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