#include <stdlib.h>
#include <iostream>
#include <iomanip>
//
#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMatrixD.h>
#include <TMatrixF.h>
#include <TArrayF.h>
#include <TArrayI.h>
//
#include <PamLevel2.h>
#include <CaloFranzini.h>
//
using namespace std;
//
extern Bool_t MATRIX;
extern Bool_t FULL;
extern Bool_t SIMU;
extern Bool_t CRIG;
extern Bool_t SRIG;
CaloFranzini *cf;
Int_t nbin;
Float_t rig[18];
Float_t rmean[17];
Int_t ntot[17];
Int_t MDIM = 8213;
//Int_t MDIM = 4128;
//Float_t qqplane[17][43];
//Int_t nnqplane[17][43];

TArrayF *qplane[17];
TArrayI *nqplane[17];
TMatrixD *matrix[17];
TMatrixD *nmat[17];

//TMatrixD *fqplane;
//TMatrixD *fnqplane;
TMatrixD *fqplane[17];
TMatrixD *fnqplane[17];
TMatrixD *fmatrix[17];
//TMatrixF *fnmat[17];
//TMatrixD *fmatrix;
//TMatrixD *fnmat;
//TMatrixF *fmatrix;
//TMatrixF *fnmat;
TMatrixF *fnmat[17];
//Int_t finmat[43][191];

