3 #include "AliTRDtracker.h"
4 #include "AliTRDcluster.h"
6 #include "AliTRDgeometry.h"
7 #include "AliTRDparameter.h"
9 #include "AliTRDmcTrack.h"
10 #include "AliTRDtrack.h"
13 #include "TParticle.h"
14 #include "TStopwatch.h"
17 void AliTRDbackTrackAnalysis() {
19 // const Int_t nPrimaries = 84210/400;
20 const Int_t nPrimaries = 84210/16;
26 TH1F *hp=new TH1F("hp","PHI resolution",100,-20.,20.); hp->SetFillColor(4);
27 TH1F *hl=new TH1F("hl","LAMBDA resolution",100,-100,100);hl->SetFillColor(4);
28 TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
30 TH1F *hpta=new TH1F("hpta","norm. Pt resolution",50,-5.,5.);
31 hpta->SetFillColor(17);
32 hpta->SetXTitle("(%)");
34 TH1F *hla=new TH1F("hla","norm. Lambda resolution",50,-5.,5.);
35 hla->SetFillColor(17);
36 TH1F *hya=new TH1F("hya","norm. Y resolution",50,-5.,5.);
37 hya->SetFillColor(17);
38 TH1F *hza=new TH1F("hza","norm. Z resolution",50,-5.,5.);
39 hza->SetFillColor(17);
43 TH1F *hy=new TH1F("hy","Y resolution",50,-0.5,0.5);hy->SetLineColor(4);
45 hy->SetXTitle("(cm)");
47 TH1F *hz=new TH1F("hz","Z resolution",80,-4,4);hz->SetLineColor(2);
49 hz->SetXTitle("(cm)");
51 TH1F *hx=new TH1F("hx","X(out)",150,250.,400.); hx->SetFillColor(4);
53 const Int_t nPtSteps = 30;
54 const Float_t maxPt = 3.;
56 TH1F *hgood=new TH1F("hgood","long TRD tracks, TPC seeds",nPtSteps,0,maxPt);
57 hgood->SetYTitle("Counts");
58 hgood->SetXTitle("Pt (GeV/c)");
59 hgood->SetLineColor(3);
61 TH1F *hGoodAndSeed = new TH1F("GoodAndSeed","TPC seed and good",nPtSteps,0,maxPt);
62 hGoodAndSeed->SetLineColor(2);
64 TH2F *hRTvsMC = new TH2F("RTvsMC","RTvsMC",100,-4.5,95.5,100,-4.5,95.5);
65 hRTvsMC->SetMarkerColor(4);
66 hRTvsMC->SetMarkerSize(2);
68 TH2F *hXvsMCX = new TH2F("XvsMCX","XvsMCX",150,250,400,150,250,400);
69 hXvsMCX->SetMarkerColor(4);
70 hXvsMCX->SetMarkerSize(2);
72 TH2F *h2seed = new TH2F("2dGood","TPC seeds",60,0.,60,50,0.,3.);
73 h2seed->SetMarkerColor(4);
74 h2seed->SetMarkerSize(2);
76 TH2F *h2lost = new TH2F("2dSeedAndGood","SeedButNotGood",60,0.,60.,50,0.,3.);
77 h2lost->SetMarkerColor(2);
78 h2lost->SetMarkerSize(2);
81 TH1F *hseed=new TH1F("seed","TPC seeds",nPtSteps,0,maxPt);
82 hseed->SetLineColor(4);
85 TH1F *hfound=new TH1F("hfound","Found tracks",nPtSteps,0,maxPt);
86 hfound->SetYTitle("Counts");
87 hfound->SetXTitle("Pt (GeV/c)");
89 TH1F *heff=new TH1F("heff","Matching Efficiency",nPtSteps,0,maxPt); // efficiency, good tracks
90 heff->SetLineColor(4); heff->SetLineWidth(2);
91 heff->SetXTitle("Pt, GeV/c");
92 heff->SetYTitle("Efficiency");
94 TH1F *hSeedEff = new TH1F("hSeedEff","TPC Efficiency",nPtSteps,0,maxPt);
95 hSeedEff->SetLineColor(4); hSeedEff->SetLineWidth(2);
96 hSeedEff->SetXTitle("Pt, GeV/c");
97 hSeedEff->SetYTitle("Efficiency");
100 TH1F *hol=new TH1F("hol","Overlap fraction",105,-2.5,102.5);
101 hol->SetLineColor(4); hol->SetLineWidth(2);
102 hol->SetXTitle("Fraction,(%)");
103 hol->SetYTitle("Counts");
105 TH1F *hend=new TH1F("end","missing tail",80,-10.5,69.5);
107 TH1F *hFraction=new TH1F("fraction","Fraction of found clusters",110,0,1.1);
108 TH1F *hCorrect=new TH1F("correct","Fraction of correct clusters",110,0,1.1);
112 const Int_t maxIndex = nPrimaries;
113 Bool_t seedLabel[maxIndex];
114 Int_t mcIndex[maxIndex];
115 Int_t rtIndex[maxIndex];
117 for(Int_t i = 0; i < maxIndex; i++) {
118 seedLabel[i] = kFALSE;
123 // mark available seeds from TPC
125 Int_t nPrimarySeeds = 0;
127 printf("marking found seeds from TPC\n");
128 TDirectory *savedir=gDirectory;
130 TFile *in=TFile::Open("AliTPCBackTracks.root");
132 cerr<<"can't open file AliTPCBackTracks.root !\n"; return;
136 sprintf(tname,"seedsTPCtoTRD_%d",nEvent);
137 TTree *seedTree=(TTree*)in->Get(tname);
139 cerr<<"AliTRDtracker::PropagateBack(): ";
140 cerr<<"can't get a tree with seeds from TPC !\n";
143 AliTPCtrack *seed=new AliTPCtrack;
144 seedTree->SetBranchAddress("tracks",&seed);
146 Int_t n=(Int_t)seedTree->GetEntries();
147 for (Int_t i=0; i<n; i++) {
148 seedTree->GetEvent(i);
149 Int_t lbl = seed->GetLabel();
151 printf("negative seed label %d \n",lbl);
154 if(lbl >= maxIndex) continue;
155 seedLabel[lbl] = kTRUE;
156 if(lbl < nPrimaries) nPrimarySeeds++;
162 printf("Found %d seeds from primaries among overall %d seeds \n",
163 nPrimarySeeds, nSeeds);
166 // done with marking TPC seeds
169 TFile *tf=TFile::Open("AliTRDtracks.root");
171 if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return;}
172 TObjArray tarray(2000);
174 sprintf(tname,"TRDb_%d",nEvent);
175 TTree *tracktree=(TTree*)tf->Get(tname);
177 TBranch *tbranch=tracktree->GetBranch("tracks");
179 Int_t nRecTracks = (Int_t) tracktree->GetEntries();
180 cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;
182 for (Int_t i=0; i<nRecTracks; i++) {
183 AliTRDtrack *iotrack=new AliTRDtrack();
184 tbranch->SetAddress(&iotrack);
185 tracktree->GetEvent(i);
186 tarray.AddLast(iotrack);
187 Int_t trackLabel = iotrack->GetLabel();
189 // printf("rt with %d clusters and label %d \n",
190 // iotrack->GetNumberOfClusters(), trackLabel);
192 if(trackLabel < 0) continue;
193 if(trackLabel >= maxIndex) continue;
194 rtIndex[trackLabel] = i;
201 TFile *mctf=TFile::Open("AliTRDmcTracks.root");
202 if (!mctf->IsOpen()) {cerr<<"Can't open AliTRDmcTracks.root !\n"; return;}
203 TObjArray mctarray(2000);
204 TTree *mctracktree=(TTree*)mctf->Get("MCtracks");
205 TBranch *mctbranch=mctracktree->GetBranch("MCtracks");
206 Int_t nMCtracks = (Int_t) mctracktree->GetEntries();
207 cerr<<"Found "<<nMCtracks<<" entries in the MC tracks tree"<<endl;
208 for (Int_t i=0; i<nMCtracks; i++) {
209 AliTRDmcTrack *ioMCtrack=new AliTRDmcTrack;
210 mctbranch->SetAddress(&ioMCtrack);
211 mctracktree->GetEvent(i);
212 mctarray.AddLast(ioMCtrack);
213 Int_t mcLabel = ioMCtrack->GetTrackIndex();
214 if(mcLabel < 0) {printf("negative mc label detected!\n"); continue;}
215 if(mcLabel >= maxIndex) continue;
216 mcIndex[mcLabel] = i;
223 TFile *geofile =TFile::Open("AliTRDclusters.root");
224 AliTRDtracker *Tracker = new AliTRDtracker(geofile);
225 Tracker->SetEventNumber(nEvent);
227 AliTRDgeometry *fGeom = (AliTRDgeometry*) geofile->Get("TRDgeometry");
228 AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter");
230 Char_t *alifile = "AliTRDclusters.root";
231 TObjArray carray(2000);
232 TObjArray *ClustersArray = &carray;
233 Tracker->ReadClusters(ClustersArray,alifile);
236 // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
237 alifile = "galice.root";
239 TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
241 cout << "Open the ALIROOT-file " << alifile << endl;
242 gafl = new TFile(alifile);
245 cout << alifile << " is already open" << endl;
248 // Get AliRun object from file or create it if not on file
249 gAlice = (AliRun*) gafl->Get("gAlice");
251 cout << "AliRun object found on file" << endl;
253 gAlice = new AliRun("gAlice","Alice test program");
256 // Define the objects
257 AliTRDv1 *TRD = (AliTRDv1*) gAlice->GetDetector("TRD");
259 // Import the Trees for the event nEvent in the file
260 const Int_t nparticles = gAlice->GetEvent(nEvent);
261 if (nparticles <= 0) return;
264 Bool_t electrons[300000] = { kFALSE };
266 Bool_t mark_electrons = kFALSE;
268 printf("mark electrons\n");
270 for(Int_t i = nPrimaries; i < nparticles; i++) {
271 p = gAlice->Particle(i);
272 if(p->GetMass() > 0.01) continue;
273 if(p->GetMass() < 0.00001) continue;
274 electrons[i] = kTRUE;
278 AliTRDcluster *cl = 0;
279 Int_t nw, label, index, ti[3], tbwc;
280 Int_t det, plane, ltb, gtb, gap, max_gap, sector, mc_sector, min_tb, max_tb;
281 Double_t Pt, Px, Py, Pz;
282 Double_t mcPt, mcPx, mcPy, mcPz;
284 Int_t rtClusters, rtCorrect;
288 AliTRDmcTrack *mct = 0;
291 Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
292 Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
293 Double_t dx = (Double_t) fPar->GetTimeBinSize();
295 Int_t tbAmp = fPar->GetTimeBefore();
296 Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
297 Int_t tbDrift = fPar->GetTimeMax();
298 Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
300 tbDrift = TMath::Min(tbDrift,maxDrift);
301 tbAmp = TMath::Min(tbAmp,maxAmp);
303 const Int_t nPlanes = fGeom->Nplan();
304 const Int_t tbpp = Tracker->GetTimeBinsPerPlane();
305 const Int_t nTB = tbpp * nPlanes;
308 for(Int_t i=0; i < maxIndex; i++) {
311 if(rtIndex[i] >= 0) rt = (AliTRDtrack *) tarray.UncheckedAt(rtIndex[i]);
312 if(mcIndex[i] >= 0) mct = (AliTRDmcTrack *) mctarray.UncheckedAt(mcIndex[i]);
315 label = mct->GetTrackIndex();
317 // if(TMath::Abs(mct->GetMass()-0.136) > 0.01) continue;
319 Int_t ncl = mct->GetNumberOfClusters();
321 // check how many time bins have a cluster
323 for(nw = 0; nw < nTB; nw++) mask[nw] = -1;
325 for(Int_t j = 0; j < ncl; j++) {
326 index = mct->GetClusterIndex(j);
327 cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
329 for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
331 if((ti[0] != label) && (ti[1] != label) && (ti[2] != label)) {
332 printf("wrong track label: %d, %d, %d != %d \n",
333 ti[0], ti[1], ti[2], label);
335 det=cl->GetDetector();
336 plane = fGeom->GetPlane(det);
337 ltb = cl->GetLocalTimeBin();
338 gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
345 for(plane = 0; plane < nPlanes; plane++) {
346 for(ltb = tbDrift-1; ltb >= -tbAmp; ltb--) {
347 gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
348 if(mask[gtb] > -1) break;
350 if(ltb >= -tbAmp) break;
352 if((plane == nPlanes) && (ltb == -tbAmp-1)) {
353 printf("warning: for track %d min tb is not found and set to %d!\n",
356 // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
359 min_tb = Tracker->GetGlobalTimeBin(0,plane,ltb);
362 for(plane = nPlanes-1 ; plane>=0; plane--) {
363 for(ltb = -tbAmp; ltb < tbDrift; ltb++) {
364 gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
365 if(mask[gtb] > -1) break;
367 if(ltb < tbDrift) break;
369 if((plane == -1) && (ltb == tbDrift)) {
370 printf("warning: for track %d max tb is not found and set to 0!\n",label);
371 // for(Int_t tb = 0; tb<nTB; tb++) printf("gtb %d; cl index %d\n",tb,mask[tb]);
373 mcX = Tracker->GetX(0,0,tbDrift-1);
376 max_tb = Tracker-> GetGlobalTimeBin(0,plane,ltb);
377 mcX = Tracker->GetX(0,plane,ltb);
378 cl = (AliTRDcluster *) ClustersArray->UncheckedAt(mask[gtb]);
379 det=cl->GetDetector();
380 sector = fGeom->GetSector(det);
389 for(nw = min_tb; nw < max_tb+1; nw++) {
393 if(gap > max_gap) max_gap = gap;
398 if(tbwc < ((Int_t) (nTB * Tracker->GetMinClustersInTrack()))) continue;
399 if(gap > Tracker->GetMaxGap()) continue;
400 p = gAlice->Particle(label);
402 printf("good track %d with min_tb %d, max_tb %d, gap %d\n",
403 label,min_tb,max_tb,gap);
405 hgood->Fill(p->Pt());
409 if(rt->GetLabel() != label) {
410 printf("mct vs rt index mismatch: %d != %d \n",
411 label, rt->GetLabel());
415 if(!seedLabel[i]) continue;
417 hGoodAndSeed->Fill(p->Pt());
419 hXvsMCX->Fill(mcX,rt->GetX());
421 rtClusters = rt->GetNumberOfClusters();
423 // find number of tb with correct cluster
426 for(Int_t j = 0; j < rtClusters; j++) {
427 index = rt->GetClusterIndex(j);
428 cl = (AliTRDcluster *) ClustersArray->UncheckedAt(index);
430 for(nw = 0; nw < 3; nw++) ti[nw] = cl->GetLabel(nw);
432 if((ti[0] != label)&&(ti[1] != label)&&(ti[2] != label)) continue;
436 Float_t foundFraction = ((Float_t) rtCorrect) / (tbwc + 0.00001);
437 Float_t correctFraction = ((Float_t) rtCorrect) / ((Float_t) rtClusters);
439 if(foundFraction > 1) printf("fraction = %f/%f \n",
440 (Float_t) rtCorrect, (Float_t) tbwc);
443 if(foundFraction <= 1) hFraction->Fill(foundFraction);
444 hCorrect->Fill(correctFraction);
446 hRTvsMC->Fill((Float_t) tbwc, (Float_t) rtClusters);
448 if((foundFraction < 0.7) || (correctFraction < 0.7)) {
449 printf("not found track %d with FrctnFound %f and FrctnCorrect %f\n",
450 label, foundFraction, correctFraction);
454 hfound->Fill(p->Pt());
456 Pt = TMath::Abs(rt->GetPt());
457 Double_t cc = TMath::Abs(rt->GetSigmaC2());
458 mct->GetPxPyPzXYZ(mcPx,mcPy,mcPz,x,y,z,-1);
459 rt->GetPxPyPz(Px,Py,Pz);
461 printf("\n\ntrack %d \n", label);
462 printf("rt Px, Py, Pz: %f, %f, %f \n", Px, Py, Pz);
463 printf("mc Px, Py, Pz: %f, %f, %f \n", mcPx, mcPy, mcPz);
465 mcPt = TMath::Sqrt(mcPx * mcPx + mcPy * mcPy);
466 if(mcPt < 0.0001) mcPt = 0.0001;
468 Float_t lamg=TMath::ATan2(mcPz,mcPt);
469 Float_t lam =TMath::ATan2(Pz,Pt);
470 if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
471 else hpta->Fill((0.3*0.4/100*(1/Pt - 1/mcPt))/TMath::Sqrt(cc));
473 Float_t phig=TMath::ATan2(mcPy,mcPx);
474 Float_t phi =TMath::ATan2(Py,Px);
477 if(!(rt->PropagateTo(x))) continue;
480 if((mcPt > Pt_min) && (mcPt < Pt_max)) {
481 hl->Fill((lam - lamg)*1000.);
482 hla->Fill((lam - lamg)/TMath::Sqrt(rt->GetSigmaTgl2()));
483 if(TMath::Abs(mcPt) < 0.0001) printf("attempt to divide by mcPt = %f\n",mcPt);
484 else hpt->Fill((1/Pt - 1/mcPt)/(1/mcPt)*100.);
485 hp->Fill((phi - phig)*1000.);
486 hy->Fill(rt->GetY() - y);
487 hz->Fill(rt->GetZ() - z);
488 hya->Fill((rt->GetY() - y)/TMath::Sqrt(rt->GetSigmaY2()));
489 hza->Fill((rt->GetZ() - z)/TMath::Sqrt(rt->GetSigmaZ2()));
490 hx->Fill((Float_t) x);
495 gStyle->SetOptStat(0);
496 // gStyle->SetOptStat(111110);
497 gStyle->SetOptFit(1);
501 TCanvas *c1=new TCanvas("c1","",0,0,700,850);
503 // gStyle->SetOptStat(0);
504 // gStyle->SetOptStat(111110);
505 gStyle->SetOptFit(1);
507 TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
508 p1->cd(); p1->SetFillColor(10); p1->SetFrameFillColor(10);
509 hgood->Draw(); hGoodAndSeed->Draw("same"); c1->cd();
511 TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
512 p2->cd(); p2->SetFillColor(10); p2->SetFrameFillColor(10);
513 hFraction->SetYTitle("Counts");
514 hFraction->SetXTitle("Fraction");
515 hFraction->Draw(); c1->cd();
517 TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
518 p3->cd(); p3->SetFillColor(10); p3->SetFrameFillColor(10);
519 hfound->Draw(); c1->cd();
521 TPad *p4=new TPad("p4","",0.5,0,1.,0.3); p4->Draw();
522 p4->cd(); p4->SetFillColor(10); p4->SetFrameFillColor(10);
523 hCorrect->SetYTitle("Counts");
524 hCorrect->SetXTitle("Fraction");
525 hCorrect->Draw(); c1->cd();
527 TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
528 p5->SetFillColor(10); p5->SetFrameFillColor(10);
529 hgood->Sumw2(); hGoodAndSeed->Sumw2(); // hfake->Sumw2();
530 hSeedEff->Divide(hGoodAndSeed,hgood,1,1.,"b");
531 // hf->Divide(hfake,hgood,1,1.,"b");
532 hSeedEff->SetMaximum(1.4);
533 hSeedEff->SetYTitle("TPC efficiency");
534 hSeedEff->SetXTitle("Pt (GeV/c)");
537 TLine *line1 = new TLine(0,1.0,maxPt,1.0); line1->SetLineStyle(4);
539 TLine *line2 = new TLine(0,0.9,maxPt,0.9); line2->SetLineStyle(4);
543 TCanvas *c2=new TCanvas("c2","",0,0,700,850);
545 TPad *p12=new TPad("p12","",0,0.3,.5,.6); p12->Draw();
546 p12->cd(); p12->SetFillColor(42); p12->SetFrameFillColor(10);
547 hp->SetFillColor(4); hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c2->cd();
549 TPad *p22=new TPad("p22","",0.5,.3,1,.6); p22->Draw();
550 p22->cd(); p22->SetFillColor(42); p22->SetFrameFillColor(10);
551 hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c2->cd();
553 TPad *p32=new TPad("p32","",0,0,0.5,0.3); p32->Draw();
554 p32->cd(); p32->SetFillColor(42); p32->SetFrameFillColor(10);
555 hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c2->cd();
557 TPad *p42=new TPad("p42","",0.5,0,1.,0.3); p42->Draw();
558 p42->cd(); p42->SetFillColor(42); p42->SetFrameFillColor(10);
559 hgood->Draw(); hGoodAndSeed->Draw("same"); c2->cd();
561 TPad *p52=new TPad("p52","",0,0.6,1,1); p52->Draw(); p52->cd();
562 p52->SetFillColor(41); p52->SetFrameFillColor(10);
564 heff->Divide(hfound,hGoodAndSeed,1,1.,"b");
565 // hf->Divide(hfake,hgood,1,1.,"b");
566 heff->SetMaximum(1.4);
567 heff->SetXTitle("Pt (GeV/c)");
570 TLine *line12 = new TLine(0,1.0,maxPt,1.0); line12->SetLineStyle(4);
571 line12->Draw("same");
572 TLine *line22 = new TLine(0,0.9,maxPt,0.9); line22->SetLineStyle(4);
573 line22->Draw("same");
575 c2->Print("matching.ps","ps");
579 TCanvas *cxyz = new TCanvas("cxyz","",50,50,750,900);
581 cxyz->cd(1); hx->Draw();
582 cxyz->cd(2); hy->Draw();
583 cxyz->cd(3); hz->Draw();
584 cxyz->cd(4); hXvsMCX->Draw();
587 TCanvas *cs = new TCanvas("cs","",0,0,700,850);
590 cs->cd(1); hy->Fit("gaus");
591 cs->cd(2); hz->Fit("gaus");
592 cs->cd(3); hpta->Fit("gaus");
593 cs->cd(4); hla->Fit("gaus");
594 cs->cd(5); hya->Fit("gaus");
595 cs->cd(6); hza->Fit("gaus");
597 cs->Print("resolution.ps","ps");
600 TCanvas *cvs = new TCanvas("cvs","",0,0,700,850);
601 cvs->cd(); hRTvsMC->Draw();