Skip to content

Commit

Permalink
Merge pull request #19 from giovannimarchiori/gmarchio-main-20231215-…
Browse files Browse the repository at this point in the history 
…ecalhcalclustering

Combined ECAL+HCAL topoclustering
  • Loading branch information
@BrieucF
BrieucF authored Feb 12, 2024
2 parents 7e03832 + 0ec80e8 commit a6c89a5
Show file tree
Hide file tree
Showing 4 changed files with 351 additions and 186 deletions.
57 changes: 43 additions & 14 deletions FCCSW_ecal/neighbours_theta.py
View file
Original file line number Diff line number Diff line change
@@ -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", "")
Expand All @@ -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=[],
Expand Down
147 changes: 95 additions & 52 deletions FCCSW_ecal/noise_map_theta.py
View file
Original file line number Diff line number Diff line change
@@ -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
)
80 changes: 61 additions & 19 deletions FCCSW_ecal/printNeighbours.C
View file
Original file line number Diff line number Diff line change
@@ -1,9 +1,23 @@
#include "TTree.h"
#include "TFile.h"
#include "TRandom.h"
#include <iostream>
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<unsigned long> *neighbours = 0;
ULong64_t cIDNoise;
double noiseLevel;
double noiseOffset;

bool useHCalReadoutWithRows = false;

Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand All @@ -131,7 +173,7 @@ void printNeighboursOfCell(ULong_t cellID)
T->GetEntry(iEntry);
if (cID == cellID)
{
printCellAndNeighbours(iEntry);
printCell(iEntry, true, false);
return;
}
}
Expand All @@ -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);
}
}
Loading

0 comments on commit a6c89a5

Please sign in to comment.