//macro to summarize results on results of <<wrong>> propagation in graphs

void Limits(TString path = "/panda/myResults/RK/HIMster/mom_15/200samp_1e4/"){
  //  const int npbeam = 9;
  //  const int nq=6;
  const int nq=1;
  double Pbeam = 15;

   // const int npbeam = 9;
   // int dpbeam[npbeam] = {-80,-40,-20,-10,0,10,20,40,80};
const int npbeam = 7;
 int dpbeam[npbeam] = {-80,-60,-40,-20,-10,0,10};

  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};
 int th1=0;
  int th2=393;
  int forwAss=0; //0=false,1=true
  //      int forwAss=1; //0=false,1=true
  int xycut=1; //we are interesting only in x-y values in LMD range
  //int xycut=0; //we are interesting only in x-y values in LMD range
  const int ncut=2;
  int thphcut[ncut]={0,1}; //no cut
  double nev[ncut][npbeam];
  double th_min_sim[ncut][npbeam];
  double th_max_sim[ncut][npbeam];
  double th_min_rec[ncut][npbeam];
  double th_max_rec[ncut][npbeam];
 double ph_min_sim[ncut][npbeam];
  double ph_max_sim[ncut][npbeam];
  double ph_min_rec[ncut][npbeam];
  double ph_max_rec[ncut][npbeam];
  double x_min_rec[ncut][npbeam];
  double x_max_rec[ncut][npbeam];
  double y_min_rec[ncut][npbeam];
  double y_max_rec[ncut][npbeam];
 
  //  int thphcut=1; //with cut

      //  for (int ich=0;ich<nq;ich++){ 
      for (int ich=0;ich<ncut;ich++){ 
	for(int im=0;im<npbeam;im++){
	  th_min_sim[ich][im]=0;
	  th_max_sim[ich][im]=0;
	  th_min_rec[ich][im]=0;
	  th_max_rec[ich][im]=0;
	  ph_min_sim[ich][im]=0;
	  ph_max_sim[ich][im]=0;
	  ph_min_rec[ich][im]=0;
	  ph_max_rec[ich][im]=0;
	  x_min_rec[ich][im]=0;
	  x_max_rec[ich][im]=0;
	  
	  int icharge = q[0];
	  int imom = dpbeam[im];
	  
	  
      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[ich];
      momName+=".root";


      TFile *f_out = new TFile(momName,"READ");
      if (f_out->IsZombie()) {
	//cout << "Error opening file" << endl;
	cout<<"File "<<momName<<"corrapted!"<<endl;
	//       exit(-1);
      }
      else{
	cout<<"File "<<momName<<" is OK :)"<<endl;
      }

      TTree *tProp = (TTree*)f_out->Get("tByPropout");
      nev[ich][im] =  (tProp->GetEntries()); //% for 1e4 events
      TH1 *hth_sim = new TH1D("hth_sim","",1e4,th1,th2);
      tProp->Project("hth_sim","1e3*th_sim");//theta sim
      TH1 *hth_rec = new TH1D("hth_rec","",1e4,0,1570);
      tProp->Project("hth_rec","1e3*th_out");//theta rec
      TH1 *hph_sim = new TH1D("hph_sim","",1e4,0,2*TMath::Pi());
      tProp->Project("hph_sim","ph_sim");//phi sim
      TH1 *hph_rec = new TH1D("hph_rec","",1e4,-TMath::Pi(),TMath::Pi());
      tProp->Project("hph_rec","ph_out");//phi rec

      // TH1 *hx_rec = new TH1D("hx_rec","",1e6,-2e2,2e2);
      // tProp->Project("hx_rec","x_out");//X

      TH1 *hx_rec = new TH1D("hx_rec","",1e6,-5e2,5e2);
      tProp->Project("hx_rec","1e4*x_out");//X
      
      TH1 *hy_rec = new TH1D("hy_rec","",1e6,-2e2,2e2);
      tProp->Project("hy_rec","y_out");//Y
      //    f_out->Close();

      /// begin th limit sim
      bool th_min_sim_fl=false; bool th_max_sim_fl=false;
      for(int ibin=0;ibin<1e4;ibin++){
	double thbin = hth_sim->GetBinContent(ibin);
	if(!th_min_sim_fl){
	  if(thbin>0){ 
	    //	    cout<<"thbin["<<ibin<<"]="<<thbin<<endl;
	    th_min_sim[ich][im]=double(ibin*(th2-th1)/1e4);
	    th_min_sim_fl = true;
	  }
	}
	else{ 
	  if(thbin<1 && !th_max_sim_fl){
	    th_max_sim[ich][im]=double((ibin-1)*(th2-th1)/1e4);
	    th_max_sim_fl=true;
	  }
	  if(th_max_sim_fl && thbin>0){ //in case of holes
	    th_max_sim_fl=false;
	  }
	}
      }/// end th limit sim

      /// begin th limit rec
     
      bool th_min_rec_fl=false; bool th_max_rec_fl=false;
      for(int ibin=0;ibin<1e4;ibin++){
	double thbin = hth_rec->GetBinContent(ibin);
	if(!th_min_rec_fl){
	  if(thbin>0){ 
	    //	    cout<<"thbin["<<ibin<<"]="<<thbin<<endl;
	    th_min_rec[ich][im]=double(ibin*1570/1e4);
	    th_min_rec_fl = true;
	  }
	}
	else{ 
	  if(thbin<1 && !th_max_rec_fl){
	    th_max_rec[ich][im]=double((ibin-1)*1570/1e4);
	    th_max_rec_fl=true;
	  }
	  if(th_max_rec_fl && thbin>0){ //in case of holes
	    th_max_rec_fl=false;
	  }
	}
      }/// end th limit rec

      /// begin ph limit sim
      bool ph_min_sim_fl=false; bool ph_max_sim_fl=false;
      for(int ibin=0;ibin<1e4;ibin++){
	double phbin = hph_sim->GetBinContent(ibin);
	if(!ph_min_sim_fl){
	  if(phbin>0){ 
	    //	    cout<<"phbin["<<ibin<<"]="<<phbin<<endl;
	    //	    ph_min_sim[ich][im]=double(ibin*(2*TMath::Pi())/1e4);
	    ph_min_sim[ich][im]=double(ibin*(2*TMath::Pi())/1e4) -TMath::Pi(); //-Pi to Pi
	    ph_min_sim_fl = true;
	  }
	}
	else{ 
	  if(phbin<1 && !ph_max_sim_fl){
	    //	    ph_max_sim[ich][im]=double((ibin-1)*(2*TMath::Pi())/1e4);
	    ph_max_sim[ich][im]=double((ibin-1)*(2*TMath::Pi())/1e4) - TMath::Pi(); //-Pi to Pi
	    ph_max_sim_fl=true;
	  }
	  if(ph_max_sim_fl && phbin>0){ //in case of holes
	    ph_max_sim_fl=false;
	  }
	}
      }/// end ph limit sim

      /// begin ph limit rec
      bool ph_min_rec_fl=false; bool ph_max_rec_fl=false;
      for(int ibin=0;ibin<1e4;ibin++){
	double phbin = hph_rec->GetBinContent(ibin);
	if(!ph_min_rec_fl){
	  if(phbin>0){ 
	    //	    cout<<"phbin["<<ibin<<"]="<<phbin<<endl;
	    ph_min_rec[ich][im]=double(ibin*(2*TMath::Pi())/1e4 - TMath::Pi());
	    ph_min_rec_fl = true;
	  }
	}
	else{ 
	  if(phbin<1 && !ph_max_rec_fl){
	    ph_max_rec[ich][im]=double((ibin-1)*(2*TMath::Pi())/1e4 - TMath::Pi());
	    ph_max_rec_fl=true;
	  }
	  if(ph_max_rec_fl && phbin>0){ //in case of holes
	    ph_max_rec_fl=false;
	  }
	}
      }/// end ph limit rec



      /// begin x limit rec
      bool x_min_rec_fl=false; bool x_max_rec_fl=false;
      for(int ibin=0;ibin<1e6;ibin++){
	double xbin = hx_rec->GetBinContent(ibin);
	if(!x_min_rec_fl){
	  if(xbin>0){ 
	    //	    cout<<"xbin["<<ibin<<"]="<<xbin<<endl;
	    x_min_rec[ich][im]=double(ibin*4e2/1e6-2e2);
	    x_min_rec_fl = true;
	  }
	}
	else{ 
	  if(xbin<1 && !x_max_rec_fl){
	    x_max_rec[ich][im]=double((ibin-1)*4e2/1e6-2e2);
	    x_max_rec_fl=true;
	  }
	  if(x_max_rec_fl && xbin>0){ //in case of holes
	    x_max_rec_fl=false;
	  }
	  if(x_max_rec_fl && xbin>0){ //in case of holes
	    x_max_rec_fl=false;
	  }
	}
      }/// end x limit rec

      /// begin y limit rec
      bool y_min_rec_fl=false; bool y_max_rec_fl=false;
      for(int ibin=0;ibin<1e6;ibin++){
	double ybin = hy_rec->GetBinContent(ibin);
	if(!y_min_rec_fl){
	  if(ybin>0){ 
	    //    cout<<"ybin["<<ibin<<"]="<<ybin<<endl;
	    y_min_rec[ich][im]=double(ibin*4e2/1e6-2e2);
	    y_min_rec_fl = true;
	  }
	}
	else{ 
	  if(ybin<1 && !y_max_rec_fl){
	    y_max_rec[ich][im]=double((ibin-1)*4e2/1e6-2e2);
	    y_max_rec_fl=true;
	  }
	  if(y_max_rec_fl && ybin>0){ //in case of holes
	    y_max_rec_fl=false;
	  }
	  if(y_max_rec_fl && ybin>0){ //in case of holes
	    y_max_rec_fl=false;
	  }
	}
      }/// end y limit rec

      cout<<"x_rec_min = "<<x_min_rec[ich][im]<<" "<<"x_rec_max = "<<x_max_rec[ich][im]<<endl;
      cout<<"y_rec_min = "<<y_min_rec[ich][im]<<" "<<"y_rec_max = "<<y_max_rec[ich][im]<<endl;
      cout<<"th_rec_min = "<<th_min_rec[ich][im]<<" "<<"th_rec_max = "<<th_max_rec[ich][im]<<endl;
      cout<<"th_sim_min = "<<th_min_sim[ich][im]<<" "<<"th_sim_max = "<<th_max_sim[ich][im]<<endl;
      cout<<"ph_rec_min = "<<ph_min_rec[ich][im]<<" "<<"ph_rec_max = "<<ph_max_rec[ich][im]<<endl;
      cout<<"ph_sim_min = "<<ph_min_sim[ich][im]<<" "<<"ph_sim_max = "<<ph_max_sim[ich][im]<<endl;
    }
  }

  //manipulation with results for more clear plots (draw result band as an uncertanty band)
  double th_mn_sim[ncut][npbeam], th_un_sim[ncut][npbeam];
  double th_mn_rec[ncut][npbeam], th_un_rec[ncut][npbeam];
  double ph_mn_sim[ncut][npbeam], ph_un_sim[ncut][npbeam];
  double x_mn_rec[ncut][npbeam], x_un_rec[ncut][npbeam];
  double err_x_sysm[ncut][npbeam];

  // //q=-1, 1.5 GeV/c
  // err_x_sysm[0][0]=-79.2561;
  // err_x_sysm[0][1]=-32.0781;
  // err_x_sysm[0][2]=-24.3141;
  // err_x_sysm[0][3]=-21.6898;
  // err_x_sysm[0][4]=-19.5718;
  // err_x_sysm[0][5]=-17.8302;
  // err_x_sysm[0][6]=-16.3701;
  // err_x_sysm[0][7]=-14.0633;
  // err_x_sysm[0][8]=-10.9661;

  // //  //q = 1, 1.5 GeV/c
  // // err_x_sysm[0][0]=121.43;
  // // err_x_sysm[0][1]=31.7415;
  // // err_x_sysm[0][2]=24.1011;
  // // err_x_sysm[0][3]=21.5176;
  // // err_x_sysm[0][4]=19.4307;
  // // err_x_sysm[0][5]=17.7121;
  // // err_x_sysm[0][6]=16.2704;
  // // err_x_sysm[0][7]=13.99;
  // // err_x_sysm[0][8]=10.921;

