/////////////////////////////////////////////////////////////////////////////
//                                                                         
// EEsmdProfile                                                            
//
// Author: Jason C. Webb <jwebb@iucf.indiana.edu>                                                                      //
//
// EEsmdProfile creates histograms of the SMD distributions which are 
// beneath the seed towers of each cluster.
//
// Example of how to improve cluster position with EEsmdProfile.  This
// example neglects several things, including the overlap of
// strips in neighboring sectors and other edge effects, and is not
// meant for production.
//
// // Will need the smd geometry class
// EEmcSmdGeom *geom = EEmcSmdGeom::instance();
//
// // We have a cluster finder from somewhere 
// EEezClAnalysis *clust = ...
//
// // The profile analysis
// EEsmdProfile *prof = new EEsmdProfile();
//
// ... code which runs the cluster finder ...
//
// // Create the SMD profiles for each cluster
// prof -> Make ( clust );
//
// // Improve the first cluster's position.  First, obtain a pointer
// // to the cluster.
// EEezCluster *cluster = clust -> getClusterPtr(0);
//
// // We will need to know what sector we are in
// Int_t sector = cluster -> getSeedTower() -> getSector();
//
// // Next, obtain the SMD histograms corresponding to the cluster
// TH1F *uHisto = prof -> getProfileU(0);
// TH1F *vHisto = prof -> getProfileV(0);
//
// // Somehow obtain the mean of these histograms... Fit to the correct
// // shower shape, once we know it...
// Float_t uMean = ...
// Float_t vMean = ...
//
// // Now, find the intersection of these two strips within the
// // sector of the seed tower using the geometry class
// TVector3 intersection = geom -> getIntersection(sector,(Int_t)uMean,(Int_t)vMean);
//
// // Make intersection a unit vector
// intersection = intersection.Unit();
//
// // Recalculate cluster momentum
// TVector3 momentum = intersection * cluster -> getEnergy();
//
// // And set its momentum
// cluster -> setMomentum( momentum );
//
// **** NOTE **** **** NOTE **** **** NOTE **** **** NOTE **** **** NOTE **** 
//
// Will need to think about whether the fit to the centroid of this 
// distribution will obtain the correct position/index or not... since
// the strips are indexed from 0... i.e. if strip number 64 (index 63)
// is fit to a gaussian... I bet the mean comes out to be 63.5... so the
// fits to these distributions should probably be incremented by 0.5,
// before being passed to EEmcSmdGeom::getIntersection().
//
// **** NOTE **** **** NOTE **** **** NOTE **** **** NOTE **** **** NOTE **** 
//
// To Do List:
//
// 1. Need to handle/flag empty histograms
// 
//
/////////////////////////////////////////////////////////////////////////////

#include "EEsmdProfile.h"
#include "EEezClAnalysis.h"
#include "EEezCluster.h"
#include "EEezTower.h"
#include "EEezStrip.h"
#include "TH1F.h"
#include "TMap.h"
#ifndef __ROOT__
#include "EEezClusterQA.h"
#endif

#include <iostream>

// Shower-shape functions
#include "eeSinglePeak.h"
#include "eeDoublePeak.h"

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
ClassImp(EEsmdProfile);

 EEsmdProfile::EEsmdProfile() 
{
  // Default Constructor (not meant to be called, use the name/title version).

  m_MinEnergy = 0.;
  m_MaxEnergy = 999999.;

  m_ShowerShape = NULL;

}

 EEsmdProfile::EEsmdProfile( const Char_t *name, const Char_t *title )
  :TDirectory(name,title) 
{
  // Constructor w/ name and title
  EEsmdProfile();
}

