/////////////////////////////////////////////////////////////////////////////
//
// EEezClusterQA
//
// Author: Jason C. Webb <jwebb@iucf.indiana.edu>
//
// EEezClusterQA provides a set of QA histograms which may be filled
// using a single call to EEezClusterQA::Fill().  The arguement of fill
// may either be a pointer to a cluster or a vector of cluster pointers.
// 
// The code assumes that the user has set a pointer to the EEsmdProfile
// object, via setSmdProfiler().  An instance of EEezAnalysis should
// also be set, via setAnalysis(), but is not required.
//
/////////////////////////////////////////////////////////////////////////////
#include "EEezClusterQA.h"

#include <TH1F.h>
#include <TH2F.h>
#include <TProfile.h>
#include <TVector3.h>

#include "StEEmcUtil/EEmcGeom/EEmcGeomDefs.h"
#include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
#ifndef __ROOT__
#include "EEezSectorStats.h"
#endif

#include "EEezTower.h"
#include "EEezStrip.h"
#include "EEezCluster.h"

#include <iostream>

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

ClassImp(EEezClusterQA);

 EEezClusterQA::EEezClusterQA( const Char_t *name, const Char_t *title ) 
  : TDirectory( name, title )
{
  // Class constructor

  m_EEezAnaly = 0;
  m_EEsmdProf = 0;
  m_EEsmdGeom = EEmcSmdGeom::instance();

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEezClusterQA::Init() 
{
  // Initialize the class

  bookHistograms();

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEezClusterQA::bookHistograms() 
{
  // Book all QA histograms
  
  cd();
  m_XY           = new TH2F("m_XY",         "Y vs X of clusters",100,-175.,175.,100,-175.,175.);
  m_XY_smd       = new TH2F("m_XY_smd",     "Y vs X of clusters, from SMD response",100,-175.,175.,100,-175.,175.);
  m_Energy       = new TH1F("m_Energy",     "Energy of clusters",100,0.,20.);
  m_EcVsEtot     = new TH2F("m_EcVsEtot",   "Cluster energy vs sector energy sum", 100, 0., 20., 100, 0., 20. );
  m_HitStrips[0] = new TH1F("m_HitStripsU", "Frequency plot of hit U strips vs strip index (id-1)",288,0.,288.);
  m_HitStrips[1] = new TH1F("m_HitStripsV", "Frequency plot of hit V strips vs strip index (id-1)",288,0.,288.);
  m_MipStrips[0] = new TH1F("m_MipStripsU", "Frequency plot of mip U strips vs strip index (id-1)",288,0.,288.);
  m_MipStrips[1] = new TH1F("m_MipStripsV", "Frequency plot of mip V strips vs strip index (id-1)",288,0.,288.);
  m_Sum9VsEc[0]  = new TH2F("m_Sum9VsEcU",  "Sum of 9 strips near SMD-U centroid vs E_{cluster}",50,0.,10.,100,0.,100.);
  m_Sum9VsEc[1]  = new TH2F("m_Sum9VsEcV",  "Sum of 9 strips near SMD-V centroid vs E_{cluster}",50,0.,10.,100,0.,100.);
  m_Sum5VsEc[0]  = new TH2F("m_Sum5VsEcU",  "Sum of 9 strips near SMD-U centroid vs E_{cluster}",50,0.,10.,100,0.,100.);
  m_Sum5VsEc[1]  = new TH2F("m_Sum5VsEcV",  "Sum of 9 strips near SMD-V centroid vs E_{cluster}",50,0.,10.,100,0.,100.);

  m_SumVsEc[0]   = new TH2F("m_SumVsEcU",   "Sum of SMD-U strips under seed tower vs E_{cluster}",50,0.,10.,100,0.,100.);
  m_SumVsEc[1]   = new TH2F("m_SumVsEcV",   "Sum of SMD-V strips under seed tower vs E_{cluster}",50,0.,10.,100,0.,100.);
  m_SumVsEc_pfx[0] = new TProfile("m_SumVsEcU_pfx", "Mean sum of SMD-U strips under seed tower vs E_{cluster}",50,0.,10.,0.,20000.);
  m_SumVsEc_pfx[1] = new TProfile("m_SumVsEcV_pfx", "Mean sum of SMD-V strips under seed tower vs E_{cluster}",50,0.,10.,0.,20000.);
   

  m_Towers       = new TH1F("m_Towers",     "Tower energy response vs index", 720, 0., 720.);

  m_Sum5         = new TH2F("m_Sum5","Sum of 5 strips U vs V",25,0.,100.,25,0.,100.);
  m_Sum9         = new TH2F("m_Sum9","Sum of 9 strips U vs V",25,0.,100.,25,0.,100.);
  m_Frac5to9[0]  = new TH1F("m_Frac5to9U","Ratio of sum of 5 strips to sum of 9 strips",100,0.,1.);
  m_Frac5to9[1]  = new TH1F("m_Frac5to9V","Ratio of sum of 5 strips to sum of 9 strips",100,0.,1.);
 
  m_Yield        = new TH2F("m_Yield","Fit yield U vs V",50,0.,200.,50,0.,200.);
  m_Sigma        = new TH2F("m_Sigma","Fit sigma U vs V",100,0.,5.,100,0.,5.);
  m_RedChi2      = new TH2F("m_RedChi2","#chi^{2}/#nu U vs V",50,0.,20.,50,0.,20.);
  m_YieldVsEc[0] = new TH2F("m_YieldVsEcU","Fit yield vs E_{cluster}",100,0.,10.,100,0.,200.);
  m_YieldVsEc[1] = new TH2F("m_YieldVsEcV","Fit yield vs E_{cluster}",100,0.,10.,100,0.,200.);
  m_YieldToSum[0] = new TH1F("m_YieldToSumU","Ratio of fit yield to sum over seed tower, U",300,0.,3.);
  m_YieldToSum[1] = new TH1F("m_YieldToSumV","Ratio of fit yield to sum over seed tower, V",300,0.,3.);

  
  m_FreqBigChi2[0] = new TH1F("m_FreqBigChi2U","Frequency plot of strip with highest #chi^{2}",288,0.,288.);
  m_FreqBigChi2[1] = new TH1F("m_FreqBigChi2V","Frequency plot of strip with highest #chi^{2}",288,0.,288.);

  m_FreqEmpty[0]   = new TH1F("m_FreqEmptyU",  "Frequency with which strips show no response near fit, and yield>50 mips",288,0.,288.);
  m_FreqEmpty[1]   = new TH1F("m_FreqEmptyV",  "Frequency with which strips show no response near fit, and yield>50 mips",288,0.,288.);


  m_SeedAdc[0] = new TH1F("m_SeedTowerAdc", "Seed-tower ADC spectrum", 1024, 0., 1024. );
  m_SeedAdc[1] = new TH1F("m_SeedPre1Adc",  "Seed-preshower-1 ADC spectrum",1024,0.,1024. );
  m_SeedAdc[2] = new TH1F("m_SeedPre2Adc",  "Seed-preshower-2 ADC spectrum",1024,0.,1024. );
  m_SeedAdc[3] = new TH1F("m_SeedPostAdc",  "Seed-postshower  ADC spectrum",1024,0.,1024. );

  m_NeighborAdc[0] = new TH1F("m_NeighborTowerAdc", "Neighbor-tower ADC spectrum", 1024, 0., 2048. );
  m_NeighborAdc[1] = new TH1F("m_NeighborPre1Adc",  "Neighbor-preshower-1 ADC spectrum",1024,0.,2048. );
  m_NeighborAdc[2] = new TH1F("m_NeighborPre2Adc",  "Neighbor-preshower-2 ADC spectrum",1024,0.,2048. );
  m_NeighborAdc[3] = new TH1F("m_NeighborPostAdc",  "Neighbor-postshower  ADC spectrum",1024,0.,2048. );

  m_SumSmd = new TH2F("m_SumSmd","Sum of mips beneath seed tower, U vs V",100,0.,400.,100,0.,400.);

  //-- Fit parameters to SMD distributions
  const Float_t hmin[] = { 0.,0.,0.,0.,0. };
  const Float_t hmax[] = { 500., 288., 5., 1., 5.0 };
  for ( Int_t abc = 0; abc < 5; abc++ ) {
    TString hname = "m_Parameters["; hname += abc; hname += "]";
    TString htitle = "Fit parameter "; htitle += abc;
    m_Parameters[abc] = new TH1F(hname,htitle,100,hmin[abc],hmax[abc]);
  }


  cd(); m_UFail = new TDirectory("UFail","SMD-U profiles which failed a user cut");
  cd(); m_VFail = new TDirectory("VFail","SMD-V profiles which passed a user cut");

  cd(); m_UPass = new TDirectory("UPass","SMD-U profiles which passed a user cut");
  cd(); m_VPass = new TDirectory("VPass","SMD-V profiles which passed a user cut");

  cd(); m_UFlag = new TDirectory("UFlag","SMD-U profiles flagged for study");
  cd(); m_VFlag = new TDirectory("VFlag","SMD-V profiles flagged for study");

  cd("..");
  
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEezClusterQA::Fill ( EEezCluster *cluster ) 
{
  // Given a pointer to a cluster, fill the QA histograms

  assert(cluster);

  //-- Get the seed tower for the cluster
  EEezTower *seed = cluster -> getSeedTower();
 

  m_SumSmd -> Fill ( seed -> getSumMipV(), seed -> getSumMipU() );

  //-- Unit vector pointing towards cluster centroid
  TVector3 unit = ( cluster -> getMomentum() ).Unit();
  m_XY     -> Fill ( kEEmcZSMD * unit.X(), kEEmcZSMD * unit.Y() );

  Float_t umean = seed -> getMeanU();
  Float_t vmean = seed -> getMeanV();

  if ( umean > 0. && vmean > 0. ) {
    TVector3 post = m_EEsmdGeom -> getIntersection( seed -> getSector(), (Int_t)umean, (Int_t)vmean );
    m_XY_smd -> Fill ( kEEmcZSMD * post.Unit().X(), 
		       kEEmcZSMD * post.Unit().Y() );
  }
  
  //-- Energy of the cluster
  Float_t eCluster = cluster -> getEnergy();
  m_Energy -> Fill ( eCluster );
  
  //-- Sector-based sums (only fill if m_EEezAnaly is not NULL)  
  //Int_t sector = seed -> getSector();
  
  //----- Hit strips vs index -----

  for ( Int_t i = 0; i < seed -> getNUStrips(); i++ ) {
    EEezStrip *strip = seed -> getUStrip(i);
    if ( strip -> getEnergy() > 0. ) {
      m_HitStrips[0] -> Fill( strip -> getIndex() + 0.5 );
      m_MipStrips[0] -> Fill( strip -> getIndex() + 0.5,
			      strip -> getNumMips() );
    }
  }
  for ( Int_t i = 0; i < seed -> getNVStrips(); i++ ) {
    EEezStrip *strip = seed -> getVStrip(i);
    if ( strip -> getEnergy() > 0. ) {
      m_HitStrips[1] -> Fill( strip -> getIndex() + 0.5 );
      m_MipStrips[1] -> Fill( strip -> getIndex() + 0.5, 
			      strip -> getNumMips() );
    }
  }



  //-----<<<<< Kludge >>>>>-----

  assert(m_EEsmdProf);


  //-- Sums over 5 and 9 strips vs cluster energy
  TH1F *hsum5U = m_EEsmdProf -> getProfileSum5U(cluster);
  TH1F *hsum5V = m_EEsmdProf -> getProfileSum5V(cluster);
  TH1F *hsum9U = m_EEsmdProf -> getProfileSum9U(cluster);
  TH1F *hsum9V = m_EEsmdProf -> getProfileSum9V(cluster);

  //-- Sum over 5 strips (but trap NULL)
  Float_t sum5U = (hsum5U!=NULL) ? hsum5U -> GetMaximum() : -1.;
  Float_t sum5V = (hsum5V!=NULL) ? hsum5V -> GetMaximum() : -1.;
  Float_t sum9U = (hsum9U!=NULL) ? hsum9U -> GetMaximum() : -1.;
  Float_t sum9V = (hsum9V!=NULL) ? hsum9V -> GetMaximum() : -1.;

  m_Sum5VsEc[0] -> Fill( eCluster, sum5U );
  m_Sum5VsEc[1] -> Fill( eCluster, sum5V );
  m_Sum9VsEc[0] -> Fill( eCluster, sum9U );
  m_Sum9VsEc[1] -> Fill( eCluster, sum9V );

  m_Sum5 -> Fill ( sum5V, sum5U );
  m_Sum9 -> Fill ( sum9V, sum9U );

  m_SumVsEc[0] -> Fill ( eCluster, seed -> getSumMipU() );
  m_SumVsEc[1] -> Fill ( eCluster, seed -> getSumMipV() );
  m_SumVsEc_pfx[0] -> Fill ( eCluster, seed -> getSumMipU() );
  m_SumVsEc_pfx[1] -> Fill ( eCluster, seed -> getSumMipV() );

  if ( sum9U > 0. ) {
    m_Frac5to9[0] -> Fill ( sum5U/sum9U );
  }
  if ( sum9V > 0. ) {
    m_Frac5to9[1] -> Fill ( sum5V/sum9V );
  }

  //-- Fit parameters
  TF1 *fitU = m_EEsmdProf -> getFitU(cluster);
  TF1 *fitV = m_EEsmdProf -> getFitV(cluster);

  if ( (fitU != NULL) && (fitV != NULL) ) {
    
    m_Yield -> Fill ( fitV -> GetParameter(0), fitU -> GetParameter(0) );
    m_Sigma -> Fill ( fitV -> GetParameter(2), fitU -> GetParameter(2) );

    Float_t chi2u = fitU -> GetChisquare();
    Float_t chi2v = fitV -> GetChisquare();
    Float_t ndfu  = fitU -> GetNDF();
    Float_t ndfv  = fitV -> GetNDF();

    if ( ndfu > 0. && ndfv > 0. ) {
      m_RedChi2 -> Fill ( chi2v/ndfv, chi2u/ndfu );
    }

    m_YieldVsEc[0] -> Fill ( eCluster, fitU -> GetParameter(0) );
    m_YieldVsEc[1] -> Fill ( eCluster, fitV -> GetParameter(0) );

    for ( Int_t abc = 0; abc < 5; abc++ ) {
      m_Parameters[abc] -> Fill ( fitU -> GetParameter(abc) );
      m_Parameters[abc] -> Fill ( fitV -> GetParameter(abc) );
    }

  }

  //-- Frequency plots for U and V strip having largest chi2 

  for ( Int_t i = 0; i < 2; i++ ) {
    TH1F *chi2strip = (i==0)?
      m_EEsmdProf -> getProfileChi2U(cluster):
      m_EEsmdProf -> getProfileChi2V(cluster);
    TH1F *prof = (i==0)?
      m_EEsmdProf -> getProfileU(cluster):
      m_EEsmdProf -> getProfileV(cluster);
    Float_t yy = prof -> Integral();
    if ( !chi2strip ) continue; //------<<<<<< CHECK WITH ASSERT! >>>>>>------
    Float_t x = chi2strip -> GetBinCenter( chi2strip -> GetMaximumBin() ); 
    m_FreqBigChi2[i] -> Fill(x);

    for ( Int_t j = 1; j <= chi2strip->GetNbinsX(); j++ ) {
      Float_t xx = chi2strip->GetBinCenter(j);
      if ( chi2strip->GetBinContent(j)==0 &&yy>50.0) 
	m_FreqEmpty[i] -> Fill(xx);
    }

  }
 
  // ------- tower, preshower, postshower ADC response -------
  
  for ( Int_t i = 0; i < 4; i++ ) m_SeedAdc[i] -> Fill ( seed -> getAdc(i) );

  EEezTowerPtrVec_t neighbors = seed -> getNeighbors();
  EEezTowerPtrVecIter_t iter = neighbors.begin();
  while ( iter != neighbors.end() ) {
    for ( Int_t i = 0; i < 4; i++ ) 
      m_NeighborAdc[i] -> Fill ( (*iter)->getAdc(i) );
    iter++;
  }


}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEezClusterQA::Fill ( EEezClusterPtrVec_t *clusters ) 
{
  // Loop over all clusters in the vector, filling histograms
  // from the clusters contained therein.

  EEezClusterPtrVecIter_t iter = clusters -> begin();
  while ( iter != clusters -> end() ) {

    //-- Get a pointer to the cluster, for simplicity
    EEezCluster *c = (*iter);
   
    Fill(c);
   
    iter++;
  }

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEezClusterQA::pass( EEezCluster *cluster, const Char_t *tag )
{
  // Save the SMD profile and associated histograms in a directory
  // as an example of a cluster which passed a cut.  One should not
  // call this on every event, just a subset of the data.

  static Int_t npass = 0;
    if ( npass++ >= 100 ) return;
  
  m_UPass -> cd();
  save( cluster, 0, tag );
  m_UPass -> cd("..");

  m_VPass -> cd();
  save( cluster, 1, tag );
  m_VPass -> cd("..");

}

 void EEezClusterQA::fail( EEezCluster *cluster, const Char_t *tag )
{
  // Save the SMD profile and associated histograms in a directory
  // as an example of a cluster which failed a cut.  One should not
  // call this on every event, just a subset of the data.

  static Int_t npass = 0;
    if ( npass++ >= 100 ) return;

  m_UFail -> cd();
  save( cluster, 0, tag );
  m_UFail -> cd("..");

  m_VFail -> cd();
  save( cluster, 1, tag );
  m_VFail -> cd("..");

}

 void EEezClusterQA::flag( EEezCluster *cluster, const Char_t *tag )
{
  // Save the specified cluster in the "flag" directory, just as 
  // with pass and fail.  Here we do not limit ourselves to the 
  // first 50 clusters.  Use with care.

  m_UFlag -> cd();
  save ( cluster, 0, tag );
  m_UFlag -> cd("..");
  
  m_VFlag -> cd();
  save ( cluster, 1, tag );
  m_VFlag -> cd("..");

}

 void EEezClusterQA::save( EEezCluster *cluster, Int_t iuv, const Char_t *tag )
{
  // Saves the specified cluster's SMD profile (and other assoicated
  // histograms) in the current directory.

  //-- Get the specified profile histograms

  TH1F *profile = (iuv)?m_EEsmdProf->getProfileV(cluster):m_EEsmdProf->getProfileU(cluster);
  TH1F *profsum = (iuv)?m_EEsmdProf->getProfileSum5V(cluster):m_EEsmdProf->getProfileSum5U(cluster);
  TH1F *fitinit = (iuv)?m_EEsmdProf->getProfileInitV(cluster):m_EEsmdProf->getProfileInitU(cluster);

  //-- Clone them (create new ones in current directory)
  TString profile_name;

  if ( profile != NULL ) {

    profile_name = profile -> GetName();
    
    profile_name.ReplaceAll("hu","smd_u_prof_");
    profile_name.ReplaceAll("hv","smd_v_prof_");

    if ( tag != "" ) {
      profile_name += tag;
      profile_name += tag; 
    }

    profile -> Clone( profile_name );
  
    if ( profsum != NULL ) {
      profile_name.ReplaceAll("prof","sum5");
      profsum -> Clone( profile_name );
    }
    if ( fitinit != NULL ) {
      profile_name.ReplaceAll("sum5","init");
      profile_name.ReplaceAll("prof","init");
      fitinit -> Clone( profile_name );
    }

  }

}

 void EEezClusterQA::save( TH1F *histo, Int_t iuv, const Char_t *tag )
{
  // Saves the specified histogram
  cd();
  histo -> Clone();
  cd("..");

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEezClusterQA::bookCut ( const Char_t *name, const Char_t *title, Int_t nbin, Float_t min, Float_t max )
{
  // Book a histogram for quantities which are cut on.   Histograms are
  // stored in a TList, which is searched whenever we fill.  This is
  // slow, and should be changed to a hash table.

  if ( !m_ListOfCuts ) 
    m_ListOfCuts = new TList();

  cd();
  m_ListOfCuts -> Add ( new TH1F(name,title,nbin,min,max) );
  cd("..");

}

 void EEezClusterQA::fillCut( const Char_t *name, Float_t value )
{
  // Fills the specified cut
  TH1F *h = ((TH1F*)m_ListOfCuts -> FindObject( name )) ;
  if ( h != NULL ) 
    h -> Fill ( value );
  else
    Warning("fillCut","%s not found",name);

}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////

 void EEezClusterQA::bookCut ( const Char_t *name, const Char_t *title, Int_t nbin, Float_t min, Float_t max, Int_t nbiny, Float_t miny, Float_t maxy )
{
  // Book a histogram for quantities which are cut on.   Histograms are
  // stored in a TList, which is searched whenever we fill.  This is
  // slow, and should be changed to a hash table.

  if ( !m_ListOf2DCuts ) 
    m_ListOf2DCuts = new TList();

  cd();
  m_ListOf2DCuts -> Add ( new TH2F(name,title,nbin,min,max,nbiny,miny,maxy) );
  cd("..");

}

 void EEezClusterQA::fillCut( const Char_t *name, Float_t value, Float_t valuey )
{
  // Fills the specified cut
  //  ((TH1F*)m_ListOf2DCuts -> FindObject( name )) -> Fill(value,valuey);

  TH2F *h = ((TH2F*)m_ListOf2DCuts -> FindObject( name )) ;
  if ( h != NULL ) 
    h -> Fill ( value, valuey );
  else
    Warning("fillCut","%s not found",name);
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////