Inheritance diagram for opdet::LEDCalibrationAna:

Public Member Functions
	LEDCalibrationAna (const fhicl::ParameterSet &)

void	endJob ()

void	analyze (const art::Event &)

Private Member Functions
uint32_t	ShaperToChannel (uint32_t Shaper)

Private Attributes
std::string	fInputModule

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

bool	fMakeNonCoincTree

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

Detailed Description

Definition at line 42 of file LEDCalibrationAna_module.cc.

Constructor & Destructor Documentation

opdet::LEDCalibrationAna::LEDCalibrationAna ( const fhicl::ParameterSet & pset )

Definition at line 95 of file LEDCalibrationAna_module.cc.

     : 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 const> 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");
   }

Member Function Documentation

void opdet::LEDCalibrationAna::analyze ( const art::Event & evt )

Definition at line 215 of file LEDCalibrationAna_module.cc.

   {
     auto const clock_data =
       art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
 
     fRunID = evt.run();
     fEventID = evt.event();
 
     // 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 = clock_data.OpticalClock().Frame(TimeStamp);
       uint32_t TimeSlice = clock_data.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 = clock_data.OpticalClock().Frame(TimeStamp);
       uint32_t TimeSlice = clock_data.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();
               }
             }
           }
         }
       }
     }
   }

void opdet::LEDCalibrationAna::endJob ( )

Definition at line 145 of file LEDCalibrationAna_module.cc.

   {
     art::ServiceHandle<art::TFileService const> 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;
     }
   }

uint32_t opdet::LEDCalibrationAna::ShaperToChannel ( uint32_t Shaper )

private

Definition at line 303 of file LEDCalibrationAna_module.cc.

   {
     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];
   }