Skip to content

GitLab

Commit 664471aa authored by Sebastien Ponce's avatar Sebastien Ponce
Browse files

Moved DetDesc specific code from Boole to DeMuonChamber

parent 32c9b668
Pipeline #3335500 passed with stages
in 13 seconds
Hide whitespace changes
Inline Side-by-side
Muon/MuonAlgs/src/MuonDigitization.cpp
View file @ 664471aa
......@@ -252,30 +252,22 @@ void MuonDigitization::addChamberNoise() {
// the readout number is essentially meaningless...
// it is chamber noise.....
if ( numberOfNoiseHit == 0 ) continue;
DeMuonChamber& p_chamber = *m_muonDetector->getChamber( k, s, chamber );
auto p_chamber = m_muonDetector->getChamber( k, s, chamber );
for ( int hit = 1; hit <= numberOfNoiseHit; ++hit ) {
// define position of hit
// first gap number
const int gap = std::min( (int)( m_flatDist() * gapPerChamber ), m_muonDetector->gapsInRegion( k, s ) - 1 );
IPVolume* gasGap = p_chamber.getGasGap( gap );
IPVolume* gasGapLayer = p_chamber.getGasGapLayer( gap );
const auto* gapBox = dynamic_cast<const SolidBox*>( gasGap->lvolume()->solid() );
const double zhalfgap = gapBox->zHalfLength();
// then x&y
const double x = m_flatDist() * m_muonDetector->getSensAreaX( k, s );
const double y = m_flatDist() * m_muonDetector->getSensAreaY( k, s );
const double time = m_flatDist() * m_BXTime;
// traslate xyz e local to global..
const Gaudi::XYZPoint pointInLayer1 = gasGap->toMother( Gaudi::XYZPoint( x, y, -1 * zhalfgap + 0.1 ) );
const Gaudi::XYZPoint pointInLayer2 = gasGap->toMother( Gaudi::XYZPoint( x, y, zhalfgap - 0.1 ) );
const Gaudi::XYZPoint pointInCh1 = gasGapLayer->toMother( pointInLayer1 );
const Gaudi::XYZPoint pointInCh2 = gasGapLayer->toMother( pointInLayer2 );
const Gaudi::XYZPoint point1 = p_chamber.geometry()->toGlobal( pointInCh1 );
const Gaudi::XYZPoint point2 = p_chamber.geometry()->toGlobal( pointInCh2 );
const double z = ( point1.z() + point2.z() ) / 2.0;
const double tofOfLight = ( sqrt( x * x + y * y + z * z ) ) / 300.0;
const int gap = std::min( (int)( m_flatDist() * gapPerChamber ), m_muonDetector->gapsInRegion( k, s ) - 1 );
#ifdef USE_DD4HEP
auto [point1, point2] = p_chamber.getGapPoints( gap, x, y );
#else
auto [point1, point2] = p_chamber->getGapPoints( gap, x, y );
#endif
const double z = ( point1.z() + point2.z() ) / 2.0;
const double tofOfLight = ( sqrt( x * x + y * y + z * z ) ) / 300.0;
const auto hitMid = Gaudi::XYZPoint{( point1.x() + point2.x() ) / 2.0F, ( point1.y() + point2.y() ) / 2.0F,
( point1.z() + point2.z() ) / 2.0F};
const int sen = m_muonDetector->sensitiveVolumeID( hitMid );
......
Muon/MuonBackground/src/MuonBackground.cpp
View file @ 664471aa
......@@ -16,6 +16,7 @@
#include "Kernel/MuonTileID.h"
#include "MuonDet/DeMuonChamber.h"
#include "MuonDet/DeMuonDetector.h"
#include "MuonDet/MuonNamespace.h"
#include "Event/GenCollision.h"
#include "Event/GenHeader.h"
......@@ -90,12 +91,7 @@ public:
private:
using DeMuonChamberHandle = DeMuonChamber&;
std::optional<std::tuple<DeMuonChamberHandle, float, float, float>> getRadialPosition( int index, int station );
StatusCode calculateHitPosInGap( DeMuonChamberHandle pChamber, int gapNumber, float xpos, float ypos, float xSlope,
float ySlope, float averageZ, Gaudi::XYZPoint& entryGlobal,
Gaudi::XYZPoint& exitGlobal );
StatusCode calculateAverageGap( DeMuonChamberHandle pChamber, int gapNumberStart, int gapNumberStop, float xpos,
float ypos, float& zaverage );
int m_type;
std::string m_typeOfBackground;
......@@ -807,22 +803,7 @@ MuonBackground::getRadialPosition( int index, int station ) {
try {
auto& chamber = m_muonDetector->pos2StChamber( xpos, ypos, station );
// check n ot only that hit is inside chamber but also gap...
Gaudi::XYZPoint pToCheck = chamber.geometry()->toGlobal( Gaudi::XYZPoint( 0, 0, 0 ) );
if ( msgLevel( MSG::DEBUG ) ) debug() << "to check " << pToCheck << endmsg;
IPVolume* firstGap = chamber.getFirstGasGap();
IPVolume* firstGapLayer = chamber.getGasGapLayer( 0 );
Gaudi::XYZPoint pInMother = firstGap->toMother( Gaudi::XYZPoint( 0, 0, 0 ) );
Gaudi::XYZPoint pInChamber = firstGapLayer->toMother( pInMother );
if ( msgLevel( MSG::DEBUG ) ) debug() << pInChamber << " p in chamber " << endmsg;
Gaudi::XYZPoint pInGlobal = ( chamber.geometry() )->toGlobal( pInChamber );
Gaudi::XYZPoint point( xpos, ypos, 0 );
if ( msgLevel( MSG::DEBUG ) ) debug() << point << " " << pInGlobal.z() << endmsg;
point.SetZ( pInGlobal.z() );
Gaudi::XYZPoint pointInChamber = ( chamber.geometry() )->toLocal( point );
Gaudi::XYZPoint pointInLayer = firstGapLayer->toLocal( pointInChamber );
if ( firstGap->isInside( pointInLayer ) ) {
if ( chamber.checkHitAndGapInChamber( xpos, ypos ) ) {
// remember this is the chamber number inside a region....
if ( msgLevel( MSG::DEBUG ) )
debug() << " chamber number " << chamber.chamberNumber() << " station number " << chamber.stationNumber()
......@@ -900,18 +881,14 @@ StatusCode MuonBackground::createHit( LHCb::MCHits* hitsContainer, int station,
}
if ( msgLevel( MSG::VERBOSE ) ) verbose() << " local slope " << xSlope << " " << ySlope << endmsg;
// define the z of the average gaps position
float averageZ = 0;
StatusCode sc = calculateAverageGap( pChamber, firstGap, lastGap, xpos, ypos, averageZ );
if ( msgLevel( MSG::VERBOSE ) ) verbose() << sc << " " << averageZ << endmsg;
float averageZ = pChamber.calculateAverageGap( firstGap, lastGap, xpos, ypos );
for ( int igap = 0; igap <= multi; igap++ ) {
int gapNumber = gapHit[igap];
Gaudi::XYZPoint entryGlobal;
Gaudi::XYZPoint exitGlobal;
// DeMuonGasGap* p_Gap=NULL;
sc = calculateHitPosInGap( pChamber, gapNumber, xpos, ypos, xSlope, ySlope, averageZ, entryGlobal, exitGlobal );
if ( msgLevel( MSG::DEBUG ) ) debug() << "status code of calhitpos " << sc << endmsg;
if ( sc.isSuccess() ) {
bool success =
pChamber.calculateHitPosInGap( gapNumber, xpos, ypos, xSlope, ySlope, averageZ, entryGlobal, exitGlobal );
if ( success ) {
if ( msgLevel( MSG::DEBUG ) ) debug() << "found hit in chamber " << endmsg;
} else {
allHitsInsideCha = false;
......@@ -930,20 +907,18 @@ StatusCode MuonBackground::createHit( LHCb::MCHits* hitsContainer, int station,
} else {
timeBest = time;
}
int _firstGap = gapHit[0];
int _lastGap = gapHit[gapHit.size() - 1];
float averageZ = 0;
StatusCode sc = calculateAverageGap( pChamber, _firstGap, _lastGap, xpos, ypos, averageZ );
if ( sc.isFailure() ) return sc;
int _firstGap = gapHit[0];
int _lastGap = gapHit[gapHit.size() - 1];
float averageZ = pChamber.calculateAverageGap( _firstGap, _lastGap, xpos, ypos );
for ( int igap = 0; igap <= multi; igap++ ) {
bool correct = true;
int gapNumber = gapHit[igap];
Gaudi::XYZPoint entryGlobal;
Gaudi::XYZPoint exitGlobal;
StatusCode _sc =
calculateHitPosInGap( pChamber, gapNumber, xpos, ypos, xSlope, ySlope, averageZ, entryGlobal, exitGlobal );
if ( _sc.isFailure() ) return _sc;
bool _sc =
pChamber.calculateHitPosInGap( gapNumber, xpos, ypos, xSlope, ySlope, averageZ, entryGlobal, exitGlobal );
if ( !_sc ) return StatusCode::FAILURE;
double x = ( entryGlobal.x() + exitGlobal.x() ) / 2.0;
double y = ( entryGlobal.y() + exitGlobal.y() ) / 2.0;
double z = ( entryGlobal.z() + exitGlobal.z() ) / 2.0;
......@@ -1014,149 +989,6 @@ int MuonBackground::howManyHit( float floatHit ) {
}
}
StatusCode MuonBackground::calculateHitPosInGap( DeMuonChamberHandle pChamber, int gapNumber, float xpos, float ypos,
float xSlope, float ySlope, float zavegaps,
Gaudi::XYZPoint& entryGlobal, Gaudi::XYZPoint& exitGlobal ) {
// StatusCode sc=StatusCode::FAILURE;
// IDetectorElement::IDEContainer::iterator itGap;
// int countGap=0;
// p_Gap=NULL;
// for(itGap=(pChamber)->childBegin(); itGap<(pChamber)->childEnd(); itGap++){
// if(countGap==gapNumber){
// p_Gap=dynamic_cast<DeMuonGasGap*>(*((*itGap)->childBegin()));
// break;
// }
// countGap++;
// }
IPVolume* p_Gap = pChamber.getGasGap( gapNumber );
IPVolume* p_GapLayer = pChamber.getGasGapLayer( gapNumber );
if ( p_Gap != NULL ) {
const ISolid* gapSolid = ( p_Gap )->lvolume()->solid();
const SolidBox* gapBox = dynamic_cast<const SolidBox*>( gapSolid );
const double zhalfgap = gapBox->zHalfLength();
const double xhalfgap = gapBox->xHalfLength();
const double yhalfgap = gapBox->yHalfLength();
// verbose()<<"half gap size "<<xhalfgap <<" "<<yhalfgap<<" "<<zhalfgap
// <<endmsg;
// verbose()<<"glob gap pos "<<gapcc.x()<<" "<<gapcc.y()<<" "
// <<gapcc.z()<<endmsg;
Gaudi::XYZPoint loch = pChamber.geometry()->toLocal( Gaudi::XYZPoint( xpos, ypos, zavegaps ) );
Gaudi::XYZPoint logaslayer = p_GapLayer->toLocal( loch );
Gaudi::XYZPoint poslocal = p_Gap->toLocal( logaslayer );
const double zcenter = poslocal.z();
// verbose()<<" input data "<<xpos<<" "<<ypos<<" "<<gapNumber<<" "
// <<zavegaps<<" "<<endmsg;
// verbose()<<"center gap "<<poslocal.x()<<" "<<poslocal.y()<<" "
// <<poslocal.z()<<" "<<zcenter<<endmsg;
// Gaudi::Transform3D matrixToLocal= (p_Gap->geometry())->matrix();
// Gaudi::XYZPoint slopelocal=matrixToLocal*
Gaudi::XYZPoint slopelocal = Gaudi::XYZPoint( xSlope, ySlope, 1.0F );
const double zinf = -zhalfgap;
const double zsup = +zhalfgap;
// verbose()<<"local slopes "<<
// slopelocal.x()<<" "<<slopelocal.y()<<" "<<slopelocal.z()<<endmsg;
double xentry = poslocal.x() + ( zinf - zcenter ) * ( slopelocal.x() / slopelocal.z() );
double xexit = poslocal.x() + ( zsup - zcenter ) * ( slopelocal.x() / slopelocal.z() );
double yentry = poslocal.y() + ( zinf - zcenter ) * ( slopelocal.y() / slopelocal.z() );
double yexit = poslocal.y() + ( zsup - zcenter ) * ( slopelocal.y() / slopelocal.z() );
// check that gap boundaries have not been crossed
const double xmin = -xhalfgap;
const double ymin = -yhalfgap;
const double xmax = xhalfgap;
const double ymax = yhalfgap;
double zmin = -zhalfgap;
double zmax = zhalfgap;
bool correct = true;
if ( xentry < xmin || xentry > xmax ) correct = false;
if ( yentry < ymin || yentry > ymax ) correct = false;
if ( xexit < xmin || xexit > xmax ) correct = false;
if ( yexit < ymin || yexit > ymax ) correct = false;
// verbose()<<"partial "<<xentry<<" "<<yentry<<" "
// <<xexit<<" "<<yexit<<" "<<correct<<endmsg;
// verbose()<<" correct "<<correct<<" "<<slopelocal.x()<<" "<<
// slopelocal.y()<<endmsg;
if ( correct ) {
// then x
if ( slopelocal.x() == 0 ) return StatusCode::FAILURE;
double z1 = ( xmin - poslocal.x() ) / ( slopelocal.x() / slopelocal.z() );
double z2 = ( xmax - poslocal.x() ) / ( slopelocal.x() / slopelocal.z() );
double zz1 = std::min( z1, z2 );
double zz2 = std::max( z1, z2 );
zmin = std::max( zmin, zz1 );
zmax = std::min( zmax, zz2 );
// then y
if ( slopelocal.y() == 0 ) return StatusCode::FAILURE;
z1 = ( xmin - poslocal.y() ) / ( slopelocal.y() / slopelocal.z() );
z2 = ( xmax - poslocal.y() ) / ( slopelocal.y() / slopelocal.z() );
zz1 = std::min( z1, z2 );
zz2 = std::max( z1, z2 );
zmin = std::max( zmin, zz1 );
zmax = std::min( zmax, zz2 );
if ( zmin > zhalfgap || zmax < -zhalfgap ) return StatusCode::FAILURE;
xentry = poslocal.x() + ( slopelocal.x() / slopelocal.z() ) * ( zmin - zcenter );
xexit = poslocal.x() + ( slopelocal.x() / slopelocal.z() ) * ( zmax - zcenter );
yentry = poslocal.y() + ( slopelocal.y() / slopelocal.z() ) * ( zmin - zcenter );
yexit = poslocal.y() + ( slopelocal.y() / slopelocal.z() ) * ( zmax - zcenter );
// back to the global frame
// verbose()<<xentry<<" "<<yentry<<" "<<zmin<<endmsg;
// verbose()<<xexit<<" "<<yexit<<" "<<zmax<<endmsg;
Gaudi::XYZPoint entryInLayer = p_Gap->toMother( Gaudi::XYZPoint( xentry, yentry, zmin ) );
Gaudi::XYZPoint entryInCh = p_GapLayer->toMother( entryInLayer );
entryGlobal = ( pChamber.geometry() )->toGlobal( entryInCh );
Gaudi::XYZPoint exitInLayer = p_Gap->toMother( ( Gaudi::XYZPoint( xexit, yexit, zmax ) ) );
Gaudi::XYZPoint exitInCh = p_GapLayer->toMother( exitInLayer );
exitGlobal = ( pChamber.geometry() )->toGlobal( exitInCh );
// verbose()<<"global "<<entryGlobal.x()<<" "<<entryGlobal.y()<<" "<<entryGlobal.z()<<endmsg;
// verbose()<<"global "<<exitGlobal.x()<<" "<<exitGlobal.y()<<" "<<exitGlobal.z()<<endmsg;
return StatusCode::SUCCESS;
} else
return StatusCode::FAILURE;
} else
return StatusCode::FAILURE;
}
StatusCode MuonBackground::calculateAverageGap( DeMuonChamberHandle pChamber, int gapNumberStart, int gapNumberStop,
float xpos, float ypos, float& zaverage ) {
StatusCode sc = StatusCode::SUCCESS;
int countGap = 0;
zaverage = 0;
for ( int i = gapNumberStart; i <= gapNumberStop; i++ ) {
IPVolume* pGapLayer = pChamber.getGasGapLayer( i );
IPVolume* pGap = pChamber.getGasGap( i );
// get 3 points to build a plane for the gap center
Gaudi::XYZPoint point1 =
pChamber.geometry()->toGlobal( pGapLayer->toMother( pGap->toMother( Gaudi::XYZPoint( 0., 0., 0. ) ) ) );
Gaudi::XYZPoint point2 =
pChamber.geometry()->toGlobal( pGapLayer->toMother( pGap->toMother( Gaudi::XYZPoint( 1., 0., 0. ) ) ) );
Gaudi::XYZPoint point3 =
pChamber.geometry()->toGlobal( pGapLayer->toMother( pGap->toMother( Gaudi::XYZPoint( 0., 1., 0. ) ) ) );
Gaudi::Plane3D plane( point1, point2, point3 );
double zInterceptThisPlane = 0.;
Gaudi::XYZVector para = plane.Normal();
if ( para.Z() != 0 ) {
double a = para.X();
double b = para.Y();
double d = plane.HesseDistance();
zInterceptThisPlane = -( a * xpos + b * ypos + d ) / ( para.Z() );
} else
zInterceptThisPlane = point1.z();
zaverage = zaverage + (float)zInterceptThisPlane;
}
countGap++;
zaverage = zaverage / ( gapNumberStop - gapNumberStart + 1 );
return sc;
}
int MuonBackground::chamberOffset( int station, int region ) { return m_chamberInRegion[station * 4 + region]; }
int MuonBackground::numberOfCollision( const LHCb::MCVertex* pointVertex ) {
......
Markdown is supported
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