//===============================================================================
bool Select( PamLevel2* event ){

  //---------------------------------------------------------
  // single track
  //---------------------------------------------------------
  if( event->GetTrkLevel2()->GetNTracks()!=1 ) return false; 
  PamTrack *track = event->GetTrack(0);
  if(!track)return false;

  if ( SIMU ) return true;
  //------------------------------------------------------------------
  // tracker pre-selection
  //------------------------------------------------------------------
  TrkTrack *trk = track->GetTrkTrack();
  float rigidity = trk->GetRigidity();
  if ( CRIG ) rigidity = event->GetCaloLevel2()->qtot/260.;
  if ( SRIG ) rigidity = event->GetGPamela()->P0;
  bool TRACK__OK = false;
  if( 
     trk->chi2 >0     &&
     trk->GetNX()>=4  &&
     trk->GetNY()>=3  &&
     trk->GetLeverArmX()>=5 &&
     true ) TRACK__OK = true;

  if( !TRACK__OK )return false;

  //------------------------------------------------------------------
  // TOF pre-selection
  //------------------------------------------------------------------
  bool TOF__OK = false;
  if(
     event->GetToFLevel2()->GetNHitPaddles(0) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(1) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(2) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(3) == 1 &&
     event->GetToFLevel2()->GetNHitPaddles(4) >= 1 &&
     event->GetToFLevel2()->GetNHitPaddles(5) >= 1 &&
     event->GetToFLevel2()->npmt() <= 18           &&
     !event->GetAcLevel2()->CARDhit() &&
     !event->GetAcLevel2()->CAThit() &&
     true ) TOF__OK = true;
  if( !TOF__OK && !SIMU)return false;
  //------------------------------------------------------
  // no albedo
  //------------------------------------------------------
  if(  	 !SIMU && (track->GetToFTrack()->beta[12]<=0.2 ||
		   track->GetToFTrack()->beta[12] >= 1.5) ) return false;

  //------------------------------------------------------
  bool CUT1 = false;
  if( 
     trk->nstep<100   &&
     rigidity<400.   &&
     rigidity>0.1   &&
     trk->resx[0]<0.001 && 
     trk->resx[5]<0.001 && 
     track->IsSolved() &&
     trk->IsInsideCavity() &&
     true ) CUT1 = true;
  //------------------------------------------------------
  if( !CUT1 )return false;
  if ( trk->GetDeflection()>0. && !SIMU ) return false;

  //
  //  ELENA'S CUT
  //
  //
  // lever-arm 6
  //====================================================
  bool LX6=false;
  if(
     track->GetTrkTrack()->GetLeverArmX()==6 &&
     !track->GetTrkTrack()->IsBad(0,0) &&
     !track->GetTrkTrack()->IsBad(5,0) &&
     track->GetTrkTrack()->resx[0]<0.001 &&
     track->GetTrkTrack()->resx[5]<0.001 &&
     track->GetTrkTrack()->IsInsideCavity() &&
     true ) LX6 = true;
  
  //====================================================
  // lever-arm 5
  //====================================================
  bool LX5=false;
  if(
     track->GetTrkTrack()->GetLeverArmX()==5 &&
     true ){
    if(
       track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(4)
       ){
      
      if(
	 !track->GetTrkTrack()->IsBad(0,0) &&
	 !track->GetTrkTrack()->IsBad(4,0) &&
	 track->GetTrkTrack()->resx[0]<0.001 &&
	 track->GetTrkTrack()->resx[4]<0.001 &&
	 track->GetTrkTrack()->IsInsideCavity() &&
	 true) LX5 = true;	  
    }else if (
	      track->GetTrkTrack()->XGood(1) && track->GetTrkTrack()->XGood(5)
	      ){
      
      if(
	 !track->GetTrkTrack()->IsBad(1,0) &&
	 !track->GetTrkTrack()->IsBad(5,0) &&
	 track->GetTrkTrack()->resx[1]<0.001 &&
	 track->GetTrkTrack()->resx[5]<0.001 &&
	 track->GetTrkTrack()->IsInsideCavity() &&
	 true) LX5 = true;
    }
  }
  if ( !LX5 && !LX6 ) return false;
  Float_t defl = trk->GetDeflection();
  float p0 = 1.111588e+00;
  float p1 = 1.707656e+00;
  float p2 = 1.489693e-01;
  float chi2m025 = p0 + fabs(defl)*p1 + defl*defl*p2;
  
  float def_0      = 0.07;
  float chi2m025_0 = p0 + fabs(def_0)*p1 + def_0*def_0*p2;
  
  //  int   nchi2cut=5;
  float chi2cut=3.;
  float chi2m = pow( chi2m025-chi2m025_0+pow(chi2cut,0.25), 4.);
  bool CUT2 = false;
  if(
     track->GetTrkTrack()->chi2 < chi2m &&
     true ) CUT2 = true;  
  if ( !CUT2 ) return false;
  float dedxtrk    = trk->GetDEDX();
  //  float zcutn =  9.   + 20./(rigidity*rigidity);
  float zcut2 =  3.   + 4.3/(rigidity*rigidity);
  float zcut1 =  0.52 + 0.455/(rigidity*rigidity);
  Bool_t Z1 = false;
  if(dedxtrk > zcut1 && dedxtrk < zcut2){
    Z1=true;
  }
  if ( !Z1 && !SIMU ) return false;
  //------------------------------------------------------
  //
  // energy momentum match
  //
  Float_t qtotimp = event->GetCaloLevel2()->qtot / trk->GetRigidity();
  Float_t qcut2 = (-0.5 * trk->GetRigidity() + 150.) * 1.1;
  if ( qcut2 < 55. ) qcut2 = 55.;
  if ( qtotimp <= qcut2 ) return false;
  //
  for (Int_t i=0; i < 22; i++){
    if ( track->GetCaloTrack()->tibar[i][1] < 0 || track->GetCaloTrack()->tibar[i][0] < 0 ){
      return false;
    };
  };
  //
  if ( event->GetCaloLevel2()->qtot == 0. ) return false;
  if ( rigidity>5. && track->GetCaloTrack()->qtrack/event->GetCaloLevel2()->qtot < 0.4 ) return false;
  if ( rigidity<1. && track->GetToFTrack()->beta[12] < 0.8 ) return false;
  if ( rigidity>50. ){
    if ( trk->GetNX()<5  &&
	 trk->GetNY()<4 ) return false;
    //
    Bool_t sphit = false;
    for ( Int_t plane = 0; plane < 6; plane++){
      if ( !trk->XGood(plane) ){
	for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsx(); sing++){
	  TClonesArray &t = *(event->GetTrkLevel2()->SingletX);
	  TrkSinglet *singlet = (TrkSinglet*)t[sing];
	  if ( (singlet->plane-1) == plane ){
	    Float_t x = (singlet->coord[0]+singlet->coord[1])/2.;
	    if ( fabs(track->GetTrkTrack()->xv[plane] - x) < 1. ) sphit = true;
	  };
	};
      };
      if ( !trk->YGood(plane) ){
	for (Int_t sing = 0; sing < event->GetTrkLevel2()->nclsy(); sing++){
	  TClonesArray &t = *(event->GetTrkLevel2()->SingletY);
	  TrkSinglet *singlet = (TrkSinglet*)t[sing];
	  if ( (singlet->plane-1) == plane ){
	    Float_t x1 = (singlet->coord[0]);
	    Float_t x2 = (singlet->coord[1]);
	    if ( fabs(track->GetTrkTrack()->yv[plane] - x1) < 1. ) sphit = true;
	    if ( fabs(track->GetTrkTrack()->yv[plane] - x2) < 1. ) sphit = true;
	  };
	};
      };
    };
    if ( sphit ) return false; // spurious hit along the track
  };
  //
  Int_t ti0 = track->GetCaloTrack()->tibar[0][1]-1;
  
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Float_t mip = 0.;
  //
  for ( Int_t i=0; i<event->GetCaloLevel1()->istrip; i++ ){
    //
    mip = event->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
    if ( view == 1 && plane == 0 && strip == ti0 && mip > 4.) return false;
    if ( view == 1 && (plane >0 || strip > ti0) ) break;
  };
  //  if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false;

  //  if ( rigidity > 2.2 || rigidity < 1.5 ) return false;
  //  printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG);
  //
  return true;
}
//===============================================================
// Create histograms
//
//
//
//
//
//===============================================================
void CreateHistos( PamLevel2* event , TString file){ 

  cf = new CaloFranzini(event);
  //
  if ( FULL ){
    cf->CalculateLongTZeta(false);
    cf->CalculateFullTZeta();
  };
  //
  if ( MATRIX ){
    cf->UpdateMatrixFile(file.Data());
    if ( !FULL ){
      cf->LoadBin();
      cf->LoadLong();
    } else {
      cf->LoadBin(true);
      cf->LoadFull(); 
    };
  } else {
    cf->CreateMatrixFile(file.Data());
  };
  //
  //
  nbin = 18;
  rig[0] = 0.1;
  rig[1] = 0.5;
  rig[2] = 1.;
  rig[3] = 1.5;
  rig[4] = 2.2;
  rig[5] = 3.;
  rig[6] = 4.;
  rig[7] = 5.;
  rig[8] = 6.;
  rig[9] = 8.;
  rig[10] = 10.;
  rig[11] = 15.;
  rig[12] = 25.;
  rig[13] = 35.;
  rig[14] = 50.;
  rig[15] = 100.;
  rig[16] = 200.;
  rig[17] = 300.;
  //
  memset(rmean, 0, 17*sizeof(Float_t));
  memset(ntot, 0, 17*sizeof(Int_t));
  //  memset(finmat, 0, 43*191*sizeof(Int_t));  
//  Double_t tol = 1E-20;
  //
  for (Int_t i=0; i < 17 ; i++){
    //  for (Int_t i=3; i < 4 ; i++){
    if ( !FULL ){
      matrix[i] = new TMatrixD(43,43);
      qplane[i] = new TArrayF(43);
      nqplane[i] = new TArrayI(43);
      nmat[i] = new TMatrixD(43,43);
    } else {
      if ( MATRIX ){
	//	fmatrix = new TMatrixF(4128,4128);
	//	fnmat = new TMatrixF(4128,4128);
	//	fmatrix = new TMatrixF(8213,8213);
	//	fnmat = new TMatrixF(8213,8213);
	//	fmatrix = new TMatrixF(MDIM,MDIM);
	//	fnmat = new TMatrixF(MDIM,MDIM);
	//	fmatrix[i] = new TMatrixF(1849,1849);
	//	fnmat[i] = new TMatrixF(43,43);
	//	fmatrix[i] = new TMatrixD(1333,1333);
	//	fnmat[i] = new TMatrixF(43,31);
	fmatrix[i] = new TMatrixD(473,473);
	fnmat[i] = new TMatrixF(43,11);
//	fmatrix[i]->SetTol(tol);
	//	cf->WriteFullMatrix(fmatrix, i);
	//	cf->WriteFullNMatrix(fnmat, i);
	//	delete fmatrix;
	//	delete fnmat;
	//fnmat[i] = new TMatrixI(8213,8213);
      } else {
	//	fqplane = new TMatrixD(43,191); // 43 planes x 191 strip (= 1 + 95 x 2, one strip is the one transversed by the track that could be on the extreme right or left) 
	//	fnqplane = new TMatrixD(43,191);//
	//	fqplane[i] = new TMatrixD(43,43); // 43 planes x 43 "strip", where 43 = 50 + 14 + 5 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + [1] + ...
	//	fnqplane[i] = new TMatrixD(43,43);//                                     0    1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20   21
	//
	//	fqplane[i] = new TMatrixD(43,31); // 43 planes x 43 "strip", where 43 = 50 + 14 + 6 + 5 + 3 + 3 + 3 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + [1] + 1 + 1 + 1 + 1 + ...
	//	fnqplane[i] = new TMatrixD(43,31);//                                     0    1   2   3   4   5   6   7   8   9  10  11  12  13  14    15  16  17  18  19  ...
	fqplane[i] = new TMatrixD(43,11); // 43 planes x 11 "strip", where 43 = 84 +  6 + 3 + 1 + 1 +[1]+ 1 + 1 + 3 + 6 + 84
	fnqplane[i] = new TMatrixD(43,11);//                                     0    1   2   3   4   5   6   7   8   9  10 
	//
	//	cf->WriteFullMean(fqplane, i);
	//	cf->WriteFullNMean(fnqplane, i);
	//	delete fqplane;
	//	delete fnqplane;
	//
      };
    };
  };
  //
}

