Apply pT cuts only to correlation part of analysis. Add weights for MC runs with...
[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
42802bba 37#define FmToGeV 0.19733 // conversion of Fm to GeV
fcb16ac9 38#define kappa3 0.15 // kappa3 Edgeworth coefficient (non-Gaussian features of C2)
39#define kappa4 0.32 // kappa4 Edgeworth coefficient (non-Gaussian features of C2)
42802bba 40#define kappa3Fit 0.1 // kappa3 for c4QS fit
41#define kappa4Fit 0.5 // kappa4 for c4QS fit
be9ef9f9 42
43// Author: Dhevan Gangadharan
44
45ClassImp(AliFourPion)
46
47//________________________________________________________________________
48AliFourPion::AliFourPion():
49AliAnalysisTaskSE(),
50 fname(0),
51 fAOD(0x0),
52 fOutputList(0x0),
53 fPIDResponse(0x0),
54 fEC(0x0),
55 fEvt(0x0),
56 fTempStruct(0x0),
57 fRandomNumber(0x0),
58 fLEGO(kTRUE),
59 fMCcase(kFALSE),
60 fAODcase(kTRUE),
61 fPbPbcase(kTRUE),
62 fGenerateSignal(kFALSE),
63 fGeneratorOnly(kFALSE),
be9ef9f9 64 fTabulatePairs(kFALSE),
95567836 65 fLinearInterpolation(kTRUE),
42802bba 66 fMixedChargeCut(kFALSE),
be9ef9f9 67 fRMax(11),
68 ffcSq(0.7),
80507acf 69 ffcSqMRC(0.6),
be9ef9f9 70 fFilterBit(7),
71 fMaxChi2NDF(10),
72 fMinTPCncls(0),
73 fBfield(0),
74 fMbin(0),
75 fFSIindex(0),
76 fEDbin(0),
77 fMbins(fCentBins),
78 fMultLimit(0),
79 fCentBinLowLimit(0),
80 fCentBinHighLimit(1),
81 fEventCounter(0),
82 fEventsToMix(0),
83 fZvertexBins(0),
84 fMultLimits(),
42802bba 85 fMinPt(0.16),
86 fMaxPt(1.0),
be9ef9f9 87 fQcut(0),
88 fQLowerCut(0),
89 fNormQcutLow(0),
90 fNormQcutHigh(0),
91 fKupperBound(0),
92 fQupperBoundQ2(0),
93 fQupperBoundQ3(0),
94 fQupperBoundQ4(0),
95 fQbinsQ2(1),
96 fQbinsQ3(1),
97 fQbinsQ4(1),
98 fQupperBoundWeights(0),
99 fKstepT(),
100 fKstepY(),
101 fKmeanT(),
102 fKmeanY(),
103 fKmiddleT(),
104 fKmiddleY(),
105 fQstep(0),
106 fQstepWeights(0),
107 fQmean(),
108 fDampStart(0),
109 fDampStep(0),
110 fTPCTOFboundry(0),
111 fTOFboundry(0),
112 fSigmaCutTPC(2.0),
113 fSigmaCutTOF(2.0),
114 fMinSepPairEta(0.03),
115 fMinSepPairPhi(0.04),
116 fShareQuality(0),
117 fShareFraction(0),
118 fTrueMassP(0),
119 fTrueMassPi(0),
120 fTrueMassK(0),
121 fTrueMassKs(0),
122 fTrueMassLam(0),
123 fKtIndexL(0),
124 fKtIndexH(0),
125 fQoIndexL(0),
126 fQoIndexH(0),
127 fQsIndexL(0),
128 fQsIndexH(0),
129 fQlIndexL(0),
130 fQlIndexH(0),
131 fDummyB(0),
132 fKT3transition(0.3),
133 fKT4transition(0.3),
95567836 134 farrP1(),
135 farrP2(),
b71263d0 136 fDefaultsCharSwitch(),
137 fLowQPairSwitch_E0E0(),
138 fLowQPairSwitch_E0E1(),
139 fLowQPairSwitch_E0E2(),
140 fLowQPairSwitch_E0E3(),
6bb6954b 141 fLowQPairSwitch_E1E1(),
b71263d0 142 fLowQPairSwitch_E1E2(),
143 fLowQPairSwitch_E1E3(),
144 fLowQPairSwitch_E2E3(),
145 fNormQPairSwitch_E0E0(),
146 fNormQPairSwitch_E0E1(),
147 fNormQPairSwitch_E0E2(),
148 fNormQPairSwitch_E0E3(),
6bb6954b 149 fNormQPairSwitch_E1E1(),
b71263d0 150 fNormQPairSwitch_E1E2(),
151 fNormQPairSwitch_E1E3(),
152 fNormQPairSwitch_E2E3(),
80507acf 153 fMomResC2SC(0x0),
154 fMomResC2MC(0x0),
5591748e 155 fWeightmuonCorrection(0x0)
be9ef9f9 156{
157 // Default constructor
158 for(Int_t mb=0; mb<fMbins; mb++){
159 for(Int_t edB=0; edB<fEDbins; edB++){
160 for(Int_t c1=0; c1<2; c1++){
161 for(Int_t c2=0; c2<2; c2++){
162 for(Int_t term=0; term<2; term++){
163
164 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
165
166 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
167 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
168 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
169 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
170 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
171 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
172
173 }// term_2
174
175
176 for(Int_t c3=0; c3<2; c3++){
177 for(Int_t term=0; term<5; term++){
178
179 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
180 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
b71263d0 181 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
be9ef9f9 182 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
183
184 }// term_3
b71263d0 185
186 for(Int_t c4=0; c4<2; c4++){
6bb6954b 187 for(Int_t term=0; term<13; term++){
b71263d0 188
189 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
190 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
191 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
192 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
193
194 }// term_4
195
196 }// c4
be9ef9f9 197 }//c3
198 }//c2
199 }//c1
200 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
201 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
202 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
203 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
204 }
205 }
206
207 }// ED
208 }// Mbin
209
210 // Initialize FSI histograms
211 for(Int_t i=0; i<12; i++){
212 fFSIss[i]=0x0;
213 fFSIos[i]=0x0;
214 }
215
216
217 // Initialize fNormWeight and fNormWeightErr to 0
218 for(Int_t i=0; i<3; i++){// Kt iterator
219 for(Int_t j=0; j<10; j++){// Mbin iterator
220 fNormWeight[i][j]=0x0;
221 }
222 }
223
224
225}
226//________________________________________________________________________
227AliFourPion::AliFourPion(const Char_t *name)
228: AliAnalysisTaskSE(name),
229 fname(name),
230 fAOD(0x0),
231 fOutputList(0x0),
232 fPIDResponse(0x0),
233 fEC(0x0),
234 fEvt(0x0),
235 fTempStruct(0x0),
236 fRandomNumber(0x0),
237 fLEGO(kTRUE),
238 fMCcase(kFALSE),
239 fAODcase(kTRUE),
240 fPbPbcase(kTRUE),
241 fGenerateSignal(kFALSE),
242 fGeneratorOnly(kFALSE),
be9ef9f9 243 fTabulatePairs(kFALSE),
95567836 244 fLinearInterpolation(kTRUE),
42802bba 245 fMixedChargeCut(kFALSE),
be9ef9f9 246 fRMax(11),
247 ffcSq(0.7),
80507acf 248 ffcSqMRC(0.6),
be9ef9f9 249 fFilterBit(7),
250 fMaxChi2NDF(10),
251 fMinTPCncls(0),
252 fBfield(0),
253 fMbin(0),
254 fFSIindex(0),
255 fEDbin(0),
256 fMbins(fCentBins),
257 fMultLimit(0),
258 fCentBinLowLimit(0),
259 fCentBinHighLimit(1),
260 fEventCounter(0),
261 fEventsToMix(0),
262 fZvertexBins(0),
263 fMultLimits(),
42802bba 264 fMinPt(0.16),
265 fMaxPt(1.0),
be9ef9f9 266 fQcut(0),
267 fQLowerCut(0),
268 fNormQcutLow(0),
269 fNormQcutHigh(0),
270 fKupperBound(0),
271 fQupperBoundQ2(0),
272 fQupperBoundQ3(0),
273 fQupperBoundQ4(0),
274 fQbinsQ2(1),
275 fQbinsQ3(1),
276 fQbinsQ4(1),
277 fQupperBoundWeights(0),
278 fKstepT(),
279 fKstepY(),
280 fKmeanT(),
281 fKmeanY(),
282 fKmiddleT(),
283 fKmiddleY(),
284 fQstep(0),
285 fQstepWeights(0),
286 fQmean(),
287 fDampStart(0),
288 fDampStep(0),
289 fTPCTOFboundry(0),
290 fTOFboundry(0),
291 fSigmaCutTPC(2.0),
292 fSigmaCutTOF(2.0),
293 fMinSepPairEta(0.03),
294 fMinSepPairPhi(0.04),
295 fShareQuality(0),
296 fShareFraction(0),
297 fTrueMassP(0),
298 fTrueMassPi(0),
299 fTrueMassK(0),
300 fTrueMassKs(0),
301 fTrueMassLam(0),
302 fKtIndexL(0),
303 fKtIndexH(0),
304 fQoIndexL(0),
305 fQoIndexH(0),
306 fQsIndexL(0),
307 fQsIndexH(0),
308 fQlIndexL(0),
309 fQlIndexH(0),
310 fDummyB(0),
311 fKT3transition(0.3),
312 fKT4transition(0.3),
95567836 313 farrP1(),
314 farrP2(),
b71263d0 315 fDefaultsCharSwitch(),
316 fLowQPairSwitch_E0E0(),
317 fLowQPairSwitch_E0E1(),
318 fLowQPairSwitch_E0E2(),
319 fLowQPairSwitch_E0E3(),
6bb6954b 320 fLowQPairSwitch_E1E1(),
b71263d0 321 fLowQPairSwitch_E1E2(),
322 fLowQPairSwitch_E1E3(),
323 fLowQPairSwitch_E2E3(),
324 fNormQPairSwitch_E0E0(),
325 fNormQPairSwitch_E0E1(),
326 fNormQPairSwitch_E0E2(),
327 fNormQPairSwitch_E0E3(),
6bb6954b 328 fNormQPairSwitch_E1E1(),
b71263d0 329 fNormQPairSwitch_E1E2(),
330 fNormQPairSwitch_E1E3(),
331 fNormQPairSwitch_E2E3(),
80507acf 332 fMomResC2SC(0x0),
333 fMomResC2MC(0x0),
5591748e 334 fWeightmuonCorrection(0x0)
be9ef9f9 335{
336 // Main constructor
337 fAODcase=kTRUE;
6432c081 338
be9ef9f9 339
340
341 for(Int_t mb=0; mb<fMbins; mb++){
342 for(Int_t edB=0; edB<fEDbins; edB++){
343 for(Int_t c1=0; c1<2; c1++){
344 for(Int_t c2=0; c2<2; c2++){
345 for(Int_t term=0; term<2; term++){
346
347 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
348
349 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
350 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
351 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
352 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
353 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
354 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
355
356 }// term_2
357
358 for(Int_t c3=0; c3<2; c3++){
359 for(Int_t term=0; term<5; term++){
360
361 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
362 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
b71263d0 363 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = 0x0;
be9ef9f9 364 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
365
366 }// term_3
b71263d0 367
368 for(Int_t c4=0; c4<2; c4++){
6bb6954b 369 for(Int_t term=0; term<13; term++){
b71263d0 370
371 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = 0x0;
372 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = 0x0;
373 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = 0x0;
374 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = 0x0;
375
376 }// term_4
377 }// c4
be9ef9f9 378 }//c3
379 }//c2
380 }//c1
381
382 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
383 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
384 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
385 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
386 }
387 }
388
389 }// ED
390 }// Mbin
391
392 // Initialize FSI histograms
393 for(Int_t i=0; i<12; i++){
394 fFSIss[i]=0x0;
395 fFSIos[i]=0x0;
396 }
397
398 // Initialize fNormWeight and fNormWeightErr to 0
399 for(Int_t i=0; i<3; i++){// Kt iterator
400 for(Int_t j=0; j<10; j++){// Mbin iterator
401 fNormWeight[i][j]=0x0;
402 }
403 }
404
405
406 DefineOutput(1, TList::Class());
407}
408//________________________________________________________________________
409AliFourPion::AliFourPion(const AliFourPion &obj)
410 : AliAnalysisTaskSE(obj.fname),
411 fname(obj.fname),
412 fAOD(obj.fAOD),
413 //fESD(obj.fESD),
414 fOutputList(obj.fOutputList),
415 fPIDResponse(obj.fPIDResponse),
416 fEC(obj.fEC),
417 fEvt(obj.fEvt),
418 fTempStruct(obj.fTempStruct),
419 fRandomNumber(obj.fRandomNumber),
420 fLEGO(obj.fLEGO),
421 fMCcase(obj.fMCcase),
422 fAODcase(obj.fAODcase),
423 fPbPbcase(obj.fPbPbcase),
424 fGenerateSignal(obj.fGenerateSignal),
425 fGeneratorOnly(obj.fGeneratorOnly),
be9ef9f9 426 fTabulatePairs(obj.fTabulatePairs),
95567836 427 fLinearInterpolation(obj.fLinearInterpolation),
42802bba 428 fMixedChargeCut(obj.fMixedChargeCut),
be9ef9f9 429 fRMax(obj.fRMax),
430 ffcSq(obj.ffcSq),
80507acf 431 ffcSqMRC(obj.ffcSqMRC),
be9ef9f9 432 fFilterBit(obj.fFilterBit),
433 fMaxChi2NDF(obj.fMaxChi2NDF),
434 fMinTPCncls(obj.fMinTPCncls),
435 fBfield(obj.fBfield),
436 fMbin(obj.fMbin),
437 fFSIindex(obj.fFSIindex),
438 fEDbin(obj.fEDbin),
439 fMbins(obj.fMbins),
440 fMultLimit(obj.fMultLimit),
441 fCentBinLowLimit(obj.fCentBinLowLimit),
442 fCentBinHighLimit(obj.fCentBinHighLimit),
443 fEventCounter(obj.fEventCounter),
444 fEventsToMix(obj.fEventsToMix),
445 fZvertexBins(obj.fZvertexBins),
446 fMultLimits(),
42802bba 447 fMinPt(obj.fMinPt),
448 fMaxPt(obj.fMaxPt),
449 fQcut(obj.fQcut),
be9ef9f9 450 fQLowerCut(obj.fQLowerCut),
451 fNormQcutLow(0),
452 fNormQcutHigh(0),
453 fKupperBound(obj.fKupperBound),
454 fQupperBoundQ2(obj.fQupperBoundQ2),
455 fQupperBoundQ3(obj.fQupperBoundQ3),
456 fQupperBoundQ4(obj.fQupperBoundQ4),
457 fQbinsQ2(obj.fQbinsQ2),
458 fQbinsQ3(obj.fQbinsQ3),
459 fQbinsQ4(obj.fQbinsQ4),
460 fQupperBoundWeights(obj.fQupperBoundWeights),
461 fKstepT(),
462 fKstepY(),
463 fKmeanT(),
464 fKmeanY(),
465 fKmiddleT(),
466 fKmiddleY(),
467 fQstep(obj.fQstep),
468 fQstepWeights(obj.fQstepWeights),
469 fQmean(),
470 fDampStart(obj.fDampStart),
471 fDampStep(obj.fDampStep),
472 fTPCTOFboundry(obj.fTPCTOFboundry),
473 fTOFboundry(obj.fTOFboundry),
474 fSigmaCutTPC(obj.fSigmaCutTPC),
475 fSigmaCutTOF(obj.fSigmaCutTOF),
476 fMinSepPairEta(obj.fMinSepPairEta),
477 fMinSepPairPhi(obj.fMinSepPairPhi),
478 fShareQuality(obj.fShareQuality),
479 fShareFraction(obj.fShareFraction),
480 fTrueMassP(obj.fTrueMassP),
481 fTrueMassPi(obj.fTrueMassPi),
482 fTrueMassK(obj.fTrueMassK),
483 fTrueMassKs(obj.fTrueMassKs),
484 fTrueMassLam(obj.fTrueMassLam),
485 fKtIndexL(obj.fKtIndexL),
486 fKtIndexH(obj.fKtIndexH),
487 fQoIndexL(obj.fQoIndexL),
488 fQoIndexH(obj.fQoIndexH),
489 fQsIndexL(obj.fQsIndexL),
490 fQsIndexH(obj.fQsIndexH),
491 fQlIndexL(obj.fQlIndexL),
492 fQlIndexH(obj.fQlIndexH),
493 fDummyB(obj.fDummyB),
494 fKT3transition(obj.fKT3transition),
495 fKT4transition(obj.fKT4transition),
95567836 496 farrP1(),
497 farrP2(),
b71263d0 498 fDefaultsCharSwitch(),
499 fLowQPairSwitch_E0E0(),
500 fLowQPairSwitch_E0E1(),
501 fLowQPairSwitch_E0E2(),
502 fLowQPairSwitch_E0E3(),
6bb6954b 503 fLowQPairSwitch_E1E1(),
b71263d0 504 fLowQPairSwitch_E1E2(),
505 fLowQPairSwitch_E1E3(),
506 fLowQPairSwitch_E2E3(),
507 fNormQPairSwitch_E0E0(),
508 fNormQPairSwitch_E0E1(),
509 fNormQPairSwitch_E0E2(),
510 fNormQPairSwitch_E0E3(),
6bb6954b 511 fNormQPairSwitch_E1E1(),
b71263d0 512 fNormQPairSwitch_E1E2(),
513 fNormQPairSwitch_E1E3(),
514 fNormQPairSwitch_E2E3(),
80507acf 515 fMomResC2SC(obj.fMomResC2SC),
516 fMomResC2MC(obj.fMomResC2MC),
5591748e 517 fWeightmuonCorrection(obj.fWeightmuonCorrection)
be9ef9f9 518{
519 // Copy Constructor
520
521 for(Int_t i=0; i<12; i++){
522 fFSIss[i]=obj.fFSIss[i];
523 fFSIos[i]=obj.fFSIos[i];
524 }
525
526 // Initialize fNormWeight and fNormWeightErr to 0
527 for(Int_t i=0; i<3; i++){// Kt iterator
528 for(Int_t j=0; j<10; j++){// Mbin iterator
529 fNormWeight[i][j]=0x0;
530 }
531 }
532
533
534}
535//________________________________________________________________________
536AliFourPion &AliFourPion::operator=(const AliFourPion &obj)
537{
538 // Assignment operator
539 if (this == &obj)
540 return *this;
541
542 fname = obj.fname;
543 fAOD = obj.fAOD;
544 fOutputList = obj.fOutputList;
545 fPIDResponse = obj.fPIDResponse;
546 fEC = obj.fEC;
547 fEvt = obj.fEvt;
548 fTempStruct = obj.fTempStruct;
549 fRandomNumber = obj.fRandomNumber;
550 fLEGO = fLEGO;
551 fMCcase = obj.fMCcase;
552 fAODcase = obj.fAODcase;
553 fPbPbcase = obj.fPbPbcase;
554 fGenerateSignal = obj.fGenerateSignal;
555 fGeneratorOnly = obj.fGeneratorOnly;
be9ef9f9 556 fTabulatePairs = obj.fTabulatePairs;
95567836 557 fLinearInterpolation = obj.fLinearInterpolation;
42802bba 558 fMixedChargeCut = obj.fMixedChargeCut;
be9ef9f9 559 fRMax = obj.fRMax;
560 ffcSq = obj.ffcSq;
80507acf 561 ffcSqMRC = obj.ffcSqMRC;
be9ef9f9 562 fFilterBit = obj.fFilterBit;
563 fMaxChi2NDF = obj.fMaxChi2NDF;
564 fMinTPCncls = obj.fMinTPCncls;
565 fBfield = obj.fBfield;
566 fMbin = obj.fMbin;
567 fFSIindex = obj.fFSIindex;
568 fEDbin = obj.fEDbin;
569 fMbins = obj.fMbins;
570 fMultLimit = obj.fMultLimit;
571 fCentBinLowLimit = obj.fCentBinLowLimit;
572 fCentBinHighLimit = obj.fCentBinHighLimit;
573 fEventCounter = obj.fEventCounter;
574 fEventsToMix = obj.fEventsToMix;
575 fZvertexBins = obj.fZvertexBins;
42802bba 576 fMinPt = obj.fMinPt;
577 fMaxPt = obj.fMaxPt;
578 fQcut = obj.fQcut;
be9ef9f9 579 fQLowerCut = obj.fQLowerCut;
580 fKupperBound = obj.fKupperBound;
581 fQupperBoundQ2 = obj.fQupperBoundQ2;
582 fQupperBoundQ3 = obj.fQupperBoundQ3;
583 fQupperBoundQ4 = obj.fQupperBoundQ4;
584 fQbinsQ2 = obj.fQbinsQ2;
585 fQbinsQ3 = obj.fQbinsQ3;
586 fQbinsQ4 = obj.fQbinsQ4;
587 fQupperBoundWeights = obj.fQupperBoundWeights;
588 fQstep = obj.fQstep;
589 fQstepWeights = obj.fQstepWeights;
590 fDampStart = obj.fDampStart;
591 fDampStep = obj.fDampStep;
592 fTPCTOFboundry = obj.fTPCTOFboundry;
593 fTOFboundry = obj.fTOFboundry;
594 fSigmaCutTPC = obj.fSigmaCutTPC;
595 fSigmaCutTOF = obj.fSigmaCutTOF;
596 fMinSepPairEta = obj.fMinSepPairEta;
597 fMinSepPairPhi = obj.fMinSepPairPhi;
598 fShareQuality = obj.fShareQuality;
599 fShareFraction = obj.fShareFraction;
600 fTrueMassP = obj.fTrueMassP;
601 fTrueMassPi = obj.fTrueMassPi;
602 fTrueMassK = obj.fTrueMassK;
603 fTrueMassKs = obj.fTrueMassKs;
604 fTrueMassLam = obj.fTrueMassLam;
605 fKtIndexL = obj.fKtIndexL;
606 fKtIndexH = obj.fKtIndexH;
607 fQoIndexL = obj.fQoIndexL;
608 fQoIndexH = obj.fQoIndexH;
609 fQsIndexL = obj.fQsIndexL;
610 fQsIndexH = obj.fQsIndexH;
611 fQlIndexL = obj.fQlIndexL;
612 fQlIndexH = obj.fQlIndexH;
613 fDummyB = obj.fDummyB;
614 fKT3transition = obj.fKT3transition;
615 fKT4transition = obj.fKT4transition;
80507acf 616 fMomResC2SC = obj.fMomResC2SC;
617 fMomResC2MC = obj.fMomResC2MC;
5591748e 618 fWeightmuonCorrection = obj.fWeightmuonCorrection;
be9ef9f9 619
620 for(Int_t i=0; i<12; i++){
621 fFSIss[i]=obj.fFSIss[i];
622 fFSIos[i]=obj.fFSIos[i];
623 }
624 for(Int_t i=0; i<3; i++){// Kt iterator
625 for(Int_t j=0; j<10; j++){// Mbin iterator
626 fNormWeight[i][j]=obj.fNormWeight[i][j];
627 }
628 }
629
630 return (*this);
631}
632//________________________________________________________________________
633AliFourPion::~AliFourPion()
634{
635 // Destructor
636 if(fAOD) delete fAOD;
637 //if(fESD) delete fESD;
638 if(fOutputList) delete fOutputList;
639 if(fPIDResponse) delete fPIDResponse;
640 if(fEC) delete fEC;
641 if(fEvt) delete fEvt;
642 if(fTempStruct) delete [] fTempStruct;
643 if(fRandomNumber) delete fRandomNumber;
80507acf 644 if(fMomResC2SC) delete fMomResC2SC;
645 if(fMomResC2MC) delete fMomResC2MC;
5591748e 646 if(fWeightmuonCorrection) delete fWeightmuonCorrection;
647
b71263d0 648 for(Int_t j=0; j<kMultLimitPbPb; j++){
649 if(fLowQPairSwitch_E0E0[j]) delete [] fLowQPairSwitch_E0E0[j];
650 if(fLowQPairSwitch_E0E1[j]) delete [] fLowQPairSwitch_E0E1[j];
651 if(fLowQPairSwitch_E0E2[j]) delete [] fLowQPairSwitch_E0E2[j];
652 if(fLowQPairSwitch_E0E3[j]) delete [] fLowQPairSwitch_E0E3[j];
6bb6954b 653 if(fLowQPairSwitch_E1E1[j]) delete [] fLowQPairSwitch_E1E1[j];
b71263d0 654 if(fLowQPairSwitch_E1E2[j]) delete [] fLowQPairSwitch_E1E2[j];
655 if(fLowQPairSwitch_E1E3[j]) delete [] fLowQPairSwitch_E1E3[j];
656 if(fLowQPairSwitch_E2E3[j]) delete [] fLowQPairSwitch_E2E3[j];
657 //
658 if(fNormQPairSwitch_E0E0[j]) delete [] fNormQPairSwitch_E0E0[j];
659 if(fNormQPairSwitch_E0E1[j]) delete [] fNormQPairSwitch_E0E1[j];
660 if(fNormQPairSwitch_E0E2[j]) delete [] fNormQPairSwitch_E0E2[j];
661 if(fNormQPairSwitch_E0E3[j]) delete [] fNormQPairSwitch_E0E3[j];
6bb6954b 662 if(fNormQPairSwitch_E1E1[j]) delete [] fNormQPairSwitch_E1E1[j];
b71263d0 663 if(fNormQPairSwitch_E1E2[j]) delete [] fNormQPairSwitch_E1E2[j];
664 if(fNormQPairSwitch_E1E3[j]) delete [] fNormQPairSwitch_E1E3[j];
665 if(fNormQPairSwitch_E2E3[j]) delete [] fNormQPairSwitch_E2E3[j];
666 }
be9ef9f9 667
668 //
669 for(Int_t mb=0; mb<fMbins; mb++){
670 for(Int_t edB=0; edB<fEDbins; edB++){
671 for(Int_t c1=0; c1<2; c1++){
672 for(Int_t c2=0; c2<2; c2++){
673 for(Int_t term=0; term<2; term++){
674
675 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2;
676
677 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal;
678 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared;
679 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;
680 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;
681 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;
682 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;
683 //
684 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv;
685 if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW;
686 }// term_2
687
688 for(Int_t c3=0; c3<2; c3++){
689 for(Int_t term=0; term<5; term++){
690
691 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;
692 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 693 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;
42802bba 694 if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted;
be9ef9f9 695 //
696 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;
697
698 }// term_3
b71263d0 699
700 for(Int_t c4=0; c4<2; c4++){
6bb6954b 701 for(Int_t term=0; term<13; term++){
b71263d0 702
703 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;
704 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;
705 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;
42802bba 706 if(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted) delete Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted;
b71263d0 707 //
708 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;
709 }// term_4
710
711 }//c4
be9ef9f9 712 }//c3
713 }//c2
714 }//c1
715 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
716 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
717 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD;
718 if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD;
719 }
720 }
721
722 }// ED
723 }// Mbin
724
725
726 for(Int_t i=0; i<12; i++){
727 if(fFSIss[i]) delete fFSIss[i];
728 if(fFSIos[i]) delete fFSIos[i];
729 }
730 for(Int_t i=0; i<3; i++){// Kt iterator
731 for(Int_t j=0; j<10; j++){// Mbin iterator
732 if(fNormWeight[i][j]) delete fNormWeight[i][j];
733 }
734 }
735
736}
737//________________________________________________________________________
738void AliFourPion::ParInit()
739{
740 cout<<"AliFourPion MyInit() call"<<endl;
741 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;
742
743 fRandomNumber = new TRandom3();
744 fRandomNumber->SetSeed(0);
745
746 //
747 fEventCounter=0;
748 fEventsToMix=3;
749 fZvertexBins=2;//2
750
751 fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
752 fTOFboundry = 2.1;// TOF pid used below this momentum
753
754 ////////////////////////////////////////////////
755 // PadRow Pair Cuts
756 fShareQuality = .5;// max
757 fShareFraction = .05;// max
758 ////////////////////////////////////////////////
759
760
761 fMultLimits[0]=0, fMultLimits[1]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
762 fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
763
764
765
766 if(fPbPbcase) {// PbPb
767 fMultLimit=kMultLimitPbPb;
768 fMbins=fCentBins;
769 fQcut=0.1;
b71263d0 770 fNormQcutLow = 0.15;// 0.15
771 fNormQcutHigh = 0.2;// 0.175
be9ef9f9 772 }else {// pp
773 fMultLimit=kMultLimitpp;
774 fMbins=kMultBinspp;
775 fQcut=0.6;
776 fNormQcutLow = 1.0;
777 fNormQcutHigh = 1.5;
778 }
779
780 fQLowerCut = 0.005;// was 0.005
781 fKupperBound = 1.0;
782 //
783 fKstepY[0] = 1.6;
784 fKmeanY[0] = 0;// central y
785 fKmiddleY[0] = 0;
786
787 // 4x1 (Kt: 0-0.25, 0.25-0.35, 0.35-0.45, 0.45-1.0)
788 if(fKbinsT==4){
789 fKstepT[0] = 0.25; fKstepT[1] = 0.1; fKstepT[2] = 0.1; fKstepT[3] = 0.55;
790 fKmeanT[0] = 0.212; fKmeanT[1] = 0.299; fKmeanT[2] = 0.398; fKmeanT[3] = 0.576;
791 fKmiddleT[0] = 0.125; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.4; fKmiddleT[3] = 0.725;
792 }
793 // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
794 if(fKbinsT==3){
795 fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
796 fKmeanT[0] = 0.240; fKmeanT[1] = 0.369; fKmeanT[2] = 0.576;
797 fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
798 }
799 // 2x1 (Kt: 0-0.35, 0.35-1.0)
800 if(fKbinsT==2){
801 fKstepT[0] = 0.35; fKstepT[1] = 0.65;
802 fKmeanT[0] = 0.264; fKmeanT[1] = 0.500;
803 fKmiddleT[0] = 0.175; fKmiddleT[1] = 0.675;
804 }
805
806 //
807 fQupperBoundWeights = 0.2;
808 fQupperBoundQ2 = 2.0;
5591748e 809 fQupperBoundQ3 = 0.6;
810 fQupperBoundQ4 = 0.6;
be9ef9f9 811 fQbinsQ2 = fQupperBoundQ2/0.005;
812 fQbinsQ3 = fQupperBoundQ3/0.005;
813 fQbinsQ4 = fQupperBoundQ4/0.005;
814 fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
815 for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
816 //
817 fDampStart = 0.5;// was 0.3, then 0.5
818 fDampStep = 0.02;
819
820 //
be9ef9f9 821
822 fEC = new AliFourPionEventCollection **[fZvertexBins];
823 for(UShort_t i=0; i<fZvertexBins; i++){
824
825 fEC[i] = new AliFourPionEventCollection *[fMbins];
826
827 for(UShort_t j=0; j<fMbins; j++){
828
829 fEC[i][j] = new AliFourPionEventCollection(fEventsToMix+1, fMultLimit, kMCarrayLimit, fMCcase);
830 }
831 }
832
b71263d0 833 for(Int_t i=0; i<kMultLimitPbPb; i++) fDefaultsCharSwitch[i]='0';
834 for(Int_t i=0; i<kMultLimitPbPb; i++) {
835 fLowQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
836 fLowQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
837 fLowQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
838 fLowQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 839 fLowQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 840 fLowQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
841 fLowQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
842 fLowQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
843 //
844 fNormQPairSwitch_E0E0[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
845 fNormQPairSwitch_E0E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
846 fNormQPairSwitch_E0E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
847 fNormQPairSwitch_E0E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 848 fNormQPairSwitch_E1E1[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 849 fNormQPairSwitch_E1E2[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
850 fNormQPairSwitch_E1E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
851 fNormQPairSwitch_E2E3[i] = new TArrayC(kMultLimitPbPb,fDefaultsCharSwitch);
852 }
be9ef9f9 853
854 fTempStruct = new AliFourPionTrackStruct[fMultLimit];
b71263d0 855
856
be9ef9f9 857 fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
858
859
860
861 // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
862 if(!fLEGO) {
863 SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
864 if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
865 if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
5591748e 866 if(!fMCcase && !fTabulatePairs) SetMuonCorrections(fLEGO);// Read Muon corrections
be9ef9f9 867 }
868
869 /////////////////////////////////////////////
870 /////////////////////////////////////////////
871
872}
873//________________________________________________________________________
874void AliFourPion::UserCreateOutputObjects()
875{
876 // Create histograms
877 // Called once
878
879 ParInit();// Initialize my settings
880
881
882 fOutputList = new TList();
883 fOutputList->SetOwner();
884
885 TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
886 fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
887 fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
888 fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
889 fOutputList->Add(fVertexDist);
890
891
892 TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
893 fOutputList->Add(fDCAxyDistPlus);
894 TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
895 fOutputList->Add(fDCAzDistPlus);
896 TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
897 fOutputList->Add(fDCAxyDistMinus);
898 TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
899 fOutputList->Add(fDCAzDistMinus);
900
901
902 TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
903 fOutputList->Add(fEvents1);
904 TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
905 fOutputList->Add(fEvents2);
906
907 TH1F *fMultDist0 = new TH1F("fMultDist0","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
908 fMultDist0->GetXaxis()->SetTitle("Multiplicity");
909 fOutputList->Add(fMultDist0);
910
911 TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
912 fMultDist1->GetXaxis()->SetTitle("Multiplicity");
913 fOutputList->Add(fMultDist1);
914
915 TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
916 fMultDist2->GetXaxis()->SetTitle("Multiplicity");
917 fOutputList->Add(fMultDist2);
918
919 TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
920 fMultDist3->GetXaxis()->SetTitle("Multiplicity");
921 fOutputList->Add(fMultDist3);
922
923 TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
924 fOutputList->Add(fPtEtaDist);
925
926 TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
927 fOutputList->Add(fPhiPtDist);
928
929 TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
930 fOutputList->Add(fTOFResponse);
931 TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
932 fOutputList->Add(fTPCResponse);
933
fac2d84b 934 TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",400,0,2);
be9ef9f9 935 fOutputList->Add(fRejectedPairs);
fac2d84b 936 TH1F *fRejectedPairsWeighting = new TH1F("fAcceptedPairsWeighting","",400,0,2);
937 fOutputList->Add(fRejectedPairsWeighting);
938 TH1F *fTotalPairsWeighting = new TH1F("fTotalPairsWeighting","",400,0,2);
939 fOutputList->Add(fTotalPairsWeighting);
940 //
941 TH1F *fRejectedPairsMC = new TH1F("fRejectedPairsMC","",400,0,2);
942 fOutputList->Add(fRejectedPairsMC);
943 TH1F *fRejectedPairsWeightingMC = new TH1F("fAcceptedPairsWeightingMC","",400,0,2);
944 fOutputList->Add(fRejectedPairsWeightingMC);
945 TH1F *fTotalPairsWeightingMC = new TH1F("fTotalPairsWeightingMC","",400,0,2);
946 fOutputList->Add(fTotalPairsWeightingMC);
947
be9ef9f9 948 TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
949 fOutputList->Add(fRejectedEvents);
fac2d84b 950
be9ef9f9 951 TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
952 if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
953 TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
954 if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
955 TH3F *fPairsShareFracDPhiNum = new TH3F("fPairsShareFracDPhiNum","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
956 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiNum);
957 TH3F *fPairsShareFracDPhiDen = new TH3F("fPairsShareFracDPhiDen","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
958 if(fMCcase) fOutputList->Add(fPairsShareFracDPhiDen);
959 TH3D* fPairsPadRowNum = new TH3D("fPairsPadRowNum","", 20,0,1, 159,0,1, 40,0,0.2);
960 if(fMCcase) fOutputList->Add(fPairsPadRowNum);
961 TH3D* fPairsPadRowDen = new TH3D("fPairsPadRowDen","", 20,0,1, 159,0,1, 40,0,0.2);
962 if(fMCcase) fOutputList->Add(fPairsPadRowDen);
963
964
965
966 TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
967 if(fMCcase) fOutputList->Add(fResonanceOSPairs);
968 TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
969 if(fMCcase) fOutputList->Add(fAllOSPairs);
970
971 TH3D *fPrimarySCPionPairs = new TH3D("fPrimarySCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
972 if(fMCcase) fOutputList->Add(fPrimarySCPionPairs);
973 TH3D *fAllSCPionPairs = new TH3D("fAllSCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
974 if(fMCcase) fOutputList->Add(fAllSCPionPairs);
975 TH3D *fPrimaryMCPionPairs = new TH3D("fPrimaryMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
976 if(fMCcase) fOutputList->Add(fPrimaryMCPionPairs);
977 TH3D *fAllMCPionPairs = new TH3D("fAllMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
978 if(fMCcase) fOutputList->Add(fAllMCPionPairs);
979
980 //
981 TH1D *fMuonParents = new TH1D("fMuonParents","",500,0.5,500.5);
982 if(fMCcase) fOutputList->Add(fMuonParents);
983 TH1D *fSecondaryMuonParents = new TH1D("fSecondaryMuonParents","",500,0.5,500.5);
984 if(fMCcase) fOutputList->Add(fSecondaryMuonParents);
985 TH3D *fMuonPionDeltaQinv = new TH3D("fMuonPionDeltaQinv","",2,-0.5,1.5, 20,0,1, 100,-0.2,0.2);
986 if(fMCcase) fOutputList->Add(fMuonPionDeltaQinv);
987 TH1D *fPionCandidates = new TH1D("fPionCandidates","",500,0.5,500.5);
988 if(fMCcase) fOutputList->Add(fPionCandidates);
989 //
990
991
992 TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
993 fOutputList->Add(fAvgMult);
994
995 TH2D *fTrackChi2NDF = new TH2D("fTrackChi2NDF","",20,0,100, 100,0,10);
996 fOutputList->Add(fTrackChi2NDF);
997 TH2D *fTrackTPCncls = new TH2D("fTrackTPCncls","",20,0,100, 110,50,160);
998 fOutputList->Add(fTrackTPCncls);
999
1000
1001 TH1D *fTPNRejects3pion1 = new TH1D("fTPNRejects3pion1","",fQbinsQ3,0,fQupperBoundQ3);
1002 fOutputList->Add(fTPNRejects3pion1);
1003 TH1D *fTPNRejects3pion2 = new TH1D("fTPNRejects3pion2","",fQbinsQ3,0,fQupperBoundQ3);
1004 fOutputList->Add(fTPNRejects3pion2);
1005 TH1D *fTPNRejects4pion1 = new TH1D("fTPNRejects4pion1","",fQbinsQ4,0,fQupperBoundQ4);
1006 fOutputList->Add(fTPNRejects4pion1);
1007
1008 TH3D *fKT3DistTerm1 = new TH3D("fKT3DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
1009 TH3D *fKT3DistTerm5 = new TH3D("fKT3DistTerm5","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
1010 fOutputList->Add(fKT3DistTerm1);
1011 fOutputList->Add(fKT3DistTerm5);
1012 TH3D *fKT4DistTerm1 = new TH3D("fKT4DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
6bb6954b 1013 TH3D *fKT4DistTerm13 = new TH3D("fKT4DistTerm13","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
be9ef9f9 1014 fOutputList->Add(fKT4DistTerm1);
6bb6954b 1015 fOutputList->Add(fKT4DistTerm13);
be9ef9f9 1016
1017
1018 TProfile2D *fKT3AvgpT = new TProfile2D("fKT3AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
1019 fOutputList->Add(fKT3AvgpT);
1020 TProfile2D *fKT4AvgpT = new TProfile2D("fKT4AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
1021 fOutputList->Add(fKT4AvgpT);
9c0557fb 1022 TH3D* fQ3AvgpT = new TH3D("fQ3AvgpT","", 2,-0.5,1.5, fQbinsQ3,0,fQupperBoundQ3, 180,0.1,1.0);
1023 fOutputList->Add(fQ3AvgpT);
1024 TH3D* fQ4AvgpT = new TH3D("fQ4AvgpT","", 2,-0.5,1.5, fQbinsQ4,0,fQupperBoundQ4, 180,0.1,1.0);
1025 fOutputList->Add(fQ4AvgpT);
be9ef9f9 1026
1027
1028 TH1D *fMCWeight3DTerm1SC = new TH1D("fMCWeight3DTerm1SC","", 20,0,0.2);
1029 TH1D *fMCWeight3DTerm1SCden = new TH1D("fMCWeight3DTerm1SCden","", 20,0,0.2);
1030 TH1D *fMCWeight3DTerm2SC = new TH1D("fMCWeight3DTerm2SC","", 20,0,0.2);
1031 TH1D *fMCWeight3DTerm2SCden = new TH1D("fMCWeight3DTerm2SCden","", 20,0,0.2);
1032 TH1D *fMCWeight3DTerm1MC = new TH1D("fMCWeight3DTerm1MC","", 20,0,0.2);
1033 TH1D *fMCWeight3DTerm1MCden = new TH1D("fMCWeight3DTerm1MCden","", 20,0,0.2);
1034 TH1D *fMCWeight3DTerm2MC = new TH1D("fMCWeight3DTerm2MC","", 20,0,0.2);
1035 TH1D *fMCWeight3DTerm2MCden = new TH1D("fMCWeight3DTerm2MCden","", 20,0,0.2);
1036 TH1D *fMCWeight3DTerm3MC = new TH1D("fMCWeight3DTerm3MC","", 20,0,0.2);
1037 TH1D *fMCWeight3DTerm3MCden = new TH1D("fMCWeight3DTerm3MCden","", 20,0,0.2);
1038 TH1D *fMCWeight3DTerm4MC = new TH1D("fMCWeight3DTerm4MC","", 20,0,0.2);
1039 TH1D *fMCWeight3DTerm4MCden = new TH1D("fMCWeight3DTerm4MCden","", 20,0,0.2);
1040 fOutputList->Add(fMCWeight3DTerm1SC);
1041 fOutputList->Add(fMCWeight3DTerm1SCden);
1042 fOutputList->Add(fMCWeight3DTerm2SC);
1043 fOutputList->Add(fMCWeight3DTerm2SCden);
1044 fOutputList->Add(fMCWeight3DTerm1MC);
1045 fOutputList->Add(fMCWeight3DTerm1MCden);
1046 fOutputList->Add(fMCWeight3DTerm2MC);
1047 fOutputList->Add(fMCWeight3DTerm2MCden);
1048 fOutputList->Add(fMCWeight3DTerm3MC);
1049 fOutputList->Add(fMCWeight3DTerm3MCden);
1050 fOutputList->Add(fMCWeight3DTerm4MC);
1051 fOutputList->Add(fMCWeight3DTerm4MCden);
1052
1053
6432c081 1054
1055 for(Int_t mb=0; mb<fMbins; mb++){
1056 if((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit)) continue;
be9ef9f9 1057
6432c081 1058 for(Int_t edB=0; edB<fEDbins; edB++){
1059 for(Int_t c1=0; c1<2; c1++){
1060 for(Int_t c2=0; c2<2; c2++){
1061 for(Int_t term=0; term<2; term++){
1062
1063 TString *nameEx2 = new TString("TwoParticle_Charge1_");
1064 *nameEx2 += c1;
1065 nameEx2->Append("_Charge2_");
1066 *nameEx2 += c2;
1067 nameEx2->Append("_M_");
1068 *nameEx2 += mb;
1069 nameEx2->Append("_ED_");
1070 *nameEx2 += edB;
1071 nameEx2->Append("_Term_");
1072 *nameEx2 += term+1;
1073
1074 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
1075
1076
1077 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1078 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2);
1079 TString *nameEx2QW=new TString(nameEx2->Data());
1080 nameEx2QW->Append("_QW");
1081 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1082 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW);
1083 TString *nameAvgP=new TString(nameEx2->Data());
1084 nameAvgP->Append("_AvgP");
1085 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP = new TProfile2D(nameAvgP->Data(),"",10,0,1, fQbinsQ2,0,fQupperBoundQ2, 0.,1.0,"");
1086 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP);
1087
1088 TString *nameUnitMult=new TString(nameEx2->Data());
1089 nameUnitMult->Append("_UnitMult");
1090 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);
1091 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin);
1092
1093 if(fMCcase){
1094 // Momentum resolution histos
1095 TString *nameIdeal = new TString(nameEx2->Data());
1096 nameIdeal->Append("_Ideal");
1097 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);
1098 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal);
1099 TString *nameSmeared = new TString(nameEx2->Data());
1100 nameSmeared->Append("_Smeared");
1101 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);
1102 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared);
1103 //
1104 // Muon correction histos
1105 TString *nameMuonIdeal=new TString(nameEx2->Data());
1106 nameMuonIdeal->Append("_MuonIdeal");
1107 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1108 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal);
1109 TString *nameMuonSmeared=new TString(nameEx2->Data());
1110 nameMuonSmeared->Append("_MuonSmeared");
1111 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1112 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared);
1113 //
1114 TString *nameMuonPionK2=new TString(nameEx2->Data());
1115 nameMuonPionK2->Append("_MuonPionK2");
1116 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2 = new TH2D(nameMuonPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1117 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2);
1118 //
1119 TString *namePionPionK2=new TString(nameEx2->Data());
1120 namePionPionK2->Append("_PionPionK2");
1121 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2 = new TH2D(namePionPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
1122 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2);
1123 //
1124 //
1125 TString *nameEx2MC=new TString(nameEx2->Data());
1126 nameEx2MC->Append("_MCqinv");
1127 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1128 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
1129 TString *nameEx2MCQW=new TString(nameEx2->Data());
1130 nameEx2MCQW->Append("_MCqinvQW");
1131 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
1132 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
1133 //
1134 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
1135 nameEx2PIDpurityDen->Append("_PIDpurityDen");
1136 Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen = new TH2D(nameEx2PIDpurityDen->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
1137 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
1138 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
1139 nameEx2PIDpurityNum->Append("_PIDpurityNum");
1140 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);
1141 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
1142 }
1143 TString *nameEx2OSLB1 = new TString(nameEx2->Data());
1144 nameEx2OSLB1->Append("_osl_b1");
1145 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);
1146 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL);
1147 nameEx2OSLB1->Append("_QW");
1148 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);
1149 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW);
1150 //
1151 TString *nameEx2OSLB2 = new TString(nameEx2->Data());
1152 nameEx2OSLB2->Append("_osl_b2");
1153 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);
1154 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL);
1155 nameEx2OSLB2->Append("_QW");
1156 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);
1157 fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW);
1158
1159 }// term_2
1160
1161
1162
1163 // skip 3-particle if Tabulate6DPairs is true
1164 if(fTabulatePairs) continue;
1165
1166 for(Int_t c3=0; c3<2; c3++){
1167 for(Int_t term=0; term<5; term++){
be9ef9f9 1168
6432c081 1169 TString *namePC3 = new TString("ThreeParticle_Charge1_");
1170 *namePC3 += c1;
1171 namePC3->Append("_Charge2_");
1172 *namePC3 += c2;
1173 namePC3->Append("_Charge3_");
1174 *namePC3 += c3;
1175 namePC3->Append("_M_");
1176 *namePC3 += mb;
1177 namePC3->Append("_ED_");
1178 *namePC3 += edB;
1179 namePC3->Append("_Term_");
1180 *namePC3 += term+1;
be9ef9f9 1181
6432c081 1182 ///////////////////////////////////////
1183 // skip degenerate histograms
1184 if( (c1+c2+c3)==1) {if(c3!=1) continue;}
1185 if( (c1+c2+c3)==2) {if(c1!=0) continue;}
1186 /////////////////////////////////////////
be9ef9f9 1187
be9ef9f9 1188
6432c081 1189 TString *nameNorm=new TString(namePC3->Data());
1190 nameNorm->Append("_Norm");
1191 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
1192 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3);
1193 //
be9ef9f9 1194
6432c081 1195 TString *name1DQ=new TString(namePC3->Data());
1196 name1DQ->Append("_1D");
1197 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH1D(name1DQ->Data(),"", fQbinsQ3,0,fQupperBoundQ3);
1198 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3);
1199 //
1200 TString *nameKfactor=new TString(namePC3->Data());
1201 nameKfactor->Append("_Kfactor");
1202 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,100, "");
1203 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor);
1204 //
80507acf 1205 TString *nameKfactorW=new TString(namePC3->Data());
1206 nameKfactorW->Append("_KfactorWeighted");
1207 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted = new TProfile(nameKfactorW->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,100, "");
1208 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactorWeighted);
1209 //
6432c081 1210 TString *nameMeanQinv=new TString(namePC3->Data());
1211 nameMeanQinv->Append("_MeanQinv");
1212 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv = new TProfile(nameMeanQinv->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,.2, "");
1213 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv);
1214
1215 if(fMCcase==kTRUE){
1216 // Momentum resolution correction histos
1217 TString *nameMomResIdeal=new TString(namePC3->Data());
1218 nameMomResIdeal->Append("_Ideal");
1219 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);
1220 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal);
1221 TString *nameMomResSmeared=new TString(namePC3->Data());
1222 nameMomResSmeared->Append("_Smeared");
1223 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);
1224 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared);
be9ef9f9 1225 // Muon correction histos
6432c081 1226 TString *nameMuonIdeal=new TString(namePC3->Data());
be9ef9f9 1227 nameMuonIdeal->Append("_MuonIdeal");
6432c081 1228 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);
1229 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal);
1230 TString *nameMuonSmeared=new TString(namePC3->Data());
be9ef9f9 1231 nameMuonSmeared->Append("_MuonSmeared");
6432c081 1232 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);
1233 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared);
be9ef9f9 1234 //
6432c081 1235 TString *nameMuonPionK3=new TString(namePC3->Data());
1236 nameMuonPionK3->Append("_MuonPionK3");
1237 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);
1238 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3);
be9ef9f9 1239 //
6432c081 1240 TString *namePionPionK3=new TString(namePC3->Data());
1241 namePionPionK3->Append("_PionPionK3");
1242 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);
1243 if(mb==0 && edB==0 && term<4) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3);
1244
1245 }// MCcase
1246 //
1247 if(c1==c2 && c1==c3 && term==4 ){
1248 TString *nameTwoPartNorm=new TString(namePC3->Data());
1249 nameTwoPartNorm->Append("_TwoPartNorm");
1250 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);
1251 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm);
be9ef9f9 1252 //
6432c081 1253 TString *nameTwoPartNegNorm=new TString(namePC3->Data());
1254 nameTwoPartNegNorm->Append("_TwoPartNegNorm");
1255 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNegNorm = new TH2D(nameTwoPartNegNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
1256 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNegNorm);
be9ef9f9 1257 //
6432c081 1258 TString *nameTwoPartNormErr=new TString(namePC3->Data());
1259 nameTwoPartNormErr->Append("_TwoPartNormErr");
1260 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNormErr = new TH2D(nameTwoPartNormErr->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
1261 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNormErr);
1262 }// term=4
be9ef9f9 1263
6432c081 1264 }// term_3
be9ef9f9 1265
6432c081 1266 for(Int_t c4=0; c4<2; c4++){
1267 for(Int_t term=0; term<13; term++){
1268
1269 TString *namePC4 = new TString("FourParticle_Charge1_");
1270 *namePC4 += c1;
1271 namePC4->Append("_Charge2_");
1272 *namePC4 += c2;
1273 namePC4->Append("_Charge3_");
1274 *namePC4 += c3;
1275 namePC4->Append("_Charge4_");
1276 *namePC4 += c4;
1277 namePC4->Append("_M_");
1278 *namePC4 += mb;
1279 namePC4->Append("_ED_");
1280 *namePC4 += edB;
1281 namePC4->Append("_Term_");
1282 *namePC4 += term+1;
be9ef9f9 1283
1284 ///////////////////////////////////////
1285 // skip degenerate histograms
6432c081 1286 if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
1287 if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
1288 if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
be9ef9f9 1289 /////////////////////////////////////////
1290
6432c081 1291 TString *nameNorm=new TString(namePC4->Data());
be9ef9f9 1292 nameNorm->Append("_Norm");
6432c081 1293 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);
1294 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4);
be9ef9f9 1295 //
6432c081 1296 TString *name1DQ=new TString(namePC4->Data());
be9ef9f9 1297 name1DQ->Append("_1D");
6432c081 1298 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = new TH1D(name1DQ->Data(),"", fQbinsQ4,0,fQupperBoundQ4);
1299 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4);
be9ef9f9 1300 //
6432c081 1301 TString *nameKfactor=new TString(namePC4->Data());
be9ef9f9 1302 nameKfactor->Append("_Kfactor");
6432c081 1303 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
1304 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor);
be9ef9f9 1305 //
80507acf 1306 TString *nameKfactorW=new TString(namePC4->Data());
1307 nameKfactorW->Append("_KfactorWeighted");
1308 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted = new TProfile(nameKfactorW->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
1309 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactorWeighted);
1310 //
6432c081 1311 if(c1==c2 && c1==c3 && c1==c4 && term==12 ){
1312 TString *nameTwoPartNorm=new TString(namePC4->Data());
1313 nameTwoPartNorm->Append("_TwoPartNorm");
1314 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);
1315 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm);
1316 //
1317 TString *nameTwoPartNegNorm=new TString(namePC4->Data());
1318 nameTwoPartNegNorm->Append("_TwoPartNegNorm");
1319 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNegNorm = new TH2D(nameTwoPartNegNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
1320 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNegNorm);
1321 //
1322 TString *nameTwoPartNormErr=new TString(namePC4->Data());
1323 nameTwoPartNormErr->Append("_TwoPartNormErr");
1324 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNormErr = new TH2D(nameTwoPartNormErr->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
1325 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNormErr);
1326 }
be9ef9f9 1327
1328 if(fMCcase==kTRUE){
1329 // Momentum resolution correction histos
6432c081 1330 TString *nameMomResIdeal=new TString(namePC4->Data());
be9ef9f9 1331 nameMomResIdeal->Append("_Ideal");
6432c081 1332 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);
1333 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal);
1334 TString *nameMomResSmeared=new TString(namePC4->Data());
be9ef9f9 1335 nameMomResSmeared->Append("_Smeared");
6432c081 1336 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);
1337 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared);
be9ef9f9 1338 // Muon correction histos
6432c081 1339 TString *nameMuonIdeal=new TString(namePC4->Data());
be9ef9f9 1340 nameMuonIdeal->Append("_MuonIdeal");
6432c081 1341 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);
1342 if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal);
1343 TString *nameMuonSmeared=new TString(namePC4->Data());
be9ef9f9 1344 nameMuonSmeared->Append("_MuonSmeared");
6432c081 1345 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);
1346 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 1347 //
6432c081 1348 TString *nameMuonPionK4=new TString(namePC4->Data());
1349 nameMuonPionK4->Append("_MuonPionK4");
1350 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);
1351 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 1352 //
6432c081 1353 TString *namePionPionK4=new TString(namePC4->Data());
1354 namePionPionK4->Append("_PionPionK4");
1355 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);
1356 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 1357
1358 }// MCcase
be9ef9f9 1359
6432c081 1360
be9ef9f9 1361 }
6432c081 1362 }
1363
1364 }//c3
1365 }//c2
1366 }//c1
1367 }// ED
1368 }// mbin
be9ef9f9 1369
be9ef9f9 1370
1371
1372 if(fTabulatePairs){
1373
1374 for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
1375 for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
1376 for(Int_t mb=0; mb<fMbins; mb++){
1377 for(Int_t edB=0; edB<fEDbins; edB++){
1378
1379 TString *nameNum = new TString("TPN_num_Kt_");
1380 *nameNum += tKbin;
1381 nameNum->Append("_Ky_");
1382 *nameNum += yKbin;
1383 nameNum->Append("_M_");
1384 *nameNum += mb;
1385 nameNum->Append("_ED_");
1386 *nameNum += edB;
1387
1388 TString *nameDen = new TString("TPN_den_Kt_");
1389 *nameDen += tKbin;
1390 nameDen->Append("_Ky_");
1391 *nameDen += yKbin;
1392 nameDen->Append("_M_");
1393 *nameDen += mb;
1394 nameDen->Append("_ED_");
1395 *nameDen += edB;
1396
1397 if(edB==0){
1398 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = new TH3D(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1399 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD);
1400
1401 KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = new TH3D(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1402 fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD);
1403 }
1404
1405 }
1406 }
1407 }
1408 }
1409
1410 }
1411
1412
1413 TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1414 fOutputList->Add(fQsmearMean);
1415 TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1416 fOutputList->Add(fQsmearSq);
1417 TH2D *fQ2Res = new TH2D("fQ2Res","",20,0,1, 200,-.2,.2);
1418 fOutputList->Add(fQ2Res);
1419 TH2D *fQ3Res = new TH2D("fQ3Res","",20,0,1, 200,-.3,.3);
1420 fOutputList->Add(fQ3Res);
1421 TH2D *fQ4Res = new TH2D("fQ4Res","",20,0,1, 200,-.4,.4);
1422 fOutputList->Add(fQ4Res);
1423
1424 TH2D *DistQinv4pion = new TH2D("DistQinv4pion","",6,0.5,6.5, 20,0,0.1);
1425 fOutputList->Add(DistQinv4pion);
1426 TH2D *DistQinvMC4pion = new TH2D("DistQinvMC4pion","",6,0.5,6.5, 20,0,0.1);
1427 if(fMCcase) fOutputList->Add(DistQinvMC4pion);
1428
1429 TH2D *fAvgQ12VersusQ3 = new TH2D("fAvgQ12VersusQ3","",10,0,0.1, 20,0,0.1);
1430 fOutputList->Add(fAvgQ12VersusQ3);
1431 TH2D *fAvgQ13VersusQ3 = new TH2D("fAvgQ13VersusQ3","",10,0,0.1, 20,0,0.1);
1432 fOutputList->Add(fAvgQ13VersusQ3);
1433 TH2D *fAvgQ23VersusQ3 = new TH2D("fAvgQ23VersusQ3","",10,0,0.1, 20,0,0.1);
1434 fOutputList->Add(fAvgQ23VersusQ3);
1435
5591748e 1436 TH1D *fDistPionParents4 = new TH1D("fDistPionParents4","",4,0.5,4.5);
1437 fOutputList->Add(fDistPionParents4);
1438
6432c081 1439 TH2D *fDistTPCNclsFindable = new TH2D("fDistTPCNclsFindable","", 100,0,0.5, 201,-0.5,200.5);
1440 fDistTPCNclsFindable->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsFindable->GetYaxis()->SetTitle("Ncls Findable");
1441 fOutputList->Add(fDistTPCNclsFindable);
1442 TProfile *fProfileTPCNclsFindable = new TProfile("fProfileTPCNclsFindable","",100,0,0.5, 0,200, "");
1443 fProfileTPCNclsFindable->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsFindable->GetYaxis()->SetTitle("<Ncls Findable>");
1444 fOutputList->Add(fProfileTPCNclsFindable);
1445 //
1446 TH2D *fDistTPCNclsCrossed = new TH2D("fDistTPCNclsCrossed","",100,0,0.5, 201,-0.5,200.5);
1447 fDistTPCNclsCrossed->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsCrossed->GetYaxis()->SetTitle("Ncls Crossed");
1448 fOutputList->Add(fDistTPCNclsCrossed);
1449 TProfile *fProfileTPCNclsCrossed = new TProfile("fProfileTPCNclsCrossed","",100,0,0.5, 0,200, "");
1450 fProfileTPCNclsCrossed->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsCrossed->GetYaxis()->SetTitle("<Ncls Crossed>");
1451 fOutputList->Add(fProfileTPCNclsCrossed);
1452 //
1453 TH2D *fDistTPCNclsFindableRatio = new TH2D("fDistTPCNclsFindableRatio","",100,0,0.5, 100,0,1);
1454 fDistTPCNclsFindableRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsFindableRatio->GetYaxis()->SetTitle("Ncls / Ncls Findable");
1455 fOutputList->Add(fDistTPCNclsFindableRatio);
1456 TProfile *fProfileTPCNclsFindableRatio = new TProfile("fProfileTPCNclsFindableRatio","",100,0,0.5, 0,1, "");
1457 fProfileTPCNclsFindableRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsFindableRatio->GetYaxis()->SetTitle("<Ncls / Ncls Findable>");
1458 fOutputList->Add(fProfileTPCNclsFindableRatio);
1459 //
1460 TH2D *fDistTPCNclsCrossedRatio = new TH2D("fDistTPCNclsCrossedRatio","",100,0,0.5, 100,0,1);
1461 fDistTPCNclsCrossedRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fDistTPCNclsCrossedRatio->GetYaxis()->SetTitle("Ncls / Ncls Crossed");
1462 fOutputList->Add(fDistTPCNclsCrossedRatio);
1463 TProfile *fProfileTPCNclsCrossedRatio = new TProfile("fProfileTPCNclsCrossedRatio","",100,0,0.5, 0,1, "");
1464 fProfileTPCNclsCrossedRatio->GetXaxis()->SetTitle("pT (GeV/c)"); fProfileTPCNclsCrossedRatio->GetYaxis()->SetTitle("<Ncls / Ncls Crossed>");
1465 fOutputList->Add(fProfileTPCNclsCrossedRatio);
1466
42802bba 1467 TH2D *fc4QSFitNum = new TH2D("fc4QSFitNum","",7,0.5,7.5, fQbinsQ4,0,fQupperBoundQ4);
1468 fOutputList->Add(fc4QSFitNum);
1469 TH2D *fc4QSFitDen = new TH2D("fc4QSFitDen","",7,0.5,7.5, fQbinsQ4,0,fQupperBoundQ4);
1470 fOutputList->Add(fc4QSFitDen);
1471
be9ef9f9 1472 ////////////////////////////////////
1473 ///////////////////////////////////
1474
1475 PostData(1, fOutputList);
1476}
1477
1478//________________________________________________________________________
1479void AliFourPion::UserExec(Option_t *)
1480{
1481 // Main loop
1482 // Called for each event
1483 //cout<<"=========== Event # "<<fEventCounter+1<<" ==========="<<endl;
1484 fEventCounter++;
1485
1486 if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1487
1488 fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1489 if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
1490
1491
1492 // Trigger Cut
1493 if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1494 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1495 if(!isSelected1 && !fMCcase) {return;}
1496 }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1497 Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1498 Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1499 if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1500 }else {return;}
1501
1502 ///////////////////////////////////////////////////////////
1503 const AliAODVertex *primaryVertexAOD;
1504 AliCentrality *centrality;// for AODs and ESDs
1505
1506
1507 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1508 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1509 fPIDResponse = inputHandler->GetPIDResponse();
1510
1511
1512 TClonesArray *mcArray = 0x0;
1513 if(fMCcase){
1514 if(fAODcase){
1515 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1516 if(!mcArray || mcArray->GetEntriesFast() >= kMCarrayLimit){
1517 cout<<"No MC particle branch found or Array too large!!"<<endl;
1518 return;
1519 }
1520 }
1521 }
1522
1523
1524 UInt_t status=0;
1525 Int_t positiveTracks=0, negativeTracks=0;
1526 Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1527
1528 Double_t vertex[3]={0};
1529 Int_t zbin=0;
1530 Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1531 /////////////////////////////////////////////////
1532
1533
1534 Float_t centralityPercentile=0;
1535 Float_t cStep=5.0, cStart=0;
1536
1537
1538 if(fAODcase){// AOD case
1539
1540 if(fPbPbcase){
1541 centrality = fAOD->GetCentrality();
1542 centralityPercentile = centrality->GetCentralityPercentile("V0M");
1543 if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1544 if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range. Skipping Event"<<endl;*/ return;}
1545 cout<<"Centrality % = "<<centralityPercentile<<endl;
1546 }
1547
b71263d0 1548
be9ef9f9 1549 ((TH1F*)fOutputList->FindObject("fMultDist0"))->Fill(fAOD->GetNumberOfTracks());
1550
1551 // Pile-up rejection
1552 AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1553 if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1554 else AnaUtil->SetUseMVPlpSelection(kFALSE);
1555 Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD);
1556 if(pileUpCase) return;
1557
1558 ////////////////////////////////
1559 // Vertexing
1560 ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1561 primaryVertexAOD = fAOD->GetPrimaryVertex();
1562 vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1563
1564 if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut
1565 ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1566
1567 if(!fMCcase && primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1568
1569 ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1570
1571 fBfield = fAOD->GetMagneticField();
1572
1573 for(Int_t i=0; i<fZvertexBins; i++){
1574 if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1575 zbin=i;
1576 break;
1577 }
1578 }
1579
1580
1581
1582 /////////////////////////////
1583 // Create Shuffled index list
1584 Int_t randomIndex[fAOD->GetNumberOfTracks()];
1585 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1586 Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1587 /////////////////////////////
1588
1589
1590 // Track loop
1591 for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1592 AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1593 if (!aodtrack) continue;
1594 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1595
1596 status=aodtrack->GetStatus();
1597
1598 if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1599 ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
1600 ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
1601 if(aodtrack->GetTPCNcls() < fMinTPCncls) continue;// TPC nCluster cut
1602 if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
1603
1604 if(fFilterBit != 7){
1605 Bool_t goodTrackOtherFB = kFALSE;
1606 for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1607 AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1608 if(!aodtrack2) continue;
1609 if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1610
1611 if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1612
1613 }
1614 if(!goodTrackOtherFB) continue;
1615 }
1616
1617
86528e76 1618 if(aodtrack->Pt() < 0.16) continue;
be9ef9f9 1619 if(fabs(aodtrack->Eta()) > 0.8) continue;
1620
1621
1622 Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1623 if(!goodMomentum) continue;
1624 aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1625
1626
1627 Double_t dca2[2]={0};
1628 dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1629 dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1630 Double_t dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1631
1632 fTempStruct[myTracks].fStatus = status;
1633 fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1634 fTempStruct[myTracks].fId = aodtrack->GetID();
1635 //
1636 fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1637 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1638 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1639 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1640 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1641 fTempStruct[myTracks].fEta = aodtrack->Eta();
1642 fTempStruct[myTracks].fCharge = aodtrack->Charge();
1643 fTempStruct[myTracks].fDCAXY = dca2[0];
1644 fTempStruct[myTracks].fDCAZ = dca2[1];
1645 fTempStruct[myTracks].fDCA = dca3d;
1646 fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1647 fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1648
1649
1650
86528e76 1651 if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1652
be9ef9f9 1653
1654
1655 // PID section
1656 fTempStruct[myTracks].fElectron = kFALSE;
1657 fTempStruct[myTracks].fPion = kFALSE;
1658 fTempStruct[myTracks].fKaon = kFALSE;
1659 fTempStruct[myTracks].fProton = kFALSE;
1660
1661 Float_t nSigmaTPC[5];
1662 Float_t nSigmaTOF[5];
1663 nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1664 nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1665 fTempStruct[myTracks].fTOFhit = kFALSE;// default
1666 Float_t signalTPC=0, signalTOF=0;
1667 Double_t integratedTimesTOF[10]={0};
1668
1669
1670 Bool_t DoPIDWorkAround=kTRUE;
1671 //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
1672 if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
1673 if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
1674 nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
1675 nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
1676 nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
1677 nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
1678 nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
1679 //
1680 nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
1681 nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
1682 nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
1683 nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
1684 nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
1685 signalTPC = aodtrack->GetTPCsignal();
1686 if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1687 fTempStruct[myTracks].fTOFhit = kTRUE;
1688 signalTOF = aodtrack->GetTOFsignal();
1689 aodtrack->GetIntegratedTimes(integratedTimesTOF);
1690 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1691
1692 }else {// FilterBit 7 PID workaround
1693
1694 for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1695 AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1696 if (!aodTrack2) continue;
1697 if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1698
1699 UInt_t status2=aodTrack2->GetStatus();
1700
1701 nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
1702 nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
1703 nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
1704 nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
1705 nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
1706 //
1707 nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
1708 nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
1709 nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
1710 nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
1711 nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
1712 signalTPC = aodTrack2->GetTPCsignal();
1713
1714 if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1715 fTempStruct[myTracks].fTOFhit = kTRUE;
1716 signalTOF = aodTrack2->GetTOFsignal();
1717 aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1718 }else fTempStruct[myTracks].fTOFhit = kFALSE;
1719
6432c081 1720 //if(aodTrack2->Pt()<0.2) cout<<aodTrack2->GetTPCNclsF()<<" "<<aodTrack2->GetTPCNCrossedRows()<<" "<<aodTrack2->GetTPCNcls()<<" "<<aodTrack2->GetTPCFoundFraction()<<endl;
1721
1722
be9ef9f9 1723 }// aodTrack2
1724 }// FilterBit 7 PID workaround
1725
1726
1727 ///////////////////
1728 ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1729 if(fTempStruct[myTracks].fTOFhit) {
1730 ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1731 }
1732 ///////////////////
1733
1734 // Use TOF if good hit and above threshold
1735 if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1736 if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1737 if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1738 if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1739 if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1740 }else {// TPC info instead
1741 if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1742 if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1743 if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1744 if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1745 }
1746
1747
1748 // Ensure there is only 1 candidate per track
1749 if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1750 if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1751 if(!fTempStruct[myTracks].fPion) continue;// only pions
1752 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1753 if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1754 if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1755 //if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;// superfluous
1756 ////////////////////////
1757 //if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons// superfluous
1758
1759
1760
1761 if(fTempStruct[myTracks].fCharge==+1) {
1762 ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1763 ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1764 }else {
1765 ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1766 ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1767 }
1768
1769 ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1770 ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1771
6432c081 1772 ((TH2D*)fOutputList->FindObject("fDistTPCNclsFindable"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNclsF());
1773 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsFindable"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNclsF());
1774 //
1775 ((TH2D*)fOutputList->FindObject("fDistTPCNclsCrossed"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNCrossedRows());
1776 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsCrossed"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNCrossedRows());
1777 //
1778 if(aodtrack->GetTPCNclsF() > 0){
1779 ((TH2D*)fOutputList->FindObject("fDistTPCNclsFindableRatio"))->Fill(aodtrack->Pt(), double(aodtrack->GetTPCNcls())/double(aodtrack->GetTPCNclsF()));
1780 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsFindableRatio"))->Fill(aodtrack->Pt(), double(aodtrack->GetTPCNcls())/double(aodtrack->GetTPCNclsF()));
1781 }
1782 //
1783 ((TH2D*)fOutputList->FindObject("fDistTPCNclsCrossedRatio"))->Fill(aodtrack->Pt(), aodtrack->GetTPCFoundFraction());
1784 ((TProfile*)fOutputList->FindObject("fProfileTPCNclsCrossedRatio"))->Fill(aodtrack->Pt(), aodtrack->GetTPCFoundFraction());
1785
1786
be9ef9f9 1787 if(fTempStruct[myTracks].fPion) {// pions
1788 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1789 fTempStruct[myTracks].fKey = 1;
1790 }else if(fTempStruct[myTracks].fKaon){// kaons
1791 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1792 fTempStruct[myTracks].fKey = 10;
1793 }else{// protons
1794 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1795 fTempStruct[myTracks].fKey = 100;
1796 }
1797
1798
1799
1800 if(aodtrack->Charge() > 0) positiveTracks++;
1801 else negativeTracks++;
1802
1803 if(fTempStruct[myTracks].fPion) pionCount++;
1804 if(fTempStruct[myTracks].fKaon) kaonCount++;
1805 if(fTempStruct[myTracks].fProton) protonCount++;
1806
1807 myTracks++;
1808
1809 if(fMCcase){// muon mothers
1810 AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
1811 if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
1812 AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
1813 if(parent->IsPhysicalPrimary()){
1814 ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
1815 }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
1816 }
1817 ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
1818 }
1819 }
1820 //cout<<"kinkcount = "<<kinkcount<<" pionkinks = "<<pionkinks<<" primarypionkinks = "<<primarypionkinks<<endl;
1821 }else {// ESD tracks
1822 cout<<"ESDs not supported currently"<<endl;
1823 return;
1824 }
1825
1826 // Generator info only
1827 if(fMCcase && fGeneratorOnly){
1828 myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
1829 for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
1830 if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1831 if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
1832
1833 AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
1834 if(!mcParticle) continue;
1835 if(fabs(mcParticle->Eta())>0.8) continue;
1836 if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
86528e76 1837 if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
be9ef9f9 1838 if(!mcParticle->IsPrimary()) continue;
1839 if(!mcParticle->IsPhysicalPrimary()) continue;
1840 if(abs(mcParticle->GetPdgCode())!=211) continue;
1841
1842 fTempStruct[myTracks].fP[0] = mcParticle->Px();
1843 fTempStruct[myTracks].fP[1] = mcParticle->Py();
1844 fTempStruct[myTracks].fP[2] = mcParticle->Pz();
1845 fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
1846
1847 fTempStruct[myTracks].fId = myTracks;// use my track counter
1848 fTempStruct[myTracks].fLabel = mctrackN;
1849 fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1850 if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1851 fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1852 fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1853 fTempStruct[myTracks].fEta = mcParticle->Eta();
1854 fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
1855 fTempStruct[myTracks].fDCAXY = 0.;
1856 fTempStruct[myTracks].fDCAZ = 0.;
1857 fTempStruct[myTracks].fDCA = 0.;
1858 fTempStruct[myTracks].fPion = kTRUE;
1859 fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2));
1860 fTempStruct[myTracks].fKey = 1;
1861
1862 myTracks++;
1863 pionCount++;
1864 }
1865 }
1866
1867 if(myTracks >= 1) {
1868 ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1869 }
1870
1871
1872 //cout<<"There are "<<myTracks<<" myTracks"<<endl;
1873 //cout<<"pionCount = "<<pionCount<<" kaonCount = "<<kaonCount<<" protonCount = "<<protonCount<<endl;
6432c081 1874 //return;
be9ef9f9 1875
1876 /////////////////////////////////////////
1877 // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1878 if(myTracks < 4) {cout<<"Less than 4 tracks. Skipping Event."<<endl; return;}
1879 /////////////////////////////////////////
1880
1881
1882 ////////////////////////////////
1883 ///////////////////////////////
1884 // Mbin determination
1885 //
1886 // Mbin set to Pion Count Only for pp!!!!!!!
1887 fMbin=-1;
1888 if(!fPbPbcase){
1889 for(Int_t i=0; i<kMultBinspp; i++){
1890 if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1891 // Mbin 0 has 1 pion
1892 }
1893 }else{
1894 for(Int_t i=0; i<fCentBins; i++){
1895 if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1896 fMbin=i;// 0 = most central
1897 break;
1898 }
1899 }
1900 }
b71263d0 1901
be9ef9f9 1902 if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1903
1904 ///////////////////
1905 // can only be called after fMbin has been set
1906 // Radius parameter only matters for Monte-Carlo data
1907 SetFSIindex(fRMax);
1908 ///////////////////
1909
1910 Int_t rBinForTPNMomRes = 10;
1911 if(fMbin==0) {rBinForTPNMomRes=10;}// 10 fm with EW (fRMax should be 11 for normal running)
1912 else if(fMbin==1) {rBinForTPNMomRes=9;}
1913 else if(fMbin<=3) {rBinForTPNMomRes=8;}
1914 else if(fMbin<=5) {rBinForTPNMomRes=7;}
1915 else {rBinForTPNMomRes=6;}
1916
1917 //////////////////////////////////////////////////
1918 fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1919 //////////////////////////////////////////////////
1920
1921
1922
1923 ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1924 ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1925
1926 ////////////////////////////////////
1927 // Add event to buffer if > 0 tracks
1928 if(myTracks > 0){
1929 fEC[zbin][fMbin]->FIFOShift();
1930 (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1931 (fEvt)->fNtracks = myTracks;
1932 (fEvt)->fFillStatus = 1;
1933 for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1934 if(fMCcase){
1935 (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1936 for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1937 AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1938 (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1939 (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1940 (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1941 (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1942 (fEvt)->fMCtracks[i].fPdgCode = tempMCTrack->GetPdgCode();
1943 (fEvt)->fMCtracks[i].fMotherLabel = tempMCTrack->GetMother();
1944 }
1945 }
1946 }
1947
1948
1949
1950 Float_t qinv12=0, qinv13=0, qinv14=0, qinv23=0, qinv24=0, qinv34=0;
1951 Float_t qout=0, qside=0, qlong=0;
1952 Float_t kT12=0;
1953 Float_t q3=0, q3MC=0;
1954 Float_t q4=0, q4MC=0;
1955 Int_t ch1=0, ch2=0, ch3=0, ch4=0;
1956 Int_t bin1=0, bin2=0, bin3=0, bin4=0;
1957 Float_t pVect1[4]={0};
1958 Float_t pVect2[4]={0};
1959 Float_t pVect3[4]={0};
1960 Float_t pVect4[4]={0};
1961 Float_t pVect1MC[4]={0};
1962 Float_t pVect2MC[4]={0};
1963 Float_t pVect3MC[4]={0};
1964 Float_t pVect4MC[4]={0};
1965 Float_t Pparent1[4]={0};
1966 Float_t Pparent2[4]={0};
1967 Float_t Pparent3[4]={0};
1968 Float_t Pparent4[4]={0};
1969 Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0;
1970 Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0;
5591748e 1971 Float_t weight12CC[3]={0};
1972 Float_t weight13CC[3]={0};
1973 Float_t weight14CC[3]={0};
1974 Float_t weight23CC[3]={0};
1975 Float_t weight24CC[3]={0};
1976 Float_t weight34CC[3]={0};
ab890ec3 1977 //Float_t weight12CC_e=0, weight13CC_e=0, weight14CC_e=0, weight23CC_e=0, weight24CC_e=0, weight34CC_e=0;
1978 Float_t weightTotal=0;//, weightTotalErr=0;
be9ef9f9 1979 Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0;
1980 Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0;
1981 Float_t parentQ3=0;
1982 Float_t FSICorr12=0, FSICorr13=0, FSICorr14=0, FSICorr23=0, FSICorr24=0, FSICorr34=0;
1983 Bool_t pionParent1=kFALSE, pionParent2=kFALSE, pionParent3=kFALSE, pionParent4=kFALSE;
1984 Bool_t FilledMCpair12=kFALSE, FilledMCtriplet123=kFALSE;
6432c081 1985 Bool_t Positive1stTripletWeights=kTRUE, Positive2ndTripletWeights=kTRUE;
ab890ec3 1986 Float_t T12=0, T13=0, T14=0, T23=0, T24=0, T34=0;
80507acf 1987 Int_t momBin12=1, momBin13=1, momBin14=1, momBin23=1, momBin24=1, momBin34=1;
1988 Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr14=1.0, MomResCorr23=1.0, MomResCorr24=1.0, MomResCorr34=1.0;
be9ef9f9 1989 //
1990 AliAODMCParticle *mcParticle1=0x0;
1991 AliAODMCParticle *mcParticle2=0x0;
1992
1993
1994 ////////////////////
b71263d0 1995 //Int_t PairCount[7]={0};
1996 //Int_t NormPairCount[7]={0};
be9ef9f9 1997 Int_t KT3index=0, KT4index=0;
1998
b71263d0 1999 // reset to defaults
2000 for(Int_t i=0; i<kMultLimitPbPb; i++) {
2001 fLowQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2002 fLowQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2003 fLowQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2004 fLowQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 2005 fLowQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 2006 fLowQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2007 fLowQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2008 fLowQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2009 //
2010 fNormQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2011 fNormQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2012 fNormQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2013 fNormQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
6bb6954b 2014 fNormQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
b71263d0 2015 fNormQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2016 fNormQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2017 fNormQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2018 }
5591748e 2019
b71263d0 2020
2021 //////////////////////////////////////////
2022 // make low-q pair storage and normalization-pair storage
2023 //
2024 for(Int_t en1=0; en1<=2; en1++){// 1st event number (en1=0 is the same event as current event)
2025 for(Int_t en2=en1; en2<=3; en2++){// 2nd event number (en2=0 is the same event as current event)
6bb6954b 2026 if(en1>1 && en1==en2) continue;
2027
b71263d0 2028 for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {// 1st particle
2029 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2030
5591748e 2031
b71263d0 2032 pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2033 pVect1[1]=(fEvt+en1)->fTracks[i].fP[0]; pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2034 pVect1[2]=(fEvt+en1)->fTracks[i].fP[1]; pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2035 pVect1[3]=(fEvt+en1)->fTracks[i].fP[2]; pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2036 ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2037 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2038
2039 qinv12 = GetQinv(pVect1, pVect2);
2040 kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2041 SetFillBins2(ch1, ch2, bin1, bin2);
2042
80507acf 2043 if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
2044 if(ch1 == ch2 && !fGeneratorOnly){
fac2d84b 2045 Int_t tempChGroup[2]={0,0};
80507acf 2046 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fTotalPairsWeighting"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
b71263d0 2047 if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2048 if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
2049 continue;
be9ef9f9 2050 }
80507acf 2051 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fAcceptedPairsWeighting"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
fac2d84b 2052 }
42802bba 2053 if(fMixedChargeCut && ch1 != ch2 && !fGeneratorOnly && !fMCcase){// remove +- low-q pairs to keep balance between ++ and +- contributions to multi-particle Q3,Q4 projections
fac2d84b 2054 Int_t tempChGroup[2]={0,1};
80507acf 2055 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fTotalPairsWeightingMC"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
fac2d84b 2056 if(!AcceptPairPM((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2057 if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairsMC"))->Fill(qinv12);
2058 continue;
2059 }
80507acf 2060 if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fAcceptedPairsWeightingMC"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
42802bba 2061 }
b71263d0 2062
2063 GetQosl(pVect1, pVect2, qout, qside, qlong);
2064 if( (en1+en2==0)) {
86528e76 2065 if(!fGenerateSignal) Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12);
b71263d0 2066 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2QW->Fill(kT12, qinv12, qinv12);
2067 // osl frame
2068 if((kT12 > 0.2) && (kT12 < 0.3)){
2069 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2070 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
be9ef9f9 2071 }
b71263d0 2072 if((kT12 > 0.6) && (kT12 < 0.7)){
2073 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2074 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
be9ef9f9 2075 }
b71263d0 2076 // unit mult bins
2077 if( (fEvt+en1)->fNtracks%100==0){
2078 Int_t kTindex=0;
2079 if(kT12>0.3) kTindex=1;
2080 Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
2081 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[0].fUnitMultBin->Fill(UnitMultBin, qinv12);
be9ef9f9 2082 }
2083
b71263d0 2084 }
2085 if( (en1+en2==1)) {
86528e76 2086 if(!fGenerateSignal) Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
b71263d0 2087 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2QW->Fill(kT12, qinv12, qinv12);
2088 // osl frame
2089 if((kT12 > 0.2) && (kT12 < 0.3)){
2090 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2091 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2092 }
2093 if((kT12 > 0.6) && (kT12 < 0.7)){
2094 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2095 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
be9ef9f9 2096 }
b71263d0 2097 // unit mult bins
2098 if( (fEvt+en1)->fNtracks%100==0){
2099 Int_t kTindex=0;
2100 if(kT12>0.3) kTindex=1;
2101 Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
2102 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[1].fUnitMultBin->Fill(UnitMultBin, qinv12);
be9ef9f9 2103 }
b71263d0 2104 }
2105 //////////////////////////////////////////
6bb6954b 2106 if(fTabulatePairs && en1==0 && en2<=1 && bin1==bin2){
b71263d0 2107 Float_t kY = 0;
2108 Int_t kTbin=-1, kYbin=-1;
86528e76 2109 Bool_t PairToReject=kFALSE;
2110 if((fEvt+en1)->fTracks[i].fPt < fMinPt || (fEvt+en1)->fTracks[i].fPt > fMaxPt) PairToReject=kTRUE;
2111 if((fEvt+en2)->fTracks[j].fPt < fMinPt || (fEvt+en2)->fTracks[j].fPt > fMaxPt) PairToReject=kTRUE;
2112 if(!PairToReject){
2113 for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(kT12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {kTbin = kIt; break;}}
2114 for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(kY < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {kYbin = kIt; break;}}
2115 if((kTbin<0) || (kYbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2116 if((kTbin>=fKbinsT) || (kYbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2117 if(fGenerateSignal && en2==0) {
2118 Int_t chGroup2[2]={ch1,ch2};
2119 Float_t WInput = MCWeight(chGroup2, fRMax, ffcSqMRC, qinv12, kT12);
b71263d0 2120 KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
86528e76 2121 }else KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2122 }
b71263d0 2123 }
2124
2125 //////////////////////////////////////////////////////////////////////////////
86528e76 2126
b71263d0 2127 if(qinv12 <= fQcut) {
b71263d0 2128 if(en1==0 && en2==0) {fLowQPairSwitch_E0E0[i]->AddAt('1',j);}
2129 if(en1==0 && en2==1) {fLowQPairSwitch_E0E1[i]->AddAt('1',j);}
2130 if(en1==0 && en2==2) {fLowQPairSwitch_E0E2[i]->AddAt('1',j);}
2131 if(en1==0 && en2==3) {fLowQPairSwitch_E0E3[i]->AddAt('1',j);}
6bb6954b 2132 if(en1==1 && en2==1) {fLowQPairSwitch_E1E1[i]->AddAt('1',j);}
b71263d0 2133 if(en1==1 && en2==2) {fLowQPairSwitch_E1E2[i]->AddAt('1',j);}
2134 if(en1==1 && en2==3) {fLowQPairSwitch_E1E3[i]->AddAt('1',j);}
2135 if(en1==2 && en2==3) {fLowQPairSwitch_E2E3[i]->AddAt('1',j);}
2136 }
2137 if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) {
b71263d0 2138 if(en1==0 && en2==0) {fNormQPairSwitch_E0E0[i]->AddAt('1',j);}
2139 if(en1==0 && en2==1) {fNormQPairSwitch_E0E1[i]->AddAt('1',j);}
2140 if(en1==0 && en2==2) {fNormQPairSwitch_E0E2[i]->AddAt('1',j);}
2141 if(en1==0 && en2==3) {fNormQPairSwitch_E0E3[i]->AddAt('1',j);}
6bb6954b 2142 if(en1==1 && en2==1) {fNormQPairSwitch_E1E1[i]->AddAt('1',j);}
b71263d0 2143 if(en1==1 && en2==2) {fNormQPairSwitch_E1E2[i]->AddAt('1',j);}
2144 if(en1==1 && en2==3) {fNormQPairSwitch_E1E3[i]->AddAt('1',j);}
2145 if(en1==2 && en2==3) {fNormQPairSwitch_E2E3[i]->AddAt('1',j);}
be9ef9f9 2146 }
b71263d0 2147
be9ef9f9 2148 }
2149 }
2150 }
b71263d0 2151 }
be9ef9f9 2152
b71263d0 2153 //cout<<PairCount[0]<<" "<<PairCount[1]<<" "<<PairCount[2]<<" "<<PairCount[3]<<" "<<PairCount[4]<<" "<<PairCount[5]<<" "<<PairCount[6]<<endl;
2154 //cout<<NormPairCount[0]<<" "<<NormPairCount[1]<<" "<<NormPairCount[2]<<" "<<NormPairCount[3]<<" "<<NormPairCount[4]<<" "<<NormPairCount[5]<<" "<<NormPairCount[6]<<endl;
2155 ///////////////////////////////////////////////////
2156 // Do not use pairs from events with too many pairs
2157
2158 ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2159
2160 ///////////////////////////////////////////////////
2161
2162
2163 if(fTabulatePairs) return;
be9ef9f9 2164
fac2d84b 2165 /*TF1 *SCpairWeight = new TF1("SCpairWeight","[0] + [1]*x + [2]*exp(-[3]*x)",0,0.2);// same-charge pair weight for monte-carlo data without two-track cuts.
2166 SCpairWeight->FixParameter(0, 0.959);
2167 SCpairWeight->FixParameter(1, 0.278);
2168 SCpairWeight->FixParameter(2, -1.759);
2169 SCpairWeight->FixParameter(3, 115.107);*/
2170
b71263d0 2171 ////////////////////////////////////////////////////
2172 ////////////////////////////////////////////////////
2173 // Normalization counting of 3- and 4-particle terms
2174 for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2175 for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
6bb6954b 2176 if(en2==0 && en3>2) continue;// not needed config
b71263d0 2177 if(en2==1 && en3==en2) continue;// not needed config
2178 for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2179 if(en3==0 && en4>1) continue;// not needed config
6bb6954b 2180 if(en3==1 && en4==3) continue;// not needed configs
2181 if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
b71263d0 2182
2183 for (Int_t i=0; i<myTracks; i++) {// 1st particle
2184 pVect1[1]=(fEvt)->fTracks[i].fP[0];
2185 pVect1[2]=(fEvt)->fTracks[i].fP[1];
2186 pVect1[3]=(fEvt)->fTracks[i].fP[2];
2187 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
be9ef9f9 2188
b71263d0 2189 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
b71263d0 2190 if(en2==0) {if(fNormQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2191 else {if(fNormQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2192
2193 pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2194 pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2195 pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2196 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2197
2198 for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
b71263d0 2199 if(en3==0) {
2200 if(fNormQPairSwitch_E0E0[i]->At(k)=='0') continue;
2201 if(fNormQPairSwitch_E0E0[j]->At(k)=='0') continue;
2202 }else if(en3==1){
2203 if(fNormQPairSwitch_E0E1[i]->At(k)=='0') continue;
2204 if(fNormQPairSwitch_E0E1[j]->At(k)=='0') continue;
2205 }else{
2206 if(fNormQPairSwitch_E0E2[i]->At(k)=='0') continue;
2207 if(fNormQPairSwitch_E1E2[j]->At(k)=='0') continue;
2208 }
be9ef9f9 2209
b71263d0 2210 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2211 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2212 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2213 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2214 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2215 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2216
2217 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2218 if(KT3<=fKT3transition) KT3index=0;
2219 else KT3index=1;
2220
2221 if(en2==0 && en3==0 && en4==0) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fNorm3->Fill(0);
2222 if(en2==1 && en3==2 && en4==3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fNorm3->Fill(0);
2223 if(en2==0 && en3==1 && en4==2) {
2224 if(fill2) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fNorm3->Fill(0);
2225 if(fill3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fNorm3->Fill(0);
2226 if(fill4) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fNorm3->Fill(0);
2227 }
2228
2229
2230 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
b71263d0 2231 if(en4==0){
2232 if(fNormQPairSwitch_E0E0[i]->At(l)=='0') continue;
2233 if(fNormQPairSwitch_E0E0[j]->At(l)=='0') continue;
2234 if(fNormQPairSwitch_E0E0[k]->At(l)=='0') continue;
2235 }else if(en4==1){
6bb6954b 2236 if(en3==0){
2237 if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2238 if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2239 if(fNormQPairSwitch_E0E1[k]->At(l)=='0') continue;
2240 }else{
2241 if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2242 if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2243 if(fNormQPairSwitch_E1E1[k]->At(l)=='0') continue;
2244 }
b71263d0 2245 }else if(en4==2){
2246 if(fNormQPairSwitch_E0E2[i]->At(l)=='0') continue;
2247 if(fNormQPairSwitch_E0E2[j]->At(l)=='0') continue;
2248 if(fNormQPairSwitch_E1E2[k]->At(l)=='0') continue;
be9ef9f9 2249 }else{
b71263d0 2250 if(fNormQPairSwitch_E0E3[i]->At(l)=='0') continue;
2251 if(fNormQPairSwitch_E1E3[j]->At(l)=='0') continue;
2252 if(fNormQPairSwitch_E2E3[k]->At(l)=='0') continue;
be9ef9f9 2253 }
86528e76 2254
b71263d0 2255 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2256 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2257 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2258 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2259 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.;
2260 if(KT4<=fKT4transition) KT4index=0;
2261 else KT4index=1;
be9ef9f9 2262
6bb6954b 2263 Bool_t FillTerms[13]={kFALSE};
b71263d0 2264 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, en2+en3+en4, FillTerms);
2265 //
6bb6954b 2266 for(int ft=0; ft<13; ft++) {
b71263d0 2267 if(FillTerms[ft]) Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fNorm4->Fill(0.);
be9ef9f9 2268 }
be9ef9f9 2269
b71263d0 2270
be9ef9f9 2271 }
2272 }
b71263d0 2273 }
2274 }
2275
be9ef9f9 2276 }
2277 }
b71263d0 2278 }
be9ef9f9 2279
6bb6954b 2280
be9ef9f9 2281
2282
2283 ///////////////////////////////////////////////////////////////////////
2284 ///////////////////////////////////////////////////////////////////////
2285 ///////////////////////////////////////////////////////////////////////
2286 //
2287 //
2288 // Start the Main Correlation Analysis
2289 //
2290 //
2291 ///////////////////////////////////////////////////////////////////////
6bb6954b 2292
be9ef9f9 2293
2294
2295 ////////////////////////////////////////////////////
2296 ////////////////////////////////////////////////////
2297 for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2298 for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
6bb6954b 2299 if(en2==0 && en3>2) continue;// not needed config
be9ef9f9 2300 if(en2==1 && en3==en2) continue;// not needed config
2301 for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2302 if(en3==0 && en4>1) continue;// not needed config
6bb6954b 2303 if(en3==1 && en4==3) continue;// not needed configs
2304 if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2305
be9ef9f9 2306 Int_t ENsum=en2+en3+en4;// 0 or 1 or 3 or 6
6bb6954b 2307
be9ef9f9 2308 /////////////////////////////////////////////////////////////
2309 for (Int_t i=0; i<myTracks; i++) {// 1st particle
2310 pVect1[0]=(fEvt)->fTracks[i].fEaccepted;
2311 pVect1[1]=(fEvt)->fTracks[i].fP[0];
2312 pVect1[2]=(fEvt)->fTracks[i].fP[1];
2313 pVect1[3]=(fEvt)->fTracks[i].fP[2];
2314 ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
86528e76 2315 if((fEvt)->fTracks[i].fPt < fMinPt) continue;
2316 if((fEvt)->fTracks[i].fPt > fMaxPt) continue;
be9ef9f9 2317
2318 /////////////////////////////////////////////////////////////
2319 for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
b71263d0 2320 if(en2==0) {if(fLowQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2321 else {if(fLowQPairSwitch_E0E1[i]->At(j)=='0') continue;}
86528e76 2322 if((fEvt+en2)->fTracks[j].fPt < fMinPt) continue;
2323 if((fEvt+en2)->fTracks[j].fPt > fMaxPt) continue;
2324
be9ef9f9 2325 pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2326 pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2327 pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2328 pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2329 ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2330 qinv12 = GetQinv(pVect1, pVect2);
2331 kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2332 SetFillBins2(ch1, ch2, bin1, bin2);
2333 Int_t kTindex=0;
2334 if(kT12<=0.3) kTindex=0;
2335 else kTindex=1;
2336
2337 FSICorr12 = FSICorrelation(ch1,ch2, qinv12);
2338
2339 // two particle terms filled during tabulation of low-q pairs
2340
2341
2342 if(fMCcase){
2343 FilledMCpair12=kFALSE;
2344
6bb6954b 2345 if(ch1==ch2 && fMbin==0 && qinv12<0.2 && ENsum!=2 && ENsum!=3 && ENsum!=6){
be9ef9f9 2346 for(Int_t rstep=0; rstep<10; rstep++){
2347 Float_t coeff = (rstep)*0.2*(0.18/1.2);
2348 Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2349 if(phi1 > 2*PI) phi1 -= 2*PI;
2350 if(phi1 < 0) phi1 += 2*PI;
2351 Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2352 if(phi2 > 2*PI) phi2 -= 2*PI;
2353 if(phi2 < 0) phi2 += 2*PI;
2354 Float_t deltaphi = phi1 - phi2;
2355 if(deltaphi > PI) deltaphi -= PI;
2356 if(deltaphi < -PI) deltaphi += PI;
2357
2358 if(ENsum==0) ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2359 else ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2360 }
2361
2362 }// pair selection
2363
2364 // Check that label does not exceed stack size
2365 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
2366 if(ENsum==0 && abs((fEvt+en2)->fTracks[j].fLabel) == abs((fEvt)->fTracks[i].fLabel)) continue;
2367 pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2368 pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2369 pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2370 pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2371 pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2372 qinv12MC = GetQinv(pVect1MC, pVect2MC);
86528e76 2373 Int_t chGroup2[2]={ch1,ch2};
2374
2375 if(fGenerateSignal && (ENsum==0 || ENsum==6)){
2376 if(ENsum==0) {
2377 Float_t WInput = MCWeight(chGroup2, fRMax, ffcSqMRC, qinv12MC, 0.);
2378 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12, WInput);
2379 }else{
2380 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
2381 }
2382 }
be9ef9f9 2383
2384 if(qinv12<0.1 && ch1==ch2 && ENsum==0) {
2385 ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC);
2386 ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
2387 ((TH2D*)fOutputList->FindObject("fQ2Res"))->Fill(kT12, qinv12-qinv12MC);
2388 }
2389
2390 // secondary contamination
2391 if(ENsum==0){
2392 mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
2393 mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
2394 if(!mcParticle1 || !mcParticle2) continue;
2395 if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
2396 if(ch1==ch2) {
2397 ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2398 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2399 ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2400 }
2401 }else{
2402 ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2403 if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2404 ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2405 }
2406 }
2407 }
2408 }
2409
2410 if(ENsum==6){// all mixed events
86528e76 2411
be9ef9f9 2412 Float_t rForQW=5.0;
2413 if(fFSIindex<=1) rForQW=10;
2414 else if(fFSIindex==2) rForQW=9;
2415 else if(fFSIindex==3) rForQW=8;
2416 else if(fFSIindex==4) rForQW=7;
2417 else if(fFSIindex==5) rForQW=6;
2418 else if(fFSIindex==6) rForQW=5;
2419 else if(fFSIindex==7) rForQW=4;
2420 else if(fFSIindex==8) rForQW=3;
2421 else rForQW=2;
2422
2423
80507acf 2424 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinv->Fill(qinv12MC, MCWeight(chGroup2, rForQW, ffcSqMRC, qinv12MC, 0.));// was 4,5
2425 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(chGroup2, rForQW, ffcSqMRC, qinv12MC, 0.));// was 4,5
be9ef9f9 2426 // pion purity
2427 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityDen->Fill(kT12, qinv12);
2428 Int_t SCNumber = 1;
2429 Int_t PdgCodeSum = abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode) + abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode);
2430 if(PdgCodeSum==22) SCNumber=1;// e-e
2431 else if(PdgCodeSum==24) SCNumber=2;// e-mu
2432 else if(PdgCodeSum==222) SCNumber=3;// e-pi
2433 else if(PdgCodeSum==332) SCNumber=4;// e-k
2434 else if(PdgCodeSum==2223) SCNumber=5;// e-p
2435 else if(PdgCodeSum==26) SCNumber=6;// mu-mu
2436 else if(PdgCodeSum==224) SCNumber=7;// mu-pi
2437 else if(PdgCodeSum==334) SCNumber=8;// mu-k
2438 else if(PdgCodeSum==2225) SCNumber=9;// mu-p
2439 else if(PdgCodeSum==422) SCNumber=10;// pi-pi
2440 else if(PdgCodeSum==532) SCNumber=11;// pi-k
2441 else if(PdgCodeSum==2423) SCNumber=12;// pi-p
2442 else if(PdgCodeSum==642) SCNumber=13;// k-k
2443 else if(PdgCodeSum==2533) SCNumber=14;// k-p
2444 else if(PdgCodeSum==4424) SCNumber=15;// p-p
2445 else {SCNumber=16;}
2446
2447 Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityNum->Fill(SCNumber, kT12, qinv12);
2448
2449 ///////////////////////
2450 // muon contamination
2451 Pparent1[0]=pVect1MC[0]; Pparent1[1]=pVect1MC[1]; Pparent1[2]=pVect1MC[2]; Pparent1[3]=pVect1MC[3];
2452 Pparent2[0]=pVect2MC[0]; Pparent2[1]=pVect2MC[1]; Pparent2[2]=pVect2MC[2]; Pparent2[3]=pVect2MC[3];
2453 pionParent1=kFALSE; pionParent2=kFALSE;
2454 FilledMCpair12=kTRUE;
2455 //
2456 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
2457 Int_t MotherLabel1 = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fMotherLabel;
2458 if(abs((fEvt)->fMCtracks[MotherLabel1].fPdgCode)==211) {
2459 pionParent1=kTRUE;
2460 Pparent1[1] = (fEvt)->fMCtracks[MotherLabel1].fPx; Pparent1[2] = (fEvt)->fMCtracks[MotherLabel1].fPy; Pparent1[3] = (fEvt)->fMCtracks[MotherLabel1].fPz;
2461 Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
2462 }
2463 //
2464 if(abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13) {
2465 Int_t MotherLabel2 = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fMotherLabel;
2466 if(abs((fEvt+en2)->fMCtracks[MotherLabel2].fPdgCode)==211) {
2467 pionParent2=kTRUE;
2468 Pparent2[1] = (fEvt+en2)->fMCtracks[MotherLabel2].fPx; Pparent2[2] = (fEvt+en2)->fMCtracks[MotherLabel2].fPy; Pparent2[3] = (fEvt+en2)->fMCtracks[MotherLabel2].fPz;
2469 Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
2470 }
2471 }
2472
2473 parentQinv12 = GetQinv(Pparent1, Pparent2);
2474
2475 if(pionParent1 || pionParent2){
2476 if(parentQinv12 > 0.001 && parentQinv12 < 0.3){
2477 Float_t muonPionK12 = FSICorrelation(ch1, ch2, qinv12MC);
2478 Float_t pionPionK12 = FSICorrelation(ch1, ch2, parentQinv12);
2479 for(Int_t term=1; term<=2; term++){
2480 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2481 Float_t Rvalue = 5+Riter;
2482 Float_t WInput = 1.0;
2483 if(term==1) {
2484 WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.);
5591748e 2485 }else{
be9ef9f9 2486 muonPionK12 = 1.0; pionPionK12=1.0;
2487 }
5591748e 2488
be9ef9f9 2489 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput);
2490 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput);
2491 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12);
2492 Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fPionPionK2->Fill(Rvalue, parentQinv12, pionPionK12);
2493 }// Riter
2494 }// term loop
2495
2496 if(ch1==ch2) ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(0., kT12, qinv12MC-parentQinv12);
2497 else ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(1., kT12, qinv12MC-parentQinv12);
2498 }// parentQ check
2499 }// pion parent check
2500 }// muon check
2501
2502
2503 // momentum resolution
2504 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2505 Float_t Rvalue = 5+Riter;
80507acf 2506 Float_t WInput = MCWeight(chGroup2, Rvalue, ffcSqMRC, qinv12MC, 0.);
be9ef9f9 2507 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fIdeal->Fill(Rvalue, qinv12MC, WInput);
2508 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fIdeal->Fill(Rvalue, qinv12MC);
2509 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fSmeared->Fill(Rvalue, qinv12, WInput);
2510 Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fSmeared->Fill(Rvalue, qinv12);
2511 }
2512
2513 }// ENsum check
2514 }// MC array check
2515 }// MC case
2516
2517
6bb6954b 2518
be9ef9f9 2519 /////////////////////////////////////////////////////////////
2520 for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
b71263d0 2521 if(en3==0) {
2522 if(fLowQPairSwitch_E0E0[i]->At(k)=='0') continue;
2523 if(fLowQPairSwitch_E0E0[j]->At(k)=='0') continue;
2524 }else if(en3==1){
2525 if(fLowQPairSwitch_E0E1[i]->At(k)=='0') continue;
2526 if(fLowQPairSwitch_E0E1[j]->At(k)=='0') continue;
2527 }else{
2528 if(fLowQPairSwitch_E0E2[i]->At(k)=='0') continue;
2529 if(fLowQPairSwitch_E1E2[j]->At(k)=='0') continue;
be9ef9f9 2530 }
86528e76 2531 if((fEvt+en3)->fTracks[k].fPt < fMinPt) continue;
2532 if((fEvt+en3)->fTracks[k].fPt > fMaxPt) continue;
2533
be9ef9f9 2534 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
2535 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2536 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2537 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2538 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2539 qinv13 = GetQinv(pVect1, pVect3);
2540 qinv23 = GetQinv(pVect2, pVect3);
2541 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
86528e76 2542 Int_t chGroup3[3]={ch1,ch2,ch3};
2543 Float_t QinvMCGroup3[3]={0};
2544 Float_t kTGroup3[3]={0};
be9ef9f9 2545 FilledMCtriplet123 = kFALSE;
86528e76 2546 if(fMCcase && fGenerateSignal){
2547 if((fEvt+en3)->fTracks[k].fLabel == (fEvt+en2)->fTracks[j].fLabel) continue;
2548 if((fEvt+en3)->fTracks[k].fLabel == (fEvt)->fTracks[i].fLabel) continue;
2549
2550 pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2551 pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2552 pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2553 pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2554 qinv13MC = GetQinv(pVect1MC, pVect3MC);
2555 qinv23MC = GetQinv(pVect2MC, pVect3MC);
2556 QinvMCGroup3[0] = qinv12MC; QinvMCGroup3[1] = qinv13MC; QinvMCGroup3[2] = qinv23MC;
2557 }
2558
be9ef9f9 2559
2560 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2561 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2562
2563 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2564 if(KT3<=fKT3transition) KT3index=0;
2565 else KT3index=1;
2566
2567 FSICorr13 = FSICorrelation(ch1,ch3, qinv13);
2568 FSICorr23 = FSICorrelation(ch2,ch3, qinv23);
80507acf 2569 if(!fGenerateSignal && !fMCcase) {
2570 momBin12 = fMomResC2SC->GetYaxis()->FindBin(qinv12);
2571 momBin13 = fMomResC2SC->GetYaxis()->FindBin(qinv13);
2572 momBin23 = fMomResC2SC->GetYaxis()->FindBin(qinv23);
2573 if(momBin12 >= 20) momBin12 = 19;
2574 if(momBin13 >= 20) momBin13 = 19;
2575 if(momBin23 >= 20) momBin23 = 19;
2576 //
2577 if(ch1==ch2) MomResCorr12 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin12);
2578 else MomResCorr12 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin12);
2579 if(ch1==ch3) MomResCorr13 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin13);
2580 else MomResCorr13 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin13);
2581 if(ch2==ch3) MomResCorr23 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin23);
2582 else MomResCorr23 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin23);
2583 }
be9ef9f9 2584 if(ENsum==0) {
86528e76 2585 Float_t Winput=1.0;
2586 if(fMCcase && fGenerateSignal) Winput = MCWeight3(1, fRMax, ffcSqMRC, chGroup3, QinvMCGroup3, kTGroup3);
2587 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fTerms3->Fill(q3, Winput);
2588 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactor->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23), Winput);
2589 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactorWeighted->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23), MomResCorr12*MomResCorr13*MomResCorr23 * Winput);
be9ef9f9 2590 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv12);
2591 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv13);
2592 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv23);
2593 }
2594 if(ENsum==6) {
2595 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTerms3->Fill(q3);
2596 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv12);
2597 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv13);
2598 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv23);
2599 }
2600 if(ENsum==3){
86528e76 2601 Float_t Winput=1.0;
be9ef9f9 2602 if(fill2) {
86528e76 2603 if(fMCcase && fGenerateSignal) Winput = MCWeight3(2, fRMax, ffcSqMRC, chGroup3, QinvMCGroup3, kTGroup3);
2604 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fTerms3->Fill(q3, Winput);
2605 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactor->Fill(q3, 1/(FSICorr12), Winput);
2606 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12 * Winput);
be9ef9f9 2607 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv12);
2608 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv13);
2609 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv23);
2610 }if(fill3) {
86528e76 2611 if(fMCcase && fGenerateSignal) Winput = MCWeight3(3, fRMax, ffcSqMRC, chGroup3, QinvMCGroup3, kTGroup3);
2612 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fTerms3->Fill(q3, Winput);
2613 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactor->Fill(q3, 1/(FSICorr12), Winput);
2614 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12 * Winput);
be9ef9f9 2615 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv12);
2616 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv13);
2617 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv23);
2618 }if(fill4) {
86528e76 2619 if(fMCcase && fGenerateSignal) Winput = MCWeight3(4, fRMax, ffcSqMRC, chGroup3, QinvMCGroup3, kTGroup3);
2620 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fTerms3->Fill(q3, Winput);
2621 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactor->Fill(q3, 1/(FSICorr12), Winput);
2622 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12 * Winput);
be9ef9f9 2623 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv12);
2624 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv13);
2625 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv23);
2626 }
2627 }
2628
2629 // r3 denominator
2630 if(ENsum==6 && ch1==ch2 && ch1==ch3){
6432c081 2631 Positive1stTripletWeights = kTRUE;
be9ef9f9 2632 //
2633 GetWeight(pVect1, pVect2, weight12, weight12Err);
2634 GetWeight(pVect1, pVect3, weight13, weight13Err);
2635 GetWeight(pVect2, pVect3, weight23, weight23Err);
2636
2637 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {// weight should never be larger than 1
2638 if(fMbin==0 && bin1==0) {
2639 ((TH1D*)fOutputList->FindObject("fTPNRejects3pion1"))->Fill(q3, sqrt(fabs(weight12*weight13*weight23)));
2640 }
2641 }else{
2642
5591748e 2643 Float_t MuonCorr12=1.0, MuonCorr13=1.0, MuonCorr23=1.0;
be9ef9f9 2644 if(!fGenerateSignal && !fMCcase) {
5591748e 2645 MuonCorr12 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin12);
2646 MuonCorr13 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin13);
2647 MuonCorr23 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin23);
be9ef9f9 2648 }
fac2d84b 2649
5591748e 2650 // no MRC, no Muon Correction
2651 weight12CC[0] = ((weight12+1) - ffcSq*FSICorr12 - (1-ffcSq));
2652 weight12CC[0] /= FSICorr12*ffcSq;
2653 weight13CC[0] = ((weight13+1) - ffcSq*FSICorr13 - (1-ffcSq));
2654 weight13CC[0] /= FSICorr13*ffcSq;
2655 weight23CC[0] = ((weight23+1) - ffcSq*FSICorr23 - (1-ffcSq));
2656 weight23CC[0] /= FSICorr23*ffcSq;
2657 if(weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2658 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, sqrt(weight12CC[0]*weight13CC[0]*weight23CC[0]));
2659 }
2660 // no Muon Correction
2661 weight12CC[1] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2662 weight12CC[1] /= FSICorr12*ffcSq;
2663 weight13CC[1] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2664 weight13CC[1] /= FSICorr13*ffcSq;
2665 weight23CC[1] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2666 weight23CC[1] /= FSICorr23*ffcSq;
2667 if(weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2668 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(2, q3, sqrt(weight12CC[1]*weight13CC[1]*weight23CC[1]));
2669 }
2670 // both Corrections
2671 weight12CC[2] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2672 weight12CC[2] /= FSICorr12*ffcSq;
2673 weight12CC[2] *= MuonCorr12;
2674 weight13CC[2] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2675 weight13CC[2] /= FSICorr13*ffcSq;
2676 weight13CC[2] *= MuonCorr13;
2677 weight23CC[2] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2678 weight23CC[2] /= FSICorr23*ffcSq;
2679 weight23CC[2] *= MuonCorr23;
6432c081 2680
5591748e 2681 if(weight12CC[2] < 0 || weight13CC[2] < 0 || weight23CC[2] < 0) {// C2^QS can never be less than unity
be9ef9f9 2682 if(fMbin==0 && bin1==0) {
5591748e 2683 ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC[2]*weight13CC[2]*weight23CC[2])));
be9ef9f9 2684 }
6432c081 2685 if(weight12CC[2] < 0) weight12CC[2]=0;
2686 if(weight13CC[2] < 0) weight13CC[2]=0;
2687 if(weight23CC[2] < 0) weight23CC[2]=0;
2688 Positive1stTripletWeights = kFALSE;
2689 }
2690 /////////////////////////////////////////////////////
2691 weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]);
2692 /////////////////////////////////////////////////////
2693 if(Positive1stTripletWeights){
5591748e 2694 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(3, q3, weightTotal);
ab890ec3 2695 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(4, q3, 1);
6432c081 2696 }else{
2697 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNegNorm->Fill(4, q3, 1);
2698 }
2699 //
2700 // Full Weight reconstruction
2701
2702 for(Int_t RcohIndex=0; RcohIndex<2; RcohIndex++){// Rcoh=0, then Rcoh=Rch
9c0557fb 2703 for(Int_t GIndex=0; GIndex<50; GIndex++){
2704 Int_t FillBin = 5 + RcohIndex*50 + GIndex;
6432c081 2705 Float_t G = 0.02*GIndex;
2706 if(RcohIndex==0){
2707 T12 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight12CC[2])) / (2*pow(1-G,2));
2708 T13 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight13CC[2])) / (2*pow(1-G,2));
2709 T23 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight23CC[2])) / (2*pow(1-G,2));
2710 weightTotal = 2*G*(1-G)*(T12 + T13 + T23) + pow(1-G,2)*(T12*T12 + T13*T13 + T23*T23);
2711 weightTotal += 2*G*pow(1-G,2)*(T12*T13 + T12*T23 + T13*T23) + 2*pow(1-G,3)*T12*T13*T23;
2712 }else{
95567836 2713 T12 = sqrt(weight12CC[2] / (1-G*G));
2714 T13 = sqrt(weight13CC[2] / (1-G*G));
2715 T23 = sqrt(weight23CC[2] / (1-G*G));
2716 weightTotal = (1-G*G)*(T12*T12 + T13*T13 + T23*T23);
2717 weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3)) * T12*T13*T23;
6432c081 2718 }
2719 if(Positive1stTripletWeights){
ab890ec3 2720 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(FillBin, q3, weightTotal);
6432c081 2721 }else{
2722 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNegNorm->Fill(FillBin, q3, weightTotal);
ab890ec3 2723 }
2724 }
6432c081 2725 }
2726 //
2727 /*weight12CC_e = weight12Err*MomResCorr12 / FSICorr12 / ffcSq * MuonCorr12;
43046300 2728 weight13CC_e = weight13Err*MomResCorr13 / FSICorr13 / ffcSq * MuonCorr13;
2729 weight23CC_e = weight23Err*MomResCorr23 / FSICorr23 / ffcSq * MuonCorr23;
2730 if(weight12CC[2]*weight13CC[2]*weight23CC[2] > 0){
6432c081 2731 weightTotalErr = pow(2 * sqrt(3) * weight12CC_e*weight13CC[2]*weight23CC[2] / sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]),2);
43046300 2732 }
2733 weightTotalErr += pow(weight12CC_e,2) + pow(weight13CC_e,2) + pow(weight23CC_e,2);
ab890ec3 2734 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNormErr->Fill(4, q3, weightTotalErr);*/
5591748e 2735
be9ef9f9 2736 }// 1st r3 den check
5591748e 2737
be9ef9f9 2738 }// r3 den
2739
2740
2741 if(ch1==ch2 && ch1==ch3 && ENsum==0){
2742 ((TH3D*)fOutputList->FindObject("fKT3DistTerm1"))->Fill(fMbin+1, KT3, q3);
9c0557fb 2743 if(q3<0.1){
be9ef9f9 2744 Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2745 Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2746 Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2747 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt1);
2748 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt2);
2749 ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt3);
9c0557fb 2750 if(fMbin==0){
2751 ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt1);
2752 ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt2);
2753 ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt3);
2754 }
be9ef9f9 2755 }
9c0557fb 2756
be9ef9f9 2757 }
2758 if(ch1==ch2 && ch1==ch3 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT3DistTerm5"))->Fill(fMbin+1, KT3, q3);
2759
2760
2761
2762
2763 if(fMCcase && ENsum==6 && FilledMCpair12){// for momentum resolution and muon correction
2764 if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
2765
2766 pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2767 pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2768 pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2769 pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2770 qinv13MC = GetQinv(pVect1MC, pVect3MC);
2771 qinv23MC = GetQinv(pVect2MC, pVect3MC);
2772
2773 q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
2774 if(q3<0.1 && ch1==ch2 && ch1==ch3) ((TH2D*)fOutputList->FindObject("fQ3Res"))->Fill(KT3, q3-q3MC);
2775
fac2d84b 2776 Float_t TripletWeightTTC=1.0;// same-charge weights to mimic two-track depletion of same-charge pairs
2777 //if(ch1==ch2 && qinv12>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv12);
2778 //if(ch1==ch3 && qinv13>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv13);
2779 //if(ch2==ch3 && qinv23>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv23);
2780
86528e76 2781
be9ef9f9 2782 Pparent3[0]=pVect3MC[0]; Pparent3[1]=pVect3MC[1]; Pparent3[2]=pVect3MC[2]; Pparent3[3]=pVect3MC[3];
2783 pionParent3=kFALSE;
2784
2785 if(abs((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPdgCode)==13){// muon check
2786 Int_t MotherLabel3 = (fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fMotherLabel;
2787 if(abs((fEvt+en3)->fMCtracks[MotherLabel3].fPdgCode)==211) {
2788 pionParent3=kTRUE;
2789 Pparent3[1] = (fEvt+en3)->fMCtracks[MotherLabel3].fPx; Pparent3[2] = (fEvt+en3)->fMCtracks[MotherLabel3].fPy; Pparent3[3] = (fEvt+en3)->fMCtracks[MotherLabel3].fPz;
2790 Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
2791 }
2792 }
2793
2794 parentQinv13 = GetQinv(Pparent1, Pparent3);
2795 parentQinv23 = GetQinv(Pparent2, Pparent3);
2796 parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
2797
2798 if(parentQinv12 > 0.001 && parentQinv13 > 0.001 && parentQinv23 > 0.001 && parentQ3 < 0.5){
2799 FilledMCtriplet123=kTRUE;
2800 if(pionParent1 || pionParent2 || pionParent3) {// want at least one pion-->muon
2801
2802 Float_t parentQinvGroup3[3]={parentQinv12, parentQinv13, parentQinv23};
2803 //Float_t parentkTGroup3[3]={float(sqrt(pow(Pparent1[1]+Pparent2[1],2) + pow(Pparent1[2]+Pparent2[2],2))/2.),
2804 //float(sqrt(pow(Pparent1[1]+Pparent3[1],2) + pow(Pparent1[2]+Pparent3[2],2))/2.),
2805 //float(sqrt(pow(Pparent2[1]+Pparent3[1],2) + pow(Pparent2[2]+Pparent3[2],2))/2.)};
2806 Float_t parentkTGroup3[3]={0};
2807
2808 ((TH2D*)fOutputList->FindObject("fAvgQ12VersusQ3"))->Fill(parentQ3, parentQinv12);
2809 ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13);
2810 ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23);
2811
5591748e 2812 for(Int_t term=1; term<=4; term++){
2813 if(term==1) {}
2814 else if(term==2) {if(!pionParent1 && !pionParent2) continue;}
2815 else if(term==3) {if(!pionParent1 && !pionParent3) continue;}
2816 else {if(!pionParent2 && !pionParent3) continue;}
be9ef9f9 2817 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2818 Float_t Rvalue = 5+Riter;
2819 Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
2820 Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
2821 Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
fac2d84b 2822 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(1, Rvalue, q3MC, WInput*TripletWeightTTC);
2823 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ3, WInput*TripletWeightTTC);
2824 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(1, Rvalue, q3MC, WInputFSI*TripletWeightTTC);
2825 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(1, Rvalue, parentQ3, WInputParentFSI*TripletWeightTTC);
5591748e 2826 //
fac2d84b 2827 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(2, Rvalue, q3MC, TripletWeightTTC);
2828 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ3, TripletWeightTTC);
2829 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(2, Rvalue, q3MC, TripletWeightTTC);
2830 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(2, Rvalue, parentQ3, TripletWeightTTC);
be9ef9f9 2831 }// Riter
2832 }// term loop
2833
2834 }// pion parent check
2835 }// parentQ check (muon correction)
2836
2837 // 3-pion momentum resolution
2838 for(Int_t term=1; term<=5; term++){
2839 for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2840 Float_t Rvalue = 5+Riter;
80507acf 2841 Float_t WInput = MCWeight3(term, Rvalue, ffcSqMRC, chGroup3, QinvMCGroup3, kTGroup3);
fac2d84b 2842 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fIdeal->Fill(Rvalue, q3MC, WInput*TripletWeightTTC);
2843 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fSmeared->Fill(Rvalue, q3, WInput*TripletWeightTTC);
be9ef9f9 2844 }
2845 }
2846
2847 }// 3rd particle label check
2848 }// MCcase and ENsum==6
2849
5591748e 2850
be9ef9f9 2851
2852
5591748e 2853 /////////////////////////////////////////////////////////////
be9ef9f9 2854 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
b71263d0 2855 if(en4==0){
2856 if(fLowQPairSwitch_E0E0[i]->At(l)=='0') continue;
2857 if(fLowQPairSwitch_E0E0[j]->At(l)=='0') continue;
2858 if(fLowQPairSwitch_E0E0[k]->At(l)=='0') continue;
2859 }else if(en4==1){
6bb6954b 2860 if(en3==0){
2861 if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2862 if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2863 if(fLowQPairSwitch_E0E1[k]->At(l)=='0') continue;
2864 }else{
2865 if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2866 if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2867 if(fLowQPairSwitch_E1E1[k]->At(l)=='0') continue;
2868 }
b71263d0 2869 }else if(en4==2){
2870 if(fLowQPairSwitch_E0E2[i]->At(l)=='0') continue;
2871 if(fLowQPairSwitch_E0E2[j]->At(l)=='0') continue;
2872 if(fLowQPairSwitch_E1E2[k]->At(l)=='0') continue;
2873 }else{
2874 if(fLowQPairSwitch_E0E3[i]->At(l)=='0') continue;
2875 if(fLowQPairSwitch_E1E3[j]->At(l)=='0') continue;
2876 if(fLowQPairSwitch_E2E3[k]->At(l)=='0') continue;
be9ef9f9 2877 }
86528e76 2878 if((fEvt+en4)->fTracks[l].fPt < fMinPt) continue;
2879 if((fEvt+en4)->fTracks[l].fPt > fMaxPt) continue;
2880
be9ef9f9 2881 pVect4[0]=(fEvt+en4)->fTracks[l].fEaccepted;
2882 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2883 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2884 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2885 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2886 qinv14 = GetQinv(pVect1, pVect4);
2887 qinv24 = GetQinv(pVect2, pVect4);
2888 qinv34 = GetQinv(pVect3, pVect4);
2889 q4 = sqrt(pow(q3,2) + pow(qinv14,2) + pow(qinv24,2) + pow(qinv34,2));
86528e76 2890 Int_t chGroup4[4]={ch1,ch2,ch3,ch4};
2891 Float_t QinvMCGroup4[6]={0};
2892 Float_t kTGroup4[6]={0};
2893
2894 if(fMCcase && fGenerateSignal){// for momentum resolution and muon correction
2895 if((fEvt+en4)->fTracks[l].fLabel == (fEvt+en3)->fTracks[k].fLabel) continue;
2896 if((fEvt+en4)->fTracks[l].fLabel == (fEvt+en2)->fTracks[j].fLabel) continue;
2897 if((fEvt+en4)->fTracks[l].fLabel == (fEvt)->fTracks[i].fLabel) continue;
2898
2899 pVect4MC[0]=sqrt(pow((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2900 pVect4MC[1]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPx;
2901 pVect4MC[2]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPy;
2902 pVect4MC[3]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPz;
2903 qinv14MC = GetQinv(pVect1MC, pVect4MC);
2904 qinv24MC = GetQinv(pVect2MC, pVect4MC);
2905 qinv34MC = GetQinv(pVect3MC, pVect4MC);
2906
2907 QinvMCGroup4[0] = qinv12MC; QinvMCGroup4[1] = qinv13MC; QinvMCGroup4[2] = qinv14MC;
2908 QinvMCGroup4[3] = qinv23MC; QinvMCGroup4[4] = qinv24MC; QinvMCGroup4[5] = qinv34MC;
2909 //if(q4<0.1 && ENsum==0 && bin1==bin2 && bin1==bin3 && bin1==bin4) cout<<q4<<" "<<fRMax<<" "<<ffcSqMRC<<" "<<chGroup4[0]<<" "<<chGroup4[1]<<" "<<chGroup4[2]<<" "<<chGroup4[3]<<" "<<QinvMCGroup4[0]<<" "<<QinvMCGroup4[1]<<" "<<QinvMCGroup4[2]<<" "<<QinvMCGroup4[3]<<" "<<QinvMCGroup4[4]<<" "<<QinvMCGroup4[5]<<endl;
be9ef9f9 2910
86528e76 2911 }
be9ef9f9 2912 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2913 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(1, qinv12); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(2, qinv13);
2914 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(3, qinv14); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(4, qinv23);
2915 ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(5, qinv24); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(6, qinv34);
2916 }
2917
2918 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.;
2919 if(KT4<=fKT4transition) KT4index=0;
2920 else KT4index=1;
2921
2922 FSICorr14 = FSICorrelation(ch1,ch4, qinv14);
2923 FSICorr24 = FSICorrelation(ch2,ch4, qinv24);
2924 FSICorr34 = FSICorrelation(ch3,ch4, qinv34);
80507acf 2925
2926 if(!fGenerateSignal && !fMCcase) {
2927 momBin14 = fMomResC2SC->GetYaxis()->FindBin(qinv14);
2928 momBin24 = fMomResC2SC->GetYaxis()->FindBin(qinv24);
2929 momBin34 = fMomResC2SC->GetYaxis()->FindBin(qinv34);
2930 if(momBin14 >= 20) momBin14 = 19;
2931 if(momBin24 >= 20) momBin24 = 19;
2932 if(momBin34 >= 20) momBin34 = 19;
2933 //
2934 if(ch1==ch4) MomResCorr14 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin14);
2935 else MomResCorr14 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin14);
2936 if(ch2==ch4) MomResCorr24 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin24);
2937 else MomResCorr24 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin24);
2938 if(ch3==ch4) MomResCorr34 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin34);
2939 else MomResCorr34 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin34);
2940 }
2941
6bb6954b 2942 Bool_t FillTerms[13]={kFALSE};
be9ef9f9 2943 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, ENsum, FillTerms);
2944 //
6bb6954b 2945 for(int ft=0; ft<13; ft++) {
be9ef9f9 2946 Float_t FSIfactor = 1.0;
80507acf 2947 Float_t MomResWeight = 1.0;
86528e76 2948 Float_t WInput = 1.0;
2949 if(fMCcase && fGenerateSignal) WInput = MCWeight4(ft+1, fRMax, ffcSqMRC, chGroup4, QinvMCGroup4, kTGroup4);
80507acf 2950 if(ft==0) {
2951 FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr14 * FSICorr23 * FSICorr24 * FSICorr34);
2952 MomResWeight = MomResCorr12 * MomResCorr13 * MomResCorr14 * MomResCorr23 * MomResCorr24 * MomResCorr34;
2953 }else if(ft<=4) {
2954 FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr23);
2955 MomResWeight = MomResCorr12 * MomResCorr13 * MomResCorr23;
2956 }else if(ft<=10) {
2957 FSIfactor = 1/(FSICorr12);
2958 MomResWeight = MomResCorr12;
2959 }else if(ft==11) {
2960 FSIfactor = 1/(FSICorr12 * FSICorr34);
2961 MomResWeight = MomResCorr12 * MomResCorr34;
2962 }else {FSIfactor = 1.0; MomResWeight = 1.0;}
be9ef9f9 2963 if(FillTerms[ft]) {
86528e76 2964 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fTerms4->Fill(q4, WInput);
2965 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactor->Fill(q4, FSIfactor, WInput);
2966 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactorWeighted->Fill(q4, FSIfactor, MomResWeight*WInput);
be9ef9f9 2967 }
2968 }
2969
2970 /////////////////////////////////////////////////////////////
2971 // r4{2}
6432c081 2972 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2973 Positive2ndTripletWeights=kTRUE;
2974 //
be9ef9f9 2975 GetWeight(pVect1, pVect4, weight14, weight14Err);
2976 GetWeight(pVect2, pVect4, weight24, weight24Err);
2977 GetWeight(pVect3, pVect4, weight34, weight34Err);
2978
5591748e 2979 Float_t MuonCorr14=1.0, MuonCorr24=1.0, MuonCorr34=1.0;
be9ef9f9 2980 if(!fGenerateSignal && !fMCcase) {
5591748e 2981 MuonCorr14 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin14);
2982 MuonCorr24 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin24);
2983 MuonCorr34 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin34);
be9ef9f9 2984 }
5591748e 2985
2986 // no MRC, no Muon Correction
2987 weight14CC[0] = ((weight14+1) - ffcSq*FSICorr14 - (1-ffcSq));
2988 weight14CC[0] /= FSICorr14*ffcSq;
2989 weight24CC[0] = ((weight24+1) - ffcSq*FSICorr24 - (1-ffcSq));
2990 weight24CC[0] /= FSICorr24*ffcSq;
2991 weight34CC[0] = ((weight34+1) - ffcSq*FSICorr34 - (1-ffcSq));
2992 weight34CC[0] /= FSICorr34*ffcSq;
2993 if(weight14CC[0] > 0 && weight24CC[0] > 0 && weight34CC[0] > 0 && weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2994 weightTotal = sqrt(weight12CC[0]*weight13CC[0]*weight24CC[0]*weight34CC[0]);
2995 weightTotal += sqrt(weight12CC[0]*weight14CC[0]*weight23CC[0]*weight34CC[0]);
2996 weightTotal += sqrt(weight13CC[0]*weight14CC[0]*weight23CC[0]*weight24CC[0]);
2997 weightTotal /= 3.;
2998 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal);
2999 }
3000 // no Muon Correction
3001 weight14CC[1] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
3002 weight14CC[1] /= FSICorr14*ffcSq;
3003 weight24CC[1] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
3004 weight24CC[1] /= FSICorr24*ffcSq;
3005 weight34CC[1] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
3006 weight34CC[1] /= FSICorr34*ffcSq;
3007 if(weight14CC[1] > 0 && weight24CC[1] > 0 && weight34CC[1] > 0 && weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
3008 weightTotal = sqrt(weight12CC[1]*weight13CC[1]*weight24CC[1]*weight34CC[1]);
3009 weightTotal += sqrt(weight12CC[1]*weight14CC[1]*weight23CC[1]*weight34CC[1]);
3010 weightTotal += sqrt(weight13CC[1]*weight14CC[1]*weight23CC[1]*weight24CC[1]);
3011 weightTotal /= 3.;
3012 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(2, q4, weightTotal);
3013 }
3014 // both corrections
3015 weight14CC[2] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
3016 weight14CC[2] /= FSICorr14*ffcSq;
3017 weight14CC[2] *= MuonCorr14;
3018 weight24CC[2] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
3019 weight24CC[2] /= FSICorr24*ffcSq;
3020 weight24CC[2] *= MuonCorr24;
3021 weight34CC[2] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
3022 weight34CC[2] /= FSICorr34*ffcSq;
3023 weight34CC[2] *= MuonCorr34;
3024
3025 if(weight14CC[2] < 0 || weight24CC[2] < 0 || weight34CC[2] < 0) {// C2^QS can never be less than unity
3026 if(fMbin==0 && bin1==0 && KT4index==0) {
3027 ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC[2]*weight23CC[2]*weight34CC[2]*weight14CC[2])));
be9ef9f9 3028 }
6432c081 3029 if(weight14CC[2] < 0) weight14CC[2]=0;
3030 if(weight24CC[2] < 0) weight24CC[2]=0;
3031 if(weight34CC[2] < 0) weight34CC[2]=0;
3032 Positive2ndTripletWeights=kFALSE;
3033 }
3034 /////////////////////////////////////////////////////
3035 weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2]);
3036 weightTotal += sqrt(weight12CC[2]*weight14CC[2]*weight23CC[2]*weight34CC[2]);
3037 weightTotal += sqrt(weight13CC[2]*weight14CC[2]*weight23CC[2]*weight24CC[2]);
3038 weightTotal /= 3.;
3039 if(Positive1stTripletWeights && Positive2ndTripletWeights){
5591748e 3040 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(3, q4, weightTotal);
ab890ec3 3041 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(4, q4, 1);
6432c081 3042 }else{
3043 Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNegNorm->Fill(4, q4, 1);
3044 }
3045 // Full Weight reconstruction
3046 for(Int_t RcohIndex=0; RcohIndex<2; RcohIndex++){// Rcoh=0, then Rcoh=Rch
86528e76 3047 for(Int_t GIndex=0; GIndex<50; GIndex++){// 20 is enough
9c0557fb 3048 Int_t FillBin = 5 + RcohIndex*50 + GIndex;
6432c081 3049 Float_t G = 0.02*GIndex;
3050 if(RcohIndex==0){// Rcoh=0
3051 Float_t a = pow(1-G,2);
3052 Float_t b = 2*G*(1-G);
3053 T12 = (-b + sqrt(pow(b,2) + 4*a*weight12CC[2])) / (2*a);
3054 T13 = (-b + sqrt(pow(b,2) + 4*a*weight13CC[2])) / (2*a);
3055 T14 = (-b + sqrt(pow(b,2) + 4*a*weight14CC[2])) / (2*a);
3056 T23 = (-b + sqrt(pow(b,2) + 4*a*weight23CC[2])) / (2*a);
3057 T24 = (-b + sqrt(pow(b,2) + 4*a*weight24CC[2])) / (2*a);
3058 T34 = (-b + sqrt(pow(b,2) + 4*a*weight34CC[2])) / (2*a);
3059 weightTotal = 2*G*(1-G)*(T12 + T13 + T14 + T23 + T24 + T34) + pow(1-G,2)*(T12*T12 + T13*T13 + T14*T14 + T23*T23 + T24*T24 + T34*T34);// 2-pion
3060 weightTotal += 2*G*pow(1-G,3)*(T12*T34*T34 + T12*T12*T34 + T13*T24*T24 + T13*T13*T24 + T14*T23*T23 + T14*T14*T23);// 2-pair
3061 weightTotal += pow(1-G,4)*(pow(T12,2)*pow(T34,2) + pow(T13,2)*pow(T24,2) + pow(T14,2)*pow(T23,2));// 2-pair fully chaotic
3062 weightTotal += 2*G*pow(1-G,2)*(T12*T13 + T12*T23 + T13*T23 + T12*T14 + T12*T24 + T14*T24);// 3-pion
3063 weightTotal += 2*G*pow(1-G,2)*(T13*T14 + T13*T34 + T14*T34 + T23*T24 + T23*T34 + T24*T34);// 3-pion
3064 weightTotal += 2*pow(1-G,3)*(T12*T13*T23 + T12*T14*T24 + T13*T14*T34 + T23*T24*T34);// 3-pion fully chaotic
3065 weightTotal += 2*G*pow(1-G,3)*(T12*T14*T34 + T12*T14*T23 + T12*T23*T34 + T14*T23*T34);// 4-pion
3066 weightTotal += 2*G*pow(1-G,3)*(T12*T13*T34 + T12*T34*T24 + T12*T24*T13 + T13*T24*T34);// 4-pion
3067 weightTotal += 2*G*pow(1-G,3)*(T14*T13*T23 + T14*T13*T24 + T13*T23*T24 + T14*T24*T23);// 4-pion
3068 weightTotal += 2*pow(1-G,4)*(T12*T13*T24*T34 + T12*T14*T23*T34 + T13*T14*T23*T24);// 4-pion fully chaotic
3069 }else{// Rcoh=Rch
95567836 3070 T12 = sqrt(weight12CC[2] / (1-G*G));
3071 T13 = sqrt(weight13CC[2] / (1-G*G));
3072 T14 = sqrt(weight14CC[2] / (1-G*G));
3073 T23 = sqrt(weight23CC[2] / (1-G*G));
3074 T24 = sqrt(weight24CC[2] / (1-G*G));
3075 T34 = sqrt(weight34CC[2] / (1-G*G));
3076 weightTotal = (1-G*G)*(T12*T12 + T13*T13 + T14*T14 + T23*T23 + T24*T24 + T34*T34);// 2-pion
3077 weightTotal += (4*G*pow(1-G,3)+pow(1-G,4))*(pow(T12,2)*pow(T34,2) + pow(T13,2)*pow(T24,2) + pow(T14,2)*pow(T23,2));// 2-pair
3078 weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3))*(T12*T13*T23);// 3-pion