// macro for bkg channels analysis 
//
// author: a.karavdina
// date: 07 January 2014

void AnaBkg(TString path="/panda/myResults/DPMoutputDec2013/FTF/NewGeometry/NoCut/mom_15/"){
  const int nfiles=3;
  //No cuts
  //  TString dates1[nfiles] = {"01032014","02032014","03032014"}; //el+inel, 
  //  TString dates1[nfiles] = {"04032014","05032014","06032014"}; //el+inel, 
  //    TString dates1[nfiles] = {"01022014","04022014","07022014"}; //inel, 
  //  TString dates1[nfiles] = {"02022014","05022014","08022014"}; //el+inel, 
  //    TString dates1[nfiles] = {"03022014","06022014","09022014"}; //el 
  TString dates1[nfiles] = {"01042014","02042014","03042014"}; //el+inel, 

    //X&Y cuts
  //TString dates1[nfiles] = {"10022014","13022014","16022014"}; //inel 15 GeV
  //  TString dates1[nfiles] = {"11022014","14022014","17022014"}; //el+inel 15 GeV
  // TString dates1[nfiles] = {"12022014","15022014","18022014"}; //el 15 GeV
// const int nfiles=2;
//  TString dates1[nfiles] = {"14022014","17022014"}; //el+inel 15 GeV
  //$\sigma^{1.5,4.6,8.9,11.91,15}_{\theta}$=601,212,100,83,74 $\mu$rad
 
  // //4.06 GeV/c
  //  const double thRes = 0.212 ;//mrad, near IP
      //const double cs_tot = 25.07;// mb,  el
   //  const double cs_tot = 49.44;// mb,  inel
  // const double cs_tot = 74.51;// mb,  inel+el

  // // 8.9 GeV/c
    //       const double thRes = 0.1;//mrad, near IP
  //     const double cs_tot = 14.05;// mb,  el
   //   const double cs_tot = 43.56;// mb,  inel
    //	const double cs_tot = 57.61;// mb,  inel+el


  // 11.91 GeV/c
    //   const double thRes = 0.083;//mrad, near IP
  //   const double cs_tot = 11.21;// mb,  el
     //   const double cs_tot = 41.78;// mb,  inel
    //   const double cs_tot = 52.99;// mb,  inel+el

  // ///// 15GeV/c
  const double thRes = 7.40773e-02 ;//mrad 15 GeV/c, near IP
  const double cs_tot = 50.26;// mb, 15 GeV/c, el+inel
    //    const double cs_tot = 9.76;// mb, 15 GeV/c, el
  //   // const double cs_tot = 40.50;// mb, 15 GeV/c, inel

 ////// 1.5GeV/c

  //  const double thRes = 0.6 ;//mrad 15 GeV/c, near IP
  //     const double cs_tot = 124.63;// mb, 1.5 GeV/c, el+inel
    //   const double cs_tot = 61.92;// mb, 1.5 GeV/c, el
    //   const double cs_tot = 62.71;// mb, 1.5 GeV/c, inel

 
  double bkg_tot1[nfiles];   double bkg_tot2[nfiles]; double bkg_tot3[nfiles]; double bkg_tot4[nfiles]; 
  double sig_tot1[nfiles]; 
  double sig_sec1[nfiles];
  //  double norm = 100./1e7; //normalisation to tot number of event, gives result in %
  double norm = 1e3*cs_tot/1e7; //normalisation to tot number of event, gives result in mub
  // double norm = 1; //no normalisation
  for (int i1=0;i1<nfiles;i1++){
    TString filename1 = path;
    filename1 += dates1[i1];
    filename1 +="/TOT_Lumi_BkgSig.root";
    TFile *fin1 = new TFile(filename1,"READ");
    TTree *tbkg1 = (TTree *)fin1->Get("tBkg"); //Bkg data
    TTree *tsig1 = (TTree *)fin1->Get("tSig"); //Signal data
    // // //No mom cut ----------------------------
    bkg_tot1[i1] = (tbkg1->GetEntries("sumid!=4424 && issecond<0"));
    bkg_tot3[i1] = (tbkg1->GetEntries("sumid!=4424 && issecond>=0"));
    bkg_tot2[i1] = (tbkg1->GetEntries("sumid==4424"));
    TString thcond = "abs(thmc-threc)>";
    thcond +=3*thRes;
    //    thcond += " && sumid==4424";
    bkg_tot4[i1] = (tsig1->GetEntries(thcond.Data()));
    TString thcond2 = "abs(thmc-threc)<=";
    thcond2 +=3*thRes;
    //    thcond2 += " && sumid==4424";
    //   cout<<"  bkg_tot1 = "<< bkg_tot1[i1] <<"  bkg_tot2 = "<< bkg_tot2[i1] <<endl;

    sig_tot1[i1] = (tsig1->GetEntries(thcond2.Data()));
    cout<<"  sig_tot1 = "<< sig_tot1[i1]<<endl;
    // // //[END] No mom cut --------------------

    // // //  cout<<"  sig_tot1 = "<< sig_tot1[i1]<<endl;
    // // // //add mom cut:
    // bkg_tot1[i1] = (tbkg1->GetEntries("sumid!=4424 && issecond<0 && abs(15.0-prec)/15.0<3e-4"));
    // bkg_tot3[i1] = (tbkg1->GetEntries("sumid!=4424 && issecond>=0 && abs(15.0-prec)/15.0<3e-4"));
    // bkg_tot2[i1] = (tbkg1->GetEntries("sumid==4424 && abs(15.0-prec)/15.0<3e-4"));

    // TString thcond = "abs(15.0-prec)/15.0<5e-4 && abs(thmc-threc)>";
    // thcond +=2*thRes;
    // //    thcond += " && sumid==4424";
    // bkg_tot4[i1] = (tsig1->GetEntries(thcond.Data()));
    // TString thcond2 = "abs(15.0-prec)/15.0<3e-4 && abs(thmc-threc)<=";
    // thcond2 +=2*thRes;
    // //    thcond2 += " && sumid==4424";
    // //    cout<<"  bkg_tot1 = "<< bkg_tot1[i1] <<"  bkg_tot2 = "<< bkg_tot2[i1] <<endl;
    // sig_tot1[i1] = (tsig1->GetEntries(thcond2.Data()));
    // // // //   fin->Close();
  }
  double bkg_tot1_mean = 0,bkg_tot1_rms = 0;
  double bkg_tot2_mean = 0,bkg_tot2_rms = 0;
  double bkg_tot3_mean = 0,bkg_tot3_rms = 0;
  double bkg_tot4_mean = 0,bkg_tot4_rms = 0;
  double sig_tot1_mean = 0,sig_tot1_rms = 0;
  for(int i1=0;i1<nfiles;i1++){
    bkg_tot1_mean += bkg_tot1[i1];
    bkg_tot2_mean += bkg_tot2[i1];
    bkg_tot3_mean += bkg_tot3[i1];
    bkg_tot4_mean += bkg_tot4[i1];
    sig_tot1_mean += sig_tot1[i1];
  }
  bkg_tot1_mean /=nfiles;
  bkg_tot2_mean /=nfiles;
  bkg_tot3_mean /=nfiles;
  bkg_tot4_mean /=nfiles;
  sig_tot1_mean /=nfiles;
 for(int i1=0;i1<nfiles;i1++){
   bkg_tot1_rms += (bkg_tot1_mean-bkg_tot1[i1])*(bkg_tot1_mean-bkg_tot1[i1]);
   bkg_tot2_rms += (bkg_tot2_mean-bkg_tot2[i1])*(bkg_tot2_mean-bkg_tot2[i1]);
   bkg_tot3_rms += (bkg_tot3_mean-bkg_tot3[i1])*(bkg_tot3_mean-bkg_tot3[i1]);
   bkg_tot4_rms += (bkg_tot4_mean-bkg_tot4[i1])*(bkg_tot4_mean-bkg_tot4[i1]);
   sig_tot1_rms += (sig_tot1_mean-sig_tot1[i1])*(sig_tot1_mean-sig_tot1[i1]);
  }
 bkg_tot1_rms /=nfiles;
 bkg_tot1_rms = sqrt(bkg_tot1_rms);
 bkg_tot2_rms /=nfiles;
 bkg_tot2_rms = sqrt(bkg_tot2_rms);
 bkg_tot3_rms /=nfiles;
 bkg_tot3_rms = sqrt(bkg_tot3_rms);
 bkg_tot4_rms /=nfiles;
 bkg_tot4_rms = sqrt(bkg_tot4_rms);
 sig_tot1_rms /=nfiles;
 sig_tot1_rms = sqrt(sig_tot1_rms);

 //normalisation to tot.sim events
 bkg_tot1_mean *=norm;
 bkg_tot1_rms *=norm;
 bkg_tot2_mean *=norm;
 bkg_tot2_rms *=norm;
 bkg_tot3_mean *=norm;
 bkg_tot3_rms *=norm;
 bkg_tot4_mean *=norm;
 bkg_tot4_rms *=norm;
 sig_tot1_mean *=norm;
 sig_tot1_rms *=norm;

  // cout<<"El+Inel: "<<endl;
 cout<<"Effective El Prim cross-section (Dth<=3Sigma)"<<sig_tot1_mean<<" +/- "<<sig_tot1_rms<<" mub"<<endl;
 cout<<"Effective El Prim cross-section (Dth>3Sigma)"<<bkg_tot4_mean<<" +/- "<<bkg_tot4_rms<<" mub"<<endl;
 cout<<"Effective El Sec cross-section "<<bkg_tot2_mean<<" +/- "<<bkg_tot2_rms<<" mub"<<endl;
 cout<<"Effective InEl Prim cross-section "<<bkg_tot1_mean<<" +/- "<<bkg_tot1_rms<<" mub"<<endl;
 cout<<"Effective InEl Sec cross-section "<<bkg_tot3_mean<<" +/- "<<bkg_tot3_rms<<" mub"<<endl;
}