/// g++ $(root-config --cflags --glibs) -L /home/karavdin/pandaroot13//buildPanda/lib AllSamplesKnossosSummaryAlign.C -o /home/karavdin/pandaroot13/macro/lmd/LMD_alignment/AllSamplesKnossosSummaryAlign
// ./AllSamplesKnossosSummaryAlign
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <TString.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TGraphErrors.h>
#include <TMultiGraph.h>
#include <sstream>
#include "TStyle.h"

using namespace std;
//void TOTALSummaryAlign(TString pathG="/panda/pandaroot/macro/lmd/tmpOutputAlignTMP/", double tr_sc=200, double rt_sc=9)
int main(int __argc,char *__argv[]){
  // TStyle* m_gStyle = new TStyle();
  gStyle->SetOptStat(0);
  gStyle->SetLabelSize(0.045,"xyz");
  gStyle->SetTitleYSize(0.04);
  gStyle->SetTitleXSize(0.04);

  double tr_sc=300; double rt_sc=3; TString strrt_sc="0";
  TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/";
  // const int misAlinSamp=0;
  const int nAS=10; // number of misalignment consts samples
  //   const int nAS=5; // number of misalignment consts samples
  const  double i_nAS = 1./nAS;
  const int nS=100; // number of data samples for each misalignment consts sample
  //  const int nS=20;
  const double i_nS = 1./nS;
  const double i_nS_r = 1./(nS-1);
  const double i_nAS_r = 1./(nAS-1);
  TString pathGb=pathG+"/RECdataAUG2014/BOX/mom_15/";
  const int ntrksSample = 5e4;
  const double i_TrksSimi = 100./ntrksSample; //relative to simulated in %
  const int nParDt=11;
  double v_Dt[nParDt]={0, 10, 50, 100,200,300,400,500,600,700,800};
// const int nParDt=7;
//  double v_Dt[nParDt]={0, 10, 50, 100,200,300,400,};
  // const int nParDt=10;
  // double v_Dt[nParDt]={0, 10, 50, 100,200,300,400,500,600,800};
  // const int nParDt=9;
  // double v_Dt[nParDt]={0, 50, 100,200,300,400,500,600,800};
  // const int nParDt=1;
  // double v_Dt[nParDt]={700};
  // const int nParDa=5;
  // double v_Da[nParDa]={0,0.2,1,3,5};
  // TString strv_Da[nParDa]={"0","0_2","1","3","5"};
  const int nParDa=4;
  double v_Da[nParDa]={0,1,3,5};
  TString strv_Da[nParDa]={"0","1","3","5"};
  // const int nParDa=2;
  // double v_Da[nParDa]={1,3};
  // TString strv_Da[nParDa]={"1","3"};
  for(int iDt=0;iDt<nParDt;iDt++){
    tr_sc = v_Dt[iDt];
    for(int iDa=0;iDa<nParDa;iDa++){
      rt_sc= v_Da[iDa];
      strrt_sc= strv_Da[iDa];
      cout<<"work with dt~"<<tr_sc<<" da~"<<rt_sc<<endl;
      double misconst_av_TrksBef[nAS+1][nS+1];double misconst_av_TrksAft[nAS+1][nS+1];
      double misconst_avEr_TrksBef[nAS+1][nS+1]; double misconst_avEr_TrksAft[nAS+1][nS+1];
      double misconst_av_mis_a_mean[nAS+1][6][nS+1];
      double misconst_av_mis_b_mean[nAS+1][6][nS+1];
      double misconst_avEr_mis_a_mean[nAS+1][6][nS+1];
      double misconst_avEr_mis_b_mean[nAS+1][6][nS+1];
      //How to save data
      TString resname= pathGb+"/ResultsSamplesSummary_Av_";
      resname += nAS;
      resname += "AlignSamples_";
      resname += nS;
      resname += "DataSamples_";
      resname+=tr_sc;
      resname+="um_";
      resname+=strrt_sc;
      resname+="mrad";
      
      TString resname_pdf = resname+".pdf";
      TString resname_pdf_o = resname_pdf+"(";
      TString resname_pdf_c = resname_pdf+")";

      for(int misAlinSamp=0;misAlinSamp<nAS;misAlinSamp++){
	double TrksBef[nS+1];
	double TrksAft[nS+1];
	double mis_b[6][nS+1];
	double mis_a[6][nS+1];
	cout<<"misAlinSamp "<<misAlinSamp<<endl;
	for(int i=0;i<nS;i++){
	  //	  if (i%20 == 0)  
	  //	  cout<<"Sum of "<<i<<"sample(s)"<<endl;
	  //	  for(int ipart=0;ipart<nS;ipart++){
	    //	    cout<<"part #"<<ipart<<endl;
	    ///Misalignment constants -----------------------------------------------------------------
	    TString namemisc = pathGb+"/";
	    namemisc+=tr_sc;
	    namemisc+="mkm_";
	    namemisc+=strrt_sc;
	    namemisc+="mrad/ConstSample";
	    namemisc+=misAlinSamp;
	    namemisc+="/RECsum";
	    namemisc +=(i+1);
	    namemisc+="DataSamples/KnossosResults.root";
	    TFile *fmisc = new TFile(namemisc,"READ");
	    if (fmisc->IsZombie()) {
	      cout << "!!! Error opening file " <<namemisc<< endl;
	      continue;
	      // return;
	    }
	    TH2F *hStatKnossos = (TH2F*)fmisc->Get("haStat_out");
	    TH2F *hStatKnossos0 = (TH2F*)fmisc->Get("haStat_in");
	    TrksBef[i] = hStatKnossos0->ProjectionY()->GetMean();
	    TrksAft[i] = hStatKnossos->ProjectionY()->GetMean();
	    TH2F *hmis_b_0 = (TH2F*)fmisc->Get("mis_before_0");
	    mis_b[0][i] = 1e4*(hmis_b_0->ProjectionY()->GetRMS());
	    TH2F *hmis_b_1 = (TH2F*)fmisc->Get("mis_before_1");
	    mis_b[1][i] = 1e4*(hmis_b_1->ProjectionY()->GetRMS());
	    TH2F *hmis_b_2 = (TH2F*)fmisc->Get("mis_before_2");
	    mis_b[2][i] = 1e4*(hmis_b_2->ProjectionY()->GetRMS());
	    TH2F *hmis_b_3 = (TH2F*)fmisc->Get("mis_before_3");
	    mis_b[3][i] = 1e3*(hmis_b_3->ProjectionY()->GetRMS());
	    TH2F *hmis_b_4 = (TH2F*)fmisc->Get("mis_before_4");
	    mis_b[4][i] = 1e3*(hmis_b_4->ProjectionY()->GetRMS());
	    TH2F *hmis_b_5 = (TH2F*)fmisc->Get("mis_before_5");
	    mis_b[5][i] = 1e3*(hmis_b_5->ProjectionY()->GetRMS());
    
	    TH2F *hmis_a_0 = (TH2F*)fmisc->Get("mis_diff_0");
	    mis_a[0][i] = 1e4*(hmis_a_0->ProjectionY()->GetRMS());
	    TH2F *hmis_a_1 = (TH2F*)fmisc->Get("mis_diff_1");
	    mis_a[1][i] = 1e4*(hmis_a_1->ProjectionY()->GetRMS());
	    TH2F *hmis_a_2 = (TH2F*)fmisc->Get("mis_diff_2");
	    mis_a[2][i] = 1e4*(hmis_a_2->ProjectionY()->GetRMS());
	    TH2F *hmis_a_3 = (TH2F*)fmisc->Get("mis_diff_3");
	    mis_a[3][i] = 1e3*(hmis_a_3->ProjectionY()->GetRMS());
	    TH2F *hmis_a_4 = (TH2F*)fmisc->Get("mis_diff_4");
	    mis_a[4][i] = 1e3*(hmis_a_4->ProjectionY()->GetRMS());
	    TH2F *hmis_a_5 = (TH2F*)fmisc->Get("mis_diff_5");
	    mis_a[5][i] = 1e3*(hmis_a_5->ProjectionY()->GetRMS());
	    //	    cout<<"mis_a[0]["<<i<<"] ="<<mis_a[0][i]<<" mis_a[1]["<<i<<"] ="<<mis_a[1][i]<<" mis_a[5]["<<i<<"] ="<<mis_a[5][i]<<endl;
	    /// END (Misalignment constants) -----------------------------------------------------------
	    fmisc->Close();
	  }
	//	}
	//return;
	double av_TrksBef[nS+1];double av_TrksAft[nS+1];
	double avEr_TrksBef[nS+1]; double avEr_TrksAft[nS+1];
	double av_mis_a_mean[6][nS+1];
	double av_mis_b_mean[6][nS+1];
	double avEr_mis_a_mean[6][nS+1];
	double avEr_mis_b_mean[6][nS+1];

	for(int i=0;i<nS;i++){
	  av_TrksBef[i]=TrksBef[i];
	  av_TrksAft[i]=TrksAft[i];
	  avEr_TrksBef[i] = 0;
	  avEr_TrksAft[i] = 0;
	  for(int ialc=0;ialc<6;ialc++){
	    av_mis_a_mean[ialc][i] = mis_a[ialc][i];
	    av_mis_b_mean[ialc][i] = mis_b[ialc][i];
	    avEr_mis_a_mean[ialc][i] = 0;
	    avEr_mis_b_mean[ialc][i] = 0;
	  }
	}
  
	for(int i=0;i<nS;i++){
	  misconst_av_TrksBef[misAlinSamp][i] = av_TrksBef[i];
	  misconst_av_TrksAft[misAlinSamp][i] = av_TrksAft[i];
	  misconst_avEr_TrksBef[misAlinSamp][i] = avEr_TrksBef[i];
	  misconst_avEr_TrksAft[misAlinSamp][i] = avEr_TrksAft[i];
	  for(int ialc=0;ialc<6;ialc++){
	    misconst_av_mis_a_mean[misAlinSamp][ialc][i] = av_mis_a_mean[ialc][i];
	    misconst_av_mis_b_mean[misAlinSamp][ialc][i] = av_mis_b_mean[ialc][i];
	    misconst_avEr_mis_a_mean[misAlinSamp][ialc][i] = avEr_mis_a_mean[ialc][i];
	    misconst_avEr_mis_b_mean[misAlinSamp][ialc][i] = avEr_mis_b_mean[ialc][i];
	    //	    cout<<"["<<misAlinSamp<<" "<<ialc<<" "<<i<<"] "<<misconst_av_mis_a_mean[misAlinSamp][ialc][i]<<" "<<misconst_avEr_mis_a_mean[misAlinSamp][ialc][i]<<endl;
	  }
	}
      }//loop over diff. misalignment samples

      /// Make avarage for different data sums
      double av_misconst_av_TrksBef[nS+1];double av_misconst_av_TrksAft[nS+1];
      double avEr_misconst_avEr_TrksBef[nS+1]; double avEr_misconst_avEr_TrksAft[nS+1]; 
      double av_misconst_av_mis_a_mean[6][nS+1];
      double av_misconst_av_mis_b_mean[6][nS+1];
      double avEr_misconst_avEr_mis_a_mean[6][nS+1];
      double avEr_misconst_avEr_mis_b_mean[6][nS+1];

      for(int i=0;i<nS;i++){
	av_misconst_av_TrksBef[i] = 0;
	av_misconst_av_TrksAft[i] = 0;
	for(int ialc=0;ialc<6;ialc++){
	  av_misconst_av_mis_a_mean[ialc][i] = 0;
	  av_misconst_av_mis_b_mean[ialc][i] = 0;
	  avEr_misconst_avEr_mis_a_mean[ialc][i] = 0;
	  avEr_misconst_avEr_mis_b_mean[ialc][i] = 0;
	}
      }

      for(int i=0;i<nS;i++){
	for(int misAlinSamp=0;misAlinSamp<nAS;misAlinSamp++){
	  av_misconst_av_TrksBef[i] += misconst_av_TrksBef[misAlinSamp][i];
	  av_misconst_av_TrksAft[i] += misconst_av_TrksAft[misAlinSamp][i];
	  for(int ialc=0;ialc<6;ialc++){
	    av_misconst_av_mis_a_mean[ialc][i] += misconst_av_mis_a_mean[misAlinSamp][ialc][i];
	    av_misconst_av_mis_b_mean[ialc][i] += misconst_av_mis_b_mean[misAlinSamp][ialc][i];
	  }
	}//loop over diff. misalignment samples
	for(int ialc=0;ialc<6;ialc++){
	  av_misconst_av_mis_a_mean[ialc][i] *= i_nAS;
	  av_misconst_av_mis_b_mean[ialc][i] *= i_nAS;
	}
	av_misconst_av_TrksBef[i] *= i_nAS;
	av_misconst_av_TrksAft[i] *= i_nAS;
	//	cout<<"av_misconst_av_TrksBef["<<i<<"] = "<<av_misconst_av_TrksBef[i]<<endl;
      }

      //calculate RMS (+ error from previous avaraging)
      for(int i=0;i<nS;i++){
	for(int misAlinSamp=0;misAlinSamp<nAS;misAlinSamp++){
	  avEr_misconst_avEr_TrksBef[i] = TMath::Power((av_misconst_av_TrksBef[i]-misconst_av_TrksBef[misAlinSamp][i]),2) 
	    + TMath::Power(misconst_avEr_TrksBef[misAlinSamp][i],2);
	  avEr_misconst_avEr_TrksAft[i] = TMath::Power((av_misconst_av_TrksAft[i]-misconst_av_TrksAft[misAlinSamp][i]),2) 
	    + TMath::Power(misconst_avEr_TrksAft[misAlinSamp][i],2);
	  for(int ialc=0;ialc<6;ialc++){
	    avEr_misconst_avEr_mis_a_mean[ialc][i] += 
	      TMath::Power((av_misconst_av_mis_a_mean[ialc][i]-misconst_av_mis_a_mean[misAlinSamp][ialc][i]),2);
	    //	    if(i==49) cout<<av_misconst_av_mis_a_mean[ialc][i]<<" "<<misconst_av_mis_a_mean[misAlinSamp][ialc][i]<<" "<< avEr_misconst_avEr_mis_a_mean[ialc][i]<<endl;
	       //	       +TMath::Power(misconst_avEr_mis_a_mean[misAlinSamp][ialc][i],2));
	    avEr_misconst_avEr_mis_b_mean[ialc][i] += 
	       TMath::Power((av_misconst_av_mis_b_mean[ialc][i]-misconst_av_mis_b_mean[misAlinSamp][ialc][i]),2);
		//	       +TMath::Power(misconst_avEr_mis_b_mean[misAlinSamp][ialc][i],2));
	  }
	}
	avEr_misconst_avEr_TrksBef[i] *=i_nAS_r;
	avEr_misconst_avEr_TrksBef[i] = sqrt(avEr_misconst_avEr_TrksBef[i]);
	avEr_misconst_avEr_TrksAft[i] *=i_nAS_r;
	avEr_misconst_avEr_TrksAft[i] = sqrt(avEr_misconst_avEr_TrksAft[i]);
	for(int ialc=0;ialc<6;ialc++){
	  avEr_misconst_avEr_mis_a_mean[ialc][i] *=i_nAS_r;
	  avEr_misconst_avEr_mis_a_mean[ialc][i] = sqrt(avEr_misconst_avEr_mis_a_mean[ialc][i]);
	  avEr_misconst_avEr_mis_b_mean[ialc][i] *=i_nAS_r;
	  avEr_misconst_avEr_mis_b_mean[ialc][i] = sqrt(avEr_misconst_avEr_mis_b_mean[ialc][i]);
	}
      }

      /// END Make avarage for different misalignment sensors

      // ///////// Draw and save results
      TCanvas c1;  
      //put together missalignment consts results
      TGraphErrors *gr_mis_b[6];
      TGraphErrors *gr_mis_a[6];
      TMultiGraph *mgr_mis[6];
      for(int ialc=0;ialc<6;ialc++){
        mgr_mis[ialc] = new TMultiGraph();
        gr_mis_b[ialc] = new TGraphErrors(nS,av_misconst_av_TrksBef,av_misconst_av_mis_b_mean[ialc],avEr_misconst_avEr_TrksBef,avEr_misconst_avEr_mis_b_mean[ialc]);
        gr_mis_b[ialc]->SetMarkerStyle(20);
        gr_mis_b[ialc]->SetMarkerColor(kGreen-3);
        gr_mis_b[ialc]->SetMarkerSize(2.5);
	//	for(int kl=0;kl<nS+1;kl++)
	  //	  cout<<"["<<ialc<<" "<<kl<<"] "<<av_misconst_av_mis_a_mean[ialc][kl]<<" "<<avEr_misconst_avEr_mis_a_mean[ialc][kl]<<endl;
	
        gr_mis_a[ialc] = new TGraphErrors(nS,av_misconst_av_TrksBef,av_misconst_av_mis_a_mean[ialc],avEr_misconst_avEr_TrksBef,avEr_misconst_avEr_mis_a_mean[ialc]);
        gr_mis_a[ialc]->SetMarkerStyle(29);
        gr_mis_a[ialc]->SetMarkerColor(kOrange+7);
        gr_mis_a[ialc]->SetMarkerSize(2.5);
	//  mgr_mis[ialc]->Add(gr_mis_b[ialc]);
        mgr_mis[ialc]->Add(gr_mis_a[ialc]);
        mgr_mis[ialc]->Draw("AP");
        mgr_mis[ialc]->GetXaxis()->SetTitle("Average number of rec. trks");
	if(ialc<3){
	  mgr_mis[ialc]->GetHistogram()->SetMaximum(100.);
	  mgr_mis[ialc]->GetHistogram()->SetMinimum(-10.);
	}
	else{
	mgr_mis[ialc]->GetHistogram()->SetMaximum(5.1);
	mgr_mis[ialc]->GetHistogram()->SetMinimum(-1.1);
	}

        if(ialc==0){ 
	  mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{x}, #mum");
	  gr_mis_a[ialc]->SetName("dX_a");
	  gr_mis_b[ialc]->SetName("dX_b");
	}
        if(ialc==1){ 
	  mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{y}, #mum");
	  gr_mis_a[ialc]->SetName("dY_a");
	  gr_mis_b[ialc]->SetName("dY_b");
	}
        if(ialc==2){
	  mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{z}, #mum");
	  gr_mis_a[ialc]->SetName("dZ_a");
	  gr_mis_b[ialc]->SetName("dZ_b");
	}
        if(ialc==3){
	  mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#alpha}, mrad");
	  gr_mis_a[ialc]->SetName("dA_a");
	  gr_mis_b[ialc]->SetName("dA_b");
	}
        if(ialc==4){
	  mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#beta}, mrad");
	  gr_mis_a[ialc]->SetName("dB_a");
	  gr_mis_b[ialc]->SetName("dB_b");
	}
        if(ialc==5){ 
	  mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#gamma}, mrad");
	  gr_mis_a[ialc]->SetName("dG_a");
	  gr_mis_b[ialc]->SetName("dG_b");
	}
        if(ialc<2 || ialc==5) c1.Print(resname_pdf_o); //write canvas and keep the file open
         c1.Clear();
      }
      // TMultiGraph *mgr_stat = new TMultiGraph();
      TGraphErrors *grstat = new TGraphErrors(nS,av_misconst_av_TrksBef,av_misconst_av_TrksAft,avEr_misconst_avEr_TrksBef,avEr_misconst_avEr_TrksAft);
      grstat->SetMarkerStyle(34);
      grstat->SetMarkerSize(2.5);
      grstat->GetXaxis()->SetTitle("Average number of rec. trks (1st iter)");
      grstat->GetYaxis()->SetTitle("Average number of rec. trks (last iter)");
      grstat->Fit("pol1");
      gStyle->SetOptFit();
      grstat->Draw("AP");
      //  TGraphErrors *grstatRef = new TGraphErrors(nS,av_misconst_av_TrksBef,av_TrksRef,avEr_misconst_avEr_TrksBef,avEr_TrksRef);
      //  grstatRef->SetMarkerStyle(24);
      //  grstatRef->SetMarkerSize(2.5);
      //  TGraphErrors *grstatAft = new TGraphErrors(nS,av_misconst_av_TrksBef,av_TrksAft,avEr_misconst_avEr_TrksBef,avEr_TrksAft);
      //  grstatAft->SetMarkerStyle(22);
      //  grstatAft->SetMarkerSize(2.5);
      //  grstatAft->SetTitle("Average number of rec. trks (empty = refer; filled = after Knossos)");
      //  // mgr_stat->Add(grstatBef);
      //  mgr_stat->Add(grstatRef);
      //  mgr_stat->Add(grstatAft);
      //  mgr_stat->Draw("AP");
      //  mgr_stat->GetXaxis()->SetTitle("Number of sim. trks for alignment");
      //  mgr_stat->GetYaxis()->SetTitle("N_{REC} [empty = N^{ before}_{REC}; filled = N^{ after}_{REC}]");
      //  //mgr_stat->SetTitle("Average number of rec. trks (empty = before Knossos; filled = after Knossos)");
      //      c1.Clear();
      c1.Print(resname_pdf_c); //write canvas and close ps file
    

       TString out = resname+".root";
       TFile *f = new TFile(out,"RECREATE");
       for(int ialc=0;ialc<6;ialc++){
         gr_mis_b[ialc]->Write();
         gr_mis_a[ialc]->Write();
         mgr_mis[ialc]->Write();
       }
      //  mgr_stat->Write();
       f->Write();
       f->Close();
       //  //////// END Draw and save results
    }
  }
  return 0;
}