if(q[0]<0){
  // q = -1, 15 GeV/c
  err_x_sysm[0][0]=191.253;
  err_x_sysm[0][1]=-90.8542;
  err_x_sysm[0][2]=-132.717;
  err_x_sysm[0][3]=-123.771;
  err_x_sysm[0][4]=-116.459;
  err_x_sysm[0][5]=-109.723;
  err_x_sysm[0][6]=-102.672;
    }
 else{
   //q = 1, 15 GeV/c
   err_x_sysm[0][0]=0;
   err_x_sysm[0][1]=357.051;
   err_x_sysm[0][2]=252.462;
   err_x_sysm[0][3]=191.409;
   err_x_sysm[0][4]=170.383;
   err_x_sysm[0][5]=153.058;
   err_x_sysm[0][6]=138.528;
 }
  for (int ich=0;ich<ncut;ich++){ 

    for(int im=0;im<npbeam;im++){
      th_mn_sim[ich][im] = th_min_sim[ich][im]+ 0.5*(th_max_sim[ich][im] - th_min_sim[ich][im]);
      th_un_sim[ich][im] = 0.5*(th_max_sim[ich][im] - th_min_sim[ich][im]);

      ph_mn_sim[ich][im] = ph_min_sim[ich][im]+ 0.5*(ph_max_sim[ich][im] - ph_min_sim[ich][im]);
      ph_un_sim[ich][im] = 0.5*(ph_max_sim[ich][im] - ph_min_sim[ich][im]);
      //      cout<<"th_mn_sim[ich][im]  = "<<th_mn_sim[ich][im]<<" th_un_sim[ich][im] = "<<th_un_sim[ich][im]<<endl;

      th_mn_rec[ich][im] = th_min_rec[ich][im]+ 0.5*(th_max_rec[ich][im] - th_min_rec[ich][im]);
      th_un_rec[ich][im] = 0.5*(th_max_rec[ich][im] - th_min_rec[ich][im]);
      x_mn_rec[ich][im] = x_min_rec[ich][im]+ 0.5*(x_max_rec[ich][im] - x_min_rec[ich][im]);
      //      x_mn_rec[ich][im] = x_min_rec[ich][im]+ 0.5*(x_max_rec[ich][im] - x_min_rec[ich][im]) + 1e-4*err_x_sysm[ich][im];
      //    x_un_rec[ich][im] = 0.5*(x_max_rec[ich][im] - x_min_rec[ich][im]))+1e-4*fabs(err_x_sysm[ich][im]);
      x_un_rec[ich][im] = 0.5*(x_max_rec[ich][im] - x_min_rec[ich][im]))+fabs(err_x_sysm[0][im]);
      cout<<"x_mn_rec[ich][im]  = "<<x_mn_rec[ich][im]<<" x_un_rec[ich][im] = "<<x_un_rec[ich][im]<<endl;
    } 
  }

  TGraph *gnev[ncut];
  TMultiGraph *mgrev = new TMultiGraph();

  //  TGraph *gth_sim_min[ncut];   TGraph *gth_sim_max[ncut];
   TGraphErrors *gth_sim_mn[ncut]; 
   
   TMultiGraph *mgr_th_sim = new TMultiGraph();   

   // TGraph *gph_sim_min[ncut];   TGraph *gph_sim_max[ncut];
   TGraphErrors *gph_sim_mn[ncut]; 
  TMultiGraph *mgr_ph_sim = new TMultiGraph();   


  //  TGraph *gth_rec_min[ncut];   TGraph *gth_rec_max[ncut];
  TGraphErrors *gth_rec_mn[ncut]; 
  TMultiGraph *mgr_th_rec = new TMultiGraph();   

  // TGraph *gx_rec_min[ncut];   TGraph *gx_rec_max[ncut];
  TGraphErrors *gx_rec_mn[ncut]; 
  TMultiGraph *mgr_x_rec = new TMultiGraph();   

  TLegend *leg = new TLegend(0.64,0.65,0.98,0.98);
  leg->SetFillColor(0);
  leg->SetTextFont(42);
  leg->SetTextSize(0.05);
  for (int ich=0;ich<ncut;ich++){ 

    gnev[ich] = new TGraph(npbeam,pbeamval,nev[ich]);
    // gnev[ich]->SetMarkerStyle(24+ich);
    // gnev[ich]->SetMarkerSize(3.5);

    if(q[0]<0){
      gnev[ich]->SetLineColor(kOrange+1+ich);
      //      gnev[ich]->SetMarkerColor(kOrange+1+ich);
    }
    else{
      gnev[ich]->SetLineColor(kTeal+2+ich);
      //      gnev[ich]->SetMarkerColor(kTeal+2+ich);
    }
    gnev[ich]->SetLineWidth(2);
    mgrev->Add(gnev[ich], "LP");
    TString grname = "q_{[e]} = ";
    grname +=q[0];
    leg->SetHeader(grname);
    //  leg->AddEntry(gnev[ich],grname,"l");

    // gth_sim_min[ich] = new TGraph(npbeam,pbeamval,th_min_sim[ich]);
    // gth_sim_max[ich] = new TGraph(npbeam,pbeamval,th_max_sim[ich]);
    gth_sim_mn[ich] = new TGraphErrors(npbeam,pbeamval,th_mn_sim[ich],0,th_un_sim[ich]);
    // gph_sim_min[ich] = new TGraph(npbeam,pbeamval,ph_min_sim[ich]);
    // gph_sim_max[ich] = new TGraph(npbeam,pbeamval,ph_max_sim[ich]);
    gph_sim_mn[ich] = new TGraphErrors(npbeam,pbeamval,ph_mn_sim[ich],0,ph_un_sim[ich]);

    // gth_rec_min[ich] = new TGraph(npbeam,pbeamval,th_min_rec[ich]);
    // gth_rec_max[ich] = new TGraph(npbeam,pbeamval,th_max_rec[ich]); 
    gth_rec_mn[ich] = new TGraphErrors(npbeam,pbeamval,th_mn_rec[ich],0,th_un_rec[ich]);

    gx_rec_mn[ich] = new TGraphErrors(npbeam,pbeamval,x_mn_rec[ich],0,x_un_rec[ich]);
    gx_rec_mn[ich]->SetMarkerStyle(20);
    // gx_rec_min[ich] = new TGraph(npbeam,pbeamval,x_min_rec[ich]);
    // gx_rec_max[ich] = new TGraph(npbeam,pbeamval,x_max_rec[ich]);

    // gth_sim_min[ich]->SetMarkerStyle(24+ich);
    // gth_sim_min[ich]->SetMarkerSize(3.5);
    // gph_sim_min[ich]->SetMarkerStyle(24+ich);
    // gph_sim_min[ich]->SetMarkerSize(3.5);
    // gth_rec_min[ich]->SetMarkerStyle(24+ich);
    // gth_rec_min[ich]->SetMarkerSize(3.5);
    // gx_rec_min[ich]->SetMarkerStyle(24+ich);
    // gx_rec_min[ich]->SetMarkerSize(3.5);

    // gth_sim_max[ich]->SetMarkerStyle(24+ich);
    // gth_sim_max[ich]->SetMarkerSize(3.5);
    // gph_sim_max[ich]->SetMarkerStyle(24+ich);
    // gph_sim_max[ich]->SetMarkerSize(3.5);
    // gth_rec_max[ich]->SetMarkerStyle(24+ich);
    // gth_rec_max[ich]->SetMarkerSize(3.5);
    // gx_rec_max[ich]->SetMarkerStyle(24+ich);
    // gx_rec_max[ich]->SetMarkerSize(3.5);

    if(q[0]<0){
      gth_sim_mn[ich]->SetLineColor(kOrange+1+ich);
      gth_sim_mn[ich]->SetFillStyle(3004+ich);
      gth_sim_mn[ich]->SetFillColor(kOrange+1+ich);

      gph_sim_mn[ich]->SetLineColor(kOrange+1+ich);
      gph_sim_mn[ich]->SetFillStyle(3004+ich);
      gph_sim_mn[ich]->SetFillColor(kOrange+1+ich);
    

      gth_rec_mn[ich]->SetLineColor(kOrange+1+ich);
      gth_rec_mn[ich]->SetFillStyle(3004+ich);
      gth_rec_mn[ich]->SetFillColor(kOrange+1+ich);

      gx_rec_mn[ich]->SetLineColor(kOrange+1+ich);
      //      gx_rec_mn[ich]->SetFillStyle(3004+ich);
      gx_rec_mn[ich]->SetFillColor(kOrange+1+ich);

   
    }
    else{
      gth_sim_mn[ich]->SetLineColor(kTeal+2+ich);
      gth_sim_mn[ich]->SetFillStyle(3004+ich);
      gth_sim_mn[ich]->SetFillColor(kTeal+2+ich);

      gph_sim_mn[ich]->SetLineColor(kTeal+2+ich);
      gph_sim_mn[ich]->SetFillStyle(3004+ich);
      gph_sim_mn[ich]->SetFillColor(kTeal+2+ich);
    

      gth_rec_mn[ich]->SetLineColor(kTeal+2+ich);
      gth_rec_mn[ich]->SetFillStyle(3004+ich);
      gth_rec_mn[ich]->SetFillColor(kTeal+2+ich);

      gx_rec_mn[ich]->SetLineColor(kTeal+2+ich);
      //      gx_rec_mn[ich]->SetFillStyle(3004+ich);
      gx_rec_mn[ich]->SetFillColor(kTeal+2+ich);


      // gth_sim_min[ich]->SetMarkerColor(kTeal+2+ich);
      // gth_sim_max[ich]->SetMarkerColor(kTeal+2+ich);
      // gth_sim_min[ich]->SetLineColor(kTeal+2+ich);
      // gth_sim_max[ich]->SetLineColor(kTeal+2+ich);
      // gph_sim_min[ich]->SetMarkerColor(kTeal+2+ich);
      // gph_sim_max[ich]->SetMarkerColor(kTeal+2+ich);
      // gph_sim_min[ich]->SetLineColor(kTeal+2+ich);
      // gph_sim_max[ich]->SetLineColor(kTeal+2+ich);
      // gth_rec_min[ich]->SetMarkerColor(kTeal+2+ich);
      // gth_rec_max[ich]->SetMarkerColor(kTeal+2+ich);
      // gth_rec_min[ich]->SetLineColor(kTeal+2+ich);
      // gth_rec_max[ich]->SetLineColor(kTeal+2+ich);
      // gx_rec_min[ich]->SetMarkerColor(kTeal+2+ich);
      // gx_rec_max[ich]->SetMarkerColor(kTeal+2+ich);
      // gx_rec_min[ich]->SetLineColor(kTeal+2+ich);
      // gx_rec_max[ich]->SetLineColor(kTeal+2+ich);
    }
    // mgr_th_sim->Add(gth_sim_min[ich],"C");
    mgr_th_sim->Add(gth_sim_mn[ich],"3");
    //mgr_th_sim->Add(gth_sim_max[ich],"C");
    // mgr_ph_sim->Add(gph_sim_min[ich],"PL");
    // mgr_ph_sim->Add(gph_sim_max[ich],"PL");
    mgr_ph_sim->Add(gph_sim_mn[ich],"3");
    mgr_th_rec->Add(gth_rec_mn[ich],"3");
    // mgr_th_rec->Add(gth_rec_min[ich],"PL");
    // mgr_th_rec->Add(gth_rec_max[ich],"PL");
    mgr_x_rec->Add(gx_rec_mn[ich],"3");
    // mgr_x_rec->Add(gx_rec_min[ich],"PL");
    // mgr_x_rec->Add(gx_rec_max[ich],"PL");
    //   grname +=", ";
    TString cutn;
    if(thphcut[ich]<1) cutn ="(x,y)";
    else cutn ="(x,y) & (#hat{#theta},#hat{#phi})";
    cutn +=" cut";
    leg->AddEntry(gth_sim_mn[ich],cutn,"f");//field
    //    leg->AddEntry(gx_rec_mn[ich],cutn,"ple");//line
  }
  //  mgrx->Draw("AP");
