packages/trackreco/PrelimDistortionCorrectionAuAu.cc

0001 /*!
0002  *  \File PrelimDistortionCorrectionAuAu.cc
0003  *  \brief  Refits AuAu TPC seeds with distortion corrections and current calibrations
0004  *  \author Tony Frawley
0005  */
0006
0007 #include "PrelimDistortionCorrectionAuAu.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 PrelimDistortionCorrectionAuAu::PrelimDistortionCorrectionAuAu(const std::string& name)
0051   : SubsysReco(name)
0052 {}
0053
0054 //______________________________________________________
0055 PrelimDistortionCorrectionAuAu::~PrelimDistortionCorrectionAuAu()
0056 {
0057   if( m_own_fieldmap )
0058   { delete _field_map; }
0059 }
0060
0061 //___________________________________________________________________________________________
0062 int PrelimDistortionCorrectionAuAu::End(PHCompositeNode* /*unused*/)
0063 {
0064   return Fun4AllReturnCodes::EVENT_OK;
0065 }
0066
0067 //___________________________________________________________________________________________
0068 int PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 << "PrelimDistortionCorrectionAuAu::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 PrelimDistortionCorrectionAuAu::process_event(PHCompositeNode* /*topNode*/)
0194 {
0195   if(_pp_mode)
0196   {
0197     std::cout << "PrelimDistortionCorrectionAuAu called but in pp_mode, do nothing!" << std::endl;
0198     return Fun4AllReturnCodes::EVENT_OK;
0199   }
0200
0201   PHTimer timer("PrelimDistortionCorrectionAuAuTimer");
0202   timer.stop();
0203   timer.restart();
0204
0205   if(Verbosity()>0)
0206   {
0207     std::cout << "starting PrelimDistortionCorrectionAuAu 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 Z0 = 0;
0248       float offset_Z = 0.0 - Z0;
0249
0250       if(Verbosity() > 0)
0251     {
0252       std::cout << "  processing seed with offset_Z " << offset_Z
0253             << " eta " << track->get_eta()
0254         << " x " << TrackSeedHelper::get_x(track)
0255         << " y " << TrackSeedHelper::get_y(track)
0256         << " z " << TrackSeedHelper::get_z(track)
0257             << std::endl;
0258     }
0259
0260       // We want to make  distortion corrections to all clusters in this seed after offsetting the z values
0261       std::vector<TrkrDefs::cluskey> dumvec;
0262       for(TrackSeed::ConstClusterKeyIter iter = track->begin_cluster_keys();
0263       iter != track->end_cluster_keys();
0264       ++iter)
0265     {
0266       TrkrDefs::cluskey cluskey = *iter;
0267       TrkrCluster *cluster = _cluster_map->findCluster(cluskey);
0268
0269     // get global position
0270     const Acts::Vector3 pos = _tgeometry->getGlobalPosition(cluskey, cluster);
0271
0272     // apply z offset
0273     Acts::Vector3 offsetpos(pos.x(), pos.y(), pos.z() + offset_Z);
0274
0275     // apply distortion corrections
0276     if (m_dcc_module_edge)
0277     {
0278       offsetpos = m_distortionCorrection.get_corrected_position(offsetpos, m_dcc_module_edge);
0279     }
0280
0281     if (m_dcc_static)
0282     {
0283       offsetpos = m_distortionCorrection.get_corrected_position(offsetpos, m_dcc_static);
0284     }
0285
0286     if (m_dcc_average)
0287     {
0288       offsetpos = m_distortionCorrection.get_corrected_position(offsetpos, m_dcc_average);
0289     }
0290
0291       // now move the distortion corrected cluster back by offset_Z, to preserve the z measurement info
0292     const Acts::Vector3 corrpos(offsetpos(0), offsetpos(1), offsetpos(2) - offset_Z);
0293       correctedOffsetTrackClusPositions.emplace(cluskey,corrpos);
0294       dumvec.push_back(cluskey);
0295
0296       if(Verbosity() > 0)
0297         {
0298           std::cout << " cluskey " << cluskey << " input pos " << pos(0) << "  " << pos(1) << "  " << pos(2)
0299             << "   corr. pos " << corrpos(0) << "  " << corrpos(1) << "  " << corrpos(2) << std::endl
0300             << "distortion correction " <<  corrpos(0) - pos(0) << "  " << corrpos(1) - pos(1) << "  " << corrpos(2) - pos(2)
0301             << std::endl;
0302
0303         }
0304     }
0305
0306       /// Can't circle fit a seed with less than 3 clusters, skip it
0307       if(dumvec.size() < 3)
0308     { continue; }
0309       keylist_A.push_back(dumvec);
0310
0311       if(Verbosity() > 0)
0312     {
0313       std::cout << "Added  input seed " << track_it << "  becomes output seed " << keylist_A.size() - 1 << std::endl;
0314     }
0315
0316     } // end if TPC seed
0317
0318   }  // end loop over tracks
0319
0320   // reset the seed map for the TPC
0321   _track_map->Reset();
0322
0323   // refit the corrected clusters
0324   std::vector<float> trackChi2;
0325   auto seeds = fitter->ALICEKalmanFilter(keylist_A, true, correctedOffsetTrackClusPositions, trackChi2);
0326
0327   // update the seed parameters on the node tree
0328   // This calls circle fit and line fit (setting x, y, z, phi, eta) and sets qOverR explicitly using q from KF
0329   publishSeeds(seeds.first, correctedOffsetTrackClusPositions);
0330
0331   return Fun4AllReturnCodes::EVENT_OK;
0332 }
0333
0334 //____________________________________________________________________________________________________________
0335 void PrelimDistortionCorrectionAuAu::publishSeeds(std::vector<TrackSeed_v2>& seeds, PositionMap& positions)
0336 {
0337   int seed_index = 0;
0338   for(auto& seed: seeds )
0339   {
0340     if(seed.size_cluster_keys() < 3)
0341       {
0342     continue;   // ALICEKalmanFilter can drop clusters. Seeds require at least 3 clusters for circle fit
0343       }
0344
0345     /// The ALICEKF gives a better charge determination at high pT
0346     int q = seed.get_charge();
0347     TrackSeedHelper::circleFitByTaubin(&seed,positions, 7, 55);
0348     TrackSeedHelper::lineFit(&seed,positions, 7, 55);
0349
0350     seed.set_qOverR(fabs(seed.get_qOverR()) * q);
0351     seed.set_phi(TrackSeedHelper::get_phi(&seed,positions));
0352     _track_map->insert(&seed);
0353
0354     if(Verbosity() > 0)
0355       {
0356     std::cout << "Publishing seed " << seed_index
0357           << " q " << q
0358           << " qOverR " << fabs(seed.get_qOverR()) * q
0359           << " x " << TrackSeedHelper::get_x(&seed)
0360           << " y " << TrackSeedHelper::get_y(&seed)
0361           << " z " << TrackSeedHelper::get_z(&seed)
0362           << " pT " << seed.get_pt()
0363           << " eta " << seed.get_eta()
0364           << " phi " << seed.get_phi()
0365           << std::endl;
0366       }
0367     if(Verbosity() > 5)
0368       {
0369     TrackSeed* readseed = _track_map->get(seed_index);
0370     readseed->identify();
0371       }
0372
0373     seed_index++;
0374   }
0375 }