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