ShortPlotFromTreePtBinned.C

#include <TBrowser.h>
#include <TBufferFile.h>
#include <TCanvas.h>
#include <TClass.h>
#include <TMath.h>
#include <TLatex.h>
#include <TObject.h>
#include <TClonesArray.h>
#include <TFile.h>
#include <TTree.h>
#include <TFitResultPtr.h>
#include <TPaveText.h>
#include <TGraph.h>
#include <TMatrixD.h>
#include <TLorentzVector.h>
#include <string>
#include <vector>

//#include "TTreeReader.h"
//#include "Event.h"

# include "TF1.h"
# include "TF2.h"
# include "TProfile.h"
# include "TVirtualFitter.h"
# include "TBackCompFitter.h"
# include "TGraphErrors.h"
# include "TFitter.h"
# include "TMinuit.h"
# include "TRandom.h"
# include <iostream>
# include <fstream>
# include <string>
# include "TH1F.h"
# include "TH1D.h"
# include "TH2F.h"
# include "Scripts/AliAnalysisTaskMyMuonTree_AOD.h"

// Code d'analyse complet

// FUNCTIONS

void ShortPlotFromTreePtBinned();
Double_t FourierV2_WrtInvMass(Double_t *x,Double_t *par);
Double_t FourierV2_WrtInvMassIntegrated(Double_t *x,Double_t *par);
Double_t BackFcnV2(Double_t *x,Double_t *par);
Double_t BackFcnV2Poly(Double_t *x,Double_t *par);
Double_t SignalFcnJPsiV2(Double_t *x,Double_t *par);
Double_t FourierV2(Double_t *x,Double_t *par);
Double_t FourierV5(Double_t *x,Double_t *par);
Double_t TwoCBE2E(Double_t *x,Double_t *par);
Double_t ExpBkg(Double_t *x,Double_t *par);
Double_t JPsiCrystalBallExtended(Double_t *x,Double_t *par);
Double_t Psi2SCrystalBallExtended(Double_t *x,Double_t *par);
TFitResultPtr FittingAllInvMass(const char *histoname, TCanvas *canvas);
TFitResultPtr FittingAllInvMassBin(const char *histoname, TCanvas *canvas, int i);
void FcnCombinedAllMass(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */  );
int GetCent(int cent);
int GetCentPM(int Ntkl, float SPDTrackletsPer, float SPDClustersPer, float V0MPer, int zvtx_idx, int groupnumber);
int GetCentCvetan(float V0MPer);
bool isCentral(int centint);
bool isPeripheral(int centint);
int GetPtBin(double pt);
int GetMassBin(double mass);
Double_t myEtaLow(double x);
Double_t myEtaHigh(double x);
Double_t fitPLL(Double_t *x, Double_t *par);
Double_t fitLin(Double_t *x, Double_t *par);

Double_t mDiff = 0.589188;

// Centrality bins:
//0: 0-1%
//1: 1-3%
//2: 3-5%
//3: 5-10%
//4: 10-15%
//5: 15-20%
//6: 20-30%
//7: 30-40%
//8: 40-50%
//9: 50-60%
//10: 60-70%
//11: 70-80%
//12: 80-90%
//13: 90-100%

// Definition of centrality classes
int CentralLowBound = 0;
int CentralHighBound = 5;
int PeripheralLowBound = 8;
int PeripheralHighBound = 13;

//CentralityEstimator to be used
string centralityMethod = "V0MPercentile";

//Do you want to run an analysis to get the V2 wrt number of tracklets ?
bool isMultiplicityStudy = kFALSE;

//Useless
int NtklCentralLowBound = 6;
int NtklCentralHighBound = 12;
int NtklPeripheralLowBound = 0;
int NtklPeripheralHighBound = 6;

//Define the pT binning, number of bins, min, max
double PtBins[] = {0,2,3,6,8};
const int NbPtBins = 4;
double LowDimuPtCut = 0;
double HighDimuPtCut = 8;

//Define the invariant mass binning, number of bins, min, max, range to be used for fitting the invariant mass, number of bins in invariant mass fit
double MassBins[] = {1.0,1.5,1.9,2.3,2.7,3.0,3.3,3.7,4.1,4.5,5.0};
//double MassBins[] = {1.0,5.0};
const int NbinsInvMass = 10;
double MinInvMass = 1.;
double MaxInvMass = 5.;
double MinInvMassFit = 1.5;
double MaxInvMassFit = 4.5;

const int NbinsDimuInvMass = 400; //Manu 300 ?


//Function to say if a centrality bin falls in the central class
bool isCentral(int centint){
    if((centint>=CentralLowBound)&&(centint<=CentralHighBound))
        return kTRUE;
    else{
        return kFALSE;
    }
}

//Function to say if a centrality bin falls in the peripheral class
bool isPeripheral(int centint){
    if((centint>=PeripheralLowBound)&&(centint<=PeripheralHighBound))
        return kTRUE;
    else{
        return kFALSE;
    }
}
// Matrix used to define centrality classes bounds in terms of #SPDTracklets when estimators were not yet calibrated

TH1D* hStatic3;
Double_t mJpsi =  3.096916;
 Double_t mPsip =  3.686108;
// Double_t ratMass = 1.01;
 Double_t ratSigma = 1.05;
Int_t npfits;

TH1F* hnseg(NULL);
TH1F* V2JPsiTkl(NULL);

double V2_Ext1[NbPtBins] = {0};
double V2_Ext2[NbPtBins] = {0};
double errV2_Ext1[NbPtBins] = {0};
double errV2_Ext2[NbPtBins] = {0};

