// Macro to check Point Resolution reading the TTree 
//04/10/2010

//xfit: Point Res XZ proj.
//yfit: Point Res YZ proj.
//sigmax: Point Res Pull XZ proj.
//sigmay: Point Res Pull YZ proj.

/*************************************************************************
You have to provide to the macro:

<Input FILE LIST of ROOT files (.txt)>
<Output ROOT file (.root) with the histos> 
<DIR PATH where you want only the output plots>
 *************************************************************************/

#include <iostream>
#include <TH1F.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TBranch.h>
#include <TMath.h>
#include <TString.h>
#include <TList.h>
#include <TLine.h>


using namespace std;



void ResolvsPane(TH2F *h_PointResvsPlane, TString histoname, char *outpath);
void Mean_Sigma_Trend_Map_LASER(TH2F *h_PointRes, TString histoname, char *path, int Zbins, int Ybins);


void pointres(char *infile, char *outfile, char *outpath)
{
  gROOT->Reset();
  
  
  if(infile == 0)
  {
    printf("pointres_for_alignment_check: \n");
    printf("<Input FILE LIST of ROOT files (.txt)>, <Output ROOT file (.root) with the histos> <DIR PATH where you want only the output plots> \n"); 

    return;

   }

  // out
  TFile *file2 = new TFile(outfile, "RECREATE");

  ////TRANSVERSE
  TH2F *h_PointResYvsPlane_EndplatePositive_2Pads = new TH2F("h_PointResYvsPlane_EndplatePositive_2Pads"," Trasv Point Res vs Plane 2Pads (EndPlate +)", 216,-0.5,215.5,100,-5.,5.);
  TH2F *h_PointResYvsPlane_EndplateNegative_2Pads = new TH2F("h_PointResYvsPlane_EndplateNegative_2Pads"," Trasv Point Res vs Plane 2Pads (EndPlate -)", 216,-0.5,215.5,100,-5.,5.);

  TH2F *h_TrasvPoinResMap_LASER_EndplatePositive_TPC1 = new TH2F("h_TrasvPoinResMap_LASER_EndplatePositive_TPC1"," Trasv Point Res LASER Map EndPlane+  TPC1", 84,-0.5,83.5,100,-5.,5.);
  TH2F *h_TrasvPoinResMap_LASER_EndplatePositive_TPC2 = new TH2F("h_TrasvPoinResMap_LASER_EndplatePositive_TPC2"," Trasv Point Res LASER Map EndPlane+ TPC2", 84,-0.5,83.5,100,-5.,5.);
  TH2F *h_TrasvPoinResMap_LASER_EndplatePositive_TPC3 = new TH2F("h_TrasvPoinResMap_LASER_EndplatePositive_TPC3"," Trasv Point Res LASER Map EndPlane+ TPC3", 84,-0.5,83.5,100,-5.,5.);

  TH2F *h_TrasvPoinResMap_LASER_EndplateNegative_TPC1 = new TH2F("h_TrasvPoinResMap_LASER_EndplateNegative_TPC1"," Trasv Point Res LASER Map EndPlane- TPC1", 84,-0.5,83.5,100,-5.,5.);
  TH2F *h_TrasvPoinResMap_LASER_EndplateNegative_TPC2 = new TH2F("h_TrasvPoinResMap_LASER_EndplateNegative_TPC2"," Trasv Point Res LASER Map EndPlane- TPC2", 84,-0.5,83.5,100,-5.,5.);
  TH2F *h_TrasvPoinResMap_LASER_EndplateNegative_TPC3 = new TH2F("h_TrasvPoinResMap_LASER_EndplateNegative_TPC3"," Trasv Point Res LASER Map EndPlane- TPC3", 84,-0.5,83.5,100,-5.,5.);

  
  //LONGITUDINAL
  TH2F *h_PointResXvsPlane_EndplatePositive_2Pads = new TH2F("h_PointResXvsPlane_EndplatePositive_2Pads"," Long Point Res vs Plane 2Pads (EndPlate +)", 216,-0.5,215.5,100,-5.,5.);
  TH2F *h_PointResXvsPlane_EndplateNegative_2Pads = new TH2F("h_PointResXvsPlane_EndplateNegative_2Pads"," Long Point Res vs Plane 2Pads (EndPlate -)", 216,-0.5,215.5,100,-5.,5.);
 

  ///////////VARIABLE
  // Retrive variables
  double xtrack;
  double ytrack; 
  double ztrack; 
  double txtrack;
  double tytrack;
  double rho;
  double pt;
  double t0;
  double trasvdiff;
  int nplanes;

  const int maxplanes = 216;
  int minplanes = 60;

  double *x = new double[maxplanes];
  double *xfit= new double[maxplanes];
  double *sigmax = new double[maxplanes];
  double *y = new double[maxplanes]; 
  double *yfit = new double[maxplanes];
  double *sigmay = new double[maxplanes];
  double *z = new double[maxplanes];
  double *tx = new double[maxplanes];
  double *ty = new double[maxplanes];
  double *q = new double[maxplanes];
  double *deltaDrift = new double[maxplanes];
  int *idplane = new int[maxplanes];
  int *npads = new int[maxplanes];
  int *enabled = new int[maxplanes];

  ///////////////////////////////////////////////////////////////////////
  
   //////////////////////////
  //Map of point res

  const int nbinsZ=14;
  const int nbinsX=4;
  const int nbinsY=6;

  double Xmin=-900.0;
  double Xmax=900.0;
  double Ymin=-1000.0;
  double Ymax=1000.0;

  
  double Z_TPC1_min=-738.0;
  double Z_TPC1_max=-26.0;
  double Z_TPC2_min=620.0;
  double Z_TPC2_max=1332.0;
  double Z_TPC3_min=1979.0;
  double Z_TPC3_max=2691.0;

  //////////////////////////////////////////
  // LASER dots GEOMETRY
  float TPC1_RP0_dots_z[nbinsZ] = {-30.92,-79.87,-128.94,-187.92,-246.64,-295.57,-344.61,-410.39,-459.31,-508.36,-567.25,-625.97,-674.86,-723.82};
  float TPC1_RP1_dots_z[nbinsZ] = {-724.03,-675.05,-626.0,-567.14,-508.41,-459.48,-410.36,-344.47,-295.51,-246.5,-187.63,-128.90,-80.0,-30.96};
  float TPC2_RP0_dots_z[nbinsZ] = {1326.47,1277.96,1229.42,1171.29,1113.11,1064.61,1016.13,946.98,898.48,849.96,791.78,733.59,685.14,636.62};
  float TPC2_RP1_dots_z[nbinsZ] = {636.56,685.06,733.58,791.8,850.,898.54,947.03,1016.11,1064.61,1113.09,1171.27,1229.45,1277.95,1326.47};
  float TPC3_RP0_dots_z[nbinsZ] = {2687.13,2638.19,2589.12,2530.24,2471.50,2422.53,2373.45,2307.59,2258.64,2209.6,2150.7,2091.93,2043.,1993.95};
  float TPC3_RP1_dots_z[nbinsZ] = {1194.04,2042.97,2092.01,2150.87,2209.57,2258.54,2307.63,2373.50,2422.41,2471.45,2530.31,2589.11,2638.04,2687.06};

  float TPC1_RP0_dots_z_edgs_high[nbinsZ+1] = {-26.,-55.39,-104.40,-158.43,-217.28,-271.11,-320.09,-377.50,-434.85,-483.84,-537.80,-596.61,-650.42,-699.34,-738.};
  float TPC1_RP1_dots_z_edgs_low[nbinsZ+1] = {-738.0,-699.53,-650.52,-596.57,-537.77,-483.94,-434.92,-377.42,-319.99,-271.,-217.06,-158.26,-104.45,-55.48,-26.};
  float TPC2_RP0_dots_z_edges_low[nbinsZ+1] = {1332.,1302.21,1253.69,1200.36,1142.2,1088.86,1040.37,981.55,922.73,874.22,820.87,762.68,709.36,660.88,620.};
  float TPC2_RP1_dots_z_edges_high[nbinsZ+1] = {620.,660.81,709.33,762.69,820.9,874.27,922.78,981.57,1040.36,1088.85,1142.18,1200.36,1253.70,1302.21,1332.};
  float TPC3_RP0_dots_z_edge_low[nbinsZ+1] = {2691.,2662.66,2613.66,2559.68,2500.87,2447.01,2397.99,2340.52,2283.12,2234.12,2180.15,2121.32,2067.47,2018.48,1979.};
  float TPC3_RP1_dots_z_edge_high[nbinsZ+1] = {1979.,2018.50,2067.49,2121.44,2180.22,2234.05,2283.09,2340.57,2397.96,2446.93,2500.88,2559.71,2613.57,2662.55,2691.};

  int binx;
  int biny;
  int binz;

  double xbin;
  double ybin;
  double zbin;

  int index;
  int index_LASER;

  //////////////////////////////////////////////////////////////////////////////////
  
  //List of .root files
   int check;
   int filenum = 0;
   FILE *lr;
   char rootfilename[100];

   TList *list_files = new TList();

   if((lr=fopen(infile, "r")) != NULL)
   {
     //TList *list_files = new TList();
      while((check=feof(lr)) == 0)
      {
   fscanf(lr,"%s \n", rootfilename);
   TFile *namefile = new TFile(rootfilename);
   list_files->Add(namefile);
   filenum++;
   printf("%d \n", filenum);
       }
       fclose(lr);
       
       
    }
    else 
    {
      printf("List files not opened! \n");
      return ;
    }

    int nfile = (int)list_files->GetSize();
    printf("root files uploaded: %d \n", nfile);

  
 for(int m=0; m<nfile; m++)
 {
 
  printf("%d \n", m+1);

  TFile *fileTree = (TFile*)list_files->At(m);
  
  TTree *pointrestree = (TTree*)fileTree->Get("pointrestree"); 
 
  pointrestree->SetBranchAddress("xtrack", &xtrack);
  pointrestree->SetBranchAddress("ytrack", &ytrack);
  pointrestree->SetBranchAddress("ztrack", &ztrack);
  pointrestree->SetBranchAddress("txtrack", &txtrack);
  pointrestree->SetBranchAddress("tytrack", &tytrack);
  pointrestree->SetBranchAddress("rho", &rho);
  pointrestree->SetBranchAddress("pt", &pt);
  pointrestree->SetBranchAddress("t0", &t0);
  pointrestree->SetBranchAddress("trasvdiff", &trasvdiff);
  pointrestree->SetBranchAddress("nplanes", &nplanes);
 
  pointrestree->SetBranchAddress("x",x);
  pointrestree->SetBranchAddress("xfit",xfit);
  pointrestree->SetBranchAddress("sigmax",sigmax);
  pointrestree->SetBranchAddress("y",y);
  pointrestree->SetBranchAddress("yfit",yfit);
  pointrestree->SetBranchAddress("sigmay",sigmay);
  pointrestree->SetBranchAddress("z",z);
  pointrestree->SetBranchAddress("tx",tx);
  pointrestree->SetBranchAddress("ty",ty);
  pointrestree->SetBranchAddress("q",q);
  pointrestree->SetBranchAddress("deltaDrift",deltaDrift);
  pointrestree->SetBranchAddress("idplane",idplane);
  pointrestree->SetBranchAddress("npads",npads);
  pointrestree->SetBranchAddress("enabled",enabled);
  ////////////////////////////////////////////////////////7
  
  
  int nentries = pointrestree->GetEntries();
  printf("Tree entries: %d \n", nentries);
  
  //Loop over events..
  for(int i=0; i<nentries;i++)
  {
  //printf(" entry = %d \n", i);
    pointrestree->GetEntry(i);
    //printf("nplanes = %d \n", nplanes);
    
    if(i%1000 == 0)
    std::cout <<  "Events Processed " << i << std::endl;


    if(nplanes < minplanes) continue; // min number of planes for fit
    //printf("nplanes = %d \n", nplanes);
  
    // CATHODE OR NO CATHODE TRACKS
    bool CathodeTrack = false; 
    double xcoord = 0;

    double QMax = 0;

    bool endplate_positive = false; 
    bool endplate_negative = false; 
   
    bool disabledTrack = false;

    // Selection of NON CATHODE tracks
    for(int j=0; j<nplanes; j++)
    {
         if(xcoord*x[j] < 0) {CathodeTrack = true;break;}
         if((xcoord*x[j] > 0) && x[j] <0) {endplate_negative = true;}
         if((xcoord*x[j] > 0) && x[j] >0) {endplate_positive = true;}
     xcoord = x[j];
   
    }

    if(t0 == 0) continue; 

    for(int npl=0; npl<nplanes; npl++)
    {
       
       if(npads[npl] == 2) 
       {

    //Select EndPlate+/-of  TPC
   if(endplate_positive) 
   {
     h_PointResYvsPlane_EndplatePositive_2Pads->Fill(idplane[npl],yfit[npl]);
     h_PointResXvsPlane_EndplatePositive_2Pads->Fill(idplane[npl],xfit[npl]);

         }
   
    if(endplate_negative)
    {
      h_PointResYvsPlane_EndplateNegative_2Pads->Fill(idplane[npl],yfit[npl]);
      h_PointResXvsPlane_EndplateNegative_2Pads->Fill(idplane[npl],xfit[npl]);
          }

       }

       /////////////////////////////////////////////////////////// MAP OF POINT RES
       
       binx = (int)(((x[npl]-Xmin)/(Xmax-Xmin))*nbinsX);
       biny = (int)(((y[npl]-Ymin)/(Ymax-Ymin))*nbinsY);

  
       // TPC1
       if((z[npl]>-738) && (z[npl]<-26))
       {
     
   if(endplate_positive)
   {  
    /////////////////////////////////////////Perform binz
    for(int k=0; k<nbinsZ; k++) 
    {
       if((z[npl] > TPC1_RP1_dots_z_edgs_low[k]) && (z[npl] < TPC1_RP1_dots_z_edgs_low[k+1])) //order is increasing!!!!
       {
         binz = k;
         break;
       }
    
    }
    
    index_LASER = biny*nbinsZ + binz;
  
    h_TrasvPoinResMap_LASER_EndplatePositive_TPC1->Fill(index_LASER,yfit[npl]);
    
   }
   if(endplate_negative)
   {	 
    /////////////////////////////////////////Perform binz
    for(int k=0; k<nbinsZ; k++) 
    {
      if((z[npl] < TPC1_RP0_dots_z_edgs_high[k]) && (z[npl] > TPC1_RP0_dots_z_edgs_high[k+1])) //order is decreasing!!!!
      {
        binz = k;
        break;
      }
      
    }
    index_LASER = biny*nbinsZ + binz;
  
    h_TrasvPoinResMap_LASER_EndplateNegative_TPC1->Fill(index_LASER,yfit[npl]);
          
          
   }	
       }

       // TPC2
       if((z[npl]>620) && (z[npl]<1332))
       {
   if(endplate_positive)
   {  
     /////////////////////////////////////////Perform binz
     for(int k=0; k<nbinsZ; k++) 
     {
       if((z[npl] > TPC2_RP1_dots_z_edges_high[k]) && (z[npl] < TPC2_RP1_dots_z_edges_high[k+1])) //order is increasing!!!!
       {
         binz = k;
         break;
       }
     }
     
     index_LASER = biny*nbinsZ + binz;
     
     h_TrasvPoinResMap_LASER_EndplatePositive_TPC2->Fill(index_LASER,yfit[npl]);
   }
   if(endplate_negative)
   {   
     /////////////////////////////////////////Perform binz
     for(int k=0; k<nbinsZ; k++) 
     {
       if((z[npl] <TPC2_RP0_dots_z_edges_low [k]) && (z[npl] > TPC2_RP0_dots_z_edges_low[k+1])) //order is decreasing!!!!
         {
     binz = k;
     break;
         }
       
     }
     
     index_LASER = biny*nbinsZ + binz;
     
     h_TrasvPoinResMap_LASER_EndplateNegative_TPC2->Fill(index_LASER,yfit[npl]);
     
   }	  
         
       }

       // TPC3
       if((z[npl]>1979) && (z[npl]<2691))
       {
     if(endplate_positive)
     {  
       /////////////////////////////////////////Perform binz
       for(int k=0; k<nbinsZ; k++) 
       {
         if((z[npl] > TPC3_RP1_dots_z_edge_high[k]) && (z[npl] < TPC3_RP1_dots_z_edge_high[k+1])) //order is increasing!!!!
         {
     binz = k;
     break;
         }
       
       }
       index_LASER = biny*nbinsZ + binz;
       
       h_TrasvPoinResMap_LASER_EndplatePositive_TPC3->Fill(index_LASER,yfit[npl]);
       
     }
         
     if(endplate_negative)
     {   
       /////////////////////////////////////////Perform binz
       for(int k=0; k<nbinsZ; k++) 
       {
         if((z[npl] < TPC3_RP0_dots_z_edge_low[k]) && (z[npl] > TPC3_RP0_dots_z_edge_low[k+1])) //order is decreasing!!!!
         {
     binz = k;
     break;
         }
    
       }
          
       index_LASER = biny*nbinsZ + binz;

       h_TrasvPoinResMap_LASER_EndplateNegative_TPC3->Fill(index_LASER,yfit[npl]);
       
     }   
        
       }                  
     //////////////////////////////////////////////////////////////////////////////////////////////////////
    }

  }

 } 
  file2->Write();

  char *a = "Y_EndplatePositive";
  TString name(a);
  ResolvsPlane(h_PointResYvsPlane_EndplatePositive_2Pads, name, outpath);

  char *a = "Y_EndplateNegative";
  TString name(a);
  ResolvsPlane(h_PointResYvsPlane_EndplateNegative_2Pads, name, outpath);
 
  /*
  char *a = "X_EndplatePositive";
  TString name(a);
  ResolvsPlane(h_PointResXvsPlane_EndplatePositive_2Pads, name, outpath);

  char *a = "X_EndplateNegative";
  TString name(a);
  ResolvsPlane(h_PointResXvsPlane_EndplateNegative_2Pads, name, outpath);   
  */

   ////////////// LASER
  char *a = "PointRes_Trasv_LASER_EndplatePositive_TPC1";
  TString name(a);
  Mean_Sigma_Trend_Map_LASER(h_TrasvPoinResMap_LASER_EndplatePositive_TPC1, name, outpath, nbinsZ, nbinsY);

  char *a = "PointRes_Trasv_LASER_EndplateNegative_TPC1";
  TString name(a);
  Mean_Sigma_Trend_Map_LASER(h_TrasvPoinResMap_LASER_EndplateNegative_TPC1, name, outpath, nbinsZ, nbinsY);
  
  char *a = "PointRes_Trasv_LASER_EndplatePositive_TPC2";
  TString name(a);
  Mean_Sigma_Trend_Map_LASER(h_TrasvPoinResMap_LASER_EndplatePositive_TPC2, name, outpath, nbinsZ, nbinsY);
  
  char *a = "PointRes_Trasv_LASER_EndplateNegative_TPC2";
  TString name(a);
  Mean_Sigma_Trend_Map_LASER(h_TrasvPoinResMap_LASER_EndplateNegative_TPC2, name, outpath, nbinsZ, nbinsY);
  
  char *a = "PointRes_Trasv_LASER_EndplatePositive_TPC3";
  TString name(a);
  Mean_Sigma_Trend_Map_LASER(h_TrasvPoinResMap_LASER_EndplatePositive_TPC3, name, outpath, nbinsZ, nbinsY);

  char *a = "PointRes_Trasv_LASER_EndplateNegative_TPC3";
  TString name(a);
  Mean_Sigma_Trend_Map_LASER(h_TrasvPoinResMap_LASER_EndplateNegative_TPC3, name, outpath, nbinsZ, nbinsY);
} 


