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