]> 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
adding z vertex axis to AliTHn output (switch for using it requires coarse binning...
[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)
7768d177 62Double_t AliPerformanceRes::fgkMergeEntriesCut=5000000.; //5*10**6 tracks (small default to keep default memory foorprint low)
7cc34f08 63
64//_____________________________________________________________________________
65AliPerformanceRes::AliPerformanceRes():
66 AliPerformanceObject("AliPerformanceRes"),
67 fResolHisto(0),
68 fPullHisto(0),
69
70 // Cuts
71 fCutsRC(0),
72 fCutsMC(0),
73
74 // histogram folder
75 fAnalysisFolder(0)
76{
77 Init();
78}
79
80//_____________________________________________________________________________
81AliPerformanceRes::AliPerformanceRes(Char_t* name="AliPerformanceRes", Char_t* title="AliPerformanceRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
82 AliPerformanceObject(name,title),
83 fResolHisto(0),
84 fPullHisto(0),
85
86 // Cuts
87 fCutsRC(0),
88 fCutsMC(0),
89
90 // histogram folder
91 fAnalysisFolder(0)
92{
93 // named constructor
94 //
95 SetAnalysisMode(analysisMode);
96 SetHptGenerator(hptGenerator);
97
98 Init();
99}
100
101//_____________________________________________________________________________
102AliPerformanceRes::~AliPerformanceRes()
103{
104 // destructor
105
106 if(fResolHisto) delete fResolHisto; fResolHisto=0;
107 if(fPullHisto) delete fPullHisto; fPullHisto=0;
108
109 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
110}
111
112//_____________________________________________________________________________
113void AliPerformanceRes::Init(){
114
115 //
116 // histogram bining
117 //
118
119 // set pt bins
120 Int_t nPtBins = 50;
ef1ed531 121 Double_t ptMin = 1.e-2, ptMax = 20.;
7cc34f08 122
123 Double_t *binsPt = 0;
ef1ed531 124
7cc34f08 125 if (IsHptGenerator()) {
ef1ed531 126 ptMax = 100.;
127 }
128 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
7cc34f08 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
7cc34f08 229 // Fill TPC only resolution comparison information
e95d4942 230 const AliExternalTrackParam* tmpTrack = esdTrack->GetTPCInnerParam();
231 if(!tmpTrack) return;
232
233 AliExternalTrackParam track = *tmpTrack;
7cc34f08 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();
e95d4942 262 Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
7cc34f08 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;
e95d4942 272 double Bz = esdEvent->GetMagneticField();
273
274 Double_t mclocal[4]; //Rotated x,y,px,py mc-coordinates - the MC data should be rotated since the track is propagated best along x
275 Double_t c = TMath::Cos(track.GetAlpha());
276 Double_t s = TMath::Sin(track.GetAlpha());
277 Double_t x = particle->Vx();
278 Double_t y = particle->Vy();
279 mclocal[0] = x*c + y*s;
280 mclocal[1] =-x*s + y*c;
281 Double_t px = particle->Px();
282 Double_t py = particle->Py();
283 mclocal[2] = px*c + py*s;
284 mclocal[3] =-px*s + py*c;
285 Float_t mcsnplocal = TMath::Sin(TMath::ATan2(mclocal[3],mclocal[2]));
286
287
288 track.AliExternalTrackParam::PropagateTo(mclocal[0],Bz);
289
290 deltaYTPC= track.GetY()-mclocal[1];
291 deltaZTPC = track.GetZ()-particle->Vz();
292 deltaLambdaTPC = TMath::ATan2(track.Pz(),track.Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
293 //See comments in ProcessInnerTPC for remarks on local and global momentum coordinates for deltaPhi / pullSnp calculation
294 deltaPhiTPC = TMath::ATan2(track.Py(),track.Px())-TMath::ATan2(particle->Py(),particle->Px());
295 //delta1PtTPC = (track.OneOverPt()-1./mcpt)*mcpt;
296 deltaPtTPC = (track.Pt()-mcpt) / mcpt;
297
298 pullYTPC= deltaYTPC / TMath::Sqrt(track.GetSigmaY2());
299 pullZTPC = deltaZTPC / TMath::Sqrt(track.GetSigmaZ2());
7cc34f08 300
e95d4942 301 //Double_t sigma_lambda = 1./(1.+track.GetTgl()*track.GetTgl()) * TMath::Sqrt(track.GetSigmaTgl2());
302 //Double_t sigma_phi = 1./TMath::Sqrt(1-track.GetSnp()*track.GetSnp()) * TMath::Sqrt(track.GetSigmaSnp2());
303 pullPhiTPC = (track.GetSnp() - mcsnplocal) / TMath::Sqrt(track.GetSigmaSnp2());
304 pullLambdaTPC = (track.GetTgl() - mctgl) / TMath::Sqrt(track.GetSigmaTgl2());
7cc34f08 305
e95d4942 306 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track.GetSigmaTgl2());
307 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track.GetSigmaSnp2());
308 if (mcpt) pull1PtTPC = (track.OneOverPt()-1./mcpt) / TMath::Sqrt(track.GetSigma1Pt2());
309 else pull1PtTPC = 0.;
7cc34f08 310
311 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
312 fResolHisto->Fill(vResolHisto);
313
314 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
315 fPullHisto->Fill(vPullHisto);
316 }
317}
318
319//_____________________________________________________________________________
758320f7 320void AliPerformanceRes::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
7cc34f08 321{
322 // Fill resolution comparison information (TPC+ITS)
758320f7 323 if(!esdEvent) return;
7cc34f08 324 if(!esdTrack) return;
325
758320f7 326 if( IsUseTrackVertex() )
327 {
328 // Relate TPC inner params to prim. vertex
329 const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
330 Double_t x[3]; esdTrack->GetXYZ(x);
331 Double_t b[3]; AliTracker::GetBxByBz(x,b);
332 Bool_t isOK = esdTrack->RelateToVertexBxByBz(vtxESD, b, kVeryBig);
333 if(!isOK) return;
334
335 /*
336 // JMT -- recaluclate DCA for HLT if not present
337 if ( dca[0] == 0. && dca[1] == 0. ) {
338 track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
339 }
340 */
341 }
342
7cc34f08 343 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
344 esdTrack->GetImpactParameters(dca,cov);
345
346 //
347 // Fill rec vs MC information
348 //
349 if(!stack) return;
350
351 Int_t label = TMath::Abs(esdTrack->GetLabel());
352 TParticle* particle = stack->Particle(label);
353 if(!particle) return;
ff7b6e14 354 if(!particle->GetPDG()) return;
7cc34f08 355 if(particle->GetPDG()->Charge()==0) return;
356 //printf("charge %d \n",particle->GetPDG()->Charge());
357
358
359 // Only 5 charged particle species (e,mu,pi,K,p)
360 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
361
362 // exclude electrons
363 if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
364
365 Float_t mceta = particle->Eta();
366 Float_t mcphi = particle->Phi();
367 if(mcphi<0) mcphi += 2.*TMath::Pi();
368 Float_t mcpt = particle->Pt();
369 Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
370 Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
371
372 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
373 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
faa3b211 374 if(esdTrack->GetITSclusters(0)<fCutsRC->GetMinNClustersITS()) return; // min. nb. ITS clusters
7cc34f08 375
376 Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
377 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
378
379 // select primaries
380 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
381 {
382 if(mcpt == 0) return;
383
384 deltaYTPC= esdTrack->GetY()-particle->Vy();
385 deltaZTPC = esdTrack->GetZ()-particle->Vz();
386 deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
387 deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
388 //delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
389 deltaPtTPC = (esdTrack->Pt()-mcpt) / mcpt;
390
391 pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
392 pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
393
394 //Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
395 //Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
396 pullPhiTPC = (esdTrack->GetSnp() - mcsnp) / TMath::Sqrt(esdTrack->GetSigmaSnp2());
397 pullLambdaTPC = (esdTrack->GetTgl() - mctgl) / TMath::Sqrt(esdTrack->GetSigmaTgl2());
398
399 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
400 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
401 if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
402 else pull1PtTPC = 0.;
403
404 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
405 fResolHisto->Fill(vResolHisto);
406
407 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
408 fPullHisto->Fill(vPullHisto);
409
410
411 /*
412 deltaYTPC= esdTrack->GetY()-particle->Vy();
413 deltaZTPC = esdTrack->GetZ()-particle->Vz();
414 deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
415 deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
416 delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
417
418 pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
419 pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
420
421 Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
422 Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
423 pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
424 pullPhiTPC = deltaPhiTPC / sigma_phi;
425
426 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
427 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
428 if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
429 else pull1PtTPC = 0.;
430
431 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
432 fResolHisto->Fill(vResolHisto);
433
434 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
435 fPullHisto->Fill(vPullHisto);
436 */
437 }
438}
439
440//_____________________________________________________________________________
758320f7 441void AliPerformanceRes::ProcessConstrained(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
7cc34f08 442{
443 // Fill resolution comparison information (constarained parameters)
444 //
758320f7 445 if(!esdEvent) return;
7cc34f08 446 if(!esdTrack) return;
447
758320f7 448 if( IsUseTrackVertex() )
449 {
450 // Relate TPC inner params to prim. vertex
451 const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
452 Double_t x[3]; esdTrack->GetXYZ(x);
453 Double_t b[3]; AliTracker::GetBxByBz(x,b);
454 Bool_t isOK = esdTrack->RelateToVertexBxByBz(vtxESD, b, kVeryBig);
455 if(!isOK) return;
456
457 /*
458 // JMT -- recaluclate DCA for HLT if not present
459 if ( dca[0] == 0. && dca[1] == 0. ) {
460 track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
461 }
462 */
463 }
464
465
7cc34f08 466 const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
467 if(!track) return;
468
469 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
470 esdTrack->GetImpactParameters(dca,cov);
471
472 //
473 // Fill rec vs MC information
474 //
475 if(!stack) return;
476
477 Int_t label = TMath::Abs(esdTrack->GetLabel());
478 TParticle* particle = stack->Particle(label);
479 if(!particle) return;
ff7b6e14 480 if(!particle->GetPDG()) return;
7cc34f08 481 if(particle->GetPDG()->Charge()==0) return;
482 //printf("charge %d \n",particle->GetPDG()->Charge());
483
484 // Only 5 charged particle species (e,mu,pi,K,p)
485 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
486
487 // exclude electrons
488 if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
489
490 Float_t mceta = particle->Eta();
491 Float_t mcphi = particle->Phi();
492 if(mcphi<0) mcphi += 2.*TMath::Pi();
493 Float_t mcpt = particle->Pt();
494 Float_t mcsnp = TMath::Sin(TMath::ATan2(particle->Py(),particle->Px()));
495 Float_t mctgl = TMath::Tan(TMath::ATan2(particle->Pz(),particle->Pt()));
496
497 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
498 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
499
500 Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
501 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
502
503 // select primaries
504 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
505 {
506
507 if(mcpt == 0) return;
508
509 deltaYTPC= track->GetY()-particle->Vy();
510 deltaZTPC = track->GetZ()-particle->Vz();
511 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
512 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
513 //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
514 deltaPtTPC = (track->Pt()-mcpt) / mcpt;
515
516 pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
517 pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
518
519 //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
520 //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
521 pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
522 pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
523
524 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
525 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
526 if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
527 else pull1PtTPC = 0.;
528
529 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,particle->Vy(),particle->Vz(),mcphi,mceta,mcpt};
530 fResolHisto->Fill(vResolHisto);
531
532 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mcsnp,mctgl,1./mcpt};
533 fPullHisto->Fill(vPullHisto);
534
535 /*
536
537 deltaYTPC= track->GetY()-particle->Vy();
538 deltaZTPC = track->GetZ()-particle->Vz();
539 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
540 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
541 delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
542
543 pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
544 pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
545
546 Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
547 Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
548 pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
549 pullPhiTPC = deltaPhiTPC / sigma_phi;
550
551 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
552 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
553
554 if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
555 else pull1PtTPC = 0.;
556
557 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
558 fResolHisto->Fill(vResolHisto);
559
560 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
561 fPullHisto->Fill(vPullHisto);
562
563 */
564 }
565}
566
567//_____________________________________________________________________________
758320f7 568void AliPerformanceRes::ProcessInnerTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
7cc34f08 569{
570 //
571 // Fill resolution comparison information (inner params at TPC reference point)
572 //
758320f7 573 if(!esdEvent) return;
7cc34f08 574 if(!esdTrack) return;
575
758320f7 576 if( IsUseTrackVertex() )
577 {
578 // Relate TPC inner params to prim. vertex
579 const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
580 Double_t x[3]; esdTrack->GetXYZ(x);
581 Double_t b[3]; AliTracker::GetBxByBz(x,b);
582 Bool_t isOK = esdTrack->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
583 if(!isOK) return;
584
585 /*
586 // JMT -- recaluclate DCA for HLT if not present
587 if ( dca[0] == 0. && dca[1] == 0. ) {
588 track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
589 }
590 */
591 }
592
7cc34f08 593 const AliExternalTrackParam * innerParam = esdTrack->GetInnerParam();
594 if(!innerParam) return;
595
596 // create new AliExternalTrackParam
597 AliExternalTrackParam *track = new AliExternalTrackParam(*innerParam);
598 if(!track) return;
599
600 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
601 esdTrack->GetImpactParametersTPC(dca,cov);
602
603 //
604 // Fill rec vs MC information
605 //
606 if(!mcEvent) return;
607
608 Int_t label = TMath::Abs(esdTrack->GetLabel());
609 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(label);
610 if(!mcParticle) return;
611
612 // get the first TPC track reference
613 AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle);
614 if(!ref0) return;
615
616 // Only 5 charged particle species (e,mu,pi,K,p)
617 TParticle *particle = mcParticle->Particle();
618 if(!particle) return;
619
620 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
621
622 // exclude electrons
623 if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
624
e95d4942 625 Double_t mclocal[4]; //Rotated x,y,Px,Py mc-coordinates - the MC data should be rotated since the track is propagated best along x
626 Double_t c = TMath::Cos(track->GetAlpha());
627 Double_t s = TMath::Sin(track->GetAlpha());
628 Double_t x = ref0->X();
629 Double_t y = ref0->Y();
630 mclocal[0] = x*c + y*s;
631 mclocal[1] =-x*s + y*c;
632 Double_t px = ref0->Px();
633 Double_t py = ref0->Py();
634 mclocal[2] = px*c + py*s;
635 mclocal[3] =-px*s + py*c;
7cc34f08 636
e95d4942 637 Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
638 // propagate track to the radius of the first track reference within TPC
7cc34f08 639 Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
7cc34f08 640 Double_t field[3]; track->GetBxByBz(field);
e95d4942 641 if (TGeoGlobalMagField::Instance()->GetField() == NULL) {
642 Error("ProcessInnerTPC", "Magnetic Field not set");
643 }
644 Bool_t isOK = track->PropagateToBxByBz(mclocal[0],field);
645 if(!isOK) {return;}
646 //track->AliExternalTrackParam::PropagateTo(mclocal[0],esdEvent->GetMagneticField()); //Use different propagation since above methods returns zero B field and does not work
647
7cc34f08 648 Float_t mceta = -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
649 Float_t mcphi = ref0->Phi();
650 if(mcphi<0) mcphi += 2.*TMath::Pi();
651 Float_t mcpt = ref0->Pt();
652 Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
e95d4942 653 Float_t mcsnplocal = TMath::Sin(TMath::ATan2(mclocal[3],mclocal[2]));
7cc34f08 654 Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
655
656 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
657 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
658
659 Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
660 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
661
662 // select primaries
e95d4942 663 Bool_t isPrimary;
664 if ( IsUseTrackVertex() )
665 {
666 isPrimary = TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ();
667 }
668 else
669 {
670 //If Track vertex is not used the above check does not work, hence we use the MC reference track
671 isPrimary = label < mcEvent->Stack()->GetNprimary();
672 }
673 if(isPrimary)
7cc34f08 674 {
675 if(mcpt == 0) return;
676
e95d4942 677 deltaYTPC= track->GetY()-mclocal[1];
7cc34f08 678 deltaZTPC = track->GetZ()-ref0->Z();
679 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
e95d4942 680 //track->Px() etc returns momentum in global coordinates, hence mc momentum must not be rotated.
7cc34f08 681 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
682 //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
683 deltaPtTPC = (track->Pt()-mcpt) / mcpt;
684
e95d4942 685 pullYTPC= deltaYTPC / TMath::Sqrt(track->GetSigmaY2());
686 pullZTPC = deltaZTPC / TMath::Sqrt(track->GetSigmaZ2());
7cc34f08 687
688 //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
689 //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
e95d4942 690 //track->GetSnp returns value in local coordinates, hence here, in contrast to deltaPhiTPC, the mc momentum must be rotated
691 pullPhiTPC = (track->GetSnp() - mcsnplocal) / TMath::Sqrt(track->GetSigmaSnp2());
7cc34f08 692 pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
693
694 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
695 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
696 if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
697 else pull1PtTPC = 0.;
698
699 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
700 fResolHisto->Fill(vResolHisto);
701
702 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
703 fPullHisto->Fill(vPullHisto);
704 }
705
706 if(track) delete track;
707}
708
709//_____________________________________________________________________________
758320f7 710void AliPerformanceRes::ProcessOuterTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack, AliESDfriendTrack *const friendTrack, AliESDEvent* const esdEvent)
7cc34f08 711{
712 //
713 // Fill resolution comparison information (outer params at TPC reference point)
714 //
715 if(!friendTrack) return;
758320f7 716 if(!esdEvent) return;
717 if(!esdTrack) return;
718
719 if( IsUseTrackVertex() )
720 {
721 // Relate TPC inner params to prim. vertex
722 const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
723 Double_t x[3]; esdTrack->GetXYZ(x);
724 Double_t b[3]; AliTracker::GetBxByBz(x,b);
725 Bool_t isOK = esdTrack->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
726 if(!isOK) return;
727
728 /*
729 // JMT -- recaluclate DCA for HLT if not present
730 if ( dca[0] == 0. && dca[1] == 0. ) {
731 track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
732 }
733 */
734 }
7cc34f08 735
736 const AliExternalTrackParam * outerParam = friendTrack->GetTPCOut();
737 if(!outerParam) return;
738
739 // create new AliExternalTrackParam
740 AliExternalTrackParam *track = new AliExternalTrackParam(*outerParam);
741 if(!track) return;
742
743 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
744 esdTrack->GetImpactParametersTPC(dca,cov);
745
746 //
747 // Fill rec vs MC information
748 //
749 if(!mcEvent) return;
750
751 Int_t label = TMath::Abs(esdTrack->GetLabel());
752 AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(label);
753 if(!mcParticle) return;
754
755 // get the last TPC track reference
756 AliTrackReference *ref0 = GetLastTPCTrackRef(mcParticle);
757 if(!ref0) return;
758
759 // Only 5 charged particle species (e,mu,pi,K,p)
760 TParticle *particle = mcParticle->Particle();
761 if(!particle) return;
762
763 if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) return;
764
765 // exclude electrons
766 if (fCutsMC->GetEM()==TMath::Abs(particle->GetPdgCode())) return;
767
768 // calculate alpha angle
769 Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
770 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
771 //if (alpha<0) alpha += TMath::TwoPi();
772
773 // rotate outer track to local coordinate system
774 // and propagate to the radius of the last track reference of TPC
775 Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
776 //Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
777 Double_t field[3]; track->GetBxByBz(field);
778 Bool_t isOK = track->PropagateBxByBz(alpha,trRadius,field);
779 if(!isOK) return;
780
781 Float_t mceta = -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
782 Float_t mcphi = ref0->Phi();
783 if(mcphi<0) mcphi += 2.*TMath::Pi();
784 Float_t mcpt = ref0->Pt();
785 Float_t mcsnp = TMath::Sin(TMath::ATan2(ref0->Py(),ref0->Px()));
786 Float_t mctgl = TMath::Tan(TMath::ATan2(ref0->Pz(),ref0->Pt()));
787
788 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
789 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
790
791 Float_t deltaPtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
792 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
793
794 // select primaries
795 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
796 {
797 if(mcpt == 0) return;
798
799 deltaYTPC= track->GetY(); // already rotated
800 deltaZTPC = track->GetZ()-ref0->Z();
801 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
802 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
803 //delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
804 deltaPtTPC = (track->Pt()-mcpt) / mcpt;
805
806 pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2());
807 pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
808
809 //Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
810 //Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
811 pullPhiTPC = (track->GetSnp() - mcsnp) / TMath::Sqrt(track->GetSigmaSnp2());
812 pullLambdaTPC = (track->GetTgl() - mctgl) / TMath::Sqrt(track->GetSigmaTgl2());
813
814 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
815 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
816 if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
817 else pull1PtTPC = 0.;
818
819 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,deltaPtTPC,ref0->Y(),ref0->Z(),mcphi,mceta,mcpt};
820 fResolHisto->Fill(vResolHisto);
821
822 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mcsnp,mctgl,1./mcpt};
823 fPullHisto->Fill(vPullHisto);
824 }
825
826 if(track) delete track;
827}
828
829//_____________________________________________________________________________
830AliTrackReference * AliPerformanceRes::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
831{
832 // return first TPC track reference
833 if(!mcParticle) return 0;
834
835 // find first track reference
836 // check direction to select proper reference point for looping tracks
837 Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
838 AliTrackReference *ref = 0;
839 AliTrackReference *refIn = 0;
840 for (Int_t iref = 0; iref < nTrackRef; iref++) {
841 ref = mcParticle->GetTrackReference(iref);
842 if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
843 {
844 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
845 if(dir < 0.) break;
846
847 refIn = ref;
848 break;
849 }
850 }
851
852return refIn;
853}
854
855//_____________________________________________________________________________
856AliTrackReference * AliPerformanceRes::GetLastTPCTrackRef(AliMCParticle *mcParticle)
857{
858 // return last TPC track reference
859 if(!mcParticle) return 0;
860
861 // find last track reference
862 // check direction to select proper reference point for looping tracks
863 Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
864 AliTrackReference *ref = 0;
865 AliTrackReference *refOut = 0;
866 for (Int_t iref = 0; iref < nTrackRef; iref++) {
867 ref = mcParticle->GetTrackReference(iref);
868 if(ref && (ref->DetectorId()==AliTrackReference::kTPC)) {
869 Float_t dir=ref->X()*ref->Px()+ref->Y()*ref->Py();
870 if(dir< 0.0) break;
871 refOut = ref;
872 }
873 }
874
875return refOut;
876}
877
878//_____________________________________________________________________________
879void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
880{
881 // Process comparison information
7cc34f08 882 if(!esdEvent)
883 {
884 Error("Exec","esdEvent not available");
885 return;
886 }
887 AliHeader* header = 0;
888 AliGenEventHeader* genHeader = 0;
889 AliStack* stack = 0;
890 TArrayF vtxMC(3);
891
892 if(bUseMC)
893 {
894 if(!mcEvent) {
895 Error("Exec","mcEvent not available");
896 return;
897 }
898 // get MC event header
899 header = mcEvent->Header();
900 if (!header) {
901 Error("Exec","Header not available");
902 return;
903 }
904 // MC particle stack
905 stack = mcEvent->Stack();
906 if (!stack) {
907 Error("Exec","Stack not available");
908 return;
909 }
910 // get MC vertex
911 genHeader = header->GenEventHeader();
912 if (!genHeader) {
913 Error("Exec","Could not retrieve genHeader from Header");
914 return;
915 }
916 genHeader->PrimaryVertex(vtxMC);
917 }
918 else {
919 Error("Exec","MC information required!");
920 return;
921 }
922
923 // use ESD friends
924 if(bUseESDfriend) {
925 if(!esdFriend) {
926 Error("Exec","esdFriend not available");
927 return;
928 }
929 }
930
758320f7 931 // get event vertex
932 const AliESDVertex *vtxESD = NULL;
933 if( IsUseTrackVertex() )
934 {
935 // track vertex
936 vtxESD = esdEvent->GetPrimaryVertexTracks();
e95d4942 937 if(vtxESD && (vtxESD->GetStatus()<=0)) return;
758320f7 938 }
5a88d599 939 // Coverity - removed else branch as vtxESD is not further used in method
940 // else {
941 // // TPC track vertex
942 // vtxESD = esdEvent->GetPrimaryVertexTPC();
943 // }
e95d4942 944
758320f7 945
946
947
7cc34f08 948 // Process events
949 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
950 {
951 AliESDtrack *track = esdEvent->GetTrack(iTrack);
952 if(!track) continue;
953
954 AliESDfriendTrack *friendTrack=0;
583e3b45 955
7cc34f08 956
758320f7 957 Int_t label = TMath::Abs(track->GetLabel());
958 if ( label > stack->GetNtrack() )
959 {
960 ULong_t status = track->GetStatus();
961 printf ("Error : ESD MCLabel %d - StackSize %d - Status %lu \n",
962 track->GetLabel(), stack->GetNtrack(), status );
963 printf(" NCluster %d \n", track->GetTPCclusters(0) );
964 /*
965 if ((status&AliESDtrack::kTPCrefit)== 0 ) printf(" kTPCrefit \n");
966 if ((status&AliESDtrack::kTPCin)== 0 ) printf(" kTPCin \n");
967 if ((status&AliESDtrack::kTPCout)== 0 ) printf(" kTPCout \n");
968 if ((status&AliESDtrack::kTRDrefit)== 0 ) printf(" kTRDrefit \n");
969 if ((status&AliESDtrack::kTRDin)== 0 ) printf(" kTRDin \n");
970 if ((status&AliESDtrack::kTRDout)== 0 ) printf(" kTRDout \n");
971 if ((status&AliESDtrack::kITSrefit)== 0 ) printf(" kITSrefit \n");
972 if ((status&AliESDtrack::kITSin)== 0 ) printf(" kITSin \n");
973 if ((status&AliESDtrack::kITSout)== 0 ) printf(" kITSout \n");
974 */
975
976 continue;
977 }
978
e95d4942 979 if (label == 0) continue; //Cannot distinguish between track or fake track
980 if (track->GetLabel() < 0) continue; //Do not consider fake tracks
981
758320f7 982 if(GetAnalysisMode() == 0) ProcessTPC(stack,track,esdEvent);
983 else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track,esdEvent);
984 else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track,esdEvent);
985 else if(GetAnalysisMode() == 3) ProcessInnerTPC(mcEvent,track,esdEvent);
583e3b45 986 else if(GetAnalysisMode() == 4) {
987
988 if(bUseESDfriend) {
989 friendTrack=esdFriend->GetTrack(iTrack);
990 if(!friendTrack) continue;
991 }
992
993 ProcessOuterTPC(mcEvent,track,friendTrack,esdEvent);
994 }
7cc34f08 995 else {
996 printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
997 return;
998 }
999 }
1000}
1001
1002//_____________________________________________________________________________
1003TH1F* AliPerformanceRes::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
1004 // Create resolution histograms
1005
1006 TH1F *hisr, *hism;
1007 if (!gPad) new TCanvas;
1008 hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
1009 if (type) return hism;
1010 else
1011 return hisr;
1012}
1013
1014//_____________________________________________________________________________
1015void AliPerformanceRes::Analyse() {
1016 // Analyse comparison information and store output histograms
1017 // in the folder "folderRes"
1018 //
1019 TH1::AddDirectory(kFALSE);
1020 TH1F *h=0;
1021 TH2F *h2D=0;
1022 TObjArray *aFolderObj = new TObjArray;
18e588e9 1023 if(!aFolderObj) return;
7cc34f08 1024
1025 // write results in the folder
1026 TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
1027 c->cd();
1028
1029 char name[256];
1030 char title[256];
1031
1032 for(Int_t i=0; i<5; i++)
1033 {
1034 for(Int_t j=5; j<10; j++)
1035 {
e95d4942 1036 if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
1037 //if(j!=8) fResolHisto->GetAxis(8)->SetRangeUser(0.0,0.89); // eta window
7cc34f08 1038 else fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.49);
1039 fResolHisto->GetAxis(9)->SetRangeUser(0.16,100.); // pt threshold
1040 if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
1041
1042 h2D = (TH2F*)fResolHisto->Projection(i,j);
1043
1044 h = AliPerformanceRes::MakeResol(h2D,1,0,100);
18e588e9 1045 snprintf(name,256,"h_res_%d_vs_%d",i,j);
7cc34f08 1046 h->SetName(name);
1047
1048 h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
18e588e9 1049 snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
7cc34f08 1050 h->GetYaxis()->SetTitle(title);
18e588e9 1051 snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
7cc34f08 1052 h->SetTitle(title);
1053
1054 if(j==9) h->SetBit(TH1::kLogX);
1055 aFolderObj->Add(h);
1056
1057 h = AliPerformanceRes::MakeResol(h2D,1,1,100);
1058 //h = (TH1F*)arr->At(1);
18e588e9 1059 snprintf(name,256,"h_mean_res_%d_vs_%d",i,j);
7cc34f08 1060 h->SetName(name);
1061
1062 h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
18e588e9 1063 snprintf(title,256,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
7cc34f08 1064 h->GetYaxis()->SetTitle(title);
1065
18e588e9 1066 snprintf(title,256,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
7cc34f08 1067 h->SetTitle(title);
1068
1069 if(j==9) h->SetBit(TH1::kLogX);
1070 aFolderObj->Add(h);
1071
1072 fResolHisto->GetAxis(8)->SetRangeUser(-1.5,1.5);
1073 fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
1074
1075 //
1076 //if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(-0.9,0.89); // eta window
1077 if(j!=8) fPullHisto->GetAxis(8)->SetRangeUser(0.0,0.89); // eta window
1078 else fPullHisto->GetAxis(8)->SetRangeUser(-1.5,1.49); // eta window
1079 fPullHisto->GetAxis(9)->SetRangeUser(0.16,100.); // pt threshold
1080 if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
1081
1082 h2D = (TH2F*)fPullHisto->Projection(i,j);
1083
1084 h = AliPerformanceRes::MakeResol(h2D,1,0,100);
18e588e9 1085 snprintf(name,256,"h_pull_%d_vs_%d",i,j);
7cc34f08 1086 h->SetName(name);
1087
1088 h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
18e588e9 1089 snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
7cc34f08 1090 h->GetYaxis()->SetTitle(title);
18e588e9 1091 snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
7cc34f08 1092 h->SetTitle(title);
1093
1094 //if(j==9) h->SetBit(TH1::kLogX);
1095 aFolderObj->Add(h);
1096
1097 h = AliPerformanceRes::MakeResol(h2D,1,1,100);
18e588e9 1098 snprintf(name,256,"h_mean_pull_%d_vs_%d",i,j);
7cc34f08 1099 h->SetName(name);
1100
1101 h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
18e588e9 1102 snprintf(title,256,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
7cc34f08 1103 h->GetYaxis()->SetTitle(title);
18e588e9 1104 snprintf(title,256,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
7cc34f08 1105 h->SetTitle(title);
1106
1107 //if(j==9) h->SetBit(TH1::kLogX);
1108 aFolderObj->Add(h);
1109 }
1110 }
1111
e95d4942 1112 for (Int_t i = 0;i < fResolHisto->GetNdimensions();i++)
1113 {
1114 fResolHisto->GetAxis(i)->SetRange(1,0); //Reset Range
1115 }
1116 for (Int_t i = 0;i < fPullHisto->GetNdimensions();i++)
1117 {
1118 fPullHisto->GetAxis(i)->SetRange(1,0); //Reset Range
1119 }
1120
7cc34f08 1121 // export objects to analysis folder
1122 fAnalysisFolder = ExportToFolder(aFolderObj);
1123
1124 // delete only TObjArray
1125 if(aFolderObj) delete aFolderObj;
1126}
1127
1128//_____________________________________________________________________________
1129TFolder* AliPerformanceRes::ExportToFolder(TObjArray * array)
1130{
1131 // recreate folder avery time and export objects to new one
1132 //
1133 AliPerformanceRes * comp=this;
1134 TFolder *folder = comp->GetAnalysisFolder();
1135
1136 TString name, title;
1137 TFolder *newFolder = 0;
1138 Int_t i = 0;
1139 Int_t size = array->GetSize();
1140
1141 if(folder) {
1142 // get name and title from old folder
1143 name = folder->GetName();
1144 title = folder->GetTitle();
1145
1146 // delete old one
1147 delete folder;
1148
1149 // create new one
1150 newFolder = CreateFolder(name.Data(),title.Data());
1151 newFolder->SetOwner();
1152
1153 // add objects to folder
1154 while(i < size) {
1155 newFolder->Add(array->At(i));
1156 i++;
1157 }
1158 }
1159
1160return newFolder;
1161}
1162
1163//_____________________________________________________________________________
1164Long64_t AliPerformanceRes::Merge(TCollection* const list)
1165{
1166 // Merge list of objects (needed by PROOF)
7729608b 1167
7cc34f08 1168 if (!list)
1169 return 0;
1170
1171 if (list->IsEmpty())
1172 return 1;
1173
1174 TIterator* iter = list->MakeIterator();
1175 TObject* obj = 0;
1176
1177 // collection of generated histograms
1178 Int_t count=0;
1179 while((obj = iter->Next()) != 0)
1180 {
1181 AliPerformanceRes* entry = dynamic_cast<AliPerformanceRes*>(obj);
1182 if (entry == 0) continue;
7729608b 1183 if (fResolHisto->GetEntries()<fgkMergeEntriesCut){
1184 fResolHisto->Add(entry->fResolHisto);
1185 fPullHisto->Add(entry->fPullHisto);
1186 }
7cc34f08 1187
1188 count++;
1189 }
1190
1191return count;
1192}
1193
1194//_____________________________________________________________________________
1195TFolder* AliPerformanceRes::CreateFolder(TString name,TString title) {
1196// create folder for analysed histograms
1197//
1198TFolder *folder = 0;
1199 folder = new TFolder(name.Data(),title.Data());
1200
1201 return folder;
1202}