| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 /* 26 /*
27 * =========================================== 27 * ============================================================================
28 * 28 *
29 * Filename: CexmcChargeExchangeReconst 29 * Filename: CexmcChargeExchangeReconstructor.cc
30 * 30 *
31 * Description: charge exchange reconstruc 31 * Description: charge exchange reconstructor
32 * 32 *
33 * Version: 1.0 33 * Version: 1.0
34 * Created: 02.12.2009 15:17:13 34 * Created: 02.12.2009 15:17:13
35 * Revision: none 35 * Revision: none
36 * Compiler: gcc 36 * Compiler: gcc
37 * 37 *
38 * Author: Alexey Radkov (), 38 * Author: Alexey Radkov (),
39 * Company: PNPI 39 * Company: PNPI
40 * 40 *
41 * =========================================== 41 * ============================================================================
42 */ 42 */
43 43
44 #include <cmath> 44 #include <cmath>
45 #include <G4ThreeVector.hh> 45 #include <G4ThreeVector.hh>
46 #include <G4LorentzVector.hh> 46 #include <G4LorentzVector.hh>
47 #include "CexmcChargeExchangeReconstructor.hh" 47 #include "CexmcChargeExchangeReconstructor.hh"
48 #include "CexmcChargeExchangeReconstructorMess 48 #include "CexmcChargeExchangeReconstructorMessenger.hh"
49 #include "CexmcEnergyDepositStore.hh" 49 #include "CexmcEnergyDepositStore.hh"
50 #include "CexmcPrimaryGeneratorAction.hh" 50 #include "CexmcPrimaryGeneratorAction.hh"
51 #include "CexmcParticleGun.hh" 51 #include "CexmcParticleGun.hh"
52 #include "CexmcProductionModel.hh" 52 #include "CexmcProductionModel.hh"
53 #include "CexmcRunManager.hh" 53 #include "CexmcRunManager.hh"
54 #include "CexmcException.hh" 54 #include "CexmcException.hh"
55 55
56 56
57 CexmcChargeExchangeReconstructor::CexmcChargeE 57 CexmcChargeExchangeReconstructor::CexmcChargeExchangeReconstructor(
58 const CexmcProduct 58 const CexmcProductionModel * productionModel ) :
59 outputParticleMass( 0 ), nucleusOutputPart 59 outputParticleMass( 0 ), nucleusOutputParticleMass( 0 ),
60 useTableMass( false ), useMassCut( false ) 60 useTableMass( false ), useMassCut( false ), massCutOPCenter( 0 ),
61 massCutNOPCenter( 0 ), massCutOPWidth( 0 ) 61 massCutNOPCenter( 0 ), massCutOPWidth( 0 ), massCutNOPWidth( 0 ),
62 massCutEllipseAngle( 0 ), useAbsorbedEnerg 62 massCutEllipseAngle( 0 ), useAbsorbedEnergyCut( false ),
63 absorbedEnergyCutCLCenter( 0 ), absorbedEn 63 absorbedEnergyCutCLCenter( 0 ), absorbedEnergyCutCRCenter( 0 ),
64 absorbedEnergyCutCLWidth( 0 ), absorbedEne 64 absorbedEnergyCutCLWidth( 0 ), absorbedEnergyCutCRWidth( 0 ),
65 absorbedEnergyCutEllipseAngle( 0 ), expect 65 absorbedEnergyCutEllipseAngle( 0 ), expectedMomentumAmp( -1 ),
66 edCollectionAlgorithm( CexmcCollectEDInAll 66 edCollectionAlgorithm( CexmcCollectEDInAllCrystals ),
67 hasMassCutTriggered( false ), hasAbsorbedE 67 hasMassCutTriggered( false ), hasAbsorbedEnergyCutTriggered( false ),
68 beamParticleIsInitialized( false ), partic 68 beamParticleIsInitialized( false ), particleGun( NULL ), messenger( NULL )
69 { 69 {
70 if ( ! productionModel ) 70 if ( ! productionModel )
71 throw CexmcException( CexmcWeirdExcept 71 throw CexmcException( CexmcWeirdException );
72 72
73 productionModelData.incidentParticle = 73 productionModelData.incidentParticle =
74 production 74 productionModel->GetIncidentParticle();
75 75
76 CexmcRunManager * runManager( static_cast 76 CexmcRunManager * runManager( static_cast< CexmcRunManager * >(
77 G4 77 G4RunManager::GetRunManager() ) );
78 const CexmcPrimaryGeneratorAction * prima 78 const CexmcPrimaryGeneratorAction * primaryGeneratorAction(
79 static_cast< const Cex 79 static_cast< const CexmcPrimaryGeneratorAction * >(
80 runManager->Ge 80 runManager->GetUserPrimaryGeneratorAction() ) );
81 CexmcPrimaryGeneratorAction * thePrimaryG 81 CexmcPrimaryGeneratorAction * thePrimaryGeneratorAction(
82 const_cast< CexmcPrima 82 const_cast< CexmcPrimaryGeneratorAction * >(
83 primaryGenerat 83 primaryGeneratorAction ) );
84 particleGun = thePrimaryGeneratorAction->G 84 particleGun = thePrimaryGeneratorAction->GetParticleGun();
85 85
86 productionModelData.nucleusParticle = 86 productionModelData.nucleusParticle =
87 production 87 productionModel->GetNucleusParticle();
88 productionModelData.outputParticle = 88 productionModelData.outputParticle =
89 production 89 productionModel->GetOutputParticle();
90 productionModelData.nucleusOutputParticle 90 productionModelData.nucleusOutputParticle =
91 production 91 productionModel->GetNucleusOutputParticle();
92 92
93 messenger = new CexmcChargeExchangeReconst 93 messenger = new CexmcChargeExchangeReconstructorMessenger( this );
94 } 94 }
95 95
96 96
97 CexmcChargeExchangeReconstructor::~CexmcCharge 97 CexmcChargeExchangeReconstructor::~CexmcChargeExchangeReconstructor()
98 { 98 {
99 delete messenger; 99 delete messenger;
100 } 100 }
101 101
102 102
103 void CexmcChargeExchangeReconstructor::SetupB 103 void CexmcChargeExchangeReconstructor::SetupBeamParticle( void )
104 { 104 {
105 if ( *productionModelData.incidentParticle 105 if ( *productionModelData.incidentParticle !=
106 *parti 106 *particleGun->GetParticleDefinition() )
107 throw CexmcException( CexmcBeamAndInci 107 throw CexmcException( CexmcBeamAndIncidentParticlesMismatch );
108 108
109 beamParticleIsInitialized = true; 109 beamParticleIsInitialized = true;
110 } 110 }
111 111
112 112
113 void CexmcChargeExchangeReconstructor::Recons 113 void CexmcChargeExchangeReconstructor::Reconstruct(
114 const Cexm 114 const CexmcEnergyDepositStore * edStore )
115 { 115 {
116 if ( ! beamParticleIsInitialized ) 116 if ( ! beamParticleIsInitialized )
117 { 117 {
118 if ( *productionModelData.incidentPart 118 if ( *productionModelData.incidentParticle !=
119 *parti 119 *particleGun->GetParticleDefinition() )
120 throw CexmcException( CexmcBeamAnd 120 throw CexmcException( CexmcBeamAndIncidentParticlesMismatch );
121 121
122 beamParticleIsInitialized = true; 122 beamParticleIsInitialized = true;
123 } 123 }
124 124
125 if ( edCollectionAlgorithm == CexmcCollect 125 if ( edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
126 collectEDInAdjacentCrystals = true; 126 collectEDInAdjacentCrystals = true;
127 127
128 ReconstructEntryPoints( edStore ); 128 ReconstructEntryPoints( edStore );
129 if ( hasBasicTrigger ) 129 if ( hasBasicTrigger )
130 ReconstructTargetPoint(); 130 ReconstructTargetPoint();
131 if ( hasBasicTrigger ) 131 if ( hasBasicTrigger )
132 ReconstructAngle(); 132 ReconstructAngle();
133 133
134 G4ThreeVector epLeft( calorimeterEPLeftWo 134 G4ThreeVector epLeft( calorimeterEPLeftWorldPosition -
135 targetEPWorldPositi 135 targetEPWorldPosition );
136 G4ThreeVector epRight( calorimeterEPRight 136 G4ThreeVector epRight( calorimeterEPRightWorldPosition -
137 targetEPWorldPosit 137 targetEPWorldPosition );
138 138
139 G4double cosTheAngle( std::cos( theAngle 139 G4double cosTheAngle( std::cos( theAngle ) );
140 G4double calorimeterEDLeft( edStore->calo 140 G4double calorimeterEDLeft( edStore->calorimeterEDLeft );
141 G4double calorimeterEDRight( edStore->cal 141 G4double calorimeterEDRight( edStore->calorimeterEDRight );
142 142
143 if ( edCollectionAlgorithm == CexmcCollect 143 if ( edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
144 { 144 {
145 calorimeterEDLeft = calorimeterEDLeftA 145 calorimeterEDLeft = calorimeterEDLeftAdjacent;
146 calorimeterEDRight = calorimeterEDRigh 146 calorimeterEDRight = calorimeterEDRightAdjacent;
147 } 147 }
148 148
149 //G4double cosOutputParticleLAB( 149 //G4double cosOutputParticleLAB(
150 //( calorimeterEDLeft * cosAngleLeft + 150 //( calorimeterEDLeft * cosAngleLeft +
151 //calorimeterEDRight * cosAngleRight 151 //calorimeterEDRight * cosAngleRight ) /
152 //std::sqrt( calorimeterEDLeft * cal 152 //std::sqrt( calorimeterEDLeft * calorimeterEDLeft +
153 //calorimeterEDRight * ca 153 //calorimeterEDRight * calorimeterEDRight +
154 //calorimeterEDLeft * cal 154 //calorimeterEDLeft * calorimeterEDRight * cosTheAngle ) );
155 155
156 outputParticleMass = std::sqrt( 2 * calori 156 outputParticleMass = std::sqrt( 2 * calorimeterEDLeft *
157 calorimete 157 calorimeterEDRight * ( 1 - cosTheAngle ) );
158 158
159 G4ThreeVector opdpLeftMomentum( epLeft 159 G4ThreeVector opdpLeftMomentum( epLeft );
160 opdpLeftMomentum.setMag( calorimeterEDLeft 160 opdpLeftMomentum.setMag( calorimeterEDLeft );
161 G4ThreeVector opdpRightMomentum( epRigh 161 G4ThreeVector opdpRightMomentum( epRight );
162 opdpRightMomentum.setMag( calorimeterEDRig 162 opdpRightMomentum.setMag( calorimeterEDRight );
163 G4ThreeVector opMomentum( opdpLeftMomen 163 G4ThreeVector opMomentum( opdpLeftMomentum + opdpRightMomentum );
164 164
165 /* opMass will be used only in calculation 165 /* opMass will be used only in calculation of output particle's total
166 * energy, in other places outputParticleM 166 * energy, in other places outputParticleMass should be used instead */
167 G4double opMass( useTableMass ? 167 G4double opMass( useTableMass ?
168 productionModelDa 168 productionModelData.outputParticle->GetPDGMass() :
169 outputParticleMas 169 outputParticleMass );
170 /* the formula below is equivalent to 170 /* the formula below is equivalent to
171 * calorimeterEDLeft + calorimeterEDRight 171 * calorimeterEDLeft + calorimeterEDRight if opMass = outputParticleMass */
172 G4double opEnergy( std::sqrt( opMo 172 G4double opEnergy( std::sqrt( opMomentum.mag2() +
173 opMa 173 opMass * opMass ) );
174 productionModelData.outputParticleLAB = G4 174 productionModelData.outputParticleLAB = G4LorentzVector( opMomentum,
175 175 opEnergy );
176 176
177 G4ThreeVector incidentParticleMomentum( p 177 G4ThreeVector incidentParticleMomentum( particleGun->GetOrigDirection() );
178 G4double incidentParticleMomentumAmp 178 G4double incidentParticleMomentumAmp( expectedMomentumAmp > 0 ?
179 expectedMomentumAmp : part 179 expectedMomentumAmp : particleGun->GetOrigMomentumAmp() );
180 incidentParticleMomentum *= incidentPartic 180 incidentParticleMomentum *= incidentParticleMomentumAmp;
181 181
182 G4double incidentParticlePDGMass( 182 G4double incidentParticlePDGMass(
183 productionModelData.in 183 productionModelData.incidentParticle->GetPDGMass() );
184 G4double incidentParticlePDGMass2( i 184 G4double incidentParticlePDGMass2( incidentParticlePDGMass *
185 i 185 incidentParticlePDGMass );
186 G4double incidentParticleEnergy( 186 G4double incidentParticleEnergy(
187 std::sqrt( incidentParticleMomentumAmp 187 std::sqrt( incidentParticleMomentumAmp * incidentParticleMomentumAmp +
188 incidentParticlePDGMass2 ) 188 incidentParticlePDGMass2 ) );
189 189
190 productionModelData.incidentParticleLAB = 190 productionModelData.incidentParticleLAB = G4LorentzVector(
191 incidentParticleMoment 191 incidentParticleMomentum, incidentParticleEnergy );
192 G4double nucleusParticlePDGMass( 192 G4double nucleusParticlePDGMass(
193 productionModelData.nu 193 productionModelData.nucleusParticle->GetPDGMass() );
194 productionModelData.nucleusParticleLAB = G 194 productionModelData.nucleusParticleLAB = G4LorentzVector(
195 G4ThreeVector( 0, 0, 0 195 G4ThreeVector( 0, 0, 0 ), nucleusParticlePDGMass );
196 196
197 G4LorentzVector lVecSum( productionModelD 197 G4LorentzVector lVecSum( productionModelData.incidentParticleLAB +
198 productionModelData.nu 198 productionModelData.nucleusParticleLAB );
199 G4ThreeVector boostVec( lVecSum.boostVe 199 G4ThreeVector boostVec( lVecSum.boostVector() );
200 200
201 productionModelData.nucleusOutputParticleL 201 productionModelData.nucleusOutputParticleLAB =
202 lVecSum - productionModelData.outp 202 lVecSum - productionModelData.outputParticleLAB;
203 203
204 productionModelData.incidentParticleSCM = 204 productionModelData.incidentParticleSCM =
205 productionModelData.incidentPartic 205 productionModelData.incidentParticleLAB;
206 productionModelData.nucleusParticleSCM = 206 productionModelData.nucleusParticleSCM =
207 productionModelData.nucleusParticl 207 productionModelData.nucleusParticleLAB;
208 productionModelData.outputParticleSCM = 208 productionModelData.outputParticleSCM =
209 productionModelData.outputParticle 209 productionModelData.outputParticleLAB;
210 productionModelData.nucleusOutputParticleS 210 productionModelData.nucleusOutputParticleSCM =
211 productionModelData.nucleusOutputP 211 productionModelData.nucleusOutputParticleLAB;
212 212
213 productionModelData.incidentParticleSCM.bo 213 productionModelData.incidentParticleSCM.boost( -boostVec );
214 productionModelData.nucleusParticleSCM.boo 214 productionModelData.nucleusParticleSCM.boost( -boostVec );
215 productionModelData.outputParticleSCM.boos 215 productionModelData.outputParticleSCM.boost( -boostVec );
216 productionModelData.nucleusOutputParticleS 216 productionModelData.nucleusOutputParticleSCM.boost( -boostVec );
217 217
218 G4ThreeVector nopMomentum( incidentPartic 218 G4ThreeVector nopMomentum( incidentParticleMomentum - opMomentum );
219 G4double nopEnergy( incidentParticle 219 G4double nopEnergy( incidentParticleEnergy + nucleusParticlePDGMass -
220 opEnergy ); 220 opEnergy );
221 nucleusOutputParticleMass = std::sqrt( nop 221 nucleusOutputParticleMass = std::sqrt( nopEnergy * nopEnergy -
222 nop 222 nopMomentum.mag2() );
223 223
224 if ( useMassCut ) 224 if ( useMassCut )
225 { 225 {
226 G4double cosMassCutEllipseAngle( std: 226 G4double cosMassCutEllipseAngle( std::cos( massCutEllipseAngle ) );
227 G4double sinMassCutEllipseAngle( std: 227 G4double sinMassCutEllipseAngle( std::sin( massCutEllipseAngle ) );
228 228
229 if ( massCutOPWidth <= 0. || massCutNO 229 if ( massCutOPWidth <= 0. || massCutNOPWidth <= 0. )
230 { 230 {
231 hasMassCutTriggered = false; 231 hasMassCutTriggered = false;
232 } 232 }
233 else 233 else
234 { 234 {
235 G4double massCutOPWidth2( massCut 235 G4double massCutOPWidth2( massCutOPWidth * massCutOPWidth );
236 G4double massCutNOPWidth2( massCu 236 G4double massCutNOPWidth2( massCutNOPWidth * massCutNOPWidth );
237 237
238 hasMassCutTriggered = 238 hasMassCutTriggered =
239 std::pow( ( outputParticleMass 239 std::pow( ( outputParticleMass - massCutOPCenter ) *
240 cosMassCutEllips 240 cosMassCutEllipseAngle +
241 ( nucleusOutputParti 241 ( nucleusOutputParticleMass - massCutNOPCenter ) *
242 sinMassCutEllips 242 sinMassCutEllipseAngle, 2 ) / massCutOPWidth2 +
243 std::pow( - ( outputParticleMa 243 std::pow( - ( outputParticleMass - massCutOPCenter ) *
244 sinMassCutEllips 244 sinMassCutEllipseAngle +
245 ( nucleusOutputParti 245 ( nucleusOutputParticleMass - massCutNOPCenter ) *
246 cosMassCutEllips 246 cosMassCutEllipseAngle, 2 ) / massCutNOPWidth2 <
247 1; 247 1;
248 } 248 }
249 } 249 }
250 250
251 if ( useAbsorbedEnergyCut ) 251 if ( useAbsorbedEnergyCut )
252 { 252 {
253 G4double cosAbsorbedEnergyCutEllipseA 253 G4double cosAbsorbedEnergyCutEllipseAngle(
254 std::cos( abso 254 std::cos( absorbedEnergyCutEllipseAngle ) );
255 G4double sinAbsorbedEnergyCutEllipseA 255 G4double sinAbsorbedEnergyCutEllipseAngle(
256 std::sin( abso 256 std::sin( absorbedEnergyCutEllipseAngle ) );
257 257
258 if ( absorbedEnergyCutCLWidth <= 0. || 258 if ( absorbedEnergyCutCLWidth <= 0. || absorbedEnergyCutCRWidth <= 0. )
259 { 259 {
260 hasAbsorbedEnergyCutTriggered = fa 260 hasAbsorbedEnergyCutTriggered = false;
261 } 261 }
262 else 262 else
263 { 263 {
264 G4double absorbedEnergyCutCLWidth 264 G4double absorbedEnergyCutCLWidth2(
265 absorbedEnergyCutCLWid 265 absorbedEnergyCutCLWidth * absorbedEnergyCutCLWidth );
266 G4double absorbedEnergyCutCRWidth 266 G4double absorbedEnergyCutCRWidth2(
267 absorbedEnergyCutCRWid 267 absorbedEnergyCutCRWidth * absorbedEnergyCutCRWidth );
268 268
269 hasAbsorbedEnergyCutTriggered = 269 hasAbsorbedEnergyCutTriggered =
270 std::pow( ( calorimeterEDLeft 270 std::pow( ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) *
271 cosAbsorbedEnerg 271 cosAbsorbedEnergyCutEllipseAngle +
272 ( calorimeterEDRight 272 ( calorimeterEDRight - absorbedEnergyCutCRCenter ) *
273 sinAbsorbedEnerg 273 sinAbsorbedEnergyCutEllipseAngle, 2 ) /
274 absorbedEnergyCutCLWidth2 + 274 absorbedEnergyCutCLWidth2 +
275 std::pow( - ( calorimeterEDLef 275 std::pow( - ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) *
276 sinAbsorbedEnerg 276 sinAbsorbedEnergyCutEllipseAngle +
277 ( calorimeterEDRight 277 ( calorimeterEDRight - absorbedEnergyCutCRCenter ) *
278 cosAbsorbedEnerg 278 cosAbsorbedEnergyCutEllipseAngle, 2 ) /
279 absorbedEnergyCutCRWidth2 < 279 absorbedEnergyCutCRWidth2 <
280 1; 280 1;
281 } 281 }
282 } 282 }
283 283
284 hasBasicTrigger = true; 284 hasBasicTrigger = true;
285 } 285 }
286 286
287 287
288 G4bool CexmcChargeExchangeReconstructor::HasF 288 G4bool CexmcChargeExchangeReconstructor::HasFullTrigger( void ) const
289 { 289 {
290 if ( ! hasBasicTrigger ) 290 if ( ! hasBasicTrigger )
291 return false; 291 return false;
292 if ( useMassCut && ! hasMassCutTriggered ) 292 if ( useMassCut && ! hasMassCutTriggered )
293 return false; 293 return false;
294 if ( useAbsorbedEnergyCut && ! hasAbsorbed 294 if ( useAbsorbedEnergyCut && ! hasAbsorbedEnergyCutTriggered )
295 return false; 295 return false;
296 296
297 return true; 297 return true;
298 } 298 }
299 299
300 300
301 void CexmcChargeExchangeReconstructor::SetExp 301 void CexmcChargeExchangeReconstructor::SetExpectedMomentumAmpDiff(
302 302 G4double value )
303 { 303 {
304 expectedMomentumAmp = particleGun->GetOrig 304 expectedMomentumAmp = particleGun->GetOrigMomentumAmp() + value;
305 } 305 }
306 306
307 307