]>
Commit | Line | Data |
---|---|---|
10a5e8ec | 1 | #ifndef __CINT__ |
2 | #include <iostream.h> | |
3 | #include "AliTRDtracker.h" | |
4 | #include "AliTRDcluster.h" | |
5 | #include "AliTRDv1.h" | |
6 | #include "AliTRDgeometry.h" | |
7 | #include "AliTRDparameter.h" | |
8 | #include "alles.h" | |
9 | #include "AliTRDmcTrack.h" | |
10 | #include "AliTRDtrack.h" | |
11 | ||
12 | #include "TFile.h" | |
13 | #include "TParticle.h" | |
14 | #include "TStopwatch.h" | |
15 | #endif | |
16 | ||
17 | ||
18 | void AliTRDseedsAnalysis() { | |
19 | ||
20 | const Int_t nPrimaries = (Int_t) 3*86030/4; | |
21 | // const Int_t nPrimaries = 500; | |
22 | ||
23 | Bool_t page[6]; for(Int_t i=0; i<6; i++) page[i]=kTRUE; | |
24 | ||
25 | const Int_t nPtSteps = 30; | |
26 | const Float_t maxPt = 3.; | |
27 | const Float_t minPt = 0.; | |
28 | ||
29 | const Int_t nEtaSteps = 40; | |
30 | const Float_t maxEta = 1.; | |
31 | const Float_t minEta = -1.; | |
32 | ||
33 | // page[0] | |
34 | TH1F *hNcl=new TH1F("hNcl","Seeds No of Clusters",255,-9.5,500.5); | |
35 | hNcl->SetFillColor(4); | |
36 | hNcl->SetXTitle("No of clusters"); | |
37 | hNcl->SetYTitle("counts"); | |
38 | TH1F *hNtb=new TH1F("hNtb","Seeds No of TB with clusters",160,-9.5,150.5); | |
39 | hNtb->SetFillColor(4); | |
40 | hNtb->SetXTitle("No of timebins with clusters"); | |
41 | hNtb->SetYTitle("counts"); | |
42 | TH1F *hNep=new TH1F("hNep","Seeds end point",160, -9.5, 150.5); | |
43 | hNep->SetFillColor(4); | |
44 | hNep->SetXTitle("outermost timebin with cluster"); | |
45 | hNep->SetYTitle("counts"); | |
46 | TH1F *hNmg=new TH1F("hNmg","Seeds max gap",160, -9.5, 150.5); | |
47 | hNmg->SetFillColor(4); | |
48 | hNmg->SetXTitle("max gap (tb)"); | |
49 | hNmg->SetYTitle("counts"); | |
50 | ||
51 | // page[1] | |
52 | Float_t fMin = -5.5, fMax = 134.5; | |
53 | Int_t iChan = 140; | |
54 | TH2F *h2ep = new TH2F("h2ep","MC vs RT (end point)",iChan,fMin,fMax,iChan,fMin,fMax); | |
55 | h2ep->SetMarkerColor(4); | |
56 | h2ep->SetMarkerSize(2); | |
57 | TH2F *h2ntb = new TH2F("h2ntb","MC vs RT (TB with Clusters)",iChan,fMin,fMax,iChan,fMin,fMax); | |
58 | h2ntb->SetMarkerColor(4); | |
59 | h2ntb->SetMarkerSize(2); | |
60 | TH2F *h2mg = new TH2F("h2mg","MC vs RT (max gap)",iChan/2,fMin,fMax/2,iChan/2,fMin,fMax/2); | |
61 | h2mg->SetMarkerColor(4); | |
62 | h2mg->SetMarkerSize(2); | |
63 | ||
64 | // page[2] | |
65 | TH1F *hPt_all=new TH1F("hPt_all","Seeds Pt",nPtSteps,minPt,maxPt); | |
66 | hPt_all->SetLineColor(4); | |
67 | hPt_all->SetXTitle("Pt (GeV/c)"); | |
68 | hPt_all->SetYTitle("counts"); | |
69 | ||
70 | TH1F *hPt_short=new TH1F("hPt_short","Short seeds Pt",nPtSteps,minPt,maxPt); | |
71 | hPt_short->SetLineColor(8); | |
72 | TH1F *hPt_long=new TH1F("hPt_long","Long seeds Pt",nPtSteps,minPt,maxPt); | |
73 | hPt_long->SetLineColor(2); | |
74 | ||
75 | TH1F *hrtPt_short=new TH1F("hrtPt_short","RT from short seeds",nPtSteps,minPt,maxPt); | |
76 | hrtPt_short->SetLineColor(8); | |
77 | TH1F *hrtPt_long=new TH1F("hrtPt_long","RT from long seeds",nPtSteps,minPt,maxPt); | |
78 | hrtPt_long->SetLineColor(2); | |
79 | ||
80 | TH1F *hEta_all=new TH1F("hEta_all","Seeds Eta",nEtaSteps,minEta,maxEta); | |
81 | hEta_all->SetLineColor(4); | |
82 | hEta_all->SetXTitle("Eta"); | |
83 | hEta_all->SetYTitle("counts"); | |
84 | TH1F *hEta_short=new TH1F("hEta_short","Short seeds Eta",nEtaSteps,minEta,maxEta); | |
85 | hEta_short->SetLineColor(8); | |
86 | TH1F *hEta_long=new TH1F("hEta_long","Long seeds Eta",nEtaSteps,minEta,maxEta); | |
87 | hEta_long->SetLineColor(2); | |
88 | ||
89 | TH1F *hrtEta_short=new TH1F("hrtEta_short","RT from short seeds",nEtaSteps,minEta,maxEta); | |
90 | hrtEta_short->SetLineColor(8); | |
91 | TH1F *hrtEta_long=new TH1F("hrtEta_long","RT from long seeds",nEtaSteps,minEta,maxEta); | |
92 | hrtPt_long->SetLineColor(2); | |
93 | ||
94 | // page[3] | |
95 | TH1F *hEff_long = new TH1F("hEff_long","Efficiency vs Pt",nPtSteps,minPt,maxPt); | |
96 | hEff_long->SetLineColor(2); hEff_long->SetLineWidth(2); | |
97 | hEff_long->SetXTitle("Pt, GeV/c"); | |
98 | hEff_long->SetYTitle("Efficiency"); | |
99 | ||
100 | TH1F *hEff_short = new TH1F("hEff_short","Efficiency short",nPtSteps,minPt,maxPt); | |
101 | hEff_short->SetLineColor(8); hEff_short->SetLineWidth(2); | |
102 | hEff_short->SetXTitle("Pt, GeV/c"); | |
103 | hEff_short->SetYTitle("Efficiency"); | |
104 | ||
105 | TH1F *hEff_long_eta = new TH1F("hEff_long_eta","Efficiency vs Eta",nEtaSteps,minEta,maxEta); | |
106 | hEff_long_eta->SetLineColor(2); hEff_long_eta->SetLineWidth(2); | |
107 | hEff_long_eta->SetXTitle("Pt, GeV/c"); | |
108 | hEff_long_eta->SetYTitle("Efficiency"); | |
109 | ||
110 | TH1F *hEff_short_eta = new TH1F("hEff_short_eta","Efficiency short",nEtaSteps,minEta,maxEta); | |
111 | hEff_short_eta->SetLineColor(8); hEff_short_eta->SetLineWidth(2); | |
112 | hEff_short_eta->SetXTitle("Pt, GeV/c"); | |
113 | hEff_short_eta->SetYTitle("Efficiency"); | |
114 | ||
115 | // page[4] | |
116 | TH1F *hY=new TH1F("hY","Y resolution",50,-0.5,0.5);hY->SetLineColor(4); | |
117 | hY->SetLineWidth(2); | |
118 | hY->SetXTitle("(cm)"); | |
119 | TH1F *hZ=new TH1F("hZ","Z resolution",80,-4,4);hZ->SetLineColor(4); | |
120 | hZ->SetLineWidth(2); | |
121 | hZ->SetXTitle("(cm)"); | |
122 | TH1F *hPhi=new TH1F("hPhi","PHI resolution",100,-20.,20.); | |
123 | hPhi->SetFillColor(4); | |
124 | hPhi->SetXTitle("(mrad)"); | |
125 | TH1F *hLambda=new TH1F("hLambda","Lambda resolution",100,-100,100); | |
126 | hLambda->SetFillColor(17); | |
127 | hLambda->SetXTitle("(mrad)"); | |
128 | TH1F *hPt=new TH1F("hPt","Relative Pt resolution",30,-10.,10.); | |
129 | ||
130 | // page[5] | |
131 | TH1F *hY_n=new TH1F("hY_n","Normalized Y resolution",50,-0.5,0.5); hY_n->SetLineColor(4); | |
132 | hY_n->SetLineWidth(2); | |
133 | hY_n->SetXTitle(" "); | |
134 | TH1F *hZ_n=new TH1F("hZ_n","Normalized Z resolution",50,-0.5,0.5); hZ_n->SetLineColor(2); | |
135 | hZ_n->SetLineWidth(2); | |
136 | hZ_n->SetXTitle(" "); | |
137 | TH1F *hLambda_n=new TH1F("hLambda_n","Normalized Lambda resolution",50,-0.5,0.5); | |
138 | hLambda_n->SetFillColor(17); | |
139 | TH1F *hPt_n=new TH1F("hPt_n","Normalized Pt resolution",50,-0.5,0.5); | |
140 | ||
141 | ||
142 | Int_t nEvent = 0; | |
143 | ||
144 | // Load Seeds saved as MC tracks | |
145 | TObjArray mctarray(2000); | |
146 | TFile *mctf=TFile::Open("AliTRDtrackableSeeds.root"); | |
147 | if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDtrackableSeeds.root !\n"; return;} | |
148 | TTree *mctracktree=(TTree*)mctf->Get("MCtracks"); | |
149 | TBranch *mctbranch=mctracktree->GetBranch("MCtracks"); | |
150 | Int_t const nMCtracks = (Int_t) mctracktree->GetEntries(); | |
151 | cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl; | |
152 | for (Int_t i=0; i<nMCtracks; i++) { | |
153 | AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack(); | |
154 | mctbranch->SetAddress(&ioMCtrack); | |
155 | mctracktree->GetEvent(i); | |
156 | mctarray.AddLast(ioMCtrack); | |
157 | } | |
158 | mctf->Close(); | |
159 | ||
160 | TFile *tf=TFile::Open("AliTRDtracks.root"); | |
161 | ||
162 | if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;} | |
163 | TObjArray rtarray(2000); | |
164 | ||
165 | char tname[100]; | |
166 | sprintf(tname,"TRDb_%d",nEvent); | |
167 | TTree *tracktree=(TTree*)tf->Get(tname); | |
168 | ||
169 | TBranch *tbranch=tracktree->GetBranch("tracks"); | |
170 | ||
171 | Int_t nRecTracks = (Int_t) tracktree->GetEntries(); | |
172 | cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl; | |
173 | ||
174 | for (Int_t i=0; i<nRecTracks; i++) { | |
175 | AliTRDtrack *iotrack=new AliTRDtrack(); | |
176 | tbranch->SetAddress(&iotrack); | |
177 | tracktree->GetEvent(i); | |
178 | rtarray.AddLast(iotrack); | |
179 | } | |
180 | tf->Close(); | |
181 | ||
182 | // return; | |
183 | ||
184 | AliTRDmcTrack *mct; | |
185 | AliTRDtrack *rt; | |
186 | Int_t rtIndex[nMCtracks]; | |
187 | for(Int_t j = 0; j < nMCtracks; j++) { | |
188 | mct = (AliTRDmcTrack*) mctarray.UncheckedAt(j); | |
189 | rtIndex[j] = 0; | |
190 | for (Int_t i=0; i<nRecTracks; i++) { | |
191 | rt = (AliTRDtrack*) rtarray.UncheckedAt(i); | |
192 | Int_t label = rt->GetSeedLabel(); | |
193 | if(mct->GetTrackIndex() == label) rtIndex[j] = i; | |
194 | } | |
195 | } | |
196 | ||
197 | // Load clusters | |
198 | ||
199 | TFile *geofile =TFile::Open("AliTRDclusters.root"); | |
200 | AliTRDtracker *Tracker = new AliTRDtracker(geofile); | |
201 | Tracker->SetEventNumber(nEvent); | |
202 | ||
203 | AliTRDgeometry *fGeom = (AliTRDgeometry*) geofile->Get("TRDgeometry"); | |
204 | AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter"); | |
205 | ||
206 | Char_t *alifile = "AliTRDclusters.root"; | |
207 | TObjArray carray(2000); | |
208 | TObjArray *ClustersArray = &carray; | |
209 | Tracker->ReadClusters(ClustersArray,alifile); | |
210 | ||
211 | ||
212 | // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits | |
213 | alifile = "galice.root"; | |
214 | TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile); | |
215 | if (!gafl) { | |
216 | cout << "Open the ALIROOT-file " << alifile << endl; | |
217 | gafl = new TFile(alifile); | |
218 | } | |
219 | else { | |
220 | cout << alifile << " is already open" << endl; | |
221 | } | |
222 | gAlice = (AliRun*) gafl->Get("gAlice"); | |
223 | if (gAlice) | |
224 | cout << "AliRun object found on file" << endl; | |
225 | else | |
226 | gAlice = new AliRun("gAlice","Alice test program"); | |
227 | ||
228 | AliTRDv1 *TRD = (AliTRDv1*) gAlice->GetDetector("TRD"); | |
229 | const Int_t nparticles = gAlice->GetEvent(nEvent); | |
230 | if (nparticles <= 0) return; | |
231 | ||
232 | TParticle *p; | |
233 | Bool_t electrons[300000]; | |
234 | for(Int_t i = 0; i < 300000; i++) electrons[i] = kFALSE; | |
235 | ||
236 | Bool_t mark_electrons = kFALSE; | |
237 | if(mark_electrons) { | |
238 | printf("mark electrons\n"); | |
239 | ||
240 | for(Int_t i = nPrimaries; i < nparticles; i++) { | |
241 | p = gAlice->Particle(i); | |
242 | if(p->GetMass() > 0.01) continue; | |
243 | if(p->GetMass() < 0.00001) continue; | |
244 | electrons[i] = kTRUE; | |
245 | } | |
246 | } | |
247 | ||
248 | AliTRDcluster *cl = 0; | |
249 | Int_t nw, label, index, ti[3]; | |
250 | Int_t mct_tbwc, rt_tbwc, mct_max_gap, rt_max_gap, mct_sector, rt_sector; | |
251 | Int_t mct_max_tb, rt_max_tb, mct_min_tb, rt_min_tb; | |
252 | ||
253 | Int_t det, plane, ltb, gtb, gap, max_gap, sector; | |
254 | Double_t Pt, Px, Py, Pz; | |
255 | Double_t mcPt, mcPx, mcPy, mcPz; | |
256 | Double_t x,y,z, mcX; | |
257 | Int_t rtClusters, rtCorrect; | |
258 | ||
259 | Int_t nToFind_long, nFound_long, nToFind_short, nFound_short; | |
260 | ||
261 | printf("\n"); | |
262 | ||
263 | Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region | |
264 | Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region | |
265 | Double_t dx = (Double_t) fPar->GetTimeBinSize(); | |
266 | ||
267 | Int_t tbAmp = fPar->GetTimeBefore(); | |
268 | Int_t maxAmp = 0; | |
269 | Int_t tbDrift = fPar->GetTimeMax(); | |
270 | Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx); | |
271 | ||
272 | tbDrift = TMath::Min(tbDrift,maxDrift); | |
273 | tbAmp = TMath::Min(tbAmp,maxAmp); | |
274 | ||
275 | const Int_t nPlanes = fGeom->Nplan(); | |
276 | const Int_t tbpp = Tracker->GetTimeBinsPerPlane(); | |
277 | const Int_t nTB = tbpp * nPlanes; | |
278 | Int_t mask[nTB][2]; | |
279 | ||
280 | nToFind_long = 0; nFound_long = 0; | |
281 | nToFind_short = 0; nFound_short = 0; | |
282 | ||
283 | for(Int_t i=0; i < nMCtracks; i++) { | |
284 | mct = (AliTRDmcTrack *) mctarray.UncheckedAt(i); | |
285 | label = mct->GetTrackIndex(); | |
286 | ||
287 | // Waveform of the MC track | |
288 | for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; } | |
289 | Int_t mct_ncl = mct->GetNumberOfClusters(); | |
290 | ||
291 | for(ltb = 0; ltb < kMAX_TB; ltb++) { | |
292 | for(plane = 0; plane < nPlanes; plane++) { | |
293 | for(Int_t n = 0; n < 2; n++) { | |
294 | index = mct->GetClusterIndex(ltb,plane,n); | |
295 | if(index < 0) continue; | |
296 | cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index); | |
297 | ||
298 | for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw); | |
299 | ||
300 | if((ti[0] != label) && (ti[1] != label) && (ti[2] != label)) { | |
301 | printf("mc wrong track label: %d, %d, %d != %d \n", | |
302 | ti[0], ti[1], ti[2], label); | |
303 | } | |
304 | det=cl->GetDetector(); | |
305 | if(fGeom->GetPlane(det) != plane) | |
306 | printf("cluster plane = %d != %d expected plane\n", | |
307 | fGeom->GetPlane(det), plane); | |
308 | if(cl->GetLocalTimeBin() != ltb) | |
309 | printf("cluster ltb = %d != %d expected ltb\n", | |
310 | cl->GetLocalTimeBin(), ltb); | |
311 | gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
312 | mask[gtb][n] = index; | |
313 | } | |
314 | } | |
315 | } | |
316 | ||
317 | for(plane = 0; plane < nPlanes; plane++) { | |
318 | for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) { | |
319 | gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
320 | if(mask[gtb][0] > -1) break; | |
321 | } | |
322 | if(ltb >= -tbAmp) break; | |
323 | } | |
324 | if((plane == nPlanes) && (ltb == -tbAmp-1)) { | |
325 | // printf("warning: for track %d min tb is not found and set to %d!\n", | |
326 | // label, nTB-1); | |
327 | mct_min_tb = nTB-1; | |
328 | // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]); | |
329 | } | |
330 | else { | |
331 | mct_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
332 | } | |
333 | ||
334 | for(plane = nPlanes-1 ; plane>=0; plane--) { | |
335 | for(ltb = -tbAmp; ltb < tbDrift; ltb++) { | |
336 | gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
337 | if(mask[gtb][0] > -1) break; | |
338 | } | |
339 | if(ltb < tbDrift) break; | |
340 | } | |
341 | if((plane == -1) && (ltb == tbDrift)) { | |
342 | // printf("warning: for track %d max tb is not found and set to 0!\n",label); | |
343 | // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]); | |
344 | mct_max_tb = 0; | |
345 | // mcX = Tracker->GetX(0,0,tbDrift-1); | |
346 | } | |
347 | else { | |
348 | mct_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb); | |
349 | // mcX = Tracker->GetX(0,plane,ltb); | |
350 | cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]); | |
351 | det=cl->GetDetector(); | |
352 | sector = fGeom->GetSector(det); | |
353 | mct_sector = sector; | |
354 | } | |
355 | ||
356 | mct_tbwc = 0; | |
357 | mct_max_gap = 0; | |
358 | gap = -1; | |
359 | ||
360 | for(nw = mct_min_tb; nw < mct_max_tb+1; nw++) { | |
361 | gap++; | |
362 | if(mask[nw][0] > -1) { | |
363 | mct_tbwc++; | |
364 | if(gap > mct_max_gap) mct_max_gap = gap; | |
365 | gap = 0; | |
366 | } | |
367 | } | |
368 | ||
369 | // Waveform of the reconstructed track | |
370 | if(rtIndex[i] >= 0) rt = (AliTRDtrack *) rtarray.UncheckedAt(rtIndex[i]); | |
371 | ||
372 | for(nw = 0; nw < nTB; nw++) { mask[nw][0] = -1; mask[nw][1] = -1; } | |
373 | Int_t rt_ncl = rt->GetNumberOfClusters(); | |
374 | ||
375 | for(Int_t n = 0; n < rt_ncl; n++) { | |
376 | index = rt->GetClusterIndex(n); | |
377 | cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index); | |
378 | ||
379 | for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw); | |
380 | ||
381 | if((ti[0] != label) && (ti[1] != label) && (ti[2] != label)) { | |
382 | // printf("rt wrong track label: %d, %d, %d != %d \n", ti[0], ti[1], ti[2], label); | |
383 | continue; | |
384 | } | |
385 | ||
386 | det=cl->GetDetector(); | |
387 | plane = fGeom->GetPlane(det); | |
388 | ltb = cl->GetLocalTimeBin(); | |
389 | gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
390 | mask[gtb][0] = index; | |
391 | } | |
392 | ||
393 | for(plane = 0; plane < nPlanes; plane++) { | |
394 | for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) { | |
395 | gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
396 | if(mask[gtb][0] > -1) break; | |
397 | } | |
398 | if(ltb >= -tbAmp) break; | |
399 | } | |
400 | if((plane == nPlanes) && (ltb == -tbAmp-1)) { | |
401 | // printf("warning: for track %d min tb is not found and set to %d!\n", | |
402 | // label, nTB-1); | |
403 | rt_min_tb = nTB-1; | |
404 | // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]); | |
405 | } | |
406 | else { | |
407 | rt_min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
408 | } | |
409 | ||
410 | for(plane = nPlanes-1 ; plane>=0; plane--) { | |
411 | for(ltb = -tbAmp; ltb < tbDrift; ltb++) { | |
412 | gtb = Tracker->GetGlobalTimeBin(0,plane,ltb); | |
413 | if(mask[gtb][0] > -1) break; | |
414 | } | |
415 | if(ltb < tbDrift) break; | |
416 | } | |
417 | if((plane == -1) && (ltb == tbDrift)) { | |
418 | // printf("warning: for track %d max tb is not found and set to 0!\n",label); | |
419 | // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]); | |
420 | rt_max_tb = 0; | |
421 | // mcX = Tracker->GetX(0,0,tbDrift-1); | |
422 | } | |
423 | else { | |
424 | rt_max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb); | |
425 | // mcX = Tracker->GetX(0,plane,ltb); | |
426 | cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb][0]); | |
427 | det=cl->GetDetector(); | |
428 | sector = fGeom->GetSector(det); | |
429 | rt_sector = sector; | |
430 | } | |
431 | ||
432 | rt_tbwc = 0; | |
433 | rt_max_gap = 0; | |
434 | gap = -1; | |
435 | ||
436 | for(nw = rt_min_tb; nw < rt_max_tb+1; nw++) { | |
437 | gap++; | |
438 | if(mask[nw][0] > -1) { | |
439 | rt_tbwc++; | |
440 | if(gap > rt_max_gap) rt_max_gap = gap; | |
441 | gap = 0; | |
442 | } | |
443 | } | |
444 | ||
445 | // Fill the histoes | |
446 | ||
447 | if(page[0]) { | |
448 | hNcl->Fill((Float_t) mct_ncl); | |
449 | hNtb->Fill((Float_t) mct_tbwc); | |
450 | hNep->Fill((Float_t) mct_max_tb); | |
451 | hNmg->Fill((Float_t) mct_max_gap); | |
452 | } | |
453 | if(page[1]) { | |
454 | h2ep->Fill((Float_t) rt_max_tb, (Float_t) mct_max_tb); | |
455 | h2ntb->Fill((Float_t) rt_tbwc, (Float_t) mct_tbwc); | |
456 | h2mg->Fill((Float_t) rt_max_gap, (Float_t) mct_max_gap); | |
457 | } | |
458 | if(page[2]) { | |
459 | p = gAlice->Particle(label); | |
460 | hPt_all->Fill(p->Pt()); | |
461 | hEta_all->Fill(p->Eta()); | |
462 | if(mct_max_tb > 60) { | |
463 | nToFind_long++; | |
464 | hPt_long->Fill(p->Pt()); | |
465 | hEta_long->Fill(p->Eta()); | |
466 | if(((mct_max_tb - rt_max_tb) < 10) && | |
467 | (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) { | |
468 | nFound_long++; | |
469 | hrtPt_long->Fill(p->Pt()); | |
470 | hrtEta_long->Fill(p->Eta()); | |
471 | } | |
472 | } | |
473 | if((mct_max_tb < 60) && (mct_max_tb > 10)) { | |
474 | nToFind_short++; | |
475 | hPt_short->Fill(p->Pt()); | |
476 | hEta_short->Fill(p->Eta()); | |
477 | if(((mct_max_tb - rt_max_tb) < 10) && | |
478 | (((Float_t) rt_tbwc) / ((Float_t) mct_tbwc) > 0.7)) { | |
479 | nFound_short++; | |
480 | hrtPt_short->Fill(p->Pt()); | |
481 | hrtEta_short->Fill(p->Eta()); | |
482 | } | |
483 | } | |
484 | } | |
485 | if(page[4] && page[5]) { | |
486 | if((mct_tbwc > 50) && (rt_tbwc > 50)) { | |
487 | if(rt->GetSeedLabel() != label) { | |
488 | printf("mct vs rt index mismatch: %d != %d \n", | |
489 | label, rt->GetSeedLabel()); | |
490 | return; | |
491 | } | |
492 | Pt = TMath::Abs(rt->GetPt()); | |
493 | Double_t cc = TMath::Abs(rt->GetSigmaC2()); | |
494 | mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1); | |
495 | rt->GetPxPyPz(Px,Py,Pz); | |
496 | ||
497 | printf("\n\ntrack %d \n", label); | |
498 | printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz); | |
499 | printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz); | |
500 | ||
501 | mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy); | |
502 | if(mcPt < 0.0001) mcPt = 0.0001; | |
503 | ||
504 | Float_t lamg=TMath::ATan2(mcPz,mcPt); | |
505 | Float_t lam =TMath::ATan2(Pz,Pt); | |
506 | if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt); | |
507 | else hPt_n->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc)); | |
508 | ||
509 | Float_t phig=TMath::ATan2(mcPy,mcPx); | |
510 | Float_t phi =TMath::ATan2(Py,Px); | |
511 | // if(!(rt->PropagateTo(x))) continue; | |
512 | // if((mcPt > Pt_min) && (mcPt < Pt_max)) { | |
513 | hLambda->Fill((lam - lamg)*1000.); | |
514 | hLambda_n->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2())); | |
515 | if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt); | |
516 | else hPt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.); | |
517 | hPhi->Fill((phi - phig)*1000.); | |
518 | hY->Fill(rt->GetY() - y); | |
519 | hZ->Fill(rt->GetZ() - z); | |
520 | hY_n->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2())); | |
521 | hZ_n->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2())); | |
522 | } | |
523 | } | |
524 | } | |
525 | ||
526 | // plot pictures | |
527 | ||
528 | if(page[0]) { | |
529 | TCanvas *c0=new TCanvas("c0","",0,0,700,850); | |
530 | gStyle->SetOptStat(111110); | |
531 | c0->Divide(2,2); | |
532 | c0->cd(1); gPad->SetLogy(); hNcl->Draw(); | |
533 | c0->cd(2); hNtb->Draw(); | |
534 | c0->cd(3); hNep->Draw(); | |
535 | c0->cd(4); hNmg->Draw(); | |
536 | c0->Print("c0.ps","ps"); | |
537 | } | |
538 | if(page[1]) { | |
539 | TCanvas *c1=new TCanvas("c1","",0,0,700,850); | |
540 | gStyle->SetOptStat(0); | |
541 | c1->Divide(2,2); | |
542 | c1->cd(1); h2ep->Draw(); | |
543 | c1->cd(2); h2ntb->Draw(); | |
544 | c1->cd(3); h2mg->Draw(); | |
545 | // c1->cd(4); hNmg->Draw(); | |
546 | c1->Print("c1.ps","ps"); | |
547 | } | |
548 | if(page[2]) { | |
549 | TCanvas *c2=new TCanvas("c2","",0,0,700,850); | |
550 | gStyle->SetOptStat(0); | |
551 | c2->Divide(2,2); | |
552 | c2->cd(1); hPt_all->Draw(); hPt_long->Draw("same"); hPt_short->Draw("same"); | |
553 | c2->cd(2); hEta_all->Draw(); hEta_long->Draw("same"); hEta_short->Draw("same"); | |
554 | // c2->cd(3); h2mg->Draw(); | |
555 | // c2->cd(4); hNmg->Draw(); | |
556 | c2->Print("c2.ps","ps"); | |
557 | } | |
558 | if(page[3]) { | |
559 | TCanvas *c3=new TCanvas("c3","",0,0,700,850); | |
560 | gStyle->SetOptStat(0); | |
561 | c3->Divide(1,2); | |
562 | c3->cd(1); | |
563 | hrtPt_long->Sumw2(); hPt_long->Sumw2(); | |
564 | hEff_long->Divide(hrtPt_long,hPt_long,1,1.,"b"); | |
565 | hEff_long->SetMaximum(1.4); | |
566 | hEff_long->SetYTitle("Matching Efficiency"); | |
567 | hEff_long->SetXTitle("Pt (GeV/c)"); | |
568 | hEff_long->Draw(); | |
569 | hrtPt_short->Sumw2(); hPt_short->Sumw2(); | |
570 | hEff_short->Divide(hrtPt_short,hPt_short,1,1.,"b"); | |
571 | hEff_short->SetMaximum(1.4); | |
572 | hEff_short->SetYTitle("Matching Efficiency"); | |
573 | hEff_short->SetXTitle("Pt (GeV/c)"); | |
574 | hEff_short->Draw("same"); | |
575 | ||
576 | c3->cd(2); | |
577 | hrtEta_long->Sumw2(); hEta_long->Sumw2(); | |
578 | hEff_long_eta->Divide(hrtEta_long,hEta_long,1,1.,"b"); | |
579 | hEff_long_eta->SetMaximum(1.4); | |
580 | hEff_long_eta->SetYTitle("Matching Efficiency"); | |
581 | hEff_long_eta->SetXTitle("Eta"); | |
582 | hEff_long_eta->Draw(); | |
583 | hrtEta_short->Sumw2(); hEta_short->Sumw2(); | |
584 | hEff_short_eta->Divide(hrtEta_short,hEta_short,1,1.,"b"); | |
585 | hEff_short_eta->SetMaximum(1.4); | |
586 | hEff_short_eta->SetYTitle("Matching Efficiency"); | |
587 | hEff_short_eta->SetXTitle("Eta"); | |
588 | hEff_short_eta->Draw("same"); | |
589 | c3->Print("c3.ps","ps"); | |
590 | } | |
591 | if(page[4]) { | |
592 | TCanvas *c4=new TCanvas("c4","",0,0,700,850); | |
593 | gStyle->SetOptStat(111110); | |
594 | c4->Divide(2,3); | |
595 | c4->cd(1); hY->Draw(); | |
596 | c4->cd(2); hZ->Draw(); | |
597 | c4->cd(3); hPhi->Draw(); | |
598 | c4->cd(4); hLambda->Draw(); | |
599 | c4->cd(5); hPt->Draw(); | |
600 | c4->Print("c4.ps","ps"); | |
601 | } | |
602 | if(page[5]) { | |
603 | TCanvas *c5=new TCanvas("c5","",0,0,700,850); | |
604 | gStyle->SetOptStat(111110); | |
605 | c5->Divide(2,3); | |
606 | c5->cd(1); hY_n->Draw(); | |
607 | c5->cd(2); hZ_n->Draw(); | |
608 | c5->cd(3); hLambda_n->Draw(); | |
609 | c5->cd(4); hPt_n->Draw(); | |
610 | c5->Print("c5.ps","ps"); | |
611 | } | |
612 | printf("Efficiency (long) = %d/%d = %f\n",nFound_long,nToFind_long, | |
613 | ((Float_t) nFound_long / (Float_t) nToFind_long)); | |
614 | printf("Efficiency (shrt) = %d/%d = %f\n",nFound_short,nToFind_short, | |
615 | ((Float_t) nFound_short / (Float_t) nToFind_short)); | |
616 | } | |
617 | ||
618 |