#include <stdlib.h>
#include <iostream>
#include <iomanip>
//
#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMatrixD.h>
#include <TArrayF.h>
#include <TArrayI.h>
//
#include <PamLevel2.h>
#include <CaloFranzini.h>
//
using namespace std;
//
CaloFranzini *cf;
Int_t nbin;
Float_t rig[18];
Float_t rmean[18];
Float_t qqplane[18][43];
Int_t nnqplane[18][43];

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

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

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

  //------------------------------------------------------------------
  // tracker pre-selection
  //------------------------------------------------------------------
  TrkTrack *trk = track->GetTrkTrack();

  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 )return false;
  //------------------------------------------------------
  // no albedo
  //------------------------------------------------------
  if(  	 track->GetToFTrack()->beta[12]<=0.2 ||
	 track->GetToFTrack()->beta[12] >= 1.5 ) return false;

  //------------------------------------------------------
  bool CUT1 = false;
  if( 
     trk->nstep<100   &&
     trk->GetRigidity()<400.   &&
     trk->GetRigidity()>0.1   &&
     trk->GetDeflection()<0.   &&
     trk->resx[0]<0.001 && 
     trk->resx[5]<0.001 && 
     track->IsSolved() &&
     trk->IsInsideCavity() &&
     true ) CUT1 = true;
  //------------------------------------------------------
  if( !CUT1 )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;
  //------------------------------------------------------
  //
  // 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->GetCaloLevel1()->qtotpl(0) > 7. ) return false;
  //
  return true;
}
//===============================================================
// Create histograms
//
//
//
//
//
//===============================================================
void CreateHistos( PamLevel2* event , TString file){

  cf = new CaloFranzini(event);
  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] = 400.;
  //
//  memset(qplane, 0, 44*18*sizeof(Float_t));
  memset(rmean, 0, 18*sizeof(Float_t));
  memset(qqplane, 0, 43*18*sizeof(Float_t));
  memset(nnqplane, 0, 43*18*sizeof(Float_t));
  //
  for (Int_t i=0; i < 18 ; i++){
    matrix[i] = new TMatrixD(43,43);
    qplane[i] = new TArrayF(43);
    nqplane[i] = new TArrayI(43);
    nmat[i] = new TMatrixD(43,43);
  };
//  for (Int_t i=0; i < 18 ; i++){
//	for (Int_t j = 0; j < 43;i++){
//		(*qplane[i])[j] = 0.;
//		qplane[i]->Reset();
//		nqplane[i]->Reset();
//		for (Int_t k = 0; k < 43;k++){
//			(*matrix[i])[j][k] = 0.;
//			(*nmat[i])[j][k] = 0.;
//		};
//	};
//  };
  //
}

//=============================================================== 
void FindAverage( PamLevel2* L2, int iev ){
  //
  //  printf("SELEZIONATO \n");
  Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
  //
  Int_t rbi = 0;
  for (Int_t i = 0; i<nbin-1; i++){ // < 17
    if ( erig>=rig[i] && erig < rig[i+1] ){
      rbi = i;
      break;
    };
  };
  //
  //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  rmean[rbi] += erig;
  //
  Int_t dgf = 43;
  //
//  for (Int_t i=0; i < 22; i++){
//    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){
//      dgf = 2 * i;
 //     break;
  //  };
 //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){
 //     dgf = 1 + 2 * i;
 //     break;
 //   };
 // };
  //
//  if ( dgf < 43 && dgf > 37 ) dgf--; 
  //
  for (Int_t i=0; i<dgf; i++){
    (*nqplane[rbi])[i]++;
 	nnqplane[rbi][i]++;
  };
  //
  //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);
  //
  //
  // 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 ( nplane > 37 ) nplane--; 
    //    printf(" mip %f view %i plane %i => nplane %i \n",mip,view,plane,nplane);
    if ( nplane < dgf ){
      //      printf(" PRIMA: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);
      (*qplane[rbi])[nplane] += mip;
 	qqplane[rbi][nplane] += mip;
      //      printf(" DOPO: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]);
    };
    //
  };
}

