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