Add 4-pion femto analysis code
authordgangadh <dhevan.raja.gangadharan@cern.ch>
Thu, 8 May 2014 15:15:47 +0000 (17:15 +0200)
committerdgangadh <dhevan.raja.gangadharan@cern.ch>
Thu, 8 May 2014 15:15:47 +0000 (17:15 +0200)
PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx [new file with mode: 0755]
PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h [new file with mode: 0755]
PWGCF/FEMTOSCOPY/Chaoticity/AliFourPionEventCollection.cxx [new file with mode: 0755]
PWGCF/FEMTOSCOPY/Chaoticity/AliFourPionEventCollection.h [new file with mode: 0755]

diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.cxx
new file mode 100755 (executable)
index 0000000..a2af30a
--- /dev/null
@@ -0,0 +1,3568 @@
+#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 "AliAnalysisUtils.h"
+
+#include "AliFourPion.h"
+
+#define PI 3.1415927
+#define G_Coeff 0.006399 // 2*pi*alpha*M_pion
+#define kappa3 0.1 // kappa3 Edgeworth coefficient (non-Gaussian features of C2)
+#define kappa4 0.5 // kappa4 Edgeworth coefficient (non-Gaussian features of C2)
+
+
+// Author: Dhevan Gangadharan
+
+ClassImp(AliFourPion)
+
+//________________________________________________________________________
+AliFourPion::AliFourPion():
+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),
+  fGeneratorOnly(kFALSE),
+  fPdensityPairCut(kTRUE),
+  fTabulatePairs(kFALSE),
+  fRMax(11),
+  ffcSq(0.7),
+  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(0),
+  fQLowerCut(0),
+  fNormQcutLow(0),
+  fNormQcutHigh(0),
+  fKupperBound(0),
+  fQupperBoundQ2(0),
+  fQupperBoundQ3(0),
+  fQupperBoundQ4(0),
+  fQbinsQ2(1),
+  fQbinsQ3(1),
+  fQbinsQ4(1),
+  fQupperBoundWeights(0),
+  fKstepT(),
+  fKstepY(),
+  fKmeanT(),
+  fKmeanY(),
+  fKmiddleT(),
+  fKmiddleY(),
+  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),
+  fKtIndexL(0),
+  fKtIndexH(0),
+  fQoIndexL(0),
+  fQoIndexH(0),
+  fQsIndexL(0),
+  fQsIndexH(0),
+  fQlIndexL(0),
+  fQlIndexH(0),
+  fDummyB(0),
+  fKT3transition(0.3),
+  fKT4transition(0.3),
+  fMomResC2(0x0)
+{
+  // 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 term=0; term<2; term++){
+           
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
+           
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
+           
+         }// term_2
+         
+         
+         for(Int_t c3=0; c3<2; c3++){
+           for(Int_t term=0; term<5; term++){
+             
+             Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
+             Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
+             Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
+                     
+           }// term_3
+         }//c3
+       }//c2
+      }//c1
+      for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
+       for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
+         KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
+         KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
+       }
+      }
+      
+    }// ED
+  }// Mbin
+  
+  // Initialize FSI histograms
+  for(Int_t i=0; i<12; i++){
+    fFSIss[i]=0x0; 
+    fFSIos[i]=0x0;
+  }
+
+
+  // Initialize fNormWeight and fNormWeightErr to 0
+  for(Int_t i=0; i<3; i++){// Kt iterator
+    for(Int_t j=0; j<10; j++){// Mbin iterator
+      fNormWeight[i][j]=0x0;
+    }
+  }
+  
+
+}
+//________________________________________________________________________
+AliFourPion::AliFourPion(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),
+  fGeneratorOnly(kFALSE),
+  fPdensityPairCut(kTRUE),
+  fTabulatePairs(kFALSE),
+  fRMax(11),
+  ffcSq(0.7),
+  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(0),
+  fQLowerCut(0),
+  fNormQcutLow(0),
+  fNormQcutHigh(0),
+  fKupperBound(0),
+  fQupperBoundQ2(0),
+  fQupperBoundQ3(0),
+  fQupperBoundQ4(0),
+  fQbinsQ2(1),
+  fQbinsQ3(1),
+  fQbinsQ4(1),
+  fQupperBoundWeights(0),
+  fKstepT(),
+  fKstepY(),
+  fKmeanT(),
+  fKmeanY(),
+  fKmiddleT(),
+  fKmiddleY(),
+  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),
+  fKtIndexL(0),
+  fKtIndexH(0),
+  fQoIndexL(0),
+  fQoIndexH(0),
+  fQsIndexL(0),
+  fQsIndexH(0),
+  fQlIndexL(0),
+  fQlIndexH(0),
+  fDummyB(0),
+  fKT3transition(0.3),
+  fKT4transition(0.3),
+  fMomResC2(0x0)
+{
+  // 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 term=0; term<2; term++){
+           
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2=0x0;
+           
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = 0x0;
+           Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = 0x0;
+           
+         }// term_2
+         
+         for(Int_t c3=0; c3<2; c3++){
+           for(Int_t term=0; term<5; term++){
+             
+             Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3 = 0x0;
+             Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = 0x0;
+             Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = 0x0;
+             
+           }// term_3
+         }//c3
+       }//c2
+      }//c1
+      
+      for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
+       for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
+         KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = 0x0;
+         KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = 0x0;
+       }
+      }
+      
+    }// ED
+  }// Mbin
+  
+  // Initialize FSI histograms
+  for(Int_t i=0; i<12; i++){
+    fFSIss[i]=0x0; 
+    fFSIos[i]=0x0;
+  }
+  
+  // Initialize fNormWeight and fNormWeightErr to 0
+  for(Int_t i=0; i<3; i++){// Kt iterator
+    for(Int_t j=0; j<10; j++){// Mbin iterator
+      fNormWeight[i][j]=0x0;
+    }
+  }
+
+
+  DefineOutput(1, TList::Class());
+}
+//________________________________________________________________________
+AliFourPion::AliFourPion(const AliFourPion &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),
+    fGeneratorOnly(obj.fGeneratorOnly),
+    fPdensityPairCut(obj.fPdensityPairCut),
+    fTabulatePairs(obj.fTabulatePairs),
+    fRMax(obj.fRMax),
+    ffcSq(obj.ffcSq),
+    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(0),
+    fQLowerCut(obj.fQLowerCut),
+    fNormQcutLow(0),
+    fNormQcutHigh(0),
+    fKupperBound(obj.fKupperBound),
+    fQupperBoundQ2(obj.fQupperBoundQ2),
+    fQupperBoundQ3(obj.fQupperBoundQ3),
+    fQupperBoundQ4(obj.fQupperBoundQ4),
+    fQbinsQ2(obj.fQbinsQ2),
+    fQbinsQ3(obj.fQbinsQ3),
+    fQbinsQ4(obj.fQbinsQ4),
+    fQupperBoundWeights(obj.fQupperBoundWeights),
+    fKstepT(),
+    fKstepY(),
+    fKmeanT(),
+    fKmeanY(),
+    fKmiddleT(),
+    fKmiddleY(),
+    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),
+    fKtIndexL(obj.fKtIndexL),
+    fKtIndexH(obj.fKtIndexH),
+    fQoIndexL(obj.fQoIndexL),
+    fQoIndexH(obj.fQoIndexH),
+    fQsIndexL(obj.fQsIndexL),
+    fQsIndexH(obj.fQsIndexH),
+    fQlIndexL(obj.fQlIndexL),
+    fQlIndexH(obj.fQlIndexH),
+    fDummyB(obj.fDummyB),
+    fKT3transition(obj.fKT3transition),
+    fKT4transition(obj.fKT4transition),
+    fMomResC2(obj.fMomResC2)
+{
+  // Copy Constructor
+  
+  for(Int_t i=0; i<12; i++){
+    fFSIss[i]=obj.fFSIss[i]; 
+    fFSIos[i]=obj.fFSIos[i];
+  }
+  
+  // Initialize fNormWeight and fNormWeightErr to 0
+  for(Int_t i=0; i<3; i++){// Kt iterator
+    for(Int_t j=0; j<10; j++){// Mbin iterator
+      fNormWeight[i][j]=0x0;
+    }
+  }
+  
+
+}
+//________________________________________________________________________
+AliFourPion &AliFourPion::operator=(const AliFourPion &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;
+  fGeneratorOnly = obj.fGeneratorOnly;
+  fPdensityPairCut = obj.fPdensityPairCut;
+  fTabulatePairs = obj.fTabulatePairs;
+  fRMax = obj.fRMax;
+  ffcSq = obj.ffcSq;
+  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;
+  fKupperBound = obj.fKupperBound;
+  fQupperBoundQ2 = obj.fQupperBoundQ2;
+  fQupperBoundQ3 = obj.fQupperBoundQ3;
+  fQupperBoundQ4 = obj.fQupperBoundQ4;
+  fQbinsQ2 = obj.fQbinsQ2;
+  fQbinsQ3 = obj.fQbinsQ3;
+  fQbinsQ4 = obj.fQbinsQ4;
+  fQupperBoundWeights = obj.fQupperBoundWeights;
+  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;
+  fKtIndexL = obj.fKtIndexL;
+  fKtIndexH = obj.fKtIndexH;
+  fQoIndexL = obj.fQoIndexL;
+  fQoIndexH = obj.fQoIndexH;
+  fQsIndexL = obj.fQsIndexL;
+  fQsIndexH = obj.fQsIndexH;
+  fQlIndexL = obj.fQlIndexL;
+  fQlIndexH = obj.fQlIndexH;
+  fDummyB = obj.fDummyB;
+  fKT3transition = obj.fKT3transition;
+  fKT4transition = obj.fKT4transition;
+  fMomResC2 = obj.fMomResC2;
+  
+  for(Int_t i=0; i<12; i++){
+    fFSIss[i]=obj.fFSIss[i]; 
+    fFSIos[i]=obj.fFSIos[i];
+  }
+  for(Int_t i=0; i<3; i++){// Kt iterator
+    for(Int_t j=0; j<10; j++){// Mbin iterator
+      fNormWeight[i][j]=obj.fNormWeight[i][j];
+    }
+  }
+  
+  return (*this);
+}
+//________________________________________________________________________
+AliFourPion::~AliFourPion()
+{
+  // 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;
+  if(fMomResC2) delete fMomResC2;
+  
+  //
+  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 term=0; term<2; term++){
+           
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2;
+           
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal;
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared;
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL;
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW;
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL;
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW;
+           //
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv;
+           if(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW) delete Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW;
+         }// term_2
+         
+         for(Int_t c3=0; c3<2; c3++){
+           for(Int_t term=0; term<5; term++){
+               
+             if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fNorm3;
+             if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3;
+             //
+             if(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm) delete Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm;
+               
+           }// term_3
+         }//c3
+       }//c2
+      }//c1
+      for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
+       for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
+         if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD;
+         if(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD) delete KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD;
+       }
+      }
+      
+    }// ED
+  }// Mbin
+  
+   
+  for(Int_t i=0; i<12; i++){
+    if(fFSIss[i]) delete fFSIss[i]; 
+    if(fFSIos[i]) delete fFSIos[i];
+  }
+  for(Int_t i=0; i<3; i++){// Kt iterator
+    for(Int_t j=0; j<10; j++){// Mbin iterator
+      if(fNormWeight[i][j]) delete fNormWeight[i][j];
+    }
+  }
+  
+}
+//________________________________________________________________________
+void AliFourPion::ParInit()
+{
+  cout<<"AliFourPion MyInit() call"<<endl;
+  cout<<"lego:"<<fLEGO<<"  MCcase:"<<fMCcase<<"  PbPbcase:"<<fPbPbcase<<"  TabulatePairs:"<<fTabulatePairs<<"  GenSignal:"<<fGenerateSignal<<"  CentLow:"<<fCentBinLowLimit<<"  CentHigh:"<<fCentBinHighLimit<<"  RMax:"<<fRMax<<"  fc^2:"<<ffcSq<<"  FB:"<<fFilterBit<<"  MaxChi2/NDF:"<<fMaxChi2NDF<<"  MinTPCncls:"<<fMinTPCncls<<"  MinPairSepEta:"<<fMinSepPairEta<<"  MinPairSepPhi:"<<fMinSepPairPhi<<"  NsigTPC:"<<fSigmaCutTPC<<"  NsigTOF:"<<fSigmaCutTOF<<endl;
+
+  fRandomNumber = new TRandom3();
+  fRandomNumber->SetSeed(0);
+    
+  //
+  fEventCounter=0;
+  fEventsToMix=3;
+  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]=2, fMultLimits[2]=4, fMultLimits[3]=6, fMultLimits[4]=8, fMultLimits[5]=10;
+  fMultLimits[6]=12, fMultLimits[7]=14, fMultLimits[8]=16, fMultLimits[9]=18, fMultLimits[10]=20, fMultLimits[11]=150;
+  
+    
+  
+  if(fPbPbcase) {// PbPb
+    fMultLimit=kMultLimitPbPb;
+    fMbins=fCentBins;
+    fQcut=0.1;
+    fNormQcutLow = 0.15;
+    fNormQcutHigh = 0.175;
+  }else {// pp
+    fMultLimit=kMultLimitpp; 
+    fMbins=kMultBinspp; 
+    fQcut=0.6;
+    fNormQcutLow = 1.0;
+    fNormQcutHigh = 1.5;
+  }
+  
+  fQLowerCut = 0.005;// was 0.005
+  fKupperBound = 1.0;
+  //
+  fKstepY[0] = 1.6;
+  fKmeanY[0] = 0;// central y
+  fKmiddleY[0] = 0;
+
+  // 4x1 (Kt: 0-0.25, 0.25-0.35, 0.35-0.45, 0.45-1.0)
+  if(fKbinsT==4){
+    fKstepT[0] = 0.25; fKstepT[1] = 0.1; fKstepT[2] = 0.1; fKstepT[3] = 0.55;
+    fKmeanT[0] = 0.212; fKmeanT[1] = 0.299; fKmeanT[2] = 0.398; fKmeanT[3] = 0.576;
+    fKmiddleT[0] = 0.125; fKmiddleT[1] = 0.3; fKmiddleT[2] = 0.4; fKmiddleT[3] = 0.725;
+  }
+  // 3x1 (Kt: 0-0.3, 0.3-0.45, 0.45-1.0)
+  if(fKbinsT==3){
+    fKstepT[0] = 0.3; fKstepT[1] = 0.15; fKstepT[2] = 0.55;
+    fKmeanT[0] = 0.240; fKmeanT[1] = 0.369; fKmeanT[2] = 0.576;
+    fKmiddleT[0] = 0.15; fKmiddleT[1] = 0.375; fKmiddleT[2] = 0.725;
+  }
+  // 2x1 (Kt: 0-0.35, 0.35-1.0)
+  if(fKbinsT==2){
+    fKstepT[0] = 0.35; fKstepT[1] = 0.65;
+    fKmeanT[0] = 0.264; fKmeanT[1] = 0.500;
+    fKmiddleT[0] = 0.175; fKmiddleT[1] = 0.675;
+  }
+  //
+  fQupperBoundWeights = 0.2;
+  fQupperBoundQ2 = 2.0;
+  fQupperBoundQ3 = 0.5;
+  fQupperBoundQ4 = 0.5;
+  fQbinsQ2 = fQupperBoundQ2/0.005;
+  fQbinsQ3 = fQupperBoundQ3/0.005;
+  fQbinsQ4 = fQupperBoundQ4/0.005;
+  fQstepWeights = fQupperBoundWeights/Float_t(kQbinsWeights);
+  for(Int_t i=0; i<kQbinsWeights; i++) {fQmean[i]=(i+0.5)*fQstepWeights;}
+  //
+  fDampStart = 0.5;// was 0.3, then 0.5
+  fDampStep = 0.02;
+  
+  //
+  fKT3transition = 0.3;
+  fKT4transition = 0.3;
+  
+  fEC = new AliFourPionEventCollection **[fZvertexBins];
+  for(UShort_t i=0; i<fZvertexBins; i++){
+    
+    fEC[i] = new AliFourPionEventCollection *[fMbins];
+
+    for(UShort_t j=0; j<fMbins; j++){
+      
+      fEC[i][j] = new AliFourPionEventCollection(fEventsToMix+1, fMultLimit, kMCarrayLimit, fMCcase);
+    }
+  }
+  
+  
+  fTempStruct = new AliFourPionTrackStruct[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
+    if(!fTabulatePairs) SetWeightArrays(fLEGO);// Set Weight Array
+    if(!fMCcase && !fTabulatePairs) SetMomResCorrections(fLEGO);// Read Momentum resolution file
+  }
+  
+  /////////////////////////////////////////////
+  /////////////////////////////////////////////
+  
+}
+//________________________________________________________________________
+void AliFourPion::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 *fMultDist0 = new TH1F("fMultDist0","Multiplicity Distribution",fMultLimit,-.5,fMultLimit-.5);
+  fMultDist0->GetXaxis()->SetTitle("Multiplicity");
+  fOutputList->Add(fMultDist0);
+
+  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);
+
+
+
+  TH2D *fResonanceOSPairs = new TH2D("fResonanceOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
+  if(fMCcase) fOutputList->Add(fResonanceOSPairs);
+  TH2D *fAllOSPairs = new TH2D("fAllOSPairs","",fMbins,.5,fMbins+.5, 1000,0,2);
+  if(fMCcase) fOutputList->Add(fAllOSPairs);
+  
+  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);
+
+  //
+  TH1D *fMuonParents = new TH1D("fMuonParents","",500,0.5,500.5);
+  if(fMCcase) fOutputList->Add(fMuonParents);
+  TH1D *fSecondaryMuonParents = new TH1D("fSecondaryMuonParents","",500,0.5,500.5);
+  if(fMCcase) fOutputList->Add(fSecondaryMuonParents);
+  TH3D *fMuonPionDeltaQinv = new TH3D("fMuonPionDeltaQinv","",2,-0.5,1.5, 20,0,1, 100,-0.2,0.2);
+  if(fMCcase) fOutputList->Add(fMuonPionDeltaQinv);
+  TH1D *fPionCandidates = new TH1D("fPionCandidates","",500,0.5,500.5);
+  if(fMCcase) fOutputList->Add(fPionCandidates);
+  //  
+  
+
+  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);
+
+
+  TH1D *fTPNRejects3pion1 = new TH1D("fTPNRejects3pion1","",fQbinsQ3,0,fQupperBoundQ3);
+  fOutputList->Add(fTPNRejects3pion1);
+  TH1D *fTPNRejects3pion2 = new TH1D("fTPNRejects3pion2","",fQbinsQ3,0,fQupperBoundQ3);
+  fOutputList->Add(fTPNRejects3pion2);
+  TH1D *fTPNRejects4pion1 = new TH1D("fTPNRejects4pion1","",fQbinsQ4,0,fQupperBoundQ4);
+  fOutputList->Add(fTPNRejects4pion1);
+
+  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);
+  TH3D *fKT4DistTerm1 = new TH3D("fKT4DistTerm1","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  TH3D *fKT4DistTerm12 = new TH3D("fKT4DistTerm12","",fMbins,.5,fMbins+.5, 20,0,1, 20,0,0.2);
+  fOutputList->Add(fKT4DistTerm1);
+  fOutputList->Add(fKT4DistTerm12);
+
+
+  TProfile2D *fKT3AvgpT = new TProfile2D("fKT3AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
+  fOutputList->Add(fKT3AvgpT);
+  TProfile2D *fKT4AvgpT = new TProfile2D("fKT4AvgpT","",fMbins,.5,fMbins+.5, 2,-0.5,1.5, 0.,1.0,"");
+  fOutputList->Add(fKT4AvgpT);
+
+
+  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);
+
+
+  if(fPdensityPairCut){
+    
+    for(Int_t mb=0; mb<fMbins; mb++){
+      if((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 term=0; term<2; term++){
+               
+             TString *nameEx2 = new TString("TwoParticle_Charge1_");
+             *nameEx2 += c1;
+             nameEx2->Append("_Charge2_");
+             *nameEx2 += c2;
+             nameEx2->Append("_M_");
+             *nameEx2 += mb;
+             nameEx2->Append("_ED_");
+             *nameEx2 += edB;
+             nameEx2->Append("_Term_");
+             *nameEx2 += term+1;
+             
+             if( (c1+c2)==1 ) {if(c1!=0) continue;}// skip degenerate histogram
+             
+             
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2 = new TH2D(nameEx2->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2);
+             TString *nameEx2QW=new TString(nameEx2->Data());
+             nameEx2QW->Append("_QW");
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW = new TH2D(nameEx2QW->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fTerms2QW);
+             TString *nameAvgP=new TString(nameEx2->Data());
+             nameAvgP->Append("_AvgP");
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP = new TProfile2D(nameAvgP->Data(),"",10,0,1, fQbinsQ2,0,fQupperBoundQ2, 0.,1.0,"");
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fAvgP);
+             
+             TString *nameUnitMult=new TString(nameEx2->Data());
+             nameUnitMult->Append("_UnitMult");
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin = new TH2D(nameUnitMult->Data(),"Two Particle Distribution",21,0.5,21.5, fQbinsQ2,0,fQupperBoundQ2);
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fUnitMultBin);
+             
+             if(fMCcase){
+               // Momentum resolution histos
+               TString *nameIdeal = new TString(nameEx2->Data());
+               nameIdeal->Append("_Ideal");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal = new TH2D(nameIdeal->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
+               if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fIdeal);
+               TString *nameSmeared = new TString(nameEx2->Data());
+               nameSmeared->Append("_Smeared");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared = new TH2D(nameSmeared->Data(),"Two Particle Distribution",11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
+               if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fSmeared);
+               //
+               // Muon correction histos
+               TString *nameMuonIdeal=new TString(nameEx2->Data());
+               nameMuonIdeal->Append("_MuonIdeal");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
+               if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonIdeal);
+               TString *nameMuonSmeared=new TString(nameEx2->Data());
+               nameMuonSmeared->Append("_MuonSmeared");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
+               if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonSmeared);
+               //
+               TString *nameMuonPionK2=new TString(nameEx2->Data());
+               nameMuonPionK2->Append("_MuonPionK2");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2 = new TH2D(nameMuonPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
+               if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMuonPionK2);
+               //
+               TString *namePionPionK2=new TString(nameEx2->Data());
+               namePionPionK2->Append("_PionPionK2");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2 = new TH2D(namePionPionK2->Data(),"", 11,0.5,11.5, fQbinsQ2,0,fQupperBoundQ2);
+               if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPionPionK2);
+               //
+               //
+               TString *nameEx2MC=new TString(nameEx2->Data());
+               nameEx2MC->Append("_MCqinv");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv = new TH1D(nameEx2MC->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
+               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinv);
+               TString *nameEx2MCQW=new TString(nameEx2->Data());
+               nameEx2MCQW->Append("_MCqinvQW");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW = new TH1D(nameEx2MCQW->Data(),"", fQbinsQ2,0,fQupperBoundQ2);
+               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fMCqinvQW);
+               //
+               TString *nameEx2PIDpurityDen=new TString(nameEx2->Data());
+               nameEx2PIDpurityDen->Append("_PIDpurityDen");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen = new TH2D(nameEx2PIDpurityDen->Data(),"Two Particle Distribution",20,0,1, fQbinsQ2,0,fQupperBoundQ2);
+               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityDen);
+               TString *nameEx2PIDpurityNum=new TString(nameEx2->Data());
+               nameEx2PIDpurityNum->Append("_PIDpurityNum");
+               Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum = new TH3D(nameEx2PIDpurityNum->Data(),"Two Particle Distribution",16,0.5,16.5, 20,0,1, fQbinsQ2,0,fQupperBoundQ2);
+               fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].fPIDpurityNum);
+             }
+             TString *nameEx2OSLB1 = new TString(nameEx2->Data()); 
+             nameEx2OSLB1->Append("_osl_b1");
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSL);
+             nameEx2OSLB1->Append("_QW");
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW = new TH3D(nameEx2OSLB1->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[0].fTerms2OSLQW);
+             //
+             TString *nameEx2OSLB2 = new TString(nameEx2->Data()); 
+             nameEx2OSLB2->Append("_osl_b2");
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSL);
+             nameEx2OSLB2->Append("_QW");
+             Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW = new TH3D(nameEx2OSLB2->Data(),"Two Particle Distribution",kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
+             fOutputList->Add(Charge1[c1].Charge2[c2].MB[mb].EDB[edB].TwoPT[term].OSL_ktbin[1].fTerms2OSLQW);
+             
+           }// term_2
+           
+           
+           
+           // skip 3-particle if Tabulate6DPairs is true
+           if(fTabulatePairs) continue;
+           
+           for(Int_t c3=0; c3<2; c3++){
+             for(Int_t term=0; term<5; term++){
+                               
+               TString *namePC3 = new TString("ThreeParticle_Charge1_");
+               *namePC3 += c1;
+               namePC3->Append("_Charge2_");
+               *namePC3 += c2;
+               namePC3->Append("_Charge3_");
+               *namePC3 += c3;
+               namePC3->Append("_M_");
+               *namePC3 += mb;
+               namePC3->Append("_ED_");
+               *namePC3 += edB;
+               namePC3->Append("_Term_");
+               *namePC3 += term+1;
+               
+               ///////////////////////////////////////
+               // skip degenerate histograms
+               if( (c1+c2+c3)==1) {if(c3!=1) continue;}
+               if( (c1+c2+c3)==2) {if(c1!=0) continue;}
+               /////////////////////////////////////////
+               
+                               
+               TString *nameNorm=new TString(namePC3->Data());
+               nameNorm->Append("_Norm");
+               Charge1[c1].Charge2[c2].Charge3[c3].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].MB[mb].EDB[edB].ThreePT[term].fNorm3);
+               //
+               
+               TString *name1DQ=new TString(namePC3->Data());
+               name1DQ->Append("_1D");
+               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3 = new TH1D(name1DQ->Data(),"", fQbinsQ3,0,fQupperBoundQ3);
+               fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTerms3);
+               //
+               TString *nameKfactor=new TString(namePC3->Data());
+               nameKfactor->Append("_Kfactor");
+               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,100, "");
+               fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fKfactor);
+               //
+               TString *nameMeanQinv=new TString(namePC3->Data());
+               nameMeanQinv->Append("_MeanQinv");
+               Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv = new TProfile(nameMeanQinv->Data(),"", fQbinsQ3,0,fQupperBoundQ3, 0,.2, "");
+               fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMeanQinv);
+               
+               if(fMCcase==kTRUE){
+                 // Momentum resolution correction histos
+                 TString *nameMomResIdeal=new TString(namePC3->Data());
+                 nameMomResIdeal->Append("_Ideal");
+                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
+                 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fIdeal);
+                 TString *nameMomResSmeared=new TString(namePC3->Data());
+                 nameMomResSmeared->Append("_Smeared");
+                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
+                 if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fSmeared);
+                 // Muon correction histos
+                 TString *nameMuonIdeal=new TString(namePC3->Data());
+                 nameMuonIdeal->Append("_MuonIdeal");
+                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
+                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonIdeal);
+                 TString *nameMuonSmeared=new TString(namePC3->Data());
+                 nameMuonSmeared->Append("_MuonSmeared");
+                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
+                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonSmeared);
+                 //
+                 TString *nameMuonPionK3=new TString(namePC3->Data());
+                 nameMuonPionK3->Append("_MuonPionK3");
+                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3 = new TH2D(nameMuonPionK3->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
+                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fMuonPionK3);
+                 //
+                 TString *namePionPionK3=new TString(namePC3->Data());
+                 namePionPionK3->Append("_PionPionK3");
+                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3 = new TH2D(namePionPionK3->Data(),"", 11,0.5,11.5, fQbinsQ3,0,fQupperBoundQ3);
+                 if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fPionPionK3);
+                 
+               }// MCcase
+               //
+               if(c1==c2 && c1==c3 && term==4 ){
+                 TString *nameTwoPartNorm=new TString(namePC3->Data());
+                 nameTwoPartNorm->Append("_TwoPartNorm");
+                 Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ3,0,fQupperBoundQ3);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].MB[mb].EDB[edB].ThreePT[term].fTwoPartNorm);
+               }// term=4
+               
+             }// term_3
+             
+             for(Int_t c4=0; c4<2; c4++){
+               for(Int_t term=0; term<12; term++){
+                               
+                 TString *namePC4 = new TString("FourParticle_Charge1_");
+                 *namePC4 += c1;
+                 namePC4->Append("_Charge2_");
+                 *namePC4 += c2;
+                 namePC4->Append("_Charge3_");
+                 *namePC4 += c3;
+                 namePC4->Append("_Charge4_");
+                 *namePC4 += c4;
+                 namePC4->Append("_M_");
+                 *namePC4 += mb;
+                 namePC4->Append("_ED_");
+                 *namePC4 += edB;
+                 namePC4->Append("_Term_");
+                 *namePC4 += term+1;
+                 
+                 ///////////////////////////////////////
+                 // skip degenerate histograms
+                 if( (c1+c2+c3+c4)==1) {if(c4!=1) continue;}
+                 if( (c1+c2+c3+c4)==2) {if(c3+c4!=2) continue;}
+                 if( (c1+c2+c3+c4)==3) {if(c1!=0) continue;}
+                 /////////////////////////////////////////
+                
+                 TString *nameNorm=new TString(namePC4->Data());
+                 nameNorm->Append("_Norm");
+                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4 = new TH1D(nameNorm->Data(),"Norm",1,-0.5,0.5);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fNorm4);
+                 //
+                 TString *name1DQ=new TString(namePC4->Data());
+                 name1DQ->Append("_1D");
+                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4 = new TH1D(name1DQ->Data(),"", fQbinsQ4,0,fQupperBoundQ4);
+                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTerms4);
+                 //
+                 TString *nameKfactor=new TString(namePC4->Data());
+                 nameKfactor->Append("_Kfactor");
+                 Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor = new TProfile(nameKfactor->Data(),"", fQbinsQ4,0,fQupperBoundQ4, 0,100, "");
+                 fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fKfactor);
+                 //
+                 if(c1==c2 && c1==c3 && c1==c4 && term==11 ){
+                   TString *nameTwoPartNorm=new TString(namePC4->Data());
+                   nameTwoPartNorm->Append("_TwoPartNorm");
+                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm = new TH2D(nameTwoPartNorm->Data(),"", kDENtypes,0.5,kDENtypes+0.5, fQbinsQ4,0,fQupperBoundQ4);
+                   fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fTwoPartNorm);
+                 }
+
+                 if(fMCcase==kTRUE){
+                   // Momentum resolution correction histos
+                   TString *nameMomResIdeal=new TString(namePC4->Data());
+                   nameMomResIdeal->Append("_Ideal");
+                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal = new TH2D(nameMomResIdeal->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
+                   if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fIdeal);
+                   TString *nameMomResSmeared=new TString(namePC4->Data());
+                   nameMomResSmeared->Append("_Smeared");
+                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared = new TH2D(nameMomResSmeared->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
+                   if(mb==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fSmeared);
+                   // Muon correction histos
+                   TString *nameMuonIdeal=new TString(namePC4->Data());
+                   nameMuonIdeal->Append("_MuonIdeal");
+                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal = new TH2D(nameMuonIdeal->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
+                   if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonIdeal);
+                   TString *nameMuonSmeared=new TString(namePC4->Data());
+                   nameMuonSmeared->Append("_MuonSmeared");
+                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared = new TH2D(nameMuonSmeared->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
+                   if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonSmeared);
+                   //
+                   TString *nameMuonPionK4=new TString(namePC4->Data());
+                   nameMuonPionK4->Append("_MuonPionK4");
+                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4 = new TH2D(nameMuonPionK4->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
+                   if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fMuonPionK4);
+                   //
+                   TString *namePionPionK4=new TString(namePC4->Data());
+                   namePionPionK4->Append("_PionPionK4");
+                   Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4 = new TH2D(namePionPionK4->Data(),"", 11,0.5,11.5, fQbinsQ4,0,fQupperBoundQ4);
+                   if(mb==0 && edB==0) fOutputList->Add(Charge1[c1].Charge2[c2].Charge3[c3].Charge4[c4].MB[mb].EDB[edB].FourPT[term].fPionPionK4);
+                   
+                 }// MCcase
+                 
+                 
+               }
+             }
+
+           }//c3
+         }//c2
+       }//c1
+      }// ED
+    }// mbin
+  }// Pdensity Method
+
+  
+  if(fTabulatePairs){
+    
+    for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
+      for(Int_t yKbin=0; yKbin<fKbinsY; yKbin++){
+       for(Int_t mb=0; mb<fMbins; mb++){
+         for(Int_t edB=0; edB<fEDbins; edB++){
+           
+           TString *nameNum = new TString("TPN_num_Kt_");
+           *nameNum += tKbin;
+           nameNum->Append("_Ky_");
+           *nameNum += yKbin;
+           nameNum->Append("_M_");
+           *nameNum += mb;
+           nameNum->Append("_ED_");
+           *nameNum += edB;
+           
+           TString *nameDen = new TString("TPN_den_Kt_");
+           *nameDen += tKbin;
+           nameDen->Append("_Ky_");
+           *nameDen += yKbin;
+           nameDen->Append("_M_");
+           *nameDen += mb;
+           nameDen->Append("_ED_");
+           *nameDen += edB;
+           
+           if(edB==0){
+             KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD = new TH3D(nameNum->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
+             fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[0].fTerms2ThreeD);
+             
+             KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD = new TH3D(nameDen->Data(),"", kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights, kQbinsWeights,0,fQupperBoundWeights);
+             fOutputList->Add(KT[tKbin].KY[yKbin].MB[mb].EDB[edB].TwoPT[1].fTerms2ThreeD);
+           }   
+         
+         }
+       }
+      }
+    }
+    
+  }
+  
+  
+  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);
+  TH2D *fQ2Res = new TH2D("fQ2Res","",20,0,1, 200,-.2,.2);
+  fOutputList->Add(fQ2Res);
+  TH2D *fQ3Res = new TH2D("fQ3Res","",20,0,1, 200,-.3,.3);
+  fOutputList->Add(fQ3Res);
+  TH2D *fQ4Res = new TH2D("fQ4Res","",20,0,1, 200,-.4,.4);
+  fOutputList->Add(fQ4Res);
+  
+  TH2D *DistQinv4pion = new TH2D("DistQinv4pion","",6,0.5,6.5, 20,0,0.1);
+  fOutputList->Add(DistQinv4pion);
+  TH2D *DistQinvMC4pion = new TH2D("DistQinvMC4pion","",6,0.5,6.5, 20,0,0.1);
+  if(fMCcase) fOutputList->Add(DistQinvMC4pion);
+
+  TH2D *fAvgQ12VersusQ3 = new TH2D("fAvgQ12VersusQ3","",10,0,0.1, 20,0,0.1);
+  fOutputList->Add(fAvgQ12VersusQ3);
+  TH2D *fAvgQ13VersusQ3 = new TH2D("fAvgQ13VersusQ3","",10,0,0.1, 20,0,0.1);
+  fOutputList->Add(fAvgQ13VersusQ3);
+  TH2D *fAvgQ23VersusQ3 = new TH2D("fAvgQ23VersusQ3","",10,0,0.1, 20,0,0.1);
+  fOutputList->Add(fAvgQ23VersusQ3);
+
+  ////////////////////////////////////
+  ///////////////////////////////////  
+  
+  PostData(1, fOutputList);
+}
+
+//________________________________________________________________________
+void AliFourPion::UserExec(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(fAOD->GetRunNumber() >= 136851 && fAOD->GetRunNumber() <= 139517){// 10h data
+  Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+  if(!isSelected1 && !fMCcase) {return;}
+  }else if(fAOD->GetRunNumber() >= 167693 && fAOD->GetRunNumber() <= 170593){// 11h data
+    Bool_t isSelected1 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);
+    Bool_t isSelected2 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);
+    if(!isSelected1 && !isSelected2 && !fMCcase) {return;}
+  }else {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;
+  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;
+      }
+    }
+  }
+  
+
+  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;
+  
+  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;
+    }
+    
+
+    ((TH1F*)fOutputList->FindObject("fMultDist0"))->Fill(fAOD->GetNumberOfTracks());
+
+    // Pile-up rejection
+    AliAnalysisUtils *AnaUtil=new AliAnalysisUtils();
+    if(!fPbPbcase) AnaUtil->SetUseMVPlpSelection(kTRUE);// use Multi-Vertex tool for pp and pPb
+    else AnaUtil->SetUseMVPlpSelection(kFALSE);
+    Bool_t pileUpCase=AnaUtil->IsPileUpEvent(fAOD); 
+    if(pileUpCase) return;
+    
+    ////////////////////////////////
+    // 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(!fMCcase && 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;
+      }
+    }
+
+   
+       
+    /////////////////////////////
+    // 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
+      ((TH2D*)fOutputList->FindObject("fTrackChi2NDF"))->Fill(centralityPercentile, aodtrack->Chi2perNDF());
+      ((TH2D*)fOutputList->FindObject("fTrackTPCncls"))->Fill(centralityPercentile, aodtrack->GetTPCncls());
+      if(aodtrack->GetTPCNcls() < fMinTPCncls) continue;// TPC nCluster cut
+      if(aodtrack->Chi2perNDF() > fMaxChi2NDF) continue;
+
+      if(fFilterBit != 7){
+       Bool_t goodTrackOtherFB = kFALSE;
+       for (Int_t j = 0; j < fAOD->GetNumberOfTracks(); j++) {
+         AliAODTrack* aodtrack2 = fAOD->GetTrack(randomIndex[j]);
+         if(!aodtrack2) continue;
+         if(!aodtrack2->TestFilterBit(BIT(fFilterBit))) continue;
+         
+         if(-(aodtrack->GetID()+1)==aodtrack2->GetID()) {goodTrackOtherFB=kTRUE; break;}
+         
+       }
+       if(!goodTrackOtherFB) continue;
+      }
+      
+      
+      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);
+      
+        
+      Double_t dca2[2]={0};
+      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));
+      Double_t 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
+     
+     
+      // 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};
+
+    
+      Bool_t DoPIDWorkAround=kTRUE;
+      //if(fFilterBit == 7) DoPIDWorkAround=kTRUE;
+      if(fMCcase && !fPbPbcase) DoPIDWorkAround=kFALSE;
+      if(DoPIDWorkAround==kFALSE && fabs(fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion)) < 900) {
+       nSigmaTPC[0]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kElectron);
+       nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kMuon);
+       nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kPion);
+       nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kKaon);
+       nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodtrack,AliPID::kProton);
+       //
+       nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kElectron);
+       nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kMuon);
+       nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kPion);
+       nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodtrack,AliPID::kKaon);
+       nSigmaTOF[4]=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]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kElectron);
+         nSigmaTPC[1]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kMuon);
+         nSigmaTPC[2]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kPion);
+         nSigmaTPC[3]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kKaon);
+         nSigmaTPC[4]=fPIDResponse->NumberOfSigmasTPC(aodTrack2,AliPID::kProton);
+         //
+         nSigmaTOF[0]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kElectron);
+         nSigmaTOF[1]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kMuon);
+         nSigmaTOF[2]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kPion);
+         nSigmaTOF[3]=fPIDResponse->NumberOfSigmasTOF(aodTrack2,AliPID::kKaon);
+         nSigmaTOF[4]=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(fabs(nSigmaTOF[0])<fSigmaCutTOF) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+       if(fabs(nSigmaTOF[2])<fSigmaCutTOF) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+       if(fabs(nSigmaTOF[3])<fSigmaCutTOF) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+       if(fabs(nSigmaTOF[4])<fSigmaCutTOF) fTempStruct[myTracks].fProton = kTRUE;// Proton candidate
+      }else {// TPC info instead
+       if(fabs(nSigmaTPC[0])<fSigmaCutTPC) fTempStruct[myTracks].fElectron = kTRUE;// Electron candidate
+       if(fabs(nSigmaTPC[2])<fSigmaCutTPC) fTempStruct[myTracks].fPion = kTRUE;// Pion candidate
+       if(fabs(nSigmaTPC[3])<fSigmaCutTPC) fTempStruct[myTracks].fKaon = kTRUE;// Kaon candidate
+       if(fabs(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) continue;// only pions
+      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;// superfluous
+      ////////////////////////
+      //if(fTempStruct[myTracks].fProton && fTempStruct[myTracks].fMom < 0.25) continue;//extra cut for protons// superfluous
+
+
+   
+      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());
+      
+     
+    
+      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;
+      }
+      
+           
+
+      if(aodtrack->Charge() > 0) positiveTracks++;
+      else negativeTracks++;
+      
+      if(fTempStruct[myTracks].fPion) pionCount++;
+      if(fTempStruct[myTracks].fKaon) kaonCount++;
+      if(fTempStruct[myTracks].fProton) protonCount++;
+
+      myTracks++;
+      
+      if(fMCcase){// muon mothers
+       AliAODMCParticle *tempMCTrack=(AliAODMCParticle*)mcArray->At(abs(aodtrack->GetLabel()));
+       if(abs(tempMCTrack->GetPdgCode())==13 && tempMCTrack->GetMother()>0){// muons
+         AliAODMCParticle *parent=(AliAODMCParticle*)mcArray->At(tempMCTrack->GetMother());
+         if(parent->IsPhysicalPrimary()){
+           ((TH1D*)fOutputList->FindObject("fMuonParents"))->Fill(abs(parent->GetPdgCode()));
+         }else ((TH1D*)fOutputList->FindObject("fSecondaryMuonParents"))->Fill(abs(parent->GetPdgCode()));
+       }
+       ((TH1D*)fOutputList->FindObject("fPionCandidates"))->Fill(abs(tempMCTrack->GetPdgCode()));
+      }
+    }
+    //cout<<"kinkcount = "<<kinkcount<<"   pionkinks = "<<pionkinks<<"   primarypionkinks = "<<primarypionkinks<<endl;
+  }else {// ESD tracks
+    cout<<"ESDs not supported currently"<<endl;
+    return;
+  }
+  
+  // Generator info only
+  if(fMCcase && fGeneratorOnly){
+    myTracks=0; pionCount=0; kaonCount=0; protonCount=0;// reset track counters
+    for(Int_t mctrackN=0; mctrackN<mcArray->GetEntriesFast(); mctrackN++){
+      if(myTracks >= fMultLimit) {cout<<"More tracks than Track Limit"<<endl; return;}
+      if(myTracks >= 1300) continue;// additional cut to limit high mult events which exceed pair # limits
+      
+      AliAODMCParticle *mcParticle = (AliAODMCParticle*)mcArray->At(mctrackN);
+      if(!mcParticle) continue;
+      if(fabs(mcParticle->Eta())>0.8) continue;
+      if(mcParticle->Charge()!=-3 && mcParticle->Charge()!=+3) continue;// x3 by convention
+      if(mcParticle->Pt() < 0.16 || mcParticle->Pt() > 1.0) continue;
+      if(!mcParticle->IsPrimary()) continue;
+      if(!mcParticle->IsPhysicalPrimary()) continue;
+      if(abs(mcParticle->GetPdgCode())!=211) continue;
+      
+      fTempStruct[myTracks].fP[0] = mcParticle->Px();
+      fTempStruct[myTracks].fP[1] = mcParticle->Py();
+      fTempStruct[myTracks].fP[2] = mcParticle->Pz();
+      fTempStruct[myTracks].fX[0] = 0.; fTempStruct[myTracks].fX[1] = 0.; fTempStruct[myTracks].fX[2] = 0.;
+      
+      fTempStruct[myTracks].fId = myTracks;// use my track counter 
+      fTempStruct[myTracks].fLabel = mctrackN;
+      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 = mcParticle->Eta();
+      fTempStruct[myTracks].fCharge = int(mcParticle->Charge()/3.);
+      fTempStruct[myTracks].fDCAXY = 0.;
+      fTempStruct[myTracks].fDCAZ = 0.;
+      fTempStruct[myTracks].fDCA = 0.;
+      fTempStruct[myTracks].fPion = kTRUE;
+      fTempStruct[myTracks].fEaccepted = sqrt(pow(fTempStruct[myTracks].fMom,2) + pow(fTrueMassPi,2)); 
+      fTempStruct[myTracks].fKey = 1;
+      
+      myTracks++;
+      pionCount++;
+    }
+  }
+  
+  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 < 4) {cout<<"Less than 4 tracks. Skipping Event."<<endl; return;}
+  /////////////////////////////////////////
+
+  ////////////////////////////////
+  ///////////////////////////////
+  // Mbin determination
+  //
+  // Mbin set to Pion Count Only for pp!!!!!!!
+  fMbin=-1;
+  if(!fPbPbcase){
+    for(Int_t i=0; i<kMultBinspp; i++){
+      if( ( pionCount > fMultLimits[i]) && ( pionCount <= fMultLimits[i+1]) ) { fMbin=i; break;}
+      // Mbin 0 has 1 pion
+    }
+  }else{
+    for(Int_t i=0; i<fCentBins; i++){
+      if( (centralityPercentile >= cStart+i*cStep) && (centralityPercentile < cStart+(i+1)*cStep) ){
+       fMbin=i;// 0 = most central
+       break;
+      }
+    }
+  }
+
+  if(fMbin==-1) {cout<<"Bad Mbin+++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; return;}
+  
+  ///////////////////
+  // can only be called after fMbin has been set
+  // Radius parameter only matters for Monte-Carlo data
+  SetFSIindex(fRMax);
+  ///////////////////  
+  
+  Int_t rBinForTPNMomRes = 10;
+  if(fMbin==0) {rBinForTPNMomRes=10;}// 10 fm with EW (fRMax should be 11 for normal running)
+  else if(fMbin==1) {rBinForTPNMomRes=9;}
+  else if(fMbin<=3) {rBinForTPNMomRes=8;}
+  else if(fMbin<=5) {rBinForTPNMomRes=7;}
+  else {rBinForTPNMomRes=6;}
+
+  //////////////////////////////////////////////////
+  fEDbin=0;// Extra Dimension bin (Kt, (Kt-Psi),....)
+  //////////////////////////////////////////////////
+  
+
+  
+  ((TH1F*)fOutputList->FindObject("fEvents1"))->Fill(fMbin+1);
+  ((TProfile*)fOutputList->FindObject("fAvgMult"))->Fill(fMbin+1., pionCount);
+
+  ////////////////////////////////////
+  // 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));
+       (fEvt)->fMCtracks[i].fPdgCode = tempMCTrack->GetPdgCode();
+       (fEvt)->fMCtracks[i].fMotherLabel = tempMCTrack->GetMother();
+      }        
+    }
+  }
+    
+  
+  
+  Float_t qinv12=0, qinv13=0, qinv14=0, qinv23=0, qinv24=0, qinv34=0;
+  Float_t qout=0, qside=0, qlong=0;
+  Float_t kT12=0;
+  Float_t q3=0, q3MC=0;
+  Float_t q4=0, q4MC=0;
+  Int_t ch1=0, ch2=0, ch3=0, ch4=0;
+  Int_t bin1=0, bin2=0, bin3=0, bin4=0;
+  Float_t pVect1[4]={0}; 
+  Float_t pVect2[4]={0};
+  Float_t pVect3[4]={0};
+  Float_t pVect4[4]={0};
+  Float_t pVect1MC[4]={0}; 
+  Float_t pVect2MC[4]={0};
+  Float_t pVect3MC[4]={0};
+  Float_t pVect4MC[4]={0};
+  Float_t Pparent1[4]={0};
+  Float_t Pparent2[4]={0};
+  Float_t Pparent3[4]={0};
+  Float_t Pparent4[4]={0};
+  Float_t weight12=0, weight13=0, weight14=0, weight23=0, weight24=0, weight34=0;
+  Float_t weight12Err=0, weight13Err=0, weight14Err=0, weight23Err=0, weight24Err=0, weight34Err=0;
+  Float_t weight12CC=0, weight13CC=0, weight14CC=0, weight23CC=0, weight24CC=0, weight34CC=0;
+  Float_t weightTotal=0;//, weightTotalErr=0;
+  Float_t qinv12MC=0, qinv13MC=0, qinv14MC=0, qinv23MC=0, qinv24MC=0, qinv34MC=0; 
+  Float_t parentQinv12=0, parentQinv13=0, parentQinv14=0, parentQinv23=0, parentQinv24=0, parentQinv34=0;
+  Float_t parentQ3=0;
+  Float_t FSICorr12=0, FSICorr13=0, FSICorr14=0, FSICorr23=0, FSICorr24=0, FSICorr34=0;
+  Bool_t pionParent1=kFALSE, pionParent2=kFALSE, pionParent3=kFALSE, pionParent4=kFALSE;
+  Bool_t FilledMCpair12=kFALSE, FilledMCtriplet123=kFALSE;
+  Bool_t GoodTripletWeights=kFALSE;
+  //
+  AliAODMCParticle *mcParticle1=0x0;
+  AliAODMCParticle *mcParticle2=0x0;
+  
+
+  ////////////////////
+  Int_t PairCount[7]={0};
+  Int_t NormPairCount[7]={0};
+  Int_t KT3index=0, KT4index=0;
+
+    // reset to defaults
+    for(Int_t i=0; i<fMultLimit; i++) {
+      for(Int_t j=0; j<fMultLimit; j++) {
+       LowQPairSwitch_E0E0[i][j] = kFALSE;
+       LowQPairSwitch_E0E1[i][j] = kFALSE;
+       LowQPairSwitch_E0E2[i][j] = kFALSE;
+       LowQPairSwitch_E0E3[i][j] = kFALSE;
+       LowQPairSwitch_E1E2[i][j] = kFALSE;
+       LowQPairSwitch_E1E3[i][j] = kFALSE;
+       LowQPairSwitch_E2E3[i][j] = kFALSE;
+       //
+       NormQPairSwitch_E0E0[i][j] = kFALSE;
+       NormQPairSwitch_E0E1[i][j] = kFALSE;
+       NormQPairSwitch_E0E2[i][j] = kFALSE;
+       NormQPairSwitch_E0E3[i][j] = kFALSE;
+       NormQPairSwitch_E1E2[i][j] = kFALSE;
+       NormQPairSwitch_E1E3[i][j] = kFALSE;
+       NormQPairSwitch_E2E3[i][j] = kFALSE;
+      }
+      
+    
+    }
+
+   
+    //////////////////////////////////////////
+    // make low-q pair storage and normalization-pair storage
+    // 
+    for(Int_t en1=0; en1<=2; en1++){// 1st event number (en1=0 is the same event as current event)
+      for(Int_t en2=en1; en2<=3; en2++){// 2nd event number (en2=0 is the same event as current event)
+       if(en1>0 && en1==en2) continue;
+       for (Int_t i=0; i<(fEvt+en1)->fNtracks; i++) {// 1st particle
+         for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
+           
+           //if(pairCount[en2] >= kPairLimit) {exitCode=kTRUE; continue;}// Too many SE pairs
+           
+           
+           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];
+           ch1 = Int_t(((fEvt+en1)->fTracks[i].fCharge + 1)/2.);
+           ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
+
+           qinv12 = GetQinv(pVect1, pVect2);
+           kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
+           SetFillBins2(ch1, ch2, bin1, bin2);
+             
+           if(qinv12 < fQLowerCut) continue;// remove unwanted low-q pairs (also a type of track splitting/merging cut)
+           if(ch1 == ch2 && !fGeneratorOnly){
+             if(!AcceptPair((fEvt+en1)->fTracks[i], (fEvt+en2)->fTracks[j])) {
+               if(en1==0 && en2==0) ((TH1F*)fOutputList->FindObject("fRejectedPairs"))->Fill(qinv12);
+               continue;
+             }
+           }
+           
+           GetQosl(pVect1, pVect2, qout, qside, qlong);
+           if( (en1+en2==0)) {
+             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2->Fill(kT12, qinv12);
+             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[0].fTerms2QW->Fill(kT12, qinv12, qinv12);
+             // osl frame
+             if((kT12 > 0.2) && (kT12 < 0.3)){
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
+             }
+             if((kT12 > 0.6) && (kT12 < 0.7)){  
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
+             }
+             // unit mult bins
+             if( (fEvt+en1)->fNtracks%100==0){
+               Int_t kTindex=0;
+               if(kT12>0.3) kTindex=1;
+               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[0].fUnitMultBin->Fill(UnitMultBin, qinv12);
+             }
+
+           }
+           if( (en1+en2==1)) {
+             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2->Fill(kT12, qinv12);
+             Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fTerms2QW->Fill(kT12, qinv12, qinv12);
+             // osl frame
+             if((kT12 > 0.2) && (kT12 < 0.3)){  
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[0].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
+             }
+             if((kT12 > 0.6) && (kT12 < 0.7)){  
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSL->Fill(fabs(qout), fabs(qside), fabs(qlong));
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[fEDbin].TwoPT[en2].OSL_ktbin[1].fTerms2OSLQW->Fill(fabs(qout), fabs(qside), fabs(qlong), qinv12);
+             }
+             // unit mult bins
+             if( (fEvt+en1)->fNtracks%100==0){
+               Int_t kTindex=0;
+               if(kT12>0.3) kTindex=1;
+               Int_t UnitMultBin = int((fEvt+en1)->fNtracks / 100.) + 1;
+               Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[kTindex].TwoPT[1].fUnitMultBin->Fill(UnitMultBin, qinv12);
+             }
+           }
+           //////////////////////////////////////////
+           if(fTabulatePairs && en2<=1 && bin1==bin2){
+             Float_t kY = 0;
+             Int_t kTbin=-1, kYbin=-1;
+             
+             for(Int_t kIt=0; kIt<fKbinsT; kIt++) {if(kT12 < (fKmiddleT[kIt] + fKstepT[kIt]/2.)) {kTbin = kIt; break;}} 
+             for(Int_t kIt=0; kIt<fKbinsY; kIt++) {if(kY < (fKmiddleY[kIt] + fKstepY[kIt]/2.)) {kYbin = kIt; break;}}
+             if((kTbin<0) || (kYbin<0)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
+             if((kTbin>=fKbinsT) || (kYbin>=fKbinsY)) {cout<<"problem!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<endl; continue;}
+             if(fGenerateSignal && en2==0) {
+               Int_t chGroup2[2]={ch1,ch2};
+               Float_t WInput = MCWeight(chGroup2, fRMax, 0.7, qinv12, kT12);
+               KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong), WInput);
+               }else KT[kTbin].KY[kYbin].MB[fMbin].EDB[0].TwoPT[en2].fTerms2ThreeD->Fill(fabs(qout), fabs(qside), fabs(qlong));
+           }
+           
+           //////////////////////////////////////////////////////////////////////////////
+
+           if(qinv12 <= fQcut) {
+             if(en1==0 && en2==0) {LowQPairSwitch_E0E0[i][j] = kTRUE; PairCount[0]++;}
+             if(en1==0 && en2==1) {LowQPairSwitch_E0E1[i][j] = kTRUE; PairCount[1]++;}
+             if(en1==0 && en2==2) {LowQPairSwitch_E0E2[i][j] = kTRUE; PairCount[2]++;}
+             if(en1==0 && en2==3) {LowQPairSwitch_E0E3[i][j] = kTRUE; PairCount[3]++;}
+             if(en1==1 && en2==2) {LowQPairSwitch_E1E2[i][j] = kTRUE; PairCount[4]++;}
+             if(en1==1 && en2==3) {LowQPairSwitch_E1E3[i][j] = kTRUE; PairCount[5]++;}
+             if(en1==2 && en2==3) {LowQPairSwitch_E2E3[i][j] = kTRUE; PairCount[6]++;}
+           }
+           if((qinv12 >= fNormQcutLow) && (qinv12 < fNormQcutHigh)) {
+             if(en1==0 && en2==0) {NormQPairSwitch_E0E0[i][j] = kTRUE; NormPairCount[0]++;}
+             if(en1==0 && en2==1) {NormQPairSwitch_E0E1[i][j] = kTRUE; NormPairCount[1]++;}
+             if(en1==0 && en2==2) {NormQPairSwitch_E0E2[i][j] = kTRUE; NormPairCount[2]++;}
+             if(en1==0 && en2==3) {NormQPairSwitch_E0E3[i][j] = kTRUE; NormPairCount[3]++;}
+             if(en1==1 && en2==2) {NormQPairSwitch_E1E2[i][j] = kTRUE; NormPairCount[4]++;}
+             if(en1==1 && en2==3) {NormQPairSwitch_E1E3[i][j] = kTRUE; NormPairCount[5]++;}
+             if(en1==2 && en2==3) {NormQPairSwitch_E2E3[i][j] = kTRUE; NormPairCount[6]++;}
+           }
+           
+         }
+       }
+      }
+    }
+    
+    //cout<<PairCount[0]<<"  "<<PairCount[1]<<"  "<<PairCount[2]<<"  "<<PairCount[3]<<"  "<<PairCount[4]<<"  "<<PairCount[5]<<"  "<<PairCount[6]<<endl;
+    //cout<<NormPairCount[0]<<"  "<<NormPairCount[1]<<"  "<<NormPairCount[2]<<"  "<<NormPairCount[3]<<"  "<<NormPairCount[4]<<"  "<<NormPairCount[5]<<"  "<<NormPairCount[6]<<endl;
+    ///////////////////////////////////////////////////  
+    // Do not use pairs from events with too many pairs
+    
+    ((TH1F*)fOutputList->FindObject("fEvents2"))->Fill(fMbin+1);
+    
+    ///////////////////////////////////////////////////
+
+    
+    if(fTabulatePairs) return;
+
+
+
+    ////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////
+    // Normalization counting of 3- and 4-particle terms
+    for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
+      for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
+       if(en2==0 && en3>1) continue;// not needed config
+       if(en2==1 && en3==en2) continue;// not needed config
+       for(Int_t en4=en3; en4<=3; en4++){// 4th event number
+         if(en3==0 && en4>1) continue;// not needed config
+         if(en3==1 && en4!=2) continue;// not needed configs
+         if(en3==2 && en4!=3) continue;// not needed configs
+         
+         for (Int_t i=0; i<myTracks; i++) {// 1st particle
+           pVect1[1]=(fEvt)->fTracks[i].fP[0];
+           pVect1[2]=(fEvt)->fTracks[i].fP[1];
+           pVect1[3]=(fEvt)->fTracks[i].fP[2];
+           ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
+
+           for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
+             if(en2==0) {if(!NormQPairSwitch_E0E0[i][j]) continue;}
+             else {if(!NormQPairSwitch_E0E1[i][j]) continue;}
+             
+             pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
+             pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
+             pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
+             ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
+
+             for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
+               if(en3==0) {
+                 if(!NormQPairSwitch_E0E0[i][k]) continue;
+                 if(!NormQPairSwitch_E0E0[j][k]) continue;
+               }else if(en3==1){
+                 if(!NormQPairSwitch_E0E1[i][k]) continue;
+                 if(!NormQPairSwitch_E0E1[j][k]) continue;
+               }else{
+                 if(!NormQPairSwitch_E0E2[i][k]) continue;
+                 if(!NormQPairSwitch_E1E2[j][k]) continue;
+               }
+               pVect3[1]=(fEvt+en3)->fTracks[k].fP[0];
+               pVect3[2]=(fEvt+en3)->fTracks[k].fP[1];
+               pVect3[3]=(fEvt+en3)->fTracks[k].fP[2];
+               ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
+               Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
+               SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
+
+               Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
+               if(KT3<=fKT3transition) KT3index=0;
+               else KT3index=1;
+               
+               if(en2==0 && en3==0 && en4==0) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fNorm3->Fill(0);
+               if(en2==1 && en3==2 && en4==3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fNorm3->Fill(0);
+               if(en2==0 && en3==1 && en4==2) {
+                 if(fill2) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fNorm3->Fill(0);
+                 if(fill3) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fNorm3->Fill(0);
+                 if(fill4) Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fNorm3->Fill(0);
+               }
+                 
+               
+               for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
+                 if(en4==0){
+                   if(!NormQPairSwitch_E0E0[i][l]) continue;
+                   if(!NormQPairSwitch_E0E0[j][l]) continue;
+                   if(!NormQPairSwitch_E0E0[k][l]) continue;
+                 }else if(en4==1){
+                   if(!NormQPairSwitch_E0E1[i][l]) continue;
+                   if(!NormQPairSwitch_E0E1[j][l]) continue;
+                   if(!NormQPairSwitch_E0E1[k][l]) continue;
+                 }else if(en4==2){
+                   if(!NormQPairSwitch_E0E2[i][l]) continue;
+                   if(!NormQPairSwitch_E0E2[j][l]) continue;
+                   if(!NormQPairSwitch_E1E2[k][l]) continue;
+                 }else{
+                   if(!NormQPairSwitch_E0E3[i][l]) continue;
+                   if(!NormQPairSwitch_E1E3[j][l]) continue;
+                   if(!NormQPairSwitch_E2E3[k][l]) continue;
+                 }
+                 
+                 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
+                 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
+                 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
+                 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
+                 Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
+                 if(KT4<=fKT4transition) KT4index=0;
+                 else KT4index=1;
+
+                 Bool_t FillTerms[12]={kFALSE};
+                 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, en2+en3+en4, FillTerms);
+                 //
+                 for(int ft=0; ft<12; ft++) {
+                   if(FillTerms[ft]) Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fNorm4->Fill(0.); 
+                 }
+       
+                 
+               }
+             }
+           }
+         }  
+         
+       }
+      }
+    }
+    
+  
+   
+
+    ///////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////
+    //
+    //
+    // Start the Main Correlation Analysis
+    //
+    //
+    ///////////////////////////////////////////////////////////////////////
+
+
+
+    ////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////
+    for(Int_t en2=0; en2<=1; en2++){// 2nd event number (en2=0 is the same event as current event)
+      for(Int_t en3=en2; en3<=2; en3++){// 3rd event number
+       if(en2==0 && en3>1) continue;// not needed config
+       if(en2==1 && en3==en2) continue;// not needed config
+       for(Int_t en4=en3; en4<=3; en4++){// 4th event number
+         if(en3==0 && en4>1) continue;// not needed config
+         if(en3==1 && en4!=2) continue;// not needed configs
+         if(en3==2 && en4!=3) continue;// not needed configs
+         Int_t ENsum=en2+en3+en4;// 0 or 1 or 3 or 6
+         /////////////////////////////////////////////////////////////
+         for (Int_t i=0; i<myTracks; i++) {// 1st particle
+           pVect1[0]=(fEvt)->fTracks[i].fEaccepted;
+           pVect1[1]=(fEvt)->fTracks[i].fP[0];
+           pVect1[2]=(fEvt)->fTracks[i].fP[1];
+           pVect1[3]=(fEvt)->fTracks[i].fP[2];
+           ch1 = Int_t(((fEvt)->fTracks[i].fCharge + 1)/2.);
+
+           /////////////////////////////////////////////////////////////
+           for (Int_t j=i+1; j<(fEvt+en2)->fNtracks; j++) {// 2nd particle
+             if(en2==0) {if(!LowQPairSwitch_E0E0[i][j]) continue;}
+             else {if(!LowQPairSwitch_E0E1[i][j]) continue;}
+             
+             pVect2[0]=(fEvt+en2)->fTracks[j].fEaccepted;
+             pVect2[1]=(fEvt+en2)->fTracks[j].fP[0];
+             pVect2[2]=(fEvt+en2)->fTracks[j].fP[1];
+             pVect2[3]=(fEvt+en2)->fTracks[j].fP[2];
+             ch2 = Int_t(((fEvt+en2)->fTracks[j].fCharge + 1)/2.);
+             qinv12 = GetQinv(pVect1, pVect2);
+             kT12 = sqrt(pow(pVect1[1]+pVect2[1],2) + pow(pVect1[2]+pVect2[2],2))/2.;
+             SetFillBins2(ch1, ch2, bin1, bin2);
+             Int_t kTindex=0;
+             if(kT12<=0.3) kTindex=0;
+             else kTindex=1;
+             
+             FSICorr12 = FSICorrelation(ch1,ch2, qinv12);
+             
+             // two particle terms filled during tabulation of low-q pairs
+             
+             
+             if(fMCcase){
+               FilledMCpair12=kFALSE;
+
+               if(ch1==ch2 && fMbin==0 && qinv12<0.2 && ENsum!=3 && ENsum!=6){
+                 for(Int_t rstep=0; rstep<10; rstep++){
+                   Float_t coeff = (rstep)*0.2*(0.18/1.2);
+                   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;
+                   
+                   if(ENsum==0) ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiNum"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
+                   else ((TH3F*)fOutputList->FindObject("fPairsDetaDPhiDen"))->Fill(rstep, (fEvt)->fTracks[i].fEta-(fEvt+en2)->fTracks[j].fEta, deltaphi);
+                 }
+                 
+               }// pair selection
+
+               // Check that label does not exceed stack size
+               if((fEvt)->fTracks[i].fLabel < (fEvt)->fMCarraySize && (fEvt+en2)->fTracks[j].fLabel < (fEvt+en2)->fMCarraySize){
+                 if(ENsum==0 && abs((fEvt+en2)->fTracks[j].fLabel) == abs((fEvt)->fTracks[i].fLabel)) continue;
+                 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(pVect1MC, pVect2MC);
+                 
+                 if(qinv12<0.1 && ch1==ch2 && ENsum==0) {
+                   ((TProfile*)fOutputList->FindObject("fQsmearMean"))->Fill(1.,qinv12-qinv12MC); 
+                   ((TProfile*)fOutputList->FindObject("fQsmearSq"))->Fill(1.,1000.*pow(qinv12-qinv12MC,2));
+                   ((TH2D*)fOutputList->FindObject("fQ2Res"))->Fill(kT12, qinv12-qinv12MC);
+                 }
+                                 
+                 // secondary contamination
+                 if(ENsum==0){
+                   mcParticle1 = (AliAODMCParticle*)mcArray->At(abs((fEvt)->fTracks[i].fLabel));
+                   mcParticle2 = (AliAODMCParticle*)mcArray->At(abs((fEvt+en2)->fTracks[j].fLabel));
+                   if(!mcParticle1 || !mcParticle2) continue;
+                   if(abs(mcParticle1->GetPdgCode())==211 && abs(mcParticle2->GetPdgCode())==211){
+                     if(ch1==ch2) {
+                       ((TH3D*)fOutputList->FindObject("fAllSCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
+                       if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
+                         ((TH3D*)fOutputList->FindObject("fPrimarySCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
+                       }             
+                     }else{
+                       ((TH3D*)fOutputList->FindObject("fAllMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
+                       if(!mcParticle1->IsSecondaryFromWeakDecay() && !mcParticle2->IsSecondaryFromWeakDecay()) {
+                         ((TH3D*)fOutputList->FindObject("fPrimaryMCPionPairs"))->Fill(fMbin+1, kT12, qinv12);
+                       }
+                     }
+                   }
+                 }
+                 
+                 if(ENsum==6){// all mixed events
+                   Int_t chGroup2[2]={ch1,ch2};
+
+                   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].MB[fMbin].EDB[0].TwoPT[1].fMCqinv->Fill(qinv12MC, MCWeight(chGroup2, rForQW, 0.7, qinv12MC, 0.));// was 4,5
+                   Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fMCqinvQW->Fill(qinv12MC, qinv12MC*MCWeight(chGroup2, rForQW, 0.7, qinv12MC, 0.));// was 4,5
+                   // pion purity
+                   Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityDen->Fill(kT12, qinv12);
+                   Int_t SCNumber = 1;
+                   Int_t PdgCodeSum = abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode) + abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode);
+                   if(PdgCodeSum==22) SCNumber=1;// e-e
+                   else if(PdgCodeSum==24) SCNumber=2;// e-mu
+                   else if(PdgCodeSum==222) SCNumber=3;// e-pi
+                   else if(PdgCodeSum==332) SCNumber=4;// e-k
+                   else if(PdgCodeSum==2223) SCNumber=5;// e-p
+                   else if(PdgCodeSum==26) SCNumber=6;// mu-mu
+                   else if(PdgCodeSum==224) SCNumber=7;// mu-pi
+                   else if(PdgCodeSum==334) SCNumber=8;// mu-k
+                   else if(PdgCodeSum==2225) SCNumber=9;// mu-p
+                   else if(PdgCodeSum==422) SCNumber=10;// pi-pi
+                   else if(PdgCodeSum==532) SCNumber=11;// pi-k
+                   else if(PdgCodeSum==2423) SCNumber=12;// pi-p
+                   else if(PdgCodeSum==642) SCNumber=13;// k-k
+                   else if(PdgCodeSum==2533) SCNumber=14;// k-p
+                   else if(PdgCodeSum==4424) SCNumber=15;// p-p
+                   else {SCNumber=16;}
+                   
+                   Charge1[bin1].Charge2[bin2].MB[fMbin].EDB[0].TwoPT[1].fPIDpurityNum->Fill(SCNumber, kT12, qinv12);
+                   
+                   ///////////////////////
+                   // muon contamination
+                   Pparent1[0]=pVect1MC[0]; Pparent1[1]=pVect1MC[1]; Pparent1[2]=pVect1MC[2]; Pparent1[3]=pVect1MC[3];
+                   Pparent2[0]=pVect2MC[0]; Pparent2[1]=pVect2MC[1]; Pparent2[2]=pVect2MC[2]; Pparent2[3]=pVect2MC[3];
+                   pionParent1=kFALSE; pionParent2=kFALSE;
+                   FilledMCpair12=kTRUE;
+                   //
+                   if(abs((fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fPdgCode)==13 || abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13){// muon check
+                     Int_t MotherLabel1 = (fEvt)->fMCtracks[abs((fEvt)->fTracks[i].fLabel)].fMotherLabel;
+                     if(abs((fEvt)->fMCtracks[MotherLabel1].fPdgCode)==211) {
+                       pionParent1=kTRUE;
+                       Pparent1[1] = (fEvt)->fMCtracks[MotherLabel1].fPx; Pparent1[2] = (fEvt)->fMCtracks[MotherLabel1].fPy; Pparent1[3] = (fEvt)->fMCtracks[MotherLabel1].fPz;
+                       Pparent1[0] = sqrt(pow(Pparent1[1],2)+pow(Pparent1[2],2)+pow(Pparent1[3],2)+pow(fTrueMassPi,2));
+                     }
+                     // 
+                     if(abs((fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fPdgCode)==13) {
+                       Int_t MotherLabel2 = (fEvt+en2)->fMCtracks[abs((fEvt+en2)->fTracks[j].fLabel)].fMotherLabel;
+                       if(abs((fEvt+en2)->fMCtracks[MotherLabel2].fPdgCode)==211) {
+                         pionParent2=kTRUE;
+                         Pparent2[1] = (fEvt+en2)->fMCtracks[MotherLabel2].fPx; Pparent2[2] = (fEvt+en2)->fMCtracks[MotherLabel2].fPy; Pparent2[3] = (fEvt+en2)->fMCtracks[MotherLabel2].fPz;
+                         Pparent2[0] = sqrt(pow(Pparent2[1],2)+pow(Pparent2[2],2)+pow(Pparent2[3],2)+pow(fTrueMassPi,2));
+                       }
+                     }
+                     
+                     parentQinv12 = GetQinv(Pparent1, Pparent2);
+                     
+                     if(pionParent1 || pionParent2){
+                       if(parentQinv12 > 0.001 && parentQinv12 < 0.3){
+                         Float_t muonPionK12 = FSICorrelation(ch1, ch2, qinv12MC);
+                         Float_t pionPionK12 = FSICorrelation(ch1, ch2, parentQinv12);
+                         for(Int_t term=1; term<=2; term++){
+                           for(Int_t Riter=0; Riter<fRVALUES; Riter++){
+                             Float_t Rvalue = 5+Riter;
+                             Float_t WInput = 1.0;
+                             if(term==1) {
+                               WInput = MCWeight(chGroup2, Rvalue, 1.0, parentQinv12, 0.);
+                               muonPionK12 = 1.0; pionPionK12=1.0;
+                             }
+                             Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonSmeared->Fill(Rvalue, qinv12MC, WInput);
+                             Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonIdeal->Fill(Rvalue, parentQinv12, WInput);
+                             Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fMuonPionK2->Fill(Rvalue, qinv12MC, muonPionK12);
+                             Charge1[bin1].Charge2[bin2].MB[0].EDB[0].TwoPT[term-1].fPionPionK2->Fill(Rvalue, parentQinv12, pionPionK12);
+                           }// Riter
+                         }// term loop
+                         
+                         if(ch1==ch2) ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(0., kT12, qinv12MC-parentQinv12);
+                         else ((TH3D*)fOutputList->FindObject("fMuonPionDeltaQinv"))->Fill(1., kT12, qinv12MC-parentQinv12);
+                       }// parentQ check
+                     }// pion parent check
+                   }// muon check
+                 
+                   
+                   // momentum resolution
+                   for(Int_t Riter=0; Riter<fRVALUES; Riter++){
+                     Float_t Rvalue = 5+Riter;
+                     Float_t WInput = MCWeight(chGroup2, Rvalue, 0.7, qinv12MC, 0.);
+                     Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fIdeal->Fill(Rvalue, qinv12MC, WInput);
+                     Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fIdeal->Fill(Rvalue, qinv12MC);
+                     Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[0].fSmeared->Fill(Rvalue, qinv12, WInput);
+                     Charge1[bin1].Charge2[bin2].MB[0].EDB[kTindex].TwoPT[1].fSmeared->Fill(Rvalue, qinv12);
+                   }
+                   
+                 }// ENsum check
+               }// MC array check
+             }// MC case
+             
+               
+         
+             /////////////////////////////////////////////////////////////
+             for (Int_t k=j+1; k<(fEvt+en3)->fNtracks; k++) {// 3rd particle
+               if(en3==0) {
+                 if(!LowQPairSwitch_E0E0[i][k]) continue;
+                 if(!LowQPairSwitch_E0E0[j][k]) continue;
+               }else if(en3==1){
+                 if(!LowQPairSwitch_E0E1[i][k]) continue;
+                 if(!LowQPairSwitch_E0E1[j][k]) continue;
+               }else{
+                 if(!LowQPairSwitch_E0E2[i][k]) continue;
+                 if(!LowQPairSwitch_E1E2[j][k]) continue;
+               }
+               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];
+               ch3 = Int_t(((fEvt+en3)->fTracks[k].fCharge + 1)/2.);
+               qinv13 = GetQinv(pVect1, pVect3);
+               qinv23 = GetQinv(pVect2, pVect3);
+               q3 = sqrt(pow(qinv12,2) + pow(qinv13,2) + pow(qinv23,2));
+               
+               FilledMCtriplet123 = kFALSE;
+               
+               Bool_t fill2=kFALSE, fill3=kFALSE, fill4=kFALSE;
+               SetFillBins3(ch1, ch2, ch3, 1, bin1, bin2, bin3, fill2, fill3, fill4);
+               
+               Float_t KT3 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2],2))/3.;
+               if(KT3<=fKT3transition) KT3index=0;
+               else KT3index=1;
+               
+               FSICorr13 = FSICorrelation(ch1,ch3, qinv13);
+               FSICorr23 = FSICorrelation(ch2,ch3, qinv23);
+
+               if(ENsum==0) {
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fTerms3->Fill(q3); 
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fKfactor->Fill(q3, 1/(FSICorr12*FSICorr12*FSICorr12));
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv12);
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv13);
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[0].fMeanQinv->Fill(q3, qinv23);
+               }
+               if(ENsum==6) {
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTerms3->Fill(q3);
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv12);
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv13);
+                 Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fMeanQinv->Fill(q3, qinv23);
+               }
+               if(ENsum==3){
+                 if(fill2) {
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fTerms3->Fill(q3);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fKfactor->Fill(q3, 1/(FSICorr12));
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv12);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv13);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[1].fMeanQinv->Fill(q3, qinv23);
+                 }if(fill3) {
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fTerms3->Fill(q3);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fKfactor->Fill(q3, 1/(FSICorr12));
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv12);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv13);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[2].fMeanQinv->Fill(q3, qinv23);
+                 }if(fill4) {
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fTerms3->Fill(q3);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fKfactor->Fill(q3, 1/(FSICorr12));
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv12);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv13);
+                   Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[3].fMeanQinv->Fill(q3, qinv23);
+                 }
+               }
+               
+               // r3 denominator
+               if(ENsum==6 && ch1==ch2 && ch1==ch3){
+                 GoodTripletWeights = kFALSE;
+                 //
+                 GetWeight(pVect1, pVect2, weight12, weight12Err);
+                 GetWeight(pVect1, pVect3, weight13, weight13Err);
+                 GetWeight(pVect2, pVect3, weight23, weight23Err);
+                 
+                 if(sqrt(fabs(weight12*weight13*weight23)) > 1.0) {// weight should never be larger than 1
+                   if(fMbin==0 && bin1==0) {
+                     ((TH1D*)fOutputList->FindObject("fTPNRejects3pion1"))->Fill(q3, sqrt(fabs(weight12*weight13*weight23)));
+                   }
+                 }else{
+                   
+                   Float_t MomResCorr12=1.0, MomResCorr13=1.0, MomResCorr23=1.0;
+                   if(!fGenerateSignal && !fMCcase) {
+                     Int_t momBin12 = fMomResC2->GetYaxis()->FindBin(qinv12);
+                     Int_t momBin13 = fMomResC2->GetYaxis()->FindBin(qinv13);
+                     Int_t momBin23 = fMomResC2->GetYaxis()->FindBin(qinv23);            
+                     if(momBin12 >= 20) momBin12 = 19;
+                     if(momBin13 >= 20) momBin13 = 19;
+                     if(momBin23 >= 20) momBin23 = 19;
+                     MomResCorr12 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin12);
+                     MomResCorr13 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin13);
+                     MomResCorr23 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin23);
+                   }
+                  
+                   weight12CC = ((weight12+1)*MomResCorr12 - ffcSq*FSICorr12 - (1-ffcSq));
+                   weight12CC /= FSICorr12*ffcSq;
+                   weight13CC = ((weight13+1)*MomResCorr13 - ffcSq*FSICorr13 - (1-ffcSq));
+                   weight13CC /= FSICorr13*ffcSq;
+                   weight23CC = ((weight23+1)*MomResCorr23 - ffcSq*FSICorr23 - (1-ffcSq));
+                   weight23CC /= FSICorr23*ffcSq;
+
+                   if(weight12CC < 0 || weight13CC < 0 || weight23CC < 0) {// C2^QS can never be less than unity
+                     if(fMbin==0 && bin1==0) {
+                       ((TH1D*)fOutputList->FindObject("fTPNRejects3pion2"))->Fill(q3, sqrt(fabs(weight12CC*weight13CC*weight23CC)));
+                     }
+                   }else{
+                     GoodTripletWeights = kTRUE;
+                     /////////////////////////////////////////////////////
+                     weightTotal = sqrt(weight12CC*weight13CC*weight23CC);
+                     /////////////////////////////////////////////////////
+                     Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[fMbin].EDB[KT3index].ThreePT[4].fTwoPartNorm->Fill(1, q3, weightTotal);
+                   }// 2nd r3 den check                      
+                 }// 1st r3 den check
+
+               }// r3 den
+               
+
+               if(ch1==ch2 && ch1==ch3 && ENsum==0){
+                 ((TH3D*)fOutputList->FindObject("fKT3DistTerm1"))->Fill(fMbin+1, KT3, q3);
+                 if(q3<0.06){
+                   Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
+                   Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
+                   Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
+                   ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt1);
+                   ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt2);
+                   ((TProfile2D*)fOutputList->FindObject("fKT3AvgpT"))->Fill(fMbin+1, KT3index, pt3);
+                 }
+               }
+               if(ch1==ch2 && ch1==ch3 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT3DistTerm5"))->Fill(fMbin+1, KT3, q3);
+               
+               
+
+               
+               if(fMCcase && ENsum==6 && FilledMCpair12){// for momentum resolution and muon correction
+                 if((fEvt+en3)->fTracks[k].fLabel < (fEvt+en3)->fMCarraySize){
+                   
+                   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(pVect1MC, pVect3MC);
+                   qinv23MC = GetQinv(pVect2MC, pVect3MC);
+                   
+                   q3MC = sqrt(pow(qinv12MC,2)+pow(qinv13MC,2)+pow(qinv23MC,2));
+                   if(q3<0.1 && ch1==ch2 && ch1==ch3) ((TH2D*)fOutputList->FindObject("fQ3Res"))->Fill(KT3, q3-q3MC);
+                   
+                   Int_t chGroup3[3]={ch1,ch2,ch3};
+                   Float_t QinvMCGroup3[3]={qinv12MC,qinv13MC,qinv23MC};
+                   //Float_t kTGroup3[3]={float(sqrt(pow(pVect1MC[1]+pVect2MC[1],2) + pow(pVect1MC[2]+pVect2MC[2],2))/2.),
+                   //float(sqrt(pow(pVect1MC[1]+pVect3MC[1],2) + pow(pVect1MC[2]+pVect3MC[2],2))/2.),
+                   //float(sqrt(pow(pVect2MC[1]+pVect3MC[1],2) + pow(pVect2MC[2]+pVect3MC[2],2))/2.)};
+                   Float_t kTGroup3[3]={0};
+                   
+                   Pparent3[0]=pVect3MC[0]; Pparent3[1]=pVect3MC[1]; Pparent3[2]=pVect3MC[2]; Pparent3[3]=pVect3MC[3];
+                   pionParent3=kFALSE;
+                 
+                   if(abs((fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fPdgCode)==13){// muon check
+                     Int_t MotherLabel3 = (fEvt+en3)->fMCtracks[abs((fEvt+en3)->fTracks[k].fLabel)].fMotherLabel;
+                     if(abs((fEvt+en3)->fMCtracks[MotherLabel3].fPdgCode)==211) {
+                       pionParent3=kTRUE;
+                       Pparent3[1] = (fEvt+en3)->fMCtracks[MotherLabel3].fPx; Pparent3[2] = (fEvt+en3)->fMCtracks[MotherLabel3].fPy; Pparent3[3] = (fEvt+en3)->fMCtracks[MotherLabel3].fPz;
+                       Pparent3[0] = sqrt(pow(Pparent3[1],2)+pow(Pparent3[2],2)+pow(Pparent3[3],2)+pow(fTrueMassPi,2));
+                     }
+                   }
+                   
+                   parentQinv13 = GetQinv(Pparent1, Pparent3);
+                   parentQinv23 = GetQinv(Pparent2, Pparent3);
+                   parentQ3 = sqrt(pow(parentQinv12,2) + pow(parentQinv13,2) + pow(parentQinv23,2));
+                  
+                   if(parentQinv12 > 0.001 && parentQinv13 > 0.001 && parentQinv23 > 0.001 && parentQ3 < 0.5){
+                     FilledMCtriplet123=kTRUE;
+                     if(pionParent1 || pionParent2 || pionParent3) {// want at least one pion-->muon
+                       
+                       Float_t parentQinvGroup3[3]={parentQinv12, parentQinv13, parentQinv23};
+                       //Float_t parentkTGroup3[3]={float(sqrt(pow(Pparent1[1]+Pparent2[1],2) + pow(Pparent1[2]+Pparent2[2],2))/2.),
+                       //float(sqrt(pow(Pparent1[1]+Pparent3[1],2) + pow(Pparent1[2]+Pparent3[2],2))/2.),
+                       //float(sqrt(pow(Pparent2[1]+Pparent3[1],2) + pow(Pparent2[2]+Pparent3[2],2))/2.)};
+                       Float_t parentkTGroup3[3]={0};
+                       
+                       ((TH2D*)fOutputList->FindObject("fAvgQ12VersusQ3"))->Fill(parentQ3, parentQinv12);
+                       ((TH2D*)fOutputList->FindObject("fAvgQ13VersusQ3"))->Fill(parentQ3, parentQinv13);
+                       ((TH2D*)fOutputList->FindObject("fAvgQ23VersusQ3"))->Fill(parentQ3, parentQinv23);
+
+                       for(Int_t term=1; term<=5; term++){
+                         for(Int_t Riter=0; Riter<fRVALUES; Riter++){
+                           Float_t Rvalue = 5+Riter;
+                           Float_t WInput = MCWeight3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3, parentkTGroup3);
+                           Float_t WInputParentFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, parentQinvGroup3);
+                           Float_t WInputFSI = MCWeightFSI3(term, Rvalue, 1.0, chGroup3, QinvMCGroup3);
+                           Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonSmeared->Fill(Rvalue, q3MC, WInput);
+                           Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonIdeal->Fill(Rvalue, parentQ3, WInput);
+                           Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fMuonPionK3->Fill(Rvalue, q3MC, WInputFSI);
+                           Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[0].ThreePT[term-1].fPionPionK3->Fill(Rvalue, parentQ3, WInputParentFSI);
+                         }// Riter
+                       }// term loop
+                   
+                     }// pion parent check
+                   }// parentQ check (muon correction)
+                   
+                   // 3-pion momentum resolution
+                   for(Int_t term=1; term<=5; term++){
+                     for(Int_t Riter=0; Riter<fRVALUES; Riter++){
+                       Float_t Rvalue = 5+Riter;
+                       Float_t WInput = MCWeight3(term, Rvalue, 0.7, chGroup3, QinvMCGroup3, kTGroup3);
+                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fIdeal->Fill(Rvalue, q3MC, WInput);
+                       Charge1[bin1].Charge2[bin2].Charge3[bin3].MB[0].EDB[KT3index].ThreePT[term-1].fSmeared->Fill(Rvalue, q3, WInput);
+                     }
+                   }
+                   
+                 }// 3rd particle label check
+               }// MCcase and ENsum==6
+               
+
+               
+
+               /////////////////////////////////////////////////////////////
+               for (Int_t l=k+1; l<(fEvt+en4)->fNtracks; l++) {// 4th particle
+                 if(en4==0){
+                   if(!LowQPairSwitch_E0E0[i][l]) continue;
+                   if(!LowQPairSwitch_E0E0[j][l]) continue;
+                   if(!LowQPairSwitch_E0E0[k][l]) continue;
+                 }else if(en4==1){
+                   if(!LowQPairSwitch_E0E1[i][l]) continue;
+                   if(!LowQPairSwitch_E0E1[j][l]) continue;
+                   if(!LowQPairSwitch_E0E1[k][l]) continue;
+                 }else if(en4==2){
+                   if(!LowQPairSwitch_E0E2[i][l]) continue;
+                   if(!LowQPairSwitch_E0E2[j][l]) continue;
+                   if(!LowQPairSwitch_E1E2[k][l]) continue;
+                 }else{
+                   if(!LowQPairSwitch_E0E3[i][l]) continue;
+                   if(!LowQPairSwitch_E1E3[j][l]) continue;
+                   if(!LowQPairSwitch_E2E3[k][l]) continue;
+                 }
+                 
+                 pVect4[0]=(fEvt+en4)->fTracks[l].fEaccepted;
+                 pVect4[1]=(fEvt+en4)->fTracks[l].fP[0];
+                 pVect4[2]=(fEvt+en4)->fTracks[l].fP[1];
+                 pVect4[3]=(fEvt+en4)->fTracks[l].fP[2];
+                 ch4 = Int_t(((fEvt+en4)->fTracks[l].fCharge + 1)/2.);
+                 qinv14 = GetQinv(pVect1, pVect4);
+                 qinv24 = GetQinv(pVect2, pVect4);
+                 qinv34 = GetQinv(pVect3, pVect4);
+                 q4 = sqrt(pow(q3,2) + pow(qinv14,2) + pow(qinv24,2) + pow(qinv34,2));
+                 
+                 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6){
+                   ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(1, qinv12); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(2, qinv13); 
+                   ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(3, qinv14); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(4, qinv23); 
+                   ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(5, qinv24); ((TH2D*)fOutputList->FindObject("DistQinv4pion"))->Fill(6, qinv34);
+                 }
+                 
+                 Float_t KT4 = sqrt(pow(pVect1[1]+pVect2[1]+pVect3[1]+pVect4[1],2) + pow(pVect1[2]+pVect2[2]+pVect3[2]+pVect4[2],2))/4.;
+                 if(KT4<=fKT4transition) KT4index=0;
+                 else KT4index=1;
+                 
+                 FSICorr14 = FSICorrelation(ch1,ch4, qinv14);
+                 FSICorr24 = FSICorrelation(ch2,ch4, qinv24);
+                 FSICorr34 = FSICorrelation(ch3,ch4, qinv34);
+
+                 Bool_t FillTerms[12]={kFALSE};
+                 SetFillBins4(ch1, ch2, ch3, ch4, bin1, bin2, bin3, bin4, ENsum, FillTerms);
+                 //
+                 for(int ft=0; ft<12; ft++) {
+                   Float_t FSIfactor = 1.0;
+                   if(ft==0) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr14 * FSICorr23 * FSICorr24 * FSICorr34);
+                   else if(ft<=4) FSIfactor = 1/(FSICorr12 * FSICorr13 * FSICorr23);
+                   else if(ft<=10) FSIfactor = 1/(FSICorr12);
+                   else FSIfactor = 1.0;
+                   if(FillTerms[ft]) {
+                     Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fTerms4->Fill(q4);
+                     Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT4index].FourPT[ft].fKfactor->Fill(q4, FSIfactor);
+                   }
+                 }
+                 
+                 /////////////////////////////////////////////////////////////
+                 // r4{2}
+                 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6 && GoodTripletWeights){
+                   GetWeight(pVect1, pVect4, weight14, weight14Err);
+                   GetWeight(pVect2, pVect4, weight24, weight24Err);
+                   GetWeight(pVect3, pVect4, weight34, weight34Err);
+                   
+                   Float_t MomResCorr14=1.0, MomResCorr24=1.0, MomResCorr34=1.0;
+                   if(!fGenerateSignal && !fMCcase) {
+                     Int_t momBin14 = fMomResC2->GetYaxis()->FindBin(qinv14);
+                     Int_t momBin24 = fMomResC2->GetYaxis()->FindBin(qinv24);
+                     Int_t momBin34 = fMomResC2->GetYaxis()->FindBin(qinv34);            
+                     if(momBin14 >= 20) momBin14 = 19;
+                     if(momBin24 >= 20) momBin24 = 19;
+                     if(momBin34 >= 20) momBin34 = 19;
+                       MomResCorr14 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin14);
+                       MomResCorr24 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin24);
+                       MomResCorr34 = fMomResC2->GetBinContent(rBinForTPNMomRes, momBin34);
+                   }
+                   weight14CC = ((weight14+1)*MomResCorr14 - ffcSq*FSICorr14 - (1-ffcSq));
+                   weight14CC /= FSICorr14*ffcSq;
+                   weight24CC = ((weight24+1)*MomResCorr24 - ffcSq*FSICorr24 - (1-ffcSq));
+                   weight24CC /= FSICorr24*ffcSq;
+                   weight34CC = ((weight34+1)*MomResCorr34 - ffcSq*FSICorr34 - (1-ffcSq));
+                   weight34CC /= FSICorr34*ffcSq;
+       
+                   if(weight14CC < 0 || weight24CC < 0 || weight34CC < 0) {// C2^QS can never be less than unity
+                     if(fMbin==0 && bin1==0) {
+                       ((TH1D*)fOutputList->FindObject("fTPNRejects4pion1"))->Fill(q4, sqrt(fabs(weight12CC*weight23CC*weight34CC*weight14CC)));
+                     }
+                   }else{
+                     /////////////////////////////////////////////////////
+                     weightTotal  = sqrt(weight12CC*weight13CC*weight24CC*weight34CC) + sqrt(weight12CC*weight14CC*weight23CC*weight34CC);
+                     weightTotal += sqrt(weight13CC*weight14CC*weight23CC*weight24CC);
+                     weightTotal /= 3.;
+                     /////////////////////////////////////////////////////
+                     Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[fMbin].EDB[KT3index].FourPT[11].fTwoPartNorm->Fill(1, q4, weightTotal);
+                   }
+                 }
+                 /////////////////////////////////////////////////////////////
+
+                 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==0){
+                   ((TH3D*)fOutputList->FindObject("fKT4DistTerm1"))->Fill(fMbin+1, KT4, q4);
+                   if(q4<0.085){
+                     Float_t pt1=sqrt(pow(pVect1[1],2)+pow(pVect1[2],2));
+                     Float_t pt2=sqrt(pow(pVect2[1],2)+pow(pVect2[2],2));
+                     Float_t pt3=sqrt(pow(pVect3[1],2)+pow(pVect3[2],2));
+                     Float_t pt4=sqrt(pow(pVect4[1],2)+pow(pVect4[2],2));
+                     ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt1);
+                     ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt2);
+                     ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt3);
+                     ((TProfile2D*)fOutputList->FindObject("fKT4AvgpT"))->Fill(fMbin+1, KT4index, pt4);
+                   }
+                 }
+                 if(ch1==ch2 && ch1==ch3 && ch1==ch4 && ENsum==6) ((TH3D*)fOutputList->FindObject("fKT4DistTerm12"))->Fill(fMbin+1, KT4, q4);
+               
+
+                 // momenumtum resolution and muon corrections
+                 if(fMCcase && ENsum==6 && FilledMCtriplet123){// for momentum resolution and muon correction
+                   if((fEvt+en4)->fTracks[l].fLabel < (fEvt+en4)->fMCarraySize){
+                     
+                     pVect4MC[0]=sqrt(pow((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPtot,2)+pow(fTrueMassPi,2)); 
+                     pVect4MC[1]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPx;
+                     pVect4MC[2]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPy;
+                     pVect4MC[3]=(fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPz;
+                     qinv14MC = GetQinv(pVect1MC, pVect4MC);
+                     qinv24MC = GetQinv(pVect2MC, pVect4MC);
+                     qinv34MC = GetQinv(pVect3MC, pVect4MC);
+                     
+                     q4MC = sqrt(pow(q3MC,2) + pow(qinv14MC,2) +  pow(qinv24MC,2) +  pow(qinv34MC,2));
+                     if(q4<0.1 && ch1==ch2 && ch1==ch3 && ch1==ch4) ((TH2D*)fOutputList->FindObject("fQ4Res"))->Fill(KT4, q4-q4MC);
+                     if(ch1==ch2 && ch1==ch3 && ch1==ch4) {
+                       ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(1, qinv12MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(2, qinv13MC); 
+                       ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(3, qinv14MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(4, qinv23MC); 
+                       ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(5, qinv24MC); ((TH2D*)fOutputList->FindObject("DistQinvMC4pion"))->Fill(6, qinv34MC);
+                     }
+                     Int_t chGroup4[4]={ch1,ch2,ch3,ch4};
+                     Float_t QinvMCGroup4[6]={qinv12MC, qinv13MC, qinv14MC, qinv23MC, qinv24MC, qinv34MC};
+                     Float_t kTGroup4[6]={0};
+                     
+                     Pparent4[0]=pVect4MC[0]; Pparent4[1]=pVect4MC[1]; Pparent4[2]=pVect4MC[2]; Pparent4[3]=pVect4MC[3];
+                     pionParent4=kFALSE;
+                     if(abs((fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fPdgCode)==13){// muon check
+                       Int_t MotherLabel4 = (fEvt+en4)->fMCtracks[abs((fEvt+en4)->fTracks[l].fLabel)].fMotherLabel;
+                       if(abs((fEvt+en4)->fMCtracks[MotherLabel4].fPdgCode)==211) {
+                         pionParent4=kTRUE;
+                         Pparent4[1] = (fEvt+en4)->fMCtracks[MotherLabel4].fPx; Pparent4[2] = (fEvt+en4)->fMCtracks[MotherLabel4].fPy; Pparent4[3] = (fEvt+en4)->fMCtracks[MotherLabel4].fPz;
+                         Pparent4[0] = sqrt(pow(Pparent4[1],2)+pow(Pparent4[2],2)+pow(Pparent4[3],2)+pow(fTrueMassPi,2));
+                       }
+                     }
+
+                     parentQinv14 = GetQinv(Pparent1, Pparent4);
+                     parentQinv24 = GetQinv(Pparent2, Pparent4);
+                     parentQinv34 = GetQinv(Pparent3, Pparent4);
+                     Float_t parentQ4 = sqrt(pow(parentQ3,2) + pow(parentQinv14,2) + pow(parentQinv24,2) + pow(parentQinv34,2));
+                     
+                     if(parentQinv14 > 0.001 && parentQinv24 > 0.001 && parentQinv34 > 0.001 && parentQ4 < 0.5){
+                       if(pionParent1 || pionParent2 || pionParent3 || pionParent4) {// want at least one pion-->muon
+                        
+                         Float_t parentQinvGroup4[6]={parentQinv12, parentQinv13, parentQinv14, parentQinv23, parentQinv24, parentQinv34};
+                         Float_t parentkTGroup4[6]={0};
+                         
+                         for(Int_t term=1; term<=12; term++){
+                           for(Int_t Riter=0; Riter<fRVALUES; Riter++){
+                             Float_t Rvalue = 5+Riter;
+                             Float_t WInput = MCWeight4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4, parentkTGroup4);
+                             Float_t WInputParentFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, parentQinvGroup4);
+                             Float_t WInputFSI = MCWeightFSI4(term, Rvalue, 1.0, chGroup4, QinvMCGroup4);
+                             Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonSmeared->Fill(Rvalue, q4MC, WInput);
+                             Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonIdeal->Fill(Rvalue, parentQ4, WInput);
+                             Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fMuonPionK4->Fill(Rvalue, q4MC, WInputFSI);
+                             Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[0].FourPT[term-1].fPionPionK4->Fill(Rvalue, parentQ4, WInputParentFSI);
+                           }// Riter
+                         }// term loop
+                         
+                       }// pion parent check
+                     }// parentQ check (muon correction)
+                     
+                     // 4-pion momentum resolution
+                     for(Int_t term=1; term<=12; term++){
+                       for(Int_t Riter=0; Riter<fRVALUES; Riter++){
+                         Float_t Rvalue = 5+Riter;
+                         Float_t WInput = MCWeight4(term, Rvalue, 0.7, chGroup4, QinvMCGroup4, kTGroup4);
+                         Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fIdeal->Fill(Rvalue, q4MC, WInput);
+                         Charge1[bin1].Charge2[bin2].Charge3[bin3].Charge4[bin4].MB[0].EDB[KT4index].FourPT[term-1].fSmeared->Fill(Rvalue, q4, WInput);
+                       }
+                     }
+                     
+                   }// label check particle 4
+                 }// MCcase
+                 
+               }// 4th particle
+             }// 3rd particle
+           }// 2nd particle
+         }// 1st particle
+         
+       }// en4
+      }// en3
+    }// en2
+    
+    
+
+  
+  
+  
+  
+  // Post output data.
+  PostData(1, fOutputList);
+  
+}
+//________________________________________________________________________
+void AliFourPion::Terminate(Option_t *) 
+{
+  // Called once at the end of the query
+  cout<<"Done"<<endl;
+
+}
+//________________________________________________________________________
+Bool_t AliFourPion::AcceptPair(AliFourPionTrackStruct first, AliFourPionTrackStruct 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 AliFourPion::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 AliFourPion::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;
+  }
+}
+
+
+//________________________________________________________________________
+Float_t AliFourPion::GetQinv(Float_t track1[], Float_t track2[]){
+  
+  Float_t qinv = sqrt( fabs(pow(track1[1]-track2[1],2) + pow(track1[2]-track2[2],2) + pow(track1[3]-track2[3],2) - pow(track1[0]-track2[0],2)) );
+  return qinv;
+  
+}
+//________________________________________________________________________
+void AliFourPion::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;
+}
+//________________________________________________________________________
+void AliFourPion::SetWeightArrays(Bool_t legoCase, TH3F *histos[AliFourPion::fKbinsT][AliFourPion::fCentBins]){
+
+  if(legoCase){
+    cout<<"LEGO call to SetWeightArrays"<<endl;
+    
+    for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
+      for(Int_t mb=0; mb<fCentBins; mb++){
+       fNormWeight[tKbin][mb] = (TH3F*)histos[tKbin][mb]->Clone();
+       fNormWeight[tKbin][mb]->SetDirectory(0);
+      }
+    }
+    
+  }else{
+    
+    TFile *wFile = new TFile("WeightFile.root","READ");
+    if(!wFile->IsOpen()) {cout<<"No Weight File!!!!!!!!!!"<<endl; return;}
+    else cout<<"Good Weight File Found!"<<endl;
+    
+    for(Int_t tKbin=0; tKbin<fKbinsT; tKbin++){
+      for(Int_t mb=0; mb<fCentBins; mb++){
+                   
+       TString *name = new TString("Weight_Kt_");
+       *name += tKbin;
+       name->Append("_Ky_0");
+       name->Append("_M_");
+       *name += mb;
+       name->Append("_ED_0");
+       
+       
+       fNormWeight[tKbin][mb] = (TH3F*)wFile->Get(name->Data());
+       fNormWeight[tKbin][mb]->SetDirectory(0);
+       
+       
+      }//mb
+    }//kt
+    
+    wFile->Close();
+  }
+  
+  cout<<"Done reading weight file"<<endl;
+  
+}
+//________________________________________________________________________
+void AliFourPion::GetWeight(Float_t track1[], Float_t track2[], Float_t& wgt, Float_t& wgtErr){
+  
+  Float_t kt=sqrt( pow(track1[1]+track2[1],2) + pow(track1[2]+track2[2],2))/2.;
+  //
+  Float_t qOut=0,qSide=0,qLong=0;
+  GetQosl(track1, track2, qOut, qSide, qLong);
+  qOut = fabs(qOut);
+  qSide = fabs(qSide);
+  qLong = fabs(qLong);
+  Float_t wd=0, xd=0, yd=0, zd=0;
+  //Float_t qinvtemp=GetQinv(0,track1, track2);
+  //
+  
+  if(kt < fKmeanT[0]) {fKtIndexL=0; fKtIndexH=1;}
+  else if(kt >= fKmeanT[fKbinsT-1]) {fKtIndexL=fKbinsT-2; fKtIndexH=fKbinsT-1;}
+  else {
+    for(Int_t i=0; i<fKbinsT-1; i++){
+      if((kt >= fKmeanT[i]) && (kt < fKmeanT[i+1])) {fKtIndexL=i; fKtIndexH=i+1; break;}
+    }
+  }
+  wd = (kt-fKmeanT[fKtIndexL])/(fKmeanT[fKtIndexH]-fKmeanT[fKtIndexL]);
+  //
+  /////////
+  if(qOut < fQmean[0]) {fQoIndexL=0; fQoIndexH=0; xd=0;}
+  else if(qOut >= fQmean[kQbinsWeights-1]) {fQoIndexL=kQbinsWeights-1; fQoIndexH=kQbinsWeights-1; xd=1;}
+  else {
+    for(Int_t i=0; i<kQbinsWeights-1; i++){
+      if((qOut >= fQmean[i]) && (qOut < fQmean[i+1])) {fQoIndexL=i; fQoIndexH=i+1; break;}
+    }
+    xd = (qOut-fQmean[fQoIndexL])/(fQmean[fQoIndexH]-fQmean[fQoIndexL]);
+  }
+  //
+  if(qSide < fQmean[0]) {fQsIndexL=0; fQsIndexH=0; yd=0;}
+  else if(qSide >= fQmean[kQbinsWeights-1]) {fQsIndexL=kQbinsWeights-1; fQsIndexH=kQbinsWeights-1; yd=1;}
+  else {
+    for(Int_t i=0; i<kQbinsWeights-1; i++){
+      if((qSide >= fQmean[i]) && (qSide < fQmean[i+1])) {fQsIndexL=i; fQsIndexH=i+1; break;}
+    }
+    yd = (qSide-fQmean[fQsIndexL])/(fQmean[fQsIndexH]-fQmean[fQsIndexL]);
+  }
+  //
+  if(qLong < fQmean[0]) {fQlIndexL=0; fQlIndexH=0; zd=0;}
+  else if(qLong >= fQmean[kQbinsWeights-1]) {fQlIndexL=kQbinsWeights-1; fQlIndexH=kQbinsWeights-1; zd=1;}
+  else {
+    for(Int_t i=0; i<kQbinsWeights-1; i++){
+      if((qLong >= fQmean[i]) && (qLong < fQmean[i+1])) {fQlIndexL=i; fQlIndexH=i+1; break;}
+    }
+    zd = (qLong-fQmean[fQlIndexL])/(fQmean[fQlIndexH]-fQmean[fQlIndexL]);
+  }
+  //
+
+  
+  //Float_t min = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
+  Float_t minErr = fNormWeight[fKtIndexL][fMbin]->GetBinError(fQoIndexH+1,fQsIndexH+1,fQlIndexH+1);
+  /*
+  Float_t deltaW=0;
+  // kt
+  deltaW += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - min)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
+  // Qo 
+  deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - min)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
+  // Qs
+  deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - min)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
+  // Ql
+  deltaW += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - min)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
+  //
+  wgt = min + deltaW;
+  */
+  
+  //
+  // w interpolation (kt)
+  Float_t c000 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexL+1)*wd;
+  Float_t c100 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexL+1)*wd;
+  Float_t c010 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexL+1)*wd;
+  Float_t c001 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexL+1, fQlIndexH+1)*wd;
+  Float_t c110 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1)*wd;
+  Float_t c101 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1)*wd;
+  Float_t c011 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1)*wd;
+  Float_t c111 = fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*(1-wd) + fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1)*wd;
+  // x interpolation (qOut)
+  Float_t c00 = c000*(1-xd) + c100*xd;
+  Float_t c10 = c010*(1-xd) + c110*xd;
+  Float_t c01 = c001*(1-xd) + c101*xd;
+  Float_t c11 = c011*(1-xd) + c111*xd;
+  // y interpolation (qSide)
+  Float_t c0 = c00*(1-yd) + c10*yd;
+  Float_t c1 = c01*(1-yd) + c11*yd;
+  // z interpolation (qLong)
+  wgt = (c0*(1-zd) + c1*zd);
+  
+  ////
+  
+  // Denominator errors negligible compared to numerator so do not waste cpu time below.  
+  //Float_t deltaWErr=0;
+  // Kt
+  /*
+  deltaWErr += (fNormWeight[fKtIndexH][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(kt-fKmeanT[fKtIndexL])/((fKstepT[fKtIndexL]+fKstepT[fKtIndexH])/2.);
+  // Qo 
+  deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexL+1, fQsIndexH+1, fQlIndexH+1) - minErr)*(qOut-fQmean[fQoIndexL])/fQstepWeights;
+  // Qs
+  deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexL+1, fQlIndexH+1) - minErr)*(qSide-fQmean[fQsIndexL])/fQstepWeights;
+  // Ql
+  deltaWErr += (fNormWeight[fKtIndexL][fMbin]->GetBinContent(fQoIndexH+1, fQsIndexH+1, fQlIndexL+1) - minErr)*(qLong-fQmean[fQlIndexL])/fQstepWeights;
+  */
+  wgtErr = minErr;
+  
+}
+//________________________________________________________________________
+Float_t AliFourPion::MCWeight(Int_t c[2], Float_t R, Float_t fcSq, Float_t qinv, Float_t k12){
+  
+  Float_t radius = R/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
+  Float_t r12=radius*(1-k12/2.0);
+  SetFSIindex(R);
+  Float_t coulCorr12 = FSICorrelation(c[0], c[1], qinv);
+  if(c[0]==c[1]){
+    Float_t arg=qinv*r12;
+    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-fcSq) + fcSq*(1 + exp(-pow(qinv*r12,2))*pow(EW,2))*coulCorr12);
+  }else {
+    return ((1-fcSq) + fcSq*coulCorr12);
+  }
+    
+}
+//________________________________________________________________________
+Float_t AliFourPion::MCWeightOSL(Int_t charge1, Int_t charge2, Int_t r, Int_t dampIndex, Float_t qinv, Float_t qo, Float_t qs, Float_t ql){
+  
+  Float_t radiusOut = Float_t(r)/0.19733;// convert to GeV (starts at 5 fm, was 3 fm)
+  Float_t radiusSide = radiusOut;
+  Float_t radiusLong = radiusOut;
+  Float_t myDamp = fDampStart + (fDampStep)*dampIndex;
+  Float_t coulCorr12 = FSICorrelation(charge1, charge2, qinv);
+  if(charge1==charge2){
+    return ((1-myDamp) + myDamp*(1 + exp(-pow(qo*radiusOut,2)) * exp(-pow(qs*radiusSide,2)) * exp(-pow(ql*radiusLong,2)))*coulCorr12);
+  }else {
+    return ((1-myDamp) + myDamp*coulCorr12);
+  }
+    
+}
+
+//________________________________________________________________________
+Float_t AliFourPion::MCWeight3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3], Float_t kT[3]){
+  // FSI + QS correlations
+  if(term==5) return 1.0;
+  
+  Float_t radius=R/0.19733;
+  Float_t r12=radius*(1-kT[0]/2.0);
+  Float_t r13=radius*(1-kT[1]/2.0);
+  Float_t r23=radius*(1-kT[2]/2.0);
+  Float_t fc = sqrt(fcSq);
+  
+  SetFSIindex(R);
+  Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
+  Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
+  Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
+  
+  if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
+    Float_t arg12=qinv[0]*r12;
+    Float_t arg13=qinv[1]*r13;
+    Float_t arg23=qinv[2]*r23;
+    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(qinv[0]*r12,2))*pow(EW12,2) + exp(-pow(qinv[1]*r13,2))*pow(EW13,2) + exp(-pow(qinv[2]*r23,2))*pow(EW23,2);
+      C3QS += 2*exp(-(pow(r12,2)*pow(qinv[0],2) + pow(r13,2)*pow(qinv[1],2) + pow(r23,2)*pow(qinv[2],2))/2.)*EW12*EW13*EW23;
+      Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12;
+      C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13;
+      C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23;
+      C3 += pow(fc,3)*C3QS*Kfactor12*Kfactor13*Kfactor23;
+      return C3;
+    }else if(term==2){
+      return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[0]*r12,2))*pow(EW12,2))*Kfactor12);
+    }else if(term==3){
+      return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[1]*r13,2))*pow(EW13,2))*Kfactor13);
+    }else if(term==4){
+      return ((1-fcSq) + fcSq*(1 + exp(-pow(qinv[2]*r23,2))*pow(EW23,2))*Kfactor23);
+    }else return 1.0;
+    
+  }else{// mixed charge case
+    Float_t arg=qinv[0]*r12;
+    Float_t KfactorSC = Kfactor12;
+    Float_t KfactorMC1 = Kfactor13;
+    Float_t KfactorMC2 = Kfactor23;
+    if(c[0]==c[2]) {arg=qinv[1]*r13; KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;} 
+    if(c[1]==c[2]) {arg=qinv[2]*r23; KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
+    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);
+    if(term==1){
+      Float_t C3QS = 1 + exp(-pow(arg,2))*pow(EW,2);
+      Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      C3 += pow(fc,2)*(1-fc)*(1+exp(-pow(arg,2))*pow(EW,2))*KfactorSC;
+      C3 += pow(fc,2)*(1-fc)*KfactorMC1;
+      C3 += pow(fc,2)*(1-fc)*KfactorMC2;
+      C3 += pow(fc,3)*C3QS*KfactorSC*KfactorMC1*KfactorMC2;
+      return C3;
+    }else if(term==2){
+      if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
+      else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
+    }else if(term==3){
+      return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
+    }else if(term==4){
+      if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*(1 + exp(-pow(arg,2))*pow(EW,2))*KfactorSC);
+      else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
+    }else return 1.0;
+  }
+  
+}
+//________________________________________________________________________
+Float_t AliFourPion::MCWeightFSI3(Int_t term, Float_t R, Float_t fcSq, Int_t c[3], Float_t qinv[3]){
+  // FSI only (no QS correlations)
+  if(term==5) return 1.0;
+  
+  Float_t fc = sqrt(fcSq);
+  SetFSIindex(R);
+  Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
+  Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
+  Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[2]);// K2
+  
+  if(c[0]==c[1] && c[0]==c[2]){// all three of the same charge
+    if(term==1){
+      Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      C3 += pow(fc,2)*(1-fc)*Kfactor12;
+      C3 += pow(fc,2)*(1-fc)*Kfactor13;
+      C3 += pow(fc,2)*(1-fc)*Kfactor23;
+      C3 += pow(fc,3)*Kfactor12*Kfactor13*Kfactor23;
+      return C3;
+    }else if(term==2){
+      return ((1-fcSq) + fcSq*Kfactor12);
+    }else if(term==3){
+      return ((1-fcSq) + fcSq*Kfactor13);
+    }else if(term==4){
+      return ((1-fcSq) + fcSq*Kfactor23);
+    }else return 1.0;
+    
+  }else{// mixed charge case
+    Float_t KfactorSC = Kfactor12;
+    Float_t KfactorMC1 = Kfactor13;
+    Float_t KfactorMC2 = Kfactor23;
+    if(c[0]==c[2]) {KfactorSC = Kfactor13; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor23;} 
+    if(c[1]==c[2]) {KfactorSC = Kfactor23; KfactorMC1 = Kfactor12; KfactorMC2 = Kfactor13;}
+    if(term==1){
+      Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      C3 += pow(fc,2)*(1-fc)*KfactorSC;
+      C3 += pow(fc,2)*(1-fc)*KfactorMC1;
+      C3 += pow(fc,2)*(1-fc)*KfactorMC2;
+      C3 += pow(fc,3)*KfactorSC*KfactorMC1*KfactorMC2;
+      return C3;
+    }else if(term==2){
+      if( (c[0]+c[1]+c[2])==1) return ((1-fcSq) + fcSq*KfactorSC);
+      else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
+    }else if(term==3){
+      return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
+    }else if(term==4){
+      if( (c[0]+c[1]+c[2])==2) return ((1-fcSq) + fcSq*KfactorSC);
+      else return ((1-fcSq) + fcSq*KfactorMC1);// doesn't matter MC1 or MC2
+    }else return 1.0;
+  }
+  
+}
+//________________________________________________________________________
+Float_t AliFourPion::MCWeight4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6], Float_t kT[6]){
+  if(term==12) return 1.0;
+
+  // Charge ordering:
+  // ----, ---+, --++, -+++, ++++
+  //
+  // term ordering:
+  // Term 1: 1-2-3-4  (all same-event)
+  // Term 2: 1-2-3 4 (particle 4 from different event)
+  // Term 3: 1-2-4 3 (particle 3 from different event)
+  // Term 4: 1-3-4 2 (particle 2 from different event)
+  // Term 5: 2-3-4 1 (particle 1 from different event)
+  // Term 6: 1-2 3 4 (particle 1 and 2 from same event)
+  // Term 7: 1-3 2 4
+  // Term 8: 1-4 2 3
+  // Term 9: 2-3 1 4
+  // Term 10: 2-4 1 3
+  // Term 11: 3-4 1 2
+  // Term 12: 1 2 3 4 (all from different events)
+
+  Float_t radius = R/0.19733;
+  Float_t r[6]={0};
+  r[0]=radius*(1-kT[0]/2.0);
+  r[1]=radius*(1-kT[1]/2.0);
+  r[2]=radius*(1-kT[2]/2.0);
+  r[3]=radius*(1-kT[3]/2.0);
+  r[4]=radius*(1-kT[4]/2.0);
+  r[5]=radius*(1-kT[5]/2.0);
+    
+  Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
+  Float_t fc = sqrt(fcSq);
+  SetFSIindex(R);
+  Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
+  Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
+  Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
+  Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
+  Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
+  Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
+  Float_t arg12=qinv[0]*r[0];
+  Float_t arg13=qinv[1]*r[1];
+  Float_t arg14=qinv[2]*r[2];
+  Float_t arg23=qinv[3]*r[3];
+  Float_t arg24=qinv[4]*r[4];
+  Float_t arg34=qinv[5]*r[5];
+  // Exchange Amplitudes
+  Float_t EA12 = exp(-pow(arg12,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg12,3) - 12.*arg12) + kappa4/(24.*pow(2.,2))*(16.*pow(arg12,4) -48.*pow(arg12,2) + 12));
+  Float_t EA13 = exp(-pow(arg13,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg13,3) - 12.*arg13) + kappa4/(24.*pow(2.,2))*(16.*pow(arg13,4) -48.*pow(arg13,2) + 12));
+  Float_t EA14 = exp(-pow(arg14,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg14,3) - 12.*arg14) + kappa4/(24.*pow(2.,2))*(16.*pow(arg14,4) -48.*pow(arg14,2) + 12));
+  Float_t EA23 = exp(-pow(arg23,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg23,3) - 12.*arg23) + kappa4/(24.*pow(2.,2))*(16.*pow(arg23,4) -48.*pow(arg23,2) + 12));
+  Float_t EA24 = exp(-pow(arg24,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg24,3) - 12.*arg24) + kappa4/(24.*pow(2.,2))*(16.*pow(arg24,4) -48.*pow(arg24,2) + 12));
+  Float_t EA34 = exp(-pow(arg34,2)/2.)*(1 + kappa3/(6.*pow(2.,1.5))*(8.*pow(arg34,3) - 12.*arg34) + kappa4/(24.*pow(2.,2))*(16.*pow(arg34,4) -48.*pow(arg34,2) + 12));
+  
+  if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
+    
+    if(term==1){
+      Float_t C4QS = 1 + pow(EA12,2) + pow(EA13,2) + pow(EA14,2) + pow(EA23,2) + pow(EA24,2) + pow(EA34,2);// baseline + single pairs
+      C4QS += pow(EA12,2) * pow(EA34,2);// 2-pairs
+      C4QS += pow(EA13,2) * pow(EA24,2);// 2-pairs
+      C4QS += pow(EA14,2) * pow(EA23,2);// 2-pairs
+      C4QS += 2*EA12*EA13*EA23 + 2*EA12*EA14*EA24 + 2*EA13*EA14*EA34 + 2*EA23*EA24*EA34;// 3-particle exhange
+      C4QS += 3*EA12*EA23*EA34*EA14 + 3*EA12*EA13*EA34*EA24;// 4-particle exchange
+      Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
+      C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA12,2))*Kfactor12 + (1 + pow(EA13,2))*Kfactor13 + (1 + pow(EA14,2))*Kfactor14 );
+      C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA23,2))*Kfactor23 + (1 + pow(EA24,2))*Kfactor24 + (1 + pow(EA34,2))*Kfactor34);
+      C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA13,2) + pow(EA23,2) + 2*EA12*EA13*EA23) * Kfactor12*Kfactor13*Kfactor23;
+      C4 += pow(fc,3)*(1-fc)*(1 + pow(EA12,2) + pow(EA14,2) + pow(EA24,2) + 2*EA12*EA14*EA24) * Kfactor12*Kfactor14*Kfactor24;
+      C4 += pow(fc,3)*(1-fc)*(1 + pow(EA13,2) + pow(EA14,2) + pow(EA34,2) + 2*EA13*EA14*EA34) * Kfactor13*Kfactor14*Kfactor34;
+      C4 += pow(fc,3)*(1-fc)*(1 + pow(EA23,2) + pow(EA24,2) + pow(EA34,2) + 2*EA23*EA24*EA34) * Kfactor23*Kfactor24*Kfactor34;
+      C4 += pow(fc,4)*C4QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
+      return C4;
+    }else if(term<=5){
+      Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0;
+      if(term==2)      {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
+      else if(term==3) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
+      else if(term==4) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
+      else             {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
+      Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2);
+      C3QS += 2*EA1*EA2*EA3;
+      Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      C3 += pow(fc,2)*(1-fc)*( (1+pow(EA1,2))*Kpair1 + (1+pow(EA2,2))*Kpair2 + (1+pow(EA3,2))*Kpair3 );
+      C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
+      return C3;
+    }else if(term<=11){
+      if(term==6)       return ((1-fcSq) + fcSq*(1 + pow(EA12,2))*Kfactor12);
+      else if(term==7)  return ((1-fcSq) + fcSq*(1 + pow(EA13,2))*Kfactor13);
+      else if(term==8)  return ((1-fcSq) + fcSq*(1 + pow(EA14,2))*Kfactor14);
+      else if(term==9)  return ((1-fcSq) + fcSq*(1 + pow(EA23,2))*Kfactor23);
+      else if(term==10) return ((1-fcSq) + fcSq*(1 + pow(EA24,2))*Kfactor24);
+      else              return ((1-fcSq) + fcSq*(1 + pow(EA34,2))*Kfactor34);
+    }else return 1.0;
+    
+  }else{// mixed charge case
+    if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
+      Float_t EA1=0, EA2=0, EA3=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
+      Int_t c_OddOneOut = 1;
+      if(ChargeSum==3) c_OddOneOut = 0;
+      //
+      if(c[0]==c_OddOneOut) {EA1=EA23; EA2=EA24; EA3=EA34; Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
+      else if(c[1]==c_OddOneOut) {EA1=EA13; EA2=EA14; EA3=EA34; Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
+      else if(c[2]==c_OddOneOut) {EA1=EA12; EA2=EA14; EA3=EA24; Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;   Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
+      else {EA1=EA12; EA2=EA13; EA3=EA23; Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;   Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
+      
+      if(term==1){
+       Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
+       Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
+       C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 + (1 + pow(EA3,2))*Kpair3 );
+       C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
+       C4 += pow(fc,3)*(1-fc)*(1 + pow(EA1,2) + pow(EA2,2) +  pow(EA3,2) + 2*EA1*EA2*EA3) * Kpair1*Kpair2*Kpair3;
+       C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair4*Kpair5 + (1+pow(EA2,2))*Kpair2*Kpair4*Kpair6 + (1+pow(EA3,2))*Kpair3*Kpair5*Kpair6);// doesn't matter which MC K's
+       C4 += pow(fc,4)*C3QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
+       return C4;
+      }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
+       Float_t C3QS = 1 + pow(EA1,2) + pow(EA2,2) + pow(EA3,2) + 2*EA1*EA2*EA3;
+       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+       C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;
+       C3 += pow(fc,2)*(1-fc)*(1+pow(EA2,2))*Kpair2;
+       C3 += pow(fc,2)*(1-fc)*(1+pow(EA3,2))*Kpair3;
+       C3 += pow(fc,3)*C3QS*Kpair1*Kpair2*Kpair3;
+       return C3;
+      }else if(term<=5){// one SC pair, two MC pairs
+       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+       C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
+       C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair5;
+       return C3;
+      }else if(term==6 || term==7){
+       if(ChargeSum==1) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
+       else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
+      }else if(term==8){
+       return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
+      }else if(term==9){
+       return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
+      }else if(term==10 || term==1){
+       if(ChargeSum==3) return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
+       else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
+      }else return 1.0;
+    }else{// --++ configuration
+      Float_t EA1=0, EA2=0, Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
+      if(c[0]==c[1]) {EA1=EA12; EA2=EA34; Kpair1=Kfactor12; Kpair2=Kfactor34;    Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
+      else if(c[0]==c[2]) {EA1=EA13; EA2=EA24; Kpair1=Kfactor13; Kpair2=Kfactor24;    Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
+      else {EA1=EA14; EA2=EA23; Kpair1=Kfactor14; Kpair2=Kfactor23;    Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
+      //
+      if(term==1){
+       Float_t C2QS = 1 + pow(EA1,2)*pow(EA2,2);
+       Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
+       C4 += pow(fc,2)*pow(1-fc,2)*( (1 + pow(EA1,2))*Kpair1 + (1 + pow(EA2,2))*Kpair2 );
+       C4 += pow(fc,2)*pow(1-fc,2)*( Kpair3 + Kpair4 + Kpair5 + Kpair6 );
+       C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair3*Kpair4 + (1 + pow(EA2,2))*Kpair2*Kpair3*Kpair4);
+       C4 += pow(fc,3)*(1-fc)*( (1 + pow(EA1,2))*Kpair1*Kpair5*Kpair6 + (1 + pow(EA2,2))*Kpair2*Kpair5*Kpair6);// doesn't matter which two MC K's used
+       C4 += pow(fc,4)*C2QS*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
+       return C4;
+      }else if(term<=5){
+       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+       C3 += pow(fc,2)*(1-fc)*(1+pow(EA1,2))*Kpair1;// any SC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
+       C3 += pow(fc,3)*(1+pow(EA1,2))*Kpair1*Kpair4*Kpair6;
+       return C3;
+      }else if(term==6 || term==11){
+       return ((1-fcSq) + fcSq*(1 + pow(EA1,2))*Kpair1);// any SC pair will do
+      }else if(term !=12){
+       return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
+      }else return 1.0;
+    }
+  }
+  
+}
+//________________________________________________________________________
+Float_t AliFourPion::MCWeightFSI4(Int_t term, Float_t R, Float_t fcSq, Int_t c[4], Float_t qinv[6]){
+  if(term==12) return 1.0;
+    
+  Int_t ChargeSum=c[0]+c[1]+c[2]+c[3];
+  Float_t fc = sqrt(fcSq);
+  SetFSIindex(R);
+  Float_t Kfactor12 = FSICorrelation(c[0],c[1], qinv[0]);// K2
+  Float_t Kfactor13 = FSICorrelation(c[0],c[2], qinv[1]);// K2
+  Float_t Kfactor14 = FSICorrelation(c[0],c[3], qinv[2]);// K2
+  Float_t Kfactor23 = FSICorrelation(c[1],c[2], qinv[3]);// K2
+  Float_t Kfactor24 = FSICorrelation(c[1],c[3], qinv[4]);// K2
+  Float_t Kfactor34 = FSICorrelation(c[2],c[3], qinv[5]);// K2
+  
+  if(c[0]==c[1] && c[0]==c[2] && c[0]==c[3]){// ---- and ++++ configuration
+    
+    if(term==1){
+      Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
+      C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor12 + Kfactor13 + Kfactor14 );
+      C4 += pow(fc,2)*pow(1-fc,2)*( Kfactor23 + Kfactor24 + Kfactor34 );
+      C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor13*Kfactor23;
+      C4 += pow(fc,3)*(1-fc) * Kfactor12*Kfactor14*Kfactor24;
+      C4 += pow(fc,3)*(1-fc) * Kfactor13*Kfactor14*Kfactor34;
+      C4 += pow(fc,3)*(1-fc) * Kfactor23*Kfactor24*Kfactor34;
+      C4 += pow(fc,4) * Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
+      return C4;
+    }else if(term<=5){
+      Float_t Kpair1=0, Kpair2=0, Kpair3=0;
+      if(term==2) {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;}
+      else if(term==3) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;}
+      else if(term==4) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;}
+      else {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;}
+      Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+      C3 += pow(fc,2)*(1-fc)*( Kpair1 + Kpair2 + Kpair3 );
+      C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
+      return C3;
+    }else if(term<=11){
+      if(term==6) return ((1-fcSq) + fcSq*Kfactor12);
+      else if(term==7) return ((1-fcSq) + fcSq*Kfactor13);
+      else if(term==8) return ((1-fcSq) + fcSq*Kfactor14);
+      else if(term==9) return ((1-fcSq) + fcSq*Kfactor23);
+      else if(term==10) return ((1-fcSq) + fcSq*Kfactor24);
+      else return ((1-fcSq) + fcSq*Kfactor34);
+    }else return 1.0;
+    
+  }else{// mixed charge case
+    if( ChargeSum==1 || ChargeSum==3){// ---+ and -+++ configuration
+      Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
+      Int_t c_OddOneOut = 1;
+      if(ChargeSum==3) c_OddOneOut = 0;
+      //
+      if(c[0]==c_OddOneOut) {Kpair1=Kfactor23; Kpair2=Kfactor24; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor13; Kpair6=Kfactor14;}
+      else if(c[1]==c_OddOneOut) {Kpair1=Kfactor13; Kpair2=Kfactor14; Kpair3=Kfactor34;   Kpair4=Kfactor12; Kpair5=Kfactor23; Kpair6=Kfactor24;}
+      else if(c[2]==c_OddOneOut) {Kpair1=Kfactor12; Kpair2=Kfactor14; Kpair3=Kfactor24;   Kpair4=Kfactor13; Kpair5=Kfactor23; Kpair6=Kfactor34;}
+      else {Kpair1=Kfactor12; Kpair2=Kfactor13; Kpair3=Kfactor23;   Kpair4=Kfactor14; Kpair5=Kfactor24; Kpair6=Kfactor34;}
+      
+      if(term==1){
+       Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
+       C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 );
+       C4 += pow(fc,2)*pow(1-fc,2)*( Kpair4 + Kpair5 + Kpair6 );
+       C4 += pow(fc,3)*(1-fc)*Kpair1*Kpair2*Kpair3;
+       C4 += pow(fc,3)*(1-fc)*(Kpair1*Kpair4*Kpair5 + Kpair2*Kpair4*Kpair6 + Kpair3*Kpair5*Kpair6);// doesn't matter which two MC K's used
+       C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
+       return C4;
+      }else if( (term==2 && ChargeSum==1) || (term==5 && ChargeSum==3)){
+       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+       C3 += pow(fc,2)*(1-fc)*Kpair1;
+       C3 += pow(fc,2)*(1-fc)*Kpair2;
+       C3 += pow(fc,2)*(1-fc)*Kpair3;
+       C3 += pow(fc,3)*Kpair1*Kpair2*Kpair3;
+       return C3;
+      }else if(term<=5){// one SC pair, two MC pairs
+       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+       C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair5;// any MC pair will do
+       C3 += pow(fc,3)*Kpair1*Kpair4*Kpair5;
+       return C3;
+      }else if(term==6 || term==7){
+       if(ChargeSum==1) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
+       else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
+      }else if(term==8){
+       return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
+      }else if(term==9){
+       return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
+      }else if(term==10 || term==1){
+       if(ChargeSum==3) return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
+       else return ((1-fcSq) + fcSq*Kpair4);// any MC pair will do
+      }else return 1.0;
+    }else{// --++ configuration
+      Float_t Kpair1=0, Kpair2=0, Kpair3=0, Kpair4=0, Kpair5=0, Kpair6=0;
+      if(c[0]==c[1]) {Kpair1=Kfactor12; Kpair2=Kfactor34;    Kpair3=Kfactor13; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor24;}
+      else if(c[0]==c[2]) {Kpair1=Kfactor13; Kpair2=Kfactor24;    Kpair3=Kfactor12; Kpair4=Kfactor14; Kpair5=Kfactor23; Kpair6=Kfactor34;}
+      else {Kpair1=Kfactor14; Kpair2=Kfactor23;    Kpair3=Kfactor12; Kpair4=Kfactor13; Kpair5=Kfactor24; Kpair6=Kfactor34;}
+      //
+      if(term==1){
+       Float_t C4 = pow(1-fc,4) + 4*fc*pow(1-fc,3);
+       C4 += pow(fc,2)*pow(1-fc,2)*( Kpair1 + Kpair2 + Kpair3 + Kpair4 + Kpair5 + Kpair6);
+       C4 += pow(fc,3)*(1-fc)*( Kpair1*Kpair3*Kpair4 + Kpair2*Kpair3*Kpair4 + Kpair1*Kpair5*Kpair6 + Kpair2*Kpair5*Kpair6);
+       C4 += pow(fc,4)*Kfactor12*Kfactor13*Kfactor14*Kfactor23*Kfactor24*Kfactor34;
+       return C4;
+      }else if(term<=5){
+       Float_t C3 = pow(1-fc,3) + 3*fc*pow(1-fc,2);
+       C3 += pow(fc,2)*(1-fc)*Kpair1;// any SC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair4;// any MC pair will do
+       C3 += pow(fc,2)*(1-fc)*Kpair6;// any MC pair will do
+       C3 += pow(fc,3)*Kpair1*Kpair4*Kpair6;
+       return C3;
+      }else if(term==6 || term==11){
+       return ((1-fcSq) + fcSq*Kpair1);// any SC pair will do
+      }else if(term !=12){
+       return ((1-fcSq) + fcSq*Kpair3);// any MC pair will do
+      }else return 1.0;
+    }
+  }
+  
+}
+//________________________________________________________________________
+void AliFourPion::SetMomResCorrections(Bool_t legoCase, TH2D *temp2D){
+  
+  if(legoCase){
+    cout<<"LEGO call to SetMomResCorrections"<<endl;
+    fMomResC2 = (TH2D*)temp2D->Clone();
+    fMomResC2->SetDirectory(0);
+  }else {
+    TFile *momResFile = new TFile("MomResFile.root","READ");
+    if(!momResFile->IsOpen()) {
+      cout<<"No momentum resolution file found"<<endl;
+      AliFatal("No momentum resolution file found.  Kill process.");
+    }else {cout<<"Good Momentum Resolution File Found!"<<endl;}
+    
+    TH2D *temp2D2 = (TH2D*)momResFile->Get("MRC_C2_SC");
+    fMomResC2 = (TH2D*)temp2D2->Clone();
+    fMomResC2->SetDirectory(0);
+    
+    momResFile->Close();
+  }
+
+  
+  for(Int_t bx=1; bx<=fMomResC2->GetNbinsX(); bx++){
+    for(Int_t by=1; by<=fMomResC2->GetNbinsY(); by++){
+      if(fMomResC2->GetBinContent(bx,by) > 1.5) fMomResC2->SetBinContent(bx,by, 1.5);// Maximum is ~1.02 
+      if(fMomResC2->GetBinContent(bx,by) < 0.95) fMomResC2->SetBinContent(bx,by, 0.95);// Minimum is ~0.98
+    }
+  }
+
+  cout<<"Done reading momentum resolution file"<<endl;
+}
+//________________________________________________________________________
+void AliFourPion::SetFSICorrelations(Bool_t legoCase, TH1D *tempss[12], TH1D *tempos[12]){
+  // read in 2-particle and 3-particle FSI correlations = K2 & K3
+  // 2-particle input histo from file is binned in qinv.  3-particle in qinv of each pair
+  if(legoCase){
+    cout<<"LEGO call to SetFSICorrelations"<<endl;
+    for(Int_t MB=0; MB<12; MB++) {
+      fFSIss[MB] = (TH1D*)tempss[MB]->Clone();
+      fFSIos[MB] = (TH1D*)tempos[MB]->Clone();
+      //
+      fFSIss[MB]->SetDirectory(0);
+      fFSIos[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;}
+    
+    TH1D *temphistoSS[12];
+    TH1D *temphistoOS[12];
+    for(Int_t MB=0; MB<12; MB++) {
+      TString *nameK2SS = new TString("K2ss_");
+      *nameK2SS += MB;
+      temphistoSS[MB] = (TH1D*)fsifile->Get(nameK2SS->Data());
+      //
+      TString *nameK2OS = new TString("K2os_");
+      *nameK2OS += MB;
+      temphistoOS[MB] = (TH1D*)fsifile->Get(nameK2OS->Data());
+      //
+      fFSIss[MB] = (TH1D*)temphistoSS[MB]->Clone();
+      fFSIos[MB] = (TH1D*)temphistoOS[MB]->Clone();
+      fFSIss[MB]->SetDirectory(0);
+      fFSIos[MB]->SetDirectory(0);
+    }
+    //
+    
+    fsifile->Close();
+  }
+  
+  cout<<"Done reading FSI file"<<endl;
+}
+//________________________________________________________________________
+Float_t AliFourPion::FSICorrelation(Int_t charge1, Int_t charge2, Float_t qinv){
+  // returns 2-particle Coulomb correlations = K2
+  Int_t qbinL = fFSIss[fFSIindex]->GetXaxis()->FindBin(qinv-fFSIss[fFSIindex]->GetXaxis()->GetBinWidth(1)/2.);
+  Int_t qbinH = qbinL+1;
+  if(qbinL <= 0) return 1.0;
+  if(qbinH > fFSIss[fFSIindex]->GetNbinsX()) return 1.0;
+  
+  Float_t slope=0;
+  if(charge1==charge2){
+    slope = fFSIss[fFSIindex]->GetBinContent(qbinL) - fFSIss[fFSIindex]->GetBinContent(qbinH);
+    slope /= fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
+    return (slope*(qinv - fFSIss[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIss[fFSIindex]->GetBinContent(qbinL));
+  }else {
+    slope = fFSIos[fFSIindex]->GetBinContent(qbinL) - fFSIos[fFSIindex]->GetBinContent(qbinH);
+    slope /= fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL) - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinH);
+    return (slope*(qinv - fFSIos[fFSIindex]->GetXaxis()->GetBinCenter(qbinL)) + fFSIos[fFSIindex]->GetBinContent(qbinL));
+  }
+}
+//________________________________________________________________________
+void AliFourPion::SetFillBins2(Int_t c1, Int_t c2, Int_t &b1, Int_t &b2){
+  if((c1+c2)==1) {b1=0; b2=1;}// Re-assign to merge degenerate histos
+  else {b1=c1; b2=c2;}
+}
+//________________________________________________________________________
+void AliFourPion::SetFillBins3(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){
+    
+  // "part" specifies which pair is from the same event.  Only relevant for terms 2-4 
+  Bool_t seSS=kFALSE;
+  if(part==1) {// default case (irrelevant for term 1 and term 5)
+    if(c1==c2) seSS=kTRUE;
+  }
+  if(part==2){
+    if(c1==c3) seSS=kTRUE;
+  }
+  
+  
+  // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
+  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;
+  }
+  
+}
+//________________________________________________________________________
+void AliFourPion::SetFillBins4(Int_t c1, Int_t c2, Int_t c3, Int_t c4, Int_t &b1, Int_t &b2, Int_t &b3, Int_t &b4, Int_t ENsum, Bool_t fillTerm[12]){
+  
+  // fill2, fill3, fill4 are only used for Cumulant Terms 2,3,4
+  if( (c1+c2+c3+c4)==0 || (c1+c2+c3+c4)==4) {// all of the same charge: ---- or ++++
+    
+    b1=c1; b2=c2; b3=c3; b4=c4;
+    if(ENsum==0) fillTerm[0]=kTRUE;
+    else if(ENsum==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
+    else if(ENsum==3) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
+    else fillTerm[11]=kTRUE;
+  
+  }else if( (c1+c2+c3+c4)==1) {// one positive charge: ---+
+  
+    b1=0; b2=0; b3=0; b4=1;// Re-assign to merge degenerate histos
+    if(ENsum==0) fillTerm[0]=kTRUE;
+    else if(ENsum==1){
+      if(c4==1) fillTerm[1]=kTRUE;
+      else {fillTerm[2]=kTRUE; fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
+    }else if(ENsum==3){
+      if(c3==1 || c4==1) {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[8]=kTRUE;} 
+      else {fillTerm[7]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
+    }else fillTerm[11]=kTRUE;
+  
+  }else if( (c1+c2+c3+c4)==2) {// two positive charges: --++
+    
+    b1=0; b2=0; b3=1; b4=1;// Re-assign to merge degenerate histos
+    if(ENsum==0) fillTerm[0]=kTRUE;
+    else if(ENsum==1){
+      if(c4==1) {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE;}
+      else {fillTerm[3]=kTRUE; fillTerm[4]=kTRUE;}
+    }else if(ENsum==3){
+      if( (c1+c2)==0) fillTerm[5]=kTRUE;
+      else if( (c1+c2)==1) {fillTerm[6]=kTRUE; fillTerm[7]=kTRUE; fillTerm[8]=kTRUE; fillTerm[9]=kTRUE;}
+      else fillTerm[10]=kTRUE;
+    }else fillTerm[11]=kTRUE;
+
+  }else{// three positive charges
+    
+    b1=0; b2=1; b3=1; b4=1;// Re-assign to merge degenerate histos
+    if(ENsum==0) fillTerm[0]=kTRUE;
+    else if(ENsum==1){
+      if(c4==0) fillTerm[4]=kTRUE;
+      else {fillTerm[1]=kTRUE; fillTerm[2]=kTRUE; fillTerm[3]=kTRUE;}
+    }else if(ENsum==3){
+      if(c3==0 || c4==0) {fillTerm[8]=kTRUE; fillTerm[9]=kTRUE; fillTerm[10]=kTRUE;}
+      else {fillTerm[5]=kTRUE; fillTerm[6]=kTRUE; fillTerm[7]=kTRUE;}
+    }else fillTerm[11]=kTRUE;
+  
+  }
+  
+}
+//________________________________________________________________________
+void AliFourPion::SetFSIindex(Float_t R){
+  if(!fMCcase){
+    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
+  }else{// FSI binning for MC 
+    if(R>=10.) fFSIindex = 0;
+    else if(R>=9.) fFSIindex = 1;
+    else if(R>=8.) fFSIindex = 2;
+    else if(R>=7.) fFSIindex = 3;
+    else if(R>=6.) fFSIindex = 4;
+    else if(R>=5.) fFSIindex = 5;
+    else if(R>=4.) fFSIindex = 6;
+    else if(R>=3.) fFSIindex = 7;
+    else if(R>=2.) fFSIindex = 8;
+    else fFSIindex = 9;
+  }
+}
diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPion.h
new file mode 100755 (executable)
index 0000000..97ab36b
--- /dev/null
@@ -0,0 +1,322 @@
+#ifndef AliFourPion_cxx
+#define AliFourPion_cxx
+//
+// Class AliFourPion
+//
+// AliFourPion
+// 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 "AliFourPionEventCollection.h"
+#include "AliCentrality.h"
+
+class AliFourPion : public AliAnalysisTaskSE {
+ public:
+
+  AliFourPion();
+  AliFourPion(const Char_t *name);
+  virtual ~AliFourPion();
+  AliFourPion(const AliFourPion &obj); 
+  AliFourPion &operator=(const AliFourPion &obj);
+  
+
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+
+  enum {
+    kPairLimit = 10000,//10000
+    kNormPairLimit = 45000,
+    kMultLimitPbPb = 1800,//1800
+    kMultLimitpp = 300,
+    kMultBinspp = 11,
+    kMCarrayLimit = 150000,// 110000
+    kQbinsWeights = 40,
+    kNDampValues = 16,
+    kRmin = 5,// EW min radii 5 fm
+    kDENtypes = 3,
+  };
+
+  static const Int_t fKbinsT   = 4;// Set fKstep as well !!!!
+  static const Int_t fKbinsY   = 1;// Set fKstep as well !!!!
+  static const Int_t fEDbins   = 2;
+  static const Int_t fCentBins = 10;// 0-50%
+  static const Int_t fRVALUES  = 7;// 7 EW radii (5-11) , was 8 Gaussian radii (3-10fm)
+
+
+  Int_t GetNumKtBins() const {return AliFourPion::fKbinsT;}
+  Int_t GetNumRValues() const {return AliFourPion::fRVALUES;}
+  Int_t GetNumCentBins() const {return AliFourPion::fCentBins;}
+  Int_t GetNumEDBins() const {return AliFourPion::fEDbins;}
+  void SetWeightArrays(Bool_t legoCase=kTRUE, TH3F *histos[AliFourPion::fKbinsT][AliFourPion::fCentBins]=0x0);
+  void SetMomResCorrections(Bool_t legoCase=kTRUE, TH2D *temp2D=0x0);
+  void SetFSICorrelations(Bool_t legoCase=kTRUE, TH1D *tempss[12]=0x0, TH1D *tempos[12]=0x0);
+  //
+  void SetMCdecision(Bool_t mc) {fMCcase = mc;}
+  void SetTabulatePairs(Bool_t tabulate) {fTabulatePairs = tabulate;}
+  void SetPbPbCase(Bool_t pbpb) {fPbPbcase = pbpb;}
+  void SetGenerateSignal(Bool_t gen) {fGenerateSignal = gen;}
+  void SetGeneratorOnly(Bool_t genOnly) {fGeneratorOnly = genOnly;}
+  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;}
+  void SetfcSq(Float_t fcSq) {ffcSq = fcSq;}
+  //
+
+
+ private:
+
+  void ParInit();
+  Bool_t AcceptPair(AliFourPionTrackStruct, AliFourPionTrackStruct);
+  Float_t GamovFactor(Int_t, Int_t, Float_t);
+  void Shuffle(Int_t*, Int_t, Int_t);
+  Float_t GetQinv(Float_t[], Float_t[]);
+  void GetQosl(Float_t[], Float_t[], Float_t&, Float_t&, Float_t&);
+  void GetWeight(Float_t[], Float_t[], Float_t&, Float_t&);
+  Float_t FSICorrelation(Int_t, Int_t, Float_t);
+  Float_t MCWeight(Int_t[2], Float_t, Float_t, Float_t, Float_t);
+  Float_t MCWeightOSL(Int_t, Int_t, Int_t, Int_t, Float_t, Float_t, Float_t, Float_t);
+  Float_t MCWeight3(Int_t, Float_t, Float_t, Int_t[3], Float_t[3], Float_t[3]);
+  Float_t MCWeightFSI3(Int_t, Float_t, Float_t, Int_t[3], Float_t[3]);
+  Float_t MCWeight4(Int_t, Float_t, Float_t, Int_t[4], Float_t[6], Float_t[6]);
+  Float_t MCWeightFSI4(Int_t, Float_t, Float_t, Int_t[4], Float_t[6]);
+  //
+  void SetFillBins2(Int_t, Int_t, Int_t&, Int_t&);
+  void SetFillBins3(Int_t, Int_t, Int_t, Short_t, Int_t&, Int_t&, Int_t&, Bool_t&, Bool_t&, Bool_t&);
+  void SetFillBins4(Int_t, Int_t, Int_t, Int_t, Int_t&, Int_t&, Int_t&, Int_t&, Int_t, Bool_t[12]);
+  void SetFSIindex(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
+  
+  
+  AliFourPionEventCollection ***fEC; //!
+  AliFourPionEventStruct *fEvt; //!
+  AliFourPionTrackStruct *fTempStruct; //!
+  TRandom3* fRandomNumber; //!
+
+  
+  //////////////////////////////
+  // histogram structures
+
+  
+  struct St6 {
+    TH1D *fNorm3; //!
+    TH1D *fTerms3; //!
+    TProfile *fKfactor; //!
+    TProfile *fMeanQinv; //!
+    TH2D *fIdeal; //!
+    TH2D *fSmeared; //!
+    //
+    TH2D *fMuonSmeared; //!
+    TH2D *fMuonIdeal; //!
+    TH2D *fMuonPionK3; //!
+    TH2D *fPionPionK3; //!
+    TH2D *fTwoPartNorm; //!
+  };
+  struct St7 {
+    TH3D *fTerms2OSL; //!
+    TH3D *fTerms2OSLQW; //!
+  };
+  struct St5 {
+    TH2D *fTerms2; //!
+    TH2D *fTerms2QW; //!
+    TH3D *fTerms2ThreeD; //!
+    TProfile2D *fAvgP; //!
+    TH2D *fIdeal; //!
+    TH2D *fSmeared; //!
+    TH2D *fUnitMultBin; //!
+    //
+    TH2D *fMuonSmeared; //!
+    TH2D *fMuonIdeal; //!
+    TH2D *fMuonPionK2; //!
+    TH2D *fPionPionK2; //!
+    TH1D *fMCqinv; //!
+    TH1D *fMCqinvQW; //!
+    TH2D *fPIDpurityDen; //!
+    TH3D *fPIDpurityNum; //!
+    struct St7 OSL_ktbin[2];
+  };
+  struct StFourPT {
+    TH1D *fNorm4; //!
+    TH1D *fTerms4; //!
+    TProfile *fKfactor; //!
+    TH2D *fIdeal; //!
+    TH2D *fSmeared; //!
+    //
+    TH2D *fMuonSmeared; //!
+    TH2D *fMuonIdeal; //!
+    TH2D *fMuonPionK4; //!
+    TH2D *fPionPionK4; //!
+    TH2D *fTwoPartNorm; //!
+  };
+  struct St_EDB {
+    struct St5 TwoPT[2];
+    struct St6 ThreePT[5];
+    struct StFourPT FourPT[5];
+  };
+  struct St_M {
+    struct St_EDB EDB[fEDbins];
+  };
+  struct St4 {
+    struct St_M MB[fCentBins];
+  };
+  struct St3 {
+    struct St4 Charge4[2];
+    struct St_M MB[fCentBins];
+  };
+  struct St2 {
+    struct St3 Charge3[2];
+    struct St_M MB[fCentBins];
+  };
+  struct St1 {
+    struct St2 Charge2[2];
+  };
+  struct St1 Charge1[2];//!
+
+
+  /////////////////////
+  // 4D r3 denominator
+  struct St_Ky {
+    struct St_M MB[fCentBins];
+  };
+  struct St_Kt {
+    struct St_Ky KY[fKbinsY];
+  };
+  struct St_Kt KT[fKbinsT];//!
+  
+  Bool_t fLEGO;
+  Bool_t fMCcase;
+  Bool_t fAODcase;
+  Bool_t fPbPbcase;
+  Bool_t fGenerateSignal;
+  Bool_t fGeneratorOnly;
+  Bool_t fPdensityPairCut;
+  Bool_t fTabulatePairs;
+  Int_t fRMax;
+  Float_t ffcSq;
+  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[kMultBinspp+1];
+  Float_t fQcut;
+  Float_t fQLowerCut;
+  Float_t fNormQcutLow;
+  Float_t fNormQcutHigh;
+  Float_t fKupperBound;
+  Float_t fQupperBoundQ2;
+  Float_t fQupperBoundQ3;
+  Float_t fQupperBoundQ4;
+  Float_t fQbinsQ2;
+  Float_t fQbinsQ3;
+  Float_t fQbinsQ4;
+  Float_t fQupperBoundWeights;
+  Float_t fKstepT[fKbinsT];
+  Float_t fKstepY[fKbinsY];
+  Float_t fKmeanT[fKbinsT];
+  Float_t fKmeanY[fKbinsY];
+  Float_t fKmiddleT[fKbinsT];
+  Float_t fKmiddleY[fKbinsY];
+  Float_t fQstep;
+  Float_t fQstepWeights;
+  Float_t fQmean[kQbinsWeights];
+  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;
+  Int_t fKtIndexL,fKtIndexH;
+  //
+  Int_t fQoIndexL,fQoIndexH;
+  Int_t fQsIndexL,fQsIndexH;
+  Int_t fQlIndexL,fQlIndexH;
+
+  Bool_t fDummyB;
+
+  Float_t fKT3transition;
+  Float_t fKT4transition;
+  
+
+  bool LowQPairSwitch_E0E0[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool LowQPairSwitch_E0E1[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool LowQPairSwitch_E0E2[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool LowQPairSwitch_E0E3[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool LowQPairSwitch_E1E2[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool LowQPairSwitch_E1E3[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool LowQPairSwitch_E2E3[kMultLimitPbPb][kMultLimitPbPb];//!
+  //
+  bool NormQPairSwitch_E0E0[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool NormQPairSwitch_E0E1[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool NormQPairSwitch_E0E2[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool NormQPairSwitch_E0E3[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool NormQPairSwitch_E1E2[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool NormQPairSwitch_E1E3[kMultLimitPbPb][kMultLimitPbPb];//!
+  bool NormQPairSwitch_E2E3[kMultLimitPbPb][kMultLimitPbPb];//!
+
+ public:
+  TH2D *fMomResC2;
+  TH1D *fFSIss[12];
+  TH1D *fFSIos[12];
+  TH3F *fNormWeight[fKbinsT][fCentBins];
+  
+
+  ClassDef(AliFourPion, 1); 
+};
+
+#endif
diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPionEventCollection.cxx b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPionEventCollection.cxx
new file mode 100755 (executable)
index 0000000..8277bff
--- /dev/null
@@ -0,0 +1,284 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+//  This class provides storage for event and track information which 
+//  are used for same-event as well as mixed-event analyses in AliFourPion 
+//
+//  authors: Dhevan Gangadharan (dhevan.raja.gangadharan@cern.ch)
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#include "AliFourPionEventCollection.h"
+
+AliFourPionTrackStruct::AliFourPionTrackStruct():
+  fStatus(0),
+  fFiltermap(0),
+  fId(0),
+  fPhi(0),
+  fPt(0),
+  fMom(0),
+  fP(),
+  fCharge(0),
+  fEta(0),
+  fMass(0),
+  fDCAXY(0),
+  fDCAZ(0),
+  fDCA(0),
+  fEaccepted(0),
+  fKey(0),
+  fClusterMap(0),
+  fSharedMap(0),
+  fX(),
+  fTOFhit(0),
+  fElectron(0),
+  fMuon(0),
+  fPion(0),
+  fKaon(0),
+  fProton(0),
+  fLabel(0)// MC
+{
+  //Default constructor
+}
+AliFourPionTrackStruct::AliFourPionTrackStruct(const AliFourPionTrackStruct &obj)
+  : fStatus(obj.fStatus),
+    fFiltermap(obj.fFiltermap),
+    fId(obj.fId),
+    fPhi(obj.fPhi),
+    fPt(obj.fPt),
+    fMom(obj.fMom),
+    fP(),
+    fCharge(obj.fCharge),
+    fEta(obj.fEta),
+    fMass(obj.fMass),
+    fDCAXY(obj.fDCAXY),
+    fDCAZ(obj.fDCAZ),
+    fDCA(obj.fDCA),
+    fEaccepted(obj.fEaccepted),
+    fKey(obj.fKey),
+    fClusterMap(obj.fClusterMap),
+    fSharedMap(obj.fSharedMap),
+    fX(),
+    fTOFhit(obj.fTOFhit),
+    fElectron(obj.fElectron),
+    fMuon(obj.fMuon),
+    fPion(obj.fPion),
+    fKaon(obj.fKaon),
+    fProton(obj.fProton),
+    fLabel(obj.fLabel)// MC
+{
+  // copy constructor
+}
+AliFourPionTrackStruct &AliFourPionTrackStruct::operator=(const AliFourPionTrackStruct &obj) 
+{
+  // Assignment operator  
+  if (this == &obj)
+    return *this;
+
+  fStatus = obj.fStatus;
+  fFiltermap = obj.fFiltermap;
+  fId = obj.fId;
+  fPhi = obj.fPhi;
+  fPt = obj.fPt;
+  fMom = obj.fMom;
+  fP[0] = obj.fP[0];
+  fP[1] = obj.fP[1];
+  fP[2] = obj.fP[2];
+  fCharge = obj.fCharge;
+  fEta = obj.fEta;
+  fMass = obj.fMass;
+  fDCAXY = obj.fDCAXY;
+  fDCAZ = obj.fDCAZ;
+  fDCA = obj.fDCA;
+  fEaccepted = obj.fEaccepted;
+  fKey = obj.fKey;
+  fClusterMap = obj.fClusterMap;
+  fSharedMap = obj.fSharedMap;
+  fX[0] = obj.fX[0];
+  fX[1] = obj.fX[1];
+  fX[2] = obj.fX[2];
+  fTOFhit = obj.fTOFhit;
+  fElectron = obj.fElectron;
+  fMuon = obj.fMuon;
+  fPion = obj.fPion;
+  fKaon = obj.fKaon;
+  fProton = obj.fProton;
+  fLabel = obj.fLabel;// MC
+
+  return (*this);
+}
+AliFourPionTrackStruct::~AliFourPionTrackStruct()
+{
+  // Destructor
+}
+
+//_____________________________________________________________________________
+AliFourPionMCStruct::AliFourPionMCStruct():
+  fPx(0),
+  fPy(0),
+  fPz(0),
+  fPtot(0),
+  fPdgCode(0),
+  fMotherLabel(0)
+{
+  // Default constructor
+}
+AliFourPionMCStruct::AliFourPionMCStruct(const AliFourPionMCStruct &obj)
+  : fPx(obj.fPx),
+    fPy(obj.fPy),
+    fPz(obj.fPz),
+    fPtot(obj.fPtot),
+    fPdgCode(obj.fPdgCode),
+    fMotherLabel(obj.fMotherLabel)
+{
+  // copy constructor
+}
+AliFourPionMCStruct &AliFourPionMCStruct::operator=(const AliFourPionMCStruct &obj) 
+{
+  // Assignment operator  
+  if (this == &obj)
+    return *this;
+
+  fPx = obj.fPx;
+  fPy = obj.fPy;
+  fPz = obj.fPz;
+  fPtot = obj.fPtot;
+  fPdgCode = obj.fPdgCode;
+  fMotherLabel = obj.fMotherLabel;
+
+  return (*this);
+}
+AliFourPionMCStruct::~AliFourPionMCStruct()
+{
+  // Destructor
+}
+
+//_____________________________________________________________________________
+AliFourPionEventStruct::AliFourPionEventStruct():
+  fFillStatus(0),
+  fNtracks(0),
+  fMCarraySize(0),
+  fTracks(0),
+  fMCtracks(0)
+{
+  // Default constructor
+}
+AliFourPionEventStruct::AliFourPionEventStruct(const AliFourPionEventStruct &obj)
+  : fFillStatus(obj.fFillStatus),
+    fNtracks(obj.fNtracks),
+    fMCarraySize(obj.fMCarraySize),
+    fTracks(obj.fTracks),
+    fMCtracks(obj.fMCtracks)
+{
+  // copy constructor
+}
+AliFourPionEventStruct &AliFourPionEventStruct::operator=(const AliFourPionEventStruct &obj) 
+{
+  // Assignment operator  
+  if (this == &obj)
+    return *this;
+
+  fFillStatus = obj.fFillStatus;
+  fNtracks = obj.fNtracks;
+  fMCarraySize = obj.fMCarraySize;
+  fTracks = obj.fTracks;
+  fMCtracks = obj.fMCtracks;
+  
+  return (*this);
+}
+AliFourPionEventStruct::~AliFourPionEventStruct()
+{
+  // Destructor
+  if(fTracks) delete fTracks;
+  if(fMCtracks) delete fMCtracks;
+}
+
+//_____________________________________________________________________________
+AliFourPionEventCollection::AliFourPionEventCollection():
+  fFIFO(0),
+  fLimit(0),
+  fMCLimit(0),
+  fEvtStr(0)
+{
+  // Default constructor
+}
+AliFourPionEventCollection::AliFourPionEventCollection(Short_t a, Int_t lim, Int_t mcarraylimit, Bool_t MCcase):
+  fFIFO(0),
+  fLimit(0),
+  fMCLimit(0),
+  fEvtStr(0)
+{
+  
+  // Main constructor
+  SetBuffSize(a);
+  
+  fEvtStr = new AliFourPionEventStruct[fFIFO];  //allocate pointer array of type particle_event
+  fLimit = lim;
+  fMCLimit = mcarraylimit;
+
+  for(Int_t ii = 0; ii < fFIFO; ii++){   //Initialize particle table pointers to NULL
+    (fEvtStr + ii)->fNtracks = 0;
+    (fEvtStr + ii)->fFillStatus = 0;
+    (fEvtStr + ii)->fMCarraySize = 0;
+    //
+    (fEvtStr + ii)->fTracks = NULL;
+    (fEvtStr + ii)->fTracks = new AliFourPionTrackStruct[fLimit];
+    if(MCcase) (fEvtStr + ii)->fMCtracks = new AliFourPionMCStruct[fMCLimit];
+    
+  }
+}
+AliFourPionEventCollection::AliFourPionEventCollection(const AliFourPionEventCollection &obj)
+  : fFIFO(obj.fFIFO),
+    fLimit(obj.fLimit),
+    fMCLimit(obj.fMCLimit),
+    fEvtStr(obj.fEvtStr)
+{
+  // copy constructor
+}
+AliFourPionEventCollection &AliFourPionEventCollection::operator=(const AliFourPionEventCollection &obj) 
+{
+  // Assignment operator  
+  if (this == &obj)
+    return *this;
+
+  fFIFO = obj.fFIFO;
+  fLimit = obj.fLimit;
+  fMCLimit = obj.fMCLimit;
+  fEvtStr = obj.fEvtStr;
+  
+  return (*this);
+}
+AliFourPionEventCollection::~AliFourPionEventCollection(){
+
+    for(Int_t i = 0; i < fFIFO; i++){
+
+       if((fEvtStr + i)->fTracks != NULL){
+         delete [] (fEvtStr + i)->fTracks;
+         delete [] (fEvtStr + i)->fMCtracks;
+        }      
+       
+    }
+    
+    delete [] fEvtStr;
+    //remove histos from heap
+
+}
+
+
+//_____________________________________________________________________________
+void AliFourPionEventCollection::FIFOShift(){ //Shift elements in FIFO by one and clear last element in FIFO 
+  
+  
+  for(UShort_t i=fFIFO-1 ; i > 0; i--){
+    for(Int_t j=0; j<(fEvtStr + i-1)->fNtracks; j++) (fEvtStr + i)->fTracks[j] = (fEvtStr + i-1)->fTracks[j];
+    for(Int_t j=0; j<(fEvtStr + i-1)->fMCarraySize; j++) (fEvtStr + i)->fMCtracks[j] = (fEvtStr + i-1)->fMCtracks[j];
+    
+    (fEvtStr + i)->fFillStatus = (fEvtStr + i-1)->fFillStatus;
+    (fEvtStr + i)->fNtracks = (fEvtStr + i-1)->fNtracks;
+    (fEvtStr + i)->fMCarraySize = (fEvtStr + i-1)->fMCarraySize;
+    
+  }// fifo loop
+
+
+  (fEvtStr)->fNtracks=0;
+  (fEvtStr)->fFillStatus=0;
+  (fEvtStr)->fMCarraySize=0;
+}
diff --git a/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPionEventCollection.h b/PWGCF/FEMTOSCOPY/Chaoticity/AliFourPionEventCollection.h
new file mode 100755 (executable)
index 0000000..3da828b
--- /dev/null
@@ -0,0 +1,112 @@
+#ifndef ALIFOURPIONEVENTCOLLECTION
+#define ALIFOURPIONEVENTCOLLECTION
+
+#include <iostream>
+#include <string>
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TBits.h"
+#include "TObject.h"
+#include "TVector2.h"
+#include "AliESDtrack.h"
+
+using namespace std;
+
+
+class AliFourPionTrackStruct{// detector tracks
+  
+ public:
+  AliFourPionTrackStruct();
+  virtual ~AliFourPionTrackStruct();
+  AliFourPionTrackStruct(const AliFourPionTrackStruct &obj); 
+  AliFourPionTrackStruct &operator=(const AliFourPionTrackStruct &obj);
+
+
+  UInt_t fStatus;
+  UInt_t fFiltermap;
+  Int_t fId;
+  Double_t fPhi;
+  Float_t fPt;
+  Float_t fMom;
+  Double_t fP[3];
+  Int_t fCharge;
+  Float_t fEta;
+  Float_t fMass;
+  Float_t fDCAXY;
+  Float_t fDCAZ;
+  Float_t fDCA;
+  Float_t fEaccepted;
+  Short_t fKey;
+  TBits fClusterMap;
+  TBits fSharedMap;
+  Double_t fX[3];
+  Bool_t fTOFhit;
+  Bool_t fElectron;
+  Bool_t fMuon;
+  Bool_t fPion;
+  Bool_t fKaon;
+  Bool_t fProton;
+  Int_t fLabel;// MC
+
+  ClassDef(AliFourPionTrackStruct, 1);
+};
+
+class AliFourPionMCStruct{// MC info
+
+ public:
+  AliFourPionMCStruct();
+  virtual ~AliFourPionMCStruct();
+  AliFourPionMCStruct(const AliFourPionMCStruct &obj); 
+  AliFourPionMCStruct &operator=(const AliFourPionMCStruct &obj);
+
+  Float_t fPx;
+  Float_t fPy;
+  Float_t fPz;
+  Float_t fPtot;
+  Int_t fPdgCode;
+  Int_t fMotherLabel;
+
+  ClassDef(AliFourPionMCStruct, 1);
+};
+
+class AliFourPionEventStruct{// like particle_event
+  
+ public:
+  AliFourPionEventStruct();
+  virtual ~AliFourPionEventStruct();
+  AliFourPionEventStruct(const AliFourPionEventStruct &obj); 
+  AliFourPionEventStruct &operator=(const AliFourPionEventStruct &obj);
+
+
+  Int_t fFillStatus;
+  Int_t fNtracks;
+  Int_t fMCarraySize;
+  AliFourPionTrackStruct *fTracks;
+  AliFourPionMCStruct *fMCtracks;
+
+  ClassDef(AliFourPionEventStruct, 1);
+};
+
+
+
+class AliFourPionEventCollection {
+  
+  public:
+    AliFourPionEventCollection();
+    AliFourPionEventCollection(Short_t,Int_t,Int_t,Bool_t);
+    virtual ~AliFourPionEventCollection();
+    AliFourPionEventCollection(const AliFourPionEventCollection &obj); 
+    AliFourPionEventCollection &operator=(const AliFourPionEventCollection &obj);
+   
+    void FIFOShift();
+    void SetBuffSize(Short_t a){fFIFO = a;}
+    Short_t fFIFO; //Size of the Event Storage buffer.
+    Int_t fLimit; //Max number of tracks
+    Int_t fMCLimit; //Max number of MC tracks
+    AliFourPionEventStruct *fEvtStr;
+
+   ClassDef(AliFourPionEventCollection, 1);
+};
+#endif