Clusterer.cxx

Go to the documentation of this file.

#include<cmath>
#include "Clusterer/Clusterer.h"
#include "Clusterer/HitCluster.h"
#include "Clusterer/HitClusterFast.h"
#include "RecoBase/PlaneCluster.h"
#include "RecoBase/CellHit.h"
#include "JobControl/Exception.h"

#include "RawData/RawDigit.h"
#include "RawData/DAQHeader.h"
#include "Geometry/Geometry.h"
#include "Geometry/PlaneGeo.h"
#include "Geometry/CellGeo.h"
using namespace recobase;
using namespace cluster;

MODULE_DECL(Clusterer);
ClassImp(Clusterer)

Clusterer::Clusterer(const char* version) : jobc::Module("Clusterer")
{
  fClusterAlgType=0;
  this->SetWatch("ClustererConfig","default");
}


void Clusterer::Update(const cfg::Config& c)
{
  c("ClusterAlgType").Get(fClusterAlgType);
}

Clusterer::~Clusterer()
{
  this->RemoveAllWatches();
}

jobc::Result Clusterer::Reco(edm::EventHandle& evt)
{
  //make cell hits from rawdigits
  /*
    std::vector<const rawdata::RawDigit*> rawdigits(0);
    evt.Raw().Get("./",rawdigits);

    printf("cluster-\n");
    for (unsigned int i =0;i<rawdigits.size();i++)
    {
    printf("%d %d \n",i,rawdigits[i]->TDC(0));
    CellHit c(rawdigits[i]->Channel(),rawdigits[i]->ADC(0),rawdigits[i]->TDC(0));
    c.SetPlane(i/100);
    c.SetCell(i%100);
    evt.Reco().Put(c,"");
    }
    printf("-cluster\n");

  */



  //get cellhits from the file
  std::vector<const recobase::CellHit*> cellhit(0);
  //assert_jobc(evt.Reco().Get("Hits",cellhit),"No CellHits found in DetSim() folder!");

  try{
    evt.Reco().Get("Hits",cellhit);
  }catch(...){return jobc::kFailed;}
  

  std::vector<const rawdata::DAQHeader*> header;
  short int det = rawdata::kFar;
  try{ evt.DAQ().Get("./",header); }
  catch(edm::Exception e){
    std::cerr << "Error retrieving daq header, while looking "<<
      "in TrackReco::Reco(), using default fardet geometry" << std::endl;
  }
  if(header.size() > 0) det = header[0]->DetId();
  //  if (!fGeo) // ideally we would use DAQHeader instead of fDetGeom...
  fGeo = &geo::Geometry::Instance(evt.Header().Run(), det);
  
  evt.Reco().List();

  std::vector<recobase::PlaneCluster> planecluster(0);

  if (fClusterAlgType==0)
    {
      HitCluster a;
      planecluster = a.MakeClusters(cellhit);
    }
  if (fClusterAlgType==1)
    {
      HitClusterFast a;
      planecluster = a.MakeClusters(cellhit);
    }


  for(unsigned int i=0; i<planecluster.size(); i++)
    {
      UpdatePlaneClusterStat(planecluster[i]);
      //                printf("cluster at %f %f w %f d %f\n",planecluster[i].Zpos(),planecluster[i].Tpos(),planecluster[i].fdt,planecluster[i].fdz);

      evt.Reco().Put(planecluster[i],"");
    }

  evt.Reco().List();

  return jobc::kPassed;
}


void Clusterer::UpdatePlaneClusterStat(recobase::PlaneCluster & p)
{
  //printf("\n");

  float fQtot=0.0;
  float tmax=-1.E10, tmin=1.E10;
  p.ftpos=0.0;
  p.fzpos=0.0;
  for (unsigned int i=0;i<p.fCellHit.size();i++){
    float mip;
    float f;
    if(p.fCellHit[i]->PE(f))mip=f;else continue;
    fQtot +=mip;

    double tpos=0;
    double zpos=0;

    unsigned int cell = p.fCellHit[i]->Cell();
    unsigned int plane = p.fCellHit[i]->Plane();
    //printf("%d %d \n",plane,cell);
    //const geo::CellGeo& gcell= fGeo->Plane(plane).Cell(cell);
    double posit[3];

    fGeo->Plane(plane).Cell(cell).GetCenter(posit,0.0);

    //fGeo->Plane(plane).Cell(cell).LocalToWorld(posit,posit);

    if(i==0)
      { 
        p.fdt=fGeo->Plane(plane).Cell(cell).HalfD();
        p.fdz=fGeo->Plane(plane).Cell(cell).HalfW();

      }
    //  double posit[3];
    //  gcell.GetCenter(posit,0.0);

    //printf("   (%d,%d,%f,%f,%f,%f)",plane,cell,posit[0],posit[1],posit[2],mip);


    if(fGeo->Plane(plane).View()==geo::kY)
      tpos = posit[1];
    else
      tpos=posit[0];

    zpos = posit[2];
     
    //  tpos=zpos=1;

    p.ftpos +=tpos*mip;
    p.fzpos +=zpos*mip;

    if( tpos > tmax ) tmax = tpos;
    if( tpos < tmin ) tmin = tpos;

  }
  p.ftpos/=fQtot;
  p.fzpos/=fQtot;
  p.fdt = sqrt(pow(p.fdt,2.) + pow((tmax-tmin)/2.,2.));

  p.fW=fQtot;


  //printf("\n");


}

---