]> 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
Coverity fix
[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);
704f2481 964 TString *nameAvgP=new TString(nameEx2->Data());
965 nameAvgP->Append("_AvgP");
966 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fAvgP = new TProfile2D(nameAvgP->Data(),"",10,0,1, 400,0,2, 0.,1.0,"");
967 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fAvgP);
968
cd12341d 969 // Momentum resolution histos
970 if(fMCcase && sc==0){
971 TString *nameIdeal = new TString(nameEx2->Data());
972 nameIdeal->Append("_Ideal");
fa109294 973 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 974 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal);
975 TString *nameSmeared = new TString(nameEx2->Data());
976 nameSmeared->Append("_Smeared");
fa109294 977 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 978 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared);
fa109294 979 //
980 TString *nameEx2MC=new TString(nameEx2->Data());
981 nameEx2MC->Append("_MCqinv");
982 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"",400,0,2);
983 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
984 TString *nameEx2MCQW=new TString(nameEx2->Data());
985 nameEx2MCQW->Append("_MCqinvQW");
986 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"",400,0,2);
987 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
54d66278 988 //
989 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
990 nameEx2PIDpurityDen->Append("_PIDpurityDen");
991 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);
992 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
993 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
994 nameEx2PIDpurityNum->Append("_PIDpurityNum");
995 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);
996 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
cd12341d 997 }
fa109294 998 if(sc==0){
cd12341d 999
1000 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
1001 nameEx2OSLB1->Append("_osl_b1");
fa109294 1002 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 1003 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSL);
1004 nameEx2OSLB1->Append("_QW");
fa109294 1005 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 1006 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fExplicit2OSLQW);
1007 //
1008 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
1009 nameEx2OSLB2->Append("_osl_b2");
fa109294 1010 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 1011 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSL);
1012 nameEx2OSLB2->Append("_QW");
fa109294 1013 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 1014 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fExplicit2OSLQW);
1015
1016 }
1017
1018 }// term_2
1019 }// SC_2
1020
1021 // skip 3-particle if Tabulate6DPairs is true
1022 if(fTabulatePairs) continue;
1023
1024 for(Int_t c3=0; c3<2; c3++){
1025 for(Int_t sc=0; sc<kSCLimit3; sc++){
1026 for(Int_t term=0; term<5; term++){
1027 TString *nameEx3 = new TString("Explicit3_Charge1_");
1028 *nameEx3 += c1;
1029 nameEx3->Append("_Charge2_");
1030 *nameEx3 += c2;
1031 nameEx3->Append("_Charge3_");
1032 *nameEx3 += c3;
1033 nameEx3->Append("_SC_");
1034 *nameEx3 += sc;
1035 nameEx3->Append("_M_");
1036 *nameEx3 += mb;
1037 nameEx3->Append("_ED_");
1038 *nameEx3 += edB;
1039 nameEx3->Append("_Term_");
1040 *nameEx3 += term+1;
1041
1042 TString *namePC3 = new TString("PairCut3_Charge1_");
1043 *namePC3 += c1;
1044 namePC3->Append("_Charge2_");
1045 *namePC3 += c2;
1046 namePC3->Append("_Charge3_");
1047 *namePC3 += c3;
1048 namePC3->Append("_SC_");
1049 *namePC3 += sc;
1050 namePC3->Append("_M_");
1051 *namePC3 += mb;
1052 namePC3->Append("_ED_");
1053 *namePC3 += edB;
1054 namePC3->Append("_Term_");
1055 *namePC3 += term+1;
1056
1057 ///////////////////////////////////////
1058 // skip degenerate histograms
1059 if(sc==0 || sc==6 || sc==9){// Identical species
1060 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
1061 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
1062 }else if(sc!=5){
1063 if( (c1+c2)==1) {if(c1!=0) continue;}
1064 }else {}// do nothing for pi-k-p case
1065
1066 /////////////////////////////////////////
1067
1068 if(fPdensityExplicitLoop){
1069 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);
1070 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3);
1071 //
1072 nameEx3->Append("_Norm");
1073 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);
1074 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3);
1075 }
1076 if(fPdensityPairCut){
1077 TString *nameNorm=new TString(namePC3->Data());
1078 nameNorm->Append("_Norm");
1079 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);
1080 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1081 //
1082 if(sc<=2){
1083 TString *name3DQ=new TString(namePC3->Data());
1084 name3DQ->Append("_3D");
1085 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);
1086 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1087 //
90814457 1088 const int NEdgesPos=16;
1089 double lowEdges4vectPos[NEdgesPos]={0};
1090 lowEdges4vectPos[0]=0.0;
1091 lowEdges4vectPos[1]=0.0001;// best resolution at low Q^2
1092 for(int edge=2; edge<NEdgesPos; edge++){
1093 lowEdges4vectPos[edge] = lowEdges4vectPos[edge-1] + lowEdges4vectPos[1]*(edge);
1094 }
1095 const int NEdges=2*NEdgesPos-1;
fa109294 1096 double lowEdges4vect[NEdges]={0};
90814457 1097 for(int edge=0; edge<NEdges; edge++){
1098 if(edge<NEdgesPos-1) lowEdges4vect[edge] = -lowEdges4vectPos[NEdgesPos-1-edge];
1099 else if(edge==NEdgesPos-1) lowEdges4vect[edge] = 0;
1100 else lowEdges4vect[edge] = lowEdges4vectPos[edge-NEdgesPos+1];
fa109294 1101 }
90814457 1102
5e3e77d6 1103 if(c1==c2 && c1==c3 && sc==0 && fMCcase==kFALSE){
90814457 1104 TString *name4vect1=new TString(namePC3->Data());
1105 TString *name4vect2=new TString(namePC3->Data());
1106 name4vect1->Append("_4VectProd1");
1107 name4vect2->Append("_4VectProd2");
5e3e77d6 1108 // use 3.75e6 MeV^4 as the resolution on QprodSum
90814457 1109 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);
1110 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1Terms);
1111 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);
1112 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2Terms);
5e3e77d6 1113 }
1114 if(sc==0 && fMCcase==kTRUE){
1115 TString *name3DMomResIdeal=new TString(namePC3->Data());
1116 name3DMomResIdeal->Append("_Ideal");
1117 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);
1118 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1119 TString *name3DMomResSmeared=new TString(namePC3->Data());
1120 name3DMomResSmeared->Append("_Smeared");
1121 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);
1122 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared);
601fbb33 1123 //
1124 TString *name3DMomResQW12=new TString(namePC3->Data());
1125 name3DMomResQW12->Append("_QW12");
1126 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);
1127 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW12);
1128 TString *name3DMomResQW13=new TString(namePC3->Data());
1129 name3DMomResQW13->Append("_QW13");
1130 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);
1131 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fQW13);
90814457 1132 //
1133 if(term==0){
1134 TString *name3DSumK3=new TString(namePC3->Data());
1135 name3DSumK3->Append("_SumK3");
1136 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);
1137 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSumK3);
1138 TString *name3DEnK3=new TString(namePC3->Data());
1139 name3DEnK3->Append("_EnK3");
1140 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);
1141 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fEnK3);
1142 }
1143
1144 if(c1==c2 && c1==c3){
1145 TString *name4vect1Ideal=new TString(namePC3->Data());
1146 TString *name4vect1Smeared=new TString(namePC3->Data());
1147 TString *name4vect2Ideal=new TString(namePC3->Data());
1148 TString *name4vect2Smeared=new TString(namePC3->Data());
1149 name4vect1Ideal->Append("_4VectProd1Ideal");
1150 name4vect1Smeared->Append("_4VectProd1Smeared");
1151 name4vect2Ideal->Append("_4VectProd2Ideal");
1152 name4vect2Smeared->Append("_4VectProd2Smeared");
1153 // use 3.75e6 MeV^4 as the resolution on QprodSum
1154 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);
1155 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsIdeal);
1156 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);
1157 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSmeared);
1158 //
1159 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);
1160 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsIdeal);
1161 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);
1162 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSmeared);
1163 //
1164 if(term==0){
1165 TString *name4vect1SumK3=new TString(namePC3->Data());
1166 TString *name4vect2SumK3=new TString(namePC3->Data());
1167 TString *name4vect1EnK3=new TString(namePC3->Data());
1168 TString *name4vect2EnK3=new TString(namePC3->Data());
1169 name4vect1SumK3->Append("_4VectProd1SumK3");
1170 name4vect2SumK3->Append("_4VectProd2SumK3");
1171 name4vect1EnK3->Append("_4VectProd1EnK3");
1172 name4vect2EnK3->Append("_4VectProd2EnK3");
1173 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);
1174 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK3);
1175 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);
1176 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK3);
1177 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);
1178 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK3);
1179 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);
1180 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK3);
1181 }// term 0
1182 if(term > 0 && term < 4){
1183 TString *name4vect1SumK2=new TString(namePC3->Data());
1184 TString *name4vect2SumK2=new TString(namePC3->Data());
1185 TString *name4vect1EnK2=new TString(namePC3->Data());
1186 TString *name4vect2EnK2=new TString(namePC3->Data());
1187 name4vect1SumK2->Append("_4VectProd1SumK2");
1188 name4vect2SumK2->Append("_4VectProd2SumK2");
1189 name4vect1EnK2->Append("_4VectProd1EnK2");
1190 name4vect2EnK2->Append("_4VectProd2EnK2");
1191 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);
1192 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsSumK2);
1193 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);
1194 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsSumK2);
1195 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);
1196 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd1TermsEnK2);
1197 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);
1198 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].f4VectProd2TermsEnK2);
1199 }// terms 1,2,3
1200 }
1201 }// MCcase
5e3e77d6 1202 //
90814457 1203 if(c1==c2 && c1==c3 && term==4 && sc==0){
cd12341d 1204 for(Int_t dt=0; dt<kDENtypes; dt++){
1205 TString *nameDenType=new TString("PairCut3_Charge1_");
1206 *nameDenType += c1;
1207 nameDenType->Append("_Charge2_");
1208 *nameDenType += c2;
1209 nameDenType->Append("_Charge3_");
1210 *nameDenType += c3;
1211 nameDenType->Append("_SC_");
1212 *nameDenType += sc;
1213 nameDenType->Append("_M_");
1214 *nameDenType += mb;
1215 nameDenType->Append("_ED_");
1216 *nameDenType += edB;
1217 nameDenType->Append("_TPN_");
1218 *nameDenType += dt;
1219
1220 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);
1221 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNorm);
1222 // neglect errors for TPN
1223 //nameDenType->Append("_Err");
1224 //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);
1225 //fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].fTwoPartNormErr);
5e3e77d6 1226 //
fa109294 1227 TString *name4vect1TPN=new TString(nameDenType->Data());
1228 TString *name4vect2TPN=new TString(nameDenType->Data());
1229 name4vect1TPN->Append("_4VectProd1");
1230 name4vect2TPN->Append("_4VectProd2");
1231 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);
1232 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNorm);
1233 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);
1234 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNorm);
90814457 1235 //
1236 if(fMCcase){
1237 TString *name4vect1TPNIdeal=new TString(nameDenType->Data());
1238 TString *name4vect2TPNIdeal=new TString(nameDenType->Data());
1239 TString *name4vect1TPNSmeared=new TString(nameDenType->Data());
1240 TString *name4vect2TPNSmeared=new TString(nameDenType->Data());
1241 name4vect1TPNIdeal->Append("_4VectProd1Ideal");
1242 name4vect2TPNIdeal->Append("_4VectProd2Ideal");
1243 name4vect1TPNSmeared->Append("_4VectProd1Smeared");
1244 name4vect2TPNSmeared->Append("_4VectProd2Smeared");
1245 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);
1246 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormIdeal);
1247 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);
1248 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormIdeal);
1249 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);
1250 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd1TwoPartNormSmeared);
1251 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);
1252 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].DT[dt].f4VectProd2TwoPartNormSmeared);
1253 }
1254
cd12341d 1255 }
1256
1257 }// term=4
1258 }// c and sc exclusion
1259 }// PdensityPairCut
1260 }// term_3
1261 }// SC_3
1262 }//c3
1263 }//c2
1264 }//c1
1265 }// ED
1266 }// mbin
1267 }// Pdensity Method
1268
1269
1270 if(fTabulatePairs){
1271
1272 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
1273 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
1274 for(Int_t mb=0; mb<fMbins; mb++){
1275 for(Int_t edB=0; edB<kEDbins; edB++){
1276
1277 TString *nameNum = new TString("TwoPart_num_Kt_");
1278 *nameNum += tKbin;
1279 nameNum->Append("_Ky_");
1280 *nameNum += yKbin;
1281 nameNum->Append("_M_");
1282 *nameNum += mb;
1283 nameNum->Append("_ED_");
1284 *nameNum += edB;
1285
1286 TString *nameDen = new TString("TwoPart_den_Kt_");
1287 *nameDen += tKbin;
1288 nameDen->Append("_Ky_");
1289 *nameDen += yKbin;
1290 nameDen->Append("_M_");
1291 *nameDen += mb;
1292 nameDen->Append("_ED_");
1293 *nameDen += edB;
1294
1295
1296 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD = new TH3I(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1297 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fExplicit2ThreeD);
1298
1299 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD = new TH3I(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1300 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fExplicit2ThreeD);
1301 }
1302 }
1303 }
1304 }
1305
1306 }
1307
1308
1309 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1310 fOutputList->Add(fQsmearMean);
1311 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1312 fOutputList->Add(fQsmearSq);
1313 TH1D *fQDist = new TH1D("fQDist","",200,-.2,.2);
1314 fOutputList->Add(fQDist);
1315
1316
1317 ////////////////////////////////////
1318 ///////////////////////////////////
1319
1320 PostData(1, fOutputList);
1321}
1322
1323//________________________________________________________________________
1324void AliChaoticity::Exec(Option_t *)
1325{
1326 // Main loop
1327 // Called for each event
1ccd6f0d 1328 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
cd12341d 1329 fEventCounter++;
1330
b6e5ec54 1331 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
cd12341d 1332
b6e5ec54 1333 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1334 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
cd12341d 1335
cd12341d 1336
1337 // Trigger Cut
1338 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1339 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
fa109294 1340 if(!isSelected1 && !fMCcase) {return;}
5e3e77d6 1341 }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
cd12341d 1342 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1343 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
fa109294 1344 if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1345 }else {return;}
cd12341d 1346
1347 ///////////////////////////////////////////////////////////
1348 const AliAODVertex *primaryVertexAOD;
1349 AliCentrality *centrality;// for AODs and ESDs
1350
1351
1352 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1353 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1354 fPIDResponse = inputHandler->GetPIDResponse();
1355
1356
1357 TClonesArray *mcArray = 0x0;
1358 if(fMCcase){
1359 if(fAODcase){
1360 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1361 if(!mcArray || mcArray->GetEntriesFast() >= 110000){
1362 cout<<"No MC particle branch found or Array too large!!"<<endl;
cd12341d 1363 return;
1364 }
1365 }
1366 }
1367
1368
1369 UInt_t status=0;
1370 Int_t positiveTracks=0, negativeTracks=0;
1371 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1372
1373 Double_t vertex[3]={0};
1374 Int_t zbin=0;
1375 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1376 /////////////////////////////////////////////////
1377
1378
1379 Float_t centralityPercentile=0;
1380 Float_t cStep=5.0, cStart=0;
1381
1382
1383 if(fAODcase){// AOD case
1384
1385 if(fPbPbcase){
1386 centrality = fAOD->GetCentrality();
1387 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1388 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
5e3e77d6 1389 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range. Skipping Event"<<endl;*/ return;}
cd12341d 1390 cout<<"Centrality % = "<<centralityPercentile<<endl;
1391 }
1392
1393
1394
1395
1396 ////////////////////////////////
1397 // Vertexing
1398 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1399 primaryVertexAOD = fAOD->GetPrimaryVertex();
1400 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1401
1402 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1403 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1404
1405 if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Reject Pile-up events
1406 if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1407
1408 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1409
1410 fBfield = fAOD->GetMagneticField();
1411
1412 for(Int_t i=0; i<fZvertexBins; i++){
1413 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1414 zbin=i;
1415 break;
1416 }
1417 }
1418
1419
1420
1421 /////////////////////////////
1422 // Create Shuffled index list
1423 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1424 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1425 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1426 /////////////////////////////
1427
1428 // Track loop
1429 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1430 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1431 if (!aodtrack) continue;
1432 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1433
1434 status=aodtrack->GetStatus();
1435
1436 if((Bool_t)(((1<<7) & aodtrack->GetFilterMap()) == 0)) continue;// AOD filterBit cut
1437 // 1<<5 is for "standard cuts with tight dca cut"
1438 // 1<<7 is for TPC only tracking
1439 if(aodtrack->Pt() < 0.16) continue;
1440 if(fabs(aodtrack->Eta()) > 0.8) continue;
1441
1442
1443 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1444 if(!goodMomentum) continue;
1445 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1446
1447 Float_t dca2[2];
1448 Float_t dca3d;
1449
1450 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1451 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1452 dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1453
1454 fTempStruct[myTracks].fStatus = status;
1455 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1456 fTempStruct[myTracks].fId = aodtrack->GetID();
1457 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1458 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1459 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1460 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1461 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1462 fTempStruct[myTracks].fEta = aodtrack->Eta();
1463 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1464 fTempStruct[myTracks].fDCAXY = dca2[0];
1465 fTempStruct[myTracks].fDCAZ = dca2[1];
1466 fTempStruct[myTracks].fDCA = dca3d;
1467 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1468 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1469
1470
1471
1472 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1473 if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1474 if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1475
1476 if(fTempStruct[myTracks].fCharge==+1) {
1477 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1478 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1479 }else {
1480 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1481 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1482 }
1483
1484 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1485 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1486
1487
1488 // nSimga PID workaround
1489 fTempStruct[myTracks].fElectron = kFALSE;
1490 fTempStruct[myTracks].fPion = kFALSE;
1491 fTempStruct[myTracks].fKaon = kFALSE;
1492 fTempStruct[myTracks].fProton = kFALSE;
1493
1494 Float_t nSigmaTPC[5];
1495 Float_t nSigmaTOF[5];
1496 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1497 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1498 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1499 Float_t signalTPC=0, signalTOF=0;
1500 Double_t integratedTimesTOF[10]={0};
1501 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1502 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1503 if (!aodTrack2) continue;
1504 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1505
1506 UInt_t status2=aodTrack2->GetStatus();
1507
1508 nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
1509 nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
1510 nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
1511 nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
1512 nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
1513 //
1514 nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
1515 nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
1516 nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
1517 nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
1518 nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
1519 signalTPC = aodTrack2->GetTPCsignal();
1520
1521 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1522 fTempStruct[myTracks].fTOFhit = kTRUE;
1523 signalTOF = aodTrack2->GetTOFsignal();
1524 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1525 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1526
1527 }// aodTrack2
5e3e77d6 1528
c4980714 1529 //cout<<nSigmaTPC[2]<<endl;
cd12341d 1530 ///////////////////
1531 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1532 if(fTempStruct[myTracks].fTOFhit) {
1533 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1534 }
1535 ///////////////////
1536
1537 // Use TOF if good hit and above threshold
1538 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1539 if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1540 if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1541 if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1542 if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1543 }else {// TPC info instead
1544 if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1545 if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1546 if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1547 if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1548 }
1549
fa109294 1550 //cout<<nSigmaTPC[2]<<endl;
cd12341d 1551 //////////////////////////////////////
1552 // Bayesian PIDs for TPC only tracking
1553 //const Double_t* PID = aodtrack->PID();
1554 //fTempStruct[myTracks].fElectron = kFALSE;
1555 //fTempStruct[myTracks].Pion = kFALSE;
1556 //fTempStruct[myTracks].Kaon = kFALSE;
1557 //fTempStruct[myTracks].Proton = kFALSE;
1558
1559 // Pions
1560 //if(PID[2] > 0.2) fTempStruct[myTracks].Pion = kTRUE;// pions: 0.2 --> 0.5
1561 //
1562 // Kaons
1563 //if(PID[3] <= 0.5) fTempStruct[myTracks].Kaon = kFALSE;// kaons
1564 //
1565 // Protons
1566 //if(PID[4] <= 0.5) fTempStruct[myTracks].Proton = kFALSE;// protons
1567 //////////////////////////////////////
1568
1569
1570 // Ensure there is only 1 candidate per track
1571 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1572 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1573 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1574 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1575 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1576 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1577 ////////////////////////
1578 if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
1579
1580 if(!fTempStruct[myTracks].fPion) continue;// only pions
1581
1582
1583
1584
1585 if(fTempStruct[myTracks].fPion) {// pions
1586 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1587 fTempStruct[myTracks].fKey = 1;
1588 }else if(fTempStruct[myTracks].fKaon){// kaons
1589 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1590 fTempStruct[myTracks].fKey = 10;
1591 }else{// protons
1592 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1593 fTempStruct[myTracks].fKey = 100;
1594 }
1595
1596
1597
1598 if(aodtrack->Charge() > 0) positiveTracks++;
1599 else negativeTracks++;
1600
1601 if(fTempStruct[myTracks].fPion) pionCount++;
1602 if(fTempStruct[myTracks].fKaon) kaonCount++;
1603 if(fTempStruct[myTracks].fProton) protonCount++;
1604
1605 myTracks++;
1606
1607 }
1608 }else {// ESD tracks
1609 cout<<"ESDs not supported currently"<<endl;
1610 return;
1611 }
1612
1613
1614 if(myTracks >= 1) {
1615 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1616 }
1617
1618
b6e5ec54 1619 //cout<<"There are "<<myTracks<<" myTracks"<<endl;
1620 //cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
cd12341d 1621
1622 /////////////////////////////////////////
1623 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1624 if(myTracks < 3) {cout<<"Less than 3 tracks. Skipping Event."<<endl; return;}
1625 /////////////////////////////////////////
fa109294 1626
cd12341d 1627
1628 ////////////////////////////////
1629 ///////////////////////////////
1630 // Mbin determination
1631 //
1632 // Mbin set to Pion Count Only for pp!!!!!!!
1633 fMbin=-1;
1634 if(!fPbPbcase){
1635 for(Int_t i=0; i<kMultBinspp; i++){
1636 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1637 // Mbin 0 has 1 pion
1638 }
1639 }else{
1640 for(Int_t i=0; i<kCentBins; i++){
1641 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1642 fMbin=i;// 0 = most central
1643 break;
1644 }
1645 }
1646 }
1647
1648 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1649
704f2481 1650 fFSIbin=0;
ae9b34d1 1651 if(fMbin==0) fFSIbin = 0;//0-5%
1652 else if(fMbin==1) fFSIbin = 1;//5-10%
1653 else if(fMbin<=3) fFSIbin = 2;//10-20%
1654 else if(fMbin<=5) fFSIbin = 3;//20-30%
1655 else if(fMbin<=7) fFSIbin = 4;//30-40%
1656 else fFSIbin = 5;//40-50%
1657
90814457 1658 Int_t rIndexForTPN = 4;
1659 if(fMbin<=1) {rIndexForTPN=5;}
1660 else if(fMbin<=3) {rIndexForTPN=4;}
1661 else if(fMbin<=5) {rIndexForTPN=3;}
1662 else {rIndexForTPN=2;}
1663
cd12341d 1664 //////////////////////////////////////////////////
1665 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1666 //////////////////////////////////////////////////
1667
1668
cd12341d 1669 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1670 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1671
1672 ////////////////////////////////////
1673 // Add event to buffer if > 0 tracks
1674 if(myTracks > 0){
1675 fEC[zbin][fMbin]->FIFOShift();
1676 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1677 (fEvt)->fNtracks = myTracks;
1678 (fEvt)->fFillStatus = 1;
1679 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1680 if(fMCcase){
1681 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1682 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1683 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1684 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1685 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1686 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1687 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1688 }
1689 }
1690 }
1691
1692
1693
1694 Float_t qinv12=0, qinv13=0, qinv23=0;
1695 Float_t qout=0, qside=0, qlong=0;
1696 Float_t qoutMC=0, qsideMC=0, qlongMC=0;
1697 Float_t transK12=0, rapK12=0, transK3=0;
1698 Int_t transKbin=0, rapKbin=0;
1699 Float_t q3=0;
1700 Int_t ch1=0, ch2=0, ch3=0;
1701 Short_t key1=0, key2=0, key3=0;
1702 Int_t bin1=0, bin2=0, bin3=0;
1703 Float_t pVect1[4]={0};
1704 Float_t pVect2[4]={0};
1705 Float_t pVect3[4]={0};
1706 Float_t pVect1MC[4]={0};
1707 Float_t pVect2MC[4]={0};
5e3e77d6 1708 Float_t pVect3MC[4]={0};
cd12341d 1709 Int_t index1=0, index2=0, index3=0;
1710 Float_t weight12=0, weight13=0, weight23=0;
1711 Float_t weight12Err=0, weight13Err=0, weight23Err=0;
1712 Float_t weight12CC=0, weight13CC=0, weight23CC=0;
1713 Float_t weightTotal=0;//, weightTotalErr=0;
5e3e77d6 1714 Float_t qinv12MC=0, qinv13MC=0, qinv23MC=0;
90814457 1715 Float_t Qsum1v1=0, Qsum2=0, Qsum3v1=0, Qsum1v2=0, Qsum3v2=0;
1716 Float_t Qsum1v1MC=0, Qsum2MC=0, Qsum3v1MC=0, Qsum1v2MC=0, Qsum3v2MC=0;
1717 //
cd12341d 1718 AliAODMCParticle *mcParticle1=0x0;
1719 AliAODMCParticle *mcParticle2=0x0;
1720
1721
1722 if(fPdensityPairCut){
1723 ////////////////////
1724 Int_t pairCountSE=0, pairCountME=0;
1725 Int_t normPairCount[2]={0};
1726 Int_t numOtherPairs1[2][fMultLimit];
1727 Int_t numOtherPairs2[2][fMultLimit];
1728 Bool_t exitCode=kFALSE;
1729 Int_t tempNormFillCount[2][2][2][10][5];
1730
1731
1732 // reset to defaults
1733 for(Int_t i=0; i<fMultLimit; i++) {
1734 fPairLocationSE[i]->Set(fMultLimit,fDefaultsInt);
1735 fPairLocationME[i]->Set(fMultLimit,fDefaultsInt);
1736
1737 // Normalization Utilities
1738 fOtherPairLocation1[0][i]->Set(fMultLimit,fDefaultsInt);
1739 fOtherPairLocation1[1][i]->Set(fMultLimit,fDefaultsInt);
1740 fOtherPairLocation2[0][i]->Set(fMultLimit,fDefaultsInt);
1741 fOtherPairLocation2[1][i]->Set(fMultLimit,fDefaultsInt);
1742 fNormPairSwitch[0][i]->Set(fMultLimit,fDefaultsCharMult);
1743 fNormPairSwitch[1][i]->Set(fMultLimit,fDefaultsCharMult);
1744 fNormPairSwitch[2][i]->Set(fMultLimit,fDefaultsCharMult);
1745 numOtherPairs1[0][i]=0;
1746 numOtherPairs1[1][i]=0;
1747 numOtherPairs2[0][i]=0;
1748 numOtherPairs2[1][i]=0;
1749
1750 // Track Merging/Splitting Utilities
1751 fPairSplitCut[0][i]->Set(fMultLimit,fDefaultsCharMult);// P11
1752 fPairSplitCut[1][i]->Set(fMultLimit,fDefaultsCharMult);// P12
1753 fPairSplitCut[2][i]->Set(fMultLimit,fDefaultsCharMult);// P13
1754 fPairSplitCut[3][i]->Set(fMultLimit,fDefaultsCharMult);// P23
1755 }
1756
1757 // Reset the temp Normalization counters
1758 for(Int_t i=0; i<2; i++){// Charge1
1759 for(Int_t j=0; j<2; j++){// Charge2
1760 for(Int_t k=0; k<2; k++){// Charge3
1761 for(Int_t l=0; l<10; l++){// FillIndex (species Combination)
1762 for(Int_t m=0; m<5; m++){// Term (Cumulant term)
1763 tempNormFillCount[i][j][k][l][m] = 0;
1764 }
1765 }
1766 }
1767 }
1768 }
1769
1770
1771 ///////////////////////////////////////////////////////
1772 // Start the pairing process
1773 // P11 pairing
1774 // 1st Particle
1775
1776 for (Int_t i=0; i<myTracks; i++) {
1777
1778 Int_t en2=0;
1779
1780 // 2nd particle
1781 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {
1782
1783 key1 = (fEvt)->fTracks[i].fKey;
1784 key2 = (fEvt+en2)->fTracks[j].fKey;
1785 Short_t fillIndex2 = FillIndex2part(key1+key2);
1786 Short_t qCutBin = SetQcutBin(fillIndex2);
1787 Short_t normBin = SetNormBin(fillIndex2);
1788 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1789 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1790 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1791 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1792 //
1793
1794 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
1795 GetQosl(pVect1, pVect2, qout, qside, qlong);
1796 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1797
1798 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1799
1800 ///////////////////////////////
1801 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
1802 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1803 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
1804
1805 if(fMCcase && ch1==ch2 && fMbin==0){
1806 for(Int_t rstep=0; rstep<10; rstep++){
1807 Float_t coeff = (rstep)*0.2*(0.18/1.2);
1808 // propagate through B field to r=1.2m
1809 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
1810 if(phi1 > 2*PI) phi1 -= 2*PI;
1811 if(phi1 < 0) phi1 += 2*PI;
1812 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
1813 if(phi2 > 2*PI) phi2 -= 2*PI;
1814 if(phi2 < 0) phi2 += 2*PI;
1815 Float_t deltaphi = phi1 - phi2;
1816 if(deltaphi > PI) deltaphi -= PI;
1817 if(deltaphi < -PI) deltaphi += PI;
1818 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
1819 }
1820 }
5e3e77d6 1821
cd12341d 1822 // Pair Splitting/Merging cut
1823 if(ch1 == ch2){
1824 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1825 fPairSplitCut[0][i]->AddAt('1',j);
1826 ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1827 continue;
1828 }
1829 }
5e3e77d6 1830
cd12341d 1831 // HIJING tests
1832 if(fMCcase && fillIndex2==0){
1833
1834 // Check that label does not exceed stack size
1835 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
1836 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1837 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
1838 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1839 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1840 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1841 qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
1842 GetQosl(pVect1MC, pVect2MC, qoutMC, qsideMC, qlongMC);
1843 if(qinv12<0.1 && ch1==ch2) {
1844 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
1845 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
1846 ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
1847 }
1848
1849
5e3e77d6 1850 for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 8fm + 1 Therminator setting
cd12341d 1851 for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
1852 Int_t denIndex = rIter*kNDampValues + myDampIt;
cd12341d 1853 Float_t WInput = MCWeight(ch1,ch2, rIter, myDampIt, qinv12MC);
1854 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
1855 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
1856 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
1857 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
1858 }
1859 }
1860
fa109294 1861 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
1862 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
1863
cd12341d 1864 //HIJING resonance test
1865 if(ch1 != ch2){
cd12341d 1866 ((TH1F*)fOutputList->FindObject("fAllOSPairs"))->Fill(fMbin+1, qinv12);
fa109294 1867 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
cd12341d 1868 if(mcParticle1->GetMother() == mcParticle2->GetMother() && mcParticle1->GetMother() >=0){
1869 ((TH1F*)fOutputList->FindObject("fResonanceOSPairs"))->Fill(fMbin+1, qinv12);
1870 }
1871 }
1872 }
90814457 1873 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinv->Fill(qinv12MC, MCWeight(ch1,ch2,4,5,qinv12MC));
1874 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(ch1,ch2,4,5,qinv12MC));
54d66278 1875 // pion purity
1876 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityDen->Fill(transK12, qinv12);
1877 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
1878 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityNum->Fill(transK12, qinv12);
1879 }
1880
cd12341d 1881 }// label check 2
1882 }// MC case
1883
1884 //////////////////////////////////////////
1885 // 2-particle term
1886 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
5e3e77d6 1887 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
704f2481 1888 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt)->fTracks[i].fMom);
1889 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt+en2)->fTracks[j].fMom);
fa109294 1890
cd12341d 1891 // osl frame
fa109294 1892 if(fillIndex2==0){
cd12341d 1893 if((transK12 > 0.2) && (transK12 < 0.3)){
1894 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1895 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1896 }
1897 if((transK12 > 0.6) && (transK12 < 0.7)){
1898 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1899 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1900 }
1901 }
5e3e77d6 1902
cd12341d 1903 //////////////////////////////////////////
1904 if(fTabulatePairs){
1905 if(fillIndex2==0 && bin1==bin2){
1906 rapK12 = 0;
1907 transKbin=-1; rapKbin=-1;
1908 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
1909 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
1910 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1911 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1912 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1913 continue;
1914 }
1915 }
1916
1917
1918 // exit out of loop if there are too many pairs
1919 if(pairCountSE >= kPairLimit) {cout<<"Too many SE pairs"<<endl; exitCode=kTRUE; continue;}
1920 if(exitCode) continue;
5e3e77d6 1921
cd12341d 1922 //////////////////////////
1923 // Enforce the Qcut
1924 if(qinv12 <= fQcut[qCutBin]) {
5e3e77d6 1925
cd12341d 1926 ///////////////////////////
1927 // particle 1
1928 (fEvt)->fPairsSE[pairCountSE].fP1[0] = (fEvt)->fTracks[i].fP[0];
1929 (fEvt)->fPairsSE[pairCountSE].fP1[1] = (fEvt)->fTracks[i].fP[1];
1930 (fEvt)->fPairsSE[pairCountSE].fP1[2] = (fEvt)->fTracks[i].fP[2];
1931 (fEvt)->fPairsSE[pairCountSE].fE1 = (fEvt)->fTracks[i].fEaccepted;
1932 (fEvt)->fPairsSE[pairCountSE].fCharge1 = (fEvt)->fTracks[i].fCharge;
1933 (fEvt)->fPairsSE[pairCountSE].fIndex1 = i;
1934 (fEvt)->fPairsSE[pairCountSE].fKey1 = key1;
1935 (fEvt)->fPairsSE[pairCountSE].fLabel1 = (fEvt)->fTracks[i].fLabel;
1936 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
1937 (fEvt)->fPairsSE[pairCountSE].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
1938 (fEvt)->fPairsSE[pairCountSE].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
1939 (fEvt)->fPairsSE[pairCountSE].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
5e3e77d6 1940 }
cd12341d 1941 // particle 2
1942 (fEvt)->fPairsSE[pairCountSE].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
1943 (fEvt)->fPairsSE[pairCountSE].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
1944 (fEvt)->fPairsSE[pairCountSE].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
1945 (fEvt)->fPairsSE[pairCountSE].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
1946 (fEvt)->fPairsSE[pairCountSE].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1947 (fEvt)->fPairsSE[pairCountSE].fIndex2 = j;
1948 (fEvt)->fPairsSE[pairCountSE].fKey2 = key2;
1949 (fEvt)->fPairsSE[pairCountSE].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
1950 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
1951 (fEvt)->fPairsSE[pairCountSE].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
1952 (fEvt)->fPairsSE[pairCountSE].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
1953 (fEvt)->fPairsSE[pairCountSE].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
1954 }
5e3e77d6 1955
cd12341d 1956 (fEvt)->fPairsSE[pairCountSE].fQinv = qinv12;
1957
1958 fPairLocationSE[i]->AddAt(pairCountSE,j);
1959
1960 pairCountSE++;
1961
1962 }
1963
1964
1965 /////////////////////////////////////////////////////////
1966 // Normalization Region
1967
1968 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
1969 // particle 1
1970 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
1971 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
1972 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
1973 // particle 2
1974 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
1975 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
1976 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
1977
1978
1979 //other past pairs with particle j
1980 for(Int_t pastpair=0; pastpair<numOtherPairs2[0][j]; pastpair++){
1981 Int_t locationOtherPair = fOtherPairLocation2[0][j]->At(pastpair);
1982 if(locationOtherPair < 0) continue;// no pair there
1983 Int_t indexOther1 = i;
1984 Int_t indexOther2 = fNormPairs[0][ locationOtherPair ].fIndex1;
1985
1986 // Both possible orderings of other indexes
1987 if( (fNormPairSwitch[0][indexOther1]->At(indexOther2)=='1') || (fNormPairSwitch[0][indexOther2]->At(indexOther1)=='1')) {
1988
1989 // 1 and 2 are from SE
1990 ch3 = Int_t((fNormPairs[0][ locationOtherPair ].fCharge1 + 1)/2.);
1991 key3 = fNormPairs[0][ locationOtherPair ].fKey1;
1992 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
1993 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
1994
1995 tempNormFillCount[bin1][bin2][bin3][fillIndex3][0]++;
1996 }
1997
1998 }// pastpair P11 loop
1999
2000
2001 fNormPairSwitch[en2][i]->AddAt('1',j);
2002 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
2003 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
2004
2005 numOtherPairs1[en2][i]++;
2006 numOtherPairs2[en2][j]++;
2007
2008
2009 normPairCount[en2]++;
2010 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
2011
2012 }// Norm Region
2013
2014 }// j particle
2015 }// i particle
2016
2017
2018
2019 //////////////////////////////////////////////
2020 // P12 pairing
2021 // 1st Particle
2022 for (Int_t i=0; i<myTracks; i++) {
2023
2024 Int_t en2=1;
2025
2026 // 2nd particle
2027 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2028
2029 key1 = (fEvt)->fTracks[i].fKey;
2030 key2 = (fEvt+en2)->fTracks[j].fKey;
2031 Short_t fillIndex2 = FillIndex2part(key1+key2);
2032 Short_t qCutBin = SetQcutBin(fillIndex2);
2033 Short_t normBin = SetNormBin(fillIndex2);
2034 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2035 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2036 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2037 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2038
2039 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2040 GetQosl(pVect1, pVect2, qout, qside, qlong);
2041 transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2042
2043 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2044
2045 ///////////////////////////////
2046 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2047 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2048 SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
2049
2050 if(fMCcase && ch1==ch2 && fMbin==0){
2051 for(Int_t rstep=0; rstep<10; rstep++){
2052 Float_t coeff = (rstep)*0.2*(0.18/1.2);
2053 // propagate through B field to r=1.2m
2054 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2055 if(phi1 > 2*PI) phi1 -= 2*PI;
2056 if(phi1 < 0) phi1 += 2*PI;
2057 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2058 if(phi2 > 2*PI) phi2 -= 2*PI;
2059 if(phi2 < 0) phi2 += 2*PI;
2060 Float_t deltaphi = phi1 - phi2;
2061 if(deltaphi > PI) deltaphi -= PI;
2062 if(deltaphi < -PI) deltaphi += PI;
2063 ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2064 }
2065 }
2066
2067 if(ch1 == ch2){
2068 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2069 fPairSplitCut[1][i]->AddAt('1',j);
2070 continue;
2071 }
2072 }
2073
2074 //////////////////////////////////////////
2075 // 2-particle term
2076 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
5e3e77d6 2077 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
2078
cd12341d 2079 // osl frame
fa109294 2080 if(fillIndex2==0){
2081 if((transK12 > 0.2) && (transK12 < 0.3)){
cd12341d 2082 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2083 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2084 }
2085 if((transK12 > 0.6) && (transK12 < 0.7)){
2086 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2087 Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fExplicit2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2088 }
2089 }
2090 //////////////////////////////////////////
2091 if(fTabulatePairs){
2092 if(fillIndex2==0 && bin1==bin2){
2093 rapK12 = 0;
2094 transKbin=-1; rapKbin=-1;
2095 for(Int_t kIt=0; kIt<kKbinsT; kIt++) {if(transK12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {transKbin = kIt; break;}}
2096 for(Int_t kIt=0; kIt<kKbinsY; kIt++) {if(rapK12 < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {rapKbin = kIt; break;}}
2097 if((transKbin<0) || (rapKbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2098 if((transKbin>=kKbinsT) || (rapKbin>=kKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2099 KT[transKbin].KY[rapKbin].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2100 continue;
2101 }
2102 }
2103
2104
2105 if(pairCountME >= 2*kPairLimit) {cout<<"Too many ME pairs"<<endl; exitCode=kTRUE; continue;}
2106 if(exitCode) continue;
2107
2108 if(qinv12 <= fQcut[qCutBin]) {
2109 ///////////////////////////
2110
2111 // particle 1
2112 (fEvt)->fPairsME[pairCountME].fP1[0] = (fEvt)->fTracks[i].fP[0];
2113 (fEvt)->fPairsME[pairCountME].fP1[1] = (fEvt)->fTracks[i].fP[1];
2114 (fEvt)->fPairsME[pairCountME].fP1[2] = (fEvt)->fTracks[i].fP[2];
2115 (fEvt)->fPairsME[pairCountME].fE1 = (fEvt)->fTracks[i].fEaccepted;
2116 (fEvt)->fPairsME[pairCountME].fCharge1 = (fEvt)->fTracks[i].fCharge;
2117 (fEvt)->fPairsME[pairCountME].fIndex1 = i;
2118 (fEvt)->fPairsME[pairCountME].fKey1 = key1;
2119 (fEvt)->fPairsME[pairCountME].fLabel1 = (fEvt)->fTracks[i].fLabel;
2120 if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
2121 (fEvt)->fPairsME[pairCountME].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
2122 (fEvt)->fPairsME[pairCountME].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
2123 (fEvt)->fPairsME[pairCountME].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
2124 }
2125 // particle 2
2126 (fEvt)->fPairsME[pairCountME].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
2127 (fEvt)->fPairsME[pairCountME].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
2128 (fEvt)->fPairsME[pairCountME].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
2129 (fEvt)->fPairsME[pairCountME].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
2130 (fEvt)->fPairsME[pairCountME].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2131 (fEvt)->fPairsME[pairCountME].fIndex2 = j;
2132 (fEvt)->fPairsME[pairCountME].fKey2 = key2;
2133 (fEvt)->fPairsME[pairCountME].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
2134 if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
2135 (fEvt)->fPairsME[pairCountME].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2136 (fEvt)->fPairsME[pairCountME].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2137 (fEvt)->fPairsME[pairCountME].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2138 }
2139
2140 (fEvt)->fPairsME[pairCountME].fQinv = qinv12;
2141
2142 fPairLocationME[i]->AddAt(Int_t(pairCountME),j);
2143
2144 pairCountME++;
2145
2146 }
2147
2148 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2149 // particle 1
2150 fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
2151 fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
2152 fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
2153 // particle 2
2154 fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
2155 fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
2156 fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
2157
2158 //other past pairs in P11 with particle i
2159 for(Int_t pastpairP11=0; pastpairP11<numOtherPairs2[0][i]; pastpairP11++){// past pair in P11 with i as 1st and 2nd particle
2160 Int_t locationOtherPairP11 = fOtherPairLocation2[0][i]->At(pastpairP11);// i is 2nd particle
2161 if(locationOtherPairP11 < 0) continue;// no pair there
2162 Int_t indexOther1P11 = fNormPairs[0][ locationOtherPairP11 ].fIndex1;
2163
2164 //Check other past pairs in P12
2165 if( (fNormPairSwitch[1][indexOther1P11]->At(j)=='0')) continue;
2166
2167 // 1 and 3 are from SE
2168 ch3 = Int_t((fNormPairs[0][ locationOtherPairP11 ].fCharge1 + 1)/2.);// charge of second particle in P11
2169 key3 = fNormPairs[0][ locationOtherPairP11 ].fKey1;
2170 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2171 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2172 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 2, bin1, bin2, bin3, fill2, fill3, fill4);
2173
2174
2175 if(fill2) tempNormFillCount[bin1][bin2][bin3][fillIndex3][1]++;
2176 if(fill3) tempNormFillCount[bin1][bin2][bin3][fillIndex3][2]++;
2177 if(fill4) tempNormFillCount[bin1][bin2][bin3][fillIndex3][3]++;
2178
2179
2180 }// P11 loop
2181
2182
2183 fNormPairSwitch[en2][i]->AddAt('1',j);
2184 fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
2185 fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
2186
2187 numOtherPairs1[en2][i]++;
2188 numOtherPairs2[en2][j]++;
2189
2190 normPairCount[en2]++;
2191 if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
2192
2193 }// Norm Region
2194
2195
2196 }
2197 }
2198
2199
2200 ///////////////////////////////////////
2201 // P13 pairing (just for Norm counting of term5)
2202 for (Int_t i=0; i<myTracks; i++) {
2203
2204 // exit out of loop if there are too many pairs
2205 // dont bother with this loop if exitCode is set.
2206 if(exitCode) break;
2207
2208 // 2nd particle
2209 Int_t en2=2;
2210
2211 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2212
2213 key1 = (fEvt)->fTracks[i].fKey;
2214 key2 = (fEvt+en2)->fTracks[j].fKey;
2215 Short_t fillIndex2 = FillIndex2part(key1+key2);
2216 Short_t normBin = SetNormBin(fillIndex2);
2217 pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2218 pVect1[1]=(fEvt)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2219 pVect1[2]=(fEvt)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2220 pVect1[3]=(fEvt)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2221
2222 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2223
2224 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2225
2226 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2227 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2228
2229 if(ch1 == ch2){
2230 if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2231 fPairSplitCut[2][i]->AddAt('1',j);
2232 continue;
2233 }
2234 }
2235
2236 /////////////////////////////////////////////////////////
2237 // Normalization Region
2238
2239 if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
2240
2241 fNormPairSwitch[en2][i]->AddAt('1',j);
2242
2243 }// Norm Region
2244 }
2245 }
2246
2247
2248
2249 ///////////////////////////////////////
2250 // P23 pairing (just for Norm counting of term5)
2251 Int_t en1=1;
2252 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {
2253
2254 // exit out of loop if there are too many pairs
2255 // dont bother with this loop if exitCode is set.
2256 if(exitCode) break;
2257
2258 // 2nd event
2259 Int_t en2=2;
2260 // 2nd particle
2261 for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
2262
2263 if(exitCode) break;
2264
2265 key1 = (fEvt+en1)->fTracks[i].fKey;
2266 key2 = (fEvt+en2)->fTracks[j].fKey;
2267 Short_t fillIndex2 = FillIndex2part(key1+key2);
2268 Short_t normBin = SetNormBin(fillIndex2);
2269 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2270 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2271 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2272 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2273
2274 qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
2275
2276 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
2277
2278 ///////////////////////////////
2279 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2280 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2281
2282 if(ch1 == ch2){
2283 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2284 fPairSplitCut[3][i]->AddAt('1',j);
2285 continue;
2286 }
2287 }
2288
2289 if((qinv12 < fNormQcutLow[normBin]) || (qinv12 >= fNormQcutHigh[normBin])) continue;
2290
2291 Int_t index1P23 = i;
2292 Int_t index2P23 = j;
2293
2294 for(Int_t pastpairP12=0; pastpairP12<numOtherPairs2[1][index1P23]; pastpairP12++){// loop in P12 with i as 2nd particle
2295 Int_t locationOtherPairP12 = fOtherPairLocation2[1][index1P23]->At(pastpairP12);
2296 if(locationOtherPairP12 < 0) continue; // no pair there
2297 Int_t index1P12 = fNormPairs[1][ locationOtherPairP12 ].fIndex1;
2298
2299
2300 //Check other past pair status in P13
2301 if( (fNormPairSwitch[2][index1P12]->At(index2P23)=='0')) continue;
2302
2303 // all from different event
2304 ch3 = Int_t((fNormPairs[1][ locationOtherPairP12 ].fCharge1 + 1)/2.);// charge of first particle in P12
2305 key3 = fNormPairs[1][ locationOtherPairP12 ].fKey1;
2306 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2307 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2308
2309 tempNormFillCount[bin1][bin2][bin3][fillIndex3][4]++;
2310 }
2311 }
2312 }
2313
2314
2315
2316
2317 ///////////////////////////////////////////////////
2318 // Do not use pairs from events with too many pairs
2319 if(exitCode) {
2320 cout<<"SE or ME or Norm PairCount too large. Discarding all pairs and skipping event"<<endl;
2321 (fEvt)->fNpairsSE = 0;
2322 (fEvt)->fNpairsME = 0;
2323 ((TH1F*)fOutputList->FindObject("fRejectedEvents"))->Fill(fMbin+1);
2324 return;// Skip event
2325 }else{
2326 (fEvt)->fNpairsSE = pairCountSE;
2327 (fEvt)->fNpairsME = pairCountME;
2328 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2329 }
2330 ///////////////////////////////////////////////////
2331
cd12341d 2332
2333 //cout<<"pairCountSE = "<<pairCountSE<<" pairCountME = "<<pairCountME<<endl;
b6e5ec54 2334 //cout<<"Start Main analysis"<<endl;
cd12341d 2335
2336 ///////////////////////////////////////////////////////////////////////
2337 ///////////////////////////////////////////////////////////////////////
2338 ///////////////////////////////////////////////////////////////////////
2339 //
2340 //
5e3e77d6 2341 // Start the Main Correlation Analysis
cd12341d 2342 //
2343 //
2344 ///////////////////////////////////////////////////////////////////////
2345
2346
2347 /////////////////////////////////////////////////////////
2348 // Skip 3-particle part if Tabulate6DPairs is set to true
2349 if(fTabulatePairs) return;
2350 /////////////////////////////////////////////////////////
2351
2352 // Set the Normalization counters
2353 for(Int_t termN=0; termN<5; termN++){
2354
2355 if(termN==0){
2356 if((fEvt)->fNtracks ==0) continue;
2357 }else if(termN<4){
2358 if((fEvt)->fNtracks ==0) continue;
2359 if((fEvt+1)->fNtracks ==0) continue;
2360 }else {
2361 if((fEvt)->fNtracks ==0) continue;
2362 if((fEvt+1)->fNtracks ==0) continue;
2363 if((fEvt+2)->fNtracks ==0) continue;
2364 }
2365
2366 for(Int_t sc=0; sc<kSCLimit3; sc++){
2367
2368 for(Int_t c1=0; c1<2; c1++){
2369 for(Int_t c2=0; c2<2; c2++){
2370 for(Int_t c3=0; c3<2; c3++){
2371
2372 if(sc==0 || sc==6 || sc==9){// Identical species
2373 if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
2374 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
2375 }else if(sc!=5){
2376 if( (c1+c2)==1) {if(c1!=0) continue;}
2377 }else {}// do nothing for pi-k-p case
2378 Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[fMbin].EDB[fEDbin].ThreePT[termN].fNorm3->Fill(0.,tempNormFillCount[c1][c2][c3][sc][termN]);
2379 }
2380 }
2381 }
2382 }
2383 }
2384
2385
2386
2387 /////////////////////////////////////////////
2388 // Calculate Pair-Cut Correlations
2389 for(Int_t en1case=0; en1case<2; en1case++){// limit at 2 (normal)
2390
2391 Int_t nump1=0;
2392 if(en1case==0) nump1 = (fEvt)->fNpairsSE;
2393 if(en1case==1) nump1 = (fEvt)->fNpairsME;
2394
2395 // 1st pair
2396 for(Int_t p1=0; p1<nump1; p1++){
2397
2398 if(en1case==0){
2399 ch1 = Int_t(((fEvt)->fPairsSE[p1].fCharge1 + 1)/2.);
2400 ch2 = Int_t(((fEvt)->fPairsSE[p1].fCharge2 + 1)/2.);
2401 pVect1[0] = (fEvt)->fPairsSE[p1].fE1; pVect2[0] = (fEvt)->fPairsSE[p1].fE2;
2402 pVect1[1] = (fEvt)->fPairsSE[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsSE[p1].fP2[0];
2403 pVect1[2] = (fEvt)->fPairsSE[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsSE[p1].fP2[1];
2404 pVect1[3] = (fEvt)->fPairsSE[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsSE[p1].fP2[2];
2405 index1 = (fEvt)->fPairsSE[p1].fIndex1; index2 = (fEvt)->fPairsSE[p1].fIndex2;
2406 key1 = (fEvt)->fPairsSE[p1].fKey1; key2 = (fEvt)->fPairsSE[p1].fKey2;
2407 pVect1MC[1] = (fEvt)->fPairsSE[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsSE[p1].fP2MC[0];
2408 pVect1MC[2] = (fEvt)->fPairsSE[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsSE[p1].fP2MC[1];
2409 pVect1MC[3] = (fEvt)->fPairsSE[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsSE[p1].fP2MC[2];
2410 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2411 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2412
2413 qinv12 = (fEvt)->fPairsSE[p1].fQinv;
2414 }
2415 if(en1case==1){
2416 ch1 = Int_t(((fEvt)->fPairsME[p1].fCharge1 + 1)/2.);
2417 ch2 = Int_t(((fEvt)->fPairsME[p1].fCharge2 + 1)/2.);
2418 pVect1[0] = (fEvt)->fPairsME[p1].fE1; pVect2[0] = (fEvt)->fPairsME[p1].fE2;
2419 pVect1[1] = (fEvt)->fPairsME[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsME[p1].fP2[0];
2420 pVect1[2] = (fEvt)->fPairsME[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsME[p1].fP2[1];
2421 pVect1[3] = (fEvt)->fPairsME[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsME[p1].fP2[2];
2422 index1 = (fEvt)->fPairsME[p1].fIndex1; index2 = (fEvt)->fPairsME[p1].fIndex2;
2423 key1 = (fEvt)->fPairsME[p1].fKey1; key2 = (fEvt)->fPairsME[p1].fKey2;
2424 pVect1MC[1] = (fEvt)->fPairsME[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsME[p1].fP2MC[0];
2425 pVect1MC[2] = (fEvt)->fPairsME[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsME[p1].fP2MC[1];
2426 pVect1MC[3] = (fEvt)->fPairsME[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsME[p1].fP2MC[2];
2427 pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
2428 pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
2429
2430 qinv12 = (fEvt)->fPairsME[p1].fQinv;
2431 }
2432
cd12341d 2433
2434 // en2 buffer
2435 for(Int_t en2=0; en2<3; en2++){
2436 //////////////////////////////////////
2437
2438 Bool_t skipcase=kTRUE;
2439 Short_t config=-1, part=-1;
2440 if(en1case==0 && en2==0) {skipcase=kFALSE; config=1; part=0;}// P11T1
2441 if(en1case==0 && en2==1) {skipcase=kFALSE; config=2; part=1;}// P11T2
2442 if(en1case==1 && en2==0) {skipcase=kFALSE; config=2; part=2;}// P12T1
2443 if(en1case==1 && en2==2) {skipcase=kFALSE; config=3; part=3;}// P12T3
2444
2445 if(skipcase) continue;
2446
2447
2448 // 3-particle terms
2449 // 3rd particle
2450 for(Int_t k=0; k<(fEvt+en2)->fNtracks; k++){
2451 index3 = k;
2452
2453
2454 // remove auto-correlations and duplicate triplets
2455 if(config==1){
2456 if( index1 == index3) continue;
2457 if( index2 == index3) continue;
2458 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2459
2460 // skip the switched off triplets
2461 if(fTripletSkip1[fPairLocationSE[index1]->At(index2)]->At(index3)=='1') {
2462 fTripletSkip1[fPairLocationSE[index1]->At(index2)]->AddAt('0',index3);// Reset
2463 continue;
2464 }
2465 ///////////////////////////////
2466 // Turn off 1st possible degenerate triplet
2467 if(index1 < index3){// verify correct id ordering ( index1 < k )
2468 if(fPairLocationSE[index1]->At(index3) >= 0){
2469 fTripletSkip1[fPairLocationSE[index1]->At(index3)]->AddAt('1',index2);
2470 }
2471 if(fPairSplitCut[0][index1]->At(index3)=='1') continue;// Track splitting/merging
2472 }else {// or k < index1
2473 if(fPairLocationSE[index3]->At(index1) >= 0){
2474 fTripletSkip1[fPairLocationSE[index3]->At(index1)]->AddAt('1',index2);
2475 }
2476 if(fPairSplitCut[0][index3]->At(index1)=='1') continue;// Track splitting/merging
2477 }
2478 // turn off 2nd possible degenerate triplet
2479 if(index2 < index3){// verify correct id ordering (index2 < k)
2480 if(fPairLocationSE[index2]->At(index3) >= 0){
2481 fTripletSkip1[fPairLocationSE[index2]->At(index3)]->AddAt('1',index1);
2482 }
2483 if(fPairSplitCut[0][index2]->At(index3)=='1') continue;// Track splitting/merging
2484 }else {// or k < index2
2485 if(fPairLocationSE[index3]->At(index2) >= 0){
2486 fTripletSkip1[fPairLocationSE[index3]->At(index2)]->AddAt('1',index1);
2487 }
2488 if(fPairSplitCut[0][index3]->At(index2)=='1') continue;// Track splitting/merging
2489 }
2490
2491 }// end config 1
2492
2493 if(config==2 && part==1){// SE pair and third particle from next event. P11T2
2494 ///////////////////////////////
2495 // Turn off 1st possible degenerate triplet
2496 if(fPairLocationME[index1]->At(index3) >= 0){
2497 fTripletSkip2[fPairLocationME[index1]->At(index3)]->AddAt('1',index2);
2498 }
2499
2500 // turn off 2nd possible degenerate triplet
2501 if(fPairLocationME[index2]->At(index3) >= 0){
2502 fTripletSkip2[fPairLocationME[index2]->At(index3)]->AddAt('1',index1);
2503 }
2504
2505 if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
2506 if(fPairSplitCut[1][index1]->At(index3)=='1') continue;// Track splitting/merging
2507 if(fPairSplitCut[1][index2]->At(index3)=='1') continue;// Track splitting/merging
2508 }// end config 2 part 1
2509
2510 if(config==2 && part==2){// P12T1
2511 if( index1 == index3) continue;
2512
2513 // skip the switched off triplets
2514 if(fTripletSkip2[fPairLocationME[index1]->At(index2)]->At(index3)=='1') {
2515 fTripletSkip2[fPairLocationME[index1]->At(index2)]->AddAt('0',index3);// Reset
2516 continue;
2517 }
2518 // turn off another possible degenerate
2519 if(fPairLocationME[index3]->At(index2) >= 0){
2520 fTripletSkip2[fPairLocationME[index3]->At(index2)]->AddAt('1',index1);
2521 }// end config 2 part 2
2522
2523 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2524 if(index1 < index3) {if(fPairSplitCut[0][index1]->At(index3)=='1') continue;}// Track splitting/merging
2525 else {if(fPairSplitCut[0][index3]->At(index1)=='1') continue;}// Track splitting/merging
2526 if(fPairSplitCut[1][index3]->At(index2)=='1') continue;// Track splitting/merging
2527 }
2528 if(config==3){// P12T3
2529 if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
2530 if(fPairSplitCut[2][index1]->At(index3)=='1') continue;// Track splitting/merging
2531 if(fPairSplitCut[3][index2]->At(index3)=='1') continue;// Track splitting/merging
2532 }// end config 3
2533
2534
5e3e77d6 2535
cd12341d 2536 ch3 = Int_t(((fEvt+en2)->fTracks[k].fCharge + 1)/2.);
2537 key3 = (fEvt+en2)->fTracks[k].fKey;
2538 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2539 Short_t fillIndex13 = FillIndex2part(key1+key3);
2540 Short_t fillIndex23 = FillIndex2part(key2+key3);
2541 Short_t qCutBin13 = SetQcutBin(fillIndex13);
2542 Short_t qCutBin23 = SetQcutBin(fillIndex23);
2543 pVect3[0] = (fEvt+en2)->fTracks[k].fEaccepted;
2544 pVect3[1] = (fEvt+en2)->fTracks[k].fP[0];
2545 pVect3[2] = (fEvt+en2)->fTracks[k].fP[1];
2546 pVect3[3] = (fEvt+en2)->fTracks[k].fP[2];
5e3e77d6 2547 if(fMCcase){
2548 pVect3MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPx;
2549 pVect3MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPy;
2550 pVect3MC[3] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPz;
2551 pVect3MC[0] = sqrt(pow(pVect3MC[1],2)+pow(pVect3MC[2],2)+pow(pVect3MC[3],2)+pow(fTrueMassPi,2));
2552 qinv12MC = GetQinv(0, pVect1MC, pVect2MC);
2553 qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
2554 qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
2555 }
cd12341d 2556 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2557 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2558
5e3e77d6 2559
cd12341d 2560 if(qinv13 < fQLowerCut) continue;
2561 if(qinv23 < fQLowerCut) continue;
d60cf3cf 2562 if(qinv13 > fQcut[qCutBin13]) continue;
2563 if(qinv23 > fQcut[qCutBin23]) continue;
5e3e77d6 2564 // if all three pair cuts are the same then the case (config=2 && term=2) never reaches here.
2565
cd12341d 2566 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2567 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2568 Float_t firstQ=0, secondQ=0, thirdQ=0;
5e3e77d6 2569 //
fa109294 2570
5e3e77d6 2571 //
cd12341d 2572 if(config==1) {// 123
2573 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
fa109294 2574
cd12341d 2575 if(fillIndex3 <= 2){
2576 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 0, 1, firstQ, secondQ, thirdQ);
2577 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ);
2578 ((TH2F*)fOutputList->FindObject("fKt3Dist"))->Fill(fMbin+1, transK3);
5e3e77d6 2579 //
2580 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
90814457 2581 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2582 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2583 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
5e3e77d6 2584 }
2585 //////////////////////////////////////
90814457 2586 // Momentum resolution and <K3> calculation
5e3e77d6 2587 if(fillIndex3==0 && fMCcase){
2588 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2589 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, 0, 1, firstQMC, secondQMC, thirdQMC);
54d66278 2590 Float_t WInput = 1;
90814457 2591 Double_t K3=1.0;
5e3e77d6 2592 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2593 WInput = MCWeight3D(kTRUE, 1, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2594 K3 = FSICorrelationOmega0(kTRUE, firstQMC, secondQMC, thirdQMC);// K3
5e3e77d6 2595 }else {// mixed charge
034e970d 2596 if(bin1==bin2) {
2597 WInput = MCWeight3D(kFALSE, 1, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2598 K3 = FSICorrelationOmega0(kFALSE, firstQMC, secondQMC, thirdQMC);// K3
2599 }else {
2600 WInput = MCWeight3D(kFALSE, 1, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2601 K3 = FSICorrelationOmega0(kFALSE, thirdQMC, secondQMC, firstQMC);// K3
2602 }
5e3e77d6 2603 }
704f2481 2604
54d66278 2605 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2606 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2607 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fQW12->Fill(firstQMC, secondQMC, thirdQMC, WInput*firstQMC);
2608 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fQW13->Fill(firstQMC, secondQMC, thirdQMC, WInput*secondQMC);
90814457 2609 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fSumK3->Fill(firstQMC, secondQMC, thirdQMC, WInput/K3);
2610 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fEnK3->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2611 if(ch1==ch2 && ch1==ch3){
2612 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2613 FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2614 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2615 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2616 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2617 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2618 //
d60cf3cf 2619 if(qinv12MC > fQLowerCut && qinv13MC > fQLowerCut && qinv23MC > fQLowerCut){
2620 // does not really matter if MC or real data triplets are used to average 1/K3...but better to use umsmeared values
2621 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsSumK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K3);
2622 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsSumK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K3);
2623 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd1TermsEnK3->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2624 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].f4VectProd2TermsEnK3->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2625 }
90814457 2626 }
fa109294 2627 }// fMCcase
2628
cd12341d 2629 }
2630
2631 }else if(config==2){// 12, 13, 23
5e3e77d6 2632
cd12341d 2633 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2634 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
2635
2636 // loop over terms 2-4
5e3e77d6 2637 for(Int_t jj=2; jj<5; jj++){
2638 if(jj==2) {if(!fill2) continue;}//12
2639 if(jj==3) {if(!fill3) continue;}//13
2640 if(jj==4) {if(!fill4) continue;}//23
cd12341d 2641
2642 if(fillIndex3 <= 2){
5e3e77d6 2643 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj, firstQ, secondQ, thirdQ);
2644 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTerms3->Fill(firstQ, secondQ, thirdQ);
90814457 2645 if(fillIndex3==0 && ch1==ch2 && ch1==ch3){
2646 if(part==1){// P11T2
2647 if(jj==2) {
2648 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2649 if(fMCcase) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2650 }if(jj==3){
2651 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2652 if(fMCcase) FourVectProdTerms(pVect1MC, pVect3MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2653 }if(jj==4) {
2654 FourVectProdTerms(pVect3, pVect1, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2655 if(fMCcase) FourVectProdTerms(pVect3MC, pVect1MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2656 }
2657 }else{// P12T1
2658 if(jj==2) {
2659 FourVectProdTerms(pVect1, pVect3, pVect2, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2660 if(fMCcase) FourVectProdTerms(pVect1MC, pVect3MC, pVect2MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2661 }if(jj==3) {
2662 FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2663 if(fMCcase) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2664 }if(jj==4) {
2665 FourVectProdTerms(pVect2, pVect1, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2666 if(fMCcase) FourVectProdTerms(pVect2MC, pVect1MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);// 4-vector product sums
2667 }
2668 }
2669 if(!fMCcase){
2670 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2671 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
fa109294 2672 }
5e3e77d6 2673 }
2674 //////////////////////////////////////
2675 // Momentum resolution calculation
2676 if(fillIndex3==0 && fMCcase){
2677 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2678 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, jj, firstQMC, secondQMC, thirdQMC);
90814457 2679 Float_t WInput = 1;
5e3e77d6 2680 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2681 WInput = MCWeight3D(kTRUE, jj, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
5e3e77d6 2682 }else {// mixed charge
90814457 2683 if(bin1==bin2) WInput = MCWeight3D(kFALSE, jj, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2684 else WInput = MCWeight3D(kFALSE, 6-jj, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
2685 }
2686 //
2687 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2688 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2689 //
2690 if(ch1==ch2 && ch1==ch3){
2691 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2692 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2693 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2694 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2695 //
d60cf3cf 2696 if(qinv12MC > fQLowerCut && qinv13MC > fQLowerCut && qinv23MC > fQLowerCut){
2697 // does not really matter if MC or real data triplets are used to average 1/K3...but better to use umsmeared values
2698 Float_t InteractingQ=0;
2699 if(part==1) {InteractingQ=qinv12;}// 12 from SE
2700 else {InteractingQ=qinv13;}// 13 from SE
2701 Double_t K2 = FSICorrelationTherm2(+1,+1, InteractingQ);// K2 from Therminator source
2702 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsSumK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput/K2);
2703 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsSumK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput/K2);
2704 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd1TermsEnK2->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2705 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].f4VectProd2TermsEnK2->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2706 }
5e3e77d6 2707 }
2708 }// fMCcase
2709
cd12341d 2710 }
2711 }
5e3e77d6 2712
cd12341d 2713 }else {// config 3: All particles from different events
cd12341d 2714
704f2481 2715 // "enhancement" differs from 1.0 only when Qinv goes over fQcut
2716 //Float_t enhancement=1.0;
2717 //Int_t nUnderCut=0;
2718 //if(qinv13<fQcut[qCutBin13]) nUnderCut++;
2719 //if(qinv23<fQcut[qCutBin23]) nUnderCut++;
2720 //if(nUnderCut==0) enhancement = (1+1+1)/1.;// 1 LowQ pair
2721 //if(nUnderCut==1) enhancement = (1+2)/2.;// 2 LowQ pair
2722 //if(nUnderCut==2) enhancement = 1.;// 3 LowQ pair
cd12341d 2723
2724 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
90814457 2725
2726 if(ch1==ch2 && ch1==ch3 && fillIndex3==0) FourVectProdTerms(pVect1, pVect2, pVect3, Qsum1v1, Qsum2, Qsum3v1, Qsum1v2, Qsum3v2);// 4-vector product sums
2727 if(fMCcase && ch1==ch2 && ch1==ch3 && fillIndex3==0) FourVectProdTerms(pVect1MC, pVect2MC, pVect3MC, Qsum1v1MC, Qsum2MC, Qsum3v1MC, Qsum1v2MC, Qsum3v2MC);
cd12341d 2728
2729 if(fillIndex3 <= 2){
2730 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
704f2481 2731 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ);
5e3e77d6 2732 if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
90814457 2733 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1Terms->Fill(Qsum1v1, Qsum2, Qsum3v1);
2734 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2Terms->Fill(Qsum1v2, Qsum2, Qsum3v2);
5e3e77d6 2735 }
2736 //////////////////////////////////////
2737 // Momentum resolution calculation
2738 if(fillIndex3==0 && fMCcase){
2739 Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
2740 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, 5, firstQMC, secondQMC, thirdQMC);
90814457 2741 Float_t WInput=1;
5e3e77d6 2742 if(ch1==ch2 && ch1==ch3){// same charge
90814457 2743 WInput = MCWeight3D(kTRUE, 5, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
5e3e77d6 2744 }else {// mixed charge
90814457 2745 if(bin1==bin2) WInput = MCWeight3D(kFALSE, 5, kMCfixedLambdabin, firstQMC, secondQMC, thirdQMC);
2746 else WInput = MCWeight3D(kFALSE, 5, kMCfixedLambdabin, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss in this case. 1st Q argument is ss
2747 }
2748 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fIdeal->Fill(firstQMC, secondQMC, thirdQMC, WInput);
2749 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fSmeared->Fill(firstQ, secondQ, thirdQ, WInput);
2750 if(ch1==ch2 && ch1==ch3){
2751 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1TermsIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, WInput);
2752 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd1TermsSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, WInput);
2753 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2TermsIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, WInput);
2754 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].f4VectProd2TermsSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, WInput);
2755
5e3e77d6 2756 }
2757 }// fMCcase
2758
cd12341d 2759 }
5e3e77d6 2760
cd12341d 2761 if(fillIndex3 !=0) continue;// only calculate TPN for pi-pi-pi
2762 if(ch1!=ch2 || ch1!=ch3) continue;// only calcualte TPN for ss
5e3e77d6 2763
90814457 2764 //if(fMCcase) continue;// only calcualte TPN for real data
2765
cd12341d 2766 GetWeight(pVect1, pVect2, weight12, weight12Err);
2767 GetWeight(pVect1, pVect3, weight13, weight13Err);
2768 GetWeight(pVect2, pVect3, weight23, weight23Err);
2769 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) continue;// weight should never be larger than 1
2770
5e3e77d6 2771
2772 // Coul correlations from Wave-functions
fa109294 2773 //for(Int_t rIter=0; rIter<kRVALUES; rIter++){// 3fm to 8fm, last value for Therminator
2774 //for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){// 0.3 to 0.6
2775 //Float_t myDamp = fDampStart + (fDampStep)*myDampIt;
2776 //Int_t denIndex = (kRVALUES-1)*kNDampValues + myDampIt;
2777 //Int_t denIndex = myDampIt;
2778 Int_t myDampIt = 5;
704f2481 2779 Float_t myDamp = 0.52;
fa109294 2780 Int_t denIndex = 0;
90814457 2781 Int_t momResIndex = rIndexForTPN*kNDampValues + myDampIt;// Therminator slot uses R=7 for momentum resolution
2782
704f2481 2783 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, qinv12);
2784 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, qinv13);
2785 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, qinv23);
d60cf3cf 2786 if(coulCorr12 < 0.1 || coulCorr13 < 0.1 || coulCorr23 < 0.1) {// Safety check
2787 if(fMbin==0 && bin1==0) {
2788 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, sqrt(fabs(weight12*weight13*weight23)));
2789 }
2790 continue;
2791 }
704f2481 2792 Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
2793 if(!fMCcase) {
2794 Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
2795 Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
2796 Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);
2797 if(momBin12 >= kQbins) momBin12 = kQbins-1;
2798 if(momBin13 >= kQbins) momBin13 = kQbins-1;
2799 if(momBin23 >= kQbins) momBin23 = kQbins-1;
2800 MomResCorr12 = fMomResC2->GetBinContent(momResIndex+1, momBin12);
2801 MomResCorr13 = fMomResC2->GetBinContent(momResIndex+1, momBin13);
2802 MomResCorr23 = fMomResC2->GetBinContent(momResIndex+1, momBin23);
d60cf3cf 2803 if(MomResCorr12 > 1.2 || MomResCorr13 > 1.2 || MomResCorr23 > 1.2) {// Safety check
2804 if(fMbin==0 && bin1==0) {
2805 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, sqrt(fabs(weight12*weight13*weight23)));
2806 }
2807 continue;
2808 }
704f2481 2809 }
2810 weight12CC = ((weight12+1)*MomResCorr12 - myDamp*coulCorr12 - (1-myDamp));
2811 weight12CC /= coulCorr12*myDamp;
2812 weight13CC = ((weight13+1)*MomResCorr13 - myDamp*coulCorr13 - (1-myDamp));
2813 weight13CC /= coulCorr13*myDamp;
2814 weight23CC = ((weight23+1)*MomResCorr23 - myDamp*coulCorr23 - (1-myDamp));
2815 weight23CC /= coulCorr23*myDamp;
2816
2817 if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {
2818 if(fMbin==0 && bin1==0) {
2819 weightTotal = sqrt(fabs(weight12CC*weight13CC*weight23CC));
2820 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, weightTotal);
2821 }
2822 continue;// C2^QS can never be less than unity
2823 }
2824
2825 /////////////////////////////////////////////////////
2826 weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
2827 /////////////////////////////////////////////////////
2828
2829 if(weightTotal > 1.5) {
2830 if(fMbin==0 && bin1==0) {
2831 ((TH3F*)fOutputList->FindObject("fTPNRejects"))->Fill(qinv12, qinv13, qinv23, weightTotal);
2832 }
2833 continue;// C2^QS never be greater than 1.0 in theory. Can be slightly larger than 1.0 with fluctuations
2834 }
2835
2836
2837
2838 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].fTwoPartNorm->Fill(qinv12, qinv13, qinv23, weightTotal);
2839
2840 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNorm->Fill(Qsum1v1, Qsum2, Qsum3v1, weightTotal);
2841 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNorm->Fill(Qsum1v2, Qsum2, Qsum3v2, weightTotal);
2842 if(fMCcase){
2843 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNormIdeal->Fill(Qsum1v1MC, Qsum2MC, Qsum3v1MC, weightTotal);
2844 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd1TwoPartNormSmeared->Fill(Qsum1v1, Qsum2, Qsum3v1, weightTotal);
2845 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNormIdeal->Fill(Qsum1v2MC, Qsum2MC, Qsum3v2MC, weightTotal);
2846 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].DT[denIndex].f4VectProd2TwoPartNormSmeared->Fill(Qsum1v2, Qsum2, Qsum3v2, weightTotal);
2847 }
90814457 2848
2849
704f2481 2850 // Save cpu time and memory by skipping r3 denominator calculation below. den errors are negligible compared to num errors.
2851 /*
2852 if(weightTotal > 0.0001){// tiny numbers cause a Float_ting point exception below
2853 weightTotalErr = pow((weight12Err*coulCorr12)*weight13CC*weight23CC,2);
2854 weightTotalErr += pow(weight12CC*(weight13Err*coulCorr13)*weight23CC,2);
2855 weightTotalErr += pow(weight12CC*weight13CC*(weight23Err*coulCorr23),2);
2856 weightTotalErr /= pow(2*weightTotal,2);
2857
2858 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].TwoPartNormErr->Fill(denIndex, q3, weightTotalErr);
2859 }
2860 */
2861
2862 //}// damp iter
2863 //}// R iter
2864
cd12341d 2865
2866
2867 }
2868 }// end 3rd particle
2869 }// en2
2870
2871
fa109294 2872 }// p1
2873 }//en1
cd12341d 2874
cd12341d 2875 ///////////////////
2876 }// end of PdensityPairs
2877
fa109294 2878
cd12341d 2879
2880
2881
2882
2883
2884 ////////////////////////////////////////////////////////
2885 // Pdensity Method with Explicit Loops
2886 if(fPdensityExplicitLoop){
2887
2888 ////////////////////////////////////
2889 // 2nd, 3rd, and 4th order Correlations
2890
2891 // First Particle
2892 for (Int_t i=0; i<myTracks; i++) {
2893 ch1 = Int_t( ((fEvt)->fTracks[i].fCharge + 1)/2. );
2894 pVect1[0] = (fEvt)->fTracks[i].fEaccepted;
2895 pVect1[1] = (fEvt)->fTracks[i].fP[0];
2896 pVect1[2] = (fEvt)->fTracks[i].fP[1];
2897 pVect1[3] = (fEvt)->fTracks[i].fP[2];
2898 key1 = (fEvt)->fTracks[i].fKey;
2899
2900 // Second Event
2901 for(Int_t en2=0; en2<fEventsToMix+1; en2++){
2902 Int_t startbin2=0;
2903 if(en2==0) startbin2=i+1;
2904
2905 // Second Particle
2906 for (Int_t j=startbin2; j<(fEvt+en2)->fNtracks; j++) {
2907 ch2 = Int_t( ((fEvt+en2)->fTracks[j].fCharge + 1)/2. );
2908 pVect2[0] = (fEvt+en2)->fTracks[j].fEaccepted;
2909 pVect2[1] = (fEvt+en2)->fTracks[j].fP[0];
2910 pVect2[2] = (fEvt+en2)->fTracks[j].fP[1];
2911 pVect2[3] = (fEvt+en2)->fTracks[j].fP[2];
2912 key2 = (fEvt+en2)->fTracks[j].fKey;
2913
2914 Short_t fillIndex12 = FillIndex2part(key1+key2);
2915 qinv12 = GetQinv(fillIndex12, pVect1, pVect2);
2916
2917 if(qinv12 < fQLowerCut) continue;
2918
2919
2920 // 2-particle part is filled always during pair creator
2921
2922 // Third Event
2923 for(Int_t en3=en2; en3<fEventsToMix+1; en3++){
2924 Int_t startbin3=0;
2925 if(en3==en2) startbin3=j+1;
2926 else startbin3=0;
2927
2928
2929 // Third Particle
2930 for (Int_t k=startbin3; k<(fEvt+en3)->fNtracks; k++) {
2931 ch3 = Int_t( ((fEvt+en3)->fTracks[k].fCharge + 1)/2. );
2932 pVect3[0] = (fEvt+en3)->fTracks[k].fEaccepted;
2933 pVect3[1] = (fEvt+en3)->fTracks[k].fP[0];
2934 pVect3[2] = (fEvt+en3)->fTracks[k].fP[1];
2935 pVect3[3] = (fEvt+en3)->fTracks[k].fP[2];
2936 key3 = (fEvt+en3)->fTracks[k].fKey;
2937
2938 Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
2939 Short_t fillIndex13 = FillIndex2part(key1+key3);
2940 qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
2941 Short_t fillIndex23 = FillIndex2part(key2+key3);
2942 qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
2943
2944
2945 if(qinv13 < fQLowerCut) continue;
2946 if(qinv23 < fQLowerCut) continue;
2947
2948
2949 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2950
2951 Short_t normBin12 = SetNormBin(fillIndex12);
2952 Short_t normBin13 = SetNormBin(fillIndex13);
2953 Short_t normBin23 = SetNormBin(fillIndex23);
2954
2955
2956 if(en3==0 && en2==0) {// 123
2957 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
2958
2959 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fExplicit3->Fill(q3);// 123
2960
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[0].fNormEx3->Fill(0.);
2964 }
2965 }
2966
2967 }else if((en2==0 && en3==1) ) {// 12-3, 13-2, 23-1
2968 Float_t gFact=1;
2969
2970 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2971 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2972
2973
2974 if(fill2){
2975 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fExplicit3->Fill(q3, gFact);// 12
2976 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2977 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2978 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[1].fNormEx3->Fill(0.);
2979 }
2980 }
2981 }
2982 if(fill3){
2983 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fExplicit3->Fill(q3, gFact);// 12
2984 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2985 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2986 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[2].fNormEx3->Fill(0.);
2987 }
2988 }
2989 }
2990 if(fill4){
2991 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fExplicit3->Fill(q3, gFact);// 12
2992 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
2993 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
2994 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[3].fNormEx3->Fill(0.);
2995 }
2996 }
2997 }
2998
2999 }else if(en2==1 && en3==2){// all uncorrelated events
3000 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
3001
3002 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fExplicit3->Fill(q3);
3003 if((qinv12>=fNormQcutLow[normBin12]) && (qinv13>=fNormQcutLow[normBin13]) && (qinv23>=fNormQcutLow[normBin23])) {
3004 if((qinv12<fNormQcutHigh[normBin12]) && (qinv13<fNormQcutHigh[normBin13]) && (qinv23<fNormQcutHigh[normBin23])) {
3005 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fNormEx3->Fill(0.);
3006 }
3007 }
3008 Short_t qCutBin12 = SetQcutBin(fillIndex12);
3009 Short_t qCutBin13 = SetQcutBin(fillIndex13);
3010 Short_t qCutBin23 = SetQcutBin(fillIndex23);
3011
3012 if( (qinv12 < fQcut[qCutBin12]) || (qinv13 < fQcut[qCutBin13]) || (qinv23 < fQcut[qCutBin23])){
3013
3014 Int_t nUnderCut=0;
3015 if(qinv12<fQcut[qCutBin12]) nUnderCut++;
3016 if(qinv13<fQcut[qCutBin13]) nUnderCut++;
3017 if(qinv23<fQcut[qCutBin23]) nUnderCut++;
3018
3019 }
3020
3021 }else {}
3022
3023
3024 }// 3rd particle
3025 }// 3rd event
3026
3027 }// 2nd particle
3028 }// 2nd event
3029
3030 }// 1st particle
3031
3032
3033
3034
3035 }// End of PdensityExplicit
3036
3037
3038
3039
3040 // Post output data.
3041 PostData(1, fOutputList);
3042
3043}
3044//________________________________________________________________________
3045void AliChaoticity::Terminate(Option_t *)
3046{
3047 // Called once at the end of the query
3048
3049 cout<<"Done"<<endl;
3050
3051}
3052//________________________________________________________________________
3053Bool_t AliChaoticity::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
3054{
3055
3056 if(fabs(first.fEta-second.fEta) > fMinSepTPCEntranceEta) return kTRUE;
3057
3058 // propagate through B field to r=1m
5e3e77d6 3059 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
cd12341d 3060 if(phi1 > 2*PI) phi1 -= 2*PI;
3061 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 3062 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
cd12341d 3063 if(phi2 > 2*PI) phi2 -= 2*PI;
3064 if(phi2 < 0) phi2 += 2*PI;
3065
3066 Float_t deltaphi = phi1 - phi2;
3067 if(deltaphi > PI) deltaphi -= 2*PI;
3068 if(deltaphi < -PI) deltaphi += 2*PI;
3069 deltaphi = fabs(deltaphi);
3070
5e3e77d6 3071 //cout<<deltaphi<<" "<<fMinSepTPCEntrancePhi<<" "<<fMinSepTPCEntranceEta<<endl;
3072 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 3073
3074 // propagate through B field to r=1.6m
5e3e77d6 3075 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
cd12341d 3076 if(phi1 > 2*PI) phi1 -= 2*PI;
3077 if(phi1 < 0) phi1 += 2*PI;
5e3e77d6 3078 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
cd12341d 3079 if(phi2 > 2*PI) phi2 -= 2*PI;
3080 if(phi2 < 0) phi2 += 2*PI;
3081
3082 deltaphi = phi1 - phi2;
3083 if(deltaphi > PI) deltaphi -= 2*PI;
3084 if(deltaphi < -PI) deltaphi += 2*PI;
3085 deltaphi = fabs(deltaphi);
3086
5e3e77d6 3087 if(deltaphi < fMinSepTPCEntrancePhi) return kFALSE;// Min Separation
cd12341d 3088
3089
3090 //
f1cbec0a 3091 /*
cd12341d 3092 Int_t ncl1 = first.fClusterMap.GetNbits();
3093 Int_t ncl2 = second.fClusterMap.GetNbits();
3094 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
3095 Double_t shfrac = 0; Double_t qfactor = 0;
3096 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
3097 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
3098 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
3099 sumQ++;
3100 sumCls+=2;
3101 sumSha+=2;}
3102 else {sumQ--; sumCls+=2;}
3103 }
3104 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
3105 sumQ++;
3106 sumCls++;}
3107 }
3108 if (sumCls>0) {
3109 qfactor = sumQ*1.0/sumCls;
3110 shfrac = sumSha*1.0/sumCls;
3111 }
3112
3113 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
f1cbec0a 3114 */
cd12341d 3115
3116 return kTRUE;
5e3e77d6 3117
cd12341d 3118
3119}
3120//________________________________________________________________________
3121Float_t AliChaoticity::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
3122{
3123 Float_t arg = G_Coeff/qinv;
3124
3125 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
3126 else {return (exp(-arg)-1)/(-arg);}
3127
3128}
3129//________________________________________________________________________
3130void AliChaoticity::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
3131{
3132 Int_t j, k;
3133 Int_t a = i2 - i1;
3134 for (Int_t i = i1; i < i2+1; i++) {
3135 j = (Int_t) (gRandom->Rndm() * a);
3136 k = iarr[j];
3137 iarr[j] = iarr[i];
3138 iarr[i] = k;
3139 }
3140}
3141//________________________________________________________________________
3142Short_t AliChaoticity::FillIndex2part(Short_t key){
3143
3144 if(key==2) return 0;// pi-pi
3145 else if(key==11) return 1;// pi-k
3146 else if(key==101) return 2;// pi-p
3147 else if(key==20) return 3;// k-k
3148 else if(key==110) return 4;// k-p
3149 else return 5;// p-p
3150}
3151//________________________________________________________________________
3152Short_t AliChaoticity::FillIndex3part(Short_t key){
3153
3154 if(key==3) return 0;// pi-pi-pi
3155 else if(key==12) return 1;// pi-pi-k
3156 else if(key==21) return 2;// k-k-pi
3157 else if(key==102) return 3;// pi-pi-p
3158 else if(key==201) return 4;// p-p-pi
3159 else if(key==111) return 5;// pi-k-p
3160 else if(key==30) return 6;// k-k-k
3161 else if(key==120) return 7;// k-k-p
3162 else if(key==210) return 8;// p-p-k
3163 else return 9;// p-p-p
3164
3165}
3166//________________________________________________________________________
3167Short_t AliChaoticity::SetQcutBin(Short_t fi){// fi=FillIndex
3168 if(fi <= 2) return 0;
3169 else if(fi==3) return 1;
3170 else return 2;
3171}
3172//________________________________________________________________________
3173Short_t AliChaoticity::SetNormBin(Short_t fi){// fi=FillIndex
3174 if(fi==0) return 0;
3175 else if(fi <= 2) return 1;
3176 else return 2;
3177}
3178//________________________________________________________________________
3179void AliChaoticity::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
3180
3181 if(fi==0 || fi==3 || fi==5){// Identical species
3182 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
3183 else {b1=c1; b2=c2;}
3184 }else {// Mixed species
3185 if(key1 < key2) { b1=c1; b2=c2;}
3186 else {b1=c2; b2=c1;}
3187 }
3188
3189}
3190//________________________________________________________________________
3191void 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){
3192
3193
3194 // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
3195 Bool_t seSS=kFALSE;
3196 Bool_t seSK=kFALSE;
3197 Short_t seKeySum=0;// only used for pi-k-p case
3198 if(part==1) {// default case (irrelevant for term 1 and term 5)
3199 if(c1==c2) seSS=kTRUE;
3200 if(key1==key2) seSK=kTRUE;
3201 seKeySum = key1+key2;
3202 }
3203 if(part==2){
3204 if(c1==c3) seSS=kTRUE;
3205 if(key1==key3) seSK=kTRUE;
3206 seKeySum = key1+key3;
3207 }
3208
3209
3210 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3211
3212 if(fi==0 || fi==6 || fi==9){// Identical species
3213 if( (c1+c2+c3)==1) {
3214 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
3215 //
3216 if(seSS) fill2=kTRUE;
3217 else {fill3=kTRUE; fill4=kTRUE;}
3218 //
3219 }else if( (c1+c2+c3)==2) {
3220 b1=0; b2=1; b3=1;
3221 //
3222 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
3223 else fill4=kTRUE;
3224 //
3225 }else {
3226 b1=c1; b2=c2; b3=c3;
3227 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
3228 }
3229 }else if(fi != 5){// all the rest except pi-k-p
3230 if(key1==key2){
3231 b3=c3;
3232 if( (c1+c2)==1) {b1=0; b2=1;}
3233 else {b1=c1; b2=c2;}
3234 }else if(key1==key3){
3235 b3=c2;
3236 if( (c1+c3)==1) {b1=0; b2=1;}
3237 else {b1=c1; b2=c3;}
3238 }else {// Key2==Key3
3239 b3=c1;
3240 if( (c2+c3)==1) {b1=0; b2=1;}
3241 else {b1=c2; b2=c3;}
3242 }
3243 //////////////////////////////
3244 if(seSK) fill2=kTRUE;// Same keys from Same Event
3245 else {// Different keys from Same Event
3246 if( (c1+c2+c3)==1) {
3247 if(b3==0) {
3248 if(seSS) fill3=kTRUE;
3249 else fill4=kTRUE;
3250 }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
3251 }else if( (c1+c2+c3)==2) {
3252 if(b3==1) {
3253 if(seSS) fill4=kTRUE;
3254 else fill3=kTRUE;
3255 }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
3256 }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
3257 }
3258 /////////////////////////////
3259 }else {// pi-k-p (no charge ordering applies since all are unique)
3260 if(key1==1){
3261 if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
3262 else {b1=c1; b2=c3; b3=c2;}// pi-p-k
3263 }else if(key1==10){
3264 if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
3265 else {b1=c3; b2=c1; b3=c2;}// k-p-pi
3266 }else {// key1==100
3267 if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
3268 else {b1=c3; b2=c2; b3=c1;}// p-k-pi
3269 }
3270 ////////////////////////////////////
3271 if(seKeySum==11) fill2=kTRUE;
3272 else if(seKeySum==101) fill3=kTRUE;
3273 else fill4=kTRUE;
3274 ////////////////////////////////////
3275 }
3276
3277}
3278//________________________________________________________________________
3279void 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){
3280
3281 // for terms 2-4: start by setting q12(part 1) or q13(part 2)
3282 if(fi==0 || fi==6 || fi==9){// Identical species
3283 if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
3284 if(term==1 || term==5){
3285 if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
3286 else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
3287 else {fQ=q23; sQ=q12; tQ=q13;}
3288 }else if(term==2 && part==1){
3289 fQ=q12; sQ=q13; tQ=q23;
3290 }else if(term==2 && part==2){
3291 fQ=q13; sQ=q12; tQ=q23;
3292 }else if(term==3 && part==1){
3293 sQ=q12;
3294 if(c1==c3) {fQ=q13; tQ=q23;}
3295 else {fQ=q23; tQ=q13;}
3296 }else if(term==3 && part==2){
3297 sQ=q13;
3298 if(c1==c2) {fQ=q12; tQ=q23;}
3299 else {fQ=q23; tQ=q12;}
3300 }else if(term==4 && part==1){
3301 tQ=q12;
3302 if(c1==c3) {fQ=q13; sQ=q23;}
3303 else {fQ=q23; sQ=q13;}
3304 }else if(term==4 && part==2){
3305 tQ=q13;
3306 if(c1==c2) {fQ=q12; sQ=q23;}
3307 else {fQ=q23; sQ=q12;}
3308 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3309 }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
3310 if(term==1 || term==5){
3311 if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
3312 else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
3313 else {tQ=q23; sQ=q12; fQ=q13;}
3314 }else if(term==2 && part==1){
3315 fQ=q12;
3316 if(c1==c3) {tQ=q13; sQ=q23;}
3317 else {tQ=q23; sQ=q13;}
3318 }else if(term==2 && part==2){
3319 fQ=q13;
3320 if(c1==c2) {tQ=q12; sQ=q23;}
3321 else {tQ=q23; sQ=q12;}
3322 }else if(term==3 && part==1){
3323 sQ=q12;
3324 if(c1==c3) {tQ=q13; fQ=q23;}
3325 else {tQ=q23; fQ=q13;}
3326 }else if(term==3 && part==2){
3327 sQ=q13;
3328 if(c1==c2) {tQ=q12; fQ=q23;}
3329 else {tQ=q23; fQ=q12;}
3330 }else if(term==4 && part==1){
3331 tQ=q12; sQ=q13; fQ=q23;
3332 }else if(term==4 && part==2){
3333 tQ=q13; sQ=q12; fQ=q23;
3334 }else cout<<"problem!!!!!!!!!!!!!"<<endl;
3335 }else {// fQ=ss, sQ=ss, tQ=ss
3336 if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
3337 else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
3338 else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
3339 else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
3340 else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
3341 else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
3342 else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
3343 }
3344 }else if(fi != 5){// all the rest except pi-k-p
3345 if(key1==key2){
3346 fQ=q12;
3347 if(c1==c2){
3348 // cases not explicity shown below are not possible
3349 if(term==1 || term==5) {sQ=q13; tQ=q23;}
3350 else if(term==2 && part==1) {sQ=q13; tQ=q23;}
3351 else if(term==3 && part==2) {sQ=q13; tQ=q23;}
3352 else if(term==4 && part==2) {tQ=q13; sQ=q23;}
3353 else cout<<"problem!!!!!!!!!!!!!"<<endl;
3354 }else if(c3==0){
3355 if(c1==c3) {sQ=q13; tQ=q23;}
3356 else {sQ=q23; tQ=q13;}
3357 }else {//c3==1
3358 if(c1==c3) {tQ=q13; sQ=q23;}
3359 else {tQ=q23; sQ=q13;}
3360 }
3361 }else if(key1==key3){
3362 fQ=q13;
3363 if(c1==c3){
3364 // cases not explicity shown below are not possible
3365 if(term==1 || term==5) {sQ=q12; tQ=q23;}
3366 else if(term==2 && part==2) {sQ=q12; tQ=q23;}
3367 else if(term==3 && part==1) {sQ=q12; tQ=q23;}
3368 else if(term==4 && part==1) {tQ=q12; sQ=q23;}
3369 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3370 }else if(c2==0){
3371 if(c1==c2) {sQ=q12; tQ=q23;}
3372 else {sQ=q23; tQ=q12;}
3373 }else {//c2==1
3374 if(c1==c2) {tQ=q12; sQ=q23;}
3375 else {tQ=q23; sQ=q12;}
3376 }
3377 }else {// key2==key3
3378 fQ=q23;
3379 if(c2==c3){
3380 // cases not explicity shown below are not possible
3381 if(term==1 || term==5) {sQ=q12; tQ=q13;}
3382 else if(term==3 && part==1) {sQ=q12; tQ=q13;}
3383 else if(term==3 && part==2) {sQ=q13; tQ=q12;}
3384 else if(term==4 && part==1) {tQ=q12; sQ=q13;}
3385 else if(term==4 && part==2) {tQ=q13; sQ=q12;}
3386 else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
3387 }else if(c1==0){
3388 if(c1==c2) {sQ=q12; tQ=q13;}
3389 else {sQ=q13; tQ=q12;}
3390 }else {//c1==1
3391 if(c1==c2) {tQ=q12; sQ=q13;}
3392 else {tQ=q13; sQ=q12;}
3393 }
3394 }
3395 }else {// pi-k-p
3396 if(key1==1){
3397 if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
3398 else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
3399 }else if(key1==10){
3400 if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
3401 else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
3402 }else {// key1==100
3403 if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
3404 else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
3405 }
3406
3407 }
3408
3409
3410}
3411//________________________________________________________________________
3412Float_t AliChaoticity::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
3413
3414 Float_t qinv=1.0;
3415
3416 if(fi==0 || fi==3 || fi==5){// identical masses
3417 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));
3418 }else{// different masses
3419 Float_t px = track1[1] + track2[1];
3420 Float_t py = track1[2] + track2[2];
3421 Float_t pz = track1[3] + track2[3];
3422 Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
3423 Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
3424 deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
3425
3426 qinv = pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
3427 qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
3428 qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
3429 qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
3430 qinv = sqrt(qinv);
3431 }
3432
3433 return qinv;
3434
3435}
3436//________________________________________________________________________
3437void AliChaoticity::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
3438
3439 Float_t p0 = track1[0] + track2[0];
3440 Float_t px = track1[1] + track2[1];
3441 Float_t py = track1[2] + track2[2];
3442 Float_t pz = track1[3] + track2[3];
3443
3444 Float_t mt = sqrt(p0*p0 - pz*pz);
3445 Float_t pt = sqrt(px*px + py*py);
3446
3447 Float_t v0 = track1[0] - track2[0];
3448 Float_t vx = track1[1] - track2[1];
3449 Float_t vy = track1[2] - track2[2];
3450 Float_t vz = track1[3] - track2[3];
3451
3452 qout = (px*vx + py*vy)/pt;
3453 qside = (px*vy - py*vx)/pt;
3454 qlong = (p0*vz - pz*v0)/mt;
3455}
3456//________________________________________________________________________
90814457 3457//void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F ***histos){
e67e09ed 3458void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F *histos[3][10]){
3459//void AliChaoticity::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliChaoticity::kKbinsT][AliChaoticity::kCentBins]){
704f2481 3460
cd12341d 3461 for(Int_t mb=0; mb<fMbins; mb++){
3462 for(Int_t edB=0; edB<kEDbins; edB++){
3463 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3464 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3465 //
3466 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3467 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3468 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3469
3470 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3471 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 0;
3472 }
3473 }
3474 }
3475 }
3476 }
3477 //
3478 }
3479 }
3480
90814457 3481 if(legoCase){
3482 cout<<"LEGO call to SetWeightArrays"<<endl;
704f2481 3483 // histos[0][0]->GetBinContent(3, 3, 3) should be ~0.14
3484 if(histos[0][0]->GetBinContent(3, 3, 3) > 0.5) AliFatal("AliChaoticity: SetWeightArray Problem");// Additional test to make sure loaded correctly
3485 if(histos[0][0]->GetBinContent(3, 3, 3) < 0.05) AliFatal("AliChaoticity: SetWeightArray Problem");// Additional test to make sure loaded correctly
3486
90814457 3487 for(Int_t mb=0; mb<fMbins; mb++){
3488 for(Int_t edB=0; edB<kEDbins; edB++){
3489 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3490 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
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] = histos[tKbin][mb]->GetBinContent(qobin, qsbin, qlbin);
3497 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = histos[tKbin][mb]->GetBinError(qobin, qsbin, qlbin);
704f2481 3498 if(fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] > 2.0) {// In theory never greater than 1.0
3499 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = 2.0;
3500 }
3501 if(fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] < -0.5) {// In theory never significantly less than 0
3502 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = -0.5;
3503 }
3504
90814457 3505 }
cd12341d 3506 }
3507 }
3508 }
90814457 3509 }
3510 //
3511 }
3512 }
3513 }else{
3514
3515 TFile *wFile = new TFile("WeightFile.root","READ");
3516 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3517 else cout<<"Good Weight File Found!"<<endl;
3518
3519 for(Int_t tKbin=0; tKbin<kKbinsT; tKbin++){
3520 for(Int_t yKbin=0; yKbin<kKbinsY; yKbin++){
3521 for(Int_t mb=0; mb<fMbins; mb++){
3522 for(Int_t edB=0; edB<kEDbins; edB++){
3523
3524 TString *name = new TString("Weight_Kt_");
3525 *name += tKbin;
3526 name->Append("_Ky_");
3527 *name += yKbin;
3528 name->Append("_M_");
3529 *name += mb;
3530 name->Append("_ED_");
3531 *name += edB;
3532
3533 TH3F *fTempHisto = (TH3F*)wFile->Get(name->Data());
3534
3535 for(Int_t qobin=1; qobin<=kQbinsWeights; qobin++){
3536 for(Int_t qsbin=1; qsbin<=kQbinsWeights; qsbin++){
3537 for(Int_t qlbin=1; qlbin<=kQbinsWeights; qlbin++){
3538
3539 fNormWeight[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinContent(qobin, qsbin, qlbin);
3540 fNormWeightErr[mb][edB][tKbin][yKbin][qobin-1][qsbin-1][qlbin-1] = fTempHisto->GetBinError(qobin, qsbin, qlbin);
3541 }
3542 }
3543 }
3544 delete fTempHisto;
3545 }//ED
3546 }//mb
3547 }//ky
3548 }//kt
3549
3550 wFile->Close();
3551 }
3552
5e3e77d6 3553 cout<<"Done reading weight file"<<endl;
cd12341d 3554
3555}
3556//________________________________________________________________________
cd12341d 3557void AliChaoticity::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3558
3559 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3560 Float_t ky=0;// central rapidity
3561 //
3562 Float_t qOut=0,qSide=0,qLong=0;
3563 GetQosl(track1, track2, qOut, qSide, qLong);
3564 qOut = fabs(qOut);
3565 qSide = fabs(qSide);
3566 qLong = fabs(qLong);
3567 //
3568
3569 if(kt < fKmeanT[0]) {fKtbinL=0; fKtbinH=0;}
3570 else if(kt >= fKmeanT[kKbinsT-1]) {fKtbinL=kKbinsT-1; fKtbinH=kKbinsT-1;}
3571 else {
3572 for(Int_t i=0; i<kKbinsT-1; i++){
3573 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtbinL=i; fKtbinH=i+1; break;}
3574 }
3575 }
3576 //
3577 if(ky < fKmeanY[0]) {fKybinL=0; fKybinH=0;}
3578 else if(ky >= fKmeanY[kKbinsY-1]) {fKybinL=kKbinsY-1; fKybinH=kKbinsY-1;}
3579 else {
3580 for(Int_t i=0; i<kKbinsY-1; i++){
3581 if((ky >= fKmeanY[i]) && (ky < fKmeanY[i+1])) {fKybinL=i; fKybinH=i+1; break;}
3582 }
3583 }
3584 //
3585 /////////
3586 if(qOut < fQmean[0]) {fQobinL=0; fQobinH=0;}
3587 else if(qOut >= fQmean[kQbinsWeights-1]) {fQobinL=kQbinsWeights-1; fQobinH=kQbinsWeights-1;}
3588 else {
3589 for(Int_t i=0; i<kQbinsWeights-1; i++){
3590 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQobinL=i; fQobinH=i+1; break;}
3591 }
3592 }
3593 //
3594 if(qSide < fQmean[0]) {fQsbinL=0; fQsbinH=0;}
3595 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsbinL=kQbinsWeights-1; fQsbinH=kQbinsWeights-1;}
3596 else {
3597 for(Int_t i=0; i<kQbinsWeights-1; i++){
3598 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsbinL=i; fQsbinH=i+1; break;}
3599 }
3600 }
3601 //
3602 if(qLong < fQmean[0]) {fQlbinL=0; fQlbinH=0;}
3603 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlbinL=kQbinsWeights-1; fQlbinH=kQbinsWeights-1;}
3604 else {
3605 for(Int_t i=0; i<kQbinsWeights-1; i++){
3606 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlbinL=i; fQlbinH=i+1; break;}
3607 }
3608 }
3609 //
3610
3611
3612 Float_t min = fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3613 Float_t minErr = fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinH];
3614
3615
3616 Float_t deltaW=0;
3617 // kt
3618 deltaW += (fNormWeight[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - min)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3619 // Ky
3620 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - min)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3621 // Qo
5e3e77d6 3622 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - min)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3623 // Qs
5e3e77d6 3624 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - min)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3625 // Ql
5e3e77d6 3626 deltaW += (fNormWeight[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - min)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3627
3628 //
3629 wgt = min + deltaW;
3630
3631
3632 ////
3633
3634 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3635 Float_t deltaWErr=0;
3636 // Kt
3637 /*
3638 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinH][fKybinL][fQobinH][fQsbinH][fQlbinH] - minErr)*(kt-fKmeanT[fKtbinL])/((fKstepT[fKtbinL]+fKstepT[fKtbinH])/2.);
3639 // Ky
3640 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinH][fQobinH][fQsbinH][fQlbinH] - minErr)*(ky-fKmeanY[fKybinL])/((fKstepY[fKybinL]+fKstepY[fKybinH])/2.);
3641 // Qo
5e3e77d6 3642 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinL][fQsbinH][fQlbinH] - minErr)*(qOut-fQmean[fQobinL])/fQstepWeights;
cd12341d 3643 // Qs
5e3e77d6 3644 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinL][fQlbinH] - minErr)*(qSide-fQmean[fQsbinL])/fQstepWeights;
cd12341d 3645 // Ql
5e3e77d6 3646 deltaWErr += (fNormWeightErr[fMbin][0][fKtbinL][fKybinL][fQobinH][fQsbinH][fQlbinL] - minErr)*(qLong-fQmean[fQlbinL])/fQstepWeights;
cd12341d 3647 */
3648 wgtErr = minErr + deltaWErr;
3649
3650
3651
cd12341d 3652}
3653//________________________________________________________________________
3654Float_t AliChaoticity::MCWeight(Int_t charge1, Int_t charge2, Int_t rIndex, Int_t dampIndex, Float_t qinv){
3655
5e3e77d6 3656 Float_t radius = Float_t(rIndex+3.)/0.19733;// convert to GeV
cd12341d 3657 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
90814457 3658 Float_t coulCorr12 = FSICorrelationGaus2(charge1, charge2, rIndex, qinv);
cd12341d 3659 if(charge1==charge2){
5e3e77d6 3660 return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*radius,2)))*coulCorr12);
cd12341d 3661 }else {
5e3e77d6 3662 return ((1-myDamp) + myDamp*coulCorr12);
cd12341d 3663 }
3664
3665}
3666//________________________________________________________________________
90814457 3667Float_t AliChaoticity::MCWeight3D(Bool_t SameCharge, Int_t term, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23){
cd12341d 3668 if(term==5) return 1.0;
5e3e77d6 3669
90814457 3670 Float_t radius=5;
3671 if(fMbin<=1) {radius = 8;}
3672 else if(fMbin<=3) {radius = 7;}
3673 else if(fMbin<=5) {radius = 6;}
3674 else {radius = 5;}
3675 radius /= 0.19733;
3676
3677 //Float_t radius = (3. + rIndex)/0.19733;//starts at 3fm. convert to GeV
cd12341d 3678 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3679 Float_t fc = sqrt(myDamp);
5e3e77d6 3680 if(SameCharge){// all three of the same charge
90814457 3681 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2
3682 Float_t coulCorr13 = FSICorrelationTherm2(+1,+1, q13);// K2
3683 Float_t coulCorr23 = FSICorrelationTherm2(+1,+1, q23);// K2
cd12341d 3684
5e3e77d6 3685 if(term==1){
3686 Float_t c3QS = 1 + exp(-pow(q12*radius,2)) + exp(-pow(q13*radius,2)) + exp(-pow(q23*radius,2));
3687 c3QS += 2*exp(-pow(radius,2)*(pow(q12,2) + pow(q13,2) + pow(q23,2))/2.);
3688 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3689 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3690 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*radius,2)))*coulCorr13;
3691 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*radius,2)))*coulCorr23;
3692 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kTRUE, q12, q13, q23);
3693 return w123;
3694 }else if(term==2){
3695 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3696 }else if(term==3){
3697 return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*radius,2)))*coulCorr13);
3698 }else if(term==4){
3699 return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*radius,2)))*coulCorr23);
3700 }else return 1.0;
3701
3702 }else{// mixed charge case pair 12 always treated as ss
90814457 3703 Float_t coulCorr12 = FSICorrelationTherm2(+1,+1, q12);// K2 ss
3704 Float_t coulCorr13 = FSICorrelationTherm2(+1,-1, q13);// K2 os
3705 Float_t coulCorr23 = FSICorrelationTherm2(+1,-1, q23);// K2 os
5e3e77d6 3706 if(term==1){
3707 Float_t c3QS = 1 + exp(-pow(q12*radius,2));
3708 Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3709 w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2)))*coulCorr12;
3710 w123 += pow(fc,2)*(1-fc)*coulCorr13;
3711 w123 += pow(fc,2)*(1-fc)*coulCorr23;
3712 w123 += pow(fc,3)*c3QS*FSICorrelationOmega0(kFALSE, q12, q13, q23);
3713 return w123;
3714 }else if(term==2){
3715 return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2)))*coulCorr12);
3716 }else if(term==3){
3717 return ((1-myDamp) + myDamp*coulCorr13);
3718 }else if(term==4){
3719 return ((1-myDamp) + myDamp*coulCorr23);
3720 }else return 1.0;
3721 }
3722
cd12341d 3723}
3724//________________________________________________________________________
90814457 3725void AliChaoticity::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D){
5e3e77d6 3726
3727
3728 if(legoCase){
90814457 3729 cout<<"LEGO call to SetMomResCorrections"<<endl;
5e3e77d6 3730 fMomResC2 = (TH2D*)temp2D->Clone();
5e3e77d6 3731 fMomResC2->SetDirectory(0);
5e3e77d6 3732 }else {
3733 TFile *momResFile = new TFile("MomResFile.root","READ");
3734 if(!momResFile->IsOpen()) {
3735 cout<<"No momentum resolution file found"<<endl;
3736 AliFatal("No momentum resolution file found. Kill process.");
3737 }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
3738
3739 TH2D *temp2D2 = (TH2D*)momResFile->Get("MomResHisto_pp");
5e3e77d6 3740 fMomResC2 = (TH2D*)temp2D2->Clone();
5e3e77d6 3741 fMomResC2->SetDirectory(0);
90814457 3742
5e3e77d6 3743 momResFile->Close();
3744 }
cd12341d 3745
704f2481 3746 // fMomResC2->GetBinContent(1,5) should be ~1.007
3747 if(fMomResC2->GetBinContent(1,5) > 1.2) AliFatal("AliChaoticity: SetMomResCorrections Problem");// Additional Safety check
3748 if(fMomResC2->GetBinContent(1,5) < 0.95) AliFatal("AliChaoticity: SetMomResCorrections Problem");// Additional Safety check
3749
3750 for(Int_t bx=1; bx<=fMomResC2->GetNbinsX(); bx++){
3751 for(Int_t by=1; by<=fMomResC2->GetNbinsX(); by++){
3752 if(fMomResC2->GetBinContent(bx,by) > 1.5) fMomResC2->SetBinContent(bx,by, 1.5);// Maximum is ~1.02
3753 if(fMomResC2->GetBinContent(bx,by) < 0.95) fMomResC2->SetBinContent(bx,by, 0.95);// Minimum is ~0.98
3754 }
3755 }
3756
5e3e77d6 3757 cout<<"Done reading momentum resolution file"<<endl;
3758}
3759//________________________________________________________________________
90814457 3760void AliChaoticity::SetFSICorrelations(Bool_t legoCase, TH2D *temp2DGaus[2], TH2D *temp2DTherm[2], TH3D *temp3Dos[6], TH3D *temp3Dss[6]){
5e3e77d6 3761 // read in 2-particle and 3-particle FSI correlations = K2 & K3
3762 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
3763 if(legoCase){
c4980714 3764 cout<<"LEGO call to SetFSICorrelations"<<endl;
90814457 3765 fFSI2SS[0] = (TH2D*)temp2DGaus[0]->Clone();
3766 fFSI2OS[0] = (TH2D*)temp2DGaus[1]->Clone();
3767 fFSI2SS[1] = (TH2D*)temp2DTherm[0]->Clone();
3768 fFSI2OS[1] = (TH2D*)temp2DTherm[1]->Clone();
5e3e77d6 3769 //
90814457 3770 fFSI2SS[0]->SetDirectory(0);
3771 fFSI2OS[0]->SetDirectory(0);
3772 fFSI2SS[1]->SetDirectory(0);
3773 fFSI2OS[1]->SetDirectory(0);
3774
3775 for(Int_t CB=0; CB<6; CB++) {
3776 fFSIOmega0OS[CB] = (TH3D*)temp3Dos[CB]->Clone();
3777 fFSIOmega0SS[CB] = (TH3D*)temp3Dss[CB]->Clone();
3778 //
3779 fFSIOmega0OS[CB]->SetDirectory(0);
3780 fFSIOmega0SS[CB]->SetDirectory(0);
3781 }
5e3e77d6 3782 }else {
c4980714 3783 cout<<"non LEGO call to SetFSICorrelations"<<endl;
5e3e77d6 3784 TFile *fsifile = new TFile("KFile.root","READ");
3785 if(!fsifile->IsOpen()) {
3786 cout<<"No FSI file found"<<endl;
3787 AliFatal("No FSI file found. Kill process.");
3788 }else {cout<<"Good FSI File Found!"<<endl;}
3789
90814457 3790 TH2D *temphisto2GausSS = (TH2D*)fsifile->Get("K2ssG");
3791 TH2D *temphisto2GausOS = (TH2D*)fsifile->Get("K2osG");
3792 TH2D *temphisto2ThermSS = (TH2D*)fsifile->Get("K2ssT");
3793 TH2D *temphisto2ThermOS = (TH2D*)fsifile->Get("K2osT");
3794 TH3D *temphisto3OS[6];
ae9b34d1 3795 TH3D *temphisto3SS[6];
3796 for(Int_t CB=0; CB<6; CB++) {
90814457 3797 TString *nameK3SS = new TString("K3ss_");
3798 *nameK3SS += CB;
3799 temphisto3SS[CB] = (TH3D*)fsifile->Get(nameK3SS->Data());
3800 //
3801 TString *nameK3OS = new TString("K3os_");
3802 *nameK3OS += CB;
3803 temphisto3OS[CB] = (TH3D*)fsifile->Get(nameK3OS->Data());
ae9b34d1 3804 }
3805
90814457 3806 fFSI2SS[0] = (TH2D*)temphisto2GausSS->Clone();
3807 fFSI2OS[0] = (TH2D*)temphisto2GausOS->Clone();
3808 fFSI2SS[1] = (TH2D*)temphisto2ThermSS->Clone();
3809 fFSI2OS[1] = (TH2D*)temphisto2ThermOS->Clone();
3810 fFSI2SS[0]->SetDirectory(0);
3811 fFSI2OS[0]->SetDirectory(0);
3812 fFSI2SS[1]->SetDirectory(0);
3813 fFSI2OS[1]->SetDirectory(0);
ae9b34d1 3814
90814457 3815 for(Int_t CB=0; CB<6; CB++) {
3816 fFSIOmega0SS[CB] = (TH3D*)temphisto3SS[CB]->Clone();
3817 fFSIOmega0OS[CB] = (TH3D*)temphisto3OS[CB]->Clone();
3818 fFSIOmega0SS[CB]->SetDirectory(0);
3819 fFSIOmega0OS[CB]->SetDirectory(0);
3820 }
3821 //
3822
5e3e77d6 3823 fsifile->Close();
cd12341d 3824 }
3825
5e3e77d6 3826 // condition FSI histogram for edge effects
ae9b34d1 3827 for(Int_t CB=0; CB<6; CB++){
3828 for(Int_t ii=1; ii<=fFSIOmega0SS[CB]->GetNbinsX(); ii++){
3829 for(Int_t jj=1; jj<=fFSIOmega0SS[CB]->GetNbinsY(); jj++){
3830 for(Int_t kk=1; kk<=fFSIOmega0SS[CB]->GetNbinsZ(); kk++){
5e3e77d6 3831
ae9b34d1 3832 if(fFSIOmega0SS[CB]->GetBinContent(ii,jj,kk) <=0){
3833 Double_t Q12 = fFSIOmega0SS[CB]->GetXaxis()->GetBinCenter(ii);
3834 Double_t Q23 = fFSIOmega0SS[CB]->GetYaxis()->GetBinCenter(jj);
3835 Double_t Q13 = fFSIOmega0SS[CB]->GetZaxis()->GetBinCenter(kk);
3836 //
3837 Int_t Q12bin=ii;
3838 Int_t Q23bin=jj;
3839 Int_t Q13bin=kk;
3840 Int_t AC=0;//Adjust Counter
3841 Int_t AClimit=10;// maximum bin shift
3842 if(Q12 < sqrt(pow(Q13,2)+pow(Q23,2) - 2*Q13*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin++; AC++;}}
3843 if(Q12 > sqrt(pow(Q13,2)+pow(Q23,2) + 2*Q13*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q12bin--; AC++;}}
3844 //
3845 if(Q13 < sqrt(pow(Q12,2)+pow(Q23,2) - 2*Q12*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin++; AC++;}}
3846 if(Q13 > sqrt(pow(Q12,2)+pow(Q23,2) + 2*Q12*Q23)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q13bin--; AC++;}}
3847 //
3848 if(Q23 < sqrt(pow(Q12,2)+pow(Q13,2) - 2*Q12*Q13)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin++; AC++;}}
3849 if(Q23 > sqrt(pow(Q12,2)+pow(Q13,2) + 2*Q12*Q13)) {while(fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0 && AC<AClimit) {Q23bin--; AC++;}}
3850
3851 // Save cpu time by setting empty cell contents (edge effects) to nearest non-zero cell (these cells are not used very often anyway.)
3852 if(AC==AClimit) {
3853 fFSIOmega0SS[CB]->SetBinContent(ii,jj,kk, 1.0);
90814457 3854 fFSIOmega0OS[CB]->SetBinContent(ii,jj,kk, 1.0);
ae9b34d1 3855 }else {
3856 fFSIOmega0SS[CB]->SetBinContent(ii,jj,kk, fFSIOmega0SS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin));
90814457 3857 fFSIOmega0OS[CB]->SetBinContent(ii,jj,kk, fFSIOmega0OS[CB]->GetBinContent(Q12bin, Q23bin, Q13bin));
ae9b34d1 3858 }
5e3e77d6 3859 }
ae9b34d1 3860
5e3e77d6 3861 }
5e3e77d6 3862 }
cd12341d 3863 }
3864 }
704f2481 3865 // fFSI2SS[1]->GetBinContent(1,2) should be ~0.32
3866 if(fFSI2SS[1]->GetBinContent(1,2) > 1.0) AliFatal("AliChaoticity: SetFSICorrelations Problem");// Additional Safety check
3867 if(fFSI2SS[1]->GetBinContent(1,2) < 0.1) AliFatal("AliChaoticity: SetFSICorrelations Problem");// Additional Safety check
3868
3869 for(Int_t ii=1; ii<=fFSI2SS[0]->GetNbinsX(); ii++){
3870 for(Int_t jj=1; jj<=fFSI2SS[0]->GetNbinsY(); jj++){
3871 if(fFSI2SS[0]->GetBinContent(ii,jj) > 1.0) fFSI2SS[0]->SetBinContent(ii,jj, 1.0);
3872 if(fFSI2SS[1]->GetBinContent(ii,jj) > 1.0) fFSI2SS[1]->SetBinContent(ii,jj, 1.0);
3873 if(fFSI2OS[0]->GetBinContent(ii,jj) > 10.0) fFSI2OS[0]->SetBinContent(ii,jj, 10.0);
3874 if(fFSI2OS[1]->GetBinContent(ii,jj) > 10.0) fFSI2OS[1]->SetBinContent(ii,jj, 10.0);
3875 //
3876 if(fFSI2SS[0]->GetBinContent(ii,jj) < 0.05) fFSI2SS[0]->SetBinContent(ii,jj, 0.05);
3877 if(fFSI2SS[1]->GetBinContent(ii,jj) < 0.05) fFSI2SS[1]->SetBinContent(ii,jj, 0.05);
3878 if(fFSI2OS[0]->GetBinContent(ii,jj) < 0.9) fFSI2OS[0]->SetBinContent(ii,jj, 0.9);
3879 if(fFSI2OS[1]->GetBinContent(ii,jj) < 0.9) fFSI2OS[1]->SetBinContent(ii,jj, 0.9);
3880 }
3881 }
3882
5e3e77d6 3883 cout<<"Done reading FSI file"<<endl;
cd12341d 3884}
3885//________________________________________________________________________
90814457 3886Float_t AliChaoticity::FSICorrelationGaus2(Int_t charge1, Int_t charge2, Int_t rIndex, Float_t qinv){
5e3e77d6 3887 // returns 2-particle Coulomb correlations = K2
3888 if(rIndex >= kRVALUES) return 1.0;
90814457 3889 Int_t qbinL = fFSI2SS[0]->GetYaxis()->FindBin(qinv-fFSI2SS[0]->GetYaxis()->GetBinWidth(1)/2.);
3890 Int_t qbinH = qbinL+1;
3891 if(qbinL <= 0) return 1.0;
3892 if(qbinH > fFSI2SS[0]->GetNbinsY()) return 1.0;
3893
3894 Float_t slope=0;
3895 if(charge1==charge2){
3896 slope = fFSI2SS[0]->GetBinContent(rIndex+1, qbinL) - fFSI2SS[0]->GetBinContent(rIndex+1, qbinH);
3897 slope /= fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinH);
3898 return (slope*(qinv - fFSI2SS[0]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS[0]->GetBinContent(rIndex+1, qbinL));
3899 }else {
3900 slope = fFSI2OS[0]->GetBinContent(rIndex+1, qbinL) - fFSI2OS[0]->GetBinContent(rIndex+1, qbinH);
3901 slope /= fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinH);
3902 return (slope*(qinv - fFSI2OS[0]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS[0]->GetBinContent(rIndex+1, qbinL));
3903 }
3904}
3905//________________________________________________________________________
3906Float_t AliChaoticity::FSICorrelationTherm2(Int_t charge1, Int_t charge2, Float_t qinv){
3907 // returns 2-particle Coulomb correlations = K2
3908 Int_t qbinL = fFSI2SS[1]->GetYaxis()->FindBin(qinv-fFSI2SS[1]->GetYaxis()->GetBinWidth(1)/2.);
5e3e77d6 3909 Int_t qbinH = qbinL+1;
3910 if(qbinL <= 0) return 1.0;
90814457 3911 if(qbinH > fFSI2SS[1]->GetNbinsY()) return 1.0;
5e3e77d6 3912
3913 Float_t slope=0;
3914 if(charge1==charge2){
90814457 3915 slope = fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinL) - fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinH);
3916 slope /= fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinL) - fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinH);
3917 return (slope*(qinv - fFSI2SS[1]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2SS[1]->GetBinContent(fFSIbin+1, qbinL));
5e3e77d6 3918 }else {
90814457 3919 slope = fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinL) - fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinH);
3920 slope /= fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinL) - fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinH);
3921 return (slope*(qinv - fFSI2OS[1]->GetYaxis()->GetBinCenter(qbinL)) + fFSI2OS[1]->GetBinContent(fFSIbin+1, qbinL));
cd12341d 3922 }
3923}
3924//________________________________________________________________________
5e3e77d6 3925Double_t AliChaoticity::FSICorrelationOmega0(Bool_t SameCharge, Double_t Q12, Double_t Q13, Double_t Q23){
3926 // remember: Omega0 histogram is in the following order (Q12, Q23, Q13)
3927 // returns 3d 3-particle Coulomb Correlation = K3
ae9b34d1 3928 Int_t Q12bin = fFSIOmega0SS[fFSIbin]->GetXaxis()->FindBin(Q12);
3929 Int_t Q13bin = fFSIOmega0SS[fFSIbin]->GetZaxis()->FindBin(Q13);
3930 Int_t Q23bin = fFSIOmega0SS[fFSIbin]->GetYaxis()->FindBin(Q23);
5e3e77d6 3931
3932 if(SameCharge){
ae9b34d1 3933 if(fFSIOmega0SS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3934 else return fFSIOmega0SS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
90814457 3935 }else{// mixed charge.
3936 if(fFSIOmega0OS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin) <=0) return 1.0;
3937 else return fFSIOmega0OS[fFSIbin]->GetBinContent(Q12bin, Q23bin, Q13bin);// K3
5e3e77d6 3938 }
cd12341d 3939}
90814457 3940//________________________________________________________________________
3941void 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){
3942 QS1v1 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
3943 QS1v1 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
3944 QS1v1 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
3945 QS2 = (pV1[1]-pV2[1])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[1]-pV3[1]);
3946 QS3v1 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
3947 QS3v1 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
3948 //
3949 QS1v2 = (pV1[0]-pV2[0])*(pV2[1]-pV3[1]) - (pV1[1]-pV2[1])*(pV2[0]-pV3[0]);
3950 QS1v2 += (pV1[0]-pV2[0])*(pV2[2]-pV3[2]) - (pV1[2]-pV2[2])*(pV2[0]-pV3[0]);
3951 QS3v2 = (pV1[1]-pV2[1])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[1]-pV3[1]);
3952 QS3v2 += (pV1[0]-pV2[0])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[0]-pV3[0]);
3953 QS3v2 += (pV1[2]-pV2[2])*(pV2[3]-pV3[3]) - (pV1[3]-pV2[3])*(pV2[2]-pV3[2]);
3954}
3955//________________________________________________________________________
3956void AliChaoticity::TestAddTask(){
3957 /*
3958 TString inputFileNameFSI = "KFile.root";
3959 TFile *inputFileFSI = TFile::Open(inputFileNameFSI,"OLD");
3960 if (!inputFileFSI){
3961 cout << "Requested file:" << inputFileFSI << " was not opened. ABORT." << endl;
3962 return;
3963 }
3964 TH2D *FSI2gaus[2];
3965 TH2D *FSI2therm[2];
3966 TH3D *FSI3ss[6];
3967 TH3D *FSI3os[6];
704f2481 3968
90814457 3969 FSI2gaus[0] = (TH2D*)inputFileFSI->Get("K2ssG");
3970 FSI2gaus[1] = (TH2D*)inputFileFSI->Get("K2osG");
3971 FSI2therm[0] = (TH2D*)inputFileFSI->Get("K2ssT");
3972 FSI2therm[1] = (TH2D*)inputFileFSI->Get("K2osT");
3973 for(Int_t CB=0; CB<6; CB++) {
3974 TString *nameSS=new TString("K3ss_");
3975 *nameSS += CB;
3976 FSI3ss[CB] = (TH3D*)inputFileFSI->Get(nameSS->Data());
3977 TString *nameOS=new TString("K3os_");
3978 *nameOS += CB;
3979 FSI3os[CB] = (TH3D*)inputFileFSI->Get(nameOS->Data());
3980 }
3981 //
3982 FSI2gaus[0]->SetDirectory(0);
3983 FSI2gaus[1]->SetDirectory(0);
3984 FSI2therm[0]->SetDirectory(0);
3985 FSI2therm[1]->SetDirectory(0);
3986 for(Int_t CB=0; CB<6; CB++) {
3987 FSI3ss[CB]->SetDirectory(0);
3988 FSI3os[CB]->SetDirectory(0);
3989 }
704f2481 3990
90814457 3991 SetFSICorrelations( kTRUE, FSI2gaus, FSI2therm , FSI3os, FSI3ss);
3992 //
3993
3994 if(!fTabulatePairs) {
3995 TString inputFileNameWeight = "WeightFile.root";
3996 TFile *inputFileWeight = TFile::Open(inputFileNameWeight,"OLD");
3997 if (!inputFileWeight){
3998 cout << "Requested file:" << inputFileWeight << " was not opened. ABORT." << endl;
3999 return;
4000 }
4001
4002 ////////////////////////////////////////////////////
4003 // C2 Weight File
704f2481 4004 const Int_t ktbins = 3;
4005 const Int_t cbins = 10;
90814457 4006 TH3F *weightHisto[ktbins][cbins];
4007 for(Int_t i=0; i<ktbins; i++){
4008 for(Int_t j=0; j<cbins; j++){
4009 TString name = "Weight_Kt_";
4010 name += i;
4011 name += "_Ky_0_M_";
4012 name += j;
4013 name += "_ED_0";
4014
4015 weightHisto[i][j] = (TH3F*)inputFileWeight->Get(name);
4016 }
4017 }
4018 SetWeightArrays( kTRUE, weightHisto );
4019 }
4020
4021 //
4022 if(!fMCcase && !fTabulatePairs){
4023 TString inputFileNameMomRes = "MomResFile.root";
4024 TFile *inputFileMomRes = TFile::Open(inputFileNameMomRes,"OLD");
4025 if (!inputFileMomRes){
4026 cout << "Requested file:" << inputFileMomRes << " was not opened. ABORT." << endl;
4027 return;
4028 }
4029 ////////////////////////////////////////////////////
4030 // Momentum Resolution File
4031 TH2D *momResHisto2D = 0;
4032 momResHisto2D = (TH2D*)inputFileMomRes->Get("MomResHisto_pp");
4033 SetMomResCorrections( kTRUE, momResHisto2D);
4034 }
4035 */
4036}