| 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: 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>
47 #include "CexmcChargeExchangeReconstructor.hh"
48 #include "CexmcChargeExchangeReconstructorMessenger.hh"
49 #include "CexmcEnergyDepositStore.hh"
50 #include "CexmcPrimaryGeneratorAction.hh"
51 #include "CexmcParticleGun.hh"
52 #include "CexmcProductionModel.hh"
53 #include "CexmcRunManager.hh"
54 #include "CexmcException.hh"
55
56
57 CexmcChargeExchangeReconstructor::CexmcChargeExchangeReconstructor(
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 )
71 throw CexmcException( CexmcWeirdException );
72
73 productionModelData.incidentParticle =
74 productionModel->GetIncidentParticle();
75
76 CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
77 G4RunManager::GetRunManager() ) );
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
86 productionModelData.nucleusParticle =
87 productionModel->GetNucleusParticle();
88 productionModelData.outputParticle =
89 productionModel->GetOutputParticle();
90 productionModelData.nucleusOutputParticle =
91 productionModel->GetNucleusOutputParticle();
92
93 messenger = new CexmcChargeExchangeReconstructorMessenger( this );
94 }
95
96
97 CexmcChargeExchangeReconstructor::~CexmcChargeExchangeReconstructor()
98 {
99 delete messenger;
100 }
101
102
103 void CexmcChargeExchangeReconstructor::SetupBeamParticle( void )
104 {
105 if ( *productionModelData.incidentParticle !=
106 *particleGun->GetParticleDefinition() )
107 throw CexmcException( CexmcBeamAndIncidentParticlesMismatch );
108
109 beamParticleIsInitialized = true;
110 }
111
112
113 void CexmcChargeExchangeReconstructor::Reconstruct(
114 const CexmcEnergyDepositStore * edStore )
115 {
116 if ( ! beamParticleIsInitialized )
117 {
118 if ( *productionModelData.incidentParticle !=
119 *particleGun->GetParticleDefinition() )
120 throw CexmcException( CexmcBeamAndIncidentParticlesMismatch );
121
122 beamParticleIsInitialized = true;
123 }
124
125 if ( edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
126 collectEDInAdjacentCrystals = true;
127
128 ReconstructEntryPoints( edStore );
129 if ( hasBasicTrigger )
130 ReconstructTargetPoint();
131 if ( hasBasicTrigger )
132 ReconstructAngle();
133
134 G4ThreeVector epLeft( calorimeterEPLeftWorldPosition -
135 targetEPWorldPosition );
136 G4ThreeVector epRight( calorimeterEPRightWorldPosition -
137 targetEPWorldPosition );
138
139 G4double cosTheAngle( std::cos( theAngle ) );
140 G4double calorimeterEDLeft( edStore->calorimeterEDLeft );
141 G4double calorimeterEDRight( edStore->calorimeterEDRight );
142
143 if ( edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
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 ?
168 productionModelData.outputParticle->GetPDGMass() :
169 outputParticleMass );
170 /* the formula below is equivalent to
171 * calorimeterEDLeft + calorimeterEDRight if opMass = outputParticleMass */
172 G4double opEnergy( std::sqrt( opMomentum.mag2() +
173 opMass * opMass ) );
174 productionModelData.outputParticleLAB = G4LorentzVector( opMomentum,
175 opEnergy );
176
177 G4ThreeVector incidentParticleMomentum( particleGun->GetOrigDirection() );
178 G4double incidentParticleMomentumAmp( expectedMomentumAmp > 0 ?
179 expectedMomentumAmp : particleGun->GetOrigMomentumAmp() );
180 incidentParticleMomentum *= incidentParticleMomentumAmp;
181
182 G4double incidentParticlePDGMass(
183 productionModelData.incidentParticle->GetPDGMass() );
184 G4double incidentParticlePDGMass2( incidentParticlePDGMass *
185 incidentParticlePDGMass );
186 G4double incidentParticleEnergy(
187 std::sqrt( incidentParticleMomentumAmp * incidentParticleMomentumAmp +
188 incidentParticlePDGMass2 ) );
189
190 productionModelData.incidentParticleLAB = G4LorentzVector(
191 incidentParticleMomentum, incidentParticleEnergy );
192 G4double nucleusParticlePDGMass(
193 productionModelData.nucleusParticle->GetPDGMass() );
194 productionModelData.nucleusParticleLAB = G4LorentzVector(
195 G4ThreeVector( 0, 0, 0 ), nucleusParticlePDGMass );
196
197 G4LorentzVector lVecSum( productionModelData.incidentParticleLAB +
198 productionModelData.nucleusParticleLAB );
199 G4ThreeVector boostVec( lVecSum.boostVector() );
200
201 productionModelData.nucleusOutputParticleLAB =
202 lVecSum - productionModelData.outputParticleLAB;
203
204 productionModelData.incidentParticleSCM =
205 productionModelData.incidentParticleLAB;
206 productionModelData.nucleusParticleSCM =
207 productionModelData.nucleusParticleLAB;
208 productionModelData.outputParticleSCM =
209 productionModelData.outputParticleLAB;
210 productionModelData.nucleusOutputParticleSCM =
211 productionModelData.nucleusOutputParticleLAB;
212
213 productionModelData.incidentParticleSCM.boost( -boostVec );
214 productionModelData.nucleusParticleSCM.boost( -boostVec );
215 productionModelData.outputParticleSCM.boost( -boostVec );
216 productionModelData.nucleusOutputParticleSCM.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
238 hasMassCutTriggered =
239 std::pow( ( outputParticleMass - massCutOPCenter ) *
240 cosMassCutEllipseAngle +
241 ( nucleusOutputParticleMass - massCutNOPCenter ) *
242 sinMassCutEllipseAngle, 2 ) / massCutOPWidth2 +
243 std::pow( - ( outputParticleMass - massCutOPCenter ) *
244 sinMassCutEllipseAngle +
245 ( nucleusOutputParticleMass - massCutNOPCenter ) *
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
258 if ( absorbedEnergyCutCLWidth <= 0. || absorbedEnergyCutCRWidth <= 0. )
259 {
260 hasAbsorbedEnergyCutTriggered = false;
261 }
262 else
263 {
264 G4double absorbedEnergyCutCLWidth2(
265 absorbedEnergyCutCLWidth * absorbedEnergyCutCLWidth );
266 G4double absorbedEnergyCutCRWidth2(
267 absorbedEnergyCutCRWidth * absorbedEnergyCutCRWidth );
268
269 hasAbsorbedEnergyCutTriggered =
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
288 G4bool CexmcChargeExchangeReconstructor::HasFullTrigger( void ) const
289 {
290 if ( ! hasBasicTrigger )
291 return false;
292 if ( useMassCut && ! hasMassCutTriggered )
293 return false;
294 if ( useAbsorbedEnergyCut && ! hasAbsorbedEnergyCutTriggered )
295 return false;
296
297 return true;
298 }
299
300
301 void CexmcChargeExchangeReconstructor::SetExpectedMomentumAmpDiff(
302 G4double value )
303 {
304 expectedMomentumAmp = particleGun->GetOrigMomentumAmp() + value;
305 }
306
307