Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:18:30 

0001 #include "TrackDeadLayerTagging.h"
0002 
0003 #include <phool/phool.h>
0004 #include <phool/getClass.h>
0005 #include <fun4all/PHTFileServer.h>
0006 #include <fun4all/Fun4AllServer.h>
0007 
0008 #include <cdbobjects/CDBTTree.h>
0009 #include <ffamodules/CDBInterface.h>
0010 
0011 #include <trackreco/ActsPropagator.h>
0012 
0013 #include <trackbase/TrkrCluster.h>
0014 #include <trackbase/MvtxDefs.h>
0015 #include <trackbase/InttDefs.h>
0016 #include <trackbase/TpcDefs.h>
0017 #include <trackbase/TrackFitUtils.h>
0018 
0019 #include <g4detectors/PHG4TpcGeom.h>
0020 
0021 #include <mvtx/SegmentationAlpide.h>
0022 #include <intt/CylinderGeomIntt.h>
0023 #include <intt/InttMapping.h>
0024 #include <micromegas/MicromegasDefs.h>
0025 #include <micromegas/CylinderGeomMicromegas.h>
0026 
0027 #include <TNtuple.h>
0028 #include <TH1.h>
0029 
0030 #include <iostream>
0031 
0032 //____________________________________
0033 TrackDeadLayerTagging::TrackDeadLayerTagging(const std::string &name):
0034     SubsysReco(name)
0035 { 
0036   //initialize
0037   _event = 0;
0038   _outfile = "dead_layers.root";
0039 }
0040 
0041 //___________________________________
0042 int TrackDeadLayerTagging::InitRun(PHCompositeNode *topNode)
0043 {
0044   std::cout << PHWHERE << " Opening file " << _outfile << std::endl;
0045   PHTFileServer::get().open( _outfile, "RECREATE");
0046 
0047   _missing_sensors = new TTree("missing_sensors","Sensors missing from track");
0048   _missing_sensors->Branch("event",&_event);
0049   _missing_sensors->Branch("trackID",&_trackID);
0050   _missing_sensors->Branch("layer",&_layer);
0051   _missing_sensors->Branch("phiID",&_phiID);
0052   _missing_sensors->Branch("zID",&_zID);
0053   _missing_sensors->Branch("localindX",&_localindX);
0054   _missing_sensors->Branch("localindY",&_localindY);
0055   _missing_sensors->Branch("isDead",&_isDead);
0056   _missing_sensors->Branch("global_x",&_global_x);
0057   _missing_sensors->Branch("global_y",&_global_y);
0058   _missing_sensors->Branch("global_z",&_global_z);
0059 
0060   std::cout << "loading nodes" << std::endl;
0061   tpcglobal.set_verbosity(5);
0062   tpcglobal.loadNodes(topNode);
0063   std::cout << "nodes loaded" << std::endl;
0064   //_outtree = new TNtuple("deadlayers_test","deadlayers_test","x:y:z:pt:pz:eta:phi:charge:nhits:nmaps:nintt:quality:crossing:vtx_x:vtx_y:vtx_z:vtx_ntracks:vtxid:ndead_mvtx:ndead_intt:ndead_tpc:ndead_tpot");
0065 
0066   return 0;
0067 }
0068 
0069 //__________________________________
0070 //Call user instructions for every event
0071 int TrackDeadLayerTagging::process_event(PHCompositeNode *topNode)
0072 {
0073   _event++;
0074   if(_event%1000==0) std::cout << PHWHERE << "Events processed: " << _event << std::endl;
0075  
0076   _surface_maps.clear(); 
0077   GetNodes(topNode);
0078   for(auto& v : _dead_sensors) v.clear();
0079   get_surfaces();
0080   get_dead_sensor_maps();
0081 
0082   if(Verbosity()>0)
0083   {
0084     std::cout << "--------------------------------" << std::endl;
0085     std::cout << "event " << _event << std::endl;
0086   }
0087 
0088   find_dead_layers(topNode);
0089 
0090   return 0;
0091 }
0092 
0093 void TrackDeadLayerTagging::get_dead_sensor_maps()
0094 {
0095   std::string mvtx_hotpixel = CDBInterface::instance()->getUrl("MVTX_HotPixelMap");
0096   std::string mvtx_deadpixel = CDBInterface::instance()->getUrl("MVTX_DeadPixelMap");
0097   std::string intt_hotstrip = CDBInterface::instance()->getUrl("INTT_HotMap");
0098   std::string tpc_hotchannel = CDBInterface::instance()->getUrl("TPC_HOTCHANNELMAP");
0099   std::string tpc_deadchannel = CDBInterface::instance()->getUrl("TPC_DEADCHANNELMAP");
0100 
0101   std::unique_ptr<CDBTTree> mvtx_hottree = std::make_unique<CDBTTree>(mvtx_hotpixel);
0102   std::unique_ptr<CDBTTree> mvtx_deadtree = std::make_unique<CDBTTree>(mvtx_deadpixel);
0103   std::unique_ptr<CDBTTree> intt_hottree = std::make_unique<CDBTTree>(intt_hotstrip);
0104   std::unique_ptr<CDBTTree> tpc_hottree;
0105   std::unique_ptr<CDBTTree> tpc_deadtree;
0106   if(!tpc_hotchannel.empty())
0107   {
0108     tpc_hottree = std::make_unique<CDBTTree>(tpc_hotchannel);
0109   }
0110   if(!tpc_deadchannel.empty())
0111   {
0112     tpc_deadtree = std::make_unique<CDBTTree>(tpc_deadchannel);
0113   }
0114 
0115   std::cout << "Verbosity is " << Verbosity() << std::endl;
0116 
0117   mvtx_hottree->LoadCalibrations();
0118   mvtx_deadtree->LoadCalibrations();
0119   intt_hottree->LoadCalibrations();
0120   if(!tpc_hotchannel.empty()) tpc_hottree->LoadCalibrations();
0121   if(!tpc_deadchannel.empty()) tpc_deadtree->LoadCalibrations();
0122 
0123   int n_mvtxhot = mvtx_hottree->GetSingleIntValue("TotalHotPixels");
0124   int n_mvtxdead = mvtx_deadtree->GetSingleIntValue("TotalDeadPixels");
0125   int n_intthot = intt_hottree->GetSingleIntValue("size");
0126   int n_tpchot = (!tpc_hotchannel.empty())? tpc_hottree->GetSingleIntValue("TotalHotChannels") : 0;
0127   int n_tpcdead = (!tpc_deadchannel.empty())? tpc_deadtree->GetSingleIntValue("TotalDeadChannels") : 0;
0128 
0129   if(Verbosity()>0)
0130   {
0131     std::cout << "MVTX has " << n_mvtxhot << " hot pixels" << std::endl;
0132     std::cout << "MVTX has " << n_mvtxdead << " dead pixels" << std::endl;
0133     std::cout << "INTT has " << n_intthot << " hot pixels" << std::endl;
0134     std::cout << "TPC has " << n_tpchot << " hot pixels" << std::endl;
0135     std::cout << "TPC has " << n_tpcdead << " dead pixels" << std::endl;
0136   }
0137 
0138   if(Verbosity()>2) std::cout << "MVTX hot sensors:" << std::endl;
0139   for(int i=0; i<n_mvtxhot; i++)
0140   {
0141     int layer = mvtx_hottree->GetIntValue(i,"layer");
0142     MissingSensor hotsensor_mvtx = {
0143       .layer = layer,
0144       .phi_id = mvtx_hottree->GetIntValue(i,"stave"),
0145       .z_id = mvtx_hottree->GetIntValue(i,"chip"),
0146       .local_x = mvtx_hottree->GetIntValue(i,"row"),
0147       .local_y = mvtx_hottree->GetIntValue(i,"col")
0148     };
0149     if(Verbosity()>2)
0150     {
0151       std::cout << "layer " << hotsensor_mvtx.layer << " stave " << hotsensor_mvtx.phi_id << " chip " << hotsensor_mvtx.z_id << " row " << hotsensor_mvtx.local_x << " col " << hotsensor_mvtx.local_y << std::endl;
0152     }
0153     if(layer>=0&&layer<=2)
0154     {
0155       _dead_sensors[layer].push_back(hotsensor_mvtx);
0156     }
0157   }
0158 
0159   if(Verbosity()>2) std::cout << "MVTX dead sensors:" << std::endl;
0160   for(int i=0; i<n_mvtxdead; i++)
0161   {
0162     int layer = mvtx_deadtree->GetIntValue(i,"layer");
0163     MissingSensor deadsensor_mvtx = {
0164       .layer = layer,
0165       .phi_id = mvtx_deadtree->GetIntValue(i,"stave"),
0166       .z_id = mvtx_deadtree->GetIntValue(i,"chip"),
0167       .local_x = mvtx_deadtree->GetIntValue(i,"row"),
0168       .local_y = mvtx_deadtree->GetIntValue(i,"col")
0169     };
0170     if(layer>=0&&layer<=2)
0171     {
0172       _dead_sensors[layer].push_back(deadsensor_mvtx);
0173     }
0174     if(Verbosity()>2)
0175     {
0176       std::cout << "layer " << deadsensor_mvtx.layer << " stave " << deadsensor_mvtx.phi_id << " chip " << deadsensor_mvtx.z_id << " row " << deadsensor_mvtx.local_x << " col " << deadsensor_mvtx.local_y << std::endl;
0177     }
0178   }
0179 
0180   if(Verbosity()>2) std::cout << "INTT hot sensors:" << std::endl;
0181   for(int i=0; i<n_intthot; i++)
0182   {
0183     InttNameSpace::RawData_s raw_sensor;
0184     raw_sensor.felix_server = intt_hottree->GetIntValue(i,"felix_server");
0185     raw_sensor.felix_channel = intt_hottree->GetIntValue(i,"felix_channel");
0186     raw_sensor.chip = intt_hottree->GetIntValue(i,"chip");
0187     raw_sensor.channel = intt_hottree->GetIntValue(i,"channel");
0188     InttNameSpace::Offline_s offline_sensor = InttNameSpace::ToOffline(raw_sensor);
0189     MissingSensor hotsensor_intt = {
0190       .layer = offline_sensor.layer,
0191       .phi_id = offline_sensor.ladder_phi,
0192       .z_id = offline_sensor.ladder_z,
0193       .local_x = offline_sensor.strip_y,
0194       .local_y = offline_sensor.strip_x
0195     };
0196     if(offline_sensor.layer>=3&&offline_sensor.layer<=6)
0197     {
0198       _dead_sensors[offline_sensor.layer].push_back(hotsensor_intt);
0199     }
0200     if(Verbosity()>2)
0201     {
0202       std::cout << "layer " << offline_sensor.layer << " ladder phi " << offline_sensor.ladder_phi << " ladder z " << offline_sensor.ladder_z << " strip y " << offline_sensor.strip_y << " strip x " << offline_sensor.strip_x << std::endl;
0203     }
0204   }
0205 
0206   if(Verbosity()>2) std::cout << "TPC hot sensors:" << std::endl;
0207   for(int i=0; i<n_tpchot; i++)
0208   {
0209     int layer = tpc_hottree->GetIntValue(i,"layer");
0210     MissingSensor hotsensor_tpc = {
0211       .layer = layer,
0212       .phi_id = tpc_hottree->GetIntValue(i,"sector"),
0213       .z_id = tpc_hottree->GetIntValue(i,"side"),
0214       .local_x = tpc_hottree->GetIntValue(i,"pad")
0215     };
0216     if(layer>=7&&layer<=54)
0217     {
0218       _dead_sensors[layer].push_back(hotsensor_tpc);
0219     }
0220     if(Verbosity()>2)
0221     {
0222       std::cout << "layer " << layer << " sector " << hotsensor_tpc.phi_id << " side " << hotsensor_tpc.z_id << " phibin " << hotsensor_tpc.local_x << std::endl;
0223     }
0224   }
0225 
0226   if(Verbosity()>2) std::cout << "TPC dead sensors:" << std::endl;
0227   for(int i=0; i<n_tpcdead; i++)
0228   {
0229     int layer = tpc_deadtree->GetIntValue(i,"layer");
0230     MissingSensor deadsensor_tpc = {
0231       .layer = layer,
0232       .phi_id = tpc_deadtree->GetIntValue(i,"sector"),
0233       .z_id = tpc_deadtree->GetIntValue(i,"side"),
0234       .local_x = tpc_deadtree->GetIntValue(i,"pad")
0235     };
0236     if(layer>=7&&layer<=54)
0237     {
0238       _dead_sensors[layer].push_back(deadsensor_tpc);
0239     }
0240     if(Verbosity()>2)
0241     {
0242       std::cout << "layer " << layer << " sector " << deadsensor_tpc.phi_id << " side " << deadsensor_tpc.z_id << " phibin " << deadsensor_tpc.local_x << std::endl;
0243     }
0244   }
0245   // TODO: figure out how TPOT works for this
0246 }
0247 
0248 void TrackDeadLayerTagging::get_surfaces()
0249 {
0250   for(int layer=0; layer<57; layer++)
0251   {
0252     _surface_maps.emplace_back();
0253     if(layer<3) // mvtx
0254     {
0255       const int n_staves = 12 + 4*layer;
0256       const int n_chips = 9;
0257       for(int stave_id=0;stave_id<n_staves;stave_id++)
0258       {
0259         for(int chip_id=0;chip_id<n_chips;chip_id++)
0260         {
0261           // we only need spatial information from surface, and ActsSurfaceMaps::getSiliconSurface sets time info to 0 anyway
0262           TrkrDefs::hitsetkey hskey = MvtxDefs::genHitSetKey(layer,stave_id,chip_id,0);
0263           //std::cout << "layer " << layer << ", stave " << stave_id << ", chip " << chip_id << ", hskey " << std::to_string(hskey) << std::endl;
0264           
0265           const auto& surface = _tGeometry->maps().getSiliconSurface(hskey);
0266           if(!surface)
0267           {
0268             std::cout << PHWHERE << " cannot find MVTX surface with layer " << layer << " stave " << stave_id << " chip " << chip_id << std::endl;
0269           }
0270           else
0271           {
0272             //std::cout << "geoId: " << surface->geometryId() << std::endl;
0273             _surface_maps[layer].insert(std::make_pair(surface->geometryId(),hskey));
0274           }
0275         }
0276       }
0277     }
0278     else if(layer<7) // intt
0279     {
0280       int n_ladders;
0281       if(layer<5) n_ladders = 12;
0282       else n_ladders = 16;
0283       const int n_zsegments = 4;
0284       for(int ladder_phi_id=0;ladder_phi_id<n_ladders;ladder_phi_id++)
0285       {
0286         for(int ladder_z_id=0;ladder_z_id<n_zsegments;ladder_z_id++)
0287         {
0288           // we only need spatial information from surface, and ActsSurfaceMaps::getSiliconSurface sets time info to 0 anyway
0289           TrkrDefs::hitsetkey hskey = InttDefs::genHitSetKey(layer,ladder_z_id,ladder_phi_id,0);
0290           //std::cout << "hskey " << std::to_string(hskey) << std::endl;
0291           //std::cout << "gen hitsetkey should have layer " << layer << ", z id " << ladder_z_id << ", phi id " << ladder_phi_id << std::endl;
0292           //std::cout << "actually has layer " << int(TrkrDefs::getLayer(hskey)) << ", z id " << int(InttDefs::getLadderZId(hskey)) << ", phi id " << int(InttDefs::getLadderPhiId(hskey)) << std::endl;
0293 
0294           const auto& surface = _tGeometry->maps().getSiliconSurface(hskey);
0295           if(!surface)
0296           {
0297             std::cout << PHWHERE << " cannot find INTT surface with layer " << layer << " ladder_z " << ladder_z_id << " ladder_phi " << ladder_phi_id << std::endl;
0298           }
0299           else
0300           {
0301             //std::cout << "geoId: " << surface->geometryId() << std::endl;
0302             _surface_maps[layer].insert(std::make_pair(surface->geometryId(),hskey));
0303           }
0304         }
0305       }
0306     }
0307     else if(layer<55) // tpc
0308     {
0309       // TPC organizes itself into "fake" subsurfaces, 120 per layer, unrelated with sector positions(?)
0310       // This means that this surface map doesn't actually do anything useful...
0311       // Keeping it anyway for now, in case there turns out to be a way to grab sector and side from subsurface key
0312       const int n_subsurf = 288;
0313       // We need a dummy hitsetkey that tells the surface map which layer we're on
0314       TrkrDefs::hitsetkey hskey = TpcDefs::genHitSetKey(layer,0,0);
0315       for(int subsurf_id = 0; subsurf_id < n_subsurf; subsurf_id++)
0316       {
0317         // this only uses the layer from the hitsetkey, so this is fine
0318         const auto& surface = _tGeometry->maps().getTpcSurface(hskey,subsurf_id);
0319         if(!surface)
0320         {
0321           std::cout << PHWHERE << " cannot find TPC surface with layer " << layer << " subsurface " << subsurf_id << std::endl;
0322         }
0323         else
0324         {
0325           //std::cout << "geoId: " << surface->geometryId() << std::endl;
0326           // we store the subsurface ID instead of the hitsetkey here because most of the information of the hitsetkey is unknown (and the layer is redundant with the index)
0327           _surface_maps[layer].insert(std::make_pair(surface->geometryId(),subsurf_id));
0328         }
0329       }
0330     }
0331     else if(layer<57) // tpot
0332     {
0333       const int n_tiles = 8;
0334       for(int tile_id = 0; tile_id < n_tiles; tile_id++)
0335       {
0336         TrkrDefs::hitsetkey hskey_phi = MicromegasDefs::genHitSetKey(layer,MicromegasDefs::SegmentationType::SEGMENTATION_PHI,tile_id);
0337         TrkrDefs::hitsetkey hskey_z = MicromegasDefs::genHitSetKey(layer,MicromegasDefs::SegmentationType::SEGMENTATION_Z,tile_id);
0338 
0339         const auto& surface_phi = _tGeometry->maps().getMMSurface(hskey_phi);
0340         const auto& surface_z = _tGeometry->maps().getMMSurface(hskey_z);
0341 
0342         if(!surface_phi)
0343         {
0344           std::cout << PHWHERE << " cannot find TPOT surface with layer " << layer << " phi segmentation, tile " << tile_id << std::endl;
0345         }
0346         else
0347         {
0348           //std::cout << "geoId: " << surface_phi->geometryId() << std::endl;
0349           _surface_maps[layer].insert(std::make_pair(surface_phi->geometryId(),hskey_phi));
0350         }
0351         if(!surface_z)
0352         {
0353           std::cout << PHWHERE << " cannot find TPOT surface with layer " << layer << " z segmentation, tile " << tile_id << std::endl;
0354         }
0355         else
0356         {
0357           //std::cout << "geoId: " << surface_z->geometryId() << std::endl;
0358           _surface_maps[layer].insert(std::make_pair(surface_z->geometryId(),hskey_z));
0359         }
0360       }
0361     }
0362   }
0363 }
0364 
0365 std::vector<float> TrackDeadLayerTagging::getFitParams(TrackSeed* seed)
0366 {
0367   std::vector<float> fitpars;
0368   fitpars.push_back(1./fabs(seed->get_qOverR()));
0369   fitpars.push_back(seed->get_X0());
0370   fitpars.push_back(seed->get_Y0());
0371   fitpars.push_back(seed->get_slope());
0372   fitpars.push_back(seed->get_Z0());
0373   return fitpars;
0374 }
0375 
0376 std::vector<MissingSpacePoint> TrackDeadLayerTagging::get_missing_space_points(SvtxTrack* track, std::vector<int> missing_layers)
0377 {
0378   std::vector<MissingSpacePoint> missing_space_points;
0379   if(useActsPropagator)
0380   {
0381     ActsPropagator propagator(_tGeometry);
0382     propagator.verbosity(0);
0383     ActsPropagator::BoundTrackParamResult result_trackparams = propagator.makeTrackParams(track,_vertexmap);
0384     if(result_trackparams.ok())
0385     {
0386       Acts::BoundTrackParameters trackparams = result_trackparams.value();
0387       for(int layer : missing_layers)
0388       {
0389         ActsPropagator::BTPPairResult result_propagated_params = propagator.propagateTrack(trackparams,layer);
0390         if(result_propagated_params.ok())
0391         {
0392           // [path length (mm), BoundTrackParams]
0393           Acts::BoundTrackParameters updated_params = result_propagated_params.value().second;
0394           //Surface surf = updated_params.referenceSurface().getSharedPtr();
0395 
0396           Acts::Vector3 global_intersection = updated_params.position(_tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
0397           Acts::Vector3 tangent = updated_params.direction();
0398 
0399           int nCompatibleSurfaces = 0;
0400           for(auto& [geoId, hskey] : _surface_maps[layer])
0401           {
0402             Surface surface = _tGeometry->geometry().tGeometry->findSurface(geoId)->getSharedPtr();
0403             auto local_result = surface->globalToLocal(_tGeometry->geometry().getGeoContext(),global_intersection * Acts::UnitConstants::cm,tangent);
0404             if(local_result.ok())
0405             {
0406               Acts::Vector2 local_intersection = local_result.value() / Acts::UnitConstants::cm;
0407               // trackparams = updated_params;
0408               if(surface->insideBounds(local_intersection * Acts::UnitConstants::cm))
0409               {
0410                 MissingSpacePoint missingpoint = { layer, surface, local_intersection, global_intersection };
0411                 missing_space_points.push_back(missingpoint);
0412                 nCompatibleSurfaces++;
0413               }
0414             }
0415           }
0416           if(nCompatibleSurfaces>1 && Verbosity()>0) std::cout << PHWHERE << "NOTE: More than one compatible surface found!" << std::endl;
0417         }
0418         else
0419         {
0420           if(Verbosity()>0) std::cout << PHWHERE << "ActsPropagator returned error " << result_propagated_params.error().message() << " on propagating to layer " << layer << std::endl;
0421         }
0422       }
0423     }
0424     else
0425     {
0426       if(Verbosity()>0) std::cout << PHWHERE << "ActsPropagator returned error " << result_trackparams.error().message() << " on formation of track params" << std::endl;
0427     }
0428   }
0429   else
0430   {
0431     TrackSeed* siliconseed = track->get_silicon_seed();
0432     TrackSeed* tpcseed = track->get_tpc_seed();
0433     
0434     // get fit parameters for track; use TPC seed for this unless silicon-only track
0435     std::vector<float> fitpars = tpcseed? getFitParams(tpcseed) : getFitParams(siliconseed);
0436 
0437     // get nearest cluster phi for each layer (necessary to disambiguate multiple compatible surfaces)
0438     std::vector<std::pair<float,float>> filledlayer_pairs;
0439     if(siliconseed)
0440     {
0441       for(auto si_ckey = siliconseed->begin_cluster_keys(); si_ckey != siliconseed->end_cluster_keys(); si_ckey++)
0442       {
0443         TrkrCluster* c = _clustermap->findCluster(*si_ckey);
0444         Acts::Vector3 globalpos = _tGeometry->getGlobalPosition(*si_ckey,c);
0445         float phi = atan2(globalpos.y(),globalpos.x());
0446         filledlayer_pairs.push_back(std::make_pair(TrkrDefs::getLayer(*si_ckey),phi));
0447       }
0448     }
0449     if(tpcseed)
0450     {
0451       for(auto tpc_ckey = tpcseed->begin_cluster_keys(); tpc_ckey != tpcseed->end_cluster_keys(); tpc_ckey++)
0452       {
0453         TrkrCluster* c = _clustermap->findCluster(*tpc_ckey);
0454         Acts::Vector3 globalpos = tpcglobal.getGlobalPositionDistortionCorrected(*tpc_ckey,c,track->get_crossing());
0455         float phi = atan2(globalpos.y(),globalpos.x());
0456         filledlayer_pairs.push_back(std::make_pair(TrkrDefs::getLayer(*tpc_ckey),phi));
0457       }
0458     }
0459 
0460     // for each missing layer, find helix intersection with compatible surface in that layer
0461     for(int missing_layer : missing_layers)
0462     {
0463       if(Verbosity()>1)
0464       {
0465         std::cout << "layer " << missing_layer << ":" << std::endl;
0466       }
0467       std::vector<MissingSpacePoint> compatible_surfaces;
0468       
0469       for(auto& [geoId, hskey] : _surface_maps[missing_layer])
0470       {
0471         Surface surface = _tGeometry->geometry().tGeometry->findSurface(geoId)->getSharedPtr();
0472         Acts::Vector3 surface_center = surface->center(_tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
0473         Acts::Vector3 intersection = TrackFitUtils::get_helix_surface_intersection(surface,fitpars,surface_center,_tGeometry);
0474         Acts::Vector3 tangent = TrackFitUtils::get_helix_tangent(fitpars,intersection).second;
0475         if(Verbosity()>2)
0476         {
0477           std::cout << "surface center: (" << surface_center.x() << ", " << surface_center.y() << ", " << surface_center.z() << ")" << std::endl;
0478           std::cout << "intersection with surface: (" << intersection.x() << ", " << intersection.y() << ", " << intersection.z() << ")" << std::endl;
0479           std::cout << "tangent vector: (" << tangent.x() << ", " << tangent.y() << ", " << tangent.z() << ")" << std::endl;
0480         }
0481         auto intersection_local_result = surface->globalToLocal(_tGeometry->geometry().getGeoContext(),intersection * Acts::UnitConstants::cm,tangent * Acts::UnitConstants::cm);
0482         if(!intersection_local_result.ok())
0483         {
0484           if(Verbosity()>2)
0485           {
0486             std::cout << "local coordinate conversion failed..." << std::endl;
0487           }
0488           continue;
0489         }
0490         Acts::Vector2 intersection_local = intersection_local_result.value() / Acts::UnitConstants::cm;
0491         if(Verbosity()>2)
0492         {
0493           std::cout << "local intersection: (" << intersection_local.x() << ", " << intersection_local.y() << ")" << std::endl;
0494         }
0495         if(surface->insideBounds(intersection_local * Acts::UnitConstants::cm))
0496         {
0497           if(Verbosity()>1)
0498           {
0499             std::cout << "missing track state within bounds" << std::endl;
0500           }
0501           MissingSpacePoint missingpoint = { missing_layer, surface, intersection_local, intersection };
0502           compatible_surfaces.push_back(missingpoint);
0503         }
0504       }
0505 
0506       // disambiguate multiple compatible surfaces on helix
0507       if(compatible_surfaces.size()>1)
0508       {
0509         // find cluster phi of closest on-track layer to space point
0510         float closest_phi = 1e9;
0511         // cover edge cases
0512         if(missing_layer < filledlayer_pairs[0].first) closest_phi = filledlayer_pairs[0].second;
0513         else if(missing_layer > filledlayer_pairs.back().first) closest_phi = filledlayer_pairs.back().second;
0514         for(size_t i=0; i<filledlayer_pairs.size(); i++)
0515         {
0516           if(filledlayer_pairs[i].first < missing_layer) continue;
0517           closest_phi = filledlayer_pairs[i].second;
0518           int abovelayer_difference = (filledlayer_pairs[i].first - missing_layer);
0519           int belowlayer_difference = filledlayer_pairs[i-1].first - missing_layer;
0520           if(abovelayer_difference < belowlayer_difference) closest_phi = filledlayer_pairs[i].second;
0521           else closest_phi = filledlayer_pairs[i-1].second;
0522           break;
0523         }
0524         // determine which space point has most similar phi to closest on-track cluster phi
0525         int mostcompatible_index = -1;
0526         float best_dphi = 1e9;
0527         for(size_t i=0; i < compatible_surfaces.size(); i++)
0528         {
0529           Acts::Vector3 globalpos = compatible_surfaces[i].global_intersection;
0530           float phi = atan2(globalpos.y(),globalpos.x());
0531           if(fabs(phi-closest_phi)<best_dphi)
0532           {
0533             mostcompatible_index = i;
0534             best_dphi = fabs(phi-closest_phi);
0535           }
0536         }
0537         missing_space_points.push_back(compatible_surfaces[mostcompatible_index]);
0538       }
0539       else if(compatible_surfaces.size()==1)
0540       {
0541         missing_space_points.push_back(compatible_surfaces[0]);
0542       }
0543     }
0544   }
0545   return missing_space_points;
0546 }
0547 
0548 std::vector<int> TrackDeadLayerTagging::get_missing_layers(SvtxTrack* track)
0549 {
0550   std::vector<TrkrDefs::cluskey> ckeys;
0551   if(track->get_silicon_seed()) std::copy(track->get_silicon_seed()->begin_cluster_keys(),track->get_silicon_seed()->end_cluster_keys(),std::back_inserter(ckeys));
0552   if(track->get_tpc_seed()) std::copy(track->get_tpc_seed()->begin_cluster_keys(),track->get_tpc_seed()->end_cluster_keys(),std::back_inserter(ckeys));
0553 
0554   int start_layer = 0;
0555   int end_layer = 56;
0556   if(ignoreSilicon) start_layer = 7;
0557   if(ignoreTpc) end_layer = 6;
0558 
0559   std::vector<int> missing_layers;
0560   for(int l=start_layer; l<=end_layer; l++)
0561   {
0562     if(std::find_if(ckeys.begin(),ckeys.end(),[&](TrkrDefs::cluskey ckey){ return (TrkrDefs::getLayer(ckey) == l); }) == ckeys.end())
0563     {
0564       missing_layers.push_back(l);
0565     }
0566   }
0567 
0568   if(missing_layers.size()>0)
0569   {
0570     if(Verbosity()>2)
0571     {
0572       std::cout << "missing layers found: ";
0573       for(auto& l : missing_layers) std::cout << l << ", ";
0574       std::cout << std::endl;
0575     }
0576   }
0577   return missing_layers;
0578 }
0579 
0580 std::vector<MissingSensor> TrackDeadLayerTagging::find_sensors(std::vector<MissingSpacePoint> &spacepoints)
0581 {
0582   std::vector<MissingSensor> missing_sensors;
0583   for(auto& spacepoint : spacepoints)
0584   {
0585     if(Verbosity()>1) std::cout << "sensor layer: " << spacepoint.layer << std::endl;
0586     if(spacepoint.layer<3) // MVTX
0587     {
0588       int row;
0589       int col;
0590       SegmentationAlpide::localToDetector(spacepoint.local_intersection.x(),spacepoint.local_intersection.y(),row,col);
0591       TrkrDefs::hitsetkey hskey = _surface_maps[spacepoint.layer][spacepoint.surface->geometryId()];
0592       int phi_id = MvtxDefs::getStaveId(hskey);
0593       int z_id = MvtxDefs::getChipId(hskey); 
0594       if(Verbosity()>1)
0595       {
0596         std::cout << "stave: " << phi_id << ", chip: " << z_id << std::endl;
0597         std::cout << "pixel: " << row << ", " << col << std::endl;
0598       }
0599       if(row >= 0 && row <= 511 && col>=0 && col<=1023)
0600       {
0601         MissingSensor missing_sensor = { _event, _trackID, spacepoint.layer, phi_id, z_id, row, col };
0602         missing_sensors.push_back(missing_sensor);
0603       }
0604     }
0605     else if(spacepoint.layer<7) // INTT
0606     {
0607       int strip_y;
0608       int strip_z;
0609       CylinderGeomIntt* intt_cylinder = dynamic_cast<CylinderGeomIntt*>(_inttGeom->GetLayerGeom(spacepoint.layer));
0610       TrkrDefs::hitsetkey hskey = _surface_maps[spacepoint.layer][spacepoint.surface->geometryId()];
0611       int phi_id = InttDefs::getLadderPhiId(hskey);
0612       int z_id = InttDefs::getLadderZId(hskey);
0613       int layer = TrkrDefs::getLayer(hskey);
0614       intt_cylinder->find_strip_index_values(z_id,spacepoint.local_intersection.x(),spacepoint.local_intersection.y(),strip_y,strip_z);
0615       if(Verbosity()>1)
0616       {
0617         std::cout << "layer: " << layer << std::endl;
0618         std::cout << "ladder: " << phi_id << ", " << z_id << std::endl;
0619         std::cout << "strip: " << strip_y << ", " << strip_z << std::endl;
0620       }
0621       if(strip_y >= 0 && strip_y <= 255 && strip_z>=0 && ((z_id == 0 && strip_z<=7) || (z_id == 1 && strip_z<=4)))
0622       {
0623         MissingSensor missing_sensor = { _event, _trackID, spacepoint.layer, phi_id, z_id, strip_y, strip_z };
0624         missing_sensors.push_back(missing_sensor);
0625       }
0626     }
0627     else if(spacepoint.layer<55) // TPC
0628     {
0629       int side = (spacepoint.global_intersection.z()>0)? 1 : 0;
0630       double phi = atan2(spacepoint.global_intersection.y(),spacepoint.global_intersection.x());
0631       PHG4TpcGeom* layergeom = _tpcGeom->GetLayerCellGeom(spacepoint.layer);
0632       int phibin = layergeom->get_phibin(phi,side);
0633       int zbin = layergeom->get_zbin(spacepoint.global_intersection.z());
0634 
0635       // turns out we don't need the _surface_maps for TPC because the sector and side information isn't carried there anyway
0636       const int sector_step = layergeom->get_phibins() / 12;
0637       int sector_id = phibin / sector_step;
0638 
0639       if(Verbosity()>1)
0640       {
0641         std::cout << "layer: " << spacepoint.layer << std::endl;
0642         std::cout << "sector: " << sector_id << ", side " << side << std::endl;
0643         std::cout << "phibin: " << phibin << ", zbin: " << zbin << std::endl;
0644       }
0645       MissingSensor missing_sensor = { _event, _trackID, spacepoint.layer, sector_id, side, phibin, zbin };
0646       missing_sensors.push_back(missing_sensor);
0647     }
0648     else if(spacepoint.layer<57) // TPOT
0649     {
0650       int strip_z = 0; // no strip z index for TPOT, that info is in tile ID
0651       CylinderGeomMicromegas* tpotlayergeom = dynamic_cast<CylinderGeomMicromegas*>(_tpotGeom->GetLayerGeom(spacepoint.layer));
0652       
0653       TrkrDefs::hitsetkey hskey = _surface_maps[spacepoint.layer][spacepoint.surface->geometryId()];
0654       int tile_id = MicromegasDefs::getTileId(hskey);
0655       int seg_id = (MicromegasDefs::getSegmentationType(hskey) == MicromegasDefs::SegmentationType::SEGMENTATION_Z)? 0 : 1;
0656       TVector2 local_tvector(spacepoint.local_intersection.x(),spacepoint.local_intersection.y());
0657 
0658       int strip_id = tpotlayergeom->find_strip_from_local_coords(tile_id,_tGeometry,local_tvector);
0659 
0660       if(Verbosity()>1)
0661       {
0662         std::cout << "layer: " << spacepoint.layer << std::endl;
0663         std::cout << "tile ID: " << tile_id << ", segmentation " << seg_id << std::endl;
0664         std::cout << "strip_id " << strip_id << std::endl;
0665       }
0666       MissingSensor missing_sensor = { _event, _trackID, spacepoint.layer, tile_id, seg_id, strip_id, strip_z };
0667       missing_sensors.push_back(missing_sensor);
0668     }
0669   }
0670   return missing_sensors;
0671 }
0672 
0673 bool TrackDeadLayerTagging::is_dead(MissingSensor& sensor)
0674 {
0675   std::vector<MissingSensor> dead_thislayer = _dead_sensors.at(sensor.layer);
0676   return (std::find_if(dead_thislayer.begin(),dead_thislayer.end(),
0677             [&sensor](const MissingSensor& candidate) { return sensor.isSameSensor(candidate); })
0678           != dead_thislayer.end());
0679 }
0680 
0681 //_____________________________________
0682 void TrackDeadLayerTagging::find_dead_layers(PHCompositeNode *topNode)
0683 {
0684 
0685   for(const auto &[key, track] : *_trackmap)
0686   {
0687     if(!track) continue;
0688 
0689     _trackID = track->get_id();
0690     if(Verbosity()>0) std::cout << "track ID " << _trackID << std::endl;
0691     if(std::isnan(track->get_x()) ||
0692        std::isnan(track->get_y()) ||
0693        std::isnan(track->get_z()) ||
0694        std::isnan(track->get_px()) ||
0695        std::isnan(track->get_py()) ||
0696        std::isnan(track->get_pz()))
0697     {
0698       std::cout << "malformed track:" << std::endl;
0699       track->identify();
0700       std::cout << "skipping..." << std::endl;
0701       continue;
0702     }
0703 
0704     std::vector<int> missing_layers = get_missing_layers(track);
0705     std::vector<MissingSpacePoint> missing_space_points = get_missing_space_points(track,missing_layers);
0706     std::vector<MissingSensor> missing_sensors = find_sensors(missing_space_points);
0707  
0708     if(Verbosity()>0 && missing_sensors.size()>0)
0709     {
0710       std::cout << "missing sensors: " << std::endl;
0711       for(auto& sensor : missing_sensors)
0712       {
0713         if(!(ignoreSilicon && sensor.layer<7) && !(ignoreTpc && sensor.layer>=7))
0714         {
0715           std::cout << "layer: " << sensor.layer << " phi_id: " << sensor.phi_id << " z_id: " << sensor.z_id << " local_x: " << sensor.local_x << " local_y: " << sensor.local_y << std::endl;
0716         }
0717       }
0718     }
0719 
0720     for(size_t i=0; i<missing_sensors.size(); i++)
0721     {
0722       MissingSensor& sensor = missing_sensors[i];
0723       MissingSpacePoint& point = missing_space_points[i];
0724       if(!(ignoreSilicon && sensor.layer<7) && !(ignoreTpc && sensor.layer>=7))
0725       {
0726         _layer = sensor.layer;
0727         _phiID = sensor.phi_id;
0728         _zID = sensor.z_id;
0729         _localindX = sensor.local_x;
0730         _localindY = sensor.local_y;
0731         _isDead = is_dead(sensor)? 1 : 0;
0732         _global_x = point.global_intersection.x();
0733         _global_y = point.global_intersection.y();
0734         _global_z = point.global_intersection.z();
0735         _missing_sensors->Fill();
0736       }
0737     }
0738   }
0739 }
0740 
0741 //___________________________________
0742 void TrackDeadLayerTagging::GetNodes(PHCompositeNode *topNode)
0743 {
0744 
0745   _clustermap = findNode::getClass<TrkrClusterContainer>(topNode,"TRKR_CLUSTER");
0746   if(!_clustermap && _event<2)
0747   {
0748     std::cout << PHWHERE << " cannot find TrkrClusterContainer" << std::endl;
0749   }
0750 
0751   _trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0752   if(!_trackmap && _event<2)
0753   {
0754     std::cout << PHWHERE << " cannot find SvtxTrackMap" << std::endl;
0755   }
0756 
0757   _vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0758   if(!_vertexmap && _event<2)
0759   {
0760     std::cout << PHWHERE << " cannot find SvtxVertexMap" << std::endl;
0761   }
0762 
0763   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0764   if(!_tGeometry && _event<2)
0765   {
0766     std::cout << PHWHERE << " cannot find ActsGeometry" << std::endl;
0767   }
0768 
0769   _inttGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0770   if(!_inttGeom && _event<2)
0771   {
0772     std::cout << PHWHERE << " cannot find CYLINDERGEOM_INTT" << std::endl;
0773   }
0774 
0775   _tpcGeom = findNode::getClass<PHG4TpcGeomContainer>(topNode,"TPCGEOMCONTAINER");
0776   if(!_tpcGeom && _event<2)
0777   {
0778     std::cout << PHWHERE << " cannot find CYLINDERCELLGEOM_SVTX" << std::endl;
0779   }
0780 
0781   _tpotGeom = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0782   if(!_tpotGeom && _event<2)
0783   {
0784     std::cout << PHWHERE << " cannot find CYLINDERGEOM_MICROMEGAS_FULL" << std::endl;
0785   }
0786 
0787 }
0788 
0789 //______________________________________
0790 int TrackDeadLayerTagging::End(PHCompositeNode *topNode)
0791 {
0792  PHTFileServer::get().cd( _outfile );
0793  _missing_sensors->Write();
0794  return 0;
0795 }
0796
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!