Skip to content
Snippets Groups Projects
Commit e1c9a7e5 authored by Kristian Zarebski's avatar Kristian Zarebski
Browse files

Removed Unneeded Randomizer and Swapped to Using Tuples

parent f03bc7af
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 96 additions and 164 deletions
LHCbG4Tests/G4GammaToDiLeptonConversionTest/G4GammaToDiLeptonConversionTest.cc
+ 8
− 2
View file @ e1c9a7e5
......@@ -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");
......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorAction.hh
+ 3
− 5
View file @ e1c9a7e5
......@@ -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......
......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/PrimaryGeneratorMessenger.hh
+ 2
− 5
View file @ e1c9a7e5
......@@ -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......
......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/include/SteppingAction.hh
+ 9
− 2
View file @ e1c9a7e5
......@@ -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......
......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/scripts/G4GammaToDiLeptonConversionTest.py
+ 8
− 3
View file @ e1c9a7e5
......@@ -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...")
......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorAction.cc
+ 1
− 14
View file @ e1c9a7e5
......@@ -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);
}
......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/PrimaryGeneratorMessenger.cc
+ 1
− 19
View file @ e1c9a7e5
......@@ -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......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/RunAction.cc
+ 12
− 8
View file @ e1c9a7e5
......@@ -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++) {
......
LHCbG4Tests/G4GammaToDiLeptonConversionTest/src/SteppingAction.cc
+ 52
− 106
View file @ e1c9a7e5
......@@ -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......
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment