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