////////////////// Point Res vs Y trend
void ResolvsPlane(TH2F *h_PointResvsPlane, TString histoname, char *outpath)
{
 char filename_out[256];
 sprintf(filename_out,"%s/Mean_Sigma_Trend_%s.png", outpath, histoname.Data());

 int nbins =  h_PointResvsPlane->GetNbinsX();
 //printf(" nbins = %d \n", nbins);
 
 double *mean,*emean,*sigma,*esigma,*b,*eb,*skewness,*kurtosis,*enull;

 mean = new double[nbins]; 
 emean = new double[nbins]; 
 sigma = new double[nbins]; 
 esigma = new double[nbins]; 
 b = new double[nbins];
 eb = new double[nbins];
 skewness = new double[nbins];
 kurtosis = new double[nbins];
 enull = new double[nbins];

 // printf("I'm here \n");

 for(int i = 0; i < nbins; i++)
 {  
   char name[256];

   sprintf(name,"proy_%d",i);
   //cout << " name:  " << name << endl;
   
   b[i] = i; 
   eb[i] = 0.5;
   enull[i] = 0.0; 
   
   TH1D *h = (TH1D*) h_PointResvsPlane->ProjectionY(name,i+1,i+1);
   //TH1D *h = h_PointResvsy->ProjectionY(name,i+1,i+1);
   

   if(h->GetEntries() < 10) continue;

   
   h->Fit("gaus","","",-2.,2.); 
   
   TF1 *g = h->GetFunction("gaus");

   mean[i] = g->GetParameter(1);
   emean[i] = g->GetParError(1);
   sigma[i] = g->GetParameter(2);
   esigma[i] = g->GetParError(2);

   h->Fit("gaus","","",mean[i]-2.*sigma[i],mean[i]+2.*sigma[i]);
   
   g = h->GetFunction("gaus");

   mean[i] = g->GetParameter(1);
   emean[i] = g->GetParError(1);
   sigma[i] = g->GetParameter(2);
   esigma[i] = g->GetParError(2);

   skewness[i] = h->GetSkewness();
   kurtosis[i] = h->GetKurtosis();

   

 }

 
 cout << " Results from ResolvsPlane" << endl; 

 TGraphErrors *MEAN = new TGraphErrors(nbins,b,mean,eb,emean);
 TGraphErrors *SIGMA = new TGraphErrors(nbins,b,sigma,eb,esigma);
 TGraphErrors *SKEW = new TGraphErrors(nbins,b,skewness,eb,enull);
 TGraphErrors *KURT = new TGraphErrors(nbins,b,kurtosis,eb,enull);


 TH2F *tmpMean = new TH2F("tmpMean"," Mean ",100,-0.5,nbins+0.5,100,-1.,1.); 
 TH2F *tmpSigma = new TH2F("tmpSigma"," Sigma ",100,-0.5,nbins+0.5,100,0.,3.); 
 TH2F *tmpSkew = new TH2F("tmpSkew"," Skewness ",100,-0.5,nbins+0.5,100,-1.,1.); 
 TH2F *tmpKurt = new TH2F("tmpKurt"," Kurtosis ",100,-0.5,nbins+0.5,100,0.,3.); 
 
 TCanvas *c_resolvsplane = new TCanvas("c_resolvsplane"," "); 
 
 
 gStyle->SetOptStat(0);

 TLine *ln1 = new TLine(71.5,-3.,71.5,3.);
 TLine *ln2 = new TLine(143.5,-3.,143.5,3.);

 ln1->SetLineColor(2);
 ln2->SetLineColor(2);

 c_resolvsplane ->Divide(2,1); 

 c_resolvsplane ->cd(1); 

 tmpMean->Draw(); 

 MEAN->Draw("same"); 

 ln1->Draw();
 ln2->Draw();

 c_resolvsplane ->cd(2); 

 tmpSigma->Draw(); 

 SIGMA->Draw("same"); 

 ln1->Draw();
 ln2->Draw();

 /*
 
 c_resolvsplane ->cd(3); 

 tmpSkew->Draw(); 

 SKEW->Draw("same"); 

 //ln1->Draw();
 //ln2->Draw();

 c_resolvsplane ->cd(4); 

 tmpKurt->Draw(); 

 KURT->Draw("same"); 

 //ln1->Draw();
 //ln2->Draw();
 
 */

 c_resolvsplane ->Update();
 c_resolvsplane ->SaveAs(filename_out);

   
}

