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/Cluster.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"
#include "Utilities/EDMUtils.h"

using namespace recobase;
using namespace cluster;

MODULE_DECL(Clusterer);
ClassImp(Clusterer)

Clusterer::Clusterer(const char* version) : jobc::Module("Clusterer")
{
 fGeo=0;
  fClusterAlgType=0;
  this->SetCfgVersion(version);

}


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.Cal().Get("Hits",cellhit);
  }catch(...){return jobc::kFailed;}
  

  if (!fGeo) fGeo = &util::EDMUtils::GetGeometry(evt);
  
  evt.Reco().List();

  std::vector<recobase::Cluster> 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++)
    {
      UpdateClusterStat(planecluster[i]);
      //                printf("cluster at %f %f w %f d %f\n",planecluster[i].Zpos(),planecluster[i].Tpos(),planecluster[i].fdt,planecluster[i].fdz);
        if(planecluster[i].NCell())
      evt.Reco().Put(planecluster[i],"");
    }

  evt.Reco().List();

  return jobc::kPassed;
}


void Clusterer::UpdateClusterStat(recobase::Cluster & p)
{
 

  double fQtot=0.0;
  double tmax=-1.E10, tmin=1.E10;
  double ftpos=0.0;
  double fzpos=0.0;

  double fdt=0;
  double fdz=0;

  if(!fGeo)
  {printf("NO GEO!!!\n");return;}



  int view = fGeo->Plane(p.Plane(0)).View();

        

  int cells = (view == geo::kX) ? p.NXCell() : p.NYCell();

  for ( int i=0;i<cells;i++){
    const recobase::CellHit* ch = (view == geo::kX) ? p.XCell(i) : p.YCell(i);


    float mip=0.0;
    if(ch->PE() > 0) mip = ch->PE();else continue;
    fQtot +=mip;

    double tpos=0;
    double zpos=0;



    unsigned int cell = ch->Cell();
    unsigned int plane = ch->Plane();
    //printf("%d %d \n",plane,cell);
    
    double posit[3];

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

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

    if(i==0)
      { 
        fdt=fGeo->Plane(plane).Cell(cell).HalfD();
        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;

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

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

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

  p.SetQ(fQtot);

  double pos[4];
  double dpos[4];
  for(int i=0;i<4;i++)
  {
    pos[i]=0;
    dpos[i]=0;
  }

  
  pos[2]=fzpos;
  dpos[2]=fdz;

  if(view==geo::kX)
  {
    pos[0]=ftpos;
    dpos[0]=fdt;
  }else{
    pos[1]=ftpos;
    dpos[1]=fdt;
  }

  p.SetXYZT(pos);
  p.SetdXYZT(dpos);


  //printf("\n");


}

---