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