diff --git a/Generators/CosmicGenerator/CMakeLists.txt b/Generators/CosmicGenerator/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..16ffa7ba363eaea5590c83150e11a0498f6e7162 --- /dev/null +++ b/Generators/CosmicGenerator/CMakeLists.txt @@ -0,0 +1,38 @@ +################################################################################ +# Package: CosmicGenerator +################################################################################ + +# Declare the package name: +atlas_subdir( CosmicGenerator ) + +# Declare the package's dependencies: +atlas_depends_on_subdirs( PUBLIC + Control/AthenaKernel + Generators/GeneratorModules + PRIVATE + GaudiKernel ) + +# External dependencies: +find_package( CLHEP ) +find_package( HepMC ) + +# Component(s) in the package: +atlas_add_library( CosmicGeneratorLib + src/cosmic2.f + src/CosmicGun.cxx + src/CosmicGenerator.cxx + PUBLIC_HEADERS CosmicGenerator + INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} + DEFINITIONS ${CLHEP_DEFINITIONS} + LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaKernel GeneratorModulesLib + PRIVATE_LINK_LIBRARIES GaudiKernel ) + +atlas_add_component( CosmicGenerator + src/components/*.cxx + INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${HEPMC_INCLUDE_DIRS} + LINK_LIBRARIES ${CLHEP_LIBRARIES} ${HEPMC_LIBRARIES} AthenaKernel GeneratorModulesLib GaudiKernel CosmicGeneratorLib ) + +# Install files from the package: +atlas_install_python_modules( python/*.py ) +atlas_install_joboptions( share/*.txt share/*.py ) + diff --git a/Generators/CosmicGenerator/CosmicGenerator/CosmicGenerator.h b/Generators/CosmicGenerator/CosmicGenerator/CosmicGenerator.h index 8610e8837cc6b5ec81a33fc329b10297ab1baa97..cf701316818e796eac2876afe7ff620e83afd85e 100644 --- a/Generators/CosmicGenerator/CosmicGenerator/CosmicGenerator.h +++ b/Generators/CosmicGenerator/CosmicGenerator/CosmicGenerator.h @@ -70,7 +70,6 @@ in the vertical position (the way it is positioned in the ATLAS detector) #include <string> -class StoreGateSvc; class ActiveStoreSvc; class CosmicGun; @@ -115,7 +114,6 @@ private: float m_mm; bool m_readfile; - bool m_readTTR; std::string m_infile; std::ifstream m_ffile; @@ -125,9 +123,6 @@ private: 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); @@ -150,12 +145,6 @@ private: double m_rvertmax; double m_pixelplanemaxx; double m_pixelplanemaxz; - - double m_smearTR; - double m_smearTRp; - - std::string m_recordName; - bool m_stopParticles; }; #endif diff --git a/Generators/CosmicGenerator/cmt/requirements b/Generators/CosmicGenerator/cmt/requirements index 4730c7d592aa3b8fca45ca2721b6cf93b5241110..c5d4f763e0592ef4dc720d3a454b0b68f84420d7 100644 --- a/Generators/CosmicGenerator/cmt/requirements +++ b/Generators/CosmicGenerator/cmt/requirements @@ -12,7 +12,6 @@ use AtlasHepMC AtlasHepMC-* External private use GaudiInterface GaudiInterface-* External -use TrackRecord TrackRecord-* Simulation/G4Sim end_private public diff --git a/Generators/CosmicGenerator/python/CosmicGeneratorConfig.py b/Generators/CosmicGenerator/python/CosmicGeneratorConfig.py index 58a580783655dc66c25aa1bc5bdd14714ebdd808..797ceb94621a55cfe1de52a66fba2404ef689139 100644 --- a/Generators/CosmicGenerator/python/CosmicGeneratorConfig.py +++ b/Generators/CosmicGenerator/python/CosmicGeneratorConfig.py @@ -1,5 +1,151 @@ # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +class CavernPropertyCalculator(object): + + xvert_low = -1000.*200. # 200 m + xvert_high = 1000.*200. # 200 m + zvert_low = -1000.*200. # 200 m + zvert_high = 1000.*200. # 200 m + radius = 10000.#barrel length ~22m + emin = 10.*1000. # 10 GeV + emax = 100.*1000. # 100 GeV + EminDict = {'slice1': 10.*1000., # 10 GeV + 'slice2': 100.*1000., # 100 GeV + 'slice3': 300.*1000., # 300 GeV + 'slice4': 1000.*1000., # 1 TeV + 'NONE': 10.*1000. } + EmaxDict = {'slice1': 100.*1000., # 100 GeV + 'slice2': 300.*1000., # 300 GeV + 'slice3': 1000.*1000., # 1 TeV + 'slice4': 1000.*1000., # 5 TeV + 'NONE': 100.*1000. } + xvert_lowDict = {'slice1': -1000.*200., # 200 m + 'slice2': -1000.*600., # 600 m + 'slice3': -1000.*1000., # 1 km + 'slice4': -1000.*3000., # 3 km + 'NONE': -301700. } + xvert_highDict = {'slice1': 1000.*200., # 200 m + 'slice2': 1000.*600., # 600 m + 'slice3': 1000.*1000., # 1 km + 'slice4': 1000.*3000., # 3 km + 'NONE': 298300. } + zvert_lowDict = {'slice1': -1000.*200., # 200 m + 'slice2': -1000.*600., # 600 m + 'slice3': -1000.*1000., # 1 km + 'slice4': -1000.*3000., # 3 km + 'NONE': -300000. } + zvert_highDict = {'slice1': 1000.*200., # 200 m + 'slice2': 1000.*600., # 600 m + 'slice3': 1000.*1000., # 1 km + 'slice4': 1000.*3000., # 3 km + 'NONE': 300000. } + + def recalculate(self): + from G4AtlasApps.SimFlags import simFlags + if not simFlags.CosmicPtSlice.statusOn: + if simFlags.CosmicFilterVolumeName.statusOn and simFlags.CosmicFilterVolumeName() == 'Muon': + self.xvert_low = -301700. + self.xvert_high = 298300. + self.zvert_low = -1000.*300. + self.zvert_high = 1000.*300. + self.radius = 20000. + if (simFlags.CosmicFilterVolumeName.statusOn and simFlags.CosmicFilterVolumeName == "Pixel") or (simFlags.CosmicFilterVolumeName2.statusOn and simFlags.CosmicFilterVolumeName2 == "Pixel"): + self.radius = 2000. + else: + n_value = simFlags.CosmicPtSlice.get_Value() + self.xvert_low = self.xvert_lowDict[n_value] + self.xvert_high = self.zvert_lowDict[n_value] + self.zvert_low = self.xvert_highDict[n_value] + self.zvert_high = self.zvert_highDict[n_value] + #self.radius = + self.emin = self.EminDict[n_value] + self.emax = self.EmaxDict[n_value] + + def CosmicEmin(self): + """ + """ + self.recalculate() + return self.emin + + + def CosmicEmax(self): + """ + """ + self.recalculate() + return self.emax + + + def CosmicLowVertex_X(self): + """ + """ + self.recalculate() + return self.xvert_low + + + def CosmicHighVertex_X(self): + """ + """ + self.recalculate() + return self.xvert_high + + + def CosmicLowVertex_Z(self): + """ + """ + self.recalculate() + return self.zvert_low + + + def CosmicHighVertex_Z(self): + """ + """ + self.recalculate() + return self.zvert_high + + + def CosmicRadius(self): + """ + """ + self.recalculate() + return self.radius + + def BedrockDX(self): + """ + """ + self.recalculate() + return (self.xvert_high - self.xvert_low) + + def BedrockDZ(self): + """ + """ + self.recalculate() + return (self.zvert_high - self.zvert_low) + + def reconfigureCavernGeometry(self): + """ + Note that in this coordinate frame the y-axis points upward + such that the cosmics arrive from upward to downward in y. + """ + newSize = max( self.BedrockDX() , self.BedrockDZ() ) + if (newSize > 350000.): + 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: + from AthenaCommon.AppMgr import theApp + 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" + + ############## Input: GenericGenerator ############### def getInput_GenericCosmicGenerator(name="GenericCosmicGenerator", **kwargs): ## Configuring the Athena application for a 'generator' job @@ -14,6 +160,9 @@ def getInput_GenericCosmicGenerator(name="GenericCosmicGenerator", **kwargs): simFlags.RandomSeedList.addSeed( "COSMICS", 2040160768, 443921183 ) kwargs.setdefault('AtRndmGenSvc', simFlags.RandomSvc.get_Value()) + from CosmicGenerator.CosmicGeneratorConfig import CavernPropertyCalculator + theCavern = CavernPropertyCalculator() + ##-------------------------------------------------------------- ## CosmicGenerator parameters ##-------------------------------------------------------------- @@ -26,78 +175,42 @@ def getInput_GenericCosmicGenerator(name="GenericCosmicGenerator", **kwargs): ## ## 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. - + + kwargs.setdefault('emin', theCavern.CosmicEmin()) + kwargs.setdefault('emax', theCavern.CosmicEmax()) + #kwargs.setdefault('emin', 10000) # default =10000 #10 GeV + #kwargs.setdefault('emax', 5000*1000) # 2 TeV - FIXME?! + + kwargs.setdefault('xvert_low', theCavern.CosmicLowVertex_X() ) + kwargs.setdefault('xvert_hig', theCavern.CosmicHighVertex_X()) + kwargs.setdefault('zvert_low', theCavern.CosmicLowVertex_Z() ) + kwargs.setdefault('zvert_hig', theCavern.CosmicHighVertex_Z()) + kwargs.setdefault('Radius', theCavern.CosmicRadius()) + 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.) + + if simFlags.CosmicFilterVolumeName.statusOn: + print 'Using %s Volume setup of Cosmic Generator...' % simFlags.CosmicFilterVolumeName.get_Value() + #special settings from Juerg Beringer + if simFlags.CosmicFilterVolumeName == "Pixel" or simFlags.CosmicFilterVolumeName2 == "Pixel": + 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 + 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" + theCavern.reconfigureCavernGeometry() - ## Adding the CosmicGenerator to the list of algs - from CosmicGenerator.CosmicGeneratorConf import CosmicGenerator - algorithm = CosmicGenerator(**kwargs) + from AthenaCommon import CfgMgr + algorithm = CfgMgr.CosmicGenerator(**kwargs) from AthenaCommon.AlgSequence import AlgSequence topSequence = AlgSequence() if not hasattr(topSequence, 'CosmicGenerator'): @@ -117,45 +230,6 @@ def getInput_EvgenCosmicGenerator(name="EvgenCosmicGenerator", **kwargs): 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 @@ -164,31 +238,3 @@ def getInput_EvgenCosmicGenerator(name="EvgenCosmicGenerator", **kwargs): 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) - diff --git a/Generators/CosmicGenerator/python/CosmicGeneratorConfigDb.py b/Generators/CosmicGenerator/python/CosmicGeneratorConfigDb.py index 39676591547c6079e68a31ef9f50235e4f8526f7..4b8de87681f9aa2372d184cc225a858af43e55a1 100644 --- a/Generators/CosmicGenerator/python/CosmicGeneratorConfigDb.py +++ b/Generators/CosmicGenerator/python/CosmicGeneratorConfigDb.py @@ -2,7 +2,5 @@ from AthenaCommon.CfgGetter import addAlgorithm -addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_GenericCosmicGenerator", "GenericCosmicGenerator") -addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_EvgenCosmicGenerator", "EvgenCosmicGenerator") -addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_TrackRecordCosmicGenerator", "TrackRecordCosmicGenerator") - +addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_GenericCosmicGenerator", "GenericCosmicGenerator") +addAlgorithm("CosmicGenerator.CosmicGeneratorConfig.getInput_EvgenCosmicGenerator", "EvgenCosmicGenerator") diff --git a/Generators/CosmicGenerator/python/cosmics_flags.py b/Generators/CosmicGenerator/python/cosmics_flags.py index 6aa700ccee8992f2017a12c86d5cee9b69f1a32b..75b0e0d0aed357ae24e1751ea9f4e1ad4124dbf9 100644 --- a/Generators/CosmicGenerator/python/cosmics_flags.py +++ b/Generators/CosmicGenerator/python/cosmics_flags.py @@ -12,6 +12,7 @@ class CosmicFilterVolumeName(JobProperty): allowedTypes = [ 'str' ] StoredValue = 'CaloEntryLayer' + class CosmicFilterVolumeName2(JobProperty): """ Declare a second volume used to do cosmics filtering. @@ -22,14 +23,17 @@ class CosmicFilterVolumeName2(JobProperty): allowedTypes = [ 'str' ] StoredValue = 'NONE' + class CosmicPtSlice(JobProperty): """ Name of the pT slice used for cosmics simulation. """ statusOn = False allowedTypes = [ 'str' ] + allowedValues = [ 'slice1', 'slice2', 'slice3', 'slice4', 'NONE' ] StoredValue = 'NONE' + class CosmicFilterID(JobProperty): """ Particle ID to be filtered in cosmics processing. @@ -38,6 +42,7 @@ class CosmicFilterID(JobProperty): allowedTypes = [ 'str' ] StoredValue = '13' + class CosmicFilterPTmin(JobProperty): """ Particle minimum pT to be filtered in cosmics processing (in MeV). @@ -46,6 +51,7 @@ class CosmicFilterPTmin(JobProperty): allowedTypes = [ 'str' ] StoredValue = "5000" + class CosmicFilterPTmax(JobProperty): """ Particle maximum pT to be filtered in cosmics processing (in MeV). @@ -54,6 +60,7 @@ class CosmicFilterPTmax(JobProperty): allowedTypes = [ 'str' ] StoredValue = "6000" + class PixelEC_cosmic_SCsetup(JobProperty): """ Pixel EndCap C cosmic scintillator setup. @@ -64,3 +71,5 @@ class PixelEC_cosmic_SCsetup(JobProperty): allowedTypes = ['int'] allowedValues = [1,2] StoredValue = 1 + + diff --git a/Generators/CosmicGenerator/share/SetCosmicGenerator.py b/Generators/CosmicGenerator/share/SetCosmicGenerator.py index 4cd8c4b43fe7103c59eac00afe0e97e737862ddf..ea501fde8fa82862d589bb7d44a5b6428583c669 100644 --- a/Generators/CosmicGenerator/share/SetCosmicGenerator.py +++ b/Generators/CosmicGenerator/share/SetCosmicGenerator.py @@ -6,73 +6,24 @@ Set up cosmic generator for simulation + trigger. from G4AtlasApps.SimFlags import simFlags simFlags.load_cosmics_flags() assert hasattr(simFlags, "ReadTR") -if simFlags.ReadTR.statusOn: - import G4AtlasApps.AtlasCosmicTrackRecordJob -else: - import AthenaCommon.AtlasUnixGeneratorJob - + ## Set up standard algorithms and random seeds -from AthenaCommon.AppMgr import ServiceMgr -from PartPropSvc.PartPropSvcConf import PartPropSvc -ServiceMgr += PartPropSvc() from AthenaCommon.AlgSequence import AlgSequence topSequence = AlgSequence() -if hasattr(simFlags, 'RandomSeedList'): - ##Only for MC12 onwards - if not simFlags.RandomSeedList.checkForExistingSeed( "COSMICS"): - simFlags.RandomSeedList.addSeed( "COSMICS", 2040160768, 443921183 ) - -## Adding the CosmicGenerator to the list of algs -from CosmicGenerator.CosmicGeneratorConf import CosmicGenerator -topSequence += CosmicGenerator() +from AthenaCommon.CfgGetter import getAlgorithm +if simFlags.ReadTR.statusOn: + try: + cosmicGen = topSequence.TrackRecordCosmicGenerator + except: + cosmicGen = getAlgorithm("TrackRecordCosmicGenerator") +else: + try: + cosmicGen = topSequence.EvgenCosmicGenerator + except: + cosmicGen = getAlgorithm("EvgenCosmicGenerator") ## Adding the McEventCollection dumper to the list of algs #from TruthExamples.TruthExamplesConf import DumpMC #topSequence += DumpMC() -##-------------------------------------------------------------- -## 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! - -from AthenaCommon.BeamFlags import jobproperties -if jobproperties.Beam.beamType() == "cosmics": - if hasattr(simFlags, 'RandomSvc'): - ##Only for MC12 onwards - topSequence.CosmicGenerator.AtRndmGenSvc = simFlags.RandomSvc.get_Value() - topSequence.CosmicGenerator.emin = 10000; # 10 GeV - topSequence.CosmicGenerator.emax = 2000*1000; # 2 TeV - topSequence.CosmicGenerator.xvert_low = -301700. - topSequence.CosmicGenerator.xvert_hig = 298300. - topSequence.CosmicGenerator.zvert_low = -300000. - topSequence.CosmicGenerator.zvert_hig = 300000. - topSequence.CosmicGenerator.yvert_val = 57300 + 41000. - topSequence.CosmicGenerator.ctcut = 0.0 - topSequence.CosmicGenerator.OptimizeForCavern= True - topSequence.CosmicGenerator.IPx = 0.0 - topSequence.CosmicGenerator.IPy = 0.0 - topSequence.CosmicGenerator.IPz = 0.0 - topSequence.CosmicGenerator.Radius = 2000.0 - topSequence.CosmicGenerator.ReadTR = simFlags.ReadTR.statusOn - -if simFlags.ReadTR.statusOn: - topSequence.CosmicGenerator.TRSmearing = -1 #in millimeters, e.g. 10 - topSequence.CosmicGenerator.TRPSmearing = -1 #in radians, e.g. 0.01 - #topSequence.CosmicGenerator.OutputLevel = DEBUG # for turning up output during testing - -if simFlags.CosmicPtSlice.statusOn: - print "Configuring cosmic pT slice: %s" % simFlags.CosmicPtSlice - include("CosmicGenerator/CosmicSliceConfig.py") - -print topSequence.CosmicGenerator +print cosmicGen diff --git a/Generators/CosmicGenerator/src/CosmicGenerator.cxx b/Generators/CosmicGenerator/src/CosmicGenerator.cxx index d197c107d7de3e69d39cf3c3d1ebb3ac826721be..872674b047b3e3bde37f257fa95016afe84ce5a9 100644 --- a/Generators/CosmicGenerator/src/CosmicGenerator.cxx +++ b/Generators/CosmicGenerator/src/CosmicGenerator.cxx @@ -62,8 +62,6 @@ #include "CLHEP/Units/PhysicalConstants.h" #include "CLHEP/Random/RandFlat.h" -#include "TrackRecord/TrackRecordCollection.h" - #include <limits> #include <cmath> #include <vector> @@ -94,7 +92,6 @@ CosmicGenerator::CosmicGenerator(const std::string& name, m_mm = 10; m_readfile = false; m_activeStore = 0; - p_msgSvc = 0; m_events = 0; m_rejected = 0; @@ -112,9 +109,6 @@ CosmicGenerator::CosmicGenerator(const std::string& name, declareProperty("tmin", m_tmin =0. ); declareProperty("tmax", m_tmax =0. ); - declareProperty("stopped_tminus", m_stopped_tminus =-25. ); - declareProperty("stopped_tplus", m_stopped_tplus =25. ); - declareProperty("IPx", m_IPx =0. ); declareProperty("IPy", m_IPy =0. ); declareProperty("IPz", m_IPz =0. ); @@ -126,9 +120,6 @@ CosmicGenerator::CosmicGenerator(const std::string& name, declareProperty("SwapYZAxis", m_swapYZAxis = false); declareProperty("OptimizeForMuonEndCap", m_muonECOpt = false); declareProperty("ctcut", m_ctcut =0.35 ); - // declareProperty("ReadTTR", m_readTTR =false ); - declareProperty("ReadTR", m_readTTR =false ); - declareProperty("TRCollection" , m_recordName = "CosmicRecord" ); declareProperty("PrintEvent", m_printEvent=10); declareProperty("PrintMod", m_printMod=100); declareProperty("RMax", m_rmax = 10000000. ); @@ -148,11 +139,6 @@ CosmicGenerator::CosmicGenerator(const std::string& name, declareProperty("pixelplane_maxx",m_pixelplanemaxx = 1150); declareProperty("pixelplane_maxz",m_pixelplanemaxz = 1650); - // Smearing of the initial position for tracks - declareProperty("TRSmearing", m_smearTR=-1, "Smear the initial position of the track by up to this amount"); - declareProperty("TRPSmearing", m_smearTRp=-1, "Smear the momentum of the track by up to this amount"); - - declareProperty("StopParticles", m_stopParticles=false, "Stop the particles and make them decay within 25 ns"); } //-------------------------------------------------------------------------- @@ -337,135 +323,6 @@ StatusCode CosmicGenerator::callGenerator() { exit(1); return StatusCode::FAILURE; } - } else if(m_readTTR){ - - //const DataHandle <TimedTrackRecordCollection> coll; - const DataHandle <TrackRecordCollection> coll; - CHECK( evtStore()->retrieve(coll,m_recordName) ); - { - - ATH_MSG_INFO("retrieved "<<coll->size()<<" TTR hits; will smear position by "<< (m_smearTR>0?m_smearTR:0.) <<" mm and momentum by "<< (m_smearTRp>0?m_smearTRp:0.) <<" radians"); - - for (auto iterTTR : *coll) { - - const HepPDT::ParticleData* particle = particleData(abs(iterTTR.GetPDGCode())); - double mass = particle->mass().value(); - - double en = mass*mass+iterTTR.GetMomentum().x()*iterTTR.GetMomentum().x()+ - iterTTR.GetMomentum().y()*iterTTR.GetMomentum().y()+ - iterTTR.GetMomentum().z()*iterTTR.GetMomentum().z(); - - en=sqrt(en); - - ATH_MSG_VERBOSE("Reading back TTR:\n pos is "<<iterTTR.GetPosition() - <<"\n mom is "<<iterTTR.GetMomentum() - <<"\n pdg is "<<iterTTR.GetPDGCode() ); - - CLHEP::HepLorentzVector particle4Position( iterTTR.GetPosition().x(), - iterTTR.GetPosition().y(), - iterTTR.GetPosition().z(), - iterTTR.GetTime()); - - ATH_MSG_DEBUG( "Smearing position by up to " << m_smearTR << " mm and momentum by up to " << m_smearTRp << " radians" ); - if (m_smearTR>0){ - ATH_MSG_DEBUG( "Cosmic ray track was starting at " << particle4Position ); - - // if Z is maximal, move in X and Y; otherwise move in Z and "phi" - if ( particle4Position.z() == 22031 || particle4Position.z() == -22031 ){ //FIXME Hardcoded limits! - particle4Position.setX( particle4Position.x() + CLHEP::RandFlat::shoot(engine, -m_smearTR, m_smearTR) ); - particle4Position.setY( particle4Position.y() + CLHEP::RandFlat::shoot(engine, -m_smearTR, m_smearTR) ); - } else { - particle4Position.setZ( particle4Position.z() + CLHEP::RandFlat::shoot(engine, -m_smearTR, m_smearTR) ); - double R = sqrt( pow( particle4Position.x(),2 ) + pow(particle4Position.y(),2 ) ); - double dPhi = atan2( m_smearTR, R ); - dPhi = CLHEP::RandFlat::shoot( engine, -dPhi, dPhi ); - double theta = atan2( particle4Position.x() , particle4Position.y() ); - particle4Position.setX( R*sin( theta + dPhi ) ); - particle4Position.setY( R*cos( theta + dPhi ) ); - } - ATH_MSG_DEBUG( "Shifted cosmic ray to " << particle4Position ); - } - CLHEP::HepLorentzVector particle4Momentum( iterTTR.GetMomentum().x(), - iterTTR.GetMomentum().y(), - iterTTR.GetMomentum().z(), - en ); - if (m_smearTRp>0){ - - ATH_MSG_DEBUG( "Cosmic ray track momentum was " << particle4Momentum ); - - // Keep p - smear an angle, and then randomly spin that change in (0,2PI) - double dTheta = CLHEP::RandFlat::shoot(engine, 0, m_smearTRp); - double dPhi = CLHEP::RandFlat::shoot(engine, 0, 2*M_PI); - - // Need a perpendicular vector... - CLHEP::HepLorentzVector perpendicularMomentum( 1 , 0 , 0 , 0); - if ( particle4Momentum.x() != 0 ){ - if (particle4Momentum.y() == 0){ - perpendicularMomentum.setX( 0 ); - perpendicularMomentum.setY( 1 ); - } else { - perpendicularMomentum.setX( particle4Momentum.y() ); - perpendicularMomentum.setY( particle4Momentum.x() ); - } - } - - // Now scale it based on dTheta - double tempP = pow(particle4Momentum.x(),2)+pow(particle4Momentum.y(),2)+pow(particle4Momentum.z(),2); - if ( tempP==0 ) { - ATH_MSG_DEBUG("Our initial momentum had zero magnitude!!"); - perpendicularMomentum.setX(0); - perpendicularMomentum.setY(0); - } else if ( tan(dTheta) == 0 ){ - ATH_MSG_DEBUG("Randomly deciding to keep the vector's direction..."); - perpendicularMomentum.setX(0); - perpendicularMomentum.setY(0); - } else { - double scale = ( tempP ) * sin(dTheta) / ( pow(perpendicularMomentum.x(),2)+pow(perpendicularMomentum.y(),2) ); - perpendicularMomentum.setX( perpendicularMomentum.x() * scale ); - perpendicularMomentum.setY( perpendicularMomentum.y() * scale ); - - // Rotate perpendicularMomentum by dPhi around particle4Momentum - perpendicularMomentum.rotate( dPhi , particle4Momentum.vect() ); - } - particle4Momentum.setX( particle4Momentum.x() + perpendicularMomentum.x() ); - particle4Momentum.setY( particle4Momentum.y() + perpendicularMomentum.y() ); - particle4Momentum.setZ( particle4Momentum.z() + perpendicularMomentum.z() ); - - // Rescale (very small effect, but want to include it) - double scale2 = tempP==0? 1 : (pow(particle4Momentum.x(),2)+pow(particle4Momentum.y(),2)+pow(particle4Momentum.z(),2)) / tempP; - particle4Momentum.setX( particle4Momentum.x() * scale2 ); - particle4Momentum.setY( particle4Momentum.y() * scale2 ); - particle4Momentum.setZ( particle4Momentum.z() * scale2 ); - ATH_MSG_DEBUG( "Rotated the vector by " << perpendicularMomentum ); - ATH_MSG_DEBUG( "And resulting momentum is " << particle4Momentum ); - } - if (m_stopParticles){ - ATH_MSG_DEBUG( "Will stop the track record where it is and give it mass " << mass ); - particle4Momentum.setX(0); - particle4Momentum.setY(0); - particle4Momentum.setZ(0); - particle4Momentum.setT(mass); - - double settime=CLHEP::RandFlat::shoot(engine,m_stopped_tminus, m_stopped_tplus); - ATH_MSG_DEBUG( "Setting particle time to something uniform between "<<m_stopped_tminus<<" and "<<m_stopped_tplus<<" ns : " << settime ); - particle4Position.setT(settime*CLHEP::c_light); // ct in mm - } - - m_fourPos.push_back( particle4Position ); - m_fourMom.push_back( particle4Momentum ); - - m_pdgCode.push_back(iterTTR.GetPDGCode()); - HepMC::Polarization thePolarization; - thePolarization.set_normal3d(HepGeom::Normal3D<double>(0,0,0)); - m_polarization.push_back(thePolarization); - - if (m_stopParticles){ - ATH_MSG_DEBUG( "Only one per event!!" ); - break; - } - } - - } } else {