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