5246387983ef05f4472f93c344edcf4f83d21e75
[u/mrichter/AliRoot.git] / TOF / AliTOFComparison.C
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
35 struct GoodTrackTOF {
36   Int_t lab;
37   Int_t code;
38   Float_t px,py,pz;
39   Float_t x,y,z;
40 };
41
42 extern AliRun *gAlice;
43
44 Int_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
141    TCanvas *c1=new TCanvas("c1","",0,0,600,900);
142    c1->Divide(1,2);
143
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
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
248    TFile fc("AliTOFComparison.root","RECREATE");
249    c1->Write();
250    fc.Close();
251
252    return 0;
253 }
254
255 Int_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;
331       //      if (TMath::Abs(p->Vz())>0.1) continue;
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 }