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