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