Geant4  10.02.p02
CexmcChargeExchangeReconstructor.cc
Go to the documentation of this file.
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 
213  productionModelData.incidentParticleSCM.boost( -boostVec );
214  productionModelData.nucleusParticleSCM.boost( -boostVec );
215  productionModelData.outputParticleSCM.boost( -boostVec );
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 
const G4ParticleDefinition * nucleusParticle
G4ThreeVector calorimeterEPLeftWorldPosition
const G4VUserPrimaryGeneratorAction * GetUserPrimaryGeneratorAction() const
CLHEP::Hep3Vector G4ThreeVector
const G4ParticleDefinition * incidentParticle
G4double GetOrigMomentumAmp(void) const
G4ParticleDefinition * GetIncidentParticle(void) const
CexmcChargeExchangeReconstructorMessenger * messenger
const G4ParticleDefinition * outputParticle
G4ThreeVector targetEPWorldPosition
void ReconstructEntryPoints(const CexmcEnergyDepositStore *edStore)
CexmcChargeExchangeReconstructor(const CexmcProductionModel *productionModel)
void Reconstruct(const CexmcEnergyDepositStore *edStore)
bool G4bool
Definition: G4Types.hh:79
G4ParticleDefinition * GetOutputParticle(void) const
void ReconstructTargetPoint(void)
G4ThreeVector calorimeterEPRightWorldPosition
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
G4double GetPDGMass() const
G4ParticleDefinition * GetNucleusParticle(void) const
G4ParticleDefinition * GetParticleDefinition() const
G4ParticleDefinition * GetNucleusOutputParticle(void) const
double G4double
Definition: G4Types.hh:76
const G4ParticleDefinition * nucleusOutputParticle
const G4ThreeVector & GetOrigDirection(void) const
CLHEP::HepLorentzVector G4LorentzVector