Skip to content
Snippets Groups Projects
calorimeter_clustering.py 3.63 KiB
Newer Older
  • Learn to ignore specific revisions
  • from Gaudi.Configuration import *
    import os
      
    from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
    from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
    from GaudiKernel import SystemOfUnits as units
    
    detector_name = "topside"
    if "JUGGLER_DETECTOR" in os.environ :
      detector_name = str(os.environ["JUGGLER_DETECTOR"])
    
    # todo add checks
    input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
    output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
    n_events = str(os.environ["JUGGLER_N_EVENTS"])
    
    geo_service  = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)])
    podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file])
    
    from Configurables import PodioInput
    from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
    from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
    
    from Configurables import Jug__Digi__HadronicCalDigi as HadCalorimeterDigi
    from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
    from Configurables import Jug__Digi__CrystalEndcapsDigi as CrystalEndcapsDigi
    
    from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
    from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
    
    from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
    from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
    
    podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits"], OutputLevel=DEBUG)
    
    ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
    copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2",OutputLevel=DEBUG) 
    calcopier = CalCopier("CalCopier", inputCollection="CrystalEcalHits", outputCollection="CrystalEcalHits2",OutputLevel=DEBUG) 
    
    emcaldigi = CrystalEndcapsDigi("ecal_digi", inputHitCollection="CrystalEcalHits", outputHitCollection="RawDigiEcalHits")
    ecdigi = EMCalorimeterDigi("ec_barrel_digi", inputHitCollection="EcalBarrelHits", outputHitCollection="RawEcalBarrelHits")
    
    crystal_ec_reco = CrystalEndcapsReco("crystal_ec_reco", inputHitCollection="RawDigiEcalHits", outputHitCollection="RecoEcalHits",
                                   minModuleEdep=1.0*units.MeV,OutputLevel=DEBUG)
    ecal_reco = EMCalReconstruction("ecal_reco", inputHitCollection="RawEcalBarrelHits", outputHitCollection="RecEcalBarrelHits",
                                   minModuleEdep=5.0*units.MeV,OutputLevel=DEBUG)
    
    ec_barrel_cluster = IslandCluster("ec_barrel_cluster", 
            inputHitCollection="RecEcalBarrelHits", 
            outputClusterCollection="EcalBarrelClusters",
            splitHitCollection="splitEcalBarrelHitCollection",
            minClusterCenterEdep=1*units.MeV, 
            groupRange=2.0,
            OutputLevel=DEBUG)
    crystal_ec_cluster = IslandCluster("crystal_ec_cluster", 
            inputHitCollection="RecoEcalHits", 
            outputClusterCollection="EcalClusters",
            minClusterCenterEdep=30*units.MeV, 
            groupRange=2.0,
            OutputLevel=DEBUG)
    
    clusterreco = RecoCoG("cluster_reco", clusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalBox_z_length",OutputLevel=DEBUG)
    
    out = PodioOutput("out", filename=output_rec_file)
    
    out.outputCommands = ["keep *"]
    
    ApplicationMgr(
        TopAlg = [podioinput, copier, calcopier,
                  ecdigi, emcaldigi, 
                  crystal_ec_reco, ecal_reco, 
                  ec_barrel_cluster, crystal_ec_cluster, clusterreco, out],
        EvtSel = 'NONE',
        EvtMax   = n_events,
        ExtSvc = [podioevent],
        OutputLevel=DEBUG
     )