K0s code update (Matt Steinpreis)
[u/mrichter/AliRoot.git] / 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 }