Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-31 09:21:54

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /*
0027  * =============================================================================
0028  *
0029  *       Filename:  CexmcSetup.cc
0030  *
0031  *    Description:  physical setup
0032  *
0033  *        Version:  1.0
0034  *        Created:  10.10.2009 23:00:50
0035  *       Revision:  none
0036  *       Compiler:  gcc
0037  *
0038  *         Author:  Alexey Radkov (), 
0039  *        Company:  PNPI
0040  *
0041  * =============================================================================
0042  */
0043 
0044 #include <G4GDMLParser.hh>
0045 #include <G4MultiFunctionalDetector.hh>
0046 #include <G4SDManager.hh>
0047 #include <G4LogicalVolume.hh>
0048 #include <G4VPhysicalVolume.hh>
0049 #include <G4PhysicalVolumeStore.hh>
0050 #include <G4Box.hh>
0051 #include <G4LogicalVolumeStore.hh>
0052 #include <G4Region.hh>
0053 #include <G4RegionStore.hh>
0054 #include <G4ProductionCuts.hh>
0055 #include <G4VUserPhysicsList.hh>
0056 #include <G4SystemOfUnits.hh>
0057 #include "CexmcSetup.hh"
0058 #include "CexmcPrimitiveScorer.hh"
0059 #include "CexmcTrackPoints.hh"
0060 #include "CexmcTrackPointsInLeftRightSet.hh"
0061 #include "CexmcTrackPointsInCalorimeter.hh"
0062 #include "CexmcTrackPointsFilter.hh"
0063 #include "CexmcSimpleEnergyDeposit.hh"
0064 #include "CexmcEnergyDepositInLeftRightSet.hh"
0065 #include "CexmcEnergyDepositInCalorimeter.hh"
0066 #include "CexmcRunManager.hh"
0067 #include "CexmcPhysicsManager.hh"
0068 #include "CexmcException.hh"
0069 
0070 
0071 CexmcSetup::CexmcSetup( const G4String &  gdmlFile_,
0072                         G4bool  validateGDMLFile_ ) :
0073     world( 0 ), gdmlFile( gdmlFile_ ), validateGDMLFile( validateGDMLFile_ ),
0074     calorimeterRegionInitialized( false ),
0075     calorimeterGeometryDataInitialized( false ), monitorVolume( NULL ),
0076     vetoCounterVolume( NULL ), calorimeterVolume( NULL ), targetVolume( NULL ),
0077     rightVetoCounter( NULL ), rightCalorimeter( NULL )
0078 {
0079 }
0080 
0081 
0082 G4VPhysicalVolume *  CexmcSetup::Construct( void )
0083 {
0084     if ( world )
0085         return world;
0086 
0087     G4GDMLParser  gdmlParser;
0088 
0089     gdmlParser.Read( gdmlFile, validateGDMLFile );
0090     world = gdmlParser.GetWorldVolume();
0091 
0092     SetupSpecialVolumes( gdmlParser );
0093 
0094     ReadTransforms( gdmlParser );
0095 
0096     ReadRightDetectors();
0097 
0098     CexmcRunManager *  runManager( static_cast< CexmcRunManager * >(
0099                                             G4RunManager::GetRunManager() ) );
0100 
0101     runManager->SetupConstructionHook();
0102 
0103     const CexmcPhysicsManager *  physicsManager(
0104             dynamic_cast< const CexmcPhysicsManager * >(
0105                                         runManager->GetUserPhysicsList() ) );
0106 
0107     if ( ! physicsManager )
0108         throw CexmcException( CexmcWeirdException );
0109 
0110     CexmcPhysicsManager *        thePhysicsManager(
0111             const_cast< CexmcPhysicsManager * >( physicsManager ) );
0112     thePhysicsManager->SetupConstructionHook( this );
0113 
0114     return world;
0115 }
0116 
0117 
0118 void  CexmcSetup::SetupSpecialVolumes( const G4GDMLParser &  gdmlParser )
0119 {
0120     G4MultiFunctionalDetector *   detector[ CexmcNumberOfDetectorRoles ] =
0121                                                                     { NULL };
0122     const G4LogicalVolumeStore *  lvs( G4LogicalVolumeStore::GetInstance() );
0123 
0124     for ( std::vector< G4LogicalVolume * >::const_iterator
0125                         lvIter( lvs->begin() ); lvIter != lvs->end(); ++lvIter )
0126     {
0127         G4String           volumeName( G4String( ( *lvIter )->GetName() ) );
0128         G4GDMLAuxListType  auxInfo( gdmlParser.GetVolumeAuxiliaryInformation(
0129                                                                     *lvIter ) );
0130         CexmcDetectorRole  curDetectorRole( CexmcNumberOfDetectorRoles );
0131 
0132         for ( G4GDMLAuxListType::const_iterator  pair( auxInfo.begin() );
0133                                               pair != auxInfo.end(); ++pair )
0134         {
0135             CexmcPrimitiveScorer *  scorer( NULL );
0136             G4String                detectorName( "uninitialized" );
0137             do
0138             {
0139                 if ( pair->type == "EnergyDepositDetector" )
0140                 {
0141                     do
0142                     {
0143                         if ( pair->value == "MonitorRole" )
0144                         {
0145                             AssertAndAsignDetectorRole( curDetectorRole,
0146                                                 CexmcMonitorDetectorRole );
0147                             scorer = new CexmcSimpleEnergyDeposit(
0148                                     CexmcDetectorTypeName[ CexmcEDDetector ] );
0149                             break;
0150                         }
0151                         if ( pair->value == "VetoCounterRole" )
0152                         {
0153                             AssertAndAsignDetectorRole( curDetectorRole,
0154                                                 CexmcVetoCounterDetectorRole );
0155                             scorer = new CexmcEnergyDepositInLeftRightSet(
0156                                     CexmcDetectorTypeName[ CexmcEDDetector ],
0157                                     this );
0158                             break;
0159                         }
0160                         if ( pair->value == "CalorimeterRole" )
0161                         {
0162                             AssertAndAsignDetectorRole( curDetectorRole,
0163                                                 CexmcCalorimeterDetectorRole );
0164                             scorer = new CexmcEnergyDepositInCalorimeter(
0165                                     CexmcDetectorTypeName[ CexmcEDDetector ],
0166                                     this );
0167                             break;
0168                         }
0169                     } while ( false );
0170                     detectorName = CexmcDetectorRoleName[ curDetectorRole ];
0171                     G4cout << CEXMC_LINE_START "ED Scorer of detector role '" <<
0172                                detectorName << "' in volume '" << volumeName <<
0173                                "'" << G4endl;
0174                     break;
0175                 }
0176                 if ( pair->type == "TrackPointsDetector" )
0177                 {
0178                     do
0179                     {
0180                         if ( pair->value == "MonitorRole" )
0181                         {
0182                             AssertAndAsignDetectorRole( curDetectorRole,
0183                                                 CexmcMonitorDetectorRole );
0184                             scorer = new CexmcTrackPoints(
0185                                     CexmcDetectorTypeName[ CexmcTPDetector ] );
0186                             break;
0187                         }
0188                         if ( pair->value == "VetoCounterRole"  )
0189                         {
0190                             AssertAndAsignDetectorRole( curDetectorRole,
0191                                                 CexmcVetoCounterDetectorRole );
0192                             scorer = new CexmcTrackPointsInLeftRightSet(
0193                                     CexmcDetectorTypeName[ CexmcTPDetector ],
0194                                     this );
0195                             break;
0196                         }
0197                         if ( pair->value == "CalorimeterRole" )
0198                         {
0199                             AssertAndAsignDetectorRole( curDetectorRole,
0200                                                 CexmcCalorimeterDetectorRole );
0201                             scorer = new CexmcTrackPointsInCalorimeter(
0202                                     CexmcDetectorTypeName[ CexmcTPDetector ],
0203                                     this );
0204                             break;
0205                         }
0206                         if ( pair->value == "TargetRole" )
0207                         {
0208                             AssertAndAsignDetectorRole( curDetectorRole,
0209                                                 CexmcTargetDetectorRole );
0210                             scorer = new CexmcTrackPoints(
0211                                     CexmcDetectorTypeName[ CexmcTPDetector ] );
0212                             break;
0213                         }
0214                     } while ( false );
0215                     detectorName = CexmcDetectorRoleName[ curDetectorRole ];
0216                     G4cout << CEXMC_LINE_START "TP Scorer of detector role '" <<
0217                                detectorName << "' in volume '" << volumeName <<
0218                                "'" << G4endl;
0219                     if ( scorer )
0220                     {
0221                         CexmcTrackPointsFilter *  filter(
0222                                  new CexmcTrackPointsFilter( "trackPoints" ) );
0223                         scorer->SetFilter( filter );
0224                     }
0225                     break;
0226                 }
0227                 if ( pair->type == "SensitiveRegion" )
0228                 {
0229                     do
0230                     {
0231                         if ( pair->value == "CalorimeterRegion" )
0232                         {
0233                             G4Region *  region( NULL );
0234                             if ( calorimeterRegionInitialized )
0235                             {
0236                                 region = G4RegionStore::GetInstance()->
0237                                         GetRegion( CexmcCalorimeterRegionName );
0238                             }
0239                             else
0240                             {
0241                                 region = new G4Region(
0242                                                 CexmcCalorimeterRegionName );
0243                                 G4ProductionCuts *  cuts(
0244                                                         new G4ProductionCuts );
0245                                 G4double  defaultProductionCut( 1.0 * mm );
0246                                 const G4VUserPhysicsList *  physicsList(
0247                                             G4RunManager::GetRunManager()->
0248                                                         GetUserPhysicsList() );
0249                                 if ( physicsList )
0250                                     defaultProductionCut =
0251                                             physicsList->GetDefaultCutValue();
0252                                 cuts->SetProductionCut( defaultProductionCut );
0253                                 region->SetProductionCuts( cuts );
0254                                 calorimeterRegionInitialized = true;
0255                             }
0256                             region->AddRootLogicalVolume( *lvIter );
0257                             break;
0258                         }
0259                     } while ( false );
0260                     G4cout << CEXMC_LINE_START "Sensitive Region for logical "
0261                                "volume '" << volumeName << "' registered" <<
0262                                G4endl;
0263                     break;
0264                 }
0265                 if ( pair->type == "SpecialVolume" )
0266                 {
0267                     do
0268                     {
0269                         if ( pair->value == "Monitor" )
0270                         {
0271                             monitorVolume = *lvIter;
0272                             G4cout << CEXMC_LINE_START "Monitor volume '";
0273                             break;
0274                         }
0275                         if ( pair->value == "VetoCounter" )
0276                         {
0277                             vetoCounterVolume = *lvIter;
0278                             G4cout << CEXMC_LINE_START "VetoCounter volume '";
0279                             break;
0280                         }
0281                         if ( pair->value == "Calorimeter" )
0282                         {
0283                             calorimeterVolume = *lvIter;
0284                             G4cout << CEXMC_LINE_START "Calorimeter volume '";
0285                             ReadCalorimeterGeometryData( *lvIter );
0286                             calorimeterGeometryDataInitialized = true;
0287                             break;
0288                         }
0289                         if ( pair->value == "Target" )
0290                         {
0291                             targetVolume = *lvIter;
0292                             G4cout << CEXMC_LINE_START "Target volume '";
0293                             break;
0294                         }
0295                     } while ( false );
0296                     G4cout << volumeName << "' registered" << G4endl;
0297                     break;
0298                 }
0299             }
0300             while ( false );
0301 
0302             if ( scorer )
0303             {
0304                 /* curDetectorRole must be intact when scorer is not NULL */
0305                 if ( ! detector[ curDetectorRole ] )
0306                 {
0307                     detector[ curDetectorRole ] =
0308                                 new G4MultiFunctionalDetector( detectorName );
0309                 }
0310                 detector[ curDetectorRole ]->RegisterPrimitive( scorer );
0311                 /* now that scorer has initialized pointer to its detector, its
0312                  * messenger's path shall be properly initialized as well */
0313                 scorer->InitializeMessenger();
0314                 /* NB: logical volumes in GDML file may not have multiple
0315                  * detector roles: for example volume Monitor may have only one
0316                  * role MonitorRole. This restriction arises from that fact that
0317                  * a logical volume may contain only one sensitive detector. */
0318                 ( *lvIter )->SetSensitiveDetector(
0319                                                 detector[ curDetectorRole ] );
0320             }
0321         }
0322     }
0323 
0324     if ( ! calorimeterRegionInitialized )
0325         throw CexmcException( CexmcCalorimeterRegionNotInitialized );
0326 
0327     if ( ! calorimeterGeometryDataInitialized )
0328         throw CexmcException( CexmcCalorimeterGeometryDataNotInitialized );
0329 
0330     for ( G4int  i( 0 ); i < CexmcNumberOfDetectorRoles; ++i )
0331     {
0332         if ( detector[ i ] )
0333             G4SDManager::GetSDMpointer()->AddNewDetector( detector[ i ] );
0334     }
0335 }
0336 
0337 
0338 void  CexmcSetup::ReadTransforms( const G4GDMLParser &  gdmlParser )
0339 {
0340     G4ThreeVector     position( gdmlParser.GetPosition( "TargetPos" ) );
0341     G4ThreeVector     rotation( gdmlParser.GetRotation( "TargetRot" ) );
0342     G4RotationMatrix  rm;
0343 
0344     RotateMatrix( rotation, rm );
0345     targetTransform.SetNetTranslation( position );
0346     targetTransform.SetNetRotation( rm );
0347 
0348     position = gdmlParser.GetPosition( "CalorimeterLeftPos" );
0349     rotation = gdmlParser.GetRotation( "CalorimeterLeftRot" );
0350     rm = G4RotationMatrix();
0351     RotateMatrix( rotation, rm );
0352     calorimeterLeftTransform.SetNetTranslation( position );
0353     calorimeterLeftTransform.SetNetRotation( rm );
0354 
0355     position = gdmlParser.GetPosition( "CalorimeterRightPos" );
0356     rotation = gdmlParser.GetRotation( "CalorimeterRightRot" );
0357     rm = G4RotationMatrix();
0358     RotateMatrix( rotation, rm );
0359     calorimeterRightTransform.SetNetTranslation( position );
0360     calorimeterRightTransform.SetNetRotation( rm );
0361 }
0362 
0363 
0364 void  CexmcSetup::ReadCalorimeterGeometryData(
0365                                             const G4LogicalVolume *  lVolume )
0366 {
0367     if ( lVolume->GetNoDaughters() == 0 )
0368         throw CexmcException( CexmcIncompatibleGeometry );
0369 
0370     G4VPhysicalVolume *  pVolume( lVolume->GetDaughter( 0 ) );
0371     EAxis                axis;
0372     G4double             width;
0373     G4double             offset;
0374     G4bool               consuming;
0375 
0376     if ( ! pVolume )
0377         throw CexmcException( CexmcIncompatibleGeometry );
0378 
0379     if ( pVolume->IsReplicated() )
0380     {
0381         pVolume->GetReplicationData( axis,
0382                                      calorimeterGeometry.nCrystalsInColumn,
0383                                      width, offset, consuming );
0384     }
0385 
0386     lVolume = pVolume->GetLogicalVolume();
0387 
0388     if ( lVolume->GetNoDaughters() == 0 )
0389         throw CexmcException( CexmcIncompatibleGeometry );
0390 
0391     pVolume = lVolume->GetDaughter( 0 );
0392 
0393     if ( ! pVolume )
0394         throw CexmcException( CexmcIncompatibleGeometry );
0395 
0396     if ( pVolume->IsReplicated() )
0397     {
0398         pVolume->GetReplicationData( axis, calorimeterGeometry.nCrystalsInRow,
0399                                      width, offset, consuming );
0400     }
0401 
0402     lVolume = pVolume->GetLogicalVolume();
0403 
0404     /* NB: this is not necessarily a crystal itself as far as crystals can be
0405      * wrapped in paper and other materials, but this is what reconstructor and
0406      * digitizers really need */
0407     G4Box *  crystalBox( dynamic_cast< G4Box * >( lVolume->GetSolid() ) );
0408 
0409     if ( ! crystalBox )
0410         throw CexmcException( CexmcIncompatibleGeometry );
0411 
0412     calorimeterGeometry.crystalWidth = crystalBox->GetXHalfLength() * 2;
0413     calorimeterGeometry.crystalHeight = crystalBox->GetYHalfLength() * 2;
0414     calorimeterGeometry.crystalLength = crystalBox->GetZHalfLength() * 2;
0415 }
0416 
0417 
0418 void  CexmcSetup::ConvertToCrystalGeometry( const G4ThreeVector &  src,
0419                     G4int &  row, G4int &  column, G4ThreeVector &  dst ) const
0420 {
0421     G4int     nCrystalsInColumn( calorimeterGeometry.nCrystalsInColumn );
0422     G4int     nCrystalsInRow( calorimeterGeometry.nCrystalsInRow );
0423     G4double  crystalWidth( calorimeterGeometry.crystalWidth );
0424     G4double  crystalHeight( calorimeterGeometry.crystalHeight );
0425 
0426     row = G4int( ( src.y() + crystalHeight * nCrystalsInColumn / 2 ) /
0427                  crystalHeight );
0428     column = G4int( ( src.x() + crystalWidth * nCrystalsInRow / 2 ) /
0429                     crystalWidth );
0430     G4double   xInCalorimeterOffset(
0431                     ( G4double( column ) - G4double( nCrystalsInRow ) / 2 ) *
0432                                             crystalWidth + crystalWidth / 2 );
0433     G4double   yInCalorimeterOffset(
0434                     ( G4double( row ) - G4double( nCrystalsInColumn ) / 2 ) *
0435                                             crystalHeight + crystalHeight / 2 );
0436     dst.setX( src.x() - xInCalorimeterOffset );
0437     dst.setY( src.y() - yInCalorimeterOffset );
0438 }
0439 
0440 
0441 void  CexmcSetup::ReadRightDetectors( void )
0442 {
0443     G4PhysicalVolumeStore *  pvs( G4PhysicalVolumeStore::GetInstance() );
0444 
0445     for ( std::vector< G4VPhysicalVolume * >::const_iterator  k( pvs->begin() );
0446                                                         k != pvs->end(); ++k )
0447     {
0448         /* FIXME: it would be more reasonable to find detectors from volumes
0449          * tagged with 'EnergyDepositDetector' or 'TrackPointsDetector', and not
0450          * from volumes tagged with 'SpecialVolume' as it is done here. However
0451          * in case of calorimeters the role of detectors are played by crystal
0452          * volumes, not the calorimeters themselves, but only calorimeters can
0453          * hold information about their left or right positions! Thus, following
0454          * considerations of convenience, right detectors for all detector roles
0455          * are chosen from volumes tagged with 'SpecialVolume' */
0456         do
0457         {
0458             if ( ( *k )->GetLogicalVolume() == vetoCounterVolume )
0459             {
0460                 if ( G4StrUtil::contains(( *k )->GetName(), "Right" ) )
0461                     rightVetoCounter = *k;
0462                 break;
0463             }
0464             if ( ( *k )->GetLogicalVolume() == calorimeterVolume )
0465             {
0466                 if ( G4StrUtil::contains(( *k )->GetName(), "Right" ) )
0467                     rightCalorimeter = *k;
0468                 break;
0469             }
0470         } while ( false );
0471     }
0472 }
0473 
0474 
0475 void  CexmcSetup::AssertAndAsignDetectorRole( CexmcDetectorRole &  detectorRole,
0476                                               CexmcDetectorRole  value )
0477 {
0478     if ( detectorRole != CexmcNumberOfDetectorRoles && detectorRole != value )
0479         throw CexmcException( CexmcMultipleDetectorRoles );
0480 
0481     detectorRole = value;
0482 }
0483 
0484 
0485 void  CexmcSetup::RotateMatrix( const G4ThreeVector &  rot,
0486                                 G4RotationMatrix &  rm )
0487 {
0488     rm.rotateX( rot.x() );
0489     rm.rotateY( rot.y() );
0490     rm.rotateZ( rot.z() );
0491 }
0492