/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
 Int_t EEsmdProfile::Init() 
{
  // Initialization
  if ( !m_ShowerShape ) 
    m_ShowerShape = new TF1("gauss","gaus");
  
  m_SinglePeak = new TF1("eeSinglePeak", eeSinglePeak, 0, 288, 5 );
  m_SinglePeak -> SetParameter(4,3.5);
  m_SinglePeak -> SetParLimits(4,3.5,3.5);

  m_SinglePeak -> SetLineColor(2);
  m_SinglePeak -> SetLineStyle(2);
  m_SinglePeak -> SetLineWidth(1);

  m_DoublePeak = new TF1("eeDoublePeak", eeDoublePeak, 0, 288, 10 );
  m_DoublePeak -> SetParameter(4,3.5);
  m_DoublePeak -> SetParLimits(4,3.5,3.5);
  m_DoublePeak -> SetParameter(9,3.5);
  m_DoublePeak -> SetParLimits(9,3.5,3.5);

  m_DoublePeak -> SetLineColor(2);
  m_DoublePeak -> SetLineStyle(2);
  m_DoublePeak -> SetLineWidth(1);

  cd();

  m_Converge = new TH1F("m_Converge","Number of iterations for single-peak fit to stabilize",10,0.,10.);

#ifndef __ROOT__
  for ( Int_t i = 0; i < 12; i++ ) {

    Int_t sec = i+1;
    TString name = "QAsec"; if ( sec < 10 ) name += "0"; name += sec;
    m_EEClusterQA[i] = new EEezClusterQA(name,"EEMC Cluster QA / EEsmdProfile");
    m_EEClusterQA[i] -> setSmdProfiler( this );
    m_EEClusterQA[i] -> Init();
    m_EEClusterQA[i] -> setSmdProfiler( this );

  }
#endif

  cd("..");

  return 1;

}

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

 Int_t EEsmdProfile::Make( EEezClAnalysis *analy, 
			  EEmcEventHeader *header 
			  ) 
{
  // Takes as input a pointer to the cluster finder and to an
  // EEmcEventHeader.  Loops over all clusters and creates a
  // histogram of the energy response of SMD strips versus
  // strip number.  The histograms ordering of the histograms
  // matches the ordering of the clusters.

  m_EEezClust = analy;
  

  // Loop over all clusters found in the cluster analysis
  Int_t nc = analy -> getNClusters();

  //  std::cout << "Preparing to loop over " << nc << " clusters" << std::endl << std::flush;

  for ( Int_t ic = 0; ic < nc; ic++ ) { 
    
    // Get the cluster from the analysis
    EEezCluster *cluster = analy -> getClusterPtr( ic );

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

    // As a starting point, we only fill a histogram
    // over the range of strips associated with the
    // seed tower.  In the future, we should probably
    // expand this to all towers in the cluster.  Of 
    // course, sector boundaries present complications.
    // As do overlap strips.
    
    // Names and titles for histograms

    TString hname="hu";
    hname += seed-> getName();
    hname += ":";
    if ( header ) 
      hname += header -> getEventNumber();
    else
      hname += "1";
    
    TString htitle="SMD u-strip distribution, event number ";
    if ( header )
      htitle += header -> getEventNumber();
    else
      htitle += "1";
    htitle += ", seed tower ";
    htitle += seed -> getName();

    TString vname="hv";
    vname += seed-> getName();
    vname += ":";
    if ( header )
      vname += header -> getEventNumber();
    else
      vname += "1";
    
    TString vtitle="SMD v-strip distributions, event number ";
    if ( header )
      vtitle += header -> getEventNumber();
    else
      vtitle += "1";
    vtitle += ", seed tower ";
    vtitle += seed -> getName();


    // Get the appropriate range of U and V strips
    EEezStripPtrVec_t *ustrips = seed -> getUStrips();
    EEezStripPtrVec_t *vstrips = seed -> getVStrips();
    
    // Get the min and max U and V strips
    EEezStripPtrVecIter_t iter;
    iter = ustrips -> begin();
    EEezStrip *ufirst = *iter;
    iter = ustrips -> end() - 1;
    EEezStrip *ulast  = *iter;
    iter = vstrips -> begin();
    EEezStrip *vfirst = *iter;
    iter = vstrips -> end() - 1;
    EEezStrip *vlast  = *iter;

    Int_t umin = ufirst -> getIndex() - 2;
    Int_t vmin = vfirst -> getIndex() - 2;
    Int_t umax = ulast  -> getIndex() + 2;
    Int_t vmax = vlast  -> getIndex() + 2;
    

    // Histograms for storage, created on the heap, so
    // we need to erase them when we clear.
    TH1F *hu = new TH1F(hname,htitle,(umax-umin),umin,umax);
    TH1F *hv = new TH1F(vname,vtitle,(vmax-vmin),vmin,vmax);

    TString xtitle = "strip number";
    TString ytitle = "N_{mips} = (E_{strip}/E_{mip}) = ( ADC - ped ) / gain[mip/ch]";

    hu -> GetXaxis() -> SetTitle(xtitle);
    hv -> GetXaxis() -> SetTitle(xtitle);
    hu -> GetYaxis() -> SetTitle(ytitle);
    hv -> GetYaxis() -> SetTitle(ytitle);


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

    // Fill said histograms with energy deposit 
    // versus strip number (index)

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

      //Float_t energy = (*iter)->getEnergy();
      Int_t   index  = (*iter)->getIndex();
      Float_t npe    = (*iter)->getNumPhotoElectrons();
      Float_t nmips  = (*iter)->getNumMips();
      
      if ( nmips > 0. ) {

// 	Float_t emips  = 0.;
// 	if ( npe > 0. ) { 
// 	  emips = nmips / TMath::Sqrt(npe);
// 	}


	Float_t emips = 0.;
	
	//	Float_t e_npe  = TMath::Sqrt(npe);
	//	Float_t e_nmip = TMath::Sqrt(nmips);
	Float_t e_tot  = TMath::Sqrt(npe+nmips);
	if ( e_tot > 0. ) {
	  emips = nmips / e_tot;
	}


	hu -> Fill( index, nmips );
	hu -> SetBinError( index - umin + 1, emips );

      }
          
      iter++;

    }

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

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

      //Float_t energy = (*iter)->getEnergy();
      Int_t   index  = (*iter)->getIndex();
      Float_t npe    = (*iter)->getNumPhotoElectrons();
      Float_t nmips  = (*iter)->getNumMips();
      
      if ( nmips > 0. ) {

	Float_t emips  = 0.;
	if ( npe > 0. ) { 
	  emips = nmips / TMath::Sqrt(npe);
	}

	hv -> Fill( index, nmips );
	hv -> SetBinError( index - vmin + 1, emips );

      }
          
      iter++;

    }

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

    // Store the histograms

    m_UProfile.AddLast(hu);
    m_VProfile.AddLast(hv);

    // Store the index of the histogram, associated with
    // the name (see EEezCluster::Hash()) of the 
    m_Name2IndexU[ cluster -> Hash() ] = m_UProfile.LastIndex();
    m_Name2IndexV[ cluster -> Hash() ] = m_VProfile.LastIndex();
    
  }

  return 1;

}

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

 void EEsmdProfile::MakeFits( Option_t *opts ) 
{
  // Loop over all histograms and fit them to a single "eeSinglePeak()"
  // shower shape (double gaussian with common mean, constrained widths).


  // Loop over all profiles (aka clustes)
  for ( Int_t i = 0; i < m_UProfile.GetSize(); i++ ) {


    // Retrieve the histograms
    TH1F *uprofile = getProfileU(i);
    TH1F *vprofile = getProfileV(i);

    assert(uprofile != NULL);
    assert(vprofile != NULL);

    // Find the highest bin in both views
    Int_t umax_bin = uprofile -> GetMaximumBin();
    Int_t vmax_bin = vprofile -> GetMaximumBin();

    
    //-- 04/10/04 The "mean" parameter in fitSinglePeak has been
    //-- disabled.  The peak will now be auto-magically located
    
    //-- fitDoublePeak will need to be updated.


    // Convert this to a strip number...
    Float_t umax = (Float_t)umax_bin + uprofile -> GetBinLowEdge(0);
    Float_t vmax = (Float_t)vmax_bin + vprofile -> GetBinLowEdge(0);
    
    // Fit the two histograms to a single peak

    fitSinglePeak(i,0,umax,0.,-1.,opts);
    fitSinglePeak(i,1,vmax,0.,-1.,opts);

  }

}


