#include <iostream>
#include <fstream>

gROOT->Reset();

void mkHist_D1primeFFCorr_BpBm_test() {
  // Root macro to produce FF histogram
  int test = 0;

  // Name of output files
  TString outHistoName("hist_D1primeFFCorr_BpBm_test.root");
  TString outInfoFileName("hist_D1primeFFCorr_BpBm_test.txt");

  // binning
  enum { nbinsx = 100 };
  Float_t xmin = 0.0, xmax = 2.0;

  // Define histos
  TH1D *Hw = new TH1D("Hw", "Hw", nbinsx, 1.0, 1.5);
  TH1D *HoldFF = new TH1D("HoldFF", "HoldFF", nbinsx, xmin, 0.15);
  TH1D *HnewFF = new TH1D("HnewFF", "HnewFF", nbinsx, xmin, 0.15);
  TH1D *HcorrFF = new TH1D("HcorrFF", "HcorrFF", nbinsx, xmin, 2.0);


  // Masses
  double mB = 5.2790;
  double mD = 2.4300;

  // Get normalization factor
  double RN(1.0);
  normFactor(RN, mB, mD);
  cout << "  mB = " << mB << endl;
  cout << "  mD = " << mD << endl;
  cout << "  RN = " << RN << endl;


  //###################################
  // Chain input data
  cout << "start chain" << endl;
  TChain * chain = new TChain("ntp2");
  int nStop = 1952;
  if (test == 1) nStop = 5;
  for (Int_t i=1; i<=nStop; i++) {
    TString inFileNameBase("/hep24/khamano/hep03/R16b/Run2/SP-1235-BToDlnu-Run2-R16b");
    chain->Add(inFileNameBase+"-"+i+".root");
  }
  TTree * tree = chain;


  cout << "read tree" << endl;
  // variables to read out
  Int_t nDl_3D(0);
  Int_t isABhabha_3D(0);
  Int_t DDecayMode_3D[200];
  Int_t DlLepType_3D[200];
  Int_t tightK_3D[200];
  Float_t DMass_3D[200];
  double FitDprob_3D[200], FitBprob_3D[200];
  Float_t cosDlnonDl_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 trkEffCorrErr_3D[200], PIDCorrErr_3D[200];
  double PIDKCorr_3D[200], PIDLeptCorr_3D[200];
  double PIDKCorrErr_3D[200], PIDLeptCorrErr_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->SetBranchAddress("nDl", &nDl_3D);
  tree->SetBranchAddress("isABhabha", &isABhabha_3D);
  tree->SetBranchAddress("DDecayMode", &DDecayMode_3D);
  tree->SetBranchAddress("DlLepType", &DlLepType_3D);
  tree->SetBranchAddress("tightK", &tightK_3D);
  tree->SetBranchAddress("DMass", &DMass_3D);
  tree->SetBranchAddress("FitDprob", &FitDprob_3D);
  tree->SetBranchAddress("FitBprob", &FitBprob_3D);
  tree->SetBranchAddress("cosDlnonDl", &cosDlnonDl_3D);
  tree->SetBranchAddress("noExtraKorPi", &noExtraKorPi_3D);
  tree->SetBranchAddress("mcKaonMissID", &mcKaonMissID_3D);
  tree->SetBranchAddress("mcPiMissID", &mcPiMissID_3D);
  tree->SetBranchAddress("mcKPifromSameD", &mcKPifromSameD_3D);
  tree->SetBranchAddress("mcLepMissID", &mcLepMissID_3D);
  tree->SetBranchAddress("mcDirectLepFromB", &mcDirectLepFromB_3D);
  tree->SetBranchAddress("trueBDecMod", &trueBDecMod_3D);
  tree->SetBranchAddress("trueLepFromSameB", &trueLepFromSameB_3D);
  tree->SetBranchAddress("trkEffCorr", &trkEffCorr_3D);
  tree->SetBranchAddress("trkEffCorrErr", &trkEffCorrErr_3D);
  tree->SetBranchAddress("PIDCorr", &PIDCorr_3D);
  tree->SetBranchAddress("PIDCorrErr", &PIDCorrErr_3D);
  tree->SetBranchAddress("PIDKCorr", &PIDKCorr_3D);
  tree->SetBranchAddress("PIDKCorrErr", &PIDKCorrErr_3D);
  tree->SetBranchAddress("PIDLeptCorr", &PIDLeptCorr_3D);
  tree->SetBranchAddress("PIDLeptCorrErr", &PIDLeptCorrErr_3D);
  tree->SetBranchAddress("mcW", &mcW_3D);
  tree->SetBranchAddress("MCcosL", &MCcosL_3D);
  tree->SetBranchAddress("MCcosV", &MCcosV_3D);
  tree->SetBranchAddress("MCcosChi", &MCcosChi_3D);
  tree->SetBranchAddress("DMass", &varx_3D);
  tree->SetBranchAddress("DMom", &vary_3D);
  tree->SetBranchAddress("DLcosBY", &varz_3D);

  // Read out root-tuple entries and fill arrays
  Int_t nevent_3D =(Int_t)tree->GetEntries(); // number of events
  cout << "number of events = " << nevent_3D << endl;
  int numCand(0); //processed candidates
  int numCand_DDecMod(0); //after D decay mode cut
  int numCand_EvtSel(0); //after event selection cuts
  int numCand_selLep(0); //after lepton selection
  int numCand_Index(0); //before index selection
  int numSelCand(0); //selcted candidates

  for (Int_t i=0;i<nevent_3D;i++) { // loop over events
    if (i==1) cout << "event 1" << endl;
    if (i==100) cout << "event 100" << endl;
    if (i==1000) cout << "event 1000" << endl;
    if (i==10000) cout << "event 10000" << endl;
    if (i==100000) cout << "event 100000" << endl;
    if (i==1000000) cout << "event 1000000" << endl;

    tree->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 leptType = DlLepType_3D[j];
      int isTightK = tightK_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];
      double PidCorr = PIDCorr_3D[j];
      double PidKCorr = PIDKCorr_3D[j];
      double PidLCorr = PIDLeptCorr_3D[j];
      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];
      double probD = FitDprob_3D[j];
      double probB = FitBprob_3D[j];
      float cosThrust = cosDlnonDl_3D[j];
      if (cosThrust < 0) cosThrust = -cosDlnonDl_3D[j];

      numCand++;
      //if (DDecMod != 1) continue;
      numCand_DDecMod++;
      if (isTightK == 0) continue;
      //if (probD < 0.001) continue;
      //if (probB < 0.01) continue;
      //if (cosThrust > 0.92) continue;
      numCand_EvtSel++;
      //if (leptType != 1) continue;
      numCand_selLep++;

      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 == 1640 )  goodBDecMod=1;
      if (mcBDecMod == 2630 )  goodBDecMod=1;
      if (mcBDecMod == -1 )  goodBDecMod=1;
      if (mcBDecMod == -1 )  goodBDecMod=1;
      if (goodBDecMod != 1) continue;

      numCand_Index++;

      numSelCand++;

      // Get trueW
      Hw->Fill(trueW);
      
      double wx = 1.0;

      // Get oldFF
      double oldFF = oldDlnuFF(trueW, mB, mD, cosL, cosV, cosChi);
      HoldFF->Fill(oldFF);
      //cout << "oldfF =" << oldFF << endl;

      // Get newFF
      double newFF = newDlnuFF(trueW, mB, mD, cosL, cosV, cosChi);
      HnewFF->Fill(newFF);
      
      // Get weight
      if (oldFF > 0) wx = RN * newFF / oldFF;
      HcorrFF->Fill(wx);


    } //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 after event selection = " << numCand_EvtSel << endl;
  cout << "number of candidates after lepton selection = " << numCand_selLep << endl;
  cout << "number of candidates before index selection = " << numCand_Index << endl;
  cout << "number of selected candidates = " << numSelCand << endl;



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

  Hw->Write();
  HoldFF->Write();
  HnewFF->Write();
  HcorrFF->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 = 2" << endl;
  outFile << "  base directory = /home1x/OtherMounts/hep03/khamano/R16b" << endl;
  outFile << "  File name base = SP-1235-BToDlnu-Run2-R16b" << endl;
  outFile << "  start number = 1" << endl;
  outFile << "  stop number = " << nStop << endl;
  outFile << "################" << endl;
  outFile << "Output files" << endl;
  outFile << "  histograms : " << outHistoName << endl;
  outFile << "  information : " << outInfoFileName << endl;
  outFile << "################" << endl;
  outFile << "Signal-background selection" << endl;
  outFile << "  reqRealD = 1" << endl;
  outFile << "  reqRealL = 1"  << endl;
  outFile << "  reqDirectL = 1" << endl;
  outFile << "  reqCorr = 1" << endl;
  outFile << "  reqBMod = 1" << endl;
  outFile << "B decay mode selection" << endl;
  outFile << "  BMod1 = 1140" << endl;
  outFile << "  BMod2 = 2130" << endl;
  outFile << "  BMod3 = -1" << endl;
  outFile << "  BMod4 = -1" << 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 after event selection = " 
	  << numCand_EvtSel << endl;
  outFile << "  number of candidates after lepton selection = " 
	  << numCand_selLep << endl;
  outFile << "  number of candidates before index selection = " 
	  << numCand_Index << endl;
  outFile << "  number of candidates selected = " << numSelCand << endl;
  outFile << "################" << endl;
  outFile << "Normalization factor" << endl;
  outFile << "  RN = " << RN << endl;
  outFile << "################" << endl;
  outFile << "Binning info" << endl;
  outFile << "  nbinsx = 70" << endl;
  outFile << "  xmin = " << xmin << endl; 
  outFile << "  xmax = " << xmax  << endl;
  outFile.close();


  //#############################################################
  //Plot histos
  
  if (test == 1) {
    cout << "plotting\n";
    TCanvas * cx = new TCanvas("cx","cx",0,0,600,600);
    cx->Divide(2,2);
    gStyle->SetOptStat(000000);
    
    cx->cd(1);
    HoldFF->Draw();
    
    cx->cd(2);
    HnewFF->Draw();
    
    cx->cd(3);
    HcorrFF->Draw(); 

    cx->cd(4);
    Hw->Draw();
    
  }
  
}




