void analysis(Int_t run_i,Int_t run_f,Int_t df_run) {
  gROOT->Reset();
  gStyle->SetOptStat(0);

  TH1F *hsthe = new TH1F("hstheta","HMS theta in degree",100,0,50);
  TH1F *hsdel = new TH1F("hsdelta","HMS \delta E",100,-12,12);
  TH1F *hsen = new TH1F("hse","HMS energy in GeV",100,1.,5.);
  
  TH1F *Q2 = new TH1F("Q2","Q^{2} in GeV^{2}",100,0.,3.);
  TH1F *W = new TH1F("W","W in GeV",100,0.8,3.);
  Float_t Wmin = 1.05; //To apply dilution factor in [0.0,3.0] with 200 bins, Wmin should be divided by 0.015 GeV.
  Int_t Wbin = 42; // For 15MeV resoultion deltaW=1.05 for Wbin=70
  Float_t Wbinsize = 0.015; // 15MeV
  Float_t Wmax = Wmin + Wbin*Wbinsize;
  TH2F *Q2vsW = new TH2F("Q2vsW","Q^{2} vs W",Wbin,Wmin,Wmax,100,0.,3.);
  //TH1F *W_p = new TH1F("W_p","W of positive helicity",Wbin,Wmin,Wmax);
  //TH1F *W_n = new TH1F("W_n","W of negative helicity",Wbin,Wmin,Wmax);
  //TH1D *W_C_p = new TH1D("W_C_p","W of positive helicity, charge normalized, LT correction muliplied",Wbin,Wmin,Wmax);
  //TH1D *W_C_n = new TH1D("W_C_n","W of negative helicity, charge normalized, LT correction muliplied",Wbin,Wmin,Wmax);
  //TH1D *W_CPbPt_p = new TH1D("W_C_p","W of positive helicity, charge normalized, divided by Pb*Pt",Wbin,Wmin,Wmax);
  //TH1D *W_CPbPt_n = new TH1D("W_C_n","W of negative helicity, charge normalized, divided by Pb*Pt",Wbin,Wmin,Wmax);

  TH1D *W_asym_raw_all = new TH1D("W_asym_raw_all","(N_{+} - N_{-})/(N_{+} + N_{-}) along W for all runs",Wbin,Wmin,Wmax);
  TH1D *Error_W_asym_raw_all = new TH1D("Error_W_asym_all","Accumulated inverse of squared error of W_asym",Wbin,Wmin,Wmax);
  TH1D *W_asym_CPbPt_all = new TH1D("W_asym_CPbPt_all","#frac{1}{PbPt} (N_{+}/C_{+} - N_{-}/C_{-})/(N_{+}/C_{+} + N_{-}/C_{-}) along W for all runs",Wbin,Wmin,Wmax);
  TH1D *Error_W_asym_CPbPt_all = new TH1D("Error_W_asym_CPbPt_all","Accumulated inverse of squared error of W_asym_CPbPt",Wbin,Wmin,Wmax);
  TH1D *W_asym_CPbPtdf_all = new TH1D("W_asym_CPbPtdf_all","#frac{1}{PbPtdf} (N_{+}/C_{+} - N_{-}/C_{-})/(N_{+}/C_{+} + N_{-}/C_{-}) along W for all runs",Wbin,Wmin,Wmax);

  TChain chain("h9010");
  Long_t *id,*size,*flags,*mt;

  for (Int_t nrun=run_i; nrun<=run_f; nrun++) {
    Int_t nfile = 0;
    chain.Reset("");
    for (Int_t index=1; index<10; index++) {
      if (gSystem->GetPathInfo(filename(nrun,index),id,size,flags,mt)==0) {
	chain.Add(filename(nrun,index));
	cout << filename(nrun,index) << " loaded : " << endl;
	nfile++;
      }
    }

    if (nfile!=0) {
//chain.Add(filename(72957,2));
//chain.Add(filename(72957,3));

      
      TH1F *W_p = new TH1F("W_p","W of positive helicity",Wbin,Wmin,Wmax);
      TH1F *W_n = new TH1F("W_n","W of negative helicity",Wbin,Wmin,Wmax);

      //Kinematic variables
      Run_info *r = new Run_info(nrun);
      r->summary();
      Int_t hwplate = r->hwplate();
      Float_t Pb = r->Pb()/100.; // Pb in [0,1]
      Float_t Pt = r->Pt()/100.; // Pt in [0,1]

      Int_t hel_flip = 0;
      if(hwplate==1) hel_flip = -1;
      if(hwplate==0) hel_flip = 1;
      //      cout << "hel_flip = " << hel_flip << " hwplate = " << hwplate << endl;

      Float_t C_p_over_n = 1.;
      if(hel_flip > 0) C_p_over_n = r->Charge_help()/r->Charge_helm();
      if(hel_flip < 0) C_p_over_n = r->Charge_helm()/r->Charge_help();
      //      cout << "Charge_help = " << r->Charge_help() << " Charge_helm = " << r->Charge_helm() << " Cp/Cn = " << C_p_over_n << endl ;

      Float_t LT_cor = r->LT_cor();
      cout << "LT_correction = " << LT_cor << endl;

      Kin_var *k = new Kin_var(r->BeamE(),r->P_HMS(),r->Theta_HMS());
      //k->set(4700,2200,16);
      cout << k->Q2()/1e6 << " " << sqrt(k->hsca2pe()) << endl;
      k->summary();
      //delete r;
      //delete k;

      TObjArray Hlist(0);
      
      Int_t nentries = chain.GetEntries();
      Float_t hsdelta, hstheta, hse, hsshtrk;
      //Float_t polariza;
      Int_t helicite;
      const Int_t r2d = 180./TMath::Pi();
      
      chain.SetBranchAddress("hsdelta",&hsdelta);
      chain.SetBranchAddress("hstheta",&hstheta);
      chain.SetBranchAddress("hse",&hse);
      chain.SetBranchAddress("hsshtrk",&hsshtrk);
      chain.SetBranchAddress("helicite",&helicite);
      //chain.SetBranchAddress("polariza",&polariza);
      
      cout << "Number of entries :" << nentries << endl;
      
      for (Int_t i = 0; i<nentries; i++) {
	chain.GetEntry(i);
      
	if(abs(hsdelta)<8. && hse>0. && (hsshtrk/hse)>0.4) {
	  k->set(r->BeamE(),hse*1e3,hstheta*r2d);
	  //k->set(r->BeamE(),3200.*(1.+hsdelta/100.),hstheta*r2d);
	  Float_t sW = sqrt(k->W2())/1e3; // W in GeV
	  //Float_t Pt = polariza/100.; // Pt in [0,1]
	  hsthe->Fill(k->theta());
	  hsdel->Fill(hsdelta);
	  hsen->Fill(hse);
	  Q2->Fill(k->Q2()/1e6);
	  W->Fill(sW);
	  Q2vsW->Fill(sW,k->Q2()/1e6);
	  
	  if(TMath::Abs(Pt)>0.2) {
	    if(helicite*hel_flip>0) {
	      W_p->Fill(sW);
	      //W_C_p->Fill(sW,LT_cor); // charge, LT correction
	      //W_CPbPt_p->Fill(sW,1./Pb/Pt); // charge, Pb, Pt
	    }
	    else if(helicite*hel_flip<0) {
	      W_n->Fill(sW);
	      //W_C_n->Fill(sW,LT_cor*C_p_over_n); // charge, LT correction
	      //W_CPbPt_n->Fill(sW,C_p_over_n/Pb/Pt); // charge, Pb, Pt
	    }
	  }
	}
      }
      
      //Getting asymmetry
      TH1D *W_p_p_n = new TH1D("W_p_p_n","N_{+} + N_{-} along W",Wbin,Wmin,Wmax);
      TH1D *W_p_m_n = new TH1D("W_p_m_n","N_{+} - N_{-} along W",Wbin,Wmin,Wmax);
      TH1D *W_CPbPt_p_p_n = new TH1D("W_CPbPt_p_p_n","#frac{1}{PbPt} (N_{+}/C_{+} + N_{-}/C_{-}) along W",Wbin,Wmin,Wmax);
      TH1D *W_CPbPt_p_m_n = new TH1D("W_CPbPt_p_m_n","#frac{1}{PbPt} (N_{+}/C_{+} - N_{-}/C_{-}) along W",Wbin,Wmin,Wmax);

      TH1D *W_asym_raw = new TH1D("W_asym_raw","(N_{+} - N_{-})/(N_{+} + N_{-}) along W",Wbin,Wmin,Wmax);
      TH1D *W_asym_CPbPt = new TH1D("W_asym_CPbPt","#frac{1}{PbPt} (N_{+}/C_{+} - N_{-}/C_{-})/(N_{+}/C_{+} + N_{-}/C_{-}) along W",Wbin,Wmin,Wmax);
      
      W_p_p_n->Add(W_p,W_n,1,1);
      W_p_m_n->Add(W_p,W_n,1,-1);
      W_asym_raw->Divide(W_p_m_n,W_p_p_n,1,1,"");
      
      W_CPbPt_p_p_n->Add(W_p,W_n,1.,1.*C_p_over_n);
      W_CPbPt_p_m_n->Add(W_p,W_n,1.,-1.*C_p_over_n);
      W_CPbPt_p_m_n->Scale(LT_cor/Pb/Pt,"");
      W_asym_CPbPt->Divide(W_CPbPt_p_m_n,W_CPbPt_p_p_n,1,1,"");
      //W_CPbPt_p_m_n->Add(W_CPbPt_p,W_CPbPt_n,1,-1);
      //W_asym_CPbPt->Divide(W_CPbPt_p_m_n,W_p_p_n,1,1,"");
      //W_p->Add(W_n,-1);
      //W_asym_raw->Divide(W_p,W_p_p_n,1,1,"");

      //Error calculation
      TH1D *E2_W_asym_raw = new TH1D("E2_W_asym_raw","Inverse of squared errors of W_asym_raw",Wbin,Wmin,Wmax);
      TH1D *W_asym_raw_weighted = new TH1D("W_asym_raw_weighted","W_asym_raw divided by error^2",Wbin,Wmin,Wmax);
      TH1D *E2_W_asym_CPbPt_CPbPt = new TH1D("E2_W_asym_CPbPt","Inverse of squared errors of W_asym_CPbPt",Wbin,Wmin,Wmax);
      TH1D *W_asym_CPbPt_weighted = new TH1D("W_asym_CPbPt_weighted","W_asym_CPbPt divided by error^2",Wbin,Wmin,Wmax);

      for(Int_t ibin = 1; ibin<=Wbin; ibin++) {
	cout << ibin << " " << W_p->GetBinCenter(ibin) << " " << W_p->GetBinContent(ibin) << endl;
	Float_t N_p = W_p->GetBinContent(ibin);
	Float_t N_n = W_n->GetBinContent(ibin);
	
	Float_t sigma2_W_asym_raw = 4*N_p*N_n/pow(N_p+N_n,3);
	Float_t sigma2_W_asym_CPbPt = 4*pow(LT_cor/Pb/Pt,2)*pow(C_p_over_n,2)*N_p*N_n*(N_p+N_n)/pow(N_p+N_n*C_p_over_n,4);

	//cout << ibin << " " << W_asym_raw->GetBinContent(ibin) << " " << W_asym_raw->GetBinError(ibin) << endl;
	W_asym_raw->SetBinError(ibin,sqrt(sigma2_W_asym_raw));
	W_asym_CPbPt->SetBinError(ibin,sqrt(sigma2_W_asym_CPbPt));
	//cout << ibin << " " << W_asym_raw->GetBinContent(ibin) << " " << W_asym_raw->GetBinError(ibin) << endl;

	E2_W_asym_raw->SetBinContent(ibin,1./sigma2_W_asym_raw);
	E2_W_asym_CPbPt->SetBinContent(ibin,1./sigma2_W_asym_CPbPt);
      }

      //Accumulation of asymmetry/sigma^2 and 1/sigma^2
      W_asym_raw_weighted->Multiply(W_asym_raw,E2_W_asym_raw,1,1,"");
      W_asym_raw_all->Add(W_asym_raw_weighted,1);
      Error_W_asym_raw_all->Add(E2_W_asym_raw,1);
      W_asym_CPbPt_weighted->Multiply(W_asym_CPbPt,E2_W_asym_CPbPt,1,1,"");
      W_asym_CPbPt_all->Add(W_asym_CPbPt_weighted,1);
      Error_W_asym_CPbPt_all->Add(E2_W_asym_CPbPt,1);
      
      delete k;
      delete r;
    }
  }

  //Calculation of final asymmetry of all runs
  for(Int_t ibin = 1; ibin<=Wbin; ibin++) {
    Float_t sum_asym_raw_of_sigma2 = W_asym_raw_all->GetBinContent(ibin);
    Float_t sum_raw_of_sigma2 = Error_W_asym_raw_all->GetBinContent(ibin);
    Float_t sum_asym_CPbPt_of_sigma2 = W_asym_CPbPt_all->GetBinContent(ibin);
    Float_t sum_CPbPt_of_sigma2 = Error_W_asym_CPbPt_all->GetBinContent(ibin);

    W_asym_raw_all->SetBinContent(ibin,sum_asym_raw_of_sigma2/sum_raw_of_sigma2);
    W_asym_raw_all->SetBinError(ibin,sqrt(1./sum_raw_of_sigma2));
    W_asym_CPbPt_all->SetBinContent(ibin,sum_asym_CPbPt_of_sigma2/sum_CPbPt_of_sigma2);
    W_asym_CPbPt_all->SetBinError(ibin,sqrt(1./sum_CPbPt_of_sigma2));
  }

  //Loading the dilution factor (If df_run = 1, df = 1.0)
  Long_t *id,*size,*flags,*mt;
  
  if (gSystem->GetPathInfo(df_filename(df_run),id,size,flags,mt)==0) {
    //TFile *df_file = new TFile(df_filename(df_run));
    //TH1D *df = (TH1D*)df_file->Get("h300");
    TH1D *df = new TH1D("df","Dilution factor without error",Wbin,Wmin,Wmax);
    
    ifstream df_text;
    df_text.open(df_filename(df_run));
    Int_t iWmin = ceil(Wmin/Wbinsize);
    cout << "iWmin = " << iWmin << endl;
    cout << "index \t W \t df \t error" << endl;

    Int_t line_number = 1;
    Int_t ibin = 1;
    TString line;
    while(!df_text.eof()) {
      
      line.ReadLine(df_text);

      if(line_number>iWmin && ibin<=Wbin) {
	TObjArray *sublist = line.Tokenize("\t ");
	TObjString *sub = sublist->At(3);
	Float_t wprnh = atof(sub->GetString().Data());
	
	df->SetBinContent(ibin,wprnh);
	df->SetBinError(ibin,0);
	
	cout << line.Data() << endl;
	cout << ibin << " " << df->GetBinCenter(ibin) << " " << wprnh << " " << df->GetBinError(ibin) << endl;

	ibin++;
      }

      line_number++;
    }
    df_text.close();
    
    TCanvas *df_c = new TCanvas("df_c","DF canvas");
    df_c->cd(1); df->Draw("E");
    W_asym_CPbPtdf_all->Divide(W_asym_CPbPt_all,df,1,1,"");
  }
  else {
    cout << "No dilution factor exists!" << endl;
    W_asym_CPbPtdf_all->Add(W_asym_CPbPt_all,1);
  }

  Hlist.Add(hsthe);
  Hlist.Add(hsdel);
  Hlist.Add(hsen);
  Hlist.Add(Q2);
  Hlist.Add(W);
  Hlist.Add(Q2vsW);
  Hlist.Add(W_asym_raw_all);
  Hlist.Add(W_asym_CPbPt_all);
  Hlist.Add(W_asym_CPbPtdf_all);
  Hlist.Add(df);

  //Drawing for check
  TCanvas *c1 = new TCanvas("c1","canvas1");
  c1->Divide(2,2);
  c1->cd(1); hsthe->Draw();
  c1->cd(2); hsen->Draw();
  c1->cd(3); W->Draw();
  c1->cd(4); Q2->Draw();
  TCanvas *c2 = new TCanvas("c2","canvas2");
  c2->cd(1); Q2vsW->SetOption("colz"); Q2vsW->Draw();
  TCanvas *c3 = new TCanvas("c3","canvas3");
  c3->Divide(2,2);
  c3->cd(1); W_p->Draw("E");
  c3->cd(2); W_n->Draw("E");
  c3->cd(3); W_p_p_n->Draw("E");
  c3->cd(4); W_p_m_n->Draw("E");  
  TCanvas *c4 = new TCanvas("c4","canvas4");
  c4->cd(1); W_asym_raw->Draw("E");
  TCanvas *c5 = new TCanvas("c5","canvas5");
  c5->cd(1); W_asym_CPbPt->Draw("E");
  TCanvas *c6 = new TCanvas("c6","canvas6");
  W_asym_raw_all->SetMaximum(0.1);
  W_asym_raw_all->SetMinimum(-0.1);
  c6->cd(1); W_asym_raw_all->Draw("E");
  TCanvas *c7 = new TCanvas("c7","canvas7");
  W_asym_CPbPt_all->SetMaximum(0.2);
  W_asym_CPbPt_all->SetMinimum(-0.2);
  c7->cd(1); W_asym_CPbPt_all->Draw("E");
  TCanvas *c8 = new TCanvas("c8","canvas8");
  W_asym_CPbPtdf_all->SetMaximum(1.0);
  W_asym_CPbPtdf_all->SetMinimum(-1.0);
  c8->cd(1); W_asym_CPbPtdf_all->Draw("E");

  //Writing to file
  TFile *result = new TFile("new.root","RECREATE");
  Hlist.Write();
  result->Close();
			
}