#ifndef __ROOT__
 void EEsmdProfile::MakeQA() 
{
  // Fills QA histograms

  for ( Int_t i = 0; i < m_EEezClust -> getNClusters(); i++ ) {    

    // Fill cluster QA histograms
    EEezCluster *cluster = m_EEezClust -> getClusterPtr(i);
    Int_t sec = cluster -> getSeedTower() -> getSector();

    //-- Somehow "this" gets lost between events!?
    m_EEClusterQA[sec] -> setSmdProfiler( this );
    m_EEClusterQA[sec] -> Fill ( cluster );
    //-- All pass for now
    m_EEClusterQA[sec] -> pass( cluster );


  }

}
#endif

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

 void EEsmdProfile::Clear( Option_t *opts ) 
{  
  // Clear the list of profiles and delete the histograms

  m_UProfile.Delete();
  m_VProfile.Delete();

  m_UProfileSum9.Delete();
  m_VProfileSum9.Delete();
  
  m_UProfileSum5.Delete();
  m_VProfileSum5.Delete();

  m_UProfileChi2.Delete();
  m_VProfileChi2.Delete();

  m_UProfileInit.Delete();
  m_VProfileInit.Delete();

  m_UResidual.Delete();
  m_VResidual.Delete();

  m_Name2IndexU.clear();
  m_Name2IndexV.clear();

}

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

 void EEsmdProfile::fitProfile ( Int_t iprofile, Int_t iuv ) 
{ 
  // Fits the SMD profile ( U = 0, V = 1 ) to the functional form
  // specified in setShowerShape().  If none was specified, the 
  // method defaults to a gaussian shower shape.
  
  
  // Obtain a pointer to the requested profile
  assert( iuv>=0 && iuv <=1 ); // Only two planes
  TH1F *profile = ( iuv == 0 ) ? getProfileU(iprofile) : getProfileV(iprofile);

  // First obtain some info from the histogram to
  // come up with initial parameters
  Float_t mean = profile -> GetMean();
  Float_t rms  = profile -> GetRMS();
  Float_t sum  = profile -> Integral();

  // Next, get the bin number of the highest-responding
  // strip and the its corresponding position

  //  Int_t maxbin = profile -> GetMaximumBin();
  //Float_t xmax = profile -> GetBinLowEdge(0) + maxbin + 0.5;

  // Initial guess for gaussian parameters
  Float_t a = sum / TMath::Sqrt( TMath::TwoPi() ) / rms;
  TF1 *fit = m_ShowerShape;

  fit -> SetParameter(0,a);
  fit -> SetParameter(1,mean);
  fit -> SetParameter(2,rms);

  fit -> SetRange(mean-5.0,mean+5.0);

  // Perform the fit
  profile -> Fit(fit,"RNQ");

  // Next, increase the range to 3 sigma and refit
  Float_t min = fit -> GetParameter(1) - 3.0 * fit -> GetParameter(2);
  Float_t max = fit -> GetParameter(1) + 3.0 * fit -> GetParameter(2);

  profile -> GetXaxis() -> SetRangeUser(min,max);
  profile -> Fit(fit,"RQ");
  profile -> GetXaxis() -> SetRange(0,0);


  
}

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


 void EEsmdProfile::fitDoublePeak ( Int_t iprofile,
				   Int_t iuv,
				   Float_t mean1,
				   Float_t mean2,
				   Float_t min,
				   Float_t max,
				   Option_t *opts
				   )
{
  // Fit the specified cluster's profile (either U or V plane) to
  // a double shower-shape peak, over the specfied range in strip
  // number, centered on the specified mean values.
  //
  // ... essentially, the first steps towards the "final" pi0 finder.  
 


  // Obtain pointer to the appropriate histogram
  assert ( iuv == 0 || iuv == 1 );
  TH1F *profile = ( iuv == 0 ) ? getProfileU(iprofile) : getProfileV(iprofile);

  // Set the range to that specified by the user
  if ( min < max )
    profile -> GetXaxis() -> SetRangeUser(min,max);


  // Perform a single-peak fit to the first point 
  fitSinglePeak(iprofile,iuv,mean1,min,max,opts);


  // Copy its parameters over to the double-peak fit,
  // also copy the limits
  for ( Int_t i = 0; i < 5; i++ ) {
    m_DoublePeak -> SetParameter(i, m_SinglePeak -> GetParameter(i) );
    Double_t pmin, pmax;
    m_SinglePeak -> GetParLimits(i, pmin, pmax );
    m_DoublePeak -> SetParLimits(i, pmin, pmax );
  }


  // Perform a single-peak fit to the second point
  fitSinglePeak(iprofile,iuv,mean2,min,max,opts);


  // And copy the new parameters and limits
  for ( Int_t i = 0; i < 5; i++ ) {
    m_DoublePeak -> SetParameter(i+5, m_SinglePeak -> GetParameter(i) );
    Double_t pmin, pmax;
    m_SinglePeak -> GetParLimits(i, pmin, pmax );
    m_DoublePeak -> SetParLimits(i+5, pmin, pmax );
  }

#if 0
  std::cout << "Double peak parameters" << std::endl;
  m_DoublePeak -> Print();
#endif

  // Now fit the profile to the double-peak shape
  profile -> Fit ( m_DoublePeak, opts );



}


 void EEsmdProfile::fitSinglePeak ( Int_t iprofile, 
				   Int_t iuv, 
				   Float_t mean,
				   Float_t min, 
				   Float_t max,
				   Option_t *opts ) 
{
  // Fit the specified cluster's SMD profile (either U or V plane) to
  // a single shower-shape peak, over the specified range in strip 
  // number, centered on the specified mean value.  If the histogram
  // is empty, no fit will be performed.  The fit is stored within
  // this histogram's list of fits (see TH1F::GetListOfFunctions()).

  // Obtain pointer to the appropriate histogram
  assert ( iuv == 0 || iuv == 1 );
  TH1F *profile = ( iuv == 0 ) ? getProfileU(iprofile) : getProfileV(iprofile);


#if 0
  std::cout << "fitSinglePeak(" 
	    << iprofile << "," 
	    << iuv << ","
	    << mean << ","
	    << min << ","
	    << max << ","
	    << opts << ");"
	    << std::endl << std::flush;
#endif

  assert(profile);



  // If the histogram is empty, we have nothing to do here...
  if ( profile -> GetEntries() == 0 ) return;

  // Set the range to that specified by the user
  if ( min < max )
    profile -> GetXaxis() -> SetRangeUser(min,max);


  // Initialize fit parameters, step sizes, and allowed ranges
  EEezCluster *cluster = m_EEezClust -> getClusterPtr(iprofile);
  estimateSinglePeak( cluster, iuv );

#if 0
  std::cout << "Initial parameterization" << std::endl << std::flush;
  m_SinglePeak -> Print();
#endif

  //-- Perform the fit based on initial parameters, and refit
  //-- over a reduced range, to limit impact on chi2 from large-
  //-- angle bremstrahlung and MCS.
  profile -> Fit( m_SinglePeak, opts );
  Float_t oldMean = m_SinglePeak -> GetParameter(1);
  m_SinglePeak -> SetRange( oldMean - 4.5, oldMean + 4.5 );    
  //-- Refit over new range
  profile -> Fit( m_SinglePeak, opts );
  Float_t myMean  = m_SinglePeak -> GetParameter(1);
  m_SinglePeak -> SetRange( myMean - 4.5, myMean + 4.5 );



  //-- Now refit until the means converge to 1/4 strip or less,
  //-- while freeing up our restriction on the width.

  m_SinglePeak -> SetParLimits(2, 0.1, 5.0);
  Int_t   count = 0;  
  while ( TMath::Abs( myMean - oldMean ) > 0.25 && count++ < 10 ) {   
    oldMean = m_SinglePeak -> GetParameter(1);
    profile -> Fit( m_SinglePeak, opts );
    myMean = m_SinglePeak -> GetParameter(1);
    m_SinglePeak -> SetRange( myMean - 4.5, myMean + 4.5 );
  }
  m_Converge -> Fill( count+0.5 );


  //-- Create a new histogram which will hold the chi2 for each 
  //-- strip. 
  TString hname = profile -> GetName();
  hname += "_Chi2_vs_strip";
  Int_t   nbin = profile -> GetNbinsX();
  Float_t myMin  = profile -> GetBinLowEdge(1);
  Float_t myMax  = profile -> GetBinLowEdge(nbin)+1.0;
  TH1F *histo = new TH1F(hname,"#chi^{2}_{i} = (Y_{i} - f(x_{i})^2 / #sigma_{i}^{2} vs. i",nbin,myMin,myMax);
  
  //-- Histogram bin containing the mean of the fit
  Int_t imean = profile -> FindBin( myMean );

  //-- Loop over the four bins to the left and right of the mean.
  //-- Abort the procedure if the fit doesn't have any degrees
  //-- of freedom, though.
  if ( m_SinglePeak -> GetNDF() > 0 )
    for ( Int_t i = imean - 4; i <= imean + 4; i++ ) {
      
      //-- Get the center of the bin
      Float_t x = profile -> GetBinCenter(i);
      
      //-- Evaluate the function at this point
      Float_t f = m_SinglePeak -> Eval(x);

      //-- Get the strip response and standard deviation for this strip
      Float_t y = profile -> GetBinContent(i);
      Float_t s = profile -> GetBinError(i);
      
      if ( s>0. ) {
	histo -> Fill ( x, (y-f)*(y-f)/s/s );
      }

    }

  //-- Now add the new histogram to the list of histograms
  if ( iuv == 0 ) 
    m_UProfileChi2.Add( histo );
  else
    m_VProfileChi2.Add( histo );

}

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

 TH1F *EEsmdProfile::getResidualU( Int_t i ) 
{
  // If it doesn't exist, it calculates the residual histogram for
  // the specified U profile, stores it in m_UResidual, and returns
  // a pointer to it.  Otherwise, returns a pointer to the histogram.
  // NOTE: After looking at root's TList::AddAt() function, I do not 
  // believe the documentation!  It looks like if you attempt to add
  // an object at a position beyond the end of the list, it blindly
  // assumes that you want to add the object at the end of the list...
  // This means that, for instance, if you were to call getResidualU(5)
  // BEFORE calling getResidual(1), the residual corresponding to the
  // 5th cluster will be inserted into the list at position 2... and
  // subsequent calls to getResidualU(2) will return the 5th cluster's
  // residual!  This is why I really hate CERN code.
  
  TH1F *resid = (TH1F *)m_UResidual.At(i);
  if ( !resid ) {
    TString name = getProfileU(i) -> GetName();
    name.ReplaceAll("hu","ru");
    resid = (TH1F*)getProfileU(i) -> Clone(name);
    TF1 *fit = (TF1*)getProfileU(i)->GetListOfFunctions()->At(0);
    if ( !fit ) return NULL;
    resid -> Add(fit,-1.0);
    m_UResidual.AddAt(resid,i);
  }
  assert(resid); // Something is really messed up here
  return resid;

}

 TH1F *EEsmdProfile::getResidualV( Int_t i )
{
  // If it doesn't exist, it calculates the residual histogram for
  // the specified V profile, stores it in m_UResidual, and returns
  // a pointer to it.  Otherwise, returns a pointer to the histogram.
  // See warnings in getResidualU().

  TH1F *resid = (TH1F *)m_VResidual.At(i);
  if ( !resid ) {
    TString name = getProfileV(i) -> GetName();
    name.ReplaceAll("hv","rv");
    resid = (TH1F*)getProfileV(i) -> Clone(name);
    TF1 *fit = (TF1*)getProfileV(i)->GetListOfFunctions()->At(0);
    if ( !fit ) return NULL;
    resid -> Add(fit,-1.0);
    m_VResidual.AddAt(resid,i);
  }
  assert(resid); // Something is really messed up here
  return resid;

}

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

 TF1 *EEsmdProfile::getSinglePeak()
{
  // Returns a pointer to the single-peak function used in the call
  // to fitSinglePeak().
  return m_SinglePeak;
}

 TF1 *EEsmdProfile::getDoublePeak()
{
  // Returns a pointer to the double-peak function used in the call
  // to fitDoublePeak().
  return m_DoublePeak;
}

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

 TH1F *EEsmdProfile::getSumOverStrips ( EEezCluster *cluster, Int_t nstrip, Int_t iuv ) 
{  
  //
  // Creates a new histogram whose entries correspond to the integral
  // over nstrip strips, versus the index of the strip in the center.
  // (NOTE-- please make sure that nstrips is odd.  It won't crash, but
  // it will be harder to interpret, and its not trapped.)  
  //
  // This histogram MUST be deleted by the user!  To retrieve this 
  // histogram once it has been created, getSumStrips ( EEezCluster *cluster ), 
  // or by index or name.
  //

  //-- Name of the cluster
  TString name_cluster = cluster -> Hash();

  TString histo = "h"; histo += (iuv)?"V":"U"; histo += nstrip; histo += "sum";
  histo += name_cluster;
  TString title = "#sum_{x}^{x+9} N_{mip}(x) vs x";

  //-- Get the range of strips for this histogram
  EEezTower *seed = cluster -> getSeedTower();
  EEezStripPtrVec_t *strips = (iuv)?seed->getVStrips():seed->getUStrips();

  EEezStripPtrVecIter_t begin = strips->begin();
  EEezStripPtrVecIter_t end   = strips->end() - 1;
  Int_t myMin = (*begin) -> getIndex() - 2;
  Int_t myMax = (*end)   -> getIndex() + 2;

  // Create the new histogram
  TH1F *sum9 = new TH1F(histo,title,(myMax-myMin),myMin,myMax);
  
  // And fill it
  while ( begin < strips -> end() ) {
    Float_t sum   = 0;
    Float_t index = (*begin) -> getIndex() + nstrip/2;
    EEezStripPtrVecIter_t iter = begin;
    while ( iter < begin+nstrip && iter < strips->end() ) {
      sum += (*iter) -> getNumMips();
      iter++;
    }
    sum9 -> Fill (index, sum );
    begin++;
  }

  //$$$ Info("getSumOverStrips","returning profile sum %s", sum9 -> GetName() );

  return sum9;

}


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

 void EEsmdProfile::estimateSinglePeak ( EEezCluster *cluster, Int_t iuv )
{
  // For the specified cluster and plane, initialize starting parameters
  // for the single-peak fit

  //-- Save the first 50 events as examples...
  static Int_t nsave = 0;



  //-- Open an integration window, 5 and 9 strips wide.  Slide this
  //-- window from first strip to last strip, and fill a histogram
  //-- with the integrated number of mips inserted at the starting
  //-- point of the integration + 1/2 the window.  This distribution
  //-- should reach maximum at the position of the maximum bin.
  TH1F *sum5 = getSumOverStrips( cluster, 5, iuv );
  TH1F *sum9 = getSumOverStrips( cluster, 9, iuv );
  if ( iuv == 0 ) {
    m_UProfileSum9.Add(sum9);
    m_UProfileSum5.Add(sum5);
  }
  else {
    m_VProfileSum9.Add(sum9);
    m_VProfileSum5.Add(sum5);
  }

  Int_t max5 = sum5 -> GetMaximumBin();
  //$$$  Int_t max9 = sum9 -> GetMaximumBin();

  Float_t mean5 = sum5 -> GetBinLowEdge( max5 ) + 0.5;
  //$$$ Float_t mean9 = sum9 -> GetBinLowEdge( max9 ) + 0.5;

  //  Float_t sumAtMax5 = sum5 -> GetBinContent( max5 );
  //  Float_t sumAtMax9 = sum9 -> GetBinContent( max9 );
  
  //-- We estimate the mean as the average of the two preceding values.
  //-- This will reduce the effects of outlying strips.  In the future,
  //-- we may consider investigating a weighted average of some sort.
  //$$$  Float_t start_mean = 0.5 * ( mean5 + mean9 );
  Float_t start_mean = mean5;  
  Float_t min_mean = start_mean - 4.5;
  Float_t max_mean = start_mean + 4.5;
  Int_t   bin_mean = sum5 -> FindBin( start_mean );



  //-- We estimate the width of the distribution using the integration-
  //-- window histograms above.  When the integral is reduced by
  //-- ~32% from the maximum, we know that we are off by 1 sigma.
  //-- We will apply a sanity check to the result... if its less than
  //-- 0.5 strips or greater than 3.0 strips, we start our estimate 
  //-- at 1.0 strips wide.

  //---> Note: for now, just start with a constant
  Float_t start_sigma = 0.84;
  Float_t min_sigma   = 0.65 * start_sigma;
  Float_t max_sigma   = 2.00 * start_sigma;

  //-- Next, we estimate the integrated yield...  We use the integral
  //-- over 9 strips, evaluated at our starting mean.
  Float_t start_yield = sum9 -> GetBinContent( bin_mean );
  Float_t min_yield   = 0.1 * start_yield;
  Float_t max_yield   = 2.5 * start_yield;


  //-- m_SinglePeak is a 5 parameter fit: 2 gaussians w/ common
  //-- means + a constrained width.  Initially, we fix the relative
  //-- fraction of the two gaussians and the relative widths.
  m_SinglePeak -> SetParameter( 0, start_yield );
  m_SinglePeak -> SetParameter( 1, start_mean );
  m_SinglePeak -> SetParameter( 2, start_sigma );
  m_SinglePeak -> FixParameter( 3, 0.2 );
  //$$$m_SinglePeak -> FixParameter( 4, 3.5 );
  m_SinglePeak -> SetParameter( 3, 0.2 );
  m_SinglePeak -> SetParameter( 4, 3.5 );
 
  m_SinglePeak -> SetParLimits( 0, min_yield, max_yield );
  m_SinglePeak -> SetParLimits( 1, min_mean,  max_mean );
  m_SinglePeak -> SetParLimits( 2, min_sigma, max_sigma );

  //  m_SinglePeak -> SetParLimits( 3, 0.05, 0.45 );
  //  m_SinglePeak -> SetParLimits( 4, 1.1,  4.0 );



  //-- Clone and save initial parameterization for fits


  TH1F *profile = (iuv)?getProfileV(cluster):getProfileU(cluster);

  TString name = profile -> GetName();

  name += "_init_"; name += nsave;

  TH1F *cloned = (TH1F *)profile -> Clone(name);
  assert(cloned);
  cloned -> GetListOfFunctions() -> Add( m_SinglePeak -> Clone("init") ); 

  if ( iuv == 0 ) 
    m_UProfileInit.Add(cloned);
  else
    m_VProfileInit.Add(cloned);


}



