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