packages/micromegas/MicromegasCombinedDataEvaluation.cc

0001 /*!
0002  * \file MicromegasCombinedDataEvaluation.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "MicromegasCombinedDataEvaluation.h"
0007 #include "MicromegasDefs.h"
0008
0009 #include <ffarawobjects/MicromegasRawHit.h>
0010 #include <ffarawobjects/MicromegasRawHitContainer.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/getClass.h>
0016
0017 #include <TFile.h>
0018 #include <TH1.h>
0019 #include <TH2.h>
0020 #include <TProfile.h>
0021
0022 #include <cassert>
0023 #include <fstream>
0024 #include <list>
0025 #include <memory>
0026
0027 namespace
0028 {
0029
0030   // streamer for lists
0031   template <class T>
0032   std::ostream& operator<<(std::ostream& out, const std::list<T>& list)
0033   {
0034     if (list.empty())
0035     {
0036       out << "{}";
0037     }
0038     else
0039     {
0040       out << "{ ";
0041       bool first = true;
0042       for (const auto& value : list)
0043       {
0044         if (!first)
0045         {
0046           out << ", ";
0047         }
0048         out << value;
0049         first = false;
0050       }
0051
0052       out << " }";
0053     }
0054
0055     return out;
0056   }
0057
0058 }  // namespace
0059
0060 //_________________________________________________________
0061 void MicromegasCombinedDataEvaluation::Waveform::copy_from(const MicromegasCombinedDataEvaluation::Sample& sample)
0062 {
0063   packet_id = sample.packet_id;
0064   lvl1_bco = sample.lvl1_bco;
0065   fee_bco = sample.fee_bco;
0066   checksum = sample.checksum;
0067   checksum_error = sample.checksum_error;
0068   fee_id = sample.fee_id;
0069   layer = sample.layer;
0070   tile = sample.tile;
0071   sampa_address = sample.sampa_address;
0072   sampa_channel = sample.sampa_channel;
0073   channel = sample.channel;
0074   strip = sample.strip;
0075   sample_max = sample.sample;
0076   adc_max = sample.adc;
0077   pedestal = sample.pedestal;
0078   rms = sample.rms;
0079 }
0080
0081 //_________________________________________________________
0082 void MicromegasCombinedDataEvaluation::Container::Reset()
0083 {
0084   n_tagger.clear();
0085   n_waveform.clear();
0086   samples.clear();
0087   waveforms.clear();
0088   lvl1_bco_list.clear();
0089   lvl1_count_list.clear();
0090 }
0091
0092 //_________________________________________________________
0093 MicromegasCombinedDataEvaluation::MicromegasCombinedDataEvaluation(const std::string& name)
0094   : SubsysReco(name)
0095 {
0096 }
0097
0098 //_____________________________________________________________________
0099 int MicromegasCombinedDataEvaluation::Init(PHCompositeNode* /*topNode*/)
0100 {
0101   // read calibrations
0102   m_calibration_data.read(m_calibration_filename);
0103
0104   m_evaluation_file.reset(new TFile(m_evaluation_filename.c_str(), "RECREATE"));
0105   m_evaluation_tree = new TTree("T", "T");
0106   m_evaluation_tree->SetAutoSave(5000);
0107   m_container = new Container;
0108   m_evaluation_tree->Branch("Event", &m_container);
0109   return Fun4AllReturnCodes::EVENT_OK;
0110 }
0111
0112 //____________________________________________________________________________..
0113 int MicromegasCombinedDataEvaluation::InitRun(PHCompositeNode* /*topNode*/)
0114 {
0115   return Fun4AllReturnCodes::EVENT_OK;
0116 }
0117
0118 //___________________________________________________________________________
0119 int MicromegasCombinedDataEvaluation::process_event(PHCompositeNode* topNode)
0120 {
0121   // load raw hits container
0122   auto rawhitcontainer = findNode::getClass<MicromegasRawHitContainer>(topNode, m_rawhitnodename);
0123   assert(rawhitcontainer);
0124
0125   // reset container
0126   m_container->Reset();
0127
0128   // temporary storage for samples and waveforms, sorted by lvl1 bco
0129   std::multimap<uint64_t, Sample> sample_map;
0130   std::multimap<uint64_t, Waveform> waveform_map;
0131
0132   // map number of waveforms per packet
0133   std::map<unsigned int, size_t> packet_waveforms;
0134
0135   // loop over raw hits
0136   if (Verbosity())
0137   {
0138     std::cout << "MicromegasCombinedDataEvaluation::process_event - hits: " << rawhitcontainer->get_nhits() << std::endl;
0139   }
0140
0141   bool first = true;
0142   uint64_t first_lvl1_bco = 0;
0143
0144   for (unsigned int ihit = 0; ihit < rawhitcontainer->get_nhits(); ++ihit)
0145   {
0146     const auto rawhit = rawhitcontainer->get_hit(ihit);
0147     const auto packet_id = rawhit->get_packetid();
0148
0149     // make sure packet is valid
0150     if (std::find(std::begin(MicromegasDefs::m_packet_ids), std::end(MicromegasDefs::m_packet_ids), packet_id) == std::end(MicromegasDefs::m_packet_ids))
0151     {
0152       std::cout << "MicromegasCombinedDataEvaluation::process_event - invalid packet: " << packet_id << std::endl;
0153       continue;
0154     }
0155
0156     ++packet_waveforms[packet_id];
0157
0158     // create running sample, assign packet, fee, layer and tile id
0159     Sample sample;
0160     sample.packet_id = packet_id;
0161
0162     // get fee id, apply mapping to current fiber set, for backward compatibility
0163     sample.fee_id = m_mapping.get_new_fee_id(rawhit->get_fee());
0164
0165     const auto hitsetkey = m_mapping.get_hitsetkey(sample.fee_id);
0166     sample.layer = TrkrDefs::getLayer(hitsetkey);
0167     sample.tile = MicromegasDefs::getTileId(hitsetkey);
0168
0169     // channel and bound check.
0170     sample.channel = rawhit->get_channel();
0171     if( sample.channel >= MicromegasDefs::m_nchannels_fee )
0172     {
0173       if( Verbosity() )
0174       { std::cout << "MicromegasCombinedDataEvaluation::process_event - invalid channel: " << sample.channel << std::endl; }
0175       continue;
0176     }
0177
0178     // beam crossing
0179     sample.fee_bco = rawhit->get_bco();
0180     sample.lvl1_bco = rawhit->get_gtm_bco();
0181
0182     if (first)
0183     {
0184       first = false;
0185       first_lvl1_bco = rawhit->get_gtm_bco();
0186     }
0187
0188     // increment bco map
0189     // ++m_bco_map[sample.lvl1_bco];
0190
0191     ++m_bco_map[first_lvl1_bco];
0192
0193     //     // checksum and checksum error
0194     //     sample.checksum = rawhit->get_checksum();
0195     //     sample.checksum_error = rawhit->get_checksum_error();
0196
0197     // channel, sampa_channel, sampa address and strip
0198     sample.sampa_address = rawhit->get_sampaaddress();
0199     sample.sampa_channel = rawhit->get_sampachannel();
0200     sample.strip = m_mapping.get_physical_strip(sample.fee_id, sample.channel);
0201
0202     // get channel rms and pedestal from calibration data
0203     const double pedestal = m_calibration_data.get_pedestal(sample.fee_id, sample.channel);
0204     const double rms = m_calibration_data.get_rms(sample.fee_id, sample.channel);
0205     sample.pedestal = pedestal;
0206     sample.rms = rms;
0207
0208     // get number of samples and loop
0209     const auto sample_range = std::make_pair( rawhit->get_sample_begin(), rawhit->get_sample_end() );
0210     if (Verbosity() > 1)
0211     {
0212       std::cout << "MicromegasCombinedDataEvaluation::process_event -"
0213                 << " fee: " << sample.fee_id
0214                 << " tile: " << sample.tile
0215                 << " layer: " << sample.layer
0216                 << " tile: " << sample.tile
0217                 << " lvl1_bco: " << sample.lvl1_bco
0218                 << " fee_bco: " << sample.fee_bco
0219                 << " error: " << sample.checksum_error
0220                 << " channel: " << sample.channel
0221                 << " strip: " << sample.strip
0222                 << " samples: (" << sample_range.first << "," << sample_range.second << ")"
0223                 << std::endl;
0224     }
0225
0226     Sample sample_max;
0227     for (unsigned short is = sample_range.first; is < std::min<unsigned short>(sample_range.second, 1024); ++is)
0228     {
0229       // assign sample id and corresponding adc, save copy in container
0230       const auto adc = rawhit->get_adc(is);
0231       if (adc == MicromegasDefs::m_adc_invalid)
0232       {
0233         continue;
0234       }
0235       sample.sample = is;
0236       sample.adc = adc;
0237       sample_map.emplace(sample.lvl1_bco, sample);
0238
0239       if (sample.adc > sample_max.adc)
0240       {
0241         sample_max = sample;
0242       }
0243     }
0244
0245     // create waveform
0246     Waveform waveform(sample_max);
0247     waveform.is_signal =
0248         rms > 0 &&
0249         waveform.adc_max >= m_min_adc &&
0250         waveform.sample_max >= m_sample_min &&
0251         waveform.sample_max < m_sample_max &&
0252         waveform.adc_max > pedestal + m_n_sigma * rms;
0253
0254     waveform_map.emplace(waveform.lvl1_bco, waveform);
0255   }
0256
0257   // copy all samples and waveform to container
0258   for (auto&& [lvl_bco, sample] : sample_map)
0259   //  { m_container->samples.push_back(std::move(sample)); } // std::move has no effect for trivially copyable types
0260   {
0261     m_container->samples.push_back(sample);
0262   }
0263
0264   for (auto&& [lvl1_bco, waveform] : waveform_map)
0265   //  { m_container->waveforms.push_back(std::move(waveform)); } // std::move has no effect for trivially copyable types
0266   {
0267     m_container->waveforms.push_back(waveform);
0268   }
0269
0270   // store number of waveforms
0271   for (const auto& [packet_id, n_waveforms] : packet_waveforms)
0272   {
0273     m_container->n_waveform.push_back(n_waveforms);
0274   }
0275
0276   // fill evaluation tree
0277   m_evaluation_tree->Fill();
0278
0279   return Fun4AllReturnCodes::EVENT_OK;
0280 }
0281
0282 //_____________________________________________________________________
0283 int MicromegasCombinedDataEvaluation::End(PHCompositeNode* /*topNode*/)
0284 {
0285   if (m_evaluation_file && m_evaluation_tree)
0286   {
0287     m_evaluation_file->cd();
0288     m_evaluation_tree->Write();
0289     m_evaluation_file->Close();
0290   }
0291
0292   // print bco map
0293   if (Verbosity())
0294   {
0295     for (const auto& [bco, nwaveforms] : m_bco_map)
0296     {
0297       std::cout << "MicromegasCombinedDataEvaluation::End - bco: 0x" << std::hex << bco << std::dec << ", nwaveforms: " << nwaveforms << std::endl;
0298     }
0299   }
0300
0301   // print bco list, for offline processing
0302   if (Verbosity())
0303   {
0304     std::cout << "const std::vector<uint64_t> lvl1_bco_list = {" << std::endl;
0305     bool first = true;
0306     int count = 0;
0307     for (const auto& [bco, nwaveforms] : m_bco_map)
0308     {
0309       if (!first)
0310       {
0311         std::cout << ", ";
0312       }
0313       first = false;
0314       if (count == 10)
0315       {
0316         count = 0;
0317         std::cout << std::endl;
0318       }
0319       std::cout << " 0x" << std::hex << bco << std::dec;
0320       ++count;
0321     }
0322     std::cout << std::endl
0323               << "};" << std::endl;
0324   }
0325
0326   return Fun4AllReturnCodes::EVENT_OK;
0327 }