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