LIPatternFinderSimple.cxx

Go to the documentation of this file.

#include <cmath>
#include <numeric>
#include <cassert>
#include <algorithm>
#include <functional>

#include "MessageService/MsgService.h"

#include "LIPatternFinderSimple.h"

ClassImp(LIPatternFinderSimple)

//_____________________________________________________________________
CVSID(
"$Id: LIPatternFinderSimple.cxx,v 1.5 2006/05/27 08:58:13 rhatcher Exp $"
);
//_____________________________________________________________________

// STL Functors

struct is_east_side: public unary_function<CandDigitHandle *, bool>
{
   bool operator()(CandDigitHandle * digit) {
     
         const PlexSEIdAltL & seid = digit->GetPlexSEIdAltL();
         return (seid.GetEnd() == StripEnd::kEast);
   }
};
//.....................................................................
struct sum_q :
              public binary_function<double, CandDigitHandle *, double>
{
   double operator()(double sum, CandDigitHandle * current_digit) {
     
          return (sum + current_digit->GetCharge());
   }
};
//.....................................................................
struct is_in_plane_and_above_thr:
                       public unary_function<CandDigitHandle *, bool>
{
   int    _plane;
   double _threshold;

   is_in_plane_and_above_thr(int plane, double threshold) :
                              _plane(plane), _threshold(threshold) { }
                              
   bool operator()(CandDigitHandle * digit) {

      const PlexSEIdAltL & seid = digit->GetPlexSEIdAltL();
      
      return ( (digit->GetCharge() > _threshold) &&
                                        (seid.GetPlane() == _plane) );
   }
};
//.....................................................................
struct is_in_plane_range_and_above_thr:
                        public unary_function<CandDigitHandle *, bool>
{
   int    _pl_min, _pl_max;
   double _threshold;

   is_in_plane_range_and_above_thr(
        int pl_min, int pl_max, double threshold) :
           _pl_min(pl_min), _pl_max(pl_max), _threshold(threshold) { }

   bool operator()(CandDigitHandle * digit) {

      const PlexSEIdAltL & seid = digit->GetPlexSEIdAltL();

      return ( (digit->GetCharge() >  _threshold) &&
                             (seid.GetPlane() >= _pl_min) &&
                                        (seid.GetPlane() <= _pl_max) );
   }
};
//.....................................................................
struct min_plane:
   public binary_function<CandDigitHandle *, CandDigitHandle *, bool> {
   bool operator()(CandDigitHandle * a, CandDigitHandle * b) {
          return a->GetPlexSEIdAltL().GetPlane() <
                                       b->GetPlexSEIdAltL().GetPlane();
   }
};
//_____________________________________________________________________
LIPatternFinderSimple::LIPatternFinderSimple()
{
  MSG("LIPatternFinder", Msg::kVerbose)
                                     << "LIPatternFinderSimple ctor\n";
  DefaultConfig();
}
//_____________________________________________________________________
LIPatternFinderSimple::~LIPatternFinderSimple()
{
  MSG("LIPatternFinder", Msg::kVerbose)
                                     << "LIPatternFinderSimple dtor\n";
}
//_____________________________________________________________________
void LIPatternFinderSimple::DefaultConfig(void)
{
  // see header file for comments on variables
  
  fMaxEWAsymmetry               = 5.e+5; 
  fPulseHeightThreshold         = 200; 
  fMaxNDigits                   = 1000;
  fMaxCharge                    = 1.e+6;
  fHighActivityThreshold        = 0.90;
  fLowActivityThreshold         = 0.20;
}
//_____________________________________________________________________
void LIPatternFinderSimple::Configure(const Registry & registry)
{
  double tmpd;
  
  if (registry.Get("MaxEWAsymmetry", tmpd))    fMaxEWAsymmetry = tmpd;
  if (registry.Get("MaxNDigits", tmpd))            fMaxNDigits = tmpd;
  if (registry.Get("MaxCharge", tmpd))              fMaxCharge = tmpd;
  if (registry.Get("PulseHeightThreshold", tmpd))
                                         fPulseHeightThreshold = tmpd;
  if (registry.Get("HighActivityThreshold", tmpd))
                                        fHighActivityThreshold = tmpd;
  if (registry.Get("LowActivityThreshold", tmpd))
                                         fLowActivityThreshold = tmpd;
}
//_____________________________________________________________________
bool LIPatternFinderSimple::IsLightInjectionTrash(
                                            CandDigitListHandle * cdlh)
{
  MSG("LIPatternFinder", Msg::kVerbose)
                            << "LIPatternFinderSimple running" << endl;

  CandDigitHandleItr cdh_iter ( cdlh->GetDaughterIterator() );

  vector<CandDigitHandle *> digits;
  
  while ( CandDigitHandle * digit = cdh_iter.Ptr() ) {

    digits.push_back(digit);
        
    cdh_iter.Next();
  }

  sort(digits.begin(), digits.end(), min_plane()); // sort in pl.-nu.
  
  //PrintSnarl(digits);

  double asymmetry  = GetEWAsymmetry(digits);
  double tot_charge = GetTotalCharge(digits);

  vector<double> segment_activity =
                                ActivityBetweenCrateBoundaries(digits);

  sort(segment_activity.begin(), segment_activity.end());

  double first_max  = segment_activity[ segment_activity.size()-1 ];
  double second_max = segment_activity[ segment_activity.size()-2 ];
  
  MSG("LIPatternFinder", Msg::kVerbose)
          << "------------------------- LIPatternFinderSimple summary:"
          << endl;             
  MSG("LIPatternFinder", Msg::kVerbose)
                           << "E-W Asymmetry  = " << asymmetry << endl;
  MSG("LIPatternFinder", Msg::kVerbose)
                            << "N-Digits  = " << digits.size() << endl;
  MSG("LIPatternFinder", Msg::kVerbose)
                           << "Total charge  = " << tot_charge << endl;
  MSG("LIPatternFinder", Msg::kVerbose)
       << "Fraction of active planes between crate boundaries" << endl;
  MSG("LIPatternFinder", Msg::kVerbose)
                            << "    - 1st max = " << first_max << endl;
  MSG("LIPatternFinder", Msg::kVerbose)
                           << "    - 2nd max = " << second_max << endl;

  //-- check for large EAST-WEST asymmetry
  bool exceeds_max_EW_asymmetry = (asymmetry     > fMaxEWAsymmetry);

  //-- check if event exceeds a maximum number of digits
  bool exceeds_max_n_digits     = (digits.size() > fMaxNDigits);

  //-- check if event exceeds a maximum charge
  bool exceeds_max_charge       = (tot_charge    > fMaxCharge);

  //-- check if event activity is confined to a single segment (=
  //   area between two crate boundaries). To do this I have
  //   calculated the fraction of active planes in each segment
  //   and check whether the 1st maximum is very high while the
  //   2nd maximum is low.
  bool shows_segment_asymmetry  =
               (first_max  > fHighActivityThreshold) && 
                                 (second_max < fLowActivityThreshold);

  bool is_LI = exceeds_max_EW_asymmetry &&
                          exceeds_max_n_digits &&
                                    exceeds_max_charge &&
                                              shows_segment_asymmetry;  
  return is_LI;
}
//_____________________________________________________________________
double LIPatternFinderSimple::GetEWAsymmetry(
                        const vector<CandDigitHandle *> & digits) const
{
// returns the net |WEST side - EAST side| charge
  
  //-- copy input vector

  vector<CandDigitHandle *> digits_cp( digits.size() );
  
  copy( digits.begin(), digits.end(), digits_cp.begin() );

  //-- East-West partition
  
  vector<CandDigitHandle *>::iterator part_iter = partition(
                 digits_cp.begin(), digits_cp.end(), is_east_side() );

  //-- accumulate charge
  
  double q_east = accumulate(
                          digits_cp.begin(), part_iter, 0.0, sum_q());
  double q_west = accumulate(
                          part_iter, digits_cp.end(), 0.0, sum_q());

  MSG("LIPatternFinder", Msg::kVerbose) << "Q_W = " << q_west << endl;
  MSG("LIPatternFinder", Msg::kVerbose) << "Q_E = " << q_east << endl;
                          
  //-- calculate asymmetry

  double asymmetry = fabs(q_west-q_east);

  return asymmetry; 
}
//_____________________________________________________________________
double LIPatternFinderSimple::GetTotalCharge(
                        const vector<CandDigitHandle *> & digits) const
{
// returns the total charge

  double q = accumulate(digits.begin(), digits.end(), 0.0, sum_q());

  return q;
}
//_____________________________________________________________________
double LIPatternFinderSimple::FractionOfActivePlanesInRange(
               int plane_min, int plane_max, 
                        const vector<CandDigitHandle *> & digits) const
{
// returns the fraction of active planes in the range [pl_min, pl_max]

  int n_total_planes  = 1 + plane_max - plane_min;
  int n_active_planes = 0;
  
  assert(n_total_planes > 0);

  for(int plane = plane_min; plane < plane_max; plane++) {

    int n_digits = count_if(digits.begin(), digits.end(),
              is_in_plane_and_above_thr(plane, fPulseHeightThreshold));

    if(n_digits > 0) n_active_planes++;  
  }

  double fraction = (double) n_active_planes / (double) n_total_planes;

  return fraction;
}
//_____________________________________________________________________
vector<double>
     LIPatternFinderSimple::ActivityBetweenCrateBoundaries(
                        const vector<CandDigitHandle *> & digits) const
{
  vector<double> active_fraction;

  const int kNDetectorSegments = 8;
  
  const int kMinPlane[kNDetectorSegments] = {
                               1,  65, 129, 193, 250, 314, 378, 442 };
  const int kMaxPlane[kNDetectorSegments] = {
                              64, 128, 192, 248, 313, 377, 441, 485 };

  for(int isegment = 0; isegment < kNDetectorSegments; isegment++)
         active_fraction.push_back(
              FractionOfActivePlanesInRange(
                   kMinPlane[isegment], kMaxPlane[isegment], digits) );
                
  return active_fraction;
}
//_____________________________________________________________________
void LIPatternFinderSimple::PrintSnarl(
                        const vector<CandDigitHandle *> & digits) const
{
  vector<CandDigitHandle *>::const_iterator digit_iter;

  for(digit_iter = digits.begin();
                            digit_iter != digits.end(); ++digit_iter) {

    const PlexSEIdAltL & seid = (*digit_iter)->GetPlexSEIdAltL();
    RawChannelId rawchid = (*digit_iter)->GetChannelId();
    
    MSG("LIPatternFinder", Msg::kVerbose)
         << " PL = "    << seid.GetPlane()
         << " END = "   << seid.GetEnd()
         << " CRATE = " << rawchid.GetCrate()
         << " VARC = "  << rawchid.GetVarcId()
         << " VMM = "   << rawchid.GetVmm()    
         << " Q = "     << (*digit_iter)->GetCharge() << endl;
  }
}
//_____________________________________________________________________

---