#if 0
  //-- We save the first 50 events as examples
  if ( nsave < 50 ) {

    cd();

    //-- Get the specified profile histogram
    //$$$H1F *profile = (iuv)?getProfileV(cluster):getProfileU(cluster);
    //  TString name = profile -> GetName();
    //  name += "_par_init_"; name += nsave;
    //  TH1F *cloned = (TH1F *)profile -> Clone(name);
    //  cloned -> GetListOfFunctions() -> Add( m_SinglePeak -> Clone("init") );

    //-- Save the sum5 and sum9 histograms
    name = profile -> GetName();
    name += "_sum5_"; name += nsave;
    TH1F *clone5 = (TH1F *)sum5 -> Clone(name);

    name = profile -> GetName();
    name += "_sum9_"; name += nsave;
    TH1F *clone9 = (TH1F *)sum9 -> Clone(name);

    nsave++;

    cd("..");


  } 
#endif


#if 0


  m_DiffPar[0] = new TH1F("diffPar0","Difference between inital yield estimate and fit",40,-20.,20.);
  m_DiffPar[1] = new TH1F("diffPar1","Difference between inital mean estimate and fit",20,-5.,5.);
  m_DiffPar[2] = new TH1F("diffPar2","Difference between inital width estimate and fit",100,-2.,2.);
  m_DiffPar[3] = new TH1F("diffPar3","Difference between inital frac. yield estimate and fit",20,-1.,1.);
  m_DiffPar[4] = new TH1F("diffPar4","Difference between inital rel. width estimate and fit",20,-5.,5.);

  m_Par[0] = new TH1F("Par0","fit yields",40,0.,400.); 
  m_Par[1] = new TH1F("Par1","fit means",288,0.,288.);
  m_Par[2] = new TH1F("Par2","fit widths",20,0.,5.);
  m_Par[3] = new TH1F("Par3","fit frac. yields",20,0.,1.);
  m_Par[4] = new TH1F("Par4","fit rel. widths",20,0.,5.);

  m_DiffParVsChi2[0] = new TProfile("diffPar0VsChi2","Difference between inital yield estimate and fit vs chi2/ndf",40,0.,10.,-1000.,1000.);
  m_DiffParVsChi2[1] = new TProfile("diffPar1VsChi2","Difference between inital mean estimate and fit vs chi2/ndf",40,0.,10.,-1000.,1000.);
  m_DiffParVsChi2[2] = new TProfile("diffPar2VsChi2","Difference between inital width estimate and fit vs chi2/ndf",40,0.,10.,-1000.,1000.);
  m_DiffParVsChi2[3] = new TProfile("diffPar3VsChi2","Difference between inital frac. yield estimate and fit vs chi2/ndf",40,0.,10.,-1000.,1000.);
  m_DiffParVsChi2[4] = new TProfile("diffPar4VsChi2","Difference between inital rel. width estimate and fit vs chi2/ndf",40,0.,10.,-1000.,1000.);

  m_Chi2VsRMS = new TH2F("chi2VsRMS","#chi^2/ndf vs RMS of the histogram",50,0.,10.,50,0.,20.);

  m_NDF = new TH1F("ndf","Degrees of freedom in single shower shape fit",10,0.,10.);
  