void normFactor(double &RNM, double mB, double mD) {

  enum { n = 7 };
  double xi[n];
  double wi[n];
  xi[0] = 0.0;
  wi[0] = 0.417959183673496;
  xi[1] = 0.405845151377397;
  wi[1] = 0.381830050505119;
  xi[2] = -xi[1];
  wi[2] = wi[1];
  xi[3] = 0.741531185599394;
  wi[3] = 0.279705391489277;
  xi[4] = -xi[3];
  wi[4] = wi[3];
  xi[5] = 0.949107912342759;
  wi[5] = 0.129484966168870;
  xi[6] = -xi[5];
  wi[6] = wi[5];


  double wMax = (mB*mB + mD*mD)/(2.0*mB*mD);

  double oldIntBF(0.0);
  double newIntBF(0.0);
  for (Int_t i = 0; i<n ; i++) {
    double yi = (wMax-1.0)*xi[i]/2.0 + (wMax+1.0)/2.0;

    double oldDW = oldDlnuDW(yi, mB, mD);
    oldIntBF += wi[i]*oldDW;

    double newDW = newDlnuDW(yi, mB, mD);
    newIntBF += wi[i]*newDW;
  }
  oldIntBF = (wMax-1.0)*oldIntBF/2.0;
  newIntBF = (wMax-1.0)*newIntBF/2.0;
  cout << "oldIntBF = " << oldIntBF << endl;
  cout << "newIntBF = " << newIntBF << endl;

  if (newIntBF>0) RNM = oldIntBF / newIntBF;
  cout << "RNM = " << RNM << endl;

}




