ChopHelper.cxx

Go to the documentation of this file.

// $Id: ChopHelper.cxx,v 1.5 2010/01/06 18:34:37 rhatcher Exp $
//
// ChopHelper.cxx
//

#include <cassert>
#include <cmath>

#include "ChopHelper.h"
#include "ChopHelp.h"
#include "DigitVector.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(ChopHelper)

CVSID( " $Id: ChopHelper.cxx,v 1.5 2010/01/06 18:34:37 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.


//______________________________________________________________________
ChopHelper::ChopHelper(const MomNavigator* mom) : fMom(mom)
{
}

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

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 ChopHelper::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;    
}


//______________________________________________________________________
ChopHelp* ChopHelper::GetChopHelp( const CandHandle* ev1, 
                                   const CandHandle* ev2,
                                   const CandHandle* ev3,
                                   const CandHandle* ev4,
                                   const CandHandle* ev5 )
{
  std::vector<CandHandle> v;
  if(ev1) v.push_back(*ev1);
  if(ev2) v.push_back(*ev2);
  if(ev3) v.push_back(*ev3);
  if(ev4) v.push_back(*ev4);
  if(ev5) v.push_back(*ev5);
  
  return GetChopHelp(v);
}

ChopHelp* ChopHelper::GetChopHelp( std::vector<CandHandle> events )
{
 

  // Config.
   bool cDoRecombobulation = false;

   if (events.size()==0) return 0;
  
   const RawDigitDataBlock* rddb = DataUtil::GetRawBlock<RawDigitDataBlock>(fMom);
   if (!rddb) {
     MSG("Chop", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
     return 0;
   }
   
   const VldContext &context = *(events[0].GetVldContext());
   if(context.GetDetector() != Detector::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;
   std::vector<DigitVector> eventDigits(events.size());
   for(UInt_t i=0;i<events.size();i++) {
     eventDigits[i].FillFrom(&(events[i]));
     digits.FillFrom(eventDigits[i],0,eventDigits[i].size(), true);
   }

   UInt_t ndigits = digits.size();
    
   // 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;

   // 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
   

   // Remove empty chops.
   
   for(UInt_t i=0;i<chops.size();i++) {
     if(chops[i].size()==0) {
       chops.erase(chops.begin()+i);
       i--;
     }
   }

   // Finally, evaluate the chops vs the events.

   ChopHelp* chophelp = new ChopHelp(chops.size(), events.size());

   for(UInt_t ichop=0;ichop<chops.size();ichop++) {
     double tot = 0;
     for(UInt_t idig=0;idig<chops[ichop].size(); idig++) {
       tot += chops[ichop][idig].GetCharge(CalDigitType::kSigCorr);
     }
     chophelp->chopph[ichop]=tot;
   }

   for(UInt_t ievt=0;ievt<events.size();ievt++) {
     double tot = 0;
     for(UInt_t idig=0;idig<eventDigits[ievt].size(); idig++) {
       tot += eventDigits[ievt][idig].GetCharge(CalDigitType::kSigCorr);
     }
     chophelp->candph[ievt]=tot;
   }

   // Fill correlation matrix.

   for(UInt_t ievt=0;ievt<events.size();ievt++) {
     for(UInt_t ichop=0;ichop<chops.size();ichop++) {
       double tot = 0;
       for(UInt_t idig=0;idig<eventDigits[ievt].size(); idig++) {
         CandDigitHandle dig = eventDigits[ievt][idig];
         for(UInt_t idig2=0;idig2<chops[ichop].size(); idig2++) {
           if(dig == chops[ichop][idig2]) 
             tot+=dig.GetCharge(CalDigitType::kSigCorr);
         }
       }
       chophelp->matchmatrix[ievt*chops.size() + ichop] = tot;
     }
   } // uuuuugly again



   // Now find the best matches and fill the eff/pur values.
   for(UInt_t ievt=0; ievt<events.size(); ievt++) {
     int bestchop=0;
     Float_t bestmatch=-1;
     for(UInt_t ichop=0;ichop<chops.size();ichop++) {
       if(chophelp->matchmatrix[ievt*chops.size() + ichop] > bestmatch) {
         bestmatch = chophelp->matchmatrix[ievt*chops.size() + ichop];
         bestchop=ichop;
       }
     }
     chophelp->chopmatch[ievt] = bestchop;

     chophelp->estpurity[ievt] = chophelp->matchmatrix[ievt*chops.size() + bestchop] / chophelp->candph[ievt];
     chophelp->estcompleteness[ievt] = chophelp->matchmatrix[ievt*chops.size() + bestchop] / chophelp->chopph[bestchop];
   }
   
   
   return chophelp;
}

---