Fixes for cmake
[u/mrichter/AliRoot.git] / PWG1 / AliPerformanceRes.cxx
1 //------------------------------------------------------------------------------
2 // Implementation of AliPerformanceRes 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 AliPerformanceRes.
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   AliPerformanceRes * compObj = (AliPerformanceRes*)coutput->FindObject("AliPerformanceRes");
21  
22   // analyse comparison data
23   compObj->Analyse();
24
25   // the output histograms/graphs will be stored in the folder "folderRes" 
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_Res.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 "TF1.h"
41
42 #include "AliPerformanceRes.h" 
43 #include "AliESDEvent.h" 
44 #include "AliESDVertex.h"
45 #include "AliESDtrack.h"
46 #include "AliESDfriendTrack.h"
47 #include "AliESDfriend.h"
48 #include "AliLog.h" 
49 #include "AliMCEvent.h" 
50 #include "AliMCParticle.h" 
51 #include "AliHeader.h" 
52 #include "AliGenEventHeader.h" 
53 #include "AliStack.h" 
54 #include "AliMCInfoCuts.h" 
55 #include "AliRecInfoCuts.h" 
56 #include "AliTracker.h" 
57 #include "AliTreeDraw.h" 
58
59 using namespace std;
60
61 ClassImp(AliPerformanceRes)
62
63 //_____________________________________________________________________________
64 AliPerformanceRes::AliPerformanceRes():
65   AliPerformanceObject("AliPerformanceRes"),
66   fResolHisto(0),
67   fPullHisto(0),
68
69   // Cuts 
70   fCutsRC(0),  
71   fCutsMC(0),  
72
73   // histogram folder 
74   fAnalysisFolder(0)
75 {
76   Init();
77 }
78
79 //_____________________________________________________________________________
80 AliPerformanceRes::AliPerformanceRes(Char_t* name="AliPerformanceRes", Char_t* title="AliPerformanceRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
81   AliPerformanceObject(name,title),
82   fResolHisto(0),
83   fPullHisto(0),
84
85   // Cuts 
86   fCutsRC(0),  
87   fCutsMC(0),  
88
89   // histogram folder 
90   fAnalysisFolder(0)
91 {
92   // named constructor  
93   // 
94   SetAnalysisMode(analysisMode);
95   SetHptGenerator(hptGenerator);
96
97   Init();
98 }
99
100 //_____________________________________________________________________________
101 AliPerformanceRes::~AliPerformanceRes()
102 {
103   // destructor
104    
105   if(fResolHisto) delete fResolHisto; fResolHisto=0;     
106   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;     
107   
108   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
109 }
110
111 //_____________________________________________________________________________
112 void AliPerformanceRes::Init(){
113
114   //
115   // histogram bining
116   //
117
118   // set pt bins
119   Int_t nPtBins = 50;
120   Double_t ptMin = 1.e-2, ptMax = 10.;
121
122   Double_t *binsPt = 0;
123   if (IsHptGenerator())  { 
124     nPtBins = 100; ptMax = 100.;
125     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
126   } else {
127     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
128   }
129
130   Double_t yMin = -0.02, yMax = 0.02;
131   Double_t zMin = -12.0, zMax = 12.0;
132   if(GetAnalysisMode() == 3) { // TrackRef coordinate system
133     yMin = -100.; yMax = 100.; 
134     zMin = -100.; zMax = 100.; 
135   }
136
137   // res_y:res_z:res_phi,res_lambda:res_pt:y:z:eta:phi:pt
138   Int_t binsResolHisto[10]={100,100,100,100,100,25,50,144,30,nPtBins};
139   Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin, 0., -1.5, ptMin};
140   Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 2.*TMath::Pi(), 1.5, ptMax};
141
142   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_pt:y:z:phi:eta:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
143   fResolHisto->SetBinEdges(9,binsPt);
144
145   fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)");
146   fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)");
147   fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)");
148   fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)");
149   fResolHisto->GetAxis(4)->SetTitle("(p_{T}/p_{Tmc}-1)");
150   fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
151   fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
152   fResolHisto->GetAxis(7)->SetTitle("#phi_{mc} (rad)");
153   fResolHisto->GetAxis(8)->SetTitle("#eta_{mc}");
154   fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
155   fResolHisto->Sumw2();
156
157   ////pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt
158   //Int_t binsPullHisto[10]={100,100,100,100,100,50,50,30,144,nPtBins};
159   //Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1.5, 0., ptMin};
160   //Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1.5, 2.*TMath::Pi(),ptMax};
161   //fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
162
163   //pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt
164   Int_t binsPullHisto[10]={100,100,100,100,100,50,50,50,50,nPtBins};
165   Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1., -2.0, ptMin};
166   Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1., 2.0, ptMax};
167   fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_y:pull_z:pull_snp:pull_tgl:pull_1pt:y:z:snp:tgl:1pt",10,binsPullHisto,minPullHisto,maxPullHisto);
168
169   /*
170   if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,bins1Pt);
171   fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
172   fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
173   fPullHisto->GetAxis(2)->SetTitle("(#phi-#phi_{mc})/#sigma");
174   fPullHisto->GetAxis(3)->SetTitle("(#lambda-#lambda_{mc})/#sigma");
175   fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
176   fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
177   fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
178   fPullHisto->GetAxis(7)->SetTitle("#eta_{mc}");
179   fPullHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
180   fPullHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
181   fPullHisto->Sumw2();
182   */
183
184   fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
185   fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
186   fPullHisto->GetAxis(2)->SetTitle("(sin#phi-sin#phi_{mc})/#sigma");
187   fPullHisto->GetAxis(3)->SetTitle("(tan#lambda-tan#lambda_{mc})/#sigma");
188   fPullHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)/#sigma");
189   fPullHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
190   fPullHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
191   fPullHisto->GetAxis(7)->SetTitle("sin#phi_{mc}");
192   fPullHisto->GetAxis(8)->SetTitle("tan#lambda_{mc}");
193   fPullHisto->GetAxis(9)->SetTitle("1/p_{Tmc} (GeV/c)^{-1}");
194   fPullHisto->Sumw2();
195
196   // Init cuts 
197   if(!fCutsMC) 
198     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
199   if(!fCutsRC) 
200     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
201
202   // init folder
203   fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder");
204 }
205
206 //_____________________________________________________________________________
207 void AliPerformanceRes::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
208 {
209   if(!esdTrack) return;
210
211   // Fill TPC only resolution comparison information 
212   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
213   if(!track) return;
214
215   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
216   esdTrack->GetImpactParametersTPC(dca,cov);
217  
218   //
219   // Fill rec vs MC information
220   //
221   if(!stack) return;
222   Int_t label = TMath::Abs(esdTrack->GetLabel()); 
223   TParticle* particle = stack->Particle(label);
224   if(!particle) return;
225   if(!particle->GetPDG()) return;
226   if(particle->GetPDG()->Charge()==0) return;
227   //printf("charge %d \n",particle->GetPDG()->Charge());
228   
229   // Only 5 charged particle species (e,mu,pi,K,p)
230   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
231
232   // exclude electrons
233   if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
234
235   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
236   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
237
238   Float_t mceta =  particle->Eta();
239   Float_t mcphi =  particle->Phi();
240   if(mcphi<0) mcphi += 2.*TMath::Pi();
241   Float_t mcpt = particle->Pt();
242   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
243   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
244
245   // nb. TPC clusters cut
246   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
247
248   // select primaries
249   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
250   { 
251     if(mcpt == 0) return;
252     
253     deltaYTPC= track->GetY()-particle->Vy();
254     deltaZTPC = track->GetZ()-particle->Vz();
255     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
256     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
257     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
258     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
259
260     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
261     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
262  
263     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2()); 
264     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
265     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
266     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
267
268     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
269     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); 
270     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
271     else pull1PtTPC = 0.;
272
273     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
274     fResolHisto->Fill(vResolHisto);
275
276     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
277     fPullHisto->Fill(vPullHisto);
278   }
279 }
280
281 //_____________________________________________________________________________
282 void AliPerformanceRes::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack)
283 {
284   // Fill resolution comparison information (TPC+ITS)
285   if(!esdTrack) return;
286
287   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
288   esdTrack->GetImpactParameters(dca,cov);
289  
290   //
291   // Fill rec vs MC information
292   //
293   if(!stack) return;
294
295   Int_t label = TMath::Abs(esdTrack->GetLabel()); 
296   TParticle* particle = stack->Particle(label);
297   if(!particle) return;
298   if(particle->GetPDG()->Charge()==0) return;
299   //printf("charge %d \n",particle->GetPDG()->Charge());
300
301
302   // Only 5 charged particle species (e,mu,pi,K,p)
303   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
304
305   // exclude electrons
306   if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
307
308   Float_t mceta =  particle->Eta();
309   Float_t mcphi =  particle->Phi();
310   if(mcphi<0) mcphi += 2.*TMath::Pi();
311   Float_t mcpt = particle->Pt();
312   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
313   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
314
315   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
316   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters  
317   Int_t clusterITS[200];
318   if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return;  // min. nb. ITS clusters
319
320   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
321   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
322
323   // select primaries
324   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
325   { 
326     if(mcpt == 0) return;
327     
328     deltaYTPC= esdTrack->GetY()-particle->Vy();
329     deltaZTPC = esdTrack->GetZ()-particle->Vz();
330     deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
331     deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
332     //delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
333     deltaPtTPC = (esdTrack->Pt()-mcpt) / mcpt;
334
335     pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
336     pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
337  
338     //Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2()); 
339     //Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
340     pullPhiTPC = (esdTrack->GetSnp() - mcsnp) / TMath::Sqrt(esdTrack->GetSigmaSnp2());
341     pullLambdaTPC = (esdTrack->GetTgl() - mctgl) / TMath::Sqrt(esdTrack->GetSigmaTgl2());
342
343     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
344     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2()); 
345     if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
346     else pull1PtTPC = 0.;
347
348     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
349     fResolHisto->Fill(vResolHisto);
350
351     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
352     fPullHisto->Fill(vPullHisto);
353
354    
355     /*
356     deltaYTPC= esdTrack->GetY()-particle->Vy();
357     deltaZTPC = esdTrack->GetZ()-particle->Vz();
358     deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
359     deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
360     delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
361
362     pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
363     pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
364
365     Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2()); 
366     Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
367     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
368     pullPhiTPC = deltaPhiTPC / sigma_phi;
369
370     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
371     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2()); 
372     if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
373     else pull1PtTPC = 0.;
374
375     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
376     fResolHisto->Fill(vResolHisto);
377
378     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
379     fPullHisto->Fill(vPullHisto);
380     */
381   }
382 }
383
384 //_____________________________________________________________________________
385 void AliPerformanceRes::ProcessConstrained(AliStack* const stack, AliESDtrack *const esdTrack)
386 {
387   // Fill resolution comparison information (constarained parameters) 
388   //
389   if(!esdTrack) return;
390
391   const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
392   if(!track) return;
393
394   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
395   esdTrack->GetImpactParameters(dca,cov);
396  
397   //
398   // Fill rec vs MC information
399   //
400   if(!stack) return;
401
402   Int_t label = TMath::Abs(esdTrack->GetLabel()); 
403   TParticle* particle = stack->Particle(label);
404   if(!particle) return;
405   if(particle->GetPDG()->Charge()==0) return;
406   //printf("charge %d \n",particle->GetPDG()->Charge());
407
408   // Only 5 charged particle species (e,mu,pi,K,p)
409   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
410
411   // exclude electrons
412   if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
413
414   Float_t mceta =  particle->Eta();
415   Float_t mcphi =  particle->Phi();
416   if(mcphi<0) mcphi += 2.*TMath::Pi();
417   Float_t mcpt = particle->Pt();
418   Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
419   Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
420
421   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
422   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
423
424   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
425   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
426
427   // select primaries
428   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
429   { 
430
431     if(mcpt == 0) return;
432     
433     deltaYTPC= track->GetY()-particle->Vy();
434     deltaZTPC = track->GetZ()-particle->Vz();
435     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
436     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
437     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
438     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
439
440     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
441     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
442  
443     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2()); 
444     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
445     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
446     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
447
448     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
449     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); 
450     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
451     else pull1PtTPC = 0.;
452
453     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
454     fResolHisto->Fill(vResolHisto);
455
456     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
457     fPullHisto->Fill(vPullHisto);
458
459     /*
460
461     deltaYTPC= track->GetY()-particle->Vy();
462     deltaZTPC = track->GetZ()-particle->Vz();
463     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
464     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
465     delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
466
467     pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
468     pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
469
470     Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2()); 
471     Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
472     pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
473     pullPhiTPC = deltaPhiTPC / sigma_phi;
474
475     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
476     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); 
477
478     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
479     else pull1PtTPC = 0.;
480
481     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
482     fResolHisto->Fill(vResolHisto);
483
484     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
485     fPullHisto->Fill(vPullHisto);
486
487     */
488   }
489 }
490  
491 //_____________________________________________________________________________
492 void AliPerformanceRes::ProcessInnerTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack)
493 {
494   //
495   // Fill resolution comparison information (inner params at TPC reference point) 
496   //
497   if(!esdTrack) return;
498
499   const AliExternalTrackParam * innerParam = esdTrack->GetInnerParam();
500   if(!innerParam) return;
501
502   // create new AliExternalTrackParam
503   AliExternalTrackParam *track = new AliExternalTrackParam(*innerParam);
504   if(!track) return;
505
506   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
507   esdTrack->GetImpactParametersTPC(dca,cov);
508  
509   //
510   // Fill rec vs MC information
511   //
512   if(!mcEvent) return;
513
514   Int_t label = TMath::Abs(esdTrack->GetLabel()); 
515   AliMCParticle *mcParticle = mcEvent->GetTrack(label);
516   if(!mcParticle) return;
517
518   // get the first TPC track reference
519   AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle);
520   if(!ref0) return;
521
522   // Only 5 charged particle species (e,mu,pi,K,p)
523   TParticle *particle = mcParticle->Particle();
524   if(!particle) return;
525
526   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
527
528   // exclude electrons
529   if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
530
531   // calculate alpha angle
532   Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
533   Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
534   //if (alpha<0) alpha += TMath::TwoPi();
535
536   // rotate inner track to local coordinate system
537   // and propagate to the radius of the first track referenco of TPC
538   Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
539   Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
540   if(!isOK) return;
541
542   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
543   Float_t mcphi =  ref0->Phi();
544   if(mcphi<0) mcphi += 2.*TMath::Pi();
545   Float_t mcpt = ref0->Pt();
546   Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
547   Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
548
549   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
550   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
551
552   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
553   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
554
555   // select primaries
556   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
557   { 
558     if(mcpt == 0) return;
559     
560     deltaYTPC= track->GetY(); // already rotated
561     deltaZTPC = track->GetZ()-ref0->Z();
562     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
563     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
564     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
565     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
566
567     pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2());
568     pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
569  
570     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2()); 
571     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
572     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
573     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
574
575     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
576     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); 
577     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
578     else pull1PtTPC = 0.;
579
580     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
581     fResolHisto->Fill(vResolHisto);
582
583     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
584     fPullHisto->Fill(vPullHisto);
585   }
586
587   if(track) delete track;
588 }
589
590 //_____________________________________________________________________________
591 void AliPerformanceRes::ProcessOuterTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack, AliESDfriendTrack *const friendTrack)
592 {
593   //
594   // Fill resolution comparison information (outer params at TPC reference point) 
595   //
596   if(!friendTrack) return;
597
598   const AliExternalTrackParam * outerParam = friendTrack->GetTPCOut();
599   if(!outerParam) return;
600
601   // create new AliExternalTrackParam
602   AliExternalTrackParam *track = new AliExternalTrackParam(*outerParam);
603   if(!track) return;
604
605   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
606   esdTrack->GetImpactParametersTPC(dca,cov);
607  
608   //
609   // Fill rec vs MC information
610   //
611   if(!mcEvent) return;
612
613   Int_t label = TMath::Abs(esdTrack->GetLabel()); 
614   AliMCParticle *mcParticle = mcEvent->GetTrack(label);
615   if(!mcParticle) return;
616
617   // get the last TPC track reference
618   AliTrackReference *ref0 = GetLastTPCTrackRef(mcParticle);
619   if(!ref0) return;
620
621   // Only 5 charged particle species (e,mu,pi,K,p)
622   TParticle *particle = mcParticle->Particle();
623   if(!particle) return;
624
625   if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
626
627   // exclude electrons
628   if (fCutsMC->GetEP()==TMath::Abs(particle->GetPdgCode())) return;
629
630   // calculate alpha angle
631   Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
632   Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
633   //if (alpha<0) alpha += TMath::TwoPi();
634
635   // rotate outer track to local coordinate system
636   // and propagate to the radius of the last track reference of TPC
637   Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
638   Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
639   if(!isOK) return;
640
641   Float_t mceta =  -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
642   Float_t mcphi =  ref0->Phi();
643   if(mcphi<0) mcphi += 2.*TMath::Pi();
644   Float_t mcpt = ref0->Pt();
645   Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
646   Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
647
648   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
649   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
650
651   Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
652   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
653
654   // select primaries
655   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
656   { 
657     if(mcpt == 0) return;
658     
659     deltaYTPC= track->GetY(); // already rotated
660     deltaZTPC = track->GetZ()-ref0->Z();
661     deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
662     deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
663     //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
664     deltaPtTPC = (track->Pt()-mcpt) / mcpt;
665
666     pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2());
667     pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
668  
669     //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2()); 
670     //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
671     pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
672     pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
673
674     //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
675     //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); 
676     if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
677     else pull1PtTPC = 0.;
678
679     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
680     fResolHisto->Fill(vResolHisto);
681
682     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
683     fPullHisto->Fill(vPullHisto);
684   }
685
686   if(track) delete track;
687 }
688
689 //_____________________________________________________________________________
690 AliTrackReference * AliPerformanceRes::GetFirstTPCTrackRef(AliMCParticle *mcParticle) 
691 {
692   // return first TPC track reference 
693   if(!mcParticle) return 0;
694
695   // find first track reference 
696   // check direction to select proper reference point for looping tracks
697   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
698   AliTrackReference *ref = 0;
699   AliTrackReference *refIn = 0;
700   for (Int_t iref = 0; iref < nTrackRef; iref++) { 
701     ref = mcParticle->GetTrackReference(iref);
702     if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
703     {
704       Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
705       if(dir < 0.) break;
706
707       refIn = ref;
708       break;
709     }
710   }
711
712 return refIn;
713 }
714
715 //_____________________________________________________________________________
716 AliTrackReference * AliPerformanceRes::GetLastTPCTrackRef(AliMCParticle *mcParticle) 
717 {
718   // return last TPC track reference 
719   if(!mcParticle) return 0;
720
721   // find last track reference 
722   // check direction to select proper reference point for looping tracks
723   Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
724   AliTrackReference *ref = 0;
725   AliTrackReference *refOut = 0;
726   for (Int_t iref = 0; iref < nTrackRef; iref++) { 
727     ref = mcParticle->GetTrackReference(iref);
728     if(ref && (ref->DetectorId()==AliTrackReference::kTPC)) { 
729       Float_t dir=ref->X()*ref->Px()+ref->Y()*ref->Py();
730       if(dir< 0.0) break;
731       refOut = ref;
732     }
733   }
734
735 return refOut;
736 }
737
738 //_____________________________________________________________________________
739 void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
740 {
741   // Process comparison information 
742   //
743   if(!esdEvent) 
744   {
745       AliDebug(AliLog::kError, "esdEvent not available");
746       return;
747   }
748   AliHeader* header = 0;
749   AliGenEventHeader* genHeader = 0;
750   AliStack* stack = 0;
751   TArrayF vtxMC(3);
752   
753   if(bUseMC)
754   {
755     if(!mcEvent) {
756       AliDebug(AliLog::kError, "mcEvent not available");
757       return;
758     }
759
760     // get MC event header
761     header = mcEvent->Header();
762     if (!header) {
763       AliDebug(AliLog::kError, "Header not available");
764       return;
765     }
766     // MC particle stack
767     stack = mcEvent->Stack();
768     if (!stack) {
769       AliDebug(AliLog::kError, "Stack not available");
770       return;
771     }
772
773     // get MC vertex
774     genHeader = header->GenEventHeader();
775     if (!genHeader) {
776       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
777       return;
778     }
779     genHeader->PrimaryVertex(vtxMC);
780
781   } // end bUseMC
782   
783   // use ESD friends
784   if(bUseESDfriend) {
785     if(!esdFriend) {
786       AliDebug(AliLog::kError, "esdFriend not available");
787       return;
788     }
789   }
790
791   //  Process events
792   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
793   { 
794     AliESDtrack *track = esdEvent->GetTrack(iTrack);
795     if(!track) continue;
796
797     AliESDfriendTrack *friendTrack=0;
798     if(bUseESDfriend) {
799       friendTrack=esdFriend->GetTrack(iTrack);
800       if(!friendTrack) continue;
801     }
802
803     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
804     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
805     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
806     else if(GetAnalysisMode() == 3) ProcessInnerTPC(mcEvent,track);
807     else if(GetAnalysisMode() == 4) ProcessOuterTPC(mcEvent,track,friendTrack);
808     else {
809       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
810       return;
811     }
812   }
813 }
814
815 //_____________________________________________________________________________
816 TH1F* AliPerformanceRes::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
817   // Create resolution histograms
818  
819   TH1F *hisr, *hism;
820   if (!gPad) new TCanvas;
821   hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
822   if (type) return hism;
823   else 
824     return hisr;
825 }
826
827 //_____________________________________________________________________________
828 void AliPerformanceRes::Analyse() {
829   // Analyse comparison information and store output histograms
830   // in the folder "folderRes"
831   //
832   TH1::AddDirectory(kFALSE);
833   TH1F *h=0;
834   TH2F *h2D=0;
835   TObjArray *aFolderObj = new TObjArray;
836
837   // write results in the folder 
838   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
839   c->cd();
840
841   char name[256];
842   char title[256];
843
844   for(Int_t i=0; i<5; i++) 
845   {
846     for(Int_t j=5; j<10; j++) 
847     {
848       if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
849       fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
850       //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.,0.9); // eta window
851       //fResolHisto->GetAxis(9)->SetRangeUser(0.16,3.); // pt window
852       if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
853
854       h2D = (TH2F*)fResolHisto->Projection(i,j);
855
856       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
857       sprintf(name,"h_res_%d_vs_%d",i,j);
858       h->SetName(name);
859
860       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
861       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
862       h->GetYaxis()->SetTitle(title);
863       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
864       h->SetTitle(title);
865
866       if(j==9) h->SetBit(TH1::kLogX);    
867       aFolderObj->Add(h);
868
869       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
870       //h = (TH1F*)arr->At(1);
871       sprintf(name,"h_mean_res_%d_vs_%d",i,j);
872       h->SetName(name);
873
874       h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
875       sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
876       h->GetYaxis()->SetTitle(title);
877
878       sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
879       h->SetTitle(title);
880
881       if(j==9) h->SetBit(TH1::kLogX);    
882       aFolderObj->Add(h);
883
884       fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
885       fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
886
887       //
888       if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-1.0,0.99); // theta window
889       else  fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.5); // theta window
890       fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);  // pt threshold
891       if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
892
893       h2D = (TH2F*)fPullHisto->Projection(i,j);
894
895       h = AliPerformanceRes::MakeResol(h2D,1,0,100);
896       sprintf(name,"h_pull_%d_vs_%d",i,j);
897       h->SetName(name);
898
899       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
900       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
901       h->GetYaxis()->SetTitle(title);
902       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
903       h->SetTitle(title);
904
905       //if(j==9) h->SetBit(TH1::kLogX);    
906       aFolderObj->Add(h);
907
908       h = AliPerformanceRes::MakeResol(h2D,1,1,100);
909       sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
910       h->SetName(name);
911
912       h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
913       sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
914       h->GetYaxis()->SetTitle(title);
915       sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
916       h->SetTitle(title);
917
918       //if(j==9) h->SetBit(TH1::kLogX);    
919       aFolderObj->Add(h);
920     }
921   }
922
923   // export objects to analysis folder
924   fAnalysisFolder = ExportToFolder(aFolderObj);
925
926   // delete only TObjArray
927   if(aFolderObj) delete aFolderObj;
928 }
929
930 //_____________________________________________________________________________
931 TFolder* AliPerformanceRes::ExportToFolder(TObjArray * array) 
932 {
933   // recreate folder avery time and export objects to new one
934   //
935   AliPerformanceRes * comp=this;
936   TFolder *folder = comp->GetAnalysisFolder();
937
938   TString name, title;
939   TFolder *newFolder = 0;
940   Int_t i = 0;
941   Int_t size = array->GetSize();
942
943   if(folder) { 
944      // get name and title from old folder
945      name = folder->GetName();  
946      title = folder->GetTitle();  
947
948          // delete old one
949      delete folder;
950
951          // create new one
952      newFolder = CreateFolder(name.Data(),title.Data());
953      newFolder->SetOwner();
954
955          // add objects to folder
956      while(i < size) {
957            newFolder->Add(array->At(i));
958            i++;
959          }
960   }
961
962 return newFolder;
963 }
964
965 //_____________________________________________________________________________
966 Long64_t AliPerformanceRes::Merge(TCollection* const list) 
967 {
968   // Merge list of objects (needed by PROOF)
969
970   if (!list)
971   return 0;
972
973   if (list->IsEmpty())
974   return 1;
975
976   TIterator* iter = list->MakeIterator();
977   TObject* obj = 0;
978
979   // collection of generated histograms
980   Int_t count=0;
981   while((obj = iter->Next()) != 0) 
982   {
983   AliPerformanceRes* entry = dynamic_cast<AliPerformanceRes*>(obj);
984   if (entry == 0) continue; 
985
986   fResolHisto->Add(entry->fResolHisto);
987   fPullHisto->Add(entry->fPullHisto);
988
989   count++;
990   }
991
992 return count;
993 }
994
995 //_____________________________________________________________________________
996 TFolder* AliPerformanceRes::CreateFolder(TString name,TString title) { 
997 // create folder for analysed histograms
998 //
999 TFolder *folder = 0;
1000   folder = new TFolder(name.Data(),title.Data());
1001
1002   return folder;
1003 }