{
//int nrun = 1467;
//cout<<"run number ?"<<endl;
//cin>>nrun;
TFile *f=new TFile("/data/doug/e94107/e94107_det_v2_4288.root","UPDATE");
//if(!f){
//  cout<<"cannot locate root file.."<<endl;
//  return();
}
Double_t emom = 2039.67;
Double_t momset = 2030.; // Delta scan momentum setting
//cout<<"momtum setting ?"<<endl;
//cin>>momset:
switch (nrun){
 case 4288:
   momset = 2000.;
   emom=2000.;
   break;
 case 1458:
   momset = 2105.;
   break;
 case 1460:
   momset = 2080.;
   break;
 case 1464:
   momset = 2060.;
   break;
 case 1467:
   momset = 2030.;
   break;
 case 1471:
   momset = 2005.;
   break;
 case 1472:
   momset = 1985.;
   break;
 case 1476:
   momset = 1945.;
   break;
 case 1457:
   momset = 2105.;
   break;
 case 1461:
   momset = 2080.;
   break;
 case 1465:
   momset = 2060.;
   break;
 case 1466:
   momset = 2030.;
   break;
 case 1470:
   momset = 2005.;
   break;
 case 1473:
   momset = 1985.;
   break;
 case 1477:
   momset = 1945.;
   break;
 case 1099:
   momset = 2983.20;
   emom = 3028.43;
   break;
 case 1071:
   momset = 3079.80;
   emom = 3028.43;
   break;
 default:
   return;
}

Double_t L_th_0=0.104720; // central angel 6 degree
Double_t R_th_0=-0.104720;

//Double_t el1=0.069; //carbon
//Double_t el2=0.070;

Double_t el1=0.498+0.029; //2.5mm waterfall+ 1mil steel //before interaction
Double_t el2=0.498+0.029; //after

Double_t L_dpkin;
Double_t L_dp;
Double_t L_th_tg;
Double_t L_ph_tg;
//Double_t L_d_tg;

Double_t R_dpkin;
Double_t R_dp;
Double_t R_th_tg;
Double_t R_ph_tg;
//Double_t R_d_tg;


Double_t mom;
Double_t momc;
Double_t m_A;
//Double_t tmp_E;
Double_t p_corr;
Double_t theta;

const Double_t m_u =931.494;
const Double_t m_O = 14903.345; //15.9994u
const Double_t m_C = 11177.928; //12u
const Double_t m_H = 938.272;     
const Double_t m_Be = 8394.793; //9.012182u
const Double_t m_Al = 25133.141; //26.981538u


if(!(T->GetBranch("L.gold.dp"))) return;
T->SetBranchAddress("L.gold.dp",&L_dp);
T->SetBranchAddress("L.gold.th",&L_th_tg);
T->SetBranchAddress("L.gold.ph",&L_ph_tg);
T->SetBranchAddress("R.gold.dp",&R_dp);
T->SetBranchAddress("R.gold.th",&R_th_tg);
T->SetBranchAddress("R.gold.ph",&R_ph_tg);

if(T->GetBranch("L.dp.kin")) {
  cout<<"Tree already modified, end."<<endl;
  return();
}

else{
  // define new branch for dp.kin
  TBranch* b_ldpkin=T->Branch("L.dp.kin",&L_dpkin,"L.dp.kin/D");
  TBranch* b_rdpkin=T->Branch("R.dp.kin",&R_dpkin,"R.dp.kin/D");
}

// generate
Int_t nentries = (Int_t) T->GetEntries();
T->SetBranchStatus("*",0);
T->SetBranchStatus("L.gold.dp",1);
T->SetBranchStatus("L.gold.th",1);
T->SetBranchStatus("L.gold.ph",1);
T->SetBranchStatus("L.dp.kin",1);
T->SetBranchStatus("R.gold.dp",1);
T->SetBranchStatus("R.gold.th",1);
T->SetBranchStatus("R.gold.ph",1);
T->SetBranchStatus("R.dp.kin",1);
for (Int_t i=0; i<nentries;i++){  
  //left 
  T->GetEntry(i);
  mom = (1+L_dp)*momset;
  Double_t ediff=abs(emom-el1-mom);

  // mass defination
  if(ediff==0) m_A=m_O;
  else{
    m_A = pow(emom-el1,2)*(1-cos(L_th_0))/ediff;   
    if ((m_A-m_H)<800) m_A = m_H;
    else m_A =m_O;
  }
  // accurate dp.kin
  theta = atan(sqrt(pow(tan(L_th_0+atan(L_ph_tg)),2)+pow(L_th_tg,2)));
  p_corr = m_A*(emom-el1)/(m_A+(emom-el1)*(1-cos(L_th_0)))-m_A*(emom-el1)/(m_A+(emom-el1)*(1-cos(theta)));
  momc = mom+p_corr;

  L_dpkin = (momc - momset)/momset;
  // fill in tree
  b_ldpkin->Fill();

  //right
  mom = (1+R_dp)*momset;
  ediff=abs(emom-el1-mom);
  if(ediff==0) m_A=m_O;
  else{
    m_A = pow(emom-el1,2)*(1-cos(R_th_0))/ediff;
    if ((m_A-m_H)<800) m_A = m_H;
    else m_A =m_O;
  }
  theta = atan(sqrt(pow(tan(R_th_0+atan(R_ph_tg)),2)+pow(R_th_tg,2)));
  p_corr = m_A*(emom-el1)/(m_A+(emom-el1)*(1-cos(R_th_0)))-m_A*(emom-el1)/(m_A+(emom-el1)*(1-cos(theta)));
  momc = mom+p_corr;
  R_dpkin = (momc - momset)/momset;
  b_rdpkin->Fill();
}
T->Write("",TObject::kOverwrite);
f->Close();
}

 // nilanga's dp.kin
  //  momc = mom*(1+(mom*sin(th_0)*ph_tg)/m_A);
  
  // Yi's dp.kin
  //  d_tg = ph_tg + th_tg*th_tg/(2*th_0);
  //  momc = mom + m_A*sin(th_0)*d_tg/pow(1-cos(th_0)+m_A/(emom-el1),2);

  // judge target mass
  //  if(abs(m_A-m_O)<1500) m_A = m_O;
  //  else if(abs(m_A-m_C)<1100) m_A = m_C;
  //  else if(abs(m_A-m_Be)<800) m_A = m_Be;
  //  else if(abs(m_A-m_H)<100) m_A = m_H;