PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx

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