CosmicTrends_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       CosmicTrends
// Module Type: analyzer
// File:        CosmicTrends_module.cc
//
// Generated at Tue Aug 21 13:38:36 2012 by Brian_Rebel using artmod
// from art v1_00_11.
////////////////////////////////////////////////////////////////////////
#ifndef CosmicTrends_h
#define CosmicTrends_h

#include <map>

#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Core/EDAnalyzer.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "Geometry/Geometry.h"
#include "RecoBase/Track.h"
#include "MCCheater/BackTracker.h"
#include "RunHistory/service/RunHistoryService.h"

#include "TStopwatch.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"

namespace calib {
  class CosmicTrends;

  enum _quality{
    kAll = 0,
    kTriCell,
    kMaxQuality
  };

  enum _oneDHistNames{
    kFractionUsableX = 0,
    kFractionUsableY,
    kRecoCosTheta,
    kRecoAz,
    kPlanesCrossed,
    kDCosZ,
    kDCosX,
    kDCosY,
    kNumTracks
  };

  enum _twoDHistNames{
    kStopperCountX = 0,
    kTotalCountX,
    kStopperCountY,
    kTotalCountY,
    kCosThetaVsAz,
    kCosThetaVsAzTru,
    kCosThetaVsPlanes,
    kAzVsPlanes,
    kDeltaAzVsAz,
    kDeltaZenVsZen,
    kTotalHitCountX,
    kTotalHitCountY
  };
  
}

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

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

  virtual void beginRun(art::Run const&);
  virtual void reconfigure(fhicl::ParameterSet const & p);

private:

  void FillHistograms(rb::Track                                           const& track,
                      unsigned int                                        const& qualityLevel,
                      std::map<unsigned short, std::set<unsigned short> > const& planesAndCells,
                      bool                                                       mc);

  // useable tracks are defined as those that pass through a cell and its
  // neighbors to either side in the same view
  // useable tracks should also span a total number of planes in the detector

  double       fTriCellFrac;     ///< minimum fraction of tri-cells in a track
  unsigned int fMinPlanes;       ///< minimum number of planes a track must cross 
                                 ///< before being deemed useable            
  double       fContainmentCut;  ///< #cm in from the detector boundary to be called a stopper
  std::string  fTrackLabel;      ///< label of module creating the tracks we use               

  std::map<unsigned int, std::vector<TH1D*> > fOneDHists; ///< all 1D hists to make
  std::map<unsigned int, std::vector<TH2D*> > fTwoDHists; ///< all 2D hists to make

  art::ServiceHandle<geo::Geometry> fGeo; 

};

//---------------------------------------------------------------------
calib::CosmicTrends::CosmicTrends(fhicl::ParameterSet const & p)
  : EDAnalyzer(p)
{
  fOneDHists.clear();
  fTwoDHists.clear();

  this->reconfigure(p);
}

//---------------------------------------------------------------------
calib::CosmicTrends::~CosmicTrends()
{
}