//=============================================================== 
void FindAverage( PamLevel2* L2, int iev ){
  //
  Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
  if ( SRIG ) erig = L2->GetGPamela()->P0;
  if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
  //
  Int_t rbi = 0;
  for (Int_t i = 0; i<nbin-1; i++){ 
    if ( erig>=rig[i] && erig < rig[i+1] ){
      rbi = i;
      break;
    };
  };
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  rmean[rbi] += erig;
  ntot[rbi]++;
  //
  if (!FULL ){
    Int_t dgf = 43;
    //
    for (Int_t i=0; i<dgf; i++){
      (*nqplane[rbi])[i]++;
    };
    //
    // Fill the estrip matrix
    //
    Int_t nplane = 0;
    Int_t view = 0;
    Int_t plane = 0;
    Int_t strip = 0;
    Float_t mip = 0.;
    //
    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
      //
      nplane = 1 - view + 2 * plane;
      if ( erig > 4. && nplane == 0 && mip > 15. ) printf(" IEV %i erig %f OBT %u pkt %u file %s \n",iev,erig,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
      //printf(" IEV %i OBT %u pkt %u file %s \n",iev,L2->GetOrbitalInfo()->OBT,L2->GetOrbitalInfo()->pkt_num,L2->GetPamTree()->GetFile()->GetName());
      if ( nplane > 37 ) nplane--; 
      if ( nplane < dgf ){
	(*qplane[rbi])[nplane] += mip;
      };
      //
    };
  } else {
    //
    // FULL CALORIMETER
    //
    //    fqplane = cf->LoadFullAverage(rbi);
    //    fnqplane = cf->LoadFullNAverage(rbi);
    CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
    //
    Int_t nplane = 0;
    Int_t view = 0;
    Int_t plane = 0;
    Int_t strip = 0;
    Float_t mip = 0.;
    //
    Int_t cs = 0;
    Int_t cd = 0;
    Int_t mstrip = 0;
    //
    for (Int_t j=0; j<2; j++){
      for (Int_t i=0; i<22; i++){
	nplane = 1 - j + 2*i;
	if ( nplane > 37 ) nplane--;
	//
	cs = ct->tibar[i][j] - 1;
	//
	cd = 95 - cs;	
	//
	Int_t oldstr = -1;
	for (Int_t k=0; k<191; k++){
	  mstrip = cd + k;
	  //	  if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.;
	  //	  if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.;
	  Int_t lstr = cf->ConvertStrip(mstrip);
	  if ( oldstr != lstr ){
	    (*fnqplane[rbi])[nplane][lstr] += 1.;
	    oldstr = lstr;
	  };	
	};
      };
    };
    //
    //
    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
      //
      nplane = 1 - view + 2 * plane;
      if ( nplane > 37 ) nplane--; 
      //
      cs = ct->tibar[plane][view] - 1;
      //
      cd = 95 - cs;
      //
      mstrip = cd + strip;
      //
      Int_t lstr = cf->ConvertStrip(mstrip);
      //      (*fqplane[rbi])[nplane][mstrip] += mip;
      //      (*fqplane)[nplane][mstrip] += mip;
      (*fqplane[rbi])[nplane][lstr] += mip;
      //
    };
    //
    //    cf->WriteFullMean(fqplane, rbi);
    //    cf->WriteFullNMean(fnqplane, rbi);
    //    cf->UnLoadFullAverage(rbi);
    //    cf->UnLoadFullNAverage(rbi);
    //    delete fqplane;
    //    delete fnqplane;
    //
  };
}

