/**
 * \file CaloFranzini.cpp
 * \author Emiliano Mocchiutti (2007/07/18)
 */
//
// headers
//
#include <CaloFranzini.h>
//--------------------------------------
/**
 * Default constructor 
 */
CaloFranzini::CaloFranzini(){
  Clear();
}

CaloFranzini::CaloFranzini(PamLevel2 *l2p){  
  //
  lfile = NULL;
  ffile = NULL;
  brig = NULL;
  brigm = NULL;
  nbin = 0;
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  // Default variables
  //
  debug = false;
  dolong = true;
  dofull = false;
  sntr = 0;
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  crig = false;
  sel = true;
  cont = false;
  N = 0;
  NC = 43;
  //
  mask18b = -1;
  //
  Clear();
  //
}

void CaloFranzini::Clear(){
  //
  longtzeta = 0.;
  fulltzeta = 0.;
  degfre = 0;
  fdegfre = 0;
  memset(estrip, 0, 2*22*96*sizeof(Float_t));
  memset(qplane, 0, 43*sizeof(Float_t));
  //
  numneg = 0;
  numpos = 0;
  negfulltzeta = 0.;
  posfulltzeta = 0.;
  aveposvar = 0.;
  avenegvar = 0.;
  minsvalue =  numeric_limits<UInt_t>::max();
  maxsvalue =  numeric_limits<UInt_t>::min();
  //
}

