AlgChopListSharp2.cxx

Go to the documentation of this file.

// $Id: AlgChopListSharp2.cxx,v 1.3 2006/12/01 18:33:03 rhatcher Exp $
//
// AlgChopListSharp2.cxx
//

#include <cassert>
#include <cmath>

#include "CandChop/AlgChopListSharp2.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"

ClassImp(AlgChopListSharp2)
CVSID( " $Id: AlgChopListSharp2.cxx,v 1.3 2006/12/01 18:33:03 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(DetectorType::kNear,ElecType::kQIE,0,0,false,false);
const float k1pe = 100.; // About 1 pe.


//______________________________________________________________________
AlgChopListSharp2::AlgChopListSharp2()
{
}

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

class Chop : public DigitVector{
public:
  std::map<PlexStripEndId,float> fStripEnergy;
  std::map<int, float>           fPlaneEnergy;
  float                          fTotalEnergy;
  bool                           fDirty;

  Chop() : fDirty(false) {};

  void Add(CandDigitHandle& cdh) { 
    fDirty = true;
    push_back(cdh);
  }

  void BuildMaps() {
    if(fDirty) {
      fStripEnergy.clear();
      fPlaneEnergy.clear();
      fTotalEnergy = 0;
      for(UInt_t i=0;i<size();i++) {
        float sigcor = (*this)[i].GetCharge(CalDigitType::kSigCorr);
        PlexStripEndId seid = (*this)[i].GetPlexSEIdAltL().GetBestSEId();
        if(seid.IsValid()) {
          fStripEnergy[seid] += sigcor;
          fPlaneEnergy[seid.GetPlane()] += sigcor;
          fTotalEnergy += sigcor;
        }
      }
    }
    fDirty = false;
  }
};

//______________________________________________________________________
bool AlgChopListSharp2::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 AlgChopListSharp2::RunAlg(AlgConfig& /*algConfig*/, 
                           CandHandle &candHandle,  // thing to make
                           CandContext &candContext)
{

  // Config.
  bool cDoRecombobulation = false;

  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 AlgChopListSharp2 is not a digit list." << std::endl;
   }
   
   const CandDigitListHandle *cdlh_ptr = 
     dynamic_cast<const CandDigitListHandle*>(candContext.GetDataIn());
   
   const MomNavigator *mom = candContext.GetMom();
   RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
   if (!rr) {
     MSG("Chop", Msg::kWarning) << "No RawRecord in MOM." << endl;
     return;
   }
   const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *>
     (rr->FindRawBlock("RawDigitDataBlock"));
   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() != DetectorType::kNear) 
     MSG("Chop",Msg::kError) << "Running the Sharp2 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_Sharp2" << 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; 
     }
   }


   
   std::vector<int> cuts;
   /*
   cuts.push_back(0); // cut at start.   
   int peak_bin = 0;
   bool rising = true;
   for(UInt_t bin=0;bin<numBins-1;bin++) {
     float dq = energyVsTime[bin+1] - energyVsTime[bin];
     float max_climb = 2.5*sqrt(fabs(energyVsTime[bin])/k1pe)*k1pe;
     if(max_climb < 4.*k1pe) max_climb = max_climb*2.;

     bool falling = !rising;

     if(falling && ((bin-peak_bin)<8)) 
       if(max_climb<5.*k1pe) max_climb = 5.*k1pe;

     if(falling) {
       if(dq > max_climb){
         cuts.push_back(bin);
         rising = true;
       }
     } 

     if(rising) {
       if(dq< 0){
         rising = false;
         peak_bin = bin;
       }
     }
   }
   
   cuts.push_back(numBins-1); // cut at the end.
   */

   // Used bins:
   std::vector<char> binsUsed(numBins,0);
   binsUsed[0]=1;
   binsUsed[numBins-1]=1;

   cuts.push_back(0);
   cuts.push_back(numBins-1);

   

   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<=0.) break; // We've hit 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;
     
     if(binsUsed[bin_start-1]==0) bin_start--;
     if(binsUsed[bin_stop +1]==0) bin_stop++;
     
     bool done = false;
     while(!done) {

       // Stop if there's a rise. If so, mark as a contested bin.
       if(ShouldSplit(energyVsTime[bin_start], 
                      energyVsTime[bin_start-1],
                      bin_start-1 - biggest_bin ) ) {
         done = true;
         cuts.push_back(bin_start);
       }
       
       // Stop if we've hit another chop.
       if(binsUsed[bin_start-1]) {
         done = true;
       }
       
       if(!done) {
         bin_start--;
       }
     };
     
     // Expand forwards until the energy starts climbing.
     // But, allow small pulses in for the first 5 buckets.
     done = false;
     while(!done) {
       
       // 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) ) {
           done = true;
           cuts.push_back(bin_stop); // mark this one as contested.
         }
       }
       
       // 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++;
     }
     
     // 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); 


   // Time-order the cuts.
   std::sort(cuts.begin(), cuts.end());

   // One chop for every cut, corresponding to the cut that ends the chop.
   std::vector<Chop> chops(cuts.size());

   // Build initial cuts, excluding hits in contested bins.
   // loop cuts, 1 to n
   for(UInt_t icut = 1; icut<= cuts.size(); icut++) {
     int bin_start = cuts[icut-1]+1; // Excluding last contested bin
     int bin_stop  = cuts[icut]-1;   // Excluding current contested bin.

     int tdc_start = bin_start + tfirst;
     int tdc_stop  = bin_stop  + tfirst;
     
     for(UInt_t idig = 0; idig < ndigits; idig++ ) {
       int tdc = digit_tdc[idig];
       if((tdc>=(int)tdc_start)&&(tdc<=(int)tdc_stop)) {
         chops[icut].Add(digits[idig]);
       }
     }     
   }

   // Now loop through again and deal with any contested bins.
   // Ignore the first and last cuts, since they are at the start
   // and end of the event.
   for(UInt_t icut = 1; icut< cuts.size()-1; icut++) {
     int contested_bin = cuts[icut];
     int contested_tdc = contested_bin + tfirst;
     
     // Make a list of digits in the contested bin.
     DigitVector contested_digits;
     for(UInt_t idig = 0; idig < ndigits; idig++ ) {
       int tdc = digit_tdc[idig];
       if(tdc==contested_tdc) {
         contested_digits.push_back(digits[idig]);
       }
     }

     if(contested_digits.size()>0) {
       // only bother if there's actually some digits in the contested bin.

       int chop1 = icut;
       int chop2 = icut+1;

       // Tell the fighting chops to build maps.
       chops[chop1].BuildMaps();
       chops[chop2].  BuildMaps();

       MSG("Chop",Msg::kDebug) << "Dealing with contested bin " << cuts[icut] 
            << " with " << contested_digits.size() << " contested digits." << endl;

       int assign_strip1 = 0;
       int assign_strip2 = 0;
       int assign_plane1 = 0;
       int assign_plane2 = 0;
       int assign_total1 = 0;
       int assign_total2 = 0;

       
       // Now loop through contested digits and assign them:
       for(UInt_t ic=0; ic<contested_digits.size(); ic++) {
         // Assign to chop with most energy on same strip:
         CandDigitHandle digit = contested_digits[ic];
         PlexStripEndId seid =digit.GetPlexSEIdAltL().GetBestSEId();
         float chop1stripE = chops[chop1].fStripEnergy[seid];
         float chop2stripE  =chops[chop2  ].fStripEnergy[seid];

         /*
         MSG("Chop",Msg::kDebug) << "  Consider " << seid.AsString() << " q=" << digit.GetCharge(CalDigitType::kSigCorr)
              << "  chop1(s/p/t) = " 
              << chops[chop1].fStripEnergy[seid] << "/"
              << chops[chop1].fPlaneEnergy[seid.GetPlane()] << "/"
              << chops[chop1].fTotalEnergy 
              << "  chop2(s/p/t) = " 
              << chops[chop2].fStripEnergy[seid] << "/"
              << chops[chop2].fPlaneEnergy[seid.GetPlane()] << "/"
              << chops[chop2].fTotalEnergy 
              << endl;
         */

         if(chop1stripE > chop2stripE){
           chops[chop1].Add(digit);
           assign_strip1++;
         } else if(chop2stripE > chop1stripE){
           chops[chop2].Add(digit);
           assign_strip2++;
         } else {
           // Darn. The strip content is the same in both. (probably 0)
           // So compare on the plane level:
           
           float chop1planeE = chops[chop1].fPlaneEnergy[seid.GetPlane()];
           float chop2planeE = chops[chop2].fPlaneEnergy[seid.GetPlane()];
           if(chop1planeE > chop2planeE){
             chops[chop1].Add(digit);
             assign_plane1++;
           } else  if(chop2planeE > chop1planeE){
             chops[chop2].Add(digit);
             assign_plane2++;      
           } else {

             // Ok, so there's nothing in either plane. Assign it to the biggest chop,
             // or the earlier one if tied.
             if(chops[chop1].fTotalEnergy>=chops[chop2].fTotalEnergy) {
               chops[chop1].Add(digit);
               assign_total1++;
             } else {
               chops[chop2].Add(digit);
               assign_total2++;
             }
           }
         }
       } //uuuuuuugly!
       
      
       MSG("Chop",Msg::kDebug) << "  Assigned by strip: " << assign_strip1 << "/" << assign_strip2
            << "   by plane: "         << assign_plane1 << "/" << assign_plane2
            << "   by total: "         << assign_total1 << "/" << assign_total2
            << endl;
     }
    

   }
   
   // And now...
   // Recombobulate!
   if(cDoRecombobulation) {

   for(UInt_t ichop=1;ichop<chops.size();ichop++) chops[ichop].BuildMaps();
   
   for(UInt_t icut=1;icut<cuts.size()-1;icut++) {
     for(int ilr = 0; ilr<2; ilr++) { //loop left and right
       int contested_tdc;
       int chopfrom;
       int chopto;
       if(ilr==0) {
         contested_tdc = cuts[icut]-1 + tfirst;
         chopfrom      = icut;
         chopto        = icut+1;
       } else {
         contested_tdc = cuts[icut]+1 + tfirst;
         chopfrom      = icut+1;
         chopto        = icut;
       }

       int moved = 0;
       
       for(UInt_t idig = 0; idig<chops[chopfrom].size(); idig++) {
         CandDigitHandle digit = chops[chopfrom][idig];
         int tdc = (cal.GetTDCFromTime(digit.GetTime(CalTimeType::kNone), kQieRcid));
         if(tdc == contested_tdc) {
           PlexStripEndId seid =digit.GetPlexSEIdAltL().GetBestSEId();
           float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr);
           
           // If the chop that owns the digit has no other digits on that strip:
           if(chops[chopfrom].fStripEnergy[seid] <= sigcor) {
             // If the neigboring chop HAS got some energy on that strip
             if(chops[chopto].fStripEnergy[seid] > 0) {
               // Then move the digit.
               chops[chopto].push_back(digit);
               chops[chopfrom].erase(chops[chopfrom].begin() + idig);
               idig--; // set so we look at the next element instead of skipping one.
               moved++;
             }
           }
         }
       }
       
       MSG("Chop",Msg::kDebug) << "Recombobulated " << moved << " digits in bin " << contested_tdc-tfirst
            << " from chop " << chopfrom << " to chop " << chopto << endl;
     } // End left/right loop
   } // End loop over cuts

   } // end Recombobulation


   // Finally, store the chops.
   
   for(UInt_t ichop=1;ichop<chops.size();ichop++) {

     if(chops[ichop].size()>0) {
   
       cxx.SetDataIn(&(chops[ichop]));
       CandDigitListHandle chopHandle = CandDigitList::MakeCandidate(ah,cxx);
       chopHandle.SetName(Form("Chop %d",ichop));
       chopList.AddDaughterLink(chopHandle);
       
       MSG("Chop",Msg::kDebug) << "Creating chop. "
            << "  Start: " << cuts[ichop-1] << "   Stop: " << cuts[ichop] 
            << "  Digits: " << chops[ichop].size() 
            << endl;

       nchop++;
     }
   } 
}

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

---