void CalculateAverage(){
  //  
  if ( !FULL ){
    for (Int_t i=0; i<nbin-1; i++){
      if ( (*nqplane[i])[0] > 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0];
      for (Int_t j=0; j<43 ; j++){
	if ( (*nqplane[i])[j] > 0 ){  
	  (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j];
	} else {
	  (*qplane[i])[j] = 0.;
	};
	printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]);
      };
    };
    for (Int_t i=0; i<nbin-1; i++){
      //
      cf->WriteLongMean(qplane[i], i);
      //
    };
  } else {
    //
    for (Int_t i=0; i<nbin-1; i++){
      //      fqplane = cf->LoadFullAverage(i);
      //      fnqplane = cf->LoadFullNAverage(i);
      if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]);
      //
      for (Int_t j=0; j<43 ; j++){
	//	for (Int_t k=0; k<191; k++){
	//	for (Int_t k=0; k<43; k++){
	//	for (Int_t k=0; k<31; k++){
	for (Int_t k=0; k<11; k++){
	  // 	  if ( (*fnqplane[i])[j][k] > 0 ){  
	  // 	    (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
	  // 	  } else {
	  // 	    (*fqplane[i])[j][k] = 0.;
	  // 	  };
	  // 	  printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane[i])[j][k],(*fnqplane[i])[j][k]);
	  if ( (*fnqplane[i])[j][k] > 0 ){  
	    if ( (*fqplane[i])[j][k] == 0. ) (*fqplane[i])[j][k] = 0.7;
	    (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k];
	  } else {
	    (*fqplane[i])[j][k] = 0.;
	  };
	  //	  printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane)[j][k],(*fnqplane)[j][k]);
	};
      };
      cf->WriteFullMean(fqplane[i], i);
      cf->WriteFullNMean(fnqplane[i], i);
      //      cf->UnLoadFullAverage(i);
      //      cf->UnLoadFullNAverage(i);
      //      delete fqplane;
      //      delete fnqplane;
    };
    //
    //     for (Int_t i=0; i<nbin-1; i++){
    //       //
    //       cf->WriteFullMean(fqplane[i], i);
    //       //
    //     };
  };
  //  
  cf->WriteNumBin(nbin);
  //
  TArrayF *rigbin = new TArrayF(18, rig);
  cf->WriteRigBin(rigbin);
  //
  TArrayF *rmeanbin = new TArrayF(17, rmean);
  if ( FULL ){
    TFile *file =  cf->GetFullFile();
    file->cd();
    file->WriteObject(&(*rmeanbin), "binrigmean");
  } else {
    TFile *file =  cf->GetFile();
    file->cd();
    file->WriteObject(&(*rmeanbin), "binrigmean");
  };
  //
  //
}

