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