//---------------------------------------------------------------------
void calib::CosmicTrends::analyze(art::Event const & e)
{
 
  // get the tracks for this event
  art::Handle< std::vector<rb::Track> > th;
  e.getByLabel(fTrackLabel, th);

  int numTriCellTracks = 0;

  // loop over all the tracks, use only those that cross the minimum
  // number of planes
  for( size_t t = 0; t < th->size(); ++t){
    
    rb::Track const& track = th->at(t);

    // now loop over the hits in the track and make a 
    // map of plane and cells that are hit
    std::map<unsigned short, std::set<unsigned short> > planesAndCells;

    // loop over the planes in the map and check out the cells 
    // that were hit
    double useableCells[2] = {0.};
    double totalCells[2]   = {0.};

    for(size_t h = 0; h < track.NCell(); ++h)
      planesAndCells[track.Cell(h)->Plane()].insert(track.Cell(h)->Cell());

    useableCells[0] = 0.;
    useableCells[1] = 0.;
    totalCells[0] = 0.;
    totalCells[1] = 0.;
  
    for(auto const& pacitr : planesAndCells){
      
      unsigned short plane = pacitr.first;
      // get the set for the current plane
      // sets are presorted, go c++
      std::set<unsigned short> cellSet = pacitr.second;
      
      // put the iterator over the set to the second entry
      //auto citr = cellSet.cbegin();
      //totalCells[fGeo->Plane(pacitr.first)->View()] += 1.;
      //++citr;
      //unsigned short cell = citr;
      // loop over the cells and look to see if there is 
      // a cell in this plane that has both its neighbors also
      // hit on this track
      //while( citr != cellSet.end() ){
      for( auto citr : cellSet ){
        if( cellSet.find(citr - 1) != cellSet.end() &&
           cellSet.find(citr + 1) != cellSet.end() ){
          useableCells[fGeo->Plane(plane)->View()] += 1.;
          
          if(fGeo->Plane(plane)->View() == geo::kX){
            fTwoDHists[kAll][kTotalCountX]->Fill(plane, citr);
          }
          else if(fGeo->Plane(plane)->View() == geo::kY){
            fTwoDHists[kAll][kTotalCountY]->Fill(plane, citr);
          }
        }
        
        if(fGeo->Plane(plane)->View() == geo::kX){
          fTwoDHists[kAll][kTotalHitCountX]->Fill(plane, citr);
        }
        else if(fGeo->Plane(plane)->View() == geo::kY){
          fTwoDHists[kAll][kTotalHitCountY]->Fill(plane, citr);
        }
        
        totalCells[fGeo->Plane(plane)->View()] += 1.;   
        ++citr;
      }// end iterator over cells for this plane

    } // end iterator over plane/cell map
    
    double usefulFrac = 0.;
    if(totalCells[0] + totalCells[1] > 0)
      usefulFrac = (useableCells[0] + useableCells[1])/(totalCells[0] + totalCells[1]);
    
    // loop over the planesAndCells map again and fill histograms
    // for cells that are passed tracks
    unsigned int quality = (usefulFrac >= fTriCellFrac) ? kTriCell : kAll;

    if(totalCells[0] > 0){
      fOneDHists[kAll][kFractionUsableX]->Fill(useableCells[0]/totalCells[0]);
      if(quality == kTriCell)
        fOneDHists[quality][kFractionUsableX]->Fill(useableCells[0]/totalCells[0]);
    }
    if(totalCells[1] > 0){
      fOneDHists[kAll][kFractionUsableY]->Fill(useableCells[1]/totalCells[1]);
      if(quality == kTriCell)
        fOneDHists[quality][kFractionUsableY]->Fill(useableCells[1]/totalCells[1]);
    }
    
    this->FillHistograms(track, kAll, planesAndCells, !e.isRealData());
    if(quality == kTriCell){
      ++numTriCellTracks;
      this->FillHistograms(track, quality, planesAndCells, !e.isRealData());
    }    
    
  }// end loop over tracks in this event

  fOneDHists[kTriCell][kNumTracks]->Fill(numTriCellTracks);
  fOneDHists[kAll][kNumTracks]->Fill(th->size());

  return;
}

