From 487307ac3de7642e1d8dcc742aba90e60c1475e1 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Fri, 27 Aug 2021 16:40:44 +0200 Subject: [PATCH 01/22] Fix phi flip and abs() bugs --- TimeOfFlight/src/TOFEstimators.cc | 13 ++-- TimeOfFlight/src/TOFUtils.cc | 103 ++++++++++++++++-------------- 2 files changed, 63 insertions(+), 53 deletions(-) diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index cb66fbd1..03c84e8d 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -97,8 +97,8 @@ void TOFEstimators::init() { const auto& detector = dd4hep::Detector::getInstance(); detector.field().magneticField({0., 0., 0.}, _bField); } - else { - throw EVENT::Exception(std::string("Invalid ProcessorVersion parameter passed: \'") + _procVersion + std::string("\'\n Viable options are idr (default), dev")); + else { + throw EVENT::Exception(std::string("Invalid ProcessorVersion parameter passed: \'") + _procVersion + std::string("\'\n Viable options are idr (default), dev")); } } @@ -497,13 +497,18 @@ dd4hep::rec::Vector3D TOFEstimators::getMomAtCalo(const Track* track){ double TOFEstimators::getFlightLength(const Track* track){ const TrackState* ts = track->getTrackState(TrackState::AtIP); - double phiIP = ts->getPhi(); + double phiIp = ts->getPhi(); ts = track->getTrackState(TrackState::AtCalorimeter); double phiCalo = ts->getPhi(); double omegaCalo = ts->getOmega(); double tanLCalo = ts->getTanLambda(); - return abs((phiIP - phiCalo)/omegaCalo)*sqrt(1. + tanLCalo*tanLCalo); + double dPhi = std::abs(phiIp - phiCalo); + + // if dPhi > PI assume phi is flipped + if (dPhi > M_PI) dPhi = 2*M_PI - dPhi; + + return dPhi/std::abs(omegaCalo)*std::sqrt(1. + tanLCalo*tanLCalo); } diff --git a/TimeOfFlight/src/TOFUtils.cc b/TimeOfFlight/src/TOFUtils.cc index 95168ed9..0ef9d235 100644 --- a/TimeOfFlight/src/TOFUtils.cc +++ b/TimeOfFlight/src/TOFUtils.cc @@ -15,50 +15,55 @@ namespace TOFUtils{ - using namespace lcio ; - - - float computeFlightLength( EVENT::Track* trk){ - - const TrackState* tsIP = trk->getTrackState( TrackState::AtIP ) ; - const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ; - - float Phicalo = tscalo->getPhi() ; - float PhiIP = tsIP->getPhi() ; - - float Omega = tsIP->getOmega() ; - float tanL = tsIP->getTanLambda() ; - - float length = (PhiIP-Phicalo)*(1/Omega) * sqrt( 1 + tanL * tanL ) ; - - return length ; - } +using namespace lcio ; + + +float computeFlightLength( EVENT::Track* trk){ + + const TrackState* tsIp = trk->getTrackState( TrackState::AtIP ) ; + const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ; + + float phiIp = tsIp->getPhi() ; + float phiCalo = tscalo->getPhi() ; + + float omega = tsIp->getOmega() ; + float tanL = tsIp->getTanLambda() ; + + float dPhi = std::abs(phiIp - phiCalo); - float computeFlightLength(const TrackState* ts0, const TrackState* ts1 ){ - - float Phicalo = ts1->getPhi() ; - float PhiIP = ts0->getPhi() ; - - float Omega = ts0->getOmega() ; + // if dPhi > PI assume phi is flipped + if (dPhi > M_PI) dPhi = 2*M_PI - dPhi; + + return dPhi/std::abs(omega) * std::sqrt( 1. + tanL * tanL ) ; +} + +float computeFlightLength(const TrackState* ts0, const TrackState* ts1 ){ + float phi1 = ts1->getPhi() ; + float phi0 = ts0->getPhi() ; + + float omega = ts0->getOmega() ; float tanL = ts0->getTanLambda() ; - - float length = (PhiIP-Phicalo)*(1/Omega) * sqrt( 1 + tanL * tanL ) ; - return length ; - } + float dPhi = std::abs(phi0 - phi1); + + // if dPhi > PI assume phi is flipped + if (dPhi > M_PI) dPhi = 2*M_PI - dPhi; + + return dPhi/std::abs(omega) * std::sqrt( 1. + tanL * tanL ) ; +} int layer( EVENT::CalorimeterHit* h ) { - + return CHT( h->getType() ).layer() ; } - - - /// helper function to check if this is an Ecal hit + + + /// helper function to check if this is an Ecal hit bool isEcal( EVENT::CalorimeterHit* h ) { - - return CHT( h->getType() ).caloID() == CHT::ecal ; - } + + return CHT( h->getType() ).caloID() == CHT::ecal ; + } std::string caloTypeStr( EVENT::CalorimeterHit* h ) { @@ -72,7 +77,7 @@ namespace TOFUtils{ std::string CaloHitData::toString(){ std::stringstream s ; - s << " l= " << layer ; + s << " l= " << layer ; s << ", st= " << smearedTime ; s << ", tr= " << timeResolution ; s << ", dIP=" << distanceFromIP ; @@ -87,30 +92,30 @@ namespace TOFUtils{ float computeDistanceFromLine( EVENT::CalorimeterHit* h, const dd4hep::rec::Vector3D& point, const dd4hep::rec::Vector3D& unitDir) { - dd4hep::rec::Vector3D pos( h->getPosition()[0], - h->getPosition()[1], + dd4hep::rec::Vector3D pos( h->getPosition()[0], + h->getPosition()[1], h->getPosition()[2] ) ; - + dd4hep::rec::Vector3D diff = pos - point ; - + return diff.cross( unitDir ).r() ; - + } - + CaloHitDataVec findHitsClosestToLine( const CaloHitLayerMap& layerMap ){ - + CaloHitDataVec hitVec ; - + for( auto m : layerMap ){ - + const CaloHitDataVec& chv = m.second ; - + CaloHitData* closestHit = *std::min_element( chv.begin() , chv.end () , [](CaloHitData* c0, CaloHitData* c1 ){ return c0->distancefromStraightline < c1->distancefromStraightline ; } - ) ; + ) ; hitVec.push_back( closestHit ) ; } @@ -128,14 +133,14 @@ namespace TOFUtils{ float meansq = 0. ; for( auto ch : chv ){ - float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; + float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; mean += t ; ++nHit ; } mean /= nHit ; for( auto ch : chv ){ - float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; + float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; meansq += ( t - mean ) * ( t - mean ) ; } From 492ddaca8e1cb7d4447c31b648283318258cacf3 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Fri, 17 Sep 2021 00:54:15 +0200 Subject: [PATCH 02/22] Major update. First draft compiling commit --- TimeOfFlight/README.md | 22 +- TimeOfFlight/example/tofplots_steer.xml | 124 --- TimeOfFlight/include/TOFEstimators.h | 159 +--- TimeOfFlight/include/TOFPlots.h | 74 -- TimeOfFlight/include/TOFUtils.h | 96 +-- TimeOfFlight/scripts/draw_cluster_time.C | 99 --- TimeOfFlight/scripts/draw_tofestimators.C | 70 -- TimeOfFlight/scripts/lcio_particle_gun.py | 192 ----- TimeOfFlight/scripts/make_root_plots.sh | 26 - .../run_all_and_make_root_plots_random.sh | 21 - TimeOfFlight/scripts/run_all_marlin.sh | 11 - TimeOfFlight/scripts/run_all_marlin_10GeV.sh | 22 - TimeOfFlight/scripts/run_all_marlin_2GeV.sh | 22 - TimeOfFlight/scripts/run_all_marlin_5GeV.sh | 22 - TimeOfFlight/scripts/run_all_marlin_v1.sh | 45 -- TimeOfFlight/scripts/run_all_marlin_v2.sh | 12 - TimeOfFlight/scripts/tofestimators.xml | 129 --- TimeOfFlight/src/TOFEstimators.cc | 751 ++++-------------- TimeOfFlight/src/TOFUtils.cc | 365 ++++++--- TimeOfFlight/xml/steer.xml | 43 + 20 files changed, 532 insertions(+), 1773 deletions(-) delete mode 100644 TimeOfFlight/example/tofplots_steer.xml delete mode 100644 TimeOfFlight/include/TOFPlots.h delete mode 100644 TimeOfFlight/scripts/draw_cluster_time.C delete mode 100644 TimeOfFlight/scripts/draw_tofestimators.C delete mode 100644 TimeOfFlight/scripts/lcio_particle_gun.py delete mode 100644 TimeOfFlight/scripts/make_root_plots.sh delete mode 100644 TimeOfFlight/scripts/run_all_and_make_root_plots_random.sh delete mode 100644 TimeOfFlight/scripts/run_all_marlin.sh delete mode 100644 TimeOfFlight/scripts/run_all_marlin_10GeV.sh delete mode 100644 TimeOfFlight/scripts/run_all_marlin_2GeV.sh delete mode 100644 TimeOfFlight/scripts/run_all_marlin_5GeV.sh delete mode 100644 TimeOfFlight/scripts/run_all_marlin_v1.sh delete mode 100644 TimeOfFlight/scripts/run_all_marlin_v2.sh delete mode 100644 TimeOfFlight/scripts/tofestimators.xml create mode 100644 TimeOfFlight/xml/steer.xml diff --git a/TimeOfFlight/README.md b/TimeOfFlight/README.md index 842b2f13..2b113a0c 100644 --- a/TimeOfFlight/README.md +++ b/TimeOfFlight/README.md @@ -9,25 +9,13 @@ based on the timing information stored in CalorimeterHits. ### TOFEstimators +author: N.Weinhold (DESY, internship 2017) author: F.Gaede, DESY, 2018 - -Compute various estimators for time of flight and add these to the +author: B.Dudar, DESY, 2021 +Compute various estimators for time of flight and add these to the ReconstrucedParticles as PIDParameters. -### TOFPlots - -author: N.Weinhold (DESY, internship 2017) - -Computes various estimators for the time of flight from CalorimeterHits. -Creates ROOT histograms for these parameters. - - -## Examples - -See [./scripts/tofestimators.xml](./scripts/tofestimators.xml) for and example -steering file that computes TOF estimators for 0,10 and 50 ps single hit -time resolution. -Use [./scripts/draw_tofestimators.C](./scripts/draw_tofestimators.C) to create -plots of beta vs. momentum for the diffefrent estimators. +## Example steering file +[./xml/steer.xml](./xml/steer.xml) is an example of a steering file that runs TOFEstimators processor. It computes momentum harmonic mean, time-of-flight and track length of the PandoraPFOs and stores this information inside PIDHandler. In the example time of flight calculated using closest Ecal hit assuming perfect time resolution. Then new slcio file "test.slcio" produced with LCIOOutputProcessor processor. diff --git a/TimeOfFlight/example/tofplots_steer.xml b/TimeOfFlight/example/tofplots_steer.xml deleted file mode 100644 index eeced279..00000000 --- a/TimeOfFlight/example/tofplots_steer.xml +++ /dev/null @@ -1,124 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - simjob.slcio - - - - - - - MESSAGE - - - - - - - - - 1 - - aida_file - - root - - - - - - - - MCParticle - - DEBUG - - - - - - - - 1 0 - 2 0 - 3 0 - - - - - - - - - - - - - - - - - outputfile.slcio - - WRITE_NEW - - - - - - - - - - - - 1 - - - - - - - 10000 - - - - - - - - input.stdhep - - - - - - - - - - - - - diff --git a/TimeOfFlight/include/TOFEstimators.h b/TimeOfFlight/include/TOFEstimators.h index a71d87b1..cf669a49 100644 --- a/TimeOfFlight/include/TOFEstimators.h +++ b/TimeOfFlight/include/TOFEstimators.h @@ -1,29 +1,10 @@ #ifndef TOFEstimators_h #define TOFEstimators_h 1 -#include "marlin/Processor.h" -#include "lcio.h" - -//std stuff -#include #include #include -#include -#include -#include -#include - -#include - -#include -#include - -#include "DDRec/Vector3D.h" - -using namespace lcio ; -using namespace marlin ; - -#include "TH2F.h" +#include "marlin/Processor.h" +#include "MarlinTrk/IMarlinTrkSystem.h" /** * @section DESCRIPTION @@ -101,118 +82,36 @@ using namespace marlin ; * * @file * @author F.Gaede, DESY, April 2018 - * @author B.Dudar, DESY, February 2021 + * @author B.Dudar, DESY, September 2021 */ -class TOFEstimators : public Processor { +class TOFEstimators : public marlin::Processor { public: - - virtual Processor* newProcessor() { return new TOFEstimators ; } - - TOFEstimators(const TOFEstimators&) = delete; - TOFEstimators& operator=(const TOFEstimators&) = delete; - - TOFEstimators() ; - - /** Called at the begin of the job before anything is read. - * Use to initialize the processor, e.g. book histograms. - */ - virtual void init() ; - - /** Called for every event - the working horse. - */ - virtual void processEvent(LCEvent* evt) ; - - /** Called after data processing to plot some histograms to debug. - */ - virtual void check(LCEvent* evt) ; - - /** Called after data processing for clean up. - */ - virtual void end() ; - - /** Called in dev version to write MomAtCalo parameter. - Can be used for more precise mass calculation - */ - dd4hep::rec::Vector3D getMomAtCalo(const Track*); - - /** Track length of the helix between IP and ECAL entrance - * assuming IP at (0, 0, 0) and ECAL entry point is - * obtained from the Kalman filter by extrapolation. - * Omega/lambda track parameters are taken from the - * closest track hit to the ECAL and assumed to be - * constant along the track. - * Length is calculated by the formula that is reasonable only for delta phi < 2pi - * meaning only tracks with less than one curl. - */ - double getFlightLength(const Track*); - - /** - * TOFClosest - ToF of the CalorimeterHit closest to track entry - * point to the ECAL corrected with the distance - * between cell center and entry point assuming - * speed of light propagation - */ - double getTOFClosest(const Track*, const Cluster*); - - /** - * TOFFastest - ToF of the CalorimeterHit arrived fastest corrected - * with the distance between cell center and entry - * point assuming speed of light propagation - */ - double getTOFFastest(const Track*, const Cluster*); - - /** - * TOFCylFit - ToF estimated from extrapolation from a fit of the - * hit times vs distance to to the track entry point - * to the ECAL at distance = 0 - */ - double getTOFCylFit(const Track*, const Cluster*); - - /** - * TOFClosestFit - ToF estimated with a fit as TOFCylFit but from hits - * that are closest to the linear prolongation of the - * trajectory line of the track in MaxLayerNumber first layers - */ - double getTOFClosestFit(const Track*, const Cluster*); - - protected: - - /** Input collection name with ReconstructedParticles (PFOs). - used for ReconstructedParticleCollection steering parameter. usually PandoraPFOs - */ - std::string _colNamePFO{}; - - /** Input version idr (default) or dev (current development). - */ - std::string _procVersion{}; - - /** Select ECAL hits only from _maxLayerNum first layers of ECAL for ToF calculations - */ - int _maxLayerNum{}; - - - /** Select hits for tofCylFit method from the cylinder limited by this radius - */ - double _cylRadiusCut{}; - - /** Single ECAL hit time resolution in ps - */ - float _resolution{}; - - /** vector of names of output PID object parameters which are written as output - */ - std::vector _TOFNames{}; - - gsl_rng* _rng = nullptr; - std::vector _h{}; - - std::mt19937 _generator{}; - std::normal_distribution _smearing{}; - - /** B field from the DD4HEP geometry processor for MomAtCalo calculations - */ - double _bField[3]; + TOFEstimators(const TOFEstimators&) = delete; + TOFEstimators& operator=(const TOFEstimators&) = delete; + + marlin::Processor* newProcessor() { return new TOFEstimators; } + TOFEstimators(); + void init(); + void processEvent(EVENT::LCEvent* evt); + private: + std::string _pfoCollectionName{}; + bool _extrapolateToEcal{}; + std::string _tofMethod{}; + double _timeResolution{}; + int _maxEcalLayer{}; + + + int _nEvent{}; + /** vector of names of output PID object parameters which are written as output + */ + std::vector _outputParNames{}; + + // MarlinTrk v02-00 release notes: + // USERS SHOULD NO LONGER DELETE THE IMarlinTrkSystem POINTER IN THEIR CODE + MarlinTrk::IMarlinTrkSystem* _trkSystem = nullptr; + double _bField{}; + double _tpcOuterR{}; }; #endif diff --git a/TimeOfFlight/include/TOFPlots.h b/TimeOfFlight/include/TOFPlots.h deleted file mode 100644 index 7b674380..00000000 --- a/TimeOfFlight/include/TOFPlots.h +++ /dev/null @@ -1,74 +0,0 @@ -#ifndef MyProcessor_h -#define MyProcessor_h 1 - -#include "marlin/Processor.h" -#include "lcio.h" -#include -#include -#include - -using namespace lcio ; -using namespace marlin ; - - -class TH1 ; - -/** Computes various estimators for the time of flight from CalorimeterHits. - * Creates ROOT histograms for these parameters. - * - * @author N. Weinhold, DESY internship 2017 - */ - -class TOFPlots : public Processor { - - public: - - virtual Processor* newProcessor() { return new TOFPlots ; } - - TOFPlots(const TOFPlots&) = delete; - TOFPlots& operator=(const TOFPlots&) = delete; - - TOFPlots() ; - - /** Called at the begin of the job before anything is read. - * Use to initialize the processor, e.g. book histograms. - */ - virtual void init() ; - - /** Called for every run. - */ - virtual void processRunHeader( LCRunHeader* run ) ; - - /** Called for every event - the working horse. - */ - virtual void processEvent( LCEvent * evt ) ; - - - virtual void check( LCEvent * evt ) ; - - - /** Called after data processing for clean up. - */ - virtual void end() ; - - - protected: - - /** Input collection name. - */ - std::string _colNameMCP{}; - std::string _colNamePFO{}; - - int _nRun{}; - int _nEvt{}; - - gsl_rng* _rng = nullptr ; - - std::vector _h{}; - -} ; - -#endif - - - diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index 4647ac4d..08f5335e 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -2,90 +2,68 @@ #define TOFUtils_h 1 /****************************************************** - * Utilities for cumputing time of flight estimators. - * + * Utility functions that are used by the TOFEstimators processor + * * @author F. gaede, DESY, 2018 - * + * @author B. Dudar, DESY, 2021 ****************************************************** */ #include -#include -#include - -#include "EVENT/CalorimeterHit.h" -#include "EVENT/Track.h" +#include "EVENT/ReconstructedParticle.h" +#include "IMPL/TrackStateImpl.h" +#include "MarlinTrk/IMarlinTrkSystem.h" +#include "MarlinTrk/IMarlinTrack.h" #include "DDRec/Vector3D.h" -namespace EVENT{ - class Track ; -} - +namespace TOFUtils{ + bool sortByRho(EVENT::TrackerHit* a, EVENT::TrackerHit* b); -namespace TOFUtils{ + IMPL::TrackStateImpl getTrackStateAtHit(MarlinTrk::IMarlinTrack* marlinTrk, EVENT::TrackerHit* hit); - /// handle for calorimeter hit meta data used for TOF estimators - struct CaloHitData{ - CaloHitData() = delete ; - CaloHitData(const CaloHitData&) = default ; - CaloHitData& operator=(const CaloHitData&) = default ; - CaloHitData(EVENT::CalorimeterHit* h) : lcioHit(h) {} ; - EVENT::CalorimeterHit* lcioHit = nullptr ; - int layer = 0 ; - float smearedTime = 0. ; - float timeResolution = 0. ; - float distanceFromIP = 0. ; - float distanceFromReferencePoint = 0. ; - float distancefromStraightline = 0. ; - std::string toString() ; - } ; - - /// define a vector of unique pointers to extended calo hits - typedef std::vector< std::unique_ptr > CaloHitUPtrVec ; + // gets momentum of the helix from parameters of the trackState + dd4hep::rec::Vector3D getHelixMomAtTrackState(const EVENT::TrackState& ts, double bField); - typedef std::vector< CaloHitData* > CaloHitDataVec ; - - /// define map type for storing the hits by layer - typedef std::map< int, CaloHitDataVec > CaloHitLayerMap ; - - /* compute the distance of the hit from a straight line - * defined by l = point + lambda * unitDir - */ - float computeDistanceFromLine( EVENT::CalorimeterHit* h, const dd4hep::rec::Vector3D& point, - const dd4hep::rec::Vector3D& unitDir) ; + double getHelixArcLength(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); - - /// return vector with hits that have the shortest distancefromStraightline for every layer - CaloHitDataVec findHitsClosestToLine( const CaloHitLayerMap& layerMap ) ; + double getHelixLengthAlongZ(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); + double getHelixNRevolutions(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); - /// compute the flight length of the particle from the IP to the calorimeter - float computeFlightLength( EVENT::Track* trk) ; + double getTPCOuterR(); - /// compute the flight length of the particle between the two given track states - float computeFlightLength(const EVENT::TrackState* ts0, const EVENT::TrackState* ts1 ) ; + EVENT::TrackerHit* getSETHit(EVENT::Track* track, double tpcOuterR); - /// helper function to get the layer of the calo hit - int layer( EVENT::CalorimeterHit* h ) ; + std::vector selectFrankEcalHits( EVENT::Cluster* cluster, EVENT::Track* track, int maxEcalLayer, double bField ); + // This is not guarantied to work 100% of times, + // but I didn't find a better way to collect subTracks... + // 1) Always add main track as a first subtrack + // 2) If nTPCHits == nHits of SubTrack0 +-1 assume subTrack0 stores TPC hits + // Deviation on 1 hit may happen because of the SET and only God knows what other reasons.. + // This would mean that VXD+SIT subTrack is not stored in the track! + // So we skip initial subTrack0 with TPC hits which we have anyhow added with global Track + // and loop over remaining subTracks + // 3) If 2) is false and nTPCHits == nHits of SubTrack1 +-1 assume subTrack1 stores TPC hits + // This would mean that subTrack0 stores VXD+SIT hits + // So we skip both VXD+SIT and 1st TPC subTracks which we have anyhow added with global Track + // and loop over remaining subTracks + std::vector getSubTracks(EVENT::Track* track); - /// helper function to check if this is an Ecal hit - bool isEcal( EVENT::CalorimeterHit* h ) ; + // return vector of trackStates for the IP, every hit in the whole track and ECAL if extrapolate to the ECAL option is chosen + // track state are sorted along the helix + std::vector getTrackStatesPerHit(std::vector tracks, MarlinTrk::IMarlinTrkSystem* trkSystem, bool extrapolateToEcal, double bField); + double getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, double timeResolution); - /// string with calo type information - std::string caloTypeStr( EVENT::CalorimeterHit* h ) ; + double getTofFrankAvg( std::vector selectedHits, EVENT::Track* track, double timeResolution); + double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); - /* compute the average TOF from all given hits - correcting for time of - * flight from reference point - */ - std::pair computeTOFEstimator( const CaloHitDataVec& chv ) ; - } diff --git a/TimeOfFlight/scripts/draw_cluster_time.C b/TimeOfFlight/scripts/draw_cluster_time.C deleted file mode 100644 index c75a6c08..00000000 --- a/TimeOfFlight/scripts/draw_cluster_time.C +++ /dev/null @@ -1,99 +0,0 @@ -#include "TCanvas.h" -#include "TROOT.h" -#include "TGraphErrors.h" -#include "TF1.h" -#include "TLegend.h" -#include "TArrow.h" -#include "TLatex.h" - -#include - -void draw_cluster_time(const std::string& histName){ - - std::string treeName("MyTOFPlots/") ; - treeName += histName ; - - - TFile fka("aida_file_kaon_random_REC.root"); - TH1* tka = (TH1*) fka.Get( treeName.c_str() ) ; - - TFile fpi("aida_file_pion_random_REC.root") ; - TH1* tpi = (TH1*) fpi.Get( treeName.c_str() ) ; - - TFile fmu("aida_file_muon_random_REC.root") ; - TH1* tmu = (TH1*) fmu.Get( treeName.c_str() ) ; - - TFile fpr("aida_file_proton_random_REC.root") ; - TH1* tpr = (TH1*) fpr.Get( treeName.c_str() ) ; - - TFile fel("aida_file_electron_random_REC.root") ; - TH1* tel = (TH1*) fel.Get( treeName.c_str() ) ; - - - - tka->SetLineColor( kRed ) ; - tpi->SetLineColor( kGreen ) ; - tmu->SetLineColor( kBlue ) ; - tpr->SetLineColor( kYellow+1 ) ; - tel->SetLineColor( kOrange ) ; - - tka->SetMarkerColor( kRed ) ; - tpi->SetMarkerColor( kGreen ) ; - tmu->SetMarkerColor( kBlue ) ; - tpr->SetMarkerColor( kYellow+1 ) ; - tel->SetMarkerColor( kOrange ) ; - - - TLegend leg(.1,.75,.2,.9,""); - leg.SetFillColor(0); - tmu->SetFillColor(0); - tpi->SetFillColor(0); - tka->SetFillColor(0); - tpr->SetFillColor(0); - tel->SetFillColor(0); - - tmu->Draw() ; - tka->Draw("same") ; - tpi->Draw("same") ; - tel->Draw("same") ; - tpr->Draw("same") ; - - TLegendEntry * l = leg.AddEntry(&fka,"kaons","l") ; - l->SetLineColor( kRed ) ; - l->SetMarkerColor( kRed ) ; - - l = leg.AddEntry(&fel,"electrons","l"); - l->SetLineColor( kOrange ); - l->SetMarkerColor( kOrange ); - - l = leg.AddEntry(&fpi,"pions","l") ; - l->SetLineColor( kGreen ) ; - l->SetMarkerColor( kGreen ) ; - - l = leg.AddEntry(&fmu,"muons","l") ; - l->SetLineColor( kBlue ) ; - l->SetMarkerColor( kBlue ) ; - - l = leg.AddEntry(&fpr,"protons","l") ; - l->SetLineColor( kYellow+1 ) ; - l->SetMarkerColor( kYellow+1 ) ; - - // l = leg.AddEntry(&fph,"photons","l") ; - // l->SetLineColor( kViolet ) ; - - leg.DrawClone("Same"); - - std::string filename( histName ) ; - filename += ".png" ; - - gPad->Print( filename.c_str() ) ; - - // gPad->WaitPrimitive() ; - -// fka.Close() ; -// fpi.Close() ; -// fmu.Close() ; -// fpr.Close() ; -// fel.Close() ; -// fph.Close() ; -} diff --git a/TimeOfFlight/scripts/draw_tofestimators.C b/TimeOfFlight/scripts/draw_tofestimators.C deleted file mode 100644 index 0f69c8b0..00000000 --- a/TimeOfFlight/scripts/draw_tofestimators.C +++ /dev/null @@ -1,70 +0,0 @@ -#include "TCanvas.h" -#include "TROOT.h" -#include "TGraphErrors.h" -#include "TF1.h" -#include "TLegend.h" -#include "TArrow.h" -#include "TLatex.h" - -#include - -/* - -.x ../TimeOfFlight/scripts/draw_tofestimators.C("./tof_estimators.root", "TOFEstimators50ps" , "hbetaFirstHitsChrg" ) -.x ../TimeOfFlight/scripts/draw_tofestimators.C("./tof_estimators.root", "TOFEstimators50ps" , "hbetaCloseHitsChrg" ) -.x ../TimeOfFlight/scripts/draw_tofestimators.C("./tof_estimators.root", "TOFEstimators50ps" , "hbetaFirstHitsNeut" ) -.x ../TimeOfFlight/scripts/draw_tofestimators.C("./tof_estimators.root", "TOFEstimators50ps" , "hbetaCloseHitsNeut" ) -.x ../TimeOfFlight/scripts/draw_tofestimators.C("./tof_estimators.root", "TOFEstimators50ps" , "hbetaLastTrkHit" ) - - */ - - -void draw_tofestimators(const std::string& fileName, - std::string treeName, - const std::string& histName){ - - std::string outFile = treeName + "_" + histName + ".png" ; - - - treeName += "/" ; - treeName += histName ; - - - - TFile file( fileName.c_str() ); - - TH2* h = (TH2*) file.Get( treeName.c_str() ) ; - - std::cout << " got histo " << treeName << " : h = " << h << std::endl ; - - if( h != 0 ) - h->Draw() ; - - - TF1 *fe = new TF1("fe","x/sqrt(x*x+0.0005109989461*0.0005109989461)",0.1,10.); - TF1 *fmu = new TF1("fmu","x/sqrt(x*x+0.105658*0.105658)",0.1,10.); - TF1 *fpi = new TF1("fpi","x/sqrt(x*x+0.13957018*0.13957018)",0.1,10.); - TF1 *fk = new TF1("fk","x/sqrt(x*x+0.493677*0.493677)",0.1,10.); - TF1 *fp = new TF1("fp","x/sqrt(x*x+0.9382720813*0.9382720813)",0.1,10.); - - fe->SetLineColor( kBlack ) ; - fmu->SetLineColor( kBlack ) ; - fpi->SetLineColor( kBlack ) ; - fk->SetLineColor( kBlack ) ; - fp->SetLineColor( kBlack ) ; - - fe-> SetLineWidth( 1 ) ; - fmu->SetLineWidth( 1 ) ; - fpi->SetLineWidth( 1 ) ; - fk-> SetLineWidth( 1 ) ; - fp-> SetLineWidth( 1 ) ; - - fe->Draw("same") ; - fmu->Draw("same") ; - fpi->Draw("same") ; - fk->Draw("same") ; - fp->Draw("same") ; - - gPad->Print( outFile.c_str() ) ; - -} diff --git a/TimeOfFlight/scripts/lcio_particle_gun.py b/TimeOfFlight/scripts/lcio_particle_gun.py deleted file mode 100644 index c0d51e20..00000000 --- a/TimeOfFlight/scripts/lcio_particle_gun.py +++ /dev/null @@ -1,192 +0,0 @@ -##################################### -# -# simple script to create lcio files with single particle -# events - modify as needed -# @author F.Gaede, DESY -# @date 1/07/2014 -# -# initialize environment: -# export PYTHONPATH=${LCIO}/src/python:${ROOTSYS}/lib -# -##################################### -import math -import random -from array import array - -# --- LCIO dependencies --- -from pyLCIO import UTIL, EVENT, IMPL, IO, IOIMPL -import sys - - -if len( sys.argv ) != 4: - print "usage: python lcio_particle_gun.py particle_type momentum polar_angle" - exit(0) - -#---- read input parameters -particle_type=sys.argv[1] -#------ one of 'muon','pion','kaon' - -momentum=sys.argv[2] - -polar_angle=sys.argv[3] - - -#---- number of events per momentum bin ----- -nevt = 1000 - - - - -momenta = [ float(momentum) ] -theta = float(polar_angle)/180. * math.pi - - - -if particle_type=="muon": - pdg = 13 - mass = 0.105658 -elif particle_type=="pion": - pdg = -211 - mass = 0.13957018 -elif particle_type=="kaon": - pdg = -321 - mass = 0.493677 -elif particle_type=="electron": - pdg = 11 - mass = 0.0005109989461 -elif particle_type=="proton": - pdg = 2212 - mass = 0.9382720813 -elif particle_type=="photon": - pdg = 22 - mass = 0. -else: - print " unknown particle type " - exit(1) - - - -outfile = "mcparticles_"+particle_type+"_"+momentum+"GeV_"+polar_angle+"deg.slcio" - -#-------------------------------------------- - - -wrt = IOIMPL.LCFactory.getInstance().createLCWriter( ) - -wrt.open( outfile , EVENT.LCIO.WRITE_NEW ) - -print " opened outfile: " , outfile - -random.seed() - - -#========== particle properties =================== - -#momenta = [ 1. , 3., 5., 10., 15., 25., 50., 100. ] - - - -genstat = 1 -charge = -1. - -decayLen = 1.e32 -#================================================= - -# write a RunHeader -run = IMPL.LCRunHeaderImpl() -run.setRunNumber( 0 ) -run.setDetectorName("ILD_l4_v02") -run.parameters().setValue("Generator","${lcgeo}_DIR/examples/lcio_particle_gun.py") -run.parameters().setValue("PDG", pdg ) -run.parameters().setValue("Charge", charge ) -run.parameters().setValue("Mass", mass ) -wrt.writeRunHeader( run ) -#================================================ - - -for p in momenta: - - for j in range( 0, nevt ): - - col = IMPL.LCCollectionVec( EVENT.LCIO.MCPARTICLE ) - evt = IMPL.LCEventImpl() - - evt.setEventNumber( j ) - - evt.addCollection( col , "MCParticle" ) - -# phi = random.random() * math.pi * 2. - phi = math.pi / 4. ; - - energy = math.sqrt( mass*mass + p * p ) - - px = p * math.cos( phi ) * math.sin( theta ) - py = p * math.sin( phi ) * math.sin( theta ) - pz = p * math.cos( theta ) - - momentum = array('f',[ px, py, pz ] ) - - epx = decayLen * math.cos( phi ) * math.sin( theta ) - epy = decayLen * math.sin( phi ) * math.sin( theta ) - epz = decayLen * math.cos( theta ) - - endpoint = array('d',[ epx, epy, epz ] ) - - -#--------------- create MCParticle ------------------- - - mcp = IMPL.MCParticleImpl() - - mcp.setGeneratorStatus( genstat ) - mcp.setMass( mass ) - mcp.setPDG( pdg ) - mcp.setMomentum( momentum ) - mcp.setCharge( charge ) - - if( decayLen < 1.e9 ) : # arbitrary ... - mcp.setEndpoint( endpoint ) - - -#------------------------------------------------------- - - - -#------------------------------------------------------- - - - col.addElement( mcp ) - - wrt.writeEvent( evt ) - - -wrt.close() - - -# -# longer format: - use ".hepevt" -# - -# -# int ISTHEP; // status code -# int IDHEP; // PDG code -# int JMOHEP1; // first mother -# int JMOHEP2; // last mother -# int JDAHEP1; // first daughter -# int JDAHEP2; // last daughter -# double PHEP1; // px in GeV/c -# double PHEP2; // py in GeV/c -# double PHEP3; // pz in GeV/c -# double PHEP4; // energy in GeV -# double PHEP5; // mass in GeV/c**2 -# double VHEP1; // x vertex position in mm -# double VHEP2; // y vertex position in mm -# double VHEP3; // z vertex position in mm -# double VHEP4; // production time in mm/c -# -# inputFile >> ISTHEP >> IDHEP -# >> JMOHEP1 >> JMOHEP2 -# >> JDAHEP1 >> JDAHEP2 -# >> PHEP1 >> PHEP2 >> PHEP3 -# >> PHEP4 >> PHEP5 -# >> VHEP1 >> VHEP2 >> VHEP3 -# >> VHEP4; diff --git a/TimeOfFlight/scripts/make_root_plots.sh b/TimeOfFlight/scripts/make_root_plots.sh deleted file mode 100644 index 751620c0..00000000 --- a/TimeOfFlight/scripts/make_root_plots.sh +++ /dev/null @@ -1,26 +0,0 @@ - -# -# create root plots -# - -# root -b -q 'draw_cluster_time.C("beta_05hits_50ps")' -# root -b -q 'draw_cluster_time.C("beta_05perc_50ps")' -# root -b -q 'draw_cluster_time.C("beta_10hits_50ps")' -# root -b -q 'draw_cluster_time.C("beta_10perc_50ps")' -# root -b -q 'draw_cluster_time.C("beta_20hits_50ps")' -# root -b -q 'draw_cluster_time.C("beta_20perc_50ps")' - -root -b -q 'draw_cluster_time.C("clutime_0ps")' -root -b -q 'draw_cluster_time.C("time05perc_0ps")' -root -b -q 'draw_cluster_time.C("time10perc_0ps")' -root -b -q 'draw_cluster_time.C("time20perc_0ps")' -root -b -q 'draw_cluster_time.C("timeoffastesthit_0ps")' -root -b -q 'draw_cluster_time.C("time05hits_0ps")' -root -b -q 'draw_cluster_time.C("time10hits_0ps")' -root -b -q 'draw_cluster_time.C("time20hits_0ps")' -root -b -q 'draw_cluster_time.C("cor_ref_time05perc_0ps")' -root -b -q 'draw_cluster_time.C("cor_ref_time10perc_0ps")' -root -b -q 'draw_cluster_time.C("cor_ref_time20perc_0ps")' -root -b -q 'draw_cluster_time.C("cor_ref_time05hits_0ps")' -root -b -q 'draw_cluster_time.C("cor_ref_time10hits_0ps")' -root -b -q 'draw_cluster_time.C("cor_ref_time20hits_0ps")' diff --git a/TimeOfFlight/scripts/run_all_and_make_root_plots_random.sh b/TimeOfFlight/scripts/run_all_and_make_root_plots_random.sh deleted file mode 100644 index dbbea2c9..00000000 --- a/TimeOfFlight/scripts/run_all_and_make_root_plots_random.sh +++ /dev/null @@ -1,21 +0,0 @@ - -# -# run all Marlin random files + create root plots -# - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm11_RandomP.n010_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_electron_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm13_RandomP.n001_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_muon_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm211_RandomP.n010_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_pion_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm2212_RandomP.n010_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_proton_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm321_RandomP.n001_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_kaon_random_REC - -root -b -q 'draw_cluster_time.C("hBeta_5hitsvsMomentum_0ps")' -root -b -q 'draw_cluster_time.C("hBeta_5percvsMomentum_0ps")' -root -b -q 'draw_cluster_time.C("hBeta_10hitsvsMomentum_0ps")' -root -b -q 'draw_cluster_time.C("hBeta_10percvsMomentum_0ps")' -root -b -q 'draw_cluster_time.C("hBeta_20hitsvsMomentum_0ps")' -root -b -q 'draw_cluster_time.C("hBeta_20percvsMomentum_0ps")' diff --git a/TimeOfFlight/scripts/run_all_marlin.sh b/TimeOfFlight/scripts/run_all_marlin.sh deleted file mode 100644 index f6ab5133..00000000 --- a/TimeOfFlight/scripts/run_all_marlin.sh +++ /dev/null @@ -1,11 +0,0 @@ - -# -# create root files with time of flight plots for all particles -# - -Marlin mysteer.xml --global.LCIOInputFiles=../../data/mcparticles_proton_5GeV_60deg_REC.slcio --MyAIDAProcessor.FileName=aida_file_proton_5GeV_60deg_REC -Marlin mysteer.xml --global.LCIOInputFiles=../../data/mcparticles_kaon_5GeV_60deg_REC.slcio --MyAIDAProcessor.FileName=aida_file_kaon_5GeV_60deg_REC -Marlin mysteer.xml --global.LCIOInputFiles=../../data/mcparticles_muon_5GeV_60deg_REC.slcio --MyAIDAProcessor.FileName=aida_file_muon_5GeV_60deg_REC -Marlin mysteer.xml --global.LCIOInputFiles=../../data/mcparticles_pion_5GeV_60deg_REC.slcio --MyAIDAProcessor.FileName=aida_file_pion_5GeV_60deg_REC -Marlin mysteer.xml --global.LCIOInputFiles=../../data/mcparticles_electron_5GeV_60deg_REC.slcio --MyAIDAProcessor.FileName=aida_file_electron_5GeV_60deg_REC - diff --git a/TimeOfFlight/scripts/run_all_marlin_10GeV.sh b/TimeOfFlight/scripts/run_all_marlin_10GeV.sh deleted file mode 100644 index 111196c0..00000000 --- a/TimeOfFlight/scripts/run_all_marlin_10GeV.sh +++ /dev/null @@ -1,22 +0,0 @@ - -# -# create root files with time of flight plots for all particles with 5 GeV -# - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG13/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG13_MOM010GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_muon_10GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG321/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG321_MOM010GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_kaon_10GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG11/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG11_MOM010GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_electron_10GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG211/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG211_MOM010GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_pion_10GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG2212/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG2212_MOM010GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_proton_10GeV_REC - - -root -b -q 'draw_cluster_time.C("beta_05hits_10GeV_10ps")' -root -b -q 'draw_cluster_time.C("beta_05perc_10GeV_10ps")' -root -b -q 'draw_cluster_time.C("beta_10hits_10GeV_10ps")' -root -b -q 'draw_cluster_time.C("beta_10perc_10GeV_10ps")' -root -b -q 'draw_cluster_time.C("beta_20hits_10GeV_10ps")' -root -b -q 'draw_cluster_time.C("beta_20perc_10GeV_10ps")' diff --git a/TimeOfFlight/scripts/run_all_marlin_2GeV.sh b/TimeOfFlight/scripts/run_all_marlin_2GeV.sh deleted file mode 100644 index e5db5e96..00000000 --- a/TimeOfFlight/scripts/run_all_marlin_2GeV.sh +++ /dev/null @@ -1,22 +0,0 @@ - -# -# create root files with time of flight plots for all particles with 2 GeV -# - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG13/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG13_MOM002GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_muon_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG321/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG321_MOM002GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_kaon_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG11/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG11_MOM002GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_electron_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG211/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG211_MOM002GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_pion_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG2212/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG2212_MOM002GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_proton_02GeV_REC - - -root -b -q 'draw_cluster_time.C("beta_05hits_02GeV_0ps")' -root -b -q 'draw_cluster_time.C("beta_05perc_02GeV_0ps")' -root -b -q 'draw_cluster_time.C("beta_10hits_02GeV_0ps")' -root -b -q 'draw_cluster_time.C("beta_10perc_02GeV_0ps")' -root -b -q 'draw_cluster_time.C("beta_20hits_02GeV_0ps")' -root -b -q 'draw_cluster_time.C("beta_20perc_02GeV_0ps")' diff --git a/TimeOfFlight/scripts/run_all_marlin_5GeV.sh b/TimeOfFlight/scripts/run_all_marlin_5GeV.sh deleted file mode 100644 index 81c08f93..00000000 --- a/TimeOfFlight/scripts/run_all_marlin_5GeV.sh +++ /dev/null @@ -1,22 +0,0 @@ - -# -# create root files with time of flight plots for all particles with 5 GeV -# - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG13/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG13_MOM005GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_muon_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG321/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG321_MOM005GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_kaon_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG11/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG11_MOM005GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_electron_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG211/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG211_MOM005GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_pion_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG2212/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG2212_MOM005GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_proton_05GeV_REC - - -root -b -q 'draw_cluster_time.C("beta_05hits_05GeV_50ps")' -root -b -q 'draw_cluster_time.C("beta_05perc_05GeV_50ps")' -root -b -q 'draw_cluster_time.C("beta_10hits_05GeV_50ps")' -root -b -q 'draw_cluster_time.C("beta_10perc_05GeV_50ps")' -root -b -q 'draw_cluster_time.C("beta_20hits_05GeV_50ps")' -root -b -q 'draw_cluster_time.C("beta_20perc_05GeV_50ps")' diff --git a/TimeOfFlight/scripts/run_all_marlin_v1.sh b/TimeOfFlight/scripts/run_all_marlin_v1.sh deleted file mode 100644 index cf192050..00000000 --- a/TimeOfFlight/scripts/run_all_marlin_v1.sh +++ /dev/null @@ -1,45 +0,0 @@ - -# -# create root files with time of flight plots for all particles -# - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG13/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG13_MOM002GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_muon_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG13/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG13_MOM005GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_muon_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG13/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG13_MOM010GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_muon_10GeV_REC - - - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG321/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG321_MOM002GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_kaon_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG321/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG321_MOM005GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_kaon_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG321/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG321_MOM010GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_kaon_10GeV_REC - - - - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG11/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG11_MOM002GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_electron_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG11/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG11_MOM005GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_electron_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG11/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG11_MOM010GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_electron_10GeV_REC - - - - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG211/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG211_MOM002GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_pion_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG211/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG211_MOM005GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_pion_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG211/ILD_l4_v02/v01-19-04_lcgeo/u014/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG211_MOM010GeV.n001_1.d_rec_u014.slcio --MyAIDAProcessor.FileName=aida_file_pion_10GeV_REC - - - - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG2212/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG2212_MOM002GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_proton_02GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG2212/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG2212_MOM005GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_proton_05GeV_REC - - Marlin mysteer.xml --global.LCIOInputFiles=/pnfs/desy.de/ilc/prod/ilc/mc-opt/ild/rec/calib/single_PDG2212/ILD_l4_v02/v01-19-04_lcgeo/u017/rv01-19-04_lcgeo.sv01-19-04_lcgeo.mILD_l4_v02.Pmcparticles_PDG2212_MOM010GeV.n001_1.d_rec_u017.slcio --MyAIDAProcessor.FileName=aida_file_proton_10GeV_REC diff --git a/TimeOfFlight/scripts/run_all_marlin_v2.sh b/TimeOfFlight/scripts/run_all_marlin_v2.sh deleted file mode 100644 index d99c0c33..00000000 --- a/TimeOfFlight/scripts/run_all_marlin_v2.sh +++ /dev/null @@ -1,12 +0,0 @@ - - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm11_RandomP.n010_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_electron_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm13_RandomP.n001_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_muon_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm211_RandomP.n010_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_pion_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm2212_RandomP.n010_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_proton_random_REC - -Marlin mysteer.xml --global.LCIOInputFiles=/nfs/dust/ilc/group/ild/DD4hep_data/REC/gun/sv01-19-05.mILD_l5_v02.Pmcparticles_PDGpm321_RandomP.n001_1.d_sim_u033_REC.slcio --MyAIDAProcessor.FileName=aida_file_kaon_random_REC - diff --git a/TimeOfFlight/scripts/tofestimators.xml b/TimeOfFlight/scripts/tofestimators.xml deleted file mode 100644 index cd6af1cc..00000000 --- a/TimeOfFlight/scripts/tofestimators.xml +++ /dev/null @@ -1,129 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - simjob.slcio - - - - - - - MESSAGE DEBUG - - - - - - - - - 1 - - tof_estimators - - root - - - - - - - - - 10 - - PandoraPFOs - - 0 - - - - - idr - - - - - 10 - 0 - dev - - - - 10 - PandoraPFOs - 10. - - - - 10 - PandoraPFOs - 50 - - - - MCParticle - PandoraPFOs - - - - - - - - - - - - - - - tof_estimators_output.slcio - - WRITE_NEW - - - - - - - - - - - /cvmfs/ilc.desy.de/sw/x86_64_gcc82_centos7/v02-02/lcgeo/v00-16-06/ILD/compact/${DetectorModel}/${DetectorModel}.xml - - - - - diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index 03c84e8d..9f34c799 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -1,638 +1,187 @@ #include "TOFEstimators.h" #include "TOFUtils.h" -#include -#include - -#include -#include -#include -#include -#include -#include -#include - #include "EVENT/LCCollection.h" -#include "EVENT/ReconstructedParticle.h" -using EVENT::LCCollection, EVENT::ReconstructedParticle; - -#include "DD4hep/Detector.h" -#include "DD4hep/DD4hepUnits.h" - -#include "TGraphErrors.h" -#include "TF1.h" -#include "CLHEP/Units/PhysicalConstants.h" - -#include "HelixClass.h" -#include "marlinutil/CalorimeterHitType.h" -#include - -// ----- include for verbosity dependend logging --------- -#include "marlin/VerbosityLevels.h" +#include "UTIL/PIDHandler.h" +#include "marlin/Global.h" #include "marlin/ProcessorEventSeeder.h" -#include -#include - - -#include -#include -#include -#include - -#include -#include - -using namespace lcio ; -using namespace marlin ; -using namespace TOFUtils ; +#include "marlin/VerbosityLevels.h" +#include "marlinutil/GeometryUtil.h" +#include "MarlinTrk/Factory.h" +#include "EVENT/SimTrackerHit.h" +#include "UTIL/LCRelationNavigator.h" +#include "CLHEP/Random/Randomize.h" + +using namespace TOFUtils; +using std::vector; +using std::string; +using EVENT::LCCollection; +using EVENT::ReconstructedParticle; +using EVENT::TrackerHit; +using EVENT::Track; +using EVENT::SimTrackerHit; +using EVENT::Cluster; +using EVENT::CalorimeterHit; +using EVENT::TrackState; +using EVENT::LCObject; +using UTIL::LCRelationNavigator; +using CLHEP::RandGauss; +using dd4hep::rec::Vector3D; TOFEstimators aTOFEstimators ; -TOFEstimators::TOFEstimators() : Processor("TOFEstimators") { - _description = "TOFEstimators compute some estimators for the time of flight from calorimeter hits" ; +TOFEstimators::TOFEstimators() : marlin::Processor("TOFEstimators") { + _description = "TOFEstimators processor computes momentum harmonic mean, track length \ + and time-of-flight of the chosen ReconstructedParticle to the \ + specified end point (SET hit or Ecal surface). To be used for a further particle ID"; - // register steering parameters: name, description, class-variable, default value - registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, + registerInputCollection(LCIO::RECONSTRUCTEDPARTICLE, "ReconstructedParticleCollection", "Name of the ReconstructedParticle collection", - _colNamePFO, + _pfoCollectionName, std::string("PandoraPFOs") ); - registerProcessorParameter( "MaxLayerNumber", - "Use only calorimeter hits up to MaxLayerNumber in TOF estimators", - _maxLayerNum, - int(100) ); - - registerProcessorParameter( "CylRadius", - "Cut-off hits further away from the shower core than this raduius for TOF calculation", - _cylRadiusCut, - double(5.) ); - - - registerProcessorParameter( "TimeResolution", - "Assumed time resolution per hit in ps", - _resolution, - float(0.) ); - - registerProcessorParameter( "ProcessorVersion", - "Legacy or new TOFEstimator", - _procVersion, - std::string("idr") ); + registerProcessorParameter("ExtrapolateToEcal", + "If true, track is extrapolated to the Ecal surface for track length calculation, \ + time of flight estimated using Ecal hits. If false, track length calculated to the last tracker hit.\ + Time of flight estimated using SET hit if exists.", + _extrapolateToEcal, + bool(true)); + + registerProcessorParameter("TofMethod", + "name of the algorithm which estimates time of flight\ + to the Ecal surface based on Ecal hits time information.\ + Available options are: closest, frankAvg, frankFit.\ + In case of _extrapolateToEcal==false is ignored", + _tofMethod, + std::string("closest") ); + + registerProcessorParameter("TimeResolution", + "Time resolution of individual SET strips or Ecal hits in ps", + _timeResolution, + double(0.)); + + registerProcessorParameter("MaxEcalLayer", + "Time of flight is calculated using Ecal hits only up to MaxLayer", + _maxEcalLayer, + int(10) ); } -void TOFEstimators::init() { - Global::EVENTSEEDER->registerProcessor(this); - if ( _procVersion == "idr" ){ - // initialize gsl random generator - _rng = gsl_rng_alloc(gsl_rng_ranlxs2); - _TOFNames = {"TOFFirstHit", "TOFClosestHits", "TOFClosestHitsError", "TOFFlightLength", "TOFLastTrkHit" , "TOFLastTrkHitFlightLength"}; - } - else if (_procVersion == "dev"){ - _smearing.param( std::normal_distribution::param_type(0., _resolution/1000.) ); - _TOFNames = {"TOFClosest", "TOFFastest", "TOFCylFit", "TOFClosestFit", "FlightLength", "MomAtCalo"}; - - const auto& detector = dd4hep::Detector::getInstance(); - detector.field().magneticField({0., 0., 0.}, _bField); - } - else { - throw EVENT::Exception(std::string("Invalid ProcessorVersion parameter passed: \'") + _procVersion + std::string("\'\n Viable options are idr (default), dev")); - } -} - - -void TOFEstimators::processEvent( LCEvent * evt ) { - - if ( _procVersion == "idr" ){ - // use the global Marlin random seed for this processor - gsl_rng_set( _rng, Global::EVENTSEEDER->getSeed(this) ) ; - - streamlog_out(DEBUG ) << "seed set to " << Global::EVENTSEEDER->getSeed(this) << std::endl; - - streamlog_out(DEBUG2) << " process event: " << evt->getEventNumber() - << " in run: " << evt->getRunNumber() << std::endl ; - - - // get the PFO collection from the event if it exists - LCCollection* colPFO = nullptr ; - - try{ - colPFO = evt->getCollection( _colNamePFO ) ; - } - catch(lcio::Exception&){ - streamlog_out( DEBUG6 ) << " collection " << _colNamePFO - << " not found in event - nothing to do ... " << std::endl ; - } - - if( colPFO->getTypeName() != LCIO::RECONSTRUCTEDPARTICLE ) { - - streamlog_out( ERROR ) << " collection " << _colNamePFO - << " not of type LCIO::RECONSTRUCTEDPARTICLE " << std::endl ; - - colPFO = nullptr ; - } - - if( colPFO != nullptr ){ - - - - PIDHandler pidh( colPFO ); - int algoID = pidh.addAlgorithm( name() , _TOFNames); - - - int nPFO = colPFO->getNumberOfElements() ; - - - for(int i=0; i< nPFO ; ++i){ - - ReconstructedParticle* pfo = dynamic_cast( colPFO->getElementAt( i ) ) ; - - bool isCharged = false ; - - if( pfo->getClusters().size() != 1 ){ - - streamlog_out( DEBUG1 ) << " ignore particle w/ cluster number other than one: " << *pfo << std::endl ; - continue ; - } - - if( fabs( pfo->getCharge() ) < 0.1 && pfo->getTracks().size() == 0 ) {;} - else if ( fabs( pfo->getCharge() ) > 0.1 && pfo->getTracks().size() == 1 ) { - isCharged = true ;} - else { - streamlog_out( DEBUG1 ) << " ignore particle w/ track number other than zero or one: " << *pfo << std::endl ; - continue ; - } - - streamlog_out( DEBUG1 ) << " --- compute TOF estimators for particle : " << *pfo << std::endl ; - - - - // ======= use only Ecal hits (requires the CalorimeterHitType to be set in the digitizer ) - // with time information ( > 1 ps) and layer <= max layer - - int maxLayerNum = _maxLayerNum ; - std::function selectHits = [maxLayerNum](CalorimeterHit* h){ - - return ( isEcal( h ) && - h->getTime() > 1.e-3 && - layer( h ) <= maxLayerNum ) ; - } ; - - // ------------------------------------------------------------------------------------------- - - Cluster* clu = pfo->getClusters()[0] ; - const CalorimeterHitVec& cluhv = clu->getCalorimeterHits() ; - - // create vectors of extended handle objects for relevant calorimeter hits - // one w/ unique_ptr for memory handling - - CaloHitUPtrVec uniqueVec ; - uniqueVec.reserve( cluhv.size() ) ; - - CaloHitDataVec caloHitVec ; - caloHitVec.reserve( cluhv.size() ) ; - //------------------------------------------------------------------------------- - - CaloHitLayerMap layerMap ; - - for( auto* clh : cluhv ){ - - if( selectHits( clh ) ){ - - uniqueVec.push_back( std::unique_ptr( new CaloHitData( clh) ) ) ; - - CaloHitData* ch = uniqueVec.back().get() ; - - caloHitVec.push_back( ch ) ; - - ch->layer = layer( ch->lcioHit ) ; - ch->timeResolution = _resolution ; - - ch->smearedTime = ( _resolution > 0. ? - ch->lcioHit->getTime() + gsl_ran_gaussian( _rng, _resolution / 1000. ) : // convert ps to ns - ch->lcioHit->getTime() ) ; - - ch->distanceFromIP = dd4hep::rec::Vector3D( clh->getPosition() ).r() ; - - - layerMap[ ch->layer ].push_back( ch ) ; - } - } - - - if( layerMap.empty() ) { - - streamlog_out( DEBUG1 ) << " --- not suitable Ecal hits found for particle " << std::endl ; - continue ; - } - - // streamlog_out( DEBUG ) << " --- map with hits per layer : " << std::endl ; - // for( auto m : layerMap ){ - // streamlog_out( DEBUG ) << " ----- layer " << m.first << " : " << std::endl ; - // for( auto ch : m.second ) - // streamlog_out( DEBUG ) << " " << caloTypeStr( ch->lcioHit ) << std::endl ; - // } - - - // --- define reference point: track state at calo for charged - hit closest to IP for neutral - // and direction of straight line - either from IP or from track state at calo - - dd4hep::rec::Vector3D refPoint ; - dd4hep::rec::Vector3D unitDir ; - float flightLength = 0. ; - - TrackerHit* lastTrackerHit = nullptr ; - float flightLengthTrkHit = 0. ; - - if( isCharged ){ - - Track* trk = pfo->getTracks()[0] ; - const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ; - - refPoint = tscalo->getReferencePoint() ; - - float tanL = tscalo->getTanLambda() ; - float theta = atan( 1. / tanL ) ; - - unitDir = dd4hep::rec::Vector3D( 1. , tscalo->getPhi() , theta , dd4hep::rec::Vector3D::spherical ) ; - - flightLength = computeFlightLength( trk ) ; - - - // also store time and flight length of last tracker hit - const TrackState* tsIP = trk->getTrackState( TrackState::AtIP ) ; - const TrackState* tslh = trk->getTrackState( TrackState::AtLastHit ) ; - - lastTrackerHit = trk->getTrackerHits().back() ; - - flightLengthTrkHit = computeFlightLength( tsIP , tslh ) ; - - - #if 1 - if( lastTrackerHit->getTime() > 1e-3 ) { - dd4hep::rec::Vector3D rpLH = tslh->getReferencePoint() ; - dd4hep::rec::Vector3D lhp = lastTrackerHit->getPosition() ; - streamlog_out( DEBUG3 ) << " *************** referenece point calo : " << refPoint << std::endl ; - streamlog_out( DEBUG3 ) << " *************** referenece point last hit : " << rpLH << std::endl ; - streamlog_out( DEBUG3 ) << " *************** poistion last hit : " << lhp << std::endl ; - streamlog_out( DEBUG3 ) << " distance hit-trkstate: " << (rpLH - lhp ).r() << " -- distance calo/last hit ref points : " << (refPoint-rpLH).r() << std::endl ; - streamlog_out( DEBUG3 ) << " flight lengths: " << flightLength << " - " << flightLengthTrkHit << " -- diff " << flightLength - flightLengthTrkHit << - " time diff: " << (flightLength - flightLengthTrkHit) / CLHEP::c_light << std::endl ; - streamlog_out( DEBUG3 ) << " track state : " << *tslh << std::endl ; - streamlog_out( DEBUG3 ) << " last hit : " << *lastTrackerHit << std::endl ; - } - #endif - - - } else { // neutral particle - - CaloHitDataVec chv = layerMap.begin()->second ; // only look in first layer w/ hits - - CaloHitData* closestHit = - *min_element( chv.begin() , chv.end () , - [](CaloHitData* c0, CaloHitData* c1 ){ return c0->distanceFromIP < c1->distanceFromIP ; } - ) ; - - refPoint = closestHit->lcioHit->getPosition() ; - - flightLength = refPoint.r() ; - - - dd4hep::rec::Vector3D cluPos = clu->getPosition() ; - - unitDir = cluPos.unit() ; - - } - - streamlog_out( DEBUG2 ) << " ----- use reference point for TOF : " << refPoint << std::endl ; +void TOFEstimators::init(){ + marlin::Global::EVENTSEEDER->registerProcessor(this); - streamlog_out( DEBUG ) << " ----- calorimeter hits considered for estimators : " << std::endl ; + _outputParNames = {"momentumHM", "trackLength", "timeOfFlight"}; + _bField = MarlinUtil::getBzAtOrigin(); + _tpcOuterR = getTPCOuterR(); - - // ------ loop again over hits and fill missing data - - for( auto ch : caloHitVec ){ - - CalorimeterHit* calohit = ch->lcioHit ; - - dd4hep::rec::Vector3D pos = ch->lcioHit->getPosition() ; - - ch->distanceFromReferencePoint = ( pos - refPoint ).r() ; - - ch->distancefromStraightline = computeDistanceFromLine( calohit, refPoint, unitDir ) ; - - streamlog_out( DEBUG ) << " ----- " << caloTypeStr( calohit ) - << " -- " << ch->toString() << std::endl ; - } - - // ---- now get hits that are closest to the extrapolated line - CaloHitDataVec tofHits = findHitsClosestToLine( layerMap ) ; - - - streamlog_out( DEBUG ) << " ***** hits used for the TOF estimator : " << std::endl ; - for( auto ch : tofHits ){ - streamlog_out( DEBUG ) << " ----- " << ch->toString() << std::endl ; - } - - auto t_dt = computeTOFEstimator( tofHits ) ; - - float tof_fh = tofHits[0]->smearedTime - tofHits[0]->distanceFromReferencePoint / CLHEP::c_light ; - - streamlog_out( DEBUG2 ) << " #### tof ( first ) : " << tof_fh << " +/- " << 0 << std::endl ; - streamlog_out( DEBUG2 ) << " #### tof ( straight line ) : " << t_dt.first << " +/- " << t_dt.second << std::endl ; - - - float trkHitTime = ( lastTrackerHit ? lastTrackerHit->getTime() : 0. ) ; - - - FloatVec TOF_params = { tof_fh, - t_dt.first, t_dt.second, - flightLength , trkHitTime , flightLengthTrkHit } ; - - pidh.setParticleID( pfo , 0, 0 , 0.0 , algoID, TOF_params ); - - - //========================================================================================= - - } - - } - } - else if (_procVersion == "dev"){ - _generator.seed( Global::EVENTSEEDER->getSeed(this) ); - - LCCollection* colPFO = evt->getCollection(_colNamePFO); - PIDHandler handler( colPFO ); - int algoID = handler.addAlgorithm( name() , _TOFNames); - - for (int i=0; igetNumberOfElements(); ++i){ - ReconstructedParticle* pfo = dynamic_cast ( colPFO->getElementAt(i) ); - int nClusters = pfo->getClusters().size(); - int nTracks = pfo->getTracks().size(); - - // Only simple cases of PFOs - if( nClusters != 1 || nTracks != 1) continue; - - const Track* track = pfo->getTracks()[0]; - const Cluster* cluster = pfo->getClusters()[0]; - - float momAtCalo = getMomAtCalo(track).r(); - float flightLength = getFlightLength(track); - float tofClosest = getTOFClosest(track, cluster); - float tofFastest = getTOFFastest(track, cluster); - float tofCylFit = getTOFCylFit(track, cluster); - float tofClosestFit = getTOFClosestFit(track, cluster); - - const std::vector results{tofClosest, tofFastest, tofCylFit, tofClosestFit, flightLength, momAtCalo}; - handler.setParticleID(pfo , 0, 0, 0., algoID, results); - } - } + _trkSystem = MarlinTrk::Factory::createMarlinTrkSystem("DDKalTest", nullptr, ""); + _trkSystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useQMS, true); + _trkSystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::usedEdx, true); + _trkSystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useSmoothing, true); + _trkSystem->init(); } -void TOFEstimators::check( LCEvent *evt) { - if( _procVersion == "idr"){ - streamlog_out( DEBUG ) << " --- check called ! " << std::endl ; - - - // create some histograms with beta vs momentum for charged particles - - if( isFirstEvent() ){ - - // this creates a directory for this processor .... - AIDAProcessor::tree( this ) ; - - _h.resize(5) ; - int nBins = 100 ; - _h[0] = new TH2F( "hbetaFirstHitsChrg", "beta vs momentum - first hit - charged", nBins, .1 , 10., nBins, 0.93 , 1.03 ) ; - _h[1] = new TH2F( "hbetaCloseHitsChrg", "beta vs momentum - closest hits - charged", nBins, .1 , 10., nBins, 0.93 , 1.03 ) ; - - _h[2] = new TH2F( "hbetaFirstHitsNeut", "beta vs momentum - first hit - neutral", nBins, .1 , 10., nBins, 0.93 , 1.03 ) ; - _h[3] = new TH2F( "hbetaCloseHitsNeut", "beta vs momentum - closest hits - neutral", nBins, .1 , 10., nBins, 0.93 , 1.03 ) ; - - _h[4] = new TH2F( "hbetaLastTrkHit", "beta vs momentum - last tracker hit", nBins, .1 , 10., nBins, 0.93 , 1.03 ) ; - - - } - - // get the PFO collection from the event if it exists - LCCollection* colPFO = nullptr ; - try{ colPFO = evt->getCollection( _colNamePFO ) ; } catch(lcio::Exception&){} - - if( colPFO != nullptr ){ - - PIDHandler pidh( colPFO ); - int algoID = pidh.getAlgorithmID( name() ); - int tof_firsthit = pidh.getParameterIndex(algoID,"TOFFirstHit") ; - int tof_closest = pidh.getParameterIndex(algoID,"TOFClosestHits") ; - int tof_length = pidh.getParameterIndex(algoID,"TOFFlightLength") ; - int tof_trkhit = pidh.getParameterIndex(algoID,"TOFLastTrkHit") ; - int tof_trk_len = pidh.getParameterIndex(algoID,"TOFLastTrkHitFlightLength") ; - - int nPFO = colPFO->getNumberOfElements() ; - - for( int i=0 ; i< nPFO ; ++i){ - - ReconstructedParticle* pfo = dynamic_cast( colPFO->getElementAt( i ) ) ; - - const ParticleID& tofPID = pidh.getParticleID( pfo , algoID ) ; - - const FloatVec& tofParams = tofPID.getParameters() ; - - streamlog_out( DEBUG ) << " **** found TOF parameters for pfo w/ size " << tofParams.size() << std::endl ; - - if( !tofParams.empty() ){ +void TOFEstimators::processEvent(LCEvent * evt){ + RandGauss::setTheSeed( marlin::Global::EVENTSEEDER->getSeed(this) ); + ++_nEvent; + streamlog_out(MESSAGE)<<"******Event****** "<<_nEvent<getMomentum() ; - double momentum = sqrt( mom[0] * mom[0] + mom[1] * mom[1] + mom[2] * mom[2] ) ; - double length = tofParams[ tof_length ] ; + LCCollection* pfos = evt->getCollection(_pfoCollectionName); + LCRelationNavigator navigatorSET = LCRelationNavigator( evt->getCollection("SETSpacePointRelations") ); - double beta_fh = ( length / tofParams[ tof_firsthit] ) / CLHEP::c_light ; - double beta_ch = ( length / tofParams[ tof_closest ] ) / CLHEP::c_light ; + PIDHandler pidHandler( pfos ); + int algoID = pidHandler.addAlgorithm( name(), _outputParNames ); - if( abs( pfo->getCharge() ) > 0.5 ) { - _h[0]->Fill( momentum , beta_fh ); - _h[1]->Fill( momentum , beta_ch ); - } else { - _h[2]->Fill( momentum , beta_fh ); - _h[3]->Fill( momentum , beta_ch ); - } + for (int i=0; igetNumberOfElements(); ++i){ + ReconstructedParticle* pfo = static_cast ( pfos->getElementAt(i) ); - if( tofParams[ tof_trk_len ] > 1.e-3 ){ // if TOF from last tracker hit has been set + int nClusters = pfo->getClusters().size(); + int nTracks = pfo->getTracks().size(); - double beta = ( tofParams[ tof_trk_len ] / tofParams[ tof_trkhit ] ) / CLHEP::c_light ; - - _h[4]->Fill( momentum , beta ); - } - - } + if( nClusters != 1 || nTracks != 1){ + // Analyze only simple pfos. Otherwise write dummy zeros + vector results{0., 0., 0.}; + pidHandler.setParticleID(pfo , 0, 0, 0., algoID, results); + continue; } - } - } -} - - -void TOFEstimators::end(){ - if( _procVersion == "idr") gsl_rng_free( _rng ); -} - - -dd4hep::rec::Vector3D TOFEstimators::getMomAtCalo(const Track* track){ - const TrackState* ts = track->getTrackState(TrackState::AtCalorimeter); - double phi = ts->getPhi(); - double d0 = ts->getD0(); - double z0 = ts->getZ0(); - double omega = ts->getOmega(); - double tanL = ts->getTanLambda(); - - HelixClass helix; - helix.Initialize_Canonical(phi, d0, z0, omega, tanL, _bField[2]/dd4hep::tesla); - return helix.getMomentum(); -} - - -double TOFEstimators::getFlightLength(const Track* track){ - const TrackState* ts = track->getTrackState(TrackState::AtIP); - double phiIp = ts->getPhi(); - ts = track->getTrackState(TrackState::AtCalorimeter); - double phiCalo = ts->getPhi(); - double omegaCalo = ts->getOmega(); - double tanLCalo = ts->getTanLambda(); - - double dPhi = std::abs(phiIp - phiCalo); - - // if dPhi > PI assume phi is flipped - if (dPhi > M_PI) dPhi = 2*M_PI - dPhi; - - return dPhi/std::abs(omegaCalo)*std::sqrt(1. + tanLCalo*tanLCalo); -} - - -double TOFEstimators::getTOFClosest(const Track* track, const Cluster* cluster){ - const TrackState* ts = track->getTrackState(TrackState::AtCalorimeter); - dd4hep::rec::Vector3D trackAtCaloPos = ts->getReferencePoint(); - - double dToImpactMin = std::numeric_limits::max(); - double hitTime = 0.; - for (const auto& hit : cluster->getCalorimeterHits() ){ - CHT hitType( hit->getType() ); - bool isEcal = (hitType.caloID() == CHT::ecal); - if (!isEcal) continue; - - dd4hep::rec::Vector3D hitPos = hit->getPosition() ; - double dToImpact = (hitPos - trackAtCaloPos).r(); - - if(dToImpact < dToImpactMin){ - dToImpactMin = dToImpact; - hitTime = hit->getTime(); + Track* track = pfo->getTracks()[0]; + Cluster* cluster = pfo->getClusters()[0]; + + /////////////////////////////////////////////////////////////// + // This part calculates track length and momentum harmonic mean + /////////////////////////////////////////////////////////////// + vector subTracks = getSubTracks(track); + vector trackStates = getTrackStatesPerHit(subTracks, _trkSystem, _extrapolateToEcal, _bField); + + double trackLength = 0.; + double harmonicMom = 0.; + int nTrackStates = trackStates.size(); + for( int j=1; j < nTrackStates; ++j ){ + if (j == nTrackStates - 1){ + double nTurns = getHelixNRevolutions( trackStates[j-1], trackStates[j] ); + if( nTurns > 0.5 ){ + // we cannot calculate arc length for more than pi revolution. Use formula with only z + double arcLength = getHelixLengthAlongZ( trackStates[j-1], trackStates[j] ); + Vector3D mom = getHelixMomAtTrackState( trackStates[j-1], _bField ); + trackLength += arcLength; + harmonicMom += arcLength/mom.r2(); + continue; + } + } + + double arcLength = getHelixArcLength( trackStates[j-1], trackStates[j] ); + Vector3D mom = getHelixMomAtTrackState( trackStates[j-1], _bField ); + trackLength += arcLength; + harmonicMom += arcLength/mom.r2(); } - } - hitTime += _smearing(_generator); - return hitTime - dToImpactMin / CLHEP::c_light; -} - - -double TOFEstimators::getTOFFastest(const Track* track, const Cluster* cluster){ - const TrackState* ts = track->getTrackState(TrackState::AtCalorimeter); - dd4hep::rec::Vector3D trackAtCaloPos = ts->getReferencePoint(); - - double dToImpactMin = 0.; - double hitTimeMin = std::numeric_limits::max(); - for (const auto& hit : cluster->getCalorimeterHits() ){ - CHT hitType( hit->getType() ); - bool isEcal = (hitType.caloID() == CHT::ecal); - if (!isEcal) continue; - - dd4hep::rec::Vector3D hitPos = hit->getPosition(); - double dToImpact = (hitPos - trackAtCaloPos).r(); - double hitTime = hit->getTime() + _smearing(_generator); - if(hitTime < hitTimeMin){ - dToImpactMin = dToImpact; - hitTimeMin = hitTime; + harmonicMom = std::sqrt(trackLength/harmonicMom); + + + ////////////////////////////////////// + // This part calculates Time of flight + ////////////////////////////////////// + double timeOfFlight = 0.; + if( _extrapolateToEcal ){ + // if we extrapolate to the calorimeter check TOF method steering parameter + if (_tofMethod == "closest"){ + timeOfFlight = getTofClosest(cluster, track, _timeResolution/1000.); + } + else if (_tofMethod == "frankAvg"){ + vector frankHits = selectFrankEcalHits(cluster, track, _maxEcalLayer, _bField); + timeOfFlight = getTofFrankAvg(frankHits, track, _timeResolution/1000.); + } + else if (_tofMethod == "frankFit"){ + vector frankHits = selectFrankEcalHits(cluster, track, _maxEcalLayer, _bField); + timeOfFlight = getTofFrankFit(frankHits, track, _timeResolution/1000.); + } + else{ + //wrong parameter silly! + } } - } - return hitTimeMin - dToImpactMin / CLHEP::c_light; -} - - -double TOFEstimators::getTOFCylFit(const Track* track, const Cluster* cluster){ - const TrackState* ts = track->getTrackState(TrackState::AtCalorimeter); - dd4hep::rec::Vector3D trackAtCaloPos = ts->getReferencePoint(); - dd4hep::rec::Vector3D trackAtCaloMom = getMomAtCalo(track); - - std::vector x, xErr, y, yErr; - - for (const auto& hit : cluster->getCalorimeterHits() ){ - CHT hitType( hit->getType() ); - bool isEcal = (hitType.caloID() == CHT::ecal); - int layer = hitType.layer(); - if (!isEcal || layer >= _maxLayerNum) continue; - - dd4hep::rec::Vector3D hitPos = hit->getPosition(); - double dToLine = (hitPos - trackAtCaloPos).cross(trackAtCaloMom.unit()).r(); - // take only hits from 5 mm cylinders. 5 mm value was ibtain optimizing on photons - if (dToLine > _cylRadiusCut) continue; - - x.push_back( (hitPos - trackAtCaloPos).r() ); - y.push_back( hit->getTime() + _smearing(_generator) ); - xErr.push_back(0.); - yErr.push_back( 0.1 ); - } - - if (x.size() <= 1) return 0.; - - TGraphErrors gr(x.size(), &x[0], &y[0], &xErr[0], &yErr[0]); - gr.Fit("pol1", "Q"); - - return gr.GetFunction("pol1")->GetParameter(0); -} - - -double TOFEstimators::getTOFClosestFit(const Track* track, const Cluster* cluster){ - const TrackState* ts = track->getTrackState(TrackState::AtCalorimeter); - dd4hep::rec::Vector3D trackAtCaloPos = ts->getReferencePoint(); - dd4hep::rec::Vector3D trackAtCaloMom = getMomAtCalo(track); - - std::map dToImpact; - std::map hitTime; - std::map dToLineMin; - for(int i=0; i<_maxLayerNum; ++i) dToImpact[i] = 0.; - for(int i=0; i<_maxLayerNum; ++i) hitTime[i] = 0.; - for(int i=0; i<_maxLayerNum; ++i) dToLineMin[i] = std::numeric_limits::max(); - - for (const auto& hit : cluster->getCalorimeterHits() ){ - CHT hitType( hit->getType() ); - bool isEcal = (hitType.caloID() == CHT::ecal); - int layer = hitType.layer(); - if (!isEcal || layer >= _maxLayerNum) continue; - - dd4hep::rec::Vector3D hitPos = hit->getPosition(); - double dToLine = (hitPos - trackAtCaloPos).cross(trackAtCaloMom.unit()).r(); - if( dToLine < dToLineMin[layer] ){ - dToLineMin[layer] = dToLine; - dToImpact[layer] = (hitPos - trackAtCaloPos).r(); - hitTime[layer] = hit->getTime(); + else{ + //define tof as average time between two SET strips + //if no SET hits found, tof=0 + TrackerHit* hitSET = getSETHit(track, _tpcOuterR); + if ( hitSET != nullptr ){ + const vector& simHitsSET = navigatorSET.getRelatedToObjects( hitSET ); + if ( simHitsSET.size() == 2 ){ + //it should be 2 always.. but just to be safe + SimTrackerHit* simHitSETFront = static_cast ( simHitsSET[0] ); + SimTrackerHit* simHitSETBack = static_cast ( simHitsSET[1] ); + double timeFront = RandGauss::shoot( simHitSETFront->getTime(), _timeResolution/1000. ); + double timeBack = RandGauss::shoot( simHitSETBack->getTime(), _timeResolution/1000. ); + timeOfFlight = (timeFront + timeBack)/2.; + } + } } - } - - std::vector x, xErr, y, yErr; - for(int i=0; i<_maxLayerNum; ++i){ - if (hitTime[i] <= 0.) continue; - x.push_back(dToImpact[i]); - y.push_back( hitTime[i] + _smearing(_generator) ); - xErr.push_back(0.); - yErr.push_back(0.1); + vector results{float(harmonicMom), float(trackLength), float(timeOfFlight)}; + pidHandler.setParticleID(pfo , 0, 0, 0., algoID, results); } - - if (x.size() <= 1) return 0.; - - TGraphErrors gr(x.size(), &x[0], &y[0], &xErr[0], &yErr[0]); - gr.Fit("pol1", "Q"); - - return gr.GetFunction("pol1")->GetParameter(0); } diff --git a/TimeOfFlight/src/TOFUtils.cc b/TimeOfFlight/src/TOFUtils.cc index 0ef9d235..aae6c66d 100644 --- a/TimeOfFlight/src/TOFUtils.cc +++ b/TimeOfFlight/src/TOFUtils.cc @@ -2,150 +2,321 @@ #include #include - -#include -#include -#include -#include -#include -#include +#include #include "marlinutil/CalorimeterHitType.h" +#include "HelixClass.h" +#include "DD4hep/Detector.h" +#include "DD4hep/DD4hepUnits.h" +#include "DDRec/DetectorData.h" +#include "CLHEP/Random/Randomize.h" +#include "CLHEP/Units/PhysicalConstants.h" +#include "TGraphErrors.h" +#include "TF1.h" +#include "MarlinTrk/MarlinTrkUtils.h" +#include "UTIL/LCRelationNavigator.h" +#include "UTIL/ILDConf.h" +#include "IMPL/TrackImpl.h" + +using std::vector; +using std::numeric_limits; +using std::pair; +using EVENT::TrackerHit; +using EVENT::Track; +using EVENT::Cluster; +using EVENT::CalorimeterHit; +using EVENT::TrackState; +using dd4hep::Detector; +using dd4hep::DetElement; +using dd4hep::rec::Vector3D; +using dd4hep::rec::FixedPadSizeTPCData; +using MarlinTrk::IMarlinTrack; +using MarlinTrk::IMarlinTrkSystem; +using UTIL::LCRelationNavigator; +using UTIL::ILDDetID; +using IMPL::TrackImpl; +using IMPL::TrackStateImpl; +using CLHEP::RandGauss; + +bool TOFUtils::sortByRho(TrackerHit* a, TrackerHit* b){ + Vector3D posA( a->getPosition() ); + Vector3D posB( b->getPosition() ); + return posA.rho() < posB.rho(); +} -namespace TOFUtils{ - -using namespace lcio ; - - -float computeFlightLength( EVENT::Track* trk){ +TrackStateImpl TOFUtils::getTrackStateAtHit(IMarlinTrack* marlinTrk, TrackerHit* hit){ + TrackStateImpl ts; + double chi2Dummy; + int ndfDummy; + marlinTrk->getTrackState(hit, ts, chi2Dummy, ndfDummy); + return ts; +} - const TrackState* tsIp = trk->getTrackState( TrackState::AtIP ) ; - const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ; - float phiIp = tsIp->getPhi() ; - float phiCalo = tscalo->getPhi() ; +Vector3D TOFUtils::getHelixMomAtTrackState(const TrackState& ts, double bField){ + double phi = ts.getPhi(); + double d0 = ts.getD0(); + double z0 = ts.getZ0(); + double omega = ts.getOmega(); + double tanL = ts.getTanLambda(); - float omega = tsIp->getOmega() ; - float tanL = tsIp->getTanLambda() ; + HelixClass helix; + helix.Initialize_Canonical(phi, d0, z0, omega, tanL, bField); + return helix.getMomentum(); +} - float dPhi = std::abs(phiIp - phiCalo); - // if dPhi > PI assume phi is flipped +double TOFUtils::getHelixArcLength(const TrackState& ts1, const TrackState& ts2){ + double omega = ts1.getOmega(); + double z1 = ts1.getReferencePoint()[2] + ts1.getZ0(); + double z2 = ts2.getReferencePoint()[2] + ts2.getZ0(); + double dPhi = std::abs( ts2.getPhi() - ts1.getPhi() ); if (dPhi > M_PI) dPhi = 2*M_PI - dPhi; - return dPhi/std::abs(omega) * std::sqrt( 1. + tanL * tanL ) ; + return std::sqrt( std::pow(dPhi/omega, 2) + std::pow(z2-z1, 2) ); } -float computeFlightLength(const TrackState* ts0, const TrackState* ts1 ){ - float phi1 = ts1->getPhi() ; - float phi0 = ts0->getPhi() ; - - float omega = ts0->getOmega() ; - float tanL = ts0->getTanLambda() ; - float dPhi = std::abs(phi0 - phi1); +double TOFUtils::getHelixLengthAlongZ(const TrackState& ts1, const TrackState& ts2){ + double tanL = ts1.getTanLambda(); + double z1 = ts1.getReferencePoint()[2] + ts1.getZ0(); + double z2 = ts2.getReferencePoint()[2] + ts2.getZ0(); - // if dPhi > PI assume phi is flipped - if (dPhi > M_PI) dPhi = 2*M_PI - dPhi; - - return dPhi/std::abs(omega) * std::sqrt( 1. + tanL * tanL ) ; + return std::abs( (z2-z1)/tanL ) * std::sqrt( 1.+tanL*tanL ); } +double TOFUtils::getHelixNRevolutions(const TrackState& ts1, const TrackState& ts2){ + double omega = ts1.getOmega(); + double tanL = ts1.getTanLambda(); + double z1 = ts1.getReferencePoint()[2] + ts1.getZ0(); + double z2 = ts2.getReferencePoint()[2] + ts2.getZ0(); - int layer( EVENT::CalorimeterHit* h ) { - - return CHT( h->getType() ).layer() ; - } + // helix length projected on xy + double circHelix = std::abs( (z2-z1)/tanL ); + double circFull = 2*M_PI/std::abs(omega); + return circHelix/circFull; +} - /// helper function to check if this is an Ecal hit - bool isEcal( EVENT::CalorimeterHit* h ) { - return CHT( h->getType() ).caloID() == CHT::ecal ; - } +double TOFUtils::getTPCOuterR(){ + auto& detector = Detector::getInstance(); + DetElement tpcDet = detector.detector("TPC"); + FixedPadSizeTPCData* tpc = tpcDet.extension (); + return tpc->rMaxReadout/dd4hep::mm; +} - std::string caloTypeStr( EVENT::CalorimeterHit* h ) { +TrackerHit* TOFUtils::getSETHit(Track* track, double tpcOuterR){ + vector hits = track->getTrackerHits(); + TrackerHit* lastHit = hits.back(); - std::stringstream s ; - s << CHT( h->getType() ) ; - return s.str() ; - } + Vector3D pos (lastHit->getPosition()); + if ( pos.rho() > tpcOuterR ) return lastHit; + return nullptr; +} - std::string CaloHitData::toString(){ +vector TOFUtils::selectFrankEcalHits( Cluster* cluster, Track* track, int maxEcalLayer, double bField ){ + vector selectedHits; + + const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); + Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); + Vector3D trackMomAtEcal = TOFUtils::getHelixMomAtTrackState(*tsEcal, bField); + + for (int l=0; l::max(); + for ( auto hit : cluster->getCalorimeterHits() ){ + CHT hitType( hit->getType() ); + bool isECALHit = ( hitType.caloID() == CHT::ecal ); + if ( (!isECALHit) || ( int( hitType.layer() ) != l) ) continue; + + Vector3D hitPos( hit->getPosition() ); + double dToLine = (hitPos - trackPosAtEcal).cross(trackMomAtEcal.unit()).r(); + if (dToLine < closestDistanceToLine){ + closestDistanceToLine = dToLine; + selectedHit = hit; + } + } + if ( selectedHit != nullptr ) selectedHits.push_back(selectedHit); + } + return selectedHits; +} - std::stringstream s ; - s << " l= " << layer ; - s << ", st= " << smearedTime ; - s << ", tr= " << timeResolution ; - s << ", dIP=" << distanceFromIP ; - s << ", dRP= " << distanceFromReferencePoint ; - s << ", dSL= " << distancefromStraightline ; -// EVENT::CalorimeterHit* lcioHit = nullptr ; - return s.str() ; - } +vector TOFUtils::getSubTracks(Track* track){ + vector subTracks; + subTracks.push_back(track); + int nSubTracks = track->getTracks().size(); + if (nSubTracks <= 1) return subTracks; - float computeDistanceFromLine( EVENT::CalorimeterHit* h, const dd4hep::rec::Vector3D& point, - const dd4hep::rec::Vector3D& unitDir) { + int nTPCHits = track->getSubdetectorHitNumbers()[(ILDDetID::TPC)*2-1]; + int nSubTrack0Hits = track->getTracks()[0]->getTrackerHits().size(); + int nSubTrack1Hits = track->getTracks()[1]->getTrackerHits().size(); - dd4hep::rec::Vector3D pos( h->getPosition()[0], - h->getPosition()[1], - h->getPosition()[2] ) ; + int startIdx; + if( std::abs(nTPCHits - nSubTrack0Hits) <= 1 ) startIdx = 1; + else if ( std::abs(nTPCHits - nSubTrack1Hits) <= 1 ) startIdx = 2; + else{ + // this shouldn't happen in princinple at all... + return subTracks; + } + for(int j=startIdx; j < nSubTracks; ++j) subTracks.push_back( track->getTracks()[j] ); + return subTracks; +} - dd4hep::rec::Vector3D diff = pos - point ; +vector TOFUtils::getTrackStatesPerHit(vector tracks, IMarlinTrkSystem* trkSystem, bool extrapolateToEcal, double bField){ + vector trackStatesPerHit; + int nTracks = tracks.size(); + for(int i=0; i hits = track->getTrackerHits(); + sort(hits.begin(), hits.end(), sortByRho); + + // setup initial dummy covariance matrix + vector covMatrix(15); + // initialize variances + covMatrix[0] = 1e+06; //sigma_d0^2 + covMatrix[2] = 100.; //sigma_phi0^2 + covMatrix[5] = 0.00001; //sigma_omega^2 + covMatrix[9] = 1e+06; //sigma_z0^2 + covMatrix[14] = 100.; //sigma_tanl^2 + double maxChi2PerHit = 100.; + IMarlinTrack* marlinTrk = trkSystem->createTrack(); + TrackImpl refittedTrack; + + //Need to initialize trackState at last hit + TrackStateImpl preFit = *track->getTrackState(TrackState::AtLastHit); + preFit.setCovMatrix( covMatrix ); + int errorFit = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, hits, &refittedTrack, IMarlinTrack::backward, &preFit, bField, maxChi2PerHit); + if (errorFit != 0) continue; + + vector< pair > hitsInFit; + marlinTrk->getHitsInFit(hitsInFit); + + int nHitsInFit = hitsInFit.size(); + // if first subTrack + if (i == 0){ + trackStatesPerHit.push_back(*(static_cast (refittedTrack.getTrackState(TrackState::AtIP)) )); + + //add hits in increasing rho for the FIRST subTrack!!!!! + for( int j=nHitsInFit-1; j>=0; --j ){ + TrackStateImpl ts = getTrackStateAtHit(marlinTrk, hitsInFit[j].first); + trackStatesPerHit.push_back(ts); + } + } + else{ + // check which hit is closer to the last hit of previous fit. + // and iterate starting from the closest + Vector3D innerHit ( hitsInFit.back().first->getPosition() ); + Vector3D outerHit ( hitsInFit.front().first->getPosition() ); + Vector3D prevHit ( trackStatesPerHit.back().getReferencePoint() ); + + if ( (innerHit - prevHit).r() < (outerHit - prevHit).r() ){ + for( int j=nHitsInFit-1; j>=0; --j ){ + //iterate in increasing rho + TrackStateImpl ts = getTrackStateAtHit(marlinTrk, hitsInFit[j].first); + trackStatesPerHit.push_back(ts); + } + } + else{ + for( int j=0; j (refittedTrack.getTrackState(TrackState::AtLastHit)) ) ); + if (extrapolateToEcal) trackStatesPerHit.push_back( *(static_cast (refittedTrack.getTrackState(TrackState::AtCalorimeter) ) ) ); + } + } + // one can maybe use hits of refittedTrack, which includes hits that failed in the fit + // code would look cleaner, but using hits failed in fit probably would have worse performance.. + // needs to be checked + return trackStatesPerHit; +} - return diff.cross( unitDir ).r() ; - } +double TOFUtils::getTofClosest( Cluster* cluster, Track* track, double timeResolution){ + const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); + Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); + double hitTime = numeric_limits::max(); + double closestDistance = numeric_limits::max(); + for( auto hit : cluster->getCalorimeterHits() ){ + CHT hitType( hit->getType() ); + bool isECALHit = ( hitType.caloID() == CHT::ecal ); + if (! isECALHit) continue; - CaloHitDataVec findHitsClosestToLine( const CaloHitLayerMap& layerMap ){ + Vector3D hitPos( hit->getPosition() ); + double dToTrack = (hitPos - trackPosAtEcal).r(); + if( dToTrack < closestDistance ){ + closestDistance = dToTrack; + hitTime = hit->getTime(); + } + } - CaloHitDataVec hitVec ; + return RandGauss::shoot(hitTime, timeResolution) - closestDistance/CLHEP::c_light; +} - for( auto m : layerMap ){ - const CaloHitDataVec& chv = m.second ; +double TOFUtils::getTofFrankAvg( vector selectedHits, Track* track, double timeResolution){ + const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); + Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); - CaloHitData* closestHit = - *std::min_element( chv.begin() , chv.end () , - [](CaloHitData* c0, CaloHitData* c1 ){ return c0->distancefromStraightline < c1->distancefromStraightline ; } - ) ; + int nHits = selectedHits.size(); + if (nHits == 0) return 0.; - hitVec.push_back( closestHit ) ; + double tof = 0.; + for ( auto hit : selectedHits ){ + Vector3D hitPos( hit->getPosition() ); + double dToTrack = (hitPos - trackPosAtEcal).r(); + tof += RandGauss::shoot(hit->getTime(), timeResolution) - dToTrack/CLHEP::c_light; } + return tof/nHits; +} - return hitVec ; - } - - - std::pair computeTOFEstimator( const CaloHitDataVec& chv ){ - - const static float c_mm_per_ns = 299.792458 ; - int nHit = 0 ; - float mean = 0. ; - float meansq = 0. ; +double TOFUtils::getTofFrankFit( vector selectedHits, Track* track, double timeResolution){ + const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); + Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); - for( auto ch : chv ){ - float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; - mean += t ; - ++nHit ; + int nHits = selectedHits.size(); + if (nHits == 0) return 0.; + else if (nHits == 1){ + //we can't fit 1 point, but lets return something reasonable + Vector3D hitPos( selectedHits[0]->getPosition() ); + double dToTrack = (hitPos - trackPosAtEcal).r(); + return RandGauss::shoot(selectedHits[0]->getTime(), timeResolution) - dToTrack/CLHEP::c_light; } - mean /= nHit ; - for( auto ch : chv ){ - float t = ch->smearedTime - ch->distanceFromReferencePoint / c_mm_per_ns ; - meansq += ( t - mean ) * ( t - mean ) ; + vector x, y, xErr, yErr; + for ( auto hit : selectedHits ){ + Vector3D hitPos( hit->getPosition() ); + double dToTrack = (hitPos - trackPosAtEcal).r(); + x.push_back(dToTrack); + double time = RandGauss::shoot(hit->getTime(), timeResolution); + y.push_back(time); + xErr.push_back(0.); + yErr.push_back(0.3); + //setting this to 0 is not good for the fit.. So I put random 300ps... + //Changing this doesn't seem to affect results, although one may want to check it more carefully } - return std::make_pair( mean, std::sqrt(meansq/nHit ) ); - } - - -} //namespace + TGraphErrors gr(nHits, &x[0], &y[0], &xErr[0], &yErr[0]); + gr.Fit("pol1", "Q"); + return gr.GetFunction("pol1")->GetParameter(0); +} diff --git a/TimeOfFlight/xml/steer.xml b/TimeOfFlight/xml/steer.xml new file mode 100644 index 00000000..4b1dfbc8 --- /dev/null +++ b/TimeOfFlight/xml/steer.xml @@ -0,0 +1,43 @@ + + + + + + + + + + + + + /pnfs/desy.de/ilc/prod/ilc/mc-2020/ild/rec/250-SetA/2f_hadronic_eL_pR/ILD_l5_o1_v02/v02-02/00015161/000/rv02-02.sv02-02.mILD_l5_o1_v02.E250-SetA.I500010.P2f_z_h.eL.pR.n001_209.d_rec_00015161_259.slcio + + + + + + MESSAGE DEBUG0 + + + + + + /cvmfs/ilc.desy.de/sw/x86_64_gcc82_centos7/v02-02-02/lcgeo/v00-16-06/ILD/compact/ILD_l5_o1_v02/ILD_l5_o1_v02.xml + + + + + true + closest + 0 + + + + + test.slcio + WRITE_NEW + * + + + + From 6c25184a0c462b6fe7dfb1c7ee61e36539c1180f Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Tue, 21 Sep 2021 19:35:06 +0200 Subject: [PATCH 03/22] Add namespaces in definitions for Doxygen --- TimeOfFlight/src/TOFEstimators.cc | 2 +- TimeOfFlight/src/TOFUtils.cc | 26 +++++++++++++------------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index 9f34c799..6c7c00e8 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -85,7 +85,7 @@ void TOFEstimators::init(){ } -void TOFEstimators::processEvent(LCEvent * evt){ +void TOFEstimators::processEvent(EVENT::LCEvent * evt){ RandGauss::setTheSeed( marlin::Global::EVENTSEEDER->getSeed(this) ); ++_nEvent; streamlog_out(MESSAGE)<<"******Event****** "<<_nEvent<getPosition() ); Vector3D posB( b->getPosition() ); return posA.rho() < posB.rho(); } -TrackStateImpl TOFUtils::getTrackStateAtHit(IMarlinTrack* marlinTrk, TrackerHit* hit){ +IMPL::TrackStateImpl TOFUtils::getTrackStateAtHit(MarlinTrk::IMarlinTrack* marlinTrk, EVENT::TrackerHit* hit){ TrackStateImpl ts; double chi2Dummy; int ndfDummy; @@ -54,7 +54,7 @@ TrackStateImpl TOFUtils::getTrackStateAtHit(IMarlinTrack* marlinTrk, TrackerHit* } -Vector3D TOFUtils::getHelixMomAtTrackState(const TrackState& ts, double bField){ +dd4hep::rec::Vector3D TOFUtils::getHelixMomAtTrackState(const EVENT::TrackState& ts, double bField){ double phi = ts.getPhi(); double d0 = ts.getD0(); double z0 = ts.getZ0(); @@ -67,7 +67,7 @@ Vector3D TOFUtils::getHelixMomAtTrackState(const TrackState& ts, double bField){ } -double TOFUtils::getHelixArcLength(const TrackState& ts1, const TrackState& ts2){ +double TOFUtils::getHelixArcLength(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2){ double omega = ts1.getOmega(); double z1 = ts1.getReferencePoint()[2] + ts1.getZ0(); double z2 = ts2.getReferencePoint()[2] + ts2.getZ0(); @@ -78,7 +78,7 @@ double TOFUtils::getHelixArcLength(const TrackState& ts1, const TrackState& ts2) } -double TOFUtils::getHelixLengthAlongZ(const TrackState& ts1, const TrackState& ts2){ +double TOFUtils::getHelixLengthAlongZ(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2){ double tanL = ts1.getTanLambda(); double z1 = ts1.getReferencePoint()[2] + ts1.getZ0(); double z2 = ts2.getReferencePoint()[2] + ts2.getZ0(); @@ -86,7 +86,7 @@ double TOFUtils::getHelixLengthAlongZ(const TrackState& ts1, const TrackState& t return std::abs( (z2-z1)/tanL ) * std::sqrt( 1.+tanL*tanL ); } -double TOFUtils::getHelixNRevolutions(const TrackState& ts1, const TrackState& ts2){ +double TOFUtils::getHelixNRevolutions(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2){ double omega = ts1.getOmega(); double tanL = ts1.getTanLambda(); double z1 = ts1.getReferencePoint()[2] + ts1.getZ0(); @@ -108,7 +108,7 @@ double TOFUtils::getTPCOuterR(){ } -TrackerHit* TOFUtils::getSETHit(Track* track, double tpcOuterR){ +EVENT::TrackerHit* TOFUtils::getSETHit(EVENT::Track* track, double tpcOuterR){ vector hits = track->getTrackerHits(); TrackerHit* lastHit = hits.back(); @@ -118,7 +118,7 @@ TrackerHit* TOFUtils::getSETHit(Track* track, double tpcOuterR){ } -vector TOFUtils::selectFrankEcalHits( Cluster* cluster, Track* track, int maxEcalLayer, double bField ){ +std::vector TOFUtils::selectFrankEcalHits( EVENT::Cluster* cluster, EVENT::Track* track, int maxEcalLayer, double bField ){ vector selectedHits; const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); @@ -148,7 +148,7 @@ vector TOFUtils::selectFrankEcalHits( Cluster* cluster, Track* } -vector TOFUtils::getSubTracks(Track* track){ +std::vector TOFUtils::getSubTracks(EVENT::Track* track){ vector subTracks; subTracks.push_back(track); @@ -171,7 +171,7 @@ vector TOFUtils::getSubTracks(Track* track){ } -vector TOFUtils::getTrackStatesPerHit(vector tracks, IMarlinTrkSystem* trkSystem, bool extrapolateToEcal, double bField){ +std::vector TOFUtils::getTrackStatesPerHit(std::vector tracks, MarlinTrk::IMarlinTrkSystem* trkSystem, bool extrapolateToEcal, double bField){ vector trackStatesPerHit; int nTracks = tracks.size(); for(int i=0; i TOFUtils::getTrackStatesPerHit(vector tracks, IMa } -double TOFUtils::getTofClosest( Cluster* cluster, Track* track, double timeResolution){ +double TOFUtils::getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, double timeResolution){ const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); @@ -273,7 +273,7 @@ double TOFUtils::getTofClosest( Cluster* cluster, Track* track, double timeResol } -double TOFUtils::getTofFrankAvg( vector selectedHits, Track* track, double timeResolution){ +double TOFUtils::getTofFrankAvg( std::vector selectedHits, EVENT::Track* track, double timeResolution){ const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); @@ -290,7 +290,7 @@ double TOFUtils::getTofFrankAvg( vector selectedHits, Track* tr } -double TOFUtils::getTofFrankFit( vector selectedHits, Track* track, double timeResolution){ +double TOFUtils::getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution){ const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); From 0a7ea4bffcd6e80e014f2b46af7c3dec4e697138 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Thu, 23 Sep 2021 01:34:23 +0200 Subject: [PATCH 04/22] Updated documentation,README and steering file example --- TimeOfFlight/README.md | 190 +++++++++++++++++++++++++-- TimeOfFlight/include/TOFEstimators.h | 149 ++++++++++----------- TimeOfFlight/include/TOFUtils.h | 90 ++++++++++--- TimeOfFlight/xml/steer.xml | 76 +++++++++-- 4 files changed, 378 insertions(+), 127 deletions(-) diff --git a/TimeOfFlight/README.md b/TimeOfFlight/README.md index 2b113a0c..90cb5281 100644 --- a/TimeOfFlight/README.md +++ b/TimeOfFlight/README.md @@ -1,21 +1,187 @@ +# Description -# TimeOfFlight +TOFEstimators is the Marlin processor that calculates track length and time-of-flight for charged particles that can be used for the mass reconstruction and particle identification.
+In addition it calculates square root of the harmonic mean of the momentum squared. One can use it as a better substitute for the particle momentum in the mass reconstruction. -This package deals with computing and displaying time of flight parameters -based on the timing information stored in CalorimeterHits. +[Doxygen documentation](https://www.desy.de/~dudarboh/tof_doc/html/index.html) contains a detailed description of the source code.
+In case reader uses dark theme on the github and finds hard to read the formulas he can find this README as the main page there. +## Requirements -## Processors +TOFEstimators is dependent on InitDD4hep processor to extract detector geometry details.
+Make sure to run InitDD4hep before this processor is executed. -### TOFEstimators +In addition, input sclio file must contain a lot of various tracker and calorimeter hit collections.
+Thus, input slcio file can be e.g. *REC* format, but **cannot** be *DST* or *mini-DST* format. -author: N.Weinhold (DESY, internship 2017) -author: F.Gaede, DESY, 2018 -author: B.Dudar, DESY, 2021 -Compute various estimators for time of flight and add these to the -ReconstrucedParticles as PIDParameters. +## Steering parameters +1. **ReconstructedParticleCollection** / default: `"PandoraPFOs"`
+ Collection of `ReconstructedParticle` objects to analyze.
+ The final results will be written in the PIDHandler inside this collection. -## Example steering file +2. **ExtrapolateToEcal** / default: `true`
+ If `true` the track length is calculated to the extrapolated track position at the ECal surface. Time-of-flight then calculated using ECal hits via one of the methods chosen in *TofMethod* steering parameter (see below). -[./xml/steer.xml](./xml/steer.xml) is an example of a steering file that runs TOFEstimators processor. It computes momentum harmonic mean, time-of-flight and track length of the PandoraPFOs and stores this information inside PIDHandler. In the example time of flight calculated using closest Ecal hit assuming perfect time resolution. Then new slcio file "test.slcio" produced with LCIOOutputProcessor processor. + If `false` the track length is calculated to the last tracker hit. Time-of-flight then calculated using the last tracker hit. Which is SET hit for the barrel region and *undefined* for the endcap region.
+ If the last tracker hit is the SET hit time-of-flight is set to the average time of two SET strips.
+ If the last tracker hit is the TPC hit time-of-flight is set to `0`. + +3. **TofMethod** / options: `"closest"`, `"frankAvg"`, `"frankFit"` / default: `"closest"`
+ If *ExtrapolateToEcal* is set to `true` then it defines how to calculate the time-of-flight at the ECal surface. + + - `"closest"` uses the closest ECal hit to the extrapolated track position at the ECal surface.
+ Time-of-flight is defined as: . + + - `"frankAvg"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Define time-of-flight as an average of their time corrected for the distance to the track position at the ECal surface assuming speed of flight is the speed of light. + . + + - `"frankFit"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Use linear fit of the hits time as a function of the distance to the track position at the ECal surface to define time-of-flight as an extrapolated time at the ECal surface . + + for more illustrative explanations of the methods above see slides 10, 12, 13 from the following [LCWS2021 talk]((https://indico.cern.ch/event/995633/contributions/4259659/attachments/2209010/3738157/Bohdan_TOF_LCWS2021.pdf)). + + If *ExtrapolateToEcal* is set to `false` then this steering parameter is ignored. + + +4. **TimeResolution** / default: `0.0`
+ If *ExtrapolateToEcal* is set to `true` it defines times resolution of individual ECal hits in ps. + + If *ExtrapolateToEcal* is set to `false` it defines times resolution of individual SET strips in ps. + +5. **MaxEcalLayer** / default: `10`
+ Defines number of inner ECal layers to use for the *frankAvg* and *frankFit* algorithms. + + If *frankAvg* and *frankFit* are unused then this steering parameter is ignored. + + +## Output parameters + +1. `momentumHM` + + A squared root of the harmonic mean of the squared momentum: . + + If momentum is constant, then it is just a momentum: .
+ If momentum changes, then using instead of or is mathematically rigorous choice for the relativistic particle assumption .
+ See [Winfried A. Mitaroff paper](https://arxiv.org/abs/2107.02031) for the details. + +2. `trackLength` + + A track length calculated between the point of the closest approach (PCA) to the IP `(0,0,0)`
+ and the ECal surface if *ExtrapolateToEcal* is `true`
+ or the last tracker hit if *ExtrapolateToEcal* is `false`. + + A total track length is calculated as a sum of track segments obtained from iterating over track states defined at the IP, every tracker hit and the ECal surface if *ExtrapolateToEcal* is `true`.
+ A track segment length is calculated between neighboring track states as: + . + + If *ExtrapolateToEcal* is `true` then we additionally check for the number of helix turns between the last tracker hit and the extrapolated track position at the ECal surface. If we are unable to use the formula above to calculate the length of the last segment. Thus we use different formula that does not rely on the azimuthal angle: + + . + + +3. `timeOfFlight` + + Time-of-flight of the particle assuming . For the detailed description of different time-of-flight estimators see the **TofMethod** bullet point in the Steering parameters section. + + + + +## Steering file example + +[./xml/steer.xml](./xml/steer.xml) is a steering file example that runs three TOFEstimators processors: *MyTofClosest0ps*, *MyTofSET10ps* and *MyTofFrankAvg50ps* and then writes their output parameters in the PIDHandlers of the *PandoraPFOs* collection in the new *output.slcio* file using LCIOOutputProcessor. + +To run this example one needs to setup iLCSoft environment. +If the reader has a NAF account he/she can setup iLCSoft environment with:
+`source /cvmfs/ilc.desy.de/sw/x86_64_gcc82_centos7/v02-02-02/init_ilcsoft.sh` + +and run the example steering file with:
+`Marlin ./xml/steer.xml` + +Then one can look at the output.sclio file with, e.g.:
+`dumpevent output.slcio 1 | less` + +Scrolling carefully one can find these lines which indicate that the new file contains our parameters for all three TOFEstimators: + +``` +collection name : PandoraPFOs +parameters: + +--------------- print out of ReconstructedParticle collection --------------- + + flag: 0x0 +. . . +parameter ParameterNames_MyTofClosest0ps [string]: momentumHM, trackLength, timeOfFlight, +parameter ParameterNames_MyTofFrankAvg50ps [string]: momentumHM, trackLength, timeOfFlight, +parameter ParameterNames_MyTofSET10ps [string]: momentumHM, trackLength, timeOfFlight, +. . . +``` + +Scrolling more down, one can see calculated parameters for each individual particle:
+``` +------------ detailed PID info: --- + + algorithms : + [id: 0] BasicVariablePID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood M + [id: 3] LikelihoodPID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MVAO + [id: 4] LowMomMuID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MVAOutp + [id: 9] MyTofClosest0ps - params: momentumHM trackLength timeOfFlight + [id: 11] MyTofFrankAvg50ps - params: momentumHM trackLength timeOfFlight + [id: 10] MyTofSET10ps - params: momentumHM trackLength timeOfFlight + [id: 2] ShowerShapesPID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MV + [id: 5] TOFEstimators0ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightLen + [id: 8] TOFEstimators100ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightL + [id: 6] TOFEstimators10ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightLe + [id: 7] TOFEstimators50ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightLe + [id: 1] dEdxPID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MVAOutput_ + + + [particle] | PDG | likelihood | type | algoId | parameters : + | | | | | +. . . + [00000073] | 2212 | 2.0802e+00 | 000000 | 1 | [ electronLikelihood : -1.64e+02, muonLikelihood : -3.47e+01, pionLikelihood : -2. + | 0 | 0.0000e+00 | 000000 | 4 | [ electronLikelihood : 9.99e+02, muonLikelihood : 9.99e+02, pionLikelihood : 9.99e + | 0 | 0.0000e+00 | 000000 | 5 | [ TOFFirstHit : 9.58e+00, TOFClosestHits : 9.58e+00, TOFClosestHitsError : 2.71e-0 + | 0 | 0.0000e+00 | 000000 | 6 | [ TOFFirstHit : 9.57e+00, TOFClosestHits : 9.58e+00, TOFClosestHitsError : 7.16e-0 + | 0 | 0.0000e+00 | 000000 | 7 | [ TOFFirstHit : 9.54e+00, TOFClosestHits : 9.58e+00, TOFClosestHitsError : 6.32e-0 + | 0 | 0.0000e+00 | 000000 | 8 | [ TOFFirstHit : 9.66e+00, TOFClosestHits : 9.64e+00, TOFClosestHitsError : 8.34e-0 + | 0 | 0.0000e+00 | 000000 | 9 | [ momentumHM : 3.84e+00, trackLength : 2.85e+03, timeOfFlight : 9.58e+00,] + | 0 | 0.0000e+00 | 000000 | 10 | [ momentumHM : 3.84e+00, trackLength : 2.62e+03, timeOfFlight : 0.00e+00,] + | 0 | 0.0000e+00 | 000000 | 11 | [ momentumHM : 3.84e+00, trackLength : 2.85e+03, timeOfFlight : 9.57e+00,] +. . . +``` +## Extracting output parameters + +Here is the pseudocode example on how one could extract PIDHandler parameters of *MyTofClosest0ps* results from the slcio file in your own Marlin processor and then use them to e.g. reconstruct the mass:
+``` +void YourAmazingProcessor::processEvent(LCEvent* event){ + + LCCollection* pfos = event->getCollection("PandoraPFOs"); + PIDHandler pidHandler(pfos); + + int algorithmID = pidHandler.getAlgorithmID("MyTofClosest0ps"); + int momParIdx = pidHandler.getParameterIndex(algorithmID, "momentumHM"); + int lenParIdx = pidHandler.getParameterIndex(algorithmID, "trackLength"); + int tofParIdx = pidHandler.getParameterIndex(algorithmID, "timeOfFlight"); + + + for(int i=0; i < pfos.getNumberOfElements(); ++i){ + ReconstructedParticle* pfo = static_cast ( pfos->getElementAt(i) ); + + const ParticleID& pfoID = pidHandler.getParticleID(pfo, algorithmID); + const std::vector& pars = pfoID.getParameters(); + + double momentum = pars[momParIdx]; // in GeV + double trackLength = pars[lenParIdx]; // in mm + double tof = pars[tofParIdx]; // in ns + + //in GeV + double mass = momentum * std::sqrt( std::pow(tof*CLHEP::c_light/trackLength, 2) - 1 ); + } +} +``` + + +## Authors +- N.Weinhold, DESY, internship 2017
+- F.Gaede, DESY, 2018
+- B.Dudar, DESY, 2021
diff --git a/TimeOfFlight/include/TOFEstimators.h b/TimeOfFlight/include/TOFEstimators.h index cf669a49..635f5ac5 100644 --- a/TimeOfFlight/include/TOFEstimators.h +++ b/TimeOfFlight/include/TOFEstimators.h @@ -1,116 +1,101 @@ #ifndef TOFEstimators_h #define TOFEstimators_h 1 +/** +\author F. Gaede, DESY, April 2018 +\author B. Dudar, DESY, September 2021 +*/ + #include #include #include "marlin/Processor.h" #include "MarlinTrk/IMarlinTrkSystem.h" -/** - * @section DESCRIPTION - * Compute time of flight (ToF) from the CalorimeterHits in the ECAL Cluster. - * The ToFs are stored in a PID object with the name of the processor attached - * to the PandoraPFOs. - * Steering parameters: - * MaxLayerNumber - Estimate ToF using only MaxLayerNumber first layers in the ECAL - * ReconstructedParticleCollection - input collection (usually PandoraPFOs) - * TimeResolution - assumed single Calorimeter hit time resolution in ps - * CylRadius - radius within which hits are accepted for the fit. Relevant only - for new TOFCylFit algorithm. Default 5 mm - * ProcessorVersion - controls output of the processor between idr - * (interim design report) and dev (current development) versions - * Output for ProcessorVersion = "idr" - * The following parameters are stored with the IDR ProcessorVersion: - * TOFFirstHit - ToF of the CalorimeterHit closest to track entry - * point to the ECAL corrected with the distance - * between cell center and entry point assuming - * speed of light propagation - * TOFClosestHits - ToF estimated as the average of corrected - * CalorimeterHit times from hits that are closest to - * the linear prolongation of the trajectory line - * of the track in MaxLayerNumber first layers - * TOFClosestHitsError - error of the TOFClosestHits - * TOFFlightLength - Track length of the helix between IP and ECAL entrance - * assuming IP at (0, 0, 0) and ECAL entry point is - * obtained from Kalman filter by extrapolation. - * Omega/lambda track parameters are taken from IP hit - * and assumed to be constant along the track. - * Results are reasonable for delta phi < 2pi - * meaning not more than one curl. - * TOFLastTrkHit - (unsmeared) time of last tracker hit (in the SET) - * TOFLastTrkHitFlightLength - track length similar to the TOFFlightLength - * but to the SET hit - * Notes: - * For charged particles the TrackState at the calorimeter is used as reference point and for the - * extrapolation into the calorimeter. For neutral particles the position of the first hit (the one closest to the IP) - * is used and a straight flight path from the IP is taken for the extrapolation. - * The hit time is corrected for the flight time from the reference assuming speed of light. - * - * Output for ProcessorVersion = "dev" - * The following parameters are stored with the IDR ProcessorVersion: - * TOFClosest - ToF of the CalorimeterHit closest to track entry - * point to the ECAL corrected with the distance - * between cell center and entry point assuming - * speed of light propagation - * TOFFastest - ToF of the CalorimeterHit arrived fastest corrected - * with the distance between cell center and entry - * point assuming speed of light propagation - * TOFCylFit - ToF estimated from extrapolation from a fit of the - * hit times vs distance to to the track entry point - * to the ECAL at distance = 0 - * TOFClosestFit - ToF estimated with a fit as TOFCylFit but from hits - * that are closest to the linear prolongation of the - * trajectory line of the track in MaxLayerNumber first layers - * FlightLength - Track length of the helix between IP and ECAL entrance - * assuming IP at (0, 0, 0) and ECAL entry point is - * obtained from Kalman filter by extrapolation. - * Omega/lambda track parameters are taken from the - * closest track hit to the ECAL and assumed to be - * constant along the track. - * Results are reasonable for delta phi < 2pi - * meaning not more than one curl. - * MomAtCalo - momentum of the track estimated from the track curvature near - * ECAL entry point. - * Notes: - * FlightLength shows less bias in mass than TOFFlightLength - * MomAtCalo shows less bias in mass than taking momentum from IP as in IDR - * Every method in both idr and dev considers only Ecal hits. - * Methods work only for charged PFOs with ONE associated track and ONE associated cluster - * Therefore neutral PFOs and PFOs with multiple tracks/clusters are ignored - * - * See ../scripts/tofestimators.xml for example steering file. - * - * @file - * @author F.Gaede, DESY, April 2018 - * @author B.Dudar, DESY, September 2021 - */ - class TOFEstimators : public marlin::Processor { public: + /** + Copy constructor is removed to avoid W-effc++ warnings. + */ TOFEstimators(const TOFEstimators&) = delete; + + /** + Copy operator is removed to avoid W-effc++ warnings. + */ TOFEstimators& operator=(const TOFEstimators&) = delete; + /** + Method required by the Marlin to register processor in the global scope. + */ marlin::Processor* newProcessor() { return new TOFEstimators; } + + /** + Default constructor. + Defines steering parameters from the xml steering file. + */ TOFEstimators(); + + /** Called at the begin of the job before anything is read. + Extracts geometry details and initializes Kalman Filter System. + */ void init(); + + /** Called for every event - the working horse. + Calculates momentum, track length and time of flight and + writes them into PIDHandler of the input collection object + */ void processEvent(EVENT::LCEvent* evt); + private: + /** Steering parameter: name of the input collection of ReconstructedParticles to analyse. + Usually PandoraPFOs + */ std::string _pfoCollectionName{}; + + /** Steering parameter: an option to indicate whether to do the calculations. + to the last tracker hit or extrapolate to the ECal surface. + */ bool _extrapolateToEcal{}; + + /** Steering parameter: name of the method to calculate time of flight based on the ECal hits. + If _extrapolateToEcal == false then ignored + */ std::string _tofMethod{}; + + /** Steering parameter: Assumed time resolution of the detector elements in ps. + In case _extrapolateToEcal == false, then smearing applied for both SET strips. + In case _extrapolateToEcal == true, then smearing applied for every ECal hit. + */ double _timeResolution{}; + + /** Steering parameter: Number of inner ECal layers to consider for TOF methods. + In case _extrapolateToEcal == true and (_tofMethod == "frankAvg" or _tofMethod == "frankFit"). + This parameter indicates how many ECal layers these tof Methods will consider to estimate TOF. + */ int _maxEcalLayer{}; + /** Number of the current event. + */ int _nEvent{}; - /** vector of names of output PID object parameters which are written as output + + /** Names of the output parameters written to the PIDHandler. + These parameters are: momentumHM, trackLength and timeOfFlight */ std::vector _outputParNames{}; - // MarlinTrk v02-00 release notes: - // USERS SHOULD NO LONGER DELETE THE IMarlinTrkSystem POINTER IN THEIR CODE + + /** Kalman Filter System for refitting. + Release notes of MarlinTrk v02-00: + USERS SHOULD NO LONGER DELETE THE IMarlinTrkSystem POINTER IN THEIR CODE + */ MarlinTrk::IMarlinTrkSystem* _trkSystem = nullptr; + + /** B_{z} at the origin in Tesla. + */ double _bField{}; + + /** TPC outer radius in mm. + */ double _tpcOuterR{}; }; diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index 08f5335e..b3968a01 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -1,12 +1,11 @@ #ifndef TOFUtils_h #define TOFUtils_h 1 -/****************************************************** - * Utility functions that are used by the TOFEstimators processor +/** + * Utility functions that are used by the TOFEstimators processor. * - * @author F. gaede, DESY, 2018 - * @author B. Dudar, DESY, 2021 - ****************************************************** + * \author F. Gaede, DESY, 2018 + * \author B. Dudar, DESY, 2021 */ #include @@ -18,49 +17,100 @@ namespace TOFUtils{ + /** Comparator function for tracker hits. + + Returns `true` if \f$ \rho_{a} < \rho_{b} \f$. + With \f$ \rho \f$ being a radius of the tracker hit projected on the \f$ xy \f$ plane: \f$ \rho = \sqrt{x^2 + y^2} \f$. + + Primarily used to sort vector of tracker hits by \f$ \rho \f$ for the Kalman Filter fit. + */ bool sortByRho(EVENT::TrackerHit* a, EVENT::TrackerHit* b); + + /** Extracts track state at the given tracker hit used by Kalman Filter in the refit. + */ IMPL::TrackStateImpl getTrackStateAtHit(MarlinTrk::IMarlinTrack* marlinTrk, EVENT::TrackerHit* hit); - // gets momentum of the helix from parameters of the trackState + /** Returns momentum vector from the given track state. + */ dd4hep::rec::Vector3D getHelixMomAtTrackState(const EVENT::TrackState& ts, double bField); + /** Estimate helix arc length between two track states. + + Works only for arcs with \Delta \varphi < \pi + Here is the formula it uses + + */ double getHelixArcLength(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); + + /** Estimate helix arc length between two track states without \varphi. + + Relies on tanL. Less precise than the getHelixMomAtTrackState. + Used to calculate track length with \Delta \varphi > \pi + + Here is the formula it uses + + */ double getHelixLengthAlongZ(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); + + /** Get number of helix revolutions between two track states assuming constant momentum. + */ double getHelixNRevolutions(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); + /** Return TPC outer radius from DD4hep detector geometry. + */ double getTPCOuterR(); + /** Return SET hit. + SET hit is assumed to be a tracker hit with a radius larger than TPC outer radius + */ EVENT::TrackerHit* getSETHit(EVENT::Track* track, double tpcOuterR); + /** TEST. + SET hit is assumed to be a tracker hit with a radius larger than TPC outer radius + */ std::vector selectFrankEcalHits( EVENT::Cluster* cluster, EVENT::Track* track, int maxEcalLayer, double bField ); - // This is not guarantied to work 100% of times, - // but I didn't find a better way to collect subTracks... - // 1) Always add main track as a first subtrack - // 2) If nTPCHits == nHits of SubTrack0 +-1 assume subTrack0 stores TPC hits - // Deviation on 1 hit may happen because of the SET and only God knows what other reasons.. - // This would mean that VXD+SIT subTrack is not stored in the track! - // So we skip initial subTrack0 with TPC hits which we have anyhow added with global Track - // and loop over remaining subTracks - // 3) If 2) is false and nTPCHits == nHits of SubTrack1 +-1 assume subTrack1 stores TPC hits - // This would mean that subTrack0 stores VXD+SIT hits - // So we skip both VXD+SIT and 1st TPC subTracks which we have anyhow added with global Track - // and loop over remaining subTracks + + /** TEST. + SET hit is assumed to be a tracker hit with a radius larger than TPC outer radius + + This is not guarantied to work 100% of times, + but I didn't find a better way to collect subTracks... + 1) Always add main track as a first subtrack + 2) If nTPCHits == nHits of SubTrack0 +-1 assume subTrack0 stores TPC hits + Deviation on 1 hit may happen because of the SET and only God knows what other reasons.. + This would mean that VXD+SIT subTrack is not stored in the track! + So we skip initial subTrack0 with TPC hits which we have anyhow added with global Track + and loop over remaining subTracks + 3) If 2) is false and nTPCHits == nHits of SubTrack1 +-1 assume subTrack1 stores TPC hits + This would mean that subTrack0 stores VXD+SIT hits + So we skip both VXD+SIT and 1st TPC subTracks which we have anyhow added with global Track + and loop over remaining subTracks + + */ std::vector getSubTracks(EVENT::Track* track); - // return vector of trackStates for the IP, every hit in the whole track and ECAL if extrapolate to the ECAL option is chosen - // track state are sorted along the helix + /** test. + return vector of trackStates for the IP, every hit in the whole track and ECAL if extrapolate to the ECAL option is chosen + track state are sorted along the helix + */ std::vector getTrackStatesPerHit(std::vector tracks, MarlinTrk::IMarlinTrkSystem* trkSystem, bool extrapolateToEcal, double bField); + /** test. + */ double getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, double timeResolution); + /** test. + */ double getTofFrankAvg( std::vector selectedHits, EVENT::Track* track, double timeResolution); + /** test. + */ double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); diff --git a/TimeOfFlight/xml/steer.xml b/TimeOfFlight/xml/steer.xml index 4b1dfbc8..e909cee5 100644 --- a/TimeOfFlight/xml/steer.xml +++ b/TimeOfFlight/xml/steer.xml @@ -3,40 +3,90 @@ + + + + + + + - + /pnfs/desy.de/ilc/prod/ilc/mc-2020/ild/rec/250-SetA/2f_hadronic_eL_pR/ILD_l5_o1_v02/v02-02/00015161/000/rv02-02.sv02-02.mILD_l5_o1_v02.E250-SetA.I500010.P2f_z_h.eL.pR.n001_209.d_rec_00015161_259.slcio - - - - - MESSAGE DEBUG0 - + + + 2 + 0 + false + false + MESSAGE + 1234567890 - + /cvmfs/ilc.desy.de/sw/x86_64_gcc82_centos7/v02-02-02/lcgeo/v00-16-06/ILD/compact/ILD_l5_o1_v02/ILD_l5_o1_v02.xml - true - closest - 0 + + + PandoraPFOs + + + true + + + closest + + + 0 + + + 10 + + + + + + + false + + + + 10 + + + + + + frankAvg + + + + 49.999 + - test.slcio + + + + MCParticle SimTrackerHit SimCalorimeterHit TrackerHit TrackerHitPlane CalorimeterHit LCRelation Track LCFloatVec + + + PandoraPFOs + + output.slcio WRITE_NEW - * From c9fca9425ecea0c2ccc1aa2c410c4db06e6e95d7 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Fri, 24 Sep 2021 15:04:31 +0200 Subject: [PATCH 05/22] Finilize documentation. --- TimeOfFlight/README.md | 201 ++++++++++++--------------- TimeOfFlight/include/TOFEstimators.h | 47 +++---- TimeOfFlight/include/TOFUtils.h | 111 ++++++++++----- 3 files changed, 188 insertions(+), 171 deletions(-) diff --git a/TimeOfFlight/README.md b/TimeOfFlight/README.md index 90cb5281..c4379f65 100644 --- a/TimeOfFlight/README.md +++ b/TimeOfFlight/README.md @@ -6,7 +6,7 @@ In addition it calculates square root of the harmonic mean of the momentum squar [Doxygen documentation](https://www.desy.de/~dudarboh/tof_doc/html/index.html) contains a detailed description of the source code.
In case reader uses dark theme on the github and finds hard to read the formulas he can find this README as the main page there. -## Requirements +# Requirements TOFEstimators is dependent on InitDD4hep processor to extract detector geometry details.
Make sure to run InitDD4hep before this processor is executed. @@ -14,57 +14,62 @@ Make sure to run InitDD4hep before this processor is executed. In addition, input sclio file must contain a lot of various tracker and calorimeter hit collections.
Thus, input slcio file can be e.g. *REC* format, but **cannot** be *DST* or *mini-DST* format. -## Steering parameters +# Steering parameters -1. **ReconstructedParticleCollection** / default: `"PandoraPFOs"`
+There are five steering parameters: + +### ReconstructedParticleCollection / default: `"PandoraPFOs"` Collection of `ReconstructedParticle` objects to analyze.
The final results will be written in the PIDHandler inside this collection. -2. **ExtrapolateToEcal** / default: `true`
- If `true` the track length is calculated to the extrapolated track position at the ECal surface. Time-of-flight then calculated using ECal hits via one of the methods chosen in *TofMethod* steering parameter (see below). +### ExtrapolateToEcal / default: `true`
+ If `true` the track length is calculated to the extrapolated track position at the ECal surface. Time-of-flight then calculated using ECal hits via one of the methods chosen in *TofMethod* steering parameter. If `false` the track length is calculated to the last tracker hit. Time-of-flight then calculated using the last tracker hit. Which is SET hit for the barrel region and *undefined* for the endcap region.
If the last tracker hit is the SET hit time-of-flight is set to the average time of two SET strips.
- If the last tracker hit is the TPC hit time-of-flight is set to `0`. + If the last tracker hit is not the SET hit time-of-flight is set to `0`. + +### TofMethod / options: `"closest"`, `"frankAvg"`, `"frankFit"` / default: `"closest"` + + If *ExtrapolateToEcal* is set to `false` then this steering parameter is ignored. -3. **TofMethod** / options: `"closest"`, `"frankAvg"`, `"frankFit"` / default: `"closest"`
- If *ExtrapolateToEcal* is set to `true` then it defines how to calculate the time-of-flight at the ECal surface. + If *ExtrapolateToEcal* is set to `true` then it defines how to use ECal hits to calculate the time-of-flight at the ECal surface. - `"closest"` uses the closest ECal hit to the extrapolated track position at the ECal surface.
- Time-of-flight is defined as: . + Time-of-flight is defined as: - `"frankAvg"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Define time-of-flight as an average of their time corrected for the distance to the track position at the ECal surface assuming speed of flight is the speed of light. - . + - - `"frankFit"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Use linear fit of the hits time as a function of the distance to the track position at the ECal surface to define time-of-flight as an extrapolated time at the ECal surface . + - `"frankFit"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Use linear fit of the hits time as a function of the distance to the track position at the ECal surface to define time-of-flight as an extrapolated time at the extrapolated track position at the ECal surface for more illustrative explanations of the methods above see slides 10, 12, 13 from the following [LCWS2021 talk]((https://indico.cern.ch/event/995633/contributions/4259659/attachments/2209010/3738157/Bohdan_TOF_LCWS2021.pdf)). - If *ExtrapolateToEcal* is set to `false` then this steering parameter is ignored. +### TimeResolution / default: `0.0` (ps) + If *ExtrapolateToEcal* is set to `true` it defines times resolution of individual ECal hits in ps. -4. **TimeResolution** / default: `0.0`
- If *ExtrapolateToEcal* is set to `true` it defines times resolution of individual ECal hits in ps. + If *ExtrapolateToEcal* is set to `false` it defines times resolution of individual SET strips in ps. - If *ExtrapolateToEcal* is set to `false` it defines times resolution of individual SET strips in ps. +### MaxEcalLayer / default: `10` + Defines number of inner ECal layers to use for the *frankAvg* and *frankFit* algorithms. -5. **MaxEcalLayer** / default: `10`
- Defines number of inner ECal layers to use for the *frankAvg* and *frankFit* algorithms. + If *frankAvg* and *frankFit* are unused then this steering parameter is ignored. - If *frankAvg* and *frankFit* are unused then this steering parameter is ignored. +# Output parameters -## Output parameters +There are three output parameters: -1. `momentumHM` +### `momentumHM` (GeV/c) - A squared root of the harmonic mean of the squared momentum: . + A squared root of the harmonic mean of the squared momentum: - If momentum is constant, then it is just a momentum: .
- If momentum changes, then using instead of or is mathematically rigorous choice for the relativistic particle assumption .
+ If momentum is constant, then it is just a momentum:
+ If momentum changes, then using instead of or is mathematically rigorous choice for the relativistic particle assumption
See [Winfried A. Mitaroff paper](https://arxiv.org/abs/2107.02031) for the details. -2. `trackLength` +### `trackLength` (mm) A track length calculated between the point of the closest approach (PCA) to the IP `(0,0,0)`
and the ECal surface if *ExtrapolateToEcal* is `true`
@@ -72,21 +77,17 @@ Thus, input slcio file can be e.g. *REC* format, but **cannot** be *DST* or *min A total track length is calculated as a sum of track segments obtained from iterating over track states defined at the IP, every tracker hit and the ECal surface if *ExtrapolateToEcal* is `true`.
A track segment length is calculated between neighboring track states as: - . + If *ExtrapolateToEcal* is `true` then we additionally check for the number of helix turns between the last tracker hit and the extrapolated track position at the ECal surface. If we are unable to use the formula above to calculate the length of the last segment. Thus we use different formula that does not rely on the azimuthal angle: - . - + -3. `timeOfFlight` +### `timeOfFlight` (ns) Time-of-flight of the particle assuming . For the detailed description of different time-of-flight estimators see the **TofMethod** bullet point in the Steering parameters section. - - - -## Steering file example +# Steering file example [./xml/steer.xml](./xml/steer.xml) is a steering file example that runs three TOFEstimators processors: *MyTofClosest0ps*, *MyTofSET10ps* and *MyTofFrankAvg50ps* and then writes their output parameters in the PIDHandlers of the *PandoraPFOs* collection in the new *output.slcio* file using LCIOOutputProcessor. @@ -102,86 +103,68 @@ Then one can look at the output.sclio file with, e.g.:
Scrolling carefully one can find these lines which indicate that the new file contains our parameters for all three TOFEstimators: -``` -collection name : PandoraPFOs -parameters: - ---------------- print out of ReconstructedParticle collection --------------- - - flag: 0x0 -. . . -parameter ParameterNames_MyTofClosest0ps [string]: momentumHM, trackLength, timeOfFlight, -parameter ParameterNames_MyTofFrankAvg50ps [string]: momentumHM, trackLength, timeOfFlight, -parameter ParameterNames_MyTofSET10ps [string]: momentumHM, trackLength, timeOfFlight, -. . . -``` - -Scrolling more down, one can see calculated parameters for each individual particle:
-``` ------------- detailed PID info: --- - - algorithms : - [id: 0] BasicVariablePID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood M - [id: 3] LikelihoodPID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MVAO - [id: 4] LowMomMuID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MVAOutp - [id: 9] MyTofClosest0ps - params: momentumHM trackLength timeOfFlight - [id: 11] MyTofFrankAvg50ps - params: momentumHM trackLength timeOfFlight - [id: 10] MyTofSET10ps - params: momentumHM trackLength timeOfFlight - [id: 2] ShowerShapesPID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MV - [id: 5] TOFEstimators0ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightLen - [id: 8] TOFEstimators100ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightL - [id: 6] TOFEstimators10ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightLe - [id: 7] TOFEstimators50ps - params: TOFFirstHit TOFClosestHits TOFClosestHitsError TOFFlightLength TOFLastTrkHit TOFLastTrkHitFlightLe - [id: 1] dEdxPID - params: electronLikelihood muonLikelihood pionLikelihood kaonLikelihood protonLikelihood hadronLikelihood MVAOutput_ - - - [particle] | PDG | likelihood | type | algoId | parameters : - | | | | | -. . . - [00000073] | 2212 | 2.0802e+00 | 000000 | 1 | [ electronLikelihood : -1.64e+02, muonLikelihood : -3.47e+01, pionLikelihood : -2. - | 0 | 0.0000e+00 | 000000 | 4 | [ electronLikelihood : 9.99e+02, muonLikelihood : 9.99e+02, pionLikelihood : 9.99e - | 0 | 0.0000e+00 | 000000 | 5 | [ TOFFirstHit : 9.58e+00, TOFClosestHits : 9.58e+00, TOFClosestHitsError : 2.71e-0 - | 0 | 0.0000e+00 | 000000 | 6 | [ TOFFirstHit : 9.57e+00, TOFClosestHits : 9.58e+00, TOFClosestHitsError : 7.16e-0 - | 0 | 0.0000e+00 | 000000 | 7 | [ TOFFirstHit : 9.54e+00, TOFClosestHits : 9.58e+00, TOFClosestHitsError : 6.32e-0 - | 0 | 0.0000e+00 | 000000 | 8 | [ TOFFirstHit : 9.66e+00, TOFClosestHits : 9.64e+00, TOFClosestHitsError : 8.34e-0 - | 0 | 0.0000e+00 | 000000 | 9 | [ momentumHM : 3.84e+00, trackLength : 2.85e+03, timeOfFlight : 9.58e+00,] - | 0 | 0.0000e+00 | 000000 | 10 | [ momentumHM : 3.84e+00, trackLength : 2.62e+03, timeOfFlight : 0.00e+00,] - | 0 | 0.0000e+00 | 000000 | 11 | [ momentumHM : 3.84e+00, trackLength : 2.85e+03, timeOfFlight : 9.57e+00,] -. . . -``` -## Extracting output parameters - -Here is the pseudocode example on how one could extract PIDHandler parameters of *MyTofClosest0ps* results from the slcio file in your own Marlin processor and then use them to e.g. reconstruct the mass:
-``` -void YourAmazingProcessor::processEvent(LCEvent* event){ - - LCCollection* pfos = event->getCollection("PandoraPFOs"); - PIDHandler pidHandler(pfos); - - int algorithmID = pidHandler.getAlgorithmID("MyTofClosest0ps"); - int momParIdx = pidHandler.getParameterIndex(algorithmID, "momentumHM"); - int lenParIdx = pidHandler.getParameterIndex(algorithmID, "trackLength"); - int tofParIdx = pidHandler.getParameterIndex(algorithmID, "timeOfFlight"); - - - for(int i=0; i < pfos.getNumberOfElements(); ++i){ - ReconstructedParticle* pfo = static_cast ( pfos->getElementAt(i) ); - - const ParticleID& pfoID = pidHandler.getParticleID(pfo, algorithmID); - const std::vector& pars = pfoID.getParameters(); - - double momentum = pars[momParIdx]; // in GeV - double trackLength = pars[lenParIdx]; // in mm - double tof = pars[tofParIdx]; // in ns - - //in GeV - double mass = momentum * std::sqrt( std::pow(tof*CLHEP::c_light/trackLength, 2) - 1 ); + + collection name : PandoraPFOs + parameters: + --------------- print out of ReconstructedParticle collection --------------- + . . . + parameter ParameterNames_MyTofClosest0ps [string]: momentumHM, trackLength, timeOfFlight, + parameter ParameterNames_MyTofFrankAvg50ps [string]: momentumHM, trackLength, timeOfFlight, + parameter ParameterNames_MyTofSET10ps [string]: momentumHM, trackLength, timeOfFlight, + . . . + + +Scrolling more down, one can see calculated parameters for each individual particle: + + + ------------ detailed PID info: --- + algorithms : + . . . + [id: 9] MyTofClosest0ps - params: momentumHM trackLength timeOfFlight + [id: 11] MyTofFrankAvg50ps - params: momentumHM trackLength timeOfFlight + [id: 10] MyTofSET10ps - params: momentumHM trackLength timeOfFlight + . . . + [particle] | PDG | likelihood | type | algoId | parameters : + | | | | | + [00000073] . . . + | 0 | 0.0000e+00 | 000000 | 9 | [ momentumHM : 3.84e+00, trackLength : 2.85e+03, timeOfFlight : 9.58e+00,] + | 0 | 0.0000e+00 | 000000 | 10 | [ momentumHM : 3.84e+00, trackLength : 2.62e+03, timeOfFlight : 0.00e+00,] + | 0 | 0.0000e+00 | 000000 | 11 | [ momentumHM : 3.84e+00, trackLength : 2.85e+03, timeOfFlight : 9.57e+00,] + . . . + +# Extracting output parameters + +Here is the pseudocode example on how one could extract PIDHandler parameters of *MyTofClosest0ps* results from the slcio file in your own Marlin processor and then use them to e.g. reconstruct the mass: + + + void YourAmazingProcessor::processEvent(LCEvent* event){ + + LCCollection* pfos = event->getCollection("PandoraPFOs"); + PIDHandler pidHandler(pfos); + + int algorithmID = pidHandler.getAlgorithmID("MyTofClosest0ps"); + int momParIdx = pidHandler.getParameterIndex(algorithmID, "momentumHM"); + int lenParIdx = pidHandler.getParameterIndex(algorithmID, "trackLength"); + int tofParIdx = pidHandler.getParameterIndex(algorithmID, "timeOfFlight"); + + + for(int i=0; i < pfos.getNumberOfElements(); ++i){ + ReconstructedParticle* pfo = static_cast ( pfos->getElementAt(i) ); + + const ParticleID& pfoID = pidHandler.getParticleID(pfo, algorithmID); + const std::vector& pars = pfoID.getParameters(); + + double momentum = pars[momParIdx]; // GeV + double trackLength = pars[lenParIdx]; // mm + double tof = pars[tofParIdx]; // ns + + //in GeV + double mass = momentum * std::sqrt( std::pow(tof*CLHEP::c_light/trackLength, 2) - 1 ); + } } -} -``` -## Authors +# Authors - N.Weinhold, DESY, internship 2017
- F.Gaede, DESY, 2018
- B.Dudar, DESY, 2021
diff --git a/TimeOfFlight/include/TOFEstimators.h b/TimeOfFlight/include/TOFEstimators.h index 635f5ac5..c5e50847 100644 --- a/TimeOfFlight/include/TOFEstimators.h +++ b/TimeOfFlight/include/TOFEstimators.h @@ -2,6 +2,7 @@ #define TOFEstimators_h 1 /** +Marlin processor than calculates output parameters. \author F. Gaede, DESY, April 2018 \author B. Dudar, DESY, September 2021 */ @@ -14,12 +15,16 @@ class TOFEstimators : public marlin::Processor { public: /** - Copy constructor is removed to avoid W-effc++ warnings. + Copy constructor . + We remove it to avoid W-effc++ warnings. + Copying objects with pointer members is a bad idea. */ TOFEstimators(const TOFEstimators&) = delete; /** - Copy operator is removed to avoid W-effc++ warnings. + Copy assignment operator. + We remove it to avoid W-effc++ warnings. + Copying objects with pointer members is a bad idea. */ TOFEstimators& operator=(const TOFEstimators&) = delete; @@ -30,7 +35,7 @@ class TOFEstimators : public marlin::Processor { /** Default constructor. - Defines steering parameters from the xml steering file. + Registers steering parameters from the xml steering file. */ TOFEstimators(); @@ -41,60 +46,52 @@ class TOFEstimators : public marlin::Processor { /** Called for every event - the working horse. Calculates momentum, track length and time of flight and - writes them into PIDHandler of the input collection object + writes them into PIDHandler of the input collection object. */ void processEvent(EVENT::LCEvent* evt); private: - /** Steering parameter: name of the input collection of ReconstructedParticles to analyse. - Usually PandoraPFOs + /** Stores ReconstructedParticleCollection steering parameter. */ std::string _pfoCollectionName{}; - /** Steering parameter: an option to indicate whether to do the calculations. - to the last tracker hit or extrapolate to the ECal surface. + /** Stores ExtrapolateToEcal steering parameter. */ bool _extrapolateToEcal{}; - /** Steering parameter: name of the method to calculate time of flight based on the ECal hits. - If _extrapolateToEcal == false then ignored + /** Stores TofMethod steering parameter. */ std::string _tofMethod{}; - /** Steering parameter: Assumed time resolution of the detector elements in ps. - In case _extrapolateToEcal == false, then smearing applied for both SET strips. - In case _extrapolateToEcal == true, then smearing applied for every ECal hit. + /** Stores TimeResolution steering parameter. */ double _timeResolution{}; - /** Steering parameter: Number of inner ECal layers to consider for TOF methods. - In case _extrapolateToEcal == true and (_tofMethod == "frankAvg" or _tofMethod == "frankFit"). - This parameter indicates how many ECal layers these tof Methods will consider to estimate TOF. + /** Stores MaxEcalLayer steering parameter. */ int _maxEcalLayer{}; - - /** Number of the current event. + /** Stores current event number. */ int _nEvent{}; - /** Names of the output parameters written to the PIDHandler. - These parameters are: momentumHM, trackLength and timeOfFlight + /** Stores names of the output parameters. + These are "momentumHM", "trackLength" and "timeOfFlight". */ std::vector _outputParNames{}; - /** Kalman Filter System for refitting. - Release notes of MarlinTrk v02-00: - USERS SHOULD NO LONGER DELETE THE IMarlinTrkSystem POINTER IN THEIR CODE + /** Kalman Filter System. + \note Release notes of MarlinTrk v02-00: + \note users should no longer delete the IMarlinTrkSystem pointer in their code */ MarlinTrk::IMarlinTrkSystem* _trkSystem = nullptr; - /** B_{z} at the origin in Tesla. + /** Stores z component of the magnetic field at the origin in Tesla. */ double _bField{}; - /** TPC outer radius in mm. + /** Stores outer TPC radius in mm. */ double _tpcOuterR{}; }; diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index b3968a01..9f21b2cd 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -17,8 +17,7 @@ namespace TOFUtils{ - /** Comparator function for tracker hits. - + /** Comparator function by radius for tracker hits. Returns `true` if \f$ \rho_{a} < \rho_{b} \f$. With \f$ \rho \f$ being a radius of the tracker hit projected on the \f$ xy \f$ plane: \f$ \rho = \sqrt{x^2 + y^2} \f$. @@ -27,89 +26,127 @@ namespace TOFUtils{ bool sortByRho(EVENT::TrackerHit* a, EVENT::TrackerHit* b); - /** Extracts track state at the given tracker hit used by Kalman Filter in the refit. + /** Get track state at tracker hit. + Returns track state at tracker hit from the Kalman Filter fit. + + Tracker hit has to be used in the fit by the Kalman Filter. */ IMPL::TrackStateImpl getTrackStateAtHit(MarlinTrk::IMarlinTrack* marlinTrk, EVENT::TrackerHit* hit); - /** Returns momentum vector from the given track state. + /** Get track momentum at the track state. + Returns momentum Vector3D from the given track state. */ dd4hep::rec::Vector3D getHelixMomAtTrackState(const EVENT::TrackState& ts, double bField); - /** Estimate helix arc length between two track states. + /** Get track length. + Returns the track length between two track states estimated by the helix length formula: - Works only for arcs with \Delta \varphi < \pi - Here is the formula it uses + \f$ \ell = \sqrt{\left( \frac{\varphi_{i+1} - \varphi_{i}}{\Omega}\right)^{2} + \left( z_{i+1} - z_{i} \right)^{2} } \f$ + Note: The formula above works only for the arcs with \f$ \Delta \varphi < \pi \f$. */ double getHelixArcLength(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); - /** Estimate helix arc length between two track states without \varphi. - Relies on tanL. Less precise than the getHelixMomAtTrackState. - Used to calculate track length with \Delta \varphi > \pi + /** Get track length. + Returns the track length between two track states estimated by the helix length formula: + + \f$ \ell = \frac{\left |z_{i+1} - z_{i}\right |}{\tan{\lambda}} \sqrt{1 + \tan{\lambda}^{2} } \f$ - Here is the formula it uses + Note: The formula above works for any \f$ \Delta \varphi \f$. + However it is less precise than getHelixArcLength() due to the less precise \f$ \tan{\lambda} \f$. + Also helix formula implies constant momentum assumption which may show higher discrepancy for long tracks with low momentum. */ double getHelixLengthAlongZ(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); - /** Get number of helix revolutions between two track states assuming constant momentum. + /** Get number of helix revolutions. + Returns number of helix revolutions between two track states. + + The calculation is done with: + \f$ N_{\mathrm{turns}} = \frac{\left |z_{i+1} - z_{i}\right |}{\tan{\lambda}} \bigg / (2 \pi \frac{1}{|\Omega|}) \f$ */ double getHelixNRevolutions(const EVENT::TrackState& ts1, const EVENT::TrackState& ts2); - /** Return TPC outer radius from DD4hep detector geometry. + /** Returns TPC outer radius from the DD4hep detector geometry. */ double getTPCOuterR(); - /** Return SET hit. - SET hit is assumed to be a tracker hit with a radius larger than TPC outer radius + /** Returns SET hit. + Checks the \f$ \rho = \sqrt{x^{2}+y^{2}} \f$ of the tracker hit. + If \f$ \rho > R_{\mathrm{TPC, outer}} \f$ then it is the SET hit. */ EVENT::TrackerHit* getSETHit(EVENT::Track* track, double tpcOuterR); - /** TEST. - SET hit is assumed to be a tracker hit with a radius larger than TPC outer radius + /** Get a subset of the cluster calorimeter hits. + Returns the subset of the cluster calorimeter hits which we call "Frank" + by historical reasons as these hits were used to calculate time-of-flight as + a default method for the IDR production. + + "Frank" hits are the closest ECal hits to the linearly extrapolated + track line inside the ECal in each of the first maxEcalLayer ECal layers. */ std::vector selectFrankEcalHits( EVENT::Cluster* cluster, EVENT::Track* track, int maxEcalLayer, double bField ); - /** TEST. - SET hit is assumed to be a tracker hit with a radius larger than TPC outer radius + /** Get all subtracks of the Track. + Returns a vector of the subTracks of the main track that is passed as an argument. + + The main purpose of this function is to capture all hits of the track. + This requires iteration over the subTracks as the main Track stores hits only + of the first half turn due to the our fit procedure. - This is not guarantied to work 100% of times, - but I didn't find a better way to collect subTracks... - 1) Always add main track as a first subtrack - 2) If nTPCHits == nHits of SubTrack0 +-1 assume subTrack0 stores TPC hits - Deviation on 1 hit may happen because of the SET and only God knows what other reasons.. - This would mean that VXD+SIT subTrack is not stored in the track! - So we skip initial subTrack0 with TPC hits which we have anyhow added with global Track - and loop over remaining subTracks - 3) If 2) is false and nTPCHits == nHits of SubTrack1 +-1 assume subTrack1 stores TPC hits - This would mean that subTrack0 stores VXD+SIT hits - So we skip both VXD+SIT and 1st TPC subTracks which we have anyhow added with global Track - and loop over remaining subTracks + The first subTrack in the returned vector is always main Track itself, which + containes VXD, SIT and TPC tracker hits for the first half turn. + Then additional subTracks are added which contain hits for additional track revolutions if such exist. + + If nTPCHits \f$ \pm 1\f$ = nHits of SubTrack0 then assume subTrack0 stores TPC hits. + This would mean that VXD and SIT subTrack is not stored! + So we need skip *only first* subtrack: initial subTrack0 with TPC hits which we have anyhow added with the main Track + Else if nTPCHits \f$ \pm 1\f$ = nHits of SubTrack1 assume subTrack1 stores TPC hits + This would mean that subTrack0 stores VXD and SIT hits. + So we need to *skip both* subtracks VXD+SIT and 1st TPC half-turn subTracks which we have anyhow added with the main Track + + Note: We consider deviations for \f$\pm 1\f$ hit may happen because of the + SET and only God knows what other reasons... + + Note: This function is not guarantied to properly work 100% of times, + but I didn't find a better way to collect subTracks. */ std::vector getSubTracks(EVENT::Track* track); - /** test. - return vector of trackStates for the IP, every hit in the whole track and ECAL if extrapolate to the ECAL option is chosen - track state are sorted along the helix + /** Get list of track states. + Returns a vector of track states at the IP, track state for every tracker hit + inside all provided subTracks as tracks argument and the track state at the ECal surface if extrapolateToEcal argument is set to true. */ std::vector getTrackStatesPerHit(std::vector tracks, MarlinTrk::IMarlinTrkSystem* trkSystem, bool extrapolateToEcal, double bField); - /** test. + /** Get the time-of-flight using the closest ECal hit. + Returns time measured by the closest ECal hit to the extrapolated track position at the ECal surface. + + \f$ \mathrm{TOF} = t_{\mathrm{closest}} - \frac{\left| \vec{r}_{\mathrm{track}} - \vec{r}_{\mathrm{closest}} \right|}{c} \f$ */ double getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, double timeResolution); - /** test. + /** Get the time-of-flight using average of the Frank ECal hits. + Returns the time-of-flight as an average of the Frank hits time corrected for the distance to + the extrapolated track position at the ECal surface assuming speed of flight is the speed of light. + + \f$ \mathrm{TOF} = \frac{1}{\mathrm{MaxEcalLayer}}\sum_{i}^{\mathrm{MaxEcalLayer}} \left( t_{i} - \frac{\left|\vec{r}_{\mathrm{track}} - \vec{r}_{i} \right|}{c} \right) \f$ */ double getTofFrankAvg( std::vector selectedHits, EVENT::Track* track, double timeResolution); - /** test. + /** Get the time-of-flight using fit of the Frank ECal hits. + Returns the time-of-flight as an extrapolated time at the extrapolated track + position at the ECal surface by using a linear fit of the Frank hits time as a + function of the distance to the extrapolated track position at the ECal surface. + + \f$ \mathrm{TOF}=f(|\vec{r}_{\mathrm{track}} - \vec{r}_{\mathrm{hit}} |=0) \f$ */ double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); From b169a093e125335defbff584122139f28b435839 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Fri, 24 Sep 2021 15:55:47 +0200 Subject: [PATCH 06/22] Finalize documentation v2.0 --- TimeOfFlight/README.md | 21 +++++++++++++++------ TimeOfFlight/include/TOFEstimators.h | 2 +- TimeOfFlight/include/TOFUtils.h | 7 +++++++ TimeOfFlight/src/TOFUtils.cc | 2 ++ 4 files changed, 25 insertions(+), 7 deletions(-) diff --git a/TimeOfFlight/README.md b/TimeOfFlight/README.md index c4379f65..fdcccd64 100644 --- a/TimeOfFlight/README.md +++ b/TimeOfFlight/README.md @@ -8,12 +8,16 @@ In case reader uses dark theme on the github and finds hard to read the formulas # Requirements -TOFEstimators is dependent on InitDD4hep processor to extract detector geometry details.
-Make sure to run InitDD4hep before this processor is executed. +TOFEstimators is currently working with only simple cases of particle flow objects.
+The analyzed *ReconstructedParticle* object must have **exactly one** associated `Track` and **exactly one** associated `Cluster`.
+In other cases it still produces PIDHandler, but with `0.0` for every output parameter. -In addition, input sclio file must contain a lot of various tracker and calorimeter hit collections.
+The input sclio file must contain a lot of various tracker and calorimeter hit collections.
Thus, input slcio file can be e.g. *REC* format, but **cannot** be *DST* or *mini-DST* format. +TOFEstimators is dependent on InitDD4hep processor to extract detector geometry details.
+Make sure to run InitDD4hep before this processor is executed. + # Steering parameters There are five steering parameters: @@ -36,12 +40,17 @@ There are five steering parameters: If *ExtrapolateToEcal* is set to `true` then it defines how to use ECal hits to calculate the time-of-flight at the ECal surface. - `"closest"` uses the closest ECal hit to the extrapolated track position at the ECal surface.
- Time-of-flight is defined as: + Time-of-flight is defined as:
+ If no ECal hits are found returns `0.0`. - `"frankAvg"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Define time-of-flight as an average of their time corrected for the distance to the track position at the ECal surface assuming speed of flight is the speed of light. - +
+ If no ECal hits are found returns `0.0`. + - - `"frankFit"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Use linear fit of the hits time as a function of the distance to the track position at the ECal surface to define time-of-flight as an extrapolated time at the extrapolated track position at the ECal surface + - `"frankFit"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Use linear fit of the hits time as a function of the distance to the track position at the ECal surface to define time-of-flight as an extrapolated time at the extrapolated track position at the ECal surface
+ If no ECal hits are found returns `0.0`.
+ If only one ECal hit is found, which is not enough to perform the linear fit, then returns the same as *closest* and *frankAvg*. for more illustrative explanations of the methods above see slides 10, 12, 13 from the following [LCWS2021 talk]((https://indico.cern.ch/event/995633/contributions/4259659/attachments/2209010/3738157/Bohdan_TOF_LCWS2021.pdf)). diff --git a/TimeOfFlight/include/TOFEstimators.h b/TimeOfFlight/include/TOFEstimators.h index c5e50847..7689a0ad 100644 --- a/TimeOfFlight/include/TOFEstimators.h +++ b/TimeOfFlight/include/TOFEstimators.h @@ -2,7 +2,7 @@ #define TOFEstimators_h 1 /** -Marlin processor than calculates output parameters. +Marlin processor that calculates output parameters. \author F. Gaede, DESY, April 2018 \author B. Dudar, DESY, September 2021 */ diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index 9f21b2cd..d4318751 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -130,6 +130,8 @@ namespace TOFUtils{ Returns time measured by the closest ECal hit to the extrapolated track position at the ECal surface. \f$ \mathrm{TOF} = t_{\mathrm{closest}} - \frac{\left| \vec{r}_{\mathrm{track}} - \vec{r}_{\mathrm{closest}} \right|}{c} \f$ + + If no ECal hits are found returns `0.0`. */ double getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, double timeResolution); @@ -138,6 +140,8 @@ namespace TOFUtils{ the extrapolated track position at the ECal surface assuming speed of flight is the speed of light. \f$ \mathrm{TOF} = \frac{1}{\mathrm{MaxEcalLayer}}\sum_{i}^{\mathrm{MaxEcalLayer}} \left( t_{i} - \frac{\left|\vec{r}_{\mathrm{track}} - \vec{r}_{i} \right|}{c} \right) \f$ + + If no ECal hits are found within selectedHits, then returns `0.0`. */ double getTofFrankAvg( std::vector selectedHits, EVENT::Track* track, double timeResolution); @@ -147,6 +151,9 @@ namespace TOFUtils{ function of the distance to the extrapolated track position at the ECal surface. \f$ \mathrm{TOF}=f(|\vec{r}_{\mathrm{track}} - \vec{r}_{\mathrm{hit}} |=0) \f$ + + If no ECal hits are found within selectedHits, then returns `0.0`. + If *only one* ECal hit is found, which is not enough to perform the linear fit, then returns the same as getTofClosest() and getTofFrankAvg(). */ double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); diff --git a/TimeOfFlight/src/TOFUtils.cc b/TimeOfFlight/src/TOFUtils.cc index e2cff788..f3b3e01a 100644 --- a/TimeOfFlight/src/TOFUtils.cc +++ b/TimeOfFlight/src/TOFUtils.cc @@ -254,6 +254,7 @@ double TOFUtils::getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, do const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); + double hitTime = numeric_limits::max(); double closestDistance = numeric_limits::max(); for( auto hit : cluster->getCalorimeterHits() ){ @@ -268,6 +269,7 @@ double TOFUtils::getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, do hitTime = hit->getTime(); } } + if ( hitTime == numeric_limits::max() ) return 0.; return RandGauss::shoot(hitTime, timeResolution) - closestDistance/CLHEP::c_light; } From 2a8af491f07de8b5740e164aa8b4d2e54d4f7711 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Sat, 25 Sep 2021 00:06:09 +0200 Subject: [PATCH 07/22] Added performance testing function and fixed mem leak --- CMakeLists.txt | 4 +++- TimeOfFlight/include/TOFUtils.h | 2 ++ TimeOfFlight/src/TOFEstimators.cc | 9 ++++++- TimeOfFlight/src/TOFUtils.cc | 40 +++++++++++++++++++++++++++++++ TimeOfFlight/xml/steer.xml | 6 ++--- 5 files changed, 56 insertions(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0e1c5d2c..c63712e4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,7 +46,9 @@ include_directories( SYSTEM ${Boost_INCLUDE_DIRS} ) FIND_PACKAGE( Marlin 1.0 REQUIRED ) FIND_PACKAGE( MarlinUtil 1.4 REQUIRED ) FIND_PACKAGE( MarlinKinfit 0.0 REQUIRED ) -FIND_PACKAGE( MarlinTrk ) +# FIND_PACKAGE( MarlinTrk ) +find_package (MarlinTrk PATHS /afs/desy.de/user/d/dudarboh/iLCSoft/MarlinTrk/ NO_DEFAULT_PATH) + FIND_PACKAGE( GSL REQUIRED ) FIND_PACKAGE( ROOT 5.27 REQUIRED COMPONENTS MathMore TMVA) FIND_PACKAGE( DD4hep COMPONENTS DDRec ) diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index d4318751..c9889550 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -158,6 +158,8 @@ namespace TOFUtils{ double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); + void debugPrint(); + } diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index 6c7c00e8..6d09e1bb 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -1,6 +1,8 @@ #include "TOFEstimators.h" #include "TOFUtils.h" +#include + #include "EVENT/LCCollection.h" #include "UTIL/PIDHandler.h" @@ -86,6 +88,7 @@ void TOFEstimators::init(){ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ + auto startTime = std::chrono::steady_clock::now(); RandGauss::setTheSeed( marlin::Global::EVENTSEEDER->getSeed(this) ); ++_nEvent; streamlog_out(MESSAGE)<<"******Event****** "<<_nEvent< results{float(harmonicMom), float(trackLength), float(timeOfFlight)}; pidHandler.setParticleID(pfo , 0, 0, 0., algoID, results); + } + auto timeNow = std::chrono::steady_clock::now(); + std::chrono::duration duration = timeNow - startTime; + streamlog_out(MESSAGE)<<"Time spent (sec): "< #include +#include "marlin/VerbosityLevels.h" #include "marlinutil/CalorimeterHitType.h" #include "HelixClass.h" #include "DD4hep/Detector.h" @@ -242,6 +243,7 @@ std::vector TOFUtils::getTrackStatesPerHit(std::vector (refittedTrack.getTrackState(TrackState::AtLastHit)) ) ); if (extrapolateToEcal) trackStatesPerHit.push_back( *(static_cast (refittedTrack.getTrackState(TrackState::AtCalorimeter) ) ) ); } + delete marlinTrk; } // one can maybe use hits of refittedTrack, which includes hits that failed in the fit // code would look cleaner, but using hits failed in fit probably would have worse performance.. @@ -322,3 +324,41 @@ double TOFUtils::getTofFrankFit( std::vector selectedHit gr.Fit("pol1", "Q"); return gr.GetFunction("pol1")->GetParameter(0); } + + +void TOFUtils::debugPrint(){ + + auto parseLine = [](char* line){ + // This assumes that a digit will be found and the line ends in " Kb". + int i = strlen(line); + const char* p = line; + while (*p <'0' || *p > '9') p++; + line[i-3] = '\0'; + i = atoi(p); + return i; + }; + + auto getValue = [&parseLine](std::string memType="VmSize:"){ + //Note: this value is in KB! + int memTypeIdx; + if(memType == "VmSize:") memTypeIdx = 7; + if(memType == "VmRSS:") memTypeIdx = 6; + + FILE* file = fopen("/proc/self/status", "r"); + int result = -1; + char line[128]; + + while (fgets(line, 128, file) != NULL){ + if (strncmp(line, memType.c_str(), memTypeIdx) == 0){ + result = parseLine(line); + break; + } + } + fclose(file); + return result; + }; + + streamlog_out(MESSAGE)<<"*************Performcance info*************"< - - + @@ -21,7 +21,7 @@ - 2 + 400 0 false false From bd097be9e720bf36d745b3efb48b78b05967242e Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Sat, 25 Sep 2021 23:03:05 +0200 Subject: [PATCH 08/22] Fix mem leak with unique_ptr --- TimeOfFlight/src/TOFUtils.cc | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/TimeOfFlight/src/TOFUtils.cc b/TimeOfFlight/src/TOFUtils.cc index 363fa18f..a87a2623 100644 --- a/TimeOfFlight/src/TOFUtils.cc +++ b/TimeOfFlight/src/TOFUtils.cc @@ -189,13 +189,13 @@ std::vector TOFUtils::getTrackStatesPerHit(std::vectorcreateTrack(); + std::unique_ptr marlinTrk( trkSystem->createTrack() ); TrackImpl refittedTrack; //Need to initialize trackState at last hit TrackStateImpl preFit = *track->getTrackState(TrackState::AtLastHit); preFit.setCovMatrix( covMatrix ); - int errorFit = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, hits, &refittedTrack, IMarlinTrack::backward, &preFit, bField, maxChi2PerHit); + int errorFit = MarlinTrk::createFinalisedLCIOTrack(marlinTrk.get(), hits, &refittedTrack, IMarlinTrack::backward, &preFit, bField, maxChi2PerHit); if (errorFit != 0) continue; vector< pair > hitsInFit; @@ -208,7 +208,7 @@ std::vector TOFUtils::getTrackStatesPerHit(std::vector=0; --j ){ - TrackStateImpl ts = getTrackStateAtHit(marlinTrk, hitsInFit[j].first); + TrackStateImpl ts = getTrackStateAtHit(marlinTrk.get(), hitsInFit[j].first); trackStatesPerHit.push_back(ts); } } @@ -222,14 +222,14 @@ std::vector TOFUtils::getTrackStatesPerHit(std::vector=0; --j ){ //iterate in increasing rho - TrackStateImpl ts = getTrackStateAtHit(marlinTrk, hitsInFit[j].first); + TrackStateImpl ts = getTrackStateAtHit(marlinTrk.get(), hitsInFit[j].first); trackStatesPerHit.push_back(ts); } } else{ for( int j=0; j TOFUtils::getTrackStatesPerHit(std::vector (refittedTrack.getTrackState(TrackState::AtLastHit)) ) ); if (extrapolateToEcal) trackStatesPerHit.push_back( *(static_cast (refittedTrack.getTrackState(TrackState::AtCalorimeter) ) ) ); } - delete marlinTrk; } // one can maybe use hits of refittedTrack, which includes hits that failed in the fit // code would look cleaner, but using hits failed in fit probably would have worse performance.. From d7d00920d4a9fcb676e7c9f21e1d5df0818e8035 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Sat, 25 Sep 2021 23:12:06 +0200 Subject: [PATCH 09/22] rm detele after unique_ptr --- TimeOfFlight/xml/steer.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TimeOfFlight/xml/steer.xml b/TimeOfFlight/xml/steer.xml index a1f75e48..9fdb5a87 100644 --- a/TimeOfFlight/xml/steer.xml +++ b/TimeOfFlight/xml/steer.xml @@ -21,7 +21,7 @@ - 400 + 0 0 false false From 027b418bb6c4f256ced0b29a618b1af719c8d599 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Sat, 25 Sep 2021 23:32:23 +0200 Subject: [PATCH 10/22] restore cmake --- CMakeLists.txt | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index c63712e4..0e1c5d2c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,9 +46,7 @@ include_directories( SYSTEM ${Boost_INCLUDE_DIRS} ) FIND_PACKAGE( Marlin 1.0 REQUIRED ) FIND_PACKAGE( MarlinUtil 1.4 REQUIRED ) FIND_PACKAGE( MarlinKinfit 0.0 REQUIRED ) -# FIND_PACKAGE( MarlinTrk ) -find_package (MarlinTrk PATHS /afs/desy.de/user/d/dudarboh/iLCSoft/MarlinTrk/ NO_DEFAULT_PATH) - +FIND_PACKAGE( MarlinTrk ) FIND_PACKAGE( GSL REQUIRED ) FIND_PACKAGE( ROOT 5.27 REQUIRED COMPONENTS MathMore TMVA) FIND_PACKAGE( DD4hep COMPONENTS DDRec ) From 61a5741079f2d2b448a4b929974ad89f78475976 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Sun, 26 Sep 2021 00:00:05 +0200 Subject: [PATCH 11/22] Add debugging info --- TimeOfFlight/include/TOFUtils.h | 4 +++- TimeOfFlight/src/TOFEstimators.cc | 13 +++++++++---- TimeOfFlight/src/TOFUtils.cc | 5 ++--- TimeOfFlight/xml/steer.xml | 3 +++ 4 files changed, 17 insertions(+), 8 deletions(-) diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index c9889550..e6129cda 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -157,7 +157,9 @@ namespace TOFUtils{ */ double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); - + /** Print debug info. + Prints a time spent per event, current usage of the virtual memory (VM) and resident set size (RSS). + */ void debugPrint(); } diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index 6d09e1bb..2e6a3f59 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -89,9 +89,10 @@ void TOFEstimators::init(){ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ auto startTime = std::chrono::steady_clock::now(); + RandGauss::setTheSeed( marlin::Global::EVENTSEEDER->getSeed(this) ); ++_nEvent; - streamlog_out(MESSAGE)<<"******Event****** "<<_nEvent<getCollection(_pfoCollectionName); @@ -185,10 +186,14 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ } vector results{float(harmonicMom), float(trackLength), float(timeOfFlight)}; pidHandler.setParticleID(pfo , 0, 0, 0., algoID, results); + streamlog_out(DEBUG6)<<"*****PFO***** "<< i+1< duration = timeNow - startTime; - streamlog_out(MESSAGE)<<"Time spent (sec): "< duration = endTime - startTime; + streamlog_out(DEBUG7)<<"Time spent (sec): "< 10 + + DEBUG6 + From acb229acd1f00cacd191a5db44b650454f4c9287 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Sun, 26 Sep 2021 01:35:53 +0200 Subject: [PATCH 12/22] Add exception handeling and warnings. --- TimeOfFlight/src/TOFEstimators.cc | 56 +++++++++++++++++++++++-------- TimeOfFlight/src/TOFUtils.cc | 22 ++++++------ TimeOfFlight/xml/steer.xml | 16 +++++++-- 3 files changed, 67 insertions(+), 27 deletions(-) diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index 2e6a3f59..7e678d31 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -69,10 +69,16 @@ TOFEstimators::TOFEstimators() : marlin::Processor("TOFEstimators") { "Time of flight is calculated using Ecal hits only up to MaxLayer", _maxEcalLayer, int(10) ); + } void TOFEstimators::init(){ + + if(_tofMethod != "closest" && _tofMethod != "frankAvg" && _tofMethod != "frankFit"){ + throw EVENT::Exception( "Invalid steering parameter for TofMethod is passed: " + _tofMethod + "\n Available options are: closest, frankAvg, frankFit" ); + } + marlin::Global::EVENTSEEDER->registerProcessor(this); _outputParNames = {"momentumHM", "trackLength", "timeOfFlight"}; @@ -88,15 +94,30 @@ void TOFEstimators::init(){ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ - auto startTime = std::chrono::steady_clock::now(); - RandGauss::setTheSeed( marlin::Global::EVENTSEEDER->getSeed(this) ); ++_nEvent; - streamlog_out(DEBUG7)<<"************Event************ "<<_nEvent<getCollection(_pfoCollectionName); + } + catch(lcio::DataNotAvailableException& e){ + streamlog_out(WARNING)<<"Input collection: "<<_pfoCollectionName<<" is not found - skipping event "<<_nEvent<getCollection(_pfoCollectionName); - LCRelationNavigator navigatorSET = LCRelationNavigator( evt->getCollection("SETSpacePointRelations") ); + LCCollection* setRelations = nullptr; + try{ + setRelations = evt->getCollection("SETSpacePointRelations"); + } + catch(lcio::DataNotAvailableException& e){ + streamlog_out(WARNING)<<"Collection: SETSpacePointRelations is not found - skipping event "<<_nEvent< 0.5 ){ @@ -146,13 +168,11 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ } harmonicMom = std::sqrt(trackLength/harmonicMom); - ////////////////////////////////////// // This part calculates Time of flight ////////////////////////////////////// double timeOfFlight = 0.; if( _extrapolateToEcal ){ - // if we extrapolate to the calorimeter check TOF method steering parameter if (_tofMethod == "closest"){ timeOfFlight = getTofClosest(cluster, track, _timeResolution/1000.); } @@ -164,24 +184,32 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ vector frankHits = selectFrankEcalHits(cluster, track, _maxEcalLayer, _bField); timeOfFlight = getTofFrankFit(frankHits, track, _timeResolution/1000.); } - else{ - //wrong parameter silly! - } } else{ //define tof as average time between two SET strips - //if no SET hits found, tof=0 + //if no SET hits found, tof alreasy is 0, just skip TrackerHit* hitSET = getSETHit(track, _tpcOuterR); if ( hitSET != nullptr ){ const vector& simHitsSET = navigatorSET.getRelatedToObjects( hitSET ); - if ( simHitsSET.size() == 2 ){ - //it should be 2 always.. but just to be safe + if ( simHitsSET.size() >= 2 ){ + //It must be always 2, but just in case... + if (simHitsSET.size() > 2) streamlog_out(WARNING)<<"Found more than two SET strip hits! Writing TOF as an average of the first two elements in the array."<( simHitsSET[0] ); SimTrackerHit* simHitSETBack = static_cast ( simHitsSET[1] ); double timeFront = RandGauss::shoot( simHitSETFront->getTime(), _timeResolution/1000. ); double timeBack = RandGauss::shoot( simHitSETBack->getTime(), _timeResolution/1000. ); timeOfFlight = (timeFront + timeBack)/2.; } + else if (simHitsSET.size() == 1){ + streamlog_out(WARNING)<<"Found only one SET strip hit! Writing TOF from a single strip."<( simHitsSET[0] ); + timeOfFlight = RandGauss::shoot( simHitSET->getTime(), _timeResolution/1000. ); + } + else{ + // this probably never happens... + streamlog_out(WARNING)<<"Found NO simHits associated with the found SET hit! Writing TOF as 0."< results{float(harmonicMom), float(trackLength), float(timeOfFlight)}; @@ -190,8 +218,8 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ streamlog_out(DEBUG6)<<"momentum: "<< float(harmonicMom)<<" Gev"< duration = endTime - startTime; streamlog_out(DEBUG7)<<"Time spent (sec): "< hits = track->getTrackerHits(); TrackerHit* lastHit = hits.back(); - Vector3D pos (lastHit->getPosition()); + if ( pos.rho() > tpcOuterR ) return lastHit; return nullptr; } @@ -128,7 +128,7 @@ std::vector TOFUtils::selectFrankEcalHits( EVENT::Cluste for (int l=0; l::max(); for ( auto hit : cluster->getCalorimeterHits() ){ @@ -160,11 +160,14 @@ std::vector TOFUtils::getSubTracks(EVENT::Track* track){ int nSubTrack0Hits = track->getTracks()[0]->getTrackerHits().size(); int nSubTrack1Hits = track->getTracks()[1]->getTrackerHits().size(); + //OPTIMIZE: this is not reliable, but seems no other way is possible. + //Read documentation in the header file for details. int startIdx; if( std::abs(nTPCHits - nSubTrack0Hits) <= 1 ) startIdx = 1; else if ( std::abs(nTPCHits - nSubTrack1Hits) <= 1 ) startIdx = 2; else{ // this shouldn't happen in princinple at all... + streamlog_out(WARNING)<<"Can't understand which subTrack is responsible for the first TPC hits! Skip adding subTracks."<getTracks()[j] ); @@ -178,7 +181,7 @@ std::vector TOFUtils::getTrackStatesPerHit(std::vector hits = track->getTrackerHits(); - sort(hits.begin(), hits.end(), sortByRho); + std::sort(hits.begin(), hits.end(), sortByRho); // setup initial dummy covariance matrix vector covMatrix(15); @@ -239,14 +242,14 @@ std::vector TOFUtils::getTrackStatesPerHit(std::vector (refittedTrack.getTrackState(TrackState::AtLastHit)) ) ); if (extrapolateToEcal) trackStatesPerHit.push_back( *(static_cast (refittedTrack.getTrackState(TrackState::AtCalorimeter) ) ) ); } } - // one can maybe use hits of refittedTrack, which includes hits that failed in the fit - // code would look cleaner, but using hits failed in fit probably would have worse performance.. - // needs to be checked + // one can maybe use hits of refittedTrack, but they include also hits that had failed in the fit + // code would look cleaner, but using hits that are failed in fit probably would have worse performance.. + // needs to be checked. return trackStatesPerHit; } @@ -255,7 +258,6 @@ double TOFUtils::getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, do const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); - double hitTime = numeric_limits::max(); double closestDistance = numeric_limits::max(); for( auto hit : cluster->getCalorimeterHits() ){ @@ -270,8 +272,8 @@ double TOFUtils::getTofClosest( EVENT::Cluster* cluster, EVENT::Track* track, do hitTime = hit->getTime(); } } - if ( hitTime == numeric_limits::max() ) return 0.; + if ( hitTime == numeric_limits::max() ) return 0.; return RandGauss::shoot(hitTime, timeResolution) - closestDistance/CLHEP::c_light; } @@ -315,7 +317,7 @@ double TOFUtils::getTofFrankFit( std::vector selectedHit y.push_back(time); xErr.push_back(0.); yErr.push_back(0.3); - //setting this to 0 is not good for the fit.. So I put random 300ps... + //OPTIMIZE: setting this to 0 is not good for the fit.. So I put random 300ps... //Changing this doesn't seem to affect results, although one may want to check it more carefully } diff --git a/TimeOfFlight/xml/steer.xml b/TimeOfFlight/xml/steer.xml index a73cae37..bfae8738 100644 --- a/TimeOfFlight/xml/steer.xml +++ b/TimeOfFlight/xml/steer.xml @@ -8,8 +8,9 @@ - + + + @@ -53,7 +54,7 @@ 10 - DEBUG6 + MESSAGE @@ -78,6 +79,15 @@ + + + frankFit + + + + 100 + + From 6a6bc59daf44a15ef46fc6e44785974d0d4670e5 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Wed, 6 Oct 2021 11:54:18 +0200 Subject: [PATCH 13/22] Resolve Thomas comments --- TimeOfFlight/include/TOFEstimators.h | 12 +++--- TimeOfFlight/include/TOFUtils.h | 2 +- TimeOfFlight/src/TOFEstimators.cc | 56 +++++++++------------------- TimeOfFlight/src/TOFUtils.cc | 45 ++++++++++------------ 4 files changed, 44 insertions(+), 71 deletions(-) diff --git a/TimeOfFlight/include/TOFEstimators.h b/TimeOfFlight/include/TOFEstimators.h index 7689a0ad..9750ef10 100644 --- a/TimeOfFlight/include/TOFEstimators.h +++ b/TimeOfFlight/include/TOFEstimators.h @@ -1,17 +1,17 @@ #ifndef TOFEstimators_h #define TOFEstimators_h 1 -/** -Marlin processor that calculates output parameters. -\author F. Gaede, DESY, April 2018 -\author B. Dudar, DESY, September 2021 -*/ - #include #include #include "marlin/Processor.h" #include "MarlinTrk/IMarlinTrkSystem.h" + +/** +Marlin processor that calculates output parameters. +\author F. Gaede, DESY, April 2018 +\author B. Dudar, DESY, September 2021 +*/ class TOFEstimators : public marlin::Processor { public: /** diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index e6129cda..a1fab649 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -158,7 +158,7 @@ namespace TOFUtils{ double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); /** Print debug info. - Prints a time spent per event, current usage of the virtual memory (VM) and resident set size (RSS). + Prints a current usage of the virtual memory (VM) and resident set size (RSS). */ void debugPrint(); diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index 7e678d31..afcf21b1 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -84,6 +84,8 @@ void TOFEstimators::init(){ _outputParNames = {"momentumHM", "trackLength", "timeOfFlight"}; _bField = MarlinUtil::getBzAtOrigin(); _tpcOuterR = getTPCOuterR(); + // internally we use time resolution in nanoseconds + _timeResolution = _timeResolution/1000.; _trkSystem = MarlinTrk::Factory::createMarlinTrkSystem("DDKalTest", nullptr, ""); _trkSystem->setOption( MarlinTrk::IMarlinTrkSystem::CFG::useQMS, true); @@ -99,23 +101,8 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ streamlog_out(MESSAGE)<<"************Event************ "<<_nEvent<getCollection(_pfoCollectionName); - } - catch(lcio::DataNotAvailableException& e){ - streamlog_out(WARNING)<<"Input collection: "<<_pfoCollectionName<<" is not found - skipping event "<<_nEvent<getCollection("SETSpacePointRelations"); - } - catch(lcio::DataNotAvailableException& e){ - streamlog_out(WARNING)<<"Collection: SETSpacePointRelations is not found - skipping event "<<_nEvent<getCollection(_pfoCollectionName); + LCCollection* setRelations = evt->getCollection("SETSpacePointRelations"); LCRelationNavigator navigatorSET = LCRelationNavigator( setRelations ); @@ -148,20 +135,13 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ double harmonicMom = 0.; int nTrackStates = trackStates.size(); for( int j=1; j < nTrackStates; ++j ){ - //for the last hit we check which track length formula to use - if (j == nTrackStates - 1){ - double nTurns = getHelixNRevolutions( trackStates[j-1], trackStates[j] ); - if( nTurns > 0.5 ){ - // we cannot calculate arc length for more than pi revolution. Use formula with only z - double arcLength = getHelixLengthAlongZ( trackStates[j-1], trackStates[j] ); - Vector3D mom = getHelixMomAtTrackState( trackStates[j-1], _bField ); - trackLength += arcLength; - harmonicMom += arcLength/mom.r2(); - continue; - } - } + //we check which track length formula to use + double nTurns = getHelixNRevolutions( trackStates[j-1], trackStates[j] ); + double arcLength; + // we cannot calculate arc length for more than pi revolution. Use formula with only z + if ( nTurns <= 0.5 ) arcLength = getHelixArcLength( trackStates[j-1], trackStates[j] ); + else arcLength = getHelixLengthAlongZ( trackStates[j-1], trackStates[j] ); - double arcLength = getHelixArcLength( trackStates[j-1], trackStates[j] ); Vector3D mom = getHelixMomAtTrackState( trackStates[j-1], _bField ); trackLength += arcLength; harmonicMom += arcLength/mom.r2(); @@ -174,19 +154,19 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ double timeOfFlight = 0.; if( _extrapolateToEcal ){ if (_tofMethod == "closest"){ - timeOfFlight = getTofClosest(cluster, track, _timeResolution/1000.); + timeOfFlight = getTofClosest(cluster, track, _timeResolution); } else if (_tofMethod == "frankAvg"){ vector frankHits = selectFrankEcalHits(cluster, track, _maxEcalLayer, _bField); - timeOfFlight = getTofFrankAvg(frankHits, track, _timeResolution/1000.); + timeOfFlight = getTofFrankAvg(frankHits, track, _timeResolution); } else if (_tofMethod == "frankFit"){ vector frankHits = selectFrankEcalHits(cluster, track, _maxEcalLayer, _bField); - timeOfFlight = getTofFrankFit(frankHits, track, _timeResolution/1000.); + timeOfFlight = getTofFrankFit(frankHits, track, _timeResolution); } } else{ - //define tof as average time between two SET strips + //define tof as an average time between two SET strips //if no SET hits found, tof alreasy is 0, just skip TrackerHit* hitSET = getSETHit(track, _tpcOuterR); if ( hitSET != nullptr ){ @@ -197,14 +177,14 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ SimTrackerHit* simHitSETFront = static_cast ( simHitsSET[0] ); SimTrackerHit* simHitSETBack = static_cast ( simHitsSET[1] ); - double timeFront = RandGauss::shoot( simHitSETFront->getTime(), _timeResolution/1000. ); - double timeBack = RandGauss::shoot( simHitSETBack->getTime(), _timeResolution/1000. ); + double timeFront = RandGauss::shoot(simHitSETFront->getTime(), _timeResolution); + double timeBack = RandGauss::shoot(simHitSETBack->getTime(), _timeResolution); timeOfFlight = (timeFront + timeBack)/2.; } else if (simHitsSET.size() == 1){ streamlog_out(WARNING)<<"Found only one SET strip hit! Writing TOF from a single strip."<( simHitsSET[0] ); - timeOfFlight = RandGauss::shoot( simHitSET->getTime(), _timeResolution/1000. ); + SimTrackerHit* simHitSET = static_cast (simHitsSET[0]); + timeOfFlight = RandGauss::shoot(simHitSET->getTime(), _timeResolution); } else{ // this probably never happens... diff --git a/TimeOfFlight/src/TOFUtils.cc b/TimeOfFlight/src/TOFUtils.cc index 93fcf60d..8ffd4a00 100644 --- a/TimeOfFlight/src/TOFUtils.cc +++ b/TimeOfFlight/src/TOFUtils.cc @@ -120,31 +120,28 @@ EVENT::TrackerHit* TOFUtils::getSETHit(EVENT::Track* track, double tpcOuterR){ std::vector TOFUtils::selectFrankEcalHits( EVENT::Cluster* cluster, EVENT::Track* track, int maxEcalLayer, double bField ){ - vector selectedHits; - const TrackState* tsEcal = track->getTrackState(TrackState::AtCalorimeter); Vector3D trackPosAtEcal ( tsEcal->getReferencePoint() ); Vector3D trackMomAtEcal = TOFUtils::getHelixMomAtTrackState(*tsEcal, bField); - for (int l=0; l::max(); - for ( auto hit : cluster->getCalorimeterHits() ){ - CHT hitType( hit->getType() ); - bool isECALHit = ( hitType.caloID() == CHT::ecal ); - if ( (!isECALHit) || ( int( hitType.layer() ) != l) ) continue; - - Vector3D hitPos( hit->getPosition() ); - double dToLine = (hitPos - trackPosAtEcal).cross(trackMomAtEcal.unit()).r(); - if (dToLine < closestDistanceToLine){ - closestDistanceToLine = dToLine; - selectedHit = hit; - } + vector selectedHits(maxEcalLayer, nullptr); + vector minDistances(maxEcalLayer, numeric_limits::max()); + + for ( auto hit : cluster->getCalorimeterHits() ){ + CHT hitType( hit->getType() ); + bool isECALHit = ( hitType.caloID() == CHT::ecal ); + int layer = hitType.layer(); + if ( (!isECALHit) || ( layer >= maxEcalLayer) ) continue; + + Vector3D hitPos( hit->getPosition() ); + double dToLine = (hitPos - trackPosAtEcal).cross(trackMomAtEcal.unit()).r(); + if ( dToLine < minDistances[layer] ){ + minDistances[layer] = dToLine; + selectedHits[layer] = hit; } - if ( selectedHit != nullptr ) selectedHits.push_back(selectedHit); } + selectedHits.erase( std::remove_if( selectedHits.begin(), selectedHits.end(), [](CalorimeterHit* h) { return h == nullptr; } ), selectedHits.end() ); + return selectedHits; } @@ -167,7 +164,7 @@ std::vector TOFUtils::getSubTracks(EVENT::Track* track){ else if ( std::abs(nTPCHits - nSubTrack1Hits) <= 1 ) startIdx = 2; else{ // this shouldn't happen in princinple at all... - streamlog_out(WARNING)<<"Can't understand which subTrack is responsible for the first TPC hits! Skip adding subTracks."<getTracks()[j] ); @@ -308,20 +305,16 @@ double TOFUtils::getTofFrankFit( std::vector selectedHit return RandGauss::shoot(selectedHits[0]->getTime(), timeResolution) - dToTrack/CLHEP::c_light; } - vector x, y, xErr, yErr; + vector x, y; for ( auto hit : selectedHits ){ Vector3D hitPos( hit->getPosition() ); double dToTrack = (hitPos - trackPosAtEcal).r(); x.push_back(dToTrack); double time = RandGauss::shoot(hit->getTime(), timeResolution); y.push_back(time); - xErr.push_back(0.); - yErr.push_back(0.3); - //OPTIMIZE: setting this to 0 is not good for the fit.. So I put random 300ps... - //Changing this doesn't seem to affect results, although one may want to check it more carefully } - TGraphErrors gr(nHits, &x[0], &y[0], &xErr[0], &yErr[0]); + TGraph gr(nHits, x.data(), y.data()); gr.Fit("pol1", "Q"); return gr.GetFunction("pol1")->GetParameter(0); } From 49690d8b5fa8118f8aa76c83abf9616c048df40f Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Wed, 6 Oct 2021 13:34:05 +0200 Subject: [PATCH 14/22] minor comment changes --- TimeOfFlight/include/TOFUtils.h | 13 ++++++------- TimeOfFlight/src/TOFEstimators.cc | 2 +- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index a1fab649..0ce46ff5 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -1,13 +1,6 @@ #ifndef TOFUtils_h #define TOFUtils_h 1 -/** - * Utility functions that are used by the TOFEstimators processor. - * - * \author F. Gaede, DESY, 2018 - * \author B. Dudar, DESY, 2021 - */ - #include #include "EVENT/ReconstructedParticle.h" #include "IMPL/TrackStateImpl.h" @@ -15,6 +8,12 @@ #include "MarlinTrk/IMarlinTrack.h" #include "DDRec/Vector3D.h" +/** + * Utility functions that are used by the TOFEstimators processor. + * + * \author F. Gaede, DESY, 2018 + * \author B. Dudar, DESY, 2021 +*/ namespace TOFUtils{ /** Comparator function by radius for tracker hits. diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index afcf21b1..ce99ef92 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -138,7 +138,7 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ //we check which track length formula to use double nTurns = getHelixNRevolutions( trackStates[j-1], trackStates[j] ); double arcLength; - // we cannot calculate arc length for more than pi revolution. Use formula with only z + // we cannot calculate arc length for more than pi revolution using delta phi. Use formula with only z if ( nTurns <= 0.5 ) arcLength = getHelixArcLength( trackStates[j-1], trackStates[j] ); else arcLength = getHelixLengthAlongZ( trackStates[j-1], trackStates[j] ); From 40d6ea06587cbaef83e5f79582a80f37d3c46590 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Wed, 6 Oct 2021 18:49:37 +0200 Subject: [PATCH 15/22] Make an example agree with the README --- TimeOfFlight/xml/steer.xml | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/TimeOfFlight/xml/steer.xml b/TimeOfFlight/xml/steer.xml index bfae8738..ec8d7517 100644 --- a/TimeOfFlight/xml/steer.xml +++ b/TimeOfFlight/xml/steer.xml @@ -10,7 +10,6 @@ - @@ -79,16 +78,6 @@ - - - frankFit - - - - 100 - - - From ab9970f12632f71338153bba218fed6a39b6493e Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Mon, 11 Oct 2021 09:39:42 +0200 Subject: [PATCH 16/22] Finilize Thomas comments --- TimeOfFlight/include/TOFUtils.h | 4 --- TimeOfFlight/src/TOFEstimators.cc | 2 +- TimeOfFlight/src/TOFUtils.cc | 43 +++---------------------------- 3 files changed, 4 insertions(+), 45 deletions(-) diff --git a/TimeOfFlight/include/TOFUtils.h b/TimeOfFlight/include/TOFUtils.h index 0ce46ff5..bc07372b 100644 --- a/TimeOfFlight/include/TOFUtils.h +++ b/TimeOfFlight/include/TOFUtils.h @@ -156,10 +156,6 @@ namespace TOFUtils{ */ double getTofFrankFit( std::vector selectedHits, EVENT::Track* track, double timeResolution); - /** Print debug info. - Prints a current usage of the virtual memory (VM) and resident set size (RSS). - */ - void debugPrint(); } diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index ce99ef92..ea5ced78 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -187,7 +187,7 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ timeOfFlight = RandGauss::shoot(simHitSET->getTime(), _timeResolution); } else{ - // this probably never happens... + // this happens very rarily (0.1%). When >1 simHits associated with a single strip none simHits are written by the DDSpacePointBuilder. streamlog_out(WARNING)<<"Found NO simHits associated with the found SET hit! Writing TOF as 0."< TOFUtils::getSubTracks(EVENT::Track* track){ int nSubTrack0Hits = track->getTracks()[0]->getTrackerHits().size(); int nSubTrack1Hits = track->getTracks()[1]->getTrackerHits().size(); - //OPTIMIZE: this is not reliable, but seems no other way is possible. + //OPTIMIZE: this is not reliable, but I don't see any other way at the moment. //Read documentation in the header file for details. int startIdx; if( std::abs(nTPCHits - nSubTrack0Hits) <= 1 ) startIdx = 1; else if ( std::abs(nTPCHits - nSubTrack1Hits) <= 1 ) startIdx = 2; else{ - // this shouldn't happen in princinple at all... - streamlog_out(ERROR)<<"Can't understand which subTrack is responsible for the first TPC hits! Skip adding subTracks."<getTracks()[j] ); @@ -318,40 +318,3 @@ double TOFUtils::getTofFrankFit( std::vector selectedHit gr.Fit("pol1", "Q"); return gr.GetFunction("pol1")->GetParameter(0); } - - -void TOFUtils::debugPrint(){ - - auto parseLine = [](char* line){ - // This assumes that a digit will be found and the line ends in " Kb". - int i = strlen(line); - const char* p = line; - while (*p <'0' || *p > '9') p++; - line[i-3] = '\0'; - i = atoi(p); - return i; - }; - - auto getValue = [&parseLine](std::string memType="VmSize:"){ - //Note: this value is in KB! - int memTypeIdx; - if(memType == "VmSize:") memTypeIdx = 7; - if(memType == "VmRSS:") memTypeIdx = 6; - - FILE* file = fopen("/proc/self/status", "r"); - int result = -1; - char line[128]; - - while (fgets(line, 128, file) != NULL){ - if (strncmp(line, memType.c_str(), memTypeIdx) == 0){ - result = parseLine(line); - break; - } - } - fclose(file); - return result; - }; - - streamlog_out(DEBUG7)<<"Using VM (MB): "< Date: Mon, 11 Oct 2021 10:17:54 +0200 Subject: [PATCH 17/22] Finilize readme, remove debugpring --- TimeOfFlight/README.md | 63 ++++++++++++++++--------------- TimeOfFlight/src/TOFEstimators.cc | 1 - 2 files changed, 32 insertions(+), 32 deletions(-) diff --git a/TimeOfFlight/README.md b/TimeOfFlight/README.md index fdcccd64..e210442e 100644 --- a/TimeOfFlight/README.md +++ b/TimeOfFlight/README.md @@ -1,16 +1,17 @@ -# Description +## Description TOFEstimators is the Marlin processor that calculates track length and time-of-flight for charged particles that can be used for the mass reconstruction and particle identification.
In addition it calculates square root of the harmonic mean of the momentum squared. One can use it as a better substitute for the particle momentum in the mass reconstruction. -[Doxygen documentation](https://www.desy.de/~dudarboh/tof_doc/html/index.html) contains a detailed description of the source code.
-In case reader uses dark theme on the github and finds hard to read the formulas he can find this README as the main page there. +The detailed description of the source code is available in the [doxygen documentation](https://www.desy.de/~dudarboh/tof_doc/html/index.html).
-# Requirements +If formulas are barely visible please switch to the light theme of the github or check the main page in the [doxygen documentation](https://www.desy.de/~dudarboh/tof_doc/html/index.html). + +## Requirements TOFEstimators is currently working with only simple cases of particle flow objects.
-The analyzed *ReconstructedParticle* object must have **exactly one** associated `Track` and **exactly one** associated `Cluster`.
-In other cases it still produces PIDHandler, but with `0.0` for every output parameter. +The analyzed `ReconstructedParticle` object must have **exactly one** associated `Track` and **exactly one** associated `Cluster`.
+In other cases it still produces PIDHandler, but with `0.0` as output for every output parameter. The input sclio file must contain a lot of various tracker and calorimeter hit collections.
Thus, input slcio file can be e.g. *REC* format, but **cannot** be *DST* or *mini-DST* format. @@ -18,85 +19,85 @@ Thus, input slcio file can be e.g. *REC* format, but **cannot** be *DST* or *min TOFEstimators is dependent on InitDD4hep processor to extract detector geometry details.
Make sure to run InitDD4hep before this processor is executed. -# Steering parameters +## Steering parameters -There are five steering parameters: +This processor has five steering parameters: -### ReconstructedParticleCollection / default: `"PandoraPFOs"` ++ "ReconstructedParticleCollection" | default: "PandoraPFOs"
Collection of `ReconstructedParticle` objects to analyze.
The final results will be written in the PIDHandler inside this collection. -### ExtrapolateToEcal / default: `true`
++ "ExtrapolateToEcal" | default: true
If `true` the track length is calculated to the extrapolated track position at the ECal surface. Time-of-flight then calculated using ECal hits via one of the methods chosen in *TofMethod* steering parameter. If `false` the track length is calculated to the last tracker hit. Time-of-flight then calculated using the last tracker hit. Which is SET hit for the barrel region and *undefined* for the endcap region.
If the last tracker hit is the SET hit time-of-flight is set to the average time of two SET strips.
If the last tracker hit is not the SET hit time-of-flight is set to `0`. -### TofMethod / options: `"closest"`, `"frankAvg"`, `"frankFit"` / default: `"closest"` - ++ TofMethod | options: "closest", "frankAvg", "frankFit" | default: "closest"
If *ExtrapolateToEcal* is set to `false` then this steering parameter is ignored. If *ExtrapolateToEcal* is set to `true` then it defines how to use ECal hits to calculate the time-of-flight at the ECal surface. - - `"closest"` uses the closest ECal hit to the extrapolated track position at the ECal surface.
+ - "closest" uses the closest ECal hit to the extrapolated track position at the ECal surface.
Time-of-flight is defined as:
If no ECal hits are found returns `0.0`. - - `"frankAvg"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Define time-of-flight as an average of their time corrected for the distance to the track position at the ECal surface assuming speed of flight is the speed of light. + - "frankAvg" Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Define time-of-flight as an average of their time corrected for the distance to the track position at the ECal surface assuming speed of flight is the speed of light.
If no ECal hits are found returns `0.0`. - - - `"frankFit"` Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Use linear fit of the hits time as a function of the distance to the track position at the ECal surface to define time-of-flight as an extrapolated time at the extrapolated track position at the ECal surface
+ - "frankFit" Select the closest ECal hit to the linearly extrapolated track line inside the ECal in each of the first *MaxEcalLayer* ECal layers. Use linear fit of the hits time as a function of the distance to the track position at the ECal surface to define time-of-flight as an extrapolated time at the extrapolated track position at the ECal surface
If no ECal hits are found returns `0.0`.
If only one ECal hit is found, which is not enough to perform the linear fit, then returns the same as *closest* and *frankAvg*. for more illustrative explanations of the methods above see slides 10, 12, 13 from the following [LCWS2021 talk]((https://indico.cern.ch/event/995633/contributions/4259659/attachments/2209010/3738157/Bohdan_TOF_LCWS2021.pdf)). -### TimeResolution / default: `0.0` (ps) ++ TimeResolution | default: 0.0 (ps)
If *ExtrapolateToEcal* is set to `true` it defines times resolution of individual ECal hits in ps. If *ExtrapolateToEcal* is set to `false` it defines times resolution of individual SET strips in ps. -### MaxEcalLayer / default: `10` ++ MaxEcalLayer | default: 10
Defines number of inner ECal layers to use for the *frankAvg* and *frankFit* algorithms. If *frankAvg* and *frankFit* are unused then this steering parameter is ignored. -# Output parameters +## Output parameters -There are three output parameters: - -### `momentumHM` (GeV/c) +This processor has three output parameters: ++ "momentumHM" (GeV/c)
A squared root of the harmonic mean of the squared momentum: If momentum is constant, then it is just a momentum:
If momentum changes, then using instead of or is mathematically rigorous choice for the relativistic particle assumption
See [Winfried A. Mitaroff paper](https://arxiv.org/abs/2107.02031) for the details. -### `trackLength` (mm) ++ "trackLength" (mm)
A track length calculated between the point of the closest approach (PCA) to the IP `(0,0,0)`
and the ECal surface if *ExtrapolateToEcal* is `true`
or the last tracker hit if *ExtrapolateToEcal* is `false`. A total track length is calculated as a sum of track segments obtained from iterating over track states defined at the IP, every tracker hit and the ECal surface if *ExtrapolateToEcal* is `true`.
- A track segment length is calculated between neighboring track states as: - - If *ExtrapolateToEcal* is `true` then we additionally check for the number of helix turns between the last tracker hit and the extrapolated track position at the ECal surface. If we are unable to use the formula above to calculate the length of the last segment. Thus we use different formula that does not rely on the azimuthal angle: + For every pair of track states we check an approximate number of revolutions of the helix with . + + In case the track segment length is calculated between neighboring track states as: + .
+ In case the track segment length is calculated between neighboring track states as: + .
- + If helix has more than half revolution between two tracker states then we are unable to use the formula with phi angle, thus we use different formula that does not rely on the azimuthal angle but only on z coordinate and the dip of the helix. Later formula is less precise, so we use it only in these exception cases. More than half turn situation generally should not happen between TPC hits as neighboring TPC hits usually close to each other. But more than half turn often can happen between last TPC hit and the extrapolated track state at the ECal endcap. -### `timeOfFlight` (ns) ++ "timeOfFlight" (ns)
Time-of-flight of the particle assuming . For the detailed description of different time-of-flight estimators see the **TofMethod** bullet point in the Steering parameters section. -# Steering file example +## Steering file example [./xml/steer.xml](./xml/steer.xml) is a steering file example that runs three TOFEstimators processors: *MyTofClosest0ps*, *MyTofSET10ps* and *MyTofFrankAvg50ps* and then writes their output parameters in the PIDHandlers of the *PandoraPFOs* collection in the new *output.slcio* file using LCIOOutputProcessor. @@ -141,7 +142,7 @@ Scrolling more down, one can see calculated parameters for each individual parti | 0 | 0.0000e+00 | 000000 | 11 | [ momentumHM : 3.84e+00, trackLength : 2.85e+03, timeOfFlight : 9.57e+00,] . . . -# Extracting output parameters +## Extracting output parameters Here is the pseudocode example on how one could extract PIDHandler parameters of *MyTofClosest0ps* results from the slcio file in your own Marlin processor and then use them to e.g. reconstruct the mass: @@ -173,7 +174,7 @@ Here is the pseudocode example on how one could extract PIDHandler parameters of } -# Authors +## Authors - N.Weinhold, DESY, internship 2017
- F.Gaede, DESY, 2018
- B.Dudar, DESY, 2021
diff --git a/TimeOfFlight/src/TOFEstimators.cc b/TimeOfFlight/src/TOFEstimators.cc index ea5ced78..16223cfe 100644 --- a/TimeOfFlight/src/TOFEstimators.cc +++ b/TimeOfFlight/src/TOFEstimators.cc @@ -203,5 +203,4 @@ void TOFEstimators::processEvent(EVENT::LCEvent * evt){ auto endTime = std::chrono::steady_clock::now(); std::chrono::duration duration = endTime - startTime; streamlog_out(DEBUG7)<<"Time spent (sec): "< Date: Mon, 11 Oct 2021 10:19:19 +0200 Subject: [PATCH 18/22] minor readme style --- TimeOfFlight/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/TimeOfFlight/README.md b/TimeOfFlight/README.md index e210442e..21c59f07 100644 --- a/TimeOfFlight/README.md +++ b/TimeOfFlight/README.md @@ -23,11 +23,11 @@ Make sure to run InitDD4hep before this processor is executed. This processor has five steering parameters: -+ "ReconstructedParticleCollection" | default: "PandoraPFOs"
++ ReconstructedParticleCollection | default: "PandoraPFOs"
Collection of `ReconstructedParticle` objects to analyze.
The final results will be written in the PIDHandler inside this collection. -+ "ExtrapolateToEcal" | default: true
++ ExtrapolateToEcal | default: true
If `true` the track length is calculated to the extrapolated track position at the ECal surface. Time-of-flight then calculated using ECal hits via one of the methods chosen in *TofMethod* steering parameter. If `false` the track length is calculated to the last tracker hit. Time-of-flight then calculated using the last tracker hit. Which is SET hit for the barrel region and *undefined* for the endcap region.
From e333e308cabedca772eb5656724bfe9410e429d2 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Mon, 11 Oct 2021 15:32:52 +0200 Subject: [PATCH 19/22] readme typo --- TimeOfFlight/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/TimeOfFlight/README.md b/TimeOfFlight/README.md index 21c59f07..bdee0f3a 100644 --- a/TimeOfFlight/README.md +++ b/TimeOfFlight/README.md @@ -158,7 +158,7 @@ Here is the pseudocode example on how one could extract PIDHandler parameters of int tofParIdx = pidHandler.getParameterIndex(algorithmID, "timeOfFlight"); - for(int i=0; i < pfos.getNumberOfElements(); ++i){ + for(int i=0; i < pfos->getNumberOfElements(); ++i){ ReconstructedParticle* pfo = static_cast ( pfos->getElementAt(i) ); const ParticleID& pfoID = pidHandler.getParticleID(pfo, algorithmID); From 3d98701f521817e7969eab0ebd47774429e645c4 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Mon, 11 Oct 2021 15:36:50 +0200 Subject: [PATCH 20/22] update docstring --- TimeOfFlight/include/TOFEstimators.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/TimeOfFlight/include/TOFEstimators.h b/TimeOfFlight/include/TOFEstimators.h index 9750ef10..5fb9d4d6 100644 --- a/TimeOfFlight/include/TOFEstimators.h +++ b/TimeOfFlight/include/TOFEstimators.h @@ -6,9 +6,8 @@ #include "marlin/Processor.h" #include "MarlinTrk/IMarlinTrkSystem.h" - /** -Marlin processor that calculates output parameters. +Marlin processor that calculates harmonic mean momentum, track length and time-of-flight for charged particles. \author F. Gaede, DESY, April 2018 \author B. Dudar, DESY, September 2021 */ From f8c905985cceec7633cdc77604a92346cbadae03 Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Mon, 11 Oct 2021 15:48:00 +0200 Subject: [PATCH 21/22] Docstring comment cleanup --- TimeOfFlight/include/TOFEstimators.h | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/TimeOfFlight/include/TOFEstimators.h b/TimeOfFlight/include/TOFEstimators.h index 5fb9d4d6..1b5a82dd 100644 --- a/TimeOfFlight/include/TOFEstimators.h +++ b/TimeOfFlight/include/TOFEstimators.h @@ -14,16 +14,14 @@ Marlin processor that calculates harmonic mean momentum, track length and time-o class TOFEstimators : public marlin::Processor { public: /** - Copy constructor . + Copy constructor. We remove it to avoid W-effc++ warnings. - Copying objects with pointer members is a bad idea. */ TOFEstimators(const TOFEstimators&) = delete; /** Copy assignment operator. We remove it to avoid W-effc++ warnings. - Copying objects with pointer members is a bad idea. */ TOFEstimators& operator=(const TOFEstimators&) = delete; From a073e8f405b8f6662649486cc53703eb3af20c0d Mon Sep 17 00:00:00 2001 From: Bohdan Dudar Date: Thu, 4 Nov 2021 20:25:11 +0100 Subject: [PATCH 22/22] Removing TOFPlots --- TimeOfFlight/src/TOFPlots.cc | 771 ----------------------------------- 1 file changed, 771 deletions(-) delete mode 100644 TimeOfFlight/src/TOFPlots.cc diff --git a/TimeOfFlight/src/TOFPlots.cc b/TimeOfFlight/src/TOFPlots.cc deleted file mode 100644 index 80e5d589..00000000 --- a/TimeOfFlight/src/TOFPlots.cc +++ /dev/null @@ -1,771 +0,0 @@ -#include "TOFPlots.h" -#include -#include -#include - -#include -#include -#include -#include -#include - -// ----- include for verbosity dependend logging --------- -#include "marlin/VerbosityLevels.h" - -#include "marlin/ProcessorEventSeeder.h" -#include -#include - - -#include -#include -#include -#include - -#include -#include - -//---- ROOT ----- -#include "TH1F.h" - -using namespace lcio ; -using namespace marlin ; - - -TOFPlots aTOFPlots ; - - -// struct with cluster timing parameters -struct CluTime{ - - float clutime ; - float time05perc ; - float time10perc ; - float time20perc ; - float timeoffastesthit ; - float time05hits ; - float time10hits ; - float time20hits ; - float cor_ref_time05perc ; - float cor_ref_time10perc ; - float cor_ref_time20perc ; - float cor_ref_time05hits ; - float cor_ref_time10hits ; - float cor_ref_time20hits ; - -} ; - - -struct CorRefTime : lcrtrel::LCFloatExtension {} ; -struct SmearedTime : lcrtrel::LCFloatExtension {} ; - - -CluTime computeClusterTimes(EVENT::Cluster* clu, const float* refPoint, float timeResolution, gsl_rng* rng ) ; - - -// compute the path length along the track from the IP to the entry point into the calorimeter -float computeFlightLength( EVENT::Track* trk ) ; - - -// helper enum defining histogram index in vector -namespace TOFHistos { - enum index{ - hMCPEnergy, hPFOEnergy, hCluEnergy, hHitEnergy, hCorrectedTime, hCorrectedTimeFromRefPoint, - htimerefpoint100ps, htimerefpoint50ps, htimerefpoint10ps, haveragetime, haveragecorrectedtime, - haveragecorrectedreftime, - clutime , - time05perc , - time10perc , - time20perc ,timeoffastesthit, - time05hits , - time10hits , - time20hits , - cor_ref_time05perc , - cor_ref_time10perc , - cor_ref_time20perc , - cor_ref_time05hits , - cor_ref_time10hits , - cor_ref_time20hits , - flightlength , - h_beta_05hits , - h_beta_05perc , - h_beta_10hits , - h_beta_10perc , - h_beta_20hits , - h_beta_20perc , - hBeta_5hitsvsMomentum , - hBeta_5percvsMomentum , - hBeta_10hitsvsMomentum , - hBeta_10percvsMomentum , - hBeta_20hitsvsMomentum , - hBeta_20percvsMomentum , - //----- keep Size as last : - Size - }; - - - class Histograms{ - public: - Histograms(std::vector& v) : _h(&v) {} - - void create(int idx, const char* n, int nBin=100, double min=0., double max=0. ){ - create( idx , n , n , nBin, min , max ) ; - } - - void create(int idx, const char* n, const char* t, int nBin=100, double min=0., double max=0. ){ - - _h->at( idx ) = new TH1D( n, t , nBin , min, max ) ; - - streamlog_out( DEBUG ) << " create histo " << n << " at index " << idx << std::endl ; - } - - void create(int idx, const char* n, const char* t, int nBin , double* bins ){ - - _h->at( idx ) = new TH1D( n, t , nBin , bins ) ; - - streamlog_out( DEBUG ) << " create histo " << n << " at index " << idx << std::endl ; - } - - void fill( int idx , double val, double weight=1.0 ){ _h->at( idx )->Fill( val , weight ) ; } - - protected: - - std::vector* _h; - - }; -} - -using namespace TOFHistos ; - -TOFPlots::TOFPlots() : Processor("TOFPlots") { - - // modify processor description - _description = "TOFPlots creates plots related to potential time of flight measurements with the calorimter" ; - - - // register steering parameters: name, description, class-variable, default value - registerInputCollection( LCIO::MCPARTICLE, - "MCPCollectionName" , - "Name of the MCParticle collection" , - _colNameMCP , - std::string("MCParticle") - ); - - // register steering parameters: name, description, class-variable, default value - registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, - "PFOCollectionName" , - "Name of the ReconstructedParticle collection" , - _colNamePFO , - std::string("PandoraPFOs") - ); -} - - - -void TOFPlots::init() { - - streamlog_out(DEBUG) << " init called " << std::endl ; - - // usually a good idea to - printParameters() ; - - _nRun = 0 ; - _nEvt = 0 ; - - // initialize gsl random generator - _rng = gsl_rng_alloc(gsl_rng_ranlxs2); - - Global::EVENTSEEDER->registerProcessor(this); -} - - -void TOFPlots::processRunHeader( LCRunHeader* ) { - - _nRun++ ; -} - - - -void TOFPlots::processEvent( LCEvent * evt ) { - - - gsl_rng_set( _rng, Global::EVENTSEEDER->getSeed(this) ) ; - streamlog_out( DEBUG4 ) << "seed set to " << Global::EVENTSEEDER->getSeed(this) << std::endl; - - // this gets called for every event - // usually the working horse ... - - - // define a histogram pointer - - static AIDA::ICloud2D* hHitDistvsTime ; - static AIDA::ICloud2D* hHitDistRefvsTime ; - static AIDA::IHistogram2D* hBeta_5hitsvsMomentum ; - static AIDA::IHistogram2D* hBeta_5percvsMomentum ; - static AIDA::IHistogram2D* hBeta_10hitsvsMomentum ; - static AIDA::IHistogram2D* hBeta_10percvsMomentum ; - static AIDA::IHistogram2D* hBeta_20hitsvsMomentum ; - static AIDA::IHistogram2D* hBeta_20percvsMomentum ; - - Histograms h(_h) ; - - if( isFirstEvent() ){ - - // this creates a directory for this processor .... - AIDAProcessor::histogramFactory( this ) ; - - _h.resize( TOFHistos::Size ) ; - - const int nBins = 100 ; - - h.create(hMCPEnergy , "hMCPEnergy ", "energy of the MCParticles", nBins ) ; - h.create(hPFOEnergy , "hPFOEnergy ", "energy of the ReconstructedParticles", nBins ) ; - h.create(hCluEnergy , "hCluEnergy ", "energy of the cluster", nBins ) ; - h.create(hHitEnergy , "hHitEnergy ", "energy of the Hit", nBins ) ; - h.create(hCorrectedTime , "hCorrectedTime", "hittime compared to c", nBins ) ; - h.create(hCorrectedTimeFromRefPoint , "hCorrectedTimeFromRefPoint", "hittime from Refpoint compared to c", nBins ) ; - h.create(htimerefpoint100ps , "htimerefpoint100ps", "hittime from Refpoint compared to c with gaus100ps", nBins ) ; - h.create(htimerefpoint50ps , "htimerefpoint50ps", "hittime from Refpoint compared to c with gaus50ps", nBins ) ; - h.create(htimerefpoint10ps , "htimerefpoint10ps", "hittime from Refpoint compared to c with gaus10ps", nBins ) ; - h.create(haveragetime , "haveragetime", "average hittime per cluster", nBins , 7.2 , 8.5 ) ; - - h.create(clutime , "clutime_0ps" , "average cluster time (all hits) smeared with 0ps", nBins , 7.3 , 8. ) ; - h.create(time05perc ,"time05perc_0ps" , "average cluster time ( 5% of hits) smeared with 0ps", nBins , 7.0 , 7.6 ) ; - h.create(time10perc ,"time10perc_0ps" , "average cluster time (10% of hits) smeared with 0ps", nBins , 7.1 , 7.7 ) ; - h.create(time20perc ,"time20perc_0ps" , "average cluster time (20% of hits) smeared with 0ps", nBins , 7.1 , 7.7 ) ; - h.create(timeoffastesthit ,"timeoffastesthit_0ps" , "time of fastest hit smeared with 0ps", nBins , 7. ,7.5 ) ; - - h.create(time05hits ,"time05hits_0ps" , "average cluster time ( 5 fastest hits) smeared with 0ps", nBins , 7.1 , 7.6 ) ; - h.create(time10hits ,"time10hits_0ps" , "average cluster time (10 fastest hits) smeared with 0ps", nBins , 7.1 , 7.6 ) ; - h.create(time20hits ,"time20hits_0ps" , "average cluster time (20 fastest hits) smeared with 0ps", nBins , 7.2 , 7.7 ) ; - - h.create(cor_ref_time05perc ,"cor_ref_time05perc_0ps" , "average cluster cor_ref_time ( 5% of hits) smeared with 0ps", nBins , 7. ,7.5 ) ; - h.create(cor_ref_time10perc ,"cor_ref_time10perc_0ps" , "average cluster cor_ref_time (10% of hits) smeared with 0ps", nBins , 7. ,7.5 ) ; - h.create(cor_ref_time20perc ,"cor_ref_time20perc_0ps" , "average cluster cor_ref_time (20% of hits) smeared with 0ps", nBins , 7. ,7.5 ) ; - - h.create(cor_ref_time05hits ,"cor_ref_time05hits_0ps" , "average cluster cor_ref_time ( 5 fastest hits) smeared with 0ps", nBins , 7. ,7.5 ) ; - h.create(cor_ref_time10hits ,"cor_ref_time10hits_0ps" , "average cluster cor_ref_time (10 fastest hits) smeared with 0ps", nBins , 7. ,7.5 ) ; - h.create(cor_ref_time20hits ,"cor_ref_time20hits_0ps" , "average cluster cor_ref_time (20 fastest hits) smeared with 0ps", nBins , 7. ,7.5 ) ; - - h.create(haveragecorrectedtime , "haveragecorrectedtime", "average corrected hittime per cluster", nBins , 0. , 1. ) ; - h.create(haveragecorrectedreftime , "haveragecorrectedreftime", "average corrected hittime from refpoint per cluster", - nBins , 7. , 8. ) ; - - h.create(flightlength , "flightlength" , "total lenght of flight from the IP to the RefPoint", nBins , 1800 , 2200 ) ; - - h.create(h_beta_05hits , "beta_05hits_10ps" , "beta factor with cor_ref_time05hits smeared with 10ps", nBins , 0.98 , 1.05 ) ; - h.create(h_beta_05perc , "beta_05perc_10ps" , "beta factor with cor_ref_time05perc smeared with 10ps", nBins , 0.98 , 1.05 ) ; - h.create(h_beta_10hits , "beta_10hits_10ps" , "beta factor with cor_ref_time10hits smeared with 10ps", nBins , 0.98 , 1.05 ) ; - h.create(h_beta_10perc , "beta_10perc_10ps" , "beta factor with cor_ref_time10perc smeared with 10ps", nBins , 0.98 , 1.05 ) ; - h.create(h_beta_20hits , "beta_20hits_10ps" , "beta factor with cor_ref_time20hits smeared with 10ps", nBins , 0.98 , 1.05 ) ; - h.create(h_beta_20perc , "beta_20perc_10ps" , "beta factor with cor_ref_time20perc smeared with 10ps", nBins , 0.98 , 1.05 ) ; - - hHitDistvsTime = AIDAProcessor::histogramFactory(this)-> - createCloud2D( "hHitDistvsTime", "global r of Hit vs time", nBins ) ; - - hHitDistRefvsTime = AIDAProcessor::histogramFactory(this)-> - createCloud2D( "hHitDistRefvsTime", "dist of Hit from calo entry point vs time", nBins ) ; - - hBeta_5hitsvsMomentum = AIDAProcessor::histogramFactory(this)-> - createHistogram2D( "hBeta_5hitsvsMomentum_0ps", "Particles speed in c (beta_05hits) vs momentum smeared with 0ps", nBins, 1. , 10., nBins, 0.93 , 1.03 ) ; - - hBeta_5percvsMomentum = AIDAProcessor::histogramFactory(this)-> - createHistogram2D( "hBeta_5percvsMomentum_0ps", "Particles speed in c (beta_05perc) vs momentum smeared with 0ps", nBins, 1. , 10., nBins, 0.93 , 1.03 ) ; - - hBeta_10hitsvsMomentum = AIDAProcessor::histogramFactory(this)-> - createHistogram2D( "hBeta_10hitsvsMomentum_0ps", "Particles speed in c (beta_10hits) vs momentum smeared with 0ps", nBins, 1. , 10., nBins, 0.93 , 1.03 ) ; - - hBeta_10percvsMomentum = AIDAProcessor::histogramFactory(this)-> - createHistogram2D( "hBeta_10percvsMomentum_0ps", "Particles speed in c (beta_10perc) vs momentum smeared with 0ps", nBins, 1. , 10., nBins, 0.93 , 1.03 ) ; - - hBeta_20hitsvsMomentum = AIDAProcessor::histogramFactory(this)-> - createHistogram2D( "hBeta_20hitsvsMomentum_0ps", "Particles speed in c (beta_20hits) vs momentum smeared with 0ps", nBins, 1. , 10., nBins, 0.93 , 1.03 ) ; - - hBeta_20percvsMomentum = AIDAProcessor::histogramFactory(this)-> - createHistogram2D( "hBeta_20percvsMomentum_0ps", "Particles speed in c (beta_20perc) vs momentum smeared with 0ps", nBins, 1. , 10., nBins, 0.93 , 1.03 ) ; - - } - - LCCollection* colMCP=0 ; - LCCollection* colPFO=0 ; - - try{ colMCP = evt->getCollection( _colNameMCP ) ; } catch(lcio::Exception&){} - try{ colPFO = evt->getCollection( _colNamePFO ) ; } catch(lcio::Exception&){} - - - if( colMCP != NULL ){ - int nMCP = colMCP->getNumberOfElements() ; - - for(int i=0; i< nMCP ; i++){ - MCParticle* p = dynamic_cast( colMCP->getElementAt( i ) ) ; - - // only look at true generator particles - if( p->getGeneratorStatus() != 1 ) - continue ; - - // fill histogram from LCIO data : - h.fill( hMCPEnergy, p->getEnergy() ) ; - } - } - - - - if( colPFO != NULL ){ - - int nPFO = colPFO->getNumberOfElements() ; - - - // loop over all PFOs ---------------------------------------- - for(int i=0; i< nPFO ; i++){ - ReconstructedParticle* p = dynamic_cast( colPFO->getElementAt( i ) ) ; - - - - const double* mom = p->getMomentum() ; - - double momentum = sqrt( mom[0] * mom[0] + mom[1] * mom[1] + mom[2] * mom[2] ) ; - - // only interested in charged particles - if( std::fabs( p->getCharge() ) < 0.5 ) - continue ; - - // fill histogram from LCIO data : - h.fill( hPFOEnergy, p->getEnergy() ) ; - - // get clusters - const ClusterVec& cv = p->getClusters() ; - - // take only the frst - Cluster* clu = ( !cv.empty() ? cv[0] : 0 ) ; - - if( clu == 0 ){ - streamlog_out(DEBUG7) << " no cluster in PFO ...." << std::endl ; - continue ; - } - - - // get tracks - const TrackVec& tv = p->getTracks() ; - - // take only the frst - Track* trk = ( !tv.empty() ? tv[0] : 0 ) ; - - if( trk == 0 ){ - streamlog_out(DEBUG7) << " no track in PFO ...." << std::endl ; - continue ; - } - - const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ; - float x_ref = tscalo->getReferencePoint()[0] ; - float y_ref = tscalo->getReferencePoint()[1] ; - float z_ref = tscalo->getReferencePoint()[2] ; - - - - - h.fill( hCluEnergy, clu->getEnergy() ) ; - - // getCalorimeterHits - const CalorimeterHitVec& chv = clu->getCalorimeterHits() ; - - float total_time = 0 ; - float total_correctedtime = 0 ; - float total_correctedreftime = 0 ; - int number_of_hits = 0 ; - - //========================= loop over all hits in this particle's cluster ================ - for( CalorimeterHit* calohit : chv){ - - h.fill( hHitEnergy, calohit->getEnergy() ) ; - - float time = calohit->getTime() ; - - if( time < 1.e-3 ){ // ignore hits if time 0 or less than 1 ps - continue ; - } - - - long cellID = calohit->getCellID0() ; - UTIL::BitField64 bf("system:5") ; - bf.setValue(cellID) ; - int systemID = bf["system"] ; - - if( systemID != UTIL::ILDDetID::ECAL ){ - continue ; - } - - // compute distance of hit from the IP - float x = calohit->getPosition()[0] ; - float y = calohit->getPosition()[1] ; - float z = calohit->getPosition()[2] ; - - float distFromIP = std::sqrt( x*x + y*y + z*z ) ; - - float distFromRef = std::sqrt( (x-x_ref)*(x-x_ref) - + (y-y_ref)*(y-y_ref) - + (z-z_ref)*(z-z_ref) ) ; - - - if( time < 10. ){ - hHitDistvsTime->fill( time, distFromIP ) ; - hHitDistRefvsTime->fill( time, distFromRef ) ; - } - - // compute calohittime from Refpoint compared to a particle travelling with speed of light - - float hittimewithcfromrefpoint = distFromRef / 299.8 ; - - float correctedtimefromrefpoint = time - hittimewithcfromrefpoint ; - - - // get hit times assuming different resolutions - // time is given in ns - assume 10, 50 and 100 ps time resolution - - float smear_100ps = gsl_ran_gaussian( _rng, .1 ) ; - float timerefpoint100ps = correctedtimefromrefpoint + smear_100ps ; - float smear_50ps = gsl_ran_gaussian( _rng, .05 ) ; - float timerefpoint50ps = correctedtimefromrefpoint + smear_50ps ; - float smear_10ps = gsl_ran_gaussian( _rng, .01 ) ; - float timerefpoint10ps = correctedtimefromrefpoint + smear_10ps ; - - // compute calohittime compared to a particle travelling with speed of light - - float hittimewithc = std::sqrt( x*x + y*y + z*z )/299.8 ; - - float correctedtime = time - hittimewithc ; - - total_time += time ; - total_correctedtime += correctedtime ; - total_correctedreftime += correctedtimefromrefpoint ; - ++number_of_hits ; - - if( correctedtime < 0.3 ){ - h.fill( hCorrectedTime, correctedtime ) ; - h.fill( hCorrectedTimeFromRefPoint, correctedtimefromrefpoint ) ; - } - - if( timerefpoint100ps < 8.5 ){ - h.fill( htimerefpoint100ps, timerefpoint100ps ) ; - } - - if( timerefpoint50ps < 8.5 ){ - h.fill( htimerefpoint50ps, timerefpoint50ps ) ; - } - - if( timerefpoint10ps < 8.5 ){ - h.fill( htimerefpoint10ps, timerefpoint10ps ) ; - } - } - - float averagetime = total_time / (float) number_of_hits ; - float averagecorrectedtime = total_correctedtime / (float) number_of_hits ; - float averagecorrectedreftime = total_correctedreftime / (float) number_of_hits ; - - - double timeRes = 0.0 ; // time resolution in ns - CluTime ct = computeClusterTimes( clu, tscalo->getReferencePoint(), timeRes , _rng ) ; - - h.fill( clutime, ct.clutime ) ; - h.fill( time05perc, ct.time05perc ) ; - h.fill( time10perc, ct.time10perc ) ; - h.fill( time20perc, ct.time20perc ) ; - h.fill( time05hits, ct.time05hits ) ; - h.fill( time10hits, ct.time10hits ) ; - h.fill( time20hits, ct.time20hits ) ; - - h.fill( cor_ref_time05perc, ct.cor_ref_time05perc ) ; - h.fill( cor_ref_time10perc, ct.cor_ref_time10perc ) ; - h.fill( cor_ref_time20perc, ct.cor_ref_time20perc ) ; - h.fill( cor_ref_time05hits, ct.cor_ref_time05hits ) ; - h.fill( cor_ref_time10hits, ct.cor_ref_time10hits ) ; - h.fill( cor_ref_time20hits, ct.cor_ref_time20hits ) ; - h.fill( timeoffastesthit, ct.timeoffastesthit ) ; - - - float length = computeFlightLength( trk) ; - - h.fill( flightlength , length ) ; - - float beta_05hits = ( length / ct.cor_ref_time05hits ) / 299.8 ; - float beta_05perc = ( length / ct.cor_ref_time05perc ) / 299.8 ; - float beta_10hits = ( length / ct.cor_ref_time10hits ) / 299.8 ; - float beta_10perc = ( length / ct.cor_ref_time10perc ) / 299.8 ; - float beta_20hits = ( length / ct.cor_ref_time20hits ) / 299.8 ; - float beta_20perc = ( length / ct.cor_ref_time20perc ) / 299.8 ; - - h.fill( h_beta_05hits , beta_05hits ) ; - h.fill( h_beta_05perc , beta_05perc ) ; - h.fill( h_beta_10hits , beta_10hits ) ; - h.fill( h_beta_10perc , beta_10perc ) ; - h.fill( h_beta_20hits , beta_20hits ) ; - h.fill( h_beta_20perc , beta_20perc ) ; - - hBeta_5hitsvsMomentum->fill( momentum, beta_05hits ) ; - hBeta_5percvsMomentum->fill( momentum, beta_05perc ) ; - hBeta_10hitsvsMomentum->fill( momentum, beta_10hits ) ; - hBeta_10percvsMomentum->fill( momentum, beta_10perc ) ; - hBeta_20hitsvsMomentum->fill( momentum, beta_20hits ) ; - hBeta_20percvsMomentum->fill( momentum, beta_20perc ) ; - - - streamlog_out(DEBUG) << " averagetime : " << averagetime - << " and from computeClusterTimes: " << ct.cor_ref_time05perc << std::endl ; - - - h.fill( haveragetime, averagetime ) ; - h.fill( haveragecorrectedtime, averagecorrectedtime ) ; - h.fill( haveragecorrectedreftime, averagecorrectedreftime ) ; - - //========================================================================================= - } - } - - - - //-- note: this will not be printed if compiled w/o MARLINDEBUG=1 ! - - streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() - << " in run: " << evt->getRunNumber() << std::endl ; - - - - _nEvt ++ ; -} - - - -void TOFPlots::check( LCEvent *) { - // nothing to check here - could be used to fill checkplots in reconstruction processor -} - - - -void TOFPlots::end(){ - - gsl_rng_free( _rng ); - - - streamlog_out(MESSAGE) << "TOFPlots::end() " << name() - << " processed " << _nEvt << " events in " << _nRun << " runs " - << std::endl ; - -} - - -//******************************************************************************************************* - - -CluTime computeClusterTimes(EVENT::Cluster* clu, const float* refPoint, float timeResolution, gsl_rng* rng){ - - struct CluTime ct ; - - // getCalorimeterHits - const CalorimeterHitVec& cluHv = clu->getCalorimeterHits() ; - std::vector< EVENT::CalorimeterHit* > chv ; - chv.reserve( cluHv.size() ) ; - - float total_time = 0 ; - int number_of_hits = 0 ; - - - - //========================= loop over all hits in this particle's cluster ================ - for( CalorimeterHit* calohit : cluHv){ - - float time = calohit->getTime() ; - - // here we apply a random time shift to simulate the resolution - float rndTimeShift = gsl_ran_gaussian( rng, timeResolution ) ; - time += rndTimeShift ; - - calohit->ext() = time ; - - - if( time < 1.e-3 ){ // ignore hits if time 0 or less than 1 ps - continue ; - } - long cellID = calohit->getCellID0() ; - UTIL::BitField64 bf("system:5") ; - bf.setValue(cellID) ; - int systemID = bf["system"] ; - - if( systemID != UTIL::ILDDetID::ECAL ){ // only look at Ecal (barrel hits) - continue ; - } - - ++number_of_hits ; - - // copy the hit into the hit vector - chv.push_back( calohit ) ; - - total_time += time ; - } - ct.clutime = total_time / number_of_hits ; - - // sort the hit vector wrt time using a lambda expression - std::sort( chv.begin(), chv.end(), []( EVENT::CalorimeterHit* h1, EVENT::CalorimeterHit* h2) { - return h1->ext() < h2->ext() ; - }); - - - float total_time05perc = 0 ; - int number_of_hits05perc = 0 ; - - float total_time10perc = 0 ; - int number_of_hits10perc = 0 ; - - float total_time20perc = 0 ; - int number_of_hits20perc = 0 ; - - float total_time05hits = 0 ; - float total_time10hits = 0 ; - float total_time20hits = 0 ; - - // for( CalorimeterHit* calohit : chv){ - - for(size_t i=0 ; i < chv.size() ; ++i){ - - EVENT::CalorimeterHit* calohit = chv.at(i); - - float time = calohit->ext() ; - - // compute time corrected for time of flight from reference point - float x = calohit->getPosition()[0] ; - float y = calohit->getPosition()[1] ; - float z = calohit->getPosition()[2] ; - - float x_ref = refPoint[0] ; - float y_ref = refPoint[1] ; - float z_ref = refPoint[2] ; - - float distFromRef = std::sqrt( (x-x_ref)*(x-x_ref) - + (y-y_ref)*(y-y_ref) - + (z-z_ref)*(z-z_ref) ) ; - - calohit->ext() = time - distFromRef / 299.8 ; - - if( i < (0.05*chv.size () ) ){ - ++number_of_hits05perc ; - total_time05perc += time ; - } - - if( i < (0.10*chv.size () ) ){ - ++number_of_hits10perc ; - total_time10perc += time ; - } - - if( i < (0.20*chv.size () ) ){ - ++number_of_hits20perc ; - total_time20perc += time ; - } - - if( i < 5 ){ - total_time05hits += time ; - } - - if( i < 10 ){ - total_time10hits += time ; - } - - if( i < 20 ){ - total_time20hits += time ; - } - - } - - ct.time05perc = total_time05perc / number_of_hits05perc ; - ct.time10perc = total_time10perc / number_of_hits10perc ; - ct.time20perc = total_time20perc / number_of_hits20perc ; - ct.time05hits = total_time05hits / 5 ; - ct.time10hits = total_time10hits / 10 ; - ct.time20hits = total_time20hits / 20 ; - - ct.timeoffastesthit = ( chv.empty() ? 0. : chv[0]->getTime() ) ; - - - // sort the hit vector again - now wrt to corrected time using a lambda expression - std::sort( chv.begin(), chv.end(), []( EVENT::CalorimeterHit* h1, EVENT::CalorimeterHit* h2) { - return h1->ext() < h2->ext() ; - }); - - float total_cor_ref_time05perc = 0 ; - int number_of_cor_ref_hits05perc = 0 ; - - float total_cor_ref_time10perc = 0 ; - int number_of_cor_ref_hits10perc = 0 ; - - float total_cor_ref_time20perc = 0 ; - int number_of_cor_ref_hits20perc = 0 ; - - float total_cor_ref_time05hits = 0 ; - float total_cor_ref_time10hits = 0 ; - float total_cor_ref_time20hits = 0 ; - - for(size_t i=0 ; i < chv.size() ; ++i){ - - EVENT::CalorimeterHit* calohit = chv.at(i); - - double time = calohit->ext() ; - - - if( i < (0.05*chv.size () ) ){ - ++number_of_cor_ref_hits05perc ; - total_cor_ref_time05perc += time ; - } - - if( i < (0.10*chv.size () ) ){ - ++number_of_cor_ref_hits10perc ; - total_cor_ref_time10perc += time ; - } - - if( i < (0.20*chv.size () ) ){ - ++number_of_cor_ref_hits20perc ; - total_cor_ref_time20perc += time ; - } - - if( i < 5 ){ - total_cor_ref_time05hits += time ; - } - - if( i < 10 ){ - total_cor_ref_time10hits += time ; - } - - if( i < 20 ){ - total_cor_ref_time20hits += time ; - } - - } - - ct.cor_ref_time05perc = total_cor_ref_time05perc / number_of_cor_ref_hits05perc ; - ct.cor_ref_time10perc = total_cor_ref_time10perc / number_of_cor_ref_hits10perc ; - ct.cor_ref_time20perc = total_cor_ref_time20perc / number_of_cor_ref_hits20perc ; - ct.cor_ref_time05hits = total_cor_ref_time05hits / 5 ; - ct.cor_ref_time10hits = total_cor_ref_time10hits / 10 ; - ct.cor_ref_time20hits = total_cor_ref_time20hits / 20 ; - - - - - return ct ; -} - - - -float computeFlightLength( EVENT::Track* trk){ - - const TrackState* tsIP = trk->getTrackState( TrackState::AtIP ) ; - const TrackState* tscalo = trk->getTrackState( TrackState::AtCalorimeter ) ; - - float Phicalo = tscalo->getPhi() ; - float PhiIP = tsIP->getPhi() ; - - float Omega = tsIP->getOmega() ; - float tanL = tsIP->getTanLambda() ; - - - float length = (PhiIP-Phicalo)*(1/Omega) * sqrt( 1 + tanL * tanL ) ; - - - return length ; -}