]> 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
fix because trunk does not compile
[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
034e970d 2592 if(bin1==bin2) {
2593 WInput = MCWeight3D(kFALSE, 1, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2594 K3 = FSICorrelationOmega0(kFALSE, firstQMC, secondQMC, thirdQMC);// K3
2595 }else {
2596 WInput = MCWeight3D(kFALSE, 1, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2597 K3 = FSICorrelationOmega0(kFALSE, thirdQMC, secondQMC, firstQMC);// K3
2598 }
5e3e77d6 2599 }
54d66278 2600 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2601 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2602 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fQW12->Fill(firstQMC, secondQMC, thirdQMC, WInput*firstQMC);
2603 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fQW13->Fill(firstQMC, secondQMC, thirdQMC, WInput*secondQMC);
90814457 2604 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSumK3->Fill(firstQMC, secondQMC, thirdQMC, WInput/K3);
2605 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fEnK3->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2606 if(ch1==ch2 && ch1==ch3){
2607 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2608 FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2609 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2610 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2611 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2612 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2613 //
2614 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsSumK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K3);
2615 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsSumK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K3);
2616 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsEnK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2617 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsEnK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2618 }
2619
fa109294 2620 }// fMCcase
2621
cd12341d 2622 }
2623
2624 }else if(config==2){// 12, 13, 23
5e3e77d6 2625
cd12341d 2626 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2627 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
2628
2629 // loop over terms 2-4
5e3e77d6 2630 for(Int_t jj=2; jj<5; jj++){
2631 if(jj==2) {if(!fill2) continue;}//12
2632 if(jj==3) {if(!fill3) continue;}//13
2633 if(jj==4) {if(!fill4) continue;}//23
cd12341d 2634
2635 if(fillIndex3 <= 2){
5e3e77d6 2636 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj, firstQ, secondQ, thirdQ);
2637 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTerms3->Fill(firstQ, secondQ, thirdQ);
90814457 2638 if(fillIndex3==0 && ch1==ch2 && ch1==ch3){
2639 if(part==1){// P11T2
2640 if(jj==2) {
2641 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2642 if(fMCcase) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2643 }if(jj==3){
2644 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2645 if(fMCcase) FourVectProdTerms(pVect1MC, pVect3MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2646 }if(jj==4) {
2647 FourVectProdTerms(pVect3, pVect1, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2648 if(fMCcase) FourVectProdTerms(pVect3MC, pVect1MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2649 }
2650 }else{// P12T1
2651 if(jj==2) {
2652 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2653 if(fMCcase) FourVectProdTerms(pVect1MC, pVect3MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2654 }if(jj==3) {
2655 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2656 if(fMCcase) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2657 }if(jj==4) {
2658 FourVectProdTerms(pVect2, pVect1, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2659 if(fMCcase) FourVectProdTerms(pVect2MC, pVect1MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2660 }
2661 }
2662 if(!fMCcase){
2663 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2664 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
fa109294 2665 }
5e3e77d6 2666 }
2667 //////////////////////////////////////
2668 // Momentum resolution calculation
2669 if(fillIndex3==0 && fMCcase){
2670 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2671 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, jj, firstQMC, secondQMC, thirdQMC);
90814457 2672 Float_t WInput = 1;
5e3e77d6 2673 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2674 WInput = MCWeight3D(kTRUE, jj, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
5e3e77d6 2675 }else {// mixed charge
90814457 2676 if(bin1==bin2) WInput = MCWeight3D(kFALSE, jj, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2677 else WInput = MCWeight3D(kFALSE, 6-jj, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2678 }
2679 //
2680 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2681 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2682 //
2683 if(ch1==ch2 && ch1==ch3){
2684 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2685 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2686 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2687 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2688 //
2689 Float_t InteractingQ=0;
2690 if(part==1) {InteractingQ=qinv12;}// 12 from SE
2691 else {InteractingQ=qinv13;}// 13 from SE
2692 Double_t K2 = FSICorrelationTherm2(+1,+1, InteractingQ);// K2 from Therminator source
2693 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSumK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K2);
2694 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSumK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K2);
2695 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsEnK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2696 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsEnK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
5e3e77d6 2697 }
2698 }// fMCcase
2699
cd12341d 2700 }
2701 }
5e3e77d6 2702
cd12341d 2703 }else {// config 3: All particles from different events
2704 //Simulate Filling of other event-mixed lowQ pairs
2705
2706 Float_t enhancement=1.0;
2707 Int_t nUnderCut=0;
2708 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2709 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2710
2711 if(nUnderCut==0) enhancement = (1+1+1)/1.;// 1 LowQ pair
2712 if(nUnderCut==1) enhancement = (1+2)/2.;// 2 LowQ pair
2713 if(nUnderCut==2) enhancement = 1.;// 3 LowQ pair
2714
2715 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
90814457 2716
2717 if(ch1==ch2 && ch1==ch3 && fillIndex3==0) FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2718 if(fMCcase && ch1==ch2 && ch1==ch3 && fillIndex3==0) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);
cd12341d 2719
2720 if(fillIndex3 <= 2){
2721 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
2722 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ, enhancement);
5e3e77d6 2723 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
90814457 2724 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2725 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
5e3e77d6 2726 }
2727 //////////////////////////////////////
2728 // Momentum resolution calculation
2729 if(fillIndex3==0 && fMCcase){
2730 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2731 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, 5, firstQMC, secondQMC, thirdQMC);
90814457 2732 Float_t WInput=1;
5e3e77d6 2733 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2734 WInput = MCWeight3D(kTRUE, 5, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
5e3e77d6 2735 }else {// mixed charge
90814457 2736 if(bin1==bin2) WInput = MCWeight3D(kFALSE, 5, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2737 else WInput = MCWeight3D(kFALSE, 5, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss in this case. 1st Q argument is ss
2738 }
2739 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2740 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2741 if(ch1==ch2 && ch1==ch3){
2742 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2743 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2744 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2745 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2746
5e3e77d6 2747 }
2748 }// fMCcase
2749
cd12341d 2750 }
5e3e77d6 2751
cd12341d 2752 if(fillIndex3 !=0) continue;// only calculate TPN for pi-pi-pi
2753 if(ch1!=ch2 || ch1!=ch3) continue;// only calcualte TPN for ss
5e3e77d6 2754
90814457 2755 //if(fMCcase) continue;// only calcualte TPN for real data
2756
cd12341d 2757 GetWeight(pVect1, pVect2, weight12, weight12Err);
2758 GetWeight(pVect1, pVect3, weight13, weight13Err);
2759 GetWeight(pVect2, pVect3, weight23, weight23Err);
2760 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) continue;// weight should never be larger than 1
2761
5e3e77d6 2762
2763 // Coul correlations from Wave-functions
fa109294 2764 //for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 8fm, last value for Therminator
2765 //for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){// 0.3 to 0.6
2766 //Float_t myDamp = fDampStart + (fDampStep)*myDampIt;
2767 //Int_t denIndex = (kRVALUES-1)*kNDampValues + myDampIt;
2768 //Int_t denIndex = myDampIt;
2769 Int_t myDampIt = 5;
2770 Float_t myDamp = 0.4;
2771 Int_t denIndex = 0;
90814457 2772 Int_t momResIndex = rIndexForTPN*kNDampValues + myDampIt;// Therminator slot uses R=7 for momentum resolution
2773
2774 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, qinv12);
2775 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, qinv13);
2776 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, qinv23);
cd12341d 2777
90814457 2778 Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
2779 if(!fMCcase) {
2780 Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
2781 Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
2782 Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);
2783 if(momBin12 >= kQbins) momBin12 = kQbins-1;
2784 if(momBin13 >= kQbins) momBin13 = kQbins-1;
2785 if(momBin23 >= kQbins) momBin23 = kQbins-1;
2786 MomResCorr12 = fMomResC2->GetBinContent(momResIndex+1, momBin12);
2787 MomResCorr13 = fMomResC2->GetBinContent(momResIndex+1, momBin13);
2788 MomResCorr23 = fMomResC2->GetBinContent(momResIndex+1, momBin23);
2789 }
2790 weight12CC = ((weight12+1)*MomResCorr12 - myDamp*coulCorr12 - (1-myDamp));
5e3e77d6 2791 weight12CC /= coulCorr12*myDamp;
90814457 2792 weight13CC = ((weight13+1)*MomResCorr13 - myDamp*coulCorr13 - (1-myDamp));
5e3e77d6 2793 weight13CC /= coulCorr13*myDamp;
90814457 2794 weight23CC = ((weight23+1)*MomResCorr23 - myDamp*coulCorr23 - (1-myDamp));
5e3e77d6 2795 weight23CC /= coulCorr23*myDamp;
cd12341d 2796
5e3e77d6 2797 if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {// testing purposes only
cd12341d 2798 if(denIndex==71 && fMbin==0 && ch1==-1) {
2799 weightTotal = sqrt(fabs(weight12CC*weight13CC*weight23CC));
2800 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
2801 }
2802 continue;// C2^QS can never be less than unity
2803 }
2804
2805 weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
90814457 2806
cd12341d 2807 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].fTwoPartNorm->Fill(qinv12, qinv13, qinv23, enhancement*weightTotal);
90814457 2808
2809 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNorm->Fill(Qsum1v1, Qsum2, Qsum3v1, enhancement*weightTotal);
2810 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNorm->Fill(Qsum1v2, Qsum2, Qsum3v2, enhancement*weightTotal);
2811 if(fMCcase){
2812 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNormIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, enhancement*weightTotal);
2813 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNormSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, enhancement*weightTotal);
2814 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNormIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, enhancement*weightTotal);
2815 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNormSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, enhancement*weightTotal);
fa109294 2816 }
90814457 2817
2818
cd12341d 2819 // Save cpu time and memory by skipping r3 denominator calculation below. den errors are negligible compared to num errors.
2820 /*
2821 if(weightTotal > 0.0001){// tiny numbers cause a Float_ting point exception below
2822 weightTotalErr = pow((weight12Err*coulCorr12)*weight13CC*weight23CC,2);
2823 weightTotalErr += pow(weight12CC*(weight13Err*coulCorr13)*weight23CC,2);
2824 weightTotalErr += pow(weight12CC*weight13CC*(weight23Err*coulCorr23),2);
2825 weightTotalErr /= pow(2*weightTotal,2);
2826
2827 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].TwoPartNormErr->Fill(denIndex, q3, enhancement*weightTotalErr);
2828 }
2829 */
2830
fa109294 2831 //}// damp iter
2832 //}// R iter
cd12341d 2833
2834
2835
2836 }
2837 }// end 3rd particle
2838 }// en2
2839
2840
fa109294 2841 }// p1
2842 }//en1
cd12341d 2843
cd12341d 2844 ///////////////////
2845 }// end of PdensityPairs
2846
fa109294 2847
cd12341d 2848
2849
2850
2851
2852
2853 ////////////////////////////////////////////////////////
2854 // Pdensity Method with Explicit Loops
2855 if(fPdensityExplicitLoop){
2856
2857 ////////////////////////////////////
2858 // 2nd, 3rd, and 4th order Correlations
2859
2860 // First Particle
2861 for (Int_t i=0; i<myTracks; i++) {
2862 ch1 = Int_t( ((fEvt)->fTracks[i].fCharge + 1)/2. );
2863 pVect1[0] = (fEvt)->fTracks[i].fEaccepted;
2864 pVect1[1] = (fEvt)->fTracks[i].fP[0];
2865 pVect1[2] = (fEvt)->fTracks[i].fP[1];
2866 pVect1[3] = (fEvt)->fTracks[i].fP[2];
2867 key1 = (fEvt)->fTracks[i].fKey;
2868
2869 // Second Event
2870 for(Int_t en2=0; en2<fEventsToMix+1; en2++){
2871 Int_t startbin2=0;
2872 if(en2==0) startbin2=i+1;
2873
2874 // Second Particle
2875 for (Int_t j=startbin2; j<(fEvt+en2)->fNtracks; j++) {
2876 ch2 = Int_t( ((fEvt+en2)->fTracks[j].fCharge + 1)/2. );
2877 pVect2[0] = (fEvt+en2)->fTracks[j].fEaccepted;
2878 pVect2[1] = (fEvt+en2)->fTracks[j].fP[0];
2879 pVect2[2] = (fEvt+en2)->fTracks[j].fP[1];
2880 pVect2[3] = (fEvt+en2)->fTracks[j].fP[2];
2881 key2 = (fEvt+en2)->fTracks[j].fKey;
2882
2883 Short_t fillIndex12 = FillIndex2part(key1+key2);
2884 qinv12 = GetQinv(fillIndex12, pVect1, pVect2);
2885
2886 if(qinv12 < fQLowerCut) continue;
2887
2888
2889 // 2-particle part is filled always during pair creator
2890
2891 // Third Event
2892 for(Int_t en3=en2; en3<fEventsToMix+1; en3++){
2893 Int_t startbin3=0;
2894 if(en3==en2) startbin3=j+1;
2895 else startbin3=0;
2896
2897
2898 // Third Particle
2899 for (Int_t k=startbin3; k<(fEvt+en3)->fNtracks; k++) {
2900 ch3 = Int_t( ((fEvt+en3)->fTracks[k].fCharge + 1)/2. );
2901 pVect3[0] = (fEvt+en3)->fTracks[k].fEaccepted;
2902 pVect3[1] = (fEvt+en3)->fTracks[k].fP[0];
2903 pVect3[2] = (fEvt+en3)->fTracks[k].fP[1];
2904 pVect3[3] = (fEvt+en3)->fTracks[k].fP[2];
2905 key3 = (fEvt+en3)->fTracks[k].fKey;
2906
2907 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2908 Short_t fillIndex13 = FillIndex2part(key1+key3);
2909 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2910 Short_t fillIndex23 = FillIndex2part(key2+key3);
2911 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2912
2913
2914 if(qinv13 < fQLowerCut) continue;
2915 if(qinv23 < fQLowerCut) continue;
2916
2917
2918 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2919
2920 Short_t normBin12 = SetNormBin(fillIndex12);
2921 Short_t normBin13 = SetNormBin(fillIndex13);
2922 Short_t normBin23 = SetNormBin(fillIndex23);
2923
2924
2925 if(en3==0 && en2==0) {// 123
2926 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2927
2928 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fExplicit3->Fill(q3);// 123
2929
2930 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2931 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2932 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fNormEx3->Fill(0.);
2933 }
2934 }
2935
2936 }else if((en2==0 && en3==1) ) {// 12-3, 13-2, 23-1
2937 Float_t gFact=1;
2938
2939 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2940 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2941
2942
2943 if(fill2){
2944 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fExplicit3->Fill(q3, gFact);// 12
2945 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2946 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2947 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fNormEx3->Fill(0.);
2948 }
2949 }
2950 }
2951 if(fill3){
2952 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fExplicit3->Fill(q3, gFact);// 12
2953 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2954 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2955 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fNormEx3->Fill(0.);
2956 }
2957 }
2958 }
2959 if(fill4){
2960 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fExplicit3->Fill(q3, gFact);// 12
2961 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2962 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2963 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fNormEx3->Fill(0.);
2964 }
2965 }
2966 }
2967
2968 }else if(en2==1 && en3==2){// all uncorrelated events
2969 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2970
2971 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fExplicit3->Fill(q3);
2972 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2973 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2974 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fNormEx3->Fill(0.);
2975 }
2976 }
2977 Short_t qCutBin12 = SetQcutBin(fillIndex12);
2978 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2979 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2980
2981 if( (qinv12 < fQcut[qCutBin12]) || (qinv13 < fQcut[qCutBin13]) || (qinv23 < fQcut[qCutBin23])){
2982
2983 Int_t nUnderCut=0;
2984 if(qinv12<fQcut[qCutBin12]) nUnderCut++;
2985 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2986 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2987
2988 }
2989
2990 }else {}
2991
2992
2993 }// 3rd particle
2994 }// 3rd event
2995
2996 }// 2nd particle
2997 }// 2nd event
2998
2999 }// 1st particle
3000
3001
3002
3003
3004 }// End of PdensityExplicit
3005
3006
3007
3008
3009 // Post output data.
3010 PostData(1, fOutputList);
3011
3012}
3013//________________________________________________________________________
3014void AliChaoticity::Terminate(Option_t *)
3015{
3016 // Called once at the end of the query
3017
3018 cout<<"Done"<<endl;
3019
3020}
3021//________________________________________________________________________
3022Bool_t AliChaoticity::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
3023{
3024
3025 if(fabs(first.fEta-second.fEta) > fMinSepTPCEntranceEta) return kTRUE;
3026
3027 // propagate through B field to r=1m
5e3e77d6 3028 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
cd12341d 3029 if(phi1 > 2*PI) phi1 -= 2*PI;
3030 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 3031 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
cd12341d 3032 if(phi2 > 2*PI) phi2 -= 2*PI;
3033 if(phi2 < 0) phi2 += 2*PI;
3034
3035 Float_t deltaphi = phi1 - phi2;
3036 if(deltaphi > PI) deltaphi -= 2*PI;
3037 if(deltaphi < -PI) deltaphi += 2*PI;
3038 deltaphi = fabs(deltaphi);
3039
5e3e77d6 3040 //cout<<deltaphi<<" "<<fMinSepTPCEntrancePhi<<" "<<fMinSepTPCEntranceEta<<endl;
3041 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 3042
3043 // propagate through B field to r=1.6m
5e3e77d6 3044 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
cd12341d 3045 if(phi1 > 2*PI) phi1 -= 2*PI;
3046 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 3047 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
cd12341d 3048 if(phi2 > 2*PI) phi2 -= 2*PI;
3049 if(phi2 < 0) phi2 += 2*PI;
3050
3051 deltaphi = phi1 - phi2;
3052 if(deltaphi > PI) deltaphi -= 2*PI;
3053 if(deltaphi < -PI) deltaphi += 2*PI;
3054 deltaphi = fabs(deltaphi);
3055
5e3e77d6 3056 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 3057
3058
3059 //
3060
3061 Int_t ncl1 = first.fClusterMap.GetNbits();
3062 Int_t ncl2 = second.fClusterMap.GetNbits();
3063 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
3064 Double_t shfrac = 0; Double_t qfactor = 0;
3065 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
3066 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
3067 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
3068 sumQ++;
3069 sumCls+=2;
3070 sumSha+=2;}
3071 else {sumQ--; sumCls+=2;}
3072 }
3073 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
3074 sumQ++;
3075 sumCls++;}
3076 }
3077 if (sumCls>0) {
3078 qfactor = sumQ*1.0/sumCls;
3079 shfrac = sumSha*1.0/sumCls;
3080 }
3081
3082 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
3083
3084
3085 return kTRUE;
5e3e77d6 3086
cd12341d 3087
3088}
3089//________________________________________________________________________
3090Float_t AliChaoticity::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
3091{
3092 Float_t arg = G_Coeff/qinv;
3093
3094 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
3095 else {return (exp(-arg)-1)/(-arg);}
3096
3097}
3098//________________________________________________________________________
3099void AliChaoticity::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
3100{
3101 Int_t j, k;
3102 Int_t a = i2 - i1;
3103 for (Int_t i = i1; i < i2+1; i++) {
3104 j = (Int_t) (gRandom->Rndm() * a);
3105 k = iarr[j];
3106 iarr[j] = iarr[i];
3107 iarr[i] = k;
3108 }
3109}
3110//________________________________________________________________________
3111Short_t AliChaoticity::FillIndex2part(Short_t key){
3112
3113 if(key==2) return 0;// pi-pi
3114 else if(key==11) return 1;// pi-k
3115 else if(key==101) return 2;// pi-p
3116 else if(key==20) return 3;// k-k
3117 else if(key==110) return 4;// k-p
3118 else return 5;// p-p
3119}
3120//________________________________________________________________________
3121Short_t AliChaoticity::FillIndex3part(Short_t key){
3122
3123 if(key==3) return 0;// pi-pi-pi
3124 else if(key==12) return 1;// pi-pi-k
3125 else if(key==21) return 2;// k-k-pi
3126 else if(key==102) return 3;// pi-pi-p
3127 else if(key==201) return 4;// p-p-pi
3128 else if(key==111) return 5;// pi-k-p
3129 else if(key==30) return 6;// k-k-k
3130 else if(key==120) return 7;// k-k-p
3131 else if(key==210) return 8;// p-p-k
3132 else return 9;// p-p-p
3133
3134}
3135//________________________________________________________________________
3136Short_t AliChaoticity::SetQcutBin(Short_t fi){// fi=FillIndex
3137 if(fi <= 2) return 0;
3138 else if(fi==3) return 1;
3139 else return 2;
3140}
3141//________________________________________________________________________
3142Short_t AliChaoticity::SetNormBin(Short_t fi){// fi=FillIndex
3143 if(fi==0) return 0;
3144 else if(fi <= 2) return 1;
3145 else return 2;
3146}
3147//________________________________________________________________________
3148void AliChaoticity::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
3149
3150 if(fi==0 || fi==3 || fi==5){// Identical species
3151 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
3152 else {b1=c1; b2=c2;}
3153 }else {// Mixed species
3154 if(key1 < key2) { b1=c1; b2=c2;}
3155 else {b1=c2; b2=c1;}
3156 }
3157
3158}
3159//________________________________________________________________________
3160void 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){
3161
3162
3163 // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
3164 Bool_t seSS=kFALSE;
3165 Bool_t seSK=kFALSE;
3166 Short_t seKeySum=0;// only used for pi-k-p case
3167 if(part==1) {// default case (irrelevant for term 1 and term 5)
3168 if(c1==c2) seSS=kTRUE;
3169 if(key1==key2) seSK=kTRUE;
3170 seKeySum = key1+key2;
3171 }
3172 if(part==2){
3173 if(c1==c3) seSS=kTRUE;
3174 if(key1==key3) seSK=kTRUE;
3175 seKeySum = key1+key3;
3176 }
3177
3178
3179 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3180
3181 if(fi==0 || fi==6 || fi==9){// Identical species
3182 if( (c1+c2+c3)==1) {
3183 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
3184 //
3185 if(seSS) fill2=kTRUE;
3186 else {fill3=kTRUE; fill4=kTRUE;}
3187 //
3188 }else if( (c1+c2+c3)==2) {
3189 b1=0; b2=1; b3=1;
3190 //
3191 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
3192 else fill4=kTRUE;
3193 //
3194 }else {
3195 b1=c1; b2=c2; b3=c3;
3196 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
3197 }
3198 }else if(fi != 5){// all the rest except pi-k-p
3199 if(key1==key2){
3200 b3=c3;
3201 if( (c1+c2)==1) {b1=0; b2=1;}
3202 else {b1=c1; b2=c2;}
3203 }else if(key1==key3){
3204 b3=c2;
3205 if( (c1+c3)==1) {b1=0; b2=1;}
3206 else {b1=c1; b2=c3;}
3207 }else {// Key2==Key3
3208 b3=c1;
3209 if( (c2+c3)==1) {b1=0; b2=1;}
3210 else {b1=c2; b2=c3;}
3211 }
3212 //////////////////////////////
3213 if(seSK) fill2=kTRUE;// Same keys from Same Event
3214 else {// Different keys from Same Event
3215 if( (c1+c2+c3)==1) {
3216 if(b3==0) {
3217 if(seSS) fill3=kTRUE;
3218 else fill4=kTRUE;
3219 }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
3220 }else if( (c1+c2+c3)==2) {
3221 if(b3==1) {
3222 if(seSS) fill4=kTRUE;
3223 else fill3=kTRUE;
3224 }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
3225 }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
3226 }
3227 /////////////////////////////
3228 }else {// pi-k-p (no charge ordering applies since all are unique)
3229 if(key1==1){
3230 if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
3231 else {b1=c1; b2=c3; b3=c2;}// pi-p-k
3232 }else if(key1==10){
3233 if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
3234 else {b1=c3; b2=c1; b3=c2;}// k-p-pi
3235 }else {// key1==100
3236 if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
3237 else {b1=c3; b2=c2; b3=c1;}// p-k-pi
3238 }
3239 ////////////////////////////////////
3240 if(seKeySum==11) fill2=kTRUE;
3241 else if(seKeySum==101) fill3=kTRUE;
3242 else fill4=kTRUE;
3243 ////////////////////////////////////
3244 }
3245
3246}
3247//________________________________________________________________________
3248void 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){
3249
3250 // for terms 2-4: start by setting q12(part 1) or q13(part 2)
3251 if(fi==0 || fi==6 || fi==9){// Identical species
3252 if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
3253 if(term==1 || term==5){
3254 if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
3255 else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
3256 else {fQ=q23; sQ=q12; tQ=q13;}
3257 }else if(term==2 && part==1){
3258 fQ=q12; sQ=q13; tQ=q23;
3259 }else if(term==2 && part==2){
3260 fQ=q13; sQ=q12; tQ=q23;
3261 }else if(term==3 && part==1){
3262 sQ=q12;
3263 if(c1==c3) {fQ=q13; tQ=q23;}
3264 else {fQ=q23; tQ=q13;}
3265 }else if(term==3 && part==2){
3266 sQ=q13;
3267 if(c1==c2) {fQ=q12; tQ=q23;}
3268 else {fQ=q23; tQ=q12;}
3269 }else if(term==4 && part==1){
3270 tQ=q12;
3271 if(c1==c3) {fQ=q13; sQ=q23;}
3272 else {fQ=q23; sQ=q13;}
3273 }else if(term==4 && part==2){
3274 tQ=q13;
3275 if(c1==c2) {fQ=q12; sQ=q23;}
3276 else {fQ=q23; sQ=q12;}
3277 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3278 }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
3279 if(term==1 || term==5){
3280 if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
3281 else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
3282 else {tQ=q23; sQ=q12; fQ=q13;}
3283 }else if(term==2 && part==1){
3284 fQ=q12;
3285 if(c1==c3) {tQ=q13; sQ=q23;}
3286 else {tQ=q23; sQ=q13;}
3287 }else if(term==2 && part==2){
3288 fQ=q13;
3289 if(c1==c2) {tQ=q12; sQ=q23;}
3290 else {tQ=q23; sQ=q12;}
3291 }else if(term==3 && part==1){
3292 sQ=q12;
3293 if(c1==c3) {tQ=q13; fQ=q23;}
3294 else {tQ=q23; fQ=q13;}
3295 }else if(term==3 && part==2){
3296 sQ=q13;
3297 if(c1==c2) {tQ=q12; fQ=q23;}
3298 else {tQ=q23; fQ=q12;}
3299 }else if(term==4 && part==1){
3300 tQ=q12; sQ=q13; fQ=q23;
3301 }else if(term==4 && part==2){
3302 tQ=q13; sQ=q12; fQ=q23;
3303 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3304 }else {// fQ=ss, sQ=ss, tQ=ss
3305 if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
3306 else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
3307 else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
3308 else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
3309 else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
3310 else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
3311 else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
3312 }
3313 }else if(fi != 5){// all the rest except pi-k-p
3314 if(key1==key2){
3315 fQ=q12;
3316 if(c1==c2){
3317 // cases not explicity shown below are not possible
3318 if(term==1 || term==5) {sQ=q13; tQ=q23;}
3319 else if(term==2 && part==1) {sQ=q13; tQ=q23;}
3320 else if(term==3 && part==2) {sQ=q13; tQ=q23;}
3321 else if(term==4 && part==2) {tQ=q13; sQ=q23;}
3322 else cout<<"problem!!!!!!!!!!!!!"<<endl;
3323 }else if(c3==0){
3324 if(c1==c3) {sQ=q13; tQ=q23;}
3325 else {sQ=q23; tQ=q13;}
3326 }else {//c3==1
3327 if(c1==c3) {tQ=q13; sQ=q23;}
3328 else {tQ=q23; sQ=q13;}
3329 }
3330 }else if(key1==key3){
3331 fQ=q13;
3332 if(c1==c3){
3333 // cases not explicity shown below are not possible
3334 if(term==1 || term==5) {sQ=q12; tQ=q23;}
3335 else if(term==2 && part==2) {sQ=q12; tQ=q23;}
3336 else if(term==3 && part==1) {sQ=q12; tQ=q23;}
3337 else if(term==4 && part==1) {tQ=q12; sQ=q23;}
3338 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3339 }else if(c2==0){
3340 if(c1==c2) {sQ=q12; tQ=q23;}
3341 else {sQ=q23; tQ=q12;}
3342 }else {//c2==1
3343 if(c1==c2) {tQ=q12; sQ=q23;}
3344 else {tQ=q23; sQ=q12;}
3345 }
3346 }else {// key2==key3
3347 fQ=q23;
3348 if(c2==c3){
3349 // cases not explicity shown below are not possible
3350 if(term==1 || term==5) {sQ=q12; tQ=q13;}
3351 else if(term==3 && part==1) {sQ=q12; tQ=q13;}
3352 else if(term==3 && part==2) {sQ=q13; tQ=q12;}
3353 else if(term==4 && part==1) {tQ=q12; sQ=q13;}
3354 else if(term==4 && part==2) {tQ=q13; sQ=q12;}
3355 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3356 }else if(c1==0){
3357 if(c1==c2) {sQ=q12; tQ=q13;}
3358 else {sQ=q13; tQ=q12;}
3359 }else {//c1==1
3360 if(c1==c2) {tQ=q12; sQ=q13;}
3361 else {tQ=q13; sQ=q12;}
3362 }
3363 }
3364 }else {// pi-k-p
3365 if(key1==1){
3366 if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
3367 else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
3368 }else if(key1==10){
3369 if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
3370 else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
3371 }else {// key1==100
3372 if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
3373 else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
3374 }
3375
3376 }
3377
3378
3379}
3380//________________________________________________________________________
3381Float_t AliChaoticity::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
3382
3383 Float_t qinv=1.0;
3384
3385 if(fi==0 || fi==3 || fi==5){// identical masses
3386 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));
3387 }else{// different masses
3388 Float_t px = track1[1] + track2[1];
3389 Float_t py = track1[2] + track2[2];
3390 Float_t pz = track1[3] + track2[3];
3391 Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
3392 Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
3393 deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
3394
3395 qinv = pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
3396 qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
3397 qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
3398 qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
3399 qinv = sqrt(qinv);
3400 }
3401
3402 return qinv;
3403
3404}
3405//________________________________________________________________________
3406void AliChaoticity::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
3407
3408 Float_t p0 = track1[0] + track2[0];
3409 Float_t px = track1[1] + track2[1];
3410 Float_t py = track1[2] + track2[2];
3411 Float_t pz = track1[3] + track2[3];
3412
3413 Float_t mt = sqrt(p0*p0 - pz*pz);
3414 Float_t pt = sqrt(px*px + py*py);
3415
3416 Float_t v0 = track1[0] - track2[0];
3417 Float_t vx = track1[1] - track2[1];
3418 Float_t vy = track1[2] - track2[2];
3419 Float_t vz = track1[3] - track2[3];
3420
3421 qout = (px*vx + py*vy)/pt;
3422 qside = (px*vy - py*vx)/pt;
3423 qlong = (p0*vz - pz*v0)/mt;
3424}
3425//________________________________________________________________________
90814457 3426//void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F ***histos){
e67e09ed 3427void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F *histos[3][10]){
3428//void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliChaoticity::kKbinsT][AliChaoticity::kCentBins]){
cd12341d 3429 for(Int_t mb=0; mb<fMbins; mb++){
3430 for(Int_t edB=0; edB<kEDbins; edB++){
3431 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3432 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3433 //
3434 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3435 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3436 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3437
3438 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3439 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3440 }
3441 }
3442 }
3443 }
3444 }
3445 //
3446 }
3447 }
3448
90814457 3449 if(legoCase){
3450 cout<<"LEGO call to SetWeightArrays"<<endl;
3451 for(Int_t mb=0; mb<fMbins; mb++){
3452 for(Int_t edB=0; edB<kEDbins; edB++){
3453 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3454 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3455 //
3456 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3457 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3458 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3459
3460 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinContent(qobin, qsbin, qlbin);
3461 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinError(qobin, qsbin, qlbin);
3462 }
cd12341d 3463 }
3464 }
3465 }
90814457 3466 }
3467 //
3468 }
3469 }
3470 }else{
3471
3472 TFile *wFile = new TFile("WeightFile.root","READ");
3473 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3474 else cout<<"Good Weight File Found!"<<endl;
3475
3476 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3477 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3478 for(Int_t mb=0; mb<fMbins; mb++){
3479 for(Int_t edB=0; edB<kEDbins; edB++){
3480
3481 TString *name = new TString("Weight_Kt_");
3482 *name += tKbin;
3483 name->Append("_Ky_");
3484 *name += yKbin;
3485 name->Append("_M_");
3486 *name += mb;
3487 name->Append("_ED_");
3488 *name += edB;
3489
3490 TH3F *fTempHisto = (TH3F*)wFile->Get(name->Data());
3491
3492 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3493 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3494 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3495
3496 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinContent(qobin, qsbin, qlbin);
3497 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinError(qobin, qsbin, qlbin);
3498 }
3499 }
3500 }
3501 delete fTempHisto;
3502 }//ED
3503 }//mb
3504 }//ky
3505 }//kt
3506
3507 wFile->Close();
3508 }
3509
5e3e77d6 3510 cout<<"Done reading weight file"<<endl;
cd12341d 3511
3512}
3513//________________________________________________________________________
cd12341d 3514void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3515
3516 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3517 Float_t ky=0;// central rapidity
3518 //
3519 Float_t qOut=0,qSide=0,qLong=0;
3520 GetQosl(track1, track2, qOut, qSide, qLong);
3521 qOut = fabs(qOut);
3522 qSide = fabs(qSide);
3523 qLong = fabs(qLong);
3524 //
3525
3526 if(kt < fKmeanT[0]) {fKtbinL=0; fKtbinH=0;}
3527 else if(kt >= fKmeanT[kKbinsT-1]) {fKtbinL=kKbinsT-1; fKtbinH=kKbinsT-1;}
3528 else {
3529 for(Int_t i=0; i<kKbinsT-1; i++){
3530 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtbinL=i; fKtbinH=i+1; break;}
3531 }
3532 }
3533 //
3534 if(ky < fKmeanY[0]) {fKybinL=0; fKybinH=0;}
3535 else if(ky >= fKmeanY[kKbinsY-1]) {fKybinL=kKbinsY-1; fKybinH=kKbinsY-1;}
3536 else {
3537 for(Int_t i=0; i<kKbinsY-1; i++){
3538 if((ky >= fKmeanY[i]) && (ky < fKmeanY[i+1])) {fKybinL=i; fKybinH=i+1; break;}
3539 }
3540 }
3541 //
3542 /////////
3543 if(qOut < fQmean[0]) {fQobinL=0; fQobinH=0;}
3544 else if(qOut >= fQmean[kQbinsWeights-1]) {fQobinL=kQbinsWeights-1; fQobinH=kQbinsWeights-1;}
3545 else {
3546 for(Int_t i=0; i<kQbinsWeights-1; i++){
3547 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQobinL=i; fQobinH=i+1; break;}
3548 }
3549 }
3550 //
3551 if(qSide < fQmean[0]) {fQsbinL=0; fQsbinH=0;}
3552 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsbinL=kQbinsWeights-1; fQsbinH=kQbinsWeights-1;}
3553 else {
3554 for(Int_t i=0; i<kQbinsWeights-1; i++){
3555 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsbinL=i; fQsbinH=i+1; break;}
3556 }
3557 }
3558 //
3559 if(qLong < fQmean[0]) {fQlbinL=0; fQlbinH=0;}
3560 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlbinL=kQbinsWeights-1; fQlbinH=kQbinsWeights-1;}
3561 else {
3562 for(Int_t i=0; i<kQbinsWeights-1; i++){
3563 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlbinL=i; fQlbinH=i+1; break;}
3564 }
3565 }
3566 //
3567
3568
3569 Float_t min = fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3570 Float_t minErr = fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3571
3572
3573 Float_t deltaW=0;
3574 // kt
3575 deltaW += (fNormWeight[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - min)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3576 // Ky
3577 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - min)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3578 // Qo
5e3e77d6 3579 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - min)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3580 // Qs
5e3e77d6 3581 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - min)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3582 // Ql
5e3e77d6 3583 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - min)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3584
3585 //
3586 wgt = min + deltaW;
3587
3588
3589 ////
3590
3591 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3592 Float_t deltaWErr=0;
3593 // Kt
3594 /*
3595 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - minErr)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3596 // Ky
3597 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - minErr)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3598 // Qo
5e3e77d6 3599 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - minErr)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3600 // Qs
5e3e77d6 3601 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - minErr)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3602 // Ql
5e3e77d6 3603 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - minErr)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3604 */
3605 wgtErr = minErr + deltaWErr;
3606
3607
3608
cd12341d 3609}
3610//________________________________________________________________________
3611Float_t AliChaoticity::MCWeight(Int_t charge1, Int_t charge2, Int_t rIndex, Int_t dampIndex, Float_t qinv){
3612
5e3e77d6 3613 Float_t radius = Float_t(rIndex+3.)/0.19733;// convert to GeV
cd12341d 3614 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
90814457 3615 Float_t coulCorr12 = FSICorrelationGaus2(charge1, charge2, rIndex, qinv);
cd12341d 3616 if(charge1==charge2){
5e3e77d6 3617 return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*radius,2)))*coulCorr12);
cd12341d 3618 }else {
5e3e77d6 3619 return ((1-myDamp) + myDamp*coulCorr12);
cd12341d 3620 }
3621
3622}
3623//________________________________________________________________________
90814457 3624Float_t AliChaoticity::MCWeight3D(Bool_t SameCharge, Int_t term, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23){
cd12341d 3625 if(term==5) return 1.0;
5e3e77d6 3626
90814457 3627 Float_t radius=5;
3628 if(fMbin<=1) {radius = 8;}
3629 else if(fMbin<=3) {radius = 7;}
3630 else if(fMbin<=5) {radius = 6;}
3631 else {radius = 5;}
3632 radius /= 0.19733;
3633
3634 //Float_t radius = (3. + rIndex)/0.19733;//starts at 3fm. convert to GeV
cd12341d 3635 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3636 Float_t fc = sqrt(myDamp);
5e3e77d6 3637 if(SameCharge){// all three of the same charge
90814457 3638 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2
3639 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, q13);// K2
3640 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, q23);// K2
cd12341d 3641
5e3e77d6 3642 if(term==1){
3643 Float_t c3QS = 1 + exp(-pow(q12*radius,2)) + exp(-pow(q13*radius,2)) + exp(-pow(q23*radius,2));
3644 c3QS += 2*exp(-pow(radius,2)*(pow(q12,2) + pow(q13,2) + pow(q23,2))/2.);
3645 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3646 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3647 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*radius,2)))*coulCorr13;
3648 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*radius,2)))*coulCorr23;
3649 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kTRUE, q12, q13, q23);
3650 return w123;
3651 }else if(term==2){
3652 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3653 }else if(term==3){
3654 return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*radius,2)))*coulCorr13);
3655 }else if(term==4){
3656 return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*radius,2)))*coulCorr23);
3657 }else return 1.0;
3658
3659 }else{// mixed charge case pair 12 always treated as ss
90814457 3660 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2 ss
3661 Float_t coulCorr13 = FSICorrelationTherm2(+1,-1, q13);// K2 os
3662 Float_t coulCorr23 = FSICorrelationTherm2(+1,-1, q23);// K2 os
5e3e77d6 3663 if(term==1){
3664 Float_t c3QS = 1 + exp(-pow(q12*radius,2));
3665 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3666 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3667 w123 += pow(fc,2)*(1-fc)*coulCorr13;
3668 w123 += pow(fc,2)*(1-fc)*coulCorr23;
3669 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kFALSE, q12, q13, q23);
3670 return w123;
3671 }else if(term==2){
3672 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3673 }else if(term==3){
3674 return ((1-myDamp) + myDamp*coulCorr13);
3675 }else if(term==4){
3676 return ((1-myDamp) + myDamp*coulCorr23);
3677 }else return 1.0;
3678 }
3679
cd12341d 3680}
3681//________________________________________________________________________
90814457 3682void AliChaoticity::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D){
5e3e77d6 3683
3684
3685 if(legoCase){
90814457 3686 cout<<"LEGO call to SetMomResCorrections"<<endl;
5e3e77d6 3687 fMomResC2 = (TH2D*)temp2D->Clone();
5e3e77d6 3688 fMomResC2->SetDirectory(0);
5e3e77d6 3689 }else {
3690 TFile *momResFile = new TFile("MomResFile.root","READ");
3691 if(!momResFile->IsOpen()) {
3692 cout<<"No momentum resolution file found"<<endl;
3693 AliFatal("No momentum resolution file found. Kill process.");
3694 }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
3695
3696 TH2D *temp2D2 = (TH2D*)momResFile->Get("MomResHisto_pp");
5e3e77d6 3697 fMomResC2 = (TH2D*)temp2D2->Clone();
5e3e77d6 3698 fMomResC2->SetDirectory(0);
90814457 3699
5e3e77d6 3700 momResFile->Close();
3701 }
cd12341d 3702
5e3e77d6 3703 cout<<"Done reading momentum resolution file"<<endl;
3704}
3705//________________________________________________________________________
90814457 3706void AliChaoticity::SetFSICorrelations(Bool_t legoCase, TH2D *temp2DGaus[2], TH2D *temp2DTherm[2], TH3D *temp3Dos[6], TH3D *temp3Dss[6]){
5e3e77d6 3707 // read in 2-particle and 3-particle FSI correlations = K2 & K3
3708 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
3709 if(legoCase){
c4980714 3710 cout<<"LEGO call to SetFSICorrelations"<<endl;
90814457 3711 fFSI2SS[0] = (TH2D*)temp2DGaus[0]->Clone();
3712 fFSI2OS[0] = (TH2D*)temp2DGaus[1]->Clone();
3713 fFSI2SS[1] = (TH2D*)temp2DTherm[0]->Clone();
3714 fFSI2OS[1] = (TH2D*)temp2DTherm[1]->Clone();
5e3e77d6 3715 //
90814457 3716 fFSI2SS[0]->SetDirectory(0);
3717 fFSI2OS[0]->SetDirectory(0);
3718 fFSI2SS[1]->SetDirectory(0);
3719 fFSI2OS[1]->SetDirectory(0);
3720
3721 for(Int_t CB=0; CB<6; CB++) {
3722 fFSIOmega0OS[CB] = (TH3D*)temp3Dos[CB]->Clone();
3723 fFSIOmega0SS[CB] = (TH3D*)temp3Dss[CB]->Clone();
3724 //
3725 fFSIOmega0OS[CB]->SetDirectory(0);
3726 fFSIOmega0SS[CB]->SetDirectory(0);
3727 }
5e3e77d6 3728 }else {
c4980714 3729 cout<<"non LEGO call to SetFSICorrelations"<<endl;
5e3e77d6 3730 TFile *fsifile = new TFile("KFile.root","READ");
3731 if(!fsifile->IsOpen()) {
3732 cout<<"No FSI file found"<<endl;
3733 AliFatal("No FSI file found. Kill process.");
3734 }else {cout<<"Good FSI File Found!"<<endl;}
3735
90814457 3736 TH2D *temphisto2GausSS = (TH2D*)fsifile->Get("K2ssG");
3737 TH2D *temphisto2GausOS = (TH2D*)fsifile->Get("K2osG");
3738 TH2D *temphisto2ThermSS = (TH2D*)fsifile->Get("K2ssT");
3739 TH2D *temphisto2ThermOS = (TH2D*)fsifile->Get("K2osT");
3740 TH3D *temphisto3OS[6];
ae9b34d1 3741 TH3D *temphisto3SS[6];
3742 for(Int_t CB=0; CB<6; CB++) {
90814457 3743 TString *nameK3SS = new TString("K3ss_");
3744 *nameK3SS += CB;
3745 temphisto3SS[CB] = (TH3D*)fsifile->Get(nameK3SS->Data());
3746 //
3747 TString *nameK3OS = new TString("K3os_");
3748 *nameK3OS += CB;
3749 temphisto3OS[CB] = (TH3D*)fsifile->Get(nameK3OS->Data());
ae9b34d1 3750 }
3751
90814457 3752 fFSI2SS[0] = (TH2D*)temphisto2GausSS->Clone();
3753 fFSI2OS[0] = (TH2D*)temphisto2GausOS->Clone();
3754 fFSI2SS[1] = (TH2D*)temphisto2ThermSS->Clone();
3755 fFSI2OS[1] = (TH2D*)temphisto2ThermOS->Clone();
3756 fFSI2SS[0]->SetDirectory(0);
3757 fFSI2OS[0]->SetDirectory(0);
3758 fFSI2SS[1]->SetDirectory(0);
3759 fFSI2OS[1]->SetDirectory(0);
ae9b34d1 3760
90814457 3761 for(Int_t CB=0; CB<6; CB++) {
3762 fFSIOmega0SS[CB] = (TH3D*)temphisto3SS[CB]->Clone();
3763 fFSIOmega0OS[CB] = (TH3D*)temphisto3OS[CB]->Clone();
3764 fFSIOmega0SS[CB]->SetDirectory(0);
3765 fFSIOmega0OS[CB]->SetDirectory(0);
3766 }
3767 //
3768
5e3e77d6 3769 fsifile->Close();
cd12341d 3770 }
3771
5e3e77d6 3772 // condition FSI histogram for edge effects
ae9b34d1 3773 for(Int_t CB=0; CB<6; CB++){
3774 for(Int_t ii=1; ii<=fFSIOmega0SS[CB]->GetNbinsX(); ii++){
3775 for(Int_t jj=1; jj<=fFSIOmega0SS[CB]->GetNbinsY(); jj++){
3776 for(Int_t kk=1; kk<=fFSIOmega0SS[CB]->GetNbinsZ(); kk++){
5e3e77d6 3777
ae9b34d1 3778 if(fFSIOmega0SS[CB]->GetBinContent(ii,jj,kk) <=0){
3779 Double_t Q12 = fFSIOmega0SS[CB]->GetXaxis()->GetBinCenter(ii);
3780 Double_t Q23 = fFSIOmega0SS[CB]->GetYaxis()->GetBinCenter(jj);
3781 Double_t Q13 = fFSIOmega0SS[CB]->GetZaxis()->GetBinCenter(kk);
3782 //
3783 Int_t Q12bin=ii;
3784 Int_t Q23bin=jj;
3785 Int_t Q13bin=kk;
3786 Int_t AC=0;//Adjust Counter
3787 Int_t AClimit=10;// maximum bin shift
3788 if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin++; AC++;}}
3789 if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin--; AC++;}}
3790 //
3791 if(Q13 < sqrt(pow(Q12,2)+pow(Q23,2) - 2*Q12*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin++; AC++;}}
3792 if(Q13 > sqrt(pow(Q12,2)+pow(Q23,2) + 2*Q12*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin--; AC++;}}
3793 //
3794 if(Q23 < sqrt(pow(Q12,2)+pow(Q13,2) - 2*Q12*Q13)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin++; AC++;}}
3795 if(Q23 > sqrt(pow(Q12,2)+pow(Q13,2) + 2*Q12*Q13)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin--; AC++;}}
3796
3797 // Save cpu time by setting empty cell contents (edge effects) to nearest non-zero cell (these cells are not used very often anyway.)
3798 if(AC==AClimit) {
3799 fFSIOmega0SS[CB]->SetBinContent(ii,jj,kk, 1.0);
90814457 3800 fFSIOmega0OS[CB]->SetBinContent(ii,jj,kk, 1.0);
ae9b34d1 3801 }else {
3802 fFSIOmega0SS[CB]->SetBinContent(ii,jj,kk, fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin));
90814457 3803 fFSIOmega0OS[CB]->SetBinContent(ii,jj,kk, fFSIOmega0OS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin));
ae9b34d1 3804 }
5e3e77d6 3805 }
ae9b34d1 3806
5e3e77d6 3807 }
5e3e77d6 3808 }
cd12341d 3809 }
3810 }
5e3e77d6 3811 cout<<"Done reading FSI file"<<endl;
cd12341d 3812}
3813//________________________________________________________________________
90814457 3814Float_t AliChaoticity::FSICorrelationGaus2(Int_t charge1, Int_t charge2, Int_t rIndex, Float_t qinv){
5e3e77d6 3815 // returns 2-particle Coulomb correlations = K2
3816 if(rIndex >= kRVALUES) return 1.0;
90814457 3817 Int_t qbinL = fFSI2SS[0]->GetYaxis()->FindBin(qinv-fFSI2SS[0]->GetYaxis()->GetBinWidth(1)/2.);
3818 Int_t qbinH = qbinL+1;
3819 if(qbinL <= 0) return 1.0;
3820 if(qbinH > fFSI2SS[0]->GetNbinsY()) return 1.0;
3821
3822 Float_t slope=0;
3823 if(charge1==charge2){
3824 slope = fFSI2SS[0]->GetBinContent(rIndex+1, qbinL) - fFSI2SS[0]->GetBinContent(rIndex+1, qbinH);
3825 slope /= fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinH);
3826 return (slope*(qinv - fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS[0]->GetBinContent(rIndex+1, qbinL));
3827 }else {
3828 slope = fFSI2OS[0]->GetBinContent(rIndex+1, qbinL) - fFSI2OS[0]->GetBinContent(rIndex+1, qbinH);
3829 slope /= fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinH);
3830 return (slope*(qinv - fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS[0]->GetBinContent(rIndex+1, qbinL));
3831 }
3832}
3833//________________________________________________________________________
3834Float_t AliChaoticity::FSICorrelationTherm2(Int_t charge1, Int_t charge2, Float_t qinv){
3835 // returns 2-particle Coulomb correlations = K2
3836 Int_t qbinL = fFSI2SS[1]->GetYaxis()->FindBin(qinv-fFSI2SS[1]->GetYaxis()->GetBinWidth(1)/2.);
5e3e77d6 3837 Int_t qbinH = qbinL+1;
3838 if(qbinL <= 0) return 1.0;
90814457 3839 if(qbinH > fFSI2SS[1]->GetNbinsY()) return 1.0;
5e3e77d6 3840
3841 Float_t slope=0;
3842 if(charge1==charge2){
90814457 3843 slope = fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinL) - fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinH);
3844 slope /= fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinH);
3845 return (slope*(qinv - fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinL));
5e3e77d6 3846 }else {
90814457 3847 slope = fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinL) - fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinH);
3848 slope /= fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinH);
3849 return (slope*(qinv - fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinL));
cd12341d 3850 }
3851}
3852//________________________________________________________________________
5e3e77d6 3853Double_t AliChaoticity::FSICorrelationOmega0(Bool_t SameCharge, Double_t Q12, Double_t Q13, Double_t Q23){
3854 // remember: Omega0 histogram is in the following order (Q12, Q23, Q13)
3855 // returns 3d 3-particle Coulomb Correlation = K3
ae9b34d1 3856 Int_t Q12bin = fFSIOmega0SS[fFSIbin]->GetXaxis()->FindBin(Q12);
3857 Int_t Q13bin = fFSIOmega0SS[fFSIbin]->GetZaxis()->FindBin(Q13);
3858 Int_t Q23bin = fFSIOmega0SS[fFSIbin]->GetYaxis()->FindBin(Q23);
5e3e77d6 3859
3860 if(SameCharge){
ae9b34d1 3861 if(fFSIOmega0SS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3862 else return fFSIOmega0SS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
90814457 3863 }else{// mixed charge.
3864 if(fFSIOmega0OS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3865 else return fFSIOmega0OS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
5e3e77d6 3866 }
cd12341d 3867}
90814457 3868//________________________________________________________________________
3869void 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){
3870 QS1v1 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
3871 QS1v1 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
3872 QS1v1 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
3873 QS2 = (pV1[1]-pV2[1])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[1]-pV3[1]);
3874 QS3v1 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
3875 QS3v1 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
3876 //
3877 QS1v2 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
3878 QS1v2 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
3879 QS3v2 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
3880 QS3v2 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
3881 QS3v2 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
3882}
3883//________________________________________________________________________
3884void AliChaoticity::TestAddTask(){
3885 /*
3886 TString inputFileNameFSI = "KFile.root";
3887 TFile *inputFileFSI = TFile::Open(inputFileNameFSI,"OLD");
3888 if (!inputFileFSI){
3889 cout << "Requested file:" << inputFileFSI << " was not opened. ABORT." << endl;
3890 return;
3891 }
3892 TH2D *FSI2gaus[2];
3893 TH2D *FSI2therm[2];
3894 TH3D *FSI3ss[6];
3895 TH3D *FSI3os[6];
3896 FSI2gaus[0] = (TH2D*)inputFileFSI->Get("K2ssG");
3897 FSI2gaus[1] = (TH2D*)inputFileFSI->Get("K2osG");
3898 FSI2therm[0] = (TH2D*)inputFileFSI->Get("K2ssT");
3899 FSI2therm[1] = (TH2D*)inputFileFSI->Get("K2osT");
3900 for(Int_t CB=0; CB<6; CB++) {
3901 TString *nameSS=new TString("K3ss_");
3902 *nameSS += CB;
3903 FSI3ss[CB] = (TH3D*)inputFileFSI->Get(nameSS->Data());
3904 TString *nameOS=new TString("K3os_");
3905 *nameOS += CB;
3906 FSI3os[CB] = (TH3D*)inputFileFSI->Get(nameOS->Data());
3907 }
3908 //
3909 FSI2gaus[0]->SetDirectory(0);
3910 FSI2gaus[1]->SetDirectory(0);
3911 FSI2therm[0]->SetDirectory(0);
3912 FSI2therm[1]->SetDirectory(0);
3913 for(Int_t CB=0; CB<6; CB++) {
3914 FSI3ss[CB]->SetDirectory(0);
3915 FSI3os[CB]->SetDirectory(0);
3916 }
3917 SetFSICorrelations( kTRUE, FSI2gaus, FSI2therm , FSI3os, FSI3ss);
3918 //
3919
3920 if(!fTabulatePairs) {
3921 TString inputFileNameWeight = "WeightFile.root";
3922 TFile *inputFileWeight = TFile::Open(inputFileNameWeight,"OLD");
3923 if (!inputFileWeight){
3924 cout << "Requested file:" << inputFileWeight << " was not opened. ABORT." << endl;
3925 return;
3926 }
3927
3928 ////////////////////////////////////////////////////
3929 // C2 Weight File
3930 Int_t ktbins = GetNumKtbins();
3931 Int_t cbins = GetNumCentBins();
3932 TH3F *weightHisto[ktbins][cbins];
3933 for(Int_t i=0; i<ktbins; i++){
3934 for(Int_t j=0; j<cbins; j++){
3935 TString name = "Weight_Kt_";
3936 name += i;
3937 name += "_Ky_0_M_";
3938 name += j;
3939 name += "_ED_0";
3940
3941 weightHisto[i][j] = (TH3F*)inputFileWeight->Get(name);
3942 }
3943 }
3944 SetWeightArrays( kTRUE, weightHisto );
3945 }
3946
3947 //
3948 if(!fMCcase && !fTabulatePairs){
3949 TString inputFileNameMomRes = "MomResFile.root";
3950 TFile *inputFileMomRes = TFile::Open(inputFileNameMomRes,"OLD");
3951 if (!inputFileMomRes){
3952 cout << "Requested file:" << inputFileMomRes << " was not opened. ABORT." << endl;
3953 return;
3954 }
3955 ////////////////////////////////////////////////////
3956 // Momentum Resolution File
3957 TH2D *momResHisto2D = 0;
3958 momResHisto2D = (TH2D*)inputFileMomRes->Get("MomResHisto_pp");
3959 SetMomResCorrections( kTRUE, momResHisto2D);
3960 }
3961 */
3962}