PWGCF/Correlations/DPhi/TriggerPID/AliTwoParticlePIDCorr.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include "AliTwoParticlePIDCorr.h"
  17 #include "AliVParticle.h"
  18 #include "TFormula.h"
  19 #include "TAxis.h"
  20 #include "TChain.h"
  21 #include "TTree.h"
  22 #include "TH1F.h"
  23 #include "TH2F.h"
  24 #include "TH3F.h"
  25 #include "TList.h"
  26 #include "TFile.h"
  27 #include "TGrid.h"
  28
  29 #include "AliCentrality.h"
  30 #include "Riostream.h"
  31
  32 #include <TSpline.h>
  33 #include <AliPID.h>
  34 #include "AliESDpid.h"
  35 #include "AliAODpidUtil.h"
  36 #include <AliPIDResponse.h>
  37
  38 #include <AliAnalysisManager.h>
  39 #include <AliInputEventHandler.h>
  40 #include "AliAODInputHandler.h"
  41
  42 #include "AliAnalysisTaskSE.h"
  43 #include "AliAnalysisManager.h"
  44 #include "AliCentrality.h"
  45
  46 #include "AliVEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODTrack.h"
  49 #include "AliVTrack.h"
  50
  51 #include "THnSparse.h"
  52
  53 #include "AliAODMCHeader.h"
  54 #include "AliAODMCParticle.h"
  55 #include "AliMCEventHandler.h"
  56 #include "AliMCEvent.h"
  57 #include "AliMCParticle.h"
  58 #include "TParticle.h"
  59 #include <TDatabasePDG.h>
  60 #include <TParticlePDG.h>
  61
  62 #include "AliGenCocktailEventHeader.h"
  63 #include "AliGenEventHeader.h"
  64
  65 #include "AliEventPoolManager.h"
  66 //#include "AliAnalysisUtils.h"
  67 using namespace AliPIDNameSpace;
  68 using namespace std;
  69
  70 ClassImp(AliTwoParticlePIDCorr)
  71 ClassImp(LRCParticlePID)
  72 //________________________________________________________________________
  73 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr() // All data members should be initialised here
  74 :AliAnalysisTaskSE(),
  75   fOutput(0),
  76   fCentralityMethod("V0A"),
  77   fSampleType("pPb"),
  78 fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default)
  79   trkVtx(0),
  80   zvtx(0),
  81   fFilterBit(768),
  82   fzvtxcut(10.0),
  83   ffilltrigassoUNID(kFALSE),
  84   ffilltrigUNIDassoID(kFALSE),
  85   ffilltrigIDassoUNID(kTRUE),
  86   ffilltrigIDassoID(kFALSE),
  87   ffilltrigIDassoIDMCTRUTH(kFALSE),
  88   fPtOrderMCTruth(kFALSE),
  89   fTriggerSpeciesSelection(kFALSE),
  90   fAssociatedSpeciesSelection(kFALSE),
  91   fTriggerSpecies(SpPion),
  92   fAssociatedSpecies(SpPion),
  93   fCustomBinning(""),
  94   fBinningString(""),
  95   fcontainPIDtrig(kTRUE),
  96   fcontainPIDasso(kFALSE),
  97 //frejectPileUp(kFALSE),
  98   fminPt(0.2),
  99   fmaxPt(10.0),
 100   fmineta(-0.8),
 101   fmaxeta(0.8),
 102   fminprotonsigmacut(-6.0),
 103   fmaxprotonsigmacut(-3.0),
 104   fminpionsigmacut(0.0),
 105   fmaxpionsigmacut(4.0),
 106   fselectprimaryTruth(kTRUE),
 107   fonlyprimarydatareco(kFALSE),
 108   fdcacut(kFALSE),
 109   fdcacutvalue(3.0),
 110   ffillhistQAReco(kFALSE),
 111   ffillhistQATruth(kFALSE),
 112   kTrackVariablesPair(0),
 113   fminPtTrig(0),
 114   fmaxPtTrig(0),
 115   fminPtComboeff(2.0),
 116   fmaxPtComboeff(4.0),
 117   fminPtAsso(0),
 118   fmaxPtAsso(0),
 119   fhistcentrality(0),
 120   fEventCounter(0),
 121   fEtaSpectrasso(0),
 122   fphiSpectraasso(0),
 123   MCtruthpt(0),
 124   MCtrutheta(0),
 125   MCtruthphi(0),
 126   MCtruthpionpt(0),
 127   MCtruthpioneta(0),
 128   MCtruthpionphi(0),
 129   MCtruthkaonpt(0),
 130   MCtruthkaoneta(0),
 131   MCtruthkaonphi(0),
 132   MCtruthprotonpt(0),
 133   MCtruthprotoneta(0),
 134   MCtruthprotonphi(0),
 135   fPioncont(0),
 136   fKaoncont(0),
 137   fProtoncont(0),
 138   fCentralityCorrelation(0x0),
 139   fHistoTPCdEdx(0x0),
 140   fHistoTOFbeta(0x0),
 141   fTPCTOFPion3d(0),
 142   fTPCTOFKaon3d(0),
 143   fTPCTOFProton3d(0),
 144   fPionPt(0),
 145   fPionEta(0),
 146   fPionPhi(0),
 147   fKaonPt(0),
 148   fKaonEta(0),
 149   fKaonPhi(0),
 150   fProtonPt(0),
 151   fProtonEta(0),
 152   fProtonPhi(0),
 153   fCorrelatonTruthPrimary(0),
 154   fCorrelatonTruthPrimarymix(0),
 155   fTHnCorrUNID(0),
 156   fTHnCorrUNIDmix(0),
 157   fTHnCorrID(0),
 158   fTHnCorrIDmix(0),
 159   fTHnCorrIDUNID(0),
 160   fTHnCorrIDUNIDmix(0),
 161   fTHnTrigcount(0),
 162   fTHnTrigcountMCTruthPrim(0),
 163   fPoolMgr(0x0),
 164   fArrayMC(0),
 165   fAnalysisType("AOD"),
 166   fefffilename(""),
 167   twoTrackEfficiencyCutValue(0.02),
 168 //fControlConvResoncances(0),
 169   fPID(NULL),
 170  eventno(0),
 171   fPtTOFPIDmin(0.6),
 172   fPtTOFPIDmax(4.0),
 173   fRequestTOFPID(kTRUE),
 174   fPIDType(NSigmaTPCTOF),
 175   fNSigmaPID(3),
 176   fUseExclusiveNSigma(kFALSE),
 177   fRemoveTracksT0Fill(kFALSE),
 178 fSelectCharge(0),
 179 fTriggerSelectCharge(0),
 180 fAssociatedSelectCharge(0),
 181 fTriggerRestrictEta(-1),
 182 fEtaOrdering(kFALSE),
 183 fCutConversions(kFALSE),
 184 fCutResonances(kFALSE),
 185 fRejectResonanceDaughters(-1),
 186   fOnlyOneEtaSide(0),
 187 fInjectedSignals(kFALSE),
 188   fRemoveWeakDecays(kFALSE),
 189 fRemoveDuplicates(kFALSE),
 190   fapplyTrigefficiency(kFALSE),
 191   fapplyAssoefficiency(kFALSE),
 192   ffillefficiency(kFALSE),
 193   fmesoneffrequired(kFALSE),
 194   fkaonprotoneffrequired(kFALSE),
 195 //fAnalysisUtils(0x0),
 196   fDCAXYCut(0)
 197
 198 {
 199  for ( Int_t i = 0; i < 16; i++) {
 200     fHistQA[i] = NULL;
 201   }
 202
 203  for ( Int_t i = 0; i < 6; i++ ){
 204     fTHnrecomatchedallPid[i] = NULL;
 205     fTHngenprimPidTruth[i] = NULL;
 206     effcorection[i]=NULL;
 207     //effmap[i]=NULL;
 208
 209   }
 210
 211  for(Int_t ipart=0;ipart<NSpecies;ipart++)
 212     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
 213       fnsigmas[ipart][ipid]=999.;
 214
 215  for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
 216
 217   }
 218 //________________________________________________________________________
 219 AliTwoParticlePIDCorr::AliTwoParticlePIDCorr(const char *name) // All data members should be initialised here
 220   :AliAnalysisTaskSE(name),
 221  fOutput(0),
 222  fCentralityMethod("V0A"),
 223   fSampleType("pPb"),
 224   fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default)
 225   trkVtx(0),
 226   zvtx(0),
 227   fFilterBit(768),
 228   fzvtxcut(10.0),
 229   ffilltrigassoUNID(kFALSE),
 230   ffilltrigUNIDassoID(kFALSE),
 231   ffilltrigIDassoUNID(kTRUE),
 232   ffilltrigIDassoID(kFALSE),
 233   ffilltrigIDassoIDMCTRUTH(kFALSE),
 234   fPtOrderMCTruth(kFALSE),
 235   fTriggerSpeciesSelection(kFALSE),
 236   fAssociatedSpeciesSelection(kFALSE),
 237   fTriggerSpecies(SpPion),
 238   fAssociatedSpecies(SpPion),
 239   fCustomBinning(""),
 240   fBinningString(""),
 241   fcontainPIDtrig(kTRUE),
 242   fcontainPIDasso(kFALSE),
 243    // frejectPileUp(kFALSE),
 244   fminPt(0.2),
 245   fmaxPt(10.0),
 246   fmineta(-0.8),
 247   fmaxeta(0.8),
 248    fminprotonsigmacut(-6.0),
 249    fmaxprotonsigmacut(-3.0),
 250    fminpionsigmacut(0.0),
 251    fmaxpionsigmacut(4.0),
 252   fselectprimaryTruth(kTRUE),
 253   fonlyprimarydatareco(kFALSE),
 254    fdcacut(kFALSE),
 255   fdcacutvalue(3.0),
 256   ffillhistQAReco(kFALSE),
 257   ffillhistQATruth(kFALSE),
 258  kTrackVariablesPair(0),
 259   fminPtTrig(0),
 260   fmaxPtTrig(0),
 261   fminPtComboeff(2.0),
 262   fmaxPtComboeff(4.0),
 263   fminPtAsso(0),
 264   fmaxPtAsso(0),
 265   fhistcentrality(0),
 266   fEventCounter(0),
 267   fEtaSpectrasso(0),
 268   fphiSpectraasso(0),
 269   MCtruthpt(0),
 270   MCtrutheta(0),
 271   MCtruthphi(0),
 272   MCtruthpionpt(0),
 273   MCtruthpioneta(0),
 274   MCtruthpionphi(0),
 275   MCtruthkaonpt(0),
 276   MCtruthkaoneta(0),
 277   MCtruthkaonphi(0),
 278   MCtruthprotonpt(0),
 279   MCtruthprotoneta(0),
 280   MCtruthprotonphi(0),
 281   fPioncont(0),
 282   fKaoncont(0),
 283   fProtoncont(0),
 284   fCentralityCorrelation(0x0),
 285   fHistoTPCdEdx(0x0),
 286   fHistoTOFbeta(0x0),
 287   fTPCTOFPion3d(0),
 288   fTPCTOFKaon3d(0),
 289   fTPCTOFProton3d(0),
 290   fPionPt(0),
 291   fPionEta(0),
 292   fPionPhi(0),
 293   fKaonPt(0),
 294   fKaonEta(0),
 295   fKaonPhi(0),
 296   fProtonPt(0),
 297   fProtonEta(0),
 298   fProtonPhi(0),
 299   fCorrelatonTruthPrimary(0),
 300  fCorrelatonTruthPrimarymix(0),
 301   fTHnCorrUNID(0),
 302   fTHnCorrUNIDmix(0),
 303   fTHnCorrID(0),
 304   fTHnCorrIDmix(0),
 305   fTHnCorrIDUNID(0),
 306   fTHnCorrIDUNIDmix(0),
 307   fTHnTrigcount(0),
 308   fTHnTrigcountMCTruthPrim(0),
 309    fPoolMgr(0x0),
 310   fArrayMC(0),
 311   fAnalysisType("AOD"),
 312    fefffilename(""),
 313   twoTrackEfficiencyCutValue(0.02),
 314 //fControlConvResoncances(0),
 315   fPID(NULL),
 316  eventno(0),
 317  fPtTOFPIDmin(0.6),
 318   fPtTOFPIDmax(4.0),
 319   fRequestTOFPID(kTRUE),
 320   fPIDType(NSigmaTPCTOF),
 321   fNSigmaPID(3),
 322   fUseExclusiveNSigma(kFALSE),
 323   fRemoveTracksT0Fill(kFALSE),
 324 fSelectCharge(0),
 325 fTriggerSelectCharge(0),
 326 fAssociatedSelectCharge(0),
 327 fTriggerRestrictEta(-1),
 328 fEtaOrdering(kFALSE),
 329 fCutConversions(kFALSE),
 330 fCutResonances(kFALSE),
 331 fRejectResonanceDaughters(-1),
 332   fOnlyOneEtaSide(0),
 333 fInjectedSignals(kFALSE),
 334   fRemoveWeakDecays(kFALSE),
 335 fRemoveDuplicates(kFALSE),
 336   fapplyTrigefficiency(kFALSE),
 337   fapplyAssoefficiency(kFALSE),
 338   ffillefficiency(kFALSE),
 339  fmesoneffrequired(kFALSE),
 340  fkaonprotoneffrequired(kFALSE),
 341    //fAnalysisUtils(0x0),
 342   fDCAXYCut(0)
 343 {
 344
 345    for ( Int_t i = 0; i < 16; i++) {
 346     fHistQA[i] = NULL;
 347   }
 348
 349 for ( Int_t i = 0; i < 6; i++ ){
 350     fTHnrecomatchedallPid[i] = NULL;
 351     fTHngenprimPidTruth[i] = NULL;
 352     effcorection[i]=NULL;
 353     //effmap[i]=NULL;
 354
 355   }
 356
 357  for(Int_t ipart=0;ipart<NSpecies;ipart++)
 358     for(Int_t ipid=0;ipid<=NSigmaPIDType;ipid++)
 359       fnsigmas[ipart][ipid]=999.;
 360
 361    for(Int_t ipart=0;ipart<NSpecies;ipart++) {fHasDoubleCounting[ipart]=kFALSE;}
 362
 363   // The last in the above list should not have a comma after it
 364
 365   // Constructor
 366   // Define input and output slots here (never in the dummy constructor)
 367   // Input slot #0 works with a TChain - it is connected to the default input container
 368   // Output slot #1 writes into a TH1 container
 369
 370   DefineOutput(1, TList::Class());                                        // for output list
 371
 372 }
 373
 374 //________________________________________________________________________
 375 AliTwoParticlePIDCorr::~AliTwoParticlePIDCorr()
 376 {
 377   // Destructor. Clean-up the output list, but not the histograms that are put inside
 378   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 379   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 380     delete fOutput;
 381
 382   }
 383   if (fPID) delete fPID;
 384
 385   }
 386 //________________________________________________________________________
 387 Float_t AliTwoParticlePIDCorr::PhiRange(Float_t DPhi)
 388
 389 {
 390         //
 391         // Puts the argument in the range [-pi/2,3 pi/2].
 392         //
 393
 394         if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
 395         if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
 396
 397         return DPhi;
 398
 399 }
 400 //________________________________________________________________________
 401 void AliTwoParticlePIDCorr::UserCreateOutputObjects()
 402 {
 403   // Create histograms
 404   // Called once (on the worker node)
 405   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 406   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 407   fPID = inputHandler->GetPIDResponse();
 408
 409   //AliAnalysisUtils *fUtils = new AliAnalysisUtils();
 410
 411 //get the efficiency correction map
 412
 413
 414   fOutput = new TList();
 415   fOutput->SetOwner();  // IMPORTANT!
 416
 417     Int_t centmultbins=10;
 418     Double_t centmultmin=0.0;
 419     Double_t centmultmax=100.0;
 420   if(fSampleType=="pPb" || fSampleType=="PbPb")
 421    {
 422       centmultbins=10;
 423       centmultmin=0.0;
 424       centmultmax=100.0;
 425    }
 426  if(fSampleType=="pp")
 427    {
 428      centmultbins=10;
 429      centmultmin=0.0;
 430      centmultmax=200.0;
 431    }
 432
 433 fhistcentrality=new TH1F("fhistcentrality","fhistcentrality",centmultbins*4,centmultmin,centmultmax);
 434 fOutput->Add(fhistcentrality);
 435
 436   fEventCounter = new TH1F("fEventCounter","EventCounter", 10, 0.5,10.5);
 437   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
 438   fEventCounter->GetXaxis()->SetBinLabel(2,"Have a vertex");
 439   fEventCounter->GetXaxis()->SetBinLabel(5,"After vertex Cut");
 440   fEventCounter->GetXaxis()->SetBinLabel(6,"Within 0-100% centrality");
 441   fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
 442   //fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
 443   fOutput->Add(fEventCounter);
 444
 445 fEtaSpectrasso=new TH2F("fEtaSpectraasso","fEtaSpectraasso",180,-0.9,0.9,100,0.,20. );
 446 fOutput->Add(fEtaSpectrasso);
 447
 448 fphiSpectraasso=new TH2F("fphiSpectraasso","fphiSpectraasso",72,0,2*TMath::Pi(),100,0.,20.);
 449 fOutput->Add(fphiSpectraasso);
 450
 451
 452  if(fSampleType=="pPb" || fSampleType=="PbPb"){ fCentralityCorrelation = new TH2D("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 20000, 0,40000);
 453       fOutput->Add(fCentralityCorrelation);
 454  }
 455
 456 fHistoTPCdEdx = new TH2F("hHistoTPCdEdx", ";p_{T} (GeV/c);dE/dx (au.)",200,0.0,10.0,500, 0., 500.);
 457 fOutput->Add(fHistoTPCdEdx);
 458 fHistoTOFbeta = new TH2F(Form("hHistoTOFbeta"), ";p_{T} (GeV/c);v/c",70, 0., 7., 500, 0.1, 1.1);
 459   fOutput->Add(fHistoTOFbeta);
 460
 461    fTPCTOFPion3d=new TH3F ("fTPCTOFpion3d", "fTPCTOFpion3d",100,0., 10., 120,-60.,60.,120,-60.,60);
 462    fOutput->Add(fTPCTOFPion3d);
 463
 464    fTPCTOFKaon3d=new TH3F ("fTPCTOFKaon3d", "fTPCTOFKaon3d",100,0., 10., 120,-60.,60.,120,-60.,60);
 465    fOutput->Add(fTPCTOFKaon3d);
 466
 467    fTPCTOFProton3d=new TH3F ("fTPCTOFProton3d", "fTPCTOFProton3d",100,0., 10., 120,-60.,60.,120,-60.,60);
 468    fOutput->Add(fTPCTOFProton3d);
 469
 470 if(ffillhistQAReco)
 471     {
 472     fPionPt = new TH1F("fHistQAPionPt","p_{T} distribution",200,0.,10.);
 473  fOutput->Add(fPionPt);
 474     fPionEta= new TH1F("fHistQAPionEta","#eta distribution",360,-1.8,1.8);
 475  fOutput->Add(fPionEta);
 476     fPionPhi = new TH1F("fHistQAPionPhi","#phi distribution",340,0,6.8);
 477  fOutput->Add(fPionPhi);
 478
 479     fKaonPt = new TH1F("fHistQAKaonPt","p_{T} distribution",200,0.,10.);
 480  fOutput->Add(fKaonPt);
 481     fKaonEta= new TH1F("fHistQAKaonEta","#eta distribution",360,-1.8,1.8);
 482  fOutput->Add(fKaonEta);
 483     fKaonPhi = new TH1F("fHistQAKaonPhi","#phi distribution",340,0,6.8);
 484  fOutput->Add(fKaonPhi);
 485
 486     fProtonPt = new TH1F("fHistQAProtonPt","p_{T} distribution",200,0.,10.);
 487  fOutput->Add(fProtonPt);
 488     fProtonEta= new TH1F("fHistQAProtonEta","#eta distribution",360,-1.8,1.8);
 489  fOutput->Add(fProtonEta);
 490     fProtonPhi= new TH1F("fHistQAProtonPhi","#phi distribution",340,0,6.8);
 491  fOutput->Add(fProtonPhi);
 492     }
 493
 494   fHistQA[0] = new TH1F("fHistQAvx", "Histo Vx All ", 50, -5., 5.);
 495   fHistQA[1] = new TH1F("fHistQAvy", "Histo Vy All", 50, -5., 5.);
 496   fHistQA[2] = new TH1F("fHistQAvz", "Histo Vz All", 50, -25., 25.);
 497   fHistQA[3] = new TH1F("fHistQAvxA", "Histo Vx  After Cut ", 50, -5., 5.);
 498   fHistQA[4] = new TH1F("fHistQAvyA", "Histo Vy After Cut", 50, -5., 5.);
 499   fHistQA[5] = new TH1F("fHistQAvzA", "Histo Vz After Cut", 50, -25., 25.);
 500   fHistQA[6] = new TH1F("fHistQADcaXyC", "Histo DCAxy after cut", 50, -5., 5.);
 501   fHistQA[7] = new TH1F("fHistQADcaZC", "Histo DCAz after cut", 50, -5., 5.);
 502   fHistQA[8] = new TH1F("fHistQAPt","p_{T} distribution",200,0.,10.);
 503   fHistQA[9] = new TH1F("fHistQAEta","#eta distribution",360,-1.8,1.8);
 504   fHistQA[10] = new TH1F("fHistQAPhi","#phi distribution",340,0,6.8);
 505   fHistQA[11] = new TH1F("fHistQANCls","Number of TPC cluster",200,0,200);
 506   fHistQA[13] = new TH1F("fHistQAChi2","Chi2 per NDF",100,0,10);
 507  fHistQA[12] = new TH1F("fHistQANCls1","Number of TPC cluster1",200,0,200);
 508  fHistQA[14] = new TH1F("nCrossedRowsTPC","Number of TPC ccrossed rows",200,0,200);
 509  fHistQA[15] = new TH1F("ratioCrossedRowsOverFindableClustersTPC","Number of TPC ccrossed rows find clusters",200,0,2);
 510
 511 for(Int_t i = 0; i < 16; i++)
 512     {
 513       fOutput->Add(fHistQA[i]);
 514     }
 515
 516    kTrackVariablesPair=6 ;
 517
 518  if(fcontainPIDtrig && !fcontainPIDasso) kTrackVariablesPair=7;
 519
 520  if(!fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=7;
 521
 522  if(fcontainPIDtrig && fcontainPIDasso) kTrackVariablesPair=8;
 523
 524
 525 // two particle histograms
 526   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
 527   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables
 528   TString* axisTitlePair;  // axis titles for track variables
 529   axisTitlePair=new TString[kTrackVariablesPair];
 530
 531  TString defaultBinningStr;
 532   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"
 533     "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"
 534     "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0,10.0\n"
 535     "p_t_eff: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"
 536     "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
 537   "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
 538         "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"
 539       "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
 540
 541   if(fcontainPIDtrig){
 542       defaultBinningStr += "PIDTrig: -0.5,0.5,1.5,2.5,3.5\n"; // course
 543   }
 544   if(fcontainPIDasso){
 545       defaultBinningStr += "PIDAsso: -0.5,0.5,1.5,2.5,3.5\n"; // course
 546   }
 547  // =========================================================
 548   // Customization (adopted from AliUEHistograms)
 549   // =========================================================
 550
 551   TObjArray* lines = defaultBinningStr.Tokenize("\n");
 552   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
 553   {
 554     TString line(lines->At(i)->GetName());
 555     TString tag = line(0, line.Index(":")+1);
 556     if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
 557       fBinningString += line + "\n";
 558     else
 559       AliInfo(Form("Using custom binning for %s", tag.Data()));
 560   }
 561   delete lines;
 562   fBinningString += fCustomBinning;
 563
 564   AliInfo(Form("Used THn Binning:\n%s",fBinningString.Data()));
 565
 566  // =========================================================
 567   // Now set the bins
 568   // =========================================================
 569
 570     dBinsPair[0]       = GetBinning(fBinningString, "multiplicity", iBinPair[0]);
 571     axisTitlePair[0]   = " multiplicity";
 572
 573     dBinsPair[1]     = GetBinning(fBinningString, "vertex", iBinPair[1]);
 574     axisTitlePair[1]  = "v_{Z} (cm)";
 575
 576     dBinsPair[2]     = GetBinning(fBinningString, "p_t_leading_course", iBinPair[2]);
 577     axisTitlePair[2]    = "p_{T,trig.} (GeV/c)";
 578
 579     dBinsPair[3]     = GetBinning(fBinningString, "p_t_assoc", iBinPair[3]);
 580     axisTitlePair[3]    = "p_{T,assoc.} (GeV/c)";
 581
 582     dBinsPair[4]       = GetBinning(fBinningString, "delta_eta", iBinPair[4]);
 583     axisTitlePair[4]   = "#Delta#eta";
 584
 585     dBinsPair[5]       = GetBinning(fBinningString, "delta_phi", iBinPair[5]);
 586     axisTitlePair[5]   = "#Delta#varphi (rad)";
 587
 588  if(fcontainPIDtrig && !fcontainPIDasso){
 589     dBinsPair[6]       = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
 590     axisTitlePair[6]   = "PIDTrig";
 591  }
 592
 593  if(!fcontainPIDtrig && fcontainPIDasso){
 594     dBinsPair[6]       = GetBinning(fBinningString, "PIDAsso", iBinPair[6]);
 595     axisTitlePair[6]   = "PIDAsso";
 596  }
 597
 598 if(fcontainPIDtrig && fcontainPIDasso){
 599
 600     dBinsPair[6]       = GetBinning(fBinningString, "PIDTrig", iBinPair[6]);
 601     axisTitlePair[6]   = "PIDTrig";
 602
 603     dBinsPair[7]       = GetBinning(fBinningString, "PIDAsso", iBinPair[7]);
 604     axisTitlePair[7]   = "PIDAsso";
 605  }
 606
 607         Int_t nEtaBin = -1;
 608         Double_t* EtaBin = GetBinning(fBinningString, "eta", nEtaBin);
 609
 610         Int_t nPteffbin = -1;
 611         Double_t* Pteff = GetBinning(fBinningString, "p_t_eff", nPteffbin);
 612
 613
 614         fminPtTrig=dBinsPair[2][0];
 615         fmaxPtTrig=dBinsPair[2][iBinPair[2]];
 616         fminPtAsso=dBinsPair[3][0];
 617         fmaxPtAsso=dBinsPair[3][iBinPair[3]];
 618
 619         //fminPtComboeff=fminPtTrig;***then this value will be fixed ,even Setter can't change it's value
 620         //fmaxPtComboeff=fmaxPtTrig;
 621 //THnSparses for calculation of efficiency
 622
 623  if((fAnalysisType =="MCAOD") && ffillefficiency) {
 624      const Int_t nDim = 4;//       cent zvtx  pt   eta
 625      Int_t fBinsCh[nDim] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//**********************change it
 626      Double_t fMinCh[nDim] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
 627      Double_t fMaxCh[nDim] = { dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
 628
 629 TString Histrename;
 630 for(Int_t jj=0;jj<6;jj++)//PID type binning
 631     {
 632   Histrename="fTHnrecomatchedallPid";Histrename+=jj;
 633   fTHnrecomatchedallPid[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
 634   fTHnrecomatchedallPid[jj]->Sumw2();
 635   fTHnrecomatchedallPid[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
 636   fTHnrecomatchedallPid[jj]->GetAxis(0)->SetTitle("Centrality");
 637   fTHnrecomatchedallPid[jj]->GetAxis(1)->Set(iBinPair[1],dBinsPair[1]);
 638   fTHnrecomatchedallPid[jj]->GetAxis(1)->SetTitle("zvtx");
 639   fTHnrecomatchedallPid[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
 640   fTHnrecomatchedallPid[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 641   fTHnrecomatchedallPid[jj]->GetAxis(3)->Set(nEtaBin,EtaBin);
 642   fTHnrecomatchedallPid[jj]->GetAxis(3)->SetTitle("#eta");
 643   fOutput->Add(fTHnrecomatchedallPid[jj]);
 644
 645 Histrename="fTHngenprimPidTruth";Histrename+=jj;
 646   fTHngenprimPidTruth[jj] = new THnSparseF(Histrename.Data(),"cent:zvtx::Pt:eta", nDim, fBinsCh, fMinCh, fMaxCh);
 647   fTHngenprimPidTruth[jj]->Sumw2();
 648   fTHngenprimPidTruth[jj]->GetAxis(0)->Set(iBinPair[0],dBinsPair[0]);
 649   fTHngenprimPidTruth[jj]->GetAxis(0)->SetTitle("Centrality");
 650   fTHngenprimPidTruth[jj]->GetAxis(1)->Set(iBinPair[1], dBinsPair[1]);
 651   fTHngenprimPidTruth[jj]->GetAxis(1)->SetTitle("zvtx");
 652   fTHngenprimPidTruth[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
 653   fTHngenprimPidTruth[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 654   fTHngenprimPidTruth[jj]->GetAxis(3)->Set(nEtaBin, EtaBin);
 655   fTHngenprimPidTruth[jj]->GetAxis(3)->SetTitle("#eta");
 656   fOutput->Add(fTHngenprimPidTruth[jj]);
 657     }
 658  }
 659
 660  Int_t fBins[kTrackVariablesPair];
 661  Double_t fMin[kTrackVariablesPair];
 662  Double_t fMax[kTrackVariablesPair];
 663
 664 //ThnSparses for Correlation plots(data &  MC)
 665  for(Int_t i=0;i<kTrackVariablesPair;i++)
 666    {
 667      fBins[i] =iBinPair[i];
 668      fMin[i]= dBinsPair[i][0];
 669      fMax[i] = dBinsPair[i][iBinPair[i]];
 670    }
 671      if(ffilltrigassoUNID)
 672        {
 673   fTHnCorrUNID = new THnSparseF("fTHnCorrUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 674   fTHnCorrUNID->Sumw2();
 675   for(Int_t i=0; i<kTrackVariablesPair;i++){
 676     SetAsymmetricBin(fTHnCorrUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 677   }
 678   fOutput->Add(fTHnCorrUNID);
 679
 680 fTHnCorrUNIDmix = new THnSparseF("fTHnCorrUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 681   fTHnCorrUNIDmix->Sumw2();
 682 for(Int_t i=0; i<kTrackVariablesPair;i++){
 683     SetAsymmetricBin(fTHnCorrUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 684   }
 685   fOutput->Add(fTHnCorrUNIDmix);
 686        }
 687
 688      if(ffilltrigIDassoID)
 689        {
 690     fTHnCorrID = new THnSparseF("fTHnCorrID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 691   fTHnCorrID->Sumw2();
 692 for(Int_t i=0; i<kTrackVariablesPair;i++){
 693     SetAsymmetricBin(fTHnCorrID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 694   }
 695   fOutput->Add(fTHnCorrID);
 696
 697 fTHnCorrIDmix = new THnSparseF("fTHnCorrIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 698   fTHnCorrIDmix->Sumw2();
 699 for(Int_t i=0; i<kTrackVariablesPair;i++){
 700     SetAsymmetricBin(fTHnCorrIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 701   }
 702   fOutput->Add(fTHnCorrIDmix);
 703        }
 704
 705      if(ffilltrigUNIDassoID || ffilltrigIDassoUNID)//***********a bit tricky, be careful
 706        {
 707  fTHnCorrIDUNID = new THnSparseF("fTHnCorrIDUNID","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 708   fTHnCorrIDUNID->Sumw2();
 709 for(Int_t i=0; i<kTrackVariablesPair;i++){
 710     SetAsymmetricBin(fTHnCorrIDUNID,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 711   }
 712   fOutput->Add(fTHnCorrIDUNID);
 713
 714  fTHnCorrIDUNIDmix = new THnSparseF("fTHnCorrIDUNIDmix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 715   fTHnCorrIDUNIDmix->Sumw2();
 716 for(Int_t i=0; i<kTrackVariablesPair;i++){
 717     SetAsymmetricBin(fTHnCorrIDUNIDmix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 718   }
 719   fOutput->Add(fTHnCorrIDUNIDmix);
 720        }
 721
 722
 723
 724   //ThnSparse for Correlation plots(truth MC)
 725      if((fAnalysisType == "MCAOD") && ffilltrigIDassoIDMCTRUTH) {//remember that in this case uidentified means other than pions, kaons, protons
 726       fCorrelatonTruthPrimary = new THnSparseF("fCorrelatonTruthPrimary","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 727   fCorrelatonTruthPrimary->Sumw2();
 728 for(Int_t i=0; i<kTrackVariablesPair;i++){
 729     SetAsymmetricBin(fCorrelatonTruthPrimary,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 730   }
 731   fOutput->Add(fCorrelatonTruthPrimary);
 732
 733   fCorrelatonTruthPrimarymix = new THnSparseF("fCorrelatonTruthPrimarymix","cent:zvtx:pttrig:ptasso:deltaeta:deltaphi", kTrackVariablesPair, fBins, fMin, fMax);
 734   fCorrelatonTruthPrimarymix->Sumw2();
 735 for(Int_t i=0; i<kTrackVariablesPair;i++){
 736     SetAsymmetricBin(fCorrelatonTruthPrimarymix,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 737   }
 738   fOutput->Add(fCorrelatonTruthPrimarymix);
 739  }
 740
 741
 742     //binning for trigger no. counting
 743
 744         Double_t* fMint;
 745         Double_t* fMaxt;
 746         Int_t* fBinst;
 747         Int_t dims=3;
 748         if(fcontainPIDtrig) dims=4;
 749         fMint= new Double_t[dims];
 750         fMaxt= new Double_t[dims];
 751         fBinst= new Int_t[dims];
 752         for(Int_t i=0; i<3;i++)
 753           {
 754             fBinst[i]=iBinPair[i];
 755             fMint[i]=dBinsPair[i][0];
 756             fMaxt[i]=dBinsPair[i][iBinPair[i]];
 757           }
 758         if(fcontainPIDtrig){
 759             fBinst[3]=iBinPair[6];
 760             fMint[3]=dBinsPair[6][0];
 761             fMaxt[3]=dBinsPair[6][iBinPair[6]];
 762         }
 763
 764   //ThSparse for trigger counting(data & reco MC)
 765   if(ffilltrigassoUNID || ffilltrigUNIDassoID || ffilltrigIDassoUNID || ffilltrigIDassoID)
 766           {
 767   fTHnTrigcount = new THnSparseF("fTHnTrigcount","cent:zvtx:pt", dims, fBinst, fMint, fMaxt);
 768   fTHnTrigcount->Sumw2();
 769 for(Int_t i=0; i<3;i++){
 770     SetAsymmetricBin(fTHnTrigcount,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 771   }
 772  if(fcontainPIDtrig)   SetAsymmetricBin(fTHnTrigcount,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);
 773   fOutput->Add(fTHnTrigcount);
 774           }
 775
 776   if((fAnalysisType =="MCAOD") && ffilltrigIDassoIDMCTRUTH) {
 777   //ThSparse for trigger counting(truth MC)
 778 fTHnTrigcountMCTruthPrim = new THnSparseF("fTHnTrigcountMCTruthPrim","cent:zvtx:pt", dims, fBinst, fMint, fMaxt);
 779   fTHnTrigcountMCTruthPrim->Sumw2();
 780 for(Int_t i=0; i<3;i++){
 781     SetAsymmetricBin(fTHnTrigcountMCTruthPrim,i,dBinsPair[i],iBinPair[i],axisTitlePair[i]);
 782   }
 783 if(fcontainPIDtrig)    SetAsymmetricBin(fTHnTrigcountMCTruthPrim,3,dBinsPair[6],iBinPair[6],axisTitlePair[6]);
 784   fOutput->Add(fTHnTrigcountMCTruthPrim);
 785  }
 786
 787 if(fAnalysisType=="MCAOD"){
 788   if(ffillhistQATruth)
 789     {
 790   MCtruthpt=new TH1F ("MCtruthpt","ptdistributiontruthprim",100,0.,10.);
 791   fOutput->Add(MCtruthpt);
 792
 793   MCtrutheta=new TH1F ("MCtrutheta","etadistributiontruthprim",360,-1.8,1.8);
 794   fOutput->Add(MCtrutheta);
 795
 796   MCtruthphi=new TH1F ("MCtruthphi","phidisttruthprim",340,0,6.8);
 797   fOutput->Add(MCtruthphi);
 798
 799   MCtruthpionpt=new TH1F ("MCtruthpionpt","MCtruthpionpt",100,0.,10.);
 800   fOutput->Add(MCtruthpionpt);
 801
 802   MCtruthpioneta=new TH1F ("MCtruthpioneta","MCtruthpioneta",360,-1.8,1.8);
 803   fOutput->Add(MCtruthpioneta);
 804
 805   MCtruthpionphi=new TH1F ("MCtruthpionphi","MCtruthpionphi",340,0,6.8);
 806   fOutput->Add(MCtruthpionphi);
 807
 808   MCtruthkaonpt=new TH1F ("MCtruthkaonpt","MCtruthkaonpt",100,0.,10.);
 809   fOutput->Add(MCtruthkaonpt);
 810
 811   MCtruthkaoneta=new TH1F ("MCtruthkaoneta","MCtruthkaoneta",360,-1.8,1.8);
 812   fOutput->Add(MCtruthkaoneta);
 813
 814   MCtruthkaonphi=new TH1F ("MCtruthkaonphi","MCtruthkaonphi",340,0,6.8);
 815   fOutput->Add(MCtruthkaonphi);
 816
 817   MCtruthprotonpt=new TH1F ("MCtruthprotonpt","MCtruthprotonpt",100,0.,10.);
 818   fOutput->Add(MCtruthprotonpt);
 819
 820   MCtruthprotoneta=new TH1F ("MCtruthprotoneta","MCtruthprotoneta",360,-1.8,1.8);
 821   fOutput->Add(MCtruthprotoneta);
 822
 823   MCtruthprotonphi=new TH1F ("MCtruthprotonphi","MCtruthprotonphi",340,0,6.8);
 824   fOutput->Add(MCtruthprotonphi);
 825     }
 826  fPioncont=new TH2F("fPioncont", "fPioncont",10,-0.5,9.5,100,0.,10.);
 827   fOutput->Add(fPioncont);
 828
 829  fKaoncont=new TH2F("fKaoncont","fKaoncont",10,-0.5,9.5,100,0.,10.);
 830   fOutput->Add(fKaoncont);
 831
 832  fProtoncont=new TH2F("fProtoncont","fProtoncont",10,-0.5,9.5,100,0.,10.);
 833   fOutput->Add(fProtoncont);
 834   }
 835
 836 //Mixing
 837 DefineEventPool();
 838
 839   if(fapplyTrigefficiency || fapplyAssoefficiency)
 840    {
 841      const Int_t nDimt = 4;//       cent zvtx  pt   eta
 842      Int_t fBinsCht[nDimt] = {iBinPair[0], iBinPair[1], nPteffbin ,nEtaBin};//*************change it
 843      Double_t fMinCht[nDimt] = { dBinsPair[0][0],dBinsPair[1][0], Pteff[0], EtaBin[0] };
 844      Double_t fMaxCht[nDimt] = {dBinsPair[0][iBinPair[0]], dBinsPair[1][iBinPair[1]], Pteff[nPteffbin], EtaBin[nEtaBin]};
 845
 846   TString Histrexname;
 847 for(Int_t jj=0;jj<6;jj++)// PID type binning
 848     {
 849   Histrexname="effcorection";Histrexname+=jj;
 850   effcorection[jj] = new THnSparseF(Histrexname.Data(),"cent:zvtx::Pt:eta", nDimt, fBinsCht, fMinCht, fMaxCht);
 851   effcorection[jj]->Sumw2();
 852   effcorection[jj]->GetAxis(0)->Set(iBinPair[0], dBinsPair[0]);
 853   effcorection[jj]->GetAxis(0)->SetTitle("Centrality");
 854   effcorection[jj]->GetAxis(1)->Set( iBinPair[1],dBinsPair[1]);
 855   effcorection[jj]->GetAxis(1)->SetTitle("zvtx");
 856   effcorection[jj]->GetAxis(2)->Set(nPteffbin, Pteff);
 857   effcorection[jj]->GetAxis(2)->SetTitle("p_{T} (GeV/c)");
 858   effcorection[jj]->GetAxis(3)->Set( nEtaBin,EtaBin);
 859   effcorection[jj]->GetAxis(3)->SetTitle("#eta");
 860   fOutput->Add(effcorection[jj]);
 861     }
 862 // TFile *fsifile = new TFile(fefffilename,"READ");
 863
 864  if (TString(fefffilename).BeginsWith("alien:"))
 865     TGrid::Connect("alien:");
 866  TFile *fileT=TFile::Open(fefffilename);
 867  TString Nameg;
 868 for(Int_t jj=0;jj<6;jj++)//type binning
 869     {
 870 Nameg="effmap";Nameg+=jj;
 871 //effcorection[jj] = (THnSparseF*)fsifile->Get(Nameg.Data());
 872 effcorection[jj] = (THnSparseF*)fileT->Get(Nameg.Data());
 873
 874 //effcorection[jj]->SetDirectory(0);//****************************not present in case oh THnF
 875     }
 876 //fsifile->Close();
 877 fileT->Close();
 878
 879    }
 880
 881 //fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
 882 // fOutput->Add(fControlConvResoncances);
 883
 884
 885   PostData(1, fOutput);              // Post data for ALL output slots >0 here, to get at least an empty histogram
 886 }
 887 //-------------------------------------------------------------------------------
 888 void AliTwoParticlePIDCorr::UserExec( Option_t * ){
 889
 890
 891   if(fAnalysisType == "AOD") {
 892
 893     doAODevent();
 894
 895   }//AOD--analysis-----
 896
 897   else if(fAnalysisType == "MCAOD") {
 898
 899     doMCAODevent();
 900
 901   }
 902
 903   else return;
 904
 905 }
 906 //-------------------------------------------------------------------------
 907 void AliTwoParticlePIDCorr::doMCAODevent()
 908 {
 909   AliVEvent *event = InputEvent();
 910   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
 911   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
 912   if (!aod) {
 913     AliError("Cannot get the AOD event");
 914     return;
 915   }
 916
 917 // count all events(physics triggered)
 918   fEventCounter->Fill(1);
 919
 920  // get centrality object and check quality(valid for p-Pb and Pb-Pb)
 921   Double_t cent_v0=0.0;
 922
 923   if(fSampleType=="pPb" || fSampleType=="PbPb")
 924     {
 925   AliCentrality *centrality=0;
 926   if(aod)
 927   centrality = aod->GetHeader()->GetCentralityP();
 928   // if (centrality->GetQuality() != 0) return ;
 929
 930   if(centrality)
 931   {
 932   cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
 933   }
 934   else
 935     {
 936   cent_v0= -1;
 937      }
 938     }
 939
 940 //check the PIDResponse handler
 941   if (!fPID) return;
 942
 943 // get mag. field required for twotrack efficiency cut
 944  Float_t bSign = 0;
 945  bSign = (aod->GetMagneticField() > 0) ? 1 : -1;
 946
 947  //check for TClonesArray(truth track MC information)
 948 fArrayMC = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 949   if (!fArrayMC) {
 950     AliFatal("Error: MC particles branch not found!\n");
 951     return;
 952   }
 953
 954   //check for AliAODMCHeader(truth event MC information)
 955   AliAODMCHeader *header=NULL;
 956   header=(AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 957   if(!header) {
 958     printf("MC header branch not found!\n");
 959     return;
 960   }
 961
 962 //Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
 963 Float_t zVtxmc =header->GetVtxZ();
 964  if(TMath::Abs(zVtxmc)>fzvtxcut) return;
 965
 966  // For productions with injected signals, figure out above which label to skip particles/tracks
 967   Int_t skipParticlesAbove = 0;
 968
 969  if (fInjectedSignals)
 970   {
 971     AliGenEventHeader* eventHeader = 0;
 972     Int_t headers = 0;
 973
 974 // AOD
 975       if (!header)
 976       AliFatal("fInjectedSignals set but no MC header found");
 977
 978       headers = header->GetNCocktailHeaders();
 979       eventHeader = header->GetCocktailHeader(0);
 980
 981  if (!eventHeader)
 982     {
 983       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
 984       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
 985       AliError("First event header not found. Skipping this event.");
 986       //fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
 987       return;
 988     }
 989 skipParticlesAbove = eventHeader->NProduced();
 990     AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
 991   }
 992
 993  //vertex selection(is it fine for PP?)
 994   if ( fSampleType=="pPb"){
 995   trkVtx = aod->GetPrimaryVertex();
 996   if (!trkVtx || trkVtx->GetNContributors()<=0) return;
 997   TString vtxTtl = trkVtx->GetTitle();
 998   if (!vtxTtl.Contains("VertexerTracks")) return;
 999    zvtx = trkVtx->GetZ();
1000   const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
1001   if (!spdVtx || spdVtx->GetNContributors()<=0) return;
1002   TString vtxTyp = spdVtx->GetTitle();
1003   Double_t cov[6]={0};
1004   spdVtx->GetCovarianceMatrix(cov);
1005   Double_t zRes = TMath::Sqrt(cov[5]);
1006   if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1007    if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1008   }
1009   else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case
1010         Int_t nVertex = aod->GetNumberOfVertices();
1011         if( nVertex > 0 ) {
1012     trkVtx = aod->GetPrimaryVertex();
1013                 Int_t nTracksPrim = trkVtx->GetNContributors();
1014                  zvtx = trkVtx->GetZ();
1015                 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1016                 // Reject TPC only vertex
1017                 TString name(trkVtx->GetName());
1018                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
1019
1020                 // Select a quality vertex by number of tracks?
1021                 if( nTracksPrim < fnTracksVertex ) {
1022                   //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1023                         return ;
1024                         }
1025                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1026                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1027                 //  return kFALSE;
1028                 //      if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1029         }
1030         else return;
1031
1032   }
1033   else return;//as there is no proper sample type
1034
1035
1036   fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ()));   //for trkVtx only before vertex cut |zvtx|<10 cm
1037
1038   //count events having a proper vertex
1039    fEventCounter->Fill(2);
1040
1041  if (TMath::Abs(zvtx) > fzvtxcut) return;
1042
1043 fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut for trkVtx only
1044
1045 //now we have events passed physics trigger, centrality,eventzvtx cut
1046
1047 //count events after vertex cut
1048   fEventCounter->Fill(5);
1049
1050 if(!aod) return;  //for safety
1051
1052 if (fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0)  return;//for pp case it is the multiplicity
1053
1054    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)
1055
1056    Double_t cent_v0_truth=0.0;
1057
1058    //TObjArray* tracksMCtruth=0;
1059 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
1060  tracksMCtruth->SetOwner(kTRUE);  //***********************************IMPORTANT!
1061
1062   eventno++;
1063
1064   //There is a small difference between zvtx and zVtxmc??????
1065   //cout<<"***********************************************zVtxmc="<<zVtxmc<<endl;
1066   //cout<<"***********************************************zvtx="<<zvtx<<endl;
1067
1068 //now process the truth particles(for both efficiency & correlation function)
1069 Int_t nMCTrack = fArrayMC->GetEntriesFast();
1070
1071 for (Int_t iMC = 0; iMC < nMCTrack; iMC++)
1072 {      //MC truth track loop starts
1073
1074 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
1075
1076     if(!partMC){
1077       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
1078       continue;
1079     }
1080
1081 //consider only charged particles
1082     if(partMC->Charge() == 0) continue;
1083
1084 //consider only primary particles; neglect all secondary particles including from weak decays
1085  if(fselectprimaryTruth && !partMC->IsPhysicalPrimary()) continue;
1086
1087
1088 //remove injected signals(primaries above <maxLabel>)
1089  if (fInjectedSignals && partMC->GetLabel() >= skipParticlesAbove) continue;
1090
1091 //remove duplicates
1092   Bool_t isduplicate=kFALSE;
1093  if (fRemoveDuplicates)
1094    {
1095  for (Int_t j=iMC+1; j<nMCTrack; ++j)
1096    {//2nd trutuh loop starts
1097 AliAODMCParticle *partMC2 = (AliAODMCParticle*) fArrayMC->At(j);
1098    if(!partMC2){
1099       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",j));
1100       continue;
1101     }
1102  if (partMC->GetLabel() == partMC2->GetLabel())
1103    {
1104 isduplicate=kTRUE;
1105  break;
1106    }
1107    }//2nd truth loop ends
1108    }
1109  if(fRemoveDuplicates && isduplicate) continue;//remove duplicates
1110
1111 //give only kinematic cuts at the truth level
1112  if (partMC->Eta() < fmineta || partMC->Eta() > fmaxeta) continue;
1113  if (partMC->Pt() < fminPt ||  partMC->Pt() > fmaxPt) continue;
1114
1115  if(!partMC) continue;//for safety
1116
1117  //To determine multiplicity in case of PP
1118  nooftrackstruth++;
1119  //cout<<"**************************************"<<TMath::Abs(partMC->GetLabel())<<endl;
1120 //only physical primary(all/unidentified)
1121 if(ffillhistQATruth)
1122     {
1123  MCtruthpt->Fill(partMC->Pt());
1124  MCtrutheta->Fill(partMC->Eta());
1125  MCtruthphi->Fill(partMC->Phi());
1126     }
1127  //get particle ID
1128 Int_t pdgtruth=((AliAODMCParticle*)partMC)->GetPdgCode();
1129 Int_t particletypeTruth=-999;
1130  if (TMath::Abs(pdgtruth)==211)
1131    {
1132  particletypeTruth=SpPion;
1133 if(ffillhistQATruth)
1134     {
1135  MCtruthpionpt->Fill(partMC->Pt());
1136  MCtruthpioneta->Fill(partMC->Eta());
1137  MCtruthpionphi->Fill(partMC->Phi());
1138     }
1139       }
1140  if (TMath::Abs(pdgtruth)==321)
1141    {
1142  particletypeTruth=SpKaon;
1143 if(ffillhistQATruth)
1144     {
1145  MCtruthkaonpt->Fill(partMC->Pt());
1146  MCtruthkaoneta->Fill(partMC->Eta());
1147  MCtruthkaonphi->Fill(partMC->Phi());
1148   }
1149     }
1150 if(TMath::Abs(pdgtruth)==2212)
1151   {
1152  particletypeTruth=SpProton;
1153 if(ffillhistQATruth)
1154     {
1155  MCtruthprotonpt->Fill(partMC->Pt());
1156  MCtruthprotoneta->Fill(partMC->Eta());
1157  MCtruthprotonphi->Fill(partMC->Phi());
1158     }
1159      }
1160  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)
1161
1162  // -- Fill THnSparse for efficiency and contamination calculation
1163  if (fSampleType=="pp") cent_v0=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
1164
1165  Double_t primmctruth[4] = {cent_v0, zVtxmc,partMC->Pt(), partMC->Eta()};
1166  if(ffillefficiency)
1167   {
1168   fTHngenprimPidTruth[5]->Fill(primmctruth);//for all primary truth particles(4)
1169 if (TMath::Abs(pdgtruth)==211 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[3]->Fill(primmctruth);//for  primary truth mesons(3)
1170 if (TMath::Abs(pdgtruth)==2212 || TMath::Abs(pdgtruth)==321) fTHngenprimPidTruth[4]->Fill(primmctruth);//for  primary truth kaons+protons(4)
1171   if (TMath::Abs(pdgtruth)==211)  fTHngenprimPidTruth[0]->Fill(primmctruth);//for pions
1172   if (TMath::Abs(pdgtruth)==321)  fTHngenprimPidTruth[1]->Fill(primmctruth);//for kaons
1173   if(TMath::Abs(pdgtruth)==2212)  fTHngenprimPidTruth[2]->Fill(primmctruth);//for protons
1174   }
1175
1176  Float_t effmatrixtruth=1.0;//In Truth MC, no case of efficiency correction so it should be always 1.0
1177 if(partMC->Pt()>=fminPtAsso || partMC->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
1178   {
1179 LRCParticlePID* copy6 = new LRCParticlePID(particletypeTruth,partMC->Charge(),partMC->Pt(),partMC->Eta(), partMC->Phi(),effmatrixtruth);
1180 //copy6->SetUniqueID(eventno * 100000 + TMath::Abs(partMC->GetLabel()));
1181  copy6->SetUniqueID(eventno * 100000 + (Int_t)nooftrackstruth);
1182  tracksMCtruth->Add(copy6);//************** TObjArray used for truth correlation function calculation
1183   }
1184   }//MC truth track loop ends
1185
1186 //*********************still in event loop
1187
1188  if (fSampleType=="pp") cent_v0_truth=nooftrackstruth;
1189  else cent_v0_truth=cent_v0;//the notation cent_v0 is reserved for reco/data case
1190
1191  //now cent_v0_truth should be used for all correlation function calculation
1192
1193 if(nooftrackstruth>0.0 && ffilltrigIDassoIDMCTRUTH)
1194   {
1195  //Fill Correlations for MC truth particles(same event)
1196 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)//hadron triggered correlation
1197  Fillcorrelation(tracksMCtruth,0,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kFALSE,"trigIDassoIDMCTRUTH");//mixcase=kFALSE for same event case
1198
1199 //start mixing
1200 AliEventPool* pool2 = fPoolMgr->GetEventPool(cent_v0_truth, zVtxmc+200);
1201 if (pool2 && pool2->IsReady())
1202   {//start mixing only when pool->IsReady
1203 if(tracksMCtruth && tracksMCtruth->GetEntriesFast()>0)
1204   {//proceed only when no. of trigger particles >0 in current event
1205 for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
1206   { //pool event loop start
1207  TObjArray* bgTracks6 = pool2->GetEvent(jMix);
1208   if(!bgTracks6) continue;
1209   Fillcorrelation(tracksMCtruth,bgTracks6,cent_v0_truth,zVtxmc,bSign,fPtOrderMCTruth,kFALSE,kTRUE,"trigIDassoIDMCTRUTH");//mixcase=kTRUE for mixing case
1210
1211    }// pool event loop ends mixing case
1212  }//if(trackstrig && trackstrig->GetEntriesFast()>0) condition ends mixing case
1213 } //if pool->IsReady() condition ends mixing case
1214
1215  //still in main event loop
1216
1217  if(tracksMCtruth){
1218 if(pool2)  pool2->UpdatePool(CloneAndReduceTrackList(tracksMCtruth));//ownership of tracksasso is with pool now, don't delete it
1219  }
1220   }
1221
1222  //still in main event loop
1223
1224 if(tracksMCtruth) delete tracksMCtruth;
1225
1226 //now deal with reco tracks
1227    //TObjArray* tracksUNID=0;
1228    TObjArray* tracksUNID = new TObjArray;
1229    tracksUNID->SetOwner(kTRUE);
1230
1231    //TObjArray* tracksID=0;
1232    TObjArray* tracksID = new TObjArray;
1233    tracksID->SetOwner(kTRUE);
1234
1235
1236    Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//used for reconstructed track dca cut
1237
1238    Double_t trackscount=0.0;
1239
1240 // loop over reconstructed tracks
1241   for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1242 { // reconstructed track loop starts
1243   AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1244   if (!track) continue;
1245  //get the corresponding MC track at the truth level (doing reco matching)
1246   AliAODMCParticle* recomatched = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track->GetLabel())));
1247   if(!recomatched) continue;//if a reco track doesn't have corresponding truth track at generated level is a fake track(label==0), ignore it
1248
1249 //remove injected signals
1250  if(fInjectedSignals)
1251    {
1252     AliAODMCParticle* mother = recomatched;
1253
1254       while (!mother->IsPhysicalPrimary())
1255       {// find the primary mother;the first stable mother is searched and checked if it is <= <maxLabel>
1256         if (mother->GetMother() < 0)
1257         {
1258           mother = 0;
1259           break;
1260         }
1261
1262    mother =(AliAODMCParticle*) fArrayMC->At(((AliAODMCParticle*)mother)->GetMother());
1263         if (!mother)
1264           break;
1265       }
1266  if (!mother)
1267     {
1268       Printf("WARNING: No mother found for particle %d:", recomatched->GetLabel());
1269       continue;
1270     }
1271  if (mother->GetLabel() >= skipParticlesAbove) continue;//remove injected signals(primaries above <maxLabel>)
1272    }//remove injected signals
1273
1274  if (fRemoveWeakDecays && ((AliAODMCParticle*) recomatched)->IsSecondaryFromWeakDecay()) continue;//remove weak decays
1275
1276   Bool_t isduplicate2=kFALSE;
1277 if (fRemoveDuplicates)
1278    {
1279   for (Int_t j =itrk+1; j < aod->GetNumberOfTracks(); j++)
1280     {//2nd loop starts
1281  AliAODTrack* track2 = dynamic_cast<AliAODTrack*>(aod->GetTrack(j));
1282  if (!track2) continue;
1283  AliAODMCParticle* recomatched2 = static_cast<AliAODMCParticle*>(fArrayMC->At(TMath::Abs(track2->GetLabel())));
1284 if(!recomatched2) continue;
1285
1286 if (track->GetLabel() == track2->GetLabel())
1287    {
1288 isduplicate2=kTRUE;
1289  break;
1290    }
1291     }//2nd loop ends
1292    }
1293  if(fRemoveDuplicates && isduplicate2) continue;//remove duplicates
1294
1295   fHistQA[11]->Fill(track->GetTPCNcls());
1296   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
1297
1298  if(tracktype==0) continue;
1299  if(tracktype==1)//tracks "not" passed AliAODTrack::kPrimary at reconstructed level & have proper TPC PID response(?)
1300 {
1301   if(!track) continue;//for safety
1302  //accepted all(primaries+secondary) reconstructed tracks(pt 0.2 to 10.0,,eta -0.8 to 0.8)
1303   trackscount++;
1304
1305 //check for eta , phi holes
1306  fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1307  fphiSpectraasso->Fill(track->Phi(),track->Pt());
1308
1309   Int_t particletypeMC=-9999;
1310
1311 //tag all particles as unidentified
1312  particletypeMC=unidentified;
1313
1314  Float_t effmatrix=1.;
1315
1316 // -- Fill THnSparse for efficiency calculation
1317  if (fSampleType=="pp") cent_v0=15.0;//integrated over multiplicity(so put any fixed value for each track so that practically means there is only one bin in multiplicityi.e multiplicity intregated out )**************Important
1318  //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)
1319
1320  //Clone & Reduce track list(TObjArray) for unidentified particles
1321 if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
1322   {
1323  if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
1324    effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);
1325    LRCParticlePID* copy = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1326    copy->SetUniqueID(eventno * 100000 +(Int_t)trackscount);
1327    tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
1328   }
1329  Double_t allrecomatchedpid[4] = {cent_v0, zVtxmc,recomatched->Pt(), recomatched->Eta()};
1330 if(ffillefficiency) fTHnrecomatchedallPid[5]->Fill(allrecomatchedpid);//for all
1331
1332  //now start the particle identification process:)
1333
1334 //get the pdg code of the corresponding truth particle
1335  Int_t pdgCode = ((AliAODMCParticle*)recomatched)->GetPdgCode();
1336
1337 Float_t dEdx = track->GetTPCsignal();
1338  fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1339
1340  if(HasTOFPID(track))
1341 {
1342 Float_t beta = GetBeta(track);
1343 fHistoTOFbeta->Fill(track->Pt(), beta);
1344  }
1345
1346  //do track identification(nsigma method)
1347   particletypeMC=GetParticle(track);//******************************problem is here
1348
1349 //2-d TPCTOF map(for each Pt interval)
1350   if(HasTOFPID(track)){
1351  fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
1352  fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
1353  fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
1354   }
1355  if(particletypeMC==SpUndefined) continue;
1356
1357  //Pt, Eta , Phi distribution of the reconstructed identified particles
1358 if(ffillhistQAReco)
1359     {
1360 if (particletypeMC==SpPion)
1361   {
1362     fPionPt->Fill(track->Pt());
1363     fPionEta->Fill(track->Eta());
1364     fPionPhi->Fill(track->Phi());
1365   }
1366 if (particletypeMC==SpKaon)
1367   {
1368     fKaonPt->Fill(track->Pt());
1369     fKaonEta->Fill(track->Eta());
1370     fKaonPhi->Fill(track->Phi());
1371   }
1372 if (particletypeMC==SpProton)
1373   {
1374     fProtonPt->Fill(track->Pt());
1375     fProtonEta->Fill(track->Eta());
1376     fProtonPhi->Fill(track->Phi());
1377   }
1378     }
1379
1380  //for misidentification fraction calculation(do it with tuneonPID)
1381  if(particletypeMC==SpPion )
1382    {
1383      if(TMath::Abs(pdgCode)==211) fPioncont->Fill(1.,track->Pt());
1384      if(TMath::Abs(pdgCode)==321) fPioncont->Fill(3.,track->Pt());
1385      if(TMath::Abs(pdgCode)==2212) fPioncont->Fill(5.,track->Pt());
1386      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fPioncont->Fill(7.,track->Pt());
1387    }
1388 if(particletypeMC==SpKaon )
1389    {
1390      if(TMath::Abs(pdgCode)==211) fKaoncont->Fill(1.,track->Pt());
1391      if(TMath::Abs(pdgCode)==321) fKaoncont->Fill(3.,track->Pt());
1392      if(TMath::Abs(pdgCode)==2212) fKaoncont->Fill(5.,track->Pt());
1393      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fKaoncont->Fill(7.,track->Pt());
1394    }
1395  if(particletypeMC==SpProton )
1396    {
1397      if(TMath::Abs(pdgCode)==211) fProtoncont->Fill(1.,track->Pt());
1398      if(TMath::Abs(pdgCode)==321) fProtoncont->Fill(3.,track->Pt());
1399      if(TMath::Abs(pdgCode)==2212) fProtoncont->Fill(5.,track->Pt());
1400      if(TMath::Abs(pdgCode)!=211 && TMath::Abs(pdgCode)!=321 && TMath::Abs(pdgCode)!=2212) fProtoncont->Fill(7.,track->Pt());
1401    }
1402
1403  //fill tracking efficiency
1404  if(ffillefficiency)
1405    {
1406  if(particletypeMC==SpPion || particletypeMC==SpKaon)
1407    {
1408 if(TMath::Abs(pdgCode)==211 ||  TMath::Abs(pdgCode)==321)   fTHnrecomatchedallPid[3]->Fill(allrecomatchedpid);//for mesons
1409    }
1410  if(particletypeMC==SpKaon || particletypeMC==SpProton)
1411    {
1412 if(TMath::Abs(pdgCode)==321 ||  TMath::Abs(pdgCode)==2212)   fTHnrecomatchedallPid[4]->Fill(allrecomatchedpid);//for kaons+protons
1413    }
1414   if(particletypeMC==SpPion && TMath::Abs(pdgCode)==211)  fTHnrecomatchedallPid[0]->Fill(allrecomatchedpid);//for pions
1415   if(particletypeMC==SpKaon && TMath::Abs(pdgCode)==321)  fTHnrecomatchedallPid[1]->Fill(allrecomatchedpid);//for kaons
1416   if(particletypeMC==SpProton && TMath::Abs(pdgCode)==2212) fTHnrecomatchedallPid[2]->Fill(allrecomatchedpid);//for protons
1417    }
1418
1419 if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)//to reduce memory consumption in pool
1420   {
1421 if (fapplyTrigefficiency || fapplyAssoefficiency)
1422     effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletypeMC);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles
1423     LRCParticlePID* copy1 = new LRCParticlePID(particletypeMC,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1424     copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1425     tracksID->Add(copy1);
1426   }
1427   }// if(tracktype==1) condition structure ands
1428
1429 }//reco track loop ends
1430
1431   //*************************************************************still in event loop
1432
1433 //same event delta-eta-deltaphi plot
1434 if(fSampleType=="pp")  cent_v0=trackscount;//multiplicity
1435
1436 if(trackscount>0.0)
1437   {
1438 //fill the centrality/multiplicity distribution of the selected events
1439  fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
1440
1441  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??????
1442
1443 //count selected events having centrality betn 0-100%
1444  fEventCounter->Fill(6);
1445
1446  //same event delte-eta, delta-phi plot
1447 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
1448   {//same event calculation starts
1449     if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
1450     if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID)  Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
1451   }
1452
1453 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
1454   {//same event calculation starts
1455     if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID)  Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
1456     if(ffilltrigIDassoID)   Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
1457   }
1458
1459 //still in  main event loop
1460 //start mixing
1461  if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
1462 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
1463 if (pool && pool->IsReady())
1464   {//start mixing only when pool->IsReady
1465  for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
1466   { //pool event loop start
1467  TObjArray* bgTracks = pool->GetEvent(jMix);
1468   if(!bgTracks) continue;
1469   if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
1470   Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
1471  if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
1472   Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
1473    }// pool event loop ends mixing case
1474
1475 } //if pool->IsReady() condition ends mixing case
1476  if(tracksUNID) {
1477 if(pool)
1478   pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
1479  }
1480  }//mixing with unidentified particles
1481
1482  if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
1483 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
1484 if (pool1 && pool1->IsReady())
1485   {//start mixing only when pool->IsReady
1486 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
1487   { //pool event loop start
1488  TObjArray* bgTracks2 = pool1->GetEvent(jMix);
1489   if(!bgTracks2) continue;
1490 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
1491   Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
1492 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
1493   Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
1494
1495    }// pool event loop ends mixing case
1496 } //if pool1->IsReady() condition ends mixing case
1497
1498 if(tracksID) {
1499 if(pool1)
1500   pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
1501  }
1502  }//mixing with identified particles
1503
1504   //no. of events analyzed
1505 fEventCounter->Fill(7);
1506   }
1507
1508 if(tracksUNID)  delete tracksUNID;
1509
1510 if(tracksID) delete tracksID;
1511
1512
1513 PostData(1, fOutput);
1514
1515 }
1516 //________________________________________________________________________
1517 void AliTwoParticlePIDCorr::doAODevent()
1518 {
1519
1520   //get AOD
1521   AliVEvent *event = InputEvent();
1522   if (!event) { Printf("ERROR: Could not retrieve event"); return; }
1523   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
1524   if (!aod) {
1525     AliError("Cannot get the AOD event");
1526     return;
1527   }
1528
1529 // count all events
1530   fEventCounter->Fill(1);
1531
1532   // get centrality object and check quality
1533   Double_t cent_v0=0;
1534
1535   if(fSampleType=="pPb" || fSampleType=="PbPb")
1536     {
1537   AliCentrality *centrality=0;
1538   if(aod)
1539   centrality = aod->GetHeader()->GetCentralityP();
1540   // if (centrality->GetQuality() != 0) return ;
1541
1542   if(centrality)
1543   {
1544   cent_v0 = centrality->GetCentralityPercentile(fCentralityMethod);
1545   }
1546   else
1547     {
1548   cent_v0= -1;
1549      }
1550     }
1551
1552  Float_t bSign = (aod->GetMagneticField() > 0) ? 1 : -1;//for two track efficiency cut in correlation function calculation
1553  Float_t bSign1=aod->GetHeader()->GetMagneticField() ;//for dca cut in ClassifyTrack(), i.e in track loop
1554
1555 // Pileup selection ************************************************
1556 // if(frejectPileUp)  //applicable only for TPC only tracks,not for hybrid tracks?.
1557      // {
1558  //if (fAnalysisUtils && fAnalysisUtils->IsPileUpEvent(aod)) return;
1559  // }
1560
1561  if (!fPID) return;//this should be available with each event even if we don't do PID selection
1562
1563   //vertex selection(is it fine for PP?)
1564   if ( fSampleType=="pPb"){
1565   trkVtx = aod->GetPrimaryVertex();
1566   if (!trkVtx || trkVtx->GetNContributors()<=0) return;
1567   TString vtxTtl = trkVtx->GetTitle();
1568   if (!vtxTtl.Contains("VertexerTracks")) return;
1569    zvtx = trkVtx->GetZ();
1570   const AliAODVertex* spdVtx = aod->GetPrimaryVertexSPD();
1571   if (!spdVtx || spdVtx->GetNContributors()<=0) return;
1572   TString vtxTyp = spdVtx->GetTitle();
1573   Double_t cov[6]={0};
1574   spdVtx->GetCovarianceMatrix(cov);
1575   Double_t zRes = TMath::Sqrt(cov[5]);
1576   if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) return;
1577    if (TMath::Abs(spdVtx->GetZ() - trkVtx->GetZ())>0.5) return;
1578   }
1579   else if(fSampleType=="PbPb" || fSampleType=="pp") {//for pp and pb-pb case
1580         Int_t nVertex = aod->GetNumberOfVertices();
1581         if( nVertex > 0 ) {
1582      trkVtx = aod->GetPrimaryVertex();
1583                 Int_t nTracksPrim = trkVtx->GetNContributors();
1584                  zvtx = trkVtx->GetZ();
1585                 //if (fDebug > 1)AliInfo(Form(" Vertex in = %f with %d particles by  %s data ...",zVertex,nTracksPrim,vertex->GetName()));
1586                 // Reject TPC only vertex
1587                 TString name(trkVtx->GetName());
1588                 if (name.CompareTo("PrimaryVertex") && name.CompareTo("SPDVertex"))return;
1589
1590                 // Select a quality vertex by number of tracks?
1591                 if( nTracksPrim < fnTracksVertex ) {
1592                   //if (fDebug > 1) AliInfo(" Primary-vertex Selection: event REJECTED ...");
1593                         return ;
1594                         }
1595                 // TODO remove vertexer Z events with dispersion > 0.02: Doesn't work for AOD at present
1596                 //if (strcmp(vertex->GetTitle(), "AliVertexerZ") == 0 && vertex->GetDispersion() > 0.02)
1597                 //  return kFALSE;
1598                 //      if (fDebug > 1) AliInfo(" Primary-vertex Selection: event ACCEPTED...");
1599         }
1600         else return;
1601
1602   }
1603   else return;//as there is no proper sample type
1604
1605
1606   fHistQA[0]->Fill((trkVtx->GetX()));fHistQA[1]->Fill((trkVtx->GetY()));fHistQA[2]->Fill((trkVtx->GetZ())); //for trkVtx only before vertex cut |zvtx|<10 cm
1607
1608 //count events having a proper vertex
1609    fEventCounter->Fill(2);
1610
1611  if (TMath::Abs(zvtx) > fzvtxcut) return;
1612
1613 //count events after vertex cut
1614   fEventCounter->Fill(5);
1615
1616
1617  //if(!fAnalysisUtils->IsVertexSelected2013pA(aod)) return;
1618
1619  fHistQA[3]->Fill((trkVtx->GetX()));fHistQA[4]->Fill((trkVtx->GetY()));fHistQA[5]->Fill((trkVtx->GetZ()));//after vertex cut,for trkVtx only
1620
1621
1622  if(!aod) return; //for safety
1623
1624 if(fSampleType=="pPb" || fSampleType=="PbPb") if (cent_v0 < 0)  return;//for pp case it is the multiplicity
1625
1626    TObjArray*  tracksUNID= new TObjArray;//track info before doing PID
1627    tracksUNID->SetOwner(kTRUE);  // IMPORTANT!
1628
1629    TObjArray* tracksID= new TObjArray;//only pions, kaons,protonsi.e. after doing the PID selection
1630    tracksID->SetOwner(kTRUE);  // IMPORTANT!
1631
1632      Double_t trackscount=0.0;//in case of pp this will give the multiplicity after the track loop(only for unidentified particles that pass the filterbit and kinematic cuts)
1633
1634     eventno++;
1635
1636  for (Int_t itrk = 0; itrk < aod->GetNumberOfTracks(); itrk++)
1637 { //track loop starts for TObjArray(containing track and event information) filling; used for correlation function calculation
1638   AliAODTrack* track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrk));
1639   if (!track) continue;
1640   fHistQA[11]->Fill(track->GetTPCNcls());
1641   Int_t particletype=-9999;//required for PID filling
1642   Int_t tracktype=ClassifyTrack(track,trkVtx,bSign1);//dcacut=kFALSE,onlyprimary=kFALSE
1643   if(tracktype!=1) continue;
1644
1645   if(!track) continue;//for safety
1646
1647 //check for eta , phi holes
1648  fEtaSpectrasso->Fill(track->Eta(),track->Pt());
1649  fphiSpectraasso->Fill(track->Phi(),track->Pt());
1650
1651  trackscount++;
1652
1653  //if no applyefficiency , set the eff factor=1.0
1654  Float_t effmatrix=1.0;
1655
1656  //tag all particles as unidentified that passed filterbit & kinematic cuts
1657  particletype=unidentified;
1658
1659
1660  if (fSampleType=="pp") cent_v0=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
1661
1662
1663  //to reduce memory consumption in pool
1664   if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)
1665   {
1666  //Clone & Reduce track list(TObjArray) for unidentified particles
1667  if (fapplyTrigefficiency || fapplyAssoefficiency)//get the trackingefficiency x contamination factor for unidentified particles
1668    effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);
1669  LRCParticlePID* copy = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1670   copy->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1671   tracksUNID->Add(copy);//track information Storage for UNID correlation function(tracks that pass the filterbit & kinematic cuts only)
1672   }
1673
1674 //now start the particle identificaion process:)
1675
1676 //track passing filterbit 768 have proper TPC response,or need to be checked explicitly before doing PID????
1677
1678   Float_t dEdx = track->GetTPCsignal();
1679   fHistoTPCdEdx->Fill(track->Pt(), dEdx);
1680
1681   //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)
1682  if(HasTOFPID(track))
1683 {
1684   Float_t beta = GetBeta(track);
1685   fHistoTOFbeta->Fill(track->Pt(), beta);
1686  }
1687
1688
1689 //track identification(using nsigma method)
1690   particletype=GetParticle(track);//*******************************change may be required(It should return only pion,kaon, proton and Spundefined; NOT unidentifed***************be careful)
1691
1692
1693 //2-d TPCTOF map(for each Pt interval)
1694   if(HasTOFPID(track)){
1695  fTPCTOFPion3d->Fill(track->Pt(),fnsigmas[SpPion][NSigmaTOF],fnsigmas[SpPion][NSigmaTPC]);
1696  fTPCTOFKaon3d->Fill(track->Pt(),fnsigmas[SpKaon][NSigmaTOF],fnsigmas[SpKaon][NSigmaTPC]);
1697  fTPCTOFProton3d->Fill(track->Pt(),fnsigmas[SpProton][NSigmaTOF],fnsigmas[SpProton][NSigmaTPC]);
1698   }
1699
1700 //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!!!!!
1701   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)!!!!!!!!!!!
1702
1703  //Pt, Eta , Phi distribution of the reconstructed identified particles
1704 if(ffillhistQAReco)
1705     {
1706 if (particletype==SpPion)
1707   {
1708     fPionPt->Fill(track->Pt());
1709     fPionEta->Fill(track->Eta());
1710     fPionPhi->Fill(track->Phi());
1711   }
1712 if (particletype==SpKaon)
1713   {
1714     fKaonPt->Fill(track->Pt());
1715     fKaonEta->Fill(track->Eta());
1716     fKaonPhi->Fill(track->Phi());
1717   }
1718 if (particletype==SpProton)
1719   {
1720     fProtonPt->Fill(track->Pt());
1721     fProtonEta->Fill(track->Eta());
1722     fProtonPhi->Fill(track->Phi());
1723   }
1724     }
1725
1726 if(track->Pt()>=fminPtAsso || track->Pt()<=fmaxPtTrig)
1727   {
1728 if (fapplyTrigefficiency || fapplyAssoefficiency)
1729   effmatrix=GetTrackbyTrackeffvalue(track,cent_v0,zvtx,particletype);//get the tracking eff x TOF matching eff x PID eff x contamination factor for identified particles; Bool_t mesoneffrequired=kFALSE
1730  LRCParticlePID* copy1 = new LRCParticlePID(particletype,track->Charge(),track->Pt(),track->Eta(), track->Phi(),effmatrix);
1731     copy1->SetUniqueID(eventno * 100000 + (Int_t)trackscount);
1732     tracksID->Add(copy1);
1733   }
1734 } //track loop ends but still in event loop
1735
1736 if(trackscount<1.0){
1737   if(tracksUNID) delete tracksUNID;
1738   if(tracksID) delete tracksID;
1739   return;
1740  }
1741
1742 // cout<<fminPtAsso<<"***"<<fmaxPtAsso<<"*********************"<<fminPtTrig<<"***"<<fmaxPtTrig<<"*****************"<<kTrackVariablesPair<<endl;
1743 if(fSampleType=="pp")  cent_v0=trackscount;//multiplicity
1744
1745 //fill the centrality/multiplicity distribution of the selected events
1746  fhistcentrality->Fill(cent_v0);//*********************************WARNING::binning of cent_v0 is different for pp and pPb/PbPb case
1747
1748 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??????
1749
1750 //count selected events having centrality betn 0-100%
1751  fEventCounter->Fill(6);
1752
1753 //same event delta-eta-deltaphi plot
1754
1755 if(tracksUNID && tracksUNID->GetEntriesFast()>0)//hadron triggered correlation
1756   {//same event calculation starts
1757     if(ffilltrigassoUNID) Fillcorrelation(tracksUNID,0,cent_v0,zvtx,bSign,kTRUE,kTRUE,kFALSE,"trigassoUNID");//mixcase=kFALSE (hadron-hadron correlation)
1758     if(tracksID && tracksID->GetEntriesFast()>0 && ffilltrigUNIDassoID)  Fillcorrelation(tracksUNID,tracksID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigUNIDassoID");//mixcase=kFALSE (hadron-ID correlation)
1759   }
1760
1761 if(tracksID && tracksID->GetEntriesFast()>0)//ID triggered correlation
1762   {//same event calculation starts
1763     if(tracksUNID && tracksUNID->GetEntriesFast()>0 && ffilltrigIDassoUNID)  Fillcorrelation(tracksID,tracksUNID,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoUNID");//mixcase=kFALSE (ID-hadron correlation)
1764     if(ffilltrigIDassoID)   Fillcorrelation(tracksID,0,cent_v0,zvtx,bSign,kFALSE,kTRUE,kFALSE,"trigIDassoID");//mixcase=kFALSE (ID-ID correlation)
1765   }
1766
1767 //still in  main event loop
1768
1769
1770 //start mixing
1771  if(ffilltrigassoUNID || ffilltrigIDassoUNID){//mixing with unidentified particles
1772 AliEventPool* pool = fPoolMgr->GetEventPool(cent_v0, zvtx);//In the pool there is tracksUNID(i.e associateds are unidentified)
1773 if (pool && pool->IsReady())
1774   {//start mixing only when pool->IsReady
1775  for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
1776   { //pool event loop start
1777  TObjArray* bgTracks = pool->GetEvent(jMix);
1778   if(!bgTracks) continue;
1779   if(ffilltrigassoUNID && tracksUNID && tracksUNID->GetEntriesFast()>0)//*******************************hadron trggered mixing
1780   Fillcorrelation(tracksUNID,bgTracks,cent_v0,zvtx,bSign,kTRUE,kTRUE,kTRUE,"trigassoUNID");//mixcase=kTRUE
1781  if(ffilltrigIDassoUNID && tracksID && tracksID->GetEntriesFast()>0)//***********************************ID trggered mixing
1782   Fillcorrelation(tracksID,bgTracks,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoUNID");//mixcase=kTRUE
1783    }// pool event loop ends mixing case
1784
1785 } //if pool->IsReady() condition ends mixing case
1786  if(tracksUNID) {
1787 if(pool)
1788   pool->UpdatePool(CloneAndReduceTrackList(tracksUNID));
1789  }
1790  }//mixing with unidentified particles
1791
1792
1793  if(ffilltrigUNIDassoID || ffilltrigIDassoID){//mixing with identified particles
1794 AliEventPool* pool1 = fPoolMgr->GetEventPool(cent_v0, zvtx+100);//In the pool1 there is tracksID(i.e associateds are identified)
1795 if (pool1 && pool1->IsReady())
1796   {//start mixing only when pool->IsReady
1797 for (Int_t jMix=0; jMix<pool1->GetCurrentNEvents(); jMix++)
1798   { //pool event loop start
1799  TObjArray* bgTracks2 = pool1->GetEvent(jMix);
1800   if(!bgTracks2) continue;
1801 if(ffilltrigUNIDassoID && tracksUNID && tracksUNID->GetEntriesFast()>0)
1802   Fillcorrelation(tracksUNID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigUNIDassoID");//mixcase=kTRUE
1803 if(ffilltrigIDassoID && tracksID && tracksID->GetEntriesFast()>0)
1804   Fillcorrelation(tracksID,bgTracks2,cent_v0,zvtx,bSign,kFALSE,kTRUE,kTRUE,"trigIDassoID");//mixcase=kTRUE
1805
1806    }// pool event loop ends mixing case
1807 } //if pool1->IsReady() condition ends mixing case
1808
1809 if(tracksID) {
1810 if(pool1)
1811   pool1->UpdatePool(CloneAndReduceTrackList(tracksID));//ownership of tracksasso is with pool now, don't delete it(tracksUNID is with pool)
1812  }
1813  }//mixing with identified particles
1814
1815
1816   //no. of events analyzed
1817 fEventCounter->Fill(7);
1818
1819 //still in main event loop
1820
1821
1822 if(tracksUNID)  delete tracksUNID;
1823
1824 if(tracksID) delete tracksID;
1825
1826
1827 PostData(1, fOutput);
1828
1829 } // *************************event loop ends******************************************//_______________________________________________________________________
1830 TObjArray* AliTwoParticlePIDCorr::CloneAndReduceTrackList(TObjArray* tracks)
1831 {
1832   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1833
1834   TObjArray* tracksClone = new TObjArray;
1835   tracksClone->SetOwner(kTRUE);
1836
1837   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1838   {
1839     LRCParticlePID* particle = (LRCParticlePID*) tracks->UncheckedAt(i);
1840     LRCParticlePID* copy100 = new LRCParticlePID(particle->getparticle(),particle->Charge(), particle->Pt(),particle->Eta(), particle->Phi(), particle->geteffcorrectionval());
1841     copy100->SetUniqueID(particle->GetUniqueID());
1842     tracksClone->Add(copy100);
1843   }
1844
1845   return tracksClone;
1846 }
1847
1848 //--------------------------------------------------------------------------------
1849 void AliTwoParticlePIDCorr::Fillcorrelation(TObjArray *trackstrig,TObjArray *tracksasso,Double_t cent,Float_t vtx,Float_t bSign,Bool_t fPtOrder,Bool_t twoTrackEfficiencyCut,Bool_t mixcase,TString fillup)
1850 {
1851
1852   //before calling this function check that either trackstrig & tracksasso are available
1853
1854  // Eta() is extremely time consuming, therefore cache it for the inner loop here:
1855   TObjArray* input = (tracksasso) ? tracksasso : trackstrig;
1856   TArrayF eta(input->GetEntriesFast());
1857   for (Int_t i=0; i<input->GetEntriesFast(); i++)
1858     eta[i] = ((LRCParticlePID*) input->UncheckedAt(i))->Eta();
1859
1860   //if(trackstrig)
1861     Int_t jmax=trackstrig->GetEntriesFast();
1862   if(tracksasso)
1863      jmax=tracksasso->GetEntriesFast();
1864
1865 // identify K, Lambda candidates and flag those particles
1866     // a TObject bit is used for this
1867 const UInt_t kResonanceDaughterFlag = 1 << 14;
1868     if (fRejectResonanceDaughters > 0)
1869     {
1870       Double_t resonanceMass = -1;
1871       Double_t massDaughter1 = -1;
1872       Double_t massDaughter2 = -1;
1873       const Double_t interval = 0.02;
1874  switch (fRejectResonanceDaughters)
1875       {
1876         case 1: resonanceMass = 0.9; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
1877         case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
1878         case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
1879         default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
1880       }
1881
1882 for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
1883         trackstrig->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
1884  for (Int_t i=0; tracksasso->GetEntriesFast(); i++)
1885           tracksasso->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
1886
1887  for (Int_t i=0; i<trackstrig->GetEntriesFast(); i++)
1888       {
1889       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
1890 for (Int_t j=0; tracksasso->GetEntriesFast(); j++)
1891         {
1892         LRCParticlePID *asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
1893
1894  // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
1895 if (trig->IsEqual(asso)) continue;
1896
1897 if (trig->Charge() * asso->Charge() > 0) continue;
1898
1899  Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
1900
1901 if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
1902           {
1903             mass = GetInvMassSquared(trig->Pt(), trig->Eta(), trig->Phi(), asso->Pt(), asso->Eta(), asso->Phi(), massDaughter1, massDaughter2);
1904
1905             if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
1906             {
1907               trig->SetBit(kResonanceDaughterFlag);
1908               asso->SetBit(kResonanceDaughterFlag);
1909
1910 //            Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
1911             }
1912           }
1913         }
1914       }
1915     }
1916
1917     //two particle correlation filling
1918
1919 for(Int_t i=0;i<trackstrig->GetEntriesFast();i++)
1920     {  //trigger loop starts
1921       LRCParticlePID *trig=(LRCParticlePID*)(trackstrig->UncheckedAt(i));
1922       if(!trig) continue;
1923       Float_t trigpt=trig->Pt();
1924     //to avoid overflow qnd underflow
1925       if(trigpt<fminPtTrig || trigpt>fmaxPtTrig) continue;
1926       Int_t particlepidtrig=trig->getparticle();
1927       if(fTriggerSpeciesSelection){ if (particlepidtrig!=fTriggerSpecies) continue;}
1928
1929       Float_t trigeta=trig->Eta();
1930
1931       // some optimization
1932  if (fTriggerRestrictEta > 0 && TMath::Abs(trigeta) > fTriggerRestrictEta)
1933         continue;
1934
1935 if (fOnlyOneEtaSide != 0)
1936       {
1937         if (fOnlyOneEtaSide * trigeta < 0)
1938           continue;
1939       }
1940   if (fTriggerSelectCharge != 0)
1941         if (trig->Charge() * fTriggerSelectCharge < 0)
1942           continue;
1943
1944       if (fRejectResonanceDaughters > 0)
1945         if (trig->TestBit(kResonanceDaughterFlag)) continue;
1946
1947       Float_t trigphi=trig->Phi();
1948       Float_t trackefftrig=1.0;
1949       if(fapplyTrigefficiency) trackefftrig=trig->geteffcorrectionval();
1950       //cout<<"*******************trackefftrig="<<trackefftrig<<endl;
1951         Double_t* trigval;
1952         Int_t dim=3;
1953         if(fcontainPIDtrig) dim=4;
1954         trigval= new Double_t[dim];
1955       trigval[0] = cent;
1956       trigval[1] = vtx;
1957       trigval[2] = trigpt;
1958 if(fcontainPIDtrig)      trigval[3] = particlepidtrig;
1959
1960       //filling only for same event case(mixcase=kFALSE)
1961
1962       //trigger particle counting only when mixcase=kFALSE i.e for same event correlation function calculation
1963 if(mixcase==kFALSE)
1964   {
1965     if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kTRUE){
1966       if(fillup=="trigassoUNID" )  fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1967     }
1968  if(ffilltrigassoUNID==kTRUE && ffilltrigUNIDassoID==kFALSE){
1969       if(fillup=="trigassoUNID" )  fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1970     }
1971 if(ffilltrigassoUNID==kFALSE && ffilltrigUNIDassoID==kTRUE){
1972       if(fillup=="trigUNIDassoID")  fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1973     }
1974  //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
1975 if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kTRUE){
1976       if(fillup=="trigIDassoID")  fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1977     }
1978  if(ffilltrigIDassoUNID==kTRUE && ffilltrigIDassoID==kFALSE){
1979       if(fillup=="trigIDassoUNID" )  fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1980     }
1981 if(ffilltrigIDassoUNID==kFALSE && ffilltrigIDassoID==kTRUE){
1982       if(fillup=="trigIDassoID")  fTHnTrigcount->Fill(trigval,1.0/trackefftrig);
1983     }
1984
1985  if(fillup=="trigIDassoIDMCTRUTH") fTHnTrigcountMCTruthPrim->Fill(trigval,1.0/trackefftrig); //In truth MC case "Unidentified" means any particle other than pion,kaon or proton and no efficiency correction(default value 1.0)************************be careful!!!!
1986   }
1987
1988     //asso loop starts within trigger loop
1989    for(Int_t j=0;j<jmax;j++)
1990              {
1991     LRCParticlePID *asso=0;
1992     if(!tracksasso)
1993     asso=(LRCParticlePID*)(trackstrig->UncheckedAt(j));
1994     else
1995     asso=(LRCParticlePID*)(tracksasso->UncheckedAt(j));
1996
1997     if(!asso) continue;
1998
1999     //to avoid overflow qnd underflow
2000  if(asso->Pt()<fminPtAsso || asso->Pt()>fmaxPtAsso) continue;//***********************Important
2001
2002  if(fmaxPtAsso==fminPtTrig) {if(asso->Pt()==fminPtTrig) continue;}//******************Think about it!
2003
2004   if(!tracksasso && i==j) continue;
2005
2006    // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
2007    // if (tracksasso && trig->IsEqual(asso))  continue;
2008
2009   if (tracksasso && (trig->GetUniqueID()==asso->GetUniqueID())) continue;
2010
2011  if (fPtOrder)
2012  if (asso->Pt() >= trig->Pt()) continue;
2013
2014   Int_t particlepidasso=asso->getparticle();
2015   if(fAssociatedSpeciesSelection){ if (particlepidasso!=fAssociatedSpecies) continue;}
2016
2017
2018 if (fAssociatedSelectCharge != 0)
2019 if (asso->Charge() * fAssociatedSelectCharge < 0) continue;
2020
2021  if (fSelectCharge > 0)
2022         {
2023           // skip like sign
2024           if (fSelectCharge == 1 && asso->Charge() * trig->Charge() > 0)
2025             continue;
2026
2027           // skip unlike sign
2028           if (fSelectCharge == 2 && asso->Charge() * trig->Charge() < 0)
2029             continue;
2030         }
2031
2032 if (fEtaOrdering)
2033         {
2034           if (trigeta < 0 && asso->Eta() < trigeta)
2035             continue;
2036           if (trigeta > 0 && asso->Eta() > trigeta)
2037             continue;
2038         }
2039
2040 if (fRejectResonanceDaughters > 0)
2041           if (asso->TestBit(kResonanceDaughterFlag))
2042           {
2043 //          Printf("Skipped j=%d", j);
2044             continue;
2045           }
2046
2047         // conversions
2048         if (fCutConversions && asso->Charge() * trig->Charge() < 0)
2049         {
2050           Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2051
2052           if (mass < 0.1)
2053           {
2054             mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.510e-3, 0.510e-3);
2055
2056             //fControlConvResoncances->Fill(0.0, mass);
2057
2058             if (mass < 0.04*0.04)
2059               continue;
2060           }
2061         }
2062
2063         // K0s
2064         if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2065         {
2066           Float_t mass = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.1396);
2067
2068           const Float_t kK0smass = 0.4976;
2069
2070           if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
2071           {
2072             mass = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.1396);
2073
2074             //fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
2075
2076             if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
2077               continue;
2078           }
2079         }
2080
2081         // Lambda
2082         if (fCutResonances && asso->Charge() * trig->Charge() < 0)
2083         {
2084           Float_t mass1 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(), eta[j], asso->Phi(), 0.1396, 0.9383);
2085           Float_t mass2 = GetInvMassSquaredCheap(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2086
2087           const Float_t kLambdaMass = 1.115;
2088
2089           if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
2090           {
2091             mass1 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j], asso->Phi(), 0.1396, 0.9383);
2092
2093             //fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
2094
2095             if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2096               continue;
2097           }
2098           if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
2099           {
2100             mass2 = GetInvMassSquared(trig->Pt(), trigeta, trig->Phi(), asso->Pt(),eta[j] , asso->Phi(), 0.9383, 0.1396);
2101
2102             //fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
2103
2104             if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
2105               continue;
2106           }
2107         }
2108
2109         if (twoTrackEfficiencyCut)
2110         {
2111           // the variables & cuthave been developed by the HBT group
2112           // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
2113           Float_t phi1 = trig->Phi();
2114           Float_t pt1 = trig->Pt();
2115           Float_t charge1 = trig->Charge();
2116           Float_t phi2 = asso->Phi();
2117           Float_t pt2 = asso->Pt();
2118           Float_t charge2 = asso->Charge();
2119
2120           Float_t deta= trigeta - eta[j];
2121
2122  // optimization
2123           if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
2124           {
2125
2126   // check first boundaries to see if is worth to loop and find the minimum
2127             Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
2128             Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
2129
2130  const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
2131
2132             Float_t dphistarminabs = 1e5;
2133             Float_t dphistarmin = 1e5;
2134
2135  if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
2136             {
2137               for (Double_t rad=0.8; rad<2.51; rad+=0.01)
2138               {
2139                 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
2140
2141                 Float_t dphistarabs = TMath::Abs(dphistar);
2142
2143         if (dphistarabs < dphistarminabs)
2144                 {
2145                   dphistarmin = dphistar;
2146                   dphistarminabs = dphistarabs;
2147                 }
2148               }
2149
2150 if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
2151               {
2152 //              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);
2153                 continue;
2154               }
2155 //fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
2156
2157             }
2158           }
2159         }
2160
2161         Float_t trackeffasso=1.0;
2162         if(fapplyAssoefficiency) trackeffasso=asso->geteffcorrectionval();
2163         //cout<<"*******************trackeffasso="<<trackeffasso<<endl;
2164         Float_t deleta=trigeta-eta[j];
2165         Float_t delphi=PhiRange(trigphi-asso->Phi());
2166
2167  //here get the two particle efficiency correction factor
2168         Float_t effweight=trackefftrig*trackeffasso;
2169         //cout<<"*******************effweight="<<effweight<<endl;
2170         Double_t* vars;
2171         vars= new Double_t[kTrackVariablesPair];
2172         vars[0]=cent;
2173         vars[1]=vtx;
2174         vars[2]=trigpt;
2175         vars[3]=asso->Pt();
2176         vars[4]=deleta;
2177         vars[5]=delphi;
2178 if(fcontainPIDtrig && !fcontainPIDasso) vars[6]=particlepidtrig;
2179 if(!fcontainPIDtrig && fcontainPIDasso) vars[6]=particlepidasso;
2180  if(fcontainPIDtrig && fcontainPIDasso){
2181         vars[6]=particlepidtrig;
2182         vars[7]=particlepidasso;
2183  }
2184
2185         //Fill Correlation ThnSparses
2186     if(fillup=="trigassoUNID")
2187       {
2188     if(mixcase==kFALSE)  fTHnCorrUNID->Fill(vars,1.0/effweight);
2189     if(mixcase==kTRUE)   fTHnCorrUNIDmix->Fill(vars,1.0/effweight);
2190       }
2191     if(fillup=="trigIDassoID")
2192       {
2193     if(mixcase==kFALSE)  fTHnCorrID->Fill(vars,1.0/effweight);
2194     if(mixcase==kTRUE)  fTHnCorrIDmix->Fill(vars,1.0/effweight);
2195       }
2196     if(fillup=="trigIDassoIDMCTRUTH")//******************************************for TRUTH MC particles
2197       {//in this case effweight should be 1 always
2198      if(mixcase==kFALSE)  fCorrelatonTruthPrimary->Fill(vars,1.0/effweight);
2199      if(mixcase==kTRUE) fCorrelatonTruthPrimarymix->Fill(vars,1.0/effweight);
2200     }
2201     if(fillup=="trigIDassoUNID" || fillup=="trigUNIDassoID")//****************************be careful
2202       {
2203     if(mixcase==kFALSE)  fTHnCorrIDUNID->Fill(vars,1.0/effweight);
2204     if(mixcase==kTRUE)   fTHnCorrIDUNIDmix->Fill(vars,1.0/effweight);
2205        }
2206
2207 delete[] vars;
2208    }//asso loop ends
2209 delete[] trigval;
2210  }//trigger loop ends
2211
2212 }
2213
2214 //--------------------------------------------------------------------------------
2215 Float_t AliTwoParticlePIDCorr::GetTrackbyTrackeffvalue(AliAODTrack* track,Double_t cent,Float_t evzvtx, Int_t parpid)
2216 {
2217   //This function is called only when applyefficiency=kTRUE; also ensure that "track" is present before calling that function
2218  Int_t effVars[4];
2219  Float_t effvalue=1.;
2220
2221   if(parpid==unidentified)
2222             {
2223             effVars[0] = effcorection[5]->GetAxis(0)->FindBin(cent);
2224             effVars[1] = effcorection[5]->GetAxis(1)->FindBin(evzvtx);
2225             effVars[2] = effcorection[5]->GetAxis(2)->FindBin(track->Pt());
2226             effVars[3] = effcorection[5]->GetAxis(3)->FindBin(track->Eta());
2227             effvalue=effcorection[5]->GetBinContent(effVars);
2228             }
2229 if(parpid==SpPion || parpid==SpKaon)
2230             {
2231               if(fmesoneffrequired && !fkaonprotoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2232                 {
2233             effVars[0] = effcorection[3]->GetAxis(0)->FindBin(cent);
2234             effVars[1] = effcorection[3]->GetAxis(1)->FindBin(evzvtx);
2235             effVars[2] = effcorection[3]->GetAxis(2)->FindBin(track->Pt());
2236             effVars[3] = effcorection[3]->GetAxis(3)->FindBin(track->Eta());
2237             effvalue=effcorection[3]->GetBinContent(effVars);
2238                 }
2239               else{
2240  if(parpid==SpPion)
2241             {
2242             effVars[0] = effcorection[0]->GetAxis(0)->FindBin(cent);
2243             effVars[1] = effcorection[0]->GetAxis(1)->FindBin(evzvtx);
2244             effVars[2] = effcorection[0]->GetAxis(2)->FindBin(track->Pt());
2245             effVars[3] = effcorection[0]->GetAxis(3)->FindBin(track->Eta());
2246             effvalue=effcorection[0]->GetBinContent(effVars);
2247             }
2248
2249  if(parpid==SpKaon)
2250             {
2251             effVars[0] = effcorection[1]->GetAxis(0)->FindBin(cent);
2252             effVars[1] =  effcorection[1]->GetAxis(1)->FindBin(evzvtx);
2253             effVars[2] =  effcorection[1]->GetAxis(2)->FindBin(track->Pt());
2254             effVars[3] =  effcorection[1]->GetAxis(3)->FindBin(track->Eta());
2255             effvalue=effcorection[1]->GetBinContent(effVars);
2256             }
2257               }
2258             }
2259
2260  if(parpid==SpProton)
2261             {
2262             effVars[0] =  effcorection[2]->GetAxis(0)->FindBin(cent);
2263             effVars[1] = effcorection[2]->GetAxis(1)->FindBin(evzvtx);
2264             effVars[2] = effcorection[2]->GetAxis(2)->FindBin(track->Pt());
2265             effVars[3] = effcorection[2]->GetAxis(3)->FindBin(track->Eta());
2266             effvalue=effcorection[2]->GetBinContent(effVars);
2267             }
2268
2269  if(fkaonprotoneffrequired && !fmesoneffrequired && track->Pt()>=fminPtComboeff && track->Pt()<=fmaxPtComboeff)
2270                 {
2271   if(parpid==SpProton || parpid==SpKaon)
2272             {
2273             effVars[0] = effcorection[4]->GetAxis(0)->FindBin(cent);
2274             effVars[1] = effcorection[4]->GetAxis(1)->FindBin(evzvtx);
2275             effVars[2] = effcorection[4]->GetAxis(2)->FindBin(track->Pt());
2276             effVars[3] = effcorection[4]->GetAxis(3)->FindBin(track->Eta());
2277             effvalue=effcorection[4]->GetBinContent(effVars);
2278             }
2279                 }
2280             //    Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
2281      if(effvalue==0.) effvalue=1.;
2282
2283      return effvalue;
2284
2285 }
2286 //-----------------------------------------------------------------------
2287
2288 Int_t AliTwoParticlePIDCorr::ClassifyTrack(AliAODTrack* track,AliAODVertex* vertex,Float_t magfield)
2289 {
2290
2291   if(!track) return 0;
2292   Bool_t trackOK = track->TestFilterBit(fFilterBit);
2293   if(!trackOK) return 0;
2294   //select only primary traks(for data & reco MC tracks)
2295   if(fonlyprimarydatareco && track->GetType()!=AliAODTrack::kPrimary) return 0;
2296   if(track->Charge()==0) return 0;
2297   fHistQA[12]->Fill(track->GetTPCNcls());
2298   Float_t dxy, dz;
2299   dxy = track->DCA();
2300   dz = track->ZAtDCA();
2301   fHistQA[6]->Fill(dxy);
2302   fHistQA[7]->Fill(dz);
2303   Float_t chi2ndf = track->Chi2perNDF();
2304   fHistQA[13]->Fill(chi2ndf);
2305   Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
2306   fHistQA[14]->Fill(nCrossedRowsTPC);
2307   //Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
2308   if (track->GetTPCNclsF()>0) {
2309    Float_t  ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
2310    fHistQA[15]->Fill(ratioCrossedRowsOverFindableClustersTPC);
2311     }
2312 //accepted tracks
2313      Float_t pt=track->Pt();
2314      if(pt< fminPt || pt> fmaxPt) return 0;
2315      if(TMath::Abs(track->Eta())> fmaxeta) return 0;
2316      if(track->Phi()<0. || track->Phi()>2*TMath::Pi()) return 0;
2317      //if (!HasTPCPID(track)) return 0;//trigger & associated particles must have TPC PID,Is it required
2318 // DCA XY
2319         if (fdcacut && fDCAXYCut)
2320         {
2321           if (!vertex)
2322             return 0;
2323
2324           Double_t pos[2];
2325           Double_t covar[3];
2326           AliAODTrack* clone =(AliAODTrack*) track->Clone();
2327           Bool_t success = clone->PropagateToDCA(vertex, magfield, fdcacutvalue, pos, covar);
2328           delete clone;
2329           if (!success)
2330             return 0;
2331
2332 //        Printf("%f", ((AliAODTrack*)part)->DCA());
2333 //        Printf("%f", pos[0]);
2334           if (TMath::Abs(pos[0]) > fDCAXYCut->Eval(track->Pt()))
2335             return 0;
2336         }
2337
2338      fHistQA[8]->Fill(pt);
2339      fHistQA[9]->Fill(track->Eta());
2340      fHistQA[10]->Fill(track->Phi());
2341      return 1;
2342   }
2343   //________________________________________________________________________________
2344 void AliTwoParticlePIDCorr::CalculateNSigmas(AliAODTrack *track)
2345 {
2346 //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 )
2347 Float_t pt=track->Pt();
2348
2349 //it is assumed that every track that passed the filterbit have proper TPC response(!!)
2350 Float_t nsigmaTPCkPion =fPID->NumberOfSigmasTPC(track, AliPID::kPion);
2351 Float_t nsigmaTPCkKaon =fPID->NumberOfSigmasTPC(track, AliPID::kKaon);
2352 Float_t nsigmaTPCkProton =fPID->NumberOfSigmasTPC(track, AliPID::kProton);
2353
2354 Float_t nsigmaTOFkProton=999.,nsigmaTOFkKaon=999.,nsigmaTOFkPion=999.;
2355 Float_t nsigmaTPCTOFkProton=999.,nsigmaTPCTOFkKaon=999.,nsigmaTPCTOFkPion=999.;
2356
2357  if(HasTOFPID(track) && pt>fPtTOFPIDmin)
2358    {
2359
2360 nsigmaTOFkPion =fPID->NumberOfSigmasTOF(track, AliPID::kPion);
2361 nsigmaTOFkKaon =fPID->NumberOfSigmasTOF(track, AliPID::kKaon);
2362 nsigmaTOFkProton =fPID->NumberOfSigmasTOF(track, AliPID::kProton);
2363 //---combined
2364 nsigmaTPCTOFkPion   = TMath::Sqrt(nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion);
2365 nsigmaTPCTOFkKaon   = TMath::Sqrt(nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon);
2366 nsigmaTPCTOFkProton = TMath::Sqrt(nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton);
2367
2368
2369    }
2370 else{
2371     // --- combined
2372     // if TOF is missing and below fPtTOFPID only the TPC information is used
2373     nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
2374     nsigmaTPCTOFkKaon   = TMath::Abs(nsigmaTPCkKaon);
2375     nsigmaTPCTOFkPion   = TMath::Abs(nsigmaTPCkPion);
2376
2377   }
2378
2379 //set data member fnsigmas
2380   fnsigmas[SpPion][NSigmaTPC]=nsigmaTPCkPion;
2381   fnsigmas[SpKaon][NSigmaTPC]=nsigmaTPCkKaon;
2382   fnsigmas[SpProton][NSigmaTPC]=nsigmaTPCkProton;
2383
2384   //for all tracks below fPtTOFPIDmin  and also for tracks above fPtTOFPIDmin without proper TOF response these TOF nsigma values will be 999.
2385   fnsigmas[SpPion][NSigmaTOF]=nsigmaTOFkPion;
2386   fnsigmas[SpKaon][NSigmaTOF]=nsigmaTOFkKaon;
2387   fnsigmas[SpProton][NSigmaTOF]=nsigmaTOFkProton;
2388
2389  //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)
2390   fnsigmas[SpPion][NSigmaTPCTOF]=nsigmaTPCTOFkPion;
2391   fnsigmas[SpKaon][NSigmaTPCTOF]=nsigmaTPCTOFkKaon;
2392   fnsigmas[SpProton][NSigmaTPCTOF]=nsigmaTPCTOFkProton;
2393
2394
2395 }
2396 //----------------------------------------------------------------------------
2397 Int_t AliTwoParticlePIDCorr::FindMinNSigma(AliAODTrack *track)
2398 {
2399   //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2400 if(fRequestTOFPID && track->Pt()>fPtTOFPIDmin && track->Pt()<=fPtTOFPIDmax && (!HasTOFPID(track)) )return SpUndefined;//so any track having Pt>0.6 && Pt<=4.0 Gev withot having proper TOF response will be defined as SpUndefined
2401 //get the identity of the particle with the minimum Nsigma
2402   Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
2403   switch (fPIDType){
2404   case NSigmaTPC:
2405     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2406     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPC])  ;
2407     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPC])  ;
2408     break;
2409   case NSigmaTOF:
2410     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2411     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTOF])  ;
2412     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTOF])  ;
2413     break;
2414   case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
2415     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
2416     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
2417     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
2418     break;
2419   }
2420
2421   if(track->Pt()<=fPtTOFPIDmax)  {
2422  // guess the particle based on the smaller nsigma (within fNSigmaPID)
2423   if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton )) return SpUndefined;//it is the default value for the three
2424 if( ( nsigmaKaon   < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID)) return SpKaon;
2425 if( ( nsigmaPion   < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID)) return SpPion;
2426 if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID)) return SpProton;
2427
2428 // else, return undefined
2429   return SpUndefined;
2430   }
2431   else  {//asymmetric nsigma cut around pion's mean position for tracks having Pt>4.0 Gev
2432     if(fminprotonsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxprotonsigmacut) return SpProton;
2433     if(fminpionsigmacut<fnsigmas[SpPion][NSigmaTPC] && fnsigmas[SpPion][NSigmaTPC]<fmaxpionsigmacut) return SpPion;
2434 // else, return undefined(here we don't detect kaons)
2435   return SpUndefined;
2436   }
2437
2438 }
2439
2440 //------------------------------------------------------------------------------------------
2441 Bool_t* AliTwoParticlePIDCorr::GetDoubleCounting(AliAODTrack * trk){
2442   //this function is always called after calling the function CalculateNSigmas(AliAODTrack *track)
2443
2444   //if a particle has double counting set fHasDoubleCounting[ipart]=kTRUE
2445   //fill DC histos
2446   for(Int_t ipart=0;ipart<NSpecies;ipart++)fHasDoubleCounting[ipart]=kFALSE;//array with kTRUE for second (or third) identity of the track
2447
2448   Int_t MinNSigma=FindMinNSigma(trk);//not filling the NSigmaRec histos
2449
2450
2451   if(MinNSigma==SpUndefined)return fHasDoubleCounting;//in case of undefined no Double counting
2452
2453   Float_t nsigmaPion=999., nsigmaKaon=999., nsigmaProton=999.;
2454   switch (fPIDType) {
2455   case NSigmaTPC:
2456     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPC]);
2457     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPC])  ;
2458     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPC])  ;
2459     break;
2460   case NSigmaTOF:
2461     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTOF]);
2462     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTOF])  ;
2463     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTOF])  ;
2464     break;
2465   case NSigmaTPCTOF://In case of no TOF matching the combined nsigma is the TPC one
2466     nsigmaProton  =  TMath::Abs(fnsigmas[SpProton][NSigmaTPCTOF]);
2467     nsigmaKaon    =  TMath::Abs(fnsigmas[SpKaon][NSigmaTPCTOF])  ;
2468     nsigmaPion    =  TMath::Abs(fnsigmas[SpPion][NSigmaTPCTOF])  ;
2469     break;
2470   }
2471
2472   //there is chance of overlapping only for particles having pt below 4.0 GEv
2473   if(trk->Pt()<=4.0){
2474   if(nsigmaPion<fNSigmaPID && MinNSigma!=SpPion)fHasDoubleCounting[SpPion]=kTRUE;
2475   if(nsigmaKaon<fNSigmaPID && MinNSigma!=SpKaon)fHasDoubleCounting[SpKaon]=kTRUE;
2476   if(nsigmaProton<fNSigmaPID && MinNSigma!=SpProton)fHasDoubleCounting[SpProton]=kTRUE;
2477
2478   }
2479
2480
2481   return fHasDoubleCounting;
2482 }
2483
2484 //////////////////////////////////////////////////////////////////////////////////////////////////
2485 Int_t AliTwoParticlePIDCorr::GetParticle(AliAODTrack * trk){
2486   //return the specie according to the minimum nsigma value
2487   //no double counting, this has to be evaluated using CheckDoubleCounting()
2488   //Printf("fPtTOFPID %.1f, fRequestTOFPID %d, fNSigmaPID %.1f, fPIDType %d",fPtTOFPID,fRequestTOFPID,fNSigmaPID,fPIDType);
2489
2490   CalculateNSigmas(trk);//fill the data member fnsigmas with the nsigmas value [ipart][iPID]
2491
2492   if(fUseExclusiveNSigma){
2493     Bool_t *HasDC;
2494     HasDC=GetDoubleCounting(trk);
2495     for(Int_t ipart=0;ipart<NSpecies;ipart++){
2496       if(HasDC[ipart]==kTRUE)  return SpUndefined;
2497     }
2498     return FindMinNSigma(trk);//NSigmaRec distr filled here
2499   }
2500   else return FindMinNSigma(trk);//NSigmaRec distr filled here
2501
2502 }
2503
2504 //-------------------------------------------------------------------------------------
2505 Bool_t
2506 AliTwoParticlePIDCorr::HasTPCPID(AliAODTrack *track) const
2507 {
2508   // check PID signal
2509    AliPIDResponse::EDetPidStatus statustpc = fPID->CheckPIDStatus(AliPIDResponse::kTPC,track);
2510    if(statustpc!=AliPIDResponse::kDetPidOk) return kFALSE;
2511    //ULong_t status=track->GetStatus();
2512    //if  (!( (status & AliAODTrack::kTPCpid  ) == AliAODTrack::kTPCpid  )) return kFALSE;//remove light nuclei
2513    //if (track->GetTPCsignal() <= 0.) return kFALSE;
2514    // 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.
2515
2516   return kTRUE;
2517 }
2518 //___________________________________________________________
2519
2520 Bool_t
2521 AliTwoParticlePIDCorr::HasTOFPID(AliAODTrack *track) const
2522 {
2523   // check TOF matched track
2524   //ULong_t status=track->GetStatus();
2525   //if  (!( (status & AliAODTrack::kITSin  ) == AliAODTrack::kITSin  )) return kFALSE;
2526  AliPIDResponse::EDetPidStatus statustof = fPID->CheckPIDStatus(AliPIDResponse::kTOF,track);
2527  if(statustof!= AliPIDResponse::kDetPidOk) return kFALSE;
2528   if(track->Pt()<=fPtTOFPIDmin) return kFALSE;
2529  //if(!((status & AliAODTrack::kTOFpid  ) == AliAODTrack::kTOFpid  )) return kFALSE;
2530  //Float_t probMis = fPIDresponse->GetTOFMismatchProbability(track);
2531  // if (probMis > 0.01) return kFALSE;
2532 if(fRemoveTracksT0Fill)
2533     {
2534 Int_t startTimeMask = fPID->GetTOFResponse().GetStartTimeMask(track->P());
2535       if (startTimeMask < 0)return kFALSE;
2536     }
2537   return kTRUE;
2538 }
2539
2540 //________________________________________________________________________
2541 Float_t AliTwoParticlePIDCorr :: GetBeta(AliAODTrack *track)
2542 {
2543   //it is called only when TOF PID is available
2544   Double_t p = track->P();
2545   Double_t time=track->GetTOFsignal()-fPID->GetTOFResponse().GetStartTime(p);
2546   Double_t timei[5];
2547   track->GetIntegratedTimes(timei);
2548   return timei[0]/time;
2549 }
2550 //------------------------------------------------------------------------------------------------------
2551
2552 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)
2553 {
2554   // calculate inv mass squared
2555   // same can be achieved, but with more computing time with
2556   /*TLorentzVector photon, p1, p2;
2557   p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
2558   p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
2559   photon = p1+p2;
2560   photon.M()*/
2561
2562   Float_t tantheta1 = 1e10;
2563
2564   if (eta1 < -1e-10 || eta1 > 1e-10)
2565     tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
2566
2567   Float_t tantheta2 = 1e10;
2568   if (eta2 < -1e-10 || eta2 > 1e-10)
2569     tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
2570
2571   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
2572   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
2573
2574   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 ) ) );
2575
2576   return mass2;
2577 }
2578 //---------------------------------------------------------------------------------
2579
2580 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)
2581 {
2582   // calculate inv mass squared approximately
2583
2584   Float_t tantheta1 = 1e10;
2585
2586   if (eta1 < -1e-10 || eta1 > 1e-10)
2587   {
2588     Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
2589     tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
2590   }
2591
2592   Float_t tantheta2 = 1e10;
2593   if (eta2 < -1e-10 || eta2 > 1e-10)
2594   {
2595     Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
2596     tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
2597   }
2598
2599   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
2600   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
2601
2602   // fold onto 0...pi
2603   Float_t deltaPhi = TMath::Abs(phi1 - phi2);
2604   while (deltaPhi > TMath::TwoPi())
2605     deltaPhi -= TMath::TwoPi();
2606   if (deltaPhi > TMath::Pi())
2607     deltaPhi = TMath::TwoPi() - deltaPhi;
2608
2609   Float_t cosDeltaPhi = 0;
2610   if (deltaPhi < TMath::Pi()/3)
2611     cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
2612   else if (deltaPhi < 2*TMath::Pi()/3)
2613     cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
2614   else
2615     cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
2616
2617   Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
2618
2619 //   Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
2620
2621   return mass2;
2622 }
2623 //--------------------------------------------------------------------------------
2624 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)
2625 {
2626   //
2627   // calculates dphistar
2628   //
2629
2630   Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2631
2632   static const Double_t kPi = TMath::Pi();
2633
2634   // circularity
2635 //   if (dphistar > 2 * kPi)
2636 //     dphistar -= 2 * kPi;
2637 //   if (dphistar < -2 * kPi)
2638 //     dphistar += 2 * kPi;
2639
2640   if (dphistar > kPi)
2641     dphistar = kPi * 2 - dphistar;
2642   if (dphistar < -kPi)
2643     dphistar = -kPi * 2 - dphistar;
2644   if (dphistar > kPi) // might look funny but is needed
2645     dphistar = kPi * 2 - dphistar;
2646
2647   return dphistar;
2648 }
2649 //_________________________________________________________________________
2650 void AliTwoParticlePIDCorr ::DefineEventPool()
2651 {
2652 const Int_t MaxNofEvents=1000;
2653 const Int_t MaxNofTracks=50000;
2654 const Int_t NofVrtxBins=10+(1+10)*2;
2655 Double_t ZvrtxBins[NofVrtxBins+1]={ -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10,
2656                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110,
2657                                       190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210
2658                                       };
2659  if (fSampleType=="pp"){
2660 const Int_t  NofCentBins=5;
2661 Double_t CentralityBins[NofCentBins+1]={0.,10., 20., 40.,80.,200.1};
2662 fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
2663  }
2664 if (fSampleType=="pPb" || fSampleType=="PbPb")
2665   {
2666 const Int_t  NofCentBins=15;
2667 Double_t CentralityBins[NofCentBins+1]={0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
2668 fPoolMgr = new AliEventPoolManager(MaxNofEvents,MaxNofTracks,NofCentBins,CentralityBins,NofVrtxBins,ZvrtxBins);
2669   }
2670 fPoolMgr->SetTargetValues(MaxNofTracks, 0.1, 5);
2671
2672 //if(!fPoolMgr) return kFALSE;
2673 //return kTRUE;
2674
2675 }
2676 //------------------------------------------------------------------------
2677 Double_t* AliTwoParticlePIDCorr::GetBinning(const char* configuration, const char* tag, Int_t& nBins)
2678 {
2679   // This method is a copy from AliUEHist::GetBinning
2680   // takes the binning from <configuration> identified by <tag>
2681   // configuration syntax example:
2682   // 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
2683   // phi: .....
2684   //
2685   // returns bin edges which have to be deleted by the caller
2686
2687   TString config(configuration);
2688   TObjArray* lines = config.Tokenize("\n");
2689   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2690   {
2691     TString line(lines->At(i)->GetName());
2692     if (line.BeginsWith(TString(tag) + ":"))
2693     {
2694       line.Remove(0, strlen(tag) + 1);
2695       line.ReplaceAll(" ", "");
2696       TObjArray* binning = line.Tokenize(",");
2697       Double_t* bins = new Double_t[binning->GetEntriesFast()];
2698       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2699         bins[j] = TString(binning->At(j)->GetName()).Atof();
2700
2701       nBins = binning->GetEntriesFast() - 1;
2702
2703       delete binning;
2704       delete lines;
2705       return bins;
2706     }
2707   }
2708
2709   delete lines;
2710   AliFatal(Form("Tag %s not found in %s", tag, configuration));
2711   return 0;
2712 }
2713 //_______________________________________________________________________________
2714 void AliTwoParticlePIDCorr::SetAsymmetricBin(THnSparse *h,Int_t dim,Double_t *arraybin,Int_t arraybinsize,TString axisTitle)
2715 {
2716   TAxis *axis = 0x0;
2717   axis = h->GetAxis(dim);
2718   axis->Set(arraybinsize,arraybin);
2719   axis->SetTitle(axisTitle);
2720 }
2721
2722 //________________________________________________________________________
2723 void AliTwoParticlePIDCorr::Terminate(Option_t *)
2724 {
2725   // Draw result to screen, or perform fitting, normalizations
2726   // Called once at the end of the query
2727   fOutput = dynamic_cast<TList*> (GetOutputData(1));
2728   if(!fOutput) { Printf("ERROR: could not retrieve TList fOutput"); return; }
2729
2730
2731 }
2732 //------------------------------------------------------------------
2733