Commit d3e4d319 authored by Zach Marshall's avatar Zach Marshall Committed by Graeme Stewart
Browse files

Whole lot of cleaning - no functional changes really (CosmicGenerator-00-01-01)

parent 80234f33
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifndef COSMICGENERATOR_COSMICEVENTPARSER_H
#define COSMICGENERATOR_COSMICEVENTPARSER_H
#include <iostream>
#include "CLHEP/Vector/ThreeVector.h"
class CosmicEventParser {
public:
CosmicEventParser() : m_eventNumber(0), m_pdgId(0) {};
const CLHEP::HepLorentzVector& Vertex(void){return m_vertex;}
const CLHEP::HepLorentzVector& Momentum(void){return m_momentum;}
int pdgID(void){return m_pdgId;}
private:
CLHEP::HepLorentzVector m_vertex;
CLHEP::HepLorentzVector m_momentum;
int m_eventNumber;
int m_pdgId;
std::istream& read(std::istream& is);
std::ostream& write(std::ostream& os) const;
friend std::istream& operator >> (std::istream& is,CosmicEventParser& ev);
friend std::ostream& operator << (std::ostream& os,const CosmicEventParser& ev);
};
std::istream& operator >> (std::istream& is,CosmicEventParser& ev) {
return ev.read(is);
}
std::ostream& operator << (std::ostream& os,const CosmicEventParser& ev) {
return ev.write(os);
}
inline
std::istream& CosmicEventParser::read(std::istream& is)
{
int dummy;
int id;
CLHEP::Hep3Vector vert,mom;
double v_x,v_y,v_z;
is >> m_eventNumber >> dummy >> id >> v_x >> v_y >> v_z >> mom;
//
// rotate over pi in x-z plane
//
vert.setX(-v_x);
vert.setY( v_y);
vert.setZ(-v_z);
//
// convert to MeV's and mm units
//
mom = 1000*mom;
vert = 10*vert;
m_vertex.setVect(vert);
m_vertex.setE(0.);
m_momentum.setVect(mom);
double energy = sqrt(pow(105.66,2)+mom.mag2());
m_momentum.setE(energy);
if(id == 5) m_pdgId = 13;
else m_pdgId = -13;
return is;
}
inline
std::ostream& CosmicEventParser::write(std::ostream& os) const
{
int dummy(1);
int id(5);
if(m_pdgId == -13) id = 6;
os << m_eventNumber << " " << dummy << " " << id << " "
<< m_vertex.x() << " " << m_vertex.y() << " " << m_vertex.z() << " "
<< m_momentum.x() << " " << m_momentum.y() << " " << m_momentum.z();
return os;
}
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
/**
@class CosmicGenerator
@brief Cosmic generator. The output will be stored in the transient event store so it can be passed to the simulation.
@author
W. Seligman: Initial Code 08-Nov-2002,
based on work by M. Shapiro and I. Hinchliffe
Modification for increasing efficiency of muon hitting the detector:
H. Ma. March 17, 2006
Property: ExzCut:
if true, the method exzCut(...) will be called to apply a
energy dependent position cut on the surface.
This rejects low energy muons at large distance.
Property: RMax
Used by exzCut to reject non-projective muons, which are
too far out on the surface
Modifications to accomodate Pixel EndCap C Cosmic Test needs
Marian Zdrazil June 7, 2006 mzdrazil@lbl.gov
Modifications to accomodate replacement of Pixel EndCap C by a Pixel EndCap A
Marian Zdrazil November 24, 2006 mzdrazil@lbl.gov
It is easier and actually more useful to leave the EndCap A
in the vertical position (the way it is positioned in the ATLAS detector)
instead of rotating it clockwise by 90deg which corresponds to the
placement during the Pixel EndCap A cosmic test in SR1 in November 2006.
This is why we will generate cosmic muons coming from the positive Z-axis
direction better than rotating the whole setup in PixelGeoModel.
Modifications July 3rd 2007, Rob McPherson
- Fix mu+/mu- bug (always present in Athena versions)
- Fix sign of Py (since tag CosmicGenerator-00-00-21, muons only upward-going)
Optimize selection of events passed to Geant4 for full simulation:
- cut on energy based on pathlength in rock
- reweighting of generated cosmic rays
- geometrical cut in plane of pixel detector
Juerg Beringer November 2007 JBeringer@lgl.gov
Robert Cahn November 2007 RNCahn@lbl.gov
*/
#ifndef GENERATORMODULESCOSMICGEN_H
#define GENERATORMODULESCOSMICGEN_H
#include "GeneratorModules/GenModule.h"
#include "AthenaKernel/IAtRndmGenSvc.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "HepMC/Polarization.h"
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
class StoreGateSvc;
class ActiveStoreSvc;
class CosmicGun;
class CosmicGenerator:public GenModule {
public:
CosmicGenerator(const std::string& name, ISvcLocator* pSvcLocator);
virtual ~CosmicGenerator();
virtual StatusCode genInitialize();
virtual StatusCode callGenerator();
virtual StatusCode genFinalize();
virtual StatusCode fillEvt(GenEvent* evt);
CLHEP::HepLorentzVector generateVertex(void);
CLHEP::HepLorentzVector generateVertexReweighted(void);
static IAtRndmGenSvc* p_AtRndmGenSvc;
private:
ActiveStoreSvc* m_activeStore;
// event counter, used for event ID
int m_events, m_rejected,m_accepted;
std::vector<int> m_pdgCode;
float m_emin, m_emax;
float m_ctcut;
float m_xlow, m_xhig, m_zlow, m_zhig, m_yval, m_IPx, m_IPy, m_IPz, m_radius, m_zpos;
float m_tmin, m_tmax;
float m_stopped_tminus, m_stopped_tplus;
bool m_cavOpt;
int m_srOneOpt;
bool m_srOnePixECOpt;
bool m_swapYZAxis;
bool m_muonECOpt;
int m_printEvent, m_printMod;
int m_selection ;
float m_thetamin, m_thetamax, m_phimin, m_phimax;
float m_GeV;
float m_mm;
bool m_readfile;
bool m_readTTR;
std::string m_infile;
std::ifstream m_ffile;
// Event scalars, three-vectors, and four-vectors:
std::vector<CLHEP::HepLorentzVector> m_fourPos;
std::vector<CLHEP::HepLorentzVector> m_fourMom;
CLHEP::Hep3Vector m_center;
std::vector<HepMC::Polarization> m_polarization;
// Pointer to Athena MessageSvc.
IMessageSvc* p_msgSvc;
// Energy dependent position cut for muons to reach the detector.
bool exzCut(const CLHEP::Hep3Vector& pos,const CLHEP::HepLorentzVector& p);
// property for calling exzCut
bool m_exzCut ;
// maximum r used in exzCut
float m_rmax ;
// Calculation of pathlength in rock and whether cosmic ray is aimed towards
// the pixel detector
double pathLengthInRock(double xgen, double ygen, double zgen, double theta, double phi);
bool pointsAtPixels(double xgen, double ygen, double zgen, double theta, double phi);
// New optimization options
bool m_doPathlengthCut;
bool m_doAimedAtPixelsCut;
bool m_doReweighting;
double m_energyCutThreshold;
double m_ysurface;
double m_rvertmax;
double m_pixelplanemaxx;
double m_pixelplanemaxz;
double m_smearTR;
double m_smearTRp;
std::string m_recordName;
bool m_stopParticles;
};
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#ifndef COSMICGENERATOR_COSMICGUN_H
#define COSMICGENERATOR_COSMICGUN_H
#include "CLHEP/Vector/LorentzVector.h"
class CosmicGun{
public:
static CosmicGun* GetCosmicGun(void);
CLHEP::HepLorentzVector GenerateEvent(void);
CLHEP::HepLorentzVector GenerateVertex(void);
void SetEnergyRange(float emin, float emax);
void SetCosCut(float ctcut);
void PrintLevel(int printevt,int printmod);
int GetMuonCharge(void);
float InitializeGenerator(); // returns flux after all cuts in cm2/s
private:
CosmicGun(void);
static CosmicGun* mpointer;
int m_event;
int m_printevt, m_printmod;
float m_emin, m_emax;
float m_coscut;
};
#endif
package CosmicGenerator
author SBentvelsen
author rmcphers
use AtlasPolicy AtlasPolicy-*
use AtlasCLHEP AtlasCLHEP-* External
use AthenaKernel AthenaKernel-* Control
use GeneratorModules GeneratorModules-* Generators
use AtlasHepMC AtlasHepMC-* External
private
use GaudiInterface GaudiInterface-* External
use TrackRecord TrackRecord-* Simulation/G4Sim
end_private
public
# Build the library (and export the headers)
apply_pattern dual_use_library files="cosmic2.f CosmicGun.cxx CosmicGenerator.cxx"
apply_pattern declare_joboptions files="*.txt *.py"
apply_pattern declare_python_modules files="*.py"
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
/**
@mainpage CosmicGenerator
@author W. Seligman, M. Shapiro, I. Hinchliffe, M. Zdrazil,
@section General information
CosmicGenerator package is a generator used for the cosmic particle production.
The output will is stored in the transient event store so it can be passed to t
he simulation.
The CosmicGenerator is used e.g. by default by the G4 cosmic simulation as a so
urce of cosmic muons at ground level.
One of the features of the CosmicGenerator is the ability to filter primary muo
ns depending on their direction and energy.
If you look in jobOptions_ConfigCosmicProd.py, you will find that the following
properties can be set:
* CosmicGenerator.emin, CosmicGenerator.emax: energy range for the primary m
uon
* CosmicGenerator.xvert_low, CosmicGenerator.xvert_high, CosmicGenerator.zv
ert_low, CosmicGenerator.zvert_high: the (x,z) surface at ground level in which
the primary vertex has to be created
* CosmicGenerator.yvert_val: the y quota at which the primary vertexes must
be created (i.e. the "ground level")
* CosmicGenerator.ctcut: angular cut (wrt to the vertical)
Another set of properties allows further optimization:
* CosmicGenerator.OptimizeForCavern: if True, muons are passed to the simul
ation only if they are pointing towards the interaction point, within a given t
olerance. In order for this to work, the CosmicGenerator must be informed on wh
ere the IP is exactly. This is what the next properties are for
* CosmicGenerator.IPx, CosmicGenerator.IPy, CosmicGenerator.IPz: the (x,y,z
) coordinates of the IP
* CosmicGenerator.Radius: the tolerance of the direction filtering. Only mu
ons pointing inside a sphere centered in the IP with the given radius will be a
ccepted.
@htmlinclude used_packages.html
@include requirements
*/
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
############## Input: GenericGenerator ###############
def getInput_GenericCosmicGenerator(name="GenericCosmicGenerator", **kwargs):
## Configuring the Athena application for a 'generator' job
from G4AtlasApps.SimFlags import simFlags
simFlags.load_cosmics_flags()
## Set up standard algorithms and random seeds
from AthenaCommon.AppMgr import ServiceMgr
from PartPropSvc.PartPropSvcConf import PartPropSvc
ServiceMgr += PartPropSvc()
if not simFlags.RandomSeedList.checkForExistingSeed( "COSMICS"):
simFlags.RandomSeedList.addSeed( "COSMICS", 2040160768, 443921183 )
kwargs.setdefault('AtRndmGenSvc', simFlags.RandomSvc.get_Value())
##--------------------------------------------------------------
## CosmicGenerator parameters
##--------------------------------------------------------------
##
## Note that in this coordinate frame the y-axis points upward
## such that the cosmics arrive from upward to downward in y.
## The production vertex of cosmics is randomly distributed (flat)
## in the x-z plane with boundaries given below.
## The energy range is given as well.
##
## The following settings are tuned to scintillators of dimensions
## 140 x 0.5 x 100 cm^3 placed at +-115.0 cm
##
## NOTE: From G4AtlasApps-04-00-00 onwards, IDET-* cosmics
## commissioning layouts are not supported!
xvert_low = -301700.
xvert_hig = 298300.
zvert_low = -300000.
zvert_hig = 300000.
if simFlags.CosmicPtSlice.statusOn:
print "Configuring cosmic pT slice: %s" % simFlags.CosmicPtSlice.get_Value()
if simFlags.CosmicPtSlice == 'slice1':
kwargs.setdefault('emin', 10.*1000.) # 10 GeV
kwargs.setdefault('emax', 100.*1000.) # 100 GeV
xvert_low = -1000.*200. # 200 m
xvert_hig = 1000.*200. # 200 m
zvert_low = -1000.*200. # 200 m
zvert_hig = 1000.*200. # 200 m
elif simFlags.CosmicPtSlice == 'slice2':
kwargs.setdefault('emin', 100.*1000.) # 100 GeV
kwargs.setdefault('emax', 300.*1000.) # 300 GeV
xvert_low = -1000.*600. # 600 m
xvert_hig = 1000.*600. # 600 m
zvert_low = -1000.*600. # 600 m
zvert_hig = 1000.*600. # 600 m
elif simFlags.CosmicPtSlice == 'slice3':
kwargs.setdefault('emin', 300.*1000.) # 300 GeV
kwargs.setdefault('emax', 1000.*1000.) # 1000 GeV
xvert_low = -1000.*1000. # 1 km
xvert_hig = 1000.*1000. # 1 km
zvert_low = -1000.*1000. # 1 km
zvert_hig = 1000.*1000. # 1 km
elif simFlags.CosmicPtSlice == 'slice4':
kwargs.setdefault('emin', 1000.*1000.) # 1 TeV
kwargs.setdefault('emax', 5000.*1000.) # 5 TeV
xvert_low = -1000.*3000. # 3 km
xvert_hig = 1000.*3000. # 3 km
zvert_low = -1000.*3000. # 3 km
zvert_hig = 1000.*3000. # 3 km
elif simFlags.CosmicPtSlice != 'NONE':
print 'Slice name incorrect!'
# TODO: theApp.exit(1)?
import sys
sys.exit(1)
kwargs.setdefault('xvert_low', xvert_low)
kwargs.setdefault('xvert_hig', xvert_hig)
kwargs.setdefault('zvert_low', zvert_low)
kwargs.setdefault('zvert_hig', zvert_hig)
bedrockDX = (xvert_hig - xvert_low)/2.
bedrockDZ = (zvert_hig - zvert_low)/2.
if (bedrockDX > 350000 or bedrockDZ > 350000) :
newSize = max( bedrockDX , bedrockDZ )
print "Resizing bedrock (mm) to fit cosmic generator:",newSize
from AthenaCommon.Configurable import Configurable
if Configurable.allConfigurables.get('GeoModelSvc'):
GeoModelSvc = Configurable.allConfigurables.get('GeoModelSvc')
else:
GeoModelSvc = theApp.service('GeoModelSvc')
if (newSize <= 500000) : GeoModelSvc.CavernInfraVersionOverride = 'CavernInfra-03-Bedrock500'
elif (newSize <= 1000000) : GeoModelSvc.CavernInfraVersionOverride = 'CavernInfra-03-Bedrock1000'
elif (newSize <= 1500000) : GeoModelSvc.CavernInfraVersionOverride = 'CavernInfra-03-Bedrock1500'
elif (newSize <= 2000000) : GeoModelSvc.CavernInfraVersionOverride = 'CavernInfra-03-Bedrock2000'
elif (newSize <= 3000000) : GeoModelSvc.CavernInfraVersionOverride = 'CavernInfra-03-Bedrock3000'
elif (newSize <= 4000000) : GeoModelSvc.CavernInfraVersionOverride = 'CavernInfra-03-Bedrock4000'
elif (newSize <= 5000000) : GeoModelSvc.CavernInfraVersionOverride = 'CavernInfra-03-Bedrock5000'
else :
print "No need to resize the bedrock for cosmic generation"
## Adding the CosmicGenerator to the list of algs
from CosmicGenerator.CosmicGeneratorConf import CosmicGenerator
algorithm = CosmicGenerator(**kwargs)
from AthenaCommon.AlgSequence import AlgSequence
topSequence = AlgSequence()
if not hasattr(topSequence, 'CosmicGenerator'):
topSequence += algorithm
return algorithm
############## Input: Creating cosmics from scratch ###############
def getInput_EvgenCosmicGenerator(name="EvgenCosmicGenerator", **kwargs):
## Configuring the Athena application for a 'generator' job
import AthenaCommon.AtlasUnixGeneratorJob
from AthenaCommon.AthenaCommonFlags import athenaCommonFlags
athenaCommonFlags.PoolEvgenInput.set_Off()
from G4AtlasApps.SimFlags import simFlags
if not simFlags.CosmicFilterVolumeName.statusOn:
print "Warning CosmicFilterVolumeName not set using default (CaloEntryLayer)"
simFlags.CosmicFilterVolumeName = "CaloEntryLayer"
kwargs.setdefault('emin', 10000) # default =10000 #10 GeV
kwargs.setdefault('emax', 5000*1000) # 2 TeV - FIXME?!
if simFlags.CosmicFilterVolumeName == "Muon":
print 'Using muon Volume setup of Cosmic Generator...'
kwargs.setdefault('xvert_low', -301700.)
kwargs.setdefault('xvert_hig', 298300.)
kwargs.setdefault('zvert_low', -300000.)
kwargs.setdefault('zvert_hig', 300000.)
kwargs.setdefault('Radius', 20000.)
else:
print 'Using Non-muon Volume setup of Cosmic Generator...'
kwargs.setdefault('xvert_low', -200000.)
kwargs.setdefault('xvert_hig', 200000.)
kwargs.setdefault('zvert_low', -200000.)
kwargs.setdefault('zvert_hig', 200000.)
kwargs.setdefault('Radius', 10000.) #barrel length ~22m
kwargs.setdefault('yvert_val', 57300+41000.)
kwargs.setdefault('ctcut', 0.)
kwargs.setdefault('OptimizeForCavern', True)
kwargs.setdefault('IPx', 0.)
kwargs.setdefault('IPy', 0.)
kwargs.setdefault('IPz', 0.)
kwargs.setdefault('ReadTR', simFlags.ReadTR.statusOn)
#special settings from Juerg Beringer
if simFlags.CosmicFilterVolumeName == "Pixel" or simFlags.CosmicFilterVolumeName2 == "Pixel":
kwargs.setdefault('Radius', 2000.)
kwargs.setdefault('doPathLengthCut', True) # Optimization based on box cut in pixel detector plane
kwargs.setdefault('energyCutThreshold', 100.) # - margin of error for energy loss calculation (in MeV)
kwargs.setdefault('doAimedAtPixelsCut', True) # Optimization based on box cut in pixel detector plane
kwargs.setdefault('pixelplane_maxx', 1150.) # - require |x| < value in mm
kwargs.setdefault('pixelplane_maxz', 1650.) # - require |y| < value in mm
kwargs.setdefault('doReweighting', True) # Whether to use reweighting for cosmic ray generation
kwargs.setdefault('rvert_max', 300000.) # - radius in mm for generating primary vertex
#fix for bug: 49362
import sys
from AthenaCommon.AppMgr import ServiceMgr
ServiceMgr.EventSelector.EventsPerRun = int(2**31 - 1) #sys.maxint on a 32-bit machine
return getInput_GenericCosmicGenerator(name, **kwargs)
############## Input: Reading Cosmics from TrackRecord Input File ###############
def getInput_TrackRecordCosmicGenerator(name="TrackRecordCosmicGenerator", **kwargs):
## Configuring the Athena application for a 'track record' job
import G4AtlasApps.AtlasCosmicTrackRecordJob
from AthenaCommon.BeamFlags import jobproperties
if jobproperties.Beam.beamType() == "cosmics":
kwargs.setdefault('emin', 10000) # 10 GeV
kwargs.setdefault('emax', 2000*1000) # 2 TeV
kwargs.setdefault('xvert_low', -301700.)
kwargs.setdefault('xvert_hig', 298300.)
kwargs.setdefault('zvert_low', -300000.)
kwargs.setdefault('zvert_hig', 300000.)
kwargs.setdefault('Radius', 2000.0)
kwargs.setdefault('yvert_val', (57300. + 41000.))
kwargs.setdefault('ctcut', 0.0)
kwargs.setdefault('OptimizeForCavern', True)
kwargs.setdefault('IPx', 0.0)
kwargs.setdefault('IPy', 0.0)
kwargs.setdefault('IPz', 0.0)
from G4AtlasApps.SimFlags import simFlags
kwargs.setdefault('ReadTR', simFlags.ReadTR.statusOn)
kwargs.setdefault('TRSmearing', -1) #in millimeters, e.g. 10
kwargs.setdefault('TRPSmearing', -1) #in radians, e.g. 0.01
#kwargs.setdefault('OutputLevel', DEBUG) # for turning up output during testing
return getInput_GenericCosmicGenerator(name, **kwargs)
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
from AthenaCommon.CfgGetter import addAlgorithm
addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_GenericCosmicGenerator", "GenericCosmicGenerator")
addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_EvgenCosmicGenerator", "EvgenCosmicGenerator")
addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_TrackRecordCosmicGenerator", "TrackRecordCosmicGenerator")
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
from AthenaCommon.JobProperties import JobProperty
class CosmicFilterVolumeName(JobProperty):
"""
Declare the volume used to do cosmics filtering.
Values are G4 volume names, e.g. TRT_Barrel, TRT_EC, SCT_Barrel, Pixel, InnerDetector, Calo, CaloEntryLayer, Muon, ...
"""
statusOn = False
allowedTypes = [ 'str' ]
StoredValue = 'CaloEntryLayer'
class CosmicFilterVolumeName2(JobProperty):
"""
Declare a second volume used to do cosmics filtering.
Values are G4 volume names, e.g. TRT_Barrel, TRT_EC, SCT_Barrel, Pixel, InnerDetector, Calo, CaloEntryLayer, Muon, ...
"""
statusOn = False
allowedTypes = [ 'str' ]
StoredValue = 'NONE'