Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/ChargeExchangeMC/src/CexmcEnergyDepositDigitizer.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 /*
 27  * ============================================================================
 28  *
 29  *       Filename:  CexmcEnergyDepositDigitizer.cc
 30  *
 31  *    Description:  digitizes of energy deposit in a single event
 32  *
 33  *        Version:  1.0
 34  *        Created:  23.11.2009 14:39:41
 35  *       Revision:  none
 36  *       Compiler:  gcc
 37  *
 38  *         Author:  Alexey Radkov (), 
 39  *        Company:  PNPI
 40  *
 41  * ============================================================================
 42  */
 43 
 44 #include <iostream>
 45 #include <iomanip>
 46 #include <G4DigiManager.hh>
 47 #include <G4String.hh>
 48 #include <Randomize.hh>
 49 #include "CexmcEnergyDepositDigitizer.hh"
 50 #include "CexmcEnergyDepositDigitizerMessenger.hh"
 51 #include "CexmcSimpleEnergyDeposit.hh"
 52 #include "CexmcEnergyDepositInLeftRightSet.hh"
 53 #include "CexmcEnergyDepositInCalorimeter.hh"
 54 #include "CexmcSetup.hh"
 55 #include "CexmcRunManager.hh"
 56 #include "CexmcSensitiveDetectorsAttributes.hh"
 57 
 58 
 59 CexmcEnergyDepositDigitizer::CexmcEnergyDepositDigitizer(
 60                                                     const G4String &  name ) :
 61     G4VDigitizerModule( name ), monitorED( 0 ),
 62     vetoCounterEDLeft( 0 ), vetoCounterEDRight( 0 ),
 63     calorimeterEDLeft( 0 ), calorimeterEDRight( 0 ),
 64     calorimeterEDLeftMaxX( 0 ), calorimeterEDLeftMaxY( 0 ),
 65     calorimeterEDRightMaxX( 0 ), calorimeterEDRightMaxY( 0 ),
 66     monitorHasTriggered( false ), hasTriggered( false ),
 67     monitorEDThreshold( 0 ),
 68     vetoCounterEDLeftThreshold( 0 ), vetoCounterEDRightThreshold( 0 ),
 69     calorimeterEDLeftThreshold( 0 ), calorimeterEDRightThreshold( 0 ),
 70     calorimeterTriggerAlgorithm( CexmcAllCrystalsMakeEDTriggerThreshold ),
 71     outerCrystalsVetoAlgorithm( CexmcNoOuterCrystalsVeto ),
 72     outerCrystalsVetoFraction( 0 ), monitorEDThresholdRef( 0 ),
 73     vetoCounterEDLeftThresholdRef( 0 ), vetoCounterEDRightThresholdRef( 0 ),
 74     calorimeterEDLeftThresholdRef( 0 ), calorimeterEDRightThresholdRef( 0 ),
 75     calorimeterTriggerAlgorithmRef( CexmcAllCrystalsMakeEDTriggerThreshold ),
 76     outerCrystalsVetoAlgorithmRef( CexmcNoOuterCrystalsVeto ),
 77     outerCrystalsVetoFractionRef( 0 ), nCrystalsInColumn( 1 ),
 78     nCrystalsInRow( 1 ), applyFiniteCrystalResolution( false ),
 79     messenger( NULL )
 80 {
 81     G4RunManager *      runManager( G4RunManager::GetRunManager() );
 82     const CexmcSetup *  setup( static_cast< const CexmcSetup * >(
 83                                 runManager->GetUserDetectorConstruction() ) );
 84     const CexmcSetup::CalorimeterGeometryData &  calorimeterGeometry(
 85                                             setup->GetCalorimeterGeometry() );
 86 
 87     nCrystalsInColumn = calorimeterGeometry.nCrystalsInColumn;
 88     nCrystalsInRow = calorimeterGeometry.nCrystalsInRow;
 89 
 90     if ( nCrystalsInColumn > 0 )
 91     {
 92         calorimeterEDLeftCollection.resize( nCrystalsInColumn );
 93         calorimeterEDRightCollection.resize( nCrystalsInColumn );
 94     }
 95 
 96     if ( nCrystalsInRow > 0 )
 97     {
 98         for ( CexmcEnergyDepositCalorimeterCollection::iterator
 99                 k( calorimeterEDLeftCollection.begin() );
100                     k != calorimeterEDLeftCollection.end(); ++k )
101         {
102             k->resize( nCrystalsInRow );
103         }
104         for ( CexmcEnergyDepositCalorimeterCollection::iterator
105                 k( calorimeterEDRightCollection.begin() );
106                     k != calorimeterEDRightCollection.end(); ++k )
107         {
108             k->resize( nCrystalsInRow );
109         }
110     }
111 
112     messenger = new CexmcEnergyDepositDigitizerMessenger( this );
113 }
114 
115 
116 CexmcEnergyDepositDigitizer::~CexmcEnergyDepositDigitizer()
117 {
118     delete messenger;
119 }
120 
121 
122 void  CexmcEnergyDepositDigitizer::InitializeData( void )
123 {
124     monitorED = 0;
125     vetoCounterEDLeft = 0;
126     vetoCounterEDRight = 0;
127     calorimeterEDLeft = 0;
128     calorimeterEDRight = 0;
129     calorimeterEDLeftMaxX = 0;
130     calorimeterEDLeftMaxY = 0;
131     calorimeterEDRightMaxX = 0;
132     calorimeterEDRightMaxY = 0;
133     monitorHasTriggered = false;
134     hasTriggered = false;
135 
136     for ( CexmcEnergyDepositCalorimeterCollection::iterator
137               k( calorimeterEDLeftCollection.begin() );
138                   k != calorimeterEDLeftCollection.end(); ++k )
139     {
140         for ( CexmcEnergyDepositCrystalRowCollection::iterator
141                 l( k->begin() ); l != k->end(); ++l )
142         {
143             *l = 0;
144         }
145     }
146     for ( CexmcEnergyDepositCalorimeterCollection::iterator
147               k( calorimeterEDRightCollection.begin() );
148                   k != calorimeterEDRightCollection.end(); ++k )
149     {
150         for ( CexmcEnergyDepositCrystalRowCollection::iterator
151                 l( k->begin() ); l != k->end(); ++l )
152         {
153             *l = 0;
154         }
155     }
156 }
157 
158 
159 void  CexmcEnergyDepositDigitizer::Digitize( void )
160 {
161     InitializeData();
162 
163     G4DigiManager *  digiManager( G4DigiManager::GetDMpointer() );
164     G4int    hcId( digiManager->GetHitsCollectionID(
165                     CexmcDetectorRoleName[ CexmcMonitorDetectorRole ] +
166                     "/" + CexmcDetectorTypeName[ CexmcEDDetector ] ) );
167     const CexmcEnergyDepositCollection *
168          hitsCollection( static_cast< const CexmcEnergyDepositCollection * >(
169                                     digiManager->GetHitsCollection( hcId ) ) );
170 
171     if ( hitsCollection )
172     {
173         /* it always must have index 0 */
174         if ( ( *hitsCollection )[ 0 ] )
175             monitorED = *( *hitsCollection )[ 0 ];
176     }
177 
178     hcId = digiManager->GetHitsCollectionID(
179                     CexmcDetectorRoleName[ CexmcVetoCounterDetectorRole ] +
180                     "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
181     hitsCollection = static_cast< const CexmcEnergyDepositCollection * >(
182                                     digiManager->GetHitsCollection( hcId ) );
183     if ( hitsCollection )
184     {
185         for ( CexmcEnergyDepositCollectionData::iterator
186                   k( hitsCollection->GetMap()->begin() );
187                       k != hitsCollection->GetMap()->end(); ++k )
188         {
189             G4int      index( k->first );
190             CexmcSide  side( CexmcEnergyDepositInLeftRightSet::GetSide(
191                                                                    index ) );
192             switch ( side )
193             {
194             case CexmcLeft :
195                 vetoCounterEDLeft = *k->second;
196                 break;
197             case CexmcRight :
198                 vetoCounterEDRight = *k->second;
199                 break;
200             default :
201                 break;
202             }
203         }
204     }
205 
206     G4double  maxEDCrystalLeft( 0 );
207     G4double  maxEDCrystalRight( 0 );
208     G4double  outerCrystalsEDLeft( 0 );
209     G4double  outerCrystalsEDRight( 0 );
210     G4double  innerCrystalsEDLeft( 0 );
211     G4double  innerCrystalsEDRight( 0 );
212 
213     CexmcRunManager *  runManager( static_cast< CexmcRunManager * >(
214                                             G4RunManager::GetRunManager() ) );
215 
216     hcId = digiManager->GetHitsCollectionID(
217                     CexmcDetectorRoleName[ CexmcCalorimeterDetectorRole ] +
218                     "/" + CexmcDetectorTypeName[ CexmcEDDetector ] );
219     hitsCollection = static_cast< const CexmcEnergyDepositCollection * >(
220                                     digiManager->GetHitsCollection( hcId ) );
221     if ( hitsCollection )
222     {
223         for ( CexmcEnergyDepositCollectionData::iterator
224                   k( hitsCollection->GetMap()->begin() );
225                       k != hitsCollection->GetMap()->end(); ++k )
226         {
227             G4int      index( k->first );
228             CexmcSide  side( CexmcEnergyDepositInLeftRightSet::GetSide(
229                                                                    index ) );
230             G4int      row( CexmcEnergyDepositInCalorimeter::GetRow( index ) );
231             G4int      column( CexmcEnergyDepositInCalorimeter::GetColumn(
232                                                                    index ) );
233             G4double   value( *k->second );
234             if ( applyFiniteCrystalResolution && value > 0. &&
235                                                  ! runManager->ProjectIsRead() )
236             {
237                 for ( CexmcEnergyRangeWithDoubleValueList::const_iterator
238                           l( crystalResolutionData.begin() );
239                           l != crystalResolutionData.end(); ++l )
240                 {
241                     if ( value < l->bottom || value >= l->top )
242                         continue;
243                     value = G4RandGauss::shoot( value,
244                                         value * l->value * CexmcFwhmToStddev );
245                     if ( value < 0. )
246                         value = 0.;
247                     break;
248                 }
249             }
250             switch ( side )
251             {
252             case CexmcLeft :
253                 if ( value > maxEDCrystalLeft )
254                 {
255                     calorimeterEDLeftMaxX = column;
256                     calorimeterEDLeftMaxY = row;
257                     maxEDCrystalLeft = value;
258                 }
259                 if ( IsOuterCrystal( column, row ) )
260                 {
261                     outerCrystalsEDLeft += value;
262                 }
263                 else
264                 {
265                     innerCrystalsEDLeft += value;
266                 }
267                 calorimeterEDLeft += value;
268                 calorimeterEDLeftCollection[ row ][ column ] = value;
269                 break;
270             case CexmcRight :
271                 if ( value > maxEDCrystalRight )
272                 {
273                     calorimeterEDRightMaxX = column;
274                     calorimeterEDRightMaxY = row;
275                     maxEDCrystalRight = value;
276                 }
277                 if ( IsOuterCrystal( column, row ) )
278                 {
279                     outerCrystalsEDRight += value;
280                 }
281                 else
282                 {
283                     innerCrystalsEDRight += value;
284                 }
285                 calorimeterEDRight += value;
286                 calorimeterEDRightCollection[ row ][ column ] = value;
287                 break;
288             default :
289                 break;
290             }
291         }
292     }
293 
294     G4double  calorimeterEDLeftEffective( calorimeterEDLeft );
295     G4double  calorimeterEDRightEffective( calorimeterEDRight );
296 
297     if ( calorimeterTriggerAlgorithm ==
298          CexmcInnerCrystalsMakeEDTriggerThreshold )
299     {
300         calorimeterEDLeftEffective = innerCrystalsEDLeft;
301         calorimeterEDRightEffective = innerCrystalsEDRight;
302     }
303 
304     monitorHasTriggered = monitorED >= monitorEDThreshold;
305 
306     hasTriggered = monitorHasTriggered &&
307                    vetoCounterEDLeft < vetoCounterEDLeftThreshold &&
308                    vetoCounterEDRight < vetoCounterEDRightThreshold &&
309                    calorimeterEDLeftEffective >= calorimeterEDLeftThreshold &&
310                    calorimeterEDRightEffective >= calorimeterEDRightThreshold;
311 
312     /* event won't trigger if outer crystals veto triggered */
313     if ( hasTriggered )
314     {
315         switch ( outerCrystalsVetoAlgorithm )
316         {
317         case CexmcNoOuterCrystalsVeto :
318             break;
319         case CexmcMaximumEDInASingleOuterCrystalVeto :
320             hasTriggered =
321                     ! IsOuterCrystal( calorimeterEDLeftMaxX,
322                                       calorimeterEDLeftMaxY ) &&
323                     ! IsOuterCrystal( calorimeterEDRightMaxX,
324                                       calorimeterEDRightMaxY );
325             break;
326         case CexmcFractionOfEDInOuterCrystalsVeto :
327             hasTriggered =
328                     ( ( outerCrystalsEDLeft / calorimeterEDLeft ) <
329                                               outerCrystalsVetoFraction ) &&
330                     ( ( outerCrystalsEDRight / calorimeterEDRight ) <
331                                               outerCrystalsVetoFraction );
332             break;
333         default :
334             break;
335         }
336     }
337 }
338 
339 
340 std::ostream &  operator<<( std::ostream &  out,
341                 const CexmcEnergyDepositCalorimeterCollection &  edCollection )
342 {
343     std::streamsize  prec( out.precision() );
344 
345     out.precision( 4 );
346 
347     out << std::endl;
348     for ( CexmcEnergyDepositCalorimeterCollection::const_reverse_iterator
349             k( edCollection.rbegin() ); k != edCollection.rend(); ++k )
350     {
351         for ( CexmcEnergyDepositCrystalRowCollection::const_reverse_iterator
352                 l( k->rbegin() ); l != k->rend(); ++l )
353             out << std::setw( 10 ) << *l;
354         out << std::endl;
355     }
356 
357     out.precision( prec );
358 
359     return out;
360 }
361 
362