CellHitMerge.cxx

Go to the documentation of this file.

// $Id: CellHitMerge.cxx,v 1.2 2009/08/28 18:35:14 vahle Exp $
//
// A module to merge CellHits if necessary
// also removes "crazy" cells
// reads in vector of CellHits from Cal() Folder
// writes out new merged vector to a subfolder in Cal called MergeHits
// originally adapted from RecoSubShower::FormStripList
//
// \author $Author: vahle $
//
#include <iostream>
#include <vector>
#include <cmath>
#include "TH2D.h"
#include "TH1D.h"
#include "Geometry/Geometry.h"
#include "JobControl/Exception.h"
#include "Header/Header.h"
#include "RecoBase/CellHit.h"
#include "CellHitMerge/CellHitMerge.h"

using namespace std;
using namespace cellhitmerge;
using namespace recobase;

// This macro declares the existence of your module to ana
MODULE_DECL(CellHitMerge); //Required!

CellHitMerge::CellHitMerge(const char* version) :
  jobc::Module("cellhitmerge::CellHitMerge"), // Required!
  nprinted(0),
  fSkipMe(false),
  nevtmerged(0),
  maxdiffcells(0),
  npruned(0),
  maxpruned(0)
{
  delTMerged = new TH1D("delTMerged","",200,0,2000);
  delPHMerged = new TH1D("delPHMerged","",300,-1500,1500);
  delPhvdelT = new TH2D("delPhvdelT","",200,0,2000,300,-1500,1500);
  celldiff = new TH1D("celldiff","",31,-0.5,30.5);
  delPhvPh = new TH2D("delPhvPh","",150,0,1500,300,-1500,1500);
  phvph = new TH2D("phvph","",150,0,1500,150,0,1500);

  this->SetCfgVersion(version); // Required!
}

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

void CellHitMerge::Update(const cfg::Config& c)
{


  //maybe we wont want to do this step sometimes
  int domergestep;
  c("DoMergeStep").Get(domergestep);
  if(domergestep==0){
    fSkipMe=true;
  }
  else{
    fSkipMe=false;
  }

}

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

CellHitMerge::~CellHitMerge()
{
  cout<<"Number of events that got pruned: "<<npruned<<" max pruning: "<<maxpruned<<" cells"<<endl;
  cout<<"Number of events that got merged or pruned: "<<nevtmerged<<" max difference: "<<maxdiffcells<<" cells"<<endl;
}

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

