]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
code updates, Add Muon Corrections, KT4index bug fix
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / Chaoticity / AliFourPion.cxx
CommitLineData
be9ef9f9 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#include "AliAnalysisUtils.h"
32
33#include "AliFourPion.h"
34
35#define PI 3.1415927
36#define G_Coeff 0.006399 // 2*pi*alpha*M_pion
37#define kappa3 0.1 // kappa3 Edgeworth coefficient (non-Gaussian features of C2)
38#define kappa4 0.5 // kappa4 Edgeworth coefficient (non-Gaussian features of C2)
39
40
41// Author: Dhevan Gangadharan
42
43ClassImp(AliFourPion)
44
45//________________________________________________________________________
46AliFourPion::AliFourPion():
47AliAnalysisTaskSE(),
48 fname(0),
49 fAOD(0x0),
50 fOutputList(0x0),
51 fPIDResponse(0x0),
52 fEC(0x0),
53 fEvt(0x0),
54 fTempStruct(0x0),
55 fRandomNumber(0x0),
56 fLEGO(kTRUE),
57 fMCcase(kFALSE),
58 fAODcase(kTRUE),
59 fPbPbcase(kTRUE),
60 fGenerateSignal(kFALSE),
61 fGeneratorOnly(kFALSE),
62 fPdensityPairCut(kTRUE),
63 fTabulatePairs(kFALSE),
64 fRMax(11),
65 ffcSq(0.7),
66 fFilterBit(7),
67 fMaxChi2NDF(10),
68 fMinTPCncls(0),
69 fBfield(0),
70 fMbin(0),
71 fFSIindex(0),
72 fEDbin(0),
73 fMbins(fCentBins),
74 fMultLimit(0),
75 fCentBinLowLimit(0),
76 fCentBinHighLimit(1),
77 fEventCounter(0),
78 fEventsToMix(0),
79 fZvertexBins(0),
80 fMultLimits(),
81 fQcut(0),
82 fQLowerCut(0),
83 fNormQcutLow(0),
84 fNormQcutHigh(0),
85 fKupperBound(0),
86 fQupperBoundQ2(0),
87 fQupperBoundQ3(0),
88 fQupperBoundQ4(0),
89 fQbinsQ2(1),
90 fQbinsQ3(1),
91 fQbinsQ4(1),
92 fQupperBoundWeights(0),
93 fKstepT(),
94 fKstepY(),
95 fKmeanT(),
96 fKmeanY(),
97 fKmiddleT(),
98 fKmiddleY(),
99 fQstep(0),
100 fQstepWeights(0),
101 fQmean(),
102 fDampStart(0),
103 fDampStep(0),
104 fTPCTOFboundry(0),
105 fTOFboundry(0),
106 fSigmaCutTPC(2.0),
107 fSigmaCutTOF(2.0),
108 fMinSepPairEta(0.03),
109 fMinSepPairPhi(0.04),
110 fShareQuality(0),
111 fShareFraction(0),
112 fTrueMassP(0),
113 fTrueMassPi(0),
114 fTrueMassK(0),
115 fTrueMassKs(0),
116 fTrueMassLam(0),
117 fKtIndexL(0),
118 fKtIndexH(0),
119 fQoIndexL(0),
120 fQoIndexH(0),
121 fQsIndexL(0),
122 fQsIndexH(0),
123 fQlIndexL(0),
124 fQlIndexH(0),
125 fDummyB(0),
126 fKT3transition(0.3),
127 fKT4transition(0.3),
b71263d0 128 fDefaultsCharSwitch(),
129 fLowQPairSwitch_E0E0(),
130 fLowQPairSwitch_E0E1(),
131 fLowQPairSwitch_E0E2(),
132 fLowQPairSwitch_E0E3(),
6bb6954b 133 fLowQPairSwitch_E1E1(),
b71263d0 134 fLowQPairSwitch_E1E2(),
135 fLowQPairSwitch_E1E3(),
136 fLowQPairSwitch_E2E3(),
137 fNormQPairSwitch_E0E0(),
138 fNormQPairSwitch_E0E1(),
139 fNormQPairSwitch_E0E2(),
140 fNormQPairSwitch_E0E3(),
6bb6954b 141 fNormQPairSwitch_E1E1(),
b71263d0 142 fNormQPairSwitch_E1E2(),
143 fNormQPairSwitch_E1E3(),
144 fNormQPairSwitch_E2E3(),
5591748e 145 fMomResC2(0x0),
146 fWeightmuonCorrection(0x0)
be9ef9f9 147{
148 // Default constructor
149 for(Int_t mb=0; mb<fMbins; mb++){
150 for(Int_t edB=0; edB<fEDbins; edB++){
151 for(Int_t c1=0; c1<2; c1++){
152 for(Int_t c2=0; c2<2; c2++){
153 for(Int_t term=0; term<2; term++){
154
155 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
156
157 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
158 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
159 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
160 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
161 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
162 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
163
164 }// term_2
165
166
167 for(Int_t c3=0; c3<2; c3++){
168 for(Int_t term=0; term<5; term++){
169
170 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
171 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
b71263d0 172 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
be9ef9f9 173 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
174
175 }// term_3
b71263d0 176
177 for(Int_t c4=0; c4<2; c4++){
6bb6954b 178 for(Int_t term=0; term<13; term++){
b71263d0 179
180 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
181 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
182 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
183 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
184
185 }// term_4
186
187 }// c4
be9ef9f9 188 }//c3
189 }//c2
190 }//c1
191 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
192 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
193 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
194 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
195 }
196 }
197
198 }// ED
199 }// Mbin
200
201 // Initialize FSI histograms
202 for(Int_t i=0; i<12; i++){
203 fFSIss[i]=0x0;
204 fFSIos[i]=0x0;
205 }
206
207
208 // Initialize fNormWeight and fNormWeightErr to 0
209 for(Int_t i=0; i<3; i++){// Kt iterator
210 for(Int_t j=0; j<10; j++){// Mbin iterator
211 fNormWeight[i][j]=0x0;
212 }
213 }
214
215
216}
217//________________________________________________________________________
218AliFourPion::AliFourPion(const Char_t *name)
219: AliAnalysisTaskSE(name),
220 fname(name),
221 fAOD(0x0),
222 fOutputList(0x0),
223 fPIDResponse(0x0),
224 fEC(0x0),
225 fEvt(0x0),
226 fTempStruct(0x0),
227 fRandomNumber(0x0),
228 fLEGO(kTRUE),
229 fMCcase(kFALSE),
230 fAODcase(kTRUE),
231 fPbPbcase(kTRUE),
232 fGenerateSignal(kFALSE),
233 fGeneratorOnly(kFALSE),
234 fPdensityPairCut(kTRUE),
235 fTabulatePairs(kFALSE),
236 fRMax(11),
237 ffcSq(0.7),
238 fFilterBit(7),
239 fMaxChi2NDF(10),
240 fMinTPCncls(0),
241 fBfield(0),
242 fMbin(0),
243 fFSIindex(0),
244 fEDbin(0),
245 fMbins(fCentBins),
246 fMultLimit(0),
247 fCentBinLowLimit(0),
248 fCentBinHighLimit(1),
249 fEventCounter(0),
250 fEventsToMix(0),
251 fZvertexBins(0),
252 fMultLimits(),
253 fQcut(0),
254 fQLowerCut(0),
255 fNormQcutLow(0),
256 fNormQcutHigh(0),
257 fKupperBound(0),
258 fQupperBoundQ2(0),
259 fQupperBoundQ3(0),
260 fQupperBoundQ4(0),
261 fQbinsQ2(1),
262 fQbinsQ3(1),
263 fQbinsQ4(1),
264 fQupperBoundWeights(0),
265 fKstepT(),
266 fKstepY(),
267 fKmeanT(),
268 fKmeanY(),
269 fKmiddleT(),
270 fKmiddleY(),
271 fQstep(0),
272 fQstepWeights(0),
273 fQmean(),
274 fDampStart(0),
275 fDampStep(0),
276 fTPCTOFboundry(0),
277 fTOFboundry(0),
278 fSigmaCutTPC(2.0),
279 fSigmaCutTOF(2.0),
280 fMinSepPairEta(0.03),
281 fMinSepPairPhi(0.04),
282 fShareQuality(0),
283 fShareFraction(0),
284 fTrueMassP(0),
285 fTrueMassPi(0),
286 fTrueMassK(0),
287 fTrueMassKs(0),
288 fTrueMassLam(0),
289 fKtIndexL(0),
290 fKtIndexH(0),
291 fQoIndexL(0),
292 fQoIndexH(0),
293 fQsIndexL(0),
294 fQsIndexH(0),
295 fQlIndexL(0),
296 fQlIndexH(0),
297 fDummyB(0),
298 fKT3transition(0.3),
299 fKT4transition(0.3),
b71263d0 300 fDefaultsCharSwitch(),
301 fLowQPairSwitch_E0E0(),
302 fLowQPairSwitch_E0E1(),
303 fLowQPairSwitch_E0E2(),
304 fLowQPairSwitch_E0E3(),
6bb6954b 305 fLowQPairSwitch_E1E1(),
b71263d0 306 fLowQPairSwitch_E1E2(),
307 fLowQPairSwitch_E1E3(),
308 fLowQPairSwitch_E2E3(),
309 fNormQPairSwitch_E0E0(),
310 fNormQPairSwitch_E0E1(),
311 fNormQPairSwitch_E0E2(),
312 fNormQPairSwitch_E0E3(),
6bb6954b 313 fNormQPairSwitch_E1E1(),
b71263d0 314 fNormQPairSwitch_E1E2(),
315 fNormQPairSwitch_E1E3(),
316 fNormQPairSwitch_E2E3(),
5591748e 317 fMomResC2(0x0),
318 fWeightmuonCorrection(0x0)
be9ef9f9 319{
320 // Main constructor
321 fAODcase=kTRUE;
322 fPdensityPairCut = kTRUE;
323
324
325 for(Int_t mb=0; mb<fMbins; mb++){
326 for(Int_t edB=0; edB<fEDbins; edB++){
327 for(Int_t c1=0; c1<2; c1++){
328 for(Int_t c2=0; c2<2; c2++){
329 for(Int_t term=0; term<2; term++){
330
331 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
332
333 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
334 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
335 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
336 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
337 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
338 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
339
340 }// term_2
341
342 for(Int_t c3=0; c3<2; c3++){
343 for(Int_t term=0; term<5; term++){
344
345 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
346 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
b71263d0 347 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
be9ef9f9 348 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
349
350 }// term_3
b71263d0 351
352 for(Int_t c4=0; c4<2; c4++){
6bb6954b 353 for(Int_t term=0; term<13; term++){
b71263d0 354
355 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
356 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
357 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
358 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
359
360 }// term_4
361 }// c4
be9ef9f9 362 }//c3
363 }//c2
364 }//c1
365
366 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
367 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
368 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
369 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
370 }
371 }
372
373 }// ED
374 }// Mbin
375
376 // Initialize FSI histograms
377 for(Int_t i=0; i<12; i++){
378 fFSIss[i]=0x0;
379 fFSIos[i]=0x0;
380 }
381
382 // Initialize fNormWeight and fNormWeightErr to 0
383 for(Int_t i=0; i<3; i++){// Kt iterator
384 for(Int_t j=0; j<10; j++){// Mbin iterator
385 fNormWeight[i][j]=0x0;
386 }
387 }
388
389
390 DefineOutput(1, TList::Class());
391}
392//________________________________________________________________________
393AliFourPion::AliFourPion(const AliFourPion &obj)
394 : AliAnalysisTaskSE(obj.fname),
395 fname(obj.fname),
396 fAOD(obj.fAOD),
397 //fESD(obj.fESD),
398 fOutputList(obj.fOutputList),
399 fPIDResponse(obj.fPIDResponse),
400 fEC(obj.fEC),
401 fEvt(obj.fEvt),
402 fTempStruct(obj.fTempStruct),
403 fRandomNumber(obj.fRandomNumber),
404 fLEGO(obj.fLEGO),
405 fMCcase(obj.fMCcase),
406 fAODcase(obj.fAODcase),
407 fPbPbcase(obj.fPbPbcase),
408 fGenerateSignal(obj.fGenerateSignal),
409 fGeneratorOnly(obj.fGeneratorOnly),
410 fPdensityPairCut(obj.fPdensityPairCut),
411 fTabulatePairs(obj.fTabulatePairs),
412 fRMax(obj.fRMax),
413 ffcSq(obj.ffcSq),
414 fFilterBit(obj.fFilterBit),
415 fMaxChi2NDF(obj.fMaxChi2NDF),
416 fMinTPCncls(obj.fMinTPCncls),
417 fBfield(obj.fBfield),
418 fMbin(obj.fMbin),
419 fFSIindex(obj.fFSIindex),
420 fEDbin(obj.fEDbin),
421 fMbins(obj.fMbins),
422 fMultLimit(obj.fMultLimit),
423 fCentBinLowLimit(obj.fCentBinLowLimit),
424 fCentBinHighLimit(obj.fCentBinHighLimit),
425 fEventCounter(obj.fEventCounter),
426 fEventsToMix(obj.fEventsToMix),
427 fZvertexBins(obj.fZvertexBins),
428 fMultLimits(),
429 fQcut(0),
430 fQLowerCut(obj.fQLowerCut),
431 fNormQcutLow(0),
432 fNormQcutHigh(0),
433 fKupperBound(obj.fKupperBound),
434 fQupperBoundQ2(obj.fQupperBoundQ2),
435 fQupperBoundQ3(obj.fQupperBoundQ3),
436 fQupperBoundQ4(obj.fQupperBoundQ4),
437 fQbinsQ2(obj.fQbinsQ2),
438 fQbinsQ3(obj.fQbinsQ3),
439 fQbinsQ4(obj.fQbinsQ4),
440 fQupperBoundWeights(obj.fQupperBoundWeights),
441 fKstepT(),
442 fKstepY(),
443 fKmeanT(),
444 fKmeanY(),
445 fKmiddleT(),
446 fKmiddleY(),
447 fQstep(obj.fQstep),
448 fQstepWeights(obj.fQstepWeights),
449 fQmean(),
450 fDampStart(obj.fDampStart),
451 fDampStep(obj.fDampStep),
452 fTPCTOFboundry(obj.fTPCTOFboundry),
453 fTOFboundry(obj.fTOFboundry),
454 fSigmaCutTPC(obj.fSigmaCutTPC),
455 fSigmaCutTOF(obj.fSigmaCutTOF),
456 fMinSepPairEta(obj.fMinSepPairEta),
457 fMinSepPairPhi(obj.fMinSepPairPhi),
458 fShareQuality(obj.fShareQuality),
459 fShareFraction(obj.fShareFraction),
460 fTrueMassP(obj.fTrueMassP),
461 fTrueMassPi(obj.fTrueMassPi),
462 fTrueMassK(obj.fTrueMassK),
463 fTrueMassKs(obj.fTrueMassKs),
464 fTrueMassLam(obj.fTrueMassLam),
465 fKtIndexL(obj.fKtIndexL),
466 fKtIndexH(obj.fKtIndexH),
467 fQoIndexL(obj.fQoIndexL),
468 fQoIndexH(obj.fQoIndexH),
469 fQsIndexL(obj.fQsIndexL),
470 fQsIndexH(obj.fQsIndexH),
471 fQlIndexL(obj.fQlIndexL),
472 fQlIndexH(obj.fQlIndexH),
473 fDummyB(obj.fDummyB),
474 fKT3transition(obj.fKT3transition),
475 fKT4transition(obj.fKT4transition),
b71263d0 476 fDefaultsCharSwitch(),
477 fLowQPairSwitch_E0E0(),
478 fLowQPairSwitch_E0E1(),
479 fLowQPairSwitch_E0E2(),
480 fLowQPairSwitch_E0E3(),
6bb6954b 481 fLowQPairSwitch_E1E1(),
b71263d0 482 fLowQPairSwitch_E1E2(),
483 fLowQPairSwitch_E1E3(),
484 fLowQPairSwitch_E2E3(),
485 fNormQPairSwitch_E0E0(),
486 fNormQPairSwitch_E0E1(),
487 fNormQPairSwitch_E0E2(),
488 fNormQPairSwitch_E0E3(),
6bb6954b 489 fNormQPairSwitch_E1E1(),
b71263d0 490 fNormQPairSwitch_E1E2(),
491 fNormQPairSwitch_E1E3(),
492 fNormQPairSwitch_E2E3(),
5591748e 493 fMomResC2(obj.fMomResC2),
494 fWeightmuonCorrection(obj.fWeightmuonCorrection)
be9ef9f9 495{
496 // Copy Constructor
497
498 for(Int_t i=0; i<12; i++){
499 fFSIss[i]=obj.fFSIss[i];
500 fFSIos[i]=obj.fFSIos[i];
501 }
502
503 // Initialize fNormWeight and fNormWeightErr to 0
504 for(Int_t i=0; i<3; i++){// Kt iterator
505 for(Int_t j=0; j<10; j++){// Mbin iterator
506 fNormWeight[i][j]=0x0;
507 }
508 }
509
510
511}
512//________________________________________________________________________
513AliFourPion &AliFourPion::operator=(const AliFourPion &obj)
514{
515 // Assignment operator
516 if (this == &obj)
517 return *this;
518
519 fname = obj.fname;
520 fAOD = obj.fAOD;
521 fOutputList = obj.fOutputList;
522 fPIDResponse = obj.fPIDResponse;
523 fEC = obj.fEC;
524 fEvt = obj.fEvt;
525 fTempStruct = obj.fTempStruct;
526 fRandomNumber = obj.fRandomNumber;
527 fLEGO = fLEGO;
528 fMCcase = obj.fMCcase;
529 fAODcase = obj.fAODcase;
530 fPbPbcase = obj.fPbPbcase;
531 fGenerateSignal = obj.fGenerateSignal;
532 fGeneratorOnly = obj.fGeneratorOnly;
533 fPdensityPairCut = obj.fPdensityPairCut;
534 fTabulatePairs = obj.fTabulatePairs;
535 fRMax = obj.fRMax;
536 ffcSq = obj.ffcSq;
537 fFilterBit = obj.fFilterBit;
538 fMaxChi2NDF = obj.fMaxChi2NDF;
539 fMinTPCncls = obj.fMinTPCncls;
540 fBfield = obj.fBfield;
541 fMbin = obj.fMbin;
542 fFSIindex = obj.fFSIindex;
543 fEDbin = obj.fEDbin;
544 fMbins = obj.fMbins;
545 fMultLimit = obj.fMultLimit;
546 fCentBinLowLimit = obj.fCentBinLowLimit;
547 fCentBinHighLimit = obj.fCentBinHighLimit;
548 fEventCounter = obj.fEventCounter;
549 fEventsToMix = obj.fEventsToMix;
550 fZvertexBins = obj.fZvertexBins;
551 fQLowerCut = obj.fQLowerCut;
552 fKupperBound = obj.fKupperBound;
553 fQupperBoundQ2 = obj.fQupperBoundQ2;
554 fQupperBoundQ3 = obj.fQupperBoundQ3;
555 fQupperBoundQ4 = obj.fQupperBoundQ4;
556 fQbinsQ2 = obj.fQbinsQ2;
557 fQbinsQ3 = obj.fQbinsQ3;
558 fQbinsQ4 = obj.fQbinsQ4;
559 fQupperBoundWeights = obj.fQupperBoundWeights;
560 fQstep = obj.fQstep;
561 fQstepWeights = obj.fQstepWeights;
562 fDampStart = obj.fDampStart;
563 fDampStep = obj.fDampStep;
564 fTPCTOFboundry = obj.fTPCTOFboundry;
565 fTOFboundry = obj.fTOFboundry;
566 fSigmaCutTPC = obj.fSigmaCutTPC;
567 fSigmaCutTOF = obj.fSigmaCutTOF;
568 fMinSepPairEta = obj.fMinSepPairEta;
569 fMinSepPairPhi = obj.fMinSepPairPhi;
570 fShareQuality = obj.fShareQuality;
571 fShareFraction = obj.fShareFraction;
572 fTrueMassP = obj.fTrueMassP;
573 fTrueMassPi = obj.fTrueMassPi;
574 fTrueMassK = obj.fTrueMassK;
575 fTrueMassKs = obj.fTrueMassKs;
576 fTrueMassLam = obj.fTrueMassLam;
577 fKtIndexL = obj.fKtIndexL;
578 fKtIndexH = obj.fKtIndexH;
579 fQoIndexL = obj.fQoIndexL;
580 fQoIndexH = obj.fQoIndexH;
581 fQsIndexL = obj.fQsIndexL;
582 fQsIndexH = obj.fQsIndexH;
583 fQlIndexL = obj.fQlIndexL;
584 fQlIndexH = obj.fQlIndexH;
585 fDummyB = obj.fDummyB;
586 fKT3transition = obj.fKT3transition;
587 fKT4transition = obj.fKT4transition;
588 fMomResC2 = obj.fMomResC2;
5591748e 589 fWeightmuonCorrection = obj.fWeightmuonCorrection;
be9ef9f9 590
591 for(Int_t i=0; i<12; i++){
592 fFSIss[i]=obj.fFSIss[i];
593 fFSIos[i]=obj.fFSIos[i];
594 }
595 for(Int_t i=0; i<3; i++){// Kt iterator
596 for(Int_t j=0; j<10; j++){// Mbin iterator
597 fNormWeight[i][j]=obj.fNormWeight[i][j];
598 }
599 }
600
601 return (*this);
602}
603//________________________________________________________________________
604AliFourPion::~AliFourPion()
605{
606 // Destructor
607 if(fAOD) delete fAOD;
608 //if(fESD) delete fESD;
609 if(fOutputList) delete fOutputList;
610 if(fPIDResponse) delete fPIDResponse;
611 if(fEC) delete fEC;
612 if(fEvt) delete fEvt;
613 if(fTempStruct) delete [] fTempStruct;
614 if(fRandomNumber) delete fRandomNumber;
615 if(fMomResC2) delete fMomResC2;
5591748e 616 if(fWeightmuonCorrection) delete fWeightmuonCorrection;
617
b71263d0 618 for(Int_t j=0; j<kMultLimitPbPb; j++){
619 if(fLowQPairSwitch_E0E0[j]) delete [] fLowQPairSwitch_E0E0[j];
620 if(fLowQPairSwitch_E0E1[j]) delete [] fLowQPairSwitch_E0E1[j];
621 if(fLowQPairSwitch_E0E2[j]) delete [] fLowQPairSwitch_E0E2[j];
622 if(fLowQPairSwitch_E0E3[j]) delete [] fLowQPairSwitch_E0E3[j];
6bb6954b 623 if(fLowQPairSwitch_E1E1[j]) delete [] fLowQPairSwitch_E1E1[j];
b71263d0 624 if(fLowQPairSwitch_E1E2[j]) delete [] fLowQPairSwitch_E1E2[j];
625 if(fLowQPairSwitch_E1E3[j]) delete [] fLowQPairSwitch_E1E3[j];
626 if(fLowQPairSwitch_E2E3[j]) delete [] fLowQPairSwitch_E2E3[j];
627 //
628 if(fNormQPairSwitch_E0E0[j]) delete [] fNormQPairSwitch_E0E0[j];
629 if(fNormQPairSwitch_E0E1[j]) delete [] fNormQPairSwitch_E0E1[j];
630 if(fNormQPairSwitch_E0E2[j]) delete [] fNormQPairSwitch_E0E2[j];
631 if(fNormQPairSwitch_E0E3[j]) delete [] fNormQPairSwitch_E0E3[j];
6bb6954b 632 if(fNormQPairSwitch_E1E1[j]) delete [] fNormQPairSwitch_E1E1[j];
b71263d0 633 if(fNormQPairSwitch_E1E2[j]) delete [] fNormQPairSwitch_E1E2[j];
634 if(fNormQPairSwitch_E1E3[j]) delete [] fNormQPairSwitch_E1E3[j];
635 if(fNormQPairSwitch_E2E3[j]) delete [] fNormQPairSwitch_E2E3[j];
636 }
be9ef9f9 637
638 //
639 for(Int_t mb=0; mb<fMbins; mb++){
640 for(Int_t edB=0; edB<fEDbins; edB++){
641 for(Int_t c1=0; c1<2; c1++){
642 for(Int_t c2=0; c2<2; c2++){
643 for(Int_t term=0; term<2; term++){
644
645 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2;
646
647 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal;
648 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared;
649 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL;
650 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW;
651 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL;
652 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW;
653 //
654 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv;
655 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW;
656 }// term_2
657
658 for(Int_t c3=0; c3<2; c3++){
659 for(Int_t term=0; term<5; term++){
660
661 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3;
662 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3;
b71263d0 663 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor;
be9ef9f9 664 //
665 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm;
666
667 }// term_3
b71263d0 668
669 for(Int_t c4=0; c4<2; c4++){
6bb6954b 670 for(Int_t term=0; term<13; term++){
b71263d0 671
672 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4;
673 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4;
674 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor;
675 //
676 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm;
677 }// term_4
678
679 }//c4
be9ef9f9 680 }//c3
681 }//c2
682 }//c1
683 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
684 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
685 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD;
686 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD;
687 }
688 }
689
690 }// ED
691 }// Mbin
692
693
694 for(Int_t i=0; i<12; i++){
695 if(fFSIss[i]) delete fFSIss[i];
696 if(fFSIos[i]) delete fFSIos[i];
697 }
698 for(Int_t i=0; i<3; i++){// Kt iterator
699 for(Int_t j=0; j<10; j++){// Mbin iterator
700 if(fNormWeight[i][j]) delete fNormWeight[i][j];
701 }
702 }
703
704}
705//________________________________________________________________________
706void AliFourPion::ParInit()
707{
708 cout<<"AliFourPion MyInit() call"<<endl;
709 cout<<"lego:"<<fLEGO<<" MCcase:"<<fMCcase<<" PbPbcase:"<<fPbPbcase<<" TabulatePairs:"<<fTabulatePairs<<" GenSignal:"<<fGenerateSignal<<" CentLow:"<<fCentBinLowLimit<<" CentHigh:"<<fCentBinHighLimit<<" RMax:"<<fRMax<<" fc^2:"<<ffcSq<<" FB:"<<fFilterBit<<" MaxChi2/NDF:"<<fMaxChi2NDF<<" MinTPCncls:"<<fMinTPCncls<<" MinPairSepEta:"<<fMinSepPairEta<<" MinPairSepPhi:"<<fMinSepPairPhi<<" NsigTPC:"<<fSigmaCutTPC<<" NsigTOF:"<<fSigmaCutTOF<<endl;
710
711 fRandomNumber = new TRandom3();
712 fRandomNumber->SetSeed(0);
713
714 //
715 fEventCounter=0;
716 fEventsToMix=3;
717 fZvertexBins=2;//2
718
719 fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
720 fTOFboundry = 2.1;// TOF pid used below this momentum
721
722 ////////////////////////////////////////////////
723 // PadRow Pair Cuts
724 fShareQuality = .5;// max
725 fShareFraction = .05;// max
726 ////////////////////////////////////////////////
727
728
729 fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
730 fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
731
732
733
734 if(fPbPbcase) {// PbPb
735 fMultLimit=kMultLimitPbPb;
736 fMbins=fCentBins;
737 fQcut=0.1;
b71263d0 738 fNormQcutLow = 0.15;// 0.15
739 fNormQcutHigh = 0.2;// 0.175
be9ef9f9 740 }else {// pp
741 fMultLimit=kMultLimitpp;
742 fMbins=kMultBinspp;
743 fQcut=0.6;
744 fNormQcutLow = 1.0;
745 fNormQcutHigh = 1.5;
746 }
747
748 fQLowerCut = 0.005;// was 0.005
749 fKupperBound = 1.0;
750 //
751 fKstepY[0] = 1.6;
752 fKmeanY[0] = 0;// central y
753 fKmiddleY[0] = 0;
754
755 // 4x1 (Kt: 0-0.25, 0.25-0.35, 0.35-0.45, 0.45-1.0)
756 if(fKbinsT==4){
757 fKstepT[0] = 0.25; fKstepT[1] = 0.1; fKstepT[2] = 0.1; fKstepT[3] = 0.55;
758 fKmeanT[0] = 0.212; fKmeanT[1] = 0.299; fKmeanT[2] = 0.398; fKmeanT[3] = 0.576;
759 fKmiddleT[0] = 0.125; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.4; fKmiddleT[3] = 0.725;
760 }
761 // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
762 if(fKbinsT==3){
763 fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
764 fKmeanT[0] = 0.240; fKmeanT[1] = 0.369; fKmeanT[2] = 0.576;
765 fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
766 }
767 // 2x1 (Kt: 0-0.35, 0.35-1.0)
768 if(fKbinsT==2){
769 fKstepT[0] = 0.35; fKstepT[1] = 0.65;
770 fKmeanT[0] = 0.264; fKmeanT[1] = 0.500;
771 fKmiddleT[0] = 0.175; fKmiddleT[1] = 0.675;
772 }
773
774 //
775 fQupperBoundWeights = 0.2;
776 fQupperBoundQ2 = 2.0;
5591748e 777 fQupperBoundQ3 = 0.6;
778 fQupperBoundQ4 = 0.6;
be9ef9f9 779 fQbinsQ2 = fQupperBoundQ2/0.005;
780 fQbinsQ3 = fQupperBoundQ3/0.005;
781 fQbinsQ4 = fQupperBoundQ4/0.005;
782 fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
783 for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
784 //
785 fDampStart = 0.5;// was 0.3, then 0.5
786 fDampStep = 0.02;
787
788 //
789 fKT3transition = 0.3;
790 fKT4transition = 0.3;
791
792 fEC = new AliFourPionEventCollection **[fZvertexBins];
793 for(UShort_t i=0; i<fZvertexBins; i++){
794
795 fEC[i] = new AliFourPionEventCollection *[fMbins];
796
797 for(UShort_t j=0; j<fMbins; j++){
798
799 fEC[i][j] = new AliFourPionEventCollection(fEventsToMix+1, fMultLimit, kMCarrayLimit, fMCcase);
800 }
801 }
802
b71263d0 803 for(Int_t i=0; i<kMultLimitPbPb; i++) fDefaultsCharSwitch[i]='0';
804 for(Int_t i=0; i<kMultLimitPbPb; i++) {
805 fLowQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
806 fLowQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
807 fLowQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
808 fLowQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 809 fLowQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 810 fLowQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
811 fLowQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
812 fLowQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
813 //
814 fNormQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
815 fNormQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
816 fNormQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
817 fNormQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 818 fNormQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 819 fNormQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
820 fNormQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
821 fNormQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
822 }
be9ef9f9 823
824 fTempStruct = new AliFourPionTrackStruct[fMultLimit];
b71263d0 825
826
be9ef9f9 827 fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
828
829
830
831 // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
832 if(!fLEGO) {
833 SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
834 if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
835 if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
5591748e 836 if(!fMCcase && !fTabulatePairs) SetMuonCorrections(fLEGO);// Read Muon corrections
be9ef9f9 837 }
838
839 /////////////////////////////////////////////
840 /////////////////////////////////////////////
841
842}
843//________________________________________________________________________
844void AliFourPion::UserCreateOutputObjects()
845{
846 // Create histograms
847 // Called once
848
849 ParInit();// Initialize my settings
850
851
852 fOutputList = new TList();
853 fOutputList->SetOwner();
854
855 TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
856 fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
857 fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
858 fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
859 fOutputList->Add(fVertexDist);
860
861
862 TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
863 fOutputList->Add(fDCAxyDistPlus);
864 TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
865 fOutputList->Add(fDCAzDistPlus);
866 TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
867 fOutputList->Add(fDCAxyDistMinus);
868 TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
869 fOutputList->Add(fDCAzDistMinus);
870
871
872 TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
873 fOutputList->Add(fEvents1);
874 TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
875 fOutputList->Add(fEvents2);
876
877 TH1F *fMultDist0 = new TH1F("fMultDist0","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
878 fMultDist0->GetXaxis()->SetTitle("Multiplicity");
879 fOutputList->Add(fMultDist0);
880
881 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
882 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
883 fOutputList->Add(fMultDist1);
884
885 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
886 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
887 fOutputList->Add(fMultDist2);
888
889 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
890 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
891 fOutputList->Add(fMultDist3);
892
893 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
894 fOutputList->Add(fPtEtaDist);
895
896 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
897 fOutputList->Add(fPhiPtDist);
898
899 TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
900 fOutputList->Add(fTOFResponse);
901 TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
902 fOutputList->Add(fTPCResponse);
903
904 TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",200,0,2);
905 fOutputList->Add(fRejectedPairs);
906 TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
907 fOutputList->Add(fRejectedEvents);
908
909 TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
910 if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
911 TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
912 if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
913 TH3F *fPairsShareFracDPhiNum = new TH3F("fPairsShareFracDPhiNum","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
914 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiNum);
915 TH3F *fPairsShareFracDPhiDen = new TH3F("fPairsShareFracDPhiDen","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
916 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiDen);
917 TH3D* fPairsPadRowNum = new TH3D("fPairsPadRowNum","", 20,0,1, 159,0,1, 40,0,0.2);
918 if(fMCcase) fOutputList->Add(fPairsPadRowNum);
919 TH3D* fPairsPadRowDen = new TH3D("fPairsPadRowDen","", 20,0,1, 159,0,1, 40,0,0.2);
920 if(fMCcase) fOutputList->Add(fPairsPadRowDen);
921
922
923
924 TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
925 if(fMCcase) fOutputList->Add(fResonanceOSPairs);
926 TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
927 if(fMCcase) fOutputList->Add(fAllOSPairs);
928
929 TH3D *fPrimarySCPionPairs = new TH3D("fPrimarySCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
930 if(fMCcase) fOutputList->Add(fPrimarySCPionPairs);
931 TH3D *fAllSCPionPairs = new TH3D("fAllSCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
932 if(fMCcase) fOutputList->Add(fAllSCPionPairs);
933 TH3D *fPrimaryMCPionPairs = new TH3D("fPrimaryMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
934 if(fMCcase) fOutputList->Add(fPrimaryMCPionPairs);
935 TH3D *fAllMCPionPairs = new TH3D("fAllMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
936 if(fMCcase) fOutputList->Add(fAllMCPionPairs);
937
938 //
939 TH1D *fMuonParents = new TH1D("fMuonParents","",500,0.5,500.5);
940 if(fMCcase) fOutputList->Add(fMuonParents);
941 TH1D *fSecondaryMuonParents = new TH1D("fSecondaryMuonParents","",500,0.5,500.5);
942 if(fMCcase) fOutputList->Add(fSecondaryMuonParents);
943 TH3D *fMuonPionDeltaQinv = new TH3D("fMuonPionDeltaQinv","",2,-0.5,1.5, 20,0,1, 100,-0.2,0.2);
944 if(fMCcase) fOutputList->Add(fMuonPionDeltaQinv);
945 TH1D *fPionCandidates = new TH1D("fPionCandidates","",500,0.5,500.5);
946 if(fMCcase) fOutputList->Add(fPionCandidates);
947 //
948
949
950 TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
951 fOutputList->Add(fAvgMult);
952
953 TH2D *fTrackChi2NDF = new TH2D("fTrackChi2NDF","",20,0,100, 100,0,10);
954 fOutputList->Add(fTrackChi2NDF);
955 TH2D *fTrackTPCncls = new TH2D("fTrackTPCncls","",20,0,100, 110,50,160);
956 fOutputList->Add(fTrackTPCncls);
957
958
959 TH1D *fTPNRejects3pion1 = new TH1D("fTPNRejects3pion1","",fQbinsQ3,0,fQupperBoundQ3);
960 fOutputList->Add(fTPNRejects3pion1);
961 TH1D *fTPNRejects3pion2 = new TH1D("fTPNRejects3pion2","",fQbinsQ3,0,fQupperBoundQ3);
962 fOutputList->Add(fTPNRejects3pion2);
963 TH1D *fTPNRejects4pion1 = new TH1D("fTPNRejects4pion1","",fQbinsQ4,0,fQupperBoundQ4);
964 fOutputList->Add(fTPNRejects4pion1);
965
966 TH3D *fKT3DistTerm1 = new TH3D("fKT3DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
967 TH3D *fKT3DistTerm5 = new TH3D("fKT3DistTerm5","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
968 fOutputList->Add(fKT3DistTerm1);
969 fOutputList->Add(fKT3DistTerm5);
970 TH3D *fKT4DistTerm1 = new TH3D("fKT4DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
6bb6954b 971 TH3D *fKT4DistTerm13 = new TH3D("fKT4DistTerm13","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
be9ef9f9 972 fOutputList->Add(fKT4DistTerm1);
6bb6954b 973 fOutputList->Add(fKT4DistTerm13);
be9ef9f9 974
975
976 TProfile2D *fKT3AvgpT = new TProfile2D("fKT3AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
977 fOutputList->Add(fKT3AvgpT);
978 TProfile2D *fKT4AvgpT = new TProfile2D("fKT4AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
979 fOutputList->Add(fKT4AvgpT);
980
981
982 TH1D *fMCWeight3DTerm1SC = new TH1D("fMCWeight3DTerm1SC","", 20,0,0.2);
983 TH1D *fMCWeight3DTerm1SCden = new TH1D("fMCWeight3DTerm1SCden","", 20,0,0.2);
984 TH1D *fMCWeight3DTerm2SC = new TH1D("fMCWeight3DTerm2SC","", 20,0,0.2);
985 TH1D *fMCWeight3DTerm2SCden = new TH1D("fMCWeight3DTerm2SCden","", 20,0,0.2);
986 TH1D *fMCWeight3DTerm1MC = new TH1D("fMCWeight3DTerm1MC","", 20,0,0.2);
987 TH1D *fMCWeight3DTerm1MCden = new TH1D("fMCWeight3DTerm1MCden","", 20,0,0.2);
988 TH1D *fMCWeight3DTerm2MC = new TH1D("fMCWeight3DTerm2MC","", 20,0,0.2);
989 TH1D *fMCWeight3DTerm2MCden = new TH1D("fMCWeight3DTerm2MCden","", 20,0,0.2);
990 TH1D *fMCWeight3DTerm3MC = new TH1D("fMCWeight3DTerm3MC","", 20,0,0.2);
991 TH1D *fMCWeight3DTerm3MCden = new TH1D("fMCWeight3DTerm3MCden","", 20,0,0.2);
992 TH1D *fMCWeight3DTerm4MC = new TH1D("fMCWeight3DTerm4MC","", 20,0,0.2);
993 TH1D *fMCWeight3DTerm4MCden = new TH1D("fMCWeight3DTerm4MCden","", 20,0,0.2);
994 fOutputList->Add(fMCWeight3DTerm1SC);
995 fOutputList->Add(fMCWeight3DTerm1SCden);
996 fOutputList->Add(fMCWeight3DTerm2SC);
997 fOutputList->Add(fMCWeight3DTerm2SCden);
998 fOutputList->Add(fMCWeight3DTerm1MC);
999 fOutputList->Add(fMCWeight3DTerm1MCden);
1000 fOutputList->Add(fMCWeight3DTerm2MC);
1001 fOutputList->Add(fMCWeight3DTerm2MCden);
1002 fOutputList->Add(fMCWeight3DTerm3MC);
1003 fOutputList->Add(fMCWeight3DTerm3MCden);
1004 fOutputList->Add(fMCWeight3DTerm4MC);
1005 fOutputList->Add(fMCWeight3DTerm4MCden);
1006
1007
1008 if(fPdensityPairCut){
1009
1010 for(Int_t mb=0; mb<fMbins; mb++){
1011 if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
1012
1013 for(Int_t edB=0; edB<fEDbins; edB++){
1014 for(Int_t c1=0; c1<2; c1++){
1015 for(Int_t c2=0; c2<2; c2++){
1016 for(Int_t term=0; term<2; term++){
1017
1018 TString *nameEx2 = new TString("TwoParticle_Charge1_");
1019 *nameEx2 += c1;
1020 nameEx2->Append("_Charge2_");
1021 *nameEx2 += c2;
1022 nameEx2->Append("_M_");
1023 *nameEx2 += mb;
1024 nameEx2->Append("_ED_");
1025 *nameEx2 += edB;
1026 nameEx2->Append("_Term_");
1027 *nameEx2 += term+1;
1028
1029 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
1030
1031
1032 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1033 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2);
1034 TString *nameEx2QW=new TString(nameEx2->Data());
1035 nameEx2QW->Append("_QW");
1036 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1037 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW);
1038 TString *nameAvgP=new TString(nameEx2->Data());
1039 nameAvgP->Append("_AvgP");
1040 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP = new TProfile2D(nameAvgP->Data(),"",10,0,1, fQbinsQ2,0,fQupperBoundQ2, 0.,1.0,"");
1041 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP);
1042
1043 TString *nameUnitMult=new TString(nameEx2->Data());
1044 nameUnitMult->Append("_UnitMult");
1045 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin = new TH2D(nameUnitMult->Data(),"Two Particle Distribution",21,0.5,21.5, fQbinsQ2,0,fQupperBoundQ2);
1046 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin);
1047
1048 if(fMCcase){
1049 // Momentum resolution histos
1050 TString *nameIdeal = new TString(nameEx2->Data());
1051 nameIdeal->Append("_Ideal");
1052 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1053 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal);
1054 TString *nameSmeared = new TString(nameEx2->Data());
1055 nameSmeared->Append("_Smeared");
1056 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1057 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared);
1058 //
1059 // Muon correction histos
1060 TString *nameMuonIdeal=new TString(nameEx2->Data());
1061 nameMuonIdeal->Append("_MuonIdeal");
1062 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1063 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal);
1064 TString *nameMuonSmeared=new TString(nameEx2->Data());
1065 nameMuonSmeared->Append("_MuonSmeared");
1066 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1067 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared);
1068 //
1069 TString *nameMuonPionK2=new TString(nameEx2->Data());
1070 nameMuonPionK2->Append("_MuonPionK2");
1071 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2 = new TH2D(nameMuonPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1072 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2);
1073 //
1074 TString *namePionPionK2=new TString(nameEx2->Data());
1075 namePionPionK2->Append("_PionPionK2");
1076 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2 = new TH2D(namePionPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1077 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2);
1078 //
1079 //
1080 TString *nameEx2MC=new TString(nameEx2->Data());
1081 nameEx2MC->Append("_MCqinv");
1082 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1083 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
1084 TString *nameEx2MCQW=new TString(nameEx2->Data());
1085 nameEx2MCQW->Append("_MCqinvQW");
1086 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1087 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
1088 //
1089 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
1090 nameEx2PIDpurityDen->Append("_PIDpurityDen");
1091 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen = new TH2D(nameEx2PIDpurityDen->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1092 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
1093 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
1094 nameEx2PIDpurityNum->Append("_PIDpurityNum");
1095 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH3D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",16,0.5,16.5, 20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1096 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
1097 }
1098 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
1099 nameEx2OSLB1->Append("_osl_b1");
1100 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1101 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL);
1102 nameEx2OSLB1->Append("_QW");
1103 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1104 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW);
1105 //
1106 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
1107 nameEx2OSLB2->Append("_osl_b2");
1108 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1109 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL);
1110 nameEx2OSLB2->Append("_QW");
1111 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1112 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW);
1113
1114 }// term_2
1115
1116
1117
1118 // skip 3-particle if Tabulate6DPairs is true
1119 if(fTabulatePairs) continue;
1120
1121 for(Int_t c3=0; c3<2; c3++){
1122 for(Int_t term=0; term<5; term++){
1123
1124 TString *namePC3 = new TString("ThreeParticle_Charge1_");
1125 *namePC3 += c1;
1126 namePC3->Append("_Charge2_");
1127 *namePC3 += c2;
1128 namePC3->Append("_Charge3_");
1129 *namePC3 += c3;
1130 namePC3->Append("_M_");
1131 *namePC3 += mb;
1132 namePC3->Append("_ED_");
1133 *namePC3 += edB;
1134 namePC3->Append("_Term_");
1135 *namePC3 += term+1;
1136
1137 ///////////////////////////////////////
1138 // skip degenerate histograms
1139 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
1140 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
1141 /////////////////////////////////////////
1142
1143
1144 TString *nameNorm=new TString(namePC3->Data());
1145 nameNorm->Append("_Norm");
1146 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1147 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1148 //
1149
1150 TString *name1DQ=new TString(namePC3->Data());
1151 name1DQ->Append("_1D");
1152 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH1D(name1DQ->Data(),"", fQbinsQ3,0,fQupperBoundQ3);
1153 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1154 //
1155 TString *nameKfactor=new TString(namePC3->Data());
1156 nameKfactor->Append("_Kfactor");
1157 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,100, "");
1158 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor);
1159 //
1160 TString *nameMeanQinv=new TString(namePC3->Data());
1161 nameMeanQinv->Append("_MeanQinv");
1162 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv = new TProfile(nameMeanQinv->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,.2, "");
1163 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv);
1164
1165 if(fMCcase==kTRUE){
1166 // Momentum resolution correction histos
1167 TString *nameMomResIdeal=new TString(namePC3->Data());
1168 nameMomResIdeal->Append("_Ideal");
1169 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1170 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1171 TString *nameMomResSmeared=new TString(namePC3->Data());
1172 nameMomResSmeared->Append("_Smeared");
1173 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1174 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared);
1175 // Muon correction histos
1176 TString *nameMuonIdeal=new TString(namePC3->Data());
1177 nameMuonIdeal->Append("_MuonIdeal");
5591748e 1178 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1179 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal);
be9ef9f9 1180 TString *nameMuonSmeared=new TString(namePC3->Data());
1181 nameMuonSmeared->Append("_MuonSmeared");
5591748e 1182 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1183 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared);
be9ef9f9 1184 //
1185 TString *nameMuonPionK3=new TString(namePC3->Data());
1186 nameMuonPionK3->Append("_MuonPionK3");
5591748e 1187 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3 = new TH3D(nameMuonPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1188 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3);
be9ef9f9 1189 //
1190 TString *namePionPionK3=new TString(namePC3->Data());
1191 namePionPionK3->Append("_PionPionK3");
5591748e 1192 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3 = new TH3D(namePionPionK3->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
1193 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3);
be9ef9f9 1194
1195 }// MCcase
1196 //
1197 if(c1==c2 && c1==c3 && term==4 ){
1198 TString *nameTwoPartNorm=new TString(namePC3->Data());
1199 nameTwoPartNorm->Append("_TwoPartNorm");
1200 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
1201 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm);
1202 }// term=4
1203
1204 }// term_3
1205
1206 for(Int_t c4=0; c4<2; c4++){
6bb6954b 1207 for(Int_t term=0; term<13; term++){
be9ef9f9 1208
1209 TString *namePC4 = new TString("FourParticle_Charge1_");
1210 *namePC4 += c1;
1211 namePC4->Append("_Charge2_");
1212 *namePC4 += c2;
1213 namePC4->Append("_Charge3_");
1214 *namePC4 += c3;
1215 namePC4->Append("_Charge4_");
1216 *namePC4 += c4;
1217 namePC4->Append("_M_");
1218 *namePC4 += mb;
1219 namePC4->Append("_ED_");
1220 *namePC4 += edB;
1221 namePC4->Append("_Term_");
1222 *namePC4 += term+1;
1223
1224 ///////////////////////////////////////
1225 // skip degenerate histograms
1226 if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
1227 if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
1228 if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
1229 /////////////////////////////////////////
1230
1231 TString *nameNorm=new TString(namePC4->Data());
1232 nameNorm->Append("_Norm");
1233 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1234 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4);
1235 //
1236 TString *name1DQ=new TString(namePC4->Data());
1237 name1DQ->Append("_1D");
1238 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = new TH1D(name1DQ->Data(),"", fQbinsQ4,0,fQupperBoundQ4);
1239 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4);
1240 //
1241 TString *nameKfactor=new TString(namePC4->Data());
1242 nameKfactor->Append("_Kfactor");
1243 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
1244 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor);
1245 //
6bb6954b 1246 if(c1==c2 && c1==c3 && c1==c4 && term==12 ){
be9ef9f9 1247 TString *nameTwoPartNorm=new TString(namePC4->Data());
1248 nameTwoPartNorm->Append("_TwoPartNorm");
1249 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
1250 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm);
1251 }
5591748e 1252
be9ef9f9 1253 if(fMCcase==kTRUE){
1254 // Momentum resolution correction histos
1255 TString *nameMomResIdeal=new TString(namePC4->Data());
1256 nameMomResIdeal->Append("_Ideal");
1257 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1258 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal);
1259 TString *nameMomResSmeared=new TString(namePC4->Data());
1260 nameMomResSmeared->Append("_Smeared");
1261 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1262 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared);
1263 // Muon correction histos
1264 TString *nameMuonIdeal=new TString(namePC4->Data());
1265 nameMuonIdeal->Append("_MuonIdeal");
5591748e 1266 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal = new TH3D(nameMuonIdeal->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1267 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal);
be9ef9f9 1268 TString *nameMuonSmeared=new TString(namePC4->Data());
1269 nameMuonSmeared->Append("_MuonSmeared");
5591748e 1270 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared = new TH3D(nameMuonSmeared->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1271 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared);
be9ef9f9 1272 //
1273 TString *nameMuonPionK4=new TString(namePC4->Data());
1274 nameMuonPionK4->Append("_MuonPionK4");
5591748e 1275 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4 = new TH3D(nameMuonPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1276 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4);
be9ef9f9 1277 //
1278 TString *namePionPionK4=new TString(namePC4->Data());
1279 namePionPionK4->Append("_PionPionK4");
5591748e 1280 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4 = new TH3D(namePionPionK4->Data(),"", 2,0.5,2.5, 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
1281 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4);
be9ef9f9 1282
1283 }// MCcase
1284
1285
1286 }
1287 }
1288
1289 }//c3
1290 }//c2
1291 }//c1
1292 }// ED
1293 }// mbin
1294 }// Pdensity Method
1295
1296
1297 if(fTabulatePairs){
1298
1299 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
1300 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
1301 for(Int_t mb=0; mb<fMbins; mb++){
1302 for(Int_t edB=0; edB<fEDbins; edB++){
1303
1304 TString *nameNum = new TString("TPN_num_Kt_");
1305 *nameNum += tKbin;
1306 nameNum->Append("_Ky_");
1307 *nameNum += yKbin;
1308 nameNum->Append("_M_");
1309 *nameNum += mb;
1310 nameNum->Append("_ED_");
1311 *nameNum += edB;
1312
1313 TString *nameDen = new TString("TPN_den_Kt_");
1314 *nameDen += tKbin;
1315 nameDen->Append("_Ky_");
1316 *nameDen += yKbin;
1317 nameDen->Append("_M_");
1318 *nameDen += mb;
1319 nameDen->Append("_ED_");
1320 *nameDen += edB;
1321
1322 if(edB==0){
1323 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = new TH3D(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1324 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD);
1325
1326 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = new TH3D(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1327 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD);
1328 }
1329
1330 }
1331 }
1332 }
1333 }
1334
1335 }
1336
1337
1338 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1339 fOutputList->Add(fQsmearMean);
1340 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1341 fOutputList->Add(fQsmearSq);
1342 TH2D *fQ2Res = new TH2D("fQ2Res","",20,0,1, 200,-.2,.2);
1343 fOutputList->Add(fQ2Res);
1344 TH2D *fQ3Res = new TH2D("fQ3Res","",20,0,1, 200,-.3,.3);
1345 fOutputList->Add(fQ3Res);
1346 TH2D *fQ4Res = new TH2D("fQ4Res","",20,0,1, 200,-.4,.4);
1347 fOutputList->Add(fQ4Res);
1348
1349 TH2D *DistQinv4pion = new TH2D("DistQinv4pion","",6,0.5,6.5, 20,0,0.1);
1350 fOutputList->Add(DistQinv4pion);
1351 TH2D *DistQinvMC4pion = new TH2D("DistQinvMC4pion","",6,0.5,6.5, 20,0,0.1);
1352 if(fMCcase) fOutputList->Add(DistQinvMC4pion);
1353
1354 TH2D *fAvgQ12VersusQ3 = new TH2D("fAvgQ12VersusQ3","",10,0,0.1, 20,0,0.1);
1355 fOutputList->Add(fAvgQ12VersusQ3);
1356 TH2D *fAvgQ13VersusQ3 = new TH2D("fAvgQ13VersusQ3","",10,0,0.1, 20,0,0.1);
1357 fOutputList->Add(fAvgQ13VersusQ3);
1358 TH2D *fAvgQ23VersusQ3 = new TH2D("fAvgQ23VersusQ3","",10,0,0.1, 20,0,0.1);
1359 fOutputList->Add(fAvgQ23VersusQ3);
1360
5591748e 1361 TH1D *fDistPionParents4 = new TH1D("fDistPionParents4","",4,0.5,4.5);
1362 fOutputList->Add(fDistPionParents4);
1363
be9ef9f9 1364 ////////////////////////////////////
1365 ///////////////////////////////////
1366
1367 PostData(1, fOutputList);
1368}
1369
1370//________________________________________________________________________
1371void AliFourPion::UserExec(Option_t *)
1372{
1373 // Main loop
1374 // Called for each event
1375 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
1376 fEventCounter++;
1377
1378 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1379
1380 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1381 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
1382
1383
1384 // Trigger Cut
1385 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1386 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1387 if(!isSelected1 && !fMCcase) {return;}
1388 }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1389 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1390 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1391 if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1392 }else {return;}
1393
1394 ///////////////////////////////////////////////////////////
1395 const AliAODVertex *primaryVertexAOD;
1396 AliCentrality *centrality;// for AODs and ESDs
1397
1398
1399 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1400 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1401 fPIDResponse = inputHandler->GetPIDResponse();
1402
1403
1404 TClonesArray *mcArray = 0x0;
1405 if(fMCcase){
1406 if(fAODcase){
1407 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1408 if(!mcArray || mcArray->GetEntriesFast() >= kMCarrayLimit){
1409 cout<<"No MC particle branch found or Array too large!!"<<endl;
1410 return;
1411 }
1412 }
1413 }
1414
1415
1416 UInt_t status=0;
1417 Int_t positiveTracks=0, negativeTracks=0;
1418 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1419
1420 Double_t vertex[3]={0};
1421 Int_t zbin=0;
1422 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1423 /////////////////////////////////////////////////
1424
1425
1426 Float_t centralityPercentile=0;
1427 Float_t cStep=5.0, cStart=0;
1428
1429
1430 if(fAODcase){// AOD case
1431
1432 if(fPbPbcase){
1433 centrality = fAOD->GetCentrality();
1434 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1435 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1436 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range. Skipping Event"<<endl;*/ return;}
1437 cout<<"Centrality % = "<<centralityPercentile<<endl;
1438 }
1439
b71263d0 1440
be9ef9f9 1441 ((TH1F*)fOutputList->FindObject("fMultDist0"))->Fill(fAOD->GetNumberOfTracks());
1442
1443 // Pile-up rejection
1444 AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1445 if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1446 else AnaUtil->SetUseMVPlpSelection(kFALSE);
1447 Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
1448 if(pileUpCase) return;
1449
1450 ////////////////////////////////
1451 // Vertexing
1452 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1453 primaryVertexAOD = fAOD->GetPrimaryVertex();
1454 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1455
1456 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1457 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1458
1459 if(!fMCcase && primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1460
1461 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1462
1463 fBfield = fAOD->GetMagneticField();
1464
1465 for(Int_t i=0; i<fZvertexBins; i++){
1466 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1467 zbin=i;
1468 break;
1469 }
1470 }
1471
1472
1473
1474 /////////////////////////////
1475 // Create Shuffled index list
1476 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1477 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1478 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1479 /////////////////////////////
1480
1481
1482 // Track loop
1483 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1484 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1485 if (!aodtrack) continue;
1486 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1487
1488 status=aodtrack->GetStatus();
1489
1490 if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1491 ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
1492 ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
1493 if(aodtrack->GetTPCNcls() < fMinTPCncls) continue;// TPC nCluster cut
1494 if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
1495
1496 if(fFilterBit != 7){
1497 Bool_t goodTrackOtherFB = kFALSE;
1498 for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1499 AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1500 if(!aodtrack2) continue;
1501 if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1502
1503 if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1504
1505 }
1506 if(!goodTrackOtherFB) continue;
1507 }
1508
1509
1510 if(aodtrack->Pt() < 0.16) continue;
1511 if(fabs(aodtrack->Eta()) > 0.8) continue;
1512
1513
1514 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1515 if(!goodMomentum) continue;
1516 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1517
1518
1519 Double_t dca2[2]={0};
1520 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1521 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1522 Double_t dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1523
1524 fTempStruct[myTracks].fStatus = status;
1525 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1526 fTempStruct[myTracks].fId = aodtrack->GetID();
1527 //
1528 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1529 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1530 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1531 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1532 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1533 fTempStruct[myTracks].fEta = aodtrack->Eta();
1534 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1535 fTempStruct[myTracks].fDCAXY = dca2[0];
1536 fTempStruct[myTracks].fDCAZ = dca2[1];
1537 fTempStruct[myTracks].fDCA = dca3d;
1538 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1539 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1540
1541
1542
1543 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1544 //if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1545 //if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1546
1547
1548 // PID section
1549 fTempStruct[myTracks].fElectron = kFALSE;
1550 fTempStruct[myTracks].fPion = kFALSE;
1551 fTempStruct[myTracks].fKaon = kFALSE;
1552 fTempStruct[myTracks].fProton = kFALSE;
1553
1554 Float_t nSigmaTPC[5];
1555 Float_t nSigmaTOF[5];
1556 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1557 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1558 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1559 Float_t signalTPC=0, signalTOF=0;
1560 Double_t integratedTimesTOF[10]={0};
1561
1562
1563 Bool_t DoPIDWorkAround=kTRUE;
1564 //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
1565 if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
1566 if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
1567 nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
1568 nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
1569 nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
1570 nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
1571 nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
1572 //
1573 nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
1574 nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
1575 nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
1576 nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
1577 nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
1578 signalTPC = aodtrack->GetTPCsignal();
1579 if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1580 fTempStruct[myTracks].fTOFhit = kTRUE;
1581 signalTOF = aodtrack->GetTOFsignal();
1582 aodtrack->GetIntegratedTimes(integratedTimesTOF);
1583 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1584
1585 }else {// FilterBit 7 PID workaround
1586
1587 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1588 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1589 if (!aodTrack2) continue;
1590 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1591
1592 UInt_t status2=aodTrack2->GetStatus();
1593
1594 nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
1595 nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
1596 nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
1597 nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
1598 nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
1599 //
1600 nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
1601 nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
1602 nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
1603 nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
1604 nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
1605 signalTPC = aodTrack2->GetTPCsignal();
1606
1607 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1608 fTempStruct[myTracks].fTOFhit = kTRUE;
1609 signalTOF = aodTrack2->GetTOFsignal();
1610 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1611 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1612
1613 }// aodTrack2
1614 }// FilterBit 7 PID workaround
1615
1616
1617 ///////////////////
1618 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1619 if(fTempStruct[myTracks].fTOFhit) {
1620 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1621 }
1622 ///////////////////
1623
1624 // Use TOF if good hit and above threshold
1625 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1626 if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1627 if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1628 if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1629 if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1630 }else {// TPC info instead
1631 if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1632 if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1633 if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1634 if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1635 }
1636
1637
1638 // Ensure there is only 1 candidate per track
1639 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1640 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1641 if(!fTempStruct[myTracks].fPion) continue;// only pions
1642 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1643 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1644 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1645 //if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;// superfluous
1646 ////////////////////////
1647 //if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons// superfluous
1648
1649
1650
1651 if(fTempStruct[myTracks].fCharge==+1) {
1652 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1653 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1654 }else {
1655 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1656 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1657 }
1658
1659 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1660 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1661
1662
1663
1664 if(fTempStruct[myTracks].fPion) {// pions
1665 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1666 fTempStruct[myTracks].fKey = 1;
1667 }else if(fTempStruct[myTracks].fKaon){// kaons
1668 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1669 fTempStruct[myTracks].fKey = 10;
1670 }else{// protons
1671 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1672 fTempStruct[myTracks].fKey = 100;
1673 }
1674
1675
1676
1677 if(aodtrack->Charge() > 0) positiveTracks++;
1678 else negativeTracks++;
1679
1680 if(fTempStruct[myTracks].fPion) pionCount++;
1681 if(fTempStruct[myTracks].fKaon) kaonCount++;
1682 if(fTempStruct[myTracks].fProton) protonCount++;
1683
1684 myTracks++;
1685
1686 if(fMCcase){// muon mothers
1687 AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
1688 if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
1689 AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
1690 if(parent->IsPhysicalPrimary()){
1691 ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
1692 }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
1693 }
1694 ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
1695 }
1696 }
1697 //cout<<"kinkcount = "<<kinkcount<<" pionkinks = "<<pionkinks<<" primarypionkinks = "<<primarypionkinks<<endl;
1698 }else {// ESD tracks
1699 cout<<"ESDs not supported currently"<<endl;
1700 return;
1701 }
1702
1703 // Generator info only
1704 if(fMCcase && fGeneratorOnly){
1705 myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
1706 for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
1707 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1708 if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
1709
1710 AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
1711 if(!mcParticle) continue;
1712 if(fabs(mcParticle->Eta())>0.8) continue;
1713 if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
1714 if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
1715 if(!mcParticle->IsPrimary()) continue;
1716 if(!mcParticle->IsPhysicalPrimary()) continue;
1717 if(abs(mcParticle->GetPdgCode())!=211) continue;
1718
1719 fTempStruct[myTracks].fP[0] = mcParticle->Px();
1720 fTempStruct[myTracks].fP[1] = mcParticle->Py();
1721 fTempStruct[myTracks].fP[2] = mcParticle->Pz();
1722 fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
1723
1724 fTempStruct[myTracks].fId = myTracks;// use my track counter
1725 fTempStruct[myTracks].fLabel = mctrackN;
1726 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1727 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1728 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1729 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1730 fTempStruct[myTracks].fEta = mcParticle->Eta();
1731 fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
1732 fTempStruct[myTracks].fDCAXY = 0.;
1733 fTempStruct[myTracks].fDCAZ = 0.;
1734 fTempStruct[myTracks].fDCA = 0.;
1735 fTempStruct[myTracks].fPion = kTRUE;
1736 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1737 fTempStruct[myTracks].fKey = 1;
1738
1739 myTracks++;
1740 pionCount++;
1741 }
1742 }
1743
1744 if(myTracks >= 1) {
1745 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1746 }
1747
1748
1749 //cout<<"There are "<<myTracks<<" myTracks"<<endl;
1750 //cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
1751
1752
1753 /////////////////////////////////////////
1754 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1755 if(myTracks < 4) {cout<<"Less than 4 tracks. Skipping Event."<<endl; return;}
1756 /////////////////////////////////////////
1757
1758
1759 ////////////////////////////////
1760 ///////////////////////////////
1761 // Mbin determination
1762 //
1763 // Mbin set to Pion Count Only for pp!!!!!!!
1764 fMbin=-1;
1765 if(!fPbPbcase){
1766 for(Int_t i=0; i<kMultBinspp; i++){
1767 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1768 // Mbin 0 has 1 pion
1769 }
1770 }else{
1771 for(Int_t i=0; i<fCentBins; i++){
1772 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1773 fMbin=i;// 0 = most central
1774 break;
1775 }
1776 }
1777 }
b71263d0 1778
be9ef9f9 1779 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1780
1781 ///////////////////
1782 // can only be called after fMbin has been set
1783 // Radius parameter only matters for Monte-Carlo data
1784 SetFSIindex(fRMax);
1785 ///////////////////
1786
1787 Int_t rBinForTPNMomRes = 10;
1788 if(fMbin==0) {rBinForTPNMomRes=10;}// 10 fm with EW (fRMax should be 11 for normal running)
1789 else if(fMbin==1) {rBinForTPNMomRes=9;}
1790 else if(fMbin<=3) {rBinForTPNMomRes=8;}
1791 else if(fMbin<=5) {rBinForTPNMomRes=7;}
1792 else {rBinForTPNMomRes=6;}
1793
1794 //////////////////////////////////////////////////
1795 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1796 //////////////////////////////////////////////////
1797
1798
1799
1800 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1801 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1802
1803 ////////////////////////////////////
1804 // Add event to buffer if > 0 tracks
1805 if(myTracks > 0){
1806 fEC[zbin][fMbin]->FIFOShift();
1807 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1808 (fEvt)->fNtracks = myTracks;
1809 (fEvt)->fFillStatus = 1;
1810 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1811 if(fMCcase){
1812 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1813 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1814 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1815 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1816 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1817 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1818 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1819 (fEvt)->fMCtracks[i].fPdgCode = tempMCTrack->GetPdgCode();
1820 (fEvt)->fMCtracks[i].fMotherLabel = tempMCTrack->GetMother();
1821 }
1822 }
1823 }
1824
1825
1826
1827 Float_t qinv12=0, qinv13=0, qinv14=0, qinv23=0, qinv24=0, qinv34=0;
1828 Float_t qout=0, qside=0, qlong=0;
1829 Float_t kT12=0;
1830 Float_t q3=0, q3MC=0;
1831 Float_t q4=0, q4MC=0;
1832 Int_t ch1=0, ch2=0, ch3=0, ch4=0;
1833 Int_t bin1=0, bin2=0, bin3=0, bin4=0;
1834 Float_t pVect1[4]={0};
1835 Float_t pVect2[4]={0};
1836 Float_t pVect3[4]={0};
1837 Float_t pVect4[4]={0};
1838 Float_t pVect1MC[4]={0};
1839 Float_t pVect2MC[4]={0};
1840 Float_t pVect3MC[4]={0};
1841 Float_t pVect4MC[4]={0};
1842 Float_t Pparent1[4]={0};
1843 Float_t Pparent2[4]={0};
1844 Float_t Pparent3[4]={0};
1845 Float_t Pparent4[4]={0};
1846 Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0;
1847 Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0;
5591748e 1848 Float_t weight12CC[3]={0};
1849 Float_t weight13CC[3]={0};
1850 Float_t weight14CC[3]={0};
1851 Float_t weight23CC[3]={0};
1852 Float_t weight24CC[3]={0};
1853 Float_t weight34CC[3]={0};
be9ef9f9 1854 Float_t weightTotal=0;//, weightTotalErr=0;
1855 Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0;
1856 Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0;
1857 Float_t parentQ3=0;
1858 Float_t FSICorr12=0, FSICorr13=0, FSICorr14=0, FSICorr23=0, FSICorr24=0, FSICorr34=0;
1859 Bool_t pionParent1=kFALSE, pionParent2=kFALSE, pionParent3=kFALSE, pionParent4=kFALSE;
1860 Bool_t FilledMCpair12=kFALSE, FilledMCtriplet123=kFALSE;
1861 Bool_t GoodTripletWeights=kFALSE;
1862 //
1863 AliAODMCParticle *mcParticle1=0x0;
1864 AliAODMCParticle *mcParticle2=0x0;
1865
1866
1867 ////////////////////
b71263d0 1868 //Int_t PairCount[7]={0};
1869 //Int_t NormPairCount[7]={0};
be9ef9f9 1870 Int_t KT3index=0, KT4index=0;
1871
b71263d0 1872 // reset to defaults
1873 for(Int_t i=0; i<kMultLimitPbPb; i++) {
1874 fLowQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1875 fLowQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1876 fLowQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1877 fLowQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 1878 fLowQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 1879 fLowQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1880 fLowQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1881 fLowQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1882 //
1883 fNormQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1884 fNormQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1885 fNormQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1886 fNormQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 1887 fNormQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 1888 fNormQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1889 fNormQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1890 fNormQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1891 }
5591748e 1892
b71263d0 1893
1894 //////////////////////////////////////////
1895 // make low-q pair storage and normalization-pair storage
1896 //
1897 for(Int_t en1=0; en1<=2; en1++){// 1st event number (en1=0 is the same event as current event)
1898 for(Int_t en2=en1; en2<=3; en2++){// 2nd event number (en2=0 is the same event as current event)
6bb6954b 1899 if(en1>1 && en1==en2) continue;
1900
b71263d0 1901 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {// 1st particle
1902 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
1903
5591748e 1904
b71263d0 1905 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1906 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1907 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1908 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1909 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
1910 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1911
1912 qinv12 = GetQinv(pVect1, pVect2);
1913 kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1914 SetFillBins2(ch1, ch2, bin1, bin2);
1915
1916 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1917 if(ch1 == ch2 && !fGeneratorOnly){
1918 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1919 if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1920 continue;
be9ef9f9 1921 }
b71263d0 1922 }
1923
1924 GetQosl(pVect1, pVect2, qout, qside, qlong);
1925 if( (en1+en2==0)) {
1926 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12);
1927 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2QW->Fill(kT12, qinv12, qinv12);
1928 // osl frame
1929 if((kT12 > 0.2) && (kT12 < 0.3)){
1930 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1931 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
be9ef9f9 1932 }
b71263d0 1933 if((kT12 > 0.6) && (kT12 < 0.7)){
1934 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1935 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
be9ef9f9 1936 }
b71263d0 1937 // unit mult bins
1938 if( (fEvt+en1)->fNtracks%100==0){
1939 Int_t kTindex=0;
1940 if(kT12>0.3) kTindex=1;
1941 Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
1942 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[0].fUnitMultBin->Fill(UnitMultBin, qinv12);
be9ef9f9 1943 }
1944
b71263d0 1945 }
1946 if( (en1+en2==1)) {
1947 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
1948 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2QW->Fill(kT12, qinv12, qinv12);
1949 // osl frame
1950 if((kT12 > 0.2) && (kT12 < 0.3)){
1951 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1952 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1953 }
1954 if((kT12 > 0.6) && (kT12 < 0.7)){
1955 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1956 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
be9ef9f9 1957 }
b71263d0 1958 // unit mult bins
1959 if( (fEvt+en1)->fNtracks%100==0){
1960 Int_t kTindex=0;
1961 if(kT12>0.3) kTindex=1;
1962 Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
1963 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[1].fUnitMultBin->Fill(UnitMultBin, qinv12);
be9ef9f9 1964 }
b71263d0 1965 }
1966 //////////////////////////////////////////
6bb6954b 1967 if(fTabulatePairs && en1==0 && en2<=1 && bin1==bin2){
b71263d0 1968 Float_t kY = 0;
1969 Int_t kTbin=-1, kYbin=-1;
be9ef9f9 1970
b71263d0 1971 for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(kT12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {kTbin = kIt; break;}}
1972 for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(kY < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {kYbin = kIt; break;}}
1973 if((kTbin<0) || (kYbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1974 if((kTbin>=fKbinsT) || (kYbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1975 if(fGenerateSignal && en2==0) {
1976 Int_t chGroup2[2]={ch1,ch2};
1977 Float_t WInput = MCWeight(chGroup2, fRMax, 0.7, qinv12, kT12);
1978 KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
1979 }else KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1980 }
1981
1982 //////////////////////////////////////////////////////////////////////////////
1983
1984 if(qinv12 <= fQcut) {
b71263d0 1985 if(en1==0 && en2==0) {fLowQPairSwitch_E0E0[i]->AddAt('1',j);}
1986 if(en1==0 && en2==1) {fLowQPairSwitch_E0E1[i]->AddAt('1',j);}
1987 if(en1==0 && en2==2) {fLowQPairSwitch_E0E2[i]->AddAt('1',j);}
1988 if(en1==0 && en2==3) {fLowQPairSwitch_E0E3[i]->AddAt('1',j);}
6bb6954b 1989 if(en1==1 && en2==1) {fLowQPairSwitch_E1E1[i]->AddAt('1',j);}
b71263d0 1990 if(en1==1 && en2==2) {fLowQPairSwitch_E1E2[i]->AddAt('1',j);}
1991 if(en1==1 && en2==3) {fLowQPairSwitch_E1E3[i]->AddAt('1',j);}
1992 if(en1==2 && en2==3) {fLowQPairSwitch_E2E3[i]->AddAt('1',j);}
1993 }
1994 if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) {
b71263d0 1995 if(en1==0 && en2==0) {fNormQPairSwitch_E0E0[i]->AddAt('1',j);}
1996 if(en1==0 && en2==1) {fNormQPairSwitch_E0E1[i]->AddAt('1',j);}
1997 if(en1==0 && en2==2) {fNormQPairSwitch_E0E2[i]->AddAt('1',j);}
1998 if(en1==0 && en2==3) {fNormQPairSwitch_E0E3[i]->AddAt('1',j);}
6bb6954b 1999 if(en1==1 && en2==1) {fNormQPairSwitch_E1E1[i]->AddAt('1',j);}
b71263d0 2000 if(en1==1 && en2==2) {fNormQPairSwitch_E1E2[i]->AddAt('1',j);}
2001 if(en1==1 && en2==3) {fNormQPairSwitch_E1E3[i]->AddAt('1',j);}
2002 if(en1==2 && en2==3) {fNormQPairSwitch_E2E3[i]->AddAt('1',j);}
be9ef9f9 2003 }
b71263d0 2004
be9ef9f9 2005 }
2006 }
2007 }
b71263d0 2008 }
be9ef9f9 2009
b71263d0 2010 //cout<<PairCount[0]<<" "<<PairCount[1]<<" "<<PairCount[2]<<" "<<PairCount[3]<<" "<<PairCount[4]<<" "<<PairCount[5]<<" "<<PairCount[6]<<endl;
2011 //cout<<NormPairCount[0]<<" "<<NormPairCount[1]<<" "<<NormPairCount[2]<<" "<<NormPairCount[3]<<" "<<NormPairCount[4]<<" "<<NormPairCount[5]<<" "<<NormPairCount[6]<<endl;
2012 ///////////////////////////////////////////////////
2013 // Do not use pairs from events with too many pairs
2014
2015 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2016
2017 ///////////////////////////////////////////////////
2018
2019
2020 if(fTabulatePairs) return;
be9ef9f9 2021
2022
b71263d0 2023
2024 ////////////////////////////////////////////////////
2025 ////////////////////////////////////////////////////
2026 // Normalization counting of 3- and 4-particle terms
2027 for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2028 for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
6bb6954b 2029 if(en2==0 && en3>2) continue;// not needed config
b71263d0 2030 if(en2==1 && en3==en2) continue;// not needed config
2031 for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2032 if(en3==0 && en4>1) continue;// not needed config
6bb6954b 2033 if(en3==1 && en4==3) continue;// not needed configs
2034 if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
b71263d0 2035
2036 for (Int_t i=0; i<myTracks; i++) {// 1st particle
2037 pVect1[1]=(fEvt)->fTracks[i].fP[0];
2038 pVect1[2]=(fEvt)->fTracks[i].fP[1];
2039 pVect1[3]=(fEvt)->fTracks[i].fP[2];
2040 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
be9ef9f9 2041
b71263d0 2042 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
b71263d0 2043 if(en2==0) {if(fNormQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2044 else {if(fNormQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2045
2046 pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2047 pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2048 pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2049 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2050
2051 for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
b71263d0 2052 if(en3==0) {
2053 if(fNormQPairSwitch_E0E0[i]->At(k)=='0') continue;
2054 if(fNormQPairSwitch_E0E0[j]->At(k)=='0') continue;
2055 }else if(en3==1){
2056 if(fNormQPairSwitch_E0E1[i]->At(k)=='0') continue;
2057 if(fNormQPairSwitch_E0E1[j]->At(k)=='0') continue;
2058 }else{
2059 if(fNormQPairSwitch_E0E2[i]->At(k)=='0') continue;
2060 if(fNormQPairSwitch_E1E2[j]->At(k)=='0') continue;
2061 }
be9ef9f9 2062
b71263d0 2063 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2064 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2065 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2066 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2067 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2068 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2069
2070 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2071 if(KT3<=fKT3transition) KT3index=0;
2072 else KT3index=1;
2073
2074 if(en2==0 && en3==0 && en4==0) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fNorm3->Fill(0);
2075 if(en2==1 && en3==2 && en4==3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fNorm3->Fill(0);
2076 if(en2==0 && en3==1 && en4==2) {
2077 if(fill2) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fNorm3->Fill(0);
2078 if(fill3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fNorm3->Fill(0);
2079 if(fill4) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fNorm3->Fill(0);
2080 }
2081
2082
2083 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
b71263d0 2084 if(en4==0){
2085 if(fNormQPairSwitch_E0E0[i]->At(l)=='0') continue;
2086 if(fNormQPairSwitch_E0E0[j]->At(l)=='0') continue;
2087 if(fNormQPairSwitch_E0E0[k]->At(l)=='0') continue;
2088 }else if(en4==1){
6bb6954b 2089 if(en3==0){
2090 if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2091 if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2092 if(fNormQPairSwitch_E0E1[k]->At(l)=='0') continue;
2093 }else{
2094 if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2095 if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2096 if(fNormQPairSwitch_E1E1[k]->At(l)=='0') continue;
2097 }
b71263d0 2098 }else if(en4==2){
2099 if(fNormQPairSwitch_E0E2[i]->At(l)=='0') continue;
2100 if(fNormQPairSwitch_E0E2[j]->At(l)=='0') continue;
2101 if(fNormQPairSwitch_E1E2[k]->At(l)=='0') continue;
be9ef9f9 2102 }else{
b71263d0 2103 if(fNormQPairSwitch_E0E3[i]->At(l)=='0') continue;
2104 if(fNormQPairSwitch_E1E3[j]->At(l)=='0') continue;
2105 if(fNormQPairSwitch_E2E3[k]->At(l)=='0') continue;
be9ef9f9 2106 }
be9ef9f9 2107
b71263d0 2108 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2109 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2110 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2111 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2112 Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
2113 if(KT4<=fKT4transition) KT4index=0;
2114 else KT4index=1;
be9ef9f9 2115
6bb6954b 2116 Bool_t FillTerms[13]={kFALSE};
b71263d0 2117 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, en2+en3+en4, FillTerms);
2118 //
6bb6954b 2119 for(int ft=0; ft<13; ft++) {
b71263d0 2120 if(FillTerms[ft]) Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fNorm4->Fill(0.);
be9ef9f9 2121 }
be9ef9f9 2122
b71263d0 2123
be9ef9f9 2124 }
2125 }
b71263d0 2126 }
2127 }
2128
be9ef9f9 2129 }
2130 }
b71263d0 2131 }
be9ef9f9 2132
6bb6954b 2133
be9ef9f9 2134
2135
2136 ///////////////////////////////////////////////////////////////////////
2137 ///////////////////////////////////////////////////////////////////////
2138 ///////////////////////////////////////////////////////////////////////
2139 //
2140 //
2141 // Start the Main Correlation Analysis
2142 //
2143 //
2144 ///////////////////////////////////////////////////////////////////////
6bb6954b 2145
be9ef9f9 2146
2147
2148 ////////////////////////////////////////////////////
2149 ////////////////////////////////////////////////////
2150 for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2151 for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
6bb6954b 2152 if(en2==0 && en3>2) continue;// not needed config
be9ef9f9 2153 if(en2==1 && en3==en2) continue;// not needed config
2154 for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2155 if(en3==0 && en4>1) continue;// not needed config
6bb6954b 2156 if(en3==1 && en4==3) continue;// not needed configs
2157 if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2158
be9ef9f9 2159 Int_t ENsum=en2+en3+en4;// 0 or 1 or 3 or 6
6bb6954b 2160
be9ef9f9 2161 /////////////////////////////////////////////////////////////
2162 for (Int_t i=0; i<myTracks; i++) {// 1st particle
2163 pVect1[0]=(fEvt)->fTracks[i].fEaccepted;
2164 pVect1[1]=(fEvt)->fTracks[i].fP[0];
2165 pVect1[2]=(fEvt)->fTracks[i].fP[1];
2166 pVect1[3]=(fEvt)->fTracks[i].fP[2];
2167 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2168
2169 /////////////////////////////////////////////////////////////
2170 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
b71263d0 2171 if(en2==0) {if(fLowQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2172 else {if(fLowQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2173
be9ef9f9 2174 pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2175 pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2176 pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2177 pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2178 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2179 qinv12 = GetQinv(pVect1, pVect2);
2180 kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2181 SetFillBins2(ch1, ch2, bin1, bin2);
2182 Int_t kTindex=0;
2183 if(kT12<=0.3) kTindex=0;
2184 else kTindex=1;
2185
2186 FSICorr12 = FSICorrelation(ch1,ch2, qinv12);
2187
2188 // two particle terms filled during tabulation of low-q pairs
2189
2190
2191 if(fMCcase){
2192 FilledMCpair12=kFALSE;
2193
6bb6954b 2194 if(ch1==ch2 && fMbin==0 && qinv12<0.2 && ENsum!=2 && ENsum!=3 && ENsum!=6){
be9ef9f9 2195 for(Int_t rstep=0; rstep<10; rstep++){
2196 Float_t coeff = (rstep)*0.2*(0.18/1.2);
2197 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2198 if(phi1 > 2*PI) phi1 -= 2*PI;
2199 if(phi1 < 0) phi1 += 2*PI;
2200 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2201 if(phi2 > 2*PI) phi2 -= 2*PI;
2202 if(phi2 < 0) phi2 += 2*PI;
2203 Float_t deltaphi = phi1 - phi2;
2204 if(deltaphi > PI) deltaphi -= PI;
2205 if(deltaphi < -PI) deltaphi += PI;
2206
2207 if(ENsum==0) ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2208 else ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2209 }
2210
2211 }// pair selection
2212
2213 // Check that label does not exceed stack size
2214 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
2215 if(ENsum==0 && abs((fEvt+en2)->fTracks[j].fLabel) == abs((fEvt)->fTracks[i].fLabel)) continue;
2216 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2217 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2218 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2219 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2220 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2221 qinv12MC = GetQinv(pVect1MC, pVect2MC);
2222
2223 if(qinv12<0.1 && ch1==ch2 && ENsum==0) {
2224 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
2225 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
2226 ((TH2D*)fOutputList->FindObject("fQ2Res"))->Fill(kT12, qinv12-qinv12MC);
2227 }
2228
2229 // secondary contamination
2230 if(ENsum==0){
2231 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
2232 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
2233 if(!mcParticle1 || !mcParticle2) continue;
2234 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
2235 if(ch1==ch2) {
2236 ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2237 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2238 ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2239 }
2240 }else{
2241 ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2242 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2243 ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2244 }
2245 }
2246 }
2247 }
2248
2249 if(ENsum==6){// all mixed events
2250 Int_t chGroup2[2]={ch1,ch2};
2251
2252 Float_t rForQW=5.0;
2253 if(fFSIindex<=1) rForQW=10;
2254 else if(fFSIindex==2) rForQW=9;
2255 else if(fFSIindex==3) rForQW=8;
2256 else if(fFSIindex==4) rForQW=7;
2257 else if(fFSIindex==5) rForQW=6;
2258 else if(fFSIindex==6) rForQW=5;
2259 else if(fFSIindex==7) rForQW=4;
2260 else if(fFSIindex==8) rForQW=3;
2261 else rForQW=2;
2262
2263
2264 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinv->Fill(qinv12MC, MCWeight(chGroup2, rForQW, 0.7, qinv12MC, 0.));// was 4,5
2265 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(chGroup2, rForQW, 0.7, qinv12MC, 0.));// was 4,5
2266 // pion purity
2267 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityDen->Fill(kT12, qinv12);
2268 Int_t SCNumber = 1;
2269 Int_t PdgCodeSum = abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode) + abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode);
2270 if(PdgCodeSum==22) SCNumber=1;// e-e
2271 else if(PdgCodeSum==24) SCNumber=2;// e-mu
2272 else if(PdgCodeSum==222) SCNumber=3;// e-pi
2273 else if(PdgCodeSum==332) SCNumber=4;// e-k
2274 else if(PdgCodeSum==2223) SCNumber=5;// e-p
2275 else if(PdgCodeSum==26) SCNumber=6;// mu-mu
2276 else if(PdgCodeSum==224) SCNumber=7;// mu-pi
2277 else if(PdgCodeSum==334) SCNumber=8;// mu-k
2278 else if(PdgCodeSum==2225) SCNumber=9;// mu-p
2279 else if(PdgCodeSum==422) SCNumber=10;// pi-pi
2280 else if(PdgCodeSum==532) SCNumber=11;// pi-k
2281 else if(PdgCodeSum==2423) SCNumber=12;// pi-p
2282 else if(PdgCodeSum==642) SCNumber=13;// k-k
2283 else if(PdgCodeSum==2533) SCNumber=14;// k-p
2284 else if(PdgCodeSum==4424) SCNumber=15;// p-p
2285 else {SCNumber=16;}
2286
2287 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityNum->Fill(SCNumber, kT12, qinv12);
2288
2289 ///////////////////////
2290 // muon contamination
2291 Pparent1[0]=pVect1MC[0]; Pparent1[1]=pVect1MC[1]; Pparent1[2]=pVect1MC[2]; Pparent1[3]=pVect1MC[3];
2292 Pparent2[0]=pVect2MC[0]; Pparent2[1]=pVect2MC[1]; Pparent2[2]=pVect2MC[2]; Pparent2[3]=pVect2MC[3];
2293 pionParent1=kFALSE; pionParent2=kFALSE;
2294 FilledMCpair12=kTRUE;
2295 //
2296 if(abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode)==13 || abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13){// muon check
2297 Int_t MotherLabel1 = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fMotherLabel;
2298 if(abs((fEvt)->fMCtracks[MotherLabel1].fPdgCode)==211) {
2299 pionParent1=kTRUE;
2300 Pparent1[1] = (fEvt)->fMCtracks[MotherLabel1].fPx; Pparent1[2] = (fEvt)->fMCtracks[MotherLabel1].fPy; Pparent1[3] = (fEvt)->fMCtracks[MotherLabel1].fPz;
2301 Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
2302 }
2303 //
2304 if(abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13) {
2305 Int_t MotherLabel2 = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fMotherLabel;
2306 if(abs((fEvt+en2)->fMCtracks[MotherLabel2].fPdgCode)==211) {
2307 pionParent2=kTRUE;
2308 Pparent2[1] = (fEvt+en2)->fMCtracks[MotherLabel2].fPx; Pparent2[2] = (fEvt+en2)->fMCtracks[MotherLabel2].fPy; Pparent2[3] = (fEvt+en2)->fMCtracks[MotherLabel2].fPz;
2309 Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
2310 }
2311 }
2312
2313 parentQinv12 = GetQinv(Pparent1, Pparent2);
2314
2315 if(pionParent1 || pionParent2){
2316 if(parentQinv12 > 0.001 && parentQinv12 < 0.3){
2317 Float_t muonPionK12 = FSICorrelation(ch1, ch2, qinv12MC);
2318 Float_t pionPionK12 = FSICorrelation(ch1, ch2, parentQinv12);
2319 for(Int_t term=1; term<=2; term++){
2320 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2321 Float_t Rvalue = 5+Riter;
2322 Float_t WInput = 1.0;
2323 if(term==1) {
2324 WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.);
5591748e 2325 }else{
be9ef9f9 2326 muonPionK12 = 1.0; pionPionK12=1.0;
2327 }
5591748e 2328
be9ef9f9 2329 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput);
2330 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput);
2331 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12);
2332 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fPionPionK2->Fill(Rvalue, parentQinv12, pionPionK12);
2333 }// Riter
2334 }// term loop
2335
2336 if(ch1==ch2) ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(0., kT12, qinv12MC-parentQinv12);
2337 else ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(1., kT12, qinv12MC-parentQinv12);
2338 }// parentQ check
2339 }// pion parent check
2340 }// muon check
2341
2342
2343 // momentum resolution
2344 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2345 Float_t Rvalue = 5+Riter;
2346 Float_t WInput = MCWeight(chGroup2, Rvalue, 0.7, qinv12MC, 0.);
2347 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fIdeal->Fill(Rvalue, qinv12MC, WInput);
2348 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fIdeal->Fill(Rvalue, qinv12MC);
2349 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fSmeared->Fill(Rvalue, qinv12, WInput);
2350 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fSmeared->Fill(Rvalue, qinv12);
2351 }
2352
2353 }// ENsum check
2354 }// MC array check
2355 }// MC case
2356
2357
6bb6954b 2358
be9ef9f9 2359 /////////////////////////////////////////////////////////////
2360 for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
b71263d0 2361 if(en3==0) {
2362 if(fLowQPairSwitch_E0E0[i]->At(k)=='0') continue;
2363 if(fLowQPairSwitch_E0E0[j]->At(k)=='0') continue;
2364 }else if(en3==1){
2365 if(fLowQPairSwitch_E0E1[i]->At(k)=='0') continue;
2366 if(fLowQPairSwitch_E0E1[j]->At(k)=='0') continue;
2367 }else{
2368 if(fLowQPairSwitch_E0E2[i]->At(k)=='0') continue;
2369 if(fLowQPairSwitch_E1E2[j]->At(k)=='0') continue;
be9ef9f9 2370 }
b71263d0 2371
be9ef9f9 2372 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
2373 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2374 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2375 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2376 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2377 qinv13 = GetQinv(pVect1, pVect3);
2378 qinv23 = GetQinv(pVect2, pVect3);
2379 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2380
2381 FilledMCtriplet123 = kFALSE;
2382
2383 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2384 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2385
2386 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2387 if(KT3<=fKT3transition) KT3index=0;
2388 else KT3index=1;
2389
2390 FSICorr13 = FSICorrelation(ch1,ch3, qinv13);
2391 FSICorr23 = FSICorrelation(ch2,ch3, qinv23);
2392
2393 if(ENsum==0) {
2394 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fTerms3->Fill(q3);
b71263d0 2395 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactor->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23));
be9ef9f9 2396 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv12);
2397 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv13);
2398 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv23);
2399 }
2400 if(ENsum==6) {
2401 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTerms3->Fill(q3);
2402 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv12);
2403 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv13);
2404 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv23);
2405 }
2406 if(ENsum==3){
2407 if(fill2) {
2408 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fTerms3->Fill(q3);
2409 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactor->Fill(q3, 1/(FSICorr12));
2410 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv12);
2411 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv13);
2412 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv23);
2413 }if(fill3) {
2414 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fTerms3->Fill(q3);
2415 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactor->Fill(q3, 1/(FSICorr12));
2416 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv12);
2417 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv13);
2418 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv23);
2419 }if(fill4) {
2420 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fTerms3->Fill(q3);
2421 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactor->Fill(q3, 1/(FSICorr12));
2422 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv12);
2423 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv13);
2424 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv23);
2425 }
2426 }
2427
2428 // r3 denominator
2429 if(ENsum==6 && ch1==ch2 && ch1==ch3){
2430 GoodTripletWeights = kFALSE;
2431 //
2432 GetWeight(pVect1, pVect2, weight12, weight12Err);
2433 GetWeight(pVect1, pVect3, weight13, weight13Err);
2434 GetWeight(pVect2, pVect3, weight23, weight23Err);
2435
2436 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {// weight should never be larger than 1
2437 if(fMbin==0 && bin1==0) {
2438 ((TH1D*)fOutputList->FindObject("fTPNRejects3pion1"))->Fill(q3, sqrt(fabs(weight12*weight13*weight23)));
2439 }
2440 }else{
2441
2442 Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
5591748e 2443 Float_t MuonCorr12=1.0, MuonCorr13=1.0, MuonCorr23=1.0;
be9ef9f9 2444 if(!fGenerateSignal && !fMCcase) {
2445 Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
2446 Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
2447 Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);
2448 if(momBin12 >= 20) momBin12 = 19;
2449 if(momBin13 >= 20) momBin13 = 19;
2450 if(momBin23 >= 20) momBin23 = 19;
2451 MomResCorr12 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin12);
2452 MomResCorr13 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin13);
2453 MomResCorr23 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin23);
5591748e 2454 MuonCorr12 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin12);
2455 MuonCorr13 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin13);
2456 MuonCorr23 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin23);
be9ef9f9 2457 }
5591748e 2458 // no MRC, no Muon Correction
2459 weight12CC[0] = ((weight12+1) - ffcSq*FSICorr12 - (1-ffcSq));
2460 weight12CC[0] /= FSICorr12*ffcSq;
2461 weight13CC[0] = ((weight13+1) - ffcSq*FSICorr13 - (1-ffcSq));
2462 weight13CC[0] /= FSICorr13*ffcSq;
2463 weight23CC[0] = ((weight23+1) - ffcSq*FSICorr23 - (1-ffcSq));
2464 weight23CC[0] /= FSICorr23*ffcSq;
2465 if(weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2466 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, sqrt(weight12CC[0]*weight13CC[0]*weight23CC[0]));
2467 }
2468 // no Muon Correction
2469 weight12CC[1] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2470 weight12CC[1] /= FSICorr12*ffcSq;
2471 weight13CC[1] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2472 weight13CC[1] /= FSICorr13*ffcSq;
2473 weight23CC[1] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2474 weight23CC[1] /= FSICorr23*ffcSq;
2475 if(weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2476 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(2, q3, sqrt(weight12CC[1]*weight13CC[1]*weight23CC[1]));
2477 }
2478 // both Corrections
2479 weight12CC[2] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2480 weight12CC[2] /= FSICorr12*ffcSq;
2481 weight12CC[2] *= MuonCorr12;
2482 weight13CC[2] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2483 weight13CC[2] /= FSICorr13*ffcSq;
2484 weight13CC[2] *= MuonCorr13;
2485 weight23CC[2] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2486 weight23CC[2] /= FSICorr23*ffcSq;
2487 weight23CC[2] *= MuonCorr23;
2488 if(weight12CC[2] < 0 || weight13CC[2] < 0 || weight23CC[2] < 0) {// C2^QS can never be less than unity
be9ef9f9 2489 if(fMbin==0 && bin1==0) {
5591748e 2490 ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC[2]*weight13CC[2]*weight23CC[2])));
be9ef9f9 2491 }
2492 }else{
2493 GoodTripletWeights = kTRUE;
2494 /////////////////////////////////////////////////////
5591748e 2495 weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]);
be9ef9f9 2496 /////////////////////////////////////////////////////
5591748e 2497 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(3, q3, weightTotal);
2498 }// 2nd r3 den check
2499
2500
be9ef9f9 2501 }// 1st r3 den check
5591748e 2502
be9ef9f9 2503 }// r3 den
2504
2505
2506 if(ch1==ch2 && ch1==ch3 && ENsum==0){
2507 ((TH3D*)fOutputList->FindObject("fKT3DistTerm1"))->Fill(fMbin+1, KT3, q3);
2508 if(q3<0.06){
2509 Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2510 Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2511 Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2512 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt1);
2513 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt2);
2514 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt3);
2515 }
2516 }
2517 if(ch1==ch2 && ch1==ch3 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT3DistTerm5"))->Fill(fMbin+1, KT3, q3);
2518
2519
2520
2521
2522 if(fMCcase && ENsum==6 && FilledMCpair12){// for momentum resolution and muon correction
2523 if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
2524
2525 pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2526 pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2527 pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2528 pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2529 qinv13MC = GetQinv(pVect1MC, pVect3MC);
2530 qinv23MC = GetQinv(pVect2MC, pVect3MC);
2531
2532 q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
2533 if(q3<0.1 && ch1==ch2 && ch1==ch3) ((TH2D*)fOutputList->FindObject("fQ3Res"))->Fill(KT3, q3-q3MC);
2534
2535 Int_t chGroup3[3]={ch1,ch2,ch3};
2536 Float_t QinvMCGroup3[3]={qinv12MC,qinv13MC,qinv23MC};
2537 //Float_t kTGroup3[3]={float(sqrt(pow(pVect1MC[1]+pVect2MC[1],2) + pow(pVect1MC[2]+pVect2MC[2],2))/2.),
2538 //float(sqrt(pow(pVect1MC[1]+pVect3MC[1],2) + pow(pVect1MC[2]+pVect3MC[2],2))/2.),
2539 //float(sqrt(pow(pVect2MC[1]+pVect3MC[1],2) + pow(pVect2MC[2]+pVect3MC[2],2))/2.)};
2540 Float_t kTGroup3[3]={0};
2541
2542 Pparent3[0]=pVect3MC[0]; Pparent3[1]=pVect3MC[1]; Pparent3[2]=pVect3MC[2]; Pparent3[3]=pVect3MC[3];
2543 pionParent3=kFALSE;
2544
2545 if(abs((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPdgCode)==13){// muon check
2546 Int_t MotherLabel3 = (fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fMotherLabel;
2547 if(abs((fEvt+en3)->fMCtracks[MotherLabel3].fPdgCode)==211) {
2548 pionParent3=kTRUE;
2549 Pparent3[1] = (fEvt+en3)->fMCtracks[MotherLabel3].fPx; Pparent3[2] = (fEvt+en3)->fMCtracks[MotherLabel3].fPy; Pparent3[3] = (fEvt+en3)->fMCtracks[MotherLabel3].fPz;
2550 Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
2551 }
2552 }
2553
2554 parentQinv13 = GetQinv(Pparent1, Pparent3);
2555 parentQinv23 = GetQinv(Pparent2, Pparent3);
2556 parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
2557
2558 if(parentQinv12 > 0.001 && parentQinv13 > 0.001 && parentQinv23 > 0.001 && parentQ3 < 0.5){
2559 FilledMCtriplet123=kTRUE;
2560 if(pionParent1 || pionParent2 || pionParent3) {// want at least one pion-->muon
2561
2562 Float_t parentQinvGroup3[3]={parentQinv12, parentQinv13, parentQinv23};
2563 //Float_t parentkTGroup3[3]={float(sqrt(pow(Pparent1[1]+Pparent2[1],2) + pow(Pparent1[2]+Pparent2[2],2))/2.),
2564 //float(sqrt(pow(Pparent1[1]+Pparent3[1],2) + pow(Pparent1[2]+Pparent3[2],2))/2.),
2565 //float(sqrt(pow(Pparent2[1]+Pparent3[1],2) + pow(Pparent2[2]+Pparent3[2],2))/2.)};
2566 Float_t parentkTGroup3[3]={0};
2567
2568 ((TH2D*)fOutputList->FindObject("fAvgQ12VersusQ3"))->Fill(parentQ3, parentQinv12);
2569 ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13);
2570 ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23);
2571
5591748e 2572 for(Int_t term=1; term<=4; term++){
2573 if(term==1) {}
2574 else if(term==2) {if(!pionParent1 && !pionParent2) continue;}
2575 else if(term==3) {if(!pionParent1 && !pionParent3) continue;}
2576 else {if(!pionParent2 && !pionParent3) continue;}
be9ef9f9 2577 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2578 Float_t Rvalue = 5+Riter;
2579 Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
2580 Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
2581 Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
5591748e 2582 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(1, Rvalue, q3MC, WInput);
2583 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ3, WInput);
2584 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(1, Rvalue, q3MC, WInputFSI);
2585 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(1, Rvalue, parentQ3, WInputParentFSI);
2586 //
2587 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(2, Rvalue, q3MC);
2588 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ3);
2589 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(2, Rvalue, q3MC);
2590 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(2, Rvalue, parentQ3);
be9ef9f9 2591 }// Riter
2592 }// term loop
2593
2594 }// pion parent check
2595 }// parentQ check (muon correction)
2596
2597 // 3-pion momentum resolution
2598 for(Int_t term=1; term<=5; term++){
2599 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2600 Float_t Rvalue = 5+Riter;
2601 Float_t WInput = MCWeight3(term, Rvalue, 0.7, chGroup3, QinvMCGroup3, kTGroup3);
2602 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fIdeal->Fill(Rvalue, q3MC, WInput);
2603 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fSmeared->Fill(Rvalue, q3, WInput);
2604 }
2605 }
2606
2607 }// 3rd particle label check
2608 }// MCcase and ENsum==6
2609
5591748e 2610
be9ef9f9 2611
2612
5591748e 2613 /////////////////////////////////////////////////////////////
be9ef9f9 2614 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
b71263d0 2615 if(en4==0){
2616 if(fLowQPairSwitch_E0E0[i]->At(l)=='0') continue;
2617 if(fLowQPairSwitch_E0E0[j]->At(l)=='0') continue;
2618 if(fLowQPairSwitch_E0E0[k]->At(l)=='0') continue;
2619 }else if(en4==1){
6bb6954b 2620 if(en3==0){
2621 if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2622 if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2623 if(fLowQPairSwitch_E0E1[k]->At(l)=='0') continue;
2624 }else{
2625 if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2626 if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2627 if(fLowQPairSwitch_E1E1[k]->At(l)=='0') continue;
2628 }
b71263d0 2629 }else if(en4==2){
2630 if(fLowQPairSwitch_E0E2[i]->At(l)=='0') continue;
2631 if(fLowQPairSwitch_E0E2[j]->At(l)=='0') continue;
2632 if(fLowQPairSwitch_E1E2[k]->At(l)=='0') continue;
2633 }else{
2634 if(fLowQPairSwitch_E0E3[i]->At(l)=='0') continue;
2635 if(fLowQPairSwitch_E1E3[j]->At(l)=='0') continue;
2636 if(fLowQPairSwitch_E2E3[k]->At(l)=='0') continue;
be9ef9f9 2637 }
6bb6954b 2638
be9ef9f9 2639 pVect4[0]=(fEvt+en4)->fTracks[l].fEaccepted;
2640 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2641 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2642 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2643 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2644 qinv14 = GetQinv(pVect1, pVect4);
2645 qinv24 = GetQinv(pVect2, pVect4);
2646 qinv34 = GetQinv(pVect3, pVect4);
2647 q4 = sqrt(pow(q3,2) + pow(qinv14,2) + pow(qinv24,2) + pow(qinv34,2));
2648
2649 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2650 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(1, qinv12); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(2, qinv13);
2651 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(3, qinv14); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(4, qinv23);
2652 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(5, qinv24); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(6, qinv34);
2653 }
2654
2655 Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
2656 if(KT4<=fKT4transition) KT4index=0;
2657 else KT4index=1;
2658
2659 FSICorr14 = FSICorrelation(ch1,ch4, qinv14);
2660 FSICorr24 = FSICorrelation(ch2,ch4, qinv24);
2661 FSICorr34 = FSICorrelation(ch3,ch4, qinv34);
2662
6bb6954b 2663 Bool_t FillTerms[13]={kFALSE};
be9ef9f9 2664 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, ENsum, FillTerms);
2665 //
6bb6954b 2666 for(int ft=0; ft<13; ft++) {
be9ef9f9 2667 Float_t FSIfactor = 1.0;
2668 if(ft==0) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr14 * FSICorr23 * FSICorr24 * FSICorr34);
2669 else if(ft<=4) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr23);
2670 else if(ft<=10) FSIfactor = 1/(FSICorr12);
6bb6954b 2671 else if(ft==11) FSIfactor = 1/(FSICorr12 * FSICorr34);
be9ef9f9 2672 else FSIfactor = 1.0;
2673 if(FillTerms[ft]) {
2674 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fTerms4->Fill(q4);
2675 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactor->Fill(q4, FSIfactor);
2676 }
2677 }
2678
2679 /////////////////////////////////////////////////////////////
2680 // r4{2}
2681 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6 && GoodTripletWeights){
2682 GetWeight(pVect1, pVect4, weight14, weight14Err);
2683 GetWeight(pVect2, pVect4, weight24, weight24Err);
2684 GetWeight(pVect3, pVect4, weight34, weight34Err);
2685
2686 Float_t MomResCorr14=1.0, MomResCorr24=1.0, MomResCorr34=1.0;
5591748e 2687 Float_t MuonCorr14=1.0, MuonCorr24=1.0, MuonCorr34=1.0;
be9ef9f9 2688 if(!fGenerateSignal && !fMCcase) {
2689 Int_t momBin14 = fMomResC2->GetYaxis()->FindBin(qinv14);
2690 Int_t momBin24 = fMomResC2->GetYaxis()->FindBin(qinv24);
2691 Int_t momBin34 = fMomResC2->GetYaxis()->FindBin(qinv34);
2692 if(momBin14 >= 20) momBin14 = 19;
2693 if(momBin24 >= 20) momBin24 = 19;
2694 if(momBin34 >= 20) momBin34 = 19;
5591748e 2695 MomResCorr14 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin14);
2696 MomResCorr24 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin24);
2697 MomResCorr34 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin34);
2698 MuonCorr14 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin14);
2699 MuonCorr24 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin24);
2700 MuonCorr34 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin34);
be9ef9f9 2701 }
5591748e 2702
2703 // no MRC, no Muon Correction
2704 weight14CC[0] = ((weight14+1) - ffcSq*FSICorr14 - (1-ffcSq));
2705 weight14CC[0] /= FSICorr14*ffcSq;
2706 weight24CC[0] = ((weight24+1) - ffcSq*FSICorr24 - (1-ffcSq));
2707 weight24CC[0] /= FSICorr24*ffcSq;
2708 weight34CC[0] = ((weight34+1) - ffcSq*FSICorr34 - (1-ffcSq));
2709 weight34CC[0] /= FSICorr34*ffcSq;
2710 if(weight14CC[0] > 0 && weight24CC[0] > 0 && weight34CC[0] > 0 && weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2711 weightTotal = sqrt(weight12CC[0]*weight13CC[0]*weight24CC[0]*weight34CC[0]);
2712 weightTotal += sqrt(weight12CC[0]*weight14CC[0]*weight23CC[0]*weight34CC[0]);
2713 weightTotal += sqrt(weight13CC[0]*weight14CC[0]*weight23CC[0]*weight24CC[0]);
2714 weightTotal /= 3.;
2715 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal);
2716 }
2717 // no Muon Correction
2718 weight14CC[1] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2719 weight14CC[1] /= FSICorr14*ffcSq;
2720 weight24CC[1] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2721 weight24CC[1] /= FSICorr24*ffcSq;
2722 weight34CC[1] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2723 weight34CC[1] /= FSICorr34*ffcSq;
2724 if(weight14CC[1] > 0 && weight24CC[1] > 0 && weight34CC[1] > 0 && weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2725 weightTotal = sqrt(weight12CC[1]*weight13CC[1]*weight24CC[1]*weight34CC[1]);
2726 weightTotal += sqrt(weight12CC[1]*weight14CC[1]*weight23CC[1]*weight34CC[1]);
2727 weightTotal += sqrt(weight13CC[1]*weight14CC[1]*weight23CC[1]*weight24CC[1]);
2728 weightTotal /= 3.;
2729 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(2, q4, weightTotal);
2730 }
2731 // both corrections
2732 weight14CC[2] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2733 weight14CC[2] /= FSICorr14*ffcSq;
2734 weight14CC[2] *= MuonCorr14;
2735 weight24CC[2] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2736 weight24CC[2] /= FSICorr24*ffcSq;
2737 weight24CC[2] *= MuonCorr24;
2738 weight34CC[2] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2739 weight34CC[2] /= FSICorr34*ffcSq;
2740 weight34CC[2] *= MuonCorr34;
2741
2742 if(weight14CC[2] < 0 || weight24CC[2] < 0 || weight34CC[2] < 0) {// C2^QS can never be less than unity
2743 if(fMbin==0 && bin1==0 && KT4index==0) {
2744 ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC[2]*weight23CC[2]*weight34CC[2]*weight14CC[2])));
be9ef9f9 2745 }
2746 }else{
2747 /////////////////////////////////////////////////////
5591748e 2748 weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2]);
2749 weightTotal += sqrt(weight12CC[2]*weight14CC[2]*weight23CC[2]*weight34CC[2]);
2750 weightTotal += sqrt(weight13CC[2]*weight14CC[2]*weight23CC[2]*weight24CC[2]);
be9ef9f9 2751 weightTotal /= 3.;
2752 /////////////////////////////////////////////////////
5591748e 2753 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(3, q4, weightTotal);
be9ef9f9 2754 }
2755 }
2756 /////////////////////////////////////////////////////////////
6bb6954b 2757
be9ef9f9 2758 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==0){
2759 ((TH3D*)fOutputList->FindObject("fKT4DistTerm1"))->Fill(fMbin+1, KT4, q4);
2760 if(q4<0.085){
2761 Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2762 Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2763 Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2764 Float_t pt4=sqrt(pow(pVect4[1],2)+pow(pVect4[2],2));
2765 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt1);
2766 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt2);
2767 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt3);
2768 ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt4);
2769 }
2770 }
6bb6954b 2771 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT4DistTerm13"))->Fill(fMbin+1, KT4, q4);
2772
be9ef9f9 2773
2774 // momenumtum resolution and muon corrections
2775 if(fMCcase && ENsum==6 && FilledMCtriplet123){// for momentum resolution and muon correction
2776 if((fEvt+en4)->fTracks[l].fLabel < (fEvt+en4)->fMCarraySize){
2777
2778 pVect4MC[0]=sqrt(pow((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2779 pVect4MC[1]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPx;
2780 pVect4MC[2]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPy;
2781 pVect4MC[3]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPz;
2782 qinv14MC = GetQinv(pVect1MC, pVect4MC);
2783 qinv24MC = GetQinv(pVect2MC, pVect4MC);
2784 qinv34MC = GetQinv(pVect3MC, pVect4MC);
2785
2786 q4MC = sqrt(pow(q3MC,2) + pow(qinv14MC,2) + pow(qinv24MC,2) + pow(qinv34MC,2));
2787 if(q4<0.1 && ch1==ch2 && ch1==ch3 && ch1==ch4) ((TH2D*)fOutputList->FindObject("fQ4Res"))->Fill(KT4, q4-q4MC);
2788 if(ch1==ch2 && ch1==ch3 && ch1==ch4) {
5591748e 2789 ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(1, qinv12MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(2, qinv13MC);
2790 ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(3, qinv14MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(4, qinv23MC);
be9ef9f9 2791 ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(5, qinv24MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(6, qinv34MC);
2792 }
2793 Int_t chGroup4[4]={ch1,ch2,ch3,ch4};
2794 Float_t QinvMCGroup4[6]={qinv12MC, qinv13MC, qinv14MC, qinv23MC, qinv24MC, qinv34MC};
2795 Float_t kTGroup4[6]={0};
2796
2797 Pparent4[0]=pVect4MC[0]; Pparent4[1]=pVect4MC[1]; Pparent4[2]=pVect4MC[2]; Pparent4[3]=pVect4MC[3];
2798 pionParent4=kFALSE;
2799 if(abs((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPdgCode)==13){// muon check
2800 Int_t MotherLabel4 = (fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fMotherLabel;
2801 if(abs((fEvt+en4)->fMCtracks[MotherLabel4].fPdgCode)==211) {
2802 pionParent4=kTRUE;
2803 Pparent4[1] = (fEvt+en4)->fMCtracks[MotherLabel4].fPx; Pparent4[2] = (fEvt+en4)->fMCtracks[MotherLabel4].fPy; Pparent4[3] = (fEvt+en4)->fMCtracks[MotherLabel4].fPz;
2804 Pparent4[0] = sqrt(pow(Pparent4[1],2)+pow(Pparent4[2],2)+pow(Pparent4[3],2)+pow(fTrueMassPi,2));
2805 }
2806 }
2807
2808 parentQinv14 = GetQinv(Pparent1, Pparent4);
2809 parentQinv24 = GetQinv(Pparent2, Pparent4);
2810 parentQinv34 = GetQinv(Pparent3, Pparent4);
2811 Float_t parentQ4 = sqrt(pow(parentQ3,2) + pow(parentQinv14,2) + pow(parentQinv24,2) + pow(parentQinv34,2));
2812
2813 if(parentQinv14 > 0.001 && parentQinv24 > 0.001 && parentQinv34 > 0.001 && parentQ4 < 0.5){
2814 if(pionParent1 || pionParent2 || pionParent3 || pionParent4) {// want at least one pion-->muon
2815
5591748e 2816 if(pionParent1) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(1);
2817 if(pionParent2) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(2);
2818 if(pionParent3) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(3);
2819 if(pionParent4) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(4);
be9ef9f9 2820 Float_t parentQinvGroup4[6]={parentQinv12, parentQinv13, parentQinv14, parentQinv23, parentQinv24, parentQinv34};
2821 Float_t parentkTGroup4[6]={0};
2822
5591748e 2823 for(Int_t term=1; term<=12; term++){
2824 if(term==1) {}
2825 else if(term==2) {if(!pionParent1 && !pionParent2 && !pionParent3) continue;}
2826 else if(term==3) {if(!pionParent1 && !pionParent2 && !pionParent4) continue;}
2827 else if(term==4) {if(!pionParent1 && !pionParent3 && !pionParent4) continue;}
2828 else if(term==5) {if(!pionParent2 && !pionParent3 && !pionParent4) continue;}
2829 else if(term==6) {if(!pionParent1 && !pionParent2) continue;}
2830 else if(term==7) {if(!pionParent1 && !pionParent3) continue;}
2831 else if(term==8) {if(!pionParent1 && !pionParent4) continue;}
2832 else if(term==9) {if(!pionParent2 && !pionParent3) continue;}
2833 else if(term==10) {if(!pionParent2 && !pionParent4) continue;}
2834 else if(term==11) {if(!pionParent3 && !pionParent4) continue;}
2835 else {}
be9ef9f9 2836 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2837 Float_t Rvalue = 5+Riter;
2838 Float_t WInput = MCWeight4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4, parentkTGroup4);
2839 Float_t WInputParentFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4);
2840 Float_t WInputFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, QinvMCGroup4);
5591748e 2841 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(1, Rvalue, q4MC, WInput);
2842 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ4, WInput);
2843 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(1, Rvalue, q4MC, WInputFSI);
2844 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(1, Rvalue, parentQ4, WInputParentFSI);
2845 //
2846 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(2, Rvalue, q4MC);
2847 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ4);
2848 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(2, Rvalue, q4MC);
2849 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(2, Rvalue, parentQ4);
be9ef9f9 2850 }// Riter
2851 }// term loop
2852
2853 }// pion parent check
2854 }// parentQ check (muon correction)
2855
2856 // 4-pion momentum resolution
6bb6954b 2857 for(Int_t term=1; term<=13; term++){
be9ef9f9 2858 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2859 Float_t Rvalue = 5+Riter;
2860 Float_t WInput = MCWeight4(term, Rvalue, 0.7, chGroup4, QinvMCGroup4, kTGroup4);
2861 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fIdeal->Fill(Rvalue, q4MC, WInput);
2862 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fSmeared->Fill(Rvalue, q4, WInput);
2863 }
2864 }
2865
2866 }// label check particle 4
2867 }// MCcase
2868
2869 }// 4th particle
2870 }// 3rd particle
2871 }// 2nd particle
2872 }// 1st particle
2873
2874 }// en4
2875 }// en3
2876 }// en2
2877
2878
2879
2880
2881
2882
2883
2884 // Post output data.
2885 PostData(1, fOutputList);
2886
2887}
2888//________________________________________________________________________
2889void AliFourPion::Terminate(Option_t *)
2890{
2891 // Called once at the end of the query
2892
2893 cout<<"Done"<<endl;
2894
2895}
2896//________________________________________________________________________
2897Bool_t AliFourPion::AcceptPair(AliFourPionTrackStruct first, AliFourPionTrackStruct second)
2898{
2899
2900 if(fabs(first.fEta-second.fEta) > fMinSepPairEta) return kTRUE;
2901
2902 // propagate through B field to r=1m
2903 Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
2904 if(phi1 > 2*PI) phi1 -= 2*PI;
2905 if(phi1 < 0) phi1 += 2*PI;
2906 Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m
2907 if(phi2 > 2*PI) phi2 -= 2*PI;
2908 if(phi2 < 0) phi2 += 2*PI;
2909
2910 Float_t deltaphi = phi1 - phi2;
2911 if(deltaphi > PI) deltaphi -= 2*PI;
2912 if(deltaphi < -PI) deltaphi += 2*PI;
2913 deltaphi = fabs(deltaphi);
2914
2915 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
2916
2917
2918 // propagate through B field to r=1.6m
2919 phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
2920 if(phi1 > 2*PI) phi1 -= 2*PI;
2921 if(phi1 < 0) phi1 += 2*PI;
2922 phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m
2923 if(phi2 > 2*PI) phi2 -= 2*PI;
2924 if(phi2 < 0) phi2 += 2*PI;
2925
2926 deltaphi = phi1 - phi2;
2927 if(deltaphi > PI) deltaphi -= 2*PI;
2928 if(deltaphi < -PI) deltaphi += 2*PI;
2929 deltaphi = fabs(deltaphi);
2930
2931 if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
2932
2933
2934
2935 //
2936
2937 Int_t ncl1 = first.fClusterMap.GetNbits();
2938 Int_t ncl2 = second.fClusterMap.GetNbits();
2939 Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
2940 Double_t shfrac = 0; Double_t qfactor = 0;
2941 for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
2942 if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
2943 if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
2944 sumQ++;
2945 sumCls+=2;
2946 sumSha+=2;}
2947 else {sumQ--; sumCls+=2;}
2948 }
2949 else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
2950 sumQ++;
2951 sumCls++;}
2952 }
2953 if (sumCls>0) {
2954 qfactor = sumQ*1.0/sumCls;
2955 shfrac = sumSha*1.0/sumCls;
2956 }
2957
2958 if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
2959
2960
2961 return kTRUE;
2962
2963
2964}
2965//________________________________________________________________________
2966Float_t AliFourPion::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
2967{
2968 Float_t arg = G_Coeff/qinv;
2969
2970 if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
2971 else {return (exp(-arg)-1)/(-arg);}
2972
2973}
2974//________________________________________________________________________
2975void AliFourPion::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
2976{
2977 Int_t j, k;
2978 Int_t a = i2 - i1;
2979 for (Int_t i = i1; i < i2+1; i++) {
2980 j = (Int_t) (gRandom->Rndm() * a);
2981 k = iarr[j];
2982 iarr[j] = iarr[i];
2983 iarr[i] = k;
2984 }
2985}
2986
2987
2988//________________________________________________________________________
2989Float_t AliFourPion::GetQinv(Float_t track1[], Float_t track2[]){
2990
2991 Float_t qinv = sqrt( fabs(pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2)) );
2992 return qinv;
2993
2994}
2995//________________________________________________________________________
2996void AliFourPion::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
2997
2998 Float_t p0 = track1[0] + track2[0];
2999 Float_t px = track1[1] + track2[1];
3000 Float_t py = track1[2] + track2[2];
3001 Float_t pz = track1[3] + track2[3];
3002
3003 Float_t mt = sqrt(p0*p0 - pz*pz);
3004 Float_t pt = sqrt(px*px + py*py);
3005
3006 Float_t v0 = track1[0] - track2[0];
3007 Float_t vx = track1[1] - track2[1];
3008 Float_t vy = track1[2] - track2[2];
3009 Float_t vz = track1[3] - track2[3];
3010
3011 qout = (px*vx + py*vy)/pt;
3012 qside = (px*vy - py*vx)/pt;
3013 qlong = (p0*vz - pz*v0)/mt;
3014}
3015//________________________________________________________________________
3016void AliFourPion::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliFourPion::fKbinsT][AliFourPion::fCentBins]){
3017
3018 if(legoCase){
3019 cout<<"LEGO call to SetWeightArrays"<<endl;
3020
3021 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3022 for(Int_t mb=0; mb<fCentBins; mb++){
3023 fNormWeight[tKbin][mb] = (TH3F*)histos[tKbin][mb]->Clone();
3024 fNormWeight[tKbin][mb]->SetDirectory(0);
3025 }
3026 }
3027
3028 }else{
3029
3030 TFile *wFile = new TFile("WeightFile.root","READ");
3031 if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3032 else cout<<"Good Weight File Found!"<<endl;
3033
3034 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3035 for(Int_t mb=0; mb<fCentBins; mb++){
3036
3037 TString *name = new TString("Weight_Kt_");
3038 *name += tKbin;
3039 name->Append("_Ky_0");
3040 name->Append("_M_");
3041 *name += mb;
3042 name->Append("_ED_0");
3043
3044
3045 fNormWeight[tKbin][mb] = (TH3F*)wFile->Get(name->Data());
3046 fNormWeight[tKbin][mb]->SetDirectory(0);
3047
3048
3049 }//mb
3050 }//kt
3051
3052 wFile->Close();
3053 }
3054
3055 cout<<"Done reading weight file"<<endl;
3056
3057}
3058//________________________________________________________________________
3059void AliFourPion::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3060
3061 Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3062 //
3063 Float_t qOut=0,qSide=0,qLong=0;
3064 GetQosl(track1, track2, qOut, qSide, qLong);
3065 qOut = fabs(qOut);
3066 qSide = fabs(qSide);
3067 qLong = fabs(qLong);
3068 Float_t wd=0, xd=0, yd=0, zd=0;
3069 //Float_t qinvtemp=GetQinv(0,track1, track2);
3070 //
3071
3072 if(kt < fKmeanT[0]) {fKtIndexL=0; fKtIndexH=1;}
3073 else if(kt >= fKmeanT[fKbinsT-1]) {fKtIndexL=fKbinsT-2; fKtIndexH=fKbinsT-1;}
3074 else {
3075 for(Int_t i=0; i<fKbinsT-1; i++){
3076 if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtIndexL=i; fKtIndexH=i+1; break;}
3077 }
3078 }
3079 wd = (kt-fKmeanT[fKtIndexL])/(fKmeanT[fKtIndexH]-fKmeanT[fKtIndexL]);
3080 //
3081 /////////
3082 if(qOut < fQmean[0]) {fQoIndexL=0; fQoIndexH=0; xd=0;}
3083 else if(qOut >= fQmean[kQbinsWeights-1]) {fQoIndexL=kQbinsWeights-1; fQoIndexH=kQbinsWeights-1; xd=1;}
3084 else {
3085 for(Int_t i=0; i<kQbinsWeights-1; i++){
3086 if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQoIndexL=i; fQoIndexH=i+1; break;}
3087 }
3088 xd = (qOut-fQmean[fQoIndexL])/(fQmean[fQoIndexH]-fQmean[fQoIndexL]);
3089 }
3090 //
3091 if(qSide < fQmean[0]) {fQsIndexL=0; fQsIndexH=0; yd=0;}
3092 else if(qSide >= fQmean[kQbinsWeights-1]) {fQsIndexL=kQbinsWeights-1; fQsIndexH=kQbinsWeights-1; yd=1;}
3093 else {
3094 for(Int_t i=0; i<kQbinsWeights-1; i++){
3095 if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsIndexL=i; fQsIndexH=i+1; break;}
3096 }
3097 yd = (qSide-fQmean[fQsIndexL])/(fQmean[fQsIndexH]-fQmean[fQsIndexL]);
3098 }
3099 //
3100 if(qLong < fQmean[0]) {fQlIndexL=0; fQlIndexH=0; zd=0;}
3101 else if(qLong >= fQmean[kQbinsWeights-1]) {fQlIndexL=kQbinsWeights-1; fQlIndexH=kQbinsWeights-1; zd=1;}
3102 else {
3103 for(Int_t i=0; i<kQbinsWeights-1; i++){
3104 if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlIndexL=i; fQlIndexH=i+1; break;}
3105 }
3106 zd = (qLong-fQmean[fQlIndexL])/(fQmean[fQlIndexH]-fQmean[fQlIndexL]);
3107 }
3108 //
3109
3110
3111 //Float_t min = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3112 Float_t minErr = fNormWeight[fKtIndexL][fMbin]->GetBinError(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3113 /*
3114 Float_t deltaW=0;
3115 // kt
3116 deltaW += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - min)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3117 // Qo
3118 deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - min)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3119 // Qs
3120 deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - min)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3121 // Ql
3122 deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - min)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3123 //
3124 wgt = min + deltaW;
3125 */
3126
3127
3128 //
3129 // w interpolation (kt)
3130 Float_t c000 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*wd;
3131 Float_t c100 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*wd;
3132 Float_t c010 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*wd;
3133 Float_t c001 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*wd;
3134 Float_t c110 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*wd;
3135 Float_t c101 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*wd;
3136 Float_t c011 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*wd;
3137 Float_t c111 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*wd;
3138 // x interpolation (qOut)
3139 Float_t c00 = c000*(1-xd) + c100*xd;
3140 Float_t c10 = c010*(1-xd) + c110*xd;
3141 Float_t c01 = c001*(1-xd) + c101*xd;
3142 Float_t c11 = c011*(1-xd) + c111*xd;
3143 // y interpolation (qSide)
3144 Float_t c0 = c00*(1-yd) + c10*yd;
3145 Float_t c1 = c01*(1-yd) + c11*yd;
3146 // z interpolation (qLong)
3147 wgt = (c0*(1-zd) + c1*zd);
3148
3149 ////
3150
3151 // Denominator errors negligible compared to numerator so do not waste cpu time below.
3152 //Float_t deltaWErr=0;
3153 // Kt
3154 /*
3155 deltaWErr += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3156 // Qo
3157 deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3158 // Qs
3159 deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - minErr)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3160 // Ql
3161 deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - minErr)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3162 */
3163 wgtErr = minErr;
3164
3165
3166}
3167//________________________________________________________________________
3168Float_t AliFourPion::MCWeight(Int_t c[2], Float_t R, Float_t fcSq, Float_t qinv, Float_t k12){
3169
3170 Float_t radius = R/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3171 Float_t r12=radius*(1-k12/2.0);
3172 SetFSIindex(R);
3173 Float_t coulCorr12 = FSICorrelation(c[0], c[1], qinv);
3174 if(c[0]==c[1]){
3175 Float_t arg=qinv*r12;
3176 Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3177 EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3178 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv*r12,2))*pow(EW,2))*coulCorr12);
3179 }else {
3180 return ((1-fcSq) + fcSq*coulCorr12);
3181 }
3182
3183}
3184//________________________________________________________________________
3185Float_t AliFourPion::MCWeightOSL(Int_t charge1, Int_t charge2, Int_t r, Int_t dampIndex, Float_t qinv, Float_t qo, Float_t qs, Float_t ql){
3186
3187 Float_t radiusOut = Float_t(r)/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3188 Float_t radiusSide = radiusOut;
3189 Float_t radiusLong = radiusOut;
3190 Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3191 Float_t coulCorr12 = FSICorrelation(charge1, charge2, qinv);
3192 if(charge1==charge2){
3193 return ((1-myDamp) + myDamp*(1 + exp(-pow(qo*radiusOut,2)) * exp(-pow(qs*radiusSide,2)) * exp(-pow(ql*radiusLong,2)))*coulCorr12);
3194 }else {
3195 return ((1-myDamp) + myDamp*coulCorr12);
3196 }
3197
3198}
3199
3200//________________________________________________________________________
3201Float_t AliFourPion::MCWeight3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3], Float_t kT[3]){
3202 // FSI + QS correlations
3203 if(term==5) return 1.0;
3204
3205 Float_t radius=R/0.19733;
3206 Float_t r12=radius*(1-kT[0]/2.0);
3207 Float_t r13=radius*(1-kT[1]/2.0);
3208 Float_t r23=radius*(1-kT[2]/2.0);
3209
3210 Float_t fc = sqrt(fcSq);
3211
3212 SetFSIindex(R);
3213 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3214 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3215 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3216
3217 if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3218 Float_t arg12=qinv[0]*r12;
3219 Float_t arg13=qinv[1]*r13;
3220 Float_t arg23=qinv[2]*r23;
3221 Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
3222 EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
3223 Float_t EW13 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
3224 EW13 += kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
3225 Float_t EW23 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
3226 EW23 += kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
3227 if(term==1){
3228 Float_t C3QS = 1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2) + exp(-pow(qinv[1]*r13,2))*pow(EW13,2) + exp(-pow(qinv[2]*r23,2))*pow(EW23,2);
3229 C3QS += 2*exp(-(pow(r12,2)*pow(qinv[0],2) + pow(r13,2)*pow(qinv[1],2) + pow(r23,2)*pow(qinv[2],2))/2.)*EW12*EW13*EW23;
3230 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3231 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12;
3232 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13;
3233 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23;
3234 C3 += pow(fc,3)*C3QS*Kfactor12*Kfactor13*Kfactor23;
3235 return C3;
3236 }else if(term==2){
3237 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12);
3238 }else if(term==3){
3239 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13);
3240 }else if(term==4){
3241 return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23);
3242 }else return 1.0;
3243
3244 }else{// mixed charge case
3245 Float_t arg=qinv[0]*r12;
3246 Float_t KfactorSC = Kfactor12;
3247 Float_t KfactorMC1 = Kfactor13;
3248 Float_t KfactorMC2 = Kfactor23;
3249 if(c[0]==c[2]) {arg=qinv[1]*r13; KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;}
3250 if(c[1]==c[2]) {arg=qinv[2]*r23; KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3251 Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3252 EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3253 if(term==1){
3254 Float_t C3QS = 1 + exp(-pow(arg,2))*pow(EW,2);
3255 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3256 C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(arg,2))*pow(EW,2))*KfactorSC;
3257 C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3258 C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3259 C3 += pow(fc,3)*C3QS*KfactorSC*KfactorMC1*KfactorMC2;
3260 return C3;
3261 }else if(term==2){
3262 if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3263 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3264 }else if(term==3){
3265 return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3266 }else if(term==4){
3267 if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3268 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3269 }else return 1.0;
3270 }
3271
3272}
3273//________________________________________________________________________
3274Float_t AliFourPion::MCWeightFSI3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3]){
3275 // FSI only (no QS correlations)
3276 if(term==5) return 1.0;
3277
3278 Float_t fc = sqrt(fcSq);
3279 SetFSIindex(R);
3280 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3281 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3282 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3283
3284 if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3285 if(term==1){
3286 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3287 C3 += pow(fc,2)*(1-fc)*Kfactor12;
3288 C3 += pow(fc,2)*(1-fc)*Kfactor13;
3289 C3 += pow(fc,2)*(1-fc)*Kfactor23;
3290 C3 += pow(fc,3)*Kfactor12*Kfactor13*Kfactor23;
3291 return C3;
3292 }else if(term==2){
3293 return ((1-fcSq) + fcSq*Kfactor12);
3294 }else if(term==3){
3295 return ((1-fcSq) + fcSq*Kfactor13);
3296 }else if(term==4){
3297 return ((1-fcSq) + fcSq*Kfactor23);
3298 }else return 1.0;
3299
3300 }else{// mixed charge case
3301 Float_t KfactorSC = Kfactor12;
3302 Float_t KfactorMC1 = Kfactor13;
3303 Float_t KfactorMC2 = Kfactor23;
3304 if(c[0]==c[2]) {KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;}
3305 if(c[1]==c[2]) {KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3306 if(term==1){
3307 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3308 C3 += pow(fc,2)*(1-fc)*KfactorSC;
3309 C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3310 C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3311 C3 += pow(fc,3)*KfactorSC*KfactorMC1*KfactorMC2;
3312 return C3;
3313 }else if(term==2){
3314 if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*KfactorSC);
3315 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3316 }else if(term==3){
3317 return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3318 }else if(term==4){
3319 if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*KfactorSC);
3320 else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3321 }else return 1.0;
3322 }
3323
3324}
3325//________________________________________________________________________
3326Float_t AliFourPion::MCWeight4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6], Float_t kT[6]){
6bb6954b 3327 if(term==13) return 1.0;
be9ef9f9 3328
3329 // Charge ordering:
3330 // ----, ---+, --++, -+++, ++++
3331 //
3332 // term ordering:
3333 // Term 1: 1-2-3-4 (all same-event)
3334 // Term 2: 1-2-3 4 (particle 4 from different event)
3335 // Term 3: 1-2-4 3 (particle 3 from different event)
3336 // Term 4: 1-3-4 2 (particle 2 from different event)
3337 // Term 5: 2-3-4 1 (particle 1 from different event)
3338 // Term 6: 1-2 3 4 (particle 1 and 2 from same event)
3339 // Term 7: 1-3 2 4
3340 // Term 8: 1-4 2 3
3341 // Term 9: 2-3 1 4
3342 // Term 10: 2-4 1 3
3343 // Term 11: 3-4 1 2
3344 // Term 12: 1 2 3 4 (all from different events)
3345
3346 Float_t radius = R/0.19733;
3347 Float_t r[6]={0};
3348 r[0]=radius*(1-kT[0]/2.0);
3349 r[1]=radius*(1-kT[1]/2.0);
3350 r[2]=radius*(1-kT[2]/2.0);
3351 r[3]=radius*(1-kT[3]/2.0);
3352 r[4]=radius*(1-kT[4]/2.0);
3353 r[5]=radius*(1-kT[5]/2.0);
3354
3355 Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
3356
3357 Float_t fc = sqrt(fcSq);
3358 SetFSIindex(R);
3359 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3360 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3361 Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
3362 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
3363 Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
3364 Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
3365 Float_t arg12=qinv[0]*r[0];
3366 Float_t arg13=qinv[1]*r[1];
3367 Float_t arg14=qinv[2]*r[2];
3368 Float_t arg23=qinv[3]*r[3];
3369 Float_t arg24=qinv[4]*r[4];
3370 Float_t arg34=qinv[5]*r[5];
3371 // Exchange Amplitudes
3372 Float_t EA12 = exp(-pow(arg12,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12) + kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12));
3373 Float_t EA13 = exp(-pow(arg13,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13) + kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12));
3374 Float_t EA14 = exp(-pow(arg14,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg14,3) - 12.*arg14) + kappa4/(24.*pow(2.,2))*(16.*pow(arg14,4) -48.*pow(arg14,2) + 12));
3375 Float_t EA23 = exp(-pow(arg23,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23) + kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12));
3376 Float_t EA24 = exp(-pow(arg24,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg24,3) - 12.*arg24) + kappa4/(24.*pow(2.,2))*(16.*pow(arg24,4) -48.*pow(arg24,2) + 12));
3377 Float_t EA34 = exp(-pow(arg34,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg34,3) - 12.*arg34) + kappa4/(24.*pow(2.,2))*(16.*pow(arg34,4) -48.*pow(arg34,2) + 12));
3378
3379 if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
3380
3381 if(term==1){
3382 Float_t C4QS = 1 + pow(EA12,2) + pow(EA13,2) + pow(EA14,2) + pow(EA23,2) + pow(EA24,2) + pow(EA34,2);// baseline + single pairs
3383 C4QS += pow(EA12,2) * pow(EA34,2);// 2-pairs
3384 C4QS += pow(EA13,2) * pow(EA24,2);// 2-pairs
3385 C4QS += pow(EA14,2) * pow(EA23,2);// 2-pairs
3386 C4QS += 2*EA12*EA13*EA23 + 2*EA12*EA14*EA24 + 2*EA13*EA14*EA34 + 2*EA23*EA24*EA34;// 3-particle exhange
3387 C4QS += 3*EA12*EA23*EA34*EA14 + 3*EA12*EA13*EA34*EA24;// 4-particle exchange
3388 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3389 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA12,2))*Kfactor12 + (1 + pow(EA13,2))*Kfactor13 + (1 + pow(EA14,2))*Kfactor14 );
3390 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA23,2))*Kfactor23 + (1 + pow(EA24,2))*Kfactor24 + (1 + pow(EA34,2))*Kfactor34);
3391 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA13,2) + pow(EA23,2) + 2*EA12*EA13*EA23) * Kfactor12*Kfactor13*Kfactor23;
3392 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA14,2) + pow(EA24,2) + 2*EA12*EA14*EA24) * Kfactor12*Kfactor14*Kfactor24;
3393 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA13,2) + pow(EA14,2) + pow(EA34,2) + 2*EA13*EA14*EA34) * Kfactor13*Kfactor14*Kfactor34;
3394 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA23,2) + pow(EA24,2) + pow(EA34,2) + 2*EA23*EA24*EA34) * Kfactor23*Kfactor24*Kfactor34;
3395 C4 += pow(fc,4)*C4QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3396 return C4;
3397 }else if(term<=5){
3398 Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0;
3399 if(term==2) {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
3400 else if(term==3) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
3401 else if(term==4) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
3402 else {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
3403 Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2);
3404 C3QS += 2*EA1*EA2*EA3;
3405 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3406 C3 += pow(fc,2)*(1-fc)*( (1+pow(EA1,2))*Kpair1 + (1+pow(EA2,2))*Kpair2 + (1+pow(EA3,2))*Kpair3 );
3407 C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
3408 return C3;
3409 }else if(term<=11){
3410 if(term==6) return ((1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12);
3411 else if(term==7) return ((1-fcSq) + fcSq*(1 + pow(EA13,2))*Kfactor13);
3412 else if(term==8) return ((1-fcSq) + fcSq*(1 + pow(EA14,2))*Kfactor14);
3413 else if(term==9) return ((1-fcSq) + fcSq*(1 + pow(EA23,2))*Kfactor23);
3414 else if(term==10) return ((1-fcSq) + fcSq*(1 + pow(EA24,2))*Kfactor24);
3415 else return ((1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34);
6bb6954b 3416 }else if(term==12){
5591748e 3417 Float_t C22 = (1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12;
3418 C22 *= (1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34;
6bb6954b 3419 return C22;
be9ef9f9 3420 }else return 1.0;
3421
3422 }else{// mixed charge case
3423 if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
3424 Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3425 Int_t c_OddOneOut = 1;
3426 if(ChargeSum==3) c_OddOneOut = 0;
3427 //
3428 if(c[0]==c_OddOneOut) {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
3429 else if(c[1]==c_OddOneOut) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3430 else if(c[2]==c_OddOneOut) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24; Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3431 else {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23; Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3432
3433 if(term==1){
3434 Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
3435 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3436 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 + (1 + pow(EA3,2))*Kpair3 );
3437 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
3438 C4 += pow(fc,3)*(1-fc)*(1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3) * Kpair1*Kpair2*Kpair3;
3439 C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair4*Kpair5 + (1+pow(EA2,2))*Kpair2*Kpair4*Kpair6 + (1+pow(EA3,2))*Kpair3*Kpair5*Kpair6);// doesn't matter which MC K's
3440 C4 += pow(fc,4)*C3QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3441 return C4;
3442 }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
3443 Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
3444 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3445 C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;
3446 C3 += pow(fc,2)*(1-fc)*(1+pow(EA2,2))*Kpair2;
3447 C3 += pow(fc,2)*(1-fc)*(1+pow(EA3,2))*Kpair3;
3448 C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
3449 return C3;
3450 }else if(term<=5){// one SC pair, two MC pairs
3451 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3452 C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
3453 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3454 C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
3455 C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair5;
3456 return C3;
3457 }else if(term==6 || term==7){
3458 if(ChargeSum==1) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3459 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3460 }else if(term==8){
3461 return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3462 }else if(term==9){
3463 return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
6bb6954b 3464 }else if(term==10 || term==11){
be9ef9f9 3465 if(ChargeSum==3) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
3466 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
6bb6954b 3467 }else return 1.0;// for 12 and 13
be9ef9f9 3468 }else{// --++ configuration
3469 Float_t EA1=0, EA2=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3470 if(c[0]==c[1]) {EA1=EA12; EA2=EA34; Kpair1=Kfactor12; Kpair2=Kfactor34; Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3471 else if(c[0]==c[2]) {EA1=EA13; EA2=EA24; Kpair1=Kfactor13; Kpair2=Kfactor24; Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3472 else {EA1=EA14; EA2=EA23; Kpair1=Kfactor14; Kpair2=Kfactor23; Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3473 //
3474 if(term==1){
3475 Float_t C2QS = 1 + pow(EA1,2)*pow(EA2,2);
3476 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3477 C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 );
3478 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair3 + Kpair4 + Kpair5 + Kpair6 );
3479 C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair3*Kpair4 + (1 + pow(EA2,2))*Kpair2*Kpair3*Kpair4);
3480 C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair5*Kpair6 + (1 + pow(EA2,2))*Kpair2*Kpair5*Kpair6);// doesn't matter which two MC K's used
3481 C4 += pow(fc,4)*C2QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3482 return C4;
3483 }else if(term<=5){
3484 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3485 C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
3486 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3487 C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
3488 C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair6;
3489 return C3;
3490 }else if(term==6 || term==11){
3491 return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
6bb6954b 3492 }else if(term !=12 && term !=13){
be9ef9f9 3493 return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
3494 }else return 1.0;
3495 }
3496 }
3497
3498}
3499//________________________________________________________________________
3500Float_t AliFourPion::MCWeightFSI4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6]){
6bb6954b 3501 if(term==13) return 1.0;
be9ef9f9 3502
3503 Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
3504
3505 Float_t fc = sqrt(fcSq);
3506 SetFSIindex(R);
3507 Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3508 Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3509 Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
3510 Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
3511 Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
3512 Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
3513
3514 if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
3515
3516 if(term==1){
3517 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3518 C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor12 + Kfactor13 + Kfactor14 );
3519 C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor23 + Kfactor24 + Kfactor34 );
3520 C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor13*Kfactor23;
3521 C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor14*Kfactor24;
3522 C4 += pow(fc,3)*(1-fc) * Kfactor13*Kfactor14*Kfactor34;
3523 C4 += pow(fc,3)*(1-fc) * Kfactor23*Kfactor24*Kfactor34;
3524 C4 += pow(fc,4) * Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3525 return C4;
3526 }else if(term<=5){
3527 Float_t Kpair1=0, Kpair2=0, Kpair3=0;
3528 if(term==2) {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
3529 else if(term==3) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
3530 else if(term==4) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
3531 else {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
3532 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3533 C3 += pow(fc,2)*(1-fc)*( Kpair1 + Kpair2 + Kpair3 );
3534 C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
3535 return C3;
3536 }else if(term<=11){
3537 if(term==6) return ((1-fcSq) + fcSq*Kfactor12);
3538 else if(term==7) return ((1-fcSq) + fcSq*Kfactor13);
3539 else if(term==8) return ((1-fcSq) + fcSq*Kfactor14);
3540 else if(term==9) return ((1-fcSq) + fcSq*Kfactor23);
3541 else if(term==10) return ((1-fcSq) + fcSq*Kfactor24);
3542 else return ((1-fcSq) + fcSq*Kfactor34);
6bb6954b 3543 }else if(term==12){
5591748e 3544 Float_t C22 = (1-fcSq) + fcSq*Kfactor12;
3545 C22 *= (1-fcSq) + fcSq*Kfactor34;
6bb6954b 3546 return C22;
be9ef9f9 3547 }else return 1.0;
3548
3549 }else{// mixed charge case
3550 if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
3551 Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3552 Int_t c_OddOneOut = 1;
3553 if(ChargeSum==3) c_OddOneOut = 0;
3554 //
3555 if(c[0]==c_OddOneOut) {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
3556 else if(c[1]==c_OddOneOut) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34; Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3557 else if(c[2]==c_OddOneOut) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24; Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3558 else {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23; Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3559
3560 if(term==1){
3561 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3562 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 );
3563 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
3564 C4 += pow(fc,3)*(1-fc)*Kpair1*Kpair2*Kpair3;
3565 C4 += pow(fc,3)*(1-fc)*(Kpair1*Kpair4*Kpair5 + Kpair2*Kpair4*Kpair6 + Kpair3*Kpair5*Kpair6);// doesn't matter which two MC K's used
3566 C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3567 return C4;
3568 }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
3569 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3570 C3 += pow(fc,2)*(1-fc)*Kpair1;
3571 C3 += pow(fc,2)*(1-fc)*Kpair2;
3572 C3 += pow(fc,2)*(1-fc)*Kpair3;
3573 C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
3574 return C3;
3575 }else if(term<=5){// one SC pair, two MC pairs
3576 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3577 C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
3578 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3579 C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
3580 C3 += pow(fc,3)*Kpair1*Kpair4*Kpair5;
3581 return C3;
3582 }else if(term==6 || term==7){
3583 if(ChargeSum==1) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3584 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3585 }else if(term==8){
3586 return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
3587 }else if(term==9){
3588 return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
6bb6954b 3589 }else if(term==10 || term==11){
be9ef9f9 3590 if(ChargeSum==3) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
3591 else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
6bb6954b 3592 }else return 1.0;// 12 and 13
be9ef9f9 3593 }else{// --++ configuration
3594 Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
3595 if(c[0]==c[1]) {Kpair1=Kfactor12; Kpair2=Kfactor34; Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
3596 else if(c[0]==c[2]) {Kpair1=Kfactor13; Kpair2=Kfactor24; Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
3597 else {Kpair1=Kfactor14; Kpair2=Kfactor23; Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
3598 //
3599 if(term==1){
3600 Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
3601 C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 + Kpair4 + Kpair5 + Kpair6);
3602 C4 += pow(fc,3)*(1-fc)*( Kpair1*Kpair3*Kpair4 + Kpair2*Kpair3*Kpair4 + Kpair1*Kpair5*Kpair6 + Kpair2*Kpair5*Kpair6);
3603 C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
3604 return C4;
3605 }else if(term<=5){
3606 Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3607 C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
3608 C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
3609 C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
3610 C3 += pow(fc,3)*Kpair1*Kpair4*Kpair6;
3611 return C3;
3612 }else if(term==6 || term==11){
3613 return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
6bb6954b 3614 }else if(term !=12 && term !=13){
be9ef9f9 3615 return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
3616 }else return 1.0;
3617 }
3618 }
3619
3620}
3621//________________________________________________________________________
3622void AliFourPion::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D){
3623
3624
3625 if(legoCase){
3626 cout<<"LEGO call to SetMomResCorrections"<<endl;
3627 fMomResC2 = (TH2D*)temp2D->Clone();
3628 fMomResC2->SetDirectory(0);
3629 }else {
3630 TFile *momResFile = new TFile("MomResFile.root","READ");
3631 if(!momResFile->IsOpen()) {
3632 cout<<"No momentum resolution file found"<<endl;
3633 AliFatal("No momentum resolution file found. Kill process.");
3634 }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
3635
3636 TH2D *temp2D2 = (TH2D*)momResFile->Get("MRC_C2_SC");
3637 fMomResC2 = (TH2D*)temp2D2->Clone();
3638 fMomResC2->SetDirectory(0);
3639
3640 momResFile->Close();
3641 }
3642
3643
3644 for(Int_t bx=1; bx<=fMomResC2->GetNbinsX(); bx++){
3645 for(Int_t by=1; by<=fMomResC2->GetNbinsY(); by++){
3646 if(fMomResC2->GetBinContent(bx,by) > 1.5) fMomResC2->SetBinContent(bx,by, 1.5);// Maximum is ~1.02
3647 if(fMomResC2->GetBinContent(bx,by) < 0.95) fMomResC2->SetBinContent(bx,by, 0.95);// Minimum is ~0.98
3648 }
3649 }
3650
3651 cout<<"Done reading momentum resolution file"<<endl;
3652}
3653//________________________________________________________________________
3654void AliFourPion::SetFSICorrelations(Bool_t legoCase, TH1D *tempss[12], TH1D *tempos[12]){
3655 // read in 2-particle and 3-particle FSI correlations = K2 & K3
3656 // 2-particle input histo from file is binned in qinv. 3-particle in qinv of each pair
3657 if(legoCase){
3658 cout<<"LEGO call to SetFSICorrelations"<<endl;
3659 for(Int_t MB=0; MB<12; MB++) {
3660 fFSIss[MB] = (TH1D*)tempss[MB]->Clone();
3661 fFSIos[MB] = (TH1D*)tempos[MB]->Clone();
3662 //
3663 fFSIss[MB]->SetDirectory(0);
3664 fFSIos[MB]->SetDirectory(0);
3665 }
3666 }else {
3667 cout<<"non LEGO call to SetFSICorrelations"<<endl;
3668 TFile *fsifile = new TFile("KFile.root","READ");
3669 if(!fsifile->IsOpen()) {
3670 cout<<"No FSI file found"<<endl;
3671 AliFatal("No FSI file found. Kill process.");
3672 }else {cout<<"Good FSI File Found!"<<endl;}
3673
3674 TH1D *temphistoSS[12];
3675 TH1D *temphistoOS[12];
3676 for(Int_t MB=0; MB<12; MB++) {
3677 TString *nameK2SS = new TString("K2ss_");
3678 *nameK2SS += MB;
3679 temphistoSS[MB] = (TH1D*)fsifile->Get(nameK2SS->Data());
3680 //
3681 TString *nameK2OS = new TString("K2os_");
3682 *nameK2OS += MB;
3683 temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data());
3684 //
3685 fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone();
3686 fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone();
3687 fFSIss[MB]->SetDirectory(0);
3688 fFSIos[MB]->SetDirectory(0);
3689 }
3690 //
3691
3692 fsifile->Close();
3693 }
3694
3695 cout<<"Done reading FSI file"<<endl;
3696}
3697//________________________________________________________________________
3698Float_t AliFourPion::FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){
3699 // returns 2-particle Coulomb correlations = K2
3700 Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.);
3701 Int_t qbinH = qbinL+1;
3702 if(qbinL <= 0) return 1.0;
3703 if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) return 1.0;
3704
3705 Float_t slope=0;
3706 if(charge1==charge2){
3707 slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH);
3708 slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
3709 return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL));
3710 }else {
3711 slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH);
3712 slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
3713 return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL));
3714 }
3715}
3716//________________________________________________________________________
3717void AliFourPion::SetFillBins2(Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
3718 if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
3719 else {b1=c1; b2=c2;}
3720}
3721//________________________________________________________________________
3722void AliFourPion::SetFillBins3(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){
3723
3724 // "part" specifies which pair is from the same event. Only relevant for terms 2-4
3725 Bool_t seSS=kFALSE;
3726 if(part==1) {// default case (irrelevant for term 1 and term 5)
3727 if(c1==c2) seSS=kTRUE;
3728 }
3729 if(part==2){
3730 if(c1==c3) seSS=kTRUE;
3731 }
3732
3733
3734 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3735 if( (c1+c2+c3)==1) {
3736 b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
3737 //
3738 if(seSS) fill2=kTRUE;
3739 else {fill3=kTRUE; fill4=kTRUE;}
3740 //
3741 }else if( (c1+c2+c3)==2) {
3742 b1=0; b2=1; b3=1;
3743 //
3744 if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
3745 else fill4=kTRUE;
3746 //
3747 }else {
3748 b1=c1; b2=c2; b3=c3;
3749 fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
3750 }
3751
3752}
3753//________________________________________________________________________
6bb6954b 3754void AliFourPion::SetFillBins4(Int_t c1, Int_t c2, Int_t c3, Int_t c4, Int_t &b1, Int_t &b2, Int_t &b3, Int_t &b4, Int_t ENsum, Bool_t fillTerm[13]){
be9ef9f9 3755
3756 // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
3757 if( (c1+c2+c3+c4)==0 || (c1+c2+c3+c4)==4) {// all of the same charge: ---- or ++++
3758
3759 b1=c1; b2=c2; b3=c3; b4=c4;
3760 if(ENsum==0) fillTerm[0]=kTRUE;
3761 else if(ENsum==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
6bb6954b 3762 else if(ENsum==2) {fillTerm[11]=kTRUE;}
be9ef9f9 3763 else if(ENsum==3) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
6bb6954b 3764 else fillTerm[12]=kTRUE;
be9ef9f9 3765
3766 }else if( (c1+c2+c3+c4)==1) {// one positive charge: ---+
3767
3768 b1=0; b2=0; b3=0; b4=1;// Re-assign to merge degenerate histos
3769 if(ENsum==0) fillTerm[0]=kTRUE;
3770 else if(ENsum==1){
3771 if(c4==1) fillTerm[1]=kTRUE;
3772 else {fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
6bb6954b 3773 }else if(ENsum==2){
3774 fillTerm[11]=kTRUE;
be9ef9f9 3775 }else if(ENsum==3){
3776 if(c3==1 || c4==1) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[8]=kTRUE;}
3777 else {fillTerm[7]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
6bb6954b 3778 }else fillTerm[12]=kTRUE;
be9ef9f9 3779
3780 }else if( (c1+c2+c3+c4)==2) {// two positive charges: --++
3781
3782 b1=0; b2=0; b3=1; b4=1;// Re-assign to merge degenerate histos
3783 if(ENsum==0) fillTerm[0]=kTRUE;
3784 else if(ENsum==1){
3785 if(c4==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE;}
3786 else {fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
6bb6954b 3787 }else if(ENsum==2){
3788 if( (c1+c2)==0) fillTerm[11]=kTRUE;
be9ef9f9 3789 }else if(ENsum==3){
3790 if( (c1+c2)==0) fillTerm[5]=kTRUE;
3791 else if( (c1+c2)==1) {fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE;}
3792 else fillTerm[10]=kTRUE;
6bb6954b 3793 }else fillTerm[12]=kTRUE;
be9ef9f9 3794
3795 }else{// three positive charges
3796
3797 b1=0; b2=1; b3=1; b4=1;// Re-assign to merge degenerate histos
3798 if(ENsum==0) fillTerm[0]=kTRUE;
3799 else if(ENsum==1){
3800 if(c4==0) fillTerm[4]=kTRUE;
3801 else {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE;}
6bb6954b 3802 }else if(ENsum==2){
3803 fillTerm[11]=kTRUE;
be9ef9f9 3804 }else if(ENsum==3){
3805 if(c3==0 || c4==0) {fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
3806 else {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE;}
6bb6954b 3807 }else fillTerm[12]=kTRUE;
be9ef9f9 3808
3809 }
3810
3811}
3812//________________________________________________________________________
3813void AliFourPion::SetFSIindex(Float_t R){
3814 if(!fMCcase){
3815 if(fPbPbcase){
3816 if(fMbin==0) fFSIindex = 0;//0-5%
3817 else if(fMbin==1) fFSIindex = 1;//5-10%
3818 else if(fMbin<=3) fFSIindex = 2;//10-20%
3819 else if(fMbin<=5) fFSIindex = 3;//20-30%
3820 else if(fMbin<=7) fFSIindex = 4;//30-40%
3821 else if(fMbin<=9) fFSIindex = 5;//40-50%
3822 else if(fMbin<=12) fFSIindex = 6;//40-50%
3823 else if(fMbin<=15) fFSIindex = 7;//40-50%
3824 else if(fMbin<=18) fFSIindex = 8;//40-50%
3825 else fFSIindex = 8;//90-100%
3826 }else fFSIindex = 9;// pp and pPb
3827 }else{// FSI binning for MC
3828 if(R>=10.) fFSIindex = 0;
3829 else if(R>=9.) fFSIindex = 1;
3830 else if(R>=8.) fFSIindex = 2;
3831 else if(R>=7.) fFSIindex = 3;
3832 else if(R>=6.) fFSIindex = 4;
3833 else if(R>=5.) fFSIindex = 5;
3834 else if(R>=4.) fFSIindex = 6;
3835 else if(R>=3.) fFSIindex = 7;
3836 else if(R>=2.) fFSIindex = 8;
3837 else fFSIindex = 9;
3838 }
3839}
5591748e 3840//________________________________________________________________________
3841void AliFourPion::SetMuonCorrections(Bool_t legoCase, TH2D *tempMuon){
3842 if(legoCase){
3843 cout<<"LEGO call to SetMuonCorrections"<<endl;
3844 fWeightmuonCorrection = (TH2D*)tempMuon->Clone();
3845 fWeightmuonCorrection->SetDirectory(0);
3846 }else {
3847 cout<<"non LEGO call to SetMuonCorrections"<<endl;
3848 TFile *MuonFile=new TFile("MuonCorrection.root","READ");
3849 if(!MuonFile->IsOpen()) {
3850 cout<<"No Muon file found"<<endl;
3851 AliFatal("No Muon file found. Kill process.");
3852 }else {cout<<"Good Muon File Found!"<<endl;}
3853
3854 fWeightmuonCorrection = (TH2D*)MuonFile->Get("WeightmuonCorrection");
3855 fWeightmuonCorrection->SetDirectory(0);
3856 //
3857 MuonFile->Close();
3858 }
3859 cout<<"Done reading Muon file"<<endl;
3860}