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