Code for Three Pion Radii Analysis
authordgangadh <dgangadh@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Aug 2013 12:39:49 +0000 (12:39 +0000)
committerdgangadh <dgangadh@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Aug 2013 12:39:49 +0000 (12:39 +0000)
PWGCF/FEMTOSCOPY/Chaoticity/AliThreePionRadii.cxx [new file with mode: 0644]
PWGCF/FEMTOSCOPY/Chaoticity/AliThreePionRadii.h [new file with mode: 0644]

diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliThreePionRadii.cxx b/PWGCF/FEMTOSCOPY/Chaoticity/AliThreePionRadii.cxx
new file mode 100644 (file)
index 0000000..cfe5988
--- /dev/null
@@ -0,0 +1,3143 @@
+#include <iostream>
+#include <math.h>
+#include "TChain.h"
+#include "TFile.h"
+#include "TKey.h"
+#include "TObject.h"
+#include "TObjString.h"
+#include "TList.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TH3D.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TCanvas.h"
+#include "TRandom3.h"
+#include "TF1.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAODMCParticle.h"
+#include "AliAODTracklets.h"
+#include "AliThreePionRadii.h"
+
+#define PI 3.1415927
+#define G_Coeff 0.006399 // 2*pi*alpha*M_pion
+#define kappa3 0.2 // kappa3 Edgeworth coefficient (non-Gaussian features of C2)
+#define kappa4 0.45 // kappa4 Edgeworth coefficient (non-Gaussian features of C2)
+
+
+// Author: Dhevan Gangadharan
+
+ClassImp(AliThreePionRadii)
+
+//________________________________________________________________________
+AliThreePionRadii::AliThreePionRadii():
+AliAnalysisTaskSE(),
+  fname(0),
+  fAOD(0x0), 
+  fOutputList(0x0),
+  fPIDResponse(0x0),
+  fEC(0x0),
+  fEvt(0x0),
+  fTempStruct(0x0),
+  fRandomNumber(0x0),
+  fLEGO(kTRUE),
+  fMCcase(kFALSE),
+  fAODcase(kTRUE),
+  fPbPbcase(kTRUE),
+  fGenerateSignal(kFALSE),
+  fPdensityPairCut(kTRUE),
+  fRMax(11),
+  fFilterBit(7),
+  fMaxChi2NDF(10),
+  fMinTPCncls(0),
+  fBfield(0),
+  fMbin(0),
+  fFSIindex(0),
+  fEDbin(0),
+  fMbins(fCentBins),
+  fMultLimit(0),
+  fCentBinLowLimit(0),
+  fCentBinHighLimit(1),
+  fEventCounter(0),
+  fEventsToMix(0),
+  fZvertexBins(0),
+  fMultLimits(),
+  fQcut(),
+  fQLowerCut(0),
+  fQlimitC2(2.0),
+  fQbinsC2(400),
+  fNormQcutLow(),
+  fNormQcutHigh(),
+  fKupperBound(0),
+  fQupperBound(0),
+  fQupperBoundWeights(0),
+  fQbins(0),
+  fQbinsWeights(0),
+  fQstep(0),
+  fQstepWeights(0),
+  fQmean(),
+  fDampStart(0),
+  fDampStep(0),
+  fTPCTOFboundry(0),
+  fTOFboundry(0),
+  fSigmaCutTPC(2.0),
+  fSigmaCutTOF(2.0),
+  fMinSepPairEta(0.03),
+  fMinSepPairPhi(0.04),
+  fShareQuality(0),
+  fShareFraction(0),
+  fTrueMassP(0), 
+  fTrueMassPi(0), 
+  fTrueMassK(0), 
+  fTrueMassKs(0), 
+  fTrueMassLam(0),
+  fDummyB(0),
+  fDefaultsCharMult(),
+  fDefaultsCharSE(),
+  fDefaultsCharME(),
+  fDefaultsInt(),
+  fPairLocationSE(),
+  fPairLocationME(),
+  fTripletSkip1(),
+  fTripletSkip2(),
+  fOtherPairLocation1(),
+  fOtherPairLocation2(),
+  fNormPairSwitch(),
+  fPairSplitCut(),
+  fNormPairs()
+{
+  // Default constructor
+  for(Int_t mb=0; mb<fMbins; mb++){
+    for(Int_t edB=0; edB<fEDbins; edB++){
+      for(Int_t c1=0; c1<2; c1++){
+       for(Int_t c2=0; c2<2; c2++){
+         for(Int_t sc=0; sc<kSCLimit2; sc++){
+           for(Int_t term=0; term<2; term++){
+             
+             Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
+             
+             Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
+             Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
+             
+           }// term_2
+         }// SC_2
+         
+         for(Int_t c3=0; c3<2; c3++){
+           for(Int_t sc=0; sc<kSCLimit3; sc++){
+             for(Int_t term=0; term<5; term++){
+               
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
+               
+             }// term_3
+           }// SC_3
+         }//c3
+       }//c2
+      }//c1
+      
+    }// ED
+  }// Mbin
+  
+  // Initialize 3-pion FSI histograms
+  for(Int_t i=0; i<10; i++){
+    fFSI2SS[i]=0x0; 
+    fFSI2OS[i]=0x0;
+  }
+
+}
+//________________________________________________________________________
+AliThreePionRadii::AliThreePionRadii(const Char_t *name) 
+: AliAnalysisTaskSE(name), 
+  fname(name),
+  fAOD(0x0), 
+  fOutputList(0x0),
+  fPIDResponse(0x0),
+  fEC(0x0),
+  fEvt(0x0),
+  fTempStruct(0x0),
+  fRandomNumber(0x0),
+  fLEGO(kTRUE),
+  fMCcase(kFALSE),
+  fAODcase(kTRUE),
+  fPbPbcase(kTRUE),
+  fGenerateSignal(kFALSE),
+  fPdensityPairCut(kTRUE),
+  fRMax(11),
+  fFilterBit(7),
+  fMaxChi2NDF(10),
+  fMinTPCncls(0),
+  fBfield(0),
+  fMbin(0),
+  fFSIindex(0),
+  fEDbin(0),
+  fMbins(fCentBins),
+  fMultLimit(0),
+  fCentBinLowLimit(0),
+  fCentBinHighLimit(1),
+  fEventCounter(0),
+  fEventsToMix(0),
+  fZvertexBins(0),
+  fMultLimits(),
+  fQcut(),
+  fQLowerCut(0),
+  fQlimitC2(2.0),
+  fQbinsC2(400),
+  fNormQcutLow(),
+  fNormQcutHigh(),
+  fKupperBound(0),
+  fQupperBound(0),
+  fQupperBoundWeights(0),
+  fQbins(0),
+  fQbinsWeights(0),
+  fQstep(0),
+  fQstepWeights(0),
+  fQmean(),
+  fDampStart(0),
+  fDampStep(0),
+  fTPCTOFboundry(0),
+  fTOFboundry(0),
+  fSigmaCutTPC(2.0),
+  fSigmaCutTOF(2.0),
+  fMinSepPairEta(0.03),
+  fMinSepPairPhi(0.04),
+  fShareQuality(0),
+  fShareFraction(0),
+  fTrueMassP(0), 
+  fTrueMassPi(0), 
+  fTrueMassK(0), 
+  fTrueMassKs(0), 
+  fTrueMassLam(0),
+  fDummyB(0),
+  fDefaultsCharMult(),
+  fDefaultsCharSE(),
+  fDefaultsCharME(),
+  fDefaultsInt(),
+  fPairLocationSE(),
+  fPairLocationME(),
+  fTripletSkip1(),
+  fTripletSkip2(),
+  fOtherPairLocation1(),
+  fOtherPairLocation2(),
+  fNormPairSwitch(),
+  fPairSplitCut(),
+  fNormPairs()
+{
+  // Main constructor
+  fAODcase=kTRUE;
+  fPdensityPairCut = kTRUE;
+  
+
+  for(Int_t mb=0; mb<fMbins; mb++){
+    for(Int_t edB=0; edB<fEDbins; edB++){
+      for(Int_t c1=0; c1<2; c1++){
+       for(Int_t c2=0; c2<2; c2++){
+         for(Int_t sc=0; sc<kSCLimit2; sc++){
+           for(Int_t term=0; term<2; term++){
+             
+             Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2=0x0;
+             
+             Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
+             Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
+             
+           }// term_2
+         }// SC_2
+         
+         for(Int_t c3=0; c3<2; c3++){
+           for(Int_t sc=0; sc<kSCLimit3; sc++){
+             for(Int_t term=0; term<5; term++){
+               
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3 = 0x0;
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3 = 0x0;
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
+               Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
+               
+             }// term_3
+           }// SC_3
+         }//c3
+       }//c2
+      }//c1
+            
+    }// ED
+  }// Mbin
+  
+  // Initialize 3-pion FSI histograms
+  for(Int_t i=0; i<10; i++){
+    fFSI2SS[i]=0x0; 
+    fFSI2OS[i]=0x0;
+  }
+
+
+  DefineOutput(1, TList::Class());
+}
+//________________________________________________________________________
+AliThreePionRadii::AliThreePionRadii(const AliThreePionRadii &obj) 
+  : AliAnalysisTaskSE(obj.fname),
+    fname(obj.fname),
+    fAOD(obj.fAOD), 
+    //fESD(obj.fESD), 
+    fOutputList(obj.fOutputList),
+    fPIDResponse(obj.fPIDResponse),
+    fEC(obj.fEC),
+    fEvt(obj.fEvt),
+    fTempStruct(obj.fTempStruct),
+    fRandomNumber(obj.fRandomNumber),
+    fLEGO(obj.fLEGO),
+    fMCcase(obj.fMCcase),
+    fAODcase(obj.fAODcase),
+    fPbPbcase(obj.fPbPbcase),
+    fGenerateSignal(obj.fGenerateSignal),
+    fPdensityPairCut(obj.fPdensityPairCut),
+    fRMax(obj.fRMax),
+    fFilterBit(obj.fFilterBit),
+    fMaxChi2NDF(obj.fMaxChi2NDF),
+    fMinTPCncls(obj.fMinTPCncls),
+    fBfield(obj.fBfield),
+    fMbin(obj.fMbin),
+    fFSIindex(obj.fFSIindex),
+    fEDbin(obj.fEDbin),
+    fMbins(obj.fMbins),
+    fMultLimit(obj.fMultLimit),
+    fCentBinLowLimit(obj.fCentBinLowLimit),
+    fCentBinHighLimit(obj.fCentBinHighLimit),
+    fEventCounter(obj.fEventCounter),
+    fEventsToMix(obj.fEventsToMix),
+    fZvertexBins(obj.fZvertexBins),
+    fMultLimits(),
+    fQcut(),
+    fQLowerCut(obj.fQLowerCut),
+    fQlimitC2(obj.fQlimitC2),
+    fQbinsC2(obj.fQbinsC2),
+    fNormQcutLow(),
+    fNormQcutHigh(),
+    fKupperBound(obj.fKupperBound),
+    fQupperBound(obj.fQupperBound),
+    fQupperBoundWeights(obj.fQupperBoundWeights),
+    fQbins(obj.fQbins),
+    fQbinsWeights(obj.fQbinsWeights),
+    fQstep(obj.fQstep),
+    fQstepWeights(obj.fQstepWeights),
+    fQmean(),
+    fDampStart(obj.fDampStart),
+    fDampStep(obj.fDampStep),
+    fTPCTOFboundry(obj.fTPCTOFboundry),
+    fTOFboundry(obj.fTOFboundry),
+    fSigmaCutTPC(obj.fSigmaCutTPC),
+    fSigmaCutTOF(obj.fSigmaCutTOF),
+    fMinSepPairEta(obj.fMinSepPairEta),
+    fMinSepPairPhi(obj.fMinSepPairPhi),
+    fShareQuality(obj.fShareQuality),
+    fShareFraction(obj.fShareFraction),
+    fTrueMassP(obj.fTrueMassP), 
+    fTrueMassPi(obj.fTrueMassPi), 
+    fTrueMassK(obj.fTrueMassK), 
+    fTrueMassKs(obj.fTrueMassKs), 
+    fTrueMassLam(obj.fTrueMassLam),
+    fDummyB(obj.fDummyB),
+    fDefaultsCharMult(),
+    fDefaultsCharSE(),
+    fDefaultsCharME(),
+    fDefaultsInt(),
+    fPairLocationSE(),
+    fPairLocationME(),
+    fTripletSkip1(),
+    fTripletSkip2(),
+    fOtherPairLocation1(),
+    fOtherPairLocation2(),
+    fNormPairSwitch(),
+    fPairSplitCut(),
+    fNormPairs()
+{
+  // Copy Constructor
+  
+  for(Int_t i=0; i<10; i++){
+    fFSI2SS[i]=obj.fFSI2SS[i]; 
+    fFSI2OS[i]=obj.fFSI2OS[i];
+  }
+
+}
+//________________________________________________________________________
+AliThreePionRadii &AliThreePionRadii::operator=(const AliThreePionRadii &obj) 
+{
+  // Assignment operator  
+  if (this == &obj)
+    return *this;
+
+  fname = obj.fname;
+  fAOD = obj.fAOD; 
+  fOutputList = obj.fOutputList;
+  fPIDResponse = obj.fPIDResponse;
+  fEC = obj.fEC;
+  fEvt = obj.fEvt;
+  fTempStruct = obj.fTempStruct;
+  fRandomNumber = obj.fRandomNumber;
+  fLEGO = fLEGO;
+  fMCcase = obj.fMCcase;
+  fAODcase = obj.fAODcase;
+  fPbPbcase = obj.fPbPbcase; 
+  fGenerateSignal = obj.fGenerateSignal;
+  fPdensityPairCut = obj.fPdensityPairCut;
+  fRMax = obj.fRMax;
+  fFilterBit = obj.fFilterBit;
+  fMaxChi2NDF = obj.fMaxChi2NDF;
+  fMinTPCncls = obj.fMinTPCncls;
+  fBfield = obj.fBfield;
+  fMbin = obj.fMbin;
+  fFSIindex = obj.fFSIindex;
+  fEDbin = obj.fEDbin;
+  fMbins = obj.fMbins;
+  fMultLimit = obj.fMultLimit;
+  fCentBinLowLimit = obj.fCentBinLowLimit;
+  fCentBinHighLimit = obj.fCentBinHighLimit;
+  fEventCounter = obj.fEventCounter;
+  fEventsToMix = obj.fEventsToMix;
+  fZvertexBins = obj.fZvertexBins;
+  fQLowerCut = obj.fQLowerCut;
+  fQlimitC2 = obj.fQlimitC2;
+  fQbinsC2 = obj.fQbinsC2;
+  fKupperBound = obj.fKupperBound;
+  fQupperBound = obj.fQupperBound;
+  fQupperBoundWeights = obj.fQupperBoundWeights;
+  fQbins = obj.fQbins;
+  fQbinsWeights = obj.fQbinsWeights;
+  fQstep = obj.fQstep;
+  fQstepWeights = obj.fQstepWeights;
+  fDampStart = obj.fDampStart;
+  fDampStep = obj.fDampStep;
+  fTPCTOFboundry = obj.fTPCTOFboundry;
+  fTOFboundry = obj.fTOFboundry;
+  fSigmaCutTPC = obj.fSigmaCutTPC;
+  fSigmaCutTOF = obj.fSigmaCutTOF;
+  fMinSepPairEta = obj.fMinSepPairEta;
+  fMinSepPairPhi = obj.fMinSepPairPhi;
+  fShareQuality = obj.fShareQuality;
+  fShareFraction = obj.fShareFraction;
+  fTrueMassP = obj.fTrueMassP; 
+  fTrueMassPi = obj.fTrueMassPi; 
+  fTrueMassK = obj.fTrueMassK; 
+  fTrueMassKs = obj.fTrueMassKs; 
+  fTrueMassLam = obj.fTrueMassLam;
+  fDummyB = obj.fDummyB;
+  for(Int_t i=0; i<10; i++){
+    fFSI2SS[i]=obj.fFSI2SS[i]; 
+    fFSI2OS[i]=obj.fFSI2OS[i];
+  }
+  
+  return (*this);
+}
+//________________________________________________________________________
+AliThreePionRadii::~AliThreePionRadii()
+{
+  // Destructor
+  if(fAOD) delete fAOD; 
+  //if(fESD) delete fESD; 
+  if(fOutputList) delete fOutputList;
+  if(fPIDResponse) delete fPIDResponse;
+  if(fEC) delete fEC;
+  if(fEvt) delete fEvt;
+  if(fTempStruct) delete [] fTempStruct;
+  if(fRandomNumber) delete fRandomNumber;
+  for(Int_t i=0; i<fMultLimit; i++){
+    if(fPairLocationSE[i]) delete [] fPairLocationSE[i];
+    if(fPairLocationME[i]) delete [] fPairLocationME[i];
+    for(Int_t j=0; j<2; j++){
+      if(fOtherPairLocation1[j][i]) delete [] fOtherPairLocation1[j][i];
+      if(fOtherPairLocation2[j][i]) delete [] fOtherPairLocation2[j][i];
+    }
+    for(Int_t j=0; j<3; j++) if(fNormPairSwitch[j][i]) delete [] fNormPairSwitch[j][i];
+    for(Int_t j=0; j<4; j++) if(fPairSplitCut[j][i]) delete [] fPairSplitCut[j][i];
+  }
+  for(Int_t i=0; i<kPairLimit; i++) if(fTripletSkip1[i]) delete [] fTripletSkip1[i];
+  for(Int_t i=0; i<2*kPairLimit; i++) if(fTripletSkip2[i]) delete [] fTripletSkip2[i];
+  for(Int_t i=0; i<3; i++) if(fNormPairs[i]) delete [] fNormPairs[i];
+  //
+  for(Int_t mb=0; mb<fMbins; mb++){
+    if(fPbPbcase && ((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit))) continue;
+    for(Int_t edB=0; edB<fEDbins; edB++){
+      for(Int_t c1=0; c1<2; c1++){
+       for(Int_t c2=0; c2<2; c2++){
+         for(Int_t sc=0; sc<kSCLimit2; sc++){
+           for(Int_t term=0; term<2; term++){
+             
+             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;
+             
+             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;
+             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;
+             //
+             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;
+             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;
+           }// term_2
+         }// SC_2
+         
+         for(Int_t c3=0; c3<2; c3++){
+           for(Int_t sc=0; sc<kSCLimit3; sc++){
+             for(Int_t term=0; term<5; term++){
+               
+               if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fExplicit3;
+               if(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3) delete Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNormEx3;
+               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;
+               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;
+                               
+             }// term_3
+           }// SC_3
+         }//c3
+       }//c2
+      }//c1
+      
+    }// ED
+  }// Mbin
+  
+   
+  for(Int_t i=0; i<10; i++){
+    if(fFSI2SS[i]) delete fFSI2SS[i]; 
+    if(fFSI2OS[i]) delete fFSI2OS[i];
+  }
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::ParInit()
+{
+  cout<<"AliThreePionRadii MyInit() call"<<endl;
+  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;
+
+  fRandomNumber = new TRandom3();
+  fRandomNumber->SetSeed(0);
+    
+  //
+  fEventCounter=0;
+  if(fPdensityPairCut) fEventsToMix=2;
+  else fEventsToMix=0;
+  fZvertexBins=2;//2
+  
+  fTPCTOFboundry = 0.6;// TPC pid used below this momentum, TOF above but below TOF_boundry
+  fTOFboundry = 2.1;// TOF pid used below this momentum
+  
+  ////////////////////////////////////////////////
+  // PadRow Pair Cuts
+  fShareQuality = .5;// max
+  fShareFraction = .05;// max
+  ////////////////////////////////////////////////
+  
+  
+  //fMultLimits[0]=0, fMultLimits[1]=5, fMultLimits[2]=10, fMultLimits[3]=15, fMultLimits[4]=20, fMultLimits[5]=25;
+  //fMultLimits[6]=30, fMultLimits[7]=35, fMultLimits[8]=40, fMultLimits[9]=45, fMultLimits[10]=kMultLimitPP;
+  for(Int_t index=0; index<fCentBins+1; index++){// (dNtracklet/deta)^(1/3)
+    fMultLimits[index]=index;
+  }
+
+  if(fPbPbcase && fCentBinLowLimit < 10) {// PbPb 0-50%
+    fMultLimit=kMultLimitPbPb; 
+    fMbins=fCentBins; 
+    fQcut[0]=0.1;//pi-pi, pi-k, pi-p
+    fQcut[1]=0.1;//k-k
+    fQcut[2]=0.6;//the rest
+    fNormQcutLow[0] = 0.15;//0.15
+    fNormQcutHigh[0] = 0.175;//0.175
+    fNormQcutLow[1] = 1.34;//1.34
+    fNormQcutHigh[1] = 1.4;//1.4
+    fNormQcutLow[2] = 1.1;//1.1
+    fNormQcutHigh[2] = 1.4;//1.4
+    //
+    fQlimitC2 = 2.0;
+    fQbinsC2 = 400;
+    fQupperBoundWeights = 2*fQcut[0];
+    fQupperBound = fQcut[0];
+    fQbins = kQbins;
+    fQbinsWeights = kQbinsWeights;
+    fQstep = fQupperBound/Float_t(fQbins);
+    fQstepWeights = fQupperBoundWeights/Float_t(fQbinsWeights);
+    for(Int_t i=0; i<fQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
+    //
+    fDampStart = 0.5;// was 0.3
+    fDampStep = 0.02;
+  }else if(fPbPbcase && fCentBinLowLimit >= 10) {// PbPb 50-100%
+    fMultLimit=kMultLimitPbPb;
+    fMbins=fCentBins;
+    fQcut[0]=0.2;//pi-pi, pi-k, pi-p
+    fQcut[1]=0.2;//k-k
+    fQcut[2]=1.2;//the rest
+    fNormQcutLow[0] = 0.3;//0.15
+    fNormQcutHigh[0] = 0.35;//0.175
+    fNormQcutLow[1] = 1.34;//1.34
+    fNormQcutHigh[1] = 1.4;//1.4
+    fNormQcutLow[2] = 1.1;//1.1
+    fNormQcutHigh[2] = 1.4;//1.4
+    //
+    fQlimitC2 = 2.0;
+    fQbinsC2 = 400;
+    fQupperBoundWeights = fQcut[0];
+    fQupperBound = fQcut[0];
+    fQbins = 2*kQbins;
+    fQbinsWeights = 2*kQbinsWeights;
+    fQstep = fQupperBound/Float_t(fQbins);
+    fQstepWeights = fQupperBoundWeights/Float_t(fQbinsWeights);
+    for(Int_t i=0; i<fQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
+    //
+    fDampStart = 0.5;// was 0.3
+    fDampStep = 0.02;
+  }else {// pp or pPb
+    fMultLimit=kMultLimitPP;
+    fMbins=5;
+    fQcut[0]=2.0;// 0.4
+    fQcut[1]=2.0;
+    fQcut[2]=2.0;
+    fNormQcutLow[0] = 1.0;
+    fNormQcutHigh[0] = 1.2;// 1.5
+    fNormQcutLow[1] = 1.0;
+    fNormQcutHigh[1] = 1.2;
+    fNormQcutLow[2] = 1.0;
+    fNormQcutHigh[2] = 1.2;
+    //
+    fQlimitC2 = 2.0;
+    fQbinsC2 = 200;
+    fQupperBoundWeights = 0.4;
+    fQupperBound = 0.4;
+    fQbins = kQbinsPP;
+    fQbinsWeights = kQbinsWeightsPP;
+    fQstep = fQupperBound/Float_t(kQbinsPP);
+    fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeightsPP);
+    for(Int_t i=0; i<kQbinsWeightsPP; i++) {fQmeanPP[i]=(i+0.5)*fQstepWeights;}
+    //
+    fDampStart = 0.5;// was 0.3
+    fDampStep = 0.02;
+  }
+
+  fQLowerCut = 0.005;// was 0.005
+  fKupperBound = 1.0;
+  //
+  
+  fEC = new AliChaoticityEventCollection **[fZvertexBins];
+  for(UShort_t i=0; i<fZvertexBins; i++){
+    
+    fEC[i] = new AliChaoticityEventCollection *[fMbins];
+
+    for(UShort_t j=0; j<fMbins; j++){
+      
+      fEC[i][j] = new AliChaoticityEventCollection(fEventsToMix+1, fMultLimit, kPairLimit, kMCarrayLimit, fMCcase);
+    }
+  }
+  
+    
+  for(Int_t i=0; i<fMultLimit; i++) fDefaultsCharMult[i]='0';
+  for(Int_t i=0; i<kPairLimit; i++) fDefaultsCharSE[i]='0';
+  for(Int_t i=0; i<2*kPairLimit; i++) fDefaultsCharME[i]='0';
+  for(Int_t i=0; i<fMultLimit; i++) fDefaultsInt[i]=-1;
+  for(Int_t i=0; i<fMultLimit; i++) fPairLocationSE[i] = new TArrayI(fMultLimit,fDefaultsInt);
+  for(Int_t i=0; i<fMultLimit; i++) fPairLocationME[i] = new TArrayI(fMultLimit,fDefaultsInt);
+  for(Int_t i=0; i<kPairLimit; i++) fTripletSkip1[i] = new TArrayC(fMultLimit,fDefaultsCharSE);
+  for(Int_t i=0; i<2*kPairLimit; i++) fTripletSkip2[i] = new TArrayC(fMultLimit,fDefaultsCharME);
+
+  // Normalization utilities
+  for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
+  for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation1[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
+  for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[0][i] = new TArrayI(fMultLimit,fDefaultsInt);
+  for(Int_t i=0; i<fMultLimit; i++) fOtherPairLocation2[1][i] = new TArrayI(fMultLimit,fDefaultsInt);
+  for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
+  for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
+  for(Int_t i=0; i<fMultLimit; i++) fNormPairSwitch[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);
+  // Track Merging/Splitting utilities
+  for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[0][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P11
+  for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[1][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P12
+  for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[2][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P13
+  for(Int_t i=0; i<fMultLimit; i++) fPairSplitCut[3][i] = new TArrayC(fMultLimit,fDefaultsCharMult);// P23
+
+  
+  fNormPairs[0] = new AliChaoticityNormPairStruct[kNormPairLimit];
+  fNormPairs[1] = new AliChaoticityNormPairStruct[kNormPairLimit];
+  
+
+  fTempStruct = new AliChaoticityTrackStruct[fMultLimit];
+   
+   
+  fTrueMassP=0.93827, fTrueMassPi=0.13957, fTrueMassK=0.493677, fTrueMassKs=0.497614, fTrueMassLam=1.11568;
+
+  
+
+  // Set weights, Coulomb corrections, and Momentum resolution corrections manually if not on LEGO
+  if(!fLEGO) {
+    SetFSICorrelations(fLEGO);// Read in 2-particle and 3-particle FSI correlations
+  }
+  
+  /////////////////////////////////////////////
+  /////////////////////////////////////////////
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+  
+  ParInit();// Initialize my settings
+
+
+  fOutputList = new TList();
+  fOutputList->SetOwner();
+  
+  TH3F *fVertexDist = new TH3F("fVertexDist","Vertex Distribution",20,-1,1, 20,-1,1, 600,-30,30);
+  fVertexDist->GetXaxis()->SetTitle("X Vertex (cm)");
+  fVertexDist->GetYaxis()->SetTitle("Y Vertex (cm)");
+  fVertexDist->GetZaxis()->SetTitle("Z Vertex (cm)");
+  fOutputList->Add(fVertexDist);
+  
+  
+  TH2F *fDCAxyDistPlus = new TH2F("fDCAxyDistPlus","DCA distribution",300,0,3., 50,0,5);
+  fOutputList->Add(fDCAxyDistPlus);
+  TH2F *fDCAzDistPlus = new TH2F("fDCAzDistPlus","DCA distribution",300,0,3., 50,0,5);
+  fOutputList->Add(fDCAzDistPlus);
+  TH2F *fDCAxyDistMinus = new TH2F("fDCAxyDistMinus","DCA distribution",300,0,3., 50,0,5);
+  fOutputList->Add(fDCAxyDistMinus);
+  TH2F *fDCAzDistMinus = new TH2F("fDCAzDistMinus","DCA distribution",300,0,3., 50,0,5);
+  fOutputList->Add(fDCAzDistMinus);
+  
+  
+  TH1F *fEvents1 = new TH1F("fEvents1","Events vs. fMbin",fMbins,.5,fMbins+.5);
+  fOutputList->Add(fEvents1);
+  TH1F *fEvents2 = new TH1F("fEvents2","Events vs. fMbin",fMbins,.5,fMbins+.5);
+  fOutputList->Add(fEvents2);
+  
+  TH1F *fMultDist1 = new TH1F("fMultDist1","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
+  fMultDist1->GetXaxis()->SetTitle("Multiplicity");
+  fOutputList->Add(fMultDist1);
+  
+  TH1F *fMultDist2 = new TH1F("fMultDist2","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
+  fMultDist2->GetXaxis()->SetTitle("Multiplicity");
+  fOutputList->Add(fMultDist2);
+  
+  TH1F *fMultDist3 = new TH1F("fMultDist3","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
+  fMultDist3->GetXaxis()->SetTitle("Multiplicity");
+  fOutputList->Add(fMultDist3);
+  
+  TH3F *fPtEtaDist = new TH3F("fPtEtaDist","fPtEtaDist",2,-1.1,1.1, 300,0,3., 28,-1.4,1.4);
+  fOutputList->Add(fPtEtaDist);
+
+  TH3F *fPhiPtDist = new TH3F("fPhiPtDist","fPhiPtDist",2,-1.1,1.1, 120,0,2*PI, 300,0,3.);
+  fOutputList->Add(fPhiPtDist);
+  
+  TH3F *fTOFResponse = new TH3F("fTOFResponse","TOF relative time",20,0,100, 200,0,2, 4000,-20000,20000);
+  fOutputList->Add(fTOFResponse);
+  TH3F *fTPCResponse = new TH3F("fTPCResponse","TPCsignal",20,0,100, 200,0,2, 1000,0,1000);
+  fOutputList->Add(fTPCResponse);
+  TH1F *fRejectedPairs = new TH1F("fRejectedPairs","",200,0,2);
+  fOutputList->Add(fRejectedPairs);
+  TH1I *fRejectedEvents = new TH1I("fRejectedEvents","",fMbins,0.5,fMbins+.5);
+  fOutputList->Add(fRejectedEvents);
+  
+  TH3F *fPairsDetaDPhiNum = new TH3F("fPairsDetaDPhiNum","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
+  if(fMCcase) fOutputList->Add(fPairsDetaDPhiNum);
+  TH3F *fPairsDetaDPhiDen = new TH3F("fPairsDetaDPhiDen","",10,-.5,9.5, 200,-0.2,0.2, 600,-0.3,0.3);
+  if(fMCcase) fOutputList->Add(fPairsDetaDPhiDen);
+  TH3F *fPairsShareFracDPhiNum = new TH3F("fPairsShareFracDPhiNum","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
+  if(fMCcase) fOutputList->Add(fPairsShareFracDPhiNum);
+  TH3F *fPairsShareFracDPhiDen = new TH3F("fPairsShareFracDPhiDen","",10,-.5,9.5, 159,0,1, 600,-0.3,0.3);
+  if(fMCcase) fOutputList->Add(fPairsShareFracDPhiDen);
+  TH3D* fPairsPadRowNum = new TH3D("fPairsPadRowNum","", 20,0,1, 159,0,1, 40,0,0.2);
+  if(fMCcase) fOutputList->Add(fPairsPadRowNum);
+  TH3D* fPairsPadRowDen = new TH3D("fPairsPadRowDen","", 20,0,1, 159,0,1, 40,0,0.2);
+  if(fMCcase) fOutputList->Add(fPairsPadRowDen);
+
+
+  
+  TH3D *fPrimarySCPionPairs = new TH3D("fPrimarySCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  if(fMCcase) fOutputList->Add(fPrimarySCPionPairs);
+  TH3D *fAllSCPionPairs = new TH3D("fAllSCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  if(fMCcase) fOutputList->Add(fAllSCPionPairs);
+  TH3D *fPrimaryMCPionPairs = new TH3D("fPrimaryMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  if(fMCcase) fOutputList->Add(fPrimaryMCPionPairs);
+  TH3D *fAllMCPionPairs = new TH3D("fAllMCPionPairs","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  if(fMCcase) fOutputList->Add(fAllMCPionPairs);
+
+  TProfile *fAvgMult = new TProfile("fAvgMult","",fMbins,.5,fMbins+.5, 0,1500,"");
+  fOutputList->Add(fAvgMult);
+
+  TH2D *fTrackChi2NDF = new TH2D("fTrackChi2NDF","",20,0,100, 100,0,10);
+  fOutputList->Add(fTrackChi2NDF);
+  TH2D *fTrackTPCncls = new TH2D("fTrackTPCncls","",20,0,100, 110,50,160);
+  fOutputList->Add(fTrackTPCncls);
+
+
+  TH3D *fTPNRejects1 = new TH3D("fTPNRejects1","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
+  fOutputList->Add(fTPNRejects1);
+  TH3D *fTPNRejects2 = new TH3D("fTPNRejects2","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
+  fOutputList->Add(fTPNRejects2);
+  TH3D *fTPNRejects3 = new TH3D("fTPNRejects3","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
+  fOutputList->Add(fTPNRejects3);
+  TH3D *fTPNRejects4 = new TH3D("fTPNRejects4","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
+  fOutputList->Add(fTPNRejects4);
+  TH3D *fTPNRejects5 = new TH3D("fTPNRejects5","",fQbins,0,fQupperBound, fQbins,0,fQupperBound, fQbins,0,fQupperBound);
+  fOutputList->Add(fTPNRejects5);
+
+
+  TH3D *fKt3DistTerm1 = new TH3D("fKt3DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  TH3D *fKt3DistTerm5 = new TH3D("fKt3DistTerm5","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  fOutputList->Add(fKt3DistTerm1);
+  fOutputList->Add(fKt3DistTerm5);
+
+  TH1D *fMCWeight3DTerm1SC = new TH1D("fMCWeight3DTerm1SC","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm1SCden = new TH1D("fMCWeight3DTerm1SCden","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm2SC = new TH1D("fMCWeight3DTerm2SC","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm2SCden = new TH1D("fMCWeight3DTerm2SCden","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm1MC = new TH1D("fMCWeight3DTerm1MC","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm1MCden = new TH1D("fMCWeight3DTerm1MCden","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm2MC = new TH1D("fMCWeight3DTerm2MC","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm2MCden = new TH1D("fMCWeight3DTerm2MCden","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm3MC = new TH1D("fMCWeight3DTerm3MC","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm3MCden = new TH1D("fMCWeight3DTerm3MCden","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm4MC = new TH1D("fMCWeight3DTerm4MC","", 20,0,0.2);
+  TH1D *fMCWeight3DTerm4MCden = new TH1D("fMCWeight3DTerm4MCden","", 20,0,0.2);
+  fOutputList->Add(fMCWeight3DTerm1SC);
+  fOutputList->Add(fMCWeight3DTerm1SCden);
+  fOutputList->Add(fMCWeight3DTerm2SC);
+  fOutputList->Add(fMCWeight3DTerm2SCden);
+  fOutputList->Add(fMCWeight3DTerm1MC);
+  fOutputList->Add(fMCWeight3DTerm1MCden);
+  fOutputList->Add(fMCWeight3DTerm2MC);
+  fOutputList->Add(fMCWeight3DTerm2MCden);
+  fOutputList->Add(fMCWeight3DTerm3MC);
+  fOutputList->Add(fMCWeight3DTerm3MCden);
+  fOutputList->Add(fMCWeight3DTerm4MC);
+  fOutputList->Add(fMCWeight3DTerm4MCden);
+
+  TH2D *fdNchdEtaResponse = new TH2D("fdNchdEtaResponse","",15,0,15, 15,0,15);
+  if(fMCcase) fOutputList->Add(fdNchdEtaResponse);
+  TH2D *fdCentVsNchdEta = new TH2D("fdCentVsNchdEta","",fMbins,.5,fMbins+.5, 15,0,15);
+  if(fPbPbcase) fOutputList->Add(fdCentVsNchdEta);
+  
+  if(fPdensityPairCut){
+    
+    for(Int_t mb=0; mb<fMbins; mb++){
+      if(fPbPbcase && ((mb < fCentBinLowLimit) || (mb > fCentBinHighLimit))) continue;
+      
+      for(Int_t edB=0; edB<fEDbins; edB++){
+       for(Int_t c1=0; c1<2; c1++){
+         for(Int_t c2=0; c2<2; c2++){
+           for(Int_t sc=0; sc<kSCLimit2; sc++){
+             for(Int_t term=0; term<2; term++){
+               
+               TString *nameEx2 = new TString("Explicit2_Charge1_");
+               *nameEx2 += c1;
+               nameEx2->Append("_Charge2_");
+               *nameEx2 += c2;
+               nameEx2->Append("_SC_");
+               *nameEx2 += sc;
+               nameEx2->Append("_M_");
+               *nameEx2 += mb;
+               nameEx2->Append("_ED_");
+               *nameEx2 += edB;
+               nameEx2->Append("_Term_");
+               *nameEx2 += term+1;
+               
+               if(sc==0 || sc==3 || sc==5){
+                 if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
+               }
+               
+               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);
+               fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2);
+               TString *nameEx2QW=new TString(nameEx2->Data());
+               nameEx2QW->Append("_QW");
+               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);
+               fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fExplicit2QW);
+               TString *nameAvgP=new TString(nameEx2->Data());
+               nameAvgP->Append("_AvgP");
+               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,"");
+               fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fAvgP);
+               
+               // Momentum resolution histos
+               if(fMCcase && sc==0){
+                 TString *nameIdeal = new TString(nameEx2->Data());
+                 nameIdeal->Append("_Ideal");
+                 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, fQbinsWeights,0,fQupperBoundWeights);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fIdeal);
+                 TString *nameSmeared = new TString(nameEx2->Data());
+                 nameSmeared->Append("_Smeared");
+                 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, fQbinsWeights,0,fQupperBoundWeights);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fSmeared);
+                 //
+                 TString *nameEx2MC=new TString(nameEx2->Data());
+                 nameEx2MC->Append("_MCqinv");
+                 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"",fQbinsC2,0,fQlimitC2);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
+                 TString *nameEx2MCQW=new TString(nameEx2->Data());
+                 nameEx2MCQW->Append("_MCqinvQW");
+                 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"",fQbinsC2,0,fQlimitC2);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
+                 //
+                 TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
+                 nameEx2PIDpurityDen->Append("_PIDpurityDen");
+                 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);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
+                 TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
+                 nameEx2PIDpurityNum->Append("_PIDpurityNum");
+                 Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH2D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",20,0,1, fQbinsC2,0,fQlimitC2);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].SC[sc].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
+               }
+               
+             }// term_2
+           }// SC_2
+           
+           
+           for(Int_t c3=0; c3<2; c3++){
+             for(Int_t sc=0; sc<kSCLimit3; sc++){
+               for(Int_t term=0; term<5; term++){
+                 TString *nameEx3 = new TString("Explicit3_Charge1_");
+                 *nameEx3 += c1;
+                 nameEx3->Append("_Charge2_");
+                 *nameEx3 += c2;
+                 nameEx3->Append("_Charge3_");
+                 *nameEx3 += c3;
+                 nameEx3->Append("_SC_");
+                 *nameEx3 += sc;
+                 nameEx3->Append("_M_");
+                 *nameEx3 += mb;
+                 nameEx3->Append("_ED_");
+                 *nameEx3 += edB;
+                 nameEx3->Append("_Term_");
+                 *nameEx3 += term+1;
+                 
+                 TString *namePC3 = new TString("PairCut3_Charge1_");
+                 *namePC3 += c1;
+                 namePC3->Append("_Charge2_");
+                 *namePC3 += c2;
+                 namePC3->Append("_Charge3_");
+                 *namePC3 += c3;
+                 namePC3->Append("_SC_");
+                 *namePC3 += sc;
+                 namePC3->Append("_M_");
+                 *namePC3 += mb;
+                 namePC3->Append("_ED_");
+                 *namePC3 += edB;
+                 namePC3->Append("_Term_");
+                 *namePC3 += term+1;
+             
+                 ///////////////////////////////////////
+                 // skip degenerate histograms
+                 if(sc==0 || sc==6 || sc==9){// Identical species
+                   if( (c1+c2+c3)==1) {if(c3!=1) continue;}
+                   if( (c1+c2+c3)==2) {if(c1!=0) continue;}
+                 }else if(sc!=5){
+                   if( (c1+c2)==1) {if(c1!=0) continue;}
+                 }else {}// do nothing for pi-k-p case
+                 
+                 /////////////////////////////////////////
+             
+                 
+
+                 
+                 if(fPdensityPairCut){
+                   TString *nameNorm=new TString(namePC3->Data());
+                   nameNorm->Append("_Norm");
+                   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);
+                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fNorm3);
+                   //
+                   if(sc<=2){
+                     TString *nameQ3=new TString(namePC3->Data());
+                     nameQ3->Append("_Q3");
+                     Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTermsQ3 = new TH1D(nameQ3->Data(),"", 200,0,2.0);
+                     fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTermsQ3);
+                     //
+                     TString *name3DQ=new TString(namePC3->Data());
+                     name3DQ->Append("_3D");
+                     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);
+                     fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fTerms3);
+                     //
+                     
+                     if(sc==0 && fMCcase==kTRUE){
+                       TString *name3DMomResIdeal=new TString(namePC3->Data());
+                       name3DMomResIdeal->Append("_Ideal");
+                       Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH1D(name3DMomResIdeal->Data(),"", 200,0,2.0);
+                       fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fIdeal);
+                       TString *name3DMomResSmeared=new TString(namePC3->Data());
+                       name3DMomResSmeared->Append("_Smeared");
+                       Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH1D(name3DMomResSmeared->Data(),"", 200,0,2.0);
+                       fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[mb].EDB[edB].ThreePT[term].fSmeared);
+                     }// MCcase
+                     
+                     
+                   }// sc exclusion
+                 }// PdensityPairCut
+               }// term_3
+             }// SC_3
+           }//c3
+         }//c2
+       }//c1
+      }// ED
+    }// mbin
+  }// Pdensity Method
+
+  
+    
+  
+  
+  
+  TProfile *fQsmearMean = new TProfile("fQsmearMean","",2,0.5,2.5, -0.2,0.2,"");
+  fOutputList->Add(fQsmearMean);
+  TProfile *fQsmearSq = new TProfile("fQsmearSq","",2,0.5,2.5, -2,2,"");
+  fOutputList->Add(fQsmearSq);
+  TH1D *fQDist = new TH1D("fQDist","",200,-.2,.2);
+  fOutputList->Add(fQDist);
+  
+  
+
+  ////////////////////////////////////
+  ///////////////////////////////////  
+  
+  PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliThreePionRadii::Exec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+  //cout<<"===========  Event # "<<fEventCounter+1<<"  ==========="<<endl;
+  fEventCounter++;
+  
+  if(!fAODcase) {cout<<"ESDs not supported"<<endl; return;}
+  
+  fAOD = dynamic_cast<AliAODEvent*> (InputEvent());
+  if (!fAOD) {Printf("ERROR: fAOD not available"); return;}
+  
+  
+  // Trigger Cut
+  if(fPbPbcase){
+    if(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
+      Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+      if(!isSelected1 && !fMCcase) {return;}
+    }
+    if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
+      Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+      Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
+      Bool_t isSelected3 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
+      if(!isSelected1 && !isSelected2 && !isSelected3 && !fMCcase) {return;}
+    }
+  }
+
+  ///////////////////////////////////////////////////////////
+  const AliAODVertex *primaryVertexAOD;
+  AliCentrality *centrality;// for AODs and ESDs
+
+  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+  fPIDResponse = inputHandler->GetPIDResponse();
+
+  
+  TClonesArray *mcArray = 0x0;
+  Int_t mcdNchdEta=0;
+  if(fMCcase){
+    if(fAODcase){ 
+      mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
+      if(!mcArray || mcArray->GetEntriesFast() >= kMCarrayLimit){
+       cout<<"No MC particle branch found or Array too large!!"<<endl;
+       return;
+      }
+      
+      // Count true Nch at mid-rapidity
+      for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
+       AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
+       if(!mcParticle) continue;
+       if(fabs(mcParticle->Eta())>0.5) continue;
+       if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
+       if(!mcParticle->IsPrimary()) continue;
+       if(!mcParticle->IsPhysicalPrimary()) continue;
+       mcdNchdEta++;
+      }
+      
+    }
+  }// fMCcase
+  
+  UInt_t status=0;
+  Int_t positiveTracks=0, negativeTracks=0;
+  Int_t myTracks=0, pionCount=0, kaonCount=0, protonCount=0;
+   
+  Double_t vertex[3]={0};
+  Int_t zbin=0;
+  Double_t zstep=2*10/Double_t(fZvertexBins), zstart=-10.;
+  /////////////////////////////////////////////////
+
+  
+  Float_t centralityPercentile=0;
+  //Float_t cStep=5.0, cStart=0;
+  Int_t trackletMult = 0;
+
+  if(fAODcase){// AOD case
+    
+    if(fPbPbcase){
+      centrality = fAOD->GetCentrality();
+      centralityPercentile = centrality->GetCentralityPercentile("V0M");
+      if(centralityPercentile == 0) {cout<<"Centrality = 0, skipping event"<<endl; return;}
+      if((centralityPercentile < 5*fCentBinLowLimit) || (centralityPercentile>= 5*(fCentBinHighLimit+1))) {/*cout<<"Centrality out of Range.  Skipping Event"<<endl;*/ return;}
+      cout<<"Centrality % = "<<centralityPercentile<<endl;
+    }else{
+      //cout<<"AOD multiplicity = "<<fAOD->GetNumberOfTracks()<<endl;
+    }
+    
+
+    
+    
+    ////////////////////////////////
+    // Vertexing
+    ((TH1F*)fOutputList->FindObject("fMultDist1"))->Fill(fAOD->GetNumberOfTracks());
+    primaryVertexAOD = fAOD->GetPrimaryVertex();
+    vertex[0]=primaryVertexAOD->GetX(); vertex[1]=primaryVertexAOD->GetY(); vertex[2]=primaryVertexAOD->GetZ();
+    
+    if(fabs(vertex[2]) > 10) {cout<<"Zvertex Out of Range. Skip Event"<<endl; return;} // Z-Vertex Cut 
+    ((TH3F*)fOutputList->FindObject("fVertexDist"))->Fill(vertex[0], vertex[1], vertex[2]);
+    
+    if(fAOD->IsPileupFromSPD()) {cout<<"PileUpEvent. Skip Event"<<endl; return;} // Reject Pile-up events
+    if(primaryVertexAOD->GetNContributors() < 1) {cout<<"Bad Vertex. Skip Event"<<endl; return;}
+   
+    ((TH1F*)fOutputList->FindObject("fMultDist2"))->Fill(fAOD->GetNumberOfTracks());
+    fBfield = fAOD->GetMagneticField();
+    
+    for(Int_t i=0; i<fZvertexBins; i++){
+      if( (vertex[2] >= zstart+i*zstep) && (vertex[2] < zstart+(i+1)*zstep) ){
+       zbin=i;
+       break;
+      }
+    }
+    
+    AliAODTracklets *tracklets = (AliAODTracklets*)fAOD->GetTracklets();
+    for(Int_t trackletN=0; trackletN<tracklets->GetNumberOfTracklets(); trackletN++){
+      if(tracklets->GetTheta(trackletN) > 1.0904 && tracklets->GetTheta(trackletN) < 2.0512) trackletMult++;// |eta|<0.5 tracklets
+    }
+   
+    
+    /////////////////////////////
+    // Create Shuffled index list
+    Int_t randomIndex[fAOD->GetNumberOfTracks()];
+    for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) randomIndex[i]=i;
+    Shuffle(randomIndex,0,fAOD->GetNumberOfTracks()-1);
+    /////////////////////////////
+  
+    // Track loop
+    for (Int_t i = 0; i < fAOD->GetNumberOfTracks(); i++) {
+      AliAODTrack* aodtrack = fAOD->GetTrack(randomIndex[i]);
+      if (!aodtrack) continue;
+      if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
+    
+      status=aodtrack->GetStatus();
+      
+          
+      if(!aodtrack->TestFilterBit(BIT(fFilterBit))) continue;// AOD filterBit cut
+      
+      if(aodtrack->Pt() < 0.16) continue;
+      if(fabs(aodtrack->Eta()) > 0.8) continue;
+      
+     
+      Bool_t goodMomentum = aodtrack->GetPxPyPz( fTempStruct[myTracks].fP);
+      if(!goodMomentum) continue; 
+      aodtrack->GetXYZ( fTempStruct[myTracks].fX);
+         
+      Float_t dca2[2];
+      Float_t dca3d;
+
+      dca2[0] = sqrt( pow(fTempStruct[myTracks].fX[0] - vertex[0],2) + pow(fTempStruct[myTracks].fX[1] - vertex[1],2));
+      dca2[1] = sqrt( pow(fTempStruct[myTracks].fX[2] - vertex[2],2));
+      dca3d = sqrt( pow(dca2[0],2) + pow(dca2[1],2));
+             
+      fTempStruct[myTracks].fStatus = status;
+      fTempStruct[myTracks].fFiltermap = aodtrack->GetFilterMap();
+      fTempStruct[myTracks].fId = aodtrack->GetID();
+      fTempStruct[myTracks].fLabel = aodtrack->GetLabel();
+      fTempStruct[myTracks].fPhi = atan2(fTempStruct[myTracks].fP[1], fTempStruct[myTracks].fP[0]);
+      if(fTempStruct[myTracks].fPhi < 0) fTempStruct[myTracks].fPhi += 2*PI;
+      fTempStruct[myTracks].fPt = sqrt(pow(fTempStruct[myTracks].fP[0],2) + pow(fTempStruct[myTracks].fP[1],2));
+      fTempStruct[myTracks].fMom = sqrt( pow(fTempStruct[myTracks].fPt,2) + pow(fTempStruct[myTracks].fP[2],2) );
+      fTempStruct[myTracks].fEta = aodtrack->Eta();
+      fTempStruct[myTracks].fCharge = aodtrack->Charge();
+      fTempStruct[myTracks].fDCAXY = dca2[0];
+      fTempStruct[myTracks].fDCAZ = dca2[1];
+      fTempStruct[myTracks].fDCA = dca3d;
+      fTempStruct[myTracks].fClusterMap = aodtrack->GetTPCClusterMap();
+      fTempStruct[myTracks].fSharedMap = aodtrack->GetTPCSharedMap();
+      
+    
+      
+      if(fTempStruct[myTracks].fMom > 0.9999) continue;// upper P bound
+      if(fTempStruct[myTracks].fPt > 0.9999) continue;// upper P bound
+      if(fTempStruct[myTracks].fP[2] > 0.9999) continue;// upper P bound
+
+      if(fTempStruct[myTracks].fCharge==+1) {
+       ((TH2F*)fOutputList->FindObject("fDCAxyDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
+       ((TH2F*)fOutputList->FindObject("fDCAzDistPlus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
+      }else {
+       ((TH2F*)fOutputList->FindObject("fDCAxyDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[0]);
+       ((TH2F*)fOutputList->FindObject("fDCAzDistMinus"))->Fill(fTempStruct[myTracks].fPt, dca2[1]);
+      }
+     
+      ((TH3F*)fOutputList->FindObject("fPhiPtDist"))->Fill(aodtrack->Charge(), aodtrack->Phi(), aodtrack->Pt());
+      ((TH3F*)fOutputList->FindObject("fPtEtaDist"))->Fill(aodtrack->Charge(), aodtrack->Pt(), aodtrack->Eta());
+
+            
+      // PID section
+      fTempStruct[myTracks].fElectron = kFALSE;
+      fTempStruct[myTracks].fPion = kFALSE;
+      fTempStruct[myTracks].fKaon = kFALSE;
+      fTempStruct[myTracks].fProton = kFALSE;
+      
+      Float_t nSigmaTPC[5];
+      Float_t nSigmaTOF[5];
+      nSigmaTPC[0]=10; nSigmaTPC[1]=10; nSigmaTPC[2]=10; nSigmaTPC[3]=10; nSigmaTPC[4]=10;
+      nSigmaTOF[0]=10; nSigmaTOF[1]=10; nSigmaTOF[2]=10; nSigmaTOF[3]=10; nSigmaTOF[4]=10;
+      fTempStruct[myTracks].fTOFhit = kFALSE;// default
+      Float_t signalTPC=0, signalTOF=0;
+      Double_t integratedTimesTOF[10]={0};
+
+      
+      if(fFilterBit != 7 || (fMCcase && !fPbPbcase)) {
+       nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron));
+       nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon));
+       nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion));
+       nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon));
+       nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton));
+       //
+       nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron));
+       nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon));
+       nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion));
+       nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon));
+       nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kProton));
+       signalTPC = aodtrack->GetTPCsignal();
+       if( (status&AliESDtrack::kTOFpid)!=0 && (status&AliESDtrack::kTIME)!=0 && (status&AliESDtrack::kTOFout)!=0 && (status&AliESDtrack::kTOFmismatch)<=0){// good tof hit
+         fTempStruct[myTracks].fTOFhit = kTRUE;
+         signalTOF = aodtrack->GetTOFsignal();
+         aodtrack->GetIntegratedTimes(integratedTimesTOF);
+       }else fTempStruct[myTracks].fTOFhit = kFALSE;
+
+      }else {// FilterBit 7 PID workaround
+    
+       for(Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
+         AliAODTrack* aodTrack2 = fAOD->GetTrack(j);
+         if (!aodTrack2) continue;
+         if(aodtrack->GetID() != (-aodTrack2->GetID() - 1)) continue;// (-aodTrack2->GetID() - 1)
+         
+         UInt_t status2=aodTrack2->GetStatus();
+         
+         nSigmaTPC[0]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron));
+         nSigmaTPC[1]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon));
+         nSigmaTPC[2]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion));
+         nSigmaTPC[3]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon));
+         nSigmaTPC[4]=fabs(fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton));
+         //
+         nSigmaTOF[0]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron));
+         nSigmaTOF[1]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon));
+         nSigmaTOF[2]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion));
+         nSigmaTOF[3]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon));
+         nSigmaTOF[4]=fabs(fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kProton));
+         signalTPC = aodTrack2->GetTPCsignal();
+         
+         if( (status2&AliESDtrack::kTOFpid)!=0 && (status2&AliESDtrack::kTIME)!=0 && (status2&AliESDtrack::kTOFout)!=0 && (status2&AliESDtrack::kTOFmismatch)<=0){// good tof hit
+           fTempStruct[myTracks].fTOFhit = kTRUE;
+           signalTOF = aodTrack2->GetTOFsignal();
+           aodTrack2->GetIntegratedTimes(integratedTimesTOF);
+         }else fTempStruct[myTracks].fTOFhit = kFALSE;
+         
+       }// aodTrack2
+      }// FilterBit 7 PID workaround
+
+     
+      ///////////////////
+      ((TH3F*)fOutputList->FindObject("fTPCResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTPC);
+      if(fTempStruct[myTracks].fTOFhit) {
+       ((TH3F*)fOutputList->FindObject("fTOFResponse"))->Fill(centralityPercentile, fTempStruct[myTracks].fMom, signalTOF - integratedTimesTOF[3]);
+      }
+      ///////////////////
+
+      // Use TOF if good hit and above threshold
+      if(fTempStruct[myTracks].fTOFhit && fTempStruct[myTracks].fMom > fTPCTOFboundry){
+       if(nSigmaTOF[0]<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+       if(nSigmaTOF[2]<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+       if(nSigmaTOF[3]<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+       if(nSigmaTOF[4]<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+      }else {// TPC info instead
+       if(nSigmaTPC[0]<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+       if(nSigmaTPC[2]<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+       if(nSigmaTPC[3]<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+       if(nSigmaTPC[4]<fSigmaCutTPC) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+      }
+               
+      
+      // Ensure there is only 1 candidate per track
+      if(fTempStruct[myTracks].fElectron && fTempStruct[myTracks].fMom < 0.45) continue;// Remove electron band
+      if(!fTempStruct[myTracks].fPion && !fTempStruct[myTracks].fKaon && !fTempStruct[myTracks].fProton) continue;
+      if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon) continue;
+      if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fProton) continue;
+      if(fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
+      if(fTempStruct[myTracks].fPion && fTempStruct[myTracks].fKaon && fTempStruct[myTracks].fProton) continue;
+      ////////////////////////
+      if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons
+
+      if(!fTempStruct[myTracks].fPion) continue;// only pions
+      
+          
+           
+    
+      if(fTempStruct[myTracks].fPion) {// pions
+       fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2)); 
+       fTempStruct[myTracks].fKey = 1;
+      }else if(fTempStruct[myTracks].fKaon){// kaons
+       fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassK,2));;
+       fTempStruct[myTracks].fKey = 10;
+      }else{// protons
+       fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassP,2));;
+       fTempStruct[myTracks].fKey = 100;
+      }
+      
+     
+      ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
+      ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
+      if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
+      if(aodtrack->GetTPCncls() < fMinTPCncls) continue;
+      
+
+      if(aodtrack->Charge() > 0) positiveTracks++;
+      else negativeTracks++;
+      
+      if(fTempStruct[myTracks].fPion) pionCount++;
+      if(fTempStruct[myTracks].fKaon) kaonCount++;
+      if(fTempStruct[myTracks].fProton) protonCount++;
+
+      myTracks++;
+      
+    }
+  }else {// ESD tracks
+    cout<<"ESDs not supported currently"<<endl;
+    return;
+  }
+  
+  
+  if(myTracks >= 1) {
+    ((TH1F*)fOutputList->FindObject("fMultDist3"))->Fill(myTracks);
+  }
+  //cout<<"There are "<<myTracks<<"  myTracks"<<endl;
+  //cout<<"pionCount = "<<pionCount<<"   kaonCount = "<<kaonCount<<"   protonCount = "<<protonCount<<endl;
+
+  /////////////////////////////////////////
+  // Pion Multiplicity Cut (To ensure all Correlation orders are present in each event)
+  if(myTracks < 3) {cout<<"Less than 3 tracks. Skipping Event."<<endl; return;}
+  /////////////////////////////////////////
+
+  ////////////////////////////////
+  ///////////////////////////////
+  // Mbin determination
+  //
+  fMbin=-1;
+  for(Int_t i=0; i<fMbins; i++){
+    if(fPbPbcase){
+      if(centralityPercentile >= 5*i && centralityPercentile < 5*(i+1)) {fMbin=i; break;}
+    }else{
+     if( (pow(trackletMult,1/3.) >= fMultLimits[i]) && (pow(trackletMult,1/3.) < fMultLimits[i+1])) {fMbin=i; break;}
+    } 
+    //if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}// Mbin 0 has 0-5 pions
+  }
+
+    
+
+  if(fMbin==-1) {cout<<pionCount<<"  Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
+  if(fMbin > fCentBinHighLimit) {cout<<"Mult out of range"<<endl; return;}
+  
+
+  fFSIindex=0;
+  if(fPbPbcase){
+    if(fMbin==0) fFSIindex = 0;//0-5%
+    else if(fMbin==1) fFSIindex = 1;//5-10%
+    else if(fMbin<=3) fFSIindex = 2;//10-20%
+    else if(fMbin<=5) fFSIindex = 3;//20-30%
+    else if(fMbin<=7) fFSIindex = 4;//30-40%
+    else if(fMbin<=9) fFSIindex = 5;//40-50%
+    else if(fMbin<=12) fFSIindex = 6;//40-50%
+    else if(fMbin<=15) fFSIindex = 7;//40-50%
+    else if(fMbin<=18) fFSIindex = 8;//40-50%
+    else fFSIindex = 8;//90-100%
+  }else fFSIindex = 9;// pp and pPb
+  
+  //////////////////////////////////////////////////
+  fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
+  //////////////////////////////////////////////////
+  
+  
+  
+  ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
+  ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
+  if(fMCcase){
+    ((TH2D*)fOutputList->FindObject("fdNchdEtaResponse"))->Fill(pow(trackletMult,1/3.), pow(mcdNchdEta,1/3.));
+  }
+  if(fPbPbcase){
+    ((TH2D*)fOutputList->FindObject("fdCentVsNchdEta"))->Fill(fMbin+1, pow(trackletMult,1/3.));
+  }
+  //cout<<trackletMult<<"  "<<mcdNchdEta<<endl;
+  
+  ////////////////////////////////////
+  // Add event to buffer if > 0 tracks
+  if(myTracks > 0){
+    fEC[zbin][fMbin]->FIFOShift();
+    (fEvt) = fEC[zbin][fMbin]->fEvtStr;
+    (fEvt)->fNtracks = myTracks;
+    (fEvt)->fFillStatus = 1;
+    for(Int_t i=0; i<myTracks; i++) (fEvt)->fTracks[i] = fTempStruct[i];
+    if(fMCcase){
+      (fEvt)->fMCarraySize = mcArray->GetEntriesFast();
+      for(Int_t i=0; i<mcArray->GetEntriesFast(); i++) {
+       AliAODMCParticle *tempMCTrack = (AliAODMCParticle*)mcArray->At(i);
+       (fEvt)->fMCtracks[i].fPx = tempMCTrack->Px();
+       (fEvt)->fMCtracks[i].fPy = tempMCTrack->Py();
+       (fEvt)->fMCtracks[i].fPz = tempMCTrack->Pz();
+       (fEvt)->fMCtracks[i].fPtot = sqrt(pow(tempMCTrack->Px(),2)+pow(tempMCTrack->Py(),2)+pow(tempMCTrack->Pz(),2));
+      }        
+    }
+  }
+    
+  
+  
+  Float_t qinv12=0, qinv13=0, qinv23=0;
+  Float_t qout=0, qside=0, qlong=0;
+  Float_t qoutMC=0, qsideMC=0, qlongMC=0;
+  Float_t firstQ=0, secondQ=0, thirdQ=0;
+  Float_t firstQMC=0, secondQMC=0, thirdQMC=0;
+  Float_t transK12=0, transK3=0;
+  Float_t q3=0, q3MC=0;
+  Int_t ch1=0, ch2=0, ch3=0;
+  Short_t key1=0, key2=0, key3=0;
+  Int_t bin1=0, bin2=0, bin3=0;
+  Float_t pVect1[4]={0}; 
+  Float_t pVect2[4]={0};
+  Float_t pVect3[4]={0}; 
+  Float_t pVect1MC[4]={0}; 
+  Float_t pVect2MC[4]={0};
+  Float_t pVect3MC[4]={0};
+  Int_t index1=0, index2=0, index3=0;
+  Float_t qinv12MC=0, qinv13MC=0, qinv23MC=0;
+  //
+  AliAODMCParticle *mcParticle1=0x0;
+  AliAODMCParticle *mcParticle2=0x0;
+  
+
+  if(fPdensityPairCut){
+    ////////////////////
+    Int_t pairCountSE=0, pairCountME=0;
+    Int_t normPairCount[2]={0};
+    Int_t numOtherPairs1[2][fMultLimit];
+    Int_t numOtherPairs2[2][fMultLimit];
+    Bool_t exitCode=kFALSE;
+    Int_t tempNormFillCount[2][2][2][10][5];
+    
+
+    // reset to defaults
+    for(Int_t i=0; i<fMultLimit; i++) {
+      fPairLocationSE[i]->Set(fMultLimit,fDefaultsInt);
+      fPairLocationME[i]->Set(fMultLimit,fDefaultsInt);
+           
+      // Normalization Utilities
+      fOtherPairLocation1[0][i]->Set(fMultLimit,fDefaultsInt);
+      fOtherPairLocation1[1][i]->Set(fMultLimit,fDefaultsInt);
+      fOtherPairLocation2[0][i]->Set(fMultLimit,fDefaultsInt);
+      fOtherPairLocation2[1][i]->Set(fMultLimit,fDefaultsInt);
+      fNormPairSwitch[0][i]->Set(fMultLimit,fDefaultsCharMult);
+      fNormPairSwitch[1][i]->Set(fMultLimit,fDefaultsCharMult);
+      fNormPairSwitch[2][i]->Set(fMultLimit,fDefaultsCharMult);
+      numOtherPairs1[0][i]=0;
+      numOtherPairs1[1][i]=0;
+      numOtherPairs2[0][i]=0;
+      numOtherPairs2[1][i]=0;
+      
+      // Track Merging/Splitting Utilities
+      fPairSplitCut[0][i]->Set(fMultLimit,fDefaultsCharMult);// P11
+      fPairSplitCut[1][i]->Set(fMultLimit,fDefaultsCharMult);// P12
+      fPairSplitCut[2][i]->Set(fMultLimit,fDefaultsCharMult);// P13
+      fPairSplitCut[3][i]->Set(fMultLimit,fDefaultsCharMult);// P23
+    }
+
+    // Reset the temp Normalization counters
+    for(Int_t i=0; i<2; i++){// Charge1
+      for(Int_t j=0; j<2; j++){// Charge2
+       for(Int_t k=0; k<2; k++){// Charge3
+         for(Int_t l=0; l<10; l++){// FillIndex (species Combination)
+           for(Int_t m=0; m<5; m++){// Term (Cumulant term)
+             tempNormFillCount[i][j][k][l][m] = 0;
+           }
+         }
+       }
+      }
+    }
+             
+    ///////////////////////////////////////////////////////
+    // Start the pairing process
+    // P11 pairing
+    // 1st Particle
+   
+    for (Int_t i=0; i<myTracks; i++) {
+         
+      Int_t en2=0;
+   
+      // 2nd particle
+      for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {
+       
+       key1 = (fEvt)->fTracks[i].fKey;
+       key2 = (fEvt+en2)->fTracks[j].fKey;
+       Short_t fillIndex2 = FillIndex2part(key1+key2);
+       Short_t qCutBin = SetQcutBin(fillIndex2);
+       Short_t normBin = SetNormBin(fillIndex2);
+       pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
+       pVect1[1]=(fEvt)->fTracks[i].fP[0];      pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
+       pVect1[2]=(fEvt)->fTracks[i].fP[1];      pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
+       pVect1[3]=(fEvt)->fTracks[i].fP[2];      pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
+       
+       //
+       
+       qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
+       GetQosl(pVect1, pVect2, qout, qside, qlong);
+       transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
+
+
+       //
+
+       ///////////////////////////////
+       ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
+       ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
+       SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
+       
+       if(fMCcase && ch1==ch2 && fMbin==0 && qinv12<0.2){
+         //////////////////////////
+         // pad-row method testing
+         Float_t coeff = (5)*0.2*(0.18/1.2);// 5 to evaluate at 1.0m in TPC
+         Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
+         if(phi1 > 2*PI) phi1 -= 2*PI;
+         if(phi1 < 0) phi1 += 2*PI;
+         Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
+         if(phi2 > 2*PI) phi2 -= 2*PI;
+         if(phi2 < 0) phi2 += 2*PI;
+         Float_t deltaphi = phi1 - phi2;
+         if(deltaphi > PI) deltaphi -= PI;
+         if(deltaphi < -PI) deltaphi += PI;
+         
+         Int_t ncl1 = (fEvt)->fTracks[i].fClusterMap.GetNbits();
+         Int_t ncl2 = (fEvt+en2)->fTracks[j].fClusterMap.GetNbits();
+         Float_t sumCls = 0; Float_t sumSha = 0; Float_t sumQ = 0;
+         Double_t shfrac = 0; //Double_t qfactor = 0;
+         for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
+           if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Both clusters
+             if ((fEvt)->fTracks[i].fSharedMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fSharedMap.TestBitNumber(imap)) { // Shared
+               sumQ++;
+               sumCls+=2;
+               sumSha+=2;}
+             else {sumQ--; sumCls+=2;}
+           }
+           else if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) || (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Non shared
+             sumQ++;
+             sumCls++;}
+         }
+         if (sumCls>0) {
+           //qfactor = sumQ*1.0/sumCls;
+           shfrac = sumSha*1.0/sumCls;
+         }
+         if(fabs(deltaphi)<0.07 && fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
+           ((TH3D*)fOutputList->FindObject("fPairsPadRowNum"))->Fill(transK12, shfrac, qinv12);
+         }
+         
+         for(Int_t rstep=0; rstep<10; rstep++){
+           coeff = (rstep)*0.2*(0.18/1.2);
+           phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
+           if(phi1 > 2*PI) phi1 -= 2*PI;
+           if(phi1 < 0) phi1 += 2*PI;
+           phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
+           if(phi2 > 2*PI) phi2 -= 2*PI;
+           if(phi2 < 0) phi2 += 2*PI;
+           deltaphi = phi1 - phi2;
+           if(deltaphi > PI) deltaphi -= PI;
+           if(deltaphi < -PI) deltaphi += PI;
+
+           if(fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
+             ((TH3F*)fOutputList->FindObject("fPairsShareFracDPhiNum"))->Fill(rstep, shfrac, deltaphi);
+           }
+           //if(shfrac < 0.05){
+           ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
+           //}
+         }
+         
+         
+       }// MCcase and pair selection
+       
+       // Pair Splitting/Merging cut
+       if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
+       if(ch1 == ch2){
+         if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
+           fPairSplitCut[0][i]->AddAt('1',j);
+           ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
+           continue;
+         }
+       }
+       
+       // HIJING tests 
+       if(fMCcase && fillIndex2==0){
+         
+         // Check that label does not exceed stack size
+         if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
+           pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
+           pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
+           pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
+           pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
+           pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
+           qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
+           GetQosl(pVect1MC, pVect2MC, qoutMC, qsideMC, qlongMC);
+           if(qinv12<0.1 && ch1==ch2) {
+             ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC); 
+             ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
+             ((TH1D*)fOutputList->FindObject("fQDist"))->Fill(qinv12-qinv12MC);
+           }
+           
+           
+          
+           mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
+           mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
+           
+          
+           // secondary contamination
+           if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
+             if(ch1==ch2) {
+               ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
+               if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
+                 ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
+               }             
+             }else{
+               ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
+               if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
+                 ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, transK12, qinv12);
+               }
+             }
+           }
+
+           Float_t rForQW=5.0;
+           if(fFSIindex<=1) rForQW=10;
+           else if(fFSIindex==2) rForQW=9;
+           else if(fFSIindex==3) rForQW=8;
+           else if(fFSIindex==4) rForQW=7;
+           else if(fFSIindex==5) rForQW=6;
+           else if(fFSIindex==6) rForQW=5;
+           else if(fFSIindex==7) rForQW=4;
+           else if(fFSIindex==8) rForQW=3;
+           else rForQW=2;
+           
+           Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fMCqinv->Fill(qinv12MC, MCWeight(ch1,ch2,rForQW,10,qinv12MC));// was 4,5
+           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
+           // pion purity          
+           Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityDen->Fill(transK12, qinv12);
+           if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){// Pions
+             Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fPIDpurityNum->Fill(transK12, qinv12);
+           }
+
+         }// label check 2
+       }// MC case
+       
+       //////////////////////////////////////////
+       // 2-particle term
+       Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
+       Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
+       Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt)->fTracks[i].fMom);
+       Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fAvgP->Fill(transK12, qinv12, (fEvt+en2)->fTracks[j].fMom);
+       
+       
+       //////////////////////////////////////////
+       
+       
+
+       // exit out of loop if there are too many pairs  
+       if(pairCountSE >= kPairLimit) {exitCode=kTRUE; continue;}// Too many SE pairs
+       if(exitCode) continue;
+
+       //////////////////////////
+       // Enforce the Qcut
+       if(qinv12 <= fQcut[qCutBin]) {
+        
+         ///////////////////////////
+         // particle 1
+         (fEvt)->fPairsSE[pairCountSE].fP1[0] = (fEvt)->fTracks[i].fP[0];
+         (fEvt)->fPairsSE[pairCountSE].fP1[1] = (fEvt)->fTracks[i].fP[1];
+         (fEvt)->fPairsSE[pairCountSE].fP1[2] = (fEvt)->fTracks[i].fP[2];
+         (fEvt)->fPairsSE[pairCountSE].fE1 = (fEvt)->fTracks[i].fEaccepted;
+         (fEvt)->fPairsSE[pairCountSE].fCharge1 = (fEvt)->fTracks[i].fCharge;
+         (fEvt)->fPairsSE[pairCountSE].fIndex1 = i;
+         (fEvt)->fPairsSE[pairCountSE].fKey1 = key1;
+         (fEvt)->fPairsSE[pairCountSE].fLabel1 = (fEvt)->fTracks[i].fLabel;
+         if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
+           (fEvt)->fPairsSE[pairCountSE].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
+           (fEvt)->fPairsSE[pairCountSE].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
+           (fEvt)->fPairsSE[pairCountSE].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
+         }
+         // particle 2
+         (fEvt)->fPairsSE[pairCountSE].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
+         (fEvt)->fPairsSE[pairCountSE].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
+         (fEvt)->fPairsSE[pairCountSE].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
+         (fEvt)->fPairsSE[pairCountSE].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
+         (fEvt)->fPairsSE[pairCountSE].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
+         (fEvt)->fPairsSE[pairCountSE].fIndex2 = j;
+         (fEvt)->fPairsSE[pairCountSE].fKey2 = key2;
+         (fEvt)->fPairsSE[pairCountSE].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
+         if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
+           (fEvt)->fPairsSE[pairCountSE].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
+           (fEvt)->fPairsSE[pairCountSE].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
+           (fEvt)->fPairsSE[pairCountSE].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
+         }
+       
+         (fEvt)->fPairsSE[pairCountSE].fQinv = qinv12;
+         
+         fPairLocationSE[i]->AddAt(pairCountSE,j);
+         
+         pairCountSE++;
+         
+       }
+
+       
+       /////////////////////////////////////////////////////////
+       // Normalization Region
+       
+       if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
+         // particle 1
+         fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
+         fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
+         fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
+         // particle 2
+         fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
+         fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
+         fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
+         
+         
+         //other past pairs with particle j
+         for(Int_t pastpair=0; pastpair<numOtherPairs2[0][j]; pastpair++){
+           Int_t locationOtherPair = fOtherPairLocation2[0][j]->At(pastpair);
+           if(locationOtherPair < 0) continue;// no pair there
+           Int_t indexOther1 = i;
+           Int_t indexOther2 = fNormPairs[0][ locationOtherPair ].fIndex1;
+           
+           // Both possible orderings of other indexes
+           if( (fNormPairSwitch[0][indexOther1]->At(indexOther2)=='1') || (fNormPairSwitch[0][indexOther2]->At(indexOther1)=='1')) {
+             
+             // 1 and 2 are from SE
+             ch3 = Int_t((fNormPairs[0][ locationOtherPair ].fCharge1 + 1)/2.);
+             key3 = fNormPairs[0][ locationOtherPair ].fKey1;
+             Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
+             SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
+             
+             tempNormFillCount[bin1][bin2][bin3][fillIndex3][0]++;
+           }
+           
+         }// pastpair P11 loop
+         
+         
+         fNormPairSwitch[en2][i]->AddAt('1',j);            
+         fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
+         fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
+         
+         numOtherPairs1[en2][i]++;
+         numOtherPairs2[en2][j]++;
+         
+         
+         normPairCount[en2]++;
+         if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
+         
+       }// Norm Region
+       
+      }// j particle
+    }// i particle
+    
+
+    
+    //////////////////////////////////////////////
+    // P12 pairing
+    // 1st Particle
+    for (Int_t i=0; i<myTracks; i++) {
+         
+      Int_t en2=1;
+
+      // 2nd particle
+      for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
+               
+       key1 = (fEvt)->fTracks[i].fKey;
+       key2 = (fEvt+en2)->fTracks[j].fKey;
+       Short_t fillIndex2 = FillIndex2part(key1+key2);
+       Short_t qCutBin = SetQcutBin(fillIndex2);
+       Short_t normBin = SetNormBin(fillIndex2);
+       pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
+       pVect1[1]=(fEvt)->fTracks[i].fP[0];      pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
+       pVect1[2]=(fEvt)->fTracks[i].fP[1];      pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
+       pVect1[3]=(fEvt)->fTracks[i].fP[2];      pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
+       
+       qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
+       GetQosl(pVect1, pVect2, qout, qside, qlong);
+       transK12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
+       //if(transK12 <= 0.35) fEDbin=0;
+       //else fEDbin=1;
+
+       
+       
+       ///////////////////////////////
+       ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
+       ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
+       SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
+       
+       if(fMCcase){
+         if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
+           pVect1MC[0]=sqrt(pow((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
+           pVect2MC[0]=sqrt(pow((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPtot,2)+pow(fTrueMassPi,2));
+           pVect1MC[1]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx; pVect2MC[1]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
+           pVect1MC[2]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy; pVect2MC[2]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
+           pVect1MC[3]=(fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz; pVect2MC[3]=(fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
+           qinv12MC = GetQinv(fillIndex2, pVect1MC, pVect2MC);
+           //
+         
+           for(Int_t rIter=0; rIter<fRVALUES; rIter++){
+             for(Int_t myDampIt=0; myDampIt<kNDampValues; myDampIt++){
+               Int_t denIndex = rIter*kNDampValues + myDampIt;
+               Float_t WInput = MCWeight(ch1,ch2, rIter+kRmin, myDampIt, qinv12MC);
+               Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fIdeal->Fill(denIndex, qinv12MC, WInput);
+               Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fIdeal->Fill(denIndex, qinv12MC);
+               Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[0].fSmeared->Fill(denIndex, qinv12, WInput);
+               Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[1].fSmeared->Fill(denIndex, qinv12);
+             }
+           }
+         
+           
+           /////////////////////////////////////////////////////
+           if(qinv12 <= fQcut[qCutBin]) {// 3-particle MRC
+             
+             // 3-particle MRC
+             Short_t fillIndex3 = 0;
+             key1=1; key2=1; key3=1;
+             Int_t en3 = 2;
+             
+             for (Int_t k=0; k<(fEvt+en3)->fNtracks; k++) {
+               if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
+                 ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
+                 pVect3[0]=(fEvt+en3)->fTracks[k].fEaccepted;
+                 pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
+                 pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
+                 pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
+                 qinv13 = GetQinv(0, pVect1, pVect3);
+                 qinv23 = GetQinv(0, pVect2, pVect3);
+                 q3 = sqrt(pow(qinv12,2)+pow(qinv13,2)+pow(qinv23,2));
+                 
+                 if(qinv13 > fQcut[qCutBin] || qinv23 > fQcut[qCutBin]) continue;
+                 
+                 
+                 pVect3MC[0]=sqrt(pow((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
+                 pVect3MC[1]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPx;
+                 pVect3MC[2]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPy;
+                 pVect3MC[3]=(fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPz;
+                 qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
+                 qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
+                 
+                 
+                 q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
+                 transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
+                 
+
+                 //
+                 // 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.
+                 Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
+                 SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
+                 
+                 
+                 for(Int_t jj=1; jj<=5; jj++){// term loop
+                   
+                   if(jj==2) {if(!fill2) continue;}//12
+                   if(jj==3) {if(!fill3) continue;}//13
+                   if(jj==4) {if(!fill4) continue;}//23
+                   
+                   Float_t WInput=1.0;
+                   ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 1, jj, firstQ, secondQ, thirdQ);
+                   ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, 1, jj, firstQMC, secondQMC, thirdQMC);
+                   
+                   if(ch1==ch2 && ch1==ch3){// same charge
+                     WInput = MCWeight3D(kTRUE, jj, 9, firstQMC, secondQMC, thirdQMC);
+                     if(jj==1) {
+                       ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1SC"))->Fill(q3MC, WInput);
+                       ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1SCden"))->Fill(q3MC);
+                     }else if(jj!=5){
+                       ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2SC"))->Fill(q3MC, WInput);
+                       ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2SCden"))->Fill(q3MC);
+                     }
+                   }else {// mixed charge
+                     if(bin1==bin2) {
+                       WInput = MCWeight3D(kFALSE, jj, 9, firstQMC, secondQMC, thirdQMC);
+                     }else {
+                       if(jj==1 || jj==5) WInput = MCWeight3D(kFALSE, jj, 9, thirdQMC, secondQMC, firstQMC);// thirdQMC is ss
+                       else WInput = MCWeight3D(kFALSE, 6-jj, 9, thirdQMC, secondQMC, firstQMC);
+                     }
+                     if(jj==1){
+                       ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1MC"))->Fill(q3MC, WInput);
+                       ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm1MCden"))->Fill(q3MC);
+                     }else{
+                       if(bin1==bin2){
+                         if(jj==2){
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MC"))->Fill(q3MC, WInput);
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MCden"))->Fill(q3MC);
+                         }else if(jj==3){
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MC"))->Fill(q3MC, WInput);
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MCden"))->Fill(q3MC);
+                         }else if(jj==4){
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MC"))->Fill(q3MC, WInput);
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MCden"))->Fill(q3MC);
+                         }else{}
+                       }else{
+                         if(jj==2){
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MC"))->Fill(q3MC, WInput);
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm4MCden"))->Fill(q3MC);
+                         }else if(jj==3){
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MC"))->Fill(q3MC, WInput);
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm3MCden"))->Fill(q3MC);
+                         }else if(jj==4){
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MC"))->Fill(q3MC, WInput);
+                           ((TH1D*)fOutputList->FindObject("fMCWeight3DTerm2MCden"))->Fill(q3MC);
+                         }else{}
+                       }
+                       
+                     }
+                   }
+                   
+                   
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fIdeal->Fill(q3MC, WInput);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fSmeared->Fill(q3, WInput);
+                   
+                                   
+                 }// jj
+               }// MCarray check, 3rd particle
+             }// 3rd particle
+             
+           }// TabulatePairs Check
+           
+         }// MCarray check, 1st and 2nd particle
+         
+         // reset key's and fill bins (they were altered for 3 particle MRC calculation)
+         key1 = (fEvt)->fTracks[i].fKey;
+         key2 = (fEvt+en2)->fTracks[j].fKey;
+         SetFillBins2(fillIndex2, key1, key2, ch1, ch2, bin1, bin2);
+
+         
+         if(ch1==ch2 && fMbin==0 && qinv12<0.2){
+           //////////////////////////
+           // pad-row method testing
+           Float_t coeff = (5)*0.2*(0.18/1.2);// 5 to evaluate at 1.0m in TPC
+           Float_t phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
+           if(phi1 > 2*PI) phi1 -= 2*PI;
+           if(phi1 < 0) phi1 += 2*PI;
+           Float_t phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
+           if(phi2 > 2*PI) phi2 -= 2*PI;
+           if(phi2 < 0) phi2 += 2*PI;
+           Float_t deltaphi = phi1 - phi2;
+           if(deltaphi > PI) deltaphi -= PI;
+           if(deltaphi < -PI) deltaphi += PI;
+           
+           Int_t ncl1 = (fEvt)->fTracks[i].fClusterMap.GetNbits();
+           Int_t ncl2 = (fEvt+en2)->fTracks[j].fClusterMap.GetNbits();
+           Float_t sumCls = 0; Float_t sumSha = 0; Float_t sumQ = 0;
+           Double_t shfrac = 0; //Double_t qfactor = 0;
+           for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
+             if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Both clusters
+               if ((fEvt)->fTracks[i].fSharedMap.TestBitNumber(imap) && (fEvt+en2)->fTracks[j].fSharedMap.TestBitNumber(imap)) { // Shared
+                 sumQ++;
+                 sumCls+=2;
+                 sumSha+=2;}
+               else {sumQ--; sumCls+=2;}
+             }
+             else if ((fEvt)->fTracks[i].fClusterMap.TestBitNumber(imap) || (fEvt+en2)->fTracks[j].fClusterMap.TestBitNumber(imap)) {// Non shared
+               sumQ++;
+               sumCls++;}
+           }
+           if (sumCls>0) {
+             //qfactor = sumQ*1.0/sumCls;
+             shfrac = sumSha*1.0/sumCls;
+           }
+           if(fabs(deltaphi)<0.07 && fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
+             ((TH3D*)fOutputList->FindObject("fPairsPadRowDen"))->Fill(transK12, shfrac, qinv12);
+           }
+           
+           for(Int_t rstep=0; rstep<10; rstep++){
+             coeff = (rstep)*0.2*(0.18/1.2);
+             // propagate through B field to r=1.2m
+             phi1 = (fEvt)->fTracks[i].fPhi - asin((fEvt)->fTracks[i].fCharge*(0.1*fBfield)*coeff/(fEvt)->fTracks[i].fPt);
+             if(phi1 > 2*PI) phi1 -= 2*PI;
+             if(phi1 < 0) phi1 += 2*PI;
+             phi2 = (fEvt+en2)->fTracks[j].fPhi - asin((fEvt+en2)->fTracks[j].fCharge*(0.1*fBfield)*coeff/(fEvt+en2)->fTracks[j].fPt);
+             if(phi2 > 2*PI) phi2 -= 2*PI;
+             if(phi2 < 0) phi2 += 2*PI;
+             deltaphi = phi1 - phi2;
+             if(deltaphi > PI) deltaphi -= PI;
+             if(deltaphi < -PI) deltaphi += PI;
+             
+             if(fabs((fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta) < 0.03){
+               ((TH3F*)fOutputList->FindObject("fPairsShareFracDPhiDen"))->Fill(rstep, shfrac, deltaphi);
+             }
+             //if(shfrac < 0.05){
+             ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
+             //}
+           }
+           
+          
+           
+           
+         }// desired pair selection
+         
+       
+  
+       }// fMCcase
+       
+       
+
+       if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
+       if(ch1 == ch2){
+         if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
+           fPairSplitCut[1][i]->AddAt('1',j);
+           continue;
+         }
+       }
+       
+       //////////////////////////////////////////
+       // 2-particle term
+       Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2->Fill(transK12, qinv12);
+       Charge1[bin1].Charge2[bin2].SC[fillIndex2].MB[fMbin].EDB[fEDbin].TwoPT[en2].fExplicit2QW->Fill(transK12, qinv12, qinv12);
+       
+       
+       //////////////////////////////////////////
+
+       
+       
+       if(pairCountME >= 2*kPairLimit) {exitCode=kTRUE; continue;}// Too many SE pairs
+       if(exitCode) continue;
+
+       if(qinv12 <= fQcut[qCutBin]) {
+         ///////////////////////////
+         
+         // particle 1
+         (fEvt)->fPairsME[pairCountME].fP1[0] = (fEvt)->fTracks[i].fP[0];
+         (fEvt)->fPairsME[pairCountME].fP1[1] = (fEvt)->fTracks[i].fP[1];
+         (fEvt)->fPairsME[pairCountME].fP1[2] = (fEvt)->fTracks[i].fP[2];
+         (fEvt)->fPairsME[pairCountME].fE1 = (fEvt)->fTracks[i].fEaccepted;
+         (fEvt)->fPairsME[pairCountME].fCharge1 = (fEvt)->fTracks[i].fCharge;
+         (fEvt)->fPairsME[pairCountME].fIndex1 = i;
+         (fEvt)->fPairsME[pairCountME].fKey1 = key1;
+         (fEvt)->fPairsME[pairCountME].fLabel1 = (fEvt)->fTracks[i].fLabel;
+         if(fMCcase && ((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize)){
+           (fEvt)->fPairsME[pairCountME].fP1MC[0] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPx;
+           (fEvt)->fPairsME[pairCountME].fP1MC[1] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPy;
+           (fEvt)->fPairsME[pairCountME].fP1MC[2] = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPz;
+         }
+         // particle 2
+         (fEvt)->fPairsME[pairCountME].fP2[0] = (fEvt+en2)->fTracks[j].fP[0];
+         (fEvt)->fPairsME[pairCountME].fP2[1] = (fEvt+en2)->fTracks[j].fP[1];
+         (fEvt)->fPairsME[pairCountME].fP2[2] = (fEvt+en2)->fTracks[j].fP[2];
+         (fEvt)->fPairsME[pairCountME].fE2 = (fEvt+en2)->fTracks[j].fEaccepted;
+         (fEvt)->fPairsME[pairCountME].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
+         (fEvt)->fPairsME[pairCountME].fIndex2 = j;
+         (fEvt)->fPairsME[pairCountME].fKey2 = key2;
+         (fEvt)->fPairsME[pairCountME].fLabel2 = (fEvt+en2)->fTracks[j].fLabel;
+         if(fMCcase && ((fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize)){
+           (fEvt)->fPairsME[pairCountME].fP2MC[0] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPx;
+           (fEvt)->fPairsME[pairCountME].fP2MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPy;
+           (fEvt)->fPairsME[pairCountME].fP2MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPz;
+         }
+         
+         (fEvt)->fPairsME[pairCountME].fQinv = qinv12;
+         
+         fPairLocationME[i]->AddAt(Int_t(pairCountME),j);
+         
+         pairCountME++;
+         
+       }
+       
+       if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
+         // particle 1
+         fNormPairs[en2][normPairCount[en2]].fCharge1 = (fEvt)->fTracks[i].fCharge;
+         fNormPairs[en2][normPairCount[en2]].fIndex1 = i;
+         fNormPairs[en2][normPairCount[en2]].fKey1 = (fEvt)->fTracks[i].fKey;
+         // particle 2
+         fNormPairs[en2][normPairCount[en2]].fCharge2 = (fEvt+en2)->fTracks[j].fCharge;
+         fNormPairs[en2][normPairCount[en2]].fIndex2 = j;
+         fNormPairs[en2][normPairCount[en2]].fKey2 = (fEvt+en2)->fTracks[j].fKey;
+         
+         //other past pairs in P11 with particle i
+         for(Int_t pastpairP11=0; pastpairP11<numOtherPairs2[0][i]; pastpairP11++){// past pair in P11 with i as 1st and 2nd particle
+           Int_t locationOtherPairP11 = fOtherPairLocation2[0][i]->At(pastpairP11);// i is 2nd particle
+           if(locationOtherPairP11 < 0) continue;// no pair there
+           Int_t indexOther1P11 = fNormPairs[0][ locationOtherPairP11 ].fIndex1; 
+                   
+           //Check other past pairs in P12
+           if( (fNormPairSwitch[1][indexOther1P11]->At(j)=='0')) continue;
+           
+           // 1 and 3 are from SE
+           ch3 = Int_t((fNormPairs[0][ locationOtherPairP11 ].fCharge1 + 1)/2.);// charge of second particle in P11
+           key3 = fNormPairs[0][ locationOtherPairP11 ].fKey1;
+           Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
+           Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
+           SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 2, bin1, bin2, bin3, fill2, fill3, fill4);
+           
+                   
+           if(fill2) tempNormFillCount[bin1][bin2][bin3][fillIndex3][1]++;
+           if(fill3) tempNormFillCount[bin1][bin2][bin3][fillIndex3][2]++;
+           if(fill4) tempNormFillCount[bin1][bin2][bin3][fillIndex3][3]++;
+           
+
+         }// P11 loop
+         
+         
+         fNormPairSwitch[en2][i]->AddAt('1',j);            
+         fOtherPairLocation1[en2][i]->AddAt(normPairCount[en2], numOtherPairs1[en2][i]);// location of otherpair with i as 1st particle
+         fOtherPairLocation2[en2][j]->AddAt(normPairCount[en2], numOtherPairs2[en2][j]);// location of otherpair with j as 2nd particle
+         
+         numOtherPairs1[en2][i]++;
+         numOtherPairs2[en2][j]++;
+         
+         normPairCount[en2]++;
+         if(normPairCount[en2] >= kNormPairLimit) exitCode=kTRUE;
+
+       }// Norm Region
+       
+
+      }
+    }
+    
+    ///////////////////////////////////////
+    // P13 pairing (just for Norm counting of term5)
+    for (Int_t i=0; i<myTracks; i++) {
+      
+      // exit out of loop if there are too many pairs
+      // dont bother with this loop if exitCode is set.
+      if(exitCode) break;
+      
+      // 2nd particle
+      Int_t en2=2;
+      
+      for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
+       
+       key1 = (fEvt)->fTracks[i].fKey;
+       key2 = (fEvt+en2)->fTracks[j].fKey;
+       Short_t fillIndex2 = FillIndex2part(key1+key2);
+       Short_t normBin = SetNormBin(fillIndex2);
+       pVect1[0]=(fEvt)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
+       pVect1[1]=(fEvt)->fTracks[i].fP[0];      pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
+       pVect1[2]=(fEvt)->fTracks[i].fP[1];      pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
+       pVect1[3]=(fEvt)->fTracks[i].fP[2];      pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
+
+       qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
+       
+       if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
+       
+       ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
+       ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
+       
+       if(ch1 == ch2){
+         if(!AcceptPair((fEvt)->fTracks[i], (fEvt+en2)->fTracks[j])) {
+           fPairSplitCut[2][i]->AddAt('1',j);
+           continue;
+         }
+       }
+       
+       /////////////////////////////////////////////////////////
+       // Normalization Region
+       
+       if((qinv12 >= fNormQcutLow[normBin]) && (qinv12 < fNormQcutHigh[normBin])){
+       
+         fNormPairSwitch[en2][i]->AddAt('1',j);            
+       
+       }// Norm Region
+      }
+    }
+
+
+   
+    ///////////////////////////////////////
+    // P23 pairing (just for Norm counting of term5)
+    Int_t en1=1;
+    for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {
+      
+      // exit out of loop if there are too many pairs
+      // dont bother with this loop if exitCode is set.
+      if(exitCode) break;
+      
+      // 2nd event
+      Int_t en2=2;
+      // 2nd particle
+      for (Int_t j=0; j<(fEvt+en2)->fNtracks; j++) {
+       
+       if(exitCode) break;
+
+       key1 = (fEvt+en1)->fTracks[i].fKey;
+       key2 = (fEvt+en2)->fTracks[j].fKey;
+       Short_t fillIndex2 = FillIndex2part(key1+key2);
+       Short_t normBin = SetNormBin(fillIndex2);
+       pVect1[0]=(fEvt+en1)->fTracks[i].fEaccepted; pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
+       pVect1[1]=(fEvt+en1)->fTracks[i].fP[0];      pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
+       pVect1[2]=(fEvt+en1)->fTracks[i].fP[1];      pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
+       pVect1[3]=(fEvt+en1)->fTracks[i].fP[2];      pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
+
+       qinv12 = GetQinv(fillIndex2, pVect1, pVect2);
+
+       if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting cut)
+       
+       ///////////////////////////////
+       ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
+       ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
+       
+       if(ch1 == ch2){
+         if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
+           fPairSplitCut[3][i]->AddAt('1',j);
+           continue;
+         }
+       }
+
+       if((qinv12 < fNormQcutLow[normBin]) || (qinv12 >= fNormQcutHigh[normBin])) continue;
+       
+       Int_t index1P23 = i;
+       Int_t index2P23 = j;
+       
+       for(Int_t pastpairP12=0; pastpairP12<numOtherPairs2[1][index1P23]; pastpairP12++){// loop in P12 with i as 2nd particle
+         Int_t locationOtherPairP12 = fOtherPairLocation2[1][index1P23]->At(pastpairP12);
+         if(locationOtherPairP12 < 0) continue; // no pair there
+         Int_t index1P12 = fNormPairs[1][ locationOtherPairP12 ].fIndex1;
+         
+                 
+         //Check other past pair status in P13
+         if( (fNormPairSwitch[2][index1P12]->At(index2P23)=='0')) continue;
+         
+         // all from different event
+         ch3 = Int_t((fNormPairs[1][ locationOtherPairP12 ].fCharge1 + 1)/2.);// charge of first particle in P12
+         key3 = fNormPairs[1][ locationOtherPairP12 ].fKey1;
+         Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
+         SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
+         
+         tempNormFillCount[bin1][bin2][bin3][fillIndex3][4]++;
+       }
+      }
+    }
+    
+    
+  
+    
+    ///////////////////////////////////////////////////  
+    // Do not use pairs from events with too many pairs
+    if(exitCode) {
+      cout<<"SE or ME or Norm PairCount too large.  Discarding all pairs and skipping event"<<endl;
+      (fEvt)->fNpairsSE = 0;
+      (fEvt)->fNpairsME = 0;
+      ((TH1F*)fOutputList->FindObject("fRejectedEvents"))->Fill(fMbin+1);
+      return;// Skip event
+    }else{
+      (fEvt)->fNpairsSE = pairCountSE;
+      (fEvt)->fNpairsME = pairCountME;  
+      ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
+    }
+    ///////////////////////////////////////////////////
+
+
+      
+    ///////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////
+    //
+    //
+    // Start the Main Correlation Analysis
+    //
+    //
+    ///////////////////////////////////////////////////////////////////////
+    
+    
+    
+  
+    /////////////////////////////////////////////////////////
+
+    // Set the Normalization counters
+    for(Int_t termN=0; termN<5; termN++){
+      
+      if(termN==0){
+       if((fEvt)->fNtracks ==0) continue;
+      }else if(termN<4){
+       if((fEvt)->fNtracks ==0) continue;
+       if((fEvt+1)->fNtracks ==0) continue;
+      }else {
+       if((fEvt)->fNtracks ==0) continue;
+       if((fEvt+1)->fNtracks ==0) continue;
+       if((fEvt+2)->fNtracks ==0) continue;
+      }
+     
+      for(Int_t sc=0; sc<kSCLimit3; sc++){
+       
+       for(Int_t c1=0; c1<2; c1++){
+         for(Int_t c2=0; c2<2; c2++){
+           for(Int_t c3=0; c3<2; c3++){
+             
+             if(sc==0 || sc==6 || sc==9){// Identical species
+               if( (c1+c2+c3)==1) {if(c1!=0 || c2!=0 || c3!=1) continue;}
+               if( (c1+c2+c3)==2) {if(c1!=0) continue;}
+             }else if(sc!=5){
+               if( (c1+c2)==1) {if(c1!=0) continue;}
+             }else {}// do nothing for pi-k-p case
+             Charge1[c1].Charge2[c2].Charge3[c3].SC[sc].MB[fMbin].EDB[fEDbin].ThreePT[termN].fNorm3->Fill(0.,tempNormFillCount[c1][c2][c3][sc][termN]);
+           }
+         }
+       }
+      }
+    }
+    
+    
+    
+    /////////////////////////////////////////////
+    // Calculate Pair-Cut Correlations
+    for(Int_t en1case=0; en1case<2; en1case++){// limit at 2 (normal)
+      
+      Int_t nump1=0;
+      if(en1case==0) nump1 = (fEvt)->fNpairsSE;
+      if(en1case==1) nump1 = (fEvt)->fNpairsME;
+     
+      // 1st pair
+      for(Int_t p1=0; p1<nump1; p1++){
+       
+       if(en1case==0){
+         ch1 = Int_t(((fEvt)->fPairsSE[p1].fCharge1 + 1)/2.);
+         ch2 = Int_t(((fEvt)->fPairsSE[p1].fCharge2 + 1)/2.);
+         pVect1[0] = (fEvt)->fPairsSE[p1].fE1; pVect2[0] = (fEvt)->fPairsSE[p1].fE2;
+         pVect1[1] = (fEvt)->fPairsSE[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsSE[p1].fP2[0]; 
+         pVect1[2] = (fEvt)->fPairsSE[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsSE[p1].fP2[1];
+         pVect1[3] = (fEvt)->fPairsSE[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsSE[p1].fP2[2];
+         index1 = (fEvt)->fPairsSE[p1].fIndex1; index2 = (fEvt)->fPairsSE[p1].fIndex2;
+         key1 = (fEvt)->fPairsSE[p1].fKey1; key2 = (fEvt)->fPairsSE[p1].fKey2;
+         qinv12 = (fEvt)->fPairsSE[p1].fQinv;
+         //
+         pVect1MC[1] = (fEvt)->fPairsSE[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsSE[p1].fP2MC[0];
+          pVect1MC[2] = (fEvt)->fPairsSE[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsSE[p1].fP2MC[1];
+          pVect1MC[3] = (fEvt)->fPairsSE[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsSE[p1].fP2MC[2];
+          pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
+          pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
+       }
+       if(en1case==1){
+         ch1 = Int_t(((fEvt)->fPairsME[p1].fCharge1 + 1)/2.);
+         ch2 = Int_t(((fEvt)->fPairsME[p1].fCharge2 + 1)/2.);
+         pVect1[0] = (fEvt)->fPairsME[p1].fE1; pVect2[0] = (fEvt)->fPairsME[p1].fE2; 
+         pVect1[1] = (fEvt)->fPairsME[p1].fP1[0]; pVect2[1] = (fEvt)->fPairsME[p1].fP2[0]; 
+         pVect1[2] = (fEvt)->fPairsME[p1].fP1[1]; pVect2[2] = (fEvt)->fPairsME[p1].fP2[1];
+         pVect1[3] = (fEvt)->fPairsME[p1].fP1[2]; pVect2[3] = (fEvt)->fPairsME[p1].fP2[2];
+         index1 = (fEvt)->fPairsME[p1].fIndex1; index2 = (fEvt)->fPairsME[p1].fIndex2;
+         key1 = (fEvt)->fPairsME[p1].fKey1; key2 = (fEvt)->fPairsME[p1].fKey2;
+         qinv12 = (fEvt)->fPairsME[p1].fQinv;
+         //
+         pVect1MC[1] = (fEvt)->fPairsME[p1].fP1MC[0]; pVect2MC[1] = (fEvt)->fPairsME[p1].fP2MC[0];
+          pVect1MC[2] = (fEvt)->fPairsME[p1].fP1MC[1]; pVect2MC[2] = (fEvt)->fPairsME[p1].fP2MC[1];
+          pVect1MC[3] = (fEvt)->fPairsME[p1].fP1MC[2]; pVect2MC[3] = (fEvt)->fPairsME[p1].fP2MC[2];
+          pVect1MC[0] = sqrt(pow(pVect1MC[1],2)+pow(pVect1MC[2],2)+pow(pVect1MC[3],2)+pow(fTrueMassPi,2));
+          pVect2MC[0] = sqrt(pow(pVect2MC[1],2)+pow(pVect2MC[2],2)+pow(pVect2MC[3],2)+pow(fTrueMassPi,2));
+       }
+       
+       
+       // en2 buffer
+       for(Int_t en2=0; en2<3; en2++){
+         //////////////////////////////////////
+
+         Bool_t skipcase=kTRUE;
+         Short_t config=-1, part=-1;
+         if(en1case==0 && en2==0) {skipcase=kFALSE; config=1; part=0;}// P11T1
+         if(en1case==0 && en2==1) {skipcase=kFALSE; config=2; part=1;}// P11T2
+         if(en1case==1 && en2==0) {skipcase=kFALSE; config=2; part=2;}// P12T1
+         if(en1case==1 && en2==2) {skipcase=kFALSE; config=3; part=3;}// P12T3
+                
+         if(skipcase) continue;
+       
+         
+         // 3-particle terms
+         // 3rd particle
+         for(Int_t k=0; k<(fEvt+en2)->fNtracks; k++){
+           index3 = k;
+           
+
+           // remove auto-correlations and duplicate triplets
+           if(config==1){
+             if( index1 == index3) continue;
+             if( index2 == index3) continue;
+             if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
+            
+             // skip the switched off triplets
+             if(fTripletSkip1[fPairLocationSE[index1]->At(index2)]->At(index3)=='1') {
+               fTripletSkip1[fPairLocationSE[index1]->At(index2)]->AddAt('0',index3);// Reset
+               continue;
+             }
+             ///////////////////////////////
+             // Turn off 1st possible degenerate triplet
+             if(index1 < index3){// verify correct id ordering ( index1 < k )
+               if(fPairLocationSE[index1]->At(index3) >= 0){
+                 fTripletSkip1[fPairLocationSE[index1]->At(index3)]->AddAt('1',index2);
+               }
+               if(fPairSplitCut[0][index1]->At(index3)=='1') continue;// Track splitting/merging
+             }else {// or k < index1
+               if(fPairLocationSE[index3]->At(index1) >= 0){
+                 fTripletSkip1[fPairLocationSE[index3]->At(index1)]->AddAt('1',index2);
+               }
+               if(fPairSplitCut[0][index3]->At(index1)=='1') continue;// Track splitting/merging
+             }
+             // turn off 2nd possible degenerate triplet
+             if(index2 < index3){// verify correct id ordering (index2 < k)
+               if(fPairLocationSE[index2]->At(index3) >= 0){
+                 fTripletSkip1[fPairLocationSE[index2]->At(index3)]->AddAt('1',index1);
+               }
+               if(fPairSplitCut[0][index2]->At(index3)=='1') continue;// Track splitting/merging
+             }else {// or k < index2
+               if(fPairLocationSE[index3]->At(index2) >= 0){
+                 fTripletSkip1[fPairLocationSE[index3]->At(index2)]->AddAt('1',index1);
+               }
+               if(fPairSplitCut[0][index3]->At(index2)=='1') continue;// Track splitting/merging
+             }
+
+           }// end config 1
+           
+           if(config==2 && part==1){// SE pair and third particle from next event. P11T2
+             ///////////////////////////////
+             // Turn off 1st possible degenerate triplet
+             if(fPairLocationME[index1]->At(index3) >= 0){
+               fTripletSkip2[fPairLocationME[index1]->At(index3)]->AddAt('1',index2);
+             }
+             
+             // turn off 2nd possible degenerate triplet
+             if(fPairLocationME[index2]->At(index3) >= 0){
+               fTripletSkip2[fPairLocationME[index2]->At(index3)]->AddAt('1',index1);
+             }
+             
+             if(fPairSplitCut[0][index1]->At(index2)=='1') continue;// Track splitting/merging
+             if(fPairSplitCut[1][index1]->At(index3)=='1') continue;// Track splitting/merging
+             if(fPairSplitCut[1][index2]->At(index3)=='1') continue;// Track splitting/merging
+           }// end config 2 part 1
+
+           if(config==2 && part==2){// P12T1
+             if( index1 == index3) continue;
+             
+             // skip the switched off triplets
+             if(fTripletSkip2[fPairLocationME[index1]->At(index2)]->At(index3)=='1') {
+               fTripletSkip2[fPairLocationME[index1]->At(index2)]->AddAt('0',index3);// Reset
+               continue;
+             }
+             // turn off another possible degenerate
+             if(fPairLocationME[index3]->At(index2) >= 0){
+               fTripletSkip2[fPairLocationME[index3]->At(index2)]->AddAt('1',index1);
+             }// end config 2 part 2
+
+             if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
+             if(index1 < index3) {if(fPairSplitCut[0][index1]->At(index3)=='1') continue;}// Track splitting/merging
+             else {if(fPairSplitCut[0][index3]->At(index1)=='1') continue;}// Track splitting/merging
+             if(fPairSplitCut[1][index3]->At(index2)=='1') continue;// Track splitting/merging
+           }
+           if(config==3){// P12T3
+             if(fPairSplitCut[1][index1]->At(index2)=='1') continue;// Track splitting/merging
+             if(fPairSplitCut[2][index1]->At(index3)=='1') continue;// Track splitting/merging
+             if(fPairSplitCut[3][index2]->At(index3)=='1') continue;// Track splitting/merging
+           }// end config 3
+           
+           
+
+           ch3 = Int_t(((fEvt+en2)->fTracks[k].fCharge + 1)/2.);
+           key3 = (fEvt+en2)->fTracks[k].fKey;
+           Short_t fillIndex3 = FillIndex3part(key1+key2+key3);
+           Short_t fillIndex13 = FillIndex2part(key1+key3);
+           Short_t fillIndex23 = FillIndex2part(key2+key3);
+           Short_t qCutBin13 = SetQcutBin(fillIndex13);
+           Short_t qCutBin23 = SetQcutBin(fillIndex23);
+           pVect3[0] = (fEvt+en2)->fTracks[k].fEaccepted;
+           pVect3[1] = (fEvt+en2)->fTracks[k].fP[0];
+           pVect3[2] = (fEvt+en2)->fTracks[k].fP[1];
+           pVect3[3] = (fEvt+en2)->fTracks[k].fP[2];
+           qinv13 = GetQinv(fillIndex13, pVect1, pVect3);
+           qinv23 = GetQinv(fillIndex23, pVect2, pVect3);
+
+           if(qinv13 < fQLowerCut) continue;
+           if(qinv23 < fQLowerCut) continue;
+           if(qinv13 > fQcut[qCutBin13]) continue;
+           if(qinv23 > fQcut[qCutBin23]) continue;
+
+          
+           
+           if(fMCcase){
+              pVect3MC[1] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPx;
+              pVect3MC[2] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPy;
+              pVect3MC[3] = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[k].fLabel)].fPz;
+              pVect3MC[0] = sqrt(pow(pVect3MC[1],2)+pow(pVect3MC[2],2)+pow(pVect3MC[3],2)+pow(fTrueMassPi,2));
+              qinv12MC = GetQinv(0, pVect1MC, pVect2MC);
+              qinv13MC = GetQinv(0, pVect1MC, pVect3MC);
+              qinv23MC = GetQinv(0, pVect2MC, pVect3MC);
+            }
+
+           
+           
+           // if all three pair cuts are the same then the case (config=2 && term=2) never reaches here.
+           
+           q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
+           transK3 = sqrt( pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
+           if(fEDbins>1){
+             if(transK3<0.3) fEDbin=0;
+             else fEDbin=1;
+           }
+           firstQ=0; secondQ=0; thirdQ=0;
+           
+           
+           //
+           
+           if(config==1) {// 123
+             SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 0, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
+             
+             if(fillIndex3 <= 2){
+               ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, 0, 1, firstQ, secondQ, thirdQ);
+               if(fillIndex3==0 && fMCcase) ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, 0, 1, firstQMC, secondQMC, thirdQMC);
+               Float_t WInput = 1.0;
+               if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, 1, 10, firstQ, secondQ, thirdQ);
+               //if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, 1, fFixedLambdaBinr3, firstQMC, secondQMC, thirdQMC);
+               ////
+               Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTermsQ3->Fill(q3, WInput);
+               Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[0].fTerms3->Fill(firstQ, secondQ, thirdQ, WInput);
+               ////
+               //
+               if(fillIndex3==0 && ch1==ch2 && ch1==ch3 && fMCcase==kFALSE){
+                 ((TH3D*)fOutputList->FindObject("fKt3DistTerm1"))->Fill(fMbin+1, transK3, q3);
+               }               
+               
+             }
+             
+           }else if(config==2){// 12, 13, 23
+             
+             Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
+             SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, part, bin1, bin2, bin3, fill2, fill3, fill4);
+         
+             // loop over terms 2-4
+             for(Int_t jj=2; jj<5; jj++){
+               if(jj==2) {if(!fill2) continue;}//12
+               if(jj==3) {if(!fill3) continue;}//13
+               if(jj==4) {if(!fill4) continue;}//23
+       
+               if(fillIndex3 <= 2){
+                 ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, jj, firstQ, secondQ, thirdQ);
+                 if(fillIndex3==0 && fMCcase) ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12MC, qinv13MC, qinv23MC, part, jj, firstQMC, secondQMC, thirdQMC);
+                 Float_t WInput = 1.0;
+                 if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, jj, 10, firstQ, secondQ, thirdQ);
+                 //if(fGenerateSignal && ch1==ch2 && ch1==ch3) WInput = MCWeight3D(kTRUE, jj, fFixedLambdaBinr3, firstQMC, secondQMC, thirdQMC);
+                 ////
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTermsQ3->Fill(q3, WInput);
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[jj-1].fTerms3->Fill(firstQ, secondQ, thirdQ, WInput);
+                 
+               }
+             }
+             
+           }else {// config 3: All particles from different events
+             
+             SetFillBins3(fillIndex3, key1, key2, key3, ch1, ch2, ch3, 3, bin1, bin2, bin3, fDummyB, fDummyB, fDummyB);
+             
+             if(ch1==ch2 && ch1==ch3 && fillIndex3==0) {
+               if(!fMCcase) ((TH3D*)fOutputList->FindObject("fKt3DistTerm5"))->Fill(fMbin+1, transK3, q3);
+             }       
+             
+             if(fillIndex3 <= 2){
+               ArrangeQs(fillIndex3, key1, key2, key3, ch1, ch2, ch3, qinv12, qinv13, qinv23, part, 5, firstQ, secondQ, thirdQ);
+               Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTermsQ3->Fill(q3);
+               Charge1[bin1].Charge2[bin2].Charge3[bin3].SC[fillIndex3].MB[fMbin].EDB[fEDbin].ThreePT[4].fTerms3->Fill(firstQ, secondQ, thirdQ);
+             }
+             
+             
+           }// config 3
+         }// end 3rd particle
+       }// en2
+       
+       
+      }// p1
+    }//en1
+    
+    ///////////////////
+  }// end of PdensityPairs
+
+  
+  // Post output data.
+  PostData(1, fOutputList);
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::Terminate(Option_t *) 
+{
+  // Called once at the end of the query
+  cout<<"Done"<<endl;
+
+}
+//________________________________________________________________________
+Bool_t AliThreePionRadii::AcceptPair(AliChaoticityTrackStruct first, AliChaoticityTrackStruct second)
+{
+  
+  if(fabs(first.fEta-second.fEta) > fMinSepPairEta) return kTRUE;
+  
+  // propagate through B field to r=1m
+  Float_t phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.15/first.fPt);// 0.15 for D=1m
+  if(phi1 > 2*PI) phi1 -= 2*PI;
+  if(phi1 < 0) phi1 += 2*PI;
+  Float_t phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.15/second.fPt);// 0.15 for D=1m 
+  if(phi2 > 2*PI) phi2 -= 2*PI;
+  if(phi2 < 0) phi2 += 2*PI;
+  
+  Float_t deltaphi = phi1 - phi2;
+  if(deltaphi > PI) deltaphi -= 2*PI;
+  if(deltaphi < -PI) deltaphi += 2*PI;
+  deltaphi = fabs(deltaphi);
+
+  if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
+    
+  
+  // propagate through B field to r=1.6m
+  phi1 = first.fPhi - asin(first.fCharge*(0.1*fBfield)*0.24/first.fPt);// mine. 0.24 for D=1.6m
+  if(phi1 > 2*PI) phi1 -= 2*PI;
+  if(phi1 < 0) phi1 += 2*PI;
+  phi2 = second.fPhi - asin(second.fCharge*(0.1*fBfield)*0.24/second.fPt);// mine. 0.24 for D=1.6m 
+  if(phi2 > 2*PI) phi2 -= 2*PI;
+  if(phi2 < 0) phi2 += 2*PI;
+  
+  deltaphi = phi1 - phi2;
+  if(deltaphi > PI) deltaphi -= 2*PI;
+  if(deltaphi < -PI) deltaphi += 2*PI;
+  deltaphi = fabs(deltaphi);
+
+  if(deltaphi < fMinSepPairPhi) return kFALSE;// Min Separation
+  
+  
+   
+  //
+  
+  Int_t ncl1 = first.fClusterMap.GetNbits();
+  Int_t ncl2 = second.fClusterMap.GetNbits();
+  Int_t sumCls = 0; Int_t sumSha = 0; Int_t sumQ = 0;
+  Double_t shfrac = 0; Double_t qfactor = 0;
+  for(Int_t imap = 0; imap < ncl1 && imap < ncl2; imap++) {
+    if (first.fClusterMap.TestBitNumber(imap) && second.fClusterMap.TestBitNumber(imap)) {// Both clusters
+      if (first.fSharedMap.TestBitNumber(imap) && second.fSharedMap.TestBitNumber(imap)) { // Shared
+       sumQ++;
+       sumCls+=2;
+       sumSha+=2;}
+      else {sumQ--; sumCls+=2;}
+    }
+    else if (first.fClusterMap.TestBitNumber(imap) || second.fClusterMap.TestBitNumber(imap)) {// Non shared
+      sumQ++;
+      sumCls++;}
+  }
+  if (sumCls>0) {
+    qfactor = sumQ*1.0/sumCls;
+    shfrac = sumSha*1.0/sumCls;
+  }
+  
+  if(qfactor > fShareQuality || shfrac > fShareFraction) return kFALSE;
+  
+  
+  return kTRUE;
+  
+
+}
+//________________________________________________________________________
+Float_t AliThreePionRadii::GamovFactor(Int_t chargeBin1, Int_t chargeBin2, Float_t qinv)
+{
+  Float_t arg = G_Coeff/qinv;
+  
+  if(chargeBin1==chargeBin2) return (exp(arg)-1)/(arg);
+  else {return (exp(-arg)-1)/(-arg);}
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::Shuffle(Int_t *iarr, Int_t i1, Int_t i2)
+{
+  Int_t j, k;
+  Int_t a = i2 - i1;
+  for (Int_t i = i1; i < i2+1; i++) {
+    j = (Int_t) (gRandom->Rndm() * a);
+    k = iarr[j];
+    iarr[j] = iarr[i];
+    iarr[i] = k;
+  }
+}
+//________________________________________________________________________
+Short_t AliThreePionRadii::FillIndex2part(Short_t key){
+
+  if(key==2) return 0;// pi-pi
+  else if(key==11) return 1;// pi-k
+  else if(key==101) return 2;// pi-p
+  else if(key==20) return 3;// k-k
+  else if(key==110) return 4;// k-p
+  else return 5;// p-p
+}
+//________________________________________________________________________
+Short_t AliThreePionRadii::FillIndex3part(Short_t key){
+  
+  if(key==3) return 0;// pi-pi-pi
+  else if(key==12) return 1;// pi-pi-k
+  else if(key==21) return 2;// k-k-pi
+  else if(key==102) return 3;// pi-pi-p
+  else if(key==201) return 4;// p-p-pi
+  else if(key==111) return 5;// pi-k-p
+  else if(key==30) return 6;// k-k-k
+  else if(key==120) return 7;// k-k-p
+  else if(key==210) return 8;// p-p-k
+  else return 9;// p-p-p
+  
+}
+//________________________________________________________________________
+Short_t AliThreePionRadii::SetQcutBin(Short_t fi){// fi=FillIndex
+  if(fi <= 2) return 0;
+  else if(fi==3) return 1;
+  else return 2;
+}
+//________________________________________________________________________
+Short_t AliThreePionRadii::SetNormBin(Short_t fi){// fi=FillIndex
+  if(fi==0) return 0;
+  else if(fi <= 2) return 1;
+  else return 2;
+}
+//________________________________________________________________________
+void AliThreePionRadii::SetFillBins2(Short_t fi, Short_t key1, Short_t key2, Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
+  
+  if(fi==0 || fi==3 || fi==5){// Identical species
+    if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
+    else {b1=c1; b2=c2;}
+  }else {// Mixed species
+    if(key1 < key2) { b1=c1; b2=c2;}
+    else {b1=c2; b2=c1;}
+  }
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::SetFillBins3(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Short_t part, Int_t &b1, Int_t &b2, Int_t &b3, Bool_t &fill2, Bool_t &fill3, Bool_t &fill4){
+  
+  
+  // seSS, seSK, SE_keysum only used to determine which terms to fill (only used for terms 2-4)
+  // part only matters for terms 2-4
+  Bool_t seSS=kFALSE;
+  Bool_t seSK=kFALSE;
+  Short_t seKeySum=0;// only used for pi-k-p case
+  if(part==1) {// default case (irrelevant for term 1 and term 5)
+    if(c1==c2) seSS=kTRUE;
+    if(key1==key2) seSK=kTRUE;
+    seKeySum = key1+key2;
+  }
+  if(part==2){
+    if(c1==c3) seSS=kTRUE;
+    if(key1==key3) seSK=kTRUE;
+    seKeySum = key1+key3;
+  }
+  
+  
+  // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
+  
+  if(fi==0 || fi==6 || fi==9){// Identical species
+    if( (c1+c2+c3)==1) {
+      b1=0; b2=0; b3=1;// Re-assign to merge degenerate histos
+      //
+      if(seSS) fill2=kTRUE;
+      else {fill3=kTRUE; fill4=kTRUE;}
+      //
+    }else if( (c1+c2+c3)==2) {
+      b1=0; b2=1; b3=1;
+      //
+      if(!seSS) {fill2=kTRUE; fill3=kTRUE;}
+      else fill4=kTRUE;
+      //
+    }else {
+      b1=c1; b2=c2; b3=c3;
+      fill2=kTRUE; fill3=kTRUE; fill4=kTRUE;
+    }
+  }else if(fi != 5){// all the rest except pi-k-p
+    if(key1==key2){
+      b3=c3;
+      if( (c1+c2)==1) {b1=0; b2=1;}
+      else {b1=c1; b2=c2;}
+    }else if(key1==key3){
+      b3=c2;
+      if( (c1+c3)==1) {b1=0; b2=1;}
+      else {b1=c1; b2=c3;}
+    }else {// Key2==Key3
+      b3=c1;
+      if( (c2+c3)==1) {b1=0; b2=1;}
+      else {b1=c2; b2=c3;}
+    }
+    //////////////////////////////
+    if(seSK) fill2=kTRUE;// Same keys from Same Event
+    else {// Different keys from Same Event
+      if( (c1+c2+c3)==1) {
+       if(b3==0) {
+         if(seSS) fill3=kTRUE;
+         else fill4=kTRUE;
+       }else{fill3=kTRUE; fill4=kTRUE;}// b3=1 so fill both
+      }else if( (c1+c2+c3)==2) {
+       if(b3==1) {
+         if(seSS) fill4=kTRUE;
+         else fill3=kTRUE;
+       }else{fill3=kTRUE; fill4=kTRUE;}// b3=0 so fill both
+      }else{fill3=kTRUE; fill4=kTRUE;}// all same charge so fill both
+    }
+    /////////////////////////////
+  }else {// pi-k-p  (no charge ordering applies since all are unique)
+    if(key1==1){
+      if(key2==10) {b1=c1; b2=c2; b3=c3;}// pi-k-p
+      else {b1=c1; b2=c3; b3=c2;}// pi-p-k
+    }else if(key1==10){
+      if(key2==1) {b1=c2; b2=c1; b3=c3;}// k-pi-p
+      else {b1=c3; b2=c1; b3=c2;}// k-p-pi
+    }else {// key1==100
+      if(key2==1) {b1=c2; b2=c3; b3=c1;}// p-pi-k
+      else {b1=c3; b2=c2; b3=c1;}// p-k-pi
+    }
+    ////////////////////////////////////
+    if(seKeySum==11) fill2=kTRUE;
+    else if(seKeySum==101) fill3=kTRUE;
+    else fill4=kTRUE;
+    ////////////////////////////////////
+  }
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::ArrangeQs(Short_t fi, Short_t key1, Short_t key2, Short_t key3, Int_t c1, Int_t c2, Int_t c3, Float_t q12, Float_t q13, Float_t q23, Short_t part, Short_t term, Float_t &fQ, Float_t &sQ, Float_t &tQ){
+  // for terms 2-4: start by setting q12(part 1) or q13(part 2)
+  if(fi==0 || fi==6 || fi==9){// Identical species
+    if( (c1+c2+c3)==1) {// fQ=ss, sQ=os, tQ=os
+      if(term==1 || term==5){
+       if(c1==c2) {fQ=q12; sQ=q13; tQ=q23;}
+       else if(c1==c3) {fQ=q13; sQ=q12; tQ=q23;}
+       else {fQ=q23; sQ=q12; tQ=q13;}
+      }else if(term==2 && part==1){
+       fQ=q12; sQ=q13; tQ=q23;
+      }else if(term==2 && part==2){
+       fQ=q13; sQ=q12; tQ=q23;
+      }else if(term==3 && part==1){
+       sQ=q12; 
+       if(c1==c3) {fQ=q13; tQ=q23;}
+       else {fQ=q23; tQ=q13;}
+      }else if(term==3 && part==2){
+       sQ=q13;
+       if(c1==c2) {fQ=q12; tQ=q23;}
+       else {fQ=q23; tQ=q12;}
+      }else if(term==4 && part==1){
+       tQ=q12;
+       if(c1==c3) {fQ=q13; sQ=q23;}
+       else {fQ=q23; sQ=q13;}
+      }else if(term==4 && part==2){
+       tQ=q13;
+       if(c1==c2) {fQ=q12; sQ=q23;}
+       else {fQ=q23; sQ=q12;}
+      }else cout<<"problem!!!!!!!!!!!!!"<<endl;
+    }else if( (c1+c2+c3)==2) {// fQ=os, sQ=os, tQ=ss
+      if(term==1 || term==5){
+       if(c1==c2) {tQ=q12; sQ=q13; fQ=q23;}
+       else if(c1==c3) {tQ=q13; sQ=q12; fQ=q23;}
+       else {tQ=q23; sQ=q12; fQ=q13;}
+      }else if(term==2 && part==1){
+       fQ=q12; 
+       if(c1==c3) {tQ=q13; sQ=q23;}
+       else {tQ=q23; sQ=q13;}
+      }else if(term==2 && part==2){
+       fQ=q13; 
+       if(c1==c2) {tQ=q12; sQ=q23;}
+       else {tQ=q23; sQ=q12;}
+      }else if(term==3 && part==1){
+       sQ=q12; 
+       if(c1==c3) {tQ=q13; fQ=q23;}
+       else {tQ=q23; fQ=q13;}
+      }else if(term==3 && part==2){
+       sQ=q13; 
+       if(c1==c2) {tQ=q12; fQ=q23;}
+       else {tQ=q23; fQ=q12;}
+      }else if(term==4 && part==1){
+       tQ=q12; sQ=q13; fQ=q23;
+      }else if(term==4 && part==2){
+       tQ=q13; sQ=q12; fQ=q23;
+      }else cout<<"problem!!!!!!!!!!!!!"<<endl;
+    }else {// fQ=ss, sQ=ss, tQ=ss
+      if(term==1 || term==5) {fQ=q12; sQ=q13; tQ=q23;}
+      else if(term==2 && part==1) {fQ=q12; sQ=q13; tQ=q23;}
+      else if(term==2 && part==2) {fQ=q13; sQ=q12; tQ=q23;}
+      else if(term==3 && part==1) {sQ=q12; fQ=q13; tQ=q23;}
+      else if(term==3 && part==2) {sQ=q13; fQ=q12; tQ=q23;}
+      else if(term==4 && part==1) {tQ=q12; fQ=q13; sQ=q23;}
+      else if(term==4 && part==2) {tQ=q13; fQ=q12; sQ=q23;}
+    }
+  }else if(fi != 5){// all the rest except pi-k-p      
+    if(key1==key2){
+      fQ=q12;
+      if(c1==c2){
+       // cases not explicity shown below are not possible
+       if(term==1 || term==5) {sQ=q13; tQ=q23;}
+       else if(term==2 && part==1) {sQ=q13; tQ=q23;}
+       else if(term==3 && part==2) {sQ=q13; tQ=q23;}
+       else if(term==4 && part==2) {tQ=q13; sQ=q23;}
+       else cout<<"problem!!!!!!!!!!!!!"<<endl;
+      }else if(c3==0){
+       if(c1==c3) {sQ=q13; tQ=q23;}
+       else {sQ=q23; tQ=q13;}
+      }else {//c3==1
+       if(c1==c3) {tQ=q13; sQ=q23;}
+       else {tQ=q23; sQ=q13;}
+      }
+    }else if(key1==key3){
+      fQ=q13;
+      if(c1==c3){
+       // cases not explicity shown below are not possible
+       if(term==1 || term==5) {sQ=q12; tQ=q23;}
+       else if(term==2 && part==2) {sQ=q12; tQ=q23;}
+       else if(term==3 && part==1) {sQ=q12; tQ=q23;}
+       else if(term==4 && part==1) {tQ=q12; sQ=q23;}
+       else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
+      }else if(c2==0){
+       if(c1==c2) {sQ=q12; tQ=q23;}
+       else {sQ=q23; tQ=q12;}
+      }else {//c2==1
+       if(c1==c2) {tQ=q12; sQ=q23;}
+       else {tQ=q23; sQ=q12;}
+      }
+    }else {// key2==key3
+      fQ=q23;
+      if(c2==c3){
+       // cases not explicity shown below are not possible
+       if(term==1 || term==5) {sQ=q12; tQ=q13;}
+       else if(term==3 && part==1) {sQ=q12; tQ=q13;}
+       else if(term==3 && part==2) {sQ=q13; tQ=q12;}
+       else if(term==4 && part==1) {tQ=q12; sQ=q13;}
+       else if(term==4 && part==2) {tQ=q13; sQ=q12;}
+       else cout<<"problem!!!!!!!!!!!!!!!!!!!!!!"<<endl;
+      }else if(c1==0){
+       if(c1==c2) {sQ=q12; tQ=q13;}
+       else {sQ=q13; tQ=q12;}
+      }else {//c1==1
+       if(c1==c2) {tQ=q12; sQ=q13;}
+       else {tQ=q13; sQ=q12;}
+      }
+    }
+  }else {// pi-k-p
+    if(key1==1){
+      if(key2==10) {fQ=q12; sQ=q13; tQ=q23;}// pi-k-p
+      else {fQ=q13; sQ=q12; tQ=q23;}// pi-p-k
+    }else if(key1==10){
+      if(key2==1) {fQ=q12; sQ=q23; tQ=q13;}// k-pi-p
+      else {fQ=q13; sQ=q23; tQ=q12;}// k-p-pi
+    }else {// key1==100
+      if(key2==1) {fQ=q23; sQ=q12; tQ=q13;}// p-pi-k
+      else {fQ=q23; sQ=q13; tQ=q12;}// p-k-pi
+    }
+    
+  }
+
+
+}
+//________________________________________________________________________
+Float_t AliThreePionRadii::GetQinv(Short_t fi, Float_t track1[], Float_t track2[]){
+  
+  Float_t qinv=1.0;
+  
+  if(fi==0 || fi==3 || fi==5){// identical masses
+    qinv = sqrt( pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2));
+  }else{// different masses
+    Float_t px = track1[1] + track2[1]; 
+    Float_t py = track1[2] + track2[2]; 
+    Float_t pz = track1[3] + track2[3];    
+    Float_t pSquared = pow(track1[0]+track2[0],2) - px*px - py*py - pz*pz;
+    Float_t deltaDOTsum = (track1[0]-track2[0])*(track1[0]+track2[0]);
+    deltaDOTsum -= (track1[1]-track2[1])*px + (track1[2]-track2[2])*py + (track1[3]-track2[3])*pz;
+    
+    qinv =  pow( (track1[1]-track2[1]) - deltaDOTsum*px/(pSquared),2);
+    qinv += pow( (track1[2]-track2[2]) - deltaDOTsum*py/(pSquared),2);
+    qinv += pow( (track1[3]-track2[3]) - deltaDOTsum*pz/(pSquared),2);
+    qinv -= pow( (track1[0]-track2[0]) - deltaDOTsum*(track1[0]+track2[0])/(pSquared),2);
+    qinv = sqrt(qinv);
+  }
+  
+  return qinv;
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::GetQosl(Float_t track1[], Float_t track2[], Float_t& qout, Float_t& qside, Float_t& qlong){
+  Float_t p0 = track1[0] + track2[0];
+  Float_t px = track1[1] + track2[1];
+  Float_t py = track1[2] + track2[2];
+  Float_t pz = track1[3] + track2[3];
+  
+  Float_t mt = sqrt(p0*p0 - pz*pz);
+  Float_t pt = sqrt(px*px + py*py);
+  
+  Float_t v0 = track1[0] - track2[0];
+  Float_t vx = track1[1] - track2[1];
+  Float_t vy = track1[2] - track2[2];
+  Float_t vz = track1[3] - track2[3];
+  
+  qout = (px*vx + py*vy)/pt;
+  qside = (px*vy - py*vx)/pt;
+  qlong = (p0*vz - pz*v0)/mt;
+}
+//________________________________________________________________________
+Float_t AliThreePionRadii::MCWeight(Int_t charge1, Int_t charge2, Float_t r, Int_t dampIndex, Float_t qinv){
+  
+  Float_t radius = r/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
+  
+  Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
+  Float_t coulCorr12 = FSICorrelation2(charge1, charge2, qinv);
+  if(charge1==charge2){
+    Float_t arg=qinv*radius;
+    Float_t EW = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg,3) - 12.*arg);
+    EW += kappa4/(24.*pow(2.,2))*(16.*pow(arg,4) -48.*pow(arg,2) + 12);
+    return ((1-myDamp) + myDamp*(1 + exp(-pow(qinv*radius,2))*pow(EW,2))*coulCorr12);
+  }else {
+    return ((1-myDamp) + myDamp*coulCorr12);
+  }
+    
+}
+//________________________________________________________________________
+Float_t AliThreePionRadii::MCWeight3D(Bool_t SameCharge, Int_t term, Int_t dampIndex, Float_t q12, Float_t q13, Float_t q23){
+  if(term==5) return 1.0;
+  
+  Float_t radius=fRMax;
+  radius /= 0.19733;
+
+  Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
+  Float_t fc = sqrt(myDamp);
+  
+  if(SameCharge){// all three of the same charge
+    Float_t coulCorr12 = FSICorrelation2(+1,+1, q12);// K2
+    Float_t coulCorr13 = FSICorrelation2(+1,+1, q13);// K2
+    Float_t coulCorr23 = FSICorrelation2(+1,+1, q23);// K2
+    Float_t arg12=q12*radius;
+    Float_t arg13=q13*radius;
+    Float_t arg23=q23*radius;
+    Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
+    EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
+    Float_t EW13 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13);
+    EW13 += kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12);
+    Float_t EW23 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23);
+    EW23 += kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12);
+    if(term==1){
+      Float_t c3QS = 1 + exp(-pow(q12*radius,2))*pow(EW12,2) + exp(-pow(q13*radius,2))*pow(EW13,2) + exp(-pow(q23*radius,2))*pow(EW23,2);
+      c3QS += 2*exp(-pow(radius,2)*(pow(q12,2) + pow(q13,2) + pow(q23,2))/2.)*EW12*EW13*EW23;
+      Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2))*pow(EW12,2))*coulCorr12;
+      w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q13*radius,2))*pow(EW13,2))*coulCorr13;
+      w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q23*radius,2))*pow(EW23,2))*coulCorr23;
+      w123 += pow(fc,3)*c3QS*coulCorr12*coulCorr13*coulCorr23;
+      return w123;
+    }else if(term==2){
+      return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2))*pow(EW12,2))*coulCorr12);
+    }else if(term==3){
+      return ((1-myDamp) + myDamp*(1 + exp(-pow(q13*radius,2))*pow(EW13,2))*coulCorr13);
+    }else if(term==4){
+      return ((1-myDamp) + myDamp*(1 + exp(-pow(q23*radius,2))*pow(EW23,2))*coulCorr23);
+    }else return 1.0;
+  
+  }else{// mixed charge case pair 12 always treated as ss
+    Float_t coulCorr12 = FSICorrelation2(+1,+1, q12);// K2 ss
+    Float_t coulCorr13 = FSICorrelation2(+1,-1, q13);// K2 os
+    Float_t coulCorr23 = FSICorrelation2(+1,-1, q23);// K2 os
+    Float_t arg12=q12*radius;
+    Float_t EW12 = 1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12);
+    EW12 += kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12);
+    if(term==1){
+      Float_t c3QS = 1 + exp(-pow(q12*radius,2))*pow(EW12,2);
+      Float_t w123 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      w123 += pow(fc,2)*(1-fc)*(1+exp(-pow(q12*radius,2))*pow(EW12,2))*coulCorr12;
+      w123 += pow(fc,2)*(1-fc)*coulCorr13;
+      w123 += pow(fc,2)*(1-fc)*coulCorr23;
+      w123 += pow(fc,3)*c3QS*coulCorr12*coulCorr13*coulCorr23;
+      return w123;
+    }else if(term==2){
+      return ((1-myDamp) + myDamp*(1 + exp(-pow(q12*radius,2))*pow(EW12,2))*coulCorr12);
+    }else if(term==3){
+      return ((1-myDamp) + myDamp*coulCorr13);
+    }else if(term==4){
+      return ((1-myDamp) + myDamp*coulCorr23);
+    }else return 1.0;
+  }
+  
+}
+//________________________________________________________________________
+void AliThreePionRadii::SetFSICorrelations(Bool_t legoCase, TH1D *temp1DSS[10], TH1D *temp1DOS[10]){
+  // read in 2-particle FSI correlations = K2 
+  // 2-particle input histo from file is binned in qinv.
+  if(legoCase){
+    cout<<"LEGO call to SetFSICorrelations"<<endl;
+    for(int mb=0; mb<10; mb++){
+      fFSI2SS[mb] = (TH1D*)temp1DSS[mb]->Clone();
+      fFSI2OS[mb] = (TH1D*)temp1DOS[mb]->Clone();
+      //
+      fFSI2SS[mb]->SetDirectory(0);
+      fFSI2OS[mb]->SetDirectory(0);
+    }
+  }else {
+    cout<<"non LEGO call to SetFSICorrelations"<<endl;
+    TFile *fsifile = new TFile("KFile.root","READ");
+    if(!fsifile->IsOpen()) {
+      cout<<"No FSI file found"<<endl;
+      AliFatal("No FSI file found.  Kill process.");
+    }else {cout<<"Good FSI File Found!"<<endl;}
+    for(int mb=0; mb<10; mb++){
+      TString *nameK2ss = new TString("K2ss_");
+      TString *nameK2os = new TString("K2os_");
+      *nameK2ss += mb;
+      *nameK2os += mb;
+      TH1D *temphistoSS = (TH1D*)fsifile->Get(nameK2ss->Data());
+      TH1D *temphistoOS = (TH1D*)fsifile->Get(nameK2os->Data());
+      //
+      fFSI2SS[mb] = (TH1D*)temphistoSS->Clone();
+      fFSI2OS[mb] = (TH1D*)temphistoOS->Clone();
+      fFSI2SS[mb]->SetDirectory(0);
+      fFSI2OS[mb]->SetDirectory(0);
+    }
+    
+    fsifile->Close();
+  }
+
+  for(int mb=0; mb<10; mb++){
+    for(Int_t ii=1; ii<=fFSI2SS[mb]->GetNbinsX(); ii++){
+      if(fFSI2SS[mb]->GetBinContent(ii) > 1.0) fFSI2SS[mb]->SetBinContent(ii, 1.0);
+      if(fFSI2OS[mb]->GetBinContent(ii) > 10.0) fFSI2OS[mb]->SetBinContent(ii, 10.0);
+      //
+      if(fFSI2SS[mb]->GetBinContent(ii) < 0.05) fFSI2SS[mb]->SetBinContent(ii, 0.05);
+      if(fFSI2OS[mb]->GetBinContent(ii) < 0.9) fFSI2OS[mb]->SetBinContent(ii, 0.9);
+    }
+  }
+  
+  cout<<"Done reading FSI file"<<endl;
+}
+//________________________________________________________________________
+Float_t AliThreePionRadii::FSICorrelation2(Int_t charge1, Int_t charge2, Float_t qinv){
+  // returns 2-particle Coulomb correlations = K2
+  Int_t qbinL = fFSI2SS[fFSIindex]->GetYaxis()->FindBin(qinv-fFSI2SS[fFSIindex]->GetYaxis()->GetBinWidth(1)/2.);
+  Int_t qbinH = qbinL+1;
+  if(qbinL <= 0) return 1.0;
+  if(qbinH > fFSI2SS[fFSIindex]->GetNbinsY()) return 1.0;
+  
+  Float_t slope=0;
+  if(charge1==charge2){
+    slope = fFSI2SS[fFSIindex]->GetBinContent(qbinL) - fFSI2SS[fFSIindex]->GetBinContent(qbinH);
+    slope /= fFSI2SS[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSI2SS[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
+    return (slope*(qinv - fFSI2SS[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSI2SS[fFSIindex]->GetBinContent(qbinL));
+  }else {
+    slope = fFSI2OS[fFSIindex]->GetBinContent(qbinL) - fFSI2OS[fFSIindex]->GetBinContent(qbinH);
+    slope /= fFSI2OS[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSI2OS[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
+    return (slope*(qinv - fFSI2OS[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSI2OS[fFSIindex]->GetBinContent(qbinL));
+  }
+}
+//________________________________________________________________________
diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliThreePionRadii.h b/PWGCF/FEMTOSCOPY/Chaoticity/AliThreePionRadii.h
new file mode 100644 (file)
index 0000000..1162bea
--- /dev/null
@@ -0,0 +1,253 @@
+#ifndef AliThreePionRadii_cxx
+#define AliThreePionRadii_cxx
+//
+// Class AliThreePionRadii
+//
+// AliThreePionRadii
+// author:
+//        Dhevan Gangadharan (dhevan.raja.gangadharan@cern.ch)
+//
+
+
+
+class TH1F;
+class TH3F;
+class TH1D;
+class TH2D;
+class TH3D;
+
+class TProfile;
+class TProfile2D;
+class TRandom3;
+
+class AliESDEvent;
+class AliAODEvent;
+class AliESDtrackCuts;
+class AliESDpid;
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliESDpid.h"
+#include "AliAODPid.h"
+#include "AliChaoticityEventCollection.h"
+#include "AliThreePionRadii.h"
+
+class AliThreePionRadii : public AliAnalysisTaskSE {
+ public:
+
+  AliThreePionRadii();
+  AliThreePionRadii(const Char_t *name);
+  virtual ~AliThreePionRadii();
+  AliThreePionRadii(const AliThreePionRadii &obj); 
+  AliThreePionRadii &operator=(const AliThreePionRadii &obj);
+  
+
+  virtual void   UserCreateOutputObjects();
+  virtual void   Exec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+
+  enum {
+    kPairLimit = 15000,//15000
+    kNormPairLimit = 45000,
+    kMultLimitPbPb = 2000,//2000
+    kMultLimitPP = 300,
+    kMCarrayLimit = 110000,
+    kQbins = 20,
+    kQbinsWeights = 40,
+    kQbinsPP = 40,
+    kQbinsWeightsPP = 40,
+    kNDampValues = 16,
+    kRmin = 2,// min radii for Momentum resolution calculations
+    kDENtypes = 1,// was (kRVALUES)*kNDampValues
+    kSCLimit2 = 1,// 1, 6
+    kSCLimit3 = 1// 1, 10
+  };
+
+  static const Int_t fEDbins   = 1;
+  static const Int_t fCentBins = 20;// 0-50%
+  static const Int_t fRVALUES  = 10;// 
+
+
+  Int_t GetNumRValues() const {return AliThreePionRadii::fRVALUES;}
+  Int_t GetNumCentBins() const {return AliThreePionRadii::fCentBins;}
+  Int_t GetNumEDBins() const {return AliThreePionRadii::fEDbins;}
+  void SetFSICorrelations(Bool_t legoCase=kTRUE, TH1D *temp1DSS[10]=0x0, TH1D *temp1DOS[10]=0x0);
+  //
+  void SetMCdecision(Bool_t mc) {fMCcase = mc;}
+  void SetPbPbCase(Bool_t pbpb) {fPbPbcase = pbpb;}
+  void SetGenerateSignal(Bool_t gen) {fGenerateSignal = gen;}
+  void SetCentBinRange(Int_t low, Int_t high) {fCentBinLowLimit = low; fCentBinHighLimit = high;}
+  void SetLEGOCase(Bool_t lego) {fLEGO = lego;}
+  void SetFilterBit(UInt_t filterbit) {fFilterBit = filterbit;}
+  void SetMaxChi2NDF(Float_t MaxChi2NDF) {fMaxChi2NDF = MaxChi2NDF;}
+  void SetMinTPCncls(Int_t MinTPCncls) {fMinTPCncls = MinTPCncls;}
+  void SetPairSeparationCutEta(Float_t pairsep) {fMinSepPairEta = pairsep;}
+  void SetPairSeparationCutPhi(Float_t pairsep) {fMinSepPairPhi = pairsep;}
+  void SetNsigmaTPC(Float_t nsig) {fSigmaCutTPC = nsig;}
+  void SetNsigmaTOF(Float_t nsig) {fSigmaCutTOF = nsig;}
+  void SetRMax(Int_t rbin) {fRMax = rbin;}
+  //
+
+
+ private:
+
+  void ParInit();
+  Bool_t AcceptPair(AliChaoticityTrackStruct, AliChaoticityTrackStruct);
+  Float_t GamovFactor(Int_t, Int_t, Float_t);
+  void Shuffle(Int_t*, Int_t, Int_t);
+  Short_t FillIndex2part(Short_t);
+  Short_t FillIndex3part(Short_t);
+  Short_t SetQcutBin(Short_t);
+  Short_t SetNormBin(Short_t);
+  void SetFillBins2(Short_t, Short_t, Short_t, Int_t, Int_t, Int_t&, Int_t&);
+  void SetFillBins3(Short_t, Short_t, Short_t, Short_t, Int_t, Int_t, Int_t, Short_t, Int_t&, Int_t&, Int_t&, Bool_t&, Bool_t&, Bool_t&);
+  void ArrangeQs(Short_t, Short_t, Short_t, Short_t, Int_t, Int_t, Int_t, Float_t, Float_t, Float_t, Short_t, Short_t, Float_t&, Float_t&, Float_t&);
+  Float_t GetQinv(Short_t, Float_t[], Float_t[]);
+  void GetQosl(Float_t[], Float_t[], Float_t&, Float_t&, Float_t&);
+  Float_t FSICorrelation2(Int_t, Int_t, Float_t);
+  Float_t MCWeight(Int_t, Int_t, Float_t, Int_t, Float_t);
+  Float_t MCWeight3D(Bool_t, Int_t, Int_t, Float_t, Float_t, Float_t);
+  //
+  
+  
+  const char* fname;// name of class
+  AliAODEvent            *fAOD; //!    // AOD object
+  TList                  *fOutputList; //! Compact Output list
+  AliPIDResponse         *fPIDResponse; //! PID response object; equivalent to AliAODpidUtil
+  
+  
+  AliChaoticityEventCollection ***fEC; //!
+  AliChaoticityEventStruct *fEvt; //!
+  AliChaoticityTrackStruct *fTempStruct; //!
+  TRandom3* fRandomNumber; //!
+
+  
+  //////////////////////////////
+  // histogram structures
+
+  struct St6 {
+    TH1D *fExplicit3; //!
+    TH1D *fNormEx3; //!
+    //
+    TH1D *fNorm3; //!
+    TH3D *fTerms3; //!
+    TH1D *fTermsQ3; //!
+    TH1D *fIdeal; //!
+    TH1D *fSmeared; //!
+  };
+  struct St5 {
+    TH2D *fExplicit2; //!
+    TH2D *fExplicit2QW; //!
+    TH3D *fExplicit2ThreeD; //!
+    TProfile2D *fAvgP; //!
+    TH2D *fIdeal; //!
+    TH2D *fSmeared; //!
+    //
+    TH1D *fMCqinv; //!
+    TH1D *fMCqinvQW; //!
+    TH2D *fPIDpurityDen; //!
+    TH2D *fPIDpurityNum; //!
+  };
+  struct St_EDB {// SC structure
+    struct St5 TwoPT[2];
+    struct St6 ThreePT[5];
+  };
+  struct St_M {
+    struct St_EDB EDB[fEDbins];
+  };
+  struct St4 {
+    struct St_M MB[fCentBins];
+  };
+  struct St3 {
+    struct St4 SC[kSCLimit3];
+  };
+  struct St2 {
+    struct St3 Charge3[2];
+    struct St4 SC[kSCLimit2];
+  };
+  struct St1 {
+    struct St2 Charge2[2];
+  };
+  struct St1 Charge1[2];//!
+
+
+  Bool_t fLEGO;
+  Bool_t fMCcase;
+  Bool_t fAODcase;
+  Bool_t fPbPbcase;
+  Bool_t fGenerateSignal;
+  Bool_t fPdensityPairCut;
+  Int_t fRMax;
+  UInt_t fFilterBit;
+  Float_t fMaxChi2NDF;
+  Int_t fMinTPCncls;
+  Double_t fBfield;
+  Int_t fMbin;
+  Int_t fFSIindex;
+  Int_t fEDbin;
+  Int_t fMbins;
+  Int_t fMultLimit;      
+  Int_t fCentBinLowLimit;
+  Int_t fCentBinHighLimit;
+  Int_t fEventCounter;
+  Int_t fEventsToMix;
+  Int_t fZvertexBins;
+  Int_t fMultLimits[fCentBins+1];
+  Float_t fQcut[3];
+  Float_t fQLowerCut;
+  Float_t fQlimitC2;
+  Int_t fQbinsC2;
+  Float_t fNormQcutLow[3];
+  Float_t fNormQcutHigh[3];
+  Float_t fKupperBound;
+  Float_t fQupperBound;
+  Float_t fQupperBoundWeights;
+  Int_t   fQbins;
+  Int_t   fQbinsWeights;
+  Float_t fQstep;
+  Float_t fQstepWeights;
+  Float_t fQmean[kQbinsWeights];
+  Float_t fQmeanPP[kQbinsWeightsPP];
+  Float_t fDampStart;
+  Float_t fDampStep;
+  
+  Float_t fTPCTOFboundry;
+  Float_t fTOFboundry;
+  Float_t fSigmaCutTPC;
+  Float_t fSigmaCutTOF;
+  
+  Float_t fMinSepPairEta;
+  Float_t fMinSepPairPhi;
+  Float_t fShareQuality;
+  Float_t fShareFraction;
+  
+  Float_t fTrueMassP, fTrueMassPi, fTrueMassK, fTrueMassKs, fTrueMassLam;
+
+  Bool_t fDummyB;
+
+  
+  Char_t fDefaultsCharMult[kMultLimitPbPb];//!
+  Char_t fDefaultsCharSE[kPairLimit];//!
+  Char_t fDefaultsCharME[2*kPairLimit];//!
+  Int_t fDefaultsInt[kMultLimitPbPb];//!
+  TArrayI *fPairLocationSE[kMultLimitPbPb];//!
+  TArrayI *fPairLocationME[kMultLimitPbPb];//!
+  TArrayC *fTripletSkip1[kPairLimit];//!
+  TArrayC *fTripletSkip2[2*kPairLimit];//!
+  TArrayI *fOtherPairLocation1[2][kMultLimitPbPb];//!
+  TArrayI *fOtherPairLocation2[2][kMultLimitPbPb];//!
+  TArrayC *fNormPairSwitch[3][kMultLimitPbPb];//!
+  TArrayC *fPairSplitCut[4][kMultLimitPbPb];//!
+  
+  AliChaoticityNormPairStruct *fNormPairs[3];//!
+  
+ public:
+  TH1D *fFSI2SS[10];
+  TH1D *fFSI2OS[10];
+
+  ClassDef(AliThreePionRadii, 1); 
+};
+
+#endif