diff --git a/MuonSpectrometer/MuonCnv/MuonByteStreamCnvTest/src/MdtDigitToMdtRDO.cxx b/MuonSpectrometer/MuonCnv/MuonByteStreamCnvTest/src/MdtDigitToMdtRDO.cxx index 689dfbe4ed13bbc0e31f9b9f28ed442974023c77..39af787e434c22de80ca27f6a33b793b17e0d802 100644 --- a/MuonSpectrometer/MuonCnv/MuonByteStreamCnvTest/src/MdtDigitToMdtRDO.cxx +++ b/MuonSpectrometer/MuonCnv/MuonByteStreamCnvTest/src/MdtDigitToMdtRDO.cxx @@ -16,6 +16,7 @@ #include <algorithm> #include <cmath> +#include <atomic> ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// @@ -183,8 +184,8 @@ StatusCode MdtDigitToMdtRDO::fill_MDTdata(const EventContext& ctx) const { // Iterate on the digits of the collection digit_iterator it_dig = mdtCollection->begin(); - for ( ; it_dig != mdtCollection->end() ; ++it_dig) - { + static std::atomic<bool> bisWarningPrinted = false; + for ( ; it_dig != mdtCollection->end() ; ++it_dig) { const MdtDigit* mdtDigit = *it_dig; Identifier channelId = mdtDigit->identify(); @@ -201,12 +202,13 @@ StatusCode MdtDigitToMdtRDO::fill_MDTdata(const EventContext& ctx) const { tdc, channel); if (!cabling) { - // as long as there is no BIS78 cabling, to avoid a hard crash, replace the tubeNumber + // as long as there is no BIS sMDT cabling, to avoid a hard crash, replace the tubeNumber // of tubes not covered in the cabling by 1 if (m_idHelperSvc->mdtIdHelper().stationName(channelId)==1 && std::abs(m_idHelperSvc->mdtIdHelper().stationEta(channelId))>6 && m_idHelperSvc->issMdt(channelId)) { - ATH_MSG_WARNING("Found BIS78 sMDT with tubeLayer="<<layer<<" and tubeNumber="<<tube<<". Setting to 1,1 for now..."); + if (!bisWarningPrinted) ATH_MSG_WARNING("Found BIS sMDT with tubeLayer="<<layer<<" and tubeNumber="<<tube<<". Setting to 1,1 until a proper cabling is implemented..."); + bisWarningPrinted=true; cabling = readCdo->getOnlineId(name, eta, phi, multilayer, 1, 1,subsystem, mrod, link, tdc, channel); } if (!cabling) { diff --git a/MuonSpectrometer/MuonDigitization/RPC_Digitization/src/RpcDigitizationTool.cxx b/MuonSpectrometer/MuonDigitization/RPC_Digitization/src/RpcDigitizationTool.cxx index 89afc5dfb54cf2036f9e66daedef3448aa8bc79e..7e38e167dd67ecde6f82ad7eb47a2a8758fef79a 100644 --- a/MuonSpectrometer/MuonDigitization/RPC_Digitization/src/RpcDigitizationTool.cxx +++ b/MuonSpectrometer/MuonDigitization/RPC_Digitization/src/RpcDigitizationTool.cxx @@ -51,6 +51,7 @@ #include <sstream> #include <iostream> #include <fstream> +#include <atomic> //12 charge points, 15 BetaGamma points, 180 efficiency points for fcp search namespace { @@ -587,16 +588,15 @@ StatusCode RpcDigitizationTool::doDigitization(const EventContext& ctx, RpcDigit double corrtimejitter = 0 ; if(m_CorrJitter>0.01) corrtimejitter = CLHEP::RandGaussZiggurat::shoot(rndmEngine,0.,m_CorrJitter); //correlated jitter - // std::cout<<" Correlated eta / phi jitter = "<<corrtimejitter<<std::endl; // handle here the special case where eta panel is dead => phi strip status (dead or eff.) cannot be resolved; // measured panel eff. will be used in that case and no phi strip killing will happen bool undefPhiStripStat = false; - std::vector<int> pcseta = PhysicalClusterSize(ctx, &idpaneleta, &hit, rndmEngine); //set to one for new algorithms + std::vector<int> pcseta = PhysicalClusterSize(ctx, &idpaneleta, &hit, rndmEngine); //set to one for new algorithms ATH_MSG_DEBUG ( "Simulated cluster on eta panel: size/first/last= "<<pcseta[0] <<"/"<<pcseta[1] <<"/"<<pcseta[2] ); std::vector<int> pcsphi = PhysicalClusterSize(ctx, &idpanelphi, &hit, rndmEngine); //set to one for new algorithms ATH_MSG_DEBUG ( "Simulated cluster on phi panel: size/first/last= "<<pcsphi[0]<<"/"<<pcsphi[1]<<"/"<<pcsphi[2] ); - + // create Identifiers Identifier atlasRpcIdeta = m_idHelper->channelID(stationName, stationEta, stationPhi, doubletR, doubletZ, doubletPhi,gasGap, 0, pcseta[1] ); @@ -610,8 +610,10 @@ StatusCode RpcDigitizationTool::doDigitization(const EventContext& ctx, RpcDigit // as long there is no proper implementation of the BI RPC cabling and digitisation, // skip this method to avoid hard crash of digitisation + static std::atomic<bool> biWarningPrinted = false; if (stationName.find("BI")!=std::string::npos) { - ATH_MSG_WARNING("skipping RPC DetectionEfficiency for BI"); + if (!biWarningPrinted) ATH_MSG_WARNING("skipping call of RPC DetectionEfficiency for BI as long as no proper cabling is implemented"); + biWarningPrinted=true; } else { ATH_CHECK(DetectionEfficiency(ctx, &atlasRpcIdeta,&atlasRpcIdphi, undefPhiStripStat, rndmEngine, particleLink)); } @@ -1061,7 +1063,7 @@ std::vector<int> RpcDigitizationTool::PhysicalClusterSize(const EventContext& ct else{ float xstripnorm=xstrip/30.; - result[0] = ClusterSizeEvaluation(ctx, id, xstripnorm, rndmEngine ); + result[0] = ClusterSizeEvaluation(ctx, id, xstripnorm, rndmEngine); int nstrips = ele->Nstrips(measuresPhi); // @@ -2095,11 +2097,16 @@ int RpcDigitizationTool::ClusterSizeEvaluation(const EventContext& ctx, const Id else if(stationName>50) index = index - 44 ; // BME and BOE 53 and 54 are at indices 7 and 8 + static std::atomic<bool> clusterIndexAPrinted = false; + static std::atomic<bool> clusterIndexCPrinted = false; if( !m_ClusterSize_fromCOOL ){ index += m_FracClusterSize1_A.size()/2*measuresPhi ; if( index>m_FracClusterSize1_A.size() || index>m_FracClusterSize2_A.size() || - index>m_FracClusterSizeTail_A.size() || index>m_MeanClusterSizeTail_A.size() ) { - ATH_MSG_WARNING ( "Index out of array in ClusterSizeEvaluation SideA " << index <<" statName "<<stationName) ; + index>m_FracClusterSizeTail_A.size() || index>m_MeanClusterSizeTail_A.size() ) { + if (stationName==0 || stationName==1) { + if (!clusterIndexAPrinted) ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideA " << index <<" statName "<<stationName<<" (ClusterSize not from COOL)"); + clusterIndexAPrinted=true; + } else ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideA " << index <<" statName "<<stationName<<" (ClusterSize not from COOL)"); return 1; } FracClusterSize1 = m_FracClusterSize1_A [index]; @@ -2110,9 +2117,12 @@ int RpcDigitizationTool::ClusterSizeEvaluation(const EventContext& ctx, const Id if(stationEta<0){ index += m_FracClusterSize1_C.size()/2*measuresPhi - m_FracClusterSize1_A.size()/2*measuresPhi ; if( index>m_FracClusterSize1_C.size() || index>m_FracClusterSize2_C.size() || - index>m_FracClusterSizeTail_C.size() || index>m_MeanClusterSizeTail_C.size() ) { - ATH_MSG_WARNING ( "Index out of array in ClusterSizeEvaluation SideC " << index <<" statName "<<stationName) ; - return 1; + index>m_FracClusterSizeTail_C.size() || index>m_MeanClusterSizeTail_C.size() ) { + if (stationName==0 || stationName==1) { + if (!clusterIndexCPrinted) ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideC " << index <<" statName "<<stationName<<" (ClusterSize not from COOL)"); + clusterIndexCPrinted=true; + } else ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideC " << index <<" statName "<<stationName<<" (ClusterSize not from COOL)"); + return 1; } FracClusterSize1 = m_FracClusterSize1_C [index]; @@ -2130,18 +2140,30 @@ int RpcDigitizationTool::ClusterSizeEvaluation(const EventContext& ctx, const Id int RPC_ProjectedTracks = 0; if( readCdo->getProjectedTracksMap().find(Id) != readCdo->getProjectedTracksMap().end()) RPC_ProjectedTracks = readCdo->getProjectedTracksMap().find(Id)->second ; - - if(readCdo->getFracClusterSize1Map().find(Id) != readCdo->getFracClusterSize1Map().end()) FracClusterSize1 = float(readCdo->getFracClusterSize1Map().find(Id)->second) ; - else ATH_MSG_INFO ( "FracClusterSize1 entry not found for id = " <<m_idHelper->show_to_string(*IdRpcStrip)<<" default will be used") ; + // the FracClusterSize maps are + static std::atomic<bool> fracCluster1Printed = false; + if(readCdo->getFracClusterSize1Map().find(Id) != readCdo->getFracClusterSize1Map().end()) FracClusterSize1 = float(readCdo->getFracClusterSize1Map().find(Id)->second) ; + else { + if (!fracCluster1Printed) ATH_MSG_WARNING("FracClusterSize1 entry not found for id = " <<m_idHelper->show_to_string(*IdRpcStrip)<<" (stationName="<<stationName<<") default will be used"); + fracCluster1Printed=true; + } + static std::atomic<bool> fracCluster2Printed = false; if(readCdo->getFracClusterSize2Map().find(Id) != readCdo->getFracClusterSize2Map().end()) FracClusterSize2 = float(readCdo->getFracClusterSize2Map().find(Id)->second) ; - else ATH_MSG_INFO ( "FracClusterSize2 entry not found for id = " <<m_idHelper->show_to_string(*IdRpcStrip)<<" default will be used") ; + else { + if (!fracCluster2Printed) ATH_MSG_WARNING("FracClusterSize2 entry not found for id = " <<m_idHelper->show_to_string(*IdRpcStrip)<<" (stationName="<<stationName<<") default will be used"); + fracCluster2Printed=true; + } ATH_MSG_DEBUG ( "FracClusterSize1 and 2 "<< FracClusterSize1 << " " << FracClusterSize2 ); FracClusterSizeTail = 1. - FracClusterSize1 - FracClusterSize2 ; + static std::atomic<bool> meanClusterPrinted = false; if(readCdo->getMeanClusterSizeMap().find(Id) != readCdo->getMeanClusterSizeMap().end()) MeanClusterSize = float(readCdo->getMeanClusterSizeMap().find(Id)->second) ; - else ATH_MSG_INFO ( "MeanClusterSize entry not found for id = " <<m_idHelper->show_to_string(*IdRpcStrip)<<" default will be used") ; + else { + if (!meanClusterPrinted) ATH_MSG_WARNING ("MeanClusterSize entry not found for id = " <<m_idHelper->show_to_string(*IdRpcStrip)<<" (stationName="<<stationName<<") default will be used"); + meanClusterPrinted=true; + } MeanClusterSizeTail = MeanClusterSize - FracClusterSize1 - 2*FracClusterSize2 ; @@ -2151,9 +2173,12 @@ int RpcDigitizationTool::ClusterSizeEvaluation(const EventContext& ctx, const Id if(RPC_ProjectedTracks<m_CutProjectedTracks || RPC_ProjectedTracks>10000000 || MeanClusterSize>m_CutMaxClusterSize || MeanClusterSize<=1 || FracClusterSizeTail < 0 || FracClusterSize1 < 0 || FracClusterSize2 < 0 || FracClusterSizeTail > 1 || FracClusterSize1 > 1 || FracClusterSize2 >1){ index += m_FracClusterSize1_A.size()/2*measuresPhi ; if( index>m_FracClusterSize1_A.size() || index>m_FracClusterSize2_A.size() || - index>m_FracClusterSizeTail_A.size() || index>m_MeanClusterSizeTail_A.size() ) { - ATH_MSG_WARNING ( "Index out of array in ClusterSizeEvaluation SideA " << index << " statName "<<stationName) ; - return 1; + index>m_FracClusterSizeTail_A.size() || index>m_MeanClusterSizeTail_A.size() ) { + if (stationName==0 || stationName==1) { + if (!clusterIndexAPrinted) ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideA " << index << " statName "<<stationName); + clusterIndexAPrinted=true; + } else ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideA " << index << " statName "<<stationName); + return 1; } FracClusterSize1 = m_FracClusterSize1_A [index]; FracClusterSize2 = m_FracClusterSize2_A [index]; @@ -2161,20 +2186,21 @@ int RpcDigitizationTool::ClusterSizeEvaluation(const EventContext& ctx, const Id MeanClusterSizeTail = m_MeanClusterSizeTail_A [index]; if(stationEta<0){ - index += m_FracClusterSize1_C.size()/2*measuresPhi - m_FracClusterSize1_A.size()/2*measuresPhi ; - if( index>m_FracClusterSize1_C.size() || index>m_FracClusterSize2_C.size() || - index>m_FracClusterSizeTail_C.size() || index>m_MeanClusterSizeTail_C.size() ) { - ATH_MSG_WARNING ( "Index out of array in ClusterSizeEvaluation SideC " << index << " statName "<<stationName ) ; - return 1; - } - - FracClusterSize1 = m_FracClusterSize1_C [index]; - FracClusterSize2 = m_FracClusterSize2_C [index]; - FracClusterSizeTail = m_FracClusterSizeTail_C [index]; - MeanClusterSizeTail = m_MeanClusterSizeTail_C [index]; + index += m_FracClusterSize1_C.size()/2*measuresPhi - m_FracClusterSize1_A.size()/2*measuresPhi ; + if( index>m_FracClusterSize1_C.size() || index>m_FracClusterSize2_C.size() || + index>m_FracClusterSizeTail_C.size() || index>m_MeanClusterSizeTail_C.size() ) { + if (stationName==0 || stationName==1) { + if (!clusterIndexCPrinted) ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideC " << index << " statName "<<stationName); + clusterIndexCPrinted=true; + } else ATH_MSG_WARNING("Index out of array in ClusterSizeEvaluation SideC " << index << " statName "<<stationName); + return 1; + } + FracClusterSize1 = m_FracClusterSize1_C [index]; + FracClusterSize2 = m_FracClusterSize2_C [index]; + FracClusterSizeTail = m_FracClusterSizeTail_C [index]; + MeanClusterSizeTail = m_MeanClusterSizeTail_C [index]; } } - } if(FracClusterSize1>1 )FracClusterSize1 =1.; diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/MuonStationIntersectSvc/MdtIntersectGeometry.h b/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/MuonStationIntersectSvc/MdtIntersectGeometry.h index 1cb2f0457dc63d1cca8d9237010e58aa15cdbb31..5abe556761ac52f4acc1984290589e7a2c24b117 100755 --- a/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/MuonStationIntersectSvc/MdtIntersectGeometry.h +++ b/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/MuonStationIntersectSvc/MdtIntersectGeometry.h @@ -38,7 +38,7 @@ namespace Muon { const Identifier& chamberId() const { return m_chid; } private: - double tubeLength( int ml, int layer, int tube ) const; + double tubeLength(const int ml, const int layer, const int tube) const; void init(MsgStream* msg); void fillDeadTubes(const MuonGM::MdtReadoutElement* mydetEl, MsgStream* msg); diff --git a/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/src/MdtIntersectGeometry.cxx b/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/src/MdtIntersectGeometry.cxx index de72e259b0dcb0952d2a68d3437bd275876337fe..1758729faf6f46574f51e83a98207b36e2d7b447 100755 --- a/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/src/MdtIntersectGeometry.cxx +++ b/MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonStationIntersectSvc/src/MdtIntersectGeometry.cxx @@ -12,6 +12,8 @@ #include "MuonIdHelpers/IMuonIdHelperSvc.h" #include "GeoModelUtilities/GeoGetIds.h" #include "MuonCondData/MdtCondDbData.h" +#include "TrkDriftCircleMath/MdtId.h" +#include <TString.h> // for Form namespace Muon{ @@ -74,7 +76,6 @@ namespace Muon{ double lineAngle = std::atan2(ldir.z(),ldir.y()); TrkDriftCircleMath::LocPos linePos( lpos.y(),lpos.z() ); - TrkDriftCircleMath::Line line( linePos, lineAngle ); const TrkDriftCircleMath::DCVec& dcs = m_mdtGeometry->tubesPassedByLine( line ); @@ -84,30 +85,34 @@ namespace Muon{ TrkDriftCircleMath::DCCit dit_end = dcs.end(); for( ; dit!=dit_end;++dit ){ + TrkDriftCircleMath::MdtId mdtId = dit->id(); + double xint = dxdy*( dit->position().x() - lpos.y() ) + lpos.x(); - Identifier tubeid = m_idHelperSvc->mdtIdHelper().channelID( m_chid, dit->id().ml()+1, dit->id().lay()+1, dit->id().tube()+1 ); + Identifier tubeid = m_idHelperSvc->mdtIdHelper().channelID( m_chid, mdtId.ml()+1, mdtId.lay()+1, mdtId.tube()+1 ); if( m_deadTubesML.find( m_idHelperSvc->mdtIdHelper().multilayerID(tubeid) ) != m_deadTubesML.end() ) { if( std::find( m_deadTubes.begin(), m_deadTubes.end(), tubeid ) != m_deadTubes.end() ) continue; } - double distWall = std::abs(xint) - 0.5*tubeLength( dit->id().ml(), dit->id().lay(), dit->id().tube() ); + double distWall = std::abs(xint) - 0.5*tubeLength( mdtId.ml(), mdtId.lay(), mdtId.tube() ); intersects.push_back( MuonTubeIntersect( tubeid, dit->dr(), distWall ) ); - } intersect.setTubeIntersects( intersects ); return intersect; } - double MdtIntersectGeometry::tubeLength( int ml, int layer, int tube ) const + double MdtIntersectGeometry::tubeLength(const int ml, const int layer, const int tube) const { + if (ml<0 || ml>1) throw std::runtime_error(Form("File: %s, Line: %d\nMdtIntersectGeometry::tubeLength() - got called with ml=%d which is definitely out of range", __FILE__, __LINE__,ml)); + if (layer<0 || layer>3) throw std::runtime_error(Form("File: %s, Line: %d\nMdtIntersectGeometry::tubeLength() - got called with layer=%d which is definitely out of range", __FILE__, __LINE__,layer)); + if (tube<0 || tube>119) throw std::runtime_error(Form("File: %s, Line: %d\nMdtIntersectGeometry::tubeLength() - got called with tube=%d which is definitely out of range", __FILE__, __LINE__,tube)); // shift by one to account for MuonGeoModel scheme - tube += 1; - layer += 1; + int theTube = tube+1; + int theLayer = layer+1; // handle case where first ml is dead - if( ml == 1 && !m_detElMl1 ) return m_detElMl0->getActiveTubeLength( layer, tube ); - if( ml == 0 ) return m_detElMl0->getActiveTubeLength( layer, tube ); - else return m_detElMl1->getActiveTubeLength( layer, tube ); + if( ml == 1 && !m_detElMl1 ) return m_detElMl0->getActiveTubeLength( theLayer, theTube ); + if( ml == 0 ) return m_detElMl0->getActiveTubeLength( theLayer, theTube ); + else return m_detElMl1->getActiveTubeLength( theLayer, theTube ); } diff --git a/Tracking/TrkUtilityPackages/TrkDriftCircleMath/TrkDriftCircleMath/MdtChamberGeometry.h b/Tracking/TrkUtilityPackages/TrkDriftCircleMath/TrkDriftCircleMath/MdtChamberGeometry.h index 44938082b8e319b01431609c36e0308de636e390..1cf1a79b2c237be07dbd58ee3c5731e5569dd648 100644 --- a/Tracking/TrkUtilityPackages/TrkDriftCircleMath/TrkDriftCircleMath/MdtChamberGeometry.h +++ b/Tracking/TrkUtilityPackages/TrkDriftCircleMath/TrkDriftCircleMath/MdtChamberGeometry.h @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ #ifndef DCMATH_MDTCHAMBERGEOMETRY_H @@ -25,7 +25,7 @@ namespace TrkDriftCircleMath { LocPos tube0ml0, LocPos tube0ml1, double tubeDist, double tubeStage, double layDist, double stationTheta ); - virtual ~MdtChamberGeometry(); + virtual ~MdtChamberGeometry()=default; void init(); void setGeometry(unsigned int nml, unsigned int nlay, diff --git a/Tracking/TrkUtilityPackages/TrkDriftCircleMath/src/MdtChamberGeometry.cxx b/Tracking/TrkUtilityPackages/TrkDriftCircleMath/src/MdtChamberGeometry.cxx index a5259b1252f6bd7d8a7d0d4a6930c63d2fcf64b2..6dce8c6487a8a9bfbabd8d267cc15e3c0c8766b5 100644 --- a/Tracking/TrkUtilityPackages/TrkDriftCircleMath/src/MdtChamberGeometry.cxx +++ b/Tracking/TrkUtilityPackages/TrkDriftCircleMath/src/MdtChamberGeometry.cxx @@ -1,5 +1,5 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ #include "TrkDriftCircleMath/MdtChamberGeometry.h" @@ -7,6 +7,7 @@ #include <iostream> #include <algorithm> +#include <TString.h> // for Form namespace TrkDriftCircleMath { @@ -57,8 +58,6 @@ namespace TrkDriftCircleMath { m_wasInit.push_back(1); m_crossedTubes.reserve(50); } - - MdtChamberGeometry::~MdtChamberGeometry() {} void MdtChamberGeometry::setGeometry(unsigned int nml, unsigned int nlay, unsigned int ntubesml0, unsigned int ntubesml1, @@ -78,19 +77,11 @@ namespace TrkDriftCircleMath { m_firstTube[1] = tube0ml1; m_allTubes.clear(); - if( m_nml > 2 ){ - std::cout << " invalid geometry: nml " << m_nml << std::endl; - m_validGeometry = false; - } - if( m_nlay > 4 ){ - std::cout << " invalid geometry: nlay " << m_nlay << std::endl; - m_validGeometry = false; - } - if( ntubesml0 > 78 || ntubesml1 > 78 ) { - std::cout << " invalid geometry: ntubes ml0 " << ntubesml0 - << " ntubes ml1 " << ntubesml1 << std::endl; - m_validGeometry = false; - } + if(m_nml<1 || m_nml>2) throw std::runtime_error(Form("File: %s, Line: %d\nMdtChamberGeometry::setGeometry() - got called with nml=%d which is definitely out of range", __FILE__, __LINE__,m_nml)); + if(m_nlay<1 || m_nlay>4) throw std::runtime_error(Form("File: %s, Line: %d\nMdtChamberGeometry::setGeometry() - got called with nlay=%d which is definitely out of range", __FILE__, __LINE__,m_nlay)); + if(ntubesml0<1 || ntubesml0>120) throw std::runtime_error(Form("File: %s, Line: %d\nMdtChamberGeometry::setGeometry() - got called with ntubesml0=%d which is definitely out of range", __FILE__, __LINE__,ntubesml0)); + // there can be chambers with only 1 multilayer. Then, the second multilayer will have ntubesml1=0 + if(ntubesml1>120) throw std::runtime_error(Form("File: %s, Line: %d\nMdtChamberGeometry::setGeometry() - got called with ntubesml1=%d which is definitely out of range", __FILE__, __LINE__,ntubesml1)); m_stationTheta = stationTheta; } @@ -99,7 +90,7 @@ namespace TrkDriftCircleMath { bool MdtChamberGeometry::validId( unsigned int ml, unsigned int lay,unsigned int tube) const { - if( !m_validGeometry) return false; + if(!m_validGeometry) return false; if( ml > 1 ) { std::cout << " ERROR wrong index: ml " << ml << " max " << m_nml << std::endl; @@ -122,16 +113,13 @@ namespace TrkDriftCircleMath { { if( !m_validGeometry) return m_allTubes; if( !m_allTubes.empty() ) return m_allTubes; - -// std::cout << " building chamber geometry mls " << m_nml << " lays " << m_nlay << " tubes ml 0 " << m_ntubeslay[0] -// << " ml 1 " << m_ntubeslay[1] << std::endl; + // create vector with all tubes in chamber for( unsigned int ml=0;ml<m_nml;++ml){ if( !m_wasInit[ml] ) continue; for( unsigned int lay=0;lay<m_nlay;++lay){ for( unsigned int tube=0;tube<m_ntubesml[ml];++tube ){ const LocPos& lp = tubePosition(ml,lay,tube); -// std::cout << " ml " << ml << " lay " << lay << " tube " << tube << " " << lp << std::endl; m_allTubes.push_back( lp ); } } @@ -142,9 +130,7 @@ namespace TrkDriftCircleMath { void MdtChamberGeometry::tubesPassedByLine( const Line& line, int inMultilayer, DCVec& crossedTubes ) const { m_resLine.set(line); - -// if( inMultilayer != -1 ) -// std::cout << " ----- tubesPassedByLine for ml " << inMultilayer << std::endl; + const LocPos& linepos = line.position(); const LocDir& linedir = line.direction(); double dxdy = fabs(linedir.y())> 0.0001 ? linedir.x()/linedir.y() : linedir.x()/0.0001; @@ -155,7 +141,6 @@ namespace TrkDriftCircleMath { // if indicated only scan single multilayer if( inMultilayer != -1 && inMultilayer != (int)ml ){ - //std::cout << " skippin ml " << ml << std::endl; continue; } for( unsigned int lay=0;lay<m_nlay;++lay){ @@ -169,66 +154,48 @@ namespace TrkDriftCircleMath { if( ctube >= (int)m_ntubesml[ml] ) ctube = m_ntubesml[ml]-1; if( inMultilayer != -1 ) - //std::cout << " y position layer " << ylay << " xpos intersect " << xintersect << " relpos " - //<< relpos << " tube " << ctube << " in ml " << ml << " lay " << lay << std::endl; - - //std::cout << " backward loop " << std::endl; for( int i = ctube-1;i>=0;--i ){ const LocPos& lp = tubePosition(ml,lay,i); double res = m_resLine.residual(lp); - //std::cout << " tube " << i << " pos " << lp << " res " << res << std::endl; if( std::abs(res) > m_tubeRad ) { - //std::cout << " tube has larger radius " << std::endl; if( m_tubeDist > 0 ){ if( res > m_tubeRad ) break; }else{ if( res < -m_tubeRad ) break; } - - //std::cout << " continue search " << std::endl; }else{ // if this is a chamber with only the second ml, set the ml index accordingly unsigned int actualMl = m_isSecondMultiLayer ? 1 : ml; crossedTubes.push_back( DriftCircle(lp,m_tubeRad,res,DriftCircle::EmptyTube,MdtId(m_id.isBarrel(),actualMl,lay,i),0) ); - //std::cout << " passed tube " << lp << " res " << res << " in ml " << ml << " lay " << lay << " tube " << i << std::endl; } } - //std::cout << " forward loop " << std::endl; for( int i =ctube;i < (int)m_ntubesml[ml];++i ){ const LocPos& lp = tubePosition(ml,lay,i); double res = m_resLine.residual(lp); - //std::cout << " tube " << i << " pos " << lp << " res " << res << std::endl; if( std::abs(res) > m_tubeRad ) { - //std::cout << " tube has larger radius " << std::endl; if( m_tubeDist > 0 ){ if( res < - m_tubeRad ) break; }else{ if( res > m_tubeRad ) break; } - //std::cout << " continue search " << std::endl; }else{ unsigned int actualMl = m_isSecondMultiLayer ? 1 : ml; crossedTubes.push_back( DriftCircle(lp,m_tubeRad,res,DriftCircle::EmptyTube,MdtId(m_id.isBarrel(),actualMl,lay,i),0) ); - //std::cout << " passed tube " << lp << " res " << res << " in ml " << ml << " lay " << lay << " tube " << i << std::endl; } } } } - //std::cout << " total number of passed tubes " << m_crossedTubes.size() << std::endl; - } + } const DCVec& MdtChamberGeometry::tubesPassedByLine( const Line& line, int inMultilayer ) const { m_crossedTubes.clear(); - if( !m_validGeometry){ - std::cout << " >>>>> invalid geometry <<<<< " << std::endl; - return m_crossedTubes; - } + if( !m_validGeometry) return m_crossedTubes; tubesPassedByLine(line,inMultilayer,m_crossedTubes);