9ae86e49a0ed8d39889e1f7137a1428dadeda952
[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() < fMinPt) 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 > fMaxPt) continue;// upper P bound
1652       //if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1653       //if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1654      
1655      
1656       // PID section
1657       fTempStruct[myTracks].fElectron = kFALSE;
1658       fTempStruct[myTracks].fPion = kFALSE;
1659       fTempStruct[myTracks].fKaon = kFALSE;
1660       fTempStruct[myTracks].fProton = kFALSE;
1661       
1662       Float_t nSigmaTPC[5];
1663       Float_t nSigmaTOF[5];
1664       nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1665       nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1666       fTempStruct[myTracks].fTOFhit = kFALSE;// default
1667       Float_t signalTPC=0, signalTOF=0;
1668       Double_t integratedTimesTOF[10]={0};
1669
1670     
1671       Bool_t DoPIDWorkAround=kTRUE;
1672       //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
1673       if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
1674       if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
1675         nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
1676         nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
1677         nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
1678         nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
1679         nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
1680         //
1681         nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
1682         nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
1683         nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
1684         nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
1685         nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
1686         signalTPC = aodtrack->GetTPCsignal();
1687         if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1688           fTempStruct[myTracks].fTOFhit = kTRUE;
1689           signalTOF = aodtrack->GetTOFsignal();
1690           aodtrack->GetIntegratedTimes(integratedTimesTOF);
1691         }else fTempStruct[myTracks].fTOFhit = kFALSE;
1692         
1693       }else {// FilterBit 7 PID workaround
1694         
1695         for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1696           AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1697           if (!aodTrack2) continue;
1698           if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1699           
1700           UInt_t status2=aodTrack2->GetStatus();
1701           
1702           nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
1703           nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
1704           nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
1705           nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
1706           nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
1707           //
1708           nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
1709           nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
1710           nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
1711           nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
1712           nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
1713           signalTPC = aodTrack2->GetTPCsignal();
1714           
1715           if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1716             fTempStruct[myTracks].fTOFhit = kTRUE;
1717             signalTOF = aodTrack2->GetTOFsignal();
1718             aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1719           }else fTempStruct[myTracks].fTOFhit = kFALSE;
1720           
1721           //if(aodTrack2->Pt()<0.2) cout<<aodTrack2->GetTPCNclsF()<<"  "<<aodTrack2->GetTPCNCrossedRows()<<"  "<<aodTrack2->GetTPCNcls()<<"  "<<aodTrack2->GetTPCFoundFraction()<<endl;
1722
1723
1724         }// aodTrack2
1725       }// FilterBit 7 PID workaround
1726      
1727      
1728       ///////////////////
1729       ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1730       if(fTempStruct[myTracks].fTOFhit) {
1731         ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1732       }
1733       ///////////////////
1734       
1735       // Use TOF if good hit and above threshold
1736       if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1737         if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1738         if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1739         if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1740         if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1741       }else {// TPC info instead
1742         if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1743         if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1744         if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1745         if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1746       }
1747       
1748     
1749       // Ensure there is only 1 candidate per track
1750       if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1751       if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1752       if(!fTempStruct[myTracks].fPion) continue;// only pions
1753       if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1754       if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1755       if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1756       //if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;// superfluous
1757       ////////////////////////
1758       //if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons// superfluous
1759
1760
1761    
1762       if(fTempStruct[myTracks].fCharge==+1) {
1763         ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1764         ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1765       }else {
1766         ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1767         ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1768       }
1769      
1770       ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1771       ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1772       
1773       ((TH2D*)fOutputList->FindObject("fDistTPCNclsFindable"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNclsF());
1774       ((TProfile*)fOutputList->FindObject("fProfileTPCNclsFindable"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNclsF());
1775       //
1776       ((TH2D*)fOutputList->FindObject("fDistTPCNclsCrossed"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNCrossedRows());
1777       ((TProfile*)fOutputList->FindObject("fProfileTPCNclsCrossed"))->Fill(aodtrack->Pt(), aodtrack->GetTPCNCrossedRows());
1778       //
1779       if(aodtrack->GetTPCNclsF() > 0){
1780         ((TH2D*)fOutputList->FindObject("fDistTPCNclsFindableRatio"))->Fill(aodtrack->Pt(), double(aodtrack->GetTPCNcls())/double(aodtrack->GetTPCNclsF()));
1781         ((TProfile*)fOutputList->FindObject("fProfileTPCNclsFindableRatio"))->Fill(aodtrack->Pt(), double(aodtrack->GetTPCNcls())/double(aodtrack->GetTPCNclsF()));
1782       }      
1783       // 
1784       ((TH2D*)fOutputList->FindObject("fDistTPCNclsCrossedRatio"))->Fill(aodtrack->Pt(), aodtrack->GetTPCFoundFraction());
1785       ((TProfile*)fOutputList->FindObject("fProfileTPCNclsCrossedRatio"))->Fill(aodtrack->Pt(), aodtrack->GetTPCFoundFraction());
1786       
1787
1788       if(fTempStruct[myTracks].fPion) {// pions
1789         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2)); 
1790         fTempStruct[myTracks].fKey = 1;
1791       }else if(fTempStruct[myTracks].fKaon){// kaons
1792         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1793         fTempStruct[myTracks].fKey = 10;
1794       }else{// protons
1795         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1796         fTempStruct[myTracks].fKey = 100;
1797       }
1798       
1799            
1800
1801       if(aodtrack->Charge() > 0) positiveTracks++;
1802       else negativeTracks++;
1803       
1804       if(fTempStruct[myTracks].fPion) pionCount++;
1805       if(fTempStruct[myTracks].fKaon) kaonCount++;
1806       if(fTempStruct[myTracks].fProton) protonCount++;
1807
1808       myTracks++;
1809       
1810       if(fMCcase){// muon mothers
1811         AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
1812         if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
1813           AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
1814           if(parent->IsPhysicalPrimary()){
1815             ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
1816           }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
1817         }
1818         ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
1819       }
1820     }
1821     //cout<<"kinkcount = "<<kinkcount<<"   pionkinks = "<<pionkinks<<"   primarypionkinks = "<<primarypionkinks<<endl;
1822   }else {// ESD tracks
1823     cout<<"ESDs not supported currently"<<endl;
1824     return;
1825   }
1826   
1827   // Generator info only
1828   if(fMCcase && fGeneratorOnly){
1829     myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
1830     for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
1831       if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1832       if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
1833       
1834       AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
1835       if(!mcParticle) continue;
1836       if(fabs(mcParticle->Eta())>0.8) continue;
1837       if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
1838       if(mcParticle->Pt() < fMinPt || mcParticle->Pt() > fMaxPt) continue;
1839       if(!mcParticle->IsPrimary()) continue;
1840       if(!mcParticle->IsPhysicalPrimary()) continue;
1841       if(abs(mcParticle->GetPdgCode())!=211) continue;
1842       
1843       fTempStruct[myTracks].fP[0] = mcParticle->Px();
1844       fTempStruct[myTracks].fP[1] = mcParticle->Py();
1845       fTempStruct[myTracks].fP[2] = mcParticle->Pz();
1846       fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
1847       
1848       fTempStruct[myTracks].fId = myTracks;// use my track counter 
1849       fTempStruct[myTracks].fLabel = mctrackN;
1850       fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1851       if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1852       fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1853       fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1854       fTempStruct[myTracks].fEta = mcParticle->Eta();
1855       fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
1856       fTempStruct[myTracks].fDCAXY = 0.;
1857       fTempStruct[myTracks].fDCAZ = 0.;
1858       fTempStruct[myTracks].fDCA = 0.;
1859       fTempStruct[myTracks].fPion = kTRUE;
1860       fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2)); 
1861       fTempStruct[myTracks].fKey = 1;
1862       
1863       myTracks++;
1864       pionCount++;
1865     }
1866   }
1867   
1868   if(myTracks >= 1) {
1869     ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1870   }
1871  
1872  
1873   //cout<<"There are "<<myTracks<<"  myTracks"<<endl;
1874   //cout<<"pionCount = "<<pionCount<<"   kaonCount = "<<kaonCount<<"   protonCount = "<<protonCount<<endl;
1875   //return;
1876
1877   /////////////////////////////////////////
1878   // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1879   if(myTracks < 4) {cout<<"Less than 4 tracks. Skipping Event."<<endl; return;}
1880   /////////////////////////////////////////
1881  
1882
1883   ////////////////////////////////
1884   ///////////////////////////////
1885   // Mbin determination
1886   //
1887   // Mbin set to Pion Count Only for pp!!!!!!!
1888   fMbin=-1;
1889   if(!fPbPbcase){
1890     for(Int_t i=0; i<kMultBinspp; i++){
1891       if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1892       // Mbin 0 has 1 pion
1893     }
1894   }else{
1895     for(Int_t i=0; i<fCentBins; i++){
1896       if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1897         fMbin=i;// 0 = most central
1898         break;
1899       }
1900     }
1901   }
1902   
1903   if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1904   
1905   ///////////////////
1906   // can only be called after fMbin has been set
1907   // Radius parameter only matters for Monte-Carlo data
1908   SetFSIindex(fRMax);
1909   ///////////////////  
1910   
1911   Int_t rBinForTPNMomRes = 10;
1912   if(fMbin==0) {rBinForTPNMomRes=10;}// 10 fm with EW (fRMax should be 11 for normal running)
1913   else if(fMbin==1) {rBinForTPNMomRes=9;}
1914   else if(fMbin<=3) {rBinForTPNMomRes=8;}
1915   else if(fMbin<=5) {rBinForTPNMomRes=7;}
1916   else {rBinForTPNMomRes=6;}
1917
1918   //////////////////////////////////////////////////
1919   fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1920   //////////////////////////////////////////////////
1921   
1922
1923   
1924   ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1925   ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1926
1927   ////////////////////////////////////
1928   // Add event to buffer if > 0 tracks
1929   if(myTracks > 0){
1930     fEC[zbin][fMbin]->FIFOShift();
1931     (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1932     (fEvt)->fNtracks = myTracks;
1933     (fEvt)->fFillStatus = 1;
1934     for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1935     if(fMCcase){
1936       (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1937       for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1938         AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1939         (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1940         (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1941         (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1942         (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1943         (fEvt)->fMCtracks[i].fPdgCode = tempMCTrack->GetPdgCode();
1944         (fEvt)->fMCtracks[i].fMotherLabel = tempMCTrack->GetMother();
1945       } 
1946     }
1947   }
1948     
1949   
1950   
1951   Float_t qinv12=0, qinv13=0, qinv14=0, qinv23=0, qinv24=0, qinv34=0;
1952   Float_t qout=0, qside=0, qlong=0;
1953   Float_t kT12=0;
1954   Float_t q3=0, q3MC=0;
1955   Float_t q4=0, q4MC=0;
1956   Int_t ch1=0, ch2=0, ch3=0, ch4=0;
1957   Int_t bin1=0, bin2=0, bin3=0, bin4=0;
1958   Float_t pVect1[4]={0}; 
1959   Float_t pVect2[4]={0};
1960   Float_t pVect3[4]={0};
1961   Float_t pVect4[4]={0};
1962   Float_t pVect1MC[4]={0}; 
1963   Float_t pVect2MC[4]={0};
1964   Float_t pVect3MC[4]={0};
1965   Float_t pVect4MC[4]={0};
1966   Float_t Pparent1[4]={0};
1967   Float_t Pparent2[4]={0};
1968   Float_t Pparent3[4]={0};
1969   Float_t Pparent4[4]={0};
1970   Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0;
1971   Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0;
1972   Float_t weight12CC[3]={0};
1973   Float_t weight13CC[3]={0};
1974   Float_t weight14CC[3]={0};
1975   Float_t weight23CC[3]={0};
1976   Float_t weight24CC[3]={0};
1977   Float_t weight34CC[3]={0};
1978   //Float_t weight12CC_e=0, weight13CC_e=0, weight14CC_e=0, weight23CC_e=0, weight24CC_e=0, weight34CC_e=0;
1979   Float_t weightTotal=0;//, weightTotalErr=0;
1980   Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0; 
1981   Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0;
1982   Float_t parentQ3=0;
1983   Float_t FSICorr12=0, FSICorr13=0, FSICorr14=0, FSICorr23=0, FSICorr24=0, FSICorr34=0;
1984   Bool_t pionParent1=kFALSE, pionParent2=kFALSE, pionParent3=kFALSE, pionParent4=kFALSE;
1985   Bool_t FilledMCpair12=kFALSE, FilledMCtriplet123=kFALSE;
1986   Bool_t Positive1stTripletWeights=kTRUE, Positive2ndTripletWeights=kTRUE;
1987   Float_t T12=0, T13=0, T14=0, T23=0, T24=0, T34=0;
1988   Int_t momBin12=1, momBin13=1, momBin14=1, momBin23=1, momBin24=1, momBin34=1;
1989   Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr14=1.0, MomResCorr23=1.0, MomResCorr24=1.0, MomResCorr34=1.0;
1990   //
1991   AliAODMCParticle *mcParticle1=0x0;
1992   AliAODMCParticle *mcParticle2=0x0;
1993   
1994
1995   ////////////////////
1996   //Int_t PairCount[7]={0};
1997   //Int_t NormPairCount[7]={0};
1998   Int_t KT3index=0, KT4index=0;
1999
2000   // reset to defaults
2001   for(Int_t i=0; i<kMultLimitPbPb; i++) {
2002     fLowQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2003     fLowQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2004     fLowQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2005     fLowQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2006     fLowQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2007     fLowQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2008     fLowQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2009     fLowQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2010     //
2011     fNormQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2012     fNormQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2013     fNormQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2014     fNormQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2015     fNormQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2016     fNormQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2017     fNormQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2018     fNormQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
2019   }
2020  
2021   
2022   //////////////////////////////////////////
2023   // make low-q pair storage and normalization-pair storage
2024   // 
2025   for(Int_t en1=0; en1<=2; en1++){// 1st event number (en1=0 is the same event as current event)
2026     for(Int_t en2=en1; en2<=3; en2++){// 2nd event number (en2=0 is the same event as current event)
2027       if(en1>1 && en1==en2) continue;
2028       
2029       for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {// 1st particle
2030         for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2031           
2032           
2033           pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2034           pVect1[1]=(fEvt+en1)->fTracks[i].fP[0];      pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2035           pVect1[2]=(fEvt+en1)->fTracks[i].fP[1];      pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2036           pVect1[3]=(fEvt+en1)->fTracks[i].fP[2];      pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2037           ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
2038           ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2039           
2040           qinv12 = GetQinv(pVect1, pVect2);
2041           kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2042           SetFillBins2(ch1, ch2, bin1, bin2);
2043           
2044           if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
2045           if(ch1 == ch2 && !fGeneratorOnly){
2046             Int_t tempChGroup[2]={0,0};
2047             if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fTotalPairsWeighting"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2048             if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2049               if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
2050               continue;
2051             }
2052             if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fAcceptedPairsWeighting"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2053           }
2054           if(fMixedChargeCut && ch1 != ch2 && !fGeneratorOnly && !fMCcase){// remove +- low-q pairs to keep balance between ++ and +- contributions to multi-particle Q3,Q4 projections
2055             Int_t tempChGroup[2]={0,1};
2056             if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fTotalPairsWeightingMC"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2057             if(!AcceptPairPM((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
2058               if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairsMC"))->Fill(qinv12);
2059               continue;
2060             }
2061             if(en1==0 && en2==1) ((TH1F*)fOutputList->FindObject("fAcceptedPairsWeightingMC"))->Fill(qinv12, MCWeight(tempChGroup, 10, ffcSqMRC, qinv12, 0.));
2062           }
2063           
2064           GetQosl(pVect1, pVect2, qout, qside, qlong);
2065           if( (en1+en2==0)) {
2066             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12);
2067             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2QW->Fill(kT12, qinv12, qinv12);
2068             // osl frame
2069             if((kT12 > 0.2) && (kT12 < 0.3)){
2070               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2071               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2072             }
2073             if((kT12 > 0.6) && (kT12 < 0.7)){  
2074               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2075               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2076             }
2077             // unit mult bins
2078             if( (fEvt+en1)->fNtracks%100==0){
2079               Int_t kTindex=0;
2080               if(kT12>0.3) kTindex=1;
2081               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
2082               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[0].fUnitMultBin->Fill(UnitMultBin, qinv12);
2083             }
2084             
2085           }
2086           if( (en1+en2==1)) {
2087             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
2088             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2QW->Fill(kT12, qinv12, qinv12);
2089             // osl frame
2090             if((kT12 > 0.2) && (kT12 < 0.3)){  
2091               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2092               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2093             }
2094             if((kT12 > 0.6) && (kT12 < 0.7)){  
2095               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
2096               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
2097             }
2098             // unit mult bins
2099             if( (fEvt+en1)->fNtracks%100==0){
2100               Int_t kTindex=0;
2101               if(kT12>0.3) kTindex=1;
2102               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
2103               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[1].fUnitMultBin->Fill(UnitMultBin, qinv12);
2104             }
2105           }
2106           //////////////////////////////////////////
2107           if(fTabulatePairs && en1==0 && en2<=1 && bin1==bin2){
2108             Float_t kY = 0;
2109             Int_t kTbin=-1, kYbin=-1;
2110             
2111             for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(kT12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {kTbin = kIt; break;}} 
2112             for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(kY < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {kYbin = kIt; break;}}
2113             if((kTbin<0) || (kYbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2114             if((kTbin>=fKbinsT) || (kYbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
2115             if(fGenerateSignal && en2==0) {
2116               Int_t chGroup2[2]={ch1,ch2};
2117               Float_t WInput = MCWeight(chGroup2, fRMax, ffcSqMRC, qinv12, kT12);
2118               KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
2119             }else KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
2120           }
2121           
2122           //////////////////////////////////////////////////////////////////////////////
2123           
2124           if(qinv12 <= fQcut) {
2125             if(en1==0 && en2==0) {fLowQPairSwitch_E0E0[i]->AddAt('1',j);}
2126             if(en1==0 && en2==1) {fLowQPairSwitch_E0E1[i]->AddAt('1',j);}
2127             if(en1==0 && en2==2) {fLowQPairSwitch_E0E2[i]->AddAt('1',j);}
2128             if(en1==0 && en2==3) {fLowQPairSwitch_E0E3[i]->AddAt('1',j);}
2129             if(en1==1 && en2==1) {fLowQPairSwitch_E1E1[i]->AddAt('1',j);}
2130             if(en1==1 && en2==2) {fLowQPairSwitch_E1E2[i]->AddAt('1',j);}
2131             if(en1==1 && en2==3) {fLowQPairSwitch_E1E3[i]->AddAt('1',j);}
2132             if(en1==2 && en2==3) {fLowQPairSwitch_E2E3[i]->AddAt('1',j);}
2133           }
2134           if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) {
2135             if(en1==0 && en2==0) {fNormQPairSwitch_E0E0[i]->AddAt('1',j);}
2136             if(en1==0 && en2==1) {fNormQPairSwitch_E0E1[i]->AddAt('1',j);}
2137             if(en1==0 && en2==2) {fNormQPairSwitch_E0E2[i]->AddAt('1',j);}
2138             if(en1==0 && en2==3) {fNormQPairSwitch_E0E3[i]->AddAt('1',j);}
2139             if(en1==1 && en2==1) {fNormQPairSwitch_E1E1[i]->AddAt('1',j);}
2140             if(en1==1 && en2==2) {fNormQPairSwitch_E1E2[i]->AddAt('1',j);}
2141             if(en1==1 && en2==3) {fNormQPairSwitch_E1E3[i]->AddAt('1',j);}
2142             if(en1==2 && en2==3) {fNormQPairSwitch_E2E3[i]->AddAt('1',j);}
2143           }
2144           
2145         }
2146       }
2147     }
2148   }
2149     
2150   //cout<<PairCount[0]<<"  "<<PairCount[1]<<"  "<<PairCount[2]<<"  "<<PairCount[3]<<"  "<<PairCount[4]<<"  "<<PairCount[5]<<"  "<<PairCount[6]<<endl;
2151   //cout<<NormPairCount[0]<<"  "<<NormPairCount[1]<<"  "<<NormPairCount[2]<<"  "<<NormPairCount[3]<<"  "<<NormPairCount[4]<<"  "<<NormPairCount[5]<<"  "<<NormPairCount[6]<<endl;
2152   ///////////////////////////////////////////////////  
2153   // Do not use pairs from events with too many pairs
2154   
2155   ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2156   
2157   ///////////////////////////////////////////////////
2158   
2159   
2160   if(fTabulatePairs) return;
2161
2162   /*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.
2163   SCpairWeight->FixParameter(0, 0.959);
2164   SCpairWeight->FixParameter(1, 0.278);
2165   SCpairWeight->FixParameter(2, -1.759);
2166   SCpairWeight->FixParameter(3, 115.107);*/
2167
2168   ////////////////////////////////////////////////////
2169   ////////////////////////////////////////////////////
2170   // Normalization counting of 3- and 4-particle terms
2171   for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2172     for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2173       if(en2==0 && en3>2) continue;// not needed config
2174       if(en2==1 && en3==en2) continue;// not needed config
2175       for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2176         if(en3==0 && en4>1) continue;// not needed config
2177         if(en3==1 && en4==3) continue;// not needed configs
2178         if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2179         
2180         for (Int_t i=0; i<myTracks; i++) {// 1st particle
2181           pVect1[1]=(fEvt)->fTracks[i].fP[0];
2182           pVect1[2]=(fEvt)->fTracks[i].fP[1];
2183           pVect1[3]=(fEvt)->fTracks[i].fP[2];
2184           ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2185           
2186           for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2187             if(en2==0) {if(fNormQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2188             else {if(fNormQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2189             
2190             pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2191             pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2192             pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2193             ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2194             
2195             for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2196               if(en3==0) {
2197                 if(fNormQPairSwitch_E0E0[i]->At(k)=='0') continue;
2198                 if(fNormQPairSwitch_E0E0[j]->At(k)=='0') continue;
2199               }else if(en3==1){
2200                 if(fNormQPairSwitch_E0E1[i]->At(k)=='0') continue;
2201                 if(fNormQPairSwitch_E0E1[j]->At(k)=='0') continue;
2202               }else{
2203                 if(fNormQPairSwitch_E0E2[i]->At(k)=='0') continue;
2204                 if(fNormQPairSwitch_E1E2[j]->At(k)=='0') continue;
2205               }
2206               
2207               pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2208               pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2209               pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2210               ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2211               Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2212               SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2213               
2214               Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2215               if(KT3<=fKT3transition) KT3index=0;
2216               else KT3index=1;
2217               
2218               if(en2==0 && en3==0 && en4==0) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fNorm3->Fill(0);
2219               if(en2==1 && en3==2 && en4==3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fNorm3->Fill(0);
2220               if(en2==0 && en3==1 && en4==2) {
2221                 if(fill2) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fNorm3->Fill(0);
2222                 if(fill3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fNorm3->Fill(0);
2223                 if(fill4) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fNorm3->Fill(0);
2224               }
2225               
2226               
2227               for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2228                 if(en4==0){
2229                   if(fNormQPairSwitch_E0E0[i]->At(l)=='0') continue;
2230                   if(fNormQPairSwitch_E0E0[j]->At(l)=='0') continue;
2231                   if(fNormQPairSwitch_E0E0[k]->At(l)=='0') continue;
2232                 }else if(en4==1){
2233                   if(en3==0){
2234                     if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2235                     if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2236                     if(fNormQPairSwitch_E0E1[k]->At(l)=='0') continue;
2237                   }else{
2238                     if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2239                     if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2240                     if(fNormQPairSwitch_E1E1[k]->At(l)=='0') continue;
2241                   }
2242                 }else if(en4==2){
2243                   if(fNormQPairSwitch_E0E2[i]->At(l)=='0') continue;
2244                   if(fNormQPairSwitch_E0E2[j]->At(l)=='0') continue;
2245                   if(fNormQPairSwitch_E1E2[k]->At(l)=='0') continue;
2246                 }else{
2247                   if(fNormQPairSwitch_E0E3[i]->At(l)=='0') continue;
2248                   if(fNormQPairSwitch_E1E3[j]->At(l)=='0') continue;
2249                   if(fNormQPairSwitch_E2E3[k]->At(l)=='0') continue;
2250                 }
2251
2252                 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2253                 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2254                 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2255                 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2256                 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.;
2257                 if(KT4<=fKT4transition) KT4index=0;
2258                 else KT4index=1;
2259                 
2260                 Bool_t FillTerms[13]={kFALSE};
2261                 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, en2+en3+en4, FillTerms);
2262                 //
2263                 for(int ft=0; ft<13; ft++) {
2264                   if(FillTerms[ft]) Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fNorm4->Fill(0.); 
2265                 }
2266                 
2267                 
2268               }
2269             }
2270           }
2271         }  
2272         
2273       }
2274     }
2275   }
2276     
2277
2278    
2279
2280     ///////////////////////////////////////////////////////////////////////
2281     ///////////////////////////////////////////////////////////////////////
2282     ///////////////////////////////////////////////////////////////////////
2283     //
2284     //
2285     // Start the Main Correlation Analysis
2286     //
2287     //
2288     ///////////////////////////////////////////////////////////////////////
2289   
2290
2291
2292     ////////////////////////////////////////////////////
2293     ////////////////////////////////////////////////////
2294     for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2295       for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2296         if(en2==0 && en3>2) continue;// not needed config
2297         if(en2==1 && en3==en2) continue;// not needed config
2298         for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2299           if(en3==0 && en4>1) continue;// not needed config
2300           if(en3==1 && en4==3) continue;// not needed configs
2301           if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2302           
2303           Int_t ENsum=en2+en3+en4;// 0 or 1 or 3 or 6
2304           
2305           /////////////////////////////////////////////////////////////
2306           for (Int_t i=0; i<myTracks; i++) {// 1st particle
2307             pVect1[0]=(fEvt)->fTracks[i].fEaccepted;
2308             pVect1[1]=(fEvt)->fTracks[i].fP[0];
2309             pVect1[2]=(fEvt)->fTracks[i].fP[1];
2310             pVect1[3]=(fEvt)->fTracks[i].fP[2];
2311             ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2312
2313             /////////////////////////////////////////////////////////////
2314             for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2315               if(en2==0) {if(fLowQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2316               else {if(fLowQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2317
2318               pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2319               pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2320               pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2321               pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2322               ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2323               qinv12 = GetQinv(pVect1, pVect2);
2324               kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2325               SetFillBins2(ch1, ch2, bin1, bin2);
2326               Int_t kTindex=0;
2327               if(kT12<=0.3) kTindex=0;
2328               else kTindex=1;
2329               
2330               FSICorr12 = FSICorrelation(ch1,ch2, qinv12);
2331               
2332               // two particle terms filled during tabulation of low-q pairs
2333               
2334               
2335               if(fMCcase){
2336                 FilledMCpair12=kFALSE;
2337
2338                 if(ch1==ch2 && fMbin==0 && qinv12<0.2 && ENsum!=2 && ENsum!=3 && ENsum!=6){
2339                   for(Int_t rstep=0; rstep<10; rstep++){
2340                     Float_t coeff = (rstep)*0.2*(0.18/1.2);
2341                     Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2342                     if(phi1 > 2*PI) phi1 -= 2*PI;
2343                     if(phi1 < 0) phi1 += 2*PI;
2344                     Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2345                     if(phi2 > 2*PI) phi2 -= 2*PI;
2346                     if(phi2 < 0) phi2 += 2*PI;
2347                     Float_t deltaphi = phi1 - phi2;
2348                     if(deltaphi > PI) deltaphi -= PI;
2349                     if(deltaphi < -PI) deltaphi += PI;
2350                     
2351                     if(ENsum==0) ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2352                     else ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2353                   }
2354                   
2355                 }// pair selection
2356
2357                 // Check that label does not exceed stack size
2358                 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
2359                   if(ENsum==0 && abs((fEvt+en2)->fTracks[j].fLabel) == abs((fEvt)->fTracks[i].fLabel)) continue;
2360                   pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
2361                   pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2362                   pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2363                   pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2364                   pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2365                   qinv12MC = GetQinv(pVect1MC, pVect2MC);
2366                   
2367                   if(qinv12<0.1 && ch1==ch2 && ENsum==0) {
2368                     ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC); 
2369                     ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
2370                     ((TH2D*)fOutputList->FindObject("fQ2Res"))->Fill(kT12, qinv12-qinv12MC);
2371                   }
2372                                   
2373                   // secondary contamination
2374                   if(ENsum==0){
2375                     mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
2376                     mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
2377                     if(!mcParticle1 || !mcParticle2) continue;
2378                     if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
2379                       if(ch1==ch2) {
2380                         ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2381                         if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2382                           ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2383                         }             
2384                       }else{
2385                         ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2386                         if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2387                           ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2388                         }
2389                       }
2390                     }
2391                   }
2392                   
2393                   if(ENsum==6){// all mixed events
2394                     Int_t chGroup2[2]={ch1,ch2};
2395
2396                     Float_t rForQW=5.0;
2397                     if(fFSIindex<=1) rForQW=10;
2398                     else if(fFSIindex==2) rForQW=9;
2399                     else if(fFSIindex==3) rForQW=8;
2400                     else if(fFSIindex==4) rForQW=7;
2401                     else if(fFSIindex==5) rForQW=6;
2402                     else if(fFSIindex==6) rForQW=5;
2403                     else if(fFSIindex==7) rForQW=4;
2404                     else if(fFSIindex==8) rForQW=3;
2405                     else rForQW=2;
2406                     
2407                     
2408                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinv->Fill(qinv12MC, MCWeight(chGroup2, rForQW, ffcSqMRC, qinv12MC, 0.));// was 4,5
2409                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(chGroup2, rForQW, ffcSqMRC, qinv12MC, 0.));// was 4,5
2410                     // pion purity
2411                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityDen->Fill(kT12, qinv12);
2412                     Int_t SCNumber = 1;
2413                     Int_t PdgCodeSum = abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode) + abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode);
2414                     if(PdgCodeSum==22) SCNumber=1;// e-e
2415                     else if(PdgCodeSum==24) SCNumber=2;// e-mu
2416                     else if(PdgCodeSum==222) SCNumber=3;// e-pi
2417                     else if(PdgCodeSum==332) SCNumber=4;// e-k
2418                     else if(PdgCodeSum==2223) SCNumber=5;// e-p
2419                     else if(PdgCodeSum==26) SCNumber=6;// mu-mu
2420                     else if(PdgCodeSum==224) SCNumber=7;// mu-pi
2421                     else if(PdgCodeSum==334) SCNumber=8;// mu-k
2422                     else if(PdgCodeSum==2225) SCNumber=9;// mu-p
2423                     else if(PdgCodeSum==422) SCNumber=10;// pi-pi
2424                     else if(PdgCodeSum==532) SCNumber=11;// pi-k
2425                     else if(PdgCodeSum==2423) SCNumber=12;// pi-p
2426                     else if(PdgCodeSum==642) SCNumber=13;// k-k
2427                     else if(PdgCodeSum==2533) SCNumber=14;// k-p
2428                     else if(PdgCodeSum==4424) SCNumber=15;// p-p
2429                     else {SCNumber=16;}
2430                     
2431                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityNum->Fill(SCNumber, kT12, qinv12);
2432                     
2433                     ///////////////////////
2434                     // muon contamination
2435                     Pparent1[0]=pVect1MC[0]; Pparent1[1]=pVect1MC[1]; Pparent1[2]=pVect1MC[2]; Pparent1[3]=pVect1MC[3];
2436                     Pparent2[0]=pVect2MC[0]; Pparent2[1]=pVect2MC[1]; Pparent2[2]=pVect2MC[2]; Pparent2[3]=pVect2MC[3];
2437                     pionParent1=kFALSE; pionParent2=kFALSE;
2438                     FilledMCpair12=kTRUE;
2439                     //
2440                     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
2441                       Int_t MotherLabel1 = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fMotherLabel;
2442                       if(abs((fEvt)->fMCtracks[MotherLabel1].fPdgCode)==211) {
2443                         pionParent1=kTRUE;
2444                         Pparent1[1] = (fEvt)->fMCtracks[MotherLabel1].fPx; Pparent1[2] = (fEvt)->fMCtracks[MotherLabel1].fPy; Pparent1[3] = (fEvt)->fMCtracks[MotherLabel1].fPz;
2445                         Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
2446                       }
2447                       // 
2448                       if(abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13) {
2449                         Int_t MotherLabel2 = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fMotherLabel;
2450                         if(abs((fEvt+en2)->fMCtracks[MotherLabel2].fPdgCode)==211) {
2451                           pionParent2=kTRUE;
2452                           Pparent2[1] = (fEvt+en2)->fMCtracks[MotherLabel2].fPx; Pparent2[2] = (fEvt+en2)->fMCtracks[MotherLabel2].fPy; Pparent2[3] = (fEvt+en2)->fMCtracks[MotherLabel2].fPz;
2453                           Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
2454                         }
2455                       }
2456                       
2457                       parentQinv12 = GetQinv(Pparent1, Pparent2);
2458                       
2459                       if(pionParent1 || pionParent2){
2460                         if(parentQinv12 > 0.001 && parentQinv12 < 0.3){
2461                           Float_t muonPionK12 = FSICorrelation(ch1, ch2, qinv12MC);
2462                           Float_t pionPionK12 = FSICorrelation(ch1, ch2, parentQinv12);
2463                           for(Int_t term=1; term<=2; term++){
2464                             for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2465                               Float_t Rvalue = 5+Riter;
2466                               Float_t WInput = 1.0;
2467                               if(term==1) {
2468                                 WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.);
2469                               }else{
2470                                 muonPionK12 = 1.0; pionPionK12=1.0;
2471                               }
2472                               
2473                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput);
2474                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput);
2475                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12);
2476                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fPionPionK2->Fill(Rvalue, parentQinv12, pionPionK12);
2477                             }// Riter
2478                           }// term loop
2479                           
2480                           if(ch1==ch2) ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(0., kT12, qinv12MC-parentQinv12);
2481                           else ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(1., kT12, qinv12MC-parentQinv12);
2482                         }// parentQ check
2483                       }// pion parent check
2484                     }// muon check
2485                   
2486                     
2487                     // momentum resolution
2488                     for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2489                       Float_t Rvalue = 5+Riter;
2490                       Float_t WInput = MCWeight(chGroup2, Rvalue, ffcSqMRC, qinv12MC, 0.);
2491                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fIdeal->Fill(Rvalue, qinv12MC, WInput);
2492                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fIdeal->Fill(Rvalue, qinv12MC);
2493                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fSmeared->Fill(Rvalue, qinv12, WInput);
2494                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fSmeared->Fill(Rvalue, qinv12);
2495                     }
2496                     
2497                   }// ENsum check
2498                 }// MC array check
2499               }// MC case
2500               
2501                 
2502
2503               /////////////////////////////////////////////////////////////
2504               for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2505                 if(en3==0) {
2506                   if(fLowQPairSwitch_E0E0[i]->At(k)=='0') continue;
2507                   if(fLowQPairSwitch_E0E0[j]->At(k)=='0') continue;
2508                 }else if(en3==1){
2509                   if(fLowQPairSwitch_E0E1[i]->At(k)=='0') continue;
2510                   if(fLowQPairSwitch_E0E1[j]->At(k)=='0') continue;
2511                 }else{
2512                   if(fLowQPairSwitch_E0E2[i]->At(k)=='0') continue;
2513                   if(fLowQPairSwitch_E1E2[j]->At(k)=='0') continue;
2514                 }
2515                 
2516                 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
2517                 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2518                 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2519                 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2520                 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2521                 qinv13 = GetQinv(pVect1, pVect3);
2522                 qinv23 = GetQinv(pVect2, pVect3);
2523                 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2524                 
2525                 FilledMCtriplet123 = kFALSE;
2526                 
2527                 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2528                 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2529                 
2530                 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2531                 if(KT3<=fKT3transition) KT3index=0;
2532                 else KT3index=1;
2533                 
2534                 FSICorr13 = FSICorrelation(ch1,ch3, qinv13);
2535                 FSICorr23 = FSICorrelation(ch2,ch3, qinv23);
2536                 if(!fGenerateSignal && !fMCcase) {
2537                   momBin12 = fMomResC2SC->GetYaxis()->FindBin(qinv12);
2538                   momBin13 = fMomResC2SC->GetYaxis()->FindBin(qinv13);
2539                   momBin23 = fMomResC2SC->GetYaxis()->FindBin(qinv23);            
2540                   if(momBin12 >= 20) momBin12 = 19;
2541                   if(momBin13 >= 20) momBin13 = 19;
2542                   if(momBin23 >= 20) momBin23 = 19;
2543                   //
2544                   if(ch1==ch2) MomResCorr12 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin12);
2545                   else MomResCorr12 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin12);
2546                   if(ch1==ch3) MomResCorr13 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin13);
2547                   else MomResCorr13 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin13);
2548                   if(ch2==ch3) MomResCorr23 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin23);
2549                   else MomResCorr23 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin23);
2550                 }
2551                 
2552                 if(ENsum==0) {
2553                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fTerms3->Fill(q3); 
2554                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactor->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23));
2555                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactorWeighted->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23), MomResCorr12*MomResCorr13*MomResCorr23);
2556                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv12);
2557                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv13);
2558                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv23);
2559                 }
2560                 if(ENsum==6) {
2561                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTerms3->Fill(q3);
2562                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv12);
2563                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv13);
2564                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv23);
2565                 }
2566                 if(ENsum==3){
2567                   if(fill2) {
2568                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fTerms3->Fill(q3);
2569                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactor->Fill(q3, 1/(FSICorr12));
2570                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12);
2571                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv12);
2572                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv13);
2573                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv23);
2574                   }if(fill3) {
2575                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fTerms3->Fill(q3);
2576                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactor->Fill(q3, 1/(FSICorr12));
2577                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12);
2578                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv12);
2579                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv13);
2580                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv23);
2581                   }if(fill4) {
2582                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fTerms3->Fill(q3);
2583                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactor->Fill(q3, 1/(FSICorr12));
2584                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactorWeighted->Fill(q3, 1/(FSICorr12), MomResCorr12);
2585                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv12);
2586                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv13);
2587                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv23);
2588                   }
2589                 }
2590                 
2591                 // r3 denominator
2592                 if(ENsum==6 && ch1==ch2 && ch1==ch3){
2593                   Positive1stTripletWeights = kTRUE;
2594                   //
2595                   GetWeight(pVect1, pVect2, weight12, weight12Err);
2596                   GetWeight(pVect1, pVect3, weight13, weight13Err);
2597                   GetWeight(pVect2, pVect3, weight23, weight23Err);
2598                   
2599                   if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {// weight should never be larger than 1
2600                     if(fMbin==0 && bin1==0) {
2601                       ((TH1D*)fOutputList->FindObject("fTPNRejects3pion1"))->Fill(q3, sqrt(fabs(weight12*weight13*weight23)));
2602                     }
2603                   }else{
2604                     
2605                     Float_t MuonCorr12=1.0, MuonCorr13=1.0, MuonCorr23=1.0;
2606                     if(!fGenerateSignal && !fMCcase) {
2607                       MuonCorr12 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin12);
2608                       MuonCorr13 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin13);
2609                       MuonCorr23 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin23);
2610                     }
2611                     
2612                     // no MRC, no Muon Correction
2613                     weight12CC[0] = ((weight12+1) - ffcSq*FSICorr12 - (1-ffcSq));
2614                     weight12CC[0] /= FSICorr12*ffcSq;
2615                     weight13CC[0] = ((weight13+1) - ffcSq*FSICorr13 - (1-ffcSq));
2616                     weight13CC[0] /= FSICorr13*ffcSq;
2617                     weight23CC[0] = ((weight23+1) - ffcSq*FSICorr23 - (1-ffcSq));
2618                     weight23CC[0] /= FSICorr23*ffcSq;
2619                     if(weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2620                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, sqrt(weight12CC[0]*weight13CC[0]*weight23CC[0]));
2621                     }
2622                     // no Muon Correction
2623                     weight12CC[1] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2624                     weight12CC[1] /= FSICorr12*ffcSq;
2625                     weight13CC[1] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2626                     weight13CC[1] /= FSICorr13*ffcSq;
2627                     weight23CC[1] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2628                     weight23CC[1] /= FSICorr23*ffcSq;
2629                     if(weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2630                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(2, q3, sqrt(weight12CC[1]*weight13CC[1]*weight23CC[1]));
2631                     }
2632                     // both Corrections
2633                     weight12CC[2] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2634                     weight12CC[2] /= FSICorr12*ffcSq;
2635                     weight12CC[2] *= MuonCorr12;
2636                     weight13CC[2] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2637                     weight13CC[2] /= FSICorr13*ffcSq;
2638                     weight13CC[2] *= MuonCorr13;
2639                     weight23CC[2] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2640                     weight23CC[2] /= FSICorr23*ffcSq;
2641                     weight23CC[2] *= MuonCorr23;
2642                     
2643                     if(weight12CC[2] < 0 || weight13CC[2] < 0 || weight23CC[2] < 0) {// C2^QS can never be less than unity
2644                       if(fMbin==0 && bin1==0) {
2645                         ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC[2]*weight13CC[2]*weight23CC[2])));
2646                       }
2647                       if(weight12CC[2] < 0) weight12CC[2]=0;
2648                       if(weight13CC[2] < 0) weight13CC[2]=0;
2649                       if(weight23CC[2] < 0) weight23CC[2]=0;
2650                       Positive1stTripletWeights = kFALSE;
2651                     }
2652                     /////////////////////////////////////////////////////
2653                     weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]);
2654                     /////////////////////////////////////////////////////
2655                     if(Positive1stTripletWeights){
2656                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(3, q3, weightTotal);
2657                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(4, q3, 1);
2658                     }else{
2659                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNegNorm->Fill(4, q3, 1);
2660                     }
2661                     //
2662                     // Full Weight reconstruction
2663                     
2664                     for(Int_t RcohIndex=0; RcohIndex<2; RcohIndex++){// Rcoh=0, then Rcoh=Rch
2665                       for(Int_t GIndex=0; GIndex<50; GIndex++){
2666                         Int_t FillBin = 5 + RcohIndex*50 + GIndex;
2667                         Float_t G = 0.02*GIndex;
2668                         if(RcohIndex==0){
2669                           T12 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight12CC[2])) / (2*pow(1-G,2));
2670                           T13 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight13CC[2])) / (2*pow(1-G,2));
2671                           T23 = (-2*G*(1-G) + sqrt(pow(2*G*(1-G),2) + 4*pow(1-G,2)*weight23CC[2])) / (2*pow(1-G,2));
2672                           weightTotal = 2*G*(1-G)*(T12 + T13 + T23) + pow(1-G,2)*(T12*T12 + T13*T13 + T23*T23);
2673                           weightTotal += 2*G*pow(1-G,2)*(T12*T13 + T12*T23 + T13*T23) + 2*pow(1-G,3)*T12*T13*T23;
2674                         }else{
2675                           T12 = sqrt(weight12CC[2] / (1-G*G));
2676                           T13 = sqrt(weight13CC[2] / (1-G*G));
2677                           T23 = sqrt(weight23CC[2] / (1-G*G));
2678                           weightTotal = (1-G*G)*(T12*T12 + T13*T13 + T23*T23);
2679                           weightTotal += (6*G*pow(1-G,2) + 2*pow(1-G,3)) * T12*T13*T23;
2680                         }
2681                         if(Positive1stTripletWeights){
2682                           Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(FillBin, q3, weightTotal);
2683                         }else{
2684                           Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNegNorm->Fill(FillBin, q3, weightTotal);
2685                         }
2686                       }
2687                     }
2688                     //
2689                     /*weight12CC_e = weight12Err*MomResCorr12 / FSICorr12 / ffcSq * MuonCorr12;
2690                       weight13CC_e = weight13Err*MomResCorr13 / FSICorr13 / ffcSq * MuonCorr13;
2691                       weight23CC_e = weight23Err*MomResCorr23 / FSICorr23 / ffcSq * MuonCorr23;
2692                       if(weight12CC[2]*weight13CC[2]*weight23CC[2] > 0){
2693                       weightTotalErr = pow(2 * sqrt(3) * weight12CC_e*weight13CC[2]*weight23CC[2] / sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]),2);
2694                       }
2695                       weightTotalErr += pow(weight12CC_e,2) + pow(weight13CC_e,2) + pow(weight23CC_e,2);
2696                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNormErr->Fill(4, q3, weightTotalErr);*/
2697                     
2698                   }// 1st r3 den check
2699                   
2700                 }// r3 den
2701                 
2702
2703                 if(ch1==ch2 && ch1==ch3 && ENsum==0){
2704                   ((TH3D*)fOutputList->FindObject("fKT3DistTerm1"))->Fill(fMbin+1, KT3, q3);
2705                   if(q3<0.1){
2706                     Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2707                     Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2708                     Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2709                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt1);
2710                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt2);
2711                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt3);
2712                     if(fMbin==0){
2713                       ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt1);
2714                       ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt2);
2715                       ((TH3D*)fOutputList->FindObject("fQ3AvgpT"))->Fill(KT3index, q3, pt3);
2716                     }
2717                   }
2718                   
2719                 }
2720                 if(ch1==ch2 && ch1==ch3 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT3DistTerm5"))->Fill(fMbin+1, KT3, q3);
2721                 
2722                 
2723
2724                 
2725                 if(fMCcase && ENsum==6 && FilledMCpair12){// for momentum resolution and muon correction
2726                   if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
2727                     
2728                     pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
2729                     pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2730                     pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2731                     pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2732                     qinv13MC = GetQinv(pVect1MC, pVect3MC);
2733                     qinv23MC = GetQinv(pVect2MC, pVect3MC);
2734                     
2735                     q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
2736                     if(q3<0.1 && ch1==ch2 && ch1==ch3) ((TH2D*)fOutputList->FindObject("fQ3Res"))->Fill(KT3, q3-q3MC);
2737                     
2738                     Float_t TripletWeightTTC=1.0;// same-charge weights to mimic two-track depletion of same-charge pairs
2739                     //if(ch1==ch2 && qinv12>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv12);
2740                     //if(ch1==ch3 && qinv13>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv13);
2741                     //if(ch2==ch3 && qinv23>0.006) TripletWeightTTC *= SCpairWeight->Eval(qinv23);
2742                     
2743                     Int_t chGroup3[3]={ch1,ch2,ch3};
2744                     Float_t QinvMCGroup3[3]={qinv12MC,qinv13MC,qinv23MC};
2745                     //Float_t kTGroup3[3]={float(sqrt(pow(pVect1MC[1]+pVect2MC[1],2) + pow(pVect1MC[2]+pVect2MC[2],2))/2.),
2746                     //float(sqrt(pow(pVect1MC[1]+pVect3MC[1],2) + pow(pVect1MC[2]+pVect3MC[2],2))/2.),
2747                     //float(sqrt(pow(pVect2MC[1]+pVect3MC[1],2) + pow(pVect2MC[2]+pVect3MC[2],2))/2.)};
2748                     Float_t kTGroup3[3]={0};
2749                     
2750                     Pparent3[0]=pVect3MC[0]; Pparent3[1]=pVect3MC[1]; Pparent3[2]=pVect3MC[2]; Pparent3[3]=pVect3MC[3];
2751                     pionParent3=kFALSE;
2752                   
2753                     if(abs((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPdgCode)==13){// muon check
2754                       Int_t MotherLabel3 = (fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fMotherLabel;
2755                       if(abs((fEvt+en3)->fMCtracks[MotherLabel3].fPdgCode)==211) {
2756                         pionParent3=kTRUE;
2757                         Pparent3[1] = (fEvt+en3)->fMCtracks[MotherLabel3].fPx; Pparent3[2] = (fEvt+en3)->fMCtracks[MotherLabel3].fPy; Pparent3[3] = (fEvt+en3)->fMCtracks[MotherLabel3].fPz;
2758                         Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
2759                       }
2760                     }
2761                     
2762                     parentQinv13 = GetQinv(Pparent1, Pparent3);
2763                     parentQinv23 = GetQinv(Pparent2, Pparent3);
2764                     parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
2765                    
2766                     if(parentQinv12 > 0.001 && parentQinv13 > 0.001 && parentQinv23 > 0.001 && parentQ3 < 0.5){
2767                       FilledMCtriplet123=kTRUE;
2768                       if(pionParent1 || pionParent2 || pionParent3) {// want at least one pion-->muon
2769                         
2770                         Float_t parentQinvGroup3[3]={parentQinv12, parentQinv13, parentQinv23};
2771                         //Float_t parentkTGroup3[3]={float(sqrt(pow(Pparent1[1]+Pparent2[1],2) + pow(Pparent1[2]+Pparent2[2],2))/2.),
2772                         //float(sqrt(pow(Pparent1[1]+Pparent3[1],2) + pow(Pparent1[2]+Pparent3[2],2))/2.),
2773                         //float(sqrt(pow(Pparent2[1]+Pparent3[1],2) + pow(Pparent2[2]+Pparent3[2],2))/2.)};
2774                         Float_t parentkTGroup3[3]={0};
2775                         
2776                         ((TH2D*)fOutputList->FindObject("fAvgQ12VersusQ3"))->Fill(parentQ3, parentQinv12);
2777                         ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13);
2778                         ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23);
2779
2780                         for(Int_t term=1; term<=4; term++){
2781                           if(term==1) {}
2782                           else if(term==2) {if(!pionParent1 && !pionParent2) continue;}
2783                           else if(term==3) {if(!pionParent1 && !pionParent3) continue;}
2784                           else {if(!pionParent2 && !pionParent3) continue;}
2785                           for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2786                             Float_t Rvalue = 5+Riter;
2787                             Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
2788                             Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
2789                             Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
2790                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(1, Rvalue, q3MC, WInput*TripletWeightTTC);
2791                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ3, WInput*TripletWeightTTC);
2792                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(1, Rvalue, q3MC, WInputFSI*TripletWeightTTC);
2793                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(1, Rvalue, parentQ3, WInputParentFSI*TripletWeightTTC);
2794                             //
2795                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(2, Rvalue, q3MC, TripletWeightTTC);
2796                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ3, TripletWeightTTC);
2797                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(2, Rvalue, q3MC, TripletWeightTTC);
2798                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(2, Rvalue, parentQ3, TripletWeightTTC);
2799                           }// Riter
2800                         }// term loop
2801                     
2802                       }// pion parent check
2803                     }// parentQ check (muon correction)
2804                     
2805                     // 3-pion momentum resolution
2806                     for(Int_t term=1; term<=5; term++){
2807                       for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2808                         Float_t Rvalue = 5+Riter;
2809                         Float_t WInput = MCWeight3(term, Rvalue, ffcSqMRC, chGroup3, QinvMCGroup3, kTGroup3);
2810                         Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fIdeal->Fill(Rvalue, q3MC, WInput*TripletWeightTTC);
2811                         Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fSmeared->Fill(Rvalue, q3, WInput*TripletWeightTTC);
2812                       }
2813                     }
2814                     
2815                   }// 3rd particle label check
2816                 }// MCcase and ENsum==6
2817                 
2818                 
2819                 
2820
2821                 /////////////////////////////////////////////////////////////
2822                 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2823                   if(en4==0){
2824                     if(fLowQPairSwitch_E0E0[i]->At(l)=='0') continue;
2825                     if(fLowQPairSwitch_E0E0[j]->At(l)=='0') continue;
2826                     if(fLowQPairSwitch_E0E0[k]->At(l)=='0') continue;
2827                   }else if(en4==1){
2828                     if(en3==0){
2829                       if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2830                       if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2831                       if(fLowQPairSwitch_E0E1[k]->At(l)=='0') continue;
2832                     }else{ 
2833                       if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2834                       if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2835                       if(fLowQPairSwitch_E1E1[k]->At(l)=='0') continue;
2836                     }
2837                   }else if(en4==2){
2838                     if(fLowQPairSwitch_E0E2[i]->At(l)=='0') continue;
2839                     if(fLowQPairSwitch_E0E2[j]->At(l)=='0') continue;
2840                     if(fLowQPairSwitch_E1E2[k]->At(l)=='0') continue;
2841                   }else{
2842                     if(fLowQPairSwitch_E0E3[i]->At(l)=='0') continue;
2843                     if(fLowQPairSwitch_E1E3[j]->At(l)=='0') continue;
2844                     if(fLowQPairSwitch_E2E3[k]->At(l)=='0') continue;
2845                   }
2846
2847                   pVect4[0]=(fEvt+en4)->fTracks[l].fEaccepted;
2848                   pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2849                   pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2850                   pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2851                   ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2852                   qinv14 = GetQinv(pVect1, pVect4);
2853                   qinv24 = GetQinv(pVect2, pVect4);
2854                   qinv34 = GetQinv(pVect3, pVect4);
2855                   q4 = sqrt(pow(q3,2) + pow(qinv14,2) + pow(qinv24,2) + pow(qinv34,2));
2856                   
2857                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2858                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(1, qinv12); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(2, qinv13); 
2859                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(3, qinv14); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(4, qinv23); 
2860                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(5, qinv24); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(6, qinv34);
2861                   }
2862                   
2863                   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.;
2864                   if(KT4<=fKT4transition) KT4index=0;
2865                   else KT4index=1;
2866                   
2867                   FSICorr14 = FSICorrelation(ch1,ch4, qinv14);
2868                   FSICorr24 = FSICorrelation(ch2,ch4, qinv24);
2869                   FSICorr34 = FSICorrelation(ch3,ch4, qinv34);
2870                   
2871                   if(!fGenerateSignal && !fMCcase) {
2872                     momBin14 = fMomResC2SC->GetYaxis()->FindBin(qinv14);
2873                     momBin24 = fMomResC2SC->GetYaxis()->FindBin(qinv24);
2874                     momBin34 = fMomResC2SC->GetYaxis()->FindBin(qinv34);                  
2875                     if(momBin14 >= 20) momBin14 = 19;
2876                     if(momBin24 >= 20) momBin24 = 19;
2877                     if(momBin34 >= 20) momBin34 = 19;
2878                     //
2879                     if(ch1==ch4) MomResCorr14 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin14);
2880                     else MomResCorr14 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin14);
2881                     if(ch2==ch4) MomResCorr24 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin24);
2882                     else MomResCorr24 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin24);
2883                     if(ch3==ch4) MomResCorr34 = fMomResC2SC->GetBinContent(rBinForTPNMomRes, momBin34);
2884                     else MomResCorr34 = fMomResC2MC->GetBinContent(rBinForTPNMomRes, momBin34);
2885                   }
2886                   
2887                   Bool_t FillTerms[13]={kFALSE};
2888                   SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, ENsum, FillTerms);
2889                   //
2890                   for(int ft=0; ft<13; ft++) {
2891                     Float_t FSIfactor = 1.0;
2892                     Float_t MomResWeight = 1.0;
2893                     if(ft==0) {
2894                       FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr14 * FSICorr23 * FSICorr24 * FSICorr34);
2895                       MomResWeight = MomResCorr12 * MomResCorr13 * MomResCorr14 * MomResCorr23 * MomResCorr24 * MomResCorr34;
2896                     }else if(ft<=4) {
2897                       FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr23);
2898                       MomResWeight = MomResCorr12 * MomResCorr13 * MomResCorr23;
2899                     }else if(ft<=10) {
2900                       FSIfactor = 1/(FSICorr12);
2901                       MomResWeight = MomResCorr12;
2902                     }else if(ft==11) {
2903                       FSIfactor = 1/(FSICorr12 * FSICorr34);
2904                       MomResWeight = MomResCorr12 * MomResCorr34;
2905                     }else {FSIfactor = 1.0; MomResWeight = 1.0;}
2906                     if(FillTerms[ft]) {
2907                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fTerms4->Fill(q4);
2908                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactor->Fill(q4, FSIfactor);
2909                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactorWeighted->Fill(q4, FSIfactor, MomResWeight);
2910                     }
2911                   }
2912                   
2913                   /////////////////////////////////////////////////////////////
2914                   // r4{2}
2915                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2916                     Positive2ndTripletWeights=kTRUE;
2917                     //
2918                     GetWeight(pVect1, pVect4, weight14, weight14Err);
2919                     GetWeight(pVect2, pVect4, weight24, weight24Err);
2920                     GetWeight(pVect3, pVect4, weight34, weight34Err);
2921                     
2922                     Float_t MuonCorr14=1.0, MuonCorr24=1.0, MuonCorr34=1.0;
2923                     if(!fGenerateSignal && !fMCcase) {
2924                       MuonCorr14 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin14);
2925                       MuonCorr24 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin24);
2926                       MuonCorr34 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin34);
2927                     }
2928                     
2929                     // no MRC, no Muon Correction
2930                     weight14CC[0] = ((weight14+1) - ffcSq*FSICorr14 - (1-ffcSq));
2931                     weight14CC[0] /= FSICorr14*ffcSq;
2932                     weight24CC[0] = ((weight24+1) - ffcSq*FSICorr24 - (1-ffcSq));
2933                     weight24CC[0] /= FSICorr24*ffcSq;
2934                     weight34CC[0] = ((weight34+1) - ffcSq*FSICorr34 - (1-ffcSq));
2935                     weight34CC[0] /= FSICorr34*ffcSq;
2936                     if(weight14CC[0] > 0 && weight24CC[0] > 0 && weight34CC[0] > 0 && weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2937                       weightTotal  = sqrt(weight12CC[0]*weight13CC[0]*weight24CC[0]*weight34CC[0]);
2938                       weightTotal += sqrt(weight12CC[0]*weight14CC[0]*weight23CC[0]*weight34CC[0]);
2939                       weightTotal += sqrt(weight13CC[0]*weight14CC[0]*weight23CC[0]*weight24CC[0]);
2940                       weightTotal /= 3.;
2941                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal);
2942                     }
2943                     // no Muon Correction
2944                     weight14CC[1] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2945                     weight14CC[1] /= FSICorr14*ffcSq;
2946             &n