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    for (Int_t i=0; i<100; i++){
 173      fThetaBins[i] =  0;
 174      fPhiBins[i] =  0;
 175    }
 176
 177    Init();
 178 }
 179
 180 //________________________________________________________________________
 181 AliPerformancePtCalibMC::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 //________________________________________________________________________
 280 AliPerformancePtCalibMC::~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 //________________________________________________________________________
 327 void 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 //________________________________________________________________________
 428 void 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 //________________________________________________________________________
 436 void 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
 627 fHistTrackMultiplicityCuts->Fill(count);
 628
 629 }
 630
 631 //______________________________________________________________________________________________________________________
 632
 633 void 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 //______________________________________________________________________________________________________________________
 683 TFolder* 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 //______________________________________________________________________________________________________________________
 718 Long64_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 //______________________________________________________________________________________________________________________
 775 TFolder* 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 //______________________________________________________________________________________________________________________
 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
 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 //____________________________________________________________________________________________________________________________________________
 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
 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 //____________________________________________________________________________________________________________________________________________
 835 void 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 }