414a7def91d13a61620d603b3e1552766b2df3ef
[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                 }// term=4
1203                 
1204               }// term_3
1205               
1206               for(Int_t c4=0; c4<2; c4++){
1207                 for(Int_t term=0; term<13; term++){
1208                                 
1209                   TString *namePC4 = new TString("FourParticle_Charge1_");
1210                   *namePC4 += c1;
1211                   namePC4->Append("_Charge2_");
1212                   *namePC4 += c2;
1213                   namePC4->Append("_Charge3_");
1214                   *namePC4 += c3;
1215                   namePC4->Append("_Charge4_");
1216                   *namePC4 += c4;
1217                   namePC4->Append("_M_");
1218                   *namePC4 += mb;
1219                   namePC4->Append("_ED_");
1220                   *namePC4 += edB;
1221                   namePC4->Append("_Term_");
1222                   *namePC4 += term+1;
1223                   
1224                   ///////////////////////////////////////
1225                   // skip degenerate histograms
1226                   if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
1227                   if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
1228                   if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
1229                   /////////////////////////////////////////
1230                  
1231                   TString *nameNorm=new TString(namePC4->Data());
1232                   nameNorm->Append("_Norm");
1233                   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);
1234                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4);
1235                   //
1236                   TString *name1DQ=new TString(namePC4->Data());
1237                   name1DQ->Append("_1D");
1238                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = new TH1D(name1DQ->Data(),"", fQbinsQ4,0,fQupperBoundQ4);
1239                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4);
1240                   //
1241                   TString *nameKfactor=new TString(namePC4->Data());
1242                   nameKfactor->Append("_Kfactor");
1243                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
1244                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor);
1245                   //
1246                   if(c1==c2 && c1==c3 && c1==c4 && term==12 ){
1247                     TString *nameTwoPartNorm=new TString(namePC4->Data());
1248                     nameTwoPartNorm->Append("_TwoPartNorm");
1249                     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);
1250                     fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm);
1251                   }
1252                   
1253                   if(fMCcase==kTRUE){
1254                     // Momentum resolution correction histos
1255                     TString *nameMomResIdeal=new TString(namePC4->Data());
1256                     nameMomResIdeal->Append("_Ideal");
1257                     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);
1258                     if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal);
1259                     TString *nameMomResSmeared=new TString(namePC4->Data());
1260                     nameMomResSmeared->Append("_Smeared");
1261                     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);
1262                     if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared);
1263                     // Muon correction histos
1264                     TString *nameMuonIdeal=new TString(namePC4->Data());
1265                     nameMuonIdeal->Append("_MuonIdeal");
1266                     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);
1267                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal);
1268                     TString *nameMuonSmeared=new TString(namePC4->Data());
1269                     nameMuonSmeared->Append("_MuonSmeared");
1270                     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);
1271                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared);
1272                     //
1273                     TString *nameMuonPionK4=new TString(namePC4->Data());
1274                     nameMuonPionK4->Append("_MuonPionK4");
1275                     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);
1276                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4);
1277                     //
1278                     TString *namePionPionK4=new TString(namePC4->Data());
1279                     namePionPionK4->Append("_PionPionK4");
1280                     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);
1281                     if(mb==0 && edB==0 && term<12) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4);
1282                     
1283                   }// MCcase
1284                   
1285                   
1286                 }
1287               }
1288
1289             }//c3
1290           }//c2
1291         }//c1
1292       }// ED
1293     }// mbin
1294   }// Pdensity Method
1295
1296   
1297   if(fTabulatePairs){
1298     
1299     for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
1300       for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
1301         for(Int_t mb=0; mb<fMbins; mb++){
1302           for(Int_t edB=0; edB<fEDbins; edB++){
1303             
1304             TString *nameNum = new TString("TPN_num_Kt_");
1305             *nameNum += tKbin;
1306             nameNum->Append("_Ky_");
1307             *nameNum += yKbin;
1308             nameNum->Append("_M_");
1309             *nameNum += mb;
1310             nameNum->Append("_ED_");
1311             *nameNum += edB;
1312             
1313             TString *nameDen = new TString("TPN_den_Kt_");
1314             *nameDen += tKbin;
1315             nameDen->Append("_Ky_");
1316             *nameDen += yKbin;
1317             nameDen->Append("_M_");
1318             *nameDen += mb;
1319             nameDen->Append("_ED_");
1320             *nameDen += edB;
1321             
1322             if(edB==0){
1323               KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = new TH3D(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1324               fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD);
1325               
1326               KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = new TH3D(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
1327               fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD);
1328             }   
1329           
1330           }
1331         }
1332       }
1333     }
1334     
1335   }
1336   
1337   
1338   TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
1339   fOutputList->Add(fQsmearMean);
1340   TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
1341   fOutputList->Add(fQsmearSq);
1342   TH2D *fQ2Res = new TH2D("fQ2Res","",20,0,1, 200,-.2,.2);
1343   fOutputList->Add(fQ2Res);
1344   TH2D *fQ3Res = new TH2D("fQ3Res","",20,0,1, 200,-.3,.3);
1345   fOutputList->Add(fQ3Res);
1346   TH2D *fQ4Res = new TH2D("fQ4Res","",20,0,1, 200,-.4,.4);
1347   fOutputList->Add(fQ4Res);
1348   
1349   TH2D *DistQinv4pion = new TH2D("DistQinv4pion","",6,0.5,6.5, 20,0,0.1);
1350   fOutputList->Add(DistQinv4pion);
1351   TH2D *DistQinvMC4pion = new TH2D("DistQinvMC4pion","",6,0.5,6.5, 20,0,0.1);
1352   if(fMCcase) fOutputList->Add(DistQinvMC4pion);
1353
1354   TH2D *fAvgQ12VersusQ3 = new TH2D("fAvgQ12VersusQ3","",10,0,0.1, 20,0,0.1);
1355   fOutputList->Add(fAvgQ12VersusQ3);
1356   TH2D *fAvgQ13VersusQ3 = new TH2D("fAvgQ13VersusQ3","",10,0,0.1, 20,0,0.1);
1357   fOutputList->Add(fAvgQ13VersusQ3);
1358   TH2D *fAvgQ23VersusQ3 = new TH2D("fAvgQ23VersusQ3","",10,0,0.1, 20,0,0.1);
1359   fOutputList->Add(fAvgQ23VersusQ3);
1360
1361   TH1D *fDistPionParents4 = new TH1D("fDistPionParents4","",4,0.5,4.5);
1362   fOutputList->Add(fDistPionParents4);
1363
1364   ////////////////////////////////////
1365   ///////////////////////////////////  
1366   
1367   PostData(1, fOutputList);
1368 }
1369
1370 //________________________________________________________________________
1371 void AliFourPion::UserExec(Option_t *) 
1372 {
1373   // Main loop
1374   // Called for each event
1375   //cout<<"===========  Event # "<<fEventCounter+1<<"  ==========="<<endl;
1376   fEventCounter++;
1377   
1378   if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
1379   
1380   fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
1381   if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
1382   
1383   
1384   // Trigger Cut
1385   if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
1386   Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
1387   if(!isSelected1 && !fMCcase) {return;}
1388   }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
1389     Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
1390     Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
1391     if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
1392   }else {return;}
1393
1394   ///////////////////////////////////////////////////////////
1395   const AliAODVertex *primaryVertexAOD;
1396   AliCentrality *centrality;// for AODs and ESDs
1397
1398  
1399   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1400   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1401   fPIDResponse = inputHandler->GetPIDResponse();
1402
1403   
1404   TClonesArray *mcArray = 0x0;
1405   if(fMCcase){
1406     if(fAODcase){ 
1407       mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
1408       if(!mcArray || mcArray->GetEntriesFast() >= kMCarrayLimit){
1409         cout<<"No MC particle branch found or Array too large!!"<<endl;
1410         return;
1411       }
1412     }
1413   }
1414   
1415
1416   UInt_t status=0;
1417   Int_t positiveTracks=0, negativeTracks=0;
1418   Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
1419    
1420   Double_t vertex[3]={0};
1421   Int_t zbin=0;
1422   Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
1423   /////////////////////////////////////////////////
1424
1425   
1426   Float_t centralityPercentile=0;
1427   Float_t cStep=5.0, cStart=0;
1428   
1429  
1430   if(fAODcase){// AOD case
1431     
1432     if(fPbPbcase){
1433       centrality = fAOD->GetCentrality();
1434       centralityPercentile = centrality->GetCentralityPercentile("V0M");
1435       if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
1436       if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range.  Skipping Event"<<endl;*/ return;}
1437       cout<<"Centrality % = "<<centralityPercentile<<endl;
1438     }
1439     
1440     
1441     ((TH1F*)fOutputList->FindObject("fMultDist0"))->Fill(fAOD->GetNumberOfTracks());
1442
1443     // Pile-up rejection
1444     AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
1445     if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
1446     else AnaUtil->SetUseMVPlpSelection(kFALSE);
1447     Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD); 
1448     if(pileUpCase) return;
1449     
1450     ////////////////////////////////
1451     // Vertexing
1452     ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
1453     primaryVertexAOD = fAOD->GetPrimaryVertex();
1454     vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
1455     
1456     if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut 
1457     ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
1458     
1459     if(!fMCcase && primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
1460    
1461     ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
1462  
1463     fBfield = fAOD->GetMagneticField();
1464     
1465     for(Int_t i=0; i<fZvertexBins; i++){
1466       if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
1467         zbin=i;
1468         break;
1469       }
1470     }
1471
1472    
1473        
1474     /////////////////////////////
1475     // Create Shuffled index list
1476     Int_t randomIndex[fAOD->GetNumberOfTracks()];
1477     for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
1478     Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
1479     /////////////////////////////
1480   
1481    
1482     // Track loop
1483     for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
1484       AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
1485       if (!aodtrack) continue;
1486       if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1487     
1488       status=aodtrack->GetStatus();
1489       
1490       if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
1491       ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
1492       ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
1493       if(aodtrack->GetTPCNcls() < fMinTPCncls) continue;// TPC nCluster cut
1494       if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
1495
1496       if(fFilterBit != 7){
1497         Bool_t goodTrackOtherFB = kFALSE;
1498         for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1499           AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
1500           if(!aodtrack2) continue;
1501           if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
1502           
1503           if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
1504           
1505         }
1506         if(!goodTrackOtherFB) continue;
1507       }
1508       
1509       
1510       if(aodtrack->Pt() < 0.16) continue;
1511       if(fabs(aodtrack->Eta()) > 0.8) continue;
1512      
1513       
1514       Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
1515       if(!goodMomentum) continue; 
1516       aodtrack->GetXYZ( fTempStruct[myTracks].fX);
1517       
1518          
1519       Double_t dca2[2]={0};
1520       dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
1521       dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
1522       Double_t dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
1523              
1524       fTempStruct[myTracks].fStatus = status;
1525       fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
1526       fTempStruct[myTracks].fId = aodtrack->GetID();
1527       //
1528       fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
1529       fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1530       if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1531       fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1532       fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1533       fTempStruct[myTracks].fEta = aodtrack->Eta();
1534       fTempStruct[myTracks].fCharge = aodtrack->Charge();
1535       fTempStruct[myTracks].fDCAXY = dca2[0];
1536       fTempStruct[myTracks].fDCAZ = dca2[1];
1537       fTempStruct[myTracks].fDCA = dca3d;
1538       fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
1539       fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
1540       
1541     
1542       
1543       if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
1544       //if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
1545       //if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
1546      
1547      
1548       // PID section
1549       fTempStruct[myTracks].fElectron = kFALSE;
1550       fTempStruct[myTracks].fPion = kFALSE;
1551       fTempStruct[myTracks].fKaon = kFALSE;
1552       fTempStruct[myTracks].fProton = kFALSE;
1553       
1554       Float_t nSigmaTPC[5];
1555       Float_t nSigmaTOF[5];
1556       nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
1557       nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
1558       fTempStruct[myTracks].fTOFhit = kFALSE;// default
1559       Float_t signalTPC=0, signalTOF=0;
1560       Double_t integratedTimesTOF[10]={0};
1561
1562     
1563       Bool_t DoPIDWorkAround=kTRUE;
1564       //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
1565       if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
1566       if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
1567         nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
1568         nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
1569         nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
1570         nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
1571         nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
1572         //
1573         nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
1574         nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
1575         nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
1576         nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
1577         nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton);
1578         signalTPC = aodtrack->GetTPCsignal();
1579         if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1580           fTempStruct[myTracks].fTOFhit = kTRUE;
1581           signalTOF = aodtrack->GetTOFsignal();
1582           aodtrack->GetIntegratedTimes(integratedTimesTOF);
1583         }else fTempStruct[myTracks].fTOFhit = kFALSE;
1584         
1585       }else {// FilterBit 7 PID workaround
1586         
1587         for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
1588           AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
1589           if (!aodTrack2) continue;
1590           if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
1591           
1592           UInt_t status2=aodTrack2->GetStatus();
1593           
1594           nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
1595           nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
1596           nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
1597           nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
1598           nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
1599           //
1600           nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
1601           nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
1602           nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
1603           nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
1604           nSigmaTOF[4]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton);
1605           signalTPC = aodTrack2->GetTPCsignal();
1606           
1607           if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
1608             fTempStruct[myTracks].fTOFhit = kTRUE;
1609             signalTOF = aodTrack2->GetTOFsignal();
1610             aodTrack2->GetIntegratedTimes(integratedTimesTOF);
1611           }else fTempStruct[myTracks].fTOFhit = kFALSE;
1612           
1613         }// aodTrack2
1614       }// FilterBit 7 PID workaround
1615      
1616      
1617       ///////////////////
1618       ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
1619       if(fTempStruct[myTracks].fTOFhit) {
1620         ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
1621       }
1622       ///////////////////
1623       
1624       // Use TOF if good hit and above threshold
1625       if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
1626         if(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1627         if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1628         if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1629         if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1630       }else {// TPC info instead
1631         if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
1632         if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
1633         if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
1634         if(fabs(nSigmaTPC[4])<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
1635       }
1636       
1637     
1638       // Ensure there is only 1 candidate per track
1639       if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
1640       if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
1641       if(!fTempStruct[myTracks].fPion) continue;// only pions
1642       if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
1643       if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
1644       if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
1645       //if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;// superfluous
1646       ////////////////////////
1647       //if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons// superfluous
1648
1649
1650    
1651       if(fTempStruct[myTracks].fCharge==+1) {
1652         ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1653         ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1654       }else {
1655         ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
1656         ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
1657       }
1658      
1659       ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
1660       ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
1661       
1662      
1663     
1664       if(fTempStruct[myTracks].fPion) {// pions
1665         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2)); 
1666         fTempStruct[myTracks].fKey = 1;
1667       }else if(fTempStruct[myTracks].fKaon){// kaons
1668         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
1669         fTempStruct[myTracks].fKey = 10;
1670       }else{// protons
1671         fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
1672         fTempStruct[myTracks].fKey = 100;
1673       }
1674       
1675            
1676
1677       if(aodtrack->Charge() > 0) positiveTracks++;
1678       else negativeTracks++;
1679       
1680       if(fTempStruct[myTracks].fPion) pionCount++;
1681       if(fTempStruct[myTracks].fKaon) kaonCount++;
1682       if(fTempStruct[myTracks].fProton) protonCount++;
1683
1684       myTracks++;
1685       
1686       if(fMCcase){// muon mothers
1687         AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
1688         if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
1689           AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
1690           if(parent->IsPhysicalPrimary()){
1691             ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
1692           }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
1693         }
1694         ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
1695       }
1696     }
1697     //cout<<"kinkcount = "<<kinkcount<<"   pionkinks = "<<pionkinks<<"   primarypionkinks = "<<primarypionkinks<<endl;
1698   }else {// ESD tracks
1699     cout<<"ESDs not supported currently"<<endl;
1700     return;
1701   }
1702   
1703   // Generator info only
1704   if(fMCcase && fGeneratorOnly){
1705     myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
1706     for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
1707       if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
1708       if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
1709       
1710       AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
1711       if(!mcParticle) continue;
1712       if(fabs(mcParticle->Eta())>0.8) continue;
1713       if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
1714       if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
1715       if(!mcParticle->IsPrimary()) continue;
1716       if(!mcParticle->IsPhysicalPrimary()) continue;
1717       if(abs(mcParticle->GetPdgCode())!=211) continue;
1718       
1719       fTempStruct[myTracks].fP[0] = mcParticle->Px();
1720       fTempStruct[myTracks].fP[1] = mcParticle->Py();
1721       fTempStruct[myTracks].fP[2] = mcParticle->Pz();
1722       fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
1723       
1724       fTempStruct[myTracks].fId = myTracks;// use my track counter 
1725       fTempStruct[myTracks].fLabel = mctrackN;
1726       fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
1727       if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
1728       fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
1729       fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
1730       fTempStruct[myTracks].fEta = mcParticle->Eta();
1731       fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
1732       fTempStruct[myTracks].fDCAXY = 0.;
1733       fTempStruct[myTracks].fDCAZ = 0.;
1734       fTempStruct[myTracks].fDCA = 0.;
1735       fTempStruct[myTracks].fPion = kTRUE;
1736       fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2)); 
1737       fTempStruct[myTracks].fKey = 1;
1738       
1739       myTracks++;
1740       pionCount++;
1741     }
1742   }
1743   
1744   if(myTracks >= 1) {
1745     ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
1746   }
1747  
1748  
1749   //cout<<"There are "<<myTracks<<"  myTracks"<<endl;
1750   //cout<<"pionCount = "<<pionCount<<"   kaonCount = "<<kaonCount<<"   protonCount = "<<protonCount<<endl;
1751   
1752
1753   /////////////////////////////////////////
1754   // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
1755   if(myTracks < 4) {cout<<"Less than 4 tracks. Skipping Event."<<endl; return;}
1756   /////////////////////////////////////////
1757  
1758
1759   ////////////////////////////////
1760   ///////////////////////////////
1761   // Mbin determination
1762   //
1763   // Mbin set to Pion Count Only for pp!!!!!!!
1764   fMbin=-1;
1765   if(!fPbPbcase){
1766     for(Int_t i=0; i<kMultBinspp; i++){
1767       if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
1768       // Mbin 0 has 1 pion
1769     }
1770   }else{
1771     for(Int_t i=0; i<fCentBins; i++){
1772       if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
1773         fMbin=i;// 0 = most central
1774         break;
1775       }
1776     }
1777   }
1778   
1779   if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
1780   
1781   ///////////////////
1782   // can only be called after fMbin has been set
1783   // Radius parameter only matters for Monte-Carlo data
1784   SetFSIindex(fRMax);
1785   ///////////////////  
1786   
1787   Int_t rBinForTPNMomRes = 10;
1788   if(fMbin==0) {rBinForTPNMomRes=10;}// 10 fm with EW (fRMax should be 11 for normal running)
1789   else if(fMbin==1) {rBinForTPNMomRes=9;}
1790   else if(fMbin<=3) {rBinForTPNMomRes=8;}
1791   else if(fMbin<=5) {rBinForTPNMomRes=7;}
1792   else {rBinForTPNMomRes=6;}
1793
1794   //////////////////////////////////////////////////
1795   fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
1796   //////////////////////////////////////////////////
1797   
1798
1799   
1800   ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
1801   ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
1802
1803   ////////////////////////////////////
1804   // Add event to buffer if > 0 tracks
1805   if(myTracks > 0){
1806     fEC[zbin][fMbin]->FIFOShift();
1807     (fEvt) = fEC[zbin][fMbin]->fEvtStr;
1808     (fEvt)->fNtracks = myTracks;
1809     (fEvt)->fFillStatus = 1;
1810     for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
1811     if(fMCcase){
1812       (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
1813       for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
1814         AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
1815         (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
1816         (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
1817         (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
1818         (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
1819         (fEvt)->fMCtracks[i].fPdgCode = tempMCTrack->GetPdgCode();
1820         (fEvt)->fMCtracks[i].fMotherLabel = tempMCTrack->GetMother();
1821       } 
1822     }
1823   }
1824     
1825   
1826   
1827   Float_t qinv12=0, qinv13=0, qinv14=0, qinv23=0, qinv24=0, qinv34=0;
1828   Float_t qout=0, qside=0, qlong=0;
1829   Float_t kT12=0;
1830   Float_t q3=0, q3MC=0;
1831   Float_t q4=0, q4MC=0;
1832   Int_t ch1=0, ch2=0, ch3=0, ch4=0;
1833   Int_t bin1=0, bin2=0, bin3=0, bin4=0;
1834   Float_t pVect1[4]={0}; 
1835   Float_t pVect2[4]={0};
1836   Float_t pVect3[4]={0};
1837   Float_t pVect4[4]={0};
1838   Float_t pVect1MC[4]={0}; 
1839   Float_t pVect2MC[4]={0};
1840   Float_t pVect3MC[4]={0};
1841   Float_t pVect4MC[4]={0};
1842   Float_t Pparent1[4]={0};
1843   Float_t Pparent2[4]={0};
1844   Float_t Pparent3[4]={0};
1845   Float_t Pparent4[4]={0};
1846   Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0;
1847   Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0;
1848   Float_t weight12CC[3]={0};
1849   Float_t weight13CC[3]={0};
1850   Float_t weight14CC[3]={0};
1851   Float_t weight23CC[3]={0};
1852   Float_t weight24CC[3]={0};
1853   Float_t weight34CC[3]={0};
1854   Float_t weightTotal=0;//, weightTotalErr=0;
1855   Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0; 
1856   Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0;
1857   Float_t parentQ3=0;
1858   Float_t FSICorr12=0, FSICorr13=0, FSICorr14=0, FSICorr23=0, FSICorr24=0, FSICorr34=0;
1859   Bool_t pionParent1=kFALSE, pionParent2=kFALSE, pionParent3=kFALSE, pionParent4=kFALSE;
1860   Bool_t FilledMCpair12=kFALSE, FilledMCtriplet123=kFALSE;
1861   Bool_t GoodTripletWeights=kFALSE;
1862   //
1863   AliAODMCParticle *mcParticle1=0x0;
1864   AliAODMCParticle *mcParticle2=0x0;
1865   
1866
1867   ////////////////////
1868   //Int_t PairCount[7]={0};
1869   //Int_t NormPairCount[7]={0};
1870   Int_t KT3index=0, KT4index=0;
1871
1872   // reset to defaults
1873   for(Int_t i=0; i<kMultLimitPbPb; i++) {
1874     fLowQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1875     fLowQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1876     fLowQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1877     fLowQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1878     fLowQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1879     fLowQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1880     fLowQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1881     fLowQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1882     //
1883     fNormQPairSwitch_E0E0[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1884     fNormQPairSwitch_E0E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1885     fNormQPairSwitch_E0E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1886     fNormQPairSwitch_E0E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1887     fNormQPairSwitch_E1E1[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1888     fNormQPairSwitch_E1E2[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1889     fNormQPairSwitch_E1E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1890     fNormQPairSwitch_E2E3[i]->Set(kMultLimitPbPb,fDefaultsCharSwitch);
1891   }
1892  
1893   
1894   //////////////////////////////////////////
1895   // make low-q pair storage and normalization-pair storage
1896   // 
1897   for(Int_t en1=0; en1<=2; en1++){// 1st event number (en1=0 is the same event as current event)
1898     for(Int_t en2=en1; en2<=3; en2++){// 2nd event number (en2=0 is the same event as current event)
1899       if(en1>1 && en1==en2) continue;
1900       
1901       for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {// 1st particle
1902         for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
1903           
1904           
1905           pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
1906           pVect1[1]=(fEvt+en1)->fTracks[i].fP[0];      pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
1907           pVect1[2]=(fEvt+en1)->fTracks[i].fP[1];      pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
1908           pVect1[3]=(fEvt+en1)->fTracks[i].fP[2];      pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
1909           ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
1910           ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
1911           
1912           qinv12 = GetQinv(pVect1, pVect2);
1913           kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
1914           SetFillBins2(ch1, ch2, bin1, bin2);
1915           
1916           if(qinv12 < fQLowerCut && !fMCcase) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
1917           if(ch1 == ch2 && !fGeneratorOnly && !fMCcase){
1918             if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
1919               if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
1920               continue;
1921             }
1922           }
1923           
1924           GetQosl(pVect1, pVect2, qout, qside, qlong);
1925           if( (en1+en2==0)) {
1926             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12);
1927             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2QW->Fill(kT12, qinv12, qinv12);
1928             // osl frame
1929             if((kT12 > 0.2) && (kT12 < 0.3)){
1930               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1931               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1932             }
1933             if((kT12 > 0.6) && (kT12 < 0.7)){  
1934               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1935               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1936             }
1937             // unit mult bins
1938             if( (fEvt+en1)->fNtracks%100==0){
1939               Int_t kTindex=0;
1940               if(kT12>0.3) kTindex=1;
1941               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
1942               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[0].fUnitMultBin->Fill(UnitMultBin, qinv12);
1943             }
1944             
1945           }
1946           if( (en1+en2==1)) {
1947             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
1948             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2QW->Fill(kT12, qinv12, qinv12);
1949             // osl frame
1950             if((kT12 > 0.2) && (kT12 < 0.3)){  
1951               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1952               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1953             }
1954             if((kT12 > 0.6) && (kT12 < 0.7)){  
1955               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
1956               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
1957             }
1958             // unit mult bins
1959             if( (fEvt+en1)->fNtracks%100==0){
1960               Int_t kTindex=0;
1961               if(kT12>0.3) kTindex=1;
1962               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
1963               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[1].fUnitMultBin->Fill(UnitMultBin, qinv12);
1964             }
1965           }
1966           //////////////////////////////////////////
1967           if(fTabulatePairs && en1==0 && en2<=1 && bin1==bin2){
1968             Float_t kY = 0;
1969             Int_t kTbin=-1, kYbin=-1;
1970             
1971             for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(kT12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {kTbin = kIt; break;}} 
1972             for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(kY < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {kYbin = kIt; break;}}
1973             if((kTbin<0) || (kYbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1974             if((kTbin>=fKbinsT) || (kYbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
1975             if(fGenerateSignal && en2==0) {
1976               Int_t chGroup2[2]={ch1,ch2};
1977               Float_t WInput = MCWeight(chGroup2, fRMax, 0.65, qinv12, kT12);
1978               KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
1979             }else KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
1980           }
1981           
1982           //////////////////////////////////////////////////////////////////////////////
1983           
1984           if(qinv12 <= fQcut) {
1985             if(en1==0 && en2==0) {fLowQPairSwitch_E0E0[i]->AddAt('1',j);}
1986             if(en1==0 && en2==1) {fLowQPairSwitch_E0E1[i]->AddAt('1',j);}
1987             if(en1==0 && en2==2) {fLowQPairSwitch_E0E2[i]->AddAt('1',j);}
1988             if(en1==0 && en2==3) {fLowQPairSwitch_E0E3[i]->AddAt('1',j);}
1989             if(en1==1 && en2==1) {fLowQPairSwitch_E1E1[i]->AddAt('1',j);}
1990             if(en1==1 && en2==2) {fLowQPairSwitch_E1E2[i]->AddAt('1',j);}
1991             if(en1==1 && en2==3) {fLowQPairSwitch_E1E3[i]->AddAt('1',j);}
1992             if(en1==2 && en2==3) {fLowQPairSwitch_E2E3[i]->AddAt('1',j);}
1993           }
1994           if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) {
1995             if(en1==0 && en2==0) {fNormQPairSwitch_E0E0[i]->AddAt('1',j);}
1996             if(en1==0 && en2==1) {fNormQPairSwitch_E0E1[i]->AddAt('1',j);}
1997             if(en1==0 && en2==2) {fNormQPairSwitch_E0E2[i]->AddAt('1',j);}
1998             if(en1==0 && en2==3) {fNormQPairSwitch_E0E3[i]->AddAt('1',j);}
1999             if(en1==1 && en2==1) {fNormQPairSwitch_E1E1[i]->AddAt('1',j);}
2000             if(en1==1 && en2==2) {fNormQPairSwitch_E1E2[i]->AddAt('1',j);}
2001             if(en1==1 && en2==3) {fNormQPairSwitch_E1E3[i]->AddAt('1',j);}
2002             if(en1==2 && en2==3) {fNormQPairSwitch_E2E3[i]->AddAt('1',j);}
2003           }
2004           
2005         }
2006       }
2007     }
2008   }
2009     
2010   //cout<<PairCount[0]<<"  "<<PairCount[1]<<"  "<<PairCount[2]<<"  "<<PairCount[3]<<"  "<<PairCount[4]<<"  "<<PairCount[5]<<"  "<<PairCount[6]<<endl;
2011   //cout<<NormPairCount[0]<<"  "<<NormPairCount[1]<<"  "<<NormPairCount[2]<<"  "<<NormPairCount[3]<<"  "<<NormPairCount[4]<<"  "<<NormPairCount[5]<<"  "<<NormPairCount[6]<<endl;
2012   ///////////////////////////////////////////////////  
2013   // Do not use pairs from events with too many pairs
2014   
2015   ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
2016   
2017   ///////////////////////////////////////////////////
2018   
2019   
2020   if(fTabulatePairs) return;
2021
2022     
2023   
2024   ////////////////////////////////////////////////////
2025   ////////////////////////////////////////////////////
2026   // Normalization counting of 3- and 4-particle terms
2027   for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2028     for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2029       if(en2==0 && en3>2) continue;// not needed config
2030       if(en2==1 && en3==en2) continue;// not needed config
2031       for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2032         if(en3==0 && en4>1) continue;// not needed config
2033         if(en3==1 && en4==3) continue;// not needed configs
2034         if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2035         
2036         for (Int_t i=0; i<myTracks; i++) {// 1st particle
2037           pVect1[1]=(fEvt)->fTracks[i].fP[0];
2038           pVect1[2]=(fEvt)->fTracks[i].fP[1];
2039           pVect1[3]=(fEvt)->fTracks[i].fP[2];
2040           ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2041           
2042           for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2043             if(en2==0) {if(fNormQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2044             else {if(fNormQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2045             
2046             pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2047             pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2048             pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2049             ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2050             
2051             for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2052               if(en3==0) {
2053                 if(fNormQPairSwitch_E0E0[i]->At(k)=='0') continue;
2054                 if(fNormQPairSwitch_E0E0[j]->At(k)=='0') continue;
2055               }else if(en3==1){
2056                 if(fNormQPairSwitch_E0E1[i]->At(k)=='0') continue;
2057                 if(fNormQPairSwitch_E0E1[j]->At(k)=='0') continue;
2058               }else{
2059                 if(fNormQPairSwitch_E0E2[i]->At(k)=='0') continue;
2060                 if(fNormQPairSwitch_E1E2[j]->At(k)=='0') continue;
2061               }
2062               
2063               pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2064               pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2065               pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2066               ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2067               Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2068               SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2069               
2070               Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2071               if(KT3<=fKT3transition) KT3index=0;
2072               else KT3index=1;
2073               
2074               if(en2==0 && en3==0 && en4==0) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fNorm3->Fill(0);
2075               if(en2==1 && en3==2 && en4==3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fNorm3->Fill(0);
2076               if(en2==0 && en3==1 && en4==2) {
2077                 if(fill2) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fNorm3->Fill(0);
2078                 if(fill3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fNorm3->Fill(0);
2079                 if(fill4) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fNorm3->Fill(0);
2080               }
2081               
2082               
2083               for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2084                 if(en4==0){
2085                   if(fNormQPairSwitch_E0E0[i]->At(l)=='0') continue;
2086                   if(fNormQPairSwitch_E0E0[j]->At(l)=='0') continue;
2087                   if(fNormQPairSwitch_E0E0[k]->At(l)=='0') continue;
2088                 }else if(en4==1){
2089                   if(en3==0){
2090                     if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2091                     if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2092                     if(fNormQPairSwitch_E0E1[k]->At(l)=='0') continue;
2093                   }else{
2094                     if(fNormQPairSwitch_E0E1[i]->At(l)=='0') continue;
2095                     if(fNormQPairSwitch_E0E1[j]->At(l)=='0') continue;
2096                     if(fNormQPairSwitch_E1E1[k]->At(l)=='0') continue;
2097                   }
2098                 }else if(en4==2){
2099                   if(fNormQPairSwitch_E0E2[i]->At(l)=='0') continue;
2100                   if(fNormQPairSwitch_E0E2[j]->At(l)=='0') continue;
2101                   if(fNormQPairSwitch_E1E2[k]->At(l)=='0') continue;
2102                 }else{
2103                   if(fNormQPairSwitch_E0E3[i]->At(l)=='0') continue;
2104                   if(fNormQPairSwitch_E1E3[j]->At(l)=='0') continue;
2105                   if(fNormQPairSwitch_E2E3[k]->At(l)=='0') continue;
2106                 }
2107
2108                 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2109                 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2110                 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2111                 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2112                 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.;
2113                 if(KT4<=fKT4transition) KT4index=0;
2114                 else KT4index=1;
2115                 
2116                 Bool_t FillTerms[13]={kFALSE};
2117                 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, en2+en3+en4, FillTerms);
2118                 //
2119                 for(int ft=0; ft<13; ft++) {
2120                   if(FillTerms[ft]) Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fNorm4->Fill(0.); 
2121                 }
2122                 
2123                 
2124               }
2125             }
2126           }
2127         }  
2128         
2129       }
2130     }
2131   }
2132     
2133
2134    
2135
2136     ///////////////////////////////////////////////////////////////////////
2137     ///////////////////////////////////////////////////////////////////////
2138     ///////////////////////////////////////////////////////////////////////
2139     //
2140     //
2141     // Start the Main Correlation Analysis
2142     //
2143     //
2144     ///////////////////////////////////////////////////////////////////////
2145   
2146
2147
2148     ////////////////////////////////////////////////////
2149     ////////////////////////////////////////////////////
2150     for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
2151       for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
2152         if(en2==0 && en3>2) continue;// not needed config
2153         if(en2==1 && en3==en2) continue;// not needed config
2154         for(Int_t en4=en3; en4<=3; en4++){// 4th event number
2155           if(en3==0 && en4>1) continue;// not needed config
2156           if(en3==1 && en4==3) continue;// not needed configs
2157           if(en3==2 && (en2+en3+en4)!=6) continue;// not needed configs
2158           
2159           Int_t ENsum=en2+en3+en4;// 0 or 1 or 3 or 6
2160           
2161           /////////////////////////////////////////////////////////////
2162           for (Int_t i=0; i<myTracks; i++) {// 1st particle
2163             pVect1[0]=(fEvt)->fTracks[i].fEaccepted;
2164             pVect1[1]=(fEvt)->fTracks[i].fP[0];
2165             pVect1[2]=(fEvt)->fTracks[i].fP[1];
2166             pVect1[3]=(fEvt)->fTracks[i].fP[2];
2167             ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
2168
2169             /////////////////////////////////////////////////////////////
2170             for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
2171               if(en2==0) {if(fLowQPairSwitch_E0E0[i]->At(j)=='0') continue;}
2172               else {if(fLowQPairSwitch_E0E1[i]->At(j)=='0') continue;}
2173
2174               pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
2175               pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
2176               pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
2177               pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
2178               ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
2179               qinv12 = GetQinv(pVect1, pVect2);
2180               kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
2181               SetFillBins2(ch1, ch2, bin1, bin2);
2182               Int_t kTindex=0;
2183               if(kT12<=0.3) kTindex=0;
2184               else kTindex=1;
2185               
2186               FSICorr12 = FSICorrelation(ch1,ch2, qinv12);
2187               
2188               // two particle terms filled during tabulation of low-q pairs
2189               
2190               
2191               if(fMCcase){
2192                 FilledMCpair12=kFALSE;
2193
2194                 if(ch1==ch2 && fMbin==0 && qinv12<0.2 && ENsum!=2 && ENsum!=3 && ENsum!=6){
2195                   for(Int_t rstep=0; rstep<10; rstep++){
2196                     Float_t coeff = (rstep)*0.2*(0.18/1.2);
2197                     Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
2198                     if(phi1 > 2*PI) phi1 -= 2*PI;
2199                     if(phi1 < 0) phi1 += 2*PI;
2200                     Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
2201                     if(phi2 > 2*PI) phi2 -= 2*PI;
2202                     if(phi2 < 0) phi2 += 2*PI;
2203                     Float_t deltaphi = phi1 - phi2;
2204                     if(deltaphi > PI) deltaphi -= PI;
2205                     if(deltaphi < -PI) deltaphi += PI;
2206                     
2207                     if(ENsum==0) ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2208                     else ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
2209                   }
2210                   
2211                 }// pair selection
2212
2213                 // Check that label does not exceed stack size
2214                 if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
2215                   if(ENsum==0 && abs((fEvt+en2)->fTracks[j].fLabel) == abs((fEvt)->fTracks[i].fLabel)) continue;
2216                   pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
2217                   pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
2218                   pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
2219                   pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
2220                   pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
2221                   qinv12MC = GetQinv(pVect1MC, pVect2MC);
2222                   
2223                   if(qinv12<0.1 && ch1==ch2 && ENsum==0) {
2224                     ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC); 
2225                     ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
2226                     ((TH2D*)fOutputList->FindObject("fQ2Res"))->Fill(kT12, qinv12-qinv12MC);
2227                   }
2228                                   
2229                   // secondary contamination
2230                   if(ENsum==0){
2231                     mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
2232                     mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
2233                     if(!mcParticle1 || !mcParticle2) continue;
2234                     if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
2235                       if(ch1==ch2) {
2236                         ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2237                         if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2238                           ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2239                         }             
2240                       }else{
2241                         ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2242                         if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
2243                           ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
2244                         }
2245                       }
2246                     }
2247                   }
2248                   
2249                   if(ENsum==6){// all mixed events
2250                     Int_t chGroup2[2]={ch1,ch2};
2251
2252                     Float_t rForQW=5.0;
2253                     if(fFSIindex<=1) rForQW=10;
2254                     else if(fFSIindex==2) rForQW=9;
2255                     else if(fFSIindex==3) rForQW=8;
2256                     else if(fFSIindex==4) rForQW=7;
2257                     else if(fFSIindex==5) rForQW=6;
2258                     else if(fFSIindex==6) rForQW=5;
2259                     else if(fFSIindex==7) rForQW=4;
2260                     else if(fFSIindex==8) rForQW=3;
2261                     else rForQW=2;
2262                     
2263                     
2264                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinv->Fill(qinv12MC, MCWeight(chGroup2, rForQW, 0.65, qinv12MC, 0.));// was 4,5
2265                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(chGroup2, rForQW, 0.65, qinv12MC, 0.));// was 4,5
2266                     // pion purity
2267                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityDen->Fill(kT12, qinv12);
2268                     Int_t SCNumber = 1;
2269                     Int_t PdgCodeSum = abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode) + abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode);
2270                     if(PdgCodeSum==22) SCNumber=1;// e-e
2271                     else if(PdgCodeSum==24) SCNumber=2;// e-mu
2272                     else if(PdgCodeSum==222) SCNumber=3;// e-pi
2273                     else if(PdgCodeSum==332) SCNumber=4;// e-k
2274                     else if(PdgCodeSum==2223) SCNumber=5;// e-p
2275                     else if(PdgCodeSum==26) SCNumber=6;// mu-mu
2276                     else if(PdgCodeSum==224) SCNumber=7;// mu-pi
2277                     else if(PdgCodeSum==334) SCNumber=8;// mu-k
2278                     else if(PdgCodeSum==2225) SCNumber=9;// mu-p
2279                     else if(PdgCodeSum==422) SCNumber=10;// pi-pi
2280                     else if(PdgCodeSum==532) SCNumber=11;// pi-k
2281                     else if(PdgCodeSum==2423) SCNumber=12;// pi-p
2282                     else if(PdgCodeSum==642) SCNumber=13;// k-k
2283                     else if(PdgCodeSum==2533) SCNumber=14;// k-p
2284                     else if(PdgCodeSum==4424) SCNumber=15;// p-p
2285                     else {SCNumber=16;}
2286                     
2287                     Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityNum->Fill(SCNumber, kT12, qinv12);
2288                     
2289                     ///////////////////////
2290                     // muon contamination
2291                     Pparent1[0]=pVect1MC[0]; Pparent1[1]=pVect1MC[1]; Pparent1[2]=pVect1MC[2]; Pparent1[3]=pVect1MC[3];
2292                     Pparent2[0]=pVect2MC[0]; Pparent2[1]=pVect2MC[1]; Pparent2[2]=pVect2MC[2]; Pparent2[3]=pVect2MC[3];
2293                     pionParent1=kFALSE; pionParent2=kFALSE;
2294                     FilledMCpair12=kTRUE;
2295                     //
2296                     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
2297                       Int_t MotherLabel1 = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fMotherLabel;
2298                       if(abs((fEvt)->fMCtracks[MotherLabel1].fPdgCode)==211) {
2299                         pionParent1=kTRUE;
2300                         Pparent1[1] = (fEvt)->fMCtracks[MotherLabel1].fPx; Pparent1[2] = (fEvt)->fMCtracks[MotherLabel1].fPy; Pparent1[3] = (fEvt)->fMCtracks[MotherLabel1].fPz;
2301                         Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
2302                       }
2303                       // 
2304                       if(abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13) {
2305                         Int_t MotherLabel2 = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fMotherLabel;
2306                         if(abs((fEvt+en2)->fMCtracks[MotherLabel2].fPdgCode)==211) {
2307                           pionParent2=kTRUE;
2308                           Pparent2[1] = (fEvt+en2)->fMCtracks[MotherLabel2].fPx; Pparent2[2] = (fEvt+en2)->fMCtracks[MotherLabel2].fPy; Pparent2[3] = (fEvt+en2)->fMCtracks[MotherLabel2].fPz;
2309                           Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
2310                         }
2311                       }
2312                       
2313                       parentQinv12 = GetQinv(Pparent1, Pparent2);
2314                       
2315                       if(pionParent1 || pionParent2){
2316                         if(parentQinv12 > 0.001 && parentQinv12 < 0.3){
2317                           Float_t muonPionK12 = FSICorrelation(ch1, ch2, qinv12MC);
2318                           Float_t pionPionK12 = FSICorrelation(ch1, ch2, parentQinv12);
2319                           for(Int_t term=1; term<=2; term++){
2320                             for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2321                               Float_t Rvalue = 5+Riter;
2322                               Float_t WInput = 1.0;
2323                               if(term==1) {
2324                                 WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.);
2325                               }else{
2326                                 muonPionK12 = 1.0; pionPionK12=1.0;
2327                               }
2328                               
2329                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput);
2330                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput);
2331                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12);
2332                               Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fPionPionK2->Fill(Rvalue, parentQinv12, pionPionK12);
2333                             }// Riter
2334                           }// term loop
2335                           
2336                           if(ch1==ch2) ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(0., kT12, qinv12MC-parentQinv12);
2337                           else ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(1., kT12, qinv12MC-parentQinv12);
2338                         }// parentQ check
2339                       }// pion parent check
2340                     }// muon check
2341                   
2342                     
2343                     // momentum resolution
2344                     for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2345                       Float_t Rvalue = 5+Riter;
2346                       Float_t WInput = MCWeight(chGroup2, Rvalue, 0.65, qinv12MC, 0.);
2347                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fIdeal->Fill(Rvalue, qinv12MC, WInput);
2348                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fIdeal->Fill(Rvalue, qinv12MC);
2349                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fSmeared->Fill(Rvalue, qinv12, WInput);
2350                       Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fSmeared->Fill(Rvalue, qinv12);
2351                     }
2352                     
2353                   }// ENsum check
2354                 }// MC array check
2355               }// MC case
2356               
2357                 
2358
2359               /////////////////////////////////////////////////////////////
2360               for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
2361                 if(en3==0) {
2362                   if(fLowQPairSwitch_E0E0[i]->At(k)=='0') continue;
2363                   if(fLowQPairSwitch_E0E0[j]->At(k)=='0') continue;
2364                 }else if(en3==1){
2365                   if(fLowQPairSwitch_E0E1[i]->At(k)=='0') continue;
2366                   if(fLowQPairSwitch_E0E1[j]->At(k)=='0') continue;
2367                 }else{
2368                   if(fLowQPairSwitch_E0E2[i]->At(k)=='0') continue;
2369                   if(fLowQPairSwitch_E1E2[j]->At(k)=='0') continue;
2370                 }
2371                 
2372                 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
2373                 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
2374                 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
2375                 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
2376                 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
2377                 qinv13 = GetQinv(pVect1, pVect3);
2378                 qinv23 = GetQinv(pVect2, pVect3);
2379                 q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
2380                 
2381                 FilledMCtriplet123 = kFALSE;
2382                 
2383                 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
2384                 SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
2385                 
2386                 Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
2387                 if(KT3<=fKT3transition) KT3index=0;
2388                 else KT3index=1;
2389                 
2390                 FSICorr13 = FSICorrelation(ch1,ch3, qinv13);
2391                 FSICorr23 = FSICorrelation(ch2,ch3, qinv23);
2392
2393                 if(ENsum==0) {
2394                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fTerms3->Fill(q3); 
2395                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactor->Fill(q3, 1/(FSICorr12*FSICorr13*FSICorr23));
2396                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv12);
2397                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv13);
2398                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv23);
2399                 }
2400                 if(ENsum==6) {
2401                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTerms3->Fill(q3);
2402                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv12);
2403                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv13);
2404                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv23);
2405                 }
2406                 if(ENsum==3){
2407                   if(fill2) {
2408                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fTerms3->Fill(q3);
2409                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactor->Fill(q3, 1/(FSICorr12));
2410                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv12);
2411                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv13);
2412                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv23);
2413                   }if(fill3) {
2414                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fTerms3->Fill(q3);
2415                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactor->Fill(q3, 1/(FSICorr12));
2416                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv12);
2417                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv13);
2418                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv23);
2419                   }if(fill4) {
2420                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fTerms3->Fill(q3);
2421                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactor->Fill(q3, 1/(FSICorr12));
2422                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv12);
2423                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv13);
2424                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv23);
2425                   }
2426                 }
2427                 
2428                 // r3 denominator
2429                 if(ENsum==6 && ch1==ch2 && ch1==ch3){
2430                   GoodTripletWeights = kFALSE;
2431                   //
2432                   GetWeight(pVect1, pVect2, weight12, weight12Err);
2433                   GetWeight(pVect1, pVect3, weight13, weight13Err);
2434                   GetWeight(pVect2, pVect3, weight23, weight23Err);
2435                   
2436                   if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {// weight should never be larger than 1
2437                     if(fMbin==0 && bin1==0) {
2438                       ((TH1D*)fOutputList->FindObject("fTPNRejects3pion1"))->Fill(q3, sqrt(fabs(weight12*weight13*weight23)));
2439                     }
2440                   }else{
2441                     
2442                     Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
2443                     Float_t MuonCorr12=1.0, MuonCorr13=1.0, MuonCorr23=1.0;
2444                     if(!fGenerateSignal && !fMCcase) {
2445                       Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
2446                       Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
2447                       Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);            
2448                       if(momBin12 >= 20) momBin12 = 19;
2449                       if(momBin13 >= 20) momBin13 = 19;
2450                       if(momBin23 >= 20) momBin23 = 19;
2451                       MomResCorr12 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin12);
2452                       MomResCorr13 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin13);
2453                       MomResCorr23 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin23);
2454                       MuonCorr12 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin12);
2455                       MuonCorr13 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin13);
2456                       MuonCorr23 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin23);
2457                     }
2458                     // no MRC, no Muon Correction
2459                     weight12CC[0] = ((weight12+1) - ffcSq*FSICorr12 - (1-ffcSq));
2460                     weight12CC[0] /= FSICorr12*ffcSq;
2461                     weight13CC[0] = ((weight13+1) - ffcSq*FSICorr13 - (1-ffcSq));
2462                     weight13CC[0] /= FSICorr13*ffcSq;
2463                     weight23CC[0] = ((weight23+1) - ffcSq*FSICorr23 - (1-ffcSq));
2464                     weight23CC[0] /= FSICorr23*ffcSq;
2465                     if(weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2466                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, sqrt(weight12CC[0]*weight13CC[0]*weight23CC[0]));
2467                     }
2468                     // no Muon Correction
2469                     weight12CC[1] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2470                     weight12CC[1] /= FSICorr12*ffcSq;
2471                     weight13CC[1] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2472                     weight13CC[1] /= FSICorr13*ffcSq;
2473                     weight23CC[1] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2474                     weight23CC[1] /= FSICorr23*ffcSq;
2475                     if(weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2476                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(2, q3, sqrt(weight12CC[1]*weight13CC[1]*weight23CC[1]));
2477                     }
2478                     // both Corrections
2479                     weight12CC[2] = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
2480                     weight12CC[2] /= FSICorr12*ffcSq;
2481                     weight12CC[2] *= MuonCorr12;
2482                     weight13CC[2] = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
2483                     weight13CC[2] /= FSICorr13*ffcSq;
2484                     weight13CC[2] *= MuonCorr13;
2485                     weight23CC[2] = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
2486                     weight23CC[2] /= FSICorr23*ffcSq;
2487                     weight23CC[2] *= MuonCorr23;
2488                     if(weight12CC[2] < 0 || weight13CC[2] < 0 || weight23CC[2] < 0) {// C2^QS can never be less than unity
2489                       if(fMbin==0 && bin1==0) {
2490                         ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC[2]*weight13CC[2]*weight23CC[2])));
2491                       }
2492                     }else{
2493                       GoodTripletWeights = kTRUE;
2494                       /////////////////////////////////////////////////////
2495                       weightTotal = sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]);
2496                       /////////////////////////////////////////////////////
2497                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(3, q3, weightTotal);
2498                       //
2499                       // Full Weight reconstruction
2500                       weightTotal *= 2.0;
2501                       weightTotal += weight12CC[2] + weight13CC[2] + weight23CC[2];
2502                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(4, q3, weightTotal);
2503                     }// 2nd r3 den check
2504                     
2505                    
2506                   }// 1st r3 den check
2507                   
2508                 }// r3 den
2509                 
2510
2511                 if(ch1==ch2 && ch1==ch3 && ENsum==0){
2512                   ((TH3D*)fOutputList->FindObject("fKT3DistTerm1"))->Fill(fMbin+1, KT3, q3);
2513                   if(q3<0.06){
2514                     Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2515                     Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2516                     Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2517                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt1);
2518                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt2);
2519                     ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt3);
2520                   }
2521                 }
2522                 if(ch1==ch2 && ch1==ch3 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT3DistTerm5"))->Fill(fMbin+1, KT3, q3);
2523                 
2524                 
2525
2526                 
2527                 if(fMCcase && ENsum==6 && FilledMCpair12){// for momentum resolution and muon correction
2528                   if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
2529                     
2530                     pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
2531                     pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
2532                     pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
2533                     pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
2534                     qinv13MC = GetQinv(pVect1MC, pVect3MC);
2535                     qinv23MC = GetQinv(pVect2MC, pVect3MC);
2536                     
2537                     q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
2538                     if(q3<0.1 && ch1==ch2 && ch1==ch3) ((TH2D*)fOutputList->FindObject("fQ3Res"))->Fill(KT3, q3-q3MC);
2539                     
2540                     Int_t chGroup3[3]={ch1,ch2,ch3};
2541                     Float_t QinvMCGroup3[3]={qinv12MC,qinv13MC,qinv23MC};
2542                     //Float_t kTGroup3[3]={float(sqrt(pow(pVect1MC[1]+pVect2MC[1],2) + pow(pVect1MC[2]+pVect2MC[2],2))/2.),
2543                     //float(sqrt(pow(pVect1MC[1]+pVect3MC[1],2) + pow(pVect1MC[2]+pVect3MC[2],2))/2.),
2544                     //float(sqrt(pow(pVect2MC[1]+pVect3MC[1],2) + pow(pVect2MC[2]+pVect3MC[2],2))/2.)};
2545                     Float_t kTGroup3[3]={0};
2546                     
2547                     Pparent3[0]=pVect3MC[0]; Pparent3[1]=pVect3MC[1]; Pparent3[2]=pVect3MC[2]; Pparent3[3]=pVect3MC[3];
2548                     pionParent3=kFALSE;
2549                   
2550                     if(abs((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPdgCode)==13){// muon check
2551                       Int_t MotherLabel3 = (fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fMotherLabel;
2552                       if(abs((fEvt+en3)->fMCtracks[MotherLabel3].fPdgCode)==211) {
2553                         pionParent3=kTRUE;
2554                         Pparent3[1] = (fEvt+en3)->fMCtracks[MotherLabel3].fPx; Pparent3[2] = (fEvt+en3)->fMCtracks[MotherLabel3].fPy; Pparent3[3] = (fEvt+en3)->fMCtracks[MotherLabel3].fPz;
2555                         Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
2556                       }
2557                     }
2558                     
2559                     parentQinv13 = GetQinv(Pparent1, Pparent3);
2560                     parentQinv23 = GetQinv(Pparent2, Pparent3);
2561                     parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
2562                    
2563                     if(parentQinv12 > 0.001 && parentQinv13 > 0.001 && parentQinv23 > 0.001 && parentQ3 < 0.5){
2564                       FilledMCtriplet123=kTRUE;
2565                       if(pionParent1 || pionParent2 || pionParent3) {// want at least one pion-->muon
2566                         
2567                         Float_t parentQinvGroup3[3]={parentQinv12, parentQinv13, parentQinv23};
2568                         //Float_t parentkTGroup3[3]={float(sqrt(pow(Pparent1[1]+Pparent2[1],2) + pow(Pparent1[2]+Pparent2[2],2))/2.),
2569                         //float(sqrt(pow(Pparent1[1]+Pparent3[1],2) + pow(Pparent1[2]+Pparent3[2],2))/2.),
2570                         //float(sqrt(pow(Pparent2[1]+Pparent3[1],2) + pow(Pparent2[2]+Pparent3[2],2))/2.)};
2571                         Float_t parentkTGroup3[3]={0};
2572                         
2573                         ((TH2D*)fOutputList->FindObject("fAvgQ12VersusQ3"))->Fill(parentQ3, parentQinv12);
2574                         ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13);
2575                         ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23);
2576
2577                         for(Int_t term=1; term<=4; term++){
2578                           if(term==1) {}
2579                           else if(term==2) {if(!pionParent1 && !pionParent2) continue;}
2580                           else if(term==3) {if(!pionParent1 && !pionParent3) continue;}
2581                           else {if(!pionParent2 && !pionParent3) continue;}
2582                           for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2583                             Float_t Rvalue = 5+Riter;
2584                             Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
2585                             Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
2586                             Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
2587                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(1, Rvalue, q3MC, WInput);
2588                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ3, WInput);
2589                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(1, Rvalue, q3MC, WInputFSI);
2590                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(1, Rvalue, parentQ3, WInputParentFSI);
2591                             //
2592                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(2, Rvalue, q3MC);
2593                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ3);
2594                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(2, Rvalue, q3MC);
2595                             Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(2, Rvalue, parentQ3);
2596                           }// Riter
2597                         }// term loop
2598                     
2599                       }// pion parent check
2600                     }// parentQ check (muon correction)
2601                     
2602                     // 3-pion momentum resolution
2603                     for(Int_t term=1; term<=5; term++){
2604                       for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2605                         Float_t Rvalue = 5+Riter;
2606                         Float_t WInput = MCWeight3(term, Rvalue, 0.65, chGroup3, QinvMCGroup3, kTGroup3);
2607                         Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fIdeal->Fill(Rvalue, q3MC, WInput);
2608                         Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fSmeared->Fill(Rvalue, q3, WInput);
2609                       }
2610                     }
2611                     
2612                   }// 3rd particle label check
2613                 }// MCcase and ENsum==6
2614                 
2615                 
2616                 
2617
2618                 /////////////////////////////////////////////////////////////
2619                 for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
2620                   if(en4==0){
2621                     if(fLowQPairSwitch_E0E0[i]->At(l)=='0') continue;
2622                     if(fLowQPairSwitch_E0E0[j]->At(l)=='0') continue;
2623                     if(fLowQPairSwitch_E0E0[k]->At(l)=='0') continue;
2624                   }else if(en4==1){
2625                     if(en3==0){
2626                       if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2627                       if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2628                       if(fLowQPairSwitch_E0E1[k]->At(l)=='0') continue;
2629                     }else{ 
2630                       if(fLowQPairSwitch_E0E1[i]->At(l)=='0') continue;
2631                       if(fLowQPairSwitch_E0E1[j]->At(l)=='0') continue;
2632                       if(fLowQPairSwitch_E1E1[k]->At(l)=='0') continue;
2633                     }
2634                   }else if(en4==2){
2635                     if(fLowQPairSwitch_E0E2[i]->At(l)=='0') continue;
2636                     if(fLowQPairSwitch_E0E2[j]->At(l)=='0') continue;
2637                     if(fLowQPairSwitch_E1E2[k]->At(l)=='0') continue;
2638                   }else{
2639                     if(fLowQPairSwitch_E0E3[i]->At(l)=='0') continue;
2640                     if(fLowQPairSwitch_E1E3[j]->At(l)=='0') continue;
2641                     if(fLowQPairSwitch_E2E3[k]->At(l)=='0') continue;
2642                   }
2643
2644                   pVect4[0]=(fEvt+en4)->fTracks[l].fEaccepted;
2645                   pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
2646                   pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
2647                   pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
2648                   ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
2649                   qinv14 = GetQinv(pVect1, pVect4);
2650                   qinv24 = GetQinv(pVect2, pVect4);
2651                   qinv34 = GetQinv(pVect3, pVect4);
2652                   q4 = sqrt(pow(q3,2) + pow(qinv14,2) + pow(qinv24,2) + pow(qinv34,2));
2653                   
2654                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
2655                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(1, qinv12); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(2, qinv13); 
2656                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(3, qinv14); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(4, qinv23); 
2657                     ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(5, qinv24); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(6, qinv34);
2658                   }
2659                   
2660                   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.;
2661                   if(KT4<=fKT4transition) KT4index=0;
2662                   else KT4index=1;
2663                   
2664                   FSICorr14 = FSICorrelation(ch1,ch4, qinv14);
2665                   FSICorr24 = FSICorrelation(ch2,ch4, qinv24);
2666                   FSICorr34 = FSICorrelation(ch3,ch4, qinv34);
2667
2668                   Bool_t FillTerms[13]={kFALSE};
2669                   SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, ENsum, FillTerms);
2670                   //
2671                   for(int ft=0; ft<13; ft++) {
2672                     Float_t FSIfactor = 1.0;
2673                     if(ft==0) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr14 * FSICorr23 * FSICorr24 * FSICorr34);
2674                     else if(ft<=4) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr23);
2675                     else if(ft<=10) FSIfactor = 1/(FSICorr12);
2676                     else if(ft==11) FSIfactor = 1/(FSICorr12 * FSICorr34);
2677                     else FSIfactor = 1.0;
2678                     if(FillTerms[ft]) {
2679                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fTerms4->Fill(q4);
2680                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactor->Fill(q4, FSIfactor);
2681                     }
2682                   }
2683                   
2684                   /////////////////////////////////////////////////////////////
2685                   // r4{2}
2686                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6 && GoodTripletWeights){
2687                     GetWeight(pVect1, pVect4, weight14, weight14Err);
2688                     GetWeight(pVect2, pVect4, weight24, weight24Err);
2689                     GetWeight(pVect3, pVect4, weight34, weight34Err);
2690                     
2691                     Float_t MomResCorr14=1.0, MomResCorr24=1.0, MomResCorr34=1.0;
2692                     Float_t MuonCorr14=1.0, MuonCorr24=1.0, MuonCorr34=1.0;
2693                     if(!fGenerateSignal && !fMCcase) {
2694                       Int_t momBin14 = fMomResC2->GetYaxis()->FindBin(qinv14);
2695                       Int_t momBin24 = fMomResC2->GetYaxis()->FindBin(qinv24);
2696                       Int_t momBin34 = fMomResC2->GetYaxis()->FindBin(qinv34);            
2697                       if(momBin14 >= 20) momBin14 = 19;
2698                       if(momBin24 >= 20) momBin24 = 19;
2699                       if(momBin34 >= 20) momBin34 = 19;
2700                       MomResCorr14 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin14);
2701                       MomResCorr24 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin24);
2702                       MomResCorr34 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin34);
2703                       MuonCorr14 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin14);
2704                       MuonCorr24 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin24);
2705                       MuonCorr34 = fWeightmuonCorrection->GetBinContent(rBinForTPNMomRes, momBin34);
2706                     }
2707                     
2708                     // no MRC, no Muon Correction
2709                     weight14CC[0] = ((weight14+1) - ffcSq*FSICorr14 - (1-ffcSq));
2710                     weight14CC[0] /= FSICorr14*ffcSq;
2711                     weight24CC[0] = ((weight24+1) - ffcSq*FSICorr24 - (1-ffcSq));
2712                     weight24CC[0] /= FSICorr24*ffcSq;
2713                     weight34CC[0] = ((weight34+1) - ffcSq*FSICorr34 - (1-ffcSq));
2714                     weight34CC[0] /= FSICorr34*ffcSq;
2715                     if(weight14CC[0] > 0 && weight24CC[0] > 0 && weight34CC[0] > 0 && weight12CC[0] > 0 && weight13CC[0] > 0 && weight23CC[0] > 0){
2716                       weightTotal  = sqrt(weight12CC[0]*weight13CC[0]*weight24CC[0]*weight34CC[0]);
2717                       weightTotal += sqrt(weight12CC[0]*weight14CC[0]*weight23CC[0]*weight34CC[0]);
2718                       weightTotal += sqrt(weight13CC[0]*weight14CC[0]*weight23CC[0]*weight24CC[0]);
2719                       weightTotal /= 3.;
2720                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(1, q4, weightTotal);
2721                     }
2722                     // no Muon Correction
2723                     weight14CC[1] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2724                     weight14CC[1] /= FSICorr14*ffcSq;
2725                     weight24CC[1] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2726                     weight24CC[1] /= FSICorr24*ffcSq;
2727                     weight34CC[1] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2728                     weight34CC[1] /= FSICorr34*ffcSq;
2729                     if(weight14CC[1] > 0 && weight24CC[1] > 0 && weight34CC[1] > 0 && weight12CC[1] > 0 && weight13CC[1] > 0 && weight23CC[1] > 0){
2730                       weightTotal  = sqrt(weight12CC[1]*weight13CC[1]*weight24CC[1]*weight34CC[1]);
2731                       weightTotal += sqrt(weight12CC[1]*weight14CC[1]*weight23CC[1]*weight34CC[1]);
2732                       weightTotal += sqrt(weight13CC[1]*weight14CC[1]*weight23CC[1]*weight24CC[1]);
2733                       weightTotal /= 3.;
2734                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(2, q4, weightTotal);
2735                     }
2736                     // both corrections
2737                     weight14CC[2] = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
2738                     weight14CC[2] /= FSICorr14*ffcSq;
2739                     weight14CC[2] *= MuonCorr14;
2740                     weight24CC[2] = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
2741                     weight24CC[2] /= FSICorr24*ffcSq;
2742                     weight24CC[2] *= MuonCorr24;
2743                     weight34CC[2] = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
2744                     weight34CC[2] /= FSICorr34*ffcSq;
2745                     weight34CC[2] *= MuonCorr34;
2746  
2747                     if(weight14CC[2] < 0 || weight24CC[2] < 0 || weight34CC[2] < 0) {// C2^QS can never be less than unity
2748                       if(fMbin==0 && bin1==0 && KT4index==0) {
2749                         ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC[2]*weight23CC[2]*weight34CC[2]*weight14CC[2])));
2750                       }
2751                     }else{
2752                       /////////////////////////////////////////////////////
2753                       weightTotal  = sqrt(weight12CC[2]*weight13CC[2]*weight24CC[2]*weight34CC[2]);
2754                       weightTotal += sqrt(weight12CC[2]*weight14CC[2]*weight23CC[2]*weight34CC[2]);
2755                       weightTotal += sqrt(weight13CC[2]*weight14CC[2]*weight23CC[2]*weight24CC[2]);
2756                       weightTotal /= 3.;
2757                       /////////////////////////////////////////////////////
2758                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(3, q4, weightTotal);
2759                       // Full Weight reconstruction
2760                       weightTotal *= 6.0;
2761                       weightTotal += weight12CC[2] + weight13CC[2] + weight14CC[2] + weight23CC[2] + weight24CC[2] + weight34CC[2];
2762                       weightTotal += 2*sqrt(weight12CC[2]*weight13CC[2]*weight23CC[2]);
2763                       weightTotal += 2*sqrt(weight12CC[2]*weight14CC[2]*weight24CC[2]);
2764                       weightTotal += 2*sqrt(weight13CC[2]*weight14CC[2]*weight34CC[2]);
2765                       weightTotal += 2*sqrt(weight23CC[2]*weight24CC[2]*weight34CC[2]);
2766                       weightTotal += weight12CC[2]*weight34CC[2] + weight13CC[2]*weight24CC[2] + weight14CC[2]*weight23CC[2];
2767                       Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[12].fTwoPartNorm->Fill(4, q4, weightTotal);
2768                     }
2769                   }
2770                   /////////////////////////////////////////////////////////////
2771                   
2772                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==0){
2773                     ((TH3D*)fOutputList->FindObject("fKT4DistTerm1"))->Fill(fMbin+1, KT4, q4);
2774                     if(q4<0.085){
2775                       Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
2776                       Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
2777                       Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
2778                       Float_t pt4=sqrt(pow(pVect4[1],2)+pow(pVect4[2],2));
2779                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt1);
2780                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt2);
2781                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt3);
2782                       ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt4);
2783                     }
2784                   }
2785                   if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT4DistTerm13"))->Fill(fMbin+1, KT4, q4);
2786
2787
2788                   // momenumtum resolution and muon corrections
2789                   if(fMCcase && ENsum==6 && FilledMCtriplet123){// for momentum resolution and muon correction
2790                     if((fEvt+en4)->fTracks[l].fLabel < (fEvt+en4)->fMCarraySize){
2791                       
2792                       pVect4MC[0]=sqrt(pow((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
2793                       pVect4MC[1]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPx;
2794                       pVect4MC[2]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPy;
2795                       pVect4MC[3]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPz;
2796                       qinv14MC = GetQinv(pVect1MC, pVect4MC);
2797                       qinv24MC = GetQinv(pVect2MC, pVect4MC);
2798                       qinv34MC = GetQinv(pVect3MC, pVect4MC);
2799                       
2800                       q4MC = sqrt(pow(q3MC,2) + pow(qinv14MC,2) +  pow(qinv24MC,2) +  pow(qinv34MC,2));
2801                       if(q4<0.1 && ch1==ch2 && ch1==ch3 && ch1==ch4) ((TH2D*)fOutputList->FindObject("fQ4Res"))->Fill(KT4, q4-q4MC);
2802                       if(ch1==ch2 && ch1==ch3 && ch1==ch4) {
2803                         ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(1, qinv12MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(2, qinv13MC);
2804                         ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(3, qinv14MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(4, qinv23MC);
2805                         ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(5, qinv24MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(6, qinv34MC);
2806                       }
2807                       Int_t chGroup4[4]={ch1,ch2,ch3,ch4};
2808                       Float_t QinvMCGroup4[6]={qinv12MC, qinv13MC, qinv14MC, qinv23MC, qinv24MC, qinv34MC};
2809                       Float_t kTGroup4[6]={0};
2810                       
2811                       Pparent4[0]=pVect4MC[0]; Pparent4[1]=pVect4MC[1]; Pparent4[2]=pVect4MC[2]; Pparent4[3]=pVect4MC[3];
2812                       pionParent4=kFALSE;
2813                       if(abs((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPdgCode)==13){// muon check
2814                         Int_t MotherLabel4 = (fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fMotherLabel;
2815                         if(abs((fEvt+en4)->fMCtracks[MotherLabel4].fPdgCode)==211) {
2816                           pionParent4=kTRUE;
2817                           Pparent4[1] = (fEvt+en4)->fMCtracks[MotherLabel4].fPx; Pparent4[2] = (fEvt+en4)->fMCtracks[MotherLabel4].fPy; Pparent4[3] = (fEvt+en4)->fMCtracks[MotherLabel4].fPz;
2818                           Pparent4[0] = sqrt(pow(Pparent4[1],2)+pow(Pparent4[2],2)+pow(Pparent4[3],2)+pow(fTrueMassPi,2));
2819                         }
2820                       }
2821
2822                       parentQinv14 = GetQinv(Pparent1, Pparent4);
2823                       parentQinv24 = GetQinv(Pparent2, Pparent4);
2824                       parentQinv34 = GetQinv(Pparent3, Pparent4);
2825                       Float_t parentQ4 = sqrt(pow(parentQ3,2) + pow(parentQinv14,2) + pow(parentQinv24,2) + pow(parentQinv34,2));
2826                       
2827                       if(parentQinv14 > 0.001 && parentQinv24 > 0.001 && parentQinv34 > 0.001 && parentQ4 < 0.5){
2828                         if(pionParent1 || pionParent2 || pionParent3 || pionParent4) {// want at least one pion-->muon
2829                          
2830                           if(pionParent1) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(1);
2831                           if(pionParent2) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(2);
2832                           if(pionParent3) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(3);
2833                           if(pionParent4) ((TH1D*)fOutputList->FindObject("fDistPionParents4"))->Fill(4);
2834                           Float_t parentQinvGroup4[6]={parentQinv12, parentQinv13, parentQinv14, parentQinv23, parentQinv24, parentQinv34};
2835                           Float_t parentkTGroup4[6]={0};
2836                           
2837                           for(Int_t term=1; term<=12; term++){
2838                             if(term==1) {}
2839                             else if(term==2) {if(!pionParent1 && !pionParent2 && !pionParent3) continue;}
2840                             else if(term==3) {if(!pionParent1 && !pionParent2 && !pionParent4) continue;}
2841                             else if(term==4) {if(!pionParent1 && !pionParent3 && !pionParent4) continue;}
2842                             else if(term==5) {if(!pionParent2 && !pionParent3 && !pionParent4) continue;}
2843                             else if(term==6) {if(!pionParent1 && !pionParent2) continue;}
2844                             else if(term==7) {if(!pionParent1 && !pionParent3) continue;}
2845                             else if(term==8) {if(!pionParent1 && !pionParent4) continue;}
2846                             else if(term==9) {if(!pionParent2 && !pionParent3) continue;}
2847                             else if(term==10) {if(!pionParent2 && !pionParent4) continue;}
2848                             else if(term==11) {if(!pionParent3 && !pionParent4) continue;}
2849                             else {} 
2850                             for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2851                               Float_t Rvalue = 5+Riter;
2852                               Float_t WInput = MCWeight4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4, parentkTGroup4);
2853                               Float_t WInputParentFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4);
2854                               Float_t WInputFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, QinvMCGroup4);
2855                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(1, Rvalue, q4MC, WInput);
2856                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(1, Rvalue, parentQ4, WInput);
2857                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(1, Rvalue, q4MC, WInputFSI);
2858                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(1, Rvalue, parentQ4, WInputParentFSI);
2859                               //
2860                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(2, Rvalue, q4MC);
2861                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(2, Rvalue, parentQ4);
2862                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(2, Rvalue, q4MC);
2863                               Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(2, Rvalue, parentQ4);
2864                             }// Riter
2865                           }// term loop
2866                           
2867                         }// pion parent check
2868                       }// parentQ check (muon correction)
2869                       
2870                       // 4-pion momentum resolution
2871                       for(Int_t term=1; term<=13; term++){
2872                         for(Int_t Riter=0; Riter<fRVALUES; Riter++){
2873                           Float_t Rvalue = 5+Riter;
2874                           Float_t WInput = MCWeight4(term, Rvalue, 0.65, chGroup4, QinvMCGroup4, kTGroup4);
2875                           Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fIdeal->Fill(Rvalue, q4MC, WInput);
2876                           Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fSmeared->Fill(Rvalue, q4, WInput);
2877                         }
2878                       }
2879                       
2880                     }// label check particle 4
2881                   }// MCcase
2882                   
2883                 }// 4th particle
2884               }// 3rd particle
2885             }// 2nd particle
2886           }// 1st particle
2887           
2888         }// en4
2889       }// en3
2890     }// en2
2891     
2892     
2893
2894   
2895   
2896   
2897   
2898   // Post output data.
2899   PostData(1, fOutputList);
2900   
2901 }
2902 //________________________________________________________________________
2903 void AliFourPion::Terminate(Option_t *) 
2904 {
2905   // Called once at the end of the query
2906  
2907   cout<<"Done"<<endl;
2908
2909 }
2910 //________________________________________________________________________
2911 Bool_t AliFourPion::AcceptPair(AliFourPionTrackStruct first, AliFourPionTrackStruct second)
2912 {
2913   
2914   if(fabs(first.fEta-second.fEta) > fMinSepPairEta) return kTRUE;
2915   
2916   // propagate through B field to r=1m
2917   Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
2918   if(phi1 > 2*PI) phi1 -= 2*PI;
2919   if(phi1 < 0) phi1 += 2*PI;
2920   Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m 
2921   if(phi2 > 2*PI) phi2 -= 2*PI;
2922   if(phi2 < 0) phi2 += 2*PI;
2923   
2924   Float_t deltaphi = phi1 - phi2;
2925   if(deltaphi > PI) deltaphi -= 2*PI;
2926   if(deltaphi < -PI) deltaphi += 2*PI;
2927   deltaphi = fabs(deltaphi);
2928
2929   if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
2930     
2931   
2932   // propagate through B field to r=1.6m
2933   phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
2934   if(phi1 > 2*PI) phi1 -= 2*PI;
2935   if(phi1 < 0) phi1 += 2*PI;
2936   phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m 
2937   if(phi2 > 2*PI) phi2 -= 2*PI;
2938   if(phi2 < 0) phi2 += 2*PI;
2939   
2940   deltaphi = phi1 - phi2;
2941   if(deltaphi > PI) deltaphi -= 2*PI;
2942   if(deltaphi < -PI) deltaphi += 2*PI;
2943   deltaphi = fabs(deltaphi);
2944
2945   if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
2946   
2947   
2948    
2949   //
2950   
2951   Int_t ncl1 = first.fClusterMap.GetNbits();
2952   Int_t ncl2 = second.fClusterMap.GetNbits();
2953   Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
2954   Double_t shfrac = 0; Double_t qfactor = 0;
2955   for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
2956     if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
2957       if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
2958         sumQ++;
2959         sumCls+=2;
2960         sumSha+=2;}
2961       else {sumQ--; sumCls+=2;}
2962     }
2963     else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
2964       sumQ++;
2965       sumCls++;}
2966   }
2967   if (sumCls>0) {
2968     qfactor = sumQ*1.0/sumCls;
2969     shfrac = sumSha*1.0/sumCls;
2970   }
2971   
2972   if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
2973   
2974   
2975   return kTRUE;
2976   
2977
2978 }
2979 //________________________________________________________________________
2980 Float_t AliFourPion::Gamov(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
2981 {
2982   Float_t arg = G_Coeff/qinv;
2983   
2984   if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
2985   else {return (exp(-arg)-1)/(-arg);}
2986   
2987 }
2988 //________________________________________________________________________
2989 void AliFourPion::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
2990 {
2991   Int_t j, k;
2992   Int_t a = i2 - i1;
2993   for (Int_t i = i1; i < i2+1; i++) {
2994     j = (Int_t) (gRandom->Rndm() * a);
2995     k = iarr[j];
2996     iarr[j] = iarr[i];
2997     iarr[i] = k;
2998   }
2999 }
3000
3001
3002 //________________________________________________________________________
3003 Float_t AliFourPion::GetQinv(Float_t track1[], Float_t track2[]){
3004   
3005   Float_t qinv = sqrt( fabs(pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2)) );
3006   return qinv;
3007   
3008 }
3009 //________________________________________________________________________
3010 void AliFourPion::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
3011  
3012   Float_t p0 = track1[0] + track2[0];
3013   Float_t px = track1[1] + track2[1];
3014   Float_t py = track1[2] + track2[2];
3015   Float_t pz = track1[3] + track2[3];
3016   
3017   Float_t mt = sqrt(p0*p0 - pz*pz);
3018   Float_t pt = sqrt(px*px + py*py);
3019   
3020   Float_t v0 = track1[0] - track2[0];
3021   Float_t vx = track1[1] - track2[1];
3022   Float_t vy = track1[2] - track2[2];
3023   Float_t vz = track1[3] - track2[3];
3024   
3025   qout = (px*vx + py*vy)/pt;
3026   qside = (px*vy - py*vx)/pt;
3027   qlong = (p0*vz - pz*v0)/mt;
3028 }
3029 //________________________________________________________________________
3030 void AliFourPion::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliFourPion::fKbinsT][AliFourPion::fCentBins]){
3031
3032   if(legoCase){
3033     cout<<"LEGO call to SetWeightArrays"<<endl;
3034     
3035     for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3036       for(Int_t mb=0; mb<fCentBins; mb++){
3037         fNormWeight[tKbin][mb] = (TH3F*)histos[tKbin][mb]->Clone();
3038         fNormWeight[tKbin][mb]->SetDirectory(0);
3039       }
3040     }
3041     
3042   }else{
3043     
3044     TFile *wFile = new TFile("WeightFile.root","READ");
3045     if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
3046     else cout<<"Good Weight File Found!"<<endl;
3047     
3048     for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
3049       for(Int_t mb=0; mb<fCentBins; mb++){
3050                     
3051         TString *name = new TString("Weight_Kt_");
3052         *name += tKbin;
3053         name->Append("_Ky_0");
3054         name->Append("_M_");
3055         *name += mb;
3056         name->Append("_ED_0");
3057         
3058         
3059         fNormWeight[tKbin][mb] = (TH3F*)wFile->Get(name->Data());
3060         fNormWeight[tKbin][mb]->SetDirectory(0);
3061         
3062         
3063       }//mb
3064     }//kt
3065     
3066     wFile->Close();
3067   }
3068   
3069   cout<<"Done reading weight file"<<endl;
3070   
3071 }
3072 //________________________________________________________________________
3073 void AliFourPion::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
3074   
3075   Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
3076   //
3077   Float_t qOut=0,qSide=0,qLong=0;
3078   GetQosl(track1, track2, qOut, qSide, qLong);
3079   qOut = fabs(qOut);
3080   qSide = fabs(qSide);
3081   qLong = fabs(qLong);
3082   Float_t wd=0, xd=0, yd=0, zd=0;
3083   //Float_t qinvtemp=GetQinv(0,track1, track2);
3084   //
3085   
3086   if(kt < fKmeanT[0]) {fKtIndexL=0; fKtIndexH=1;}
3087   else if(kt >= fKmeanT[fKbinsT-1]) {fKtIndexL=fKbinsT-2; fKtIndexH=fKbinsT-1;}
3088   else {
3089     for(Int_t i=0; i<fKbinsT-1; i++){
3090       if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtIndexL=i; fKtIndexH=i+1; break;}
3091     }
3092   }
3093   wd = (kt-fKmeanT[fKtIndexL])/(fKmeanT[fKtIndexH]-fKmeanT[fKtIndexL]);
3094   //
3095   /////////
3096   if(qOut < fQmean[0]) {fQoIndexL=0; fQoIndexH=0; xd=0;}
3097   else if(qOut >= fQmean[kQbinsWeights-1]) {fQoIndexL=kQbinsWeights-1; fQoIndexH=kQbinsWeights-1; xd=1;}
3098   else {
3099     for(Int_t i=0; i<kQbinsWeights-1; i++){
3100       if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQoIndexL=i; fQoIndexH=i+1; break;}
3101     }
3102     xd = (qOut-fQmean[fQoIndexL])/(fQmean[fQoIndexH]-fQmean[fQoIndexL]);
3103   }
3104   //
3105   if(qSide < fQmean[0]) {fQsIndexL=0; fQsIndexH=0; yd=0;}
3106   else if(qSide >= fQmean[kQbinsWeights-1]) {fQsIndexL=kQbinsWeights-1; fQsIndexH=kQbinsWeights-1; yd=1;}
3107   else {
3108     for(Int_t i=0; i<kQbinsWeights-1; i++){
3109       if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsIndexL=i; fQsIndexH=i+1; break;}
3110     }
3111     yd = (qSide-fQmean[fQsIndexL])/(fQmean[fQsIndexH]-fQmean[fQsIndexL]);
3112   }
3113   //
3114   if(qLong < fQmean[0]) {fQlIndexL=0; fQlIndexH=0; zd=0;}
3115   else if(qLong >= fQmean[kQbinsWeights-1]) {fQlIndexL=kQbinsWeights-1; fQlIndexH=kQbinsWeights-1; zd=1;}
3116   else {
3117     for(Int_t i=0; i<kQbinsWeights-1; i++){
3118       if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlIndexL=i; fQlIndexH=i+1; break;}
3119     }
3120     zd = (qLong-fQmean[fQlIndexL])/(fQmean[fQlIndexH]-fQmean[fQlIndexL]);
3121   }
3122   //
3123
3124   
3125   //Float_t min = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3126   Float_t minErr = fNormWeight[fKtIndexL][fMbin]->GetBinError(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
3127   /*
3128   Float_t deltaW=0;
3129   // kt
3130   deltaW += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - min)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3131   // Qo 
3132   deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - min)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3133   // Qs
3134   deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - min)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3135   // Ql
3136   deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - min)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3137   //
3138   wgt = min + deltaW;
3139   */
3140   
3141  
3142   //
3143   // w interpolation (kt)
3144   Float_t c000 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*wd;
3145   Float_t c100 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*wd;
3146   Float_t c010 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*wd;
3147   Float_t c001 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*wd;
3148   Float_t c110 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*wd;
3149   Float_t c101 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*wd;
3150   Float_t c011 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*wd;
3151   Float_t c111 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*wd;
3152   // x interpolation (qOut)
3153   Float_t c00 = c000*(1-xd) + c100*xd;
3154   Float_t c10 = c010*(1-xd) + c110*xd;
3155   Float_t c01 = c001*(1-xd) + c101*xd;
3156   Float_t c11 = c011*(1-xd) + c111*xd;
3157   // y interpolation (qSide)
3158   Float_t c0 = c00*(1-yd) + c10*yd;
3159   Float_t c1 = c01*(1-yd) + c11*yd;
3160   // z interpolation (qLong)
3161   wgt = (c0*(1-zd) + c1*zd);
3162   
3163   ////
3164   
3165   // Denominator errors negligible compared to numerator so do not waste cpu time below.  
3166   //Float_t deltaWErr=0;
3167   // Kt
3168   /*
3169   deltaWErr += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
3170   // Qo 
3171   deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
3172   // Qs
3173   deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - minErr)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
3174   // Ql
3175   deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - minErr)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
3176   */
3177   wgtErr = minErr;
3178   
3179  
3180 }
3181 //________________________________________________________________________
3182 Float_t AliFourPion::MCWeight(Int_t c[2], Float_t R, Float_t fcSq, Float_t qinv, Float_t k12){
3183   
3184   Float_t radius = R/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3185   Float_t r12=radius*(1-k12/2.0);
3186   SetFSIindex(R);
3187   Float_t coulCorr12 = FSICorrelation(c[0], c[1], qinv);
3188   if(c[0]==c[1]){
3189     Float_t arg=qinv*r12;
3190     Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3191     EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3192     return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv*r12,2))*pow(EW,2))*coulCorr12);
3193   }else {
3194     return ((1-fcSq) + fcSq*coulCorr12);
3195   }
3196     
3197 }
3198 //________________________________________________________________________
3199 Float_t AliFourPion::MCWeightOSL(Int_t charge1, Int_t charge2, Int_t r, Int_t dampIndex, Float_t qinv, Float_t qo, Float_t qs, Float_t ql){
3200   
3201   Float_t radiusOut = Float_t(r)/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
3202   Float_t radiusSide = radiusOut;
3203   Float_t radiusLong = radiusOut;
3204   Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
3205   Float_t coulCorr12 = FSICorrelation(charge1, charge2, qinv);
3206   if(charge1==charge2){
3207     return ((1-myDamp) + myDamp*(1 + exp(-pow(qo*radiusOut,2)) * exp(-pow(qs*radiusSide,2)) * exp(-pow(ql*radiusLong,2)))*coulCorr12);
3208   }else {
3209     return ((1-myDamp) + myDamp*coulCorr12);
3210   }
3211     
3212 }
3213
3214 //________________________________________________________________________
3215 Float_t AliFourPion::MCWeight3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3], Float_t kT[3]){
3216   // FSI + QS correlations
3217   if(term==5) return 1.0;
3218   
3219   Float_t radius=R/0.19733;
3220   Float_t r12=radius*(1-kT[0]/2.0);
3221   Float_t r13=radius*(1-kT[1]/2.0);
3222   Float_t r23=radius*(1-kT[2]/2.0);
3223  
3224   Float_t fc = sqrt(fcSq);
3225   
3226   SetFSIindex(R);
3227   Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3228   Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3229   Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3230   
3231   if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3232     Float_t arg12=qinv[0]*r12;
3233     Float_t arg13=qinv[1]*r13;
3234     Float_t arg23=qinv[2]*r23;
3235     Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
3236     EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
3237     Float_t EW13 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
3238     EW13 += kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
3239     Float_t EW23 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
3240     EW23 += kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
3241     if(term==1){
3242       Float_t C3QS = 1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2) + exp(-pow(qinv[1]*r13,2))*pow(EW13,2) + exp(-pow(qinv[2]*r23,2))*pow(EW23,2);
3243       C3QS += 2*exp(-(pow(r12,2)*pow(qinv[0],2) + pow(r13,2)*pow(qinv[1],2) + pow(r23,2)*pow(qinv[2],2))/2.)*EW12*EW13*EW23;
3244       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3245       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12;
3246       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13;
3247       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23;
3248       C3 += pow(fc,3)*C3QS*Kfactor12*Kfactor13*Kfactor23;
3249       return C3;
3250     }else if(term==2){
3251       return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12);
3252     }else if(term==3){
3253       return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13);
3254     }else if(term==4){
3255       return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23);
3256     }else return 1.0;
3257     
3258   }else{// mixed charge case
3259     Float_t arg=qinv[0]*r12;
3260     Float_t KfactorSC = Kfactor12;
3261     Float_t KfactorMC1 = Kfactor13;
3262     Float_t KfactorMC2 = Kfactor23;
3263     if(c[0]==c[2]) {arg=qinv[1]*r13; KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;} 
3264     if(c[1]==c[2]) {arg=qinv[2]*r23; KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3265     Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
3266     EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
3267     if(term==1){
3268       Float_t C3QS = 1 + exp(-pow(arg,2))*pow(EW,2);
3269       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3270       C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(arg,2))*pow(EW,2))*KfactorSC;
3271       C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3272       C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3273       C3 += pow(fc,3)*C3QS*KfactorSC*KfactorMC1*KfactorMC2;
3274       return C3;
3275     }else if(term==2){
3276       if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3277       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3278     }else if(term==3){
3279       return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3280     }else if(term==4){
3281       if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
3282       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3283     }else return 1.0;
3284   }
3285   
3286 }
3287 //________________________________________________________________________
3288 Float_t AliFourPion::MCWeightFSI3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3]){
3289   // FSI only (no QS correlations)
3290   if(term==5) return 1.0;
3291   
3292   Float_t fc = sqrt(fcSq);
3293   SetFSIindex(R);
3294   Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3295   Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3296   Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
3297   
3298   if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
3299     if(term==1){
3300       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3301       C3 += pow(fc,2)*(1-fc)*Kfactor12;
3302       C3 += pow(fc,2)*(1-fc)*Kfactor13;
3303       C3 += pow(fc,2)*(1-fc)*Kfactor23;
3304       C3 += pow(fc,3)*Kfactor12*Kfactor13*Kfactor23;
3305       return C3;
3306     }else if(term==2){
3307       return ((1-fcSq) + fcSq*Kfactor12);
3308     }else if(term==3){
3309       return ((1-fcSq) + fcSq*Kfactor13);
3310     }else if(term==4){
3311       return ((1-fcSq) + fcSq*Kfactor23);
3312     }else return 1.0;
3313     
3314   }else{// mixed charge case
3315     Float_t KfactorSC = Kfactor12;
3316     Float_t KfactorMC1 = Kfactor13;
3317     Float_t KfactorMC2 = Kfactor23;
3318     if(c[0]==c[2]) {KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;} 
3319     if(c[1]==c[2]) {KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
3320     if(term==1){
3321       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
3322       C3 += pow(fc,2)*(1-fc)*KfactorSC;
3323       C3 += pow(fc,2)*(1-fc)*KfactorMC1;
3324       C3 += pow(fc,2)*(1-fc)*KfactorMC2;
3325       C3 += pow(fc,3)*KfactorSC*KfactorMC1*KfactorMC2;
3326       return C3;
3327     }else if(term==2){
3328       if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*KfactorSC);
3329       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3330     }else if(term==3){
3331       return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3332     }else if(term==4){
3333       if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*KfactorSC);
3334       else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
3335     }else return 1.0;
3336   }
3337   
3338 }
3339 //________________________________________________________________________
3340 Float_t AliFourPion::MCWeight4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6], Float_t kT[6]){
3341   if(term==13) return 1.0;
3342
3343   // Charge ordering:
3344   // ----, ---+, --++, -+++, ++++
3345   //
3346   // term ordering:
3347   // Term 1: 1-2-3-4  (all same-event)
3348   // Term 2: 1-2-3 4 (particle 4 from different event)
3349   // Term 3: 1-2-4 3 (particle 3 from different event)
3350   // Term 4: 1-3-4 2 (particle 2 from different event)
3351   // Term 5: 2-3-4 1 (particle 1 from different event)
3352   // Term 6: 1-2 3 4 (particle 1 and 2 from same event)
3353   // Term 7: 1-3 2 4
3354   // Term 8: 1-4 2 3
3355   // Term 9: 2-3 1 4
3356   // Term 10: 2-4 1 3
3357   // Term 11: 3-4 1 2
3358   // Term 12: 1 2 3 4 (all from different events)
3359
3360   Float_t radius = R/0.19733;
3361   Float_t r[6]={0};
3362   r[0]=radius*(1-kT[0]/2.0);
3363   r[1]=radius*(1-kT[1]/2.0);
3364   r[2]=radius*(1-kT[2]/2.0);
3365   r[3]=radius*(1-kT[3]/2.0);
3366   r[4]=radius*(1-kT[4]/2.0);
3367   r[5]=radius*(1-kT[5]/2.0);
3368     
3369   Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
3370  
3371   Float_t fc = sqrt(fcSq);
3372   SetFSIindex(R);
3373   Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
3374   Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
3375   Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
3376   Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
3377   Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
3378   Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
3379   Float_t arg12=qinv[0]*r[0];
3380   Float_t arg13=qinv[1]*r[1];
3381   Float_t arg14=qinv[2]*r[2];
3382   Float_t arg23=qinv[3]*r[3];
3383   Float_t arg24=qinv[4]*r[4];
3384   Float_t arg34=qinv[5]*r[5];
3385   // Exchange Amplitudes
3386   Float_t EA12 = exp(-pow(arg12,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12) + kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12));
3387   Float_t EA13 = exp(-pow(arg13,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13) + kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12));
3388   Float_t EA14 = exp(-pow(arg14,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg14,3) - 12.*arg14) + kappa4/(24.*pow(2.,2))*(16.*pow(arg14,4) -48.*pow(arg14,2) + 12));
3389   Float_t EA23 = exp(-pow(arg23,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23) + kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12));
3390   Float_t EA24 = exp(-pow(arg24,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg24,3) - 12.*arg24) + kappa4/(24.*pow(2.,2))*(16.*pow(arg24,4) -48.*pow(arg24,2) + 12));
3391   Float_t EA34 = exp(-pow(arg34,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg34,3) - 12.*arg34) + kappa4/(24.*pow(2.,2))*(16.*pow(arg34,4) -48.*pow(arg34,2) + 12));