PWGCF/FEMTOSCOPY/V0LamAnalysis/AliAnalysisV0Lam.cxx

   1 #include <iostream>
   2 #include <math.h>
   3 #include "TChain.h"
   4 #include "TFile.h"
   5 #include "TKey.h"
   6 #include "TObject.h"
   7 #include "TObjString.h"
   8 #include "TList.h"
   9 #include "TTree.h"
  10 #include "TH1F.h"
  11 #include "TH1D.h"
  12 #include "TH2D.h"
  13 #include "TH3D.h"
  14 #include "TCanvas.h"
  15 #include "TRandom3.h"
  16 #include "AliAnalysisTask.h"
  17 #include "AliAnalysisManager.h"
  18 #include "AliESDEvent.h"
  19 #include "AliESDInputHandler.h"
  20 #include "AliESDtrackCuts.h"
  21 #include "AliAODEvent.h"
  22 #include "AliAODInputHandler.h"
  23 #include "AliAODMCParticle.h"
  24 #include "AliAODv0.h"
  25 #include "AliAODRecoDecay.h"
  26 #include "AliAnalysisV0Lam.h"
  27 #include "AliAODMCHeader.h"
  28 #include "AliGenHijingEventHeader.h"
  29 #define PI 3.1415927
  30
  31 // Analysis task for studying lambda-lambda femtoscopic correlations
  32 // Author: Jai Salzwedel, jai.salzwedel@cern.ch
  33
  34
  35 ClassImp(AliAnalysisV0Lam)
  36
  37 //________________________________________________________________________
  38 AliAnalysisV0Lam::AliAnalysisV0Lam():
  39 AliAnalysisTaskSE(),
  40   fAOD(0),
  41   fOutputList(0),
  42   fpidAOD(0)
  43 {
  44 }
  45 //________________________________________________________________________
  46 AliAnalysisV0Lam::AliAnalysisV0Lam(const char *name)
  47 : AliAnalysisTaskSE(name),
  48   fAOD(0),
  49   fOutputList(0),
  50   fpidAOD(0)
  51 {
  52   // Define output slots here
  53   // Output slot #1
  54   DefineOutput(1, TList::Class());
  55 }
  56 //________________________________________________________________________
  57 AliAnalysisV0Lam::~AliAnalysisV0Lam()
  58 {
  59   // Destructor
  60
  61   for(unsigned short i=0; i<zVertexBins; i++)
  62   {
  63     for(unsigned short j=0; j<nCentBins; j++)
  64     {
  65       delete fEC[i][j];
  66     }
  67     delete[] fEC[i];
  68   }
  69   delete[] fEC;
  70   delete fCutProcessor;
  71   if(fOutputList) delete fOutputList; //This cleans up all output hists
  72 }
  73 //________________________________________________________________________
  74 void AliAnalysisV0Lam::MyInit()
  75 {
  76   // Set global variables
  77   fEventCount=           0;
  78   fPDGLambda =    1.115683;
  79   fPDGK0     =     .497614;
  80   fPDGProton =     .938272;
  81   fPDGPion   =    .1395702;
  82   fEtaDaughter =       0.8;    //max eta of daughter particles - default 0.8
  83   fMassWindowK0 =    0.018;    //accept .480-.514
  84   fMassWindowLam = 0.00568;    //accept 1.11003-1.121363
  85   fTOFLow =            0.8;    // Lower |P| limit
  86   fSigmaCutTPCProton = 3.0;    // max Nsigma allowed
  87   fSigmaCutTPCPion   = 3.0;
  88   fSigmaCutTOFProton = 4.0;
  89   fSigmaCutTOFPion   = 4.0;
  90   fIsUsingVariableAvgSepCut = kFALSE; //Relevant for two-track cuts in DoPairStudies().
  91   fMaxV0Mult = 700;
  92   int numberVariableAvgSepCuts = 12;
  93
  94
  95   // Setup V0 cut processor
  96   fCutProcessor = new AliAnalysisV0LamCutProcessing(fOutputList);
  97   fNumberOfVariableCutValues = fCutProcessor->GetNumberOfVariableCutValues();
  98   fDefaultVariableCutIndex = fCutProcessor->GetVariableCutIndex();
  99   if(fIsUsingVariableAvgSepCut){
 100     fNumberOfCfVariableCutValues = numberVariableAvgSepCuts;
 101   }
 102   else fNumberOfCfVariableCutValues = fNumberOfVariableCutValues;
 103   cout<<"Number of variable cf cut values: "<<fNumberOfCfVariableCutValues<<endl;
 104   fTotalLambda =         0; //tabulates number of v0s found for locally run jobs
 105   fTotalAntiLambda =     0;
 106
 107   //setup event collection for event mixing
 108   fEC = new AliAnalysisV0LamEventCollection **[zVertexBins];
 109   for(unsigned short i=0; i<zVertexBins; i++)
 110   {
 111     fEC[i] = new AliAnalysisV0LamEventCollection *[nCentBins];
 112     for(unsigned short j=0; j<nCentBins; j++)
 113     {
 114       fEC[i][j] = new AliAnalysisV0LamEventCollection(nEventsToMix+1, fMaxV0Mult);
 115     }
 116   }
 117   AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 118   fpidAOD = aodH->GetAODpidUtil();
 119 }
 120 //________________________________________________________________________
 121 void AliAnalysisV0Lam::UserCreateOutputObjects()
 122 {
 123   // Create output histograms.
 124   // Histograms are added to fOutputList
 125   // When fOutputList is deleted, it automatically cleans up all
 126   // associated histograms
 127
 128   fOutputList = new TList();
 129   fOutputList->SetOwner();
 130   MyInit();// Initialize my settings
 131
 132   fMultDistRough = new TH1F("fMultDistRough","Multiplicity Distribution",301,-.5,3001-.5);
 133   fMultDistRough->GetXaxis()->SetTitle("Event Multiplicity (pions)");
 134   fMultDistRough->GetYaxis()->SetTitle("# of events");
 135   fOutputList->Add(fMultDistRough);
 136
 137   fCentrality = new TH1F("fCentrality", "Centrality Percentage of Event", 100, 0., 100.);
 138   fOutputList->Add(fCentrality);
 139
 140   //TPC signal vs track momentum
 141   fTPCVsPPosLam = new TH2F("fTPCVsPPosLam","",50,0.,2.,250,0.,500.);
 142   fOutputList->Add(fTPCVsPPosLam);
 143     //TPC signal vs track momentum
 144   fTPCVsPNegLam = new TH2F("fTPCVsPNegLam","",50,0.,2.,250,0.,500.);
 145   fOutputList->Add(fTPCVsPNegLam);
 146     //TPC signal vs track momentum
 147   fTPCVsPPosALam = new TH2F("fTPCVsPPosALam","",50,0.,2.,250,0.,500.);
 148   fOutputList->Add(fTPCVsPPosALam);
 149     //TPC signal vs track momentum
 150   fTPCVsPNegALam = new TH2F("fTPCVsPNegALam","",50,0.,2.,250,0.,500.);
 151   fOutputList->Add(fTPCVsPNegALam);
 152
 153   //V0 Shared daughter culling statistics
 154   fDataCompeted = new TH1F("fDataCompeted","",26, -0.5, 25.5);
 155   fOutputList->Add(fDataCompeted);
 156
 157   fDataCulled = new TH1F("fDataCulled","",26, -0.5, 25.5);
 158   fOutputList->Add(fDataCulled);
 159
 160   fRemainingV0s = new TH1F("fRemainingV0s","",26, -0.5, 25.5);
 161   fOutputList->Add(fRemainingV0s);
 162
 163   fRemainingFrac = new TH1F("fRemainingFrac","",101, -.005, 1.005);
 164   fOutputList->Add(fRemainingFrac);
 165
 166   fRemainingFromBeginningToV0Finder = new TH2F("fRemainingFromBeginningToV0Finder","Fraction Remaining At V0Finder Stage", AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1, 21, -0.025, 1.025);
 167   fOutputList->Add(fRemainingFromBeginningToV0Finder);
 168
 169   fMCTruthOfOriginalParticles = new TH1F("fMCTruthOfOriginalParticles","MC Truth of Original Particles in Event", AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1);
 170   fOutputList->Add(fMCTruthOfOriginalParticles);
 171
 172   fMCTruthOfV0FinderParticles = new TH1F("fMCTruthOfV0FinderParticles","MC Truth of V0Finder Particles in Event", AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1);
 173   fOutputList->Add(fMCTruthOfV0FinderParticles);
 174
 175   SetBinsOnOriginHists(fRemainingFromBeginningToV0Finder);
 176   SetBinsOnOriginHists(fMCTruthOfOriginalParticles);
 177   SetBinsOnOriginHists(fMCTruthOfV0FinderParticles);
 178   //The first dimension is the index of the variable cut value
 179   fRemainingFromBeginningToRecon = new TH3F("fRemainingFromBeginningToRecon", "Fraction Remaining After Reconstruction", fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5,  AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1, 21, -0.025, 1.025);
 180   fRemainingFromV0FinderToRecon = new TH3F("fRemainingFromV0FinderToRecon", "Fraction From V0Finder That Remain After Reconstruction Stage", fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5,  AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1, 21, -0.025, 1.025);
 181   fMCTruthOfReconstructedParticles = new TH2F("fMCTruthOfReconstructedParticles", "MC Truth of Reconstructed Particles in Event", AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1, fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5);
 182
 183   // Particle multiplicities
 184   fMultDistLambda = new TH2F("fMultDistLambda", "Lambda multiplicity", fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5,  21, -0.5, 21-0.5);
 185   fMultDistAntiLambda = new TH2F("fMultDistAntiLambda", "AntiLambda multiplicity", fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5,  21, -0.5, 21-0.5);
 186   fMultCentLambda = new TH3F("fMultCentLambda", "Lambda multiplicity vs centrality", fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5,  nCentBins, .5, nCentBins+.5, 21, -0.5, 21-0.5);
 187   fMultCentAntiLambda =  new TH3F("fMultCentAntiLambda", "AntiLambda multiplicity vs centrality", fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5,  nCentBins, .5, nCentBins+.5, 21, -0.5, 21-0.5);
 188   SetBinsOnOriginHists(fRemainingFromBeginningToRecon);
 189   SetBinsOnOriginHists(fRemainingFromV0FinderToRecon);
 190   SetBinsOnOriginHists(fMCTruthOfReconstructedParticles);
 191   fMultDistLambda->GetXaxis()->SetTitle("Var Bin");
 192   fMultDistLambda->GetYaxis()->SetTitle("Event Multiplicity (Lambdas)");
 193   fMultDistAntiLambda->GetXaxis()->SetTitle("Var Bin");
 194   fMultDistAntiLambda->GetXaxis()->SetTitle("Event Multiplicity (AntiLambdas)");
 195   fMultCentLambda->GetXaxis()->SetTitle("Var Bin");
 196   fMultCentLambda->GetYaxis()->SetTitle("Centrality");
 197   fMultCentLambda->GetZaxis()->SetTitle("Event Multiplicity (Lambdas)");
 198   fMultCentAntiLambda->GetXaxis()->SetTitle("Var Bin");
 199   fMultCentAntiLambda->GetYaxis()->SetTitle("Centrality");
 200   fMultCentAntiLambda->GetZaxis()->SetTitle("Event Multiplicity (AntiLambdas)");
 201   fOutputList->Add(fRemainingFromBeginningToRecon);
 202   fOutputList->Add(fRemainingFromV0FinderToRecon);
 203   fOutputList->Add(fMCTruthOfReconstructedParticles);
 204   fOutputList->Add(fMultDistLambda);
 205   fOutputList->Add(fMultDistAntiLambda);
 206   fOutputList->Add(fMultCentLambda);
 207   fOutputList->Add(fMultCentAntiLambda);
 208
 209
 210   fMCFakeParticleIdentity = new TH1F("fMCFakeParticleIdentity", "Breakdown of fake particles", 3, 0, 3);
 211   fOutputList->Add(fMCFakeParticleIdentity);
 212   fMCFakeParticleIdentity->GetXaxis()->SetBinLabel(1,"Fake Lambda");
 213   fMCFakeParticleIdentity->GetXaxis()->SetBinLabel(2,"Fake AntiLambda");
 214   fMCFakeParticleIdentity->GetXaxis()->SetBinLabel(3,"Fake K0Short");
 215
 216   fMCOtherV0Identity = new TH1F("fMCOtherV0Identity", "Breakdown of otherV0s particles", 3, 0, 3);
 217   fOutputList->Add(fMCOtherV0Identity);
 218   fMCOtherV0Identity->GetXaxis()->SetBinLabel(1,"Fake Lambda");
 219   fMCOtherV0Identity->GetXaxis()->SetBinLabel(2,"Fake AntiLambda");
 220   fMCOtherV0Identity->GetXaxis()->SetBinLabel(3,"Fake K0Short");
 221   int kTBins = 200;
 222   int maxKtBin = 10.;
 223   int kStarBins = 800;
 224   int maxKStar = 2.;
 225   //Pair kT Tracking: Centrality bins, kT bins, k* bins
 226   fKtLamLamSig = new TH3F("fKtLamLamSig", "LamLam Pair Kt Same Event", nCentBins, .5, nCentBins +.5, kTBins, 0., maxKtBin, kStarBins, 0., maxKStar);
 227   fOutputList->Add(fKtLamLamSig);
 228   fKtALamALamSig = new TH3F("fKtALamALamSig", "ALamALam Pair Kt Same Event", nCentBins, .5, nCentBins +.5, kTBins, 0., maxKtBin, kStarBins, 0., maxKStar);
 229   fOutputList->Add(fKtALamALamSig);
 230   fKtLamALamSig = new TH3F("fKtLamALamSig", "LamALam Pair Kt Same Event", nCentBins, .5, nCentBins +.5, kTBins, 0., maxKtBin, kStarBins, 0., maxKStar);
 231   fOutputList->Add(fKtLamALamSig);
 232     fKtLamLamBkg = new TH3F("fKtLamLamBkg", "LamLam Pair Kt Mixed Event", nCentBins, .5, nCentBins +.5, kTBins, 0., maxKtBin, kStarBins, 0., maxKStar);
 233   fOutputList->Add(fKtLamLamBkg);
 234   fKtALamALamBkg = new TH3F("fKtALamALamBkg", "ALamALam Pair Kt Mixed Event", nCentBins, .5, nCentBins +.5, kTBins, 0., maxKtBin, kStarBins, 0., maxKStar);
 235   fOutputList->Add(fKtALamALamBkg);
 236   fKtLamALamBkg = new TH3F("fKtLamALamBkg", "LamALam Pair Kt Mixed Event", nCentBins, .5, nCentBins +.5, kTBins, 0., maxKtBin, kStarBins, 0., maxKStar);
 237   fOutputList->Add(fKtLamALamBkg);
 238
 239   //Momentum resolution (pre-)correction analysis
 240   fHistPsmearingKreconVsKtruthLL = new TH2F ("fHistPsmearingKreconVsKtruthLL","Relative momentum resolution, recon vs truth LL", kStarBins, 0., maxKStar, kStarBins, 0., maxKStar);
 241   fOutputList->Add(fHistPsmearingKreconVsKtruthLL);
 242   fHistPsmearingKreconMinusKtruthLL = new TH1F ("fHistPsmearingKreconMinusKtruthLL","Relative momentum resolution, recon - truth LL", kStarBins, -0.5, 0.5);
 243   fOutputList->Add(fHistPsmearingKreconMinusKtruthLL);
 244
 245     fHistPsmearingKreconVsKtruthAA = new TH2F ("fHistPsmearingKreconVsKtruthAA","Relative momentum resolution, recon vs truth AA", kStarBins, 0., maxKStar, kStarBins, 0., maxKStar);
 246   fOutputList->Add(fHistPsmearingKreconVsKtruthAA);
 247   fHistPsmearingKreconMinusKtruthAA = new TH1F ("fHistPsmearingKreconMinusKtruthAA","Relative momentum resolution, recon - truth AA", kStarBins, -0.5, 0.5);
 248   fOutputList->Add(fHistPsmearingKreconMinusKtruthAA);
 249
 250     fHistPsmearingKreconVsKtruthLA = new TH2F ("fHistPsmearingKreconVsKtruthLA","Relative momentum resolution, recon vs truth LA", kStarBins, 0., maxKStar, kStarBins, 0., maxKStar);
 251   fOutputList->Add(fHistPsmearingKreconVsKtruthLA);
 252   fHistPsmearingKreconMinusKtruthLA = new TH1F ("fHistPsmearingKreconMinusKtruthLA","Relative momentum resolution, recon - truth LA", kStarBins, -0.5, 0.5);
 253   fOutputList->Add(fHistPsmearingKreconMinusKtruthLA);
 254
 255   /////////Signal Distributions///////////////////
 256   //First bin is variable cut value, second bin is centrality, third bin is Kstar
 257   fSignalLamLam = new TH3F("fSignalLamLam","Same Event Pair Distribution", fNumberOfCfVariableCutValues, -0.5, fNumberOfCfVariableCutValues -0.5, nCentBins, .5, nCentBins+.5, kStarBins, 0., maxKStar);
 258   fOutputList->Add(fSignalLamLam);
 259
 260   fBkgLamLam = new TH3F("fBkgLamLam","Mixed Event Pair Distribution", fNumberOfCfVariableCutValues, -0.5, fNumberOfCfVariableCutValues -0.5, nCentBins, .5, nCentBins+.5, kStarBins, 0., maxKStar);
 261   fOutputList->Add(fBkgLamLam);
 262
 263   fSignalALamALam = new TH3F("fSignalALamALam","Same Event Pair Distribution", fNumberOfCfVariableCutValues, -0.5, fNumberOfCfVariableCutValues -0.5, nCentBins, .5, nCentBins+.5, kStarBins, 0., maxKStar);
 264   fOutputList->Add(fSignalALamALam);
 265
 266   fBkgALamALam = new TH3F("fBkgALamALam","Mixed Event Pair Distribution", fNumberOfCfVariableCutValues, -0.5, fNumberOfCfVariableCutValues -0.5, nCentBins, .5, nCentBins+.5, kStarBins, 0., maxKStar);
 267   fOutputList->Add(fBkgALamALam);
 268
 269   fSignalLamALam = new TH3F("fSignalLamALam","Same Event Pair Distribution", fNumberOfCfVariableCutValues, -0.5, fNumberOfCfVariableCutValues -0.5, nCentBins, .5, nCentBins+.5, kStarBins, 0., maxKStar);
 270   fOutputList->Add(fSignalLamALam);
 271
 272   fBkgLamALam = new TH3F("fBkgLamALam","Mixed Event Pair Distribution", fNumberOfCfVariableCutValues, -0.5, fNumberOfCfVariableCutValues -0.5, nCentBins, .5, nCentBins+.5, kStarBins, 0., maxKStar);
 273   fOutputList->Add(fBkgLamALam);
 274
 275   // Daughter pair separation distributions
 276   // Binned according to (average separation, qinv)
 277   int avgSepBins = 200;
 278   double avgSepMaxValue = 20.;
 279   fSignalLamLamProtSep = new TH2F ("fSignalLamLamProtSep","Proton pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 280   fOutputList->Add(fSignalLamLamProtSep);
 281
 282   fSignalLamLamPiMinusSep = new TH2F ("fSignalLamLamPiMinusSep","PiMinus pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 283   fOutputList->Add(fSignalLamLamPiMinusSep);
 284
 285   fSignalALamALamAntiProtSep = new TH2F ("fSignalALamALamAntiProtSep","AntiProton pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 286   fOutputList->Add(fSignalALamALamAntiProtSep);
 287
 288   fSignalALamALamPiPlusSep = new TH2F ("fSignalALamALamPiPlusSep","PiPlus pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 289   fOutputList->Add(fSignalALamALamPiPlusSep);
 290
 291   fSignalLamALamAntiProtPiMinusSep = new TH2F ("fSignalLamALamAntiProtPiMinusSep","Neg particle pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 292   fOutputList->Add(fSignalLamALamAntiProtPiMinusSep);
 293
 294   fSignalLamALamProtPiPlusSep = new TH2F ("fSignalLamALamProtPiPlusSep","Pos pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 295   fOutputList->Add(fSignalLamALamProtPiPlusSep);
 296
 297   fBkgLamLamProtSep = new TH2F ("fBkgLamLamProtSep","Proton pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 298   fOutputList->Add(fBkgLamLamProtSep);
 299
 300   fBkgLamLamPiMinusSep = new TH2F ("fBkgLamLamPiMinusSep","PiMinus pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 301   fOutputList->Add(fBkgLamLamPiMinusSep);
 302
 303   fBkgALamALamAntiProtSep = new TH2F ("fBkgALamALamAntiProtSep","AntiProton pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 304   fOutputList->Add(fBkgALamALamAntiProtSep);
 305
 306   fBkgALamALamPiPlusSep = new TH2F ("fBkgALamALamPiPlusSep","PiPlus pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 307   fOutputList->Add(fBkgALamALamPiPlusSep);
 308
 309   fBkgLamALamAntiProtPiMinusSep = new TH2F ("fBkgLamALamAntiProtPiMinusSep","Neg particle pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 310   fOutputList->Add(fBkgLamALamAntiProtPiMinusSep);
 311
 312   fBkgLamALamProtPiPlusSep = new TH2F ("fBkgLamALamProtPiPlusSep","Pos pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 313   fOutputList->Add(fBkgLamALamProtPiPlusSep);
 314
 315
 316
 317   //opposite charged pair separation
 318   fSignalLamLamPlusMinusSep = new TH2F ("fSignalLamLamPlusMinusSep","Proton Pion pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 319   fOutputList->Add(fSignalLamLamPlusMinusSep);
 320
 321   fSignalALamALamPlusMinusSep = new TH2F ("fSignalALamALamPlusMinusSep","Proton Pion pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 322   fOutputList->Add(fSignalALamALamPlusMinusSep);
 323
 324   fSignalLamALamProtSep = new TH2F ("fSignalLamALamProtSep","Proton pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 325   fOutputList->Add(fSignalLamALamProtSep);
 326
 327   fSignalLamALamPionSep = new TH2F ("fSignalLamALamPionSep","Pion pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 328   fOutputList->Add(fSignalLamALamPionSep);
 329
 330
 331   fBkgLamLamPlusMinusSep = new TH2F ("fBkgLamLamPlusMinusSep","Proton Pion pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 332   fOutputList->Add(fBkgLamLamPlusMinusSep);
 333
 334   fBkgALamALamPlusMinusSep = new TH2F ("fBkgALamALamPlusMinusSep","Proton Pion pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 335   fOutputList->Add(fBkgALamALamPlusMinusSep);
 336
 337   fBkgLamALamProtSep = new TH2F ("fBkgLamALamProtSep","Proton pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 338   fOutputList->Add(fBkgLamALamProtSep);
 339
 340   fBkgLamALamPionSep = new TH2F ("fBkgLamALamPionSep","Pion pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 341   fOutputList->Add(fBkgLamALamPionSep);
 342
 343
 344
 345
 346   //primary-vertex corrected pair separation
 347
 348   fSignalLamLamProtSepCorrected = new TH2F ("fSignalLamLamProtSepCorrected","Proton pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 349   fOutputList->Add(fSignalLamLamProtSepCorrected);
 350
 351   fSignalLamLamPiMinusSepCorrected = new TH2F ("fSignalLamLamPiMinusSepCorrected","PiMinus pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 352   fOutputList->Add(fSignalLamLamPiMinusSepCorrected);
 353
 354   fSignalALamALamAntiProtSepCorrected = new TH2F ("fSignalALamALamAntiProtSepCorrected","AntiProton pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 355   fOutputList->Add(fSignalALamALamAntiProtSepCorrected);
 356
 357   fSignalALamALamPiPlusSepCorrected = new TH2F ("fSignalALamALamPiPlusSepCorrected","PiPlus pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 358   fOutputList->Add(fSignalALamALamPiPlusSepCorrected);
 359
 360   fSignalLamALamAntiProtPiMinusSepCorrected = new TH2F ("fSignalLamALamAntiProtPiMinusSepCorrected","Neg particle pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 361   fOutputList->Add(fSignalLamALamAntiProtPiMinusSepCorrected);
 362
 363   fSignalLamALamProtPiPlusSepCorrected = new TH2F ("fSignalLamALamProtPiPlusSepCorrected","Pos pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 364   fOutputList->Add(fSignalLamALamProtPiPlusSepCorrected);
 365
 366   fBkgLamLamProtSepCorrected = new TH2F ("fBkgLamLamProtSepCorrected","Proton pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 367   fOutputList->Add(fBkgLamLamProtSepCorrected);
 368
 369   fBkgLamLamPiMinusSepCorrected = new TH2F ("fBkgLamLamPiMinusSepCorrected","PiMinus pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 370   fOutputList->Add(fBkgLamLamPiMinusSepCorrected);
 371
 372   fBkgALamALamAntiProtSepCorrected = new TH2F ("fBkgALamALamAntiProtSepCorrected","AntiProton pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 373   fOutputList->Add(fBkgALamALamAntiProtSepCorrected);
 374
 375   fBkgALamALamPiPlusSepCorrected = new TH2F ("fBkgALamALamPiPlusSepCorrected","PiPlus pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 376   fOutputList->Add(fBkgALamALamPiPlusSepCorrected);
 377
 378   fBkgLamALamAntiProtPiMinusSepCorrected = new TH2F ("fBkgLamALamAntiProtPiMinusSepCorrected","Neg particle pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 379   fOutputList->Add(fBkgLamALamAntiProtPiMinusSepCorrected);
 380
 381   fBkgLamALamProtPiPlusSepCorrected = new TH2F ("fBkgLamALamProtPiPlusSepCorrected","Pos pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 382   fOutputList->Add(fBkgLamALamProtPiPlusSepCorrected);
 383
 384   //opposite charged pair separation corrected for primary vertex
 385   fSignalLamLamPlusMinusSepCorrected = new TH2F ("fSignalLamLamPlusMinusSepCorrected","Proton Pion pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 386   fOutputList->Add(fSignalLamLamPlusMinusSepCorrected);
 387
 388   fSignalALamALamPlusMinusSepCorrected = new TH2F ("fSignalALamALamPlusMinusSepCorrected","Proton Pion pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 389   fOutputList->Add(fSignalALamALamPlusMinusSepCorrected);
 390
 391   fSignalLamALamProtSepCorrected = new TH2F ("fSignalLamALamProtSepCorrected","Proton pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 392   fOutputList->Add(fSignalLamALamProtSepCorrected);
 393
 394   fSignalLamALamPionSepCorrected = new TH2F ("fSignalLamALamPionSepCorrected","Pion pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 395   fOutputList->Add(fSignalLamALamPionSepCorrected);
 396
 397
 398   fBkgLamLamPlusMinusSepCorrected = new TH2F ("fBkgLamLamPlusMinusSepCorrected","Proton Pion pair sep for Lam-Lam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 399   fOutputList->Add(fBkgLamLamPlusMinusSepCorrected);
 400
 401   fBkgALamALamPlusMinusSepCorrected = new TH2F ("fBkgALamALamPlusMinusSepCorrected","Proton Pion pair sep for ALam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 402   fOutputList->Add(fBkgALamALamPlusMinusSepCorrected);
 403
 404   fBkgLamALamProtSepCorrected = new TH2F ("fBkgLamALamProtSepCorrected","Proton pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 405   fOutputList->Add(fBkgLamALamProtSepCorrected);
 406
 407   fBkgLamALamPionSepCorrected = new TH2F ("fBkgLamALamPionSepCorrected","Pion pair sep for Lam-ALam", avgSepBins, 0., avgSepMaxValue, kStarBins, 0., maxKStar);
 408   fOutputList->Add(fBkgLamALamPionSepCorrected);
 409
 410
 411   PostData(1, fOutputList);
 412 }
 413
 414 //________________________________________________________________________
 415 void AliAnalysisV0Lam::Exec(Option_t *)
 416 {
 417   // Main loop
 418   // Called for each event
 419
 420   //Make sure we are using the correct event triggers
 421   if(!IsCorrectEventTrigger()) return;
 422   //  cout<<"===========  Event # "<<fEventCount+1<<"  ==========="<<endl;
 423   fEventCount++;
 424   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
 425   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
 426   AliAODVertex *primaryVertexAOD;
 427   double vertex[3]={0};
 428   int zBin=0;
 429   double zStep=2*10/double(zVertexBins), zStart=-10.;
 430   //Get Monte Carlo data if available
 431   TClonesArray *mcArray = 0x0;
 432   mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 433   //Check for MC event headers and get # of original Hijing particles
 434   //This is used later for rejecting injected signals
 435   Int_t numberOfLastHijingLabel = 0;
 436   AliAODMCHeader *mcHeader = (AliAODMCHeader*)fAOD->FindListObject(AliAODMCHeader::StdBranchName());
 437   if(mcHeader){
 438     // Get the iterator on the list of cocktail headers
 439     TIter next(mcHeader->GetCocktailHeaders());
 440     // Loop over the cocktail headers
 441     while (const TObject *obj=next()){
 442       // Check whether it's a Hijing header
 443       const AliGenHijingEventHeader* hijingHeader = dynamic_cast<const AliGenHijingEventHeader*>(obj);
 444       if(hijingHeader) {
 445         numberOfLastHijingLabel=hijingHeader->NProduced()-1;
 446         // We're done!
 447         //printf("Good! Last Hijing particle label %d\n",numberOfLastHijingLabel);
 448       } // End of found the hijing header
 449     } // End of loop over cocktail headers
 450   } // End of MC header exists
 451   else if(mcArray) cerr<<"Could not find mcHeader!"<<endl;
 452   if(mcArray) fIsMCEvent = kTRUE;
 453   else fIsMCEvent = kFALSE;
 454
 455   //Centrality selection
 456   AliCentrality *centrality = fAOD->GetCentrality();
 457   float centralityPercentile = centrality->GetCentralityPercentile("V0M");
 458   int centralityBin=0;
 459   //Printf("Centrality percent = %f", centralityPercentile);
 460   fCentrality->Fill(centralityPercentile);
 461   AliAODVZERO *aodV0 = fAOD->GetVZEROData();
 462   Float_t multV0A=aodV0->GetMTotV0A();
 463   Float_t multV0C=aodV0->GetMTotV0C();
 464   if(centralityPercentile < 0) {
 465     //Printf("No centrality info");
 466     return;}
 467   else if(centralityPercentile == 0 && (multV0A + multV0C < 19500)) {
 468     //Printf("No centrality info");
 469     return;}
 470   // Centrality info is good.  Now find the correct 5% centrality bin
 471   else if(centralityPercentile <= 5.) centralityBin=19;
 472   else if(centralityPercentile <= 10.) centralityBin=18;
 473   else if(centralityPercentile <= 15.) centralityBin=17;
 474   else if(centralityPercentile <= 20.) centralityBin=16;
 475   else if(centralityPercentile <= 25.) centralityBin=15;
 476   else if(centralityPercentile <= 30.) centralityBin=14;
 477   else if(centralityPercentile <= 35.) centralityBin=13;
 478   else if(centralityPercentile <= 40.) centralityBin=12;
 479   else if(centralityPercentile <= 45.) centralityBin=11;
 480   else if(centralityPercentile <= 50.) centralityBin=10;
 481   else if(centralityPercentile <= 55.) centralityBin=9;
 482   else if(centralityPercentile <= 60.) centralityBin=8;
 483   else if(centralityPercentile <= 65.) centralityBin=7;
 484   else if(centralityPercentile <= 70.) centralityBin=6;
 485   else if(centralityPercentile <= 75.) centralityBin=5;
 486   else if(centralityPercentile <= 80.) centralityBin=4;
 487   else if(centralityPercentile <= 85.) centralityBin=3;
 488   else if(centralityPercentile <= 90.) centralityBin=2;
 489   else if(centralityPercentile <= 95.) centralityBin=1;
 490   else if(centralityPercentile <= 100.) centralityBin=0;
 491   else {Printf("Skipping Peripheral Event"); return;}
 492
 493   //Vertexing
 494   primaryVertexAOD = fAOD->GetPrimaryVertex();
 495   vertex[0]=primaryVertexAOD->GetX();
 496   vertex[1]=primaryVertexAOD->GetY();
 497   vertex[2]=primaryVertexAOD->GetZ();
 498   if(vertex[0]<10e-5 && vertex[1]<10e-5 &&  vertex[2]<10e-5) return;
 499   if(fabs(vertex[2]) > 10) return; // Z-Vertex Cut
 500   for(int i=0; i<zVertexBins; i++)
 501   {
 502     if((vertex[2] > zStart+i*zStep) && (vertex[2] < zStart+(i+1)*zStep))
 503     {
 504       zBin=i;
 505       break;
 506     }
 507   }
 508   double bfield = fAOD->GetMagneticField();
 509   //Printf("Rough multiplicity = %d", fAOD->GetNumberOfTracks());
 510   fMultDistRough->Fill(fAOD->GetNumberOfTracks());
 511   /////////////////////////////////////////////////////////
 512   //Add Event to buffer - this is for event mixing
 513   fEC[zBin][centralityBin]->FifoShift();
 514   fEvt = fEC[zBin][centralityBin]->fEvt;
 515   fCutProcessor->SetCentralityBin(centralityBin+1);
 516 //////////////////////////////////////////////////////////////////
 517   //v0 tester
 518 ////////////////////////////////////////////////////////////////
 519
 520   int v0Count = 0;
 521   vector<int> lambdaCount(fNumberOfVariableCutValues,0);
 522   vector<int> antiLambdaCount(fNumberOfVariableCutValues,0);
 523   TH1F *mcTruthOriginHist = nullptr;
 524   if(fIsMCEvent) mcTruthOriginHist = CreateLambdaOriginHist(mcArray,numberOfLastHijingLabel); //Find MC truths of all the MC particles (before detector effects)
 525   TH1F *v0OriginHist = new TH1F("v0OriginHist", "Lambda Origins", AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1);
 526   TH2F *v0PassedCutsOriginHist = new TH2F("v0PassedCutsOriginHist", "Lambda Origins", AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1, fNumberOfVariableCutValues, -0.5, fNumberOfVariableCutValues -0.5);
 527   for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
 528   {
 529     //Loop over all the V0 candidates to look for (anti)Lambdas
 530     bool hasPiPlusDaughter     = kFALSE;
 531     bool hasPiMinusDaughter    = kFALSE;
 532     bool hasProtonDaughter     = kFALSE;
 533     bool hasAntiProtonDaughter = kFALSE;
 534     AliAODv0* v0 = fAOD->GetV0(i);
 535     if(!v0) continue;
 536     //Make sure the V0 satisifies a few basic criteria
 537     if(v0->GetNDaughters() > 2)                          continue;
 538     if(v0->GetNProngs() > 2)                             continue;
 539     if(v0->GetCharge() != 0)                             continue;
 540     if(v0->ChargeProng(0) == v0->ChargeProng(1))         continue;
 541     if(v0->GetOnFlyStatus())                             continue;
 542     //Now look at daughter tracks
 543     AliAODTrack* daughterTrackPos = (AliAODTrack*)v0->GetDaughter(0);
 544     AliAODTrack* daughterTrackNeg = (AliAODTrack*)v0->GetDaughter(1);
 545     if(!daughterTrackPos) continue; //Daughter tracks must exist
 546     if(!daughterTrackNeg) continue;
 547     daughterTrackPos->SetAODEvent(fAOD); //Need to set this for PID purposes
 548     daughterTrackNeg->SetAODEvent(fAOD);
 549     //Need to manually apply ITS refit cut for hybrid tracks in AOD115
 550     if((daughterTrackPos->GetStatus() & AliVTrack::kTPCrefit)==0) continue;
 551     if((daughterTrackNeg->GetStatus() & AliVTrack::kTPCrefit)==0) continue;
 552     AliPIDResponse::EDetPidStatus statusPosTPC = fpidAOD->CheckPIDStatus(AliPIDResponse::kTPC,daughterTrackPos);
 553     AliPIDResponse::EDetPidStatus statusNegTPC = fpidAOD->CheckPIDStatus(AliPIDResponse::kTPC,daughterTrackNeg);
 554     if(AliPIDResponse::kDetPidOk != statusPosTPC) continue;
 555     if(AliPIDResponse::kDetPidOk != statusNegTPC) continue;
 556     if(daughterTrackPos->GetTPCNcls() < 80) continue;
 557     if(daughterTrackNeg->GetTPCNcls() < 80) continue;
 558     //Need to manually apply shared cluster cut for hybrid tracks in AOD115 and 124
 559     Double_t fracSharedClustersPos = Double_t(daughterTrackPos->GetTPCnclsS()) / Double_t(daughterTrackPos->GetTPCncls());
 560     Double_t fracSharedClustersNeg = Double_t(daughterTrackNeg->GetTPCnclsS()) / Double_t(daughterTrackNeg->GetTPCncls());
 561     if(fracSharedClustersPos  > 0.4) continue;
 562     if(fracSharedClustersNeg  > 0.4) continue;
 563     if(daughterTrackPos->Pt() < .16) continue;
 564     if(daughterTrackNeg->Pt() < .16) continue;
 565     if(fabs(daughterTrackPos->Eta()) > fEtaDaughter) continue;
 566     if(fabs(daughterTrackNeg->Eta()) > fEtaDaughter) continue;
 567
 568     //Now we'll get particle origin and momentum truth for MC particles
 569     AliReconstructedV0::MCV0Origin_t mcV0Origin = AliReconstructedV0::kUnassigned;
 570     double v0MomentumTruth[3] = {0.,0.,0.};
 571     if(fIsMCEvent){
 572       //first reject injected particles
 573       if(IsInjectedParticle(v0,mcArray,numberOfLastHijingLabel)) continue;
 574       mcV0Origin = DetermineV0Origin(v0, mcArray);
 575       v0OriginHist->Fill(mcV0Origin);
 576       GetMCParticleMomentumTruth(v0MomentumTruth, v0, mcArray);
 577     }
 578     //Now perform daughter track PID using TPC
 579     if (daughterTrackPos->Pt() > 0.5) {// min-pt cut to fix proton PID issues in LHC10h data
 580       if(fabs(fpidAOD->NumberOfSigmasTPC(daughterTrackPos,AliPID::kProton))
 581          < fSigmaCutTPCProton)   hasProtonDaughter = kTRUE;
 582     }
 583     if(fabs(fpidAOD->NumberOfSigmasTPC(daughterTrackPos,AliPID::kPion)) < fSigmaCutTPCPion)  hasPiPlusDaughter = kTRUE;
 584     if (daughterTrackNeg->Pt() > 0.3){
 585       if(fabs(fpidAOD->NumberOfSigmasTPC(daughterTrackNeg,AliPID::kProton))
 586          < fSigmaCutTPCProton) hasAntiProtonDaughter = kTRUE;
 587     }
 588     if(fabs(fpidAOD->NumberOfSigmasTPC(daughterTrackNeg,AliPID::kPion)) < fSigmaCutTPCPion) hasPiMinusDaughter = kTRUE;
 589     //Use TOF PID info if available.  This overrides TPC PID results.
 590     if(daughterTrackPos->P() > fTOFLow)
 591     { // positive daughter PID
 592       AliPIDResponse::EDetPidStatus statusPosTOF = fpidAOD->CheckPIDStatus(AliPIDResponse::kTOF,daughterTrackPos);
 593       if (AliPIDResponse::kDetPidOk == statusPosTOF) { // TOF signal is available for PID
 594         Float_t probMis = fpidAOD->GetTOFMismatchProbability(daughterTrackPos);
 595         if (probMis < 0.01) { // avoid TOF-TPC mismatch
 596           //PiPlus
 597           if(fabs(fpidAOD->NumberOfSigmasTOF(daughterTrackPos,AliPID::kPion)) < fSigmaCutTOFPion) hasPiPlusDaughter = kTRUE;
 598           else hasPiPlusDaughter = kFALSE;
 599           //Proton
 600           if(fabs(fpidAOD->NumberOfSigmasTOF(daughterTrackPos,AliPID::kProton)) < fSigmaCutTOFProton) hasProtonDaughter = kTRUE;
 601           else hasProtonDaughter = kFALSE;
 602         }
 603       }
 604     }
 605     if(daughterTrackNeg->P() > fTOFLow)
 606     { //negative daughter PID
 607       AliPIDResponse::EDetPidStatus statusNegTOF = fpidAOD->CheckPIDStatus(AliPIDResponse::kTOF,daughterTrackNeg);
 608       if (AliPIDResponse::kDetPidOk == statusNegTOF) { // TOF signal is available for PID
 609         Float_t probMis = fpidAOD->GetTOFMismatchProbability(daughterTrackNeg);
 610         if (probMis < 0.01) {
 611           //PiMinus
 612           if(fabs(fpidAOD->NumberOfSigmasTOF(daughterTrackNeg,AliPID::kPion)) < fSigmaCutTOFPion) hasPiMinusDaughter = kTRUE;
 613           else hasPiMinusDaughter = kFALSE;
 614           //AntiProton
 615           if(fabs(fpidAOD->NumberOfSigmasTOF(daughterTrackNeg,AliPID::kProton)) < fSigmaCutTOFProton) hasAntiProtonDaughter = kTRUE;
 616           else hasAntiProtonDaughter = kFALSE;
 617         }
 618       }
 619     }
 620     //If V0 doesn't have the right daughter combinations,
 621     //move on to the next candidate
 622     if(!((hasProtonDaughter && hasPiMinusDaughter) || (hasAntiProtonDaughter && hasPiPlusDaughter))) continue;
 623
 624     //Save V0 information
 625     fEvt->fReconstructedV0[v0Count].v0Momentum[0]  = v0->Px();
 626     fEvt->fReconstructedV0[v0Count].v0Momentum[1]  = v0->Py();
 627     fEvt->fReconstructedV0[v0Count].v0Momentum[2]  = v0->Pz();
 628     for(int pIndex = 0; pIndex <3; pIndex++)
 629     {
 630       fEvt->fReconstructedV0[v0Count].v0MomentumTruth[pIndex]  = v0MomentumTruth[pIndex];
 631     }
 632     fEvt->fReconstructedV0[v0Count].v0Pt     = v0->Pt();
 633     fEvt->fReconstructedV0[v0Count].v0Eta    = v0->Eta();
 634     fEvt->fReconstructedV0[v0Count].v0Phi    = v0->Phi();
 635     fEvt->fReconstructedV0[v0Count].massLam  = v0->MassLambda();
 636     fEvt->fReconstructedV0[v0Count].massALam = v0->MassAntiLambda();
 637     fEvt->fReconstructedV0[v0Count].massLamDifference  = fabs(v0->MassLambda() - fPDGLambda);
 638     fEvt->fReconstructedV0[v0Count].massALamDifference = fabs(v0->MassAntiLambda() - fPDGLambda);
 639     fEvt->fReconstructedV0[v0Count].massK0  = v0->MassK0Short();
 640     fEvt->fReconstructedV0[v0Count].lorentzGammaLam = v0->ELambda()/fPDGLambda;
 641     fEvt->fReconstructedV0[v0Count].v0DCA   = v0->DcaV0ToPrimVertex();
 642     fEvt->fReconstructedV0[v0Count].decayLength = v0->DecayLength(primaryVertexAOD);
 643     fEvt->fReconstructedV0[v0Count].decayVertexPosition[0] = v0->DecayVertexV0X();
 644     fEvt->fReconstructedV0[v0Count].decayVertexPosition[1] = v0->DecayVertexV0Y();
 645     fEvt->fReconstructedV0[v0Count].decayVertexPosition[2] = v0->DecayVertexV0Z();
 646     fEvt->fReconstructedV0[v0Count].cosPointing = v0->CosPointingAngle(primaryVertexAOD);
 647     fEvt->fReconstructedV0[v0Count].hasProtonDaughter     = hasProtonDaughter;
 648     fEvt->fReconstructedV0[v0Count].hasAntiProtonDaughter = hasAntiProtonDaughter;
 649     fEvt->fReconstructedV0[v0Count].hasPiPlusDaughter     = hasPiPlusDaughter;
 650     fEvt->fReconstructedV0[v0Count].hasPiMinusDaughter    = hasPiMinusDaughter;
 651     fEvt->fReconstructedV0[v0Count].mcOriginType = mcV0Origin;
 652     //Save Daughter information
 653     fEvt->fReconstructedV0[v0Count].daughter1ID = v0->GetNegID();
 654     fEvt->fReconstructedV0[v0Count].daughter2ID = v0->GetPosID();
 655     fEvt->fReconstructedV0[v0Count].daughterPosMomentum[0] = v0->MomPosX();
 656     fEvt->fReconstructedV0[v0Count].daughterPosMomentum[1] = v0->MomPosY();
 657     fEvt->fReconstructedV0[v0Count].daughterPosMomentum[2] = v0->MomPosZ();
 658     fEvt->fReconstructedV0[v0Count].daughterPosProtonE = v0->EPosProton();
 659     fEvt->fReconstructedV0[v0Count].daughterPosPionE = v0->EPosPion();
 660     fEvt->fReconstructedV0[v0Count].daughterNegMomentum[0] = v0->MomNegX();
 661     fEvt->fReconstructedV0[v0Count].daughterNegMomentum[1] = v0->MomNegY();
 662     fEvt->fReconstructedV0[v0Count].daughterNegMomentum[2] = v0->MomNegZ();
 663     fEvt->fReconstructedV0[v0Count].daughterNegProtonE = v0->ENegProton();
 664     fEvt->fReconstructedV0[v0Count].daughterNegPionE = v0->ENegPion();
 665     fEvt->fReconstructedV0[v0Count].daughtersDCA = v0->DcaV0Daughters();
 666     fEvt->fReconstructedV0[v0Count].daughterPosDCAPrimaryVertex = v0->DcaPosToPrimVertex();
 667     fEvt->fReconstructedV0[v0Count].daughterNegDCAPrimaryVertex = v0->DcaNegToPrimVertex();
 668     GetGlobalPositionAtGlobalRadiiThroughTPC(daughterTrackPos, bfield, fEvt->fReconstructedV0[v0Count].daughterPosGlobalPositions);// used for merging cuts later
 669     GetGlobalPositionAtGlobalRadiiThroughTPC(daughterTrackNeg, bfield, fEvt->fReconstructedV0[v0Count].daughterNegGlobalPositions);
 670     daughterTrackPos->GetXYZ(fEvt->fReconstructedV0[v0Count].daughterPosPositionDCA);
 671     daughterTrackNeg->GetXYZ(fEvt->fReconstructedV0[v0Count].daughterNegPositionDCA);
 672     daughterTrackPos->GetPxPyPz(fEvt->fReconstructedV0[v0Count].daughterPosMomentumDCA);
 673     daughterTrackNeg->GetPxPyPz(fEvt->fReconstructedV0[v0Count].daughterNegMomentumDCA);
 674     daughterTrackPos->GetCovarianceXYZPxPyPz(fEvt->fReconstructedV0[v0Count].daughterPosCovariance);
 675     daughterTrackNeg->GetCovarianceXYZPxPyPz(fEvt->fReconstructedV0[v0Count].daughterNegCovariance);
 676     for(int coord = 0; coord <3; coord++){
 677       fEvt->fPrimaryVertex[coord] = vertex[coord];
 678       for(int location = 0; location <9; location++) {
 679         //find the track locations relative to the primary vertex location.
 680         fEvt->fReconstructedV0[v0Count].daughterPosCorrectedGlobalPositions[location][coord] = fEvt->fReconstructedV0[v0Count].daughterPosGlobalPositions[location][coord] - vertex[coord];
 681         fEvt->fReconstructedV0[v0Count].daughterNegCorrectedGlobalPositions[location][coord] = fEvt->fReconstructedV0[v0Count].daughterNegGlobalPositions[location][coord] - vertex[coord];
 682       }
 683     }
 684
 685     //Now analyze and histogram the V0
 686     fCutProcessor->CheckIfV0PassesCuts(& fEvt->fReconstructedV0[v0Count]);
 687     fCutProcessor->DoV0Histogramming(& fEvt->fReconstructedV0[v0Count]);
 688     FillReconstructedV0MCOrigin(& fEvt->fReconstructedV0[v0Count], v0PassedCutsOriginHist);
 689     AddV0ToMultiplicityCounts(& fEvt->fReconstructedV0[v0Count], lambdaCount, antiLambdaCount);
 690     FillTPCSignalHists(& fEvt->fReconstructedV0[v0Count], daughterTrackPos->P(), daughterTrackPos->GetTPCsignal(), daughterTrackNeg->P(), daughterTrackNeg->GetTPCsignal());
 691     if(fIsMCEvent) CheckForFakeV0s(& fEvt->fReconstructedV0[v0Count], fMCFakeParticleIdentity, fMCOtherV0Identity, mcV0Origin);
 692
 693     //Increment V0 count and check that we don't exceed size of V0 array
 694     v0Count++;
 695     if(fMaxV0Mult <= v0Count){
 696       cerr<<"V0 Count has exceeded"<<fMaxV0Mult<<"!"<<endl;
 697       break;
 698     }
 699   } //End of V0 loop
 700   //cout<<"Finished with V0 storage.  V0 candidate count is "<<v0Count<<endl;
 701
 702   fEvt->fNumberCandidateV0 = v0Count;
 703   if(fIsMCEvent) BinOriginInformationForMCParticles(mcTruthOriginHist, v0OriginHist, v0PassedCutsOriginHist);
 704
 705   //The following histograms don't get used again, so clean them up
 706   if(mcTruthOriginHist){
 707     delete mcTruthOriginHist;
 708     mcTruthOriginHist = nullptr;
 709   }
 710   if(v0OriginHist){
 711     delete v0OriginHist;
 712     v0OriginHist = nullptr;
 713   }
 714   if(v0PassedCutsOriginHist){
 715     delete v0PassedCutsOriginHist;
 716     v0PassedCutsOriginHist = nullptr;
 717   }
 718
 719
 720   DoV0JudgmentCuts(fEvt, v0Count);
 721   HistogramEventMultiplicities(lambdaCount, antiLambdaCount, centralityBin);
 722   fTotalLambda += lambdaCount[fDefaultVariableCutIndex];
 723   fTotalAntiLambda += antiLambdaCount[fDefaultVariableCutIndex];
 724   //Printf("Reconstruction Finished. Starting pair studies.");
 725
 726   //Now look at pairs for correlation function binning
 727   DoPairStudies(fEvt, centralityBin);
 728   //cout<<"Pair studies completed.  Event finished"<<endl;
 729
 730   // Post output data.
 731   PostData(1, fOutputList);
 732 }
 733 //________________________________________________________________________
 734 void AliAnalysisV0Lam::Terminate(Option_t *)
 735 {
 736   // Called once at the end of the query
 737   cout<<"Total Lambdas found:\t"<<fTotalLambda<<"."<<endl
 738       <<"Total AntiLambdas found:\t"<<fTotalAntiLambda<<"."<<endl
 739       <<"Done"<<endl;
 740 }
 741
 742
 743
 744
 745
 746
 747 //________________________________________________________________________
 748 void AliAnalysisV0Lam::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3])
 749 {
 750   // Gets the global position of the track at nine different radii in the TPC
 751   // track is the track you want to propagate
 752   // bfield is the magnetic field of your event
 753   //globalPositionsAtRadii is the array of global positions in the radii and xyz
 754   // Initialize the array to something indicating there was no propagation
 755   for(Int_t i=0;i<9;i++){
 756     for(Int_t j=0;j<3;j++){
 757       globalPositionsAtRadii[i][j]=-9999.;
 758     }
 759   }
 760   // Make a copy of the track to not change parameters of the track
 761   AliExternalTrackParam etp;
 762   etp.CopyFromVTrack(track);
 763   //printf("\nAfter CopyFromVTrack\n");
 764   //etp.Print();
 765   // The global position of the the track
 766   Double_t xyz[3]={-9999.,-9999.,-9999.};
 767   // Counter for which radius we want
 768   Int_t iR=0;
 769   // The radii at which we get the global positions
 770   // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
 771   Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
 772   // The global radius we are at
 773   Float_t globalRadius=0;
 774   // Propagation is done in local x of the track
 775   for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
 776     // Starts at the tracks fX and goes outwards. x = 245 is the outer radial
 777     // limit of the TPC when the track is straight, i.e. has inifinite pt
 778     // and doesn't get bent. If the track's momentum is smaller than infinite,
 779     // it will develop a y-component, which adds to the global radius
 780     // Stop if the propagation was not succesful. This can happen for low pt
 781     // tracks that don't reach outer radii
 782     if(!etp.PropagateTo(x,bfield))break;
 783     etp.GetXYZ(xyz); // GetXYZ returns global coordinates
 784     globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
 785     // Roughly reached the radius we want
 786     if(globalRadius > Rwanted[iR]){
 787       // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
 788       while (globalRadius>Rwanted[iR]){
 789         x-=.1;
 790         //      printf("propagating to x %5.2f\n",x);
 791         if(!etp.PropagateTo(x,bfield))break;
 792         etp.GetXYZ(xyz); // GetXYZ returns global coordinates
 793         globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
 794       }
 795       //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
 796       globalPositionsAtRadii[iR][0]=xyz[0];
 797       globalPositionsAtRadii[iR][1]=xyz[1];
 798       globalPositionsAtRadii[iR][2]=xyz[2];
 799       // Indicate we want the next radius
 800       iR+=1;
 801     }
 802     if(iR>=8){
 803       // TPC edge reached
 804       return;
 805     }
 806   }
 807 }
 808
 809
 810
 811 //________________________________________________________________________
 812 double AliAnalysisV0Lam::GetAverageSeparation(Float_t globalPositions1st[9][3], Float_t globalPositions2nd[9][3])
 813 {
 814   //Compute the separation of two daughter tracks, averaged over 9 different positions
 815   double avgSeparation = 0;
 816   for(int RadiusNumber = 0; RadiusNumber <9; RadiusNumber++){
 817     double sumsquare=0;
 818     for(int Component = 0; Component <3; Component++){
 819       sumsquare+= pow(globalPositions1st[RadiusNumber][Component]-globalPositions2nd[RadiusNumber][Component],2);
 820     }
 821     avgSeparation+=sqrt(sumsquare);
 822   }
 823   return avgSeparation/9.;
 824 }
 825
 826 //________________________________________________________________________
 827 void AliAnalysisV0Lam::DoV0JudgmentCuts(const AliAnalysisV0LamEvent * const event, const int totalV0s)
 828 {
 829   // Looks at all V0s in a given event, and selectively removes V0s such
 830   // that each daughter track is claimed by no more than one V0.  This is
 831   // done by making judgment cuts.  The judgment cut compares a
 832   // characteristic (e.g. cosine of pointing angle) of two V0s that share a
 833   // daughter.  The V0 closer to the ideal value (e.g. cos(pointing) = 1) is
 834   // kept, while the other V0 is removed.
 835
 836   // Occasionally several V0s will share a single daughter, or several V0s
 837   // will share several daughters.  Because of this, the possibility
 838   // exists for a V0 to be removed, and subsequently the V0 that it
 839   // competed with is also removed.  In that case the original V0 should be
 840   // "restored". This DoV0JudgmentCuts method includes an iterative process
 841   // which first removes V0s that fail the judgment cuts, and then
 842   // subsequently restores V0s which no longer compete with any other V0s.
 843   // This process of removing and restoring V0s continues until the
 844   // process stabilizes or 20 iterations occurs.
 845
 846   //"Removed" V0s have a boolean flag "isDeemedUnworthy" which is set to
 847   // true.  Those V0s do not get used in correlation function pairs.
 848
 849   // The selectionCriterion is used to set which of V0 DCA, daughter DCA to
 850   // each other, V0 cosine of pointing angle, or V0 mass is used as the
 851   // judgment cut.
 852   const int selectionCriterion = 0;
 853   // Start by looping over variable reconstruction cuts.  There will be
 854   // different lists of V0s for each reconstruction cut value, which may
 855   // lead to different sets of V0s competing over daughters.
 856   for (int cutIndex = 0; cutIndex < fNumberOfVariableCutValues; cutIndex++){
 857     bool converged;
 858     int iterations = 0;
 859     do { //Loop until the judgment cuts converge or 20 iterations pass
 860       converged = kTRUE;
 861       iterations++;
 862       for (int currentV0Number = 0; currentV0Number <totalV0s; currentV0Number++) { // Loop over each V0 in event
 863         // Mark V0s as bad if they are worse than other V0s with shared
 864         // daughters
 865         if(!(event->fReconstructedV0[currentV0Number].isLamCenter[cutIndex] || event->fReconstructedV0[currentV0Number].isALamCenter[cutIndex])) continue;
 866         // Don't bother if the current v0 isn't a center V0 (a center V0 has an m_inv inside the accepted mass window)
 867         if(event->fReconstructedV0[currentV0Number].isDeemedUnworthy[cutIndex]) continue;
 868         for (int comparisonV0Number = 0; comparisonV0Number <totalV0s; comparisonV0Number++)
 869         { // Loop over all other V0s in event.
 870           if(comparisonV0Number == currentV0Number) continue;
 871           if(!(event->fReconstructedV0[comparisonV0Number].isLamCenter[cutIndex] || event->fReconstructedV0[comparisonV0Number].isALamCenter[cutIndex])) continue; //Don't bother if the comparison V0 isn't a center v0
 872           if(!((event->fReconstructedV0[currentV0Number].daughter1ID == event->fReconstructedV0[comparisonV0Number].daughter1ID) || (event->fReconstructedV0[currentV0Number].daughter2ID == event->fReconstructedV0[comparisonV0Number].daughter2ID))) continue; //Don't bother if they don't share daughters
 873           if(event->fReconstructedV0[comparisonV0Number].isDeemedUnworthy[cutIndex]) continue;
 874           // If we reach this point in the loop, then these V0s compete
 875           // over a daughter.  Compare them and determine which V0 needs
 876           // to be removed.
 877           int worseV0 = DetermineWhichV0IsWorse(event, currentV0Number, comparisonV0Number, selectionCriterion, cutIndex);
 878           if (worseV0 != -1) event->fReconstructedV0[worseV0].isDeemedUnworthy[cutIndex] = kTRUE;
 879           // A V0 has been removed, so process has not converged.
 880           converged = kFALSE;
 881           if(event->fReconstructedV0[currentV0Number].isDeemedUnworthy[cutIndex]) break;
 882         } //End loop over comparison V0s
 883       }// End V0 removal
 884       // Now Restore V0s if they no long compete over a daughter OR if they
 885       // compete and are judged to be good.
 886       for (int currentV0Number = 0; currentV0Number <totalV0s; currentV0Number++) { // Loop over each V0 in event
 887         if(!event->fReconstructedV0[currentV0Number].isDeemedUnworthy[cutIndex]) continue; //Only look at V0s that have been removed
 888         bool stillCompeting = kFALSE; // Assume that they no longer compete
 889         for (int comparisonV0Number = 0; comparisonV0Number < totalV0s; comparisonV0Number++)
 890         { // Loop over all other V0s
 891           if(comparisonV0Number == currentV0Number) continue;
 892           if(!(event->fReconstructedV0[comparisonV0Number].isLamCenter[cutIndex] || event->fReconstructedV0[comparisonV0Number].isALamCenter[cutIndex])) continue; //Don't bother if the V0 is a center V0
 893           if(event->fReconstructedV0[comparisonV0Number].isDeemedUnworthy[cutIndex]) continue; //Only compare with V0s that HAVE NOT been removed
 894           if(!((event->fReconstructedV0[currentV0Number].daughter1ID == event->fReconstructedV0[comparisonV0Number].daughter1ID) || (event->fReconstructedV0[currentV0Number].daughter2ID == event->fReconstructedV0[comparisonV0Number].daughter2ID))) continue; //Don't bother if they don't share daughters
 895           // If we reach this point in the loop, then these V0s compete
 896           // over a daughter.  Compare them and determine which V0 needs
 897           // to be removed.
 898           int worseV0 = DetermineWhichV0IsWorse(event, currentV0Number, comparisonV0Number, selectionCriterion, cutIndex);
 899           if(worseV0 == -1)
 900           {
 901             //Something has gone wrong.
 902             cerr<<"Could not determine which V0 is worse"<<endl;
 903           }
 904           else if (worseV0 == currentV0Number)
 905           {
 906             // The V0 is still competing with another V0, and it has failed
 907             // the judgment cut, so it stays removed.
 908             stillCompeting = kTRUE;
 909             break; //No need to keep comparing with other V0s
 910           }
 911           else {
 912             //The comparison V0 is worse.  However, do nothing here.
 913             //The comparison V0 will be removed in the removal loop
 914             //if it still competes at that time.
 915           }
 916         } //end comparison V0 loop
 917         if(!stillCompeting)
 918         {
 919           //The V0 no longer fails any judgment cuts.  Restore it.
 920           event->fReconstructedV0[currentV0Number].isDeemedUnworthy[cutIndex] = kFALSE;
 921           converged = kFALSE;
 922         }
 923       }//End V0 restoration.
 924     } while ((converged == kFALSE) && (iterations < 20));
 925   } //End loop over variable reconstruction cuts.
 926   return;
 927 }
 928
 929 //________________________________________________________________________
 930 int AliAnalysisV0Lam::DetermineWhichV0IsWorse(const AliAnalysisV0LamEvent * const event, const int V01, const int V02, const int Criterion, const int cutIndex)
 931 {
 932   // Performs a judgment cut on two V0 by comparing characteristics of those V0
 933   // and looking to see which of those V0 is further from the ideal value.
 934   // Cut is only performed on two V0 that claim the same daughter track.
 935   int worseV0 = -1;
 936   if (Criterion == 0)//compare using DCA to primary vertex
 937   {
 938     if (event->fReconstructedV0[V01].v0DCA <
 939         event->fReconstructedV0[V02].v0DCA) worseV0=V02;
 940     else worseV0=V01;
 941   }
 942   else if (Criterion == 1)//compare using DCA of daughters
 943   {
 944     if(event->fReconstructedV0[V01].daughtersDCA <
 945        event->fReconstructedV0[V02].daughtersDCA) worseV0 = V02;
 946     else worseV0 = V01;
 947   }
 948   else if (Criterion == 2)//compare using cos(pointing) of V0s
 949   {
 950     if(event->fReconstructedV0[V01].cosPointing >
 951        event->fReconstructedV0[V02].cosPointing) worseV0 = V02;
 952     else worseV0 = V01;
 953   }
 954
 955   else if (Criterion == 3)//compare using Minv
 956   {
 957     double deltaM1=500.;
 958     double deltaM2=500.;
 959     if(event->fReconstructedV0[V01].isLamCenter[cutIndex]){
 960       deltaM1 = event->fReconstructedV0[V01].massLamDifference;
 961     }
 962     else if(event->fReconstructedV0[V01].isALamCenter[cutIndex]){
 963       deltaM1 = event->fReconstructedV0[V01].massALamDifference;
 964     }
 965     if(event->fReconstructedV0[V02].isLamCenter[cutIndex]){
 966       deltaM2 = event->fReconstructedV0[V02].massLamDifference;
 967     }
 968     else if(event->fReconstructedV0[V02].isALamCenter[cutIndex]){
 969       deltaM2 = event->fReconstructedV0[V02].massALamDifference;
 970     }
 971     if(deltaM1 <= deltaM2) worseV0 = V02;
 972     else worseV0 = V01;
 973   }
 974   else cerr<<"Invalid judgment cut criterion selected: "<<Criterion<<endl;
 975   return worseV0;
 976 }
 977
 978 //________________________________________________________________________
 979 TH1F *AliAnalysisV0Lam::CreateLambdaOriginHist(TClonesArray *mcArray, Int_t numberOfLastHijingLabel)
 980 {
 981   //Create a histogram of the MC truth origin of each (anti)Lambda in the event
 982   //This allows us to see how many primary and secondary lambda there are in the
 983   //event.  We'll count again after all reconstruction is done to get an idea
 984   //of our reconstruction efficiency.
 985   TH1F *mcTruthOriginHist = new TH1F("mcTruthOriginHist", "Lambda Origins", AliReconstructedV0::kOriginTypeMax+1, 0, AliReconstructedV0::kOriginTypeMax+1);
 986   for (int i=0; i < mcArray->GetEntriesFast(); i++){
 987     AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(i);
 988     if(mcParticle->GetNDaughters() != 2) continue;
 989     //Reject injected particles.  Injected particles have a label greater
 990     //than numberOfLastHijingLabel.  Many secondary particles also have a label
 991     //greater than numberOfLastHijingLabel.  So first check if a particle has a
 992     //parent.  If it does, check if that parent is "original".  If not
 993     //original, reject that particle (e.g. this will reject lambdas with
 994     //injected cascade parents).  If a particle has no parent and it has
 995     //a label greater than numberOfLastHijingLabel, reject it (because it is
 996     //injected)
 997     if(mcParticle->GetMother() > -1){ // This MCParticle has a mother
 998       AliAODMCParticle *mcMother = (AliAODMCParticle*)mcArray->At(mcParticle->GetMother());
 999       if(!mcMother) continue;
1000       //Reject this MCParticle if its mother is injected
1001       if(mcMother->GetLabel() > numberOfLastHijingLabel) continue;
1002     }
1003     //If this MCParticle has no mother but it has a label > LastHijingLabel,
1004     //reject it
1005     else if(mcParticle->GetLabel() > numberOfLastHijingLabel) continue;
1006     AliAODMCParticle *mcDaughter1 = (AliAODMCParticle*)mcArray->At(mcParticle->GetDaughter(0));
1007     AliAODMCParticle *mcDaughter2 = (AliAODMCParticle*)mcArray->At(mcParticle->GetDaughter(1));
1008     //We won't count any MC Particles that have daughters outside the acceptance
1009     //region
1010     if(fabs(mcDaughter1->Eta()) > fEtaDaughter) continue;
1011     if(fabs(mcDaughter2->Eta()) > fEtaDaughter) continue;
1012     if(mcDaughter1->Pt() < 0.16) continue;
1013     if(mcDaughter2->Pt() < 0.16) continue;
1014     //Finally, get the PDG code of the MCParticle (or of its parent in the case
1015     //of secondary particles)
1016     AliReconstructedV0::MCV0Origin_t mcParticleOrigin = DeterminePdgCodeOfMcParticle(mcParticle,mcArray);
1017     mcTruthOriginHist->Fill(mcParticleOrigin);
1018   }
1019   return mcTruthOriginHist;
1020 }
1021
1022 //________________________________________________________________________
1023 void AliAnalysisV0Lam::FillReconstructedV0MCOrigin(const AliReconstructedV0 * v0, TH2F *histPassedCutsOrigin)
1024 {
1025   //Make a histogram showing the MCTruth particle type of reconstructed V0s
1026   //(or the type of the mother particle if the V0 is secondary).
1027   for(int i = 0; i < fNumberOfVariableCutValues; i++){
1028     if(v0->isLamCenter[i]
1029        && (AliReconstructedV0::kFake == v0->mcOriginType))
1030     {
1031       histPassedCutsOrigin->Fill(AliReconstructedV0::kFakeLambda,i);
1032     }
1033     else if(v0->isLamCenter[i]){
1034       histPassedCutsOrigin->Fill(v0->mcOriginType,i);
1035     }
1036     if(v0->isALamCenter[i]
1037        && (AliReconstructedV0::kFake == v0->mcOriginType))
1038     {
1039       histPassedCutsOrigin->Fill(AliReconstructedV0::kFakeAntiLambda,i);
1040     }
1041     else if(v0->isALamCenter[i]){
1042       histPassedCutsOrigin->Fill(v0->mcOriginType,i);
1043     }
1044   }
1045 }
1046
1047 //________________________________________________________________________
1048 AliReconstructedV0::MCV0Origin_t AliAnalysisV0Lam::DetermineV0Origin(AliAODv0 *v0, TClonesArray *mcArray)
1049 {
1050   // Determines the particle type (identity of it or of its parent particle)
1051   // from MC truth information
1052   AliReconstructedV0::MCV0Origin_t mcV0Origin = AliReconstructedV0::kUnassigned;
1053   //Get the MCParticle index for the V0
1054   int v0Id = GetV0MCParticleID(v0,mcArray);
1055   if (v0Id > 0){ // A real MC particle exists for this V0
1056     AliAODMCParticle* mcV0 = (AliAODMCParticle*)mcArray->At(v0Id);
1057     // Get the PDG code for this particle (or get the PDG code of its
1058     // mother in the case of secondary particles).
1059     // The PDG code gets converted into an MCV0Origin_t object
1060     mcV0Origin = DeterminePdgCodeOfMcParticle(mcV0,mcArray);
1061   }
1062   else{ // No MC truth exists for this V0.  It is fake.
1063     mcV0Origin = AliReconstructedV0::kFake;
1064   }
1065   return mcV0Origin;
1066 }
1067
1068 //________________________________________________________________________
1069 void AliAnalysisV0Lam::GetMCParticleMomentumTruth(double *v0Momentum, AliAODv0 *v0, TClonesArray *mcArray)
1070 {
1071   // Get the MC truth of the 3D momentum of a V0.  That info is copied into
1072   // v0Momentum
1073   int v0Id = GetV0MCParticleID(v0,mcArray);
1074   if (v0Id > 0){
1075     AliAODMCParticle* mcV0 = (AliAODMCParticle*)mcArray->At(v0Id);
1076     if(!mcV0->PxPyPz(v0Momentum)) cout<<"Err copying momentum truth"<<endl;
1077   }
1078 }
1079
1080 //________________________________________________________________________
1081 int AliAnalysisV0Lam::GetV0MCParticleID(AliAODv0 *v0, TClonesArray *mcArray)
1082 {
1083   // Returns the MCParticle index of a V0. Do this by finding the
1084   // corresponding MCParticles of the daughter tracks. If those daughter
1085   // MCParticles have the same mother, return the MCParticle index of that
1086   // mother.  If they don't have the same mother (or if both tracks are
1087   // primary) the V0 is a fake.  In that case, return -1.
1088   AliAODTrack* daughterTrackPos = (AliAODTrack*)v0->GetDaughter(0);
1089   AliAODTrack* daughterTrackNeg = (AliAODTrack*)v0->GetDaughter(1);
1090   daughterTrackPos->SetAODEvent(fAOD);
1091   daughterTrackNeg->SetAODEvent(fAOD);
1092   AliAODMCParticle* mcParticlePos = (AliAODMCParticle*)mcArray->At(abs(daughterTrackPos->GetLabel()));
1093   AliAODMCParticle* mcParticleNeg = (AliAODMCParticle*)mcArray->At(abs(daughterTrackNeg->GetLabel()));
1094   if(!(mcParticlePos) || !(mcParticleNeg)){
1095     //if either of these does not exist, V0 is fake.
1096     return -1;
1097   }
1098   //mcparticle->GetMother() will return a "-1" if the particle doesn't have a true mother (i.e. it's a fake track or primary)
1099   int motherOfPosID = mcParticlePos->GetMother();
1100   int motherOfNegID = mcParticleNeg->GetMother();
1101   if ((motherOfPosID > 0) && (motherOfPosID == motherOfNegID)){
1102     // Both daughter tracks refer to the same mother.  Return the MCParticle
1103     // index of that mother.
1104     return motherOfPosID;
1105   }
1106   else return -1; //Mother does not exist, or they refer to different
1107   // mothers. So this V0 is a fake
1108 }
1109
1110 //________________________________________________________________________
1111 AliReconstructedV0::MCV0Origin_t AliAnalysisV0Lam::DeterminePdgCodeOfMcParticle(AliAODMCParticle *mcParticle, TClonesArray *mcArray)
1112 {
1113   // Get the PDG code for this particle (or get the PDG code of its
1114   // mother in the case of secondary particles)
1115   // The PDG code gets converted into an MCV0Origin_t object
1116   AliReconstructedV0::MCV0Origin_t mcParticleOrigin = AliReconstructedV0::kUnassigned;
1117   int v0PDG = mcParticle->GetPdgCode();
1118   //find if it has a parent and note the parent's pdg code
1119   int motherOfV0ID = mcParticle->GetMother();
1120
1121   if(3122 == v0PDG){ //V0 is a Lambda
1122     if (motherOfV0ID <= 0) mcParticleOrigin = AliReconstructedV0::kPrimaryLambda;
1123     else { //V0 has a mother
1124       AliAODMCParticle* mcMotherOfV0 = (AliAODMCParticle*)mcArray->At(motherOfV0ID);
1125       int motherOfV0PDG = mcMotherOfV0->GetPdgCode();
1126       if(3212 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimarySigmaZero;
1127       else if(3322 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimaryCascadeZero;
1128       else if(3312 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimaryCascadeMinus;
1129       else if(3334 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimaryOmega;
1130       else if(3224 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kExcitedSigma;
1131       else if(3214 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kExcitedSigma;
1132       else if(3114 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kExcitedSigma;
1133       else {
1134         mcParticleOrigin = AliReconstructedV0::kOtherOriginLambda;
1135       }
1136     }
1137   }
1138   else if(-3122 == v0PDG){ // V0 is an Antilambda
1139     if (motherOfV0ID <= 0) mcParticleOrigin = AliReconstructedV0::kPrimaryAntiLambda;
1140     else { // V0 has a mother
1141       AliAODMCParticle* mcMotherOfV0 = (AliAODMCParticle*)mcArray->At(motherOfV0ID);
1142       int motherOfV0PDG = mcMotherOfV0->GetPdgCode();
1143       if(-3212 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimaryAntiSigmaZero;
1144       else if(-3322 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimaryAntiCascadeZero;
1145       else if(-3312 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimaryAntiCascadePlus;
1146       else if(-3334 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kPrimaryAntiOmega;
1147       else if(-3224 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kExcitedAntiSigma;
1148       else if(-3214 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kExcitedAntiSigma;
1149       else if(-3114 == motherOfV0PDG) mcParticleOrigin = AliReconstructedV0::kExcitedAntiSigma;
1150       else {
1151         mcParticleOrigin = AliReconstructedV0::kOtherOriginAntiLambda;
1152       }
1153     }
1154   }
1155   else if(310 == mcParticle->GetPdgCode()) mcParticleOrigin = AliReconstructedV0::kKZeroShort;
1156   return mcParticleOrigin;
1157 }
1158
1159 //________________________________________________________________________
1160 void AliAnalysisV0Lam::BinOriginInformationForMCParticles(TH1F *mcOriginalV0Hist, TH1F *mcV0FinderHist, TH2F *mcV0PassedCutsHist)
1161 {
1162   //Bin number of particles of each V0 type (e.g. primary lambda, lambda from Xi decay, etc.)
1163   //Also bin fraction of particles of each type still remaining at different stages
1164   SetBinsOnOriginHists(mcOriginalV0Hist); //Put labels on these histograms
1165   SetBinsOnOriginHists(mcV0FinderHist);
1166   SetBinsOnOriginHists(mcV0PassedCutsHist);
1167   TH1F *originalToV0Ratio = (TH1F*)mcV0FinderHist->Clone("originalToV0Ratio");
1168   originalToV0Ratio->Divide(mcOriginalV0Hist);
1169   for(int i = 0; i < AliReconstructedV0::kOriginTypeMax+1; i++){
1170     if(mcOriginalV0Hist->GetBinContent(i+1) >= 1){
1171       //Bin fraction of particles remaining at the V0 Finder stage
1172       //Only fill bin if original V0 distribution had content in that bin.
1173       //This avoids divide by zero problems
1174       fRemainingFromBeginningToV0Finder->Fill(i,originalToV0Ratio->GetBinContent(i+1));
1175     }
1176   }
1177   //Add yield results into output histograms
1178   fMCTruthOfOriginalParticles->Add(mcOriginalV0Hist);
1179   fMCTruthOfV0FinderParticles->Add(mcV0FinderHist);
1180   fMCTruthOfReconstructedParticles->Add(mcV0PassedCutsHist);
1181   delete originalToV0Ratio;
1182   for(int i = 0; i < fNumberOfVariableCutValues; i++){
1183     //Need to use a loop here because different variable cuts lead
1184     //to different distributions of reconstructed particles
1185     TH1F *originalToReconstructedRatio = (TH1F*)mcV0PassedCutsHist->ProjectionX("originalToReconstructedRatio",i+1,i+1);
1186     TH1F *v0FinderToReconstructedRatio = (TH1F*)mcV0PassedCutsHist->ProjectionX("v0FinderToReconstructedRatio",i+1,i+1);
1187     originalToReconstructedRatio->Divide(mcOriginalV0Hist);
1188     v0FinderToReconstructedRatio->Divide(mcV0FinderHist);
1189     for(int j = 0; j < AliReconstructedV0::kOriginTypeMax+1; j++){
1190       // Bin fraction of particles remaining after all reconstruction cuts
1191       // have been applied.  Only report a fraction remaining as zero if
1192       // the original event had one or more particles of that type
1193       if(mcOriginalV0Hist->GetBinContent(j+1) >= 1){
1194         fRemainingFromBeginningToRecon->Fill(i,j,originalToReconstructedRatio->GetBinContent(j+1));
1195       }
1196       if(mcV0FinderHist->GetBinContent(j+1) >= 1){
1197         fRemainingFromV0FinderToRecon->Fill(i,j,v0FinderToReconstructedRatio->GetBinContent(j+1));
1198       }
1199     }
1200     delete originalToReconstructedRatio;
1201     delete v0FinderToReconstructedRatio;
1202   }
1203 }
1204
1205 //________________________________________________________________________
1206 void AliAnalysisV0Lam::SetBinsOnOriginHists(TH1 *mcHist)
1207 {
1208   mcHist->GetXaxis()->SetBinLabel(1,"OtherV0");
1209   mcHist->GetXaxis()->SetBinLabel(2,"Fake");
1210   mcHist->GetXaxis()->SetBinLabel(3,"Fake #Lambda");
1211   mcHist->GetXaxis()->SetBinLabel(4,"#Lambda");
1212   mcHist->GetXaxis()->SetBinLabel(5,"#Sigma0");
1213   mcHist->GetXaxis()->SetBinLabel(6,"#Sigma*");
1214   mcHist->GetXaxis()->SetBinLabel(7,"#Xi0");
1215   mcHist->GetXaxis()->SetBinLabel(8,"#Xi-");
1216   mcHist->GetXaxis()->SetBinLabel(9,"#Omega");
1217   mcHist->GetXaxis()->SetBinLabel(10,"Other #Lambda");
1218   mcHist->GetXaxis()->SetBinLabel(11,"Fake #bar{#Lambda}");
1219   mcHist->GetXaxis()->SetBinLabel(12,"#bar{#Lambda}");
1220   mcHist->GetXaxis()->SetBinLabel(13,"#bar{#Sigma}0");
1221   mcHist->GetXaxis()->SetBinLabel(14,"#bar{#Sigma}*");
1222   mcHist->GetXaxis()->SetBinLabel(15,"#bar{#Xi}0");
1223   mcHist->GetXaxis()->SetBinLabel(16,"#bar{#Xi}+");
1224   mcHist->GetXaxis()->SetBinLabel(17,"#bar{#Omega}");
1225   mcHist->GetXaxis()->SetBinLabel(18,"Other #bar{#Lambda}");
1226   mcHist->GetXaxis()->SetBinLabel(19,"K0s");
1227   return;
1228 }
1229
1230 //________________________________________________________________________
1231 void AliAnalysisV0Lam::SetBinsOnOriginHists(TH3 *mcHist)
1232 {
1233   mcHist->GetYaxis()->SetBinLabel(1,"OtherV0");
1234   mcHist->GetYaxis()->SetBinLabel(2,"Fake");
1235   mcHist->GetYaxis()->SetBinLabel(3,"Fake #Lambda");
1236   mcHist->GetYaxis()->SetBinLabel(4,"#Lambda");
1237   mcHist->GetYaxis()->SetBinLabel(5,"#Sigma0");
1238   mcHist->GetYaxis()->SetBinLabel(6,"#Sigma*");
1239   mcHist->GetYaxis()->SetBinLabel(7,"#Xi0");
1240   mcHist->GetYaxis()->SetBinLabel(8,"#Xi-");
1241   mcHist->GetYaxis()->SetBinLabel(9,"#Omega");
1242   mcHist->GetYaxis()->SetBinLabel(10,"Other #Lambda");
1243   mcHist->GetYaxis()->SetBinLabel(11,"Fake #bar{#Lambda}");
1244   mcHist->GetYaxis()->SetBinLabel(12,"#bar{#Lambda}");
1245   mcHist->GetYaxis()->SetBinLabel(13,"#bar{#Sigma}0");
1246   mcHist->GetYaxis()->SetBinLabel(14,"#bar{#Sigma}*");
1247   mcHist->GetYaxis()->SetBinLabel(15,"#bar{#Xi}0");
1248   mcHist->GetYaxis()->SetBinLabel(16,"#bar{#Xi}+");
1249   mcHist->GetYaxis()->SetBinLabel(17,"#bar{#Omega}");
1250   mcHist->GetYaxis()->SetBinLabel(18,"Other #bar{#Lambda}");
1251   mcHist->GetYaxis()->SetBinLabel(19,"K0s");
1252   return;
1253 }
1254
1255 //________________________________________________________________________
1256 bool AliAnalysisV0Lam::IsInjectedParticle(AliAODv0 *v0, TClonesArray *mcArray, Int_t numberOfLastHijingLabel)
1257 {
1258   //Check if a Monte Carlo particle comes from the base MC event, or if it
1259   //was injected.  Primary particles are injected if they have an
1260   //AliAODMCParticle::GetLabel() greater than AliGenHijingEventHeader::NProduced()-1
1261   bool isInjected = false;
1262   Int_t v0ID = GetV0MCParticleID(v0,mcArray);
1263   if(v0ID > -1){ //if the v0 comes from an actual MC particle V0
1264     AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(v0ID);
1265     if(mcParticle->GetMother() > -1){ //if it has a mother
1266       AliAODMCParticle *mcMother = (AliAODMCParticle*)mcArray->At(mcParticle->GetMother());
1267       if(!mcMother) return true; // if this doesn't exist, there was an error
1268       if(mcMother->GetLabel() > numberOfLastHijingLabel) isInjected = true;
1269     }
1270     else if(mcParticle->GetLabel() > numberOfLastHijingLabel) isInjected = true;
1271   }
1272   return isInjected;
1273 }
1274
1275 //________________________________________________________________________
1276 bool AliAnalysisV0Lam::IsCorrectEventTrigger()
1277 {
1278   //Pick out Central, SemiCentral, and MinBias events.  False if not using one of those event triggers.
1279   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral));
1280   return isSelected;
1281 }
1282
1283 //________________________________________________________________________
1284 void AliAnalysisV0Lam::AddV0ToMultiplicityCounts(AliReconstructedV0 *v0, vector<int> & lambdaCount, vector<int> & antiLambdaCount)
1285 {
1286   //If a V0 is a Lambda or an Antilambda, add one to the respective V0
1287   //yields for this event.  Depending on the variable cut value, the V0
1288   //may or may not get categorized as a (anti)Lambda.
1289   //This information is used for histogramming event multiplicities.
1290   for(int i = 0; i < fNumberOfVariableCutValues; i++){
1291     if(v0->isLamCenter[i]) lambdaCount[i]++;
1292     if(v0->isALamCenter[i]) antiLambdaCount[i]++;
1293   }
1294 }
1295
1296 //________________________________________________________________________
1297 void AliAnalysisV0Lam::HistogramEventMultiplicities(vector<int> & lambdaCount, vector<int> & antiLambdaCount, int centralityBin)
1298 {
1299   //Add the event yields to the output yield histograms
1300   for(int i = 0; i < fNumberOfVariableCutValues; i++){
1301     //Centrality integrated histograms
1302     fMultDistLambda->Fill(i, lambdaCount[i]);
1303     fMultDistAntiLambda->Fill(i, antiLambdaCount[i]);
1304     //Centrality differential histograms
1305     fMultCentLambda->Fill(i, centralityBin+1, lambdaCount[i]);
1306     fMultCentAntiLambda->Fill(i, centralityBin+1, antiLambdaCount[i]);
1307   }
1308 }
1309
1310 //________________________________________________________________________
1311 void AliAnalysisV0Lam::FillTPCSignalHists(const AliReconstructedV0 *v0, const double posDaughterP, const double posDaughterTPCSignal, const double negDaughterP, const double negDaughterTPCSignal)
1312 {
1313   //Histogram P vs TPCsignal for V0 daughters
1314   if(v0->isLamCenter[fDefaultVariableCutIndex])
1315   {
1316     fTPCVsPPosLam->Fill(posDaughterP,posDaughterTPCSignal);
1317     fTPCVsPNegLam->Fill(negDaughterP,negDaughterTPCSignal);
1318   }
1319   if(v0->isALamCenter[fDefaultVariableCutIndex])
1320   {
1321     fTPCVsPPosALam->Fill(posDaughterP,posDaughterTPCSignal);
1322     fTPCVsPNegALam->Fill(negDaughterP,negDaughterTPCSignal);
1323   }
1324 }
1325
1326 //________________________________________________________________________
1327 void AliAnalysisV0Lam::CheckForFakeV0s(const AliReconstructedV0 *v0, TH1F *mcFakeParticleIdentity, TH1F *mcOtherV0Identity, const AliReconstructedV0::MCV0Origin_t mcV0Origin)
1328 {
1329   // Used in MC studies to determine how many reconstructed Lambdas and
1330   // Antilambdas are actually fake.  For simplicity, this is only done for
1331   // the default value of the variable reconstruction cut.
1332   if(v0->isLamCenter[fDefaultVariableCutIndex]
1333      || v0->isALamCenter[fDefaultVariableCutIndex])
1334   {
1335     if(AliReconstructedV0::kFake == mcV0Origin){
1336       //These V0s are fake
1337       if(v0->isLamCenter[fDefaultVariableCutIndex]) mcFakeParticleIdentity->Fill(0);
1338       if(v0->isALamCenter[fDefaultVariableCutIndex]) mcFakeParticleIdentity->Fill(1);
1339     }
1340     if(AliReconstructedV0::kUnassigned == mcV0Origin){
1341       // These V0s correspond to actual MC particles, but they aren't
1342       // Lambdas, Antilambdas, or K0s.
1343       if(v0->isLamCenter[fDefaultVariableCutIndex]) mcOtherV0Identity->Fill(0);
1344       if(v0->isALamCenter[fDefaultVariableCutIndex]) mcOtherV0Identity->Fill(1);
1345     }
1346   }
1347 }
1348
1349 //________________________________________________________________________
1350 double AliAnalysisV0Lam::CalculateKstar(double momentum1[3], double momentum2[3], double mass1, double mass2)
1351 {
1352   // Calculate k* for any pair of particles, regardless of whether the
1353   // particles have the same mass.
1354   double kstar = 0.;
1355   double e1 = 0.;
1356   double e2 = 0.;
1357   for(int i = 0; i < 3; i++){
1358     kstar -= pow(momentum1[i]-momentum2[i],2);
1359     e1 += pow(momentum1[i],2);
1360     e2 += pow(momentum2[i],2);
1361   }
1362   e1 += pow(mass1,2);
1363   e1 = sqrt(e1);
1364   e2 += pow(mass2,2);
1365   e2 = sqrt(e2);
1366
1367   kstar += pow(e1-e2,2);
1368
1369   double totalMomentumSquared = 0;
1370   for(int i = 0; i < 3; i++){
1371     totalMomentumSquared -= pow(momentum1[i]+momentum2[i],2);
1372   }
1373   totalMomentumSquared += pow(e1+e2,2);
1374   kstar -= pow((pow(mass1,2)-pow(mass2,2)),2)/totalMomentumSquared;
1375
1376   kstar *= -1.;
1377   kstar = sqrt(kstar); //At this point, we've actually calculated Qinv
1378   kstar *= 0.5; // kstar is 0.5*Qinv
1379   return kstar;
1380 }
1381
1382 //________________________________________________________________________
1383 void AliAnalysisV0Lam::BinMomentumSmearing(double *v01MomentumRecon, double *v01MomentumTruth, double *v02MomentumRecon, double *v02MomentumTruth, int pairType)
1384 {
1385   // Used in MC studies to determine the amount of pair-wise relative-
1386   // momentum smearing.  For a given pair of V0s, fills output histograms
1387   // with reconstructed k* vs true k*.  Separate histograms are filled
1388   // for each pair type
1389   double v01PMag = sqrt( pow(v01MomentumTruth[0],2) + pow(v01MomentumTruth[1],2) + pow(v01MomentumTruth[2],2) );
1390   if(v01PMag < 0.0001) return;  //Not a real V0
1391   double v02PMag = sqrt( pow(v02MomentumTruth[0],2) + pow(v02MomentumTruth[1],2) + pow(v02MomentumTruth[2],2) );
1392   if(v02PMag < 0.0001) return;  //Not a real V0
1393   double reconKstar = CalculateKstar(v01MomentumRecon,v02MomentumRecon, fPDGLambda, fPDGLambda);
1394   double truthKstar = CalculateKstar(v01MomentumTruth,v02MomentumTruth, fPDGLambda, fPDGLambda);
1395   if(0 == pairType){ //Lambda-Lambda
1396     fHistPsmearingKreconVsKtruthLL->Fill(reconKstar,truthKstar);
1397     fHistPsmearingKreconMinusKtruthLL->Fill(reconKstar-truthKstar);
1398   }
1399   if(1 == pairType){ //Antilambda-Antilambda
1400     fHistPsmearingKreconVsKtruthAA->Fill(reconKstar,truthKstar);
1401     fHistPsmearingKreconMinusKtruthAA->Fill(reconKstar-truthKstar);
1402   }
1403   if(2 == pairType){ //Lambda-Antilambda
1404     fHistPsmearingKreconVsKtruthLA->Fill(reconKstar,truthKstar);
1405     fHistPsmearingKreconMinusKtruthLA->Fill(reconKstar-truthKstar);
1406   }
1407 }
1408
1409 //________________________________________________________________________
1410 void AliAnalysisV0Lam::DoPairStudies(const AliAnalysisV0LamEvent *event, const int centralityBin)
1411 {
1412   //Loop over same- and mixed-event pairs and bin correlation function
1413   //numerators and denominators.
1414
1415   //Default cut values used for the average separation cuts
1416   double avgSepIdenticalProtonCut = 3.; // used for same-charge prot avg sep
1417   double avgSepIdenticalPionCut = 4.; // used for same-charge pi-pi
1418   double avgSepNonIdenticalCut = 3.5; //used for same-charge pion-proton avg sep
1419   double &variableAvgSepValue = avgSepIdenticalPionCut; //the identical
1420   //pion cut will be the variable cut in this run.  This needs to be changed
1421   //by hand if one wants to study a different cut
1422
1423   // These values are used for the average separation cuts of the
1424   // daughter tracks. When the average separation cut is being studied,
1425   // we can loop over many different cut values to see how those different
1426   // cuts affect the correlation function results
1427   double avgSepCutArray[12] = {0.,1.,2.,2.5,3.,3.5,4.,4.5,5.,6.,7.,10.};
1428
1429   int numberOfAvgSepCutIndices; //this is the size of the avgSepCutArray
1430   if(fIsUsingVariableAvgSepCut) numberOfAvgSepCutIndices = sizeof(avgSepCutArray) / sizeof (avgSepCutArray[0]);
1431   else numberOfAvgSepCutIndices = 1;
1432
1433   int defaultVariableAvgSepCutIndex = 0;
1434   if(fIsUsingVariableAvgSepCut) {
1435     for(int i = 0; i < numberOfAvgSepCutIndices; i++){
1436       // If variable avg sep cuts are being used, this find the default
1437       // value for the avg sep cut currently being studied.
1438       if (fabs(variableAvgSepValue - avgSepCutArray[i]) < 0.1) defaultVariableAvgSepCutIndex = i;
1439     }
1440   }
1441   for(int cutIndex = 0; cutIndex < fNumberOfVariableCutValues; cutIndex++)
1442   { // Start looping over all variable cut values
1443     for(int i=0; i < event->fNumberCandidateV0; i++)
1444     { //Start looping over reconstructed V0s in this event
1445       bool center1Lam  = event->fReconstructedV0[i].isLamCenter[cutIndex];
1446       bool center1ALam = event->fReconstructedV0[i].isALamCenter[cutIndex];
1447       // Disregard V0 if it wasn't reconstructed as a center (anti)Lambda
1448       if(!(center1Lam || center1ALam)) continue;
1449       // Disregard V0 if it was removed via the judgment cuts
1450       if(event->fReconstructedV0[i].isDeemedUnworthy[cutIndex]) continue;
1451       for(int eventNumber=0; eventNumber<nEventsToMix+1; eventNumber++)
1452       { // Event buffer loop: eventNumber=0 is the current event, all other eventNumbers are past events
1453         int startBin=0;
1454         // For same event pairs, start 2nd V0 loop at i+1 V0 to avoid
1455         // double counting
1456         if(eventNumber==0) startBin=i+1;
1457         for(int j=startBin; j<(event+eventNumber)->fNumberCandidateV0; j++)
1458         { // Second V0 loop (from past or current event)
1459           if(eventNumber==0)
1460           { // Don't make pairs of V0s if they shared daughter tracks.
1461             // This is redundant if the judgment cut is already employed
1462             if(event->fReconstructedV0[i].daughter1ID
1463                == (event+eventNumber)->fReconstructedV0[j].daughter1ID) continue;
1464             if(event->fReconstructedV0[i].daughter1ID
1465                == (event+eventNumber)->fReconstructedV0[j].daughter2ID) continue;
1466             if(event->fReconstructedV0[i].daughter2ID
1467                == (event+eventNumber)->fReconstructedV0[j].daughter1ID) continue;
1468             if(event->fReconstructedV0[i].daughter2ID
1469                == (event+eventNumber)->fReconstructedV0[j].daughter2ID) continue;
1470           }
1471           //Disregard second V0 if it was removed via judgment cuts
1472           if((event+eventNumber)->fReconstructedV0[j].isDeemedUnworthy[cutIndex]) continue;
1473           // A central V0 has a mass that falls within the accepted inv
1474           // mass range.  Only make pairs with central V0s
1475           bool center2Lam  = (event+eventNumber)->fReconstructedV0[j].isLamCenter[cutIndex];
1476           bool center2ALam = (event+eventNumber)->fReconstructedV0[j].isALamCenter[cutIndex];
1477           if(!(center2Lam || center2ALam)) continue;
1478           // Now we calculate a bunch of values that are used later during
1479           // histogramming.
1480           double pairKt = pow(event->fReconstructedV0[i].v0Momentum[0] + (event+eventNumber)->fReconstructedV0[j].v0Momentum[0],2.);
1481           pairKt+= pow(event->fReconstructedV0[i].v0Momentum[1] + (event+eventNumber)->fReconstructedV0[j].v0Momentum[1],2.);
1482           pairKt = sqrt(pairKt)/2.;
1483           //Calculate k* for V0s and daughters using different mass assumptions
1484           double pairKstarLam = CalculateKstar(event->fReconstructedV0[i].v0Momentum, (event+eventNumber)->fReconstructedV0[j].v0Momentum, fPDGLambda,fPDGLambda);
1485           double pairKstarProtPlus = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGProton,fPDGProton);
1486           double pairKstarProtMinus = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGProton,fPDGProton);
1487           double pairKstarPiPlus = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGPion,fPDGPion);
1488           double pairKstarPiMinus = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGPion,fPDGPion);
1489           //used for lambda-antilambda daughter kstar
1490           double pairKstarProtPlusPiPlus = 0;
1491           double pairKstarProtMinusPiMinus = 0;
1492           double pairKstarProtPlusProtMinus = 0;
1493           double pairKstarPiPlusPiMinus = 0;
1494
1495           // Need to be careful when calculating the k* of daughter tracks
1496           // of non-identical V0s
1497           if(center1Lam && center2ALam){
1498             pairKstarProtPlusPiPlus = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGProton,fPDGPion);
1499             pairKstarProtMinusPiMinus = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGPion,fPDGProton);
1500             pairKstarProtPlusProtMinus = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGProton,fPDGProton);
1501             pairKstarPiPlusPiMinus = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGPion,fPDGPion);
1502           }
1503           else if(center1ALam && center2Lam){
1504             pairKstarProtPlusPiPlus = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGPion,fPDGProton);
1505             pairKstarProtMinusPiMinus = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGProton,fPDGPion);
1506             pairKstarProtPlusProtMinus = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGProton,fPDGProton);
1507             pairKstarPiPlusPiMinus = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGPion,fPDGPion);
1508           }
1509           double pairKstarProtPlusPiMinus1 = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGProton,fPDGPion); // only relevant for LamLam
1510           double pairKstarProtPlusPiMinus2 = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGPion,fPDGProton); // only relevant for LamLam
1511
1512           double pairKstarProtMinusPiPlus1 = CalculateKstar(event->fReconstructedV0[i].daughterNegMomentum,(event+eventNumber)->fReconstructedV0[j].daughterPosMomentum, fPDGProton,fPDGPion); // only relevant for ALamALam
1513           double pairKstarProtMinusPiPlus2 = CalculateKstar(event->fReconstructedV0[i].daughterPosMomentum,(event+eventNumber)->fReconstructedV0[j].daughterNegMomentum, fPDGPion, fPDGProton); // only relevant for ALamALam
1514           //Now find the average separation distance between daughter pairs.  Used to make a merging/splitting cut.
1515           double avgSepPos = GetAverageSeparation(event->fReconstructedV0[i].daughterPosGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterPosGlobalPositions);
1516           double avgSepNeg = GetAverageSeparation(event->fReconstructedV0[i].daughterNegGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterNegGlobalPositions);
1517           double avgSepNegPos = GetAverageSeparation(event->fReconstructedV0[i].daughterNegGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterPosGlobalPositions);
1518           double avgSepPosNeg = GetAverageSeparation(event->fReconstructedV0[i].daughterPosGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterNegGlobalPositions);
1519           double correctedAvgSepPos = GetAverageSeparation(event->fReconstructedV0[i].daughterPosCorrectedGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterPosCorrectedGlobalPositions);
1520           double correctedAvgSepNeg = GetAverageSeparation(event->fReconstructedV0[i].daughterNegCorrectedGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterNegCorrectedGlobalPositions);
1521           double correctedAvgSepNegPos = GetAverageSeparation(event->fReconstructedV0[i].daughterNegCorrectedGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterPosCorrectedGlobalPositions);
1522           double correctedAvgSepPosNeg = GetAverageSeparation(event->fReconstructedV0[i].daughterPosCorrectedGlobalPositions, (event+eventNumber)->fReconstructedV0[j].daughterNegCorrectedGlobalPositions);
1523
1524           //Now we get to the actual pair histogramming
1525           if(eventNumber==0) //Same event pair histogramming
1526           {
1527             //We do separate binning for each pair type
1528             if(center1Lam && center2Lam){
1529               // Some histograms we only fill when default variable cut
1530               // values have been used
1531               if(cutIndex == fDefaultVariableCutIndex){
1532                 //same sign tracks
1533                 fSignalLamLamProtSep->Fill(avgSepPos, pairKstarProtPlus);
1534                 fSignalLamLamPiMinusSep->Fill(avgSepNeg, pairKstarPiMinus);
1535                 fSignalLamLamProtSepCorrected->Fill(correctedAvgSepPos, pairKstarProtPlus);
1536                 fSignalLamLamPiMinusSepCorrected->Fill(correctedAvgSepNeg, pairKstarPiMinus);
1537                 //opposite sign tracks
1538                 fSignalLamLamPlusMinusSep->Fill(avgSepPosNeg,pairKstarProtPlusPiMinus1);
1539                 fSignalLamLamPlusMinusSep->Fill(avgSepNegPos,pairKstarProtPlusPiMinus2);
1540                 fSignalLamLamPlusMinusSepCorrected->Fill(correctedAvgSepPosNeg,pairKstarProtPlusPiMinus1);
1541                 fSignalLamLamPlusMinusSepCorrected->Fill(correctedAvgSepNegPos,pairKstarProtPlusPiMinus2);
1542               }
1543               for(int sepCutIndex = 0; sepCutIndex < numberOfAvgSepCutIndices; sepCutIndex++){ //looping over different avg sep cut values (if applicable)
1544                 if(fIsUsingVariableAvgSepCut) variableAvgSepValue = avgSepCutArray[sepCutIndex];
1545                 if((avgSepIdenticalProtonCut < correctedAvgSepPos)
1546                    && (avgSepIdenticalPionCut < correctedAvgSepNeg))
1547                 {
1548                   fSignalLamLam->Fill(sepCutIndex, centralityBin+1, pairKstarLam);
1549                   if(defaultVariableAvgSepCutIndex == sepCutIndex){
1550                     //This implementation doesn't work properly
1551                     fKtLamLamSig->Fill(centralityBin+1,pairKt,pairKstarLam);
1552                   }
1553                 }
1554               }
1555             }
1556             if(center1ALam && center2ALam){
1557               if(cutIndex == fDefaultVariableCutIndex){
1558                 //same sign tracks
1559                 fSignalALamALamAntiProtSep->Fill(avgSepNeg, pairKstarProtMinus);
1560                 fSignalALamALamPiPlusSep->Fill(avgSepPos, pairKstarPiPlus);
1561                 fSignalALamALamAntiProtSepCorrected->Fill(correctedAvgSepNeg, pairKstarProtMinus);
1562                 fSignalALamALamPiPlusSepCorrected->Fill(correctedAvgSepPos, pairKstarPiPlus);
1563
1564                 //opposite sign tracks
1565                 fSignalALamALamPlusMinusSep->Fill(avgSepPosNeg,pairKstarProtMinusPiPlus1);
1566                 fSignalALamALamPlusMinusSep->Fill(avgSepNegPos,pairKstarProtMinusPiPlus2);
1567                 fSignalALamALamPlusMinusSepCorrected->Fill(correctedAvgSepPosNeg,pairKstarProtMinusPiPlus1);
1568                 fSignalALamALamPlusMinusSepCorrected->Fill(correctedAvgSepNegPos,pairKstarProtMinusPiPlus2);
1569               }
1570               for(int sepCutIndex = 0; sepCutIndex < numberOfAvgSepCutIndices; sepCutIndex++){
1571                 if(fIsUsingVariableAvgSepCut) variableAvgSepValue = avgSepCutArray[sepCutIndex];
1572                 if((avgSepIdenticalPionCut < correctedAvgSepPos)
1573                    && (avgSepIdenticalProtonCut < correctedAvgSepNeg))
1574                 {
1575                   fSignalALamALam->Fill(sepCutIndex, centralityBin+1, pairKstarLam);
1576                   if(defaultVariableAvgSepCutIndex == sepCutIndex){
1577                     //This implementation doesn't work properly
1578                     fKtALamALamSig->Fill(centralityBin+1,pairKt,pairKstarLam);
1579                   }
1580                 }
1581               }
1582             }
1583             if((center1Lam && center2ALam) || (center1ALam && center2Lam)){
1584               if(cutIndex == fDefaultVariableCutIndex){
1585                 fSignalLamALamProtPiPlusSep->Fill(avgSepPos, pairKstarProtPlusPiPlus);
1586                 fSignalLamALamAntiProtPiMinusSep->Fill(avgSepNeg, pairKstarProtMinusPiMinus);
1587                 fSignalLamALamProtPiPlusSepCorrected->Fill(correctedAvgSepPos, pairKstarProtPlusPiPlus);
1588                 fSignalLamALamAntiProtPiMinusSepCorrected->Fill(correctedAvgSepNeg, pairKstarProtMinusPiMinus);
1589                 //opposite charge tracks
1590                 if(center1Lam)
1591                 {
1592                   fSignalLamALamProtSep->Fill(avgSepPosNeg, pairKstarProtPlusProtMinus);
1593                   fSignalLamALamPionSep->Fill(avgSepNegPos, pairKstarPiPlusPiMinus);
1594                   fSignalLamALamProtSepCorrected->Fill(correctedAvgSepPosNeg, pairKstarProtPlusProtMinus);
1595                   fSignalLamALamPionSepCorrected->Fill(correctedAvgSepNegPos, pairKstarPiPlusPiMinus);
1596                 }
1597                 else
1598                 {
1599                   fSignalLamALamProtSep->Fill(avgSepNegPos, pairKstarProtPlusProtMinus);
1600                   fSignalLamALamPionSep->Fill(avgSepPosNeg, pairKstarPiPlusPiMinus);
1601                   fSignalLamALamProtSepCorrected->Fill(correctedAvgSepNegPos, pairKstarProtPlusProtMinus);
1602                   fSignalLamALamPionSepCorrected->Fill(correctedAvgSepPosNeg, pairKstarPiPlusPiMinus);
1603                 }
1604               }
1605               for(int sepCutIndex = 0; sepCutIndex < numberOfAvgSepCutIndices; sepCutIndex++){
1606                 if(fIsUsingVariableAvgSepCut) variableAvgSepValue = avgSepCutArray[sepCutIndex];
1607                 if((avgSepNonIdenticalCut < correctedAvgSepPos)
1608                    && (avgSepNonIdenticalCut < correctedAvgSepNeg))
1609                 {
1610                   fSignalLamALam->Fill(sepCutIndex, centralityBin+1, pairKstarLam);
1611                   if(defaultVariableAvgSepCutIndex == sepCutIndex){
1612                    //This implementation doesn't work properly
1613                     fKtLamALamSig->Fill(centralityBin+1,pairKt,pairKstarLam);
1614                   }
1615                 }
1616               }
1617             }
1618           } //end same event pair histogramming
1619           else //Mixed event pair histogramming
1620           {
1621             if(center1Lam && center2Lam){
1622               if(cutIndex == fDefaultVariableCutIndex){
1623                 fBkgLamLamProtSep->Fill(avgSepPos, pairKstarProtPlus);
1624                 fBkgLamLamPiMinusSep->Fill(avgSepNeg, pairKstarPiMinus);
1625                 fBkgLamLamProtSepCorrected->Fill(correctedAvgSepPos, pairKstarProtPlus);
1626                 fBkgLamLamPiMinusSepCorrected->Fill(correctedAvgSepNeg, pairKstarPiMinus);
1627                 //opposite sign tracks
1628                 fBkgLamLamPlusMinusSep->Fill(avgSepPosNeg,pairKstarProtPlusPiMinus1);
1629                 fBkgLamLamPlusMinusSep->Fill(avgSepNegPos,pairKstarProtPlusPiMinus2);
1630                 fBkgLamLamPlusMinusSepCorrected->Fill(correctedAvgSepPosNeg,pairKstarProtPlusPiMinus1);
1631                 fBkgLamLamPlusMinusSepCorrected->Fill(correctedAvgSepNegPos,pairKstarProtPlusPiMinus2);
1632                 if(fIsMCEvent)
1633                 { //collect momentum smearing information
1634                   int pairType = 0;
1635                   BinMomentumSmearing(event->fReconstructedV0[i].v0Momentum, event->fReconstructedV0[i].v0MomentumTruth, (event+eventNumber)->fReconstructedV0[j].v0Momentum, (event+eventNumber)->fReconstructedV0[j].v0MomentumTruth,pairType);
1636                 }
1637               }
1638               for(int sepCutIndex = 0; sepCutIndex < numberOfAvgSepCutIndices; sepCutIndex++){
1639                 if(fIsUsingVariableAvgSepCut) variableAvgSepValue = avgSepCutArray[sepCutIndex];
1640                 if((avgSepIdenticalProtonCut < correctedAvgSepPos)
1641                    && (avgSepIdenticalPionCut < correctedAvgSepNeg))
1642                 {
1643                   fBkgLamLam->Fill(sepCutIndex, centralityBin+1, pairKstarLam);
1644                   if(defaultVariableAvgSepCutIndex == sepCutIndex){
1645                     //This implementation doesn't work properly
1646                     fKtLamLamBkg->Fill(centralityBin+1,pairKt,pairKstarLam);
1647                   }
1648                 }
1649               }
1650             }
1651             if(center1ALam && center2ALam){
1652               if(cutIndex == fDefaultVariableCutIndex){
1653                 fBkgALamALamAntiProtSep->Fill(avgSepNeg, pairKstarProtMinus);
1654                 fBkgALamALamPiPlusSep->Fill(avgSepPos, pairKstarPiPlus);
1655                 fBkgALamALamAntiProtSepCorrected->Fill(correctedAvgSepNeg, pairKstarProtMinus);
1656                 fBkgALamALamPiPlusSepCorrected->Fill(correctedAvgSepPos, pairKstarPiPlus);
1657                 //opposite sign tracks
1658                 fBkgALamALamPlusMinusSep->Fill(avgSepPosNeg,pairKstarProtMinusPiPlus1);
1659                 fBkgALamALamPlusMinusSep->Fill(avgSepNegPos,pairKstarProtMinusPiPlus2);
1660                 fBkgALamALamPlusMinusSepCorrected->Fill(correctedAvgSepPosNeg,pairKstarProtMinusPiPlus1);
1661                 fBkgALamALamPlusMinusSepCorrected->Fill(correctedAvgSepNegPos,pairKstarProtMinusPiPlus2);
1662                 if(fIsMCEvent)
1663                 { //collect momentum smearing information
1664                   int pairType = 1;
1665                   BinMomentumSmearing(event->fReconstructedV0[i].v0Momentum, event->fReconstructedV0[i].v0MomentumTruth, (event+eventNumber)->fReconstructedV0[j].v0Momentum, (event+eventNumber)->fReconstructedV0[j].v0MomentumTruth,pairType);
1666                 }
1667               }
1668               for(int sepCutIndex = 0; sepCutIndex < numberOfAvgSepCutIndices; sepCutIndex++){
1669                 if(fIsUsingVariableAvgSepCut) variableAvgSepValue = avgSepCutArray[sepCutIndex];
1670                 if((avgSepIdenticalPionCut < correctedAvgSepPos)
1671                    && (avgSepIdenticalProtonCut < correctedAvgSepNeg))
1672                 {
1673                   fBkgALamALam->Fill(sepCutIndex, centralityBin+1, pairKstarLam);
1674                   if(defaultVariableAvgSepCutIndex == sepCutIndex){
1675                     //This implementation doesn't work properly
1676                     fKtALamALamBkg->Fill(centralityBin+1,pairKt,pairKstarLam);
1677                   }
1678                 }
1679               }
1680             }
1681             if((center1Lam && center2ALam) || (center1ALam && center2Lam)){
1682               if(cutIndex == fDefaultVariableCutIndex){
1683                 fBkgLamALamProtPiPlusSep->Fill(avgSepPos, pairKstarProtPlusPiPlus);
1684                 fBkgLamALamAntiProtPiMinusSep->Fill(avgSepNeg, pairKstarProtMinusPiMinus);
1685                 fBkgLamALamProtPiPlusSepCorrected->Fill(correctedAvgSepPos, pairKstarProtPlusPiPlus);
1686                 fBkgLamALamAntiProtPiMinusSepCorrected->Fill(correctedAvgSepNeg, pairKstarProtMinusPiMinus);
1687                 //opposite charge tracks
1688                 if(center1Lam)
1689                 {
1690                   fBkgLamALamProtSep->Fill(avgSepPosNeg, pairKstarProtPlusProtMinus);
1691                   fBkgLamALamPionSep->Fill(avgSepNegPos, pairKstarPiPlusPiMinus);
1692                   fBkgLamALamProtSepCorrected->Fill(correctedAvgSepPosNeg, pairKstarProtPlusProtMinus);
1693                   fBkgLamALamPionSepCorrected->Fill(correctedAvgSepNegPos, pairKstarPiPlusPiMinus);
1694                 }
1695                 else
1696                 {
1697                   fBkgLamALamProtSep->Fill(avgSepNegPos, pairKstarProtPlusProtMinus);
1698                   fBkgLamALamPionSep->Fill(avgSepPosNeg, pairKstarPiPlusPiMinus);
1699                   fBkgLamALamProtSepCorrected->Fill(correctedAvgSepNegPos, pairKstarProtPlusProtMinus);
1700                   fBkgLamALamPionSepCorrected->Fill(correctedAvgSepPosNeg, pairKstarPiPlusPiMinus);
1701                 }
1702                 if(fIsMCEvent)
1703                 { //collect momentum smearing information
1704                   int pairType = 2;
1705                   BinMomentumSmearing(event->fReconstructedV0[i].v0Momentum, event->fReconstructedV0[i].v0MomentumTruth, (event+eventNumber)->fReconstructedV0[j].v0Momentum, (event+eventNumber)->fReconstructedV0[j].v0MomentumTruth,pairType);
1706                 }
1707               }
1708               for(int sepCutIndex = 0; sepCutIndex < numberOfAvgSepCutIndices; sepCutIndex++){
1709                 if(fIsUsingVariableAvgSepCut) variableAvgSepValue = avgSepCutArray[sepCutIndex];
1710                 if((avgSepNonIdenticalCut < correctedAvgSepPos)
1711                    && (avgSepNonIdenticalCut < correctedAvgSepNeg))
1712                 {
1713                   fBkgLamALam->Fill(sepCutIndex, centralityBin+1, pairKstarLam);
1714                   if(defaultVariableAvgSepCutIndex == sepCutIndex){
1715                     //This implementation doesn't work properly
1716                     fKtLamALamBkg->Fill(centralityBin+1,pairKt,pairKstarLam);
1717                   }
1718                 }
1719               }//end loop over sepCutIndex
1720             }//end Lam-ALam pair binning
1721           }//end mixed event pair histogramming
1722         }//end past event
1723       }//end event buffer
1724     }//end current event
1725   }//end variable cut loop
1726 }//end DoPairStudies()
1727
1728