Possiblity to run over several events in a file (Yu.Belikov)
[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 <TMath.h>
12   #include <TError.h>
13   #include <Riostream.h>
14   #include <TH1F.h>
15   #include <TTree.h>
16   #include <TParticle.h>
17   #include <TCanvas.h>
18   #include <TLine.h>
19   #include <TText.h>
20   #include <TBenchmark.h>
21   #include <TStyle.h>
22   #include <TKey.h>
23   #include <TROOT.h>
24
25   #include "AliStack.h"
26   #include "AliHeader.h"
27   #include "AliTrackReference.h"
28   #include "AliRunLoader.h"
29   #include "AliRun.h"
30   #include "AliESD.h"
31
32   #include "AliTOFdigit.h"
33   #include "AliLoader.h"
34
35   #include "TClonesArray.h"
36 #endif
37
38 Int_t GoodTracksTOF(const Char_t *dir=".");
39
40 extern AliRun *gAlice;
41 extern TBenchmark *gBenchmark;
42 extern TROOT *gROOT;
43
44 static Int_t allgood=0;
45 static Int_t allmatched=0;
46 static Int_t allmismatched=0;
47
48 Int_t AliTOFComparison(const Char_t *dir=".") {
49    gBenchmark->Start("AliTOFComparison");
50
51    ::Info("AliTOFComparison.C","Doing comparison...");
52    
53
54    Double_t pmin=0.2;
55    Double_t pmax=3.0;
56
57    TH1F *hgood=(TH1F*)gROOT->FindObject("hgood");    
58    if (!hgood) hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
59     
60    TH1F *hfound=(TH1F*)gROOT->FindObject("hfound");
61    if (!hfound) hfound=new TH1F("hfound","Matched tracks",30,pmin,pmax);
62
63    TH1F *hfake=(TH1F*)gROOT->FindObject("hfake");
64    if (!hfake) hfake=new TH1F("hfake","Mismatched tracks",30,pmin,pmax);
65
66    TH1F *hgp=(TH1F*)gROOT->FindObject("hgp");
67    if (!hgp) hgp=new TH1F("hgp","",30,pmin,pmax);
68    hgp->SetLineColor(4); hgp->SetLineWidth(2);
69
70    TH1F *hfp=(TH1F*)gROOT->FindObject("hfp");
71    if (!hfp) hfp=new TH1F("hfp","Probability of mismatching",30,pmin,pmax);
72    hfp->SetFillColor(1); hfp->SetFillStyle(3013); hfp->SetLineWidth(2);
73
74    TH1F *hgoo=(TH1F*)gROOT->FindObject("hgoo");    
75    if (!hgoo) hgoo=new TH1F("hgoo","Good tracks",30,-1,1);
76     
77    TH1F *hfoun=(TH1F*)gROOT->FindObject("hfoun");
78    if (!hfoun) hfoun=new TH1F("hfoun","Matched tracks",30,-1,1);
79
80    TH1F *hfak=(TH1F*)gROOT->FindObject("hfak");
81    if (!hfak) hfak=new TH1F("hfak","Mismatched tracks",30,-1,1);
82
83    TH1F *hgl=(TH1F*)gROOT->FindObject("hgl");
84    if (!hgl) hgl=new TH1F("hgl","",30,-1,1);
85    hgl->SetLineColor(4); hgl->SetLineWidth(2);
86
87    TH1F *hfl=(TH1F*)gROOT->FindObject("hfl");
88    if (!hfl) hfl=new TH1F("hfl","Probability of mismatching",30,-1,1);
89    hfl->SetFillColor(1); hfl->SetFillStyle(3013); hfl->SetLineWidth(2);
90
91
92
93    Char_t fname[100];
94    sprintf(fname,"%s/GoodTracksTOF.root",dir);
95
96    TFile *refFile=TFile::Open(fname,"old");
97    if (!refFile || !refFile->IsOpen()) {
98    ::Info("AliTOFComparison.C","Marking good tracks (will take a while)...");
99      if (GoodTracksTOF(dir)) {
100         ::Error("AliTOFComparison.C","Can't generate the reference file !");
101         return 1;
102      }
103    }
104    refFile=TFile::Open(fname,"old");
105    if (!refFile || !refFile->IsOpen()) {
106      ::Error("AliTOFComparison.C","Can't open the reference file !");
107      return 1;
108    }   
109
110    TTree *tofTree=(TTree*)refFile->Get("tofTree");
111    if (!tofTree) {
112      ::Error("AliTOFComparison.C","Can't get the reference tree !");
113      return 2;
114    }
115    TBranch *branch=tofTree->GetBranch("TOF");
116    if (!branch) {
117      ::Error("AliTOFComparison.C","Can't get the TOF branch !");
118      return 3;
119    }
120    TClonesArray dummy("AliTrackReference",1000), *refs=&dummy;
121    branch->SetAddress(&refs);
122
123
124    if (gAlice) { 
125      delete gAlice->GetRunLoader();
126      delete gAlice;//if everything was OK here it is already NULL
127      gAlice = 0x0;
128    }
129    sprintf(fname,"%s/galice.root",dir);
130    AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
131    if (rl == 0x0) {
132       cerr<<"Can not open session"<<endl;
133       return 1;
134    }
135    AliLoader* tofl = rl->GetLoader("TOFLoader");
136    if (tofl == 0x0) {
137       cerr<<"Can not get the TOF loader"<<endl;
138       return 2;
139    }
140    tofl->LoadDigits("read");
141
142
143    sprintf(fname,"%s/AliESDs.root",dir);
144    TFile *ef=TFile::Open(fname);
145    if ((!ef)||(!ef->IsOpen())) {
146       ::Error("AliTOFComparison.C","Can't open AliESDs.root !");
147       delete rl;
148       return 4;
149    }
150    TKey *key=0;
151    TIter next(ef->GetListOfKeys());
152
153
154
155    //******* Loop over events *********
156    Int_t e=0;
157    while ((key=(TKey*)next())!=0) {
158      cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
159
160      rl->GetEvent(e); ef->cd();
161
162      TTree *digTree=tofl->TreeD();
163      if (!digTree) {
164         cerr<<"Can't get the TOF cluster tree !\n";
165         return 3;
166      } 
167      TBranch *branch=digTree->GetBranch("TOF");
168      if (!branch) { 
169         cerr<<"Can't get the branch with the TOF digits !\n";
170         return 4;
171      }
172      TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;   
173      branch->SetAddress(&digits);
174
175      digTree->GetEvent(0);
176      Int_t nd=digits->GetEntriesFast();
177      cerr<<"Number of the TOF digits: "<<nd<<endl;
178
179
180
181      AliESD *event=(AliESD*)key->ReadObj();
182      Int_t ntrk=event->GetNumberOfTracks();
183      cerr<<"Number of ESD tracks : "<<ntrk<<endl; 
184
185
186      if (tofTree->GetEvent(e++)==0) {
187         cerr<<"No reconstructable tracks !\n";
188         continue;
189      }
190
191      Int_t ngood=refs->GetEntriesFast(); 
192
193      Int_t matched=0;
194      Int_t mismatched=0;
195      for (Int_t i=0; i<ngood; i++) {
196         AliTrackReference *ref=(AliTrackReference*)refs->UncheckedAt(i); 
197         Int_t lab=ref->Label();
198         Float_t ptg=TMath::Sqrt(ref->Px()*ref->Px() + ref->Py()*ref->Py());
199
200         Double_t tgl=ref->Pz()/ptg; //tan(lambda)
201
202          if (ptg>pmin) { hgood->Fill(ptg); hgoo->Fill(tgl); }
203
204          Int_t j;
205          AliESDtrack *t=0;
206          for (j=0; j<ntrk; j++) {
207              AliESDtrack *tt=event->GetTrack(j);
208              if (lab!=TMath::Abs(tt->GetLabel())) continue;
209              t=tt;
210              //if ((tt->GetStatus()&AliESDtrack::kTOFpid) == 0) continue;
211              if (tt->GetTOFsignal() < 0) continue;
212              UInt_t idx=tt->GetTOFcluster();
213              if ((Int_t)idx>=nd) {
214                cerr<<"Wrong digit index ! "<<idx<<endl;
215                return 5;
216              }
217              AliTOFdigit *dig=(AliTOFdigit*)digits->UncheckedAt(idx);
218              Int_t *label=dig->GetTracks();
219              if (label[0]!=lab)
220              if (label[1]!=lab)
221                if (label[2]!=lab) {
222                  mismatched++; 
223                  if (ptg>pmin) { hfake->Fill(ptg); hfak->Fill(tgl); } 
224                  break;
225                }
226              if (ptg>pmin) { hfound->Fill(ptg); hfoun->Fill(tgl); }
227              matched++;
228              break;
229          }
230          if (j==ntrk) {
231             cerr<<"Not matched: "<<lab<<"   ";
232             if (t) {
233                cerr<<(t->GetStatus()&AliESDtrack::kITSout)<<' '
234                    <<(t->GetStatus()&AliESDtrack::kTPCout)<<' '
235                    <<(t->GetStatus()&AliESDtrack::kTRDout)<<' '
236                    <<(t->GetStatus()&AliESDtrack::kTIME);
237             } else cerr<<"No ESD track !";
238             cerr<<endl;
239          }
240      }
241      cout<<"Number of good tracks: "<<ngood<<endl;
242      cout<<"Number of matched tracks: "<<matched<<endl;
243      cout<<"Number of mismatched tracks: "<<mismatched<<endl;
244
245      allgood+=ngood; allmatched+=matched; allmismatched+=mismatched;
246
247      refs->Clear();
248      delete event;
249    } //***** End of the loop over events
250
251    ef->Close();
252    
253    delete tofTree;
254    refFile->Close();
255
256    if (allgood!=0) cerr<<"\n\nEfficiency: "<<Float_t(allmatched)/allgood<<endl;
257    cout<<"Total number of good tracks: "<<allgood<<endl;
258    cout<<"Total number of matched tracks: "<<allmatched<<endl;
259    cout<<"Total number of mismatched tracks: "<<allmismatched<<endl;
260
261    TCanvas *c1=(TCanvas*)gROOT->FindObject("c1");
262    if (!c1) {
263       c1=new TCanvas("c1","",0,0,600,900);
264       c1->Divide(1,2);
265    }
266    hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2();
267    hgp->Divide(hfound,hgood,1,1.,"b");
268    hfp->Divide(hfake,hgood,1,1.,"b");
269    hgp->SetMaximum(1.4);
270    hgp->SetYTitle("Matching efficiency");
271    hgp->SetXTitle("Pt (GeV/c)");
272
273    hfoun->Sumw2(); hgoo->Sumw2(); hfak->Sumw2();
274    hgl->Divide(hfoun,hgoo,1,1.,"b");
275    hfl->Divide(hfak,hgoo,1,1.,"b");
276    hgl->SetMaximum(1.4);
277    hgl->SetYTitle("Matching efficiency");
278    hgl->SetXTitle("Tan(lambda)");
279
280    c1->cd(1);
281
282    hgp->Draw();
283    hfp->Draw("histsame");
284    TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4);
285    line1->Draw("same");
286    TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4);
287    line2->Draw("same");
288
289    c1->cd(2);
290
291    hgl->Draw();
292    hfl->Draw("histsame");
293    TLine *line3 = new TLine(-1,1.0,1,1.0); line3->SetLineStyle(4);
294    line3->Draw("same");
295    TLine *line4 = new TLine(-1,0.9,1,0.9); line4->SetLineStyle(4);
296    line4->Draw("same");
297
298    c1->Update();
299    
300    TFile fc("AliTOFComparison.root","RECREATE");
301    c1->Write();
302    fc.Close();
303
304    gBenchmark->Stop("AliTOFComparison");
305    gBenchmark->Show("AliTOFComparison");
306
307    delete rl;
308    return 0;
309 }
310
311 Int_t GoodTracksTOF(const Char_t *dir) {
312    if (gAlice) { 
313        delete gAlice->GetRunLoader();
314        delete gAlice;//if everything was OK here it is already NULL
315        gAlice = 0x0;
316    }
317
318    Char_t fname[100];
319    sprintf(fname,"%s/galice.root",dir);
320
321    AliRunLoader *rl = AliRunLoader::Open(fname,"COMPARISON");
322    if (!rl) {
323       ::Error("GoodTracksTOF","Can't start session !");
324       return 1;
325    }
326
327    rl->LoadgAlice();
328    rl->LoadHeader();
329    rl->LoadKinematics();
330
331
332    AliLoader* tofl = rl->GetLoader("TOFLoader");
333    if (tofl == 0x0) {
334       ::Error("GoodTracksTOF","Can not get the TOF loader !");
335       delete rl;
336       return 2;
337    }
338    tofl->LoadDigits("read");
339
340    Int_t nev=rl->GetNumberOfEvents();
341    ::Info("GoodTracksTOF","Number of events : %d\n",nev);  
342
343    sprintf(fname,"%s/GoodTracksITS.root",dir);
344    TFile *itsFile=TFile::Open(fname);
345    if ((!itsFile)||(!itsFile->IsOpen())) {
346        ::Error("GoodTracksTOF","Can't open the GoodTracksITS.root !");
347        delete rl;
348        return 5; 
349    }
350    TClonesArray dum("AliTrackReference",1000), *itsRefs=&dum;
351    TTree *itsTree=(TTree*)itsFile->Get("itsTree");
352    if (!itsTree) {
353        ::Error("GoodTracksTOF","Can't get the ITS reference tree !");
354        delete rl;
355        return 6;
356    }
357    TBranch *itsBranch=itsTree->GetBranch("ITS");
358    if (!itsBranch) {
359       ::Error("GoodTracksTOF","Can't get the ITS reference branch !");
360       delete rl;
361       return 7;
362    }
363    itsBranch->SetAddress(&itsRefs);
364
365
366    sprintf(fname,"%s/GoodTracksTOF.root",dir);
367    TFile *tofFile=TFile::Open(fname,"recreate");
368    TClonesArray dummy("AliTrackReference",1000), *tofRefs=&dummy;
369    TTree tofTree("tofTree","Tree with info about the reconstructable TOF tracks");
370    tofTree.Branch("TOF",&tofRefs);
371
372
373    //********  Loop over generated events 
374    for (Int_t e=0; e<nev; e++) {
375
376      rl->GetEvent(e); tofFile->cd();
377
378      Int_t np = rl->GetHeader()->GetNtrack();
379      cout<<"Event "<<e<<" Number of particles: "<<np<<endl;
380
381      //******** Fill the "good" masks
382      Int_t *good=new Int_t[np]; Int_t k; for (k=0; k<np; k++) good[k]=0;
383
384      TTree *dTree=tofl->TreeD();
385      if (!dTree) {
386         ::Error("GoodTracksTOF","Can't get the TOF cluster tree !");
387         delete rl;
388         return 8;
389      } 
390      TBranch *branch=dTree->GetBranch("TOF");
391      if (!branch) { 
392         ::Error("GoodTracksTOF","Can't get the branch with the TOF digits !");
393         return 9;
394      }
395      TClonesArray dummy("AliTOFdigit",10000), *digits=&dummy;
396      branch->SetAddress(&digits);
397    
398      dTree->GetEvent(0);
399      Int_t nd=digits->GetEntriesFast();
400
401      for (Int_t i=0; i<nd; i++) {
402        AliTOFdigit *d=(AliTOFdigit*)digits->UncheckedAt(i);
403        Int_t l0=d->GetTrack(0);
404           if (l0>=np) {cerr<<"Wrong label: "<<l0<<endl; continue;}
405        Int_t l1=d->GetTrack(1);
406           if (l1>=np) {cerr<<"Wrong label: "<<l1<<endl; continue;}
407        Int_t l2=d->GetTrack(2);
408           if (l2>=np) {cerr<<"Wrong label: "<<l2<<endl; continue;}
409        if (l0>=0) good[l0]++; 
410        if (l1>=0) good[l1]++; 
411        if (l2>=0) good[l2]++;
412      }
413      digits->Clear();
414
415
416      //****** select tracks which are "good" enough
417      AliStack* stack = rl->Stack();
418
419      itsTree->GetEvent(e);
420      Int_t nk=itsRefs->GetEntriesFast();
421
422      Int_t nt=0;
423      for (k=0; k<nk; k++) {
424         AliTrackReference *itsRef=(AliTrackReference *)itsRefs->UncheckedAt(k);
425         Int_t lab=itsRef->Label();
426         if (good[lab] == 0) continue;
427         TParticle *p = (TParticle*)stack->Particle(lab);
428         if (p == 0x0) {
429            cerr<<"Can not get particle "<<lab<<endl;
430            continue;
431         }
432
433         if (TMath::Abs(p->Vx())>0.1) continue;
434         if (TMath::Abs(p->Vy())>0.1) continue;
435         //if (TMath::Abs(p->Vz())>0.1) continue;
436
437         new((*tofRefs)[nt]) AliTrackReference(*itsRef);
438         nt++;
439      }
440      itsRefs->Clear();
441
442      tofTree.Fill();
443      tofRefs->Clear();
444
445      delete[] good;
446
447    } //*** end of the loop over generated events
448
449    tofTree.Write();
450    tofFile->Close();
451
452    delete itsTree;
453    itsFile->Close();
454
455    delete rl;
456    return 0;
457 }