//=============================================================== 
void FindMatrix( PamLevel2* L2, int iev ){
  //
  Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
  if ( SRIG ) erig = L2->GetGPamela()->P0;
  if ( CRIG ) erig = L2->GetCaloLevel2()->qtot/260.;
  //
  Int_t rbi = 0;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( erig>=rig[i] && erig < rig[i+1] ){
      rbi = i;
      break;
    };
  };
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  if ( !FULL ){
    Int_t dgf = 43;
    //
    for (Int_t i=0; i<dgf; i++){
      for (Int_t j=0; j<dgf; j++){
	(*nmat[rbi])[i][j] += 1.;
      };
    };
    //
    // Fill the estrip matrix
    //
    Int_t nplane = 0;
    Int_t view = 0;
    Int_t plane = 0;
    Int_t strip = 0;
    Float_t mip = 0.;
    Float_t hpl[43];
    memset(hpl,0,43*sizeof(Float_t));
    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
      //
      nplane = 1 - view + 2 * plane;
      if ( nplane > 37 ) nplane--; 
      if ( nplane < dgf ){
	hpl[nplane] += mip;
      };
      //
    };
    //
    for (Int_t i=0; i<dgf; i++){
      for (Int_t j=0; j<dgf; j++){
	(*matrix[rbi])[i][j] += (hpl[i] - cf->GetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig));
      };
    };
  } else {
    //
    // FULL CALORIMETER
    //
    //    if ( rbi != 3 ) return;
    printf(" matrix %i IEV %i \n",rbi,iev);
    //    fmatrix = cf->LoadFullMatrix(rbi);
    //    cf->LoadFullMatrix(rbi,fmatrix);
    //    fnmat = cf->LoadFullNMatrix(rbi);
    //    printf(" done \n");
    //    printf(" start loop \n");
    //
    CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
    //
    Int_t nplane = 0;
    Int_t view = 0;
    Int_t plane = 0;
    Int_t strip = 0;
    Float_t mip = 0.;
    //
    //    Int_t mindgf = 48;
    //    Int_t dgf = 143;
    //       Int_t mindgf = 0; //tutto
    //        Int_t dgf = 191; //tutto
    //    Int_t mindgf = 94;
    //    Int_t dgf = 96;
    //    Int_t mindgf = 84;
    //    Int_t dgf = 106;
    Int_t mindgf = 0;
    //    Int_t dgf = 43;
    Int_t dgf = 11;
    Int_t cs = 0;
    Int_t cd = 0;
    Int_t mstrip = 0;
    //
    //    Float_t mipv[43][43];
    //    memset(mipv,0,43*43*sizeof(Float_t));
    //    Float_t mipv[43][31];
    //    memset(mipv,0,43*31*sizeof(Float_t));
    Float_t mipv[43][11];
    memset(mipv,0,43*11*sizeof(Float_t));
    //
    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
      //
      nplane = 1 - view + 2 * plane;
      if ( nplane > 37 ) nplane--; 
      //
      cs = ct->tibar[plane][view] - 1;
      //
      cd = 95 - cs;
      //
      mstrip = cd + strip;
      //
      Int_t lstr = cf->ConvertStrip(mstrip);
      mipv[nplane][lstr] += mip;
      //
    };
    //
    Float_t mip1 = 1.;
    Int_t cs1;
    Int_t cd1;
    Float_t mip2 = 1.;
    Int_t cs2;
    Int_t cd2;
    Int_t mi = -1;
    Int_t mj = -1;
    Int_t nn1 = 0;
    Int_t pl1 = 0;
    Int_t vi1 = 0;
    Int_t nn2 = 0;
    Int_t pl2 = 0;
    Int_t vi2 = 0;
    Int_t mstrip1min = 0;
    Int_t mstrip1max = 0;
    Int_t mstrip2min = 0;
    Int_t mstrip2max = 0;
    //
    Int_t toto = 0;
    //
    //
    //
    Int_t therigb = rbi;
    //
    if ( erig < rig[0] ){
      therigb = 0;
    };
    if ( erig >= rig[nbin-2] ){
      therigb = nbin - 3;
    };
    //
    for (Int_t nplane1=0; nplane1<43; nplane1++){
      if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
      vi1 = 1;
      if ( nn1%2 ) vi1 = 0;
      pl1 = (nn1 - 1 + vi1)/2;
      //
      cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1
      //
      cd1 = 95 - cs1;
      //
      Int_t at1 = TMath::Max(0,(cd1+0));
      Int_t at2 = TMath::Min(190,(cd1+95));
      mstrip1min = cf->ConvertStrip(at1);
      mstrip1max = cf->ConvertStrip(at2) + 1;
      //      mstrip1min = cf->ConvertStrip(TMath::Max(mindgf,(cd1+0)));
      //      mstrip1max = cf->ConvertStrip(TMath::Min(dgf,(cd1+95))) + 1;
      //
      if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i mindgf %i dgf %i cd1 %i\n",nplane1,mstrip1min,mstrip1max,mindgf,dgf,cd1);
      //
      for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
	//	printf(".\n");
	//
	mj = -1;
	//
	mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,therigb);
	//
	//	mi = (nplane1 * 191) + mstrip1;
	//	mi = (nplane1 * 43) + mstrip1;
	//	mi = (nplane1 * 31) + mstrip1;
	mi = (nplane1 * 11) + mstrip1;
	//
	//	if ( mstrip1 > mstrip1min ) break;
	//	if ( mstrip1 > dgf ) break;
	//	if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){
	//
	//	finmat[nplane1][mstrip1]++;
	(*fnmat[rbi])[nplane1][mstrip1] += 1.;
	//
	if ( mip1 != 0. ){
	  //
	  for (Int_t nplane2=0; nplane2<43; nplane2++){
	    //
	    if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
	    vi2 = 1;
	    if ( nn2%2 ) vi2 = 0;
	    pl1 = (nn2 - 1 + vi2)/2;
	    //
	    cs2 = ct->tibar[pl2][vi2] - 1;
	    //
	    cd2 = 95 - cs2;
	    //
	    //	  mstrip2min = cd2 + 0;
	    //	  mstrip2max = cd2 + 95;
	    Int_t t1 = TMath::Max(0,(cd2+0));
	    Int_t t2 = TMath::Min(190,(cd2+95));
	    mstrip2min = cf->ConvertStrip(t1);
	    mstrip2max = cf->ConvertStrip(t2) + 1;
	    //
	    if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max);
	    //
	    for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
	      //	      
	      mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,therigb);
	      //
	      if ( mip2 != 0. ){
		//
		//		mj = (nplane2 * 191) + mstrip2;
		//		mj = (nplane2 * 43) + mstrip2;
//		mj = (nplane2 * 31) + mstrip2;
//		Int_t sh = -15 + nplane1;
//		if ( sh > 15 ) sh -= 31*nplane1;
//		//
//		mj = (nplane2 * 31) + mstrip2 + sh;
//		//
//		if ( mj < 0 ) mj += 1333;
//		if ( mj >= 1333 ) mj -= 1333;
//
//
//
//
// 		Int_t sh = -5 + nplane1;
// 		if ( sh > 5 ) sh -= 11*nplane1;
// 		//
// 		mj = (nplane2 * 11) + mstrip2 + sh;
// 		//
// 		if ( mj < 0 ) mj += 473;
// 		if ( mj >= 473 ) mj -= 473;
// 		//
 		//
 		mj = (nplane2 * 11) + mstrip2;
		//
//		printf(" mi %i mj %i sh %i \n",mi,mj,sh);
	        //
		//	      mj++;
		//
		//	      if ( mstrip2 > mstrip2min ) break;
		//	      if ( mstrip2 > dgf ) break;
		//	      if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){
		//		if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){
		//		(*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig));
		//	    	  (*fnmat[rbi])[mi][mj] += 1.;
		(*fmatrix[rbi])[mi][mj] += (mip1 * mip2); // giusto
		//		(*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.;
		toto++;
		//	      (*fmatrix)[mi][mj] += 1.;
		//		cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi);
		//		cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi);
		//	    	  (*fnmat)[mi][mj] += 1.;
		//		};
	      };
	    };
	  };
	};
      };
    };
    //  
    printf(" toto = %i \n",toto);
    printf("\n done \n");
    //    printf(" write matrix \n");
    //    cf->WriteFullMatrix(fmatrix, rbi);
    //    cf->WriteFullNMatrix(fnmat, rbi);
    //    printf(" done \n");
    //    printf(" unload matrix \n");
    //   cf->UnLoadFullMatrix(rbi);
    //    cf->UnLoadFullNMatrix(rbi);
    //    printf(" done \n");
    //    printf(" delete matrix \n");
    //    delete fmatrix;
    //    delete fnmat;
    //    printf(" done \n");
  };
}

