#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TString.h>
#include <TCanvas.h>
#include <TNtuple.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <TLegend.h>
#include <TMath.h>
#include <sstream>
using namespace std;
int ModulesStat(TString pathG="/home/karavdina/soft/pandaroot/macro/lmd/testPixelAlignSIM/mom_15/")
{

  //  const int nDt = 8;
  // double v_tr_sc[nDt]={0,50,100,200,300,400,500,600};//misalign, trans
   const int nDt = 1;
  double v_tr_sc[nDt]={100};//misalign, trans


  //  const int nMod = 10;
  const int nMod = 1;//TEST
  // const int nDt = 2;
  // double v_tr_sc[nDt]={0,600};//misalign, trans

  // const int nDa = 6;
  // double v_rt_sc[nDa]={0,1,3,5,7,9};//misalign, trans
  // const int nDa = 4; //TEST
  // double v_rt_sc[nDa]={0,1,3,5};//misalign, trans
  const int nDa = 1; //TEST
  double v_rt_sc[nDa]={0};//misalign, trans
  const int nSdGL=100;//number of data samples
  const int nSd=10;//number of data samples sim in step C
  const int nSdstep = 10; //step size for "number of used trks in Knossos"
  const int nuSd = nSdGL/nSdstep;
  //  cout<<"nuSd = "<<nuSd<<endl;
  // const int nSc=10;//number of align.consts samples
  const int nSc=1;//number of align.consts samples TEST
  for(int iSc=0;iSc<nSc;iSc++){//diff output files
    for(int iDa=0;iDa<nDa;iDa++){//diff output files
  // for(int iSc=0;iSc<1;iSc++){//diff output files //TEST
    //    for(int iDa=0;iDa<1;iDa++){//diff output files //TEST

      double ev_mod_b[10][nSd][nuSd+1][nDt]; //module, data sample, sum of data samples, dt
      double ev_mod_a[10][nSd][nuSd+1][nDt]; //module, data sample, sum of data samples, dt
      //mean values of residuals
      double resx_mod_b[10][nSd][nuSd+1][nDt]; //module, data sample, sum of data samples, dt
      double resx_mod_a[10][nSd][nuSd+1][nDt]; //module, data sample, sum of data samples, dt
      double resy_mod_b[10][nSd][nuSd+1][nDt]; //module, data sample, sum of data samples, dt
      double resy_mod_a[10][nSd][nuSd+1][nDt]; //module, data sample, sum of data samples, dt
      //How to save data
      TString resname= pathG+"/AlignmentLimits_EvModStat";
      resname+=v_rt_sc[iDa];
      resname+="mrad_";
      resname+=iSc;
      resname+="constsSample";
      TString resname_pdf = resname+".pdf";
      TString resname_pdf_o = resname_pdf+"(";
      TString resname_pdf_c = resname_pdf+")"; 

	for(int idt=0;idt<nDt;idt++){//loop over dt values (x-axis of plots)
	  for(int imod=0;imod<nMod;imod++){// loop over modules (diff plots)
	    //   for(int isum=1;isum<nSd+1;isum++){ //loop over diff stat.samples (plot with diff markers)
	    for(int ijsum=1;ijsum<nuSd+1;ijsum++){ //loop over diff stat.samples (plot with diff markers)
	      int isum = 1+nSdstep*(ijsum-1);
	      //  cout<<"isum = "<<isum<<endl;
		for(int ipart=0;ipart<nSd;ipart++){//loop over samples (take average)
		  //	  cout<<"Mod#"<<imod<<endl;
	    double tr_sc = v_tr_sc[idt];
	    double rt_sc = v_rt_sc[iDa];
	    //Where search for files
	    TString path = pathG+"/";
	    path+=tr_sc;
	    path+="mkm_";
	    path+=rt_sc;
	    path+="mrad/";
	    path+=iSc;
	    path+="Sum";
	    path+=isum; //sample for alignment with sum of samples
	    path+="Samples/Sample";
	    path +=ipart;
	    path+="/results/";

	    /// Hits hists --------------------------------------------------------------------------------------  
	    TString name1 = path + "/hitsRes_before_alignment_sector";
	    name1+=imod;
	    name1+=".root";
	    TString name2 = path + "/hitsRes_after_alignment_sector"; 
	    name2+=imod;
	    name2+=".root";
	    TFile *f1 = new TFile(name1,"READ");
	    TNtuple *nb = (TNtuple*)f1->Get("nhits");
	    // nb->Print();
	    TH1I *ntrksb = (TH1I*)f1->Get("trks");
	    ev_mod_b[imod][ipart][ijsum][idt] = 1e-3*double(ntrksb->GetEntries());
	    TH1D *hresx_b = new TH1D("hresx_b","",1e3,-1.,1.);
	    nb->Project("hresx_b","xrec-xmc");
	    resx_mod_b[imod][ipart][ijsum][idt] = 1e4*(hresx_b->GetMean());
	    TH1D *hresy_b = new TH1D("hresy_b","",1e3,-1.,1.);
	    nb->Project("hresy_b","yrec-ymc");
	    resy_mod_b[imod][ipart][ijsum][idt] = 1e4*(hresy_b->GetMean());
	    TFile *f2 = new TFile(name2,"READ");
	    TNtuple *na = (TNtuple*)f2->Get("nhits");
	    //	    ev_mod_a[imod][ipart][ijsum][idt] = 1e-3*double(na->GetEntries());
	    TH1I *ntrksa = (TH1I*)f2->Get("trks");
	    ev_mod_a[imod][ipart][ijsum][idt] = 1e-3*double(ntrksa->GetEntries());
	    TH1D *hresx_a = new TH1D("hresx_a","",1e3,-1.,1.);
	    na->Project("hresx_a","xrec-xmc");
	    resx_mod_a[imod][ipart][ijsum][idt] = 1e4*(hresx_a->GetMean());
	    TH1D *hresy_a = new TH1D("hresy_a","",1e3,-1.,1.);
	    na->Project("hresy_a","yrec-ymc");
	    resy_mod_a[imod][ipart][ijsum][idt] = 1e4*(hresy_a->GetMean());
	    cout<<"imod:"<<imod<<" dt:"<<idt<<" isum:"<<isum<<" ipart:"<<ipart<<" nEv:"<<ev_mod_b[imod][ipart][ijsum][idt]<<" "<<ev_mod_a[imod][ipart][ijsum][idt]<<" "<<endl;
	    /// [end] Hits hists ------------------------------------------------------------------------------
	    f1->Close();
	    f2->Close();
	    }// end loop over modules (diff plots)
	  }//end
	}//end //loop over dt values 
      }//end loop over diff sum values
    
      double av_evMod_b[10][nuSd+1][nDt];
      double av_evMod_a[10][nuSd+1][nDt];
      double err_av_evMod_b[10][nuSd+1][nDt];
      double err_av_evMod_a[10][nuSd+1][nDt];
      double av_resxMod_b[10][nuSd+1][nDt];
      double av_resxMod_a[10][nuSd+1][nDt];
      double err_av_resxMod_b[10][nuSd+1][nDt];
      double err_av_resxMod_a[10][nuSd+1][nDt];
      double av_resyMod_b[10][nuSd+1][nDt];
      double av_resyMod_a[10][nuSd+1][nDt];
      double err_av_resyMod_b[10][nuSd+1][nDt];
      double err_av_resyMod_a[10][nuSd+1][nDt];
      //    for(int isum=1;isum<nSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
      //      for(int isum=1;isum<nSd+1;isum=isum+10){ //loop over diff stat.samples (plot with diff markers)
    for(int ijsum=1;ijsum<nuSd+1;ijsum++){ //loop over diff stat.samples (plot with diff markers)
      //	      int isum = nSdstep*ijsum;
      int isum = 1+nSdstep*(ijsum-1);
	for(int idt=0;idt<nDt;idt++){//loop over dt values (x-axis of plots)
	  for(int imod=0;imod<nMod;imod++){// loop over modules (diff plots)
	    av_evMod_b[imod][ijsum][idt]=0;
	    av_evMod_a[imod][ijsum][idt]=0;
	    err_av_evMod_b[imod][ijsum][idt]=0;
	    err_av_evMod_a[imod][ijsum][idt]=0;

	    av_resxMod_b[imod][ijsum][idt]=0;
	    av_resxMod_a[imod][ijsum][idt]=0;
	    err_av_resxMod_b[imod][ijsum][idt]=0;
	    err_av_resxMod_a[imod][ijsum][idt]=0;
	    av_resyMod_b[imod][ijsum][idt]=0;
	    av_resyMod_a[imod][ijsum][idt]=0;
	    err_av_resyMod_b[imod][ijsum][idt]=0;
	    err_av_resyMod_a[imod][ijsum][idt]=0;
	  }//end loop over modules (diff plots)
	}//end loop over dt
      }//end loop over diif sum values

      double i_nS = 1./nSd;
      double i_nS_r = 1./(nSd-1);
      /// average calculatuon -------------------------------------------
      //      for(int isum=1;isum<nSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
      //      for(int isum=1;isum<nSd+1;isum=isum+10){ //loop over diff stat.samples (plot with diff markers)
    for(int ijsum=1;ijsum<nuSd+1;ijsum++){ //loop over diff stat.samples (plot with diff markers)
      //	      int isum = nSdstep*ijsum;
      int isum = 1+nSdstep*(ijsum-1);
	for(int idt=0;idt<nDt;idt++){//loop over dt values (x-axis of plots)
	  for(int imod=0;imod<nMod;imod++){// loop over modules (diff plots)
	    for(int ipart=0;ipart<isum;ipart++){//loop over samples
	      av_evMod_b[imod][ijsum][idt]+= ev_mod_b[imod][ipart][ijsum][idt];
	      av_evMod_a[imod][ijsum][idt]+= ev_mod_a[imod][ipart][ijsum][idt];

	      av_resxMod_b[imod][ijsum][idt]+= resx_mod_b[imod][ipart][ijsum][idt];
	      av_resxMod_a[imod][ijsum][idt]+= resx_mod_a[imod][ipart][ijsum][idt];
	      av_resyMod_b[imod][ijsum][idt]+= resy_mod_b[imod][ipart][ijsum][idt];
	      av_resyMod_a[imod][ijsum][idt]+= resy_mod_a[imod][ipart][ijsum][idt];
	      //     cout<<" "<<ev_mod_b[imod][ipart][ijsum][idt]<<" "<<ev_mod_a[imod][ipart][ijsum][idt]<<" "<<endl;
	    }//end over samples
	    av_resxMod_b[imod][ijsum][idt] *= 1./isum;
	    av_resxMod_a[imod][ijsum][idt] *= 1./isum; 
	    av_resyMod_b[imod][ijsum][idt] *= 1./isum;
	    av_resyMod_a[imod][ijsum][idt] *= 1./isum; 
	    //    cout<<" av for sum of "<<isum<<" :"<<av_evMod_b[imod][ijsum][idt]<<" "<<av_evMod_a[imod][ijsum][idt]<<" "<<endl;
	  }//end loop over modules (diff plots)
	}//end loop over dt
      }//end loop over diif sum values 
      /// END average calculatuon --------------------------------------
      
      /// average err. calculatuon -------------------------------------------
      //      for(int isum=1;isum<nSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
    //      for(int isum=1;isum<nSd+1;isum=isum+10){ //loop over diff stat.samples (plot with diff markers)
    for(int ijsum=1;ijsum<nuSd+1;ijsum++){ //loop over diff stat.samples (plot with diff markers)
      //	      int isum = nSdstep*ijsum;
      int isum = 1+nSdstep*(ijsum-1);
      	for(int idt=0;idt<nDt;idt++){//loop over dt values (x-axis of plots)
      	  for(int imod=0;imod<nMod;imod++){// loop over modules (diff plots)
      	    for(int ipart=0;ipart<isum;ipart++){//loop over samples
      	      err_av_resxMod_b[imod][ijsum][idt] += TMath::Power((resx_mod_b[imod][ipart][ijsum][idt]-av_resxMod_b[imod][ijsum][idt]),2);
      	      err_av_resxMod_a[imod][ijsum][idt] += TMath::Power((resx_mod_a[imod][ipart][ijsum][idt]-av_resxMod_a[imod][ijsum][idt]),2);
	      err_av_resyMod_b[imod][ijsum][idt] += TMath::Power((resy_mod_b[imod][ipart][ijsum][idt]-av_resyMod_b[imod][ijsum][idt]),2);
      	      err_av_resyMod_a[imod][ijsum][idt] += TMath::Power((resy_mod_a[imod][ipart][ijsum][idt]-av_resyMod_a[imod][ijsum][idt]),2);
      	    }//end over samples
      	    err_av_resyMod_b[imod][ijsum][idt] *=1./isum;
      	    err_av_resyMod_b[imod][ijsum][idt] = sqrt(err_av_resyMod_b[imod][ijsum][idt]);
      	    err_av_resyMod_a[imod][ijsum][idt] *=1./isum;
      	    err_av_resyMod_a[imod][ijsum][idt] = sqrt(err_av_resyMod_a[imod][ijsum][idt]);
      	    // cout<<" ipart="<<ipart<<": err_av_resxMod_a["<<imod<<"]["<<isum<<"]["<<idt<<"] = "<<err_av_resxMod_a[imod][ijsum][idt]<<endl;
      	    // cout<<" ipart="<<ipart<<": err_av_resyMod_a["<<imod<<"]["<<isum<<"]["<<idt<<"] = "<<err_av_resyMod_a[imod][ijsum][idt]<<endl;
      	  }//end loop over modules (diff plots)
      	}//end loop over dt
      }//end loop over diif sum values 
      /// END average err. calculatuon -------------------------------------


      // // PRINT ---------------------------
      // for(int isum=1;isum<nSd+1;isum++){//loop over diff stat.samples (plot with diff markers)
      // 	for(int idt=0;idt<nDt;idt++){//loop over dt values (x-axis of plots)
      // 	  for(int imod=0;imod<10;imod++){// loop over modules (diff plots)
      // 	    cout<<"Module #"<<imod<<"@ isum="<<isum<<", idt="<<idt<<" has "<<av_evMod_b[imod][ijsum][idt]<<"+/-"<<err_av_evMod_b[imod][ijsum][idt]<<endl;
      // 	    }
      // 	  }
      // 	}
      // // END PRINT ---------------------

      TString out = resname+".root";
      TFile *f = new TFile(out,"RECREATE");

      /// Let's draw now! ---------------------------------------------------------------------
      TCanvas cRes;  
      for(int imod=0;imod<nMod;imod++){// loop over modules (diff plots)
	TGraphErrors *gr_av_b[nuSd+1];
	TGraphErrors *gr_av_a[nuSd+1];
	TMultiGraph *mgr_av_b = new TMultiGraph();
	TMultiGraph *mgr_av_a = new TMultiGraph();
	TGraphErrors *gr_avResx_b[nuSd+1];
	TGraphErrors *gr_avResx_a[nuSd+1];
	TMultiGraph *mgr_avResx_b = new TMultiGraph();
	TMultiGraph *mgr_avResx_a = new TMultiGraph();
	TGraphErrors *gr_avResy_b[nuSd+1];
	TGraphErrors *gr_avResy_a[nuSd+1];
	TMultiGraph *mgr_avResy_b = new TMultiGraph();
	TMultiGraph *mgr_avResy_a = new TMultiGraph();
	TLegend *leg_b = new TLegend(0.91,0.65,0.99,0.98);
	leg_b->SetHeader("BEFORE");
	TLegend *leg_a = new TLegend(0.91,0.65,0.99,0.98);
	leg_a->SetHeader("AFTER");
	//	for(int isum=nSd;isum>0;isum--){//loop over diff stat.samples (plot with diff markers)
	//	for(int isum=nSd;isum>0;isum=isum-10){//loop over diff stat.samples (plot with diff markers)
	for(int ijsum=nuSd;ijsum>0;ijsum--){ //loop over diff stat.samples (plot with diff markers)
	  //	  int isum = nSdstep*ijsum;
	  int isum = 1+nSdstep*(ijsum-1);
	  TString sname = " ";
	  sname +=" sum of ";
	  sname += isum;
	  gr_av_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_evMod_b[imod][ijsum],0,err_av_evMod_b[imod][ijsum]);
	  gr_av_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_evMod_a[imod][ijsum],0,err_av_evMod_a[imod][ijsum]);
	  gr_avResx_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_b[imod][ijsum],0,err_av_resxMod_b[imod][ijsum]);
	  gr_avResx_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_a[imod][ijsum],0,err_av_resxMod_a[imod][ijsum]);
	  gr_avResx_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_b[imod][ijsum],0,err_av_resxMod_b[imod][ijsum]);
	  gr_avResx_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resxMod_a[imod][ijsum],0,err_av_resxMod_a[imod][ijsum]);
	  gr_avResy_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_b[imod][ijsum],0,err_av_resyMod_b[imod][ijsum]);
	  gr_avResy_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_a[imod][ijsum],0,err_av_resyMod_a[imod][ijsum]);
	  gr_avResy_b[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_b[imod][ijsum],0,err_av_resyMod_b[imod][ijsum]);
	  gr_avResy_a[ijsum] = new TGraphErrors(nDt,v_tr_sc, av_resyMod_a[imod][ijsum],0,err_av_resyMod_a[imod][ijsum]);
	  TString name_ev = "Ev_";
	  name_ev +=imod;
	  name_ev +="mod_sum";
	  name_ev += isum;
	  gr_av_b[ijsum]->SetName(name_ev+"_b");
	  gr_av_b[ijsum]->SetMarkerStyle(20+ijsum); 
	  gr_av_b[ijsum]->SetMarkerColor(kGreen-3); 
	  gr_av_b[ijsum]->SetMarkerSize(1.5); 
	  gr_av_a[ijsum]->SetName(name_ev+"_a");
	  gr_av_a[ijsum]->SetMarkerStyle(20+ijsum); 
	  gr_av_a[ijsum]->SetMarkerColor(kOrange+7); 
	  gr_av_a[ijsum]->SetMarkerSize(1.5); 
	  mgr_av_b->Add(gr_av_b[ijsum]);
	  mgr_av_a->Add(gr_av_a[ijsum]);
	  gr_av_b[ijsum]->Write();
	  gr_av_a[ijsum]->Write();
	  leg_b->AddEntry(gr_av_b[ijsum],sname,"P");
	  leg_a->AddEntry(gr_av_a[ijsum],sname,"P");

	  gr_avResx_b[ijsum]->SetName(name_ev+"_b");
	  gr_avResx_b[ijsum]->SetMarkerStyle(20+ijsum); 
	  gr_avResx_b[ijsum]->SetMarkerColor(kGreen-3); 
	  gr_avResx_b[ijsum]->SetMarkerSize(1.5); 
	  gr_avResx_a[ijsum]->SetName(name_ev+"_a");
	  gr_avResx_a[ijsum]->SetMarkerStyle(20+ijsum); 
	  gr_avResx_a[ijsum]->SetMarkerColor(kOrange+7); 
	  gr_avResx_a[ijsum]->SetMarkerSize(1.5); 
	  mgr_avResx_b->Add(gr_avResx_b[ijsum]);
	  mgr_avResx_a->Add(gr_avResx_a[ijsum]);
	  gr_avResx_b[ijsum]->Write();
	  gr_avResx_a[ijsum]->Write();
	  gr_avResy_b[ijsum]->SetName(name_ev+"_b");
	  gr_avResy_b[ijsum]->SetMarkerStyle(20+ijsum); 
	  gr_avResy_b[ijsum]->SetMarkerColor(kGreen-3); 
	  gr_avResy_b[ijsum]->SetMarkerSize(1.5); 
	  gr_avResy_a[ijsum]->SetName(name_ev+"_a");
	  gr_avResy_a[ijsum]->SetMarkerStyle(20+ijsum); 
	  gr_avResy_a[ijsum]->SetMarkerColor(kOrange+7); 
	  gr_avResy_a[ijsum]->SetMarkerSize(1.5); 
	  mgr_avResy_b->Add(gr_avResy_b[ijsum]);
	  mgr_avResy_a->Add(gr_avResy_a[ijsum]);
	  gr_avResy_b[ijsum]->Write();
	  gr_avResy_a[ijsum]->Write();

	}//end over diff stat.samples (plot with diff markers)
  
      TString smname = "Module_";
      smname += imod;
      mgr_av_b->SetTitle(smname+"_b");
      mgr_av_a->SetTitle(smname+"_a");
      leg_a->SetFillColor(0);
      leg_b->SetFillColor(0);
      mgr_av_b->SetName(smname+"_b");
      mgr_av_a->SetName(smname+"_a");
      TString smnamex,smnamey;
      smnamex =smname;
      smnamex+="_resX";
      smnamey =smname;
      smnamey+="_resY";
      mgr_avResx_b->SetName(smnamex+"_b");
      mgr_avResx_a->SetName(smnamex+"_a");
      mgr_avResy_b->SetName(smnamey+"_b");
      mgr_avResy_a->SetName(smnamey+"_a");
      mgr_av_b->Draw("AP");
      mgr_av_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
      mgr_av_b->GetYaxis()->SetTitle("events, 10^{3}");
      leg_b->Draw();
      cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
      cRes.Clear();
      mgr_av_a->Draw("AP");
      mgr_av_a->GetXaxis()->SetTitle("#delta_{t}, #mum");
      mgr_av_a->GetYaxis()->SetTitle("events, 10^{3}");
      leg_a->Draw();
      cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
      cRes.Clear();
      mgr_av_a->Write();
      mgr_av_b->Write();
      if(imod<1){
	leg_b->SetName("LegBefore");
	leg_a->SetName("LegAfter");
	leg_b->Write();
	leg_a->Write();
      }

      mgr_avResx_b->Draw("AP");
      mgr_avResx_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
      mgr_avResx_b->GetYaxis()->SetTitle("res_{X}, mean [#mum]");
      leg_b->Draw();
      cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
      cRes.Clear();
      mgr_avResx_a->Draw("AP");
      mgr_avResx_a->GetXaxis()->SetTitle("#delta_{t}, #mum");
      mgr_avResx_a->GetYaxis()->SetTitle("res_{X}, mean [#mum]");
      leg_a->Draw();
      cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
      cRes.Clear();
      mgr_avResx_a->Write();
      mgr_avResx_b->Write();

      mgr_avResy_b->Draw("AP");
      mgr_avResy_b->GetXaxis()->SetTitle("#delta_{t}, #mum");
      mgr_avResy_b->GetYaxis()->SetTitle("res_{Y}, mean [#mum]");
      leg_b->Draw();
      cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
      cRes.Clear();
      mgr_avResy_a->Draw("AP");
      mgr_avResy_a->GetXaxis()->SetTitle("#delta_{t}, #mum");
      mgr_avResy_a->GetYaxis()->SetTitle("res_{Y}, mean [#mum]");
      leg_a->Draw();
      cRes.Print(resname_pdf_o); //write canvas and keep the ps file open
      cRes.Clear();
      mgr_avResy_a->Write();
      mgr_avResy_b->Write();
     
      }//end loop over modules (diff plots)
      cRes.Print(resname_pdf_c); //write canvas and close the ps file 
      ///--------------------------------------------------------------------------------------------
      f->Write();
      f->Close();
    }
  }
  return 0;
}