1 void AliTPCTestTracking() {
2 const char *pname="Param1";
3 const char *tname="TreeD0_Param1";
5 // Dynamically link some shared libs
6 if (gClassTable->GetID("AliRun") < 0) {
7 gROOT->LoadMacro("loadlibs.C");
14 // Connect the Root Galice file containing Geometry, Kine and Hits
15 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
16 if (!file) file = new TFile("galice.root");
18 // Get AliRun object from file or create it if not on file
20 gAlice = (AliRun*)file->Get("gAlice");
21 if (gAlice) printf("AliRun object found on file\n");
22 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
27 TClonesArray *particles=gAlice->Particles();
28 int np=particles->GetEntriesFast();
30 AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
31 int ver=TPC->IsVersion();
32 cerr<<"TPC version "<<ver<<" has been found !\n";
34 AliTPCD *digp= (AliTPCD*)file->Get(pname);
35 if (digp!=0) TPC->SetDigParam(digp);
36 else cerr<<"Warning: default TPC parameters will be used !\n";
38 int nrows=TPC->GetDigParam()->GetParam().GetNRowUp()-1;
41 cerr<<"Making clusters...\n";
43 TClonesArray *clusters=TPC->Clusters();
44 if (!clusters) {cerr<<"No clusters found !\n"; return;}
45 int n=clusters->GetEntriesFast();
46 cerr<<"Number of clusters "<<n<<" \n";
48 cerr<<"Marking \"good\" tracks... \n";
49 for (int i=0; i<n; i++) {
50 AliTPCcluster *c=(AliTPCcluster*)clusters->UncheckedAt(i);
51 int lab=c->fTracks[0];
52 if (lab<0) continue; //noise cluster
54 int sector=c->fSector-1, row=c->fPadRow-1;
55 GParticle *p=(GParticle*)particles->UncheckedAt(lab);
57 if (row==nrows) {ks=p->GetKS()|0x1000; p->SetKS(ks);}
58 if (row==nrows-8) {ks=p->GetKS()|0x800; p->SetKS(ks);}
59 ks=p->GetKS()+1; p->SetKS(ks);
64 cerr<<"Looking for clusters...\n";
65 TPC->Digits2Clusters();
66 TClonesArray *clusters=TPC->Clusters();
67 if (!clusters) {cerr<<"No clusters found !\n"; return;}
68 int n=clusters->GetEntriesFast();
69 cerr<<"Number of clusters "<<n<<" \n";
71 cerr<<"Marking \"good\" tracks... \n";
72 TTree *TD=gDirectory->Get(tname);
73 TClonesArray *digits=TPC->Digits();
74 TD->GetBranch("Digits")->SetAddress(&digits);
76 int *count = new int[np];
78 for (i=0; i<np; i++) count[i]=0;
79 int sectors_by_rows=(int)TD->GetEntries();
80 for (i=0; i<sectors_by_rows; i++) {
81 if (!TD->GetEvent(i)) continue;
83 int ndigits=digits->GetEntriesFast();
85 for (j=0; j<ndigits; j++) {
86 AliTPCdigit *dig = (AliTPCdigit*)digits->UncheckedAt(j);
87 int idx0=dig->fTracks[0];
88 int idx1=dig->fTracks[1];
89 int idx2=dig->fTracks[2];
90 sec=dig->fSector-1; row=dig->fPadRow-1;
91 if (idx0>=0 && dig->fSignal>0) count[idx0]+=1;
92 if (idx1>=0 && dig->fSignal>0) count[idx1]+=1;
93 if (idx2>=0 && dig->fSignal>0) count[idx2]+=1;
95 for (j=0; j<np; j++) {
96 GParticle *p=(GParticle*)particles->UncheckedAt(j);
99 if (row==nrows ) {ks=p->GetKS()|0x1000; p->SetKS(ks);}
100 if (row==nrows-8 ) {ks=p->GetKS()|0x800; p->SetKS(ks);}
101 ks=p->GetKS()+1; p->SetKS(ks);
110 cerr<<"Invalid TPC version !\n";
114 cerr<<"Looking for tracks...\n";
115 TPC->Clusters2Tracks();
117 TClonesArray *tracks=TPC->Tracks();
118 if (tracks) nt=tracks->GetEntriesFast();
119 cerr<<"Number of found tracks "<<nt<<endl;
121 /////////////////////////////////////////////////////////////////////////
122 cerr<<"Doing comparison...\n";
123 TH1F *hp=new TH1F("hp","PHI resolution",50,-100.,100.); hp->SetFillColor(4);
124 TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-100,100); hl->SetFillColor(4);
125 TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
126 hpt->SetFillColor(2);
127 TH1F *hd=new TH1F("hd","Impact parameter distribution ",30,0,25);
130 TH1F *hgood=new TH1F("hgood","Good tracks",20,0,2);
131 TH1F *hfound=new TH1F("hfound","Found tracks",20,0,2);
132 TH1F *hfake=new TH1F("hfake","Fake tracks",20,0,2);
133 TH1F *hg=new TH1F("hg","",20,0,2); //efficiency for good tracks
134 hg->SetLineColor(4); hg->SetLineWidth(2);
135 TH1F *hf=new TH1F("hf","Efficiency for fake tracks",20,0,2);
136 hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
138 for (int i=0; i<np; i++) {
139 GParticle *p = (GParticle*)particles->UncheckedAt(i);
140 if (p->GetParent()>=0) continue; //secondary particle
141 if (p->GetKS()<0x1000+0x800+2+30) continue;
142 Double_t ptg=p->GetPT(),pxg=p->GetPx(),pyg=p->GetPy(),pzg=p->GetPz();
143 if (ptg<0.100) continue;
144 if (fabs(pzg/ptg)>0.999) continue;
150 for (int j=0; j<nt; j++) {
151 AliTPCtrack *track=(AliTPCtrack*)tracks->UncheckedAt(j);
152 int lab=track->GetLab();
153 if (fabs(lab)!=i) continue;
154 //if (lab!=i) continue;
158 track->PropagateTo(xk);
160 track->PropagateTo(xk,42.7,2.27); //C
162 track->PropagateTo(xk,36.2,1.98e-3); //C02
164 track->PropagateTo(xk,42.7,2.27); //C
167 track->PropagateTo(xk,21.82,2.33); //ITS+beam_pipe+etc (approximately)
169 track->PropagateToVertex(); //comparison should be done at the vertex
171 if (lab==i) hfound->Fill(ptg);
172 else { hfake->Fill(ptg); cerr<<lab<<" fake\n";}
173 Double_t px,py,pz,pt=fabs(track->GetPt());track->GetPxPyPz(px,py,pz);
174 Double_t phig=TMath::ATan(pyg/pxg);
175 Double_t phi =TMath::ATan(py /px );
176 hp->Fill((phi - phig)*1000.);
177 Double_t lamg=TMath::ATan(pzg/ptg);
178 Double_t lam =TMath::ATan(pz /pt );
179 hl->Fill((lam - lamg)*1000.);
180 hpt->Fill((pt - ptg)/ptg*100.);
181 Double_t x,y,z; track->GetXYZ(x,y,z);
182 hd->Fill(sqrt(x*x + y*y + z*z)*10.);
185 if (!found) cerr<<"Track number "<<i<<" was not found !\n";
187 Stat_t ngood =hgood->GetEntries(); cerr<<"Good tracks "<<ngood<<endl;
188 Stat_t nfound=hfound->GetEntries();
190 cerr<<"Integral efficiency is about "<<nfound/ngood*100.<<" %\n";
192 gStyle->SetOptStat(111110);
193 gStyle->SetOptFit(1);
195 TCanvas *c1=new TCanvas("c1","",0,0,700,850);
197 TPad *p1=new TPad("p1","",0,0.3,.5,.6); p1->Draw();
198 p1->cd(); p1->SetFillColor(42); p1->SetFrameFillColor(10);
199 hp->SetFillColor(4); hp->SetXTitle("(mrad)"); hp->Fit("gaus"); c1->cd();
201 TPad *p2=new TPad("p2","",0.5,.3,1,.6); p2->Draw();
202 p2->cd(); p2->SetFillColor(42); p2->SetFrameFillColor(10);
203 hl->SetXTitle("(mrad)"); hl->Fit("gaus"); c1->cd();
205 TPad *p3=new TPad("p3","",0,0,0.5,0.3); p3->Draw();
206 p3->cd(); p3->SetFillColor(42); p3->SetFrameFillColor(10);
207 hpt->SetXTitle("(%)"); hpt->Fit("gaus"); c1->cd();
209 TPad *p4=new TPad("p4","",0.5,0,1,0.3); p4->Draw();
210 p4->cd(); p4->SetFillColor(42); p4->SetFrameFillColor(10);
211 hd->SetXTitle("(mm)"); hd->Draw(); c1->cd();
213 TPad *p5=new TPad("p5","",0,0.6,1,1); p5->Draw(); p5->cd();
214 /*p5->SetTopMargin(0.25);*/ p5->SetFillColor(41); p5->SetFrameFillColor(10);
215 hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
216 hg->Divide(hfound,hgood,1,1.,"b");
217 hf->Divide(hfake,hgood,1,1.,"b");
219 hg->SetYTitle("Tracking efficiency");
220 hg->SetXTitle("Pt (GeV/c)");
223 TLine *line1 = new TLine(0,1.0,2,1.0); line1->SetLineStyle(4);
225 TLine *line2 = new TLine(0,0.9,2,0.9); line2->SetLineStyle(4);
229 hf->SetFillStyle(3013);
232 hf->Draw("histsame");
233 TText *text = new TText(0.461176,0.248448,"Fake tracks");
234 text->SetTextSize(0.05);
236 text = new TText(0.453919,1.11408,"Good tracks");
237 text->SetTextSize(0.05);