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