ClusterAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       ClusterAna
// Module Type: analyzer
// File:        ClusterAna_module.cc
//
// Generated at Tue Apr  2 09:40:02 2013 by Jan Zirnstein using artmod
// from art v1_02_06.
////////////////////////////////////////////////////////////////////////

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/FindOneP.h"

// NOvADDT includes
#include "DDTBaseDataProducts/Cluster.h"
#include "DDTBaseDataProducts/BaseProducts.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/CompareDAQHit.h"
#include "DDTUtilities/func/MathUtil.h"

// ROOT includes
#include "TH1F.h"

#include <algorithm>
#include <cmath>
#include <vector>

namespace novaddt {
  class ClusterAna;
}

class novaddt::ClusterAna : public art::EDAnalyzer {
public:
  explicit ClusterAna(fhicl::ParameterSet const & p);
  virtual ~ClusterAna();

  void analyze(art::Event const & e) override;

  void beginJob() override;
  void endJob() override;

private:

  TH1F* fNClustersHist;
  TH1F* fNHitsPerClusterHist;
  TH1F* fExtentPlanesHist;
  TH1F* fClusterEnergyHist;
  TH1F* fWidthXHist;
  TH1F* fWidthYHist;
  TH1F* fClusterEnDenHist;
  TH1F* fXViewTrackGood;
  TH1F* fYViewTrackGood;
  TH1F* fCombViewTrackGood;

  std::string fClusterModule;

};


novaddt::ClusterAna::ClusterAna(fhicl::ParameterSet const & p) :
  EDAnalyzer(p),
  fClusterModule(p.get<std::string>("ClusterModule"))
{
}

novaddt::ClusterAna::~ClusterAna()
{
  // Clean up dynamic memory and other resources here.
}

