| Geant4 Cross Reference |
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: CexmcReconstructor.cc
30 *
31 * Description: reconstructor base class
32 *
33 * Version: 1.0
34 * Created: 02.12.2009 16:21:15
35 * Revision: none
36 * Compiler: gcc
37 *
38 * Author: Alexey Radkov (),
39 * Company: PNPI
40 *
41 * ============================================================================
42 */
43
44 #include "CexmcReconstructor.hh"
45 #include "CexmcReconstructorMessenger.hh"
46 #include "CexmcEnergyDepositStore.hh"
47 #include "CexmcRunManager.hh"
48
49
50 CexmcReconstructor::CexmcReconstructor() : hasBasicTrigger( false ),
51 epDefinitionAlgorithm( CexmcEntryPointBySqrtEDWeights ),
52 epDepthDefinitionAlgorithm( CexmcEntryPointDepthPlain ),
53 csAlgorithm( CexmcSelectAllCrystals ), useInnerRefCrystal( false ),
54 epDepth( 0 ), theAngle( 0 ), calorimeterEDLeftAdjacent( 0 ),
55 calorimeterEDRightAdjacent( 0 ), collectEDInAdjacentCrystals( false ),
56 targetEPInitialized( false ), messenger( NULL )
57 {
58 G4RunManager * runManager( G4RunManager::GetRunManager() );
59 const CexmcSetup * setup( static_cast< const CexmcSetup * >(
60 runManager->GetUserDetectorConstruction() ) );
61 calorimeterGeometry = setup->GetCalorimeterGeometry();
62 targetTransform = setup->GetTargetTransform();
63 calorimeterLeftTransform = setup->GetCalorimeterLeftTransform();
64 calorimeterRightTransform = setup->GetCalorimeterRightTransform();
65
66 messenger = new CexmcReconstructorMessenger( this );
67 }
68
69
70 CexmcReconstructor::~CexmcReconstructor()
71 {
72 delete messenger;
73 }
74
75
76 void CexmcReconstructor::Reconstruct(
77 const CexmcEnergyDepositStore * edStore )
78 {
79 ReconstructEntryPoints( edStore );
80 if ( hasBasicTrigger )
81 ReconstructTargetPoint();
82 if ( hasBasicTrigger )
83 ReconstructAngle();
84 }
85
86
87 G4bool CexmcReconstructor::HasFullTrigger( void ) const
88 {
89 return hasBasicTrigger;
90 }
91
92
93 void CexmcReconstructor::ReconstructEntryPoints(
94 const CexmcEnergyDepositStore * edStore )
95 {
96 G4int columnLeft( edStore->calorimeterEDLeftMaxX );
97 G4int rowLeft( edStore->calorimeterEDLeftMaxY );
98 G4int columnRight( edStore->calorimeterEDRightMaxX );
99 G4int rowRight( edStore->calorimeterEDRightMaxY );
100 G4double crystalLength( calorimeterGeometry.crystalLength );
101
102 if ( useInnerRefCrystal )
103 {
104 TransformToAdjacentInnerCrystal( columnLeft, rowLeft );
105 TransformToAdjacentInnerCrystal( columnRight, rowRight );
106 }
107
108 calorimeterEPLeftPosition.setX( 0 );
109 calorimeterEPLeftPosition.setY( 0 );
110 calorimeterEPLeftPosition.setZ( -crystalLength / 2 + epDepth );
111 calorimeterEPLeftDirection.setX( 0 );
112 calorimeterEPLeftDirection.setY( 0 );
113 calorimeterEPLeftDirection.setZ( 0 );
114 calorimeterEPRightPosition.setX( 0 );
115 calorimeterEPRightPosition.setY( 0 );
116 calorimeterEPRightPosition.setZ( -crystalLength / 2 + epDepth );
117 calorimeterEPRightDirection.setX( 0 );
118 calorimeterEPRightDirection.setY( 0 );
119 calorimeterEPRightDirection.setZ( 0 );
120
121 G4bool edInAdjacentCrystalsCollected( false );
122
123 switch ( epDefinitionAlgorithm )
124 {
125 case CexmcEntryPointInTheCenter :
126 break;
127 case CexmcEntryPointInTheCenterOfCrystalWithMaxED :
128 {
129 G4int nCrystalsInColumn(
130 calorimeterGeometry.nCrystalsInColumn );
131 G4int nCrystalsInRow( calorimeterGeometry.nCrystalsInRow );
132 G4double crystalWidth( calorimeterGeometry.crystalWidth );
133 G4double crystalHeight( calorimeterGeometry.crystalHeight );
134
135 calorimeterEPLeftPosition.setX(
136 ( G4double( columnLeft ) - G4double( nCrystalsInRow ) / 2 ) *
137 crystalWidth + crystalWidth / 2 );
138 calorimeterEPLeftPosition.setY(
139 ( G4double( rowLeft ) - G4double( nCrystalsInColumn ) / 2 ) *
140 crystalHeight + crystalHeight / 2 );
141 calorimeterEPRightPosition.setX(
142 ( G4double( columnRight ) - G4double( nCrystalsInRow ) / 2 ) *
143 crystalWidth + crystalWidth / 2 );
144 calorimeterEPRightPosition.setY(
145 ( G4double( rowRight ) - G4double( nCrystalsInColumn ) / 2 ) *
146 crystalHeight + crystalHeight / 2 );
147 }
148 break;
149 case CexmcEntryPointByLinearEDWeights :
150 case CexmcEntryPointBySqrtEDWeights :
151 {
152 G4double x( 0 );
153 G4double y( 0 );
154
155 CalculateWeightedEPPosition( edStore->calorimeterEDLeftCollection,
156 rowLeft, columnLeft, x, y,
157 calorimeterEDLeftAdjacent );
158 calorimeterEPLeftPosition.setX( x );
159 calorimeterEPLeftPosition.setY( y );
160 CalculateWeightedEPPosition( edStore->calorimeterEDRightCollection,
161 rowRight, columnRight, x, y,
162 calorimeterEDRightAdjacent );
163 calorimeterEPRightPosition.setX( x );
164 calorimeterEPRightPosition.setY( y );
165
166 if ( csAlgorithm == CexmcSelectAdjacentCrystals )
167 edInAdjacentCrystalsCollected = true;
168 }
169 break;
170 default :
171 break;
172 }
173
174 switch ( epDepthDefinitionAlgorithm )
175 {
176 case CexmcEntryPointDepthPlain :
177 break;
178 case CexmcEntryPointDepthSphere :
179 {
180 G4double calorimeterEPLeftRadiusOfTheSphere(
181 calorimeterLeftTransform.NetTranslation().mag() +
182 calorimeterEPLeftPosition.z() );
183 G4double calorimeterEPLeftRadiusOfTheSphere2(
184 calorimeterEPLeftRadiusOfTheSphere *
185 calorimeterEPLeftRadiusOfTheSphere );
186 G4double calorimeterEPLeftPositionX2(
187 calorimeterEPLeftPosition.x() *
188 calorimeterEPLeftPosition.x() );
189 G4double calorimeterEPLeftPositionY2(
190 calorimeterEPLeftPosition.y() *
191 calorimeterEPLeftPosition.y() );
192 G4double calorimeterEPLeftPositionZOffset(
193 calorimeterEPLeftRadiusOfTheSphere - std::sqrt(
194 calorimeterEPLeftRadiusOfTheSphere2 -
195 calorimeterEPLeftPositionX2 -
196 calorimeterEPLeftPositionY2 ) );
197 G4double calorimeterEPRightRadiusOfTheSphere(
198 calorimeterRightTransform.NetTranslation().mag() +
199 calorimeterEPRightPosition.z() );
200 G4double calorimeterEPRightRadiusOfTheSphere2(
201 calorimeterEPRightRadiusOfTheSphere *
202 calorimeterEPRightRadiusOfTheSphere );
203 G4double calorimeterEPRightPositionX2(
204 calorimeterEPRightPosition.x() *
205 calorimeterEPRightPosition.x() );
206 G4double calorimeterEPRightPositionY2(
207 calorimeterEPRightPosition.y() *
208 calorimeterEPRightPosition.y() );
209 G4double calorimeterEPRightPositionZOffset(
210 calorimeterEPRightRadiusOfTheSphere - std::sqrt(
211 calorimeterEPRightRadiusOfTheSphere2 -
212 calorimeterEPRightPositionX2 -
213 calorimeterEPRightPositionY2 ) );
214 calorimeterEPLeftPosition.setZ( calorimeterEPLeftPosition.z() -
215 calorimeterEPLeftPositionZOffset );
216 calorimeterEPRightPosition.setZ( calorimeterEPRightPosition.z() -
217 calorimeterEPRightPositionZOffset );
218 }
219 break;
220 default :
221 break;
222 }
223
224 calorimeterEPLeftWorldPosition = calorimeterLeftTransform.TransformPoint(
225 calorimeterEPLeftPosition );
226 calorimeterEPLeftWorldDirection = calorimeterLeftTransform.TransformAxis(
227 calorimeterEPLeftDirection );
228 calorimeterEPRightWorldPosition = calorimeterRightTransform.TransformPoint(
229 calorimeterEPRightPosition );
230 calorimeterEPRightWorldDirection = calorimeterRightTransform.TransformAxis(
231 calorimeterEPRightDirection );
232
233 if ( collectEDInAdjacentCrystals && ! edInAdjacentCrystalsCollected )
234 {
235 CollectEDInAdjacentCrystals( edStore->calorimeterEDLeftCollection,
236 rowLeft, columnLeft,
237 calorimeterEDLeftAdjacent );
238 CollectEDInAdjacentCrystals( edStore->calorimeterEDRightCollection,
239 rowRight, columnRight,
240 calorimeterEDRightAdjacent );
241 }
242
243 hasBasicTrigger = true;
244 }
245
246
247 void CexmcReconstructor::ReconstructTargetPoint( void )
248 {
249 if ( ! targetEPInitialized )
250 {
251 targetEPWorldPosition = targetTransform.TransformPoint(
252 G4ThreeVector( 0, 0, 0 ) );
253 targetEPWorldDirection.setX( 0 );
254 targetEPWorldDirection.setY( 0 );
255 targetEPWorldDirection.setZ( 1 );
256
257 targetEPPosition = targetTransform.Inverse().TransformPoint(
258 targetEPWorldPosition );
259 targetEPDirection = targetTransform.Inverse().TransformAxis(
260 targetEPWorldDirection );
261 targetEPInitialized = true;
262 }
263
264 hasBasicTrigger = true;
265 }
266
267
268 void CexmcReconstructor::ReconstructAngle( void )
269 {
270 G4ThreeVector epLeft( calorimeterEPLeftWorldPosition -
271 targetEPWorldPosition );
272 G4ThreeVector epRight( calorimeterEPRightWorldPosition -
273 targetEPWorldPosition );
274 theAngle = epLeft.angle( epRight );
275
276 hasBasicTrigger = true;
277 }
278
279
280 void CexmcReconstructor::CollectEDInAdjacentCrystals(
281 const CexmcEnergyDepositCalorimeterCollection & edHits,
282 G4int row, G4int column, G4double & ed )
283 {
284 G4int i( 0 );
285
286 for ( CexmcEnergyDepositCalorimeterCollection::const_iterator
287 k( edHits.begin() ); k != edHits.end(); ++k )
288 {
289 if ( i - row > 1 || i - row < -1 )
290 {
291 ++i;
292 continue;
293 }
294
295 G4int j( 0 );
296 for ( CexmcEnergyDepositCrystalRowCollection::const_iterator
297 l( k->begin() ); l != k->end(); ++l )
298 {
299 if ( j - column > 1 || j - column < -1 )
300 {
301 ++j;
302 continue;
303 }
304 ed += *l;
305 ++j;
306 }
307 ++i;
308 }
309 }
310
311
312 void CexmcReconstructor::CalculateWeightedEPPosition(
313 const CexmcEnergyDepositCalorimeterCollection & edHits,
314 G4int row, G4int column, G4double & x, G4double & y,
315 G4double & ed )
316 {
317 G4int nCrystalsInColumn( calorimeterGeometry.nCrystalsInColumn );
318 G4int nCrystalsInRow( calorimeterGeometry.nCrystalsInRow );
319 G4double crystalWidth( calorimeterGeometry.crystalWidth );
320 G4double crystalHeight( calorimeterGeometry.crystalHeight );
321
322 G4int i( 0 );
323 G4double xWeightsSum( 0 );
324 G4double yWeightsSum( 0 );
325 G4double energyWeightsSum( 0 );
326
327 if ( csAlgorithm == CexmcSelectAdjacentCrystals )
328 ed = 0.;
329
330 for ( CexmcEnergyDepositCalorimeterCollection::const_iterator
331 k( edHits.begin() ); k != edHits.end(); ++k )
332 {
333 if ( csAlgorithm == CexmcSelectAdjacentCrystals &&
334 ( i - row > 1 || i - row < -1 ) )
335 {
336 ++i;
337 continue;
338 }
339
340 G4int j( 0 );
341 for ( CexmcEnergyDepositCrystalRowCollection::const_iterator
342 l( k->begin() ); l != k->end(); ++l )
343 {
344 if ( csAlgorithm == CexmcSelectAdjacentCrystals &&
345 ( j - column > 1 || j - column < -1 ) )
346 {
347 ++j;
348 continue;
349 }
350
351 if ( csAlgorithm == CexmcSelectAdjacentCrystals )
352 ed += *l;
353
354 G4double xInCalorimeterOffset(
355 ( G4double( j ) - G4double( nCrystalsInRow ) / 2 ) *
356 crystalWidth + crystalWidth / 2 );
357 G4double energyWeight(
358 epDefinitionAlgorithm ==
359 CexmcEntryPointBySqrtEDWeights ?
360 std::sqrt( *l ) : *l );
361 xWeightsSum += energyWeight * xInCalorimeterOffset;
362 G4double yInCalorimeterOffset(
363 ( G4double( i ) - G4double( nCrystalsInColumn ) / 2 ) *
364 crystalHeight + crystalHeight / 2 );
365 yWeightsSum += energyWeight * yInCalorimeterOffset;
366 energyWeightsSum += energyWeight;
367 ++j;
368 }
369 ++i;
370 }
371
372 x = xWeightsSum / energyWeightsSum;
373 y = yWeightsSum / energyWeightsSum;
374 }
375
376