1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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.
27 // Author: S.Schuchmann 11/13/2009
28 //------------------------------------------------------------------------------
32 // after running the performance task, read the file, and get component
34 TFile f("Output.root");
35 AliPerformancePtCalibMC *compObj = (AliPerformancePtCalibMC*)coutput->FindObject("AliPerformancePtCalibMC");
37 // set phi and theta bins for fitting and analyse comparison data
38 compObj->SetProjBinsTheta(thetaBins,nThetaBins,minPhi, maxPhi);
39 compObj->SetProjBinsPhi(phiBins,nPhiBins,minTheta,maxTheta);
40 compObj->SetMakeFitOption(kFALSE,exclRange,fitRange);
41 compObj->SetDoRebin(rebin);
42 //compObj->SetAnaMCOff();
44 //details see functions of this class
46 // the output histograms/graphs will be stored in the folder "folderRes"
47 compObj->GetAnalysisFolder()->ls("*");
49 // user can save whole comparison object (or only folder with anlysed histograms)
50 // in the seperate output file (e.g.)
51 TFile fout("Analysed_InvPt.root","recreate");
52 compObj->Write(); // compObj->GetAnalysisFolder()->Write();
60 #include "THnSparse.h"
65 #include "AliESDEvent.h"
66 #include "AliESDtrack.h"
67 #include "AliESDtrackCuts.h"
68 #include "AliMCEvent.h"
70 #include "AliESDfriendTrack.h"
71 #include "AliESDfriend.h"
73 #include "AliPerformancePtCalibMC.h"
74 #include "AliPerfAnalyzeInvPt.h"
79 ClassImp(AliPerformancePtCalibMC)
81 //________________________________________________________________________
82 AliPerformancePtCalibMC::AliPerformancePtCalibMC() :
83 AliPerformanceObject("AliPerformancePtCalibMC"),
84 // option parameter for AliPerformancePtCalibMC::Analyse()
97 // option parameter for user defined charge/pt shift
110 fHistInvPtPtThetaPhi(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),
128 fHistTPCMomentaPosInvPtMC(0),
129 fHistTPCMomentaNegInvPtMC(0),
130 fHistTPCMomentaPosPtMC(0),
131 fHistTPCMomentaNegPtMC(0),
132 fHistESDMomentaPosInvPtMC(0),
133 fHistESDMomentaNegInvPtMC(0),
134 fHistESDMomentaPosPtMC(0),
135 fHistESDMomentaNegPtMC(0),
146 fShift = kFALSE; // shift in charge/pt yes/no
147 fDeltaInvP = 0.00; // shift value
150 fOptTPC = kTRUE; // read TPC tracks yes/no
156 fEtaAcceptance = 0.8;
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
172 for (Int_t i=0; i<100; i++){
180 //________________________________________________________________________
181 AliPerformancePtCalibMC::AliPerformancePtCalibMC(const char *name= "AliPerformancePtCalibMC", const char *title="AliPerformancePtCalibMC"):
182 AliPerformanceObject(name,title),
183 // option parameter for AliPerformancePtCalibMC::Analyse()
196 // option parameter for user defined charge/pt shift
209 fHistInvPtPtThetaPhi(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),
227 fHistTPCMomentaPosInvPtMC(0),
228 fHistTPCMomentaNegInvPtMC(0),
229 fHistTPCMomentaPosPtMC(0),
230 fHistTPCMomentaNegPtMC(0),
231 fHistESDMomentaPosInvPtMC(0),
232 fHistESDMomentaNegInvPtMC(0),
233 fHistESDMomentaPosPtMC(0),
234 fHistESDMomentaNegPtMC(0),
245 fShift = kFALSE; // shift in charge/pt yes/no
246 fDeltaInvP = 0.00; // shift value
249 fOptTPC = kTRUE; // read TPC tracks yes/no
255 fEtaAcceptance = 0.8;
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
271 for (Int_t i=0; i<100; i++){
279 //________________________________________________________________________
280 AliPerformancePtCalibMC::~AliPerformancePtCalibMC() {
285 if(fList) delete fList;
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;
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;
319 if(fESDTrackCuts) delete fESDTrackCuts;
320 if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
326 //________________________________________________________________________
327 void AliPerformancePtCalibMC::Init()
335 fAnalysisFolder = CreateFolder("folderPt_TPC","Analysis Pt Resolution Folder");
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);
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);
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;
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};
363 fHistInvPtPtThetaPhi = new THnSparseF("fHistInvPtPtThetaPhi","1/pt vs pt vs #theta vs #phi ",invPtDims,binsInvPt,xminInvPt,xmaxInvPt);
364 fList->Add(fHistInvPtPtThetaPhi);
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);
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);
395 fHistInvPtPtThetaPhiMC = new THnSparseF("fHistInvPtPtThetaPhiMC","MC 1/pt vs pt vs #theta vs #phi ",invPtDims,binsInvPt,xminInvPt,xmaxInvPt);
396 fList->Add(fHistInvPtPtThetaPhiMC);
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);
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);
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);
427 //________________________________________________________________________
428 void AliPerformancePtCalibMC::SetPtShift(const Double_t shiftVal ) {
430 //set user defined shift in charge/pt
432 if(shiftVal) { fShift=kTRUE; fDeltaInvP = shiftVal; }
435 //________________________________________________________________________
436 void AliPerformancePtCalibMC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const /*esdFriend*/, const Bool_t /*bUseMC*/, const Bool_t /*bUseESDfriend*/)
438 //exec: read MC and esd or tpc tracks
440 AliStack* stack = NULL;
443 Printf("ERROR: Event not available");
448 fHistTrackMultiplicity->Fill(esdEvent->GetNumberOfTracks());
451 Printf("ERROR: Could not retrieve MC event");
454 stack = mcEvent->Stack();
456 Printf("ERROR: Could not retrieve stack");
461 //vertex info for cut
462 //const AliESDVertex *vtx = esdEvent->GetPrimaryVertex();
463 //if (!vtx->GetStatus()) return ;
465 if(fShift) fHistUserPtShift->Fill(fDeltaInvP);
467 // read primary vertex info
468 Double_t tPrimaryVtxPosition[3];
469 // Double_t tPrimaryVtxCov[3];
470 const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexTPC();
472 tPrimaryVtxPosition[0] = primaryVtx->GetXv();
473 tPrimaryVtxPosition[1] = primaryVtx->GetYv();
474 tPrimaryVtxPosition[2] = primaryVtx->GetZv();
476 fHistPrimaryVertexPosX->Fill(tPrimaryVtxPosition[0]);
477 fHistPrimaryVertexPosY->Fill(tPrimaryVtxPosition[1]);
478 fHistPrimaryVertexPosZ->Fill(tPrimaryVtxPosition[2]);
481 //fill histos for pt spectra and shift of transverse momentum
484 for(Int_t j = 0;j<esdEvent->GetNumberOfTracks();j++){
485 AliESDtrack *esdTrack = esdEvent->GetTrack(j);
486 if(!esdTrack) continue;
490 if(!fESDTrackCuts->AcceptTrack(esdTrack)) continue;
494 Int_t esdLabel = esdTrack->GetLabel();
495 if(esdLabel<0) continue;
496 TParticle * partMC = stack->Particle(esdLabel);
497 if (!partMC) continue;
499 // fill correlation histos MC ESD
500 Double_t pESD = esdTrack->GetP();
501 Double_t ptESD = esdTrack->GetSignedPt();
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();
509 if(fPions && !(fabs(signMC)-211<0.001)) continue;// pions only
512 if(signMC>0) signMC = 1;
515 //fill MC information
516 Double_t thetaMC = partMC->Theta();
517 Double_t phiMC = partMC->Phi();
519 Double_t momAngMC[4] = {signMC*(fabs(invPtMC)),signMC*(fabs(mcPt)),thetaMC,phiMC};
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);
525 //correlation histos MC ESD
526 fHistInvPtMCESD->Fill(signMC*(fabs(invPtMC)),1.0/ptESD);
527 fHistPtMCESD->Fill(fabs(mcPt),fabs(ptESD));
530 // fill histos TPC or ESD
532 //TPC tracks and MC tracks
533 const AliExternalTrackParam *tpcTrack = esdTrack->GetTPCInnerParam();
534 if(!tpcTrack) continue;
535 if(fabs(tpcTrack->Eta())>= fEtaAcceptance) continue;
537 Double_t signedPt = tpcTrack->GetSignedPt();
538 Double_t invPt = 0.0;
540 invPt = 1.0/signedPt;
542 fHistPtShift0->Fill(signedPt);//changed
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;
551 Double_t theta = tpcTrack->Theta();
552 Double_t phi = tpcTrack->Phi();
554 Double_t momAng[4] = {invPt,signedPt,theta,phi};
555 fHistInvPtPtThetaPhi->Fill(momAng);
557 //correlation histos MC TPC
558 fHistInvPtMCTPC->Fill(signMC*(fabs(invPtMC)),invPt);
559 fHistPtMCTPC->Fill(fabs(mcPt),fabs(signedPt));
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();
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);
575 fHistTPCMomentaNegP->Fill(fabs(pESD),fabs(pTPC));
576 fHistTPCMomentaNegPt->Fill(fabs(ptESD),fabs(signedPt));
577 fHistTPCMomentaNegInvPtMC->Fill(invPtDiffESD,invPtDiffTPC);
578 fHistTPCMomentaNegPtMC->Fill(ptDiffESD,ptDiffTPC);
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));
589 // ESD tracks and MC tracks
590 if(fabs(esdTrack->Eta())>= fEtaAcceptance) continue;
591 Double_t invPt = 0.0;
595 fHistPtShift0->Fill(ptESD);
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;
603 Double_t theta = esdTrack->Theta();
604 Double_t phi = esdTrack->Phi();
606 Double_t momAng[4] = {invPt,ptESD,theta,phi};
607 fHistInvPtPtThetaPhi->Fill(momAng);
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);
617 fHistESDMomentaNegInvPtMC->Fill(invPtdiffESD);
618 fHistESDMomentaNegPtMC->Fill(ptDiffESD);
620 fHistdedxPions->Fill(ptESD,esdTrack->GetTPCsignal());
621 fHistMomresMCESD->Fill(fabs(mcPt),(fabs(ptESD)-fabs(mcPt))/fabs(mcPt));
627 fHistTrackMultiplicityCuts->Fill(count);
631 //______________________________________________________________________________________________________________________
633 void AliPerformancePtCalibMC::Analyse()
636 // analyse charge/pt spectra in bins of theta and phi. Bins can be set by user
639 THnSparseF *copyTHnSparseTheta;
640 THnSparseF *copyTHnSparsePhi;
643 Printf("AliPerformancePtCalibMC::Analyse: analysing MC!");
644 copyTHnSparseTheta = (THnSparseF*)fHistInvPtPtThetaPhiMC->Clone("copyTHnSparseThetaMC");
645 copyTHnSparsePhi = (THnSparseF*)fHistInvPtPtThetaPhiMC->Clone("copyTHnSparsePhiMC");
648 Printf("AliPerformancePtCalibMC::Analyse: analysing ESD (reco)!");
649 copyTHnSparseTheta = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparseTheta");
650 copyTHnSparsePhi = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparsePhi");
653 copyTHnSparseTheta->GetAxis(3)->SetRangeUser(fMinPhi,fMaxPhi);
654 copyTHnSparsePhi->GetAxis(2)->SetRangeUser(fMinTheta,fMaxTheta);
656 TH2F *histInvPtTheta = (TH2F*)copyTHnSparseTheta->Projection(2,0);
657 TH2F *histInvPtPhi = (TH2F*)copyTHnSparsePhi->Projection(3,0);
659 AliPerfAnalyzeInvPt *ana = new AliPerfAnalyzeInvPt("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt");
662 TH1::AddDirectory(kFALSE);
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;
671 ana->StartAnalysis(histInvPtTheta,histInvPtPhi, aFolderObj);
673 // export objects to analysis folder
674 fAnalysisFolder = ExportToFolder(aFolderObj);
676 // delete only TObjArray
677 if(aFolderObj) delete aFolderObj;
682 //______________________________________________________________________________________________________________________
683 TFolder* AliPerformancePtCalibMC::ExportToFolder(TObjArray * array)
685 // recreate folder avery time and export objects to new one
687 AliPerformancePtCalibMC * comp=this;
688 TFolder *folder = comp->GetAnalysisFolder();
691 TFolder *newFolder = 0;
693 Int_t size = array->GetSize();
696 // get name and title from old folder
697 name = folder->GetName();
698 title = folder->GetTitle();
704 newFolder = CreateFolder(name.Data(),title.Data());
705 newFolder->SetOwner();
707 // add objects to folder
709 newFolder->Add(array->At(i));
717 //______________________________________________________________________________________________________________________
718 Long64_t AliPerformancePtCalibMC::Merge(TCollection* const list)
720 // Merge list of objects (needed by PROOF)
728 TIterator* iter = list->MakeIterator();
731 // collection of generated histograms
733 while((obj = iter->Next()) != 0)
735 AliPerformancePtCalibMC* entry = dynamic_cast<AliPerformancePtCalibMC*>(obj);
736 if (!entry) continue;
737 fHistInvPtPtThetaPhi->Add(entry->fHistInvPtPtThetaPhi);
739 fHistInvPtPtThetaPhiMC->Add(entry->fHistInvPtPtThetaPhiMC);
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);
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);
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);
774 //______________________________________________________________________________________________________________________
775 TFolder* AliPerformancePtCalibMC::CreateFolder(TString name,TString title) {
776 // create folder for analysed histograms
779 folder = new TFolder(name.Data(),title.Data());
785 // set variables for Analyse()
786 //______________________________________________________________________________________________________________________
787 void 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
794 for(Int_t k = 0;k<fNPhiBins;k++){
795 fPhiBins[k] = phiBinArray[k];
797 Printf("AliPerformancePtCalibMC::SetProjBinsPhi: number of bins in phi set to %i",fNPhiBins);
799 else Printf("Warning AliPerformancePtCalibMC::SetProjBinsPhi: number of bins in phi NOT set!!! Default values are taken.");
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.");
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;
810 //____________________________________________________________________________________________________________________________________________
811 void 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
816 fNThetaBins = nthBins;
817 for(Int_t k = 0;k<fNThetaBins;k++){
818 fThetaBins[k] = thetaBinArray[k];
820 Printf("AliPerformancePtCalibMC::SetProjBinsTheta: number of bins in theta set to %i",fNThetaBins);
822 else Printf("Warning AliPerformancePtCalibMC::SetProjBinsTheta: number of bins in theta NOT set!!! Default values are taken.");
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.");
828 Printf("AliPerformancePtCalibMC::SetProjBinsTheta: phi range for projection on theta and charge/pt is selected by user: %1.3f - %1.3f rad",minPhi,maxPhi);
834 //____________________________________________________________________________________________________________________________________________
835 void AliPerformancePtCalibMC::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
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
842 fFitGaus = setGausFit;
843 fExclRange = exclusionR;
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);