TString filename(Int_t nrun, Int_t index) {

  TString dir, prefix, suffix, s_number;

  dir = "ntup/";
  prefix = "hms";
  suffix = ".root";
  s_number = Form("%d.%d",nrun,index);

  return dir + prefix + s_number + suffix;
}

TString df_filename(Int_t nrun) {

  TString dir, prefix, suffix, s_number;

  dir = "df/";
  prefix = "df_";
  suffix = ".out"; //".root";
  s_number = Form("%d",nrun);

  return dir + prefix + s_number + suffix;
}

// Kinematics from report files
class Run_info {
private:
  Int_t Run;
  Float_t BeamE;
  Float_t HMS_Deadtime;
  Float_t Current;
  Float_t Runtime;
  Float_t P_HMS;
  Float_t Theta_HMS;
  Float_t Charge_help;
  Float_t Charge_helm;
  Int_t hwplate;
  Float_t Pb;
  Float_t Pt;
  Float_t LT_cor;
  Float_t line_elem(TString*, Int_t);
public:
  Run_info(Int_t);
  ~Run_info();
  void read(Int_t);
  void summary();
  Int_t Run() { return Run; }
  Float_t BeamE() { return BeamE; }
  Float_t HMS_Deadtime() { return HMS_Deadtime; }
  Float_t Current() { return Current; }
  Float_t Runtime() { return Runtime; }
  Float_t P_HMS() { return P_HMS; }
  Float_t Theta_HMS() { return Theta_HMS; }
  Int_t hwplate() { return hwplate; }
  Float_t Pb() { return Pb; }
  Float_t Pt() { return Pt; }
  Float_t LT_cor() { return LT_cor; }
  Float_t Charge_help() { return Charge_help; }
  Float_t Charge_helm() { return Charge_helm; }
};

