possiblity to read ESD friends (Jacek)
[u/mrichter/AliRoot.git] / PWG1 / 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   fTPCHisto(0),
64
65   // Cuts 
66   fCutsRC(0),  
67   fCutsMC(0),  
68
69   // histogram folder 
70   fAnalysisFolder(0)
71 {
72   Init();
73 }
74
75 //_____________________________________________________________________________
76 AliPerformanceTPC::AliPerformanceTPC(Char_t* name="AliPerformanceTPC", Char_t* title="AliPerformanceTPC",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
77   AliPerformanceObject(name,title),
78   fTPCHisto(0),
79
80   // Cuts 
81   fCutsRC(0),  
82   fCutsMC(0),  
83
84   // histogram folder 
85   fAnalysisFolder(0)
86 {
87   // named constructor  
88   // 
89   SetAnalysisMode(analysisMode);
90   SetHptGenerator(hptGenerator);
91
92   Init();
93 }
94
95 //_____________________________________________________________________________
96 AliPerformanceTPC::~AliPerformanceTPC()
97 {
98   // destructor
99    
100   if(fTPCHisto) delete fTPCHisto; fTPCHisto=0;     
101   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
102 }
103
104 //_____________________________________________________________________________
105 void AliPerformanceTPC::Init(){
106   //
107   // histogram bining
108   //
109
110   // set pt bins
111   Int_t nPtBins = 50;
112   Double_t ptMin = 1.e-2, ptMax = 10.;
113
114   Double_t *binsPt = 0;
115   if (IsHptGenerator())  { 
116     nPtBins = 100; ptMax = 100.;
117     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
118   } else {
119     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
120   }
121
122   /*
123   Int_t nPtBins = 31;
124   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.};
125   Double_t ptMin = 0., ptMax = 10.; 
126
127   if(IsHptGenerator() == kTRUE) {
128     nPtBins = 100;
129     ptMin = 0.; ptMax = 100.; 
130   }
131   */
132
133   // nTPCClust:chi2PerTPCClust:nTPCClustFindRatio:eta:phi:pt
134   Int_t binsTPCHisto[6]={160,100,100,30,144,nPtBins};
135   Double_t minTPCHisto[6]={0., 0., 0., -1.5, 0., ptMin};
136   Double_t maxTPCHisto[6]={160.,10.,1.2, 1.5, 2.*TMath::Pi(), ptMax};
137
138   fTPCHisto = new THnSparseF("fTPCHisto","nClust:chi2PerClust:nClust/nFindableClust:eta:phi:pt",6,binsTPCHisto,minTPCHisto,maxTPCHisto);
139   fTPCHisto->SetBinEdges(5,binsPt);
140
141   fTPCHisto->GetAxis(0)->SetTitle("nClust");
142   fTPCHisto->GetAxis(1)->SetTitle("chi2PerClust");
143   fTPCHisto->GetAxis(2)->SetTitle("nClust/nFindableClust");
144   fTPCHisto->GetAxis(3)->SetTitle("#eta");
145   fTPCHisto->GetAxis(4)->SetTitle("#phi (rad)");
146   fTPCHisto->GetAxis(5)->SetTitle("p_{T} (GeV/c)");
147   fTPCHisto->Sumw2();
148
149   // Init cuts 
150   if(!fCutsMC) 
151     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
152   if(!fCutsRC) 
153     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
154
155   // init folder
156   fAnalysisFolder = CreateFolder("folderTPC","Analysis Resolution Folder");
157 }
158
159 //_____________________________________________________________________________
160 void AliPerformanceTPC::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
161 {
162   if(!esdTrack) return;
163
164   // Fill TPC only resolution comparison information 
165   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
166   if(!track) return;
167
168   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
169   esdTrack->GetImpactParametersTPC(dca,cov);
170
171   //Float_t q = esdTrack->Charge();
172   Float_t pt = track->Pt();
173   Float_t eta = track->Eta();
174   Float_t phi = track->Phi();
175   Int_t nClust = esdTrack->GetTPCclusters(0);
176   Int_t nFindableClust = esdTrack->GetTPCNclsF();
177
178   Float_t chi2PerCluster = 0.;
179   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
180
181   Float_t clustPerFindClust = 0.;
182   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
183   
184   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) {
185     Double_t vTPCHisto[6] = {nClust,chi2PerCluster,clustPerFindClust,eta,phi,pt};
186     fTPCHisto->Fill(vTPCHisto); 
187   }
188  
189   //
190   // Fill rec vs MC information
191   //
192   if(!stack) return;
193
194 }
195
196 //_____________________________________________________________________________
197 void AliPerformanceTPC::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
198 {
199   // Fill comparison information (TPC+ITS) 
200   AliDebug(AliLog::kWarning, "Warning: Not implemented");
201 }
202  
203 //_____________________________________________________________________________
204 void AliPerformanceTPC::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
205 {
206   // Fill comparison information (constarained parameters) 
207   AliDebug(AliLog::kWarning, "Warning: Not implemented");
208 }
209  
210 //_____________________________________________________________________________
211 void AliPerformanceTPC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
212 {
213   // Process comparison information 
214   //
215   if(!esdEvent) 
216   {
217       AliDebug(AliLog::kError, "esdEvent not available");
218       return;
219   }
220   AliHeader* header = 0;
221   AliGenEventHeader* genHeader = 0;
222   AliStack* stack = 0;
223   TArrayF vtxMC(3);
224   
225   if(bUseMC)
226   {
227     if(!mcEvent) {
228       AliDebug(AliLog::kError, "mcEvent not available");
229       return;
230     }
231
232     // get MC event header
233     header = mcEvent->Header();
234     if (!header) {
235       AliDebug(AliLog::kError, "Header not available");
236       return;
237     }
238     // MC particle stack
239     stack = mcEvent->Stack();
240     if (!stack) {
241       AliDebug(AliLog::kError, "Stack not available");
242       return;
243     }
244
245     // get MC vertex
246     genHeader = header->GenEventHeader();
247     if (!genHeader) {
248       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
249       return;
250     }
251     genHeader->PrimaryVertex(vtxMC);
252
253   } // end bUseMC
254
255   // use ESD friends
256   if(bUseESDfriend) {
257     if(!esdFriend) {
258       AliDebug(AliLog::kError, "esdFriend not available");
259       return;
260     }
261   }
262
263
264
265
266   //  Process events
267   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
268   { 
269     AliESDtrack*track = esdEvent->GetTrack(iTrack);
270     if(!track) continue;
271
272     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
273     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
274     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
275     else {
276       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
277       return;
278     }
279   }
280 }
281
282 //_____________________________________________________________________________
283 void AliPerformanceTPC::Analyse() {
284   //
285   // Analyse comparison information and store output histograms
286   // in the folder "folderTPC"
287   //
288   TH1::AddDirectory(kFALSE);
289   //TH1F *h=0;
290   TH2F *h2D=0;
291   TObjArray *aFolderObj = new TObjArray;
292
293   char name[256];
294   char title[256];
295   for(Int_t i=0; i<5; i++) 
296   {
297     for(Int_t j=i+1; j<6; j++) 
298     {
299       if(j==5) fTPCHisto->GetAxis(5)->SetRangeUser(0.1,10.);
300       h2D = (TH2F*)fTPCHisto->Projection(i,j);
301       sprintf(name,"h_tpc_%d_vs_%d",i,j);
302       h2D->SetName(name);
303       h2D->GetXaxis()->SetTitle(fTPCHisto->GetAxis(j)->GetTitle());
304       h2D->GetYaxis()->SetTitle(fTPCHisto->GetAxis(i)->GetTitle());
305       sprintf(title,"%s vs %s",fTPCHisto->GetAxis(j)->GetTitle(),fTPCHisto->GetAxis(i)->GetTitle());
306       h2D->SetTitle(title);
307
308       if(j==5) h2D->SetBit(TH1::kLogX);
309       aFolderObj->Add(h2D);
310     }  
311   }
312
313   // export objects to analysis folder
314   fAnalysisFolder = ExportToFolder(aFolderObj);
315
316   // delete only TObjArray
317   if(aFolderObj) delete aFolderObj;
318 }
319
320 //_____________________________________________________________________________
321 TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array) 
322 {
323   // recreate folder avery time and export objects to new one
324   //
325   AliPerformanceTPC * comp=this;
326   TFolder *folder = comp->GetAnalysisFolder();
327
328   TString name, title;
329   TFolder *newFolder = 0;
330   Int_t i = 0;
331   Int_t size = array->GetSize();
332
333   if(folder) { 
334      // get name and title from old folder
335      name = folder->GetName();  
336      title = folder->GetTitle();  
337
338          // delete old one
339      delete folder;
340
341          // create new one
342      newFolder = CreateFolder(name.Data(),title.Data());
343      newFolder->SetOwner();
344
345          // add objects to folder
346      while(i < size) {
347            newFolder->Add(array->At(i));
348            i++;
349          }
350   }
351
352 return newFolder;
353 }
354
355 //_____________________________________________________________________________
356 Long64_t AliPerformanceTPC::Merge(TCollection* const list) 
357 {
358   // Merge list of objects (needed by PROOF)
359
360   if (!list)
361   return 0;
362
363   if (list->IsEmpty())
364   return 1;
365
366   TIterator* iter = list->MakeIterator();
367   TObject* obj = 0;
368
369   // collection of generated histograms
370   Int_t count=0;
371   while((obj = iter->Next()) != 0) 
372   {
373   AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
374   if (entry == 0) continue; 
375
376   fTPCHisto->Add(entry->fTPCHisto);
377
378   count++;
379   }
380
381 return count;
382 }
383
384 //_____________________________________________________________________________
385 TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) { 
386 // create folder for analysed histograms
387 //
388 TFolder *folder = 0;
389   folder = new TFolder(name.Data(),title.Data());
390
391   return folder;
392 }