]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGPP/TPC/AliPerformancePtCalibMC.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Update T0 (Alla).
[u/mrichter/AliRoot.git] / PWGPP / TPC / AliPerformancePtCalibMC.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15//------------------------------------------------------------------------------
16// Implementation of AliPerformancePtCalibMC class. It compares ESD, TPC track
17// momenta with MC information to study charge/pt spectra.
18// The output can be analysed with AliPerfAnalyzeInvPt via AliPerformancePtCalibMC::Analyse():
19// Projection of charge/pt vs theta and vs phi resp. histoprams will be fitted with either
20// polynomial or gaussian fit function to extract minimum position of charge/pt.
21// Fit options and theta, phi bins can be set by user.
22// Attention: use the Set functions of AliPerformancePtCalibMC when running
23// AliPerformancePtCalibMC::Analyse()
24// The result of the analysis (histograms/graphs) are stored in the folder which is
25// a data member of AliPerformancePtCalibMC.
26//
27// Author: S.Schuchmann 11/13/2009
28//------------------------------------------------------------------------------
29
30/*
31
32// after running the performance task, read the file, and get component
33
34TFile f("Output.root");
35AliPerformancePtCalibMC *compObj = (AliPerformancePtCalibMC*)coutput->FindObject("AliPerformancePtCalibMC");
36
37// set phi and theta bins for fitting and analyse comparison data
38compObj->SetProjBinsTheta(thetaBins,nThetaBins,minPhi, maxPhi);
39compObj->SetProjBinsPhi(phiBins,nPhiBins,minTheta,maxTheta);
40compObj->SetMakeFitOption(kFALSE,exclRange,fitRange);
41compObj->SetDoRebin(rebin);
42//compObj->SetAnaMCOff();
43compObj->Analyse();
44//details see functions of this class
45
46// the output histograms/graphs will be stored in the folder "folderRes"
47compObj->GetAnalysisFolder()->ls("*");
48
49// user can save whole comparison object (or only folder with anlysed histograms)
50// in the seperate output file (e.g.)
51TFile fout("Analysed_InvPt.root","recreate");
52compObj->Write(); // compObj->GetAnalysisFolder()->Write();
53fout.Close();
54
55*/
56
57
58#include "TH1F.h"
59#include "TH2F.h"
60#include "THnSparse.h"
61#include "TList.h"
62#include "TMath.h"
63#include "TFolder.h"
64
65#include "AliESDEvent.h"
66#include "AliESDtrack.h"
67#include "AliESDtrackCuts.h"
68#include "AliMCEvent.h"
69#include "AliStack.h"
70#include "AliESDfriendTrack.h"
71#include "AliESDfriend.h"
72
73#include "AliPerformancePtCalibMC.h"
74#include "AliPerfAnalyzeInvPt.h"
75
76
77using namespace std;
78
79ClassImp(AliPerformancePtCalibMC)
80
81//________________________________________________________________________
82 AliPerformancePtCalibMC::AliPerformancePtCalibMC() :
83 AliPerformanceObject("AliPerformancePtCalibMC"),
84 // option parameter for AliPerformancePtCalibMC::Analyse()
85 fNThetaBins(0),
86 fNPhiBins(0),
87 fMaxPhi(0),
88 fMinPhi(0),
89 fMaxTheta(0),
90 fMinTheta(0),
91 fRange(0),
92 fExclRange(0),
93 fFitGaus(0) ,
94 fDoRebin(0),
95 fRebin(0),
96 fAnaMC(0),
97 // option parameter for user defined charge/pt shift
98 fShift(0),
99 fDeltaInvP(0),
100 //options for cuts
101 fOptTPC(0),
102 fESDcuts(0),
103 fPions(0),
104 fEtaAcceptance(0),
105 fCutsRC(0),
106 fCutsMC(0),
107 fList(0),
108
109 // histograms
110 fHistInvPtPtThetaPhi(0),
111 fHistPtShift0(0),
112 fHistPrimaryVertexPosX(0),
113 fHistPrimaryVertexPosY(0),
114 fHistPrimaryVertexPosZ(0),
115 fHistTrackMultiplicity(0),
116 fHistTrackMultiplicityCuts(0),
117 fHistTPCMomentaPosP(0),
118 fHistTPCMomentaNegP(0),
119 fHistTPCMomentaPosPt(0),
120 fHistTPCMomentaNegPt(0),
121 fHistInvPtPtThetaPhiMC(0),
122 fHistInvPtMCESD(0),
123 fHistInvPtMCTPC(0),
124 fHistPtMCESD(0),
125 fHistPtMCTPC(0),
126 fHistMomresMCESD(0),
127 fHistMomresMCTPC(0),
128 fHistTPCMomentaPosInvPtMC(0),
129 fHistTPCMomentaNegInvPtMC(0),
130 fHistTPCMomentaPosPtMC(0),
131 fHistTPCMomentaNegPtMC(0),
132 fHistESDMomentaPosInvPtMC(0),
133 fHistESDMomentaNegInvPtMC(0),
134 fHistESDMomentaPosPtMC(0),
135 fHistESDMomentaNegPtMC(0),
136 fHistUserPtShift(0),
137 fHistdedxPions(0),
138 //esd track cuts
139 fESDTrackCuts(0),
140 // analysis folder
141 fAnalysisFolder(0)
142{
143 // Dummy constructor
144
145
146 fShift = kFALSE; // shift in charge/pt yes/no
147 fDeltaInvP = 0.00; // shift value
148 //options for cuts
149
150 fOptTPC = kTRUE; // read TPC tracks yes/no
151 fESDcuts = kFALSE;
152 fPions = kFALSE;
153 fCutsRC = NULL;
154 fCutsMC = NULL;
155
156 fEtaAcceptance = 0.8;
157
158 // options for function AliPerformancePtCalibMC::Analyse()
159 fFitGaus = kFALSE;// use gaussian function for fitting charge/pt yes/no
160 fNThetaBins = 0; //number of theta bins
161 fNPhiBins = 0; //number of phi bins
162 fMaxPhi = 6.5;// max phi for 2D projection on theta and charge/pt axis
163 fMinPhi = 0.0;// min phi for 2D projection on theta and charge/pt axis
164 fMaxTheta = 3.0;// max theta for 2D projection on phi and charge/pt axis
165 fMinTheta = 0.0;// min theta for 2D projection on phi and charge/pt axis
166 fRange = 0; //fit range around 0
167 fExclRange =0; //range of rejection of points around 0
168 fAnaMC = kTRUE; // analyse MC tracks yes/no
169 fDoRebin = kFALSE;
170 fRebin = 0;
171
172 for (Int_t i=0; i<100; i++){
173 fThetaBins[i] = 0;
174 fPhiBins[i] = 0;
175 }
176
177 Init();
178}
179
180//________________________________________________________________________
181AliPerformancePtCalibMC::AliPerformancePtCalibMC(const char *name= "AliPerformancePtCalibMC", const char *title="AliPerformancePtCalibMC"):
182 AliPerformanceObject(name,title),
183 // option parameter for AliPerformancePtCalibMC::Analyse()
184 fNThetaBins(0),
185 fNPhiBins(0),
186 fMaxPhi(0),
187 fMinPhi(0),
188 fMaxTheta(0),
189 fMinTheta(0),
190 fRange(0),
191 fExclRange(0),
192 fFitGaus(0) ,
193 fDoRebin(0),
194 fRebin(0),
195 fAnaMC(0),
196 // option parameter for user defined charge/pt shift
197 fShift(0),
198 fDeltaInvP(0),
199 //options for cuts
200 fOptTPC(0),
201 fESDcuts(0),
202 fPions(0),
203 fEtaAcceptance(0),
204 fCutsRC(0),
205 fCutsMC(0),
206 fList(0),
207
208 // histograms
209 fHistInvPtPtThetaPhi(0),
210 fHistPtShift0(0),
211 fHistPrimaryVertexPosX(0),
212 fHistPrimaryVertexPosY(0),
213 fHistPrimaryVertexPosZ(0),
214 fHistTrackMultiplicity(0),
215 fHistTrackMultiplicityCuts(0),
216 fHistTPCMomentaPosP(0),
217 fHistTPCMomentaNegP(0),
218 fHistTPCMomentaPosPt(0),
219 fHistTPCMomentaNegPt(0),
220 fHistInvPtPtThetaPhiMC(0),
221 fHistInvPtMCESD(0),
222 fHistInvPtMCTPC(0),
223 fHistPtMCESD(0),
224 fHistPtMCTPC(0),
225 fHistMomresMCESD(0),
226 fHistMomresMCTPC(0),
227 fHistTPCMomentaPosInvPtMC(0),
228 fHistTPCMomentaNegInvPtMC(0),
229 fHistTPCMomentaPosPtMC(0),
230 fHistTPCMomentaNegPtMC(0),
231 fHistESDMomentaPosInvPtMC(0),
232 fHistESDMomentaNegInvPtMC(0),
233 fHistESDMomentaPosPtMC(0),
234 fHistESDMomentaNegPtMC(0),
235 fHistUserPtShift(0),
236 fHistdedxPions(0),
237 //esd track cuts
238 fESDTrackCuts(0),
239 // analysis folder
240 fAnalysisFolder(0)
241{
242 // Dummy constructor
243
244
245 fShift = kFALSE; // shift in charge/pt yes/no
246 fDeltaInvP = 0.00; // shift value
247 //options for cuts
248
249 fOptTPC = kTRUE; // read TPC tracks yes/no
250 fESDcuts = kFALSE;
251 fPions = kFALSE;
252 fCutsRC = NULL;
253 fCutsMC = NULL;
254
255 fEtaAcceptance = 0.8;
256
257 // options for function AliPerformancePtCalibMC::Analyse()
258 fFitGaus = kFALSE;// use gaussian function for fitting charge/pt yes/no
259 fNThetaBins = 0; //number of theta bins
260 fNPhiBins = 0; //number of phi bins
261 fMaxPhi = 6.5;// max phi for 2D projection on theta and charge/pt axis
262 fMinPhi = 0.0;// min phi for 2D projection on theta and charge/pt axis
263 fMaxTheta = 3.0;// max theta for 2D projection on phi and charge/pt axis
264 fMinTheta = 0.0;// min theta for 2D projection on phi and charge/pt axis
265 fRange = 0; //fit range around 0
266 fExclRange =0; //range of rejection of points around 0
267 fAnaMC = kTRUE; // analyse MC tracks yes/no
268 fDoRebin = kFALSE;// flag for rebin
269 fRebin = 0;// bins for rebin
270
271 for (Int_t i=0; i<100; i++){
272 fThetaBins[i] = 0;
273 fPhiBins[i] = 0;
274 }
275
276 Init();
277}
278
279//________________________________________________________________________
280AliPerformancePtCalibMC::~AliPerformancePtCalibMC() {
281 //
282 // destructor
283 //
284
285 if(fList) delete fList;
286 // histograms
287 if(fHistInvPtPtThetaPhi) delete fHistInvPtPtThetaPhi;fHistInvPtPtThetaPhi=0;
288 if(fHistPtShift0) delete fHistPtShift0;fHistPtShift0=0;
289 if(fHistPrimaryVertexPosX) delete fHistPrimaryVertexPosX;fHistPrimaryVertexPosX=0;
290 if(fHistPrimaryVertexPosY) delete fHistPrimaryVertexPosY;fHistPrimaryVertexPosY=0;
291 if(fHistPrimaryVertexPosZ) delete fHistPrimaryVertexPosZ;fHistPrimaryVertexPosZ=0;
292 if(fHistTrackMultiplicity) delete fHistTrackMultiplicity;fHistTrackMultiplicity=0;
293 if(fHistTrackMultiplicityCuts) delete fHistTrackMultiplicityCuts;fHistTrackMultiplicityCuts=0;
294
295 if(fHistTPCMomentaPosP) delete fHistTPCMomentaPosP;fHistTPCMomentaPosP=0;
296 if(fHistTPCMomentaNegP) delete fHistTPCMomentaNegP;fHistTPCMomentaNegP=0;
297 if(fHistTPCMomentaPosPt) delete fHistTPCMomentaPosPt;fHistTPCMomentaPosPt=0;
298 if(fHistTPCMomentaNegPt) delete fHistTPCMomentaNegPt ;fHistTPCMomentaNegPt=0;
299 if(fHistUserPtShift) delete fHistUserPtShift;fHistUserPtShift=0;
300 if(fHistInvPtPtThetaPhiMC) delete fHistInvPtPtThetaPhiMC;fHistInvPtPtThetaPhiMC=0;
301 if(fHistInvPtMCESD) delete fHistInvPtMCESD;fHistInvPtMCESD=0;
302 if(fHistInvPtMCTPC) delete fHistInvPtMCTPC;fHistInvPtMCTPC=0;
303 if(fHistPtMCESD) delete fHistPtMCESD;fHistPtMCESD=0;
304 if(fHistPtMCTPC) delete fHistPtMCTPC;fHistPtMCTPC=0;
305 if(fHistMomresMCESD) delete fHistMomresMCESD;fHistMomresMCESD=0;
306 if(fHistMomresMCTPC) delete fHistMomresMCTPC;fHistMomresMCTPC=0;
307 if(fHistTPCMomentaPosInvPtMC) delete fHistTPCMomentaPosInvPtMC;fHistTPCMomentaPosInvPtMC=0;
308 if(fHistTPCMomentaNegInvPtMC) delete fHistTPCMomentaNegInvPtMC;fHistTPCMomentaNegInvPtMC=0;
309 if(fHistTPCMomentaPosPtMC) delete fHistTPCMomentaPosPtMC;fHistTPCMomentaPosPtMC=0;
310 if(fHistTPCMomentaNegPtMC) delete fHistTPCMomentaNegPtMC;fHistTPCMomentaNegPtMC=0;
311 if(fHistESDMomentaPosInvPtMC) delete fHistESDMomentaPosInvPtMC;fHistESDMomentaPosInvPtMC=0;
312 if(fHistESDMomentaNegInvPtMC) delete fHistESDMomentaNegInvPtMC;fHistESDMomentaNegInvPtMC=0;
313 if(fHistESDMomentaPosPtMC) delete fHistESDMomentaPosPtMC;fHistESDMomentaPosPtMC=0;
314 if(fHistESDMomentaNegPtMC) delete fHistESDMomentaNegPtMC;fHistESDMomentaNegPtMC=0;
315 if(fHistdedxPions) delete fHistdedxPions;fHistdedxPions=0;
316
317
318 //esd track cuts
319 if(fESDTrackCuts) delete fESDTrackCuts;
320 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
321
322
323
324}
325
326//________________________________________________________________________
327void AliPerformancePtCalibMC::Init()
328{
329 // Create histograms
330 // Called once
331
332 fList = new TList();
333
334 // init folder
335 fAnalysisFolder = CreateFolder("folderPt_TPC","Analysis Pt Resolution Folder");
336
337 // Primary Vertex:
338 fHistPrimaryVertexPosX = new TH1F("fHistPrimaryVertexPosX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
339 fList->Add(fHistPrimaryVertexPosX);
340 fHistPrimaryVertexPosY = new TH1F("fHistPrimaryVertexPosY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
341 fList->Add(fHistPrimaryVertexPosY);
342 fHistPrimaryVertexPosZ = new TH1F("fHistPrimaryVertexPosZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-2.0,2.0);
343 fList->Add(fHistPrimaryVertexPosZ);
344
345 // Multiplicity:
346 fHistTrackMultiplicity = new TH1F("fHistTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 250, 0, 250);
347 fList->Add(fHistTrackMultiplicity);
348 fHistTrackMultiplicityCuts = new TH1F("fHistTrackMultiplicityCuts", "Multiplicity distribution;Number of tracks after cuts;Events", 250, 0, 250);
349 fList->Add(fHistTrackMultiplicityCuts);
350
351 // momentum histos
352 fHistPtShift0 = new TH1F("fHistPtShift0","1/pt dN/pt vs. pt of ESD track ",800,-20.0,20.0);
353 fList->Add(fHistPtShift0);
354 const Int_t invPtDims = 4;
355 fMaxPhi=6.52;
356 fMinPhi=0.0;
357 fMaxTheta=3.0;
358 fMinTheta=0.0;
359 Double_t xminInvPt[invPtDims] = {-4.5,-20.0,fMinTheta,fMinPhi};
360 Double_t xmaxInvPt[invPtDims] = {4.5,20.0,fMaxTheta,fMaxPhi};
361 Int_t binsInvPt[invPtDims] = {450,400,150,163};
362
363 fHistInvPtPtThetaPhi = new THnSparseF("fHistInvPtPtThetaPhi","1/pt vs pt vs #theta vs #phi ",invPtDims,binsInvPt,xminInvPt,xmaxInvPt);
364 fList->Add(fHistInvPtPtThetaPhi);
365
366 // momentum test histos
367 fHistTPCMomentaPosP = new TH2F("fHistTPCMomentaPosP","TPC p vs global esd track p pos",300,0.0,15.0,300,0.0,15.0);
368 fList->Add(fHistTPCMomentaPosP);
369 fHistTPCMomentaNegP = new TH2F("fHistTPCMomentaNegP","TPC p vs global esd track p neg",300,0.0,15.0,300,0.0,15.0);
370 fList->Add(fHistTPCMomentaNegP);
371 fHistTPCMomentaPosPt = new TH2F("fHistTPCMomentaPosPt","TPC pt vs global esd track pt pos",300,0.0,15.0,300,0.0,15.0);
372 fList->Add(fHistTPCMomentaPosPt);
373 fHistTPCMomentaNegPt = new TH2F("fHistTPCMomentaNegPt","TPC pt vs global esd track pt neg",300,0.0,15.0,300,0.0,15.0);
374 fList->Add(fHistTPCMomentaNegPt);
375
376 // momentum test histos MC
377 fHistTPCMomentaPosInvPtMC = new TH2F("fHistTPCMomentaPosInvPtMC","TPC-MC of 1/pt vs global ESD-MC of 1/pt pos",250, -5.0, 5.0,250, -5.0,5.0);
378 fList->Add(fHistTPCMomentaPosInvPtMC);
379 fHistTPCMomentaNegInvPtMC = new TH2F("fHistTPCMomentaNegInvPtMC","TPC-MC of 1/pt vs global ESD-MC 1/pt neg",250, -5.0, 5.0,250, -5.0, 5.0);
380 fList->Add(fHistTPCMomentaNegInvPtMC);
381 fHistTPCMomentaPosPtMC = new TH2F("fHistTPCMomentaPosPtMC","(TPC-MC)/MC^2 of pt vs global (ESD-MC)/MC^2 of pt pos",200,-2.0,2.0,200,-2.0,2.0);
382 fList->Add(fHistTPCMomentaPosPtMC);
383 fHistTPCMomentaNegPtMC = new TH2F("fHistTPCMomentaNegPtMC","(TPC-MC/)MC^2 of pt vs global (ESD-MC)/MC^2 of pt neg",200,-2.0,2.0,200,-2.0,2.0);
384 fList->Add(fHistTPCMomentaNegPtMC);
385 fHistESDMomentaPosInvPtMC = new TH1F("fHistESDMomentaPosInvPtMC","ESD-MC of 1/pt pos ",200, -2.0, 2.0);
386 fList->Add(fHistESDMomentaPosInvPtMC);
387 fHistESDMomentaNegInvPtMC = new TH1F("fHistESDMomentaNegInvPtMC","ESD-MC of 1/pt neg",200, -2.0, 2.0);
388 fList->Add(fHistESDMomentaNegInvPtMC);
389 fHistESDMomentaPosPtMC = new TH1F("fHistESDMomentaPosPtMC","(ESD-MC)/MC^2 of pt pos",200,-2.0,2.0);
390 fList->Add(fHistESDMomentaPosPtMC);
391 fHistESDMomentaNegPtMC = new TH1F("fHistESDMomentaNegPtMC","(ESD-MC)/MC^2 of pt neg",200,-2.0,2.0);
392 fList->Add(fHistESDMomentaNegPtMC);
393
394 // MC only info
395 fHistInvPtPtThetaPhiMC = new THnSparseF("fHistInvPtPtThetaPhiMC","MC 1/pt vs pt vs #theta vs #phi ",invPtDims,binsInvPt,xminInvPt,xmaxInvPt);
396 fList->Add(fHistInvPtPtThetaPhiMC);
397
398
399 //correlation histos MC ESD or TPC
400 fHistInvPtMCESD = new TH2F("fHistInvPtMCESD","inv pt ESD vs MC",450, -4.5, 4.5,450, -4.5, 4.5);
401 fList->Add(fHistInvPtMCESD);
402 fHistPtMCESD = new TH2F("fHistPtMCESD"," pt ESD vs MC",500, 0.0, 50.0,500, 0.0, 50.0);
403 fList->Add(fHistPtMCESD);
404 fHistInvPtMCTPC = new TH2F("fHistInvPtMCTPC","inv pt TPC vs MC",450, -4.5, 4.5,450, -4.5, 4.5);
405 fList->Add(fHistInvPtMCTPC);
406 fHistPtMCTPC = new TH2F("fHistPtMCTPC"," pt TPC vs MC",500, 0.0, 50.0,500, 0.0,50.0);
407 fList->Add(fHistPtMCTPC);
408 fHistMomresMCESD = new TH2F("fHistMomresMCESD"," (pt ESD - pt MC)/ptMC vs pt MC",500, 0.0, 50.0,200, -2.0, 2.0);
409 fList->Add(fHistMomresMCESD);
410 fHistMomresMCTPC = new TH2F("fHistMomresMCTPC"," (pt TPC - pt MC)/ptMC vs pt MC",500, 0.0, 50.0,200, -2.0, 2.0);
411 fList->Add(fHistMomresMCTPC);
412
413
414 //user pt shift check
415 fHistUserPtShift = new TH1F("fHistUserPtShift","user defined shift in 1/pt",100,-0.5,1.5);
416 fList->Add(fHistUserPtShift);
417 // pid
418 fHistdedxPions = new TH2F ("fHistdedxPions","dEdx of pions ident. via PDG code vs signed Pt",300,-15.05,15.05,200,0.0,400.0);
419 fList->Add(fHistdedxPions);
420
421 // esd track cuts
422 fESDTrackCuts =NULL;
423
424
425}
426
427//________________________________________________________________________
428void AliPerformancePtCalibMC::SetPtShift(const Double_t shiftVal ) {
429
430 //set user defined shift in charge/pt
431
432 if(shiftVal) { fShift=kTRUE; fDeltaInvP = shiftVal; }
433}
434
435//________________________________________________________________________
436void AliPerformancePtCalibMC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const /*esdFriend*/, const Bool_t /*bUseMC*/, const Bool_t /*bUseESDfriend*/)
437{
438 //exec: read MC and esd or tpc tracks
439
440 AliStack* stack = NULL;
441
442 if (!esdEvent) {
443 Printf("ERROR: Event not available");
444 return;
445 }
446
447
448 fHistTrackMultiplicity->Fill(esdEvent->GetNumberOfTracks());
449
450 if (!mcEvent) {
451 Printf("ERROR: Could not retrieve MC event");
452 return;
453 }
454 stack = mcEvent->Stack();
455 if (!stack) {
456 Printf("ERROR: Could not retrieve stack");
457 return;
458 }
459
460
461 //vertex info for cut
462 //const AliESDVertex *vtx = esdEvent->GetPrimaryVertex();
463 //if (!vtx->GetStatus()) return ;
464
465 if(fShift) fHistUserPtShift->Fill(fDeltaInvP);
466
467 // read primary vertex info
468 Double_t tPrimaryVtxPosition[3];
469 // Double_t tPrimaryVtxCov[3];
470 const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexTPC();
471
472 tPrimaryVtxPosition[0] = primaryVtx->GetXv();
473 tPrimaryVtxPosition[1] = primaryVtx->GetYv();
474 tPrimaryVtxPosition[2] = primaryVtx->GetZv();
475
476 fHistPrimaryVertexPosX->Fill(tPrimaryVtxPosition[0]);
477 fHistPrimaryVertexPosY->Fill(tPrimaryVtxPosition[1]);
478 fHistPrimaryVertexPosZ->Fill(tPrimaryVtxPosition[2]);
479
480
481 //fill histos for pt spectra and shift of transverse momentum
482 Int_t count=0;
483
484 for(Int_t j = 0;j<esdEvent->GetNumberOfTracks();j++){
485 AliESDtrack *esdTrack = esdEvent->GetTrack(j);
486 if(!esdTrack) continue;
487
488 //esd track cuts
489 if(fESDcuts){
490 if(!fESDTrackCuts->AcceptTrack(esdTrack)) continue;
491 }
492
493 // get MC info
494 Int_t esdLabel = esdTrack->GetLabel();
495 if(esdLabel<0) continue;
496 TParticle * partMC = stack->Particle(esdLabel);
497 if (!partMC) continue;
498
499 // fill correlation histos MC ESD
500 Double_t pESD = esdTrack->GetP();
501 Double_t ptESD = esdTrack->GetSignedPt();
502
503 if(!ptESD || !(partMC->Pt()) ) continue;
504 Double_t mcPt = partMC->Pt();
505 Double_t invPtMC = 1.0/mcPt;
506 Int_t signMC = partMC->GetPdgCode();
507
508 //pid
509 if(fPions && !(fabs(signMC)-211<0.001)) continue;// pions only
510
511 //MC only
512 if(signMC>0) signMC = 1;
513 else signMC = -1;
514
515 //fill MC information
516 Double_t thetaMC = partMC->Theta();
517 Double_t phiMC = partMC->Phi();
518
519 Double_t momAngMC[4] = {signMC*(fabs(invPtMC)),signMC*(fabs(mcPt)),thetaMC,phiMC};
520
521 // fill only if MC track is in eta acceptance of TPC in order to be compareable to TPC tracks
522 if(fabs( partMC->Eta())< fEtaAcceptance) {
523 fHistInvPtPtThetaPhiMC->Fill(momAngMC);
524
525 //correlation histos MC ESD
526 fHistInvPtMCESD->Fill(signMC*(fabs(invPtMC)),1.0/ptESD);
527 fHistPtMCESD->Fill(fabs(mcPt),fabs(ptESD));
528 }
529
530 // fill histos TPC or ESD
531 if(fOptTPC){
532 //TPC tracks and MC tracks
533 const AliExternalTrackParam *tpcTrack = esdTrack->GetTPCInnerParam();
534 if(!tpcTrack) continue;
535 if(fabs(tpcTrack->Eta())>= fEtaAcceptance) continue;
536
537 Double_t signedPt = tpcTrack->GetSignedPt();
538 Double_t invPt = 0.0;
539 if(signedPt) {
540 invPt = 1.0/signedPt;
541
542 fHistPtShift0->Fill(signedPt);//changed
543
544 if(fShift){Printf("user shift of momentum SET to non zero value!");
545 invPt += fDeltaInvP; //shift momentum for tests
546 if(invPt) signedPt = 1.0/invPt;
547 else continue;
548 }
549
550
551 Double_t theta = tpcTrack->Theta();
552 Double_t phi = tpcTrack->Phi();
553
554 Double_t momAng[4] = {invPt,signedPt,theta,phi};
555 fHistInvPtPtThetaPhi->Fill(momAng);
556
557 //correlation histos MC TPC
558 fHistInvPtMCTPC->Fill(signMC*(fabs(invPtMC)),invPt);
559 fHistPtMCTPC->Fill(fabs(mcPt),fabs(signedPt));
560
561 //compare to MC info
562 Double_t ptDiffESD = (fabs(ptESD)-fabs(mcPt))/pow(mcPt,2);
563 Double_t ptDiffTPC = (fabs(signedPt)-fabs(mcPt))/pow(mcPt,2);
564 Double_t invPtDiffESD = fabs(1.0/ptESD)-1.0/fabs(mcPt);
565 Double_t invPtDiffTPC = fabs(invPt)-1.0/fabs(mcPt);
566 Double_t pTPC = tpcTrack->GetP();
567
568 if(esdTrack->GetSign()>0){//compare momenta ESD track and TPC track
569 fHistTPCMomentaPosP->Fill(fabs(pESD),fabs(pTPC));
570 fHistTPCMomentaPosPt->Fill(fabs(ptESD),fabs(signedPt));
571 fHistTPCMomentaPosInvPtMC->Fill(invPtDiffESD,invPtDiffTPC);
572 fHistTPCMomentaPosPtMC->Fill(ptDiffESD,ptDiffTPC);
573 }
574 else{
575 fHistTPCMomentaNegP->Fill(fabs(pESD),fabs(pTPC));
576 fHistTPCMomentaNegPt->Fill(fabs(ptESD),fabs(signedPt));
577 fHistTPCMomentaNegInvPtMC->Fill(invPtDiffESD,invPtDiffTPC);
578 fHistTPCMomentaNegPtMC->Fill(ptDiffESD,ptDiffTPC);
579 }
580 fHistdedxPions->Fill(signedPt,esdTrack->GetTPCsignal());
581 fHistMomresMCESD->Fill(fabs(mcPt),(fabs(ptESD)-fabs(mcPt))/fabs(mcPt));
582 fHistMomresMCTPC->Fill(fabs(mcPt),(fabs(signedPt)-fabs(mcPt))/fabs(mcPt));
583 count++;
584 }
585 else continue;
586 }
587
588 else{
589 // ESD tracks and MC tracks
590 if(fabs(esdTrack->Eta())>= fEtaAcceptance) continue;
591 Double_t invPt = 0.0;
592
593 if(ptESD) {
594 invPt = 1.0/ptESD;
595 fHistPtShift0->Fill(ptESD);
596
597 if(fShift){Printf("user shift of momentum SET to non zero value!");
598 invPt += fDeltaInvP; //shift momentum for tests
599 if(invPt) ptESD = 1.0/invPt;
600 else continue;
601 }
602
603 Double_t theta = esdTrack->Theta();
604 Double_t phi = esdTrack->Phi();
605
606 Double_t momAng[4] = {invPt,ptESD,theta,phi};
607 fHistInvPtPtThetaPhi->Fill(momAng);
608
609 //differences MC ESD tracks
610 Double_t ptDiffESD = (fabs(ptESD)-fabs(mcPt))/pow(mcPt,2);
611 Double_t invPtdiffESD = fabs(1.0/ptESD)-1.0/fabs(mcPt);
612 if(esdTrack->GetSign()>0){
613 fHistESDMomentaPosInvPtMC->Fill(invPtdiffESD);
614 fHistESDMomentaPosPtMC->Fill(ptDiffESD);
615 }
616 else{
617 fHistESDMomentaNegInvPtMC->Fill(invPtdiffESD);
618 fHistESDMomentaNegPtMC->Fill(ptDiffESD);
619 }
620 fHistdedxPions->Fill(ptESD,esdTrack->GetTPCsignal());
621 fHistMomresMCESD->Fill(fabs(mcPt),(fabs(ptESD)-fabs(mcPt))/fabs(mcPt));
622 count++;
623 }
624 }
625}
626
627fHistTrackMultiplicityCuts->Fill(count);
628
629}
630
631//______________________________________________________________________________________________________________________
632
633void AliPerformancePtCalibMC::Analyse()
634{
635
636 // analyse charge/pt spectra in bins of theta and phi. Bins can be set by user
637
638
639 THnSparseF *copyTHnSparseTheta;
640 THnSparseF *copyTHnSparsePhi;
641
642 if(fAnaMC){
643 Printf("AliPerformancePtCalibMC::Analyse: analysing MC!");
644 copyTHnSparseTheta = (THnSparseF*)fHistInvPtPtThetaPhiMC->Clone("copyTHnSparseThetaMC");
645 copyTHnSparsePhi = (THnSparseF*)fHistInvPtPtThetaPhiMC->Clone("copyTHnSparsePhiMC");
646 }
647 else {
648 Printf("AliPerformancePtCalibMC::Analyse: analysing ESD (reco)!");
649 copyTHnSparseTheta = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparseTheta");
650 copyTHnSparsePhi = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparsePhi");
651 }
652
653 copyTHnSparseTheta->GetAxis(3)->SetRangeUser(fMinPhi,fMaxPhi);
654 copyTHnSparsePhi->GetAxis(2)->SetRangeUser(fMinTheta,fMaxTheta);
655
656 TH2F *histInvPtTheta = (TH2F*)copyTHnSparseTheta->Projection(2,0);
657 TH2F *histInvPtPhi = (TH2F*)copyTHnSparsePhi->Projection(3,0);
658
659 AliPerfAnalyzeInvPt *ana = new AliPerfAnalyzeInvPt("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt");
660 if(!ana) return;
661
662 TH1::AddDirectory(kFALSE);
663
664 ana->SetProjBinsTheta(fThetaBins,fNThetaBins);
665 ana->SetProjBinsPhi(fPhiBins,fNPhiBins);
666 ana->SetMakeFitOption(fFitGaus,fExclRange,fRange);
667 if(fDoRebin) ana->SetDoRebin(fRebin);
668 TObjArray *aFolderObj = new TObjArray;
669 if(!aFolderObj) return;
670
671 ana->StartAnalysis(histInvPtTheta,histInvPtPhi, aFolderObj);
672
673 // export objects to analysis folder
674 fAnalysisFolder = ExportToFolder(aFolderObj);
675
676 // delete only TObjArray
677 if(aFolderObj) delete aFolderObj;
678 if(ana) delete ana;
679
680}
681
682//______________________________________________________________________________________________________________________
683TFolder* AliPerformancePtCalibMC::ExportToFolder(TObjArray * array)
684{
685 // recreate folder avery time and export objects to new one
686 //
687 AliPerformancePtCalibMC * comp=this;
688 TFolder *folder = comp->GetAnalysisFolder();
689
690 TString name, title;
691 TFolder *newFolder = 0;
692 Int_t i = 0;
693 Int_t size = array->GetSize();
694
695 if(folder) {
696 // get name and title from old folder
697 name = folder->GetName();
698 title = folder->GetTitle();
699
700 // delete old one
701 delete folder;
702
703 // create new one
704 newFolder = CreateFolder(name.Data(),title.Data());
705 newFolder->SetOwner();
706
707 // add objects to folder
708 while(i < size) {
709 newFolder->Add(array->At(i));
710 i++;
711 }
712 }
713
714 return newFolder;
715}
716
717//______________________________________________________________________________________________________________________
718Long64_t AliPerformancePtCalibMC::Merge(TCollection* const list)
719{
720 // Merge list of objects (needed by PROOF)
721
722 if (!list)
723 return 0;
724
725 if (list->IsEmpty())
726 return 1;
727
728 TIterator* iter = list->MakeIterator();
729 TObject* obj = 0;
730
731 // collection of generated histograms
732 Int_t count=0;
733 while((obj = iter->Next()) != 0)
734 {
735 AliPerformancePtCalibMC* entry = dynamic_cast<AliPerformancePtCalibMC*>(obj);
736 if (!entry) continue;
737 fHistInvPtPtThetaPhi->Add(entry->fHistInvPtPtThetaPhi);
738
739 fHistInvPtPtThetaPhiMC->Add(entry->fHistInvPtPtThetaPhiMC);
740
741 fHistInvPtMCESD->Add(entry->fHistInvPtMCESD);
742 fHistPtMCESD->Add(entry->fHistPtMCESD);
743 fHistInvPtMCTPC->Add(entry->fHistInvPtMCTPC);
744 fHistPtMCTPC->Add(entry->fHistPtMCTPC);
745 fHistMomresMCESD->Add(entry->fHistMomresMCESD);
746 fHistMomresMCTPC->Add(entry->fHistMomresMCTPC);
747
748 fHistPtShift0->Add(entry->fHistPtShift0);
749 fHistPrimaryVertexPosX->Add(entry->fHistPrimaryVertexPosX);
750 fHistPrimaryVertexPosY->Add(entry->fHistPrimaryVertexPosY);
751 fHistPrimaryVertexPosZ->Add(entry->fHistPrimaryVertexPosZ);
752 fHistTrackMultiplicity->Add(entry->fHistTrackMultiplicity);
753 fHistTrackMultiplicityCuts->Add(entry->fHistTrackMultiplicityCuts);
754
755 fHistTPCMomentaPosP->Add(entry->fHistTPCMomentaPosP);
756 fHistTPCMomentaNegP->Add(entry->fHistTPCMomentaNegP);
757 fHistTPCMomentaPosPt->Add(entry->fHistTPCMomentaPosPt);
758 fHistTPCMomentaNegPt->Add(entry->fHistTPCMomentaNegPt);
759 fHistTPCMomentaPosInvPtMC->Add(entry->fHistTPCMomentaPosInvPtMC);
760 fHistTPCMomentaNegInvPtMC->Add(entry->fHistTPCMomentaNegInvPtMC);
761 fHistTPCMomentaPosPtMC->Add(entry->fHistTPCMomentaPosPtMC);
762 fHistTPCMomentaNegPtMC->Add(entry->fHistTPCMomentaNegPtMC);
763 fHistESDMomentaPosInvPtMC->Add(entry->fHistESDMomentaPosInvPtMC);
764 fHistESDMomentaNegInvPtMC->Add(entry->fHistESDMomentaNegInvPtMC);
765 fHistESDMomentaPosPtMC->Add(entry->fHistESDMomentaPosPtMC);
766 fHistESDMomentaNegPtMC->Add(entry->fHistESDMomentaNegPtMC);
767 fHistdedxPions->Add(entry->fHistdedxPions);
768 count++;
769 }
770
771 return count;
772}
773
774//______________________________________________________________________________________________________________________
775TFolder* AliPerformancePtCalibMC::CreateFolder(TString name,TString title) {
776 // create folder for analysed histograms
777 //
778 TFolder *folder = 0;
779 folder = new TFolder(name.Data(),title.Data());
780
781 return folder;
782}
783
784
785// set variables for Analyse()
786//______________________________________________________________________________________________________________________
787void AliPerformancePtCalibMC::SetProjBinsPhi(const Double_t *phiBinArray,const Int_t nphBins, const Double_t minTheta, const Double_t maxTheta){
788 // set phi bins for Analyse()
789 // set phi bins as array and set number of this array which is equal to number of bins analysed
790 // the last analysed bin will always be the projection from first to last bin in the array
791 if(nphBins){
792 fNPhiBins = nphBins;
793
794 for(Int_t k = 0;k<fNPhiBins;k++){
795 fPhiBins[k] = phiBinArray[k];
796 }
797 Printf("AliPerformancePtCalibMC::SetProjBinsPhi: number of bins in phi set to %i",fNPhiBins);
798 }
799 else Printf("Warning AliPerformancePtCalibMC::SetProjBinsPhi: number of bins in phi NOT set!!! Default values are taken.");
800
801 if(fabs(minTheta-maxTheta)<0.001){
802 Printf("AliPerformancePtCalibMC::SetProjBinsPhi: theta range for projection for projection on phi and charge/pt is too small. whole range of theta selected.");
803 }
804 else{
805 Printf("AliPerformancePtCalibMC::SetProjBinsPhi: theta range for projection on phi and charge/pt is selected by user: %1.3f - %1.3f rad",minTheta,maxTheta);
806 fMinTheta = minTheta;
807 fMaxTheta = maxTheta;
808 }
809}
810//____________________________________________________________________________________________________________________________________________
811void AliPerformancePtCalibMC::SetProjBinsTheta(const Double_t *thetaBinArray, const Int_t nthBins, const Double_t minPhi, const Double_t maxPhi){
812 // set theta bins for Analyse()
813 // set theta bins as array and set number of this array which is equal to number of bins analysed
814 // the last analysed bin will always be the projection from first to last bin in the array
815 if(nthBins){
816 fNThetaBins = nthBins;
817 for(Int_t k = 0;k<fNThetaBins;k++){
818 fThetaBins[k] = thetaBinArray[k];
819 }
820 Printf("AliPerformancePtCalibMC::SetProjBinsTheta: number of bins in theta set to %i",fNThetaBins);
821 }
822 else Printf("Warning AliPerformancePtCalibMC::SetProjBinsTheta: number of bins in theta NOT set!!! Default values are taken.");
823
824 if(fabs(minPhi-maxPhi)<0.001){
825 Printf("AliPerformancePtCalibMC::SetProjBinsTheta: phi range for projection for projection on theta and charge/pt is too small. whole range of phi selected.");
826 }
827 else{
828 Printf("AliPerformancePtCalibMC::SetProjBinsTheta: phi range for projection on theta and charge/pt is selected by user: %1.3f - %1.3f rad",minPhi,maxPhi);
829 fMinPhi = minPhi;
830 fMaxPhi = maxPhi;
831 }
832}
833
834//____________________________________________________________________________________________________________________________________________
835void AliPerformancePtCalibMC::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
836
837 // set the fit options:
838 // for usage of gaussian function instead of polynomial (default) set setGausFit=kTRUE
839 // set the range of rejection of points around 0 via exclusionR
840 // set the fit range around 0 with fitR
841
842 fFitGaus = setGausFit;
843 fExclRange = exclusionR;
844 fRange = fitR;
845
846 if(fFitGaus) Printf("AliPerformancePtCalibMC::SetMakeFitOption: set MakeGausFit with fit range %2.3f and exclusion range in fabs(1/pt): %2.3f",fRange,fExclRange);
847 else Printf("AliPerformancePtCalibMC::SetMakeFitOption: set standard polynomial fit with fit range %2.3f and exclusion range in fabs(1/pt): %2.3f",fRange,fExclRange);
848
849}