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