Fixing coding convention violations
[u/mrichter/AliRoot.git] / TOF / AliTOFComparison.C
CommitLineData
3f83f224 1/****************************************************************************
2 * This macro estimates efficiency of matching with the TOF. *
3 * TOF "Good" tracks are those originating from the primary vertex, *
4 * being "good" in the ITS and having at least one digit in the TOF. *
5 * (To get the list of "good" tracks one should first run *
6 * AliTPCComparison.C and AliITSComparisonV2.C macros) *
7 * Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch *
8 ****************************************************************************/
9
10#if !defined(__CINT__) || defined(__MAKECINT__)
11 #include <Riostream.h>
12 #include <fstream.h>
13
14 #include "AliRun.h"
15 #include "AliHeader.h"
16 #include "AliStack.h"
17 #include "AliRunLoader.h"
18 #include "AliLoader.h"
19
20 #include "AliESD.h"
21 #include "AliTOFdigit.h"
22
23 #include "TKey.h"
24 #include "TFile.h"
25 #include "TTree.h"
26 #include "TH1.h"
27 #include "TClonesArray.h"
28 #include "TStyle.h"
29 #include "TCanvas.h"
30 #include "TLine.h"
31 #include "TText.h"
32 #include "TParticle.h"
33#endif
34
35struct GoodTrackTOF {
36 Int_t lab;
37 Int_t code;
38 Float_t px,py,pz;
39 Float_t x,y,z;
40};
41
42extern AliRun *gAlice;
43
44Int_t AliTOFComparison() {
45 Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max);
46
47 cerr<<"Doing comparison...\n";
48
49 if (gAlice) {
50 delete gAlice->GetRunLoader();
51 delete gAlice;//if everything was OK here it is already NULL
52 gAlice = 0x0;
53 }
54
55 TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
56
57 {
58 AliRunLoader *rl = AliRunLoader::Open("galice.root");
59 if (rl == 0x0) {
60 cerr<<"Can not open session"<<endl;
61 return 1;
62 }
63 AliLoader* tofl = rl->GetLoader("TOFLoader");
64 if (tofl == 0x0) {
65 cerr<<"Can not get the TOF loader"<<endl;
66 return 2;
67 }
68 tofl->LoadDigits("read");
69
70 rl->GetEvent(0);
71
72 TTree *tofTree=tofl->TreeD();
73 if (!tofTree) {
74 cerr<<"Can't get the TOF cluster tree !\n";
75 return 3;
76 }
77 TBranch *branch=tofTree->GetBranch("TOF");
78 if (!branch) {
79 cerr<<"Can't get the branch with the TOF digits !\n";
80 return 4;
81 }
82 branch->SetAddress(&digits);
83
84 tofTree->GetEvent(0);
85
86 delete rl;
87 }
88
89 Int_t nd=digits->GetEntriesFast();
90 cerr<<"Number of digits: "<<nd<<endl;
91
92 const Int_t MAX=15000;
93 GoodTrackTOF gt[MAX];
94 Int_t ngood=0;
95 ifstream in("good_tracks_tof");
96 if (in) {
97 cerr<<"Reading good tracks...\n";
98 while (in>>gt[ngood].lab>>gt[ngood].code>>
99 gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
100 gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
101 ngood++;
102 cerr<<ngood<<'\r';
103 if (ngood==MAX) {
104 cerr<<"Too many good tracks !\n";
105 break;
106 }
107 }
108 if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
109 } else {
110 cerr<<"Marking good tracks (this will take a while)...\n";
111 ngood=good_tracks_tof(gt,MAX);
112 ofstream out("good_tracks_tof");
113 if (out) {
114 for (Int_t ngd=0; ngd<ngood; ngd++)
115 out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
116 <<gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
117 <<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
118 } else cerr<<"Can not open file (good_tracks_tof) !\n";
119 out.close();
120 }
121
122 Double_t pmin=0.2;
123 Double_t pmax=4.0;
124
125 TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
126 TH1F *hfound=new TH1F("hfound","Matched tracks",30,pmin,pmax);
127 TH1F *hfake=new TH1F("hfake","Mismatched tracks",30,pmin,pmax);
128 TH1F *hgp=new TH1F("hgp","",30,pmin,pmax); //efficiency for good tracks
129 hgp->SetLineColor(4); hgp->SetLineWidth(2);
130 TH1F *hfp=new TH1F("hfp","Probability of mismatching",30,pmin,pmax);
131 hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
132
133 TH1F *hgoo=new TH1F("hgoo","Good tracks",30,-1,1);
134 TH1F *hfoun=new TH1F("hfoun","Matched tracks",30,-1,1);
135 TH1F *hfak=new TH1F("hfak","Mismatched tracks",30,-1,1);
136 TH1F *hgl=new TH1F("hgl","",30,-1,1); //efficiency for good tracks
137 hgl->SetLineColor(4); hgl->SetLineWidth(2);
138 TH1F *hfl=new TH1F("hfl","Probability of mismatching",30,-1,1);
139 hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
140
1c3e95d8 141 TCanvas *c1=new TCanvas("c1","",0,0,600,900);
142 c1->Divide(1,2);
143
3f83f224 144 TFile *ef=TFile::Open("AliESDs.root");
145 if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
146
147 TIter next(ef->GetListOfKeys());
148 TKey *key=0;
149 Int_t nev=0;
150 while ((key=(TKey*)next())!=0) {
151 cerr<<"Processing event number : "<<nev++<<endl;
152
153 AliESD *event=(AliESD*)key->ReadObj();
154 Int_t ntrk=event->GetNumberOfTracks();
155 cerr<<"Number of ESD tracks : "<<ntrk<<endl;
156
157 Int_t matched=0;
158 Int_t mismatched=0;
159 for (Int_t i=0; i<ngood; i++) {
160 Int_t lab=gt[i].lab;
161 Double_t pxg=gt[i].px, pyg=gt[i].py, pzg=gt[i].pz;
162 Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
163
164 if (ptg<0.1) continue;
165
166 Double_t tgl=pzg/ptg; //tan(lambda)
167
168 if (ptg>pmin) { hgood->Fill(ptg); hgoo->Fill(tgl); }
169
170 Int_t j;
171 AliESDtrack *t=0;
172 for (j=0; j<ntrk; j++) {
173 AliESDtrack *tt=event->GetTrack(j);
174 if (lab!=TMath::Abs(tt->GetLabel())) continue;
175 t=tt;
176 //if ((tt->GetStatus()&AliESDtrack::kTOFpid) == 0) continue;
177 if (tt->GetTOFsignal() < 0) continue;
178 UInt_t idx=tt->GetTOFcluster();
179 if ((Int_t)idx>=nd) {
180 cerr<<"Wrong digit index ! "<<idx<<endl;
181 return 5;
182 }
183 AliTOFdigit *dig=(AliTOFdigit*)digits->UncheckedAt(idx);
184 Int_t *label=dig->GetTracks();
185 if (label[0]!=lab)
186 if (label[1]!=lab)
187 if (label[2]!=lab) {
188 mismatched++;
189 if (ptg>pmin) { hfake->Fill(ptg); hfak->Fill(tgl); }
190 break;
191 }
192 if (ptg>pmin) { hfound->Fill(ptg); hfoun->Fill(tgl); }
193 matched++;
194 break;
195 }
196 if (j==ntrk) {
197 cerr<<"Not matched: "<<lab<<" ";
198 if (t) {
199 cerr<<(t->GetStatus()&AliESDtrack::kITSout)<<' '
200 <<(t->GetStatus()&AliESDtrack::kTPCout)<<' '
201 <<(t->GetStatus()&AliESDtrack::kTRDout)<<' '
202 <<(t->GetStatus()&AliESDtrack::kTIME);
203 } else cerr<<"No ESD track !";
204 cerr<<endl;
205 }
206 }
207
208 cerr<<"Number of good tracks: "<<ngood<<endl;
209 cerr<<"Number of matched tracks: "<<matched<<endl;
210 cerr<<"Number of mismatched tracks: "<<mismatched<<endl;
211 if (ngood!=0) cerr<<"Efficiency: "<<Float_t(matched)/ngood<<endl;
212
213 hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
214 hgp->Divide(hfound,hgood,1,1.,"b");
215 hfp->Divide(hfake,hgood,1,1.,"b");
216 hgp->SetMaximum(1.4);
217 hgp->SetYTitle("Matching efficiency");
218 hgp->SetXTitle("Pt (GeV/c)");
219
220 hfoun->Sumw2(); hgoo->Sumw2(); hfak->Sumw2();
221 hgl->Divide(hfoun,hgoo,1,1.,"b");
222 hfl->Divide(hfak,hgoo,1,1.,"b");
223 hgl->SetMaximum(1.4);
224 hgl->SetYTitle("Matching efficiency");
225 hgl->SetXTitle("Tan(lambda)");
226
3f83f224 227 c1->cd(1);
228
229 hgp->Draw();
230 hfp->Draw("histsame");
231 TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
232 line1->Draw("same");
233 TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
234 line2->Draw("same");
235
236 c1->cd(2);
237
238 hgl->Draw();
239 hfl->Draw("histsame");
240 TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
241 line3->Draw("same");
242 TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
243 line4->Draw("same");
244
245 break;
246 }
247
3f7a17bb 248 TFile fc("AliTOFComparison.root","RECREATE");
249 c1->Write();
250 fc.Close();
251
3f83f224 252 return 0;
253}
254
255Int_t good_tracks_tof(GoodTrackTOF *gt, const Int_t max) {
256 ifstream in("good_tracks_its");
257 if (!in) {
258 cerr<<"Can't get good_tracks_its !\n"; exit(11);
259 }
260
261 AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON");
262 if (!rl) {
263 cerr<<"Can't start session !\n";
264 exit(1);
265 }
266
267 rl->GetEvent(0);
268
269
270 rl->LoadgAlice();
271 rl->LoadHeader();
272 Int_t np = rl->GetHeader()->GetNtrack();
273
274 Int_t *good=new Int_t[np];
275 Int_t k;
276 for (k=0; k<np; k++) good[k]=0;
277
278 AliLoader* tofl = rl->GetLoader("TOFLoader");
279 if (tofl == 0x0) {
280 cerr<<"Can not get the TOF loader"<<endl;
281 exit(2);
282 }
283 tofl->LoadDigits("read");
284
285 TTree *dTree=tofl->TreeD();
286 if (!dTree) {
287 cerr<<"Can't get the TOF cluster tree !\n";
288 exit(3);
289 }
290
291 TBranch *branch=dTree->GetBranch("TOF");
292 if (!branch) {
293 cerr<<"Can't get the branch with the TOF digits !\n";
294 exit(4);
295 }
296 TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
297 branch->SetAddress(&digits);
298
299 dTree->GetEvent(0);
300 Int_t nd=digits->GetEntriesFast();
301 cerr<<"Number of digits: "<<nd<<endl;
302
303 for (Int_t i=0; i<nd; i++) {
304 AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
305 Int_t l0=d->GetTrack(0);
306 if (l0>=np) {cerr<<"Wrong label: "<<l0<<endl; continue;}
307 Int_t l1=d->GetTrack(1);
308 if (l1>=np) {cerr<<"Wrong label: "<<l1<<endl; continue;}
309 Int_t l2=d->GetTrack(2);
310 if (l2>=np) {cerr<<"Wrong label: "<<l2<<endl; continue;}
311 if (l0>=0) good[l0]++;
312 if (l1>=0) good[l1]++;
313 if (l2>=0) good[l2]++;
314 }
315
316
317 rl->LoadKinematics();
318 AliStack* stack = rl->Stack();
319 Int_t nt=0;
320 Double_t px,py,pz,x,y,z;
321 Int_t code,lab;
322 while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
323 if (good[lab] == 0) continue;
324 TParticle *p = (TParticle*)stack->Particle(lab);
325 if (p == 0x0) {
326 cerr<<"Can not get particle "<<lab<<endl;
327 exit(1);
328 }
329 if (TMath::Abs(p->Vx())>0.1) continue;
330 if (TMath::Abs(p->Vy())>0.1) continue;
2b85d556 331 // if (TMath::Abs(p->Vz())>0.1) continue;
3f83f224 332
333 gt[nt].lab=lab;
334 gt[nt].code=p->GetPdgCode();
335//**** px py pz - in global coordinate system
336 gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
337 gt[nt].x=p->Vx(); gt[nt].y=p->Vy(); gt[nt].z=p->Vz();
338 nt++;
339 if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
340 }
341
342 delete[] good;
343
344 rl->UnloadKinematics();
345 rl->UnloadHeader();
346 rl->UnloadgAlice();
347 delete rl;
348
349 return nt;
350}