]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/Chaoticity/AliChaoticity.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
r3 methodology change
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / Chaoticity / AliChaoticity.cxx
CommitLineData
cd12341d 1#include <iostream>
2#include <math.h>
3#include "TChain.h"
4#include "TFile.h"
5#include "TKey.h"
6#include "TObject.h"
7#include "TObjString.h"
8#include "TList.h"
9#include "TTree.h"
10#include "TH1F.h"
11#include "TH1D.h"
12#include "TH2D.h"
13#include "TH3D.h"
14#include "TProfile.h"
15#include "TProfile2D.h"
16#include "TCanvas.h"
17#include "TRandom3.h"
18#include "TF1.h"
19
20#include "AliAnalysisTask.h"
21#include "AliAnalysisManager.h"
22
23
24#include "AliESDEvent.h"
25#include "AliESDInputHandler.h"
26#include "AliESDtrackCuts.h"
27
28#include "AliAODEvent.h"
29#include "AliAODInputHandler.h"
30#include "AliAODMCParticle.h"
31
32#include "AliChaoticity.h"
33
34#define PI 3.1415927
35#define G_Coeff 0.006399 // 2*pi*alpha*M_pion
36
37
38// Author: Dhevan Gangadharan
39
40ClassImp(AliChaoticity)
41
42//________________________________________________________________________
43AliChaoticity::AliChaoticity():
44AliAnalysisTaskSE(),
45 fname(0),
46 fAOD(0x0),
cd12341d 47 fOutputList(0x0),
48 fPIDResponse(0x0),
49 fEC(0x0),
50 fEvt(0x0),
51 fTempStruct(0x0),
52 fRandomNumber(0x0),
53 fLEGO(kTRUE),
54 fMCcase(kFALSE),
55 fAODcase(kTRUE),
56 fPbPbcase(kTRUE),
57 fPdensityExplicitLoop(kFALSE),
58 fPdensityPairCut(kTRUE),
59 fTabulatePairs(kFALSE),
60 fBfield(0),
61 fMbin(0),
ae9b34d1 62 fFSIbin(0),
cd12341d 63 fEDbin(0),
64 fMbins(0),
65 fMultLimit(0),
66 fCentBinLowLimit(0),
67 fCentBinHighLimit(1),
68 fEventCounter(0),
69 fEventsToMix(0),
70 fZvertexBins(0),
71 fMultLimits(),
72 fQcut(),
73 fQLowerCut(0),
74 fNormQcutLow(),
75 fNormQcutHigh(),
76 fKupperBound(0),
77 fQupperBound(0),
78 fQupperBoundWeights(0),
79 fKstepT(),
80 fKstepY(),
81 fKmeanT(),
82 fKmeanY(),
83 fKmiddleT(),
84 fKmiddleY(),
85 fQstep(0),
5e3e77d6 86 fQstepWeights(0),
cd12341d 87 fQmean(),
88 fDampStart(0),
89 fDampStep(0),
cd12341d 90 fTPCTOFboundry(0),
91 fTOFboundry(0),
92 fSigmaCutTPC(0),
93 fSigmaCutTOF(0),
94 fMinSepTPCEntrancePhi(0),
95 fMinSepTPCEntranceEta(0),
96 fShareQuality(0),
97 fShareFraction(0),
98 fTrueMassP(0),
99 fTrueMassPi(0),
100 fTrueMassK(0),
101 fTrueMassKs(0),
102 fTrueMassLam(0),
103 fKtbinL(0),
104 fKtbinH(0),
105 fKybinL(0),
106 fKybinH(0),
107 fQobinL(0),
108 fQobinH(0),
109 fQsbinL(0),
110 fQsbinH(0),
111 fQlbinL(0),
112 fQlbinH(0),
113 fDummyB(0),
cd12341d 114 fDefaultsCharMult(),
115 fDefaultsCharSE(),
116 fDefaultsCharME(),
117 fDefaultsInt(),
118 fPairLocationSE(),
119 fPairLocationME(),
120 fTripletSkip1(),
121 fTripletSkip2(),
122 fOtherPairLocation1(),
123 fOtherPairLocation2(),
124 fNormPairSwitch(),
125 fPairSplitCut(),
5e3e77d6 126 fNormPairs(),
127 fFSI2SS(),
128 fFSI2OS(),
129 fFSIOmega0SS(),
bdf4d668 130 fFSIOmega0OS(),
131 fMomResC2(),
bdf4d668 132 fNormWeight(0x0),
133 fNormWeightErr(0x0)
cd12341d 134{
135 // Default constructor
136 for(Int_t mb=0; mb<fMbins; mb++){
137 for(Int_t edB=0; edB<kEDbins; edB++){
138 for(Int_t c1=0; c1<2; c1++){
139 for(Int_t c2=0; c2<2; c2++){
140 for(Int_t sc=0; sc<kSCLimit2; sc++){
141 for(Int_t term=0; term<2; term++){
142
143 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
144
145 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
146 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
147 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
148 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
149 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
150 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
151
152 }// term_2
153 }// SC_2
154
155 for(Int_t c3=0; c3<2; c3++){
156 for(Int_t sc=0; sc<kSCLimit3; sc++){
157 for(Int_t term=0; term<5; term++){
158
159 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
160 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
161 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
162 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
163 for(Int_t dt=0; dt<kDENtypes; dt++){
164 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
90814457 165 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm = 0x0;
166 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm = 0x0;
167 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal = 0x0;
168 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal = 0x0;
169 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared = 0x0;
170 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared = 0x0;
cd12341d 171 }//dt
172
173 }// term_3
174 }// SC_3
175 }//c3
176 }//c2
177 }//c1
178 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
179 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
180 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
181 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
182 }
183 }
184
185 }// ED
186 }// Mbin
187
188}
189//________________________________________________________________________
190AliChaoticity::AliChaoticity(const Char_t *name, Bool_t MCdecision, Bool_t Tabulatedecision, Bool_t PbPbdecision, Int_t lowCentBin, Int_t highCentBin, Bool_t lego)
191: AliAnalysisTaskSE(name),
192 fname(name),
193 fAOD(0x0),
b6e5ec54 194 //fESD(0x0),
cd12341d 195 fOutputList(0x0),
196 fPIDResponse(0x0),
197 fEC(0x0),
198 fEvt(0x0),
199 fTempStruct(0x0),
200 fRandomNumber(0x0),
201 fLEGO(lego),
202 fMCcase(MCdecision),
203 fAODcase(kTRUE),
204 fPbPbcase(PbPbdecision),
205 fPdensityExplicitLoop(kFALSE),
206 fPdensityPairCut(kTRUE),
207 fTabulatePairs(Tabulatedecision),
208 fBfield(0),
209 fMbin(0),
ae9b34d1 210 fFSIbin(0),
cd12341d 211 fEDbin(0),
212 fMbins(0),
213 fMultLimit(0),
214 fCentBinLowLimit(lowCentBin),
215 fCentBinHighLimit(highCentBin),
216 fEventCounter(0),
217 fEventsToMix(0),
218 fZvertexBins(0),
219 fMultLimits(),
220 fQcut(),
221 fQLowerCut(0),
222 fNormQcutLow(),
223 fNormQcutHigh(),
224 fKupperBound(0),
225 fQupperBound(0),
226 fQupperBoundWeights(0),
227 fKstepT(),
228 fKstepY(),
229 fKmeanT(),
230 fKmeanY(),
231 fKmiddleT(),
232 fKmiddleY(),
233 fQstep(0),
5e3e77d6 234 fQstepWeights(0),
cd12341d 235 fQmean(),
236 fDampStart(0),
237 fDampStep(0),
cd12341d 238 fTPCTOFboundry(0),
239 fTOFboundry(0),
240 fSigmaCutTPC(0),
241 fSigmaCutTOF(0),
242 fMinSepTPCEntrancePhi(0),
243 fMinSepTPCEntranceEta(0),
244 fShareQuality(0),
245 fShareFraction(0),
246 fTrueMassP(0),
247 fTrueMassPi(0),
248 fTrueMassK(0),
249 fTrueMassKs(0),
250 fTrueMassLam(0),
251 fKtbinL(0),
252 fKtbinH(0),
253 fKybinL(0),
254 fKybinH(0),
255 fQobinL(0),
256 fQobinH(0),
257 fQsbinL(0),
258 fQsbinH(0),
259 fQlbinL(0),
260 fQlbinH(0),
261 fDummyB(0),
cd12341d 262 fDefaultsCharMult(),
263 fDefaultsCharSE(),
264 fDefaultsCharME(),
265 fDefaultsInt(),
266 fPairLocationSE(),
267 fPairLocationME(),
268 fTripletSkip1(),
269 fTripletSkip2(),
270 fOtherPairLocation1(),
271 fOtherPairLocation2(),
272 fNormPairSwitch(),
273 fPairSplitCut(),
5e3e77d6 274 fNormPairs(),
275 fFSI2SS(),
276 fFSI2OS(),
277 fFSIOmega0SS(),
bdf4d668 278 fFSIOmega0OS(),
279 fMomResC2(),
bdf4d668 280 fNormWeight(0x0),
281 fNormWeightErr(0x0)
cd12341d 282{
283 // Main constructor
284 fLEGO = lego;
285 fAODcase=kTRUE;
286 fMCcase=MCdecision;
287 fTabulatePairs=Tabulatedecision;
288 fPbPbcase=PbPbdecision;
289 fPdensityExplicitLoop = kFALSE;
290 fPdensityPairCut = kTRUE;
291 fCentBinLowLimit = lowCentBin;
292 fCentBinHighLimit = highCentBin;
293
294 for(Int_t mb=0; mb<fMbins; mb++){
295 for(Int_t edB=0; edB<kEDbins; edB++){
296 for(Int_t c1=0; c1<2; c1++){
297 for(Int_t c2=0; c2<2; c2++){
298 for(Int_t sc=0; sc<kSCLimit2; sc++){
299 for(Int_t term=0; term<2; term++){
300
301 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
302
303 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
304 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
305 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = 0x0;
306 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = 0x0;
307 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = 0x0;
308 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = 0x0;
309
310 }// term_2
311 }// SC_2
312
313 for(Int_t c3=0; c3<2; c3++){
314 for(Int_t sc=0; sc<kSCLimit3; sc++){
315 for(Int_t term=0; term<5; term++){
316
317 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
318 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
319 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
320 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
321 for(Int_t dt=0; dt<kDENtypes; dt++){
322 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = 0x0;
90814457 323 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm = 0x0;
324 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm = 0x0;
325 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal = 0x0;
326 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal = 0x0;
327 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared = 0x0;
328 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared = 0x0;
cd12341d 329 }//dt
330
331 }// term_3
332 }// SC_3
333 }//c3
334 }//c2
335 }//c1
336 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
337 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
338 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = 0x0;
339 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = 0x0;
340 }
341 }
342
343 }// ED
344 }// Mbin
345
346
347 DefineOutput(1, TList::Class());
348}
349//________________________________________________________________________
350AliChaoticity::AliChaoticity(const AliChaoticity &obj)
351 : AliAnalysisTaskSE(obj.fname),
352 fname(obj.fname),
353 fAOD(obj.fAOD),
b6e5ec54 354 //fESD(obj.fESD),
cd12341d 355 fOutputList(obj.fOutputList),
356 fPIDResponse(obj.fPIDResponse),
357 fEC(obj.fEC),
358 fEvt(obj.fEvt),
359 fTempStruct(obj.fTempStruct),
360 fRandomNumber(obj.fRandomNumber),
361 fLEGO(obj.fLEGO),
362 fMCcase(obj.fMCcase),
363 fAODcase(obj.fAODcase),
364 fPbPbcase(obj.fPbPbcase),
365 fPdensityExplicitLoop(obj.fPdensityExplicitLoop),
366 fPdensityPairCut(obj.fPdensityPairCut),
367 fTabulatePairs(obj.fTabulatePairs),
368 fBfield(obj.fBfield),
369 fMbin(obj.fMbin),
ae9b34d1 370 fFSIbin(obj.fFSIbin),
cd12341d 371 fEDbin(obj.fEDbin),
372 fMbins(obj.fMbins),
373 fMultLimit(obj.fMultLimit),
374 fCentBinLowLimit(obj.fCentBinLowLimit),
375 fCentBinHighLimit(obj.fCentBinHighLimit),
376 fEventCounter(obj.fEventCounter),
377 fEventsToMix(obj.fEventsToMix),
378 fZvertexBins(obj.fZvertexBins),
379 fMultLimits(),
380 fQcut(),
381 fQLowerCut(obj.fQLowerCut),
382 fNormQcutLow(),
383 fNormQcutHigh(),
384 fKupperBound(obj.fKupperBound),
385 fQupperBound(obj.fQupperBound),
386 fQupperBoundWeights(obj.fQupperBoundWeights),
387 fKstepT(),
388 fKstepY(),
389 fKmeanT(),
390 fKmeanY(),
391 fKmiddleT(),
392 fKmiddleY(),
393 fQstep(obj.fQstep),
5e3e77d6 394 fQstepWeights(obj.fQstepWeights),
cd12341d 395 fQmean(),
396 fDampStart(obj.fDampStart),
397 fDampStep(obj.fDampStep),
cd12341d 398 fTPCTOFboundry(obj.fTPCTOFboundry),
399 fTOFboundry(obj.fTOFboundry),
400 fSigmaCutTPC(obj.fSigmaCutTPC),
401 fSigmaCutTOF(obj.fSigmaCutTOF),
402 fMinSepTPCEntrancePhi(obj.fMinSepTPCEntrancePhi),
403 fMinSepTPCEntranceEta(obj.fMinSepTPCEntranceEta),
404 fShareQuality(obj.fShareQuality),
405 fShareFraction(obj.fShareFraction),
406 fTrueMassP(obj.fTrueMassP),
407 fTrueMassPi(obj.fTrueMassPi),
408 fTrueMassK(obj.fTrueMassK),
409 fTrueMassKs(obj.fTrueMassKs),
410 fTrueMassLam(obj.fTrueMassLam),
411 fKtbinL(obj.fKtbinL),
412 fKtbinH(obj.fKtbinH),
413 fKybinL(obj.fKybinL),
414 fKybinH(obj.fKybinH),
415 fQobinL(obj.fQobinL),
416 fQobinH(obj.fQobinH),
417 fQsbinL(obj.fQsbinL),
418 fQsbinH(obj.fQsbinH),
419 fQlbinL(obj.fQlbinL),
420 fQlbinH(obj.fQlbinH),
421 fDummyB(obj.fDummyB),
cd12341d 422 fDefaultsCharMult(),
423 fDefaultsCharSE(),
424 fDefaultsCharME(),
425 fDefaultsInt(),
426 fPairLocationSE(),
427 fPairLocationME(),
428 fTripletSkip1(),
429 fTripletSkip2(),
430 fOtherPairLocation1(),
431 fOtherPairLocation2(),
432 fNormPairSwitch(),
433 fPairSplitCut(),
5e3e77d6 434 fNormPairs(),
435 fFSI2SS(),
436 fFSI2OS(),
437 fFSIOmega0SS(),
bdf4d668 438 fFSIOmega0OS(),
439 fMomResC2(),
bdf4d668 440 fNormWeight(),
441 fNormWeightErr()
cd12341d 442{
443 // Copy constructor
cd12341d 444}
445//________________________________________________________________________
446AliChaoticity &AliChaoticity::operator=(const AliChaoticity &obj)
447{
448 // Assignment operator
449 if (this == &obj)
450 return *this;
451
452 fname = obj.fname;
453 fAOD = obj.fAOD;
b6e5ec54 454 //fESD = obj.fESD;
cd12341d 455 fOutputList = obj.fOutputList;
456 fPIDResponse = obj.fPIDResponse;
457 fEC = obj.fEC;
458 fEvt = obj.fEvt;
459 fTempStruct = obj.fTempStruct;
460 fRandomNumber = obj.fRandomNumber;
461 fLEGO = fLEGO;
462 fMCcase = obj.fMCcase;
463 fAODcase = obj.fAODcase;
464 fPbPbcase = obj.fPbPbcase;
465 fPdensityExplicitLoop = obj.fPdensityExplicitLoop;
466 fPdensityPairCut = obj.fPdensityPairCut;
467 fTabulatePairs = obj.fTabulatePairs;
468 fBfield = obj.fBfield;
469 fMbin = obj.fMbin;
ae9b34d1 470 fFSIbin = obj.fFSIbin;
cd12341d 471 fEDbin = obj.fEDbin;
472 fMbins = obj.fMbins;
473 fMultLimit = obj.fMultLimit;
474 fCentBinLowLimit = obj.fCentBinLowLimit;
475 fCentBinHighLimit = obj.fCentBinHighLimit;
476 fEventCounter = obj.fEventCounter;
477 fEventsToMix = obj.fEventsToMix;
478 fZvertexBins = obj.fZvertexBins;
479 //fMultLimits = ;
480 //fQcut = ;
481 fQLowerCut = obj.fQLowerCut;
482 //fNormQcutLow = ;
483 //fNormQcutHigh = ;
484 fKupperBound = obj.fKupperBound;
485 fQupperBound = obj.fQupperBound;
486 fQupperBoundWeights = obj.fQupperBoundWeights;
487 //fKstepT = ;
488 //fKstepY = ;
489 //fKmeanT = ;
490 //fKmeanY = ;
491 //fKmiddleT = ;
492 //fKmiddleY = ;
493 fQstep = obj.fQstep;
5e3e77d6 494 fQstepWeights = obj.fQstepWeights;
cd12341d 495 //fQmean = ;
496 fDampStart = obj.fDampStart;
497 fDampStep = obj.fDampStep;
cd12341d 498 fTPCTOFboundry = obj.fTPCTOFboundry;
499 fTOFboundry = obj.fTOFboundry;
500 fSigmaCutTPC = obj.fSigmaCutTPC;
501 fSigmaCutTOF = obj.fSigmaCutTOF;
502 fMinSepTPCEntrancePhi = obj.fMinSepTPCEntrancePhi;
503 fMinSepTPCEntranceEta = obj.fMinSepTPCEntranceEta;
504 fShareQuality = obj.fShareQuality;
505 fShareFraction = obj.fShareFraction;
506 fTrueMassP = obj.fTrueMassP;
507 fTrueMassPi = obj.fTrueMassPi;
508 fTrueMassK = obj.fTrueMassK;
509 fTrueMassKs = obj.fTrueMassKs;
510 fTrueMassLam = obj.fTrueMassLam;
511 fKtbinL = obj.fKtbinL;
512 fKtbinH = obj.fKtbinH;
513 fKybinL = obj.fKybinL;
514 fKybinH = obj.fKybinH;
515 fQobinL = obj.fQobinL;
516 fQobinH = obj.fQobinH;
517 fQsbinL = obj.fQsbinL;
518 fQsbinH = obj.fQsbinH;
519 fQlbinL = obj.fQlbinL;
520 fQlbinH = obj.fQlbinH;
521 fDummyB = obj.fDummyB;
522 fNormWeight = obj.fNormWeight;
523 fNormWeightErr = obj.fNormWeightErr;
524
525 return (*this);
526}
527//________________________________________________________________________
528AliChaoticity::~AliChaoticity()
529{
530 // Destructor
531 if(fAOD) delete fAOD;
b6e5ec54 532 //if(fESD) delete fESD;
cd12341d 533 if(fOutputList) delete fOutputList;
534 if(fPIDResponse) delete fPIDResponse;
535 if(fEC) delete fEC;
536 if(fEvt) delete fEvt;
68b2d0b9 537 if(fTempStruct) delete [] fTempStruct;
cd12341d 538 if(fRandomNumber) delete fRandomNumber;
90814457 539 if(fFSI2SS[0]) delete fFSI2SS[0];
540 if(fFSI2OS[0]) delete fFSI2OS[0];
541 if(fFSI2SS[1]) delete fFSI2SS[1];
542 if(fFSI2OS[1]) delete fFSI2OS[1];
ae9b34d1 543 if(fFSIOmega0SS[0]) delete fFSIOmega0SS[0];
544 if(fFSIOmega0SS[1]) delete fFSIOmega0SS[1];
545 if(fFSIOmega0SS[2]) delete fFSIOmega0SS[2];
546 if(fFSIOmega0SS[3]) delete fFSIOmega0SS[3];
547 if(fFSIOmega0SS[4]) delete fFSIOmega0SS[4];
548 if(fFSIOmega0SS[5]) delete fFSIOmega0SS[5];
90814457 549 if(fFSIOmega0OS[0]) delete fFSIOmega0OS[0];
550 if(fFSIOmega0OS[1]) delete fFSIOmega0OS[1];
551 if(fFSIOmega0OS[2]) delete fFSIOmega0OS[2];
552 if(fFSIOmega0OS[3]) delete fFSIOmega0OS[3];
553 if(fFSIOmega0OS[4]) delete fFSIOmega0OS[4];
554 if(fFSIOmega0OS[5]) delete fFSIOmega0OS[5];
5e3e77d6 555 if(fMomResC2) delete fMomResC2;
bdf4d668 556 if(fNormWeight) delete fNormWeight;
557 if(fNormWeightErr) delete fNormWeightErr;
5e3e77d6 558
68b2d0b9 559 for(int i=0; i<fMultLimit; i++){
560 if(fPairLocationSE[i]) delete [] fPairLocationSE[i];
561 if(fPairLocationME[i]) delete [] fPairLocationME[i];
562 for(int j=0; j<2; j++){
563 if(fOtherPairLocation1[j][i]) delete [] fOtherPairLocation1[j][i];
564 if(fOtherPairLocation2[j][i]) delete [] fOtherPairLocation2[j][i];
565 }
566 for(int j=0; j<3; j++) if(fNormPairSwitch[j][i]) delete [] fNormPairSwitch[j][i];
567 for(int j=0; j<4; j++) if(fPairSplitCut[j][i]) delete [] fPairSplitCut[j][i];
568 }
569 for(int i=0; i<kPairLimit; i++) if(fTripletSkip1[i]) delete [] fTripletSkip1[i];
570 for(int i=0; i<2*kPairLimit; i++) if(fTripletSkip2[i]) delete [] fTripletSkip2[i];
571 for(int i=0; i<3; i++) if(fNormPairs[i]) delete [] fNormPairs[i];
cd12341d 572 //
573 for(Int_t mb=0; mb<fMbins; mb++){
574 for(Int_t edB=0; edB<kEDbins; edB++){
575 for(Int_t c1=0; c1<2; c1++){
576 for(Int_t c2=0; c2<2; c2++){
577 for(Int_t sc=0; sc<kSCLimit2; sc++){
578 for(Int_t term=0; term<2; term++){
579
580 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2;
581
582 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal;
583 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared;
584 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL;
585 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW;
586 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL;
587 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW;
fa109294 588 //
589 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv;
590 if(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW) delete Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW;
cd12341d 591 }// term_2
592 }// SC_2
593
594 for(Int_t c3=0; c3<2; c3++){
595 for(Int_t sc=0; sc<kSCLimit3; sc++){
596 for(Int_t term=0; term<5; term++){
597
598 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3;
599 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3;
600 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3;
601 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3;
90814457 602 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms;
603 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms;
604 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal;
605 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal;
606 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared;
607 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared;
fa109294 608 //
90814457 609 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3;
610 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3;
611 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3;
612 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3;
613 //
614 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2;
615 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2;
616 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2;
617 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2;
618
fa109294 619 //
cd12341d 620 for(Int_t dt=0; dt<kDENtypes; dt++){
621 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm;
fa109294 622 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm;
623 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm;
90814457 624 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal;
625 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal;
626 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared;
627 if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared;
628
cd12341d 629 }//dt
630
631 }// term_3
632 }// SC_3
633 }//c3
634 }//c2
635 }//c1
636 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
637 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
638 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD;
639 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD;
640 }
641 }
642
643 }// ED
644 }// Mbin
645
646
647}
648//________________________________________________________________________
649void AliChaoticity::ParInit()
650{
651 cout<<"AliChaoticity MyInit() call"<<endl;
652
653 fRandomNumber = new TRandom3();
654 fRandomNumber->SetSeed(0);
cd12341d 655
656 //
657 fEventCounter=0;
658 if(fPdensityExplicitLoop) fEventsToMix=3;
659 else if(fPdensityPairCut && !fPdensityExplicitLoop) fEventsToMix=2;
660 else fEventsToMix=0;
661 fZvertexBins=2;//2
662
663 fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
664 fTOFboundry = 2.1;// TOF pid used below this momentum
665 fSigmaCutTPC = 2.0;// 2.0 (norm). 1.5 Sys Deviation
666 fSigmaCutTOF = 2.0;// 2.0 (norm). 1.5 Sys Deviation
667
668 ////////////////////////////////////////////////
669 // Proton Pair Cuts
670 fMinSepTPCEntrancePhi = 0.035;// 0.035(norm). 0.04 Sys Deviation
671 fMinSepTPCEntranceEta = 0.035;// 0.035(norm). 0.04 Sys Deviation
672 fShareQuality = .5;// max
673 fShareFraction = .05;// max
674 ////////////////////////////////////////////////
5e3e77d6 675
cd12341d 676
677 fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
678 fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
679
680
90814457 681
cd12341d 682 if(fPbPbcase) {// PbPb
683 fMultLimit=kMultLimitPbPb;
684 fMbins=kCentBins;
685 fQcut[0]=0.1;
686 fQcut[1]=0.1;
687 fQcut[2]=0.6;
688 fNormQcutLow[0] = 0.15;//1.06 (test also at 0.15)
689 fNormQcutHigh[0] = 0.175;//1.1 (test also at 0.175)
690 fNormQcutLow[1] = 1.34;//1.34
691 fNormQcutHigh[1] = 1.4;//1.4
692 fNormQcutLow[2] = 1.1;//1.1
693 fNormQcutHigh[2] = 1.4;//1.4
694 }
695 else {// pp
696 fMultLimit=kMultLimitpp;
697 fMbins=kMultBinspp;
698 fQcut[0]=0.6;
699 fQcut[1]=0.6;
700 fQcut[2]=0.6;
701 fNormQcutLow[0] = 1.0;
702 fNormQcutHigh[0] = 1.5;
703 fNormQcutLow[1] = 1.0;
704 fNormQcutHigh[1] = 1.5;
705 fNormQcutLow[2] = 1.0;
706 fNormQcutHigh[2] = 1.5;
707 }
708
90814457 709 fQLowerCut = 0.005;// was 0.005
cd12341d 710 fKupperBound = 1.0;
711 //
712 // 4x1 (Kt: 0-0.2, 0.2-0.4, 0.4-0.6, 0.6-1.0)
713 //fKstepT[0] = 0.2; fKstepT[1] = 0.2; fKstepT[2] = 0.2; fKstepT[3] = 0.4;
714 //fKstepY[0] = 1.6;// -0.8 to +0.8
715 //fKmeanT[0] = 0.1; fKmeanT[1] = 0.3; fKmeanT[2] = 0.5; fKmeanT[3] = 0.8;
716 //fKmeanY[0] = 0;// central y
717 //fKmiddleT[0] = 0.1; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.5; fKmiddleT[3] = 0.8;
718 //fKmiddleY[0] = 0;
719 // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
720 fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
721 fKstepY[0] = 1.6;// -0.8 to +0.8
722 fKmeanT[0] = 0.241; fKmeanT[1] = 0.369; fKmeanT[2] = 0.573;
723 fKmeanY[0] = 0;// central y
724 fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
725 fKmiddleY[0] = 0;
726 // 2x1 (Kt: 0-0.4, 0.4-1.0)
727 //fKstepT[0] = 0.4; fKstepT[1] = 0.6;
728 //fKstepY[0] = 1.6;// -0.8 to +0.8
729 //fKmeanT[0] = 0.255; fKmeanT[1] = 0.505;
730 //fKmiddleT[0] = 0.2; fKmiddleT[1] = 0.7;
731 //fKmeanY[0] = 0;// central y
732 //fKmiddleY[0] = 0;
733 // 1x1 (Kt: 0-1.0)
734 //fKstepT[0] = 1.0;
735 //fKstepY[0] = 1.6;// -0.8 to +0.8
736 //fKmeanT[0] = 0.306;
737 //fKmiddleT[0] = 0.5;
738 //fKmeanY[0] = 0;// central y
739 //fKmiddleY[0] = 0;
740 //
741 fQupperBoundWeights = 0.2;
742 fQupperBound = 0.1;
5e3e77d6 743 fQstep = fQupperBound/Float_t(kQbins);
744 fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
745 for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
cd12341d 746 //
747 fDampStart = 0.3;
748 fDampStep = 0.02;
749
750
751
752 fEC = new AliChaoticityEventCollection **[fZvertexBins];
753 for(UShort_t i=0; i<fZvertexBins; i++){
754
755 fEC[i] = new AliChaoticityEventCollection *[fMbins];
756
757 for(UShort_t j=0; j<fMbins; j++){
758
759 fEC[i][j] = new AliChaoticityEventCollection(fEventsToMix+1, fMultLimit, kPairLimit, fMCcase);
760 }
761 }
762
763
764 for(Int_t i=0; i<fMultLimit; i++) fDefaultsCharMult[i]='0';
765 for(Int_t i=0; i<kPairLimit; i++) fDefaultsCharSE[i]='0';
766 for(Int_t i=0; i<2*kPairLimit; i++) fDefaultsCharME[i]='0';
767 for(Int_t i=0; i<fMultLimit; i++) fDefaultsInt[i]=-1;
768 for(Int_t i=0; i<fMultLimit; i++) fPairLocationSE[i] = new TArrayI(fMultLimit,fDefaultsInt);
769 for(Int_t i=0; i<fMultLimit; i++) fPairLocationME[i] = new TArrayI(fMultLimit,fDefaultsInt);
770 for(Int_t i=0; i<kPairLimit; i++) fTripletSkip1[i] = new TArrayC(fMultLimit,fDefaultsCharSE);
771 for(Int_t i=0; i<2*kPairLimit; i++) fTripletSkip2[i] = new TArrayC(fMultLimit,fDefaultsCharME);
772
773
774 // Normalization utilities
775 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
776 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
777 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
778 for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
779 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
780 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
781 for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
782
783 // Track Merging/Splitting utilities
784 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P11
785 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P12
786 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P13
787 for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[3][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P23
788
789
790 fNormPairs[0] = new AliChaoticityNormPairStruct[kNormPairLimit];
791 fNormPairs[1] = new AliChaoticityNormPairStruct[kNormPairLimit];
792
793
794 fTempStruct = new AliChaoticityTrackStruct[fMultLimit];
795
796
797 fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
798
799
800 // Initialize Weight Array
801 fNormWeight = new Float_t******[fMbins];// fMbin
802 fNormWeightErr = new Float_t******[fMbins];// fMbin
803 for(Int_t i=0; i<fMbins; i++){// Mbin iterator
804 fNormWeight[i] = new Float_t*****[kEDbins];// create ED bins
805 fNormWeightErr[i] = new Float_t*****[kEDbins];// create ED bins
806
807 for(Int_t j=0; j<kEDbins; j++){// ED iterator
808 fNormWeight[i][j] = new Float_t****[kKbinsT];// create Kt bins
809 fNormWeightErr[i][j] = new Float_t****[kKbinsT];// create Kt bins
810
811 for(Int_t k=0; k<kKbinsT; k++){// Kt iterator
812 fNormWeight[i][j][k] = new Float_t***[kKbinsY];// create Ky bins
813 fNormWeightErr[i][j][k] = new Float_t***[kKbinsY];// create Ky bins
814
815 for(Int_t l=0; l<kKbinsY; l++){// Ky iterator
816 fNormWeight[i][j][k][l] = new Float_t**[kQbinsWeights];// create Qout bins
817 fNormWeightErr[i][j][k][l] = new Float_t**[kQbinsWeights];// create Qout bins
818
819 for(Int_t m=0; m<kQbinsWeights; m++){// Qout iterator
820 fNormWeight[i][j][k][l][m] = new Float_t*[kQbinsWeights];// create Qside bins
821 fNormWeightErr[i][j][k][l][m] = new Float_t*[kQbinsWeights];// create Qside bins
822
823 for(Int_t n=0; n<kQbinsWeights; n++){// Qside iterator
824 fNormWeight[i][j][k][l][m][n] = new Float_t[kQbinsWeights];// create Qlong bins
825 fNormWeightErr[i][j][k][l][m][n] = new Float_t[kQbinsWeights];// create Qlong bins
826 }
827 }
828 }
829 }
830 }
831 }
832
833
834
835 // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
90814457 836 if(!fLEGO) {
5e3e77d6 837 SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
90814457 838 if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
839 if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
cd12341d 840 }
90814457 841
842 /////////////////////////////////////////////
843 // AddTaskChaoticity function call testing
844 //TestAddTask();
845 ////////////////////////////////////////////////
846
cd12341d 847}
848//________________________________________________________________________
849void AliChaoticity::UserCreateOutputObjects()
850{
851 // Create histograms
852 // Called once
853
854 ParInit();// Initialize my settings
5e3e77d6 855
cd12341d 856
857 fOutputList = new TList();
858 fOutputList->SetOwner();
859
860 TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
861 fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
862 fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
863 fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
864 fOutputList->Add(fVertexDist);
865
866
867 TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
868 fOutputList->Add(fDCAxyDistPlus);
869 TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
870 fOutputList->Add(fDCAzDistPlus);
871 TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
872 fOutputList->Add(fDCAxyDistMinus);
873 TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
874 fOutputList->Add(fDCAzDistMinus);
875
876
877 TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
878 fOutputList->Add(fEvents1);
879 TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
880 fOutputList->Add(fEvents2);
881
882 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
883 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
884 fOutputList->Add(fMultDist1);
885
886 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
887 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
888 fOutputList->Add(fMultDist2);
889
890 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
891 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
892 fOutputList->Add(fMultDist3);
893
894 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
895 fOutputList->Add(fPtEtaDist);
896
897 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
898 fOutputList->Add(fPhiPtDist);
899
900 TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
901 fOutputList->Add(fTOFResponse);
902 TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
903 fOutputList->Add(fTPCResponse);
904
905 TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",200,0,2);
906 fOutputList->Add(fRejectedPairs);
907 TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
908 fOutputList->Add(fRejectedEvents);
909
910 TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
911 if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
912 TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
913 if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
914
915 TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
916 if(fMCcase) fOutputList->Add(fResonanceOSPairs);
917 TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
918 if(fMCcase) fOutputList->Add(fAllOSPairs);
919
920 TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
921 fOutputList->Add(fAvgMult);
922
923 TH3D *fTPNRejects = new TH3D("fTPNRejects","",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
924 fOutputList->Add(fTPNRejects);
925
926 TH2D *fKt3Dist = new TH2D("fKt3Dist","",fMbins,.5,fMbins+.5, 10,0,1);
927 fOutputList->Add(fKt3Dist);
928
929
cd12341d 930 if(fPdensityExplicitLoop || fPdensityPairCut){
931
932 for(Int_t mb=0; mb<fMbins; mb++){
933 if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
934
935 for(Int_t edB=0; edB<kEDbins; edB++){
936 for(Int_t c1=0; c1<2; c1++){
937 for(Int_t c2=0; c2<2; c2++){
938 for(Int_t sc=0; sc<kSCLimit2; sc++){
939 for(Int_t term=0; term<2; term++){
940
941 TString *nameEx2 = new TString("Explicit2_Charge1_");
942 *nameEx2 += c1;
943 nameEx2->Append("_Charge2_");
944 *nameEx2 += c2;
945 nameEx2->Append("_SC_");
946 *nameEx2 += sc;
947 nameEx2->Append("_M_");
948 *nameEx2 += mb;
949 nameEx2->Append("_ED_");
950 *nameEx2 += edB;
951 nameEx2->Append("_Term_");
952 *nameEx2 += term+1;
953
954 if(sc==0 || sc==3 || sc==5){
955 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
956 }
957
958 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
959 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2);
5e3e77d6 960 TString *nameEx2QW=new TString(nameEx2->Data());
961 nameEx2QW->Append("_QW");
962 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
963 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW);
cd12341d 964 // Momentum resolution histos
965 if(fMCcase && sc==0){
966 TString *nameIdeal = new TString(nameEx2->Data());
967 nameIdeal->Append("_Ideal");
fa109294 968 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",kRVALUES*kNDampValues,-0.5,kRVALUES*kNDampValues-0.5, kQbinsWeights,0,fQupperBoundWeights);
cd12341d 969 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal);
970 TString *nameSmeared = new TString(nameEx2->Data());
971 nameSmeared->Append("_Smeared");
fa109294 972 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",kRVALUES*kNDampValues,-0.5,kRVALUES*kNDampValues-0.5, kQbinsWeights,0,fQupperBoundWeights);
cd12341d 973 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared);
fa109294 974 //
975 TString *nameEx2MC=new TString(nameEx2->Data());
976 nameEx2MC->Append("_MCqinv");
977 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"",400,0,2);
978 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
979 TString *nameEx2MCQW=new TString(nameEx2->Data());
980 nameEx2MCQW->Append("_MCqinvQW");
981 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"",400,0,2);
982 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
54d66278 983 //
984 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
985 nameEx2PIDpurityDen->Append("_PIDpurityDen");
986 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen = new TH2D(nameEx2PIDpurityDen->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
987 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
988 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
989 nameEx2PIDpurityNum->Append("_PIDpurityNum");
990 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH2D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",20,0,1, 400,0,2);
991 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
cd12341d 992 }
fa109294 993 if(sc==0){
cd12341d 994
995 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
996 nameEx2OSLB1->Append("_osl_b1");
fa109294 997 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
cd12341d 998 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL);
999 nameEx2OSLB1->Append("_QW");
fa109294 1000 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
cd12341d 1001 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW);
1002 //
1003 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
1004 nameEx2OSLB2->Append("_osl_b2");
fa109294 1005 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
cd12341d 1006 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL);
1007 nameEx2OSLB2->Append("_QW");
fa109294 1008 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
cd12341d 1009 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW);
1010
1011 }
1012
1013 }// term_2
1014 }// SC_2
1015
1016 // skip 3-particle if Tabulate6DPairs is true
1017 if(fTabulatePairs) continue;
1018
1019 for(Int_t c3=0; c3<2; c3++){
1020 for(Int_t sc=0; sc<kSCLimit3; sc++){
1021 for(Int_t term=0; term<5; term++){
1022 TString *nameEx3 = new TString("Explicit3_Charge1_");
1023 *nameEx3 += c1;
1024 nameEx3->Append("_Charge2_");
1025 *nameEx3 += c2;
1026 nameEx3->Append("_Charge3_");
1027 *nameEx3 += c3;
1028 nameEx3->Append("_SC_");
1029 *nameEx3 += sc;
1030 nameEx3->Append("_M_");
1031 *nameEx3 += mb;
1032 nameEx3->Append("_ED_");
1033 *nameEx3 += edB;
1034 nameEx3->Append("_Term_");
1035 *nameEx3 += term+1;
1036
1037 TString *namePC3 = new TString("PairCut3_Charge1_");
1038 *namePC3 += c1;
1039 namePC3->Append("_Charge2_");
1040 *namePC3 += c2;
1041 namePC3->Append("_Charge3_");
1042 *namePC3 += c3;
1043 namePC3->Append("_SC_");
1044 *namePC3 += sc;
1045 namePC3->Append("_M_");
1046 *namePC3 += mb;
1047 namePC3->Append("_ED_");
1048 *namePC3 += edB;
1049 namePC3->Append("_Term_");
1050 *namePC3 += term+1;
1051
1052 ///////////////////////////////////////
1053 // skip degenerate histograms
1054 if(sc==0 || sc==6 || sc==9){// Identical species
1055 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
1056 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
1057 }else if(sc!=5){
1058 if( (c1+c2)==1) {if(c1!=0) continue;}
1059 }else {}// do nothing for pi-k-p case
1060
1061 /////////////////////////////////////////
1062
1063 if(fPdensityExplicitLoop){
1064 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = new TH1D(nameEx3->Data(),"Three Particle Distribution",200,0,2);
1065 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3);
1066 //
1067 nameEx3->Append("_Norm");
1068 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = new TH1D(nameEx3->Data(),"Explicit_3 Norm",1,-0.5,0.5);
1069 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3);
1070 }
1071 if(fPdensityPairCut){
1072 TString *nameNorm=new TString(namePC3->Data());
1073 nameNorm->Append("_Norm");
1074 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1075 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1076 //
1077 if(sc<=2){
1078 TString *name3DQ=new TString(namePC3->Data());
1079 name3DQ->Append("_3D");
1080 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH3D(name3DQ->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1081 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1082 //
90814457 1083 const int NEdgesPos=16;
1084 double lowEdges4vectPos[NEdgesPos]={0};
1085 lowEdges4vectPos[0]=0.0;
1086 lowEdges4vectPos[1]=0.0001;// best resolution at low Q^2
1087 for(int edge=2; edge<NEdgesPos; edge++){
1088 lowEdges4vectPos[edge] = lowEdges4vectPos[edge-1] + lowEdges4vectPos[1]*(edge);
1089 }
1090 const int NEdges=2*NEdgesPos-1;
fa109294 1091 double lowEdges4vect[NEdges]={0};
90814457 1092 for(int edge=0; edge<NEdges; edge++){
1093 if(edge<NEdgesPos-1) lowEdges4vect[edge] = -lowEdges4vectPos[NEdgesPos-1-edge];
1094 else if(edge==NEdgesPos-1) lowEdges4vect[edge] = 0;
1095 else lowEdges4vect[edge] = lowEdges4vectPos[edge-NEdgesPos+1];
fa109294 1096 }
90814457 1097
5e3e77d6 1098 if(c1==c2 && c1==c3 && sc==0 && fMCcase==kFALSE){
90814457 1099 TString *name4vect1=new TString(namePC3->Data());
1100 TString *name4vect2=new TString(namePC3->Data());
1101 name4vect1->Append("_4VectProd1");
1102 name4vect2->Append("_4VectProd2");
5e3e77d6 1103 // use 3.75e6 MeV^4 as the resolution on QprodSum
90814457 1104 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms = new TH3D(name4vect1->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1105 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms);
1106 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms = new TH3D(name4vect2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1107 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms);
5e3e77d6 1108 }
1109 if(sc==0 && fMCcase==kTRUE){
1110 TString *name3DMomResIdeal=new TString(namePC3->Data());
1111 name3DMomResIdeal->Append("_Ideal");
1112 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH3D(name3DMomResIdeal->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1113 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1114 TString *name3DMomResSmeared=new TString(namePC3->Data());
1115 name3DMomResSmeared->Append("_Smeared");
1116 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH3D(name3DMomResSmeared->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1117 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared);
601fbb33 1118 //
1119 TString *name3DMomResQW12=new TString(namePC3->Data());
1120 name3DMomResQW12->Append("_QW12");
1121 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW12 = new TH3D(name3DMomResQW12->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1122 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW12);
1123 TString *name3DMomResQW13=new TString(namePC3->Data());
1124 name3DMomResQW13->Append("_QW13");
1125 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW13 = new TH3D(name3DMomResQW13->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1126 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW13);
90814457 1127 //
1128 if(term==0){
1129 TString *name3DSumK3=new TString(namePC3->Data());
1130 name3DSumK3->Append("_SumK3");
1131 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSumK3 = new TH3D(name3DSumK3->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1132 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSumK3);
1133 TString *name3DEnK3=new TString(namePC3->Data());
1134 name3DEnK3->Append("_EnK3");
1135 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fEnK3 = new TH3D(name3DEnK3->Data(),"", kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1136 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fEnK3);
1137 }
1138
1139 if(c1==c2 && c1==c3){
1140 TString *name4vect1Ideal=new TString(namePC3->Data());
1141 TString *name4vect1Smeared=new TString(namePC3->Data());
1142 TString *name4vect2Ideal=new TString(namePC3->Data());
1143 TString *name4vect2Smeared=new TString(namePC3->Data());
1144 name4vect1Ideal->Append("_4VectProd1Ideal");
1145 name4vect1Smeared->Append("_4VectProd1Smeared");
1146 name4vect2Ideal->Append("_4VectProd2Ideal");
1147 name4vect2Smeared->Append("_4VectProd2Smeared");
1148 // use 3.75e6 MeV^4 as the resolution on QprodSum
1149 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal = new TH3D(name4vect1Ideal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1150 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal);
1151 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared = new TH3D(name4vect1Smeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1152 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared);
1153 //
1154 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal = new TH3D(name4vect2Ideal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1155 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal);
1156 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared = new TH3D(name4vect2Smeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1157 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared);
1158 //
1159 if(term==0){
1160 TString *name4vect1SumK3=new TString(namePC3->Data());
1161 TString *name4vect2SumK3=new TString(namePC3->Data());
1162 TString *name4vect1EnK3=new TString(namePC3->Data());
1163 TString *name4vect2EnK3=new TString(namePC3->Data());
1164 name4vect1SumK3->Append("_4VectProd1SumK3");
1165 name4vect2SumK3->Append("_4VectProd2SumK3");
1166 name4vect1EnK3->Append("_4VectProd1EnK3");
1167 name4vect2EnK3->Append("_4VectProd2EnK3");
1168 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3 = new TH3D(name4vect1SumK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1169 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3);
1170 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3 = new TH3D(name4vect2SumK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1171 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3);
1172 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3 = new TH3D(name4vect1EnK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1173 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3);
1174 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3 = new TH3D(name4vect2EnK3->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1175 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3);
1176 }// term 0
1177 if(term > 0 && term < 4){
1178 TString *name4vect1SumK2=new TString(namePC3->Data());
1179 TString *name4vect2SumK2=new TString(namePC3->Data());
1180 TString *name4vect1EnK2=new TString(namePC3->Data());
1181 TString *name4vect2EnK2=new TString(namePC3->Data());
1182 name4vect1SumK2->Append("_4VectProd1SumK2");
1183 name4vect2SumK2->Append("_4VectProd2SumK2");
1184 name4vect1EnK2->Append("_4VectProd1EnK2");
1185 name4vect2EnK2->Append("_4VectProd2EnK2");
1186 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2 = new TH3D(name4vect1SumK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1187 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2);
1188 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2 = new TH3D(name4vect2SumK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1189 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2);
1190 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2 = new TH3D(name4vect1EnK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1191 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2);
1192 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2 = new TH3D(name4vect2EnK2->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1193 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2);
1194 }// terms 1,2,3
1195 }
1196 }// MCcase
5e3e77d6 1197 //
90814457 1198 if(c1==c2 && c1==c3 && term==4 && sc==0){
cd12341d 1199 for(Int_t dt=0; dt<kDENtypes; dt++){
1200 TString *nameDenType=new TString("PairCut3_Charge1_");
1201 *nameDenType += c1;
1202 nameDenType->Append("_Charge2_");
1203 *nameDenType += c2;
1204 nameDenType->Append("_Charge3_");
1205 *nameDenType += c3;
1206 nameDenType->Append("_SC_");
1207 *nameDenType += sc;
1208 nameDenType->Append("_M_");
1209 *nameDenType += mb;
1210 nameDenType->Append("_ED_");
1211 *nameDenType += edB;
1212 nameDenType->Append("_TPN_");
1213 *nameDenType += dt;
1214
1215 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm = new TH3D(nameDenType->Data(),"",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1216 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm);
1217 // neglect errors for TPN
1218 //nameDenType->Append("_Err");
1219 //Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr = new TH3D(nameDenType->Data(),"",kQbins,0,fQupperBound, kQbins,0,fQupperBound, kQbins,0,fQupperBound);
1220 //fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr);
5e3e77d6 1221 //
fa109294 1222 TString *name4vect1TPN=new TString(nameDenType->Data());
1223 TString *name4vect2TPN=new TString(nameDenType->Data());
1224 name4vect1TPN->Append("_4VectProd1");
1225 name4vect2TPN->Append("_4VectProd2");
1226 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm = new TH3D(name4vect1TPN->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1227 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm);
1228 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm = new TH3D(name4vect2TPN->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1229 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm);
90814457 1230 //
1231 if(fMCcase){
1232 TString *name4vect1TPNIdeal=new TString(nameDenType->Data());
1233 TString *name4vect2TPNIdeal=new TString(nameDenType->Data());
1234 TString *name4vect1TPNSmeared=new TString(nameDenType->Data());
1235 TString *name4vect2TPNSmeared=new TString(nameDenType->Data());
1236 name4vect1TPNIdeal->Append("_4VectProd1Ideal");
1237 name4vect2TPNIdeal->Append("_4VectProd2Ideal");
1238 name4vect1TPNSmeared->Append("_4VectProd1Smeared");
1239 name4vect2TPNSmeared->Append("_4VectProd2Smeared");
1240 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal = new TH3D(name4vect1TPNIdeal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1241 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal);
1242 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal = new TH3D(name4vect2TPNIdeal->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1243 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal);
1244 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared = new TH3D(name4vect1TPNSmeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1245 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared);
1246 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared = new TH3D(name4vect2TPNSmeared->Data(),"",NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect, NEdges-1,lowEdges4vect);
1247 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared);
1248 }
1249
cd12341d 1250 }
1251
1252 }// term=4
1253 }// c and sc exclusion
1254 }// PdensityPairCut
1255 }// term_3
1256 }// SC_3
1257 }//c3
1258 }//c2
1259 }//c1
1260 }// ED
1261 }// mbin
1262 }// Pdensity Method
1263
1264
1265 if(fTabulatePairs){
1266
1267 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
1268 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
1269 for(Int_t mb=0; mb<fMbins; mb++){
1270 for(Int_t edB=0; edB<kEDbins; edB++){
1271
1272 TString *nameNum = new TString("TwoPart_num_Kt_");
1273 *nameNum += tKbin;
1274 nameNum->Append("_Ky_");
1275 *nameNum += yKbin;
1276 nameNum->Append("_M_");
1277 *nameNum += mb;
1278 nameNum->Append("_ED_");
1279 *nameNum += edB;
1280
1281 TString *nameDen = new TString("TwoPart_den_Kt_");
1282 *nameDen += tKbin;
1283 nameDen->Append("_Ky_");
1284 *nameDen += yKbin;
1285 nameDen->Append("_M_");
1286 *nameDen += mb;
1287 nameDen->Append("_ED_");
1288 *nameDen += edB;
1289
1290
1291 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = new TH3I(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1292 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD);
1293
1294 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = new TH3I(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1295 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD);
1296 }
1297 }
1298 }
1299 }
1300
1301 }
1302
1303
1304 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1305 fOutputList->Add(fQsmearMean);
1306 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1307 fOutputList->Add(fQsmearSq);
1308 TH1D *fQDist = new TH1D("fQDist","",200,-.2,.2);
1309 fOutputList->Add(fQDist);
1310
1311
1312 ////////////////////////////////////
1313 ///////////////////////////////////
1314
1315 PostData(1, fOutputList);
1316}
1317
1318//________________________________________________________________________
1319void AliChaoticity::Exec(Option_t *)
1320{
1321 // Main loop
1322 // Called for each event
1ccd6f0d 1323 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
cd12341d 1324 fEventCounter++;
1325
b6e5ec54 1326 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
cd12341d 1327
b6e5ec54 1328 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1329 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
cd12341d 1330
cd12341d 1331
1332 // Trigger Cut
1333 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1334 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
fa109294 1335 if(!isSelected1 && !fMCcase) {return;}
5e3e77d6 1336 }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
cd12341d 1337 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1338 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
fa109294 1339 if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1340 }else {return;}
cd12341d 1341
1342 ///////////////////////////////////////////////////////////
1343 const AliAODVertex *primaryVertexAOD;
1344 AliCentrality *centrality;// for AODs and ESDs
1345
1346
1347 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1348 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1349 fPIDResponse = inputHandler->GetPIDResponse();
1350
1351
1352 TClonesArray *mcArray = 0x0;
1353 if(fMCcase){
1354 if(fAODcase){
1355 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1356 if(!mcArray || mcArray->GetEntriesFast() >= 110000){
1357 cout<<"No MC particle branch found or Array too large!!"<<endl;
cd12341d 1358 return;
1359 }
1360 }
1361 }
1362
1363
1364 UInt_t status=0;
1365 Int_t positiveTracks=0, negativeTracks=0;
1366 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1367
1368 Double_t vertex[3]={0};
1369 Int_t zbin=0;
1370 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1371 /////////////////////////////////////////////////
1372
1373
1374 Float_t centralityPercentile=0;
1375 Float_t cStep=5.0, cStart=0;
1376
1377
1378 if(fAODcase){// AOD case
1379
1380 if(fPbPbcase){
1381 centrality = fAOD->GetCentrality();
1382 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1383 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
5e3e77d6 1384 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range. Skipping Event"<<endl;*/ return;}
cd12341d 1385 cout<<"Centrality % = "<<centralityPercentile<<endl;
1386 }
1387
1388
1389
1390
1391 ////////////////////////////////
1392 // Vertexing
1393 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1394 primaryVertexAOD = fAOD->GetPrimaryVertex();
1395 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1396
1397 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1398 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1399
1400 if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Reject Pile-up events
1401 if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1402
1403 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1404
1405 fBfield = fAOD->GetMagneticField();
1406
1407 for(Int_t i=0; i<fZvertexBins; i++){
1408 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1409 zbin=i;
1410 break;
1411 }
1412 }
1413
1414
1415
1416 /////////////////////////////
1417 // Create Shuffled index list
1418 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1419 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1420 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1421 /////////////////////////////
1422
1423 // Track loop
1424 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1425 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1426 if (!aodtrack) continue;
1427 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1428
1429 status=aodtrack->GetStatus();
1430
1431 if((Bool_t)(((1<<7) & aodtrack->GetFilterMap()) == 0)) continue;// AOD filterBit cut
1432 // 1<<5 is for "standard cuts with tight dca cut"
1433 // 1<<7 is for TPC only tracking
1434 if(aodtrack->Pt() < 0.16) continue;
1435 if(fabs(aodtrack->Eta()) > 0.8) continue;
1436
1437
1438 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1439 if(!goodMomentum) continue;
1440 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1441
1442 Float_t dca2[2];
1443 Float_t dca3d;
1444
1445 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1446 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1447 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1448
1449 fTempStruct[myTracks].fStatus = status;
1450 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1451 fTempStruct[myTracks].fId = aodtrack->GetID();
1452 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1453 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1454 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1455 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1456 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1457 fTempStruct[myTracks].fEta = aodtrack->Eta();
1458 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1459 fTempStruct[myTracks].fDCAXY = dca2[0];
1460 fTempStruct[myTracks].fDCAZ = dca2[1];
1461 fTempStruct[myTracks].fDCA = dca3d;
1462 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1463 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1464
1465
1466
1467 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1468 if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1469 if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1470
1471 if(fTempStruct[myTracks].fCharge==+1) {
1472 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1473 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1474 }else {
1475 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1476 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1477 }
1478
1479 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1480 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1481
1482
1483 // nSimga PID workaround
1484 fTempStruct[myTracks].fElectron = kFALSE;
1485 fTempStruct[myTracks].fPion = kFALSE;
1486 fTempStruct[myTracks].fKaon = kFALSE;
1487 fTempStruct[myTracks].fProton = kFALSE;
1488
1489 Float_t nSigmaTPC[5];
1490 Float_t nSigmaTOF[5];
1491 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1492 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1493 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1494 Float_t signalTPC=0, signalTOF=0;
1495 Double_t integratedTimesTOF[10]={0};
1496 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1497 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1498 if (!aodTrack2) continue;
1499 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1500
1501 UInt_t status2=aodTrack2->GetStatus();
1502
1503 nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
1504 nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
1505 nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
1506 nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
1507 nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
1508 //
1509 nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
1510 nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
1511 nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
1512 nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
1513 nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
1514 signalTPC = aodTrack2->GetTPCsignal();
1515
1516 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1517 fTempStruct[myTracks].fTOFhit = kTRUE;
1518 signalTOF = aodTrack2->GetTOFsignal();
1519 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1520 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1521
1522 }// aodTrack2
5e3e77d6 1523
c4980714 1524 //cout<<nSigmaTPC[2]<<endl;
cd12341d 1525 ///////////////////
1526 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1527 if(fTempStruct[myTracks].fTOFhit) {
1528 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1529 }
1530 ///////////////////
1531
1532 // Use TOF if good hit and above threshold
1533 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1534 if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1535 if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1536 if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1537 if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1538 }else {// TPC info instead
1539 if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1540 if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1541 if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1542 if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1543 }
1544
fa109294 1545 //cout<<nSigmaTPC[2]<<endl;
cd12341d 1546 //////////////////////////////////////
1547 // Bayesian PIDs for TPC only tracking
1548 //const Double_t* PID = aodtrack->PID();
1549 //fTempStruct[myTracks].fElectron = kFALSE;
1550 //fTempStruct[myTracks].Pion = kFALSE;
1551 //fTempStruct[myTracks].Kaon = kFALSE;
1552 //fTempStruct[myTracks].Proton = kFALSE;
1553
1554 // Pions
1555 //if(PID[2] > 0.2) fTempStruct[myTracks].Pion = kTRUE;// pions: 0.2 --> 0.5
1556 //
1557 // Kaons
1558 //if(PID[3] <= 0.5) fTempStruct[myTracks].Kaon = kFALSE;// kaons
1559 //
1560 // Protons
1561 //if(PID[4] <= 0.5) fTempStruct[myTracks].Proton = kFALSE;// protons
1562 //////////////////////////////////////
1563
1564
1565 // Ensure there is only 1 candidate per track
1566 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1567 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1568 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1569 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1570 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1571 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1572 ////////////////////////
1573 if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
1574
1575 if(!fTempStruct[myTracks].fPion) continue;// only pions
1576
1577
1578
1579
1580 if(fTempStruct[myTracks].fPion) {// pions
1581 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1582 fTempStruct[myTracks].fKey = 1;
1583 }else if(fTempStruct[myTracks].fKaon){// kaons
1584 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1585 fTempStruct[myTracks].fKey = 10;
1586 }else{// protons
1587 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1588 fTempStruct[myTracks].fKey = 100;
1589 }
1590
1591
1592
1593 if(aodtrack->Charge() > 0) positiveTracks++;
1594 else negativeTracks++;
1595
1596 if(fTempStruct[myTracks].fPion) pionCount++;
1597 if(fTempStruct[myTracks].fKaon) kaonCount++;
1598 if(fTempStruct[myTracks].fProton) protonCount++;
1599
1600 myTracks++;
1601
1602 }
1603 }else {// ESD tracks
1604 cout<<"ESDs not supported currently"<<endl;
1605 return;
1606 }
1607
1608
1609 if(myTracks >= 1) {
1610 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1611 }
1612
1613
b6e5ec54 1614 //cout<<"There are "<<myTracks<<" myTracks"<<endl;
1615 //cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
cd12341d 1616
1617 /////////////////////////////////////////
1618 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1619 if(myTracks < 3) {cout<<"Less than 3 tracks. Skipping Event."<<endl; return;}
1620 /////////////////////////////////////////
fa109294 1621
cd12341d 1622
1623 ////////////////////////////////
1624 ///////////////////////////////
1625 // Mbin determination
1626 //
1627 // Mbin set to Pion Count Only for pp!!!!!!!
1628 fMbin=-1;
1629 if(!fPbPbcase){
1630 for(Int_t i=0; i<kMultBinspp; i++){
1631 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1632 // Mbin 0 has 1 pion
1633 }
1634 }else{
1635 for(Int_t i=0; i<kCentBins; i++){
1636 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1637 fMbin=i;// 0 = most central
1638 break;
1639 }
1640 }
1641 }
1642
1643 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1644
ae9b34d1 1645 if(fMbin==0) fFSIbin = 0;//0-5%
1646 else if(fMbin==1) fFSIbin = 1;//5-10%
1647 else if(fMbin<=3) fFSIbin = 2;//10-20%
1648 else if(fMbin<=5) fFSIbin = 3;//20-30%
1649 else if(fMbin<=7) fFSIbin = 4;//30-40%
1650 else fFSIbin = 5;//40-50%
1651
90814457 1652 Int_t rIndexForTPN = 4;
1653 if(fMbin<=1) {rIndexForTPN=5;}
1654 else if(fMbin<=3) {rIndexForTPN=4;}
1655 else if(fMbin<=5) {rIndexForTPN=3;}
1656 else {rIndexForTPN=2;}
1657
cd12341d 1658 //////////////////////////////////////////////////
1659 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1660 //////////////////////////////////////////////////
1661
1662
cd12341d 1663 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1664 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1665
1666 ////////////////////////////////////
1667 // Add event to buffer if > 0 tracks
1668 if(myTracks > 0){
1669 fEC[zbin][fMbin]->FIFOShift();
1670 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1671 (fEvt)->fNtracks = myTracks;
1672 (fEvt)->fFillStatus = 1;
1673 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1674 if(fMCcase){
1675 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1676 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1677 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1678 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1679 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1680 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1681 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1682 }
1683 }
1684 }
1685
1686
1687
1688 Float_t qinv12=0, qinv13=0, qinv23=0;
1689 Float_t qout=0, qside=0, qlong=0;
1690 Float_t qoutMC=0, qsideMC=0, qlongMC=0;
1691 Float_t transK12=0, rapK12=0, transK3=0;
1692 Int_t transKbin=0, rapKbin=0;
1693 Float_t q3=0;
1694 Int_t ch1=0, ch2=0, ch3=0;
1695 Short_t key1=0, key2=0, key3=0;
1696 Int_t bin1=0, bin2=0, bin3=0;
1697 Float_t pVect1[4]={0};
1698 Float_t pVect2[4]={0};
1699 Float_t pVect3[4]={0};
1700 Float_t pVect1MC[4]={0};
1701 Float_t pVect2MC[4]={0};
5e3e77d6 1702 Float_t pVect3MC[4]={0};
cd12341d 1703 Int_t index1=0, index2=0, index3=0;
1704 Float_t weight12=0, weight13=0, weight23=0;
1705 Float_t weight12Err=0, weight13Err=0, weight23Err=0;
1706 Float_t weight12CC=0, weight13CC=0, weight23CC=0;
1707 Float_t weightTotal=0;//, weightTotalErr=0;
5e3e77d6 1708 Float_t qinv12MC=0, qinv13MC=0, qinv23MC=0;
90814457 1709 Float_t Qsum1v1=0, Qsum2=0, Qsum3v1=0, Qsum1v2=0, Qsum3v2=0;
1710 Float_t Qsum1v1MC=0, Qsum2MC=0, Qsum3v1MC=0, Qsum1v2MC=0, Qsum3v2MC=0;
1711 //
cd12341d 1712 AliAODMCParticle *mcParticle1=0x0;
1713 AliAODMCParticle *mcParticle2=0x0;
1714
1715
1716 if(fPdensityPairCut){
1717 ////////////////////
1718 Int_t pairCountSE=0, pairCountME=0;
1719 Int_t normPairCount[2]={0};
1720 Int_t numOtherPairs1[2][fMultLimit];
1721 Int_t numOtherPairs2[2][fMultLimit];
1722 Bool_t exitCode=kFALSE;
1723 Int_t tempNormFillCount[2][2][2][10][5];
1724
1725
1726 // reset to defaults
1727 for(Int_t i=0; i<fMultLimit; i++) {
1728 fPairLocationSE[i]->Set(fMultLimit,fDefaultsInt);
1729 fPairLocationME[i]->Set(fMultLimit,fDefaultsInt);
1730
1731 // Normalization Utilities
1732 fOtherPairLocation1[0][i]->Set(fMultLimit,fDefaultsInt);
1733 fOtherPairLocation1[1][i]->Set(fMultLimit,fDefaultsInt);
1734 fOtherPairLocation2[0][i]->Set(fMultLimit,fDefaultsInt);
1735 fOtherPairLocation2[1][i]->Set(fMultLimit,fDefaultsInt);
1736 fNormPairSwitch[0][i]->Set(fMultLimit,fDefaultsCharMult);
1737 fNormPairSwitch[1][i]->Set(fMultLimit,fDefaultsCharMult);
1738 fNormPairSwitch[2][i]->Set(fMultLimit,fDefaultsCharMult);
1739 numOtherPairs1[0][i]=0;
1740 numOtherPairs1[1][i]=0;
1741 numOtherPairs2[0][i]=0;
1742 numOtherPairs2[1][i]=0;
1743
1744 // Track Merging/Splitting Utilities
1745 fPairSplitCut[0][i]->Set(fMultLimit,fDefaultsCharMult);// P11
1746 fPairSplitCut[1][i]->Set(fMultLimit,fDefaultsCharMult);// P12
1747 fPairSplitCut[2][i]->Set(fMultLimit,fDefaultsCharMult);// P13
1748 fPairSplitCut[3][i]->Set(fMultLimit,fDefaultsCharMult);// P23
1749 }
1750
1751 // Reset the temp Normalization counters
1752 for(Int_t i=0; i<2; i++){// Charge1
1753 for(Int_t j=0; j<2; j++){// Charge2
1754 for(Int_t k=0; k<2; k++){// Charge3
1755 for(Int_t l=0; l<10; l++){// FillIndex (species Combination)
1756 for(Int_t m=0; m<5; m++){// Term (Cumulant term)
1757 tempNormFillCount[i][j][k][l][m] = 0;
1758 }
1759 }
1760 }
1761 }
1762 }
1763
1764
1765 ///////////////////////////////////////////////////////
1766 // Start the pairing process
1767 // P11 pairing
1768 // 1st Particle
1769
1770 for (Int_t i=0; i<myTracks; i++) {
1771
1772 Int_t en2=0;
1773
1774 // 2nd particle
1775 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {
1776
1777 key1 = (fEvt)->fTracks[i].fKey;
1778 key2 = (fEvt+en2)->fTracks[j].fKey;
1779 Short_t fillIndex2 = FillIndex2part(key1+key2);
1780 Short_t qCutBin = SetQcutBin(fillIndex2);
1781 Short_t normBin = SetNormBin(fillIndex2);
1782 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1783 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1784 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1785 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1786 //
1787
1788 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1789 GetQosl(pVect1, pVect2, qout, qside, qlong);
1790 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1791
1792 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1793
1794 ///////////////////////////////
1795 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1796 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1797 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1798
1799 if(fMCcase && ch1==ch2 && fMbin==0){
1800 for(Int_t rstep=0; rstep<10; rstep++){
1801 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1802 // propagate through B field to r=1.2m
1803 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1804 if(phi1 > 2*PI) phi1 -= 2*PI;
1805 if(phi1 < 0) phi1 += 2*PI;
1806 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1807 if(phi2 > 2*PI) phi2 -= 2*PI;
1808 if(phi2 < 0) phi2 += 2*PI;
1809 Float_t deltaphi = phi1 - phi2;
1810 if(deltaphi > PI) deltaphi -= PI;
1811 if(deltaphi < -PI) deltaphi += PI;
1812 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1813 }
1814 }
5e3e77d6 1815
cd12341d 1816 // Pair Splitting/Merging cut
1817 if(ch1 == ch2){
1818 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1819 fPairSplitCut[0][i]->AddAt('1',j);
1820 ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1821 continue;
1822 }
1823 }
5e3e77d6 1824
cd12341d 1825 // HIJING tests
1826 if(fMCcase && fillIndex2==0){
1827
1828 // Check that label does not exceed stack size
1829 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
1830 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1831 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1832 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1833 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1834 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1835 qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
1836 GetQosl(pVect1MC, pVect2MC, qoutMC, qsideMC, qlongMC);
1837 if(qinv12<0.1 && ch1==ch2) {
1838 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
1839 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
1840 ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
1841 }
1842
1843
5e3e77d6 1844 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 8fm + 1 Therminator setting
cd12341d 1845 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
1846 Int_t denIndex = rIter*kNDampValues + myDampIt;
cd12341d 1847 Float_t WInput = MCWeight(ch1,ch2, rIter, myDampIt, qinv12MC);
1848 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
1849 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
1850 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
1851 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
1852 }
1853 }
1854
fa109294 1855 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
1856 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
1857
cd12341d 1858 //HIJING resonance test
1859 if(ch1 != ch2){
cd12341d 1860 ((TH1F*)fOutputList->FindObject("fAllOSPairs"))->Fill(fMbin+1, qinv12);
fa109294 1861 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
cd12341d 1862 if(mcParticle1->GetMother() == mcParticle2->GetMother() && mcParticle1->GetMother() >=0){
1863 ((TH1F*)fOutputList->FindObject("fResonanceOSPairs"))->Fill(fMbin+1, qinv12);
1864 }
1865 }
1866 }
90814457 1867 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinv->Fill(qinv12MC, MCWeight(ch1,ch2,4,5,qinv12MC));
1868 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(ch1,ch2,4,5,qinv12MC));
54d66278 1869 // pion purity
1870 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityDen->Fill(transK12, qinv12);
1871 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
1872 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityNum->Fill(transK12, qinv12);
1873 }
1874
cd12341d 1875 }// label check 2
1876 }// MC case
1877
1878 //////////////////////////////////////////
1879 // 2-particle term
1880 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
5e3e77d6 1881 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
fa109294 1882
1883
cd12341d 1884 // osl frame
fa109294 1885 if(fillIndex2==0){
cd12341d 1886 if((transK12 > 0.2) && (transK12 < 0.3)){
1887 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1888 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1889 }
1890 if((transK12 > 0.6) && (transK12 < 0.7)){
1891 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1892 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1893 }
1894 }
5e3e77d6 1895
cd12341d 1896 //////////////////////////////////////////
1897 if(fTabulatePairs){
1898 if(fillIndex2==0 && bin1==bin2){
1899 rapK12 = 0;
1900 transKbin=-1; rapKbin=-1;
1901 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1902 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1903 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1904 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1905 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1906 continue;
1907 }
1908 }
1909
1910
1911 // exit out of loop if there are too many pairs
1912 if(pairCountSE >= kPairLimit) {cout<<"Too many SE pairs"<<endl; exitCode=kTRUE; continue;}
1913 if(exitCode) continue;
5e3e77d6 1914
cd12341d 1915 //////////////////////////
1916 // Enforce the Qcut
1917 if(qinv12 <= fQcut[qCutBin]) {
5e3e77d6 1918
cd12341d 1919 ///////////////////////////
1920 // particle 1
1921 (fEvt)->fPairsSE[pairCountSE].fP1[0] = (fEvt)->fTracks[i].fP[0];
1922 (fEvt)->fPairsSE[pairCountSE].fP1[1] = (fEvt)->fTracks[i].fP[1];
1923 (fEvt)->fPairsSE[pairCountSE].fP1[2] = (fEvt)->fTracks[i].fP[2];
1924 (fEvt)->fPairsSE[pairCountSE].fE1 = (fEvt)->fTracks[i].fEaccepted;
1925 (fEvt)->fPairsSE[pairCountSE].fCharge1 = (fEvt)->fTracks[i].fCharge;
1926 (fEvt)->fPairsSE[pairCountSE].fIndex1 = i;
1927 (fEvt)->fPairsSE[pairCountSE].fKey1 = key1;
1928 (fEvt)->fPairsSE[pairCountSE].fLabel1 = (fEvt)->fTracks[i].fLabel;
1929 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1930 (fEvt)->fPairsSE[pairCountSE].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1931 (fEvt)->fPairsSE[pairCountSE].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1932 (fEvt)->fPairsSE[pairCountSE].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
5e3e77d6 1933 }
cd12341d 1934 // particle 2
1935 (fEvt)->fPairsSE[pairCountSE].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1936 (fEvt)->fPairsSE[pairCountSE].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1937 (fEvt)->fPairsSE[pairCountSE].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1938 (fEvt)->fPairsSE[pairCountSE].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1939 (fEvt)->fPairsSE[pairCountSE].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1940 (fEvt)->fPairsSE[pairCountSE].fIndex2 = j;
1941 (fEvt)->fPairsSE[pairCountSE].fKey2 = key2;
1942 (fEvt)->fPairsSE[pairCountSE].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1943 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1944 (fEvt)->fPairsSE[pairCountSE].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1945 (fEvt)->fPairsSE[pairCountSE].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1946 (fEvt)->fPairsSE[pairCountSE].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1947 }
5e3e77d6 1948
cd12341d 1949 (fEvt)->fPairsSE[pairCountSE].fQinv = qinv12;
1950
1951 fPairLocationSE[i]->AddAt(pairCountSE,j);
1952
1953 pairCountSE++;
1954
1955 }
1956
1957
1958 /////////////////////////////////////////////////////////
1959 // Normalization Region
1960
1961 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1962 // particle 1
1963 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1964 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1965 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1966 // particle 2
1967 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1968 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1969 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1970
1971
1972 //other past pairs with particle j
1973 for(Int_t pastpair=0; pastpair<numOtherPairs2[0][j]; pastpair++){
1974 Int_t locationOtherPair = fOtherPairLocation2[0][j]->At(pastpair);
1975 if(locationOtherPair < 0) continue;// no pair there
1976 Int_t indexOther1 = i;
1977 Int_t indexOther2 = fNormPairs[0][ locationOtherPair ].fIndex1;
1978
1979 // Both possible orderings of other indexes
1980 if( (fNormPairSwitch[0][indexOther1]->At(indexOther2)=='1') || (fNormPairSwitch[0][indexOther2]->At(indexOther1)=='1')) {
1981
1982 // 1 and 2 are from SE
1983 ch3 = Int_t((fNormPairs[0][ locationOtherPair ].fCharge1 + 1)/2.);
1984 key3 = fNormPairs[0][ locationOtherPair ].fKey1;
1985 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1986 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
1987
1988 tempNormFillCount[bin1][bin2][bin3][fillIndex3][0]++;
1989 }
1990
1991 }// pastpair P11 loop
1992
1993
1994 fNormPairSwitch[en2][i]->AddAt('1',j);
1995 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
1996 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
1997
1998 numOtherPairs1[en2][i]++;
1999 numOtherPairs2[en2][j]++;
2000
2001
2002 normPairCount[en2]++;
2003 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
2004
2005 }// Norm Region
2006
2007 }// j particle
2008 }// i particle
2009
2010
2011
2012 //////////////////////////////////////////////
2013 // P12 pairing
2014 // 1st Particle
2015 for (Int_t i=0; i<myTracks; i++) {
2016
2017 Int_t en2=1;
2018
2019 // 2nd particle
2020 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2021
2022 key1 = (fEvt)->fTracks[i].fKey;
2023 key2 = (fEvt+en2)->fTracks[j].fKey;
2024 Short_t fillIndex2 = FillIndex2part(key1+key2);
2025 Short_t qCutBin = SetQcutBin(fillIndex2);
2026 Short_t normBin = SetNormBin(fillIndex2);
2027 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2028 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2029 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2030 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2031
2032 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2033 GetQosl(pVect1, pVect2, qout, qside, qlong);
2034 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2035
2036 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2037
2038 ///////////////////////////////
2039 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2040 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2041 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
2042
2043 if(fMCcase && ch1==ch2 && fMbin==0){
2044 for(Int_t rstep=0; rstep<10; rstep++){
2045 Float_t coeff = (rstep)*0.2*(0.18/1.2);
2046 // propagate through B field to r=1.2m
2047 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2048 if(phi1 > 2*PI) phi1 -= 2*PI;
2049 if(phi1 < 0) phi1 += 2*PI;
2050 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2051 if(phi2 > 2*PI) phi2 -= 2*PI;
2052 if(phi2 < 0) phi2 += 2*PI;
2053 Float_t deltaphi = phi1 - phi2;
2054 if(deltaphi > PI) deltaphi -= PI;
2055 if(deltaphi < -PI) deltaphi += PI;
2056 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2057 }
2058 }
2059
2060 if(ch1 == ch2){
2061 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2062 fPairSplitCut[1][i]->AddAt('1',j);
2063 continue;
2064 }
2065 }
2066
2067 //////////////////////////////////////////
2068 // 2-particle term
2069 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
5e3e77d6 2070 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
2071
cd12341d 2072 // osl frame
fa109294 2073 if(fillIndex2==0){
2074 if((transK12 > 0.2) && (transK12 < 0.3)){
cd12341d 2075 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2076 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2077 }
2078 if((transK12 > 0.6) && (transK12 < 0.7)){
2079 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2080 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2081 }
2082 }
2083 //////////////////////////////////////////
2084 if(fTabulatePairs){
2085 if(fillIndex2==0 && bin1==bin2){
2086 rapK12 = 0;
2087 transKbin=-1; rapKbin=-1;
2088 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
2089 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
2090 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2091 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2092 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2093 continue;
2094 }
2095 }
2096
2097
2098 if(pairCountME >= 2*kPairLimit) {cout<<"Too many ME pairs"<<endl; exitCode=kTRUE; continue;}
2099 if(exitCode) continue;
2100
2101 if(qinv12 <= fQcut[qCutBin]) {
2102 ///////////////////////////
2103
2104 // particle 1
2105 (fEvt)->fPairsME[pairCountME].fP1[0] = (fEvt)->fTracks[i].fP[0];
2106 (fEvt)->fPairsME[pairCountME].fP1[1] = (fEvt)->fTracks[i].fP[1];
2107 (fEvt)->fPairsME[pairCountME].fP1[2] = (fEvt)->fTracks[i].fP[2];
2108 (fEvt)->fPairsME[pairCountME].fE1 = (fEvt)->fTracks[i].fEaccepted;
2109 (fEvt)->fPairsME[pairCountME].fCharge1 = (fEvt)->fTracks[i].fCharge;
2110 (fEvt)->fPairsME[pairCountME].fIndex1 = i;
2111 (fEvt)->fPairsME[pairCountME].fKey1 = key1;
2112 (fEvt)->fPairsME[pairCountME].fLabel1 = (fEvt)->fTracks[i].fLabel;
2113 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
2114 (fEvt)->fPairsME[pairCountME].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
2115 (fEvt)->fPairsME[pairCountME].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
2116 (fEvt)->fPairsME[pairCountME].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
2117 }
2118 // particle 2
2119 (fEvt)->fPairsME[pairCountME].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
2120 (fEvt)->fPairsME[pairCountME].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
2121 (fEvt)->fPairsME[pairCountME].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
2122 (fEvt)->fPairsME[pairCountME].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
2123 (fEvt)->fPairsME[pairCountME].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2124 (fEvt)->fPairsME[pairCountME].fIndex2 = j;
2125 (fEvt)->fPairsME[pairCountME].fKey2 = key2;
2126 (fEvt)->fPairsME[pairCountME].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
2127 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
2128 (fEvt)->fPairsME[pairCountME].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2129 (fEvt)->fPairsME[pairCountME].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2130 (fEvt)->fPairsME[pairCountME].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2131 }
2132
2133 (fEvt)->fPairsME[pairCountME].fQinv = qinv12;
2134
2135 fPairLocationME[i]->AddAt(Int_t(pairCountME),j);
2136
2137 pairCountME++;
2138
2139 }
2140
2141 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2142 // particle 1
2143 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
2144 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
2145 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
2146 // particle 2
2147 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2148 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
2149 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
2150
2151 //other past pairs in P11 with particle i
2152 for(Int_t pastpairP11=0; pastpairP11<numOtherPairs2[0][i]; pastpairP11++){// past pair in P11 with i as 1st and 2nd particle
2153 Int_t locationOtherPairP11 = fOtherPairLocation2[0][i]->At(pastpairP11);// i is 2nd particle
2154 if(locationOtherPairP11 < 0) continue;// no pair there
2155 Int_t indexOther1P11 = fNormPairs[0][ locationOtherPairP11 ].fIndex1;
2156
2157 //Check other past pairs in P12
2158 if( (fNormPairSwitch[1][indexOther1P11]->At(j)=='0')) continue;
2159
2160 // 1 and 3 are from SE
2161 ch3 = Int_t((fNormPairs[0][ locationOtherPairP11 ].fCharge1 + 1)/2.);// charge of second particle in P11
2162 key3 = fNormPairs[0][ locationOtherPairP11 ].fKey1;
2163 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2164 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2165 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 2, bin1, bin2, bin3, fill2, fill3, fill4);
2166
2167
2168 if(fill2) tempNormFillCount[bin1][bin2][bin3][fillIndex3][1]++;
2169 if(fill3) tempNormFillCount[bin1][bin2][bin3][fillIndex3][2]++;
2170 if(fill4) tempNormFillCount[bin1][bin2][bin3][fillIndex3][3]++;
2171
2172
2173 }// P11 loop
2174
2175
2176 fNormPairSwitch[en2][i]->AddAt('1',j);
2177 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
2178 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
2179
2180 numOtherPairs1[en2][i]++;
2181 numOtherPairs2[en2][j]++;
2182
2183 normPairCount[en2]++;
2184 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
2185
2186 }// Norm Region
2187
2188
2189 }
2190 }
2191
2192
2193 ///////////////////////////////////////
2194 // P13 pairing (just for Norm counting of term5)
2195 for (Int_t i=0; i<myTracks; i++) {
2196
2197 // exit out of loop if there are too many pairs
2198 // dont bother with this loop if exitCode is set.
2199 if(exitCode) break;
2200
2201 // 2nd particle
2202 Int_t en2=2;
2203
2204 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2205
2206 key1 = (fEvt)->fTracks[i].fKey;
2207 key2 = (fEvt+en2)->fTracks[j].fKey;
2208 Short_t fillIndex2 = FillIndex2part(key1+key2);
2209 Short_t normBin = SetNormBin(fillIndex2);
2210 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2211 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2212 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2213 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2214
2215 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2216
2217 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2218
2219 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2220 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2221
2222 if(ch1 == ch2){
2223 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2224 fPairSplitCut[2][i]->AddAt('1',j);
2225 continue;
2226 }
2227 }
2228
2229 /////////////////////////////////////////////////////////
2230 // Normalization Region
2231
2232 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2233
2234 fNormPairSwitch[en2][i]->AddAt('1',j);
2235
2236 }// Norm Region
2237 }
2238 }
2239
2240
2241
2242 ///////////////////////////////////////
2243 // P23 pairing (just for Norm counting of term5)
2244 Int_t en1=1;
2245 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {
2246
2247 // exit out of loop if there are too many pairs
2248 // dont bother with this loop if exitCode is set.
2249 if(exitCode) break;
2250
2251 // 2nd event
2252 Int_t en2=2;
2253 // 2nd particle
2254 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2255
2256 if(exitCode) break;
2257
2258 key1 = (fEvt+en1)->fTracks[i].fKey;
2259 key2 = (fEvt+en2)->fTracks[j].fKey;
2260 Short_t fillIndex2 = FillIndex2part(key1+key2);
2261 Short_t normBin = SetNormBin(fillIndex2);
2262 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2263 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2264 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2265 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2266
2267 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2268
2269 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2270
2271 ///////////////////////////////
2272 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2273 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2274
2275 if(ch1 == ch2){
2276 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2277 fPairSplitCut[3][i]->AddAt('1',j);
2278 continue;
2279 }
2280 }
2281
2282 if((qinv12 < fNormQcutLow[normBin]) || (qinv12 >= fNormQcutHigh[normBin])) continue;
2283
2284 Int_t index1P23 = i;
2285 Int_t index2P23 = j;
2286
2287 for(Int_t pastpairP12=0; pastpairP12<numOtherPairs2[1][index1P23]; pastpairP12++){// loop in P12 with i as 2nd particle
2288 Int_t locationOtherPairP12 = fOtherPairLocation2[1][index1P23]->At(pastpairP12);
2289 if(locationOtherPairP12 < 0) continue; // no pair there
2290 Int_t index1P12 = fNormPairs[1][ locationOtherPairP12 ].fIndex1;
2291
2292
2293 //Check other past pair status in P13
2294 if( (fNormPairSwitch[2][index1P12]->At(index2P23)=='0')) continue;
2295
2296 // all from different event
2297 ch3 = Int_t((fNormPairs[1][ locationOtherPairP12 ].fCharge1 + 1)/2.);// charge of first particle in P12
2298 key3 = fNormPairs[1][ locationOtherPairP12 ].fKey1;
2299 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2300 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2301
2302 tempNormFillCount[bin1][bin2][bin3][fillIndex3][4]++;
2303 }
2304 }
2305 }
2306
2307
2308
2309
2310 ///////////////////////////////////////////////////
2311 // Do not use pairs from events with too many pairs
2312 if(exitCode) {
2313 cout<<"SE or ME or Norm PairCount too large. Discarding all pairs and skipping event"<<endl;
2314 (fEvt)->fNpairsSE = 0;
2315 (fEvt)->fNpairsME = 0;
2316 ((TH1F*)fOutputList->FindObject("fRejectedEvents"))->Fill(fMbin+1);
2317 return;// Skip event
2318 }else{
2319 (fEvt)->fNpairsSE = pairCountSE;
2320 (fEvt)->fNpairsME = pairCountME;
2321 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2322 }
2323 ///////////////////////////////////////////////////
2324
cd12341d 2325
2326 //cout<<"pairCountSE = "<<pairCountSE<<" pairCountME = "<<pairCountME<<endl;
b6e5ec54 2327 //cout<<"Start Main analysis"<<endl;
cd12341d 2328
2329 ///////////////////////////////////////////////////////////////////////
2330 ///////////////////////////////////////////////////////////////////////
2331 ///////////////////////////////////////////////////////////////////////
2332 //
2333 //
5e3e77d6 2334 // Start the Main Correlation Analysis
cd12341d 2335 //
2336 //
2337 ///////////////////////////////////////////////////////////////////////
2338
2339
2340 /////////////////////////////////////////////////////////
2341 // Skip 3-particle part if Tabulate6DPairs is set to true
2342 if(fTabulatePairs) return;
2343 /////////////////////////////////////////////////////////
2344
2345 // Set the Normalization counters
2346 for(Int_t termN=0; termN<5; termN++){
2347
2348 if(termN==0){
2349 if((fEvt)->fNtracks ==0) continue;
2350 }else if(termN<4){
2351 if((fEvt)->fNtracks ==0) continue;
2352 if((fEvt+1)->fNtracks ==0) continue;
2353 }else {
2354 if((fEvt)->fNtracks ==0) continue;
2355 if((fEvt+1)->fNtracks ==0) continue;
2356 if((fEvt+2)->fNtracks ==0) continue;
2357 }
2358
2359 for(Int_t sc=0; sc<kSCLimit3; sc++){
2360
2361 for(Int_t c1=0; c1<2; c1++){
2362 for(Int_t c2=0; c2<2; c2++){
2363 for(Int_t c3=0; c3<2; c3++){
2364
2365 if(sc==0 || sc==6 || sc==9){// Identical species
2366 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
2367 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
2368 }else if(sc!=5){
2369 if( (c1+c2)==1) {if(c1!=0) continue;}
2370 }else {}// do nothing for pi-k-p case
2371 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[fMbin].EDB[fEDbin].ThreePT[termN].fNorm3->Fill(0.,tempNormFillCount[c1][c2][c3][sc][termN]);
2372 }
2373 }
2374 }
2375 }
2376 }
2377
2378
2379
2380 /////////////////////////////////////////////
2381 // Calculate Pair-Cut Correlations
2382 for(Int_t en1case=0; en1case<2; en1case++){// limit at 2 (normal)
2383
2384 Int_t nump1=0;
2385 if(en1case==0) nump1 = (fEvt)->fNpairsSE;
2386 if(en1case==1) nump1 = (fEvt)->fNpairsME;
2387
2388 // 1st pair
2389 for(Int_t p1=0; p1<nump1; p1++){
2390
2391 if(en1case==0){
2392 ch1 = Int_t(((fEvt)->fPairsSE[p1].fCharge1 + 1)/2.);
2393 ch2 = Int_t(((fEvt)->fPairsSE[p1].fCharge2 + 1)/2.);
2394 pVect1[0] = (fEvt)->fPairsSE[p1].fE1; pVect2[0] = (fEvt)->fPairsSE[p1].fE2;
2395 pVect1[1] = (fEvt)->fPairsSE[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsSE[p1].fP2[0];
2396 pVect1[2] = (fEvt)->fPairsSE[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsSE[p1].fP2[1];
2397 pVect1[3] = (fEvt)->fPairsSE[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsSE[p1].fP2[2];
2398 index1 = (fEvt)->fPairsSE[p1].fIndex1; index2 = (fEvt)->fPairsSE[p1].fIndex2;
2399 key1 = (fEvt)->fPairsSE[p1].fKey1; key2 = (fEvt)->fPairsSE[p1].fKey2;
2400 pVect1MC[1] = (fEvt)->fPairsSE[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsSE[p1].fP2MC[0];
2401 pVect1MC[2] = (fEvt)->fPairsSE[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsSE[p1].fP2MC[1];
2402 pVect1MC[3] = (fEvt)->fPairsSE[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsSE[p1].fP2MC[2];
2403 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2404 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2405
2406 qinv12 = (fEvt)->fPairsSE[p1].fQinv;
2407 }
2408 if(en1case==1){
2409 ch1 = Int_t(((fEvt)->fPairsME[p1].fCharge1 + 1)/2.);
2410 ch2 = Int_t(((fEvt)->fPairsME[p1].fCharge2 + 1)/2.);
2411 pVect1[0] = (fEvt)->fPairsME[p1].fE1; pVect2[0] = (fEvt)->fPairsME[p1].fE2;
2412 pVect1[1] = (fEvt)->fPairsME[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsME[p1].fP2[0];
2413 pVect1[2] = (fEvt)->fPairsME[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsME[p1].fP2[1];
2414 pVect1[3] = (fEvt)->fPairsME[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsME[p1].fP2[2];
2415 index1 = (fEvt)->fPairsME[p1].fIndex1; index2 = (fEvt)->fPairsME[p1].fIndex2;
2416 key1 = (fEvt)->fPairsME[p1].fKey1; key2 = (fEvt)->fPairsME[p1].fKey2;
2417 pVect1MC[1] = (fEvt)->fPairsME[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsME[p1].fP2MC[0];
2418 pVect1MC[2] = (fEvt)->fPairsME[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsME[p1].fP2MC[1];
2419 pVect1MC[3] = (fEvt)->fPairsME[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsME[p1].fP2MC[2];
2420 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2421 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2422
2423 qinv12 = (fEvt)->fPairsME[p1].fQinv;
2424 }
2425
cd12341d 2426
2427 // en2 buffer
2428 for(Int_t en2=0; en2<3; en2++){
2429 //////////////////////////////////////
2430
2431 Bool_t skipcase=kTRUE;
2432 Short_t config=-1, part=-1;
2433 if(en1case==0 && en2==0) {skipcase=kFALSE; config=1; part=0;}// P11T1
2434 if(en1case==0 && en2==1) {skipcase=kFALSE; config=2; part=1;}// P11T2
2435 if(en1case==1 && en2==0) {skipcase=kFALSE; config=2; part=2;}// P12T1
2436 if(en1case==1 && en2==2) {skipcase=kFALSE; config=3; part=3;}// P12T3
2437
2438 if(skipcase) continue;
2439
2440
2441 // 3-particle terms
2442 // 3rd particle
2443 for(Int_t k=0; k<(fEvt+en2)->fNtracks; k++){
2444 index3 = k;
2445
2446
2447 // remove auto-correlations and duplicate triplets
2448 if(config==1){
2449 if( index1 == index3) continue;
2450 if( index2 == index3) continue;
2451 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2452
2453 // skip the switched off triplets
2454 if(fTripletSkip1[fPairLocationSE[index1]->At(index2)]->At(index3)=='1') {
2455 fTripletSkip1[fPairLocationSE[index1]->At(index2)]->AddAt('0',index3);// Reset
2456 continue;
2457 }
2458 ///////////////////////////////
2459 // Turn off 1st possible degenerate triplet
2460 if(index1 < index3){// verify correct id ordering ( index1 < k )
2461 if(fPairLocationSE[index1]->At(index3) >= 0){
2462 fTripletSkip1[fPairLocationSE[index1]->At(index3)]->AddAt('1',index2);
2463 }
2464 if(fPairSplitCut[0][index1]->At(index3)=='1') continue;// Track splitting/merging
2465 }else {// or k < index1
2466 if(fPairLocationSE[index3]->At(index1) >= 0){
2467 fTripletSkip1[fPairLocationSE[index3]->At(index1)]->AddAt('1',index2);
2468 }
2469 if(fPairSplitCut[0][index3]->At(index1)=='1') continue;// Track splitting/merging
2470 }
2471 // turn off 2nd possible degenerate triplet
2472 if(index2 < index3){// verify correct id ordering (index2 < k)
2473 if(fPairLocationSE[index2]->At(index3) >= 0){
2474 fTripletSkip1[fPairLocationSE[index2]->At(index3)]->AddAt('1',index1);
2475 }
2476 if(fPairSplitCut[0][index2]->At(index3)=='1') continue;// Track splitting/merging
2477 }else {// or k < index2
2478 if(fPairLocationSE[index3]->At(index2) >= 0){
2479 fTripletSkip1[fPairLocationSE[index3]->At(index2)]->AddAt('1',index1);
2480 }
2481 if(fPairSplitCut[0][index3]->At(index2)=='1') continue;// Track splitting/merging
2482 }
2483
2484 }// end config 1
2485
2486 if(config==2 && part==1){// SE pair and third particle from next event. P11T2
2487 ///////////////////////////////
2488 // Turn off 1st possible degenerate triplet
2489 if(fPairLocationME[index1]->At(index3) >= 0){
2490 fTripletSkip2[fPairLocationME[index1]->At(index3)]->AddAt('1',index2);
2491 }
2492
2493 // turn off 2nd possible degenerate triplet
2494 if(fPairLocationME[index2]->At(index3) >= 0){
2495 fTripletSkip2[fPairLocationME[index2]->At(index3)]->AddAt('1',index1);
2496 }
2497
2498 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2499 if(fPairSplitCut[1][index1]->At(index3)=='1') continue;// Track splitting/merging
2500 if(fPairSplitCut[1][index2]->At(index3)=='1') continue;// Track splitting/merging
2501 }// end config 2 part 1
2502
2503 if(config==2 && part==2){// P12T1
2504 if( index1 == index3) continue;
2505
2506 // skip the switched off triplets
2507 if(fTripletSkip2[fPairLocationME[index1]->At(index2)]->At(index3)=='1') {
2508 fTripletSkip2[fPairLocationME[index1]->At(index2)]->AddAt('0',index3);// Reset
2509 continue;
2510 }
2511 // turn off another possible degenerate
2512 if(fPairLocationME[index3]->At(index2) >= 0){
2513 fTripletSkip2[fPairLocationME[index3]->At(index2)]->AddAt('1',index1);
2514 }// end config 2 part 2
2515
2516 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2517 if(index1 < index3) {if(fPairSplitCut[0][index1]->At(index3)=='1') continue;}// Track splitting/merging
2518 else {if(fPairSplitCut[0][index3]->At(index1)=='1') continue;}// Track splitting/merging
2519 if(fPairSplitCut[1][index3]->At(index2)=='1') continue;// Track splitting/merging
2520 }
2521 if(config==3){// P12T3
2522 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2523 if(fPairSplitCut[2][index1]->At(index3)=='1') continue;// Track splitting/merging
2524 if(fPairSplitCut[3][index2]->At(index3)=='1') continue;// Track splitting/merging
2525 }// end config 3
2526
2527
5e3e77d6 2528
cd12341d 2529 ch3 = Int_t(((fEvt+en2)->fTracks[k].fCharge + 1)/2.);
2530 key3 = (fEvt+en2)->fTracks[k].fKey;
2531 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2532 Short_t fillIndex13 = FillIndex2part(key1+key3);
2533 Short_t fillIndex23 = FillIndex2part(key2+key3);
2534 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2535 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2536 pVect3[0] = (fEvt+en2)->fTracks[k].fEaccepted;
2537 pVect3[1] = (fEvt+en2)->fTracks[k].fP[0];
2538 pVect3[2] = (fEvt+en2)->fTracks[k].fP[1];
2539 pVect3[3] = (fEvt+en2)->fTracks[k].fP[2];
5e3e77d6 2540 if(fMCcase){
2541 pVect3MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPx;
2542 pVect3MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPy;
2543 pVect3MC[3] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPz;
2544 pVect3MC[0] = sqrt(pow(pVect3MC[1],2)+pow(pVect3MC[2],2)+pow(pVect3MC[3],2)+pow(fTrueMassPi,2));
2545 qinv12MC = GetQinv(0, pVect1MC, pVect2MC);
2546 qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
2547 qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
2548 }
cd12341d 2549 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2550 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2551
5e3e77d6 2552
cd12341d 2553 if(qinv13 < fQLowerCut) continue;
2554 if(qinv23 < fQLowerCut) continue;
5e3e77d6 2555 if(qinv13 > fQcut[qCutBin13]) continue;// cut value was 3x higher before
2556 if(qinv23 > fQcut[qCutBin23]) continue;// cut value was 3x higher before
2557 // if all three pair cuts are the same then the case (config=2 && term=2) never reaches here.
2558
cd12341d 2559 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2560 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2561 Float_t firstQ=0, secondQ=0, thirdQ=0;
5e3e77d6 2562 //
90814457 2563
2564
2565
fa109294 2566
5e3e77d6 2567 //
cd12341d 2568 if(config==1) {// 123
2569 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
fa109294 2570
cd12341d 2571 if(fillIndex3 <= 2){
2572 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 0, 1, firstQ, secondQ, thirdQ);
2573 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ);
2574 ((TH2F*)fOutputList->FindObject("fKt3Dist"))->Fill(fMbin+1, transK3);
5e3e77d6 2575 //
2576 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
90814457 2577 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2578 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2579 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
5e3e77d6 2580 }
2581 //////////////////////////////////////
90814457 2582 // Momentum resolution and <K3> calculation
5e3e77d6 2583 if(fillIndex3==0 && fMCcase){
2584 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2585 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, 0, 1, firstQMC, secondQMC, thirdQMC);
54d66278 2586 Float_t WInput = 1;
90814457 2587 Double_t K3=1.0;
5e3e77d6 2588 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2589 WInput = MCWeight3D(kTRUE, 1, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2590 K3 = FSICorrelationOmega0(kTRUE, firstQMC, secondQMC, thirdQMC);// K3
5e3e77d6 2591 }else {// mixed charge
90814457 2592 if(bin1==bin2) WInput = MCWeight3D(kFALSE, 1, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2593 else WInput = MCWeight3D(kFALSE, 1, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2594 K3 = FSICorrelationOmega0(kFALSE, firstQMC, secondQMC, thirdQMC);// K3
5e3e77d6 2595 }
54d66278 2596 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2597 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2598 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fQW12->Fill(firstQMC, secondQMC, thirdQMC, WInput*firstQMC);
2599 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fQW13->Fill(firstQMC, secondQMC, thirdQMC, WInput*secondQMC);
90814457 2600 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSumK3->Fill(firstQMC, secondQMC, thirdQMC, WInput/K3);
2601 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fEnK3->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2602 if(ch1==ch2 && ch1==ch3){
2603 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2604 FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2605 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2606 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2607 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2608 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2609 //
2610 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsSumK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K3);
2611 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsSumK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K3);
2612 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsEnK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2613 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsEnK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2614 }
2615
fa109294 2616 }// fMCcase
2617
cd12341d 2618 }
2619
2620 }else if(config==2){// 12, 13, 23
5e3e77d6 2621
cd12341d 2622 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2623 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
2624
2625 // loop over terms 2-4
5e3e77d6 2626 for(Int_t jj=2; jj<5; jj++){
2627 if(jj==2) {if(!fill2) continue;}//12
2628 if(jj==3) {if(!fill3) continue;}//13
2629 if(jj==4) {if(!fill4) continue;}//23
cd12341d 2630
2631 if(fillIndex3 <= 2){
5e3e77d6 2632 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj, firstQ, secondQ, thirdQ);
2633 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTerms3->Fill(firstQ, secondQ, thirdQ);
90814457 2634 if(fillIndex3==0 && ch1==ch2 && ch1==ch3){
2635 if(part==1){// P11T2
2636 if(jj==2) {
2637 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2638 if(fMCcase) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2639 }if(jj==3){
2640 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2641 if(fMCcase) FourVectProdTerms(pVect1MC, pVect3MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2642 }if(jj==4) {
2643 FourVectProdTerms(pVect3, pVect1, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2644 if(fMCcase) FourVectProdTerms(pVect3MC, pVect1MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2645 }
2646 }else{// P12T1
2647 if(jj==2) {
2648 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2649 if(fMCcase) FourVectProdTerms(pVect1MC, pVect3MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2650 }if(jj==3) {
2651 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2652 if(fMCcase) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2653 }if(jj==4) {
2654 FourVectProdTerms(pVect2, pVect1, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2655 if(fMCcase) FourVectProdTerms(pVect2MC, pVect1MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2656 }
2657 }
2658 if(!fMCcase){
2659 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2660 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
fa109294 2661 }
5e3e77d6 2662 }
2663 //////////////////////////////////////
2664 // Momentum resolution calculation
2665 if(fillIndex3==0 && fMCcase){
2666 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2667 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, jj, firstQMC, secondQMC, thirdQMC);
90814457 2668 Float_t WInput = 1;
5e3e77d6 2669 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2670 WInput = MCWeight3D(kTRUE, jj, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
5e3e77d6 2671 }else {// mixed charge
90814457 2672 if(bin1==bin2) WInput = MCWeight3D(kFALSE, jj, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2673 else WInput = MCWeight3D(kFALSE, 6-jj, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2674 }
2675 //
2676 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2677 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2678 //
2679 if(ch1==ch2 && ch1==ch3){
2680 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2681 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2682 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2683 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2684 //
2685 Float_t InteractingQ=0;
2686 if(part==1) {InteractingQ=qinv12;}// 12 from SE
2687 else {InteractingQ=qinv13;}// 13 from SE
2688 Double_t K2 = FSICorrelationTherm2(+1,+1, InteractingQ);// K2 from Therminator source
2689 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSumK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K2);
2690 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSumK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K2);
2691 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsEnK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2692 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsEnK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
5e3e77d6 2693 }
2694 }// fMCcase
2695
cd12341d 2696 }
2697 }
5e3e77d6 2698
cd12341d 2699 }else {// config 3: All particles from different events
2700 //Simulate Filling of other event-mixed lowQ pairs
2701
2702 Float_t enhancement=1.0;
2703 Int_t nUnderCut=0;
2704 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2705 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2706
2707 if(nUnderCut==0) enhancement = (1+1+1)/1.;// 1 LowQ pair
2708 if(nUnderCut==1) enhancement = (1+2)/2.;// 2 LowQ pair
2709 if(nUnderCut==2) enhancement = 1.;// 3 LowQ pair
2710
2711 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
90814457 2712
2713 if(ch1==ch2 && ch1==ch3 && fillIndex3==0) FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2714 if(fMCcase && ch1==ch2 && ch1==ch3 && fillIndex3==0) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);
cd12341d 2715
2716 if(fillIndex3 <= 2){
2717 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
2718 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ, enhancement);
5e3e77d6 2719 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
90814457 2720 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2721 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
5e3e77d6 2722 }
2723 //////////////////////////////////////
2724 // Momentum resolution calculation
2725 if(fillIndex3==0 && fMCcase){
2726 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2727 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, 5, firstQMC, secondQMC, thirdQMC);
90814457 2728 Float_t WInput=1;
5e3e77d6 2729 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2730 WInput = MCWeight3D(kTRUE, 5, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
5e3e77d6 2731 }else {// mixed charge
90814457 2732 if(bin1==bin2) WInput = MCWeight3D(kFALSE, 5, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2733 else WInput = MCWeight3D(kFALSE, 5, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss in this case. 1st Q argument is ss
2734 }
2735 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2736 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2737 if(ch1==ch2 && ch1==ch3){
2738 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2739 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2740 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2741 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2742
5e3e77d6 2743 }
2744 }// fMCcase
2745
cd12341d 2746 }
5e3e77d6 2747
cd12341d 2748 if(fillIndex3 !=0) continue;// only calculate TPN for pi-pi-pi
2749 if(ch1!=ch2 || ch1!=ch3) continue;// only calcualte TPN for ss
5e3e77d6 2750
90814457 2751 //if(fMCcase) continue;// only calcualte TPN for real data
2752
cd12341d 2753 GetWeight(pVect1, pVect2, weight12, weight12Err);
2754 GetWeight(pVect1, pVect3, weight13, weight13Err);
2755 GetWeight(pVect2, pVect3, weight23, weight23Err);
2756 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) continue;// weight should never be larger than 1
2757
5e3e77d6 2758
2759 // Coul correlations from Wave-functions
fa109294 2760 //for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 8fm, last value for Therminator
2761 //for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){// 0.3 to 0.6
2762 //Float_t myDamp = fDampStart + (fDampStep)*myDampIt;
2763 //Int_t denIndex = (kRVALUES-1)*kNDampValues + myDampIt;
2764 //Int_t denIndex = myDampIt;
2765 Int_t myDampIt = 5;
2766 Float_t myDamp = 0.4;
2767 Int_t denIndex = 0;
90814457 2768 Int_t momResIndex = rIndexForTPN*kNDampValues + myDampIt;// Therminator slot uses R=7 for momentum resolution
2769
2770 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, qinv12);
2771 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, qinv13);
2772 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, qinv23);
cd12341d 2773
90814457 2774 Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
2775 if(!fMCcase) {
2776 Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
2777 Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
2778 Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);
2779 if(momBin12 >= kQbins) momBin12 = kQbins-1;
2780 if(momBin13 >= kQbins) momBin13 = kQbins-1;
2781 if(momBin23 >= kQbins) momBin23 = kQbins-1;
2782 MomResCorr12 = fMomResC2->GetBinContent(momResIndex+1, momBin12);
2783 MomResCorr13 = fMomResC2->GetBinContent(momResIndex+1, momBin13);
2784 MomResCorr23 = fMomResC2->GetBinContent(momResIndex+1, momBin23);
2785 }
2786 weight12CC = ((weight12+1)*MomResCorr12 - myDamp*coulCorr12 - (1-myDamp));
5e3e77d6 2787 weight12CC /= coulCorr12*myDamp;
90814457 2788 weight13CC = ((weight13+1)*MomResCorr13 - myDamp*coulCorr13 - (1-myDamp));
5e3e77d6 2789 weight13CC /= coulCorr13*myDamp;
90814457 2790 weight23CC = ((weight23+1)*MomResCorr23 - myDamp*coulCorr23 - (1-myDamp));
5e3e77d6 2791 weight23CC /= coulCorr23*myDamp;
cd12341d 2792
5e3e77d6 2793 if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {// testing purposes only
cd12341d 2794 if(denIndex==71 && fMbin==0 && ch1==-1) {
2795 weightTotal = sqrt(fabs(weight12CC*weight13CC*weight23CC));
2796 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
2797 }
2798 continue;// C2^QS can never be less than unity
2799 }
2800
2801 weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
90814457 2802
cd12341d 2803 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].fTwoPartNorm->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
90814457 2804
2805 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNorm->Fill(Qsum1v1, Qsum2, Qsum3v1, enhancement*weightTotal);
2806 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNorm->Fill(Qsum1v2, Qsum2, Qsum3v2, enhancement*weightTotal);
2807 if(fMCcase){
2808 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNormIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, enhancement*weightTotal);
2809 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNormSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, enhancement*weightTotal);
2810 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNormIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, enhancement*weightTotal);
2811 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNormSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, enhancement*weightTotal);
fa109294 2812 }
90814457 2813
2814
cd12341d 2815 // Save cpu time and memory by skipping r3 denominator calculation below. den errors are negligible compared to num errors.
2816 /*
2817 if(weightTotal > 0.0001){// tiny numbers cause a Float_ting point exception below
2818 weightTotalErr = pow((weight12Err*coulCorr12)*weight13CC*weight23CC,2);
2819 weightTotalErr += pow(weight12CC*(weight13Err*coulCorr13)*weight23CC,2);
2820 weightTotalErr += pow(weight12CC*weight13CC*(weight23Err*coulCorr23),2);
2821 weightTotalErr /= pow(2*weightTotal,2);
2822
2823 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].TwoPartNormErr->Fill(denIndex, q3, enhancement*weightTotalErr);
2824 }
2825 */
2826
fa109294 2827 //}// damp iter
2828 //}// R iter
cd12341d 2829
2830
2831
2832 }
2833 }// end 3rd particle
2834 }// en2
2835
2836
fa109294 2837 }// p1
2838 }//en1
cd12341d 2839
cd12341d 2840 ///////////////////
2841 }// end of PdensityPairs
2842
fa109294 2843
cd12341d 2844
2845
2846
2847
2848
2849 ////////////////////////////////////////////////////////
2850 // Pdensity Method with Explicit Loops
2851 if(fPdensityExplicitLoop){
2852
2853 ////////////////////////////////////
2854 // 2nd, 3rd, and 4th order Correlations
2855
2856 // First Particle
2857 for (Int_t i=0; i<myTracks; i++) {
2858 ch1 = Int_t( ((fEvt)->fTracks[i].fCharge + 1)/2. );
2859 pVect1[0] = (fEvt)->fTracks[i].fEaccepted;
2860 pVect1[1] = (fEvt)->fTracks[i].fP[0];
2861 pVect1[2] = (fEvt)->fTracks[i].fP[1];
2862 pVect1[3] = (fEvt)->fTracks[i].fP[2];
2863 key1 = (fEvt)->fTracks[i].fKey;
2864
2865 // Second Event
2866 for(Int_t en2=0; en2<fEventsToMix+1; en2++){
2867 Int_t startbin2=0;
2868 if(en2==0) startbin2=i+1;
2869
2870 // Second Particle
2871 for (Int_t j=startbin2; j<(fEvt+en2)->fNtracks; j++) {
2872 ch2 = Int_t( ((fEvt+en2)->fTracks[j].fCharge + 1)/2. );
2873 pVect2[0] = (fEvt+en2)->fTracks[j].fEaccepted;
2874 pVect2[1] = (fEvt+en2)->fTracks[j].fP[0];
2875 pVect2[2] = (fEvt+en2)->fTracks[j].fP[1];
2876 pVect2[3] = (fEvt+en2)->fTracks[j].fP[2];
2877 key2 = (fEvt+en2)->fTracks[j].fKey;
2878
2879 Short_t fillIndex12 = FillIndex2part(key1+key2);
2880 qinv12 = GetQinv(fillIndex12, pVect1, pVect2);
2881
2882 if(qinv12 < fQLowerCut) continue;
2883
2884
2885 // 2-particle part is filled always during pair creator
2886
2887 // Third Event
2888 for(Int_t en3=en2; en3<fEventsToMix+1; en3++){
2889 Int_t startbin3=0;
2890 if(en3==en2) startbin3=j+1;
2891 else startbin3=0;
2892
2893
2894 // Third Particle
2895 for (Int_t k=startbin3; k<(fEvt+en3)->fNtracks; k++) {
2896 ch3 = Int_t( ((fEvt+en3)->fTracks[k].fCharge + 1)/2. );
2897 pVect3[0] = (fEvt+en3)->fTracks[k].fEaccepted;
2898 pVect3[1] = (fEvt+en3)->fTracks[k].fP[0];
2899 pVect3[2] = (fEvt+en3)->fTracks[k].fP[1];
2900 pVect3[3] = (fEvt+en3)->fTracks[k].fP[2];
2901 key3 = (fEvt+en3)->fTracks[k].fKey;
2902
2903 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2904 Short_t fillIndex13 = FillIndex2part(key1+key3);
2905 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2906 Short_t fillIndex23 = FillIndex2part(key2+key3);
2907 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2908
2909
2910 if(qinv13 < fQLowerCut) continue;
2911 if(qinv23 < fQLowerCut) continue;
2912
2913
2914 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2915
2916 Short_t normBin12 = SetNormBin(fillIndex12);
2917 Short_t normBin13 = SetNormBin(fillIndex13);
2918 Short_t normBin23 = SetNormBin(fillIndex23);
2919
2920
2921 if(en3==0 && en2==0) {// 123
2922 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2923
2924 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fExplicit3->Fill(q3);// 123
2925
2926 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2927 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2928 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fNormEx3->Fill(0.);
2929 }
2930 }
2931
2932 }else if((en2==0 && en3==1) ) {// 12-3, 13-2, 23-1
2933 Float_t gFact=1;
2934
2935 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2936 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2937
2938
2939 if(fill2){
2940 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fExplicit3->Fill(q3, gFact);// 12
2941 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2942 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2943 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fNormEx3->Fill(0.);
2944 }
2945 }
2946 }
2947 if(fill3){
2948 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fExplicit3->Fill(q3, gFact);// 12
2949 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2950 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2951 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fNormEx3->Fill(0.);
2952 }
2953 }
2954 }
2955 if(fill4){
2956 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fExplicit3->Fill(q3, gFact);// 12
2957 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2958 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2959 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fNormEx3->Fill(0.);
2960 }
2961 }
2962 }
2963
2964 }else if(en2==1 && en3==2){// all uncorrelated events
2965 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2966
2967 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fExplicit3->Fill(q3);
2968 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2969 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2970 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fNormEx3->Fill(0.);
2971 }
2972 }
2973 Short_t qCutBin12 = SetQcutBin(fillIndex12);
2974 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2975 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2976
2977 if( (qinv12 < fQcut[qCutBin12]) || (qinv13 < fQcut[qCutBin13]) || (qinv23 < fQcut[qCutBin23])){
2978
2979 Int_t nUnderCut=0;
2980 if(qinv12<fQcut[qCutBin12]) nUnderCut++;
2981 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2982 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2983
2984 }
2985
2986 }else {}
2987
2988
2989 }// 3rd particle
2990 }// 3rd event
2991
2992 }// 2nd particle
2993 }// 2nd event
2994
2995 }// 1st particle
2996
2997
2998
2999
3000 }// End of PdensityExplicit
3001
3002
3003
3004
3005 // Post output data.
3006 PostData(1, fOutputList);
3007
3008}
3009//________________________________________________________________________
3010void AliChaoticity::Terminate(Option_t *)
3011{
3012 // Called once at the end of the query
3013
3014 cout<<"Done"<<endl;
3015
3016}
3017//________________________________________________________________________
3018Bool_t AliChaoticity::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
3019{
3020
3021 if(fabs(first.fEta-second.fEta) > fMinSepTPCEntranceEta) return kTRUE;
3022
3023 // propagate through B field to r=1m
5e3e77d6 3024 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
cd12341d 3025 if(phi1 > 2*PI) phi1 -= 2*PI;
3026 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 3027 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
cd12341d 3028 if(phi2 > 2*PI) phi2 -= 2*PI;
3029 if(phi2 < 0) phi2 += 2*PI;
3030
3031 Float_t deltaphi = phi1 - phi2;
3032 if(deltaphi > PI) deltaphi -= 2*PI;
3033 if(deltaphi < -PI) deltaphi += 2*PI;
3034 deltaphi = fabs(deltaphi);
3035
5e3e77d6 3036 //cout<<deltaphi<<" "<<fMinSepTPCEntrancePhi<<" "<<fMinSepTPCEntranceEta<<endl;
3037 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 3038
3039 // propagate through B field to r=1.6m
5e3e77d6 3040 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
cd12341d 3041 if(phi1 > 2*PI) phi1 -= 2*PI;
3042 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 3043 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
cd12341d 3044 if(phi2 > 2*PI) phi2 -= 2*PI;
3045 if(phi2 < 0) phi2 += 2*PI;
3046
3047 deltaphi = phi1 - phi2;
3048 if(deltaphi > PI) deltaphi -= 2*PI;
3049 if(deltaphi < -PI) deltaphi += 2*PI;
3050 deltaphi = fabs(deltaphi);
3051
5e3e77d6 3052 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 3053
3054
3055 //
3056
3057 Int_t ncl1 = first.fClusterMap.GetNbits();
3058 Int_t ncl2 = second.fClusterMap.GetNbits();
3059 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
3060 Double_t shfrac = 0; Double_t qfactor = 0;
3061 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
3062 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
3063 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
3064 sumQ++;
3065 sumCls+=2;
3066 sumSha+=2;}
3067 else {sumQ--; sumCls+=2;}
3068 }
3069 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
3070 sumQ++;
3071 sumCls++;}
3072 }
3073 if (sumCls>0) {
3074 qfactor = sumQ*1.0/sumCls;
3075 shfrac = sumSha*1.0/sumCls;
3076 }
3077
3078 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
3079
3080
3081 return kTRUE;
5e3e77d6 3082
cd12341d 3083
3084}
3085//________________________________________________________________________
3086Float_t AliChaoticity::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
3087{
3088 Float_t arg = G_Coeff/qinv;
3089
3090 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
3091 else {return (exp(-arg)-1)/(-arg);}
3092
3093}
3094//________________________________________________________________________
3095void AliChaoticity::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
3096{
3097 Int_t j, k;
3098 Int_t a = i2 - i1;
3099 for (Int_t i = i1; i < i2+1; i++) {
3100 j = (Int_t) (gRandom->Rndm() * a);
3101 k = iarr[j];
3102 iarr[j] = iarr[i];
3103 iarr[i] = k;
3104 }
3105}
3106//________________________________________________________________________
3107Short_t AliChaoticity::FillIndex2part(Short_t key){
3108
3109 if(key==2) return 0;// pi-pi
3110 else if(key==11) return 1;// pi-k
3111 else if(key==101) return 2;// pi-p
3112 else if(key==20) return 3;// k-k
3113 else if(key==110) return 4;// k-p
3114 else return 5;// p-p
3115}
3116//________________________________________________________________________
3117Short_t AliChaoticity::FillIndex3part(Short_t key){
3118
3119 if(key==3) return 0;// pi-pi-pi
3120 else if(key==12) return 1;// pi-pi-k
3121 else if(key==21) return 2;// k-k-pi
3122 else if(key==102) return 3;// pi-pi-p
3123 else if(key==201) return 4;// p-p-pi
3124 else if(key==111) return 5;// pi-k-p
3125 else if(key==30) return 6;// k-k-k
3126 else if(key==120) return 7;// k-k-p
3127 else if(key==210) return 8;// p-p-k
3128 else return 9;// p-p-p
3129
3130}
3131//________________________________________________________________________
3132Short_t AliChaoticity::SetQcutBin(Short_t fi){// fi=FillIndex
3133 if(fi <= 2) return 0;
3134 else if(fi==3) return 1;
3135 else return 2;
3136}
3137//________________________________________________________________________
3138Short_t AliChaoticity::SetNormBin(Short_t fi){// fi=FillIndex
3139 if(fi==0) return 0;
3140 else if(fi <= 2) return 1;
3141 else return 2;
3142}
3143//________________________________________________________________________
3144void AliChaoticity::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
3145
3146 if(fi==0 || fi==3 || fi==5){// Identical species
3147 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
3148 else {b1=c1; b2=c2;}
3149 }else {// Mixed species
3150 if(key1 < key2) { b1=c1; b2=c2;}
3151 else {b1=c2; b2=c1;}
3152 }
3153
3154}
3155//________________________________________________________________________
3156void AliChaoticity::SetFillBins3(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Short_t part, Int_t &b1, Int_t &b2, Int_t &b3, Bool_t &fill2, Bool_t &fill3, Bool_t &fill4){
3157
3158
3159 // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
3160 Bool_t seSS=kFALSE;
3161 Bool_t seSK=kFALSE;
3162 Short_t seKeySum=0;// only used for pi-k-p case
3163 if(part==1) {// default case (irrelevant for term 1 and term 5)
3164 if(c1==c2) seSS=kTRUE;
3165 if(key1==key2) seSK=kTRUE;
3166 seKeySum = key1+key2;
3167 }
3168 if(part==2){
3169 if(c1==c3) seSS=kTRUE;
3170 if(key1==key3) seSK=kTRUE;
3171 seKeySum = key1+key3;
3172 }
3173
3174
3175 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3176
3177 if(fi==0 || fi==6 || fi==9){// Identical species
3178 if( (c1+c2+c3)==1) {
3179 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
3180 //
3181 if(seSS) fill2=kTRUE;
3182 else {fill3=kTRUE; fill4=kTRUE;}
3183 //
3184 }else if( (c1+c2+c3)==2) {
3185 b1=0; b2=1; b3=1;
3186 //
3187 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
3188 else fill4=kTRUE;
3189 //
3190 }else {
3191 b1=c1; b2=c2; b3=c3;
3192 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
3193 }
3194 }else if(fi != 5){// all the rest except pi-k-p
3195 if(key1==key2){
3196 b3=c3;
3197 if( (c1+c2)==1) {b1=0; b2=1;}
3198 else {b1=c1; b2=c2;}
3199 }else if(key1==key3){
3200 b3=c2;
3201 if( (c1+c3)==1) {b1=0; b2=1;}
3202 else {b1=c1; b2=c3;}
3203 }else {// Key2==Key3
3204 b3=c1;
3205 if( (c2+c3)==1) {b1=0; b2=1;}
3206 else {b1=c2; b2=c3;}
3207 }
3208 //////////////////////////////
3209 if(seSK) fill2=kTRUE;// Same keys from Same Event
3210 else {// Different keys from Same Event
3211 if( (c1+c2+c3)==1) {
3212 if(b3==0) {
3213 if(seSS) fill3=kTRUE;
3214 else fill4=kTRUE;
3215 }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
3216 }else if( (c1+c2+c3)==2) {
3217 if(b3==1) {
3218 if(seSS) fill4=kTRUE;
3219 else fill3=kTRUE;
3220 }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
3221 }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
3222 }
3223 /////////////////////////////
3224 }else {// pi-k-p (no charge ordering applies since all are unique)
3225 if(key1==1){
3226 if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
3227 else {b1=c1; b2=c3; b3=c2;}// pi-p-k
3228 }else if(key1==10){
3229 if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
3230 else {b1=c3; b2=c1; b3=c2;}// k-p-pi
3231 }else {// key1==100
3232 if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
3233 else {b1=c3; b2=c2; b3=c1;}// p-k-pi
3234 }
3235 ////////////////////////////////////
3236 if(seKeySum==11) fill2=kTRUE;
3237 else if(seKeySum==101) fill3=kTRUE;
3238 else fill4=kTRUE;
3239 ////////////////////////////////////
3240 }
3241
3242}
3243//________________________________________________________________________
3244void AliChaoticity::ArrangeQs(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Float_t q12, Float_t q13, Float_t q23, Short_t part, Short_t term, Float_t &fQ, Float_t &sQ, Float_t &tQ){
3245
3246 // for terms 2-4: start by setting q12(part 1) or q13(part 2)
3247 if(fi==0 || fi==6 || fi==9){// Identical species
3248 if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
3249 if(term==1 || term==5){
3250 if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
3251 else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
3252 else {fQ=q23; sQ=q12; tQ=q13;}
3253 }else if(term==2 && part==1){
3254 fQ=q12; sQ=q13; tQ=q23;
3255 }else if(term==2 && part==2){
3256 fQ=q13; sQ=q12; tQ=q23;
3257 }else if(term==3 && part==1){
3258 sQ=q12;
3259 if(c1==c3) {fQ=q13; tQ=q23;}
3260 else {fQ=q23; tQ=q13;}
3261 }else if(term==3 && part==2){
3262 sQ=q13;
3263 if(c1==c2) {fQ=q12; tQ=q23;}
3264 else {fQ=q23; tQ=q12;}
3265 }else if(term==4 && part==1){
3266 tQ=q12;
3267 if(c1==c3) {fQ=q13; sQ=q23;}
3268 else {fQ=q23; sQ=q13;}
3269 }else if(term==4 && part==2){
3270 tQ=q13;
3271 if(c1==c2) {fQ=q12; sQ=q23;}
3272 else {fQ=q23; sQ=q12;}
3273 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3274 }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
3275 if(term==1 || term==5){
3276 if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
3277 else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
3278 else {tQ=q23; sQ=q12; fQ=q13;}
3279 }else if(term==2 && part==1){
3280 fQ=q12;
3281 if(c1==c3) {tQ=q13; sQ=q23;}
3282 else {tQ=q23; sQ=q13;}
3283 }else if(term==2 && part==2){
3284 fQ=q13;
3285 if(c1==c2) {tQ=q12; sQ=q23;}
3286 else {tQ=q23; sQ=q12;}
3287 }else if(term==3 && part==1){
3288 sQ=q12;
3289 if(c1==c3) {tQ=q13; fQ=q23;}
3290 else {tQ=q23; fQ=q13;}
3291 }else if(term==3 && part==2){
3292 sQ=q13;
3293 if(c1==c2) {tQ=q12; fQ=q23;}
3294 else {tQ=q23; fQ=q12;}
3295 }else if(term==4 && part==1){
3296 tQ=q12; sQ=q13; fQ=q23;
3297 }else if(term==4 && part==2){
3298 tQ=q13; sQ=q12; fQ=q23;
3299 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3300 }else {// fQ=ss, sQ=ss, tQ=ss
3301 if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
3302 else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
3303 else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
3304 else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
3305 else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
3306 else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
3307 else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
3308 }
3309 }else if(fi != 5){// all the rest except pi-k-p
3310 if(key1==key2){
3311 fQ=q12;
3312 if(c1==c2){
3313 // cases not explicity shown below are not possible
3314 if(term==1 || term==5) {sQ=q13; tQ=q23;}
3315 else if(term==2 && part==1) {sQ=q13; tQ=q23;}
3316 else if(term==3 && part==2) {sQ=q13; tQ=q23;}
3317 else if(term==4 && part==2) {tQ=q13; sQ=q23;}
3318 else cout<<"problem!!!!!!!!!!!!!"<<endl;
3319 }else if(c3==0){
3320 if(c1==c3) {sQ=q13; tQ=q23;}
3321 else {sQ=q23; tQ=q13;}
3322 }else {//c3==1
3323 if(c1==c3) {tQ=q13; sQ=q23;}
3324 else {tQ=q23; sQ=q13;}
3325 }
3326 }else if(key1==key3){
3327 fQ=q13;
3328 if(c1==c3){
3329 // cases not explicity shown below are not possible
3330 if(term==1 || term==5) {sQ=q12; tQ=q23;}
3331 else if(term==2 && part==2) {sQ=q12; tQ=q23;}
3332 else if(term==3 && part==1) {sQ=q12; tQ=q23;}
3333 else if(term==4 && part==1) {tQ=q12; sQ=q23;}
3334 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3335 }else if(c2==0){
3336 if(c1==c2) {sQ=q12; tQ=q23;}
3337 else {sQ=q23; tQ=q12;}
3338 }else {//c2==1
3339 if(c1==c2) {tQ=q12; sQ=q23;}
3340 else {tQ=q23; sQ=q12;}
3341 }
3342 }else {// key2==key3
3343 fQ=q23;
3344 if(c2==c3){
3345 // cases not explicity shown below are not possible
3346 if(term==1 || term==5) {sQ=q12; tQ=q13;}
3347 else if(term==3 && part==1) {sQ=q12; tQ=q13;}
3348 else if(term==3 && part==2) {sQ=q13; tQ=q12;}
3349 else if(term==4 && part==1) {tQ=q12; sQ=q13;}
3350 else if(term==4 && part==2) {tQ=q13; sQ=q12;}
3351 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3352 }else if(c1==0){
3353 if(c1==c2) {sQ=q12; tQ=q13;}
3354 else {sQ=q13; tQ=q12;}
3355 }else {//c1==1
3356 if(c1==c2) {tQ=q12; sQ=q13;}
3357 else {tQ=q13; sQ=q12;}
3358 }
3359 }
3360 }else {// pi-k-p
3361 if(key1==1){
3362 if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
3363 else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
3364 }else if(key1==10){
3365 if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
3366 else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
3367 }else {// key1==100
3368 if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
3369 else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
3370 }
3371
3372 }
3373
3374
3375}
3376//________________________________________________________________________
3377Float_t AliChaoticity::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
3378
3379 Float_t qinv=1.0;
3380
3381 if(fi==0 || fi==3 || fi==5){// identical masses
3382 qinv = sqrt( pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2));
3383 }else{// different masses
3384 Float_t px = track1[1] + track2[1];
3385 Float_t py = track1[2] + track2[2];
3386 Float_t pz = track1[3] + track2[3];
3387 Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
3388 Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
3389 deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
3390
3391 qinv = pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
3392 qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
3393 qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
3394 qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
3395 qinv = sqrt(qinv);
3396 }
3397
3398 return qinv;
3399
3400}
3401//________________________________________________________________________
3402void AliChaoticity::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
3403
3404 Float_t p0 = track1[0] + track2[0];
3405 Float_t px = track1[1] + track2[1];
3406 Float_t py = track1[2] + track2[2];
3407 Float_t pz = track1[3] + track2[3];
3408
3409 Float_t mt = sqrt(p0*p0 - pz*pz);
3410 Float_t pt = sqrt(px*px + py*py);
3411
3412 Float_t v0 = track1[0] - track2[0];
3413 Float_t vx = track1[1] - track2[1];
3414 Float_t vy = track1[2] - track2[2];
3415 Float_t vz = track1[3] - track2[3];
3416
3417 qout = (px*vx + py*vy)/pt;
3418 qside = (px*vy - py*vx)/pt;
3419 qlong = (p0*vz - pz*v0)/mt;
3420}
3421//________________________________________________________________________
90814457 3422//void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F ***histos){
3423void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F *histos[3][10]){
3424
cd12341d 3425 for(Int_t mb=0; mb<fMbins; mb++){
3426 for(Int_t edB=0; edB<kEDbins; edB++){
3427 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3428 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3429 //
3430 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3431 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3432 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3433
3434 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3435 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3436 }
3437 }
3438 }
3439 }
3440 }
3441 //
3442 }
3443 }
3444
90814457 3445 if(legoCase){
3446 cout<<"LEGO call to SetWeightArrays"<<endl;
3447 for(Int_t mb=0; mb<fMbins; mb++){
3448 for(Int_t edB=0; edB<kEDbins; edB++){
3449 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3450 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3451 //
3452 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3453 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3454 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3455
3456 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinContent(qobin, qsbin, qlbin);
3457 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinError(qobin, qsbin, qlbin);
3458 }
cd12341d 3459 }
3460 }
3461 }
90814457 3462 }
3463 //
3464 }
3465 }
3466 }else{
3467
3468 TFile *wFile = new TFile("WeightFile.root","READ");
3469 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3470 else cout<<"Good Weight File Found!"<<endl;
3471
3472 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3473 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3474 for(Int_t mb=0; mb<fMbins; mb++){
3475 for(Int_t edB=0; edB<kEDbins; edB++){
3476
3477 TString *name = new TString("Weight_Kt_");
3478 *name += tKbin;
3479 name->Append("_Ky_");
3480 *name += yKbin;
3481 name->Append("_M_");
3482 *name += mb;
3483 name->Append("_ED_");
3484 *name += edB;
3485
3486 TH3F *fTempHisto = (TH3F*)wFile->Get(name->Data());
3487
3488 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3489 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3490 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3491
3492 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinContent(qobin, qsbin, qlbin);
3493 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinError(qobin, qsbin, qlbin);
3494 }
3495 }
3496 }
3497 delete fTempHisto;
3498 }//ED
3499 }//mb
3500 }//ky
3501 }//kt
3502
3503 wFile->Close();
3504 }
3505
5e3e77d6 3506 cout<<"Done reading weight file"<<endl;
cd12341d 3507
3508}
3509//________________________________________________________________________
cd12341d 3510void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3511
3512 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3513 Float_t ky=0;// central rapidity
3514 //
3515 Float_t qOut=0,qSide=0,qLong=0;
3516 GetQosl(track1, track2, qOut, qSide, qLong);
3517 qOut = fabs(qOut);
3518 qSide = fabs(qSide);
3519 qLong = fabs(qLong);
3520 //
3521
3522 if(kt < fKmeanT[0]) {fKtbinL=0; fKtbinH=0;}
3523 else if(kt >= fKmeanT[kKbinsT-1]) {fKtbinL=kKbinsT-1; fKtbinH=kKbinsT-1;}
3524 else {
3525 for(Int_t i=0; i<kKbinsT-1; i++){
3526 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtbinL=i; fKtbinH=i+1; break;}
3527 }
3528 }
3529 //
3530 if(ky < fKmeanY[0]) {fKybinL=0; fKybinH=0;}
3531 else if(ky >= fKmeanY[kKbinsY-1]) {fKybinL=kKbinsY-1; fKybinH=kKbinsY-1;}
3532 else {
3533 for(Int_t i=0; i<kKbinsY-1; i++){
3534 if((ky >= fKmeanY[i]) && (ky < fKmeanY[i+1])) {fKybinL=i; fKybinH=i+1; break;}
3535 }
3536 }
3537 //
3538 /////////
3539 if(qOut < fQmean[0]) {fQobinL=0; fQobinH=0;}
3540 else if(qOut >= fQmean[kQbinsWeights-1]) {fQobinL=kQbinsWeights-1; fQobinH=kQbinsWeights-1;}
3541 else {
3542 for(Int_t i=0; i<kQbinsWeights-1; i++){
3543 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQobinL=i; fQobinH=i+1; break;}
3544 }
3545 }
3546 //
3547 if(qSide < fQmean[0]) {fQsbinL=0; fQsbinH=0;}
3548 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsbinL=kQbinsWeights-1; fQsbinH=kQbinsWeights-1;}
3549 else {
3550 for(Int_t i=0; i<kQbinsWeights-1; i++){
3551 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsbinL=i; fQsbinH=i+1; break;}
3552 }
3553 }
3554 //
3555 if(qLong < fQmean[0]) {fQlbinL=0; fQlbinH=0;}
3556 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlbinL=kQbinsWeights-1; fQlbinH=kQbinsWeights-1;}
3557 else {
3558 for(Int_t i=0; i<kQbinsWeights-1; i++){
3559 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlbinL=i; fQlbinH=i+1; break;}
3560 }
3561 }
3562 //
3563
3564
3565 Float_t min = fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3566 Float_t minErr = fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3567
3568
3569 Float_t deltaW=0;
3570 // kt
3571 deltaW += (fNormWeight[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - min)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3572 // Ky
3573 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - min)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3574 // Qo
5e3e77d6 3575 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - min)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3576 // Qs
5e3e77d6 3577 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - min)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3578 // Ql
5e3e77d6 3579 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - min)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3580
3581 //
3582 wgt = min + deltaW;
3583
3584
3585 ////
3586
3587 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3588 Float_t deltaWErr=0;
3589 // Kt
3590 /*
3591 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - minErr)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3592 // Ky
3593 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - minErr)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3594 // Qo
5e3e77d6 3595 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - minErr)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3596 // Qs
5e3e77d6 3597 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - minErr)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3598 // Ql
5e3e77d6 3599 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - minErr)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3600 */
3601 wgtErr = minErr + deltaWErr;
3602
3603
3604
cd12341d 3605}
3606//________________________________________________________________________
3607Float_t AliChaoticity::MCWeight(Int_t charge1, Int_t charge2, Int_t rIndex, Int_t dampIndex, Float_t qinv){
3608
5e3e77d6 3609 Float_t radius = Float_t(rIndex+3.)/0.19733;// convert to GeV
cd12341d 3610 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
90814457 3611 Float_t coulCorr12 = FSICorrelationGaus2(charge1, charge2, rIndex, qinv);
cd12341d 3612 if(charge1==charge2){
5e3e77d6 3613 return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*radius,2)))*coulCorr12);
cd12341d 3614 }else {
5e3e77d6 3615 return ((1-myDamp) + myDamp*coulCorr12);
cd12341d 3616 }
3617
3618}
3619//________________________________________________________________________
90814457 3620Float_t AliChaoticity::MCWeight3D(Bool_t SameCharge, Int_t term, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23){
cd12341d 3621 if(term==5) return 1.0;
5e3e77d6 3622
90814457 3623 Float_t radius=5;
3624 if(fMbin<=1) {radius = 8;}
3625 else if(fMbin<=3) {radius = 7;}
3626 else if(fMbin<=5) {radius = 6;}
3627 else {radius = 5;}
3628 radius /= 0.19733;
3629
3630 //Float_t radius = (3. + rIndex)/0.19733;//starts at 3fm. convert to GeV
cd12341d 3631 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3632 Float_t fc = sqrt(myDamp);
5e3e77d6 3633 if(SameCharge){// all three of the same charge
90814457 3634 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2
3635 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, q13);// K2
3636 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, q23);// K2
cd12341d 3637
5e3e77d6 3638 if(term==1){
3639 Float_t c3QS = 1 + exp(-pow(q12*radius,2)) + exp(-pow(q13*radius,2)) + exp(-pow(q23*radius,2));
3640 c3QS += 2*exp(-pow(radius,2)*(pow(q12,2) + pow(q13,2) + pow(q23,2))/2.);
3641 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3642 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3643 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*radius,2)))*coulCorr13;
3644 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*radius,2)))*coulCorr23;
3645 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kTRUE, q12, q13, q23);
3646 return w123;
3647 }else if(term==2){
3648 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3649 }else if(term==3){
3650 return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*radius,2)))*coulCorr13);
3651 }else if(term==4){
3652 return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*radius,2)))*coulCorr23);
3653 }else return 1.0;
3654
3655 }else{// mixed charge case pair 12 always treated as ss
90814457 3656 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2 ss
3657 Float_t coulCorr13 = FSICorrelationTherm2(+1,-1, q13);// K2 os
3658 Float_t coulCorr23 = FSICorrelationTherm2(+1,-1, q23);// K2 os
5e3e77d6 3659 if(term==1){
3660 Float_t c3QS = 1 + exp(-pow(q12*radius,2));
3661 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3662 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3663 w123 += pow(fc,2)*(1-fc)*coulCorr13;
3664 w123 += pow(fc,2)*(1-fc)*coulCorr23;
3665 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kFALSE, q12, q13, q23);
3666 return w123;
3667 }else if(term==2){
3668 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3669 }else if(term==3){
3670 return ((1-myDamp) + myDamp*coulCorr13);
3671 }else if(term==4){
3672 return ((1-myDamp) + myDamp*coulCorr23);
3673 }else return 1.0;
3674 }
3675
cd12341d 3676}
3677//________________________________________________________________________
90814457 3678void AliChaoticity::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D){
5e3e77d6 3679
3680
3681 if(legoCase){
90814457 3682 cout<<"LEGO call to SetMomResCorrections"<<endl;
5e3e77d6 3683 fMomResC2 = (TH2D*)temp2D->Clone();
5e3e77d6 3684 fMomResC2->SetDirectory(0);
5e3e77d6 3685 }else {
3686 TFile *momResFile = new TFile("MomResFile.root","READ");
3687 if(!momResFile->IsOpen()) {
3688 cout<<"No momentum resolution file found"<<endl;
3689 AliFatal("No momentum resolution file found. Kill process.");
3690 }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
3691
3692 TH2D *temp2D2 = (TH2D*)momResFile->Get("MomResHisto_pp");
5e3e77d6 3693 fMomResC2 = (TH2D*)temp2D2->Clone();
5e3e77d6 3694 fMomResC2->SetDirectory(0);
90814457 3695
5e3e77d6 3696 momResFile->Close();
3697 }
cd12341d 3698
5e3e77d6 3699 cout<<"Done reading momentum resolution file"<<endl;
3700}
3701//________________________________________________________________________
90814457 3702void AliChaoticity::SetFSICorrelations(Bool_t legoCase, TH2D *temp2DGaus[2], TH2D *temp2DTherm[2], TH3D *temp3Dos[6], TH3D *temp3Dss[6]){
5e3e77d6 3703 // read in 2-particle and 3-particle FSI correlations = K2 & K3
3704 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
3705 if(legoCase){
c4980714 3706 cout<<"LEGO call to SetFSICorrelations"<<endl;
90814457 3707 fFSI2SS[0] = (TH2D*)temp2DGaus[0]->Clone();
3708 fFSI2OS[0] = (TH2D*)temp2DGaus[1]->Clone();
3709 fFSI2SS[1] = (TH2D*)temp2DTherm[0]->Clone();
3710 fFSI2OS[1] = (TH2D*)temp2DTherm[1]->Clone();
5e3e77d6 3711 //
90814457 3712 fFSI2SS[0]->SetDirectory(0);
3713 fFSI2OS[0]->SetDirectory(0);
3714 fFSI2SS[1]->SetDirectory(0);
3715 fFSI2OS[1]->SetDirectory(0);
3716
3717 for(Int_t CB=0; CB<6; CB++) {
3718 fFSIOmega0OS[CB] = (TH3D*)temp3Dos[CB]->Clone();
3719 fFSIOmega0SS[CB] = (TH3D*)temp3Dss[CB]->Clone();
3720 //
3721 fFSIOmega0OS[CB]->SetDirectory(0);
3722 fFSIOmega0SS[CB]->SetDirectory(0);
3723 }
5e3e77d6 3724 }else {
c4980714 3725 cout<<"non LEGO call to SetFSICorrelations"<<endl;
5e3e77d6 3726 TFile *fsifile = new TFile("KFile.root","READ");
3727 if(!fsifile->IsOpen()) {
3728 cout<<"No FSI file found"<<endl;
3729 AliFatal("No FSI file found. Kill process.");
3730 }else {cout<<"Good FSI File Found!"<<endl;}
3731
90814457 3732 TH2D *temphisto2GausSS = (TH2D*)fsifile->Get("K2ssG");
3733 TH2D *temphisto2GausOS = (TH2D*)fsifile->Get("K2osG");
3734 TH2D *temphisto2ThermSS = (TH2D*)fsifile->Get("K2ssT");
3735 TH2D *temphisto2ThermOS = (TH2D*)fsifile->Get("K2osT");
3736 TH3D *temphisto3OS[6];
ae9b34d1 3737 TH3D *temphisto3SS[6];
3738 for(Int_t CB=0; CB<6; CB++) {
90814457 3739 TString *nameK3SS = new TString("K3ss_");
3740 *nameK3SS += CB;
3741 temphisto3SS[CB] = (TH3D*)fsifile->Get(nameK3SS->Data());
3742 //
3743 TString *nameK3OS = new TString("K3os_");
3744 *nameK3OS += CB;
3745 temphisto3OS[CB] = (TH3D*)fsifile->Get(nameK3OS->Data());
ae9b34d1 3746 }
3747
90814457 3748 fFSI2SS[0] = (TH2D*)temphisto2GausSS->Clone();
3749 fFSI2OS[0] = (TH2D*)temphisto2GausOS->Clone();
3750 fFSI2SS[1] = (TH2D*)temphisto2ThermSS->Clone();
3751 fFSI2OS[1] = (TH2D*)temphisto2ThermOS->Clone();
3752 fFSI2SS[0]->SetDirectory(0);
3753 fFSI2OS[0]->SetDirectory(0);
3754 fFSI2SS[1]->SetDirectory(0);
3755 fFSI2OS[1]->SetDirectory(0);
ae9b34d1 3756
90814457 3757 for(Int_t CB=0; CB<6; CB++) {
3758 fFSIOmega0SS[CB] = (TH3D*)temphisto3SS[CB]->Clone();
3759 fFSIOmega0OS[CB] = (TH3D*)temphisto3OS[CB]->Clone();
3760 fFSIOmega0SS[CB]->SetDirectory(0);
3761 fFSIOmega0OS[CB]->SetDirectory(0);
3762 }
3763 //
3764
5e3e77d6 3765 fsifile->Close();
cd12341d 3766 }
3767
5e3e77d6 3768 // condition FSI histogram for edge effects
ae9b34d1 3769 for(Int_t CB=0; CB<6; CB++){
3770 for(Int_t ii=1; ii<=fFSIOmega0SS[CB]->GetNbinsX(); ii++){
3771 for(Int_t jj=1; jj<=fFSIOmega0SS[CB]->GetNbinsY(); jj++){
3772 for(Int_t kk=1; kk<=fFSIOmega0SS[CB]->GetNbinsZ(); kk++){
5e3e77d6 3773
ae9b34d1 3774 if(fFSIOmega0SS[CB]->GetBinContent(ii,jj,kk) <=0){
3775 Double_t Q12 = fFSIOmega0SS[CB]->GetXaxis()->GetBinCenter(ii);
3776 Double_t Q23 = fFSIOmega0SS[CB]->GetYaxis()->GetBinCenter(jj);
3777 Double_t Q13 = fFSIOmega0SS[CB]->GetZaxis()->GetBinCenter(kk);
3778 //
3779 Int_t Q12bin=ii;
3780 Int_t Q23bin=jj;
3781 Int_t Q13bin=kk;
3782 Int_t AC=0;//Adjust Counter
3783 Int_t AClimit=10;// maximum bin shift
3784 if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin++; AC++;}}
3785 if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin--; AC++;}}
3786 //
3787 if(Q13 < sqrt(pow(Q12,2)+pow(Q23,2) - 2*Q12*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin++; AC++;}}
3788 if(Q13 > sqrt(pow(Q12,2)+pow(Q23,2) + 2*Q12*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin--; AC++;}}
3789 //
3790 if(Q23 < sqrt(pow(Q12,2)+pow(Q13,2) - 2*Q12*Q13)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin++; AC++;}}
3791 if(Q23 > sqrt(pow(Q12,2)+pow(Q13,2) + 2*Q12*Q13)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin--; AC++;}}
3792
3793 // Save cpu time by setting empty cell contents (edge effects) to nearest non-zero cell (these cells are not used very often anyway.)
3794 if(AC==AClimit) {
3795 fFSIOmega0SS[CB]->SetBinContent(ii,jj,kk, 1.0);
90814457 3796 fFSIOmega0OS[CB]->SetBinContent(ii,jj,kk, 1.0);
ae9b34d1 3797 }else {
3798 fFSIOmega0SS[CB]->SetBinContent(ii,jj,kk, fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin));
90814457 3799 fFSIOmega0OS[CB]->SetBinContent(ii,jj,kk, fFSIOmega0OS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin));
ae9b34d1 3800 }
5e3e77d6 3801 }
ae9b34d1 3802
5e3e77d6 3803 }
5e3e77d6 3804 }
cd12341d 3805 }
3806 }
5e3e77d6 3807 cout<<"Done reading FSI file"<<endl;
cd12341d 3808}
3809//________________________________________________________________________
90814457 3810Float_t AliChaoticity::FSICorrelationGaus2(Int_t charge1, Int_t charge2, Int_t rIndex, Float_t qinv){
5e3e77d6 3811 // returns 2-particle Coulomb correlations = K2
3812 if(rIndex >= kRVALUES) return 1.0;
90814457 3813 Int_t qbinL = fFSI2SS[0]->GetYaxis()->FindBin(qinv-fFSI2SS[0]->GetYaxis()->GetBinWidth(1)/2.);
3814 Int_t qbinH = qbinL+1;
3815 if(qbinL <= 0) return 1.0;
3816 if(qbinH > fFSI2SS[0]->GetNbinsY()) return 1.0;
3817
3818 Float_t slope=0;
3819 if(charge1==charge2){
3820 slope = fFSI2SS[0]->GetBinContent(rIndex+1, qbinL) - fFSI2SS[0]->GetBinContent(rIndex+1, qbinH);
3821 slope /= fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinH);
3822 return (slope*(qinv - fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS[0]->GetBinContent(rIndex+1, qbinL));
3823 }else {
3824 slope = fFSI2OS[0]->GetBinContent(rIndex+1, qbinL) - fFSI2OS[0]->GetBinContent(rIndex+1, qbinH);
3825 slope /= fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinH);
3826 return (slope*(qinv - fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS[0]->GetBinContent(rIndex+1, qbinL));
3827 }
3828}
3829//________________________________________________________________________
3830Float_t AliChaoticity::FSICorrelationTherm2(Int_t charge1, Int_t charge2, Float_t qinv){
3831 // returns 2-particle Coulomb correlations = K2
3832 Int_t qbinL = fFSI2SS[1]->GetYaxis()->FindBin(qinv-fFSI2SS[1]->GetYaxis()->GetBinWidth(1)/2.);
5e3e77d6 3833 Int_t qbinH = qbinL+1;
3834 if(qbinL <= 0) return 1.0;
90814457 3835 if(qbinH > fFSI2SS[1]->GetNbinsY()) return 1.0;
5e3e77d6 3836
3837 Float_t slope=0;
3838 if(charge1==charge2){
90814457 3839 slope = fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinL) - fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinH);
3840 slope /= fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinH);
3841 return (slope*(qinv - fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinL));
5e3e77d6 3842 }else {
90814457 3843 slope = fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinL) - fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinH);
3844 slope /= fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinH);
3845 return (slope*(qinv - fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinL));
cd12341d 3846 }
3847}
3848//________________________________________________________________________
5e3e77d6 3849Double_t AliChaoticity::FSICorrelationOmega0(Bool_t SameCharge, Double_t Q12, Double_t Q13, Double_t Q23){
3850 // remember: Omega0 histogram is in the following order (Q12, Q23, Q13)
3851 // returns 3d 3-particle Coulomb Correlation = K3
ae9b34d1 3852 Int_t Q12bin = fFSIOmega0SS[fFSIbin]->GetXaxis()->FindBin(Q12);
3853 Int_t Q13bin = fFSIOmega0SS[fFSIbin]->GetZaxis()->FindBin(Q13);
3854 Int_t Q23bin = fFSIOmega0SS[fFSIbin]->GetYaxis()->FindBin(Q23);
5e3e77d6 3855
3856 if(SameCharge){
ae9b34d1 3857 if(fFSIOmega0SS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3858 else return fFSIOmega0SS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
90814457 3859 }else{// mixed charge.
3860 if(fFSIOmega0OS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3861 else return fFSIOmega0OS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
5e3e77d6 3862 }
cd12341d 3863}
90814457 3864//________________________________________________________________________
3865void AliChaoticity::FourVectProdTerms(Float_t pV1[], Float_t pV2[], Float_t pV3[], Float_t& QS1v1, Float_t& QS2, Float_t& QS3v1, Float_t& QS1v2, Float_t& QS3v2){
3866 QS1v1 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
3867 QS1v1 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
3868 QS1v1 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
3869 QS2 = (pV1[1]-pV2[1])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[1]-pV3[1]);
3870 QS3v1 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
3871 QS3v1 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
3872 //
3873 QS1v2 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
3874 QS1v2 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
3875 QS3v2 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
3876 QS3v2 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
3877 QS3v2 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
3878}
3879//________________________________________________________________________
3880void AliChaoticity::TestAddTask(){
3881 /*
3882 TString inputFileNameFSI = "KFile.root";
3883 TFile *inputFileFSI = TFile::Open(inputFileNameFSI,"OLD");
3884 if (!inputFileFSI){
3885 cout << "Requested file:" << inputFileFSI << " was not opened. ABORT." << endl;
3886 return;
3887 }
3888 TH2D *FSI2gaus[2];
3889 TH2D *FSI2therm[2];
3890 TH3D *FSI3ss[6];
3891 TH3D *FSI3os[6];
3892 FSI2gaus[0] = (TH2D*)inputFileFSI->Get("K2ssG");
3893 FSI2gaus[1] = (TH2D*)inputFileFSI->Get("K2osG");
3894 FSI2therm[0] = (TH2D*)inputFileFSI->Get("K2ssT");
3895 FSI2therm[1] = (TH2D*)inputFileFSI->Get("K2osT");
3896 for(Int_t CB=0; CB<6; CB++) {
3897 TString *nameSS=new TString("K3ss_");
3898 *nameSS += CB;
3899 FSI3ss[CB] = (TH3D*)inputFileFSI->Get(nameSS->Data());
3900 TString *nameOS=new TString("K3os_");
3901 *nameOS += CB;
3902 FSI3os[CB] = (TH3D*)inputFileFSI->Get(nameOS->Data());
3903 }
3904 //
3905 FSI2gaus[0]->SetDirectory(0);
3906 FSI2gaus[1]->SetDirectory(0);
3907 FSI2therm[0]->SetDirectory(0);
3908 FSI2therm[1]->SetDirectory(0);
3909 for(Int_t CB=0; CB<6; CB++) {
3910 FSI3ss[CB]->SetDirectory(0);
3911 FSI3os[CB]->SetDirectory(0);
3912 }
3913 SetFSICorrelations( kTRUE, FSI2gaus, FSI2therm , FSI3os, FSI3ss);
3914 //
3915
3916 if(!fTabulatePairs) {
3917 TString inputFileNameWeight = "WeightFile.root";
3918 TFile *inputFileWeight = TFile::Open(inputFileNameWeight,"OLD");
3919 if (!inputFileWeight){
3920 cout << "Requested file:" << inputFileWeight << " was not opened. ABORT." << endl;
3921 return;
3922 }
3923
3924 ////////////////////////////////////////////////////
3925 // C2 Weight File
3926 Int_t ktbins = GetNumKtbins();
3927 Int_t cbins = GetNumCentBins();
3928 TH3F *weightHisto[ktbins][cbins];
3929 for(Int_t i=0; i<ktbins; i++){
3930 for(Int_t j=0; j<cbins; j++){
3931 TString name = "Weight_Kt_";
3932 name += i;
3933 name += "_Ky_0_M_";
3934 name += j;
3935 name += "_ED_0";
3936
3937 weightHisto[i][j] = (TH3F*)inputFileWeight->Get(name);
3938 }
3939 }
3940 SetWeightArrays( kTRUE, weightHisto );
3941 }
3942
3943 //
3944 if(!fMCcase && !fTabulatePairs){
3945 TString inputFileNameMomRes = "MomResFile.root";
3946 TFile *inputFileMomRes = TFile::Open(inputFileNameMomRes,"OLD");
3947 if (!inputFileMomRes){
3948 cout << "Requested file:" << inputFileMomRes << " was not opened. ABORT." << endl;
3949 return;
3950 }
3951 ////////////////////////////////////////////////////
3952 // Momentum Resolution File
3953 TH2D *momResHisto2D = 0;
3954 momResHisto2D = (TH2D*)inputFileMomRes->Get("MomResHisto_pp");
3955 SetMomResCorrections( kTRUE, momResHisto2D);
3956 }
3957 */
3958}