CosmicMetrics_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \brief Module that plots metrics from reconstructed cosmic ray data
3 /// \author gsdavies@iastate.edu
4 ////////////////////////////////////////////////////////////////////////////
5 
6 #include <string>
7 
8 // ROOT includes
9 #include "TH1D.h"
10 #include "TVector3.h"
11 
12 // Framework includes
20 
21 // NOvA includes
23 #include "Geometry/Geometry.h"
25 #include "NovaDatabase/Table.h"
26 #include "RawData/RawTrigger.h"
27 #include "RecoBase/CellHit.h"
28 #include "RecoBase/Track.h"
29 #include "Utilities/func/MathUtil.h"
30 
31 static const double ns = 1.0E-9;
32 
33 namespace comi {
34  // A module to plot cosmic metrics from reconstructed cosmic ray muons
35  class CosmicMetrics : public art::EDAnalyzer {
36  public:
37  explicit CosmicMetrics(fhicl::ParameterSet const& pset);
38  virtual ~CosmicMetrics();
39 
40  void analyze(art::Event const& evt);
41  void reconfigure(const fhicl::ParameterSet& pset);
42 
43  void beginJob();
44  void endJob();
45 
46  private:
47 
48  bool IsEntering(double x, double y, double z); ///< yes/no is "track" entering within a buffer window
49  std::string fRawDataLabel; ///< Where to find RawData
50  std::string fCellHitsListLabel; ///< Where to find CellHits
51  std::string fTrackLabel; ///< Where to find Tracks
52  std::string fSoftwareRelease; ///< Software release
53 
54  double fTrkLengthCut; ///< Minimum length of track to accept [cm]
55  double fMinTDC; ///< Minimum TDC value to accept track
56  double fMaxTDC; ///< Maximum TDC value to accept track
57  unsigned int fMuonNhit; ///< Required # hits on muon track
58  unsigned int fMuonNhitX; ///< # required in x view
59  unsigned int fMuonNhitY; ///< # required in y view
60  double fDetEdgeBuffer; ///< Detector edge buffer zone for entering cosmic ray track decision [cm]
61 
62  unsigned int fNevents; ///< total number of events counter
63  int fNTracks; ///< total number of reconstructed tracks counter, using time cuts
64  int fRun; ///< Run number
65  int fSubRun; ///< Subrun number
66  unsigned int fTriggerType;
67  int fWriteToDB;
68 
69  //......................................................................
70 
71  TH1D *fTrackLength; ///< length of track [cm]
72  TH1D *fCosCosmic; ///< cosine theta of cosmic track with respect to zenith
73  TH1D *fCosmicRate; ///< histogram holding single value for cosmic rate at end of job [Hz]
74  TH1D *fNEvents; ///< total number of events
75  TH1D *fTNS[2]; ///< distribution of XZ/YZ CellHit times in nano seconds
76 
77  };
78 }
79 
80 ////////////////////////////////////////////////////////////////////////
81 
82 namespace comi
83 {
84  //......................................................................
86  EDAnalyzer(pset),
87  fRawDataLabel (pset.get< std::string >("RawDataLabel") ),
88  fCellHitsListLabel(pset.get< std::string >("CellHitsListLabel") ),
89  fTrackLabel (pset.get< std::string >("TrackLabel") ),
90  fTrkLengthCut (pset.get< double >("TrkLengthCut") ),
91  fMinTDC (pset.get< double >("MinTDC") ),
92  fMaxTDC (pset.get< double >("MaxTDC") ),
93  fMuonNhit (pset.get< unsigned int>("MuonNhit") ),
94  fMuonNhitX (pset.get< unsigned int>("MuonNhitX") ),
95  fMuonNhitY (pset.get< unsigned int>("MuonNhitY") ),
96  fDetEdgeBuffer (pset.get< double >("DetEdgeBuffer") ),
97  fNevents(0), fNTracks(0), fRun(-999), fSubRun(-999), fTriggerType(-999)
98  {
99  reconfigure(pset);
100  }
101 
102  //..................................................................................
103 
105  {
106  }
107 
108  //......................................................................
109 
111  {
113 
114  fNEvents = tfs->make<TH1D>("fNEvents",";;Number of events",2,0,2);
115  fTrackLength = tfs->make<TH1D>("trackLength", ";Track Length (cm);" , 100, 0., 1000.);
116  fCosCosmic = tfs->make<TH1D>("coscosmic",";Cos #theta wrt zenith",1000,0.,1.);
117  fCosmicRate = tfs->make<TH1D>("cosmicrate",";Cosmic Rate (Hz)",1,0.,100000);
118 
119  for(int view = geo::kX; view <= geo::kY; ++view){
120  const TString viewStr = (view == geo::kX) ? "X" : "Y";
121  fTNS[view] = tfs->make<TH1D>("HitTNS"+viewStr, ";Time (ns);hits", 70, -100000.0, 600000.);
122  }
123  }
124 
125  //......................................................................
126 
127  ///
128  /// \brief The analysis of CosmicTrack reconstructed tracks in roder to
129  /// calculate cosmic ray muon rates, track lengths and more
130  /// cosmic ray metrics to be monitored for commissioning
131  ///
132  /// \param evt - Read-only access to the event data
133  ///
134 
136  {
137 
139 
140  ++fNevents; //Updating counter of number of events in sub-run
141  if(fRun < 0) fRun = evt.run();
142  if(fSubRun < 0) fSubRun = evt.subRun();
143  fNEvents->Fill(1);
144 
146  evt.getByLabel(fRawDataLabel, trigv);
147  const rawdata::RawTrigger& trig = (*trigv)[0];
149 
150  // now get the tracks that were made by CosmicTrack
151  // and check that there are tracks in the event
153  evt.getByLabel(fTrackLabel,trks);
154  if (trks->size()==0) {
155  mf::LogWarning("CosmicMetricsWarn") << "There are " << trks->size()
156  << " Tracks in this event, skip";
157  return;
158  }
159  LOG_DEBUG("CosmicMetrics") << "There are " << trks->size() << " Tracks in this event";
160 
161  const unsigned int trkMax = trks->size(); // maximum number of tracks in event
162  double cosCosmic = 0.;
163  int numtracks = 0 ;
164 
165  // Loop over tracks
166  for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
167 
168  bool isgoodXZtrack = false;
169  bool isgoodYZtrack = false;
170 
171  art::Ptr<rb::Track> track(trks, trkIdx); // ART-pointer to each track in event
172 
173  // Get vectors for the end points of the track
174  TVector3 start(track->Start());
175  TVector3 end(track->Stop());
176  TVector3 dir(track->Dir());
177  TVector3 length = end - start;
178 
179  if(end.Y() > start.Y()) std::swap(start, end); // force start point in Y to be greater than end point in Y for cosmics
180 
181  // Does this look like an entering track?
182  if (this->IsEntering(start.X(),start.Y(),start.Z()) == false) continue;
183 
184  const unsigned int hitMax = track->NCell(); // maximum number of CellHits on track
185  const unsigned int hitX = track->NXCell(); // maximum number of CellHits on track in XZ
186  const unsigned int hitY = track->NYCell(); // maximum number of CellHits on track in YZ
187 
188  // Do NOT select tracks if max number of hits or XZ- YZ-hits less than configurable values
189  if ( (hitMax < fMuonNhit) || (hitX < fMuonNhitX) || (hitY < fMuonNhitY) ) continue;
190  assert((track->NXCell() > 0) && (track->NYCell() > 0));
191 
192  double TNS[2] = {0.}; // TNS() value in each view
193  double TNStotal[2] = {0.}; // total TNS() in each view
194  double TNSaverage[2] = {0.}; // average TNS() in each view
195 
196  // Only select tracks that cross more than fTrkLengthCut [cm] in Z
197  if(fabs(end.Z() - start.Z()) > fTrkLengthCut){
198 
199  fTrackLength->Fill(length.Mag());
200 
201  // Loop through both views and grab every CellHit in this track
202  for(int view = geo::kX; view <= geo::kY; ++view){
203  // C++ is a bit strict about the conversions. Save typing.
204  const geo::View_t geoview = geo::View_t(view);
205  const int vhitMax = track->NCell(geoview);
206 
207  for(int hitIdx = 0; hitIdx < vhitMax; ++hitIdx){
208 
209  art::Ptr<rb::CellHit> chit = track->Cell(geoview, hitIdx);
210  fTNS[view]->Fill(chit->TNS());
211  TNS[view] = chit->TNS();
212  TNStotal[view]+=TNS[view];
213  } // hitIdx
214 
215  TNSaverage[view] = TNStotal[view]/track->NCell(geoview); // average TNS() CellHits in each view
216  } // view
217  } // track length cut
218 
219  if(TNSaverage[geo::kX] >= fMinTDC && TNSaverage[geo::kX] <= fMaxTDC){
220  isgoodXZtrack = true;
221  }
222  if(TNSaverage[geo::kY] >= fMinTDC && TNSaverage[geo::kY] <= fMaxTDC){
223  isgoodYZtrack = true;
224  }
225 
226  // Only count 3D tracks within configurable window size and with matching XZ/YZ 2D tracks
227  if(isgoodXZtrack == true && isgoodYZtrack == true){
228  numtracks++;
229  }
230 
231  // Cosmic track angular distribution with respect to zenith
232  cosCosmic = TMath::Abs(-dir.Y()/util::pythag(dir.X(),dir.Y(),dir.Z()));
233  fCosCosmic->Fill(cosCosmic);
234 
235  } // loop over tracks
236  fNTracks += numtracks;
237 
238  return;
239  }
240 
241  //......................................................................
242 
244  {
245  fTrkLengthCut = pset.get< double >("TrkLengthCut" );
246  fMinTDC = pset.get< double >("MinTDC" );
247  fMaxTDC = pset.get< double >("MaxTDC" );
248  fDetEdgeBuffer= pset.get< double >("DetEdgeBuffer");
249  fWriteToDB = pset.get< int >("WriteToDB");
250  }
251 
252  //......................................................................
253 
255  {
256 
257 
259 
260  float fLivetime = fNevents*((fMaxTDC-fMinTDC)*ns);
261  float fCosmicrate = fNTracks/(fNevents*((fMaxTDC-fMinTDC)*ns));
262  std::cout<< "There are: " << fNevents << " events in this job.\n";
263  std::cout<< "There are: " << fNTracks << " tracks selected in this job.\n";
264  std::cout<< "The livetime is " << fLivetime << "!\n";
265  std::cout<< "The cosmic rate is therefore " << fCosmicrate << " Hz \n";
266  std::cout<< "tracklengthmean="<<fTrackLength->GetMean(1)<<std::endl;
267  std::cout<< "tracklengthrms=" <<fTrackLength->GetRMS(1) <<std::endl;
268  fCosmicRate->Fill(fCosmicrate);
269 
270  //Get software release
271  fSoftwareRelease = getenv("SRT_BASE_RELEASE");
272 
273  if (fWriteToDB) {
274 
275  nova::database::Table* t = new nova::database::Table(std::string("/Commissioning/tables/CosmicMetrics.xml"));
276 
277  t->SetDetector(geom->DetId());
278  t->SetValidityRange("detector", novadaq::cnv::DetInfo::GetName(geom->DetId()));
279  t->SetValidityRange("run",fRun);
280  t->SetValidityRange("subrun",fSubRun);
281  t->SetValidityRange("triggertype",fTriggerType);
282  t->SetValidityRange("softwarerelease", fSoftwareRelease);
283  t->LoadFromDB();
284  if ( t->NRow() ){
285  std::cerr<< " Subrun metrics are already in the database "<<std::endl;
286  }
287  else{
288  std::cout<< " Writting subrun metrics to the database "<<std::endl;
289 
290  nova::database::Row* r = t->NewRow();
291 
292  r->Set("detector", novadaq::cnv::DetInfo::GetName(geom->DetId()));
293  r->Set("run", fRun);
294  r->Set("subrun", fSubRun);
295  r->Set("triggertype", fTriggerType);
296  r->Set("nevents", fNevents);
297  r->Set("ntracks", fNTracks);
298  r->Set("livetime", fLivetime);
299  r->Set("cosmicrate", fCosmicrate);
300  r->Set("tracklengthmean", fTrackLength->GetMean(1));
301  r->Set("tracklengthrms", fTrackLength->GetRMS(1));
302  r->Set("softwarerelease", fSoftwareRelease);
303 
304  t->AddRow(r);
305  delete r;
306  t->WriteToDB();
307  }
308  }
309  }
310 
311  //......................................................................
312 
313  bool CosmicMetrics::IsEntering(double x, double y, double z)
314  {
316  if (y > geom->DetHalfHeight()-fDetEdgeBuffer) return true;
317  if (fabs(x) > geom->DetHalfWidth()-fDetEdgeBuffer) return true;
318  if (z < fDetEdgeBuffer) return true;
319  if (z > geom->DetLength()-fDetEdgeBuffer) return true;
320  return false;
321  }
322 } // end namespace comi
323 
324 //////////////////////////////////////////////////////////////////////////
325 namespace comi
326 {
328 }
329 ////////////////////////////////////////////////////////////////////////
static std::string GetName(int id)
float TNS() const
Definition: CellHit.h:46
bool Set(std::string cname, T value)
Definition: Row.h:30
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
SubRunNumber_t subRun() const
Definition: Event.h:72
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::string fTrackLabel
Where to find Tracks.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH1D * fNEvents
total number of events
bool IsEntering(double x, double y, double z)
yes/no is "track" entering within a buffer window
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double DetLength() const
OStream cerr
Definition: OStream.cxx:7
std::string fCellHitsListLabel
Where to find CellHits.
TH1D * fTrackLength
length of track [cm]
unsigned int fMuonNhit
Required # hits on muon track.
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
CosmicMetrics(fhicl::ParameterSet const &pset)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
bool SetValidityRange(std::string cname, T start, T end)
Definition: Table.h:206
unsigned int fMuonNhitY
required in y view
uint8_t fTriggerMask_TriggerType
Definition: RawTrigger.h:43
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
double fMaxTDC
Maximum TDC value to accept track.
unsigned int fMuonNhitX
required in x view
TH1D * fTNS[2]
distribution of XZ/YZ CellHit times in nano seconds
bool WriteToDB(bool commit=true)
use commit=false if just testing
Definition: Table.cpp:1810
length
Definition: demo0.py:21
std::string getenv(std::string const &name)
Commissioning files to look at the quality of our data.
Definition: Cana_module.cc:39
T get(std::string const &key) const
Definition: ParameterSet.h:231
unsigned int fNevents
total number of events counter
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
static const double ns
Module that plots metrics from reconstructed cosmic ray data.
z
Definition: test.py:28
void AddRow(const Row *row)
Definition: Table.cpp:727
bool SetDetector(std::string det)
Definition: Table.cpp:513
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
TH1D * fCosmicRate
histogram holding single value for cosmic rate at end of job [Hz]
std::string fSoftwareRelease
Software release.
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
int fNTracks
total number of reconstructed tracks counter, using time cuts
std::string fRawDataLabel
Where to find RawData.
T * make(ARGS...args) const
double fTrkLengthCut
Minimum length of track to accept [cm].
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
struct Table Table
Definition: TexBuilder.h:2
double fMinTDC
Minimum TDC value to accept track.
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
TH1D * fCosCosmic
cosine theta of cosmic track with respect to zenith
double fDetEdgeBuffer
Detector edge buffer zone for entering cosmic ray track decision [cm].
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
RunNumber_t run() const
Definition: Event.h:77
void reconfigure(const fhicl::ParameterSet &pset)
void analyze(art::Event const &evt)
The analysis of CosmicTrack reconstructed tracks in roder to calculate cosmic ray muon rates...
Float_t track
Definition: plot.C:35
int fSubRun
Subrun number.
Encapsulate the geometry of one entire detector (near, far, ndos)