PWG1/TPC/AliPerformancePtCalibMC.cxx

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