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