Run_info::Run_info(Int_t nrun) {
  read(nrun);
}

Run_info::~Run_info() { }

Run_info::read(Int_t nrun) {
  // BeamE in MeV
  // HMS_Deadtime as ratio between 0 to 1 (BeamOn)
  // Current in microA (BCM1 BeamOn)
  // Runtime in sec (BCM1 BeamOn)
  // P_HMS in MeV/c
  // Theta_HMS in degrees
  // Charge_help/m in microC
  // hwplate 1/0
  // Pb, Pt in percent
  TString genfile = Form("scalers/gen%d.txt", nrun);
  TString run = Form("%d", nrun);
  TString line;

  Run=nrun;
  // Run infomation from E_Pb file
  ifstream Beaminfo;
  Beaminfo.open("E_Pb.txt");
  while(!Beaminfo.eof()) {
    line.ReadLine(Beaminfo);
    if(line.BeginsWith(run.Data())) {
      cout << "E_Pb : " << line.Data() << endl;
      BeamE = line_elem(line, 1);
      hwplate = (Int_t)line_elem(line, 5);
      Pb = line_elem(line,6); 
    }
  }
  Beaminfo.close();

  // Live time correction from livetime_list file
  ifstream LTinfo;
  LTinfo.open("livetime_list.txt");
  while(!LTinfo.eof()) {
    line.ReadLine(LTinfo);
    if(line.BeginsWith(run.Data())) {
      cout << "LT_correction : " << line.Data() << endl;
      LT_cor = line_elem(line, 3);
    }
  }
  LTinfo.close();

  // Target polarization using 2011/02/09 J. Maxwell charge averaged polarization with 50 nA cut or no cut
  ifstream Ptinfo;
  Ptinfo.open("charge_averaged_runs.nocut.txt");
  while(!Ptinfo.eof()) {
    line.ReadLine(Ptinfo);
    if(line.BeginsWith(run.Data())) {
      cout << "Pt : " << line.Data() << endl;
      Pt = line_elem(line, 1);
    }
  }
  Ptinfo.close();

  ifstream Report;
  Report.open(genfile.Data());
  while(!Report.eof()) {
    line.ReadLine(Report);
    if(line.BeginsWith("HMS_Deadtime")) HMS_Deadtime = line_elem(line,2);
    if(line.BeginsWith("Current")) Current = line_elem(line,2);
    if(line.BeginsWith("Runtime")) Runtime = line_elem(line,2);
    if(line.BeginsWith("P_HMS")) P_HMS = line_elem(line,2)*1e3;
    if(line.BeginsWith("Theta_HMS")) Theta_HMS = line_elem(line,2);
    if(line.BeginsWith("Charge_help")) Charge_help = line_elem(line,2);
    if(line.BeginsWith("Charge_helm")) Charge_helm = line_elem(line,2);
  }
  Report.close();
}