void Mean_Sigma_Trend_Map_LASER(TH2F *h_PointRes, TString histoname, char *path, int Zbins, int Ybins)
{
 char filename_out[256];

 sprintf(filename_out,"%s/Mean_Sigma_Trend_%s.png", path, histoname.Data());
 
 int nbins =  h_PointRes->GetNbinsX();
 //printf(" nbins = %d \n", nbins);
 
 double *mean,*emean,*sigma,*esigma,*b,*eb,*skewness,*kurtosis,*enull;

 mean = new double[nbins]; 
 emean = new double[nbins]; 
 sigma = new double[nbins]; 
 esigma = new double[nbins]; 
 b = new double[nbins];
 eb = new double[nbins];
 skewness = new double[nbins];
 kurtosis = new double[nbins];
 enull = new double[nbins];

 // printf("I'm here \n");

 for(int i = 0; i < nbins; i++)
 {  
   char name[256];

   sprintf(name,"proy_%d",i);
   //cout << " name:  " << name << endl;
   
   b[i] = i; 
   eb[i] = 0.5;
   enull[i] = 0.0; 
   
   TH1D *h = (TH1D*) h_PointRes->ProjectionY(name,i+1,i+1);
   //TH1D *h = h_PointResvsy->ProjectionY(name,i+1,i+1);
   

   if(h->GetEntries() < 10) continue;

   
   h->Fit("gaus","","",-2.,2.); 
   
   TF1 *g = h->GetFunction("gaus");

   mean[i] = g->GetParameter(1);
   emean[i] = g->GetParError(1);
   sigma[i] = g->GetParameter(2);
   esigma[i] = g->GetParError(2);

   h->Fit("gaus","","",mean[i]-2.*sigma[i],mean[i]+2.*sigma[i]);
   
   g = h->GetFunction("gaus");

   mean[i] = g->GetParameter(1);
   emean[i] = g->GetParError(1);
   sigma[i] = g->GetParameter(2);
   esigma[i] = g->GetParError(2);

   skewness[i] = h->GetSkewness();
   kurtosis[i] = h->GetKurtosis();

   h->Write();

 }

 
 cout << " Results from ResolvsPlane" << endl; 

 TGraphErrors *MEAN = new TGraphErrors(nbins,b,mean,eb,emean);
 TGraphErrors *SIGMA = new TGraphErrors(nbins,b,sigma,eb,esigma);
 TGraphErrors *SKEW = new TGraphErrors(nbins,b,skewness,eb,enull);
 TGraphErrors *KURT = new TGraphErrors(nbins,b,kurtosis,eb,enull);


 TH2F *tmpMean = new TH2F("tmpMean"," Mean ",100,-0.5,nbins+0.5,100,-1.,1.); 
 TH2F *tmpSigma = new TH2F("tmpSigma"," Sigma ",100,-0.5,nbins+0.5,100,0.,3.); 
 TH2F *tmpSkew = new TH2F("tmpSkew"," Skewness ",100,-0.5,nbins+0.5,100,-1.,1.); 
 TH2F *tmpKurt = new TH2F("tmpKurt"," Kurtosis ",100,-0.5,nbins+0.5,100,0.,3.); 
 
 TCanvas *c_resolvsplane = new TCanvas("c_resolvsplane"," "); 
 
 
 gStyle->SetOptStat(0);

 c_resolvsplane ->Divide(1,2); 
 //c_resolvsplane ->Divide(2,1); 

 c_resolvsplane ->cd(1); 

 tmpMean->Draw(); 

 MEAN->Draw("same"); 

 for(int j=0; j<Ybins; j++)
 {  
  for(int i=0; i<Zbins; i++)
  { 
    TLine *ln = new TLine(0.5+i+Zbins*j,-3.,0.5+i+Zbins*j,3.);
    if(i == Zbins-1) //only green line!
    { 
      ln->SetLineColor(3);
       ln->Draw();
    }
    //ln->Draw(); //all lines, black and green
  }

 }

 c_resolvsplane ->cd(2); 

 tmpSigma->Draw(); 

 SIGMA->Draw("same"); 

 for(int j=0; j<Ybins; j++)
 {  
   for(int i=0; i<Zbins; i++)
   {
     TLine *ln = new TLine(0.5+i+Zbins*j,-3.,0.5+i+Zbins*j,3.);
     if(i == Zbins-1) //only green line!
     { 
       ln->SetLineColor(3);
       ln->Draw();
     }
     //ln->Draw(); //all lines, black and green
   }

 }

 /*
 c_resolvsplane ->cd(3); 

 tmpSkew->Draw(); 

 SKEW->Draw("same"); 

 //ln1->Draw();
 //ln2->Draw();

 c_resolvsplane ->cd(4); 

 tmpKurt->Draw(); 

 KURT->Draw("same"); 

 //ln1->Draw();
 //ln2->Draw();
 
 */

 c_resolvsplane ->Update();
 c_resolvsplane ->SaveAs(filename_out);

}