void CaloFranzini::Print(){
  //
  Process();
  //
  printf("========================================================================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
  printf(" debug                [debug flag]:.. %i\n",debug);
  printf(" long  degree of freedom          :.. %i\n",this->GetLongDegreeOfFreedom());  
  printf(" longtzeta                        :.. %f\n",longtzeta);  
  printf(" longtzeta normalized             :.. %f\n",this->GetNormLongTZeta());  
  printf(" full degree of freedom           :.. %i\n",this->GetFullDegreeOfFreedom());  
  printf(" fulltzeta                        :.. %f\n",fulltzeta);  
  printf(" fulltzeta normalized             :.. %f\n",this->GetNormFullTZeta());  
  printf(" fulltzeta negative contribute    :.. %f\n",negfulltzeta);   
  printf(" fulltzeta positive contribute    :.. %f\n",posfulltzeta);   
  printf(" fulltzeta number of negatives    :.. %i\n",numneg);   
  printf(" fulltzeta number of positives    :.. %i\n",numpos);   
  printf(" fulltzeta minimum variance       :.. %f\n",minsvalue);   
  printf(" fulltzeta maximum variance       :.. %f\n",maxsvalue);   
  printf(" fulltzeta average positive var   :.. %f\n",aveposvar);   
  printf(" fulltzeta average negative var   :.. %f\n",avenegvar);   
  printf("========================================================================\n");
  //
}

void CaloFranzini::Delete(){
  //
  if ( ffile ) ffile->Close();
  if ( lfile ) lfile->Close();
  //
  Clear();
  //
}

void CaloFranzini::SetNoWpreSampler(Int_t n){
  Int_t nc2 = NC/2;
  if ( NC >= 37 ) nc2 = (NC+1)/2;
  if ( nc2+n < 23 ){
    N = n;
  } else {
    printf(" ERROR! Calorimeter is made of 22 W planes\n");
    printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
    printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
    NC = 43;
    N = 0;
  };
}

void CaloFranzini::SetNoWcalo(Int_t n){
  if ( N+n < 23 ){
    NC = n*2;
    if ( NC >37 ) NC--;
  } else {
    printf(" ERROR! Calorimeter is made of 22 W planes\n");
    printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
    printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
    NC = 43;
    N = 0;
  };
}

void CaloFranzini::Process(){
  this->Process(0);
}

void CaloFranzini::Process(Int_t itr){
  //  
  if ( !L2 ){
    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
    printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
    return;
  };
  //
  if ( !nbin || !brig || (!lfile && !ffile) ){
    printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n");
    printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
    return;
  };
  //
  Bool_t newentry = false;
  //
  if ( L2->IsORB() ){
    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime || itr != sntr ){
      newentry = true;
      OBT = L2->GetOrbitalInfo()->OBT;
      PKT = L2->GetOrbitalInfo()->pkt_num;
      atime = L2->GetOrbitalInfo()->absTime;
      sntr = itr;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  // Some variables
  //
  if ( debug ) printf(" Processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  this->Clear();
  //
  Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity();
  if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.;
  //
  // Fill the estrip matrix
  //
  Int_t nplane = 0;
  Int_t view = 0;
  Int_t plane = 0;
  Int_t strip = 0;
  Float_t mip = 0.;
  //
  //
  //
  if ( dolong ){
    //
    for ( Int_t i=0; i<L2->GetCaloLevel1()->istrip; i++ ){
      //
      mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
      //
      estrip[view][plane][strip] = mip;
      //    
      nplane = 1 - view + 2 * (plane - N);
      if ( nplane > (37-(2*N)) ) nplane--; 
      //
      //    if ( plane == (18+N) ) mip = 0.;
      if ( nplane > -1 ) qplane[nplane] += mip;
      //
    };
    //
    if ( cont ){
      for (Int_t i=0; i<22; i++){
 	if ( i == (18+N) ){
	  mask18b =  1 + 2 * (i - N);
	  break;
	};
      };
    };  
    //  
    if ( sel ){
      for (Int_t i=0; i<22; i++){
 	if ( i == (18-N) ){
	  mask18b =  1 + 2 * (i - N);
	  break;
	};
      };
    };
    //
    if ( mask18b == 37 ) mask18b = -1;
    //
    Int_t dgf = 43;
    //
    for (Int_t i=0; i < 22; i++){
      if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){
	dgf = 2 * i;
	break;
      };
      if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){
      dgf = 1 + 2 * i;
      break;
      };
    };
    //
    if ( dgf < 43 && dgf > 37 ) dgf--; 
    // 
    degfre = TMath::Min(dgf,NC);
    //
    //    Float_t longzdiag = 0.;
    //    Float_t longzout = 0.;
    //
    if ( degfre > 0 ){
      for (Int_t i = 0; i < degfre; i++){
	if ( i != mask18b ){
	  for (Int_t j = 0; j < degfre; j++){	
	    if ( j != mask18b ){
// 	      if ( i == j ){
// 		longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
// 		if ( debug ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
// 	      } else {
// 		longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
// 		if ( debug && i == (j+1) ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig));
// 	      };
	      //
	      longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig));
	      //
	    };
	  };
	};
      };
      //      if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout);
    };
  };
  if ( dofull ){
    //    printf(" ERROR: NOT IMPLEMENTED YET\n");
    //    printf(" ERROR: CaloFranzini variables _NOT_ filled \n");
    //
    // FULL CALORIMETER
    //
    CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack();
    //
    Int_t dgf = 0;
    Int_t cs = 0;
    Int_t cd = 0;
    Int_t mstrip = 0;
    //
    Int_t maxnpl = 0;
    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--; 
      nplane = 1 - view + 2 * (plane - N);
      if ( nplane > (37-(2*N)) ) nplane--; 
      //
      cs = ct->tibar[plane][view] - 1;
      //
      if ( ct->tibar[plane][view] > -1 ){
	//
	cd = 95 - cs;
	//
	mstrip = cd + strip;
	//
	Int_t lstr = this->ConvertStrip(mstrip);
	//
	if ( nplane > maxnpl ) maxnpl = nplane;
	if ( nplane > -1 ) mipv[nplane][lstr] += mip;
	//
      };
      //      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 tonpl = TMath::Min(maxnpl,NC);
    //
    Int_t rbi = 0;
    for (Int_t i = 0; i<nbin-1; i++){
      if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
	rbi = i;
	break;
      };
    };
    //
    //
    Int_t therigb = rbi;
    //
    if ( rig < brigm->At(2) ){
//      therigb = 0;
      therigb = 2;
    };
    //    if ( rig >= brigm->At(nbin-5) ){
    if ( rig >= brigm->At(nbin-2) ){
      therigb = nbin - 3;
      //      therigb = nbin - 5;
    };
    Int_t mtherig = therigb;
    if ( mtherig >= 13 ) mtherig = 12;
    //
    if ( debug ) printf(" rig %f brigm0 %f brigm n2 %f nbin %i \n",rig,brigm->At(0),brigm->At(nbin-2),nbin);
    //
    if ( cont ){
      for (Int_t i=0; i<22; i++){
 	if ( i == (18+N) ){
	  mask18b =  1 + 2 * (i - N);
	  break;
	};
      };
    };  
    //  
    if ( sel ){
      for (Int_t i=0; i<22; i++){
 	if ( i == (18-N) ){
	  mask18b =  1 + 2 * (i - N);
	  break;
	};
      };
    };
    //
    if ( mask18b == 37 ) mask18b = -1;
    //
    if ( debug ) printf("Start main loop \n");
    //
    for (Int_t nplane1=0; nplane1<tonpl; nplane1++){
      //
      if ( nplane1 != mask18b ){
        //
	if ( nplane1 >= 37 ) nn1 = nplane1 + 1;
	vi1 = 1;
	if ( nn1%2 ) vi1 = 0;
	pl1 = (nn1 - 1 + vi1)/2;
	//
	cs1 = ct->tibar[pl1][vi1] - 1; 
	//
	if ( ct->tibar[pl1][vi1] > -1 ){
          //
          dgf++;
	  //
	  cd1 = 95 - cs1;
	  //
	  Int_t at1 = TMath::Max(0,(cd1+0));
	  Int_t at2 = TMath::Min(190,(cd1+95));
	  mstrip1min = this->ConvertStrip(at1);
	  mstrip1max = this->ConvertStrip(at2) + 1;
	  //
	  for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){
	    //
	    mj = -1;
	    //
	    //	    if ( debug ) printf(" get mip1 %i %i %f %i \n",nplane1,mstrip1,rig,therigb);
	    mip1 = mipv[nplane1][mstrip1] - this->GetFullAverageAt(nplane1,mstrip1,rig,therigb);
	    if ( mip1 < 0 ){
	      numneg++;
	      avenegvar += mip1;
	    };
	    if ( mip1 >= 0 ){
	      numpos++;
	      aveposvar += mip1;
	    };
	    if ( mip1 < minsvalue ) minsvalue = mip1;
	    if ( mip1 >= maxsvalue ) maxsvalue = mip1;
	    //
	    mi = (nplane1 * 11) + mstrip1;
	    //
	    if ( mip1 != 0. ){
	      //
	      for (Int_t nplane2=0; nplane2<tonpl; nplane2++){
		//
		if ( nplane2 != mask18b ){
		  //
		  if ( nplane2 >= 37 ) nn2 = nplane2 + 1;
		  vi2 = 1;
		  if ( nn2%2 ) vi2 = 0;
		  pl1 = (nn2 - 1 + vi2)/2;
		  //
		  cs2 = ct->tibar[pl2][vi2] - 1;
		  //
		  if ( ct->tibar[pl2][vi2] > -1 ){
		    //
		    cd2 = 95 - cs2;
		    //
		    Int_t t1 = TMath::Max(0,(cd2+0));
		    Int_t t2 = TMath::Min(190,(cd2+95));
		    mstrip2min = this->ConvertStrip(t1);
		    mstrip2max = this->ConvertStrip(t2) + 1;
		    //
		    for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){
		      //	      
		      mip2 = mipv[nplane2][mstrip2] - this->GetFullAverageAt(nplane2,mstrip2,rig,therigb);
		      //
		      if ( mip2 != 0. ){
			//
// 			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;
			//
			Float_t svalue =  mip1 * this->GetFullHmatrixAt(mi,mj,rig,mtherig,therigb) * mip2;
			//			fulltzeta += mip1 * this->GetFullHmatrixAt(mi,mj,rig,therigb) * mip2;
			fulltzeta += svalue;
			if ( svalue < 0. ) negfulltzeta += svalue;
			if ( svalue > 0. ) posfulltzeta += svalue;
			//		(*fmatrix[rbi])[mi][mj] += (mip1 * mip2); 
			//			if ( mstrip1 == mstrip2 ) printf("\n\n=> nplane1 %i nplane2 %i mstrip1 %i mstrip2 %i \n => mip1 %f H %f mip2 %f \n => mipv1 %f ave1 %f mipv2 %f ave2 %f \n => rig %f rigbin %i mtherigb %i\n",nplane1,nplane2,mstrip1,mstrip2,mip1,this->GetFullHmatrixAt(mi,mj,rig,mtherig),mip2,mipv[nplane1][mstrip1],this->GetFullAverageAt(nplane1,mstrip1,rig,therigb),mipv[nplane2][mstrip2],this->GetFullAverageAt(nplane2,mstrip2,rig,therigb),rig,therigb,mtherig);
			//
		      };
		    };
		  };
		};
	      };
	    };
	  };
	};
      };
    };
    //
    fdegfre = dgf*11;
    if ( numpos ) aveposvar /= numpos;
    if ( numneg ) avenegvar /= numneg;
    //
  };
  //
  //  if ( debug ) this->Print();
  if ( debug ) printf(" exit \n");
}