void Run_info::summary() {
  cout << "Summary of Run information" << endl;
  cout << "(Run#, BeamE, HMS_Deadtime, Current, Runtime, P_HMS, Theta_HMS, hwplate, Pb, Pt, LT_correction)" << endl;
  cout << Run << ", ";
  cout << BeamE << ", ";
  cout << HMS_Deadtime << ", ";
  cout << Current << ", ";
  cout << Runtime << ", ";
  cout << P_HMS << ", ";
  cout << Theta_HMS << ", ";
  cout << hwplate << ", ";
  cout << Pb << ", ";
  cout << Pt << ", ";
  cout << LT_cor << endl;
}

Float_t Run_info::line_elem(TString line, Int_t index) {
  TObjArray *sublist = line.Tokenize("\t ");
  TObjString *sub = sublist->At(index);
  return atof(sub->GetString().Data());
}

// Kinematic variables
class Kin_var {
private:
  Float_t sEi;
  Float_t sEf;
  Float_t stheta;
  Float_t stheta_degree;
  Float_t sPhi;
  Float_t sQ2;
  Float_t M;
public:
  Kin_var(Float_t, Float_t, Float_t);
  ~Kin_var();
  void set(Float_t, Float_t, Float_t);
  void summary();
  Float_t Ei() { return sEi; }
  Float_t Ef() { return sEf; }
  Float_t theta() { return stheta_degree; }
  Float_t nu() { return sEi-sEf; }
  Float_t Q2() { return sQ2; }
  Float_t W2();
  Float_t x_bj();
  Float_t gamma();
  Float_t hsca1pa();
  Float_t hsca1pe();
  Float_t hsca2pe();
  Float_t hsna1();
  Float_t hsna2();
};

