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