BNBInfo_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   This module reads IFDB DB and then stores BNB spill info
/// \author ranjan@fnal.gov
////////////////////////////////////////////////////////////////////////

//C++ includes 
#include <iostream>
#include <fstream>
#include <math.h>
#include <string>
#include <time.h> 
#include <fstream>
#include <sstream>
#include <iomanip>
#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include "ifbeam.h" 

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "art/Framework/Services/Optional/TFileService.h"
// NOvASoft includes
#include "NovaTimingUtilities/TimingUtilities.h"
#include "RawData/RawTrigger.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "SummaryData/POTSum.h"
#include "SummaryData/SpillData.h"

//IFDB includes
#include "IFBeam_service.h"
#include "Munits.h"
#include "ifbeam.h" 
//ROOT includes

#include "TF1.h"
#include "TGraph.h"
#include "TMinuit.h"
// ROOT includes
#include "TTree.h"
#include "TH1F.h"
#include "TFile.h"
#include "TH2F.h"
#include "TGraphAsymmErrors.h"
#include "TGraph.h"
#include "TVector.h"




/// Information about the IFDB
namespace ifdb 
{
  class MIN : public art::EDProducer 
  {
  public:
    explicit MIN(fhicl::ParameterSet const& pset); // Required! explicit tag tells the compiler this is not a copy constructor      
    virtual ~MIN();                        
    
    void produce(art::Event& evt);
    
    void reconfigure(const fhicl::ParameterSet& pset);
    
    void beginSubRun(art::SubRun &sr);
    void endSubRun(art::SubRun &sr);
          
    
  private:
    double fMinPOTCut;
    double fMinHornICut;
    double fMaxHornICut;
    double fMaxDeltaTCut;
    double fMinPosXCut;
    double fMaxPosXCut;
    double fMinPosYCut;
    double fMaxPosYCut;
    double fBNBOffset;
    std::string pRawDataLabel;
    
    art::ServiceHandle<ifbeam_ns::IFBeam> ifbeam_handle;
    std::unique_ptr<BeamFolder> bfp;
    sumdata::POTSum *totpots;
    std::string pPotLabel;
    //bool useBPMS = false;
    //bool useBeamP = false;
    bool isMC;

 
      TTree* fOutTree;
      int channel;
      double response;


    simb::Origin_t fMCOrigin;
    
    static constexpr double
      PI          = 3.14159265,
    // Positions on BPM's along the beamline (in meters).
      z_mw875     = 201.820648,
      z_hp875     = 202.078003, // // HP875
      z_vp875     = 202.281219, // // VP875
      z_mw876     = 202.966766,
      
      z_vptg1     = 204.629608,
      z_hptg1     = 204.833267,  
      z_vptg2     = 205.036835,
      z_hptg2     = 205.240662,
      z_mwtgt     = 205.867645, // multiwire (profile monitor),
      z_mtgt      = 206.870895, // target,
      z_hrndn     = 208.317703,
        hptg1station = 2.76,
        hp875station = 0.0,
        vptg1station = 2.55,
        vp875station = 0.20,
        tgtstation   = 4.80,

        hptg1offset = 0.46,
        hp875offset = -3.40,
        vptg1offset = 0.39,
        vp875offset = 1.48,
//target LE and actrn1 are not actually used...
      z_target_le = 0.0*Munits::foot,       // LE location
      actrn1      = 1.2467192*Munits::foot, // ACTRN1 OLD       
      wirespacing = 0.5*Munits::mm;         // wire spacing on SWIC
                                
  };
}


////////////////////////////////////////////////////////////////////////
namespace ifdb
{
  //..............................................................
  MIN::MIN(fhicl::ParameterSet const& pset) :
    fMinPOTCut     (pset.get< double      >("MinPOTCut")     ),
    fMinHornICut   (pset.get< double      >("MinHornICut")   ),
    fMaxHornICut   (pset.get< double      >("MaxHornICut")   ),
    fMaxDeltaTCut  (pset.get< double      >("MaxDeltaTCut")  ),
    fMinPosXCut    (pset.get< double      >("MinPosXCut")    ),
    fMaxPosXCut    (pset.get< double      >("MaxPosXCut")    ),
    fMinPosYCut    (pset.get< double      >("MinPosYCut")    ),
    fMaxPosYCut    (pset.get< double      >("MaxPosYCut")    ),
    fBNBOffset     (pset.get< double      >("BNBOffset")     ),
    pRawDataLabel  (pset.get< std::string >("RawDataLabel")  ),
    bfp( ifbeam_handle->getBeamFolder(pset.get< std::string >("Bundle"), pset.get< std::string >("URL"), pset.get< double >("TimeWindow"))),
    pPotLabel      (pset.get< std::string >("PotLabel")      ),
    isMC(false)

  { bfp->set_epsilon(0.01); // how close in time does the spill time have to be from the DAQ time (in seconds). 
   // bfp->setValidWindow(0.01);
    // tell the module what it is making as this is a EDProducer
    produces< sumdata::SpillData >();
    produces< sumdata::POTSum, art::InSubRun >();
    // Be noisy to demonstrate what's happening
    mf::LogInfo("MIN") << "In MIN Constructor\n";
  
  }
  
