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