Protection against zerro pointer
[u/mrichter/AliRoot.git] / PWG1 / TPC / AliPerformanceTPC.cxx
1 //------------------------------------------------------------------------------
2 // Implementation of AliPerformanceTPC class. It keeps information from 
3 // comparison of reconstructed and MC particle tracks. In addtion, 
4 // it keeps selection cuts used during comparison. The comparison 
5 // information is stored in the ROOT histograms. Analysis of these 
6 // histograms can be done by using Analyse() class function. The result of 
7 // the analysis (histograms/graphs) are stored in the folder which is
8 // a data member of AliPerformanceTPC.
9 //
10 // Author: J.Otwinowski 04/02/2008 
11 //------------------------------------------------------------------------------
12
13 /*
14  
15   // after running comparison task, read the file, and get component
16   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
17   LoadMyLibs();
18
19   TFile f("Output.root");
20   AliPerformanceTPC * compObj = (AliPerformanceTPC*)coutput->FindObject("AliPerformanceTPC");
21  
22   // analyse comparison data
23   compObj->Analyse();
24
25   // the output histograms/graphs will be stored in the folder "folderTPC" 
26   compObj->GetAnalysisFolder()->ls("*");
27
28   // user can save whole comparison object (or only folder with anlysed histograms) 
29   // in the seperate output file (e.g.)
30   TFile fout("Analysed_TPC.root","recreate");
31   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
32   fout.Close();
33
34 */
35
36 #include "TCanvas.h"
37 #include "TH1.h"
38 #include "TH2.h"
39 #include "TAxis.h"
40 #include "TPostScript.h"
41
42 #include "AliPerformanceTPC.h" 
43 #include "AliESDEvent.h" 
44 #include "AliESDVertex.h"
45 #include "AliESDtrack.h"
46 #include "AliLog.h" 
47 #include "AliMCEvent.h" 
48 #include "AliHeader.h" 
49 #include "AliGenEventHeader.h" 
50 #include "AliStack.h" 
51 #include "AliMCInfoCuts.h" 
52 #include "AliRecInfoCuts.h" 
53 #include "AliTracker.h" 
54 #include "AliTreeDraw.h" 
55 #include "AliTPCTransform.h" 
56 #include "AliTPCseed.h" 
57 #include "AliTPCcalibDB.h" 
58 #include "AliESDfriend.h" 
59 #include "AliESDfriendTrack.h" 
60 #include "AliTPCclusterMI.h" 
61
62 using namespace std;
63
64 ClassImp(AliPerformanceTPC)
65
66 //_____________________________________________________________________________
67 AliPerformanceTPC::AliPerformanceTPC():
68   AliPerformanceObject("AliPerformanceTPC"),
69   fTPCClustHisto(0),
70   fTPCEventHisto(0),
71   fTPCTrackHisto(0),
72
73   // Cuts 
74   fCutsRC(0),  
75   fCutsMC(0),  
76
77   // histogram folder 
78   fAnalysisFolder(0)
79 {
80   Init();
81 }
82
83 //_____________________________________________________________________________
84 AliPerformanceTPC::AliPerformanceTPC(Char_t* name="AliPerformanceTPC", Char_t* title="AliPerformanceTPC",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
85   AliPerformanceObject(name,title),
86   fTPCClustHisto(0),
87   fTPCEventHisto(0),
88   fTPCTrackHisto(0),
89
90   // Cuts 
91   fCutsRC(0),  
92   fCutsMC(0),  
93
94   // histogram folder 
95   fAnalysisFolder(0)
96 {
97   // named constructor  
98   // 
99   SetAnalysisMode(analysisMode);
100   SetHptGenerator(hptGenerator);
101
102   Init();
103 }
104
105 //_____________________________________________________________________________
106 AliPerformanceTPC::~AliPerformanceTPC()
107 {
108   // destructor
109    
110   if(fTPCClustHisto) delete fTPCClustHisto; fTPCClustHisto=0;     
111   if(fTPCEventHisto) delete fTPCEventHisto; fTPCEventHisto=0;     
112   if(fTPCTrackHisto) delete fTPCTrackHisto; fTPCTrackHisto=0;     
113   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
114 }
115
116 //_____________________________________________________________________________
117 void AliPerformanceTPC::Init(){
118   //
119   // histogram bining
120   //
121
122   // set pt bins
123   Int_t nPtBins = 50;
124   Double_t ptMin = 1.e-2, ptMax = 10.;
125
126   Double_t *binsPt = 0;
127   if (IsHptGenerator())  { 
128     nPtBins = 100; ptMax = 100.;
129     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
130   } else {
131     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
132   }
133
134   /*
135   Int_t nPtBins = 31;
136   Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
137   Double_t ptMin = 0., ptMax = 10.; 
138
139   if(IsHptGenerator() == kTRUE) {
140     nPtBins = 100;
141     ptMin = 0.; ptMax = 100.; 
142   }
143   */
144   // 
145   //gclX:gclY:TPCSide
146   Int_t binsTPCClustHisto[3]=  { 500,   500,  2 };
147   Double_t minTPCClustHisto[3]={-250., -250., 0.};
148   Double_t maxTPCClustHisto[3]={ 250.,  250., 2.};
149
150   fTPCClustHisto = new THnSparseF("fTPCClustHisto","gclX:gclY:TPCSide",3,binsTPCClustHisto,minTPCClustHisto,maxTPCClustHisto);
151   fTPCClustHisto->GetAxis(0)->SetTitle("gclX (cm)");
152   fTPCClustHisto->GetAxis(1)->SetTitle("gclY (cm)");
153   fTPCClustHisto->GetAxis(2)->SetTitle("TPCSide");
154   fTPCClustHisto->Sumw2();
155  
156
157   // Xv:Yv:Zv:mult:multP:multN:vertStatus
158   Int_t binsTPCEventHisto[7]=  {100,  100,   100,  151,   151,   151, 2   };
159   Double_t minTPCEventHisto[7]={-10., -10., -30.,  -0.5,  -0.5,  -0.5, 0.  };
160   Double_t maxTPCEventHisto[7]={ 10.,  10.,  30.,  150.5, 150.5, 150.5, 2. };
161
162   fTPCEventHisto = new THnSparseF("fTPCEventHisto","Xv:Yv:Zv:mult:multP:multN:vertStatus",7,binsTPCEventHisto,minTPCEventHisto,maxTPCEventHisto);
163   fTPCEventHisto->GetAxis(0)->SetTitle("Xv (cm)");
164   fTPCEventHisto->GetAxis(1)->SetTitle("Yv (cm)");
165   fTPCEventHisto->GetAxis(2)->SetTitle("Zv (cm)");
166   fTPCEventHisto->GetAxis(3)->SetTitle("mult");
167   fTPCEventHisto->GetAxis(4)->SetTitle("multP");
168   fTPCEventHisto->GetAxis(5)->SetTitle("multN");
169   fTPCEventHisto->GetAxis(6)->SetTitle("vertStatus");
170   fTPCEventHisto->Sumw2();
171
172
173   // nTPCClust:chi2PerTPCClust:nTPCClustFindRatio:DCAr:DCAz:eta:phi:pt:charge
174   Int_t binsTPCTrackHisto[9]=  { 160,  50,  60,  100, 100,  30,   144,             nPtBins,  3 };
175   Double_t minTPCTrackHisto[9]={ 0.,   0.,  0., -10,  -10., -1.5, 0.,             ptMin,   -1.5};
176   Double_t maxTPCTrackHisto[9]={ 160., 10., 1.2, 10,   10.,  1.5, 2.*TMath::Pi(), ptMax,    1.5};
177
178   fTPCTrackHisto = new THnSparseF("fTPCTrackHisto","nClust:chi2PerClust:nClust/nFindableClust:DCAr:DCAz:eta:phi:pt:charge",9,binsTPCTrackHisto,minTPCTrackHisto,maxTPCTrackHisto);
179   fTPCTrackHisto->SetBinEdges(7,binsPt);
180
181   fTPCTrackHisto->GetAxis(0)->SetTitle("nClust");
182   fTPCTrackHisto->GetAxis(1)->SetTitle("chi2PerClust");
183   fTPCTrackHisto->GetAxis(2)->SetTitle("nClust/nFindableClust");
184   fTPCTrackHisto->GetAxis(3)->SetTitle("DCAr (cm)");
185   fTPCTrackHisto->GetAxis(4)->SetTitle("DCAz (cm)");
186   fTPCTrackHisto->GetAxis(5)->SetTitle("#eta");
187   fTPCTrackHisto->GetAxis(6)->SetTitle("#phi (rad)");
188   fTPCTrackHisto->GetAxis(7)->SetTitle("p_{T} (GeV/c)");
189   fTPCTrackHisto->GetAxis(8)->SetTitle("charge");
190   fTPCTrackHisto->Sumw2();
191
192   // Init cuts 
193   if(!fCutsMC) 
194     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
195   if(!fCutsRC) 
196     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
197
198   // init folder
199   fAnalysisFolder = CreateFolder("folderTPC","Analysis Resolution Folder");
200 }
201
202 //_____________________________________________________________________________
203 void AliPerformanceTPC::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
204 {
205   if(!esdTrack) return;
206
207   // Fill TPC only resolution comparison information 
208   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
209   if(!track) return;
210
211   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
212   esdTrack->GetImpactParametersTPC(dca,cov);
213
214   Float_t q = esdTrack->Charge();
215   Float_t pt = track->Pt();
216   Float_t eta = track->Eta();
217   Float_t phi = track->Phi();
218   Int_t nClust = esdTrack->GetTPCclusters(0);
219   Int_t nFindableClust = esdTrack->GetTPCNclsF();
220
221   Float_t chi2PerCluster = 0.;
222   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
223
224   Float_t clustPerFindClust = 0.;
225   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
226   
227   //
228   // select primaries
229   //
230   Double_t dcaToVertex = -1;
231   if( fCutsRC->GetDCAToVertex2D() ) 
232   {
233       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()                    + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ()); 
234   }
235   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
236   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
237   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
238
239   Double_t vTPCTrackHisto[9] = {nClust,chi2PerCluster,clustPerFindClust,dca[0],dca[1],eta,phi,pt,q};
240   fTPCTrackHisto->Fill(vTPCTrackHisto); 
241  
242   //
243   // Fill rec vs MC information
244   //
245   if(!stack) return;
246
247 }
248
249 //_____________________________________________________________________________
250 void AliPerformanceTPC::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
251 {
252   // Fill comparison information (TPC+ITS) 
253   AliDebug(AliLog::kWarning, "Warning: Not implemented");
254 }
255  
256 //_____________________________________________________________________________
257 void AliPerformanceTPC::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
258 {
259   // Fill comparison information (constarained parameters) 
260   AliDebug(AliLog::kWarning, "Warning: Not implemented");
261 }
262  
263 //_____________________________________________________________________________
264 void AliPerformanceTPC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
265 {
266   // Process comparison information 
267   //
268   if(!esdEvent) 
269   {
270     Error("Exec","esdEvent not available");
271     return;
272   }
273   AliHeader* header = 0;
274   AliGenEventHeader* genHeader = 0;
275   AliStack* stack = 0;
276   TArrayF vtxMC(3);
277   
278   if(bUseMC)
279   {
280     if(!mcEvent) {
281       Error("Exec","mcEvent not available");
282       return;
283     }
284     // get MC event header
285     header = mcEvent->Header();
286     if (!header) {
287       Error("Exec","Header not available");
288       return;
289     }
290     // MC particle stack
291     stack = mcEvent->Stack();
292     if (!stack) {
293       Error("Exec","Stack not available");
294       return;
295     }
296     // get MC vertex
297     genHeader = header->GenEventHeader();
298     if (!genHeader) {
299       Error("Exec","Could not retrieve genHeader from Header");
300       return;
301     }
302     genHeader->PrimaryVertex(vtxMC);
303   } 
304   
305   // use ESD friends
306   if(bUseESDfriend) {
307     if(!esdFriend) {
308       Error("Exec","esdFriend not available");
309       return;
310     }
311   }
312
313   // trigger
314   if(!bUseMC &&GetTriggerClass()) {
315     Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
316     if(!isEventTriggered) return; 
317   }
318
319   // get TPC event vertex
320   const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
321
322   //  events with rec. vertex
323   Int_t mult=0; Int_t multP=0; Int_t multN=0;
324   if(vtxESD->GetStatus() >0)
325   {
326   //  Process ESD events
327   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
328   { 
329     AliESDtrack *track = esdEvent->GetTrack(iTrack);
330     if(!track) continue;
331
332     if(bUseESDfriend) {
333     AliESDfriendTrack *friendTrack=esdFriend->GetTrack(iTrack);
334     if(!friendTrack) continue;
335
336       TObject *calibObject=0;
337       AliTPCseed *seed=0;
338       if (!friendTrack) continue;
339         for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j) {
340           if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
341           break;
342         }
343
344         //AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
345         for (Int_t irow=0;irow<159;irow++) {
346         if(!seed) continue;
347           
348           AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
349           if (!cluster) continue;
350
351               Float_t gclf[3];
352               cluster->GetGlobalXYZ(gclf);
353
354               //Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
355               //Int_t i[1]={cluster->GetDetector()};
356               //transform->Transform(x,i,0,1);
357               //printf("gx %f gy  %f  gz %f \n", cluster->GetX(), cluster->GetY(),cluster->GetZ());
358               //printf("gclf[0] %f gclf[1]  %f  gclf[2] %f \n", gclf[0], gclf[1],  gclf[2]);
359      
360              Int_t TPCside; 
361              if(gclf[2]>0.) TPCside=0; // A side 
362              else TPCside=1;
363
364              Double_t vTPCClust[3] = { gclf[0], gclf[1],  TPCside };
365              fTPCClustHisto->Fill(vTPCClust);
366          }
367     }
368     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
369     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
370     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
371     else {
372       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
373       return;
374     }
375
376    // TPC only
377    AliESDtrack *tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
378    if(!tpcTrack) continue;
379
380    // track selection
381    if( fCutsRC->AcceptTrack(tpcTrack) ) { 
382      mult++;
383      if(tpcTrack->Charge()>0.) multP++;
384      if(tpcTrack->Charge()<0.) multN++;
385    }
386
387    if(tpcTrack) delete tpcTrack;
388   }
389   }
390   //
391   
392   Double_t vTPCEvent[7] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),mult,multP,multN,vtxESD->GetStatus()};
393   fTPCEventHisto->Fill(vTPCEvent);
394 }
395
396 //_____________________________________________________________________________
397 void AliPerformanceTPC::Analyse() {
398   //
399   // Analyse comparison information and store output histograms
400   // in the folder "folderTPC"
401   //
402   TH1::AddDirectory(kFALSE);
403   TH1F *h=0;
404   TH2D *h2D=0;
405   TObjArray *aFolderObj = new TObjArray;
406   char name[256];
407   char title[256];
408
409   //
410   // Cluster histograms
411   //
412   fTPCClustHisto->GetAxis(2)->SetRange(1,1); // A-side
413   h2D = fTPCClustHisto->Projection(1,0);
414   h2D->SetName("h_clust_A_side");
415   h2D->SetTitle("gclX:gclY - A_side");
416   aFolderObj->Add(h2D);
417
418   fTPCClustHisto->GetAxis(2)->SetRange(2,2); // C-side
419   h2D = fTPCClustHisto->Projection(1,0);
420   h2D->SetName("h_clust_C_side");
421   h2D->SetTitle("gclX:gclY - C_side");
422   aFolderObj->Add(h2D);
423
424   //
425   // event histograms
426   //
427   for(Int_t i=0; i<6; i++) 
428   {
429       h = (TH1F*)fTPCEventHisto->Projection(i);
430       sprintf(name,"h_tpc_event_%d",i);
431       h->SetName(name);
432       h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
433       h->GetYaxis()->SetTitle("events");
434       sprintf(title,"%s",fTPCEventHisto->GetAxis(i)->GetTitle());
435       h->SetTitle(title);
436
437       aFolderObj->Add(h);
438   }
439
440   // reconstructed vertex status > 0
441   fTPCEventHisto->GetAxis(6)->SetRange(2,2);
442   for(Int_t i=0; i<6; i++) 
443   {
444       h = (TH1F*)fTPCEventHisto->Projection(i);
445       sprintf(name,"h_tpc_event_recVertex%d",i);
446       h->SetName(name);
447       h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
448       h->GetYaxis()->SetTitle("events");
449       sprintf(title,"%s rec. vertex",fTPCEventHisto->GetAxis(i)->GetTitle());
450       h->SetTitle(title);
451
452       aFolderObj->Add(h);
453   }
454
455   //
456   // Track histograms
457   //
458   for(Int_t i=0; i<9; i++) 
459   {
460       h = (TH1F*)fTPCTrackHisto->Projection(i);
461       sprintf(name,"h_tpc_track_%d",i);
462       h->SetName(name);
463       h->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
464       h->GetYaxis()->SetTitle("tracks");
465       sprintf(title,"%s",fTPCTrackHisto->GetAxis(i)->GetTitle());
466       h->SetTitle(title);
467
468       if(i==7) h->Scale(1,"width");
469       aFolderObj->Add(h);
470   }
471
472   //
473   for(Int_t i=0; i<8; i++) 
474   {
475     for(Int_t j=i+1; j<9; j++) 
476     {
477       h2D = fTPCTrackHisto->Projection(i,j);
478       sprintf(name,"h_tpc_track_%d_vs_%d",i,j);
479       h2D->SetName(name);
480       h2D->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(j)->GetTitle());
481       h2D->GetYaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
482       sprintf(title,"%s vs %s",fTPCTrackHisto->GetAxis(j)->GetTitle(),fTPCTrackHisto->GetAxis(i)->GetTitle());
483       h2D->SetTitle(title);
484
485       if(j==7) h2D->SetBit(TH1::kLogX);
486       aFolderObj->Add(h2D);
487     }  
488   }
489
490   // export objects to analysis folder
491   fAnalysisFolder = ExportToFolder(aFolderObj);
492
493   // delete only TObjArray
494   if(aFolderObj) delete aFolderObj;
495 }
496
497 //_____________________________________________________________________________
498 TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array) 
499 {
500   // recreate folder avery time and export objects to new one
501   //
502   AliPerformanceTPC * comp=this;
503   TFolder *folder = comp->GetAnalysisFolder();
504
505   TString name, title;
506   TFolder *newFolder = 0;
507   Int_t i = 0;
508   Int_t size = array->GetSize();
509
510   if(folder) { 
511      // get name and title from old folder
512      name = folder->GetName();  
513      title = folder->GetTitle();  
514
515          // delete old one
516      delete folder;
517
518          // create new one
519      newFolder = CreateFolder(name.Data(),title.Data());
520      newFolder->SetOwner();
521
522          // add objects to folder
523      while(i < size) {
524            newFolder->Add(array->At(i));
525            i++;
526          }
527   }
528
529 return newFolder;
530 }
531
532 //_____________________________________________________________________________
533 Long64_t AliPerformanceTPC::Merge(TCollection* const list) 
534 {
535   // Merge list of objects (needed by PROOF)
536
537   if (!list)
538   return 0;
539
540   if (list->IsEmpty())
541   return 1;
542
543   TIterator* iter = list->MakeIterator();
544   TObject* obj = 0;
545
546   // collection of generated histograms
547   Int_t count=0;
548   while((obj = iter->Next()) != 0) 
549   {
550     AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
551     if (entry == 0) continue; 
552
553     fTPCClustHisto->Add(entry->fTPCClustHisto);
554     fTPCEventHisto->Add(entry->fTPCEventHisto);
555     fTPCTrackHisto->Add(entry->fTPCTrackHisto);
556
557     count++;
558   }
559
560 return count;
561 }
562
563 //_____________________________________________________________________________
564 TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) { 
565 // create folder for analysed histograms
566 //
567 TFolder *folder = 0;
568   folder = new TFolder(name.Data(),title.Data());
569
570   return folder;
571 }