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