ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CexmcReconstructor.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CexmcReconstructor.cc
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