#include <iostream>
#include <fstream>

gROOT->Reset();

// Setup file including parameters
#include "histMk_3D_setup.h"

void histMk_3D_B0B0_Dlnu_2nd_DchKPiPi() {

  // Root macro to produce 3D histogram of B0B0 MC.
  // B0B0 MC is luminosity nomalized to OnPeak data.
  //###########################################################
  // Change these values to meet your condition
  //###########################################################
 
  TString RunNum("RunAll");
  if (run == 1) RunNum = "Run1";
  if (run == 2) RunNum = "Run2";
  if (run == 3) RunNum = "Run3";
  if (run == 4) RunNum = "Run4";

  TString MCCorr("noCorr");
  if (correctTrk) MCCorr ="effCorr";

  // Name base of output files
  TString outFileNameBase("histo_3D_"+RunNum+"_"+MCCorr);

  if (newSlope) {
    TString Slope("New");
  } else {
    TString Slope("Old");
  }

  // Info of output files
  TString out3DHistoName(outFileNameBase+"_B0B0_Dlnu_2nd"+Slope+"_DchKPiPi.root");
  TString outInfoFileName(outFileNameBase+"_B0B0_Dlnu_2nd"+Slope+"_DchKPiPi.txt");
  TString histTitle_peak("B0B0 to Dlnu (DchKPipi) D mass peak");
  TString histTitle_side("B0B0 to Dlnu (DchKPipi) D mass sideband");
  TString histTitle_sub("B0B0 to Dlnu (DchKPipi) sideband subtracted");

  // B decay mode selection
  int selBDecMod1 = 3120;
  int selBDecMod2 = 4110;
  int selBDecMod3 = -1;
  int selBDecMod4 = -1;

  // D decay mode selection
  int DDecMod_3D = DDecMod_DchKPiPi;

  // D mass range
  float DMass1 = DMass_DchKPiPi1; // lower edge of sideband
  float DMass2 = DMass_DchKPiPi2; // lower edge of D mass peak
  float DMass3 = DMass_DchKPiPi3; // upper edge of D mass peak
  float DMass4 = DMass_DchKPiPi4; // upper edge of sideband

  // Input files
  TString inFileName_Run1(baseDir_Run1+"/"+inFile_Run1_B0B0);
  int nStart_Run1 = nStart_Run1_B0B0;
  int nStop_Run1 = nStop_Run1_B0B0;
  TString inFileName_Run2(baseDir_Run2+"/"+inFile_Run2_B0B0);
  int nStart_Run2 = nStart_Run2_B0B0;
  int nStop_Run2 = nStop_Run2_B0B0;
  TString inFileName_Run3(baseDir_Run3+"/"+inFile_Run3_B0B0);
  int nStart_Run3 = nStart_Run3_B0B0;
  int nStop_Run3 = nStop_Run3_B0B0;
  TString inFileName_Run4(baseDir_Run4+"/"+inFile_Run4_B0B0);
  int nStart_Run4 = nStart_Run4_B0B0;
  int nStop_Run4 = nStop_Run4_B0B0;

  // Luminosity normalization factors
  double c_Run1 = c_Run1_B0B0;
  double c_Run2 = c_Run2_B0B0;
  double c_Run3 = c_Run3_B0B0;
  double c_Run4 = c_Run4_B0B0;
  double c_All = c_All_B0B0;


  //#######################################################################
  // Create D mass peak and sideband arrays, readin data and fill the arrays.

  cout << "starat B0B0ToDlnu\n";
  cout << "start chain\n";
  if (run == 1) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run1; i<=nStop_Run1; i++) {
      chain_3D->Add(inFileName_Run1+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run1;
  } else if (run == 2) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run2; i<=nStop_Run2; i++) {
      chain_3D->Add(inFileName_Run2+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run2;
  } else if (run == 3) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run3; i<=nStop_Run3; i++) {
      chain_3D->Add(inFileName_Run3+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run3;
  } else if (run == 4) {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run4; i<=nStop_Run4; i++) {
      chain_3D->Add(inFileName_Run4+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_Run4;
  } else {
    TChain * chain_3D = new TChain(ntpName);
    for (Int_t i=nStart_Run1; i<=nStop_Run1; i++) {
      chain_3D->Add(inFileName_Run1+"-"+i+".root");
    }
    for (Int_t i=nStart_Run2; i<=nStop_Run2; i++) {
      chain_3D->Add(inFileName_Run2+"-"+i+".root");
    }
    for (Int_t i=nStart_Run3; i<=nStop_Run3; i++) {
      chain_3D->Add(inFileName_Run3+"-"+i+".root");
    }
    for (Int_t i=nStart_Run4; i<=nStop_Run4; i++) {
      chain_3D->Add(inFileName_Run4+"-"+i+".root");
    }
    TTree * tree_3D = chain_3D;
    double conv = c_All;
  }





  // Define and initialise arrays for 3D histo
  double num_peak[nbinsx][nbinsy][nbinsz];
  double num_side[nbinsx][nbinsy][nbinsz];
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	num_peak[i][j][k] = 0;
	num_side[i][j][k] = 0;
      }
    }
  }

  cout << "read tree\n";
  // variables to read out
  Int_t nDl_3D(0);
  Int_t DDecayMode_3D[200];
  Float_t DMass_3D[200];
  Int_t noExtraKorPi_3D[200], mcKaonMissID_3D[200];
  Int_t mcPiMissID_3D[200], mcKPifromSameD_3D[200];
  Int_t mcLepMissID_3D[200], mcDirectLepFromB_3D[200];
  Int_t trueLepFromSameB_3D[200];
  Int_t trueBDecMod_3D[200];
  double trkEffCorr_3D[200], PIDCorr_3D[200];
  double mcW_3D[200];
  Float_t varx_3D[200], vary_3D[200], varz_3D[200];
  
  // Set reference between the above variables and root-tuple entries
  tree_3D->SetBranchAddress("nDl", &nDl_3D);
  tree_3D->SetBranchAddress("DDecayMode", &DDecayMode_3D);
  tree_3D->SetBranchAddress("DMass", &DMass_3D);
  tree_3D->SetBranchAddress("noExtraKorPi", &noExtraKorPi_3D);
  tree_3D->SetBranchAddress("mcKaonMissID", &mcKaonMissID_3D);
  tree_3D->SetBranchAddress("mcPiMissID", &mcPiMissID_3D);
  tree_3D->SetBranchAddress("mcKPifromSameD", &mcKPifromSameD_3D);
  tree_3D->SetBranchAddress("mcLepMissID", &mcLepMissID_3D);
  tree_3D->SetBranchAddress("mcDirectLepFromB", &mcDirectLepFromB_3D);
  tree_3D->SetBranchAddress("trueBDecMod", &trueBDecMod_3D);
  tree_3D->SetBranchAddress("trueLepFromSameB", &trueLepFromSameB_3D);
  tree_3D->SetBranchAddress("trkEffCorr", &trkEffCorr_3D);
  tree_3D->SetBranchAddress("PIDCorr", &PIDCorr_3D);
  tree_3D->SetBranchAddress("mcW", &mcW_3D);
  tree_3D->SetBranchAddress(varx, &varx_3D);
  tree_3D->SetBranchAddress(vary, &vary_3D);
  tree_3D->SetBranchAddress(varz, &varz_3D);

  // Read out root-tuple entries and fill arrays
  Int_t nevent_3D =(Int_t)tree_3D->GetEntries(); // number of events
  cout << "number of 3D events = " << nevent_3D << endl;
  int numCand(0), numSelCand(0);
  int numCand_DDecMod(0), numCand_Index(0);
  for (Int_t i=0;i<nevent_3D;i++) { // loop over events
    if (i==1) cout << "3D event 1" << endl;
    if (i==100) cout << "3D event 100" << endl;
    if (i==1000) cout << "3D event 1000" << endl;
    if (i==10000) cout << "3D event 10000" << endl;
    if (i==100000) cout << "3D event 100000" << endl;
    if (i==1000000) cout << "3D event 1000000" << endl;
    if (i==2000000) cout << "3D event 2000000" << endl;
    if (i==3000000) cout << "3D event 3000000" << endl;
    if (i==4000000) cout << "3D event 4000000" << endl;
    
    tree_3D->GetEntry(i); //read complete event in memory

    for (Int_t j=0;j<nDl_3D;j++) { // loop over Dl candidates
      int DDecMod = DDecayMode_3D[j];
      float DhadMass = DMass_3D[j];
      int noExtKorPi = noExtraKorPi_3D[j];
      int kaonMissID = mcKaonMissID_3D[j];
      int piMissID = mcPiMissID_3D[j];
      int KPiFromSameD = mcKPifromSameD_3D[j];
      int lepMissID = mcLepMissID_3D[j];
      int directLepFromB = mcDirectLepFromB_3D[j];
      int lepFromSameB = trueLepFromSameB_3D[j];
      int mcBDecMod = trueBDecMod_3D[j];
      double trkCorr = trkEffCorr_3D[j] * trkCorrModif;
      double PidCorr = PIDCorr_3D[j] * PidCorrModif;
      double trueW = mcW_3D[j];
      float xvar = varx_3D[j];
      float yvar = vary_3D[j];
      float zvar = varz_3D[j];

      numCand++;
      if (DDecMod != DDecMod_3D) continue;
      numCand_DDecMod++;
      if (noExtKorPi == 0) continue;
      if (kaonMissID == 1) continue;
      if (piMissID == 1) continue;
      if (KPiFromSameD == 0) continue;
      if (lepMissID == 1) continue;
      if (directLepFromB == 0) continue;
      if (lepFromSameB == 0) continue;
      int goodBDecMod(0);
      if (mcBDecMod == selBDecMod1 )  goodBDecMod=1;
      if (mcBDecMod == selBDecMod2 )  goodBDecMod=1;
      if (mcBDecMod == selBDecMod3 )  goodBDecMod=1;
      if (mcBDecMod == selBDecMod4 )  goodBDecMod=1;
      if (goodBDecMod != 1) continue;

      // Get indeces of this candidate
      //cout << "getting indecies" << endl;
      int xIndex(-1);
      for (Int_t k=0;k<nbinsx;k++) {
	if (xbins[k]<  xvar &&  xvar <= xbins[k+1]) {
	  xIndex = k;
	}
      }
      int yIndex(-1);
      for (Int_t k=0;k<nbinsy;k++) {
	if (ybins[k]< yvar && yvar <= ybins[k+1]) {
	  yIndex = k;
	}
      }
      int zIndex(-1);
      for (Int_t k=0;k<nbinsz;k++) {
	if (zbins[k]< zvar && zvar <= zbins[k+1]) {
	  zIndex = k;
	}
      }

      numCand_Index++;

      if (xIndex>=0 && yIndex>=0 && zIndex>=0) {
	numSelCand++;

	// Get FF correction factor
	// Get q-square
	double mB = 5.2790;
	double mD = 1.8646;
	double q2 = mB*mB + mD*mD - 2*mB*mD*trueW;
	//cout << "q2 = " << q2 << endl;

	// Get Fplus
	double mass=mD;
	double mtb, mbb;
	double msd, mx,mb,nf,nfp; 
	double msq,bx2,mbx,mtx;
	double zji,cji,w,gammaji,chiji,betaji_fppfm;
	double rfppfm,rfpmfm,f3fppfm,f3fpmfm,fppfm,fpmfm,al,ai,rcji,f3; 
	double mqm,msb,bb2,mup,mum,bbx2,tm,wt,r2,betaji_fpmfm;
	
	double t=q2;
	
	msb=5.2;
	msd=0.33;
	bb2=0.431*0.431;
	mbb=5.31;
	nf = 4.0;
  
	msq=1.82;
	bx2=0.45*0.45;
	mbx=0.75*2.01+0.25*1.87;
	nfp = 3.0;
	
	mtb = msb + msd;
	mtx = msq + msd;
	mb=5.279;
	mx=mass;
  
	mup=1.0/(1.0/msq+1.0/msb);
	mum=1.0/(1.0/msq-1.0/msb);
	bbx2=0.5*(bb2+bx2);
	tm=(mb-mx)*(mb-mx);
	if ( t>tm ) t=0.99*tm;
	wt=1.0+(tm-t)/(2.0*mbb*mbx);
	
	mqm = 0.1;
  
	double PI = 3.1415;
	double lqcd2 = 0.04;
	double nflav = 4;
	double As0=0.6;
	if ( mqm > 0.6 ) {
	  if ( mqm < 1.85 ) {
	    nflav = 3.0;}
	  
	  As0 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( mqm*mqm/lqcd2);
	}
	
	double As2=0.6;
	if ( msq > 0.6 ) {
	  if ( msq < 1.85 ) {
	    nflav = 3.0;}
	  
	  As2 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( msq*msq/lqcd2);
	}
	r2=3.0/(4.0*msb*msq)+3*msd*msd/(2*mbb*mbx*bbx2) + 
	  (16.0/(mbb*mbx*(33.0-2.0*nfp)))*
	  log(As0/As2);
  
	f3 = sqrt(mtx/mtb)*pow(sqrt(bx2*bb2)/bbx2,1.5) /
	  ((1.0+r2*(tm-t)/12.0)*(1.0+r2*(tm-t)/12.0));
  
	w = 1.0 + (( tm - t ) / ( 2.0* mb * mx ));
	rcji = ( 1/sqrt(w*w -1 ))*log( w + sqrt( w*w -1 ));
	al = (8.0 / ( 33.0 - 2.0*nfp ))*(w*rcji -1.0 );
	ai = -1.0* ( 6.0/( 33.0 - 2.0*nf));  
	
	double As1=0.6;
	if ( msb > 0.6 ) {
	  if ( msb < 1.85 ) {
	    nflav = 3.0;}
	  
	  As1 = 12.0*PI / ( 33.0 - 2.0*nflav) / log( msb*msb/lqcd2);
	}

  
	cji = pow(( As1 / As2 ),ai);
  
	zji = msq / msb;
  
	gammaji = 2+((2.0*zji)/(1-zji))*log(zji);
	gammaji = -1.0*gammaji;
  
	chiji = -1.0 - ( gammaji / ( 1- zji ));
	betaji_fppfm = gammaji - (2.0/3.0)*chiji;
	betaji_fpmfm = gammaji + (2.0/3.0)*chiji;
  
	double As3=0.6;
	double m3=sqrt(msb*msq);
  
	if ( m3 > 0.6 ) {
	  if ( msq < 1.85 ) {
	    nflav = 3.0;}
	  
	  As3 = 12.0*PI / ( 33.0 - 2.0*nflav) /
	    log( m3*m3/lqcd2);
	}
	rfppfm = cji *(1.0 + betaji_fppfm*As3/PI);
  
	rfpmfm = cji *(1.0 + betaji_fpmfm*As3/PI);
  
	f3fppfm = f3*pow(( mbb / mtb ),-0.5)*pow((mbx/mtx),0.5);
	f3fpmfm = f3*pow(( mbb / mtb ),0.5)*pow((mbx/mtx),-0.5);
	fppfm = f3fppfm* rfppfm * ( 2.0 - 
				    ( ( mtx/msq)*(1- ( (msd*msq*bb2)
						       /(2.0*mup*mtx*bbx2)))));
	fpmfm = f3fpmfm* rfpmfm * ( mtb/msq) * ( 1 - 
						 ( ( msd*msq*bb2)/
						   ( 2.0*mup*mtx*bbx2)));
  
	double fpf = (fppfm + fpmfm)/2.0;
	//cout << "fpf = " << fpf << endl;

	// Get R
	double mcR = 2*sqrt(mB*mD) / (mB+mD);

	// Get oldFF
	double oldFF = mcR * fpf;
	double oldFF2 = oldFF * oldFF;
	//cout << "oldFF = " << oldFF << endl;

	// Get newFF
	if (newSlope) {
	  double trueZ 
	    = (sqrt(trueW+1) - sqrt(2))/(sqrt(trueW+1) + sqrt(2));
	  double N01_new = 1 - 10*trueZ*trueZ + 84*trueZ*trueZ*trueZ;
	  double N02_new = -8*trueZ + 51*trueZ*trueZ - 252*trueZ*trueZ*trueZ;
	  double newFF2 =  N02_new *  N02_new;
	} else {
	  double newFF2 = (trueW - 1) * (trueW - 1);
	}

	// Get weight
	double W_D = newFF2 / oldFF2; 

	if (DMass2<DhadMass && DhadMass<DMass3) {
	  if (correctTrk) {
	    num_peak[xIndex][yIndex][zIndex] 
	      += (trkCorr * PidCorr)*W_D;
	  } else {
	    num_peak[xIndex][yIndex][zIndex] += W_D;
	  }
	}
	if ((DMass1<DhadMass && DhadMass<DMass2) 
	    ||(DMass3<DhadMass && DhadMass<DMass4)){
	  if (correctTrk) {
	    num_side[xIndex][yIndex][zIndex] 
	      += (trkCorr * PidCorr)*W_D;
	  } else {
	    num_side[xIndex][yIndex][zIndex] += W_D;
	  }
	}
      }
      
    } //end of Dl loop
  } //end of event loop
  cout << "number of processed candidates = " << numCand << endl;
  cout << "number of candidates after D decay mode selection = " 
       << numCand_DDecMod << endl;
  cout << "number of candidates before index selection = " 
       << numCand_Index << endl;
  cout << "number of selected candidates = " << numSelCand << endl;
  

  //###################################

  // Define and fill histos
  TH3D *peak 
    = new TH3D("peak", histTitle_peak,
	       nbinsx, xbins, nbinsy, ybins, nbinsz, zbins);
  TH3D *side 
    = new TH3D("side", histTitle_side, 
	       nbinsx, xbins, nbinsy, ybins, nbinsz, zbins);
  TH3D *sub 
    = new TH3D("sub", histTitle_sub, 
	       nbinsx, xbins, nbinsy, ybins, nbinsz, zbins);
 
  // Now fill the histograms
  cout << "filling histograms\n";
  for (int i=0; i<nbinsx; i++) {
    for (int j=0; j<nbinsy; j++) {
      for (int k=0; k<nbinsz; k++) {
	peak->SetBinContent(i+1,j+1,k+1,num_peak[i][j][k]);
        side->SetBinContent(i+1,j+1,k+1,num_side[i][j][k]);

	double peak3D = conv * num_peak[i][j][k];
	double side3D = conv * num_side[i][j][k];
	double sub3D = peak3D - side3D;
	double var3D = num_peak[i][j][k] + num_side[i][j][k];
	double err3D = conv * sqrt(var3D);
        sub->SetBinContent(i+1,j+1,k+1,sub3D);
        sub->SetBinError(i+1,j+1,k+1,err3D);
      }
    }
  }

  //#############################################################
  // Open output file and write histograms in it. ######
  cout <<"saving histograms to " << out3DHistoName << endl;
  TFile *f = new TFile(out3DHistoName, "Recreate");

  peak->Write();
  side->Write();
  sub->Write();

  f->Write();
  f->Close();

  //#######################################################################
  // Open output file and save parameters. #####
  cout << "dumping parameters to " << outInfoFileName << endl;
  ofstream outFile(outInfoFileName); //open output file
  if (!outFile) {
    cout << "Cannot open file " << outFile << endl;
  }
  outFile.setf(ios::fixed);
  outFile.setf(ios::showpoint);
  outFile << "Input files" << endl;
  outFile << "  run = " << run << endl;
  if (run == 1) {
    outFile << "  File name base = " << inFileName_Run1 << endl;
    outFile << "  start number = " << nStart_Run1 << endl;
    outFile << "  stop number = " << nStop_Run1 << endl;
  } else if (run == 2) {
    outFile << "  File name base = " << inFileName_Run2 << endl;
    outFile << "  start number = " << nStart_Run2 << endl;
    outFile << "  stop number = " << nStop_Run2 << endl;
  } else if (run == 3) {
    outFile << "  File name base = " << inFileName_Run3 << endl;
    outFile << "  start number = " << nStart_Run3 << endl;
    outFile << "  stop number = " << nStop_Run3 << endl;
  } else if (run == 4) {
    outFile << "  File name base = " << inFileName_Run4 << endl;
    outFile << "  start number = " << nStart_Run4 << endl;
    outFile << "  stop number = " << nStop_Run4 << endl;
  } else {
    outFile << "  File name base = " << inFileName_Run1 << endl;
    outFile << "  start number = " << nStart_Run1 << endl;
    outFile << "  stop number = " << nStop_Run1 << endl;
    outFile << "  File name base = " << inFileName_Run2 << endl;
    outFile << "  start number = " << nStart_Run2 << endl;
    outFile << "  stop number = " << nStop_Run2 << endl;
    outFile << "  File name base = " << inFileName_Run3 << endl;
    outFile << "  start number = " << nStart_Run3 << endl;
    outFile << "  stop number = " << nStop_Run3 << endl;
    outFile << "  File name base = " << inFileName_Run4 << endl;
    outFile << "  start number = " << nStart_Run4 << endl;
    outFile << "  stop number = " << nStop_Run4 << endl;
  }
  outFile << "################" << endl;
  outFile << "Output files" << endl;
  outFile << "  histograms : " << out3DHistoName << endl;
  outFile << "  information : " << outInfoFileName << endl;
  outFile << "################" << endl;
  outFile << "Event number info" << endl;
  outFile << "  number of events processed = " << nevent_3D << endl;
  outFile << "  number of candidates processed = " << numCand << endl;
  outFile << "  number of candidates after D decay mode selection = " 
	  << numCand_DDecMod << endl;
  outFile << "  number of candidates before index selection = " 
	  << numCand_Index << endl;
  outFile << "  number of candidates selected = " << numSelCand << endl;
  outFile << "################" << endl;
  outFile << "Luminosity normalization factor" << endl;
  outFile << "  conv = " << conv << endl;
  outFile << "B decay mode selection" << endl;
  outFile << "  selBDecMod1 = " << selBDecMod1 << endl;
  outFile << "  selBDecMod2 = " << selBDecMod2 << endl;
  outFile << "  selBDecMod3 = " << selBDecMod3 << endl;
  outFile << "  selBDecMod4 = " << selBDecMod4 << endl;
  outFile << "D decay mode selection" << endl;
  outFile << "  DDecMod_3D = " << DDecMod_3D << endl;
  outFile << "D mass cut" << endl;
  outFile << "  DMass1 = " << DMass1 << endl;
  outFile << "  DMass2 = " << DMass2 << endl;
  outFile << "  DMass3 = " << DMass3 << endl;
  outFile << "  DMass4 = " << DMass4 << endl;
  outFile << "3 variables" << endl;
  outFile << "  Variable1 = " << varx << endl;
  outFile << "  Variable2 = " << vary << endl;
  outFile << "  Variable3 = " << varz << endl;
  outFile << "################" << endl;
  outFile << "Do or do not tracking and PID correction" << endl;
  outFile << "  correctTrk = " << correctTrk << endl;
  outFile << "track and PID correction modificatin factors" << endl;
  outFile << "modified correction = Corr * Modif" << endl;
  outFile << "  trkCorrModif = " << trkCorrModif << endl;
  outFile << "  PidCorrModif = " << PidCorrModif << endl;
  outFile << "FF weight W = [1 - rho_D_sq*(w-1)]^2" << endl;
  outFile << "  slope rho_D_sq = " << rho_D_sq << endl;
  outFile.close();

}