AlgChopListMitre.cxx

Go to the documentation of this file.

// $Id: AlgChopListMitre.cxx,v 1.5 2007/11/11 08:26:13 rhatcher Exp $
//
// AlgChopListMitre.cxx
//

#include <cassert>
#include <cmath>

#include "CandChop/AlgChopListMitre.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 <TMath.h>

ClassImp(AlgChopListMitre)
CVSID( " $Id: AlgChopListMitre.cxx,v 1.5 2007/11/11 08:26:13 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.


//______________________________________________________________________
AlgChopListMitre::AlgChopListMitre()
{
}

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



//______________________________________________________________________
void AlgChopListMitre::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 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() != 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(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(s).

   std::vector< float >             profAngles;  
   profAngles.push_back(0);
   profAngles.push_back(1.0/60.); // 1 bucket per 60 planes.
   profAngles.push_back(-1.0/60.); // 1 bucket per 60 planes.
   UInt_t nProf = profAngles.size();

   std::vector< std::vector<float> >profiles(nProf);
   std::vector<float> meanPlane(numBins,0);


   for(UInt_t p = 0; p<nProf; p++) {
     profiles[p].resize(numBins,0.);
   }

   for(UInt_t idig = 0; idig < ndigits; idig++ ) {
     float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr);
     float tdcbin = digit_tdc[idig]-tfirst;
     float dplane = digit_plane[idig];
     if(dplane>PlexPlaneId::LastPlaneNearCalor()) 
       dplane = PlexPlaneId::LastPlaneNearCalor();
     dplane = dplane-60.;

     std::vector<int> proj(nProf);
     for(UInt_t p = 0; p<nProf; p++) 
       proj[p] = TMath::Nint(tdcbin + profAngles[p]*dplane);

     if((tdcbin<0) || ((int)numBins<=tdcbin)) MSG("Chop",Msg::kDebug) << "Whups!" << endl;
     else if(digit_plane[idig]<=PlexPlaneId::LastPlaneNearCalor()) {
       meanPlane[proj[0]] += dplane*sigcor;

       for(UInt_t p = 0; p<nProf; p++) {
         profiles[p][proj[p]] += sigcor;
       }
     }
   }

   for(UInt_t i=0;i<numBins;i++) meanPlane[i] /= profiles[0][i];
   

   std::vector<int> cuts;
   std::vector<int> cuts_type;

   // Insert a vertical cut at beginning of spill:
   cuts.push_back(0);
   cuts_type.push_back(0);

   
   int peak_bin = 0;
   bool rising = true;
   for(UInt_t bin=1;bin<numBins-1;bin++) {
     float dq = profiles[0][bin+1] - profiles[0][bin];
     float max_climb = 2.5*sqrt(fabs(profiles[0][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);
         cuts_type.push_back(0);
         rising = true;
       }
     } 

     if(rising) {
       if(dq< 0){
         rising = false;
         peak_bin = bin;
       }
     }
   }
      
   // We have first set of trial cuts. Now re-cut by the mitre saw.
   // Loop through cuts:
   for(UInt_t icut=0;icut<cuts.size(); icut++) {
     // find the lowest-charge bin of all the mitres, going back and forth +-1 bin
     float lowest_size = profiles[0][cuts[icut]];
     if(lowest_size>0) { // we can improve
       
       MSG("Chop",Msg::kDebug) << "Evaluating cut " << icut
            << " at bin " << cuts[icut] << " mean plane " << meanPlane[cuts[icut]] << endl;

       int a = cuts[icut]-5; if(a<0) a=0;
       int b = cuts[icut]+5; if(b>=(int)numBins) b=numBins-1;
       for(int i=a; i<=b; i++) {
         MSG("Chop",Msg::kDebug) << "Bin: " << i;
         for(UInt_t p = 0; p<nProf; p++)  MSG("Chop",Msg::kDebug) << "\t" << Form("%6.0f",profiles[p][i]);
         MSG("Chop",Msg::kDebug) << endl;
       }
           
       for(UInt_t p = 0; p<nProf; p++) {
         for(Int_t bin=cuts[icut]; bin<= cuts[icut]; bin++) {
           if(profiles[p][bin] < lowest_size) {
             lowest_size = profiles[p][bin];
             cuts[icut]      = bin;
             cuts_type[icut] = p;
           }
         }
       }

       MSG("Chop",Msg::kDebug) << "  -> Cutting bin " << cuts[icut] << " on projection " << cuts_type[icut] 
            << " lowbin: " << lowest_size << endl;
       
     }
   }


   // Insert vertical cut at end:
   cuts.push_back(numBins-1);
   cuts_type.push_back(0);
   
   for(UInt_t icut=0; icut<cuts.size()-1; icut++) {
     DigitVector chop;
     
     for(UInt_t idig = 0; idig < ndigits; idig++ ) {

       float tdcbin = digit_tdc[idig]-tfirst;
       float dplane = digit_plane[idig];
       if(dplane>PlexPlaneId::LastPlaneNearCalor()) 
         dplane = PlexPlaneId::LastPlaneNearCalor();
       dplane -= 60.;

       std::vector<int> proj(nProf);
       for(UInt_t p = 0; p<nProf; p++) 
         proj[p] = TMath::Nint(tdcbin + profAngles[p]*dplane);
       
       // Is it after start cut?
       if(proj[cuts_type[icut]] >= cuts[icut]) {
         
         // Is it after end cut?
         if(proj[cuts_type[icut+1]] < cuts[icut+1] ) {
           
           // Add it to the chop.
           chop.push_back(digits[idig]);
           
         } 
       }
     }
     
     cxx.SetDataIn(&(chop));
     CandDigitListHandle chopHandle = CandDigitList::MakeCandidate(ah,cxx);
     chopHandle.SetName(Form("Chop %d",nchop++));
     chopList.AddDaughterLink(chopHandle);
     
     MSG("Chop",Msg::kDebug) << "Creating chop. "
          << "  Start: " << cuts[icut] << "   Stop: " << cuts[icut+1]
          << "  Digits: " << chop.size() 
          << endl;
     
   } 
}

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

---