FairRoot/PandaRoot
Functions
hit_noise_studies.C File Reference
#include <iostream>
#include <TCanvas.h>
#include <TApplication.h>
#include <TROOT.h>
#include <TChain.h>
#include <TStyle.h>
#include <TGaxis.h>
#include <TColor.h>
#include <TClonesArray.h>
#include <TH1.h>
#include <TH2.h>
#include <PndLmdDim.h>
#include <PndTrack.h>
#include <PndSdsMCPoint.h>
#include <PndSdsDigiPixel.h>
#include <PndSdsHit.h>
#include <PndSdsMergedHit.h>
#include <FairTrackParH.h>

Go to the source code of this file.

Functions

void Root_Appearance ()
 
double last_percent (-1.)
 
void DrawProgressBar (int len, double percent)
 
int hit_noise_studies ()
 
int main ()
 

Function Documentation

void DrawProgressBar ( int  len,
double  percent 
)
int hit_noise_studies ( )

Definition at line 45 of file hit_noise_studies.C.

References DrawProgressBar(), PndLmdDim::Get_instance(), PndLmdDim::Get_sensor_by_id(), PndTrack::GetParamFirst(), PndSdsHit::GetPosition(), PndSdsDigi::GetSensorID(), PndSdsHit::GetSensorID(), hit, hit_pos(), hits, nEvents, nHits, sqrt(), theta, track, and PndLmdDim::Transform_global_to_lmd_local().

Referenced by main().

