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