PWGCF/Correlations/DPhi/TriggerPID/AliTwoParticlePIDCorr.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 #include "AliTwoParticlePIDCorr.h"
  17 #include "AliVParticle.h"
  18 #include "TFormula.h"
  19 #include "TAxis.h"
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TH1F.h"
  23 #include "TH2F.h"
  24 #include "TH3F.h"
  25 #include "TProfile.h"
  26 #include "TList.h"
  27 #include "TFile.h"
  28 #include "TGrid.h"
  29
  30 #include "AliCentrality.h"
  31 #include "Riostream.h"
  32
  33 #include "AliTHn.h"
  34 #include "AliCFContainer.h"
  35 #include "THn.h"
  36 #include "THnSparse.h"
  37
  38 #include <TSpline.h>
  39 #include <AliPID.h>
  40 #include "AliESDpid.h"
  41 #include "AliAODpidUtil.h"
  42 #include <AliPIDResponse.h>
  43 #include "AliPIDCombined.h"
  44
  45 #include <AliAnalysisManager.h>
  46 #include <AliInputEventHandler.h>
  47 #include "AliAODInputHandler.h"
  48
  49 #include "AliAnalysisTaskSE.h"
  50 #include "AliAnalysisManager.h"
  51 #include "AliCentrality.h"
  52
  53 #include "AliVEvent.h"
  54 #include "AliAODEvent.h"
  55 #include "AliAODTrack.h"
  56 #include "AliVTrack.h"
  57
  58 #include "THnSparse.h"
  59
  60 #include "AliAODMCHeader.h"
  61 #include "AliAODMCParticle.h"
  62 #include "AliMCEventHandler.h"
  63 #include "AliMCEvent.h"
  64 #include "AliMCParticle.h"
  65 #include "TParticle.h"
  66 #include <TDatabasePDG.h>
  67 #include <TParticlePDG.h>
  68
  69 #include "AliGenCocktailEventHeader.h"
  70 #include "AliGenEventHeader.h"
  71 #include "AliCollisionGeometry.h"
  72 #include "AliOADBContainer.h"
  73
  74 #include "AliEventPoolManager.h"
  75 #include "AliAnalysisUtils.h"
  76 using namespace AliPIDNameSpace;
  77 using namespace std;
  78
  79 ClassImp(AliTwoParticlePIDCorr)
  80 ClassImp(LRCParticlePID)
  81
  82 const char * kPIDTypeName[]={"TPC","TOF","TPC-TOF"} ;
  83 const char * kDetectorName[]={"ITS","TPC","TOF"} ;
  84 const char * kParticleSpeciesName[]={"Pions","Kaons","Protons","Undefined"} ;
  85 //Source code::dphicorrelations,VnV0, TaskBFpsi, AliHelperPID,
  86
  87 //________________________________________________________________________
  88 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
  89 :AliAnalysisTaskSE(),
  90   fOutput(0),
  91    fOutputList(0),
  92   fList(0),
  93   fCentralityMethod("V0A"),
  94   fSampleType("pPb"),
  95  fRequestEventPlane(kFALSE),
  96   fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default)
  97   trkVtx(0),
  98   zvtx(0),
  99   fFilterBit(768),
 100   fTrackStatus(0),
 101   fSharedClusterCut(-1),
 102  fSharedTPCmapCut(-1),
 103   fVertextype(1),
 104  skipParticlesAbove(0),
 105   fzvtxcut(10.0),
 106   ffilltrigassoUNID(kFALSE),
 107   ffilltrigUNIDassoID(kFALSE),
 108   ffilltrigIDassoUNID(kTRUE),
 109   ffilltrigIDassoID(kFALSE),
 110   ffilltrigIDassoIDMCTRUTH(kFALSE),
 111   fMaxNofMixingTracks(50000),
 112   fPtOrderMCTruth(kTRUE),
 113  fPtOrderDataReco(kTRUE),
 114   fWeightPerEvent(kFALSE),
 115   fTriggerSpeciesSelection(kFALSE),
 116   fAssociatedSpeciesSelection(kFALSE),
 117  fRandomizeReactionPlane(kFALSE),
 118   fTriggerSpecies(SpPion),
 119   fAssociatedSpecies(SpPion),
 120   fCustomBinning(""),
 121   fBinningString(""),
 122   fSelectHighestPtTrig(kFALSE),
 123   fcontainPIDtrig(kTRUE),
 124   fcontainPIDasso(kFALSE),
 125   SetChargeAxis(0),
 126   frejectPileUp(kFALSE),
 127   fminPt(0.2),
 128   fmaxPt(20.0),
 129   fmineta(-0.8),
 130   fmaxeta(0.8),
 131   fselectprimaryTruth(kTRUE),
 132   fonlyprimarydatareco(kFALSE),
 133   fdcacut(kFALSE),
 134   fdcacutvalue(3.0),
 135   ffillhistQAReco(kFALSE),
 136   ffillhistQATruth(kFALSE),
 137   kTrackVariablesPair(0),
 138   fminPtTrig(0),
 139   fmaxPtTrig(0),
 140   fminPtComboeff(2.0),
 141   fmaxPtComboeff(4.0),
 142   fminPtAsso(0),
 143   fmaxPtAsso(0),
 144  fmincentmult(0),
 145  fmaxcentmult(0),
 146  fPriHistShare(0),
 147   fhistcentrality(0),
 148   fEventCounter(0),
 149   fEtaSpectrasso(0),
 150   fphiSpectraasso(0),
 151   MCtruthpt(0),
 152   MCtrutheta(0),
 153   MCtruthphi(0),
 154   MCtruthpionpt(0),
 155   MCtruthpioneta(0),
 156   MCtruthpionphi(0),
 157   MCtruthkaonpt(0),
 158   MCtruthkaoneta(0),
 159   MCtruthkaonphi(0),
 160   MCtruthprotonpt(0),
 161   MCtruthprotoneta(0),
 162   MCtruthprotonphi(0),
 163   fPioncont(0),
 164   fKaoncont(0),
 165   fProtoncont(0),
 166   fUNIDcont(0),
 167   fEventno(0),
 168   fEventnobaryon(0),
 169   fEventnomeson(0),
 170  fhistJetTrigestimate(0),
 171 fTwoTrackDistancePtdip(0x0),
 172 fTwoTrackDistancePtdipmix(0x0),
 173   fCentralityCorrelation(0x0),
 174  fHistVZEROAGainEqualizationMap(0),
 175   fHistVZEROCGainEqualizationMap(0),
 176  fHistVZEROChannelGainEqualizationMap(0),
 177 fCentralityWeights(0),
 178  fHistCentStats(0x0),
 179  fHistRefmult(0x0),
 180  fHistEQVZEROvsTPCmultiplicity(0x0),
 181     fHistEQVZEROAvsTPCmultiplicity(0x0),
 182     fHistEQVZEROCvsTPCmultiplicity(0x0),
 183     fHistVZEROCvsEQVZEROCmultiplicity(0x0),
 184     fHistVZEROAvsEQVZEROAmultiplicity(0x0),
 185     fHistVZEROCvsVZEROAmultiplicity(0x0),
 186     fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
 187     fHistVZEROSignal(0x0),
 188 fHistEventPlaneTruth(0x0),
 189 fHistPsiMinusPhi(0x0),
 190 fEventPlanePID(0x0),
 191 evplaneMC(999.),
 192  fgPsi2v0a(999.),
 193     fgPsi2v0c(999.),
 194     fgPsi2tpc(999.),
 195     fgPsi3v0a(999.),
 196     fgPsi3v0c(999.),
 197     fgPsi3tpc(999.),
 198     fgPsi2v0aMC(999.),
 199     fgPsi2v0cMC(999.),
 200     fgPsi2tpcMC(999.),
 201     fgPsi3v0aMC(999.),
 202     fgPsi3v0cMC(999.),
 203     fgPsi3tpcMC(999.),
 204  gReactionPlane(999.),
 205   fV2(kTRUE),
 206  fV3(kFALSE),
 207  fIsAfter2011(kTRUE),
 208   fRun(-1),
 209   fNcluster(70),
 210  fEPdet("V0A"),
 211  fMultV0(NULL),
 212   fV0Cpol(100),
 213   fV0Apol(100),
 214  fHResTPCv0A2(NULL),
 215 fHResTPCv0C2(NULL),
 216 fHResv0Cv0A2(NULL),
 217 fHResTPCv0A3(NULL),
 218 fHResTPCv0C3(NULL),
 219 fHResv0Cv0A3(NULL),
 220  fHResMA2(NULL),
 221 fHResMC2(NULL),
 222 fHResAC2(NULL),
 223 fHResMA3(NULL),
 224 fHResMC3(NULL),
 225 fHResAC3(NULL),
 226 fPhiRPTPC(NULL),
 227 fPhiRPTPCv3(NULL),
 228 fPhiRPv0A(NULL),
 229 fPhiRPv0C(NULL),
 230 fPhiRPv0Av3(NULL),
 231 fPhiRPv0Cv3(NULL),
 232  fControlConvResoncances(0),
 233   fHistoTPCdEdx(0x0),
 234   fHistoTOFbeta(0x0),
 235   fTPCTOFPion3d(0),
 236   fTPCTOFKaon3d(0),
 237   fTPCTOFProton3d(0),
 238   fPionPt(0),
 239   fPionEta(0),
 240   fPionPhi(0),
 241   fKaonPt(0),
 242   fKaonEta(0),
 243   fKaonPhi(0),
 244   fProtonPt(0),
 245   fProtonEta(0),
 246   fProtonPhi(0),
 247   fCorrelatonTruthPrimary(0),
 248   fCorrelatonTruthPrimarymix(0),
 249   fTHnCorrUNID(0),
 250   fTHnCorrUNIDmix(0),
 251   fTHnCorrID(0),
 252   fTHnCorrIDmix(0),
 253   fTHnCorrIDUNID(0),
 254   fTHnCorrIDUNIDmix(0),
 255   fTHnTrigcount(0),
 256   fTHnTrigcountMCTruthPrim(0),
 257   fPoolMgr(0x0),
 258   fArrayMC(0),
 259   fAnalysisType("AOD"),
 260   fefffilename(""),
 261  ftwoTrackEfficiencyCutDataReco(kTRUE),
 262   twoTrackEfficiencyCutValue(0.02),
 263   fPID(NULL),
 264  fPIDCombined(NULL),
 265  eventno(0),
 266   fPtTOFPIDmin(0.5),
 267   fPtTOFPIDmax(4.0),
 268   fRequestTOFPID(kTRUE),
 269   fPIDType(NSigmaTPCTOF),
 270  fFIllPIDQAHistos(kTRUE),
 271   fNSigmaPID(3),
 272   fBayesCut(0.8),
 273  fdiffPIDcutvalues(kFALSE),
 274  fPIDCutval1(0.0),
 275  fPIDCutval2(0.0),
 276  fPIDCutval3(0.0),
 277  fPIDCutval4(0.0),
 278  fHighPtKaonNSigmaPID(-1),
 279  fHighPtKaonSigma(3.5),
 280   fUseExclusiveNSigma(kFALSE),
 281   fRemoveTracksT0Fill(kFALSE),
 282 fSelectCharge(0),
 283 fTriggerSelectCharge(0),
 284 fAssociatedSelectCharge(0),
 285 fTriggerRestrictEta(-1),
 286 fEtaOrdering(kFALSE),
 287 fCutConversions(kFALSE),
 288 fCutResonances(kFALSE),
 289 fRejectResonanceDaughters(-1),
 290   fOnlyOneEtaSide(0),
 291 fInjectedSignals(kFALSE),
 292   fRemoveWeakDecays(kFALSE),
 293 fRemoveDuplicates(kFALSE),
 294   fapplyTrigefficiency(kFALSE),
 295   fapplyAssoefficiency(kFALSE),
 296   ffillefficiency(kFALSE),
 297   fmesoneffrequired(kFALSE),
 298   fkaonprotoneffrequired(kFALSE),
 299 fAnalysisUtils(0x0),
 300   fDCAXYCut(0)
 301
 302 {
 303  for ( Int_t i = 0; i < 16; i++) {
 304     fHistQA[i] = NULL;
 305   }
 306
 307  for ( Int_t i = 0; i < 6; i++ ){
 308     fTrackHistEfficiency[i] = NULL;
 309     effcorection[i]=NULL;
 310     //effmap[i]=NULL;
 311   }
 312  for ( Int_t i = 0; i < 2; i++ ){
 313    fTwoTrackDistancePt[i]=NULL;
 314    fTwoTrackDistancePtmix[i]=NULL;
 315 }
 316
 317  for(Int_t ipart=0;ipart<NSpecies;ipart++)
 318     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
 319       fnsigmas[ipart][ipid]=999.;
 320
 321  for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
 322
 323   for(Int_t i = 0; i != 2; ++i)
 324     for(Int_t j = 0; j != 2; ++j)
 325       for(Int_t iC = 0; iC < 9; iC++){
 326         fMeanQ[iC][i][j] = 0.;
 327         fWidthQ[iC][i][j] = 1.;
 328         fMeanQv3[iC][i][j] = 0.;
 329         fWidthQv3[iC][i][j] = 1.;
 330     }
 331
 332   }
 333 //________________________________________________________________________
 334 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
 335   :AliAnalysisTaskSE(name),
 336  fOutput(0),
 337    fOutputList(0),
 338    fList(0),
 339  fCentralityMethod("V0A"),
 340   fSampleType("pPb"),
 341  fRequestEventPlane(kFALSE),
 342   fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default)
 343   trkVtx(0),
 344   zvtx(0),
 345   fFilterBit(768),
 346   fTrackStatus(0),
 347   fSharedClusterCut(-1),
 348   fSharedTPCmapCut(-1),
 349   fVertextype(1),
 350    skipParticlesAbove(0),
 351   fzvtxcut(10.0),
 352   ffilltrigassoUNID(kFALSE),
 353   ffilltrigUNIDassoID(kFALSE),
 354   ffilltrigIDassoUNID(kTRUE),
 355   ffilltrigIDassoID(kFALSE),
 356   ffilltrigIDassoIDMCTRUTH(kFALSE),
 357   fMaxNofMixingTracks(50000),
 358   fPtOrderMCTruth(kTRUE),
 359   fPtOrderDataReco(kTRUE),
 360   fWeightPerEvent(kFALSE),
 361   fTriggerSpeciesSelection(kFALSE),
 362   fAssociatedSpeciesSelection(kFALSE),
 363    fRandomizeReactionPlane(kFALSE),
 364   fTriggerSpecies(SpPion),
 365   fAssociatedSpecies(SpPion),
 366   fCustomBinning(""),
 367   fBinningString(""),
 368   fSelectHighestPtTrig(kFALSE),
 369   fcontainPIDtrig(kTRUE),
 370   fcontainPIDasso(kFALSE),
 371   SetChargeAxis(0),
 372   frejectPileUp(kFALSE),
 373   fminPt(0.2),
 374   fmaxPt(20.0),
 375   fmineta(-0.8),
 376   fmaxeta(0.8),
 377   fselectprimaryTruth(kTRUE),
 378   fonlyprimarydatareco(kFALSE),
 379    fdcacut(kFALSE),
 380   fdcacutvalue(3.0),
 381   ffillhistQAReco(kFALSE),
 382   ffillhistQATruth(kFALSE),
 383  kTrackVariablesPair(0),
 384   fminPtTrig(0),
 385   fmaxPtTrig(0),
 386   fminPtComboeff(2.0),
 387   fmaxPtComboeff(4.0),
 388   fminPtAsso(0),
 389   fmaxPtAsso(0),
 390    fmincentmult(0),
 391    fmaxcentmult(0),
 392    fPriHistShare(0),
 393   fhistcentrality(0),
 394   fEventCounter(0),
 395   fEtaSpectrasso(0),
 396   fphiSpectraasso(0),
 397   MCtruthpt(0),
 398   MCtrutheta(0),
 399   MCtruthphi(0),
 400   MCtruthpionpt(0),
 401   MCtruthpioneta(0),
 402   MCtruthpionphi(0),
 403   MCtruthkaonpt(0),
 404   MCtruthkaoneta(0),
 405   MCtruthkaonphi(0),
 406   MCtruthprotonpt(0),
 407   MCtruthprotoneta(0),
 408   MCtruthprotonphi(0),
 409   fPioncont(0),
 410   fKaoncont(0),
 411   fProtoncont(0),
 412    fUNIDcont(0),
 413   fEventno(0),
 414   fEventnobaryon(0),
 415   fEventnomeson(0),
 416   fhistJetTrigestimate(0),
 417 fTwoTrackDistancePtdip(0x0),
 418 fTwoTrackDistancePtdipmix(0x0),
 419   fCentralityCorrelation(0x0),
 420  fHistVZEROAGainEqualizationMap(0),
 421   fHistVZEROCGainEqualizationMap(0),
 422    fHistVZEROChannelGainEqualizationMap(0),
 423 fCentralityWeights(0),
 424   fHistCentStats(0x0),
 425   fHistRefmult(0x0),
 426     fHistEQVZEROvsTPCmultiplicity(0x0),
 427     fHistEQVZEROAvsTPCmultiplicity(0x0),
 428     fHistEQVZEROCvsTPCmultiplicity(0x0),
 429     fHistVZEROCvsEQVZEROCmultiplicity(0x0),
 430     fHistVZEROAvsEQVZEROAmultiplicity(0x0),
 431     fHistVZEROCvsVZEROAmultiplicity(0x0),
 432     fHistEQVZEROCvsEQVZEROAmultiplicity(0x0),
 433     fHistVZEROSignal(0x0),
 434 fHistEventPlaneTruth(0x0),
 435    fHistPsiMinusPhi(0x0),
 436 fEventPlanePID(0x0),
 437 evplaneMC(999.),
 438  fgPsi2v0a(999.),
 439     fgPsi2v0c(999.),
 440     fgPsi2tpc(999.),
 441     fgPsi3v0a(999.),
 442     fgPsi3v0c(999.),
 443     fgPsi3tpc(999.),
 444     fgPsi2v0aMC(999.),
 445     fgPsi2v0cMC(999.),
 446     fgPsi2tpcMC(999.),
 447     fgPsi3v0aMC(999.),
 448     fgPsi3v0cMC(999.),
 449     fgPsi3tpcMC(999.),
 450    gReactionPlane(999.),
 451  fV2(kTRUE),
 452  fV3(kFALSE),
 453  fIsAfter2011(kTRUE),
 454   fRun(-1),
 455   fNcluster(70),
 456    fEPdet("V0A"),
 457  fMultV0(NULL),
 458   fV0Cpol(100),
 459   fV0Apol(100),
 460  fHResTPCv0A2(NULL),
 461 fHResTPCv0C2(NULL),
 462 fHResv0Cv0A2(NULL),
 463 fHResTPCv0A3(NULL),
 464 fHResTPCv0C3(NULL),
 465 fHResv0Cv0A3(NULL),
 466  fHResMA2(NULL),
 467 fHResMC2(NULL),
 468 fHResAC2(NULL),
 469 fHResMA3(NULL),
 470 fHResMC3(NULL),
 471 fHResAC3(NULL),
 472 fPhiRPTPC(NULL),
 473 fPhiRPTPCv3(NULL),
 474 fPhiRPv0A(NULL),
 475 fPhiRPv0C(NULL),
 476 fPhiRPv0Av3(NULL),
 477 fPhiRPv0Cv3(NULL),
 478   fControlConvResoncances(0),
 479   fHistoTPCdEdx(0x0),
 480   fHistoTOFbeta(0x0),
 481   fTPCTOFPion3d(0),
 482   fTPCTOFKaon3d(0),
 483   fTPCTOFProton3d(0),
 484   fPionPt(0),
 485   fPionEta(0),
 486   fPionPhi(0),
 487   fKaonPt(0),
 488   fKaonEta(0),
 489   fKaonPhi(0),
 490   fProtonPt(0),
 491   fProtonEta(0),
 492   fProtonPhi(0),
 493   fCorrelatonTruthPrimary(0),
 494  fCorrelatonTruthPrimarymix(0),
 495   fTHnCorrUNID(0),
 496   fTHnCorrUNIDmix(0),
 497   fTHnCorrID(0),
 498   fTHnCorrIDmix(0),
 499   fTHnCorrIDUNID(0),
 500   fTHnCorrIDUNIDmix(0),
 501   fTHnTrigcount(0),
 502   fTHnTrigcountMCTruthPrim(0),
 503   fPoolMgr(0x0),
 504   fArrayMC(0),
 505   fAnalysisType("AOD"),
 506   fefffilename(""),
 507   ftwoTrackEfficiencyCutDataReco(kTRUE),
 508   twoTrackEfficiencyCutValue(0.02),
 509   fPID(NULL),
 510   fPIDCombined(NULL),
 511   eventno(0),
 512  fPtTOFPIDmin(0.5),
 513   fPtTOFPIDmax(4.0),
 514   fRequestTOFPID(kTRUE),
 515   fPIDType(NSigmaTPCTOF),
 516   fFIllPIDQAHistos(kTRUE),
 517   fNSigmaPID(3),
 518   fBayesCut(0.8),
 519  fdiffPIDcutvalues(kFALSE),
 520  fPIDCutval1(0.0),
 521  fPIDCutval2(0.0),
 522  fPIDCutval3(0.0),
 523  fPIDCutval4(0.0),
 524 fHighPtKaonNSigmaPID(-1),
 525  fHighPtKaonSigma(3.5),
 526   fUseExclusiveNSigma(kFALSE),
 527   fRemoveTracksT0Fill(kFALSE),
 528 fSelectCharge(0),
 529 fTriggerSelectCharge(0),
 530 fAssociatedSelectCharge(0),
 531 fTriggerRestrictEta(-1),
 532 fEtaOrdering(kFALSE),
 533 fCutConversions(kFALSE),
 534 fCutResonances(kFALSE),
 535 fRejectResonanceDaughters(-1),
 536   fOnlyOneEtaSide(0),
 537 fInjectedSignals(kFALSE),
 538   fRemoveWeakDecays(kFALSE),
 539 fRemoveDuplicates(kFALSE),
 540   fapplyTrigefficiency(kFALSE),
 541   fapplyAssoefficiency(kFALSE),
 542   ffillefficiency(kFALSE),
 543  fmesoneffrequired(kFALSE),
 544  fkaonprotoneffrequired(kFALSE),
 545    fAnalysisUtils(0x0),
 546   fDCAXYCut(0)
 547   // The last in the above list should not have a comma after it
 548
 549 {
 550
 551    for ( Int_t i = 0; i < 16; i++) {
 552     fHistQA[i] = NULL;
 553   }
 554
 555 for ( Int_t i = 0; i < 6; i++ ){
 556     fTrackHistEfficiency[i] = NULL;
 557     effcorection[i]=NULL;
 558     //effmap[i]=NULL;
 559   }
 560
 561 for ( Int_t i = 0; i < 2; i++ ){
 562    fTwoTrackDistancePt[i]=NULL;
 563    fTwoTrackDistancePtmix[i]=NULL;
 564 }
 565
 566  for(Int_t ipart=0;ipart<NSpecies;ipart++)
 567     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
 568       fnsigmas[ipart][ipid]=999.;
 569
 570    for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
 571
 572   for(Int_t i = 0; i != 2; ++i)
 573     for(Int_t j = 0; j != 2; ++j)
 574       for(Int_t iC = 0; iC < 9; iC++){
 575         fMeanQ[iC][i][j] = 0.;
 576         fWidthQ[iC][i][j] = 1.;
 577         fMeanQv3[iC][i][j] = 0.;
 578         fWidthQv3[iC][i][j] = 1.;
 579     }
 580
 581   // Constructor
 582   // Define input and output slots here (never in the dummy constructor)
 583   // Input slot #0 works with a TChain - it is connected to the default input container
 584   // Output slot #1 writes into a TH1 container
 585
 586   DefineOutput(1, TList::Class());                                        // for output list
 587   DefineOutput(2, TList::Class());
 588
 589 }
 590
 591 //________________________________________________________________________
 592 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
 593 {
 594   // Destructor. Clean-up the output list, but not the histograms that are put inside
 595   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 596   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 597     delete fOutput;
 598
 599   }
 600
 601 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 602     delete fOutputList;
 603
 604   }
 605
 606   if (fPID) delete fPID;
 607   if (fPIDCombined) delete fPIDCombined;
 608
 609   }
 610 //________________________________________________________________________
 611
 612 //////////////////////////////////////////////////////////////////////////////////////////////////
 613
 614 TH2F* AliTwoParticlePIDCorr::GetHistogram2D(const char * name){
 615   // returns histo named name
 616   return (TH2F*) fOutputList->FindObject(name);
 617 }
 618
 619 //////////////////////////////////////////////////////////////////////////////////////////////////
 620
 621 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
 622
 623 {
 624         //
 625         // Puts the argument in the range [-pi/2,3 pi/2].
 626         //
 627
 628         if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
 629         if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
 630
 631         return DPhi;
 632
 633 }
 634 //________________________________________________________________________
 635 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
 636 {
 637   // Create histograms
 638   // Called once (on the worker node)
 639   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 640   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 641   fPID = inputHandler->GetPIDResponse();
 642
 643   //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
 644
 645 //get the efficiency correction map
 646
 647 // global switch disabling the reference
 648   // (to avoid "Replacing existing TH1" if several wagons are created in train)
 649   Bool_t oldStatus = TH1::AddDirectoryStatus();
 650   TH1::AddDirectory(kFALSE);
 651
 652   const Int_t nPsiTOF = 10;
 653   const Int_t nCentrBin = 9;
 654
 655
 656   fOutput = new TList();
 657   fOutput->SetOwner();  // IMPORTANT!
 658
 659   fOutputList = new TList;
 660   fOutputList->SetOwner();
 661   fOutputList->SetName("PIDQAList");
 662
 663   fList = new TList;
 664   fList->SetOwner();
 665   fList->SetName("EPQAList");
 666
 667   fEventCounter = new TH1F("fEventCounter","EventCounter", 19, 0.5,19.5);
 668   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
 669   fEventCounter->GetXaxis()->SetBinLabel(3,"After PileUP Cut");//only for Data
 670   fEventCounter->GetXaxis()->SetBinLabel(5,"Have A Vertex");
 671   fEventCounter->GetXaxis()->SetBinLabel(7,"After vertex Cut");
 672   fEventCounter->GetXaxis()->SetBinLabel(9,"Getting centrality");
 673   fEventCounter->GetXaxis()->SetBinLabel(11,"After centrality flattening");
 674   fEventCounter->GetXaxis()->SetBinLabel(13,"Within 0-100% centrality");
 675   fEventCounter->GetXaxis()->SetBinLabel(15,"Event Analyzed");
 676   //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
 677   fOutput->Add(fEventCounter);
 678
 679 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
 680 fOutput->Add(fEtaSpectrasso);
 681
 682 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
 683 fOutput->Add(fphiSpectraasso);
 684
 685  if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
 686       fOutput->Add(fCentralityCorrelation);
 687  }
 688
 689 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")
 690   {
 691  TString gCentName[8] = {"V0A","V0C","V0M","V0AEq","V0CEq","V0MEq","CL1","ZNA"};
 692   fHistCentStats = new TH2F("fHistCentStats",
 693                              "Centrality statistics;;Cent percentile",
 694                             8,-0.5,7.5,220,-5,105);
 695   for(Int_t i = 1; i <= 8; i++){
 696     fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
 697     //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
 698   }
 699   fOutput->Add(fHistCentStats);
 700   }
 701
 702 if(fCentralityMethod.EndsWith("_MANUAL"))
 703   {
 704 fhistcentrality=new TH1F("fhistcentrality","referencemultiplicity",30001,-0.5,30000.5);
 705 fOutput->Add(fhistcentrality);
 706   }
 707  else{
 708 fhistcentrality=new TH1F("fhistcentrality","centrality",220,-5,105);
 709 fOutput->Add(fhistcentrality);
 710  }
 711
 712 if(fCentralityMethod.EndsWith("_MANUAL"))
 713   {
 714 TString gmultName[4] = {"V0A_MANUAL","V0C_MANUAL","V0M_MANUAL","TRACKS_MANUAL"};
 715   fHistRefmult = new TH2F("fHistRefmult",
 716                              "Reference multiplicity",
 717                             4,-0.5,3.5,10000,0,20000);
 718   for(Int_t i = 1; i <= 4; i++){
 719     fHistRefmult->GetXaxis()->SetBinLabel(i,gmultName[i-1].Data());
 720     //fHistCentStatsUsed->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
 721   }
 722   fOutput->Add(fHistRefmult);
 723
 724  if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
 725  //TPC vs EQVZERO multiplicity
 726     fHistEQVZEROvsTPCmultiplicity = new TH2F("fHistEQVZEROvsTPCmultiplicity","EqVZERO vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
 727     fHistEQVZEROvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO multiplicity (a.u.)");
 728     fHistEQVZEROvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
 729     fOutput->Add(fHistEQVZEROvsTPCmultiplicity);
 730
 731
 732     fHistEQVZEROAvsTPCmultiplicity = new TH2F("fHistEQVZEROAvsTPCmultiplicity","EqVZERO_A vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
 733     fHistEQVZEROAvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
 734     fHistEQVZEROAvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
 735     fOutput->Add(fHistEQVZEROAvsTPCmultiplicity);
 736
 737
 738     fHistEQVZEROCvsTPCmultiplicity = new TH2F("fHistEQVZEROCvsTPCmultiplicity","EqVZERO_C vs TPC multiplicity",10001,-0.5,10000.5,4001,-0.5,4000.5);
 739     fHistEQVZEROCvsTPCmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
 740     fHistEQVZEROCvsTPCmultiplicity->GetYaxis()->SetTitle("TPC multiplicity (a.u.)");
 741     fOutput->Add(fHistEQVZEROCvsTPCmultiplicity);
 742
 743  //EQVZERO vs VZERO multiplicity
 744   fHistVZEROCvsEQVZEROCmultiplicity = new TH2F("fHistVZEROCvsEQVZEROCmultiplicity","EqVZERO_C vs VZERO_C multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
 745     fHistVZEROCvsEQVZEROCmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
 746     fHistVZEROCvsEQVZEROCmultiplicity->GetYaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
 747     fOutput->Add(fHistVZEROCvsEQVZEROCmultiplicity);
 748
 749
 750 fHistVZEROAvsEQVZEROAmultiplicity = new TH2F("fHistVZEROAvsEQVZEROAmultiplicity","EqVZERO_A vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
 751     fHistVZEROAvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
 752     fHistVZEROAvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
 753     fOutput->Add(fHistVZEROAvsEQVZEROAmultiplicity);
 754
 755
 756   //VZEROC vs VZEROA multiplicity
 757 fHistVZEROCvsVZEROAmultiplicity = new TH2F("fHistVZEROCvsVZEROAmultiplicity","VZERO_C vs VZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
 758     fHistVZEROCvsVZEROAmultiplicity->GetXaxis()->SetTitle("VZERO_C multiplicity (a.u.)");
 759     fHistVZEROCvsVZEROAmultiplicity->GetYaxis()->SetTitle("VZERO_A multiplicity (a.u.)");
 760     fOutput->Add(fHistVZEROCvsVZEROAmultiplicity);
 761
 762
 763
 764   //EQVZEROC vs EQVZEROA multiplicity
 765 fHistEQVZEROCvsEQVZEROAmultiplicity = new TH2F("fHistEQVZEROCvsEQVZEROAmultiplicity","EqVZERO_C vs EqVZERO_A multiplicity",10001,-0.5,10000.5,10001,-0.5,10000.5);
 766     fHistEQVZEROCvsEQVZEROAmultiplicity->GetXaxis()->SetTitle("EqVZERO_C multiplicity (a.u.)");
 767     fHistEQVZEROCvsEQVZEROAmultiplicity->GetYaxis()->SetTitle("EqVZERO_A multiplicity (a.u.)");
 768     fOutput->Add(fHistEQVZEROCvsEQVZEROAmultiplicity);
 769
 770  fHistVZEROSignal = new TH2F("fHistVZEROSignal","VZERO signal vs VZERO channel;VZERO channel; Signal (a.u.)",64,0.5,64.5,3001,-0.5,30000.5);
 771   fOutput->Add(fHistVZEROSignal);
 772  }
 773 }
 774
 775  if(fRequestEventPlane){
 776 //Event plane
 777
 778   fHistPsiMinusPhi = new TH2D("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
 779   fList->Add(fHistPsiMinusPhi);
 780
 781   fEventPlanePID = new TH3F("fEventPlanePID",";centrality;eventplane;PID",20,0.0,100.0,4,-0.5,3.5,4,-0.5,3.5);
 782   fList->Add(fEventPlanePID);
 783
 784  }
 785
 786 if(fCutConversions || fCutResonances)
 787     {
 788 fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
 789  fOutput->Add(fControlConvResoncances);
 790     }
 791
 792 fHistoTPCdEdx = new TH2F("fHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
 793 fOutputList->Add(fHistoTPCdEdx);
 794 fHistoTOFbeta = new TH2F(Form("fHistoTOFbeta"), ";p_{T} (GeV/c);v/c",100, 0., fmaxPt, 500, 0.1, 1.1);
 795   fOutputList->Add(fHistoTOFbeta);
 796
 797    fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
 798    fOutputList->Add(fTPCTOFPion3d);
 799
 800    fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
 801    fOutputList->Add(fTPCTOFKaon3d);
 802
 803    fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
 804    fOutputList->Add(fTPCTOFProton3d);
 805
 806 if(ffillhistQAReco)
 807     {
 808     fPionPt = new TH1F("fPionPt","p_{T} distribution",200,0.,10.);
 809  fOutputList->Add(fPionPt);
 810     fPionEta= new TH1F("fPionEta","#eta distribution",360,-1.8,1.8);
 811  fOutputList->Add(fPionEta);
 812     fPionPhi = new TH1F("fPionPhi","#phi distribution",340,0,6.8);
 813  fOutputList->Add(fPionPhi);
 814
 815     fKaonPt = new TH1F("fKaonPt","p_{T} distribution",200,0.,10.);
 816  fOutputList->Add(fKaonPt);
 817     fKaonEta= new TH1F("fKaonEta","#eta distribution",360,-1.8,1.8);
 818  fOutputList->Add(fKaonEta);
 819     fKaonPhi = new TH1F("fKaonPhi","#phi distribution",340,0,6.8);
 820  fOutputList->Add(fKaonPhi);
 821
 822     fProtonPt = new TH1F("fProtonPt","p_{T} distribution",200,0.,10.);
 823  fOutputList->Add(fProtonPt);
 824     fProtonEta= new TH1F("fProtonEta","#eta distribution",360,-1.8,1.8);
 825  fOutputList->Add(fProtonEta);
 826     fProtonPhi= new TH1F("fProtonPhi","#phi distribution",340,0,6.8);
 827  fOutputList->Add(fProtonPhi);
 828     }
 829
 830   fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
 831   fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
 832   fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
 833   fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx  After Cut ", 50, -5., 5.);
 834   fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
 835   fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
 836   fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
 837   fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
 838   fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
 839   fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
 840   fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
 841   fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
 842   fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
 843  fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
 844  fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
 845  fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
 846
 847 for(Int_t i = 0; i < 16; i++)
 848     {
 849       fOutput->Add(fHistQA[i]);
 850     }
 851
 852     fPriHistShare = new TH1F ("fPriHistShare","Shared clusters, primaries;#shared clusters;counts",160,0,160);
 853     fOutput->Add(fPriHistShare);
 854
 855    Int_t eventplaneaxis=0;
 856
 857    if (fRequestEventPlane) eventplaneaxis=1;
 858
 859    kTrackVariablesPair=6+SetChargeAxis+eventplaneaxis;
 860
 861    if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
 862
 863  if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7+SetChargeAxis+eventplaneaxis;
 864
 865  if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8+SetChargeAxis+eventplaneaxis;
 866
 867
 868 // two particle histograms
 869   Int_t anaSteps   = 1;       // analysis steps
 870   const char* title = "d^{2}N_{ch}/d#varphid#eta";
 871
 872   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
 873   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables
 874   TString* axisTitlePair;  // axis titles for track variables
 875   axisTitlePair=new TString[kTrackVariablesPair];
 876
 877  TString defaultBinningStr;
 878   defaultBinningStr =   "eta: -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\n"
 879     "p_t_assoc: 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0,10.0\n"
 880     "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
 881     "p_t_eff:0.0,0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0,5.5, 6.0, 7.0, 8.0,9.0,10.0\n"
 882     "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
 883   "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, 0.0, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n" // this binning starts at -pi/2 and is modulo 3
 884         "delta_eta: -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,2.1, 2.2, 2.3, 2.4\n"
 885       "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
 886
 887  if(fRequestEventPlane){
 888    defaultBinningStr += "eventPlane: -0.5,0.5,1.5,2.5,3.5\n"; // Event Plane Bins (Psi: -0.5->0.5 (in plane), 0.5->1.5 (intermediate), 1.5->2.5 (out of plane), 2.5->3.5 (rest))
 889   }
 890
 891   if(fcontainPIDtrig){
 892       defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
 893   }
 894   if(fcontainPIDasso){
 895       defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
 896   }
 897
 898   if(SetChargeAxis==2){
 899       defaultBinningStr += "TrigCharge: -2.0,0.0,2.0\n"; // course
 900       defaultBinningStr += "AssoCharge: -2.0,0.0,2.0\n"; // course
 901   }
 902  // =========================================================
 903   // Customization (adopted from AliUEHistograms)
 904   // =========================================================
 905
 906   TObjArray* lines = defaultBinningStr.Tokenize("\n");
 907   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
 908   {
 909     TString line(lines->At(i)->GetName());
 910     TString tag = line(0, line.Index(":")+1);
 911     if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
 912       fBinningString += line + "\n";
 913     else
 914       AliInfo(Form("Using custom binning for %s", tag.Data()));
 915   }
 916   delete lines;
 917   fBinningString += fCustomBinning;
 918
 919   AliInfo(Form("Used AliTHn Binning:\n%s",fBinningString.Data()));
 920
 921  //  =========================================================
 922   // Now set the bins
 923   // =========================================================
 924
 925     dBinsPair[0]       = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
 926     axisTitlePair[0]   = " multiplicity";
 927
 928     dBinsPair[1]     = GetBinning(fBinningString, "vertex", iBinPair[1]);
 929     axisTitlePair[1]  = "v_{Z} (cm)";
 930
 931     dBinsPair[2]     = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
 932     axisTitlePair[2]    = "p_{T,trig.} (GeV/c)";
 933
 934     dBinsPair[3]     = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
 935     axisTitlePair[3]    = "p_{T,assoc.} (GeV/c)";
 936
 937     dBinsPair[4]       = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
 938     axisTitlePair[4]   = "#Delta#eta";
 939
 940     dBinsPair[5]       = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
 941     axisTitlePair[5]   = "#Delta#varphi (rad)";
 942
 943     Int_t dim_val=6;
 944
 945     if(fRequestEventPlane){
 946     dBinsPair[dim_val]       = GetBinning(fBinningString, "eventPlane", iBinPair[dim_val]);
 947     axisTitlePair[dim_val]   = "#varphi - #Psi_{2} (a.u.)";
 948     dim_val=7;
 949     }
 950
 951     if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
 952     dBinsPair[dim_val]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val]);
 953     axisTitlePair[dim_val]   = "TrigCharge";
 954
 955     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+1]);
 956     axisTitlePair[dim_val+1]   = "AssoCharge";
 957     }
 958
 959  if(fcontainPIDtrig && !fcontainPIDasso){
 960     dBinsPair[dim_val]       = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
 961     axisTitlePair[dim_val]   = "PIDTrig";
 962     if(SetChargeAxis==2){
 963     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
 964     axisTitlePair[dim_val+1]   = "TrigCharge";
 965
 966     dBinsPair[dim_val+2]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
 967     axisTitlePair[dim_val+2]   = "AssoCharge";
 968     }
 969  }
 970
 971  if(!fcontainPIDtrig && fcontainPIDasso){
 972     dBinsPair[dim_val]       = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val]);
 973     axisTitlePair[dim_val]   = "PIDAsso";
 974
 975  if(SetChargeAxis==2){
 976     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+1]);
 977     axisTitlePair[dim_val+1]   = "TrigCharge";
 978
 979     dBinsPair[dim_val+2]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+2]);
 980     axisTitlePair[dim_val+2]   = "AssoCharge";
 981     }
 982  }
 983
 984 if(fcontainPIDtrig && fcontainPIDasso){
 985
 986     dBinsPair[dim_val]       = GetBinning(fBinningString, "PIDTrig", iBinPair[dim_val]);
 987     axisTitlePair[dim_val]   = "PIDTrig";
 988
 989     dBinsPair[dim_val+1]       = GetBinning(fBinningString, "PIDAsso", iBinPair[dim_val+1]);
 990     axisTitlePair[dim_val+1]   = "PIDAsso";
 991
 992     if(SetChargeAxis==2){
 993     dBinsPair[dim_val+2]       = GetBinning(fBinningString, "TrigCharge", iBinPair[dim_val+2]);
 994     axisTitlePair[dim_val+2]   = "TrigCharge";
 995
 996     dBinsPair[dim_val+3]       = GetBinning(fBinningString, "AssoCharge", iBinPair[dim_val+3]);
 997     axisTitlePair[dim_val+3]   = "AssoCharge";
 998     }
 999  }
