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> << 45 #include <G4ThreeVector.hh> 44 #include <G4ThreeVector.hh> 46 #include <G4LorentzVector.hh> 45 #include <G4LorentzVector.hh> 47 #include "CexmcChargeExchangeReconstructor.hh" 46 #include "CexmcChargeExchangeReconstructor.hh" 48 #include "CexmcChargeExchangeReconstructorMess 47 #include "CexmcChargeExchangeReconstructorMessenger.hh" 49 #include "CexmcEnergyDepositStore.hh" 48 #include "CexmcEnergyDepositStore.hh" 50 #include "CexmcPrimaryGeneratorAction.hh" 49 #include "CexmcPrimaryGeneratorAction.hh" 51 #include "CexmcParticleGun.hh" 50 #include "CexmcParticleGun.hh" 52 #include "CexmcProductionModel.hh" 51 #include "CexmcProductionModel.hh" 53 #include "CexmcRunManager.hh" 52 #include "CexmcRunManager.hh" 54 #include "CexmcException.hh" 53 #include "CexmcException.hh" >> 54 #include "CexmcCommon.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 ), hasMassCutTriggered( false ), 66 edCollectionAlgorithm( CexmcCollectEDInAll << 66 hasAbsorbedEnergyCutTriggered( false ), beamParticleIsInitialized( false ), 67 hasMassCutTriggered( false ), hasAbsorbedE << 67 particleGun( NULL ), messenger( NULL ) 68 beamParticleIsInitialized( false ), partic << 69 { 68 { 70 if ( ! productionModel ) 69 if ( ! productionModel ) 71 throw CexmcException( CexmcWeirdExcept 70 throw CexmcException( CexmcWeirdException ); 72 71 73 productionModelData.incidentParticle = 72 productionModelData.incidentParticle = 74 production 73 productionModel->GetIncidentParticle(); 75 74 76 CexmcRunManager * runManager( static_cast 75 CexmcRunManager * runManager( static_cast< CexmcRunManager * >( 77 G4 76 G4RunManager::GetRunManager() ) ); 78 const CexmcPrimaryGeneratorAction * prima 77 const CexmcPrimaryGeneratorAction * primaryGeneratorAction( 79 static_cast< const Cex 78 static_cast< const CexmcPrimaryGeneratorAction * >( 80 runManager->Ge 79 runManager->GetUserPrimaryGeneratorAction() ) ); 81 CexmcPrimaryGeneratorAction * thePrimaryG 80 CexmcPrimaryGeneratorAction * thePrimaryGeneratorAction( 82 const_cast< CexmcPrima 81 const_cast< CexmcPrimaryGeneratorAction * >( 83 primaryGenerat 82 primaryGeneratorAction ) ); 84 particleGun = thePrimaryGeneratorAction->G 83 particleGun = thePrimaryGeneratorAction->GetParticleGun(); 85 84 86 productionModelData.nucleusParticle = 85 productionModelData.nucleusParticle = 87 production 86 productionModel->GetNucleusParticle(); 88 productionModelData.outputParticle = 87 productionModelData.outputParticle = 89 production 88 productionModel->GetOutputParticle(); 90 productionModelData.nucleusOutputParticle 89 productionModelData.nucleusOutputParticle = 91 production 90 productionModel->GetNucleusOutputParticle(); 92 91 93 messenger = new CexmcChargeExchangeReconst 92 messenger = new CexmcChargeExchangeReconstructorMessenger( this ); 94 } 93 } 95 94 96 95 97 CexmcChargeExchangeReconstructor::~CexmcCharge 96 CexmcChargeExchangeReconstructor::~CexmcChargeExchangeReconstructor() 98 { 97 { 99 delete messenger; 98 delete messenger; 100 } 99 } 101 100 102 101 103 void CexmcChargeExchangeReconstructor::SetupB 102 void CexmcChargeExchangeReconstructor::SetupBeamParticle( void ) 104 { 103 { 105 if ( *productionModelData.incidentParticle 104 if ( *productionModelData.incidentParticle != 106 *parti 105 *particleGun->GetParticleDefinition() ) 107 throw CexmcException( CexmcBeamAndInci 106 throw CexmcException( CexmcBeamAndIncidentParticlesMismatch ); 108 107 109 beamParticleIsInitialized = true; 108 beamParticleIsInitialized = true; 110 } 109 } 111 110 112 111 113 void CexmcChargeExchangeReconstructor::Recons 112 void CexmcChargeExchangeReconstructor::Reconstruct( 114 const Cexm 113 const CexmcEnergyDepositStore * edStore ) 115 { 114 { 116 if ( ! beamParticleIsInitialized ) 115 if ( ! beamParticleIsInitialized ) 117 { 116 { 118 if ( *productionModelData.incidentPart 117 if ( *productionModelData.incidentParticle != 119 *parti 118 *particleGun->GetParticleDefinition() ) 120 throw CexmcException( CexmcBeamAnd 119 throw CexmcException( CexmcBeamAndIncidentParticlesMismatch ); 121 120 122 beamParticleIsInitialized = true; 121 beamParticleIsInitialized = true; 123 } 122 } 124 123 125 if ( edCollectionAlgorithm == CexmcCollect << 126 collectEDInAdjacentCrystals = true; << 127 << 128 ReconstructEntryPoints( edStore ); 124 ReconstructEntryPoints( edStore ); 129 if ( hasBasicTrigger ) 125 if ( hasBasicTrigger ) 130 ReconstructTargetPoint(); 126 ReconstructTargetPoint(); 131 if ( hasBasicTrigger ) 127 if ( hasBasicTrigger ) 132 ReconstructAngle(); 128 ReconstructAngle(); 133 129 134 G4ThreeVector epLeft( calorimeterEPLeftWo 130 G4ThreeVector epLeft( calorimeterEPLeftWorldPosition - 135 targetEPWorldPositi 131 targetEPWorldPosition ); 136 G4ThreeVector epRight( calorimeterEPRight 132 G4ThreeVector epRight( calorimeterEPRightWorldPosition - 137 targetEPWorldPosit 133 targetEPWorldPosition ); 138 134 139 G4double cosTheAngle( std::cos( theAngle 135 G4double cosTheAngle( std::cos( theAngle ) ); 140 G4double calorimeterEDLeft( edStore->calo 136 G4double calorimeterEDLeft( edStore->calorimeterEDLeft ); 141 G4double calorimeterEDRight( edStore->cal 137 G4double calorimeterEDRight( edStore->calorimeterEDRight ); 142 138 143 if ( edCollectionAlgorithm == CexmcCollect << 144 { << 145 calorimeterEDLeft = calorimeterEDLeftA << 146 calorimeterEDRight = calorimeterEDRigh << 147 } << 148 << 149 //G4double cosOutputParticleLAB( 139 //G4double cosOutputParticleLAB( 150 //( calorimeterEDLeft * cosAngleLeft + 140 //( calorimeterEDLeft * cosAngleLeft + 151 //calorimeterEDRight * cosAngleRight 141 //calorimeterEDRight * cosAngleRight ) / 152 //std::sqrt( calorimeterEDLeft * cal 142 //std::sqrt( calorimeterEDLeft * calorimeterEDLeft + 153 //calorimeterEDRight * ca 143 //calorimeterEDRight * calorimeterEDRight + 154 //calorimeterEDLeft * cal 144 //calorimeterEDLeft * calorimeterEDRight * cosTheAngle ) ); 155 145 156 outputParticleMass = std::sqrt( 2 * calori 146 outputParticleMass = std::sqrt( 2 * calorimeterEDLeft * 157 calorimete 147 calorimeterEDRight * ( 1 - cosTheAngle ) ); 158 148 159 G4ThreeVector opdpLeftMomentum( epLeft 149 G4ThreeVector opdpLeftMomentum( epLeft ); 160 opdpLeftMomentum.setMag( calorimeterEDLeft 150 opdpLeftMomentum.setMag( calorimeterEDLeft ); 161 G4ThreeVector opdpRightMomentum( epRigh 151 G4ThreeVector opdpRightMomentum( epRight ); 162 opdpRightMomentum.setMag( calorimeterEDRig 152 opdpRightMomentum.setMag( calorimeterEDRight ); 163 G4ThreeVector opMomentum( opdpLeftMomen 153 G4ThreeVector opMomentum( opdpLeftMomentum + opdpRightMomentum ); 164 154 165 /* opMass will be used only in calculation << 166 * energy, in other places outputParticleM << 167 G4double opMass( useTableMass ? 155 G4double opMass( useTableMass ? 168 productionModelDa 156 productionModelData.outputParticle->GetPDGMass() : 169 outputParticleMas 157 outputParticleMass ); 170 /* the formula below is equivalent to << 158 G4double opEnergy( std::sqrt( 171 * calorimeterEDLeft + calorimeterEDRight << 159 opMomentum.mag2() + opMass * opMass ) ); 172 G4double opEnergy( std::sqrt( opMo << 173 opMa << 174 productionModelData.outputParticleLAB = G4 160 productionModelData.outputParticleLAB = G4LorentzVector( opMomentum, 175 161 opEnergy ); 176 162 177 G4ThreeVector incidentParticleMomentum( p 163 G4ThreeVector incidentParticleMomentum( particleGun->GetOrigDirection() ); 178 G4double incidentParticleMomentumAmp << 164 G4double incidentParticleMomentumAmp( 179 expectedMomentumAmp : part << 165 particleGun->GetOrigMomentumAmp() ); 180 incidentParticleMomentum *= incidentPartic 166 incidentParticleMomentum *= incidentParticleMomentumAmp; 181 167 182 G4double incidentParticlePDGMass( 168 G4double incidentParticlePDGMass( 183 productionModelData.in 169 productionModelData.incidentParticle->GetPDGMass() ); 184 G4double incidentParticlePDGMass2( i 170 G4double incidentParticlePDGMass2( incidentParticlePDGMass * 185 i 171 incidentParticlePDGMass ); 186 G4double incidentParticleEnergy( 172 G4double incidentParticleEnergy( 187 std::sqrt( incidentParticleMomentumAmp 173 std::sqrt( incidentParticleMomentumAmp * incidentParticleMomentumAmp + 188 incidentParticlePDGMass2 ) 174 incidentParticlePDGMass2 ) ); 189 175 190 productionModelData.incidentParticleLAB = 176 productionModelData.incidentParticleLAB = G4LorentzVector( 191 incidentParticleMoment 177 incidentParticleMomentum, incidentParticleEnergy ); 192 G4double nucleusParticlePDGMass( 178 G4double nucleusParticlePDGMass( 193 productionModelData.nu 179 productionModelData.nucleusParticle->GetPDGMass() ); 194 productionModelData.nucleusParticleLAB = G 180 productionModelData.nucleusParticleLAB = G4LorentzVector( 195 G4ThreeVector( 0, 0, 0 181 G4ThreeVector( 0, 0, 0 ), nucleusParticlePDGMass ); 196 << 197 G4LorentzVector lVecSum( productionModelD 182 G4LorentzVector lVecSum( productionModelData.incidentParticleLAB + 198 productionModelData.nu 183 productionModelData.nucleusParticleLAB ); 199 G4ThreeVector boostVec( lVecSum.boostVe 184 G4ThreeVector boostVec( lVecSum.boostVector() ); 200 185 201 productionModelData.nucleusOutputParticleL 186 productionModelData.nucleusOutputParticleLAB = 202 lVecSum - productionModelData.outp 187 lVecSum - productionModelData.outputParticleLAB; 203 188 204 productionModelData.incidentParticleSCM = 189 productionModelData.incidentParticleSCM = 205 productionModelData.incidentPartic 190 productionModelData.incidentParticleLAB; 206 productionModelData.nucleusParticleSCM = 191 productionModelData.nucleusParticleSCM = 207 productionModelData.nucleusParticl 192 productionModelData.nucleusParticleLAB; 208 productionModelData.outputParticleSCM = 193 productionModelData.outputParticleSCM = 209 productionModelData.outputParticle 194 productionModelData.outputParticleLAB; 210 productionModelData.nucleusOutputParticleS 195 productionModelData.nucleusOutputParticleSCM = 211 productionModelData.nucleusOutputP 196 productionModelData.nucleusOutputParticleLAB; 212 197 213 productionModelData.incidentParticleSCM.bo 198 productionModelData.incidentParticleSCM.boost( -boostVec ); 214 productionModelData.nucleusParticleSCM.boo 199 productionModelData.nucleusParticleSCM.boost( -boostVec ); 215 productionModelData.outputParticleSCM.boos 200 productionModelData.outputParticleSCM.boost( -boostVec ); 216 productionModelData.nucleusOutputParticleS 201 productionModelData.nucleusOutputParticleSCM.boost( -boostVec ); 217 202 >> 203 G4double edDelta2( >> 204 std::pow( ( calorimeterEDLeft - calorimeterEDRight ) / >> 205 ( calorimeterEDLeft + calorimeterEDRight ), >> 206 2 ) ); >> 207 G4double outputParticleKinEnergy( >> 208 std::sqrt( 2 * opMass * opMass / ( 1 - cosTheAngle ) / >> 209 ( 1 - edDelta2 ) ) - opMass ); 218 G4ThreeVector nopMomentum( incidentPartic 210 G4ThreeVector nopMomentum( incidentParticleMomentum - opMomentum ); 219 G4double nopEnergy( incidentParticle << 211 G4double nopEnergy( 220 opEnergy ); << 212 std::sqrt( incidentParticleMomentum.mag2() + >> 213 incidentParticlePDGMass2 ) + >> 214 nucleusParticlePDGMass - >> 215 ( outputParticleKinEnergy + opMass ) ); 221 nucleusOutputParticleMass = std::sqrt( nop 216 nucleusOutputParticleMass = std::sqrt( nopEnergy * nopEnergy - 222 nop 217 nopMomentum.mag2() ); 223 218 224 if ( useMassCut ) 219 if ( useMassCut ) 225 { 220 { 226 G4double cosMassCutEllipseAngle( std: 221 G4double cosMassCutEllipseAngle( std::cos( massCutEllipseAngle ) ); 227 G4double sinMassCutEllipseAngle( std: 222 G4double sinMassCutEllipseAngle( std::sin( massCutEllipseAngle ) ); 228 223 229 if ( massCutOPWidth <= 0. || massCutNO 224 if ( massCutOPWidth <= 0. || massCutNOPWidth <= 0. ) 230 { 225 { 231 hasMassCutTriggered = false; 226 hasMassCutTriggered = false; 232 } 227 } 233 else 228 else 234 { 229 { 235 G4double massCutOPWidth2( massCut 230 G4double massCutOPWidth2( massCutOPWidth * massCutOPWidth ); 236 G4double massCutNOPWidth2( massCu 231 G4double massCutNOPWidth2( massCutNOPWidth * massCutNOPWidth ); 237 232 238 hasMassCutTriggered = 233 hasMassCutTriggered = 239 std::pow( ( outputParticleMass 234 std::pow( ( outputParticleMass - massCutOPCenter ) * 240 cosMassCutEllips 235 cosMassCutEllipseAngle + 241 ( nucleusOutputParti 236 ( nucleusOutputParticleMass - massCutNOPCenter ) * 242 sinMassCutEllips 237 sinMassCutEllipseAngle, 2 ) / massCutOPWidth2 + 243 std::pow( - ( outputParticleMa 238 std::pow( - ( outputParticleMass - massCutOPCenter ) * 244 sinMassCutEllips 239 sinMassCutEllipseAngle + 245 ( nucleusOutputParti 240 ( nucleusOutputParticleMass - massCutNOPCenter ) * 246 cosMassCutEllips 241 cosMassCutEllipseAngle, 2 ) / massCutNOPWidth2 < 247 1; 242 1; 248 } 243 } 249 } 244 } 250 245 251 if ( useAbsorbedEnergyCut ) 246 if ( useAbsorbedEnergyCut ) 252 { 247 { 253 G4double cosAbsorbedEnergyCutEllipseA 248 G4double cosAbsorbedEnergyCutEllipseAngle( 254 std::cos( abso 249 std::cos( absorbedEnergyCutEllipseAngle ) ); 255 G4double sinAbsorbedEnergyCutEllipseA 250 G4double sinAbsorbedEnergyCutEllipseAngle( 256 std::sin( abso 251 std::sin( absorbedEnergyCutEllipseAngle ) ); 257 252 258 if ( absorbedEnergyCutCLWidth <= 0. || 253 if ( absorbedEnergyCutCLWidth <= 0. || absorbedEnergyCutCRWidth <= 0. ) 259 { 254 { 260 hasAbsorbedEnergyCutTriggered = fa 255 hasAbsorbedEnergyCutTriggered = false; 261 } 256 } 262 else 257 else 263 { 258 { 264 G4double absorbedEnergyCutCLWidth 259 G4double absorbedEnergyCutCLWidth2( 265 absorbedEnergyCutCLWid 260 absorbedEnergyCutCLWidth * absorbedEnergyCutCLWidth ); 266 G4double absorbedEnergyCutCRWidth 261 G4double absorbedEnergyCutCRWidth2( 267 absorbedEnergyCutCRWid 262 absorbedEnergyCutCRWidth * absorbedEnergyCutCRWidth ); 268 263 269 hasAbsorbedEnergyCutTriggered = 264 hasAbsorbedEnergyCutTriggered = 270 std::pow( ( calorimeterEDLeft 265 std::pow( ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) * 271 cosAbsorbedEnerg 266 cosAbsorbedEnergyCutEllipseAngle + 272 ( calorimeterEDRight 267 ( calorimeterEDRight - absorbedEnergyCutCRCenter ) * 273 sinAbsorbedEnerg 268 sinAbsorbedEnergyCutEllipseAngle, 2 ) / 274 absorbedEnergyCutCLWidth2 + 269 absorbedEnergyCutCLWidth2 + 275 std::pow( - ( calorimeterEDLef 270 std::pow( - ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) * 276 sinAbsorbedEnerg 271 sinAbsorbedEnergyCutEllipseAngle + 277 ( calorimeterEDRight 272 ( calorimeterEDRight - absorbedEnergyCutCRCenter ) * 278 cosAbsorbedEnerg 273 cosAbsorbedEnergyCutEllipseAngle, 2 ) / 279 absorbedEnergyCutCRWidth2 < 274 absorbedEnergyCutCRWidth2 < 280 1; 275 1; 281 } 276 } 282 } 277 } 283 278 284 hasBasicTrigger = true; 279 hasBasicTrigger = true; 285 } 280 } 286 281 287 282 288 G4bool CexmcChargeExchangeReconstructor::HasF 283 G4bool CexmcChargeExchangeReconstructor::HasFullTrigger( void ) const 289 { 284 { 290 if ( ! hasBasicTrigger ) 285 if ( ! hasBasicTrigger ) 291 return false; 286 return false; 292 if ( useMassCut && ! hasMassCutTriggered ) 287 if ( useMassCut && ! hasMassCutTriggered ) 293 return false; 288 return false; 294 if ( useAbsorbedEnergyCut && ! hasAbsorbed 289 if ( useAbsorbedEnergyCut && ! hasAbsorbedEnergyCutTriggered ) 295 return false; 290 return false; 296 291 297 return true; 292 return true; 298 } << 299 << 300 << 301 void CexmcChargeExchangeReconstructor::SetExp << 302 << 303 { << 304 expectedMomentumAmp = particleGun->GetOrig << 305 } 293 } 306 294 307 295