LEDCalibrationAna_module.cc

Go to the documentation of this file.

// -*- mode: c++; c-basic-offset: 2; -*-
// Ben Jones, MIT, 2013
//
//  This ana module extracts pedestals and gains from 
//  LED calibration run data (incomplete)
//


#ifndef LEDCalibrationAna_H
#define LEDCalibrationAna_H 1

// LArSoft includes
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RawData/OpDetWaveform.h"
#include "larana/OpticalDetector/OpHitFinder/AlgoThreshold.h"
#include "larana/OpticalDetector/OpHitFinder/PedAlgoEdges.h"
#include "larana/OpticalDetector/OpHitFinder/PulseRecoManager.h"
#include "lardata/DetectorInfoServices/DetectorClocksService.h"

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// ROOT includes
#include "TH1.h"
#include "THStack.h"
#include "TF1.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include <TTree.h>

// C++ Includes
#include <map>
#include <vector>
#include <iostream>
#include <cstring>
#include <sstream>
#include "math.h"
#include <climits>

namespace opdet {
 
  class LEDCalibrationAna : public art::EDAnalyzer{
  public:
 
    // Standard constructor and destructor for an ART module.
    LEDCalibrationAna(const fhicl::ParameterSet&);
    virtual ~LEDCalibrationAna();

    void beginJob();

    void endJob();

    // The analyzer routine, called once per event. 
    void analyze (const art::Event&); 

  private:

    // The parameters we'll read from the .fcl file.
    std::string fInputModule;              // Input tag for OpDet collection
    uint32_t    fTriggerChannel;
    float       fCoincThreshold;
    float       fMaxTimeMean;
    float       fMaxTimeThresh;
    float       fAreaDivs;
    float       fAreaMin;
    float       fAreaMax;
    float       fTriggerDelay;
    
    pmtana::PulseRecoManager  fPulseRecoMgr;
    pmtana::AlgoThreshold     fThreshAlg;
    pmtana::PedAlgoEdges      fPedAlg;
    TTree * fPulseTree;
    TTree * fPulseTreeNonCoinc;
    Float_t fPeak;
    Float_t fPedRMS;
    Float_t fPedMean;
    Float_t fArea;
    Float_t fTBegin;
    Float_t fTEnd;
    Float_t fTMax;
    Float_t fOffset;
    Int_t   fEventID;
    Int_t   fRunID;
    Int_t   fChannel;
    Int_t   fShaper;

    uint32_t ShaperToChannel(uint32_t Shaper);
   
    bool    fMakeNonCoincTree;

    std::map<uint32_t, std::vector<double> > fAreas;

  
    
  };


} 

#endif 

namespace opdet {

  //-----------------------------------------------------------------------
  // Constructor
  LEDCalibrationAna::LEDCalibrationAna(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset),
      fPulseRecoMgr(),
      fThreshAlg(),
      fPedAlg()
  {
    // Indicate that the Input Module comes from .fcl
    fInputModule    = pset.get<std::string>("InputModule");
    fTriggerChannel = pset.get<uint32_t>   ("TriggerChannel");
    fTriggerDelay   = pset.get<uint32_t>   ("TriggerDelay");
    fCoincThreshold = pset.get<float>      ("CoincThreshold");
    fMaxTimeThresh  = pset.get<float>      ("MaxTimeThresh");
    fMaxTimeMean    = pset.get<float>      ("MaxTimeMean");
    fAreaMin        = pset.get<float>      ("AreaMin");
    fAreaMax        = pset.get<float>      ("AreaMax");
    fAreaDivs       = pset.get<float>      ("AreaDivs");
    fMakeNonCoincTree=pset.get<bool>       ("MakeNonCoincTree");
    
    fPulseRecoMgr.AddRecoAlgo(&fThreshAlg);
    fPulseRecoMgr.SetDefaultPedAlgo(&fPedAlg);

    art::ServiceHandle<art::TFileService> tfs;

    fPulseTree = tfs->make<TTree>("fPulseTree","fPulseTree");
    fPulseTree->Branch("Area",     &fArea,     "Area/F");
    fPulseTree->Branch("Peak",     &fPeak,     "Peak/F");
    fPulseTree->Branch("TBegin",   &fTBegin,   "TBegin/F");
    fPulseTree->Branch("TEnd",     &fTEnd,     "TEnd/F");
    fPulseTree->Branch("TMax",     &fTMax,     "TMax/F");
    fPulseTree->Branch("Channel",  &fChannel,  "Channel/I");
    fPulseTree->Branch("Shaper",   &fShaper,   "Shaper/I");
    fPulseTree->Branch("PedMean",  &fPedMean,  "PedMean/F");
    fPulseTree->Branch("PedRMS",   &fPedRMS,   "PedRMS/F");
    fPulseTree->Branch("Offset",   &fOffset,   "Offset/F");
    fPulseTree->Branch("EventID",  &fEventID,  "EventID/I");
    fPulseTree->Branch("RunID",    &fRunID,    "RunID/I");

    fPulseTreeNonCoinc = tfs->make<TTree>("fPulseTreeNonCoinc","fPulseTreeNonCoinc");
    fPulseTreeNonCoinc->Branch("Area",     &fArea,     "Area/F");
    fPulseTreeNonCoinc->Branch("Peak",     &fPeak,     "Peak/F");
    fPulseTreeNonCoinc->Branch("TBegin",   &fTBegin,   "TBegin/F");
    fPulseTreeNonCoinc->Branch("TEnd",     &fTEnd,     "TEnd/F");
    fPulseTreeNonCoinc->Branch("TMax",     &fTMax,     "TMax/F");
    fPulseTreeNonCoinc->Branch("Channel",  &fChannel,  "Channel/I");
    fPulseTreeNonCoinc->Branch("Shaper",   &fShaper,   "Shaper/I");
    fPulseTreeNonCoinc->Branch("PedMean",  &fPedMean,  "PedMean/F");
    fPulseTreeNonCoinc->Branch("PedRMS",   &fPedRMS,   "PedRMS/F");
    fPulseTreeNonCoinc->Branch("EventID",  &fEventID,  "EventID/I");
    fPulseTreeNonCoinc->Branch("RunID",    &fRunID,    "RunID/I");

  }

