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