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