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