new histograms added
[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
56 using namespace std;
57
58 ClassImp(AliPerformanceTPC)
59
60 //_____________________________________________________________________________
61 AliPerformanceTPC::AliPerformanceTPC():
62   AliPerformanceObject("AliPerformanceTPC"),
63   fTPCEventHisto(0),
64   fTPCTrackHisto(0),
65
66   // Cuts 
67   fCutsRC(0),  
68   fCutsMC(0),  
69
70   // histogram folder 
71   fAnalysisFolder(0)
72 {
73   Init();
74 }
75
76 //_____________________________________________________________________________
77 AliPerformanceTPC::AliPerformanceTPC(Char_t* name="AliPerformanceTPC", Char_t* title="AliPerformanceTPC",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
78   AliPerformanceObject(name,title),
79   fTPCEventHisto(0),
80   fTPCTrackHisto(0),
81
82   // Cuts 
83   fCutsRC(0),  
84   fCutsMC(0),  
85
86   // histogram folder 
87   fAnalysisFolder(0)
88 {
89   // named constructor  
90   // 
91   SetAnalysisMode(analysisMode);
92   SetHptGenerator(hptGenerator);
93
94   Init();
95 }
96
97 //_____________________________________________________________________________
98 AliPerformanceTPC::~AliPerformanceTPC()
99 {
100   // destructor
101    
102   if(fTPCEventHisto) delete fTPCEventHisto; fTPCEventHisto=0;     
103   if(fTPCTrackHisto) delete fTPCTrackHisto; fTPCTrackHisto=0;     
104   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
105 }
106
107 //_____________________________________________________________________________
108 void AliPerformanceTPC::Init(){
109   //
110   // histogram bining
111   //
112
113   // set pt bins
114   Int_t nPtBins = 50;
115   Double_t ptMin = 1.e-2, ptMax = 10.;
116
117   Double_t *binsPt = 0;
118   if (IsHptGenerator())  { 
119     nPtBins = 100; ptMax = 100.;
120     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
121   } else {
122     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
123   }
124
125   /*
126   Int_t nPtBins = 31;
127   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.};
128   Double_t ptMin = 0., ptMax = 10.; 
129
130   if(IsHptGenerator() == kTRUE) {
131     nPtBins = 100;
132     ptMin = 0.; ptMax = 100.; 
133   }
134   */
135
136   // Xv:Yv:Zv:mult:multP:multN:vertStatus
137   Int_t binsTPCEventHisto[7]=  {100,  100,   100,  151,   151,   151, 2   };
138   Double_t minTPCEventHisto[7]={-10., -10., -30.,  -0.5,  -0.5,  -0.5, 0.  };
139   Double_t maxTPCEventHisto[7]={ 10.,  10.,  30.,  150.5, 150.5, 150.5, 2. };
140
141   fTPCEventHisto = new THnSparseF("fTPCEventHisto","Xv:Yv:Zv:mult:multP:multN:vertStatus",7,binsTPCEventHisto,minTPCEventHisto,maxTPCEventHisto);
142   fTPCEventHisto->GetAxis(0)->SetTitle("Xv (cm)");
143   fTPCEventHisto->GetAxis(1)->SetTitle("Yv (cm)");
144   fTPCEventHisto->GetAxis(2)->SetTitle("Zv (cm)");
145   fTPCEventHisto->GetAxis(3)->SetTitle("mult");
146   fTPCEventHisto->GetAxis(4)->SetTitle("multP");
147   fTPCEventHisto->GetAxis(5)->SetTitle("multN");
148   fTPCEventHisto->GetAxis(6)->SetTitle("vertStatus");
149   fTPCEventHisto->Sumw2();
150
151
152   // nTPCClust:chi2PerTPCClust:nTPCClustFindRatio:DCAr:DCAz:eta:phi:pt:charge
153   Int_t binsTPCTrackHisto[9]=  { 160,  50,  60,  100, 100,  30,   90,             nPtBins,  3 };
154   Double_t minTPCTrackHisto[9]={ 0.,   0.,  0., -10,  -10., -1.5, 0.,             ptMin,   -1.5};
155   Double_t maxTPCTrackHisto[9]={ 160., 10., 1.2, 10,   10.,  1.5, 2.*TMath::Pi(), ptMax,    1.5};
156
157   fTPCTrackHisto = new THnSparseF("fTPCTrackHisto","nClust:chi2PerClust:nClust/nFindableClust:DCAr:DCAz:eta:phi:pt:charge",9,binsTPCTrackHisto,minTPCTrackHisto,maxTPCTrackHisto);
158   fTPCTrackHisto->SetBinEdges(7,binsPt);
159
160   fTPCTrackHisto->GetAxis(0)->SetTitle("nClust");
161   fTPCTrackHisto->GetAxis(1)->SetTitle("chi2PerClust");
162   fTPCTrackHisto->GetAxis(2)->SetTitle("nClust/nFindableClust");
163   fTPCTrackHisto->GetAxis(3)->SetTitle("DCAr (cm)");
164   fTPCTrackHisto->GetAxis(4)->SetTitle("DCAz (cm)");
165   fTPCTrackHisto->GetAxis(5)->SetTitle("#eta");
166   fTPCTrackHisto->GetAxis(6)->SetTitle("#phi (rad)");
167   fTPCTrackHisto->GetAxis(7)->SetTitle("p_{T} (GeV/c)");
168   fTPCTrackHisto->GetAxis(8)->SetTitle("charge");
169   fTPCTrackHisto->Sumw2();
170
171   // Init cuts 
172   if(!fCutsMC) 
173     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
174   if(!fCutsRC) 
175     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
176
177   // init folder
178   fAnalysisFolder = CreateFolder("folderTPC","Analysis Resolution Folder");
179 }
180
181 //_____________________________________________________________________________
182 void AliPerformanceTPC::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
183 {
184   if(!esdTrack) return;
185
186   // Fill TPC only resolution comparison information 
187   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
188   if(!track) return;
189
190   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
191   esdTrack->GetImpactParametersTPC(dca,cov);
192
193   Float_t q = esdTrack->Charge();
194   Float_t pt = track->Pt();
195   Float_t eta = track->Eta();
196   Float_t phi = track->Phi();
197   Int_t nClust = esdTrack->GetTPCclusters(0);
198   Int_t nFindableClust = esdTrack->GetTPCNclsF();
199
200   Float_t chi2PerCluster = 0.;
201   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
202
203   Float_t clustPerFindClust = 0.;
204   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
205   
206   //
207   // select primaries
208   //
209   Double_t dcaToVertex = -1;
210   if( fCutsRC->GetDCAToVertex2D() ) 
211   {
212       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()                    + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ()); 
213   }
214   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
215   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
216   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
217
218   Double_t vTPCTrackHisto[9] = {nClust,chi2PerCluster,clustPerFindClust,dca[0],dca[1],eta,phi,pt,q};
219   fTPCTrackHisto->Fill(vTPCTrackHisto); 
220  
221   //
222   // Fill rec vs MC information
223   //
224   if(!stack) return;
225
226 }
227
228 //_____________________________________________________________________________
229 void AliPerformanceTPC::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
230 {
231   // Fill comparison information (TPC+ITS) 
232   AliDebug(AliLog::kWarning, "Warning: Not implemented");
233 }
234  
235 //_____________________________________________________________________________
236 void AliPerformanceTPC::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
237 {
238   // Fill comparison information (constarained parameters) 
239   AliDebug(AliLog::kWarning, "Warning: Not implemented");
240 }
241  
242 //_____________________________________________________________________________
243 void AliPerformanceTPC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
244 {
245   // Process comparison information 
246   //
247   if(!esdEvent) 
248   {
249     Error("Exec","esdEvent not available");
250     return;
251   }
252   AliHeader* header = 0;
253   AliGenEventHeader* genHeader = 0;
254   AliStack* stack = 0;
255   TArrayF vtxMC(3);
256   
257   if(bUseMC)
258   {
259     if(!mcEvent) {
260       Error("Exec","mcEvent not available");
261       return;
262     }
263     // get MC event header
264     header = mcEvent->Header();
265     if (!header) {
266       Error("Exec","Header not available");
267       return;
268     }
269     // MC particle stack
270     stack = mcEvent->Stack();
271     if (!stack) {
272       Error("Exec","Stack not available");
273       return;
274     }
275     // get MC vertex
276     genHeader = header->GenEventHeader();
277     if (!genHeader) {
278       Error("Exec","Could not retrieve genHeader from Header");
279       return;
280     }
281     genHeader->PrimaryVertex(vtxMC);
282   } 
283   
284   // use ESD friends
285   if(bUseESDfriend) {
286     if(!esdFriend) {
287       Error("Exec","esdFriend not available");
288       return;
289     }
290   }
291
292   //  Process events
293   Int_t mult=0; Int_t multP=0; Int_t multN=0;
294   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
295   { 
296     AliESDtrack *track = esdEvent->GetTrack(iTrack);
297     if(!track) continue;
298
299     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
300     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
301     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
302     else {
303       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
304       return;
305     }
306
307    // TPC only
308    AliESDtrack *tpcTrack = AliESDtrackCuts::GetTPCOnlyTrack(esdEvent,iTrack);
309    if(!tpcTrack) continue;
310
311    // track selection
312    if( fCutsRC->AcceptTrack(tpcTrack) ) { 
313      mult++;
314      if(tpcTrack->Charge()>0.) multP++;
315      if(tpcTrack->Charge()<0.) multN++;
316    }
317
318    if(tpcTrack) delete tpcTrack;
319   }
320   //
321   const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
322   Double_t vTPCEvent[7] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),mult,multP,multN,vtxESD->GetStatus()};
323   fTPCEventHisto->Fill(vTPCEvent);
324 }
325
326 //_____________________________________________________________________________
327 void AliPerformanceTPC::Analyse() {
328   //
329   // Analyse comparison information and store output histograms
330   // in the folder "folderTPC"
331   //
332   TH1::AddDirectory(kFALSE);
333   TH1F *h=0;
334   TH2F *h2D=0;
335   TObjArray *aFolderObj = new TObjArray;
336   char name[256];
337   char title[256];
338
339   //
340   // event histograms
341   //
342   for(Int_t i=0; i<6; i++) 
343   {
344       h = (TH1F*)fTPCEventHisto->Projection(i);
345       sprintf(name,"h_tpc_event_%d",i);
346       h->SetName(name);
347       h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
348       h->GetYaxis()->SetTitle("events");
349       sprintf(title,"%s",fTPCEventHisto->GetAxis(i)->GetTitle());
350       h->SetTitle(title);
351
352       aFolderObj->Add(h);
353   }
354
355   // reconstructed vertex status > 0
356   fTPCEventHisto->GetAxis(6)->SetRange(2,2);
357   for(Int_t i=0; i<6; i++) 
358   {
359       h = (TH1F*)fTPCEventHisto->Projection(i);
360       sprintf(name,"h_tpc_event_recVertex%d",i);
361       h->SetName(name);
362       h->GetXaxis()->SetTitle(fTPCEventHisto->GetAxis(i)->GetTitle());
363       h->GetYaxis()->SetTitle("events");
364       sprintf(title,"%s rec. vertex",fTPCEventHisto->GetAxis(i)->GetTitle());
365       h->SetTitle(title);
366
367       aFolderObj->Add(h);
368   }
369
370   //
371   // Track histograms
372   //
373   for(Int_t i=0; i<9; i++) 
374   {
375       h = (TH1F*)fTPCTrackHisto->Projection(i);
376       sprintf(name,"h_tpc_track_%d",i);
377       h->SetName(name);
378       h->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
379       h->GetYaxis()->SetTitle("tracks");
380       sprintf(title,"%s",fTPCTrackHisto->GetAxis(i)->GetTitle());
381       h->SetTitle(title);
382
383       if(i==7) h->Scale(1,"width");
384       aFolderObj->Add(h);
385   }
386
387   //
388   for(Int_t i=0; i<8; i++) 
389   {
390     for(Int_t j=i+1; j<9; j++) 
391     {
392       h2D = (TH2F*)fTPCTrackHisto->Projection(i,j);
393       sprintf(name,"h_tpc_track_%d_vs_%d",i,j);
394       h2D->SetName(name);
395       h2D->GetXaxis()->SetTitle(fTPCTrackHisto->GetAxis(j)->GetTitle());
396       h2D->GetYaxis()->SetTitle(fTPCTrackHisto->GetAxis(i)->GetTitle());
397       sprintf(title,"%s vs %s",fTPCTrackHisto->GetAxis(j)->GetTitle(),fTPCTrackHisto->GetAxis(i)->GetTitle());
398       h2D->SetTitle(title);
399
400       if(j==7) h2D->SetBit(TH1::kLogX);
401       aFolderObj->Add(h2D);
402     }  
403   }
404
405   // export objects to analysis folder
406   fAnalysisFolder = ExportToFolder(aFolderObj);
407
408   // delete only TObjArray
409   if(aFolderObj) delete aFolderObj;
410 }
411
412 //_____________________________________________________________________________
413 TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array) 
414 {
415   // recreate folder avery time and export objects to new one
416   //
417   AliPerformanceTPC * comp=this;
418   TFolder *folder = comp->GetAnalysisFolder();
419
420   TString name, title;
421   TFolder *newFolder = 0;
422   Int_t i = 0;
423   Int_t size = array->GetSize();
424
425   if(folder) { 
426      // get name and title from old folder
427      name = folder->GetName();  
428      title = folder->GetTitle();  
429
430          // delete old one
431      delete folder;
432
433          // create new one
434      newFolder = CreateFolder(name.Data(),title.Data());
435      newFolder->SetOwner();
436
437          // add objects to folder
438      while(i < size) {
439            newFolder->Add(array->At(i));
440            i++;
441          }
442   }
443
444 return newFolder;
445 }
446
447 //_____________________________________________________________________________
448 Long64_t AliPerformanceTPC::Merge(TCollection* const list) 
449 {
450   // Merge list of objects (needed by PROOF)
451
452   if (!list)
453   return 0;
454
455   if (list->IsEmpty())
456   return 1;
457
458   TIterator* iter = list->MakeIterator();
459   TObject* obj = 0;
460
461   // collection of generated histograms
462   Int_t count=0;
463   while((obj = iter->Next()) != 0) 
464   {
465     AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
466     if (entry == 0) continue; 
467
468     fTPCEventHisto->Add(entry->fTPCEventHisto);
469     fTPCTrackHisto->Add(entry->fTPCTrackHisto);
470
471     count++;
472   }
473
474 return count;
475 }
476
477 //_____________________________________________________________________________
478 TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) { 
479 // create folder for analysed histograms
480 //
481 TFolder *folder = 0;
482   folder = new TFolder(name.Data(),title.Data());
483
484   return folder;
485 }