diff --git a/Detector/DetCommon/compact/HcalBarrelAir.xml b/Detector/DetCommon/compact/HcalBarrelAir.xml
new file mode 100644
index 000000000..6dc41bd68
--- /dev/null
+++ b/Detector/DetCommon/compact/HcalBarrelAir.xml
@@ -0,0 +1,51 @@
+
+
+
+
+ HCal Place-Holder
+
+
+
+
+
+
+
+
+
+
+
+
+ Envelope for HCal barrel
+
+
+
+
+
+
+
+
diff --git a/Detector/DetCommon/compact/TrackerBarrelAir.xml b/Detector/DetCommon/compact/TrackerBarrelAir.xml
new file mode 100644
index 000000000..b0e8614f9
--- /dev/null
+++ b/Detector/DetCommon/compact/TrackerBarrelAir.xml
@@ -0,0 +1,37 @@
+
+
+
+
+ Tracker Place-Holder
+
+
+
+
+
+
+
+
+
+ Envelope for Tracker
+
+
+
+
+
+
diff --git a/Detector/DetComponents/src/RedoSegmentation.cpp b/Detector/DetComponents/src/RedoSegmentation.cpp
index b22163226..97c76691b 100644
--- a/Detector/DetComponents/src/RedoSegmentation.cpp
+++ b/Detector/DetComponents/src/RedoSegmentation.cpp
@@ -17,6 +17,7 @@ DECLARE_ALGORITHM_FACTORY(RedoSegmentation)
RedoSegmentation::RedoSegmentation(const std::string& aName, ISvcLocator* aSvcLoc) : GaudiAlgorithm(aName, aSvcLoc) {
declareProperty("inhits", m_inHits, "Hit collection with old segmentation (input)");
declareProperty("outhits", m_outHits, "Hit collection with modified segmentation (output)");
+
}
RedoSegmentation::~RedoSegmentation() {}
@@ -83,11 +84,17 @@ StatusCode RedoSegmentation::execute() {
// cellID contains the volumeID that needs to be copied to the new id
uint64_t oldid = 0;
uint debugIter = 0;
+
+ int nhits = 0;
+
for (const auto& hit : *inHits) {
+
fcc::CaloHit newHit = outHits->create();
newHit.energy(hit.energy());
newHit.time(hit.time());
- m_oldDecoder->setValue(hit.cellId());
+
+ nhits++;
+
if (debugIter < m_debugPrint) {
debug() << "OLD: " << m_oldDecoder->valueString() << endmsg;
}
@@ -107,6 +114,7 @@ StatusCode RedoSegmentation::execute() {
debugIter++;
}
}
+ info() << "nhits = " << nhits << endmsg;
return StatusCode::SUCCESS;
}
diff --git a/Detector/DetComponents/src/RedoSegmentation.h b/Detector/DetComponents/src/RedoSegmentation.h
index cf1238ba1..382b4598c 100644
--- a/Detector/DetComponents/src/RedoSegmentation.h
+++ b/Detector/DetComponents/src/RedoSegmentation.h
@@ -83,5 +83,6 @@ class RedoSegmentation : public GaudiAlgorithm {
std::vector m_detectorIdentifiers;
/// Limit of debug printing
Gaudi::Property m_debugPrint{this, "debugPrint", 10, "Limit of debug printing"};
+
};
#endif /* DETCOMPONENTS_REDOSEGMENTATION_H */
diff --git a/Detector/DetComponents/tests/options/redoSegmentationXYZ_DEcal.py b/Detector/DetComponents/tests/options/redoSegmentationXYZ_DEcal.py
new file mode 100644
index 000000000..2638a7db8
--- /dev/null
+++ b/Detector/DetComponents/tests/options/redoSegmentationXYZ_DEcal.py
@@ -0,0 +1,61 @@
+from Gaudi.Configuration import *
+from Configurables import ApplicationMgr
+
+from Configurables import FCCDataSvc
+inputfile = "/afs/cern.ch/user/t/toprice/private/FCC/FCCSW/batch_eos/50Layers_2.1mmW_50umPixels_18umThick_FCCSW0.8/100GeV_BFIELD4T_ETAMIN-0.001_ETAMAX0.001/output_100GeV_BFIELD4T_1.root"
+podiosvc = FCCDataSvc("EventDataSvc", input=inputfile)
+
+from Configurables import PodioInput
+podioinput = PodioInput("PodioReader", collections=["positionedCaloHits"], OutputLevel=DEBUG)
+
+from Configurables import GeoSvc
+geoservice = GeoSvc("GeoSvc", detectors=[ 'file:/afs/cern.ch/user/t/toprice/private/FCC/FCCSW/Detector/DetFCChhBaseline1/compact/FCChh_DectEmptyMaster.xml',
+ 'file:/afs/cern.ch/user/t/toprice/private/FCC/FCCSW/Detector/DetFCChhECalDigital/compact/FCChh_DECalBarrel_50Layers_2.1mmW_50umPixels_18umThick.xml'
+ ],
+ OutputLevel = DEBUG)
+
+from Configurables import RedoSegmentation
+resegment = RedoSegmentation("ReSegmentation",
+ # old bitfield (readout)
+ oldReadoutName = "BarDECal_Readout",
+ # specify which fields are going to be deleted
+ oldSegmentationIds = ["x","y","z"],
+ # new bitfield (readout), with new segmentation
+ newReadoutName="BarDECal_Pads",
+ OutputLevel = DEBUG)
+# clusters are needed, with deposit position and cellID in bits
+resegment.inhits.Path = "positionedCaloHits"
+resegment.outhits.Path = "newCaloHits"
+
+from Configurables import CreateCaloCells
+createcells = CreateCaloCells("CreateCaloCells",
+ doCellCalibration = False,
+ addCellNoise = False, filterCellNoise = False,
+ OutputLevel = DEBUG)
+createcells.hits.Path="newCaloHits"
+createcells.cells.Path="newCaloCells"
+
+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
+
+from Configurables import FCCDataSvc, PodioOutput
+#podiosvc = FCCDataSvc("EventDataSvc")
+podioout = PodioOutput("out", filename="/afs/cern.ch/user/t/toprice/private/FCC/FCCSW/batch_eos/50Layers_2.1mmW_50umPixels_18umThick_FCCSW0.8/100GeV_BFIELD4T_ETAMIN-0.001_ETAMAX0.001/testResegmentationXYZ_100GeV.root")
+podioout.outputCommands = ["keep *"]
+
+ApplicationMgr(EvtSel='NONE',
+ EvtMax=5000,
+ TopAlg=[podioinput, resegment, createcells,hist, podioout],
+ ExtSvc = [podiosvc, geoservice],
+ OutputLevel=DEBUG)
diff --git a/Detector/DetFCChhECalDigital/CMakeLists.txt b/Detector/DetFCChhECalDigital/CMakeLists.txt
new file mode 100644
index 000000000..6404a80c8
--- /dev/null
+++ b/Detector/DetFCChhECalDigital/CMakeLists.txt
@@ -0,0 +1,27 @@
+################################################################################
+# Package: DetFCChhECalDigital
+################################################################################
+gaudi_subdir(DetFCChhECalDigital v1r0)
+
+gaudi_depends_on_subdirs(GaudiKernel
+ Detector/DetExtensions)
+
+
+find_package(DD4hep)
+find_package(Geant4)
+include(${Geant4_USE_FILE})
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${DD4hep_ROOT}/cmake )
+include( DD4hep )
+
+find_package(ROOT COMPONENTS MathCore GenVector Geom REQUIRED)
+
+gaudi_add_module(DetFCChhECalDigital
+ src/*.cpp
+ INCLUDE_DIRS DD4hep ROOT DetExtensions Geant4
+ LINK_LIBRARIES GaudiKernel DD4hep ROOT Geant4)
+
+set(LIBRARY_OUTPUT_PATH ${CMAKE_LIBRARY_OUTPUT_DIRECTORY})
+message(STATUS "LIBRARY_OUTPUT_PATH -> ${LIBRARY_OUTPUT_PATH}")
+dd4hep_generate_rootmap(DetFCChhECalDigital)
+
diff --git a/Detector/DetFCChhECalDigital/compact/FCChh_DECalBarrel_50Layers_2.1mmW_50umPixels_18umThick_RRB.xml b/Detector/DetFCChhECalDigital/compact/FCChh_DECalBarrel_50Layers_2.1mmW_50umPixels_18umThick_RRB.xml
new file mode 100644
index 000000000..c6048ab9e
--- /dev/null
+++ b/Detector/DetFCChhECalDigital/compact/FCChh_DECalBarrel_50Layers_2.1mmW_50umPixels_18umThick_RRB.xml
@@ -0,0 +1,84 @@
+
+
+
+
+
+ DECal very conceptual design
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ system:3,active:1,EM_barrel:1,layer:6,digital:1,x:-17,y:-17,z:-18
+
+
+
+
+
+
+
+
+ system:3,active:1,EM_barrel:1,layer:6,digital:1,x:-15,y:-15,z:-18
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Detector/DetFCChhECalDigital/src/ECalBarrel_geo.cpp b/Detector/DetFCChhECalDigital/src/ECalBarrel_geo.cpp
new file mode 100644
index 000000000..50b374caa
--- /dev/null
+++ b/Detector/DetFCChhECalDigital/src/ECalBarrel_geo.cpp
@@ -0,0 +1,356 @@
+// DD4hep includes
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/Shapes.h"
+#include "DD4hep/Objects.h"
+
+#include "XML/XMLElements.h"
+
+// Gaudi
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IMessageSvc.h"
+#include "GaudiKernel/MsgStream.h"
+
+#include "cmath"
+
+
+namespace det {
+
+static dd4hep::detail::Ref_t createECal (dd4hep::Detector& lcdd,dd4hep::xml::Handle_t xmlElement,
+ dd4hep::SensitiveDetector sensDet)
+{
+
+ ServiceHandle msgSvc("MessageSvc", "DECalConstruction");
+ MsgStream lLog(&(*msgSvc), "DECalConstruction");
+
+ lLog << MSG::DEBUG << "++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endmsg;
+ lLog << MSG::DEBUG << "\t\t\t Building DECAL " << endmsg;
+ lLog << MSG::DEBUG << " +++++++++++++++++++++++++++++++++++++++++++++++++++++" << endmsg;
+
+ // Getting the detectors various values from the xml file
+
+ dd4hep::xml::DetElement xmlDet = xmlElement;
+ std::string detName = xmlDet.nameStr();
+ //Make DetElement
+ dd4hep::DetElement eCal(detName, xmlDet.id());
+
+
+
+ dd4hep::xml::DetElement calo = xmlDet.child(_Unicode(calorimeter));
+ dd4hep::xml::Dimension calo_dims(calo.dimensions());
+ std::string calo_name=calo.nameStr();
+ double calo_id=calo.id();
+
+ //dd4hep::xml::DetElement trkr = calo.child(_Unicode(tracker));
+ std::string trkr_mat = "Silicon";//trkr.materialStr();
+ double trkr_tck = 0;//trkr.thickness();
+
+ dd4hep::xml::DetElement active = calo.child(_Unicode(active));
+ std::string active_mat=active.materialStr();
+ double active_tck=active.thickness();
+
+
+ dd4hep::xml::DetElement samples_xml = calo.child(_Unicode(layers));
+ int active_samples = samples_xml.repeat();
+
+ dd4hep::xml::DetElement substrate = calo.child(_Unicode(substrate));
+ std::string substrate_mat=substrate.materialStr();
+ double substrate_tck=substrate.thickness();
+
+ dd4hep::xml::DetElement passive = calo.child(_Unicode(passive));
+ std::string passive_mat=passive.materialStr();
+ double passive_tck=passive.thickness();
+
+ dd4hep::xml::DetElement padding = calo.child(_Unicode(spacing));
+ std::string padding_mat=padding.materialStr();
+ double padding_tck=padding.thickness();
+
+
+
+
+ // here we take the fabs value of padding as if it is negative it goes before the W volume
+ double module_tck = active_tck+substrate_tck+passive_tck+std::fabs(padding_tck);
+ double calo_tck=(active_samples*module_tck) + trkr_tck;
+
+ //Introduce number of staves, modules from the xml file:
+ int nStaves = calo_dims.numsides();
+ int nModules = calo_dims.nmodules();
+ int testing = 1;
+
+ lLog << MSG::INFO << "Tracker Material Thickness = " << trkr_tck << endmsg;
+
+ lLog << MSG::INFO << "ECAL Epitaxial Thickness = " << active_tck << endmsg;
+ lLog << MSG::INFO << "ECAL Substrate Thickness = " << substrate_tck << endmsg;
+ lLog << MSG::INFO << "ECAL Passive Thickness = " << passive_tck << endmsg;
+ lLog << MSG::INFO << "ECAL Between Layers Thickness = " << padding_tck << endmsg;
+ lLog << MSG::INFO << "ECAL Module Thickness = " << module_tck << endmsg;
+
+ lLog << MSG::INFO << "nSamplings Layers = " << active_samples << endmsg;
+ lLog << MSG::INFO << "Total Calorimeter thickness = " << calo_tck << endmsg;
+
+ //Log number of staves/modules in the log file:
+ lLog << MSG::INFO << "nStaves = " << nStaves << endmsg;
+ lLog << MSG::INFO << "nModules = " << nModules << endmsg;
+
+ // Make volume that envelopes the whole barrel; set material to air
+ dd4hep::xml::Dimension dimensions(xmlDet.dimensions());
+ dd4hep::Tube envelopeShape(dimensions.rmin(), dimensions.rmax(), dimensions.dz());
+ dd4hep::Volume envelopeVolume(detName, envelopeShape, lcdd.air());
+ // Invisibility seems to be broken in visualisation tags, have to hardcode that
+ std::cout << "dimensions.visStr() = " << dimensions.visStr()<< std::endl;
+ envelopeVolume.setVisAttributes(lcdd, dimensions.visStr());
+
+ // create a lump of air the entire size of the barrel ECAL (Cylindrical)
+ // place it within the ECAL envelope
+ dd4hep::DetElement caloDet(calo_name, calo_id);
+ dd4hep::Tube caloShape(calo_dims.rmin() , calo_dims.rmax(), calo_dims.dz());
+ lLog << MSG::INFO << "ECAL: Building the calorimeter from " << calo_dims.rmin() << " to " << calo_dims.rmax() << endmsg;
+ dd4hep::Volume caloVolume(padding_mat, caloShape, lcdd.material(padding_mat));
+ lLog << MSG::DEBUG << "ECAL: Filling the calorimeter volume with " << padding_mat << endmsg;
+ dd4hep::PlacedVolume placedCalo = envelopeVolume.placeVolume(caloVolume);
+ placedCalo.addPhysVolID("EM_barrel", calo_id);
+ caloDet.setPlacement(placedCalo);
+
+ // set the sensitive detector type to the DD4hep calorimeter
+ sensDet.setType("DigitalCalorimeterSD");
+
+ // if we add material for the tracker then it needs to be placed first
+ if(trkr_tck > 0.0) {
+ dd4hep::DetElement trkrLayer(trkr_mat, 0);
+ dd4hep::Tube trkrShape(calo_dims.rmin(), calo_dims.rmin()+trkr_tck, calo_dims.dz());
+ lLog << MSG::DEBUG << "TRKR: Tracker from " << calo_dims.rmin() << " to " << calo_dims.rmin()+trkr_tck << endmsg;
+ dd4hep::Volume trkrVolume(trkr_mat, trkrShape, lcdd.material(trkr_mat));
+ dd4hep::PlacedVolume placedTrkrLayer = caloVolume.placeVolume(trkrVolume);
+ trkrLayer.setPlacement(placedTrkrLayer);
+ }
+
+
+ // START OF TONY's STUFF
+ // Calorimeter runs
+ // 1) epitaxial
+ // 2) substrate
+ // 3) passive material (absorber)
+ // 4) air gap between layers (this is the caleVolume made of air
+
+ // STANDARD CONFIGURATION WITH PADDING (air) >= 0 mm
+ // | MODULE 1 | MODULE 2 | ... | MODULE N |
+ // | | | | | | | | | ... | | | | |
+ // | E | S | A | A | E | S | A | A | ... | E | S | A | A |
+ // | P | U | B | I | P | U | B | I | ... | P | U | B | I |
+ // | I | B | S | R | I | B | S | R | ... | I | B | S | R |
+ // | | | | | | | | | ... | | | | |
+
+ // ALTERNATIVE CONFIGURATION WITH PADDING (air) < 0 mm
+ // | MODULE 1 | MODULE 2 | ... | MODULE N |
+ // | | | | | | | | | ... | | | | |
+ // | E | S | A | A | E | S | A | A | ... | E | S | A | A |
+ // | P | U | I | B | P | U | I | B | ... | P | U | I | B |
+ // | I | B | R | S | I | B | R | S | ... | I | B | R | S |
+ // | | | | | | | | | ... | | | | |
+
+ //Leave option for circular geometry if nStaves == 1:
+ if (nStaves == 1)
+ {
+ // loop on the sensitive and substrate layers and place within the caloVolume
+ int l=1;
+ for (int i=0;i
+#include
+
+#include "TFile.h"
+#include "TH1F.h"
+
+/** DigitalCalorimeterSD DetectorDescription/DetSensitive/src/DigitalCalorimeterSD.h DigitalCalorimeterSD.h
+ *
+ * Simple sensitive detector for calorimeter.
+ * It is based on DD4hep::Simulation::Geant4GenericSD (but it is not identical).
+ * In particular, the position of the hit is set to G4Step::GetPreStepPoint() position.
+ * New hit is created for each energy deposit.
+ * No timing information is saved.
+ *
+ * @author Anna Zaborowska
+ */
+
+struct Strixel
+{
+ int layer;
+ int x;
+ int y;
+ int z;
+
+ bool operator<(const Strixel &other) const
+ {
+ if(layer
+ struct hash
+ {
+ std::size_t operator()(const Strixel& s) const
+ {
+ using std::size_t;
+ using std::hash;
+ using std::string;
+
+ // Compute individual hash values for first,
+ // second and third and combine them using XOR
+ // and bit shifting:
+
+ return hash()(
+ std::to_string(s.layer) +","+
+ std::to_string(s.x) +","+
+ std::to_string(s.y) +","+
+ std::to_string(s.z)
+ );
+ }
+ };
+
+}
+
+struct Pad
+{
+ int layer;
+ int r;
+ int z;
+
+ bool operator<(const Pad &other) const
+ {
+ if(layer
+ struct hash
+ {
+ std::size_t operator()(const Pad& p) const
+ {
+ using std::size_t;
+ using std::hash;
+ using std::string;
+
+ // Compute individual hash values for first,
+ // second and third and combine them using XOR
+ // and bit shifting:
+
+ return hash()(
+ std::to_string(p.layer) +","+
+ std::to_string(p.r) +","+
+ std::to_string(p.z)
+ );
+ }
+ };
+
+}
+
+namespace det {
+class DigitalCalorimeterSD : public G4VSensitiveDetector
+{
+ public:
+ /** Constructor.
+ * @param aDetectorName Name of the detector
+ * @param aReadoutName Name of the readout (used to name the collection)
+ * @param aSeg Segmentation of the detector (used to retrieve the cell ID)
+ */
+ DigitalCalorimeterSD(const std::string& aDetectorName,
+ const std::string& aReadoutName,
+ const dd4hep::Segmentation& aSeg);
+ /// Destructor
+ virtual ~DigitalCalorimeterSD();
+ /** Initialization.
+ * Creates the hit collection with the name passed in the constructor.
+ * The hit collection is registered in Geant.
+ * @param aHitsCollections Geant hits collection.
+ */
+ virtual void Initialize(G4HCofThisEvent* aHitsCollections) final;
+ /** Process hit once the particle hit the sensitive volume.
+ * Checks if the energy deposit is larger than 0, calculates the position and cellID,
+ * saves that into the hit collection.
+ * New hit is created for each energy deposit.
+ * @param aStep Step in which particle deposited the energy.
+ */
+ virtual bool ProcessHits(G4Step* aStep, G4TouchableHistory*) final;
+
+ virtual void EndOfEvent(G4HCofThisEvent* aHitsCollections) final;
+
+ virtual bool IsAllowedCellID(unsigned long cid) final;
+ virtual void AddCellIDMask(unsigned long cid) final;
+ virtual void UpdateCellIDMask() final;
+
+private:
+ /// Collection of calorimeter hits
+ G4THitsCollection* m_calorimeterCollection;
+ G4THitsCollection* m_tempCollection;
+ /// Segmentation of the detector used to retrieve the cell Ids
+ dd4hep::Segmentation m_seg;
+
+ std::string m_pixelsOverThresholdFileName;
+ std::ofstream m_pixelsOverThresholdFile;
+
+ std::string m_pixelsOverThresholdPerLayerFileName;
+ std::ofstream m_pixelsOverThresholdPerLayerFile;
+
+ std::string m_mipsPerPixelFileName;
+ std::ofstream m_mipsPerPixelFile;
+
+ std::string m_padMultiplicityFileName;
+ std::ofstream m_padMultiplicityFile;
+
+ bool m_headerPrinted;
+ bool m_countParticlesPerPixel;
+ int m_incidentParticles;
+
+ std::map m_CellIDsMaskedFromPreviousEvent;
+ std::map m_TrackIDsPerEvent;
+
+ //TFile* m_rootFile;
+ //TH1F* m_padMultiplicity;
+};
+}
+
+#endif /* DETSENSITIVE_DIGITALCALORIMETERSD_H */
diff --git a/Detector/DetSensitive/src/DigitalCalorimeterSD.cpp b/Detector/DetSensitive/src/DigitalCalorimeterSD.cpp
new file mode 100644
index 000000000..4787a8aaf
--- /dev/null
+++ b/Detector/DetSensitive/src/DigitalCalorimeterSD.cpp
@@ -0,0 +1,400 @@
+#include "DetSensitive/DigitalCalorimeterSD.h"
+
+// FCCSW
+#include "DDSegmentation/BitField64.h"
+#include "DetCommon/DetUtils.h"
+#include "DDSegmentation/Segmentation.h"
+#include "DDSegmentation/CartesianGridXYZ.h"
+
+// DD4hep
+#include "DDG4/Geant4Mapping.h"
+#include "DDG4/Geant4VolumeManager.h"
+
+
+// CLHEP
+#include "CLHEP/Vector/ThreeVector.h"
+
+// Geant4
+#include "G4SDManager.hh"
+#include "G4ThreeVector.hh"
+#include "G4SystemOfUnits.hh"
+
+#include "iostream"
+#include "map"
+#include "unordered_map"
+
+namespace det {
+DigitalCalorimeterSD::DigitalCalorimeterSD(const std::string& aDetectorName,
+ const std::string& aReadoutName,
+ const dd4hep::Segmentation& aSeg)
+ : G4VSensitiveDetector(aDetectorName), m_seg(aSeg), m_calorimeterCollection(nullptr) {
+ // name of the collection of hits is determined byt the readout name (from XML)
+ collectionName.insert(aReadoutName);
+
+ m_headerPrinted = false;
+ m_countParticlesPerPixel = true;
+
+}
+
+DigitalCalorimeterSD::~DigitalCalorimeterSD(){
+}
+
+void DigitalCalorimeterSD::Initialize(G4HCofThisEvent* aHitsCollections)
+{
+ // create a collection of hits and add it to G4HCofThisEvent
+ // deleted in ~G4Event
+ m_calorimeterCollection = new G4THitsCollection
+ (SensitiveDetectorName,collectionName[0]);
+ aHitsCollections->AddHitsCollection(G4SDManager::GetSDMpointer()->GetCollectionID(m_calorimeterCollection), m_calorimeterCollection);
+
+ //this collection is just a temp and not added to the HitsCollection just used to store raw hits before summing together
+ m_tempCollection = new G4THitsCollection
+ (SensitiveDetectorName,"temp");
+
+ // m_pixelsOverThresholdFile.open("pixelsOverThreshold.txt", std::ios_base::app);.
+
+ if(!m_headerPrinted){
+ if(m_countParticlesPerPixel) {
+ std::cout << "\n\n\n###############################" << std::endl;
+ std::cout << "#\n#" << std::endl;
+ std::cout << "## Counting Pixels in DECAL Mode " << std::endl;
+ std::cout << "#\n#\n###############################" << std::endl;
+ }
+
+ std::cout << "incident particles\tunique cellIDs\tcells over thresh\tunique strixels\tstrixel pixels (3)\tstrixel pixels (7)\tstrixel pixels (15)\tstrixel pixels (31)" << std::endl;
+ m_headerPrinted=true;
+ }
+ m_incidentParticles = 0;
+
+ m_TrackIDsPerEvent.clear();
+
+
+}
+
+bool DigitalCalorimeterSD::ProcessHits(G4Step* aStep, G4TouchableHistory*)
+{
+ // check if energy was deposited
+ G4double edep = aStep->GetTotalEnergyDeposit();
+ if(edep==0.)
+ return false;
+
+ CLHEP::Hep3Vector prePos = aStep->GetPreStepPoint()->GetPosition();
+ CLHEP::Hep3Vector postPos = aStep->GetPostStepPoint()->GetPosition();
+
+
+ bool isIncident = false;
+
+ // count up the number of particles which enter the epitaxial layer
+ if(aStep->GetPreStepPoint()->GetStepStatus() == fGeomBoundary){
+ m_incidentParticles++;
+ isIncident = true;
+
+ G4int trackID = aStep->GetTrack()->GetTrackID();
+ std::map::iterator trackID_it = m_TrackIDsPerEvent.find(trackID);
+ if(m_TrackIDsPerEvent.find(trackID) != m_TrackIDsPerEvent.end()) {
+ (*trackID_it).second++;
+ } else {
+ m_TrackIDsPerEvent[trackID] = 1;
+ }
+ }
+
+ dd4hep::Position pos(prePos.x(), prePos.y(), prePos.z());
+ auto hit = new dd4hep::sim::Geant4CalorimeterHit(pos);
+ //hit->cellID = segmentation::cellID(m_seg, *aStep);
+ hit->cellID = utils::cellID(m_seg, *aStep);
+
+ // if we have a step which is not incident make energy negative
+ // this allows me to count incident particles later on
+ if(!isIncident) {
+ edep *= -1;
+ }
+ hit->energyDeposit = edep;
+ m_tempCollection->insert(hit);
+
+ return true;
+}
+
+
+
+void DigitalCalorimeterSD::EndOfEvent(G4HCofThisEvent*) {
+
+ int nHits = m_tempCollection->entries();
+
+ std::cout << "Total number of TrackIDs in SD = " << m_TrackIDsPerEvent.size() << std::endl;
+
+ // this map stored the energy per cell
+ // map.first = cellID
+ // map.second.first = number of particles into the cell
+ // map.second.second = total energy deposited in the cell
+ std::map > counter;
+ dd4hep::sim::Geant4CalorimeterHit* hit = nullptr;
+
+ //std::cout << "####" << m_seg->type() << std::endl;
+
+ // loop through all of the hits in an event and sum the energy deposited per cell
+ for(int i=0; i (m_tempCollection->GetHit(i));
+
+ // this bit of code allows us to check that is incident and that we should count as new particles
+ // we then make the value of energy +ve again for ease of adding
+ bool isIncident = true;
+ if(hit->energyDeposit < 0.0) {
+ isIncident = false;
+ hit->energyDeposit *= -1;
+ //std::cout << "NON INCIDENT STEP" << std::endl;
+ }
+
+ // ensure that there are no longer energy less than zero
+ if(hit->energyDeposit < 0.0) {
+ std::cout << "STILL A HIT WITH NEGATIVE EDEP MORON!" << std::endl;
+ }
+
+ uint64_t cellID = hit->cellID;
+
+ // for the DEcal chip being designed we have a dead time. This bool ensures the pixel is not dead
+ // if useMask=false then IsAllowedCellID always returns true
+ if(IsAllowedCellID(cellID)){
+
+ std::map >::iterator it = counter.find(cellID);
+
+ if(counter.find(cellID) != counter.end()) {
+
+ // incrimnet the number of incident particles
+ if(isIncident) {
+ (*it).second.first++;
+ }
+ (*it).second.second->energyDeposit += hit->energyDeposit;
+ } else {
+ // this ensures that we only count incident particles
+ // if not incident set to zero so we do count incident later on but do not double count
+ int nparts = 0;
+ if(isIncident) {
+ nparts = 1;
+ }
+ std::pair info = std::make_pair(nparts, hit);
+ counter[cellID] = info;
+ }
+ }
+ }
+
+ // now loop through all the hits and check the Mask. Cannot do it before as may cut out multiple hits to cells
+ bool useMask = false;
+ for(int i=0; i (m_tempCollection->GetHit(i));
+ uint64_t cellID = hit->cellID;
+ AddCellIDMask(cellID);
+ }
+
+ // apply a threshold of 480 electrons to each cell
+ double threshold = 0.0; //480*0.0000036; // 480 electrons in MeV which should be default value
+ int pixelsOverThreshold = 0;
+ int pixelsOverThresholdPerLayer[50];
+ std::vector mipsPerPixel[50];
+ for(int n(0); n<50; ++n){
+ pixelsOverThresholdPerLayer[n]=0;
+ mipsPerPixel[n].clear();
+ }
+
+ dd4hep::sim::Geant4CalorimeterHit* digi = nullptr;
+ std::map >::iterator it;
+ std::map StrixelHitCounter;
+ for(it=counter.begin(); it!=counter.end(); it++)
+ {
+ if(useMask) {
+ AddCellIDMask((*it).second.second->cellID);
+ }
+
+ // if the energy deposited in the cellID is gtreather than threshold then do something
+ if((*it).second.second->energyDeposit>=threshold) {
+ pixelsOverThreshold++;
+ // insert the collection to G4 so can be passed on later
+ digi = dynamic_cast ((*it).second.second);
+ // this line overwrites the time parameter with the number of particles per pixel
+ std::vector conts;
+ dd4hep::sim::Geant4Hit::MonteCarloContrib cont(0,0,0,0);
+
+ if(m_countParticlesPerPixel) {
+ //digi->energyDeposit = (*it).second.first; // overwrite the energy deposit with the number of mips in the cell
+ cont.time = (*it).second.first;
+ }
+
+ conts.push_back(cont);
+ digi->truth = conts;
+ m_calorimeterCollection->insert(digi); // this passes to ROOT the readout... interesting!
+
+ }
+ }
+
+ std::cout << m_TrackIDsPerEvent.size() << " " << m_tempCollection->entries() << " " << m_incidentParticles << " " << m_calorimeterCollection->entries() << std::endl;
+
+ UpdateCellIDMask();
+
+ // below this point is DECal chip specific where the readout is reconfidured
+ // not needed for the analogue SiW
+
+ int strixelHits = StrixelHitCounter.size();
+ int strixelPixelHits_3 = 0;
+ int strixelPixelHits_7 = 0;
+ int strixelPixelHits_15 = 0;
+ int strixelPixelHits_31 = 0;
+ //int strixelMaxPixels = 3;
+
+ std::map > PadHitCounter;
+ std::map::iterator it_strixel;
+ for(it_strixel=StrixelHitCounter.begin(); it_strixel!=StrixelHitCounter.end(); ++it_strixel) {
+
+ // convert strixel to pad
+ double x1 = (*it_strixel).first.x;
+ double y1 = (*it_strixel).first.y;
+ double r1 = pow((x1*x1+y1*y1),0.5);
+
+ double x0 = r1;
+ double y0 = 0.0; // set r0 to be at y=0 x=radius
+ double cord = pow(pow((x1-x0),2)+pow((y1-y0),2),0.5);
+ double theta = 2*asin(cord/(2*r1));
+ double pos_round_circ = theta*r1;
+
+ // for pad counting
+ std::vector counter;
+ counter.resize(6);
+ counter.at(0) = 1;
+
+ int npixels = (*it_strixel).second;
+
+ if(npixels<3) {
+ strixelPixelHits_3 += npixels;
+ counter.at(1) = npixels;
+ }
+ else {
+ strixelPixelHits_3 += 3;
+ counter.at(1) = 3;
+ }
+
+ if(npixels<7) {
+ strixelPixelHits_7 += npixels;
+ counter.at(2) = npixels;
+ }
+ else {
+ strixelPixelHits_7 += 7;
+ counter.at(2) = 7;
+ }
+ if(npixels<15) {
+ strixelPixelHits_15 += npixels;
+ counter.at(3) = npixels;
+ }
+ else {
+ strixelPixelHits_15 += 15;
+ counter.at(3) = 15;
+ }
+
+ if(npixels<31) {
+ strixelPixelHits_31 += npixels;
+ counter.at(4) = npixels;
+ }
+ else {
+ strixelPixelHits_31 += 31;
+ counter.at(4) = 31;
+ }
+
+ counter.at(5) = npixels;
+
+ //-------------------------------------------------
+ Pad pad;
+ pad.layer = (*it_strixel).first.layer;
+ pad.r = floor(pos_round_circ/5.0);
+ pad.z = (*it_strixel).first.z;
+
+ std::map >::iterator it_pad = PadHitCounter.find(pad);
+ if( PadHitCounter.find(pad) != PadHitCounter.end() ) {
+ // add to the relevant int part
+ (*it_pad).second.at(0) += counter.at(0); // count hits in pad limited by 1 hit per column
+ (*it_pad).second.at(1) += counter.at(1); // count hits in pad limited by 3 hits per column
+ (*it_pad).second.at(2) += counter.at(2); // count hits in pad limited by 7 hits per column
+ (*it_pad).second.at(3) += counter.at(3); // count hits in pad limited by 15 hits per column
+ (*it_pad).second.at(4) += counter.at(4); // count hits in pad limited by 31 hits per column
+ (*it_pad).second.at(5) += counter.at(5); // count all hits in pad
+ } else {
+ PadHitCounter[pad] = counter;
+ }
+
+ }
+
+ // now we loop through the map of pads to get total hits
+ if(!m_padMultiplicityFile.is_open()) {
+ std::map >::iterator it_pad;
+ m_padMultiplicityFile.open(m_padMultiplicityFileName, std::ios_base::app);
+ for(it_pad=PadHitCounter.begin(); it_pad!=PadHitCounter.end(); ++it_pad) {
+
+ //std::cout << (*it_pad).second.at(1) << " ";
+ m_padMultiplicityFile << (*it_pad).second.at(1) << " ";
+
+ }
+ m_padMultiplicityFile << "\n";
+ m_padMultiplicityFile.close();
+ }
+
+// std::cout << std::endl;
+// std::cout << nHits << " " << pixelsOverThreshold << std::endl;
+// std::cout << "Total number of unique cellIDs = " << pixelsOverThreshold << std::endl;
+// std::cout << "Total number of unique strixels = " << strixelHits << std::endl;
+// std::cout << "Total number of strixel pixels (3) = " << strixelPixelHits_3 << std::endl ;
+// std::cout << "Total number of strixel pixels (7) = " << strixelPixelHits_7 << std::endl;
+// std::cout << "Total number of strixel pixels (15) = " << strixelPixelHits_15 << std::endl;
+// std::cout << "Total number of strixel pixels (31) = " << strixelPixelHits_31 << std::endl;
+// std::cout << "Total number of unique pads = " << PadHitCounter.size() << std::endl << std::endl;
+
+
+ std::cout << "Counters\t"<< m_incidentParticles << "\t" << counter.size() << "\t" << pixelsOverThreshold << "\t" << strixelHits << "\t" << strixelPixelHits_3 << "\t" << strixelPixelHits_7 << "\t" << strixelPixelHits_15 << "\t" << strixelPixelHits_31 << "\t" << PadHitCounter.size() << std::endl;
+
+
+
+
+
+}
+
+
+bool DigitalCalorimeterSD::IsAllowedCellID(unsigned long cid) {
+ // check if the cellID occurred in the previous event
+ // if the counter is zero the cell is allowed and can remove cell from the mask
+ std::map::iterator valid_it = m_CellIDsMaskedFromPreviousEvent.find(cid);
+ if(m_CellIDsMaskedFromPreviousEvent.find(cid) != m_CellIDsMaskedFromPreviousEvent.end()) {
+ if((*valid_it).second > 0) {
+ //std::cout << "CellID " << cid << " is inactive for the next " << (*valid_it).second << " events" << std::endl;
+ return false;
+ } else {
+ m_CellIDsMaskedFromPreviousEvent.erase(valid_it);
+ return true;
+ }
+ } else {
+ return true;
+ }
+}
+
+void DigitalCalorimeterSD::AddCellIDMask(unsigned long cid) {
+ // look for the cellID in the std::map. If it exists then reset the dead time
+ // if it does not exist then add it to the map
+ int nEventsMasked = 40;
+ std::map::iterator add_it = m_CellIDsMaskedFromPreviousEvent.find(cid);
+ if(m_CellIDsMaskedFromPreviousEvent.find(cid) != m_CellIDsMaskedFromPreviousEvent.end()) {
+ (*add_it).second = nEventsMasked;
+ } else {
+ m_CellIDsMaskedFromPreviousEvent[cid] = nEventsMasked;
+ }
+}
+
+void DigitalCalorimeterSD::UpdateCellIDMask() {
+ // loop through the mask and remove one from each mask ID.
+ // if it reaches zero then remove from the list.
+ std::map::iterator update_it;
+ for (update_it = m_CellIDsMaskedFromPreviousEvent.begin(); update_it != m_CellIDsMaskedFromPreviousEvent.end(); update_it++){
+ (*update_it).second -= 1;
+ }
+ for (update_it = m_CellIDsMaskedFromPreviousEvent.begin(); update_it != m_CellIDsMaskedFromPreviousEvent.end(); update_it++){
+ if((*update_it).second == 0) m_CellIDsMaskedFromPreviousEvent.erase(update_it);
+ }
+ std::cout << "There are currently " << m_CellIDsMaskedFromPreviousEvent.size() << " pixels masked out" << std::endl;
+}
+
+
+}
+
diff --git a/Detector/DetSensitive/src/SDWrapper.cpp b/Detector/DetSensitive/src/SDWrapper.cpp
index 2a76e12fc..53da0c58e 100644
--- a/Detector/DetSensitive/src/SDWrapper.cpp
+++ b/Detector/DetSensitive/src/SDWrapper.cpp
@@ -5,6 +5,7 @@
#include "DetSensitive/BirksLawCalorimeterSD.h"
#include "DetSensitive/FullParticleAbsorptionSD.h"
#include "DetSensitive/GflashCalorimeterSD.h"
+#include "DetSensitive/DigitalCalorimeterSD.h"
#include "DetSensitive/MiddleStepTrackerSD.h"
#include "DetSensitive/SimpleCalorimeterSD.h"
#include "DetSensitive/SimpleTrackerSD.h"
@@ -47,6 +48,13 @@ static G4VSensitiveDetector* create_aggregate_calorimeter_sd(const std::string&
return new det::AggregateCalorimeterSD(
aDetectorName, readoutName, aLcdd.sensitiveDetector(aDetectorName).readout().segmentation());
}
+// Factory method to create an instance of DigitalCalorimeterSD
+static G4VSensitiveDetector* create_digital_calorimeter_sd(const std::string& aDetectorName,
+ dd4hep::Detector& aLcdd) {
+ std::string readoutName = aLcdd.sensitiveDetector(aDetectorName).readout().name();
+ return new det::DigitalCalorimeterSD(
+ aDetectorName,readoutName,aLcdd.sensitiveDetector(aDetectorName).readout().segmentation());
+}
// Factory method to create an instance of GflashCalorimeterSD
static G4VSensitiveDetector* create_gflash_calorimeter_sd(const std::string& aDetectorName,
dd4hep::Detector& aLcdd) {
@@ -70,3 +78,4 @@ DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(BirksLawCalorimeterSD, dd4hep::sim::cre
DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(AggregateCalorimeterSD, dd4hep::sim::create_aggregate_calorimeter_sd)
DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(GflashCalorimeterSD, dd4hep::sim::create_gflash_calorimeter_sd)
DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(FullParticleAbsorptionSD, dd4hep::sim::create_full_particle_absorbtion_sd)
+DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(DigitalCalorimeterSD, dd4hep::sim::create_digital_calorimeter_sd)
diff --git a/Detector/DetStudies/src/components/DECalAnalysis.cpp b/Detector/DetStudies/src/components/DECalAnalysis.cpp
new file mode 100644
index 000000000..334e21c18
--- /dev/null
+++ b/Detector/DetStudies/src/components/DECalAnalysis.cpp
@@ -0,0 +1,295 @@
+
+#include "DECalAnalysis.h"
+
+// FCCSW
+#include "DetInterface/IGeoSvc.h"
+
+// datamodel
+#include "datamodel/PositionedCaloHitCollection.h"
+#include "datamodel/CaloHitCollection.h"
+#include "datamodel/MCParticleCollection.h"
+
+#include "CLHEP/Vector/ThreeVector.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "TH1F.h"
+#include "TVector2.h"
+#include "TMath.h"
+#include "TString.h"
+#include "string"
+
+// DD4hep
+#include "DD4hep/Detector.h"
+#include "DD4hep/Readout.h"
+
+DECLARE_ALGORITHM_FACTORY(DECalAnalysis)
+
+DECalAnalysis::DECalAnalysis(const std::string& aName, ISvcLocator* aSvcLoc)
+ : GaudiAlgorithm(aName, aSvcLoc),
+ m_histSvc("THistSvc", "DECalAnalysis"),
+ m_geoSvc("GeoSvc", "DECalAnalysis")
+{
+ declareProperty("pixels", m_deposits, "Energy deposits in sampling calorimeter (input)");
+ declareProperty("pads", m_pads, "Energy deposits in Pad mode (input)");
+ declareProperty("truth", m_truth, "Generated particle truth");
+
+ m_calP0 = 239.74;
+ m_calP1 = 99.1602;
+ m_calP2 = -0.0143;
+
+}
+DECalAnalysis::~DECalAnalysis() {}
+
+StatusCode DECalAnalysis::initialize() {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
+ }
+ // check if readouts exist
+ if (m_geoSvc->lcdd()->readouts().find(m_pixelReadoutName) == m_geoSvc->lcdd()->readouts().end()) {
+ error() << "Readout <<" << m_pixelReadoutName << ">> does not exist." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_sumPixelsLayers.assign(m_numLayers+1,0);
+ m_sumPadsLayers.assign(m_numLayers+1,0);
+ m_sumParticlesLayers.assign(m_numLayers+1,0);
+
+ m_treeInfo = new TTree("decal_info", "for Kostas");
+ m_treeInfo->Branch("event_number", &m_treeEvtNumber);
+ m_treeInfo->Branch("truth_energy", &m_truthEnergy);
+ m_treeInfo->Branch("edep_tot", &m_sumEnergyDep);
+ m_treeInfo->Branch("pixels_cal", &m_calPixels);
+ m_treeInfo->Branch("pixels_tot", &m_sumPixels);
+ m_treeInfo->Branch("particles_tot", &m_sumParticles);
+ for(uint l(0); lBranch(Form("pixels_l%i", l), &m_sumPixelsLayers[l]);
+ }
+ m_treeInfo->Branch("pads_tot", &m_sumPads);
+ for(uint l(0); lBranch(Form("pads_l%i", l), &m_sumPadsLayers[l]);
+ }
+ for(uint l(0); lBranch(Form("particles_l%i", l), &m_sumParticlesLayers[l]);
+ }
+ if(m_histSvc->regTree("/rec/decalInfo", m_treeInfo).isFailure()) {
+ error() << "Couldn't register tree" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ // create histograms
+ m_particlesPerPixel = new TProfile("decal_ParticlesPerPixel", "Particles Per Pixels; Layer Number; Particles / Pixel", m_numLayers+1, -0.5, m_numLayers+0.5);
+ if (m_histSvc->regHist("/rec/profile_test", m_particlesPerPixel).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+
+ m_totalHits = new TH1F("decal_totalHits", "Total hits in an event", 15000, 0, 150000);
+ if (m_histSvc->regHist("/rec/decal_hits", m_totalHits).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ m_totalHitsVsLayer = new TH1F("decal_totalHitsVsLayer", "Total hits per layer", m_numLayers+1, -0.5, m_numLayers+0.5);
+ m_totalHitsVsLayer->Sumw2(kFALSE);
+ if(m_histSvc->regHist("/rec/decal_hitsVsLayer", m_totalHitsVsLayer).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_totalPadsVsLayer = new TH1F("decal_totalPadsVsLayer", "Total pads per layer", m_numLayers+1, -0.5, m_numLayers+0.5);
+ m_totalPadsVsLayer->Sumw2(kFALSE);
+ if(m_histSvc->regHist("/rec/decal_padsVsLayer", m_totalPadsVsLayer).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ for(uint n(0); n<100; n++) {
+ m_hitsVsLayerEvent[n] = new TH1F(("decal_hitsVsLayerEvent"+std::to_string(n)).c_str(), "", m_numLayers+1, -0.5, m_numLayers+0.5);
+ m_hitsVsLayerEvent[n]->Sumw2(kFALSE);
+ if(m_histSvc->regHist("/rec/decal_hitsVsLayerEvent"+std::to_string(n), m_hitsVsLayerEvent[n]).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+
+ m_hitsInMaxLayer = new TH1F("decal_hitsInMaxLayer", "Total hits in max layer", 150, 0, 15000);
+ if (m_histSvc->regHist("/rec/decal_hitsInMaxLayer", m_hitsInMaxLayer).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_maxLayer = new TH1F("decal_maxLayer", "Layer with the maximum hits in an event",m_numLayers+1, -0.5, m_numLayers+0.5 );
+ if (m_histSvc->regHist("/rec/decal_maxLayer", m_maxLayer).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_totalPadsAbovePixelThreshold = new TH1F("decal_totalPadsAbovePixelThreshold","", 1000,0,10000);
+ if (m_histSvc->regHist("/rec/decal_totalPadsAbovePixelThreshold", m_totalPadsAbovePixelThreshold).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_totalPads = new TH1F("decal_totalPads","", 1000,0,10000);
+ if (m_histSvc->regHist("/rec/decal_totalPads", m_totalPads).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ for (uint i = 0; i < m_numLayers+1; i++) {
+ m_totalHitsLayers.push_back(new TH1F(("decal_Hits_layer" + std::to_string(i)).c_str(),
+ ("Total pixel hits in layer " + std::to_string(i)).c_str(), 1000, 0,
+ 10000));
+ if (m_histSvc->regHist("/rec/decal_total_hits_layer" + std::to_string(i), m_totalHitsLayers.back()).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_totalPadsLayers.push_back(new TH1F(("decal_Pads_layer" + std::to_string(i)).c_str(),
+ ("Total pads with hits in layer " + std::to_string(i)).c_str(), 100, 0,
+ 1000));
+ if (m_histSvc->regHist("/rec/decal_total_pads_layer" + std::to_string(i), m_totalPadsLayers.back()).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+
+ m_totalPixelsPerPadLayers.push_back(new TH1F(("decal_Pixels_Per_Pad_layer" + std::to_string(i)).c_str(),
+ ("Total pixels per pad in layer " + std::to_string(i)).c_str(), 200, 0,
+ 200));
+ if (m_histSvc->regHist("/rec/decal_Pixels_Per_Pad_layer" + std::to_string(i), m_totalPixelsPerPadLayers.back()).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ m_XYEvent0 = new TH2F("decal_XYEvent0", "XY Positions for Event 0", 5000, -2500, 2500, 5000, -2500, 2500);
+ if (m_histSvc->regHist("/rec/decal_XYEvent0", m_XYEvent0).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_EtaPhiEvent0 = new TH2F("decal_EtaPhiEvent0", "Eta Phi Positions for Event 0", 222, -0.004, 0.002, 222, -0.8305,-0.8245);
+ if (m_histSvc->regHist("/rec/decal_EtaPhiEvent0", m_EtaPhiEvent0).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_EtaPhiSeparation = new TH1F("decal_EtaSeparation", "Eta Phi Positions", 1200, 0, 1);
+ if (m_histSvc->regHist("/rec/decal_EtaPhiSeparation", m_EtaPhiSeparation).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_nEventsFilled = 0;
+
+ m_calFit = new TF1("calFit", "pol2", 0,3000);
+ m_calFit->SetParameter("p0", m_calP0);
+ m_calFit->SetParameter("p1", m_calP1);
+ m_calFit->SetParameter("p2", m_calP2);
+
+ // if it gets to here then it has registered all histograms
+ return StatusCode::SUCCESS;
+}
+
+StatusCode DECalAnalysis::execute() {
+ // get the decoder for the readout
+ // this will allow to extract layers etc.
+ auto decoder = m_geoSvc->lcdd()->readout(m_pixelReadoutName).idSpec().decoder();
+
+ // calculate the MC truth energy of the incident particle
+ const auto truth = m_truth.get();
+ for (const auto& t : *truth) {
+ double px = t.core().p4.px;
+ double py = t.core().p4.py;
+ double pz = t.core().p4.pz;
+ double m = t.core().p4.mass;
+
+ m_truthEnergy = TMath::Sqrt(m*m + px*px + py*py + pz*pz);
+ }
+
+ // added a little print out to ensure it is still working
+ if(m_treeEvtNumber < 10 or m_treeEvtNumber % 100 == 0) {
+ info() << "Event Number " << m_treeEvtNumber << endmsg;
+ info() << "Truth energy of single particle = " << m_truthEnergy << endmsg;
+ }
+
+
+ // set some variables to hold information for an event
+ m_sumPixels = 0.;
+ m_sumParticles = 0.;
+ m_sumPixelsLayers.assign(m_numLayers+1, 0);
+ m_sumParticlesLayers.assign(m_numLayers+1, 0);
+ m_meanEta.assign(m_numLayers+1, 0);
+ m_meanPhi.assign(m_numLayers+1, 0);
+
+ m_sumEnergyDep = 0;
+
+ bool fillXY = true;
+ if (m_XYEvent0->GetEntries()>0) {
+ fillXY = false;
+ }
+
+ int isdigital = 0, notdigital = 0;
+
+ const auto deposits = m_deposits.get();
+ for (const auto& hit : *deposits) {
+ decoder->setValue(hit.core().cellId);
+
+ int layer = (*decoder)[m_layerFieldName];
+ m_sumPixelsLayers[layer] += 1;
+
+ // check if energy was deposited in the calorimeter (active/passive material)
+ // layers are numbered starting from 1, layer == 0 is cryostat/bath
+ if (layer > 0) {
+ m_sumPixels++;
+ decoder->setValue(hit.core().cellId);
+ int layer = (*decoder)[m_layerFieldName];
+ // get the energy deposited in the active area
+ double edep = hit.core().energy;
+ m_sumEnergyDep += edep;
+ double particles = hit.core().time; // I overwrite now the time for number of particles
+ m_sumParticles += particles;
+ m_sumParticlesLayers[layer] += particles; // call this particles for DEcal as changed SD readout to pass nparticles / pixel not enerrgy dep / pixel
+ m_particlesPerPixel->Fill(layer, particles);
+ // calculate eta and phi
+ double r = TMath::Sqrt(hit.position().x*hit.position().x + hit.position().y*hit.position().y + hit.position().z*hit.position().z);
+ double theta = TMath::ACos(hit.position().z/r);
+ double eta = -TMath::Log2(TMath::Tan(theta/2.0));
+ double phi = TMath::ATan(hit.position().y / hit.position().x);
+
+ m_meanEta[layer] += eta;
+ m_meanPhi[layer] += phi;
+
+ }
+ }
+
+ std::cout << "isdigital = " << isdigital << std::endl;
+ std::cout << "notdigital = " << notdigital << std::endl;
+
+ if(m_sumPixels>0 || m_sumPads>0) {
+ m_treeInfo->Fill();
+ }
+
+ m_treeEvtNumber++;
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode DECalAnalysis::finalize() {
+ debug() << "StatusCode DECalAnalysis::finalize()" << endmsg;
+
+ // count from 1 to avoid the hits in the Trkr volume
+ for (uint i = 1; i < m_numLayers+1; i++) {
+ m_totalHitsVsLayer->SetBinContent(i+1,m_totalHitsLayers[i]->GetMean());
+ m_totalPadsVsLayer->SetBinContent(i+1,m_totalPadsLayers[i]->GetMean());
+ }
+
+return GaudiAlgorithm::finalize(); }
+
+Double_t DECalAnalysis::FitLongProfile( Double_t *x, Double_t *par)
+{
+ double E_0 = par[0];
+ double a = par[1];
+ double b = par[2];
+ double led = E_0*b * (pow( (b*x[0]), a-1 )* exp(-(b*x[0])) ) / TMath::Gamma(a);
+
+ return led;
+}
diff --git a/Detector/DetStudies/src/components/DECalAnalysis.h b/Detector/DetStudies/src/components/DECalAnalysis.h
new file mode 100644
index 000000000..9c1ca574f
--- /dev/null
+++ b/Detector/DetStudies/src/components/DECalAnalysis.h
@@ -0,0 +1,140 @@
+#ifndef DETSTUDIES_DECalAnalysis_H
+#define DETSTUDIES_DECalAnalysis_H
+
+// GAUDI
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// FCCSW
+#include "FWCore/DataHandle.h"
+#include "TH2F.h"
+#include "TTree.h"
+#include "TF1.h"
+#include "TProfile.h"
+
+class IGeoSvc;
+
+// datamodel
+namespace fcc {
+class PositionedCaloHitCollection;
+class CaloHitCollection;
+class MCParticleCollection;
+}
+
+class TH1F;
+class ITHistSvc;
+/** @class DECalAnalysis DECalAnalysis.h
+ *
+ * Histograms of energy deposited in active material and total energy deposited in the calorimeter.
+ * 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 DECalAnalysis : public GaudiAlgorithm {
+public:
+ explicit DECalAnalysis(const std::string&, ISvcLocator*);
+ virtual ~DECalAnalysis();
+ /** 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;
+
+ Double_t FitLongProfile( Double_t *x, Double_t *par);
+
+private:
+ /// Pointer to the interface of histogram service
+ ServiceHandle m_histSvc;
+ /// Pointer to the geometry service
+ ServiceHandle m_geoSvc;
+ /// Number of layers/cells
+ Gaudi::Property m_numLayers{this, "numLayers", 50, "Number of layers"};
+ /// Name of the layer/cell field
+ Gaudi::Property m_layerFieldName{this, "layerFieldName", "", "Identifier of layers"};
+
+ /// The Pixel analysis
+ /// Handle for the energy deposits
+ DataHandle m_deposits{"rec/pixels", Gaudi::DataHandle::Reader, this};
+ /// Name of the detector readout
+ Gaudi::Property m_pixelReadoutName{this, "pixelReadoutName", "", "Name of the detector readout"};
+
+ /// The Pad analysis
+ /// Handle for the energy deposits
+ DataHandle m_pads{"rec/pads", Gaudi::DataHandle::Reader, this};
+ /// Name of the detector readout
+ Gaudi::Property m_padReadoutName{this, "padReadoutName", "", "Name of the detector readout"};
+
+ DataHandle m_truth{"truth", Gaudi::DataHandle::Reader, this};
+
+
+ // Histograms of total deposited energy within layer
+ // Histograms for the number of hits per layer
+ std::vector m_totalHitsLayers;
+
+ // Histogram of total hits in the calorimeter
+ TH1F* m_totalHits;
+ // Histogram of number of hits per layer in single histogram
+ TH1F* m_totalHitsVsLayer;
+ TH1F* m_hitsVsLayerEvent[100];
+
+ TH1F* m_totalPadsVsLayer;
+ std::vector m_totalPadsLayers;
+ std::vector m_totalPixelsPerPadLayers;
+ TH1F* m_totalPixelsPerPad;
+ TH1F* m_totalPads;
+ TH1F* m_totalPadsAbovePixelThreshold;
+
+ TH2F* m_XYEvent0;
+ TH2F* m_EtaPhiEvent0;
+ TH1F* m_EtaPhiSeparation;
+
+ TProfile* m_particlesPerPixel;
+
+ // Histogram of the number of hits in the layer with most hits
+ TH1F* m_hitsInMaxLayer;
+ // Histogram of the layer with maximum hits
+ TH1F* m_maxLayer;
+
+ // calibration factors to compensate non linearity
+ float m_calP0, m_calP1, m_calP2;
+ TF1* m_calFit;
+
+ TTree* m_treeInfo;
+ int m_treeEvtNumber;
+
+ int m_nEventsFilled;
+ int m_truthEnergy;
+ float m_calPixels;
+ int m_sumPixels;
+ std::vector m_sumPixelsLayers;
+ std::vector m_sumParticlesLayers;
+ int m_sumPads;
+ std::vector m_sumPadsLayers;
+
+ float m_sumEnergyDep;
+ float m_sumParticles;
+
+ std::vector m_meanEta;
+ std::vector m_meanPhi;
+
+ double m_fitE0;
+ double m_fita;
+ double m_fitb;
+ double m_rising50;
+ double m_rising90;
+ double m_fwhm;
+
+ TF1* m_fitLongProfile;
+};
+#endif /* DETSTUDIES_DECalAnalysis_H */
diff --git a/Detector/DetStudies/src/components/DECalLongitudinalTest.cpp b/Detector/DetStudies/src/components/DECalLongitudinalTest.cpp
new file mode 100644
index 000000000..6f73960d5
--- /dev/null
+++ b/Detector/DetStudies/src/components/DECalLongitudinalTest.cpp
@@ -0,0 +1,142 @@
+
+#include "DECalLongitudinalTest.h"
+
+// FCCSW
+#include "DetInterface/IGeoSvc.h"
+
+// datamodel
+#include "datamodel/PositionedCaloHitCollection.h"
+#include "datamodel/CaloHitCollection.h"
+
+#include "CLHEP/Vector/ThreeVector.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "TH1F.h"
+#include "TVector2.h"
+#include "TMath.h"
+
+// DD4hep
+#include "DD4hep/Detector.h"
+#include "DD4hep/Readout.h"
+
+DECLARE_ALGORITHM_FACTORY(DECalLongitudinalTest)
+
+DECalLongitudinalTest::DECalLongitudinalTest(const std::string& aName, ISvcLocator* aSvcLoc)
+ : GaudiAlgorithm(aName, aSvcLoc),
+ m_histSvc("THistSvc", "DECalLongitudinalTest"),
+ m_geoSvc("GeoSvc", "DECalLongitudinalTest")
+{
+ declareProperty("deposits", m_deposits, "Energy deposits in sampling calorimeter (input)");
+}
+DECalLongitudinalTest::~DECalLongitudinalTest() {}
+
+StatusCode DECalLongitudinalTest::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
+ m_totalHits = new TH1F("decal_totalHits", "Total hits in an event", 15000, 0, 150000);
+ if (m_histSvc->regHist("/rec/decal_hits", m_totalHits).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ m_totalHitsVsLayer = new TH1F("decal_totalHitsVsLayer", "Total hits per layer", m_numLayers+1, -0.5, m_numLayers+0.5);
+ m_totalHitsVsLayer->Sumw2(kFALSE);
+ if(m_histSvc->regHist("/rec/decal_hitsVsLayer", m_totalHitsVsLayer).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ for (uint i = 0; i < m_numLayers; i++) {
+ m_totalHitsLayers.push_back(new TH1F(("decal_Hits_layer" + std::to_string(i)).c_str(),
+ ("Total pixel hits in layer " + std::to_string(i)).c_str(), 1000, 0,
+ 10000));
+ if (m_histSvc->regHist("/rec/decal_total_hits_layer" + std::to_string(i), m_totalHitsLayers.back()).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ m_XYEvent0 = new TH2F("decal_XYEvent0", "XY Positions for Event 0", 5000, -2500, 2500, 5000, -2500, 2500);
+ if (m_histSvc->regHist("/rec/decal_XYEvent0", m_XYEvent0).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_EtaPhiEvent0 = new TH2F("decal_EtaPhiEvent0", "Eta Phi Positions for Event 0", 2000, -0.04, 0.02, 2000, -4,4);
+ if (m_histSvc->regHist("/rec/decal_EtaPhiEvent0", m_EtaPhiEvent0).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_EtaPhiSeparation = new TH1F("decal_EtaSeparation", "Eta Phi Positions", 1200, 0, 1);
+ if (m_histSvc->regHist("/rec/decal_EtaPhiSeparation", m_EtaPhiSeparation).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+
+
+ // if it gets to here then it has registered all histograms
+ return StatusCode::SUCCESS;
+}
+
+StatusCode DECalLongitudinalTest::execute() {
+ // get the decoder for the readout
+ // this will allow to extract layers etc.
+ auto decoder = m_geoSvc->lcdd()->readout(m_readoutName).idSpec().decoder();
+
+ // set some variables to hold information for an event
+ double sumHits = 0.;
+ std::vector sumHitslayers;
+ sumHitslayers.assign(m_numLayers, 0);
+
+ bool fillXY = true;
+ if (m_XYEvent0->GetEntries()>0) {
+ fillXY = false;
+ }
+
+ const auto deposits = m_deposits.get();
+ for (const auto& hit : *deposits) {
+ decoder->setValue(hit.core().cellId);
+ sumHitslayers[(*decoder)[m_layerFieldName]] += 1;
+ // check if energy was deposited in the calorimeter (active/passive material)
+ // layers are numbered starting from 1, layer == 0 is cryostat/bath
+ std::cout << (*decoder)[m_layerFieldName] << std::endl;
+ if ((*decoder)[m_layerFieldName] > 0) {
+ sumHits++;
+ decoder->setValue(hit.core().cellId);
+ // calculate eta and phi
+ double r = TMath::Sqrt(hit.position().x*hit.position().x + hit.position().y*hit.position().y + hit.position().z*hit.position().z);
+ double theta = TMath::ACos(hit.position().z/r);
+ double eta = -TMath::Log2(TMath::Tan(theta/2.0));
+ double phi = TMath::ATan(hit.position().y / hit.position().x);
+ if (fillXY) {
+ m_XYEvent0->Fill((*decoder)["x"],(*decoder)["y"]);
+ m_XYEvent0->Fill(hit.position().x, hit.position().y);
+ if((*decoder)[m_layerFieldName] == 18) {
+ m_EtaPhiEvent0->Fill(eta,phi);
+ }
+ }
+ }
+ }
+
+ // Fill histograms
+ m_totalHits->Fill(sumHits);
+ for (uint i = 1; i < m_numLayers; i++) {
+ m_totalHitsLayers[i]->Fill(sumHitslayers[i]);
+ }
+ return StatusCode::SUCCESS;
+}
+
+StatusCode DECalLongitudinalTest::finalize() {
+ debug() << "StatusCode DECalLongitudinalTest::finalize()" << endmsg;
+
+ // count from 1 to avoid the hits in the Trkr volume
+ for (uint i = 1; i < m_numLayers; i++) {
+ m_totalHitsVsLayer->SetBinContent(i+1,m_totalHitsLayers[i]->GetMean());
+ }
+
+return GaudiAlgorithm::finalize(); }
diff --git a/Detector/DetStudies/src/components/DECalLongitudinalTest.h b/Detector/DetStudies/src/components/DECalLongitudinalTest.h
new file mode 100644
index 000000000..ac2f5ff36
--- /dev/null
+++ b/Detector/DetStudies/src/components/DECalLongitudinalTest.h
@@ -0,0 +1,75 @@
+#ifndef DETSTUDIES_DECALLONGITUDINALTEST_H
+#define DETSTUDIES_DECALLONGITUDINALTEST_H
+
+// GAUDI
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// FCCSW
+#include "FWCore/DataHandle.h"
+#include "TH2F.h"
+class IGeoSvc;
+
+// datamodel
+namespace fcc {
+class PositionedCaloHitCollection;
+class CaloHitCollection;
+}
+
+class TH1F;
+class ITHistSvc;
+/** @class DECalLongitudinalTest DECalLongitudinalTest.h
+ *
+ * Histograms of energy deposited in active material and total energy deposited in the calorimeter.
+ * 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 DECalLongitudinalTest : public GaudiAlgorithm {
+public:
+ explicit DECalLongitudinalTest(const std::string&, ISvcLocator*);
+ virtual ~DECalLongitudinalTest();
+ /** 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 m_histSvc;
+ /// Pointer to the geometry service
+ ServiceHandle m_geoSvc;
+ /// Handle for the energy deposits
+ DataHandle m_deposits{"rec/caloHits", Gaudi::DataHandle::Reader, this};
+ /// Name of the layer/cell field
+ Gaudi::Property m_layerFieldName{this, "layerFieldName", "", "Identifier of layers"};
+ /// Number of layers/cells
+ Gaudi::Property m_numLayers{this, "numLayers", 8, "Number of layers"};
+ /// Name of the detector readout
+ Gaudi::Property m_readoutName{this, "readoutName", "", "Name of the detector readout"};
+
+ // Histograms of total deposited energy within layer
+ // Histograms for the number of hits per layer
+ std::vector m_totalHitsLayers;
+ // Histogram of total hits in the calorimeter
+ TH1F* m_totalHits;
+ // Histogram of number of hits per layer in single histogram
+ TH1F* m_totalHitsVsLayer;
+
+ TH2F* m_XYEvent0;
+ TH2F* m_EtaPhiEvent0;
+ TH1F* m_EtaPhiSeparation;
+};
+#endif /* DETSTUDIES_DECALLONGITUDINALTEST_H */
diff --git a/Detector/DetStudies/src/components/FilterSiliconEcalHits.cpp b/Detector/DetStudies/src/components/FilterSiliconEcalHits.cpp
new file mode 100644
index 000000000..dc2d4dc8e
--- /dev/null
+++ b/Detector/DetStudies/src/components/FilterSiliconEcalHits.cpp
@@ -0,0 +1,107 @@
+
+#include "FilterSiliconEcalHits.h"
+
+// FCCSW
+#include "DetInterface/IGeoSvc.h"
+
+// datamodel
+#include "datamodel/PositionedCaloHitCollection.h"
+#include "datamodel/CaloHitCollection.h"
+#include "datamodel/MCParticleCollection.h"
+
+#include "CLHEP/Vector/ThreeVector.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "TH1F.h"
+#include "TVector2.h"
+#include "TMath.h"
+#include "TString.h"
+#include "string"
+
+// DD4hep
+#include "DD4hep/Detector.h"
+#include "DD4hep/Readout.h"
+
+DECLARE_ALGORITHM_FACTORY(FilterSiliconEcalHits)
+
+FilterSiliconEcalHits::FilterSiliconEcalHits(const std::string& aName, ISvcLocator* aSvcLoc)
+ : GaudiAlgorithm(aName, aSvcLoc),
+ m_histSvc("THistSvc", "FilterSiliconEcalHits"),
+ m_geoSvc("GeoSvc", "FilterSiliconEcalHits")
+{
+ declareProperty("deposits", m_deposits, "Energy deposits in sampling calorimeter (input)");
+ declareProperty("filtered", m_filtered, "Filtered energy deposits (output)");
+
+}
+FilterSiliconEcalHits::~FilterSiliconEcalHits() {}
+
+StatusCode FilterSiliconEcalHits::initialize() {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
+ }
+ // check if readouts exist
+ if (m_geoSvc->lcdd()->readouts().find(m_padReadoutName) == m_geoSvc->lcdd()->readouts().end()) {
+ error() << "Readout <<" << m_padReadoutName << ">> does not exist." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ // if it gets to here then it has registered all histograms
+ return StatusCode::SUCCESS;
+}
+
+StatusCode FilterSiliconEcalHits::execute() {
+ // get the decoder for the readout
+ // this will allow to extract layers etc.
+ auto decoder = m_geoSvc->lcdd()->readout(m_padReadoutName).idSpec().decoder();
+
+ // added a little print out to ensure it is still working
+ if(m_EvtNumber < 10 or m_EvtNumber % 100 == 0) {
+ info() << "Event Number " << m_EvtNumber << endmsg;
+ }
+
+
+ // set some variables to hold information for an event
+ int nFilteredHits = 0;
+
+ const auto deposits = m_deposits.get();
+ auto filtered = m_filtered.createAndPut();
+
+ // if looking at digital need to threshold at this point so the merging of cells takes correct deposits
+ // the analgoue we can keep all hits and sum at a later date
+ // eventually this will be passed as a parameter but coding quickly now
+ double m_threshold = 0;
+ if(m_digitalFlag==1) {
+ m_threshold = 480*0.0000000036;
+ }
+
+ for (const auto& hit : *deposits) {
+ decoder->setValue(hit.core().cellId);
+ int digital = (*decoder)["digital"];
+ if(digital==m_digitalFlag) {
+ nFilteredHits++;
+
+ if(hit.energy()>m_threshold) {
+ fcc::PositionedCaloHit filteredHit = filtered->create();
+ filteredHit.energy(hit.energy());
+ filteredHit.time(hit.time());
+ filteredHit.cellId(hit.cellId());
+ filteredHit.bits(hit.bits());
+ filteredHit.position(hit.position());
+ }
+ }
+
+ }
+
+ info() << "nFilteredHits = " << nFilteredHits << endmsg;
+
+ m_EvtNumber++;
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode FilterSiliconEcalHits::finalize() {
+ debug() << "StatusCode FilterSiliconEcalHits::finalize()" << endmsg;
+
+
+return GaudiAlgorithm::finalize(); }
+
+
diff --git a/Detector/DetStudies/src/components/FilterSiliconEcalHits.h b/Detector/DetStudies/src/components/FilterSiliconEcalHits.h
new file mode 100644
index 000000000..7d7bc6f93
--- /dev/null
+++ b/Detector/DetStudies/src/components/FilterSiliconEcalHits.h
@@ -0,0 +1,70 @@
+#ifndef DETSTUDIES_FilterSiliconEcalHits_H
+#define DETSTUDIES_FilterSiliconEcalHits_H
+
+// GAUDI
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// FCCSW
+#include "FWCore/DataHandle.h"
+#include "TH2F.h"
+#include "TTree.h"
+#include "TF1.h"
+#include "TProfile.h"
+
+class IGeoSvc;
+
+// datamodel
+namespace fcc {
+class PositionedCaloHitCollection;
+}
+
+class ITHistSvc;
+/** @class DECalAnalysis DECalAnalysis.h
+ *
+ * Histograms of energy deposited in active material and total energy deposited in the calorimeter.
+ * 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 FilterSiliconEcalHits : public GaudiAlgorithm {
+public:
+ explicit FilterSiliconEcalHits(const std::string&, ISvcLocator*);
+ virtual ~FilterSiliconEcalHits();
+ /** 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 m_histSvc;
+ /// Pointer to the geometry service
+ ServiceHandle m_geoSvc;
+
+ /// The Pixel analysis
+ /// Handle for the energy deposits
+ DataHandle m_deposits{"rec/deposits", Gaudi::DataHandle::Reader, this};
+ DataHandle m_filtered{"rec/filtered", Gaudi::DataHandle::Writer, this};
+ /// Name of the detector readout
+ Gaudi::Property m_padReadoutName{this, "readoutName", "", "Name of the detector readout"};
+
+ Gaudi::Property m_digitalFlag{this, "digitalFlag", 0, "If 1 then filter for digital hits, 0 filter for analgoue hits"};
+
+ int m_EvtNumber;
+
+};
+#endif /* DETSTUDIES_FilterSiliconEcalHits_H */
diff --git a/Detector/DetStudies/src/components/SiWAnalysis.cpp b/Detector/DetStudies/src/components/SiWAnalysis.cpp
new file mode 100644
index 000000000..34312c6c2
--- /dev/null
+++ b/Detector/DetStudies/src/components/SiWAnalysis.cpp
@@ -0,0 +1,154 @@
+
+#include "SiWAnalysis.h"
+
+// FCCSW
+#include "DetInterface/IGeoSvc.h"
+
+// datamodel
+#include "datamodel/PositionedCaloHitCollection.h"
+#include "datamodel/CaloHitCollection.h"
+#include "datamodel/MCParticleCollection.h"
+
+#include "CLHEP/Vector/ThreeVector.h"
+#include "GaudiKernel/ITHistSvc.h"
+#include "TH1F.h"
+#include "TVector2.h"
+#include "TMath.h"
+#include "TString.h"
+#include "string"
+
+// DD4hep
+#include "DD4hep/Detector.h"
+#include "DD4hep/Readout.h"
+
+DECLARE_ALGORITHM_FACTORY(SiWAnalysis)
+
+SiWAnalysis::SiWAnalysis(const std::string& aName, ISvcLocator* aSvcLoc)
+ : GaudiAlgorithm(aName, aSvcLoc),
+ m_histSvc("THistSvc", "SiWAnalysis"),
+ m_geoSvc("GeoSvc", "SiWAnalysis")
+{
+ declareProperty("deposits", m_deposits, "Energy deposits in sampling calorimeter (input)");
+ declareProperty("truth", m_truth, "Generated particle truth");
+
+ m_calP0 = 239.74;
+ m_calP1 = 99.1602;
+ m_calP2 = -0.0143;
+
+}
+SiWAnalysis::~SiWAnalysis() {}
+
+StatusCode SiWAnalysis::initialize() {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
+ }
+ // check if readouts exist
+ if (m_geoSvc->lcdd()->readouts().find(m_padReadoutName) == m_geoSvc->lcdd()->readouts().end()) {
+ error() << "Readout <<" << m_padReadoutName << ">> does not exist." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+
+ m_treeInfo = new TTree("info", "Analysis Tree");
+ m_treeInfo->Branch("event_number", &m_treeEvtNumber);
+ m_treeInfo->Branch("truth_energy", &m_truthEnergy);
+ m_treeInfo->Branch("edep_tot", &m_sumEnergyDep);
+
+
+ if(m_histSvc->regTree("/rec/Info", m_treeInfo).isFailure()) {
+ error() << "Couldn't register tree" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ // create histograms
+
+ m_totalEdep = new TH1F("SiW_totalEdep", "Total edep in an event", 15000, 0, 15);
+ if (m_histSvc->regHist("/rec/SiW_Edep", m_totalEdep).isFailure()) {
+ error() << "Couldn't register histogram" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_calFit = new TF1("calFit", "pol2", 0,3000);
+ m_calFit->SetParameter("p0", m_calP0);
+ m_calFit->SetParameter("p1", m_calP1);
+ m_calFit->SetParameter("p2", m_calP2);
+
+ // if it gets to here then it has registered all histograms
+ return StatusCode::SUCCESS;
+}
+
+StatusCode SiWAnalysis::execute() {
+ // get the decoder for the readout
+ // this will allow to extract layers etc.
+ auto decoder = m_geoSvc->lcdd()->readout(m_padReadoutName).idSpec().decoder();
+
+ // calculate the MC truth energy of the incident particle
+ const auto truth = m_truth.get();
+ for (const auto& t : *truth) {
+ double px = t.core().p4.px;
+ double py = t.core().p4.py;
+ double pz = t.core().p4.pz;
+ double m = t.core().p4.mass;
+
+ m_truthEnergy = TMath::Sqrt(m*m + px*px + py*py + pz*pz);
+ }
+
+ // added a little print out to ensure it is still working
+ if(m_treeEvtNumber < 10 or m_treeEvtNumber % 100 == 0) {
+ info() << "Event Number " << m_treeEvtNumber << endmsg;
+ info() << "Truth energy of single particle = " << m_truthEnergy << endmsg;
+ }
+
+
+ // set some variables to hold information for an event
+ m_sumEdepLayers.assign(m_numLayers+1, 0);
+ m_sumEnergyDep = 0;
+
+ int isdigital = 0, notdigital = 0;
+
+ const auto deposits = m_deposits.get();
+ for (const auto& hit : *deposits) {
+ decoder->setValue(hit.core().cellId);
+
+
+ int layer = (*decoder)[m_layerFieldName];
+
+ // check if energy was deposited in the calorimeter (active/passive material)
+ // layers are numbered starting from 1, layer == 0 is cryostat/bath
+ if (layer > 0) {
+ decoder->setValue(hit.core().cellId);
+ int layer = (*decoder)[m_layerFieldName];
+ // get the energy deposited in the active area
+ double edep = hit.core().energy;
+ m_sumEnergyDep += edep;
+
+ }
+ }
+
+// Fill histograms
+ m_totalEdep->Fill(m_sumEnergyDep);
+
+ if(m_sumEnergyDep>0) {
+ m_treeInfo->Fill();
+ }
+
+ m_treeEvtNumber++;
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode SiWAnalysis::finalize() {
+ debug() << "StatusCode SiWAnalysis::finalize()" << endmsg;
+
+
+return GaudiAlgorithm::finalize(); }
+
+Double_t SiWAnalysis::FitLongProfile( Double_t *x, Double_t *par)
+{
+ double E_0 = par[0];
+ double a = par[1];
+ double b = par[2];
+ double led = E_0*b * (pow( (b*x[0]), a-1 )* exp(-(b*x[0])) ) / TMath::Gamma(a);
+
+ return led;
+}
diff --git a/Detector/DetStudies/src/components/SiWAnalysis.h b/Detector/DetStudies/src/components/SiWAnalysis.h
new file mode 100644
index 000000000..6ee703ef4
--- /dev/null
+++ b/Detector/DetStudies/src/components/SiWAnalysis.h
@@ -0,0 +1,98 @@
+#ifndef DETSTUDIES_SiWAnalysis_H
+#define DETSTUDIES_SiWAnalysis_H
+
+// GAUDI
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// FCCSW
+#include "FWCore/DataHandle.h"
+#include "TH2F.h"
+#include "TTree.h"
+#include "TF1.h"
+#include "TProfile.h"
+
+class IGeoSvc;
+
+// datamodel
+namespace fcc {
+class PositionedCaloHitCollection;
+class CaloHitCollection;
+class MCParticleCollection;
+}
+
+class TH1F;
+class ITHistSvc;
+/** @class DECalAnalysis DECalAnalysis.h
+ *
+ * Histograms of energy deposited in active material and total energy deposited in the calorimeter.
+ * 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 SiWAnalysis : public GaudiAlgorithm {
+public:
+ explicit SiWAnalysis(const std::string&, ISvcLocator*);
+ virtual ~SiWAnalysis();
+ /** 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;
+
+ Double_t FitLongProfile( Double_t *x, Double_t *par);
+
+private:
+ /// Pointer to the interface of histogram service
+ ServiceHandle m_histSvc;
+ /// Pointer to the geometry service
+ ServiceHandle m_geoSvc;
+ /// Number of layers/cells
+ Gaudi::Property m_numLayers{this, "numLayers", 50, "Number of layers"};
+ /// Name of the layer/cell field
+ Gaudi::Property m_layerFieldName{this, "layerFieldName", "", "Identifier of layers"};
+
+ /// The Pixel analysis
+ /// Handle for the energy deposits
+ DataHandle m_deposits{"rec/deposits", Gaudi::DataHandle::Reader, this};
+ /// Name of the detector readout
+ Gaudi::Property m_padReadoutName{this, "padReadoutName", "", "Name of the detector readout"};
+
+
+ DataHandle m_truth{"truth", Gaudi::DataHandle::Reader, this};
+
+
+
+ // calibration factors to compensate non linearity
+ float m_calP0, m_calP1, m_calP2;
+ TF1* m_calFit;
+
+ TTree* m_treeInfo;
+ int m_treeEvtNumber;
+
+ int m_truthEnergy;
+ float m_sumEnergyDep;
+ std::vector m_sumEdepLayers;
+
+ double m_fitE0;
+ double m_fita;
+ double m_fitb;
+ double m_rising50;
+ double m_rising90;
+ double m_fwhm;
+
+ TF1* m_fitLongProfile;
+ TH1F* m_totalEdep;
+};
+#endif /* DETSTUDIES_SiWAnalysis_H */
diff --git a/Detector/DetStudies/tests/options/SiWAnalysis_batch.py b/Detector/DetStudies/tests/options/SiWAnalysis_batch.py
new file mode 100644
index 000000000..d7446bf47
--- /dev/null
+++ b/Detector/DetStudies/tests/options/SiWAnalysis_batch.py
@@ -0,0 +1,96 @@
+from Gaudi.Configuration import *
+
+# Data service
+from Configurables import FCCDataSvc
+batch_dir = "/eos/user/t/toprice/private/FCC/FCCSW/v8.0/XYZ/"
+fccsw_version = "FCCSW0.8"
+det_config = ""
+run_config = ""
+file = ""
+file = file[7:]
+#inputfile = batch_dir+"/"+det_config+"_"+fccsw_version+"/"+run_config+"/"+file
+inputfile = batch_dir+"/"+det_config+"/"+run_config+"/output_"+file
+podiosvc = FCCDataSvc("EventDataSvc", input=inputfile)
+
+
+from Configurables import PodioInput
+podioinput = PodioInput("PodioReader", collections=["positionedCaloHits", "GenParticles"], OutputLevel=DEBUG)
+
+# DD4hep geometry service
+from Configurables import GeoSvc
+geoservice = GeoSvc("GeoSvc", detectors=[ 'file:/afs/cern.ch/user/t/toprice/private/FCC/FCCSW/Detector/DetFCChhBaseline1/compact/FCChh_DectEmptyMaster.xml',
+ 'file:/afs/cern.ch/user/t/toprice/private/FCC/FCCSW/Detector/DetFCChhECalDigital/compact/FCChh_DECalBarrel_'+det_config[:-9]+'.xml'
+],
+ OutputLevel = INFO)
+
+
+from Configurables import FilterSiliconEcalHits
+filtered = FilterSiliconEcalHits("FilterSiEcal",
+ readoutName = "BarDECal_Readout",
+ digitalFlag = 0)
+filtered.deposits.Path="positionedCaloHits"
+filtered.filtered.Path="filteredCaloHits"
+
+from Configurables import RedoSegmentation
+resegment = RedoSegmentation("ReSegmentation",
+ # old bitfield (readout)
+ oldReadoutName = "BarDECal_Readout",
+ # specify which fields are going to be deleted
+ oldSegmentationIds = ["x","y","z"],
+ # new bitfield (readout), with new segmentation
+ newReadoutName="BarDECal_Pads",
+
+ OutputLevel = INFO)
+# clusters are needed, with deposit position and cellID in bits
+resegment.inhits.Path = "filteredCaloHits"
+resegment.outhits.Path = "newCaloHits"
+
+from Configurables import CreateCaloCells
+createcells = CreateCaloCells("CreateCaloCells",
+ doCellCalibration = False,
+ addCellNoise = False, filterCellNoise = False,
+ OutputLevel = INFO)
+createcells.hits.Path="newCaloHits"
+createcells.cells.Path="newCaloCells"
+
+from Configurables import SiWAnalysis
+hist = SiWAnalysis("SiWAnalysis",
+ padReadoutName = "BarDECal_Pads",
+ layerFieldName = "layer",
+ numLayers = 50, # one more because index starts at 1 - layer 0 will be always empty
+ OutputLevel = INFO)
+
+hist.deposits.Path="newCaloCells"
+hist.truth.Path="GenParticles"
+
+
+
+#THistSvc().Output = ["rec DATAFILE='"+batch_dir+"/"+det_config+"_"+fccsw_version+"/"+run_config+"/