Production point stored in AliFastGlauber.
[u/mrichter/AliRoot.git] / PWG1 / AliPerformanceRes.cxx
CommitLineData
777a0ba6 1//------------------------------------------------------------------------------
2// Implementation of AliPerformanceRes class. It keeps information from
3// comparison of reconstructed and MC particle tracks. In addtion,
4// it keeps selection cuts used during comparison. The comparison
5// information is stored in the ROOT histograms. Analysis of these
6// histograms can be done by using Analyse() class function. The result of
7// the analysis (histograms/graphs) are stored in the folder which is
8// a data member of AliPerformanceRes.
9//
10// Author: J.Otwinowski 04/02/2008
11//------------------------------------------------------------------------------
12
13/*
14
15 // after running comparison task, read the file, and get component
16 gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
17 LoadMyLibs();
18
19 TFile f("Output.root");
20 //AliPerformanceRes * compObj = (AliPerformanceRes*)f.Get("AliPerformanceRes");
21 AliPerformanceRes * compObj = (AliPerformanceRes*)cOutput->FindObject("AliPerformanceRes");
22
23 // analyse comparison data
24 compObj->Analyse();
25
26 // the output histograms/graphs will be stored in the folder "folderRes"
27 compObj->GetAnalysisFolder()->ls("*");
28
29 // user can save whole comparison object (or only folder with anlysed histograms)
30 // in the seperate output file (e.g.)
31 TFile fout("Analysed_Res.root","recreate");
32 compObj->Write(); // compObj->GetAnalysisFolder()->Write();
33 fout.Close();
34
35*/
36
37#include "TCanvas.h"
38#include "TH1.h"
39#include "TH2.h"
40#include "TAxis.h"
41#include "TF1.h"
42
43#include "AliPerformanceRes.h"
44#include "AliESDEvent.h"
45#include "AliESDVertex.h"
46#include "AliESDtrack.h"
47#include "AliLog.h"
48#include "AliMCEvent.h"
49#include "AliMCParticle.h"
50#include "AliHeader.h"
51#include "AliGenEventHeader.h"
52#include "AliStack.h"
53#include "AliMCInfoCuts.h"
54#include "AliRecInfoCuts.h"
55#include "AliTracker.h"
56#include "AliTreeDraw.h"
57
58using namespace std;
59
60ClassImp(AliPerformanceRes)
61
62//_____________________________________________________________________________
63AliPerformanceRes::AliPerformanceRes():
64 AliPerformanceObject("AliPerformanceRes"),
65 fResolHisto(0),
66 fPullHisto(0),
67
68 // Cuts
69 fCutsRC(0),
70 fCutsMC(0),
71
72 // histogram folder
73 fAnalysisFolder(0)
74{
75 Init();
76}
77
78//_____________________________________________________________________________
79AliPerformanceRes::AliPerformanceRes(Char_t* name="AliPerformanceRes", Char_t* title="AliPerformanceRes",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
80 AliPerformanceObject(name,title),
81 fResolHisto(0),
82 fPullHisto(0),
83
84 // Cuts
85 fCutsRC(0),
86 fCutsMC(0),
87
88 // histogram folder
89 fAnalysisFolder(0)
90{
91 // named constructor
92 //
93 SetAnalysisMode(analysisMode);
94 SetHptGenerator(hptGenerator);
95
96 Init();
97}
98
99//_____________________________________________________________________________
100AliPerformanceRes::~AliPerformanceRes()
101{
102 // destructor
103
104 if(fResolHisto) delete fResolHisto; fResolHisto=0;
105 if(fPullHisto) delete fPullHisto; fPullHisto=0;
106
107 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
108}
109
110//_____________________________________________________________________________
111void AliPerformanceRes::Init(){
112
113 //
114 // histogram bining
115 //
116
117 // set pt bins
118 Int_t nPtBins = 50;
119 Double_t ptMin = 1.e-2, ptMax = 10.;
120
121 Double_t *binsPt = 0;
122 if (IsHptGenerator()) {
123 nPtBins = 100; ptMax = 100.;
124 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
125 } else {
126 binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
127 }
128
129 /*
130 Int_t nPtBins = 31;
131 Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
132 Double_t ptMin = 0., ptMax = 10.;
133
134 if(IsHptGenerator() == kTRUE) {
135 nPtBins = 100;
136 ptMin = 0.; ptMax = 100.;
137 }
138 */
139
140 Double_t yMin = -0.02, yMax = 0.02;
141 Double_t zMin = -12.0, zMax = 12.0;
142 if(GetAnalysisMode() == 3) { // TrackRef coordinate system
143 yMin = -100.; yMax = 100.;
144 zMin = -100.; zMax = 100.;
145 }
146
147 // res_y:res_z:res_phi,res_lambda:res_1pt:y:z:eta:phi:pt
148 Int_t binsResolHisto[10]={100,100,100,100,100,25,50,30,144,nPtBins};
149 Double_t minResolHisto[10]={-1.,-1.,-0.03,-0.03,-0.2, yMin, zMin,-1.5,0.,ptMin};
150 Double_t maxResolHisto[10]={ 1., 1., 0.03, 0.03, 0.2, yMax, zMax, 1.5,2.*TMath::Pi(), ptMax};
151
152 fResolHisto = new THnSparseF("fResolHisto","res_y:res_z:res_phi:res_lambda:res_1pt:y:z:eta:phi:pt",10,binsResolHisto,minResolHisto,maxResolHisto);
153 fResolHisto->SetBinEdges(9,binsPt);
154
155 fResolHisto->GetAxis(0)->SetTitle("y-y_{mc} (cm)");
156 fResolHisto->GetAxis(1)->SetTitle("z-z_{mc} (cm)");
157 fResolHisto->GetAxis(2)->SetTitle("#phi-#phi_{mc} (rad)");
158 fResolHisto->GetAxis(3)->SetTitle("#lambda-#lambda_{mc} (rad)");
159 fResolHisto->GetAxis(4)->SetTitle("(p_{Tmc}/p_{T}-1)");
160 fResolHisto->GetAxis(5)->SetTitle("y_{mc} (cm)");
161 fResolHisto->GetAxis(6)->SetTitle("z_{mc} (cm)");
162 fResolHisto->GetAxis(7)->SetTitle("#eta_{mc}");
163 fResolHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
164 fResolHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
165 fResolHisto->Sumw2();
166
167 //pull_y:pull_z:pull_lambda:pull_phi:pull_1pt:y:z:eta:phi:pt
168 Int_t binsPullHisto[10]={100,100,100,100,100,50,50,30,144,nPtBins};
169 Double_t minPullHisto[10]={-5.,-5.,-5.,-5.,-5.,yMin, zMin,-1.5, 0., ptMin};
170 Double_t maxPullHisto[10]={ 5., 5., 5., 5., 5., yMax, zMax, 1.5, 2.*TMath::Pi(),ptMax};
171
172 fPullHisto = new THnSparseF("fPullHisto","pull_y:pull_z:pull_phi:pull_lambda:pull_1pt:y:z:eta:phi:pt",10,binsPullHisto,minPullHisto,maxPullHisto);
173 if(!IsHptGenerator()) fPullHisto->SetBinEdges(9,binsPt);
174
175 fPullHisto->GetAxis(0)->SetTitle("(y-y_{mc})/#sigma");
176 fPullHisto->GetAxis(1)->SetTitle("(z-z_{mc})/#sigma");
177 fPullHisto->GetAxis(2)->SetTitle("(#phi-#phi_{mc})/#sigma");
178 fPullHisto->GetAxis(3)->SetTitle("(#lambda-#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("#eta_{mc}");
183 fPullHisto->GetAxis(8)->SetTitle("#phi_{mc} (rad)");
184 fPullHisto->GetAxis(9)->SetTitle("p_{Tmc} (GeV/c)");
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//_____________________________________________________________________________
198void AliPerformanceRes::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
199{
200 if(!esdTrack) return;
201
202 // Fill TPC only resolution comparison information
203 const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
204 if(!track) return;
205
206 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
207 esdTrack->GetImpactParametersTPC(dca,cov);
208
209 //
210 // Fill rec vs MC information
211 //
212 if(!stack) return;
213 Int_t label = TMath::Abs(esdTrack->GetLabel());
214 if(label == 0) return;
215
216 TParticle* particle = stack->Particle(label);
217 if(!particle) return;
218 if(!particle->GetPDG()) return;
219 if(particle->GetPDG()->Charge()==0) return;
220 //printf("charge %d \n",particle->GetPDG()->Charge());
221
222 Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
223 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
224
225 Float_t mceta = particle->Eta();
226 Float_t mcphi = particle->Phi();
227 if(mcphi<0) mcphi += 2.*TMath::Pi();
228 Float_t mcpt = particle->Pt();
229
230 // nb. TPC clusters cut
231 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return;
232
233 // select primaries
234 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
235 {
236 deltaYTPC= track->GetY()-particle->Vy();
237 deltaZTPC = track->GetZ()-particle->Vz();
238 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
239 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
240 delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
241
242 pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
243 pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
244
245 Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
246 Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
247 pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
248 pullPhiTPC = deltaPhiTPC / sigma_phi;
249
250 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
251 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
252 if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
253 else pull1PtTPC = 0.;
254
255 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
256 fResolHisto->Fill(vResolHisto);
257
258 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
259 fPullHisto->Fill(vPullHisto);
260 }
261}
262
263//_____________________________________________________________________________
264void AliPerformanceRes::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack)
265{
266 // Fill resolution comparison information (TPC+ITS)
267 if(!esdTrack) return;
268
269 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
270 esdTrack->GetImpactParameters(dca,cov);
271
272 //
273 // Fill rec vs MC information
274 //
275 if(!stack) return;
276
277 Int_t label = TMath::Abs(esdTrack->GetLabel());
278 if(label == 0) return;
279
280 TParticle* particle = stack->Particle(label);
281 if(!particle) return;
282 if(particle->GetPDG()->Charge()==0) return;
283 //printf("charge %d \n",particle->GetPDG()->Charge());
284
285 Float_t mceta = particle->Eta();
286 Float_t mcphi = particle->Phi();
287 if(mcphi<0) mcphi += 2.*TMath::Pi();
288 Float_t mcpt = particle->Pt();
289
290 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
291 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
292 Int_t clusterITS[200];
293 if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return; // min. nb. ITS clusters
294
295 Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
296 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
297
298 // select primaries
299 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
300 {
301 deltaYTPC= esdTrack->GetY()-particle->Vy();
302 deltaZTPC = esdTrack->GetZ()-particle->Vz();
303 deltaLambdaTPC = TMath::ATan2(esdTrack->Pz(),esdTrack->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
304 deltaPhiTPC = TMath::ATan2(esdTrack->Py(),esdTrack->Px())-TMath::ATan2(particle->Py(),particle->Px());
305 delta1PtTPC = (esdTrack->OneOverPt()-1./mcpt)*mcpt;
306
307 pullYTPC= (esdTrack->GetY()-particle->Vy()) / TMath::Sqrt(esdTrack->GetSigmaY2());
308 pullZTPC = (esdTrack->GetZ()-particle->Vz()) / TMath::Sqrt(esdTrack->GetSigmaZ2());
309
310 Double_t sigma_lambda = 1./(1.+esdTrack->GetTgl()*esdTrack->GetTgl()) * TMath::Sqrt(esdTrack->GetSigmaTgl2());
311 Double_t sigma_phi = 1./TMath::Sqrt(1-esdTrack->GetSnp()*esdTrack->GetSnp()) * TMath::Sqrt(esdTrack->GetSigmaSnp2());
312 pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
313 pullPhiTPC = deltaPhiTPC / sigma_phi;
314
315 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(esdTrack->GetSigmaTgl2());
316 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(esdTrack->GetSigmaSnp2());
317 if (mcpt) pull1PtTPC = (esdTrack->OneOverPt()-1./mcpt) / TMath::Sqrt(esdTrack->GetSigma1Pt2());
318 else pull1PtTPC = 0.;
319
320 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
321 fResolHisto->Fill(vResolHisto);
322
323 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
324 fPullHisto->Fill(vPullHisto);
325 }
326}
327
328//_____________________________________________________________________________
329void AliPerformanceRes::ProcessConstrained(AliStack* const stack, AliESDtrack *const esdTrack)
330{
331 // Fill resolution comparison information (constarained parameters)
332 //
333 if(!esdTrack) return;
334
335 const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
336 if(!track) return;
337
338 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
339 esdTrack->GetImpactParameters(dca,cov);
340
341 //
342 // Fill rec vs MC information
343 //
344 if(!stack) return;
345
346 Int_t label = TMath::Abs(esdTrack->GetLabel());
347 if(label == 0) return;
348
349 TParticle* particle = stack->Particle(label);
350 if(!particle) return;
351 if(particle->GetPDG()->Charge()==0) return;
352 //printf("charge %d \n",particle->GetPDG()->Charge());
353
354 Float_t mceta = particle->Eta();
355 Float_t mcphi = particle->Phi();
356 if(mcphi<0) mcphi += 2.*TMath::Pi();
357 Float_t mcpt = particle->Pt();
358
359 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
360 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
361
362 Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
363 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
364
365 // select primaries
366 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
367 {
368 deltaYTPC= track->GetY()-particle->Vy();
369 deltaZTPC = track->GetZ()-particle->Vz();
370 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(particle->Pz(),particle->Pt());
371 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(particle->Py(),particle->Px());
372 delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
373
374 pullYTPC= (track->GetY()-particle->Vy()) / TMath::Sqrt(track->GetSigmaY2());
375 pullZTPC = (track->GetZ()-particle->Vz()) / TMath::Sqrt(track->GetSigmaZ2());
376
377 Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
378 Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
379 pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
380 pullPhiTPC = deltaPhiTPC / sigma_phi;
381
382 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
383 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
384
385 if (mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
386 else pull1PtTPC = 0.;
387
388 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
389 fResolHisto->Fill(vResolHisto);
390
391 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,particle->Vy(),particle->Vz(),mceta,mcphi,mcpt};
392 fPullHisto->Fill(vPullHisto);
393 }
394}
395
396//_____________________________________________________________________________
397void AliPerformanceRes::ProcessInnerTPC(AliMCEvent *const mcEvent, AliESDtrack *const esdTrack)
398{
399 //
400 // Fill resolution comparison information (inner params at TPC reference point)
401 //
402 if(!esdTrack) return;
403
404 const AliExternalTrackParam * innerParam = esdTrack->GetInnerParam();
405 if(!innerParam) return;
406
407 // create new AliExternalTrackParam
408 AliExternalTrackParam *track = new AliExternalTrackParam(*innerParam);
409 if(!track) return;
410
411 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
412 esdTrack->GetImpactParameters(dca,cov);
413
414 //
415 // Fill rec vs MC information
416 //
417 if(!mcEvent) return;
418
419 Int_t label = TMath::Abs(esdTrack->GetLabel());
420 if(label == 0) return;
421 AliMCParticle *mcParticle = mcEvent->GetTrack(label);
422 if(!mcParticle) return;
423
424 // get the first TPC track reference
425 AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle);
426 if(!ref0) return;
427
428 // calculate alpha angle
429 Double_t xyz[3] = {ref0->X(),ref0->Y(),ref0->Z()};
430 Double_t alpha = TMath::ATan2(xyz[1],xyz[0]);
431 //if (alpha<0) alpha += TMath::TwoPi();
432
433 // rotate inner track to local coordinate system
434 // and propagate to the radius of the first track referenco of TPC
435 Double_t trRadius = TMath::Sqrt(xyz[1] * xyz[1] + xyz[0] * xyz[0]);
436 Bool_t isOK = track->Propagate(alpha,trRadius,AliTracker::GetBz());
437 if(!isOK) return;
438
439 Float_t mceta = -TMath::Log(TMath::Tan(0.5 * ref0->Theta()));
440 Float_t mcphi = ref0->Phi();
441 if(mcphi<0) mcphi += 2.*TMath::Pi();
442 Float_t mcpt = ref0->Pt();
443
444 if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
445 if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
446
447 Float_t delta1PtTPC, deltaYTPC, deltaZTPC, deltaPhiTPC, deltaLambdaTPC;
448 Float_t pull1PtTPC, pullYTPC, pullZTPC, pullPhiTPC, pullLambdaTPC;
449
450 // select primaries
451 if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
452 {
453 //deltaYTPC= track->GetY()-ref0->Y();
454 deltaYTPC= track->GetY(); // it corresponds to deltaY after alpha rotation
455 deltaZTPC = track->GetZ()-ref0->Z();
456 deltaLambdaTPC = TMath::ATan2(track->Pz(),track->Pt())-TMath::ATan2(ref0->Pz(),ref0->Pt());
457 deltaPhiTPC = TMath::ATan2(track->Py(),track->Px())-TMath::ATan2(ref0->Py(),ref0->Px());
458 if(mcpt) delta1PtTPC = (track->OneOverPt()-1./mcpt)*mcpt;
459 else delta1PtTPC = 0.;
460
461 //pullYTPC= (track->GetY()-ref0->Y()) / TMath::Sqrt(track->GetSigmaY2());
462 pullYTPC= track->GetY() / TMath::Sqrt(track->GetSigmaY2()); //it corresponds to deltaY after alpha rotation
463 pullZTPC = (track->GetZ()-ref0->Z()) / TMath::Sqrt(track->GetSigmaZ2());
464
465 //pullLambdaTPC = deltaLambdaTPC / TMath::Sqrt(track->GetSigmaTgl2());
466 //pullPhiTPC = deltaPhiTPC / TMath::Sqrt(track->GetSigmaSnp2());
467
468 Double_t sigma_lambda = 1./(1.+track->GetTgl()*track->GetTgl()) * TMath::Sqrt(track->GetSigmaTgl2());
469 Double_t sigma_phi = 1./TMath::Sqrt(1-track->GetSnp()*track->GetSnp()) * TMath::Sqrt(track->GetSigmaSnp2());
470 pullLambdaTPC = deltaLambdaTPC / sigma_lambda;
471 pullPhiTPC = deltaPhiTPC / sigma_phi;
472
473 if(mcpt) pull1PtTPC = (track->OneOverPt()-1./mcpt) / TMath::Sqrt(track->GetSigma1Pt2());
474 else pull1PtTPC = 0.;
475
476 Double_t vResolHisto[10] = {deltaYTPC,deltaZTPC,deltaPhiTPC,deltaLambdaTPC,delta1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
477 fResolHisto->Fill(vResolHisto);
478
479 Double_t vPullHisto[10] = {pullYTPC,pullZTPC,pullPhiTPC,pullLambdaTPC,pull1PtTPC,ref0->Y(),ref0->Z(),mceta,mcphi,mcpt};
480 fPullHisto->Fill(vPullHisto);
481 }
482
483 if(track) delete track;
484}
485
486//_____________________________________________________________________________
487AliTrackReference * AliPerformanceRes::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
488{
489 // return first TPC track reference
490 if(!mcParticle) return 0;
491
492 Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
493 AliTrackReference *ref0 = 0;
494 for (Int_t iref = 0; iref < nTrackRef; iref++) {
495 ref0 = mcParticle->GetTrackReference(iref);
496 if(ref0 && (ref0->DetectorId()==AliTrackReference::kTPC))
497 break;
498 }
499
500return ref0;
501}
502
503//_____________________________________________________________________________
504void AliPerformanceRes::Exec(AliMCEvent* const mcEvent, AliESDEvent* const esdEvent, const Bool_t bUseMC)
505{
506 // Process comparison information
507 //
508 if(!esdEvent)
509 {
510 AliDebug(AliLog::kError, "esdEvent not available");
511 return;
512 }
513 AliHeader* header = 0;
514 AliGenEventHeader* genHeader = 0;
515 AliStack* stack = 0;
516 TArrayF vtxMC(3);
517
518 if(bUseMC)
519 {
520 if(!mcEvent) {
521 AliDebug(AliLog::kError, "mcEvent not available");
522 return;
523 }
524
525 // get MC event header
526 header = mcEvent->Header();
527 if (!header) {
528 AliDebug(AliLog::kError, "Header not available");
529 return;
530 }
531 // MC particle stack
532 stack = mcEvent->Stack();
533 if (!stack) {
534 AliDebug(AliLog::kError, "Stack not available");
535 return;
536 }
537
538 // get MC vertex
539 genHeader = header->GenEventHeader();
540 if (!genHeader) {
541 AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
542 return;
543 }
544 genHeader->PrimaryVertex(vtxMC);
545
546 } // end bUseMC
547
548 // Process events
549 for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
550 {
551 AliESDtrack*track = esdEvent->GetTrack(iTrack);
552 if(!track) continue;
553
554 if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
555 else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
556 else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
557 else if(GetAnalysisMode() == 3) ProcessInnerTPC(mcEvent,track);
558 else {
559 printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
560 return;
561 }
562 }
563}
564
565//_____________________________________________________________________________
566TH1F* AliPerformanceRes::MakeResol(TH2F * his, Int_t integ, Bool_t type, Int_t cut){
567 // Create resolution histograms
568
569 TH1F *hisr, *hism;
570 if (!gPad) new TCanvas;
571 hisr = AliTreeDraw::CreateResHistoII(his,&hism,integ,kTRUE,cut);
572 if (type) return hism;
573 else
574 return hisr;
575}
576
577//_____________________________________________________________________________
578void AliPerformanceRes::Analyse(){
579 // Analyse comparison information and store output histograms
580 // in the folder "folderRes"
581 //
582 TH1::AddDirectory(kFALSE);
583 TH1F *h=0;
584 TH2F *h2D=0;
585 TObjArray *aFolderObj = new TObjArray;
586
587 // write results in the folder
588 TCanvas * c = new TCanvas("Phi resol Tan","Phi resol Tan");
589 c->cd();
590
591 char name[256];
592 char title[256];
593
594 for(Int_t i=0; i<5; i++)
595 {
596 for(Int_t j=5; j<10; j++)
597 {
598 if(j!=7) fResolHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
599 fResolHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
600 if(GetAnalysisMode() == 3) fResolHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
601
602 h2D = (TH2F*)fResolHisto->Projection(i,j);
603
604 h = AliPerformanceRes::MakeResol(h2D,1,0,100);
605 sprintf(name,"h_res_%d_vs_%d",i,j);
606
607 /*
608 if(i<2) {
609 h->SetMinimum(0.);
610 h->SetMaximum(1.);
611 } else {
612 h->SetMinimum(0.);
613 h->SetMaximum(0.2);
614 }
615 */
616 h->SetName(name);
617
618 h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
619 sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(resolution)");
620 h->GetYaxis()->SetTitle(title);
621 sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
622 h->SetTitle(title);
623
624 if(j==9) h->SetBit(TH1::kLogX);
625 aFolderObj->Add(h);
626
627 h = AliPerformanceRes::MakeResol(h2D,1,1,100);
628 //h = (TH1F*)arr->At(1);
629 sprintf(name,"h_mean_res_%d_vs_%d",i,j);
630 /*
631 h->SetMinimum(-0.05);
632 h->SetMaximum(0.05);
633 if(i<2) {
634 h->SetMinimum(-0.05);
635 h->SetMaximum(0.05);
636 } else {
637 h->SetMinimum(-0.02);
638 h->SetMaximum(0.02);
639 }
640 */
641 h->SetName(name);
642
643 h->GetXaxis()->SetTitle(fResolHisto->GetAxis(j)->GetTitle());
644 sprintf(title,"%s %s",fResolHisto->GetAxis(i)->GetTitle(),"(mean)");
645 h->GetYaxis()->SetTitle(title);
646
647 sprintf(title,"%s vs %s",title,fResolHisto->GetAxis(j)->GetTitle());
648 h->SetTitle(title);
649
650 if(j==9) h->SetBit(TH1::kLogX);
651 aFolderObj->Add(h);
652
653 fResolHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
654 fResolHisto->GetAxis(9)->SetRangeUser(0.0,10.);
655
656 //
657 if(j!=7) fPullHisto->GetAxis(7)->SetRangeUser(-0.9,0.9); // eta window
658 fPullHisto->GetAxis(9)->SetRangeUser(0.16,10.); // pt threshold
659 if(GetAnalysisMode() == 3) fPullHisto->GetAxis(5)->SetRangeUser(-80.,80.); // y range
660
661 h2D = (TH2F*)fPullHisto->Projection(i,j);
662
663 h = AliPerformanceRes::MakeResol(h2D,1,0,100);
664 sprintf(name,"h_pull_%d_vs_%d",i,j);
665 //h->SetMinimum(0.);
666 //h->SetMaximum(2.);
667 h->SetName(name);
668
669 h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
670 sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(resolution)");
671 h->GetYaxis()->SetTitle(title);
672 sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
673 h->SetTitle(title);
674
675 if(j==9) h->SetBit(TH1::kLogX);
676 aFolderObj->Add(h);
677
678 h = AliPerformanceRes::MakeResol(h2D,1,1,100);
679 sprintf(name,"h_mean_pull_%d_vs_%d",i,j);
680 //h->SetMinimum(-0.2);
681 //h->SetMaximum(0.2);
682 h->SetName(name);
683
684 h->GetXaxis()->SetTitle(fPullHisto->GetAxis(j)->GetTitle());
685 sprintf(title,"%s %s",fPullHisto->GetAxis(i)->GetTitle(),"(mean)");
686 h->GetYaxis()->SetTitle(title);
687 sprintf(title,"%s vs %s",title,fPullHisto->GetAxis(j)->GetTitle());
688 h->SetTitle(title);
689
690 if(j==9) h->SetBit(TH1::kLogX);
691 aFolderObj->Add(h);
692
693 fPullHisto->GetAxis(7)->SetRangeUser(-1.5,1.5);
694 fPullHisto->GetAxis(9)->SetRangeUser(0.0,10.);
695 }
696 }
697
698 // export objects to analysis folder
699 fAnalysisFolder = ExportToFolder(aFolderObj);
700
701 // delete only TObjArray
702 if(aFolderObj) delete aFolderObj;
703}
704
705//_____________________________________________________________________________
706TFolder* AliPerformanceRes::ExportToFolder(TObjArray * array)
707{
708 // recreate folder avery time and export objects to new one
709 //
710 AliPerformanceRes * comp=this;
711 TFolder *folder = comp->GetAnalysisFolder();
712
713 TString name, title;
714 TFolder *newFolder = 0;
715 Int_t i = 0;
716 Int_t size = array->GetSize();
717
718 if(folder) {
719 // get name and title from old folder
720 name = folder->GetName();
721 title = folder->GetTitle();
722
723 // delete old one
724 delete folder;
725
726 // create new one
727 newFolder = CreateFolder(name.Data(),title.Data());
728 newFolder->SetOwner();
729
730 // add objects to folder
731 while(i < size) {
732 newFolder->Add(array->At(i));
733 i++;
734 }
735 }
736
737return newFolder;
738}
739
740//_____________________________________________________________________________
741Long64_t AliPerformanceRes::Merge(TCollection* const list)
742{
743 // Merge list of objects (needed by PROOF)
744
745 if (!list)
746 return 0;
747
748 if (list->IsEmpty())
749 return 1;
750
751 TIterator* iter = list->MakeIterator();
752 TObject* obj = 0;
753
754 // collection of generated histograms
755 Int_t count=0;
756 while((obj = iter->Next()) != 0)
757 {
758 AliPerformanceRes* entry = dynamic_cast<AliPerformanceRes*>(obj);
759 if (entry == 0) continue;
760
761 fResolHisto->Add(entry->fResolHisto);
762 fPullHisto->Add(entry->fPullHisto);
763
764 count++;
765 }
766
767return count;
768}
769
770//_____________________________________________________________________________
771TFolder* AliPerformanceRes::CreateFolder(TString name,TString title) {
772// create folder for analysed histograms
773//
774TFolder *folder = 0;
775 folder = new TFolder(name.Data(),title.Data());
776
777 return folder;
778}