diff --git a/FCCSW_ecal/neighbours_theta.py b/FCCSW_ecal/neighbours_theta.py index ca0094f..0e204c0 100644 --- a/FCCSW_ecal/neighbours_theta.py +++ b/FCCSW_ecal/neighbours_theta.py @@ -1,10 +1,17 @@ +from Configurables import GeoSvc from Configurables import ApplicationMgr from Configurables import CreateFCCeeCaloNeighbours import os -from Gaudi.Configuration import * +from Gaudi.Configuration import INFO, DEBUG + +# produce neighbour map also for HCal +doHCal = True + +# only relevant if doHCal is True, set to True +# for combined ECal+HCal topoclusters +linkECalHCalBarrels = True # Detector geometry -from Configurables import GeoSvc geoservice = GeoSvc("GeoSvc") # if K4GEO is empty, this should use relative path to working directory path_to_detector = os.environ.get("K4GEO", "") @@ -18,18 +25,40 @@ path_to_detector, _det) for _det in detectors_to_use] geoservice.OutputLevel = INFO -# Geant4 service -# Configures the Geant simulation: geometry, physics list and user actions -neighbours = CreateFCCeeCaloNeighbours("neighbours", - outputFileName = "neighbours_map_barrel_thetamodulemerged.root", - readoutNames=["ECalBarrelModuleThetaMerged"], - systemNames=["system"], - systemValues=[4], - activeFieldNames=["layer"], - activeVolumesNumbers=[12], - includeDiagonalCells=False, - connectBarrels=False, - OutputLevel=DEBUG) +if doHCal: + # create the neighbour file for ECAL+HCAL barrel cells + neighbours = CreateFCCeeCaloNeighbours("neighbours", + outputFileName="neighbours_map_ecalB_thetamodulemerged_hcalB_thetaphi.root", + readoutNames=[ + "ECalBarrelModuleThetaMerged", "BarHCal_Readout_phitheta"], + systemNames=["system", "system"], + systemValues=[4, 8], + activeFieldNames=["layer", "layer"], + activeVolumesNumbers=[12, 13], + activeVolumesTheta=[ + [], + [ + 0.788969, 0.797785, 0.806444, 0.814950, 0.823304, + 0.839573, 0.855273, 0.870425, 0.885051, 0.899172, + 0.912809, 0.938708, 0.962896 + ] + ], + includeDiagonalCells=False, + connectBarrels=linkECalHCalBarrels, + OutputLevel=DEBUG) +else: + # create the neighbour file for ECAL+HCAL barrel cells + neighbours = CreateFCCeeCaloNeighbours("neighbours", + outputFileName="neighbours_map_ecalB_thetamodulemerged.root", + readoutNames=[ + "ECalBarrelModuleThetaMerged"], + systemNames=["system"], + systemValues=[4], + activeFieldNames=["layer"], + activeVolumesNumbers=[12], + includeDiagonalCells=False, + connectBarrels=False, + OutputLevel=DEBUG) # ApplicationMgr ApplicationMgr(TopAlg=[], diff --git a/FCCSW_ecal/noise_map_theta.py b/FCCSW_ecal/noise_map_theta.py index fa9393c..76553cf 100644 --- a/FCCSW_ecal/noise_map_theta.py +++ b/FCCSW_ecal/noise_map_theta.py @@ -1,74 +1,117 @@ -from Gaudi.Configuration import * -# Detector geometry from Configurables import GeoSvc -geoservice = GeoSvc("GeoSvc") +from Configurables import ApplicationMgr +from Configurables import CreateFCCeeCaloNoiseLevelMap +# from Configurables import ReadNoiseFromFileTool +from Configurables import ReadNoiseVsThetaFromFileTool +from Configurables import ConstNoiseTool +from Configurables import CellPositionsECalBarrelModuleThetaSegTool import os +from Gaudi.Configuration import INFO, DEBUG + +doHCal = True + +# Detector geometry +geoservice = GeoSvc("GeoSvc") # if K4GEO is empty, this should use relative path to working directory path_to_detector = os.environ.get("K4GEO", "") print(path_to_detector) -detectors_to_use=[ +detectors_to_use = [ 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml' ] # prefix all xmls with path_to_detector -geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use] +geoservice.detectors = [os.path.join( + path_to_detector, _det) for _det in detectors_to_use] geoservice.OutputLevel = INFO +# readout names for ECAL and HCAL (latter is ignored if doHCal is False) ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" hcalBarrelReadoutName = "BarHCal_Readout_phitheta" -BarrelNoisePath = os.environ['FCCBASEDIR']+"/LAr_scripts/data/elecNoise_ecalBarrelFCCee_theta.root" +# names of file and histograms with noise per layer vs theta for barrel ECAL +BarrelNoisePath = os.environ['FCCBASEDIR'] + \ + "/LAr_scripts/data/elecNoise_ecalBarrelFCCee_theta.root" ecalBarrelNoiseHistName = "h_elecNoise_fcc_" -from Configurables import CellPositionsECalBarrelModuleThetaSegTool +# cell positioning and noise tool for the ecal barrel ECalBcells = CellPositionsECalBarrelModuleThetaSegTool("CellPositionsECalBarrel", - readoutName = ecalBarrelReadoutName) + readoutName=ecalBarrelReadoutName) -from Configurables import CreateFCCeeCaloNoiseLevelMap, ReadNoiseFromFileTool, ReadNoiseVsThetaFromFileTool ECalNoiseTool = ReadNoiseVsThetaFromFileTool("ReadNoiseFromFileToolECal", - useSegmentation = False, - cellPositionsTool = ECalBcells, - readoutName = ecalBarrelReadoutName, - noiseFileName = BarrelNoisePath, - elecNoiseHistoName = ecalBarrelNoiseHistName, - setNoiseOffset = False, - activeFieldName = "layer", - addPileup = False, - numRadialLayers = 12, - scaleFactor = 1/1000., #MeV to GeV - OutputLevel = INFO) + useSegmentation=False, + cellPositionsTool=ECalBcells, + readoutName=ecalBarrelReadoutName, + noiseFileName=BarrelNoisePath, + elecNoiseHistoName=ecalBarrelNoiseHistName, + setNoiseOffset=False, + activeFieldName="layer", + addPileup=False, + numRadialLayers=12, + scaleFactor=1 / 1000., # MeV to GeV + OutputLevel=INFO) -# HCAL noise file has yet to be created/implemented -# HCalNoiseTool = ReadNoiseFromFileTool("ReadNoiseFromFileToolHCal", -# readoutName = hcalBarrelReadoutName, -# noiseFileName = BarrelNoisePath, -# elecNoiseHistoName = ecalBarrelNoiseHistName, -# setNoiseOffset = False, -# activeFieldName = "layer", -# addPileup = False, -# numRadialLayers = 12, -# scaleFactor = 1/1000., #MeV to GeV -# OutputLevel = INFO) -HCalNoiseTool = None +if doHCal: + # noise tool for the HCAL barrel + # HCAL noise file has yet to be created/implemented + # HCalNoiseTool = ReadNoiseFromFileTool("ReadNoiseFromFileToolHCal", + # readoutName = hcalBarrelReadoutName, + # noiseFileName = BarrelNoisePath, + # elecNoiseHistoName = ecalBarrelNoiseHistName, + # setNoiseOffset = False, + # activeFieldName = "layer", + # addPileup = False, + # numRadialLayers = 12, + # scaleFactor = 1/1000., #MeV to GeV + # OutputLevel = INFO) + # ConstNoiseTool provides constant noise for all calo subsystems + # here we are going to use it only for hcal barrel + HCalNoiseTool = ConstNoiseTool("ConstNoiseTool") + + # create the noise file for ECAL+HCAL barrel cells + noisePerCell = CreateFCCeeCaloNoiseLevelMap("noisePerCell", + ECalBarrelNoiseTool=ECalNoiseTool, + ecalBarrelSysId=4, + HCalBarrelNoiseTool=HCalNoiseTool, + hcalBarrelSysId=8, + readoutNames=[ + ecalBarrelReadoutName, hcalBarrelReadoutName], + systemNames=[ + "system", "system"], + systemValues=[4, 8], + activeFieldNames=[ + "layer", "layer"], + activeVolumesNumbers=[12, 13], + # activeVolumesEta = [1.2524, 1.2234, 1.1956, 1.1561, 1.1189, 1.0839, 1.0509, 0.9999, 0.9534, 0.91072], + activeVolumesTheta=[ + [], + [ + 0.788969, 0.797785, 0.806444, 0.81495, 0.823304, + 0.839573, 0.855273, 0.870425, 0.885051, 0.899172, + 0.912809, 0.938708, 0.962896 + ] + ], + outputFileName="cellNoise_map_electronicsNoiseLevel_ecalB_thetamodulemerged_hcalB_thetaphi.root", + OutputLevel=DEBUG) +else: + # create the noise file for ECAL barrel cells + noisePerCell = CreateFCCeeCaloNoiseLevelMap("noisePerCell", + ECalBarrelNoiseTool=ECalNoiseTool, + ecalBarrelSysId=4, + HCalBarrelNoiseTool=None, + hcalBarrelSysId=8, + readoutNames=[ecalBarrelReadoutName], + systemNames=["system"], + systemValues=[4], + activeFieldNames=["layer"], + activeVolumesNumbers=[12], + activeVolumesTheta=[[]], + outputFileName="cellNoise_map_electronicsNoiseLevel_ecalB_thetamodulemerged.root", + OutputLevel=DEBUG) -# when HCal noise is ready, add also hcal readout info -noisePerCell = CreateFCCeeCaloNoiseLevelMap("noisePerCell", - ECalBarrelNoiseTool = ECalNoiseTool, - ecalBarrelSysId = 4, - HCalBarrelNoiseTool = HCalNoiseTool, - readoutNames=[ecalBarrelReadoutName], - systemNames=["system"], - systemValues=[4], - activeFieldNames=["layer"], - activeVolumesNumbers = [12], - #activeVolumesEta = [1.2524, 1.2234, 1.1956, 1.1561, 1.1189, 1.0839, 1.0509, 0.9999, 0.9534, 0.91072], - outputFileName = "cellNoise_map_electronicsNoiseLevel_thetamodulemerged.root", - OutputLevel = DEBUG) # ApplicationMgr -from Configurables import ApplicationMgr -ApplicationMgr( TopAlg = [], - EvtSel = 'NONE', - EvtMax = 1, - # order is important, as GeoSvc is needed by G4SimSvc - ExtSvc = [geoservice, noisePerCell], - OutputLevel=INFO +ApplicationMgr(TopAlg=[], + EvtSel='NONE', + EvtMax=1, + # order is important, as GeoSvc is needed by G4SimSvc + ExtSvc=[geoservice, noisePerCell], + OutputLevel=INFO ) diff --git a/FCCSW_ecal/printNeighbours.C b/FCCSW_ecal/printNeighbours.C index 9f0a417..8ba5f42 100644 --- a/FCCSW_ecal/printNeighbours.C +++ b/FCCSW_ecal/printNeighbours.C @@ -1,9 +1,23 @@ +#include "TTree.h" +#include "TFile.h" +#include "TRandom.h" +#include +using namespace std; + TTree *T = nullptr; -// const std::string filename = "neighbours_map_barrel_thetamodulemerged.root"; +//const std::string filename = "neighbours_map_barrel_thetamodulemerged.root"; const std::string filename = "neighbours_map_HCalBarrel.root"; const std::string treename = "neighbours"; + +TTree *Tnoise = nullptr; +const std::string filenameNoise = "cellNoise_map_electronicsNoiseLevel_ecalB_thetamodulemerged_hcalB_thetaphi.root"; +const std::string treenameNoise = "noisyCells"; + ULong64_t cID; std::vector *neighbours = 0; +ULong64_t cIDNoise; +double noiseLevel; +double noiseOffset; bool useHCalReadoutWithRows = false; @@ -79,6 +93,29 @@ ULong_t HCalBarrelPhiBin(ULong_t cellID) return 999999999; } +void LoadNeighboursMap() +{ + if (T == nullptr && filename!="") + { + TFile *f = TFile::Open(filename.c_str(), "READ"); + T = (TTree *)f->Get(treename.c_str()); + T->SetBranchAddress("cellId", &cID); + T->SetBranchAddress("neighbours", &neighbours); + } +} + +void LoadNoiseMap() +{ + if (Tnoise == nullptr && filenameNoise!="") + { + TFile *fNoise = TFile::Open(filenameNoise.c_str(), "READ"); + Tnoise = (TTree *)fNoise->Get(treenameNoise.c_str()); + Tnoise->SetBranchAddress("cellId", &cIDNoise); + Tnoise->SetBranchAddress("noiseLevel", &noiseLevel); + Tnoise->SetBranchAddress("noiseOffset", &noiseOffset); + } +} + void printCell(ULong_t cellID) { cout << "cellID: " << cellID << endl; @@ -97,28 +134,33 @@ void printCell(ULong_t cellID) cout << endl; } -void LoadNeighboursMap() -{ - if (T == nullptr) - { - TFile *f = TFile::Open(filename.c_str(), "READ"); - T = (TTree *)f->Get(treename.c_str()); - T->SetBranchAddress("cellId", &cID); - T->SetBranchAddress("neighbours", &neighbours); - } -} - -void printCellAndNeighbours(ULong64_t iEntry) +void printCell(ULong64_t iEntry, bool showNeighbours, bool showNoise) { T->GetEntry(iEntry); cout << "=================================================" << endl; cout << endl; printCell(cID); - cout << "Neighbours: " << endl - << endl; - for (unsigned int i = 0; i < neighbours->size(); i++) + if (showNeighbours) + { + cout << "Neighbours: " << endl + << endl; + for (unsigned int i = 0; i < neighbours->size(); i++) + { + printCell(neighbours->at(i)); + } + } + if (showNoise) { - printCell(neighbours->at(i)); + for (ULong64_t jEntry = 0; jEntry < Tnoise->GetEntries(); jEntry++) + { + Tnoise->GetEntry(jEntry); + if (cIDNoise == cID) + { + cout << "Noise: level = " << noiseLevel << " , offset = " << noiseOffset + << endl << endl; + break; + } + } } cout << "=================================================" << endl; } @@ -131,7 +173,7 @@ void printNeighboursOfCell(ULong_t cellID) T->GetEntry(iEntry); if (cID == cellID) { - printCellAndNeighbours(iEntry); + printCell(iEntry, true, false); return; } } @@ -144,6 +186,6 @@ void printNeighbours(int n = 10) for (int i = 0; i < n; i++) { int entry = (int)gRandom->Uniform(T->GetEntries()); - printCellAndNeighbours(entry); + printCell(entry, true, false); } } diff --git a/FCCSW_ecal/run_thetamodulemerged.py b/FCCSW_ecal/run_thetamodulemerged.py index 4e3f7a3..98be5eb 100644 --- a/FCCSW_ecal/run_thetamodulemerged.py +++ b/FCCSW_ecal/run_thetamodulemerged.py @@ -27,11 +27,16 @@ from Configurables import TopoCaloNeighbours from Configurables import TopoCaloNoisyCells from Configurables import CaloTopoClusterFCCee +from Configurables import RewriteBitfield +from Gaudi.Configuration import INFO +# , VERBOSE, DEBUG +# from Gaudi.Configuration import * -from Gaudi.Configuration import * import os -from GaudiKernel.SystemOfUnits import GeV, tesla +from GaudiKernel.SystemOfUnits import GeV, tesla, mm +from GaudiKernel.PhysicalConstants import pi, halfpi, twopi +from math import cos, sin, tan use_pythia = False addNoise = False @@ -46,23 +51,30 @@ # (in strips: 0.5625/4=0.14) # Nevts = 20000 -Nevts = 100 +Nevts = 10 # Nevts = 1 # Nevts=1000 -# momentum = 100 # in GeV -# momentum = 50 # in GeV -momentum = 10 # in GeV -_pi = 3.14159 -thetaMin = 40 # degrees -thetaMax = 140 # degrees + +# particle momentum and direction +# momentum = 100 # in GeV +momentum = 50 # in GeV +# momentum = 10 # in GeV +thetaMin = 40 # degrees +thetaMax = 140 # degrees # thetaMin = 89 # thetaMax = 91 -# thetaMin = 90 # degrees -# thetaMax = 90 # degrees -# phiMin = _pi/2. -# phiMax = _pi/2. +# thetaMin = 90 # degrees +# thetaMax = 90 # degrees +# phiMin = halfpi +# phiMax = halfpi phiMin = 0 -phiMax = 2 * _pi +phiMax = twopi + +# particle origin +# origR = 1000.0*mm +origR = 0.0 * mm +origTheta = halfpi +origPhi = 0.0 # particle type: 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ pdgCode = 11 @@ -95,14 +107,6 @@ pythia8gentool.printPythiaStatistics = False pythia8gentool.pythiaExtraSettings = [""] genAlg.SignalProvider = pythia8gentool - # to smear the primary vertex position: - # from Configurables import GaussSmearVertex - # smeartool = GaussSmearVertex() - # smeartool.xVertexSigma = 0.5*units.mm - # smeartool.yVertexSigma = 0.5*units.mm - # smeartool.zVertexSigma = 40.0*units.mm - # smeartool.tVertexSigma = 180.0*units.picosecond - # genAlg.VertexSmearingTool = smeartool else: from Configurables import MomentumRangeParticleGun pgun = MomentumRangeParticleGun("ParticleGun") @@ -111,12 +115,30 @@ pgun.MomentumMax = momentum * GeV pgun.PhiMin = phiMin pgun.PhiMax = phiMax - pgun.ThetaMin = thetaMin * _pi / 180. - pgun.ThetaMax = thetaMax * _pi / 180. + pgun.ThetaMin = thetaMin * pi / 180. + pgun.ThetaMax = thetaMax * pi / 180. genAlg.SignalProvider = pgun genAlg.hepmc.Path = "hepmc" +# smear/shift vertex +if origR > 0.0: + origX = origR * cos(origPhi) + origY = origR * sin(origPhi) + origZ = origR / tan(origTheta) + print("Particle gun will be moved to %f %f %f" % (origX, origY, origZ)) + from Configurables import GaussSmearVertex + vertexSmearAndShiftTool = GaussSmearVertex() + vertexSmearAndShiftTool.xVertexSigma = 0. + vertexSmearAndShiftTool.yVertexSigma = 0. + vertexSmearAndShiftTool.tVertexSigma = 0. + vertexSmearAndShiftTool.xVertexMean = origX + vertexSmearAndShiftTool.yVertexMean = origY + vertexSmearAndShiftTool.zVertexMean = origZ + vertexSmearAndShiftTool.tVertexMean = 0. + genAlg.VertexSmearingTool = vertexSmearAndShiftTool + +# hepMC -> EDM converter hepmc_converter = HepMCToEDMConverter() hepmc_converter.hepmc.Path = "hepmc" genParticlesOutputName = "genParticles" @@ -131,7 +153,6 @@ path_to_detector = os.environ.get("K4GEO", "") print(path_to_detector) detectors_to_use = [ - # 'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml', 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml' ] # prefix all xmls with path_to_detector @@ -175,7 +196,7 @@ if magneticField == 1: field = SimG4ConstantMagneticFieldTool( "SimG4ConstantMagneticFieldTool", - FieldComponentZ=-2*tesla, + FieldComponentZ=-2 * tesla, FieldOn=True, IntegratorStepper="ClassicalRK4" ) @@ -197,9 +218,11 @@ # HCAL if runHCal: hcalBarrelReadoutName = "HCalBarrelReadout" + hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" hcalEndcapReadoutName = "HCalEndcapReadout" else: hcalBarrelReadoutName = "" + hcalBarrelReadoutName2 = "" hcalEndcapReadoutName = "" # Configure saving of calorimeter hits @@ -357,7 +380,7 @@ if runHCal: # Create cells in HCal - # 1. step - merge hits into cells with the default readout + # 1 - merge hits into cells with the default readout hcalBarrelCellsName = "HCalBarrelCells" createHcalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", doCellCalibration=True, @@ -369,22 +392,14 @@ cells=hcalBarrelCellsName, OutputLevel=INFO) - # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", - # doCellCalibration=True, - # calibTool=calibHcalEndcap, - # addCellNoise=False, - # filterCellNoise=False, - # OutputLevel=INFO) - # createHcalEndcapCells.hits.Path="HCalEndcapHits" - # createHcalEndcapCells.cells.Path="HCalEndcapCells" - + # 2 - attach positions to the cells from Configurables import CellPositionsHCalBarrelPhiThetaSegTool cellPositionHcalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( "CellPositionsHCalBarrel", readoutName=hcalBarrelReadoutName, OutputLevel=INFO ) - hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" + hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" createHcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( "CreateHcalBarrelPositionedCells", OutputLevel=INFO @@ -392,31 +407,75 @@ createHcalBarrelPositionedCells.positionsTool = cellPositionHcalBarrelTool createHcalBarrelPositionedCells.hits.Path = hcalBarrelCellsName createHcalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName + + # 3 - compute new cellID of cells based on new readout - removing row information + hcalBarrelCellsName2 = "HCalBarrelCells2" + rewriteHCalBarrel = RewriteBitfield("RewriteHCalBarrel", + # old bitfield (readout) + oldReadoutName=hcalBarrelReadoutName, + # specify which fields are going to be deleted + removeIds=["row"], + # new bitfield (readout), with new segmentation + newReadoutName=hcalBarrelReadoutName2, + debugPrint=10, + OutputLevel=INFO) + # clusters are needed, with deposit position and cellID in bits + rewriteHCalBarrel.inhits.Path = hcalBarrelCellsName + rewriteHCalBarrel.outhits.Path = hcalBarrelCellsName2 + + # 4 - attach positions to the new cells + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool + hcalBarrelPositionedCellsName2 = "HCalBarrelPositionedCells2" + cellPositionHcalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( + "CellPositionsHCalBarrel2", + readoutName=hcalBarrelReadoutName2, + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( + "CreateHCalBarrelPositionedCells2", + OutputLevel=INFO + ) + createHcalBarrelPositionedCells2.positionsTool = cellPositionHcalBarrelTool2 + createHcalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 + createHcalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 + + # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", + # doCellCalibration=True, + # calibTool=calibHcalEndcap, + # addCellNoise=False, + # filterCellNoise=False, + # OutputLevel=INFO) + # createHcalEndcapCells.hits.Path="HCalEndcapHits" + # createHcalEndcapCells.cells.Path="HCalEndcapCells" + else: hcalBarrelCellsName = "emptyCaloCells" hcalBarrelPositionedCellsName = "emptyCaloCells" + hcalBarrelCellsName2 = "emptyCaloCells" + hcalBarrelPositionedCellsName2 = "emptyCaloCells" cellPositionHcalBarrelTool = None - + cellPositionHcalBarrelTool2 = None + # Empty cells for parts of calorimeter not implemented yet createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") createemptycells.cells.Path = "emptyCaloCells" -# PRODUCE SLIDING WINDOW CLUSTERS +# Produce sliding window clusters (ECAL only) towers = CaloTowerToolFCCee("towers", - deltaThetaTower = 0.009817477, deltaPhiTower = 2*2*_pi/1536., - ecalBarrelReadoutName=ecalBarrelReadoutName, - ecalEndcapReadoutName=ecalEndcapReadoutName, - ecalFwdReadoutName="", - hcalBarrelReadoutName=hcalBarrelReadoutName, - hcalExtBarrelReadoutName="", - hcalEndcapReadoutName="", - hcalFwdReadoutName="", - OutputLevel=INFO) + deltaThetaTower = 0.009817477, deltaPhiTower = 2*2*pi/1536., + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalFwdReadoutName="", + hcalBarrelReadoutName="", + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName towers.ecalEndcapCells.Path = "ECalEndcapCells" towers.ecalFwdCells.Path = "emptyCaloCells" -towers.hcalBarrelCells.Path = hcalBarrelCellsName +towers.hcalBarrelCells.Path = "emptyCaloCells" towers.hcalExtBarrelCells.Path = "emptyCaloCells" towers.hcalEndcapCells.Path = "emptyCaloCells" towers.hcalFwdCells.Path = "emptyCaloCells" @@ -434,41 +493,34 @@ threshold = 0.040 createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", - towerTool=towers, - nThetaWindow=windT, nPhiWindow=windP, - nThetaPosition=posT, nPhiPosition=posP, - nThetaDuplicates=dupT, nPhiDuplicates=dupP, - nThetaFinal=finT, nPhiFinal=finP, - energyThreshold=threshold, - energySharingCorrection=False, - attachCells=True, - OutputLevel=INFO - ) + towerTool=towers, + nThetaWindow=windT, nPhiWindow=windP, + nThetaPosition=posT, nPhiPosition=posP, + nThetaDuplicates=dupT, nPhiDuplicates=dupP, + nThetaFinal=finT, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) createClusters.clusters.Path = "CaloClusters" createClusters.clusterCells.Path = "CaloClusterCells" -createEcalBarrelPositionedCaloClusterCells = CreateCaloCellPositionsFCCee( - "ECalBarrelPositionedCaloClusterCells", - OutputLevel=INFO -) -createEcalBarrelPositionedCaloClusterCells.positionsTool = cellPositionEcalBarrelTool -createEcalBarrelPositionedCaloClusterCells.hits.Path = "CaloClusterCells" -createEcalBarrelPositionedCaloClusterCells.positionedHits.Path = "PositionedCaloClusterCells" - correctCaloClusters = CorrectCaloClusters("correctCaloClusters", inClusters=createClusters.clusters.Path, - outClusters="Corrected"+createClusters.clusters.Path, + outClusters="Corrected" + createClusters.clusters.Path, numLayers=[12], firstLayerIDs=[0], lastLayerIDs=[11], - # readoutNames = [ecalBarrelReadoutNamePhiEta], readoutNames=[ecalBarrelReadoutName], - # upstreamParameters = [[0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + # upstreamParameters = [ + # [0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], upstreamParameters=[ [0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], upstreamFormulas=[ ['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], - # downstreamParameters = [[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + # downstreamParameters = [ + # [-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], downstreamParameters=[ [-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], downstreamFormulas=[ @@ -481,33 +533,35 @@ ecalBarrelReadoutName=ecalBarrelReadoutName, ecalEndcapReadoutName="", ecalFwdReadoutName="", - hcalBarrelReadoutName=hcalBarrelReadoutName, + hcalBarrelReadoutName=hcalBarrelReadoutName2, hcalExtBarrelReadoutName="", hcalEndcapReadoutName="", hcalFwdReadoutName="", OutputLevel=INFO) createTopoInput.ecalBarrelCells.Path = ecalBarrelPositionedCellsName -# createTopoInput.ecalBarrelCells.Path = "ECalBarrelPositionedCells2" createTopoInput.ecalEndcapCells.Path = "emptyCaloCells" createTopoInput.ecalFwdCells.Path = "emptyCaloCells" -createTopoInput.hcalBarrelCells.Path = hcalBarrelPositionedCellsName +createTopoInput.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 createTopoInput.hcalExtBarrelCells.Path = "emptyCaloCells" createTopoInput.hcalEndcapCells.Path = "emptyCaloCells" createTopoInput.hcalFwdCells.Path = "emptyCaloCells" cellPositionHcalBarrelNoSegTool = None cellPositionHcalExtBarrelTool = None +neighboursMap = "/LAr_scripts/data/neighbours_map_barrel_thetamodulemerged.root" +noiseMap = "/LAr_scripts/data/cellNoise_map_electronicsNoiseLevel_thetamodulemerged.root" +if runHCal: + neighboursMap = "/LAr_scripts/data/neighbours_map_ecalB_thetamodulemerged_hcalB_thetaphi.root" + noiseMap = "/LAr_scripts/data/cellNoise_map_electronicsNoiseLevel_ecalB_thetamodulemerged_hcalB_thetaphi.root" + readNeighboursMap = TopoCaloNeighbours("ReadNeighboursMap", - fileName=os.environ['FCCBASEDIR'] + - "/LAr_scripts/data/neighbours_map_barrel_thetamodulemerged.root", - #"/LAr_scripts/data/neighbours_map_HCalBarrel.root", + fileName=os.environ['FCCBASEDIR'] + neighboursMap, OutputLevel=INFO) # Noise levels per cell readNoisyCellsMap = TopoCaloNoisyCells("ReadNoisyCellsMap", - fileName=os.environ['FCCBASEDIR'] + - "/LAr_scripts/data/cellNoise_map_electronicsNoiseLevel_thetamodulemerged.root", + fileName=os.environ['FCCBASEDIR'] + noiseMap, OutputLevel=INFO) createTopoClusters = CaloTopoClusterFCCee("CreateTopoClusters", @@ -518,9 +572,9 @@ noiseTool=readNoisyCellsMap, # cell positions tools for all sub - systems positionsECalBarrelTool=cellPositionEcalBarrelTool, - positionsHCalBarrelTool=cellPositionHcalBarrelTool, - positionsHCalBarrelNoSegTool=cellPositionHcalBarrelNoSegTool, - positionsHCalExtBarrelTool=cellPositionHcalExtBarrelTool, + positionsHCalBarrelTool=cellPositionHcalBarrelTool2, + # positionsHCalBarrelNoSegTool=cellPositionHcalBarrelNoSegTool, + # positionsHCalExtBarrelTool=cellPositionHcalExtBarrelTool, # positionsHCalExtBarrelTool = HCalExtBcells, # positionsEMECTool = EMECcells, # positionsHECTool = HECcells, @@ -534,31 +588,25 @@ createTopoClusters.clusters.Path = "CaloTopoClusters" createTopoClusters.clusterCells.Path = "CaloTopoClusterCells" -createEcalBarrelPositionedCaloTopoClusterCells = CreateCaloCellPositionsFCCee( - "ECalBarrelPositionedCaloTopoClusterCells", - OutputLevel=INFO -) -# createEcalBarrelPositionedCaloTopoClusterCells.positionsTool = cellPositionEcalBarrelTool2 -createEcalBarrelPositionedCaloTopoClusterCells.positionsTool = cellPositionEcalBarrelTool -createEcalBarrelPositionedCaloTopoClusterCells.hits.Path = "CaloTopoClusterCells" -createEcalBarrelPositionedCaloTopoClusterCells.positionedHits.Path = "PositionedCaloTopoClusterCells" - +# Correction below is for EM-only clusters +# What to do for EM+HAD topoclusters? correctCaloTopoClusters = CorrectCaloClusters( "correctCaloTopoClusters", inClusters=createTopoClusters.clusters.Path, - outClusters="Corrected"+createTopoClusters.clusters.Path, + outClusters="Corrected" + createTopoClusters.clusters.Path, numLayers=[12], firstLayerIDs=[0], lastLayerIDs=[11], - # readoutNames = [ecalBarrelReadoutNamePhiEta], readoutNames=[ecalBarrelReadoutName], - # upstreamParameters = [[0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + # upstreamParameters = [[0.02729094887360858, -1.378665489864182, -68.40424543618059, + # 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], upstreamParameters=[[0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], upstreamFormulas=[['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], - # downstreamParameters = [[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], - downstreamParameters=[[-0.0032351643028483354, 0.006597484738888312, - 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + # downstreamParameters = [[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, + # -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + downstreamParameters=[[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, + -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], downstreamFormulas=[['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], OutputLevel=INFO ) @@ -568,8 +616,11 @@ OutputLevel=INFO) # out.outputCommands = ["keep *"] -# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells"] -# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells", "drop %s"%ecalBarrelCellsName, "drop %s"%createEcalBarrelPositionedCells.positionedHits.Path] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", +# "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells"] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", +# "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells", +# "drop %s" % ecalBarrelCellsName, "drop %s" % createEcalBarrelPositionedCells.positionedHits.Path] # out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop *ells*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells"] # out.outputCommands = ["keep *", "drop HCal*", "drop ECalBarrel*", "drop emptyCaloCells"] if runHCal: @@ -618,15 +669,15 @@ TopAlg += [ createHcalBarrelCells, createHcalBarrelPositionedCells, + rewriteHCalBarrel, + createHcalBarrelPositionedCells2, # createHcalEndcapCells ] TopAlg += [ createemptycells, createClusters, - createEcalBarrelPositionedCaloClusterCells, correctCaloClusters, createTopoClusters, - createEcalBarrelPositionedCaloTopoClusterCells, correctCaloTopoClusters, out ]