Kin_var::Kin_var(Float_t Ei,Float_t Ef, Float_t theta) {
  set(Ei,Ef,theta);
}

Kin_var::~Kin_var() { }

void Kin_var::set(Float_t Ei,Float_t Ef, Float_t theta) {
  // Energy in MeV, theta in degree, but stheta in radian
  sEi = Ei;
  sEf = Ef;
  stheta = theta*TMath::Pi()/180;
  stheta_degree = theta;
  sPhi = 0.0;
  sQ2 = 4*sEi*sEf*TMath::Sin(stheta/2)*TMath::Sin(stheta/2);
  M = 955.74;
}

void Kin_var::summary() {
  cout << "Summary of Kinematic variables" << endl;
  cout << "Ei = " << sEi << endl;
  cout << "Ef = " << sEf << endl;
  cout << "theta = " << stheta << endl;
  cout << "nu = " << nu() << endl;
  cout << "Q2 = " << Q2() << endl;
  cout << "W2 = " << W2() << endl;
  cout << "x_bj = " << x_bj() << endl;
  cout << "gamma = " << gamma() << endl;
}

Float_t Kin_var::W2() {
  return M*M + 2*M*nu() - sQ2;
}

Float_t Kin_var::x_bj() {
  return sQ2/(2*M*nu());
}

Float_t Kin_var::gamma() {
  return 2*M*x_bj()/sqrt(sQ2);
}

//# Col 55: HSCA1PA = E-EP*COS(TH)                
//# Col 56: HSCA1PE = EP*SIN(TH)/COS(PHI)         
//# Col 57: HSCA2PE = (E-EP*COS(TH) )/(EP*SIN(TH)COS(PHI)  
//# Col 58: HSNA1   = 1/(E+EP)                             
//# Col 59: HSNA2   = SQRT(Q2)/(2E)

Float_t Kin_var::hsca1pa() {
  return sEi - sEf*TMath::Cos(stheta);
}

Float_t Kin_var::hsca1pe() {
  return sEf*TMath::Sin(stheta)/TMath::Cos(sPhi);
}

Float_t Kin_var::hsca2pe() {
  return (sEi - sEf*TMath::Cos(stheta))/(sEf*TMath::Sin(stheta)*TMath::Cos(sPhi));
}

Float_t Kin_var::hsna1() {
  return 1/(sEi + sEf);
}

Float_t Kin_var::hsna2() {
  return sqrt(sQ2)/(2*sEi);
}