double oldDlnuFF(double trueW, double mB, double mD, double cosL, double cosV, double cosChi) {
  double gA(0);
  double gV1(0);
  double gV3rgV2(0);
  ISGW2D1prime(trueW, mB, mD, gA, gV1, gV3rgV2);

  double oldFF = D1primeFF(trueW, mB, mD, cosL, gA, gV1, gV3rgV2);
  return oldFF;
}


double newDlnuFF(double trueW, double mB, double mD, double cosL, double cosV, double cosChi) {
  double gA(0);
  double gV1(0);
  double gV3rgV2(0);
  LLSWD1primeB1(trueW, mB, mD, gA, gV1, gV3rgV2);

  double newFF = D1primeFF(trueW, mB, mD, cosL, gA, gV1, gV3rgV2);
  return newFF;
}


double oldDlnuDW(double trueW, double mB, double mD) {
  double gA(0);
  double gV1(0);
  double gV3rgV2(0);
  ISGW2D1prime(trueW, mB, mD, gA, gV1, gV3rgV2);

  double oldFF = D1primeDW(trueW, mB, mD, gA, gV1, gV3rgV2);
  double oldDW = sqrt(trueW*trueW - 1.0)*oldFF;
  return oldDW;
}


double newDlnuDW(double trueW, double mB, double mD) {
  double gA(0);
  double gV1(0);
  double gV3rgV2(0);
  LLSWD1primeB1(trueW, mB, mD, gA, gV1, gV3rgV2);

  double newFF = D1primeDW(trueW, mB, mD, gA, gV1, gV3rgV2);
  double newDW = sqrt(trueW*trueW - 1.0)*newFF;
  return newDW;
}



