changes in Analyse function
[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   Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
315   if(!isEventTriggered) return; 
316
317   // get TPC event vertex
318   const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
319
320   //  events with rec. vertex
321   Int_t mult=0; Int_t multP=0; Int_t multN=0;
322   if(vtxESD->GetStatus() >0)
323   {
324   //  Process ESD events
325   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
326   { 
327     AliESDtrack *track = esdEvent->GetTrack(iTrack);
328     if(!track) continue;
329
330     if(bUseESDfriend) {
331     AliESDfriendTrack *friendTrack=esdFriend->GetTrack(iTrack);
332     if(!friendTrack) continue;
333
334       TObject *calibObject=0;
335       AliTPCseed *seed=0;
336       if (!friendTrack) continue;
337         for (Int_t j=0;(calibObject=friendTrack->GetCalibObject(j));++j) {
338           if ((seed=dynamic_cast<AliTPCseed*>(calibObject)))
339           break;
340         }
341
342         //AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
343         for (Int_t irow=0;irow<159;irow++) {
344         if(!seed) continue;
345           
346           AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
347           if (!cluster) continue;
348
349               Float_t gclf[3];
350               cluster->GetGlobalXYZ(gclf);
351
352               //Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
353               //Int_t i[1]={cluster->GetDetector()};
354               //transform->Transform(x,i,0,1);
355               //printf("gx %f gy  %f  gz %f \n", cluster->GetX(), cluster->GetY(),cluster->GetZ());
356               //printf("gclf[0] %f gclf[1]  %f  gclf[2] %f \n", gclf[0], gclf[1],  gclf[2]);
357      
358              Int_t TPCside; 
359              if(gclf[2]>0.) TPCside=0; // A side 
360              else TPCside=1;
361
362              Double_t vTPCClust[3] = { gclf[0], gclf[1],  TPCside };
363              fTPCClustHisto->Fill(vTPCClust);
364          }
365     }
366     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
367     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
368     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
369     else {
370       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
371       return;
372     }
373
374    // TPC only
375    AliESDtrack *tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
376    if(!tpcTrack) continue;
377
378    // track selection
379    if( fCutsRC->AcceptTrack(tpcTrack) ) { 
380      mult++;
381      if(tpcTrack->Charge()>0.) multP++;
382      if(tpcTrack->Charge()<0.) multN++;
383    }
384
385    if(tpcTrack) delete tpcTrack;
386   }
387   }
388   //
389   
390   Double_t vTPCEvent[7] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),mult,multP,multN,vtxESD->GetStatus()};
391   fTPCEventHisto->Fill(vTPCEvent);
392 }
393
394 //_____________________________________________________________________________
395 void AliPerformanceTPC::Analyse() {
396   //
397   // Analyse comparison information and store output histograms
398   // in the folder "folderTPC"
399   //
400   TH1::AddDirectory(kFALSE);
401   TH1F *h=0;
402   TH2D *h2D=0;
403   TObjArray *aFolderObj = new TObjArray;
404   char name[256];
405   char title[256];
406
407   //
408   // Cluster histograms
409   //
410   fTPCClustHisto->GetAxis(2)->SetRange(1,1); // A-side
411   h2D = fTPCClustHisto->Projection(1,0);
412   h2D->SetName("h_clust_A_side");
413   h2D->SetTitle("gclX:gclY - A_side");
414   aFolderObj->Add(h2D);
415
416   fTPCClustHisto->GetAxis(2)->SetRange(2,2); // C-side
417   h2D = fTPCClustHisto->Projection(1,0);
418   h2D->SetName("h_clust_C_side");
419   h2D->SetTitle("gclX:gclY - C_side");
420   aFolderObj->Add(h2D);
421
422   //
423   // event histograms
424   //
425   for(Int_t i=0; i<6; i++) 
426   {
427       h = (TH1F*)fTPCEventHisto->Projection(i);
428       sprintf(name,"h_tpc_event_%d",i);
429       h->SetName(name);
430       h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
431       h->GetYaxis()->SetTitle("events");
432       sprintf(title,"%s",fTPCEventHisto->GetAxis(i)->GetTitle());
433       h->SetTitle(title);
434
435       aFolderObj->Add(h);
436   }
437
438   // reconstructed vertex status > 0
439   fTPCEventHisto->GetAxis(6)->SetRange(2,2);
440   for(Int_t i=0; i<6; i++) 
441   {
442       h = (TH1F*)fTPCEventHisto->Projection(i);
443       sprintf(name,"h_tpc_event_recVertex%d",i);
444       h->SetName(name);
445       h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
446       h->GetYaxis()->SetTitle("events");
447       sprintf(title,"%s rec. vertex",fTPCEventHisto->GetAxis(i)->GetTitle());
448       h->SetTitle(title);
449
450       aFolderObj->Add(h);
451   }
452
453   //
454   // Track histograms
455   //
456   for(Int_t i=0; i<9; i++) 
457   {
458       h = (TH1F*)fTPCTrackHisto->Projection(i);
459       sprintf(name,"h_tpc_track_%d",i);
460       h->SetName(name);
461       h->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
462       h->GetYaxis()->SetTitle("tracks");
463       sprintf(title,"%s",fTPCTrackHisto->GetAxis(i)->GetTitle());
464       h->SetTitle(title);
465
466       if(i==7) h->Scale(1,"width");
467       aFolderObj->Add(h);
468   }
469
470   //
471   for(Int_t i=0; i<8; i++) 
472   {
473     for(Int_t j=i+1; j<9; j++) 
474     {
475       h2D = fTPCTrackHisto->Projection(i,j);
476       sprintf(name,"h_tpc_track_%d_vs_%d",i,j);
477       h2D->SetName(name);
478       h2D->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(j)->GetTitle());
479       h2D->GetYaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
480       sprintf(title,"%s vs %s",fTPCTrackHisto->GetAxis(j)->GetTitle(),fTPCTrackHisto->GetAxis(i)->GetTitle());
481       h2D->SetTitle(title);
482
483       if(j==7) h2D->SetBit(TH1::kLogX);
484       aFolderObj->Add(h2D);
485     }  
486   }
487
488   // export objects to analysis folder
489   fAnalysisFolder = ExportToFolder(aFolderObj);
490
491   // delete only TObjArray
492   if(aFolderObj) delete aFolderObj;
493 }
494
495 //_____________________________________________________________________________
496 TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array) 
497 {
498   // recreate folder avery time and export objects to new one
499   //
500   AliPerformanceTPC * comp=this;
501   TFolder *folder = comp->GetAnalysisFolder();
502
503   TString name, title;
504   TFolder *newFolder = 0;
505   Int_t i = 0;
506   Int_t size = array->GetSize();
507
508   if(folder) { 
509      // get name and title from old folder
510      name = folder->GetName();  
511      title = folder->GetTitle();  
512
513          // delete old one
514      delete folder;
515
516          // create new one
517      newFolder = CreateFolder(name.Data(),title.Data());
518      newFolder->SetOwner();
519
520          // add objects to folder
521      while(i < size) {
522            newFolder->Add(array->At(i));
523            i++;
524          }
525   }
526
527 return newFolder;
528 }
529
530 //_____________________________________________________________________________
531 Long64_t AliPerformanceTPC::Merge(TCollection* const list) 
532 {
533   // Merge list of objects (needed by PROOF)
534
535   if (!list)
536   return 0;
537
538   if (list->IsEmpty())
539   return 1;
540
541   TIterator* iter = list->MakeIterator();
542   TObject* obj = 0;
543
544   // collection of generated histograms
545   Int_t count=0;
546   while((obj = iter->Next()) != 0) 
547   {
548     AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
549     if (entry == 0) continue; 
550
551     fTPCClustHisto->Add(entry->fTPCClustHisto);
552     fTPCEventHisto->Add(entry->fTPCEventHisto);
553     fTPCTrackHisto->Add(entry->fTPCTrackHisto);
554
555     count++;
556   }
557
558 return count;
559 }
560
561 //_____________________________________________________________________________
562 TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) { 
563 // create folder for analysed histograms
564 //
565 TFolder *folder = 0;
566   folder = new TFolder(name.Data(),title.Data());
567
568   return folder;
569 }