diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml index a070d543de61ff6e44f8952b24be492cea927c9f..bfd81e73dea5721de9e4b9bbb16f0d25512b6d3b 100644 --- a/Detector/DetDriftChamber/compact/det.xml +++ b/Detector/DetDriftChamber/compact/det.xml @@ -1,5 +1,7 @@ <?xml version="1.0" encoding="UTF-8"?> -<lccdd> +<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> <info name="DriftChamber" title="Test with Drift Chamber" @@ -15,6 +17,7 @@ <gdmlFile ref="materials.xml"/> </includes> + <define> <constant name="world_size" value="2226*mm"/> <constant name="world_x" value="world_size"/> @@ -22,6 +25,7 @@ <constant name="world_z" value="world_size"/> <!-- SDT --> + <constant name="DetID_DC" value="7"/> <constant name="SDT_radius_min" value="234*mm"/> <constant name="SDT_radius_max" value="1720*mm"/> <constant name="SDT_half_length" value="2225*mm"/> diff --git a/Examples/options/sim_fit_CRD.py b/Examples/options/sim_fit_CRD.py deleted file mode 100644 index e877aeb4f5eef6c197cc7a939baba95e76ea46e0..0000000000000000000000000000000000000000 --- a/Examples/options/sim_fit_CRD.py +++ /dev/null @@ -1,278 +0,0 @@ -#!/usr/bin/env python -from Gaudi.Configuration import * - -############################################################################## -# Event service -############################################################################## -from Configurables import k4DataSvc -dsvc = k4DataSvc("EventDataSvc") - -############################################################################## -# Random seed -############################################################################## -from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_ -# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi -rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4 -rndmengine.SetSingleton = True -rndmengine.Seeds = [10] - -from Configurables import MarlinEvtSeeder -evtseeder = MarlinEvtSeeder("EventSeeder") - -############################################################################## -# Detector geometry -############################################################################## -geometry_option = "CRD_o1_v01/CRD_o1_v01_noECAL.xml" - -if not os.getenv("DETCRDROOT"): - print("Can't find the geometry. Please setup envvar DETCRDROOT." ) - sys.exit(-1) - -geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option) -if not os.path.exists(geometry_path): - print("Can't find the compact geometry file: %s"%geometry_path) - sys.exit(-1) - -from Configurables import GeomSvc -geosvc = GeomSvc("GeomSvc") -geosvc.compact = geometry_path - -from Configurables import GearSvc -gearsvc = GearSvc("GearSvc") - -############################################################################## -# Physics Generator -############################################################################## -from Configurables import GenAlgo -from Configurables import GtGunTool -from Configurables import StdHepRdr -from Configurables import SLCIORdr -from Configurables import HepMCRdr -from Configurables import GenPrinter -gun = GtGunTool("GtGunTool") -gun.Particles = ["mu-"] -gun.EnergyMins = [100.] # GeV -gun.EnergyMaxs = [100.] # GeV -gun.ThetaMins = [0] # deg -gun.ThetaMaxs = [180] # deg -gun.PhiMins = [0] # deg -gun.PhiMaxs = [360] # deg -# stdheprdr = StdHepRdr("StdHepRdr") -# stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep" -# lciordr = SLCIORdr("SLCIORdr") -# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio" -# hepmcrdr = HepMCRdr("HepMCRdr") -# hepmcrdr.Input = "example_UsingIterators.txt" - -genprinter = GenPrinter("GenPrinter") - -genalg = GenAlgo("GenAlgo") -genalg.GenTools = ["GtGunTool"] -#genalg.GenTools = ["StdHepRdr"] -# genalg.GenTools = ["StdHepRdr", "GenPrinter"] -# genalg.GenTools = ["SLCIORdr", "GenPrinter"] -# genalg.GenTools = ["HepMCRdr", "GenPrinter"] - -############################################################################## -# Detector Simulation -############################################################################## -from Configurables import DetSimSvc -detsimsvc = DetSimSvc("DetSimSvc") - -from Configurables import DetSimAlg -detsimalg = DetSimAlg("DetSimAlg") -#detsimalg.VisMacs = ["vis.mac"] -detsimalg.RunCmds = [ -# "/tracking/verbose 1", -] -detsimalg.AnaElems = [ - # example_anatool.name() - # "ExampleAnaElemTool" - "Edm4hepWriterAnaElemTool" -] -detsimalg.RootDetElem = "WorldDetElemTool" - -from Configurables import AnExampleDetElemTool -example_dettool = AnExampleDetElemTool("AnExampleDetElemTool") - -from Configurables import CalorimeterSensDetTool -from Configurables import DriftChamberSensDetTool - -calo_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool") -driftchamber_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool") - -############################################################################## -# DedxSimTool -############################################################################## -# dedxoption = "DummyDedxSimTool" -dedxoption = "BetheBlochEquationDedxSimTool" - -driftchamber_sensdettool.DedxSimTool = dedxoption - -from Configurables import DummyDedxSimTool -from Configurables import BetheBlochEquationDedxSimTool - -if dedxoption == "DummyDedxSimTool": - dedx_simtool = DummyDedxSimTool("DummyDedxSimTool") -elif dedxoption == "BetheBlochEquationDedxSimTool": - dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool") - dedx_simtool.material_Z = 2 - dedx_simtool.material_A = 4 - dedx_simtool.scale = 10 - dedx_simtool.resolution = 0.0001 - -############################################################################## -# DCHDigiAlg -############################################################################## -from Configurables import DCHDigiAlg -dCHDigiAlg = DCHDigiAlg("DCHDigiAlg") -dCHDigiAlg.readout = "DriftChamberHitsCollection" -dCHDigiAlg.drift_velocity = 40#um/ns -dCHDigiAlg.mom_threshold = 0 #GeV -dCHDigiAlg.WriteAna = True - -############################################################################## -# PODIO -############################################################################## -from Configurables import PodioOutput -out = PodioOutput("outputalg") -out.filename = "CRD-o1-v01-sim00.root" -out.outputCommands = ["keep *"] - - -############################################################################## -# Silicon digitization -############################################################################## -vxdhitname = "VXDTrackerHits" -sithitname = "SITTrackerHits" -sitspname = "SITSpacePoints" -tpchitname = "TPCTrackerHits" -sethitname = "SETTrackerHits" -setspname = "SETSpacePoints" -ftdspname = "FTDSpacePoints" -ftdhitname = "FTDTrackerHits" -from Configurables import PlanarDigiAlg -digiVXD = PlanarDigiAlg("VXDDigi") -digiVXD.SimTrackHitCollection = "VXDCollection" -digiVXD.TrackerHitCollection = vxdhitname -digiVXD.ResolutionU = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004] -digiVXD.ResolutionV = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004] - -digiSIT = PlanarDigiAlg("SITDigi") -digiSIT.IsStrip = 1 -digiSIT.SimTrackHitCollection = "SITCollection" -digiSIT.TrackerHitCollection = sithitname -digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation" -digiSIT.ResolutionU = [0.007] -digiSIT.ResolutionV = [0.000] - -digiSET = PlanarDigiAlg("SETDigi") -digiSET.IsStrip = 1 -digiSET.SimTrackHitCollection = "SETCollection" -digiSET.TrackerHitCollection = sethitname -digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation" -digiSET.ResolutionU = [0.007] -digiSET.ResolutionV = [0.000] - -digiFTD = PlanarDigiAlg("FTDDigi") -digiFTD.SimTrackHitCollection = "FTDCollection" -digiFTD.TrackerHitCollection = ftdhitname -digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" -digiFTD.ResolutionU = [0.003, 0.003, 0.007, 0.007, 0.007, 0.007, 0.007, 0.007] -digiFTD.ResolutionV = [0.003, 0.003, 0, 0, 0, 0, 0, 0 ] -#digiFTD.OutputLevel = DEBUG - -from Configurables import SpacePointBuilderAlg -spSIT = SpacePointBuilderAlg("SITBuilder") -spSIT.TrackerHitCollection = sithitname -spSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation" -spSIT.SpacePointCollection = sitspname -spSIT.SpacePointAssociationCollection = "SITSpacePointAssociation" -#spSIT.OutputLevel = DEBUG - -spFTD = SpacePointBuilderAlg("FTDBuilder") -spFTD.TrackerHitCollection = ftdhitname -spFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" -spFTD.SpacePointCollection = ftdspname -spFTD.SpacePointAssociationCollection = "FTDSpacePointAssociation" -#spFTD.OutputLevel = DEBUG - - -############################################################################## -# Silicon reconstruction -############################################################################## -from Configurables import TrackSystemSvc -tracksystemsvc = TrackSystemSvc("TrackSystemSvc") - -from Configurables import SiliconTrackingAlg -tracking = SiliconTrackingAlg("SiliconTracking") -tracking.HeaderCol = "EventHeader" -tracking.VTXHitCollection = vxdhitname -tracking.SITHitCollection = sitspname -tracking.FTDPixelHitCollection = ftdhitname -tracking.FTDSpacePointCollection = ftdspname -tracking.SITRawHitCollection = sithitname -tracking.FTDRawHitCollection = ftdhitname -tracking.UseSIT = 1 -tracking.SmoothOn = 0 -#tracking.OutputLevel = DEBUG - -from Configurables import ForwardTrackingAlg -forward = ForwardTrackingAlg("ForwardTracking") -forward.FTDPixelHitCollection = ftdhitname -forward.FTDSpacePointCollection = ftdspname -forward.FTDRawHitCollection = ftdhitname -forward.Chi2ProbCut = 0.0 -forward.HitsPerTrackMin = 3 -forward.BestSubsetFinder = "SubsetSimple" -forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation", - "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"] -forward.CriteriaMin = [0, 0.9, 0, 0.995, 0, 0.8, 0, 20, 1.002, 0.1, 0, 0.99, 0, 0.999, 0, 0.99, 0] -forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08, 99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05] -#forward.OutputLevel = DEBUG - -from Configurables import TrackSubsetAlg -subset = TrackSubsetAlg("TrackSubset") -subset.TrackInputCollections = ["ForwardTracks", "SiTracks"] -subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, sitspname, ftdspname] -subset.TrackSubsetCollection = "SubsetTracks" -#subset.OutputLevel = DEBUG - -############################################################################## -# TruthTrackerAlg -############################################################################## -from Configurables import TruthTrackerAlg -truthTrackerAlg = TruthTrackerAlg("TruthTrackerAlg") -truthTrackerAlg.SiSubsetTrackCollection = "SubsetTracks" - -############################################################################## -# RecGenfitAlgSDT -############################################################################## -from Configurables import RecGenfitAlgSDT -recGenfitAlgSDT = RecGenfitAlgSDT("RecGenfitAlgSDT") -recGenfitAlgSDT.debug=10 - -############################################################################## -# NTupleSvc -############################################################################## -from Configurables import NTupleSvc -ntsvc = NTupleSvc("NTupleSvc") -ntsvc.Output = [ -#"MyTuples DATAFILE='DCH_digi_ana.root' OPT='NEW' TYP='ROOT'", - "RecGenfitAlgSDT DATAFILE='fit_SDT.root' OPT='NEW' TYP='ROOT'"] - - -############################################################################## -# ApplicationMgr -############################################################################## -from Configurables import ApplicationMgr -ApplicationMgr( - TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSIT, - spFTD, tracking, forward, subset, dCHDigiAlg, truthTrackerAlg, - recGenfitAlgSDT, out], - EvtSel = 'NONE', - EvtMax = 1, - ExtSvc = [rndmengine, dsvc, ntsvc, evtseeder, geosvc, gearsvc, tracksystemsvc], - HistogramPersistency = "ROOT", - OutputLevel=DEBUG -) diff --git a/Examples/options/tut_detsim_digi_fit_DC.py b/Examples/options/tut_detsim_digi_fit_DC.py index 52a7c026596dcee1047590b74bf08c513b1966a3..770bf4267256f540decd991517a4afb37634e57d 100644 --- a/Examples/options/tut_detsim_digi_fit_DC.py +++ b/Examples/options/tut_detsim_digi_fit_DC.py @@ -70,8 +70,8 @@ gun.Particles = ["e-"] # gun.PositionZs = [0.] # mm -gun.EnergyMins = [10.] # GeV -gun.EnergyMaxs = [10.] # GeV +gun.EnergyMins = [100.] # GeV +gun.EnergyMaxs = [100.] # GeV gun.ThetaMins = [90] # rad; 45deg gun.ThetaMaxs = [90] # rad; 45deg @@ -195,5 +195,5 @@ ApplicationMgr( TopAlg = [genalg, detsimalg, dCHDigiAlg, truthTrackerAlg, EvtMax = 10, ExtSvc = [rndmengine, dsvc, geosvc, ntsvc], HistogramPersistency = "ROOT", - OutputLevel=DEBUG + OutputLevel=ERROR ) diff --git a/Reconstruction/RecGenfitAlg/CMakeLists.txt b/Reconstruction/RecGenfitAlg/CMakeLists.txt index 630086dce7a22cef058903d2859aaed2d469959c..866156036919ed186102e28a6abfbd02938588fa 100644 --- a/Reconstruction/RecGenfitAlg/CMakeLists.txt +++ b/Reconstruction/RecGenfitAlg/CMakeLists.txt @@ -2,13 +2,11 @@ if (GenFit_FOUND) gaudi_add_module(RecGenfitAlg - SOURCES src/RecGenfitAlgSDT.cpp - src/RecGenfitAlgDC.cpp + SOURCES src/RecGenfitAlgDC.cpp src/GenfitTrack.cpp src/GenfitField.cpp src/GenfitFitter.cpp src/GenfitMaterialInterface.cpp - src/GenfitMsg.cpp LINK GearSvc Gaudi::GaudiAlgLib Gaudi::GaudiKernel diff --git a/Reconstruction/RecGenfitAlg/README.md b/Reconstruction/RecGenfitAlg/README.md index c9cddb4c8b0b9285c1a42a72140c953b9f79797a..0870fb5dd8766a3201b19c3af73949907ce0c483 100644 --- a/Reconstruction/RecGenfitAlg/README.md +++ b/Reconstruction/RecGenfitAlg/README.md @@ -7,7 +7,6 @@ It's incudingg the interface to genfit fitter, track and measurement: * GenfitTrack: class for create and access of genfit Track and Measurement * GenfitField: implementation of genfit AbsField * GenfitGeoMaterialInterface: implementation of genfit AbsMaterialInterface -* GenfitMsg: a class to get Gaudi message instance in Genfit classes. * GenfitAlgDC: a class to fit drift chamber * GenfitAlgSDT: a class to fit drift chamber+silicon diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp index 341e52f3d9c1a0e6539b62bec6225410dbbad7d4..744084254143d4480d2cb096970944aa2b6f0a58 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp +++ b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp @@ -1,5 +1,4 @@ #include "GenfitField.h" -#include "GenfitMsg.h" //External #include "DD4hep/DD4hepUnits.h" diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp index 15d40830f9d3fa8deb71e0ca3e306ee7e5ad7c37..4850cf4d360e4ed71c705a32d5a18e4cdc8f319c 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp +++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp @@ -1,6 +1,5 @@ #include "GenfitFitter.h" #include "GenfitTrack.h" -#include "GenfitMsg.h" #include "GenfitField.h" #include "GenfitMaterialInterface.h" @@ -92,8 +91,8 @@ int GenfitFitter::init(bool deleteOldFitter) { if(deleteOldFitter && m_absKalman) delete m_absKalman; - GenfitMsg::get()<<MSG::DEBUG<<"Initialize GenfitFitter with " - <<m_fitterType<<endmsg; + if(m_debug>=2)std::cout<<"Initialize GenfitFitter with " + <<m_fitterType<<std::endl; if (m_fitterType=="DAFRef") { m_absKalman = new genfit::DAF(true,getDeltaPval(), @@ -111,11 +110,11 @@ int GenfitFitter::init(bool deleteOldFitter) } else { m_absKalman = nullptr; - GenfitMsg::get()<< MSG::DEBUG<<"Fitter type is invalid:" - <<m_fitterType<<endmsg; + if(m_debug>=2)std::cout<<"Fitter type is invalid:" + <<m_fitterType<<std::endl; return -1; } - GenfitMsg::get()<<MSG::DEBUG<<"Fitter type is "<<m_fitterType<<endmsg; + if(m_debug>=2)std::cout<<"Fitter type is "<<m_fitterType<<std::endl; m_absKalman->setDebugLvl(m_debug); return 0; @@ -124,11 +123,11 @@ int GenfitFitter::init(bool deleteOldFitter) /// Fit a track from a candidate track int GenfitFitter::processTrackWithRep(GenfitTrack* track,int repID,bool resort) { - GenfitMsg::get()<<MSG::DEBUG<< "In ProcessTrackWithRep rep "<<repID<<endmsg; + if(m_debug>=2)std::cout<< "In ProcessTrackWithRep rep "<<repID<<std::endl; if(getDebug()>2) print(""); if(getDebug()>0){ - GenfitMsg::get()<<MSG::DEBUG<<"Print track seed "<<endmsg; + if(m_debug>=2)std::cout<<"Print track seed "<<std::endl; track->getTrack()->getStateSeed().Print(); } /// Do the fitting @@ -136,29 +135,29 @@ int GenfitFitter::processTrackWithRep(GenfitTrack* track,int repID,bool resort) m_absKalman->processTrackWithRep(track->getTrack(), track->getRep(repID), resort); }catch(genfit::Exception& e){ - GenfitMsg::get()<<MSG::DEBUG<<"Genfit exception caught "<<endmsg; + if(m_debug>=2)std::cout<<"Genfit exception caught "<<std::endl; e.what(); return false; } - GenfitMsg::get()<<MSG::DEBUG<<"End of ProcessTrackWithRep"<<endmsg; + if(m_debug>=2)std::cout<<"End of ProcessTrackWithRep"<<std::endl; return true; } // End of ProcessTrackWithRep /// Fit a track from a candidate track int GenfitFitter::processTrack(GenfitTrack* track, bool resort) { - GenfitMsg::get()<<MSG::DEBUG<<"In ProcessTrack"<<endmsg; + if(m_debug>=2)std::cout<<"In ProcessTrack"<<std::endl; if(getDebug()>2) print(""); /// Do the fitting try{ m_absKalman->processTrack(track->getTrack(),resort); }catch(genfit::Exception& e){ - GenfitMsg::get()<<MSG::DEBUG<<"Genfit exception caught "<<endmsg; + if(m_debug>=2)std::cout<<"Genfit exception caught "<<std::endl; e.what(); return false; } - GenfitMsg::get()<<MSG::DEBUG<<"End of ProcessTrack"<<endmsg; + if(m_debug>=2)std::cout<<"End of ProcessTrack"<<std::endl; return true; } // End of ProcessTrack @@ -167,7 +166,7 @@ void GenfitFitter::setFitterType(const char* val) { std::string oldFitterType=m_fitterType; m_fitterType = val; - GenfitMsg::get()<<MSG::DEBUG<<"Fitter type is "<<m_fitterType<<endmsg; + if(m_debug>=2)std::cout<<"Fitter type is "<<m_fitterType<<std::endl; init(oldFitterType==val); } @@ -210,8 +209,8 @@ GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom, // get track rep genfit::AbsTrackRep* rep = track->getRep(repID); if(nullptr == rep) { - GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToCylinder rep is null" - <<endmsg; + if(m_debug>=2)std::cout<<"In ExtrapolateToCylinder rep is null" + <<std::endl; return trackLength*dd4hep::cm; } @@ -219,8 +218,8 @@ GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom, genfit::TrackPoint* tp = track->getTrack()->getPointWithFitterInfo(hitID,rep); if(nullptr == tp) { - GenfitMsg::get()<<MSG::DEBUG<< - "In ExtrapolateToCylinder TrackPoint is null"<<endmsg; + if(m_debug>=2)std::cout<< + "In ExtrapolateToCylinder TrackPoint is null"<<std::endl; return trackLength*dd4hep::cm; } @@ -230,8 +229,8 @@ GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom, tp->getFitterInfo(rep))->getBackwardUpdate(); if(nullptr == state){ - GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToCylinder " - <<"no KalmanFittedStateOnPlane in backwardUpdate"<<endmsg; + if(m_debug>=2)std::cout<<"In ExtrapolateToCylinder " + <<"no KalmanFittedStateOnPlane in backwardUpdate"<<std::endl; return trackLength*dd4hep::cm; } rep->setPosMom(*state, pos*(1/dd4hep::cm), mom*(1/dd4hep::GeV)); @@ -245,9 +244,9 @@ GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom, pos = pos*dd4hep::cm; mom = mom*dd4hep::GeV; } catch(genfit::Exception& e){ - GenfitMsg::get() << MSG::ERROR + if(m_debug>=3)std::cout <<"Exception in GenfitFitter::extrapolateToCylinder " - << e.what()<<endmsg; + << e.what()<<std::endl; trackLength = 1e9*dd4hep::cm; } return trackLength*dd4hep::cm; @@ -266,8 +265,8 @@ double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom, // get track rep genfit::AbsTrackRep* rep = track->getRep(repID); if(nullptr == rep) { - GenfitMsg::get()<<MSG::DEBUG<<"In ExtrapolateToPoint rep " - <<repID<<" not exist!"<<endmsg; + if(m_debug>=2)std::cout<<"In ExtrapolateToPoint rep " + <<repID<<" not exist!"<<std::endl; return trackLength*dd4hep::cm; } @@ -278,8 +277,8 @@ double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom, genfit::TrackPoint* tp = track->getTrack()->getPointWithFitterInfo(0,rep); if(nullptr == tp) { - GenfitMsg::get()<<MSG::DEBUG<< - "In ExtrapolateToPoint TrackPoint is null"<<endmsg; + if(m_debug>=2)std::cout<< + "In ExtrapolateToPoint TrackPoint is null"<<std::endl; return trackLength*dd4hep::cm; } @@ -289,8 +288,8 @@ double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom, tp->getFitterInfo(rep))->getBackwardUpdate(); if(nullptr == state) { - GenfitMsg::get()<<MSG::DEBUG<< - "In ExtrapolateToPoint KalmanFittedStateOnPlane is null"<<endmsg; + if(m_debug>=2)std::cout<< + "In ExtrapolateToPoint KalmanFittedStateOnPlane is null"<<std::endl; return trackLength*dd4hep::cm; } trackLength = rep->extrapolateToPoint(*state, @@ -299,9 +298,9 @@ double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom, pos = pos*dd4hep::cm; mom = mom*dd4hep::GeV; } catch(genfit::Exception& e){ - GenfitMsg::get() << MSG::ERROR + if(m_debug>=3)std::cout <<"Exception in GenfitFitter::extrapolateToPoint" - << e.what()<<endmsg; + << e.what()<<std::endl; trackLength = 1e9*dd4hep::cm; } return trackLength*dd4hep::cm; @@ -343,41 +342,41 @@ void GenfitFitter::print(const char* name) if(0==strlen(name)) name=m_name; //TODO print all fitting parameters - GenfitMsg::get()<<MSG::DEBUG<<name - <<" GenfitFitter Global Parameters:"<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" Fitter type = " << m_fitterType<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" Fitting Parameters:"<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" MinIteration = " << m_minIterations<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" MaxIteration = " << m_maxIterations<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_deltaPval = " << m_deltaPval<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" Material Effect Parameters:"<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_noEffects = " << m_noEffects<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_energyLossBetheBloch= " << m_energyLossBetheBloch<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_noiseBetheBloch = " << m_noiseBetheBloch<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_noiseCoulomb = " << m_noiseCoulomb<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_energyLossBrems = " << m_energyLossBrems<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_noiseBrems = " << m_noiseBrems<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name + if(m_debug>=2)std::cout<<name + <<" GenfitFitter Global Parameters:"<<std::endl; + if(m_debug>=2)std::cout<<name + <<" Fitter type = " << m_fitterType<<std::endl; + if(m_debug>=2)std::cout<<name + <<" Fitting Parameters:"<<std::endl; + if(m_debug>=2)std::cout<<name + <<" MinIteration = " << m_minIterations<<std::endl; + if(m_debug>=2)std::cout<<name + <<" MaxIteration = " << m_maxIterations<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_deltaPval = " << m_deltaPval<<std::endl; + if(m_debug>=2)std::cout<<name + <<" Material Effect Parameters:"<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_noEffects = " << m_noEffects<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_energyLossBetheBloch= " << m_energyLossBetheBloch<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_noiseBetheBloch = " << m_noiseBetheBloch<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_noiseCoulomb = " << m_noiseCoulomb<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_energyLossBrems = " << m_energyLossBrems<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_noiseBrems = " << m_noiseBrems<<std::endl; + if(m_debug>=2)std::cout<<name <<" m_ignoreBoundariesBetweenEqualMaterials= " - << m_ignoreBoundariesBetweenEqualMaterials<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_mscModelName = " << m_mscModelName<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_debug = " << m_debug<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<name - <<" m_hist = " << m_hist<<endmsg; + << m_ignoreBoundariesBetweenEqualMaterials<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_mscModelName = " << m_mscModelName<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_debug = " << m_debug<<std::endl; + if(m_debug>=2)std::cout<<name + <<" m_hist = " << m_hist<<std::endl; if(m_fitterType=="DAF"||m_fitterType=="DAFRef"){ genfit::DAF* daf = getDAF(); @@ -388,27 +387,27 @@ void GenfitFitter::print(const char* name) for (unsigned int i=0; i<betas.size(); ++i) { sBetas<<betas.at(i)<<","; } - GenfitMsg::get()<<MSG::DEBUG<<" print "<<betas.size()<< - " betas: %s "<<sBetas.str()<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<m_name - << " m_deltaWeight = " << m_deltaWeight<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<" DAF maxIterations= " - << daf->getMaxIterations()<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<" DAF minIterations= " - << daf->getMinIterations()<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<" DAF deltaPval= " - << daf->getDeltaPval()<<endmsg; + if(m_debug>=2)std::cout<<" print "<<betas.size()<< + " betas: %s "<<sBetas.str()<<std::endl; + if(m_debug>=2)std::cout<<m_name + << " m_deltaWeight = " << m_deltaWeight<<std::endl; + if(m_debug>=2)std::cout<<" DAF maxIterations= " + << daf->getMaxIterations()<<std::endl; + if(m_debug>=2)std::cout<<" DAF minIterations= " + << daf->getMinIterations()<<std::endl; + if(m_debug>=2)std::cout<<" DAF deltaPval= " + << daf->getDeltaPval()<<std::endl; } } genfit::KalmanFitterRefTrack* ref = getKalRef(); if(nullptr != ref){ std::ostringstream sBetas; - GenfitMsg::get()<<MSG::DEBUG<<" DAF maxIterations= " - << ref->getMaxIterations()<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<" DAF minIterations= " - << ref->getMinIterations()<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<" DAF deltaPval= " - << ref->getDeltaPval()<<endmsg; + if(m_debug>=2)std::cout<<" DAF maxIterations= " + << ref->getMaxIterations()<<std::endl; + if(m_debug>=2)std::cout<<" DAF minIterations= " + << ref->getMinIterations()<<std::endl; + if(m_debug>=2)std::cout<<" DAF deltaPval= " + << ref->getDeltaPval()<<std::endl; } } @@ -422,7 +421,7 @@ void GenfitFitter::setMinIterations(unsigned int val) void GenfitFitter::setMaxIterations(unsigned int val) { if(val<=4) { - GenfitMsg::get()<<MSG::ERROR<<"Could NOT set MaxIteration<=4!"<<endmsg; + if(m_debug>=3)std::cout<<"Could NOT set MaxIteration<=4!"<<std::endl; } m_absKalman->setMaxIterations(val); m_maxIterations = val; @@ -434,10 +433,10 @@ void GenfitFitter::setMaxIterationsBetas(double bStart,double bFinal, m_absKalman->setMaxIterations(val); m_maxIterations = val; if(val<=4) { - GenfitMsg::get()<<MSG::ERROR<<"Could NOT set MaxIteration<=4!"<<endmsg; + if(m_debug>=3)std::cout<<"Could NOT set MaxIteration<=4!"<<std::endl; } - GenfitMsg::get()<<MSG::DEBUG<<"bStart "<<bStart<<" bFinal "<<bFinal - <<" MaxIteration "<<val<<endmsg; + if(m_debug>=2)std::cout<<"bStart "<<bStart<<" bFinal "<<bFinal + <<" MaxIteration "<<val<<std::endl; // genfit version less than 2.2.0 genfit::DAF* daf = dynamic_cast<genfit::DAF*> (m_absKalman); daf->setAnnealingScheme(bStart,bFinal,val); @@ -562,7 +561,7 @@ void GenfitFitter::setMaterialDebugLvl(unsigned int val) ///Set GenfitFitter parameters void GenfitFitter::setDebug(unsigned int val) { - GenfitMsg::get()<<MSG::DEBUG<<"set fitter debugLvl "<<val<<endmsg; + if(m_debug>=2)std::cout<<"set fitter debugLvl "<<val<<std::endl; m_absKalman->setDebugLvl(val); if(nullptr!=getDAF()) getDAF()->setDebugLvl(val); if(val>10) genfit::MaterialEffects::getInstance()->setDebugLvl(val); @@ -577,8 +576,8 @@ genfit::DAF* GenfitFitter::getDAF() try{ daf = dynamic_cast<genfit::DAF*> (m_absKalman); }catch(...){ - GenfitMsg::get()<< MSG::ERROR - << "dynamic_cast m_rom AbsFitter to DAF m_ailed!"<<endmsg; + if(m_debug>=3)std::cout + << "dynamic_cast m_rom AbsFitter to DAF m_ailed!"<<std::endl; return nullptr; } return daf; @@ -590,22 +589,22 @@ genfit::KalmanFitterRefTrack* GenfitFitter::getKalRef() try{ ref = dynamic_cast<genfit::KalmanFitterRefTrack*> (m_absKalman); }catch(...){ - GenfitMsg::get()<< MSG::ERROR + if(m_debug>=3)std::cout << "dynamic_cast m_rom AbsFitter to KalmanFitterRefTrack m_ailed!" - <<endmsg; + <<std::endl; } return ref; } void GenfitFitter::initHist(std::string name) { - GenfitMsg::get()<<MSG::DEBUG<<"GenfitFitter::initHist "<<name<<endmsg; + if(m_debug>=2)std::cout<<"GenfitFitter::initHist "<<name<<std::endl; //getDAF()->initHist(name); } void GenfitFitter::writeHist() { - GenfitMsg::get()<<MSG::DEBUG<<"GenfitFitter::writeHist "<<endmsg; + if(m_debug>=2)std::cout<<"GenfitFitter::writeHist "<<std::endl; //getDAF()->writeHist(); } diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp b/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp deleted file mode 100644 index 6799d68e0ed8120274c7b5a0d6ac24eb8fb4b078..0000000000000000000000000000000000000000 --- a/Reconstruction/RecGenfitAlg/src/GenfitMsg.cpp +++ /dev/null @@ -1,19 +0,0 @@ -#include "GenfitMsg.h" -#include "GaudiKernel/ServiceHandle.h" -#include "GaudiKernel/IMessageSvc.h" - -MsgStream* GenfitMsg::m_log = nullptr; - -MsgStream GenfitMsg::get(){ - if(nullptr == m_log){ - SmartIF<IMessageSvc> msgSvc( - Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc")); - m_log = new MsgStream(msgSvc, "GenfitMsg"); - (*m_log) << MSG::DEBUG << "initialize GenfitMsg" << endmsg; - } - return (*m_log); -} - -GenfitMsg::~GenfitMsg(){ - delete m_log; -} diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMsg.h b/Reconstruction/RecGenfitAlg/src/GenfitMsg.h deleted file mode 100644 index 43808b42b425ef29245445c2144c8a255facda4b..0000000000000000000000000000000000000000 --- a/Reconstruction/RecGenfitAlg/src/GenfitMsg.h +++ /dev/null @@ -1,26 +0,0 @@ -////////////////////////////////////////////////////////////////////// -/// -/// This is a class for genfit interfaces classes to access Gaudi log -/// system -/// -/// Authors: -/// Yao ZHANG(zhangyao@ihep.ac.cn) -/// -////////////////////////////////////////////////////////////////////// -#ifndef RECGENFITALG_GENFITMSG_H -#define RECGENFITALG_GENFITMSG_H - -#include "GaudiKernel/MsgStream.h" - -class GenfitMsg { - public: - static MsgStream get(); - private: - GenfitMsg(){}; - virtual ~GenfitMsg(); - static MsgStream* m_log; - -}; - -#endif - diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp index f210dc9f69c88a0cbc6407b18911c181c700accf..7b361351a17c987464727d72a02cee69859e942b 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp +++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp @@ -1,5 +1,4 @@ #include "GenfitTrack.h" -#include "GenfitMsg.h" #include "GenfitField.h" //CEPCSW @@ -105,11 +104,11 @@ bool GenfitTrack::createGenfitTrack(int pdgType,int charge, int chargeId=0; charge>0 ? chargeId=0 : chargeId=1; - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" CreateGenfitTrack seed pos(" + if(m_debug>=2)std::cout<<m_name<<" CreateGenfitTrack seed pos(" <<seedState[0]<<" "<<seedState[1]<<" "<<seedState[2]<<")cm (" <<seedState[3]<<" "<<seedState[4]<<" "<<seedState[5]<<")GeV charge " - <<charge<<" pdg "<<s_PDG[chargeId][pdgType]<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<"seedCov "<<endmsg; + <<charge<<" pdg "<<s_PDG[chargeId][pdgType]<<std::endl; + if(m_debug>=2)std::cout<<"seedCov "<<std::endl; addTrackRep(s_PDG[chargeId][pdgType]); @@ -124,10 +123,10 @@ bool GenfitTrack::createGenfitTrackFromMCParticle(int pidType, ///get track parameters from McParticle edm4hep::Vector3d mcPocaPos = mcParticle.getVertex();//mm edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV - GenfitMsg::get()<<MSG::DEBUG<<"seedPos poca "<< mcPocaPos.x - <<" "<<mcPocaPos.y<<" "<<mcPocaPos.z<<" mm "<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<"seedMom poca "<< mcPocaMom.x - <<" "<<mcPocaMom.y<<" "<<mcPocaMom.z<<" GeV "<<endmsg; + if(m_debug>=2)std::cout<<"seedPos poca "<< mcPocaPos.x + <<" "<<mcPocaPos.y<<" "<<mcPocaPos.z<<" mm "<<std::endl; + if(m_debug>=2)std::cout<<"seedMom poca "<< mcPocaMom.x + <<" "<<mcPocaMom.y<<" "<<mcPocaMom.z<<" GeV "<<std::endl; ///Pivot to first layer to avoid correction of beam pipe edm4hep::Vector3d firstLayerPos(1e9,1e9,1e9); @@ -139,24 +138,24 @@ bool GenfitTrack::createGenfitTrackFromMCParticle(int pidType, TLorentzVector seedPos(firstLayerPos.x,firstLayerPos.y,firstLayerPos.z, eventStartTime); TVector3 seedMom(firstLayerMom.x,firstLayerMom.y,firstLayerMom.z); - GenfitMsg::get()<<MSG::DEBUG<<"seedPos "<< firstLayerPos.x - <<" "<<firstLayerPos.y<<" "<<firstLayerPos.z<<endmsg; - GenfitMsg::get()<<MSG::DEBUG<<"seedMom "<< firstLayerMom.x - <<" "<<firstLayerMom.y<<" "<<firstLayerMom.z<<endmsg; + if(m_debug>=2)std::cout<<"seedPos "<< firstLayerPos.x + <<" "<<firstLayerPos.y<<" "<<firstLayerPos.z<<std::endl; + if(m_debug>=2)std::cout<<"seedMom "<< firstLayerMom.x + <<" "<<firstLayerMom.y<<" "<<firstLayerMom.z<<std::endl; ///Get error matrix of seed track TMatrixDSym covM(5);//FIXME, TODO ///Create a genfit track with seed - GenfitMsg::get()<<MSG::DEBUG<<"createGenfitTrack " ; + if(m_debug>=2)std::cout<<"createGenfitTrack " ; if(!GenfitTrack::createGenfitTrack(pidType,mcParticle.getCharge(), seedPos,seedMom,covM)){ - GenfitMsg::get()<<MSG::ERROR<<"GenfitTrack" - <<" Error in createGenfitTrack" <<endmsg; + if(m_debug>=2)std::cout<<"GenfitTrack" + <<" Error in createGenfitTrack" <<std::endl; return false; }else{ - GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack " - <<"createGenfitTrackFromMCParticle track created" <<endmsg; + if(m_debug>=2)std::cout<<"GenfitTrack " + <<"createGenfitTrackFromMCParticle track created" <<std::endl; } return true; }//end of createGenfitTrackFromMCParticle @@ -206,13 +205,13 @@ bool GenfitTrack::addSpacePointTrakerHit(edm4hep::ConstTrackerHit& hit, p[1]=pos.y*dd4hep::mm; p[2]=pos.z*dd4hep::mm; - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" addSpacePointTrakerHit"<<hitID - <<"pos "<<p[0]<<" "<<p[1]<<" "<<p[2]<<" cm"<<endmsg; + if(m_debug>=2)std::cout<<m_name<<" addSpacePointTrakerHit"<<hitID + <<"pos "<<p[0]<<" "<<p[1]<<" "<<p[2]<<" cm"<<std::endl; /// New a SpacepointMeasurement double cov[6]; for(int i=0;i<6;i++) { cov[i]=hit.getCovMatrix(i); - GenfitMsg::get()<<MSG::DEBUG<<"cov "<<cov[i]<<endmsg; + if(m_debug>=2)std::cout<<"cov "<<cov[i]<<std::endl; } TMatrixDSym hitCov(3);//FIXME? //in SimpleDigi/src/PlanarDigiAlg.cpp @@ -232,7 +231,7 @@ bool GenfitTrack::addSpacePointTrakerHit(edm4hep::ConstTrackerHit& hit, genfit::TrackPoint* trackPoint = new genfit::TrackPoint(sMeas,m_track); m_track->insertPoint(trackPoint); - GenfitMsg::get()<<MSG::DEBUG<<"end of addSpacePointTrakerHit"<<endmsg; + if(m_debug>=2)std::cout<<"end of addSpacePointTrakerHit"<<std::endl; return true; } @@ -260,10 +259,10 @@ bool GenfitTrack::addSpacePointMeasurement(const TVectorD& pos, hitCov(1,1)=sigma_t*sigma_t; hitCov(2,2)=sigma_t*sigma_t; - GenfitMsg::get()<< MSG::DEBUG<<m_name<<" addSpacePointMeasurement detID " + if(m_debug>=2)std::cout<<m_name<<" addSpacePointMeasurement detID " <<detID<<" hitId "<<hitID<<" " <<pos_t[0]<<" "<<pos_t[1]<<" "<<pos_t[2] <<" cm smeared "<<pos_smeared[0]<<" "<<pos_smeared[1]<<" " - <<pos_smeared[2]<<" sigma_t "<<sigma_t<<" cm"<<endmsg; + <<pos_smeared[2]<<" sigma_t "<<sigma_t<<" cm"<<std::endl; genfit::SpacepointMeasurement* sMeas = new genfit::SpacepointMeasurement(pos_smeared,hitCov,detID,hitID,nullptr); @@ -287,12 +286,12 @@ void GenfitTrack::addWireMeasurement(double driftDistance, wireMeas->setMaxDistance(0.5);//0.5 cm FIXME wireMeas->setLeftRightResolution(lrAmbig); - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" Add wire measurement(cm) "<<hitID + if(m_debug>=2)std::cout<<m_name<<" Add wire measurement(cm) "<<hitID <<" ep1("<<endPoint1[0]<<" "<<endPoint1[1]<<" "<<endPoint1[2] <<") ep2("<<endPoint2[0]<<" "<<endPoint2[1]<<" "<<endPoint2[2] <<") drift "<<driftDistance<<" driftErr "<<driftDistanceError <<" lr "<<lrAmbig<<" detId "<<detID << " hitId "<< hitID - <<endmsg; + <<std::endl; ///New a TrackPoint,create connection between meas. and trackPoint genfit::TrackPoint* trackPoint=new genfit::TrackPoint(wireMeas,m_track); @@ -301,8 +300,8 @@ void GenfitTrack::addWireMeasurement(double driftDistance, m_track->insertPoint(trackPoint); }catch(genfit::Exception& e){ - GenfitMsg::get()<<MSG::DEBUG<<m_name - <<"Add wire measurement exception"<<endmsg; + if(m_debug>=2)std::cout<<m_name + <<"Add wire measurement exception"<<std::endl; e.what(); } }//end of addWireMeasurementOnTrack @@ -320,9 +319,9 @@ bool GenfitTrack::addWireMeasurementOnTrack(edm4hep::Track& track,double sigma) m_gridDriftChamber->cellposition(hit.getCellID(),endPointStart, endPointEnd); int lrAmbig=0; - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" time "<<hit.getTime() + if(m_debug>=2)std::cout<<m_name<<" time "<<hit.getTime() <<" driftVelocity " <<driftVelocity<<std::endl; - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" wire pos " <<endPointStart.X() + if(m_debug>=2)std::cout<<m_name<<" wire pos " <<endPointStart.X() <<" "<<endPointStart.Y()<<" " <<endPointStart.Z()<<" " <<endPointEnd.X()<<" " <<endPointEnd.Y()<<" " <<endPointEnd.Z()<<" " <<std::endl; @@ -363,8 +362,8 @@ genfit::RKTrackRep* GenfitTrack::addTrackRep(int pdg) // genfit::MeasuredStateOnPlane stateInit(rep); // rep->setPosMomCov(stateInit, pos, mom, covM); //}catch(genfit::Exception e){ - // GenfitMsg::get() << MSG::DEBUG<<m_name<<" Exception in set track status" - // <<endmsg ; + // if(m_debug>=2)std::cout<<m_name<<" Exception in set track status" + // <<std::endl ; // std::cout<<e.what()<<std::endl; // return false; //} @@ -414,22 +413,22 @@ bool GenfitTrack::fitSuccess(int repID) const /// Test fitting converged if (!fitStatus->isFitted()||!fitStatus->isFitConverged() ||fitStatus->isFitConvergedFully()) { - GenfitMsg::get() << MSG::DEBUG<<m_name<< "Fitting is failed... isFitted" + if(m_debug>=2)std::cout<<m_name<< "Fitting is failed... isFitted" <<fitStatus->isFitted()<<" , isFitConverged " <<fitStatus->isFitConverged()<<", isFitConvergedFully " - <<fitStatus->isFitConvergedFully()<<endmsg; + <<fitStatus->isFitConvergedFully()<<std::endl; return false; } double chi2 = fitStatus->getChi2(); double ndf = fitStatus->getNdf(); - GenfitMsg::get() << MSG::INFO << "Fit result: chi2 "<<chi2 <<" ndf "<<ndf - << " chi2/ndf = " << chi2/ndf<<endmsg; + if(m_debug>=2)std::cout<< "Fit result: chi2 "<<chi2 <<" ndf "<<ndf + << " chi2/ndf = " << chi2/ndf<<std::endl; /// Test fitting chi2 if (chi2<= 0) { - GenfitMsg::get() << MSG::DEBUG<<m_name<< "Fit chi2<0 (chi2,ndf) = (" << - chi2 << "," << ndf << ")"<<endmsg; + if(m_debug>=2)std::cout<<m_name<< "Fit chi2<0 (chi2,ndf) = (" << + chi2 << "," << ndf << ")"<<std::endl; return false; } return true; @@ -446,7 +445,7 @@ void GenfitTrack::printSeed() const TLorentzVector pos = getSeedStatePos(); TVector3 mom = getSeedStateMom(); print(pos,mom); - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" NumPoints "<<getNumPoints()<<endmsg; + if(m_debug>=2)std::cout<<m_name<<" NumPoints "<<getNumPoints()<<std::endl; } void GenfitTrack::printFitted(int repID) const @@ -455,14 +454,14 @@ void GenfitTrack::printFitted(int repID) const TVector3 fittedMom; TMatrixDSym cov; - GenfitMsg::get()<<MSG::DEBUG<<m_name<< "printFitted nHit=" - <<m_track->getNumPoints()<<endmsg; + if(m_debug>=2)std::cout<<m_name<< "printFitted nHit=" + <<m_track->getNumPoints()<<std::endl; for(unsigned int iHit=0; iHit<m_track->getNumPoints(); iHit++){ if (getPosMomCovMOP((int) iHit, fittedPos, fittedMom, cov, repID)){ //print(fittedPos,fittedMom,to_string(iHit).c_str());//TODO }else{ - GenfitMsg::get()<<MSG::DEBUG<<m_name<<"Hit "<<iHit - <<" have no fitter info"<<endmsg; + if(m_debug>=2)std::cout<<m_name<<"Hit "<<iHit + <<" have no fitter info"<<std::endl; } } } @@ -475,7 +474,7 @@ void GenfitTrack::print( TLorentzVector pos, TVector3 mom, TVector3 mglo = mom; // TODO - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" "<<comment<<endmsg; + if(m_debug>=2)std::cout<<m_name<<" "<<comment<<std::endl; if(m_debug>1){ for(unsigned int i=0;i<m_reps.size();i++){ m_reps[i]->Print(); } @@ -531,8 +530,8 @@ int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom, std::cout<<" getNumPointsWithFittedInfo " <<getNumPointsWithFittedInfo(repID) <<" no TrackPoint with fitted info "<<std::endl; - GenfitMsg::get()<<MSG::DEBUG<<m_name - <<"Exception in getFittedState"<<endmsg; + if(m_debug>=2)std::cout<<m_name + <<"Exception in getFittedState"<<std::endl; std::cout<<e.what()<<std::endl; return 3; } @@ -625,8 +624,8 @@ double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir, extrapoLen = rep->extrapolateToLine(state, end0, end1, poca, pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise); }catch(genfit::Exception& e) { - GenfitMsg::get() << MSG::ERROR - <<"Exception in GenfigFitter::ExtrapolateToHit"<<e.what()<<endmsg; + if(m_debug>=3)std::cout<< + "Exception in GenfigFitter::ExtrapolateToHit"<<e.what()<<std::endl; return extrapoLen; } @@ -635,10 +634,10 @@ double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir, pocaOnWire = pocaOnWire; doca = (pocaOnWire-poca).Mag(); //TODO: debug pocaOnWire - GenfitMsg::get() << MSG::DEBUG << " poca "<<poca.x()<<","<<poca.y() - <<" "<<poca.z()<<" doca "<<doca<<endmsg; - GenfitMsg::get() << MSG::DEBUG << " pocaOnWire "<<pocaOnWire.x() - <<","<<pocaOnWire.y()<<" "<<pocaOnWire.z()<<" doca "<<doca<<endmsg; + if(m_debug>=2)std::cout<< " poca "<<poca.x()<<","<<poca.y() + <<" "<<poca.z()<<" doca "<<doca<<std::endl; + if(m_debug>=2)std::cout<< " pocaOnWire "<<pocaOnWire.x() + <<","<<pocaOnWire.y()<<" "<<pocaOnWire.z()<<" doca "<<doca<<std::endl; return extrapoLen*(dd4hep::cm); }//end of ExtrapolateToHit @@ -651,11 +650,11 @@ int GenfitTrack::addSimTrackerHits(const edm4hep::Track& track, //A TrakerHit collection std::vector<edm4hep::ConstSimTrackerHit> sortedDCTrackHitCol; - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" VXD " + if(m_debug>=2)std::cout<<m_name<<" VXD " <<lcio::ILDDetID::VXD<<" SIT " <<lcio::ILDDetID::SIT<<" SET " <<lcio::ILDDetID::SET<<" FTD " - <<lcio::ILDDetID::FTD<<" "<<endmsg; + <<lcio::ILDDetID::FTD<<" "<<std::endl; ///Get TrackerHit on Track int hitID=0; for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){ @@ -664,25 +663,25 @@ int GenfitTrack::addSimTrackerHits(const edm4hep::Track& track, UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string); encoder.setValue(hit.getCellID()); int detID=encoder[lcio::ILDCellID0::subdet]; - GenfitMsg::get()<<MSG::DEBUG<<m_name<<" "<<iHit<<" hit "<<hit - <<" detID "<<detID<<endmsg; + if(m_debug>=2)std::cout<<m_name<<" "<<iHit<<" hit "<<hit + <<" detID "<<detID<<std::endl; if(detID==lcio::ILDDetID::VXD || detID==lcio::ILDDetID::SIT || detID==lcio::ILDDetID::SET || detID==lcio::ILDDetID::FTD){ if(addSpacePointTrakerHit(hit,hitID)){ - GenfitMsg::get()<<MSG::DEBUG<<"add slicon space point"<<endmsg; + if(m_debug>=2)std::cout<<"add slicon space point"<<std::endl; hitID++; }else{ - GenfitMsg::get()<<MSG::ERROR<<"silicon addSpacePointTrakerHit" - <<detID<<" "<<hit.getCellID()<<" faieled"<<endmsg; + if(m_debug>=2)std::cout<<"silicon addSpacePointTrakerHit" + <<detID<<" "<<hit.getCellID()<<" faieled"<<std::endl; } }else if(detID==7){ //if(addSpacePointMeasurement(p,sigma,hit.getCellID(),hitID)){ - // GenfitMsg::get()<<MSG::DEBUG<<"add DC space point"<<endmsg; + // if(m_debug>=2)std::cout<<"add DC space point"<<std::endl; // hitID++; //}else{ - // GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurement" - // <<detID<<" faieled" <<endmsg; + // if(m_debug>=2)std::cout<<"addSpacePointMeasurement" + // <<detID<<" faieled" <<std::endl; //} float minTime=FLT_MAX; edm4hep::ConstSimTrackerHit minTimeSimHit; @@ -697,12 +696,12 @@ int GenfitTrack::addSimTrackerHits(const edm4hep::Track& track, //std::cout<<"minTimeSimHit "<<minTimeSimHit<<std::endl; sortedDCTrackHitCol.push_back(minTimeSimHit); }else{ - GenfitMsg::get()<<MSG::DEBUG<<" Skip add this hit!"<<endmsg; + if(m_debug>=2)std::cout<<" Skip add this hit!"<<std::endl; } }//end loop over hit on track - GenfitMsg::get()<<MSG::DEBUG<<" addSimTrakerHits trackerHits_size=" - <<track.trackerHits_size()<<endmsg; + if(m_debug>=2)std::cout<<" addSimTrakerHits trackerHits_size=" + <<track.trackerHits_size()<<std::endl; ///Add DC hits to track //Sort sim DC hits by time @@ -715,15 +714,15 @@ int GenfitTrack::addSimTrackerHits(const edm4hep::Track& track, p[2]=pos.z; unsigned long long detID = dCTrackerHit.getCellID(); if(addSpacePointMeasurement(p,sigma,detID,hitID,smear)){ - GenfitMsg::get()<<MSG::DEBUG<<"add DC space point"<<endmsg; + if(m_debug>=2)std::cout<<"add DC space point"<<std::endl; hitID++; }else{ - GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurement" - <<detID<<" faieled" <<endmsg; + if(m_debug>=2)std::cout<<"addSpacePointMeasurement" + <<detID<<" faieled" <<std::endl; } } - GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack nHitAdded "<<hitID<<endmsg; + if(m_debug>=2)std::cout<<"GenfitTrack nHitAdded "<<hitID<<std::endl; return hitID; } @@ -745,13 +744,13 @@ bool GenfitTrack::storeTrack(edm4hep::ReconstructedParticle& recParticle, TLorentzVector fittedPos; TVector3 fittedMom; int fittedState=getFittedState(fittedPos,fittedMom,fittedCov); - GenfitMsg::get()<<MSG::DEBUG<<"fit result: get status OK? pidType " + if(m_debug>=2)std::cout<<"fit result: get status OK? pidType " <<pidType<<" fittedState==0 " <<(0==fittedState)<<" isFitted "<<isFitted - <<" isConverged "<<isConverged<<" ndf "<<ndf<<endmsg; + <<" isConverged "<<isConverged<<" ndf "<<ndf<<std::endl; if((0!=fittedState) || (!isFitted) || !isConvergedFully || (ndf<ndfCut)){ - GenfitMsg::get()<<MSG::DEBUG<<" fitting failed"<<endmsg; + if(m_debug>=2)std::cout<<" fitting failed"<<std::endl; }else{ - GenfitMsg::get()<<MSG::DEBUG<<"fit result: Pos("<< + if(m_debug>=2)std::cout<<"fit result: Pos("<< fittedPos.X()<<" "<< fittedPos.Y()<<" "<< fittedPos.Z()<<") mom("<< @@ -762,7 +761,7 @@ bool GenfitTrack::storeTrack(edm4hep::ReconstructedParticle& recParticle, fittedMom.Perp()<< " chi2 "<<chi2<< " ndf "<<ndf - <<endmsg; + <<std::endl; } float pos[3]={float(fittedPos.X()),float(fittedPos.Y()),float(fittedPos.Z())}; float mom[3]={float(fittedMom.X()),float(fittedMom.Y()),float(fittedMom.Z())}; diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp deleted file mode 100644 index 18f19c8eac116e44c259725e39951ab97b481bd4..0000000000000000000000000000000000000000 --- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp +++ /dev/null @@ -1,404 +0,0 @@ -#include "RecGenfitAlgSDT.h" -#include "GenfitTrack.h" -#include "GenfitFitter.h" -#include "GenfitField.h" - -//genfit -#include "EventDisplay.h" - -//cepcsw -#include "DetInterface/IGeomSvc.h" -#include "DataHelper/HelixClass.h" -#include "DetSegmentation/GridDriftChamber.h" - -//externals -#include "edm4hep/EventHeaderCollection.h" -#include "edm4hep/MCParticle.h" -#include "edm4hep/MCParticleCollection.h" -#include "edm4hep/SimTrackerHitCollection.h" -#include "edm4hep/TrackerHitCollection.h" -#include "edm4hep/TrackCollection.h" -#include "edm4hep/MCRecoTrackerAssociationCollection.h" -#include "edm4hep/ReconstructedParticle.h" -#include "edm4hep/ReconstructedParticleCollection.h" -#include "edm4hep/Track.h" -#include "edm4hep/TrackCollection.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/Detector.h" -#include "UTIL/BitField64.h" -#include "DDSegmentation/Segmentation.h" -#include "TRandom.h" -#include "TLorentzVector.h" - -//stl -#include "time.h" - - -DECLARE_COMPONENT( RecGenfitAlgSDT ) - -///////////////////////////////////////////////////////////////////// -RecGenfitAlgSDT::RecGenfitAlgSDT(const std::string& name, - ISvcLocator* pSvcLocator): - GaudiAlgorithm(name, pSvcLocator),m_nPDG(5),m_dd4hep(nullptr), - m_gridDriftChamber(nullptr),m_decoder(nullptr) -{ - declareProperty("EventHeaderCollection", m_headerCol); - declareProperty("MCParticleCollection", m_mcParticleCol, - "Handle of the input MCParticle collection"); - declareProperty("DriftChamberHitCollection", m_DCSimHitCol, - "Handle of the input DC SimTrakerHit collection"); - declareProperty("DriftChamberDigiCollection", m_DCDigiCol, - "Handle of DC digi(TrakerHit) collection"); - declareProperty("DCTrackCollection", m_DCTrackCol, - "Handle of DC track collection"); - declareProperty("DCHitAssociationCollection", m_DCHitAssociationCol, - "Handle of simTrackerHit and TrackerHit association collection"); - declareProperty("SDTTrackCollection", m_SDTTrackCol, - "Handle of input silicon track collection"); - declareProperty("SDTRecParticleCollection", m_SDTRecParticleCol, - "Handle of silicon+drift chamber rec. particle collection"); -} - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * - -StatusCode RecGenfitAlgSDT::initialize() -{ - MsgStream log(msgSvc(), name()); - info()<<" RecGenfitAlgSDT initialize()"<<endmsg; - - ///Get GeomSvc - m_geomSvc=Gaudi::svcLocator()->service("GeomSvc"); - if (nullptr==m_geomSvc) { - std::cout<<"Failed to find GeomSvc"<<std::endl; - return StatusCode::FAILURE; - } - ///Get Detector - m_dd4hep = m_geomSvc->lcdd(); - ///Get Field - m_dd4hepField=m_geomSvc->lcdd()->field(); - - /// New a genfit fitter - m_genfitFitter=new GenfitFitter(m_fitterType.toString().c_str()); - m_genfitField=new GenfitField(m_dd4hepField); - m_genfitFitter->setField(m_genfitField); - m_genfitFitter->setGeoMaterial(m_geomSvc->lcdd()); - m_genfitFitter->setEnergyLossBrems(m_correctBremsstrahlung); - m_genfitFitter->setNoiseBrems(m_correctBremsstrahlung); - if(m_debug>10) m_genfitFitter->setDebug(m_debug-10); - if(m_noMaterialEffects) m_genfitFitter->setNoEffects(true); - if(-1==m_debugPid) m_genfitFitter->setNoEffects(true); - if(-1==m_debugPid) m_debugPid=0;//charged geantino with electron pid - if(m_fitterType=="DAF"||m_fitterType=="DafRef"){ - m_genfitFitter->setMaxIterationsBetas(m_bStart,m_bFinal,m_maxIteration); - } else { - m_genfitFitter->setMaxIterations(m_maxIteration); - } - //print genfit parameters - if(m_debug) m_genfitFitter->print(); - if(""!=m_genfitHistRootName) m_genfitFitter->initHist(m_genfitHistRootName); - - //initialize member vairables - for(int i=0;i<5;i++) m_fitSuccess[i]=0; - m_nRecTrack=0; - ///Get Readout - dd4hep::Readout readout=m_dd4hep->readout(m_readout_name); - ///Get Segmentation - m_gridDriftChamber=dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*> - (readout.segmentation().segmentation()); - if(nullptr==m_gridDriftChamber){ - error() << "Failed to get the GridDriftChamber" << endmsg; - return StatusCode::FAILURE; - } - ///Get Decoder - m_decoder = m_geomSvc->getDecoder(m_readout_name); - if (nullptr==m_decoder) { - error() << "Failed to get the decoder" << endmsg; - return StatusCode::FAILURE; - } - - - ///book tuple - NTuplePtr nt(ntupleSvc(), "RecGenfitAlgSDT/recGenfitAlgSDT"); - if(nt){ - m_tuple=nt; - }else{ - m_tuple=ntupleSvc()->book("RecGenfitAlgSDT/recGenfitAlgSDT", - CLID_ColumnWiseTuple,"RecGenfitAlgSDT"); - if(m_tuple){ - StatusCode sc; - sc=m_tuple->addItem("run",m_run); - sc=m_tuple->addItem("evt",m_evt); - sc=m_tuple->addItem("tkId",m_tkId); - sc=m_tuple->addItem("mcIndex",m_mcIndex,0,100);//max. 100 particles - sc=m_tuple->addItem("truthPocaMc",m_mcIndex,m_truthPocaMc,3); - sc=m_tuple->addItem("pocaPosMc",m_mcIndex,m_pocaPosMc,3); - sc=m_tuple->addItem("pocaMomMc",m_mcIndex,m_pocaMomMc,3); - sc=m_tuple->addItem("pocaMomMcP",m_mcIndex,m_pocaMomMcP); - sc=m_tuple->addItem("pocaMomMcPt",m_mcIndex,m_pocaMomMcPt); - sc=m_tuple->addItem("pocaPosMdc",3,m_pocaPosMdc); - sc=m_tuple->addItem("pocaMomMdc",3,m_pocaMomMdc); - sc=m_tuple->addItem("index",m_pidIndex, 0, 5); - sc=m_tuple->addItem("firstPosKalP",5,3,m_firstPosKal); - sc=m_tuple->addItem("firstMomKalP",5,m_firstMomKalP); - sc=m_tuple->addItem("firstMomKalPt",5,m_firstMomKalPt); - sc=m_tuple->addItem("pocaPosKal",5,3,m_pocaPosKal); - sc=m_tuple->addItem("pocaMomKal",5,3,m_pocaMomKal); - sc=m_tuple->addItem("pocaMomKalP",5,m_pocaMomKalP); - sc=m_tuple->addItem("pocaMomKalPt",5,m_pocaMomKalPt); - sc=m_tuple->addItem("chargeKal",5,m_chargeKal); - sc=m_tuple->addItem("nDofKal",5,m_nDofKal); - sc=m_tuple->addItem("chi2Kal",5,m_chi2Kal); - sc=m_tuple->addItem("isFitted",5,m_isFitted); - sc=m_tuple->addItem("isFitConverged",5,m_isFitConverged); - sc=m_tuple->addItem("isFitConvergedFully",5, - m_isFitConvergedFully); - sc=m_tuple->addItem("nHitFailedKal",5,m_nHitFailedKal); - sc=m_tuple->addItem("nHitFitted",5,m_nHitFitted); - sc=m_tuple->addItem("nDigi",m_nDigi); - sc=m_tuple->addItem("nHitMc",m_nHitMc); - sc=m_tuple->addItem("nHitKalInput",m_nHitKalInput); - sc=m_tuple->addItem("nHitWithFitInfo",5,m_nHitWithFitInfo); - sc=m_tuple->addItem("nSimDCHit",m_nSimDCHit,0,50000); - sc=m_tuple->addItem("mdcHitDriftT",m_nSimDCHit,m_mdcHitDriftT); - sc=m_tuple->addItem("mdcHitDriftDl",m_nSimDCHit,m_mdcHitDriftDl); - sc=m_tuple->addItem("mdcHitDriftDr",m_nSimDCHit,m_mdcHitDriftDr); - sc=m_tuple->addItem("mdcHitLr",m_nSimDCHit,m_mdcHitLr); - sc=m_tuple->addItem("mdcHitLayer",m_nSimDCHit,m_mdcHitLayer); - sc=m_tuple->addItem("mdcHitWire",m_nSimDCHit,m_mdcHitWire); - sc=m_tuple->addItem("mdcHitExpDoca",m_nSimDCHit,m_mdcHitExpDoca); - sc=m_tuple->addItem("mdcHitExpMcDoca",m_nSimDCHit,m_mdcHitExpMcDoca); - sc=m_tuple->addItem("mdcHitErr",m_nSimDCHit,m_mdcHitErr); - sc=m_tuple->addItem("time",5,m_time); - sc=m_tuple->addItem("mdcHitMcTkId",m_nSimDCHit,m_mdcHitMcTkId); - sc=m_tuple->addItem("mdcHitMcLr",m_nSimDCHit,m_mdcHitMcLr); - sc=m_tuple->addItem("mdcHitMcDrift",m_nSimDCHit,m_mdcHitMcDrift); - sc=m_tuple->addItem("mdcHitMcX",m_nSimDCHit,m_mdcHitMcX); - sc=m_tuple->addItem("mdcHitMcY",m_nSimDCHit,m_mdcHitMcY); - sc=m_tuple->addItem("mdcHitMcZ",m_nSimDCHit,m_mdcHitMcZ); - sc=m_tuple->addItem("mcPocaX",m_nSimDCHit,m_mdcHitExpMcPocaX); - sc=m_tuple->addItem("mcPocaY",m_nSimDCHit,m_mdcHitExpMcPocaY); - sc=m_tuple->addItem("mcPocaZ",m_nSimDCHit,m_mdcHitExpMcPocaZ); - sc=m_tuple->addItem("mcPocaWireX",m_nSimDCHit,m_mdcHitExpMcPocaWireX); - sc=m_tuple->addItem("mcPocaWireY",m_nSimDCHit,m_mdcHitExpMcPocaWireY); - sc=m_tuple->addItem("mcPocaWireZ",m_nSimDCHit,m_mdcHitExpMcPocaWireZ); - debug()<< "Book tuple RecGenfitAlgSDT/genfit" << endmsg; - }else{ - error()<< "Cannot book tuple RecGenfitAlgSDT/genfit" << endmsg; - } - }//end of book tuple - - //init genfit event display - //if(m_showDisplay) m_genfitDisplay = genfit::EventDisplay::getInstance(); - - return StatusCode::SUCCESS; -} - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * - -StatusCode RecGenfitAlgSDT::execute() -{ - StatusCode sc=StatusCode::SUCCESS; - m_timer=clock(); - info()<<"RecGenfitAlgSDT in execute()"<<endmsg; - - /////retrieve EventHeader - //auto header = _headerCol.get()->at(0); - //int evtNo = header.getEventNumber(); - //int runNo = header.getRunNumber(); - - //std::cout<<"run "<<header.getEventNumber() - // <<" "<<header.getRunNumber()<<std::endl; - - ///retrieve silicon Track and TrackHits - const edm4hep::TrackCollection* sdtTrackCol=nullptr; - sdtTrackCol=m_SDTTrackCol.get(); - if(nullptr==sdtTrackCol) { - debug()<<"TrackCollection not found"<<endmsg; - return StatusCode::SUCCESS; - } - //const edm4hep::TrackerHitCollection* dcDigiCol=nullptr; - //dcDigiCol=m_DCDigiCol.get(); - //if(nullptr==dcDigiCol) { - // debug()<<"DigiDCHitCollection not found"<<endmsg; - // return StatusCode::SUCCESS; - //} - ///retrieve DC Hit Association - auto dcHitAssociationCol=m_DCHitAssociationCol.get(); - - - edm4hep::ReconstructedParticleCollection* sdtRecParticleCol= - m_SDTRecParticleCol.createAndPut(); - - ///---------------------------------------------------- - ///Loop over Track and do fitting for each track - ///---------------------------------------------------- - debug()<<"SDTTrackCol size="<<sdtTrackCol->size()<<endmsg; - for(auto sdtTrack: *sdtTrackCol){ - ///Loop over 5 particle hypothesis(0-4): e,mu,pi,K,p - ///-1 for chargedgeantino - for(unsigned int pidType=0;pidType<m_nPDG;pidType++){ - if((m_debugPid>=0) && (m_debugPid!=pidType)) continue; - ///----------------------------------- - ///Create a GenFit track - ///----------------------------------- - GenfitTrack* genfitTrack=new GenfitTrack(m_genfitField, - m_gridDriftChamber); - genfitTrack->setDebug(m_debug); - double eventStartTime=0; - if(!genfitTrack->createGenfitTrackFromEDM4HepTrack(pidType,sdtTrack, - eventStartTime)){ - debug()<<"createGenfitTrack failed!"<<endmsg; - return StatusCode::SUCCESS; - } - if(0==genfitTrack->addSimTrackerHits(sdtTrack,dcHitAssociationCol, - m_sigmaHit.value())){ - debug()<<"addSimTrackerHits failed!"<<endmsg; - return StatusCode::SUCCESS; - } - if(m_debug) genfitTrack->printSeed(); - - ///----------------------------------- - ///call genfit fitting procedure - ///----------------------------------- - m_genfitFitter->processTrack(genfitTrack,m_resortHits); - m_genfitFitter->setDebug(m_debug); - - ///----------------------------------- - ///Store track - ///----------------------------------- - auto dcRecParticle=sdtRecParticleCol->create(); - genfitTrack->storeTrack(dcRecParticle,pidType,m_ndfCut, - m_chi2Cut); - if(m_debug) genfitTrack->printSeed(); - - if(m_tuple) debugTrack(pidType,genfitTrack); - if(m_showDisplay) { - //m_genfitDisplay->addEvent(genfitTrack->getTrack()); - //m_genfitDisplay->open(); - }else{ - delete genfitTrack; - } - ++m_fitSuccess[pidType]; - }//end loop over particle type - }//end loop over a track - m_nRecTrack++; - - if(m_tuple) debugEvent(); - - - //if(m_genfitDisplay) while(1){ - // std::cout<<"Press any key to finish..."<<std::endl; - // //system ("pause"); - //} - - if(m_tuple) sc=m_tuple->write(); - - return StatusCode::SUCCESS; -} - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * - -StatusCode RecGenfitAlgSDT::finalize() -{ - MsgStream log(msgSvc(), name()); - info()<< " RecGenfitAlgSDT in finalize()" << endmsg; - - m_genfitFitter->writeHist(); - delete m_genfitFitter; - info()<<"RecGenfitAlgSDT nRecTrack="<<m_nRecTrack<<" success e " - <<m_fitSuccess[0]<<" mu "<<m_fitSuccess[1]<<" pi "<<m_fitSuccess[2] - <<" K "<<m_fitSuccess[3]<<" p "<<m_fitSuccess[4]<<std::endl; - if(m_nRecTrack>0){ - std::cout<<"RecGenfitAlgSDT Success rate = "<<std::endl; - for (int i=0;i<5;i++){ - std::cout<<Form("%d %2.2f",i,((float) m_fitSuccess[i])/m_nRecTrack) - <<std::endl; - } - } - return StatusCode::SUCCESS; -} - -void RecGenfitAlgSDT::debugTrack(int pidType,const GenfitTrack* genfitTrack) -{ - /// Get fit status - const genfit::FitStatus* fitState = genfitTrack->getFitStatus(); - int charge= fitState->getCharge(); - - m_chargeKal[pidType]= charge; - m_nHitWithFitInfo[pidType]=genfitTrack->getNumPointsWithFittedInfo(); - m_chi2Kal[pidType]=fitState->getChi2(); - m_nDofKal[pidType]=fitState->getNdf(); - m_isFitted[pidType]=(int)fitState->isFitted(); - m_isFitConverged[pidType]=(int) fitState->isFitConverged(); - m_isFitConvergedFully[pidType]=(int) fitState->isFitConvergedFully(); - - ///get fitted state of track - TMatrixDSym fittedCov; - TLorentzVector fittedPos; - TVector3 fittedMom; - int fittedState=genfitTrack->getFittedState(fittedPos,fittedMom,fittedCov); - HelixClass helix;//mm and GeV - float pos[3]={float(fittedPos.X()/dd4hep::mm),float(fittedPos.Y()/dd4hep::mm), - float(fittedPos.Z()/dd4hep::mm)}; - float mom[3]={float(fittedMom.X()),float(fittedMom.Y()),float(fittedMom.Z())}; - helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect())); - m_pocaMomKalP[pidType]=fittedMom.Mag(); - - if(m_debug>0){ - /// Get fit status - debug()<<"evt "<<m_evt<<" fit result: get status OK? pidType " - <<pidType<<" fittedState "<<fittedState<<" isFitted " - <<m_isFitted[pidType]<<" isConverged "<<m_isFitConverged[pidType] - <<" isFitConvergedFully "<<m_isFitConvergedFully[pidType] - <<" ndf "<<m_nDofKal[pidType] - <<" chi2 "<<m_chi2Kal[pidType]<<endmsg; - if((0!=fittedState)||(!m_isFitted[pidType])||(m_nDofKal[pidType]<m_ndfCut)){ - debug()<<"fitting failed"<<endmsg; - }else{ - debug()<<"evt "<<m_evt<<" fit result: Pos("<< - fittedPos.X()<<" "<< - fittedPos.Y()<<" "<< - fittedPos.Z()<<") mom("<< - fittedMom.X()<<" "<< - fittedMom.Y()<<" "<< - fittedMom.Z()<<") p_tot "<< - fittedMom.Mag()<<" pt "<< - fittedMom.Perp()<<endmsg; - } - } -} - -void RecGenfitAlgSDT::debugEvent() -{ - const edm4hep::MCParticleCollection* mcParticleCol = nullptr; - const edm4hep::SimTrackerHitCollection* simDCHitCol=nullptr; - - m_pidIndex=5; - - mcParticleCol=m_mcParticleCol.get(); - simDCHitCol=m_simDCHitCol.get(); - m_nSimDCHit=simDCHitCol->size(); - int iMcParticle=0; - int iHit=0; - for(auto mcParticle : *mcParticleCol){ - for(auto simDCHit: *simDCHitCol){ - edm4hep::Vector3d pos=simDCHit.position(); - TVectorD p(3); - p[0]=pos.x;//no unit conversion here - p[1]=pos.y; - p[2]=pos.z; - m_mdcHitMcX[iHit]=pos.x; - m_mdcHitMcY[iHit]=pos.y; - m_mdcHitMcZ[iHit]=pos.z; - iHit++; - } - edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV - float px=mcPocaMom.x; - float py=mcPocaMom.y; - float pz=mcPocaMom.z; - debug()<<" "<<px<<" "<<py<<" "<<pz<<endmsg; - m_pocaMomMcP[iMcParticle]=sqrt(px*px+py*py+pz*pz); - iMcParticle++; - } - m_mcIndex=iHit; -} diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h deleted file mode 100644 index 1a1031e1b3db3b17d5c509ebbed76304f09a1600..0000000000000000000000000000000000000000 --- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h +++ /dev/null @@ -1,192 +0,0 @@ -////////////////////////////////////////////////////////////////////// -/// -/// This is an algorithm for track fitting for CEPC track with genfit. -/// -/// In this file, including: -/// An algorithm for combined silicon and drift chamber track fitting -/// with genfit for 5 particle hypothesis -/// -/// Units are following DD4hepUnits -/// -/// Authors: -/// Yao ZHANG(zhangyao@ihep.ac.cn) -/// -///////////////////////////////////////////////////////////////////// - -#ifndef RECGENFITALG_RECGENFITALGSDT_H -#define RECGENFITALG_RECGENFITALGSDT_H - -#include "GaudiAlg/GaudiAlgorithm.h" -#include "GaudiKernel/NTuple.h" -#include "k4FWCore/DataHandle.h" -#include "DD4hep/Fields.h" -#include <string> - -class GenfitFitter; -class GenfitField; -class GenfitTrack; -class IGeomSvc; -class time; -namespace genfit{ - class EventDisplay; -} -namespace dd4hep { - class Detector; - //class rec::CellIDPositionConverter; - namespace DDSegmentation{ - class GridDriftChamber; - class BitFieldCoder; - } -} -namespace edm4hep{ - class EventHeaderCollection; - class MCParticleCollection; - class SimTrackerHitCollection; - class TrackCollection; - class TrackerHitCollection; - class MCRecoTrackerAssociationCollection; - class ReconstructedParticle; - class ReconstructedParticleCollection; -} - -///////////////////////////////////////////////////////////////////////////// - -class RecGenfitAlgSDT:public GaudiAlgorithm { - public: - RecGenfitAlgSDT (const std::string& name, ISvcLocator* pSvcLocator); - StatusCode initialize() override; StatusCode execute() override; StatusCode finalize() override; - - private: - GenfitFitter* m_genfitFitter;//The pointer to a GenfitFitter - const GenfitField* m_genfitField;//The pointer to a GenfitField - - void debugTrack(int pidType,const GenfitTrack* genfitTrack); - void debugEvent(); - - DataHandle<edm4hep::EventHeaderCollection> m_headerCol{ - "EventHeaderCol", Gaudi::DataHandle::Reader, this}; - //Drift chamber rec hit and trac - DataHandle<edm4hep::TrackerHitCollection> m_DCDigiCol{ - "DigiDCHitCollection", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::SimTrackerHitCollection> m_DCSimHitCol{ - "DCSimHitsCollection", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::TrackCollection> m_DCTrackCol{ - "DCTrackCollection", Gaudi::DataHandle::Reader, this}; - //Mc truth - DataHandle<edm4hep::MCParticleCollection> m_mcParticleCol{ - "MCParticle", Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::SimTrackerHitCollection> m_simDCHitCol{ - "DriftChamberHitsCollection" , Gaudi::DataHandle::Reader, this}; - DataHandle<edm4hep::MCRecoTrackerAssociationCollection> - m_DCHitAssociationCol{"DCHitAssociationCollection", - Gaudi::DataHandle::Reader, this}; - - //Track from silicon detectors - DataHandle<edm4hep::TrackCollection> m_SDTTrackCol{"SDTTrackCollection", - Gaudi::DataHandle::Reader, this}; - - //Output hits and particles - DataHandle<edm4hep::ReconstructedParticleCollection> m_SDTRecParticleCol{ - "SDTRecParticleCollection", Gaudi::DataHandle::Writer, this}; - - const unsigned int m_nPDG;//5:e,mu,pi,K,proton - SmartIF<IGeomSvc> m_geomSvc; - dd4hep::OverlayedField m_dd4hepField; - dd4hep::Detector* m_dd4hep; - dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber; - dd4hep::DDSegmentation::BitFieldCoder* m_decoder; - Gaudi::Property<std::string> m_readout_name{this, - "readout", "DriftChamberHitsCollection"}; - Gaudi::Property<int> m_debug{this,"debug",0}; - Gaudi::Property<float> m_sigmaHit{this,"sigmaHit",0.11};//mm - Gaudi::Property<float> m_nSigmaHit{this,"nSigmaHit",5}; - Gaudi::Property<double> m_initCovResPos{this,"initCovResPos",1}; - Gaudi::Property<double> m_initCovResMom{this,"initCovResMom",0.1}; - //Fitter type default is DAFRef. - //Candidates are DAF,DAFRef,KalmanFitter and KalmanFitterRefTrack. - Gaudi::Property<std::string> m_fitterType{this,"fitterTyep","DAFRef"}; - Gaudi::Property<bool> m_correctBremsstrahlung{this, - "correctBremsstrahlung",false}; - Gaudi::Property<bool> m_noMaterialEffects{this, - "noMaterialEffects",false}; - Gaudi::Property<int> m_maxIteration{this,"maxIteration",20}; - Gaudi::Property<int> m_resortHits{this,"resortHits",true}; - Gaudi::Property<double> m_bStart{this,"bStart",100}; - Gaudi::Property<double> m_bFinal{this,"bFinal",0.01}; - Gaudi::Property<double> m_DCCornerCuts{this,"dcCornerCuts",-999}; - Gaudi::Property<double> m_ndfCut{this,"ndfCut",1}; - Gaudi::Property<double> m_chi2Cut{this,"chi2Cut",1000}; - //-1,chargedGeantino;0,1,2,3,4:e,mu,pi,K,proton - Gaudi::Property<int> m_debugPid{this,"debugPid",-99}; - Gaudi::Property<bool> m_useTruthTrack{this,"useTruthTrack",true}; - Gaudi::Property<bool> m_useTruthHit{this,"useTruthHit",true}; - Gaudi::Property<std::string> m_genfitHistRootName{this, - "genfitHistRootName",""}; - Gaudi::Property<bool> m_showDisplay{this,"showDisplay",false}; - int m_fitSuccess[5]; - int m_nRecTrack; - //bool m_useRecLRAmbig; - - genfit::EventDisplay* m_genfitDisplay; - clock_t m_timer; - - /// tuples - NTuple::Tuple* m_tuple; - NTuple::Item<int> m_run; - NTuple::Item<int> m_evt; - NTuple::Item<int> m_tkId; - NTuple::Item<int> m_mcIndex;//number of navigated mcParicle - NTuple::Matrix<double> m_truthPocaMc;//2 dim matched particle and 3 pos. - NTuple::Matrix<double> m_pocaPosMc;//2 dim matched particle and 3 pos. - NTuple::Matrix<double> m_pocaMomMc;//2 dim matched particle and 3 mom. - NTuple::Array<double> m_pocaMomMcP;//2 dim matched particle and p - NTuple::Array<double> m_pocaMomMcPt;//2 dim matched particle and pt - NTuple::Array<double> m_pocaPosMdc;//pos 0:x,1:y,2:z - NTuple::Array<double> m_pocaMomMdc;//mom. 0:px,1:py,2:pz - NTuple::Item<int> m_pidIndex; - NTuple::Matrix<double> m_firstPosKal;//5 hyposis and pos. at first - NTuple::Array<double> m_firstMomKalP;//5 hyposis and mom. at first - NTuple::Array<double> m_firstMomKalPt;//5 hyposis and mom. at first - NTuple::Matrix<double> m_pocaPosKal;//5 hyposis and 3 mom. - NTuple::Matrix<double> m_pocaMomKal;//5 hyposis and 3 mom. - NTuple::Array<double> m_pocaMomKalP;//5 hyposis and p - NTuple::Array<double> m_pocaMomKalPt;//5 hyposis and pt - NTuple::Array<int> m_chargeKal; - NTuple::Array<double> m_chi2Kal; - NTuple::Array<double> m_nDofKal; - NTuple::Array<int> m_isFitConverged; - NTuple::Array<int> m_isFitConvergedFully; - NTuple::Array<int> m_isFitted; - NTuple::Item<int> m_nDigi; - NTuple::Item<int> m_nHitMc; - NTuple::Item<int> m_nSimDCHit; - NTuple::Array<int> m_nHitWithFitInfo; - NTuple::Item<int> m_nHitKalInput; - NTuple::Array<double> m_mdcHitDriftT; - NTuple::Array<double> m_mdcHitDriftDl; - NTuple::Array<double> m_mdcHitDriftDr; - NTuple::Array<int> m_mdcHitLr; - NTuple::Array<int> m_mdcHitLayer; - NTuple::Array<int> m_mdcHitWire; - NTuple::Array<double> m_mdcHitExpDoca; - NTuple::Array<double> m_mdcHitExpMcDoca; - NTuple::Array<double> m_mdcHitErr; - NTuple::Array<int> m_nHitFailedKal; - NTuple::Array<int> m_nHitFitted; - NTuple::Array<double> m_time; - //truth - NTuple::Array<int> m_mdcHitMcLr; - NTuple::Array<int> m_mdcHitMcTkId; - NTuple::Array<double> m_mdcHitMcDrift; - NTuple::Array<double> m_mdcHitMcX; - NTuple::Array<double> m_mdcHitMcY; - NTuple::Array<double> m_mdcHitMcZ; - NTuple::Array<double> m_mdcHitExpMcPocaX; - NTuple::Array<double> m_mdcHitExpMcPocaY; - NTuple::Array<double> m_mdcHitExpMcPocaZ; - NTuple::Array<double> m_mdcHitExpMcPocaWireX; - NTuple::Array<double> m_mdcHitExpMcPocaWireY; - NTuple::Array<double> m_mdcHitExpMcPocaWireZ; - -}; -#endif diff --git a/Utilities/KalDet/CMakeLists.txt b/Utilities/KalDet/CMakeLists.txt index faacdc5b6132c6762b8098a79a6b90184f59526c..12fbb25371c389078d6787aff43cbda32f07b8cd 100644 --- a/Utilities/KalDet/CMakeLists.txt +++ b/Utilities/KalDet/CMakeLists.txt @@ -12,7 +12,7 @@ endif() set( DICT_CINT_DEFINITIONS "HANDLE_DICT_EXCEPTIONS=IGNORED_FOR_CINT" ) set( DICT_INPUT_DIRS gen kern lctpc/gearTPC ) -set( DICT_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}/rootdict ) +set( DICT_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}) foreach( DICT_DIR ${DICT_INPUT_DIRS} ) list( APPEND DICT_INCLUDE_DIRS ${CMAKE_CURRENT_SOURCE_DIR}/src/${DICT_DIR} ) diff --git a/Utilities/KalTest/CMakeLists.txt b/Utilities/KalTest/CMakeLists.txt index d1a29b90d61522b3804673af3636b3ba62641694..be765815acad2802a421885f49ba1d5fa339f4c7 100644 --- a/Utilities/KalTest/CMakeLists.txt +++ b/Utilities/KalTest/CMakeLists.txt @@ -5,7 +5,7 @@ set( DICT_CINT_DEFINITIONS "HANDLE_DICT_EXCEPTIONS=IGNORED_FOR_CINT" ) set( DICT_INPUT_DIRS geomlib kallib kaltracklib utils ) -set( DICT_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}/rootdict ) +set( DICT_OUTPUT_DIR ${CMAKE_CURRENT_BINARY_DIR}) foreach( DICT ${DICT_INPUT_DIRS} ) list( APPEND DICT_INCLUDE_DIRS ${CMAKE_CURRENT_SOURCE_DIR}/src/${DICT} )