PWGCF/Correlations/DPhi/TriggerPID/AliTwoParticlePIDCorr.cxx

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