//---------------------------------------------------------------------
void calib::CosmicTrends::beginRun(art::Run const&)
{
  if(fOneDHists.size() > 0) return;

  art::ServiceHandle<art::TFileService> tfs;

  // use the geometry to get the number of cells for planes in each view
  // also get the total number of planes

  unsigned int planeInView = *(fGeo->GetPlanesByView(geo::kX).begin());  
  double xCells = fGeo->Plane(planeInView)->Ncells();
  
  planeInView = *(fGeo->GetPlanesByView(geo::kY).begin());  
  double yCells = fGeo->Plane(planeInView)->Ncells();
  
  double totalPlanes = fGeo->NPlanes();
  double maxCells = (yCells > xCells) ? yCells : xCells;

  for(int q = kAll; q < kMaxQuality; ++q){

    std::string qualString("AllTracks");
    if(q == kTriCell) qualString = "TriCell";

    // # of useable stopping muons passing through each cell
    fTwoDHists[q].push_back(tfs->make<TH2D>(("stopperCountX"+qualString).c_str(), 
                                            ";Plane;Cell",
                                            (int)totalPlanes, 0., totalPlanes,
                                            (int)maxCells,    0., maxCells));

    // # of useable muons passing through each cell         
    fTwoDHists[q].push_back(tfs->make<TH2D>(("totalCountX"+qualString).c_str(), 
                                            ";Plane;Cell",
                                            (int)totalPlanes, 0., totalPlanes,
                                            (int)maxCells,    0., maxCells));
  
    // # of useable stopping muons passing through each cell
    fTwoDHists[q].push_back(tfs->make<TH2D>(("stopperCountY"+qualString).c_str(), 
                                            ";Plane;Cell",
                                            (int)totalPlanes, 0., totalPlanes,
                                            (int)maxCells,    0., maxCells));
    // # of useable muons passing through each cell         
    fTwoDHists[q].push_back(tfs->make<TH2D>(("totalCountY"+qualString).c_str(), 
                                            ";Plane;Cell",
                                            (int)totalPlanes, 0., totalPlanes,
                                            (int)maxCells,    0., maxCells));
  
    // reconstructed cos(theta) vs reconstructed azimuth
    fTwoDHists[q].push_back(tfs->make<TH2D>(("cosThetaVsAz"+qualString).c_str(), 
                                            ";cos(#theta);#phi (degrees)",
                                            100,              0., 1.,
                                            36,               0., 360.));
    
    // reconstructed cos(theta) vs reconstructed azimuth
    fTwoDHists[q].push_back(tfs->make<TH2D>(("cosThetaVsAzTru"+qualString).c_str(), 
                                            ";True cos(#theta);True #phi (degrees)",
                                            100,              0., 1.,
                                            36,               0., 360.));
    
    // reconstructed cos(theta) vs planes crossed
    fTwoDHists[q].push_back(tfs->make<TH2D>(("cosThetaVsPlanes"+qualString).c_str(), 
                                            ";cos(#theta);Planes Crossed",
                                            100,              0., 1.,
                                            (int)totalPlanes, 0., totalPlanes));
                           
    // reconstructed azimuth vs planes crossed
    fTwoDHists[q].push_back(tfs->make<TH2D>(("azimuthVsPlanes"+qualString).c_str(), 
                                            ";#phi (degrees);Planes Crossed",
                                            36,               0., 360.,
                                            (int)totalPlanes, 0., totalPlanes));
    
    // reco - true vs reconstructed azimuth
    fTwoDHists[q].push_back(tfs->make<TH2D>(("deltaAzVsAz"+qualString).c_str(), 
                                            ";#phi (degrees);#Delta#phi (degrees)",
                                            36,               0., 360.,
                                            360,           -180., 180.));
    
    // reco - true vs reconstructed cos(theta)
    fTwoDHists[q].push_back(tfs->make<TH2D>(("deltaZenVsZen"+qualString).c_str(), 
                                            ";cos(#theta);#Deltacos(#theta)",
                                            100,              0., 1.,
                                            200,             -1., 1.));
    
    // fraction of tri-cells in a muon track
    fOneDHists[q].push_back(tfs->make<TH1D>(("fractionUsableX"+qualString).c_str(), 
                                            ";Usable Fraction;Number of Muons",
                                            100, 0., 1.));

    // fraction of tri-cells in a muon track
    fOneDHists[q].push_back(tfs->make<TH1D>(("fractionUsableY"+qualString).c_str(), 
                                            ";Usable Fraction;Number of Muons",
                                            100, 0., 1.));
    
    // reconstructed cos(theta)
    fOneDHists[q].push_back(tfs->make<TH1D>(("recoCosTheta"+qualString).c_str(),   
                                            ";cos(#theta);Number of Muons",
                                            100, 0., 1.   ));
    
    // reconstructed azimuth
    fOneDHists[q].push_back(tfs->make<TH1D>(("recoAzimuth"+qualString).c_str(),   
                                            ";#phi (degrees);Number of Muons",
                                            36,  0., 360. ));
    
    // reconstructed number of planes crossed by track
    fOneDHists[q].push_back(tfs->make<TH1D>(("planesCrossed"+qualString).c_str(),   
                                            ";Planes Crossed;Number of Muons",
                                            (int)totalPlanes*0.2, 0., (int)totalPlanes));
    
    // direction cosine in z 
    fOneDHists[q].push_back(tfs->make<TH1D>(("dcosZ"+qualString).c_str(),   
                                            ";dz/ds;Number of Tracks",
                                            40, -1., 1.));
    
    // direction cosine in x 
    fOneDHists[q].push_back(tfs->make<TH1D>(("dcosX"+qualString).c_str(),   
                                            ";dx/ds;Number of Tracks",
                                            40, -1., 1.));
    
    // direction cosine in y
    fOneDHists[q].push_back(tfs->make<TH1D>(("dcosY"+qualString).c_str(),   
                                            ";dy/ds;Number of TriCell Tracks",
                                            40, -1., 1.));
    
    // tracks in event
    fOneDHists[q].push_back(tfs->make<TH1D>(("numTracks"+qualString).c_str(),   
                                            ";Number of TriCell Tracks;Events",
                                            10, 0., 10.));

    // # of muons passing through each cell         
    fTwoDHists[q].push_back(tfs->make<TH2D>(("totalHitCountX"+qualString).c_str(), 
                                            ";Plane;Cell",
                                            (int)totalPlanes, 0., totalPlanes,
                                            (int)maxCells,    0., maxCells));
    // # of muons passing through each cell         
    fTwoDHists[q].push_back(tfs->make<TH2D>(("totalHitCountY"+qualString).c_str(), 
                                            ";Plane;Cell",
                                            (int)totalPlanes, 0., totalPlanes,
                                            (int)maxCells,    0., maxCells));
    

    
  }

  return;
}

//---------------------------------------------------------------------
void calib::CosmicTrends::reconfigure(fhicl::ParameterSet const & p)
{
  fTriCellFrac    = p.get< double       >("TriCellFraction",       0.6);
  fMinPlanes      = p.get< unsigned int >("MinimumPlanesPerTrack",   8);
  fTrackLabel     = p.get< std::string  >("TrackLabel"                );
  fContainmentCut = p.get< double       >("ContainmentCut",        20.);

  return;
}