45  {
46  cout << " lmd hit noise studies " << endl;
47  TCanvas* canvas = new TCanvas("canvas", "canvas", 900, 600);
48  canvas->Divide(3,2);
49 
50  // load the MC File
51  TChain tMC("pndsim");
52  tMC.Add("./Lumi_digi_0.root");
53 
54  TChain thits("pndsim");
55  thits.Add("./Lumi_recoMerged_0.root");
56 
57  TChain ttrk("pndsim");
58  ttrk.Add("./Lumi_TrackNotFiltered_0.root");
59 
60  TChain ttrk_filt("pndsim");
61  ttrk_filt.Add("./Lumi_Track_0.root");
62 
63  bool usefiltered(true);
64  TChain ttrk_prop("pndsim");
65  if (usefiltered)
66  ttrk_prop.Add("./Lumi_GeaneFiltered_0.root");
67  else
68  ttrk_prop.Add("./Lumi_Geane_0.root");
69 
70 
71  //--- assign MC info -----------------------------------------------------
72  //TClonesArray* true_tracks = new TClonesArray("PndSdsDigiPixel");
73  //tMC.SetBranchAddress("PndSdsDigiPixel", &true_tracks); //True Track to compare
74 
75  TClonesArray* digis = new TClonesArray("PndSdsDigiPixel");
76  tMC.SetBranchAddress("LMDPixelDigis", &digis); //hits
77 
78  TClonesArray* hits = new TClonesArray("PndSdsMergedHit");
79  thits.SetBranchAddress("LMDHitsMerged", &hits); //True Track to compare
80 
81  TClonesArray* tracks = new TClonesArray("PndTrack");
82  ttrk.SetBranchAddress("LMDPndTrack", &tracks);
83 
84  TClonesArray* tracks_filt = new TClonesArray("PndTrack");
85  ttrk_filt.SetBranchAddress("LMDPndTrack", &tracks_filt);
86 
87  TClonesArray* tracks_prop = new TClonesArray("FairTrackParH");
88  if (usefiltered){
89  ttrk_prop.SetBranchAddress("LMDCleanTrack", &tracks_prop);
90  cout << " using geanie filtered tracks " << endl;
91  } else
92  ttrk_prop.SetBranchAddress("GeaneTrackFinal", &tracks_prop);
93 
94 
95  //TClonesArray* tracks = new TClonesArray("LMDPndTrack");
96  //tMC.SetBranchAddress("LMDPndTrack", &tracks);
97 
98  // get the lmd dim help class
99  PndLmdDim& lmddim = PndLmdDim::Get_instance();
100 
101  // setup histograms
102  /*
103  TH1F* hist_eloss_total = new TH1F("hist_eloss_total", "energy deposit in 50#mum x 2cm x2cm sensors", 200, 0, 60);
104  TH1F* hist_eloss = new TH1F("hist_eloss", "energy deposit scaled to 1#mum", 200, 0, 1200);
105  hist_eloss_total->SetXTitle("E_{dep} [keV]");
106  hist_eloss->SetXTitle("E_{dep} [eV]");
107 
108  TH2* plane3_hits = new TH2F("plane3_hits", "hits on plane 3", 200, -10, 10, 200, -10, 10); //lmddim.Get_histogram_Plane(0, 0);
109  plane3_hits->SetXTitle("X [cm]");
110  plane3_hits->SetYTitle("Y [cm]");
111  plane3_hits->SetZTitle("# hits");
112  TH2Poly* plane3_Edep = lmddim.Get_histogram_Plane(3, 0);
113  plane3_Edep->SetNameTitle("plane3_Edep", "E dep on plane 3");
114  plane3_Edep->SetZTitle("E_{dep} [J]");
115  TH2Poly* plane3_rate = lmddim.Get_histogram_Plane(3, 0);
116  plane3_rate->SetNameTitle("plane3_rate", "rate on plane 3");
117  plane3_rate->SetZTitle("rate [kHz]");
118  TH2Poly* plane3_rad = lmddim.Get_histogram_Plane(3, 0);
119  plane3_rad->SetNameTitle("plane3_rad", "IEL dose on plane 3");
120  plane3_rad->SetZTitle("IEL dose [Gy/a]");
121 */
122  TH2D* hist_track_dir = new TH2D("hist_track_dir", "track dir", 500, -200, 200, 500, -200, 200);
123  TH2D* hist_track_dir_filt = new TH2D("hist_track_dir_filt", "track dir filtered", 500, -200, 200, 500, -200, 200);
124  TH2D* hist_track_dir_prop = new TH2D("hist_track_dir_prop", "track dir propagated", 50, -20, 20, 50, -20, 20);
125  TH1D* hist_track_theta = new TH1D("hist_track_theta", "#theta distr", 50, 0, 15);
126  // TH2D* hist_hit_mult = new TH2D("hist_hit_mult", "hit multiplicity", 100, 0, 100, 4, 0, 4);
127  TH2D* hist_hit_mult = new TH2D("hist_hit_mult", "hit multiplicity", 1000, 0, 1000, 4, 0, 4);
128  TH1D* hist_trk_mult = new TH1D("hist_trk_mult", "track multiplicity", 100, 0, 100);
129  TH2D* hist_hit_distr = new TH2D("hist_hit_distr", "hit distribution", 100, -10, 10, 100, -10, 10);
130  // read events and analyze
131  int nEvents = tMC.GetEntries();
132  int nEvents_counted(0);
133  int nhitstotal(0);
134  int ntrackstotal(0);
135  int ntrackstotal_filt(0);
136  int ntrackstotal_prop(0);
137  int ntrackstotal_rel(0);
138 
139  cout << " reading " << nEvents << " Events " << endl;
140 
141  int real_hits(0);
142  int noise_hits(0);
143  for (Int_t j = 0; j < nEvents ; j++) {
144  nEvents_counted++;
145  DrawProgressBar(50, (j + 1) / ((double) nEvents));
146  tMC.GetEntry(j);
147  ttrk.GetEntry(j);
148  ttrk_filt.GetEntry(j);
149  ttrk_prop.GetEntry(j);
150  thits.GetEntry(j);
151  ttrk_prop.GetEntry(j);
152 
153  int nTracks = tracks->GetEntriesFast();
154  hist_trk_mult->Fill(nTracks);
155  //cout << " reading " << nTracks << " tracks" << endl;
156  for (auto iTrack = 0; iTrack < nTracks; iTrack++){
157  ntrackstotal++;
158  PndTrack* track = (PndTrack*) tracks->At(iTrack);
159  if (track){
160  double pxpz = track->GetParamFirst().GetPx()/track->GetParamFirst().GetPz()*1e3;
161  pxpz -= 40;
162  double pypz = track->GetParamFirst().GetPy()/track->GetParamFirst().GetPz()*1e3;
163  hist_track_dir->Fill(pxpz, pypz);
164  double theta2 = pxpz*pxpz + pypz*pypz;
165  if ((2*2 < theta2 && theta2 < 12*12)){
166  ntrackstotal_rel++;
167  //hist_track_dir->Fill(pxpz, pypz);
168  }
169  }
170  }
171 
172  nTracks = tracks_filt->GetEntriesFast();
173  //cout << " reading " << nTracks << " tracks" << endl;
174  for (auto iTrack = 0; iTrack < nTracks; iTrack++){
175  ntrackstotal_filt++;
176  PndTrack* track = (PndTrack*) tracks->At(iTrack);
177  if (track){
178  double pxpz = track->GetParamFirst().GetPx()/track->GetParamFirst().GetPz()*1e3;
179  pxpz -= 40;
180  double pypz = track->GetParamFirst().GetPy()/track->GetParamFirst().GetPz()*1e3;
181  hist_track_dir_filt->Fill(pxpz, pypz);
182  }
183  }
184 
185  nTracks = tracks_prop->GetEntriesFast();
186  //cout << " reading " << nTracks << " tracks" << endl;
187  for (auto iTrack = 0; iTrack < nTracks; iTrack++){
188  FairTrackParH* track = (FairTrackParH*) tracks_prop->At(iTrack);
189  if (track){
190  if(track->GetLambda()!=0){ //ignore not propagated tracks
191  ntrackstotal_prop++;
192  double pxpz = track->GetPx()/track->GetPz()*1e3;
193  double pypz = track->GetPy()/track->GetPz()*1e3;
194  hist_track_dir_prop->Fill(pxpz, pypz);
195  double theta = sqrt(pxpz*pxpz+pypz*pypz);
196  hist_track_theta->Fill(theta);
197  }
198  }
199  }
200 
201  int nHits = digis->GetEntriesFast();
202  int nSdsHits = 0;
203  //cout << " reading " << nHits << " digis " << endl;
204  int ihits[4] = {0,0,0,0};
205  for (Int_t iHit = 0; iHit < nHits; iHit++) {
206  nhitstotal++;
207  PndSdsDigiPixel* mcpoint = (PndSdsDigiPixel*) digis->At(iHit);
208  //hist_hit_distr->Fill(mcpoint->Get);
209  //if (mcpoint-> < 0)
210  noise_hits++;
211  //else {
212  // real_hits++;
213  //}
214  int ihalf, iplane, imodule, iside, idie, isens;
215  int sensid = mcpoint->GetSensorID();
216  lmddim.Get_sensor_by_id(sensid, ihalf, iplane, imodule, iside, idie, isens);
217 
218  ihits[iplane]++;
219  // cout << mcpoint->GetEnergyLoss() << endl;
220  //hist_eloss_total->Fill(mcpoint->GetEnergyLoss()*1.e6); // in keV
221  //hist_eloss->Fill(mcpoint->GetEnergyLoss()*1.e9/50.); // in eV
222 /*
223  int sensID = mcpoint->GetSensorID();
224  int ihalf, iplane, imodule, iside, idie, isensor;
225  lmddim.Get_sensor_by_id(sensID, ihalf, iplane, imodule, iside, idie, isensor);
226  TVector3 hit_pos = mcpoint->GetPosition();
227  TVector3 hit_pos_lmd = lmddim.Transform_global_to_lmd_local(hit_pos, false);
228 
229  cout << mcpoint->GetTime() << endl;
230  //cout << mcpoint->Get << endl << endl;
231 
232  if (iplane == 3 && iside == 0){
233  //plane3_hits->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y());
234  //plane3_rate->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y());
235  //plane3_Edep->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y(), mcpoint->GetEnergyLoss()*1.e9*1.602e-19); // in J
236  //plane3_rad->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y(), mcpoint->GetEnergyLoss()*1.e9*1.602e-19); // in J
237  }
238 
239  if (mcpoint->GetTrackID() < 0) continue;
240  PndMCTrack *mctrk = (PndMCTrack*) true_tracks->At(
241  mcpoint->GetTrackID());
242  if (mctrk->IsGeneratorCreated()) { // take only anti-protons from the generator
243  if (mcpoint) {
244  TVector3 _mcpoint((mcpoint->GetPosition())); // position at the sensor entrance
245  TVector3 _mctrack(mcpoint->GetPx(), mcpoint->GetPy(),
246  mcpoint->GetPz()); // momentum of the track at the entrance
247 
248  TVector3 momMC = mctrk->GetMomentum(); // momentum in the primary vertex
249  TVector3 posMC = mctrk->GetStartVertex(); // position of the primary vertex
250  }
251  } */
252  }
253  hist_hit_mult->Fill(ihits[0], 0);
254  hist_hit_mult->Fill(ihits[1], 1);
255  hist_hit_mult->Fill(ihits[2], 2);
256  hist_hit_mult->Fill(ihits[3], 3);
257 
258  nHits = hits->GetEntriesFast();
259  for (auto iHit = 0; iHit < nHits; iHit++){
260  PndSdsMergedHit* hit = (PndSdsMergedHit*) hits->At(iHit);
261  int sensID = hit->GetSensorID();
262  int ihalf, iplane, imodule, iside, idie, isensor;
263  lmddim.Get_sensor_by_id(sensID, ihalf, iplane, imodule, iside, idie, isensor);
264  TVector3 hit_pos = hit->GetPosition();
265  TVector3 hit_pos_lmd = lmddim.Transform_global_to_lmd_local(hit_pos, false);
266  hist_hit_distr->Fill(hit_pos_lmd.X(), hit_pos_lmd.Y());
267  }
268  }
269  // cross section for N generated events at 1.5 GeV/c momentum was
270  // sigm_coul = 36.683 mbarn
271  // sigm_had_el = 25.3021 mbarn
272  // sigm_inel = 62.71 mbarn according to Anastasias note
273  // -> sig_total = 124.63 mbarn or 124.63 x 10-27 cm^-2
274  //double sig_total = 124.63e-27;//62e-27;//124.63e-27; // cm^2
275  // the mean instantanious luminosity at high luminosity mode is 1.6 x 10^32
276  // but we use here for reference 1 x 10^32 which can be scaled easily to any other number
277  //double lumi_ref = 1e32; // cm^-2 s^-1
278  // therefore we simulated an interaction rate of
279  //double rate_total = sig_total * lumi_ref; // 1/s or Hz
280  //cout << "\n The total rate at 1.5 GeV/c and " << lumi_ref << " mean luminosity to normalize to is " << rate_total << " Hz" << endl;
281  // The simulated time scale is then
282  //double tsim = nEvents_counted / rate_total; // s
283  cout << nEvents_counted << " events " << endl;// correspond to a real time of " << tsim << " s " << endl;
284  cout << " in " << nhitstotal << " hits found " << real_hits << " real hits and " << noise_hits << " noise hits " << endl;
285  cout << " the ratio over nevents is " << (double) noise_hits / (double) nEvents_counted << endl;
286  cout << " found " << ntrackstotal << " tracks " << (double) ntrackstotal/(double) nEvents_counted << endl;
287  cout << " found " << ntrackstotal_rel << " tracks relevant for analysis " << (double) ntrackstotal_rel/(double) nEvents_counted << endl;
288  cout << " found " << ntrackstotal_filt << " filtered tracks " << (double) ntrackstotal_filt/(double) nEvents_counted << endl;
289  cout << " found " << ntrackstotal_prop << " propagated tracks " << (double) ntrackstotal_prop/(double) nEvents_counted << endl;
290 
291  canvas->cd(1);
292  hist_track_dir->Draw("COLZ");
293  hist_track_dir->SetXTitle("p_{x}/p_{z}x10^{-3}");
294  hist_track_dir->SetYTitle("p_{y}/p_{z}x10^{-3}");
295 
296  canvas->cd(2);
297  hist_hit_distr->Draw("COLZ");
298  hist_hit_distr->SetXTitle("X [cm]");
299  hist_hit_distr->SetYTitle("Y [cm]");
300 
301  //hist_track_dir_filt->Draw("COLZ");
302  //hist_track_dir_filt->SetXTitle("p_{x}/p_{z}x10^{-3}");
303  //hist_track_dir_filt->SetYTitle("p_{y}/p_{z}x10^{-3}");
304 
305  canvas->cd(3);
306  hist_trk_mult->Draw();
307  hist_trk_mult->SetXTitle("# tracks / event");
308 
309  canvas->cd(4);
310  hist_hit_mult->Draw("colz");
311  hist_hit_mult->SetXTitle("# hits / event");
312  hist_hit_mult->SetYTitle("plane");
313 
314  canvas->cd(5);
315  hist_track_dir_prop->Draw("COLZ");
316  hist_track_dir_prop->SetXTitle("p_{x}/p_{z}x10^{-3}");
317  hist_track_dir_prop->SetYTitle("p_{y}/p_{z}x10^{-3}");
318 
319  canvas->cd(6);
320  hist_track_theta->Draw();
321  hist_track_theta->SetXTitle("#Theta [mrad]");
322  /*
323  plane3_rate->Scale(1./tsim*1e-3); // in kHz
324  //plane3_rate->SetMaximum(150);
325  cout << " total rate of " << nhitstotal << " hits is " << nhitstotal/tsim*1e-6 << "MHz" << endl;
326  plane3_rad->Scale(1./tsim*60.*60.*24.*365.*0.5/(2.33e-3*(2*2*50.e-4))); // 50% duty cycle over a year and normalized to the mass in kg of one sensor
327  //canvas->cd(1);
328  plane3_Edep->Draw("colz");
329  //hist_eloss_total->Draw();
330  //canvas->cd(2);
331  plane3_hits->Draw("colz");
332  canvas->Print("MC_hit_distr_1_5_plane_3_side_0.pdf");
333  canvas->Print("MC_hit_distr_1_5_plane_3_side_0.root");
334  hist_eloss->Draw();
335  canvas->Print("MC_IEL_distrib_1_5.pdf");
336  canvas->Print("MC_IEL_distrib_1_5.root");
337  //canvas->cd(3);
338  plane3_rate->Draw("colz");
339  cout << " max rate is " << plane3_rate->GetBinContent(plane3_rate->GetMaximumBin()) << " kHz" << endl;
340  canvas->Print("MC_hit_rate_1_5_plane_3_side_0.pdf");
341  canvas->Print("MC_hit_rate_1_5_plane_3_side_0.root");
342  //canvas->cd(4);
343  plane3_rad->Draw("colz");
344  cout << " max dose is " << plane3_rad->GetBinContent(plane3_rad->GetMaximumBin()) << " Gy" << endl;
345  canvas->Print("MC_IEL_dose_1_5_plane_3_side_0.pdf");
346  canvas->Print("MC_IEL_does_1_5_plane_3_side_0.root");
347  */
348  return 0;
349 }
theta
double theta
Definition: macro/detectors/emc/dedicated/reco_analys.C:24
hit_pos
TVector3 hit_pos(hit->GetX(), hit->GetY(), hit->GetZ())
PndSdsDigi::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsDigi.h:59
PndSdsHit::GetPosition
TVector3 GetPosition() const
Definition: PndSdsHit.h:93
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
cd
f cd()
PndLmdDim::Transform_global_to_lmd_local
TVector3 Transform_global_to_lmd_local(const TVector3 &point, bool isvector=false, bool aligned=true)
Definition: PndLmdDim.cxx:3113
DrawProgressBar
void DrawProgressBar(int len, double percent)
Definition: geo/Check_particle_path.C:447
cm
Double_t const cm
Definition: createRootGeoFile.C:31
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
SetXTitle
hist4 SetXTitle("p/(GeV/c)")
Y
double Y
Definition: anaLmdDigi.C:68
for
for(int j=0;j< ncounts;j++)
Definition: anal_hit_digi_cluster_fwendcap.C:156
Draw
tree Draw("fELoss:TMath::Sqrt(fPx_out*fPx_out+ fPy_out*fPy_out+ fPz_out*fPz_out)", ppos &&"fELoss < 0.04")
analysis
int analysis()
Definition: analysis.C:8
nHits
int nHits
Definition: RiemannTest.C:16
PndLmdDim::Get_sensor_by_id
void Get_sensor_by_id(const int sensor_id, int &ihalf, int &iplane, int &imodule, int &iside, int &idie, int &isensor)
Definition: PndLmdDim.h:526
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
track
PndMCTrack * track
Definition: anaLmdCluster.C:89
z
Double_t z
Definition: createRootFscGeometryFile.C:146
total
TF1 * total
Definition: anaLmdCluster.C:147
is
static int is
Definition: ranlxd.cxx:374
X
double X
Definition: anaLmdDigi.C:68
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndLmdDim::Get_instance
static PndLmdDim & Get_instance()
Definition: PndLmdDim.cxx:242
PndSdsDigiPixel
Data class to store the digi output of a pixel module.
Definition: PndSdsDigiPixel.h:25
hit
PndSdsMCPoint * hit
Definition: anasim.C:70
hits
CbmHit * hits[nHits]
Definition: RiemannTest.C:19
y
Double_t y
Definition: createRootFscGeometryFile.C:159
PndSdsHit::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
PndTrack::GetParamFirst
FairTrackParP GetParamFirst()
Definition: PndTrack.h:49
double last_percent ( -  1.)
int main ( void  )

