MuondEdx_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       MuondEdx
// Module Type: filter
// File:        MuondEdx_module.cc
//
// Generated at Tue May  7 16:23:55 2013 by Abigail Waldron using artmod
// from cetpkgsupport v1_02_00.
////////////////////////////////////////////////////////////////////////


// framework includes....
#include "art/Framework/Core/EDFilter.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/FindManyP.h" 

// novasoft stuff
#include "CalibrationDataProducts/PCHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/CellHit.h"
#include "Calibrator/Calibrator.h"
#include "GeometryObjects/GeometryBase.h"

// C++
#include <memory>
#include <string>
#include <iostream>
#include <vector>

// ROOT
#include "TH2D.h"
#include "TFile.h"
#include "TVector3.h"

namespace caldp {
  class PCHit;
}


namespace calib {
  class MuondEdx;
  
  class MuondEdx : public art::EDFilter {
  public:
    explicit MuondEdx(fhicl::ParameterSet const & p);
    virtual ~MuondEdx();
    
    bool filter(art::Event & e) override;
    void reconfigure(fhicl::ParameterSet const & p);
    void beginJob() override;
    void endJob() override;
    bool beginRun(art::Run & r) override;
    bool endRun(art::Run & r) override;
    
  private:
    double getPECorr(art::Event & e,caldp::PCHit const &pchit);

    TH2D* fdEdx;
    TH2D* fdEdx_minus;
    TH2D* fdEdx_plus;
    TH2D* fdEdx_true; 

    std::string fPCHitLabel;
    std::string fQualXYName;
    std::string fStopperLabel;

    double fSysShift;

  };
  
  
  MuondEdx::MuondEdx(fhicl::ParameterSet const & p)
  {
    this->reconfigure(p);
    produces< TH2D, art::InRun >("dEdx");
    produces< TH2D, art::InRun >("dEdxMinus");
    produces< TH2D, art::InRun >("dEdxPlus");
    produces< TH2D, art::InRun >("dEdxTrue");
  }
  
  MuondEdx::~MuondEdx()
  {
    // Clean up dynamic memory and other resources here.
  }
  
  bool MuondEdx::filter(art::Event & e)
  {
    /// Get the reconstructed tracks                                           
    art::Handle<std::vector<rb::Track> > tracks;
    e.getByLabel(fStopperLabel, tracks);
    if( tracks->size() == 0 ) return true;

    art::FindManyP<caldp::PCHit> PCHitGetter(tracks, e, art::InputTag(fPCHitLabel,fQualXYName)); 
    if( !PCHitGetter.isValid() ) return false;    

    // loop over stopping tracks in event
    for( unsigned int i=0; i<tracks->size(); ++i ){
      
      // needed for track window, but good generally to reject protons
      if( (*tracks)[i].TotalLength() < 200. ) continue;  

      TVector3 stopXYZ;
      art::ServiceHandle<geo::Geometry> fGeo;
      geo::CellUniqueId cid;
      int stopPlaneID = 10000;
      int firstMuCPlane = -1;
      bool isNearDet = false;

      //continue if track ends in the muon catcher
      if(fGeo->DetId() == novadaq::cnv::kNEARDET){
        isNearDet = true;
        stopXYZ = (*tracks)[i].Start();
        cid = fGeo->CellId(stopXYZ[0], stopXYZ[1], stopXYZ[2], 1.5, 1.5, 6., 1.0);
        if(cid) stopPlaneID = fGeo->getPlaneID(cid);
        firstMuCPlane = fGeo->FirstPlaneInMuonCatcher();
        if(stopPlaneID >= firstMuCPlane) continue;
      }

      // loop over PCHits associated to stopping track
      std::vector<art::Ptr<caldp::PCHit> > hits = PCHitGetter.at(i);
      if( hits.size() == 0 ) continue;      
      for( unsigned int j=0; j<hits.size(); ++j ){
        
        //cut if hit is inside muon catcher
        if(isNearDet)
          if(hits[j]->Plane() >= firstMuCPlane) continue;

        double x = (*tracks)[i].TotalLength() - hits[j]->FlightLen();
        double path = hits[j]->Path();
        double trueE = hits[j]->TrueMeV();
        double PECorr = this->getPECorr(e,*hits[j]);
        fdEdx->Fill(x,PECorr/path);
        fdEdx_minus->Fill(x-fSysShift,PECorr/path);
        fdEdx_plus->Fill(x+fSysShift,PECorr/path);
        fdEdx_true->Fill(x,trueE/path);
      }
    }


    return true;


  }
  double MuondEdx::getPECorr(art::Event & e, caldp::PCHit const &pchit){
   
    //Calibrator necessary to get Atten information                             
    art::ServiceHandle<calib::Calibrator> calibrator;

    //Make mock cellhit
    rb::CellHit cellhit;
    cellhit.SetPlane(pchit.Plane());
    cellhit.SetCell(pchit.Cell());
    cellhit.SetView(pchit.View());
    cellhit.SetTNS(pchit.TNS(),pchit.GoodTime());
    if(!e.isRealData())cellhit.SetMC(); // assume only filled for MC...
    cellhit.SetPE(pchit.PE());

    double pecorr = calibrator->GetPECorr(cellhit, pchit.W());
    return pecorr;
  }

  void MuondEdx::reconfigure(fhicl::ParameterSet const & pset)
  {
    fPCHitLabel = pset.get< std::string >("PCHitLabel"); //Label of PCHitList module            
    fQualXYName = pset.get< std::string >("QualXYName"); //Instance label, "XY" quality hits    
    fStopperLabel = pset.get<std::string>("StopperLabel");
    fSysShift = pset.get<double>("SysShift"); // track end point systematic shift
  }
  
  void MuondEdx::beginJob()
  {
    // histograms for output
    art::ServiceHandle<art::TFileService> tfs;
    fdEdx = tfs->make<TH2D>("dEdx",";Distance to track end (cm);PECorr / cm;",
                            50,0.,500.,200,0.,200.);
    fdEdx_minus = tfs->make<TH2D>("dEdx_minus",";Distance to track end (cm);PECorr / cm;",
                            50,0.,500.,200,0.,200.);
    fdEdx_plus = tfs->make<TH2D>("dEdx_plus",";Distance to track end (cm);PECorr / cm;",
                            50,0.,500.,200,0.,200.);
    fdEdx_true = tfs->make<TH2D>("dEdx_true",";Distance to track end (cm);True MeV / cm;",
                                 50,0.,500.,200,0.,5.);
  }
  
  bool MuondEdx::beginRun(art::Run &r){
    return true;
  }

  bool MuondEdx::endRun(art::Run &r){

    art::ServiceHandle<art::TFileService> tfs;
    std::unique_ptr<TH2D> pdEdx(tfs->make<TH2D>(*fdEdx));
    std::unique_ptr<TH2D> pdEdx_minus(tfs->make<TH2D>(*fdEdx_minus));
    std::unique_ptr<TH2D> pdEdx_plus(tfs->make<TH2D>(*fdEdx_plus));
    std::unique_ptr<TH2D> pdEdx_true(tfs->make<TH2D>(*fdEdx_true));

    r.put(std::move(pdEdx),"dEdx");
    r.put(std::move(pdEdx_minus),"dEdxMinus");
    r.put(std::move(pdEdx_plus),"dEdxPlus");
    r.put(std::move(pdEdx_true),"dEdxTrue");
    return true;
  }

  void MuondEdx::endJob()
  {

  }
}// end namespace
DEFINE_ART_MODULE(calib::MuondEdx)