void novaddt::ClusterAna::analyze(art::Event const & e)
{
  art::Handle< std::vector<novaddt::Cluster> > clusters;
  e.getByLabel(fClusterModule, clusters);

  // put in a mf::LogWarning why this statement returns if true
  if(clusters->empty()){
    mf::LogWarning ("No Clusters")<<"No Clusters in the input file";
    return;
  }

  // Fill Number of Clusters in Event Histogram
  fNClustersHist->Fill(clusters->size());

  // Grab the Hitlists associated with the clusters
  art::FindOneP<novaddt::HitList> hls(clusters, e, fClusterModule);

  //std::cout<< "Sanity check, there are: "<<clusters->size()<<" clusters and " <<hls.size()<<" hitlists."<<std::endl;

  unsigned int i(0);
  for(auto const& clust: *clusters){
    fExtentPlanesHist->Fill(clust.Extent());
    fWidthXHist->Fill(clust.WidthX());
    fWidthYHist->Fill(clust.WidthY());
    fClusterEnergyHist->Fill(clust.TotalADC());
    fClusterEnDenHist->Fill(clust.TotalADC()/(1.0*clust.Volume()));

    art::Ptr<novaddt::HitList> hl_ptr(hls.at(i));
    ++i;

    novaddt::HitList hlist(*hl_ptr);

    fNHitsPerClusterHist->Fill(hl_ptr->size());
    //std::cout<<"Hits in cluster: "<<hlist.size()<<std::endl;

    // Split the hits in x and y to do 2D fitting
    auto bound = std::partition(hlist.begin(),hlist.end(),CompareDAQHit<View>());
    //std::cout<<"Boundary is at: "<<bound-hlist.begin()<<std::endl;

    unsigned int xsize = bound-hlist.begin();
    unsigned int ysize = hlist.end()-bound;

    //std::cout<<"X Size: "<<xsize<<", Y size: "<<ysize<<std::endl;

    // Smarter way of vector managment 

    std::vector<double> XXpoints;
    std::vector<double> XZpoints;
    std::vector<double> YYpoints;
    std::vector<double> YZpoints;

    XXpoints.reserve(xsize);
    XZpoints.reserve(xsize);
    YYpoints.reserve(ysize);
    YZpoints.reserve(ysize);

    {
      auto it = hlist.begin();
      for(it=hlist.begin();it!=bound;++it){
        XXpoints.emplace_back((*it).Cell().val);
        XZpoints.emplace_back((*it).Plane().val);
      }

      for(it=bound;it!=hlist.end();++it){
        YYpoints.emplace_back((*it).Cell().val);
        YZpoints.emplace_back((*it).Plane().val);
      }
    }

    // Fit to a straight line
    double xslope;
    double xint;
    unsigned int xdof = xsize-1;
    double xsigm=0.0;
    double pred;

    if(xsize > 1){
      util::LinFitUnweighted(XXpoints, XZpoints, xslope, xint);
      //std::cout<<"XView fit done: y = "<<xslope<<"*x + "<<xint<<std::endl;
  
      // Evaluate goodness of fit
      for(unsigned int xi=0;xi<xsize;++xi){
        pred = xslope*XXpoints[xi] + xint;
        xsigm += sqrt(util::sqr(XZpoints[xi]-pred));
      }

      fXViewTrackGood->Fill(xsigm/xdof);
    }

    else{
      fXViewTrackGood->Fill(0.0);
      xdof = 1;
    }
    //std::cout<<"Goodness of XView fit: "<<xsigm/xdof<<std::endl;

    double yslope;
    double yint;
    unsigned int ydof = ysize-1;
    double ysigm=0.0;

    if(ysize > 1){
      util::LinFitUnweighted(YYpoints, YZpoints, yslope, yint);
      //std::cout<<"YView fit done: y = "<<yslope<<"*x + "<<yint<<std::endl;

      for(unsigned int yi=0;yi<ysize;++yi){
        pred = yslope*YYpoints[yi] + yint;
        ysigm += sqrt(util::sqr(YZpoints[yi]-pred));
      }

      fYViewTrackGood->Fill(ysigm/ydof);
    }

    else{
      fYViewTrackGood->Fill(0.0);
      ydof = 1;
    }
    //std::cout<<"Goodness of YView fit: "<<ysigm/ydof<<std::endl;

    fCombViewTrackGood->Fill(sqrt(util::sqr(xsigm/xdof)+util::sqr(ysigm/ydof)));

    //std::cout<<"Analysis of Cluster finished\n\n";

  } // end of cluster

} // end of analyze function

void novaddt::ClusterAna::beginJob()
{
  art::ServiceHandle<art::TFileService> f;

  fNClustersHist       = f->make<TH1F>("fNClusters",";NClusters;Events",100,0.,100.);
  fNHitsPerClusterHist = f->make<TH1F>("fNHitsPerCluster",";NHits;Clusters",100,0.,100.);
  fExtentPlanesHist    = f->make<TH1F>("fExtentPlanes",";NPlanes;Clusters",896,0.,896.);
  fWidthXHist          = f->make<TH1F>("fWidthX",";NCells;Clusters",384,0.,384.);
  fWidthYHist          = f->make<TH1F>("fWidthY",";NCells;Clusters",384,0.,384.);
  fClusterEnergyHist   = f->make<TH1F>("fClusterEnergy",";TotADC;Clusters",1000,0.,100000.);
  fClusterEnDenHist    = f->make<TH1F>("fClusterEnergyDensity",";ADC per plane*cell^2;Clusters",10000,0.,0.5);
  fXViewTrackGood      = f->make<TH1F>("fXViewTrackGood",";Goodness of Lin Fit;Tracks",1000,0.,100.);
  fYViewTrackGood      = f->make<TH1F>("fYViewTrackGood",";Goodness of Lin Fit;Tracks",1000,0.,100.);
  fCombViewTrackGood   = f->make<TH1F>("fCombViewTrackGood",";Goodness of Lin Fit;Tracks",1000,0.,100.);

}

void novaddt::ClusterAna::endJob()
{
  // Implementation of optional member function here.
}

DEFINE_ART_MODULE(novaddt::ClusterAna)