ChargeExchangeMC/src/CexmcChargeExchangeReconstructor.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /*
0027  * ============================================================================
0028  *
0029  *       Filename:  CexmcChargeExchangeReconstructor.cc
0030  *
0031  *    Description:  charge exchange reconstructor
0032  *
0033  *        Version:  1.0
0034  *        Created:  02.12.2009 15:17:13
0035  *       Revision:  none
0036  *       Compiler:  gcc
0037  *
0038  *         Author:  Alexey Radkov (),
0039  *        Company:  PNPI
0040  *
0041  * ============================================================================
0042  */
0043
0044 #include <cmath>
0045 #include <G4ThreeVector.hh>
0046 #include <G4LorentzVector.hh>
0047 #include "CexmcChargeExchangeReconstructor.hh"
0048 #include "CexmcChargeExchangeReconstructorMessenger.hh"
0049 #include "CexmcEnergyDepositStore.hh"
0050 #include "CexmcPrimaryGeneratorAction.hh"
0051 #include "CexmcParticleGun.hh"
0052 #include "CexmcProductionModel.hh"
0053 #include "CexmcRunManager.hh"
0054 #include "CexmcException.hh"
0055
0056
0057 CexmcChargeExchangeReconstructor::CexmcChargeExchangeReconstructor(
0058                             const CexmcProductionModel *  productionModel ) :
0059     outputParticleMass( 0 ), nucleusOutputParticleMass( 0 ),
0060     useTableMass( false ), useMassCut( false ), massCutOPCenter( 0 ),
0061     massCutNOPCenter( 0 ), massCutOPWidth( 0 ), massCutNOPWidth( 0 ),
0062     massCutEllipseAngle( 0 ), useAbsorbedEnergyCut( false ),
0063     absorbedEnergyCutCLCenter( 0 ), absorbedEnergyCutCRCenter( 0 ),
0064     absorbedEnergyCutCLWidth( 0 ), absorbedEnergyCutCRWidth( 0 ),
0065     absorbedEnergyCutEllipseAngle( 0 ), expectedMomentumAmp( -1 ),
0066     edCollectionAlgorithm( CexmcCollectEDInAllCrystals ),
0067     hasMassCutTriggered( false ), hasAbsorbedEnergyCutTriggered( false ),
0068     beamParticleIsInitialized( false ), particleGun( NULL ), messenger( NULL )
0069 {
0070     if ( ! productionModel )
0071         throw CexmcException( CexmcWeirdException );
0072
0073     productionModelData.incidentParticle =
0074                                     productionModel->GetIncidentParticle();
0075
0076     CexmcRunManager *  runManager( static_cast< CexmcRunManager * >(
0077                                             G4RunManager::GetRunManager() ) );
0078     const CexmcPrimaryGeneratorAction *  primaryGeneratorAction(
0079                         static_cast< const CexmcPrimaryGeneratorAction * >(
0080                                 runManager->GetUserPrimaryGeneratorAction() ) );
0081     CexmcPrimaryGeneratorAction *  thePrimaryGeneratorAction(
0082                         const_cast< CexmcPrimaryGeneratorAction * >(
0083                                 primaryGeneratorAction ) );
0084     particleGun = thePrimaryGeneratorAction->GetParticleGun();
0085
0086     productionModelData.nucleusParticle =
0087                                     productionModel->GetNucleusParticle();
0088     productionModelData.outputParticle =
0089                                     productionModel->GetOutputParticle();
0090     productionModelData.nucleusOutputParticle =
0091                                     productionModel->GetNucleusOutputParticle();
0092
0093     messenger = new CexmcChargeExchangeReconstructorMessenger( this );
0094 }
0095
0096
0097 CexmcChargeExchangeReconstructor::~CexmcChargeExchangeReconstructor()
0098 {
0099     delete messenger;
0100 }
0101
0102
0103 void  CexmcChargeExchangeReconstructor::SetupBeamParticle( void )
0104 {
0105     if ( *productionModelData.incidentParticle !=
0106                                         *particleGun->GetParticleDefinition() )
0107         throw CexmcException( CexmcBeamAndIncidentParticlesMismatch );
0108
0109     beamParticleIsInitialized = true;
0110 }
0111
0112
0113 void  CexmcChargeExchangeReconstructor::Reconstruct(
0114                                     const CexmcEnergyDepositStore *  edStore )
0115 {
0116     if ( ! beamParticleIsInitialized )
0117     {
0118         if ( *productionModelData.incidentParticle !=
0119                                         *particleGun->GetParticleDefinition() )
0120             throw CexmcException( CexmcBeamAndIncidentParticlesMismatch );
0121
0122         beamParticleIsInitialized = true;
0123     }
0124
0125     if ( edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
0126         collectEDInAdjacentCrystals = true;
0127
0128     ReconstructEntryPoints( edStore );
0129     if ( hasBasicTrigger )
0130         ReconstructTargetPoint();
0131     if ( hasBasicTrigger )
0132         ReconstructAngle();
0133
0134     G4ThreeVector  epLeft( calorimeterEPLeftWorldPosition -
0135                            targetEPWorldPosition );
0136     G4ThreeVector  epRight( calorimeterEPRightWorldPosition -
0137                             targetEPWorldPosition );
0138
0139     G4double  cosTheAngle( std::cos( theAngle ) );
0140     G4double  calorimeterEDLeft( edStore->calorimeterEDLeft );
0141     G4double  calorimeterEDRight( edStore->calorimeterEDRight );
0142
0143     if ( edCollectionAlgorithm == CexmcCollectEDInAdjacentCrystals )
0144     {
0145         calorimeterEDLeft = calorimeterEDLeftAdjacent;
0146         calorimeterEDRight = calorimeterEDRightAdjacent;
0147     }
0148
0149     //G4double  cosOutputParticleLAB(
0150         //( calorimeterEDLeft * cosAngleLeft +
0151           //calorimeterEDRight * cosAngleRight ) /
0152           //std::sqrt( calorimeterEDLeft * calorimeterEDLeft +
0153                      //calorimeterEDRight * calorimeterEDRight +
0154                      //calorimeterEDLeft * calorimeterEDRight * cosTheAngle ) );
0155
0156     outputParticleMass = std::sqrt( 2 * calorimeterEDLeft *
0157                                     calorimeterEDRight * ( 1 - cosTheAngle ) );
0158
0159     G4ThreeVector    opdpLeftMomentum( epLeft );
0160     opdpLeftMomentum.setMag( calorimeterEDLeft );
0161     G4ThreeVector    opdpRightMomentum( epRight );
0162     opdpRightMomentum.setMag( calorimeterEDRight );
0163     G4ThreeVector    opMomentum( opdpLeftMomentum + opdpRightMomentum );
0164
0165     /* opMass will be used only in calculation of output particle's total
0166      * energy, in other places outputParticleMass should be used instead */
0167     G4double         opMass( useTableMass ?
0168                              productionModelData.outputParticle->GetPDGMass() :
0169                              outputParticleMass );
0170     /* the formula below is equivalent to
0171      * calorimeterEDLeft + calorimeterEDRight if opMass = outputParticleMass */
0172     G4double         opEnergy( std::sqrt( opMomentum.mag2() +
0173                                           opMass * opMass ) );
0174     productionModelData.outputParticleLAB = G4LorentzVector( opMomentum,
0175                                                              opEnergy );
0176
0177     G4ThreeVector  incidentParticleMomentum( particleGun->GetOrigDirection() );
0178     G4double       incidentParticleMomentumAmp( expectedMomentumAmp > 0 ?
0179                     expectedMomentumAmp : particleGun->GetOrigMomentumAmp() );
0180     incidentParticleMomentum *= incidentParticleMomentumAmp;
0181
0182     G4double       incidentParticlePDGMass(
0183                         productionModelData.incidentParticle->GetPDGMass() );
0184     G4double       incidentParticlePDGMass2( incidentParticlePDGMass *
0185                                              incidentParticlePDGMass );
0186     G4double       incidentParticleEnergy(
0187         std::sqrt( incidentParticleMomentumAmp * incidentParticleMomentumAmp +
0188                    incidentParticlePDGMass2 ) );
0189
0190     productionModelData.incidentParticleLAB = G4LorentzVector(
0191                         incidentParticleMomentum, incidentParticleEnergy );
0192     G4double       nucleusParticlePDGMass(
0193                         productionModelData.nucleusParticle->GetPDGMass() );
0194     productionModelData.nucleusParticleLAB = G4LorentzVector(
0195                         G4ThreeVector( 0, 0, 0 ), nucleusParticlePDGMass );
0196
0197     G4LorentzVector  lVecSum( productionModelData.incidentParticleLAB +
0198                         productionModelData.nucleusParticleLAB );
0199     G4ThreeVector    boostVec( lVecSum.boostVector() );
0200
0201     productionModelData.nucleusOutputParticleLAB =
0202             lVecSum - productionModelData.outputParticleLAB;
0203
0204     productionModelData.incidentParticleSCM =
0205             productionModelData.incidentParticleLAB;
0206     productionModelData.nucleusParticleSCM =
0207             productionModelData.nucleusParticleLAB;
0208     productionModelData.outputParticleSCM =
0209             productionModelData.outputParticleLAB;
0210     productionModelData.nucleusOutputParticleSCM =
0211             productionModelData.nucleusOutputParticleLAB;
0212
0213     productionModelData.incidentParticleSCM.boost( -boostVec );
0214     productionModelData.nucleusParticleSCM.boost( -boostVec );
0215     productionModelData.outputParticleSCM.boost( -boostVec );
0216     productionModelData.nucleusOutputParticleSCM.boost( -boostVec );
0217
0218     G4ThreeVector  nopMomentum( incidentParticleMomentum - opMomentum );
0219     G4double       nopEnergy( incidentParticleEnergy + nucleusParticlePDGMass -
0220                               opEnergy );
0221     nucleusOutputParticleMass = std::sqrt( nopEnergy * nopEnergy -
0222                                            nopMomentum.mag2() );
0223
0224     if ( useMassCut )
0225     {
0226         G4double  cosMassCutEllipseAngle( std::cos( massCutEllipseAngle ) );
0227         G4double  sinMassCutEllipseAngle( std::sin( massCutEllipseAngle ) );
0228
0229         if ( massCutOPWidth <= 0. || massCutNOPWidth <= 0. )
0230         {
0231             hasMassCutTriggered = false;
0232         }
0233         else
0234         {
0235             G4double  massCutOPWidth2( massCutOPWidth * massCutOPWidth );
0236             G4double  massCutNOPWidth2( massCutNOPWidth * massCutNOPWidth );
0237
0238             hasMassCutTriggered =
0239                 std::pow( ( outputParticleMass - massCutOPCenter ) *
0240                               cosMassCutEllipseAngle +
0241                           ( nucleusOutputParticleMass - massCutNOPCenter ) *
0242                               sinMassCutEllipseAngle, 2 ) / massCutOPWidth2 +
0243                 std::pow( - ( outputParticleMass - massCutOPCenter ) *
0244                               sinMassCutEllipseAngle +
0245                           ( nucleusOutputParticleMass - massCutNOPCenter ) *
0246                               cosMassCutEllipseAngle, 2 ) / massCutNOPWidth2 <
0247                 1;
0248         }
0249     }
0250
0251     if ( useAbsorbedEnergyCut )
0252     {
0253         G4double  cosAbsorbedEnergyCutEllipseAngle(
0254                                 std::cos( absorbedEnergyCutEllipseAngle ) );
0255         G4double  sinAbsorbedEnergyCutEllipseAngle(
0256                                 std::sin( absorbedEnergyCutEllipseAngle ) );
0257
0258         if ( absorbedEnergyCutCLWidth <= 0. || absorbedEnergyCutCRWidth <= 0. )
0259         {
0260             hasAbsorbedEnergyCutTriggered = false;
0261         }
0262         else
0263         {
0264             G4double  absorbedEnergyCutCLWidth2(
0265                         absorbedEnergyCutCLWidth * absorbedEnergyCutCLWidth );
0266             G4double  absorbedEnergyCutCRWidth2(
0267                         absorbedEnergyCutCRWidth * absorbedEnergyCutCRWidth );
0268
0269             hasAbsorbedEnergyCutTriggered =
0270                 std::pow( ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) *
0271                               cosAbsorbedEnergyCutEllipseAngle +
0272                           ( calorimeterEDRight - absorbedEnergyCutCRCenter ) *
0273                               sinAbsorbedEnergyCutEllipseAngle, 2 ) /
0274                 absorbedEnergyCutCLWidth2 +
0275                 std::pow( - ( calorimeterEDLeft - absorbedEnergyCutCLCenter ) *
0276                               sinAbsorbedEnergyCutEllipseAngle +
0277                           ( calorimeterEDRight - absorbedEnergyCutCRCenter ) *
0278                               cosAbsorbedEnergyCutEllipseAngle, 2 ) /
0279                 absorbedEnergyCutCRWidth2 <
0280                 1;
0281         }
0282     }
0283
0284     hasBasicTrigger = true;
0285 }
0286
0287
0288 G4bool  CexmcChargeExchangeReconstructor::HasFullTrigger( void ) const
0289 {
0290     if ( ! hasBasicTrigger )
0291         return false;
0292     if ( useMassCut && ! hasMassCutTriggered )
0293         return false;
0294     if ( useAbsorbedEnergyCut && ! hasAbsorbedEnergyCutTriggered )
0295         return false;
0296
0297     return true;
0298 }
0299
0300
0301 void  CexmcChargeExchangeReconstructor::SetExpectedMomentumAmpDiff(
0302                                                             G4double  value )
0303 {
0304     expectedMomentumAmp = particleGun->GetOrigMomentumAmp() + value;
0305 }
0306