Float_t CaloFranzini::GetNormLongTZeta(){
  Process(); 
  Float_t normz = 0.;
  if ( degfre != 0 ) normz = longtzeta/(Float_t)degfre;
  return normz; 
}

Float_t CaloFranzini::GetNormFullTZeta(){
  Process(); 
  Float_t normz = 0.;
  if ( fdegfre != 0 ) normz = fulltzeta/(Float_t)fdegfre;
  return normz; 
}

Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){  
  //
  lfile = new TFile(matrixfile.Data(),"READ");
  //
  if ( !lfile || lfile->IsZombie() ){
    if ( dolong ){
      lfile = new TFile(matrixfile.Data(),"RECREATE"); 
      printf(" Create file %s \n",lfile->GetName());
    };
    if ( dofull ){
      ffile = new TFile(matrixfile.Data(),"RECREATE"); 
      printf(" Create file %s \n",ffile->GetName());
    };
  } else {
    lfile->Close();
    printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
    return(false);
  };
  //
  return(true);
  //
}

Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){  
  //
  if ( dolong ){
    lfile = new TFile(matrixfile.Data(),"UPDATE");
    //
    if ( !lfile || lfile->IsZombie() ){
      printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
      return(false);
    };
  };
  //
  if ( dofull ){
    ffile = new TFile(matrixfile.Data(),"UPDATE");
    //
    if ( !ffile || ffile->IsZombie() ){
      printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data());
      return(false);
    };
  };
  //
  return(true);
  //
}

