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