  //......................................................................
  MIN::~MIN()
  {
  }
  
  //......................................................................
  void MIN::beginSubRun(art::SubRun &sr)
  {
    totpots = new sumdata::POTSum(); 

    art::ServiceHandle<art::TFileService> tfs;
 fOutTree = tfs->make<TTree>("TimingCalib","Timing Calibaration");

    fOutTree->Branch("channel", &channel, "channel/I");
    fOutTree->Branch("response", &response, "respone/I");

  }
  
  //......................................................................
  void MIN::endSubRun(art::SubRun &sr)
  {
    if(isMC && fMCOrigin == simb::kBeamNeutrino){
      art::Handle< sumdata::POTSum > mcpots;
      
      // If MC, get the old POTSum object that was filled in GENIEGen, 
      // and copy it to my current one
      sr.getByLabel(pPotLabel,mcpots);
      totpots->totpot     = mcpots->totpot;
      totpots->totgoodpot = mcpots->totgoodpot;
      totpots->totspills  = mcpots->totspills;
      totpots->goodspills = mcpots->goodspills;
     }
          
    std::unique_ptr< sumdata::POTSum > p(totpots);
    sr.put(std::move(p));
  }
  
  //......................................................................
  void MIN::reconfigure(const fhicl::ParameterSet& pset)
  {
    fMinPOTCut      = pset.get< double >("MinPOTCut"     );
    fMinHornICut    = pset.get< double >("MinHornICut"   );
    fMaxHornICut    = pset.get< double >("MaxHornICut"   );
    fMinPosXCut     = pset.get< double >("MinPosXCut"    );
    fMaxPosXCut     = pset.get< double >("MaxPosXCut"    );
    fMinPosYCut     = pset.get< double >("MinPosYCut"    );
    fMaxPosYCut     = pset.get< double >("MaxPosYCut"    );
    fMaxDeltaTCut   = pset.get< double >("MaxDeltaTCut"  );
    fBNBOffset      = pset.get< double >("BNBOffset"  );
  }

        


  