Definition at line 37 of file hit_noise_studies.C.

References hit_noise_studies(), and Root_Appearance().

37  {
38  Root_Appearance();
39  TApplication myapp("myapp", 0, 0);
40  hit_noise_studies();
41  myapp.Run();
42  return 0;
43 }
hit_noise_studies
int hit_noise_studies()
Definition: hit_noise_studies.C:45
Root_Appearance
void Root_Appearance()
Definition: hit_noise_studies.C:351
void Root_Appearance ( )

Definition at line 351 of file hit_noise_studies.C.

References Double_t.

Referenced by main(), and online_monitoring_studies().

351  {
352  TPad foo; // never remove this line :-)))
353  if(1){
354  gROOT->SetStyle("Plain");
355  const Int_t NRGBs = 5;
356  const Int_t NCont = 255;
357  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
358  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
359  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
360  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
361  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
362  gStyle->SetNumberContours(NCont);
363  gStyle->SetTitleFont(10*13+2,"xyz");
364  gStyle->SetTitleSize(0.06, "xyz");
365  gStyle->SetTitleOffset(1,"xy");
366  gStyle->SetTitleOffset(1.3,"z");
367  gStyle->SetLabelFont(10*13+2,"xyz");
368  gStyle->SetLabelSize(0.06,"xyz");
369  gStyle->SetLabelOffset(0.009,"xyz");
370  gStyle->SetPadBottomMargin(0.2);
371  gStyle->SetPadTopMargin(0.15);
372  gStyle->SetPadLeftMargin(0.15);
373  gStyle->SetPadRightMargin(0.20);
374  gStyle->SetOptTitle(0);
375  gStyle->SetOptStat(0);
376  gROOT->ForceStyle();
377  gStyle->SetFrameFillColor(0);
378  gStyle->SetFrameFillStyle(0);
379  TGaxis::SetMaxDigits(3);
380  }
381 }
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
Generated on Wed Apr 3 2019 08:44:18 for FairRoot/PandaRoot by   doxygen 1.8.5