ChargeExchangeMC/src/CexmcEnergyDepositDigitizer.cc

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:  CexmcEnergyDepositDigitizer.cc
0030  *
0031  *    Description:  digitizes of energy deposit in a single event
0032  *
0033  *        Version:  1.0
0034  *        Created:  23.11.2009 14:39:41
0035  *       Revision:  none
0036  *       Compiler:  gcc
0037  *
0038  *         Author:  Alexey Radkov (),
0039  *        Company:  PNPI
0040  *
0041  * ============================================================================
0042  */
0043
0044 #include <iostream>
0045 #include <iomanip>
0046 #include <G4DigiManager.hh>
0047 #include <G4String.hh>
0048 #include <Randomize.hh>
0049 #include "CexmcEnergyDepositDigitizer.hh"
0050 #include "CexmcEnergyDepositDigitizerMessenger.hh"
0051 #include "CexmcSimpleEnergyDeposit.hh"
0052 #include "CexmcEnergyDepositInLeftRightSet.hh"
0053 #include "CexmcEnergyDepositInCalorimeter.hh"
0054 #include "CexmcSetup.hh"
0055 #include "CexmcRunManager.hh"
0056 #include "CexmcSensitiveDetectorsAttributes.hh"
0057
0058
0059 CexmcEnergyDepositDigitizer::CexmcEnergyDepositDigitizer(
0060                                                     const G4String &  name ) :
0061     G4VDigitizerModule( name ), monitorED( 0 ),
0062     vetoCounterEDLeft( 0 ), vetoCounterEDRight( 0 ),
0063     calorimeterEDLeft( 0 ), calorimeterEDRight( 0 ),
0064     calorimeterEDLeftMaxX( 0 ), calorimeterEDLeftMaxY( 0 ),
0065     calorimeterEDRightMaxX( 0 ), calorimeterEDRightMaxY( 0 ),
0066     monitorHasTriggered( false ), hasTriggered( false ),
0067     monitorEDThreshold( 0 ),
0068     vetoCounterEDLeftThreshold( 0 ), vetoCounterEDRightThreshold( 0 ),
0069     calorimeterEDLeftThreshold( 0 ), calorimeterEDRightThreshold( 0 ),
0070     calorimeterTriggerAlgorithm( CexmcAllCrystalsMakeEDTriggerThreshold ),
0071     outerCrystalsVetoAlgorithm( CexmcNoOuterCrystalsVeto ),
0072     outerCrystalsVetoFraction( 0 ), monitorEDThresholdRef( 0 ),
0073     vetoCounterEDLeftThresholdRef( 0 ), vetoCounterEDRightThresholdRef( 0 ),
0074     calorimeterEDLeftThresholdRef( 0 ), calorimeterEDRightThresholdRef( 0 ),
0075     calorimeterTriggerAlgorithmRef( CexmcAllCrystalsMakeEDTriggerThreshold ),
0076     outerCrystalsVetoAlgorithmRef( CexmcNoOuterCrystalsVeto ),
0077     outerCrystalsVetoFractionRef( 0 ), nCrystalsInColumn( 1 ),
0078     nCrystalsInRow( 1 ), applyFiniteCrystalResolution( false ),
0079     messenger( NULL )
0080 {
0081     G4RunManager *      runManager( G4RunManager::GetRunManager() );
0082     const CexmcSetup *  setup( static_cast< const CexmcSetup * >(
0083                                 runManager->GetUserDetectorConstruction() ) );
0084     const CexmcSetup::CalorimeterGeometryData &  calorimeterGeometry(
0085                                             setup->GetCalorimeterGeometry() );
0086
0087     nCrystalsInColumn = calorimeterGeometry.nCrystalsInColumn;
0088     nCrystalsInRow = calorimeterGeometry.nCrystalsInRow;
0089
0090     if ( nCrystalsInColumn > 0 )
0091     {
0092         calorimeterEDLeftCollection.resize( nCrystalsInColumn );
0093         calorimeterEDRightCollection.resize( nCrystalsInColumn );
0094     }
0095
0096     if ( nCrystalsInRow > 0 )
0097     {
0098         for ( CexmcEnergyDepositCalorimeterCollection::iterator
0099                 k( calorimeterEDLeftCollection.begin() );
0100                     k != calorimeterEDLeftCollection.end(); ++k )
0101         {
0102             k->resize( nCrystalsInRow );
0103         }
0104         for ( CexmcEnergyDepositCalorimeterCollection::iterator
0105                 k( calorimeterEDRightCollection.begin() );
0106                     k != calorimeterEDRightCollection.end(); ++k )
0107         {
0108             k->resize( nCrystalsInRow );
0109         }
0110     }
0111
0112     messenger = new CexmcEnergyDepositDigitizerMessenger( this );
0113 }
0114
0115
0116 CexmcEnergyDepositDigitizer::~CexmcEnergyDepositDigitizer()
0117 {
0118     delete messenger;
0119 }
0120
0121
0122 void  CexmcEnergyDepositDigitizer::InitializeData( void )
0123 {
0124     monitorED = 0;
0125     vetoCounterEDLeft = 0;
0126     vetoCounterEDRight = 0;
0127     calorimeterEDLeft = 0;
0128     calorimeterEDRight = 0;
0129     calorimeterEDLeftMaxX = 0;
0130     calorimeterEDLeftMaxY = 0;
0131     calorimeterEDRightMaxX = 0;
0132     calorimeterEDRightMaxY = 0;
0133     monitorHasTriggered = false;
0134     hasTriggered = false;
0135
0136     for ( CexmcEnergyDepositCalorimeterCollection::iterator
0137               k( calorimeterEDLeftCollection.begin() );
0138                   k != calorimeterEDLeftCollection.end(); ++k )
0139     {
0140         for ( CexmcEnergyDepositCrystalRowCollection::iterator
0141                 l( k->begin() ); l != k->end(); ++l )
0142         {
0143             *l = 0;
0144         }
0145     }
0146     for ( CexmcEnergyDepositCalorimeterCollection::iterator
0147               k( calorimeterEDRightCollection.begin() );
0148                   k != calorimeterEDRightCollection.end(); ++k )
0149     {
0150         for ( CexmcEnergyDepositCrystalRowCollection::iterator
0151                 l( k->begin() ); l != k->end(); ++l )
0152         {
0153             *l = 0;
0154         }
0155     }
0156 }
0157
0158
0159 void  CexmcEnergyDepositDigitizer::Digitize( void )
0160 {
0161     InitializeData();
0162
0163     G4DigiManager *  digiManager( G4DigiManager::GetDMpointer() );
0164     G4int    hcId( digiManager->GetHitsCollectionID(
0165                     CexmcDetectorRoleName[ CexmcMonitorDetectorRole ] +
0166                     "/" + CexmcDetectorTypeName[ CexmcEDDetector ] ) );
0167     const CexmcEnergyDepositCollection *
0168          hitsCollection( static_cast< const CexmcEnergyDepositCollection * >(
0169                                     digiManager->GetHitsCollection( hcId ) ) );
0170
0171     if ( hitsCollection )
0172     {
0173         /* it always must have index 0 */
0174         if ( ( *hitsCollection )[ 0 ] )
0175             monitorED = *( *hitsCollection )[ 0 ];
0176     }
0177
0178     hcId = digiManager->GetHitsCollectionID(
0179                     CexmcDetectorRoleName[ CexmcVetoCounterDetectorRole ] +
0180                     "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
0181     hitsCollection = static_cast< const CexmcEnergyDepositCollection * >(
0182                                     digiManager->GetHitsCollection( hcId ) );
0183     if ( hitsCollection )
0184     {
0185         for ( CexmcEnergyDepositCollectionData::iterator
0186                   k( hitsCollection->GetMap()->begin() );
0187                       k != hitsCollection->GetMap()->end(); ++k )
0188         {
0189             G4int      index( k->first );
0190             CexmcSide  side( CexmcEnergyDepositInLeftRightSet::GetSide(
0191                                                                    index ) );
0192             switch ( side )
0193             {
0194             case CexmcLeft :
0195                 vetoCounterEDLeft = *k->second;
0196                 break;
0197             case CexmcRight :
0198                 vetoCounterEDRight = *k->second;
0199                 break;
0200             default :
0201                 break;
0202             }
0203         }
0204     }
0205
0206     G4double  maxEDCrystalLeft( 0 );
0207     G4double  maxEDCrystalRight( 0 );
0208     G4double  outerCrystalsEDLeft( 0 );
0209     G4double  outerCrystalsEDRight( 0 );
0210     G4double  innerCrystalsEDLeft( 0 );
0211     G4double  innerCrystalsEDRight( 0 );
0212
0213     CexmcRunManager *  runManager( static_cast< CexmcRunManager * >(
0214                                             G4RunManager::GetRunManager() ) );
0215
0216     hcId = digiManager->GetHitsCollectionID(
0217                     CexmcDetectorRoleName[ CexmcCalorimeterDetectorRole ] +
0218                     "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
0219     hitsCollection = static_cast< const CexmcEnergyDepositCollection * >(
0220                                     digiManager->GetHitsCollection( hcId ) );
0221     if ( hitsCollection )
0222     {
0223         for ( CexmcEnergyDepositCollectionData::iterator
0224                   k( hitsCollection->GetMap()->begin() );
0225                       k != hitsCollection->GetMap()->end(); ++k )
0226         {
0227             G4int      index( k->first );
0228             CexmcSide  side( CexmcEnergyDepositInLeftRightSet::GetSide(
0229                                                                    index ) );
0230             G4int      row( CexmcEnergyDepositInCalorimeter::GetRow( index ) );
0231             G4int      column( CexmcEnergyDepositInCalorimeter::GetColumn(
0232                                                                    index ) );
0233             G4double   value( *k->second );
0234             if ( applyFiniteCrystalResolution && value > 0. &&
0235                                                  ! runManager->ProjectIsRead() )
0236             {
0237                 for ( CexmcEnergyRangeWithDoubleValueList::const_iterator
0238                           l( crystalResolutionData.begin() );
0239                           l != crystalResolutionData.end(); ++l )
0240                 {
0241                     if ( value < l->bottom || value >= l->top )
0242                         continue;
0243                     value = G4RandGauss::shoot( value,
0244                                         value * l->value * CexmcFwhmToStddev );
0245                     if ( value < 0. )
0246                         value = 0.;
0247                     break;
0248                 }
0249             }
0250             switch ( side )
0251             {
0252             case CexmcLeft :
0253                 if ( value > maxEDCrystalLeft )
0254                 {
0255                     calorimeterEDLeftMaxX = column;
0256                     calorimeterEDLeftMaxY = row;
0257                     maxEDCrystalLeft = value;
0258                 }
0259                 if ( IsOuterCrystal( column, row ) )
0260                 {
0261                     outerCrystalsEDLeft += value;
0262                 }
0263                 else
0264                 {
0265                     innerCrystalsEDLeft += value;
0266                 }
0267                 calorimeterEDLeft += value;
0268                 calorimeterEDLeftCollection[ row ][ column ] = value;
0269                 break;
0270             case CexmcRight :
0271                 if ( value > maxEDCrystalRight )
0272                 {
0273                     calorimeterEDRightMaxX = column;
0274                     calorimeterEDRightMaxY = row;
0275                     maxEDCrystalRight = value;
0276                 }
0277                 if ( IsOuterCrystal( column, row ) )
0278                 {
0279                     outerCrystalsEDRight += value;
0280                 }
0281                 else
0282                 {
0283                     innerCrystalsEDRight += value;
0284                 }
0285                 calorimeterEDRight += value;
0286                 calorimeterEDRightCollection[ row ][ column ] = value;
0287                 break;
0288             default :
0289                 break;
0290             }
0291         }
0292     }
0293
0294     G4double  calorimeterEDLeftEffective( calorimeterEDLeft );
0295     G4double  calorimeterEDRightEffective( calorimeterEDRight );
0296
0297     if ( calorimeterTriggerAlgorithm ==
0298          CexmcInnerCrystalsMakeEDTriggerThreshold )
0299     {
0300         calorimeterEDLeftEffective = innerCrystalsEDLeft;
0301         calorimeterEDRightEffective = innerCrystalsEDRight;
0302     }
0303
0304     monitorHasTriggered = monitorED >= monitorEDThreshold;
0305
0306     hasTriggered = monitorHasTriggered &&
0307                    vetoCounterEDLeft < vetoCounterEDLeftThreshold &&
0308                    vetoCounterEDRight < vetoCounterEDRightThreshold &&
0309                    calorimeterEDLeftEffective >= calorimeterEDLeftThreshold &&
0310                    calorimeterEDRightEffective >= calorimeterEDRightThreshold;
0311
0312     /* event won't trigger if outer crystals veto triggered */
0313     if ( hasTriggered )
0314     {
0315         switch ( outerCrystalsVetoAlgorithm )
0316         {
0317         case CexmcNoOuterCrystalsVeto :
0318             break;
0319         case CexmcMaximumEDInASingleOuterCrystalVeto :
0320             hasTriggered =
0321                     ! IsOuterCrystal( calorimeterEDLeftMaxX,
0322                                       calorimeterEDLeftMaxY ) &&
0323                     ! IsOuterCrystal( calorimeterEDRightMaxX,
0324                                       calorimeterEDRightMaxY );
0325             break;
0326         case CexmcFractionOfEDInOuterCrystalsVeto :
0327             hasTriggered =
0328                     ( ( outerCrystalsEDLeft / calorimeterEDLeft ) <
0329                                               outerCrystalsVetoFraction ) &&
0330                     ( ( outerCrystalsEDRight / calorimeterEDRight ) <
0331                                               outerCrystalsVetoFraction );
0332             break;
0333         default :
0334             break;
0335         }
0336     }
0337 }
0338
0339
0340 std::ostream &  operator<<( std::ostream &  out,
0341                 const CexmcEnergyDepositCalorimeterCollection &  edCollection )
0342 {
0343     std::streamsize  prec( out.precision() );
0344
0345     out.precision( 4 );
0346
0347     out << std::endl;
0348     for ( CexmcEnergyDepositCalorimeterCollection::const_reverse_iterator
0349             k( edCollection.rbegin() ); k != edCollection.rend(); ++k )
0350     {
0351         for ( CexmcEnergyDepositCrystalRowCollection::const_reverse_iterator
0352                 l( k->rbegin() ); l != k->rend(); ++l )
0353             out << std::setw( 10 ) << *l;
0354         out << std::endl;
0355     }
0356
0357     out.precision( prec );
0358
0359     return out;
0360 }
0361