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