double D1primeFF(double trueW, double mB, double mD, double cosL, double gA, double gV1, double gV3rgV2) {
  double cosT = -cosL;
  double r = mD/mB;
  double r_tm1 = (trueW - r)*gV1 + (trueW*trueW - 1.0)*gV3rgV2;
  double tm1 = (1.0 - cosT*cosT) * r_tm1 * r_tm1;
  double tm2a = (1.0 + cosT*cosT)* (gV1*gV1 + (trueW*trueW-1.0)*gA*gA);
  double tm2b = -4*cosT*sqrt(trueW*trueW-1.0)*gV1*gA;
  double tm2 = (1.0 - 2.0*r*trueW + r*r) * (tm2a + tm2b);
  double DW = tm1 + tm2;
  return DW;
}

double D1primeDW(double trueW, double mB, double mD, double gA, double gV1, double gV3rgV2) {
  double r = mD/mB;
  double r_tm1 = (trueW - r)*gV1 + (trueW*trueW - 1.0)*gV3rgV2;
  double tm1 = r_tm1 * r_tm1;
  double tm2 
    = 2.0 * (1.0 - 2.0*r*trueW + r*r) * (gV1*gV1 + (trueW*trueW-1.0)*gA*gA);
  double DW = tm1 + tm2;
  return DW;
}