//===============================================================
// Save histograms
//
//
//
//
//
//===============================================================
void SaveHistos(){
  //
  if ( MATRIX ){
    //
    printf("Finished, calculating average and inverting matrices\n");
    //
    if ( !FULL ){
      for (Int_t i=0; i<nbin-1; i++){
	//
	// determine the average matrix
	//    
	for (Int_t ii=0; ii<43; ii++){
	  for (Int_t j=0; j<43; j++){
	    if ( (*nmat[i])[ii][j] > 0. ){
	      (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
	    } else {
	      (*matrix[i])[ii][j] = 0.;
	    };
	  };
	};
	//
	cf->WriteLongMatrix(matrix[i],i);
	//
	if ( matrix[i]->Determinant() == 0. ){
	  printf("\n");
	  for (Int_t ii=0; ii<43; ii++){
	    for (Int_t j=0; j<43; j++){
	      printf(" %.f",(*matrix[i])[ii][j]);	  
	    };
	    printf("\n");
	  };
	  printf("\n");
	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
	} else {
	  Double_t det = 0.;
	  TMatrixD invmatrix = (TMatrixD)(matrix[i]->Invert(&det));
	  printf(" Bin %i determinant is %f \n",i,det);
	  cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i);
	};
      };
    } else {
      //
      // FULL
      //
      for (Int_t i=0; i<nbin-1; i++){
	//      for (Int_t i=3; i<5; i++){
	//      for (Int_t i=3; i<4; i++){
	//
	// determine the average matrix
	//    
	//	fmatrix = cf->LoadFullMatrix(i);
	//	fnmat = cf->LoadFullNMatrix(i);	
	//
	// 	for (Int_t ii=0; ii<MDIM; ii++){
	// 	  for (Int_t j=0; j<MDIM; j++){
	// // 	    if ( (*fnmat[i])[ii][j] > 0. ){
	// // 	      (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
	// // 	    } else {
	// // 	      (*fmatrix[i])[ii][j] = 0.;
	// // 	    };
	// 	    if ( (*fnmat)[ii][j] > 0. ){
	// 	      (*fmatrix)[ii][j] /= (*fnmat)[ii][j];
	// 	    } else {
	// 	      (*fmatrix)[ii][j] = 0.;
	// 	    };
	// 	  };
	// 	};
	//
// 	TMatrixD *mymat3 = new TMatrixD(129,129);
// 	TMatrixD *mymat5 = new TMatrixD(215,215);
// 	TMatrixD *mymat7 = new TMatrixD(301,301);
// 	TMatrixD *mymat9 = new TMatrixD(387,387);
// 	TMatrixD *mymat11 = new TMatrixD(473,473);
// 	TMatrixD *mymat17 = new TMatrixD(731,731);
	//	TMatrixF *mymat = new TMatrixF(129,129);
	//	TMatrixF *mymat = new TMatrixF(989,989);
	Int_t i1 = -1;
	Int_t j1 = -1;
	//	int mi,mj;
	Int_t nonzero = 0;
	Int_t nonzero1 = 0;
	for (Int_t ii=0; ii<43; ii++){
	  //	  for (Int_t j=0; j<191; j++){
	  //	  for (Int_t j=0; j<43; j++){
	  //	  for (Int_t j=0; j<31; j++){
	  for (Int_t j=0; j<11; j++){
	    // 	    if ( (*fnmat[i])[ii][j] > 0. ){
	    // 	      (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j];
	    // 	    } else {
	    // 	      (*fmatrix[i])[ii][j] = 0.;
	    // 	    };
	    //	    i1 =  (ii * 191) + j;
	    //	    i1 =  (ii * 43) + j;
	    //	    i1 =  (ii * 31) + j;
	    i1 =  (ii * 11) + j;
	    //	    j1 = -1;
	    for (Int_t iij=0; iij<43; iij++){
	      //	      for (Int_t jj=0; jj<191; jj++){
	      //	      for (Int_t jj=0; jj<43; jj++){
	      //	      for (Int_t jj=0; jj<31; jj++){
	      for (Int_t jj=0; jj<11; jj++){
		//
		//		j1 = (iij * 191) + jj;
		//		j1 = (iij * 43) + jj;
// 		Int_t sh = -15 + ii;
// 		if ( sh > 15 ) sh -= 31*ii;
// 		//
// 		j1 = (iij * 31) + jj + sh;
// 		//
// 		if ( j1 < 0 ) j1 += 1333;
// 		if ( j1 >= 1333 ) j1 -= 1333;
		//
		//
		//
// 		Int_t sh = -5 + ii;
// 		if ( sh > 5 ) sh -= 11*ii;
// 		//
// 		j1 = (iij * 11) + jj + sh;
// 		//
// 		if ( j1 < 0 ) j1 += 473;
// 		if ( j1 >= 473 ) j1 -= 473;
//
		j1 = (iij * 11) + jj;
		//		j1++;
		//		if ( finmat[ii][j] > 0 ){
		//		  (*fmatrix)[i1][j1] /= finmat[ii][j];
		if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix[i])[i1][j1] == 0. || !((*fmatrix[i])[i1][j1] == (*fmatrix[i])[i1][j1]) ){
		  (*fmatrix[i])[i1][j1] = 1.;
		} else {
		  (*fmatrix[i])[i1][j1] /= (*fnmat[i])[ii][j];
		  nonzero++;
		  if ( i1 == 0 ) nonzero1++;
		};
		//
// 		if ( j>=7 && j <=23 && jj >=7 && jj<=23 ){
// 		  Int_t mi17 = (ii*3) + j -7;
// 		  Int_t mj17 = (iij*3) + jj -7;
// 		  (*mymat17)[mi17][mj17] = (*fmatrix[i])[i1][j1];
// 		};
// 		if ( j>=10 && j <=20 && jj >=10 && jj<=20 ){
// 		  Int_t mi11 = (ii*3) + j -10;
// 		  Int_t mj11 = (iij*3) + jj -10;
// 		  (*mymat11)[mi11][mj11] = (*fmatrix[i])[i1][j1];
// 		};
// 		if ( j>=11 && j <=19 && jj >=11 && jj<=19 ){
// 		  Int_t mi9 = (ii*3) + j -11;
// 		  Int_t mj9 = (iij*3) + jj -11;
// 		  (*mymat9)[mi9][mj9] = (*fmatrix[i])[i1][j1];
// 		};
// 		if ( j>=12 && j <=18 && jj >=12 && jj<=18 ){
// 		  Int_t mi7 = (ii*3) + j -12;
// 		  Int_t mj7 = (iij*3) + jj -12;
// 		  (*mymat7)[mi7][mj7] = (*fmatrix[i])[i1][j1];
// 		};
// 		if ( j>=13 && j <=17 && jj >=13 && jj<=17 ){
// 		  Int_t mi5 = (ii*3) + j -13;
// 		  Int_t mj5 = (iij*3) + jj -13;
// 		  (*mymat5)[mi5][mj5] = (*fmatrix[i])[i1][j1];
// 		};
// 		if ( j>=14 && j <=16 && jj >=14 && jj<=16 ){
// 		  Int_t mi3 = (ii*3) + j -14;
// 		  Int_t mj3 = (iij*3) + jj -14;
// 		  (*mymat3)[mi3][mj3] = (*fmatrix[i])[i1][j1];
// 		};


		//		if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){
		//			mi = (ii*3) + j -94;
		//			mj = (iij*3) + jj -94;
		//			(*mymat)[mi][mj] = (*fmatrix)[i1][j1];
		//		};


		//		if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){
		//			mi = (ii*3) + j -84;
		//			mj = (iij*3) + jj -84;
		//			(*mymat)[mi][mj] = (*fmatrix)[i1][j1];
		//		};

	      };
	    };
	  };
	};
	//
	printf(" Matrix has %i non-zero elements \n",nonzero);	
	//	printf(" Matrix has %i non-zero elements on the first row\n",nonzero1);	
	//
	// 	Bool_t BAD = false;
	// 	for (Int_t ii=0; ii<43; ii++){
	// 	  for (Int_t j=0; j<191; j++){
	// 	    //
	// 	    i1 =  (ii * 191) + j;
	// 	    //
	// 	    for (Int_t iij=0; iij<43; iij++){
	// 	      for (Int_t jj=0; jj<191; jj++){
	// 		//
	// 		j1 = (iij * 191) + jj;
	// 		//
	// 		//		printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
	// 		if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){
	// 		  printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]);
	// 		  printf(" che schifo! \n");
	// 		  BAD = true;
	// 		};
	// 		//
	// 	      };
	// 	    };
	// 	  };
	// 	};
	// 	//
	// 	if ( BAD ) printf(" questa matrice fa cagare \n");
	//
	//
	//	cf->WriteFullMatrix(fmatrix, i);
	//	cf->WriteFullNMatrix(fnmat, i);
	//cf->WriteFullMatrix(fmatrix[i],i);// OK
	//	cf->WriteFullNMatrix(fnmat[i], i); //OK
	//
	//	TDecompSVD svd(*fmatrix[i]); 
	//	Bool_t ok = svd.Decompose(); 
	//
	Double_t zero = (Double_t)0.0;
	//
	if ( fmatrix[i]->Determinant() == zero ){
	  //if ( fmatrix->Determinant() == 0. ){
	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
	  cf->WriteFullMatrix(fmatrix[i],i);// OK
	  cf->WriteFullNMatrix(fnmat[i], i); //OK
	} else {
	  //	};
	  //	if ( i == 3 ){
	  //	if ( ok ){
	  //	  Double_t tol = 1E-20;
	  //	  TDecompSVD svd((*fmatrix)[i],tol);
	  //	  svd.Decompose(); 
	  //	  TMatrixD svdInv = svd.Invert();
	  //	  svdInv.Print("svdInv"); 
	  //	  cout << "condition: " << svd.Condition() << endl; 
	  //	  cf->WriteInvertedFullMatrix((TMatrixD)svdInv,999);
	  
 	  Double_t det = 0.;
 	  TMatrixD invmatrix = (TMatrixD)(fmatrix[i]->Invert(&det));
 	  printf(" Bin %i determinant is %f \n",i,det);
          cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,i);
	};

// 	if ( mymat3->Determinant() == 0. ){
// 	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
// 	} else {
// 	  Double_t det = 0.;
// 	  TMatrixD invmatrix = (TMatrixD)(mymat3->Invert(&det));
// 	  printf(" Mymat3 determinant is %f \n",det);
// 	  cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1103);
// 	};
// 	cf->WriteFullMatrix(mymat3, 103);
// 	if ( mymat5->Determinant() == 0. ){
// 	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
// 	} else {
// 	  Double_t det = 0.;
// 	  TMatrixD invmatrix = (TMatrixD)(mymat5->Invert(&det));
// 	  printf(" Mymat5 determinant is %f \n",det);
// 	  cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1105);
// 	};
// 	cf->WriteFullMatrix(mymat5, 105);
// 	if ( mymat7->Determinant() == 0. ){
// 	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
// 	} else {
// 	  Double_t det = 0.;
// 	  TMatrixD invmatrix = (TMatrixD)(mymat7->Invert(&det));
// 	  printf(" Mymat7 determinant is %f \n",det);
// 	  cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1107);
// 	};
// 	cf->WriteFullMatrix(mymat7, 107);
// 	if ( mymat9->Determinant() == 0. ){
// 	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
// 	} else {
// 	  Double_t det = 0.;
// 	  TMatrixD invmatrix = (TMatrixD)(mymat9->Invert(&det));
// 	  printf(" Mymat3 determinant is %f \n",det);
// 	  cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1109);
// 	};
// 	cf->WriteFullMatrix(mymat9, 109);
// 	if ( mymat11->Determinant() == 0. ){
// 	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
// 	} else {
// 	  Double_t det = 0.;
// 	  TMatrixD invmatrix = (TMatrixD)(mymat11->Invert(&det));
// 	  printf(" Mymat11 determinant is %f \n",det);
// 	  cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1111);
// 	};
// 	cf->WriteFullMatrix(mymat11, 111);
// 	if ( mymat17->Determinant() == 0. ){
// 	  printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i);
// 	} else {
// 	  Double_t det = 0.;
// 	  TMatrixD invmatrix = (TMatrixD)(mymat17->Invert(&det));
// 	  printf(" Mymat3 determinant is %f \n",det);
// 	  cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1117);
// 	};
// 	cf->WriteFullMatrix(mymat17, 117);


	//
	//	cf->UnLoadFullMatrix(i);
	//	cf->UnLoadFullNMatrix(i);
	//	delete fmatrix;
	//	delete fnmat;
	//
      };
    };
    //
    printf(" done, closing file and exiting\n");
    //
  };
  //
  cf->CloseMatrixFile();
  //
  cf->Delete();
  //
}