PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ////////////////////////////////////////////////////////////////////////////
  17 //
  18 //  This class is used to perform femtoscopic analysis on K0s particles,
  19 //  which are reconstructed using the AliAODv0 class.
  20 //
  21 //  authors: Matthew Steinpreis (matthew.steinpreis@cern.ch)
  22 //
  23 //  Change log:
  24 //      - TOF mismatch function calls changed (4/18/13)
  25 //      - added minimum decay length cut (rarely used though) (3/28/13)
  26 //      - K0 multiplicity histogram now filled with "unskippedCount" instead
  27 //        of k0Count (which included skipped k0s with shared daughters)
  28 //        (3/25/13)
  29 //      - added hists for 3D mom. in LF and PRF (3/28/13)
  30 //      - changed calling of PIDResponse (should be same actions) (3/28/13)
  31 //      - keep "side" K0s for mass plot (4/18)
  32 //              - tweaked loading and skipping appropriately
  33 //              - use merit test to skip sides (against good and side K0s)
  34 //              - a good K0 cant be skipped by a side
  35 //      - moved TPC propagation (via Hans' method) up to v0 level, which now
  36 //        uses an AliAODTrack(AliVTrack) instead of AliESDtrack (5/31/13)
  37 //      - added primary vertex subtraction in TPC propagation   (5/31/13)
  38 //      - removed all instances of AliESDtrack usage    (5/31/13)
  39 //      - removed old separation method/histograms      (5/31/13)
  40 //      - tidied up LCMS boost                                                  (6/10/13)
  41 //      - added new boosting prescription, get q out-side-long for LCMS and PRF (6/24/13)
  42 //              - added histograms and values for LCMS momenta (for simulation)
  43 //      - added random particle order switch in correlations (9/09/13)
  44 //      - added more bins for 3D OSL analysis (9/19/13)
  45 //      - added merit cut choice, pass as argument (10/16/13)
  46 //              - 1-mass, 2-v0dca, 3-dddca, 4-combination (used to be v0dca)
  47 //      - added passable argument for two-track minimum separation (10/16/13)
  48 //      - added boolean to turn off field-sign dependence for train (10/30/13)
  49 //      - changed destructors to minimize lost memory (11/27/13)
  50 //      - added Case3D to switch off all 3D objects (11/27/13)
  51 //      - added centrality flattening routine (and switch) (12/04/13)
  52 //      - added event plane stuff (12/11/13)
  53 //      - added NPsiBins argument (1/8/14)
  54 ////////////////////////////////////////////////////////////////////////////////
  55
  56
  57
  58 #include <iostream>
  59 #include <math.h>
  60 #include "TMath.h"
  61 #include "TChain.h"
  62 #include "TFile.h"
  63 #include "TKey.h"
  64 #include "TObject.h"
  65 #include "TObjString.h"
  66 #include "TList.h"
  67 #include "TTree.h"
  68 #include "TH1F.h"
  69 #include "TH1D.h"
  70 #include "TH2D.h"
  71 #include "TH3D.h"
  72 #include "TProfile.h"
  73 #include "TCanvas.h"
  74 #include "TRandom3.h"
  75
  76 #include "AliAnalysisTask.h"
  77 #include "AliAnalysisManager.h"
  78
  79 #include "AliAODEvent.h"
  80 #include "AliAODInputHandler.h"
  81 #include "AliAODMCParticle.h"
  82 #include "AliAODv0.h"
  83 #include "AliAODRecoDecay.h"
  84 #include "AliCentrality.h"
  85
  86 #include "AliFemtoK0Analysis.h"
  87
  88 #define PI 3.1415927
  89
  90
  91 // Author: Matt Steinpreis, adapted from Dhevan Gangadharan
  92
  93 ClassImp(AliFemtoK0Analysis)
  94
  95 //________________________________________________________________________
  96 AliFemtoK0Analysis::AliFemtoK0Analysis():
  97 AliAnalysisTaskSE(),
  98   fSignDep(kFALSE),
  99   fFieldPos(kTRUE),
 100   fOnlineCase(kTRUE),
 101   fMeritCase(kTRUE),
 102   fCase3D(kFALSE),
 103   fMinDecayLength(0.0),
 104   fMeritCutChoice(0),
 105   fMinSep(0.0),
 106   fFlatCent(kFALSE),
 107   fPsiBinning(kFALSE),
 108   fNPsiBins(0),
 109   fEventCount(0),
 110   fEC(0x0),
 111   fEvt(0X0),
 112   fRandomNumber(0x0),
 113   fName(0x0),
 114   fAOD(0x0),
 115   fOutputList(0x0),
 116   fPidAOD(0x0)
 117 {
 118 }
 119 //________________________________________________________________________
 120 AliFemtoK0Analysis::AliFemtoK0Analysis(const char *name, bool SignDep, bool FieldPositive, bool OnlineCase, bool MeritCase, bool Case3D, float MinDL, int MeritCutChoice, float MinSep, bool FlatCent, bool PsiBinning, int NPsiBins)
 121 : AliAnalysisTaskSE(name),
 122   fSignDep(SignDep),
 123   fFieldPos(FieldPositive),
 124   fOnlineCase(OnlineCase),
 125   fMeritCase(MeritCase),
 126   fCase3D(Case3D),
 127   fMinDecayLength(MinDL),
 128   fMeritCutChoice(MeritCutChoice),
 129   fMinSep(MinSep),
 130   fFlatCent(FlatCent),
 131   fPsiBinning(PsiBinning),
 132   fNPsiBins(0),
 133   fEventCount(0),
 134   fEC(0x0),
 135   fEvt(0X0),
 136   fRandomNumber(0x0),
 137   fName(name),
 138   fAOD(0x0),
 139   fOutputList(0x0),
 140   fPidAOD(0x0)
 141 {
 142   //main constructor
 143   fSignDep              = SignDep;
 144   fFieldPos     = FieldPositive;
 145   fOnlineCase   = OnlineCase;
 146   fMeritCase    = MeritCase;
 147   fCase3D               = Case3D;
 148   fMinDecayLength = MinDL;
 149   fMeritCutChoice = MeritCutChoice;
 150   fMinSep               = MinSep;
 151   fFlatCent             = FlatCent;
 152   fPsiBinning   = PsiBinning;
 153   fNPsiBins             = NPsiBins;
 154
 155   // Define output slots here
 156   // Output slot #1
 157   DefineOutput(1, TList::Class());
 158
 159 }
 160 //________________________________________________________________________
 161 AliFemtoK0Analysis::AliFemtoK0Analysis(const AliFemtoK0Analysis &obj)
 162 : AliAnalysisTaskSE(obj.fName),
 163   fSignDep(obj.fSignDep),
 164   fFieldPos(obj.fFieldPos),
 165   fOnlineCase(obj.fOnlineCase),
 166   fMeritCase(obj.fMeritCase),
 167   fCase3D(obj.fCase3D),
 168   fMinDecayLength(obj.fMinDecayLength),
 169   fMeritCutChoice(obj.fMeritCutChoice),
 170   fMinSep(obj.fMinSep),
 171   fFlatCent(obj.fFlatCent),
 172   fPsiBinning(obj.fPsiBinning),
 173   fNPsiBins(obj.fNPsiBins),
 174   fEventCount(obj.fEventCount),
 175   fEC(obj.fEC),
 176   fEvt(obj.fEvt),
 177   fRandomNumber(obj.fRandomNumber),
 178   fName(obj.fName),
 179   fAOD(obj.fAOD),
 180   fOutputList(obj.fOutputList),
 181   fPidAOD(obj.fPidAOD)
 182 {
 183 }
 184 //________________________________________________________________________
 185 AliFemtoK0Analysis &AliFemtoK0Analysis::operator=(const AliFemtoK0Analysis &obj)
 186 {
 187  //Assignment operator
 188  if (this == &obj) return *this;
 189
 190  fSignDep   = obj.fSignDep;
 191  fFieldPos      = obj.fFieldPos;
 192  fOnlineCase    = obj.fOnlineCase;
 193  fMeritCase     = obj.fMeritCase;
 194  fCase3D                = obj.fCase3D;
 195  fMinDecayLength= obj.fMinDecayLength;
 196  fMeritCutChoice= obj.fMeritCutChoice;
 197  fMinSep                = obj.fMinSep;
 198  fFlatCent              = obj.fFlatCent;
 199  fPsiBinning    = obj.fPsiBinning;
 200  fNPsiBins              = obj.fNPsiBins;
 201  fEventCount    = obj.fEventCount;
 202  fEC            = obj.fEC;
 203  fEvt           = obj.fEvt;
 204  fRandomNumber  = obj.fRandomNumber;
 205  fName          = obj.fName;
 206  fAOD           = obj.fAOD;
 207  fOutputList    = obj.fOutputList;
 208  fPidAOD        = obj.fPidAOD;
 209
 210  return (*this);
 211 }
 212 //________________________________________________________________________
 213 AliFemtoK0Analysis::~AliFemtoK0Analysis()
 214 {
 215   // Destructor
 216   for(unsigned short i=0; i<kZVertexBins; i++)
 217   {
 218     for(unsigned short j=0; j<kCentBins; j++)
 219     {
 220           for(unsigned short k=0; k<fNPsiBins; k++)
 221           {
 222         fEC[i][j][k]->~AliFemtoK0EventCollection();
 223         fEC[i][j][k] = NULL;
 224       }
 225       delete [] fEC[i][j]; fEC[i][j] = NULL;
 226     }
 227     delete[] fEC[i]; fEC[i] = NULL;
 228   }
 229   delete[] fEC; fEC = NULL;
 230
 231   if(fEC){ delete fEC; fEC = NULL;}
 232   if(fRandomNumber){ delete fRandomNumber; fRandomNumber = NULL;}
 233   if(fAOD){ delete fAOD; fAOD = NULL;}
 234   if(fOutputList){ delete fOutputList; fOutputList = NULL;}
 235   if(fPidAOD){ delete fPidAOD; fPidAOD = NULL;}
 236 }
 237 //________________________________________________________________________
 238 void AliFemtoK0Analysis::MyInit()
 239 {
 240
 241   // One can set global variables here
 242   fEventCount = 0;
 243
 244   fEC = new AliFemtoK0EventCollection ***[kZVertexBins];
 245   for(unsigned short i=0; i<kZVertexBins; i++)
 246   {
 247     fEC[i] = new AliFemtoK0EventCollection **[kCentBins];
 248
 249     for(unsigned short j=0; j<kCentBins; j++)
 250     {
 251       fEC[i][j] = new AliFemtoK0EventCollection *[fNPsiBins];
 252
 253       for(unsigned short k=0; k<fNPsiBins; k++)
 254           {
 255             fEC[i][j][k] = new AliFemtoK0EventCollection(kEventsToMix+1, kMultLimit);
 256       }
 257     }
 258   }
 259
 260   AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 261   fPidAOD = aodH->GetPIDResponse();
 262
 263   fRandomNumber = new TRandom3();  //for 3D, random sign switching
 264   fRandomNumber->SetSeed(0);
 265
 266 }
 267 //________________________________________________________________________
 268 void AliFemtoK0Analysis::UserCreateOutputObjects()
 269 {
 270   // Create histograms
 271   // Called once
 272
 273   MyInit();// Initialize my settings
 274
 275   fOutputList = new TList();
 276   fOutputList->SetOwner();
 277
 278   TH1F *fHistCent = new TH1F("fHistCent","",100,0,100);
 279   fOutputList->Add(fHistCent);
 280   TH1F *fHistCentFlat = new TH1F("fHistCentFlat","",100,0,100);
 281   fOutputList->Add(fHistCentFlat);
 282   TH1F *fHistCentUsed = new TH1F("fHistCentUsed","",100,0,100);
 283   fOutputList->Add(fHistCentUsed);
 284
 285   //pion parameters
 286   TH1F* fHistDCAPiPlus = new TH1F("fHistDCAPiPlus","",100,0,10);
 287   fOutputList->Add(fHistDCAPiPlus);
 288   TH1F* fHistDCAPiMinus = new TH1F("fHistDCAPiMinus","",100,0,10);
 289   fOutputList->Add(fHistDCAPiMinus);
 290   TH1F* fHistDCADaughters = new TH1F("fHistDCADaughters", "DCA of pions to each other", 50, 0., 0.5);
 291   fOutputList->Add(fHistDCADaughters);
 292   TH2F* fHistK0PiPlusPt = new TH2F("fHistK0PiPlusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
 293   fOutputList->Add(fHistK0PiPlusPt);
 294   TH2F* fHistK0PiMinusPt = new TH2F("fHistK0PiMinusPt", "", kCentBins, .5,kCentBins+.5, 40,0.,4.);
 295   fOutputList->Add(fHistK0PiMinusPt);
 296   TH1F* fHistDaughterPhi = new TH1F("fHistDaughterPhi","",180,-PI,PI);
 297   fOutputList->Add(fHistDaughterPhi);
 298
 299   //K0 parameters
 300   TH1F* fHistMultK0 = new TH1F("fHistMultK0", "K0 multiplicity", 51, -0.5, 51-0.5);
 301   fOutputList->Add(fHistMultK0);
 302   TH2F* fHistPtK0 = new TH2F("fHistPtK0", "K0 pt distribution",kCentBins,.5,kCentBins+.5, 100, 0., 10.);
 303   fOutputList->Add(fHistPtK0);
 304   TH1F* fHistDecayLengthK0 = new TH1F("fHistDecayLengthK0", "K0 decay length", 100, 0., 100.);
 305   fOutputList->Add(fHistDecayLengthK0);
 306   TH1F* fHistDCAK0 = new TH1F("fHistDCAK0", "DCA of K0 to primary vertex", 40, 0., 0.4);
 307   fOutputList->Add(fHistDCAK0);
 308   TH2F* fHistKtK0 = new TH2F("fHistKtK0", "Kt distribution of K0 pairs", kCentBins, .5, kCentBins+.5, 300, 0., 3.);
 309   fOutputList->Add(fHistKtK0);
 310
 311   TH1F* fHistPx = new TH1F("fHistPx","",200,0,2);
 312   TH1F* fHistPy = new TH1F("fHistPy","",200,0,2);
 313   TH1F* fHistPz = new TH1F("fHistPz","",200,0,2);
 314   TH1F* fHistPxCM = new TH1F("fHistPxCM","",200,0,2);
 315   TH1F* fHistPyCM = new TH1F("fHistPyCM","",200,0,2);
 316   TH1F* fHistPzCM = new TH1F("fHistPzCM","",200,0,2);
 317   TH1F* fHistKsCM = new TH1F("fHistKsCM","",200,0,2);
 318   fOutputList->Add(fHistPx);
 319   fOutputList->Add(fHistPy);
 320   fOutputList->Add(fHistPz);
 321   fOutputList->Add(fHistPxCM);
 322   fOutputList->Add(fHistPyCM);
 323   fOutputList->Add(fHistPzCM);
 324   fOutputList->Add(fHistKsCM);
 325
 326   TH1F* fHistPOutLCMS = new TH1F("fHistPOutLCMS","",200,0,2);
 327   TH1F* fHistPSideLCMS = new TH1F("fHistPSideLCMS","",200,0,2);
 328   TH1F* fHistPLongLCMS = new TH1F("fHistPLongLCMS","",200,0,2);
 329   fOutputList->Add(fHistPOutLCMS);
 330   fOutputList->Add(fHistPSideLCMS);
 331   fOutputList->Add(fHistPLongLCMS);
 332
 333   //pair gamma (LCMS to PRF, OSL)
 334   TH2F* fHistGamma = new TH2F("fHistGamma","Gamma from LCMS to PRF",500,1,5,100,0,1);
 335   fOutputList->Add(fHistGamma);
 336
 337   //invariant mass distributions
 338   TH3F* fHistMass = new TH3F("fHistMass","",kCentBins,.5,kCentBins+.5,50,0.,5.,400,.3,.7);
 339   fOutputList->Add(fHistMass);
 340   //TH3F* fHistMassPtCFK0 = new TH3F("fHistMassPtCFK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
 341   //fOutputList->Add(fHistMassPtCFK0);
 342   //TH3F* fHistMassPtCFBkgK0 = new TH3F("fHistMassPtCFBkgK0","",kCentBins,.5,kCentBins+.5,50,0.,5.,200,.4,.6);
 343   //fOutputList->Add(fHistMassPtCFBkgK0);
 344   //TH3F* fHistMassQKt = new TH3F("fHistMassQKt","",100,0,1,200,0,2,200,.4,.6);
 345   //fOutputList->Add(fHistMassQKt);
 346   //TH3F* fHistMassKtK0 = new TH3F("fHistMassKtK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
 347   //fOutputList->Add(fHistMassKtK0);
 348   //TH3F* fHistMassKtBkgK0 = new TH3F("fHistMassKtBkgK0","",kCentBins,.5,kCentBins+.5,300,0.,3.,200,.4,.6);
 349   //fOutputList->Add(fHistMassKtBkgK0);
 350
 351   //separation studies
 352   TH1F* fHistSepNumPos = new TH1F("fHistSepNumPos","",200,0,20);
 353   fOutputList->Add(fHistSepNumPos);
 354   TH1F* fHistSepDenPos = new TH1F("fHistSepDenPos","",200,0,20);
 355   fOutputList->Add(fHistSepDenPos);
 356   TH1F* fHistSepNumNeg = new TH1F("fHistSepNumNeg","",200,0,20);
 357   fOutputList->Add(fHistSepNumNeg);
 358   TH1F* fHistSepDenNeg = new TH1F("fHistSepDenNeg","",200,0,20);
 359   fOutputList->Add(fHistSepDenNeg);
 360
 361   TH2F* fHistSepNumPos2 = new TH2F("fHistSepNumPos2","",100,0,20,100,0,20);
 362   TH2F* fHistSepDenPos2 = new TH2F("fHistSepDenPos2","",100,0,20,100,0,20);
 363   TH2F* fHistSepNumNeg2 = new TH2F("fHistSepNumNeg2","",100,0,20,100,0,20);
 364   TH2F* fHistSepDenNeg2 = new TH2F("fHistSepDenNeg2","",100,0,20,100,0,20);
 365   fOutputList->Add(fHistSepNumPos2);
 366   fOutputList->Add(fHistSepDenPos2);
 367   fOutputList->Add(fHistSepNumNeg2);
 368   fOutputList->Add(fHistSepDenNeg2);
 369
 370   TH2F* fHistSepDPC = new TH2F("fHistSepDPC","",200,-1,1,50,0,10);
 371   TH2F* fHistSepDPCBkg = new TH2F("fHistSepDPCBkg","",200,-1,1,50,0,10);
 372   fOutputList->Add(fHistSepDPC);
 373   fOutputList->Add(fHistSepDPCBkg);
 374
 375   TH1F *fHistPsi = new TH1F("fHistPsi","",90,-PI/2,PI/2);
 376   fOutputList->Add(fHistPsi);
 377   TH2F *fHistPhi = new TH2F("fHistPhi","",kCentBins,.5,kCentBins+.5,180,0,2*PI);
 378   fOutputList->Add(fHistPhi);
 379   TH2F *fHistPhiPsi = new TH2F("fHistPhiPsi","",kCentBins,.5,kCentBins+.5,180,0,2*PI);
 380   fOutputList->Add(fHistPhiPsi);
 381
 382
 383   TH3F *fHistDPhi = new TH3F("fHistDPhi","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
 384   TH3F *fHistDPhiBkg = new TH3F("fHistDPhiBkg","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
 385   TH3F *fHistDPhiPsi = new TH3F("fHistDPhiPsi","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
 386   TH3F *fHistDPhiPsiBkg = new TH3F("fHistDPhiPsiBkg","",kCentBins,.5,kCentBins+.5,200,0.,2.,90,0,PI);
 387   fOutputList->Add(fHistDPhi);
 388   fOutputList->Add(fHistDPhiBkg);
 389   fOutputList->Add(fHistDPhiPsi);
 390   fOutputList->Add(fHistDPhiPsiBkg);
 391
 392 /////////Pair Distributions///////////////////
 393
 394   //1D Q invariant
 395   TH3F* fHistQinvSignal = new TH3F("fHistQinvSignal","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 396   fOutputList->Add(fHistQinvSignal);
 397   TH3F* fHistQinvBkg = new TH3F("fHistQinvBkg","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
 398   fOutputList->Add(fHistQinvBkg);
 399
 400   //event plane
 401   TH3F* fHistQinvSignalEPIn = new TH3F("fHistQinvSignalEPIn","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 402   fOutputList->Add(fHistQinvSignalEPIn);
 403   TH3F* fHistQinvBkgEPIn = new TH3F("fHistQinvBkgEPIn","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
 404   fOutputList->Add(fHistQinvBkgEPIn);
 405   TH3F* fHistQinvSignalEPOut = new TH3F("fHistQinvSignalEPOut","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 406   fOutputList->Add(fHistQinvSignalEPOut);
 407   TH3F* fHistQinvBkgEPOut = new TH3F("fHistQinvBkgEPOut","Mixed Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1.);
 408   fOutputList->Add(fHistQinvBkgEPOut);
 409
 410   //mass bins within peak
 411   //TH3F* fHistCLCLSignal = new TH3F("fHistCLCLSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 412   //TH3F* fHistCLCLBkg = new TH3F("fHistCLCLBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 413   //TH3F* fHistCLCRSignal = new TH3F("fHistCLCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 414   //TH3F* fHistCLCRBkg = new TH3F("fHistCLCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 415   //TH3F* fHistCRCRSignal = new TH3F("fHistCRCRSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 416   //TH3F* fHistCRCRBkg = new TH3F("fHistCRCRBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 417   //fOutputList->Add(fHistCLCLSignal);
 418   //fOutputList->Add(fHistCLCLBkg);
 419   //fOutputList->Add(fHistCLCRSignal);
 420   //fOutputList->Add(fHistCLCRBkg);
 421   //fOutputList->Add(fHistCRCRSignal);
 422   //fOutputList->Add(fHistCRCRBkg);
 423
 424   //3D out-side-long
 425   TH3F *fHist3DOSLSignal[10][4];
 426   TH3F *fHist3DOSLBkg[10][4];
 427
 428   if(fCase3D){
 429    for(int i3D=0;i3D<10;i3D++){
 430     for(int j3D=0;j3D<4;j3D++){
 431      TString *histname = new TString("fHist3DOSL");
 432      *histname += i3D;
 433      *histname += j3D;
 434      histname->Append("Signal");
 435      fHist3DOSLSignal[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
 436      fOutputList->Add(fHist3DOSLSignal[i3D][j3D]);
 437      histname->Replace(12,6,"Bkg");
 438      fHist3DOSLBkg[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
 439      fOutputList->Add(fHist3DOSLBkg[i3D][j3D]);
 440     }
 441    }
 442   }
 443
 444   //3D Spherical Harmonics
 445   //TH3F* fHistSHCentLowKt = new TH3F("fHistSHCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 446   //TH3F* fHistSHCentHighKt = new TH3F("fHistSHCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 447   //TH3F* fHistSHSemiCentLowKt = new TH3F("fHistSHSemiCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 448   //TH3F* fHistSHSemiCentHighKt = new TH3F("fHistSHSemiCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 449   //TH3F* fHistSHCentLowKtBkg = new TH3F("fHistSHCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 450   //TH3F* fHistSHCentHighKtBkg = new TH3F("fHistSHCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 451   //TH3F* fHistSHSemiCentLowKtBkg = new TH3F("fHistSHSemiCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 452   //TH3F* fHistSHSemiCentHighKtBkg = new TH3F("fHistSHSemiCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 453   //fOutputList->Add(fHistSHCentLowKt);
 454   //fOutputList->Add(fHistSHCentHighKt);
 455   //fOutputList->Add(fHistSHSemiCentLowKt);
 456   //fOutputList->Add(fHistSHSemiCentHighKt);
 457   //fOutputList->Add(fHistSHCentLowKtBkg);
 458   //fOutputList->Add(fHistSHCentHighKtBkg);
 459   //fOutputList->Add(fHistSHSemiCentLowKtBkg);
 460   //fOutputList->Add(fHistSHSemiCentHighKtBkg);
 461
 462   //side-side
 463   //TH3F* fHistLeftLeftSignal = new TH3F("fHistLeftLeftSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 464   //TH3F* fHistLeftRightSignal = new TH3F("fHistLeftRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 465   //TH3F* fHistRightRightSignal = new TH3F("fHistRightRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 466   //TH3F* fHistLeftLeftBkg = new TH3F("fHistLeftLeftBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 467   //TH3F* fHistLeftRightBkg = new TH3F("fHistLeftRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 468   //TH3F* fHistRightRightBkg = new TH3F("fHistRightRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 469   //fOutputList->Add(fHistLeftLeftSignal);
 470   //fOutputList->Add(fHistLeftRightSignal);
 471   //fOutputList->Add(fHistRightRightSignal);
 472   //fOutputList->Add(fHistLeftLeftBkg);
 473   //fOutputList->Add(fHistLeftRightBkg);
 474   //fOutputList->Add(fHistRightRightBkg);
 475
 476   //TH3F* fHistSplitK0Sides = new TH3F("fHistSplitK0Sides","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 477   //fOutputList->Add(fHistSplitK0Sides);
 478   //TH3F* fHistSplitK0Centers = new TH3F("fHistSplitK0Centers","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 479   //fOutputList->Add(fHistSplitK0Centers);
 480   //TH3F* fHistQinvSignalNoSplit = new TH3F("fHistQinvSignalNoSplit","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 481   //fOutputList->Add(fHistQinvSignalNoSplit);
 482
 483   PostData(1, fOutputList);
 484
 485 }
 486
 487 //________________________________________________________________________
 488 void AliFemtoK0Analysis::Exec(Option_t *)
 489 {
 490   // Main loop
 491   // Called for each event
 492   //cout<<"===========  Event # "<<fEventCount+1<<"  ==========="<<endl;
 493   fEventCount++;
 494   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
 495   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
 496
 497   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral));
 498   bool isCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
 499   //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 500   if(!isSelected) {
 501    //cout << "Failed trigger selection." << endl;
 502    return;
 503   }
 504
 505   ///////////////////////////////////////////////////////////
 506
 507   unsigned int statusPos=0;
 508   unsigned int statusNeg=0;
 509
 510   float bField=0;
 511   bField = fAOD->GetMagneticField();
 512   if(bField == 0) return;
 513   if(fSignDep){
 514    if(fFieldPos && bField < 0) return;
 515    if(!fFieldPos && bField > 0) return;
 516   }
 517
 518
 519   /////////////////////////////////////////////////
 520
 521
 522   //Centrality selection
 523
 524   AliCentrality *centrality = fAOD->GetCentrality();
 525   float percent = centrality->GetCentralityPercentile("V0M");
 526   int centBin=0;
 527   //Printf("Centrality percent = %f", percent);
 528
 529   AliAODVZERO *aodV0 = fAOD->GetVZEROData();
 530   float multV0A=aodV0->GetMTotV0A();
 531   float multV0C=aodV0->GetMTotV0C();
 532
 533   if(percent < 0) {
 534    //Printf("No centrality info");
 535    return;
 536   }
 537   if(percent < 0.1 && (multV0A + multV0C < 19500)){
 538    //Printf("No centrality info");
 539    return;
 540   }
 541   else if(percent <= 5)   centBin=15;
 542   else if(percent <= 10)  centBin=14;
 543   else if(percent <= 15)  centBin=13;
 544   else if(percent <= 20)  centBin=12;
 545   else if(percent <= 25)  centBin=11;
 546   else if(percent <= 30)  centBin=10;
 547   else if(percent <= 35)  centBin=9;
 548   else if(percent <= 40)  centBin=8;
 549   else if(percent <= 45)  centBin=7;
 550   else if(percent <= 50)  centBin=6;
 551   else if(percent <= 55)  centBin=5;
 552   else if(percent <= 60)  centBin=4;
 553   else if(percent <= 65)  centBin=3;
 554   else if(percent <= 70)  centBin=2;
 555   else if(percent <= 75)  centBin=1;
 556   else if(percent <= 80)  centBin=0;
 557   else {
 558    //Printf("Skipping Peripheral Event");
 559    return;
 560   }
 561   if(percent > 10 && isCentral) return;
 562   ((TH1F*)fOutputList->FindObject("fHistCent"))->Fill(percent);
 563
 564   //flatten centrality dist.
 565   if(percent < 9){
 566    if(fFlatCent){
 567     if(RejectEventCentFlat(bField,percent)) return;
 568    }
 569   }
 570   ((TH1F*)fOutputList->FindObject("fHistCentFlat"))->Fill(percent);
 571
 572   //Vertexing
 573   AliAODVertex *primaryVertex;
 574   double vertex[3]={0};
 575   primaryVertex = fAOD->GetPrimaryVertex();
 576   vertex[0]=primaryVertex->GetX();
 577   vertex[1]=primaryVertex->GetY();
 578   vertex[2]=primaryVertex->GetZ();
 579   if(vertex[0]<10e-5 && vertex[1]<10e-5 &&  vertex[2]<10e-5) return;
 580   if(fabs(vertex[2]) > 10) return; // Z-vertex Cut
 581
 582   int zBin=0;
 583   double zStep=2*10/double(kZVertexBins), zStart=-10.;
 584   for(int i=0; i<kZVertexBins; i++)
 585   {
 586    if((vertex[2] > zStart+i*zStep) && (vertex[2] < zStart+(i+1)*zStep))
 587    {
 588     zBin=i;
 589     break;
 590    }
 591   }
 592
 593   //Event plane
 594   int psiBin = 0;
 595   AliEventplane *eventplane = fAOD->GetEventplane();
 596   if(fPsiBinning && !eventplane) return;
 597   double psiEP = eventplane->GetEventplane("V0",fAOD,2); //[-PI/2,PI/2]
 598   ((TH1F*)fOutputList->FindObject("fHistPsi"))->Fill(psiEP);
 599
 600   double psiStep = PI/double(fNPsiBins);
 601   double psiStart = -0.5*PI;
 602   for(int i=0; i<fNPsiBins; i++)
 603   {
 604    if((psiEP > psiStart+i*psiStep) && (psiEP < psiStart+(i+1)*psiStep))
 605    {
 606     psiBin = i;
 607     break;
 608    }
 609   }
 610   if(!fPsiBinning) psiBin = 0;
 611
 612 ////////////////////////////////////////////////////////////////
 613 //Cut Values and constants
 614
 615   //const bool kMCCase = kFALSE;                     //switch for MC analysis
 616   const int kMaxNumK0 = 300;                       //maximum number of K0s, array size
 617   const float kMinDCAPrimaryPion = 0.4;            //minimum dca of pions to primary
 618   const float kMaxDCADaughtersK0 = 0.3;            //maximum dca of pions to each other - 3D
 619   const float kMaxDCAK0 = 0.3;                     //maximum dca of K0 to primary
 620   const float kMaxDLK0 = 30.0;                     //maximum decay length of K0
 621   const float kMinDLK0 = fMinDecayLength;          //minimum decay length of K0
 622   const float kEtaCut = 0.8;                       //maximum |pseudorapidity|
 623   const float kMinCosAngle = 0.99;                 //minimum cosine of K0 pointing angle
 624
 625   const float kMinSeparation = fMinSep;                //minimum daughter (pair) separation
 626
 627   const float kTOFLow = 0.8;                       //boundary for TOF usage
 628   const float kMaxTOFSigmaPion = 3.0;              //TOF # of sigmas
 629   const float kMaxTPCSigmaPion = 3.0;              //TPC # of sigmas
 630
 631   //const float kMassPion = .13957;
 632   const float kMassK0Short = .497614;       //true PDG masses
 633
 634 ////////////////////////////////////////////////////////////////
 635   //v0 tester
 636 ////////////////////////////////////////////////////////////////
 637   int v0Count = 0;      //number of v0s (entries in array)
 638   int k0Count = 0;      //number of good K0s
 639
 640   AliFemtoK0Particle *tempK0 = new AliFemtoK0Particle[kMultLimit];
 641
 642   //for daughter sharing studies
 643   //int idArray[100] = {0};
 644   //int idCount = 0;
 645
 646   //for MC
 647   //TClonesArray *mcArray = 0x0;
 648   //if(kMCCase){
 649   //mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 650   //if(!mcArray){cout<<"No MC particle branch found"<<endl;return;}}
 651
 652   for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
 653   {
 654     bool goodK0 = kFALSE;
 655     bool goodPiPlus = kFALSE;
 656     bool goodPiMinus = kFALSE;
 657
 658     //load v0 track
 659     AliAODv0* v0 = fAOD->GetV0(i);
 660     if(!v0) continue;
 661     if(fOnlineCase){
 662      if(!(v0->GetOnFlyStatus())) continue;
 663     } //for online
 664     else{
 665      if((v0->GetOnFlyStatus())) continue; //for offline
 666     }
 667
 668     //for on-the-fly ordering
 669     AliAODTrack* tempTrack = (AliAODTrack*)v0->GetDaughter(0);
 670     short int pos0or1;
 671     short int neg0or1;
 672     bool orderswitch = kFALSE;
 673     if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}
 674     else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}
 675
 676     //load daughter tracks
 677     AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);
 678     AliAODTrack* prongTrackNeg = (AliAODTrack*)v0->GetDaughter(neg0or1);
 679     if(!prongTrackPos) continue;
 680     if(!prongTrackNeg) continue;
 681
 682     //daughter cuts
 683     if(v0->PtProng(pos0or1) < .15) continue;
 684     if(v0->PtProng(neg0or1) < .15) continue;
 685     if(fabs(v0->EtaProng(pos0or1)) > .8) continue;
 686     if(fabs(v0->EtaProng(neg0or1)) > .8) continue;
 687
 688     //load status for PID
 689     statusPos=prongTrackPos->GetStatus();
 690     if((statusPos&AliESDtrack::kTPCrefit)==0) continue;
 691     prongTrackPos->SetAODEvent(fAOD);
 692     statusNeg=prongTrackNeg->GetStatus();
 693     if((statusNeg&AliESDtrack::kTPCrefit)==0) continue;
 694     prongTrackNeg->SetAODEvent(fAOD);
 695
 696     //TPC PID
 697     if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackPos,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiPlus = kTRUE;
 698     if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;
 699
 700     //Positive daughter identification TOF
 701     float probMis;
 702     AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);
 703     double Ppos = v0->PProng(pos0or1);
 704     if(Ppos > kTOFLow) //PiPlus
 705     {
 706      //if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
 707      if(AliPIDResponse::kDetPidOk == statusPosTOF)
 708      {
 709       probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
 710       if(probMis < 0.01) //avoid TOF-TPC mismatch
 711       {
 712        if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;
 713        else goodPiPlus = kFALSE;
 714       }
 715      }
 716     }
 717     //Negative daughter identification TOF
 718     AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);
 719     double Pneg = v0->PProng(neg0or1);
 720     if(Pneg > kTOFLow) //PiMinus
 721     {
 722      //if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
 723      if(AliPIDResponse::kDetPidOk == statusNegTOF)
 724      {
 725       probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
 726       if(probMis < 0.01) //avoid TOF-TPC mismatch
 727       {
 728        if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;
 729        else goodPiMinus = kFALSE;
 730       }
 731      }
 732     }
 733
 734     //K0 cuts
 735     if(v0->Eta() > kEtaCut)                                     continue;
 736     if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle)      continue;
 737     if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8)        continue;
 738     if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion)           continue;
 739     if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion)           continue;
 740     if(v0->DecayLength(primaryVertex) > kMaxDLK0)               continue;
 741     if(v0->DecayLength(primaryVertex) < kMinDLK0)                       continue;
 742     if(v0->DcaV0Daughters() > kMaxDCADaughtersK0)               continue;
 743     double v0Dca = v0->DcaV0ToPrimVertex();
 744     if(v0Dca > kMaxDCAK0)                                                       continue;
 745     if(!goodPiMinus || !goodPiPlus)                             continue;
 746
 747     //EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT
 748
 749     //for MC
 750     //bool MCgood = kFALSE;
 751     //if(kMCCase){
 752     //AliAODMCParticle* mck0dp = (AliAODMCParticle*)mcArray->At(abs(prongTrackPos->GetLabel()));
 753     //AliAODMCParticle* mck0dn = (AliAODMCParticle*)mcArray->At(abs(prongTrackNeg->GetLabel()));
 754     //if(mck0dp->GetMother() >= 0){
 755      //if(mck0dp->GetMother() == mck0dn->GetMother()){
 756       //if(abs(mck0dp->GetPdgCode()) == 211 && abs(mck0dn->GetPdgCode()) == 211){
 757        //AliAODMCParticle* mck0 = (AliAODMCParticle*)mcArray->At(mck0dp->GetMother());
 758        //if(abs(mck0->GetPdgCode()) == 310){
 759         //MCgood = kTRUE;
 760        //}
 761       //}
 762      //}
 763     //}
 764     //}// if kMCCase
 765
 766     if(v0->MassK0Short() > .48 && v0->MassK0Short() < .515) goodK0 = kTRUE;
 767     //else continue; //removed, Apr 18
 768
 769     //Check for shared daughters, using v0 DCA to judge
 770     bool v0JudgeNew; //true if new v0 beats old
 771     tempK0[v0Count].fSkipShared = kFALSE;
 772     double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut
 773     if(fMeritCase){
 774      for(int iID = 0; iID<v0Count; iID++){
 775       if(tempK0[iID].fSkipShared == kFALSE){            //if old is already skipped, go to next old
 776        if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){
 777         double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca};
 778         v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins
 779         if(!v0JudgeNew){                //if old beats new...
 780          if(!tempK0[iID].fK0 && goodK0) continue;       //if bad old beats new good, do nothing...
 781          else{                                          //but if bad old beats new bad, or good old beats anything, skip new
 782           tempK0[v0Count].fSkipShared = kTRUE;          //skip new
 783           break;                                        //no need to keep checking others
 784          }
 785         }
 786         else{                                           //if new beats old...
 787          if(tempK0[iID].fK0 && !goodK0) continue;       //if bad new beats good old, do nothing...
 788          else{                                          //but if bad new beats bad old, or good new beats anything, skip old
 789               tempK0[iID].fSkipShared = kTRUE;          //skip old
 790               if(tempK0[iID].fK0) k0Count--;            //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)
 791          }
 792         }
 793        }
 794       }
 795      }
 796      if(tempK0[v0Count].fSkipShared) continue;          //if new K0 is skipped, don't load; go to next v0
 797     }//if MeritCase
 798
 799     //load parameters into temporary class instance
 800     if(v0Count < kMaxNumK0)
 801     {
 802         if(goodK0){
 803          tempK0[v0Count].fK0 = kTRUE;
 804          k0Count++;
 805         }
 806         else tempK0[v0Count].fK0 = kFALSE;
 807
 808         //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;
 809         //else tempK0[v0Count].fSideLeft = kFALSE;
 810         //if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;
 811         //else tempK0[v0Count].fSideRight = kFALSE;
 812             //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)
 813
 814         tempK0[v0Count].fDaughterID1    = prongTrackPos->GetID();
 815         tempK0[v0Count].fDaughterID2    = prongTrackNeg->GetID();
 816         tempK0[v0Count].fMomentum[0]    = v0->Px();
 817         tempK0[v0Count].fMomentum[1]    = v0->Py();
 818         tempK0[v0Count].fMomentum[2]    = v0->Pz();
 819         tempK0[v0Count].fPt             = v0->Pt();
 820         tempK0[v0Count].fMass           = v0->MassK0Short();
 821         tempK0[v0Count].fV0Dca          = v0Dca;
 822
 823         //for hists
 824         tempK0[v0Count].fDDDca          = v0->DcaV0Daughters();
 825             tempK0[v0Count].fDecayLength        = v0->DecayLength(primaryVertex);
 826         tempK0[v0Count].fPosPt          = v0->PtProng(pos0or1);
 827         tempK0[v0Count].fNegPt          = v0->PtProng(neg0or1);
 828         tempK0[v0Count].fPosPhi         = v0->PhiProng(pos0or1);
 829         tempK0[v0Count].fNegPhi         = v0->PhiProng(neg0or1);
 830             if(!orderswitch){
 831          tempK0[v0Count].fPosDca        = v0->DcaPosToPrimVertex();
 832          tempK0[v0Count].fNegDca        = v0->DcaNegToPrimVertex();
 833             }
 834         else{
 835          tempK0[v0Count].fPosDca        = v0->DcaNegToPrimVertex();
 836          tempK0[v0Count].fNegDca        = v0->DcaPosToPrimVertex();
 837         }
 838
 839                 //for psi studies
 840         double v0Phi = v0->Phi(); //between [0,2pi]
 841         double v0PhiPsi = v0Phi-psiEP;
 842         if(v0PhiPsi < 0) v0PhiPsi += 2.*PI;
 843         else if (v0PhiPsi > 2.*PI) v0PhiPsi -= 2.*PI;
 844                 else{};
 845         tempK0[v0Count].fPhi    = v0Phi;
 846         tempK0[v0Count].fPhiPsi = v0PhiPsi;
 847
 848         //for separation
 849         GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);
 850             GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);
 851         //for DPC
 852         prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);
 853         prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);
 854
 855         //if(idCount < 50){
 856         // if(goodK0){
 857         //  idArray[idCount*2]   = prongTrackPos->GetID();
 858         //  idArray[idCount*2+1] = prongTrackNeg->GetID();
 859         //  idCount++;
 860         //}}
 861
 862         v0Count++;
 863     }
 864
 865   }//v0
 866
 867   if(k0Count<2) return;  //only keep events with more than 1 good K0
 868
 869   //Add Event to buffer - this is for event mixing
 870   fEC[zBin][centBin][psiBin]->FIFOShift();
 871   (fEvt) = fEC[zBin][centBin][psiBin]->fEvt;
 872   (fEvt)->fFillStatus = 1;
 873   int unskippedCount = 0;
 874   for(int i=0;i<v0Count;i++)
 875   {
 876    if(!tempK0[i].fSkipShared)                           //don't include skipped v0s (from shared daughters)
 877    {
 878     ((TH3F*)fOutputList->FindObject("fHistMass"))->Fill(centBin+1,tempK0[i].fPt,tempK0[i].fMass);
 879     if(tempK0[i].fK0)                                   //make sure particle is good (mass)
 880     {
 881      (fEvt)->fK0Particle[unskippedCount] = tempK0[i];   //load good, unskipped K0s
 882      unskippedCount++;                                  //count good, unskipped K0s
 883     }
 884    }
 885   }
 886   (fEvt)->fNumV0s = unskippedCount;
 887   //Printf("Number of v0s: %d", v0Count);
 888   //Printf("Number of K0s: %d", k0Count);
 889   delete [] tempK0; tempK0 = NULL;
 890
 891   ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount);        // changed 3/25, used to be "k0Count"
 892   ((TH1F*)fOutputList->FindObject("fHistCentUsed"))->Fill(percent);
 893
 894   //Printf("Reconstruction Finished. Starting pair studies.");
 895
 896   //////////////////////////////////////////////////////////////////////
 897   // Correlations
 898   //////////////////////////////////////////////////////////////////////
 899
 900   float px1, py1, pz1, px2, py2, pz2;                   //single kaon values
 901   float en1, en2;                                                               //single kaon values
 902   //float pt1, pt2;                                                             //single kaon values
 903   float pairPx, pairPy, pairPz, pairP0;                 //pair momentum values
 904   float pairPt, pairMt, pairKt;                         //pair momentum values
 905   float pairMInv, pairPDotQ;
 906   float qinv, q0, qx, qy, qz;                           //pair q values
 907   //float qLength, thetaSH, thetaSHCos, phiSH;            //Spherical Harmonics values
 908   float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks;    //PRF
 909   //float qOutLCMS;
 910   float qOutPRF, qSide, qLong;                          //relative momentum in LCMS/PRF frame
 911   float betasq, gamma;
 912   float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;
 913   float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;
 914
 915
 916   for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0
 917   {
 918     //single particle histograms (done here to avoid "skipped" v0s
 919     ((TH1F*)fOutputList->FindObject("fHistDCADaughters"))       ->Fill((fEvt)->fK0Particle[i].fDDDca);
 920     ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0"))      ->Fill((fEvt)->fK0Particle[i].fDecayLength);
 921     ((TH1F*)fOutputList->FindObject("fHistDCAK0"))              ->Fill((fEvt)->fK0Particle[i].fV0Dca);
 922     ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus")) ->Fill((fEvt)->fK0Particle[i].fNegDca);
 923     ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus"))          ->Fill((fEvt)->fK0Particle[i].fPosDca);
 924     ((TH2F*)fOutputList->FindObject("fHistPtK0"))               ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);
 925     ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);
 926     ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt"))        ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);
 927     ((TH1F*)fOutputList->FindObject("fHistDaughterPhi"))        ->Fill((fEvt)->fK0Particle[i].fPosPhi);
 928     ((TH1F*)fOutputList->FindObject("fHistDaughterPhi"))        ->Fill((fEvt)->fK0Particle[i].fNegPhi);
 929
 930     ((TH1F*)fOutputList->FindObject("fHistPx"))         ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);
 931     ((TH1F*)fOutputList->FindObject("fHistPy"))         ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);
 932     ((TH1F*)fOutputList->FindObject("fHistPz"))         ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);
 933
 934     ((TH2F*)fOutputList->FindObject("fHistPhi"))        ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhi);
 935     ((TH2F*)fOutputList->FindObject("fHistPhiPsi"))     ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhiPsi);
 936
 937     for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events
 938     {
 939       int startbin=0;
 940       if(evnum==0) startbin=i+1;
 941
 942       for(int j=startbin; j<(fEvt+evnum)->fNumV0s; j++) // Past event V0
 943       {
 944         if(evnum==0)  // Get rid of shared tracks
 945         {
 946           if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
 947                   if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
 948               if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
 949               if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
 950             }
 951
 952         px1 = (fEvt)->fK0Particle[i].fMomentum[0];
 953                 py1 = (fEvt)->fK0Particle[i].fMomentum[1];
 954                 pz1 = (fEvt)->fK0Particle[i].fMomentum[2];
 955         //pt1 = (fEvt)->fK0Particle[i].fPt;
 956                 px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];
 957                 py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];
 958                 pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];
 959         //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;
 960         if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias
 961                  double tempvar;
 962          tempvar = px1; px1 = px2; px2 = tempvar;
 963          tempvar = py1; py1 = py2; py2 = tempvar;
 964          tempvar = pz1; pz1 = pz2; pz2 = tempvar;
 965                 }
 966
 967                 en1  = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));
 968                 en2  = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));
 969
 970         q0 = en1 - en2;
 971         qx = px1 - px2;
 972         qy = py1 - py2;
 973         qz = pz1 - pz2;
 974         qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));
 975
 976         pairPx = px1 + px2;
 977                 pairPy = py1 + py2;
 978                 pairPz = pz1 + pz2;
 979         pairP0 = en1 + en2;
 980                 pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);
 981         pairKt = pairPt/2.;                                                                     //used for KT binning
 982         pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz);           //used for LCMS (not plots)
 983                 pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF
 984         pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz;    //used for PRF
 985
 986             //PRF (this section will probably be removed in favor of later boosting section)
 987         p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2));  //pair momentum length
 988         am12 = sqrt(pow(en1+en2,2)-p12*p12);                                    //sqrt(s)=|p1+p2|(4vec)
 989         epm = en1+en2+am12;                                                                             //"energy plus mass"
 990         p112 = px1*pairPx+py1*pairPy+pz1*pairPz;                                //proj. of p1 on pairP
 991         if(am12 == 0) continue;
 992         h1 = (p112/epm - en1)/am12;
 993         ppx = px1+pairPx*h1;                                                                    //px in PRF
 994         ppy = py1+pairPy*h1;                                                                    //py in PRF
 995         ppz = pz1+pairPz*h1;                                                                    //pz in PRF
 996         ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz);                                             //k*
 997         ((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);
 998         ((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);
 999         ((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);
1000         ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);
1001
1002         //relative momentum in out-side-long for LCMS and PRF
1003         if(pairMt == 0 || pairPt == 0) continue;
1004         qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames
1005         qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames
1006         //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;
1007                 qOutPRF  = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;
1008
1009                 //finding gamma for gamma binning/hists (likely will be removed after tests)
1010                 p1LCMSOut  = (pairPx*px1+pairPy*py1)/pairPt;
1011                 p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;
1012                 p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;
1013                 p2LCMSOut  = (pairPx*px2+pairPy*py2)/pairPt;
1014                 p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;
1015                 p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;
1016                 en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));
1017                 en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2));
1018                 betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);
1019                 gamma = 1./sqrt(1-betasq);
1020                 ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);
1021                 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);
1022                 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);
1023                 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);
1024                 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);
1025                 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);
1026                 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);
1027                 //getting bin numbers and names for 3D histogram
1028         TString *histname3D = new TString("fHist3DOSL");
1029         int ktBin;
1030         if(pairKt < 0.6) ktBin = 0;
1031                 else if(pairKt < 0.8) ktBin = 1;
1032                 else if(pairKt < 1.0) ktBin = 2;
1033                 else ktBin = 3;
1034                 *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]
1035                 *histname3D += ktBin;
1036
1037         //Spherical harmonics
1038         //qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);
1039         //thetaSHCos = qLong/qLength;
1040         //thetaSH = acos(thetaSHCos);
1041         //phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));
1042
1043         //Finding average separation of daughters throughout TPC - two-track cut
1044         float posPositions1[9][3] = {{0}};
1045         float negPositions1[9][3] = {{0}};
1046         float posPositions2[9][3] = {{0}};
1047         float negPositions2[9][3] = {{0}};
1048         for(int iPos = 0; iPos < 9; iPos++){
1049          for(int jPos = 0; jPos < 3; jPos++){
1050            posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];
1051            negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];
1052            posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];
1053            negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];
1054          }
1055         }
1056         float pMean = 0.;       //average separation for positive daughters
1057         float nMean = 0.;       //average separation for negative daughters
1058         float pDiff;
1059         float nDiff;
1060         float pMin = 9999.;     //minimum separation (updates) - pos
1061         float nMin = 9999.;     //minimum separation (updates) - neg
1062         double pCount=0;        //counter for number of points used - pos
1063         double nCount=0;        //counter for number of points used - neg
1064         for(int ss=0;ss<9;ss++){
1065          if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){
1066           pCount++;
1067           pDiff = sqrt(pow(posPositions1[ss][0]-posPositions2[ss][0],2)+pow(posPositions1[ss][1]-posPositions2[ss][1],2)+pow(posPositions1[ss][2]-posPositions2[ss][2],2));
1068           pMean = pMean + pDiff;
1069           if(pDiff < pMin) pMin = pDiff;
1070          }
1071          if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){
1072           nCount++;
1073           nDiff = sqrt(pow(negPositions1[ss][0]-negPositions2[ss][0],2)+pow(negPositions1[ss][1]-negPositions2[ss][1],2)+pow(negPositions1[ss][2]-negPositions2[ss][2],2));
1074           nMean = nMean + nDiff;
1075           if(nDiff < nMin) nMin = nDiff;
1076          }
1077         }
1078         pMean = pMean/pCount;
1079         nMean = nMean/nCount;
1080
1081         if(evnum==0){
1082          ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean);
1083          ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);
1084          ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);
1085          ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);
1086         }
1087         else{
1088          ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean);
1089          ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);
1090          ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);
1091          ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);
1092         }
1093
1094         //Decay plane coincidence
1095         //daughter momenta
1096         float a1 = (fEvt)->fK0Particle[i].fPPos[0];
1097         float b1 = (fEvt)->fK0Particle[i].fPPos[1];
1098         float c1 = (fEvt)->fK0Particle[i].fPPos[2];
1099         float d1 = (fEvt)->fK0Particle[i].fPNeg[0];
1100         float e1 = (fEvt)->fK0Particle[i].fPNeg[1];
1101         float f1 = (fEvt)->fK0Particle[i].fPNeg[2];
1102         float a2 = (fEvt+evnum)->fK0Particle[j].fPPos[0];
1103         float b2 = (fEvt+evnum)->fK0Particle[j].fPPos[1];
1104         float c2 = (fEvt+evnum)->fK0Particle[j].fPPos[2];
1105         float d2 = (fEvt+evnum)->fK0Particle[j].fPNeg[0];
1106         float e2 = (fEvt+evnum)->fK0Particle[j].fPNeg[1];
1107         float f2 = (fEvt+evnum)->fK0Particle[j].fPNeg[2];
1108
1109         float cross1[3];
1110         float cross2[3];
1111         cross1[0] = b1*f1-c1*e1;
1112         cross1[1] = c1*d1-a1*f1;
1113         cross1[2] = a1*e1-b1*d1;
1114         cross2[0] = b2*f2-c2*e2;
1115         cross2[1] = c2*d2-a2*f2;
1116         cross2[2] = a2*e2-b2*d2;
1117         float crosslength1 = sqrt(pow(cross1[0],2)+pow(cross1[1],2)+pow(cross1[2],2));
1118         float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));
1119         float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);
1120
1121         if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);
1122         else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);
1123
1124         if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)
1125         //end separation studies
1126
1127         //Fill Histograms
1128         bool center1K0   = kFALSE;  //accepted mass K0
1129             bool center2K0   = kFALSE;
1130         if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;
1131         if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;
1132         //bool CL1 = kFALSE;
1133         //bool CL2 = kFALSE;
1134         //bool CR1 = kFALSE;
1135         //bool CR2 = kFALSE;
1136         //if(center1K0 && center2K0){
1137         // if((fEvt)->fK0Particle[i].fMass < kMassK0Short) CL1 = kTRUE;
1138         // else CR1 = kTRUE;
1139         // if((fEvt+evnum)->fK0Particle[j].fMass < kMassK0Short) CL2 = kTRUE;
1140         // else CR2 = kTRUE;
1141         //}
1142
1143         //bool SideLeft1 = kFALSE;
1144         //bool SideLeft2 = kFALSE;
1145         //bool SideRight1 = kFALSE;
1146         //bool SideRight2 = kFALSE;
1147         //if((fEvt)->fK0Particle[i].fSideLeft) SideLeft1 = kTRUE;
1148         //else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;
1149         //if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;
1150         //else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;
1151
1152         //for psi binning
1153                 float phipsi1 = (fEvt)->fK0Particle[i].fPhiPsi;
1154         float phipsi2 = (fEvt+evnum)->fK0Particle[j].fPhiPsi;
1155         bool inPlane1 = kFALSE;
1156         bool inPlane2 = kFALSE;
1157         if(phipsi1 > PI) phipsi1 = phipsi1-PI;
1158         if(phipsi2 > PI) phipsi2 = phipsi2-PI;
1159         if(phipsi1 < 0.25*PI || phipsi1 > 0.75*PI) inPlane1 = kTRUE;
1160         if(phipsi2 < 0.25*PI || phipsi2 > 0.75*PI) inPlane2 = kTRUE;
1161
1162                 //for dphi, dphipsi
1163                 float dPhi = fabs((fEvt)->fK0Particle[i].fPhi - (fEvt+evnum)->fK0Particle[j].fPhi);
1164         if(dPhi > PI) dPhi = 2*PI-dPhi;
1165         float dPhiPsi = fabs((fEvt)->fK0Particle[i].fPhiPsi - (fEvt+evnum)->fK0Particle[j].fPhiPsi);
1166         if(dPhiPsi > PI) dPhiPsi = 2*PI-dPhiPsi;
1167
1168         if(evnum==0) //Same Event
1169         {
1170           //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt)->fK0Particle[i].fMass);
1171           //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1172           //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1173           //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1174           //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1175           //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1176
1177           if(center1K0 && center2K0){
1178            //1D
1179            ((TH3F*)fOutputList->FindObject("fHistQinvSignal"))->Fill(centBin+1, pairKt, qinv);
1180            //if(!splitK0centers)((TH3F*)fOutputList->FindObject("fHistQinvSignalNoSplit"))->Fill(centBin+1, pairKt, qinv);
1181            ((TH2F*)fOutputList->FindObject("fHistKtK0"))->Fill(centBin+1, pairKt);
1182
1183                    //eventplane
1184            if(inPlane1 && inPlane2)
1185              ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPIn"))->Fill(centBin+1, pairKt, qinv);
1186            else if(!inPlane1 && !inPlane2)
1187              ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPOut"))->Fill(centBin+1, pairKt, qinv);
1188
1189            //dPhi,dPhiPsi
1190                    ((TH3F*)fOutputList->FindObject("fHistDPhi"))->Fill(centBin+1,pairKt,dPhi);
1191                    ((TH3F*)fOutputList->FindObject("fHistDPhiPsi"))->Fill(centBin+1,pairKt,dPhiPsi);
1192
1193
1194            //for mass bin study
1195            //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLSignal"))->Fill(centBin+1, pairKt, qinv);
1196            //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRSignal"))->Fill(centBin+1, pairKt, qinv);
1197            //else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);
1198
1199            //3D
1200            if(fCase3D){
1201             if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1202                      histname3D->Append("Signal");
1203                          ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1204                     }
1205            }
1206            /*if(pairKt < 1.0){
1207             if(centBin > 13){
1208              ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1209              ((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}
1210             else if(centBin > 9){
1211              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1212              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}}
1213            else if(pairKt < 2.0){
1214             if(centBin > 13){
1215              ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1216              ((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}
1217             else if(centBin > 9){
1218              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1219
1220              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1221
1222           }//centercenter
1223
1224           //side-side correlations
1225           //if(!splitK0sides){
1226           // if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftSignal"))->Fill(centBin+1, pairKt, qinv);
1227            //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightSignal"))->Fill(centBin+1, pairKt, qinv);
1228            //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightSignal"))->Fill(centBin+1, pairKt, qinv);
1229          //}
1230
1231         }//same event
1232
1233         else //Mixed Events
1234         {
1235           //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1236           //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1237           //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1238           //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1239
1240           if(center1K0 && center2K0){
1241            //1D
1242            ((TH3F*)fOutputList->FindObject("fHistQinvBkg"))->Fill(centBin+1, pairKt, qinv);
1243
1244                    //eventplane
1245            if(inPlane1 && inPlane2)
1246              ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPIn"))->Fill(centBin+1, pairKt, qinv);
1247            else if(!inPlane1 && !inPlane2)
1248              ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPOut"))->Fill(centBin+1, pairKt, qinv);
1249
1250                    ((TH3F*)fOutputList->FindObject("fHistDPhiBkg"))->Fill(centBin+1,pairKt,dPhi);
1251                    ((TH3F*)fOutputList->FindObject("fHistDPhiPsiBkg"))->Fill(centBin+1,pairKt,dPhiPsi);
1252
1253            //for mass bin study
1254            //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLBkg"))->Fill(centBin+1, pairKt, qinv);
1255            //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRBkg"))->Fill(centBin+1, pairKt, qinv);
1256            //else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);
1257
1258            //3D
1259            if(fCase3D){
1260             if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1261                      histname3D->Replace(12,6,"Bkg");
1262                          ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1263                     }
1264            }
1265            /*if(pairKt < 1.0){
1266             if(centBin > 13){
1267              ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1268              ((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}
1269             else if(centBin > 9){
1270              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1271              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}
1272            else if(pairKt < 2.0){
1273             if(centBin > 13){
1274              ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1275              ((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}
1276             else if(centBin > 9){
1277              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1278              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1279           }
1280
1281           //side-side correlations
1282           //if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftBkg"))->Fill(centBin+1, pairKt, qinv);
1283           //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightBkg"))->Fill(centBin+1, pairKt, qinv);
1284           //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightBkg"))->Fill(centBin+1, pairKt, qinv);
1285
1286         }//Mixed Events
1287
1288       }//past event
1289     }//event buffer
1290   }//current event
1291
1292   // Post output data.
1293   PostData(1, fOutputList);
1294
1295   }
1296 //________________________________________________________________________
1297 void AliFemtoK0Analysis::Terminate(Option_t *)
1298 {
1299   // Called once at the end of the query
1300   cout<<"Done"<<endl;
1301
1302 }
1303
1304 //_________________________________________________________________________
1305 void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){
1306   // Gets the global position of the track at nine different radii in the TPC
1307   // track is the track you want to propagate
1308   // bfield is the magnetic field of your event
1309   // globalPositionsAtRadii is the array of global positions in the radii and xyz
1310
1311   // Initialize the array to something indicating there was no propagation
1312   for(Int_t i=0;i<9;i++){
1313     for(Int_t j=0;j<3;j++){
1314       globalPositionsAtRadii[i][j]=-9999.;
1315     }
1316   }
1317
1318    // Make a copy of the track to not change parameters of the track
1319   AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1320   //printf("\nAfter CopyFromVTrack\n");
1321   //etp.Print();
1322
1323   // The global position of the the track
1324   Double_t xyz[3]={-9999.,-9999.,-9999.};
1325
1326   // Counter for which radius we want
1327   Int_t iR=0;
1328   // The radii at which we get the global positions
1329   // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1330   Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1331   // The global radius we are at
1332   Float_t globalRadius=0;
1333
1334   // Propagation is done in local x of the track
1335   for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1336     // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1337     // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1338     // If the track's momentum is smaller than infinite, it will develop a y-component, which
1339     // adds to the global radius
1340
1341     // Stop if the propagation was not succesful. This can happen for low pt tracks
1342     // that don't reach outer radii
1343     if(!etp.PropagateTo(x,bfield))break;
1344     etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1345     globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1346
1347     // Roughly reached the radius we want
1348     if(globalRadius > Rwanted[iR]){
1349
1350       // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1351       while (globalRadius>Rwanted[iR]){
1352         x-=.1;
1353         //      printf("propagating to x %5.2f\n",x);
1354         if(!etp.PropagateTo(x,bfield))break;
1355         etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1356         globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1357       }
1358       //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]);
1359       globalPositionsAtRadii[iR][0]=xyz[0];
1360       globalPositionsAtRadii[iR][1]=xyz[1];
1361       globalPositionsAtRadii[iR][2]=xyz[2];
1362       //subtract primary vertex, "zero" track for correct mixed-event comparison
1363       globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];
1364       globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];
1365       globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];
1366
1367       // Indicate we want the next radius
1368       iR+=1;
1369     }
1370     if(iR>=8){
1371       // TPC edge reached
1372       return;
1373     }
1374   }
1375 }
1376
1377 bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){
1378  //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.
1379  //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'
1380
1381  bool newV0Wins = kFALSE;
1382  double pardiff[3] = {newPars[0]-oldPars[0],
1383                       newPars[1]-oldPars[1],
1384                       newPars[2]-oldPars[2]};
1385  if(cutChoice > 0 && cutChoice < 4){
1386   if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;
1387  }
1388  else if(cutChoice == 4){
1389   int newWinCount = 0;
1390   for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}
1391   if(newWinCount > 1) newV0Wins = kTRUE;
1392  }
1393  else{};
1394  return newV0Wins;
1395 }
1396
1397 bool AliFemtoK0Analysis::RejectEventCentFlat(float MagField, float CentPercent)
1398 { // to flatten centrality distribution
1399  bool RejectEvent = kFALSE;
1400  int weightBinSign;
1401  if(MagField > 0) weightBinSign = 0;
1402  else weightBinSign = 1;
1403  float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1404                                                  {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1405  int weightBinCent = (int) CentPercent;
1406  if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
1407
1408  return RejectEvent;
1409 }