DEdx_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 // Class: DEdx
3 // Module Type: analyzer
4 // File: DEdx_module.cc
5 //
6 // This is a simple module that makes plot of dE/dx using 3D reconstructed tracks.
7 // It cab be run before each MC production on the sample of generated muon tracks.
8 // It can be also used on the data with reconstructed 3D tracks tagged as cosmic muons.
9 //
10 // Authors: Casandra Hazel Morris (University of Houston) and Dorota Stefan (CERN/NCBJ)
11 //
12 ////////////////////////////////////////////////////////////////////////
13 
23 
25 
28 #include "canvas/Persistency/Common/FindManyP.h"
34 #include "art_root_io/TFileService.h"
36 #include "fhiclcpp/ParameterSet.h"
38 
42 
43 #include "TH1.h"
44 #include "TTree.h"
45 
46 namespace proto
47 {
48  class DEdx;
49 }
50 
51 class proto::DEdx : public art::EDAnalyzer {
52 public:
53  explicit DEdx(fhicl::ParameterSet const & p);
54 
55  //plugins should not be copied or assigned.
56  DEdx(DEdx const &) = delete;
57  DEdx(DEdx &&) = delete;
58  DEdx & operator = (DEdx const &) = delete;
59  DEdx & operator = (DEdx &&) = delete;
60 
61  // Required functions.
62  void analyze(art::Event const & e) override;
63  void beginJob() override;
64  void endJob() override;
65 
66  void reconfigure(fhicl::ParameterSet const& p) ;
67 
68 private:
69 
70  // Declare member data here.
71 
72  // void CountdEdx(const std::vector < art::Ptr< recob::Hit > > & hits, const std::vector < recob::TrackHitMeta const* > & data); // MeV/cm
73  void CountdEdx(detinfo::DetectorClocksData const& clockData,
74  detinfo::DetectorPropertiesData const& detProp,
75  const std::vector < art::Ptr< recob::Hit > > & hits, const std::vector< recob::TrackHitMeta const* > & data);
76  void ResetVars();
77 
78  bool fCosmics;
79 
80  TTree *fTree;
81  TTree *fTreere;
82  int fRun;
83  int fEvent;
87 
88  double fdQdx;
89  double fdEdx;
90  double fdQ;
91  double fdx;
92 
93  // Module labels to get data products
98 
99 };
100 
102  :
103  EDAnalyzer(p),
104  fCalorimetryAlg(p.get<fhicl::ParameterSet>("CalorimetryAlg"))
105  // More initializers here.
106 {
107  reconfigure(p);
108 }
109 
111 {
113  fTreere = tfs->make<TTree>("entries","entries tree");
114  fTree = tfs->make<TTree>("dedx","dedx tree");
115 
116  fTreere->Branch("fRun", &fRun, "fRun/I");
117  fTreere->Branch("fEvent", &fEvent, "fEvent/I");
118  fTreere->Branch("fChosenView", &fChosenView, "fChosenView/I");
119  fTreere->Branch("fNumberOfTracks", &fNumberOfTracks, "fNumberOfTracks/I");
120  fTreere->Branch("fNumberOfTaggedTracks", &fNumberOfTaggedTracks, "fNumberOfTaggedTracks/I");
121 
122  fTree->Branch("fdQdx", &fdQdx, "fdQdx/D");
123  fTree->Branch("fdEdx", &fdEdx, "fdEdx/D");
124  fTree->Branch("fdQ", &fdQ, "fdQ/D");
125  fTree->Branch("fdx", &fdx, "fdx/D");
126 }
127 
129 {
130 }
131 
133 {
134  fChosenView = p.get<int>("ChosenView");
135  fCosmics = p.get<bool>("Cosmics");
136  fHitModuleLabel = p.get< std::string >("HitModuleLabel");
137  fTrackModuleLabel = p.get< std::string >("TrackModuleLabel");
138 }
139 
141 {
142  // Implementation of required member function here.
143  ResetVars();
144  fRun = e.run();
145  fEvent = e.id().event();
146 
147  // tracks
148  auto trkHandle = e.getValidHandle< std::vector<recob::Track> >(fTrackModuleLabel);
149  fNumberOfTracks = trkHandle->size();
150 
151  const art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trkHandle, e, fTrackModuleLabel);
152 
153  // Find the tagged tracks as cosmic muons
154  const art::FindManyP<anab::CosmicTag> ct(trkHandle, e, fTrackModuleLabel);
155 
156  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(e);
157  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService>()->DataFor(e, clockData);
158  if (fCosmics)
159  {
160  if (ct.isValid())
161  {
162  // loop over tracks
163  for (size_t t = 0; t < trkHandle->size(); ++t)
164  {
165  if (ct.at(t).size())
166  {
168 
169  auto vhit = fmthm.at(t);
170  auto vmeta = fmthm.data(t);
171 
172  CountdEdx(clockData, detProp, vhit, vmeta);
173  }
174  }
175  }
176  }
177  else if (fmthm.isValid())
178  {
179 
180  // loop over tracks
181  for (size_t t = 0; t < trkHandle->size(); ++t)
182  {
183  auto vhit = fmthm.at(t);
184  auto vmeta = fmthm.data(t);
185 
186  CountdEdx(clockData, detProp, vhit, vmeta);
187  }
188  }
189 
190  fTreere->Fill();
191 }
192 
194  detinfo::DetectorPropertiesData const& detProp,
195  const std::vector < art::Ptr< recob::Hit > > & hits,
196  const std::vector< recob::TrackHitMeta const* > & data) // MeV/cm
197 {
198  for (size_t h = 0; h < hits.size(); ++h)
199  {
200  unsigned short plane = hits[h]->WireID().Plane;
201  unsigned short time = hits[h]->PeakTime();
202 
203  if (plane == fChosenView)
204  {
205  double dqadc = hits[h]->Integral();
206  if (!std::isnormal(dqadc) || (dqadc < 0)) dqadc = 0.0;
207 
208  double t0 = 0;
209 
210  fdQdx = 0.0;
211  fdQ = dqadc;
212  fdx = data[h]->Dx();
213  if ((fdx > 0) && (fdQ > 0))
214  {
215  fdQdx = fdQ/fdx;
216  fdEdx = fCalorimetryAlg.dEdx_AREA(clockData, detProp, fdQdx, time, plane, t0);
217  if (fdEdx > 35) fdEdx = 35;
218  }
219  else if ((fdx == 0) && (fdQ > 0))
220  {
221 // std::cout << " the charge is positive but dx is equal zero " << std::endl;
222  }
223 
224  fTree->Fill();
225  }
226  }
227 }
228 
230 {
231  fdQdx = 0.0;
232  fdEdx = 0.0;
233  fdQ = 0.0;
234  fdx = 0.0;
235  fNumberOfTracks = 0;
237 }
238 
239 
code to link reconstructed objects back to the MC truth information
void ResetVars()
Definition: DEdx_module.cc:229
void analyze(art::Event const &e) override
Definition: DEdx_module.cc:140
double fdQ
Definition: DEdx_module.cc:90
std::string string
Definition: nybbler.cc:12
DEdx(fhicl::ParameterSet const &p)
Definition: DEdx_module.cc:101
struct vector vector
Class to keep data related to recob::Hit associated with recob::Track.
DEdx & operator=(DEdx const &)=delete
TTree * fTreere
Definition: DEdx_module.cc:81
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
void reconfigure(fhicl::ParameterSet const &p)
Definition: DEdx_module.cc:132
TTree * fTree
Definition: DEdx_module.cc:80
art framework interface to geometry description
const double e
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
double fdx
Definition: DEdx_module.cc:91
void CountdEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit > > &hits, const std::vector< recob::TrackHitMeta const * > &data)
Definition: DEdx_module.cc:193
double fdQdx
Definition: DEdx_module.cc:88
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::string fHitModuleLabel
Definition: DEdx_module.cc:94
void beginJob() override
Definition: DEdx_module.cc:110
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
p
Definition: test.py:223
size_t fNumberOfTracks
Definition: DEdx_module.cc:85
General LArSoft Utilities.
RunNumber_t run() const
Definition: DataViewImpl.cc:71
std::string fCalorimetryModuleLabel
Definition: DEdx_module.cc:96
double fdEdx
Definition: DEdx_module.cc:89
void endJob() override
Definition: DEdx_module.cc:128
Declaration of signal hit object.
Contains all timing reference information for the detector.
std::string fTrackModuleLabel
Definition: DEdx_module.cc:95
Provides recob::Track data product.
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
calo::CalorimetryAlg fCalorimetryAlg
Definition: DEdx_module.cc:97
EventNumber_t event() const
Definition: EventID.h:116
Access the description of detector geometry.
size_t fNumberOfTaggedTracks
Definition: DEdx_module.cc:86
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
Collection of Physical constants used in LArSoft.
EventID id() const
Definition: Event.cc:34