CalHit.cxx

Go to the documentation of this file.

#include <sstream>
#include "RawData/RawDigit.h"
#include "RecoBase/CellHit.h"
#include "CalHit/CalHit.h"
#include "JobControl/Exception.h"
#include "Header/Header.h"
#include "Utilities/EDMUtils.h"

using namespace calhit;

MODULE_DECL(CalHit);

CalHit::CalHit(const char* version) : 
  jobc::Module("CalHit"),
  fCal(0),
  fGeo(0),
  fRawDigitInputFolder("."),
  fCellHitOutputFolder("Hits"),
  fRecoHitOutputFolder("Rough"),
  fRunRoughReco(0)
{
  this->SetCfgVersion(version);
}

//......................................................................
CalHit::~CalHit() { }

//......................................................................

void CalHit::Update(const cfg::Config& c)
{
  c("RawDigitInputFolder").Get(fRawDigitInputFolder);
  c("CellHitOutputFolder").Get(fCellHitOutputFolder);
  c("RecoHitOutputFolder").Get(fRecoHitOutputFolder);
  c("RunRoughReco")       .Get(fRunRoughReco);
}

//......................................................................

jobc::Result CalHit::Reco(edm::EventHandle& evt)
{

  // Load the RawDigit list
  std::vector<const rawdata::RawDigit*> digitlist(0);
  try {
    evt.Raw().Get(fRawDigitInputFolder.c_str(),digitlist);
  } 
  catch (edm::Exception e) {
    return jobc::kFailed;
  }

  // Get a Calibrator instance for this event
  if (!fCal) fCal = &util::EDMUtils::GetCalibrator(evt);

  // Get a Geometry instance for this event, if needed (RoughReco)
  if (fRunRoughReco&&!fGeo) fGeo = &util::EDMUtils::GetGeometry(evt);

  // Make a CellHit for every RawDigit and write it out
  UInt_t Ndigits = digitlist.size();

  if (Ndigits) {

    // Make the output folder
    if(evt.Cal().GetFolder(fCellHitOutputFolder.c_str()) == 0)
      evt.Cal().MakeFolder(fCellHitOutputFolder.c_str());

    // Make the hits and write them out
    for (UInt_t iD=0;iD<Ndigits;iD++) {

      CellHit curr_cellhit;
      fCal->MakeCellHit(digitlist[iD],&curr_cellhit);

      evt.Cal().Put(curr_cellhit,fCellHitOutputFolder.c_str());

    }

    // Should we make rough RecoHits based on neighboring planes' hits?
    if (fRunRoughReco) RoughReco(evt);

  }


  return jobc::kPassed;

}


//-------
// RoughReco probably shouldn't be in CalHit.  But, it's more of an
// example of how to use the RecoHit class.  I don't expect any
// physics to necessarily get done with the output of RoughReco.
//
// This code is simply an adaptation of P. Shanahan's pre-Oct-2009
// version of CalHit.cxx.

void CalHit::RoughReco(edm::EventHandle& evt) {

  // Load the CellHit list (just made)
  std::vector<const recobase::CellHit*> cellhitlist(0);
  evt.Cal().Get(fCellHitOutputFolder.c_str(),cellhitlist);

  // Make the output folder
  if(evt.Reco().GetFolder(fRecoHitOutputFolder.c_str()) == 0)
    evt.Reco().MakeFolder(fRecoHitOutputFolder.c_str());

  // to handle mean positions
  std::map<unsigned int, std::pair<double, double> > planeMean;
  planeMean.clear();

  // Get Q-weighted mean positions in each plane
  UInt_t Ncellhits = cellhitlist.size();

  for (UInt_t iC=0;iC<Ncellhits;iC++) {
    const CellHit *cellhit = cellhitlist[iC];
    UInt_t plane = cellhit->Plane();
    UInt_t cell  = cellhit->Cell();

    Double_t xyz[3];
    fGeo->Plane(plane).Cell(cell).GetCenter(xyz);
    Double_t viewpos = 0;
    if      (cellhit->View()==geo::kX) viewpos = xyz[0];
    else if (cellhit->View()==geo::kY) viewpos = xyz[1];

    if (planeMean.find(plane)==planeMean.end()) {
      // init
      planeMean[plane].first  = 0;
      planeMean[plane].second = 0;
    }
    planeMean[plane].first  += cellhit->PECorr();
    planeMean[plane].second += cellhit->PECorr()*viewpos;
  }
  
  // Determine ortho pos for each CellHit, and make the corresponding RecoHit
  for (UInt_t iC=0;iC<Ncellhits;iC++) {
    const CellHit *cellhit = cellhitlist[iC];
    UInt_t plane = cellhit->Plane();

    // first, just accumulate existing plane numbers in other view,
    // within 3 planes of hit planes
    // use a fixed array so we can skip specific elements later if necessary
    unsigned int planes[6]={9999,9999,9999,9999,9999,9999};

    // attempt to walk back three planes.  
    unsigned int pl=fGeo->NextPlaneOtherView(plane,-1);
    if (pl!=geo::kPLANE_NOT_FOUND) {
      planes[0]=pl;
      pl=fGeo->NextPlaneInView(pl,-1);// have to protect against Plane
      // not found, since otherwise the second
      // try will give the end of the detector
      if (pl!=geo::kPLANE_NOT_FOUND) { 
        planes[1]=pl;
        pl=fGeo->NextPlaneInView(pl,-1);
        if (pl!=geo::kPLANE_NOT_FOUND)  planes[2]=pl;
      }
    }

    // and forward three planes
    pl=fGeo->NextPlaneOtherView(plane,1);
    if (pl!=geo::kPLANE_NOT_FOUND) {
      planes[3]=pl;
      pl=fGeo->NextPlaneInView(pl,1);// have to protect against Plane
      // not found, since otherwise the second
      // try will give plane 1 (I think)
      if (pl!=geo::kPLANE_NOT_FOUND) {
        planes[4]=pl;
        pl=fGeo->NextPlaneInView(pl,1);
        if (pl!=geo::kPLANE_NOT_FOUND)  planes[5]=pl;
      }
    }

    double sPh=0; // sum of pulseheight and (ph-weighted) position 
    double sTPos=0;

    for (unsigned int ip=0; ip<6; ip++) {
      // only uses the -3 and +3 planes if we don't already have signal
      // from +/-1 and 2.
      if (sPh>0 && (ip==2 || ip==5) ) continue;
      
      unsigned int planeno=planes[ip];
      if (planeMean.find(planeno)!=planeMean.end()) {//9999 won't be found...
        sPh+=planeMean[planeno].first;
        sTPos+=planeMean[planeno].second;
      } 
    }

    // use middle of cell if position not calculable
    Double_t w = 0;
    if (sPh>0) w = sTPos/sPh;

    // Okay, got the ortho position (w).  Make the RecoHit.
    RecoHit recohit;
    fCal->MakeRecoHit(cellhit,&recohit,w);
    evt.Reco().Put(recohit,fRecoHitOutputFolder.c_str());

  } // hit loop
  
}

---