| Report problems to ATLAS LXR Team (with time and IP address indicated) |
|
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ] |
||||
|
||||||
| Links to LXR source navigation pages for stable releases | [ 12.*.* ] [ 13.*.* ] [ 14.*.* ] [ 15.*.* ] | |||||
001 #include "GaudiKernel/ToolFactory.h" 002 #include "GaudiKernel/SmartDataPtr.h" 003 #include "GaudiKernel/IDataProviderSvc.h" 004 005 //----------------------------------------------------------------// 006 #include "CLHEP/Matrix/SymMatrix.h" 007 008 #include "TrkParameters/MeasuredAtaDisc.h" 009 #include "TrkParameters/MeasuredAtaCylinder.h" 010 #include "TrkParameters/MeasuredAtaPlane.h" 011 012 #include "TrkEventPrimitives/JacobianThetaPToCotThetaPt.h" 013 #include "TrkEventPrimitives/ParamDefs.h" 014 #include "TrkEventPrimitives/LocalDirection.h" 015 #include "TrkEventPrimitives/LocalPosition.h" 016 017 #include "TrkParameters/AtaPlane.h" 018 019 #include "TrkSurfaces/CylinderSurface.h" 020 #include "TrkSurfaces/DiscSurface.h" 021 022 #include "TrkEventPrimitives/DefinedParameter.h" 023 #include "TrkEventPrimitives/LocalParameters.h" 024 025 #include "Identifier/Identifier.h" 026 027 #include "TrkMeasurementBase/MeasurementBase.h" 028 #include "TrkParameters/MeasuredTrackParameters.h" 029 030 #include "MuonSegment/MuonSegment.h" 031 #include "MuonSegment/MuonSegmentQuality.h" 032 033 //----------------------------------------------------------------// 034 #include "AtlasDetDescr/AtlasDetectorID.h" 035 #include "MuonIdHelpers/MdtIdHelper.h" 036 #include "MuonIdHelpers/RpcIdHelper.h" 037 #include "MuonIdHelpers/TgcIdHelper.h" 038 #include "MuonIdHelpers/CscIdHelper.h" 039 040 #include "Particle/TrackParticle.h" 041 #include "Particle/TrackParticleContainer.h" 042 043 #include "TrkParameters/MeasuredPerigee.h" 044 #include "TrkRIO_OnTrack/RIO_OnTrack.h" 045 #include "TrkExInterfaces/IExtrapolator.h" 046 #include "TrkParameters/AtaDisc.h" 047 #include "TrkParameters/AtaCylinder.h" 048 #include "TrkSegment/Segment.h" 049 #include "InDetBeamSpotService/IBeamCondSvc.h" 050 051 ///////////////////////////////////////////////////////// 052 #include "STACOTools/StacoEDMHelper.h" 053 054 StacoEDMHelper::StacoEDMHelper(const std::string& t, 055 const std::string& n, 056 const IInterface* p ):AthAlgTool(t,n,p), 057 p_IExtrapolator("Trk::Extrapolator/AtlasExtrapolator"), 058 p_IBeamCondSvc ( "BeamCondSvc",t ) 059 { 060 061 declareInterface<IStacoEDMHelper>(this); 062 063 declareProperty("IExtrapolator" , p_IExtrapolator ) ; 064 065 declareProperty("IBeamCondSvc", p_IBeamCondSvc); 066 067 } 068 069 StacoEDMHelper::~StacoEDMHelper(){} 070 071 // Initialize 072 StatusCode StacoEDMHelper::initialize(){ 073 074 StatusCode sc = StatusCode::SUCCESS; 075 076 msg(MSG::INFO) << "Initialisation started " << endreq; 077 078 sc = AthAlgTool::initialize(); 079 if ( sc.isFailure() ) { 080 msg(MSG::FATAL) << " AthAlgTool::initialize() failed" << endreq; 081 return( StatusCode::FAILURE ); 082 } 083 084 //Set pointer on BeamCondSvc 085 m_IsOnIBeamCondSvc = 1 ; 086 if ( p_IBeamCondSvc.retrieve().isFailure() || p_IBeamCondSvc == 0 ){ 087 m_IsOnIBeamCondSvc = 0 ; 088 msg(MSG::WARNING) << "Could not find BeamCondSvc." << endreq; 089 msg(MSG::WARNING) << "MS track can not be propgated down to beam spot" << endreq; 090 msg(MSG::WARNING) << "z axis will used as perigee axis" << endreq; 091 }else{ 092 msg(MSG::INFO) << "Retrieved service " << p_IBeamCondSvc << endreq; 093 } 094 095 //Set pointer on AtlasDetectorID 096 sc = detStore()->retrieve(m_detID, "AtlasID" ); 097 if (sc.isFailure()) { 098 msg(MSG::FATAL) <<"Could not get AtlasDetectorID "<<endreq; 099 return( StatusCode::FAILURE ); 100 } 101 msg(MSG::INFO) << "Found AtlasDetectorID "<<endreq; 102 103 sc = detStore()->retrieve(m_mdtId,"MDTIDHELPER"); 104 if (sc.isFailure()){ 105 msg(MSG::FATAL) << "Cannot retrieve MDTIDHELPER" << endreq; 106 return( StatusCode::FAILURE ); 107 } 108 sc = detStore()->retrieve(m_cscId,"CSCIDHELPER"); 109 if (sc.isFailure()){ 110 msg(MSG::FATAL) << "Cannot retrieve CSCIDHELPER" << endreq; 111 return( StatusCode::FAILURE ); 112 } 113 sc = detStore()->retrieve(m_rpcId,"RPCIDHELPER"); 114 if (sc.isFailure()){ 115 msg(MSG::FATAL) << "Cannot retrieve RPCIDHELPER" << endreq; 116 return( StatusCode::FAILURE ); 117 } 118 sc = detStore()->retrieve(m_tgcId,"TGCIDHELPER"); 119 if (sc.isFailure()){ 120 msg(MSG::FATAL) << "Cannot retrieve TGCIDHELPER" << endreq; 121 return( StatusCode::FAILURE ); 122 } 123 124 //Retrieve p_IExtrapolator 125 if ( p_IExtrapolator.retrieve().isFailure() ) { 126 msg(MSG::WARNING) << "Failed to retrieve tool " << p_IExtrapolator << endreq; 127 p_IExtrapolator = 0 ; 128 }else{ 129 msg(MSG::INFO) << "Retrieved tool " << p_IExtrapolator << endreq; 130 } 131 132 msg(MSG::INFO) << "Initialisation ended " << endreq; 133 134 return StatusCode::SUCCESS; 135 136 } 137 138 StatusCode StacoEDMHelper::finalize(){return StatusCode::SUCCESS;} 139 140 std::vector<const Rec::TrackParticle*> StacoEDMHelper::AlignTrackParticleContainers( 141 const Rec::TrackParticleContainer* pTrackParticleContainer , 142 const Rec::TrackParticleContainer* pTrackParticleContainer_Aux 143 ){ 144 145 //Minimal distance for identification 146 double DistanceSmall = 0.001; 147 148 //Collect in AuxVecColl_Aligned the TP from pTrackParticleContainer_Aux corresponding to those in pTrackParticleContainer 149 std::vector<const Rec::TrackParticle*> AuxVecColl_Aligned; 150 151 152 //Copy pTrackParticleContainer_Aux in AuxVecColl 153 std::vector<const Rec::TrackParticle*> AuxVecColl; 154 Rec::TrackParticleContainer::const_iterator it_Aux = pTrackParticleContainer_Aux->begin() ; 155 Rec::TrackParticleContainer::const_iterator it_Aux_End = pTrackParticleContainer_Aux->end() ; 156 for ( ; it_Aux!=it_Aux_End; ++it_Aux ){ 157 AuxVecColl.push_back(*it_Aux); 158 } 159 160 161 //Loop on pTrackParticleContainer 162 Rec::TrackParticleContainer::const_iterator it = pTrackParticleContainer->begin() ; 163 Rec::TrackParticleContainer::const_iterator it_End = pTrackParticleContainer->end() ; 164 for ( ; it!=it_End; ++it ){ 165 166 // std::cout 167 // << "----------" 168 // << "----------" 169 // << "----------" 170 // << "----------" 171 // << "----------" 172 // << "----------" 173 // << "----------" 174 // << "----------" 175 // << std::endl ; 176 177 std::vector<const Rec::TrackParticle*>::iterator it_AuxVecColl_Associated = AuxVecColl.end(); 178 179 // Loop on AuxVecColl 180 std::vector<const Rec::TrackParticle*>::iterator it_AuxVecColl = AuxVecColl.begin(); 181 std::vector<const Rec::TrackParticle*>::iterator it_AuxVecColl_End = AuxVecColl.end(); 182 for ( ; it_AuxVecColl!= it_AuxVecColl_End ;++it_AuxVecColl ) { 183 if (it_AuxVecColl_Associated != it_AuxVecColl_End ) break ; 184 185 // std::cout 186 // << "----------" 187 // << "----------" 188 // << "----------" 189 // << "----------" 190 // << std::endl ; 191 192 // Get the TrkTrack 193 const Trk::Track* pTrack = (*it)->originalTrack(); 194 const Trk::Track* pTrack_Aux = (*it_AuxVecColl)->originalTrack(); 195 196 // Loop on the TrackParameters to decide if the tracks match 197 std::vector<const Trk::TrackParameters*>::const_iterator TrkPar = pTrack->trackParameters()->begin(); 198 std::vector<const Trk::TrackParameters*>::const_iterator TrkPar_End = pTrack->trackParameters()->end(); 199 for ( ; TrkPar!=TrkPar_End; ++TrkPar){ 200 if (it_AuxVecColl_Associated != it_AuxVecColl_End ) break ; 201 const Trk::MeasuredTrackParameters* MEAS_TrkPar = dynamic_cast<const Trk::MeasuredTrackParameters*>(*TrkPar); 202 if ( !MEAS_TrkPar ) continue ; 203 204 std::vector<const Trk::TrackParameters*>::const_iterator TrkPar_Aux = pTrack_Aux->trackParameters()->begin(); 205 std::vector<const Trk::TrackParameters*>::const_iterator TrkPar_Aux_End = pTrack_Aux->trackParameters()->end(); 206 for ( ; TrkPar_Aux!=TrkPar_Aux_End; ++TrkPar_Aux){ 207 if (it_AuxVecColl_Associated != it_AuxVecColl_End ) break ; 208 const Trk::MeasuredTrackParameters* MEAS_TrkPar_Aux = dynamic_cast<const Trk::MeasuredTrackParameters*>(*TrkPar_Aux); 209 if ( !MEAS_TrkPar_Aux ) continue ; 210 211 double Dx = ((*TrkPar)->position()).x() - ((*TrkPar_Aux)->position()).x() ; 212 double Dy = ((*TrkPar)->position()).y() - ((*TrkPar_Aux)->position()).y() ; 213 double Dz = ((*TrkPar)->position()).z() - ((*TrkPar_Aux)->position()).z() ; 214 double Distance = sqrt ( Dx * Dx + Dy * Dy + Dz * Dz ) ; 215 if ( Distance <= DistanceSmall ) it_AuxVecColl_Associated = it_AuxVecColl ; 216 217 // std::cout 218 // << " " << Distance 219 // << " " << Dx 220 // << " " << Dy 221 // << " " << Dz 222 // << " " << ((*TrkPar)->position()).x() 223 // << " " << ((*TrkPar)->position()).y() 224 // << " " << ((*TrkPar)->position()).z() 225 // << " " << ((*TrkPar_Aux)->position()).x() 226 // << " " << ((*TrkPar_Aux)->position()).y() 227 // << " " << ((*TrkPar_Aux)->position()).z() 228 // << std::endl ; 229 230 } 231 } 232 233 } 234 235 // Collect the corresponding menber of AuxVecColl or 0 236 if (it_AuxVecColl_Associated != it_AuxVecColl_End ) { 237 AuxVecColl_Aligned.push_back(*it_AuxVecColl_Associated); 238 AuxVecColl.erase(it_AuxVecColl_Associated); 239 }else{ 240 AuxVecColl_Aligned.push_back(0); 241 } 242 243 } 244 245 return AuxVecColl_Aligned; 246 247 } 248 249 const Trk::TrackParameters* StacoEDMHelper::GetSecuredAtaPlane( 250 const Trk::GlobalPosition & aGlobalPosition , 251 const Trk::GlobalMomentum & aGlobalMomentum , 252 double TheCharge , 253 const Trk::PlaneSurface & aPlaneSurface 254 ){ 255 256 // return ( new Trk::AtaPlane (aGlobalPosition, aGlobalMomentum , TheCharge , aPlaneSurface) ); 257 258 double locx = 0. ; 259 double locy = 0. ; 260 261 const Trk::LocalPosition* pLocalPosition = aPlaneSurface.globalToLocal( aGlobalPosition ); 262 if (pLocalPosition) { 263 locx = (*pLocalPosition)[Trk::locX] ; 264 locy = (*pLocalPosition)[Trk::locY] ; 265 }else{ 266 msg(MSG::DEBUG) << " globalToLocal failed in GetSecuredAtaPlane " << endreq; 267 Trk::GlobalPosition aLocalPosition = aPlaneSurface.transform().inverse()*aGlobalPosition ; 268 locx = aLocalPosition[Trk::locX] ; 269 locy = aLocalPosition[Trk::locY] ; 270 } 271 delete pLocalPosition; 272 double qoverp = TheCharge/fabs(aGlobalMomentum.mag()); 273 return ( new Trk::AtaPlane(locx,locy,aGlobalMomentum.phi(),aGlobalMomentum.theta(),qoverp, aPlaneSurface) ); 274 275 } 276 int StacoEDMHelper::HitsPerML( const Trk::Segment* pSegment, int ML ){ 277 if ( IsValid(pSegment) == 0 ) return 0; 278 279 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 280 if( !pMuonSegment ) return 0; 281 282 int hitsOnML(0); 283 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 284 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 285 for (; rotIter!=pRIOSet.end(); ++rotIter){ 286 if ( (*rotIter)!=0 ){ 287 Identifier id = (*rotIter)->identify() ; 288 if ( m_mdtId->is_mdt( id ) ) 289 if( ML == m_mdtId->multilayer( id ) ) ++hitsOnML; 290 } 291 } 292 293 return hitsOnML; 294 } 295 296 std::string StacoEDMHelper::EtaStationString( const Trk::Segment* pSegment ){ 297 std::string s("0"); 298 if( IsValid(pSegment) == 0 ) return s; 299 300 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 301 if( !pMuonSegment ) return s; 302 303 std::stringstream ss; 304 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 305 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 306 for (; rotIter!=pRIOSet.end(); ++rotIter){ 307 if ( (*rotIter)!=0 ){ 308 Identifier id = (*rotIter)->identify() ; 309 if( m_mdtId->is_mdt(id) ) ss << m_mdtId->stationEta(id); 310 if( m_cscId->is_csc(id) ) ss << m_cscId->stationEta(id); 311 ss >> s; 312 return s; 313 } 314 } 315 return s; 316 317 } 318 // Dumps TrkTrack 319 StatusCode StacoEDMHelper::dump_TrkTrack(const Trk::Track* pTrack,std::ofstream* pOutCurrent){ 320 321 StatusCode sc = StatusCode::SUCCESS; 322 323 std::vector<const Trk::TrackParameters*>::const_iterator TrkPar = pTrack->trackParameters()->begin(); 324 std::vector<const Trk::TrackParameters*>::const_iterator TrkParEnd = pTrack->trackParameters()->end(); 325 for ( ; TrkPar!=TrkParEnd; ++TrkPar){ 326 const Trk::Perigee* pPerigee = dynamic_cast<const Trk::Perigee*>(*TrkPar); 327 if ( pPerigee !=0 ) { 328 sc = dump_Perigee(pPerigee,pOutCurrent) ; 329 if (sc.isFailure()) msg(MSG::WARNING) << "dump_Perigee failed" << endreq; 330 } 331 const Trk::AtaCylinder* pAtaCylinder = dynamic_cast<const Trk::AtaCylinder*>(*TrkPar); 332 if ( pAtaCylinder !=0 ) { 333 sc = dump_AtaCylinder(pAtaCylinder,pOutCurrent) ; 334 if (sc.isFailure()) msg(MSG::WARNING) << "dump_AtaCylinder failed" << endreq; 335 } 336 const Trk::AtaDisc* pAtaDisc = dynamic_cast<const Trk::AtaDisc*>(*TrkPar); 337 if ( pAtaDisc !=0 ) { 338 sc = dump_AtaDisc(pAtaDisc,pOutCurrent) ; 339 if (sc.isFailure()) msg(MSG::WARNING) << "dump_AtaDisc failed" << endreq; 340 } 341 if ( pPerigee ==0 && pAtaCylinder ==0 && pAtaDisc ==0 ) { 342 *pOutCurrent << "* NOT Perigee, AtaCylinder, AtaDisc"<<std::endl; 343 *pOutCurrent << *(*TrkPar) << std::endl; 344 } 345 } 346 347 *pOutCurrent << *pTrack << std::endl; 348 349 return StatusCode::SUCCESS; 350 351 } 352 StatusCode StacoEDMHelper::dump_Perigee(const Trk::Perigee* pPerigee,std::ofstream* pOutCurrent){ 353 354 if (pPerigee == 0) return StatusCode::SUCCESS; 355 *pOutCurrent << "* Its parameters at Perigee are:"<<std::endl; 356 if (IsThisPositionBeyondThelimitCaloSpectro( 357 (pPerigee->position()).x(), 358 (pPerigee->position()).y(), 359 (pPerigee->position()).z()) == 1 ) { 360 *pOutCurrent << "* At Muon Spectro Entrance ?"<<std::endl; 361 }else{ 362 *pOutCurrent << "* At beam?"<<std::endl; 363 } 364 365 *pOutCurrent << "*parameters()[Trk::d0] " << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::d0] << std::endl; 366 *pOutCurrent << " parameters()[Trk::z0] " << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::z0] << std::endl; 367 *pOutCurrent << " parameters()[Trk::phi0] " << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::phi0] << std::endl; 368 *pOutCurrent << " parameters()[Trk::theta] " << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::theta] << std::endl; 369 *pOutCurrent << " parameters()[Trk::qOverP] " << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::qOverP] << std::endl; 370 *pOutCurrent << " cotTheta() " << std::setw(12)<<std::setprecision(5)<<pPerigee->cotTheta() << std::endl; 371 *pOutCurrent << " charge()/peri->pT() " << std::setw(12)<<std::setprecision(5)<<pPerigee->charge()/pPerigee->pT() << std::endl; 372 373 const Trk::MeasuredPerigee* pMeasuredPerigee = dynamic_cast<const Trk::MeasuredPerigee*>(pPerigee); 374 if (pMeasuredPerigee != 0){ 375 double mtheta = (pPerigee->parameters())[Trk::theta]; 376 double mqp = (pPerigee->parameters())[Trk::qOverP]; 377 Trk::JacobianThetaPToCotThetaPt TheJac(mtheta,mqp); 378 Trk::CovarianceMatrix covVert(5); 379 covVert = pMeasuredPerigee->localErrorMatrix().covariance().similarity(TheJac); 380 *pOutCurrent 381 << "* at Perigee : " 382 << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::d0] 383 << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::z0] 384 << std::setw(12)<<std::setprecision(5)<<pPerigee->parameters()[Trk::phi0] 385 << std::setw(12)<<std::setprecision(5)<<pPerigee->cotTheta() 386 << std::setw(12)<<std::setprecision(5)<<pPerigee->charge()/pPerigee->pT() 387 << std::setw(12)<<std::setprecision(5)<<pPerigee->pT()/pPerigee->charge() 388 <<std::endl; 389 *pOutCurrent 390 << "* +/- " 391 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[0][0]) 392 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[1][1]) 393 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[2][2]) 394 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[3][3]) 395 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[4][4]) 396 <<std::endl; 397 *pOutCurrent 398 << "* Covariance matrice " 399 <<std::endl; 400 *pOutCurrent 401 << "* : " 402 << std::setw(12)<<std::setprecision(2)<<covVert[0][0]*100./sqrt(covVert[0][0]*covVert[0][0]) 403 << std::setw(12)<<std::setprecision(2)<<covVert[0][1]*100./sqrt(covVert[0][0]*covVert[1][1]) 404 << std::setw(12)<<std::setprecision(2)<<covVert[0][2]*100./sqrt(covVert[0][0]*covVert[2][2]) 405 << std::setw(12)<<std::setprecision(2)<<covVert[0][3]*100./sqrt(covVert[0][0]*covVert[3][3]) 406 << std::setw(12)<<std::setprecision(2)<<covVert[0][4]*100./sqrt(covVert[0][0]*covVert[4][4]) 407 <<std::endl; 408 for (int i1=1; i1<=4; i1++){ 409 *pOutCurrent 410 << "* " 411 << std::setw(12)<<std::setprecision(2)<<covVert[i1][0]*100./sqrt(covVert[i1][i1]*covVert[0][0]) 412 << std::setw(12)<<std::setprecision(2)<<covVert[i1][1]*100./sqrt(covVert[i1][i1]*covVert[1][1]) 413 << std::setw(12)<<std::setprecision(2)<<covVert[i1][2]*100./sqrt(covVert[i1][i1]*covVert[2][2]) 414 << std::setw(12)<<std::setprecision(2)<<covVert[i1][3]*100./sqrt(covVert[i1][i1]*covVert[3][3]) 415 << std::setw(12)<<std::setprecision(2)<<covVert[i1][4]*100./sqrt(covVert[i1][i1]*covVert[4][4]) 416 <<std::endl; 417 } 418 } 419 *pOutCurrent 420 << "* position " 421 << std::setw(12)<<std::setprecision(3)<<(pPerigee->position()).x() 422 <<std::endl; 423 *pOutCurrent 424 << "* " 425 << std::setw(12)<<std::setprecision(3)<<(pPerigee->position()).y() 426 <<std::endl; 427 *pOutCurrent 428 << "* " 429 << std::setw(12)<<std::setprecision(3)<<(pPerigee->position()).z() 430 <<std::endl; 431 *pOutCurrent 432 << "* momentum " 433 << std::setw(12)<<std::setprecision(3)<<(pPerigee->momentum()).x() 434 <<std::endl; 435 *pOutCurrent 436 << "* " 437 << std::setw(12)<<std::setprecision(3)<<(pPerigee->momentum()).y() 438 <<std::endl; 439 *pOutCurrent 440 << "* " 441 << std::setw(12)<<std::setprecision(3)<<(pPerigee->momentum()).z() 442 <<std::endl; 443 444 return StatusCode::SUCCESS; 445 446 } 447 StatusCode StacoEDMHelper::dump_AtaCylinder(const Trk::AtaCylinder* pAtaCylinder,std::ofstream* pOutCurrent){ 448 449 if (pAtaCylinder == 0) return StatusCode::SUCCESS; 450 *pOutCurrent << "* Its parameters at AtaCylinder are:"<<std::endl; 451 if (IsThisPositionBeyondThelimitCaloSpectro( 452 (pAtaCylinder->position()).x(), 453 (pAtaCylinder->position()).y(), 454 (pAtaCylinder->position()).z()) == 1 ) { 455 *pOutCurrent << "* At Muon Spectro Entrance ?"<<std::endl; 456 }else{ 457 *pOutCurrent << "* At Calo Entrance ?"<<std::endl; 458 } 459 460 *pOutCurrent << "*parameters()[Trk::locRPhi] " << std::setw(12)<<std::setprecision(5)<<pAtaCylinder->parameters()[Trk::locRPhi] << std::endl; 461 *pOutCurrent << " parameters()[Trk::locZ] " << std::setw(12)<<std::setprecision(5)<<pAtaCylinder->parameters()[Trk::locZ] << std::endl; 462 *pOutCurrent << " parameters()[Trk::phi] " << std::setw(12)<<std::setprecision(5)<<pAtaCylinder->parameters()[Trk::phi] << std::endl; 463 *pOutCurrent << " parameters()[Trk::theta] " << std::setw(12)<<std::setprecision(5)<<pAtaCylinder->parameters()[Trk::theta] << std::endl; 464 *pOutCurrent << " parameters()[Trk::qOverP] " << std::setw(12)<<std::setprecision(5)<<pAtaCylinder->parameters()[Trk::qOverP] << std::endl; 465 466 double RCrossing = sqrt ( (pAtaCylinder->position()).x()*(pAtaCylinder->position()).x() + (pAtaCylinder->position()).y()*(pAtaCylinder->position()).y() ); 467 //double ZCrossing = (pAtaCylinder->position()).z() ; 468 double ThetaX = atan( RCrossing/pAtaCylinder->parameters()[Trk::locZ] ) ; 469 if ( ThetaX < 0. ) ThetaX = ThetaX + 3.1415927 ; 470 double PhiX = pAtaCylinder->parameters()[Trk::locRPhi]/RCrossing ; 471 double ThetaV = pAtaCylinder->parameters()[Trk::theta] ; 472 double PhiV = pAtaCylinder->parameters()[Trk::phi] ; 473 double qOnp = pAtaCylinder->parameters()[Trk::qOverP] ; 474 475 const Trk::MeasuredAtaCylinder* pMeasuredAtaCylinder = dynamic_cast<const Trk::MeasuredAtaCylinder*>(pAtaCylinder); 476 if (pMeasuredAtaCylinder != 0){ 477 HepMatrix TheJac(5,5,0); 478 TheJac(1,2) = -std::pow( sin(ThetaX) , 2 )/RCrossing ; 479 TheJac(2,1) = 1./RCrossing ; 480 TheJac(3,4) = 1. ; 481 TheJac(4,3) = 1. ; 482 TheJac(5,5) = 1. ; 483 Trk::CovarianceMatrix covVert(5); 484 covVert = pMeasuredAtaCylinder->localErrorMatrix().covariance().similarity(TheJac); 485 *pOutCurrent 486 << "* at AtaCylinder : " 487 << std::setw(12)<<std::setprecision(5)<<ThetaX 488 << std::setw(12)<<std::setprecision(5)<<PhiX 489 << std::setw(12)<<std::setprecision(5)<<ThetaV 490 << std::setw(12)<<std::setprecision(5)<<PhiV 491 << std::setw(12)<<std::setprecision(5)<<qOnp 492 << std::setw(12)<<std::setprecision(5)<<1./qOnp 493 <<std::endl; 494 *pOutCurrent 495 << "* +/- " 496 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[0][0]) 497 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[1][1]) 498 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[2][2]) 499 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[3][3]) 500 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[4][4]) 501 <<std::endl; 502 *pOutCurrent 503 << "* Covariance matrice " 504 <<std::endl; 505 *pOutCurrent 506 << "* : " 507 << std::setw(12)<<std::setprecision(2)<<covVert[0][0]*100./sqrt(covVert[0][0]*covVert[0][0]) 508 << std::setw(12)<<std::setprecision(2)<<covVert[0][1]*100./sqrt(covVert[0][0]*covVert[1][1]) 509 << std::setw(12)<<std::setprecision(2)<<covVert[0][2]*100./sqrt(covVert[0][0]*covVert[2][2]) 510 << std::setw(12)<<std::setprecision(2)<<covVert[0][3]*100./sqrt(covVert[0][0]*covVert[3][3]) 511 << std::setw(12)<<std::setprecision(2)<<covVert[0][4]*100./sqrt(covVert[0][0]*covVert[4][4]) 512 <<std::endl; 513 for (int i1=1; i1<=4; i1++){ 514 *pOutCurrent 515 << "* " 516 << std::setw(12)<<std::setprecision(2)<<covVert[i1][0]*100./sqrt(covVert[i1][i1]*covVert[0][0]) 517 << std::setw(12)<<std::setprecision(2)<<covVert[i1][1]*100./sqrt(covVert[i1][i1]*covVert[1][1]) 518 << std::setw(12)<<std::setprecision(2)<<covVert[i1][2]*100./sqrt(covVert[i1][i1]*covVert[2][2]) 519 << std::setw(12)<<std::setprecision(2)<<covVert[i1][3]*100./sqrt(covVert[i1][i1]*covVert[3][3]) 520 << std::setw(12)<<std::setprecision(2)<<covVert[i1][4]*100./sqrt(covVert[i1][i1]*covVert[4][4]) 521 <<std::endl; 522 } 523 } 524 *pOutCurrent 525 << "* position " 526 << std::setw(12)<<std::setprecision(3)<<(pAtaCylinder->position()).x() 527 <<std::endl; 528 *pOutCurrent 529 << "* " 530 << std::setw(12)<<std::setprecision(3)<<(pAtaCylinder->position()).y() 531 <<std::endl; 532 *pOutCurrent 533 << "* " 534 << std::setw(12)<<std::setprecision(3)<<(pAtaCylinder->position()).z() 535 <<std::endl; 536 *pOutCurrent 537 << "* momentum " 538 << std::setw(12)<<std::setprecision(3)<<(pAtaCylinder->momentum()).x() 539 <<std::endl; 540 *pOutCurrent 541 << "* " 542 << std::setw(12)<<std::setprecision(3)<<(pAtaCylinder->momentum()).y() 543 <<std::endl; 544 *pOutCurrent 545 << "* " 546 << std::setw(12)<<std::setprecision(3)<<(pAtaCylinder->momentum()).z() 547 <<std::endl; 548 549 return StatusCode::SUCCESS; 550 551 } 552 StatusCode StacoEDMHelper::dump_AtaDisc(const Trk::AtaDisc* pAtaDisc,std::ofstream* pOutCurrent){ 553 554 if (pAtaDisc == 0) return StatusCode::SUCCESS; 555 *pOutCurrent << "* Its parameters at AtaDisc are:"<<std::endl; 556 if (IsThisPositionBeyondThelimitCaloSpectro( 557 (pAtaDisc->position()).x(), 558 (pAtaDisc->position()).y(), 559 (pAtaDisc->position()).z()) == 1 ) { 560 *pOutCurrent << "* At Muon Spectro Entrance ?"<<std::endl; 561 }else{ 562 *pOutCurrent << "* At Calo Entrance ?"<<std::endl; 563 } 564 565 *pOutCurrent << "*parameters()[Trk::locR] " <<std::setw(12)<<std::setprecision(5)<< pAtaDisc->parameters()[Trk::locR] << std::endl; 566 *pOutCurrent << " parameters()[Trk::locPhi] " << std::setw(12)<<std::setprecision(5)<<pAtaDisc->parameters()[Trk::locPhi] << std::endl; 567 *pOutCurrent << " parameters()[Trk::phi] " << std::setw(12)<<std::setprecision(5)<<pAtaDisc->parameters()[Trk::phi] << std::endl; 568 *pOutCurrent << " parameters()[Trk::theta] " << std::setw(12)<<std::setprecision(5)<<pAtaDisc->parameters()[Trk::theta] << std::endl; 569 *pOutCurrent << " parameters()[Trk::qOverP] " << std::setw(12)<<std::setprecision(5)<<pAtaDisc->parameters()[Trk::qOverP] << std::endl; 570 571 //double RCrossing = sqrt ( (pAtaDisc->position()).x()*(pAtaDisc->position()).x() + (pAtaDisc->position()).y()*(pAtaDisc->position()).y() ); 572 double ZCrossing = (pAtaDisc->position()).z() ; 573 double ThetaX = atan( pAtaDisc->parameters()[Trk::locR]/ZCrossing ) ; 574 if ( ThetaX < 0. ) ThetaX = ThetaX + 3.1415927 ; 575 double PhiX = pAtaDisc->parameters()[Trk::locPhi] ; 576 double ThetaV = pAtaDisc->parameters()[Trk::theta] ; 577 double PhiV = pAtaDisc->parameters()[Trk::phi] ; 578 double qOnp = pAtaDisc->parameters()[Trk::qOverP] ; 579 580 const Trk::MeasuredAtaDisc* pMeasuredAtaDisc = dynamic_cast<const Trk::MeasuredAtaDisc*>(pAtaDisc); 581 if (pMeasuredAtaDisc != 0){ 582 HepMatrix TheJac(5,5,0); 583 TheJac(1,1) = std::pow( cos(ThetaX) , 2 )/ZCrossing; 584 TheJac(2,2) = 1. ; 585 TheJac(3,4) = 1. ; 586 TheJac(4,3) = 1. ; 587 TheJac(5,5) = 1. ; 588 Trk::CovarianceMatrix covVert(5); 589 covVert = pMeasuredAtaDisc->localErrorMatrix().covariance().similarity(TheJac); 590 *pOutCurrent 591 << "* at AtaDisc : " 592 << std::setw(12)<<std::setprecision(5)<<ThetaX 593 << std::setw(12)<<std::setprecision(5)<<PhiX 594 << std::setw(12)<<std::setprecision(5)<<ThetaV 595 << std::setw(12)<<std::setprecision(5)<<PhiV 596 << std::setw(12)<<std::setprecision(5)<<qOnp 597 << std::setw(12)<<std::setprecision(5)<<1./qOnp 598 <<std::endl; 599 *pOutCurrent 600 << "* +/- " 601 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[0][0]) 602 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[1][1]) 603 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[2][2]) 604 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[3][3]) 605 << std::setw(12)<<std::setprecision(5)<<sqrt(covVert[4][4]) 606 <<std::endl; 607 *pOutCurrent 608 << "* Covariance matrice " 609 <<std::endl; 610 *pOutCurrent 611 << "* : " 612 << std::setw(12)<<std::setprecision(2)<<covVert[0][0]*100./sqrt(covVert[0][0]*covVert[0][0]) 613 << std::setw(12)<<std::setprecision(2)<<covVert[0][1]*100./sqrt(covVert[0][0]*covVert[1][1]) 614 << std::setw(12)<<std::setprecision(2)<<covVert[0][2]*100./sqrt(covVert[0][0]*covVert[2][2]) 615 << std::setw(12)<<std::setprecision(2)<<covVert[0][3]*100./sqrt(covVert[0][0]*covVert[3][3]) 616 << std::setw(12)<<std::setprecision(2)<<covVert[0][4]*100./sqrt(covVert[0][0]*covVert[4][4]) 617 <<std::endl; 618 for (int i1=1; i1<=4; i1++){ 619 *pOutCurrent 620 << "* " 621 << std::setw(12)<<std::setprecision(2)<<covVert[i1][0]*100./sqrt(covVert[i1][i1]*covVert[0][0]) 622 << std::setw(12)<<std::setprecision(2)<<covVert[i1][1]*100./sqrt(covVert[i1][i1]*covVert[1][1]) 623 << std::setw(12)<<std::setprecision(2)<<covVert[i1][2]*100./sqrt(covVert[i1][i1]*covVert[2][2]) 624 << std::setw(12)<<std::setprecision(2)<<covVert[i1][3]*100./sqrt(covVert[i1][i1]*covVert[3][3]) 625 << std::setw(12)<<std::setprecision(2)<<covVert[i1][4]*100./sqrt(covVert[i1][i1]*covVert[4][4]) 626 <<std::endl; 627 } 628 } 629 *pOutCurrent 630 << "* position " 631 << std::setw(12)<<std::setprecision(3)<<(pAtaDisc->position()).x() 632 <<std::endl; 633 *pOutCurrent 634 << "* " 635 << std::setw(12)<<std::setprecision(3)<<(pAtaDisc->position()).y() 636 <<std::endl; 637 *pOutCurrent 638 << "* " 639 << std::setw(12)<<std::setprecision(3)<<(pAtaDisc->position()).z() 640 <<std::endl; 641 *pOutCurrent 642 << "* momentum " 643 << std::setw(12)<<std::setprecision(3)<<(pAtaDisc->momentum()).x() 644 <<std::endl; 645 *pOutCurrent 646 << "* " 647 << std::setw(12)<<std::setprecision(3)<<(pAtaDisc->momentum()).y() 648 <<std::endl; 649 *pOutCurrent 650 << "* " 651 << std::setw(12)<<std::setprecision(3)<<(pAtaDisc->momentum()).z() 652 <<std::endl; 653 654 return StatusCode::SUCCESS; 655 656 } 657 658 // Dumps TrkSegment 659 StatusCode StacoEDMHelper::dump_TrkSegment( 660 const Trk::Segment* pSegment, 661 std::ofstream* pOutCurrent){ 662 663 if ( IsValid(pSegment) == 0 ) return StatusCode::SUCCESS; 664 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 665 666 //Get ErrorMatrix 667 //Change of interface of MuonSegment: JFL Fri Jan 19 10:00:33 CET 2007 668 //const Trk::ErrorMatrix TheErrorMatrix = pMuonSegment->fullLocalErrorMatrix() ; 669 const Trk::ErrorMatrix TheErrorMatrix = pMuonSegment->localErrorMatrix(); 670 const Trk::WeightMatrix TheWeightMatrix = TheErrorMatrix.weight(); 671 const Trk::CovarianceMatrix TheCovarianceMatrix = TheErrorMatrix.covariance(); 672 673 // *pOutCurrent 674 // << "Station Index/Name " 675 // << " " << SegmentAssociatedStationIndex(pSegment) 676 // << " " << SegmentAssociatedStationName(pSegment) 677 // <<std::endl; 678 679 int Ihas2ndCoord = 0 ; 680 if (SegmentHas2ndCoordMeasurement(pSegment)) Ihas2ndCoord = 1 ; 681 *pOutCurrent 682 << "Has2ndCoord " 683 << std::setw(20)<< Ihas2ndCoord 684 << std::endl; 685 686 double X_Xaxis,Y_Xaxis,Z_Xaxis; 687 GetXaxis(pSegment,X_Xaxis,Y_Xaxis,Z_Xaxis); 688 *pOutCurrent 689 << "Xaxis " 690 << " " << std::setw(12)<<std::setprecision(5) << X_Xaxis 691 << " " << std::setw(12)<<std::setprecision(5) << Y_Xaxis 692 << " " << std::setw(12)<<std::setprecision(5) << Z_Xaxis 693 <<std::endl; 694 695 double X_Yaxis,Y_Yaxis,Z_Yaxis; 696 GetYaxis(pSegment,X_Yaxis,Y_Yaxis,Z_Yaxis); 697 *pOutCurrent 698 << "Yaxis " 699 << " " << std::setw(12)<<std::setprecision(5) << X_Yaxis 700 << " " << std::setw(12)<<std::setprecision(5) << Y_Yaxis 701 << " " << std::setw(12)<<std::setprecision(5) << Z_Yaxis 702 <<std::endl; 703 704 double ThePositionX,ThePositionY,ThePositionZ; 705 GetPoint(pSegment,ThePositionX,ThePositionY,ThePositionZ); 706 *pOutCurrent 707 << "Position " 708 << " " << std::setw(12)<<std::setprecision(5) << ThePositionX 709 << " " << std::setw(12)<<std::setprecision(5) << ThePositionY 710 << " " << std::setw(12)<<std::setprecision(5) << ThePositionZ 711 <<std::endl; 712 713 double TheDirectionX,TheDirectionY,TheDirectionZ; 714 GetDirection(pSegment,TheDirectionX,TheDirectionY,TheDirectionZ); 715 *pOutCurrent 716 << "Direction " 717 << " " << std::setw(12)<<std::setprecision(5) << TheDirectionX 718 << " " << std::setw(12)<<std::setprecision(5) << TheDirectionY 719 << " " << std::setw(12)<<std::setprecision(5) << TheDirectionZ 720 <<std::endl; 721 722 *pOutCurrent 723 << "Errors " 724 <<std::endl; 725 for (int Ind1=0; Ind1<4; Ind1++){ 726 *pOutCurrent 727 << " " << std::setw(12)<<std::setprecision(5) << TheWeightMatrix[Ind1][0] 728 << " " << std::setw(12)<<std::setprecision(5) << TheWeightMatrix[Ind1][1] 729 << " " << std::setw(12)<<std::setprecision(5) << TheWeightMatrix[Ind1][2] 730 << " " << std::setw(12)<<std::setprecision(5) << TheWeightMatrix[Ind1][3] 731 <<std::endl; 732 } 733 *pOutCurrent 734 << "dx Axz mm/mrad " 735 << std::setw(12)<<std::setprecision(4)<< 1./sqrt(TheWeightMatrix[0][0]) 736 << std::setw(12)<<std::setprecision(4)<< 1000./sqrt(TheWeightMatrix[1][1]) 737 << std::setw(12)<<std::setprecision(4)<< TheWeightMatrix[0][1]/(sqrt(TheWeightMatrix[0][0]*TheWeightMatrix[1][1])) 738 << std::endl; 739 *pOutCurrent 740 << "dy Axz micron/microrad " 741 << std::setw(12)<<std::setprecision(4)<< 1000./sqrt(TheWeightMatrix[2][2]) 742 << std::setw(12)<<std::setprecision(4)<< 1000000./sqrt(TheWeightMatrix[3][3]) 743 << std::setw(12)<<std::setprecision(4)<< TheWeightMatrix[2][3]/(sqrt(TheWeightMatrix[2][2]*TheWeightMatrix[3][3])) 744 << std::endl; 745 746 *pOutCurrent 747 << "CErrors " 748 <<std::endl; 749 for (int Ind1=0; Ind1<4; Ind1++){ 750 *pOutCurrent 751 << " " << std::setw(12)<<std::setprecision(5) << TheCovarianceMatrix[Ind1][0] 752 << " " << std::setw(12)<<std::setprecision(5) << TheCovarianceMatrix[Ind1][1] 753 << " " << std::setw(12)<<std::setprecision(5) << TheCovarianceMatrix[Ind1][2] 754 << " " << std::setw(12)<<std::setprecision(5) << TheCovarianceMatrix[Ind1][3] 755 <<std::endl; 756 } 757 *pOutCurrent 758 << "Cdx Axz mm/mrad " 759 << std::setw(12)<<std::setprecision(4)<< sqrt(TheCovarianceMatrix[0][0]) 760 << std::setw(12)<<std::setprecision(4)<< 1000.*sqrt(TheCovarianceMatrix[1][1]) 761 << std::setw(12)<<std::setprecision(4)<< TheCovarianceMatrix[0][1]/(sqrt(TheCovarianceMatrix[0][0]*TheCovarianceMatrix[1][1])) 762 << std::endl; 763 *pOutCurrent 764 << "Cdy Axz micron/microrad" 765 << std::setw(12)<<std::setprecision(4)<< 1000.*sqrt(TheCovarianceMatrix[2][2]) 766 << std::setw(12)<<std::setprecision(4)<< 1000000.*sqrt(TheCovarianceMatrix[3][3]) 767 << std::setw(12)<<std::setprecision(4)<< TheCovarianceMatrix[2][3]/(sqrt(TheCovarianceMatrix[2][2]*TheCovarianceMatrix[3][3])) 768 << std::endl; 769 770 return StatusCode::SUCCESS; 771 772 } 773 774 // Dumps Rots of Trk Objects 775 StatusCode StacoEDMHelper::dump_TrkTrackDigits( 776 const Trk::Track* pTrack, 777 std::ofstream* pOutCurrent){ 778 779 if ( pTrack == 0 ) return StatusCode::SUCCESS; 780 if ( pTrack->measurementsOnTrack() == 0 ) return StatusCode::SUCCESS; 781 782 std::vector<const Trk::MeasurementBase*>::const_iterator imEnd = pTrack->measurementsOnTrack()->end(); 783 784 std::vector<const Trk::MeasurementBase*>::const_iterator im = pTrack->measurementsOnTrack()->begin(); 785 for ( ; im!=imEnd; ++im){ 786 const Trk::RIO_OnTrack* pRIO_OnTrack = dynamic_cast<const Trk::RIO_OnTrack*>(*im); 787 if ( pRIO_OnTrack!=0 ){ 788 const Identifier& channelId = pRIO_OnTrack->identify(); 789 if (m_detID->is_muon(channelId)){ 790 StatusCode sc = dump_rot(pRIO_OnTrack,pOutCurrent); 791 if (sc.isFailure()) msg(MSG::WARNING) << "dump_rot failed" << endreq; 792 } 793 }else{ 794 msg(MSG::WARNING) <<"Null pointer to ROT found on Track " <<endreq; 795 } 796 } 797 798 im = pTrack->measurementsOnTrack()->begin(); 799 for ( ; im!=imEnd; ++im){ 800 const Trk::RIO_OnTrack* pRIO_OnTrack = dynamic_cast<const Trk::RIO_OnTrack*>(*im); 801 if ( pRIO_OnTrack!=0 ){ 802 const Identifier& channelId = pRIO_OnTrack->identify(); 803 if (!m_detID->is_muon(channelId)) *pOutCurrent << *pRIO_OnTrack << std::endl; 804 }else{ 805 msg(MSG::WARNING) <<"Null pointer to ROT found on Track " <<endreq; 806 } 807 } 808 809 return StatusCode::SUCCESS; 810 811 } 812 StatusCode StacoEDMHelper::dump_TrkSegmentDigits( 813 const Trk::Segment* pSegment, 814 std::ofstream* pOutCurrent){ 815 816 if ( IsValid(pSegment) == 0 ) return StatusCode::SUCCESS; 817 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 818 819 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 820 821 std::vector<const Trk::RIO_OnTrack*>::const_iterator pRIO_OnTrackIter = pRIOSet.begin(); 822 for (; pRIO_OnTrackIter!=pRIOSet.end(); ++pRIO_OnTrackIter){ 823 if ( (*pRIO_OnTrackIter)!=0 ){ 824 const Identifier& channelId = (*pRIO_OnTrackIter)->identify(); 825 if (m_detID->is_muon(channelId)){ 826 StatusCode sc= dump_rot((*pRIO_OnTrackIter),pOutCurrent); 827 if (sc.isFailure()) msg(MSG::WARNING) << "dump_rot failed" << endreq; 828 } 829 }else{ 830 msg(MSG::WARNING) <<"Null pointer to ROT found in pRIOSet " <<endreq; 831 } 832 } 833 834 pRIO_OnTrackIter = pRIOSet.begin(); 835 for (; pRIO_OnTrackIter!=pRIOSet.end(); ++pRIO_OnTrackIter){ 836 if ( (*pRIO_OnTrackIter)!=0 ){ 837 const Identifier& channelId = (*pRIO_OnTrackIter)->identify(); 838 if (!m_detID->is_muon(channelId)) *pOutCurrent << *(*pRIO_OnTrackIter) << std::endl; 839 }else{ 840 msg(MSG::WARNING) <<"Null pointer to ROT found in pRIOSet " <<endreq; 841 } 842 } 843 844 return StatusCode::SUCCESS; 845 846 } 847 StatusCode StacoEDMHelper::dump_TrkSegmentHoles( 848 const Trk::Segment* pSegment, 849 std::ofstream* pOutCurrent){ 850 851 if ( IsValid(pSegment) == 0 ) return StatusCode::SUCCESS; 852 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 853 854 const Muon::MuonSegmentQuality* pMuonSegmentQuality = dynamic_cast<const Muon::MuonSegmentQuality*>(pMuonSegment->fitQuality()); 855 if( pMuonSegmentQuality ){ 856 if (pMuonSegmentQuality->channelsWithoutHit().size() != 0 ){ 857 std::vector<Identifier>::const_iterator IdCurrent = pMuonSegmentQuality->channelsWithoutHit().begin(); 858 std::vector<Identifier>::const_iterator IdCurrent_end = pMuonSegmentQuality->channelsWithoutHit().end(); 859 for( ;IdCurrent!=IdCurrent_end;++IdCurrent ){ 860 const Identifier& channelId = *IdCurrent; 861 StatusCode sc = dump_id(channelId,pOutCurrent); 862 if (sc.isFailure()) msg(MSG::WARNING) << "dump_id failed" << endreq; 863 } 864 }else{ 865 *pOutCurrent << "This Segment does not contain holes " << std::endl; 866 } 867 }else{ 868 *pOutCurrent << "This Segment does not contain a pMuonSegmentQuality" << std::endl; 869 } 870 871 return StatusCode::SUCCESS; 872 873 } 874 875 // Dumps Rot 876 StatusCode StacoEDMHelper::dump_rot(const Trk::RIO_OnTrack* rot,std::ofstream* pOutCurrent){ 877 878 if ( rot == 0 ) return StatusCode::SUCCESS; 879 880 const Identifier& channelId = rot->identify(); 881 StatusCode sc = dump_id(channelId,pOutCurrent); 882 *pOutCurrent << *rot << std::endl; 883 if (sc.isFailure()) msg(MSG::WARNING) << "dump_id failed" << endreq; 884 885 return StatusCode::SUCCESS; 886 887 } 888 // Dumps Ids 889 StatusCode StacoEDMHelper::dump_id(Identifier channelId,std::ofstream* pOutCurrent){ 890 891 if (m_detID->is_muon(channelId)){ 892 if(m_mdtId->is_mdt(channelId)){ 893 StatusCode sc = dump_mdt(channelId,pOutCurrent); 894 if (sc.isFailure()) msg(MSG::WARNING) << "dump_mdt failed" << endreq; 895 } 896 if(m_rpcId->is_rpc(channelId)){ 897 StatusCode sc = dump_rpc(channelId,pOutCurrent); 898 if (sc.isFailure()) msg(MSG::WARNING) << "dump_rpc failed" << endreq; 899 } 900 if(m_tgcId->is_tgc(channelId)){ 901 StatusCode sc = dump_tgc(channelId,pOutCurrent); 902 if (sc.isFailure()) msg(MSG::WARNING) << "dump_tgc failed" << endreq; 903 } 904 if(m_cscId->is_csc(channelId)){ 905 StatusCode sc = dump_csc(channelId,pOutCurrent); 906 if (sc.isFailure()) msg(MSG::WARNING) << "dump_csc failed" << endreq; 907 } 908 } 909 910 return StatusCode::SUCCESS; 911 912 } 913 StatusCode StacoEDMHelper::dump_mdt(Identifier channelId,std::ofstream* pOutCurrent){ 914 915 if(m_mdtId->is_mdt(channelId)){ 916 int SGStationNber = m_mdtId->stationName(channelId); 917 std::string SGStationName = m_mdtId->stationNameString(SGStationNber); 918 int SGStationEta = m_mdtId->stationEta(channelId); 919 int SGStationPhi = m_mdtId->stationPhi(channelId); 920 int SGMultilayer = m_mdtId->multilayer(channelId) ; 921 int SGTubeLayer = m_mdtId->tubeLayer(channelId) ; 922 int SGTube = m_mdtId->tube(channelId) ; 923 *pOutCurrent 924 << "MDT " 925 << " " 926 << std::setw(3)<< SGStationName 927 << std::setw(4)<< SGStationEta 928 << std::setw(4)<< SGStationPhi 929 << std::setw(4)<< SGMultilayer 930 << std::setw(4)<< SGTubeLayer 931 << std::setw(4)<< SGTube 932 << std::endl; 933 } 934 935 return StatusCode::SUCCESS; 936 937 } 938 StatusCode StacoEDMHelper::dump_rpc(Identifier channelId,std::ofstream* pOutCurrent){ 939 940 if(m_rpcId->is_rpc(channelId)){ 941 int SGStationNber = m_rpcId->stationName(channelId); 942 std::string SGStationName = m_rpcId->stationNameString(SGStationNber); 943 int SGStationEta = m_rpcId->stationEta(channelId); 944 int SGStationPhi = m_rpcId->stationPhi(channelId); 945 int SGDoubletR = m_rpcId->doubletR(channelId) ; 946 int SGDoubletZ = m_rpcId->doubletZ(channelId) ; 947 int SGDoubletPhi = m_rpcId->doubletPhi(channelId) ; 948 int SGGasGap = m_rpcId->gasGap(channelId) ; 949 int SGMeasuresPhi = m_rpcId->measuresPhi(channelId) ; 950 int SGStrip = m_rpcId->strip(channelId) ; 951 *pOutCurrent 952 << "RPC " 953 << " " 954 << std::setw(3)<< SGStationName 955 << std::setw(4)<< SGStationEta 956 << std::setw(4)<< SGStationPhi 957 << std::setw(4)<< SGDoubletR 958 << std::setw(4)<< SGDoubletZ 959 << std::setw(4)<< SGDoubletPhi 960 << std::setw(4)<< SGGasGap 961 << std::setw(4)<< SGMeasuresPhi 962 << std::setw(4)<< SGStrip 963 << std::endl; 964 } 965 966 return StatusCode::SUCCESS; 967 968 } 969 StatusCode StacoEDMHelper::dump_tgc(Identifier channelId,std::ofstream* pOutCurrent){ 970 971 if(m_tgcId->is_tgc(channelId)){ 972 int SGStationNber = m_tgcId->stationName(channelId); 973 std::string SGStationName = m_tgcId->stationNameString(SGStationNber); 974 int SGStationEta = m_tgcId->stationEta(channelId); 975 int SGStationPhi = m_tgcId->stationPhi(channelId); 976 int SGGasGap = m_tgcId->gasGap(channelId) ; 977 int SGIsStrip = m_tgcId->isStrip(channelId) ; 978 int SGChannel = m_tgcId->channel(channelId) ; 979 *pOutCurrent 980 << "TGC " 981 << " " 982 << std::setw(3)<< SGStationName 983 << std::setw(4)<< SGStationEta 984 << std::setw(4)<< SGStationPhi 985 << std::setw(4)<< SGGasGap 986 << std::setw(4)<< SGIsStrip 987 << std::setw(4)<< SGChannel 988 << std::endl; 989 } 990 991 return StatusCode::SUCCESS; 992 993 } 994 StatusCode StacoEDMHelper::dump_csc(Identifier channelId,std::ofstream* pOutCurrent){ 995 996 if(m_cscId->is_csc(channelId)){ 997 int SGStationNber = m_cscId->stationName(channelId); 998 std::string SGStationName = m_cscId->stationNameString(SGStationNber); 999 int SGStationEta = m_cscId->stationEta(channelId); 1000 int SGStationPhi = m_cscId->stationPhi(channelId); 1001 int SGChamberLayer = m_cscId->chamberLayer(channelId) ; 1002 int SGWireLayer = m_cscId->wireLayer(channelId) ; 1003 int SGMeasuresPhi = m_cscId->measuresPhi(channelId) ; 1004 int SGStrip = m_cscId->strip(channelId) ; 1005 *pOutCurrent 1006 << "CSC " 1007 << " " 1008 << std::setw(3)<< SGStationName 1009 << std::setw(4)<< SGStationEta 1010 << std::setw(4)<< SGStationPhi 1011 << std::setw(4)<< SGChamberLayer 1012 << std::setw(4)<< SGWireLayer 1013 << std::setw(4)<< SGMeasuresPhi 1014 << std::setw(4)<< SGStrip 1015 << std::endl; 1016 } 1017 1018 return StatusCode::SUCCESS; 1019 1020 } 1021 // Compute Eloss 1022 StatusCode StacoEDMHelper::BackTrackingEloss( 1023 const Rec::TrackParticle* pTrackParticle, 1024 double& Eloss){ 1025 1026 Eloss = 0. ; 1027 1028 if ( pTrackParticle == 0 ) return( StatusCode::FAILURE ); 1029 1030 StatusCode scEloss = BackTrackingEloss( (pTrackParticle->originalTrack()) , Eloss) ; 1031 if ( scEloss.isFailure() ) return( StatusCode::FAILURE ); 1032 1033 return StatusCode::SUCCESS; 1034 1035 } 1036 // StatusCode StacoEDMHelper::BackTrackingEloss(const Trk::Track* pTrack, double& Eloss){ 1037 // 1038 // 1039 // Eloss = 0. ; 1040 // 1041 // if ( pTrack == 0) return StatusCode::FAILURE; 1042 // if ( pTrack->trackParameters() == 0) return StatusCode::FAILURE; 1043 // 1044 // std::vector<const Trk::TrackParameters*>::const_iterator TrkPar = pTrack->trackParameters()->begin(); 1045 // std::vector<const Trk::TrackParameters*>::const_iterator TrkParEnd = pTrack->trackParameters()->end(); 1046 // 1047 // //Find the MS entrance track 1048 // const Trk::AtaCylinder* pAtaCylinder_MS = 0 ; 1049 // const Trk::AtaDisc* pAtaDisc_MS = 0 ; 1050 // TrkPar = pTrack->trackParameters()->begin(); 1051 // for ( ; TrkPar!=TrkParEnd; ++TrkPar){ 1052 // pAtaCylinder_MS = dynamic_cast<const Trk::AtaCylinder*>(*TrkPar); 1053 // if ( pAtaCylinder_MS == 0 ) continue ; 1054 // if (IsThisPositionBeyondThelimitCaloSpectro( 1055 // (pAtaCylinder_MS->position()).x(), 1056 // (pAtaCylinder_MS->position()).y(), 1057 // (pAtaCylinder_MS->position()).z()) != 1 ) pAtaCylinder_MS = 0 ; 1058 // if ( pAtaCylinder_MS !=0 ) break; 1059 // } 1060 // if ( pAtaCylinder_MS == 0 ){ 1061 // TrkPar = pTrack->trackParameters()->begin(); 1062 // for ( ; TrkPar!=TrkParEnd; ++TrkPar){ 1063 // pAtaDisc_MS = dynamic_cast<const Trk::AtaDisc*>(*TrkPar); 1064 // if ( pAtaDisc_MS == 0 ) continue ; 1065 // if (IsThisPositionBeyondThelimitCaloSpectro( 1066 // (pAtaDisc_MS->position()).x(), 1067 // (pAtaDisc_MS->position()).y(), 1068 // (pAtaDisc_MS->position()).z()) != 1 ) pAtaDisc_MS = 0 ; 1069 // if ( pAtaDisc_MS !=0 ) break; 1070 // } 1071 // } 1072 // if ( pAtaCylinder_MS == 0 && pAtaDisc_MS == 0) return StatusCode::FAILURE; 1073 // 1074 // double E_MS = 0. ; 1075 // if ( pAtaCylinder_MS != 0 ) { 1076 // E_MS = sqrt ( (pAtaCylinder_MS->momentum()).x()*(pAtaCylinder_MS->momentum()).x() 1077 // + (pAtaCylinder_MS->momentum()).y()*(pAtaCylinder_MS->momentum()).y() 1078 // + (pAtaCylinder_MS->momentum()).z()*(pAtaCylinder_MS->momentum()).z() ) ; 1079 // }else{ 1080 // E_MS = sqrt ( (pAtaDisc_MS->momentum()).x()*(pAtaDisc_MS->momentum()).x() 1081 // + (pAtaDisc_MS->momentum()).y()*(pAtaDisc_MS->momentum()).y() 1082 // + (pAtaDisc_MS->momentum()).z()*(pAtaDisc_MS->momentum()).z() ) ; 1083 // } 1084 // 1085 // //Find the Kalo entrance track 1086 // const Trk::AtaCylinder* pAtaCylinder_KA = 0 ; 1087 // const Trk::AtaDisc* pAtaDisc_KA = 0 ; 1088 // TrkPar = pTrack->trackParameters()->begin(); 1089 // for ( ; TrkPar!=TrkParEnd; ++TrkPar){ 1090 // pAtaCylinder_KA = dynamic_cast<const Trk::AtaCylinder*>(*TrkPar); 1091 // if ( pAtaCylinder_KA == 0 ) continue ; 1092 // if (IsThisPositionBeyondThelimitCaloSpectro( 1093 // (pAtaCylinder_KA->position()).x(), 1094 // (pAtaCylinder_KA->position()).y(), 1095 // (pAtaCylinder_KA->position()).z()) == 1 ) pAtaCylinder_KA = 0 ; 1096 // if ( pAtaCylinder_KA !=0 ) break; 1097 // } 1098 // if ( pAtaCylinder_KA == 0 ){ 1099 // TrkPar = pTrack->trackParameters()->begin(); 1100 // for ( ; TrkPar!=TrkParEnd; ++TrkPar){ 1101 // pAtaDisc_KA = dynamic_cast<const Trk::AtaDisc*>(*TrkPar); 1102 // if ( pAtaDisc_KA == 0 ) continue ; 1103 // if (IsThisPositionBeyondThelimitCaloSpectro( 1104 // (pAtaDisc_KA->position()).x(), 1105 // (pAtaDisc_KA->position()).y(), 1106 // (pAtaDisc_KA->position()).z()) == 1 ) pAtaDisc_KA = 0 ; 1107 // if ( pAtaDisc_KA !=0 ) break; 1108 // } 1109 // } 1110 // if ( pAtaCylinder_KA == 0 && pAtaDisc_KA == 0) return StatusCode::FAILURE; 1111 // 1112 // double E_KA = 0. ; 1113 // if ( pAtaCylinder_KA != 0 ) { 1114 // E_KA = sqrt ( (pAtaCylinder_KA->momentum()).x()*(pAtaCylinder_KA->momentum()).x() 1115 // + (pAtaCylinder_KA->momentum()).y()*(pAtaCylinder_KA->momentum()).y() 1116 // + (pAtaCylinder_KA->momentum()).z()*(pAtaCylinder_KA->momentum()).z() ) ; 1117 // }else{ 1118 // E_KA = sqrt ( (pAtaDisc_KA->momentum()).x()*(pAtaDisc_KA->momentum()).x() 1119 // + (pAtaDisc_KA->momentum()).y()*(pAtaDisc_KA->momentum()).y() 1120 // + (pAtaDisc_KA->momentum()).z()*(pAtaDisc_KA->momentum()).z() ) ; 1121 // } 1122 // 1123 // Eloss = E_KA - E_MS ; 1124 // 1125 // return StatusCode::SUCCESS; 1126 // 1127 // } 1128 StatusCode StacoEDMHelper::BackTrackingEloss(const Trk::Track* pTrack, double& Eloss){ 1129 1130 Eloss = 0. ; 1131 1132 if ( pTrack == 0) return StatusCode::FAILURE; 1133 if ( pTrack->trackParameters() == 0) return StatusCode::FAILURE; 1134 1135 std::vector<const Trk::TrackParameters*>::const_iterator TrkPar = pTrack->trackParameters()->begin(); 1136 std::vector<const Trk::TrackParameters*>::const_iterator TrkParEnd = pTrack->trackParameters()->end(); 1137 1138 double ToCaloMax = 10000000. ; 1139 double EatCalo = 0. ; 1140 double ToSpectroMax = 10000000. ; 1141 double EatSpectro = 0. ; 1142 int FoundCalo = 0 ; 1143 int FoundSpectro = 0 ; 1144 TrkPar = pTrack->trackParameters()->begin(); 1145 for ( ; TrkPar!=TrkParEnd; ++TrkPar){ 1146 const Trk::MeasuredTrackParameters* pMeasuredTrackParameters = dynamic_cast<const Trk::MeasuredTrackParameters*>(*TrkPar); 1147 if ( pMeasuredTrackParameters ) { 1148 if (IsThisPositionBeyondThelimitCaloSpectro( 1149 ((*TrkPar)->position()).x(), 1150 ((*TrkPar)->position()).y(), 1151 ((*TrkPar)->position()).z()) != 1 ) { 1152 FoundCalo = 1 ; 1153 double ToCalo = DistanceToCaloCylinder( 1154 ((*TrkPar)->position()).x(), 1155 ((*TrkPar)->position()).y(), 1156 ((*TrkPar)->position()).z() 1157 ); 1158 if (ToCalo < ToCaloMax ){ 1159 ToCaloMax = ToCalo ; 1160 EatCalo = sqrt ( ((*TrkPar)->momentum()).x()*((*TrkPar)->momentum()).x() 1161 + ((*TrkPar)->momentum()).y()*((*TrkPar)->momentum()).y() 1162 + ((*TrkPar)->momentum()).z()*((*TrkPar)->momentum()).z() ) ; 1163 } 1164 }else{ 1165 FoundSpectro = 1 ; 1166 double ToSpectro = DistanceToSpectroCylinder( 1167 ((*TrkPar)->position()).x(), 1168 ((*TrkPar)->position()).y(), 1169 ((*TrkPar)->position()).z() 1170 ); 1171 if (ToSpectro < ToSpectroMax ){ 1172 ToSpectroMax = ToSpectro ; 1173 EatSpectro = sqrt ( ((*TrkPar)->momentum()).x()*((*TrkPar)->momentum()).x() 1174 + ((*TrkPar)->momentum()).y()*((*TrkPar)->momentum()).y() 1175 + ((*TrkPar)->momentum()).z()*((*TrkPar)->momentum()).z() ) ; 1176 } 1177 } 1178 } 1179 } 1180 1181 if (( FoundCalo !=1 ) ||( FoundSpectro !=1 )) 1182 return StatusCode::FAILURE; 1183 1184 Eloss = EatCalo - EatSpectro ; 1185 1186 return StatusCode::SUCCESS; 1187 1188 } 1189 double StacoEDMHelper::DistanceToCaloCylinder(double x, double y, double z){ 1190 1191 double TheR = 1050. ; 1192 double TheHalfZ = 3200. ; 1193 return DistanceToCylinder(TheHalfZ,TheR,x,y,z) ; 1194 1195 } 1196 double StacoEDMHelper::DistanceToSpectroCylinder(double x, double y, double z){ 1197 1198 double TheR = 4250. ; 1199 double TheHalfZ = 6820. ; 1200 return DistanceToCylinder(TheHalfZ,TheR,x,y,z) ; 1201 1202 } 1203 double StacoEDMHelper::DistanceToCylinder( double TheHalfZ, double TheR, double x, double y, double z){ 1204 1205 double DeltaR = fabs( sqrt ( x*x + y*y ) - TheR ) ; 1206 double DeltaZ = fabs( fabs( z ) - TheHalfZ ) ; 1207 1208 double ToBeReturned = DeltaR ; 1209 if ( DeltaZ < ToBeReturned ) ToBeReturned = DeltaZ ; 1210 1211 return ToBeReturned ; 1212 1213 } 1214 int StacoEDMHelper::IsThisPositionBeyondThelimitCaloSpectro(double x, double y, double z){ 1215 1216 double RLimitCaloSpectro = 3000. ; 1217 double ZLimitCaloSpectro = 5000. ; 1218 1219 double R = sqrt ( x*x + y*y ) ; 1220 if (R >= RLimitCaloSpectro ) return 1 ; 1221 if (fabs(z) >= fabs(ZLimitCaloSpectro) ) return 1; 1222 return 0 ; 1223 1224 } 1225 1226 // Compute Eloss 1227 StatusCode StacoEDMHelper::BackTrackingEloss( 1228 const Rec::TrackParticle* muonSpecTP, 1229 const Rec::TrackParticle* extrapTP, 1230 double& Eloss){ 1231 1232 Eloss = 0. ; 1233 1234 if ( muonSpecTP == 0 ) return( StatusCode::FAILURE ); 1235 if ( extrapTP == 0 ) return( StatusCode::FAILURE ); 1236 1237 StatusCode scEloss = BackTrackingEloss( (muonSpecTP->originalTrack()), (extrapTP->originalTrack()) ,Eloss) ; 1238 if ( scEloss.isFailure() ) return( StatusCode::FAILURE ); 1239 1240 return StatusCode::SUCCESS; 1241 1242 } 1243 StatusCode StacoEDMHelper::BackTrackingEloss( 1244 const Trk::Track* muonSpecTT, 1245 const Trk::Track* extrapTT, 1246 double& Eloss){ 1247 1248 Eloss = 0. ; 1249 1250 if ( muonSpecTT == 0 ) return( StatusCode::FAILURE ); 1251 if ( extrapTT == 0 ) return( StatusCode::FAILURE ); 1252 1253 double E_muonSpecTT = sqrt ( 1254 ((muonSpecTT->perigeeParameters())->momentum()).x()*((muonSpecTT->perigeeParameters())->momentum()).x() 1255 + ((muonSpecTT->perigeeParameters())->momentum()).y()*((muonSpecTT->perigeeParameters())->momentum()).y() 1256 + ((muonSpecTT->perigeeParameters())->momentum()).z()*((muonSpecTT->perigeeParameters())->momentum()).z() ) ; 1257 1258 double E_extrapTT = sqrt ( 1259 ((extrapTT->perigeeParameters())->momentum()).x()*((extrapTT->perigeeParameters())->momentum()).x() 1260 + ((extrapTT->perigeeParameters())->momentum()).y()*((extrapTT->perigeeParameters())->momentum()).y() 1261 + ((extrapTT->perigeeParameters())->momentum()).z()*((extrapTT->perigeeParameters())->momentum()).z() ) ; 1262 1263 Eloss = E_extrapTT - E_muonSpecTT ; 1264 1265 return StatusCode::SUCCESS; 1266 1267 } 1268 1269 // Segment Associated Station Name 1270 int StacoEDMHelper::SegmentAssociatedStationIndex(const Trk::Segment* pSegment){ 1271 1272 int ToBeReturned = -1 ; 1273 1274 std::string TheStationName = SegmentAssociatedStationName(pSegment); 1275 if (TheStationName != "XXX" ) ToBeReturned = m_mdtId->stationNameIndex(TheStationName); 1276 1277 return ToBeReturned ; 1278 1279 } 1280 1281 std::string StacoEDMHelper::SegmentAssociatedStationName(const Trk::Segment* pSegment){ 1282 1283 std::string ToBeReturned = "XXX" ; 1284 1285 if ( IsValid(pSegment) == 0 ) return ToBeReturned; 1286 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1287 1288 //JFL 06/09/16 Let the segment with no RDO be named CSC 1289 //JFL 06/09/16 NB: no station has this name 1290 ToBeReturned = "CSC" ; 1291 1292 std::vector<int> NberOfRot ; 1293 std::vector<std::string> NameOfStation ; 1294 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1295 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1296 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1297 if ( (*rotIter)!=0 ){ 1298 const Identifier& id = (*rotIter)->identify(); 1299 if (m_detID->is_muon(id)){ 1300 std::string SGStationName ; 1301 int ThereIsARotToSort = 0 ; 1302 if(m_mdtId->is_mdt(id)){ 1303 SGStationName = m_mdtId->stationNameString(m_mdtId->stationName(id)); 1304 ThereIsARotToSort = 1 ; 1305 } 1306 if(m_cscId->is_csc(id)){ 1307 if(m_cscId->measuresPhi(id) == 0 ){ 1308 SGStationName = m_cscId->stationNameString(m_cscId->stationName(id)); 1309 ThereIsARotToSort = 1 ; 1310 } 1311 } 1312 if (ThereIsARotToSort == 1){ 1313 int NameOfStationSize = NameOfStation.size(); 1314 int IsFound = 0 ; 1315 for (int ItemNber=0; ItemNber<NameOfStationSize ; ItemNber++){ 1316 if (IsFound == 0 ){ 1317 if ( NameOfStation[ItemNber] == SGStationName ){ 1318 IsFound = 1 ; 1319 NberOfRot[ItemNber] = NberOfRot[ItemNber] + 1; 1320 } 1321 } 1322 } 1323 if ( IsFound == 0 ){ 1324 NameOfStation.push_back(SGStationName) ; 1325 int IntToAdd = 1 ; 1326 NberOfRot.push_back(IntToAdd) ; 1327 } 1328 } 1329 } 1330 } 1331 } 1332 1333 int IsFound = 0 ; 1334 std::string NameSelected = "XXX" ; 1335 int NameOfStationSize = NameOfStation.size(); 1336 int NberOfRotMax = -1 ; 1337 for (int ItemNber=0; ItemNber<NameOfStationSize ; ItemNber++){ 1338 if (NberOfRot[ItemNber] > NberOfRotMax){ 1339 NberOfRotMax = NberOfRot[ItemNber] ; 1340 NameSelected = NameOfStation[ItemNber] ; 1341 IsFound = 1 ; 1342 } 1343 } 1344 if ( IsFound == 1 ) ToBeReturned = NameSelected ; 1345 1346 return ToBeReturned ; 1347 1348 } 1349 1350 int StacoEDMHelper::StationLevel( const std::string& TheSegmentStationName ) 1351 { 1352 char StationNameFragment[1]; 1353 StationNameFragment[0]= TheSegmentStationName[1]; 1354 if ( StationNameFragment[0] == 'O' ) return 2; 1355 else if ( StationNameFragment[0] == 'M' ) return 1; 1356 else if ( StationNameFragment[0] == 'I' ) return 0; 1357 return -1; 1358 } 1359 1360 bool StacoEDMHelper::IsEndCap( const std::string& TheSegmentStationName ) 1361 { 1362 char StationNameFragment[1]; 1363 StationNameFragment[0]= TheSegmentStationName[0]; 1364 if ( StationNameFragment[0] == 'B' ) return 0; 1365 else return 1; 1366 } 1367 1368 bool StacoEDMHelper::IsLargeSector( const std::string& TheSegmentStationName ) 1369 { 1370 char StationNameFragment[1]; 1371 StationNameFragment[0]= TheSegmentStationName[2]; 1372 if ( StationNameFragment[0] == 'S' ) return 0; 1373 else return 1; 1374 } 1375 1376 // 1377 int StacoEDMHelper::GetNberOfPrecisionDigitsInThisTrack( 1378 const Trk::Segment* pSegment, 1379 const Trk::Track* pTrack){ 1380 1381 int ToBeReturned = 0 ; 1382 1383 if ( IsValid(pSegment) == 0 ) return ToBeReturned; 1384 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1385 1386 if ( pTrack == 0 ) return 0 ; 1387 if ( pTrack->measurementsOnTrack() == 0 ) return 0 ; 1388 1389 int NberOfRotShared = 0; 1390 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1391 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1392 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1393 if ( (*rotIter)!=0 ){ 1394 int Is1stCoord = 0 ; 1395 if ( m_mdtId->is_mdt( (*rotIter)->identify() ) ) Is1stCoord = 1; 1396 if ( m_cscId->is_csc( (*rotIter)->identify() ) ){ 1397 if(m_cscId->measuresPhi( (*rotIter)->identify() ) == 0 ) Is1stCoord = 1; 1398 } 1399 if ( Is1stCoord == 1 ){ 1400 if (IsThisRotInThisTrack( (*rotIter) , pTrack ) == 1 ) NberOfRotShared = NberOfRotShared + 1; 1401 } 1402 } 1403 } 1404 1405 ToBeReturned = NberOfRotShared ; 1406 1407 return ToBeReturned; 1408 1409 } 1410 1411 // Is this Segment in this Track? 1412 int StacoEDMHelper::IsThisSegmentInThisTrack( 1413 const Trk::Segment* pSegment, 1414 const Trk::Track* pTrack){ 1415 1416 int ToBeReturned = 0 ; 1417 1418 if ( IsValid(pSegment) == 0 ) return ToBeReturned; 1419 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1420 1421 if ( pTrack == 0 ) return 0 ; 1422 if ( pTrack->measurementsOnTrack() == 0 ) return 0 ; 1423 1424 int NberOfRot = 0; 1425 int NberOfRotShared = 0; 1426 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1427 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1428 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1429 if ( (*rotIter)!=0 ){ 1430 int Is1stCoord = 0 ; 1431 if ( m_mdtId->is_mdt( (*rotIter)->identify() ) ) Is1stCoord = 1; 1432 if ( m_cscId->is_csc( (*rotIter)->identify() ) ){ 1433 if(m_cscId->measuresPhi( (*rotIter)->identify() ) == 0 ) Is1stCoord = 1; 1434 } 1435 if ( Is1stCoord == 1 ){ 1436 NberOfRot = NberOfRot + 1; 1437 if (IsThisRotInThisTrack( (*rotIter) , pTrack ) == 1 ) NberOfRotShared = NberOfRotShared + 1; 1438 } 1439 } 1440 } 1441 int NberOfRotNotShared = NberOfRot - NberOfRotShared ; 1442 if (NberOfRotNotShared <= 2) ToBeReturned = 1 ; 1443 if (NberOfRot == 0) ToBeReturned = 0 ; 1444 1445 return ToBeReturned; 1446 1447 } 1448 1449 // Is this Rot in this Track? 1450 int StacoEDMHelper::IsThisRotInThisTrack( 1451 const Trk::RIO_OnTrack* rot, 1452 const Trk::Track* pTrack){ 1453 1454 if ( rot == 0 ) return 0 ; 1455 if ( pTrack == 0 ) return 0 ; 1456 if ( pTrack->measurementsOnTrack() == 0 ) return 0 ; 1457 1458 int ToBeReturned = 0 ; 1459 std::vector<const Trk::MeasurementBase*>::const_iterator im = pTrack->measurementsOnTrack()->begin(); 1460 std::vector<const Trk::MeasurementBase*>::const_iterator imEnd = pTrack->measurementsOnTrack()->end(); 1461 for ( ; im!=imEnd; ++im){ 1462 const Trk::RIO_OnTrack* RotTrack = dynamic_cast<const Trk::RIO_OnTrack*>(*im); 1463 if ( RotTrack != 0 ){ 1464 int Test = AreTheseRotsIdentical(rot,RotTrack) ; 1465 if ( Test == 1 ) { 1466 ToBeReturned = 1 ; 1467 break ; 1468 } 1469 } 1470 } 1471 1472 return ToBeReturned; 1473 1474 } 1475 1476 // Are these Rots Indentical? 1477 int StacoEDMHelper::AreTheseRotsIdentical( 1478 const Trk::RIO_OnTrack* Rot1, 1479 const Trk::RIO_OnTrack* Rot2){ 1480 1481 if ( Rot1 == 0 ) return 0 ; 1482 if ( Rot2 == 0 ) return 0 ; 1483 1484 const Identifier& Id1 = Rot1->identify() ; 1485 const Identifier& Id2 = Rot2->identify() ; 1486 if ( Id1 == Id2 ) return 1 ; 1487 1488 return 0; 1489 1490 } 1491 1492 // Reconstructed Objects Validity 1493 int StacoEDMHelper::IsValid( 1494 const Trk::Segment* pSegment){ 1495 1496 int ToBeReturned = 0 ; 1497 1498 if ( pSegment == 0 ) return ToBeReturned ; 1499 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1500 if ( pMuonSegment == 0 ) return ToBeReturned; 1501 1502 ToBeReturned = 1 ; 1503 1504 return ToBeReturned ; 1505 1506 } 1507 1508 // Utilities for Segment 1509 Hep3Vector StacoEDMHelper::GetXaxis(const Trk::Segment* pSegment) 1510 { 1511 double X,Y,Z; 1512 GetXaxis(pSegment,X,Y,Z); 1513 Hep3Vector v(X,Y,Z); 1514 return v; 1515 } 1516 void StacoEDMHelper::GetXaxis( 1517 const Trk::Segment* pSegment, 1518 double& X_Xaxis, 1519 double& Y_Xaxis, 1520 double& Z_Xaxis 1521 ){ 1522 1523 X_Xaxis = 0. ; 1524 Y_Xaxis = 0. ; 1525 Z_Xaxis = 0. ; 1526 1527 if ( IsValid(pSegment) != 0 ){ 1528 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1529 const Trk::PlaneSurface ThePlaneSurface = pMuonSegment->associatedSurface() ; 1530 double dx = 1.; 1531 double dy = 0.; 1532 double dz = 0.; 1533 HepVector3D LocXaxis(dx,dy,dz); 1534 const Trk::GlobalPosition GlobXaxis( (ThePlaneSurface.transform())*LocXaxis ); 1535 X_Xaxis = GlobXaxis.x(); 1536 Y_Xaxis = GlobXaxis.y(); 1537 Z_Xaxis = GlobXaxis.z(); 1538 } 1539 1540 } 1541 Hep3Vector StacoEDMHelper::GetYaxis(const Trk::Segment* pSegment) 1542 { 1543 double X,Y,Z; 1544 GetYaxis(pSegment,X,Y,Z); 1545 Hep3Vector v(X,Y,Z); 1546 return v; 1547 } 1548 void StacoEDMHelper::GetYaxis( 1549 const Trk::Segment* pSegment, 1550 double& X_Yaxis, 1551 double& Y_Yaxis, 1552 double& Z_Yaxis 1553 ){ 1554 1555 X_Yaxis = 0. ; 1556 Y_Yaxis = 0. ; 1557 Z_Yaxis = 0. ; 1558 1559 if ( IsValid(pSegment) != 0 ){ 1560 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1561 const Trk::PlaneSurface ThePlaneSurface = pMuonSegment->associatedSurface() ; 1562 double dx = 0.; 1563 double dy = 1.; 1564 double dz = 0.; 1565 HepVector3D LocYaxis(dx,dy,dz); 1566 const Trk::GlobalPosition GlobYaxis( (ThePlaneSurface.transform())*LocYaxis ); 1567 X_Yaxis = GlobYaxis.x(); 1568 Y_Yaxis = GlobYaxis.y(); 1569 Z_Yaxis = GlobYaxis.z(); 1570 } 1571 1572 } 1573 Hep3Vector StacoEDMHelper::GetPoint(const Trk::Segment* pSegment) 1574 { 1575 double X,Y,Z; 1576 GetPoint(pSegment,X,Y,Z); 1577 Hep3Vector v(X,Y,Z); 1578 return v; 1579 } 1580 void StacoEDMHelper::GetPoint( 1581 const Trk::Segment* pSegment, 1582 double& Xpt, 1583 double& Ypt, 1584 double& Zpt 1585 ){ 1586 1587 Xpt = 0. ; 1588 Ypt = 0. ; 1589 Zpt = 0. ; 1590 1591 if ( IsValid(pSegment) != 0 ){ 1592 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1593 const Trk::GlobalPosition TheGlobalPosition = pMuonSegment->globalPosition() ; 1594 Xpt = TheGlobalPosition.x(); 1595 Ypt = TheGlobalPosition.y(); 1596 Zpt = TheGlobalPosition.z(); 1597 } 1598 1599 } 1600 Hep3Vector StacoEDMHelper::GetDirection(const Trk::Segment* pSegment) 1601 { 1602 double X,Y,Z; 1603 GetDirection(pSegment,X,Y,Z); 1604 Hep3Vector v(X,Y,Z); 1605 return v; 1606 } 1607 void StacoEDMHelper::GetDirection( 1608 const Trk::Segment* pSegment, 1609 double& Xdir, 1610 double& Ydir, 1611 double& Zdir 1612 ){ 1613 1614 Xdir = 0. ; 1615 Ydir = 0. ; 1616 Zdir = 0. ; 1617 1618 if ( IsValid(pSegment) != 0 ){ 1619 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1620 const Trk::GlobalDirection TheGlobalDirection = pMuonSegment->globalDirection() ; 1621 Xdir = TheGlobalDirection.x(); 1622 Ydir = TheGlobalDirection.y(); 1623 Zdir = TheGlobalDirection.z(); 1624 } 1625 1626 } 1627 double StacoEDMHelper::GetfitQuality( 1628 const Trk::Segment* pSegment 1629 ){ 1630 1631 double ToBeReturned = 9999. ; 1632 1633 if ( IsValid(pSegment) != 0 ){ 1634 const Trk::FitQuality* pFitQuality = pSegment->fitQuality() ; 1635 ToBeReturned = pFitQuality->chiSquared () ; 1636 } 1637 1638 return ToBeReturned ; 1639 1640 } 1641 // nd = total number of digits 1642 // ndMDT1 = number of digits on 1 MDT multilayer 1643 // ndMDT2 = number of digits on the other MDT multilayer 1644 void StacoEDMHelper::GetNberOfDigi( 1645 const Trk::Segment* pSegment, 1646 int& nd, int& ndMDT1, int& ndMDT2 1647 ){ 1648 1649 nd = 0 ; 1650 ndMDT1 = 0 ; 1651 ndMDT2 = 0 ; 1652 1653 if ( IsValid(pSegment) != 0 ){ 1654 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1655 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1656 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1657 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1658 if ( (*rotIter)!=0 ){ 1659 const Identifier& id = (*rotIter)->identify(); 1660 if (m_detID->is_muon(id)){ 1661 if( m_mdtId->is_mdt(id) 1662 || m_rpcId->is_rpc(id) 1663 || m_tgcId->is_tgc(id) 1664 || m_cscId->is_csc(id)) 1665 { 1666 ++nd; 1667 if(m_mdtId->is_mdt(id)) 1668 { 1669 int idML = m_mdtId->multilayer(id); 1670 if(1==idML)++ndMDT1; 1671 else ++ndMDT2; 1672 } 1673 } 1674 } 1675 } 1676 } 1677 } 1678 1679 } 1680 1681 int StacoEDMHelper::GetNberOfDigi( 1682 const Trk::Segment* pSegment 1683 ){ 1684 1685 int ToBeReturned = 0 ; 1686 1687 if ( IsValid(pSegment) != 0 ){ 1688 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1689 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1690 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1691 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1692 if ( (*rotIter)!=0 ){ 1693 const Identifier& id = (*rotIter)->identify(); 1694 if (m_detID->is_muon(id)){ 1695 if( m_mdtId->is_mdt(id) 1696 || m_rpcId->is_rpc(id) 1697 || m_tgcId->is_tgc(id) 1698 || m_cscId->is_csc(id)) 1699 ToBeReturned = ToBeReturned + 1 ; 1700 } 1701 } 1702 } 1703 } 1704 1705 return ToBeReturned ; 1706 1707 } 1708 1709 int StacoEDMHelper::GetNberOfDigiMDT( 1710 const Trk::Segment* pSegment 1711 ){ 1712 1713 int ToBeReturned = 0 ; 1714 1715 if ( IsValid(pSegment) != 0 ){ 1716 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1717 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1718 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1719 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1720 if ( (*rotIter)!=0 ){ 1721 const Identifier& id = (*rotIter)->identify(); 1722 if (m_detID->is_muon(id)){ 1723 if( m_mdtId->is_mdt(id)) 1724 ToBeReturned = ToBeReturned + 1 ; 1725 } 1726 } 1727 } 1728 } 1729 1730 return ToBeReturned ; 1731 1732 } 1733 1734 int StacoEDMHelper::GetNberOfDigiRPC( 1735 const Trk::Segment* pSegment 1736 ){ 1737 1738 int ToBeReturned = 0 ; 1739 1740 if ( IsValid(pSegment) != 0 ){ 1741 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1742 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1743 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1744 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1745 if ( (*rotIter)!=0 ){ 1746 const Identifier& id = (*rotIter)->identify(); 1747 if (m_detID->is_muon(id)){ 1748 if( m_rpcId->is_rpc(id)) 1749 ToBeReturned = ToBeReturned + 1 ; 1750 } 1751 } 1752 } 1753 } 1754 1755 return ToBeReturned ; 1756 1757 } 1758 1759 int StacoEDMHelper::GetNberOfDigiTGC( 1760 const Trk::Segment* pSegment 1761 ){ 1762 1763 int ToBeReturned = 0 ; 1764 1765 if ( IsValid(pSegment) != 0 ){ 1766 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1767 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1768 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1769 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1770 if ( (*rotIter)!=0 ){ 1771 const Identifier& id = (*rotIter)->identify(); 1772 if (m_detID->is_muon(id)){ 1773 if( m_tgcId->is_tgc(id)) 1774 ToBeReturned = ToBeReturned + 1 ; 1775 } 1776 } 1777 } 1778 } 1779 1780 return ToBeReturned ; 1781 1782 } 1783 1784 int StacoEDMHelper::GetNberOfDigiCSC( 1785 const Trk::Segment* pSegment 1786 ){ 1787 1788 int ToBeReturned = 0 ; 1789 1790 if ( IsValid(pSegment) != 0 ){ 1791 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1792 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1793 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1794 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1795 if ( (*rotIter)!=0 ){ 1796 const Identifier& id = (*rotIter)->identify(); 1797 if (m_detID->is_muon(id)){ 1798 if( m_cscId->is_csc(id)){ 1799 ToBeReturned = ToBeReturned + 1 ; 1800 } 1801 } 1802 } 1803 } 1804 } 1805 1806 return ToBeReturned ; 1807 1808 } 1809 1810 int StacoEDMHelper::GetNberOfDigiCSC1( 1811 const Trk::Segment* pSegment 1812 ){ 1813 1814 int ToBeReturned = 0 ; 1815 1816 if ( IsValid(pSegment) != 0 ){ 1817 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1818 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1819 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1820 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1821 if ( (*rotIter)!=0 ){ 1822 const Identifier& id = (*rotIter)->identify(); 1823 if (m_detID->is_muon(id)){ 1824 if( m_cscId->is_csc(id)){ 1825 if( !m_cscId->measuresPhi(id)) 1826 ToBeReturned = ToBeReturned + 1 ; 1827 } 1828 } 1829 } 1830 } 1831 } 1832 1833 return ToBeReturned ; 1834 1835 } 1836 1837 int StacoEDMHelper::GetNberOfDigiCSC2( 1838 const Trk::Segment* pSegment 1839 ){ 1840 1841 int ToBeReturned = 0 ; 1842 1843 if ( IsValid(pSegment) != 0 ){ 1844 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1845 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1846 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1847 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1848 if ( (*rotIter)!=0 ){ 1849 const Identifier& id = (*rotIter)->identify(); 1850 if (m_detID->is_muon(id)){ 1851 if( m_cscId->is_csc(id)){ 1852 if( m_cscId->measuresPhi(id)) 1853 ToBeReturned = ToBeReturned + 1 ; 1854 } 1855 } 1856 } 1857 } 1858 } 1859 1860 return ToBeReturned ; 1861 1862 } 1863 1864 int StacoEDMHelper::GetNberOfPrecisionDigitsInCommon( 1865 const Trk::Segment* pSegment, 1866 const Trk::Segment* pSegment2 1867 ){ 1868 1869 int ToBeReturned = 0 ; 1870 1871 if ( IsValid(pSegment) != 0 ){ 1872 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1873 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1874 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1875 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1876 if ( (*rotIter)!=0 ){ 1877 const Identifier& id = (*rotIter)->identify(); 1878 if (m_detID->is_muon(id)){ 1879 if( m_mdtId->is_mdt(id) || (m_cscId->is_csc(id) && m_cscId->measuresPhi(id)) ){ 1880 // 1881 if ( IsValid(pSegment2) != 0 ){ 1882 const Muon::MuonSegment* pMuonSegment2 = dynamic_cast<const Muon::MuonSegment*>(pSegment2); 1883 const std::vector<const Trk::RIO_OnTrack*> pRIOSet2 = pMuonSegment2->containedROTs(); 1884 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter2 = pRIOSet2.begin(); 1885 for (; rotIter2!=pRIOSet2.end(); ++rotIter2){ 1886 if ( (*rotIter2)!=0 ){ 1887 const Identifier& id2 = (*rotIter2)->identify(); 1888 if (id == id2){ 1889 ToBeReturned = ToBeReturned + 1 ; 1890 } 1891 } 1892 } 1893 } 1894 } 1895 } 1896 } 1897 } 1898 } 1899 1900 return ToBeReturned ; 1901 1902 } 1903 1904 void StacoEDMHelper::GetNberOfDigiVec( 1905 const Trk::Segment* pSegment, 1906 std::vector<int>& etaDigitVec, 1907 std::vector<int>& phiDigitVec 1908 ){ 1909 int nmdt(0); 1910 int nrpcEta(0); 1911 int ntgcEta(0); 1912 int ncscEta(0); 1913 int nrpcPhi(0); 1914 int ntgcPhi(0); 1915 int ncscPhi(0); 1916 if ( IsValid(pSegment) != 0 ){ 1917 const Muon::MuonSegment* pMuonSegment = dynamic_cast<const Muon::MuonSegment*>(pSegment); 1918 const std::vector<const Trk::RIO_OnTrack*> pRIOSet = pMuonSegment->containedROTs(); 1919 std::vector<const Trk::RIO_OnTrack*>::const_iterator rotIter = pRIOSet.begin(); 1920 for (; rotIter!=pRIOSet.end(); ++rotIter){ 1921 if ( (*rotIter)!=0 ){ 1922 const Identifier& id = (*rotIter)->identify(); 1923 if (m_detID->is_muon(id)) { 1924 if( m_mdtId->is_mdt(id) ) ++nmdt; 1925 if( m_rpcId->is_rpc(id) ) 1926 { 1927 if(m_rpcId->measuresPhi(id))++nrpcPhi; 1928 else ++nrpcEta; 1929 } 1930 if( m_tgcId->is_tgc(id) ) 1931 { 1932 if(m_tgcId->isStrip(id))++ntgcPhi; 1933 else ++ntgcEta; 1934 } 1935 if( m_cscId->is_csc(id) ) 1936 { 1937 if(m_cscId->measuresPhi(id))++ncscPhi; 1938 else ++ncscEta; 1939 } 1940 } 1941 } 1942 } 1943 } 1944 etaDigitVec.clear(); 1945 etaDigitVec.push_back(nmdt); // MDTs provide Eta measurements only 1946 etaDigitVec.push_back(ncscEta); 1947 etaDigitVec.push_back(nrpcEta); 1948 etaDigitVec.push_back(ntgcEta); 1949 phiDigitVec.clear(); 1950 phiDigitVec.push_back(0); 1951 phiDigitVec.push_back(ncscPhi); 1952 phiDigitVec.