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