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