AlgChopListSharp.cxx

Go to the documentation of this file.

// $Id: AlgChopListSharp.cxx,v 1.5 2010/01/06 17:15:44 rhatcher Exp $
//
// AlgChopListSharp.cxx
//

#include <cassert>
#include <cmath>

#include "CandChop/AlgChopListSharp.h"
#include "CandChop/CandChopListHandle.h"
#include "CandChop/DigitVector.h"


#include "Algorithm/AlgConfig.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "CandData/CandHeader.h"
#include "CandData/CandRecord.h"
#include "CandDigit/CandDigitHandle.h"
#include "CandDigit/CandDigitListHandle.h"
#include "CandDigit/CandDigitList.h"
#include "Candidate/CandContext.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "RawData/RawHeader.h"
#include "RawData/RawRecord.h"
#include "RawData/RawDigitDataBlock.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliStripHandle.h"
#include "Validity/VldContext.h"
#include "Calibrator/Calibrator.h"
#include "DataUtil/GetRawBlock.h"

ClassImp(AlgChopListSharp)
CVSID( " $Id: AlgChopListSharp.cxx,v 1.5 2010/01/06 17:15:44 rhatcher Exp $ ");

struct compareDigitTimes : public binary_function<const CandDigitHandle&, const CandDigitHandle&, bool> {
  bool operator()(const CandDigitHandle& d1, const CandDigitHandle& d2) {
    return (d1.GetTime() < d2.GetTime());
  }
};

const RawChannelId kQieRcid(Detector::kNear,ElecType::kQIE,0,0,false,false);
const float k1pe = 100.; // About 1 pe.


//______________________________________________________________________
AlgChopListSharp::AlgChopListSharp()
{
}

//______________________________________________________________________
AlgChopListSharp::~AlgChopListSharp()
{
}

//______________________________________________________________________
bool AlgChopListSharp::ShouldSplit( float this_ph, // sigcorrs in current bin
                                      float next_ph, // sigcorrs in next/prev bin we want to expand into
                                    float /*d_tmax*/  // distance from next bin to peak of event.
                                      )
{
  float climb = next_ph - this_ph;
  
  // the maximum delta that the algorithm will climb before making a new chop
  float max_climb = 2.5*sqrt(fabs(this_ph)/k1pe)*k1pe;
  if(max_climb < (6 * k1pe)) max_climb=max_climb*2;

  // the maximum pulse height in this bin if we're making a new chop.
  //const float size_limit = 20 * k1pe; 
  //const float min_time = 2;

  //if(d_tmax  < min_time) return false;
  //if(this_ph < size_limit) return false;

  if( climb >= max_climb ) return true;
  else return false;    
}


