packages/micromegas/MicromegasCombinedDataCalibration.cc

0001 /*!
0002  * \file MicromegasCombinedDataCalibration.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "MicromegasCombinedDataCalibration.h"
0007 #include "MicromegasCalibrationData.h"
0008 #include "MicromegasDefs.h"
0009
0010 #include <ffarawobjects/MicromegasRawHit.h>
0011 #include <ffarawobjects/MicromegasRawHitContainer.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <phool/getClass.h>
0016 #include <phool/PHCompositeNode.h>
0017
0018 #include <TFile.h>
0019 #include <TProfile.h>
0020
0021 #include <boost/format.hpp>
0022
0023 #include <cassert>
0024 #include <fstream>
0025 #include <memory>
0026
0027 //_________________________________________________________
0028 MicromegasCombinedDataCalibration::MicromegasCombinedDataCalibration( const std::string& name ):
0029   SubsysReco( name )
0030 {}
0031
0032 //_____________________________________________________________________
0033 int MicromegasCombinedDataCalibration::Init(PHCompositeNode* /*topNode*/ )
0034 {
0035   // histogram evaluation
0036   return Fun4AllReturnCodes::EVENT_OK;
0037 }
0038
0039 //____________________________________________________________________________..
0040 int MicromegasCombinedDataCalibration::InitRun(PHCompositeNode* /*topNode*/)
0041 { return Fun4AllReturnCodes::EVENT_OK; }
0042
0043 //___________________________________________________________________________
0044 int MicromegasCombinedDataCalibration::process_event(PHCompositeNode *topNode)
0045 {
0046
0047   // load raw hits container
0048   auto rawhitcontainer = findNode::getClass<MicromegasRawHitContainer>(topNode, m_rawhitnodename);
0049   assert(rawhitcontainer);
0050
0051   // loop over TPOT packets
0052   for (unsigned int ihit = 0; ihit < rawhitcontainer->get_nhits(); ++ihit)
0053   {
0054     const auto rawhit = rawhitcontainer->get_hit(ihit);
0055     const auto fee_id = rawhit->get_fee();
0056     const auto channel = rawhit->get_channel();
0057     const auto sample_range = std::make_pair(rawhit->get_sample_begin(), rawhit->get_sample_end());
0058     if( Verbosity() )
0059     {
0060       std::cout
0061         << "MicromegasCombinedDataCalibration::process_event -"
0062         << " waveform: " << ihit
0063         << " fee_id: " << fee_id
0064         << " channel: " << channel
0065         << " samples: (" << sample_range.first << "," << sample_range.second << ")"
0066         << std::endl;
0067     }
0068
0069     // find relevant profile histogram
0070     TProfile* profile = nullptr;
0071     auto piter = m_profile_map.lower_bound( fee_id );
0072     if( piter == m_profile_map.end() || fee_id < piter->first )
0073     {
0074       // create and insert
0075       const auto hname = (boost::format("h_adc_channel_%i") % fee_id ).str();
0076       profile = new TProfile( hname.c_str(), "ADC vs channel;channel;adc", MicromegasDefs::m_nchannels_fee, 0, MicromegasDefs::m_nchannels_fee );
0077       profile->SetErrorOption( "s" );
0078       m_profile_map.insert(  piter, std::make_pair( fee_id, profile ) );
0079     } else {
0080       profile = piter->second;
0081     }
0082
0083     // fill
0084     for( auto is = std::max(m_sample_min,sample_range.first); is < std::min(m_sample_max,sample_range.second); ++ is )
0085     {
0086         const uint16_t adc =  rawhit->get_adc(is);
0087         if( adc != MicromegasDefs::m_adc_invalid )
0088         { profile->Fill( channel, adc); }
0089     }
0090   }
0091   return Fun4AllReturnCodes::EVENT_OK;
0092 }
0093
0094 //_____________________________________________________________________
0095 int MicromegasCombinedDataCalibration::End(PHCompositeNode* /*topNode*/ )
0096 {
0097
0098   // write calibration data to ouput file
0099   if( m_profile_map.empty() )
0100   {
0101     std::cout << "MicromegasCombinedDataCalibration::End - no data" << std::endl;
0102   } else {
0103
0104     // create calibration data object
0105     MicromegasCalibrationData calibration_data;
0106     for( const auto& [fee_id, profile]:m_profile_map )
0107     {
0108       for( int i = 0; i < profile->GetNbinsX(); ++ i )
0109       {
0110         const auto pedestal = profile->GetBinContent(i+1);
0111         const auto rms = profile->GetBinError(i+1);
0112         calibration_data.set_pedestal( fee_id, i, pedestal );
0113         calibration_data.set_rms( fee_id, i, rms );
0114       }
0115     }
0116     calibration_data.write( m_calibration_filename );
0117   }
0118
0119   return Fun4AllReturnCodes::EVENT_OK;
0120 }