  //-----------------------------------------------------------------------
  // Destructor
  LEDCalibrationAna::~LEDCalibrationAna() 
  {}
   
  //-----------------------------------------------------------------------
  void LEDCalibrationAna::beginJob()
  {
  }

  //-----------------------------------------------------------------------
  void LEDCalibrationAna::endJob()
  { 
    art::ServiceHandle<art::TFileService> tfs;

    for(auto it = fAreas.begin(); it!=fAreas.end(); ++it)
      {
        uint32_t Channel = it->first;
        
        std::stringstream histname;
        histname.flush();
        histname<<"ch"<<Channel<<"area";

        TH1D * HistArea = tfs->make<TH1D>(histname.str().c_str(), histname.str().c_str(), fAreaDivs, fAreaMin, fAreaMax);

        for(size_t j=0; j!=it->second.size(); ++j)
          {
            HistArea->Fill(it->second.at(j));
          }
        
        std::stringstream fitname;
        fitname.flush();
        fitname<<"ch"<<Channel<<"fit";

        double Max = HistArea->GetMaximum();
        double Mid = HistArea->GetBinContent(fAreaDivs/2.);
        
        TF1 * GausFit = new TF1(fitname.str().c_str(),
                            "gaus(0)+gaus(3)+gaus(6)",
                            fAreaMin,
                            fAreaMax);
        
        GausFit->SetParameters(Mid, (fAreaMin+fAreaMax)/2., (fAreaMax-fAreaMin)/2.,
                               Max, 0, (fAreaMin+fAreaMax)/8., 
                               Max/5., 0, (fAreaMin+fAreaMax)/4.);


        GausFit->SetParLimits(0,0,1.1*Max);
        GausFit->SetParLimits(1,0,fAreaMax);
        GausFit->SetParLimits(2,0,fAreaMax);

        GausFit->SetParLimits(3,0,1.1*Max);
        GausFit->FixParameter(4,0);
        GausFit->SetParLimits(5,0,(fAreaMin+fAreaMax)/2.);

        GausFit->SetParLimits(6,0,1.1*Max);
        GausFit->FixParameter(7,0);
        GausFit->SetParLimits(8,0,(fAreaMin+fAreaMax)/2.);


        HistArea->Fit(GausFit);
        
        double Mean     = GausFit->GetParameter(1);
        double Width    = GausFit->GetParameter(2);
        
        double MeanErr  = GausFit->GetParError(1);
        double WidthErr = GausFit->GetParError(2);
        
        double NPE      = pow(Mean,2)/pow(Width,2);
        double SPEScale = Mean / NPE;

        double NPEError = NPE      * pow(2.*(pow(MeanErr/Mean,2) + pow(WidthErr/Width,2)),0.5);
        double SPEError = SPEScale * pow(2.*pow(WidthErr/Width,2) + pow(MeanErr/Mean,2),0.5);
        
        std::cout <<"Channel " << Channel<< ":\tSPE Scale \t"<<SPEScale<<"\t +/- \t" << SPEError<<",\t NPE \t" << NPE<<"\t +/- \t"<< NPEError <<std::endl;
      }
    
  }



