packages/micromegas/MicromegasClusterizer.cc

0001 /*!
0002  * \file MicromegasClusterizer.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "MicromegasClusterizer.h"
0007 #include "MicromegasDefs.h"
0008 #include "CylinderGeomMicromegas.h"
0009
0010 #include <g4detectors/PHG4CylinderGeomContainer.h>
0011 #include <g4detectors/PHG4CylinderGeom.h>           // for PHG4CylinderGeom
0012
0013 #include <trackbase/ActsGeometry.h>
0014 #include <trackbase/TrkrClusterContainerv4.h>        // for TrkrCluster
0015 #include <trackbase/TrkrClusterv5.h>
0016 #include <trackbase/TrkrDefs.h>
0017 #include <trackbase/TrkrHitSet.h>
0018 #include <trackbase/TrkrHit.h>
0019 #include <trackbase/TrkrHitSetContainer.h>
0020 #include <trackbase/TrkrClusterHitAssocv3.h>
0021
0022 #include <Acts/Definitions/Units.hpp>
0023 #include <Acts/Surfaces/Surface.hpp>
0024
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026 #include <fun4all/SubsysReco.h>                     // for SubsysReco
0027
0028 #include <phool/getClass.h>
0029 #include <phool/PHCompositeNode.h>
0030 #include <phool/PHIODataNode.h>                     // for PHIODataNode
0031 #include <phool/PHNode.h>                           // for PHNode
0032 #include <phool/PHNodeIterator.h>                   // for PHNodeIterator
0033 #include <phool/PHObject.h>                         // for PHObject
0034
0035 #include <Eigen/Dense>
0036
0037 #include <TVector3.h>
0038
0039 #include <cassert>
0040 #include <cmath>
0041 #include <cstdint>                                 // for uint16_t
0042 #include <iterator>                                 // for distance
0043 #include <map>                                      // for _Rb_tree_const_it...
0044 #include <utility>                                  // for pair, make_pair
0045 #include <vector>
0046
0047
0048 namespace
0049 {
0050   //! convenience square method
0051   template<class T>
0052     inline constexpr T square( const T& x ) { return x*x; }
0053
0054   // streamers
0055   [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const Acts::Vector3& vector )
0056   {
0057     out << "( " << vector[0] << "," << vector[1] << "," << vector[2] << ")";
0058     return out;
0059   }
0060
0061   // streamers
0062   [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const Acts::Vector2& vector )
0063   {
0064     out << "( " << vector[0] << "," << vector[1] << ")";
0065     return out;
0066   }
0067
0068   // streamers
0069   [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const TVector3& vector )
0070   {
0071     out << "( " << vector.x() << "," << vector.y() << "," << vector.z() << ")";
0072     return out;
0073   }
0074
0075   // streamers
0076   [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const TVector2& vector )
0077   {
0078     out << "( " << vector.X() << "," << vector.Y() << ")";
0079     return out;
0080   }
0081
0082 }
0083
0084 //_______________________________________________________________________________
0085 MicromegasClusterizer::MicromegasClusterizer(const std::string &name )
0086   : SubsysReco(name)
0087 {}
0088
0089 //_____________________________________________________________________
0090 int MicromegasClusterizer::Init(PHCompositeNode* /*topNode*/ )
0091 {
0092   // print configuration
0093   std::cout << "MicromegasClusterizer::Init - m_use_default_pedestal: " << m_use_default_pedestal << std::endl;
0094   std::cout << "MicromegasClusterizer::Init - m_default_pedestal: " << m_default_pedestal << std::endl;
0095   std::cout
0096     << "MicromegasClusterizer::Init -"
0097     << " m_calibration_filename: "
0098     << (m_calibration_filename.empty() ? "unspecified":m_calibration_filename )
0099     << std::endl;
0100
0101   // read calibrations
0102   if( !m_calibration_filename.empty() )
0103   { m_calibration_data.read( m_calibration_filename ); }
0104
0105   return Fun4AllReturnCodes::EVENT_OK;
0106 }
0107
0108 //_______________________________________________________________________________
0109 int MicromegasClusterizer::InitRun(PHCompositeNode *topNode)
0110 {
0111   PHNodeIterator iter(topNode);
0112
0113   // Looking for the DST node
0114   auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0115   assert( dstNode );
0116
0117   // Create the Cluster node if missing
0118   auto trkrClusterContainer = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
0119   if (!trkrClusterContainer)
0120   {
0121     PHNodeIterator dstiter(dstNode);
0122     auto trkrNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0123     if(!trkrNode)
0124     {
0125       trkrNode = new PHCompositeNode("TRKR");
0126       dstNode->addNode(trkrNode);
0127     }
0128
0129     trkrClusterContainer = new TrkrClusterContainerv4;
0130     auto TrkrClusterContainerNode = new PHIODataNode<PHObject>(trkrClusterContainer, "TRKR_CLUSTER", "PHObject");
0131     trkrNode->addNode(TrkrClusterContainerNode);
0132   }
0133
0134   // create cluster to hit association node, if missing
0135   auto trkrClusterHitAssoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,"TRKR_CLUSTERHITASSOC");
0136   if(!trkrClusterHitAssoc)
0137   {
0138     PHNodeIterator dstiter(dstNode);
0139     auto trkrNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0140     if(!trkrNode)
0141     {
0142       trkrNode = new PHCompositeNode("TRKR");
0143       dstNode->addNode(trkrNode);
0144     }
0145
0146     trkrClusterHitAssoc = new TrkrClusterHitAssocv3;
0147     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(trkrClusterHitAssoc, "TRKR_CLUSTERHITASSOC", "PHObject");
0148     trkrNode->addNode(newNode);
0149   }
0150   return Fun4AllReturnCodes::EVENT_OK;
0151 }
0152
0153 //_______________________________________________________________________________
0154 int MicromegasClusterizer::process_event(PHCompositeNode *topNode)
0155 {
0156
0157   // geometry
0158   PHG4CylinderGeomContainer* geonode = nullptr;
0159   for( std::string geonodename: {"CYLINDERGEOM_MICROMEGAS_FULL", "CYLINDERGEOM_MICROMEGAS" } )
0160   { if(( geonode =  findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str()) )) { break;
0161 }}
0162   assert(geonode);
0163
0164   // hitset container
0165   auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0166   assert( trkrhitsetcontainer );
0167
0168   // cluster container
0169   auto trkrClusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0170   assert( trkrClusterContainer );
0171
0172   // cluster-hit association
0173   auto trkrClusterHitAssoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0174   assert( trkrClusterHitAssoc );
0175
0176   // geometry
0177   auto acts_geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0178   assert( acts_geometry );
0179
0180   // loop over micromegas hitsets
0181   const auto hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::micromegasId);
0182   for( auto hitset_it = hitset_range.first; hitset_it != hitset_range.second; ++hitset_it )
0183   {
0184     // get hitset, key and layer
0185     TrkrHitSet* hitset = hitset_it->second;
0186     const TrkrDefs::hitsetkey hitsetkey = hitset_it->first;
0187     const auto layer = TrkrDefs::getLayer(hitsetkey);
0188     const auto tileid = MicromegasDefs::getTileId(hitsetkey);
0189
0190     // get micromegas geometry object
0191     const auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(geonode->GetLayerGeom(layer));
0192     assert(layergeom);
0193
0194     // get micromegas acts surface
0195     const auto acts_surface = acts_geometry->maps().getMMSurface( hitsetkey);
0196     if( !acts_surface )
0197     {
0198       std::cout
0199         << "MicromegasClusterizer::process_event -"
0200         << " could not find surface for layer " << (int) layer << " tile: " << (int) tileid
0201         << " skipping hitset"
0202         << std::endl;
0203       continue;
0204     }
0205
0206     /*
0207      * get segmentation type, layer thickness, strip length and pitch.
0208      * They are used to calculate cluster errors
0209      */
0210     const auto segmentation_type = layergeom->get_segmentation_type();
0211     const double pitch = layergeom->get_pitch();
0212     const double strip_length = layergeom->get_strip_length( tileid, acts_geometry );
0213
0214     // keep a list of ranges corresponding to each cluster
0215     using range_list_t = std::vector<TrkrHitSet::ConstRange>;
0216     range_list_t ranges;
0217
0218     // loop over hits
0219     const auto hit_range = hitset->getHits();
0220
0221     // keep track of first iterator of runing cluster
0222     auto begin = hit_range.first;
0223
0224     // keep track of previous strip
0225     uint16_t previous_strip = 0;
0226     bool first = true;
0227
0228     for( auto hit_it = hit_range.first; hit_it != hit_range.second; ++hit_it )
0229     {
0230
0231       // get hit key
0232       const auto hitkey = hit_it->first;
0233
0234       // get strip number
0235       const auto strip = MicromegasDefs::getStrip( hitkey );
0236
0237       if( first )
0238       {
0239
0240         previous_strip = strip;
0241         first = false;
0242         continue;
0243
0244       } else if( strip - previous_strip > 1 ) {
0245
0246         // store current cluster range
0247         ranges.push_back( std::make_pair( begin, hit_it ) );
0248
0249         // reinitialize begin of next cluster range
0250         begin = hit_it;
0251
0252       }
0253
0254       // update previous strip
0255       previous_strip = strip;
0256
0257     }
0258
0259     // store last cluster
0260     if( begin != hit_range.second ) { ranges.push_back( std::make_pair( begin, hit_range.second ) );
0261 }
0262
0263     // initialize cluster count
0264     int cluster_count = 0;
0265
0266     // loop over found hit ranges and create clusters
0267     for( const auto& range : ranges )
0268     {
0269       // create cluster key and corresponding cluster
0270       const auto ckey = TrkrDefs::genClusKey( hitsetkey, cluster_count++ );
0271
0272       TVector2 local_coordinates;
0273       double weight_sum = 0;
0274
0275       // needed for proper error calculation
0276       // it is either the sum over z, or phi, depending on segmentation
0277       double coord_sum = 0;
0278       double coordsquare_sum = 0;
0279
0280       // also store adc value
0281       unsigned int adc_sum = 0;
0282       unsigned int max_adc = 0;
0283       const unsigned int strip_count = std::distance(range.first,range.second);
0284       if(m_drop_single_strips && strip_count < 2)
0285       { continue; }
0286
0287       // loop over constituting hits
0288       for( auto hit_it = range.first; hit_it != range.second; ++hit_it )
0289       {
0290         // get hit key
0291         const auto hitkey = hit_it->first;
0292         const auto hit = hit_it->second;
0293
0294         // associate cluster key to hit key
0295         trkrClusterHitAssoc->addAssoc(ckey, hitkey );
0296
0297         // get strip number
0298         const auto strip = MicromegasDefs::getStrip( hitkey );
0299
0300         // get adc, remove pedestal
0301         const double pedestal = m_use_default_pedestal ?
0302           m_default_pedestal:
0303           m_calibration_data.get_pedestal_mapped(hitsetkey, strip);
0304         const double weight = double(hit->getAdc()) - pedestal;
0305
0306         // increment cluster adc
0307         const auto hit_adc = hit->getAdc();
0308         if( hit_adc > max_adc) { max_adc = hit_adc; }
0309         adc_sum += hit_adc;
0310
0311         // get strip local coordinate and update relevant sums
0312         const auto strip_local_coordinate = layergeom->get_local_coordinates( tileid, acts_geometry, strip );
0313         local_coordinates += strip_local_coordinate*weight;
0314         switch( segmentation_type )
0315         {
0316           case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0317           {
0318
0319             coord_sum += strip_local_coordinate.X()*weight;
0320             coordsquare_sum += square(strip_local_coordinate.X())*weight;
0321             break;
0322           }
0323
0324           case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0325           {
0326             coord_sum += strip_local_coordinate.Y()*weight;
0327             coordsquare_sum += square(strip_local_coordinate.Y())*weight;
0328             break;
0329           }
0330         }
0331
0332         weight_sum += weight;
0333
0334       }
0335
0336       local_coordinates *= (1./weight_sum);
0337
0338       // dimension and error in r, rphi and z coordinates
0339       static const float invsqrt12 = 1./std::sqrt(12);
0340       static constexpr float error_scale_phi = 1.6;
0341       static constexpr float error_scale_z = 0.8;
0342
0343       auto coord_cov = coordsquare_sum/weight_sum - square( coord_sum/weight_sum );
0344       auto coord_error_sq = coord_cov/weight_sum;
0345
0346       // local errors (x is along rphi, y is along z)
0347       double error_sq_x = 0;
0348       double error_sq_y = 0;
0349       switch( segmentation_type )
0350       {
0351         case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0352         {
0353           if( coord_error_sq == 0 ) { coord_error_sq = square(pitch)/12;
0354           } else { coord_error_sq *= square(error_scale_phi);
0355 }
0356           error_sq_x = coord_error_sq;
0357           error_sq_y = square(strip_length*invsqrt12);
0358           break;
0359         }
0360
0361         case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0362         {
0363           if( coord_error_sq == 0 ) { coord_error_sq = square(pitch)/12;
0364           } else { coord_error_sq *= square(error_scale_z);
0365 }
0366           error_sq_x = square(strip_length*invsqrt12);
0367           error_sq_y = coord_error_sq;
0368           break;
0369         }
0370       }
0371
0372       auto cluster = std::make_unique<TrkrClusterv5>();
0373       cluster->setAdc( adc_sum );
0374       cluster->setMaxAdc( max_adc );
0375       cluster->setLocalX(local_coordinates.X());
0376       cluster->setLocalY(local_coordinates.Y());
0377       cluster->setPhiError(sqrt(error_sq_x));
0378       cluster->setZError(sqrt(error_sq_y));
0379
0380       // store cluster size
0381       switch( segmentation_type )
0382       {
0383         case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0384         {
0385           cluster->setPhiSize(strip_count);
0386           cluster->setZSize(1);
0387           break;
0388         }
0389
0390         case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0391         {
0392           cluster->setPhiSize(1);
0393           cluster->setZSize(strip_count);
0394           break;
0395         }
0396       }
0397
0398       trkrClusterContainer->addClusterSpecifyKey( ckey, cluster.release() );
0399
0400       // increment counter
0401       ++m_clustercounts[hitsetkey];
0402
0403     }
0404
0405   }
0406   // done
0407   return Fun4AllReturnCodes::EVENT_OK;
0408 }
0409
0410 //_____________________________________________________________________
0411 int MicromegasClusterizer::End(PHCompositeNode* /*topNode*/)
0412 {
0413   // if( Verbosity() )
0414   {
0415     for (const auto& [hitsetkey, count] : m_clustercounts)
0416     {
0417       std::cout << "MicromegasClusterizer::End - hitsetkey: " << hitsetkey << ", cluster count: " << count << std::endl;
0418     }
0419   }
0420
0421   return Fun4AllReturnCodes::EVENT_OK;
0422 }