/////////////////////////////////////////////////////////////////////////////
//
// EEezPi0nCandidate
//
// Base class for pi0 --> 2 gamma candidates.  Candidates can be assembled
// by adding two clusters or by adding individual towers and SMD strips.
// As individual towers and strips are added, they are checked against
// the pi0 candiate's seed thresholds.  Towers which exceed the tower
// threshold are added to the list of seed towers in the candidate.  
//
// Strips which exceed a seed threshold are added to a list of potential
// seed strips.  The numberSmdSeeds() method can be used to count the
// number of SMD seeds, which satisfy a condition on the number and
// sum of adjacent strips.  In other words, we generally consider 3 adjacent 
// strips summing to > 20.0 MIPs to be an SMD seed.
//
/////////////////////////////////////////////////////////////////////////////

#include "EEezPi0nCandidate.h"
#include <iostream>

ClassImp(EEezPi0nCandidate);

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

 EEezPi0nCandidate::EEezPi0nCandidate()
{
  // Class Constructor
  mTowerThreshold = 0.7;
  mStripThreshold = 1.0;
}

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

 void EEezPi0nCandidate::setTowerThreshold( Float_t t )
{
  // Sets the seed threshold for towers.
  mTowerThreshold = t;
}

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

 void EEezPi0nCandidate::setStripThreshold( Float_t s )
{
  // Sets the seed threshold for strips.
  mStripThreshold = s;
}

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

 void EEezPi0nCandidate::addTower( EEezTower *tower )
{
  // Adds the tower to the list of towers.  If it exceeds the
  // seed threshold, it is also added to the list of seed towers.
  // Rejects duplicate entries.
  
  if ( hasTower(tower) ) return;

  //-- Add the tower to seeds if above threshold
  if ( tower -> getEnergy() > mTowerThreshold ) {
    mSeedTowers.push_back( tower );
  }
  
  //-- Add the tower to the tower list
  EEezPatch::addTower( tower );

}

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

 void EEezPi0nCandidate::addStrip( EEezStrip *strip )
{
  // Adds the strip to the list of strips.  If it exceeds
  // the seed strip threshold, it is also added to the list 
  // of potential seed strips.  Rejects duplicate entries.

  if ( hasStrip(strip) )return;

  if ( strip -> getNumMips() > mStripThreshold ) {
    if ( strip -> getPlane() == 0 ) 
      mSeedUStrips.push_back( strip );
    else
      mSeedVStrips.push_back( strip );
  }

  EEezPatch::addStrip( strip );

}

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

 void EEezPi0nCandidate::addCluster( EEezCluster *cluster )
{
  // Add the towers and underlying SMD strips in the cluster
  // to the patch.  Identify seed towers and potential seed
  // strips as well.
  // Need to increase speed.  Too many redundant
  // calls to hasTower() and hasStrip() are possible here.

  //-- Get the seed tower of the cluster and add it to
  //-- the patch.
  EEezTower *seed = cluster -> getSeedTower();
  EEezPatch::addTower(seed);
  if ( seed -> getEnergy() > mTowerThreshold ) addTower(seed);
  Int_t nu = seed -> getNUStrips();
  Int_t nv = seed -> getNVStrips();
  for ( Int_t j = 0; j < nu; j++ ) {
    EEezStrip *s = seed->getUStrip(j);
    if ( s -> getNumMips() > mStripThreshold ) addStrip(s);
  }
  for ( Int_t j = 0; j < nv; j++ ) {
    EEezStrip *s = seed->getVStrip(j);
    if ( s -> getNumMips() > mStripThreshold ) addStrip(s);
  }
 
  //-- Loop over all neighbors, and add them to the patch
  for ( Int_t i = 0; i < seed -> getNNeighbors(); i++ ) {

    EEezTower *tower = seed -> getNeighbor(i);
    EEezPatch::addTower( tower );

    nu = tower -> getNUStrips();
    nv = tower -> getNVStrips();
    for ( Int_t j = 0; j < nu; j++ ) {
      EEezStrip *s = tower->getUStrip(j);
      EEezPatch::addStrip(s);
      if ( s -> getNumMips() > mStripThreshold ) addStrip(s);
    }
    for ( Int_t j = 0; j < nv; j++ ) {
      EEezStrip *s = tower->getVStrip(j);
      EEezPatch::addStrip(s);
      if ( s -> getNumMips() > mStripThreshold ) addStrip(s);
    }
  }

}

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

 EEezPi0nCandidate::~EEezPi0nCandidate()
{
  // The Destructor
}

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

 Int_t EEezPi0nCandidate::numberSmdSeeds ( Int_t iuv, Float_t threshold )
{
  // Count SMD seeds.  To be defined...

  //-- Clone the profile histogram for this candidate
  TH1F *profile = (TH1F *)((iuv==0)?profileU()->Clone("profile"):profileV()->Clone("profile"));

  //-- Some histogram info... remember, root is braindead
  //-- and indexes bins from 1.
  Int_t nbins  = profile -> GetNbinsX();

  Int_t nseeds = 0;
  Int_t ntries = 0; // Will abort after 10 tries with no seeds
  while (1) {

    if ( ntries++ >= 10 ) break;

    //-- Get the (current) maximum bin, and trap it if it looks
    //-- like a questionable pedestal or MIP
    Int_t maxbin = profile -> GetMaximumBin();
    if ( profile -> GetBinContent(maxbin+1) *
	 profile -> GetBinContent(maxbin-1) == 0 ) {
      profile -> SetBinContent(maxbin,0.);
      continue;
    }

    //-- Restrict the range of the histogram to +/- 5 strips
    //-- around the maximum bin.
    Int_t min = TMath::Max(1,maxbin-5);
    Int_t max = TMath::Min(nbins,maxbin+5);
    if ( profile -> Integral(min+1,max-1) > threshold ) {
      nseeds++;
      ntries--;
    }
    //-- Wipe these strips off for the next search
    for ( Int_t i = min; i <= max; i++ ) 
      profile -> SetBinContent(i,0.);
 
    if ( profile -> Integral() < threshold ) break;

  }
  
  //-- Remove the histogram, we're done with it.
  delete profile;

  return nseeds;

}


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

#if 0
 Int_t EEezPi0nCandidate::numberSmdSeeds ( Int_t iuv, Float_t threshold )
{
  // Count the number of SMD seeds-- three adjacent strips which sum above
  // the specified threshold (default 25 mips).  iuv: U=0, V=1.
 

  //-- Get a pointer to the right plane of strips
  EEezStripPtrVec_t *strips = (iuv==0)?&mUStrips:&mVStrips;

  Int_t nstrips = strips -> size();
  if ( nstrips < 3 ) return 0;

  Int_t nseed = 0;
  for ( Int_t i = 1; i < nstrips-1; i++ ) {

    //-- Adjacency
    EEezStrip *s1 = strips -> at(i-1);
    EEezStrip *s2 = strips -> at(i);
    EEezStrip *s3 = strips -> at(i+1);

    if ( (s1->getIndex() != s2->getIndex() - 1) ||
	 (s3->getIndex() != s2->getIndex() + 1) ) continue;

    Float_t sum = s1->getNumMips()+s2->getNumMips()+s3->getNumMips();
    if ( sum < threshold ) continue;

    i = i + 2;
    nseed++;

  }

  return nseed;

}
#endif

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