]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG1/TPC/AliComparisonRes.cxx
Summary and trend extraction for TPC with modified AliPerformanceTPC
[u/mrichter/AliRoot.git] / PWG1 / TPC / AliComparisonRes.cxx
1 //------------------------------------------------------------------------------
2 // Implementation of AliComparisonRes 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 AliComparisonRes.
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   //AliComparisonRes * compObj = (AliComparisonRes*)f.Get("AliComparisonRes");
21   AliComparisonRes * compObj = (AliComparisonRes*)cOutput->FindObject("AliComparisonRes");
22  
23   // analyse comparison data
24   compObj->Analyse();
25
26   // the output histograms/graphs will be stored in the folder "folderRes" 
27   compObj->GetAnalysisFolder()->ls("*");
28
29   // user can save whole comparison object (or only folder with anlysed histograms) 
30   // in the seperate output file (e.g.)
31   TFile fout("Analysed_Res.root","recreate");
32   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
33   fout.Close();
34
35 */
36
37 #include "TCanvas.h"
38 #include "TH1.h"
39 #include "TAxis.h"
40
41 #include "AliComparisonRes.h" 
42 #include "AliESDRecInfo.h" 
43 #include "AliESDVertex.h"
44 #include "AliESDtrack.h"
45 #include "AliLog.h" 
46 #include "AliMCInfo.h" 
47 #include "AliMCInfoCuts.h" 
48 #include "AliRecInfoCuts.h" 
49 #include "AliTracker.h" 
50 #include "AliTreeDraw.h" 
51
52 using namespace std;
53
54 ClassImp(AliComparisonRes)
55
56 //_____________________________________________________________________________
57 AliComparisonRes::AliComparisonRes():
58   AliComparisonObject("AliComparisonRes"),
59   fResolHisto(0),
60   fPullHisto(0),
61
62   // Cuts 
63   fCutsRC(0),  
64   fCutsMC(0),  
65
66   // histogram folder 
67   fAnalysisFolder(0)
68 {
69   //Init();
70 }
71
72 //_____________________________________________________________________________
73 AliComparisonRes::AliComparisonRes(Char_t* name="AliComparisonRes", Char_t* title="AliComparisonRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
74 //AliComparisonRes::AliComparisonRes(Char_t* name, Char_t* title,Int_t analysisMode,Bool_t hptGenerator):
75   AliComparisonObject(name,title),
76   fResolHisto(0),
77   fPullHisto(0),
78
79   // Cuts 
80   fCutsRC(0),  
81   fCutsMC(0),  
82
83   // histogram folder 
84   fAnalysisFolder(0)
85 {
86   // named constructor  
87   // 
88   SetAnalysisMode(analysisMode);
89   SetHptGenerator(hptGenerator);
90
91   Init();
92 }
93
94 //_____________________________________________________________________________
95 AliComparisonRes::~AliComparisonRes()
96 {
97   // destructor
98    
99   if(fResolHisto) delete fResolHisto; fResolHisto=0;     
100   if(fPullHisto)  delete fPullHisto;  fPullHisto=0;     
101
102   if(fCutsRC) delete fCutsRC; fCutsRC=0;
103   if(fCutsMC) delete fCutsMC; fCutsMC=0;
104
105   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
106 }
107
108 //_____________________________________________________________________________
109 void AliComparisonRes::Init(){
110
111   //
112   // histogram bining
113   //
114
115   // set pt bins
116    Int_t nPtBins = 31;
117    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.};
118   Double_t ptMin = 0., ptMax = 10.; 
119
120
121   if(IsHptGenerator() == kTRUE) {
122     nPtBins = 100;
123     ptMin = 0.; ptMax = 100.; 
124   }
125
126   // res_y:res_z:res_phi,res_lambda:res_1pt:y:z:eta:phi:pt
127   Int_t binsResolHisto[10]={100,100,100,100,100,50,100,30,90,nPtBins};
128   Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2,-0.03,-20.,-1.5,0.,ptMin};
129   Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, 0.03, 20., 1.5,2.*TMath::Pi(), ptMax};
130
131
132   fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_1pt:y:z:eta:phi:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
133   if(!IsHptGenerator()) fResolHisto->SetBinEdges(9,binsPt);
134
135   fResolHisto->GetAxis(0)->SetTitle("res_y (cm)");
136   fResolHisto->GetAxis(1)->SetTitle("res_z (cm)");
137   fResolHisto->GetAxis(2)->SetTitle("res_phi (rad)");
138   fResolHisto->GetAxis(3)->SetTitle("res_lambda (rad)");
139   fResolHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}");
140   fResolHisto->GetAxis(5)->SetTitle("y (cm)");
141   fResolHisto->GetAxis(6)->SetTitle("z (cm)");
142   fResolHisto->GetAxis(7)->SetTitle("eta");
143   fResolHisto->GetAxis(8)->SetTitle("phi (rad)");
144   fResolHisto->GetAxis(9)->SetTitle("pt (GeV/c)");
145   fResolHisto->Sumw2();
146
147   //pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt
148   Int_t binsPullHisto[10]={100,100,100,100,100,50,100,30,90,nPtBins};
149   Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,-0.03, -20.,-1.5, 0., ptMin};
150   Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., 0.03, 20., 1.5, 2.*TMath::Pi(),ptMax};
151
152   fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
153   if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,binsPt);
154   
155   fPullHisto->GetAxis(0)->SetTitle("res_y (cm)");
156   fPullHisto->GetAxis(1)->SetTitle("res_z (cm)");
157   fPullHisto->GetAxis(2)->SetTitle("res_phi (rad)");
158   fPullHisto->GetAxis(3)->SetTitle("res_lambda (rad)");
159   fPullHisto->GetAxis(4)->SetTitle("res_1pt (GeV/c)^{-1}");
160   fPullHisto->GetAxis(5)->SetTitle("y (rad)");
161   fPullHisto->GetAxis(6)->SetTitle("z (rad)");
162   fPullHisto->GetAxis(7)->SetTitle("eta");
163   fPullHisto->GetAxis(8)->SetTitle("phi (rad)");
164   fPullHisto->GetAxis(9)->SetTitle("pt (GeV/c)");
165   fPullHisto->Sumw2();
166
167   // Init cuts 
168   if(!fCutsMC) 
169     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
170   if(!fCutsRC) 
171     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
172
173   // init folder
174   fAnalysisFolder = CreateFolder("folderRes","Analysis Resolution Folder");
175 }
176
177 //_____________________________________________________________________________
178 void AliComparisonRes::ProcessTPC(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
179 {
180   // Fill resolution comparison information 
181   AliExternalTrackParam *track = 0;
182   //Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
183   Double_t kMaxD      = 123456.0; // max distance
184   AliESDVertex vertexMC;  // MC primary vertex
185
186   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
187
188   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
189   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
190   //Float_t delta1Pt2TPC, deltaY1PtTPC, deltaZ1PtTPC, deltaPhi1PtTPC, deltaTheta1PtTPC; 
191
192   Float_t mceta =  infoMC->GetParticle().Eta();
193   Float_t mcphi =  infoMC->GetParticle().Phi();
194   if(mcphi<0) mcphi += 2.*TMath::Pi();
195   Float_t mcpt = infoMC->GetParticle().Pt();
196
197   // distance to Prim. vertex 
198   const Double_t* dv = infoMC->GetVDist(); 
199   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
200
201   // Only 5 charged particle species (e,mu,pi,K,p)
202   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
203   if (!isPrim) return;
204   //if (infoRC->GetStatus(1)!=3) return; // TPC refit
205   if (!infoRC->GetESDtrack()) return;  
206   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
207
208   // calculate and set prim. vertex
209   vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] );
210   vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] );
211   vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] );
212
213   // calculate track parameters at vertex
214   const AliExternalTrackParam *innerTPC =  0;
215   if ((innerTPC = infoRC->GetESDtrack()->GetTPCInnerParam()) != 0)
216   {
217     if ((track = new AliExternalTrackParam(*infoRC->GetESDtrack()->GetTPCInnerParam())) != 0 )
218     {
219       //Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov);
220       Double_t field[3];  track->GetBxByBz(field); 
221       Bool_t bDCAStatus = track->PropagateToDCABxByBz(&vertexMC,field,kMaxD,dca,cov);
222
223       // Fill parametrisation histograms (only TPC track)
224       if(bDCAStatus) 
225       {
226           // select primaries
227           if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
228           { 
229
230             deltaYTPC= track->GetY()-infoMC->GetParticle().Vy();
231             deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz();
232             deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
233             deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
234             delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
235
236             pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2());
237             pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2());
238             pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
239             pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); 
240             pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
241
242             Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
243             fResolHisto->Fill(vResolHisto);
244
245             Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
246             fPullHisto->Fill(vPullHisto);
247           }
248
249           /*
250             delta1Pt2TPC = (1/mcpt-innerTPC->OneOverPt())/TMath::Power(1+1/mcpt,2);       
251             deltaY1PtTPC= (infoMC->GetParticle().Vy()-innerTPC->GetY()) / (0.2+1/mcpt);
252             deltaZ1PtTPC = (infoMC->GetParticle().Vz()-innerTPC->GetZ()) / (0.2+1/mcpt);
253             deltaPhi1PtTPC = deltaPhiTPC   / (0.1+1/mcpt);
254             deltaTheta1PtTPC = deltaLambdaTPC / (0.1+1/mcpt);
255
256             fPhiEtaPtTPC->Fill(TMath::ATan2(innerTPC->Py(),innerTPC->Px()), innerTPC->Eta(), innerTPC->Pt()); 
257
258             f1Pt2ResolS1PtTPC->Fill(s1mcpt,delta1Pt2TPC);
259             fYResolS1PtTPC->Fill(s1mcpt,deltaY1PtTPC);
260             fZResolS1PtTPC->Fill(s1mcpt,deltaZ1PtTPC);
261             fPhiResolS1PtTPC->Fill(s1mcpt,deltaPhi1PtTPC);
262             fThetaResolS1PtTPC->Fill(s1mcpt,deltaTheta1PtTPC);
263
264             fPtResolPtTPC->Fill(mcpt,deltaPtTPC);
265             fPtPullPtTPC->Fill(mcpt,pullPtTPC);
266             fPhiResolEtaTPC->Fill(eta,deltaPhiTPC);
267             fPhiPullEtaTPC->Fill(eta,pullPhiTPC);
268             fLambdaResolEtaTPC->Fill(eta,deltaLambdaTPC);
269             fLambdaPullEtaTPC->Fill(eta,pullLambdaTPC);
270           */
271         }
272     delete track;
273     }
274   }
275 }
276
277 //_____________________________________________________________________________
278 void AliComparisonRes::ProcessTPCITS(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
279 {
280   Int_t clusterITS[200];
281   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
282
283   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
284   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
285
286   Float_t mceta =  infoMC->GetParticle().Eta();
287   Float_t mcphi =  infoMC->GetParticle().Phi();
288   if(mcphi<0) mcphi += 2.*TMath::Pi();
289   Float_t mcpt = infoMC->GetParticle().Pt();
290
291   // distance to Prim. vertex 
292   const Double_t* dv = infoMC->GetVDist(); 
293   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
294
295   // Only 5 charged particle species (e,mu,pi,K,p)
296   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
297   if (!isPrim) return;
298   if (infoRC->GetStatus(1)!=3) return; // TPC refit
299   if (!infoRC->GetESDtrack()) return;  
300   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
301
302   infoRC->GetESDtrack()->GetImpactParameters(dca,cov);
303
304   // select primaries
305   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) 
306   { 
307     deltaYTPC= infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy();
308     deltaZTPC = infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz();
309     deltaLambdaTPC = TMath::ATan2(infoRC->GetESDtrack()->Pz(),infoRC->GetESDtrack()->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
310     deltaPhiTPC = TMath::ATan2(infoRC->GetESDtrack()->Py(),infoRC->GetESDtrack()->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
311     delta1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt)*mcpt;
312
313     pullYTPC= (infoRC->GetESDtrack()->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaY2());
314     pullZTPC = (infoRC->GetESDtrack()->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaZ2());
315     pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaTgl2());
316     pullPhiTPC = deltaPhiTPC / TMath::Sqrt(infoRC->GetESDtrack()->GetSigmaSnp2()); 
317     pull1PtTPC = (infoRC->GetESDtrack()->OneOverPt()-1./mcpt) / TMath::Sqrt(infoRC->GetESDtrack()->GetSigma1Pt2());
318
319     Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
320     Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
321
322     // TPC and ITS clusters in the system
323     if(infoRC->GetESDtrack()->GetITSclusters(clusterITS)>fCutsRC->GetMinNClustersITS()) 
324     {
325       fResolHisto->Fill(vResolHisto);
326       fPullHisto->Fill(vPullHisto);
327     }
328   }
329 }
330
331 //_____________________________________________________________________________
332 void AliComparisonRes::ProcessConstrained(AliMCInfo* const infoMC, AliESDRecInfo *const infoRC)
333 {
334   // Fill resolution comparison information (constarained parameters) 
335   //
336   AliExternalTrackParam *track = 0;
337   //Double_t field      = AliTracker::GetBz(); // nominal Bz field [kG]
338   Double_t kMaxD      = 123456.0; // max distance
339   AliESDVertex vertexMC;  // MC primary vertex
340
341   Double_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
342
343   Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC; 
344   Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC; 
345
346   Float_t mceta =  infoMC->GetParticle().Eta();
347   Float_t mcphi =  infoMC->GetParticle().Phi();
348   if(mcphi<0) mcphi += 2.*TMath::Pi();
349   Float_t mcpt = infoMC->GetParticle().Pt();
350
351   // distance to Prim. vertex 
352   const Double_t* dv = infoMC->GetVDist(); 
353   Bool_t isPrim = TMath::Sqrt(dv[0]*dv[0] + dv[1]*dv[1])<fCutsMC->GetMaxR() && TMath::Abs(dv[2])<fCutsMC->GetMaxVz();
354
355   // calculate and set prim. vertex
356   vertexMC.SetXv( infoMC->GetParticle().Vx() - dv[0] );
357   vertexMC.SetYv( infoMC->GetParticle().Vy() - dv[1] );
358   vertexMC.SetZv( infoMC->GetParticle().Vz() - dv[2] );
359
360   // Check selection cuts
361   if (fCutsMC->IsPdgParticle(TMath::Abs(infoMC->GetParticle().GetPdgCode())) == kFALSE) return; 
362   if (!isPrim) return;
363   if (infoRC->GetStatus(1)!=3) return; // TPC refit
364   if (!infoRC->GetESDtrack()) return;  
365   if (infoRC->GetESDtrack()->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
366
367   // constrained parameters resolution
368   const AliExternalTrackParam * cparam = infoRC->GetESDtrack()->GetConstrainedParam();
369   if(!cparam) return;
370
371   if ((track = new AliExternalTrackParam(*cparam)) != 0)
372   {
373     //Bool_t bDCAStatus = track->PropagateToDCA(&vertexMC,field,kMaxD,dca,cov);
374     Double_t field[3];  track->GetBxByBz(field); 
375     Bool_t bDCAStatus = track->PropagateToDCABxByBz(&vertexMC,field,kMaxD,dca,cov);
376     if(bDCAStatus) {
377       if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
378       {
379         deltaYTPC= track->GetY()-infoMC->GetParticle().Vy();
380         deltaZTPC = track->GetZ()-infoMC->GetParticle().Vz();
381         deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(infoMC->GetParticle().Pz(),infoMC->GetParticle().Pt());
382         deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(infoMC->GetParticle().Py(),infoMC->GetParticle().Px());
383         delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
384
385         pullYTPC= (track->GetY()-infoMC->GetParticle().Vy()) / TMath::Sqrt(track->GetSigmaY2());
386         pullZTPC = (track->GetZ()-infoMC->GetParticle().Vz()) / TMath::Sqrt(track->GetSigmaZ2());
387         pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
388         pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2()); 
389         pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
390
391         Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
392         fResolHisto->Fill(vResolHisto);
393
394         Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,infoMC->GetParticle().Vy(),infoMC->GetParticle().Vz(),mceta,mcphi,mcpt};
395         fPullHisto->Fill(vPullHisto);
396       }
397     }
398   delete track;
399   }
400 }
401  
402 //_____________________________________________________________________________
403 void AliComparisonRes::Exec(AliMCInfo* const infoMC, AliESDRecInfo* const infoRC){
404   
405   // Process comparison information 
406   if(GetAnalysisMode() == 0) ProcessTPC(infoMC,infoRC);
407   else if(GetAnalysisMode() == 1) ProcessTPCITS(infoMC,infoRC);
408   else if(GetAnalysisMode() == 2) ProcessConstrained(infoMC,infoRC);
409   else {
410     printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
411     return;
412   }
413 }
414
415 //_____________________________________________________________________________
416 TH1F* AliComparisonRes::MakeResol(TH2F * his, Int_t integ, Bool_t type){
417   // Create resolution histograms
418  
419   TH1F *hisr, *hism;
420   if (!gPad) new TCanvas;
421   hisr = AliTreeDraw::CreateResHistoI(his,&hism,integ);
422   if (type) return hism;
423   else 
424     return hisr;
425 }
426
427 //_____________________________________________________________________________
428 void AliComparisonRes::Analyse(){
429   // Analyse comparison information and store output histograms
430   // in the folder "folderRes"
431   //
432   TH1::AddDirectory(kFALSE);
433   TH1F *h=0;
434   TH2F *h2D=0;
435   TObjArray *aFolderObj = new TObjArray;
436
437   // write results in the folder 
438   TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
439   c->cd();
440
441   char name[256];
442   char res_title[256] = {"res_y:res_z:res_lambda:res_phi:res_1pt:y:z:eta:phi:pt"} ;
443   char pull_title[256] = {"pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt"};
444
445   for(Int_t i=0; i<5; i++) 
446   {
447     for(Int_t j=5; j<10; j++) 
448     {
449       if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.2,10.); // pt threshold
450       if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
451
452       h2D = (TH2F*)fResolHisto->Projection(i,j);
453       h = AliComparisonRes::MakeResol(h2D,1,0);
454       sprintf(name,"h_res_%d_vs_%d",i,j);
455       h->SetName(name);
456       h->SetTitle(res_title);
457
458       aFolderObj->Add(h);
459
460       h = AliComparisonRes::MakeResol(h2D,1,1);
461       sprintf(name,"h_mean_res_%d_vs_%d",i,j);
462       h->SetName(name);
463       h->SetTitle(res_title);
464
465       aFolderObj->Add(h);
466
467       if(j==7) fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
468       if(j==9) fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
469
470       //
471       if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.2,10.);
472       if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9);
473
474       h2D = (TH2F*)fPullHisto->Projection(i,j);
475       h = AliComparisonRes::MakeResol(h2D,1,0);
476       sprintf(name,"h_pull_%d_vs_%d",i,j);
477       h->SetName(name);
478       h->SetTitle(pull_title);
479
480       aFolderObj->Add(h);
481
482       h = AliComparisonRes::MakeResol(h2D,1,1);
483       sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
484       h->SetName(name);
485       h->SetTitle(pull_title);
486
487       aFolderObj->Add(h);
488
489       if(j==7) fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);
490       if(j==9) fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
491     }
492   }
493
494   // export objects to analysis folder
495   fAnalysisFolder = ExportToFolder(aFolderObj);
496
497   // delete only TObjArray
498   if(aFolderObj) delete aFolderObj;
499 }
500
501 //_____________________________________________________________________________
502 TFolder* AliComparisonRes::ExportToFolder(TObjArray * array) 
503 {
504   // recreate folder avery time and export objects to new one
505   //
506   AliComparisonRes * comp=this;
507   TFolder *folder = comp->GetAnalysisFolder();
508
509   TString name, title;
510   TFolder *newFolder = 0;
511   Int_t i = 0;
512   Int_t size = array->GetSize();
513
514   if(folder) { 
515      // get name and title from old folder
516      name = folder->GetName();  
517      title = folder->GetTitle();  
518
519          // delete old one
520      delete folder;
521
522          // create new one
523      newFolder = CreateFolder(name.Data(),title.Data());
524      newFolder->SetOwner();
525
526          // add objects to folder
527      while(i < size) {
528            newFolder->Add(array->At(i));
529            i++;
530          }
531   }
532
533 return newFolder;
534 }
535
536 //_____________________________________________________________________________
537 Long64_t AliComparisonRes::Merge(TCollection* const list) 
538 {
539   // Merge list of objects (needed by PROOF)
540
541   if (!list)
542   return 0;
543
544   if (list->IsEmpty())
545   return 1;
546
547   TIterator* iter = list->MakeIterator();
548   TObject* obj = 0;
549
550   // collection of generated histograms
551   Int_t count=0;
552   while((obj = iter->Next()) != 0) 
553   {
554   AliComparisonRes* entry = dynamic_cast<AliComparisonRes*>(obj);
555   if (entry == 0) continue; 
556
557   fResolHisto->Add(entry->fResolHisto);
558   fPullHisto->Add(entry->fPullHisto);
559
560   count++;
561   }
562
563 return count;
564 }
565
566 //_____________________________________________________________________________
567 TFolder* AliComparisonRes::CreateFolder(TString name,TString title) { 
568 // create folder for analysed histograms
569 //
570 TFolder *folder = 0;
571   folder = new TFolder(name.Data(),title.Data());
572
573   return folder;
574 }