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