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