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