ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CexmcChargeExchangeReconstructor.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CexmcChargeExchangeReconstructor.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: CexmcChargeExchangeReconstructor.cc
30  *
31  * Description: charge exchange reconstructor
32  *
33  * Version: 1.0
34  * Created: 02.12.2009 15:17:13
35  * Revision: none
36  * Compiler: gcc
37  *
38  * Author: Alexey Radkov (),
39  * Company: PNPI
40  *
41  * ============================================================================
42  */
43 
44 #include <cmath>
45 #include <G4ThreeVector.hh>
46 #include <G4LorentzVector.hh>
51 #include "CexmcParticleGun.hh"
52 #include "CexmcProductionModel.hh"
53 #include "CexmcRunManager.hh"
54 #include "CexmcException.hh"
55 
56 
58  const CexmcProductionModel * productionModel ) :
59  outputParticleMass( 0 ), nucleusOutputParticleMass( 0 ),
60  useTableMass( false ), useMassCut( false ), massCutOPCenter( 0 ),
61  massCutNOPCenter( 0 ), massCutOPWidth( 0 ), massCutNOPWidth( 0 ),
62  massCutEllipseAngle( 0 ), useAbsorbedEnergyCut( false ),
63  absorbedEnergyCutCLCenter( 0 ), absorbedEnergyCutCRCenter( 0 ),
64  absorbedEnergyCutCLWidth( 0 ), absorbedEnergyCutCRWidth( 0 ),
65  absorbedEnergyCutEllipseAngle( 0 ), expectedMomentumAmp( -1 ),
66  edCollectionAlgorithm( CexmcCollectEDInAllCrystals ),
67  hasMassCutTriggered( false ), hasAbsorbedEnergyCutTriggered( false ),
68  beamParticleIsInitialized( false ), particleGun( NULL ), messenger( NULL )
69 {
70  if ( ! productionModel )
72 
74  productionModel->GetIncidentParticle();
75 
76  CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
78  const CexmcPrimaryGeneratorAction * primaryGeneratorAction(
79  static_cast< const CexmcPrimaryGeneratorAction * >(
80  runManager->GetUserPrimaryGeneratorAction() ) );
81  CexmcPrimaryGeneratorAction * thePrimaryGeneratorAction(
82  const_cast< CexmcPrimaryGeneratorAction * >(
83  primaryGeneratorAction ) );
84  particleGun = thePrimaryGeneratorAction->GetParticleGun();
85 
87  productionModel->GetNucleusParticle();
89  productionModel->GetOutputParticle();
91  productionModel->GetNucleusOutputParticle();
92 
94 }
95 
96 
98 {
99  delete messenger;
100 }
101 
102 
104 {
108 
110 }
111 
112 
114  const CexmcEnergyDepositStore * edStore )
115 {
117  {
121 
123  }
124 
127 
128  ReconstructEntryPoints( edStore );
129  if ( hasBasicTrigger )
131  if ( hasBasicTrigger )
133 
138 
139  G4double cosTheAngle( std::cos( theAngle ) );
140  G4double calorimeterEDLeft( edStore->calorimeterEDLeft );
141  G4double calorimeterEDRight( edStore->calorimeterEDRight );
142 
144  {
145  calorimeterEDLeft = calorimeterEDLeftAdjacent;
146  calorimeterEDRight = calorimeterEDRightAdjacent;
147  }
148 
149  //G4double cosOutputParticleLAB(
150  //( calorimeterEDLeft * cosAngleLeft +
151  //calorimeterEDRight * cosAngleRight ) /
152  //std::sqrt( calorimeterEDLeft * calorimeterEDLeft +
153  //calorimeterEDRight * calorimeterEDRight +
154  //calorimeterEDLeft * calorimeterEDRight * cosTheAngle ) );
155 
156  outputParticleMass = std::sqrt( 2 * calorimeterEDLeft *
157  calorimeterEDRight * ( 1 - cosTheAngle ) );
158 
159  G4ThreeVector opdpLeftMomentum( epLeft );
160  opdpLeftMomentum.setMag( calorimeterEDLeft );
161  G4ThreeVector opdpRightMomentum( epRight );
162  opdpRightMomentum.setMag( calorimeterEDRight );
163  G4ThreeVector opMomentum( opdpLeftMomentum + opdpRightMomentum );
164 
165  /* opMass will be used only in calculation of output particle's total
166  * energy, in other places outputParticleMass should be used instead */
167  G4double opMass( useTableMass ?
170  /* the formula below is equivalent to
171  * calorimeterEDLeft + calorimeterEDRight if opMass = outputParticleMass */
172  G4double opEnergy( std::sqrt( opMomentum.mag2() +
173  opMass * opMass ) );
175  opEnergy );
176 
177  G4ThreeVector incidentParticleMomentum( particleGun->GetOrigDirection() );
178  G4double incidentParticleMomentumAmp( expectedMomentumAmp > 0 ?
180  incidentParticleMomentum *= incidentParticleMomentumAmp;
181 
182  G4double incidentParticlePDGMass(
184  G4double incidentParticlePDGMass2( incidentParticlePDGMass *
185  incidentParticlePDGMass );
186  G4double incidentParticleEnergy(
187  std::sqrt( incidentParticleMomentumAmp * incidentParticleMomentumAmp +
188  incidentParticlePDGMass2 ) );
189 
191  incidentParticleMomentum, incidentParticleEnergy );
192  G4double nucleusParticlePDGMass(
195  G4ThreeVector( 0, 0, 0 ), nucleusParticlePDGMass );
196 
199  G4ThreeVector boostVec( lVecSum.boostVector() );
200 
203 
212 
217 
218  G4ThreeVector nopMomentum( incidentParticleMomentum - opMomentum );
219  G4double nopEnergy( incidentParticleEnergy + nucleusParticlePDGMass -
220  opEnergy );
221  nucleusOutputParticleMass = std::sqrt( nopEnergy * nopEnergy -
222  nopMomentum.mag2() );
223 
224  if ( useMassCut )
225  {
226  G4double cosMassCutEllipseAngle( std::cos( massCutEllipseAngle ) );
227  G4double sinMassCutEllipseAngle( std::sin( massCutEllipseAngle ) );
228 
229  if ( massCutOPWidth <= 0. || massCutNOPWidth <= 0. )
230  {
231  hasMassCutTriggered = false;
232  }
233  else
234  {
235  G4double massCutOPWidth2( massCutOPWidth * massCutOPWidth );
236  G4double massCutNOPWidth2( massCutNOPWidth * massCutNOPWidth );
237 
239  std::pow( ( outputParticleMass - massCutOPCenter ) *
240  cosMassCutEllipseAngle +
242  sinMassCutEllipseAngle, 2 ) / massCutOPWidth2 +
243  std::pow( - ( outputParticleMass - massCutOPCenter ) *
244  sinMassCutEllipseAngle +
246  cosMassCutEllipseAngle, 2 ) / massCutNOPWidth2 <
247  1;
248  }
249  }
250 
251  if ( useAbsorbedEnergyCut )
252  {
253  G4double cosAbsorbedEnergyCutEllipseAngle(
254  std::cos( absorbedEnergyCutEllipseAngle ) );
255  G4double sinAbsorbedEnergyCutEllipseAngle(
256  std::sin( absorbedEnergyCutEllipseAngle ) );
257 
259  {
261  }
262  else
263  {
264  G4double absorbedEnergyCutCLWidth2(
266  G4double absorbedEnergyCutCRWidth2(
268 
270  std::pow( ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) *
271  cosAbsorbedEnergyCutEllipseAngle +
272  ( calorimeterEDRight - absorbedEnergyCutCRCenter ) *
273  sinAbsorbedEnergyCutEllipseAngle, 2 ) /
274  absorbedEnergyCutCLWidth2 +
275  std::pow( - ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) *
276  sinAbsorbedEnergyCutEllipseAngle +
277  ( calorimeterEDRight - absorbedEnergyCutCRCenter ) *
278  cosAbsorbedEnergyCutEllipseAngle, 2 ) /
279  absorbedEnergyCutCRWidth2 <
280  1;
281  }
282  }
283 
284  hasBasicTrigger = true;
285 }
286 
287 
289 {
290  if ( ! hasBasicTrigger )
291  return false;
292  if ( useMassCut && ! hasMassCutTriggered )
293  return false;
295  return false;
296 
297  return true;
298 }
299 
300 
302  G4double value )
303 {
305 }
306