#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;
//
extern Bool_t MATRIX;
CaloFranzini *cf;
Int_t nbin;
Float_t rig[18];
Float_t rmean[17];
//Float_t qqplane[17][43];
//Int_t nnqplane[17][43];

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

//===============================================================================
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);
  //
  if ( MATRIX ){
    cf->UpdateMatrixFile(file.Data());
    cf->LoadBin();
  } 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] = 400.;
  //
  memset(rmean, 0, 17*sizeof(Float_t));
  //
  for (Int_t i=0; i < 17 ; i++){
    matrix[i] = new TMatrixD(43,43);
    qplane[i] = new TArrayF(43);
    nqplane[i] = new TArrayI(43);
    nmat[i] = new TMatrixD(43,43);
  };
  //
}

//=============================================================== 
void FindAverage( PamLevel2* L2, int iev ){
  //
  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;
    };
  };
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  rmean[rbi] += erig;
  //
  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 ( nplane > 37 ) nplane--; 
    if ( nplane < dgf ){
      (*qplane[rbi])[nplane] += mip;
    };
    //
  };
}

void CalculateAverage(){
  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]);
    };
  };
  
  cf->WriteNumBin(nbin);

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

  TArrayF *rmeanbin = new TArrayF(17, rmean);
  TFile *file =  cf->GetFile();
  file->cd();
  file->WriteObject(&(*rmeanbin), "binrigmean");
  //
  for (Int_t i=0; i<nbin-1; i++){
    //
    cf->WriteLongMean(qplane[i], i);
    //
  };
  //
}

//=============================================================== 
void FindMatrix( PamLevel2* L2, int iev ){
  //
  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;
    };
  };
  //
  if ( erig < rig[0] ) return;
  if ( erig >= rig[nbin-1] ) return;
  //
  Int_t dgf = 43;
//  printf("1qui \n");
  //
  for (Int_t i=0; i<dgf; i++){
    for (Int_t j=0; j<dgf; j++){
      (*nmat[rbi])[i][j] += 1.;
    };
  };
//  printf("2qui \n");
  //
  // 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;
    };
    //
  };
  //
//  printf("3qui \n");
  //
  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));
    };
  };
//  printf("4qui \n");  
  //  gObjectTable->Print();
}

//===============================================================
// Save histograms
//
//
//
//
//
//===============================================================
void SaveHistos(){
  //
  if ( MATRIX ){
    //
    printf("Finished, calculating average and inverting matrices\n");
    //
    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);
      };
    };
    //
    printf(" done, closing file and exiting\n");
    //
  };
  //
  cf->CloseMatrixFile();
  //
  cf->Delete();
  //
}