//=============================================================== 
void FindMatrix( PamLevel2* L2, int iev ){
  //
  //  printf("SELEZIONATO \n");
  Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity();
  //
  Int_t rbi = 0;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( erig>=rig[i] && erig < rig[i+1] ){
      rbi = i;
      break;
    };
  };
  //
  //  printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]);
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  Int_t dgf = 43;
  //
//  for (Int_t i=0; i < 22; i++){
//    if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){
//      dgf = 2 * i;
 //     break;
 //   };
 //   if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){
 //     dgf = 1 + 2 * i;
 //     break;
 //   };
//  };
  //
//  if ( dgf < 43 && dgf > 37 ) dgf--; 
  //
  for (Int_t i=0; i<dgf; i++){
    for (Int_t j=0; j<dgf; j++){
      (*nmat[rbi])[i][j] += 1.;
    };
  };
  //
  //  printf(" SELEZIONATO rig %f dgf %i \n",erig,dgf);
  //
  //
  // 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++){
      // if ( rbi == 1 ) printf(" PRIMA: matrix %f qpli %f qplj %f hpli %f hplj %f \n",(*matrix[rbi])[i][j],(*qplane[rbi])[i],(*qplane[rbi])[j],hpl[i],hpl[j]);
      (*matrix[rbi])[i][j] += (hpl[i] - (*qplane[rbi])[i]) * (hpl[j] - (*qplane[rbi])[j]);
      //      if ( rbi == 1 ) printf(" DOPO: matrix %f qpli %f qplj %f hpli %f hplj %f \n",(*matrix[rbi])[i][j],(*qplane[rbi])[i],(*qplane[rbi])[j],hpl[i],hpl[j]);
    };
  };
}

void CalculateAverage(){
  for (Int_t i=0; i<nbin; 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]);
//      if ( nnqplane[i][j] > 0 ){  
//	qqplane[i][j] /= (Float_t)nnqplane[i][j];
//      } else {
//	qqplane[i][j] = 0.;
//      };
//      printf(" BBIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],qqplane[i][j],nnqplane[i][j]);
    };
  };
  
  cf->WriteNumBin(nbin);

  TArrayF *rigbin = new TArrayF(18, rig);
  cf->WriteRigBin(rigbin);

  TArrayF *rmeanbin = new TArrayF(18, rmean);
  TFile *file =  cf->GetFile();
  file->cd();
  //  rigbin->Write("binrig");
  file->WriteObject(&(*rmeanbin), "binrigmean");

//  cf->WriteRigBin(rigbin);
  
  for (Int_t i=0; i<nbin; i++){
    
    cf->WriteLongMean(qplane[i], i);

  };
}
//===============================================================
// Save histograms
//
//
//
//
//
//===============================================================
void SaveHistos(){
  
  printf("Finished, calculating average and inverting matrices\n");


  
//  cf->WriteNumBin(nbin);

//  TArrayF *rigbin = new TArrayF(18, rig);
//  cf->WriteRigBin(rigbin);
  
  for (Int_t i=0; i<nbin; i++){
    
//    cf->WriteLongMean(qplane[i], i);

    //
    // determine the average matrix
    //    
    for (Int_t ii=0; ii<43; ii++){
      for (Int_t j=0; j<43; j++){
	//	if ( i == 1 ) printf(" num %f \n",(*nmat[i])[ii][j]);
	if ( (*nmat[i])[ii][j] > 0. ){
	  //	  if ( i == 1 ) printf(" value %f \n",(*matrix[i])[ii][j]);
	  (*matrix[i])[ii][j] /= (*nmat[i])[ii][j];
	} else {
	  (*matrix[i])[ii][j] = 0.;
	};
      };
    };
    
//     if ( i == 2 ){
//       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");
//     };
    

    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);
    };
  };
  
  printf(" done, closing file and exiting\n");

  cf->CloseMatrixFile();

  cf->Delete();
}