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