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