PWG1/TPC/AliPerformancePtCalibMC.cxx

   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
  34 TFile f("Output.root");
  35 AliPerformancePtCalibMC *compObj = (AliPerformancePtCalibMC*)coutput->FindObject("AliPerformancePtCalibMC");
  36
  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();
  43 compObj->Analyse();
  44 //details see functions of this class
  45
  46 // the output histograms/graphs will be stored in the folder "folderRes"
  47 compObj->GetAnalysisFolder()->ls("*");
  48
  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();
  53 fout.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
  77 using namespace std;
  78
  79 ClassImp(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    Init();
 173 }
 174
 175 //________________________________________________________________________
 176 AliPerformancePtCalibMC::AliPerformancePtCalibMC(const char *name= "AliPerformancePtCalibMC", const char *title="AliPerformancePtCalibMC"):
 177    AliPerformanceObject(name,title),
 178    // option parameter for AliPerformancePtCalibMC::Analyse()
 179    fNThetaBins(0),
 180    fNPhiBins(0),
 181    fMaxPhi(0),
 182    fMinPhi(0),
 183    fMaxTheta(0),
 184    fMinTheta(0),
 185    fRange(0),
 186    fExclRange(0),
 187    fFitGaus(0) ,
 188    fDoRebin(0),
 189    fRebin(0),
 190    fAnaMC(0),
 191    // option parameter for user defined charge/pt shift
 192    fShift(0),
 193    fDeltaInvP(0),
 194    //options for cuts
 195    fOptTPC(0),
 196    fESDcuts(0),
 197    fPions(0),
 198    fEtaAcceptance(0),
 199    fCutsRC(0),
 200    fCutsMC(0),
 201    fList(0),
 202
 203    // histograms
 204    fHistInvPtPtThetaPhi(0),
 205    fHistPtShift0(0),
 206    fHistPrimaryVertexPosX(0),
 207    fHistPrimaryVertexPosY(0),
 208    fHistPrimaryVertexPosZ(0),
 209    fHistTrackMultiplicity(0),
 210    fHistTrackMultiplicityCuts(0),
 211    fHistTPCMomentaPosP(0),
 212    fHistTPCMomentaNegP(0),
 213    fHistTPCMomentaPosPt(0),
 214    fHistTPCMomentaNegPt(0),
 215    fHistInvPtPtThetaPhiMC(0),
 216    fHistInvPtMCESD(0),
 217    fHistInvPtMCTPC(0),
 218    fHistPtMCESD(0),
 219    fHistPtMCTPC(0),
 220    fHistMomresMCESD(0),
 221    fHistMomresMCTPC(0),
 222    fHistTPCMomentaPosInvPtMC(0),
 223    fHistTPCMomentaNegInvPtMC(0),
 224    fHistTPCMomentaPosPtMC(0),
 225    fHistTPCMomentaNegPtMC(0),
 226    fHistESDMomentaPosInvPtMC(0),
 227    fHistESDMomentaNegInvPtMC(0),
 228    fHistESDMomentaPosPtMC(0),
 229    fHistESDMomentaNegPtMC(0),
 230    fHistUserPtShift(0),
 231    fHistdedxPions(0),
 232    //esd track cuts
 233    fESDTrackCuts(0),
 234    // analysis folder
 235    fAnalysisFolder(0)
 236 {
 237    // Dummy constructor
 238
 239
 240    fShift = kFALSE;                       // shift in charge/pt yes/no
 241    fDeltaInvP = 0.00;                     // shift value
 242    //options for cuts
 243
 244    fOptTPC =  kTRUE;                      // read TPC tracks yes/no
 245    fESDcuts = kFALSE;
 246    fPions = kFALSE;
 247    fCutsRC = NULL;
 248    fCutsMC = NULL;
 249
 250    fEtaAcceptance = 0.8;
 251
 252    // options for function AliPerformancePtCalibMC::Analyse()
 253    fFitGaus = kFALSE;// use gaussian function for fitting charge/pt yes/no
 254    fNThetaBins = 0; //number of theta bins
 255    fNPhiBins = 0; //number of phi bins
 256    fMaxPhi = 6.5;// max phi for 2D projection on theta and charge/pt axis
 257    fMinPhi = 0.0;// min phi for 2D projection on theta and charge/pt axis
 258    fMaxTheta = 3.0;// max theta for 2D projection on phi and charge/pt axis
 259    fMinTheta = 0.0;// min theta for 2D projection on phi and charge/pt axis
 260    fRange = 0; //fit range around 0
 261    fExclRange =0; //range of rejection of points around 0
 262    fAnaMC = kTRUE; // analyse MC tracks yes/no
 263    fDoRebin = kFALSE;// flag for rebin
 264    fRebin = 0;// bins for rebin
 265
 266    Init();
 267 }
 268
 269 //________________________________________________________________________
 270 AliPerformancePtCalibMC::~AliPerformancePtCalibMC() {
 271    //
 272    // destructor
 273    //
 274
 275    if(fList) delete fList;
 276    // histograms
 277    if(fHistInvPtPtThetaPhi)       delete fHistInvPtPtThetaPhi;fHistInvPtPtThetaPhi=0;
 278    if(fHistPtShift0)              delete fHistPtShift0;fHistPtShift0=0;
 279    if(fHistPrimaryVertexPosX)     delete fHistPrimaryVertexPosX;fHistPrimaryVertexPosX=0;
 280    if(fHistPrimaryVertexPosY)     delete fHistPrimaryVertexPosY;fHistPrimaryVertexPosY=0;
 281    if(fHistPrimaryVertexPosZ)     delete fHistPrimaryVertexPosZ;fHistPrimaryVertexPosZ=0;
 282    if(fHistTrackMultiplicity)     delete fHistTrackMultiplicity;fHistTrackMultiplicity=0;
 283    if(fHistTrackMultiplicityCuts) delete fHistTrackMultiplicityCuts;fHistTrackMultiplicityCuts=0;
 284
 285    if(fHistTPCMomentaPosP)    delete fHistTPCMomentaPosP;fHistTPCMomentaPosP=0;
 286    if(fHistTPCMomentaNegP)    delete fHistTPCMomentaNegP;fHistTPCMomentaNegP=0;
 287    if(fHistTPCMomentaPosPt)   delete fHistTPCMomentaPosPt;fHistTPCMomentaPosPt=0;
 288    if(fHistTPCMomentaNegPt)   delete fHistTPCMomentaNegPt ;fHistTPCMomentaNegPt=0;
 289    if(fHistUserPtShift)       delete fHistUserPtShift;fHistUserPtShift=0;
 290    if(fHistInvPtPtThetaPhiMC) delete fHistInvPtPtThetaPhiMC;fHistInvPtPtThetaPhiMC=0;
 291    if(fHistInvPtMCESD)  delete fHistInvPtMCESD;fHistInvPtMCESD=0;
 292    if(fHistInvPtMCTPC)  delete fHistInvPtMCTPC;fHistInvPtMCTPC=0;
 293    if(fHistPtMCESD)     delete fHistPtMCESD;fHistPtMCESD=0;
 294    if(fHistPtMCTPC)     delete fHistPtMCTPC;fHistPtMCTPC=0;
 295    if(fHistMomresMCESD) delete fHistMomresMCESD;fHistMomresMCESD=0;
 296    if(fHistMomresMCTPC) delete fHistMomresMCTPC;fHistMomresMCTPC=0;
 297    if(fHistTPCMomentaPosInvPtMC) delete fHistTPCMomentaPosInvPtMC;fHistTPCMomentaPosInvPtMC=0;
 298    if(fHistTPCMomentaNegInvPtMC) delete fHistTPCMomentaNegInvPtMC;fHistTPCMomentaNegInvPtMC=0;
 299    if(fHistTPCMomentaPosPtMC)    delete fHistTPCMomentaPosPtMC;fHistTPCMomentaPosPtMC=0;
 300    if(fHistTPCMomentaNegPtMC)    delete fHistTPCMomentaNegPtMC;fHistTPCMomentaNegPtMC=0;
 301    if(fHistESDMomentaPosInvPtMC) delete fHistESDMomentaPosInvPtMC;fHistESDMomentaPosInvPtMC=0;
 302    if(fHistESDMomentaNegInvPtMC) delete fHistESDMomentaNegInvPtMC;fHistESDMomentaNegInvPtMC=0;
 303    if(fHistESDMomentaPosPtMC)    delete fHistESDMomentaPosPtMC;fHistESDMomentaPosPtMC=0;
 304    if(fHistESDMomentaNegPtMC)    delete fHistESDMomentaNegPtMC;fHistESDMomentaNegPtMC=0;
 305    if(fHistdedxPions)    delete fHistdedxPions;fHistdedxPions=0;
 306
 307
 308    //esd track cuts
 309    if(fESDTrackCuts)   delete fESDTrackCuts;
 310    if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 311
 312
 313
 314 }
 315
 316 //________________________________________________________________________
 317 void AliPerformancePtCalibMC::Init()
 318 {
 319    // Create histograms
 320    // Called once
 321
 322    fList = new TList();
 323
 324    // init folder
 325    fAnalysisFolder = CreateFolder("folderPt_TPC","Analysis Pt Resolution Folder");
 326
 327    // Primary Vertex:
 328    fHistPrimaryVertexPosX       = new TH1F("fHistPrimaryVertexPosX", "Primary Vertex Position X;Primary Vertex Position X (cm);Events",100,-0.5,0.5);
 329    fList->Add(fHistPrimaryVertexPosX);
 330    fHistPrimaryVertexPosY       = new TH1F("fHistPrimaryVertexPosY", "Primary Vertex Position Y;Primary Vertex Position Y (cm);Events",100,-0.5,0.5);
 331    fList->Add(fHistPrimaryVertexPosY);
 332    fHistPrimaryVertexPosZ       = new TH1F("fHistPrimaryVertexPosZ", "Primary Vertex Position Z;Primary Vertex Position Z (cm);Events",200,-2.0,2.0);
 333    fList->Add(fHistPrimaryVertexPosZ);
 334
 335    // Multiplicity:
 336    fHistTrackMultiplicity     = new TH1F("fHistTrackMultiplicity", "Multiplicity distribution;Number of tracks;Events", 250, 0, 250);
 337    fList->Add(fHistTrackMultiplicity);
 338    fHistTrackMultiplicityCuts = new TH1F("fHistTrackMultiplicityCuts", "Multiplicity distribution;Number of tracks after cuts;Events", 250, 0, 250);
 339    fList->Add(fHistTrackMultiplicityCuts);
 340
 341    // momentum histos
 342    fHistPtShift0   = new TH1F("fHistPtShift0","1/pt dN/pt vs. pt of ESD track  ",800,-20.0,20.0);
 343    fList->Add(fHistPtShift0);
 344    const   Int_t invPtDims = 4;
 345    fMaxPhi=6.52;
 346    fMinPhi=0.0;
 347    fMaxTheta=3.0;
 348    fMinTheta=0.0;
 349    Double_t xminInvPt[invPtDims] = {-4.5,-20.0,fMinTheta,fMinPhi};
 350    Double_t xmaxInvPt[invPtDims] = {4.5,20.0,fMaxTheta,fMaxPhi};
 351    Int_t  binsInvPt[invPtDims] = {450,400,150,163};
 352
 353    fHistInvPtPtThetaPhi = new THnSparseF("fHistInvPtPtThetaPhi","1/pt vs pt vs #theta vs #phi ",invPtDims,binsInvPt,xminInvPt,xmaxInvPt);
 354    fList->Add(fHistInvPtPtThetaPhi);
 355
 356    // momentum test histos
 357    fHistTPCMomentaPosP  =  new TH2F("fHistTPCMomentaPosP","TPC p vs global esd track p pos",300,0.0,15.0,300,0.0,15.0);
 358    fList->Add(fHistTPCMomentaPosP);
 359    fHistTPCMomentaNegP  =  new TH2F("fHistTPCMomentaNegP","TPC p vs global esd track p neg",300,0.0,15.0,300,0.0,15.0);
 360    fList->Add(fHistTPCMomentaNegP);
 361    fHistTPCMomentaPosPt =  new TH2F("fHistTPCMomentaPosPt","TPC pt vs global esd track pt pos",300,0.0,15.0,300,0.0,15.0);
 362    fList->Add(fHistTPCMomentaPosPt);
 363    fHistTPCMomentaNegPt =  new TH2F("fHistTPCMomentaNegPt","TPC pt vs global esd track pt neg",300,0.0,15.0,300,0.0,15.0);
 364    fList->Add(fHistTPCMomentaNegPt);
 365
 366    // momentum test histos MC
 367    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);
 368    fList->Add(fHistTPCMomentaPosInvPtMC);
 369    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);
 370    fList->Add(fHistTPCMomentaNegInvPtMC);
 371    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);
 372    fList->Add(fHistTPCMomentaPosPtMC);
 373    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);
 374    fList->Add(fHistTPCMomentaNegPtMC);
 375    fHistESDMomentaPosInvPtMC = new TH1F("fHistESDMomentaPosInvPtMC","ESD-MC of 1/pt pos ",200, -2.0, 2.0);
 376    fList->Add(fHistESDMomentaPosInvPtMC);
 377    fHistESDMomentaNegInvPtMC = new TH1F("fHistESDMomentaNegInvPtMC","ESD-MC of 1/pt neg",200, -2.0, 2.0);
 378    fList->Add(fHistESDMomentaNegInvPtMC);
 379    fHistESDMomentaPosPtMC    = new TH1F("fHistESDMomentaPosPtMC","(ESD-MC)/MC^2 of pt pos",200,-2.0,2.0);
 380    fList->Add(fHistESDMomentaPosPtMC);
 381    fHistESDMomentaNegPtMC    = new TH1F("fHistESDMomentaNegPtMC","(ESD-MC)/MC^2 of pt neg",200,-2.0,2.0);
 382    fList->Add(fHistESDMomentaNegPtMC);
 383
 384    // MC only info
 385    fHistInvPtPtThetaPhiMC = new THnSparseF("fHistInvPtPtThetaPhiMC","MC 1/pt vs pt vs #theta vs #phi ",invPtDims,binsInvPt,xminInvPt,xmaxInvPt);
 386    fList->Add(fHistInvPtPtThetaPhiMC);
 387
 388
 389    //correlation histos MC ESD or TPC
 390    fHistInvPtMCESD  = new TH2F("fHistInvPtMCESD","inv pt ESD vs MC",450, -4.5, 4.5,450, -4.5, 4.5);
 391    fList->Add(fHistInvPtMCESD);
 392    fHistPtMCESD     = new TH2F("fHistPtMCESD"," pt ESD vs MC",500, 0.0, 50.0,500, 0.0, 50.0);
 393    fList->Add(fHistPtMCESD);
 394    fHistInvPtMCTPC  = new TH2F("fHistInvPtMCTPC","inv pt TPC vs MC",450, -4.5, 4.5,450, -4.5, 4.5);
 395    fList->Add(fHistInvPtMCTPC);
 396    fHistPtMCTPC     = new TH2F("fHistPtMCTPC"," pt TPC vs MC",500, 0.0, 50.0,500, 0.0,50.0);
 397    fList->Add(fHistPtMCTPC);
 398    fHistMomresMCESD = new TH2F("fHistMomresMCESD"," (pt ESD - pt MC)/ptMC vs pt MC",500, 0.0, 50.0,200, -2.0, 2.0);
 399    fList->Add(fHistMomresMCESD);
 400    fHistMomresMCTPC = new TH2F("fHistMomresMCTPC"," (pt TPC - pt MC)/ptMC vs pt MC",500, 0.0, 50.0,200, -2.0, 2.0);
 401    fList->Add(fHistMomresMCTPC);
 402
 403
 404    //user pt shift check
 405    fHistUserPtShift = new TH1F("fHistUserPtShift","user defined shift in 1/pt",100,-0.5,1.5);
 406    fList->Add(fHistUserPtShift);
 407    // pid
 408    fHistdedxPions = new TH2F ("fHistdedxPions","dEdx of pions ident. via PDG code vs signed Pt",300,-15.05,15.05,200,0.0,400.0);
 409    fList->Add(fHistdedxPions);
 410
 411    // esd track cuts
 412    fESDTrackCuts =NULL;
 413
 414
 415 }
 416
 417 //________________________________________________________________________
 418 void AliPerformancePtCalibMC::SetPtShift(const Double_t shiftVal ) {
 419
 420    //set user defined shift in charge/pt
 421
 422    if(shiftVal) { fShift=kTRUE; fDeltaInvP = shiftVal; }
 423 }
 424
 425 //________________________________________________________________________
 426 void AliPerformancePtCalibMC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const /*esdFriend*/, const Bool_t /*bUseMC*/, const Bool_t /*bUseESDfriend*/)
 427 {
 428    //exec: read MC and esd or tpc tracks
 429
 430    AliStack* stack = NULL;
 431
 432    if (!esdEvent) {
 433       Printf("ERROR: Event not available");
 434       return;
 435    }
 436
 437
 438    fHistTrackMultiplicity->Fill(esdEvent->GetNumberOfTracks());
 439
 440    if (!mcEvent) {
 441       Printf("ERROR: Could not retrieve MC event");
 442       return;
 443    }
 444    stack = mcEvent->Stack();
 445    if (!stack) {
 446       Printf("ERROR: Could not retrieve stack");
 447       return;
 448    }
 449
 450
 451    //vertex info for cut
 452    //const AliESDVertex *vtx = esdEvent->GetPrimaryVertex();
 453    //if (!vtx->GetStatus()) return ;
 454
 455    if(fShift) fHistUserPtShift->Fill(fDeltaInvP);
 456
 457    // read primary vertex info
 458    Double_t tPrimaryVtxPosition[3];
 459    // Double_t tPrimaryVtxCov[3];
 460    const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexTPC();
 461
 462    tPrimaryVtxPosition[0] = primaryVtx->GetXv();
 463    tPrimaryVtxPosition[1] = primaryVtx->GetYv();
 464    tPrimaryVtxPosition[2] = primaryVtx->GetZv();
 465
 466    fHistPrimaryVertexPosX->Fill(tPrimaryVtxPosition[0]);
 467    fHistPrimaryVertexPosY->Fill(tPrimaryVtxPosition[1]);
 468    fHistPrimaryVertexPosZ->Fill(tPrimaryVtxPosition[2]);
 469
 470
 471    //fill histos for pt spectra and shift of transverse momentum
 472    Int_t count=0;
 473
 474    for(Int_t j = 0;j<esdEvent->GetNumberOfTracks();j++){
 475       AliESDtrack *esdTrack = esdEvent->GetTrack(j);
 476       if(!esdTrack) continue;
 477
 478       //esd track cuts
 479       if(fESDcuts){
 480          if(!fESDTrackCuts->AcceptTrack(esdTrack)) continue;
 481       }
 482
 483       // get MC info
 484       Int_t esdLabel = esdTrack->GetLabel();
 485       if(esdLabel<0) continue;
 486       TParticle *  partMC = stack->Particle(esdLabel);
 487       if (!partMC) continue;
 488
 489       // fill correlation histos MC ESD
 490       Double_t pESD  = esdTrack->GetP();
 491       Double_t ptESD = esdTrack->GetSignedPt();
 492
 493       if(!ptESD || !(partMC->Pt()) ) continue;
 494       Double_t mcPt = partMC->Pt();
 495       Double_t invPtMC = 1.0/mcPt;
 496       Int_t signMC = partMC->GetPdgCode();
 497
 498       //pid
 499       if(fPions && !(fabs(signMC)-211<0.001)) continue;// pions only
 500
 501       //MC only
 502       if(signMC>0) signMC = 1;
 503       else signMC = -1;
 504
 505       //fill MC information
 506       Double_t thetaMC = partMC->Theta();
 507       Double_t phiMC = partMC->Phi();
 508
 509       Double_t momAngMC[4] = {signMC*(fabs(invPtMC)),signMC*(fabs(mcPt)),thetaMC,phiMC};
 510
 511       // fill only if MC track is in eta acceptance of TPC in order to be compareable to TPC tracks
 512       if(fabs( partMC->Eta())< fEtaAcceptance) {
 513          fHistInvPtPtThetaPhiMC->Fill(momAngMC);
 514
 515          //correlation histos MC ESD
 516          fHistInvPtMCESD->Fill(signMC*(fabs(invPtMC)),1.0/ptESD);
 517          fHistPtMCESD->Fill(fabs(mcPt),fabs(ptESD));
 518       }
 519
 520       // fill histos TPC or ESD
 521       if(fOptTPC){
 522          //TPC tracks and MC tracks
 523          const AliExternalTrackParam *tpcTrack = esdTrack->GetTPCInnerParam();
 524          if(!tpcTrack) continue;
 525          if(fabs(tpcTrack->Eta())>=  fEtaAcceptance) continue;
 526
 527          Double_t signedPt = tpcTrack->GetSignedPt();
 528          Double_t invPt = 0.0;
 529          if(signedPt) {
 530             invPt = 1.0/signedPt;
 531
 532             fHistPtShift0->Fill(signedPt);//changed
 533
 534             if(fShift){Printf("user shift of momentum SET to non zero value!");
 535                invPt += fDeltaInvP; //shift momentum for tests
 536                if(invPt) signedPt = 1.0/invPt;
 537                else continue;
 538             }
 539
 540
 541             Double_t theta = tpcTrack->Theta();
 542             Double_t phi = tpcTrack->Phi();
 543
 544             Double_t momAng[4] = {invPt,signedPt,theta,phi};
 545             fHistInvPtPtThetaPhi->Fill(momAng);
 546
 547             //correlation histos MC TPC
 548             fHistInvPtMCTPC->Fill(signMC*(fabs(invPtMC)),invPt);
 549             fHistPtMCTPC->Fill(fabs(mcPt),fabs(signedPt));
 550
 551             //compare to MC info
 552             Double_t  ptDiffESD = (fabs(ptESD)-fabs(mcPt))/pow(mcPt,2);
 553             Double_t  ptDiffTPC = (fabs(signedPt)-fabs(mcPt))/pow(mcPt,2);
 554             Double_t  invPtDiffESD = fabs(1.0/ptESD)-1.0/fabs(mcPt);
 555             Double_t  invPtDiffTPC = fabs(invPt)-1.0/fabs(mcPt);
 556             Double_t  pTPC  = tpcTrack->GetP();
 557
 558             if(esdTrack->GetSign()>0){//compare momenta ESD track and TPC track
 559                fHistTPCMomentaPosP->Fill(fabs(pESD),fabs(pTPC));
 560                fHistTPCMomentaPosPt->Fill(fabs(ptESD),fabs(signedPt));
 561                fHistTPCMomentaPosInvPtMC->Fill(invPtDiffESD,invPtDiffTPC);
 562                fHistTPCMomentaPosPtMC->Fill(ptDiffESD,ptDiffTPC);
 563             }
 564             else{
 565                fHistTPCMomentaNegP->Fill(fabs(pESD),fabs(pTPC));
 566                fHistTPCMomentaNegPt->Fill(fabs(ptESD),fabs(signedPt));
 567                fHistTPCMomentaNegInvPtMC->Fill(invPtDiffESD,invPtDiffTPC);
 568                fHistTPCMomentaNegPtMC->Fill(ptDiffESD,ptDiffTPC);
 569             }
 570             fHistdedxPions->Fill(signedPt,esdTrack->GetTPCsignal());
 571             fHistMomresMCESD->Fill(fabs(mcPt),(fabs(ptESD)-fabs(mcPt))/fabs(mcPt));
 572             fHistMomresMCTPC->Fill(fabs(mcPt),(fabs(signedPt)-fabs(mcPt))/fabs(mcPt));
 573             count++;
 574          }
 575          else continue;
 576       }
 577
 578  else{
 579     // ESD tracks and MC tracks
 580     if(fabs(esdTrack->Eta())>= fEtaAcceptance) continue;
 581     Double_t invPt = 0.0;
 582
 583     if(ptESD) {
 584        invPt = 1.0/ptESD;
 585        fHistPtShift0->Fill(ptESD);
 586
 587        if(fShift){Printf("user shift of momentum SET to non zero value!");
 588           invPt += fDeltaInvP; //shift momentum for tests
 589           if(invPt) ptESD = 1.0/invPt;
 590           else continue;
 591        }
 592
 593        Double_t theta = esdTrack->Theta();
 594        Double_t phi = esdTrack->Phi();
 595
 596        Double_t momAng[4] = {invPt,ptESD,theta,phi};
 597        fHistInvPtPtThetaPhi->Fill(momAng);
 598
 599        //differences MC ESD tracks
 600        Double_t ptDiffESD = (fabs(ptESD)-fabs(mcPt))/pow(mcPt,2);
 601        Double_t invPtdiffESD = fabs(1.0/ptESD)-1.0/fabs(mcPt);
 602        if(esdTrack->GetSign()>0){
 603           fHistESDMomentaPosInvPtMC->Fill(invPtdiffESD);
 604           fHistESDMomentaPosPtMC->Fill(ptDiffESD);
 605        }
 606        else{
 607           fHistESDMomentaNegInvPtMC->Fill(invPtdiffESD);
 608           fHistESDMomentaNegPtMC->Fill(ptDiffESD);
 609        }
 610        fHistdedxPions->Fill(ptESD,esdTrack->GetTPCsignal());
 611        fHistMomresMCESD->Fill(fabs(mcPt),(fabs(ptESD)-fabs(mcPt))/fabs(mcPt));
 612        count++;
 613     }
 614  }
 615 }
 616
 617 fHistTrackMultiplicityCuts->Fill(count);
 618
 619 }
 620
 621 //______________________________________________________________________________________________________________________
 622
 623 void AliPerformancePtCalibMC::Analyse()
 624 {
 625
 626    // analyse charge/pt spectra in bins of theta and phi. Bins can be set by user
 627
 628
 629    THnSparseF *copyTHnSparseTheta;
 630    THnSparseF *copyTHnSparsePhi;
 631
 632    if(fAnaMC){
 633       Printf("AliPerformancePtCalibMC::Analyse: analysing MC!");
 634       copyTHnSparseTheta = (THnSparseF*)fHistInvPtPtThetaPhiMC->Clone("copyTHnSparseThetaMC");
 635       copyTHnSparsePhi   = (THnSparseF*)fHistInvPtPtThetaPhiMC->Clone("copyTHnSparsePhiMC");
 636    }
 637    else {
 638       Printf("AliPerformancePtCalibMC::Analyse: analysing ESD (reco)!");
 639       copyTHnSparseTheta = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparseTheta");
 640       copyTHnSparsePhi   = (THnSparseF*)fHistInvPtPtThetaPhi->Clone("copyTHnSparsePhi");
 641    }
 642
 643    copyTHnSparseTheta->GetAxis(3)->SetRangeUser(fMinPhi,fMaxPhi);
 644    copyTHnSparsePhi->GetAxis(2)->SetRangeUser(fMinTheta,fMaxTheta);
 645
 646    TH2F *histInvPtTheta = (TH2F*)copyTHnSparseTheta->Projection(2,0);
 647    TH2F *histInvPtPhi   = (TH2F*)copyTHnSparsePhi->Projection(3,0);
 648
 649    AliPerfAnalyzeInvPt *ana = new  AliPerfAnalyzeInvPt("AliPerfAnalyzeInvPt","AliPerfAnalyzeInvPt");
 650
 651    TH1::AddDirectory(kFALSE);
 652
 653    ana->SetProjBinsTheta(fThetaBins,fNThetaBins);
 654    ana->SetProjBinsPhi(fPhiBins,fNPhiBins);
 655    ana->SetMakeFitOption(fFitGaus,fExclRange,fRange);
 656    if(fDoRebin) ana->SetDoRebin(fRebin);
 657    TObjArray *aFolderObj = new TObjArray;
 658
 659    ana->StartAnalysis(histInvPtTheta,histInvPtPhi, aFolderObj);
 660
 661    // export objects to analysis folder
 662    fAnalysisFolder = ExportToFolder(aFolderObj);
 663
 664    // delete only TObjArray
 665    if(aFolderObj) delete aFolderObj;
 666    if(ana) delete ana;
 667
 668 }
 669
 670 //______________________________________________________________________________________________________________________
 671 TFolder* AliPerformancePtCalibMC::ExportToFolder(TObjArray * array)
 672 {
 673    // recreate folder avery time and export objects to new one
 674    //
 675    AliPerformancePtCalibMC * comp=this;
 676    TFolder *folder = comp->GetAnalysisFolder();
 677
 678    TString name, title;
 679    TFolder *newFolder = 0;
 680    Int_t i = 0;
 681    Int_t size = array->GetSize();
 682
 683    if(folder) {
 684       // get name and title from old folder
 685       name = folder->GetName();
 686       title = folder->GetTitle();
 687
 688       // delete old one
 689       delete folder;
 690
 691       // create new one
 692       newFolder = CreateFolder(name.Data(),title.Data());
 693       newFolder->SetOwner();
 694
 695       // add objects to folder
 696       while(i < size) {
 697          newFolder->Add(array->At(i));
 698          i++;
 699       }
 700    }
 701
 702    return newFolder;
 703 }
 704
 705 //______________________________________________________________________________________________________________________
 706 Long64_t AliPerformancePtCalibMC::Merge(TCollection* const list)
 707 {
 708    // Merge list of objects (needed by PROOF)
 709
 710    if (!list)
 711       return 0;
 712
 713    if (list->IsEmpty())
 714       return 1;
 715
 716    TIterator* iter = list->MakeIterator();
 717    TObject* obj = 0;
 718
 719    // collection of generated histograms
 720    Int_t count=0;
 721    while((obj = iter->Next()) != 0)
 722       {
 723          AliPerformancePtCalibMC* entry = dynamic_cast<AliPerformancePtCalibMC*>(obj);
 724          if (!entry) continue;
 725          fHistInvPtPtThetaPhi->Add(entry->fHistInvPtPtThetaPhi);
 726
 727          fHistInvPtPtThetaPhiMC->Add(entry->fHistInvPtPtThetaPhiMC);
 728
 729          fHistInvPtMCESD->Add(entry->fHistInvPtMCESD);
 730          fHistPtMCESD->Add(entry->fHistPtMCESD);
 731          fHistInvPtMCTPC->Add(entry->fHistInvPtMCTPC);
 732          fHistPtMCTPC->Add(entry->fHistPtMCTPC);
 733          fHistMomresMCESD->Add(entry->fHistMomresMCESD);
 734          fHistMomresMCTPC->Add(entry->fHistMomresMCTPC);
 735
 736          fHistPtShift0->Add(entry->fHistPtShift0);
 737          fHistPrimaryVertexPosX->Add(entry->fHistPrimaryVertexPosX);
 738          fHistPrimaryVertexPosY->Add(entry->fHistPrimaryVertexPosY);
 739          fHistPrimaryVertexPosZ->Add(entry->fHistPrimaryVertexPosZ);
 740          fHistTrackMultiplicity->Add(entry->fHistTrackMultiplicity);
 741          fHistTrackMultiplicityCuts->Add(entry->fHistTrackMultiplicityCuts);
 742
 743          fHistTPCMomentaPosP->Add(entry->fHistTPCMomentaPosP);
 744          fHistTPCMomentaNegP->Add(entry->fHistTPCMomentaNegP);
 745          fHistTPCMomentaPosPt->Add(entry->fHistTPCMomentaPosPt);
 746          fHistTPCMomentaNegPt->Add(entry->fHistTPCMomentaNegPt);
 747          fHistTPCMomentaPosInvPtMC->Add(entry->fHistTPCMomentaPosInvPtMC);
 748          fHistTPCMomentaNegInvPtMC->Add(entry->fHistTPCMomentaNegInvPtMC);
 749          fHistTPCMomentaPosPtMC->Add(entry->fHistTPCMomentaPosPtMC);
 750          fHistTPCMomentaNegPtMC->Add(entry->fHistTPCMomentaNegPtMC);
 751          fHistESDMomentaPosInvPtMC->Add(entry->fHistESDMomentaPosInvPtMC);
 752          fHistESDMomentaNegInvPtMC->Add(entry->fHistESDMomentaNegInvPtMC);
 753          fHistESDMomentaPosPtMC->Add(entry->fHistESDMomentaPosPtMC);
 754          fHistESDMomentaNegPtMC->Add(entry->fHistESDMomentaNegPtMC);
 755          fHistdedxPions->Add(entry->fHistdedxPions);
 756          count++;
 757       }
 758
 759    return count;
 760 }
 761
 762 //______________________________________________________________________________________________________________________
 763 TFolder* AliPerformancePtCalibMC::CreateFolder(TString name,TString title) {
 764    // create folder for analysed histograms
 765    //
 766    TFolder *folder = 0;
 767    folder = new TFolder(name.Data(),title.Data());
 768
 769    return folder;
 770 }
 771
 772
 773 // set variables for Analyse()
 774 //______________________________________________________________________________________________________________________
 775 void AliPerformancePtCalibMC::SetProjBinsPhi(const Double_t *phiBinArray,const Int_t nphBins, const  Double_t minTheta, const  Double_t maxTheta){
 776    // set phi bins for Analyse()
 777    // set phi bins as array and set number of this array which is equal to number of bins analysed
 778    // the last analysed bin will always be the projection from first to last bin in the array
 779    if(nphBins){
 780       fNPhiBins = nphBins;
 781
 782       for(Int_t k = 0;k<fNPhiBins;k++){
 783          fPhiBins[k] = phiBinArray[k];
 784       }
 785       Printf("AliPerformancePtCalibMC::SetProjBinsPhi: number of bins in phi set to %i",fNPhiBins);
 786    }
 787    else  Printf("Warning AliPerformancePtCalibMC::SetProjBinsPhi: number of bins in phi NOT set!!! Default values are taken.");
 788
 789    if(fabs(minTheta-maxTheta)<0.001){
 790       Printf("AliPerformancePtCalibMC::SetProjBinsPhi: theta range for projection for projection on phi and charge/pt is too small. whole range of theta selected.");
 791    }
 792    else{
 793       Printf("AliPerformancePtCalibMC::SetProjBinsPhi: theta range for projection on phi and charge/pt is selected by user: %1.3f - %1.3f rad",minTheta,maxTheta);
 794       fMinTheta = minTheta;
 795       fMaxTheta = maxTheta;
 796    }
 797 }
 798 //____________________________________________________________________________________________________________________________________________
 799 void AliPerformancePtCalibMC::SetProjBinsTheta(const Double_t *thetaBinArray, const Int_t nthBins, const Double_t minPhi, const Double_t maxPhi){
 800    // set theta bins for Analyse()
 801    // set theta bins as array and set number of this array which is equal to number of bins analysed
 802    // the last analysed bin will always be the projection from first to last bin in the array
 803    if(nthBins){
 804       fNThetaBins = nthBins;
 805       for(Int_t k = 0;k<fNThetaBins;k++){
 806          fThetaBins[k] = thetaBinArray[k];
 807       }
 808       Printf("AliPerformancePtCalibMC::SetProjBinsTheta: number of bins in theta set to %i",fNThetaBins);
 809    }
 810    else  Printf("Warning AliPerformancePtCalibMC::SetProjBinsTheta: number of bins in theta NOT set!!! Default values are taken.");
 811
 812    if(fabs(minPhi-maxPhi)<0.001){
 813       Printf("AliPerformancePtCalibMC::SetProjBinsTheta: phi range for projection for projection on theta and charge/pt is too small. whole range of phi selected.");
 814    }
 815    else{
 816       Printf("AliPerformancePtCalibMC::SetProjBinsTheta: phi range for projection on theta and charge/pt is selected by user: %1.3f - %1.3f rad",minPhi,maxPhi);
 817       fMinPhi = minPhi;
 818       fMaxPhi = maxPhi;
 819    }
 820 }
 821
 822 //____________________________________________________________________________________________________________________________________________
 823 void AliPerformancePtCalibMC::SetMakeFitOption(const Bool_t setGausFit, const Double_t exclusionR,const Double_t fitR ){
 824
 825    // set the fit options:
 826    // for usage of gaussian function instead of polynomial (default) set setGausFit=kTRUE
 827    // set the range of rejection of points around 0 via exclusionR
 828    // set the fit range around 0 with fitR
 829
 830    fFitGaus = setGausFit;
 831    fExclRange  = exclusionR;
 832    fRange = fitR;
 833
 834    if(fFitGaus) Printf("AliPerformancePtCalibMC::SetMakeFitOption: set MakeGausFit with fit range %2.3f and exclusion range in fabs(1/pt): %2.3f",fRange,fExclRange);
 835    else  Printf("AliPerformancePtCalibMC::SetMakeFitOption: set standard polynomial fit with fit range %2.3f and exclusion range in fabs(1/pt): %2.3f",fRange,fExclRange);
 836
 837 }