PWGCF/FEMTOSCOPY/K0Analysis/AliFemtoK0Analysis.cxx

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