IFDBSpillInfo_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief This module reads IFDB DB and then stores NuMI spill info
3 /// \authors lgoodenough@anl.gov and sphanbudd@winona.edu
4 /// \modif brunetti@fnal.gov
5 ////////////////////////////////////////////////////////////////////////
6 
7 //C++ includes
8 #include <iostream>
9 #include <fstream>
10 #include <math.h>
11 #include <string>
12 #include <time.h>
13 
14 // Framework includes
21 #include "fhiclcpp/ParameterSet.h"
23 
24 // NOvASoft includes
25 #include "NovaTimingUtilities/TimingUtilities.h"
26 #include "RawData/RawTrigger.h"
28 #include "SummaryData/POTSum.h"
29 #include "SummaryData/SpillData.h"
31 
32 //IFDB includes
33 #include "IFBeam_service.h"
34 #include "Munits.h"
35 
36 //ROOT includes
37 #include "TF1.h"
38 #include "TGraph.h"
39 #include "TMinuit.h"
40 
41 
42 
43 /// Information about the IFDB
44 namespace ifdb
45 {
47  {
48  public:
49  explicit IFDBSpillInfo(fhicl::ParameterSet const& pset); // Required! explicit tag tells the compiler this is not a copy constructor
50  virtual ~IFDBSpillInfo();
51 
52  void produce(art::Event& evt);
53 
54  void reconfigure(const fhicl::ParameterSet& pset);
55 
56  void beginSubRun(art::SubRun &sr);
57  void endSubRun(art::SubRun &sr);
58 
59  void BpmProjection(std::vector<double> &, std::vector<double> &,
60  std::vector<double> &, std::vector<double> &,
61  std::vector<std::vector<double>>);
62 
63 
64  // Project Profile centers and widths to target.
65  // Return int giving status 0 if no device data
66  // Units are Munits.
67  int ProfileProjection(double &, double &, double &, double &,
68  std::vector<double>);
69 
70  double extrapolate_position(double, double, double, double, double);
71  double BpmAtTarget(double &xpmean, double &ypmean, double &xpstd, double &ypstd,
72  std::vector<double>TargBpmX, std::vector<double>TargBpmY,
73  std::vector<double>BpmIntX, std::vector<double>BpmIntY);
74 
75 
76  double GetStats(TGraph *, double &, double &);
77 
78  double GetGaussFit(double &, double &, std::vector<double>);
79 
80 
81  private:
82  double fMinPOTCut;
83  double fMinHornICut;
84  double fMaxHornICut;
85  double fHornICut0HC;
86  double fMinPosXCut;
87  double fMaxPosXCut;
88  std::vector<double> fMinPosYCutVec;
89  std::vector<double> fMaxPosYCutVec;
90  double fMinPosYCut;
91  double fMaxPosYCut;
92  double fMinWidthXCut;
93  double fMaxWidthXCut;
94  double fMinWidthYCut;
95  double fMaxWidthYCut;
96  double fMaxDeltaTCut;
98 
100  std::unique_ptr<BeamFolder> bfp;
103  bool useBPMS = false;
104  bool useBeamP = false;
105  bool useOnlyNuMI = true;
106  bool isMC;
107  bool fDataIsCalledMC; // For Data mixed into primary-stream MC..
108  // evt think it's not real data
110 
112 
113 
114  static constexpr double
115  PI = 3.14159265,
116  //new beam monitor positions after survey. Origin is between two budal monitors of the NOVA target
117  z_hp121 = -68.04458*Munits::foot,//-72.2309*Munits::foot, // HP121
118  z_vp121 = -66.99283*Munits::foot,//-71.3142*Munits::foot, // VP121
119  //pm121 is not used...
120  z_pm121 = -70.5267*Munits::foot, // PM121 OLD
121  z_hptgt = -31.25508*Munits::foot,//-33.1564*Munits::foot, // HPTGT
122  z_vptgt = -30.16533*Munits::foot,//-32.2397*Munits::foot, // VPTGT
123  z_pmtgt = -28.356*Munits::foot,//-31.5980*Munits::foot, // PMTGT
124  //target LE and actrn1 are not actually used...
125  z_target_le = 0.0*Munits::foot, // LE location
126  actrn1 = 1.2467192*Munits::foot, // ACTRN1 OLD
127  wirespacing = 0.5*Munits::mm; // wire spacing on SWIC
128 
129  };
130 }
131 
132 
133 ////////////////////////////////////////////////////////////////////////
134 namespace ifdb
135 {
136  //..............................................................
138  EDProducer(pset),
139  fMinPOTCut (pset.get< double >("MinPOTCut") ),
140  fMinHornICut (pset.get< double >("MinHornICut") ),
141  fMaxHornICut (pset.get< double >("MaxHornICut") ),
142  fHornICut0HC (pset.get< double >("HornICut0HC") ),
143  fMinPosXCut (pset.get< double >("MinPosXCut") ),
144  fMaxPosXCut (pset.get< double >("MaxPosXCut") ),
145  fMinPosYCutVec (pset.get< std::vector<double> >("MinPosYCut") ),
146  fMaxPosYCutVec (pset.get< std::vector<double> >("MaxPosYCut") ),
147  fMinWidthXCut (pset.get< double >("MinWidthXCut") ),
148  fMaxWidthXCut (pset.get< double >("MaxWidthXCut") ),
149  fMinWidthYCut (pset.get< double >("MinWidthYCut") ),
150  fMaxWidthYCut (pset.get< double >("MaxWidthYCut") ),
151  fMaxDeltaTCut (pset.get< double >("MaxDeltaTCut") ),
152  pRawDataLabel (pset.get< std::string >("RawDataLabel") ),
153  bfp( ifbeam_handle->getBeamFolder(pset.get< std::string >("Bundle"), pset.get< std::string >("URL"), pset.get< double >("TimeWindow"))),
154  pPotLabel (pset.get< std::string >("PotLabel") ),
155  useBPMS (pset.get< bool >("UseBPMS") ),
156  useBeamP (pset.get< bool >("UseBeamP") ),
157  useOnlyNuMI (pset.get< bool >("UseOnlyNuMITriggers")),
158  isMC(false),
159  fDataIsCalledMC(pset.get< bool >("DataIsCalledMC")),
160  fCosmicIsRHC (pset.get< bool >("CosmicIsRHC"))
161  {
162  // tell the module what it is making as this is a EDProducer
163  produces< sumdata::SpillData >();
164  produces< sumdata::POTSum, art::InSubRun >();
165  // Be noisy to demonstrate what's happening
166  mf::LogInfo("IFDBSpillInfo") << "In IFDBSpillInfo Constructor\n";
167  }
168 
169  //......................................................................
171  {
172  }
173 
174  //......................................................................
176  {
177  totpots = new sumdata::POTSum();
178  }
179 
180  //......................................................................
182  {
185 
186  // If MC, get the old POTSum object that was filled in GENIEGen,
187  // and copy it to my current one
188  sr.getByLabel(pPotLabel,mcpots);
189  totpots->totpot = mcpots->totpot;
190  totpots->totgoodpot = mcpots->totgoodpot;
191  totpots->totspills = mcpots->totspills;
192  totpots->goodspills = mcpots->goodspills;
193  }
194 
195  std::unique_ptr< sumdata::POTSum > p(totpots);
196  sr.put(std::move(p));
197  }
198 
199  //......................................................................
201  {
202  fMinPOTCut = pset.get< double >("MinPOTCut" );
203  fMinHornICut = pset.get< double >("MinHornICut" );
204  fMaxHornICut = pset.get< double >("MaxHornICut" );
205  fHornICut0HC = pset.get< double >("HornICut0HC" );
206  fMinPosXCut = pset.get< double >("MinPosXCut" );
207  fMaxPosXCut = pset.get< double >("MaxPosXCut" );
208  fMinPosYCut = pset.get< double >("MinPosYCut" );
209  fMaxPosYCut = pset.get< double >("MaxPosYCut" );
210  fMinWidthXCut = pset.get< double >("MinWidthXCut" );
211  fMaxWidthXCut = pset.get< double >("MaxWidthXCut" );
212  fMinWidthYCut = pset.get< double >("MinWidthYCut" );
213  fMaxWidthYCut = pset.get< double >("MaxWidthYCut" );
214  fMaxDeltaTCut = pset.get< double >("MaxDeltaTCut" );
215  }
216 
217  // Linear extrapolation to location z3
218  double IFDBSpillInfo::extrapolate_position(double t1, double z1, double t2, double z2, double z3)
219  {
220  return t1+(t2-t1)*(z3-z1)/(z2-z1);
221  }
222 
223 
224  // Calculate mean and standard deviation for points in TGraph to use as
225  // starting parameters values in Fit
226  double IFDBSpillInfo::GetStats(TGraph *prof, double &mean, double &stdev)
227  {
228 
229  double qx = 0, qtot = 0;
230 
231  double X = 0.0;
232  double Qx = 0.0;
233 
234  double minprof = 1E6;
235  double maxprof = -1E6;
236 
237  // Find min and max values of profile. Profile is upside down Gaussian
238  // with offset. All profile values have already been inverted to make
239  // Gaussian right-side up - offset is now negative.
240  // To get good starting param values, shift all profile values up by the
241  // most negative of all profile values, ie minprof, as found above.
242  for(int i = 0; i < prof->GetN(); ++i) {
243 
244  Qx = prof->GetY()[i];
245 
246  if (Qx <= minprof) minprof = Qx;
247  if (Qx >= maxprof) maxprof = Qx;
248 
249  }
250 
251  for(int i = 0; i < prof->GetN(); ++i) {
252  // For each point in profile, get wire position and charge values
253  X = prof->GetX()[i];
254  Qx = prof->GetY()[i];
255 
256  prof->SetPoint(i,X,Qx-minprof);
257  Qx = prof->GetY()[i];
258 
259  // Total charge on all wires
260  qtot += Qx;
261 
262  // Total charge-weighted position value
263  qx += Qx*X;
264  }
265 
266  // If total charge is 0, mean and stdev are 0
267  if(qtot == 0.0) {
268  mean = 0.0;
269  stdev = 0.0;
270  }
271 
272  // Otherwise, mean position is charge-weighted mean
273  else{
274 
275  mean = qx/qtot;
276 
277  // Charge-weighted variance
278  double var = 0.0;
279 
280  for(int i = 0; i < prof->GetN(); ++i) {
281 
282  X = prof->GetX()[i];
283  Qx = prof->GetY()[i];
284 
285  var += Qx*(X-mean)*(X-mean);
286  }
287 
288  var = var/qtot;
289 
290  if(var > 0.0){
291  stdev = std::sqrt(var);
292  }
293  else{
294  stdev = 0.0;
295  }
296  }
297 
298  return qtot;
299  } // GetStats
300 
301 
302  void IFDBSpillInfo::BpmProjection(std::vector<double> &xp, std::vector<double> &yp,
303  std::vector<double> &xi, std::vector<double> &yi,
304  std::vector<std::vector<double>> BPMS)
305  {
306  xp.clear();
307  yp.clear();
308  xi.clear();
309  yi.clear();
310 
311  // In NOvA era (for now) we will assume the following
312  // In MINOS era, target location was 3937 after time period 1155564632
313  double z_targ = 0;
314 
315  // Return shortest array length of all BPM positions and intensities used
316  int size = BPMS.size();
317  int n = BPMS[0].size();
318  for(int ind = 1; ind < size; ++ind) {
319  int n2 = BPMS[ind].size();
320  if(n2 < n) n = n2;
321  }
322 
323  int start = 0;
324  // Skip first BPMS value (average used in auto-tune) but only if
325  // there are multiple batches
326  if(n > 1) ++start;
327 
328  // Check to see if elements of BPMS[4] and BPMS[5] (BPM intensities) are 0.
329  // If they are, set position values to 0. Otherwise, use corresponding values of
330  // vertical and horizontal position monitors to extrapolate vertical and
331  // horizontal position at target. Simply transfer intensities over.
332  for(int ind = start; ind < n; ++ind) {
333 
334  if(BPMS[4][ind] == 0.0 || BPMS[5][ind]== 0.0) {
335  yp.push_back(0.0);
336  xp.push_back(0.0);
337  yi.push_back(0.0);
338  xi.push_back(0.0);
339  continue;
340  }
341 
342  // Returns transverse position (x or y) at z_targ assuming a linear
343  // extrapolation from z_vp121 to z_vptgt
344  yp.push_back
345  (extrapolate_position(BPMS[0][ind]*Munits::mm,z_vp121,BPMS[2][ind]*Munits::mm,z_vptgt,z_targ));
346  xp.push_back
347  (extrapolate_position(BPMS[1][ind]*Munits::mm,z_hp121,BPMS[3][ind]*Munits::mm,z_hptgt,z_targ));
348  // Returns the horizontal or vertical intensity of beam
349  yi.push_back(BPMS[4][ind]);
350  xi.push_back(BPMS[5][ind]);
351  }
352 
353  return;
354 
355  } // BpmProjection
356 
357 
358  int IFDBSpillInfo::ProfileProjection(double &x, double &y,
359  double &xstdev, double &ystdev,
360  std::vector<double> PMTGT)
361  {
362 
363  double xtgt = 0., ytgt = 0., xstdevtgt = 0., ystdevtgt = 0.;
364 
365  std::vector<double> hchannelstgt, vchannelstgt;
366 
367  // Need to add line here to return 0 if no data from devices
368 
369  // Get relevant channels from the device
370  for(int ind = 0; ind <= 47; ++ind) {
371 
372  hchannelstgt.push_back(PMTGT[103+ind]);
373  vchannelstgt.push_back(PMTGT[103+48+ind]);
374  }
375 
376  GetGaussFit(xtgt,xstdevtgt,hchannelstgt);
377  GetGaussFit(ytgt,ystdevtgt,vchannelstgt);
378 
379  // If values at MTGTDS are reasonable, set them to values at target
380  if(std::abs(xtgt) < 11*Munits::mm && std::abs(ytgt) < 11*Munits::mm) {
381  x = xtgt;
382  y = ytgt;
383  }
384 
385  // Standard deviation values are not extrapolated to target
386  xstdev = xstdevtgt;
387  ystdev = ystdevtgt;
388 
389  return 1;
390  }// ProfileProjection
391 
392 
393  double IFDBSpillInfo::GetGaussFit(double &mean, double &sigma,
394  std::vector<double> profile)
395  {
396  TGraph Prof;
397  Prof.Set(0);
398  //Prof.Clear("");
399 
400  int npoints = 0;
401 
402  // For each of 48 channels of profile monitor
403  for(int ch = 1; ch <= 48; ++ch) {
404 
405  // Get position of wire. Origin is halfway between channels 24 and 25
406  double X = (ch-24.5)*wirespacing;
407 
408  // Get the voltage value.
409  double Qx = profile[ch-1];
410 
411  // ensure that they are greater than or equal to zero
412  // NOTE: FOR NOW DON'T DO THIS - INSTEAD FIT TO AN OFFSET
413  // Qx = Qx > 0 ? Qx : 0;
414 
415  // Set values of points in TGraph, x=X/Munits::mm and y=-1*Qx/Munits::millivolt
416  // Invert voltage to make it positive in Gaussian peak - may lead to
417  // better fitting
418  Prof.SetPoint(npoints,X/Munits::mm,-1*Qx/Munits::millivolt);
419 
420  ++npoints;
421 
422  } // for channels 1 to 48
423 
424 
425  // Now that all of channel information has been loaded...
426 
427  // Remove dead and hot channels from profiles being fitted.
428  // I DON'T DO THIS BUT I WILL REMOVE DEAD CHANNELS (IE THOSE WITH ZERO READING)
429  // SuppressDeadHot(&xProf);
430 
431  for(int i = 0; i < Prof.GetN(); ++i) {
432  if (Prof.GetY()[i] == 0.0){
433  Prof.RemovePoint(i);
434  --i;
435  }
436  }
437 
438  // For some reason this is faster than if I use a TMath::Gaus:
439  TF1 gausF("gausF","[3]+([0]*exp((-1*(x-[1])*(x-[1]))/(2*[2]*[2]))/(sqrt(2*3.142)*[2]))");
440  //TF1 gausF("gausF","[3]+[0]*TMath::Gaus(x,[1],[2],1)");
441 
442  gausF.SetParName(0,"Area");
443  gausF.SetParName(1,"Mean");
444  gausF.SetParName(2,"Sigma");
445  gausF.SetParName(3,"Offset");
446 
447  // NOT SURE WE NEED THIS
448  //gausF.SetParLimits(0,0.00,100000000);
449  //gausF.SetParLimits(2,0,1000);
450 
451  // Gaussian fit parameters
452  double area = -9.;
453  double offset = -9.;
454  mean = -9., sigma = -9.;
455 
456  // Get statistical mean and standard deviation of profile distribution
457  double StatMean = 0.0, stdev = 0.0;
458  area = GetStats(&Prof,StatMean,stdev);
459  area*= Munits::millivolt;
460 
461 
462  // Check we have enough points and enough charge for a reasonable fit
463  if(Prof.GetN() > 5 && std::abs(area) > 100*Munits::millivolt &&
464  std::abs(StatMean) < 20*wirespacing/Munits::mm &&
466 
467  gausF.SetParameter(0,area/Munits::millivolt);
468  gausF.SetParameter(1,StatMean);
469  gausF.SetParameter(2,stdev/2);
470  //offset value should be fixed to the baseline, get avg value of the neagtive readings
471  int nbneg=0;
472  double avgneg=0.;
473  for(int i = 0; i < Prof.GetN(); ++i) {
474  if (Prof.GetY()[i] < 0.0){
475  avgneg=avgneg+Prof.GetY()[i];
476  nbneg++;
477  }
478  }
479  if(nbneg!=0){
480  offset=avgneg/nbneg;
481  }else{
482  // Starting offset value is a guess based on looking at data
483  offset = -3.0;
484  }
485  gausF.SetParameter(3,offset);
486 
487  if(Prof.Fit("gausF","q","",-24*wirespacing/Munits::mm,24*wirespacing/Munits::mm) == 0) {
488 
489  bool goodfit = true;
490 
491  // Are fit results physical?
492  if(gausF.GetParameter(0) < 10){
493  goodfit = false;
494  }
495 
496  if(std::abs(gausF.GetParameter(1)) > 24*wirespacing/Munits::mm){
497  goodfit = false;
498  }
499 
500  if(gausF.GetParameter(2) > 24*wirespacing/Munits::mm || gausF.GetParameter(2) < 0.0){
501  goodfit = false;
502  }
503 
504  if(goodfit) {
505  // Remember to correct fit result for area by dividing by "bin size" and correcting units
506  area = gausF.GetParameter(0)/wirespacing*Munits::mm*Munits::millivolt;
507  mean = gausF.GetParameter(1)*Munits::mm;
508  sigma = gausF.GetParameter(2)*Munits::mm;
509  offset = gausF.GetParameter(3)*Munits::millivolt;
510  } // if goodfit
511 
512  }
513  } // if we can make a reasonable fit
514 
515 // else{
516 // //mf::LogInfo("IFDBSpillInfo") << "POINTS DID NOT PASS BASIC CUTS in profile fit to Gaussian" << std::endl;
517 // }
518 
519  return area;
520  } // GetGaussFit
521 
522 
523  //......................................................................
525  {
526  std::unique_ptr<sumdata::SpillData> spilldata(new sumdata::SpillData());
527 
528  //count the total number of spills
529  // totpots->totspills++;
530 
531 
532  if(!evt.isRealData() && !fDataIsCalledMC){
533  // remove empty spill info write from here, move to GENIEGen and fill
534  isMC = true;
535  // Ask MCTruth what I am
536 
537  // TODO: It would be ideal to get this from the Trigger information
538  // but currently it all reads back as TriggerType 0 ...
540  evt.getByLabel(pPotLabel, mcTruths);
541  if(mcTruths.failedToGet() || mcTruths->empty()) return;
542 
543  const simb::MCTruth& mcTruth = (*mcTruths)[0];
544 
545  fMCOrigin = mcTruth.Origin();
547  spilldata->spillpot = 0;
548  spilldata->goodbeam = true;
549  spilldata->deltaspilltimensec = 0;
550  spilldata->hornI = 0;
551  spilldata->posx = 0;
552  spilldata->posy = 0;
553  spilldata->widthx = 0;
554  spilldata->widthy = 0;
555  spilldata->isRHC = fCosmicIsRHC;
556 
557  evt.put(std::move(spilldata));
558  }
559 
560  return;
561  }
562 
563  //this will be the CORRECT UTC time we want ultimately
564  struct timespec unixTime;
565 
566  // Get trigger information to extract trigger time
567  auto trigv = evt.getValidHandle<std::vector<rawdata::RawTrigger> >(pRawDataLabel);
568  if( trigv.failedToGet() )
569  throw cet::exception("IFDBSpillInfo") << "failed to get valid trigger handle";
570 
571  const rawdata::RawTrigger& trig = (*trigv)[0];
572 
573  if(trig.fTriggerMask_TriggerType != 0 && useOnlyNuMI){
574  // This isn't a NuMI trigger. Don't try to get beam info
575  // make a product to fix a horn current if these are cosmics
576  if (trig.fTriggerMask_TriggerType == 2){
577  spilldata->spillpot = 0;
578  spilldata->goodbeam = true;
579  spilldata->deltaspilltimensec = 0;
580  spilldata->hornI = 0;
581  spilldata->posx = 0;
582  spilldata->posy = 0;
583  spilldata->widthx = 0;
584  spilldata->widthy = 0;
585  spilldata->isRHC = fCosmicIsRHC;
586 
587  evt.put(std::move(spilldata));
588  }
589 
590  return;
591  }
592 
594 
595  unsigned long int uevt_sec = unixTime.tv_sec; //utval >> 32 & 0xFFFFFFFFF;
596  unsigned long int uevt_nsec = unixTime.tv_nsec; //utval & 0xFFFFFFFFF;
597 
598  //DAQ event time
599  double DAQtime = uevt_sec + (0.000000001)*uevt_nsec;
600 
601  //torroid values and timestamps
602  double trtgtd=0.,tr101d=0.;
603  double trtgtdtime=0.;
604 
605  //for horn current
606  double nslina=0.,nslinb=0.,nslinc=0.,nslind=0.;
607  //for horn polarity
608  double horndir=0.;
609 
611  // Set posy cut based on run number.
612  // Beam vertical position changed starting period 11.
613  if ( (rh->Detector() == "NearDet" && evt.run() < 13887) ||
614  (rh->Detector() == "FarDet" && evt.run() < 38301) ) {
615  fMinPosYCut = fMinPosYCutVec.at(0);
616  fMaxPosYCut = fMaxPosYCutVec.at(0);
617  }
618  else {
619  fMinPosYCut = fMinPosYCutVec.at(1);
620  fMaxPosYCut = fMaxPosYCutVec.at(1);
621  }
622 
623  try{
624 
625  bfp->GetNamedData(DAQtime,"E:TRTGTD@,E:TR101D", &trtgtd, &trtgtdtime, &tr101d);
626 
627  // This checks that the event time is close to the database time
628  // to ensure that we are grabbing the right data
629  double diff = 1e6; // a large number
630  diff = DAQtime - trtgtdtime;
631 
632  spilldata->deltaspilltimensec = diff*1e9;
633 
634  unsigned long int time_closest_int = (int) trtgtdtime;
635  double time_closest_ns = (trtgtdtime - time_closest_int)*1e9;
636 
637  spilldata->spilltimesec = time_closest_int;
638  spilldata->spilltimensec = time_closest_ns;
639 
640  // If spill meets time requirement, calculate all quantities and
641  // fill appropriate fields
642  if(std::abs(spilldata->deltaspilltimensec) < fMaxDeltaTCut){
643 
644  double temppot = trtgtd;
645 
646  if(temppot < 0.02) temppot = tr101d;
647 
648  // This guarantees that when toroid values go negative, total POT is not reduced
649  if(temppot < 0.0) temppot = 0.0;
650 
651  spilldata->spillpot = temppot;
652 
653  // Total POT
654  totpots->totpot += spilldata->spillpot;
655  totpots->totspills += 1;
656 
657  //These were 0 in old code, so set them to 0 here.
658  spilldata->gpsspilltimesec = 0.;
659  spilldata->gpsspilltimensec = 0.;
660 
661  // Horn polarity: note, this is a fairly new variable, so it is not
662  // available in older files
663  bfp->GetNamedData(DAQtime,"E:HRNDIR",&horndir);
664 
665  spilldata->isRHC = rh->IsRHC();
666  /* if (spilldata->isRHC)
667  std::cout << "********** IsRHC **********" << std::endl;
668  else
669  std::cout << "********** IsFHC **********" << std::endl;
670  */
671 
672  /*
673  spilldata->isRHC = false;
674  if(horndir > 0){
675  spilldata->isRHC = true;
676  }
677  */
678 
679  // Horn current
680  bfp->GetNamedData(DAQtime,"E:NSLINA,E:NSLINB,E:NSLINC,E:NSLIND",&nslina,
681  &nslinb, &nslinc, &nslind);
682 
683  double ihorn = 0.;
684 
685  // Updated numbers from Jim Hylen, email Oct 4, 2013
686  ihorn += (nslina - (+0.01)) / 0.9951;
687  ihorn += (nslinb - (-0.14)) / 0.9957;
688  ihorn += (nslinc - (-0.05)) / 0.9965;
689  ihorn += (nslind - (-0.07)) / 0.9945;
690 
691  spilldata->hornI = ihorn;
692 
693  // Horn Off, flagged when abs(ihorn) < 1
694  spilldata->is0HC = false;
695  if(std::abs(spilldata->hornI) < fHornICut0HC){
696  spilldata->isRHC = rh->IsRHC(); //override this just in case horndir > 0 but horn current still within 0HC cuts
697  spilldata->is0HC = true;
698  }
699 
700 
701  // Beam position
702  if(useBPMS == true){
703  std::vector<double> VP121, HP121, VPTGT, HPTGT, VITGT, HITGT;
704  std::vector<std::vector<double>> BPMS;
705 
706  // BPM Positions and Intensities - they each have 7 elements
707  // First is an average value of a few batches (often 2,3,4)
708  // used for auto-tuning, so we should disregard it
709  VP121 = bfp->GetNamedVector(DAQtime,"E:VP121[]");
710  HP121 = bfp->GetNamedVector(DAQtime,"E:HP121[]");
711  VPTGT = bfp->GetNamedVector(DAQtime,"E:VPTGT[]");
712  HPTGT = bfp->GetNamedVector(DAQtime,"E:HPTGT[]");
713  VITGT = bfp->GetNamedVector(DAQtime,"E:VITGT[]");
714  HITGT = bfp->GetNamedVector(DAQtime,"E:HITGT[]");
715 
716  BPMS.push_back(VP121);
717  BPMS.push_back(HP121);
718  BPMS.push_back(VPTGT);
719  BPMS.push_back(HPTGT);
720  BPMS.push_back(VITGT);
721  BPMS.push_back(HITGT);
722 
723  std::vector<double> xp, yp, xi, yi;
724  std::vector<double> xpmm, ypmm;
725  double totxi=0;
726  double totyi=0;
727  double convx=1;
728  double convy=1;
729  // Use Beam Position Monitors to estimate x and y position of
730  // beam (xp,yp) and intensities (xi,yi) at target.
731  // Do this for all six readings of positions monitors
732  BpmProjection(xp,yp,xi,yi,BPMS);
733  //positions in mm
734  for(int i=0;i<6;i++){
735  xpmm.push_back(xp.at(i)*1000.);
736  ypmm.push_back(yp.at(i)*1000.);
737  }
738  spilldata->bposx=xpmm;
739  spilldata->bposy=ypmm;
740  //convert batch by batch intensities in POT units
741  for(int i=0;i<6;i++){
742  totxi+=xi[i];
743  totyi+=yi[i];
744  }
745  convx=totxi/temppot;
746  convy=totyi/temppot;
747  for(int i=0;i<6;i++){
748  xi[i]=xi[i]/convx;
749  yi[i]=yi[i]/convy;
750  }
751  //units will be e12
752  spilldata->intx=xi;
753  spilldata->inty=yi;
754  double xpmean = -9., ypmean = -9., xpstdev = -9., ypstdev=-9.;
755  // Use 6 position and intensity values to get mean and stdev
756  // values of positions
757  BpmAtTarget(xpmean,ypmean,xpstdev,ypstdev,xp,yp,xi,yi);
758  spilldata->posx = xpmean*1000.; //make it in mm
759  spilldata->posy = ypmean*1000.; //make it in mm
760 
761  } //if useBPMS
762 
763 
764  if(useBeamP == true){
765 
766  //BPM Positions
767  std::vector<double> MTGTDS;
768  MTGTDS = bfp->GetNamedVector(DAQtime,"E:MTGTDS[]");
769 
770  double x = -9., y = -9., xstdev= -9., ystdev= -9.;
771  // Use Beam Profile Monitors to estimate x and y widths and positions
772  // of beam at target
773  ProfileProjection(x,y,xstdev,ystdev,MTGTDS);
774 
775  spilldata->widthx = xstdev*1000.; //make it in mm
776  spilldata->widthy = ystdev*1000.; //make it in mm
777 
778  } //if useBeamP
779 
780  spilldata->goodbeam = 0;
781 
782  // Apply good beam cuts
783  if((spilldata->spillpot > fMinPOTCut) // POT
784  && (spilldata->hornI > fMinHornICut
785  && spilldata->hornI < fMaxHornICut) // Horn I
786  && (spilldata->posx > fMinPosXCut
787  && spilldata->posx < fMaxPosXCut) // x position
788  && (spilldata->posy > fMinPosYCut
789  && spilldata->posy < fMaxPosYCut) // y position
790  && (spilldata->widthx > fMinWidthXCut
791  && spilldata->widthx < fMaxWidthXCut) // x width
792  && (spilldata->widthy > fMinWidthYCut
793  && spilldata->widthy < fMaxWidthYCut) // y width
794  && (std::abs(spilldata->deltaspilltimensec) < fMaxDeltaTCut)) // delta t
795  {
796  spilldata->goodbeam = 1;
797  } // Good beam cuts
798 
799 
800  // Total good POT
801  if(spilldata->goodbeam == 1){
802  totpots->totgoodpot += spilldata->spillpot;
803  totpots->goodspills += 1;
804  }
805 
806  /*
807  if(!(spilldata->spillpot>fMinPOTCut)){
808  mf::LogInfo("IFDBSpillInfo") << "Failed POT cut" << std::endl;
809  }
810  if(!(spilldata->hornI>fMinHornICut && spilldata->hornI<fMaxHornICut)){
811  mf::LogInfo("IFDBSpillInfo") << "Failed horn I cut" << std::endl;
812  }
813 
814  if(!(spilldata->posx>fMinPosXCut && spilldata->posx<fMaxPosXCut)){
815  mf::LogInfo("IFDBSpillInfo") << "Failed pos x cut" << std::endl;
816  }
817 
818  if(!(spilldata->posy>fMinPosYCut && spilldata->posy<fMaxPosYCut)){
819  mf::LogInfo("IFDBSpillInfo") << "Failed pos y cut" << std::endl;
820  }
821 
822  if(!(spilldata->widthx>fMinWidthXCut && spilldata->widthx<fMaxWidthXCut)){
823  mf::LogInfo("IFDBSpillInfo") << "Failed width x cut" << std::endl;
824  }
825 
826  if(!(spilldata->widthy>fMinWidthYCut && spilldata->widthy<fMaxWidthYCut)){
827  mf::LogInfo("IFDBSpillInfo") << "Failed width y cut" << std::endl;
828  }
829 
830  if(!(std::abs(spilldata->deltaspilltimensec)<1e9)){
831  mf::LogInfo("IFDBSpillInfo") << "Failed time cut" << std::endl;
832  }
833  */
834 
835  } // if deltaspilltimensec < fMaxDeltaTCut
836 
837  else mf::LogInfo("IFDBSpillInfo") << "Spill did not meet MaxDeltaT cut\n";
838 
839  }
840  catch (const WebAPIException &we){
841  // Eventually put in better error catching. What if each device isn't filled?
842  mf::LogError("IFDBSpillInfo") << "Exception: " << we.what() << "\n";
843 
844  //There are at least two types of exceptions. One, relatively infrequent case, which is
845  //from database flakes. This should be aborted on and retried. However, much more frequently,
846  //some of the variable information is missing from the record. This is unfortunate.
847  //We have implemented a nasty way to try to tell the difference between the two; it isn't
848  //super robust and there should be a better method. If you ever find one,
849  //please implement it.
850 
851  std::string exceptionString(we.what());
852 
853  std::size_t missVar = exceptionString.find("variable not found");
854  if(missVar != std::string::npos){
855  mf::LogInfo("IFDBSpillInfo") << "***************We have a case of missing data\n"<<
856  "in the beam database! This isn't awesome, but \n"<<
857  "is probably ok and we have decided to not abort. \n"<<
858  "This decision was made on string comparisions and \n"<<
859  "is a bad hack.In the future, we should actually \n"<<
860  "throw different types of exceptions for different \n"<<
861  "types of problems!\n"<<
862  "***************"<<std::endl;
863  }
864  else if(useOnlyNuMI){
865  std::abort();
866  }
867 
868  } // end try/catch
869 
870  evt.put(std::move(spilldata));
871 
872  } //end event loop
873 
874 
875  double IFDBSpillInfo::BpmAtTarget(double &xpmean, double &ypmean,
876  double &xpstdev, double &ypstdev,
877  std::vector<double>TargBpmX,
878  std::vector<double>TargBpmY,
879  std::vector<double>BpmIntX,
880  std::vector<double>BpmIntY)
881  {
882  //xpmean = ypmean = xpstdev = ypstdev = 0.0;
883  double IXtot = 0.0, IYtot = 0.0;
884  double ix = 0.0, iy = 0.0, ix2 = 0.0, iy2 = 0.0;
885 
886  // For intensity-weighted position values...
887  // Remember, TargBpmX, TargBpmY, BpmIntX and BpmIntY have one fewer element than
888  // the devices because first device value was an average used for beam tuning
889  for(size_t ind = 0; ind < BpmIntX.size(); ++ind) {
890 
891  // If average intensity is 0, skip to next value
892  if(BpmIntX[ind] == 0.0 || BpmIntY[ind] == 0.0) continue;
893 
894  double X = TargBpmX[ind]; //xp[ind]
895  double Y = TargBpmY[ind]; //yp[ind]
896  double IX = BpmIntX[ind]; //xi[ind]
897  double IY = BpmIntX[ind]; //yi[ind]
898 
899  // if I>0, set I=I, else set I=0
900  IX = IX > 0 ? IX : 0;
901  IY = IY > 0 ? IY : 0;
902 
903  // sum intensities
904  IXtot += IX;
905  IYtot += IY;
906 
907  // intensity-weighted position value
908  ix += IX*X;
909  iy += IY*Y;
910 
911  ix2 += IX*X*X;
912  iy2 += IY*Y*Y;
913 
914  } // end for-loop over elements of xp, yp, xi and yi
915 
916  if(IXtot <= 0.0) return 0.0;
917  if(IYtot <= 0.0) return 0.0;
918 
919  // intensity-weighted average positions
920  xpmean = ix/IXtot;
921  ypmean = iy/IYtot;
922 
923  // standard deviation
924  double xvar = ix2/IXtot-xpmean*xpmean;
925  if(xvar > 0.0) xpstdev = sqrt(xvar);
926  double yvar = iy2/IYtot-ypmean*ypmean;
927  if(yvar > 0.0) ypstdev = sqrt(yvar);
928 
929  return xpmean;
930 
931  } // BpmAtTarget
932 
933 
935 
936 } // end namespace ifdb
937 ////////////////////////////////////////////////////////////////////////
std::unique_ptr< BeamFolder > bfp
IFDBSpillInfo(fhicl::ParameterSet const &pset)
def stdev(lst)
Definition: HandyFuncs.py:286
TH2 * rh
Definition: drawXsec.C:5
void BpmProjection(std::vector< double > &, std::vector< double > &, std::vector< double > &, std::vector< double > &, std::vector< std::vector< double >>)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
simb::Origin_t Origin() const
Definition: MCTruth.h:74
enum simb::_ev_origin Origin_t
event origin types
static constexpr double z_vp121
static constexpr double z_vptgt
var_value< double > var
Definition: StanTypedefs.h:14
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
This module reads IFDB DB and then stores BNB spill info.
int ProfileProjection(double &, double &, double &, double &, std::vector< double >)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
float abs(float number)
Definition: d0nt_math.hpp:39
static constexpr double z_pmtgt
double extrapolate_position(double, double, double, double, double)
Float_t Y
Definition: plot.C:38
sumdata::POTSum * totpots
uint8_t fTriggerMask_TriggerType
Definition: RawTrigger.h:43
static constexpr Double_t foot
Definition: Munits.h:124
static constexpr double wirespacing
bool isRealData() const
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
std::void_t< T > n
single particles thrown at the detector
Definition: MCTruth.h:26
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
static constexpr double actrn1
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:491
caf::StandardRecord * sr
static constexpr double z_hptgt
std::vector< double > fMinPosYCutVec
RunNumber_t run() const
double t2
double sigma(TH1F *hist, double percentile)
static constexpr Double_t mm
Definition: Munits.h:136
double GetStats(TGraph *, double &, double &)
double GetGaussFit(double &, double &, std::vector< double >)
static constexpr double PI
void endSubRun(art::SubRun &sr)
unsigned long long fTriggerTimingMarker_TimeStart
Definition: RawTrigger.h:38
int goodspills
Definition: POTSum.h:31
static constexpr double z_target_le
double BpmAtTarget(double &xpmean, double &ypmean, double &xpstd, double &ypstd, std::vector< double >TargBpmX, std::vector< double >TargBpmY, std::vector< double >BpmIntX, std::vector< double >BpmIntY)
void reconfigure(const fhicl::ParameterSet &pset)
int totspills
Definition: POTSum.h:30
void beginSubRun(art::SubRun &sr)
void produce(art::Event &evt)
Event generator information.
Definition: MCTruth.h:32
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
Float_t X
Definition: plot.C:38
std::vector< double > fMaxPosYCutVec
static constexpr Double_t millivolt
Definition: Munits.h:332
double totpot
normalized by 10^12 POT
Definition: POTSum.h:27
static constexpr double z_hp121
bool failedToGet() const
Definition: Handle.h:190
art::ServiceHandle< ifbeam_ns::IFBeam > ifbeam_handle
Cosmic rays.
Definition: MCTruth.h:24
static constexpr double z_pm121
Beam neutrinos.
Definition: MCTruth.h:23
std::string Detector() const
Definition: RunHistory.h:376
enum BeamMode string