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 #include <TH2.h>
0021
0022 #include <boost/format.hpp>
0023
0024 #include <cassert>
0025 #include <fstream>
0026 #include <memory>
0027
0028 //_________________________________________________________
0029 MicromegasCombinedDataCalibration::MicromegasCombinedDataCalibration( const std::string& name ):
0030   SubsysReco( name )
0031 {}
0032
0033 //_____________________________________________________________________
0034 int MicromegasCombinedDataCalibration::Init(PHCompositeNode* /*topNode*/ )
0035 {
0036
0037   std::cout << "MicromegasCombinedDataCalibration::Init - m_do_evaluation: " << m_do_evaluation << std::endl;
0038   std::cout << "MicromegasCombinedDataCalibration::Init - m_evaluation_filename: " << m_evaluation_filename << std::endl;
0039
0040   // histogram evaluation
0041   if( m_do_evaluation )
0042   { m_evaluation_file.reset( TFile::Open( m_evaluation_filename.c_str(), "RECREATE" ) ); }
0043
0044   return Fun4AllReturnCodes::EVENT_OK;
0045 }
0046
0047 //____________________________________________________________________________..
0048 int MicromegasCombinedDataCalibration::InitRun(PHCompositeNode* /*topNode*/)
0049 { return Fun4AllReturnCodes::EVENT_OK; }
0050
0051 //___________________________________________________________________________
0052 int MicromegasCombinedDataCalibration::process_event(PHCompositeNode *topNode)
0053 {
0054
0055   // load raw hits container
0056   auto rawhitcontainer = findNode::getClass<MicromegasRawHitContainer>(topNode, m_rawhitnodename);
0057   assert(rawhitcontainer);
0058
0059   // loop over TPOT packets
0060   for (unsigned int ihit = 0; ihit < rawhitcontainer->get_nhits(); ++ihit)
0061   {
0062     const auto rawhit = rawhitcontainer->get_hit(ihit);
0063     const auto fee_id = rawhit->get_fee();
0064     const auto channel = rawhit->get_channel();
0065     const auto sample_range = std::make_pair(rawhit->get_sample_begin(), rawhit->get_sample_end());
0066     if( Verbosity() )
0067     {
0068       std::cout
0069         << "MicromegasCombinedDataCalibration::process_event -"
0070         << " waveform: " << ihit
0071         << " fee_id: " << fee_id
0072         << " channel: " << channel
0073         << " samples: (" << sample_range.first << "," << sample_range.second << ")"
0074         << std::endl;
0075     }
0076
0077     // find relevant profile histogram
0078     TProfile* profile = nullptr;
0079     auto piter = m_profile_map.lower_bound( fee_id );
0080     if( piter == m_profile_map.end() || fee_id < piter->first )
0081     {
0082       // create and insert
0083       const auto pname = (boost::format("p_adc_channel_%i") % fee_id ).str();
0084       profile = new TProfile( pname.c_str(), "ADC vs channel;channel;adc", MicromegasDefs::m_nchannels_fee, 0, MicromegasDefs::m_nchannels_fee );
0085       profile->SetErrorOption( "s" );
0086       m_profile_map.insert(  piter, std::make_pair( fee_id, profile ) );
0087     } else {
0088       profile = piter->second;
0089     }
0090
0091     // find relevant 2D histogram
0092     TH2* histogram = nullptr;
0093     if( m_do_evaluation )
0094     {
0095       auto hiter = m_histogram_map.lower_bound( fee_id );
0096       if( hiter == m_histogram_map.end() || fee_id < hiter->first )
0097       {
0098         static constexpr int max_adc = 1100;
0099
0100         // create and insert
0101         const auto hname = (boost::format("h_adc_channel_%i") % fee_id ).str();
0102         histogram = new TH2F( hname.c_str(), "ADC vs channel;channel;adc",
0103           MicromegasDefs::m_nchannels_fee, 0, MicromegasDefs::m_nchannels_fee,
0104           max_adc, 0, max_adc );
0105         m_histogram_map.insert(  hiter, std::make_pair( fee_id, histogram ) );
0106       } else {
0107         histogram = hiter->second;
0108       }
0109     }
0110
0111     // fill
0112     for( auto is = std::max(m_sample_min,sample_range.first); is < std::min(m_sample_max,sample_range.second); ++ is )
0113     {
0114         const uint16_t adc =  rawhit->get_adc(is);
0115         if( adc != MicromegasDefs::m_adc_invalid )
0116         {
0117           if( profile ) { profile->Fill( channel, adc); }
0118           if( histogram ) { histogram->Fill( channel, adc); }
0119         }
0120     }
0121   }
0122   return Fun4AllReturnCodes::EVENT_OK;
0123 }
0124
0125 //_____________________________________________________________________
0126 int MicromegasCombinedDataCalibration::End(PHCompositeNode* /*topNode*/ )
0127 {
0128
0129   // write calibration data to ouput file
0130   if( m_profile_map.empty() )
0131   {
0132     std::cout << "MicromegasCombinedDataCalibration::End - no data" << std::endl;
0133   } else {
0134
0135     // create calibration data object
0136     MicromegasCalibrationData calibration_data;
0137     for( const auto& [fee_id, profile]:m_profile_map )
0138     {
0139       for( int i = 0; i < profile->GetNbinsX(); ++ i )
0140       {
0141         const auto pedestal = profile->GetBinContent(i+1);
0142         const auto rms = profile->GetBinError(i+1);
0143         calibration_data.set_pedestal( fee_id, i, pedestal );
0144         calibration_data.set_rms( fee_id, i, rms );
0145       }
0146     }
0147     calibration_data.write( m_calibration_filename );
0148   }
0149
0150   if( m_do_evaluation && m_evaluation_file )
0151   {
0152     m_evaluation_file->cd();
0153     for( const auto& [feeid,h]:m_histogram_map ) { h->Write(); }
0154     for( const auto& [feeid,p]:m_profile_map ) { p->Write(); }
0155   }
0156
0157   return Fun4AllReturnCodes::EVENT_OK;
0158 }