/// g++ $(root-config --cflags --glibs) -L /home/karavdin/pandaroot13//buildPanda/lib AllSamplesSummaryAlign.C -o /home/karavdin/pandaroot13/macro/lmd/LMD_alignment/AllSamplesSummaryAlign

#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>

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[]){
  double tr_sc=300; double rt_sc=3;
  TString pathG="/home/karavdin/datastorage/AlignmentLMDpixel/";
  // const int misAlinSamp=0;
  const int nAS=2;
  const  double i_nAS = 1./nAS;
  const int nS=10;
  //  const int nS=3;
  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+"/RECdataJuly2014/BOX/mom_15/";
  const int ntrksSample = 2e4;
  const double i_TrksSimi = 100./ntrksSample; //relative to simulated in %
  // const int nParDt=10;
  // double v_Dt[nParDt]={0, 50,100,200,300,400,500,600,800,1000};
  const int nParDt=1;
  double v_Dt[nParDt]={100};
  const int nParDa=1;
  double v_Da[nParDa]={0};
  // const int nParDt=1;
  // double v_Dt[nParDt]={600};
  // const int nParDa=1;
  //  double v_Da[nParDa]={0};
  for(int iDt=0;iDt<nParDt;iDt++){
    tr_sc = v_Dt[iDt];
    for(int iDa=0;iDa<nParDa;iDa++){
      rt_sc= v_Da[iDa];
      cout<<"work with dt~"<<tr_sc<<" da~"<<rt_sc<<endl;
     
      double misconst_av_theta_mean_b[nAS][nS+1],misconst_av_theta_rms_b[nAS][nS+1];
      double misconst_av_theta_mean_a[nAS][nS+1],misconst_av_theta_rms_a[nAS][nS+1];
      double misconst_av_theta_mean_r[nAS][nS+1],misconst_av_theta_rms_r[nAS][nS+1];
      double misconst_av_TrksBef[nAS][nS+1];double misconst_av_TrksAft[nAS][nS+1];double misconst_av_TrksDiff[nAS][nS+1];
      double misconst_av_TrksRef[nAS][nS+1];
      double misconst_avEr_theta_mean_b[nAS][nS+1],misconst_avEr_theta_rms_b[nAS][nS+1];
      double misconst_avEr_theta_mean_a[nAS][nS+1],misconst_avEr_theta_rms_a[nAS][nS+1];
      double misconst_avEr_theta_mean_r[nAS][nS+1],misconst_avEr_theta_rms_r[nAS][nS+1];
      double misconst_avEr_TrksBef[nAS][nS+1]; double misconst_avEr_TrksAft[nAS][nS+1];double misconst_avEr_TrksDiff[nAS][nS+1]; double misconst_avEr_TrksRef[nAS][nS+1]; 
      double misconst_av_mis_a_mean[nAS][6][nS+1];
      double misconst_av_mis_b_mean[nAS][6][nS+1];
      double misconst_avEr_mis_a_mean[nAS][6][nS+1];
      double misconst_avEr_mis_b_mean[nAS][6][nS+1];
      //How to save data
      TString resname= pathGb+"/ResultsSamplesSummary_AvAlignSamples";
      // resname+=misAlinSamp;
      resname+="_";
      resname+=tr_sc;
      resname+="um_";
      resname+=rt_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++){
	//   const int nS=2;
	double TrksSim[nS+1];
	double TrksBef[nS+1][nS+1];
	double TrksAft[nS+1][nS+1];
	double TrksRef[nS+1][nS+1];
	double TrksDiff[nS+1][nS+1];

	double mis_b[nS+1][6][nS+1];
	double mis_a[nS+1][6][nS+1];

	double theta_mean_b[nS+1][nS+1], theta_mean_a[nS+1][nS+1], theta_mean_r[nS+1][nS+1];
	double theta_rms_b[nS+1][nS+1], theta_rms_a[nS+1][nS+1], theta_rms_r[nS+1][nS+1];
	TH1D *hthetabefore = new TH1D("hthetabefore","#theta_{MC}-#theta_{rec};#delta#theta,rad",1e2,-1e-3,1e-3);
	TH1D *hthetaafter = new TH1D("hthetaafter","#theta_{MC}-#theta_{rec};#delta#theta,rad",1e2,-1e-3,1e-3);
	TH1D *hthetaref = new TH1D("hthetaref","#theta_{MC}-#theta_{rec};#delta#theta,rad",1e2,-1e-3,1e-3);

	for(int ipart=0;ipart<nS;ipart++){
	  for(int i=0;i<nS;i++){
	    TrksSim[i]=(i+1)*ntrksSample ; //each sample contains ntrksSample events with 1trk/event
	    theta_mean_b[ipart][i] = 0;
	    theta_mean_a[ipart][i] = 0;
	    theta_mean_r[ipart][i] = 0;
	    theta_rms_b[ipart][i] = 0;
	    theta_rms_a[ipart][i] = 0;
	    theta_rms_r[ipart][i] = 0;
	  }
	}

	for(int i=0;i<nS;i++){
	  //  cout<<"Sum of "<<i<<"sample(s)"<<endl;
	  for(int ipart=0;ipart<nS;ipart++){
	    //  cout<<"part #"<<ipart<<endl;
	    int nTrks =  TrksSim[i];
	    //Where search for files
	    TString pathb = pathGb+"/";
	    pathb+=tr_sc;
	    pathb+="mkm_";
	    pathb+=rt_sc;
	    pathb+="mrad/ConstSample";
	    pathb+=misAlinSamp;
	    pathb+="/RECrawDataSample"; //ToDo: check for different misalignment consts
	    pathb +=ipart;

	    TString patha = pathGb+"/";
	    patha+=tr_sc;
	    patha+="mkm_";
	    patha+=rt_sc;
	    patha+="mrad/ConstSample";
	    patha+=misAlinSamp;
	    patha+="/RECsum"; //ToDo: check for different misalignment consts
	    patha +=(i+1);
	    patha+="DataSamples/RECalignDataSample";
	    patha+=ipart;

	    /// Theta resolution--------------------------------------------------------------------------------------  
	    TString name1 = pathb + "/Lumi_QA_before.root";
	    TString name2 = patha + "/Lumi_QA_after.root";
	    TString name3 = pathGb+"/0mkm_0mrad/ConstSample0/RECsum"; //ToDo: check for different misalignment consts
	    name3+=nS;
	    name3+="DataSamples/RECalignDataSample";
	    name3+=ipart;
	    name3+="/Lumi_QA_after.root";


	    TFile *f1 = new TFile(name1,"READ");
	    if (f1->IsZombie()) {
	      cout << "!!! Error opening file " <<name1<< endl;
	      continue;
	      //   return;
	    }
	    TH1 *hbefore = (TH1*)f1->Get("NearIP/hResTheta");
	    //cout<<"read QA before"<<endl;
	    //hbefore->Print();
	    hbefore->SetName("hResThetaBefore");
	    hbefore->SetTitle("#theta resolution");
	    //  hbefore->Print();
	    TrksBef[ipart][i] = hbefore->GetEntries();
	    theta_mean_b[ipart][i]  = 1e3*(hbefore->GetMean());
	    theta_rms_b[ipart][i]  = 1e3*(hbefore->GetRMS());
	    f1->Close();
	    // nsectors_b = (TNtuple*)f1->Get("nsectors");
	    TFile *f2 = new TFile(name2,"READ");
	    if (f2->IsZombie()) {
	      cout << "!!! Error opening file " <<name2<< endl;
	      continue;
	      // return;
	    }
	    //  cout<<" "<<name2<<endl;
	    TH1 *hafter = (TH1F*)f2->Get("NearIP/hResTheta");
	    //  cout<<"read QA after"<<endl;
	    hafter->SetName("hResThetaAfter");
	    hafter->SetTitle("#theta resolution");

	    TrksAft[ipart][i] = hafter->GetEntries();
	    TrksDiff[ipart][i] = (TrksAft[ipart][i] - TrksBef[ipart][i]);
	    theta_mean_a[ipart][i]  = 1e3*(hafter->GetMean());
	    theta_rms_a[ipart][i]  = 1e3*(hafter->GetRMS());
	    f2->Close();
  
	    TFile *f3 = new TFile(name3,"READ");
  
	    if (f3->IsZombie()) {
	      cout << "!!! Error opening file " <<name3<< endl;
	      continue;
	      // return;
	    }
	    TH1 *hrefer = (TH1F*)f3->Get("NearIP/hResTheta");
	    //  cout<<"read QA ref"<<endl;
	    hrefer->SetName("hResThetaRef");
	    hrefer->SetTitle("#theta resolution");
	    TrksRef[ipart][i] = hrefer->GetEntries();
	    theta_mean_r[ipart][i]  = 1e3*(hrefer->GetMean());
	    theta_rms_r[ipart][i]  = 1e3*(hrefer->GetRMS());
	    /// END (Theta resolution) ---------------------------------------------------------------

 
	    f3->Close();
	    ///Misalignment constants -----------------------------------------------------------------

	    TString namemisc = pathGb+"/";
	    namemisc+=tr_sc;
	    namemisc+="mkm_";
	    namemisc+=rt_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 *hmis_b_0 = (TH2F*)fmisc->Get("mis_before_0");
	    mis_b[ipart][0][i] = 1e4*(hmis_b_0->ProjectionY()->GetRMS());
	    TH2F *hmis_b_1 = (TH2F*)fmisc->Get("mis_before_1");
	    mis_b[ipart][1][i] = 1e4*(hmis_b_1->ProjectionY()->GetRMS());
	    TH2F *hmis_b_2 = (TH2F*)fmisc->Get("mis_before_2");
	    mis_b[ipart][2][i] = 1e4*(hmis_b_2->ProjectionY()->GetRMS());
	    TH2F *hmis_b_3 = (TH2F*)fmisc->Get("mis_before_3");
	    mis_b[ipart][3][i] = 1e3*(hmis_b_3->ProjectionY()->GetRMS());
	    TH2F *hmis_b_4 = (TH2F*)fmisc->Get("mis_before_4");
	    mis_b[ipart][4][i] = 1e3*(hmis_b_4->ProjectionY()->GetRMS());
	    TH2F *hmis_b_5 = (TH2F*)fmisc->Get("mis_before_5");
	    mis_b[ipart][5][i] = 1e3*(hmis_b_5->ProjectionY()->GetRMS());
    
	    TH2F *hmis_a_0 = (TH2F*)fmisc->Get("mis_diff_0");
	    mis_a[ipart][0][i] = 1e4*(hmis_a_0->ProjectionY()->GetRMS());
	    TH2F *hmis_a_1 = (TH2F*)fmisc->Get("mis_diff_1");
	    mis_a[ipart][1][i] = 1e4*(hmis_a_1->ProjectionY()->GetRMS());
	    TH2F *hmis_a_2 = (TH2F*)fmisc->Get("mis_diff_2");
	    mis_a[ipart][2][i] = 1e4*(hmis_a_2->ProjectionY()->GetRMS());
	    TH2F *hmis_a_3 = (TH2F*)fmisc->Get("mis_diff_3");
	    mis_a[ipart][3][i] = 1e3*(hmis_a_3->ProjectionY()->GetRMS());
	    TH2F *hmis_a_4 = (TH2F*)fmisc->Get("mis_diff_4");
	    mis_a[ipart][4][i] = 1e3*(hmis_a_4->ProjectionY()->GetRMS());
	    TH2F *hmis_a_5 = (TH2F*)fmisc->Get("mis_diff_5");
	    mis_a[ipart][5][i] = 1e3*(hmis_a_5->ProjectionY()->GetRMS());
	    //	    cout<<"mis_a["<<ipart<<"][5]["<<i<<"] ="<<mis_a[ipart][5][i]<<endl;
	    /// END (Misalignment constants) -----------------------------------------------------------
	    fmisc->Close();
	  }
	}
	//return;
	double av_theta_mean_b[nS+1],av_theta_rms_b[nS+1];
	double av_theta_mean_a[nS+1],av_theta_rms_a[nS+1];
	double av_theta_mean_r[nS+1],av_theta_rms_r[nS+1];
	double av_TrksBef[nS+1];double av_TrksAft[nS+1];double av_TrksDiff[nS+1];double av_TrksRef[nS+1];
	double avEr_theta_mean_b[nS+1],avEr_theta_rms_b[nS+1];
	double avEr_theta_mean_a[nS+1],avEr_theta_rms_a[nS+1];
	double avEr_theta_mean_r[nS+1],avEr_theta_rms_r[nS+1];
	double avEr_TrksBef[nS+1]; double avEr_TrksAft[nS+1];double avEr_TrksDiff[nS+1]; double avEr_TrksRef[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_theta_mean_b[i] = 0;
	  av_theta_rms_b[i] = 0;
	  av_theta_mean_a[i] = 0;
	  av_theta_rms_a[i] = 0;
	  av_theta_mean_r[i] = 0;
	  av_theta_rms_r[i] = 0;
	  av_TrksBef[i]=0;
	  av_TrksRef[i]=0;
	  av_TrksAft[i]=0;
	  av_TrksDiff[i]=0;
	  avEr_theta_mean_b[i] = 0;
	  avEr_theta_rms_b[i] =0;
	  avEr_theta_mean_a[i] = 0;
	  avEr_theta_rms_a[i] = 0;
	  avEr_theta_mean_r[i] = 0;
	  avEr_theta_rms_r[i] = 0;
	  avEr_TrksBef[i] = 0;
	  avEr_TrksRef[i] = 0;
	  avEr_TrksAft[i] = 0;
  
	  for(int ialc=0;ialc<6;ialc++){
	    av_mis_a_mean[ialc][i] = 0;
	    av_mis_b_mean[ialc][i] = 0;
	    avEr_mis_a_mean[ialc][i] = 0;
	    avEr_mis_b_mean[ialc][i] = 0;
	  }
	}

	for(int i=0;i<nS;i++){
	  for(int ipart=0;ipart<nS;ipart++){
   	    av_theta_mean_b[i]+=theta_mean_b[ipart][i];
	    av_theta_mean_a[i]+=theta_mean_a[ipart][i];
	    av_theta_mean_r[i]+=theta_mean_r[ipart][i];
	    av_theta_rms_b[i]+=theta_rms_b[ipart][i];
	    av_theta_rms_a[i]+=theta_rms_a[ipart][i];
	    av_theta_rms_r[i]+=theta_rms_r[ipart][i];
	    av_TrksBef[i]+=TrksBef[ipart][i];
	    av_TrksRef[i]+=TrksRef[ipart][i];
	    av_TrksAft[i]+=TrksAft[ipart][i];
	    av_TrksDiff[i]+=TrksDiff[ipart][i];
     
	    for(int ialc=0;ialc<6;ialc++){
	      av_mis_a_mean[ialc][i] += mis_a[ipart][ialc][i];
	      av_mis_b_mean[ialc][i] += mis_b[ipart][ialc][i];
	    }
	  }
   
	  av_theta_mean_b[i] *=i_nS;
	  av_theta_mean_a[i] *=i_nS;
	  //   cout<<" av_theta_mean_a["<<i<<"] = "<< av_theta_mean_a[i]<<endl;
	  av_theta_mean_r[i] *=i_nS;
	  av_theta_rms_b[i] *=i_nS;
	  av_theta_rms_a[i] *=i_nS;
	  av_theta_rms_r[i] *=i_nS;
	  av_TrksBef[i]*=i_nS;
	  av_TrksRef[i]*=i_nS;
	  av_TrksAft[i]*=i_nS;
	  if(i>0) av_TrksBef[i] +=av_TrksBef[i-1];
	  //	  cout<<"av_TrksBef["<<i<<"] = "<<av_TrksBef[i]<<endl;
	  // av_TrksDiff[i]*=i_nS;
   
	  for(int ialc=0;ialc<6;ialc++){
	    av_mis_a_mean[ialc][i] *= i_nS;
	    av_mis_b_mean[ialc][i] *= i_nS;
	  }
	}



	for(int i=0;i<nS;i++){
	  for(int ipart=0;ipart<nS;ipart++){
	    avEr_theta_mean_b[i]+=TMath::Power((theta_mean_b[ipart][i]-av_theta_mean_b[i]),2);
	    avEr_theta_mean_a[i]+=TMath::Power((theta_mean_a[ipart][i]-av_theta_mean_a[i]),2);
	    avEr_theta_mean_r[i]+=TMath::Power((theta_mean_r[ipart][i]-av_theta_mean_r[i]),2);
	    avEr_theta_rms_b[i]+=TMath::Power((theta_rms_b[ipart][i]-av_theta_rms_b[i]),2);
	    avEr_theta_rms_a[i]+=TMath::Power((theta_rms_a[ipart][i]-av_theta_rms_a[i]),2);
	    avEr_theta_rms_r[i]+=TMath::Power((theta_rms_r[ipart][i]-av_theta_rms_r[i]),2);
	    // avEr_TrksBef[i]+=TMath::Power((TrksBef[ipart][i]-av_TrksBef[i]),2);
	    // avEr_TrksRef[i]+=TMath::Power((TrksRef[ipart][i]-av_TrksRef[i]),2);
	    // avEr_TrksAft[i]+=TMath::Power((TrksAft[ipart][i]-av_TrksAft[i]),2);
	    // avEr_TrksDiff[i]+=TMath::Power((TrksDiff[ipart][i]-av_TrksDiff[i]),2);
 
	    for(int ialc=0;ialc<6;ialc++){
	      avEr_mis_a_mean[ialc][i] +=TMath::Power((mis_a[ipart][ialc][i]-av_mis_a_mean[ialc][i]),2);
	      avEr_mis_b_mean[ialc][i] +=TMath::Power((mis_b[ipart][ialc][i]-av_mis_b_mean[ialc][i]),2);
	    }
	  }
	  avEr_theta_mean_b[i] *=i_nS_r;
	  avEr_theta_mean_a[i] *=i_nS_r;
	  avEr_theta_mean_r[i] *=i_nS_r;
	  avEr_theta_rms_b[i] *=i_nS_r;
	  avEr_theta_rms_a[i] *=i_nS_r;
	  avEr_theta_rms_r[i] *=i_nS_r;
	  // avEr_TrksBef[i] *= i_nS_r;
	  // avEr_TrksRef[i] *= i_nS_r;
	  // avEr_TrksAft[i] *= i_nS_r;
	  // avEr_TrksDiff[i] *= i_nS_r;
	  avEr_theta_mean_b[i]  = sqrt(avEr_theta_mean_b[i]);
	  avEr_theta_mean_a[i]  = sqrt(avEr_theta_mean_a[i]);
	  avEr_theta_mean_r[i] = sqrt(avEr_theta_mean_r[i]);
	  avEr_theta_rms_b[i] = sqrt(avEr_theta_rms_b[i]);
	  avEr_theta_rms_a[i] = sqrt(avEr_theta_rms_a[i]);
	  avEr_theta_rms_r[i] = sqrt(avEr_theta_rms_r[i]);
	  // avEr_TrksBef[i] = sqrt(avEr_TrksBef[i]);
	  // avEr_TrksRef[i] = sqrt(avEr_TrksRef[i]);
	  // avEr_TrksAft[i] = sqrt(avEr_TrksAft[i]);
	  // avEr_TrksDiff[i] = sqrt(avEr_TrksDiff[i]);
  
	  for(int ialc=0;ialc<6;ialc++){
	    avEr_mis_a_mean[ialc][i] *=i_nS_r;
	    avEr_mis_b_mean[ialc][i] *=i_nS_r;
	    avEr_mis_a_mean[ialc][i] = sqrt(avEr_mis_a_mean[ialc][i]);
	    avEr_mis_b_mean[ialc][i] = sqrt(avEr_mis_b_mean[ialc][i]);
	  }
	}//loop over sum of data samples

  
	for(int i=0;i<nS;i++){
	  misconst_av_theta_mean_b[misAlinSamp][i] = av_theta_mean_b[i];
	  misconst_av_theta_rms_b[misAlinSamp][i] = av_theta_rms_b[i];
	  misconst_av_theta_mean_a[misAlinSamp][i] = av_theta_mean_a[i];
	  misconst_av_theta_rms_a[misAlinSamp][i] = av_theta_rms_a[i];
	  //	  cout<<"av_theta_mean_a[i] = "<<av_theta_mean_a[i]<<endl;
	  //	  cout<<"av_theta_rms_a[i] = "<<av_theta_rms_a[i]<<endl;
	  misconst_av_theta_mean_r[misAlinSamp][i] = av_theta_mean_r[i];
	  misconst_av_theta_rms_r[misAlinSamp][i] = av_theta_rms_r[i];
	  misconst_av_TrksBef[misAlinSamp][i] = av_TrksBef[i];
	  misconst_av_TrksAft[misAlinSamp][i] = av_TrksAft[i];
	  misconst_av_TrksRef[misAlinSamp][i] = av_TrksRef[i];
	  misconst_avEr_theta_mean_b[misAlinSamp][i] = avEr_theta_mean_b[i];
	  misconst_avEr_theta_mean_a[misAlinSamp][i] = avEr_theta_mean_a[i];
	  misconst_avEr_theta_mean_r[misAlinSamp][i] = avEr_theta_mean_r[i];
	  misconst_avEr_theta_rms_b[misAlinSamp][i] = avEr_theta_rms_b[i];
	  misconst_avEr_theta_rms_a[misAlinSamp][i] = avEr_theta_rms_a[i];
	  misconst_avEr_theta_rms_r[misAlinSamp][i] = avEr_theta_rms_r[i];
	  misconst_avEr_TrksBef[misAlinSamp][i] = avEr_TrksBef[i];
	  misconst_avEr_TrksRef[misAlinSamp][i] = avEr_TrksRef[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];
	  }
	}
      }//loop over diff. misalignment samples

      /// Make avarage for different data sums
      double av_misconst_av_theta_mean_b[nS+1],av_misconst_av_theta_rms_b[nS+1];
      double av_misconst_av_theta_mean_a[nS+1],av_misconst_av_theta_rms_a[nS+1];
      double av_misconst_av_theta_mean_r[nS+1],av_misconst_av_theta_rms_r[nS+1];
      double av_misconst_av_TrksBef[nS+1];double av_misconst_av_TrksAft[nS+1];
      double av_misconst_av_TrksRef[nS+1];
      double avEr_misconst_avEr_theta_mean_b[nS+1],avEr_misconst_avEr_theta_rms_b[nS+1];
      double avEr_misconst_avEr_theta_mean_a[nS+1],avEr_misconst_avEr_theta_rms_a[nS+1];
      double avEr_misconst_avEr_theta_mean_r[nS+1],avEr_misconst_avEr_theta_rms_r[nS+1];
      double avEr_misconst_avEr_TrksBef[nS+1]; double avEr_misconst_avEr_TrksAft[nS+1]; double avEr_misconst_avEr_TrksDiff[nS+1]; 
      double avEr_misconst_avEr_TrksRef[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++){
	for(int misAlinSamp=0;misAlinSamp<nAS;misAlinSamp++){
	  av_misconst_av_theta_mean_b[i] += misconst_av_theta_mean_b[misAlinSamp][i];
	  av_misconst_av_theta_mean_a[i] += misconst_av_theta_mean_a[misAlinSamp][i];
	  av_misconst_av_theta_mean_r[i] += misconst_av_theta_mean_r[misAlinSamp][i];
	  av_misconst_av_theta_rms_b[i] += misconst_av_theta_rms_b[misAlinSamp][i];
	  av_misconst_av_theta_rms_a[i] += misconst_av_theta_rms_a[misAlinSamp][i];
	  av_misconst_av_theta_rms_r[i] += misconst_av_theta_rms_r[misAlinSamp][i];
	  av_misconst_av_TrksBef[i] += misconst_av_TrksBef[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;
	//	cout<<"av_misconst_av_TrksBef["<<i<<"] = "<<av_misconst_av_TrksBef[i]<<endl;
	av_misconst_av_theta_mean_b[i] *= i_nAS;
	av_misconst_av_theta_mean_r[i] *= i_nAS;
	av_misconst_av_theta_mean_a[i] *= i_nAS;
	//       	cout<<"av_misconst_av_theta_mean_a["<<i<<"] = "<<av_misconst_av_theta_mean_a[i]<<endl;
	av_misconst_av_theta_rms_b[i] *= i_nAS;
	av_misconst_av_theta_rms_r[i] *= i_nAS;
	av_misconst_av_theta_rms_a[i] *= i_nAS;
       	// cout<<"av_misconst_av_theta_rms_a["<<i<<"] = "<<av_misconst_av_theta_rms_a[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_theta_mean_b[i] += (TMath::Power((misconst_av_theta_mean_b[misAlinSamp][i]-av_misconst_av_theta_mean_b[i]),2)
						 + TMath::Power(misconst_avEr_theta_mean_b[misAlinSamp][i],2));
	  avEr_misconst_avEr_theta_mean_r[i] += (TMath::Power((misconst_av_theta_mean_r[misAlinSamp][i]-av_misconst_av_theta_mean_r[i]),2) 
						 + TMath::Power(misconst_avEr_theta_mean_r[misAlinSamp][i],2));
	  avEr_misconst_avEr_theta_mean_a[i] += (TMath::Power((misconst_av_theta_mean_a[misAlinSamp][i]-av_misconst_av_theta_mean_a[i]),2) 
						 + TMath::Power(misconst_avEr_theta_mean_a[misAlinSamp][i],2));
	  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);

	  //cout<<"misconst_av_theta_mean_a["<<misAlinSamp<<"]["<<i<<"] = "<<misconst_av_theta_mean_a[misAlinSamp][i]<<endl;
	  //cout<<"av_misconst_av_theta_mean_a["<<i<<"] = "<<av_misconst_av_theta_mean_a[i]<<endl;
	  //cout<<"diff = "<<(misconst_av_theta_mean_a[misAlinSamp][i]-av_misconst_av_theta_mean_a[i])<<endl;
	  //cout<<"s1 = "<<TMath::Power((misconst_av_theta_mean_a[misAlinSamp][i]-av_misconst_av_theta_mean_a[i]),2)<<" s2 = "<<TMath::Power(misconst_avEr_theta_mean_a[misAlinSamp][i],2)<<endl;
	  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)
	       +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_theta_mean_b[i] *=i_nAS_r;
	avEr_misconst_avEr_theta_mean_b[i] = sqrt(avEr_misconst_avEr_theta_mean_b[i]);
	avEr_misconst_avEr_theta_mean_a[i] *=i_nAS_r;
	avEr_misconst_avEr_theta_mean_a[i] = sqrt(avEr_misconst_avEr_theta_mean_a[i]);
	//	cout<<" avEr_misconst_avEr_theta_mean_a["<<i<<"] = "<<avEr_misconst_avEr_theta_mean_a[i]<<endl;
	avEr_misconst_avEr_theta_mean_r[i] *=i_nAS_r;
	avEr_misconst_avEr_theta_mean_r[i] = sqrt(avEr_misconst_avEr_theta_mean_r[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;  
      //cout<<"av_Trks = "<<av_Trks<<endl;
      //put together theta resolution results
      TGraphErrors *gr_theta_b = new TGraphErrors(nS,av_misconst_av_TrksBef, av_misconst_av_theta_mean_b,avEr_misconst_avEr_TrksBef, avEr_misconst_avEr_theta_mean_b);
      gr_theta_b->SetMarkerStyle(20); 
      gr_theta_b->SetMarkerColor(kGreen-3); 
      gr_theta_b->SetMarkerSize(2.5); 
      TGraphErrors *gr_theta_a = new TGraphErrors(nS,av_misconst_av_TrksBef, av_misconst_av_theta_mean_a,avEr_misconst_avEr_TrksBef, avEr_misconst_avEr_theta_mean_a);
      gr_theta_a->SetMarkerStyle(29); 
      gr_theta_a->SetMarkerColor(kOrange+7); 
      gr_theta_a->SetMarkerSize(2.5); 
      TGraphErrors *gr_theta_r = new TGraphErrors(nS,av_misconst_av_TrksBef, av_misconst_av_theta_mean_r,avEr_misconst_avEr_TrksBef, avEr_misconst_avEr_theta_mean_r);
      gr_theta_r->SetMarkerStyle(21); 
      gr_theta_r->SetMarkerColor(15); 
      gr_theta_r->SetMarkerSize(2.5); 


      TMultiGraph *mgr_theta = new TMultiGraph();
      TMultiGraph *mgr_theta_b = new TMultiGraph();
      mgr_theta_b->Add(gr_theta_b);
      mgr_theta_b->Add(gr_theta_r);
      mgr_theta_b->Draw("AP");
      mgr_theta_b->GetXaxis()->SetTitle("Number of rec. trks used for alignment");
      mgr_theta_b->GetYaxis()->SetTitle("#theta_{mean}, #murad");
      c1.Print(resname_pdf_o); //write canvas and keep the ps file open
      c1.Clear();
      mgr_theta->Add(gr_theta_r);
      mgr_theta->Add(gr_theta_a);
      mgr_theta->Draw("AP");
      mgr_theta->GetXaxis()->SetTitle("Number of rec. trks used for alignment");
      mgr_theta->GetYaxis()->SetTitle("#theta_{mean}, #murad");
      c1.Print(resname_pdf_o); //write canvas and keep the ps file open
      c1.Clear();

      TGraphErrors *gr_theta_rms_b = new TGraphErrors(nS,av_misconst_av_TrksBef,av_misconst_av_theta_rms_b,avEr_misconst_avEr_TrksBef,avEr_misconst_avEr_theta_rms_b);
      gr_theta_rms_b->SetMarkerStyle(20); 
      gr_theta_rms_b->SetMarkerColor(kGreen-3); 
      gr_theta_rms_b->SetMarkerSize(2.5); 
      TGraphErrors *gr_theta_rms_a = new TGraphErrors(nS,av_misconst_av_TrksBef,av_misconst_av_theta_rms_a,avEr_misconst_avEr_TrksBef,avEr_misconst_avEr_theta_rms_a);
      gr_theta_rms_a->SetMarkerStyle(29); 
      gr_theta_rms_a->SetMarkerColor(kOrange+7); 
      gr_theta_rms_a->SetMarkerSize(2.5); 
      TGraphErrors *gr_theta_rms_r = new TGraphErrors(nS,av_misconst_av_TrksBef,av_misconst_av_theta_rms_r,avEr_misconst_avEr_TrksBef,avEr_misconst_avEr_theta_rms_r);
      gr_theta_rms_r->SetMarkerStyle(21); 
      gr_theta_rms_r->SetMarkerColor(15); 
      gr_theta_rms_r->SetMarkerSize(2.5); 
      TMultiGraph *mgr_theta_rms = new TMultiGraph();
      TMultiGraph *mgr_theta_rms_b = new TMultiGraph();
      mgr_theta_rms_b->Add(gr_theta_rms_b);
      mgr_theta_rms_b->Add(gr_theta_rms_r);
      mgr_theta_rms->Add(gr_theta_rms_r);
      mgr_theta_rms->Add(gr_theta_rms_a);
      mgr_theta_rms_b->Draw("AP");
      mgr_theta_rms_b->GetXaxis()->SetTitle("Number of rec. trks used for alignment");
      mgr_theta_rms_b->GetYaxis()->SetTitle("#theta_{#sigma}, #murad");
      c1.Print(resname_pdf_o); //write canvas and keep the ps file open
      c1.Clear();
      mgr_theta_rms->Draw("AP");
      mgr_theta_rms->GetXaxis()->SetTitle("Number of rec. trks used for alignment");
      mgr_theta_rms->GetYaxis()->SetTitle("#theta_{#sigma}, #murad");
      c1.Print(resname_pdf_o); //write canvas and keep the ps file open
      c1.Clear();

      //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++){
        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);
        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] = new TMultiGraph();
        //  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 ");
	mgr_mis[ialc]->GetHistogram()->SetMaximum(3.);
	mgr_mis[ialc]->GetHistogram()->SetMinimum(-3.);
        if(ialc==0) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{x}, #mum");
        if(ialc==1) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{y}, #mum");
        if(ialc==2) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{z}, #mum");
        if(ialc==3) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#alpha}, mrad");
        if(ialc==4) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#beta}, mrad");
        if(ialc==5) mgr_mis[ialc]->GetYaxis()->SetTitle("#Delta_{#gamma}, mrad");
        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 *grstatBef = new TGraphErrors(nS,av_misconst_av_TrksBef,av_TrksBef,avEr_misconst_avEr_TrksBef,avEr_TrksBef);
      //  grstatBef->SetMarkerStyle(24);
      //  grstatBef->SetMarkerSize(2.5);
      //  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");

      //  gr_theta_b->Write();
      //  gr_theta_r->Write();
      //  gr_theta_a->Write();
      //  mgr_theta->Write();
      //  gr_theta_rms_b->Write();
      //  gr_theta_rms_r->Write();
      //  gr_theta_rms_a->Write();
      //  mgr_theta_rms->Write();

      //  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;
}