BNBInfo_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief This module reads IFDB DB and then stores BNB spill info
3 /// \author ranjan@fnal.gov
4 ////////////////////////////////////////////////////////////////////////
5 
6 //C++ includes
7 #include <iostream>
8 #include <fstream>
9 #include <math.h>
10 #include <string>
11 #include <time.h>
12 #include <fstream>
13 #include <sstream>
14 #include <iomanip>
15 #include <stdarg.h>
16 #include <stdlib.h>
17 #include <stdio.h>
18 #include <string.h>
19 #include <errno.h>
20 #include "ifbeam.h"
21 
22 // Framework includes
29 #include "fhiclcpp/ParameterSet.h"
31 
33 // NOvASoft includes
34 #include "NovaTimingUtilities/TimingUtilities.h"
35 #include "RawData/RawTrigger.h"
37 #include "SummaryData/POTSum.h"
38 #include "SummaryData/SpillData.h"
39 
40 //IFDB includes
41 #include "IFBeam_service.h"
42 #include "Munits.h"
43 #include "ifbeam.h"
44 //ROOT includes
45 
46 #include "TF1.h"
47 #include "TGraph.h"
48 #include "TMinuit.h"
49 // ROOT includes
50 #include "TTree.h"
51 #include "TH1F.h"
52 #include "TFile.h"
53 #include "TH2F.h"
54 #include "TGraphAsymmErrors.h"
55 #include "TGraph.h"
56 #include "TVector.h"
57 
58 
59 
60 
61 /// Information about the IFDB
62 namespace ifdb
63 {
64  class MIN : public art::EDProducer
65  {
66  public:
67  explicit MIN(fhicl::ParameterSet const& pset); // Required! explicit tag tells the compiler this is not a copy constructor
68  virtual ~MIN();
69 
70  void produce(art::Event& evt);
71 
72  void reconfigure(const fhicl::ParameterSet& pset);
73 
74  void beginSubRun(art::SubRun &sr);
75  void endSubRun(art::SubRun &sr);
76 
77 
78  private:
79  double fMinPOTCut;
80  double fMinHornICut;
81  double fMaxHornICut;
82  double fMaxDeltaTCut;
83  double fMinPosXCut;
84  double fMaxPosXCut;
85  double fMinPosYCut;
86  double fMaxPosYCut;
87  double fBNBOffset;
89 
91  std::unique_ptr<BeamFolder> bfp;
94  //bool useBPMS = false;
95  //bool useBeamP = false;
96  bool isMC;
97 
98 
99  TTree* fOutTree;
100  int channel;
101  double response;
102 
103 
105 
106  static constexpr double
107  PI = 3.14159265,
108  // Positions on BPM's along the beamline (in meters).
109  z_mw875 = 201.820648,
110  z_hp875 = 202.078003, // // HP875
111  z_vp875 = 202.281219, // // VP875
112  z_mw876 = 202.966766,
113 
114  z_vptg1 = 204.629608,
115  z_hptg1 = 204.833267,
116  z_vptg2 = 205.036835,
117  z_hptg2 = 205.240662,
118  z_mwtgt = 205.867645, // multiwire (profile monitor),
119  z_mtgt = 206.870895, // target,
120  z_hrndn = 208.317703,
121  hptg1station = 2.76,
123  vptg1station = 2.55,
124  vp875station = 0.20,
125  tgtstation = 4.80,
126 
127  hptg1offset = 0.46,
128  hp875offset = -3.40,
129  vptg1offset = 0.39,
130  vp875offset = 1.48,
131 //target LE and actrn1 are not actually used...
132  z_target_le = 0.0*Munits::foot, // LE location
133  actrn1 = 1.2467192*Munits::foot, // ACTRN1 OLD
134  wirespacing = 0.5*Munits::mm; // wire spacing on SWIC
135 
136  };
137 }
138 
139 
140 ////////////////////////////////////////////////////////////////////////
141 namespace ifdb
142 {
143  //..............................................................
145  fMinPOTCut (pset.get< double >("MinPOTCut") ),
146  fMinHornICut (pset.get< double >("MinHornICut") ),
147  fMaxHornICut (pset.get< double >("MaxHornICut") ),
148  fMaxDeltaTCut (pset.get< double >("MaxDeltaTCut") ),
149  fMinPosXCut (pset.get< double >("MinPosXCut") ),
150  fMaxPosXCut (pset.get< double >("MaxPosXCut") ),
151  fMinPosYCut (pset.get< double >("MinPosYCut") ),
152  fMaxPosYCut (pset.get< double >("MaxPosYCut") ),
153  fBNBOffset (pset.get< double >("BNBOffset") ),
154  pRawDataLabel (pset.get< std::string >("RawDataLabel") ),
155  bfp( ifbeam_handle->getBeamFolder(pset.get< std::string >("Bundle"), pset.get< std::string >("URL"), pset.get< double >("TimeWindow"))),
156  pPotLabel (pset.get< std::string >("PotLabel") ),
157  isMC(false)
158 
159  { bfp->set_epsilon(0.01); // how close in time does the spill time have to be from the DAQ time (in seconds).
160  // bfp->setValidWindow(0.01);
161  // tell the module what it is making as this is a EDProducer
162  produces< sumdata::SpillData >();
163  produces< sumdata::POTSum, art::InSubRun >();
164  // Be noisy to demonstrate what's happening
165  mf::LogInfo("MIN") << "In MIN Constructor\n";
166 
167  }
168 
169  //......................................................................
171  {
172  }
173 
174  //......................................................................
176  {
177  totpots = new sumdata::POTSum();
178 
180  fOutTree = tfs->make<TTree>("TimingCalib","Timing Calibaration");
181 
182  fOutTree->Branch("channel", &channel, "channel/I");
183  fOutTree->Branch("response", &response, "respone/I");
184 
185  }
186 
187  //......................................................................
189  {
192 
193  // If MC, get the old POTSum object that was filled in GENIEGen,
194  // and copy it to my current one
195  sr.getByLabel(pPotLabel,mcpots);
196  totpots->totpot = mcpots->totpot;
197  totpots->totgoodpot = mcpots->totgoodpot;
198  totpots->totspills = mcpots->totspills;
199  totpots->goodspills = mcpots->goodspills;
200  }
201 
202  std::unique_ptr< sumdata::POTSum > p(totpots);
203  sr.put(std::move(p));
204  }
205 
206  //......................................................................
208  {
209  fMinPOTCut = pset.get< double >("MinPOTCut" );
210  fMinHornICut = pset.get< double >("MinHornICut" );
211  fMaxHornICut = pset.get< double >("MaxHornICut" );
212  fMinPosXCut = pset.get< double >("MinPosXCut" );
213  fMaxPosXCut = pset.get< double >("MaxPosXCut" );
214  fMinPosYCut = pset.get< double >("MinPosYCut" );
215  fMaxPosYCut = pset.get< double >("MaxPosYCut" );
216  fMaxDeltaTCut = pset.get< double >("MaxDeltaTCut" );
217  fBNBOffset = pset.get< double >("BNBOffset" );
218  }
219 
220 
221 
222 
223 
224 
225  //......................................................................
227  {
228  std::unique_ptr<sumdata::SpillData> spilldata(new sumdata::SpillData());
229 
230  //count the total number of spills
231  // totpots->totspills++;
232 
233 
234  if(!evt.isRealData()){
235  // remove empty spill info write from here, move to GENIEGen and fill
236  isMC = true;
237 
238  // Ask MCTruth what I am
239 
240  // TODO: It would be ideal to get this from the Trigger information
241  // but currently it all reads back as TriggerType 0 ...
243  evt.getByLabel(pPotLabel, mcTruths);
244  if(mcTruths.failedToGet() || mcTruths->empty()) return;
245 
246  const simb::MCTruth& mcTruth = (*mcTruths)[0];
247 
248  fMCOrigin = mcTruth.Origin();
249 
250  return;
251  }
252 
253  //this will be the CORRECT UTC time we want ultimately
254  struct timespec unixTime;
255 
256  // Get trigger information to extract trigger time
257  auto trigv = evt.getValidHandle<std::vector<rawdata::RawTrigger> >(pRawDataLabel);
258  if( trigv.failedToGet() )
259  throw cet::exception("MIN") << "failed to get valid trigger handle";
260 
261  const rawdata::RawTrigger& trig = (*trigv)[0];
262 
263  if(trig.fTriggerMask_TriggerType != 1){
264  // This isn't a BNB trigger. Don't try to get beam info
265  return;
266  }
267 
269 
270  unsigned long int uevt_sec = unixTime.tv_sec; //utval >> 32 & 0xFFFFFFFFF;
271  unsigned long int uevt_nsec = unixTime.tv_nsec; //utval & 0xFFFFFFFFF;
272 
273  //DAQ event time
274  double DAQtime = uevt_sec + (0.000000001)*uevt_nsec;
275 
276 
277 
278  //torroid values and timestamps
279  double tor860=0.;
280  double tor875=0.;
281  double tor860time=0.;
282 
283  //for horn current
284  double thcurr=0.;
285 
286 
287  try{
288 
289  bfp->GetNamedData((DAQtime)-fBNBOffset,"E:TOR860@,E:TOR875", &tor860, &tor860time, &tor875);
290 
291  // This checks that the event time is close to the database time
292  // to ensure that we are grabbing the right data
293  double diff = 1e6; // a large number
294  diff = (DAQtime-fBNBOffset) - tor860time;
295 
296  spilldata->deltaspilltimensec = diff*1e9;
297 
298  unsigned long int time_closest_int = (int) tor860time;
299  double time_closest_ns = (tor860time - time_closest_int)*1e9;
300 
301  spilldata->spilltimesec = time_closest_int;
302  spilldata->spilltimensec = time_closest_ns;
303 
304  // If spill meets time requirement, calculate all quantities and
305  // fill appropriate fields
306  if(std::abs(spilldata->deltaspilltimensec) < fMaxDeltaTCut){
307 
308  double temppot = tor860;
309 
310  if(temppot < 0.02) temppot = tor875;
311 
312  // This guarantees that when toroid values go negative, total POT is not reduced
313  if(temppot < 0.0) temppot = 0.0;
314 
315  spilldata->spillpot = temppot;
316 
317  // Total POT
318  totpots->totpot += spilldata->spillpot;
319  totpots->totspills += 1;
320 
321  //These were 0 in old code, so set them to 0 here.
322  spilldata->gpsspilltimesec = 0.;
323  spilldata->gpsspilltimensec = 0.;
324 
325  // Horn current
326  bfp->GetNamedData(DAQtime-fBNBOffset,"E:THCURR", &thcurr);
327 
328  double ihorn = 0.;
329 
330  // Unlike NuMI, BNB gets horn current from a single input
331  ihorn += (thcurr - (+0.01));
332 
333  spilldata->hornI = ihorn;
334 
335 
336  // Horn polarity
337  spilldata->isRHC = false;
338  if(ihorn > 0){
339  spilldata->isRHC = true;
340  }
341 
342 
343 //===========Beam position==================
344 double HPTG1, VPTG1, HP875, VP875;
345 double xp, yp;
346 
347 
348  bfp->GetNamedData((DAQtime)-fBNBOffset,"E:HPTG1,E:VPTG1,E:HP875,E:VP875", &HPTG1, &VPTG1, &HP875, &VP875);
349 
350 
351  double x2 = HPTG1-hptg1offset;
352  double x1 = HP875-hp875offset;
353  double dx = x2-x1;
355  double y2 = VPTG1-vptg1offset;
356  double y1 = VP875-vp875offset;
357  double dy= y2-y1;
359 
360 
361 
362 
363  spilldata->posx = xp; //make it in mm
364  spilldata->posy = yp; //make it in mm
365 
366 
367 
368 
369 //==========================================
370 
371  spilldata->goodbeam = 0;
372 
373  // Apply good beam cuts
374  if((spilldata->spillpot > fMinPOTCut) // POT
375  && (spilldata->hornI > fMinHornICut
376  && spilldata->hornI < fMaxHornICut) // Horn I
377  && (spilldata->posx > fMinPosXCut
378  && spilldata->posx < fMaxPosXCut) // x position
379  && (spilldata->posy > fMinPosYCut
380  && spilldata->posy < fMaxPosYCut) // y position
381  && (std::abs(spilldata->deltaspilltimensec) < fMaxDeltaTCut)) // delta t
382  {
383  spilldata->goodbeam = 1;
384  } // Good beam cuts
385 
386 
387  // Total good POT
388  if(spilldata->goodbeam == 1){
389  totpots->totgoodpot += spilldata->spillpot;
390  totpots->goodspills += 1;
391  }
392 
393 
394  /* if(!(spilldata->spillpot>fMinPOTCut)){
395  mf::LogInfo("MIN") << "Failed POT cut" << std::endl;
396  }
397  if(!(spilldata->hornI>fMinHornICut && spilldata->hornI<fMaxHornICut)){
398  mf::LogInfo("MIN") << "Failed horn I cut" << std::endl;
399  }
400 
401  if(!(std::abs(spilldata->deltaspilltimensec)<1e9)){
402  mf::LogInfo("MIN") << "Failed time cut" << std::endl;
403  }
404  */
405 
406  } // if deltaspilltimensec < fMaxDeltaTCut
407 
408  else mf::LogInfo("MIN") << "Spill did not meet MaxDeltaT cut\n";
409 
410  }
411  catch (WebAPIException we){
412  // Eventually put in better error catching. What if each device isn't filled?
413  mf::LogError("MIN") << "Exception: " << we.what() << "\n";
414 
415 
416  //There are at least two types of exceptions. One, relatively infrequent case, which is
417  //from database flakes. This should be aborted on and retried. However, much more frequently,
418  //some of the variable information is missing from the record. This is unfortunate.
419  //We have implemented a nasty way to try to tell the difference between the two; it isn't
420  //super robust and there should be a better method. If you ever find one,
421  //please implement it.
422 
423  std::string exceptionString(we.what());
424 
425  std::size_t missVar = exceptionString.find("variable not found");
426  if(missVar != std::string::npos){
427  mf::LogInfo("MIN") << "***************We have a case of missing data\n"<<
428  "in the beam database! This isn't awesome, but \n"<<
429  "is probably ok and we have decided to not abort. \n"<<
430  "This decision was made on string comparisions and \n"<<
431  "is a bad hack.In the future, we should actually \n"<<
432  "throw different types of exceptions for different \n"<<
433  "types of problems!\n"<<
434  "***************"<<std::endl;
435  }
436  else{
437  std::abort();
438  }
439 
440  } // end try/catch
441 
442  evt.put(std::move(spilldata));
443 
444  } //end event loop
445 
446 
447 
449 
450 } // end namespace ifdb
451 ////////////////////////////////////////////////////////////////////////
virtual ~MIN()
void endSubRun(art::SubRun &sr)
static constexpr double vp875station
std::unique_ptr< BeamFolder > bfp
static constexpr double z_mwtgt
double fMinHornICut
Float_t y1[n_points_granero]
Definition: compare.C:5
double fMinPOTCut
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
static constexpr double hp875station
static constexpr double z_vp875
Float_t x1[n_points_granero]
Definition: compare.C:5
std::string pRawDataLabel
const char * p
Definition: xmltok.h:285
static constexpr double PI
simb::Origin_t Origin() const
Definition: MCTruth.h:73
enum simb::_ev_origin Origin_t
event origin types
std::string pPotLabel
void reconfigure(const fhicl::ParameterSet &pset)
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double fBNBOffset
static constexpr double z_mw876
This module reads IFDB DB and then stores BNB spill info.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
static constexpr double hptg1station
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
static constexpr double z_target_le
static constexpr double z_mw875
static constexpr double z_mtgt
float abs(float number)
Definition: d0nt_math.hpp:39
static constexpr double z_hptg1
static constexpr double vp875offset
void beginSubRun(art::SubRun &sr)
uint8_t fTriggerMask_TriggerType
Definition: RawTrigger.h:43
static constexpr Double_t foot
Definition: Munits.h:124
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void produce(art::Event &evt)
double response
double dy[NP][NC]
double dx[NP][NC]
double fMinPosYCut
double fMinPosXCut
static constexpr double z_hrndn
T get(std::string const &key) const
Definition: ParameterSet.h:231
double fMaxPosXCut
simb::Origin_t fMCOrigin
sumdata::POTSum * totpots
static constexpr double z_vptg1
static constexpr double hptg1offset
MIN(fhicl::ParameterSet const &pset)
static constexpr double vptg1offset
static constexpr double wirespacing
static constexpr Double_t mm
Definition: Munits.h:136
TTree * fOutTree
double fMaxDeltaTCut
art::ServiceHandle< ifbeam_ns::IFBeam > ifbeam_handle
double fMaxHornICut
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static constexpr double z_hptg2
double fMaxPosYCut
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
unsigned long long fTriggerTimingMarker_TimeStart
Definition: RawTrigger.h:38
int goodspills
Definition: POTSum.h:31
ProductID put(std::unique_ptr< PROD > &&)
static constexpr double hp875offset
static constexpr double actrn1
static constexpr double z_hp875
int totspills
Definition: POTSum.h:30
static constexpr double tgtstation
static constexpr double vptg1station
static constexpr double z_vptg2
Event generator information.
Definition: MCTruth.h:32
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
double totpot
normalized by 10^12 POT
Definition: POTSum.h:27
bool failedToGet() const
Definition: Handle.h:196
static constexpr Double_t sr
Definition: Munits.h:164
Beam neutrinos.
Definition: MCTruth.h:23