void CaloFranzini::WriteNumBin(Int_t numbin){
  if ( dolong ){
    lfile->cd();
    TArrayI nbi(1, &numbin);
    lfile->WriteObject(&nbi, "nbinenergy");
  };
  if ( dofull ){
    ffile->cd();
    TArrayI nbi(1, &numbin);
    ffile->WriteObject(&nbi, "nbinenergy");
  };
}

void CaloFranzini::WriteRigBin(TArrayF *rigbin){
  if ( dolong ){
    lfile->cd();
    //  rigbin->Write("binrig");
    lfile->WriteObject(&(*rigbin), "binrig");
  };
  if ( dofull ){
    ffile->cd();
    //  rigbin->Write("binrig");
    ffile->WriteObject(&(*rigbin), "binrig");
  };
}

void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){
  lfile->cd();
  TString name = Form("qplmeann%i",bin);
  lfile->WriteObject(&(*qpl),name.Data());
}

void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){
  ffile->cd();
  TString name = Form("fqplmeann%i",bin);
  ffile->WriteObject(&(*qpl),name.Data());
  ffile->Purge();
}

void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){
  ffile->cd();
  TString name = Form("fnqplmeann%i",bin);
  ffile->WriteObject(&(*qpl),name.Data());
  ffile->Purge();
}

void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){
  lfile->cd();
  TString name = Form("matrixn%i",bin);
  //  mat.Write(name.Data());
  lfile->WriteObject(&mat,name.Data());
}

void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){
  ffile->cd();
  TString name = Form("fmatrixn%i",bin);
  //  mat.Write(name.Data());
  ffile->WriteObject(&mat,name.Data());
}

void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){
  lfile->cd();
  TString name = Form("origmatrixn%i",bin);
  //  mat.Write(name.Data());
  lfile->WriteObject(&(*mat),name.Data());
}

void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){
  ffile->cd();
  TString name = Form("origfmatrixn%i",bin);
  //  mat.Write(name.Data());
  ffile->WriteObject(&(*mat),name.Data());
  ffile->Purge();
}

void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){
  ffile->cd();
  TString name = Form("origfnmatrixn%i",bin);
  //  mat.Write(name.Data());
  ffile->WriteObject(&(*mat),name.Data());
  ffile->Purge();
}

void CaloFranzini::CloseMatrixFile(){
  if ( dolong ){
    lfile->cd();
    lfile->Close();
  };
  if ( dofull ){
    ffile->cd();
    ffile->Close();
  };
}

Bool_t CaloFranzini::Open(TString matrixfile){  
  return(this->Open(matrixfile,""));
}

Bool_t CaloFranzini::Open(TString longmatrixfile, TString fullmatrixfile){  
  //
  // find matrix file
  //
  if ( !strcmp(longmatrixfile.Data(),"") ){
    if (dolong) longmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root";    
  };
  if ( !strcmp(fullmatrixfile.Data(),"") ){
    if (dofull) fullmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root";    
  };
  //
  if ( dolong ){
    //
    lfile = new TFile(longmatrixfile.Data(),"READ");
    //
    if ( !lfile || lfile->IsZombie() ){
      printf(" ERROR: cannot open file %s \n",longmatrixfile.Data());
      return(false);
    };
    //
    if ( !this->LoadBin() ){
      printf(" %s \n",longmatrixfile.Data());
      return(false);
    };
    //
    if ( !this->LoadLong() ){
      printf(" %s \n",longmatrixfile.Data());
      return(false);
    };
    //
    if ( !this->LoadMatrices() ){
      printf(" %s \n",longmatrixfile.Data());
      return(false);
    };
  };
  //
  if ( dofull ){
    //
    ffile = new TFile(fullmatrixfile.Data(),"READ");
    //
    if ( !ffile || ffile->IsZombie() ){
      printf(" ERROR: cannot open file %s \n",fullmatrixfile.Data());
      return(false);
    };
    //
    if ( !dolong ){
      //
      if ( !this->LoadBin(true) ){
	printf(" %s \n",fullmatrixfile.Data());
	return(false);
      };      
      //
    };
    if ( !this->LoadFull() ){
      printf(" %s \n",fullmatrixfile.Data());
      return(false);
    };
    //
    if ( !this->LoadFullMatrices() ){
      printf(" %s \n",fullmatrixfile.Data());
      return(false);
    };
  };
  //
  //
  return(true);
  //
}


