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