AlgChopListPerfectMC.cxx

Go to the documentation of this file.

// $Id: AlgChopListPerfectMC.cxx,v 1.3 2007/11/11 08:26:13 rhatcher Exp $
//
// AlgChopListPerfectMC.cxx
//
// A perfect slicer, that cheats by looking at MC truth info
//

#include <cassert>

#include "CandChop/AlgChopListPerfectMC.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/Truthifier.h"

ClassImp(AlgChopListPerfectMC)
CVSID( " $Id: AlgChopListPerfectMC.cxx,v 1.3 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);


//______________________________________________________________________
AlgChopListPerfectMC::AlgChopListPerfectMC()
{

}

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



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

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

   assert(candContext.GetDataIn());
   // Check for CandDigitListHandle input
   if (!(candContext.GetDataIn()->InheritsFrom("CandDigitListHandle"))) {
     MSG("Chop",Msg::kWarning ) << "Data into AlgChopListPerfectMC 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 PerfectMC 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 nslice = 0;

   const Truthifier& truth = Truthifier::Instance(candContext.GetMom());
   std::vector<Int_t> neutrinos = truth.GetListOfNeutrinoIndecies();

   for(UInt_t inu=0; inu < neutrinos.size(); inu++) {
     // For each true neutrino in the event, make a slice:
     DigitVector slc;

     // Loop through digits and put the relevant ones in.
     for(UInt_t idig = 0; idig<ndigits; idig++) {
       float frac = truth.IsDigitFromNeutrino(neutrinos[inu],digits[idig]);
       if(frac>0.20) // Include the digit if more than 20% of it belongs to this neutrino.
         slc.push_back(digits[idig]);
     }
     if(slc.size()>0) { // there was at least one digit for this neutrino:
       // Store the slice.
       cxx.SetDataIn(&(slc));
       CandDigitListHandle sliceHandle = CandDigitList::MakeCandidate(ah,cxx);
       sliceHandle.SetName(Form("Chop %d",nslice++));
       sliceHandle.SetTitle(Form("Neutrino %d StdHepIndex %d",inu,neutrinos[inu]));
       sliceList.AddDaughterLink(sliceHandle);       
     }
   }

}

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

---