subsystems/tpot/TpotMon.cc

0001 #include "TpotMon.h"
0002 #include "TpotMonDefs.h"
0003
0004 #include <onlmon/OnlMon.h>  // for OnlMon
0005 #include <onlmon/OnlMonServer.h>
0006
0007 #include <Event/msg_profile.h>
0008 #include <Event/Event.h>
0009
0010 #include <micromegas/MicromegasDefs.h>
0011
0012 #include <TH1.h>
0013 #include <TH2.h>
0014 #include <TH2Poly.h>
0015
0016 #include <cmath>
0017 #include <cstdio>  // for printf
0018 #include <cstdlib>
0019 #include <fstream>
0020 #include <iostream>
0021 #include <memory>
0022 #include <sstream>
0023 #include <string>  // for allocator, string, char_traits
0024
0025 namespace
0026 {
0027
0028   enum
0029   {
0030     TRGMESSAGE = 1,
0031     FILLMESSAGE = 2
0032   };
0033
0034   // get first member of pairs into a list
0035   std::vector<double> get_x( const MicromegasGeometry::point_list_t& point_list )
0036   {
0037     std::vector<double> out;
0038     std::transform( point_list.begin(), point_list.end(), std::back_inserter( out ), []( const auto& p ) { return p.first; } );
0039     return out;
0040   }
0041
0042   // get second member of pairs into a list
0043   std::vector<double> get_y( const MicromegasGeometry::point_list_t& point_list )
0044   {
0045     std::vector<double> out;
0046     std::transform( point_list.begin(), point_list.end(), std::back_inserter( out ), []( const auto& p ) { return p.second; } );
0047     return out;
0048   }
0049
0050   // streamer for sample window
0051   std::ostream& operator << ( std::ostream&o, const TpotMon::sample_window_t& window )
0052   {
0053     o << "{ " << window.first << ", " << window.second << "}";
0054     return o;
0055   }
0056
0057   // number of sampa chips per fee
0058   static constexpr int m_nsampa_fee = 8;
0059
0060   /* see: https://git.racf.bnl.gov/gitea/Instrumentation/sampa_data/src/branch/fmtv2/README.md */
0061   enum ModeBitType
0062   {
0063     BX_COUNTER_SYNC_T = 0,
0064     ELINK_HEARTBEAT_T = 1,
0065     SAMPA_EVENT_TRIGGER_T = 2,
0066     CLEAR_LV1_LAST_T = 6,
0067     CLEAR_LV1_ENDAT_T = 7
0068   };
0069
0070   /*
0071    * returns true if a given channel for a given FEE is permanently masked
0072    * for now all channels from 128 to 255, for FEE 8 (SCOZ) are masked
0073    */
0074   bool is_masked( int fee_id, int channel )
0075   { return fee_id==8 && channel>=128; }
0076
0077 }
0078
0079 //__________________________________________________
0080 TpotMon::TpotMon(const std::string &name)
0081   : OnlMon(name)
0082 {
0083   // setup default calibration filename
0084   // note: this can be overriden by calling set_calibration_filename from the parent macro
0085   const auto tpotcalib = getenv("TPOTCALIB");
0086   if (!tpotcalib)
0087   {
0088     std::cout << "TpotMon::TpotMon - TPOTCALIB environment variable not set" << std::endl;
0089     exit(1);
0090   }
0091
0092   m_calibration_filename = std::string(tpotcalib) + "/" + "TPOT_Pedestal-000.root";
0093 }
0094
0095 //__________________________________________________
0096 int TpotMon::Init()
0097 {
0098
0099   if( Verbosity() )
0100   {
0101     std::cout << "TpotMon::Init - m_calibration_filename: " << m_calibration_filename << std::endl;
0102     std::cout << "TpotMon::Init - m_max_sample: " << m_max_sample << std::endl;
0103     std::cout << "TpotMon::Init - m_sample_window_signal: " << m_sample_window_signal << std::endl;
0104     std::cout << "TpotMon::Init - m_n_sigma: " << m_n_sigma << std::endl;
0105   }
0106
0107   // setup calibrations
0108   if( std::ifstream( m_calibration_filename.c_str() ).good() )
0109   {
0110     m_calibration_data.read( m_calibration_filename );
0111   } else {
0112     std::cout << "TpotMon::Init -"
0113       << " file " << m_calibration_filename << " cannot be opened."
0114       << " No calibration loaded"
0115       << std::endl;
0116   }
0117
0118   // server instance
0119   auto se = OnlMonServer::instance();
0120
0121   // map tile centers to fee id
0122   const auto fee_id_list = m_mapping.get_fee_id_list();
0123   for( const auto& fee_id:fee_id_list )
0124   {
0125     const auto tile_id = MicromegasDefs::getTileId( m_mapping.get_hitsetkey(fee_id));
0126     m_tile_centers.emplace(fee_id, m_geometry.get_tile_center(tile_id));
0127   }
0128
0129   // counters
0130   /* arbitrary counters. First bin is number of events */
0131   m_counters = new TH1F( "m_counters", "counters", 10, 0, 10 );
0132   m_counters->GetXaxis()->SetBinLabel(TpotMonDefs::kEventCounter, "RCDAQ frame" );
0133   m_counters->GetXaxis()->SetBinLabel(TpotMonDefs::kValidEventCounter, "valid RCDAQ frames" );
0134   m_counters->GetXaxis()->SetBinLabel(TpotMonDefs::kTriggerCounter, "triggers" );
0135   m_counters->GetXaxis()->SetBinLabel(TpotMonDefs::kHeartBeatCounter, "heartbeats" );
0136   se->registerHisto(this, m_counters);
0137
0138   // global occupancy
0139   m_detector_multiplicity_phi = new TH2Poly( "m_detector_multiplicity_phi", "multiplicity (#phi); z (cm); x (cm)", -120, 120, -60, 60 );
0140   m_detector_occupancy_phi = new TH2Poly( "m_detector_occupancy_phi", "occupancy (#phi); z (cm); x (cm);occupancy (%)", -120, 120, -60, 60 );
0141   se->registerHisto(this, m_detector_occupancy_phi);
0142
0143   m_detector_multiplicity_z = new TH2Poly( "m_detector_multiplicity_z", "multiplicity (z); z (cm); x (cm)", -120, 120, -60, 60 );
0144   m_detector_occupancy_z = new TH2Poly( "m_detector_occupancy_z", "occupancy (z); z (cm); x(cm);occupancy (%)", -120, 120, -60, 60 );
0145   se->registerHisto(this, m_detector_occupancy_z);
0146
0147   // setup bins
0148   for( auto&& h:{m_detector_multiplicity_phi, m_detector_occupancy_phi, m_detector_multiplicity_z, m_detector_occupancy_z } )
0149   {
0150     setup_detector_bins( h );
0151     h->SetMinimum(0);
0152   }
0153
0154   // resist region occupancy
0155   m_resist_multiplicity_phi = new TH2Poly( "m_resist_multiplicity_phi", "multiplicity (#phi); z (cm); x (cm)", -120, 120, -60, 60 );
0156   m_resist_occupancy_phi = new TH2Poly( "m_resist_occupancy_phi", "occupancy (#phi); z (cm); x (cm);occupancy (%)", -120, 120, -60, 60 );
0157   se->registerHisto(this, m_resist_occupancy_phi);
0158
0159   m_resist_multiplicity_z = new TH2Poly( "m_resist_multiplicity_z", "multiplicity (z); z (cm); x (cm)", -120, 120, -60, 60 );
0160   m_resist_occupancy_z = new TH2Poly( "m_resist_occupancy_z", "occupancy (z); z (cm); x(cm);occupancy (%)", -120, 120, -60, 60 );
0161   se->registerHisto(this, m_resist_occupancy_z);
0162
0163   // setup bins
0164   for( auto&& h:{m_resist_multiplicity_z, m_resist_occupancy_z } )
0165   {
0166     setup_resist_bins( h, MicromegasDefs::SegmentationType::SEGMENTATION_Z );
0167     h->SetMinimum(0);
0168   }
0169
0170   for( auto&& h:{m_resist_multiplicity_phi, m_resist_occupancy_phi } )
0171   {
0172     setup_resist_bins( h, MicromegasDefs::SegmentationType::SEGMENTATION_PHI );
0173     h->SetMinimum(0);
0174   }
0175
0176   for( const auto& fee_id:fee_id_list )
0177   {
0178     // local copy of detector name
0179     const auto& detector_name=m_mapping.get_detname_sphenix(fee_id);
0180
0181     detector_histograms_t detector_histograms;
0182
0183     detector_histograms.m_counts_sample = new TH1I(
0184       Form( "m_counts_sample_%s", detector_name.c_str() ),
0185       Form( "hit count vs sample (%s);sample;counts", detector_name.c_str() ),
0186       m_max_sample, 0, m_max_sample );
0187     se->registerHisto(this, detector_histograms.m_counts_sample);
0188
0189     static constexpr int max_adc = 1100;
0190     detector_histograms.m_adc_sample = new TH2I(
0191       Form( "m_adc_sample_%s", detector_name.c_str() ),
0192       Form( "adc count vs sample (%s);sample;adc", detector_name.c_str() ),
0193       m_max_sample, 0, m_max_sample,
0194       max_adc, 0, max_adc );
0195     se->registerHisto(this, detector_histograms.m_adc_sample);
0196
0197     detector_histograms.m_adc_channel = new TH2I(
0198       Form( "m_adc_channel_%s", detector_name.c_str() ),
0199       Form( "adc count vs strip (%s);strip;adc", detector_name.c_str() ),
0200       MicromegasDefs::m_nchannels_fee, 0, MicromegasDefs::m_nchannels_fee,
0201       max_adc, 0, max_adc );
0202     se->registerHisto(this, detector_histograms.m_adc_channel);
0203
0204
0205     static constexpr int max_sample = 25;
0206     detector_histograms.m_sample_channel = new TH2I(
0207       Form("m_sample_channel_%s", detector_name.c_str() ),
0208       Form("strip vs sample (%s);strip;sample", detector_name.c_str()),
0209       MicromegasDefs::m_nchannels_fee, 0, MicromegasDefs::m_nchannels_fee,
0210       max_sample, 0 , max_sample );
0211     se->registerHisto(this, detector_histograms.m_sample_channel);
0212
0213     // hit charge
0214     static constexpr double max_hit_charge = 1100;
0215     detector_histograms.m_hit_charge = new TH1I(
0216       Form( "m_hit_charge_%s", detector_name.c_str() ),
0217       Form( "hit charge distribution (%s);adc", detector_name.c_str() ),
0218       100, 0, max_hit_charge );
0219     se->registerHisto(this, detector_histograms.m_hit_charge);
0220
0221     // hit multiplicity
0222     detector_histograms.m_hit_multiplicity = new TH1I(
0223       Form( "m_hit_multiplicity_%s", detector_name.c_str() ),
0224       Form( "hit multiplicity (%s);#hits", detector_name.c_str() ),
0225       MicromegasDefs::m_nchannels_fee+5, 0, MicromegasDefs::m_nchannels_fee+5 );
0226     se->registerHisto(this, detector_histograms.m_hit_multiplicity);
0227
0228     // waveform per channel
0229     detector_histograms.m_wf_vs_channel = new TH1F(
0230       Form( "m_wf_vs_channel_%s", detector_name.c_str() ),
0231       Form( "waveform profile (%s);strip", detector_name.c_str() ),
0232       MicromegasDefs::m_nchannels_fee, 0, MicromegasDefs::m_nchannels_fee );
0233     se->registerHisto(this, detector_histograms.m_wf_vs_channel);
0234
0235     // hit per channel
0236     detector_histograms.m_hit_vs_channel = new TH1F(
0237       Form( "m_hit_vs_channel_%s", detector_name.c_str() ),
0238       Form( "hit profile (%s);strip", detector_name.c_str() ),
0239       MicromegasDefs::m_nchannels_fee, 0, MicromegasDefs::m_nchannels_fee );
0240     se->registerHisto(this, detector_histograms.m_hit_vs_channel);
0241
0242     // heartbeat hit per channel
0243     detector_histograms.m_heartbeat_vs_channel = new TH1F(
0244       Form( "m_heartbeat_vs_channel_%s", detector_name.c_str() ),
0245       Form( "heartbeat profile (%s);channel", detector_name.c_str() ),
0246       m_nsampa_fee, 0, MicromegasDefs::m_nchannels_fee );
0247     se->registerHisto(this, detector_histograms.m_heartbeat_vs_channel);
0248
0249     // store in map
0250     m_detector_histograms.emplace( fee_id, std::move( detector_histograms ) );
0251
0252   }
0253
0254   // use monitor name for db table name
0255   Reset();
0256   return 0;
0257 }
0258
0259 //________________________________
0260 int TpotMon::BeginRun(const int /* runno */)
0261 {
0262   // if you need to read calibrations on a run by run basis
0263   // this is the place to do it
0264   return 0;
0265 }
0266
0267 //________________________________
0268 int TpotMon::process_event(Event* event)
0269 {
0270
0271   // increment a given bin number by weight
0272   auto increment = []( TH1* h, int bin, double weight = 1.0 )
0273   { h->SetBinContent(bin, h->GetBinContent(bin)+weight ); };
0274
0275   // increment total number of event
0276   increment( m_counters, TpotMonDefs::kEventCounter );
0277
0278   // check event and event type
0279   if( !event ) { return 0; }
0280   if(event->getEvtType() >= 8) { return 0; }
0281
0282   // increment total number of valid events
0283   ++m_evtcnt;
0284
0285   increment( m_counters, TpotMonDefs::kValidEventCounter );
0286
0287   // hit multiplicity vs fee id
0288   std::map<int, int> multiplicity;
0289
0290   // read the data
0291   for( const auto& packet_id:MicromegasDefs::m_packet_ids )
0292   {
0293     std::unique_ptr<Packet> packet(event->getPacket(packet_id));
0294     if( !packet )
0295     {
0296       // no data
0297       if( Verbosity() )
0298       { std::cout << "TpotMon::process_event - packet " << packet_id << " not found" << std::endl; }
0299       continue;
0300     }
0301
0302     // get number of lvl1 tagger
0303     const int n_tagger = packet->lValue(0, "N_TAGGER");
0304     int n_lvl1_tagger = 0;
0305     int n_heartbeat_tagger = 0;
0306     for (int t = 0; t < n_tagger; t++)
0307     {
0308       const bool is_lvl1_tagger( static_cast<uint8_t>(packet->lValue(t, "IS_LEVEL1_TRIGGER" )));
0309       if( is_lvl1_tagger )
0310       { ++n_lvl1_tagger; }
0311
0312       // also save hearbeat bco
0313       const bool is_modebit = static_cast<uint8_t>(packet->lValue(t, "IS_MODEBIT"));
0314       if( is_modebit )
0315       {
0316         // get modebits
0317         uint64_t modebits = static_cast<uint8_t>(packet->lValue(t, "MODEBITS"));
0318         if( modebits&(1<<ELINK_HEARTBEAT_T) )
0319         { ++n_heartbeat_tagger; }
0320       }
0321     }
0322
0323     // increment counters
0324     increment( m_counters, TpotMonDefs::kTriggerCounter, double(n_lvl1_tagger)/MicromegasDefs::m_npackets_active );
0325     increment( m_counters, TpotMonDefs::kHeartBeatCounter, double(n_heartbeat_tagger)/MicromegasDefs::m_npackets_active );
0326
0327     m_triggercnt += double(n_lvl1_tagger)/MicromegasDefs::m_npackets_active;
0328
0329     // get number of datasets (also call waveforms)
0330     const auto n_waveforms = packet->iValue(0, "NR_WF" );
0331
0332     if( Verbosity() )
0333     { std::cout << "TpotMon::process_event - n_waveforms: " << n_waveforms << std::endl; }
0334     for( int i=0; i<n_waveforms; ++i )
0335     {
0336
0337       // get waveform type
0338       const int type = packet->iValue(i, "TYPE");
0339
0340       // get channel
0341       const auto channel = packet->iValue( i, "CHANNEL" );
0342
0343       // bound check channel
0344       if( channel < 0 || channel >= MicromegasDefs::m_nchannels_fee )
0345       {
0346         if( Verbosity() )
0347         { std::cout << "TpotMon::process_event - invalid channel: " << channel << std::endl; }
0348         continue;
0349       }
0350
0351       // account for fiber swapping
0352       const int fee_id = packet->iValue(i, "FEE");
0353
0354       // get detector index from fee id
0355       const auto iter = m_detector_histograms.find( fee_id );
0356       if( iter == m_detector_histograms.end() )
0357       {
0358         std::cout << "TpotMon::process_event - invalid fee_id: " << fee_id << std::endl;
0359         continue;
0360       }
0361       const auto& detector_histograms = iter->second;
0362
0363       // strip
0364       const auto strip_index = m_mapping.get_physical_strip(fee_id, channel );
0365
0366       // heartbeat hits
0367       if( type == MicromegasDefs::HEARTBEAT_T )
0368       {
0369         // fill dedicated histogram, ignore the rest
0370         detector_histograms.m_heartbeat_vs_channel->Fill( channel );
0371         continue;
0372       }
0373
0374       // get channel rms and pedestal from calibration data
0375       const double pedestal = m_calibration_data.get_pedestal( fee_id, channel );
0376       const double rms = m_calibration_data.get_rms( fee_id, channel );
0377
0378       // get tile center, segmentation
0379       const auto& [tile_x, tile_y]  = m_tile_centers.at(fee_id);
0380       const auto segmentation = MicromegasDefs::getSegmentationType( m_mapping.get_hitsetkey(fee_id));
0381
0382       // fill 2D histograms ADC vs sample and hit charge vs sample
0383       const int samples = packet->iValue( i, "SAMPLES" );
0384       for( int is = 0; is < samples; ++is )
0385       {
0386         const uint16_t adc =  packet->iValue( i, is );
0387         if( adc == MicromegasDefs::m_adc_invalid ) continue;
0388         const bool is_signal = rms>0 && (adc > m_min_adc) && (adc> pedestal+m_n_sigma*rms);
0389         if( is_signal )
0390         {
0391           if( !is_masked(fee_id,channel) )
0392           { detector_histograms.m_counts_sample->Fill( is ); }
0393
0394           detector_histograms.m_sample_channel->Fill( strip_index, is);
0395         }
0396
0397         if( !is_masked(fee_id, channel) )
0398         {
0399           detector_histograms.m_adc_sample->Fill( is, adc );
0400           detector_histograms.m_hit_charge->Fill( adc );
0401         }
0402
0403         detector_histograms.m_adc_channel->Fill( strip_index, adc );
0404
0405       }
0406
0407       // fill waveform profile for this channel
0408       detector_histograms.m_wf_vs_channel->Fill( strip_index );
0409
0410       // define if hit is signal
0411       bool is_signal = false;
0412       for( int is = std::max<int>(0,m_sample_window_signal.first); is < std::min<int>(samples,m_sample_window_signal.second); ++is )
0413       {
0414         const uint16_t adc =  packet->iValue( i, is );
0415         if( adc == MicromegasDefs::m_adc_invalid ) continue;
0416         if( rms>0 && (adc > m_min_adc) && (adc>pedestal + m_n_sigma*rms) )
0417         {
0418           is_signal = true;
0419           break;
0420         }
0421       }
0422
0423       // fill hit profile for this channel
0424       if( is_signal )
0425       {
0426         detector_histograms.m_hit_vs_channel->Fill( strip_index );
0427
0428         // fill detector multiplicity
0429         if( !is_masked( fee_id, channel ) )
0430         {
0431           ++multiplicity[fee_id];
0432
0433           switch( segmentation )
0434           {
0435             case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0436             m_detector_multiplicity_z->Fill( tile_x, tile_y );
0437             m_resist_multiplicity_z->Fill( tile_x + MicromegasGeometry::m_tile_length*( double(strip_index)/MicromegasDefs::m_nchannels_fee - 0.5), tile_y );
0438             break;
0439
0440             case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0441             m_detector_multiplicity_phi->Fill( tile_x, tile_y );
0442             m_resist_multiplicity_phi->Fill( tile_x, tile_y + MicromegasGeometry::m_tile_width*( double(strip_index)/MicromegasDefs::m_nchannels_fee - 0.5) );
0443             break;
0444           }
0445         }
0446
0447       }
0448     }
0449   }
0450
0451   // fill hit multiplicities
0452   for( auto&& [fee_id, detector_histograms]:m_detector_histograms )
0453   { detector_histograms.m_hit_multiplicity->Fill( multiplicity[fee_id] ); }
0454
0455   // convert multiplicity histogram into occupancy
0456   auto copy_content = []( TH2Poly* source, TH2Poly* destination, double scale )
0457   {
0458     for( int bin = 0; bin < source->GetNumberOfBins(); ++bin )
0459     { destination->SetBinContent( bin+1, source->GetBinContent( bin+1 )*scale ); }
0460   };
0461
0462   copy_content( m_detector_multiplicity_z, m_detector_occupancy_z, 100./(m_triggercnt*MicromegasDefs::m_nchannels_fee) );
0463   copy_content( m_detector_multiplicity_phi, m_detector_occupancy_phi, 100./(m_triggercnt*MicromegasDefs::m_nchannels_fee) );
0464   copy_content( m_resist_multiplicity_z, m_resist_occupancy_z, 400./(m_triggercnt*MicromegasDefs::m_nchannels_fee) );
0465   copy_content( m_resist_multiplicity_phi, m_resist_occupancy_phi, 400./(m_triggercnt*MicromegasDefs::m_nchannels_fee) );
0466
0467   return 0;
0468 }
0469
0470 //________________________________
0471 int TpotMon::Reset()
0472 {
0473   // reset our internal counters
0474   m_evtcnt = 0;
0475   m_triggercnt = 0;
0476
0477   // reset all histograms
0478   for( TH1* h:{
0479     m_detector_multiplicity_z, m_detector_multiplicity_phi,
0480     m_detector_occupancy_z, m_detector_occupancy_phi,
0481     m_resist_multiplicity_z, m_resist_multiplicity_phi,
0482     m_resist_occupancy_z, m_resist_occupancy_phi } )
0483   { h->Reset(); }
0484
0485   for( const auto& [fee_id,hlist]:m_detector_histograms )
0486   {
0487     for( TH1* h:std::initializer_list<TH1*>{
0488       hlist.m_counts_sample,
0489       hlist.m_adc_sample,
0490       hlist.m_adc_channel,
0491       hlist.m_sample_channel,
0492       hlist.m_hit_charge,
0493       hlist.m_hit_multiplicity,
0494       hlist.m_wf_vs_channel,
0495       hlist.m_hit_vs_channel } )
0496     { h->Reset(); }
0497   }
0498
0499   return 0;
0500 }
0501
0502 //________________________________
0503 void TpotMon::setup_detector_bins(TH2Poly* h2)
0504 {
0505   // loop over tile centers
0506   for( size_t i = 0; i < m_geometry.get_ntiles(); ++i )
0507   {
0508     const auto boundaries = m_geometry.get_tile_boundaries( i );
0509     h2->AddBin( boundaries.size(), &get_x( boundaries )[0], &get_y( boundaries )[0] );
0510   }
0511 }
0512
0513 //________________________________
0514 void TpotMon::setup_resist_bins(TH2Poly* h2, MicromegasDefs::SegmentationType segmentation)
0515 {
0516   // loop over tile centers
0517   for( size_t itile = 0; itile < m_geometry.get_ntiles(); ++itile )
0518   {
0519     for( size_t iresist = 0; iresist < MicromegasGeometry::m_nresist; ++iresist )
0520     {
0521       const auto boundaries = m_geometry.get_resist_boundaries( itile, iresist, segmentation );
0522       h2->AddBin( boundaries.size(), &get_x( boundaries )[0], &get_y( boundaries )[0] );
0523     }
0524   }
0525 }