FairRoot/PandaRoot
Functions
QAmacro_stt_4.C File Reference

Go to the source code of this file.

Functions

int QAmacro_stt_4 ()
 

Function Documentation

int QAmacro_stt_4 ( )

Definition at line 2 of file QAmacro_stt_4.C.

References Bool_t, ctime, Double_t, evt, fgaus, filedigi(), filereco(), filerun(), PndTrack::GetFlag(), PndSttHit::GetIsochrone(), PndMCTrack::GetMotherID(), PndTrack::GetParamLast(), PndMCTrack::GetPdgCode(), hit(), i, mctrack, mean, pnt, printf(), rtime, sigma, sigma1, t, timer, track, treedigi, treepnt, treereco, and xmax.

3 {
4  TStopwatch timer;
5  timer.Start();
6 
7  Bool_t fTest = kTRUE;
8  Int_t kindOftest[9];
9  for(int t = 0; t < 9; t++) kindOftest[t] = 0;
10 
11  // MCpoints
12  TFile filerun("testrun.root");
13  TTree *treepnt = (TTree*) filerun.Get("pndsim");
14  TClonesArray *pnt = new TClonesArray("PndSttPoint");
15  treepnt->SetBranchAddress("STTPoint",&pnt);
16 
17  // MCTracks
18  TClonesArray *mctrackarray = new TClonesArray("PndMCTrack");
19  treepnt->SetBranchAddress("MCTrack",&mctrackarray);
20 
21  // Hits
22  TFile filedigi("testdigi.root");
23  TTree *treedigi = (TTree*) filedigi.Get("pndsim");
24  TClonesArray *digiarray = new TClonesArray("PndSttHit");
25  treedigi->SetBranchAddress("STTHit",&digiarray);
26 
27  // Helix Tracks
28  TFile filereco("testreco.root");
29  TTree *treereco = (TTree*) filereco.Get("pndsim");
30  TClonesArray *trackarray = new TClonesArray("PndTrack");
31  treereco->SetBranchAddress("SttMvdTrack",&trackarray);
32 
33  // Track Match
34 // TClonesArray *trackIDarray = new TClonesArray("PndTrackID");
35 // treereco->SetBranchAddress("SttMvdTrackID",&trackIDarray);
36 
37 
38  // histograms
39  TH1F *hisochrone = new TH1F("hisochrone","Simulated drift radius", 100, 0., 0.6);
40  TH1F *hptot_mum = new TH1F("hptot_mum","Total momentum for #mu^{-}", 100, 0.6, 1.4);
41  TH1F *hptot_mup = new TH1F("hptot_mup","Total momentum for #mu^{+}", 100, 0.6, 1.4);
42  TH1F *hptot = new TH1F("hptot","Total momentum", 100, 0.6, 1.4);
43 
44  // limits
45  int nevents = treereco->GetEntriesFast();
46  int mum_tracks = 3 * nevents;
47  int mup_tracks = 3 * nevents;
48  // mu - HELIX
49  double mum_hel_res = 0.03;
50  double mum_hel_res_tol = 0.01;
51  double mum_hel_eff = 0.87;
52  double mum_hel_eff_tol = 0.1;
53  double mum_hel_peak_eff = 0.81;
54  double mum_hel_peak_eff_tol = 0.1;
55 
56  // mu + HELIX
57  double mup_hel_res = 0.03;
58  double mup_hel_res_tol = 0.01;
59  double mup_hel_eff = 0.87;
60  double mup_hel_eff_tol = 0.1;
61  double mup_hel_peak_eff = 0.81;
62  double mup_hel_peak_eff_tol = 0.1;
63 
64  // cout << "nevents " << nevents << endl;
65 
66  for (Int_t evt = 0; evt < nevents; evt++) {
67 
68  treepnt->GetEntry(evt);
69  treedigi->GetEntry(evt);
70  treereco->GetEntry(evt);
71 
72  // =================== DIGI ====================
73  for(Int_t i = 0; i < digiarray->GetEntriesFast(); i++) {
74 
75  PndSttHit *hit = (PndSttHit*) digiarray->At(i);
76  if(!hit) continue;
77  hisochrone->Fill(hit->GetIsochrone());
78  // if the reconstructed radius is > 0.5 (0.506 because of binning) -> kFALSE
79  if(hisochrone->GetBinCenter(hisochrone->FindLastBinAbove(0)) > 0.506) {
80  fTest = kFALSE; kindOftest[0] = 1;
81  cout << "TEST 0: hisocrone " << hisochrone->GetBinCenter(hisochrone->FindLastBinAbove(0)) << " (limit = 0.506)" << endl;
82  }
83  }
84 
85 
86  // =================== RECO ====================
87  // tracks loop
88  for (Int_t k = 0; k < trackarray->GetEntriesFast(); k++) {
89 
90  PndTrack *track = (PndTrack*) trackarray->At(k);
91  if(!track) continue;
92 // PndTrackID *trackID = (PndTrackID*) trackIDarray->At(k);
93 // if(!trackID) continue;
94 
95  if(track->GetFlag() < 0) continue;
96  TVector3 lastmom = track->GetParamLast().GetMomentum();
97 
98  // mu - or mu + ?
99  Int_t MCTrackID = -1;
100  if (track->GetNLinks() > 0){
101  FairLink MCLink = track->GetLink(1);
102  if (MCLink.GetType() == 0)
103  MCTrackID = MCLink.GetIndex();
104  }
105  if(MCTrackID == -1) continue;
106  PndMCTrack *mctrack = (PndMCTrack*) mctrackarray->At(MCTrackID);
107  if(!mctrack) continue;
108  if(mctrack->GetMotherID() != -1) continue;
109 
110  if(mctrack->GetPdgCode() == 13) hptot_mum->Fill(lastmom.Mag());
111  else if(mctrack->GetPdgCode() == -13) hptot_mup->Fill(lastmom.Mag());
112  hptot->Fill(lastmom.Mag());
113 
114  // else nomutracks++;
115  }
116 
117 
118  }
119 
120  // get limits on sum histo mu+ + mu-
121  int binmax = hptot->GetMaximumBin();
122  double ymax = hptot->GetMaximum();
123  double xmax = hptot->GetBinCenter(binmax);
124  int bininf = hptot->FindFirstBinAbove(ymax/2., 1);
125  bininf -= 1;
126  int binsup = hptot->FindLastBinAbove(ymax/2., 1);
127  binsup += 1;
128  double xinf = hptot->GetBinCenter(bininf);
129  double xsup = hptot->GetBinCenter(binsup);
130  double gamma = xsup - xinf;
131  double sigma1 = gamma / 2.34;
132  xinf = xmax - 3 * sigma1;
133  xsup = xmax + 3 * sigma1;
134 
135 
136 
137  // negative muon ...............................................
138  // HELIX
139  TF1 *fgaus = new TF1("fgaus", "gaus", xinf, xsup);
140  hptot_mum->Fit("fgaus","R0QN");
141  Double_t gauspar[3];
142  fgaus->GetParameters(gauspar);
143  Double_t mean = gauspar[1];
144  Double_t sigma = gauspar[2];
145 
146 
147  // mean in [1 GeV/c - 3 sigma, 1 GeV/c + 3 sigma]
148  if(mean < (mean - 3 * sigma) || mean > (mean + 3 * sigma)) {
149  fTest = kFALSE; kindOftest[1] = 1;
150  cout << "TEST 1: helix mu- mean " << mean
151  << " (limits [" << (mean - 3 * sigma) << ", " << (mean + 3 * sigma) << "])" << endl;
152  }
153 
154 // // resolution sigma/p // CHECK reactivate it when it works
155 // if((sigma/mean) < (mum_hel_res - mum_hel_res_tol) || (sigma/mean) > (mum_hel_res + mum_hel_res_tol)) {
156 // fTest = kFALSE; kindOftest[2] = 1;
157 // cout << "TEST 2: helix mu- resolution " << sigma/mean
158 // << " (limits [" << (mum_hel_res - mum_hel_res_tol) << ", " << (mum_hel_res + mum_hel_res_tol) << "])" << endl;
159 // }
160 
161 // // efficiency integral in 0.6, 1.4 / generated tracks // CHECK reactivate it when it works
162 // if((hptot_mum->Integral() / mum_tracks) < (mum_hel_eff - mum_hel_eff_tol) || (hptot_mum->Integral() / mum_tracks) > (mum_hel_eff + mum_hel_eff_tol)) {
163 // fTest = kFALSE; kindOftest[3] = 1;
164 // cout << "TEST 3: helix mu- efficiency " << (hptot_mum->Integral() / mum_tracks)
165 // << " (limits [" << (mum_hel_eff - mum_hel_eff_tol) << ", " << (mum_hel_eff + mum_hel_eff_tol) << "])" << endl;
166 // }
167 
168 // // efficiency integral under the peak/generated tracks // CHECK reactivate it when it works
169 // hptot_mum->GetXaxis()->SetRangeUser(mean - 3 * sigma, mean + 3 * sigma);
170 // if((hptot_mum->Integral() / mum_tracks) < (mum_hel_peak_eff - mum_hel_peak_eff_tol) || (hptot_mum->Integral() / mum_tracks) > (mum_hel_peak_eff + mum_hel_peak_eff_tol)) {
171 // fTest = kFALSE; kindOftest[4] = 1;
172 // cout << "TEST 4: helix mu- peak efficiency " << (hptot_mum->Integral() / mum_tracks)
173 // << " (limits [" << (mum_hel_peak_eff - mum_hel_peak_eff_tol) << ", " << (mum_hel_peak_eff + mum_hel_peak_eff_tol) << "])" << endl;
174 // }
175 
176 
177  // positive muon ...............................................
178  // HELIX
179  hptot_mup->Fit("fgaus","R0QN");
180  fgaus->GetParameters(gauspar);
181  mean = gauspar[1];
182  sigma = gauspar[2];
183 
184  // mean in [1 GeV/c - 3 sigma, 1 GeV/c + 3 sigma]
185  if(mean < (mean - 3 * sigma) || mean > (mean + 3 * sigma)) {
186  fTest = kFALSE; kindOftest[5] = 1;
187  cout << "TEST 9: helix mu+ mean " << mean
188  << " (limits [" << (mean - 3 * sigma) << ", " << (mean + 3 * sigma) << "])" << endl;
189  }
190 
191 // // resolution sigma/p p // CHECK reactivate it when it works
192 // if((sigma/mean) < (mup_hel_res - mup_hel_res_tol) || (sigma/mean) > (mup_hel_res + mup_hel_res_tol)) { fTest = kFALSE; kindOftest[6] = 1;
193 // cout << "TEST 10: helix mu+ resolution " << sigma/mean
194 // << " (limits [" << (mup_hel_res - mup_hel_res_tol) << ", " << (mup_hel_res + mup_hel_res_tol) << "])" << endl;
195 // }
196 
197 // // efficiency integral in 0.6, 1.4 / generated tracks // CHECK reactivate it when it works
198 // if((hptot_mup->Integral() / mup_tracks) < (mup_hel_eff - mup_hel_eff_tol) || (hptot_mup->Integral() / mup_tracks) > (mup_hel_eff + mup_hel_eff_tol)) {
199 // fTest = kFALSE; kindOftest[7] = 1;
200 // cout << "TEST 11: helix mu+ efficiency " << (hptot_mup->Integral() / mup_tracks)
201 // << " (limits [" << (mup_hel_eff - mup_hel_eff_tol) << ", " << (mup_hel_eff + mup_hel_eff_tol) << "])" << endl;
202 // }
203 
204  // // efficiency integral under the peak/generated tracks // CHECK reactivate it when it works
205 // hptot_mup->GetXaxis()->SetRangeUser(mean - 3 * sigma, mean + 3 * sigma);
206 // if((hptot_mup->Integral() / mup_tracks) < (mup_hel_peak_eff - mup_hel_peak_eff_tol) || (hptot_mup->Integral() / mup_tracks) > (mup_hel_peak_eff + mup_hel_peak_eff_tol)) { fTest = kFALSE; kindOftest[8] = 1;
207 // cout << "TEST 12: helix mu+ peak efficiency " << (hptot_mup->Integral() / mup_tracks)
208 // << " (limits [" << (mup_hel_peak_eff - mup_hel_peak_eff_tol) << ", " << (mup_hel_peak_eff + mup_hel_peak_eff_tol) << "])" << endl;
209 // }
210 
211 
212 
213 
214  if (fTest == kTRUE){
215  cout << " Test Passed" << endl;
216  cout << " All Ok " << endl;
217  }else{
218  cout << " Test Failed" << endl;
219  cout << " Not Ok " << endl;
220 
221  for(int t = 0; t < 9; t++) {
222 
223  if(kindOftest[t] == 0) continue;
224 
225  switch(t) {
226 
227  case 0:
228  cout << "isochrone out of 0.5 cm" << endl;
229  break;
230 
231  case 1:
232  cout << "negative muon ptot out of 3 sigma" << endl;
233  break;
234 
235  case 2:
236  cout << "negative muon ptot resolution wrong" << endl;
237  break;
238 
239  case 3:
240  cout << "negative muon ptot efficiency in [0.6, 1.4] wrong" << endl;
241  break;
242 
243  case 4:
244  cout << "negative muon ptot efficiency in peak wrong" << endl;
245  break;
246 
247  case 5:
248  cout << "positive muon ptot out of 3 sigma" << endl;
249  break;
250 
251  case 6:
252  cout << "positive muon ptot resolution wrong" << endl;
253  break;
254 
255  case 7:
256  cout << "positive muon ptot efficiency in [0.6, 1.4] wrong" << endl;
257  break;
258 
259  case 8:
260  cout << "positive muon ptot efficiency in peak wrong" << endl;
261  break;
262 
263  }
264  }
265  }
266  timer.Stop();
267  Double_t rtime = timer.RealTime();
268  Double_t ctime = timer.CpuTime();
269  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
270  return 0;
271 }
mctrack
PndMCTrack * mctrack
Definition: analysis_point_hit_fwendcap.C:250
filereco
TFile filereco("MvdStt_Test_reco.root")
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
i
Int_t i
Definition: run_full.C:25
track
PndRiemannTrack track
Definition: RiemannTest.C:33
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
xmax
Double_t xmax
Definition: convolutionAnalysis.C:11
PndTrack::GetFlag
Int_t GetFlag() const
Definition: PndTrack.h:33
pnt
TClonesArray * pnt
Definition: checkdedx_helixhit.C:12
treepnt
TTree * treepnt
Definition: checkhelixhit.C:10
evt
int evt
Definition: checkhelixhit.C:36
sigma
Double_t sigma[nsteps]
Definition: dedx_bands.C:65
PndSttHit::GetIsochrone
Double_t GetIsochrone() const
Definition: PndSttHit.h:62
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
treedigi
TTree * treedigi
Definition: checkdedx_helixhit.C:19
treereco
TTree * treereco
Definition: checkdedx_helixhit.C:25
timer
TStopwatch timer
Definition: hit_dirc.C:51
filedigi
TFile filedigi("testdigi.root")
PndTrack::GetParamLast
FairTrackParP GetParamLast()
Definition: PndTrack.h:50
ctime
Double_t ctime
Definition: hit_dirc.C:114
filerun
TFile filerun("testrun.root")
sigma1
double sigma1
Definition: reco_analys2.C:57
hit
int hit(Int_t nEvents=0, TString inFile="sim.root", TString parFile="par.root", TString inDigi="digi.root", TString outFile="hit.root", Int_t timeBased=0)
Definition: hit.C:1
t
TTree * t
Definition: bump_analys.C:13
mean
Double_t mean[nsteps]
Definition: dedx_bands.C:65
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
rtime
Double_t rtime
Definition: hit_dirc.C:113
fgaus
TF1 * fgaus[nsteps]
Definition: dedx_bands.C:62
Generated on Sun Aug 18 2019 03:02:46 for FairRoot/PandaRoot by   doxygen 1.8.8