#endif

#if 0

  m_DiffPar[0] -> Fill( nmip - m_SinglePeak -> GetParameter(0) );
  m_DiffPar[1] -> Fill( mean - m_SinglePeak -> GetParameter(1) );
  m_DiffPar[2] -> Fill( width - m_SinglePeak -> GetParameter(2) );
  m_DiffPar[3] -> Fill( frac - m_SinglePeak -> GetParameter(3) );
  //$$$  m_DiffPar[4] -> Fill( nmip - m_SinglePeak -> GetParameter(4) );

  Float_t chi2 = m_SinglePeak -> GetChisquare();
  m_Chi2VsRMS -> Fill( rms, chi2 );

  m_DiffParVsChi2[0] -> Fill( chi2, nmip - m_SinglePeak -> GetParameter(0) );
  m_DiffParVsChi2[1] -> Fill( chi2, mean - m_SinglePeak -> GetParameter(1) );
  m_DiffParVsChi2[2] -> Fill( chi2, width - m_SinglePeak -> GetParameter(2) );
  m_DiffParVsChi2[3] -> Fill( chi2, frac - m_SinglePeak -> GetParameter(3) );
  //$$  m_DiffParVsChi2[4] -> Fill( chi2, nmip - m_SinglePeak -> GetParameter(4) );
  
  for ( Int_t i = 0; i < 5; i++ ) m_Par[i] -> Fill( m_SinglePeak -> GetParameter(i) );

  m_NDF -> Fill ( m_SinglePeak -> GetNDF() );

  // Calculate the residual and save in the appropriate list
  TH1F *resid = (TH1F*)profile -> Clone();
  TString name = profile -> GetName();
  if ( iuv == 0 ) 
    name.ReplaceAll("hu","ru");
  else
    name.ReplaceAll("hv","rv");
  resid -> SetName(name);
  resid -> Add(m_SinglePeak,-1.0);

  if ( iuv == 0 ) 
    m_UResidual.Add( resid );
  else 
    m_VResidual.Add( resid );
  
#endif