  //......................................................................
  void MIN::produce(art::Event& evt)
  {
       std::unique_ptr<sumdata::SpillData> spilldata(new sumdata::SpillData());
    
    //count the total number of spills
    //          totpots->totspills++;
    
    
    if(!evt.isRealData()){
      // remove empty spill info write from here, move to GENIEGen and fill
      isMC = true;
      
      // Ask MCTruth what I am

      // TODO: It would be ideal to get this from the Trigger information
      //       but currently it all reads back as TriggerType 0 ...
      art::Handle< std::vector<simb::MCTruth> > mcTruths;
      evt.getByLabel(pPotLabel, mcTruths);
      if(mcTruths.failedToGet() || mcTruths->empty()) return;
     
      const simb::MCTruth& mcTruth = (*mcTruths)[0];

      fMCOrigin = mcTruth.Origin();
      
      return;
    }
    
    //this will be the CORRECT UTC time we want ultimately
    struct timespec unixTime;
    
    // Get trigger information to extract trigger time
    auto trigv = evt.getValidHandle<std::vector<rawdata::RawTrigger> >(pRawDataLabel);
    if( trigv.failedToGet() )
      throw cet::exception("MIN") << "failed to get valid trigger handle";
    
    const rawdata::RawTrigger& trig = (*trigv)[0];
    
    if(trig.fTriggerMask_TriggerType != 1){
      // This isn't a BNB trigger. Don't try to get beam info
      return;
    }
    
    novadaq::timeutils::convertNovaTimeToUnixTime(trig.fTriggerTimingMarker_TimeStart, unixTime);
    
    unsigned long int uevt_sec = unixTime.tv_sec; //utval >> 32 & 0xFFFFFFFFF;
    unsigned long int uevt_nsec = unixTime.tv_nsec; //utval & 0xFFFFFFFFF;
    
    //DAQ event time 
    double DAQtime = uevt_sec + (0.000000001)*uevt_nsec;
   
            

   //torroid values and timestamps
    double tor860=0.;
    double tor875=0.;
    double tor860time=0.;
    
    //for horn current
    double thcurr=0.;
      
 
    try{
    
      bfp->GetNamedData((DAQtime)-fBNBOffset,"E:TOR860@,E:TOR875", &tor860, &tor860time, &tor875);
   
    // This checks that the event time is close to the database time 
   // to ensure that we are grabbing the right data
      double diff = 1e6; // a large number
      diff = (DAQtime-fBNBOffset) - tor860time;
      
      spilldata->deltaspilltimensec = diff*1e9;
      
      unsigned long int time_closest_int = (int) tor860time;
      double time_closest_ns = (tor860time - time_closest_int)*1e9;
      
      spilldata->spilltimesec = time_closest_int;
      spilldata->spilltimensec = time_closest_ns;
      
      // If spill meets time requirement, calculate all quantities and 
      // fill appropriate fields
      if(std::abs(spilldata->deltaspilltimensec) < fMaxDeltaTCut){
        
        double temppot = tor860;
        
        if(temppot < 0.02) temppot = tor875;
        
        // This guarantees that when toroid values go negative, total POT is not reduced
        if(temppot < 0.0) temppot = 0.0;
        
        spilldata->spillpot = temppot;
        
        // Total POT
        totpots->totpot += spilldata->spillpot;
        totpots->totspills += 1;
        
        //These were 0 in old code, so set them to 0 here.
        spilldata->gpsspilltimesec  = 0.;
        spilldata->gpsspilltimensec = 0.;
  
        // Horn current
        bfp->GetNamedData(DAQtime-fBNBOffset,"E:THCURR", &thcurr);

        double ihorn = 0.;
        
        // Unlike NuMI, BNB gets horn current from a single input
        ihorn += (thcurr - (+0.01));
               
        spilldata->hornI = ihorn;
        

      // Horn polarity
        spilldata->isRHC = false;
        if(ihorn > 0){
          spilldata->isRHC = true;
        }

        
//===========Beam position==================
double HPTG1, VPTG1, HP875, VP875;
double xp, yp;


      bfp->GetNamedData((DAQtime)-fBNBOffset,"E:HPTG1,E:VPTG1,E:HP875,E:VP875", &HPTG1, &VPTG1, &HP875, &VP875);

  
       double x2 = HPTG1-hptg1offset;
       double x1 = HP875-hp875offset;
       double dx = x2-x1;
       xp = ((x2+dx)/(hptg1station-hp875station))*(tgtstation-hptg1station); 
       double y2 = VPTG1-vptg1offset;
       double y1 = VP875-vp875offset;
       double dy= y2-y1;
       yp = ((y2+dy)/(vptg1station-vp875station))*(tgtstation-vptg1station); 
 


  
       spilldata->posx = xp; //make it in mm
       spilldata->posy = yp; //make it in mm




//==========================================

        spilldata->goodbeam = 0;
        
        // Apply good beam cuts
        if((spilldata->spillpot > fMinPOTCut)        // POT
           && (spilldata->hornI > fMinHornICut
               && spilldata->hornI < fMaxHornICut)   // Horn I
           && (spilldata->posx > fMinPosXCut
               && spilldata->posx < fMaxPosXCut)     // x position
           && (spilldata->posy > fMinPosYCut
               && spilldata->posy < fMaxPosYCut)     // y position
           && (std::abs(spilldata->deltaspilltimensec) < fMaxDeltaTCut)) // delta t
        {
          spilldata->goodbeam = 1;
        } // Good beam cuts
        
       
        // Total good POT
        if(spilldata->goodbeam == 1){
          totpots->totgoodpot += spilldata->spillpot;
          totpots->goodspills += 1;
        }
        
        
        /* if(!(spilldata->spillpot>fMinPOTCut)){
         mf::LogInfo("MIN") << "Failed POT cut" << std::endl;
         }
         if(!(spilldata->hornI>fMinHornICut && spilldata->hornI<fMaxHornICut)){
         mf::LogInfo("MIN") << "Failed horn I cut" << std::endl;
         }
         
         if(!(std::abs(spilldata->deltaspilltimensec)<1e9)){
         mf::LogInfo("MIN") << "Failed time cut" << std::endl;
         }
        */
        
      } // if deltaspilltimensec < fMaxDeltaTCut
      
      else mf::LogInfo("MIN") << "Spill did not meet MaxDeltaT cut\n";
      
    }
    catch (WebAPIException we){
      // Eventually put in better error catching. What if each device isn't filled?
      mf::LogError("MIN") << "Exception: " << we.what() << "\n";
      
      
      //There are at least two types of exceptions. One, relatively infrequent case, which is 
      //from database flakes. This should be aborted on and retried. However, much more frequently,
      //some of the variable information is missing from the record. This is unfortunate. 
      //We have implemented a nasty way to try to tell the difference between the two; it isn't 
      //super robust and there should be a better method. If you ever find one,
      //please implement it.

      std::string exceptionString(we.what());
      
      std::size_t missVar = exceptionString.find("variable not found");
      if(missVar != std::string::npos){
        mf::LogInfo("MIN") << "***************We have a case of missing data\n"<<
          "in the beam database! This isn't awesome, but \n"<<
          "is probably ok and we have decided to not abort.  \n"<<
          "This decision was made on string comparisions and \n"<<
          "is a bad hack.In the future, we should actually \n"<<
          "throw different types of exceptions for different \n"<<
          "types of problems!\n"<<
          "***************"<<std::endl;
      }
      else{
        std::abort();
      }
      
    } // end try/catch

    evt.put(std::move(spilldata));
    
  } //end event loop
  
  

  DEFINE_ART_MODULE(MIN)
  
} // end namespace ifdb
////////////////////////////////////////////////////////////////////////