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