Bool_t CaloFranzini::LoadBin(){
  return(this->LoadBin(false));
}

Bool_t CaloFranzini::LoadBin(Bool_t full){
  //
  if ( !full ) {
    //
    TArrayI *numbin = (TArrayI*)lfile->Get("nbinenergy");
    if ( !numbin ){
      printf(" ERROR: cannot read number of bins from file ");
      return(false);
    };
    nbin = (Int_t)numbin->At(0);
    if ( nbin <= 0 ){
      printf(" ERROR: cannot work with 0 energy bins from file ");
      return(false);
    };
    //
    brig = (TArrayF*)lfile->Get("binrig");
    if ( !brig ){
      printf(" ERROR: cannot read rigidity binning from file ");
      return(false);
    };
    //
    brigm=(TArrayF*)lfile->Get("binrigmean");
    if ( !brigm ){
      printf(" ERROR: cannot read mean rigidity binning from file ");
      return(false);
    };
    //
  } else {
    //
    TArrayI *numbin = (TArrayI*)ffile->Get("nbinenergy");
    if ( !numbin ){
      printf(" ERROR: cannot read number of bins from file ");
      return(false);
    };
    nbin = (Int_t)numbin->At(0);
    if ( nbin <= 0 ){
      printf(" ERROR: cannot work with 0 energy bins from file ");
      return(false);
    };
    //
    brig = (TArrayF*)ffile->Get("binrig");
    if ( !brig ){
      printf(" ERROR: cannot read rigidity binning from file ");
      return(false);
    };
    //
    brigm=(TArrayF*)ffile->Get("binrigmean");
    if ( !brigm ){
      printf(" ERROR: cannot read mean rigidity binning from file ");
      return(false);
    };
  };
  //
  return(true);
}

Bool_t CaloFranzini::LoadLong(){
  //
  for (Int_t i=0;i<17;i++){
    TString name = Form("qplmeann%i",i);
    qplmean[i] = (TArrayF*)lfile->Get(name.Data());
    if ( !qplmean[i] ){
      printf(" ERROR: cannot read average from file ");
      return(false);
    };
  };
  //
  return(true);
}

Bool_t CaloFranzini::LoadFull(){
  //
  for (Int_t i=0;i<17;i++){
    TString name = Form("fqplmeann%i",i);
    fqplmean[i] = (TMatrixD*)ffile->Get(name.Data());
    if ( !fqplmean[i] ){
      printf(" ERROR: cannot read average from file ");
      return(false);
    };
  };
  //
  return(true);
}

Bool_t CaloFranzini::LoadMatrices(){
  //
  for (Int_t i=0;i<17;i++){
        TString name1 = Form("matrixn%i",i);
        hmat[i] = (TMatrixD*)lfile->Get(name1.Data());
  };
  //
  return(true);
}

Bool_t CaloFranzini::LoadFullMatrices(){
  //
  for (Int_t i=0;i<17;i++){
        TString name1 = Form("fmatrixn%i",i);
        hfmat[i] = (TMatrixD*)ffile->Get(name1.Data());
  };
  //
  return(true);
}

TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){
  //
  Int_t mv = 0;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
      mv = i;
      break;
    };
  };
  if ( rig < brig->At(0) ){
    printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
    mv = 0;
  };
  if ( rig >= brig->At(nbin-1) ){
    printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
    mv = nbin-2;
  };
  //
  return(hmat[mv]);
  //
}

TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){
  //
  TString name = Form("fqplmeann%i",rigbin);
  TMatrixD *fmean=(TMatrixD*)ffile->Get(name.Data());
  //
  return(fmean);
  //
}

TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){
  //
  TString name = Form("origfmatrixn%i",rigbin);
  TMatrixF *fmatri=(TMatrixF*)ffile->Get(name.Data());
  //
  return(fmatri);
  //
}

void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *&fmatri){
  //
  TString name = Form("origfmatrixn%i",rigbin);
  fmatri=(TMatrixF*)ffile->Get(name.Data());
  //
}

void CaloFranzini::UnLoadFullAverage(Int_t rigbin){
  //
  TString name = Form("fqplmeann%i",rigbin);
  ffile->Delete(name.Data());
  //
}

