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[5][5];
 463   TH3F *fHist3DOSLCutsBkg[5][5];
 464   if(fCutCheck){
 465    for(int i3D=0;i3D<5;i3D++){
 466     for(int j3D=0;j3D<5;j3D++){
 467      TString *histname = new TString("fHist3DOSLCuts");
 468      *histname += i3D;
 469      *histname += j3D;
 470      histname->Append("Signal");
 471      fHist3DOSLCutsSignal[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
 472      fOutputList->Add(fHist3DOSLCutsSignal[i3D][j3D]);
 473      histname->Replace(16,6,"Bkg");
 474      fHist3DOSLCutsBkg[i3D][j3D] = new TH3F(histname->Data(),"",100,-.5,.5,100,-.5,.5,100,-.5,.5);
 475      fOutputList->Add(fHist3DOSLCutsBkg[i3D][j3D]);
 476     }
 477    }
 478   }
 479
 480   //3D Spherical Harmonics
 481   //TH3F* fHistSHCentLowKt = new TH3F("fHistSHCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 482   //TH3F* fHistSHCentHighKt = new TH3F("fHistSHCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 483   //TH3F* fHistSHSemiCentLowKt = new TH3F("fHistSHSemiCentLowKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 484   //TH3F* fHistSHSemiCentHighKt = new TH3F("fHistSHSemiCentHighKt","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 485   //TH3F* fHistSHCentLowKtBkg = new TH3F("fHistSHCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 486   //TH3F* fHistSHCentHighKtBkg = new TH3F("fHistSHCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 487   //TH3F* fHistSHSemiCentLowKtBkg = new TH3F("fHistSHSemiCentLowKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 488   //TH3F* fHistSHSemiCentHighKtBkg = new TH3F("fHistSHSemiCentHighKtBkg","",50,0,.5,ncthetabins,-1,1,nphibins,0,2*PI);
 489   //fOutputList->Add(fHistSHCentLowKt);
 490   //fOutputList->Add(fHistSHCentHighKt);
 491   //fOutputList->Add(fHistSHSemiCentLowKt);
 492   //fOutputList->Add(fHistSHSemiCentHighKt);
 493   //fOutputList->Add(fHistSHCentLowKtBkg);
 494   //fOutputList->Add(fHistSHCentHighKtBkg);
 495   //fOutputList->Add(fHistSHSemiCentLowKtBkg);
 496   //fOutputList->Add(fHistSHSemiCentHighKtBkg);
 497
 498   //side-side
 499   //TH3F* fHistLeftLeftSignal = new TH3F("fHistLeftLeftSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 500   //TH3F* fHistLeftRightSignal = new TH3F("fHistLeftRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 501   //TH3F* fHistRightRightSignal = new TH3F("fHistRightRightSignal","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 502   //TH3F* fHistLeftLeftBkg = new TH3F("fHistLeftLeftBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 503   //TH3F* fHistLeftRightBkg = new TH3F("fHistLeftRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 504   //TH3F* fHistRightRightBkg = new TH3F("fHistRightRightBkg","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 505   //fOutputList->Add(fHistLeftLeftSignal);
 506   //fOutputList->Add(fHistLeftRightSignal);
 507   //fOutputList->Add(fHistRightRightSignal);
 508   //fOutputList->Add(fHistLeftLeftBkg);
 509   //fOutputList->Add(fHistLeftRightBkg);
 510   //fOutputList->Add(fHistRightRightBkg);
 511
 512   //TH3F* fHistSplitK0Sides = new TH3F("fHistSplitK0Sides","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 513   //fOutputList->Add(fHistSplitK0Sides);
 514   //TH3F* fHistSplitK0Centers = new TH3F("fHistSplitK0Centers","", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 515   //fOutputList->Add(fHistSplitK0Centers);
 516   //TH3F* fHistQinvSignalNoSplit = new TH3F("fHistQinvSignalNoSplit","Same Event Pair Distribution", kCentBins, .5, kCentBins+.5, 300, 0., 3., 100, 0., 1);
 517   //fOutputList->Add(fHistQinvSignalNoSplit);
 518
 519   PostData(1, fOutputList);
 520
 521 }
 522
 523 //________________________________________________________________________
 524 void AliFemtoK0Analysis::Exec(Option_t *)
 525 {
 526   // Main loop
 527   // Called for each event
 528   //cout<<"===========  Event # "<<fEventCount+1<<"  ==========="<<endl;
 529   fEventCount++;
 530   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
 531   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
 532
 533   Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral));
 534   bool isCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
 535   //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 536   if(!isSelected) {
 537    //cout << "Failed trigger selection." << endl;
 538    return;
 539   }
 540
 541   ///////////////////////////////////////////////////////////
 542
 543   unsigned int statusPos=0;
 544   unsigned int statusNeg=0;
 545
 546   float bField=0;
 547   bField = fAOD->GetMagneticField();
 548   if(bField == 0) return;
 549   if(fSignDep){
 550    if(fFieldPos && bField < 0) return;
 551    if(!fFieldPos && bField > 0) return;
 552   }
 553
 554
 555   /////////////////////////////////////////////////
 556
 557
 558   //Centrality selection
 559
 560   AliCentrality *centrality = fAOD->GetCentrality();
 561   float percent = centrality->GetCentralityPercentile("V0M");
 562   int centBin=0;
 563   //Printf("Centrality percent = %f", percent);
 564
 565   AliAODVZERO *aodV0 = fAOD->GetVZEROData();
 566   float multV0A=aodV0->GetMTotV0A();
 567   float multV0C=aodV0->GetMTotV0C();
 568
 569   if(percent < 0) {
 570    //Printf("No centrality info");
 571    return;
 572   }
 573   if(percent < 0.1 && (multV0A + multV0C < 19500)){
 574    //Printf("No centrality info");
 575    return;
 576   }
 577   else if(percent <= 5)   centBin=15;
 578   else if(percent <= 10)  centBin=14;
 579   else if(percent <= 15)  centBin=13;
 580   else if(percent <= 20)  centBin=12;
 581   else if(percent <= 25)  centBin=11;
 582   else if(percent <= 30)  centBin=10;
 583   else if(percent <= 35)  centBin=9;
 584   else if(percent <= 40)  centBin=8;
 585   else if(percent <= 45)  centBin=7;
 586   else if(percent <= 50)  centBin=6;
 587   else if(percent <= 55)  centBin=5;
 588   else if(percent <= 60)  centBin=4;
 589   else if(percent <= 65)  centBin=3;
 590   else if(percent <= 70)  centBin=2;
 591   else if(percent <= 75)  centBin=1;
 592   else if(percent <= 80)  centBin=0;
 593   else {
 594    //Printf("Skipping Peripheral Event");
 595    return;
 596   }
 597   if(percent > 10 && isCentral) return;
 598   if(fCutCheck && percent > 5) return;  //only looking at 0-5% for Cut Check
 599   ((TH1F*)fOutputList->FindObject("fHistCent"))->Fill(percent);
 600
 601   //flatten centrality dist.
 602   if(percent < 9){
 603    if(fFlatCent){
 604     if(RejectEventCentFlat(bField,percent)) return;
 605    }
 606   }
 607   ((TH1F*)fOutputList->FindObject("fHistCentFlat"))->Fill(percent);
 608
 609   //Vertexing
 610   AliAODVertex *primaryVertex;
 611   double vertex[3]={0};
 612   primaryVertex = fAOD->GetPrimaryVertex();
 613   vertex[0]=primaryVertex->GetX();
 614   vertex[1]=primaryVertex->GetY();
 615   vertex[2]=primaryVertex->GetZ();
 616   if(vertex[0]<10e-5 && vertex[1]<10e-5 &&  vertex[2]<10e-5) return;
 617   if(fabs(vertex[2]) > 10) return; // Z-vertex Cut
 618
 619   int zBin=0;
 620   double zStep=2*10/double(kZVertexBins), zStart=-10.;
 621   for(int i=0; i<kZVertexBins; i++)
 622   {
 623    if((vertex[2] > zStart+i*zStep) && (vertex[2] < zStart+(i+1)*zStep))
 624    {
 625     zBin=i;
 626     break;
 627    }
 628   }
 629
 630   //Event plane
 631   int psiBin = 0;
 632   AliEventplane *eventplane = fAOD->GetEventplane();
 633   if(fPsiBinning && !eventplane) return;
 634   double psiEP = eventplane->GetEventplane("V0",fAOD,2); //[-PI/2,PI/2]
 635   ((TH1F*)fOutputList->FindObject("fHistPsi"))->Fill(psiEP);
 636
 637   double psiStep = PI/double(fNPsiBins);
 638   double psiStart = -0.5*PI;
 639   for(int i=0; i<fNPsiBins; i++)
 640   {
 641    if((psiEP > psiStart+i*psiStep) && (psiEP < psiStart+(i+1)*psiStep))
 642    {
 643     psiBin = i;
 644     break;
 645    }
 646   }
 647   if(!fPsiBinning) psiBin = 0;
 648
 649 ////////////////////////////////////////////////////////////////
 650 //Cut Values and constants
 651
 652   //const bool kMCCase = kFALSE;                     //switch for MC analysis
 653   const int kMaxNumK0 = 300;                       //maximum number of K0s, array size
 654   const float kMinDCAPrimaryPion = 0.4;            //minimum dca of pions to primary
 655   const float kMaxDCADaughtersK0 = 0.3;            //maximum dca of pions to each other - 3D
 656   const float kMaxDCAK0 = 0.3;                     //maximum dca of K0 to primary
 657   const float kMaxDLK0 = 30.0;                     //maximum decay length of K0
 658   const float kMinDLK0 = fMinDecayLength;          //minimum decay length of K0
 659   const float kEtaCut = 0.8;                       //maximum |pseudorapidity|
 660   const float kMinCosAngle = 0.99;                 //minimum cosine of K0 pointing angle
 661
 662   const float kMinSeparation = fMinSep;                //minimum daughter (pair) separation
 663
 664   const float kTOFLow = 0.8;                       //boundary for TOF usage
 665   const float kMaxTOFSigmaPion = 3.0;              //TOF # of sigmas
 666   const float kMaxTPCSigmaPion = 3.0;              //TPC # of sigmas
 667
 668   //const float kMassPion = .13957;
 669   const float kMassK0Short = .497614;       //true PDG masses
 670
 671   //for cut checks
 672   double kCheckMassLow  [5] = {0.49,0.485,0.48,0.47,0.45};
 673   double kCheckMassHigh [5] = {0.505,0.510,0.515,0.525,0.550};
 674   double kCheckDCAK0    [5] = {0.1,0.2,0.3,0.5,1.0};
 675   double kCheckDCAPi    [5] = {1.0,0.6,0.4,0.2,0.1};
 676   double kCheckDCAPiPi  [5] = {0.1,0.2,0.3,0.5,1.0};
 677   double kCheckAvgSep   [5] = {0.0,2.5,5.0,7.5,10.0};
 678
 679 ////////////////////////////////////////////////////////////////
 680   //v0 tester
 681 ////////////////////////////////////////////////////////////////
 682   int v0Count = 0;      //number of v0s (entries in array)
 683   int k0Count = 0;      //number of good K0s
 684
 685   AliFemtoK0Particle *tempK0 = new AliFemtoK0Particle[kMultLimit];
 686
 687   //for daughter sharing studies
 688   //int idArray[100] = {0};
 689   //int idCount = 0;
 690
 691   //for MC
 692   //TClonesArray *mcArray = 0x0;
 693   //if(kMCCase){
 694   //mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
 695   //if(!mcArray){cout<<"No MC particle branch found"<<endl;return;}}
 696
 697   for(int i = 0; i < fAOD->GetNumberOfV0s(); i++)
 698   {
 699     bool goodK0 = kFALSE;
 700     bool goodPiPlus = kFALSE;
 701     bool goodPiMinus = kFALSE;
 702
 703     //load v0 track
 704     AliAODv0* v0 = fAOD->GetV0(i);
 705     if(!v0) continue;
 706     if(fOnlineCase){
 707      if(!(v0->GetOnFlyStatus())) continue;
 708     } //for online
 709     else{
 710      if((v0->GetOnFlyStatus())) continue; //for offline
 711     }
 712
 713     //for on-the-fly ordering
 714     AliAODTrack* tempTrack = (AliAODTrack*)v0->GetDaughter(0);
 715     short int pos0or1;
 716     short int neg0or1;
 717     bool orderswitch = kFALSE;
 718     if(tempTrack->Charge() > 0) {pos0or1 = 0; neg0or1 = 1;}
 719     else {pos0or1 = 1; neg0or1 = 0; orderswitch = kTRUE;}
 720
 721     //load daughter tracks
 722     AliAODTrack* prongTrackPos = (AliAODTrack*)v0->GetDaughter(pos0or1);
 723     AliAODTrack* prongTrackNeg = (AliAODTrack*)v0->GetDaughter(neg0or1);
 724     if(!prongTrackPos) continue;
 725     if(!prongTrackNeg) continue;
 726
 727     //daughter cuts
 728     if(v0->PtProng(pos0or1) < .15) continue;
 729     if(v0->PtProng(neg0or1) < .15) continue;
 730     if(fabs(v0->EtaProng(pos0or1)) > .8) continue;
 731     if(fabs(v0->EtaProng(neg0or1)) > .8) continue;
 732
 733     //load status for PID
 734     statusPos=prongTrackPos->GetStatus();
 735     if((statusPos&AliESDtrack::kTPCrefit)==0) continue;
 736     prongTrackPos->SetAODEvent(fAOD);
 737     statusNeg=prongTrackNeg->GetStatus();
 738     if((statusNeg&AliESDtrack::kTPCrefit)==0) continue;
 739     prongTrackNeg->SetAODEvent(fAOD);
 740
 741     //TPC PID
 742     if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackPos,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiPlus = kTRUE;
 743     if(fabs(fPidAOD->NumberOfSigmasTPC(prongTrackNeg,AliPID::kPion)) < kMaxTPCSigmaPion) goodPiMinus = kTRUE;
 744
 745     //Positive daughter identification TOF
 746     float probMis;
 747     AliPIDResponse::EDetPidStatus statusPosTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackPos);
 748     double Ppos = v0->PProng(pos0or1);
 749     if(Ppos > kTOFLow) //PiPlus
 750     {
 751      //if( (statusPos&AliESDtrack::kTOFpid)!=0 && (statusPos&AliESDtrack::kTIME)!=0 && (statusPos&AliESDtrack::kTOFout)!=0 && (statusPos&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
 752      if(AliPIDResponse::kDetPidOk == statusPosTOF)
 753      {
 754       probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
 755       if(probMis < 0.01) //avoid TOF-TPC mismatch
 756       {
 757        if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackPos,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiPlus = kTRUE;
 758        else goodPiPlus = kFALSE;
 759       }
 760      }
 761     }
 762     //Negative daughter identification TOF
 763     AliPIDResponse::EDetPidStatus statusNegTOF = fPidAOD->CheckPIDStatus(AliPIDResponse::kTOF, prongTrackNeg);
 764     double Pneg = v0->PProng(neg0or1);
 765     if(Pneg > kTOFLow) //PiMinus
 766     {
 767      //if( (statusNeg&AliESDtrack::kTOFpid)!=0 && (statusNeg&AliESDtrack::kTIME)!=0 && (statusNeg&AliESDtrack::kTOFout)!=0 && (statusNeg&AliESDtrack::kTOFmismatch)<=0) (OBSOLETE; NEW CALL BELOW)
 768      if(AliPIDResponse::kDetPidOk == statusNegTOF)
 769      {
 770       probMis = fPidAOD->GetTOFMismatchProbability(prongTrackPos);
 771       if(probMis < 0.01) //avoid TOF-TPC mismatch
 772       {
 773        if(fabs(fPidAOD->NumberOfSigmasTOF(prongTrackNeg,AliPID::kPion)) < kMaxTOFSigmaPion) goodPiMinus = kTRUE;
 774        else goodPiMinus = kFALSE;
 775       }
 776      }
 777     }
 778
 779     //K0 cuts
 780     if(!goodPiMinus || !goodPiPlus)                             continue;
 781     if(v0->Eta() > kEtaCut)                                     continue;
 782     if(v0->CosPointingAngle(primaryVertex) < kMinCosAngle)      continue;
 783     if(v0->MassK0Short() < .2 || v0->MassK0Short() > .8)        continue;
 784     if(v0->DecayLength(primaryVertex) > kMaxDLK0)               continue;
 785     if(v0->DecayLength(primaryVertex) < kMinDLK0)                       continue;
 786
 787     double v0Dca = v0->DcaV0ToPrimVertex();
 788     if(!fCutCheck){
 789      if(v0->DcaNegToPrimVertex() < kMinDCAPrimaryPion)      continue;
 790      if(v0->DcaPosToPrimVertex() < kMinDCAPrimaryPion)      continue;
 791      if(v0->DcaV0Daughters() > kMaxDCADaughtersK0)          continue;
 792      if(v0Dca > kMaxDCAK0)                                                      continue;
 793     }
 794     else{
 795      if(v0->DcaNegToPrimVertex() < kCheckDCAPi[4])                  continue;
 796      if(v0->DcaPosToPrimVertex() < kCheckDCAPi[4])              continue;
 797      if(v0->DcaV0Daughters() > kCheckDCAPiPi[4])                continue;
 798      if(v0Dca > kCheckDCAK0[4])                                         continue;
 799     }
 800
 801     //EVERYTHING BELOW HERE PASSES SINGLE PARTICLE CUTS, PION PID, and LOOSE MASS CUT
 802
 803     //for MC
 804     //bool MCgood = kFALSE;
 805     //if(kMCCase){
 806     //AliAODMCParticle* mck0dp = (AliAODMCParticle*)mcArray->At(abs(prongTrackPos->GetLabel()));
 807     //AliAODMCParticle* mck0dn = (AliAODMCParticle*)mcArray->At(abs(prongTrackNeg->GetLabel()));
 808     //if(mck0dp->GetMother() >= 0){
 809      //if(mck0dp->GetMother() == mck0dn->GetMother()){
 810       //if(abs(mck0dp->GetPdgCode()) == 211 && abs(mck0dn->GetPdgCode()) == 211){
 811        //AliAODMCParticle* mck0 = (AliAODMCParticle*)mcArray->At(mck0dp->GetMother());
 812        //if(abs(mck0->GetPdgCode()) == 310){
 813         //MCgood = kTRUE;
 814        //}
 815       //}
 816      //}
 817     //}
 818     //}// if kMCCase
 819
 820     if(!fCutCheck){
 821      if(v0->MassK0Short() > .48 && v0->MassK0Short() < .515) goodK0 = kTRUE;
 822     }
 823     else{
 824      if(v0->MassK0Short() > kCheckMassLow[4] && v0->MassK0Short() < kCheckMassHigh[4]) goodK0 = kTRUE;
 825     }
 826
 827     //Check for shared daughters, using v0 DCA to judge
 828     bool v0JudgeNew; //true if new v0 beats old
 829     tempK0[v0Count].fSkipShared = kFALSE;
 830     double newV0Pars[3] = {fabs(v0->MassK0Short()-kMassK0Short),v0Dca,v0->DcaV0Daughters()}; //parameters used in merit cut
 831     if(fMeritCase){
 832      for(int iID = 0; iID<v0Count; iID++){
 833       if(tempK0[iID].fSkipShared == kFALSE){            //if old is already skipped, go to next old
 834        if(tempK0[iID].fDaughterID1 == prongTrackPos->GetID() || tempK0[iID].fDaughterID2 == prongTrackNeg->GetID()){
 835         double oldV0Pars[3] = {fabs(tempK0[iID].fMass-kMassK0Short), tempK0[iID].fV0Dca, tempK0[iID].fDDDca};
 836         v0JudgeNew = CheckMeritCutWinner(fMeritCutChoice, oldV0Pars, newV0Pars); //true if new wins
 837         if(!v0JudgeNew){                //if old beats new...
 838          if(!tempK0[iID].fK0 && goodK0) continue;       //if bad old beats new good, do nothing...
 839          else{                                          //but if bad old beats new bad, or good old beats anything, skip new
 840           tempK0[v0Count].fSkipShared = kTRUE;          //skip new
 841           break;                                        //no need to keep checking others
 842          }
 843         }
 844         else{                                           //if new beats old...
 845          if(tempK0[iID].fK0 && !goodK0) continue;       //if bad new beats good old, do nothing...
 846          else{                                          //but if bad new beats bad old, or good new beats anything, skip old
 847               tempK0[iID].fSkipShared = kTRUE;          //skip old
 848               if(tempK0[iID].fK0) k0Count--;            //if good old gets skipped, subtract from number of K0s (new one will be added later, if it succeeds)
 849          }
 850         }
 851        }
 852       }
 853      }
 854      if(tempK0[v0Count].fSkipShared) continue;          //if new K0 is skipped, don't load; go to next v0
 855     }//if MeritCase
 856
 857     //for cut check
 858     if(fCutCheck){
 859      for(int iSet=0;iSet<4;iSet++){for(int jSet=0;jSet<4;jSet++){tempK0[v0Count].fCutPass[iSet][jSet] = kFALSE;}}
 860      for(int jCut = 0;jCut<5;jCut++){
 861       if(v0->MassK0Short() >    kCheckMassLow[jCut] && v0->MassK0Short() < kCheckMassHigh[jCut])
 862          tempK0[v0Count].fCutPass[0][jCut] = kTRUE;
 863       if(v0Dca < kCheckDCAK0[jCut]) tempK0[v0Count].fCutPass[1][jCut] = kTRUE;
 864       if(v0->DcaPosToPrimVertex() > kCheckDCAPi[jCut] && v0->DcaNegToPrimVertex() > kCheckDCAPi[jCut])
 865          tempK0[v0Count].fCutPass[2][jCut] = kTRUE;
 866       if(v0->DcaV0Daughters() < kCheckDCAPiPi[jCut]) tempK0[v0Count].fCutPass[3][jCut] = kTRUE;
 867      }
 868         }
 869
 870     //load parameters into temporary class instance
 871     if(v0Count < kMaxNumK0)
 872     {
 873         if(goodK0){
 874          tempK0[v0Count].fK0 = kTRUE;
 875          k0Count++;
 876         }
 877         else tempK0[v0Count].fK0 = kFALSE;
 878
 879         //if(v0->MassK0Short() > .45 && v0->MassK0Short() < .48) tempK0[v0Count].fSideLeft = kTRUE;
 880         //else tempK0[v0Count].fSideLeft = kFALSE;
 881         //if(v0->MassK0Short() > .515 && v0->MassK0Short() < .545) tempK0[v0Count].fSideRight = kTRUE;
 882         //else tempK0[v0Count].fSideRight = kFALSE;
 883             //if(!goodK0) continue; //no sides, speed up analysis (REDUNDANT RIGHT NOW)
 884
 885         tempK0[v0Count].fDaughterID1    = prongTrackPos->GetID();
 886         tempK0[v0Count].fDaughterID2    = prongTrackNeg->GetID();
 887         tempK0[v0Count].fMomentum[0]    = v0->Px();
 888         tempK0[v0Count].fMomentum[1]    = v0->Py();
 889         tempK0[v0Count].fMomentum[2]    = v0->Pz();
 890         tempK0[v0Count].fPt             = v0->Pt();
 891         tempK0[v0Count].fMass           = v0->MassK0Short();
 892         tempK0[v0Count].fV0Dca          = v0Dca;
 893
 894         //for hists
 895         tempK0[v0Count].fDDDca          = v0->DcaV0Daughters();
 896             tempK0[v0Count].fDecayLength        = v0->DecayLength(primaryVertex);
 897         tempK0[v0Count].fPosPt          = v0->PtProng(pos0or1);
 898         tempK0[v0Count].fNegPt          = v0->PtProng(neg0or1);
 899         tempK0[v0Count].fPosPhi         = v0->PhiProng(pos0or1);
 900         tempK0[v0Count].fNegPhi         = v0->PhiProng(neg0or1);
 901             if(!orderswitch){
 902          tempK0[v0Count].fPosDca        = v0->DcaPosToPrimVertex();
 903          tempK0[v0Count].fNegDca        = v0->DcaNegToPrimVertex();
 904             }
 905         else{
 906          tempK0[v0Count].fPosDca        = v0->DcaNegToPrimVertex();
 907          tempK0[v0Count].fNegDca        = v0->DcaPosToPrimVertex();
 908         }
 909
 910                 //for psi studies
 911         double v0Phi = v0->Phi(); //between [0,2pi]
 912         double v0PhiPsi = v0Phi-psiEP;
 913         if(v0PhiPsi < 0) v0PhiPsi += 2.*PI;
 914         else if (v0PhiPsi > 2.*PI) v0PhiPsi -= 2.*PI;
 915                 else{};
 916         tempK0[v0Count].fPhi    = v0Phi;
 917         tempK0[v0Count].fPhiPsi = v0PhiPsi;
 918
 919         //for separation
 920         GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackPos, bField, tempK0[v0Count].fPosXYZ, vertex);
 921             GetGlobalPositionAtGlobalRadiiThroughTPC(prongTrackNeg, bField, tempK0[v0Count].fNegXYZ, vertex);
 922         //for DPC
 923         prongTrackPos->GetPxPyPz(tempK0[v0Count].fPPos);
 924         prongTrackNeg->GetPxPyPz(tempK0[v0Count].fPNeg);
 925
 926         //if(idCount < 50){
 927         // if(goodK0){
 928         //  idArray[idCount*2]   = prongTrackPos->GetID();
 929         //  idArray[idCount*2+1] = prongTrackNeg->GetID();
 930         //  idCount++;
 931         //}}
 932
 933         v0Count++;
 934     }
 935
 936   }//v0
 937   if(k0Count<2) return;  //only keep events with more than 1 good K0
 938
 939   //Add Event to buffer - this is for event mixing
 940   fEC[zBin][centBin][psiBin]->FIFOShift();
 941   (fEvt) = fEC[zBin][centBin][psiBin]->fEvt;
 942   (fEvt)->fFillStatus = 1;
 943   int unskippedCount = 0;
 944   for(int i=0;i<v0Count;i++)
 945   {
 946    if(!tempK0[i].fSkipShared)                           //don't include skipped v0s (from shared daughters)
 947    {
 948     ((TH3F*)fOutputList->FindObject("fHistMass"))->Fill(centBin+1,tempK0[i].fPt,tempK0[i].fMass);
 949
 950     if(fCutCheck){
 951      for(int iCut=1;iCut<4;iCut++){
 952       for(int jCut=0;jCut<5;jCut++){
 953        TString *histname = new TString("fHistMassCuts");
 954        *histname += iCut;
 955        *histname += jCut;
 956        if(tempK0[i].fCutPass[iCut][jCut]) ((TH1F*)fOutputList->FindObject(histname->Data()))->Fill(tempK0[i].fMass);
 957       }
 958      }
 959     }
 960
 961     if(tempK0[i].fK0)                                   //make sure particle is good (mass)
 962     {
 963      (fEvt)->fK0Particle[unskippedCount] = tempK0[i];   //load good, unskipped K0s
 964      unskippedCount++;                                  //count good, unskipped K0s
 965     }
 966    }
 967   }
 968   (fEvt)->fNumV0s = unskippedCount;
 969   //Printf("Number of v0s: %d", v0Count);
 970   //Printf("Number of K0s: %d", k0Count);
 971   delete [] tempK0; tempK0 = NULL;
 972
 973   ((TH1F*)fOutputList->FindObject("fHistMultK0"))->Fill(unskippedCount);        // changed 3/25, used to be "k0Count"
 974   ((TH1F*)fOutputList->FindObject("fHistCentUsed"))->Fill(percent);
 975
 976   //Printf("Reconstruction Finished. Starting pair studies.");
 977
 978   //////////////////////////////////////////////////////////////////////
 979   // Correlations
 980   //////////////////////////////////////////////////////////////////////
 981
 982   float px1, py1, pz1, px2, py2, pz2;                   //single kaon values
 983   float en1, en2;                                                               //single kaon values
 984   //float pt1, pt2;                                                             //single kaon values
 985   float pairPx, pairPy, pairPz, pairP0;                 //pair momentum values
 986   float pairPt, pairMt, pairKt;                         //pair momentum values
 987   float pairMInv, pairPDotQ;
 988   float qinv, q0, qx, qy, qz;                           //pair q values
 989   //float qLength, thetaSH, thetaSHCos, phiSH;            //Spherical Harmonics values
 990   float am12, epm, h1, p12, p112, ppx, ppy, ppz, ks;    //PRF
 991   //float qOutLCMS;
 992   float qOutPRF, qSide, qLong;                          //relative momentum in LCMS/PRF frame
 993   float betasq, gamma;
 994   float p1LCMSOut, p1LCMSSide, p1LCMSLong, en1LCMS;
 995   float p2LCMSOut, p2LCMSSide, p2LCMSLong, en2LCMS;
 996
 997
 998   for(int i=0; i<(fEvt)->fNumV0s; i++) // Current event V0
 999   {
1000     //single particle histograms (done here to avoid "skipped" v0s
1001     ((TH1F*)fOutputList->FindObject("fHistDCADaughters"))       ->Fill((fEvt)->fK0Particle[i].fDDDca);
1002     ((TH1F*)fOutputList->FindObject("fHistDecayLengthK0"))      ->Fill((fEvt)->fK0Particle[i].fDecayLength);
1003     ((TH1F*)fOutputList->FindObject("fHistDCAK0"))              ->Fill((fEvt)->fK0Particle[i].fV0Dca);
1004     ((TH1F*)fOutputList->FindObject("fHistDCAPiMinus")) ->Fill((fEvt)->fK0Particle[i].fNegDca);
1005     ((TH1F*)fOutputList->FindObject("fHistDCAPiPlus"))          ->Fill((fEvt)->fK0Particle[i].fPosDca);
1006     ((TH2F*)fOutputList->FindObject("fHistPtK0"))               ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPt);
1007     ((TH2F*)fOutputList->FindObject("fHistK0PiPlusPt")) ->Fill(centBin+1, (fEvt)->fK0Particle[i].fPosPt);
1008     ((TH2F*)fOutputList->FindObject("fHistK0PiMinusPt"))        ->Fill(centBin+1, (fEvt)->fK0Particle[i].fNegPt);
1009     ((TH1F*)fOutputList->FindObject("fHistDaughterPhi"))        ->Fill((fEvt)->fK0Particle[i].fPosPhi);
1010     ((TH1F*)fOutputList->FindObject("fHistDaughterPhi"))        ->Fill((fEvt)->fK0Particle[i].fNegPhi);
1011
1012     ((TH1F*)fOutputList->FindObject("fHistPx"))         ->Fill((fEvt)->fK0Particle[i].fMomentum[0]);
1013     ((TH1F*)fOutputList->FindObject("fHistPy"))         ->Fill((fEvt)->fK0Particle[i].fMomentum[1]);
1014     ((TH1F*)fOutputList->FindObject("fHistPz"))         ->Fill((fEvt)->fK0Particle[i].fMomentum[2]);
1015
1016     ((TH2F*)fOutputList->FindObject("fHistPhi"))        ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhi);
1017     ((TH2F*)fOutputList->FindObject("fHistPhiPsi"))     ->Fill(centBin+1,(fEvt)->fK0Particle[i].fPhiPsi);
1018
1019     for(int evnum=0; evnum<kEventsToMix+1; evnum++)// Event buffer loop: evnum=0 is the current event, all other evnum's are past events
1020     {
1021       int startbin=0;
1022       if(evnum==0) startbin=i+1;
1023
1024       for(int j=startbin; j<(fEvt+evnum)->fNumV0s; j++) // Past event V0
1025       {
1026         if(evnum==0)  // Get rid of shared tracks
1027         {
1028           if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
1029                   if((fEvt)->fK0Particle[i].fDaughterID1 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
1030               if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID1) continue;
1031               if((fEvt)->fK0Particle[i].fDaughterID2 == (fEvt+evnum)->fK0Particle[j].fDaughterID2) continue;
1032             }
1033
1034         px1 = (fEvt)->fK0Particle[i].fMomentum[0];
1035                 py1 = (fEvt)->fK0Particle[i].fMomentum[1];
1036                 pz1 = (fEvt)->fK0Particle[i].fMomentum[2];
1037         //pt1 = (fEvt)->fK0Particle[i].fPt;
1038                 px2 = (fEvt+evnum)->fK0Particle[j].fMomentum[0];
1039                 py2 = (fEvt+evnum)->fK0Particle[j].fMomentum[1];
1040                 pz2 = (fEvt+evnum)->fK0Particle[j].fMomentum[2];
1041         //pt2 = (fEvt+evnum)->fK0Particle[j].fPt;
1042         if(fRandomNumber->Rndm() < .5){ //switch particle order for 3D qout bias
1043                  double tempvar;
1044          tempvar = px1; px1 = px2; px2 = tempvar;
1045          tempvar = py1; py1 = py2; py2 = tempvar;
1046          tempvar = pz1; pz1 = pz2; pz2 = tempvar;
1047                 }
1048
1049                 en1  = sqrt(pow(px1,2)+pow(py1,2)+pow(pz1,2)+pow(kMassK0Short,2));
1050                 en2  = sqrt(pow(px2,2)+pow(py2,2)+pow(pz2,2)+pow(kMassK0Short,2));
1051
1052         q0 = en1 - en2;
1053         qx = px1 - px2;
1054         qy = py1 - py2;
1055         qz = pz1 - pz2;
1056         qinv = sqrt(pow(qx,2) + pow(qy,2) + pow(qz,2) - pow(q0,2));
1057
1058         pairPx = px1 + px2;
1059                 pairPy = py1 + py2;
1060                 pairPz = pz1 + pz2;
1061         pairP0 = en1 + en2;
1062                 pairPt = sqrt(pairPx*pairPx + pairPy*pairPy);
1063         pairKt = pairPt/2.;                                                                     //used for KT binning
1064         pairMt = sqrt(pairP0*pairP0 - pairPz*pairPz);           //used for LCMS (not plots)
1065                 pairMInv = sqrt(pow(pairP0,2)-pow(pairPx,2)-pow(pairPy,2)-pow(pairPz,2));//used for PRF
1066         pairPDotQ = pairP0*q0-pairPx*qx-pairPy*qy-pairPz*qz;    //used for PRF
1067
1068             //PRF (this section will probably be removed in favor of later boosting section)
1069         p12 = sqrt(pow(pairPx,2)+pow(pairPy,2)+pow(pairPz,2));  //pair momentum length
1070         am12 = sqrt(pow(en1+en2,2)-p12*p12);                                    //sqrt(s)=|p1+p2|(4vec)
1071         epm = en1+en2+am12;                                                                             //"energy plus mass"
1072         p112 = px1*pairPx+py1*pairPy+pz1*pairPz;                                //proj. of p1 on pairP
1073         if(am12 == 0) continue;
1074         h1 = (p112/epm - en1)/am12;
1075         ppx = px1+pairPx*h1;                                                                    //px in PRF
1076         ppy = py1+pairPy*h1;                                                                    //py in PRF
1077         ppz = pz1+pairPz*h1;                                                                    //pz in PRF
1078         ks = sqrt(ppx*ppx+ppy*ppy+ppz*ppz);                                             //k*
1079         ((TH1F*)fOutputList->FindObject("fHistPxCM"))->Fill(ppx);
1080         ((TH1F*)fOutputList->FindObject("fHistPyCM"))->Fill(ppy);
1081         ((TH1F*)fOutputList->FindObject("fHistPzCM"))->Fill(ppz);
1082         ((TH1F*)fOutputList->FindObject("fHistKsCM"))->Fill(ks);
1083
1084         //relative momentum in out-side-long for LCMS and PRF
1085         if(pairMt == 0 || pairPt == 0) continue;
1086         qLong = (pairP0*qz - pairPz*q0)/pairMt; //same for both frames
1087         qSide = (pairPx*qy - pairPy*qx)/pairPt; //same for both frames
1088         //qOutLCMS = (pairPx*qx + pairPy*qy)/pairPt;
1089                 qOutPRF  = pairMInv*(pairPx*qx+pairPy*qy)/pairMt/pairPt - pairPt*pairPDotQ/pairMt/pairMInv;
1090
1091                 //finding gamma for gamma binning/hists (likely will be removed after tests)
1092                 p1LCMSOut  = (pairPx*px1+pairPy*py1)/pairPt;
1093                 p1LCMSSide = (pairPx*py1-pairPy*px1)/pairPt;
1094                 p1LCMSLong = (pairP0*pz1-pairPz*en1)/pairMt;
1095                 p2LCMSOut  = (pairPx*px2+pairPy*py2)/pairPt;
1096                 p2LCMSSide = (pairPx*py2-pairPy*px2)/pairPt;
1097                 p2LCMSLong = (pairP0*pz2-pairPz*en2)/pairMt;
1098                 en1LCMS = sqrt(pow(p1LCMSOut,2)+pow(p1LCMSSide,2)+pow(p1LCMSLong,2)+pow(kMassK0Short,2));
1099                 en2LCMS = sqrt(pow(p2LCMSOut,2)+pow(p2LCMSSide,2)+pow(p2LCMSLong,2)+pow(kMassK0Short,2));
1100                 betasq = pow((p1LCMSOut+p2LCMSOut)/(en1LCMS+en2LCMS),2);
1101                 gamma = 1./sqrt(1-betasq);
1102                 ((TH2F*)fOutputList->FindObject("fHistGamma"))->Fill(gamma,qinv);
1103                 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p1LCMSOut);
1104                 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p1LCMSSide);
1105                 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p1LCMSLong);
1106                 ((TH1F*)fOutputList->FindObject("fHistPOutLCMS"))->Fill(p2LCMSOut);
1107                 ((TH1F*)fOutputList->FindObject("fHistPSideLCMS"))->Fill(p2LCMSSide);
1108                 ((TH1F*)fOutputList->FindObject("fHistPLongLCMS"))->Fill(p2LCMSLong);
1109                 //getting bin numbers and names for 3D histogram
1110         TString *histname3D = new TString("fHist3DOSL");
1111         int ktBin;
1112         if(pairKt < 0.6) ktBin = 0;
1113                 else if(pairKt < 0.8) ktBin = 1;
1114                 else if(pairKt < 1.0) ktBin = 2;
1115                 else ktBin = 3;
1116                 *histname3D += centBin-6; //centBins: [6,15] -> array bins: [0,9]
1117                 *histname3D += ktBin;
1118
1119         //Spherical harmonics
1120         //qLength = sqrt(qLong*qLong + qSide*qSide + qOutPRF*qOutPRF);
1121         //thetaSHCos = qLong/qLength;
1122         //thetaSH = acos(thetaSHCos);
1123         //phiSH = acos(qOutPRF/(qLength*sin(thetaSH)));
1124
1125         //Finding average separation of daughters throughout TPC - two-track cut
1126         float posPositions1[9][3] = {{0}};
1127         float negPositions1[9][3] = {{0}};
1128         float posPositions2[9][3] = {{0}};
1129         float negPositions2[9][3] = {{0}};
1130         for(int iPos = 0; iPos < 9; iPos++){
1131          for(int jPos = 0; jPos < 3; jPos++){
1132            posPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fPosXYZ[iPos][jPos];
1133            negPositions1[iPos][jPos] = (fEvt)->fK0Particle[i].fNegXYZ[iPos][jPos];
1134            posPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fPosXYZ[iPos][jPos];
1135            negPositions2[iPos][jPos] = (fEvt+evnum)->fK0Particle[j].fNegXYZ[iPos][jPos];
1136          }
1137         }
1138         float pMean = 0.;       //average separation for positive daughters
1139         float nMean = 0.;       //average separation for negative daughters
1140         float pDiff;
1141         float nDiff;
1142         float pMin = 9999.;     //minimum separation (updates) - pos
1143         float nMin = 9999.;     //minimum separation (updates) - neg
1144         double pCount=0;        //counter for number of points used - pos
1145         double nCount=0;        //counter for number of points used - neg
1146         for(int ss=0;ss<9;ss++){
1147          if(posPositions1[ss][0] != -9999 && posPositions2[ss][0] != -9999){
1148           pCount++;
1149           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));
1150           pMean = pMean + pDiff;
1151           if(pDiff < pMin) pMin = pDiff;
1152          }
1153          if(negPositions1[ss][0] != -9999 && negPositions1[ss][0] != -9999){
1154           nCount++;
1155           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));
1156           nMean = nMean + nDiff;
1157           if(nDiff < nMin) nMin = nDiff;
1158          }
1159         }
1160         pMean = pMean/pCount;
1161         nMean = nMean/nCount;
1162
1163         if(evnum==0){
1164          ((TH1F*)fOutputList->FindObject("fHistSepNumPos"))->Fill(pMean);
1165          ((TH1F*)fOutputList->FindObject("fHistSepNumNeg"))->Fill(nMean);
1166          ((TH2F*)fOutputList->FindObject("fHistSepNumPos2"))->Fill(pMean,pMin);
1167          ((TH2F*)fOutputList->FindObject("fHistSepNumNeg2"))->Fill(nMean,nMin);
1168         }
1169         else{
1170          ((TH1F*)fOutputList->FindObject("fHistSepDenPos"))->Fill(pMean);
1171          ((TH1F*)fOutputList->FindObject("fHistSepDenNeg"))->Fill(nMean);
1172          ((TH2F*)fOutputList->FindObject("fHistSepDenPos2"))->Fill(pMean,pMin);
1173          ((TH2F*)fOutputList->FindObject("fHistSepDenNeg2"))->Fill(nMean,nMin);
1174         }
1175
1176         //Decay plane coincidence
1177         //daughter momenta
1178         float a1 = (fEvt)->fK0Particle[i].fPPos[0];
1179         float b1 = (fEvt)->fK0Particle[i].fPPos[1];
1180         float c1 = (fEvt)->fK0Particle[i].fPPos[2];
1181         float d1 = (fEvt)->fK0Particle[i].fPNeg[0];
1182         float e1 = (fEvt)->fK0Particle[i].fPNeg[1];
1183         float f1 = (fEvt)->fK0Particle[i].fPNeg[2];
1184         float a2 = (fEvt+evnum)->fK0Particle[j].fPPos[0];
1185         float b2 = (fEvt+evnum)->fK0Particle[j].fPPos[1];
1186         float c2 = (fEvt+evnum)->fK0Particle[j].fPPos[2];
1187         float d2 = (fEvt+evnum)->fK0Particle[j].fPNeg[0];
1188         float e2 = (fEvt+evnum)->fK0Particle[j].fPNeg[1];
1189         float f2 = (fEvt+evnum)->fK0Particle[j].fPNeg[2];
1190
1191         float cross1[3];
1192         float cross2[3];
1193         cross1[0] = b1*f1-c1*e1;
1194         cross1[1] = c1*d1-a1*f1;
1195         cross1[2] = a1*e1-b1*d1;
1196         cross2[0] = b2*f2-c2*e2;
1197         cross2[1] = c2*d2-a2*f2;
1198         cross2[2] = a2*e2-b2*d2;
1199         float crosslength1 = sqrt(pow(cross1[0],2)+pow(cross1[1],2)+pow(cross1[2],2));
1200         float crosslength2 = sqrt(pow(cross2[0],2)+pow(cross2[1],2)+pow(cross2[2],2));
1201         float dpc = (cross1[0]*cross2[0]+cross1[1]*cross2[1]+cross1[2]*cross2[2])/(crosslength1*crosslength2);
1202
1203         if(evnum==0)((TH2F*)fOutputList->FindObject("fHistSepDPC"))->Fill(dpc,pMean);
1204         else ((TH2F*)fOutputList->FindObject("fHistSepDPCBkg"))->Fill(dpc,pMean);
1205
1206         bool SepPass[5] = {0};
1207         if(!fCutCheck){
1208          if(pMean < kMinSeparation || nMean < kMinSeparation) continue; //using the "new" method (ala Hans)
1209         }
1210         else{
1211          if(pMean < kCheckAvgSep[4] || nMean < kCheckAvgSep[4]) continue;
1212          for(int jCut=0;jCut<5;jCut++){
1213           if(pMean > kCheckAvgSep[4] && nMean > kCheckAvgSep[4]) SepPass[jCut] = kTRUE;
1214          }
1215         }
1216
1217         //end separation studies
1218
1219         //Fill Histograms
1220         bool center1K0   = kFALSE;  //accepted mass K0
1221             bool center2K0   = kFALSE;
1222         if((fEvt)->fK0Particle[i].fK0) center1K0=kTRUE;
1223         if((fEvt+evnum)->fK0Particle[j].fK0) center2K0=kTRUE;
1224         //bool CL1 = kFALSE;
1225         //bool CL2 = kFALSE;
1226         //bool CR1 = kFALSE;
1227         //bool CR2 = kFALSE;
1228         //if(center1K0 && center2K0){
1229         // if((fEvt)->fK0Particle[i].fMass < kMassK0Short) CL1 = kTRUE;
1230         // else CR1 = kTRUE;
1231         // if((fEvt+evnum)->fK0Particle[j].fMass < kMassK0Short) CL2 = kTRUE;
1232         // else CR2 = kTRUE;
1233         //}
1234
1235         //bool SideLeft1 = kFALSE;
1236         //bool SideLeft2 = kFALSE;
1237         //bool SideRight1 = kFALSE;
1238         //bool SideRight2 = kFALSE;
1239         //if((fEvt)->fK0Particle[i].fSideLeft) SideLeft1 = kTRUE;
1240         //else if((fEvt)->fK0Particle[i].fSideRight) SideRight1 = kTRUE;
1241         //if((fEvt+evnum)->fK0Particle[j].fSideLeft) SideLeft2 = kTRUE;
1242         //else if((fEvt+evnum)->fK0Particle[j].fSideRight) SideRight2 = kTRUE;
1243
1244         //for psi binning
1245                 float phipsi1 = (fEvt)->fK0Particle[i].fPhiPsi;
1246         float phipsi2 = (fEvt+evnum)->fK0Particle[j].fPhiPsi;
1247         bool inPlane1 = kFALSE;
1248         bool inPlane2 = kFALSE;
1249         if(phipsi1 > PI) phipsi1 = phipsi1-PI;
1250         if(phipsi2 > PI) phipsi2 = phipsi2-PI;
1251         if(phipsi1 < 0.25*PI || phipsi1 > 0.75*PI) inPlane1 = kTRUE;
1252         if(phipsi2 < 0.25*PI || phipsi2 > 0.75*PI) inPlane2 = kTRUE;
1253
1254                 //for dphi, dphipsi
1255                 float dPhi = fabs((fEvt)->fK0Particle[i].fPhi - (fEvt+evnum)->fK0Particle[j].fPhi);
1256         if(dPhi > PI) dPhi = 2*PI-dPhi;
1257         float dPhiPsi = fabs((fEvt)->fK0Particle[i].fPhiPsi - (fEvt+evnum)->fK0Particle[j].fPhiPsi);
1258         if(dPhiPsi > PI) dPhiPsi = 2*PI-dPhiPsi;
1259
1260         if(evnum==0) //Same Event
1261         {
1262           //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt)->fK0Particle[i].fMass);
1263           //((TH3F*)fOutputList->FindObject("fHistMassQKt"))->Fill(qinv, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1264           //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1265           //((TH3F*)fOutputList->FindObject("fHistMassKtK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1266           //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1267           //((TH3F*)fOutputList->FindObject("fHistMassPtCFK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1268
1269           if(center1K0 && center2K0){
1270            //1D
1271            ((TH3F*)fOutputList->FindObject("fHistQinvSignal"))->Fill(centBin+1, pairKt, qinv);
1272            //if(!splitK0centers)((TH3F*)fOutputList->FindObject("fHistQinvSignalNoSplit"))->Fill(centBin+1, pairKt, qinv);
1273            ((TH2F*)fOutputList->FindObject("fHistKtK0"))->Fill(centBin+1, pairKt);
1274
1275                    //eventplane
1276            if(inPlane1 && inPlane2)
1277              ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPIn"))->Fill(centBin+1, pairKt, qinv);
1278            else if(!inPlane1 && !inPlane2)
1279              ((TH3F*)fOutputList->FindObject("fHistQinvSignalEPOut"))->Fill(centBin+1, pairKt, qinv);
1280
1281            //dPhi,dPhiPsi
1282                    ((TH3F*)fOutputList->FindObject("fHistDPhi"))->Fill(centBin+1,pairKt,dPhi);
1283                    ((TH3F*)fOutputList->FindObject("fHistDPhiPsi"))->Fill(centBin+1,pairKt,dPhiPsi);
1284
1285
1286            //for mass bin study
1287            //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLSignal"))->Fill(centBin+1, pairKt, qinv);
1288            //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRSignal"))->Fill(centBin+1, pairKt, qinv);
1289            //else ((TH3F*)fOutputList->FindObject("fHistCRCRSignal"))->Fill(centBin+1, pairKt, qinv);
1290
1291            //3D
1292            if(fCase3D){
1293             if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1294                      histname3D->Append("Signal");
1295                          ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1296                     }
1297            }
1298
1299            //for cut check (3D, but fCase3D set to false)
1300            if(fCutCheck){
1301             for(int iCut=0;iCut<4;iCut++){//different cuts (4 + AvgSep)
1302              bool Skip = kFALSE;
1303              for(int iCut2=0;iCut2<4;iCut2++){//for setting other cuts to "3"
1304               if(iCut2 != iCut){
1305                if(!(fEvt)->fK0Particle[i].fCutPass[iCut2][2] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut2][2]) Skip = kTRUE;
1306               }
1307              }
1308              if(!SepPass[2]) Skip = kTRUE;
1309              if(Skip) continue;
1310              for(int jCut=0;jCut<5;jCut++){//different cut values
1311               TString *histname = new TString("fHist3DOSLCuts");
1312               *histname += iCut;
1313               *histname += jCut;
1314               histname->Append("Signal");
1315               if((fEvt)->fK0Particle[i].fCutPass[iCut][jCut] && (fEvt+evnum)->fK0Particle[j].fCutPass[iCut][jCut])
1316                ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1317              }//jcut
1318             }//icut
1319
1320             //for avg sep
1321             bool asSkip = kFALSE;
1322             for(int iCut=0;iCut<4;iCut++){
1323              if(!(fEvt)->fK0Particle[i].fCutPass[iCut][2] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut][2]) asSkip=kTRUE;
1324             }
1325             if(asSkip) continue;
1326             for(int jCut=0;jCut<5;jCut++){
1327              TString *histname = new TString("fHist3DOSLCuts");
1328              *histname += 4; //4 for AvgSep
1329              *histname += jCut;
1330              histname->Append("Signal");
1331              if(SepPass[jCut]) ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1332             }
1333            }//cutCheck
1334
1335
1336            /*if(pairKt < 1.0){
1337             if(centBin > 13){
1338              ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1339              ((TH3F*)fOutputList->FindObject("fHistSHCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}
1340             else if(centBin > 9){
1341              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKt"))->Fill(qOutPRF,qSide,qLong);
1342              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKt"))->Fill(qLength,thetaSHCos,phiSH);}}
1343            else if(pairKt < 2.0){
1344             if(centBin > 13){
1345              ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1346              ((TH3F*)fOutputList->FindObject("fHistSHCentHighKt"))->Fill(qLength,thetaSHCos, phiSH);}
1347             else if(centBin > 9){
1348              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKt"))->Fill(qOutPRF,qSide,qLong);
1349
1350              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKt"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1351
1352           }//centercenter
1353
1354           //side-side correlations
1355           //if(!splitK0sides){
1356           // if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftSignal"))->Fill(centBin+1, pairKt, qinv);
1357            //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightSignal"))->Fill(centBin+1, pairKt, qinv);
1358            //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightSignal"))->Fill(centBin+1, pairKt, qinv);
1359          //}
1360
1361         }//same event
1362
1363         else //Mixed Events
1364         {
1365           //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt)->fK0Particle[i].fMass);
1366           //((TH3F*)fOutputList->FindObject("fHistMassKtBkgK0"))->Fill(centBin+1, pairKt, (fEvt+evnum)->fK0Particle[j].fMass);
1367           //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt1, (fEvt)->fK0Particle[i].fMass);
1368           //((TH3F*)fOutputList->FindObject("fHistMassPtCFBkgK0"))->Fill(centBin+1, pt2, (fEvt+evnum)->fK0Particle[j].fMass);
1369
1370           if(center1K0 && center2K0){
1371            //1D
1372            ((TH3F*)fOutputList->FindObject("fHistQinvBkg"))->Fill(centBin+1, pairKt, qinv);
1373
1374                    //eventplane
1375            if(inPlane1 && inPlane2)
1376              ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPIn"))->Fill(centBin+1, pairKt, qinv);
1377            else if(!inPlane1 && !inPlane2)
1378              ((TH3F*)fOutputList->FindObject("fHistQinvBkgEPOut"))->Fill(centBin+1, pairKt, qinv);
1379
1380                    ((TH3F*)fOutputList->FindObject("fHistDPhiBkg"))->Fill(centBin+1,pairKt,dPhi);
1381                    ((TH3F*)fOutputList->FindObject("fHistDPhiPsiBkg"))->Fill(centBin+1,pairKt,dPhiPsi);
1382
1383            //for mass bin study
1384            //if(CL1 && CL2) ((TH3F*)fOutputList->FindObject("fHistCLCLBkg"))->Fill(centBin+1, pairKt, qinv);
1385            //else if ((CL1 && CR2) || (CR1 && CL2)) ((TH3F*)fOutputList->FindObject("fHistCLCRBkg"))->Fill(centBin+1, pairKt, qinv);
1386            //else ((TH3F*)fOutputList->FindObject("fHistCRCRBkg"))->Fill(centBin+1, pairKt, qinv);
1387
1388            //3D
1389            if(fCase3D){
1390             if(pairKt > 0.2 && pairKt < 1.5 && centBin > 5){
1391                      histname3D->Replace(12,6,"Bkg");
1392                          ((TH3F*)fOutputList->FindObject(histname3D->Data()))->Fill(qOutPRF,qSide,qLong);
1393                     }
1394            }
1395
1396            //for cut check (3D, but fCase3D set to false)
1397            if(fCutCheck){
1398             for(int iCut=0;iCut<4;iCut++){//different cuts (4 + AvgSep)
1399              bool Skip = kFALSE;
1400              for(int iCut2=0;iCut2<4;iCut2++){//for setting other cuts to "3"
1401               if(iCut2 != iCut){
1402                if(!(fEvt)->fK0Particle[i].fCutPass[iCut2][2] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut2][2]) Skip = kTRUE;
1403               }
1404              }
1405              if(!SepPass[2]) Skip = kTRUE;
1406              if(Skip) continue;
1407              for(int jCut=0;jCut<5;jCut++){//different cut values
1408               TString *histname = new TString("fHist3DOSLCuts");
1409               *histname += iCut;
1410               *histname += jCut;
1411               histname->Append("Bkg");
1412               if((fEvt)->fK0Particle[i].fCutPass[iCut][jCut] && (fEvt+evnum)->fK0Particle[j].fCutPass[iCut][jCut])
1413                ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1414              }//jcut
1415             }//icut
1416
1417             //for avg sep
1418             bool asSkip = kFALSE;
1419             for(int iCut=0;iCut<4;iCut++){
1420              if(!(fEvt)->fK0Particle[i].fCutPass[iCut][2] || !(fEvt+evnum)->fK0Particle[j].fCutPass[iCut][2]) asSkip=kTRUE;
1421             }
1422             if(asSkip) continue;
1423             for(int jCut=0;jCut<5;jCut++){
1424              TString *histname = new TString("fHist3DOSLCuts");
1425              *histname += 4; //4 for AvgSep
1426              *histname += jCut;
1427              histname->Append("Bkg");
1428              if(SepPass[jCut]) ((TH3F*)fOutputList->FindObject(histname->Data()))->Fill(qOutPRF,qSide,qLong);
1429             }
1430            }//cutCheck
1431
1432
1433            /*if(pairKt < 1.0){
1434             if(centBin > 13){
1435              ((TH3F*)fOutputList->FindObject("fHistOSLCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1436              ((TH3F*)fOutputList->FindObject("fHistSHCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}
1437             else if(centBin > 9){
1438              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentLowKtBkg"))->Fill(qOutPRF,qSide,qLong);
1439              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentLowKtBkg"))->Fill(qLength,thetaSHCos,phiSH);}}
1440            else if(pairKt < 2.0){
1441             if(centBin > 13){
1442              ((TH3F*)fOutputList->FindObject("fHistOSLCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1443              ((TH3F*)fOutputList->FindObject("fHistSHCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}
1444             else if(centBin > 9){
1445              ((TH3F*)fOutputList->FindObject("fHistOSLSemiCentHighKtBkg"))->Fill(qOutPRF,qSide,qLong);
1446              ((TH3F*)fOutputList->FindObject("fHistSHSemiCentHighKtBkg"))->Fill(qLength, thetaSHCos, phiSH);}}*/
1447           }
1448
1449           //side-side correlations
1450           //if(SideLeft1 && SideLeft2) ((TH3F*)fOutputList->FindObject("fHistLeftLeftBkg"))->Fill(centBin+1, pairKt, qinv);
1451           //else if((SideLeft1 && SideRight2) || (SideRight1 && SideLeft2)) ((TH3F*)fOutputList->FindObject("fHistLeftRightBkg"))->Fill(centBin+1, pairKt, qinv);
1452           //else if(SideRight1 && SideRight2) ((TH3F*)fOutputList->FindObject("fHistRightRightBkg"))->Fill(centBin+1, pairKt, qinv);
1453
1454         }//Mixed Events
1455
1456       }//past event
1457     }//event buffer
1458   }//current event
1459
1460   // Post output data.
1461   PostData(1, fOutputList);
1462
1463   }
1464 //________________________________________________________________________
1465 void AliFemtoK0Analysis::Terminate(Option_t *)
1466 {
1467   // Called once at the end of the query
1468   cout<<"Done"<<endl;
1469
1470 }
1471
1472 //_________________________________________________________________________
1473 void AliFemtoK0Analysis::GetGlobalPositionAtGlobalRadiiThroughTPC(const AliAODTrack *track, const Float_t bfield, Float_t globalPositionsAtRadii[9][3], double PrimaryVertex[3]){
1474   // Gets the global position of the track at nine different radii in the TPC
1475   // track is the track you want to propagate
1476   // bfield is the magnetic field of your event
1477   // globalPositionsAtRadii is the array of global positions in the radii and xyz
1478
1479   // Initialize the array to something indicating there was no propagation
1480   for(Int_t i=0;i<9;i++){
1481     for(Int_t j=0;j<3;j++){
1482       globalPositionsAtRadii[i][j]=-9999.;
1483     }
1484   }
1485
1486    // Make a copy of the track to not change parameters of the track
1487   AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1488   //printf("\nAfter CopyFromVTrack\n");
1489   //etp.Print();
1490
1491   // The global position of the the track
1492   Double_t xyz[3]={-9999.,-9999.,-9999.};
1493
1494   // Counter for which radius we want
1495   Int_t iR=0;
1496   // The radii at which we get the global positions
1497   // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1498   Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1499   // The global radius we are at
1500   Float_t globalRadius=0;
1501
1502   // Propagation is done in local x of the track
1503   for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1504     // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1505     // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1506     // If the track's momentum is smaller than infinite, it will develop a y-component, which
1507     // adds to the global radius
1508
1509     // Stop if the propagation was not succesful. This can happen for low pt tracks
1510     // that don't reach outer radii
1511     if(!etp.PropagateTo(x,bfield))break;
1512     etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1513     globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1514
1515     // Roughly reached the radius we want
1516     if(globalRadius > Rwanted[iR]){
1517
1518       // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1519       while (globalRadius>Rwanted[iR]){
1520         x-=.1;
1521         //      printf("propagating to x %5.2f\n",x);
1522         if(!etp.PropagateTo(x,bfield))break;
1523         etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1524         globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1525       }
1526       //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]);
1527       globalPositionsAtRadii[iR][0]=xyz[0];
1528       globalPositionsAtRadii[iR][1]=xyz[1];
1529       globalPositionsAtRadii[iR][2]=xyz[2];
1530       //subtract primary vertex, "zero" track for correct mixed-event comparison
1531       globalPositionsAtRadii[iR][0] -= PrimaryVertex[0];
1532       globalPositionsAtRadii[iR][1] -= PrimaryVertex[1];
1533       globalPositionsAtRadii[iR][2] -= PrimaryVertex[2];
1534
1535       // Indicate we want the next radius
1536       iR+=1;
1537     }
1538     if(iR>=8){
1539       // TPC edge reached
1540       return;
1541     }
1542   }
1543 }
1544
1545 bool AliFemtoK0Analysis::CheckMeritCutWinner(int cutChoice, double oldPars[3], double newPars[3]){
1546  //performs "merit cut" judgement check on v0s with shared daughters, using one of three criteria.
1547  //if cutChoice = 4, it uses all three criteria, needed 2 of 3 'points'
1548
1549  bool newV0Wins = kFALSE;
1550  double pardiff[3] = {newPars[0]-oldPars[0],
1551                       newPars[1]-oldPars[1],
1552                       newPars[2]-oldPars[2]};
1553  if(cutChoice > 0 && cutChoice < 4){
1554   if(pardiff[cutChoice] <= 0.) newV0Wins = kTRUE;
1555  }
1556  else if(cutChoice == 4){
1557   int newWinCount = 0;
1558   for(int i=0;i<3;i++){if(pardiff[i+1] <= 0) newWinCount++;}
1559   if(newWinCount > 1) newV0Wins = kTRUE;
1560  }
1561  else{};
1562  return newV0Wins;
1563 }
1564
1565 bool AliFemtoK0Analysis::RejectEventCentFlat(float MagField, float CentPercent)
1566 { // to flatten centrality distribution
1567  bool RejectEvent = kFALSE;
1568  int weightBinSign;
1569  if(MagField > 0) weightBinSign = 0;
1570  else weightBinSign = 1;
1571  float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1572                                                  {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1573  int weightBinCent = (int) CentPercent;
1574  if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
1575
1576  return RejectEvent;
1577 }