PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx

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