//---------------------------------------------------------------------
void calib::CosmicTrends::FillHistograms(rb::Track                                           const& track,
                                         unsigned int                                        const& qualityLevel,
                                         std::map<unsigned short, std::set<unsigned short> > const& planesAndCells,
                                         bool                                                       mc)
{

  fOneDHists[qualityLevel][kPlanesCrossed]->Fill(track.ExtentPlane());
  
  if(track.ExtentPlane() < fMinPlanes) return;

  TVector3 startDir = track.Dir();
  double   cosTheta = -startDir.y()/startDir.Mag();
  double   azimuth  = std::atan2(startDir.x()/startDir.Mag(), -startDir.z()/startDir.Mag()) * 180./TMath::Pi(); 
  if(azimuth < 0) azimuth += 360.;

  fTwoDHists[qualityLevel][kCosThetaVsAz]    ->Fill(cosTheta, azimuth);
  fTwoDHists[qualityLevel][kCosThetaVsPlanes]->Fill(cosTheta, track.ExtentPlane());
  fTwoDHists[qualityLevel][kAzVsPlanes]      ->Fill(azimuth,  track.ExtentPlane());

  fOneDHists[qualityLevel][kRecoCosTheta]->Fill(cosTheta);
  fOneDHists[qualityLevel][kRecoAz]      ->Fill(azimuth);

  // add code to determine reconstruction qualityLevel of the 
  // tracking in terms of azimuth and cos(theta) when 
  // back tracker is able to do that.
  if(mc){
    art::ServiceHandle<cheat::BackTracker> bt;
    
    // get the MCParticle that contributes the most information to the
    // current track
    const sim::Particle* truth = bt->HitsToParticle(track.AllCells()).at(0);
    if(!truth) 
      LOG_WARNING("CosmicTrends") << "no simulated particle matches current track, skip it";
    else{
      float trueCosTheta = -truth->Py()/truth->P();
      float trueAzimuth  = std::atan2(truth->Px()/truth->P(), -truth->Pz()/truth->P()) * 180./TMath::Pi();
      if(trueAzimuth < 0) trueAzimuth += 360.;
      fTwoDHists[qualityLevel][kCosThetaVsAzTru]->Fill(trueCosTheta, trueAzimuth);
      fTwoDHists[qualityLevel][kDeltaAzVsAz]    ->Fill(azimuth,      azimuth-trueAzimuth);
      fTwoDHists[qualityLevel][kDeltaZenVsZen]  ->Fill(cosTheta,     cosTheta-trueCosTheta);
    }
  }

  // check if this track stops in the detector
  TVector3 startPoint(track.Start());
  TVector3 endPoint(track.Stop());
  bool stopper = false;
  if(std::abs(endPoint.X()) < fGeo->DetHalfWidth()  - fContainmentCut &&
     std::abs(endPoint.Y()) < fGeo->DetHalfHeight() - fContainmentCut &&
     endPoint.Z()  < fGeo->DetLength()     - fContainmentCut &&
     endPoint.Z()  > fContainmentCut                        )
    stopper = true;

  for(auto const& pacitr : planesAndCells){
    
    unsigned short plane = pacitr.first;
    // get the set for the current plane
    // sets are presorted, go c++
    std::set<unsigned short> cellSet = pacitr.second;
    
    // loop over the cells and look to see if there is 
    // a cell in this plane that has both its neighbors also
    // hit on this track
    for( auto citr : cellSet ){
    
//       if(fGeo->Plane(plane)->View() == geo::kX)
//      fTwoDHists[qualityLevel][kTotalCountX]->Fill(plane, citr);
//       else if(fGeo->Plane(plane)->View() == geo::kY)
//      fTwoDHists[qualityLevel][kTotalCountY]->Fill(plane, citr);
    
      if(stopper){
        if(fGeo->Plane(plane)->View() == geo::kX)
          fTwoDHists[qualityLevel][kStopperCountX]->Fill(plane, citr);
        else if(fGeo->Plane(plane)->View() == geo::kY)
          fTwoDHists[qualityLevel][kStopperCountY]->Fill(plane, citr);
      }
      
      ++citr;
    }// end iterator over cells for this plane
  } // end iterator over plane/cell map
  
  fOneDHists[qualityLevel][kDCosZ]->Fill(startDir.Z()/startDir.Mag());
  fOneDHists[qualityLevel][kDCosX]->Fill(startDir.X()/startDir.Mag());
  fOneDHists[qualityLevel][kDCosY]->Fill(startDir.Y()/startDir.Mag());

  return;
}

DEFINE_ART_MODULE(calib::CosmicTrends)

#endif /* CosmicTrends_h */