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