#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>
#include <fstream>
//
// 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 a,b,ae,be,chi2,q;
   float x[3];					  // x coordinates of hits 
   float z[3]; 					 // z coordinates of hits 
   float zLG[1]; 					 // z coordinates of hits 
   float ex[3]; // errors of x coordinates of hits
   float p101;
   float p102;
   float p201;
   float p202;
   int p1p201=2;
   int p1p202=1;
   //----------------END FACA----------------------------------------
  
  
void TBAlign::Loop()
{
	//----------------START FACA----------------------------------------
	ofstream results;
 	results.open("C:/si-calo-sep18/Fit_parameters.txt",  ios::out);
 	char s[]="Fit Parameters";
	int i;
	int n=3;
 //----------------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----------------------------------------
    
    
    //----------------START FACA----------------------------------------
     //TH1F *hxb0n12=new TH1F("hxb0n12","b parameter LG 0 !1 !2", 100,0.,0.);
     //TH1F *hxb2n1=new TH1F("hxb2n1","b parameter LG 2 !1", 100,0.,0.);
     
    /*  TH1F *hbT=new TH1F("hbT","b parameter 33-35-37 xh >-99", 100,0.,0.);
       TH1F *hbTx37=new TH1F("hbTx37","b parameter 33-35-37 xh >-99 x37 ", 100,0.,0.);
       TH1F *hbTx37LG0=new TH1F("hbTx37LG0","b parameter 33-35-37 xh >-99 x37 LG0", 100,0.,0.);
       TH1F *hbTx37H=new TH1F("hbTx37H","b parameter 33-35-37 xh >-99 x37 HORSA", 100,0.,0.);*/
     //----------------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];
    //----------------START FACA----------------------------------------
    Float_t bpar[10];
    Float_t bparT[10];
    Float_t xLG[10];
     //----------------END FACA----------------------------------------
    
   

   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");
   
     //----------------START FACA----------------------------------------
      newtree->Branch("bpar",&bpar,"bpar");
      newtree->Branch("bparT",&bparT,"bparT");
      newtree->Branch("xLG",&xLG,"xLG");
       //----------------END FACA----------------------------------------
       
       //----------------START FACA----------------------------------------
       
       			TCut cut912 = TString::Format("TMath::Abs(Calo_Time[9] - Calo_Time[12] - (3.32785)) < (3)*(1.47188) ").Data();
  				TCut cut913 = TString::Format("TMath::Abs(Calo_Time[9] - Calo_Time[13] - (4.43584)) < (3)*(0.685633) ").Data();
				TCut cut914 = TString::Format("TMath::Abs(Calo_Time[9] - Calo_Time[14] - (3.42323)) < (3)*(0.646554) ").Data();
				TCut cut915 = TString::Format("TMath::Abs(Calo_Time[9] - Calo_Time[15] - (1.71963)) < (3)*(1.78036) ").Data();
					TCut cut90 = TString::Format("TMath::Abs(Calo_Time[9] - Calo_Time[0] - (-125.007)) < (3)*(2.67618) ").Data();
					TCut xh99 = TString::Format("xh > -99 ").Data();
					TCut cutx37 = TString::Format("TMath::Abs(x37[0] - (-20.7411)) < (1)*(1.16973) ").Data();
					TCut subdet33 = TString::Format("subdet == 33 ").Data();
					TCut subdet35 = TString::Format("subdet == 35 ").Data();
					TCut subdet37 = TString::Format("subdet == 37 ").Data();
       //----------------END FACA----------------------------------------
   
   
   // ------------------//
   // 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();
	//----------------START FACA----------------------------------------
    
	  
	  			if(nx33==1 && nx35==1 && nx37==1){
	  			if(x[0]>-999 && x[1]>-999 && x[2]>-999){
			 	x[0]=x33[0];
        		x[1]=x35[0];
        		x[2]=x37[0];
				z[0]=1873.93;
				z[1]=1995.43;
				z[2]=2125.43;
				ex[0]=0.1;
				ex[1]=0.1;
				ex[2]=0.1;
				zLG[0]=2195.43; 
					fit_line(z, x, n, ex, a, b, ae, be, chi2, q); // hits x versus z
					xLG[0]=a+b*zLG[0];					 // x LG coordinates of hits 
					std::cout << "x0 = " << x[0] << "\t x1 = " << x[1]  << " \t x2 = " << x[2] << "\t xLG = " << xLG[0] << "\t z0 = " << z[0] << "\t z1 = " << z[1]  << "\t z2 = " << z[2] << "\t zLG = " << zLG[0] << std::endl;
					results << jentry << "\t " <<  a << "\t " << ae << "\t " << b << "\t " << be <<  "\t " << chi2 << "\t "<< q << "\t "<< xLG[0] << endl;
					bpar[0]=b;
					//hbT->Fill(bpar[0]);
					if (x33[0]!=NULL && x35[0]!=NULL && x37[0]!=NULL) bparT[0]=b;
				/*	if (subdet37 && xh99) hbT->Fill(bpar[0]);
					if (subdet37 && xh99 && cutx37 ) hbTx37->Fill(bpar[0]);
					if (subdet37 && xh99 && cut90 ) hbTx37LG0->Fill(bpar[0]);
					if (subdet37 && xh99 && cut912 && cut913 && cut914 && cut915 ) hbTx37H->Fill(bpar[0]);*/
					
						}
				}	
		//----------------END FACA----------------------------------------

   } // loop in jentry
   results.close(); //close the file
   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;

}