//  mgrev->Draw("AP");

  // // //draw th_sim
  // mgr_th_sim->Draw("AC");
  // mgr_th_sim->GetXaxis()->SetTitle("P, GeV/c");
  // mgr_th_sim->GetYaxis()->SetTitle("#theta_{sim}, mrad");
  // leg->Draw();
  // // /// end draw th_sim

  // // //draw ph_sim
  // mgr_ph_sim->Draw("AC");
  // mgr_ph_sim->GetXaxis()->SetTitle("P, GeV/c");
  // mgr_ph_sim->GetYaxis()->SetTitle("#phi_{sim}, rad");
  // leg->Draw();
  // // /// end draw ph_sim


  // //draw th_rec
  // mgr_th_rec->Draw("AC");
  // mgr_th_rec->GetXaxis()->SetTitle("P, GeV/c");
  // mgr_th_rec->GetYaxis()->SetTitle("#theta_{rec}, mrad");
  // leg->Draw();
  // /// end draw th_rec

  //draw x_rec
  mgr_x_rec->Draw("AC");
  mgr_x_rec->GetXaxis()->SetTitle("P, GeV/c");
//  mgr_x_rec->GetYaxis()->SetTitle("x_{rec}, cm");
  mgr_x_rec->GetYaxis()->SetTitle("x_{rec}, #mum");
  leg->Draw();
  /// end draw x_rec

  // 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*");

}