//macro to summarize results on forward-backward systematic in graphs

void Systematic(TString path = "/panda/myResults/RK/HIMster/mom_15/200samp_1e4/"){
//void Systematic(TString path = "/panda/karavdina/RungeKuttaSimplified"){
  const int nq=2;
  double Pbeam = 15;
 const int npbeam = 7;
 // int dpbeam[npbeam] = {-80,-40,-20,-10,0,10};
 int dpbeam[npbeam] = {-80,-60,-40,-20,-10,0,10};
  // const int npbeam = 9;
  //  int dpbeam[npbeam] = {-80,-40,-20,-10,0,10,20,40,80};
  double pbeamval[npbeam];
  for(int ip=0;ip<npbeam;ip++){
    pbeamval[ip]=Pbeam+dpbeam[ip]*Pbeam/100.;
  }
  //double q[nq]={-3,-2,-1,1,2,3};
  double q[nq]={-1,1};
  double x_mean[nq][npbeam];
  double x_rms[nq][npbeam];

  for (int ich=0;ich<nq;ich++){ 
    for(int im=0;im<npbeam;im++){
      int icharge = q[ich];
      int imom = dpbeam[im];
      // int th1=0;
      // int th2=785; 
      const int th1=0;   const int th2=393;
      //      const int th1=0;   const int th2=1570;
      int forwAss=1; //0=false,1=true
      int xycut=1; //we are interesting only in x-y values in LMD range
      //     int xycut=0;
      int thphcut=0;
      TString momName = path+"/backward/Prop_By_th_";
      momName+=th1;
      momName+="_";
      momName+=th2;
      momName+="_mrad_";
      momName +=Pbeam;
      momName +="_PLUS_";
      momName+=imom;
      momName+="procent_q_";
      momName+=icharge;
      momName+="_isFfB_";
      momName+=forwAss;
      momName+="_isXYcut_";
      momName+=xycut;
      momName+="_isTHPHcut_";
      momName+=thphcut;
      momName+=".root";



      TFile *f_out = new TFile(momName,"READ");
      TTree *tProp = (TTree*)f_out->Get("tByPropout");
      TH1 *hx = new TH1D("hx","",6e6,-3e6,3e6);
      //  tProp->Project("hx","1e4*x_out","x_in>15 && x_in<35");//X
      tProp->Project("hx","1e4*x_out");//X
      //      tProp->Project("hx","1e4*y_out");//Y
      //      x_mean[ich][im] = fabs(hx->GetMean());  //ABS!
      x_mean[ich][im] = hx->GetMean(); 
      x_rms[ich][im] = hx->GetRMS();
      //x_rms[ich][im] = 0;
      //  hx->Draw();
    }
  }

  TGraphErrors *grsysx[nq]; 
  TMultiGraph *mgrx = new TMultiGraph();
  TLegend *leg = new TLegend(0.8,0.65,0.98,0.98);
  leg->SetFillColor(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.05);
  for (int ich=0;ich<nq;ich++){ 
    grsysx[ich] = new TGraphErrors(npbeam,pbeamval,x_mean[ich],0,x_rms[ich]);
    grsysx[ich]->SetMarkerStyle(20+ich);
    grsysx[ich]->SetMarkerSize(2.5);
    if(q[ich]<0) grsysx[ich]->SetMarkerColor(kRed+ich);
    else grsysx[ich]->SetMarkerColor(kGreen-ich);
    mgrx->Add(grsysx[ich], "P");
    TString grname = "q = ";
    grname +=q[ich];
    leg->AddEntry(grsysx[ich],grname,"lep");
  }
  mgrx->Draw("AP");
  leg->Draw();
  cout<<"q = "<<q[0]<<endl;
  grsysx[0]->Print();
  cout<<""<<endl;
  cout<<"q = "<<q[1]<<endl;
  grsysx[1]->Print();
  // TGraphErrors *grsysx = new TGraphErrors(npbeam,pbeamval,x_mean[0],0,x_rms[0]);
  // grsysx->Draw("A*");

  // TGraphErrors *grsysx = new TGraphErrors(npbeam,pbeamval,x_mean[1],0,x_rms[1]);
  // grsysx->Draw("A*");

}