//Function called to do the Jpsi/tkl analysis
void ShortPlotFromTreePtBinned(){
 //   freopen( "logPlotFromTreeJavier16h_NoDPhiCut_NoConstraintPileUp.txt", "w", stdout );
    TH1::SetDefaultSumw2();
    bool doTracklets = kFALSE; //Do tracklet-tracklet analysis ? Obsolete now
    bool isCMSMethod = kFALSE; //Do CMS way of getting to the V2 ?
    bool doMixedEvents = kTRUE; //Do event mixing ?
    bool KeepOnlyOne = kFALSE; //Useless
    int valueOnlyOne = 3; //Useless
    bool AdditionalCutNtkl = kFALSE; //Do you add a cut on the number of tracklets ? Remove very periph events
    int valueAdditionalCutNtkl = 3; //<=
    int preciseZbinFocus = 10; //Useless
    int preciseCentFocus = 0; //Useless
    bool PtBinned = kTRUE; //Have you defined pT bins ?
    
    
// ************************************
// Définitions de paramètres          *
// ************************************
 
    double ZvtxCut = 10; //Cut on z_vtx (in cm)
    double SigmaZvtxCut = 0.25; //Cut on resolution on z_vtx
    double DPhiCut = 0.005; //0.01 Cut on DeltaPhi of tracklets
    double TklEtaCut  = 1; //Cut on Tracklet Eta (used for centrality estimation when estimators were not calibrated)
    double LowDimuYCut = -4; //Rapidity cuts on dimuon
    double HighDimuYCut = -2.5;
//    double MinMultCentral = 37;
//    double MaxMultPeriph = 23;

    double LowJPsiMass = 3.0; //Useless
    double HighJPsiMass = 3.3;

//    const int NbinsInvMass = 10; //FIXME
//    double SizeBinInvMass = (MaxInvMass-MinInvMass)/NbinsInvMass;

    //Define the binning in Deltaphi
    double MinDeltaPhi = 0;//-TMath::Pi()/2;
    double MaxDeltaPhi = TMath::Pi();//1.5*TMath::Pi();
    const int NbinsDeltaPhi = 6; //12
    double SizeBinDeltaPhi = (MaxDeltaPhi-MinDeltaPhi)/NbinsDeltaPhi;
    int BinZeroLeft = 1;//floor((0-MinDeltaPhi)*NbinsDeltaPhi/(TMath::Pi())); // 2 What is the value of the bin on the left of 0 ? Used for baseline
    
    //Define the binning in DeltaEta
    double MinDeltaEta = 0;
    double MaxDeltaEta = 5;
    const int NbinsDeltaEta = 10;
    double DeltaEtaDimuCut = 1.5;//1.5
    double SizeBinDeltaEta = (MaxDeltaEta-MinDeltaEta)/NbinsDeltaEta;
    
    //Define the binning in DeltaEta for tracklets
   double MinDeltaEtaTKL = -2.4; //1.2  //2.4
    double MaxDeltaEtaTKL = 2.4;
    const int NbinsDeltaEtaTKL = 48; //24
    double DeltaEtaTKLCut = 1.2; //1.2  //1.2
    double SizeBinDeltaEtaTKL = (MaxDeltaEtaTKL-MinDeltaEtaTKL)/NbinsDeltaEtaTKL;
    
    //Define the binning in Deltaphi for tracklets
    const int NbinsDeltaPhiTKL = 48;
    double SizeBinDeltaPhiTKL = (MaxDeltaPhi-MinDeltaPhi)/NbinsDeltaPhiTKL;
    int BinZeroLeftTKL = floor((0-MinDeltaPhi)*NbinsDeltaPhiTKL/(2*TMath::Pi()));
    
    //Define the number of centrality bins and z_vtx bins used in the analysis
    const int NbBinsCent = 14;
    const int NbBinsZvtx = 20;
    
    //Define the radii of SPD layers to get the corrective factor on tracklet Deltaphi
    double R_SPD1 = 3.9;
    double R_SPD2 = 7.6;
    double corrFactorDeltaPhi = R_SPD1/(R_SPD2-R_SPD1);
    
    //Define the list of periods to look at
    
  //  Char_t Group_Period[50] = "Group1";
 //  Char_t *arrayOfPeriods[] = {"Group1_LHC16h","Group1_LHC16j","Group1_LHC16k","Group1_LHC16o","Group1_LHC16p","Group1_LHC17i","Group1_LHC17k","Group1_LHC17l","Group2_LHC17h","Group3_LHC17h","Group4_LHC17k","Group4_LHC18l","Group4_LHC18m","Group4_LHC18o","Group4_LHC18p","Group5_LHC17l","Group5_LHC17m","Group5_LHC17o","Group5_LHC17r","Group5_LHC18c","Group5_LHC18d","Group5_LHC18e","Group5_LHC18f","Group6_LHC18m","Group7_LHC18m","Group8_LHC18m","Group9_LHC18m","Group10_LHC18m","Group11_LHC18m","Group12_LHC18m"};
  //  Char_t *arrayOfPeriods[] = {"Group1_LHC16h","Group1_LHC16j","Group1_LHC16k","Group1_LHC16o","Group1_LHC16p","Group1_LHC17i","Group1_LHC17k","Group1_LHC17l"};
   // Char_t *arrayOfPeriods[] = {"Group1_LHC16h"};
    Char_t *arrayOfPeriods[] = {"Cvetan16r_PhySelTRUEe_Vertexer"};
   // Char_t *arrayOfPeriods[] = {"Group1_LHC16h_CMUL_MULi_Match2TrigMULMLL","Group1_LHC16j_CMUL_MULi_Match2TrigMULMLL","Group1_LHC16o_CMUL_MULi_Match2TrigMULMLL","Group1_LHC16p_CMUL_MULi_Match2TrigMULMLL","Group1_LHC17i_CMUL_MULi_Match2TrigMULMLL","Group4_LHC17k_CMUL_MULi_Match2TrigMULMLL","Group4_LHC18o_CMUL_MULi_Match2TrigMULMLL"};
  //  Char_t *arrayOfPeriods[] = {"Group1_LHC16h","Group1_LHC16j"};
  // Char_t *arrayOfPeriods[] = {"Group8_LHC18m_CvetanPU_OnlyMuonTrackCutsApplied"};
    int numberOfPeriods = sizeof(arrayOfPeriods) / sizeof(arrayOfPeriods[0]);
    
    const double binsCent[6] = {0,1,10,20,40,100};//Useless
    Char_t fileInLoc[500];
    Char_t FolderName[500];
    Char_t CanvasName[500];
    Char_t FitFileName[500];
 //   Char_t AssociateFileName[200];

    //Define the location and name of the .root file storing the plots
    sprintf(FitFileName,"~/../../Volumes/Transcend2/ppAnalysis/Scripts/FitFile_NewAnalysisAllEst_Cvetan16r_TRUEe_Vertexer_NoPsi2s_TailspPb_Eta1.5-5_Q0lin_Elagued_V0MPercentile_0-20_40-100_pt0-2-3-6-8_0612e_pichanged_VWGVarChange_strictcent_rang1.5-4.5_invmassnotinboucle_restrictesSBfit_IAll.root");
//    sprintf(FitFileName,"~/../../Volumes/Transcend2/ppAnalysis/Scripts/Systematics/FitFile_TKL_DPhi10mrad_PercentileMethodSPDTracklets_EtaGap1.2_Zcut10.root");
//    sprintf(AssociateFileName,"~/../../Volumes/Transcend2/ppAnalysis/Scripts/Systematics/AssociateFile_TKL_DPhi10mrad_PercentileMethodSPDTracklets_EtaGap1.2_Zcut10.txt");
   
    //Define the location and name of the folder storing the .pdf plots
    sprintf(FolderName,"~/Desktop/ImagesJavierAnalysis/2021novembre/NewAnalysisAllEst_Cvetan16r_TRUEe_Vertexer_NoPsi2s_TailspPb_Eta1.5-5_Q0lin_Elagued_V0MPercentile_0-20_40-100_pt0-2-3-6-8_0612e_pichanged_VWGVarChange_strictcent_rang1.5-4.5_invmassnotinboucle_restrictesSBfit_IAll");
    
    
    

// *************************
// Initialiser les graphes *
// *************************
    TH1F* ScaledME1a(NULL);
    TH1F* SE1a(NULL);
    TH2F* hsingletrac(NULL);
    TH1I* hnseg;
    TH1F* hnseg2(NULL);
    TH1F* hnseg3(NULL);
    TH1F* hnseg4(NULL);
    TH1F* hnseg5(NULL);
    TH1F* hnseg6(NULL);
    TH1F* hnseg7(NULL);
    TH2F* hnseg8(NULL);
    TH2F* hnsegSigma(NULL);
    TH1F* hnseg8_proj_tampon(NULL);
    TH1F* hnseg8Tkl_proj_tampon(NULL);
    TH1F* hnseg8_proj(NULL);
    TH1F* hnseg8Tkl_proj(NULL);
    TH2F* hnseg8ME(NULL);
    TH1F* hnseg8ME_proj_tampon(NULL);
    TH1F* hnseg8TklME_proj_tampon(NULL);
    TH1D* hnseg8ME_proj(NULL);
    TH1D* hnseg8TklME_proj(NULL);
    TH1F* ME_proj_tampon(NULL);
    TH1F* ProjCopy(NULL);
    TH1F* ProjCopy2(NULL);
    TH1F* ProjCopyTkl(NULL);
    TH1F* ProjCopy2Tkl(NULL);
    TH1F* ME_proj_Tkl_tampon(NULL);
    TH1F* ME_proj(NULL);
    TH1F* SE_proj_tampon(NULL);
    TH1I* SE_proj_Tkl_tampon(NULL);
    TH1F* SE_proj(NULL);
    
    TH1F* Sijk(NULL);
    TH1F* Mijk(NULL);
    TH1F* SoverMijk(NULL);
    TH1F* SoverMik(NULL);
    
    TH1F* Sij(NULL);
    TH1F* Mij(NULL);
    TH1F* SoverMij(NULL);
    TH1F* SoverMi(NULL);
    
    TH1F* Yields[NbBinsCent][NbinsInvMass]{ NULL };
    TH1F* YieldsTkl[NbBinsCent]{ NULL };
    TH1F* YieldsTklShortCMS[NbBinsCent]{ NULL };
    TH1F* YieldsTklLongCMS[NbBinsCent]{ NULL };
    TH2F* Correlations[NbBinsCent][NbBinsZvtx][NbinsInvMass]{ NULL };
    TH2I* CorrelationsTkl[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2I* CorrelationsTklShortCMS[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2I* CorrelationsTklLongCMS[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2I* CorrelationsME[NbBinsCent][NbBinsZvtx][NbinsInvMass]{ NULL };
    TH2I* CorrelationsTklME[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2I* CorrelationsTklMEShortCMS[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2I* CorrelationsTklMELongCMS[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2F* CorrelationsMEScaled[NbBinsCent][NbBinsZvtx][NbinsInvMass]{ NULL };
    TH2F* CorrelationsTklMEScaled[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2F* CorrelationsTklMEScaledShortCMS[NbBinsCent][NbBinsZvtx]{ NULL };
    TH2F* CorrelationsTklMEScaledLongCMS[NbBinsCent][NbBinsZvtx]{ NULL };
    TH1F* Yield_tampon(NULL);
    TH1F* YieldTkl_tampon(NULL);
    TH1F* YieldWrtMass_tampon(NULL);
    TH1F* Yield_allC(NULL);
    TH1F* Yield_Central(NULL);
    TH1F* Yield_Periph(NULL);
    TH1F* Yield_Difference(NULL);
    TH1F* YieldTkl_allC(NULL);
    TH1F* YieldTkl_Central(NULL);
    TH1F* YieldTkl_Periph(NULL);
    TH1F* YieldTkl_Difference(NULL);
    TH1F* YieldTkl_allCShortCMS(NULL);
    TH1F* YieldTkl_CentralShortCMS(NULL);
    TH1F* YieldTkl_PeriphShortCMS(NULL);
    TH1F* YieldTkl_DifferenceShortCMS(NULL);
    TH1F* YieldTkl_allCLongCMS(NULL);
    TH1F* YieldTkl_CentralLongCMS(NULL);
    TH1F* YieldTkl_PeriphLongCMS(NULL);
    TH1F* YieldTkl_DifferenceLongCMS(NULL);
    TH1F* Yield_Central_MassBin[NbinsInvMass] = { NULL };
    TH1F* Yield_Periph_MassBin[NbinsInvMass] = { NULL };
    TH1F* Yield_Difference_MassBin[NbinsInvMass] = { NULL };
    TH1F* Yields_PhiBin[NbBinsCent][NbinsDeltaPhi]{ NULL };
    TH1F* YieldWrtMass_allC[NbinsDeltaPhi]{ NULL };
    TH1F* YieldWrtMass_Central[NbinsDeltaPhi]{ NULL };
    TH1F* YieldWrtMass_Periph[NbinsDeltaPhi]{ NULL };
    
    // Plots added for Pt Binning
    TH1F* YieldsPtBinned[NbBinsCent][NbinsInvMass][NbPtBins]{ NULL };//
    TH2F* YieldsPtBinned1a[NbBinsCent][NbinsInvMass][NbPtBins]{ NULL };//
    TH2I* CorrelationsMEPtBinned[NbBinsCent][NbBinsZvtx][NbinsInvMass][NbPtBins]{ NULL };//
    TH2I* CorrelationsMEPtBinnedBigClasses[3][NbBinsZvtx][NbinsInvMass][NbPtBins]{ NULL };//
    TH2F* CorrelationsPtBinned[NbBinsCent][NbBinsZvtx][NbinsInvMass][NbPtBins]{ NULL };//
    TH2F* CorrelationsMEScaledPtBinned[NbBinsCent][NbBinsZvtx][NbinsInvMass][NbPtBins]{ NULL };//
    TH2F* CorrelationsPtBinnedBigClasses[3][NbBinsZvtx][NbinsInvMass][NbPtBins]{ NULL };//
    TH2F* CorrelationsMEScaledPtBinnedBigClasses[3][NbBinsZvtx][NbinsInvMass][NbPtBins]{ NULL };//
    TH1F* Yield_Central_MassBinPtBinned[NbinsInvMass][NbPtBins] = { NULL };//
    TH1F* Yield_Periph_MassBinPtBinned[NbinsInvMass][NbPtBins] = { NULL };//
    TH1F* Yield_Difference_MassBinPtBinned[NbinsInvMass][NbPtBins] = { NULL };//
    TH1F* Yields_PhiBinPtBinned[NbBinsCent][NbinsDeltaPhi][NbPtBins]{ NULL };//
    TH1F* YieldWrtMass_allCPtBinned[NbinsDeltaPhi][NbPtBins]{ NULL };//
    TH1F* YieldWrtMass_CentralPtBinned[NbinsDeltaPhi][NbPtBins]{ NULL };//
    TH1F* YieldWrtMass_PeriphPtBinned[NbinsDeltaPhi][NbPtBins]{ NULL };//
    TH1F* baselines0PtBinned[NbPtBins]{ NULL };
    TH1F* coefficients0PtBinned[NbPtBins]{ NULL };
    TH1F* coefficients1PtBinned[NbPtBins]{ NULL };
    TH1F* coefficients2PtBinned[NbPtBins]{ NULL };
    TH1F* c2b0PtBinned[NbPtBins]{ NULL };
    TH1F* V2JPsiTklPtBinned[NbPtBins]{ NULL };
    
    TH1F* Yields_Central_1PtBinned[NbPtBins]{ NULL };
       TH1F* Yields_Periph_1PtBinned[NbPtBins]{ NULL };
       TH1F* Yields_Difference_1PtBinned[NbPtBins]{ NULL };
    
    
    // Try plots 2D DetaDphi
    
    TH2F* YCentral(NULL);
    TH2F* YPeriph(NULL);
    TH2F* YDifference(NULL);
    TH2F* YTklCentral(NULL);
    TH2I* YTklCentralME(NULL);
    TH2F* YTklPeriph(NULL);
    TH2I* YTklPeriphME(NULL);
    TH2F* YTklDifference(NULL);
    TH1D* YTklCentral_proj_tampon(NULL);
    TH1D* Correl_tampon(NULL);
    TH1D* YTklPeriph_proj_tampon(NULL);
    TH1D* YTklDifference_proj_tampon(NULL);
    TH1D* YTklCentralME_proj_tampon(NULL);
    TH1D* YTklPeriphME_proj_tampon(NULL);
    
    ScaledME1a = new TH1F("ScaledME1a",
                      "ScaledME1a",
                      NbBinsZvtx, -1*int(ZvtxCut), int(ZvtxCut));
    ScaledME1a->SetXTitle("Zvtx");
    SE1a = new TH1F("SE1a",
                      "SE1a",
                      NbBinsZvtx, -1*int(ZvtxCut), int(ZvtxCut));
    SE1a->SetXTitle("Zvtx");
    
    hsingletrac = new TH2F("hsingletrac",
                      "Single tracklet #eta,#phi distribution",
                           NbinsDeltaPhi,0,MaxDeltaPhi-MinDeltaPhi,40,-2,2);
    hsingletrac->SetXTitle("Tkl #phi (rad)");
    hsingletrac->SetYTitle("Tkl #eta");
    
    YCentral = new TH2F("YCentral",
                      "Yield dimuon-tkl #Delta#eta wrt #Delta#phi - Central",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
    YCentral->SetXTitle("Correlation #Delta#phi (rad)");
    YCentral->SetYTitle("Correlation #Delta#eta");
    YPeriph = new TH2F("YPeriph",
                      "Yield dimuon-tkl #Delta#eta wrt #Delta#phi - Periph",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
    YPeriph->SetXTitle("Correlation #Delta#phi (rad)");
    YPeriph->SetYTitle("Correlation #Delta#eta");
    YDifference = new TH2F("YDifference",
                      "Yield dimuon-tkl #Delta#eta wrt #Delta#phi - Difference",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
    YDifference->SetXTitle("Correlation #Delta#phi (rad)");
    YDifference->SetYTitle("Correlation #Delta#eta");
    YTklCentral = new TH2F("YTklCentral",
                      "Yield #Delta#eta wrt #Delta#phi - Tkl Central",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
    YTklCentral->SetXTitle("Correlation #Delta#phi (rad)");
    YTklCentral->SetYTitle("Correlation #Delta#eta");
    YTklPeriph = new TH2F("YTklPeriph",
                      "Yield #Delta#eta wrt #Delta#phi - Tkl Periph",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
    YTklPeriph->SetXTitle("Correlation #Delta#phi (rad)");
    YTklPeriph->SetYTitle("Correlation #Delta#eta");
    YTklDifference = new TH2F("YTklDifference",
                      "Yield #Delta#eta wrt #Delta#phi - Tkl Difference",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
    YTklDifference->SetXTitle("Correlation #Delta#phi (rad)");
    YTklDifference->SetYTitle("Correlation #Delta#eta");
    YTklCentralME = new TH2I("YTklCentralME",
                      "Yield #Delta#eta wrt #Delta#phi - Tkl Central ME",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
    YTklCentralME->SetXTitle("Correlation #Delta#phi (rad)");
    YTklCentralME->SetYTitle("Correlation #Delta#eta");
    YTklPeriphME = new TH2I("YTklPeriphME",
                      "Yield #Delta#eta wrt #Delta#phi - Tkl Periph ME",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
    YTklPeriphME->SetXTitle("Correlation #Delta#phi (rad)");
    YTklPeriphME->SetYTitle("Correlation #Delta#eta");
    // AE
    
    TH2F* hnseg9(NULL);
    TH2F* hnseg10(NULL);
    TH1F* hDPhi(NULL);
    TH2F* hPtWrtMassInv[3]{NULL};
    TH1D *hPtWrtMassInvSliced[3][NbPtBins]= {NULL};
    TH1D *hPtWrtMassInvSlicedRebinned[3][NbPtBins]= {NULL};
    
    TH1F* CentV0M(NULL);
    TH1F* CentTKL(NULL);
    TH1F* CentCL0(NULL);
    TH1F* CentCL1(NULL);
    TH1F* CentSPDTracklets(NULL);
    TH1F* CentSPDClusters(NULL);
    TH2F* CentV0Mwrttkl(NULL);
    TH2F* CentTKLwrttkl(NULL);
    TH2F* CentCL0wrttkl(NULL);
    TH2F* CentCL1wrttkl(NULL);
    TH2F* CentSPDTrackletswrttkl(NULL);
    TH2F* CentSPDClusterswrttkl(NULL);
    
    TH1F* Yields_Central_1(NULL);
    TH1F* Yields_Periph_1(NULL);
    TH1F* Yields_Difference_1(NULL);
    
    TH1F* coefficients0(NULL);
    TH1F* coefficients1(NULL);
    TH1F* coefficients2(NULL);
    TH1F* baselines0(NULL);
    TH1F* c2b0(NULL);
    
    ProjCopy = new TH1F("ProjCopy",
                      "ProjCopy",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    ProjCopy->SetXTitle("Correlation #Delta#phi (rad)");
    ProjCopy->SetYTitle("Correlation #Delta#eta");
    ProjCopy2 = new TH1F("ProjCopy2",
                      "ProjCopy2",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    ProjCopy2->SetXTitle("Correlation #Delta#phi (rad)");
    ProjCopy2->SetYTitle("Correlation #Delta#eta");
    
    ProjCopyTkl = new TH1F("ProjCopyTkl",
                      "ProjCopyTkl",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    ProjCopyTkl->SetXTitle("Correlation #Delta#phi (rad)");
    ProjCopyTkl->SetYTitle("Correlation #Delta#eta");
    ProjCopy2Tkl = new TH1F("ProjCopy2Tkl",
                      "ProjCopy2Tkl",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    ProjCopy2Tkl->SetXTitle("Correlation #Delta#phi (rad)");
    ProjCopy2Tkl->SetYTitle("Correlation #Delta#eta");
    
    
    Yields_Central_1 = new TH1F("Yields_Central_1",
                     "Yield of J/#psi-tkl in Central collisions wrt #Delta#phi",
                     NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Yields_Central_1->SetXTitle("#Delta#phi of the correlation (rad)");
    Yields_Central_1->SetYTitle("Yield_{Central}");
    
    Yields_Periph_1 = new TH1F("Yields_Periph_1",
                     "Yield of J/#psi-tkl in Periph collisions wrt #Delta#phi",
                     NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Yields_Periph_1->SetXTitle("#Delta#phi of the correlation (rad)");
    Yields_Periph_1->SetYTitle("Yield_{Periph}");
    
    Yields_Difference_1 = new TH1F("Yields_Difference_1",
                     "Subtracted yield of J/#psi-tkl (Central-Periph) wrt #Delta#phi",
                     NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Yields_Difference_1->SetXTitle("#Delta#phi of the correlation (rad)");
    Yields_Difference_1->SetYTitle("Yield_{Subtracted}");
    
    coefficients0 = new TH1F("coefficients0",
                        "0^{th} Fourier coefficient wrt dimuon mass",
                        NbinsInvMass,MassBins);
       coefficients0->SetXTitle("Mass of dimuon (GeV/c^{2})");
       coefficients0->SetYTitle("0^{th} Fourier coefficient");
    coefficients1 = new TH1F("coefficients1",
                           "1^{st} Fourier coefficient wrt dimuon mass",
                           NbinsInvMass,MassBins);
          coefficients1->SetXTitle("Mass of dimuon (GeV/c^{2})");
          coefficients1->SetYTitle("1^{st} Fourier coefficient");
    coefficients2 = new TH1F("coefficients2",
                           "2^{nd} Fourier coefficient wrt dimuon mass",
                           NbinsInvMass,MassBins);
          coefficients2->SetXTitle("Mass of dimuon (GeV/c^{2})");
          coefficients2->SetYTitle("2^{nd} Fourier coefficient");
    baselines0 = new TH1F("baselines0",
                     "Baseline_{Periph} wrt dimuon mass",
                     NbinsInvMass,MassBins);
    baselines0->SetXTitle("Mass of dimuon (GeV/c^{2})");
    baselines0->SetYTitle("Baseline_{Periph}");
    c2b0 = new TH1F("c2b0",
                     "2^{nd} Fourier coefficient + Baseline_{Periph} wrt dimuon mass",
                     NbinsInvMass,MassBins);
    c2b0->SetXTitle("Mass of dimuon (GeV/c^{2})");
    c2b0->SetYTitle("2^{nd} Fourier coefficient + Baseline_{Periph}");
    V2JPsiTkl = new TH1F("V2JPsiTkl",
                     "V_{2,J/#psi-tkl} wrt dimuon mass",
                     NbinsInvMass,MassBins);
    V2JPsiTkl->SetXTitle("Mass of dimuon (GeV/c^{2})");
    V2JPsiTkl->SetYTitle("V_{2,J/#psi-tkl}");
    
    char hname[200];
    char hname1[200];
    char hname2[200];
    
    for (int j=0; j <NbinsInvMass; j++){
    sprintf(hname1,"Projected yield in Mass Bin %f GeV to %f GeV - Central",MassBins[j],MassBins[j+1]);
    sprintf(hname2,"Projected yield in Mass Bin, #m_{#mu#mu} #in [%.2f,%.2f] GeV/c^{2} - Central",MassBins[j],MassBins[j+1]);
     Yield_Central_MassBin[j] = new TH1F(hname1, hname2,NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    }
       
       for (int j=0; j <NbinsInvMass; j++){
       sprintf(hname1,"Projected yield in Mass Bin %f GeV to %f GeV - Periph",MassBins[j],MassBins[j+1]);
          sprintf(hname2,"Projected yield in Mass Bin, #m_{#mu#mu} #in [%.2f,%.2f] GeV/c^{2} - Periph",MassBins[j],MassBins[j+1]);
          Yield_Periph_MassBin[j] = new TH1F(hname1, hname2,NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
          }
       
       for (int j=0; j <NbinsInvMass; j++){
           sprintf(hname1,"Projected yield in Mass Bin %f GeV to %f GeV - Difference",MassBins[j],MassBins[j+1]);
       sprintf(hname2,"Projected yield in Mass Bin, #m_{#mu#mu} #in [%.2f,%.2f] GeV/c^{2} - Difference",MassBins[j],MassBins[j+1]);
       Yield_Difference_MassBin[j] = new TH1F(hname1, hname2,NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
       }
    
    TH1I* InvMass_Central(NULL);
    TH1I* InvMass_Periph(NULL);
    TH1I* InvMass_Neither(NULL);
    
//    TH1F* YieldWrtMass_PhiBin_Central[NbinsDeltaPhi] = { NULL };
//    TH1F* YieldWrtMass_PhiBin_Periph[NbinsDeltaPhi] = { NULL };
//
//    for (int j=0; j<NbinsDeltaPhi; j++){
//        sprintf(hname,"Yield wrt Mass, Phi: %d pi/6 to %d pi/6 - Central",j-3,j-2);
//        sprintf(hname2,"YieldWrtMassPhiBin_Central_%d",j);
//        YieldWrtMass_PhiBin_Central[j] = new TH1F(hname2, hname,NbinsInvMass,MinInvMass,MaxInvMass);
//    }
//
//    for (int j=0; j<NbinsDeltaPhi; j++){
//        sprintf(hname,"Yield wrt Mass, Phi: %d pi/6 to %d pi/6 - Periph",j-3,j-2);
//        sprintf(hname2,"YieldWrtMassPhiBin_Periph_%d",j);
//        YieldWrtMass_PhiBin_Periph[j] = new TH1F(hname2, hname,NbinsInvMass,MinInvMass,MaxInvMass);
//    }
    
    
    for(int i=0; i<NbBinsCent; i++){
        for(int k=0; k<NbinsInvMass; k++){
            char hname[200];
            sprintf(hname,"Yields %d %d ",i,k);
            Yields[i][k] = new TH1F(hname,
                              "Yields Correlation wrt #Delta#phi",
                              NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
            Yields[i][k]->SetXTitle("Correlation #Delta#phi (rad)");
            for(int j=0; j<NbBinsZvtx; j++){
                sprintf(hname,"Correlations %d %d %d ",i,j,k);
                Correlations[i][j][k] = new TH2F(hname,
                                  "Correlation #Delta#eta wrt #Delta#phi",
                                  NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                Correlations[i][j][k]->SetXTitle("Correlation #Delta#phi (rad)");
                Correlations[i][j][k]->SetYTitle("Correlation #Delta#eta");
                sprintf(hname,"CorrelationsME %d %d %d ",i,j,k);
                CorrelationsME[i][j][k] = new TH2I(hname,
                                  "MixedEvent Correlation #Delta#eta wrt #Delta#phi",
                                  NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                CorrelationsME[i][j][k]->SetXTitle("Correlation #Delta#phi (rad)");
                CorrelationsME[i][j][k]->SetYTitle("Correlation #Delta#eta");
                sprintf(hname,"CorrelationsMEScaled %d %d %d ",i,j,k);
                CorrelationsMEScaled[i][j][k] = new TH2F(hname,
                                                 "MixedEvent Correlation #Delta#eta wrt #Delta#phi - Scaled",
                                                 NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                               CorrelationsMEScaled[i][j][k]->SetXTitle("Correlation #Delta#phi (rad)");
                               CorrelationsMEScaled[i][j][k]->SetYTitle("Correlation #Delta#eta");
            }
        }
        
    }
    
    for(int i=0; i<NbBinsCent; i++){
        char hname[200];
        sprintf(hname,"Yields Tkl %d ",i);
        YieldsTkl[i] = new TH1F(hname,
                          "Yields Correlation Tkl wrt #Delta#phi",
                          NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
        YieldsTkl[i]->SetXTitle("Correlation #Delta#phi (rad)");
        sprintf(hname,"Yields Tkl ShortCMS %d ",i);
        YieldsTklShortCMS[i] = new TH1F(hname,
                          "Yields Correlation Tkl ShortCMS wrt #Delta#phi",
                          NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
        YieldsTklShortCMS[i]->SetXTitle("Correlation #Delta#phi (rad)");
        sprintf(hname,"Yields Tkl LongCMS %d ",i);
        YieldsTklLongCMS[i] = new TH1F(hname,
                          "Yields Correlation Tkl LongCMS wrt #Delta#phi",
                          NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
        YieldsTklLongCMS[i]->SetXTitle("Correlation #Delta#phi (rad)");
        for(int j=0; j<NbBinsZvtx; j++){
            sprintf(hname,"Correlations Tkl %d %d ",i,j);
            CorrelationsTkl[i][j] = new TH2I(hname,
                              "Correlation Tkl #Delta#eta wrt #Delta#phi",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTkl[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTkl[i][j]->SetYTitle("Correlation #Delta#eta");
            sprintf(hname,"CorrelationsME Tkl %d %d ",i,j);
            CorrelationsTklME[i][j] = new TH2I(hname,
                              "MixedEvent Correlation Tkl #Delta#eta wrt #Delta#phi",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklME[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklME[i][j]->SetYTitle("Correlation #Delta#eta");
            sprintf(hname,"CorrelationsMEScaled Tkl %d %d ",i,j);
            CorrelationsTklMEScaled[i][j] = new TH2F(hname,
                              "MixedEvent Correlation Tkl #Delta#eta wrt #Delta#phi - Scaled",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklMEScaled[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklMEScaled[i][j]->SetYTitle("Correlation #Delta#eta");
            
            sprintf(hname,"Correlations Tkl ShortCMS %d %d ",i,j);
            CorrelationsTklShortCMS[i][j] = new TH2I(hname,
                              "Correlation Tkl ShortCMS #Delta#eta wrt #Delta#phi",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklShortCMS[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklShortCMS[i][j]->SetYTitle("Correlation #Delta#eta");
            sprintf(hname,"CorrelationsME Tkl ShortCMS %d %d ",i,j);
            CorrelationsTklMEShortCMS[i][j] = new TH2I(hname,
                              "MixedEvent Correlation Tkl ShortCMS #Delta#eta wrt #Delta#phi",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklMEShortCMS[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklMEShortCMS[i][j]->SetYTitle("Correlation #Delta#eta");
            sprintf(hname,"CorrelationsME Tkl ShortCMS Scaled %d %d ",i,j);
            CorrelationsTklMEScaledShortCMS[i][j] = new TH2F(hname,
                              "MixedEvent Correlation Tkl ShortCMS #Delta#eta wrt #Delta#phi - Scaled",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklMEScaledShortCMS[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklMEScaledShortCMS[i][j]->SetYTitle("Correlation #Delta#eta");
            
            sprintf(hname,"Correlations Tkl LongCMS %d %d ",i,j);
            CorrelationsTklLongCMS[i][j] = new TH2I(hname,
                              "Correlation Tkl LongCMS #Delta#eta wrt #Delta#phi",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklLongCMS[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklLongCMS[i][j]->SetYTitle("Correlation #Delta#eta");
            sprintf(hname,"CorrelationsME Tkl LongCMS %d %d ",i,j);
            CorrelationsTklMELongCMS[i][j] = new TH2I(hname,
                              "MixedEvent Correlation Tkl LongCMS #Delta#eta wrt #Delta#phi",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklMELongCMS[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklMELongCMS[i][j]->SetYTitle("Correlation #Delta#eta");
            sprintf(hname,"CorrelationsME Tkl LongCMS Scaled %d %d ",i,j);
            CorrelationsTklMEScaledLongCMS[i][j] = new TH2F(hname,
                              "MixedEvent Correlation Tkl LongCMS #Delta#eta wrt #Delta#phi - Scaled",
                              NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
            CorrelationsTklMEScaledLongCMS[i][j]->SetXTitle("Correlation #Delta#phi (rad)");
            CorrelationsTklMEScaledLongCMS[i][j]->SetYTitle("Correlation #Delta#eta");
        }
    }
    
    for(int i=0; i<NbBinsCent; i++){
           for(int p=0; p<NbinsDeltaPhi; p++){
                   char hname[200];
                   sprintf(hname,"Yields PhiBin %d %d ",i,p);
                   Yields_PhiBin[i][p] = new TH1F(hname,
                                     "Yields Correlation wrt dimuon mass",
                                     NbinsInvMass,MassBins);
                   Yields_PhiBin[i][p]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
           }
    }
        
    for(int p=0; p<NbinsDeltaPhi; p++){
            char hname[200];
            char hname2[200];
            sprintf(hname,"Yields PhiBin %d All C",p);
            YieldWrtMass_allC[p] = new TH1F(hname,
                              "Yields dimuon-tkl wrt dimuon mass, all C",
                              NbinsInvMass,MassBins);
            YieldWrtMass_allC[p]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
        
        sprintf(hname,"Yields PhiBin %d Periph",p);
        sprintf(hname2,"Yields dimuon-tkl wrt wrt dimuon mass, Periph, #Delta#phi #in [#frac{%d#pi}{6},#frac{%d#pi}{6}]",p-3, p-2);
        YieldWrtMass_Periph[p] = new TH1F(hname,
                          hname2,
                          NbinsInvMass,MassBins);
        YieldWrtMass_Periph[p]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
        
        sprintf(hname,"Yields PhiBin %d Central",p);
        sprintf(hname2,"Yields dimuon-tkl wrt wrt dimuon mass, Central, #Delta#phi #in [#frac{%d#pi}{6},#frac{%d#pi}{6}]",p-3, p-2);
        YieldWrtMass_Central[p] = new TH1F(hname,
                          hname2,
                          NbinsInvMass,MassBins);
        YieldWrtMass_Central[p]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
    }
        
    
    Yield_allC = new TH1F("Yield_allC",
                         "Yields dimuon-tkl wrt #Delta#phi, all C, all mass",
                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
       Yield_allC->SetXTitle("#Delta#phi (rad)");
        Yield_allC->SetYTitle("Yields");
    Yield_Central = new TH1F("Yield_Central",
                      "Yields dimuon-tkl wrt #Delta#phi, Central, 3.0-3.25 GeV",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Yield_Central->SetXTitle("#Delta#phi (rad)");
    Yield_Central->SetYTitle("Yields_{Central}");
    Yield_Periph = new TH1F("Yield_Periph",
                      "Yields dimuon-tkl wrt #Delta#phi, Periph, 3.0-3.25 GeV",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Yield_Periph->SetXTitle("#Delta#phi (rad)");
    Yield_Periph->SetYTitle("Yields_{Periph}");
    Yield_Difference = new TH1F("Yield_Difference",
                      "Yields dimuon-tkl wrt #Delta#phi, Central-Periph, all mass",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Yield_Difference->SetXTitle("#Delta#phi (rad)");
    Yield_Difference->SetXTitle("Yields_{Subtracted}");
       Yield_tampon = new TH1F("Yield_tampon",
                         "Yields dimuon-tkl wrt #Delta#phi, all C, all mass",
                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
       Yield_tampon->SetXTitle("#Delta#phi (rad)");
        Yield_tampon->SetYTitle("Yields");
    YieldWrtMass_tampon = new TH1F("YieldWrtMass_tampon",
                      "Yields dimuon-tkl wrt dimuon mass",
                      NbinsInvMass,MassBins);
    YieldWrtMass_tampon->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
    YieldWrtMass_tampon->SetYTitle("Yields");
    
     YieldTkl_allC = new TH1F("YieldTkl_allC",
                            "Yields tkl-tkl wrt #Delta#phi, all C, all mass",
                            NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
          YieldTkl_allC->SetXTitle("#Delta#phi (rad)");
       YieldTkl_allC->SetXTitle("Yields");
       YieldTkl_Central = new TH1F("YieldTkl_Central",
                         "Yields tkl-tkl wrt #Delta#phi, Central",
                         NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
       YieldTkl_Central->SetXTitle("#Delta#phi (rad)");
       YieldTkl_Central->SetYTitle("Yields_{Central}");
       YieldTkl_Periph = new TH1F("YieldTkl_Periph",
                         "Yields tkl-tkl wrt #Delta#phi, Periph",
                         NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
       YieldTkl_Periph->SetXTitle("#Delta#phi (rad)");
       YieldTkl_Periph->SetYTitle("Yields_{Periph}");
       YieldTkl_Difference = new TH1F("YieldTkl_Difference",
                         "Yields tkl-tkl wrt #Delta#phi, Central-Periph, all mass",
                         NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
       YieldTkl_Difference->SetXTitle("#Delta#phi (rad)");
       YieldTkl_Difference->SetYTitle("Yields_{Subtracted}");
       YieldTkl_tampon = new TH1F("YieldTkl_tampon",
                         "Yields tkl-tkl wrt #Delta#phi, all C, all mass",
                         NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
       YieldTkl_tampon->SetXTitle("#Delta#phi (rad)");
       YieldTkl_tampon->SetYTitle("Yields");
    
    YieldTkl_allCShortCMS = new TH1F("YieldTkl_allCShortCMS",
                         "Yields Correlation Tkl ShortCMS wrt #Delta#phi, all C, all mass",
                         NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
       YieldTkl_allCShortCMS->SetXTitle("Correlation #Delta#phi (rad)");
    YieldTkl_CentralShortCMS = new TH1F("YieldTkl_CentralShortCMS",
                      "Yields Correlation Tkl ShortCMS wrt #Delta#phi, Central",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    YieldTkl_CentralShortCMS->SetXTitle("Correlation #Delta#phi (rad)");
    YieldTkl_PeriphShortCMS = new TH1F("YieldTkl_PeriphShortCMS",
                      "Yields Correlation Tkl ShortCMS wrt #Delta#phi, Periph",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    YieldTkl_PeriphShortCMS->SetXTitle("Correlation #Delta#phi (rad)");
    YieldTkl_DifferenceShortCMS = new TH1F("YieldTkl_DifferenceShortCMS",
                      "Yields Correlation Tkl ShortCMS wrt #Delta#phi, Central-Periph, all mass",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    YieldTkl_DifferenceShortCMS->SetXTitle("Correlation #Delta#phi (rad)");
    YieldTkl_allCLongCMS = new TH1F("YieldTkl_allCLongCMS",
                         "Yields Correlation Tkl LongCMS wrt #Delta#phi, all C, all mass",
                         NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
       YieldTkl_allCLongCMS->SetXTitle("Correlation #Delta#phi (rad)");
    YieldTkl_CentralLongCMS = new TH1F("YieldTkl_CentralLongCMS",
                      "Yields Correlation Tkl LongCMS wrt #Delta#phi, Central",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    YieldTkl_CentralLongCMS->SetXTitle("Correlation #Delta#phi (rad)");
    YieldTkl_PeriphLongCMS = new TH1F("YieldTkl_PeriphLongCMS",
                      "Yields Correlation Tkl LongCMS wrt #Delta#phi, Periph",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    YieldTkl_PeriphLongCMS->SetXTitle("Correlation #Delta#phi (rad)");
    YieldTkl_DifferenceLongCMS = new TH1F("YieldTkl_DifferenceLongCMS",
                      "Yields Correlation Tkl LongCMS wrt #Delta#phi, Central-Periph, all mass",
                      NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    YieldTkl_DifferenceLongCMS->SetXTitle("Correlation #Delta#phi (rad)");
    
    
   char haxis[200];
    sprintf(haxis, "Count within bin of %d MeV/c^{2}", 1000*(MaxInvMass-MinInvMass)/NbinsDimuInvMass);
        hnseg = new TH1I("hnseg",
                         "Invariant mass of dimuon",
                         NbinsDimuInvMass,MinInvMass,MaxInvMass);
        hnseg->SetXTitle("Mass of dimuon (GeV/c^{2})");
        hnseg->SetYTitle(haxis);
    hDPhi = new TH1F("hDPhi",
                     "#Delta#Phi of Tracklet",
                     10000,-TMath::Pi(),TMath::Pi());
    hDPhi->SetXTitle("Tracklet #Delta#Phi (rad)");
    hDPhi->SetYTitle("Count");
    InvMass_Central = new TH1I("InvMass_Central",
                     "Invariant mass of dimuon - Central",
                     NbinsDimuInvMass,MinInvMass,MaxInvMass);
    InvMass_Central->SetXTitle("Mass of dimuon (GeV/c^{2})");
    InvMass_Central->SetYTitle(haxis);
    InvMass_Periph = new TH1I("InvMass_Periph",
                     "Invariant mass of dimuon - Periph",
                     NbinsDimuInvMass,MinInvMass,MaxInvMass);
    InvMass_Periph->SetXTitle("Mass of dimuon (GeV/c^{2})");
    InvMass_Periph->SetYTitle(haxis);
    InvMass_Neither = new TH1I("InvMass_Neither",
                     "Invariant mass of dimuon - Neither",
                     NbinsDimuInvMass,MinInvMass,MaxInvMass);
    InvMass_Neither->SetXTitle("Mass of dimuon (GeV/c^{2})");
    InvMass_Neither->SetYTitle(haxis);
    hnseg2 = new TH1F("hnseg2",
                     "p_{T} of dimuon",
                     100,0,10);
    hnseg2->SetXTitle("Dimuon p_{T} (GeV/c)");
    hnseg2->SetYTitle("Count");
    hnseg3 = new TH1F("hnseg3",
                     "#eta of dimuon",
                     1000,-10,0);
    hnseg3->SetXTitle("Dimuon #eta");
    hnseg3->SetYTitle("Count");
    hnseg4 = new TH1F("hnseg4",
                      "y of dimuon",
                      1000,-2.4,-4.2);
    hnseg4->SetXTitle("Dimuon y");
    hnseg4->SetYTitle("Count");
    hnseg5 = new TH1F("hnseg5",
                      "#phi of dimuon",
                      1000,0, 2*TMath::Pi());
    hnseg5->SetXTitle("Dimuon #phi (rad)");
    hnseg5->SetYTitle("Count");
    hnseg6 = new TH1F("hnseg6",
                      "Correlation #Delta#phi",
                      1000,MinDeltaPhi,MaxDeltaPhi);
    hnseg6->SetXTitle("Correlation #Delta#phi (rad)");
    hnseg6->SetYTitle("Count");
    hnseg7 = new TH1F("hnseg7",
                      "Correlation #Delta#eta",
                      1000,MinDeltaEta,MaxDeltaEta);
    hnseg7->SetXTitle("Correlation #Delta#eta");
    hnseg7->SetYTitle("Count");
    hnsegSigma = new TH2F("hnsegSigma",
                      "Z_{vtx} resolution wrt NContributors",
                      60,0,60,1000,0,10);
    hnsegSigma->SetXTitle("NContributors");
    hnsegSigma->SetYTitle("Z_{vtx} Resolution (SPD)");
    hnseg8 = new TH2F("hnseg8",
                      "Correlation #Delta#eta wrt #Delta#phi",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
    hnseg8->SetXTitle("Correlation #Delta#phi (rad)");
    hnseg8->SetYTitle("Correlation #Delta#eta");
    hnseg8ME = new TH2F("hnseg8ME",
                      "MixedEvent Correlation #Delta#eta wrt #Delta#phi",
                      NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
    hnseg8ME->SetXTitle("Correlation #Delta#phi (rad)");
    hnseg8ME->SetYTitle("Correlation #Delta#eta");
    hnseg8_proj = new TH1F("hnseg8_proj", "Correlation wrt #Delta#eta, projection",NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
    hnseg8_proj->SetXTitle("Correlation #Delta#eta");
    hnseg8_proj->SetYTitle("Correlation mean");
    hnseg8Tkl_proj = new TH1F("hnseg8Tkl_proj", "Correlation tkl-tkl wrt #Delta#eta, projection",NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
    hnseg8Tkl_proj->SetXTitle("Correlation #Delta#eta");
    hnseg8Tkl_proj->SetYTitle("Correlation mean");
    hnseg8ME_proj = new TH1D("hnseg8ME_proj", "MixedEvent Correlation wrt #Delta#eta, projection",NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
    hnseg8ME_proj->SetXTitle("Correlation #Delta#eta");
    hnseg8ME_proj->SetYTitle("Correlation mean");
    hnseg8TklME_proj = new TH1D("hnseg8TklME_proj", "MixedEvent Correlation tkl-tkl wrt #Delta#eta, projection",NbinsDeltaEtaTKL,MinDeltaEtaTKL,MaxDeltaEtaTKL);
    hnseg8TklME_proj->SetXTitle("Correlation #Delta#eta");
    hnseg8TklME_proj->SetYTitle("Correlation mean");
 
    Sijk = new TH1F("Sijk",
                                   "Sijk",
                                   NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Sijk->SetXTitle("Correlation #Delta#phi (rad)");
    Mijk = new TH1F("Mijk",
                                   "Mijk",
                                   NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    Mijk->SetXTitle("Correlation #Delta#phi (rad)");
    SoverMijk = new TH1F("SoverMijk",
                                   "SoverMijk",
                                   NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    SoverMijk->SetXTitle("Correlation #Delta#phi (rad)");
    SoverMik = new TH1F("SoverMik",
                                   "SoverMik",
                                   NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
    SoverMik->SetXTitle("Correlation #Delta#phi (rad)");
    
    Sij = new TH1F("Sij",
                                   "Sij",
                                   NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    Sij->SetXTitle("Correlation #Delta#phi (rad)");
    Mij = new TH1F("Mij",
                                   "Mij",
                                   NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    Mij->SetXTitle("Correlation #Delta#phi (rad)");
    SoverMij = new TH1F("SoverMij",
                                   "SoverMij",
                                   NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    SoverMij->SetXTitle("Correlation #Delta#phi (rad)");
    SoverMi = new TH1F("SoverMi",
                                   "SoverMi",
                                   NbinsDeltaPhiTKL,MinDeltaPhi,MaxDeltaPhi);
    SoverMi->SetXTitle("Correlation #Delta#phi (rad)");

    hnseg9 = new TH2F("hnseg9",
                      "Tracklet #eta wrt z_{vtx}",
                      100,-15,15,100,-2.5,2.5);
    hnseg9->SetXTitle("Z_{vtx} (cm)");
    hnseg9->SetYTitle("Tracklet #eta");
    hnseg10 = new TH2F("hnseg10",
                      "Number of tracklets wrt V0M percentile",
                      100,0,100,100,0,400);
    hnseg10->SetXTitle("V0M Percentile");
    hnseg10->SetYTitle("Tracklet count");
   // Double_t binsPt[7] = {0,2,3,4,6,8,12};
    hPtWrtMassInv[0] = new TH2F("hPtWrtMassInv_0",
                      "p_{T} wrt dimuon inv mass - All centralities",
                      NbinsDimuInvMass,MinInvMass,MaxInvMass,NbPtBins,PtBins);
    hPtWrtMassInv[0]->SetXTitle("Dimuon inv mass (GeV/c^{2})");
    hPtWrtMassInv[0]->SetYTitle("p_{T} (GeV/c)");
    hPtWrtMassInv[1] = new TH2F("hPtWrtMassInv_1",
                      "p_{T} wrt dimuon inv mass - Central",
                      NbinsDimuInvMass,MinInvMass,MaxInvMass,NbPtBins,PtBins);
    hPtWrtMassInv[1]->SetXTitle("Dimuon inv mass (GeV/c^{2})");
    hPtWrtMassInv[1]->SetYTitle("p_{T} (GeV/c)");
    hPtWrtMassInv[2] = new TH2F("hPtWrtMassInv_2",
                      "p_{T} wrt dimuon inv mass - Periph",
                      NbinsDimuInvMass,MinInvMass,MaxInvMass,NbPtBins,PtBins);
    hPtWrtMassInv[2]->SetXTitle("Dimuon inv mass (GeV/c^{2})");
    hPtWrtMassInv[2]->SetYTitle("p_{T} (GeV/c)");
    
    
    CentV0M = new TH1F("CentV0M",
                     "Distribution of CentV0M",
                     100,0,100);
    CentV0M->SetXTitle("CentV0M");
    CentV0M->SetYTitle("Count");
    CentTKL = new TH1F("CentTKL",
                     "Distribution of CentTKL",
                     100,0,100);
    CentTKL->SetXTitle("CentTKL");
    CentTKL->SetYTitle("Count");
    CentCL0 = new TH1F("CentCL0",
                     "Distribution of CentCL0",
                     100,0,100);
    CentCL0->SetXTitle("CentCL0");
    CentCL0->SetYTitle("Count");
    CentCL1 = new TH1F("CentCL1",
                     "Distribution of CentCL1",
                     100,0,100);
    CentCL1->SetXTitle("CentCL1");
    CentCL1->SetYTitle("Count");
    CentSPDTracklets = new TH1F("CentSPDTracklets",
                     "Distribution of CentSPDTracklets",
                     100,0,100);
    CentSPDTracklets->SetXTitle("CentSPDTracklets");
    CentSPDTracklets->SetYTitle("Count");
    CentSPDClusters = new TH1F("CentSPDClusters",
                     "Distribution of CentSPDClusters",
                     100,0,100);
    CentSPDClusters->SetXTitle("CentSPDClusters");
    CentSPDClusters->SetYTitle("Count");
    
    CentV0Mwrttkl = new TH2F("CentV0Mwrttkl",
                     "Distribution of CentV0Mwrttkl",
                     100,0,100, 100, 0, 100);
    CentV0Mwrttkl->SetXTitle("CentV0Mwrttkl");
    CentV0Mwrttkl->SetYTitle("NTracklets");
    CentTKLwrttkl = new TH2F("CentTKLwrttkl",
                     "Distribution of CentTKLwrttkl",
                     100,0,100, 100, 0, 100);
    CentTKLwrttkl->SetXTitle("CentTKLwrttkl");
    CentTKLwrttkl->SetYTitle("NTracklets");
    CentCL0wrttkl = new TH2F("CentCL0wrttkl",
                     "Distribution of CentCL0wrttkl",
                     100,0,100, 100, 0, 100);
    CentCL0wrttkl->SetXTitle("CentCL0wrttkl");
    CentCL0wrttkl->SetYTitle("NTracklets");
    CentCL1wrttkl = new TH2F("CentCL1wrttkl",
                     "Distribution of CentCL1wrttkl",
                     100,0,100, 100, 0, 100);
    CentCL1wrttkl->SetXTitle("CentCL1wrttkl");
    CentCL1wrttkl->SetYTitle("NTracklets");
    CentSPDTrackletswrttkl = new TH2F("CentSPDTrackletswrttkl",
                     "Distribution of CentSPDTrackletswrttkl",
                     100,0,100, 100, 0, 100);
    CentSPDTrackletswrttkl->SetXTitle("CentSPDTrackletswrttkl");
    CentSPDTrackletswrttkl->SetYTitle("NTracklets");
    CentSPDClusterswrttkl = new TH2F("CentSPDClusterswrttkl",
                     "Distribution of CentSPDClusterswrttkl",
                     100,0,100, 100, 0, 100);
    CentSPDClusterswrttkl->SetXTitle("CentSPDClusterswrttkl");
    CentSPDClusterswrttkl->SetYTitle("NTracklets");
    
    
    
    for(int ptbin=0;ptbin<NbPtBins;ptbin++){
        for(int i=0; i<NbBinsCent; i++){
               for(int k=0; k<NbinsInvMass; k++){
                   char hname[200];
                   sprintf(hname,"Yields %d %d %d ",i,k,ptbin);
                   YieldsPtBinned[i][k][ptbin] = new TH1F(hname,
                                     "YieldsPtBinned Correlation wrt #Delta#phi",
                                     NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
                   YieldsPtBinned[i][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                   sprintf(hname,"Yields1a %d %d %d ",i,k,ptbin);
                   YieldsPtBinned1a[i][k][ptbin] = new TH2F(hname,
                                     "Yields 1a Correlation wrt #Delta#phi and #Delta#eta",
                                     NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                   YieldsPtBinned1a[i][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                   YieldsPtBinned1a[i][k][ptbin]->SetYTitle("Correlation #Delta#eta");
                   for(int j=0; j<NbBinsZvtx; j++){
                       sprintf(hname,"CorrelationsPtBinned %d %d %d %d ",i,j,k,ptbin);
                       CorrelationsPtBinned[i][j][k][ptbin] = new TH2F(hname,
                                         "CorrelationPtBinned #Delta#eta wrt #Delta#phi",
                                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                       CorrelationsPtBinned[i][j][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                       CorrelationsPtBinned[i][j][k][ptbin]->SetYTitle("Correlation #Delta#eta");
                       sprintf(hname,"CorrelationsMEPtBinned %d %d %d %d ",i,j,k,ptbin);
                       CorrelationsMEPtBinned[i][j][k][ptbin] = new TH2I(hname,
                                         "MixedEvent CorrelationPtBinned #Delta#eta wrt #Delta#phi",
                                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                       CorrelationsMEPtBinned[i][j][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                       CorrelationsMEPtBinned[i][j][k][ptbin]->SetYTitle("Correlation #Delta#eta");
                       
                       sprintf(hname,"CorrelationsMEScaledPtBinned %d %d %d %d ",i,j,k,ptbin);
                       CorrelationsMEScaledPtBinned[i][j][k][ptbin] = new TH2F(hname,
                                         "MixedEvent CorrelationPtBinned #Delta#eta wrt #Delta#phi - Scaled",
                                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                       CorrelationsMEScaledPtBinned[i][j][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                       CorrelationsMEScaledPtBinned[i][j][k][ptbin]->SetYTitle("Correlation #Delta#eta");
                       if(i<3){
                           sprintf(hname,"CorrelationsPtBinnedBigClasses %d %d %d %d ",i,j,k,ptbin);
                           CorrelationsPtBinnedBigClasses[i][j][k][ptbin] = new TH2F(hname,
                                             "CorrelationPtBinnedBigClasses #Delta#eta wrt #Delta#phi",
                                             NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                           CorrelationsPtBinnedBigClasses[i][j][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                           CorrelationsPtBinnedBigClasses[i][j][k][ptbin]->SetYTitle("Correlation #Delta#eta");
                           
                           sprintf(hname,"CorrelationsMEPtBinnedBigClasses %d %d %d %d ",i,j,k,ptbin);
                           CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin] = new TH2I(hname,
                                             "MixedEvent CorrelationPtBinnedBigClasses #Delta#eta wrt #Delta#phi",
                                             NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                           CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                           CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin]->SetYTitle("Correlation #Delta#eta");
                           
                           sprintf(hname,"CorrelationsMEScaledPtBinnedBigClasses %d %d %d %d ",i,j,k,ptbin);
                           CorrelationsMEScaledPtBinnedBigClasses[i][j][k][ptbin] = new TH2F(hname,
                                             "MixedEvent CorrelationPtBinnedBigClasses #Delta#eta wrt #Delta#phi - Scaled",
                                             NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi,NbinsDeltaEta,MinDeltaEta,MaxDeltaEta);
                           CorrelationsMEScaledPtBinnedBigClasses[i][j][k][ptbin]->SetXTitle("Correlation #Delta#phi (rad)");
                           CorrelationsMEScaledPtBinnedBigClasses[i][j][k][ptbin]->SetYTitle("Correlation #Delta#eta");
                       }
                   }
               }
               
           }
        
        for (int j=0; j <NbinsInvMass; j++){
        char hname1[200];
        char hname2[200];
        sprintf(hname1,"Projected yield in Mass Bin %f GeV to %f GeV - Central - PtBinned %d",MassBins[j],MassBins[j+1],ptbin);
        sprintf(hname2,"Projected yield in Mass Bin, #m_{#mu#mu} #in [%.2f,%.2f] GeV/c^{2} - Central - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",MassBins[j],MassBins[j+1],PtBins[ptbin], PtBins[ptbin+1]);
        Yield_Central_MassBinPtBinned[j][ptbin] = new TH1F(hname1, hname2,NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
        }
        
        for (int j=0; j <NbinsInvMass; j++){
        char hname1[200];
        char hname2[200];
        sprintf(hname1,"Projected yield in Mass Bin %f GeV to %f GeV - Periph - PtBinned %d",MassBins[j],MassBins[j+1],ptbin);
        sprintf(hname2,"Projected yield in Mass Bin, #m_{#mu#mu} #in [%.2f,%.2f] GeV/c^{2} - Periph - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",MassBins[j],MassBins[j+1],PtBins[ptbin], PtBins[ptbin+1]);
        Yield_Periph_MassBinPtBinned[j][ptbin] = new TH1F(hname1, hname2,NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
        }
        
        for (int j=0; j <NbinsInvMass; j++){
        char hname1[200];
        char hname2[200];
        sprintf(hname1,"Projected yield in Mass Bin %f GeV to %f GeV - Difference - PtBinned %d",MassBins[j],MassBins[j+1],ptbin);
        sprintf(hname2,"Projected yield in Mass Bin, #m_{#mu#mu} #in [%.2f,%.2f] GeV/c^{2} - Difference - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",MassBins[j],MassBins[j+1],PtBins[ptbin], PtBins[ptbin+1]);
        Yield_Difference_MassBinPtBinned[j][ptbin] = new TH1F(hname1, hname2,NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
        }
        
        for(int i=0; i<NbBinsCent; i++){
                  for(int p=0; p<NbinsDeltaPhi; p++){
                          char hname[200];
                          sprintf(hname,"Yields PhiBinPtBinned %d %d %d",i,p,ptbin);
                          Yields_PhiBinPtBinned[i][p][ptbin] = new TH1F(hname,
                                            "Yields CorrelationPtBinned wrt mass",
                                            NbinsInvMass,MassBins);
                          Yields_PhiBinPtBinned[i][p][ptbin]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
                  }
           }
        
        for(int p=0; p<NbinsDeltaPhi; p++){
                char hname[200];
                char hname2[200];
                sprintf(hname,"Yields PhiBinPtBinned %d All C %d",p,ptbin);
                YieldWrtMass_allCPtBinned[p][ptbin] = new TH1F(hname,
                                  "Yields CorrelationPtBinned wrt dimuon mass, all C",
                                  NbinsInvMass,MassBins);
                YieldWrtMass_allCPtBinned[p][ptbin]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
            
            sprintf(hname,"Yields PhiBinPtBinned %d Periph %d",p,ptbin);
            sprintf(hname2,"Yields CorrelationPtBinned wrt dimuon mass, Periph, #Delta#phi #in [#frac{%d#pi}{6},#frac{%d#pi}{6}]",p-3, p-2);
            YieldWrtMass_PeriphPtBinned[p][ptbin] = new TH1F(hname,
                              hname2,
                              NbinsInvMass,MassBins);
            YieldWrtMass_PeriphPtBinned[p][ptbin]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
            
            sprintf(hname,"Yields PhiBinPtBinned %d Central %d",p,ptbin);
            sprintf(hname2,"Yields CorrelationPtBinned wrt dimuon mass, Central, #Delta#phi #in [#frac{%d#pi}{6},#frac{%d#pi}{6}]",p-3, p-2);
            YieldWrtMass_CentralPtBinned[p][ptbin] = new TH1F(hname,
                              hname2,
                             NbinsInvMass,MassBins);
            YieldWrtMass_CentralPtBinned[p][ptbin]->SetXTitle("Correlation dimuon Inv Mass (GeV/c^{2})");
        }
        
        char hname[200];
        char hname2[200];
        sprintf(hname,"Yields_Central_1PtBinned %d",ptbin);
        sprintf(hname2,"Yield of J/#psi-tkl in Central collisions wrt #Delta#phi - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        Yields_Central_1PtBinned[ptbin] = new TH1F(hname,
                         hname2,
                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
        Yields_Central_1PtBinned[ptbin]->SetXTitle("#Delta#phi (rad)");
        Yields_Central_1PtBinned[ptbin]->SetYTitle("Yield_{Central}");
        
        sprintf(hname,"Yields_Periph_1PtBinned %d",ptbin);
        sprintf(hname2,"Yield of J/#psi-tkl in Periph collisions wrt #Delta#phi - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        Yields_Periph_1PtBinned[ptbin] = new TH1F(hname,
                         hname2,
                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
        Yields_Periph_1PtBinned[ptbin]->SetXTitle("#Delta#phi (rad)");
        Yields_Periph_1PtBinned[ptbin]->SetYTitle("Yield_{Periph}");
        
        sprintf(hname,"Yields_Difference_1PtBinned %d",ptbin);
        sprintf(hname2,"Subtracted Yield of J/#psi-tkl wrt #Delta#phi - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        Yields_Difference_1PtBinned[ptbin] = new TH1F(hname,
                         hname2,
                         NbinsDeltaPhi,MinDeltaPhi,MaxDeltaPhi);
        Yields_Difference_1PtBinned[ptbin]->SetXTitle("#Delta#phi (rad)");
        Yields_Difference_1PtBinned[ptbin]->SetYTitle("Yield_{Subtracted}");
        
        sprintf(hname,"coefficients0PtBinned %d",ptbin);
        sprintf(hname2,"0^{th} Fourier coefficient - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        coefficients0PtBinned[ptbin] = new TH1F(hname,
                            hname2,
                            NbinsInvMass,MassBins);
           coefficients0PtBinned[ptbin]->SetXTitle("Mass of dimuon (GeV/c^{2})");
           coefficients0PtBinned[ptbin]->SetYTitle("0^{th} Fourier coefficient");
        
        sprintf(hname,"coefficients1PtBinned %d",ptbin);
        sprintf(hname2,"1^{st} Fourier coefficient - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        coefficients1PtBinned[ptbin] = new TH1F(hname,
                            hname2,
                            NbinsInvMass,MassBins);
           coefficients1PtBinned[ptbin]->SetXTitle("Mass of dimuon (GeV/c^{2})");
           coefficients1PtBinned[ptbin]->SetYTitle("1^{st} Fourier coefficient");
        
        sprintf(hname,"coefficients2PtBinned %d",ptbin);
        sprintf(hname2,"2^{nd} Fourier coefficient - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        coefficients2PtBinned[ptbin] = new TH1F(hname,
                            hname2,
                            NbinsInvMass,MassBins);
           coefficients2PtBinned[ptbin]->SetXTitle("Mass of dimuon (GeV/c^{2})");
           coefficients2PtBinned[ptbin]->SetYTitle("2^{nd} Fourier coefficient");
        
        sprintf(hname,"baselines0PtBinned %d",ptbin);
        sprintf(hname2,"Baseline_{Periph} - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        baselines0PtBinned[ptbin] = new TH1F(hname,
                         hname2,
                         NbinsInvMass,MassBins);
        baselines0PtBinned[ptbin]->SetXTitle("Mass of dimuon (GeV/c^{2})");
        baselines0PtBinned[ptbin]->SetYTitle("Baseline_{Periph}");
        
        sprintf(hname,"c2b0PtBinned %d",ptbin);
        sprintf(hname2,"2^{nd} Fourier coefficient + Baseline_{Periph} - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        c2b0PtBinned[ptbin] = new TH1F(hname,
                         hname2,
                         NbinsInvMass,MassBins);
        c2b0PtBinned[ptbin]->SetXTitle("Mass of dimuon (GeV/c^{2})");
        c2b0PtBinned[ptbin]->SetYTitle("2^{nd} Fourier coefficient + Baseline_{Periph}");
        
        sprintf(hname,"V2JPsiTklPtBinned %d",ptbin);
        sprintf(hname2,"V_{2,J/#psi-tkl} wrt dimuon mass - PtBinned, p_{T} #in [%.2f,%.2f] GeV/c",PtBins[ptbin], PtBins[ptbin+1]);
        V2JPsiTklPtBinned[ptbin] = new TH1F(hname,
                               hname2,
                               NbinsInvMass,MassBins);
              V2JPsiTklPtBinned[ptbin]->SetXTitle("Mass of dimuon (GeV/c^{2})");
              V2JPsiTklPtBinned[ptbin]->SetYTitle("V_{2,J/#psi-tkl}");
    }
    
    
    
    
    
// *************************
// Analyse                 *
// *************************
    
    
    int DimuCentralSeen = 0;
    int DimuPeriphSeen = 0;
    int DimuSeenMassCut = 0;
    int DimuSeenNoMassCut = 0;
    int EventRejected = 0;
    int EventNC = 0;
    int EventPileUpMult = 0;
    int EventPileUpVtx = 0;
    int RefTracklets = 0;
    int RefTrackletsCentral = 0;
    int RefTrackletsPeriph = 0;
    int cmul = 0;
    int barWidth = 50;
    int DimuonCounter[NbBinsCent][NbBinsZvtx][NbinsInvMass] = {0};
    int DimuonCounterZint[NbBinsCent][NbinsInvMass] = {0};
    int RefTklCounter[NbBinsCent][NbBinsZvtx] = {0};
    int RefTklCounterZint[NbBinsCent] = {0};
    double NormME[NbBinsCent][NbBinsZvtx][NbinsInvMass] = {0};
    double NormMETkl[NbBinsCent][NbBinsZvtx] = {0};
    
    //Nassoc
    int Nassoc_Central_TKL = 0;
    int Nassoc_Periph_TKL = 0;
    
    
    //For PtBinning
    int DimuSeenMassCutPtBinned[NbPtBins] = {0};
    int DimuonCounterPtBinned[NbBinsCent][NbBinsZvtx][NbinsInvMass][NbPtBins] = {0};
    int DimuonCounterZintPtBinned[NbBinsCent][NbinsInvMass][NbPtBins] = {0};
    int DimuonCounterPtBinnedBigClasses[3][NbBinsZvtx][NbinsInvMass][NbPtBins] = {0};
    int DimuonCounterZintPtBinnedBigClasses[3][NbinsInvMass][NbPtBins] = {0};
    double NormMEPtBinned[NbBinsCent][NbBinsZvtx][NbinsInvMass][NbPtBins] = {0};
    double NormMEPtBinnedBigClasses[3][NbBinsZvtx][NbinsInvMass][NbPtBins] = {0};
    
    int DimuC = 0;
    int DimuP = 0;
    int TklC = 0;
    int TklP = 0;
    int NormTklCentral = 0;
    int NormTklPeriph = 0;
    int countsigma = 0;
    
    MyEventLight* fEvent = 0;
    TClonesArray *fCorrelations = 0;
    TClonesArray *fTracklets = 0;
    TClonesArray *fDimuons = 0;
    CorrelationLight *correl = 0; //= new CorrelationLight(); //object must be created before
    TrackletLight *trac = 0;
    TrackletLight *tracklet1 = 0;
    TrackletLight *tracklet2 = 0;
    TrackletLight *trackletME = 0;
    DimuonLight *dimu = 0;
    TTree* theTree = NULL;
    
    
// ************************************
// Boucle sur les evenements, notés i *
// ************************************
    
    cout << "=== NOW EVENT STUDY ===" <<endl;
    cout << "Switching back to 1st TTree" <<endl;
    
    for(int tree_idx=0; tree_idx<numberOfPeriods; tree_idx++){
        
        //Define where the .root file with the TTree of the period can be found
        
   //     sprintf(fileInLoc,"~/../../Volumes/Sauvegarde /LegacySebAnalysepp/NewAnalysis/%s/muonGrid.root",arrayOfPeriods[tree_idx]);
        sprintf(fileInLoc,"~/../../Volumes/Transcend2/ppAnalysis/Scripts/NewAnalysis_AllEst/CMUL/%s/muonGrid.root",arrayOfPeriods[tree_idx]);
    //    sprintf(fileInLoc,"~/../../Volumes/Transcend2/ppAnalysis/Scripts/NewAnalysis_AllEst/CINT/%s_CINT_AllEst/muonGrid.root",arrayOfPeriods[tree_idx]);
    TFile fileIn(fileInLoc);
        
        char str [50];
        int GroupNum;
        
        std::vector <double> Pools[NbBinsCent][NbBinsZvtx]; //[12]
        std::vector <int> PoolsEventTracker[NbBinsCent][NbBinsZvtx];
        int PoolsSize[NbBinsCent][NbBinsZvtx] = {0}; //[12]
        std::vector <double> PoolsTkl[NbBinsCent][NbBinsZvtx];
        std::vector <int> PoolsTklEventTracker[NbBinsCent][NbBinsZvtx];
        int PoolsSizeTkl[NbBinsCent][NbBinsZvtx] = {0};
        
        std::vector <double> PoolsPtBinned[NbBinsCent][NbBinsZvtx][NbPtBins]; //[12]
        std::vector <int> PoolsEventTrackerPtBinned[NbBinsCent][NbBinsZvtx][NbPtBins];
        int PoolsSizePtBinned[NbBinsCent][NbBinsZvtx][NbPtBins] = {0}; //[12]

        sscanf (arrayOfPeriods[tree_idx],"Group%d_%s",&GroupNum,str);
    
    std::cout << "1" <<std::endl;
            fileIn.GetObject("MyMuonTree",theTree);
        std::cout << "2" <<std::endl;
            //setting the branch address
        std::cout << "3" <<std::endl;
      //  theTree->SetBranchAddress("event", &fEvent);
        TBranch *branch = theTree->GetBranch("event");
    //        auto bntrack = theTree->GetBranch("tracks");
            //auto branch  = theTree->GetBranch("mcparticles.fPx");
        std::cout << "3.5" <<std::endl;
            branch->SetAddress(&fEvent);
        std::cout << "4" <<std::endl;
        auto nevent = theTree->GetEntries();
        
        std::cout << "5" <<std::endl;
        fCorrelations = fEvent->fCorrelations;
        fTracklets = fEvent->fTracklets;
        fDimuons = fEvent->fDimuons;
    
    branch = theTree->GetBranch("event");
    branch->SetAddress(&fEvent);
    fCorrelations = fEvent->fCorrelations;
    fTracklets = fEvent->fTracklets;
    fDimuons = fEvent->fDimuons;
        int DimuMEcounter =0;
        int TklMEcounter =0;
        
        ofstream myfiletxt;
        myfiletxt.open("/tmp/memorycheckbin.txt");
    
        for (int i=0;i<nevent;i++) {
            if(i%100000 == 0){
                    std::cout << "[";
                    double portion = double(i)/nevent;
                    long pos = barWidth * portion;
                    for (int k = 0; k < barWidth; ++k) {
                        if (k < pos) std::cout << "=";
                        else if (k == pos) std::cout << ">";
                        else std::cout << " ";
                    }
                    std::cout << "] " << long(100 * portion) << "%     " << i << "/" << nevent << " Tree " << tree_idx+1 << "/" << numberOfPeriods;
                    std::cout.flush();
                std::cout << std::endl;
            }
            if(doTracklets && (i%4!=0)){
                continue;
            }
                theTree->GetEvent(i);

            if(fEvent->fNDimuons<1){
                continue;
            }
//            if(fEvent->fPassPhysicsSelection == 0){
//                continue;
//            }
            if(fEvent->fIsPileupFromSPDMultBins){
                EventPileUpMult++;
            }
//            if(fEvent->fNPileupVtx != 0){
//                EventPileUpVtx++;
//            }
            if(fEvent->fSPDVertexNC < 1){
                EventNC++;
            }
            if(fEvent->fIsPileupFromSPDMultBins || fEvent->fSPDVertexNC < 1){//} || fEvent->fIsPileupClustVsTkl == 1 || fEvent->fIsPileupReimplemented == 1 ){ //|| fEvent->fIsPileupClustVsTkl == 1  || fEvent->fNPileupVtx != 0 || fEvent->fNPileupVtx != 0)
                EventRejected++;
                continue; //PhysicsSelection
            }
          //  cout << "There are "<<fEvent->fNDimuons<< " dimuons, " << fEvent->fNTracklets<< " tracklets and "<< fEvent->fNCorrelations<< " correlations."<<endl;
            //For each event, compute the number of accepted tracklets see if the event is worth the analysis (also computes the number of tracklets that should be used when using an estimator calibrated by data-driven method)
            int NumberAcceptedEtaTracklets = 0;
            int NumberOfTrackletsForCentrality = 0;
            for (Int_t j=0; j<fEvent->fNTracklets; j++) {
                trac = (TrackletLight*)fTracklets->At(j);
                hsingletrac->Fill((trac->fPhi)+corrFactorDeltaPhi*(trac->fDPhi), trac->fEta);
               if((TMath::Abs(trac->fEta) < TklEtaCut) && (TMath::Abs(trac->fDPhi) < DPhiCut)){
                    NumberOfTrackletsForCentrality++;
                   //NumberOfTrackletsPassingEtaCut++;
               }
                    
                    if((trac->fEta < myEtaHigh(fEvent->fVertexZ)) && (trac->fEta > myEtaLow(fEvent->fVertexZ)) && (TMath::Abs(trac->fDPhi) < DPhiCut)){
                        NumberAcceptedEtaTracklets++;
                        hnseg9->Fill(fEvent->fVertexZ, trac->fEta);
                        hDPhi->Fill(trac->fDPhi);
                      //  NumberCloseEtaTracklets++;
                    }
                }
            
            for (Int_t j=0; j<fEvent->fNDimuons; j++) {
            dimu = (DimuonLight*)fDimuons->At(j);
            if ((TMath::Abs(fEvent->fVertexZ) < ZvtxCut) && (!fEvent->fIsFromVertexerZ || TMath::Abs(fEvent->fSPDVertexSigmaZ) < SigmaZvtxCut) && (dimu->fY < HighDimuYCut ) && (dimu->fY > LowDimuYCut) && (dimu->fCharge == 0) && (dimu->fPt > LowDimuPtCut) && (dimu->fPt < HighDimuPtCut)){
                int centint = 99;
                int ptint = 99;
                double zv = fEvent->fVertexZ;
                int zvint = floor(zv) + ZvtxCut;
               // centint = GetCentPM(NumberOfTrackletsForCentrality, fEvent->fCentralitySPDTracklets, fEvent->fCentralitySPDClusters, fEvent->fCentralityV0M, zvint, GroupNum);
                centint = GetCentCvetan(fEvent->fCentralityV0M);
                    ptint = GetPtBin(dimu->fPt);

                hnseg->Fill(dimu->fInvMass);
                  hPtWrtMassInv[0]->Fill(dimu->fInvMass, dimu->fPt);
                //  if(fEvent->fCentralitySPDTracklets<=CentSPDTrackletsCentral){
                  if(isCentral(centint)){
                      InvMass_Central->Fill(dimu->fInvMass);
                      hPtWrtMassInv[1]->Fill(dimu->fInvMass, dimu->fPt);
                  }
                //  if(fEvent->fCentralitySPDTracklets>=CentSPDTrackletsPeriph){
                  if(isPeripheral(centint)){
                      InvMass_Periph->Fill(dimu->fInvMass);
                      hPtWrtMassInv[2]->Fill(dimu->fInvMass, dimu->fPt);
                  }
                if(!isPeripheral(centint) && !isCentral(centint)){
                    InvMass_Neither->Fill(dimu->fInvMass);
                }

            }
            }
            
        
            if(NumberAcceptedEtaTracklets==0){
                continue;
            }
            
            if(KeepOnlyOne){
                if(valueOnlyOne!=NumberAcceptedEtaTracklets){
                    continue;
                }
            }
            
            if(AdditionalCutNtkl){
                if(NumberAcceptedEtaTracklets <= valueAdditionalCutNtkl){
                    continue;
                }
            }
            
            CentV0M->Fill(fEvent->fCentralityV0M);
//            CentTKL->Fill(fEvent->fCentralityTKL);
//            CentCL0->Fill(fEvent->fCentralityCL0);
//            CentCL1->Fill(fEvent->fCentralityCL1);
            CentSPDTracklets->Fill(fEvent->fCentralitySPDTracklets);
//            CentSPDClusters->Fill(fEvent->fCentralitySPDClusters);
            
            CentV0Mwrttkl->Fill(fEvent->fCentralityV0M, NumberAcceptedEtaTracklets);
//            CentTKLwrttkl->Fill(fEvent->fCentralityTKL, NumberCloseEtaTracklets);
//            CentCL0wrttkl->Fill(fEvent->fCentralityCL0, NumberCloseEtaTracklets);
//            CentCL1wrttkl->Fill(fEvent->fCentralityCL1, NumberCloseEtaTracklets);
            CentSPDTrackletswrttkl->Fill(fEvent->fCentralitySPDTracklets, NumberAcceptedEtaTracklets);
//            CentSPDClusterswrttkl->Fill(fEvent->fCentralitySPDClusters, NumberCloseEtaTracklets);

//            hnsegSigma->Fill(fEvent->fVertexNC, fEvent->fSPDVertexSigmaZ); //SIGMA
//
//            if(sqrt(fEvent->fSPDVertexSigmaZ)>0.25){
//                countsigma++;
//            }
            
          //  cout << "POUET"<<endl;
            
            bool HasAlreadyBeenAdded[NbPtBins] = {kFALSE};
            int DimuMECounterAdded[NbPtBins] = {0};
            
            for (Int_t j=0; j<fEvent->fNDimuons; j++) {
                dimu = (DimuonLight*)fDimuons->At(j);
                if ((TMath::Abs(fEvent->fVertexZ) < ZvtxCut) && (!fEvent->fIsFromVertexerZ || TMath::Abs(fEvent->fSPDVertexSigmaZ) < SigmaZvtxCut) && (dimu->fY < HighDimuYCut ) && (dimu->fY > LowDimuYCut) && (dimu->fCharge == 0) && (dimu->fPt > LowDimuPtCut) && (dimu->fPt < HighDimuPtCut)){
                    //(fEvent->fNPileupVtx == 0) &&
                    //  (TMath::Abs(fEvent->fVertexZ) < ZvtxCut) && (dimu->fEta < HighDimuEtaCut ) && (dimu->fEta > LowDimuEtaCut) && (dimu->fCharge == 0) && (dimu->fPt > LowDimuPtCut) && (dimu->fPt < HighDimuPtCut)
                    // fEvent->fNPileupVtx == 0) && (TMath::Abs(fEvent->fVertexZ) < ZvtxCut) && (dimu->fY < HighDimuYCut ) && (dimu->fY > LowDimuYCut) && (dimu->fCharge == 0) && (dimu->fPt > LowDimuPtCut) && (dimu->fPt < HighDimuPtCut)
                    DimuSeenNoMassCut++;
                    
                    int centint = 99;
                    int ptint = 99;
                    
                    if((dimu->fInvMass > MinInvMass) && (dimu->fInvMass < MaxInvMass)){
                        double mass = dimu->fInvMass;
                        int massint = GetMassBin(mass);
                       double cent = fEvent->fCentralitySPDTracklets;
                       double zv = fEvent->fVertexZ;
                       int zvint = floor(zv) + ZvtxCut;
                        
                        {
                            hnseg2->Fill(dimu->fPt );
                            hnseg3->Fill(dimu->fEta );
                            hnseg4->Fill(dimu->fY);
                            hnseg5->Fill(dimu->fPhi);
                        }
                        
                //      int centint = GetCent(cent);
                      //  centint = GetCentPM(NumberOfTrackletsForCentrality, fEvent->fCentralitySPDTracklets, fEvent->fCentralitySPDClusters, fEvent->fCentralityV0M, zvint, GroupNum);
                        centint = GetCentCvetan(fEvent->fCentralityV0M);
                        DimuonCounter[centint][zvint][massint]++;
                        if(PtBinned){
                            ptint = GetPtBin(dimu->fPt);
                            DimuonCounterPtBinned[centint][zvint][massint][ptint]++;
                            DimuonCounterPtBinnedBigClasses[2][zvint][massint][ptint]++;
                            if(isCentral(centint)){
                                  DimuonCounterPtBinnedBigClasses[0][zvint][massint][ptint]++;
                              }
                              if(isPeripheral(centint)){
                                 DimuonCounterPtBinnedBigClasses[1][zvint][massint][ptint]++;
                              }
                        }
                     //   if(cent <= CentSPDTrackletsCentral){
                        if(isCentral(centint)){
                            DimuC++;
                        }
                      //  else if(cent > CentSPDTrackletsPeriph){
                        if(isPeripheral(centint)){
                           DimuP++;
                        }
                    }
                    
                                
// Event mixing
                    
                    
                    if(doMixedEvents && PtBinned){
                        if((dimu->fInvMass > MinInvMass) && (dimu->fInvMass < MaxInvMass)){
                        double massME = dimu->fInvMass;
                            int massintME = GetMassBin(massME);
                     //  double centME = fEvent->fCentralitySPDTracklets;
                     //   int centintME = GetCent(centME);
                       double zvME = fEvent->fVertexZ;
                       int zvintME = floor(zvME) + ZvtxCut;
                            int ptintME = GetPtBin(dimu->fPt);
                       // int centintME = GetCentPM(NumberOfTrackletsForCentrality, fEvent->fCentralitySPDTracklets, fEvent->fCentralitySPDClusters, fEvent->fCentralityV0M, zvintME, GroupNum);
                            int centintME = GetCentCvetan(fEvent->fCentralityV0M);
                       int phintME = 0;
                       double phiME = dimu->fPhi;
                       if(phiME < -TMath::Pi()/2){
                        phiME += 2* TMath::Pi();
                          }
                          if(phiME > 1.5*TMath::Pi()){
                              phiME -= 2* TMath::Pi();
                          }
                       phintME = floor(phiME*6/TMath::Pi()) + 3;
                            DimuMEcounter++;

                        if(PoolsSizePtBinned[centintME][zvintME][ptintME]>=10){
                           for(int k=0; k< PoolsPtBinned[centintME][zvintME][ptintME].size(); k+=2){ //[phintME]
                               if(HasAlreadyBeenAdded[ptintME] && PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].at(k)==DimuMECounterAdded[ptintME]){
                                   continue;
                               }
                                       double correlMEPhi = PoolsPtBinned[centintME][zvintME][ptintME].at(k) - dimu->fPhi; //[phintME]
                         //      cout << "Pools[centintME][zvintME].at(k) " << k << " : " << Pools[centintME][zvintME].at(k) <<endl; //[phintME]
                               while(correlMEPhi<0 || correlMEPhi>TMath::Pi()){
                                  // cout<<"correlMEPhi "<<correlMEPhi<<endl;
                                       if(correlMEPhi <= -TMath::Pi()/2){
                                           correlMEPhi += 2* TMath::Pi();
                                       }
                                       if(correlMEPhi > 1.5*TMath::Pi()){
                                           correlMEPhi -= 2* TMath::Pi();
                                       }
                              if(correlMEPhi > 1.0*TMath::Pi() && correlMEPhi <= 1.5*TMath::Pi()){
                                   correlMEPhi = -correlMEPhi;
                               }
                                   if(correlMEPhi <= 0*TMath::Pi() && correlMEPhi > -0.5*TMath::Pi()){
                                       correlMEPhi = -correlMEPhi;
                                   }
                               }
                                       double correlMEEta = TMath::Abs(dimu->fY - PoolsPtBinned[centintME][zvintME][ptintME].at(k+1)); //FETAED //[phintME]
                                        if((correlMEEta>DeltaEtaDimuCut) && (correlMEEta<MaxDeltaEta)){
                                            CorrelationsMEPtBinned[centintME][zvintME][massintME][ptintME]->Fill(correlMEPhi,correlMEEta);
                                        }
                           }
                        }
                        if(PoolsSizePtBinned[centintME][zvintME][ptintME] == 100 && !HasAlreadyBeenAdded[ptintME]){
                            int valueDiscarded = PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].front();
                            while(PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].front()==valueDiscarded){
                                PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].erase(PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].begin(),PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].begin()+2);
                                PoolsPtBinned[centintME][zvintME][ptintME].erase(PoolsPtBinned[centintME][zvintME][ptintME].begin(),PoolsPtBinned[centintME][zvintME][ptintME].begin()+2);
                            }
                            
                            PoolsSizePtBinned[centintME][zvintME][ptintME] -= 1;
                            
                            bool hasBeenFilled = kFALSE;
                            
                            for (Int_t tk=0; tk<fEvent->fNTracklets; tk++) {
                               trackletME = (TrackletLight*)fTracklets->At(tk);
                                //FIXME ok
                               if(trackletME->fEta < myEtaHigh(zvME) && trackletME->fEta > myEtaLow(zvME) && (TMath::Abs(trackletME->fDPhi) < DPhiCut)){
                                   double trackletMEPhi = (Float_t)trackletME->fPhi + corrFactorDeltaPhi*(trackletME->fDPhi);
                                   double trackletMEEta = trackletME->fEta;
                                   PoolsPtBinned[centintME][zvintME][ptintME].push_back(trackletMEPhi);
                                   PoolsPtBinned[centintME][zvintME][ptintME].push_back(trackletMEEta);
                                   PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].push_back(DimuMEcounter);
                                   PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].push_back(DimuMEcounter);
                                   hasBeenFilled = kTRUE;
                               }

                           }
                            if(hasBeenFilled){
                                PoolsSizePtBinned[centintME][zvintME][ptintME] += 1;
                                HasAlreadyBeenAdded[ptintME] = kTRUE;
                                DimuMECounterAdded[ptintME] = DimuMEcounter;
                            }
                        }
                        else if(PoolsSizePtBinned[centintME][zvintME][ptintME] < 100 && !HasAlreadyBeenAdded[ptintME]){
                            bool hasBeenFilled = kFALSE;
                               for (Int_t tk=0; tk<fEvent->fNTracklets; tk++) {
                                   trackletME = (TrackletLight*)fTracklets->At(tk);
                                   //FIXME ok
                                   if(trackletME->fEta < myEtaHigh(zvME) && trackletME->fEta > myEtaLow(zvME) && (TMath::Abs(trackletME->fDPhi) < DPhiCut)){
                                       double trackletMEPhi = (Float_t)trackletME->fPhi + corrFactorDeltaPhi*(trackletME->fDPhi);
                                       double trackletMEEta = trackletME->fEta;
                                       PoolsPtBinned[centintME][zvintME][ptintME].push_back(trackletMEPhi);
                                       PoolsPtBinned[centintME][zvintME][ptintME].push_back(trackletMEEta);
                                       PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].push_back(DimuMEcounter);
                                       PoolsEventTrackerPtBinned[centintME][zvintME][ptintME].push_back(DimuMEcounter);
                                       hasBeenFilled = kTRUE;
                                   }

                               }
                               if(hasBeenFilled){
                                   PoolsSizePtBinned[centintME][zvintME][ptintME] += 1;
                                   HasAlreadyBeenAdded[ptintME] = kTRUE;
                                   DimuMECounterAdded[ptintME] = DimuMEcounter;
                               }
                           }
                            
                            
                            
                            
                            
                        }
                        
                    }
                                       
                                       
                }
            }
           
            
// ***********************************
// Boucle sur les correlations, notés j    *
// ***********************************
            for (Int_t j=0; j<fEvent->fNCorrelations; j++) {
                correl = (CorrelationLight*)fCorrelations->At(j);
                if ((TMath::Abs(fEvent->fVertexZ) < ZvtxCut) && (!fEvent->fIsFromVertexerZ || TMath::Abs(fEvent->fSPDVertexSigmaZ) < SigmaZvtxCut) && (correl->fDimuonY < HighDimuYCut ) && (correl->fDimuonY > LowDimuYCut) && (correl->fDimuonCharge == 0) && (correl->fDimuonPt > LowDimuPtCut) && (correl->fDimuonPt < HighDimuPtCut)  && (TMath::Abs(correl->fTrackletDPhi) < DPhiCut) && (correl->fTrackletEta < myEtaHigh(fEvent->fVertexZ)) && (correl->fTrackletEta > myEtaLow(fEvent->fVertexZ))){   //Cuts    (fEvent->fNPileupVtx == 0) && FIXME ok
                    Float_t DeltaPhi = correl->fDeltaPhi + (corrFactorDeltaPhi*(correl->fTrackletDPhi));
                    int tries = 0;
                    while((DeltaPhi<0 || DeltaPhi>TMath::Pi())&& tries<20){
                     //   cout<<"DeltaPhi is "<<DeltaPhi<<endl;
                    if(DeltaPhi <= -TMath::Pi()/2){
                      //  cout << "Below -pi/2: add 2pi"<<endl;
                        DeltaPhi += 2* TMath::Pi();
                    }
                    else if(DeltaPhi <= 0*TMath::Pi() && DeltaPhi > -TMath::Pi()/2){
                      //  cout << "-pi/2 to 0, rotate"<<endl;
                        DeltaPhi = -DeltaPhi;
                    }
                        else if(DeltaPhi > 1.0*TMath::Pi() && DeltaPhi <= 1.5*TMath::Pi()){
                       //     cout << "pi to 3pi/2, rotate"<<endl;
                        DeltaPhi = -DeltaPhi;
                    }
                    else if(DeltaPhi > 1.5*TMath::Pi()){
                       // cout << "above 3pi/2, -2pi"<<endl;
                        DeltaPhi -= 2* TMath::Pi();
                    }
                        tries++;
                        if(tries == 10){
                            cout << "TRUNCATURE HAS BEEN USED"<<endl;
                            DeltaPhi = trunc(DeltaPhi*10000)/10000.;
                        }
                    }
                 //   cout << "tries = "<<tries<<endl;

                    hnseg6->Fill(DeltaPhi);
                    hnseg7->Fill(correl->fTrackletEta-correl->fDimuonY);//FETAED
                    hnseg8->Fill(DeltaPhi, correl->fTrackletEta-correl->fDimuonY); //FETAED
                    
                    if((correl->fDimuonInvMass > MinInvMass) && (correl->fDimuonInvMass < MaxInvMass) && (correl->fTrackletEta-correl->fDimuonY>DeltaEtaDimuCut) && (correl->fTrackletEta-correl->fDimuonY<MaxDeltaEta) ){//FETAED
                    double mass = correl->fDimuonInvMass;
                        int massint = GetMassBin(mass);
                   // double cent = fEvent->fCentralitySPDTracklets;
                  //  int centint = GetCent(cent);
                    double zv = fEvent->fVertexZ;
                    int zvint = floor(zv) + ZvtxCut;
                   // int centint = GetCentPM(NumberOfTrackletsForCentrality, fEvent->fCentralitySPDTracklets, fEvent->fCentralitySPDClusters, fEvent->fCentralityV0M, zvint, GroupNum);
                        int centint = GetCentCvetan(fEvent->fCentralityV0M);
                    //    cout << centint << " " << zvint << " " << massint <<endl;
                    Correlations[centint][zvint][massint]->Fill(DeltaPhi,correl->fTrackletEta-correl->fDimuonY); //FETAED
                        if(PtBinned){
                            int ptint = GetPtBin(correl->fDimuonPt);
                            CorrelationsPtBinned[centint][zvint][massint][ptint]->Fill(DeltaPhi,correl->fTrackletEta-correl->fDimuonY); //FETAED
                        }
                        
                    }
                }
            }
            
            if((TMath::Abs(fEvent->fVertexZ) < ZvtxCut) && (!fEvent->fIsFromVertexerZ || TMath::Abs(fEvent->fSPDVertexSigmaZ) < SigmaZvtxCut)){
                hnseg10->Fill(fEvent->fCentralityV0M, NumberAcceptedEtaTracklets); //fEvent->fNTracklets
            }
            
            
            
    }
        myfiletxt.close();
    }
    cout << "Reading of events finished" <<endl;
  //  cout << "CorrelationsME, CorrelationsTklME are being set to 1" <<endl;
    

    cout << "Calculation of DimuonCOunterZint[Cent][mass]" <<endl;
    for(int i=0; i<NbBinsCent; i++){
              for(int k=0; k<NbinsInvMass; k++){
                  for(int j=0; j<NbBinsZvtx; j++){
                      DimuonCounterZint[i][k] += DimuonCounter[i][j][k];
                      if((j>8) && (j<14)){
                      cout << "DimuonCounter[cent i =" << i << "][zvint j = " << j << "][mass k =" << k << "]: " << DimuonCounter[i][j][k]<<endl;
                      }
                  }
              }
    }
    if(PtBinned){
        for(int i=0; i<3; i++){ //BigClasses
                          for(int k=0; k<NbinsInvMass; k++){
                              for(int j=0; j<NbBinsZvtx; j++){
                                  for(int ptbin=0; ptbin<NbPtBins; ptbin++){
                                      DimuonCounterZintPtBinnedBigClasses[i][k][ptbin] += DimuonCounterPtBinnedBigClasses[i][j][k][ptbin];
        //                              if((j>8) && (j<14)){
        //                              cout << "DimuonCounter[cent i =" << i << "][zvint j = " << j << "][mass k =" << k << "]: " << DimuonCounter[i][j][k]<<endl;
        //                              }
                                  }
                              }
                          }
                }
        //Summation sur Centralité BigClasses
        for(int ptbin=0; ptbin<NbPtBins; ptbin++){
            for(int j=0; j<NbBinsZvtx; j++){
                for(int k=0; k<NbinsInvMass; k++){
                    
                   for(int i=CentralLowBound; i<CentralHighBound+1; i++){ //CentPeriph
                       CorrelationsPtBinnedBigClasses[0][j][k][ptbin]->Add(CorrelationsPtBinned[i][j][k][ptbin]);
                       CorrelationsMEPtBinnedBigClasses[0][j][k][ptbin]->Add(CorrelationsMEPtBinned[i][j][k][ptbin]);
                   }
                    
                    CorrelationsPtBinnedBigClasses[0][j][k][ptbin]->Sumw2();
                    CorrelationsMEPtBinnedBigClasses[0][j][k][ptbin]->Sumw2();
                    
                    for(int i=PeripheralLowBound; i<PeripheralHighBound+1; i++){
                        CorrelationsPtBinnedBigClasses[1][j][k][ptbin]->Add(CorrelationsPtBinned[i][j][k][ptbin]);
                        CorrelationsMEPtBinnedBigClasses[1][j][k][ptbin]->Add(CorrelationsMEPtBinned[i][j][k][ptbin]);
                    }
                    
                    CorrelationsPtBinnedBigClasses[1][j][k][ptbin]->Sumw2();
                    CorrelationsMEPtBinnedBigClasses[1][j][k][ptbin]->Sumw2();
                    
                    for(int i=0; i<NbBinsCent; i++){
                        CorrelationsPtBinnedBigClasses[2][j][k][ptbin]->Add(CorrelationsPtBinned[i][j][k][ptbin]);
                        CorrelationsMEPtBinnedBigClasses[2][j][k][ptbin]->Add(CorrelationsMEPtBinned[i][j][k][ptbin]);
                    }
                    
                    CorrelationsPtBinnedBigClasses[2][j][k][ptbin]->Sumw2();
                    CorrelationsMEPtBinnedBigClasses[2][j][k][ptbin]->Sumw2();
                    
                    
                }
            }
        }
        
        if(doMixedEvents){
            for(int ptbin=0; ptbin<NbPtBins; ptbin++){
            
                cout << "Normalisation of CorrelationsME based on max of Eta projected applied" <<endl;
                    for(int i=0; i<3; i++){ //BigClasses
                        for(int j=0; j<NbBinsZvtx; j++){
                            for(int k=0; k<NbinsInvMass; k++){
                                float maxME = 0;
                                int maxindexME = 0;
                                hnseg8ME_proj_tampon = (TH1F*)(CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin]->ProjectionY("hnseg8ME_proj",1,NbinsDeltaPhi,"e")); //Change underflow
                                for(int l=1; l<hnseg8ME_proj->GetNbinsX()+1; l++){
                                    float valueME = hnseg8ME_proj_tampon->GetBinContent(l);
                                   // cout << "valueME " << valueME <<endl;
                                    if (valueME > maxME){
                                        maxME = valueME;
                                        maxindexME = l;
                                    }
                                }
                                if(maxME>0){
                                    NormMEPtBinnedBigClasses[i][j][k][ptbin] = 1/maxME;
                                }
                            }
                        }
                    }
            
                
                    //Scaling all CorrelationsME CVETANSCALE (line2987 was 0)
                    
                    for(int i=0; i<3; i++){
                           for(int j=0; j<NbBinsZvtx; j++){
                               for(int k=0; k<NbinsInvMass; k++){
                                   if(NormMEPtBinnedBigClasses[i][j][k][ptbin]>0){
                                       for(int binx=1; binx<(1+CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin]->GetNbinsX()); binx++){
                                           for(int biny=1; biny<(1+CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin]->GetNbinsY()); biny++){
                                               double error = CorrelationsMEScaledPtBinnedBigClasses[i][j][k][ptbin]->GetBinError(binx,biny);
                                               CorrelationsMEScaledPtBinnedBigClasses[i][j][k][ptbin]->SetBinContent(binx,biny, (CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin]->GetBinContent(binx,biny))*NormMEPtBinnedBigClasses[i][j][k][ptbin]);
                                               CorrelationsMEScaledPtBinnedBigClasses[i][j][k][ptbin]->SetBinError(binx,biny,error*NormMEPtBinnedBigClasses[i][j][k][ptbin]); //(CorrelationsMEPtBinnedBigClasses[i][j][k][ptbin]->GetBinError(binx,biny))*NormMEPtBinnedBigClasses[i][j][k][ptbin]); MEchange
                                           }
                                       }
                                   }
                                   
                               }
                           }
                    }
            }
                
        }
        
        
        
        
    }

    
    
    if(PtBinned){
        
        if(kTRUE){ // Summation 1c
            char hname2[200];
            for(int ptbin=0; ptbin<NbPtBins; ptbin++){
                    for(int i=0; i<3; i++){//BigClasses
                        for(int k=0; k<NbinsInvMass; k++){
                    
                        for(int binEta=1; binEta<=NbinsDeltaEta; binEta++){
                            
                        for(int binPhi=1; binPhi<=NbinsDeltaPhi; binPhi++){
                            sprintf(hname2,"htemp_pt%i_cent%i_mass%i_eta%i_phi%i",ptbin,i,k,binEta-1,binPhi-1);
                            TH1D *htemp = new TH1D(hname2,hname2,NbBinsZvtx,0,1);
                            sprintf(hname2,"htempMix_pt%i_cent%i_mass%i_eta%i_phi%i",ptbin,i,k,binEta-1,binPhi-1);
                            TH1D *htempMix = new TH1D(hname2,hname2,NbBinsZvtx,0,1);
                          for(int binZvtx=1; binZvtx<=NbBinsZvtx; binZvtx++) {
                htemp->SetBinContent(binZvtx,CorrelationsPtBinnedBigClasses[i][binZvtx-1][k][ptbin]->GetBinContent(binPhi,binEta));
                htempMix->SetBinContent(binZvtx,DimuonCounterPtBinnedBigClasses[i][binZvtx-1][k][ptbin]*CorrelationsMEScaledPtBinnedBigClasses[i][binZvtx-1][k][ptbin]->GetBinContent(binPhi,binEta));
                              
                          }
                          Double_t err1;
                          Double_t cont1 = htemp->IntegralAndError(1,NbBinsZvtx,err1);
                          Double_t err2;
                          Double_t cont2 = htempMix->IntegralAndError(1,NbBinsZvtx,err2);
                            cout << hname2<<endl;
                          if (htempMix->Integral() > 0) {
                        TF1 *fPoissonLL = new TF1("fPoissonLL",fitPLL,0,1,1);
                        hStatic3 = htempMix;
                        htemp->Fit(fPoissonLL,"LOS");
                              if(ptbin==0 && i==0 && k > NbinsInvMass-3){
                                  cout << "Potential FISHY 1 !!"<<endl;
                                  TCanvas* cAAAfish1 = new TCanvas;
                                  cAAAfish1->Divide(1,2);
                                  cAAAfish1->cd(1);
                                  htemp->DrawCopy();
                                  cAAAfish1->cd(2);
                                  htempMix->DrawCopy();
                                  cAAAfish1->Update();
                              }
        //                chi2 += fPoissonLL->GetChisquare()/fPoissonLL->GetNDF();
        //                nchi2++;
                YieldsPtBinned1a[i][k][ptbin]->SetBinContent(binPhi,binEta,fPoissonLL->GetParameter(0));
                YieldsPtBinned1a[i][k][ptbin]->SetBinError(binPhi,binEta,fPoissonLL->GetParError(0));
                        delete fPoissonLL;
                          }
                          delete htemp;
                          delete htempMix;
                        }
                        }
                            
                            if(ptbin==0 && i==0 && k > NbinsInvMass-3){
                                cout << "DRAWTIME"<<endl;
                                TCanvas* cAAAfish2 = new TCanvas;
                                cAAAfish2->cd();
                                YieldsPtBinned1a[i][k][ptbin]->DrawCopy("colz");
                                cAAAfish2->Update();
                            }
                        
                            char hname1[200];
                        bool kWasFail1a = kFALSE;
                                        for (int binPhi = 1; binPhi<=NbinsDeltaPhi; binPhi++){
                                            sprintf(hname1,"ME_proj_tampon_pt%i_cent%i_mass%i_phi%",ptbin,i,k,binPhi-1);
                                                            ME_proj_tampon = (TH1F*)(YieldsPtBinned1a[i][k][ptbin]->ProjectionY(hname1,binPhi,binPhi));//TestErr
                                        //                    for(int binEta=1; binEta<=(NbinsDeltaEtaTKL/2); binEta++){
                                        //                        ME_proj_Tkl_tamponEta->SetBinContent(NbinsDeltaEtaTKL-binEta+1, (ME_proj_Tkl_tamponEta->GetBinContent(binEta)+ME_proj_Tkl_tamponEta->GetBinContent(NbinsDeltaEtaTKL-binEta+1))/2);
                                        //                        ME_proj_Tkl_tamponEta->SetBinError(NbinsDeltaEtaTKL-binEta+1, sqrt(pow(ME_proj_Tkl_tamponEta->GetBinContent(binEta),2)+pow(ME_proj_Tkl_tamponEta->GetBinContent(NbinsDeltaEtaTKL-binEta+1),2))/2);
                                        //                    }
                                        //                    for(int binEta=(NbinsDeltaEtaTKL/2)+1; binEta<=NbinsDeltaEtaTKL; binEta++){
                                        //                        ME_proj_Tkl_tamponEta->SetBinContent(NbinsDeltaEtaTKL-binEta+1, 0);
                                        //                    }
                                            if(ME_proj_tampon->GetEntries() == 0){
                                                continue;
                                            }
                                                            TF1 *fLinear = new TF1("fLinear", fitLin, DeltaEtaDimuCut ,MaxDeltaEta,2);
                                                              Double_t paramslin[] = {0.01,1.};
                                                              fLinear->SetParameters(paramslin);
                                                               TVirtualFitter::Fitter(ME_proj_tampon)->SetMaxIterations(100000);
                                                               TVirtualFitter::Fitter(ME_proj_tampon)->SetPrecision();

                                                                       fLinear->SetParName(0,"b");
                                                                       fLinear->SetParName(1,"a");
                                                                  //  fLinear->SetParLimits(0,0,20);

                                                                      TFitResultPtr res = ME_proj_tampon->Fit("fLinear","SQ0","ep");
                                                            if(res->CovMatrixStatus() == 3){
                                                                cout<<"Linear fit for summation binPhi= "<<binPhi<< "is a success"<<endl;
                                                            }
                                                            else{
                                                                cout<<"Linear fit for summation binPhi= "<<binPhi<< "is a failure"<<endl;
                                                                kWasFail1a = kTRUE;
                                                            }
                                                                Double_t parlin[2];
                                                                fLinear->GetParameters(parlin);
                                                                fLinear->SetParameters(parlin);
                                                                cout<<"a = "<<parlin[1]<< " b = "<<parlin[0]<<endl;
                                            cout << "Chi2/ndf = " << fLinear->GetChisquare()/fLinear->GetNDF()<<endl;
                                            
                                            if(fLinear->GetChisquare()/fLinear->GetNDF()>2){
                                                cout << "FISHY !!"<<endl;
                                                TCanvas* cAAAfish = new TCanvas;
                                                cAAAfish->cd();
                                                ME_proj_tampon->DrawCopy();
                                                cAAAfish->Update();
                                            }
                                            TFile *f = new TFile(FitFileName,"UPDATE");
                                            ME_proj_tampon->Write();
                                            f->Close();
//CVETANSCALE ?
                                                                double integral = fLinear->Integral(DeltaEtaDimuCut,MaxDeltaEta);
                                                                double integralerror = fLinear->IntegralError(DeltaEtaDimuCut,MaxDeltaEta, fLinear->GetParameters(), res->GetCovarianceMatrix().GetMatrixArray());

                                                                cout << "Integral = " <<integral << " error = " <<integralerror<<endl;
                                                                YieldsPtBinned[i][k][ptbin]->SetBinContent(binPhi, integral);
                                                                YieldsPtBinned[i][k][ptbin]->SetBinError(binPhi, integralerror);

                                                        }

                                        if(kWasFail1a){
                                            cout << "Summation Method1c did not execute properly"<<endl;
                                        }

                                        TFile *f = new TFile(FitFileName,"UPDATE");
                                        YieldsPtBinned[i][k][ptbin]->Write();
                                        YieldsPtBinned1a[i][k][ptbin]->Write();
                                        f->Close();
                        
                    }
            }
                    
                    cout << "Reorganisaton of Yields[Cent][mass] to create Yields_PhiBin wrt inv mass and calculation them bor all C, periph, central" <<endl;
                    //TRUC A PLOT
                    
                        for(int i=0; i<3; i++){//BigClasses
                            for(int k=0; k<NbinsInvMass; k++){
                                Yield_tampon->Reset();
                                Yield_tampon->Add(YieldsPtBinned[i][k][ptbin]);
                              //  Yield_tampon->Scale(DimuonCounterZintPtBinnedBigClasses[i][k][ptbin]);
                                for(int deltaphibin=0; deltaphibin<NbinsDeltaPhi; deltaphibin++){
                                    Yields_PhiBinPtBinned[i][deltaphibin][ptbin]->SetBinContent(k+1,Yield_tampon->GetBinContent(deltaphibin+1));
                                    Yields_PhiBinPtBinned[i][deltaphibin][ptbin]->SetBinError(k+1,Yield_tampon->GetBinError(deltaphibin+1));
                                }
                            }
                        }
                    
                }
        }
        
     }
   
  
    cout << "DimuSeenMassCut " << DimuSeenMassCut <<endl;
    cout << "Calculation of Yields in all C, central, periph and difference" <<endl;
    
    if(PtBinned){
        for(int ptbin=0; ptbin<NbPtBins; ptbin++){
            for(int k=0; k<NbinsInvMass; k++){
 //                  DimuCnt=0;//BigClasses
//                   for(int i=CentralLowBound; i<CentralHighBound+1; i++){ //CentPeriph
//                       if(DimuonCounterZintPtBinned[i][k][ptbin]>0){
                           Yield_tampon->Reset();
                           Yield_tampon->Add(YieldsPtBinned[0][k][ptbin]);
//                           Yield_tampon->Scale(DimuonCounterZintPtBinned[i][k][ptbin]);
                           Yield_Central_MassBinPtBinned[k][ptbin]->Add(Yield_tampon);
//                           DimuCnt += DimuonCounterZintPtBinned[i][k][ptbin];
//                       }
//                   }
//                cout << "ptbin: " << ptbin << " Massbin: " << k << " DimuCnt Central: " <<DimuCnt<<endl;
//                if(DimuCnt>0){
//                    Yield_Central_MassBinPtBinned[k][ptbin]->Scale(1./DimuCnt);
//                }
//                   DimuCnt=0;
//                   for(int i=PeripheralLowBound; i<PeripheralHighBound+1; i++){
//                       if(DimuonCounterZintPtBinned[i][k][ptbin]>0){
                           Yield_tampon->Reset();
                           Yield_tampon->Add(YieldsPtBinned[1][k][ptbin]);
//                           Yield_tampon->Scale(DimuonCounterZintPtBinned[i][k][ptbin]);
                           Yield_Periph_MassBinPtBinned[k][ptbin]->Add(Yield_tampon);
//                           DimuCnt += DimuonCounterZintPtBinned[i][k][ptbin];
//
//                           if(k == NbinsInvMass-1 && ptbin ==0){
//                               cout << "Looking at YieldsPtBinned ptbin 0 massbin last centralbin "<< i <<endl;
//                               for(int index = 0; index<Yield_tampon->GetNbinsX(); index++){
//                                   cout<< "Content bin " << index << " is " << Yield_tampon->GetBinContent(index+1)<<endl;
//                                   cout<< "Error bin " << index << " is " <<Yield_tampon->GetBinError(index+1)<<endl;
//                               }
//                           }
//
//                           if(k == NbinsInvMass-1 && ptbin ==0){
//                                              cout << "Looking at progressive Periph Yield ptbin 0 massbin last"<<endl;
//                                              for(int index = 0; index<Yield_Periph_MassBinPtBinned[k][ptbin]->GetNbinsX(); index++){
//                                                  cout<< "Content bin " << index << " is " << Yield_Periph_MassBinPtBinned[k][ptbin]->GetBinContent(index+1)<<endl;
//                                                  cout<< "Error bin " << index << " is " << Yield_Periph_MassBinPtBinned[k][ptbin]->GetBinError(index+1)<<endl;
//                                              }
//                                          }
//
//                       }
//                   }
//                cout << "ptbin: " << ptbin << " Massbin: " << k << " DimuCnt Periph: " <<DimuCnt<<endl;
//                if(k == NbinsInvMass-1 && ptbin ==0){
//                    cout << "Looking at Periph Yield ptbin 0 massbin last"<<endl;
//                    for(int index = 0; index<Yield_Periph_MassBinPtBinned[k][ptbin]->GetNbinsX(); index++){
//                        cout<< "Content bin " << index << " is " << Yield_Periph_MassBinPtBinned[k][ptbin]->GetBinContent(index+1)<<endl;
//                        cout<< "Error bin " << index << " is " << Yield_Periph_MassBinPtBinned[k][ptbin]->GetBinError(index+1)<<endl;
//                    }
//                }
//                if(DimuCnt>0){
//                   Yield_Periph_MassBinPtBinned[k][ptbin]->Scale(1./DimuCnt);
//                }
                   Yield_Difference_MassBinPtBinned[k][ptbin]->Add(Yield_Central_MassBinPtBinned[k][ptbin],Yield_Periph_MassBinPtBinned[k][ptbin],1,-1);
                         
               }
        }
    }
    
    TCanvas*cTestRebinning = new TCanvas();
    cTestRebinning->Divide(1,2);
    cTestRebinning->cd(1);
    hnseg->DrawCopy();
    
    int RebinningFactor = NbinsDimuInvMass/NbinsInvMass;
    
        cout << "Rebinning of InvMass plots to prepare for division" <<endl;
    TH1I *InvMass_TotRebinned = dynamic_cast<TH1I*>(hnseg->Rebin(NbinsInvMass,"InvMass_TotRebinned",MassBins));
    TH1I *InvMass_CentralRebinned = dynamic_cast<TH1I*>(InvMass_Central->Rebin(NbinsInvMass,"InvMass_CentralRebinned",MassBins));
    TH1I *InvMass_PeriphRebinned = dynamic_cast<TH1I*>(InvMass_Periph->Rebin(NbinsInvMass,"InvMass_PeriphRebinned",MassBins));

    
    
    if(PtBinned){
        for(int ptbin=0;ptbin<NbPtBins;ptbin++){
            for(int bin_idx=1; bin_idx<=hPtWrtMassInvSlicedRebinned[0][ptbin]->GetNbinsX();bin_idx++){
                hPtWrtMassInvSlicedRebinned[0][ptbin]->SetBinError(bin_idx,0);
                hPtWrtMassInvSlicedRebinned[1][ptbin]->SetBinError(bin_idx,0);
                hPtWrtMassInvSlicedRebinned[2][ptbin]->SetBinError(bin_idx,0);
            }
            
            int DimuCntMass[NbinsInvMass] = {0};
            for(int p=0; p<NbinsDeltaPhi; p++){
                for(int k = 0; k<NbinsInvMass;k++){
                    DimuCntMass[k] = 0;
                }
//                for(int i=0; i<NbBinsCent; i++){//BigClasses
                    YieldWrtMass_tampon->Reset();
                    YieldWrtMass_tampon->Add(Yields_PhiBinPtBinned[2][p][ptbin]);
                    YieldWrtMass_allCPtBinned[p][ptbin]->Add(YieldWrtMass_tampon);
                    for(int k = 0; k<NbinsInvMass;k++){
                        DimuCntMass[k] += DimuonCounterZintPtBinnedBigClasses[2][k][ptbin];
                    }
//                }
                
//                for(int k = 0; k<NbinsInvMass;k++){
//                    YieldWrtMass_allCPtBinned[p][ptbin]->SetBinContent(k+1, YieldWrtMass_allCPtBinned[p][ptbin]->GetBinContent(k+1)/DimuCntMass[k]);
//                    YieldWrtMass_allCPtBinned[p][ptbin]->SetBinError(k+1, YieldWrtMass_allCPtBinned[p][ptbin]->GetBinError(k+1)/DimuCntMass[k]);
//                }
                
            //    YieldWrtMass_allCPtBinned[p][ptbin]->Divide(hPtWrtMassInvSlicedRebinned[0][ptbin]);
                
            }
            
             for(int p=0; p<NbinsDeltaPhi; p++){
                 for(int k = 0; k<NbinsInvMass;k++){
                     DimuCntMass[k] = 0;
                 }
//                for(int i=CentralLowBound; i<CentralHighBound+1; i++){ //CentPeriph BigClasses
                    YieldWrtMass_tampon->Reset();
                    YieldWrtMass_tampon->Add(Yields_PhiBinPtBinned[0][p][ptbin]);
                    YieldWrtMass_CentralPtBinned[p][ptbin]->Add(YieldWrtMass_tampon);
//
                        for(int k = 0; k<NbinsInvMass;k++){
                            DimuCntMass[k] += DimuonCounterZintPtBinnedBigClasses[0][k][ptbin];
                        }
//                }
                 
//                 for(int k = 0; k<NbinsInvMass;k++){
//                     YieldWrtMass_CentralPtBinned[p][ptbin]->SetBinContent(k+1, YieldWrtMass_CentralPtBinned[p][ptbin]->GetBinContent(k+1)/DimuCntMass[k]);
//                     YieldWrtMass_CentralPtBinned[p][ptbin]->SetBinError(k+1, YieldWrtMass_CentralPtBinned[p][ptbin]->GetBinError(k+1)/DimuCntMass[k]);
//                 }
                 
                 
              //  YieldWrtMass_CentralPtBinned[p][ptbin]->Divide(hPtWrtMassInvSlicedRebinned[1][ptbin]);
                 
            }
            
            for(int p=0; p<NbinsDeltaPhi; p++){
                for(int k = 0; k<NbinsInvMass;k++){
                    DimuCntMass[k] = 0;
                }
//                for(int i=PeripheralLowBound; i<PeripheralHighBound+1; i++){ //BigClasses
                    YieldWrtMass_tampon->Reset();
                    YieldWrtMass_tampon->Add(Yields_PhiBinPtBinned[1][p][ptbin]);
                    YieldWrtMass_PeriphPtBinned[p][ptbin]->Add(YieldWrtMass_tampon);
                    for(int k = 0; k<NbinsInvMass;k++){
                        DimuCntMass[k] += DimuonCounterZintPtBinnedBigClasses[1][k][ptbin];
                    }
  //              }
//                for(int k = 0; k<NbinsInvMass;k++){
//                    YieldWrtMass_PeriphPtBinned[p][ptbin]->SetBinContent(k+1, YieldWrtMass_PeriphPtBinned[p][ptbin]->GetBinContent(k+1)/DimuCntMass[k]);
//                    YieldWrtMass_PeriphPtBinned[p][ptbin]->SetBinError(k+1, YieldWrtMass_PeriphPtBinned[p][ptbin]->GetBinError(k+1)/DimuCntMass[k]);
//                }
                
                //YieldWrtMass_PeriphPtBinned[p][ptbin]->Divide(hPtWrtMassInvSlicedRebinned[2][ptbin]);
                
            }
        }
        
    }
        
    
    std::cout << "==================== Analysis Terminated - PLOTTING TIME ====================" << std::endl;
    //}
    
    
// ***********************************
// Tracer les graphes sur les canevas *
// ***********************************


    TCanvas*c1=new TCanvas();
    //Canvas avec plein de petits plots généraux
    c1->SetTitle("General Plots");
    c1->Divide(4,2);
    c1->cd(1);
    hnseg->SetFillColor(18);
    hnseg->Draw();
    c1->cd(2);
    hnseg2->SetFillColor(18);
    hnseg2->Draw();
    c1->cd(3);
    hnseg3->SetFillColor(18);
    hnseg3->Draw();
    c1->cd(4);
    hnseg4->SetFillColor(18);
    hnseg4->Draw();
    c1->cd(5);
    hnseg5->SetFillColor(18);
    hnseg5->Draw();
    c1->cd(6);
    hnseg6->SetFillColor(18);
    hnseg6->Draw();
    c1->cd(7);
    hnseg7->SetFillColor(18);
    hnseg7->Draw();
    c1->cd(8);
    hnseg8->Draw("colz");
    sprintf(CanvasName,"%s/GeneralPlots.pdf",FolderName);
    c1->SaveAs(CanvasName);
    
    TCanvas*cJav=new TCanvas();
    hsingletrac->Draw("colz");
    sprintf(CanvasName,"%s/TklEtaPhi.pdf",FolderName);
    cJav->SaveAs(CanvasName);
    
    TCanvas*cDPhi=new TCanvas();
    cDPhi->SetTitle("DPhi distribution");
    //Plot DPhi distribution
    hDPhi->SetFillColor(18);
    hDPhi->Draw();
    sprintf(CanvasName,"%s/DPhiDist.pdf",FolderName);
    cDPhi->SaveAs(CanvasName);

    TF1 *fEtaLow = new TF1("fEtaLow","myEtaLow(x)",-10,10);
    TF1 *fEtaHigh = new TF1("fEtaHigh","myEtaHigh(x)",-10,10);
    
    TCanvas*c2=new TCanvas();
    c2->SetTitle("TklEtaWrtZvtx");
    //Plot TrkEtaWrtZvtx
    hnseg9->Draw("colz");
    fEtaLow->Draw("same");
    fEtaHigh->Draw("same");
    sprintf(CanvasName,"%s/TklEtaZvtx.pdf",FolderName);
    c2->SaveAs(CanvasName);
    
    TCanvas*c3=new TCanvas();
    c3->SetTitle("Ntrk wrt V0M");
    //Plut Ntrk wrt V0M
    hnseg10->Draw("colz");
    
    TCanvas* c8zoom = new TCanvas;
    c8zoom->SetTitle("Yields in Mass Bins CtrlPeriphDiff zoom");
    //Correlation yield DeltaEta wrt DeltaPhi TH2 -> Central, Periph, Difference [MassBins]
    c8zoom->Divide(3,3);
    for(int i=1; i<=3; i++){
        c8zoom->cd(i);
        Yield_Central_MassBin[i+1]->Draw("E");
        c8zoom->cd(3+i);
        Yield_Periph_MassBin[i+1]->Draw("E");
        c8zoom->cd(2*3+i);
        Yield_Difference_MassBin[i+1]->Draw("E");
    }
    
    TFile *FitFile;
    FitFile = new TFile(FitFileName,"RECREATE");
    
    TCanvas* c8 = new TCanvas;
    // ROOT THAT Yield in bins of mass
    c8->SetTitle("Yields in Mass Bins CtrlPeriphDiff");
    //Correlation yield DeltaEta wrt DeltaPhi TH2 -> Central, Periph, Difference [MassBins]
    c8->Divide(NbinsInvMass,3);
    for(int i=1; i<=NbinsInvMass; i++){
        c8->cd(i);
        Yield_Central_MassBin[i-1]->Draw("E");
        Yield_Central_MassBin[i-1]->Write();
        c8->cd(NbinsInvMass+i);
        Yield_Periph_MassBin[i-1]->Draw("E");
        Yield_Periph_MassBin[i-1]->Write();
        c8->cd(2*NbinsInvMass+i);
        Yield_Difference_MassBin[i-1]->Draw("E");
        Yield_Difference_MassBin[i-1]->Write();
    }
    
    for(int i=1; i<NbinsInvMass+1; i++){
        baselines0->SetBinContent(i, (Yield_Periph_MassBin[i-1]->GetBinContent(BinZeroLeft)));//Pichange +Yield_Periph_MassBin[i-1]->GetBinContent(BinZeroLeft+1))/2);
        baselines0->SetBinError(i, sqrt(pow(Yield_Periph_MassBin[i-1]->GetBinError(BinZeroLeft),2)/2+pow(Yield_Periph_MassBin[i-1]->GetBinError(BinZeroLeft),2)/2));
    }
    
    if(PtBinned){
        for(int ptbin=0;ptbin<NbPtBins;ptbin++){
            for(int i=1; i<NbinsInvMass+1; i++){
                baselines0PtBinned[ptbin]->SetBinContent(i, (Yield_Periph_MassBinPtBinned[i-1][ptbin]->GetBinContent(BinZeroLeft)+Yield_Periph_MassBinPtBinned[i-1][ptbin]->GetBinContent(BinZeroLeft))/2);
                baselines0PtBinned[ptbin]->SetBinError(i, sqrt(pow(Yield_Periph_MassBinPtBinned[i-1][ptbin]->GetBinError(BinZeroLeft),2)/2+pow(Yield_Periph_MassBinPtBinned[i-1][ptbin]->GetBinError(BinZeroLeft),2)/2));
            }
        }
        
        for(int ptbin=0;ptbin<NbPtBins;ptbin++){
            for(int i=1; i<=NbinsInvMass; i++){
                Yield_Central_MassBinPtBinned[i-1][ptbin]->Write();
                Yield_Periph_MassBinPtBinned[i-1][ptbin]->Write();
                Yield_Difference_MassBinPtBinned[i-1][ptbin]->Write();

            }
        }
        
        for(int ptbin=0;ptbin<NbPtBins;ptbin++){
            for(int i=0; i<NbinsDeltaPhi; i++){
                YieldWrtMass_CentralPtBinned[i][ptbin]->Write();
                YieldWrtMass_PeriphPtBinned[i][ptbin]->Write();

            }
        }
    }
    
    TCanvas* c10 = new TCanvas;
    // ROOT THAT yields phi bins
    c10->SetTitle("Yields wrt mass Central Periph");
    // Yields wrt mass -> Periph, Central [Phi bins]
    c10->Divide(NbinsDeltaPhi/2,4);
    for(int i=1; i<NbinsDeltaPhi+1; i++){
        c10->cd(i);
        YieldWrtMass_Central[i-1]->Draw("E");
        YieldWrtMass_Central[i-1]->Write();
        c10->cd(NbinsDeltaPhi+i);
        YieldWrtMass_Periph[i-1]->Draw("E");
        YieldWrtMass_Periph[i-1]->Write();
    }
    
    FitFile->Close();
    
    TCanvas*c12=new TCanvas();
    c12->SetTitle("Centrality Estimators");
    // Estimators
    c12->Divide(3,2);
    c12->cd(1);
    CentV0M->SetFillColor(18);
    CentV0M->Draw();
    c12->cd(2);
    CentTKL->SetFillColor(18);
    CentTKL->Draw();
    c12->cd(3);
    CentCL0->SetFillColor(18);
    CentCL0->Draw();
    c12->cd(4);
    CentCL1->SetFillColor(18);
    CentCL1->Draw();
    c12->cd(5);
    CentSPDTracklets->SetFillColor(18);
    CentSPDTracklets->Draw();
    c12->cd(6);
    CentSPDClusters->SetFillColor(18);
    CentSPDClusters->Draw();
    
    TCanvas*c13=new TCanvas();
    c13->SetTitle("Estimators wrt tkl");
    //Estimators
    c13->Divide(3,2);
    c13->cd(1);
    CentV0Mwrttkl->SetFillColor(18);
    CentV0Mwrttkl->Draw();
    c13->cd(2);
    CentTKLwrttkl->SetFillColor(18);
    CentTKLwrttkl->Draw();
    c13->cd(3);
    CentCL0wrttkl->SetFillColor(18);
    CentCL0wrttkl->Draw();
    c13->cd(4);
    CentCL1wrttkl->SetFillColor(18);
    CentCL1wrttkl->Draw();
    c13->cd(5);
    CentSPDTrackletswrttkl->SetFillColor(18);
    CentSPDTrackletswrttkl->Draw();
    c13->cd(6);
    CentSPDClusterswrttkl->SetFillColor(18);
    CentSPDClusterswrttkl->Draw();
    
    
    TCanvas*c14=new TCanvas();
    c14->SetTitle("Yield wrt Phi difference fits");
    //Yield difference wrt Phi [Mass bins] fits
    c14->Divide(5,2);
    for(int i=1; i<=NbinsInvMass; i++){
        cout << "Fit of the difference to extract Fourier wrt mass pt integrated. Mass bin: " << i <<endl;
        c14->cd(i);
    // Ici on fit YieldsWrtDeltaPhiMassBin_DifferenceProj
        TH1F *histo = Yield_Difference_MassBin[i-1];
      // create a TF1 with the range from 0 to 3 and 6 parameters
      TF1 *fitFcnV2_2 = new TF1("fitFcnV2_2",FourierV2,0,1.0*TMath::Pi(),3); //Pichange
      fitFcnV2_2->SetNpx(500);
      fitFcnV2_2->SetLineWidth(4);
      fitFcnV2_2->SetLineColor(kMagenta);
      // first try without starting values for the parameters
      // This defaults to 1 for each param.
      // this results in an ok fit for the polynomial function
      // however the non-linear part (lorenzian) does not
      // respond well.
       Double_t params[3] = {1,1,1};
      fitFcnV2_2->SetParameters(params);
       TVirtualFitter::Fitter(histo)->SetMaxIterations(10000);
       TVirtualFitter::Fitter(histo)->SetPrecision();
    //  histo->Fit("fitFcn","0");
      // second try: set start values for some parameters
       
       fitFcnV2_2->SetParName(0,"a0");
       fitFcnV2_2->SetParName(1,"a1");
       fitFcnV2_2->SetParName(2,"a2");
       
      TFitResultPtr res = histo->Fit("fitFcnV2_2","SBMERI+","ep");
        gStyle->SetOptStat("n");
        gStyle->SetOptFit(1011);
        if(!doTracklets){
        TPaveStats *st = (TPaveStats*)histo->FindObject("stats");
        st->SetX1NDC(0.75); //new x start position
        st->SetY1NDC(0.75); //new x end position
        }
      // improve the pictu
    //   std::cout << "integral error: " << integralerror << std::endl;
        Double_t par[3];
        fitFcnV2_2->GetParameters(par);
        coefficients0->SetBinContent(i, par[0]);
        coefficients1->SetBinContent(i, par[1]);
        coefficients2->SetBinContent(i, par[2]);
        coefficients0->SetBinError(i,fitFcnV2_2->GetParError(0));
        coefficients1->SetBinError(i,fitFcnV2_2->GetParError(1));
        coefficients2->SetBinError(i,fitFcnV2_2->GetParError(2));
      fitFcnV2_2->Draw("same");
        cout << "STATUS COV : " << res->CovMatrixStatus() <<endl;
      // draw the legend
      TLegend *legend=new TLegend(0.12,0.80,0.60,0.90);
      legend->SetFillColorAlpha(kWhite, 0.);
      legend->SetBorderSize(0);
        legend->SetTextFont(42);
      legend->SetTextSize(0.03);
        Char_t message[80];
        sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitFcnV2_2->GetChisquare(),fitFcnV2_2->GetNDF());
        legend->AddEntry(fitFcnV2_2,message);
        if(res->CovMatrixStatus() == 3){
          //  sprintf(message,"The fit is a success");
        }
        else{
            sprintf(message,"The fit is a failure");
            legend->AddEntry(fitFcnV2_2,message);
        }
      legend->AddEntry(histo,"Data","lpe");
      legend->Draw();
    //FIXME Ajouter d'autres méthodes ok
    }
    
    if(PtBinned){
        for(int ptbin=0;ptbin<NbPtBins;ptbin++){
            for(int i=1; i<=NbinsInvMass; i++){
                
                cout << "Fit of the difference to extract Fourier wrt mass. Mass bin: " << i <<" Pt bin: "<<ptbin <<endl;
            // Ici on fit YieldsWrtDeltaPhiMassBin_DifferenceProj
                TH1F *histo = Yield_Difference_MassBinPtBinned[i-1][ptbin];
              // create a TF1 with the range from 0 to 3 and 6 parameters
              TF1 *fitFcnV2_2PtBinned = new TF1("fitFcnV2_2PtBinned",FourierV2,0,1.0*TMath::Pi(),3);//Pichange
              fitFcnV2_2PtBinned->SetNpx(500);
              fitFcnV2_2PtBinned->SetLineWidth(4);
              fitFcnV2_2PtBinned->SetLineColor(kMagenta);
              // first try without starting values for the parameters
              // This defaults to 1 for each param.
              // this results in an ok fit for the polynomial function
              // however the non-linear part (lorenzian) does not
              // respond well.
               Double_t params[3] = {1,1,1};
              fitFcnV2_2PtBinned->SetParameters(params);
               TVirtualFitter::Fitter(histo)->SetMaxIterations(10000);
               TVirtualFitter::Fitter(histo)->SetPrecision();
            //  histo->Fit("fitFcn","0");
              // second try: set start values for some parameters
               
               fitFcnV2_2PtBinned->SetParName(0,"a0");
               fitFcnV2_2PtBinned->SetParName(1,"a1");
               fitFcnV2_2PtBinned->SetParName(2,"a2");
               
              TFitResultPtr res = histo->Fit("fitFcnV2_2PtBinned","SBMERI+","ep");
                gStyle->SetOptStat("n");
                gStyle->SetOptFit(1011);
                if(!doTracklets){
                TPaveStats *st = (TPaveStats*)histo->FindObject("stats");
                st->SetX1NDC(0.8); //new x start position
                st->SetY1NDC(0.8); //new x end position
                }
              // improve the pictu
            //   std::cout << "integral error: " << integralerror << std::endl;
                Double_t par[3];
                fitFcnV2_2PtBinned->GetParameters(par);
                coefficients0PtBinned[ptbin]->SetBinContent(i, par[0]);
                coefficients1PtBinned[ptbin]->SetBinContent(i, par[1]);
                coefficients2PtBinned[ptbin]->SetBinContent(i, par[2]);
                coefficients0PtBinned[ptbin]->SetBinError(i,fitFcnV2_2PtBinned->GetParError(0));
                coefficients1PtBinned[ptbin]->SetBinError(i,fitFcnV2_2PtBinned->GetParError(1));
                coefficients2PtBinned[ptbin]->SetBinError(i,fitFcnV2_2PtBinned->GetParError(2));
            //  fitFcnV2_2->Draw("same");
                cout << "STATUS COV : " << res->CovMatrixStatus() <<endl;
              // draw the legend
//              TLegend *legend=new TLegend(0.12,0.80,0.60,0.90);
//              legend->SetTextFont(61);
//              legend->SetTextSize(0.03);
//                Char_t message[80];
//                sprintf(message,"Global Fit : #chi^{2}/NDF = %.2f / %d",fitFcnV2_2->GetChisquare(),fitFcnV2_2->GetNDF());
//                legend->AddEntry(fitFcnV2_2,message);
//                if(res->CovMatrixStatus() == 3){
//                    sprintf(message,"The fit is a success");
//                }
//                else{
//                    sprintf(message,"The fit is a failure");
//                }
//                legend->AddEntry(fitFcnV2_2,message);
//              legend->AddEntry(histo,"Data","lpe");
//              legend->Draw();
            //FIXME Ajouter d'autres méthodes ok
            }
        }
    }
    
    
    TFile *f = new TFile(FitFileName,"UPDATE");
    
    TCanvas*c15=new TCanvas();
    // ROOT THAT Fouriers
    c15->SetTitle("Fourier coefs of yields wrt Phi");
    // Fourier coefficients of Yields wrt Phi [Mass bins] -> Extraction method 2
    c15->Divide(2,2);
    c15->cd(1);
    coefficients0->Draw();
    coefficients0->Write();
    c15->cd(2);
    baselines0->Draw();
    baselines0->Write();
    c15->cd(3);
    coefficients1->Draw();
    coefficients1->Write();
    c15->cd(4);
    coefficients2->Draw();
    coefficients2->Write();
    
    sprintf(CanvasName,"%s/Fouriers.pdf",FolderName);
    c15->SaveAs(CanvasName);
            

    // Plot du V2 JPsi Tracklet
    TCanvas*c16=new TCanvas();
    // ROOT THAT Extraction 2 plot
    c16->SetTitle("Extraction method 2");
    //V2_2 Jpsi-tkl wrt Mass fit (Extraction method 2)
    c2b0->Add(coefficients0, baselines0);
    V2JPsiTkl->Divide(coefficients2, c2b0);
    
    V2JPsiTkl->Draw();
    V2JPsiTkl->Write();
    
    hnseg->Write();
    InvMass_Central->Write();
    InvMass_Periph->Write();
    InvMass_Neither->Write();
    
    for(int j=0;j<3;j++){
        for (int i=0;i<NbPtBins;i++){
            hPtWrtMassInvSliced[j][i]->Write();
        }
    }
    
    // Plot du V2 JPsi Tracklet PtBinned
    TCanvas*c16PtBinned=new TCanvas();
    // ROOT THAT Extraction 2 plot
    c16PtBinned->SetTitle("Extraction method 2 PtBinned");
    c16PtBinned->Divide(2,3);
    //V2_2 Jpsi-tkl wrt Mass fit (Extraction method 2)
    for(int ptbin=0;ptbin<NbPtBins;ptbin++){
        c16PtBinned->cd(ptbin+1);
        c2b0PtBinned[ptbin]->Add(coefficients0PtBinned[ptbin], baselines0PtBinned[ptbin]);
        V2JPsiTklPtBinned[ptbin]->Divide(coefficients2PtBinned[ptbin], c2b0PtBinned[ptbin]);
        
//        if(ptbin==2){
//            //Set valeurs Cvetan
//          //  V2JPsiTklPtBinned[ptbin]->SetBinContent(2,0.00277);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(2,0.000997);
//           // V2JPsiTklPtBinned[ptbin]->SetBinContent(3,0.00144);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(3,0.0011);
//           // V2JPsiTklPtBinned[ptbin]->SetBinContent(4,0.00183);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(4,0.00127);
//            //V2JPsiTklPtBinned[ptbin]->SetBinContent(5,0.00312);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(5,0.00138);
//            //V2JPsiTklPtBinned[ptbin]->SetBinContent(6,0.00373);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(6,0.000957);
//            //V2JPsiTklPtBinned[ptbin]->SetBinContent(7,0.000234);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(7,0.00196);
//            //V2JPsiTklPtBinned[ptbin]->SetBinContent(8,0.00156);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(8,0.00247);
//            //V2JPsiTklPtBinned[ptbin]->SetBinContent(9,-0.000989);
//            V2JPsiTklPtBinned[ptbin]->SetBinError(9,0.00293);
//        }
        
        V2JPsiTklPtBinned[ptbin]->Draw();
    }

    
    for(int ptbin=0;ptbin<NbPtBins;ptbin++){
        coefficients0PtBinned[ptbin]->Write();
        coefficients1PtBinned[ptbin]->Write();
        coefficients2PtBinned[ptbin]->Write();
        baselines0PtBinned[ptbin]->Write();
        V2JPsiTklPtBinned[ptbin]->Write();
    }
    
    f->Close();
    
    TFile *outputFile;
     outputFile = new TFile("/Volumes/SEBUSB/InvMass9.root","RECREATE");
     hnseg->Write();
    InvMass_Central->Write();
    InvMass_Periph->Write();
    InvMass_Neither->Write();
    for(int j=0;j<3;j++){
        for (int i=0;i<NbPtBins;i++){
            hPtWrtMassInvSliced[j][i]->Write();
        }
    }
    
//    for(int i=0; i<12; i++){
//        InvMass_PhiBin[i]->Write();
//    }
     outputFile->Close();
    
    
    
//    TCanvas *cinvmassphibin = new TCanvas("cinvmassphibin","Fitting PhiBins",10,10,700,500);
//    cinvmassphibin->Divide(4,3);
    TCanvas *cinvmass = new TCanvas("cinvmass","Fitting All Phi - Different Centralities",10,10,700,500);
    cinvmass->SetTitle("Inv Mass Fits");
    // ROOT THAT Inv mass fits
    cinvmass->Divide(1,4);
    
    TFitResultPtr res;
    TFitResultPtr rescent;
    TFitResultPtr resperiph;
    TFitResult resu;
    char histoname[50];
    //Values OK for 20M evts
//    double startvalues[16] = {800, 3.096916, 0.07, 0.9, 10, 2, 15, 15,   10,     0.1,          50000,   0.01,        1,         1,   1,   1};
//    double lowvalues[16] =   {1,    3.05,    0.03, 0.7, 1,  1, 5,  0.01, 0.0001,  0.000000001,  10,    0.0000000001, 0.000001, -10, -10, -10};
//    double highvalues[16] = {100000, 3.15,   0.1,  1.1, 30, 5, 25, 50,   100000, 10,           100000, 10,     10,        10,  10,  10};
    
    double startvalues[16] = {1100, 3.096916, 0.07,  0.9, 10, 2, 15, 20,    90,      0.01,     80000,      1.4,      1,         1,   1,   1};
    double lowvalues[16] =   {10,    3.05,     0.03, 0.7, 1,  1, 5,  0.001,  1,    0.00001,   1000,        0.0001, 0.00000001, -10, -10, -10};
    double highvalues[16] = {100000, 3.15,     1,  1.1, 30, 5, 25, 10000, 1000000, 100,      10000000,  100,    10,        10,  10,  10};
    
//    for(int index=0; index<16; index++){
//        if(index==0 || index==7 || index==8 || index==10){
//            startvalues[index]*=nevent1/20000000;
//            lowvalues[index]*=nevent1/20000000;
//            highvalues[index]*=nevent1/20000000;
//        }
//    }
    //Values adapted to the number of events considered
    sprintf(histoname,"hnseg");
    res = FittingAllInvMassBin(histoname, cinvmass, 0);
    double par[16];
    
    // Parameter setting (either from mass fit either start values)
    for(int i=0; i <16; i++){
            par[i] = res->Parameter(i);
            if(i>11){
                par[i] = 1;
            }
    }
    
    //Definition and design of function on V2 plot
    c16->cd();
    TF1 *fitV2_2 = new TF1("fitV2_2",FourierV2_WrtInvMass,MinInvMass,MaxInvMass,16);
    fitV2_2->SetNpx(500);
    fitV2_2->SetLineWidth(4);
    fitV2_2->SetLineColor(kMagenta);
    fitV2_2->SetParameters(par);
    TF1 *fitJustV2_2 = new TF1("fitJustV2_2",FourierV2_WrtInvMass,MinInvMass,MaxInvMass,4);
    fitJustV2_2->SetNpx(500);
    fitJustV2_2->SetLineWidth(0);
    fitJustV2_2->SetLineColor(kWhite);
    TF1 *backFcnV2_2 = new TF1("backFcnV2_2",BackFcnV2Poly,2.1,5.1,16);
    backFcnV2_2->SetLineColor(kRed);
    TF1 *signalFcnJPsiV2_2 = new TF1("signalFcnJPsiV2_2",SignalFcnJPsiV2,2.1,5.1,16);
    signalFcnJPsiV2_2->SetLineColor(kBlue);
    signalFcnJPsiV2_2->SetNpx(500);
    
       TVirtualFitter::Fitter(V2JPsiTkl)->SetMaxIterations(10000);
       TVirtualFitter::Fitter(V2JPsiTkl)->SetPrecision();
        for(int i=0; i<=11; i++){
            fitV2_2->FixParameter(i,par[i]);
        }
       
       fitV2_2->SetParName(0,"Norm_{J/#psi}");
          fitV2_2->SetParName(1,"M_{J/#psi}");
          fitV2_2->SetParName(2,"Sigma_{J/#psi}");
          fitV2_2->SetParName(3,"a_{1}");
          fitV2_2->SetParName(4,"n_{1}");
          fitV2_2->SetParName(5,"a_{2}");
          fitV2_2->SetParName(6,"n_{2}");
          fitV2_2->SetParName(7,"Norm_{#Psi(2S)}");
          fitV2_2->SetParName(8,"Norm_{TailLowM}");
          fitV2_2->SetParName(9,"Exp_{TailLowM}");
          fitV2_2->SetParName(10,"Norm_{TailHighM}");
          fitV2_2->SetParName(11,"Exp_{TailHighM}");
           fitV2_2->SetParName(12,"V2_2 J/#psi-tkl");
           fitV2_2->SetParName(13,"V2_2 Bkg M2");
       fitV2_2->SetParName(14,"V2_2 Bkg M1");
       fitV2_2->SetParName(15,"V2_2 Bkg M0");
    
    //PSI2sTEST
    fitV2_2->FixParameter(7,0);
    
    fitJustV2_2->SetParName(0,"V2_2 J/#psi-tkl");
        fitJustV2_2->SetParName(1,"V2_2 Bkg M2");
    fitJustV2_2->SetParName(2,"V2_2 Bkg M1");
    fitJustV2_2->SetParName(3,"V2_2 Bkg M0");
    
       double minParams[16];
       double parErrors[16];
    
    //Fit of V2
    cout << "Fit of V2 wrt mass"<<endl;
        res = V2JPsiTkl->Fit("fitV2_2","SBMERI+","ep");
    gStyle->SetOptFit(1011);
           gStyle->SetOptStat("n");
    if(!doTracklets){
    TPaveStats *st = (TPaveStats*)V2JPsiTkl->FindObject("stats");
           st->SetX1NDC(0.75); //new x start position
           st->SetY1NDC(0.75); //new x end position
    }
        Double_t para[16];
        fitV2_2->GetParameters(para);
        backFcnV2_2->SetParameters(para);
        signalFcnJPsiV2_2->SetParameters(para);
    
    
    fitV2_2->Draw("same");
  //  signalFcnJPsiV2_2->Draw("same");
    backFcnV2_2->Draw("same");
      // draw the legend
      TLegend *legend=new TLegend(0.12,0.80,0.60,0.90);
      legend->SetFillColorAlpha(kWhite, 0.);
      legend->SetBorderSize(0);
        legend->SetTextFont(42);
        legend->SetTextSize(0.03);
    Char_t message[80];
    sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitV2_2->GetChisquare(),fitV2_2->GetNDF());
     legend->AddEntry(fitV2_2,message);

        if(res->CovMatrixStatus() == 3){
                   // sprintf(message,"The fit is a success");
             }
             else{
                    sprintf(message,"The fit is a failure");
                 legend->AddEntry(fitV2_2,message);
             }
        // legend->AddEntry(signalFcnJPsiV2_2,"JPsi signal");
    legend->AddEntry(backFcnV2_2,"Background");
      legend->AddEntry(V2JPsiTkl,"Data","lpe");
      legend->Draw();
    
    
    
    sprintf(CanvasName,"%s/V2_Extraction2.pdf",FolderName);
    c16->SaveAs(CanvasName);
    
    
    //Ext 2 Pt Binned
    
     TCanvas *cinvmassPtBinned = new TCanvas("cinvmassPtBinned","Fitting All Phi - Different Centralities - Pt Binned",10,10,700,500);
        cinvmassPtBinned->SetTitle("Inv Mass Fits - Pt Binned");
        // ROOT THAT Inv mass fits
        cinvmassPtBinned->Divide(NbPtBins,3);
        
    if(PtBinned){
        for(int ptbin=0;ptbin<NbPtBins;ptbin++){
    
        sprintf(histoname,Form("bin%d_0",ptbin+1)); // ZZZZZZZ
        res = FittingAllInvMassBin(histoname, cinvmassPtBinned, ptbin);
        double par[16];
        
        // Parameter setting (either from mass fit either start values)
        for(int i=0; i <16; i++){
                par[i] = res->Parameter(i);
                if(i>11){
                    par[i] = 1;
                }
        }
        
        //Definition and design of function on V2 plot
        c16PtBinned->cd(ptbin+1);
        TF1 *fitV2_2PtBinned = new TF1("fitV2_2PtBinned",FourierV2_WrtInvMass,1.5,4.5,16);//Trying this
        fitV2_2PtBinned->SetNpx(500);
        fitV2_2PtBinned->SetLineWidth(4);
        fitV2_2PtBinned->SetLineColor(kMagenta);
        fitV2_2PtBinned->SetParameters(par);
        TF1 *backFcnV2_2 = new TF1("backFcnV2_2",BackFcnV2Poly,2.1,5.1,16);
        backFcnV2_2->SetLineColor(kRed);
        TF1 *signalFcnJPsiV2_2 = new TF1("signalFcnJPsiV2_2",SignalFcnJPsiV2,2.1,5.1,16);
        signalFcnJPsiV2_2->SetLineColor(kBlue);
        signalFcnJPsiV2_2->SetNpx(500);
        
           TVirtualFitter::Fitter(V2JPsiTklPtBinned[ptbin])->SetMaxIterations(10000);
           TVirtualFitter::Fitter(V2JPsiTklPtBinned[ptbin])->SetPrecision();
            for(int i=0; i<=11; i++){
                fitV2_2PtBinned->FixParameter(i,par[i]);
            }
           
           fitV2_2PtBinned->SetParName(0,"Norm_{J/#psi}");
              fitV2_2PtBinned->SetParName(1,"M_{J/#psi}");
              fitV2_2PtBinned->SetParName(2,"Sigma_{J/#psi}");
              fitV2_2PtBinned->SetParName(3,"a_{1}");
              fitV2_2PtBinned->SetParName(4,"n_{1}");
              fitV2_2PtBinned->SetParName(5,"a_{2}");
              fitV2_2PtBinned->SetParName(6,"n_{2}");
              fitV2_2PtBinned->SetParName(7,"Norm_{#Psi(2S)}");
              fitV2_2PtBinned->SetParName(8,"Norm_{TailLowM}");
              fitV2_2PtBinned->SetParName(9,"Exp_{TailLowM}");
              fitV2_2PtBinned->SetParName(10,"Norm_{TailHighM}");
              fitV2_2PtBinned->SetParName(11,"Exp_{TailHighM}");
               fitV2_2PtBinned->SetParName(12,"V2_2 J/#psi");
               fitV2_2PtBinned->SetParName(13,"V2_2 Bkg M2");
           fitV2_2PtBinned->SetParName(14,"V2_2 Bkg M1");
           fitV2_2PtBinned->SetParName(15,"V2_2 Bkg M0");
            
            //PSI2STEST
            fitV2_2PtBinned->FixParameter(7,0);
            
        
           double minParams[16];
           double parErrors[16];
        
        //Fit of V2
            cout << "Fit of V2 wrt mass ptbinned, ptbin: "<<ptbin<<endl;
            res = V2JPsiTklPtBinned[ptbin]->Fit("fitV2_2PtBinned","SBMERI+","ep");//Trying this
            gStyle->SetOptStat("n");
            gStyle->SetOptFit(1011);
            if(!doTracklets){
            TPaveStats *st = (TPaveStats*)V2JPsiTklPtBinned[ptbin]->FindObject("stats");
            st->SetX1NDC(0.75); //new x start position
            st->SetY1NDC(0.75); //new x end position
            }
            Double_t para[16];
            fitV2_2PtBinned->GetParameters(para);
            backFcnV2_2->SetParameters(para);
            signalFcnJPsiV2_2->SetParameters(para);
        
        
        fitV2_2PtBinned->Draw("same");
      //  signalFcnJPsiV2_2->Draw("same");
        backFcnV2_2->Draw("same");
          // draw the legend
          TLegend *legend=new TLegend(0.12,0.80,0.60,0.90);
          legend->SetFillColorAlpha(kWhite, 0.);
          legend->SetBorderSize(0);
            legend->SetTextFont(42);
          legend->SetTextSize(0.03);
        Char_t message[80];
        sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitV2_2PtBinned->GetChisquare(),fitV2_2PtBinned->GetNDF());
          legend->AddEntry(fitV2_2PtBinned,message);

             if(res->CovMatrixStatus() == 3){
                   // sprintf(message,"The fit is a success");
             }
             else{
                    sprintf(message,"The fit is a failure");
                 legend->AddEntry(fitV2_2PtBinned,message);
             }
        // legend->AddEntry(signalFcnJPsiV2_2,"JPsi signal");
        legend->AddEntry(backFcnV2_2,"Background");
          legend->AddEntry(V2JPsiTklPtBinned[ptbin],"Data","lpe");
          legend->Draw();
            
            V2_Ext2[ptbin] = para[12];
            errV2_Ext2[ptbin] = fitV2_2PtBinned->GetParError(12);
        }
    }
    
    
    
    
    sprintf(CanvasName,"%s/V2_Extraction2PtBinned.pdf",FolderName);
    c16PtBinned->SaveAs(CanvasName);
    
    
    
    
    
    //Combining or not fit for Central-Periph method Ext1
    
    
    int niterations = 1;
    cinvmass->cd();
    sprintf(histoname,"InvMass_Central");
   res = FittingAllInvMassBin(histoname, cinvmass, 1); //ZZZZZZZ
    
    TCanvas* chnsegsigma = new TCanvas;
    chnsegsigma->cd();
    hnsegSigma->Draw("colz");

    double parerr[12];
    TCanvas* c10_fit = new TCanvas;
       c10_fit->Divide(NbinsDeltaPhi/2,4);
    TRandom rand = 0;
    for(int j=0; j<12; j++){
        parerr[j] = res->ParError(j);
    }
    for(int j=0; j <16; j++){
        par[j] = res->Parameter(j);
        if(j>11){
            par[j] = 3;
        }
    }
    //PLOT ICI

    double par12[niterations];
    double parerr12[niterations];
    
    for(int i=0; i<NbinsDeltaPhi; i++){
        for(int t=0; t<niterations; t++){
        TVirtualPad* subpad = c10_fit->cd(i+1);
            c10_fit->cd(i+1);
            if(t==niterations-1){
                subpad->Clear();
            }
        TF1 *fitY_1Central = new TF1("fitY_1Central",FourierV2_WrtInvMass,MinInvMass,MaxInvMass,16);
          fitY_1Central->SetNpx(500);
          fitY_1Central->SetLineWidth(4);
          fitY_1Central->SetLineColor(kMagenta);
          fitY_1Central->SetParameters(par);
           TVirtualFitter::Fitter(YieldWrtMass_Central[i])->SetMaxIterations(10000);
           TVirtualFitter::Fitter(YieldWrtMass_Central[i])->SetPrecision();
        //  histo->Fit("fitFcn","0");
          // second try: set start values for some parameters
        for(int k=0; k<=11; k++){
                fitY_1Central->FixParameter(k,par[k]);
        }

           fitY_1Central->SetParName(0,"Norm_{J/#psi}");
              fitY_1Central->SetParName(1,"M_{J/#psi}");
              fitY_1Central->SetParName(2,"Sigma_{J/#psi}");
              fitY_1Central->SetParName(3,"a_{1}");
              fitY_1Central->SetParName(4,"n_{1}");
              fitY_1Central->SetParName(5,"a_{2}");
              fitY_1Central->SetParName(6,"n_{2}");
              fitY_1Central->SetParName(7,"Norm_{#Psi(2S)}");
              fitY_1Central->SetParName(8,"Norm_{TailLowM}");
              fitY_1Central->SetParName(9,"Exp_{TailLowM}");
              fitY_1Central->SetParName(10,"Norm_{TailHighM}");
              fitY_1Central->SetParName(11,"Exp_{TailHighM}");
               fitY_1Central->SetParName(12,"Y_1 J/#psi");
               fitY_1Central->SetParName(13,"Y_1 Bkg M2");
           fitY_1Central->SetParName(14,"Y_1 Bkg M1");
           fitY_1Central->SetParName(15,"Y_1 Bkg M0");
            
            //PSI2STEST
            fitY_1Central->FixParameter(7,0);
            
            cout << "YieldWrtMass_Central in phibins. PhiBin: "<<i<<endl;

          rescent = YieldWrtMass_Central[i]->Fit("fitY_1Central","SBMERIQ+","ep");
            gStyle->SetOptStat("n");
            gStyle->SetOptFit(1011);
            if(!doTracklets){
            TPaveStats *st = (TPaveStats*)YieldWrtMass_Central[i]->FindObject("stats");
            st->SetX1NDC(0.75); //new x start position
            st->SetY1NDC(0.75); //new x end position
            }
          Double_t param[16];
            Double_t paramerrs[16];
          fitY_1Central->GetParameters(param);
            for(int i=0; i<16; i++){
                paramerrs[i] = fitY_1Central->GetParError(i);
            }
            
            if(rescent->CovMatrixStatus() == 3){
            par12[t] = param[12];
            parerr12[t] = paramerrs[12];
            }
            else{
                par12[t]=parerr12[t]=0;
            }
            
            std::cout << "COV status Central t=" << t << " : " << rescent->CovMatrixStatus()<<endl;
            std::cout << "par12: " << par12[t] << ", parerr12: " << parerr12[t] <<endl;

        TF1 *backFcnY_1Central = new TF1("backFcnY_1Central",BackFcnV2Poly,2.1,5.1,16);
        backFcnY_1Central->SetLineColor(kRed);
        TF1 *signalFcnJPsiY_1Central = new TF1("signalFcnJPsiY_1Central",SignalFcnJPsiV2,2.1,5.1,16);
          // writes the fit results into the par array
        
            if(t==niterations-1){
                YieldWrtMass_Central[i]->GetListOfFunctions()->Clear();
                YieldWrtMass_Central[i]->GetListOfFunctions()->Add(fitY_1Central);
            }
        signalFcnJPsiY_1Central->SetLineColor(kBlue);
        signalFcnJPsiY_1Central->SetNpx(500);
        backFcnY_1Central->SetParameters(param);
        signalFcnJPsiY_1Central->SetParameters(param);
      //  signalFcnJPsiY_1Central->Draw("same");
        backFcnY_1Central->Draw("same");
            
                Yields_Central_1->SetBinContent(i+1, param[12]);
                Yields_Central_1->SetBinError(i+1,fitY_1Central->GetParError(12));
            
          // writes the fit results into the par array
           gStyle->SetOptFit(1011);
          // draw the legend
          TLegend *legend=new TLegend(0.15,0.65,0.3,0.85);
          legend->SetFillColorAlpha(kWhite, 0.);
          legend->SetBorderSize(0);
            legend->SetTextFont(42);
          legend->SetTextSize(0.03);
        Char_t message[80];
        sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitY_1Central->GetChisquare(),fitY_1Central->GetNDF());
           legend->AddEntry(fitY_1Central,message);
          if(rescent->CovMatrixStatus() == 3){
                    // sprintf(message,"The fit is a success");
                 }
                 else{
                     sprintf(message,"The fit is a failure");
                     legend->AddEntry(fitY_1Central,message);
                 }
        //  legend->AddEntry(signalFcnJPsiY_1Central,"JPsi signal");
          legend->AddEntry(backFcnY_1Central,"Background");
            legend->AddEntry(YieldWrtMass_Central[i],"Data","lpe");
            legend->Draw();
          }
        
    }
    cinvmass->cd();
    sprintf(histoname,"InvMass_Neither");
    res = FittingAllInvMassBin(histoname, cinvmass, 3); //ZZZZZ
    
    cinvmass->cd();
    sprintf(histoname,"InvMass_Periph");
    res = FittingAllInvMassBin(histoname, cinvmass, 2); //ZZZZZ
    
       for(int j=0; j<12; j++){
           parerr[j] = res->ParError(j);
       }
       for(int j=0; j <16; j++){
           par[j] = res->Parameter(j);
           if(j>11){
               par[j] = 3;
           }
       }
    
    double baseline = 1;
    double errbaseline = 1;
       
       for(int i=0; i<NbinsDeltaPhi; i++){
           for(int t=0; t<niterations; t++){
               TVirtualPad* subpad = c10_fit->cd(NbinsDeltaPhi+i+1);
               c10_fit->cd(NbinsDeltaPhi+i+1);
               if(t==niterations-1){
                   subpad->Clear();
               }
           TF1 *fitY_1Periph = new TF1("fitY_1Periph",FourierV2_WrtInvMass,MinInvMass,MaxInvMass,16);
             fitY_1Periph->SetNpx(500);
             fitY_1Periph->SetLineWidth(4);
             fitY_1Periph->SetLineColor(kMagenta);
             fitY_1Periph->SetParameters(par);
              TVirtualFitter::Fitter(YieldWrtMass_Periph[i])->SetMaxIterations(10000);
              TVirtualFitter::Fitter(YieldWrtMass_Periph[i])->SetPrecision();
           //  histo->Fit("fitFcn","0");
             // second try: set start values for some parameters
            for(int k=0; k<=11; k++){

                    fitY_1Periph->FixParameter(k,par[k]);
                }
            
              
              fitY_1Periph->SetParName(0,"Norm_{J/#psi}");
                 fitY_1Periph->SetParName(1,"M_{J/#psi}");
                 fitY_1Periph->SetParName(2,"Sigma_{J/#psi}");
                 fitY_1Periph->SetParName(3,"a_{1}");
                 fitY_1Periph->SetParName(4,"n_{1}");
                 fitY_1Periph->SetParName(5,"a_{2}");
                 fitY_1Periph->SetParName(6,"n_{2}");
                 fitY_1Periph->SetParName(7,"Norm_{#Psi(2S)}");
                 fitY_1Periph->SetParName(8,"Norm_{TailLowM}");
                 fitY_1Periph->SetParName(9,"Exp_{TailLowM}");
                 fitY_1Periph->SetParName(10,"Norm_{TailHighM}");
                 fitY_1Periph->SetParName(11,"Exp_{TailHighM}");
                  fitY_1Periph->SetParName(12,"Y_1 J/#psi");
                  fitY_1Periph->SetParName(13,"Y_1 Bkg M2");
              fitY_1Periph->SetParName(14,"Y_1 Bkg M1");
              fitY_1Periph->SetParName(15,"Y_1 Bkg M0");
               
               //PSI2STEST
               fitY_1Periph->FixParameter(7,0);
              
               cout << "YieldWrtMass_Periph in phibins. PhiBin: "<<i<<endl;
             resperiph = YieldWrtMass_Periph[i]->Fit("fitY_1Periph","SBMERIQ+","ep");
               gStyle->SetOptStat("n");
               gStyle->SetOptFit(1011);
               if(!doTracklets){
               TPaveStats *st = (TPaveStats*)YieldWrtMass_Periph[i]->FindObject("stats");
               st->SetX1NDC(0.75); //new x start position
               st->SetY1NDC(0.75); //new x end position
               }
               Double_t param[16];
               Double_t paramerrs[16];
             fitY_1Periph->GetParameters(param);
               
               for(int i=0; i<16; i++){
                               paramerrs[i] = fitY_1Periph->GetParError(i);
                           }
                           
                           if(resperiph->CovMatrixStatus() == 3){
                           par12[t] = param[12];
                           parerr12[t] = paramerrs[12];
                           }
                           else{
                               par12[t]=parerr12[t]=0;
                           }
               std::cout << "COV status Central t=" << t << " : " << rescent->CovMatrixStatus()<<endl;
                std::cout << "par12: " << par12[t] << ", parerr12: " << parerr12[t] <<endl;

           TF1 *backFcnY_1Periph = new TF1("backFcnY_1Periph",BackFcnV2Poly,2.1,5.1,16);
           backFcnY_1Periph->SetLineColor(kRed);
           TF1 *signalFcnJPsiY_1Periph = new TF1("signalFcnJPsiY_1Periph",SignalFcnJPsiV2,2.1,5.1,16);
             // writes the fit results into the par array
               
               if(t==niterations-1){
                   YieldWrtMass_Periph[i]->GetListOfFunctions()->Clear();
                   YieldWrtMass_Periph[i]->GetListOfFunctions()->Add(fitY_1Periph);
               }
           
           signalFcnJPsiY_1Periph->SetLineColor(kBlue);
           signalFcnJPsiY_1Periph->SetNpx(500);
           backFcnY_1Periph->SetParameters(param);
           signalFcnJPsiY_1Periph->SetParameters(param);
        //   signalFcnJPsiY_1Periph->Draw("same");
           backFcnY_1Periph->Draw("same");
               

                   if(i==0){//Pichange
                       baseline = param[12];
                       errbaseline = fitY_1Periph->GetParError(12);
                   }
//                   if(i==3){
//                       baseline += param[12];
//                       baseline /= 2;
//                       errbaseline = sqrt(pow(errbaseline,2)+pow(fitY_1Periph->GetParError(12),2));
//                       errbaseline/=2;
//                   }
                 Yields_Periph_1->SetBinContent(i+1, param[12]);
                 Yields_Periph_1->SetBinError(i+1,fitY_1Periph->GetParError(12));
               
               
             // writes the fit results into the par array
              gStyle->SetOptFit(1011);
             // draw the legend
             TLegend *legend=new TLegend(0.15,0.65,0.3,0.85);
             legend->SetFillColorAlpha(kWhite, 0.);
             legend->SetBorderSize(0);
               legend->SetTextFont(42);
             legend->SetTextSize(0.04);
           Char_t message[80];
               sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitY_1Periph->GetChisquare(),fitY_1Periph->GetNDF());
                legend->AddEntry(fitY_1Periph,message);
               if(resperiph->CovMatrixStatus() == 3){
                         // sprintf(message,"The fit is a success");
                      }
                      else{
                          sprintf(message,"The fit is a failure");
                          legend->AddEntry(fitY_1Periph,message);
                      }
           //    legend->AddEntry(signalFcnJPsiY_1Periph,"JPsi signal");
               legend->AddEntry(backFcnY_1Periph,"Background");
                 legend->AddEntry(YieldWrtMass_Periph[i],"Data","lpe");
                 legend->Draw();
               }
           
       }
    
    sprintf(CanvasName,"%s/YieldWrtMass_PhiBins.pdf",FolderName);
    c10_fit->SaveAs(CanvasName);
    
    {
        TCanvas* c17 = new TCanvas;
        c17->Divide(1,3);
        //Creer canvas pour imprimer les plots Periph yield et Central yield wrt phi et leur difference
        c17->cd(1);
        Yields_Central_1->Draw();
        c17->cd(2);
        Yields_Periph_1->Draw();
        c17->cd(3);
        Yields_Difference_1->Add(Yields_Central_1,Yields_Periph_1,1,-1);
        Yields_Difference_1->Draw();
    
             TF1 *fitFcnV2_2 = new TF1("fitFcnV2_2",FourierV2,0,1.0*TMath::Pi(),3); //Pichange
             fitFcnV2_2->SetNpx(500);
             fitFcnV2_2->SetLineWidth(4);
             fitFcnV2_2->SetLineColor(kMagenta);
             // first try without starting values for the parameters
             // This defaults to 1 for each param.
             // this results in an ok fit for the polynomial function
             // however the non-linear part (lorenzian) does not
             // respond well.
              Double_t params[3] = {1,1,1};
             fitFcnV2_2->SetParameters(params);
              TVirtualFitter::Fitter(Yields_Difference_1)->SetMaxIterations(10000);
              TVirtualFitter::Fitter(Yields_Difference_1)->SetPrecision();
           //  histo->Fit("fitFcn","0");
             // second try: set start values for some parameters
              
              fitFcnV2_2->SetParName(0,"a0");
              fitFcnV2_2->SetParName(1,"a1");
              fitFcnV2_2->SetParName(2,"a2");
        
        cout << "YieldWrtMass_Difference " <<endl;
              
             TFitResultPtr res = Yields_Difference_1->Fit("fitFcnV2_2","SBMERI+","ep");
            gStyle->SetOptStat("n");
            gStyle->SetOptFit(1011);
        if(!doTracklets){
            TPaveStats *st = (TPaveStats*)Yields_Difference_1->FindObject("stats");
            st->SetX1NDC(0.75); //new x start position
            st->SetY1NDC(0.75); //new x end position
        }
            Double_t par[3];
            fitFcnV2_2->GetParameters(par);
             // improve the pictu
           //   std::cout << "integral error: " << integralerror << std::endl;
             fitFcnV2_2->Draw("same");
             // draw the legend
             TLegend *legend=new TLegend(0.15,0.65,0.3,0.85);
             legend->SetFillColorAlpha(kWhite, 0.);
             legend->SetBorderSize(0);
               legend->SetTextFont(42);
             legend->SetTextSize(0.03);
               Char_t message[80];
               sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitFcnV2_2->GetChisquare(),fitFcnV2_2->GetNDF());
               legend->AddEntry(fitFcnV2_2,message);
        sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitFcnV2_2->GetChisquare(),fitFcnV2_2->GetNDF());
               legend->AddEntry(fitFcnV2_2,message);
        sprintf(message,"V2_1 J/#psi-tkl: #frac{%.4f}{%.4f + %.4f} = %.4f +- %.4f",par[2],par[0], baseline,par[2]/(par[0] + baseline),(par[2]/(par[0] + baseline)*sqrt(pow(fitFcnV2_2->GetParError(2)/par[2],2)+pow(sqrt(pow(fitFcnV2_2->GetParError(0),2)+pow(errbaseline,2))/(par[0] + baseline),2))));
        legend->AddEntry(fitFcnV2_2,message);
        if(res->CovMatrixStatus() == 3){
                 //  sprintf(message,"The fit is a success");
               }
               else{
                   sprintf(message,"The fit is a failure");
                   legend->AddEntry(fitFcnV2_2,message);
               }
             legend->AddEntry(Yields_Difference_1,"Data","lpe");
             legend->Draw();
           
        sprintf(CanvasName,"%s/V2_Extraction1.pdf",FolderName);
        c17->SaveAs(CanvasName);
    }
    
    
    //DimuTkl Ext1 Pt Binned
        
    if(PtBinned){
        for(int ptbin=0;ptbin<NbPtBins;ptbin++){
    
        int niterations = 1;
        cinvmassPtBinned->cd();
        sprintf(histoname,Form("bin%d_1",ptbin+1));
        res = FittingAllInvMassBin(histoname, cinvmassPtBinned, NbPtBins + ptbin);

       double parerr[12];
        TCanvas* c10_fitPtBinned = new TCanvas;
           c10_fitPtBinned->Divide(NbinsDeltaPhi/2,4);
        TRandom rand = 0;
        for(int j=0; j<12; j++){
            parerr[j] = res->ParError(j);
        }
        for(int j=0; j <16; j++){
            par[j] = res->Parameter(j);
            if(j>11){
                par[j] = 3;
            }
        }
            //PLOT ICI

        double par12[niterations];
        double parerr12[niterations];

        
        for(int i=0; i<NbinsDeltaPhi; i++){
            for(int t=0; t<niterations; t++){
                TVirtualPad* subpad = c10_fitPtBinned->cd(i+1);
                c10_fitPtBinned->cd(i+1);
                if(t==niterations-1){
                    subpad->Clear();
                }
            TF1 *fitY_1Central = new TF1("fitY_1Central",FourierV2_WrtInvMass,1.5,4.5,16); //Cvetan
              fitY_1Central->SetNpx(50000);
              fitY_1Central->SetLineWidth(4);
              fitY_1Central->SetLineColor(kMagenta);
              fitY_1Central->SetParameters(par);
               TVirtualFitter::Fitter(YieldWrtMass_CentralPtBinned[i][ptbin])->SetMaxIterations(10000);
               TVirtualFitter::Fitter(YieldWrtMass_CentralPtBinned[i][ptbin])->SetPrecision();
            //  histo->Fit("fitFcn","0");
              // second try: set start values for some parameters
            for(int k=0; k<=11; k++){

                    fitY_1Central->FixParameter(k,par[k]);
                
            }

               fitY_1Central->SetParName(0,"Norm_{J/#psi}");
                  fitY_1Central->SetParName(1,"M_{J/#psi}");
                  fitY_1Central->SetParName(2,"Sigma_{J/#psi}");
                  fitY_1Central->SetParName(3,"a_{1}");
                  fitY_1Central->SetParName(4,"n_{1}");
                  fitY_1Central->SetParName(5,"a_{2}");
                  fitY_1Central->SetParName(6,"n_{2}");
                  fitY_1Central->SetParName(7,"Norm_{#Psi(2S)}");
                  fitY_1Central->SetParName(8,"Norm_{TailLowM}");
                  fitY_1Central->SetParName(9,"Exp_{TailLowM}");
                  fitY_1Central->SetParName(10,"Norm_{TailHighM}");
                  fitY_1Central->SetParName(11,"Exp_{TailHighM}");
                   fitY_1Central->SetParName(12,"Y_1 J/#psi");
                   fitY_1Central->SetParName(13,"Y_1 Bkg M2");
               fitY_1Central->SetParName(14,"Y_1 Bkg M1");
               fitY_1Central->SetParName(15,"Y_1 Bkg M0");
                
                //PSI2STEST
                fitY_1Central->FixParameter(7,0);
                
                cout << "YieldWrtMass_CentralPtBinned in phibins. PhiBin: "<<i<<" ptbin: "<<ptbin <<endl;
                
//                if(ptbin==2 && i==1){
//                    //Set valeurs Cvetan
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(2,2.855);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(2,0.00679);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(3,2.848);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(3,0.0076);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(4,2.867);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(4,0.00856);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(5,2.842);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(5,0.00966);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(6,2.83);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(6,0.00715);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(7,2.848);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(7,0.0135);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(8,2.86);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(8,0.0164);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinContent(9,2.83);
//                    YieldWrtMass_CentralPtBinned[i][ptbin]->SetBinError(9,0.0199);
//                }
                

              rescent = YieldWrtMass_CentralPtBinned[i][ptbin]->Fit("fitY_1Central","SBMERIQ+","ep");
                gStyle->SetOptStat("n");
                gStyle->SetOptFit(1011);
                if(!doTracklets){
                TPaveStats *st = (TPaveStats*)YieldWrtMass_CentralPtBinned[i][ptbin]->FindObject("stats");
                st->SetX1NDC(0.8); //new x start position
                st->SetY1NDC(0.8); //new x end position
                }
              Double_t param[16];
                Double_t paramerrs[16];
              fitY_1Central->GetParameters(param);
                for(int i=0; i<16; i++){
                    paramerrs[i] = fitY_1Central->GetParError(i);
                }
                
                if(rescent->CovMatrixStatus() == 3){
                par12[t] = param[12];
                parerr12[t] = paramerrs[12];
                }
                else{
                    par12[t]=parerr12[t]=0;
                }
                
                std::cout << "COV status Central t=" << t << " : " << rescent->CovMatrixStatus()<<endl;
                std::cout << "par12: " << par12[t] << ", parerr12: " << parerr12[t] <<endl;

            TF1 *backFcnY_1Central = new TF1("backFcnY_1Central",BackFcnV2Poly,2.1,5.1,16);
            backFcnY_1Central->SetLineColor(kRed);
            TF1 *signalFcnJPsiY_1Central = new TF1("signalFcnJPsiY_1Central",SignalFcnJPsiV2,2.1,5.1,16);
              // writes the fit results into the par array
            
                if(t==niterations-1){
                    YieldWrtMass_CentralPtBinned[i][ptbin]->GetListOfFunctions()->Clear();
                    YieldWrtMass_CentralPtBinned[i][ptbin]->GetListOfFunctions()->Add(fitY_1Central);
                }
            signalFcnJPsiY_1Central->SetLineColor(kBlue);
            signalFcnJPsiY_1Central->SetNpx(500);
            backFcnY_1Central->SetParameters(param);
            signalFcnJPsiY_1Central->SetParameters(param);
          //  signalFcnJPsiY_1Central->Draw("same");
            backFcnY_1Central->Draw("same");
                

                    Yields_Central_1PtBinned[ptbin]->SetBinContent(i+1, param[12]); //AEUGG
                    Yields_Central_1PtBinned[ptbin]->SetBinError(i+1,fitY_1Central->GetParError(12));
                
                
              // writes the fit results into the par array
               gStyle->SetOptFit(1011);
              // draw the legend
              TLegend *legend=new TLegend(0.12,0.75,0.60,0.90);
                legend->SetBorderSize(0);
              legend->SetTextFont(61);
              legend->SetTextSize(0.03);
            Char_t message[80];
            sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitY_1Central->GetChisquare(),fitY_1Central->GetNDF());
                legend->AddEntry(fitY_1Central,message);
               if(rescent->CovMatrixStatus() == 3){
                        //  sprintf(message,"The fit is a success");
                      }
                      else{
                          sprintf(message,"The fit is a failure");
                          legend->AddEntry(fitY_1Central,message);
                      }
            //   legend->AddEntry(signalFcnJPsiY_1Central,"JPsi signal");
               legend->AddEntry(backFcnY_1Central,"Background");
                 legend->AddEntry(YieldWrtMass_CentralPtBinned[i][ptbin],"Data","lpe");
                 legend->Draw();
            }
            
        }
        
            
        cinvmassPtBinned->cd(); //aeugg
        sprintf(histoname,Form("bin%d_2",ptbin+1));
        res = FittingAllInvMassBin(histoname, cinvmassPtBinned, NbPtBins+NbPtBins+ptbin);
        
           for(int j=0; j<12; j++){
               parerr[j] = res->ParError(j);
           }
           for(int j=0; j <16; j++){
               par[j] = res->Parameter(j);
               if(j>11){
                   par[j] = 3;
               }
           }
        
        double baseline = 1;
        double errbaseline = 1;
           
           for(int i=0; i<NbinsDeltaPhi; i++){
               for(int t=0; t<niterations; t++){
                   TVirtualPad* subpad = c10_fitPtBinned->cd(NbinsDeltaPhi+i+1);
                   c10_fitPtBinned->cd(NbinsDeltaPhi+i+1);
                   if(t==niterations-1){
                       subpad->Clear();
                   }
               TF1 *fitY_1Periph = new TF1("fitY_1Periph",FourierV2_WrtInvMass,MinInvMass,MaxInvMass,16);
                 fitY_1Periph->SetNpx(500);
                 fitY_1Periph->SetLineWidth(4);
                 fitY_1Periph->SetLineColor(kMagenta);
                 fitY_1Periph->SetParameters(par);
                  TVirtualFitter::Fitter(YieldWrtMass_PeriphPtBinned[i][ptbin])->SetMaxIterations(10000);
                  TVirtualFitter::Fitter(YieldWrtMass_PeriphPtBinned[i][ptbin])->SetPrecision();
               //  histo->Fit("fitFcn","0");
                 // second try: set start values for some parameters
                for(int k=0; k<=11; k++){
                        fitY_1Periph->FixParameter(k,par[k]);
                }
                  
                  fitY_1Periph->SetParName(0,"Norm_{J/#psi}");
                     fitY_1Periph->SetParName(1,"M_{J/#psi}");
                     fitY_1Periph->SetParName(2,"Sigma_{J/#psi}");
                     fitY_1Periph->SetParName(3,"a_{1}");
                     fitY_1Periph->SetParName(4,"n_{1}");
                     fitY_1Periph->SetParName(5,"a_{2}");
                     fitY_1Periph->SetParName(6,"n_{2}");
                     fitY_1Periph->SetParName(7,"Norm_{#Psi(2S)}");
                     fitY_1Periph->SetParName(8,"Norm_{TailLowM}");
                     fitY_1Periph->SetParName(9,"Exp_{TailLowM}");
                     fitY_1Periph->SetParName(10,"Norm_{TailHighM}");
                     fitY_1Periph->SetParName(11,"Exp_{TailHighM}");
                      fitY_1Periph->SetParName(12,"Y_1 J/#psi");
                      fitY_1Periph->SetParName(13,"Y_1 Bkg M2");
                  fitY_1Periph->SetParName(14,"Y_1 Bkg M1");
                  fitY_1Periph->SetParName(15,"Y_1 Bkg M0");
                   
                   //PSI2STEST
                   fitY_1Periph->FixParameter(7,0);
                  
                   cout << "YieldWrtMass_PeriphPtBinned in phibins. PhiBin: "<<i<<" ptbin: "<<ptbin <<endl;
                 resperiph = YieldWrtMass_PeriphPtBinned[i][ptbin]->Fit("fitY_1Periph","SBMERIQ+","ep");
                   gStyle->SetOptStat("n");
                   gStyle->SetOptFit(1011);
                   if(!doTracklets){
                   TPaveStats *st = (TPaveStats*)YieldWrtMass_PeriphPtBinned[i][ptbin]->FindObject("stats");
                   st->SetX1NDC(0.8); //new x start position
                   st->SetY1NDC(0.8); //new x end position
                   }
                   Double_t param[16];
                   Double_t paramerrs[16];
                 fitY_1Periph->GetParameters(param);
                   
                   for(int i=0; i<16; i++){
                                   paramerrs[i] = fitY_1Periph->GetParError(i);
                               }
                               
                               if(resperiph->CovMatrixStatus() == 3){
                               par12[t] = param[12];
                               parerr12[t] = paramerrs[12];
                               }
                               else{
                                   par12[t]=parerr12[t]=0;
                               }
                   std::cout << "COV status Central t=" << t << " : " << rescent->CovMatrixStatus()<<endl;
                    std::cout << "par12: " << par12[t] << ", parerr12: " << parerr12[t] <<endl;

               TF1 *backFcnY_1Periph = new TF1("backFcnY_1Periph",BackFcnV2Poly,2.1,5.1,16);
               backFcnY_1Periph->SetLineColor(kRed);
               TF1 *signalFcnJPsiY_1Periph = new TF1("signalFcnJPsiY_1Periph",SignalFcnJPsiV2,2.1,5.1,16);
                 // writes the fit results into the par array
                   
                   if(t==niterations-1){
                       YieldWrtMass_PeriphPtBinned[i][ptbin]->GetListOfFunctions()->Clear();
                       YieldWrtMass_PeriphPtBinned[i][ptbin]->GetListOfFunctions()->Add(fitY_1Periph);
                   }
               
               signalFcnJPsiY_1Periph->SetLineColor(kBlue);
               signalFcnJPsiY_1Periph->SetNpx(500);
               backFcnY_1Periph->SetParameters(param);
               signalFcnJPsiY_1Periph->SetParameters(param);
            //   signalFcnJPsiY_1Periph->Draw("same");
               backFcnY_1Periph->Draw("same");
                   
                       if(i==0){//Pichange
                           baseline = param[12];
                           errbaseline = fitY_1Periph->GetParError(12);
                       }
//                       if(i==3){//Pichange
//                           baseline += param[12];
//                           baseline /= 2;
//                           errbaseline = sqrt(pow(errbaseline,2)+pow(fitY_1Periph->GetParError(12),2));
//                           errbaseline/=2;
//                       }
                     Yields_Periph_1PtBinned[ptbin]->SetBinContent(i+1, param[12]);
                     Yields_Periph_1PtBinned[ptbin]->SetBinError(i+1,fitY_1Periph->GetParError(12));
                   
                   
                 // writes the fit results into the par array
                  gStyle->SetOptFit(1011);
                 // draw the legend
                 TLegend *legend=new TLegend(0.12,0.75,0.60,0.90);
                 legend->SetFillColorAlpha(kWhite, 0.);
                 legend->SetBorderSize(0);
                   legend->SetTextFont(42);
                 legend->SetTextSize(0.03);
               Char_t message[80];
               sprintf(message,"Global Fit : #chi^{2}/NDF = %.2f / %d",fitY_1Periph->GetChisquare(),fitY_1Periph->GetNDF());
                legend->AddEntry(fitY_1Periph,message);
               if(resperiph->CovMatrixStatus() == 3){
                          sprintf(message,"The fit is a success");
                      }
                      else{
                          sprintf(message,"The fit is a failure");
                      }
                      legend->AddEntry(fitY_1Periph,message);
          //     legend->AddEntry(signalFcnJPsiY_1Periph,"JPsi signal");
               legend->AddEntry(backFcnY_1Periph,"Background");
                 legend->AddEntry(YieldWrtMass_Periph[i],"Data","lpe");
                 legend->Draw();
               }
               
           }
            
            
            
            {
                TCanvas* c17PtBinned = new TCanvas;
                c17PtBinned->Divide(1,3);
                //Creer canvas pour imprimer les plots Periph yield et Central yield wrt phi et leur difference
                c17PtBinned->cd(1);
                Yields_Central_1PtBinned[ptbin]->Draw();
                c17PtBinned->cd(2);
                Yields_Periph_1PtBinned[ptbin]->Draw();
                c17PtBinned->cd(3);
                Yields_Difference_1PtBinned[ptbin]->Add(Yields_Central_1PtBinned[ptbin],Yields_Periph_1PtBinned[ptbin],1,-1);
                Yields_Difference_1PtBinned[ptbin]->Draw();
            
                     TF1 *fitFcnV2_2 = new TF1("fitFcnV2_2",FourierV2,0,1.0*TMath::Pi(),3); //Pichange
                     fitFcnV2_2->SetNpx(500);
                     fitFcnV2_2->SetLineWidth(4);
                     fitFcnV2_2->SetLineColor(kMagenta);
                     // first try without starting values for the parameters
                     // This defaults to 1 for each param.
                     // this results in an ok fit for the polynomial function
                     // however the non-linear part (lorenzian) does not
                     // respond well.
                      Double_t params[3] = {1,1,1};
                     fitFcnV2_2->SetParameters(params);
                      TVirtualFitter::Fitter(Yields_Difference_1PtBinned[ptbin])->SetMaxIterations(10000);
                      TVirtualFitter::Fitter(Yields_Difference_1PtBinned[ptbin])->SetPrecision();
                   //  histo->Fit("fitFcn","0");
                     // second try: set start values for some parameters
                      
                      fitFcnV2_2->SetParName(0,"a0");
                      fitFcnV2_2->SetParName(1,"a1");
                      fitFcnV2_2->SetParName(2,"a2");
                      
                cout << "Yields_Difference_1PtBinned in ptbin: "<<ptbin <<endl;
                     TFitResultPtr res = Yields_Difference_1PtBinned[ptbin]->Fit("fitFcnV2_2","SBMERI+","ep");
                gStyle->SetOptStat("n");
                gStyle->SetOptFit(1011);
                if(!doTracklets){
                TPaveStats *st = (TPaveStats*)Yields_Difference_1PtBinned[ptbin]->FindObject("stats");
                st->SetX1NDC(0.8); //new x start position
                st->SetY1NDC(0.8); //new x end position
                }
                    Double_t par[3];
                    fitFcnV2_2->GetParameters(par);
                     // improve the pictu
                   //   std::cout << "integral error: " << integralerror << std::endl;
                     fitFcnV2_2->Draw("same");
                     // draw the legend
                     TLegend *legend=new TLegend(0.15,0.70,0.4,0.90);
                     legend->SetFillColorAlpha(kWhite, 0.);
                     legend->SetBorderSize(0);
                       legend->SetTextFont(42);
                     legend->SetTextSize(0.03);
                       Char_t message[80];
                       sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitFcnV2_2->GetChisquare(),fitFcnV2_2->GetNDF());
                       legend->AddEntry(fitFcnV2_2,message);
                sprintf(message,"V2_1 J/#psi-tkl: #frac{%.4f}{%.4f + %.4f} = %.4f +- %.4f",par[2],par[0], baseline,par[2]/(par[0] + baseline),(par[2]/(par[0] + baseline)*sqrt(pow(fitFcnV2_2->GetParError(2)/par[2],2)+pow(sqrt(pow(fitFcnV2_2->GetParError(0),2)+pow(errbaseline,2))/(par[0] + baseline),2))));
                legend->AddEntry(fitFcnV2_2,message);
                if(res->CovMatrixStatus() == 3){
                          // sprintf(message,"The fit is a success");
                       }
                       else{
                           sprintf(message,"The fit is a failure");
                           legend->AddEntry(fitFcnV2_2,message);
                       }
                       
                     legend->AddEntry(Yields_Difference_1,"Data","lpe");
                     legend->Draw();
                
                V2_Ext1[ptbin] = par[2]/(par[0] + baseline);
                errV2_Ext1[ptbin] = abs(par[2]/(par[0] + baseline)*sqrt(pow(fitFcnV2_2->GetParError(2)/par[2],2)+pow(sqrt(pow(fitFcnV2_2->GetParError(0),2)+pow(errbaseline,2))/(par[0] + baseline),2)));
                   
                
                sprintf(CanvasName,"%s/V2_Extraction1PtBin_%d.pdf",FolderName,ptbin);
                c17PtBinned->SaveAs(CanvasName);
            }
    
    
            sprintf(CanvasName,"%s/YieldWrtMass_PhiBins_PtBin_%d.pdf",FolderName,ptbin);
            c10_fitPtBinned->SaveAs(CanvasName);
        }
    }
    
    
    double PtMiddle[NbPtBins] = {0};
    double PtErrorSize[NbPtBins] = {0};
    
    for(int ptbin=0;ptbin<NbPtBins;ptbin++){
        PtMiddle[ptbin] = (PtBins[ptbin]+PtBins[ptbin+1])/2;
        PtErrorSize[ptbin] = (PtBins[ptbin+1]-PtBins[ptbin])/2;
    }
    
    TCanvas* cV2Pt = new TCanvas;
    cV2Pt->cd();
    
    TGraphErrors *grV2_wrt_Pt1 = new TGraphErrors(NbPtBins,PtMiddle,V2_Ext1,PtErrorSize,errV2_Ext1);
             // TGraph *gr3 = new TGraph (n, K3, chi);
              grV2_wrt_Pt1->SetTitle("V_{2,J/#psi-tkl} wrt p_{T} for different extraction methods");
              grV2_wrt_Pt1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
              grV2_wrt_Pt1->GetYaxis()->SetTitle("V_{2,J/#psi-tkl}");
            grV2_wrt_Pt1->SetMarkerColor(kAzure-3);
            grV2_wrt_Pt1->SetLineColor(kAzure-3);
            grV2_wrt_Pt1->SetLineStyle(9);
            grV2_wrt_Pt1->SetMarkerStyle(4);
            grV2_wrt_Pt1->GetYaxis()->SetRangeUser(-0.004,0.014);
            grV2_wrt_Pt1->GetXaxis()->SetRangeUser(PtBins[0],PtBins[NbPtBins]);
              grV2_wrt_Pt1->Draw("AP");
    
    TGraphErrors *grV2_wrt_Pt2 = new TGraphErrors(NbPtBins,PtMiddle,V2_Ext2,PtErrorSize,errV2_Ext2);
     // TGraph *gr3 = new TGraph (n, K3, chi);
    grV2_wrt_Pt2->SetMarkerColor(kGreen-7);
    grV2_wrt_Pt2->SetLineColor(kGreen-7);
    grV2_wrt_Pt2->SetLineStyle(9);
    grV2_wrt_Pt2->SetMarkerStyle(32);
      grV2_wrt_Pt2->Draw("P");
    
    TLegend *legendo=new TLegend(0.12,0.80,0.40,0.90);
             legendo->SetFillColorAlpha(kWhite, 0.);
             legendo->SetBorderSize(0);
              legendo->SetTextFont(42);
              legendo->SetTextSize(0.02);
           legendo->AddEntry(grV2_wrt_Pt1,"V2_Ext1 (pPb-like)");
            legendo->AddEntry(grV2_wrt_Pt2,"V2_Ext2 (pPb-like)");
             legendo->Draw();
    
    cV2Pt->Update();
    TLine *l=new TLine(cV2Pt->GetUxmin(),0.0,cV2Pt->GetUxmax(),0.0);
    l->SetLineColor(kBlack);
    l->SetLineWidth(1);
    l->SetLineStyle(9);
    l->Draw();
    
    sprintf(CanvasName,"%s/V2WrtPt.pdf",FolderName);
    cV2Pt->SaveAs(CanvasName);
    
    sprintf(CanvasName,"%s/InvMassFits.pdf",FolderName);
    cinvmass->SaveAs(CanvasName);
    sprintf(CanvasName,"%s/InvMassFitsPtBinned.pdf",FolderName);
    cinvmassPtBinned->SaveAs(CanvasName);
    sprintf(CanvasName,"%s/YieldWrtMass_PhiBins_NoFit.pdf",FolderName);
    c10->SaveAs(CanvasName);
    sprintf(CanvasName,"%s/YieldWrPhi_MassBins.pdf",FolderName);
    c14->SaveAs(CanvasName);
    sprintf(CanvasName,"%s/YieldWrtPhi_MassBin_SubZoom.pdf",FolderName);
    c8zoom->SaveAs(CanvasName);
    sprintf(CanvasName,"%s/YieldWrtPhi_MassBin_Sub.pdf",FolderName);
    c8->SaveAs(CanvasName);
    
    
    
    double baselineTKL = 1;
    double errbaselineTKL = 1;
    
    if(doTracklets){
        
        TCanvas* c6TKL = new TCanvas;
            //Tracklets yield DeltaEta wrt DeltaPhi TH2 -> Projected for Periph and Central

            c6TKL->Divide(2,2);
            c6TKL->cd(1);
            YieldTkl_allC->Draw("E");
            c6TKL->cd(2);
            YieldTkl_Difference->Draw("E");
            c6TKL->cd(3);
            YieldTkl_Central->Draw("E");
            c6TKL->cd(4);
            YieldTkl_Periph->Draw("E");
            c6TKL->Draw();
            c6TKL->Modified();
            c6TKL->ForceUpdate();
        
        sprintf(CanvasName,"%s/V2TKL_Sub.pdf",FolderName);
        c6TKL->SaveAs(CanvasName);
        
        TFile *f = new TFile(FitFileName,"UPDATE");
        YieldTkl_allC->Write();
        YieldTkl_Periph->Write();
        YieldTkl_Central->Write();
        YieldTkl_Difference->Write();
        f->Close();
        
        baselineTKL = (YieldTkl_Periph->GetBinContent(BinZeroLeftTKL) + YieldTkl_Periph->GetBinContent(BinZeroLeftTKL+1))/2;
        errbaselineTKL = sqrt(pow(YieldTkl_Periph->GetBinError(BinZeroLeftTKL),2) + pow(YieldTkl_Periph->GetBinError(BinZeroLeftTKL+1),2))/2;
            
        TCanvas*c14TKL=new TCanvas();
        //Tracklets Yield difference wrt Phi fit
        c14TKL->cd();
        // Ici on fit YieldsWrtDeltaPhiMassBin_DifferenceProj
            TH1F *histo = YieldTkl_Difference;
          // create a TF1 with the range from 0 to 3 and 6 parameters
          TF1 *fitFcnV2TKL = new TF1("fitFcnV2TKL",FourierV2,0,1.0*TMath::Pi(),3); //pichange
          fitFcnV2TKL->SetNpx(500);
          fitFcnV2TKL->SetLineWidth(4);
          fitFcnV2TKL->SetLineColor(kMagenta);
          // first try without starting values for the parameters
          // This defaults to 1 for each param.
          // this results in an ok fit for the polynomial function
          // however the non-linear part (lorenzian) does not
          // respond well.
           Double_t params[3] = {1,0.01,0.01};
          fitFcnV2TKL->SetParameters(params);
           TVirtualFitter::Fitter(histo)->SetMaxIterations(10000);
           TVirtualFitter::Fitter(histo)->SetPrecision();
        //  histo->Fit("fitFcn","0");
          // second try: set start values for some parameters
           
           fitFcnV2TKL->SetParName(0,"a0");
           fitFcnV2TKL->SetParName(1,"a1");
           fitFcnV2TKL->SetParName(2,"a2");
            fitFcnV2TKL->SetParName(3,"a3");
//            fitFcnV2TKL->SetParName(4,"a4");
//        fitFcnV2TKL->SetParName(5,"a5");
//        fitFcnV2TKL->SetParName(6,"a6");
//        fitFcnV2TKL->SetParName(7,"a7");
//        fitFcnV2TKL->SetParName(8,"a8");
//        fitFcnV2TKL->SetParName(9,"a9");
//        fitFcnV2TKL->SetParName(10,"a10");
//        fitFcnV2TKL->SetParName(11,"a11");
//        fitFcnV2TKL->SetParLimits(6,-0.1,0.1);
           
          TFitResultPtr res = histo->Fit("fitFcnV2TKL","SBMERI+","ep");
        gStyle->SetOptStat("n");
        gStyle->SetOptFit(1011);
        if(!doTracklets){
        TPaveStats *st = (TPaveStats*)histo->FindObject("stats");
        st->SetX1NDC(0.8); //new x start position
        st->SetY1NDC(0.8); //new x end position
        }
          // improve the pictu
        //   std::cout << "integral error: " << integralerror << std::endl;
            Double_t par[3];
            fitFcnV2TKL->GetParameters(par);
          fitFcnV2TKL->Draw("same");
          // draw the legend
          TLegend *legend=new TLegend(0.12,0.80,0.60,0.90);
          legend->SetFillColorAlpha(kWhite, 0.);
          legend->SetBorderSize(0);
            legend->SetTextFont(42);
          legend->SetTextSize(0.03);
            Char_t message[80];
            sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitFcnV2TKL->GetChisquare(),fitFcnV2TKL->GetNDF());
            legend->AddEntry(fitFcnV2TKL,message);
        sprintf(message,"V_{2,tkl-tkl}: #frac{%.4f}{%.4f + %.4f} = %.4f +- %.4f",par[2],par[0], baselineTKL,par[2]/(par[0] + baselineTKL),(par[2]/(par[0] + baselineTKL)*sqrt(pow(fitFcnV2TKL->GetParError(2)/par[2],2)+pow(sqrt(pow(fitFcnV2TKL->GetParError(0),2)+pow(errbaselineTKL,2))/(par[0] + baselineTKL),2))));
        legend->AddEntry(fitFcnV2TKL,message);
        if(res->CovMatrixStatus() == 3){
                 //  sprintf(message,"The fit is a success");
               }
               else{
                   sprintf(message,"The fit is a failure");
                   legend->AddEntry(fitFcnV2TKL,message);
               }
          legend->AddEntry(histo,"Data","lpe");
          legend->Draw();
        
        sprintf(CanvasName,"%s/V2TKL.pdf",FolderName);
        c14TKL->SaveAs(CanvasName);
        
    }
    
    Double_t Acoef_Periph[4];
    Double_t Acoef_Central[4];
    Double_t errAcoef_Periph[4];
    Double_t errAcoef_Central[4];
    
    if (doTracklets && isCMSMethod){
            
            TCanvas*c14TKLperiph=new TCanvas();
            //Tracklets Yield difference wrt Phi fit
            c14TKLperiph->cd();
            // Ici on fit YieldsWrtDeltaPhiMassBin_DifferenceProj
                TH1F *histo = YieldTkl_Periph;
              // create a TF1 with the range from 0 to 3 and 6 parameters
              TF1 *fitFcnV2TKL = new TF1("fitFcnV2TKL",FourierV5,0,1.0*TMath::Pi(),4);//Pichange
              fitFcnV2TKL->SetNpx(500);
              fitFcnV2TKL->SetLineWidth(4);
              fitFcnV2TKL->SetLineColor(kMagenta);
              // first try without starting values for the parameters
              // This defaults to 1 for each param.
              // this results in an ok fit for the polynomial function
              // however the non-linear part (lorenzian) does not
              // respond well.
               Double_t params[4] = {1,0.01,0.01, 0.01};
              fitFcnV2TKL->SetParameters(params);
               TVirtualFitter::Fitter(histo)->SetMaxIterations(10000);
               TVirtualFitter::Fitter(histo)->SetPrecision();
            //  histo->Fit("fitFcn","0");
              // second try: set start values for some parameters
               
               fitFcnV2TKL->SetParName(0,"a0");
               fitFcnV2TKL->SetParName(1,"a1");
               fitFcnV2TKL->SetParName(2,"a2");
                fitFcnV2TKL->SetParName(3,"a3");
    //            fitFcnV2TKL->SetParName(4,"a4");
    //        fitFcnV2TKL->SetParName(5,"a5");
    //        fitFcnV2TKL->SetParName(6,"a6");
    //        fitFcnV2TKL->SetParName(7,"a7");
    //        fitFcnV2TKL->SetParName(8,"a8");
    //        fitFcnV2TKL->SetParName(9,"a9");
    //        fitFcnV2TKL->SetParName(10,"a10");
    //        fitFcnV2TKL->SetParName(11,"a11");
    //        fitFcnV2TKL->SetParLimits(6,-0.1,0.1);
               
              TFitResultPtr res = histo->Fit("fitFcnV2TKL","SBMERI+","ep");
        gStyle->SetOptStat("n");
        gStyle->SetOptFit(1011);
        if(!doTracklets){
        TPaveStats *st = (TPaveStats*)histo->FindObject("stats");
        st->SetX1NDC(0.8); //new x start position
        st->SetY1NDC(0.8); //new x end position
        }
              // improve the pictu
            //   std::cout << "integral error: " << integralerror << std::endl;
                fitFcnV2TKL->GetParameters(Acoef_Periph);
        for(int index=0; index<4;index++){
            errAcoef_Periph[index] = fitFcnV2TKL->GetParError(index);
        }
              fitFcnV2TKL->Draw("same");
              // draw the legend
              TLegend *legend=new TLegend(0.12,0.80,0.60,0.90);
              legend->SetFillColorAlpha(kWhite, 0.);
              legend->SetBorderSize(0);
                legend->SetTextFont(42);
              legend->SetTextSize(0.03);
                Char_t message[80];
                sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitFcnV2TKL->GetChisquare(),fitFcnV2TKL->GetNDF());
                legend->AddEntry(fitFcnV2TKL,message);
            sprintf(message,"V_{2,tkl-tkl} Periph: #frac {%.4f}{%.4f} = %.4f +- %.4f",Acoef_Periph[2],Acoef_Periph[0],Acoef_Periph[2]/(Acoef_Periph[0]),(Acoef_Periph[2]/(Acoef_Periph[0])*sqrt(pow(fitFcnV2TKL->GetParError(2)/Acoef_Periph[2],2)+pow(fitFcnV2TKL->GetParError(0)/Acoef_Periph[0],2))));
            legend->AddEntry(fitFcnV2TKL,message);
            if(res->CovMatrixStatus() == 3){
                       //sprintf(message,"The fit is a success");
                   }
                   else{
                       sprintf(message,"The fit is a failure");
                       legend->AddEntry(fitFcnV2TKL,message);
                   }
              legend->AddEntry(histo,"Data","lpe");
              legend->Draw();
        
        }
    
    if (doTracklets && isCMSMethod){
            
            TCanvas*c14TKLcent=new TCanvas();
            //Tracklets Yield difference wrt Phi fit
            c14TKLcent->cd();
            // Ici on fit YieldsWrtDeltaPhiMassBin_DifferenceProj
                TH1F *histo = YieldTkl_Central;
              // create a TF1 with the range from 0 to 3 and 6 parameters
              TF1 *fitFcnV2TKL = new TF1("fitFcnV2TKL",FourierV5,0,1.0*TMath::Pi(),4);//Pichange
              fitFcnV2TKL->SetNpx(500);
              fitFcnV2TKL->SetLineWidth(4);
              fitFcnV2TKL->SetLineColor(kMagenta);
              // first try without starting values for the parameters
              // This defaults to 1 for each param.
              // this results in an ok fit for the polynomial function
              // however the non-linear part (lorenzian) does not
              // respond well.
               Double_t params[4] = {1,0.01,0.01, 0.01};
              fitFcnV2TKL->SetParameters(params);
               TVirtualFitter::Fitter(histo)->SetMaxIterations(10000);
               TVirtualFitter::Fitter(histo)->SetPrecision();
            //  histo->Fit("fitFcn","0");
              // second try: set start values for some parameters
               
               fitFcnV2TKL->SetParName(0,"a0");
               fitFcnV2TKL->SetParName(1,"a1");
               fitFcnV2TKL->SetParName(2,"a2");
                fitFcnV2TKL->SetParName(3,"a3");
    //            fitFcnV2TKL->SetParName(4,"a4");
    //        fitFcnV2TKL->SetParName(5,"a5");
    //        fitFcnV2TKL->SetParName(6,"a6");
    //        fitFcnV2TKL->SetParName(7,"a7");
    //        fitFcnV2TKL->SetParName(8,"a8");
    //        fitFcnV2TKL->SetParName(9,"a9");
    //        fitFcnV2TKL->SetParName(10,"a10");
    //        fitFcnV2TKL->SetParName(11,"a11");
    //        fitFcnV2TKL->SetParLimits(6,-0.1,0.1);
               
              TFitResultPtr res = histo->Fit("fitFcnV2TKL","SBMERI+","ep");
            gStyle->SetOptStat("n");
            gStyle->SetOptFit(1011);
        if(!doTracklets){
            TPaveStats *st = (TPaveStats*)histo->FindObject("stats");
            st->SetX1NDC(0.8); //new x start position
            st->SetY1NDC(0.8); //new x end position
        }
              // improve the pictu
            //   std::cout << "integral error: " << integralerror << std::endl;
                fitFcnV2TKL->GetParameters(Acoef_Central);
                for(int index=0; index<4;index++){
                    errAcoef_Central[index] = fitFcnV2TKL->GetParError(index);
                }
              fitFcnV2TKL->Draw("same");
              // draw the legend
              TLegend *legend=new TLegend(0.12,0.80,0.60,0.90);
              legend->SetFillColorAlpha(kWhite, 0.);
              legend->SetBorderSize(0);
                legend->SetTextFont(42);
              legend->SetTextSize(0.03);
                Char_t message[80];
                sprintf(message,"Fit : #chi^{2}/NDF = %.2f / %d",fitFcnV2TKL->GetChisquare(),fitFcnV2TKL->GetNDF());
                legend->AddEntry(fitFcnV2TKL,message);
            sprintf(message,"V_{2,tkl-tkl} Central: #frac{%.4f}{%.4f} = %.4f +- %.4f",Acoef_Central[2],Acoef_Central[0],Acoef_Central[2]/(Acoef_Central[0]),(Acoef_Central[2]/(Acoef_Central[0])*sqrt(pow(fitFcnV2TKL->GetParError(2)/Acoef_Central[2],2)+pow(fitFcnV2TKL->GetParError(0)/Acoef_Central[0],2))));
            legend->AddEntry(fitFcnV2TKL,message);
            if(res->CovMatrixStatus() == 3){
                     //  sprintf(message,"The fit is a success");
                   }
                   else{
                       sprintf(message,"The fit is a failure");
                       legend->AddEntry(fitFcnV2TKL,message);
                   }
              legend->AddEntry(histo,"Data","lpe");
              legend->Draw();
        
        }
    
    if(doTracklets && isCMSMethod){
        Double_t V2_CMS_Periph[4];
        Double_t V2_CMS_Central[4];
        Double_t V2_CMS_Difference[4];
        Double_t errV2_CMS_Periph[4];
        Double_t errV2_CMS_Central[4];
        Double_t errV2_CMS_Difference[4];
        Double_t errCompound[4];
        
        for(int index=0; index<4; index++){
            cout << "Acoef_Periph_" << index << " = " << Acoef_Periph[index] <<" +/- "<<errAcoef_Periph[index]<<endl;
        }
        for(int index=0; index<4; index++){
                   cout << "Acoef_Central_" << index << " = " << Acoef_Central[index] <<" +/- "<< errAcoef_Central[index] <<endl;
               }
        
        //Calculation V2 Periph and Central
        for(int index=0; index<4; index++){
            V2_CMS_Periph[index]=Acoef_Periph[index]/Acoef_Periph[0];
            V2_CMS_Central[index]=Acoef_Central[index]/Acoef_Central[0];
            
            errV2_CMS_Periph[index] = V2_CMS_Periph[index]*sqrt(pow(errAcoef_Periph[index]/Acoef_Periph[index],2)+pow(errAcoef_Periph[0]/Acoef_Periph[0],2));
            errV2_CMS_Central[index] = V2_CMS_Central[index]*sqrt(pow(errAcoef_Central[index]/Acoef_Central[index],2)+pow(errAcoef_Central[0]/Acoef_Central[0],2));
        }
        
        for(int index=0; index<4; index++){
                   cout << "V2_Periph_" << index << " = " << V2_CMS_Periph[index] <<" +/- "<<errV2_CMS_Periph[index]<<endl;
               }
               for(int index=0; index<4; index++){
                          cout << "V2_Central_" << index << " = " << V2_CMS_Central[index] <<" +/- "<< errV2_CMS_Central[index] <<endl;
                      }
        
        TCanvas* ctry=new TCanvas();
        ctry->Divide(1,2);
        ctry->cd(1);
        YieldTkl_CentralLongCMS->DrawCopy();
        ctry->cd(2);
        YieldTkl_CentralShortCMS->DrawCopy();
        
        TCanvas* ctry2=new TCanvas();
        ctry2->Divide(1,2);
        ctry2->cd(1);
        YieldTkl_PeriphLongCMS->DrawCopy();
        ctry2->cd(2);
        YieldTkl_PeriphShortCMS->DrawCopy();
        
        //Calculation Y_jet_Central and Y_jet_Periph
        Double_t e1;
        Double_t e2;
        Double_t e3;
        Double_t e4;
        
        Double_t Y_jet_Central = -(YieldTkl_CentralLongCMS->IntegralAndError(1,YieldTkl_CentralLongCMS->GetNbinsX(),e1) - YieldTkl_CentralShortCMS->IntegralAndError(1,YieldTkl_CentralShortCMS->GetNbinsX(),e2));
        Double_t Y_jet_Periph = -(YieldTkl_PeriphLongCMS->IntegralAndError(1,YieldTkl_PeriphLongCMS->GetNbinsX(),e3) - YieldTkl_PeriphShortCMS->IntegralAndError(1,YieldTkl_PeriphShortCMS->GetNbinsX(),e4)) ;
        Double_t errY_jet_Central = sqrt(pow(e1,2)+pow(e2,2));
        Double_t errY_jet_Periph = sqrt(pow(e3,2)+pow(e4,2));
        
        cout << "Y_jet_Central = " << Y_jet_Central << " +/- " << errY_jet_Central <<endl;
        cout << "Y_jet_Periph = " << Y_jet_Periph << " +/- " << errY_jet_Periph <<endl;
        
        cout << "Nassoc_Central_TKL = " << Nassoc_Central_TKL<<endl;
        cout << "Nassoc_Periph_TKL = " << Nassoc_Periph_TKL<<endl;
        
        Double_t factorCMS = static_cast<Double_t>(double(Nassoc_Periph_TKL)/Nassoc_Central_TKL)*(Y_jet_Central/Y_jet_Periph);
        Double_t errfactorCMS = factorCMS*sqrt(pow(1/sqrt(Nassoc_Periph_TKL),2) + pow(1/sqrt(Nassoc_Central_TKL),2) + pow(errY_jet_Central/Y_jet_Central,2) + pow(errY_jet_Periph/Y_jet_Periph,2));
        
        cout << "factorCMS = " << factorCMS<<endl;
        
        for(int index=0; index<4; index++){
            V2_CMS_Difference[index] = V2_CMS_Central[index] - V2_CMS_Periph[index]*factorCMS;
            errCompound[index] = V2_CMS_Periph[index]*factorCMS*sqrt(pow(errfactorCMS/factorCMS,2) + pow(errV2_CMS_Periph[index]/V2_CMS_Periph[index],2));
            errV2_CMS_Difference[index] = sqrt(pow(errV2_CMS_Central[index],2)+pow(errCompound[index],2));
            
            cout << "V_CMSDiff_" << index << " = " << V2_CMS_Difference[index] << " +/- " <<  errV2_CMS_Difference[index]<<endl;
        }

    }
    
    
    
    
//        for (int i=0;i<6;i++){
//               sprintf(histoname,Form("bin%d_0",i+1)); // ZZZZZZZ
//               res = FittingAllInvMassBin(histoname, c11b_0, i);
//           }
//    for (int i=0;i<6;i++){
//        sprintf(histoname,Form("bin%d_1",i+1));
//        res = FittingAllInvMassBin(histoname, c11b_1, i);
//    }
//    for (int i=0;i<6;i++){
//        sprintf(histoname,Form("bin%d_2",i+1));
//        res = FittingAllInvMassBin(histoname, c11b_2, i);
//    }
    
    
    
    
    
//    for(int i=0; i<12; i++){
//  //  histoname = "Inv Mass, Phi: " + std::to_string(0) + " pi/6 to " + std::to_string(1) + " pi/6";
//        sprintf(histoname,"InvMassPhiBin_%d",i);
//        FittingAllInvMassPhiBin(histoname, cinvmassphibin, i);
//    }
    std::cout << "===== Dimuons counted =====" <<std::endl;
    std::cout<< "DimuCentralSeen " << DimuCentralSeen <<std::endl;
    std::cout<< "DimuPeriphSeen " << DimuPeriphSeen <<std::endl;
    std::cout<< "DimuSeenMassCut " << DimuSeenMassCut <<std::endl;
    std::cout<< "DimuSeenNoMassCut " << DimuSeenNoMassCut <<std::endl;
    std::cout << "EventRejected: " << EventRejected <<std::endl;
    std::cout << "EventPileUpVtx: " << EventPileUpVtx <<std::endl;
    std::cout << "EventPileUpMult: " << EventPileUpMult <<std::endl;
    std::cout << "EventNC: " << EventNC <<std::endl;
    std::cout << "cmul  " << cmul <<std::endl;
    std::cout <<"countsigma :" << countsigma <<std::endl;
    std::cout <<"BinZeroLeft :" << BinZeroLeft <<std::endl;
}
    









// FITTING INVARIANT MASS METHODS

Double_t FourierV2_WrtInvMass(Double_t *x,Double_t *par)
// Par 0->7: signal, 8->11 Bkg, 12: v2 JPsi, 13->15: V2 bkg
// PSI2STEST { return (JPsiCrystalBallExtended(x,par)*par[12] + (ExpBkg(x,&par[8])+Psi2SCrystalBallExtended(x,par))*(par[13]*(x[0]-mJpsi)*(x[0]-mJpsi) + par[14]*(x[0]-mJpsi) + par[15]))/(JPsiCrystalBallExtended(x,par)+Psi2SCrystalBallExtended(x,par)+ExpBkg(x,&par[8])) ;}
{ return (JPsiCrystalBallExtended(x,par)*par[12] + (ExpBkg(x,&par[8]))*(par[13]*(x[0]-mJpsi)*(x[0]-mJpsi) + par[14]*(x[0]-mJpsi) + par[15]))/(JPsiCrystalBallExtended(x,par)+ExpBkg(x,&par[8])) ;}

Double_t FourierV2_WrtInvMassIntegrated(Double_t *x,Double_t *par)
// Par 0->7: signal, 8->11 Bkg, 12: v2 JPsi, 13->15: V2 bkg
{   Double_t xx = x[0];
    Double_t iBin = V2JPsiTkl->FindBin(xx);
    Double_t binLow = V2JPsiTkl->GetXaxis()->GetBinLowEdge(iBin);
    Double_t binHigh = V2JPsiTkl->GetXaxis()->GetBinUpEdge(iBin);
    
    TF1 *bloug = new TF1("bloug", FourierV2_WrtInvMass, 0,5,16);
    bloug->SetParameters(par);
    Double_t valeur = bloug->Integral(binLow,binHigh)/(binHigh-binLow);

    return valeur;}

Double_t BackFcnV2(Double_t *x,Double_t *par)
//PSI2STEST {return ((ExpBkg(x,&par[8])+Psi2SCrystalBallExtended(x,par))*(par[13]*(x[0]-mJpsi)*(x[0]-mJpsi) + par[14]*(x[0]-mJpsi) + par[15]))/(JPsiCrystalBallExtended(x,par)+Psi2SCrystalBallExtended(x,par)+ExpBkg(x,&par[8])) ;}
{return ((ExpBkg(x,&par[8]))*(par[13]*(x[0]-mJpsi)*(x[0]-mJpsi) + par[14]*(x[0]-mJpsi) + par[15]))/(JPsiCrystalBallExtended(x,par)+ExpBkg(x,&par[8])) ;}

Double_t BackFcnV2Poly(Double_t *x,Double_t *par)
{return (par[13]*(x[0]-mJpsi)*(x[0]-mJpsi) + par[14]*(x[0]-mJpsi) + par[15]) ;}

Double_t SignalFcnJPsiV2(Double_t *x,Double_t *par)
//PSI2STEST {return (JPsiCrystalBallExtended(x,par)*par[12])/(JPsiCrystalBallExtended(x,par)+Psi2SCrystalBallExtended(x,par)+ExpBkg(x,&par[8])) ;}
{return (JPsiCrystalBallExtended(x,par)*par[12])/(JPsiCrystalBallExtended(x,par)+ExpBkg(x,&par[8])) ;}

Double_t FourierV2(Double_t *x,Double_t *par)

{ return par[0] + 2*par[2]*cos(2*x[0]) + 2*par[1]*cos(x[0]);}

Double_t FourierV5(Double_t *x,Double_t *par)

{ return par[0] + 2*par[2]*cos(2*x[0]) + 2*par[1]*cos(x[0]) + 2*par[3]*cos(3*x[0]);}// + 2*par[4]*cos(4*x[0]);}// + 2*par[5]*cos(5*x[0]) + 2*par[6]*cos(6*x[0]) + 2*par[7]*cos(7*x[0]) + 2*par[8]*cos(8*x[0]) + 2*par[9]*cos(9*x[0]) + 2*par[10]*cos(10*x[0]) + 2*par[11]*cos(11*x[0]);}

Double_t TwoCBE2E(Double_t *x,Double_t *par)

//PSI2STEST { return JPsiCrystalBallExtended(x,par) + Psi2SCrystalBallExtended(x,par) + ExpBkg(x,&par[8]);}
{ return JPsiCrystalBallExtended(x,par) + ExpBkg(x,&par[8]);}

Double_t ExpBkg(Double_t *x,Double_t *par)

{ double sigma = par[2] + par[3]*((x[0]-par[1])/par[1]);
return par[0]*exp(-1.0*pow(x[0]-par[1],2)/(2*pow(sigma,2))); }
//{ return par[0]*(exp(x[0]*par[1]*(-1))) + par[2]*(exp(x[0]*par[3]*(-1))); }

Double_t JPsiCrystalBallExtended(Double_t *x,Double_t *par)
{
    
    Double_t sum = 0;
    
    Double_t t = (x[0]-par[1])/par[2];
    if (par[3] < 0) t = -t;
    
    Double_t absAlpha = fabs((Double_t)par[3]);
    Double_t absAlpha2 = fabs((Double_t)par[5]);
    
    
    if (t < -absAlpha) //left tail
    {
        Double_t a =  TMath::Power(par[4]/absAlpha,par[4])*exp(-0.5*absAlpha*absAlpha);
        Double_t b = par[4]/absAlpha - absAlpha;
        
        sum += (par[0]/(par[2]*sqrt(TMath::Pi()*2)))*(a/TMath::Power(b - t, par[4]));
    } else if (t >= -absAlpha && t < absAlpha2) // gaussian core
    {
        sum += (par[0]/(par[2]*sqrt(TMath::Pi()*2)))*(exp(-0.5*t*t));
    } else if (t >= absAlpha2) //right
    {
        Double_t c =  TMath::Power(par[6]/absAlpha2,par[6])*exp(-0.5*absAlpha2*absAlpha2);
        Double_t d = par[6]/absAlpha2 - absAlpha2;
        
        sum += (par[0]/(par[2]*sqrt(TMath::Pi()*2)))*(c/TMath::Power(d + t, par[6]));
    } else
        sum += 0;
    
    return sum ;
}

Double_t Psi2SCrystalBallExtended(Double_t *x,Double_t *par)
{
      Double_t sum = 0;
    
    Double_t t = (x[0]-(par[1]+mDiff))/(par[2]);
  if (par[3] < 0) t = -t;
  
  Double_t absAlpha = fabs((Double_t)par[3]);
  Double_t absAlpha2 = fabs((Double_t)par[5]);
  
  if (t < -absAlpha) //left tail
  {
      Double_t a =  TMath::Power(par[4]/absAlpha,par[4])*exp(-0.5*absAlpha*absAlpha);
      Double_t b = par[4]/absAlpha - absAlpha;
      
      sum += (par[7]/(par[2]*sqrt(2*TMath::Pi())))*(a/TMath::Power(b - t, par[4]));
  }
  else if (t >= -absAlpha && t < absAlpha2) // gaussian core
  {
      sum += (par[7]/(par[2]*sqrt(2*TMath::Pi())))*(exp(-0.5*t*t));
  }
  else if (t >= absAlpha2) //right tail
  {
      Double_t c =  TMath::Power(par[6]/absAlpha2,par[6])*exp(-0.5*absAlpha2*absAlpha2);
      Double_t d = par[6]/absAlpha2 - absAlpha2;
      
      sum += (par[7]/(par[2]*sqrt(2*TMath::Pi())))*(c/TMath::Power(d + t, par[6]));
  } else
      sum += 0;
    
    return sum ;

}
    
    TFitResultPtr FittingAllInvMass(const char *histoname, TCanvas *cinvmass){
        TFitResultPtr res = FittingAllInvMassBin(histoname, cinvmass, 0);
        return res;
    }

TFitResultPtr FittingAllInvMassBin(const char *histoname, TCanvas *cinvmass, int i){
 //Bevington Exercise by Peter Malzacher, modified by Rene Brun, TO BE MODIFIED FOR OUR USE
    // PARAMETERS
    /*
     ------ Signal CB2
     0: Normalisation Gauss
     1: Mean x
     2: Sig Gauss
     3: Alpha 1
     4: Power 1
     5: Alpha 2
     6: Power 2
     7: Normalisation Psi2s, JPsi
     ------ Background 2Exp
     8: Normalisation Before
     9: Exponent Before
     10: Norlmalisation After
     11: Exponent After
     */
    TFile *file0 = new TFile("/Volumes/SEBUSB/InvMass9.root");
    
    cinvmass->cd(i+1);
   cinvmass->SetFillColor(33);
   cinvmass->SetFrameFillColor(41);
   cinvmass->SetGrid();
   TH1F *histo = (TH1F*)file0->Get(histoname);
   // create a TF1 with the range from 0 to 3 and 6 parameters
   TF1 *fitFcn = new TF1("fitFcn",TwoCBE2E,MinInvMassFit,MaxInvMassFit,12); //CVETAN 2.1,5.1->2,4
   fitFcn->SetNpx(500);
   fitFcn->SetLineWidth(4);
   fitFcn->SetLineColor(kMagenta);
   // first try without starting values for the parameters
   // This defaults to 1 for each param.
   // this results in an ok fit for the polynomial function
   // however the non-linear part (lorenzian) does not
   // respond well.
    Double_t params[12] = {10,1,1,1,1,1,1,1,1,1,1,1};
   fitFcn->SetParameters(params);
    TVirtualFitter::Fitter(histo)->SetMaxIterations(100000);
    TVirtualFitter::Fitter(histo)->SetPrecision();
 //  histo->Fit("fitFcn","0");
   // First fit, fix MJPsi and Sigma
    
    fitFcn->SetParLimits(0,0.1,100000000);
    fitFcn->SetParLimits(1,3.05,3.15);
    fitFcn->SetParLimits(2,0.03,0.1);
    fitFcn->SetParLimits(3,0.7,1.1);
    fitFcn->SetParLimits(4,1,30);
    fitFcn->SetParLimits(5,1,5);
    fitFcn->SetParLimits(6,5,25);
    fitFcn->SetParLimits(7,0.001,1000000);
    fitFcn->SetParLimits(8,0.001,10000000);
    fitFcn->SetParLimits(9,0.001,100);
    fitFcn->SetParLimits(10,0.001,100);
    fitFcn->SetParLimits(11,0.001,100);
//    fitFcn->SetParLimits(8,0.1,10000000);
//    fitFcn->SetParLimits(9,-20,20);
//    fitFcn->SetParLimits(10,0.01,100000000);
//    fitFcn->SetParLimits(11,0.01,50);
    
   fitFcn->FixParameter(1,3.096916); // Mean x core
    fitFcn->FixParameter(2,0.07);
//    fitFcn->FixParameter(3,0.883);
//    fitFcn->FixParameter(4,9.940);
//    fitFcn->FixParameter(5,1.832);
//    fitFcn->FixParameter(6,15.323);
    //Tails pPb
        fitFcn->FixParameter(3,1.04);
        fitFcn->FixParameter(4,5.09);
        fitFcn->FixParameter(5,2.02);
        fitFcn->FixParameter(6,9.18);
    
    //PSI2STEST
    fitFcn->FixParameter(7,0);
    fitFcn->SetParameter(8,30000);
    fitFcn->SetParameter(9,2);
    fitFcn->SetParameter(10,0.5);
    fitFcn->SetParameter(11,0.5);
//    fitFcn->FixParameter(8,1);
//    fitFcn->SetParameter(9,0.01);
//    fitFcn->SetParameter(10,100);
//    fitFcn->SetParameter(11,10);

    
    fitFcn->SetParName(0,"Norm_{J/#psi}");
    fitFcn->SetParName(1,"M_{J/#psi}");
    fitFcn->SetParName(2,"Sigma_{J/#psi}");
    fitFcn->SetParName(3,"a_{1}");
    fitFcn->SetParName(4,"n_{1}");
    fitFcn->SetParName(5,"a_{2}");
    fitFcn->SetParName(6,"n_{2}");
    fitFcn->SetParName(7,"Norm_{#Psi(2S)}");
    fitFcn->SetParName(8,"Norm_{Bkg}");
    fitFcn->SetParName(9,"Mean_{Bkg}");
    fitFcn->SetParName(10,"Alpha");
    fitFcn->SetParName(11,"Beta");
//    fitFcn->SetParName(8,"Norm_{TailLowM}");
//    fitFcn->SetParName(9,"Exp_{TailLowM}");
//    fitFcn->SetParName(10,"Norm_{TailHighM}");
//    fitFcn->SetParName(11,"Exp_{TailHighM}");
    
    cout<< "Fit of "<<histoname<<endl;
    
   TFitResultPtr res = histo->Fit("fitFcn","SBMER","ep");
   res = histo->Fit("fitFcn","SBMERQ","ep");
      fitFcn->ReleaseParameter(1);
   res = histo->Fit("fitFcn","SBMERQ","ep");
    fitFcn->ReleaseParameter(2);
   res = histo->Fit("fitFcn","SBMERQ","ep");
   // PSI2STEST fitFcn->ReleaseParameter(7);
   res = histo->Fit("fitFcn","SBMERQ","ep");
   res = histo->Fit("fitFcn","SBMER","ep");
    gStyle->SetOptStat("ne");
    gStyle->SetOptFit(1012);
    TPaveStats *st = (TPaveStats*)histo->FindObject("stats");
    st->SetX1NDC(0.8); //new x start position
    st->SetY1NDC(0.3); //new x end position

//    fitFcn->ReleaseParameter(3);
//    fitFcn->ReleaseParameter(4);
//   res = histo->Fit("fitFcn","SBMER","ep");
//   res = histo->Fit("fitFcn","SBMER","ep");
//    fitFcn->ReleaseParameter(5);
//    fitFcn->ReleaseParameter(6);
//   res = histo->Fit("fitFcn","SBMER","ep");
//    res = histo->Fit("fitFcn","SBMER","ep");

    
   // improve the picture:
   TF1 *backFcn = new TF1("backFcn",ExpBkg,2.1,5.1,4);
   backFcn->SetLineColor(kRed);
   TF1 *signalFcnJPsi = new TF1("signalFcnJPsi",JPsiCrystalBallExtended,2.1,5.1,8);
   TF1 *signalFcnPsi2S = new TF1("signalFcnPsi2S",Psi2SCrystalBallExtended,2.1,5.1,8);
   TPaveText *pave = new TPaveText(0.5,0.5,0.7,0.6,"brNDC");
   signalFcnJPsi->SetLineColor(kBlue);
   signalFcnJPsi->SetNpx(500);
    signalFcnPsi2S->SetLineColor(kGreen);
    signalFcnPsi2S->SetNpx(500);
   Double_t par[12];
   // writes the fit results into the par array
   fitFcn->GetParameters(par);
   signalFcnJPsi->SetParameters(par);
    signalFcnPsi2S->SetParameters(par);
//   auto covMatrix = res->GetCovarianceMatrix();
//   std::cout << "Covariance matrix from the fit ";
//   covMatrix.Print();
   Double_t integral = (signalFcnJPsi->Integral(0,5))*NbinsDimuInvMass/(MaxInvMass-MinInvMass);
    auto covtot = res->GetCovarianceMatrix();
    auto covsgn = covtot.GetSub(0,7,0,7);
    std::cout << "Matrice totale" <<endl;
    covtot.Print();
    std::cout << "Matrice réduite" <<endl;
    covsgn.Print();
    std::cout << "STATUS COV " << res->CovMatrixStatus() <<endl;
   Double_t integralerror = (signalFcnJPsi->IntegralError(0,5,signalFcnJPsi->GetParameters(), res->GetCovarianceMatrix().GetSub(0,7,0,7).GetMatrixArray() ))*NbinsDimuInvMass/(MaxInvMass-MinInvMass);
    std::cout << "Erreur integrale " << integralerror <<endl;
    
    std::cout << "Fitted " << histoname << std::endl;
    std::cout << "Nb JPsi total: " << integral << std::endl;
 //   std::cout << "integral error: " << integralerror << std::endl;
   signalFcnJPsi->Draw("same");
    signalFcnPsi2S->Draw("same");
   backFcn->SetParameters(&par[8]);
   backFcn->Draw("same");
   // draw the legend
    char str[50];
    sprintf(str, "N_{J/#psi} %i +/- %i", int(integral), int(integralerror));
   pave->AddText(str);
    sprintf(str, "M_{#Psi(2S)} = %f, Sig_{#Psi(2S)} = %f", int((mJpsi+mDiff)*1000)/1000., ratSigma*par[2]);
   // pave->AddText(str);
    pave->SetTextFont(42);
    pave->SetTextSize(0.04);
    pave->SetBorderSize(0);
    pave->SetFillStyle(0);
   pave->Draw();
   TLegend *legend=new TLegend(0.4,0.6,0.6,0.8);
   legend->SetTextFont(42);
   legend->SetTextSize(0.03);
    legend->SetFillColorAlpha(kWhite, 0.);
    legend->SetBorderSize(0);
    Char_t message[80];
    sprintf(message,"Fit, #chi^{2}/NDF = %.2f / %d",fitFcn->GetChisquare(),fitFcn->GetNDF());
    legend->AddEntry(fitFcn,message);
    if(res->CovMatrixStatus() == 3){
               // sprintf(message,"The fit is a success");
            }
            else{
                sprintf(message,"The fit is a failure");
                legend->AddEntry(fitFcn,message);
            }
    legend->AddEntry(histo,"Data","lpe");
    legend->AddEntry(backFcn,"Background fit","l");
    legend->AddEntry(signalFcnJPsi,"J/#psi Signal fit","l");
     legend->AddEntry(signalFcnPsi2S,"#Psi(2S) Signal fit","l");
    legend->Draw();
    
//    cout << "Histogramme: " << histoname << endl;
//    for(int j=0; j<10; j++){
//        cout << "Bin de masse numero " << j << " - Signal : " << (signalFcn->Integral(2 + j*0.25,2 + (j+1)*0.25))/0.01 << " et Background : " << (backFcn->Integral(2 + j*0.25,2 + (j+1)*0.25))/0.01 <<endl;
//    }
    
    return res;
}

void FcnCombinedAllMass(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */  )
{
  TAxis *xaxis1  = hnseg->GetXaxis();
  TAxis *xaxis2  = V2JPsiTkl->GetXaxis();

  int nbinX1 = hnseg->GetNbinsX();
  int nbinX2 = V2JPsiTkl->GetNbinsX();

  double chi2 = 0;
  double x[1];
  double tmp;
  npfits = 0;
  for (int ix = 1; ix <= nbinX1; ++ix) {
    x[0] = xaxis1->GetBinCenter(ix);
      if ( hnseg->GetBinError(ix) > 0 ) {
        tmp = (hnseg->GetBinContent(ix) - TwoCBE2E(x,p))/hnseg->GetBinError(ix);
        chi2 += tmp*tmp;
        npfits++;
      }
  }
  for (int ix = 1; ix <= nbinX2; ++ix) {
     x[0] = xaxis2->GetBinCenter(ix);
      if ( V2JPsiTkl->GetBinError(ix) > 0 ) {
        tmp = (V2JPsiTkl->GetBinContent(ix) - FourierV2_WrtInvMass(x,p))/V2JPsiTkl->GetBinError(ix);
        chi2 += tmp*tmp;
        npfits++;
      }
  }
  fval = chi2;
}

int GetCent(float cent){
    if(cent < 1){
        return 0;
    }
    else if(cent < 3){
        return 1;
    }
    else if(cent < 5){
        return 2;
    }
    else if(cent < 10){
        return 3;
    }
    else if(cent < 15){
        return 4;
    }
    else if(cent < 20){
        return 5;
    }
    else if(cent < 30){
        return 6;
    }
    else if(cent < 40){
        return 7;
    }
    else if(cent < 50){
        return 8;
    }
    else if(cent < 60){
        return 9;
    }
    else if(cent < 70){
        return 10;
    }
    else if(cent < 80){
        return 11;
    }
    else if(cent < 90){
        return 12;
    }
    else{
        return 13;
    }
    
//    if(cent<=20){
//        return 0;
//    }
//    else if(cent<=40){
//        return 1;
//    }
//    else{
//        return 2;
 //   }
}

int GetCentPM(int Ntkl, float SPDTrackletsPer, float SPDClustersPer, float V0MPer, int zvtx_idx, int groupnumber){
    
    float estimator = 0;

    if(centralityMethod == "SPDTrackletsPercentile"){
        estimator = SPDTrackletsPer;
    }
    else if(centralityMethod == "SPDClustersPercentile"){
        estimator = SPDClustersPer;
    }
    else if(centralityMethod == "V0MPercentile"){
        estimator = V0MPer;
    }
    else{
        return -1;
    }
    return GetCent(estimator);
    
}

int GetCentCvetan(float V0MPer){
    
    float estimator = 0;

  
    if(centralityMethod == "V0MPercentile"){
        estimator = V0MPer;
    }
    else{
        return -1;
    }
    return GetCent(estimator);
    
}

int GetPtBin(double pt){

    for(int pt_index=1;pt_index<NbPtBins+1;pt_index++){
        if(PtBins[pt_index] > pt){
            return pt_index-1;
        }
    }
}

int GetMassBin(double mass){

    for(int mass_index=1;mass_index<NbinsInvMass+1;mass_index++){
        if(MassBins[mass_index] > mass){
            return mass_index-1;
        }
    }
}

Double_t myEtaLow(double x){
    return TMath::Log(TMath::Tan(0.5*TMath::ATan(7.6/(14.1+(floor(x)+0.5))))) + 0.1;}

Double_t myEtaHigh(double x){
return -(TMath::Log(TMath::Tan(0.5*TMath::ATan(7.6/(14.1-(floor(x)+0.5)))))) - 0.05;}

Double_t fitPLL(Double_t *x, Double_t *par)
{
  Double_t xx = x[0];
  Double_t iBin = hStatic3->FindBin(xx);

  return par[0]*hStatic3->GetBinContent(iBin);
}

Double_t fitLin(Double_t *x, Double_t *par)
{
    Double_t DeltaEtaTKLCut = 1.5;
  Double_t xx = x[0];
    if(abs(xx)<DeltaEtaTKLCut){
        TF1::RejectPoint();
        return 0;
    }

    return par[0]+(par[1]*abs(xx));
}