PHGenFitPkg/PHGenFit/Fitter.cc

0001 /*!
0002  *  \file       Fitter.cc
0003  *  \brief      Fitter class handles setups for the fitting.
0004  *  \details    Fitter class handles setups for the fitting like Geometry, Fields, fitter choice, etc.
0005  *  \author     Haiwang Yu <yuhw@nmsu.edu>
0006  */
0007
0008 // PHGenFit
0009 #include "Fitter.h"
0010
0011 #include "Track.h"
0012
0013 // ROOT
0014 #include <RVersion.h>  // for ROOT_VERSION, ROOT_VERSION...
0015 #include <TGeoManager.h>
0016
0017 // GenFit
0018 #include <GenFit/AbsKalmanFitter.h>
0019 #include <GenFit/DAF.h>
0020 #include <GenFit/EventDisplay.h>
0021 #include <GenFit/FieldManager.h>
0022 #include <GenFit/FitStatus.h>
0023 #include <GenFit/KalmanFitter.h>
0024 #include <GenFit/KalmanFitterRefTrack.h>
0025 #include <GenFit/MaterialEffects.h>
0026 #include <GenFit/TGeoMaterialInterface.h>
0027 #include <GenFit/Track.h>
0028
0029 // GenFitExp
0030 #include <genfitexp/Field.h>
0031
0032 #include <cassert>
0033 #include <cstddef>
0034 #include <iostream>
0035
0036 namespace genfit
0037 {
0038   class AbsTrackRep;
0039 }
0040
0041 #define LogDEBUG(exp) std::cout << "DEBUG: " << __FILE__ << ": " << __LINE__ << ": " << (exp) << std::endl
0042 #define LogERROR(exp) std::cout << "ERROR: " << __FILE__ << ": " << __LINE__ << ": " << (exp) << std::endl
0043 #define LogWARNING(exp) std::cout << "WARNING: " << __FILE__ << ": " << __LINE__ << ": " << (exp) << std::endl
0044
0045 namespace PHGenFit
0046 {
0047   Fitter::Fitter(
0048       const std::string& tgeo_file_name,
0049       const PHField* field,
0050       const std::string& fitter_choice,
0051       const std::string& /*track_rep_choice*/,
0052       const bool doEventDisplay)
0053     : verbosity(1000)
0054     , _doEventDisplay(doEventDisplay)
0055   {
0056     _tgeo_manager = new TGeoManager("Default", "Geane geometry");
0057     TGeoManager::Import(tgeo_file_name.data());
0058
0059     assert(field);
0060     genfit::Field* fieldMap = new genfit::Field(field);
0061
0062     genfit::FieldManager::getInstance()->init(fieldMap);
0063     genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
0064
0065     // init event display
0066     if (_doEventDisplay)
0067     {
0068       _display = genfit::EventDisplay::getInstance();
0069     }
0070     else
0071     {
0072       _display = nullptr;
0073     }
0074
0075     // init fitter
0076     if (fitter_choice.compare("KalmanFitterRefTrack") == 0)
0077     {
0078       _fitter = new genfit::KalmanFitterRefTrack();
0079     }
0080     else if (fitter_choice.compare("KalmanFitter") == 0)
0081 // NOLINTNEXTLINE(bugprone-branch-clone)
0082     {
0083       _fitter = new genfit::KalmanFitter();
0084     }
0085     else if (fitter_choice.compare("DafSimple") == 0)
0086     {
0087       _fitter = new genfit::DAF(false);
0088     }
0089     else if (fitter_choice.compare("DafRef") == 0)
0090     {
0091       _fitter = new genfit::DAF(true);
0092     }
0093     else
0094 // NOLINTNEXTLINE(bugprone-branch-clone)
0095     {
0096       _fitter = new genfit::KalmanFitter();
0097     }
0098
0099     genfit::Exception::quiet(true);
0100   }
0101
0102   Fitter::~Fitter()
0103   {
0104     delete _fitter;
0105     // delete _tgeo_manager;
0106     //_tgeo_manager->Delete();
0107     delete _display;
0108   }
0109
0110   int Fitter::processTrack(PHGenFit::Track* track, const bool save_to_evt_disp)
0111   {
0112     genfit::Track* fitTrack = track->getGenFitTrack();
0113
0114 #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 00, 0)
0115     try
0116     {
0117       fitTrack->checkConsistency();
0118     }
0119     catch (genfit::Exception& e)
0120     {
0121       if (verbosity >= 2)
0122       {
0123         std::cerr << "genfit::Track::checkConsistency() failed!" << std::endl;
0124         std::cerr << e.what();
0125       }
0126       return -1;
0127     }
0128 #else
0129     if (!fitTrack->checkConsistency())
0130     {
0131       if (verbosity >= 2) LogWARNING("genfit::Track::checkConsistency() failed!");
0132       return -1;
0133     }
0134 #endif
0135     try
0136     {
0137       _fitter->processTrack(fitTrack);
0138     }
0139     catch (genfit::Exception& e)
0140     {
0141       if (verbosity >= 1)
0142       {
0143         std::cerr << "PHGenFit::Fitter::processTrack::Exception: \n";
0144         std::cerr << e.what();
0145         std::cerr << "Exception, next track" << std::endl;
0146       }
0147       return -1;
0148     }
0149 #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 00, 0)
0150     try
0151     {
0152       fitTrack->checkConsistency();
0153     }
0154     catch (genfit::Exception& e)
0155     {
0156       if (verbosity >= 2)
0157       {
0158         std::cerr << "genfit::Track::checkConsistency() failed!" << std::endl;
0159         std::cerr << e.what();
0160       }
0161       return -1;
0162     }
0163 #else
0164
0165     if (!fitTrack->checkConsistency())
0166     {
0167       if (verbosity >= 2) LogWARNING("genfit::Track::checkConsistency() failed!");
0168       return -1;
0169     }
0170 #endif
0171     genfit::AbsTrackRep* rep = fitTrack->getCardinalRep();
0172     if (!fitTrack->getFitStatus(rep)->isFitConverged())
0173     {
0174       if (verbosity >= 2)
0175       {
0176         LogWARNING("Track could not be fitted successfully! Fit is not converged!");
0177       }
0178       return -1;
0179     }
0180
0181     if (_display and save_to_evt_disp)
0182     {
0183       _display->addEvent(track->getGenFitTrack());
0184     }
0185
0186     return 0;
0187   }
0188
0189   Fitter* Fitter::getInstance(const std::string& tgeo_file_name,
0190                               const PHField* field,
0191                               const std::string& fitter_choice, const std::string& track_rep_choice,
0192                               const bool doEventDisplay)
0193   {
0194     TGeoManager* tgeo_manager = TGeoManager::Import(tgeo_file_name.data(), "Default");
0195     if (!tgeo_manager)
0196     {
0197       LogERROR("No TGeoManager found!");
0198       return nullptr;
0199     }
0200
0201     assert(field);
0202     genfit::Field* fieldMap = new genfit::Field(field);
0203     return new Fitter(tgeo_manager, fieldMap, fitter_choice, track_rep_choice, doEventDisplay);
0204   }
0205
0206   Fitter::Fitter(TGeoManager* tgeo_manager, genfit::AbsBField* fieldMap,
0207                  const PHGenFit::Fitter::FitterType& fitter_choice,
0208                  const PHGenFit::Fitter::TrackRepType& /*track_rep_choice*/,
0209                  const bool doEventDisplay)
0210     : verbosity(0)
0211     , _tgeo_manager(tgeo_manager)
0212     , _doEventDisplay(doEventDisplay)
0213   {
0214     genfit::FieldManager::getInstance()->init(
0215         fieldMap);
0216     genfit::MaterialEffects::getInstance()->init(
0217         new genfit::TGeoMaterialInterface());
0218
0219     // init event display
0220     if (_doEventDisplay)
0221     {
0222       _display = genfit::EventDisplay::getInstance();
0223     }
0224     else
0225     {
0226       _display = nullptr;
0227     }
0228
0229     // init fitter
0230     if (fitter_choice == PHGenFit::Fitter::KalmanFitter)
0231     {
0232       _fitter = new genfit::KalmanFitter();
0233     }
0234     else if (fitter_choice == PHGenFit::Fitter::KalmanFitterRefTrack)
0235     {
0236       _fitter = new genfit::KalmanFitterRefTrack();
0237     }
0238     if (fitter_choice == PHGenFit::Fitter::DafSimple)
0239     {
0240       _fitter = new genfit::DAF(false);
0241     }
0242     else if (fitter_choice == PHGenFit::Fitter::DafRef)
0243     {
0244       _fitter = new genfit::DAF(true);
0245     }
0246     else
0247     {
0248       _fitter = nullptr;
0249       LogERROR("This fitter not implemented!");
0250     }
0251   }
0252
0253   Fitter* Fitter::getInstance(TGeoManager* tgeo_manager,
0254                               const PHField* field,
0255                               const std::string& fitter_choice, const std::string& track_rep_choice,
0256                               const bool doEventDisplay)
0257   {
0258     if (!tgeo_manager)
0259     {
0260       LogERROR("No TGeoManager found!");
0261       return nullptr;
0262     }
0263
0264     assert(field);
0265     genfit::Field* fieldMap = new genfit::Field(field);
0266     return new Fitter(tgeo_manager, fieldMap, fitter_choice, track_rep_choice, doEventDisplay);
0267   }
0268
0269   Fitter::Fitter(TGeoManager* tgeo_manager, genfit::AbsBField* fieldMap,
0270                  const std::string& fitter_choice, const std::string& /*track_rep_choice*/,
0271                  const bool doEventDisplay)
0272     : verbosity(0)
0273     , _tgeo_manager(tgeo_manager)
0274     , _doEventDisplay(doEventDisplay)
0275   {
0276     genfit::FieldManager::getInstance()->init(
0277         fieldMap);
0278     genfit::MaterialEffects::getInstance()->init(
0279         new genfit::TGeoMaterialInterface());
0280
0281     // init event display
0282     if (_doEventDisplay)
0283     {
0284       _display = genfit::EventDisplay::getInstance();
0285     }
0286     else
0287     {
0288       _display = nullptr;
0289     }
0290
0291     // init fitter
0292     if (fitter_choice.compare("KalmanFitterRefTrack") == 0)
0293     {
0294       _fitter = new genfit::KalmanFitterRefTrack();
0295     }
0296     else if (fitter_choice.compare("KalmanFitter") == 0)
0297     {
0298       _fitter = new genfit::KalmanFitter();
0299     }
0300     else if (fitter_choice.compare("DafSimple") == 0)
0301     {
0302       _fitter = new genfit::DAF(false);
0303     }
0304     else if (fitter_choice.compare("DafRef") == 0)
0305     {
0306       _fitter = new genfit::DAF(true);
0307     }
0308     else
0309     {
0310       _fitter = nullptr;
0311       LogERROR("This fitter not implemented!");
0312     }
0313   }
0314
0315   int Fitter::displayEvent()
0316   {
0317     if (_display)
0318     {
0319       _display->open();
0320     }
0321     else if (verbosity >= 0)
0322     {
0323       LogERROR("No genfit::EventDisplay found!");
0324     }
0325
0326     return 0;
0327   }
0328
0329   Fitter* Fitter::getInstance(TGeoManager* tgeo_manager,
0330                               const PHField* field,
0331                               const PHGenFit::Fitter::FitterType& fitter_choice,
0332                               const PHGenFit::Fitter::TrackRepType& track_rep_choice,
0333                               const bool doEventDisplay)
0334   {
0335     if (!tgeo_manager)
0336     {
0337       LogERROR("No TGeoManager found!");
0338       return nullptr;
0339     }
0340
0341     assert(field);
0342     genfit::Field* fieldMap = new genfit::Field(field);
0343
0344     return new Fitter(tgeo_manager, fieldMap, fitter_choice, track_rep_choice, doEventDisplay);
0345   }
0346
0347 }  // namespace PHGenFit