void ISGW2D1prime(double trueW, double mB, double mD, double &gA, double &gV1, double &gV3rgV2) {
  
  // Get q-square
  double q2 = mB*mB + mD*mD - 2*mB*mD*trueW;
  // Get R
  double mcR = 2.0*sqrt(mB*mD) / (mB+mD);
  // Get r
  double r = mD/mB;

  // Get FF
  double t=q2;
  double mass=mD;
  double lf, qf, cpf, cmf;

  double mtb, mbb;
  double msd, mx,mb,nf,nfp;
  double msq,bx2,mbx,mtx;
  double f5cppcm,f5cpmcm,f5,ql,ll,cppcm,cpmcm,f5q,f5l;
  double mqm,msb,bb2,mup,mum,bbx2,tm,wt,r2;

  msb=5.2;
  msd=0.33;
  bb2=0.431*0.431;
  mbb=5.31;
  nf = 4.0;

  msq=1.82;
  bx2=0.33*0.33;
  mbx=(3.0*2.49+2.40)/4.0;
  nfp = 3.0;

  mtb = msb + msd;
  mtx = msq + msd;
  mb=mB;
  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;
  r2=3.0/(4.0*msb*msq)+3*msd*msd/(2*mbb*mbx*bbx2) + 
    (16.0/(mbb*mbx*(33.0-2.0*nfp)))*
    log(EvtGetas(mqm)/EvtGetas(msq));

  f5 = sqrt(mtx/mtb)*pow(sqrt(bx2*bb2)/bbx2,5.0/2.0) /
    (pow((1.0+r2*(tm-t)/18.0),3.0));

  f5q = f5*pow(( mbb / mtb ),-0.5)*pow((mbx/mtx),-0.5);
  f5l = f5*pow(( mbb / mtb ),0.5)*pow((mbx/mtx),0.5);
  f5cppcm = f5*pow(( mbb / mtb ),-1.5)*pow((mbx/mtx),0.5);
  f5cpmcm = f5*pow(( mbb / mtb ),-0.5)*pow((mbx/mtx),-0.5);
  
  if (msq == msd) { 
    
    ql = -1.0*(msd*(5.0+wt)*f5q/(2.0*mtx*sqrt(bb2)*6.0));
    
    ll = -1.0*mtb*sqrt(bb2)*f5l*(1/mum+ ( (msd*mtx*(wt-1)/bb2)*
					  ( (5.0+wt)/(6.0*msq)-(msd*bb2)/(2.0*mum*mtx*bbx2))));
    
    cppcm = (-1.0*(msd*mtx*f5cppcm/(2.0*msq*mtb*sqrt(bb2)))*
	     (1-(msd*msq*bb2)/(2.0*mtx*mum*bbx2)));
    
    cpmcm = 1.0*(msd*mtx*f5cpmcm/(2.0*msq*mtb*sqrt(bb2)))*
      (((wt+2.0)/3.0)-(msd*msq*bb2)/(2.0*mtx*mum*bbx2))
      *(mtb/mtx);
  } else {

    ql = f5q*sqrt(1.0/6.0)*msd/(sqrt(bb2)*mtx)*
      (1.0-bb2*mtb/(4.0*msd*msq*msb));
    ll = f5l*sqrt(2.0/3.0)*mtb*sqrt(bb2)*(1.0/(2.0*msq) - 3.0/(2.0*msb) +
					  msd*mtx*(wt-1)/bb2*(1.0/msq-msd*bb2/(2.0*mum*mtx*bbx2)));  
    cppcm = msd*msd*bx2*f5cppcm/(sqrt(6.0)*mtb*msq*sqrt(bb2)*bbx2);
    cpmcm = -sqrt(2.0/3.0)*msd*f5cpmcm/(sqrt(bb2)*mtx)*
      (1+msd*bx2/(2.0*msq*bbx2));
  }

  /*
  //smooth out the mass(meson) dependence a little
  double parMass=mB;
  double q2max = parMass*parMass + mass*mass - 2.0*parMass*mass;
  double massNom= mD;
  double q2maxNom = parMass*parMass + massNom*massNom - 2.0*parMass*massNom;
  double q2maxin=sqrt(q2maxNom/q2max);
  if ( q2maxin > 1000. ) q2maxin=1000.;
  ql*=q2maxin;
  ll*=q2maxin;
  cppcm*=q2maxin;
  cpmcm*=q2maxin;
  */

  qf = ql;
  lf = ll;
  cpf = (cppcm + cpmcm)/2.0;
  cmf = (cppcm - cpmcm)/2.0;

  double conv = mcR * mB * (1.0 + r);
  //gA = -conv * qf;
  gA = conv * qf;
  gV1 = 2.0 * lf / conv;
  gV3rgV2 = conv * cpf;

}

double EvtGetas ( double mass ) {

  double pi = 3.141592654;

  double lqcd2 = 0.04;
  double nflav = 4;
  double temp = 0.6;
  
  if ( mass > 0.6 ) {
    if ( mass < 1.85 ) {
      nflav = 3.0;}
    
    temp = 12.0*pi / ( 33.0 - 2.0*nflav) /
      log( mass*mass/lqcd2);
  }
  return temp;
  
} //EvtGetas



double LLSWD1primeA(double trueW, double mB, double mD) {

  double tau_1_org = 0.60;

  double r = mD/mB;

  double m_b = 4.8;
  double m_c = 1.4;
  double e_b = 1.0/(2*m_b);
  double e_c = 1.0/(2*m_c);

  double Lambda = 0.4;
  double Lambda_p = 0.39 + Lambda;
  double Lambda_s = 0.35 + Lambda;
  
  //double Br_BD1lnu = 0.00406;
  double Br_BD1lnu = 0.0060;
  double Br_BD1lnu_org = 0.0060;
  double tau_1 = tau_1_org*sqrt(Br_BD1lnu/Br_BD1lnu_org);
  double tau_p = -1.5;
  double zeta_1 = tau_1 * 2.0 / sqrt(3.0);
  //double zeta_1 = 1.0;
  double zeta_p = tau_p + 0.5;

  double s10 = (1-r)*(1-r)*(e_c-3.0*e_b)*(e_c-3.0*e_b)
    *(Lambda_s-Lambda)*(Lambda_s-Lambda);
  double s11 = -2.0*(1-r*r)*(e_c-3.0*e_b)*(Lambda_s-Lambda);
  double s12 = (1.0+r)*(1.0+r);
  double t10 = (e_c - 3.0*e_b)*(e_c - 3.0*e_b)*(Lambda_s-Lambda)
    *(Lambda_s-Lambda);
  double t11 = 2.0 + 4.0*(e_c-3.0*e_b)*(Lambda_s-Lambda);
  double t12 = 2.0*(1 + 2.0*zeta_p);

  double newFF0 = s10 + 2.0*(1.0 - 2.0*r*trueW + r*r)*t10;
  double newFF1 = (trueW - 1.0) * (s11 + 2.0*(1.0 - 2.0*r*trueW + r*r)*t11);
  double newFF2 = (trueW - 1.0)*(trueW - 1.0) 
    * (s12 + 2.0*(1.0 - 2.0*r*trueW + r*r)*t12);
  double newFF = zeta_1*zeta_1 * (newFF0 + newFF1 + newFF2);

  return newFF;

}
 

 
void LLSWD1primeB1(double trueW, double mB, double mD, double &gA, double &gV1, double &gV3rgV2) {
  
  int opt = 1;
  double tau_1_org = 0.71;

  LLSWD1primeB(trueW, mB, mD, opt, tau_1_org, gA, gV1, gV3rgV2);

}


