ClusterCheck.cxx

Go to the documentation of this file.

// $Id: ClusterCheck.cxx,v 1.1.1.1 2009/08/21 19:14:29 vahle Exp $
//
// A Module to make clusters within 1 view (2d clusters)
// Based on RecoSubShower code 
// A Module to make clusters within 1 view (2d clusters)
// Based on RecoSubShower::FindCluster and RecoSubShower::MergeCluster code 
//
// \author $Author: vahle $
//
#include <iostream>
#include <cmath>
#include "TGraphErrors.h"
#include "TF1.h"
#include "TH2D.h"
#include "TH1D.h"
#include "TMath.h"
#include "Geometry/Geometry.h"
#include "Geometry/PlaneGeo.h"
#include "Geometry/CellGeo.h"
#include "Header/Header.h"
#include "Simulation/MCTruth.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/FLSHit.h"
#include "JobControl/Exception.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/ShowerBase.h"
#include "ClusterMakerSS/ClusterSSHelpers.h"
#include "ClusterMakerSS/MakeClusterSS.h"
#include "ClusterMakerSS/ClusterCheck.h"

using namespace std;
using namespace clusterss;

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

ClusterCheck::ClusterCheck(const char* version) :
  jobc::Module("makeclusterss::ClusterCheck")
{
  // Set the configuration version tag
  this->SetCfgVersion(version); // Required!
}

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

void ClusterCheck::Update(const cfg::Config& /*c*/)
{
}

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

ClusterCheck::~ClusterCheck()
{
}

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

jobc::Result ClusterCheck::Reco(edm::EventHandle& /*evt*/)
{
  
  return jobc::kPassed; // kFailed if you want to fail the event
}

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

jobc::Result ClusterCheck::Ana(const edm::EventHandle& evt)
{
  const edm::EventHeader &eh = evt.Header();
  int run = eh.Run();
  int event = eh.Event();
  cout<<"*************** On run "<<run<<" event "<<event<<endl;

  //now get a daq header and geometry
  std::vector<const hdr::Header*> head(0);
  try{ evt.DAQ().Get(".",head); }
  catch(edm::Exception e){
    cerr<<"Couldn't get a head"<<endl;
    return jobc::kFailed;
  }
  geo::Geometry &geom = geo::Geometry::Instance(1,head[0]->DetectorID());


  //get the list of CellHits
  std::vector<const recobase::CellHit*> hitlist(0);
  try{ evt.Cal().Get("./MergeHits",hitlist); }
  catch(edm::Exception e){
    cerr<<"Couldn't find the MergedCellHit list, trying for the original"<<endl;
    try{ evt.Cal().Get("./Hits",hitlist); }
    catch(edm::Exception e){
      cerr<<"Couldn't find the Original CellHit list either, failing."<<endl;
      return jobc::kFailed;
    }
  }
  if(hitlist.size()<5){
    //    cerr<<"Not enough activity"<<endl;
    //    return jobc::kFailed; //lets just look at things that have some real activity
  }

  std::vector<const sim::MCTruth*> truelist(0);
  try { evt.MC().Get(".", truelist); }
  catch(edm::Exception e){
    cerr<<"Couldn't get a truelist"<<endl;
    return jobc::kFailed;
  }

  //get some truthiness
  int tidx=0;
  const TParticle &nupart = truelist[tidx]->GetNeutrino()->Nu();
  double eenu = 1.*nupart.Energy();
  if(eenu<0.100){
    //    cout<<"Not enough energy"<<endl;
    //    return jobc::kFailed;//let's just look at reasonable energy nus
  }
  //  double nuvtxx=nupart.Vx();
  //  double nuvtxy=nupart.Vy();
  double nuvtxz=nupart.Vz();
  int vcont=1;
  
  //  if(nuvtxx>-1.*(geom.DetHalfWidth()-20)&&nuvtxx<1.*(geom.DetHalfWidth()-20)){//20 cm in from xedges    
  //      if(nuvtxy>-1.*(geom.DetHalfHeight()-20)&&nuvtxy<1.*(geom.DetHalfHeight()-20)){//20 cm in from yedges
  //      if(nuvtxz>20&&nuvtxz<(geom.DetLength()-300)){//20 cm from front, and 300 cm from back (muon catcher)
  //      vcont=1;
  //}
  //}
  //}
  if(nuvtxz>(geom.DetLength()-300)){// 300 cm from back (muon catcher)
    vcont=0;
  }   
  if(vcont==0){ 
    //    cout<<"Not in Fid Vol"<<endl;
    //    return jobc::kFailed;
  }

  cout<<"Found an interaction of inu: "<<nupart.GetPdgCode()<<" CC/NC: "<<truelist[tidx]->GetNeutrino()->CCNC()<<endl;
  std::vector<const sim::FLSHitList*> flslist(0);
  try { evt.DetSim().Get(".", flslist); }
  catch(edm::Exception e){
    cerr<<"Couldn't get a flslist"<<endl;
    return jobc::kFailed;
  }


  std::vector<sim::FLSHit> flshits = flslist[tidx]->fHits;
  cout<<"Size of fls "<<flshits.size()<<endl;
  for(unsigned int i=0;i<flshits.size();i++){
    cout<<"Have a FLS hit in PLane: "<<flshits[i].fPlaneId<<" Cell "<<flshits[i].fCellId<<endl;
    cout<<" edep "<<flshits[i].fEdep<<" pdg of parent: "
        <<flshits[i].fPDG<<" track id "<<flshits[i].fTrackId<<endl;
  }


  std::vector<const recobase::Cluster *>xclusterlist(0);
  bool foundx=true;
  try{evt.Reco().Get("./SSClusterX",xclusterlist);}
  catch(edm::Exception e){
    //    cerr<<"Couldn't find an xcluster"<<endl;
    foundx=false;
  }

  std::vector<const recobase::Cluster *>yclusterlist(0);
  bool foundy=true;
  try{evt.Reco().Get("./SSClusterY",yclusterlist);}
  catch(edm::Exception e){
    //    cerr<<"Couldn't find an ycluster"<<endl;
    foundy=false;
  }

  std::vector<const recobase::Cluster *> &clist=xclusterlist;
  /*
  for(int i=0;i<2;i++){
    bool fc=false;
    geo::View_t view;
    if(i==0){
      clist=xclusterlist;
      fc = foundx;
      view=geo::kX;
    }
    else{
      clist=yclusterlist;
      fc = foundy;
      view=geo::kY;
    }
 
    if(fc){
      cout<<"*************In view "<<view<<" found "<<clist.size()<<" clusters "<<endl;
      for(unsigned int i=0;i<clist.size();i++){
        //      cout<<"NEw tv cell"<<endl;
        recobase::Cluster *c = new recobase::Cluster(*clist[i]);
        cout<<"cluster "<<i<<" has "<<c->NCell(view)<<" cells"<<endl;
        for(int j=0;j<c->NCell(view);j++){
          cout<<"Found cell number "<<j<<endl;
          c->Cell(view,j)->Print("");     
        }
      }
    }    
  }
*/

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

---