1000
1001         Int_t nEtaBin = -1;
1002         Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
1003
1004         Int_t nPteffbin = -1;
1005         Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
1006
1007
1008         fminPtTrig=dBinsPair[2][0];
1009         fmaxPtTrig=dBinsPair[2][iBinPair[2]];
1010         fminPtAsso=dBinsPair[3][0];
1011         fmaxPtAsso=dBinsPair[3][iBinPair[3]];
1012         fmincentmult=dBinsPair[0][0];
1013         fmaxcentmult=dBinsPair[0][iBinPair[0]];
1014
1015         //event pool manager
1016 Int_t MaxNofEvents=1000;
1017 const Int_t NofVrtxBins=10+(1+10)*2;
1018 Double_t ZvrtxBins[NofVrtxBins+1]={ -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10,
1019                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110,
1020                                     190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210};
1021
1022 if(fCentralityMethod.EndsWith("_MANUAL"))
1023    {
1024  const Int_t NofCentBins=10;
1025  Double_t CentralityBins[NofCentBins+1]={0.,9.,14.,19.,26.,34.,44.,58.,80.,500.,1000.};//Is This binning is fine for pp, or we don't require them....
1026 if(fRequestEventPlane){
1027     // Event plane angle (Psi) bins
1028   /*
1029     Double_t* psibins = NULL;
1030     Int_t nPsiBins;
1031     psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1032   */
1033  const Int_t  nPsiBins=6;
1034  Double_t psibins[nPsiBins+1]={0.0*TMath::DegToRad(), 30.0*TMath::DegToRad(), 60.0*TMath::DegToRad(), 90.0*TMath::DegToRad(), 120.0*TMath::DegToRad(),150.0*TMath::DegToRad(),180.1*TMath::DegToRad()};
1035 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1036 // if(psibins)  delete [] psibins;
1037                                     }
1038
1039  else{
1040  const Int_t  nPsiBinsd=1;
1041  Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1042 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1043
1044 // fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1045  }
1046 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1047
1048    }
1049  else
1050    {
1051  const Int_t  NofCentBins=15;
1052 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
1053  if(fRequestEventPlane){
1054     // Event plane angle (Psi) bins
1055    /*
1056     Double_t* psibins = NULL;
1057     Int_t nPsiBins;
1058     psibins = GetBinning(fBinningString, "eventPlane", nPsiBins);
1059    */
1060  const Int_t  nPsiBins=6;
1061  Double_t psibins[nPsiBins+1]={0.0*TMath::DegToRad(), 30.0*TMath::DegToRad(), 60.0*TMath::DegToRad(), 90.0*TMath::DegToRad(), 120.0*TMath::DegToRad(),150.0*TMath::DegToRad(),180.1*TMath::DegToRad()};
1062 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBins, psibins);
1063 // if(psibins)  delete [] psibins;
1064                                     }
1065
1066  else{
1067 const Int_t  nPsiBinsd=1;
1068  Double_t psibinsd[nPsiBinsd+1]={0.0, 2000.0};
1069 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins, nPsiBinsd, psibinsd);
1070
1071 //fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
1072  }
1073 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
1074    }
1075
1076
1077    if(!fPoolMgr){
1078       AliError("Event Mixing required, but Pool Manager not initialized...");
1079       return;
1080     }
1081
1082         //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
1083         //fmaxPtComboeff=fmaxPtTrig;
1084 //THnSparses for calculation of efficiency
1085
1086  if((fAnalysisType =="MCAOD") && ffillefficiency) {
1087 TString Histrename;
1088   Int_t effbin[4];
1089   effbin[0]=iBinPair[0];
1090   effbin[1]=iBinPair[1];
1091   effbin[2]=nPteffbin;
1092   effbin[3]=nEtaBin;
1093   Int_t effsteps=5;//for each species type::primMCParticles(0),primRecoTracksMatched(1),allRecoTracksMatched(2),primRecoTracksMatchedPID(3),allRecoTracksMatchedPID(4)
1094 for(Int_t jj=0;jj<6;jj++)//PID type binning
1095     {
1096      if(jj==5) effsteps=3;//for unidentified particles
1097   Histrename="fTrackHistEfficiency";Histrename+=jj;
1098   fTrackHistEfficiency[jj] = new AliTHn(Histrename.Data(), "Tracking efficiency", effsteps, 4, effbin);
1099   fTrackHistEfficiency[jj]->SetBinLimits(0, dBinsPair[0]);
1100   fTrackHistEfficiency[jj]->SetVarTitle(0, "Centrality");
1101   fTrackHistEfficiency[jj]->SetBinLimits(1, dBinsPair[1]);
1102   fTrackHistEfficiency[jj]->SetVarTitle(1, "zvtx");
1103   fTrackHistEfficiency[jj]->SetBinLimits(2, Pteff);
1104   fTrackHistEfficiency[jj]->SetVarTitle(2, "p_{T} (GeV/c)");
1105   fTrackHistEfficiency[jj]->SetBinLimits(3, EtaBin);
1106   fTrackHistEfficiency[jj]->SetVarTitle(3, "#eta");
1107   fOutput->Add(fTrackHistEfficiency[jj]);
1108     }
1109  }
1110
1111 //AliThns for Correlation plots(data &  MC)
1112
1113      if(ffilltrigassoUNID)
1114        {
1115     fTHnCorrUNID = new AliTHn("fTHnCorrUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1116 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1117     fTHnCorrUNID->SetBinLimits(j, dBinsPair[j]);
1118     fTHnCorrUNID->SetVarTitle(j, axisTitlePair[j]);
1119   }
1120   fOutput->Add(fTHnCorrUNID);
1121
1122  fTHnCorrUNIDmix = new AliTHn("fTHnCorrUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1123 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1124     fTHnCorrUNIDmix->SetBinLimits(j, dBinsPair[j]);
1125     fTHnCorrUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1126   }
1127   fOutput->Add(fTHnCorrUNIDmix);
1128        }
1129
1130      if(ffilltrigIDassoID)
1131        {
1132 fTHnCorrID = new AliTHn("fTHnCorrID", title, anaSteps, kTrackVariablesPair, iBinPair);
1133 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1134     fTHnCorrID->SetBinLimits(j, dBinsPair[j]);
1135     fTHnCorrID->SetVarTitle(j, axisTitlePair[j]);
1136   }
1137   fOutput->Add(fTHnCorrID);
1138
1139 fTHnCorrIDmix = new AliTHn("fTHnCorrIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1140 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1141     fTHnCorrIDmix->SetBinLimits(j, dBinsPair[j]);
1142     fTHnCorrIDmix->SetVarTitle(j, axisTitlePair[j]);
1143   }
1144   fOutput->Add(fTHnCorrIDmix);
1145        }
1146
1147      if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
1148        {
1149 fTHnCorrIDUNID = new AliTHn("fTHnCorrIDUNID", title, anaSteps, kTrackVariablesPair, iBinPair);
1150 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1151     fTHnCorrIDUNID->SetBinLimits(j, dBinsPair[j]);
1152     fTHnCorrIDUNID->SetVarTitle(j, axisTitlePair[j]);
1153   }
1154   fOutput->Add(fTHnCorrIDUNID);
1155
1156
1157 fTHnCorrIDUNIDmix = new AliTHn("fTHnCorrIDUNIDmix", title, anaSteps, kTrackVariablesPair, iBinPair);
1158 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1159     fTHnCorrIDUNIDmix->SetBinLimits(j, dBinsPair[j]);
1160     fTHnCorrIDUNIDmix->SetVarTitle(j, axisTitlePair[j]);
1161   }
1162   fOutput->Add(fTHnCorrIDUNIDmix);
1163        }
1164
1165
1166
1167   //ThnSparse for Correlation plots(truth MC)
1168      if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
1169
1170 fCorrelatonTruthPrimary = new AliTHn("fCorrelatonTruthPrimary", title, anaSteps, kTrackVariablesPair, iBinPair);
1171 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1172     fCorrelatonTruthPrimary->SetBinLimits(j, dBinsPair[j]);
1173     fCorrelatonTruthPrimary->SetVarTitle(j, axisTitlePair[j]);
1174   }
1175   fOutput->Add(fCorrelatonTruthPrimary);
1176
1177
1178 fCorrelatonTruthPrimarymix = new AliTHn("fCorrelatonTruthPrimarymix", title, anaSteps, kTrackVariablesPair, iBinPair);
1179 for (Int_t j=0; j<kTrackVariablesPair; j++) {
1180     fCorrelatonTruthPrimarymix->SetBinLimits(j, dBinsPair[j]);
1181     fCorrelatonTruthPrimarymix->SetVarTitle(j, axisTitlePair[j]);
1182   }
1183   fOutput->Add(fCorrelatonTruthPrimarymix);
1184  }
1185
1186     //binning for trigger no. counting
1187
1188      Int_t ChargeAxis=0;
1189      if(SetChargeAxis==2) ChargeAxis=1;
1190
1191         Int_t* fBinst;
1192         Int_t dims=3+ChargeAxis+eventplaneaxis;
1193         if(fcontainPIDtrig) dims=4+ChargeAxis+eventplaneaxis;
1194         fBinst= new Int_t[dims];
1195    Double_t* dBinsTrig[dims];    // bins for track variables
1196    TString* axisTitleTrig;  // axis titles for track variables
1197    axisTitleTrig=new TString[dims];
1198
1199         for(Int_t i=0; i<3;i++)
1200           {
1201             fBinst[i]=iBinPair[i];
1202             dBinsTrig[i]=dBinsPair[i];
1203             axisTitleTrig[i]=axisTitlePair[i];
1204           }
1205         Int_t dim_val_trig=3;
1206     if(fRequestEventPlane){
1207       fBinst[dim_val_trig]=iBinPair[6];//if fRequestEventPlane=TRUE, dim_val already becomes 7.
1208       dBinsTrig[dim_val_trig]=dBinsPair[6];
1209       axisTitleTrig[dim_val_trig]=axisTitlePair[6];
1210       dim_val_trig=4;
1211     }
1212
1213 if(!fcontainPIDtrig && !fcontainPIDasso && ChargeAxis==1){
1214 fBinst[dim_val_trig]=iBinPair[dim_val];
1215 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1216 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1217     }
1218
1219 if(fcontainPIDtrig && !fcontainPIDasso){
1220 fBinst[dim_val_trig]=iBinPair[dim_val];
1221 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1222 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1223     if(ChargeAxis==1){
1224 fBinst[dim_val_trig+1]=iBinPair[dim_val+1];
1225 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+1];
1226 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+1];
1227     }
1228  }
1229
1230  if(!fcontainPIDtrig && fcontainPIDasso){
1231  if(ChargeAxis==1){
1232     fBinst[dim_val_trig]=iBinPair[dim_val+1];
1233 dBinsTrig[dim_val_trig]=dBinsPair[dim_val+1];
1234 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val+1];
1235     }
1236  }
1237
1238 if(fcontainPIDtrig && fcontainPIDasso){
1239   fBinst[dim_val_trig]=iBinPair[dim_val];
1240 dBinsTrig[dim_val_trig]=dBinsPair[dim_val];
1241 axisTitleTrig[dim_val_trig]=axisTitlePair[dim_val];
1242     if(ChargeAxis==1){
1243 fBinst[dim_val_trig+1]=iBinPair[dim_val+2];
1244 dBinsTrig[dim_val_trig+1]=dBinsPair[dim_val+2];
1245 axisTitleTrig[dim_val_trig+1]=axisTitlePair[dim_val+2];
1246     }
1247     }
1248
1249   //ThSparse for trigger counting(data & reco MC)
1250   if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
1251           {
1252             fTHnTrigcount = new  AliTHn("fTHnTrigcount", "fTHnTrigcount", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1253    for(Int_t i=0; i<dims;i++){
1254     fTHnTrigcount->SetBinLimits(i, dBinsTrig[i]);
1255     fTHnTrigcount->SetVarTitle(i, axisTitleTrig[i]);
1256   }
1257   fOutput->Add(fTHnTrigcount);
1258           }
1259
1260   if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
1261   //AliTHns for trigger counting(truth MC)
1262   fTHnTrigcountMCTruthPrim = new  AliTHn("fTHnTrigcountMCTruthPrim", "fTHnTrigcountMCTruthPrim", 2, dims, fBinst); //2 steps;;;;0->same event;;;;;1->mixed event
1263  for(Int_t i=0; i<dims;i++){
1264     fTHnTrigcountMCTruthPrim->SetBinLimits(i, dBinsTrig[i]);
1265     fTHnTrigcountMCTruthPrim->SetVarTitle(i, axisTitleTrig[i]);
1266   }
1267   fOutput->Add(fTHnTrigcountMCTruthPrim);
1268  }
1269
1270 if(fAnalysisType=="MCAOD"){
1271   if(ffillhistQATruth)
1272     {
1273   MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
1274   fOutputList->Add(MCtruthpt);
1275
1276   MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
1277   fOutputList->Add(MCtrutheta);
1278
1279   MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
1280   fOutputList->Add(MCtruthphi);
1281
1282   MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
1283   fOutputList->Add(MCtruthpionpt);
1284
1285   MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
1286   fOutputList->Add(MCtruthpioneta);
1287
1288   MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
1289   fOutputList->Add(MCtruthpionphi);
1290
1291   MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
1292   fOutputList->Add(MCtruthkaonpt);
1293
1294   MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
1295   fOutputList->Add(MCtruthkaoneta);
1296
1297   MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
1298   fOutputList->Add(MCtruthkaonphi);
1299
1300   MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
1301   fOutputList->Add(MCtruthprotonpt);
1302
1303   MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
1304   fOutputList->Add(MCtruthprotoneta);
1305
1306   MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
1307   fOutputList->Add(MCtruthprotonphi);
1308     }
1309  fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
1310   fOutputList->Add(fPioncont);
1311
1312  fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
1313   fOutputList->Add(fKaoncont);
1314
1315  fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
1316   fOutputList->Add(fProtoncont);
1317
1318 fUNIDcont=new TH2F("fUNIDcont","fUNIDcont",10,-0.5,9.5,100,0.,10.);
1319   fOutputList->Add(fUNIDcont);
1320   }
1321
1322 fEventno=new TH2F("fEventno","fEventno",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1323  fEventno->GetXaxis()->SetTitle("Centrality");
1324  fEventno->GetYaxis()->SetTitle("Z_Vtx");
1325 fOutput->Add(fEventno);
1326 fEventnobaryon=new TH2F("fEventnobaryon","fEventnobaryon",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1327  fEventnobaryon->GetXaxis()->SetTitle("Centrality");
1328  fEventnobaryon->GetYaxis()->SetTitle("Z_Vtx");
1329 fOutput->Add(fEventnobaryon);
1330 fEventnomeson=new TH2F("fEventnomeson","fEventnomeson",iBinPair[0], dBinsPair[0],iBinPair[1],dBinsPair[1]);
1331  fEventnomeson->GetXaxis()->SetTitle("Centrality");
1332  fEventnomeson->GetYaxis()->SetTitle("Z_Vtx");
1333 fOutput->Add(fEventnomeson);
1334
1335 fhistJetTrigestimate=new TH2F("fhistJetTrigestimate","fhistJetTrigestimate",iBinPair[0],dBinsPair[0],6,-0.5,5.5);
1336 fOutput->Add(fhistJetTrigestimate);
1337
1338    fTwoTrackDistancePtdip = new TH3F("fTwoTrackDistancePtdip", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
1339   fOutput->Add(fTwoTrackDistancePtdip);
1340
1341 fTwoTrackDistancePtdipmix = new TH3F("fTwoTrackDistancePtdipmix", ";#Delta#eta;#Delta#varphi;#Delta p_{T}", 36, -1.8, 1.8, 72,-TMath::Pi()/2, 3*TMath::Pi()/2, 40, 0, 10);
1342   fOutput->Add(fTwoTrackDistancePtdipmix);
1343
1344   TString Histttrname;
1345 for(Int_t jj=0;jj<2;jj++)// PID type binning
1346     {
1347   Histttrname="fTwoTrackDistancePt";Histttrname+=jj;
1348   fTwoTrackDistancePt[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
1349   fOutput->Add(fTwoTrackDistancePt[jj]);
1350
1351  Histttrname="fTwoTrackDistancePtmix";Histttrname+=jj;
1352   fTwoTrackDistancePtmix[jj] = new TH3F(Histttrname.Data(), ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
1353   fOutput->Add(fTwoTrackDistancePtmix[jj]);
1354     }
1355 //Mixing
1356 //DefineEventPool();
1357
1358   if(fapplyTrigefficiency || fapplyAssoefficiency)
1359    {
1360      const Int_t nDimt = 4;//       cent zvtx  pt   eta
1361      Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
1362      Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
1363      Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
1364
1365   TString Histrexname;
1366 for(Int_t jj=0;jj<6;jj++)// PID type binning
1367     {
1368   Histrexname="effcorection";Histrexname+=jj;
1369   effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
1370   effcorection[jj]->Sumw2();
1371   effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
1372   effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
1373   effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
1374   effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
1375   effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
1376   effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
1377   effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
1378   effcorection[jj]->GetAxis(3)->SetTitle("#eta");
1379   fOutput->Add(effcorection[jj]);
1380     }
1381 // TFile *fsifile = new TFile(fefffilename,"READ");
1382
1383  if (TString(fefffilename).BeginsWith("alien:"))
1384     TGrid::Connect("alien:");
1385  TFile *fileT=TFile::Open(fefffilename);
1386  TString Nameg;
1387 for(Int_t jj=0;jj<6;jj++)//type binning
1388     {
1389 Nameg="effmap";Nameg+=jj;
1390 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
1391 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
1392
1393 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
1394     }
1395 //fsifile->Close();
1396 fileT->Close();
1397
1398    }
1399
1400   //*************************************************************EP plots***********************************************//
1401   if(fRequestEventPlane){
1402   // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
1403   // v2
1404   fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
1405   fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
1406   fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
1407
1408   fList->Add(fHResTPCv0A2);
1409   fList->Add(fHResTPCv0C2);
1410   fList->Add(fHResv0Cv0A2);
1411
1412   // v3
1413   fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
1414   fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
1415   fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
1416
1417   fList->Add(fHResTPCv0A3);
1418   fList->Add(fHResTPCv0C3);
1419   fList->Add(fHResv0Cv0A3);
1420
1421   // MC as in the dataEP resolution (but using MC tracks)
1422   if(fAnalysisType == "MCAOD"  && fV2){
1423     fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
1424     fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
1425     fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
1426     fList->Add(fHResMA2);
1427     fList->Add(fHResMC2);
1428     fList->Add(fHResAC2);
1429   }
1430   if(fAnalysisType == "MCAOD" && fV3){
1431     fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
1432     fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
1433     fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
1434     fList->Add(fHResMA3);
1435     fList->Add(fHResMC3);
1436     fList->Add(fHResAC3);
1437   }
1438
1439
1440   // V0A and V0C event plane distributions
1441   //v2
1442   fPhiRPTPC = new TH2F("fPhiRPTPCv2","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1443   fPhiRPTPCv3 = new TH2F("fPhiRPTPCv3","#phi distribution of EP TPC;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1444   fList->Add(fPhiRPTPC);
1445   fList->Add(fPhiRPTPCv3);
1446
1447   fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1448   fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1449   fList->Add(fPhiRPv0A);
1450   fList->Add(fPhiRPv0C);
1451
1452   //v3
1453   fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1454   fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
1455   fList->Add(fPhiRPv0Av3);
1456   fList->Add(fPhiRPv0Cv3);
1457
1458   fHistEventPlaneTruth = new TH2F("fHistEventPlaneTruth","#phi distribution of EP MCTRUTHheader;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
1459   fList->Add(fHistEventPlaneTruth);
1460
1461   }
1462
1463   //*****************************************************PIDQA histos*****************************************************//
1464
1465
1466   //nsigma plot
1467   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1468     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1469       Double_t miny=-30;
1470       Double_t maxy=30;
1471       if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1472       TH2F *fHistoNSigma=new TH2F(Form("NSigma_%d_%d",ipart,ipid),Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1473       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1474       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1475       fOutputList->Add(fHistoNSigma);
1476     }
1477   }
1478
1479   //nsigmaRec plot
1480   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1481     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1482       Double_t miny=-10;
1483       Double_t maxy=10;
1484       if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1485       TH2F *fHistoNSigma=new TH2F(Form("NSigmaRec_%d_%d",ipart,ipid),
1486                                   Form("n#sigma for reconstructed %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1487       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1488       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1489       fOutputList->Add(fHistoNSigma);
1490     }
1491   }
1492
1493   //BayesRec plot
1494   if(fPIDType==Bayes){//use bayesianPID
1495     fPIDCombined = new AliPIDCombined();
1496     fPIDCombined->SetDefaultTPCPriors();//****************************************Need to know about it
1497
1498   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1499     Double_t miny=0.;
1500     Double_t maxy=1;
1501     TH2F *fHistoBayes=new TH2F(Form("BayesRec_%d",ipart),
1502                                Form("probability for reconstructed %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1503     fHistoBayes->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1504     fHistoBayes->GetYaxis()->SetTitle(Form("Bayes prob %s",kParticleSpeciesName[ipart]));
1505     fOutputList->Add(fHistoBayes);
1506
1507
1508    TH2F *fHistoBayesTPC=new TH2F(Form("probBayes_TPC_%d",ipart),
1509                                Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1510     fHistoBayesTPC->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1511     fHistoBayesTPC->GetYaxis()->SetTitle(Form("Bayes prob TPC %s",kParticleSpeciesName[ipart]));
1512     fOutputList->Add(fHistoBayesTPC);
1513
1514   TH2F *fHistoBayesTOF=new TH2F(Form("probBayes_TOF_%d",ipart),
1515                                Form("probability for Tracks as %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1516     fHistoBayesTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1517     fHistoBayesTOF->GetYaxis()->SetTitle(Form("Bayes prob TOF %s",kParticleSpeciesName[ipart]));
1518     fOutputList->Add(fHistoBayesTOF);
1519
1520  TH2F *fHistoBayesTPCTOF=new TH2F(Form("probBayes_TPCTOF_%d",ipart),
1521                                Form("probability for Tracks as  %s",kParticleSpeciesName[ipart]),200,0,10,500,miny,maxy);
1522     fHistoBayesTPCTOF->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1523     fHistoBayesTPCTOF->GetYaxis()->SetTitle(Form("Bayes prob TPCTOF %s",kParticleSpeciesName[ipart]));
1524     fOutputList->Add(fHistoBayesTPCTOF);
1525   }
1526   }
1527
1528   //nsigma separation power plot
1529     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1530  Double_t miny=0;
1531  Double_t maxy=10;
1532    TH2F *Pi_Ka_sep=new TH2F(Form("Pi_Ka_sep_%d",ipid),
1533                                Form("Pi_Ka separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1534     Pi_Ka_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1535     Pi_Ka_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1536     fOutputList->Add(Pi_Ka_sep);
1537
1538    TH2F *Pi_Pr_sep=new TH2F(Form("Pi_Pr_sep_%d",ipid),
1539                                Form("Pi_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1540     Pi_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1541     Pi_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1542     fOutputList->Add(Pi_Pr_sep);
1543
1544     TH2F *Ka_Pr_sep=new TH2F(Form("Ka_Pr_sep_%d",ipid),
1545                                Form("Ka_Pr separation in %s",kPIDTypeName[ipid]),50,0,10,200,miny,maxy);
1546     Ka_Pr_sep->GetXaxis()->SetTitle("P_{T} (GeV/C)");
1547     Ka_Pr_sep->GetYaxis()->SetTitle(Form("expected seaparation(n#sigma) in %s",kPIDTypeName[ipid]));
1548     fOutputList->Add(Ka_Pr_sep);
1549     }
1550
1551   //nsigmaDC plot
1552   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1553     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1554       Double_t miny=-10;
1555       Double_t maxy=10;
1556       if(ipid==NSigmaTPCTOF){miny=0;maxy=20;}
1557       TH2F *fHistoNSigma=new TH2F(Form("NSigmaDC_%d_%d",ipart,ipid),
1558                                   Form("n#sigma for double counting %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1559       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1560       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1561       fOutputList->Add(fHistoNSigma);
1562     }
1563   }
1564
1565   //nsigmaMC plot
1566  if (fAnalysisType == "MCAOD"){
1567   for(Int_t ipart=0;ipart<NSpecies;ipart++){
1568     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
1569       Double_t miny=-30;
1570       Double_t maxy=30;
1571       if(ipid==NSigmaTPCTOF){miny=0;maxy=50;}
1572       TH2F *fHistoNSigma=new TH2F(Form("NSigmaMC_%d_%d",ipart,ipid),
1573                                   Form("n#sigma for MC %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]),200,0,10,500,miny,maxy);
1574       fHistoNSigma->GetXaxis()->SetTitle("P_{T} (GeV / c)");
1575       fHistoNSigma->GetYaxis()->SetTitle(Form("n#sigma %s %s",kParticleSpeciesName[ipart],kPIDTypeName[ipid]));
1576       fOutputList->Add(fHistoNSigma);
1577     }
1578   }
1579   }
1580   //PID signal plot
1581   for(Int_t idet=0;idet<fNDetectors;idet++){
1582     for(Int_t ipart=0;ipart<NSpecies;ipart++){
1583       Double_t maxy=500;
1584       if(idet==fTOF)maxy=1.1;
1585       TH2F *fHistoPID=new TH2F(Form("PID_%d_%d",idet,ipart),Form("%s signal - %s",kDetectorName[idet],kParticleSpeciesName[ipart]),200,0,10,500,-maxy,maxy);
1586       fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1587       fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1588       fOutputList->Add(fHistoPID);
1589     }
1590   }
1591   //PID signal plot, before PID cut
1592   for(Int_t idet=0;idet<fNDetectors;idet++){
1593     Double_t maxy=500;
1594     if(idet==fTOF)maxy=1.1;
1595     TH2F *fHistoPID=new TH2F(Form("PIDAll_%d",idet),Form("%s signal",kDetectorName[idet]),200,0,10,500,-maxy,maxy);
1596     fHistoPID->GetXaxis()->SetTitle("P (GeV / c)");
1597     fHistoPID->GetYaxis()->SetTitle(Form("%s signal",kDetectorName[idet]));
1598     fOutputList->Add(fHistoPID);
1599   }
1600
1601   PostData(1, fOutput);              // Post data for ALL output slots >0 here, to get at least an empty histogram
1602   PostData(2, fOutputList);
1603   if(fRequestEventPlane) PostData(3, fList);
1604   AliInfo("Finished setting up the Output");
1605
1606   TH1::AddDirectory(oldStatus);
1607
1608
1609
1610 }
1611 //-------------------------------------------------------------------------------
1612 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
1613
1614
1615   if(fAnalysisType == "AOD") {
1616
1617     doAODevent();
1618
1619   }//AOD--analysis-----
1620
1621   else if(fAnalysisType == "MCAOD") {
1622
1623     doMCAODevent();
1624
1625   }
1626
1627   else return;
1628
1629 }
1630 //-------------------------------------------------------------------------
1631 void AliTwoParticlePIDCorr::doMCAODevent()
1632 {
1633   AliVEvent *event = InputEvent();
1634   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1635   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1636   if (!aod) {
1637     AliError("Cannot get the AOD event");
1638     return;
1639   }
1640
1641 // count all events(physics triggered)
1642   fEventCounter->Fill(1);
1643
1644     evplaneMC=999.;
1645     fgPsi2v0aMC=999.;
1646     fgPsi2v0cMC=999.;
1647     fgPsi2tpcMC=999.;
1648     fgPsi3v0aMC=999.;
1649     fgPsi3v0cMC=999.;
1650     fgPsi3tpcMC=999.;
1651     gReactionPlane = 999.;
1652
1653  // get centrality object and check quality(valid for p-Pb and Pb-Pb; coming soon for pp 7 TeV)
1654   Double_t cent_v0=-1.0;
1655   Double_t effcent=1.0;
1656   Double_t refmultReco =0.0;
1657
1658 //check the PIDResponse handler
1659   if (!fPID) return;
1660
1661 // get mag. field required for twotrack efficiency cut
1662  Float_t bSign = 0;
1663  bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
1664
1665  //check for TClonesArray(truth track MC information)
1666 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1667   if (!fArrayMC) {
1668     AliFatal("Error: MC particles branch not found!\n");
1669     return;
1670   }
1671
1672   //check for AliAODMCHeader(truth event MC information)
1673   AliAODMCHeader *header=NULL;
1674   header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1675   if(!header) {
1676     printf("MC header branch not found!\n");
1677     return;
1678   }
1679
1680 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
1681 Float_t zVtxmc =header->GetVtxZ();
1682  if(TMath::Abs(zVtxmc)>fzvtxcut) return;
1683
1684  // For productions with injected signals, figure out above which label to skip particles/tracks
1685
1686  if (fInjectedSignals)
1687   {
1688     AliGenEventHeader* eventHeader = 0;
1689     Int_t headers = 0;
1690
1691 // AOD
1692       if (!header)
1693       AliFatal("fInjectedSignals set but no MC header found");
1694
1695       headers = header->GetNCocktailHeaders();
1696       eventHeader = header->GetCocktailHeader(0);
1697
1698  if (!eventHeader)
1699     {
1700       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
1701       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
1702       AliError("First event header not found. Skipping this event.");
1703       //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
1704       return;
1705     }
1706 skipParticlesAbove = eventHeader->NProduced();
1707     AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
1708   }
1709
1710  if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL"))
1711    {
1712  //make the event selection with reco vertex cut and centrality cut and return the value of the centrality
1713      Double_t refmultTruth = GetAcceptedEventMultiplicity(aod,kTRUE);  //incase of ref multiplicity it will return the truth MC ref mullt value; need to determine the ref mult value separately for reco Mc case; in case of centrality this is final and fine
1714      refmultReco = GetAcceptedEventMultiplicity(aod,kFALSE);
1715      if(refmultTruth<=0 || refmultReco<=0) return;
1716      cent_v0=refmultTruth;
1717    }
1718  else {
1719  cent_v0=GetAcceptedEventMultiplicity(aod,kTRUE); //centrality value; 2nd argument has no meaning
1720  if(cent_v0<0.) return;
1721  }
1722
1723  effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1724
1725   //get the event plane in case of PbPb
1726    if(fRequestEventPlane){
1727    gReactionPlane=GetEventPlane(aod,kTRUE,cent_v0);//get the truth event plane
1728    if(gReactionPlane==999.) return;
1729  }
1730
1731
1732    Double_t nooftrackstruth=0.0;//in case of pp this will give the multiplicity(for truth case) after the track loop(only for unidentified particles that pass  kinematic cuts)
1733
1734    //TObjArray* tracksMCtruth=0;
1735 TObjArray* tracksMCtruth=new TObjArray;//for truth MC particles with PID,here unidentified means any particle other than pion, kaon or proton(Basicaly Spundefined of AliHelperPID)******WARNING::different from data and reco MC
1736  tracksMCtruth->SetOwner(kTRUE);  //***********************************IMPORTANT!
1737
1738   eventno++;
1739
1740   //There is a small difference between zvtx and zVtxmc??????
1741   //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1742   //cout<<"***********************************************zvtx="<<zvtx<<endl;
1743
1744 //now process the truth particles(for both efficiency & correlation function)
1745 Int_t nMCTrack = fArrayMC->GetEntriesFast();
1746
1747 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1748 {      //MC truth track loop starts
1749
1750 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1751
1752     if(!partMC){
1753       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1754       continue;
1755     }
1756
1757 //consider only charged particles
1758     if(partMC->Charge() == 0) continue;
1759
1760 //consider only primary particles; neglect all secondary particles including from weak decays
1761  if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1762
1763
1764 //remove injected signals(primaries above <maxLabel>)
1765  if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1766
1767 //remove duplicates
1768   Bool_t isduplicate=kFALSE;
1769  if (fRemoveDuplicates)
1770    {
1771  for (Int_t j=iMC+1; j<nMCTrack; ++j)
1772    {//2nd trutuh loop starts
1773 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1774    if(!partMC2){
1775       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1776       continue;
1777     }
1778  if (partMC->GetLabel() == partMC2->GetLabel())
1779    {
1780 isduplicate=kTRUE;
1781  break;
1782    }
1783    }//2nd truth loop ends
1784    }
1785  if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1786
1787 //give only kinematic cuts at the truth level
1788  if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1789  if (partMC->Pt() < fminPt ||  partMC->Pt() > fmaxPt) continue;
1790
1791  if(!partMC) continue;//for safety
1792
1793  //To determine multiplicity in case of PP
1794  nooftrackstruth++;
1795  //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1796 //only physical primary(all/unidentified)
1797 if(ffillhistQATruth)
1798     {
1799  MCtruthpt->Fill(partMC->Pt());
1800  MCtrutheta->Fill(partMC->Eta());
1801  MCtruthphi->Fill(partMC->Phi());
1802     }
1803  //get particle ID
1804 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1805 Int_t particletypeTruth=-999;
1806  if (TMath::Abs(pdgtruth)==211)
1807    {
1808  particletypeTruth=SpPion;
1809 if(ffillhistQATruth)
1810     {
1811  MCtruthpionpt->Fill(partMC->Pt());
1812  MCtruthpioneta->Fill(partMC->Eta());
1813  MCtruthpionphi->Fill(partMC->Phi());
1814     }
1815       }
1816  if (TMath::Abs(pdgtruth)==321)
1817    {
1818  particletypeTruth=SpKaon;
1819 if(ffillhistQATruth)
1820     {
1821  MCtruthkaonpt->Fill(partMC->Pt());
1822  MCtruthkaoneta->Fill(partMC->Eta());
1823  MCtruthkaonphi->Fill(partMC->Phi());
1824   }
1825     }
1826 if(TMath::Abs(pdgtruth)==2212)
1827   {
1828  particletypeTruth=SpProton;
1829 if(ffillhistQATruth)
1830     {
1831  MCtruthprotonpt->Fill(partMC->Pt());
1832  MCtruthprotoneta->Fill(partMC->Eta());
1833  MCtruthprotonphi->Fill(partMC->Phi());
1834     }
1835      }
1836  if(TMath::Abs(pdgtruth)!=211 && TMath::Abs(pdgtruth)!=321 && TMath::Abs(pdgtruth)!=2212)  particletypeTruth=unidentified;//*********************WARNING:: situation is different from reco MC and data case(here we don't have SpUndefined particles,because here unidentified=SpUndefined)
1837
1838     if(fRequestEventPlane){
1839       FillPIDEventPlane(cent_v0,particletypeTruth,partMC->Phi(),gReactionPlane);
1840     }
1841
1842  // -- Fill THnSparse for efficiency and contamination calculation
1843  if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) effcent=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important
1844
1845  Double_t primmctruth[4] = {effcent, zVtxmc,partMC->Pt(), partMC->Eta()};
1846  if(ffillefficiency)
1847   {
1848     fTrackHistEfficiency[5]->Fill(primmctruth,0);//for all primary truth particles(4)
1849     if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[3]->Fill(primmctruth,0);//for  primary truth mesons(3)
1850     if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTrackHistEfficiency[4]->Fill(primmctruth,0);//for  primary truth kaons+protons(4)
1851     if (TMath::Abs(pdgtruth)==211)  fTrackHistEfficiency[0]->Fill(primmctruth,0);//for pions
1852     if (TMath::Abs(pdgtruth)==321)  fTrackHistEfficiency[1]->Fill(primmctruth,0);//for kaons
1853     if(TMath::Abs(pdgtruth)==2212)  fTrackHistEfficiency[2]->Fill(primmctruth,0);//for protons
1854   }
1855
1856  Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1857 if((partMC->Pt()>=fminPtAsso && partMC->Pt()<=fmaxPtAsso) || (partMC->Pt()>=fminPtTrig && partMC->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
1858   {
1859     Short_t chargeval=0;
1860     if(partMC->Charge()>0)   chargeval=1;
1861     if(partMC->Charge()<0)   chargeval=-1;
1862     if(chargeval==0) continue;
1863 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,chargeval,partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
1864 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1865  copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1866  tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1867   }
1868   }//MC truth track loop ends
1869
1870 //*********************still in event loop
1871
1872  if (fSampleType=="PbPb"){
1873    if (fRandomizeReactionPlane)//only for TRuth MC??
1874   {
1875     Double_t centralityDigits = cent_v0*1000. - (Int_t)(cent_v0*1000.);
1876     Double_t angle = TMath::TwoPi() * centralityDigits;
1877     AliInfo(Form("Shifting phi of all tracks by %f (digits %f)", angle, centralityDigits));
1878     ShiftTracks(tracksMCtruth, angle);
1879   }
1880  }
1881
1882  Float_t weghtval=1.0;
1883
1884 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1885   {
1886  //Fill Correlations for MC truth particles(same event)
1887 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1888   Fillcorrelation(gReactionPlane,tracksMCtruth,0,cent_v0,zVtxmc,weghtval,kFALSE,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1889
1890 //start mixing
1891  AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0, zVtxmc+200, gReactionPlane);
1892 if (pool2 && pool2->IsReady())
1893   {//start mixing only when pool->IsReady
1894 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1895   {//proceed only when no. of trigger particles >0 in current event
1896     Float_t nmix=(Float_t)pool2->GetCurrentNEvents();
1897 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1898   { //pool event loop start
1899  TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1900   if(!bgTracks6) continue;
1901   Fillcorrelation(gReactionPlane,tracksMCtruth,bgTracks6,cent_v0,zVtxmc,nmix,(jMix == 0),bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1902
1903    }// pool event loop ends mixing case
1904  }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1905 } //if pool->IsReady() condition ends mixing case
1906
1907  //still in main event loop
1908
1909  if(tracksMCtruth){
1910 if(pool2)  pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1911  }
1912   }
1913
1914  //still in main event loop
1915
1916 if(tracksMCtruth) delete tracksMCtruth;
1917
1918 //now deal with reco tracks
1919
1920 //detrmine the ref mult in case of Reco(not required if we get centrality info from AliCentrality)
1921  if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) cent_v0=refmultReco;
1922  effcent=cent_v0;// This will be required for efficiency THn filling(specially in case of pp)
1923
1924  if(fRequestEventPlane){
1925    gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);//get the reconstructed event plane
1926     if(gReactionPlane==999.) return;
1927  }
1928
1929    //TObjArray* tracksUNID=0;
1930    TObjArray* tracksUNID = new TObjArray;
1931    tracksUNID->SetOwner(kTRUE);
1932
1933    //TObjArray* tracksID=0;
1934    TObjArray* tracksID = new TObjArray;
1935    tracksID->SetOwner(kTRUE);
1936
1937
1938    Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1939
1940
1941    Double_t trackscount=0.0;
1942 // loop over reconstructed tracks
1943   for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1944 { // reconstructed track loop starts
1945   AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1946   if (!track) continue;
1947  //get the corresponding MC track at the truth level (doing reco matching)
1948   AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1949   if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1950
1951 //remove injected signals
1952  if(fInjectedSignals)
1953    {
1954     AliAODMCParticle* mother = recomatched;
1955
1956       while (!mother->IsPhysicalPrimary())
1957       {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1958         if (mother->GetMother() < 0)
1959         {
1960           mother = 0;
1961           break;
1962         }
1963
1964    mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1965         if (!mother)
1966           break;
1967       }
1968  if (!mother)
1969     {
1970       Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
1971       continue;
1972     }
1973  if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
1974    }//remove injected signals
1975
1976  if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
1977
1978   Bool_t isduplicate2=kFALSE;
1979 if (fRemoveDuplicates)
1980    {
1981   for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
1982     {//2nd loop starts
1983  AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
1984  if (!track2) continue;
1985  AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
1986 if(!recomatched2) continue;
1987
1988 if (track->GetLabel() == track2->GetLabel())
1989    {
1990 isduplicate2=kTRUE;
1991  break;
1992    }
1993     }//2nd loop ends
1994    }
1995  if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
1996
1997   fHistQA[11]->Fill(track->GetTPCNcls());
1998   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
1999
2000  if(tracktype==0) continue;
2001  if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
2002 {
2003   if(!track) continue;//for safety
2004  //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
2005   trackscount++;
2006
2007 //check for eta , phi holes
2008  fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2009  fphiSpectraasso->Fill(track->Phi(),track->Pt());
2010
2011   Int_t particletypeMC=-9999;
2012
2013 //tag all particles as unidentified
2014  particletypeMC=unidentified;
2015
2016  Float_t effmatrix=1.;
2017
2018 // -- Fill THnSparse for efficiency calculation
2019  if (fSampleType=="pp" && fCentralityMethod.EndsWith("_MANUAL")) effcent=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicity i.e. multiplicity intregated out )**************Important
2020  //NOTE:: this will be used for fillinfg THnSparse of efficiency & also to get the the track by track eff. factor on the fly(only in case of pp)
2021
2022  //Clone & Reduce track list(TObjArray) for unidentified particles
2023 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2024   {
2025      Short_t chargeval=0;
2026     if(track->Charge()>0)   chargeval=1;
2027     if(track->Charge()<0)   chargeval=-1;
2028     if(chargeval==0) continue;
2029  if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2030    effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);
2031    LRCParticlePID* copy = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2032    copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
2033    tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2034   }
2035
2036 //get the pdg code of the corresponding truth particle
2037  Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
2038
2039  Double_t allrecomatchedpid[4] = {effcent, zVtxmc,recomatched->Pt(), recomatched->Eta()};
2040  if(ffillefficiency) {
2041 fTrackHistEfficiency[5]->Fill(allrecomatchedpid,2);//for allreco matched
2042  if(TMath::Abs(pdgCode)==211 ||  TMath::Abs(pdgCode)==321)   fTrackHistEfficiency[3]->Fill(allrecomatchedpid,2);//for mesons
2043  if(TMath::Abs(pdgCode)==321 ||  TMath::Abs(pdgCode)==2212)   fTrackHistEfficiency[4]->Fill(allrecomatchedpid,2);//for kaons+protons
2044  if(TMath::Abs(pdgCode)==211)  fTrackHistEfficiency[0]->Fill(allrecomatchedpid,2);//for pions
2045  if(TMath::Abs(pdgCode)==321)  fTrackHistEfficiency[1]->Fill(allrecomatchedpid,2);//for kaons
2046  if(TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,2);//for protons
2047
2048  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary()) {
2049  fTrackHistEfficiency[5]->Fill(allrecomatchedpid,1);//for primreco matched
2050  if(TMath::Abs(pdgCode)==211 ||  TMath::Abs(pdgCode)==321)   fTrackHistEfficiency[3]->Fill(allrecomatchedpid,1);//for mesons
2051  if(TMath::Abs(pdgCode)==321 ||  TMath::Abs(pdgCode)==2212)   fTrackHistEfficiency[4]->Fill(allrecomatchedpid,1);//for kaons+protons
2052  if( TMath::Abs(pdgCode)==211)  fTrackHistEfficiency[0]->Fill(allrecomatchedpid,1);//for pions
2053  if( TMath::Abs(pdgCode)==321)  fTrackHistEfficiency[1]->Fill(allrecomatchedpid,1);//for kaons
2054  if( TMath::Abs(pdgCode)==2212) fTrackHistEfficiency[2]->Fill(allrecomatchedpid,1);//for protons
2055  }
2056  }
2057
2058  //now start the particle identification process:)
2059
2060 Float_t dEdx = track->GetTPCsignal();
2061  fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2062
2063  if(HasTOFPID(track))
2064 {
2065 Double_t beta = GetBeta(track);
2066 fHistoTOFbeta->Fill(track->Pt(), beta);
2067  }
2068
2069 //do track identification(nsigma method)
2070   particletypeMC=GetParticle(track,fFIllPIDQAHistos);//******************************problem is here
2071   cout<<particletypeMC<<endl;
2072 switch(TMath::Abs(pdgCode)){
2073   case 2212:
2074     if(fFIllPIDQAHistos){
2075       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2076         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
2077         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpProton,ipid));
2078         h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
2079       }
2080     }
2081     break;
2082   case 321:
2083     if(fFIllPIDQAHistos){
2084       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2085         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
2086         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpKaon,ipid));
2087         h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
2088       }
2089     }
2090     break;
2091   case 211:
2092     if(fFIllPIDQAHistos){
2093       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
2094         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
2095         TH2F *h=GetHistogram2D(Form("NSigmaMC_%d_%d",SpPion,ipid));
2096         h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
2097       }
2098     }
2099     break;
2100   }
2101
2102
2103 //2-d TPCTOF map(for each Pt interval)
2104   if(HasTOFPID(track)){
2105  fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2106  fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2107  fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2108   }
2109
2110  //Pt, Eta , Phi distribution of the reconstructed identified particles
2111 if(ffillhistQAReco)
2112     {
2113 if (particletypeMC==SpPion)
2114   {
2115     fPionPt->Fill(track->Pt());
2116     fPionEta->Fill(track->Eta());
2117     fPionPhi->Fill(track->Phi());
2118   }
2119 if (particletypeMC==SpKaon)
2120   {
2121     fKaonPt->Fill(track->Pt());
2122     fKaonEta->Fill(track->Eta());
2123     fKaonPhi->Fill(track->Phi());
2124   }
2125 if (particletypeMC==SpProton)
2126   {
2127     fProtonPt->Fill(track->Pt());
2128     fProtonEta->Fill(track->Eta());
2129     fProtonPhi->Fill(track->Phi());
2130   }
2131     }
2132
2133  //for misidentification fraction calculation(do it with tuneonPID)
2134  if(particletypeMC==SpPion )
2135    {
2136      if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
2137      if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
2138      if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
2139      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
2140    }
2141 if(particletypeMC==SpKaon )
2142    {
2143      if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
2144      if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
2145      if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
2146      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
2147    }
2148  if(particletypeMC==SpProton )
2149    {
2150      if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
2151      if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
2152      if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
2153      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
2154    }
2155  if(particletypeMC==SpUndefined )
2156    {
2157      if(TMath::Abs(pdgCode)==211) fUNIDcont->Fill(1.,track->Pt());
2158      if(TMath::Abs(pdgCode)==321) fUNIDcont->Fill(3.,track->Pt());
2159      if(TMath::Abs(pdgCode)==2212) fUNIDcont->Fill(5.,track->Pt());
2160      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fUNIDcont->Fill(7.,track->Pt());
2161    }
2162
2163  if(particletypeMC==SpUndefined) continue;
2164
2165     if(fRequestEventPlane){
2166       FillPIDEventPlane(cent_v0,particletypeMC,track->Phi(),gReactionPlane);
2167     }
2168
2169  //fill tracking efficiency
2170  if(ffillefficiency)
2171    {
2172  if(particletypeMC==SpPion || particletypeMC==SpKaon)
2173    {
2174      if(TMath::Abs(pdgCode)==211 ||  TMath::Abs(pdgCode)==321) {
2175        fTrackHistEfficiency[3]->Fill(allrecomatchedpid,4);//for mesons
2176  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[3]->Fill(allrecomatchedpid,3);//for mesons
2177      }
2178    }
2179  if(particletypeMC==SpKaon || particletypeMC==SpProton)
2180    {
2181      if(TMath::Abs(pdgCode)==321 ||  TMath::Abs(pdgCode)==2212) {
2182        fTrackHistEfficiency[4]->Fill(allrecomatchedpid,4);//for kaons+protons
2183  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[4]->Fill(allrecomatchedpid,3);
2184      }
2185    }
2186  if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211)  {
2187    fTrackHistEfficiency[0]->Fill(allrecomatchedpid,4);//for pions
2188  if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[0]->Fill(allrecomatchedpid,3);
2189  }
2190  if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321) {
2191    fTrackHistEfficiency[1]->Fill(allrecomatchedpid,4);//for kaons
2192 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[1]->Fill(allrecomatchedpid,3);
2193  }
2194  if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212){
2195    fTrackHistEfficiency[2]->Fill(allrecomatchedpid,4);//for protons
2196 if (((AliAODMCParticle*)recomatched)->IsPhysicalPrimary())  fTrackHistEfficiency[2]->Fill(allrecomatchedpid,3);
2197  }
2198    }
2199
2200 if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))//to reduce memory consumption in pool
2201   {
2202     Short_t chargeval=0;
2203     if(track->Charge()>0)   chargeval=1;
2204     if(track->Charge()<0)   chargeval=-1;
2205     if(chargeval==0) continue;
2206 if (fapplyTrigefficiency || fapplyAssoefficiency)
2207     effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
2208     LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2209     copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2210     tracksID->Add(copy1);
2211   }
2212   }// if(tracktype==1) condition structure ands
2213
2214 }//reco track loop ends
2215
2216   //*************************************************************still in event loop
2217
2218
2219 if(trackscount>0.0)
2220   {
2221 //fill the centrality/multiplicity distribution of the selected events
2222  fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2223
2224  if (fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2225
2226 //count selected events having centrality betn 0-100%
2227  fEventCounter->Fill(13);
2228
2229 //***************************************event no. counting
2230 Bool_t isbaryontrig=kFALSE;
2231 Bool_t ismesontrig=kFALSE;
2232 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2233
2234 if(tracksID && tracksID->GetEntriesFast()>0)
2235   {
2236 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2237     {  //trigger loop starts
2238       LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2239       if(!trig) continue;
2240       if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2241       Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2242       if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2243       if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2244     }//trig loop ends
2245  if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2246  if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2247   }
2248
2249  //same event delte-eta, delta-phi plot
2250 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2251   {//same event calculation starts
2252     if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2253     if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID)  Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
2254   }
2255
2256 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2257   {//same event calculation starts
2258     if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID)  Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2259     if(ffilltrigIDassoID)   Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weghtval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2260   }
2261
2262 //still in  main event loop
2263 //start mixing
2264  if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2265   AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2266 if (pool && pool->IsReady())
2267   {//start mixing only when pool->IsReady
2268     Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2269  for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2270   { //pool event loop start
2271  TObjArray* bgTracks = pool->GetEvent(jMix);
2272   if(!bgTracks) continue;
2273   if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2274     Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2275  if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2276    Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2277    }// pool event loop ends mixing case
2278
2279 } //if pool->IsReady() condition ends mixing case
2280  if(tracksUNID) {
2281 if(pool)
2282   pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2283  }
2284  }//mixing with unidentified particles
2285
2286  if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2287   AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2288 if (pool1 && pool1->IsReady())
2289   {//start mixing only when pool->IsReady
2290   Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2291 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2292   { //pool event loop start
2293  TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2294   if(!bgTracks2) continue;
2295 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2296   Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2297 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2298   Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2299
2300    }// pool event loop ends mixing case
2301 } //if pool1->IsReady() condition ends mixing case
2302
2303 if(tracksID) {
2304 if(pool1)
2305   pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2306  }
2307  }//mixing with identified particles
2308
2309   //no. of events analyzed
2310 fEventCounter->Fill(15);
2311   }
2312
2313 if(tracksUNID)  delete tracksUNID;
2314
2315 if(tracksID) delete tracksID;
2316
2317
2318 PostData(1, fOutput);
2319
2320 }
2321 //________________________________________________________________________
2322 void AliTwoParticlePIDCorr::doAODevent()
2323 {
2324
2325   //get AOD
2326   AliVEvent *event = InputEvent();
2327   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
2328   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
2329   if (!aod) {
2330     AliError("Cannot get the AOD event");
2331     return;
2332   }
2333
2334 // count all events
2335   fEventCounter->Fill(1);
2336
2337 if (!fPID) return;//this should be available with each event even if we don't do PID selection
2338
2339     fgPsi2v0a=999.;
2340     fgPsi2v0c=999.;
2341     fgPsi2tpc=999.;
2342     fgPsi3v0a=999.;
2343     fgPsi3v0c=999.;
2344     fgPsi3tpc=999.;
2345     gReactionPlane = 999.;
2346
2347   Double_t cent_v0=   -999.;
2348   Double_t effcent=1.0;
2349   Float_t bSign = 0.;
2350   Double_t trackscount=0;//counts particles passed filterbit cuts and kinematic cuts used in this analysis
2351
2352
2353  bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
2354  Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
2355
2356
2357 // check event cuts and fill event histograms and return the centrality or reference multiplicity value
2358  if((cent_v0 = GetAcceptedEventMultiplicity(aod,kFALSE)) < 0){
2359     return;
2360   }
2361
2362   //get the event plane in case of PbPb
2363     if(fRequestEventPlane){
2364       gReactionPlane = GetEventPlane(aod,kFALSE,cent_v0);
2365     if(gReactionPlane==999.) return;
2366     }
2367
2368    TObjArray*  tracksUNID= new TObjArray;//track info before doing PID
2369    tracksUNID->SetOwner(kTRUE);  // IMPORTANT!
2370
2371    TObjArray* tracksID= new TObjArray;//only pions, kaons,protons i.e. after doing the PID selection
2372    tracksID->SetOwner(kTRUE);  // IMPORTANT!
2373
2374
2375     eventno++;
2376
2377     Bool_t fTrigPtmin1=kFALSE;
2378     Bool_t fTrigPtmin2=kFALSE;
2379     Bool_t fTrigPtJet=kFALSE;
2380
2381  for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
2382 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
2383   AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
2384   if (!track) continue;
2385   fHistQA[11]->Fill(track->GetTPCNcls());
2386   Int_t particletype=-9999;//required for PID filling
2387   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kTRUE);//dcacut=kFALSE,onlyprimary=kFALSE
2388   if(tracktype!=1) continue;
2389
2390   if(!track) continue;//for safety
2391
2392 //check for eta , phi holes
2393  fEtaSpectrasso->Fill(track->Eta(),track->Pt());
2394  fphiSpectraasso->Fill(track->Phi(),track->Pt());
2395
2396  trackscount++;
2397
2398  //if no applyefficiency , set the eff factor=1.0
2399  Float_t effmatrix=1.0;
2400
2401  //tag all particles as unidentified that passed filterbit & kinematic cuts
2402  particletype=unidentified;
2403
2404  //To count the no. of tracks having an accepted track in a certain PT(e.g. Jet Pt) range
2405  if(track->Pt()>=fminPtTrig) fTrigPtmin1=kTRUE;
2406  if(track->Pt()>=(fminPtTrig+0.5)) fTrigPtmin2=kTRUE;
2407  if(track->Pt()>=fmaxPtTrig) fTrigPtJet=kTRUE;
2408
2409
2410  if (fSampleType=="pp") effcent=15.0;//integrated over multiplicity [i.e each track has multiplicity 15.0](so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important for efficiency related issues
2411
2412
2413  //to reduce memory consumption in pool
2414   if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2415   {
2416  //Clone & Reduce track list(TObjArray) for unidentified particles
2417     Short_t chargeval=0;
2418     if(track->Charge()>0)   chargeval=1;
2419     if(track->Charge()<0)   chargeval=-1;
2420     if(chargeval==0) continue;
2421  if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
2422    effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);
2423  LRCParticlePID* copy = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2424   copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2425   tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
2426   }
2427
2428 //now start the particle identificaion process:)
2429
2430 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
2431
2432   Float_t dEdx = track->GetTPCsignal();
2433   fHistoTPCdEdx->Fill(track->Pt(), dEdx);
2434
2435   //fill beta vs Pt plots only for tracks having proper TOF response(much less tracks compared to the no. that pass the filterbit & kinematic cuts)
2436  if(HasTOFPID(track))
2437 {
2438   Double_t beta = GetBeta(track);
2439   fHistoTOFbeta->Fill(track->Pt(), beta);
2440  }
2441
2442
2443 //track identification(using nsigma method)
2444      particletype=GetParticle(track,fFIllPIDQAHistos);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
2445
2446
2447 //2-d TPCTOF map(for each Pt interval)
2448   if(HasTOFPID(track)){
2449  fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
2450  fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
2451  fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
2452   }
2453
2454 //ignore the Spundefined particles as they also contain pion, kaon, proton outside the nsigma cut(also if tracks don't have proper TOF PID in a certain Pt interval) & these tracks are actually counted when we do the the efficiency correction, so considering them as unidentified particles & doing the efficiency correction(i.e defining unidentified=pion+Kaon+proton+SpUndefined is right only without efficiency correction) for them will be two times wrong!!!!!
2455   if (particletype==SpUndefined) continue;//this condition creating a modulated structure in delphi projection in mixed event case(only when we are dealing with identified particles i.e. tracksID)!!!!!!!!!!!
2456
2457     if(fRequestEventPlane){
2458       FillPIDEventPlane(cent_v0,particletype,track->Phi(),gReactionPlane);
2459     }
2460
2461
2462  //Pt, Eta , Phi distribution of the reconstructed identified particles
2463 if(ffillhistQAReco)
2464     {
2465 if (particletype==SpPion)
2466   {
2467     fPionPt->Fill(track->Pt());
2468     fPionEta->Fill(track->Eta());
2469     fPionPhi->Fill(track->Phi());
2470   }
2471 if (particletype==SpKaon)
2472   {
2473     fKaonPt->Fill(track->Pt());
2474     fKaonEta->Fill(track->Eta());
2475     fKaonPhi->Fill(track->Phi());
2476   }
2477 if (particletype==SpProton)
2478   {
2479     fProtonPt->Fill(track->Pt());
2480     fProtonEta->Fill(track->Eta());
2481     fProtonPhi->Fill(track->Phi());
2482   }
2483     }
2484
2485  if((track->Pt()>=fminPtAsso && track->Pt()<=fmaxPtAsso) || (track->Pt()>=fminPtTrig && track->Pt()<=fmaxPtTrig))
2486   {
2487     Short_t chargeval=0;
2488     if(track->Charge()>0)   chargeval=1;
2489     if(track->Charge()<0)   chargeval=-1;
2490     if(chargeval==0) continue;
2491 if (fapplyTrigefficiency || fapplyAssoefficiency)
2492   effmatrix=GetTrackbyTrackeffvalue(track,effcent,zvtx,particletype);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles; Bool_t mesoneffrequired=kFALSE
2493  LRCParticlePID* copy1 = new LRCParticlePID(particletype,chargeval,track->Pt(),track->Eta(), track->Phi(),effmatrix);
2494     copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
2495     tracksID->Add(copy1);
2496   }
2497 } //track loop ends but still in event loop
2498
2499 if(trackscount<1.0){
2500   if(tracksUNID) delete tracksUNID;
2501   if(tracksID) delete tracksID;
2502   return;
2503  }
2504
2505  if (fTrigPtmin1) fhistJetTrigestimate->Fill(cent_v0,0.0);
2506  if (fTrigPtmin2) fhistJetTrigestimate->Fill(cent_v0,2.0);
2507  if (fTrigPtJet) fhistJetTrigestimate->Fill(cent_v0,4.0);
2508
2509  Float_t weightval=1.0;
2510
2511
2512
2513 //fill the centrality/multiplicity distribution of the selected events
2514  fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
2515
2516 if(fSampleType=="pPb" || fSampleType=="PbPb") fCentralityCorrelation->Fill(cent_v0, trackscount);//only with unidentified tracks(i.e before PID selection);;;;;can be used to remove centrality outliers??????
2517
2518 //count selected events having centrality betn 0-100%
2519  fEventCounter->Fill(13);
2520
2521 //***************************************event no. counting
2522 Bool_t isbaryontrig=kFALSE;
2523 Bool_t ismesontrig=kFALSE;
2524 if(tracksUNID && tracksUNID->GetEntriesFast()>0) fEventno->Fill(cent_v0,zvtx);
2525
2526 if(tracksID && tracksID->GetEntriesFast()>0)
2527   {
2528 for(Int_t i=0;i<tracksID->GetEntriesFast();i++)
2529     {  //trigger loop starts
2530       LRCParticlePID *trig=(LRCParticlePID*)(tracksID->UncheckedAt(i));
2531       if(!trig) continue;
2532       if(trig->Pt()<fminPtTrig || trig->Pt()>fmaxPtTrig) continue;
2533       Int_t particlepidtrig=trig->getparticle(); //either 1 or 2
2534       if(particlepidtrig==SpProton) isbaryontrig=kTRUE;
2535       if(particlepidtrig==SpPion) ismesontrig=kTRUE;
2536     }//trig loop ends
2537  if (isbaryontrig) fEventnobaryon->Fill(cent_v0,zvtx);
2538  if (ismesontrig) fEventnomeson->Fill(cent_v0,zvtx);
2539   }
2540
2541
2542 //same event delta-eta-deltaphi plot
2543
2544 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
2545   {//same event calculation starts
2546     if(ffilltrigassoUNID) Fillcorrelation(gReactionPlane,tracksUNID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
2547     if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID)  Fillcorrelation(gReactionPlane,tracksUNID,tracksID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
2548   }
2549
2550 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
2551   {//same event calculation starts
2552     if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID)  Fillcorrelation(gReactionPlane,tracksID,tracksUNID,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
2553     if(ffilltrigIDassoID)   Fillcorrelation(gReactionPlane,tracksID,0,cent_v0,zvtx,weightval,kFALSE,bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
2554   }
2555
2556 //still in  main event loop
2557
2558
2559 //start mixing
2560  if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
2561 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx,gReactionPlane);//In the pool there is tracksUNID(i.e associateds are unidentified)
2562 if (pool && pool->IsReady())
2563   {//start mixing only when pool->IsReady
2564   Float_t nmix1=(Float_t)pool->GetCurrentNEvents();
2565  for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
2566   { //pool event loop start
2567  TObjArray* bgTracks = pool->GetEvent(jMix);
2568   if(!bgTracks) continue;
2569   if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
2570     Fillcorrelation(gReactionPlane,tracksUNID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigassoUNID");//mixcase=kTRUE
2571  if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
2572    Fillcorrelation(gReactionPlane,tracksID,bgTracks,cent_v0,zvtx,nmix1,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
2573    }// pool event loop ends mixing case
2574 } //if pool->IsReady() condition ends mixing case
2575  if(tracksUNID) {
2576 if(pool)
2577   pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
2578  }
2579  }//mixing with unidentified particles
2580
2581
2582  if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
2583  AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100,gReactionPlane);//In the pool1 there is tracksID(i.e associateds are identified)
2584 if (pool1 && pool1->IsReady())
2585   {//start mixing only when pool->IsReady
2586   Float_t nmix2=(Float_t)pool1->GetCurrentNEvents();
2587 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
2588   { //pool event loop start
2589  TObjArray* bgTracks2 = pool1->GetEvent(jMix);
2590   if(!bgTracks2) continue;
2591 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
2592   Fillcorrelation(gReactionPlane,tracksUNID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
2593 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
2594   Fillcorrelation(gReactionPlane,tracksID,bgTracks2,cent_v0,zvtx,nmix2,(jMix == 0),bSign,fPtOrderDataReco,ftwoTrackEfficiencyCutDataReco,kTRUE,"trigIDassoID");//mixcase=kTRUE
2595
2596    }// pool event loop ends mixing case
2597 } //if pool1->IsReady() condition ends mixing case
2598
2599 if(tracksID) {
2600 if(pool1)
2601   pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
2602  }
2603  }//mixing with identified particles
2604
2605
2606   //no. of events analyzed
2607 fEventCounter->Fill(15);
2608
2609 //still in main event loop
2610
2611
2612 if(tracksUNID)  delete tracksUNID;
2613
2614 if(tracksID) delete tracksID;
2615
2616
2617 PostData(1, fOutput);
2618
2619 } // *************************event loop ends******************************************//_______________________________________________________________________
2620 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
2621 {
2622   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
2623
2624   TObjArray* tracksClone = new TObjArray;
2625   tracksClone->SetOwner(kTRUE);
2626
2627   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
2628   {
2629     LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
2630     LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());
2631     copy100->SetUniqueID(particle->GetUniqueID());
2632     tracksClone->Add(copy100);
2633   }
2634
2635   return tracksClone;
2636 }
2637
2638 //--------------------------------------------------------------------------------
2639 void AliTwoParticlePIDCorr::Fillcorrelation(Float_t ReactionPlane,TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t weight,Bool_t firstTime,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
2640 {
2641
2642   //before calling this function check that either trackstrig & tracksasso are available
2643
2644  // Eta() is extremely time consuming, therefore cache it for the inner loop here:
2645   TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
2646   TArrayF eta(input->GetEntriesFast());
2647   for (Int_t i=0; i<input->GetEntriesFast(); i++)
2648     eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
2649
2650   //if(trackstrig)
2651     Int_t jmax=trackstrig->GetEntriesFast();
2652   if(tracksasso)
2653      jmax=tracksasso->GetEntriesFast();
2654
2655 // identify K, Lambda candidates and flag those particles
2656     // a TObject bit is used for this
2657 const UInt_t kResonanceDaughterFlag = 1 << 14;
2658     if (fRejectResonanceDaughters > 0)
2659     {
2660       Double_t resonanceMass = -1;
2661       Double_t massDaughter1 = -1;
2662       Double_t massDaughter2 = -1;
2663       const Double_t interval = 0.02;
2664  switch (fRejectResonanceDaughters)
2665       {
2666         case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
2667         case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
2668         case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
2669         default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
2670       }
2671
2672 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2673         trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2674  for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
2675           tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
2676
2677  for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
2678       {
2679       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2680 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
2681         {
2682         LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2683
2684  // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2685 if (trig->IsEqual(asso)) continue;
2686
2687 if (trig->Charge() * asso->Charge() > 0) continue;
2688
2689  Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2690
2691 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
2692           {
2693             mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
2694
2695             if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
2696             {
2697               trig->SetBit(kResonanceDaughterFlag);
2698               asso->SetBit(kResonanceDaughterFlag);
2699
2700 //            Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
2701             }
2702           }
2703         }
2704       }
2705     }
2706
2707       //Select the highest Pt trigger particle in an event (within a given Pt trigger range)
2708
2709     Float_t TriggerPtMin=fminPtTrig;
2710     Float_t TriggerPtMax=fmaxPtTrig;
2711     Int_t HighestPtTriggerIndx=-99999;
2712     TH1* triggerWeighting = 0;
2713
2714 if(fSelectHighestPtTrig || fWeightPerEvent)//**************add this data member to the constructor
2715       {
2716 if (fWeightPerEvent)
2717     {
2718       TAxis* axis=0;
2719    if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID) axis = fTHnTrigcount->GetGrid(0)->GetGrid()->GetAxis(2);
2720   if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH)    axis = fTHnTrigcountMCTruthPrim->GetGrid(0)->GetGrid()->GetAxis(2);
2721       triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
2722     }
2723 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2724     {  //trigger loop starts(highest Pt trigger selection)
2725       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2726       if(!trig) continue;
2727       Float_t trigpt=trig->Pt();
2728     //to avoid overflow qnd underflow
2729       if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2730       Int_t particlepidtrig=trig->getparticle();
2731       if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2732
2733       Float_t trigeta=trig->Eta();
2734
2735       // some optimization
2736  if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2737         continue;
2738
2739 if (fOnlyOneEtaSide != 0)
2740       {
2741         if (fOnlyOneEtaSide * trigeta < 0)
2742           continue;
2743       }
2744   if (fTriggerSelectCharge != 0)
2745         if (trig->Charge() * fTriggerSelectCharge < 0)
2746           continue;
2747
2748       if (fRejectResonanceDaughters > 0)
2749         if (trig->TestBit(kResonanceDaughterFlag)) continue;
2750
2751       if(fSelectHighestPtTrig){
2752  if(trigpt>TriggerPtMin && trigpt<=TriggerPtMax)
2753           {
2754           HighestPtTriggerIndx=(Int_t)trig->GetUniqueID();
2755           TriggerPtMin=trigpt;
2756           }
2757       }
2758
2759 if (fWeightPerEvent)  triggerWeighting->Fill(trigpt);
2760
2761     }//trigger loop ends(highest Pt trigger selection)
2762
2763       }//******************(fSelectHighestPtTrig || fWeightPerEvent) condition ends
2764
2765
2766  //two particle correlation filling
2767 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
2768     {  //trigger loop starts
2769       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
2770       if(!trig) continue;
2771       Float_t trigpt=trig->Pt();
2772     //to avoid overflow qnd underflow
2773       if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
2774       Int_t particlepidtrig=trig->getparticle();
2775       if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
2776
2777       Float_t trigeta=trig->Eta();
2778
2779       // some optimization
2780  if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
2781         continue;
2782
2783 if (fOnlyOneEtaSide != 0)
2784       {
2785         if (fOnlyOneEtaSide * trigeta < 0)
2786           continue;
2787       }
2788   if (fTriggerSelectCharge != 0)
2789         if (trig->Charge() * fTriggerSelectCharge < 0)
2790           continue;
2791
2792       if (fRejectResonanceDaughters > 0)
2793         if (trig->TestBit(kResonanceDaughterFlag)) continue;
2794
2795       if(fSelectHighestPtTrig && HighestPtTriggerIndx!=-99999) {
2796         if(trig->GetUniqueID()!=(UInt_t)HighestPtTriggerIndx) continue;
2797       }
2798
2799       Float_t trigphi=trig->Phi();
2800       Float_t trackefftrig=1.0;
2801       if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
2802
2803     // Event plane (determine psi bin)
2804     Double_t gPsiMinusPhi    =   0.;
2805     Double_t gPsiMinusPhiBin = -10.;
2806 if(fRequestEventPlane){
2807     gPsiMinusPhi   = TMath::Abs(trigphi - ReactionPlane);
2808     //in-plane(Note thet event plane angle has to be defined within 0 to 180 degree(do not use this if event ), otherwise the definition of in plane and out plane particles is wrong)
2809     if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2810       (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
2811        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2812       gPsiMinusPhiBin = 0.0;
2813     /*
2814  if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
2815        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
2816       gPsiMinusPhiBin = 0.0;
2817     */
2818     //intermediate
2819     else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
2820             ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
2821             ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
2822             ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
2823       gPsiMinusPhiBin = 1.0;
2824     //out of plane
2825     else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
2826             ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
2827       gPsiMinusPhiBin = 2.0;
2828     //everything else
2829     else
2830       gPsiMinusPhiBin = 3.0;
2831
2832     fHistPsiMinusPhi->Fill(gPsiMinusPhiBin,gPsiMinusPhi);
2833  }
2834
2835       //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
2836         Double_t* trigval;
2837         Int_t dim=3;
2838         Int_t eventplaneAxis=0;
2839         if(fRequestEventPlane) eventplaneAxis=1;
2840         if(fcontainPIDtrig && SetChargeAxis==0) dim=4+eventplaneAxis;
2841         if(!fcontainPIDtrig && SetChargeAxis==2) dim=4+eventplaneAxis;
2842         if(fcontainPIDtrig && SetChargeAxis==2) dim=5+eventplaneAxis;
2843         trigval= new Double_t[dim];
2844       trigval[0] = cent;
2845       trigval[1] = vtx;
2846       trigval[2] = trigpt;
2847
2848       if(fRequestEventPlane){
2849       trigval[3] = gPsiMinusPhiBin;
2850       if(fcontainPIDtrig && SetChargeAxis==0) trigval[4] = particlepidtrig;
2851       if(!fcontainPIDtrig && SetChargeAxis==2) trigval[4] = trig->Charge();
2852       if(fcontainPIDtrig && SetChargeAxis==2) {
2853       trigval[4] = particlepidtrig;
2854       trigval[5] = trig->Charge();
2855        }
2856       }
2857
2858   if(!fRequestEventPlane){
2859       if(fcontainPIDtrig && SetChargeAxis==0) trigval[3] = particlepidtrig;
2860       if(!fcontainPIDtrig && SetChargeAxis==2) trigval[3] = trig->Charge();
2861       if(fcontainPIDtrig && SetChargeAxis==2) {
2862       trigval[3] = particlepidtrig;
2863       trigval[4] = trig->Charge();
2864        }
2865       }
2866
2867
2868
2869         if (fWeightPerEvent)
2870         {
2871           // leads effectively to a filling of one entry per filled trigger particle pT bin
2872           Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(trigval[2]);
2873 //        Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
2874           trackefftrig *= triggerWeighting->GetBinContent(weightBin);
2875         }
2876
2877
2878       //trigger particle counting for both same and mixed event case;;;;;step=0->same event case;;;;;step=1->mixed event case;;;;;;
2879 if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
2880       if(fillup=="trigassoUNID" ) {
2881 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2882 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2883       }
2884     }
2885  if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
2886    if(fillup=="trigassoUNID" )
2887      {
2888 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2889 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2890      }
2891     }
2892 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
2893   if(fillup=="trigUNIDassoID")
2894     {
2895 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2896 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2897     }
2898     }
2899  //ensure that trigIDassoID , trigassoUNID, trigIDassoUNID & trigUNIDassoID  case FillCorrelation called only once in the event loop for same event correlation function calculation, otherwise there will be multiple counting of pion, kaon,proton,unidentified
2900 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
2901   if(fillup=="trigIDassoID")
2902     {
2903 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2904 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2905     }
2906     }
2907  if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
2908    if(fillup=="trigIDassoUNID" )
2909      {
2910 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2911 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2912      }
2913     }
2914 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
2915   if(fillup=="trigIDassoID")
2916     {
2917 if(mixcase==kFALSE)   fTHnTrigcount->Fill(trigval,0,1.0/trackefftrig);
2918 if(mixcase==kTRUE && firstTime)   fTHnTrigcount->Fill(trigval,1,1.0/trackefftrig);
2919     }
2920     }
2921
2922  if(fillup=="trigIDassoIDMCTRUTH") { //In truth MC case "Unidentified" means any particle other than pion,kaon or proton and no efficiency correction(default value 1.0)************************be careful!!!!
2923 if(mixcase==kFALSE)   fTHnTrigcountMCTruthPrim->Fill(trigval,0,1.0/trackefftrig);
2924 if(mixcase==kTRUE && firstTime)   fTHnTrigcountMCTruthPrim->Fill(trigval,1,1.0/trackefftrig);
2925   }
2926
2927     //asso loop starts within trigger loop
2928    for(Int_t j=0;j<jmax;j++)
2929              {
2930     LRCParticlePID *asso=0;
2931     if(!tracksasso)
2932     asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
2933     else
2934     asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
2935
2936     if(!asso) continue;
2937
2938     //to avoid overflow and underflow
2939  if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
2940
2941  //if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
2942
2943   if(!tracksasso && j==i) continue;
2944
2945    // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event,i.e. both Trig and asso TObjArray belongs to the same Pt range but say Trig is Unidentified but asso is identified then the serial no. wise particles are not same and and j==i doesn't aplly)
2946    // if (tracksasso && trig->IsEqual(asso))  continue;
2947
2948   if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
2949
2950  if (fPtOrder)
2951  if (asso->Pt() >= trig->Pt()) continue;
2952
2953   Int_t particlepidasso=asso->getparticle();
2954   if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
2955
2956
2957 if (fAssociatedSelectCharge != 0)
2958 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
2959
2960  if (fSelectCharge > 0)
2961         {
2962           // skip like sign
2963           if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
2964             continue;
2965
2966           // skip unlike sign
2967           if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
2968             continue;
2969         }
2970
2971 if (fEtaOrdering)
2972         {
2973           if (trigeta < 0 && asso->Eta() < trigeta)
2974             continue;
2975           if (trigeta > 0 && asso->Eta() > trigeta)
2976             continue;
2977         }
2978
2979 if (fRejectResonanceDaughters > 0)
2980           if (asso->TestBit(kResonanceDaughterFlag))
2981           {
2982 //          Printf("Skipped j=%d", j);
2983             continue;
2984           }
2985
2986         // conversions
2987         if (fCutConversions && asso->Charge() * trig->Charge() < 0)
2988         {
2989           Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2990
2991           if (mass < 0.1)
2992           {
2993             mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2994
2995             fControlConvResoncances->Fill(0.0, mass);
2996
2997             if (mass < 0.04*0.04)
2998               continue;
2999           }
3000         }
3001
3002         // K0s
3003         if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3004         {
3005           Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
3006
3007           const Float_t kK0smass = 0.4976;
3008
3009           if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
3010           {
3011             mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
3012
3013             fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
3014
3015             if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
3016               continue;
3017           }
3018         }
3019
3020         // Lambda
3021         if (fCutResonances && asso->Charge() * trig->Charge() < 0)
3022         {
3023           Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
3024           Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3025
3026           const Float_t kLambdaMass = 1.115;
3027
3028           if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
3029           {
3030             mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
3031
3032             fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
3033
3034             if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3035               continue;
3036           }
3037           if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
3038           {
3039             mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
3040
3041             fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
3042
3043             if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
3044               continue;
3045           }
3046         }
3047
3048         if (twoTrackEfficiencyCut)
3049         {
3050           // the variables & cuthave been developed by the HBT group
3051           // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
3052           Float_t phi1 = trig->Phi();
3053           Float_t pt1 = trig->Pt();
3054           Float_t charge1 = trig->Charge();
3055           Float_t phi2 = asso->Phi();
3056           Float_t pt2 = asso->Pt();
3057           Float_t charge2 = asso->Charge();
3058
3059           Float_t deta= trigeta - eta[j];
3060
3061  // optimization
3062           if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
3063           {
3064
3065   // check first boundaries to see if is worth to loop and find the minimum
3066             Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
3067             Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
3068
3069  const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
3070
3071             Float_t dphistarminabs = 1e5;
3072             Float_t dphistarmin = 1e5;
3073
3074  if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
3075             {
3076               for (Double_t rad=0.8; rad<2.51; rad+=0.01)
3077               {
3078                 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
3079
3080                 Float_t dphistarabs = TMath::Abs(dphistar);
3081
3082         if (dphistarabs < dphistarminabs)
3083                 {
3084                   dphistarmin = dphistar;
3085                   dphistarminabs = dphistarabs;
3086                 }
3087               }
3088               if(mixcase==kFALSE)  fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3089               if(mixcase==kTRUE)  fTwoTrackDistancePtmix[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3090
3091 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
3092               {
3093 //              Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
3094                 continue;
3095               }
3096              if(mixcase==kFALSE) fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for same event
3097              if(mixcase==kTRUE) fTwoTrackDistancePtmix[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));//for mixed event
3098
3099             }
3100           }
3101         }
3102
3103         Float_t weightperevent=weight;
3104         Float_t trackeffasso=1.0;
3105         if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
3106         //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
3107         Float_t deleta=trigeta-eta[j];
3108         Float_t delphi=PhiRange(trigphi-asso->Phi());
3109
3110         Float_t delpt=trigpt-asso->Pt();
3111         //fill it with/without two track efficiency cut
3112         if(mixcase==kFALSE)  fTwoTrackDistancePtdip->Fill(deleta, delphi, TMath::Abs(delpt));//for same event
3113         if(mixcase==kTRUE)  fTwoTrackDistancePtdipmix->Fill(deleta, delphi, TMath::Abs(delpt));//for mixed event
3114
3115  //here get the two particle efficiency correction factor
3116         Float_t effweight=trackefftrig*trackeffasso*weightperevent;
3117         // if(mixcase==kFALSE)  cout<<"*******************effweight="<<effweight<<endl;
3118         Double_t* vars;
3119         Int_t dimused=kTrackVariablesPair+eventplaneAxis;
3120         vars= new Double_t[dimused];
3121         vars[0]=cent;
3122         vars[1]=vtx;
3123         vars[2]=trigpt;
3124         vars[3]=asso->Pt();
3125         vars[4]=deleta;
3126         vars[5]=delphi;
3127
3128         Int_t dimension=6;
3129         if(fRequestEventPlane)
3130         {
3131        vars[6]=gPsiMinusPhiBin;
3132        dimension=7;
3133         }
3134
3135 if(!fcontainPIDtrig && !fcontainPIDasso && SetChargeAxis==2){
3136         vars[dimension]=trig->Charge();
3137         vars[dimension+1]=asso->Charge();
3138  }
3139 if(fcontainPIDtrig && !fcontainPIDasso){
3140         vars[dimension]=particlepidtrig;
3141 if(SetChargeAxis==2){
3142         vars[dimension+1]=trig->Charge();
3143         vars[dimension+2]=asso->Charge();
3144  }
3145         }
3146 if(!fcontainPIDtrig && fcontainPIDasso){
3147         vars[dimension]=particlepidasso;
3148 if(SetChargeAxis==2){
3149         vars[dimension+1]=trig->Charge();
3150         vars[dimension+2]=asso->Charge();
3151    }
3152  }
3153  if(fcontainPIDtrig && fcontainPIDasso){
3154         vars[dimension]=particlepidtrig;
3155         vars[dimension+1]=particlepidasso;
3156 if(SetChargeAxis==2){
3157         vars[dimension+2]=trig->Charge();
3158         vars[dimension+3]=asso->Charge();
3159    }
3160  }
3161
3162         if (fWeightPerEvent)
3163         {
3164           Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
3165 //        Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
3166          effweight *= triggerWeighting->GetBinContent(weightBin);
3167         }
3168
3169
3170         //Fill Correlation ThnSparses
3171     if(fillup=="trigassoUNID")
3172       {
3173     if(mixcase==kFALSE)  fTHnCorrUNID->Fill(vars,0,1.0/effweight);
3174     if(mixcase==kTRUE)   fTHnCorrUNIDmix->Fill(vars,0,1.0/effweight);
3175       }
3176     if(fillup=="trigIDassoID")
3177       {
3178         if(mixcase==kFALSE)  fTHnCorrID->Fill(vars,0,1.0/effweight);
3179         if(mixcase==kTRUE)  fTHnCorrIDmix->Fill(vars,0,1.0/effweight);
3180       }
3181     if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
3182       {//in this case effweight should be 1 always
3183         if(mixcase==kFALSE)  fCorrelatonTruthPrimary->Fill(vars,0,1.0/effweight);
3184         if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,0,1.0/effweight);
3185     }
3186     if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
3187       {
3188         if(mixcase==kFALSE)  fTHnCorrIDUNID->Fill(vars,0,1.0/effweight);
3189         if(mixcase==kTRUE)   fTHnCorrIDUNIDmix->Fill(vars,0,1.0/effweight);
3190        }
3191
3192 delete[] vars;
3193    }//asso loop ends
3194 delete[] trigval;
3195  }//trigger loop ends
3196
3197  if (triggerWeighting)
3198     {
3199       delete triggerWeighting;
3200       triggerWeighting = 0;
3201     }
3202 }
3203
3204 //--------------------------------------------------------------------------------
3205 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
3206 {
3207   //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
3208  Int_t effVars[4];
3209  Float_t effvalue=1.;
3210
3211   if(parpid==unidentified)
3212             {
3213             effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
3214             effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
3215             effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
3216             effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
3217             effvalue=effcorection[5]->GetBinContent(effVars);
3218             }
3219 if(parpid==SpPion || parpid==SpKaon)
3220             {
3221               if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3222                 {
3223             effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
3224             effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
3225             effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
3226             effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
3227             effvalue=effcorection[3]->GetBinContent(effVars);
3228                 }
3229               else{
3230  if(parpid==SpPion)
3231             {
3232             effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
3233             effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
3234             effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
3235             effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
3236             effvalue=effcorection[0]->GetBinContent(effVars);
3237             }
3238
3239  if(parpid==SpKaon)
3240             {
3241             effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
3242             effVars[1] =  effcorection[1]->GetAxis(1)->FindBin(evzvtx);
3243             effVars[2] =  effcorection[1]->GetAxis(2)->FindBin(track->Pt());
3244             effVars[3] =  effcorection[1]->GetAxis(3)->FindBin(track->Eta());
3245             effvalue=effcorection[1]->GetBinContent(effVars);
3246             }
3247               }
3248             }
3249
3250  if(parpid==SpProton)
3251             {
3252             effVars[0] =  effcorection[2]->GetAxis(0)->FindBin(cent);
3253             effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
3254             effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
3255             effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
3256             effvalue=effcorection[2]->GetBinContent(effVars);
3257             }
3258
3259  if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
3260                 {
3261   if(parpid==SpProton || parpid==SpKaon)
3262             {
3263             effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
3264             effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
3265             effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
3266             effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
3267             effvalue=effcorection[4]->GetBinContent(effVars);
3268             }
3269                 }
3270             //    Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
3271      if(effvalue==0.) effvalue=1.;
3272
3273      return effvalue;
3274
3275 }
3276 //---------------------------------------------------------------------------------
3277
3278 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield, Bool_t fill)
3279 {
3280
3281   if(!track) return 0;
3282   Bool_t trackOK = track->TestFilterBit(fFilterBit);
3283   if(!trackOK) return 0;
3284   if (fTrackStatus != 0 && !CheckTrack(track)) return 0;
3285   //select only primary traks(for data & reco MC tracks)
3286   if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
3287   if(track->Charge()==0) return 0;
3288   if (fill) fHistQA[12]->Fill(track->GetTPCNcls());
3289   Float_t dxy, dz;
3290   dxy = track->DCA();
3291   dz = track->ZAtDCA();
3292   if (fill) fHistQA[6]->Fill(dxy);
3293   if (fill) fHistQA[7]->Fill(dz);
3294   Float_t chi2ndf = track->Chi2perNDF();
3295   if (fill) fHistQA[13]->Fill(chi2ndf);
3296   // Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
3297   Float_t nCrossedRowsTPC = track->GetTPCNCrossedRows();
3298   if (fill) fHistQA[14]->Fill(nCrossedRowsTPC);
3299   //Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
3300   if (track->GetTPCNclsF()>0) {
3301    Float_t  ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
3302    if (fill) fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
3303     }
3304 //accepted tracks
3305      Float_t pt=track->Pt();
3306      if(pt< fminPt || pt> fmaxPt) return 0;
3307      if(TMath::Abs(track->Eta())> fmaxeta) return 0;
3308      if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
3309      //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required???
3310 // DCA XY
3311         if (fdcacut && fDCAXYCut)
3312         {
3313           if (!vertex)
3314             return 0;
3315
3316           Double_t pos[2];
3317           Double_t covar[3];
3318           AliAODTrack* clone =(AliAODTrack*) track->Clone();
3319           Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
3320           delete clone;
3321           if (!success)
3322             return 0;
3323
3324 //        Printf("%f", ((AliAODTrack*)part)->DCA());
3325 //        Printf("%f", pos[0]);
3326           if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
3327             return 0;
3328         }
3329
3330         if (fSharedClusterCut >= 0)
3331         {
3332           Double_t frac = Double_t(((AliAODTrack*)track)->GetTPCnclsS()) / Double_t(((AliAODTrack*)track)->GetTPCncls());
3333           if (frac > fSharedClusterCut)
3334             return 0;
3335         }
3336
3337    // Rejects tracks with shared clusters after filling a control histogram
3338    // This overload is used for primaries
3339
3340      // Get the shared maps
3341       const TBits sharedMap = track->GetTPCSharedMap();
3342      // Fill a control histogram
3343       fPriHistShare->Fill(sharedMap.CountBits());
3344
3345     // Reject shared clusters
3346        if (fSharedTPCmapCut >= 0)
3347         {
3348       if((sharedMap.CountBits()) >= 1)  return 0;// Bad track, has too many shared clusters!
3349         }
3350
3351      if (fill) fHistQA[8]->Fill(pt);
3352      if (fill) fHistQA[9]->Fill(track->Eta());
3353      if (fill) fHistQA[10]->Fill(track->Phi());
3354      return 1;
3355   }
3356   //________________________________________________________________________________
3357 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track, Bool_t FIllQAHistos)
3358 {
3359 //This function is called within the func GetParticle() for accepted tracks only i.e.after call of Classifytrack() & for those tracks which have proper TPC PID response . combined nsigma(circular) cut only for particles having pt upto  4.0 Gev/c and beyond that use the asymmetric nsigma cut around pion's mean position in TPC ( while filling the  TObjArray for trig & asso )
3360 Float_t pt=track->Pt();
3361
3362 //plot the separation power
3363
3364 Double_t bethe[AliPID::kSPECIES]={0.};
3365 Double_t sigma_TPC[AliPID::kSPECIES]={0.};
3366
3367  Double_t Pi_Ka_sep[NSigmaPIDType+1]={0.};
3368  Double_t Pi_Pr_sep[NSigmaPIDType+1]={0.};
3369  Double_t Ka_Pr_sep[NSigmaPIDType+1]={0.};
3370
3371
3372     Double_t ptpc = track->GetTPCmomentum();
3373     Int_t dEdxN = track->GetTPCsignalN();
3374  for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
3375        bethe[ipart] = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart);
3376       //Double_t diff = dEdx - bethe;
3377        sigma_TPC[ipart] = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart);
3378       //nSigma[ipart] = diff / sigma;
3379     }
3380  Pi_Ka_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kKaon])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kKaon])/2.0);
3381  Pi_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kPion]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kPion]+sigma_TPC[AliPID::kProton])/2.0);
3382  Ka_Pr_sep[NSigmaTPC]=TMath::Abs(bethe[AliPID::kKaon]-bethe[AliPID::kProton])/((sigma_TPC[AliPID::kKaon]+sigma_TPC[AliPID::kProton])/2.0);
3383
3384
3385 Double_t sigma_TOF[AliPID::kSPECIES]={0.};
3386
3387 if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3388    {
3389  Double_t timei[AliPID::kSPECIES];
3390  track->GetIntegratedTimes(timei);
3391  for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {  sigma_TOF[ipart]= fPID->GetTOFResponse().GetExpectedSigma(track->P(), timei[ipart], AliPID::ParticleMass(ipart));}
3392  Pi_Ka_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kKaon])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kKaon])/2.0);
3393  Pi_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kPion]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kPion]+sigma_TOF[AliPID::kProton])/2.0);
3394  Ka_Pr_sep[NSigmaTOF]=TMath::Abs(timei[AliPID::kKaon]-timei[AliPID::kProton])/((sigma_TOF[AliPID::kKaon]+sigma_TOF[AliPID::kProton])/2.0);
3395
3396   Pi_Ka_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Ka_sep[NSigmaTPC]*Pi_Ka_sep[NSigmaTPC]+Pi_Ka_sep[NSigmaTOF]*Pi_Ka_sep[NSigmaTOF]);
3397   Pi_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Pi_Pr_sep[NSigmaTPC]*Pi_Pr_sep[NSigmaTPC]+Pi_Pr_sep[NSigmaTOF]*Pi_Pr_sep[NSigmaTOF]);
3398   Ka_Pr_sep[NSigmaTPCTOF]=TMath::Abs(Ka_Pr_sep[NSigmaTPC]*Ka_Pr_sep[NSigmaTPC]+Ka_Pr_sep[NSigmaTOF]*Ka_Pr_sep[NSigmaTOF]);
3399    }
3400
3401
3402 //fill separation power histograms
3403  for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3404    if(ipid==0){
3405         TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3406         h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3407         TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3408         h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3409         TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3410         h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3411    }
3412    if(HasTOFPID(track) && pt>fPtTOFPIDmin && ipid!=0){
3413        TH2F *h=GetHistogram2D(Form("Pi_Ka_sep_%d",ipid));
3414         h->Fill(track->Pt(),Pi_Ka_sep[ipid]);
3415         TH2F *h1=GetHistogram2D(Form("Pi_Pr_sep_%d",ipid));
3416         h1->Fill(track->Pt(),Pi_Pr_sep[ipid]);
3417         TH2F *h2=GetHistogram2D(Form("Ka_Pr_sep_%d",ipid));
3418         h2->Fill(track->Pt(),Ka_Pr_sep[ipid]);
3419    }
3420  }
3421
3422
3423 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
3424 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
3425 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
3426 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
3427
3428 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
3429 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
3430
3431  if(HasTOFPID(track) && pt>fPtTOFPIDmin)
3432    {
3433
3434 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
3435 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
3436 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
3437 //---combined
3438 nsigmaTPCTOFkPion   = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
3439 nsigmaTPCTOFkKaon   = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
3440 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
3441
3442
3443    }
3444 else{
3445     // --- combined
3446     // if TOF is missing and below fPtTOFPID only the TPC information is used
3447     nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
3448     nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
3449     nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
3450
3451   }
3452
3453 //set data member fnsigmas
3454   fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
3455   fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
3456   fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
3457
3458   //for all tracks below fPtTOFPIDmin  and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
3459   fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
3460   fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
3461   fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
3462
3463  //for all tracks below fPtTOFPIDmin  and also for tracks above fPtTOFPIDmin without proper TOF response these TPCTOF nsigma values will be TMath::Abs(TPC only nsigma)
3464   fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
3465   fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
3466   fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
3467
3468  if(FIllQAHistos){
3469     //Fill NSigma SeparationPlot
3470     for(Int_t ipart=0;ipart<NSpecies;ipart++){
3471       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3472         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3473         TH2F *h=GetHistogram2D(Form("NSigma_%d_%d",ipart,ipid));
3474         h->Fill(track->Pt(),fnsigmas[ipart][ipid]);
3475       }
3476     }
3477   }
3478
3479 }
3480 //----------------------------------------------------------------------------
3481 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track,Bool_t FillQAHistos)
3482 {
3483   //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3484 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && (!HasTOFPID(track)) )return SpUndefined;//so any track having Pt>0.6 withot having proper TOF response will be defined as SpUndefined
3485 //get the identity of the particle with the minimum Nsigma
3486   Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3487   switch (fPIDType){
3488   case NSigmaTPC:
3489     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3490     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPC])  ;
3491     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPC])  ;
3492     break;
3493   case NSigmaTOF:
3494     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3495     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTOF])  ;
3496     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTOF])  ;
3497     break;
3498   case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3499     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3500     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3501     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3502     break;
3503   case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3504     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3505     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3506     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3507     break;
3508   }
3509
3510
3511 if(fdiffPIDcutvalues){
3512   if(track->Pt()<=4) fNSigmaPID=fPIDCutval1;
3513   if(track->Pt()>4 && track->Pt()<=6) fNSigmaPID=fPIDCutval2;
3514   if(track->Pt()>6 && track->Pt()<=8) fNSigmaPID=fPIDCutval3;
3515   if(track->Pt()>8) fNSigmaPID=fPIDCutval4;
3516   }
3517
3518  // guess the particle based on the smaller nsigma (within fNSigmaPID)
3519   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
3520
3521   if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)){
3522     if((fHighPtKaonNSigmaPID>0) && (track->Pt()>fHighPtKaonSigma) && (nsigmaKaon > fHighPtKaonNSigmaPID)) return SpUndefined;//different nsigma cut for kaons above a particular Pt range(within the TPC-TOF PID range)
3523 if(FillQAHistos){
3524       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3525         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3526         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpKaon,ipid));
3527         h->Fill(track->Pt(),fnsigmas[SpKaon][ipid]);
3528       }
3529     }
3530  return SpKaon;
3531   }
3532   if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) {
3533  if(FillQAHistos){
3534       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3535         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3536         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpPion,ipid));
3537         h->Fill(track->Pt(),fnsigmas[SpPion][ipid]);
3538       }
3539     }
3540 return SpPion;
3541   }
3542   if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) {
3543 if(FillQAHistos){
3544       for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3545         if((ipid!=NSigmaTPC) && (!HasTOFPID(track)))continue;//not filling TOF and combined if no TOF PID
3546         TH2F *h=GetHistogram2D(Form("NSigmaRec_%d_%d",SpProton,ipid));
3547         h->Fill(track->Pt(),fnsigmas[SpProton][ipid]);
3548       }
3549  }
3550 return SpProton;
3551   }
3552
3553 // else, return undefined
3554   return SpUndefined;
3555
3556
3557 }
3558
3559 //------------------------------------------------------------------------------------------
3560 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk,Bool_t FIllQAHistos){
3561   //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
3562
3563   //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
3564   //fill DC histos
3565   for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
3566
3567   Int_t MinNSigma=FindMinNSigma(trk,kFALSE);//not filling the NSigmaRec histos
3568
3569
3570   if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
3571
3572   Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
3573   switch (fPIDType) {
3574   case NSigmaTPC:
3575     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
3576     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPC])  ;
3577     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPC])  ;
3578     break;
3579   case NSigmaTOF:
3580     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
3581     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTOF])  ;
3582     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTOF])  ;
3583     break;
3584   case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
3585     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3586     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3587     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3588     break;
3589   case Bayes://the nsigma in the bayesian is used to clean with a very large n-sigma value
3590     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
3591     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
3592     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
3593     break;
3594   }
3595
3596   // Actually the tracks in the overlapping region(in TPC-TOF nSigma plane) will be ignored
3597
3598   if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
3599   if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
3600   if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
3601
3602
3603
3604 if(FIllQAHistos){
3605     //fill NSigma distr for double counting
3606     for(Int_t ipart=0;ipart<NSpecies;ipart++){
3607       if(fHasDoubleCounting[ipart]){//this may be kTRUE only for particles having Pt<=4.0 GeV/C, so this histo contains all the particles having Pt<=4.0 GeV/C in the nsigma overlapping region in TPC/TPC-TOF plane
3608         for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++){
3609           if((ipid!=NSigmaTPC) && (!HasTOFPID(trk)))continue;//not filling TOF and combined if no TOF PID
3610           TH2F *h=GetHistogram2D(Form("NSigmaDC_%d_%d",ipart,ipid));
3611           h->Fill(trk->Pt(),fnsigmas[ipart][ipid]);
3612         }
3613       }
3614     }
3615   }
3616
3617
3618   return fHasDoubleCounting;
3619 }
3620
3621 //////////////////////////////////////////////////////////////////////////////////////////////////
3622
3623 Bool_t* AliTwoParticlePIDCorr::GetAllCompatibleIdentitiesNSigma(AliAODTrack * trk,Bool_t FIllQAHistos){
3624  //mainly intended to check the probability of the PID of the tracks which are in the overlapping nSigma regions and near about the middle position from the   mean position of two ID particle
3625   Bool_t *IDs=GetDoubleCounting(trk,FIllQAHistos);
3626   IDs[FindMinNSigma(trk,FIllQAHistos)]=kTRUE;
3627   return IDs;
3628
3629 }
3630 //////////////////////////////////////////////////////////////////////////////////////////////////
3631
3632 UInt_t AliTwoParticlePIDCorr::CalcPIDCombined(AliAODTrack *track, Int_t detMask, Double_t* prob) const{
3633   //
3634   // Bayesian PID calculation
3635   //
3636   for(Int_t i=0;i<AliPID::kSPECIES;i++)
3637     {
3638       prob[i]=0.;
3639     }
3640   fPIDCombined->SetDetectorMask(detMask);
3641
3642   return fPIDCombined->ComputeProbabilities((AliAODTrack*)track, fPID, prob);
3643 }
3644
3645 //////////////////////////////////////////////////////////////////////////////////////////////////
3646
3647 Int_t AliTwoParticlePIDCorr::GetIDBayes(AliAODTrack * trk, Bool_t FIllQAHistos){
3648
3649   Bool_t *IDs=GetAllCompatibleIdentitiesNSigma(trk,FIllQAHistos);
3650
3651
3652   //Filling of Probability histos
3653         Double_t probTPC[AliPID::kSPECIES]={0.};
3654         Double_t probTOF[AliPID::kSPECIES]={0.};
3655         Double_t probTPCTOF[AliPID::kSPECIES]={0.};
3656
3657         UInt_t detUsedTPC = 0;
3658         UInt_t detUsedTOF = 0;
3659         UInt_t detUsedTPCTOF = 0;
3660
3661  //get the TPC probability
3662           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
3663           detUsedTPC = fPIDCombined->ComputeProbabilities(trk, fPID, probTPC);
3664 if(detUsedTPC == AliPIDResponse::kDetTPC)
3665   {
3666 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3667
3668         TH2F *h=GetHistogram2D(Form("probBayes_TPC_%d",ipart));
3669         if(ipart==0)    h->Fill(trk->Pt(),probTPC[AliPID::kPion]);
3670         if(ipart==1)    h->Fill(trk->Pt(),probTPC[AliPID::kKaon]);
3671         if(ipart==2)    h->Fill(trk->Pt(),probTPC[AliPID::kProton]);
3672  }
3673   }
3674
3675   //get the TOF probability
3676           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
3677           detUsedTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTOF);
3678 if(detUsedTOF == AliPIDResponse::kDetTOF)
3679   {
3680 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3681         TH2F *h=GetHistogram2D(Form("probBayes_TOF_%d",ipart));
3682         if(ipart==0)    h->Fill(trk->Pt(),probTOF[AliPID::kPion]);
3683         if(ipart==1)    h->Fill(trk->Pt(),probTOF[AliPID::kKaon]);
3684         if(ipart==2)    h->Fill(trk->Pt(),probTOF[AliPID::kProton]);
3685  }
3686   }
3687
3688  //get the TPC-TOF probability
3689           fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC);
3690           detUsedTPCTOF = fPIDCombined->ComputeProbabilities(trk, fPID, probTPCTOF);
3691 if(detUsedTPCTOF == (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))
3692   {
3693 for(Int_t ipart=0;ipart<NSpecies;ipart++){
3694         TH2F *h=GetHistogram2D(Form("probBayes_TPCTOF_%d",ipart));
3695         if(ipart==0)    h->Fill(trk->Pt(),probTPCTOF[AliPID::kPion]);
3696         if(ipart==1)    h->Fill(trk->Pt(),probTPCTOF[AliPID::kKaon]);
3697         if(ipart==2)    h->Fill(trk->Pt(),probTPCTOF[AliPID::kProton]);
3698 }
3699   }
3700
3701
3702   Double_t probBayes[AliPID::kSPECIES];
3703
3704   UInt_t detUsed= 0;
3705   if(HasTOFPID(trk) && trk->Pt()>fPtTOFPIDmin){//use TOF information
3706     detUsed = CalcPIDCombined(trk, AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC, probBayes);
3707     if(detUsed != (AliPIDResponse::kDetTOF|AliPIDResponse::kDetTPC))return SpUndefined;//check that TPC and TOF are used
3708   }else{
3709     detUsed = CalcPIDCombined(trk,AliPIDResponse::kDetTPC, probBayes);
3710     if(detUsed != AliPIDResponse::kDetTPC)return SpUndefined;//check that TPC is used
3711   }
3712
3713   //the probability has to be normalized to one, we check it
3714   Double_t sump=0.;
3715   for(Int_t ipart=0;ipart<AliPID::kSPECIES;ipart++)sump+=probBayes[ipart];
3716   if(sump<.99 && sump>1.01){//FIXME precision problem in the sum, workaround
3717     AliFatal("Bayesian probability not normalized to one");
3718   }
3719
3720   if(fdiffPIDcutvalues){
3721   if(trk->Pt()<=4) fBayesCut=fPIDCutval1;
3722   if(trk->Pt()>4 && trk->Pt()<=6) fBayesCut=fPIDCutval2;
3723   if(trk->Pt()>6 && trk->Pt()<=8) fBayesCut=fPIDCutval3;
3724   if(trk->Pt()>8) fBayesCut=fPIDCutval4;
3725   }
3726
3727
3728   //probabilities are normalized to one, if the cut is above .5 there is no problem
3729   if(probBayes[AliPID::kPion]>fBayesCut && IDs[SpPion]==1){
3730     TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpPion));
3731     h->Fill(trk->Pt(),probBayes[AliPID::kPion]);
3732     return SpPion;
3733   }
3734   else if(probBayes[AliPID::kKaon]>fBayesCut && IDs[SpKaon]==1){
3735     TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpKaon));
3736     h->Fill(trk->Pt(),probBayes[AliPID::kKaon]);
3737     return SpKaon;
3738   }
3739   else if(probBayes[AliPID::kProton]>fBayesCut && IDs[SpProton]==1){
3740     TH2F *h=GetHistogram2D(Form("BayesRec_%d",SpProton));
3741     h->Fill(trk->Pt(),probBayes[AliPID::kProton]);
3742     return SpProton;
3743   }
3744   else{
3745     return SpUndefined;
3746   }
3747 }
3748
3749
3750 //////////////////////////////////////////////////////////////////////////////////////////////////
3751 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk, Bool_t FIllQAHistos){
3752   //return the specie according to the minimum nsigma value
3753   //no double counting, this has to be evaluated using CheckDoubleCounting()
3754
3755   Int_t ID=SpUndefined;
3756
3757   CalculateNSigmas(trk,FIllQAHistos);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
3758
3759
3760  //Do PID
3761   if(fPIDType==Bayes){//use bayesianPID
3762
3763     if(!fPIDCombined) {
3764       AliFatal("PIDCombined object has to be set in the steering macro");
3765     }
3766
3767     ID = GetIDBayes(trk,FIllQAHistos);
3768
3769   }else{ //use nsigma PID
3770
3771    ID=FindMinNSigma(trk,FIllQAHistos);
3772 if(fUseExclusiveNSigma){ //if one particle has double counting and exclusive nsigma is requested ID = kSpUndefined
3773       Bool_t *HasDC;
3774       HasDC=GetDoubleCounting(trk,FIllQAHistos);
3775       for(Int_t ipart=0;ipart<NSpecies;ipart++){
3776         if(HasDC[ipart]==kTRUE)  ID = SpUndefined;
3777       }
3778     }
3779   }
3780 //Fill PID signal plot
3781   if(ID != SpUndefined){
3782     for(Int_t idet=0;idet<fNDetectors;idet++){
3783       TH2F *h=GetHistogram2D(Form("PID_%d_%d",idet,ID));
3784       if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3785       if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3786       if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3787     }
3788   }
3789   //Fill PID signal plot without cuts
3790   for(Int_t idet=0;idet<fNDetectors;idet++){
3791     TH2F *h=GetHistogram2D(Form("PIDAll_%d",idet));
3792     if(idet==fITS)h->Fill(trk->P(),trk->GetITSsignal()*trk->Charge());
3793     if(idet==fTPC)h->Fill(trk->P(),trk->GetTPCsignal()*trk->Charge());
3794     if(idet==fTOF && HasTOFPID(trk))h->Fill(trk->P(),GetBeta(trk)*trk->Charge());
3795   }
3796
3797   return ID;
3798 }
3799
3800 //-------------------------------------------------------------------------------------
3801 Bool_t
3802 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
3803 {
3804   // check PID signal
3805    AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
3806    if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
3807    //ULong_t status=track->GetStatus();
3808    //if  (!( (status & AliAODTrack::kTPCpid  ) == AliAODTrack::kTPCpid  )) return kFALSE;//remove light nuclei
3809    //if (track->GetTPCsignal() <= 0.) return kFALSE;
3810    // if(track->GetTPCsignalN() < 60) return kFALSE;//tracks with TPCsignalN< 60 have questionable dEdx,cutting on TPCsignalN > 70 or > 60 shouldn't make too much difference in statistics,also  it is IMO safe to use TPC also for MIPs.
3811
3812   return kTRUE;
3813 }
3814 //___________________________________________________________
3815
3816 Bool_t
3817 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
3818 {
3819   // check TOF matched track
3820   //ULong_t status=track->GetStatus();
3821   //if  (!( (status & AliAODTrack::kITSin  ) == AliAODTrack::kITSin  )) return kFALSE;
3822  AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
3823  if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
3824   if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
3825  //if(!((status & AliAODTrack::kTOFpid  ) == AliAODTrack::kTOFpid  )) return kFALSE;
3826  //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
3827  // if (probMis > 0.01) return kFALSE;
3828 if(fRemoveTracksT0Fill)
3829     {
3830 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
3831       if (startTimeMask < 0)return kFALSE;
3832     }
3833   return kTRUE;
3834 }
3835
3836 //________________________________________________________________________
3837 Double_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
3838 {
3839   //it is called only when TOF PID is available
3840 //TOF beta calculation
3841   Double_t tofTime=track->GetTOFsignal();
3842
3843   Double_t c=TMath::C()*1.E-9;// m/ns
3844   Float_t startTime = fPID->GetTOFResponse().GetStartTime(((AliVTrack*)track)->P());//in ps
3845   Double_t length= fPID->GetTOFResponse().GetExpectedSignal(track,AliPID::kElectron)*1E-3*c;
3846   tofTime -= startTime;      // subtract startTime to the signal
3847   Double_t tof= tofTime*1E-3; // ns, average T0 fill subtracted, no info from T0detector
3848   tof=tof*c;
3849   return length/tof;
3850
3851
3852   /*
3853   Double_t p = track->P();
3854   Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
3855   Double_t timei[5];
3856   track->GetIntegratedTimes(timei);
3857   return timei[0]/time;
3858   */
3859 }
3860 //------------------------------------------------------------------------------------------------------
3861
3862 Float_t AliTwoParticlePIDCorr::GetInvMassSquared(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
3863 {
3864   // calculate inv mass squared
3865   // same can be achieved, but with more computing time with
3866   /*TLorentzVector photon, p1, p2;
3867   p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
3868   p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
3869   photon = p1+p2;
3870   photon.M()*/
3871
3872   Float_t tantheta1 = 1e10;
3873
3874   if (eta1 < -1e-10 || eta1 > 1e-10)
3875     tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
3876
3877   Float_t tantheta2 = 1e10;
3878   if (eta2 < -1e-10 || eta2 > 1e-10)
3879     tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
3880
3881   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3882   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3883
3884   Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
3885
3886   return mass2;
3887 }
3888 //---------------------------------------------------------------------------------
3889
3890 Float_t AliTwoParticlePIDCorr::GetInvMassSquaredCheap(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
3891 {
3892   // calculate inv mass squared approximately
3893
3894   Float_t tantheta1 = 1e10;
3895
3896   if (eta1 < -1e-10 || eta1 > 1e-10)
3897   {
3898     Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
3899     tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3900   }
3901
3902   Float_t tantheta2 = 1e10;
3903   if (eta2 < -1e-10 || eta2 > 1e-10)
3904   {
3905     Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
3906     tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3907   }
3908
3909   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3910   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3911
3912   // fold onto 0...pi
3913   Float_t deltaPhi = TMath::Abs(phi1 - phi2);
3914   while (deltaPhi > TMath::TwoPi())
3915     deltaPhi -= TMath::TwoPi();
3916   if (deltaPhi > TMath::Pi())
3917     deltaPhi = TMath::TwoPi() - deltaPhi;
3918
3919   Float_t cosDeltaPhi = 0;
3920   if (deltaPhi < TMath::Pi()/3)
3921     cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
3922   else if (deltaPhi < 2*TMath::Pi()/3)
3923     cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
3924   else
3925     cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
3926
3927   Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
3928
3929 //   Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
3930
3931   return mass2;
3932 }
3933 //--------------------------------------------------------------------------------
3934 Float_t  AliTwoParticlePIDCorr::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign)
3935 {
3936   //
3937   // calculates dphistar
3938   //
3939
3940   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
3941
3942   static const Double_t kPi = TMath::Pi();
3943
3944   // circularity
3945 //   if (dphistar > 2 * kPi)
3946 //     dphistar -= 2 * kPi;
3947 //   if (dphistar < -2 * kPi)
3948 //     dphistar += 2 * kPi;
3949
3950   if (dphistar > kPi)
3951     dphistar = kPi * 2 - dphistar;
3952   if (dphistar < -kPi)
3953     dphistar = -kPi * 2 - dphistar;
3954   if (dphistar > kPi) // might look funny but is needed
3955     dphistar = kPi * 2 - dphistar;
3956
3957   return dphistar;
3958 }
3959 //_________________________________________________________________________
3960 /*
3961 void AliTwoParticlePIDCorr ::DefineEventPool()
3962 {
3963 Int_t MaxNofEvents=1000;
3964 const Int_t NofVrtxBins=10+(1+10)*2;
3965 Double_t ZvrtxBins[NofVrtxBins+1]={ -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10,
3966                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110,
3967                                       190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
3968
3969 //default values are for centrality
3970 Int_t  NofCentBins=15;
3971 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
3972
3973  if(fCentralityMethod.EndsWith("_MANUAL"))
3974    {
3975  Int_t  NofCentBins=9;
3976  CentralityBins[NofCentBins+1]={0.,9.,14.,19.,26.,34.,44.,58.,80.,500.};//Is This binning is fine for pp, or we don't require them....
3977    }
3978 fPoolMgr = new AliEventPoolManager(MaxNofEvents,fMaxNofMixingTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
3979
3980
3981
3982
3983 fPoolMgr->SetTargetValues(fMaxNofMixingTracks, 0.1, 5);
3984
3985 //if(!fPoolMgr) return kFALSE;
3986 //return kTRUE;
3987
3988 }
3989 */
3990 //------------------------------------------------------------------------
3991 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
3992 {
3993   // This method is a copy from AliUEHist::GetBinning
3994   // takes the binning from <configuration> identified by <tag>
3995   // configuration syntax example:
3996   // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
3997   // phi: .....
3998   //
3999   // returns bin edges which have to be deleted by the caller
4000
4001   TString config(configuration);
4002   TObjArray* lines = config.Tokenize("\n");
4003   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
4004   {
4005     TString line(lines->At(i)->GetName());
4006     if (line.BeginsWith(TString(tag) + ":"))
4007     {
4008       line.Remove(0, strlen(tag) + 1);
4009       line.ReplaceAll(" ", "");
4010       TObjArray* binning = line.Tokenize(",");
4011       Double_t* bins = new Double_t[binning->GetEntriesFast()];
4012       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
4013         bins[j] = TString(binning->At(j)->GetName()).Atof();
4014
4015       nBins = binning->GetEntriesFast() - 1;
4016
4017       delete binning;
4018       delete lines;
4019       return bins;
4020     }
4021   }
4022
4023   delete lines;
4024   AliFatal(Form("Tag %s not found in %s", tag, configuration));
4025   return 0;
4026 }
4027
4028 //____________________________________________________________________
4029
4030 Bool_t AliTwoParticlePIDCorr::CheckTrack(AliAODTrack * part)
4031 {
4032   // check if the track status flags are set
4033
4034   UInt_t status=((AliAODTrack*)part)->GetStatus();
4035   if ((status & fTrackStatus) == fTrackStatus)
4036     return kTRUE;
4037   return kFALSE;
4038 }
4039 //________________________________________________________________________
4040
4041 Bool_t AliTwoParticlePIDCorr::AcceptEventCentralityWeight(Double_t centrality)
4042 {
4043   // rejects "randomly" events such that the centrality gets flat
4044   // uses fCentralityWeights histogram
4045
4046   // TODO code taken and adapted from AliRDHFCuts; waiting for general class AliCentralityFlattening
4047
4048   Double_t weight = fCentralityWeights->GetBinContent(fCentralityWeights->FindBin(centrality));
4049   Double_t centralityDigits = centrality*100. - (Int_t)(centrality*100.);
4050
4051   Bool_t result = kFALSE;
4052   if (centralityDigits < weight)
4053     result = kTRUE;
4054
4055   AliInfo(Form("Centrality: %f; Digits: %f; Weight: %f; Result: %d", centrality, centralityDigits, weight, result));
4056
4057   return result;
4058 }
4059
4060 //____________________________________________________________________
4061 void AliTwoParticlePIDCorr::ShiftTracks(TObjArray* tracks, Double_t angle)
4062 {
4063   // shifts the phi angle of all tracks by angle
4064   // 0 <= angle <= 2pi
4065
4066   for (Int_t i=0; i<tracks->GetEntriesFast(); ++i)
4067   {
4068    LRCParticlePID *part=(LRCParticlePID*)(tracks->UncheckedAt(i));
4069
4070     Double_t newAngle = part->Phi() + angle;
4071     if (newAngle >= TMath::TwoPi())
4072       newAngle -= TMath::TwoPi();
4073
4074     part->SetPhi(newAngle);
4075   }
4076 }
4077
4078
4079 //________________________________________________________________________
4080 void  AliTwoParticlePIDCorr::SetVZEROCalibrationFile(const char* filename,const char* lhcPeriod) {
4081   //Function to setup the VZERO gain equalization
4082     //============Get the equilization map============//
4083   TFile *calibrationFile = TFile::Open(filename);
4084   if((!calibrationFile)||(!calibrationFile->IsOpen())) {
4085     Printf("No calibration file found!!!");
4086     return;
4087   }
4088
4089   TList *list = dynamic_cast<TList *>(calibrationFile->Get(lhcPeriod));
4090   if(!list) {
4091     Printf("Calibration TList not found!!!");
4092     return;
4093   }
4094
4095   fHistVZEROAGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROAGainEqualizationMap"));
4096   if(!fHistVZEROAGainEqualizationMap) {
4097     Printf("VZERO-A calibration object not found!!!");
4098     return;
4099   }
4100   fHistVZEROCGainEqualizationMap = dynamic_cast<TH1F *>(list->FindObject("gHistVZEROCGainEqualizationMap"));
4101   if(!fHistVZEROCGainEqualizationMap) {
4102     Printf("VZERO-C calibration object not found!!!");
4103     return;
4104   }
4105
4106   fHistVZEROChannelGainEqualizationMap = dynamic_cast<TH2F *>(list->FindObject("gHistVZEROChannelGainEqualizationMap"));
4107   if(!fHistVZEROChannelGainEqualizationMap) {
4108     Printf("VZERO channel calibration object not found!!!");
4109     return;
4110   }
4111 }
4112
4113 //________________________________________________________________________
4114 Double_t AliTwoParticlePIDCorr::GetChannelEqualizationFactor(Int_t run,Int_t channel) {
4115   //
4116   if(!fHistVZEROAGainEqualizationMap) return 1.0;
4117
4118   for(Int_t iBinX = 1; iBinX <= fHistVZEROChannelGainEqualizationMap->GetNbinsX(); iBinX++) {
4119     Int_t gRunNumber = atoi(fHistVZEROChannelGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4120     if(gRunNumber == run)
4121       return fHistVZEROChannelGainEqualizationMap->GetBinContent(iBinX,channel+1);
4122   }
4123
4124   return 1.0;
4125 }
4126
4127 //________________________________________________________________________
4128 Double_t AliTwoParticlePIDCorr::GetEqualizationFactor(Int_t run, const char* side) {
4129   //
4130   if(!fHistVZEROAGainEqualizationMap) return 1.0;
4131
4132   TString gVZEROSide = side;
4133   for(Int_t iBinX = 1; iBinX < fHistVZEROAGainEqualizationMap->GetNbinsX(); iBinX++) {
4134     Int_t gRunNumber = atoi(fHistVZEROAGainEqualizationMap->GetXaxis()->GetBinLabel(iBinX));
4135     //cout<<"Looking for run "<<run<<" - current run: "<<gRunNumber<<endl;
4136     if(gRunNumber == run) {
4137       if(gVZEROSide == "A")
4138         return fHistVZEROAGainEqualizationMap->GetBinContent(iBinX);
4139       else if(gVZEROSide == "C")
4140         return fHistVZEROCGainEqualizationMap->GetBinContent(iBinX);
4141     }
4142   }
4143
4144   return 1.0;
4145 }
4146 //________________________________________________________________________
4147 Double_t AliTwoParticlePIDCorr::GetReferenceMultiplicityVZEROFromAOD(AliAODEvent *event){
4148   //Function that returns the reference multiplicity from AODs (data or reco MC, Not for Truth)
4149   //Different ref. mult. implemented: V0M, V0A, V0C, TPC
4150   Double_t gRefMultiplicity = 0., gRefMultiplicityTPC = 0.;
4151   Double_t gRefMultiplicityVZERO = 0., gRefMultiplicityVZEROA = 0., gRefMultiplicityVZEROC = 0.;
4152
4153   AliAODHeader *header = dynamic_cast<AliAODHeader *>(event->GetHeader());
4154   if(!header) {
4155     Printf("ERROR: AOD header not available");
4156     return -999;
4157   }
4158   Int_t gRunNumber = header->GetRunNumber();
4159  Float_t bSign1=header->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
4160
4161
4162  for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++)
4163 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
4164   AliAODTrack* track = dynamic_cast<AliAODTrack*>(event->GetTrack(itrk));
4165   if (!track) continue;
4166   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1,kFALSE);//don't fill the histos here
4167   if(tracktype!=1) continue;
4168
4169   if(!track) continue;//for safety
4170
4171     gRefMultiplicityTPC += 1.0;
4172
4173  }//track looop ends
4174
4175  if(fCentralityMethod == "V0A_MANUAL" || fCentralityMethod == "V0M_MANUAL" || fCentralityMethod == "V0C_MANUAL" ){
4176   //VZERO segmentation in two detectors (0-31: VZERO-C, 32-63: VZERO-A)
4177   for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
4178     fHistVZEROSignal->Fill(iChannel,event->GetVZEROEqMultiplicity(iChannel));
4179
4180     if(iChannel < 32)
4181       gRefMultiplicityVZEROC += event->GetVZEROEqMultiplicity(iChannel);
4182     else if(iChannel >= 32)
4183       gRefMultiplicityVZEROA += event->GetVZEROEqMultiplicity(iChannel);
4184   }//loop over PMTs
4185
4186   //Equalization of gain
4187   Double_t gFactorA = GetEqualizationFactor(gRunNumber,"A");
4188   if(gFactorA != 0)
4189     gRefMultiplicityVZEROA /= gFactorA;
4190   Double_t gFactorC = GetEqualizationFactor(gRunNumber,"C");
4191   if(gFactorC != 0)
4192     gRefMultiplicityVZEROC /= gFactorC;
4193   if((gFactorA != 0)&&(gFactorC != 0))
4194     gRefMultiplicityVZERO = (gRefMultiplicityVZEROA/gFactorA)+(gRefMultiplicityVZEROC/gFactorC);
4195
4196
4197   //EQVZERO vs TPC multiplicity
4198   fHistEQVZEROvsTPCmultiplicity->Fill(gRefMultiplicityVZERO,gRefMultiplicityTPC);
4199   fHistEQVZEROAvsTPCmultiplicity->Fill(gRefMultiplicityVZEROA,gRefMultiplicityTPC);
4200   fHistEQVZEROCvsTPCmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityTPC);
4201
4202   //EQVZERO vs VZERO multiplicity
4203   fHistVZEROCvsEQVZEROCmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),gRefMultiplicityVZEROC);
4204   fHistVZEROAvsEQVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0A(),gRefMultiplicityVZEROA);
4205
4206   //VZEROC vs VZEROA multiplicity
4207   fHistVZEROCvsVZEROAmultiplicity->Fill(event->GetVZEROData()->GetMTotV0C(),event->GetVZEROData()->GetMTotV0A());
4208
4209   //EQVZEROC vs EQVZEROA multiplicity
4210   fHistEQVZEROCvsEQVZEROAmultiplicity->Fill(gRefMultiplicityVZEROC,gRefMultiplicityVZEROA);
4211  }
4212
4213  if(fCentralityMethod == "TRACKS_MANUAL") {
4214     gRefMultiplicity = gRefMultiplicityTPC;
4215     fHistRefmult->Fill(3.,gRefMultiplicityTPC);
4216  }
4217  else if(fCentralityMethod == "V0M_MANUAL"){
4218     gRefMultiplicity = gRefMultiplicityVZERO;
4219     fHistRefmult->Fill(2.,gRefMultiplicityVZERO);
4220
4221  }
4222  else if(fCentralityMethod == "V0A_MANUAL"){
4223     gRefMultiplicity = gRefMultiplicityVZEROA;
4224     fHistRefmult->Fill(0.,gRefMultiplicityVZEROA);
4225
4226  }
4227  else if(fCentralityMethod == "V0C_MANUAL"){
4228     gRefMultiplicity = gRefMultiplicityVZEROC;
4229     fHistRefmult->Fill(1.,gRefMultiplicityVZEROC);
4230  }
4231
4232
4233   return gRefMultiplicity;
4234 }
4235
4236 //-------------------------------------------------------------------------------------------------------
4237 Double_t AliTwoParticlePIDCorr::GetRefMultiOrCentrality(AliAODEvent *event, Bool_t truth){
4238
4239   if(!event) return -1;
4240   // get centrality object and check quality
4241   Double_t cent_v0=-1;
4242   Double_t nooftrackstruth=0;
4243
4244 if(fCentralityMethod=="V0M" || fCentralityMethod=="V0A" || fCentralityMethod=="V0C" || fCentralityMethod=="CL1" || fCentralityMethod=="ZNA" || fCentralityMethod=="V0AEq" || fCentralityMethod=="V0CEq" || fCentralityMethod=="V0MEq")//for PbPb, pPb, pp7TeV(still to be introduced)//data or RecoMC and also for TRUTH
4245     {
4246   AliCentrality *centralityObj=0;
4247   AliAODHeader *header = (AliAODHeader*) event->GetHeader();
4248   if(!header) return -1;
4249   centralityObj = header->GetCentralityP();
4250   // if (centrality->GetQuality() != 0) return ;
4251
4252   if(centralityObj)
4253   {
4254   fHistCentStats->Fill(0.,centralityObj->GetCentralityPercentile("V0A"));
4255   fHistCentStats->Fill(1.,centralityObj->GetCentralityPercentile("V0C"));
4256   fHistCentStats->Fill(2.,centralityObj->GetCentralityPercentile("V0M"));
4257 if(fSampleType=="pp")   fHistCentStats->Fill(3.,centralityObj->GetCentralityPercentile("V0AEq"));//only available for LHC10d at present (Quantile info)
4258 if(fSampleType=="pp")   fHistCentStats->Fill(4.,centralityObj->GetCentralityPercentile("V0CEq"));//only available for LHC10d at present (Quantile info)
4259 if(fSampleType=="pp")   fHistCentStats->Fill(5.,centralityObj->GetCentralityPercentile("V0MEq"));//only available for LHC10d at present (Quantile info)
4260
4261 if(fSampleType=="pPb" || fSampleType=="PbPb")      fHistCentStats->Fill(6.,centralityObj->GetCentralityPercentile("CL1"));
4262 if(fSampleType=="pPb" || fSampleType=="PbPb")      fHistCentStats->Fill(7.,centralityObj->GetCentralityPercentile("ZNA"));
4263
4264       cent_v0 = centralityObj->GetCentralityPercentile(fCentralityMethod);
4265   }
4266   else cent_v0= -1;
4267     }//centralitymethod condition
4268
4269  else if(fCentralityMethod=="V0M_MANUAL" || fCentralityMethod=="V0A_MANUAL" || fCentralityMethod=="V0C_MANUAL" || fCentralityMethod=="TRACKS_MANUAL")//data or RecoMc and also for TRUTH
4270    {
4271      if(!truth){//for data or RecoMC
4272     cent_v0 = GetReferenceMultiplicityVZEROFromAOD(event);
4273    }//for data or RecoMC
4274
4275     if(truth){//condition for TRUTH case
4276 //check for TClonesArray(truth track MC information)
4277 fArrayMC = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
4278   if (!fArrayMC) {
4279     //AliFatal("Error: MC particles branch not found!\n");
4280     return -1;
4281   }
4282 //now process the truth particles(for both efficiency & correlation function)
4283 Int_t nMCTrack = fArrayMC->GetEntriesFast();
4284
4285 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
4286 {//MC truth track loop starts
4287
4288 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
4289
4290     if(!partMC){
4291       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
4292       continue;
4293     }
4294
4295 //consider only charged particles
4296     if(partMC->Charge() == 0) continue;
4297
4298 //consider only primary particles; neglect all secondary particles including from weak decays
4299  if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
4300
4301
4302 //remove injected signals(primaries above <maxLabel>)
4303  if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
4304
4305 //remove duplicates
4306   Bool_t isduplicate=kFALSE;
4307  if (fRemoveDuplicates)
4308    {
4309  for (Int_t j=iMC+1; j<nMCTrack; ++j)
4310    {//2nd trutuh loop starts
4311 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
4312    if(!partMC2){
4313       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
4314       continue;
4315     }
4316  if (partMC->GetLabel() == partMC2->GetLabel())
4317    {
4318 isduplicate=kTRUE;
4319  break;
4320    }
4321    }//2nd truth loop ends
4322    }
4323  if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
4324
4325
4326       if (fCentralityMethod=="V0M_MANUAL") {
4327         if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8)    continue;
4328         if (partMC->Eta() < -3.7 || partMC->Eta() > -1.7) continue;
4329 }
4330       else if (fCentralityMethod=="V0A_MANUAL") {
4331         if(partMC->Eta() > 5.1 || partMC->Eta() < 2.8)  continue;}
4332       else if (fCentralityMethod=="V0C_MANUAL") {
4333         if(partMC->Eta() > -1.7 || partMC->Eta() < -3.7)  continue;}
4334       else if (fCentralityMethod=="TRACKS_MANUAL") {
4335         if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4336         if (partMC->Pt() < fminPt ||  partMC->Pt() > fmaxPt) continue;
4337            }
4338       else{//basically returns the tracks manual case
4339 //give only kinematic cuts at the truth level
4340        if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
4341        if (partMC->Pt() < fminPt ||  partMC->Pt() > fmaxPt) continue;
4342       }
4343
4344  //To determine multiplicity in case of PP
4345  nooftrackstruth+= 1;;
4346
4347  }//truth track loop ends
4348  cent_v0=nooftrackstruth;
4349
4350     }//condition for TRUTH case
4351
4352    }//end of MANUAL method
4353
4354  else if ((fAnalysisType == "MCAOD") && (fCentralityMethod == "MC_b"))//TRUTH MC
4355     {
4356     AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4357     if (!header)
4358     return -1;
4359
4360       AliGenEventHeader* eventHeader = header->GetCocktailHeader(0);  // get first MC header from either ESD/AOD (including cocktail header if available)
4361       if (!eventHeader)
4362       {
4363         // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
4364         // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
4365         AliError("Event header not found. Skipping this event.");
4366         return -1;
4367       }
4368
4369       AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4370
4371
4372      if (collGeometry)   cent_v0 = collGeometry->ImpactParameter();
4373       else cent_v0=-1.;
4374     }//end of Impact parameter method
4375
4376 //else return -1
4377  else cent_v0=-1.;
4378
4379  return cent_v0;
4380 }
4381 //-----------------------------------------------------------------------------------------
4382 Double_t AliTwoParticlePIDCorr::GetAcceptedEventMultiplicity(AliAODEvent *aod,Bool_t truth){
4383   //do the event selection(zvtx, pileup, centrality/multiplicity cut) and then return the value of the centrality of that event
4384   if(!aod) return -1;
4385
4386   Float_t gRefMultiplicity = -1.;
4387
4388   // check first event in chunk (is not needed for new reconstructions)
4389   if(fCheckFirstEventInChunk){
4390     AliAnalysisUtils ut;
4391     if(ut.IsFirstEventInChunk(aod))
4392       return -1.;
4393   }
4394
4395  if(frejectPileUp){
4396     AliAnalysisUtils ut;
4397     ut.SetUseMVPlpSelection(kTRUE);
4398     ut.SetUseOutOfBunchPileUp(kTRUE);
4399     if(ut.IsPileUpEvent(aod))
4400       return -1.;
4401   }
4402
4403 //count events after pileup selection
4404    fEventCounter->Fill(3);
4405
4406   //vertex selection(is it fine for PP?)
4407  if ( fVertextype==1){//for pPb basically if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4408   trkVtx = aod->GetPrimaryVertex();
4409   if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;
4410   TString vtxTtl = trkVtx->GetTitle();
4411   if (!vtxTtl.Contains("VertexerTracks")) return -1;
4412    zvtx = trkVtx->GetZ();
4413   const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
4414   if (!spdVtx || spdVtx->GetNContributors()<=0) return -1;
4415   TString vtxTyp = spdVtx->GetTitle();
4416   Double_t cov[6]={0};
4417   spdVtx->GetCovarianceMatrix(cov);
4418   Double_t zRes = TMath::Sqrt(cov[5]);
4419   if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return -1;
4420    if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return -1;
4421   }
4422   else if(fVertextype==2) {//for pp and pb-pb case , taken from Jan's code
4423         Int_t nVertex = aod->GetNumberOfVertices();
4424         if( nVertex > 0 ) {
4425      trkVtx = aod->GetPrimaryVertex();
4426                 Int_t nTracksPrim = trkVtx->GetNContributors();
4427                  zvtx = trkVtx->GetZ();
4428                 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
4429                 // Reject TPC only vertex
4430                 TString name(trkVtx->GetName());
4431                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return -1;
4432
4433                 // Select a quality vertex by number of tracks?
4434                 if( nTracksPrim < fnTracksVertex ) {
4435                   //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
4436                         return -1;
4437                         }
4438                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
4439                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
4440                 //  return kFALSE;
4441                 //      if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
4442         }
4443         else return -1;
4444
4445   }
4446  else if(fVertextype==0){//default case
4447   trkVtx = aod->GetPrimaryVertex();
4448   if (!trkVtx || trkVtx->GetNContributors()<=0) return -1;//proper number of contributors
4449   zvtx = trkVtx->GetZ();
4450   Double32_t fCov[6];
4451   trkVtx->GetCovarianceMatrix(fCov);
4452   if(fCov[5] == 0) return -1;//proper vertex resolution
4453   }
4454   else {
4455    AliInfo("Wrong Vertextype set for Primary-vertex Selection: event REJECTED ...");
4456    return -1;//as there is no proper sample type
4457   }
4458
4459 fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
4460
4461 //count events having a proper vertex
4462    fEventCounter->Fill(5);
4463
4464  if (TMath::Abs(zvtx) > fzvtxcut) return -1;
4465
4466 //count events after vertex cut
4467   fEventCounter->Fill(7);
4468
4469
4470  //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
4471
4472  fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
4473
4474  //get the centrality or multiplicity
4475  if(truth)  {gRefMultiplicity = GetRefMultiOrCentrality(aod,kTRUE);}//kTRUE-->for Truth case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
4476
4477  else {gRefMultiplicity = GetRefMultiOrCentrality(aod,kFALSE);}//kFALSE-->for data and RecoMc case(only meaningful in case of ref multiplicity,in case of centrality it has no meaning)
4478
4479   if(gRefMultiplicity<0) return -1;
4480
4481  // take events only within the  multiplicity class mentioned in the custom binning
4482   if(gRefMultiplicity < fmincentmult || gRefMultiplicity > fmaxcentmult) return -1;
4483
4484 //count events having proper centrality/ref multiplicity
4485   fEventCounter->Fill(9);
4486
4487
4488 // centrality weighting (optional for 2011 if central and semicentral triggers are used);only for data and recoMC
4489  if (fCentralityWeights && !AcceptEventCentralityWeight(gRefMultiplicity))//**********************
4490   {
4491     AliInfo(Form("Rejecting event because of centrality weighting: %f", gRefMultiplicity));
4492     return -1;
4493   }
4494
4495 //count events after rejection due to centrality weighting
4496   fEventCounter->Fill(11);
4497
4498   return gRefMultiplicity;
4499
4500 }
4501 //--------------------------------------------------------------------------------------------------------
4502 Float_t AliTwoParticlePIDCorr::GetEventPlane(AliAODEvent *event,Bool_t truth, Double_t v0Centr)
4503 {
4504   // Get the event plane
4505  //reset Q vector info
4506
4507     Int_t run = event->GetRunNumber();
4508
4509     if(run != fRun){
4510         // Load the calibrations run dependent
4511       if(! fIsAfter2011) OpenInfoCalbration(run);
4512       fRun=run;
4513     }
4514
4515
4516   Int_t iC = -1;
4517 if (v0Centr < 80){ // analysis only for 0-80% centrality classes
4518  // centrality bins
4519     if(v0Centr < 5) iC = 0;
4520     else if(v0Centr < 10) iC = 1;
4521     else if(v0Centr < 20) iC = 2;
4522     else if(v0Centr < 30) iC = 3;
4523     else if(v0Centr < 40) iC = 4;
4524     else if(v0Centr < 50) iC = 5;
4525     else if(v0Centr < 60) iC = 6;
4526     else if(v0Centr < 70) iC = 7;
4527     else iC = 8;
4528
4529
4530      Int_t iCcal = iC;
4531
4532  //reset Q vector info
4533     Double_t Qxa2 = 0, Qya2 = 0;
4534     Double_t Qxc2 = 0, Qyc2 = 0;
4535     Double_t Qxa3 = 0, Qya3 = 0;
4536     Double_t Qxc3 = 0, Qyc3 = 0;
4537
4538
4539   //MC: from reaction plane
4540  if(truth)
4541 {
4542     AliAODMCHeader* header = (AliAODMCHeader*) event->GetList()->FindObject(AliAODMCHeader::StdBranchName());
4543     if (header){
4544       evplaneMC = header->GetReactionPlaneAngle();//[0, 360]
4545         //make it within [-pi/2,pi/2] to make it general
4546         if(evplaneMC > TMath::Pi()/2 && evplaneMC <=  TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
4547         else if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
4548          fHistEventPlaneTruth->Fill(iC,evplaneMC);
4549         /*
4550         AliGenEventHeader* eventHeader = header->GetCocktailHeader(0);  // get first MC header from either ESD/AOD (including cocktail header if available)
4551       if (eventHeader)
4552       {
4553
4554         AliCollisionGeometry* collGeometry = dynamic_cast<AliCollisionGeometry*> (eventHeader);
4555
4556         if (collGeometry){//get the reaction plane from MC header
4557           gReactionPlane = collGeometry->ReactionPlaneAngle();//[0,180]
4558  }
4559       }
4560         */
4561      //taken from vnv0 code(get the TPC, V0A, V0C event plane using truth tracks)
4562        TClonesArray *mcArray = NULL;
4563         mcArray = (TClonesArray*)event->GetList()->FindObject(AliAODMCParticle::StdBranchName());
4564         if(mcArray){
4565           Float_t QxMCv2[3] = {0,0,0};
4566           Float_t QyMCv2[3] = {0,0,0};
4567           Float_t QxMCv3[3] = {0,0,0};
4568           Float_t QyMCv3[3] = {0,0,0};
4569           Float_t EvPlaneMCV2[3] = {0,0,0};
4570           Float_t EvPlaneMCV3[3] = {0,0,0};
4571           Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
4572           Float_t etaMax[3] = {4.88,-1.8,0.8};
4573
4574           // analysis on MC tracks
4575           Int_t nMCtrack = mcArray->GetEntries() ;
4576
4577           // EP computation with MC tracks
4578           for(Int_t iT=0;iT < nMCtrack;iT++){
4579             AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
4580             if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
4581
4582             Float_t eta = mctr->Eta();
4583   for(Int_t iD=0;iD<3;iD++){
4584               if(eta > etaMin[iD] && eta < etaMax[iD]){
4585                 Float_t phi = mctr->Phi();
4586                 QxMCv2[iD] += TMath::Cos(2*phi);
4587                 QyMCv2[iD] += TMath::Sin(2*phi);
4588                 QxMCv3[iD] += TMath::Cos(3*phi);
4589                 QyMCv3[iD] += TMath::Sin(3*phi);
4590               }
4591             }
4592           }
4593
4594             EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
4595             EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
4596             EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
4597             fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
4598             fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
4599             fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
4600             fgPsi2v0aMC = EvPlaneMCV2[0];
4601             fgPsi2v0cMC = EvPlaneMCV2[1];
4602             fgPsi2tpcMC = EvPlaneMCV2[2];
4603
4604
4605             EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
4606             EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
4607             EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
4608             fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
4609             fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
4610             fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
4611             fgPsi3v0aMC = EvPlaneMCV3[0];
4612             fgPsi3v0cMC = EvPlaneMCV3[1];
4613             fgPsi3tpcMC = EvPlaneMCV3[2];
4614
4615         }
4616     }
4617
4618     }
4619  else{
4620     Int_t nAODTracks = event->GetNumberOfTracks();
4621
4622 // TPC EP needed for resolution studies (TPC subevent)
4623    //AliEventplane * ep = (fAOD->GetHeader())->GetEventplaneP();
4624    //Double_t psiTPC = ep->GetEventplane("Q", fAOD, 2); // in range of [0, pi]
4625     Double_t Qx2 = 0, Qy2 = 0;
4626     Double_t Qx3 = 0, Qy3 = 0;
4627
4628     for(Int_t iT = 0; iT < nAODTracks; iT++) {
4629
4630       AliAODTrack* aodTrack = event->GetTrack(iT);
4631
4632       if (!aodTrack){
4633         continue;
4634       }
4635
4636       Bool_t trkFlag = aodTrack->TestFilterBit(1);
4637
4638       if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster)  || !trkFlag)
4639         continue;
4640
4641       Double_t b[2] = {-99., -99.};
4642       Double_t bCov[3] = {-99., -99., -99.};
4643
4644
4645       AliAODTrack param(*aodTrack);
4646       if (!param.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov)){
4647         continue;
4648       }
4649
4650       if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
4651         continue;
4652
4653       Qx2 += TMath::Cos(2*aodTrack->Phi());
4654       Qy2 += TMath::Sin(2*aodTrack->Phi());
4655       Qx3 += TMath::Cos(3*aodTrack->Phi());
4656       Qy3 += TMath::Sin(3*aodTrack->Phi());
4657
4658     }
4659
4660    Float_t evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
4661    Float_t evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
4662
4663     fgPsi2tpc = evPlAng2;
4664     fgPsi3tpc = evPlAng3;
4665
4666      fPhiRPTPC->Fill(iC,evPlAng2);
4667      fPhiRPTPCv3->Fill(iC,evPlAng3);
4668
4669
4670
4671 //V0 info
4672     AliAODVZERO* aodV0 = event->GetVZEROData();
4673
4674     for (Int_t iv0 = 0; iv0 < 64; iv0++) {
4675       Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
4676       Float_t multv0 = aodV0->GetMultiplicity(iv0);
4677
4678       if(! fIsAfter2011){
4679         if(fAnalysisType == "AOD"){//not for reco MC tracks, only for real data
4680           if (iv0 < 32){ // V0C
4681             Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4682             Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4683             Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4684             Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4685           } else {       // V0A
4686             Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4687             Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4688             Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4689             Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
4690           }
4691         }
4692         else{
4693           if (iv0 < 32){ // V0C
4694             Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4695             Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4696             Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4697             Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
4698           } else {       // V0A
4699             Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4700             Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4701             Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4702             Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
4703           }
4704         }
4705       }
4706     }
4707    //grab for each centrality the proper histo with the Qx and Qy to do the recentering
4708     Double_t Qxamean2 = fMeanQ[iCcal][1][0];
4709     Double_t Qxarms2  = fWidthQ[iCcal][1][0];
4710     Double_t Qyamean2 = fMeanQ[iCcal][1][1];
4711     Double_t Qyarms2  = fWidthQ[iCcal][1][1];
4712     Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
4713     Double_t Qxarms3  = fWidthQv3[iCcal][1][0];
4714     Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
4715     Double_t Qyarms3  = fWidthQv3[iCcal][1][1];
4716
4717     Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
4718     Double_t Qxcrms2  = fWidthQ[iCcal][0][0];
4719     Double_t Qycmean2 = fMeanQ[iCcal][0][1];
4720     Double_t Qycrms2  = fWidthQ[iCcal][0][1];
4721     Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
4722     Double_t Qxcrms3  = fWidthQv3[iCcal][0][0];
4723     Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
4724     Double_t Qycrms3  = fWidthQv3[iCcal][0][1];
4725
4726     Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
4727     Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
4728     Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
4729     Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
4730     Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
4731     Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
4732     Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
4733     Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
4734     /*
4735     //to calculate 2nd order event plane with v0M
4736  Double_t QxCor2 = (Qxa2 - Qxamean2 + Qxc2 - Qxcmean2)
4737     /TMath::Sqrt(Qxarms2*Qxarms2 + Qxcrms2*Qxcrms2);
4738   Double_t QyCor2 = (Qya2 - Qyamean2 + Qyc2 - Qycmean2)
4739     /TMath::Sqrt(Qyarms2*Qyarms2 + Qycrms2*Qycrms2);
4740
4741   //here the calculated event plane is within -Pi to +Pi(delete it , no use here , only for definition)
4742   Double_t psiV0A =(TMath::Pi() + TMath::ATan2(-QyaCor2, -QxaCor2))/2.;
4743   Double_t psiV0C = (TMath::Pi() + TMath::ATan2(-QycCor2, -QxcCor2))/2.;
4744   Double_t psiVZero = (TMath::Pi() + TMath::ATan2(-QyCor2, -QxCor2))/2.;
4745
4746     */
4747
4748     Float_t evPlAngV0ACor2=999.;
4749     Float_t evPlAngV0CCor2=999.;
4750     Float_t evPlAngV0ACor3=999.;
4751     Float_t evPlAngV0CCor3=999.;
4752
4753    if(! fIsAfter2011){
4754       if(fAnalysisType == "AOD"){
4755         evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
4756         evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
4757         evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
4758         evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
4759       }
4760       else{
4761         evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
4762         evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
4763         evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
4764         evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
4765       }
4766     }
4767     else{
4768       AliEventplane *ep =  event->GetEventplane();
4769       evPlAngV0ACor2 = ep->GetEventplane("V0A", event, 2);
4770       evPlAngV0CCor2 = ep->GetEventplane("V0C", event, 2);
4771       evPlAngV0ACor3 = ep->GetEventplane("V0A", event, 3);
4772       evPlAngV0CCor3 = ep->GetEventplane("V0C", event, 3);
4773     }
4774
4775     fgPsi2v0a = evPlAngV0ACor2;
4776     fgPsi2v0c = evPlAngV0CCor2;
4777     fgPsi3v0a = evPlAngV0ACor3;
4778     fgPsi3v0c = evPlAngV0CCor3;
4779
4780  // Fill EP distribution histograms evPlAng
4781
4782      fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
4783      fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
4784
4785      fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
4786      fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
4787
4788     // Fill histograms needed for resolution evaluation
4789     fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
4790     fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
4791     fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
4792
4793     fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
4794     fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
4795     fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
4796
4797
4798     /*
4799  Float_t gVZEROEventPlane    = -10.;
4800   Float_t gReactionPlane      = -10.;
4801   Double_t qxTot = 0.0, qyTot = 0.0;
4802
4803     AliEventplane *ep = event->GetEventplane();
4804     if(ep){
4805       gVZEROEventPlane = ep->CalculateVZEROEventPlane(event,10,2,qxTot,qyTot);
4806       if(gVZEROEventPlane < 0.) gVZEROEventPlane += TMath::Pi();
4807       //gReactionPlane = gVZEROEventPlane*TMath::RadToDeg();
4808       gReactionPlane = gVZEROEventPlane;
4809     }
4810     */
4811   }//AOD,ESD,ESDMC
4812  //return gReactionPlane;
4813
4814  //make the final 2nd order event plane within 0 to Pi
4815      //using data and reco tracks only
4816       if(fgPsi2v0a!=999. && fgPsi2v0a < 0.) fgPsi2v0a += TMath::Pi();
4817       if(fgPsi2v0c!=999. && fgPsi2v0c < 0.) fgPsi2v0c += TMath::Pi();
4818       if(fgPsi2tpc!=999. && fgPsi2tpc < 0.) fgPsi2tpc += TMath::Pi();
4819       //using truth tracks only
4820       if(evplaneMC!=999. && evplaneMC < 0.) evplaneMC += TMath::Pi();
4821       if(fgPsi2v0aMC!=999. && fgPsi2v0aMC < 0.) fgPsi2v0aMC += TMath::Pi();
4822       if(fgPsi2v0cMC!=999. && fgPsi2v0cMC < 0.) fgPsi2v0cMC += TMath::Pi();
4823       if(fgPsi2tpcMC!=999. && fgPsi2tpcMC < 0.) fgPsi2tpcMC += TMath::Pi();
4824       //for the time being leave the 3rd order event planes within -pi/3 t0 +pi/3
4825
4826       if(truth){//for truth MC
4827         if(fV2 && fEPdet=="header")gReactionPlane=evplaneMC;
4828         if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0aMC;
4829         if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0cMC;
4830         if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpcMC;
4831
4832         if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0aMC;
4833         if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0cMC;
4834         if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpcMC;
4835 }
4836       else{//for data and recoMC
4837         if(fV2 && fEPdet=="V0A")gReactionPlane=fgPsi2v0a;
4838         if(fV2 && fEPdet=="V0C")gReactionPlane=fgPsi2v0c;
4839         if(fV2 && fEPdet=="TPC")gReactionPlane=fgPsi2tpc;
4840
4841         if(fV3 && fEPdet=="V0A")gReactionPlane=fgPsi3v0a;
4842         if(fV3 && fEPdet=="V0C")gReactionPlane=fgPsi3v0c;
4843         if(fV3 && fEPdet=="TPC")gReactionPlane=fgPsi3tpc;
4844
4845       }
4846
4847  } //centrality cut condition
4848
4849 return gReactionPlane;
4850 }
4851 //------------------------------------------------------------------------------------------------------------------
4852 void AliTwoParticlePIDCorr::OpenInfoCalbration(Int_t run){
4853     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
4854     TFile *foadb = TFile::Open(oadbfilename.Data());
4855
4856     if(!foadb){
4857         printf("OADB file %s cannot be opened\n",oadbfilename.Data());
4858         return;
4859     }
4860
4861     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
4862     if(!cont){
4863         printf("OADB object hMultV0BefCorr is not available in the file\n");
4864         return;
4865     }
4866
4867     if(!(cont->GetObject(run))){
4868         printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
4869         run = 137366;
4870     }
4871     fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
4872
4873     TF1 *fpol0 = new TF1("fpol0","pol0");
4874     fMultV0->Fit(fpol0,"","",0,31);
4875     fV0Cpol = fpol0->GetParameter(0);
4876     fMultV0->Fit(fpol0,"","",32,64);
4877     fV0Apol = fpol0->GetParameter(0);
4878
4879     for(Int_t iside=0;iside<2;iside++){
4880         for(Int_t icoord=0;icoord<2;icoord++){
4881             for(Int_t i=0;i  < 9;i++){
4882                 char namecont[100];
4883                 if(iside==0 && icoord==0)
4884                   snprintf(namecont,100,"hQxc2_%i",i);
4885                 else if(iside==1 && icoord==0)
4886                   snprintf(namecont,100,"hQxa2_%i",i);
4887                 else if(iside==0 && icoord==1)
4888                   snprintf(namecont,100,"hQyc2_%i",i);
4889                 else if(iside==1 && icoord==1)
4890                   snprintf(namecont,100,"hQya2_%i",i);
4891
4892                 cont = (AliOADBContainer*) foadb->Get(namecont);
4893                 if(!cont){
4894                     printf("OADB object %s is not available in the file\n",namecont);
4895                     return;
4896                 }
4897
4898                 if(!(cont->GetObject(run))){
4899                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
4900                     run = 137366;
4901                 }
4902                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
4903                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
4904
4905                 //for v3
4906                 if(iside==0 && icoord==0)
4907                   snprintf(namecont,100,"hQxc3_%i",i);
4908                 else if(iside==1 && icoord==0)
4909                   snprintf(namecont,100,"hQxa3_%i",i);
4910                 else if(iside==0 && icoord==1)
4911                   snprintf(namecont,100,"hQyc3_%i",i);
4912                 else if(iside==1 && icoord==1)
4913                   snprintf(namecont,100,"hQya3_%i",i);
4914
4915                 cont = (AliOADBContainer*) foadb->Get(namecont);
4916                 if(!cont){
4917                     printf("OADB object %s is not available in the file\n",namecont);
4918                     return;
4919                 }
4920
4921                 if(!(cont->GetObject(run))){
4922                     printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
4923                     run = 137366;
4924                 }
4925                 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
4926                 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
4927
4928             }
4929         }
4930     }
4931 }
4932 //____________________________________________________________________
4933 void AliTwoParticlePIDCorr::FillPIDEventPlane(Double_t centrality,Int_t par,Float_t trigphi,Float_t fReactionPlane)
4934 {
4935
4936  // Event plane (determine psi bin)
4937     Double_t gPsiMinusPhi    =   0.;
4938     Double_t gPsiMinusPhiBin = -10.;
4939 if(fRequestEventPlane){
4940     gPsiMinusPhi   = TMath::Abs(trigphi - fReactionPlane);
4941     //in-plane
4942     if((gPsiMinusPhi <= 7.5*TMath::DegToRad())||
4943       (gPsiMinusPhi >= 352.5*TMath::DegToRad())||
4944        ((172.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 187.5*TMath::DegToRad())))
4945       gPsiMinusPhiBin = 0.0;
4946     //intermediate
4947     else if(((37.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 52.5*TMath::DegToRad()))||
4948             ((127.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 142.5*TMath::DegToRad()))||
4949             ((217.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 232.5*TMath::DegToRad()))||
4950             ((307.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 322.5*TMath::DegToRad())))
4951       gPsiMinusPhiBin = 1.0;
4952     //out of plane
4953     else if(((82.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 97.5*TMath::DegToRad()))||
4954             ((262.5*TMath::DegToRad() <= gPsiMinusPhi)&&(gPsiMinusPhi <= 277.5*TMath::DegToRad())))
4955       gPsiMinusPhiBin = 2.0;
4956     //everything else
4957     else
4958       gPsiMinusPhiBin = 3.0;
4959
4960     fEventPlanePID->Fill(centrality,gPsiMinusPhiBin,(Float_t)par);
4961  }
4962 }
4963
4964 //----------------------------------------------------------
4965 void AliTwoParticlePIDCorr::Terminate(Option_t *)
4966 {
4967   // Draw result to screen, or perform fitting, normalizations
4968   // Called once at the end of the query
4969   fOutput = dynamic_cast<TList*> (GetOutputData(1));
4970   if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
4971
4972
4973 }
4974 //------------------------------------------------------------------