#include <iostream>
#include <fstream>

gROOT->Reset();

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

void histMk_3D_B0B0_Dstar_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);

  TString DstarFFMode("oldFF");
  if (newDstarFF) DstarFFMode = "newFF"+DlnuFFMod;

  // FF parameters
  // SP5/SP6
  double R1_old = 1.18;
  double R2_old = 0.72;
  double rhoSq_old = 0.92;

  // mass
  double m_Dstarp = 2010.0;
  double m_Dstar0 = 2006.7;
  double m_Bp = 5279.0;
  double m_B0 = 5279.4;

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

  // B decay mode selection
  int selBDecMod1 = 3220;
  int selBDecMod2 = 4210;
  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;


  //###########################################################
  // Get normalization factor
  
  TH1D *y1H = new TH1D("y1H", "y1H", 100, -1, 3);
  TH1D *y2H = new TH1D("y2H", "y2H", 100, -1, 3);
  TH1D *cosLH = new TH1D("cosLH", "cosLH", 100, -1.1, 1.1);
  TH1D *cosVH = new TH1D("cosVH", "cosVH", 100, -1.1, 1.1);
  TH1D *ChiH = new TH1D("ChiH", "ChiH", 100, -0.1, 6.3);
  TH1D *cosChiH = new TH1D("cosChiH", "cosChiH", 100, -1.1, 1.1);
  TH1D *wH = new TH1D("wH", "wH", 100, 0.5, 2);
  TH1D *FacH = new TH1D("FacH", "FacH", 100, -0.1, 7);
  TH1D *oldFFH = new TH1D("oldFFH", "oldFFH", 100, -0.1, 2);
  TH1D *newFFH = new TH1D("newFFH", "newFFH", 100, -0.1, 2);
  TH1D *oldIntH = new TH1D("oldIntH", "oldIntH", 100, -0.1, 2);
  TH1D *newIntH = new TH1D("newIntH", "newIntH", 100, -0.1, 2);


  TRandom *rdm1 = new TRandom();
  TRandom *rdm2 = new TRandom2();
  TRandom3 *rdm3 = new TRandom();
  TRandom2 *rdm4 = new TRandom();
  TRandom3 *rdm5 = new TRandom3();
  TRandom2 *rdm6 = new TRandom();
  TRandom3 *rdm7 = new TRandom3();
  int N = 1000000;
  int N_old = 0;
  int N_new = 0;
  double RN = 1;
  for (Int_t i=0;i<N;i++) {
    double y1 = (rdm1->Rndm(i))*2; 
    y1H->Fill(y1);
    double y2 = (rdm1->Rndm(i))*2; 
    y2H->Fill(y2);
    double cosL = (rdm3->Rndm(i))*2 - 1; 
    cosLH->Fill(cosL);
    double cosV = (rdm4->Rndm(i))*2 - 1; 
    cosVH->Fill(cosV);
    double Chi = (rdm5->Rndm(i))*2*3.1415; 
    ChiH->Fill(Chi);
    double cosChi = cos(Chi); 
    cosChiH->Fill(cosChi);
    double trueW = (rdm6->Rndm(i))*0.504 + 1;
    wH->Fill(trueW);
    // Get FF correction factor
    double r = m_Dstar0 / m_Bp;
    double hp1 = sqrt(1-2*trueW*r+r*r);
    double hp2 = sqrt((trueW-1)/(trueW+1));
    double F1 = 1 - cosL;
    double F2 = 1 + cosL;
    double F3 = 1 - cosL*cosL;
    double F4 = sqrt(1 - cosL*cosL);
    double F5 = 1 - cosV*cosV;
    double F6 = sqrt(1 - cosV*cosV);
    double F7 = 2 * cosChi*cosChi - 1;
    double F11 = sqrt(trueW - 1);
    double F12 = (trueW + 1)*(trueW + 1);
    FacH->Fill(F11*F12);

    double hp_old = hp1*(1 - hp2*R1_old);
    double hm_old = hp1*(1 + hp2*R1_old);
    double h0_old = 1 - r + (trueW - 1)*(1 - R2_old);
    double N0_old = 1 - rhoSq_old*(trueW -1);
    double N1_old = F1*F1 * F5 * hp_old*hp_old;
    double N2_old = F2*F2 * F5 * hm_old*hm_old;
    double N3_old = 4*F3 * cosV*cosV * h0_old*h0_old;
    double N4_old = 2*F3 * F5 * F7 * hp_old*hm_old;
    double N5_old = 4*F4 * F1 * F6 * cosV * cosChi * hp_old*h0_old;
    double N6_old = 4*F4 * F2 * F6 * cosV * cosChi * hm_old*h0_old;
    double N7_old 
      = N1_old + N2_old + N3_old - N4_old - N5_old + N6_old;
    double FF_old = N0_old * N0_old * N7_old;
    oldFFH->Fill(FF_old);
    double Int_old = F11*F12*FF_old;
    oldIntH->Fill(Int_old);
    if (y1 <= Int_old) N_old++;
 
    double hp_new = hp1*(1 - hp2*R1_new);
    double hm_new = hp1*(1 + hp2*R1_new);
    double h0_new = 1 - r + (trueW - 1)*(1 - R2_new);
    double trueZ 
      = (sqrt(trueW+1) - sqrt(2))/(sqrt(trueW+1) + sqrt(2));
    double N01_new = 8*rhoSq_new*trueZ;
    double N02_new = (53*rhoSq_new - 15)*trueZ*trueZ;
    double N03_new = (231*rhoSq_new - 91)*trueZ*trueZ*trueZ;
    double N0_new = 1 - N01_new + N02_new - N03_new;
    double N1_new = F1*F1 * F5 * hp_new*hp_new;
    double N2_new = F2*F2 * F5 * hm_new*hm_new;
    double N3_new = 4*F3 * cosV*cosV * h0_new*h0_new;
    double N4_new = 2*F3 * F5 * F7 * hp_new*hm_new;
    double N5_new = 4*F4 * F1 * F6 * cosV * cosChi * hp_new*h0_new;
    double N6_new = 4*F4 * F2 * F6 * cosV * cosChi * hm_new*h0_new;
    double N7_new 
      = N1_new + N2_new + N3_new - N4_new - N5_new + N6_new;
    //double FF_new = FF_old;
    double FF_new = N0_new * N0_new * N7_new;
    newFFH->Fill(FF_new);
    double Int_new = F11*F12*FF_new;
    newIntH->Fill(Int_new);
    if (y2 <= Int_new) N_new++;
  }

  if (N_new > 0) {
    RN = (double) N_old / N_new;
  }
  double RN_err = 0;
  if (N_new > 0 && N_old > 0) {
    RN_err = sqrt((double)(N_old + N_new)/(N_old*N_new));
  }
  cout << "N_old = " << N_old << endl;
  cout << "N_new = " << N_new << endl;
  cout << "RN = " << RN << endl;
  cout << "RN_err = " << RN_err << endl;



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

  cout << "starat B0B0ToDstar\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];
  double MCcosL_3D[200], MCcosV_3D[200], MCcosChi_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("MCcosL", &MCcosL_3D);
  tree_3D->SetBranchAddress("MCcosV", &MCcosV_3D);
  tree_3D->SetBranchAddress("MCcosChi", &MCcosChi_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];
      double cosL = MCcosL_3D[j];
      double cosV = MCcosV_3D[j];
      double cosChi = MCcosChi_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
        double r = m_Dstar0 / m_Bp;
	double hp1 = sqrt(1-2*trueW*r+r*r);
	double hp2 = sqrt((trueW-1)/(trueW+1));
	double F1 = 1 - cosL;
	double F2 = 1 + cosL;
	double F3 = 1 - cosL*cosL;
	double F4 = sqrt(1 - cosL*cosL);
	double F5 = 1 - cosV*cosV;
	double F6 = sqrt(1 - cosV*cosV);
	double F7 = 2 * cosChi*cosChi - 1;
	double hp_old = hp1*(1 - hp2*R1_old);
	double hm_old = hp1*(1 + hp2*R1_old);
	double h0_old = 1 - r + (trueW - 1)*(1 - R2_old);
        double N0_old = 1 - rhoSq_old*(trueW -1);
	double N1_old = F1*F1 * F5 * hp_old*hp_old;
        double N2_old = F2*F2 * F5 * hm_old*hm_old;
	double N3_old = 4*F3 * cosV*cosV * h0_old*h0_old;
	double N4_old = 2*F3 * F5 * F7 * hp_old*hm_old;
        double N5_old = 4*F4 * F1 * F6 * cosV * cosChi * hp_old*h0_old;
        double N6_old = 4*F4 * F2 * F6 * cosV * cosChi * hm_old*h0_old;
        double N7_old 
	  = N1_old + N2_old + N3_old - N4_old - N5_old + N6_old;
        double FF_old = N0_old * N0_old * N7_old;
	double hp_new = hp1*(1 - hp2*R1_new);
	double hm_new = hp1*(1 + hp2*R1_new);
	double h0_new = 1 - r + (trueW - 1)*(1 - R2_new);
	double trueZ 
	  = (sqrt(trueW+1) - sqrt(2))/(sqrt(trueW+1) + sqrt(2));
	double N01_new = 8*rhoSq_new*trueZ;
	double N02_new = (53*rhoSq_new - 15)*trueZ*trueZ;
	double N03_new = (231*rhoSq_new - 91)*trueZ*trueZ*trueZ;
        double N0_new = 1 - N01_new + N02_new - N03_new;
	double N1_new = F1*F1 * F5 * hp_new*hp_new;
        double N2_new = F2*F2 * F5 * hm_new*hm_new;
	double N3_new = 4*F3 * cosV*cosV * h0_new*h0_new;
	double N4_new = 2*F3 * F5 * F7 * hp_new*hm_new;
        double N5_new = 4*F4 * F1 * F6 * cosV * cosChi * hp_new*h0_new;
        double N6_new = 4*F4 * F2 * F6 * cosV * cosChi * hm_new*h0_new;
        double N7_new 
	  = N1_new + N2_new + N3_new - N4_new - N5_new + N6_new;
        double FF_new = N0_new * N0_new * N7_new;
	double W_D(1);
        if (newDstarFF && FF_old>0) W_D = (FF_new / FF_old);

	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 * RN * num_peak[i][j][k];
	double side3D = conv * RN * num_side[i][j][k];
	double sub3D = peak3D - side3D;
	double var3D = num_peak[i][j][k] + num_side[i][j][k];
	double err3D = conv * RN * 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 << "Normalization info" << endl;
  outFile << "  N_old = " << N_old << endl;
  outFile << "  N_new = " << N_new << endl;
  outFile << "  RN = " << RN << endl;
  outFile << "  RN_err = " << RN_err << 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 parameters" << endl;
  outFile << "SP5/SP6" << endl;
  outFile << "  R1_old = " << R1_old << endl;
  outFile << "  R2_old = " << R2_old << endl;
  outFile << "  rhoSq_old = " << rhoSq_old << endl;
  outFile << "New values" << endl;
  outFile << "  R1_new = " << R1_new << endl;
  outFile << "  R2_new = " << R2_new << endl;
  outFile << "  rhoSq_new = " << rhoSq_new << endl;
  outFile << "Masses" << endl;
  outFile << "  m_Dstarp = " << m_Dstarp << endl;
  outFile << "  m_Dstar0 = " << m_Dstar0 << endl;
  outFile << "  m_Bp = " << m_Bp << endl;
  outFile << "  m_B0 = " << m_B0 << endl;
  outFile.close();

}