#define TBAlign_cxx
#include "TBAlign.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TH1.h>
//#include "fit_line.C"
#include "TF1.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH3.h"
#include "TGaxis.h"
#include "TRandom.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TGraph.h"
#include "TSystem.h"
#include "TTree.h"
#include "TTreePlayer.h"
#include "TCut.h"
#include "TPaletteAxis.h"
//
#include <iostream>
//
// Size of the Si trackers
Float_t L10(1.91), L20(1.91), L30(9.50), L31(9.50);
Float_t L32(8.0),L33(8.0),L34(9.5),L35(9.5),L36(18.0),L37(18.0);
// 
//

  //----------------START FACA----------------------------------------
   float x[3];					  // x coordinates of hits 
   float z[] = {1873.93, 1995.43, 2125.43}; // z coordinates of 33, 35, 37
   float ex[] = {0.1, 0.1, 0.1}; // errors of x coordinates of hits
   //----------------END FACA----------------------------------------
  
  
void TBAlign::Loop()
{
	//----------------START FACA----------------------------------------
	int p1p233=2;							//Si33
	int p1p235=1;							//Si35	
	int p1p237=2;							//Si37
	float p133;
	float p233;
	float p135;
	float p235;
	float p137;
	float p237;
 //----------------END FACA----------------------------------------
  
   if (fChain == 0) return;

   TFile *newfile = new TFile("C:/si-calo-sep18/si-500596-trasl.root","recreate");
   TTree *newtree = fChain->CloneTree(0);

   // -----------------------------------------------------------------------------
   // If Lab==1, all coordinated x10, x20... will be trasformed in Lab frame
   // otherwise they will be in the local frame (or as will be in xh[1], yh[]
   // -----------------------------------------------------------------------------
   Int_t Lab = 1;

// Control Histogram  to be filled


  TH1F *hnx10=new TH1F("hnx10","N hits in Si10, X coord",10,0,10);
  TH1F *hnx20=new TH1F("hnx20","N hits in Si20, X coord",10,0,10);
  TH1F *hnx30=new TH1F("hnx30","N hits in Si30, X coord",10,0,10);
  TH1F *hnx31=new TH1F("hnx31","N hits in Si31, X coord",10,0,10);
  //----------------START FACA----------------------------------------
  TH1F *hnx33=new TH1F("hnx33","N hits in Si33, X coord",10,0,10);
  TH1F *hnx35=new TH1F("hnx35","N hits in Si35, X coord",10,0,10);
  TH1F *hnx37=new TH1F("hnx37","N hits in Si37, X coord",10,0,10);
  //----------------END FACA----------------------------------------

  TH1F *hx10=new TH1F("hx10","X cood, Si10", 100,0,L10);
  TH1F *hx20=new TH1F("hx20","X cood, Si20", 100,0,L20);
  TH1F *hx30=new TH1F("hx30","X cood, Si30", 100,0,L30);
  TH1F *hx31=new TH1F("hx31","X cood, Si31", 100,0,L31);
  
  //----------------START FACA----------------------------------------
  //TH1F *hx33=new TH1F("hx33","X cood, Si33", 100,0,L33);
  //TH1F *hx35=new TH1F("hx35","X cood, Si35", 100,0,L35);
  //TH1F *hx37=new TH1F("hx37","X cood, Si37", 100,0,L37);
  TH1F *hx33=new TH1F("hx33","X cood, Si33", 100,0.,0.);
  TH1F *hx35=new TH1F("hx35","X cood, Si35", 100,0.,0.);
  TH1F *hx37=new TH1F("hx37","X cood, Si37", 100,0.,0.);
  //----------------END FACA----------------------------------------

  TH2F *h2x20v10=new TH2F("h2x20v10","X, x20:x10", 100, 0, L10, 100, 0, L20);
  TH2F *h2x30v20=new TH2F("h2x30v20","X, x30:x20", 100, 0, L20, 100, 0, L30);
  TH2F *h2x31v30=new TH2F("h2x31v30","X, x31:x30", 100, 0, L30, 100, 0, L31);
  
  //----------------START FACA----------------------------------------
  //TH2F *h2x33v31=new TH2F("h2x33v31","X, x33:x31", 100, 0, L31, 100, 0, L33);
  //TH2F *h2x35v33=new TH2F("h2x35v33","X, x35:x33", 100, 0, L33, 100, 0, L35);
  //TH2F *h2x37v35=new TH2F("h2x37v35","X, x37:x35", 100, 0, L35, 100, 0, L37);
  
   TH2F *h2x33v31=new TH2F("h2x33v31","X, x33:x31", 100, 0., 0., 100, 0., 0.);
  TH2F *h2x35v33=new TH2F("h2x35v33","X, x35:x33", 100, 0., 0., 100, 0., 0.);
  TH2F *h2x37v35=new TH2F("h2x37v35","X, x37:x35", 100, 0., 0., 100, 0., 0.);
  //----------------END FACA----------------------------------------

  TH2F *h2y20v10=new TH2F("h2y20v10","Y, y20:y10", 100, 0, L10, 100, 0, L20);
  TH2F *h2y30v20=new TH2F("h2y30v20","Y, y30:y20", 100, 0, L20, 100, 0, L30);
  TH2F *h2y31v30=new TH2F("h2y31v30","Y, y31:y30", 100, 0, L30, 100, 0, L31);
  
   //----------------START FACA----------------------------------------
   TH2F *h2y33v31=new TH2F("h2y33v31","Y, y33:y31", 100, 0, L31, 100, 0, L33);
   TH2F *h2y35v33=new TH2F("h2y35v33","Y, y35:y33", 100, 0, L33, 100, 0, L35);
   TH2F *h2y37v35=new TH2F("h2y37v35","Y, y37:y35", 100, 0, L35, 100, 0, L37);  
   //----------------END FACA----------------------------------------

 
  TH1F *h20diff = new TH1F("h20diff","X: x20-x10", 100, -1., 1.);
  TH1F *h30diff = new TH1F("h30diff","X: x30-x20", 100, -8., 8.);
  TH1F *h31diff = new TH1F("h31diff","X: x31-x30", 100, -8., 8.);
  
  //----------------START FACA----------------------------------------
  //TH1F *h33diff = new TH1F("h33diff","X: x33-x31", 100, -1., 1.);
  //TH1F *h35diff = new TH1F("h35diff","X: x35-x33", 100, -1., 1.);
  //TH1F *h37diff = new TH1F("h37diff","X: x37-x35", 100, -1., 1.);
  TH1F *h33diff = new TH1F("h33diff","X: x33-x31", 100, 0., 0.);
  TH1F *h35diff = new TH1F("h35diff","X: x35-x33", 100, 0., 0.);
  TH1F *h37diff = new TH1F("h37diff","X: x37-x35", 100, 0., 0.);
  //----------------END FACA----------------------------------------

  TH1F *h20xz = new TH1F("h20xz","X: x20-x10, mrad", 100, -3., 3.);
  TH1F *h20yz = new TH1F("h20yz","Y: y20-y10, mrad", 100, -3., 3.);
  
   //----------------START FACA----------------------------------------
    //TH1F *h33xz = new TH1F("h33xz","X: x33-x31, mrad", 100, -3., 3.);
    //TH1F *h35xz = new TH1F("h35xz","X: x35-x33, mrad", 100, -3., 3.);
    //TH1F *h37xz = new TH1F("h37xz","X: x37-x35, mrad", 100, -3., 3.);
    TH1F *h33xz = new TH1F("h33xz","X: x33-x31, mrad", 100, 0., 0.);
    TH1F *h35xz = new TH1F("h35xz","X: x35-x33, mrad", 100, 0., 0.);
    TH1F *h37xz = new TH1F("h37xz","X: x37-x35, mrad", 100, 0., 0.);
    
    TH1F *h33yz = new TH1F("h33yz","Y: y33-y31, mrad", 100, -3., 3.);
    TH1F *h35yz = new TH1F("h35yz","Y: y35-y33, mrad", 100, -3., 3.);
    TH1F *h37yz = new TH1F("h37yz","Y: y37-y35, mrad", 100, -3., 3.);
    //----------------END FACA----------------------------------------
    
  
  TH2F *h2x20v10all=new TH2F("h2x20v10all","X, x20:x10", 100, 0, L10, 100, 0, L20);
  TH2F *h2x30v20all=new TH2F("h2x30v20all","X, x30:x20", 100, 0, L20, 100, 0, L30);
  TH2F *h2x31v30all=new TH2F("h2x31v30all","X, x31:x30", 100, 0, L30, 100, 0, L31);
  
    //----------------START FACA----------------------------------------
     TH2F *h2x33v31all=new TH2F("h2x33v31all","X, x33:x31", 100, 0, L31, 100, 0, L33);
     TH2F *h2x35v33all=new TH2F("h2x35v33all","X, x35:x33", 100, 0, L33, 100, 0, L35);
     TH2F *h2x37v35all=new TH2F("h2x37v35all","X, x37:x35", 100, 0, L35, 100, 0, L37);
      //----------------END FACA----------------------------------------
  
  

  TH2F *h2y20v10all=new TH2F("h2y20v10all","Y, y20:y10", 100, 0, L10, 100, 0, L20);
  TH2F *h2y30v20all=new TH2F("h2y30v20all","Y, y30:y20", 100, 0, L20, 100, 0, L30);
  TH2F *h2y31v30all=new TH2F("h2y31v30all","Y, y31:y30", 100, 0, L30, 100, 0, L31);
  
  //----------------START FACA----------------------------------------
  TH2F *h2y33v31all=new TH2F("h2y33v31all","Y, y33:y31", 100, 0, L31, 100, 0, L33);
  TH2F *h2y35v33all=new TH2F("h2y35v33all","Y, y35:y33", 100, 0, L33, 100, 0, L35);
  TH2F *h2y37v35all=new TH2F("h2y37v35all","Y, y37:y35", 100, 0, L35, 100, 0, L37);
  
    //----------------END FACA----------------------------------------

//   Float_t L10(1.91), L20(1.91), L30(9.50), L31(9.50);
//   Float_t L32(8.0),L33(8.0),L34(9.5),L35(9.5),L36(18.0),L37(18.0);
   Float_t dX10(0.),dX20(0.),dX30(0.0),dX31(0.0);
   Float_t dY10(0.),dY20(0.),dY30(0.0),dY31(0.0);
   Float_t dX32(0.),dX33(0.),dX34(0.),dX35(0.),dX36(0.),dX37(0.);
   Float_t dY32(0.),dY33(0.),dY34(0.),dY35(0.),dY36(0.),dY37(0.);

   Float_t D32(8.0),D33(8.0),D34(14.0),D35(14.0),D36(24.0),D37(24.0);

   if( Lab == 1) {
       dX10 = -L10/2.0;
       dX20 = -L20/2.0 - 0.2274;
       dX30 = -L30/2.0 + 0.4229;
       dX31 = -L31/2.0 + 2.0740;

       dY10 = -L10/2.0;
       dY20 = -L20/2.0 - 0.3533;
       dY30 = -L30/2.0 - 0.5034;
       dY31 = -L31/2.0 - 0.7914;

       dX32 = D32 - L32/2.0 + 0.0;
       dX34 = D34 - L34/2.0 + 0.0;
       dX36 = D36 - L36/2.0 + 0.0;
        
       dX33 = -D33 - L33/2.0 + 0.0;
       dX35 = -D35 - L35/2.0 + 0.0;
       dX37 = -D37 - L37/2.0 + 0.0;
       
       dY32 = -L32/2.0 +0.0;
       dY34 = -L34/2.0 +0.0;
       dY36 = -L36/2.0 +0.0;

       dY33 = -L33/2.0 +0.0;
       dY35 = -L35/2.0 +0.0;
       dY37 = -L37/2.0 +0.0;
   }


   /*
   newtree->Branch("hnx10","TH1F",&hnx10,32000,0);
   newtree->Branch("hnx20","TH1F",&hnx20,32000,0);
   newtree->Branch("hnx30","TH1F",&hnx30,32000,0);
   newtree->Branch("hnx31","TH1F",&hnx31,32000,0);

   newtree->Branch("hx10","TH1F",&hx10,32000,0);
   newtree->Branch("hx20","TH1F",&hx20,32000,0);
   newtree->Branch("hx30","TH1F",&hx30,32000,0);
   newtree->Branch("hx31","TH1F",&hx31,32000,0);

   newtree->Branch("h2x20v10","TH2F",&h2x20v10,32000,0);
   newtree->Branch("h2x30v20","TH2F",&h2x30v20,32000,0);
   newtree->Branch("h2x31v30","TH2F",&h2x31v30,32000,0);

   newtree->Branch("h20diff","TH1F",&h20diff,32000,0);
   newtree->Branch("h30diff","TH1F",&h30diff,32000,0);
   newtree->Branch("h31iff","TH1F",&h31diff,32000,0);
   */

   // Si: before MAG
   Int_t nx10,nx20,nx30,nx31;
   Int_t ny10,ny20,ny30,ny31;
   Float_t x10[10],x20[10],x30[10],x31[10];
   Float_t y10[10],y20[10],y30[10],y31[10];
   Float_t z10[10],z20[10],z30[10],z31[10];


   newtree->Branch("nx10",&nx10,"nx10/I");
   newtree->Branch("nx20",&nx20,"nx20/I");
   newtree->Branch("nx30",&nx30,"nx30/I");
   newtree->Branch("nx31",&nx31,"nx31/I");

   newtree->Branch("x10",&x10,"x10[10]/F");
   newtree->Branch("x20",&x20,"x20[10]/F");
   newtree->Branch("x30",&x30,"x30[10]/F");
   newtree->Branch("x31",&x31,"x31[10]/F");

   newtree->Branch("ny10",&ny10,"ny10/I");
   newtree->Branch("ny20",&ny20,"ny20/I");
   newtree->Branch("ny30",&ny30,"ny30/I");
   newtree->Branch("ny31",&ny31,"ny31/I");

   newtree->Branch("y10",&y10,"y10[10]/F");
   newtree->Branch("y20",&y20,"y20[10]/F");
   newtree->Branch("y30",&y30,"y30[10]/F");
   newtree->Branch("y31",&y31,"y31[10]/F");

   newtree->Branch("z10",&z10,"z10[10]/F");
   newtree->Branch("z20",&z20,"z20[10]/F");
   newtree->Branch("z30",&z30,"z30[10]/F");
   newtree->Branch("z31",&z31,"z31[10]/F");


   // arm: e-
   Int_t nx32, nx34, nx36;
   Int_t ny32, ny34, ny36;
   Float_t x32[10],x34[10],x36[10];
   Float_t y32[10],y34[10],y36[10];
   Float_t z32[10],z34[10],z36[10];

   newtree->Branch("nx32",&nx32,"nx32/I");
   newtree->Branch("nx34",&nx34,"nx34/I");
   newtree->Branch("nx36",&nx36,"nx36/I");
 
   newtree->Branch("x32",&x32,"x32[10]/F");
   newtree->Branch("x34",&x34,"x34[10]/F");
   newtree->Branch("x36",&x36,"x36[10]/F");
 
   newtree->Branch("ny32",&ny32,"ny32/I");
   newtree->Branch("ny34",&ny34,"ny34/I");
   newtree->Branch("ny36",&ny36,"ny36/I");

   newtree->Branch("y32",&y32,"y32[10]/F");
   newtree->Branch("y34",&y34,"y34[10]/F");
   newtree->Branch("y36",&y36,"y36[10]/F");
 
   newtree->Branch("z32",&z32,"z32[10]/F");
   newtree->Branch("z34",&z34,"z34[10]/F");
   newtree->Branch("z36",&z36,"z36[10]/F");
 
   // arm: e+
   Int_t nx33, nx35, nx37;
   Int_t ny33, ny35, ny37;
   Float_t x33[10],x35[10],x37[10];
   Float_t y33[10],y35[10],y37[10];
   Float_t z33[10],z35[10],z37[10];

   newtree->Branch("nx33",&nx33,"nx33/I");
   newtree->Branch("nx35",&nx35,"nx35/I");
   newtree->Branch("nx37",&nx37,"nx37/I");
 
   newtree->Branch("x33",&x33,"x33[10]/F");
   newtree->Branch("x35",&x35,"x35[10]/F");
   newtree->Branch("x37",&x37,"x37[10]/F");
 
   newtree->Branch("ny33",&ny33,"ny33/I");
   newtree->Branch("ny35",&ny35,"ny35/I");
   newtree->Branch("ny37",&ny37,"ny37/I");

   newtree->Branch("y33",&y33,"y33[10]/F");
   newtree->Branch("y35",&y35,"y35[10]/F");
   newtree->Branch("y37",&y37,"y37[10]/F");
   
   newtree->Branch("z33",&z33,"z33[10]/F");
   newtree->Branch("z35",&z35,"z35[10]/F");
   newtree->Branch("z37",&z37,"z37[10]/F");
   // ------------------//
   // Loop over hits    //
   // ------------------//
   
   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   
 //  Int_t nn10(0),nn20(0),nn30(0),nn31(0),nn32(0),nn33(0),nn34(0),nn35(0),nn36(0),nn37(0);
   for(Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if(ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

      // Temporary hit index
      
      // zero for nhit
      nx10=0;nx20=0;nx30=0;nx31=0;nx32=0;nx33=0;nx34=0;nx35=0;nx36=0;nx37=0;
      ny10=0;ny20=0;ny30=0;ny31=0;ny32=0;ny33=0;ny34=0;ny35=0;ny36=0;ny37=0;

 //     Int_t nn10(0),nn20(0),nn30(0),nn31(0),nn32(0),nn33(0),nn34(0),nn35(0),nn36(0),nn37(0);

      for(Int_t i=0; i<10; i++){
        x10[i]=-999.;x20[i]=-999.;x30[i]=-999.;x31[i]=-999.;x32[i]=-999.;x33[i]=-999.;x34[i]=-999.;x35[i]=-999.;x36[i]=-999.;x37[i]=-999.;
        y10[i]=-999.;y20[i]=-999.;y30[i]=-999.;y31[i]=-999.;y32[i]=-999.;y33[i]=-999.;y34[i]=-999.;y35[i]=-999.;y36[i]=-999.;y37[i]=-999.;
      }

      for(Int_t jhit=0; jhit<nhits; jhit++) {
          if(xh[jhit]>-999){
            if(subdet[jhit]==10) {
               x10[nx10]=xh[jhit] + dX10;
               z10[nx10]=zh[jhit];
               nx10++;
            }//if(subdet[jhit]) 
           
            if(subdet[jhit]==20) {
               x20[nx20]=xh[jhit] + dX20;
               z20[nx20]=zh[jhit];
               nx20++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==30) {
               x30[nx30]=xh[jhit] + dX30;
               z30[nx30]=zh[jhit];
               nx30++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==31) {
               x31[nx31]=xh[jhit] + dX31;
               z31[nx31]=zh[jhit];
               //cout << jhit << " " <<  subdet[jhit] << " " << xh[jhit] << " ----- " << nn31 << " " << x31[nn31]<< endl;
               nx31++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==32) {
               x32[nx32]=xh[jhit] + dX32;
               z32[nx32]=zh[jhit];
               nx32++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==33) {
               //cout << subdet[jhit] << " " << xh[jhit] << endl;
               x33[nx33]=xh[jhit] + dX33;
               z33[nx33]=zh[jhit]; 
               nx33++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==34) {
               x34[nx34]=xh[jhit] + dX34;
               z34[nx34]=zh[jhit];
               nx34++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==35) {
               x35[nx35]=xh[jhit] + dX35;
               z35[nx35]=zh[jhit];
               nx35++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==36) {
               x36[nx36]=xh[jhit] + dX36;
               z36[nx36]=zh[jhit]; 
               nx36++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==37) {
               x37[nx37]=xh[jhit] + dX37;
               z37[nx37]=zh[jhit];
               nx37++;
            }//if(subdet[jhit]) 
          }// if(xh[

        if(yh[jhit]>-999){
            if(subdet[jhit]==10) {
               y10[ny10]=yh[jhit] + dY10; ny10++;
            }//if(subdet[jhit]) 
           
            if(subdet[jhit]==20) {
               y20[ny20]=yh[jhit] + dY20; ny20++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==30) {
               y30[ny30]=yh[jhit] + dY30; ny30++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==31) {
               y31[ny31]=yh[jhit] + dY31; ny31++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==32) {
               y32[ny32]=yh[jhit] + dY32; ny32++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==33) {
               y33[ny33]=yh[jhit] + dY33; ny33++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==34) {
               y34[ny34]=yh[jhit] + dY34; ny34++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==35) {
               y35[ny35]=yh[jhit] + dY35; ny35++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==36) {
               y36[ny36]=yh[jhit] + dY36; ny36++;
            }//if(subdet[jhit]) 

            if(subdet[jhit]==37) {
               y37[ny37]=yh[jhit] + dY37; ny37++;
            }//if(subdet[jhit]) 
          }// if(yh[
      } // loop on hits 

      if(nx31>9 || nx33>9 || nx35>9 || nx37>9)  cout << nx31 << " " << nx33 << " " << nx35 << " " << nx37 << endl;
      if(nx32>9 || nx34>9 || nx36>9)  cout << nx32 << " " << nx34 << " " << nx36 << endl;


      // -----------------------------------
      // Fill Mon Histograms for calib runs
      // -----------------------------------
      hnx10->Fill(nx10);
      hnx20->Fill(nx20);
      hnx30->Fill(nx30);
      hnx31->Fill(nx31);
      
      //----------------START FACA----------------------------------------
       hnx33->Fill(nx33);
       hnx35->Fill(nx35);
       hnx37->Fill(nx37);
      //----------------END FACA----------------------------------------

      if(nx10==1) hx10->Fill(x10[0]);
      if(nx20==1) hx20->Fill(x20[0]);
      if(nx30==1) hx30->Fill(x30[0]);
      if(nx31==1) hx31->Fill(x31[0]);
      
       //----------------START FACA----------------------------------------
       if(nx33==1) hx33->Fill(x33[0]);
       if(nx35==1) hx35->Fill(x35[0]);
       if(nx37==1) hx37->Fill(x37[0]);
       
         //----------------END FACA----------------------------------------
 
      if(nx10==1 && nx20==1) h2x20v10->Fill(x10[0],x20[0]);
      if(nx20==1 && nx30==1) h2x30v20->Fill(x20[0],x30[0]);
      if(nx30==1 && nx31==1) h2x31v30->Fill(x30[0],x31[0]);

 //----------------START FACA----------------------------------------
   if(nx31==1 && nx33==1) h2x33v31->Fill(x31[0],x33[0]);
   if(nx33==1 && nx35==1) h2x35v33->Fill(x33[0],x35[0]);
   if(nx35==1 && nx37==1) h2x37v35->Fill(x35[0],x37[0]);
         //----------------END FACA----------------------------------------
 
      for(Int_t j=0; j<10; j++){
         h2x20v10all->Fill(x10[j],x20[j]);
         h2x30v20all->Fill(x20[j],x30[j]);
         h2x31v30all->Fill(x30[j],x31[j]);
         //----------------START FACA----------------------------------------
             h2x33v31all->Fill(x31[j],x33[j]);
             h2x35v33all->Fill(x33[j],x35[j]);
             h2x37v35all->Fill(x35[j],x37[j]);
            //----------------END FACA----------------------------------------
      }      

      if(ny10==1 && ny20==1) h2y20v10->Fill(y10[0],y20[0]);
      if(ny20==1 && ny30==1) h2y30v20->Fill(y20[0],y30[0]);
      if(ny30==1 && ny31==1) h2y31v30->Fill(y30[0],y31[0]);
       //----------------START FACA----------------------------------------
        if(ny31==1 && ny33==1) h2y33v31->Fill(y31[0],y33[0]);
        if(ny33==1 && ny35==1) h2y35v33->Fill(y33[0],y35[0]);
        if(ny35==1 && ny37==1) h2y37v35->Fill(y35[0],y37[0]);
           //----------------END FACA----------------------------------------
       
      for(Int_t j=0; j<10; j++){
         h2y20v10all->Fill(y10[j],y20[j]);
         h2y30v20all->Fill(y20[j],y30[j]);
         h2y31v30all->Fill(y30[j],y31[j]);
         
          //----------------START FACA----------------------------------------
           h2y33v31all->Fill(y31[j],y33[j]);
           h2y35v33all->Fill(y33[j],y35[j]);
           h2y37v35all->Fill(y35[j],y37[j]);
          
              //----------------END FACA----------------------------------------
          
      }      

      if(nx10==1 && nx20==1) h20diff->Fill(x20[0]-x10[0]);
      if(nx20==1 && nx30==1) h30diff->Fill(x30[0]-x20[0]);
      if(nx30==1 && nx31==1) h31diff->Fill(x31[0]-x30[0]);
      
      //----------------START FACA----------------------------------------
      if(nx31==1 && nx33==1) h33diff->Fill(x33[0]-x31[0]);
      if(nx33==1 && nx35==1) h35diff->Fill(x35[0]-x33[0]);
	  if(nx35==1 && nx37==1) h37diff->Fill(x37[0]-x35[0]);
      
          //----------------END FACA----------------------------------------

      if(nx10==1 && nx20==1) h20xz->Fill(1e3*(x20[0]-x10[0])/(z20[0]-z10[0]));
      if(ny10==1 && ny20==1) h20yz->Fill(1e3*(y20[0]-y10[0])/(z20[0]-z10[0]));
      
      //----------------START FACA----------------------------------------
       if(nx31==1 && nx33==1) h33xz->Fill(1e3*(x33[0]-x31[0])/(z33[0]-z31[0]));
        if(nx33==1 && nx35==1) h35xz->Fill(1e3*(x35[0]-x33[0])/(z35[0]-z33[0]));
         if(nx35==1 && nx37==1) h37xz->Fill(1e3*(x37[0]-x35[0])/(z37[0]-z35[0]));
         
       if(ny31==1 && ny33==1) h33yz->Fill(1e3*(y33[0]-y31[0])/(z33[0]-z31[0]));
        if(ny33==1 && ny35==1) h35yz->Fill(1e3*(y35[0]-y33[0])/(z35[0]-z33[0]));
         if(ny35==1 && ny37==1) h37yz->Fill(1e3*(y37[0]-y35[0])/(z37[0]-z35[0]));
    
	//----------------END FACA----------------------------------------
    
	  newtree->Fill();
   } // loop in jentry
   
    //----------------START FACA----------------------------------------
   
   TCanvas *c0 = new TCanvas("c0","Si33",1280,1024);
			gStyle->SetOptFit();
			gStyle->SetOptStat(111110);
			hx33->Draw();
 			TF1 *g133 = new TF1 ("m133", "gaus", -12, -7);
  			g133->SetLineColor(kYellow);
  			TF1 *g233 = new TF1 ("m233", "gaus", -7, -3);
  			g233->SetLineColor(kGreen);
  			TF1 *f133 = new TF1("double_gaus", "gaus(0) + gaus(3)", -12, -3);
  			f133->SetParNames("Constant 1", "Mean 1", "Sigma 1","Constant 2", "Mean 2", "Sigma 2");
  			f133->SetLineColor(kRed);  
			hx33->Fit(g133, "R");
  			hx33->Fit(g233, "R");
  			Double_t par33[6];
  			g133->GetParameters(&par33[0]);
  			g233->GetParameters(&par33[3]);
  			f133->SetParameters(par33);
  			hx33->Fit(f133, "R");
  		    hx33->Draw("e1");
  			g133->Draw("SAME");
  			g233->Draw("SAME");
  			f133->Draw("SAME");
  			if (p1p233==1) {
  				p133=f133->GetParameter(1);
  				p233=f133->GetParameter(2);	
  				x[0]=p133;
			  }
			else  if (p1p233==2) {
  				p133=f133->GetParameter(4);
  				p233=f133->GetParameter(5);	
  				x[0]=p133;
			  }
			hx33->SetTitle("si-500596, #mu^{+}_{22 GeV}- x coordinates Silicon Tracker 33 @z=1873.93 cm"); 
			hx33->GetXaxis()->SetTitle("x (cm)"); 
			hx33->GetYaxis()->SetTitle("#Counts");
			hx33->SetFillColorAlpha(kBlue,0.35);
	   		hx33->Draw();
	   		hx33->SetName("Si 33");
	   		gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			c0->Update();
   			TLegend* leg33 = new TLegend(0.15, 0.7, .25, .75);
   			leg33->SetHeader("Legend");
  			leg33->SetNColumns(1);
   			leg33->AddEntry(hx33, "Data", "l");
 			leg33->AddEntry(f133, "Fit", "l");  //TO USE FOR 2 GAUSSIAN FIT
   			leg33->Draw(); 
   			c0->Update();
   			gPad->Modified();
			gPad->Update();
			TPaveStats *stat33= (TPaveStats*)hx33->GetListOfFunctions()->FindObject("stats");
			stat33->SetTextColor(kBlue);
			stat33->SetX1NDC(0.80); stat33->SetX2NDC(0.98);
			stat33->SetY1NDC(0.77); stat33->SetY2NDC(0.92);
			stat33->AddText(TString::Format("Mean = %g", p133));
			stat33->AddText(TString::Format("Sigma	= %g", p233));
			stat33->DrawClone();
			gPad->Update();
			c0->Print("C:/si-calo-sep18/Si33.pdf");
			delete c0;
			
			TCanvas *c1 = new TCanvas("c1","Si35",1280,1024);
			gStyle->SetOptFit();
			gStyle->SetOptStat(111110);
			hx35->Draw();
 			TF1 *g135 = new TF1 ("m135", "gaus", -19, -13);
  			g135->SetLineColor(kYellow);
  			TF1 *g235 = new TF1 ("m235", "gaus", -14, -10);
  			g235->SetLineColor(kGreen);
  			TF1 *f135 = new TF1("double_gaus", "gaus(0) + gaus(3)", -19, -10);
  			f135->SetParNames("Constant 1", "Mean 1", "Sigma 1","Constant 2", "Mean 2", "Sigma 2");
  			f135->SetLineColor(kRed);  
			hx35->Fit(g135, "R");
  			hx35->Fit(g235, "R");
  			Double_t par35[6];
  			g135->GetParameters(&par35[0]);
  			g235->GetParameters(&par35[3]);
  			f135->SetParameters(par35);
  			hx35->Fit(f135, "R");
  		    hx35->Draw("e1");
  			g135->Draw("SAME");
  			g235->Draw("SAME");
  			f135->Draw("SAME");
  			if (p1p235==1) {
  				p135=f135->GetParameter(1);
  				p235=f135->GetParameter(2);	
  				x[1]=p135;
			  }
			else  if (p1p235==2) {
  				p135=f135->GetParameter(4);
  				p235=f135->GetParameter(5);	
  				x[1]=p135;
			  }
			hx35->SetTitle("si-500596, #mu^{+}_{22 GeV}- x coordinates Silicon Tracker 35 @z=1995.43 cm"); 
			hx35->GetXaxis()->SetTitle("x (cm)"); 
			hx35->GetYaxis()->SetTitle("#Counts");
			hx35->SetFillColorAlpha(kBlue,0.35);
	   		hx35->Draw();
	   		hx35->SetName("Si 35");
	   		gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			c1->Update();
   			TLegend* leg35 = new TLegend(0.15, 0.7, .25, .75);
   			leg35->SetHeader("Legend");
  			leg35->SetNColumns(1);
   			leg35->AddEntry(hx35, "Data", "l");
 			leg35->AddEntry(f135, "Fit", "l");  //TO USE FOR 2 GAUSSIAN FIT
   			leg35->Draw(); 
   			c1->Update();
   			gPad->Modified();
			gPad->Update();
			TPaveStats *stat35= (TPaveStats*)hx35->GetListOfFunctions()->FindObject("stats");
			stat35->SetTextColor(kBlue);
			stat35->SetX1NDC(0.80); stat35->SetX2NDC(0.98);
			stat35->SetY1NDC(0.77); stat35->SetY2NDC(0.92);
			stat35->AddText(TString::Format("Mean = %g", p135));
			stat35->AddText(TString::Format("Sigma	= %g", p235));
			stat35->DrawClone();
			gPad->Update();
			c1->Print("C:/si-calo-sep18/Si35.pdf");
			delete c1;
			
			TCanvas *c2 = new TCanvas("c2","Si37",1280,1024);
			gStyle->SetOptFit();
			gStyle->SetOptStat(111110);
			hx37->Draw();
 			TF1 *g137 = new TF1 ("m137", "gaus", -28, -24);
  			g137->SetLineColor(kYellow);
  			TF1 *g237 = new TF1 ("m237", "gaus", -24, -16);
  			g237->SetLineColor(kGreen);
  			TF1 *f137 = new TF1("double_gaus", "gaus(0) + gaus(3)", -28, -16);
  			f137->SetParNames("Constant 1", "Mean 1", "Sigma 1","Constant 2", "Mean 2", "Sigma 2");
  			f137->SetLineColor(kRed);  
			hx37->Fit(g137, "R");
  			hx37->Fit(g237, "R");
  			Double_t par37[6];
  			g137->GetParameters(&par37[0]);
  			g237->GetParameters(&par37[3]);
  			f137->SetParameters(par37);
  			hx37->Fit(f137, "R");
  		    hx37->Draw("e1");
  			g137->Draw("SAME");
  			g237->Draw("SAME");
  			f137->Draw("SAME");
  			if (p1p237==1) {
  				p137=f137->GetParameter(1);
  				p237=f137->GetParameter(2);	
  				x[2]=p137;
			  }
			else  if (p1p237==2) {
  				p137=f137->GetParameter(4);
  				p237=f137->GetParameter(5);	
  				x[2]=p137;
			  }
			hx37->SetTitle("si-500596, #mu^{+}_{22 GeV}- x coordinates Silicon Tracker 37 @z=2125.43 cm"); 
			hx37->GetXaxis()->SetTitle("x (cm)"); 
			hx37->GetYaxis()->SetTitle("#Counts");
			hx37->SetFillColorAlpha(kBlue,0.35);
	   		hx37->Draw();
	   		hx37->SetName("Si 37");
	   		gPad->Modified();
			gPad->Update(); // make sure it's really (re)drawn
   			c2->Update();
   			TLegend* leg37 = new TLegend(0.15, 0.7, .25, .75);
   			leg37->SetHeader("Legend");
  			leg37->SetNColumns(1);
   			leg37->AddEntry(hx37, "Data", "l");
 			leg37->AddEntry(f137, "Fit", "l");  //TO USE FOR 2 GAUSSIAN FIT
   			leg37->Draw(); 
   			c2->Update();
   			gPad->Modified();
			gPad->Update();
			TPaveStats *stat37= (TPaveStats*)hx37->GetListOfFunctions()->FindObject("stats");
			stat37->SetTextColor(kBlue);
			stat37->SetX1NDC(0.80); stat37->SetX2NDC(0.98);
			stat37->SetY1NDC(0.77); stat37->SetY2NDC(0.92);
			stat37->AddText(TString::Format("Mean = %g", p137));
			stat37->AddText(TString::Format("Sigma	= %g", p237));
			stat37->DrawClone();
			gPad->Update();
			c2->Print("C:/si-calo-sep18/Si37.pdf");
			delete c2;
			
		TCanvas *c3 = new TCanvas("c3", "track", 200, 10, 700, 500);
  		TGraphErrors *g = new TGraphErrors(3, z, x, 0, ex); // hits x versus z
  		g->SetTitle("si-500596, #mu^{+}_{22 GeV}- x-z coordinates Silicon Trackers 33-35-37;z (cm);x (cm)");
  		g->SetMarkerColor(4);
  		g->SetMarkerStyle(21);
  		const char *f_name = "pol1"; // standard ROOT linear function
  		g->Fit(f_name, ""); // e.g. "" or "Q" or "V"
		TF1 *f = g->GetFunction(f_name); // retrieve the fit function
  		if (f) {
    			std::cout << "fit function name = " << f->GetName() << std::endl;
    			std::cout << "Chi2 = " << f->GetChisquare() << std::endl;
    			std::cout << "NDF = " << f->GetNDF() << std::endl;
    			std::cout << "Prob = " << f->GetProb() << std::endl;
			    for (Int_t i = 0; i < f->GetNpar(); i++) {
      				std::cout << f->GetParName(i) << " = " << f->GetParameter(i)
                	<< " +- " << f->GetParError(i) << std::endl;
    				}
  				}
 		g->Draw("ap");
 		gPad->Modified();
		gPad->Update(); // make sure it's really (re)drawn
   		c3->Update();
   		TLegend* leg = new TLegend(0.65, 0.7, .75, .75);
   		leg->SetHeader("Legend");
  		leg->SetNColumns(1);
   		leg->AddEntry(g, "Data", "l");
 		leg->AddEntry(f_name, "Fit", "l");  //TO USE FOR 2 GAUSSIAN FIT
   		leg->Draw(); 
   		c3->Update();
   		gPad->Modified();
		gPad->Update();
 		TPaveStats *stat= (TPaveStats*)g->GetListOfFunctions()->FindObject("stats");
		stat->SetTextColor(kBlue);
		stat->SetX1NDC(0.80); stat->SetX2NDC(0.98);
		stat->SetY1NDC(0.77); stat->SetY2NDC(0.92);
		gPad->Update();
  		c3->SaveAs("C:/si-calo-sep18/Fit_x-z_33-35-37.pdf");
  // delete c; // no longer needed
  // delete g; // no longer needed
			
		
    //----------------END FACA----------------------------------------
   
   newtree->Fill();
   newfile->Write();
  /* h20xz->Fit("gaus");
   h20yz->Fit("gaus");
   //----------------START FACA----------------------------------------
   h33xz->Fit("gaus");
   h33yz->Fit("gaus");
   h35xz->Fit("gaus");
   h35yz->Fit("gaus");
   h37xz->Fit("gaus");
   h37yz->Fit("gaus");*/
   	//----------------END FACA----------------------------------------
      
       
       
//   delete newfile;

}