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