MaxBinCFAPGenerator.cpp

Go to the documentation of this file.

#include "MaxBinCFAPGenerator.h"
#include "TPulseIsland.h"
#include "TAnalysedPulse.h"
#include "TVAnalysedPulseGenerator.h"

#include <algorithm>
#include <vector>

// IsTimeOrdered()
// -- Returns tru of the first pulse is before the second
// Static so that only this file sees it
// Useful for sorting the pulses
static bool IsTimeOrdered(TAnalysedPulse* a, TAnalysedPulse* b) {
  // Want a to be before b
  return ( a->GetTime() < b->GetTime() );
}

MaxBinCFAPGenerator::MaxBinCFAPGenerator(double fraction) : fFraction(fraction) {}

// Constant Fraction
// The time the pulse reaches some fraction of its height
// is recorded as the time of the pulse.
// Max Bin
// The height of the largest bin (pedestal subtracted and
// flipped according to its polarity) is considered the
// amplitude.
void MaxBinCFAPGenerator::ProcessPulses(const TSetupData* eventSetup,
                                      const PulseIslandList_t& pulseList,
                                      AnalysedPulseList_t& analysedList){

  std::string bank = pulseList[0]->GetBankName();
  int ped = gSetup->GetPedestal(bank);
  int pol = gSetup->GetTriggerPolarity(bank);

  for (PulseIslandList_t::const_iterator iPulse = pulseList.begin(); iPulse != pulseList.end(); ++iPulse) {
    std::vector<int> samples = (*iPulse)->GetSamples();
    int t1 = 0; // Time tick right before crossing CF level
    int t2 = 1; // Time tick right after crossing CF level
    int tmax = 1; // Time tick of maximum value
    int s1 = pol * (samples[t1] - ped); // Sample value right before crossing CF level
    int s2 = pol * (samples[t2] - ped); // Sample value right after crossing CF level
    int smax = s2; // Sample value at maximum
    double cf = fFraction * (double)smax;
    for (int iSample = 2; iSample < (int)samples.size(); ++iSample) {
      if (pol * (samples[iSample] - ped) > smax) {
        tmax = iSample;
        smax = pol * (samples[tmax] - ped);
        cf = fFraction * (double)smax;
        while (s2 < (int)cf && s2 <= smax) {
          ++t1; ++t2;
          s1 = pol * (samples[t1] - ped); s2 = pol * (samples[t2] - ped);
        }
      }
    }

    // Some sanity checks
    if (s1 > (int)cf || s2 < (int)cf || s1 >= s2 || smax < 0)
      std::cout << "MaxBinCFAPGenerator: Warning: Sanity check failed.." std::endl;
    
    // Linear interpolation of constant fraction time
    // Then shift and scale by tick size
    double t0 = gSetup->GetTimeShift(bank);
    double dt = gSetup->GetClockTick(bank);
    double t = (double)(t2 - t1) / (double)(s2 - s1) * (cf - (double)s1) + (double)t1; // Units of clock tick
    t = t*dt + t0;

    // Add the pulse into the list
    analysedList.push_back(new  TAnalysedPulse((double)smax, t, 0., 0., gSetup->GetDetectorName(bank)));

    std::sort(analysedList.begin(), analysedList.end(), IsTimeOrdered);

}