Skip to content
Snippets Groups Projects
Commit a37b037f authored by FU Chengdong's avatar FU Chengdong
Browse files

move to Digitisers/SimpleDigi

parent 52ac3921
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Reconstruction/Digitisers/src/PlanarDigiAlg.cpp deleted 100644 → 0
+ 0
− 338
View file @ 52ac3921
/* -*- 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();
}
Reconstruction/Digitisers/src/PlanarDigiAlg.h deleted 100644 → 0
+ 0
− 91
View file @ 52ac3921
#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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment