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