MuondEdx_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MuondEdx
3 // Module Type: filter
4 // File: MuondEdx_module.cc
5 //
6 // Generated at Tue May 7 16:23:55 2013 by Abigail Waldron using artmod
7 // from cetpkgsupport v1_02_00.
8 ////////////////////////////////////////////////////////////////////////
9 
10 
11 // framework includes....
19 #include "fhiclcpp/ParameterSet.h"
23 
24 // novasoft stuff
26 #include "RecoBase/Track.h"
27 #include "RecoBase/CellHit.h"
28 #include "Calibrator/Calibrator.h"
30 
31 // C++
32 #include <memory>
33 #include <string>
34 #include <iostream>
35 #include <vector>
36 
37 // ROOT
38 #include "TH2D.h"
39 #include "TFile.h"
40 #include "TVector3.h"
41 
42 namespace caldp {
43  class PCHit;
44 }
45 
46 
47 namespace calib {
48  class MuondEdx;
49 
50  class MuondEdx : public art::EDFilter {
51  public:
52  explicit MuondEdx(fhicl::ParameterSet const & p);
53  virtual ~MuondEdx();
54 
55  bool filter(art::Event & e) override;
56  void reconfigure(fhicl::ParameterSet const & p);
57  void beginJob() override;
58  void endJob() override;
59  bool beginRun(art::Run & r) override;
60  bool endRun(art::Run & r) override;
61 
62  private:
63  double getPECorr(art::Event & e,caldp::PCHit const &pchit);
64 
65  TH2D* fdEdx;
66  TH2D* fdEdx_minus;
67  TH2D* fdEdx_plus;
68  TH2D* fdEdx_true;
69 
73 
74  double fSysShift;
75 
76  };
77 
78 
79  MuondEdx::MuondEdx(fhicl::ParameterSet const & p)
80  {
81  this->reconfigure(p);
82  produces< TH2D, art::InRun >("dEdx");
83  produces< TH2D, art::InRun >("dEdxMinus");
84  produces< TH2D, art::InRun >("dEdxPlus");
85  produces< TH2D, art::InRun >("dEdxTrue");
86  }
87 
88  MuondEdx::~MuondEdx()
89  {
90  // Clean up dynamic memory and other resources here.
91  }
92 
93  bool MuondEdx::filter(art::Event & e)
94  {
95  /// Get the reconstructed tracks
97  e.getByLabel(fStopperLabel, tracks);
98  if( tracks->size() == 0 ) return true;
99 
100  art::FindManyP<caldp::PCHit> PCHitGetter(tracks, e, art::InputTag(fPCHitLabel,fQualXYName));
101  if( !PCHitGetter.isValid() ) return false;
102 
103  // loop over stopping tracks in event
104  for( unsigned int i=0; i<tracks->size(); ++i ){
105 
106  // needed for track window, but good generally to reject protons
107  if( (*tracks)[i].TotalLength() < 200. ) continue;
108 
109  TVector3 stopXYZ;
112  int stopPlaneID = 10000;
113  int firstMuCPlane = -1;
114  bool isNearDet = false;
115 
116  //continue if track ends in the muon catcher
117  if(fGeo->DetId() == novadaq::cnv::kNEARDET){
118  isNearDet = true;
119  stopXYZ = (*tracks)[i].Start();
120  cid = fGeo->CellId(stopXYZ[0], stopXYZ[1], stopXYZ[2], 1.5, 1.5, 6., 1.0);
121  if(cid) stopPlaneID = fGeo->getPlaneID(cid);
122  firstMuCPlane = fGeo->FirstPlaneInMuonCatcher();
123  if(stopPlaneID >= firstMuCPlane) continue;
124  }
125 
126  // loop over PCHits associated to stopping track
127  std::vector<art::Ptr<caldp::PCHit> > hits = PCHitGetter.at(i);
128  if( hits.size() == 0 ) continue;
129  for( unsigned int j=0; j<hits.size(); ++j ){
130 
131  //cut if hit is inside muon catcher
132  if(isNearDet)
133  if(hits[j]->Plane() >= firstMuCPlane) continue;
134 
135  double x = (*tracks)[i].TotalLength() - hits[j]->FlightLen();
136  double path = hits[j]->Path();
137  double trueE = hits[j]->TrueMeV();
138  double PECorr = this->getPECorr(e,*hits[j]);
139  fdEdx->Fill(x,PECorr/path);
140  fdEdx_minus->Fill(x-fSysShift,PECorr/path);
141  fdEdx_plus->Fill(x+fSysShift,PECorr/path);
142  fdEdx_true->Fill(x,trueE/path);
143  }
144  }
145 
146 
147  return true;
148 
149 
150  }
151  double MuondEdx::getPECorr(art::Event & e, caldp::PCHit const &pchit){
152 
153  //Calibrator necessary to get Atten information
155 
156  //Make mock cellhit
157  rb::CellHit cellhit;
158  cellhit.SetPlane(pchit.Plane());
159  cellhit.SetCell(pchit.Cell());
160  cellhit.SetView(pchit.View());
161  cellhit.SetTNS(pchit.TNS(),pchit.GoodTime());
162  if(!e.isRealData())cellhit.SetMC(); // assume only filled for MC...
163  cellhit.SetPE(pchit.PE());
164 
165  double pecorr = calibrator->GetPECorr(cellhit, pchit.W());
166  return pecorr;
167  }
168 
169  void MuondEdx::reconfigure(fhicl::ParameterSet const & pset)
170  {
171  fPCHitLabel = pset.get< std::string >("PCHitLabel"); //Label of PCHitList module
172  fQualXYName = pset.get< std::string >("QualXYName"); //Instance label, "XY" quality hits
173  fStopperLabel = pset.get<std::string>("StopperLabel");
174  fSysShift = pset.get<double>("SysShift"); // track end point systematic shift
175  }
176 
178  {
179  // histograms for output
181  fdEdx = tfs->make<TH2D>("dEdx",";Distance to track end (cm);PECorr / cm;",
182  50,0.,500.,200,0.,200.);
183  fdEdx_minus = tfs->make<TH2D>("dEdx_minus",";Distance to track end (cm);PECorr / cm;",
184  50,0.,500.,200,0.,200.);
185  fdEdx_plus = tfs->make<TH2D>("dEdx_plus",";Distance to track end (cm);PECorr / cm;",
186  50,0.,500.,200,0.,200.);
187  fdEdx_true = tfs->make<TH2D>("dEdx_true",";Distance to track end (cm);True MeV / cm;",
188  50,0.,500.,200,0.,5.);
189  }
190 
191  bool MuondEdx::beginRun(art::Run &r){
192  return true;
193  }
194 
195  bool MuondEdx::endRun(art::Run &r){
196 
198  std::unique_ptr<TH2D> pdEdx(tfs->make<TH2D>(*fdEdx));
199  std::unique_ptr<TH2D> pdEdx_minus(tfs->make<TH2D>(*fdEdx_minus));
200  std::unique_ptr<TH2D> pdEdx_plus(tfs->make<TH2D>(*fdEdx_plus));
201  std::unique_ptr<TH2D> pdEdx_true(tfs->make<TH2D>(*fdEdx_true));
202 
203  r.put(std::move(pdEdx),"dEdx");
204  r.put(std::move(pdEdx_minus),"dEdxMinus");
205  r.put(std::move(pdEdx_plus),"dEdxPlus");
206  r.put(std::move(pdEdx_true),"dEdxTrue");
207  return true;
208  }
209 
210  void MuondEdx::endJob()
211  {
212 
213  }
214 }// end namespace
float TNS() const
Return uncorrected hit time.
Definition: PCHit.h:54
void SetTNS(float tns, bool good)
Definition: CellHit.h:56
const char * p
Definition: xmltok.h:285
bool GoodTime() const
Return quality of timing fit for cell.
Definition: PCHit.h:56
std::string fStopperLabel
int Cell() const
Return cell value.
Definition: PCHit.h:26
std::string fQualXYName
bool isRealData() const
Definition: Event.h:83
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
"Pre-calibration hit". Common input to calibration procedures
Definition: PCHit.h:16
Definition: Run.h:31
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
void SetView(geo::View_t view)
Definition: CellHit.h:54
void SetPlane(unsigned short plane)
Definition: CellHit.h:53
float W() const
Return W value.
Definition: PCHit.h:44
void SetCell(unsigned short cell)
Definition: CellHit.h:52
CDPStorage service.
void hits()
Definition: readHits.C:15
void SetPE(float pe)
Definition: CellHit.h:55
float PE() const
Return PE value.
Definition: PCHit.h:38
void beginJob()
Encapsulate the geometry of one entire detector (near, far, ndos)
T get(std::string const &key) const
Definition: ParameterSet.h:231
int getPlaneID(const CellUniqueId &id) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Near Detector in the NuMI cavern.
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
const double j
Definition: BetheBloch.cxx:29
int Plane() const
Return plane value.
Definition: PCHit.h:24
geo::View_t View() const
Return view.
Definition: PCHit.h:36
void SetMC(bool isMC=true)
Definition: RawDigit.h:106
Histograms used by attenuation calibration.
const std::string path
Definition: plot_BEN.C:43
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
TRandom3 r(0)
double GetPECorr(rb::CellHit const &cellhit, double w)
Float_t e
Definition: plot.C:35
std::string fPCHitLabel
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
enum BeamMode string