void LLSWD1primeB2(double trueW, double mB, double mD, double &gA, double &gV1, double &gV3rgV2) {

  int opt = 2;
  double tau_1_org = 0.75;

  LLSWD1primeB(trueW, mB, mD, opt, tau_1_org, gA, gV1, gV3rgV2);

}


void LLSWD1primeB(double trueW, double mB, double mD, int opt, double tau_1_org, double &gA, double &gV1, double &gV3rgV2) {

  // setup parameters
  double r = mD/mB;

  double m_b = 4.8;
  double m_c = 1.4;
  double e_b = 1.0/(2*m_b);
  double e_c = 1.0/(2*m_c);

  double Lambda = 0.4;
  double Lambda_p = 0.39 + Lambda;
  double Lambda_s = 0.35 + Lambda;
  
  //double Br_BD1lnu = 0.00406;
  double Br_BD1lnu = 0.0060;
  double Br_BD1lnu_org = 0.0060;
  double tau_1_org = 0.71;
  double tau_1 = tau_1_org*sqrt(Br_BD1lnu/Br_BD1lnu_org);
  double tau_p = -1.5;

  // functions tau, zeta
  double tau = tau_1 * (1.0 + tau_p * (trueW - 1.0));
  //double zeta = (trueW + 1.0)*tau / sqrt(3.0);
  double zeta_1 = tau_1 * 2.0 /  sqrt(3.0);
  //double zeta_1 = 1.0;
  double zeta_p = tau_p + 0.5;
  double zeta = zeta_1 * (1.0 + zeta_p * (trueW - 1.0));

  double zeta1 = 0.0;
  double zeta2 = 0.0;
    
  if (opt == 2) {
    zeta1 = Lambda * zeta;
    zeta2 = -Lambda_s * zeta;
  }
  
  // form factors
  double gA_t1 = zeta;
  double gA_t2 = e_c*(trueW*Lambda_s - Lambda)*zeta / (trueW+1.0);
  double gA_t3 
    = -e_b*(((2.0*trueW+1.0)*Lambda_s-(trueW+2.0)*Lambda)*zeta/(trueW+1.0)
	    - 2.0*(trueW-1.0)*zeta1);
  gA = gA_t1 + gA_t2 + gA_t3;
  
  double gV1_t1 = (trueW-1.0)*zeta;
  double gV1_t2 = e_c*(trueW*Lambda_s - Lambda)*zeta; 
  double gV1_t3 
    = -e_b*(((2.0*trueW+1.0)*Lambda_s-(trueW+2.0)*Lambda)*zeta 
	    - 2.0*(trueW*trueW-1.0)*zeta1);
  gV1 = gV1_t1 + gV1_t2 + gV1_t3;
  
  double gV2 = 2.0*e_c*zeta1;
  
  double gV3_t1 = -zeta;
  double gV3_t2  
    = -e_c*((trueW*Lambda_s - Lambda)*zeta / (trueW+1.0) + 2.0*zeta1);
  double gV3_t3 
    = e_b*(((2.0*trueW+1.0)*Lambda_s-(trueW+2.0)*Lambda)*zeta/(trueW+1.0) 
	   - 2.0*(trueW-1.0)*zeta1);
  double gV3 = gV3_t1 + gV3_t2 + gV3_t3;
    
  // calculate FF
  gV3rgV2 = gV3 + r*gV2;
  
}