void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){
  //
  TString name = Form("origfmatrixn%i",rigbin);
  ffile->Delete(name.Data());
  //
}

TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){
  //
  TString name = Form("origfnmatrixn%i",rigbin);
  TMatrixF *fnmatri=(TMatrixF*)ffile->Get(name.Data());
  //
  return(fnmatri);
  //
}

void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){
  //
  TString name = Form("origfnmatrixn%i",rigbin);
  ffile->Delete(name.Data());
  //
}

TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){
  //
  TString name = Form("fnqplmeann%i",rigbin);
  TMatrixD *fnmean=(TMatrixD*)ffile->Get(name.Data());
  //
  return(fnmean);
  //
}

void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){
  //
  TString name = Form("fnqplmeann%i",rigbin);
  ffile->Delete(name.Data());
  //
}


TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){
  //
  Int_t mv=0;
  for (Int_t i = 0; i<nbin-1; i++){
    if ( rig>=brig->At(i) && rig < brig->At(i+1) ){
	mv = i;
      break;
    };
  };
  if ( rig < brig->At(0) ){
    printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0));
    mv = 0;
  };
  if ( rig >= brig->At(nbin-1) ){
    printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1));
    mv=nbin-2;
  };
  //
  return(qplmean[mv]);
  //
}

Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){
  //
  Int_t therigb = 0;
  for (Int_t i = 0; i<nbin-2; i++){
    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
      therigb = i;
      break;
    };
  };
  //
  Float_t minrig;
  Float_t maxrig;
  //
  //
  if ( rig < brigm->At(0) ){
    if ( rig < brig->At(0) ){
//      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
    };
    therigb = 0;
  };
  if ( rig >= brigm->At(nbin-2) ){
    if ( rig >= brig->At(nbin-2) ) {
//      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
    };
    therigb = nbin - 3;
  };
  //
  minrig = brigm->At(therigb);
  maxrig = brigm->At(therigb+1);
  //
  Float_t minene = (*qplmean[therigb])[plane];
  Float_t maxene = (*qplmean[therigb+1])[plane];
  //
  if ( maxrig == minrig ){
   printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
   return(0.);
  };
  Float_t b = log(maxene/minene)/(maxrig-minrig);
  Float_t a = minene/exp(minrig*b);
  Float_t ave = a*exp(b*rig);
  if ( a == 0. || minene == 0. || ave != ave ){
    //  if ( a == 0. || minene == 0. ){
    Float_t m = (maxene-minene)/(maxrig-minrig);
    Float_t q = minene - m * minrig;
    ave = rig * m + q;
  };
  //
  return(ave);
  //
}

Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){
  //
  Int_t therigb = 0;
  for (Int_t i = 0; i<nbin-2; i++){
    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
      therigb = i;
      break;
    };
  };
  //
  //
  if ( rig < brigm->At(0) ){
    if ( rig < brig->At(0) ){
//      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
    };
    therigb = 0;
  };
  if ( rig >= brigm->At(nbin-2) ){
    if ( rig >= brig->At(nbin-2) ) {
//      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
    };
    therigb = nbin - 3;
  };
  //
  return(this->GetFullAverageAt(plane,strip,rig,therigb));
  //
}

Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig, Int_t therigb){
  //
  Float_t minrig;
  Float_t maxrig;
  //
  minrig = brigm->At(therigb);
  maxrig = brigm->At(therigb+1);
  //
  Float_t minene = (*fqplmean[therigb])[plane][strip];
  Float_t maxene = (*fqplmean[therigb+1])[plane][strip];
  //
  if ( maxrig == minrig ){
   printf("Unrecoverable ERROR! Matrix will be screwed up... \n");
   return(0.);
  };
  Float_t b = log(maxene/minene)/(maxrig-minrig);
  Float_t a = minene/exp(minrig*b);
  Float_t ave = a*exp(b*rig);
  if ( a == 0. || minene == 0. || ave != ave ){
    Float_t m = (maxene-minene)/(maxrig-minrig);
    Float_t q = minene - m * minrig;
    ave = rig * m + q; 
  };
  //  
  //  ave += (44.-plane)*strip;
  //if ( a == 0. ) ave = 0.;
  if ( !(ave == ave) ) printf("a %f b %f ave %f maxene %f minene %f maxrig %f minrig %f \n",a,b,ave,maxene,minene,maxrig,minrig);
  //
  return(ave);
  //
}

