From e1c9a7e5dda80138a0d8740c24c92be0dc35c32f Mon Sep 17 00:00:00 2001 From: Kristian Zarebski <kristian.alexander.zarebski@cern.ch> Date: Thu, 22 Mar 2018 16:45:52 +0100 Subject: [PATCH] Removed Unneeded Randomizer and Swapped to Using Tuples --- .../G4GammaToDiLeptonConversionTest.cc | 10 +- .../include/PrimaryGeneratorAction.hh | 8 +- .../include/PrimaryGeneratorMessenger.hh | 7 +- .../include/SteppingAction.hh | 11 +- .../G4GammaToDiLeptonConversionTest.py | 11 +- .../src/PrimaryGeneratorAction.cc | 15 +- .../src/PrimaryGeneratorMessenger.cc | 20 +-- .../src/RunAction.cc | 20 ++- .../src/SteppingAction.cc | 158 ++++++------------ 9 files changed, 96 insertions(+), 164 deletions(-) diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc index f098c28..527d89e 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc @@ -62,7 +62,7 @@ G4int nEvts = 100; G4double thickness = 0.3; //in mm -G4String physList = "FTFP_BERT"; +G4String dimu_xsec_fac = "1000"; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -91,6 +91,12 @@ int main(int argc,char** argv) { thickness = G4double(atof(argv[i])); ++i; } + if(G4String(argv[i]) == "--scale") + { + ++i; + dimu_xsec_fac = G4String(argv[i]); + ++i; + } } @@ -130,7 +136,7 @@ int main(int argc,char** argv) { // Increase Gamma2MuMu XSec to 1000x - UI->ApplyCommand("/testem/phys/SetGammaToMuPairFac 1000"); + UI->ApplyCommand("/testem/phys/SetGammaToMuPairFac "+dimu_xsec_fac); //Start Simulation UI->ApplyCommand("/run/beamOn"); diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh index 45b3517..4fa9549 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh @@ -47,18 +47,16 @@ class PrimaryGeneratorMessenger; class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction { public: - PrimaryGeneratorAction(DetectorConstruction*, G4double, G4int); + PrimaryGeneratorAction(DetectorConstruction*, G4double, G4int); ~PrimaryGeneratorAction(); public: - void SetRndmBeam(G4double val) {fRndmBeam = val;} - void GeneratePrimaries(G4Event*) override; + virtual void GeneratePrimaries(G4Event*); private: G4ParticleGun* fParticleGun; DetectorConstruction* fDetector; - G4double fRndmBeam; - PrimaryGeneratorMessenger* fGunMessenger; + PrimaryGeneratorMessenger* fGunMessenger; }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh index 0e08a95..0743dd3 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh @@ -48,13 +48,10 @@ class PrimaryGeneratorMessenger: public G4UImessenger public: PrimaryGeneratorMessenger(PrimaryGeneratorAction*); ~PrimaryGeneratorMessenger(); - - void SetNewValue(G4UIcommand*, G4String) override; - + private: PrimaryGeneratorAction* fAction; - G4UIdirectory* fGunDir; - G4UIcmdWithADouble* fRndmCmd; + G4UIdirectory* fGunDir; }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh index bc835af..96913fd 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh @@ -36,6 +36,8 @@ #include "Geant4/G4UserSteppingAction.hh" #include "Geant4/globals.hh" +#include "Geant4/G4ParticleDefinition.hh" +#include "Geant4/G4ParticleTypes.hh" class RunAction; @@ -47,12 +49,17 @@ public: SteppingAction(RunAction*); ~SteppingAction(); - void UserSteppingAction(const G4Step*) override; - + virtual void UserSteppingAction(const G4Step*); + private: RunAction* fRunAction; G4double fMuonMass; G4double fEMass; + + G4MuonPlus* muplus_def = G4MuonPlus::MuonPlusDefinition(); + G4Positron* eplus_def = G4Positron::PositronDefinition(); + + int tuple_index = -1; }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py index 617f8fa..9f30f0c 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py @@ -14,6 +14,7 @@ parser = argparse.ArgumentParser("DILEPTONPROD") parser.add_argument("--nevts", default=10000, help="Number of Events") parser.add_argument("--energies", default="[1000]", help="Energies List") parser.add_argument("--thickness", default="[0.3]", help="List of thicknesses") +parser.add_argument("--scale", default="100000", help="Scale DiMuon Conversion Cross Section by factor") parser.add_argument("--debug", default=False, action="store_true", help="Debug Mode") args = parser.parse_args() @@ -41,10 +42,14 @@ except TypeError: for energy in energies: for t in thickness: - _logger.info("%s\n\nRunning Gamma->l+ l- with options:\n\n\tNEvents :\t%s\n\tEnergy :\t%s GeV\n\tThickness :\t%s mm\n", - '='*40, args.nevts, energy, t) + _logger.info("%s\n\nRunning Gamma->l+ l- Conversion in Aluminium with options:\n"+ + "\n\tNumber of Events :\t%s "+ + "\n\tBeam Energy :\t%s GeV"+ + "\n\tTarget Thickness :\t%s mm"+ + "\n\tDiMuon XSec Scale Factor :\t%s\n", + '='*40, args.nevts, energy, t, args.scale) _logger.info('='*40) - subprocess.check_call("G4GammaToDiLeptonConversionTest.exe --nevts {} --energy {} --thickness {}".format(args.nevts, energy, t), shell=True) + subprocess.check_call("G4GammaToDiLeptonConversionTest.exe --nevts {} --energy {} --thickness {} --scale {}".format(args.nevts, energy, t, args.scale), shell=True) _logger.info("Creating Single Combined Plots ROOT File...") diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc index 47af5be..9141a2f 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc @@ -41,7 +41,6 @@ #include "Geant4/G4ParticleTable.hh" #include "Geant4/G4ParticleDefinition.hh" #include "Geant4/G4SystemOfUnits.hh" -#include "Geant4/Randomize.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -50,7 +49,6 @@ PrimaryGeneratorAction::PrimaryGeneratorAction( :G4VUserPrimaryGeneratorAction(), fParticleGun(0), fDetector(DC), - fRndmBeam(0), fGunMessenger(0) { fParticleGun = new G4ParticleGun(nevts); @@ -78,18 +76,7 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { //this function is called at the begining of event // - G4double x0 = -0.5*(fDetector->GetBoxDepth()); - G4double y0 = 0.*cm, z0 = 0.*cm; - - //randomize the beam, if requested. - // - if (fRndmBeam > 0.) - { - G4double rbeam = 0.5*(fDetector->GetBoxDepth())*fRndmBeam; - y0 = (2*G4UniformRand()-1.)*rbeam; - z0 = (2*G4UniformRand()-1.)*rbeam; - } - fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0)); + fParticleGun->SetParticlePosition(G4ThreeVector(-0.5*(fDetector->GetBoxDepth()),0.*cm,0.*cm)); fParticleGun->GeneratePrimaryVertex(anEvent); } diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc index 1ef2718..e778f31 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc @@ -42,36 +42,18 @@ PrimaryGeneratorMessenger::PrimaryGeneratorMessenger( PrimaryGeneratorAction* Gun) :G4UImessenger(),fAction(Gun), - fGunDir(0), - fRndmCmd(0) + fGunDir(0) { fGunDir = new G4UIdirectory("/testem/gun/"); fGunDir->SetGuidance("gun control"); - fRndmCmd = new G4UIcmdWithADouble("/testem/gun/rndm",this); - fRndmCmd->SetGuidance("random lateral extension on the beam"); - fRndmCmd->SetGuidance("in fraction of 0.5*sizeYZ"); - fRndmCmd->SetParameterName("rBeam",false); - fRndmCmd->SetRange("rBeam>=0.&&rBeam<=1."); - fRndmCmd->AvailableForStates(G4State_Idle); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... PrimaryGeneratorMessenger::~PrimaryGeneratorMessenger() { - delete fRndmCmd; delete fGunDir; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... - -void PrimaryGeneratorMessenger::SetNewValue(G4UIcommand* command, - G4String newValue) -{ - if (command == fRndmCmd) - {fAction->SetRndmBeam(fRndmCmd->GetNewDoubleValue(newValue));} -} - -//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... - diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc index 532af81..a76c886 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc @@ -47,8 +47,6 @@ #include "Geant4/G4RunManager.hh" #include "Geant4/G4MuonMinus.hh" -#include "Geant4/Randomize.hh" - #include <sstream> //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -107,6 +105,12 @@ void RunAction::BeginOfRunAction(const G4Run* aRun) fAnalysis->CreateNtupleDColumn("DiMuon_ThetaPlus"); fAnalysis->CreateNtupleDColumn("DiMuon_ThetaMinus"); fAnalysis->CreateNtupleDColumn("DiMuon_InvMass"); + fAnalysis->CreateNtupleDColumn("DiMuon_PX_Plus"); + fAnalysis->CreateNtupleDColumn("DiMuon_PY_Plus"); + fAnalysis->CreateNtupleDColumn("DiMuon_PZ_Plus"); + fAnalysis->CreateNtupleDColumn("DiMuon_PX_Minus"); + fAnalysis->CreateNtupleDColumn("DiMuon_PY_Minus"); + fAnalysis->CreateNtupleDColumn("DiMuon_PZ_Minus"); fAnalysis->FinishNtuple(); fAnalysis->CreateNtuple("DiElectron", "DiElectron Production"); @@ -115,6 +119,12 @@ void RunAction::BeginOfRunAction(const G4Run* aRun) fAnalysis->CreateNtupleDColumn("DiElectron_ThetaPlus"); fAnalysis->CreateNtupleDColumn("DiElectron_ThetaMinus"); fAnalysis->CreateNtupleDColumn("DiElectron_InvMass"); + fAnalysis->CreateNtupleDColumn("DiElectron_PX_Plus"); + fAnalysis->CreateNtupleDColumn("DiElectron_PY_Plus"); + fAnalysis->CreateNtupleDColumn("DiElectron_PZ_Plus"); + fAnalysis->CreateNtupleDColumn("DiElectron_PX_Minus"); + fAnalysis->CreateNtupleDColumn("DiElectron_PY_Minus"); + fAnalysis->CreateNtupleDColumn("DiElectron_PZ_Minus"); fAnalysis->FinishNtuple(); G4cout << "\n----> NTuple file is opened in " << fileName << G4endl; @@ -138,12 +148,6 @@ void RunAction::CountProcesses(G4String procName) void RunAction::EndOfRunAction(const G4Run*) { - // show Rndm status : not needed - // - //CLHEP::HepRandom::showEngineStatus(); - - //total number of process calls - // G4double countTot = 0.; G4cout << "\n Number of process calls --->"; for (size_t i=0; i< fProcCounter->size();i++) { diff --git a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc index df2417a..3ad44d9 100644 --- a/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc +++ b/LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc @@ -35,7 +35,6 @@ #include "RunAction.hh" #include "Geant4/G4SteppingManager.hh" #include "Geant4/G4VProcess.hh" -#include "Geant4/G4ParticleTypes.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -67,117 +66,64 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) G4ThreeVector PGamma = PrePoint->GetMomentum(); G4double Eplus(0), Eminus(0), ETot(0); - G4ThreeVector Pplus , Pminus, PTot; + G4ThreeVector Pplus , Pminus, PTot; const G4TrackVector* secondary = fpSteppingManager->GetSecondary(); - if(processName == "GammaToMuPair") - { - for (size_t lp=0; lp<(*secondary).size(); lp++) { - - if ((*secondary)[lp]->GetDefinition()==G4MuonPlus::MuonPlusDefinition()) { - Eplus = (*secondary)[lp]->GetTotalEnergy(); - Pplus = (*secondary)[lp]->GetMomentum(); - } else { - Eminus = (*secondary)[lp]->GetTotalEnergy(); - Pminus = (*secondary)[lp]->GetMomentum(); - } - - ETot += (*secondary)[lp]->GetTotalEnergy(); - PTot += (*secondary)[lp]->GetMomentum(); - - } - - G4double xPlus = Eplus/EGamma, xMinus = Eminus/EGamma; - G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus); - G4double GammaPlus = Eplus/fMuonMass; - G4double GammaMinus= Eminus/fMuonMass; - - G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); - - if(0.0 == thetaPlus || 0.0 == thetaMinus) { - G4cout << "SteppingAction for GammaToMuPair: " - << "thetaPlus= " << thetaPlus << " thetaMinus= " << thetaMinus - << " gamPlus= " << GammaPlus << " gamMinus= " << GammaMinus - << " " << thetaPlus *GammaPlus - thetaMinus*GammaMinus << G4endl; - return; - } - analysisManager->FillNtupleDColumn(0, 0, xPlus); - analysisManager->FillNtupleDColumn(0, 1, xMinus); - analysisManager->FillNtupleDColumn(0, 2, thetaPlus); - analysisManager->FillNtupleDColumn(0, 3, thetaMinus); - analysisManager->FillNtupleDColumn(0, 4, sqrt(pow(ETot,2)-pow(PTot.getR(),2))); - analysisManager->AddNtupleRow(0); - - /*analysisManager->FillH1(1,1./(1.+std::pow(thetaPlus*GammaPlus,2))); - analysisManager->FillH1(2,std::log10(thetaPlus*GammaPlus)); - - analysisManager->FillH1(3,std::log10(thetaMinus*GammaMinus)); - analysisManager->FillH1(4,std::log10(std::fabs(thetaPlus *GammaPlus - -thetaMinus*GammaMinus))); - - analysisManager->FillH1(5,xPlus); - analysisManager->FillH1(6,xMinus); - - analysisManager->FillH1(18, sqrt(pow(ETot,2)-pow(PTot.getR(),2))); - analysisManager->FillH1(19, thetaPlus); - analysisManager->FillH1(20, thetaMinus);*/ - } - else if(processName == "conv") - { - for (size_t lp=0; lp<(*secondary).size(); lp++) { - - if ((*secondary)[lp]->GetDefinition()==G4Positron::PositronDefinition()) { - Eplus = (*secondary)[lp]->GetTotalEnergy(); - Pplus = (*secondary)[lp]->GetMomentum(); - } else { - Eminus = (*secondary)[lp]->GetTotalEnergy(); - Pminus = (*secondary)[lp]->GetMomentum(); - } - - ETot += (*secondary)[lp]->GetTotalEnergy(); - PTot += (*secondary)[lp]->GetMomentum(); - - } - - G4double xPlus = Eplus/EGamma, xMinus = Eminus/EGamma; - G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus); - G4double GammaPlus = Eplus/fEMass; - G4double GammaMinus= Eminus/fEMass; - - G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); - - if(0.0 == thetaPlus || 0.0 == thetaMinus) { - G4cout << "SteppingAction for GammaToEPair: " - << "thetaPlus= " << thetaPlus << " thetaMinus= " << thetaMinus - << " gamPlus= " << GammaPlus << " gamMinus= " << GammaMinus - << " " << thetaPlus *GammaPlus - thetaMinus*GammaMinus << G4endl; - return; - } - analysisManager->FillNtupleDColumn(1, 0, xPlus); - analysisManager->FillNtupleDColumn(1, 1, xMinus); - analysisManager->FillNtupleDColumn(1, 2, thetaPlus); - analysisManager->FillNtupleDColumn(1, 3, thetaMinus); - analysisManager->FillNtupleDColumn(1, 4, sqrt(pow(ETot,2)-pow(PTot.getR(),2))); - analysisManager->AddNtupleRow(1); - - /*analysisManager->FillH1(21,1./(1.+std::pow(thetaPlus*GammaPlus,2))); - analysisManager->FillH1(22,std::log10(thetaPlus*GammaPlus)); - - analysisManager->FillH1(23,std::log10(thetaMinus*GammaMinus)); - analysisManager->FillH1(24,std::log10(std::fabs(thetaPlus *GammaPlus - -thetaMinus*GammaMinus))); - - analysisManager->FillH1(25,xPlus); - analysisManager->FillH1(26,xMinus); + + G4double particle_mass = -1; + tuple_index = -1; - analysisManager->FillH1(27, sqrt(pow(ETot,2)-pow(PTot.getR(),2))); - analysisManager->FillH1(28, thetaPlus); - analysisManager->FillH1(29, thetaMinus);*/ + if(processName == "conv") + { + particle_mass = fEMass; + tuple_index = 1; } else - { - G4cout << "Unrecognised Process!" << G4endl; //Should not get to here! + { + particle_mass = fMuonMass; + tuple_index = 0; + } + for (size_t lp=0; lp<(*secondary).size(); lp++) { + if ((*secondary)[lp]->GetDefinition()== muplus_def || (*secondary)[lp]->GetDefinition()== eplus_def) { + Eplus = (*secondary)[lp]->GetTotalEnergy(); + Pplus = (*secondary)[lp]->GetMomentum(); + } else { + Eminus = (*secondary)[lp]->GetTotalEnergy(); + Pminus = (*secondary)[lp]->GetMomentum(); + } + + ETot += (*secondary)[lp]->GetTotalEnergy(); + PTot += (*secondary)[lp]->GetMomentum(); + } + + G4double xPlus = Eplus/EGamma, xMinus = Eminus/EGamma; + G4double thetaPlus = PGamma.angle(Pplus), thetaMinus = PGamma.angle(Pminus); + G4double GammaPlus = Eplus/particle_mass; + G4double GammaMinus= Eminus/particle_mass; + + G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); + + if(0.0 == thetaPlus || 0.0 == thetaMinus) { + G4cout << "SteppingAction for GammaToMuPair: " + << "thetaPlus= " << thetaPlus << " thetaMinus= " << thetaMinus + << " gamPlus= " << GammaPlus << " gamMinus= " << GammaMinus + << " " << thetaPlus *GammaPlus - thetaMinus*GammaMinus << G4endl; + return; + } + analysisManager->FillNtupleDColumn(tuple_index, 0, xPlus); + analysisManager->FillNtupleDColumn(tuple_index, 1, xMinus); + analysisManager->FillNtupleDColumn(tuple_index, 2, thetaPlus); + analysisManager->FillNtupleDColumn(tuple_index, 3, thetaMinus); + analysisManager->FillNtupleDColumn(tuple_index, 4, sqrt(pow(ETot,2)-pow(PTot.getR(),2))); + analysisManager->FillNtupleDColumn(tuple_index, 5, Pplus.x()); + analysisManager->FillNtupleDColumn(tuple_index, 6, Pplus.y()); + analysisManager->FillNtupleDColumn(tuple_index, 7, Pplus.z()); + analysisManager->FillNtupleDColumn(tuple_index, 8, Pminus.x()); + analysisManager->FillNtupleDColumn(tuple_index, 9, Pminus.y()); + analysisManager->FillNtupleDColumn(tuple_index, 10, Pminus.z()); + analysisManager->AddNtupleRow(tuple_index); + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... -- GitLab