Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] First DECal / SiW ECal commit #294

Open
wants to merge 32 commits into
base: master
Choose a base branch
from
Open
Changes from 1 commit
Commits
Show all changes
32 commits
Select commit Hold shift + click to select a range
8479e80
toprice first commit of DECAL geometry
Nov 17, 2016
2fd40d7
first commit of the DECAL before updating master FCCSW
Nov 25, 2016
ba42ad9
Merge remote branch 'hep-fcc/master'
Nov 25, 2016
ba49172
DECAL Geometry first commit
Apr 5, 2017
e838711
SD classes for DECAL
Apr 5, 2017
acd1cad
Redo Segmentation of DECAL pixel to Pad
Apr 5, 2017
fd3674e
Merge with v0.8 SDWrapper.cpp
Apr 5, 2017
e5210b8
added back in geant_fullsim_ecal_SPG_new.py
Apr 7, 2017
f8feb3e
Merge remote branch 'hep-fcc/master'
Apr 27, 2017
2e12c47
Added basic plotting for the DECAL into Detector/DetectorStudies and …
Apr 28, 2017
77fcad1
Added DEcal analysis to count pixels and analyse variables per layer.…
Sep 11, 2017
ae03431
file to calculate the non linearity factor from total counts in an event
Sep 28, 2017
77a19a6
New xml files to run epitaxial and subtrate (analogue) Silicon ECAL
Nov 21, 2017
cc7e7a2
Analysis code to extract analogue component and analyse
Nov 21, 2017
e627926
code which creates resolution plot from SiWAnalysis
Nov 21, 2017
80030f6
code which creates resoltuion plot from DECalAnalysis
Nov 21, 2017
12bb45c
code to submit analogue analysis to lxbatch
Nov 21, 2017
5b42744
Previous commit did not take into account changes to pre-existing code
Nov 21, 2017
40b9222
Added in flags to DetSegmentation to allow extraction of Digital or n…
Nov 22, 2017
b671bef
Many many changes....
Feb 28, 2018
2f7de76
added many things....
Feb 28, 2018
f408156
Merge remote branch 'hep-fcc/master'
Feb 28, 2018
1361ac2
Modified DECal geometry to match with DD4hep-->dd4hep and removed nam…
Mar 2, 2018
47b2b43
Fixed DD4hep depenacies
Mar 4, 2018
06f4b22
modified DD4hep/LCCD.h for DD4hep/Detector.h and removed all commente…
Mar 6, 2018
d295777
Changing xml child levels
Mar 20, 2018
a3e14ae
Modified XML file to take childs of calorimeter
Mar 20, 2018
8fa4bb7
Removed all old xml files
Mar 20, 2018
b5dfa6a
Fixed xml file
Mar 20, 2018
57c968c
Actually fixed XML file this time
Mar 20, 2018
b047b93
added scripts to submit pile up events
Mar 21, 2018
e6c1161
modified pile-up files and hacked geometry to just draw a single stave
Apr 7, 2018
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Added basic plotting for the DECAL into Detector/DetectorStudies and …
…added script to run a simulation then reconstruct
tonypba committed Apr 28, 2017
commit 2e12c474e8e7fb168feb61715ae881eb5c89115c
127 changes: 127 additions & 0 deletions Detector/DetStudies/src/components/DECalLongitudinalTest.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,127 @@
#include "SamplingFractionInLayers.h"

// FCCSW
#include "DetInterface/IGeoSvc.h"

// datamodel
#include "datamodel/PositionedCaloHitCollection.h"

#include "CLHEP/Vector/ThreeVector.h"
#include "GaudiKernel/ITHistSvc.h"
#include "TH1F.h"
#include "TVector2.h"

// DD4hep
#include "DD4hep/LCDD.h"
#include "DD4hep/Readout.h"

DECLARE_ALGORITHM_FACTORY(SamplingFractionInLayers)

SamplingFractionInLayers::SamplingFractionInLayers(const std::string& aName, ISvcLocator* aSvcLoc)
: GaudiAlgorithm(aName, aSvcLoc),
m_histSvc("THistSvc", "SamplingFractionInLayers"),
m_geoSvc("GeoSvc", "SamplingFractionInLayers"),
m_totalEnergy(nullptr),
m_totalActiveEnergy(nullptr),
m_sf(nullptr) {
declareProperty("deposits", m_deposits, "Energy deposits in sampling calorimeter (input)");
}
SamplingFractionInLayers::~SamplingFractionInLayers() {}

StatusCode SamplingFractionInLayers::initialize() {
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
// check if readouts exist
if (m_geoSvc->lcdd()->readouts().find(m_readoutName) == m_geoSvc->lcdd()->readouts().end()) {
error() << "Readout <<" << m_readoutName << ">> does not exist." << endmsg;
return StatusCode::FAILURE;
}
// create histograms
for (uint i = 0; i < m_numLayers; i++) {
m_totalEnLayers.push_back(new TH1F(("ecal_totalEnergy_layer" + std::to_string(i)).c_str(),
("Total deposited energy in layer " + std::to_string(i)).c_str(), 1000, 0,
1.2 * m_energy));
if (m_histSvc->regHist("/rec/ecal_total_layer" + std::to_string(i), m_totalEnLayers.back()).isFailure()) {
error() << "Couldn't register histogram" << endmsg;
return StatusCode::FAILURE;
}
m_activeEnLayers.push_back(new TH1F(("ecal_activeEnergy_layer" + std::to_string(i)).c_str(),
("Deposited energy in active material, in layer " + std::to_string(i)).c_str(),
1000, 0, 1.2 * m_energy));
if (m_histSvc->regHist("/rec/ecal_active_layer" + std::to_string(i), m_activeEnLayers.back()).isFailure()) {
error() << "Couldn't register histogram" << endmsg;
return StatusCode::FAILURE;
}
m_sfLayers.push_back(new TH1F(("ecal_sf_layer" + std::to_string(i)).c_str(),
("SF for layer " + std::to_string(i)).c_str(), 1000, 0, 1));
if (m_histSvc->regHist("/rec/ecal_sf_layer" + std::to_string(i), m_sfLayers.back()).isFailure()) {
error() << "Couldn't register histogram" << endmsg;
return StatusCode::FAILURE;
}
}
m_totalEnergy = new TH1F("ecal_totalEnergy", "Total deposited energy", 1000, 0, 1.2 * m_energy);
if (m_histSvc->regHist("/rec/ecal_total", m_totalEnergy).isFailure()) {
error() << "Couldn't register histogram" << endmsg;
return StatusCode::FAILURE;
}
m_totalActiveEnergy = new TH1F("ecal_active", "Deposited energy in active material", 1000, 0, 1.2 * m_energy);
if (m_histSvc->regHist("/rec/ecal_active", m_totalActiveEnergy).isFailure()) {
error() << "Couldn't register histogram" << endmsg;
return StatusCode::FAILURE;
}
m_sf = new TH1F("ecal_sf", "Sampling fraction", 1000, 0, 1);
if (m_histSvc->regHist("/rec/ecal_sf", m_sf).isFailure()) {
error() << "Couldn't register histogram" << endmsg;
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}

StatusCode SamplingFractionInLayers::execute() {
auto decoder = m_geoSvc->lcdd()->readout(m_readoutName).idSpec().decoder();
double sumE = 0.;
std::vector<double> sumElayers;
double sumEactive = 0.;
std::vector<double> sumEactiveLayers;
sumElayers.assign(m_numLayers, 0);
sumEactiveLayers.assign(m_numLayers, 0);

const auto deposits = m_deposits.get();
for (const auto& hit : *deposits) {
sumElayers[(*decoder)[m_layerFieldName]] += hit.core().energy;
// check if energy was deposited in the calorimeter (active/passive material)
// layers are numbered starting from 1, layer == 0 is cryostat/bath
if ((*decoder)[m_layerFieldName] > 0) {
sumE += hit.core().energy;
decoder->setValue(hit.core().cellId);
// active material of calorimeter
if ((*decoder)[m_activeFieldName] == m_activeFieldValue) {
sumEactive += hit.core().energy;
sumEactiveLayers[(*decoder)[m_layerFieldName]] += hit.core().energy;
}
}
}
// Fill histograms
m_totalEnergy->Fill(sumE);
m_totalActiveEnergy->Fill(sumEactive);
if (sumE > 0) {
m_sf->Fill(sumEactive / sumE);
}
for (uint i = 0; i < m_numLayers; i++) {
m_totalEnLayers[i]->Fill(sumElayers[i]);
m_activeEnLayers[i]->Fill(sumEactiveLayers[i]);
if (i == 0) {
debug() << "total energy deposited in cryostat and bath = " << sumElayers[i] << endmsg;
} else {
debug() << "total energy in layer " << i << " = " << sumElayers[i] << " active = " << sumEactiveLayers[i]
<< endmsg;
}
if (sumElayers[i] > 0) {
m_sfLayers[i]->Fill(sumEactiveLayers[i] / sumElayers[i]);
}
}
return StatusCode::SUCCESS;
}

StatusCode SamplingFractionInLayers::finalize() { return GaudiAlgorithm::finalize(); }
82 changes: 82 additions & 0 deletions Detector/DetStudies/src/components/DECalLongitudinalTest.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
#ifndef DETSTUDIES_SAMPLINGFRACTIONINLAYERS_H
#define DETSTUDIES_SAMPLINGFRACTIONINLAYERS_H

// GAUDI
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ServiceHandle.h"

// FCCSW
#include "FWCore/DataHandle.h"
class IGeoSvc;

// datamodel
namespace fcc {
class PositionedCaloHitCollection;
}

class TH1F;
class ITHistSvc;
/** @class SamplingFractionInLayers SamplingFractionInLayers.h
*
* Histograms of energy deposited in active material and total energy deposited in the calorimeter.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update.

* Passive material needs to be marked as sensitive. It needs to be divided into layers (cells) as active material.
* Layers (cells) are numbered starting at 1 so that energy depoisted in cryostat and bath could be easily recognised.
* Sampling fraction is calculated for each layer as the ratio of energy deposited in active material to energy
* deposited in the layer (also in passive material).
*
* @author Anna Zaborowska
*/

class SamplingFractionInLayers : public GaudiAlgorithm {
public:
explicit SamplingFractionInLayers(const std::string&, ISvcLocator*);
virtual ~SamplingFractionInLayers();
/** Initialize.
* @return status code
*/
virtual StatusCode initialize() final;
/** Fills the histograms.
* @return status code
*/
virtual StatusCode execute() final;
/** Finalize.
* @return status code
*/
virtual StatusCode finalize() final;

private:
/// Pointer to the interface of histogram service
ServiceHandle<ITHistSvc> m_histSvc;
/// Pointer to the geometry service
ServiceHandle<IGeoSvc> m_geoSvc;
/// Handle for the energy deposits
DataHandle<fcc::PositionedCaloHitCollection> m_deposits{"rec/caloHits", Gaudi::DataHandle::Reader, this};
/// Name of the active field
Gaudi::Property<std::string> m_activeFieldName{this, "activeFieldName", "", "Identifier of active material"};
/// Value of the active material
Gaudi::Property<int> m_activeFieldValue{this, "activeFieldValue", 0, "Value of identifier for active material"};
/// Name of the layer/cell field
Gaudi::Property<std::string> m_layerFieldName{this, "layerFieldName", "", "Identifier of layers"};
/// Number of layers/cells
Gaudi::Property<uint> m_numLayers{this, "numLayers", 8, "Number of layers"};
/// Name of the detector readout
Gaudi::Property<std::string> m_readoutName{this, "readoutName", "", "Name of the detector readout"};
// Maximum energy for the axis range
Gaudi::Property<double> m_energy{this, "energyAxis", 500, "Maximum energy for axis range"};
// Histograms of total deposited energy within layer
// Layers are numbered starting at 1. Layer 0 includes total energy deposited in cryostat and bath (in front and
// behind calo)
std::vector<TH1F*> m_totalEnLayers;
// Histogram of total deposited energy in the calorimeter (in active and passive material, excluding cryostat and
// bath)
TH1F* m_totalEnergy;
// Histograms of energy deposited in the active material within layer
std::vector<TH1F*> m_activeEnLayers;
// Histogram of energy deposited in the active material of the calorimeter
TH1F* m_totalActiveEnergy;
// Histograms of sampling fraction (active/total energy) calculated within layer
std::vector<TH1F*> m_sfLayers;
// Histogram of sampling fraction (active/total energy) calculated for the calorimeter (excluding cryostat and bath)
TH1F* m_sf;
};
#endif /* DETSTUDIES_SAMPLINGFRACTIONINLAYERS_H */
83 changes: 83 additions & 0 deletions Detector/DetStudies/tests/options/longitudinalTest_DEcal.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
from Gaudi.Configuration import *
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You may add some tests to CMakeLists.


# Data service
from Configurables import FCCDataSvc
podioevent = FCCDataSvc("EventDataSvc")

# DD4hep geometry service
from Configurables import GeoSvc
geoservice = GeoSvc("GeoSvc", detectors=[ 'file:Detector/DetFCChhBaseline1/compact/FCChh_DectEmptyMaster.xml',
'file:Detector/DetFCChhECalDigital/compact/FCChh_DECalBarrel_Mockup.xml'
],
OutputLevel = INFO)

# Geant4 service
# Configures the Geant simulation: geometry, physics list and user actions
from Configurables import SimG4Svc
geantservice = SimG4Svc("SimG4Svc", detector='SimG4DD4hepDetector', physicslist="SimG4FtfpBert", actions="SimG4FullSimActions")
geantservice.G4commands += ["/run/setCut 0.1 mm"]

# Geant4 algorithm
# Translates EDM to G4Event, passes the event to G4, writes out outputs via tools
# and a tool that saves the calorimeter hits
from Configurables import SimG4Alg, SimG4SaveCalHits, SimG4SingleParticleGeneratorTool
saveecaltool = SimG4SaveCalHits("saveECalHits",readoutNames = ["BarDECal_Readout"])
#saveecaltool.DataOutputs.caloClusters.Path = "ECalClusters"
saveecaltool.positionedCaloHits.Path = "positionedCaloHits"
saveecaltool.caloHits.Path = "ECalHits"

from Configurables import SimG4SingleParticleGeneratorTool
pgun=SimG4SingleParticleGeneratorTool("SimG4SingleParticleGeneratorTool",
saveEdm=True,
particleName="e-",
energyMin=100000,energyMax=100000,
etaMin=0,etaMax=0,
# phiMin=0, phiMax=0.001,
vertexX=0,vertexY=0,vertexZ=0,
OutputLevel =DEBUG)

# next, create the G4 algorithm, giving the list of names of tools ("XX/YY")
geantsim = SimG4Alg("SimG4Alg",
outputs= ["SimG4SaveCalHits/saveECalHits"],
eventProvider = pgun,
OutputLevel = DEBUG)

from Configurables import DECalLongitudinalTest
hist = DECalLongitudinalTest("hists",
readoutName = "BarDECal_Readout",
layerFieldName = "layer",
numLayers = 50, # one more because index starts at 1 - layer 0 will be always empty
OutputLevel = DEBUG)
hist.deposits.Path="positionedCaloHits"

THistSvc().Output = ["rec DATAFILE='hist_test.root' TYP='ROOT' OPT='RECREATE'"]
THistSvc().PrintAll=True
THistSvc().AutoSave=True
THistSvc().AutoFlush=False
THistSvc().OutputLevel=INFO


# PODIO algorithm
from Configurables import PodioOutput
out = PodioOutput("out",
OutputLevel=DEBUG)
out.outputCommands = ["keep *"]
out.filename = "output_10GeV_1mmAir_0T_SiAirW.root"

#CPU information
from Configurables import AuditorSvc, ChronoAuditor
chra = ChronoAuditor()
audsvc = AuditorSvc()
audsvc.Auditors = [chra]
geantsim.AuditExecute = True
hist.AuditExecute = True

# ApplicationMgr
from Configurables import ApplicationMgr
ApplicationMgr( TopAlg = [geantsim, out, hist],
EvtSel = 'NONE',
EvtMax = 10,
# order is important, as GeoSvc is needed by G4SimSvc
ExtSvc = [podioevent, geoservice, geantservice, audsvc],
OutputLevel = DEBUG
)