jobc::Result CellHitMerge::Reco(edm::EventHandle& evt)
{
  if(fSkipMe) return jobc::kPassed;
  
  //get the list of CellHits
  std::vector<const recobase::CellHit*> hitlist(0);
  try{ evt.Cal().Get("Hits",hitlist); }
  catch(edm::Exception e){
    cerr<<"Couldn't find the CellHit list"<<endl;
    return jobc::kFailed;
  }

  //make a new hitlist to hold the pruned/merged hits
  std::vector<recobase::CellHit*> newhitlist;

  //now get a header
  std::vector<const hdr::Header*> head(0);
  try{ evt.DAQ().Get(".",head); }
  catch(edm::Exception e){
    cerr<<"Couldn't get a head"<<endl;
    return false;
  }
  geo::Geometry& geom = geo::Geometry::Instance(1,head[0]->DetectorID());
  int DETLASTPLANE = (int) geom.NPlanes();
  //now loop over CellHits in the list and look for crazies/duplicates
  std::vector<const recobase::CellHit*>::iterator chit=hitlist.begin();
  bool pruned=false;
  int ncellpruned = 0;
  while(chit!=hitlist.end()){
    const CellHit* stp = *chit;
    //    cout<<"Printing cell hit in CellHitMerge"<<endl;
    //    stp->Print();
    // if this cell has a plane number beyond the bounds of the detector
    //erase it from the list
    if(stp->Plane()<0 || stp->Plane()>DETLASTPLANE){
      hitlist.erase(chit); 
      cerr<<" CellHitMerge:: found plane crazy "<<stp->Plane()<<" strip "<<stp->Cell()<<endl;
      pruned=true;
      ncellpruned++;
      continue;
    }
    //if this cell has a cell number beyond the bounds of this plane
    //erase it from the list
    int PLANELASTCELL = (int) geom.Plane(stp->Plane()).Ncells();
    if(stp->Cell()<0 || stp->Cell()>PLANELASTCELL){
      hitlist.erase(chit);
      cerr<<" CellHitMerge:: found strip crazy "<< stp->Cell()<<" plane "<<stp->Plane()<< endl;       
      pruned=true;
      ncellpruned++;
      continue;
    }

    //make a new CellHit that we can modify if need be
    CellHit *newstp = new CellHit(*stp);

    //now check if this cell is duplicated in the list
    std::vector<const recobase::CellHit*>::iterator dhit=chit;
    dhit++;
    while(dhit!=hitlist.end()){
      const CellHit* stp2 = *dhit;
      if(newstp->Plane()==stp2->Plane() && newstp->Cell()==stp2->Cell()){
        if(nprinted<11){
          //      cerr<<"Merging Cells: "<<newstp->Plane()<<" "<<newstp->Cell()<<endl;
        }
        unsigned short tdc1, tdc2;
        newstp->TDC(tdc1);
        stp2->TDC(tdc2);
        float tdiff = fabs(1.*tdc1-1.*tdc2);
        delTMerged->Fill(tdiff);

        unsigned short adc1,adc2;
        newstp->ADC(adc1);
        stp2->ADC(adc2);
        float phdiff = 1.*adc1-1.*adc2;
        unsigned short maxph=adc1;
        if(adc2>adc1) maxph=adc2;
        delPHMerged->Fill(phdiff);
        delPhvdelT->Fill(tdiff,phdiff);
        delPhvPh->Fill(maxph,phdiff);
        phvph->Fill(adc1,adc2);
        newstp->Merge(stp2);
        hitlist.erase(dhit);
        continue;
      }      
      //      cout<<"about to put a CellHit in the newhitlist"<<endl;
      //      newstp->Print();
      dhit++;
    }
    newhitlist.push_back(newstp);
    chit++;
  }
  if(pruned){
    npruned++;
    if(ncellpruned>maxpruned){
      maxpruned=ncellpruned;
    }
  }
                      
  if(evt.Cal().GetFolder("./MergeHits")==0){
    evt.Cal().MakeFolder("./MergeHits");
  }
    
  //  cout<<"about to put vector of new hits in file "<<newhitlist.size()<<endl;
  evt.Cal().PutVector(newhitlist,"./MergeHits");

  return jobc::kPassed; // kFailed if you want to fail the event   
}

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

jobc::Result CellHitMerge::Ana(const edm::EventHandle& evt)
{
//======================================================================
// Called for every event. "Ana" implies that you are looking at the
// event, but not altering it.
//======================================================================


  std::vector<const recobase::CellHit*> hitlist(0);
  try{ evt.Cal().Get("Hits",hitlist); }
  catch(edm::Exception e){
    cerr<<"Couldn't find the Original hit list"<<endl;
    return jobc::kFailed;
  }

  std::vector<const recobase::CellHit*> mergehitlist(0);
  try{ evt.Cal().Get("./MergeHits",mergehitlist); }
  catch(edm::Exception e){
    cerr<<"Couldn't find the merged hit list"<<endl;
    return jobc::kFailed;
  }


  if(hitlist.size()!=mergehitlist.size()){
    nevtmerged++;
    int diff = hitlist.size()-mergehitlist.size();
    if(diff<0){
       diff = mergehitlist.size()-hitlist.size();
    }
    if(diff>maxdiffcells){
      maxdiffcells=diff;
    }
    celldiff->Fill(diff);
    //    cout<<"Hit List Sizes differ.  Printing both.  Size of original: "
    //  <<hitlist.size()<<" Size of Merged "<<mergehitlist.size()<<" Diff: "<<diff<<endl;
  }
  else{
    celldiff->Fill(0);
    return jobc::kPassed;
  }

  nprinted++;
  if(nprinted<11){
    cout<<"Will print 10 events that got merged.  On number: "<<nprinted<<endl;
    cout<<"Printing Original Hits: "<<hitlist.size()<<endl;
    for(unsigned int i=0;i<hitlist.size();i++){
      cout<<"** Printing Original CellHit index "<<i<<" **"<<endl;
      hitlist[i]->Print("");
      cout<<"*************"<<endl;
    }
    
    cout<<"Printing Merged Hits: "<<mergehitlist.size()<<endl;
    for(unsigned int i=0;i<mergehitlist.size();i++){
      cout<<"** Printing Merged CellHit index "<<i<<" **"<<endl;
      mergehitlist[i]->Print("");
      cout<<"*************"<<endl;
    }
  }
  return jobc::kPassed; // kFailed if you want to fail the event
  
}

---