diff --git a/Reconstruction/Digitisers/src/PlanarDigiAlg.cpp b/Reconstruction/Digitisers/src/PlanarDigiAlg.cpp deleted file mode 100644 index 203fec53b1e5c9da461da41df7718ee15d8ae3fb..0000000000000000000000000000000000000000 --- a/Reconstruction/Digitisers/src/PlanarDigiAlg.cpp +++ /dev/null @@ -1,338 +0,0 @@ -/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */ -#include "PlanarDigiAlg.h" -#include "GearSvc/IGearSvc.h" -#include "EventSeeder/IEventSeeder.h" -#include "plcio/MCParticleConst.h" - -/* -#include <EVENT/LCCollection.h> -#include <EVENT/SimTrackerHit.h> -#include <EVENT/MCParticle.h> -#include <IMPL/LCCollectionVec.h> -#include <IMPL/LCRelationImpl.h> -#include <IMPL/TrackerHitPlaneImpl.h> - -#include "UTIL/CellIDEncoder.h" -#include <UTIL/Operators.h> -*/ - -#include "UTIL/ILDConf.h" - -// STUFF needed for GEAR -#include <gear/GEAR.h> -#include "gear/gearsurf/MeasurementSurfaceStore.h" -#include "gear/gearsurf/MeasurementSurface.h" -#include "gear/gearsurf/ICoordinateSystem.h" -#include "gear/gearsurf/CartesianCoordinateSystem.h" - -#include <gsl/gsl_rng.h> -#include <gsl/gsl_randist.h> - -#include "CLHEP/Vector/TwoVector.h" - -#include <cmath> -#include <algorithm> - -DECLARE_COMPONENT( PlanarDigiAlg ) - -PlanarDigiAlg::PlanarDigiAlg(const std::string& name, ISvcLocator* svcLoc) - : GaudiAlgorithm(name, svcLoc), - _nEvt(0) -{ - //_description = "PlanarDigiAlg creates TrackerHits from SimTrackerHits, smearing them according to the input parameters." ; - - // Input collections - declareProperty("HeaderCol", _headerCol); - declareProperty("SimTrackHitCollection", _inColHdl, "Handle of the Input SimTrackerHit collection"); - - // Output collections - declareProperty("TrackerHitCollection", _outColHdl, "Handle of the TrackerHit output collection"); - declareProperty("SimTrkHitRelCollection", _outRelColHdl, "Handle of TrackerHit SimTrackHit relation collection"); -} - -StatusCode PlanarDigiAlg::initialize() -{ - - if( _resU.size() != _resV.size() ) { - fatal() << "Inconsistent number of resolutions given for U and V coordinate: " - << "ResolutionU :" << _resU.size() << " != ResolutionV : " << _resV.size() - << endmsg; - - return StatusCode::FAILURE; - } - - // get the GEAR manager - auto _gear = service<IGearSvc>("GearSvc"); - if ( !_gear ) { - error() << "Failed to find GearSvc ..." << endmsg; - return StatusCode::FAILURE; - } - _GEAR = _gear->getGearMgr(); - - // initialize gsl random generator - _rng = gsl_rng_alloc(gsl_rng_ranlxs2); - _SEEDER = service<IEventSeeder>("EventSeeder"); - _SEEDER->registerAlg(this); - - // zoujh: TODO - is this necessary? ... - //FIXME:SJA: if we want the surface store to be filled we need to create an instance of MarlinTrk implemented with KalTest/KalDet - //TODO:MarlinTrk::IMarlinTrkSystem* trksystem = MarlinTrk::Factory::createMarlinTrkSystem( "KalTest" , marlin::Global::GEAR , "" ) ; - //TODO:if( trksystem == 0 ) { - //TODO: throw EVENT::Exception( std::string(" Cannot initialize MarlinTrkSystem of Type: ") + std::string("KalTest" ) ) ; - //TODO:} - //TODO:trksystem->init() ; - //FIXME:SJA gear surface store has now been filled so we can dispose of the MarlinTrkSystem - //TODO:delete trksystem; - - return GaudiAlgorithm::initialize(); -} - -StatusCode PlanarDigiAlg::execute() -{ - auto header = _headerCol.get()->at(0); - int evtNo = header.getEventNumber(); - int runNo = header.getRunNumber(); - debug() << "Processing Run[" << runNo << "]::Event[" << evtNo << "]" << endmsg; - - unsigned int thisSeed = _SEEDER->getSeed(this, evtNo, runNo); - gsl_rng_set( _rng, thisSeed ) ; - debug() << "seed set to " << thisSeed << endmsg; - - auto STHcol = _inColHdl.get(); - if ( STHcol == nullptr ) { - debug() << "Collection " << _inColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg; - return StatusCode::SUCCESS; // or FAILURE ? - } - - // ** - // ** STHcol != 0 - // ** - - unsigned nCreatedHits=0; - unsigned nDismissedHits=0; - - auto trkhitVec = _outColHdl.createAndPut(); - auto relCol = _outRelColHdl.createAndPut(); - - //CellIDEncoder<TrackerHitPlaneImpl> cellid_encoder( lcio::ILDCellID0::encoder_string , trkhitVec ) ; - - int nSimHits = STHcol->size() ; - - int det_id = 0 ; - UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string ) ; - - if ( nSimHits>0 ) { - auto SimTHit = STHcol->at( 0 ) ; - encoder.setValue(SimTHit.getCellID0()) ; - det_id = encoder[lcio::ILDCellID0::subdet] ; - } - - if ( det_id == lcio::ILDDetID::VXD ){} - else if( det_id == lcio::ILDDetID::SIT ){} - else if( det_id == lcio::ILDDetID::SET ){} - else if( det_id == lcio::ILDDetID::FTD ){} - else { - fatal() << "unsupported detector ID = " << det_id - << ": file " << __FILE__ << " line " << __LINE__ - << endmsg; - return StatusCode::FAILURE; - } - - //smearing - debug() << " processing collection " << _inColHdl.fullKey() - << " with " << nSimHits << " hits ... " << endmsg ; - - - int i = 0; - for( auto SimTHit : *STHcol ) { - - const int celId = SimTHit->getCellID0() ; - - encoder.setValue(celId) ; - int side = encoder[lcio::ILDCellID0::side]; - int layer = encoder[lcio::ILDCellID0::layer]; - int module = encoder[lcio::ILDCellID0::module]; - int sensor = encoder[lcio::ILDCellID0::sensor]; - - debug() << "Hit = "<< i << " has celId " << celId << endmsg; - debug() << "side = " << side << endmsg; - debug() << "layerNumber = " << layer << endmsg; - debug() << "moduleNumber = " << module << endmsg; - debug() << "sensorNumber = " << sensor << endmsg; - - // //************************************************************ - // // Quick check if the MCParticle is in the list of MCParticles - // //************************************************************ - // - // if ( ! SimTHit->getMCParticle().isAvailable() ) { - // error() << " SimHit " << SimTHit << " Created by zero MCParticle which is not in the list of MCParticles: " - // << endmsg; - // continue; - // } - // - // //************************************************************ - // // Quick check if the MCParticle is of zero charge - // //************************************************************ - // - // if( abs( SimTHit->getMCParticle().getCharge()) < 0.01 ){ - // error() << " SimHit Created by zero charge particle: " - // << " Charge = " << SimTHit->getMCParticle().getCharge() - // << " EDep = " << SimTHit->getEDep() - // << " x = " << SimTHit->getPosition()[0] - // << " PDG = " << SimTHit->getMCParticle().getPDG() - // << " ID = " << SimTHit->getMCParticle().id() - // << endmsg; - // continue; - // } - - - float resU = ( _resU.size() > 1 ? _resU.value().at( layer ) : _resU.value().at(0) ) ; - float resV = ( _resV.size() > 1 ? _resV.value().at( layer ) : _resV.value().at(0) ) ; - - debug() << " --- will smear hit with resU = " << resU << " and resV = " << resV << endmsg; - - auto& pos = SimTHit->getPosition() ; - - double smearedPos[3]; - - // GearSurfaces::MeasurementSurface* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( SimTHit->getCellID0() ); - - gear::MeasurementSurface const* ms = _GEAR->getMeasurementSurfaceStore().GetMeasurementSurface( encoder.lowWord() );; - CLHEP::Hep3Vector globalPoint(pos[0],pos[1],pos[2]); - CLHEP::Hep3Vector localPoint = ms->getCoordinateSystem()->getLocalPoint(globalPoint); - CLHEP::Hep3Vector localPointSmeared = localPoint; - - - debug() << "Position of hit before smearing global: ( "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<< " ) " - << "local: ( " << localPoint.x() << " " << localPoint.y() << " " << localPoint.z() << " )" << endmsg; - - // A small check, if the hit is in the boundaries: - if ( !ms->isLocalInBoundary( localPoint ) ){ - - error() << "Hit local: ( " << localPoint.x() << " " << localPoint.y() << " " << localPoint.z() << " )" - << " is not within boundaries. Hit is skipped." << endmsg; - - nDismissedHits++; - continue; - } - - unsigned tries = 0; - - bool accept_hit = false; - - // Now try to smear the hit and make sure it is still on the surface - while( tries < 100 ) { - - if (tries > 0) { - debug() << "retry smearing for side" << side << " layer"<< layer<< " module" << module - << " sensor" << sensor << " : retries " << tries << endmsg; - } - - localPointSmeared.setX( localPoint.x() + gsl_ran_gaussian(_rng, resU) ); - localPointSmeared.setY( localPoint.y() + gsl_ran_gaussian(_rng, resV) ); - - //check if hit is in boundaries - if ( ms->isLocalInBoundary( localPointSmeared ) ) { - accept_hit = true; - break; - } - - tries++; - } - - if( accept_hit == false ) { - debug() << "hit could not be smeared within ladder after 100 tries: hit dropped" << endmsg; - continue; - } - - // for 1D strip measurements: set v to 0! Only the measurement in u counts! - if( _isStrip ) localPointSmeared[1] = 0. ; - - // convert back to global position for TrackerHitPlaneImpl - CLHEP::Hep3Vector globalPointSmeared = ms->getCoordinateSystem()->getGlobalPoint(localPointSmeared); - - debug() <<"Position of hit after smearing global: ( " - << globalPointSmeared.x() <<" "<< globalPointSmeared.y() <<" "<< globalPointSmeared.z() << " ) " - << "local: ( " - << localPointSmeared.x() << " " << localPointSmeared.y() << " " << localPointSmeared.z() << " )" - << endmsg; - - smearedPos[0] = globalPointSmeared.x(); - smearedPos[1] = globalPointSmeared.y(); - smearedPos[2] = globalPointSmeared.z(); - - //make the TrackerHitPlaneImpl - auto trkHit = trkhitVec->create(); - - trkHit->setCellID0( encoder.lowWord() ); - - trkHit->setPosition( smearedPos ) ; - - gear::CartesianCoordinateSystem* cartesian = dynamic_cast< gear::CartesianCoordinateSystem* >( ms->getCoordinateSystem() ); - CLHEP::Hep3Vector uVec = cartesian->getLocalXAxis(); - CLHEP::Hep3Vector vVec = cartesian->getLocalYAxis(); - - float u_direction[2] ; - u_direction[0] = uVec.theta(); - u_direction[1] = uVec.phi(); - - float v_direction[2] ; - v_direction[0] = vVec.theta(); - v_direction[1] = vVec.phi(); - - debug() << " U[0] = "<< u_direction[0] << " U[1] = "<< u_direction[1] - << " V[0] = "<< v_direction[0] << " V[1] = "<< v_direction[1] - << endmsg ; - - /* zoujh: TODO - generate TrackerHitPlane with podio - trkHit->setU( u_direction ) ; - trkHit->setV( v_direction ) ; - - trkHit->setdU( resU ) ; - - if( _isStrip ) trkHit->setdV( 0 ); // no error in v direction for strip hits as there is no meesurement information in v direction - else trkHit->setdV( resV ) ; - */ - - if( _isStrip ){ - trkHit->setType( UTIL::set_bit( trkHit->getType() , UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ) ) ; - } - - trkHit->setEDep( SimTHit->getEDep() ); - - // make the relation - auto rel = relCol->create(); - - float weight = 1.0; - - debug() <<" Set relation between " - << " sim hit " << SimTHit - << " to tracker hit " << trkHit - << " with a weight of " << weight - << endmsg; - - rel.setFrom( trkHit->getObjectID() ); - rel.setTo( SimTHit->getObjectID() ); - rel.setWeight( weight ); - - nCreatedHits++; - - debug() << "-------------------------------------------------------" << endmsg; - - ++i; - } - - debug() << "Created " << nCreatedHits << " hits, " - << nDismissedHits << " hits got dismissed for being out of boundary" - << endmsg; - - _nEvt ++ ; - - return StatusCode::SUCCESS; -} - -StatusCode PlanarDigiAlg::finalize() -{ - info() << "Processed " << _nEvt << " events " << endmsg; - return GaudiAlgorithm::finalize(); -} diff --git a/Reconstruction/Digitisers/src/PlanarDigiAlg.h b/Reconstruction/Digitisers/src/PlanarDigiAlg.h deleted file mode 100644 index 0b7b12b7d59a5b3ac089a4eb28a00089b2046d3c..0000000000000000000000000000000000000000 --- a/Reconstruction/Digitisers/src/PlanarDigiAlg.h +++ /dev/null @@ -1,91 +0,0 @@ -#ifndef PLANAR_DIGI_ALG_H -#define PLANAR_DIGI_ALG_H - -#include "FWCore/DataHandle.h" -#include "GaudiAlg/GaudiAlgorithm.h" -#include <gsl/gsl_rng.h> -#include "plcio/EventHeaderCollection.h" -#include "plcio/SimTrackerHitCollection.h" -#include "plcio/TrackerHitCollection.h" -#include "plcio/LCRelationCollection.h" - -/** ======= PlanarDigiProcessor / PlanarDigiAlg ========== <br> - * Creates TrackerHits from SimTrackerHits, smearing them according to the input parameters. - * The SimTrackerHits should come from a planar detector like VXD, SIT, SET or FTD. - * - * WARNING: this processor depends on correctly set CellID0s and is NOT backwards compatible to - * SimTrackerHit output with wrong CellID0s!!! - * - * The positions of "digitized" TrackerHits are obtained by gaussian smearing positions - * of SimTrackerHits in u and v direction. - * <h4>Input collections and prerequisites</h4> - * Processor requires a collection of SimTrackerHits <br> - * <h4>Output</h4> - * Processor produces collection of smeared TrackerHits<br> - * @param SimTrackHitCollectionName The name of input collection of SimTrackerHits <br> - * (default name VXDCollection) <br> - * @param TrackerHitCollectionName The name of output collection of smeared TrackerHits <br> - * (default name VTXTrackerHits) <br> - * @param SimTrkHitRelCollection The name of the TrackerHit SimTrackerHit relation collection <br> - * (default name VTXTrackerHitRelations) <br> - * @param ResolutionU resolution in direction of u (in mm) <br> - * (default value 0.004) <br> - * @param ResolutionV Resolution in direction of v (in mm) <br> - * (default value 0.004) <br> - * @param IsStrip whether the hits are 1 dimensional strip measurements <br> - * (default value false)<br> - * <br> - * - */ - -namespace gear { class GearMgr; } - -class IEventSeeder; - -class PlanarDigiAlg : public GaudiAlgorithm -{ - -public: - - PlanarDigiAlg(const std::string& name, ISvcLocator* svcLoc); - - /** Called at the begin of the job before anything is read. - * Use to initialize the processor, e.g. book histograms. - */ - virtual StatusCode initialize() ; - - /** Called for every event - the working horse. - */ - virtual StatusCode execute() ; - - /** Called after data processing for clean up. - */ - virtual StatusCode finalize() ; - - -protected: - - typedef std::vector<float> FloatVec; - - int _nEvt ; - - gear::GearMgr* _GEAR; - IEventSeeder * _SEEDER; - gsl_rng* _rng ; - - // resolution in direction of u - either one per layer or one for all layers - Gaudi::Property<FloatVec> _resU{ this, "ResolutionU", {0.0040} }; - // resolution in direction of v - either one per layer or one for all layers - Gaudi::Property<FloatVec> _resV{ this, "ResolutionV", {0.0040} }; - // whether hits are 1D strip hits - Gaudi::Property<bool> _isStrip{ this, "IsStrip", false }; - - // Input collections - DataHandle<plcio::EventHeaderCollection> _headerCol{"EventHeaderCol", Gaudi::DataHandle::Reader, this}; - DataHandle<plcio::SimTrackerHitCollection> _inColHdl{"VXDCollection", Gaudi::DataHandle::Reader, this}; - // Output collections - DataHandle<plcio::TrackerHitCollection> _outColHdl{"VXDTrackerHits", Gaudi::DataHandle::Writer, this}; - DataHandle<plcio::LCRelationCollection> _outRelColHdl{"VTXTrackerHitRelations", Gaudi::DataHandle::Reader, this}; -}; - -#endif