  //-----------------------------------------------------------------------
  void LEDCalibrationAna::analyze(const art::Event& evt) 
  {

    auto const* ts = lar::providerFrom<detinfo::DetectorClocksService>();


    fRunID=evt.run();
    fEventID=evt.event();

    art::ServiceHandle<art::TFileService> tfs;

    // Create a handle for our vector of pulses
    art::Handle< std::vector< raw::OpDetWaveform > > OpDetWaveformHandle;

    // Read in WaveformHandle
    evt.getByLabel(fInputModule, OpDetWaveformHandle);

    std::map<uint32_t, std::vector<int> > OrgOpDigitByChannel;
    
    for(size_t i=0; i!=OpDetWaveformHandle->size(); ++i)
      {
        OrgOpDigitByChannel[ShaperToChannel(OpDetWaveformHandle->at(i).ChannelNumber())].push_back(i);
      }

    std::vector<uint32_t> FrameNumbersForTrig;
    std::vector<uint32_t> TimeSlicesForTrig;
    
    for(size_t i=0; i!=OrgOpDigitByChannel[fTriggerChannel].size(); ++i)
      {
        double TimeStamp = OpDetWaveformHandle->at(OrgOpDigitByChannel[ fTriggerChannel][i] ).TimeStamp();
        uint32_t Frame     = ts->OpticalClock().Frame(TimeStamp);
        uint32_t TimeSlice = ts->OpticalClock().Sample(TimeStamp);
        FrameNumbersForTrig.push_back(Frame);
        TimeSlicesForTrig.push_back(TimeSlice);
      }

    for(size_t i=0; i!=OpDetWaveformHandle->size(); ++i)
      {
        double TimeStamp = OpDetWaveformHandle->at(i).TimeStamp();
        uint32_t Frame     = ts->OpticalClock().Frame(TimeStamp);
        uint32_t TimeSlice = ts->OpticalClock().Sample(TimeStamp);
        fShaper   = OpDetWaveformHandle->at(i).ChannelNumber();
        fChannel  = ShaperToChannel(fShaper);
        

        if(uint32_t(fChannel) != fTriggerChannel)
          {
            for(size_t j=0; j!=FrameNumbersForTrig.size(); ++j)
              {
                if( (Frame == FrameNumbersForTrig.at(j))
                    && (fabs(TimeSlice - TimeSlicesForTrig.at(j) - fTriggerDelay)<fCoincThreshold))
                  {
                    
                   
                    const raw::OpDetWaveform& wf = OpDetWaveformHandle->at(i);

                    //fPulseRecoMgr.RecoPulse(wf);
                    fPulseRecoMgr.Reconstruct(wf);
                    
                    size_t NPulses = fThreshAlg.GetNPulse();

                    fOffset  = TimeSlice-TimeSlicesForTrig.at(j);
                    //fPedMean = fThreshAlg.PedMean();
                    //fPedRMS  = fThreshAlg.PedRms();
                    
                    for(size_t k=0; k!=NPulses; ++k)
                      {
                        if(fabs(fMaxTimeMean-fThreshAlg.GetPulse(k).t_max)<fMaxTimeThresh)
                          {
                            
                            fPeak    = fThreshAlg.GetPulse(k).peak;
                            fArea    = fThreshAlg.GetPulse(k).area;
                            fTBegin  = fThreshAlg.GetPulse(k).t_start;
                            fTEnd    = fThreshAlg.GetPulse(k).t_end;
                            fTMax    = fThreshAlg.GetPulse(k).t_max;
                            fPedMean = fThreshAlg.GetPulse(k).ped_mean;
                            fPedRMS  = fThreshAlg.GetPulse(k).ped_sigma;

                            fPulseTree->Fill();
                            
                            fAreas[fChannel].push_back(fArea);
                          }
                        else if(fMakeNonCoincTree)
                          {
                            fPeak    = fThreshAlg.GetPulse(k).peak;
                            fArea    = fThreshAlg.GetPulse(k).area;
                            fTBegin  = fThreshAlg.GetPulse(k).t_start;
                            fTEnd    = fThreshAlg.GetPulse(k).t_end;
                            fTMax    = fThreshAlg.GetPulse(k).t_max;
                            
                            fPulseTreeNonCoinc->Fill();
                          }
                      }
                  }             
              }
          }
      }

  }

  //---------------------------------
  uint32_t LEDCalibrationAna::ShaperToChannel(uint32_t Shaper)
  {
    static std::map<uint32_t, uint32_t> ShaperToChannelMap;
    if(ShaperToChannelMap.size()==0)
      {
        
        // temporary
        for(size_t i=0; i!=40; ++i)
          {
            ShaperToChannelMap[i]=i;
          }
      }
    
    return ShaperToChannelMap[Shaper];
  }


} // namespace opdet

namespace opdet {
  DEFINE_ART_MODULE(LEDCalibrationAna)
}