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