//______________________________________________________________________
void AlgChopListSharp::RunAlg(AlgConfig& /*algConfig*/, 
                           CandHandle &candHandle,  // thing to make
                           CandContext &candContext)
{

  assert(candHandle.InheritsFrom("CandChopListHandle"));
  CandChopListHandle &chopList = dynamic_cast<CandChopListHandle &>(candHandle);

   assert(candContext.GetDataIn());
   // Check for CandDigitListHandle input
   if (!(candContext.GetDataIn()->InheritsFrom("CandDigitListHandle"))) {
     MSG("Chop",Msg::kWarning ) << "Data into AlgChopListSharp is not a digit list." << std::endl;
   }
   
   const CandDigitListHandle *cdlh_ptr = 
     dynamic_cast<const CandDigitListHandle*>(candContext.GetDataIn());
   
   const MomNavigator *mom = candContext.GetMom();
   const RawDigitDataBlock* rddb = DataUtil::GetRawBlock<RawDigitDataBlock>(mom);
   if (!rddb) {
     MSG("Chop", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
     return;
   }
   
   // Get setup for the DigitList maker algorithm
   AlgFactory &af = AlgFactory::GetInstance();
   AlgHandle ah = af.GetAlgHandle("AlgChop","default");
   CandContext cxx(this,candContext.GetMom());

   const VldContext &context = *(candContext.GetCandRecord()->GetVldContext());
   if(context.GetDetector() != Detector::kNear) 
     MSG("Chop",Msg::kError) << "Running the Sharp algorithm on FD data is a no-no!" << endl;

   Calibrator& cal = Calibrator::Instance();
   UgliGeomHandle ugli(context);
  
   // Now do the actual slicing.

   // First, make a nice stl vector of the digits.
   DigitVector digits(cdlh_ptr);

   UInt_t ndigits = digits.size();
   UInt_t nchop = 0;
    
   // Sort the list by time.
   // Not neccessary for this algorithm.
   // std::sort(digits.begin(), digits.end(), compareDigitTimes());


   // Also, I want some other pieces of info:
   std::vector<int>    digit_tdc(ndigits);
   std::vector<UInt_t> digit_plane(ndigits);
   //std::vector<float>  digit_tpos(ndigits);
   for(UInt_t i=0;i<ndigits;i++) {
     digit_tdc[i] = (cal.GetTDCFromTime(digits[i].GetTime(CalTimeType::kNone), kQieRcid));
     digit_plane[i] = digits[i].GetPlexSEIdAltL().GetPlane(); 
     //if(digit_plane[i]<=PlexPlaneId::LastPlaneNearCalor())
     //  digit_tpos[i]  = ugli.GetStripHandle(digits[i].GetPlexSEIdAltL().GetBestSEId()).GetTPos(); 
     //else 
     //  digit_tpos[i]  = -999;
   }

   // Find first and last times. Add some padding so sertain operations are easier to code.
   Int_t tfirst = digit_tdc[0];
   Int_t tlast  = digit_tdc[0];
   for(UInt_t i=0;i<ndigits;i++) {
     if(digit_tdc[i] < tfirst) tfirst = digit_tdc[i];
     if(digit_tdc[i] > tlast ) tlast  = digit_tdc[i];
   }
   tfirst-=5;
   tlast +=5;


   // Make the energy histogram.
   MSG("Chop",Msg::kDebug) << "Running Chop_Sharp" << endl;
      
   UInt_t numBins = tlast-tfirst;

   // Create the energy-time profile.
   std::vector<float> energyVsTime(numBins,0.);
   
   for(UInt_t idig = 0; idig < ndigits; idig++ ) {
     float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr);
     int tdcbin = digit_tdc[idig]-tfirst;
     if((tdcbin<0) || ((int)numBins<=tdcbin)) MSG("Chop",Msg::kDebug) << "Whups!" << endl;
     else if(digit_plane[idig]<=PlexPlaneId::LastPlaneNearCalor()) {
       energyVsTime[digit_tdc[idig]-tfirst] += sigcor; 
     }
   }

   // Used bins:
   std::vector<char> binsUsed(numBins,0);

   do {
     // Look for biggest peak. 
     UInt_t biggest_bin  = 99999;
     float  biggest_size = 0;
     for(UInt_t i=0;i<numBins;i++) {
       if(binsUsed[i]==0) 
         if(energyVsTime[i]>biggest_size) {
           biggest_size = energyVsTime[i];
           biggest_bin = i;
         }
     }
     
     if(biggest_bin==99999) break; // We've gone through all of them.
     if(biggest_size<100.) break; // We've hit 1 pe, rock bottom.
     
     // Collect the start and stop time for this chop.
     // Start 1 bin before the peak, and at least 1 bin after the peak.
     UInt_t bin_start = biggest_bin;
     UInt_t bin_stop =  biggest_bin;
     
     //for(Int_t i=-3;i<5;i++) {
     //  MSG("Chop",Msg::kDebug) << i << "\t" << energyVsTime[biggest_bin+i] << endl;
     // }
     
     if(binsUsed[bin_start-1]==0) bin_start--;
     if(binsUsed[bin_stop +1]==0) bin_stop++;
     
     bool done = false;
     while(!done) {
       // Stop at start of spill.
       if(bin_start==0) {
         //MSG("Chop",Msg::kDebug) << "Hit start of spill" << endl; 
         done=true;
       }
       
       if(ShouldSplit(energyVsTime[bin_start], 
                      energyVsTime[bin_start-1],
                      bin_start-1 - biggest_bin ) ) {
         //MSG("Chop",Msg::kDebug) << "Didn't move back.  dE = " << energyVsTime[bin_start-1] << "-" << energyVsTime[bin_start] << "  maxdE = " << max << endl;
         done = true;
       }
       
       // Stop if we've hit another chop.
       if(binsUsed[bin_start-1]) {
         //MSG("Chop",Msg::kDebug) << "Hit another chop." << endl;
         done = true;
       }
       
       if(!done) {
         //MSG("Chop",Msg::kDebug) << "Moving start back 1 notch. dE = " << climb << "  maxdE = " << max << endl;
         bin_start--;
       }
     };
     
     // Expand forwards until the energy starts climbing.
     // But, allow small pulses in for the first 5 buckets.
     done = false;
     while(!done) {
       // Stop at end of spill
       if(bin_stop >= numBins-1) {
         //MSG("Chop",Msg::kDebug) << "Hit end of spill" << endl;
         done = true;
       }
       
       // Allow 5 buckets worth of small stuff:
       if((energyVsTime[bin_stop+1] < 500.) && (bin_stop+1 < biggest_bin+7)) {
         // keep going
       } else {
         if(ShouldSplit(energyVsTime[bin_stop],
                        energyVsTime[bin_stop+1],
                        bin_stop+1 - biggest_bin) ) {
           //MSG("Chop",Msg::kDebug) << "Didn't move forward.  dE = " << energyVsTime[bin_stop+1] << "-" << energyVsTime[bin_stop] << "  maxdE = " << max << endl;
           done = true;
         }
       }
       
       // Stop if we hit another chop.
       if(binsUsed[bin_stop+1]) {
         //MSG("Chop",Msg::kDebug) << "Didn't move forward; hit another chop." << endl;
         done = true;
       }
       
       // If we're ok, increment and continue.
       if(!done) bin_stop++;
     }
     
     int tdc_start = bin_start+tfirst;
     int tdc_stop  = bin_stop+tfirst;
     
     // We have a chop candidate. Create it.
     DigitVector slc;

     for(UInt_t idig = 0; idig < ndigits; idig++ ) {
       int tdc = digit_tdc[idig];
       if((tdc>=(int)tdc_start)&&(tdc<=(int)tdc_stop))
         slc.push_back(digits[idig]);
     }
     cxx.SetDataIn(&(slc));
     CandDigitListHandle chopHandle = CandDigitList::MakeCandidate(ah,cxx);
     chopHandle.SetName(Form("Chop %d",nchop++));
     chopList.AddDaughterLink(chopHandle);
     
     MSG("Chop",Msg::kDebug) << "Creating chop. Big: " << biggest_bin
          << "  Start: " << bin_start << "   Stop: " << bin_stop 
          << "  Digits: " << slc.size() 
          << endl;

     

     // Zero out these buckets so they won't be caught again.
     for(UInt_t i=bin_start; i<=bin_stop; i++) {
       binsUsed[i] = 1;
     }
     
   } while(true); 
 
}

//______________________________________________________________________
void AlgChopListSharp::Trace(const char * /* c */) const
{
}

---