packages/trackreco/PrelimDistortionCorrection.cc

0001 /*!
0002  *  \File PrelimDistortionCorrection.cc
0003  *  \brief  Makes preliminary TPC distortion corrections in pp running and replaces TPC seed parameters with new fits
0004  *  \author Tony Frawley
0005  */
0006
0007 #include "PrelimDistortionCorrection.h"
0008
0009 #include "ALICEKF.h"
0010
0011 #include <ffamodules/CDBInterface.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <g4detectors/PHG4TpcCylinderGeom.h>
0016 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0017
0018 #include <phfield/PHField.h>
0019 #include <phfield/PHFieldUtility.h>
0020 #include <phfield/PHFieldConfig.h>
0021 #include <phfield/PHFieldConfigv1.h>
0022
0023 #include <phool/PHTimer.h>
0024 #include <phool/getClass.h>
0025 #include <phool/phool.h>                       // for PHWHERE
0026
0027 // tpc distortion correction
0028 #include <tpc/TpcDistortionCorrectionContainer.h>
0029
0030 #include <trackbase/TrackFitUtils.h>
0031 #include <trackbase/TrkrCluster.h>
0032 #include <trackbase/TrkrClusterContainer.h>
0033 #include <trackbase/TrkrDefs.h>
0034 #include <trackbase/ActsGeometry.h>
0035
0036 #include <trackbase_historic/ActsTransformations.h>
0037 #include <trackbase_historic/TrackSeedContainer.h>
0038 #include <trackbase_historic/TrackSeed_v2.h>
0039 #include <trackbase_historic/TrackSeedHelper.h>
0040
0041 #include <Geant4/G4SystemOfUnits.hh>
0042
0043 #include <Eigen/Core>
0044 #include <Eigen/Dense>
0045
0046 #include <iostream>                            // for operator<<, basic_ostream
0047 #include <vector>
0048
0049 //___________________________________________________________________________________________
0050 PrelimDistortionCorrection::PrelimDistortionCorrection(const std::string& name)
0051   : SubsysReco(name)
0052 {}
0053
0054 //______________________________________________________
0055 PrelimDistortionCorrection::~PrelimDistortionCorrection()
0056 {
0057   if( m_own_fieldmap )
0058   { delete _field_map; }
0059 }
0060
0061 //___________________________________________________________________________________________
0062 int PrelimDistortionCorrection::End(PHCompositeNode* /*unused*/)
0063 {
0064   return Fun4AllReturnCodes::EVENT_OK;
0065 }
0066
0067 //___________________________________________________________________________________________
0068 int PrelimDistortionCorrection::InitRun(PHCompositeNode* topNode)
0069 {
0070
0071   int ret = get_nodes(topNode);
0072   if (ret != Fun4AllReturnCodes::EVENT_OK) { return ret; }
0073
0074   PHFieldConfigv1 fcfg;
0075   fcfg.set_field_config(PHFieldConfig::FieldConfigTypes::Field3DCartesian);
0076
0077   // load magnetic field map from CDB
0078   auto magField = CDBInterface::instance()->getUrl("FIELDMAP_TRACKING");
0079   fcfg.set_filename(magField);
0080
0081   // compare field config from that on node tree
0082   /*
0083    * if the magnetic field is already on the node tree PHFieldUtility::GetFieldConfigNode returns the existing configuration.
0084    * One must then check wheter the two configurations are identical, to decide whether one must use the field from node tree or create our own.
0085    * Otherwise the configuration passed as argument is stored on the node tree.
0086    */
0087   const auto node_fcfg = PHFieldUtility::GetFieldConfigNode(&fcfg, topNode);
0088   if( fcfg == *node_fcfg )
0089   {
0090     // both configurations are identical, use field map from node tree
0091     std::cout << "PrelimDistortionCorrection::InitRun - using field map found from node tree" << std::endl;
0092     _field_map = PHFieldUtility::GetFieldMapNode(&fcfg, topNode);
0093     m_own_fieldmap = false;
0094   } else {
0095     // both configurations differ. Use our own field map
0096     std::cout << "PrelimDistortionCorrection::InitRun - using own field map" << std::endl;
0097     _field_map = PHFieldUtility::BuildFieldMap(&fcfg);
0098     m_own_fieldmap = true;
0099   }
0100
0101   fitter = std::make_unique<ALICEKF>(topNode,_cluster_map,_field_map, _fieldDir,
0102                      _min_clusters_per_track,_max_sin_phi,Verbosity());
0103   fitter->setNeonFraction(Ne_frac);
0104   fitter->setArgonFraction(Ar_frac);
0105   fitter->setCF4Fraction(CF4_frac);
0106   fitter->setNitrogenFraction(N2_frac);
0107   fitter->setIsobutaneFraction(isobutane_frac);
0108   fitter->useConstBField(_use_const_field);
0109   fitter->setConstBField(_const_field);
0110   fitter->useFixedClusterError(_use_fixed_clus_err);
0111   fitter->setFixedClusterError(0,_fixed_clus_err.at(0));
0112   fitter->setFixedClusterError(1,_fixed_clus_err.at(1));
0113   fitter->setFixedClusterError(2,_fixed_clus_err.at(2));
0114   //  _field_map = PHFieldUtility::GetFieldMapNode(nullptr,topNode);
0115   // m_Cache = magField->makeCache(m_tGeometry->magFieldContext);
0116
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 //____________________________________________________________________________________________
0121 int PrelimDistortionCorrection::get_nodes(PHCompositeNode* topNode)
0122 {
0123   //---------------------------------
0124   // Get Objects off of the Node Tree
0125   //---------------------------------
0126
0127   _tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0128   if(!_tgeometry)
0129     {
0130       std::cout << "No Acts tracking geometry, exiting." << std::endl;
0131       return Fun4AllReturnCodes::ABORTEVENT;
0132     }
0133
0134   // tpc distortion correction
0135   // input tpc distortion correction module edge
0136   m_dcc_module_edge = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerModuleEdge");
0137   if (m_dcc_module_edge)
0138   {
0139     std::cout << "PrelimDistortionCorrection::get_nodes - found TPC distortion correction container module edge" << std::endl;
0140   }
0141
0142   // input tpc distortion correction static
0143   m_dcc_static = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerStatic");
0144   if (m_dcc_static)
0145   {
0146     std::cout << "PrelimDistortionCorrection::get_nodes - found TPC distortion correction container static" << std::endl;
0147   }
0148
0149   // input tpc distortion correction average
0150   m_dcc_average = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerAverage");
0151   if (m_dcc_average)
0152   {
0153     std::cout << "PrelimDistortionCorrection::get_nodes - found TPC distortion correction container average" << std::endl;
0154   }
0155
0156
0157   if(_use_truth_clusters) {
0158     _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER_TRUTH");
0159   } else {
0160     _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0161 }
0162
0163   if (!_cluster_map)
0164   {
0165     std::cerr << PHWHERE << "PrelimDistortionCorrection::get_nodes - ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0166     return Fun4AllReturnCodes::ABORTEVENT;
0167   }
0168
0169   _track_map = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0170   if (!_track_map)
0171   {
0172     std::cerr << PHWHERE << "PrelimDistortionCorrection::get_nodes - ERROR: Can't find TrackSeedContainer " << std::endl;
0173     return Fun4AllReturnCodes::ABORTEVENT;
0174   }
0175
0176   PHG4TpcCylinderGeomContainer *geom_container =
0177       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0178   if (!geom_container)
0179   {
0180     std::cerr << PHWHERE << "PrelimDistortionCorrection::get_nodes - ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0181     return Fun4AllReturnCodes::ABORTRUN;
0182   }
0183
0184   for(int i=7;i<=54;i++)
0185   {
0186     radii.push_back(geom_container->GetLayerCellGeom(i)->get_radius());
0187   }
0188
0189   return Fun4AllReturnCodes::EVENT_OK;
0190 }
0191
0192 //________________________________________________________________________________________________________________
0193 int PrelimDistortionCorrection::process_event(PHCompositeNode* /*topNode*/)
0194 {
0195   if(!_pp_mode)
0196   {
0197     std::cout << "PrelimDistortionCorrection called but not in pp_mode, do nothing!" << std::endl;
0198     return Fun4AllReturnCodes::EVENT_OK;
0199   }
0200
0201   PHTimer timer("PrelimDistortionCorrectionTimer");
0202   timer.stop();
0203   timer.restart();
0204
0205   if(Verbosity()>0)
0206   {
0207     std::cout << "starting PrelimDistortionCorrection process_event" << std::endl;
0208   }
0209
0210   // These accumulate trackseeds as we loop over tracks, keylist is a vector of vectors of seed cluskeys
0211   // The seeds are all given to the fitter at once
0212   std::vector<std::vector<TrkrDefs::cluskey>> keylist_A;
0213   PositionMap correctedOffsetTrackClusPositions;   //  this is an std::map<TrkrDefs::cluskey, Acts::Vector3>
0214
0215   for(unsigned int track_it = 0; track_it != _track_map->size(); ++track_it )
0216   {
0217     if(Verbosity()>0) { std::cout << "TPC seed " << track_it << std::endl; }
0218     // if not a TPC track, ignore
0219     TrackSeed* track = _track_map->get(track_it);
0220     if(!track)
0221       {
0222     continue;
0223       }
0224
0225     if(Verbosity() > 0)
0226       {
0227     std::cout << "Input seed pars for " << track_it
0228           << " q " << track->get_charge()
0229           << " qOverR " << fabs(track->get_qOverR()) * track->get_charge()
0230           << " X0 " << TrackSeedHelper::get_x(track)
0231           << " Y0 " << TrackSeedHelper::get_y(track)
0232           << " Z0 " << TrackSeedHelper::get_z(track)
0233           << " eta " << track->get_eta()
0234           << " phi " << track->get_phi()
0235           << std::endl;
0236       }
0237
0238     const bool is_tpc = std::any_of(
0239       track->begin_cluster_keys(),
0240       track->end_cluster_keys(),
0241       []( const TrkrDefs::cluskey& key ) { return TrkrDefs::getTrkrId(key) == TrkrDefs::tpcId; } );
0242
0243     if(is_tpc)
0244     {
0245       // We want to move all clusters in this seed to point to Z0 = 0
0246       float Z0 = track->get_Z0();
0247       float offset_Z = 0.0 - Z0;
0248
0249       if(Verbosity() > 0)
0250     {
0251       std::cout << "  processing seed with offset_Z " << offset_Z
0252             << " eta " << track->get_eta()
0253         << " x " << TrackSeedHelper::get_x(track)
0254         << " y " << TrackSeedHelper::get_y(track)
0255         << " z " << TrackSeedHelper::get_z(track)
0256             << std::endl;
0257     }
0258
0259       // We want to make  distortion corrections to all clusters in this seed after offsetting the z values
0260       std::vector<TrkrDefs::cluskey> dumvec;
0261       for(TrackSeed::ConstClusterKeyIter iter = track->begin_cluster_keys();
0262       iter != track->end_cluster_keys();
0263       ++iter)
0264     {
0265       TrkrDefs::cluskey cluskey = *iter;
0266       TrkrCluster *cluster = _cluster_map->findCluster(cluskey);
0267
0268     // get global position
0269     const Acts::Vector3 pos = _tgeometry->getGlobalPosition(cluskey, cluster);
0270
0271     // apply z offset
0272     Acts::Vector3 offsetpos(pos.x(), pos.y(), pos.z() + offset_Z);
0273
0274     // apply distortion corrections
0275     if (m_dcc_module_edge)
0276     {
0277       offsetpos = m_distortionCorrection.get_corrected_position(offsetpos, m_dcc_module_edge);
0278     }
0279
0280     if (m_dcc_static)
0281     {
0282       offsetpos = m_distortionCorrection.get_corrected_position(offsetpos, m_dcc_static);
0283     }
0284
0285     if (m_dcc_average)
0286     {
0287       offsetpos = m_distortionCorrection.get_corrected_position(offsetpos, m_dcc_average);
0288     }
0289
0290       // now move the distortion corrected cluster back by offset_Z, to preserve the z measurement info
0291     const Acts::Vector3 corrpos(offsetpos(0), offsetpos(1), offsetpos(2) - offset_Z);
0292       correctedOffsetTrackClusPositions.emplace(cluskey,corrpos);
0293       dumvec.push_back(cluskey);
0294
0295       if(Verbosity() > 0)
0296         {
0297           std::cout << " cluskey " << cluskey << " input pos " << pos(0) << "  " << pos(1) << "  " << pos(2)
0298             << "   corr. pos " << corrpos(0) << "  " << corrpos(1) << "  " << corrpos(2) << std::endl
0299             << "distortion correction " <<  corrpos(0) - pos(0) << "  " << corrpos(1) - pos(1) << "  " << corrpos(2) - pos(2)
0300             << std::endl;
0301
0302         }
0303     }
0304
0305       /// Can't circle fit a seed with less than 3 clusters, skip it
0306       if(dumvec.size() < 3)
0307     { continue; }
0308       keylist_A.push_back(dumvec);
0309
0310       if(Verbosity() > 0)
0311     {
0312       std::cout << "Added  input seed " << track_it << "  becomes output seed " << keylist_A.size() - 1 << std::endl;
0313     }
0314
0315     } // end if TPC seed
0316
0317   }  // end loop over tracks
0318
0319   // reset the seed map for the TPC
0320   _track_map->Reset();
0321
0322   // refit the corrected clusters
0323   std::vector<float> trackChi2;
0324   auto seeds = fitter->ALICEKalmanFilter(keylist_A, true, correctedOffsetTrackClusPositions, trackChi2);
0325
0326   // update the seed parameters on the node tree
0327   // This calls circle fit and line fit (setting x, y, z, phi, eta) and sets qOverR explicitly using q from KF
0328   publishSeeds(seeds.first, correctedOffsetTrackClusPositions);
0329
0330   return Fun4AllReturnCodes::EVENT_OK;
0331 }
0332
0333 //____________________________________________________________________________________________________________
0334 void PrelimDistortionCorrection::publishSeeds(std::vector<TrackSeed_v2>& seeds, PositionMap& positions)
0335 {
0336   int seed_index = 0;
0337   for(auto& seed: seeds )
0338   {
0339     if(seed.size_cluster_keys() < 3)
0340       {
0341     continue;   // ALICEKalmanFilter can drop clusters. Seeds require at least 3 clusters for circle fit
0342       }
0343
0344     /// The ALICEKF gives a better charge determination at high pT
0345     int q = seed.get_charge();
0346     TrackSeedHelper::circleFitByTaubin(&seed,positions, 7, 55);
0347     TrackSeedHelper::lineFit(&seed,positions, 7, 55);
0348
0349     seed.set_qOverR(fabs(seed.get_qOverR()) * q);
0350     seed.set_phi(TrackSeedHelper::get_phi(&seed,positions));
0351     _track_map->insert(&seed);
0352
0353     if(Verbosity() > 0)
0354       {
0355     std::cout << "Publishing seed " << seed_index
0356           << " q " << q
0357           << " qOverR " << fabs(seed.get_qOverR()) * q
0358           << " x " << TrackSeedHelper::get_x(&seed)
0359           << " y " << TrackSeedHelper::get_y(&seed)
0360           << " z " << TrackSeedHelper::get_z(&seed)
0361           << " pT " << seed.get_pt()
0362           << " eta " << seed.get_eta()
0363           << " phi " << seed.get_phi()
0364           << std::endl;
0365       }
0366     if(Verbosity() > 5)
0367       {
0368     TrackSeed* readseed = _track_map->get(seed_index);
0369     readseed->identify();
0370       }
0371
0372     seed_index++;
0373   }
0374 }