Float_t CaloFranzini::GetHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
  Int_t therigb = 0;
  for (Int_t i = 0; i<nbin-2; i++){
    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
      therigb = i;
      break;
    };
  };
  //
  Float_t minrig;
  Float_t maxrig;
  //
  if ( rig < brigm->At(0) ){
    if ( rig < brig->At(0) ){
      //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
    };
    therigb = 0;
  };
  //  if ( rig >= brigm->At(nbin-4) ){ // -2
  if ( rig >= brigm->At(nbin-2) ){ 
    if ( rig >= brig->At(nbin-2) ) {
      //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
    };
    //    therigb = nbin - 5;// -3
    therigb = nbin - 3;
  };
  //
  if ( therigb < 2 ) therigb = 2;
  minrig = brigm->At(therigb);
  maxrig = brigm->At(therigb+1);
  //  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
  //
  Float_t minene = (*hmat[therigb])[iindex][jindex];
  Float_t maxene = (*hmat[therigb+1])[iindex][jindex];
  //  printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);
  //
//  Float_t a = 0.;
//  Float_t b = 0.;
//  Float_t ave = 0.;
//  if ( minene == 0. ){
//	
//  } else {
//  	b = log(maxene/minene)/(maxrig-minrig);
//  	a = minene/exp(minrig*b);
//  	ave = a*exp(b*rig);
//  };
  //
  Float_t m = (maxene-minene)/(maxrig-minrig);
  Float_t q = minene - m * minrig;
  Float_t ave = rig * m + q;
  if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);
  //
  //
  return(ave);
  //
}

Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){
  Int_t therigb = 0;
  for (Int_t i = 0; i<nbin-2; i++){
    if ( rig>=brigm->At(i) && rig < brigm->At(i+1) ){
      therigb = i;
      break;
    };
  };
  //
  if ( rig < brigm->At(0) ){
    if ( rig < brig->At(0) ){
      //      printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0));
    };
    therigb = 0;
  };
//  if ( rig >= brigm->At(nbin-4) ){ // -2
  if ( rig >= brigm->At(nbin-2) ){ 
    //if ( rig >= brigm->At(nbin-5) ){ 
    if ( rig >= brig->At(nbin-2) ) {
      //      printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2));
    };
    //    therigb = nbin - 5;// -3
    //  therigb = nbin - 5;// -3
    therigb = nbin - 3;
  };
  //
  if ( therigb < 2 ) therigb = 2;
  if ( therigb > 13 ) therigb = 13;
  //
  return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb));
  //
}


Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){
  //
  return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb,0));
  //
}

Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb, Int_t mtherig){  
  //
  //  Int_t lofit = 12;
  Int_t lofit = 3;
  //
  Float_t lowrig = 0.;
  Float_t minrig;
  Float_t maxrig;
  //
  if ( mtherig > lofit ) lowrig = brigm->At(therigb-1);
  minrig = brigm->At(therigb);
  maxrig = brigm->At(therigb+1);
  //if ( therigb > 10 )  printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex);
  //
  Float_t lowene = 0.;
  if ( mtherig > lofit ) lowene = (*hfmat[therigb-1])[iindex][jindex];
  Float_t minene = (*hfmat[therigb])[iindex][jindex];
  Float_t maxene = (*hfmat[therigb+1])[iindex][jindex];
  //  printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave);
  //
  Float_t ave = 0.;
  //
  //  Float_t a = 0.;
  //  Float_t b = 0.;
  //  Float_t ave = 0.;
  //  if ( minene == 0. ){
  //	
  //  } else {
  //  	b = log(maxene/minene)/(maxrig-minrig);
  //  	a = minene/exp(minrig*b);
  //  	ave = a*exp(b*rig);
  //  };
  //
  if ( mtherig > lofit ){
    //
    // FIT 
    //
//     Float_t x[3]={lowrig,minrig,maxrig};
//     Float_t y[3]={lowene,minene,maxene};
//     //
//     TGraph *tg= new TGraph(3,x,y);
//     //
//     //     gStyle->SetLabelSize(0.04);
//     //     gStyle->SetNdivisions(510,"XY");
//     //     TCanvas *c = new TCanvas();
//     //     c->Draw();
//     //     c->cd();
//     //     tg->Draw("AP");
//     //
//     TF1 *fun = new TF1("fun","pol2");
//     tg->Fit("fun","QNC");
//     //
//     ave = fun->GetParameter(0) + rig * fun->GetParameter(1) + rig * rig * fun->GetParameter(2);
//     //
//     //    printf(" therigb %i ave %f rig %f lowrig %f minrig %f maxrig %f lowene %f minene %f maxene %f \n",therigb,ave,rig,lowrig,minrig,maxrig,lowene,minene,maxene);
//     //
//     //
//     //     c->Modified();
//     //     c->Update();
//     //     gSystem->ProcessEvents();
//     //     gSystem->Sleep(6000);
//     //     c->Close();
//     //     gStyle->SetLabelSize(0);
//     //     gStyle->SetNdivisions(1,"XY");
//     //
//    tg->Delete();
//    fun->Delete();
    Float_t mrigl = (lowrig+minrig)/2.;
    Float_t mrigu = (minrig+maxrig)/2.;
    Float_t menel = (lowene+minene)/2.;
    Float_t meneu = (minene+maxene)/2.;
    Float_t m = (meneu-menel)/(mrigu-mrigl);
    Float_t q = menel - m * mrigl;
    ave = rig * m + q;
    //
  } else {
    Float_t m = (maxene-minene)/(maxrig-minrig);
    Float_t q = minene - m * minrig;
    ave = rig * m + q;
    if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave);
  };
  //
  //
  return(ave);
  //
}

