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