void CaloFranzini::DrawLongAverage(Float_t rig){
  //
  TArrayF *ll = this->LoadLongAverage(rig);
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //
  TString hid = Form("cfralongvyvx");	
  TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
  if ( tcf ){
    tcf->cd();
  } else {
    tcf = new TCanvas(hid,hid);
  };
  //
  TString thid = Form("hfralongvyvx");	
  TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
  if ( thf ) thf->Delete();
  thf = new TH1F(thid,thid,44,-0.5,43.5);
  tcf->cd();
  Float_t qpl[43];
  for (Int_t st=0;st<43;st++){
    qpl[st] = ll->At(st);
    printf("st %i qpl %f\n",st,qpl[st]);
  };
  for (Int_t st=0;st<44;st++){
    if ( st == 37 ){
      thf->Fill(st,0.);
    } else {
      Int_t ss = st;
      if ( st > 37 ) ss--;
      thf->Fill(st,qpl[ss]);
    };
  };
  thf->Draw();
  tcf->Modified();
  tcf->Update();
  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};

void CaloFranzini::DrawLongAverage(Int_t bin){
  //
  TArrayF *ll = this->LoadLongAverage(brigm->At(bin));
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //
  TString hid = Form("cfralongvyvx");	
  TCanvas *tcf  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
  if ( tcf ){
    tcf->cd();
  } else {
    tcf = new TCanvas(hid,hid);
  };
  //
  TString thid = Form("hfralongvyvx");	
  TH1F *thf  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
  if ( thf ) thf->Delete();
  thf = new TH1F(thid,thid,44,-0.5,43.5);
  tcf->cd();
  Float_t qpl[43];
  for (Int_t st=0;st<43;st++){
    qpl[st] = ll->At(st);
    printf("st %i qpl %f\n",st,qpl[st]);
  };
  for (Int_t st=0;st<44;st++){
    if ( st == 37 ){
      thf->Fill(st,0.);
    } else {
      Int_t ss = st;
      if ( st > 37 ) ss--;
      thf->Fill(st,qpl[ss]);
    };
  };
  thf->Draw();
  tcf->Modified();
  tcf->Update();
  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};

Int_t CaloFranzini::ConvertStrip(Int_t mstrip){  
  //
  Int_t lastrip = 0;
//   //
//   if ( mstrip < 50 ) lastrip = 0;
//   //
//   if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1;
//   //
//   if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2;
//   //
//   if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3;
//   //
//   if ( mstrip >= 75 && mstrip < 84 ){    
//     lastrip = (int)trunc((mstrip - 75)/3.) + 4;
//   };
//   //
//   if ( mstrip >= 84 && mstrip < 90 ){    
//     lastrip = (int)trunc((mstrip - 84)/2.) + 7;
//   };
//   //
//   if ( mstrip >= 90 && mstrip < 101 ){    
//     lastrip = mstrip - 90 + 10;
//   };
//   //
//   if ( mstrip >= 101 && mstrip < 107 ){    
//     lastrip = (int)trunc((mstrip - 101)/2.) + 21;
//   };
//   //
//   //
//   if ( mstrip >= 107 && mstrip < 116 ){    
//     lastrip = (int)trunc((mstrip - 107)/3.) + 24;
//   };
//   //
//   if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27;
//   //
//   if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28;
//   //
//   if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29;
//   //
//   if ( mstrip >= 141 ) lastrip = 30;
//   //
  //
  if ( mstrip < 83 ) lastrip = 0;
  //
  if ( mstrip >= 83 && mstrip < 90 ) lastrip = 1;
  //
  if ( mstrip >= 90 && mstrip < 93 ) lastrip = 2;
  //
  if ( mstrip >= 93 && mstrip < 98 ){    
    lastrip = mstrip - 93 + 3;
  };
  //
  if ( mstrip >= 98 && mstrip < 101 ) lastrip = 8;
  //
  if ( mstrip >= 101 && mstrip < 107 ) lastrip = 9;
  //
  if ( mstrip >= 107 ) lastrip = 10;
  //
  return(lastrip);
}