diff --git a/Detector/DetCEPCv4/CMakeLists.txt b/Detector/DetCEPCv4/CMakeLists.txt
index 19046a6924f3b444505c0389300cc71c283fafeb..03b897aece53d097dcbe979c271d06a39e007021 100644
--- a/Detector/DetCEPCv4/CMakeLists.txt
+++ b/Detector/DetCEPCv4/CMakeLists.txt
@@ -32,11 +32,13 @@ gaudi_add_module(DetCEPCv4
src/calorimeter/SHcalRpc02_Barrel.cpp
src/calorimeter/SHcalRpc01_Endcaps.cpp
src/calorimeter/SHcalRpc01_EndcapRing.cpp
- src/calorimeter/Yoke05_Barrel.cpp
- src/calorimeter/Yoke05_Endcaps.cpp
+ src/calorimeter/SHcalSc04_Barrel_v04.cpp
+ src/calorimeter/SHcalSc04_Endcaps_v01.cpp
+ src/calorimeter/Yoke05_Barrel.cpp
+ src/calorimeter/Yoke05_Endcaps.cpp
src/other/BoxSupport_o1_v01_geo.cpp
src/other/TubeSupport_o1_v01_geo.cpp
- src/other/SCoil02_geo.cpp
+ src/other/SCoil02_geo.cpp
LINK ${DD4hep_COMPONENT_LIBRARIES}
)
diff --git a/Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Barrel_v04.cpp b/Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Barrel_v04.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..334af303d765b12f008d7eccad8a8f4c02e036ee
--- /dev/null
+++ b/Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Barrel_v04.cpp
@@ -0,0 +1,723 @@
+//====================================================================
+// DD4hep Geometry driver for HcalBarrel
+//--------------------------------------------------------------------
+// S.Lu, DESY
+// F. Gaede, DESY : v04 : prepare for multi segmentation
+// 18.04.2017 - copied from SHcalSc04_Barrel_v01.cpp
+// - add optional parameter <subsegmentation key="" value=""/>
+// defines the segmentation to be used in reconstruction
+//====================================================================
+#include "DD4hep/Printout.h"
+#include "DD4hep/DetFactoryHelper.h"
+#include "XML/Utilities.h"
+#include "DDRec/DetectorData.h"
+#include "DDSegmentation/BitField64.h"
+#include "DDSegmentation/TiledLayerGridXY.h"
+#include "DDSegmentation/Segmentation.h"
+#include "DDSegmentation/MultiSegmentation.h"
+#include "LcgeoExceptions.h"
+
+#include <iostream>
+#include <vector>
+
+using namespace std;
+
+using dd4hep::BUILD_ENVELOPE;
+using dd4hep::Box;
+using dd4hep::DetElement;
+using dd4hep::Detector;
+using dd4hep::IntersectionSolid;
+using dd4hep::Material;
+using dd4hep::PlacedVolume;
+using dd4hep::Position;
+using dd4hep::Readout;
+using dd4hep::Ref_t;
+using dd4hep::Rotation3D;
+using dd4hep::RotationZ;
+using dd4hep::RotationZYX;
+using dd4hep::SensitiveDetector;
+using dd4hep::Transform3D;
+using dd4hep::Trapezoid;
+using dd4hep::Tube;
+using dd4hep::Volume;
+using dd4hep::_toString;
+
+using dd4hep::rec::LayeredCalorimeterData;
+
+// After reading in all the necessary parameters.
+// To check the radius range and the space for placing the total layers
+static bool validateEnvelope(double rInner, double rOuter, double radiatorThickness, double layerThickness, int layerNumber){
+
+ bool Error = false;
+ bool Warning = false;
+ double spaceAllowed = rOuter*cos(M_PI/16.) - rInner;
+ double spaceNeeded = (radiatorThickness + layerThickness)* layerNumber;
+ double spaceToleranted = (radiatorThickness + layerThickness)* (layerNumber+1);
+ double rOuterRecommaned = ( rInner + spaceNeeded )/cos(M_PI/16.);
+ int layerNumberRecommaned = floor( ( spaceAllowed ) / (radiatorThickness + layerThickness) );
+
+
+ if( spaceNeeded > spaceAllowed )
+ {
+ printout( dd4hep::ERROR, "SHcalSc04_Barrel_v01", " Layer number is more than it can be built! " ) ;
+ Error = true;
+ }
+ else if ( spaceToleranted < spaceAllowed )
+ {
+ printout( dd4hep::WARNING, "SHcalSc04_Barrel_v01", " Layer number is less than it is able to build!" ) ;
+ Warning = true;
+ }
+ else
+ {
+ printout( dd4hep::DEBUG, "SHcalSc04_Barrel_v01"," has been validated and start to build it." ) ;
+ Error = false;
+ Warning = false;
+ }
+
+ if( Error )
+ {
+ cout<<"\n ============> First Help Documentation <=============== \n"
+ <<" When you see this message, that means you are crashing the module. \n"
+ <<" Please take a cup of cafe, and think about what you want to do! \n"
+ <<" Here are few FirstAid# for you. Please read them carefully. \n"
+ <<" \n"
+ <<" ### FirstAid 1: ###\n"
+ <<" If you want to build HCAL within the rInner and rOuter range, \n"
+ <<" please reduce the layer number to " << layerNumberRecommaned <<" \n"
+ <<" with the current layer thickness structure. \n"
+ <<" \n"
+ <<" You may redisgn the layer structure and thickness, too. \n"
+ <<" \n"
+ <<" ### FirstAid 2: ###\n"
+ <<" If you want to build HCAL with this layer number and the layer thickness, \n"
+ <<" you have to update rOuter to "<< rOuterRecommaned*10. <<"*mm \n"
+ <<" and to inform other subdetector, you need this space for building your design. \n"
+ <<" \n"
+ <<" ### FirstAid 3: ###\n"
+ <<" Do you think that you are looking for another type of HCAL driver? \n"
+ <<" \n"
+ <<endl;
+ throw lcgeo::GeometryException( "SHcalSc04_Barrel: Error: Layer number is more than it can be built!" ) ;
+ }
+ else if( Warning )
+ {
+ cout<<"\n ============> First Help Documentation <=============== \n"
+ <<" When you see this warning message, that means you are changing the module. \n"
+ <<" Please take a cup of cafe, and think about what you want to do! \n"
+ <<" Here are few FirstAid# for you. Please read them carefully. \n"
+ <<" \n"
+ <<" ### FirstAid 1: ###\n"
+ <<" If you want to build HCAL within the rInner and rOuter range, \n"
+ <<" You could build the layer number up to " << layerNumberRecommaned <<" \n"
+ <<" with the current layer thickness structure. \n"
+ <<" \n"
+ <<" You may redisgn the layer structure and thickness, too. \n"
+ <<" \n"
+ <<" ### FirstAid 2: ###\n"
+ <<" If you want to build HCAL with this layer number and the layer thickness, \n"
+ <<" you could reduce rOuter to "<< rOuterRecommaned*10. <<"*mm \n"
+ <<" and to reduce the back plate thickness, which you may not need for placing layer. \n"
+ <<" \n"
+ <<" ### FirstAid 3: ###\n"
+ <<" Do you think that you are looking for another type of HCAL driver? \n"
+ <<" \n"
+ <<endl;
+ return Warning;
+ }
+ else { return true; }
+
+}
+
+static Ref_t create_detector(Detector& theDetector, xml_h element, SensitiveDetector sens) {
+
+ double boundarySafety = 0.0001;
+
+ xml_det_t x_det = element;
+ string det_name = x_det.nameStr();
+ int det_id = x_det.id();
+ DetElement sdet( det_name, det_id );
+
+ // --- create an envelope volume and position it into the world ---------------------
+
+ Volume envelope = dd4hep::xml::createPlacedEnvelope( theDetector, element , sdet ) ;
+
+ dd4hep::xml::setDetectorTypeFlag( element, sdet ) ;
+
+ if( theDetector.buildType() == BUILD_ENVELOPE ) return sdet ;
+
+ //-----------------------------------------------------------------------------------
+
+ xml_comp_t x_staves = x_det.staves();
+ Material stavesMaterial = theDetector.material(x_staves.materialStr());
+ Material air = theDetector.air();
+
+ PlacedVolume pv;
+
+ sens.setType("calorimeter");
+
+//====================================================================
+//
+// Read all the constant from ILD_o1_v05.xml
+// Use them to build HcalBarrel
+//
+//====================================================================
+ double Hcal_inner_radius = theDetector.constant<double>("Hcal_inner_radius");
+ double Hcal_outer_radius = theDetector.constant<double>("Hcal_outer_radius");
+ double Hcal_half_length = theDetector.constant<double>("Hcal_half_length");
+ int Hcal_inner_symmetry = theDetector.constant<int>("Hcal_inner_symmetry");
+ int Hcal_outer_symmetry = 0; // Fixed shape for Tube, and not allow to modify from compact xml.
+
+ double Hcal_radiator_thickness = theDetector.constant<double>("Hcal_radiator_thickness");
+ double Hcal_chamber_thickness = theDetector.constant<double>("Hcal_chamber_thickness");
+ double Hcal_back_plate_thickness = theDetector.constant<double>("Hcal_back_plate_thickness");
+ double Hcal_lateral_plate_thickness = theDetector.constant<double>("Hcal_lateral_structure_thickness");
+ double Hcal_stave_gaps = theDetector.constant<double>("Hcal_stave_gaps");
+ double Hcal_modules_gap = theDetector.constant<double>("Hcal_modules_gap");
+ double Hcal_middle_stave_gaps = theDetector.constant<double>("Hcal_middle_stave_gaps");
+ double Hcal_layer_air_gap = theDetector.constant<double>("Hcal_layer_air_gap");
+ //double Hcal_cells_size = theDetector.constant<double>("Hcal_cells_size");
+
+ int Hcal_nlayers = theDetector.constant<int>("Hcal_nlayers");
+
+ double TPC_outer_radius = theDetector.constant<double>("TPC_outer_radius");
+
+ double Ecal_outer_radius = theDetector.constant<double>("Ecal_outer_radius");
+
+ printout( dd4hep::DEBUG, "SHcalSc04_Barrel_v04", "TPC_outer_radius : %e - Ecal_outer_radius: %e ", TPC_outer_radius , Ecal_outer_radius) ;
+
+ validateEnvelope(Hcal_inner_radius, Hcal_outer_radius, Hcal_radiator_thickness, Hcal_chamber_thickness, Hcal_nlayers);
+
+ Readout readout = sens.readout();
+ dd4hep::Segmentation seg = readout.segmentation();
+
+ dd4hep::DDSegmentation::BitField64 encoder = seg.decoder();
+ encoder.setValue(0) ;
+
+
+ //fg: this is a bit tricky: we first have to check whether a multi segmentation is used and then, if so, we
+ // will check the <subsegmentation key="" value=""/> element, for which subsegmentation to use for filling the
+ // DDRec:LayeredCalorimeterData information.
+ // Additionally, we need to figure out if there is a TiledLayerGridXY instance defined -
+ // and in case, there is more than one defined, we need to pick the correct one specified via the
+ // <subsegmentation key="" value=""/> element.
+ // This involves a lot of casting: review the API in DD4hep !!
+
+ // check if we have a multi segmentation :
+
+ dd4hep::DDSegmentation::MultiSegmentation* multiSeg =
+ dynamic_cast< dd4hep::DDSegmentation::MultiSegmentation*>( seg.segmentation() ) ;
+
+ dd4hep::DDSegmentation::TiledLayerGridXY* tileSeg = 0 ;
+
+ int sensitive_slice_number = -1 ;
+
+ if( multiSeg ){
+
+ try{
+ // check if we have an entry for the subsegmentation to be used
+ xml_comp_t segxml = x_det.child( _Unicode( subsegmentation ) ) ;
+
+ std::string keyStr = segxml.attr<std::string>( _Unicode(key) ) ;
+ int keyVal = segxml.attr<int>( _Unicode(value) ) ;
+
+ encoder[ keyStr ] = keyVal ;
+
+ // if we have a multisegmentation that uses the slice as key, we need to know for the
+ // computation of the layer parameters in LayeredCalorimeterData::Layer below
+ if( keyStr == "slice" ){
+ sensitive_slice_number = keyVal ;
+ }
+
+
+ } catch(const std::runtime_error &) {
+ throw lcgeo::GeometryException( "SHcalSc04_Barrel: Error: MultiSegmentation specified but no "
+ " <subsegmentation key="" value=""/> element defined for detector ! " ) ;
+ }
+
+ // check if we have a TiledLayerGridXY segmentation :
+ const dd4hep::DDSegmentation::TiledLayerGridXY* ts0 =
+ dynamic_cast<const dd4hep::DDSegmentation::TiledLayerGridXY*>( &multiSeg->subsegmentation( encoder.getValue() ) ) ;
+
+ tileSeg = const_cast<dd4hep::DDSegmentation::TiledLayerGridXY*>( ts0 ) ;
+
+ if( ! tileSeg ){ // if the current segmentation is not a tileSeg, we see if there is another one
+
+ for( auto s : multiSeg->subSegmentations() ){
+ const dd4hep::DDSegmentation::TiledLayerGridXY* ts =
+ dynamic_cast<const dd4hep::DDSegmentation::TiledLayerGridXY*>( s.segmentation ) ;
+
+ if( ts ) {
+ tileSeg = const_cast<dd4hep::DDSegmentation::TiledLayerGridXY*>( ts ) ;
+ break ;
+ }
+ }
+ }
+
+ } else {
+
+ tileSeg =
+ dynamic_cast< dd4hep::DDSegmentation::TiledLayerGridXY*>( seg.segmentation() ) ;
+ }
+
+
+
+ std::vector<double> cellSizeVector = seg.cellDimensions( encoder.getValue() ); //Assume uniform cell sizes, provide dummy cellID
+ double cell_sizeX = cellSizeVector[0];
+ double cell_sizeY = cellSizeVector[1];
+
+
+ //========== fill data for reconstruction ============================
+ LayeredCalorimeterData* caloData = new LayeredCalorimeterData ;
+ caloData->layoutType = LayeredCalorimeterData::BarrelLayout ;
+ caloData->inner_symmetry = Hcal_inner_symmetry ;
+ caloData->outer_symmetry = Hcal_outer_symmetry ;
+ caloData->phi0 = 0 ; // fg: also hardcoded below
+
+ /// extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
+ caloData->extent[0] = Hcal_inner_radius ;
+ caloData->extent[1] = Hcal_outer_radius ;
+ caloData->extent[2] = 0. ; // Barrel zmin is "0" by default.
+ caloData->extent[3] = Hcal_half_length ;
+
+//====================================================================
+//
+// general calculated parameters
+//
+//====================================================================
+
+ double Hcal_total_dim_y = Hcal_nlayers * (Hcal_radiator_thickness + Hcal_chamber_thickness)
+ + Hcal_back_plate_thickness;
+
+ double Hcal_y_dim1_for_x = Hcal_outer_radius*cos(M_PI/Hcal_inner_symmetry) - Hcal_inner_radius;
+ double Hcal_bottom_dim_x = 2.*Hcal_inner_radius*tan(M_PI/Hcal_inner_symmetry)- Hcal_stave_gaps;
+ double Hcal_normal_dim_z = (2 * Hcal_half_length - Hcal_modules_gap)/2.;
+
+ //only the middle has the steel plate.
+ double Hcal_regular_chamber_dim_z = Hcal_normal_dim_z - Hcal_lateral_plate_thickness;
+
+ //double Hcal_cell_dim_x = Hcal_cells_size;
+ //double Hcal_cell_dim_z = Hcal_regular_chamber_dim_z / floor (Hcal_regular_chamber_dim_z/Hcal_cell_dim_x);
+
+
+// ========= Create Hcal Barrel stave ====================================
+// It will be the volume for palcing the Hcal Barrel Chamber(i.e. Layers).
+// Itself will be placed into the world volume.
+// ==========================================================================
+
+ double chambers_y_off_correction = 0.;
+
+ // stave modules shaper parameters
+ double BHX = (Hcal_bottom_dim_x + Hcal_stave_gaps)/2.;
+ double THX = (Hcal_total_dim_y + Hcal_inner_radius)*tan(M_PI/Hcal_inner_symmetry);
+ double YXH = Hcal_total_dim_y / 2.;
+ double DHZ = (Hcal_normal_dim_z - Hcal_lateral_plate_thickness) / 2.;
+
+ Trapezoid stave_shaper( THX, BHX, DHZ, DHZ, YXH);
+
+ Tube solidCaloTube(0, Hcal_outer_radius, DHZ+boundarySafety);
+
+ RotationZYX mrot(0,0,M_PI/2.);
+
+ Rotation3D mrot3D(mrot);
+ Position mxyzVec(0,0,(Hcal_inner_radius + Hcal_total_dim_y / 2.));
+ Transform3D mtran3D(mrot3D,mxyzVec);
+
+ IntersectionSolid barrelModuleSolid(stave_shaper, solidCaloTube, mtran3D);
+
+ Volume EnvLogHcalModuleBarrel(det_name+"_module",barrelModuleSolid,stavesMaterial);
+
+ EnvLogHcalModuleBarrel.setAttributes(theDetector,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
+
+
+
+
+
+ //stave modules lateral plate shaper parameters
+ double BHX_LP = BHX;
+ double THX_LP = THX;
+ double YXH_LP = YXH;
+
+ //build lateral palte here to simulate lateral plate in the middle of barrel.
+ double DHZ_LP = Hcal_lateral_plate_thickness/2.0;
+
+ Trapezoid stave_shaper_LP(THX_LP, BHX_LP, DHZ_LP, DHZ_LP, YXH_LP);
+
+ Tube solidCaloTube_LP(0, Hcal_outer_radius, DHZ_LP+boundarySafety);
+
+ IntersectionSolid Module_lateral_plate(stave_shaper_LP, solidCaloTube_LP, mtran3D);
+
+ Volume EnvLogHcalModuleBarrel_LP(det_name+"_Module_lateral_plate",Module_lateral_plate,stavesMaterial);
+
+ EnvLogHcalModuleBarrel_LP.setAttributes(theDetector,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
+
+
+
+#ifdef SHCALSC04_DEBUG
+ std::cout<< " ==> Hcal_outer_radius: "<<Hcal_outer_radius <<std::endl;
+#endif
+
+
+
+
+
+//====================================================================
+//
+// Chambers in the HCAL BARREL
+//
+//====================================================================
+ // Build Layer Chamber fill with air, which include the tolarance space at the two side
+ // place the slice into the Layer Chamber
+ // Place the Layer Chamber into the Stave module
+ // place the Stave module into the asembly Volume
+ // place the module middle lateral plate into asembly Volume
+
+ double x_length; //dimension of an Hcal barrel layer on the x-axis
+ double y_height; //dimension of an Hcal barrel layer on the y-axis
+ double z_width; //dimension of an Hcal barrel layer on the z-axis
+ x_length = 0.; // Each Layer the x_length has new value.
+ y_height = Hcal_chamber_thickness / 2.;
+ z_width = Hcal_regular_chamber_dim_z/2.;
+
+ double xOffset = 0.;//the x_length of a barrel layer is calculated as a
+ //barrel x-dimension plus (bottom barrel) or minus
+ //(top barrel) an x-offset, which depends on the angle M_PI/Hcal_inner_symmetry
+
+ double xShift = 0.;//Geant4 draws everything in the barrel related to the
+ //center of the bottom barrel, so we need to shift the layers to
+ //the left (or to the right) with the quantity xShift
+
+
+ //-------------------- start loop over HCAL layers ----------------------
+
+ for (int layer_id = 1; layer_id <= (2*Hcal_nlayers); layer_id++)
+ {
+
+ double TanPiDiv8 = tan(M_PI/Hcal_inner_symmetry);
+ double x_total = 0.;
+ double x_halfLength;
+ x_length = 0.;
+
+ int logical_layer_id = 0;
+
+ if ( (layer_id < Hcal_nlayers)
+ || (layer_id > Hcal_nlayers && layer_id < (2*Hcal_nlayers)) )
+ logical_layer_id = layer_id % Hcal_nlayers;
+ else if ( (layer_id == Hcal_nlayers)
+ || (layer_id == 2*Hcal_nlayers) ) logical_layer_id = Hcal_nlayers;
+
+ //---- bottom barrel------------------------------------------------------------
+ if( logical_layer_id *(Hcal_radiator_thickness + Hcal_chamber_thickness)
+ < (Hcal_outer_radius * cos(M_PI/Hcal_inner_symmetry) - Hcal_inner_radius ) ) {
+ xOffset = (logical_layer_id * Hcal_radiator_thickness
+ + (logical_layer_id -1) * Hcal_chamber_thickness) * TanPiDiv8;
+
+ x_total = Hcal_bottom_dim_x/2 - Hcal_middle_stave_gaps/2 + xOffset;
+ x_length = x_total - 2*Hcal_layer_air_gap;
+ x_halfLength = x_length/2.;
+
+ } else {//----- top barrel -------------------------------------------------
+ double y_layerID = logical_layer_id * (Hcal_radiator_thickness + Hcal_chamber_thickness) + Hcal_inner_radius;
+ double ro_layer = Hcal_outer_radius - Hcal_radiator_thickness;
+
+ x_total = sqrt( ro_layer * ro_layer - y_layerID * y_layerID);
+
+ x_length = x_total - Hcal_middle_stave_gaps;
+
+ x_halfLength = x_length/2.;
+
+ xOffset = (logical_layer_id * Hcal_radiator_thickness
+ + (logical_layer_id - 1) * Hcal_chamber_thickness - Hcal_y_dim1_for_x) / TanPiDiv8
+ + Hcal_chamber_thickness / TanPiDiv8;
+
+ }
+
+ double xAbsShift = (Hcal_middle_stave_gaps/2 + Hcal_layer_air_gap + x_halfLength);
+
+ if (layer_id <= Hcal_nlayers) xShift = - xAbsShift;
+ else if (layer_id > Hcal_nlayers) xShift = xAbsShift;
+
+ x_length = x_length/2.;
+
+
+ //calculate the size of a fractional tile
+ //-> this sets fract_cell_dim_x
+
+ //double fract_cell_dim_x = 0.;
+ //this->CalculateFractTileSize(2*x_length, Hcal_cell_dim_x, fract_cell_dim_x);
+
+ //Vector newFractCellDim(fract_cell_dim_x, Hcal_chamber_thickness, Hcal_cell_dim_z);
+ //theBarrilRegSD->SetFractCellDimPerLayer(layer_id, newFractCellDim);
+
+ encoder["layer"] = logical_layer_id ;
+ cellSizeVector = seg.segmentation()->cellDimensions( encoder.getValue() );
+ cell_sizeX = cellSizeVector[0];
+ cell_sizeY = cellSizeVector[1];
+
+ LayeredCalorimeterData::Layer caloLayer ;
+ caloLayer.cellSize0 = cell_sizeX;
+ caloLayer.cellSize1 = cell_sizeY;
+
+
+ //--------------------------------------------------------------------------------
+ // build chamber box, with the calculated dimensions
+ //-------------------------------------------------------------------------------
+ printout( dd4hep::DEBUG, "SHcalSc04_Barrel_v04", " \n Start to build layer chamber - layer_id: %d", layer_id ) ;
+ printout( dd4hep::DEBUG, "SHcalSc04_Barrel_v04"," chamber x:y:z: %e:%e:%e", x_length*2., z_width*2. , y_height*2. );
+
+ //check if x_length (scintillator length) is divisible with x_integerTileSize
+ if( layer_id <= Hcal_nlayers) {
+ double fracPart, intPart;
+ double temp = x_length*2./cell_sizeX;
+ fracPart = modf(temp, &intPart);
+ int noOfIntCells = int(temp);
+
+
+ if( tileSeg !=0 ){
+
+ tileSeg->setBoundaryLayerX(x_length);
+
+ if (fracPart == 0){ //divisible
+ if ( noOfIntCells%2 ) {
+ if( tileSeg !=0 ) tileSeg->setLayerOffsetX(0);
+ }
+ else {
+ if( tileSeg !=0 ) tileSeg->setLayerOffsetX(1);
+ }
+ tileSeg->setFractCellSizeXPerLayer(0);
+ }
+ else if (fracPart>0){
+ if ( noOfIntCells%2 ) {
+ if( tileSeg !=0 ) tileSeg->setLayerOffsetX(1);
+ }
+ else {
+ if( tileSeg !=0 ) tileSeg->setLayerOffsetX(0);
+ }
+ tileSeg->setFractCellSizeXPerLayer( (fracPart+1.0)/2.0*cell_sizeX );
+ }
+
+ if ( (int)( (z_width*2.) / cell_sizeX)%2 ){
+ if( tileSeg !=0 ) tileSeg->setLayerOffsetY(0);
+ }
+ else {
+ if( tileSeg !=0 ) tileSeg->setLayerOffsetY(1);
+ }
+
+ }
+ }
+ Box ChamberSolid((x_length + Hcal_layer_air_gap), //x + air gaps at two side, do not need to build air gaps individualy.
+ z_width, //z attention!
+ y_height); //y attention!
+
+ string ChamberLogical_name = det_name+_toString(layer_id,"_layer%d");
+
+ Volume ChamberLogical(ChamberLogical_name, ChamberSolid, air);
+
+
+
+ double layer_thickness = y_height*2.;
+
+ double nRadiationLengths=0.;
+ double nInteractionLengths=0.;
+ double thickness_sum=0;
+
+ nRadiationLengths = Hcal_radiator_thickness/(stavesMaterial.radLength());
+ nInteractionLengths = Hcal_radiator_thickness/(stavesMaterial.intLength());
+ thickness_sum = Hcal_radiator_thickness;
+
+//====================================================================
+// Create Hcal Barrel Chamber without radiator
+// Place into the Hcal Barrel stave, after each radiator
+//====================================================================
+ xml_coll_t c(x_det,_U(layer));
+ xml_comp_t x_layer = c;
+ string layer_name = det_name+_toString(layer_id,"_layer%d");
+
+ // Create the slices (sublayers) within the Hcal Barrel Chamber.
+ double slice_pos_z = layer_thickness/2.;
+ int slice_number = 0;
+
+ for(xml_coll_t k(x_layer,_U(slice)); k; ++k) {
+ xml_comp_t x_slice = k;
+ string slice_name = layer_name + _toString(slice_number,"_slice%d");
+ double slice_thickness = x_slice.thickness();
+ Material slice_material = theDetector.material(x_slice.materialStr());
+ DetElement slice(layer_name,_toString(slice_number,"slice%d"),x_det.id());
+
+ slice_pos_z -= slice_thickness/2.;
+
+ // Slice volume & box
+ Volume slice_vol(slice_name,Box(x_length,z_width,slice_thickness/2.),slice_material);
+
+ nRadiationLengths += slice_thickness/(2.*slice_material.radLength());
+ nInteractionLengths += slice_thickness/(2.*slice_material.intLength());
+ thickness_sum += slice_thickness/2;
+
+ if ( x_slice.isSensitive() ) {
+
+ slice_vol.setSensitiveDetector(sens);
+
+ // if we have a multisegmentation based on slices, we need to use the correct slice here
+ if ( sensitive_slice_number<0 || sensitive_slice_number == slice_number ) {
+
+ //Store "inner" quantities
+ caloLayer.inner_nRadiationLengths = nRadiationLengths;
+ caloLayer.inner_nInteractionLengths = nInteractionLengths;
+ caloLayer.inner_thickness = thickness_sum;
+ //Store scintillator thickness
+ caloLayer.sensitive_thickness = slice_thickness;
+
+ //Reset counters to measure "outside" quantitites
+ nRadiationLengths=0.;
+ nInteractionLengths=0.;
+ thickness_sum = 0.;
+ }
+ }
+
+ nRadiationLengths += slice_thickness/(2.*slice_material.radLength());
+ nInteractionLengths += slice_thickness/(2.*slice_material.intLength());
+ thickness_sum += slice_thickness/2;
+
+
+ // Set region, limitset, and vis.
+ slice_vol.setAttributes(theDetector,x_slice.regionStr(),x_slice.limitsStr(),x_slice.visStr());
+ // slice PlacedVolume
+ PlacedVolume slice_phv = ChamberLogical.placeVolume(slice_vol,Position(0.,0.,slice_pos_z));
+
+ slice_phv.addPhysVolID("layer",logical_layer_id).addPhysVolID("slice", slice_number );
+
+
+ if ( x_slice.isSensitive() ) {
+ int tower_id = (layer_id > Hcal_nlayers)? 1:-1;
+ slice_phv.addPhysVolID("tower",tower_id);
+ printout( dd4hep::DEBUG, "SHcalSc04_Barrel_v04", " logical_layer_id: %d tower_id: %d", logical_layer_id, tower_id ) ;
+ }
+
+ slice.setPlacement(slice_phv);
+ // Increment x position for next slice.
+ slice_pos_z -= slice_thickness/2.;
+ // Increment slice number.
+ ++slice_number;
+ }
+
+ //Store "outer" quantities
+ caloLayer.outer_nRadiationLengths = nRadiationLengths;
+ caloLayer.outer_nInteractionLengths = nInteractionLengths;
+ caloLayer.outer_thickness = thickness_sum;
+
+
+//--------------------------- Chamber Placements -----------------------------------------
+
+ double chamber_x_offset, chamber_y_offset, chamber_z_offset;
+ chamber_x_offset = xShift;
+
+ chamber_z_offset = 0;
+
+
+ chamber_y_offset = -(-Hcal_total_dim_y/2.
+ + (logical_layer_id-1) *(Hcal_chamber_thickness + Hcal_radiator_thickness)
+ + Hcal_radiator_thickness + Hcal_chamber_thickness/2.);
+
+
+ pv = EnvLogHcalModuleBarrel.placeVolume(ChamberLogical,
+ Position(chamber_x_offset,
+ chamber_z_offset,
+ chamber_y_offset + chambers_y_off_correction));
+
+
+
+ //-----------------------------------------------------------------------------------------
+ if( layer_id <= Hcal_nlayers ){ // add the first set of layers to the reconstruction data
+
+ caloLayer.distance = Hcal_inner_radius + Hcal_total_dim_y/2.0 - (chamber_y_offset + chambers_y_off_correction)
+ - caloLayer.inner_thickness ; // Will be added later at "DDMarlinPandora/DDGeometryCreator.cc:226" to get center of sensitive element
+
+ caloLayer.absorberThickness = Hcal_radiator_thickness ;
+
+ caloData->layers.push_back( caloLayer ) ;
+ }
+ //-----------------------------------------------------------------------------------------
+
+
+ }//end loop over HCAL nlayers;
+
+ if( tileSeg !=0 ){
+ // check the offsets directly in the TileSeg ...
+ std::vector<double> LOX = tileSeg->layerOffsetX();
+ std::vector<double> LOY = tileSeg->layerOffsetY();
+
+ std::stringstream sts ;
+ sts <<" layerOffsetX(): ";
+ for (std::vector<double>::const_iterator ix = LOX.begin(); ix != LOX.end(); ++ix)
+ sts << *ix << ' ';
+ printout( dd4hep::DEBUG, "SHcalSc04_Barrel_v04", "%s" , sts.str().c_str() ) ;
+ sts.clear() ; sts.str("") ;
+ sts <<" layerOffsetY(): ";
+ for (std::vector<double>::const_iterator iy = LOY.begin(); iy != LOY.end(); ++iy)
+ sts << *iy << ' ';
+ printout( dd4hep::DEBUG, "SHcalSc04_Barrel_v04", "%s" , sts.str().c_str() ) ;
+
+ }
+
+
+
+//====================================================================
+// Place HCAL Barrel stave module into the envelope
+//====================================================================
+ double stave_phi_offset, module_z_offset, lateral_plate_z_offset;
+
+ double Y = Hcal_inner_radius + Hcal_total_dim_y / 2.;
+
+ stave_phi_offset = M_PI/Hcal_inner_symmetry -M_PI/2.;
+
+
+ //-------- start loop over HCAL BARREL staves ----------------------------
+
+ for (int stave_id = 1;
+ stave_id <=Hcal_inner_symmetry;
+ stave_id++)
+ {
+ module_z_offset = - (Hcal_normal_dim_z + Hcal_modules_gap + Hcal_lateral_plate_thickness)/2.;
+ lateral_plate_z_offset = - (Hcal_lateral_plate_thickness + Hcal_modules_gap)/2.;
+
+ double phirot = stave_phi_offset;
+ RotationZYX srot(0,phirot,M_PI*0.5);
+ Rotation3D srot3D(srot);
+
+ for (int module_id = 1;
+ module_id <=2;
+ module_id++)
+ {
+ Position sxyzVec(-Y*sin(phirot), Y*cos(phirot), module_z_offset);
+ Transform3D stran3D(srot3D,sxyzVec);
+
+ // Place Hcal Barrel volume into the envelope volume
+ pv = envelope.placeVolume(EnvLogHcalModuleBarrel,stran3D);
+ pv.addPhysVolID("stave",stave_id);
+ pv.addPhysVolID("module",module_id);
+ pv.addPhysVolID("system",det_id);
+
+ const int staveCounter = (stave_id-1)*2+module_id-1;
+ DetElement stave(sdet, _toString(staveCounter,"stave%d"), staveCounter );
+ stave.setPlacement(pv);
+
+ Position xyzVec_LP(-Y*sin(phirot), Y*cos(phirot),lateral_plate_z_offset);
+ Transform3D tran3D_LP(srot3D,xyzVec_LP);
+ pv = envelope.placeVolume(EnvLogHcalModuleBarrel_LP,tran3D_LP);
+
+ module_z_offset = - module_z_offset;
+ lateral_plate_z_offset = - lateral_plate_z_offset;
+ }
+
+
+ stave_phi_offset -= M_PI*2.0/Hcal_inner_symmetry;
+ } //-------- end loop over HCAL BARREL staves ----------------------------
+
+
+ sdet.addExtension< LayeredCalorimeterData >( caloData ) ;
+
+
+ return sdet;
+
+}
+
+DECLARE_DETELEMENT(SHcalSc04_Barrel_v04, create_detector)
diff --git a/Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Endcaps_v01.cpp b/Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Endcaps_v01.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b074089bdb71a1938347b3cdd9008d3b598a2ba
--- /dev/null
+++ b/Detector/DetCEPCv4/src/calorimeter/SHcalSc04_Endcaps_v01.cpp
@@ -0,0 +1,498 @@
+//====================================================================
+// AIDA Detector description implementation
+// for LDC AHCAL Endcap
+//--------------------------------------------------------------------
+//
+// Author : S.Lu
+// F. Gaede, DESY : v01 : prepare for multi segmentation
+// 18.04.2017 - copied from SHcalSc04_Barrel_v01.cpp
+// - add optional parameter <subsegmentation key="" value=""/>
+// defines the segmentation to be used in reconstruction
+//
+// Basic idea:
+// 1. Create the Hcal Endcap module envelope (16 modules).
+// Note: with default material Steel235.
+//
+// 2. Create the Hcal Endcap Chamber(i.e. Layer) for each module.
+// Create the Layer with slices (Polystyrene,Cu,FR4,air).
+// Place each slice into the chamber with the right position,
+// And registry the IDs for slice
+//
+// 3. Place the same Layer into the endcap module envelope.
+// It will be repeated repeat 48 times.
+// And registry the IDs for layer, and endcapID.
+//
+// 4. Place the endcap module into the world volume,
+// with the right position and rotation.
+// And registry the IDs for stave,module and endcapID.
+//
+// 5. Customer material FR4 and Steel235 defined in materials.xml
+//
+//====================================================================
+#include "DD4hep/DetFactoryHelper.h"
+#include "DD4hep/DetType.h"
+#include "XML/Layering.h"
+#include "XML/Utilities.h"
+#include "DDRec/DetectorData.h"
+#include "DDSegmentation/BitField64.h"
+#include "DDSegmentation/Segmentation.h"
+#include "DDSegmentation/MultiSegmentation.h"
+#include "LcgeoExceptions.h"
+
+using namespace std;
+
+using dd4hep::BUILD_ENVELOPE;
+using dd4hep::BitField64;
+using dd4hep::Box;
+using dd4hep::DetElement;
+using dd4hep::DetType;
+using dd4hep::Detector;
+using dd4hep::Layering;
+using dd4hep::Material;
+using dd4hep::PlacedVolume;
+using dd4hep::Position;
+using dd4hep::Readout;
+using dd4hep::Ref_t;
+using dd4hep::RotationX;
+using dd4hep::Segmentation;
+using dd4hep::SensitiveDetector;
+using dd4hep::Transform3D;
+using dd4hep::Translation3D;
+using dd4hep::Volume;
+using dd4hep::_toString;
+
+using dd4hep::rec::LayeredCalorimeterData;
+
+
+static Ref_t create_detector(Detector& theDetector, xml_h element, SensitiveDetector sens) {
+ xml_det_t x_det = element;
+ Layering layering(x_det);
+ xml_dim_t dim = x_det.dimensions();
+ string det_name = x_det.nameStr();
+
+ Material air = theDetector.air();
+ Material stavesMaterial = theDetector.material(x_det.materialStr());
+ int numSides = dim.numsides();
+
+ int det_id = x_det.id();
+
+ DetElement sdet(det_name,det_id);
+
+ PlacedVolume pVol;
+
+ // --- create an envelope volume and position it into the world ---------------------
+
+ Volume envelope = dd4hep::xml::createPlacedEnvelope( theDetector, element , sdet ) ;
+
+ sdet.setTypeFlag( DetType::CALORIMETER | DetType::ENDCAP | DetType::HADRONIC ) ;
+
+ if( theDetector.buildType() == BUILD_ENVELOPE ) return sdet ;
+ //-----------------------------------------------------------------------------------
+
+ sens.setType("calorimeter");
+
+ DetElement stave_det("module0stave0",det_id);
+
+ // The way to reaad constant from XML/Detector file.
+ double Hcal_radiator_thickness = theDetector.constant<double>("Hcal_radiator_thickness");
+ double Hcal_endcap_lateral_structure_thickness = theDetector.constant<double>("Hcal_endcap_lateral_structure_thickness");
+ double Hcal_endcap_layer_air_gap = theDetector.constant<double>("Hcal_endcap_layer_air_gap");
+
+ //double Hcal_cells_size = theDetector.constant<double>("Hcal_cells_size");
+ double HcalEndcap_inner_radius = theDetector.constant<double>("Hcal_endcap_inner_radius");
+ double HcalEndcap_outer_radius = theDetector.constant<double>("Hcal_endcap_outer_radius");
+ double HcalEndcap_min_z = theDetector.constant<double>("Hcal_endcap_zmin");
+ double HcalEndcap_max_z = theDetector.constant<double>("Hcal_endcap_zmax");
+
+ double Hcal_steel_cassette_thickness = theDetector.constant<double>("Hcal_steel_cassette_thickness");
+ double HcalServices_outer_FR4_thickness = theDetector.constant<double>("Hcal_services_outer_FR4_thickness");
+ double HcalServices_outer_Cu_thickness = theDetector.constant<double>("Hcal_services_outer_Cu_thickness");
+ double Hcal_endcap_services_module_width = theDetector.constant<double>("Hcal_endcap_services_module_width");
+
+ Material stainless_steel = theDetector.material("stainless_steel");
+ Material PCB = theDetector.material("PCB");
+ Material copper = theDetector.material("Cu");
+
+ std::cout <<"\n HcalEndcap_inner_radius = "
+ <<HcalEndcap_inner_radius/dd4hep::mm <<" mm"
+ <<"\n HcalEndcap_outer_radius = "
+ <<HcalEndcap_outer_radius/dd4hep::mm <<" mm"
+ <<"\n HcalEndcap_min_z = "
+ <<HcalEndcap_min_z/dd4hep::mm <<" mm"
+ <<"\n HcalEndcap_max_z = "
+ <<HcalEndcap_max_z/dd4hep::mm <<" mm"
+ <<std::endl;
+
+ Readout readout = sens.readout();
+ Segmentation seg = readout.segmentation();
+
+
+ BitField64 encoder = seg.decoder();
+ encoder.setValue(0) ;
+
+ // we first have to check whether a multi segmentation is used and then, if so, we
+ // will check the <subsegmentation key="" value=""/> element, for which subsegmentation to use for filling the
+ // DDRec:LayeredCalorimeterData information.
+
+ // check if we have a multi segmentation :
+ dd4hep::DDSegmentation::MultiSegmentation* multiSeg =
+ dynamic_cast< dd4hep::DDSegmentation::MultiSegmentation*>( seg.segmentation() ) ;
+
+ int sensitive_slice_number = -1 ;
+
+ if( multiSeg ){
+
+ try{
+ // check if we have an entry for the subsegmentation to be used
+ xml_comp_t segxml = x_det.child( _Unicode( subsegmentation ) ) ;
+
+ std::string keyStr = segxml.attr<std::string>( _Unicode(key) ) ;
+ int keyVal = segxml.attr<int>( _Unicode(value) ) ;
+
+ encoder[ keyStr ] = keyVal ;
+
+ // if we have a multisegmentation that uses the slice as key, we need to know for the
+ // computation of the layer parameters in LayeredCalorimeterData::Layer below
+ if( keyStr == "slice" ){
+ sensitive_slice_number = keyVal ;
+ }
+
+ } catch(const std::runtime_error &) {
+ throw lcgeo::GeometryException( "SHcalSc04_Endcaps_v01: Error: MultiSegmentation specified but no "
+ " <subsegmentation key="" value=""/> element defined for detector ! " ) ;
+ }
+ }
+
+ //========== fill data for reconstruction ============================
+ LayeredCalorimeterData* caloData = new LayeredCalorimeterData ;
+ caloData->layoutType = LayeredCalorimeterData::EndcapLayout ;
+ caloData->inner_symmetry = 4 ; // hard code cernter box hole
+ caloData->outer_symmetry = 0 ; // outer tube, or 8 for Octagun
+ caloData->phi0 = 0 ;
+
+ /// extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm.
+ caloData->extent[0] = HcalEndcap_inner_radius ;
+ caloData->extent[1] = HcalEndcap_outer_radius ;
+ caloData->extent[2] = HcalEndcap_min_z ;
+ caloData->extent[3] = HcalEndcap_max_z ;
+
+
+ int endcapID = 0;
+ for(xml_coll_t c(x_det.child(_U(dimensions)),_U(dimensions)); c; ++c)
+ {
+ xml_comp_t l(c);
+
+ double dim_x = l.attr<double>(_Unicode(dim_x));
+ double dim_y = l.attr<double>(_Unicode(dim_y));
+ double dim_z = l.attr<double>(_Unicode(dim_z));
+ double pos_y = l.attr<double>(_Unicode(y_offset));
+
+ // Hcal Endcap module shape
+ double box_half_x= dim_x/2.0; // module width, all are same
+ double box_half_y= dim_y/2.0; // module length, changing
+ double box_half_z= dim_z/2.0; // total thickness, all are same
+
+ double x_offset = box_half_x*numSides-box_half_x*endcapID*2.0-box_half_x;
+ double y_offset = pos_y;
+
+ Box EndcapModule(box_half_x,box_half_y,box_half_z);
+
+ // define the name of each endcap Module
+ string envelopeVol_name = det_name+_toString(endcapID,"_EndcapModule%d");
+
+ Volume envelopeVol(envelopeVol_name,EndcapModule,stavesMaterial);
+
+ // Set envelope volume attributes.
+ envelopeVol.setAttributes(theDetector,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
+
+
+ double FEE_half_x = box_half_x-Hcal_endcap_services_module_width/2.0;
+ double FEE_half_y = Hcal_endcap_services_module_width/2.0;
+ double FEE_half_z = box_half_z;
+
+ Box FEEBox(FEE_half_x,FEE_half_y,FEE_half_z);
+ Volume FEEModule("Hcal_endcap_FEE",FEEBox,air);
+
+ double FEELayer_thickness = Hcal_steel_cassette_thickness + HcalServices_outer_FR4_thickness + HcalServices_outer_Cu_thickness;
+ Box FEELayerBox(FEE_half_x,FEE_half_y,FEELayer_thickness/2.0);
+ Volume FEELayer("FEELayer",FEELayerBox,air);
+
+ Box FEELayerSteelBox(FEE_half_x,FEE_half_y,Hcal_steel_cassette_thickness/2.0);
+ Volume FEELayerSteel("FEELayerSteel",FEELayerSteelBox,stainless_steel);
+ pVol = FEELayer.placeVolume(FEELayerSteel,
+ Position(0,
+ 0,
+ (-FEELayer_thickness/2.0
+ +Hcal_steel_cassette_thickness/2.0)));
+
+ Box FEELayerFR4Box(FEE_half_x,FEE_half_y,HcalServices_outer_FR4_thickness/2.0);
+ Volume FEELayerFR4("FEELayerFR4",FEELayerFR4Box,PCB);
+ pVol = FEELayer.placeVolume(FEELayerFR4,
+ Position(0,
+ 0,
+ (-FEELayer_thickness/2.0+Hcal_steel_cassette_thickness
+ +HcalServices_outer_FR4_thickness/2.0)));
+
+ Box FEELayerCuBox(FEE_half_x,FEE_half_y,HcalServices_outer_Cu_thickness/2.0);
+ Volume FEELayerCu("FEELayerCu",FEELayerCuBox,copper);
+ pVol = FEELayer.placeVolume(FEELayerCu,
+ Position(0,
+ 0,
+ (-FEELayer_thickness/2.0+Hcal_steel_cassette_thickness+HcalServices_outer_FR4_thickness +HcalServices_outer_Cu_thickness/2.0)));
+
+
+ // ========= Create Hcal Chamber (i.e. Layers) ==============================
+ // It will be the sub volume for placing the slices.
+ // Itself will be placed into the Hcal Endcap modules envelope.
+ // ==========================================================================
+
+ // create Layer (air) and place the slices (Polystyrene,Cu,FR4,air) into it.
+ // place the Layer into the Hcal Endcap Modules envelope (stavesMaterial).
+
+ // First Hcal Chamber position, start after first radiator
+ double layer_pos_z = - box_half_z + Hcal_radiator_thickness;
+
+ // Create Hcal Endcap Chamber without radiator
+ // Place into the Hcal Encap module envelope, after each radiator
+ int layer_num = 1;
+ for(xml_coll_t m(x_det,_U(layer)); m; ++m) {
+ xml_comp_t x_layer = m;
+ int repeat = x_layer.repeat(); // Get number of layers.
+
+ double layer_thickness = layering.layer(layer_num)->thickness();
+ string layer_name = envelopeVol_name+"_layer";
+ DetElement layer(stave_det,layer_name,det_id);
+
+ // Active Layer box & volume
+ double active_layer_dim_x = box_half_x - Hcal_endcap_lateral_structure_thickness - Hcal_endcap_layer_air_gap;
+ double active_layer_dim_y = box_half_y;
+ double active_layer_dim_z = layer_thickness/2.0;
+
+ // Build chamber including air gap
+ // The Layer will be filled with slices,
+ Volume layer_vol(layer_name, Box((active_layer_dim_x + Hcal_endcap_layer_air_gap),
+ active_layer_dim_y,active_layer_dim_z), air);
+
+
+
+ encoder["layer"] = layer_num ;
+ std::vector<double> cellSizeVector = seg.segmentation()->cellDimensions( encoder.getValue() );
+
+ LayeredCalorimeterData::Layer caloLayer ;
+ caloLayer.cellSize0 = cellSizeVector[0];
+ caloLayer.cellSize1 = cellSizeVector[1];
+
+ // ========= Create sublayer slices =========================================
+ // Create and place the slices into Layer
+ // ==========================================================================
+
+ // Create the slices (sublayers) within the Hcal Chamber.
+ double slice_pos_z = -(layer_thickness / 2.0);
+ int slice_number = 0;
+
+ double nRadiationLengths=0.;
+ double nInteractionLengths=0.;
+ double thickness_sum=0;
+
+ nRadiationLengths = Hcal_radiator_thickness/(stavesMaterial.radLength());
+ nInteractionLengths = Hcal_radiator_thickness/(stavesMaterial.intLength());
+ thickness_sum = Hcal_radiator_thickness;
+
+ for(xml_coll_t k(x_layer,_U(slice)); k; ++k) {
+ xml_comp_t x_slice = k;
+ string slice_name = layer_name + _toString(slice_number,"_slice%d");
+ double slice_thickness = x_slice.thickness();
+ Material slice_material = theDetector.material(x_slice.materialStr());
+ DetElement slice(layer,_toString(slice_number,"slice%d"),det_id);
+
+ slice_pos_z += slice_thickness / 2.0;
+
+ // Slice volume & box
+ Volume slice_vol(slice_name,Box(active_layer_dim_x,active_layer_dim_y,slice_thickness/2.0),slice_material);
+
+ nRadiationLengths += slice_thickness/(2.*slice_material.radLength());
+ nInteractionLengths += slice_thickness/(2.*slice_material.intLength());
+ thickness_sum += slice_thickness/2;
+
+
+ if ( x_slice.isSensitive() ) {
+
+ slice_vol.setSensitiveDetector(sens);
+
+ // if we have a multisegmentation based on slices, we need to use the correct slice here
+ if ( sensitive_slice_number<0 || sensitive_slice_number == slice_number ) {
+
+ //Store "inner" quantities
+ caloLayer.inner_nRadiationLengths = nRadiationLengths;
+ caloLayer.inner_nInteractionLengths = nInteractionLengths;
+ caloLayer.inner_thickness = thickness_sum;
+ //Store scintillator thickness
+ caloLayer.sensitive_thickness = slice_thickness;
+
+ //Reset counters to measure "outside" quantitites
+ nRadiationLengths=0.;
+ nInteractionLengths=0.;
+ thickness_sum = 0.;
+ }
+ }
+
+ nRadiationLengths += slice_thickness/(2.*slice_material.radLength());
+ nInteractionLengths += slice_thickness/(2.*slice_material.intLength());
+ thickness_sum += slice_thickness/2;
+
+ // Set region, limitset, and vis.
+ slice_vol.setAttributes(theDetector,x_slice.regionStr(),x_slice.limitsStr(),x_slice.visStr());
+ // slice PlacedVolume
+ PlacedVolume slice_phv = layer_vol.placeVolume(slice_vol,Position(0,0,slice_pos_z));
+ slice_phv.addPhysVolID("slice",slice_number);
+
+ slice.setPlacement(slice_phv);
+ // Increment Z position for next slice.
+ slice_pos_z += slice_thickness / 2.0;
+ // Increment slice number.
+ ++slice_number;
+ }
+ // Set region, limitset, and vis.
+ layer_vol.setAttributes(theDetector,x_layer.regionStr(),x_layer.limitsStr(),x_layer.visStr());
+
+
+ //Store "outer" quantities
+ caloLayer.outer_nRadiationLengths = nRadiationLengths;
+ caloLayer.outer_nInteractionLengths = nInteractionLengths;
+ caloLayer.outer_thickness = thickness_sum;
+
+ // ========= Place the Layer (i.e. Chamber) =================================
+ // Place the Layer into the Hcal Endcap module envelope.
+ // with the right position and rotation.
+ // Registry the IDs (layer, stave, module).
+ // Place the same layer 48 times into Endcap module
+ // ==========================================================================
+
+ for (int j = 0; j < repeat; j++) {
+
+ // Layer position in y within the Endcap Modules.
+ layer_pos_z += layer_thickness / 2.0;
+
+ PlacedVolume layer_phv = envelopeVol.placeVolume(layer_vol,
+ Position(0,0,layer_pos_z));
+ // registry the ID of Layer, stave and module
+ layer_phv.addPhysVolID("layer",layer_num);
+
+ // then setPlacement for it.
+ layer.setPlacement(layer_phv);
+
+ pVol = FEEModule.placeVolume(FEELayer,
+ Position(0,0,layer_pos_z));
+ //-----------------------------------------------------------------------------------------
+ if ( caloData->layers.size() < (unsigned int)repeat ) {
+
+ caloLayer.distance = HcalEndcap_min_z + box_half_z + layer_pos_z
+ - caloLayer.inner_thickness ; // Will be added later at "DDMarlinPandora/DDGeometryCreator.cc:226" to get center of sensitive element
+ caloLayer.absorberThickness = Hcal_radiator_thickness ;
+
+ caloData->layers.push_back( caloLayer ) ;
+ }
+ //-----------------------------------------------------------------------------------------
+
+
+ // ===== Prepare for next layer (i.e. next Chamber) =========================
+ // Prepare the chamber placement position and the chamber dimension
+ // ==========================================================================
+
+ // Increment the layer_pos_y
+ // Place Hcal Chamber after each radiator
+ layer_pos_z += layer_thickness / 2.0;
+ layer_pos_z += Hcal_radiator_thickness;
+ ++layer_num;
+ }
+
+
+ }
+
+
+ // =========== Place Hcal Endcap envelope ===================================
+ // Finally place the Hcal Endcap envelope into the world volume.
+ // Registry the stave(up/down), module(left/right) and endcapID.
+ // ==========================================================================
+
+ // Acording to the number of staves and modules,
+ // Place the same Hcal Endcap module volume into the world volume
+ // with the right position and rotation.
+ for(int stave_num=0;stave_num<2;stave_num++){
+
+ double EndcapModule_pos_x = 0;
+ double EndcapModule_pos_y = 0;
+ double EndcapModule_pos_z = 0;
+ double rot_EM = 0;
+
+ double EndcapModule_center_pos_z = HcalEndcap_min_z + box_half_z;
+
+ double FEEModule_pos_x = 0;
+ double FEEModule_pos_y = 0;
+ double FEEModule_pos_z = 0;
+ double FEEModule_center_pos_z = HcalEndcap_min_z + box_half_z;
+
+ switch (stave_num)
+ {
+ case 0 :
+ EndcapModule_pos_x = x_offset;
+ EndcapModule_pos_y = y_offset;
+ FEEModule_pos_x = x_offset;
+ FEEModule_pos_y = y_offset + box_half_y + Hcal_endcap_services_module_width/2.0;
+ break;
+ case 1 :
+ EndcapModule_pos_x = -x_offset;
+ EndcapModule_pos_y = -y_offset;
+ FEEModule_pos_x = -x_offset;
+ FEEModule_pos_y = -y_offset - box_half_y - Hcal_endcap_services_module_width/2.0;
+ break;
+ }
+
+ for(int module_num=0;module_num<2;module_num++) {
+
+ int module_id = (module_num==0)? 0:6;
+
+ rot_EM = (module_id==0)?M_PI:0;
+
+ EndcapModule_pos_z = (module_id==0)? -EndcapModule_center_pos_z:EndcapModule_center_pos_z;
+
+ PlacedVolume env_phv = envelope.placeVolume(envelopeVol,
+ Transform3D(RotationX(rot_EM),
+ Translation3D(EndcapModule_pos_x,
+ EndcapModule_pos_y,
+ EndcapModule_pos_z)));
+ env_phv.addPhysVolID("tower",endcapID);
+ env_phv.addPhysVolID("stave",stave_num); // y: up /down
+ env_phv.addPhysVolID("module",module_id); // z: -/+ 0/6
+ env_phv.addPhysVolID("system",det_id);
+
+ FEEModule_pos_z = (module_id==0)? -FEEModule_center_pos_z:FEEModule_center_pos_z;
+
+ if (!(endcapID==0))
+ env_phv = envelope.placeVolume(FEEModule,
+ Transform3D(RotationX(rot_EM),
+ Translation3D(FEEModule_pos_x,
+ FEEModule_pos_y,
+ FEEModule_pos_z)));
+
+
+ DetElement sd = (module_num==0&&stave_num==0) ? stave_det : stave_det.clone(_toString(module_id,"module%d")+_toString(stave_num,"stave%d"));
+ sd.setPlacement(env_phv);
+
+ }
+
+ }
+
+ endcapID++;
+
+ }
+
+ sdet.addExtension< LayeredCalorimeterData >( caloData ) ;
+
+ return sdet;
+}
+
+
+
+
+DECLARE_DETELEMENT(SHcalSc04_Endcaps_v01, create_detector)
diff --git a/Detector/DetCRD/compact/CRD_common_v01/SHcalSc04_Barrel_v04_01.xml b/Detector/DetCRD/compact/CRD_common_v01/SHcalSc04_Barrel_v04_01.xml
new file mode 100644
index 0000000000000000000000000000000000000000..32810ce1c0da9fdf6f602daa770d2815778ddb37
--- /dev/null
+++ b/Detector/DetCRD/compact/CRD_common_v01/SHcalSc04_Barrel_v04_01.xml
@@ -0,0 +1,71 @@
+<!-- comment>Calorimeters</comment -->
+
+<lccdd>
+ <define>
+ <constant name="Hcal_cell_size" value="10*mm"/>
+ <constant name="Hcal_inner_radius" value="Hcal_barrel_inner_radius"/>
+ <constant name="Hcal_half_length" value="Hcal_barrel_half_length"/>
+ <constant name="Hcal_inner_symmetry" value="Hcal_barrel_symmetry"/>
+ <constant name="Hcal_nlayers" value="38"/>
+ <constant name="Hcal_radiator_thickness" value="20.0*mm"/>
+ <constant name="Hcal_chamber_thickness" value="6.5*mm"/>
+ <constant name="Hcal_back_plate_thickness" value="15*mm"/>
+ <constant name="Hcal_lateral_structure_thickness" value="10*mm"/>
+ <constant name="Hcal_stave_gaps" value="0*mm"/>
+ <constant name="Hcal_middle_stave_gaps" value="0*mm"/>
+ <constant name="Hcal_modules_gap" value="2*mm"/>
+ <constant name="Hcal_layer_air_gap" value="0*mm"/>
+ <constant name="HcalSD_glass_anode_thickness" value="0.7*mm"/>
+ <constant name="HcalSD_sensitive_gas_gap" value="1.2*mm"/>
+ <constant name="HcalSD_glass_cathode_thickness" value="1.1*mm"/>
+ <constant name="Hcal_scintillator_thickness" value="3.0*mm"/>
+ <constant name="Ecal_outer_radius" value="Ecal_barrel_outer_radius"/>
+ <constant name="Hcal_readout_segmentation_slice" value="3"/>
+ </define>
+ <detectors>
+ <detector name="HcalBarrel" type="SHcalSc04_Barrel_v04" id="DetID_HCAL" readout="HcalBarrelCollection" vis="GreenVis" insideTrackingVolume="false" >
+ <comment>Hadron Calorimeter Barrel</comment>
+
+ <envelope vis="ILD_HCALVis">
+ <shape type="BooleanShape" operation="Subtraction" material="Air" >
+ <shape type="Cone" rmin1="0.0" rmax1="Hcal_outer_radius + env_safety" rmin2="0.0" rmax2="Hcal_outer_radius + env_safety" z="Hcal_half_length + env_safety/2.0"/>
+ <shape type="PolyhedraRegular" numsides="Hcal_inner_symmetry" rmin="0.0"
+ rmax="Hcal_inner_radius - env_safety" dz="2*(Hcal_half_length + env_safety)"/>
+ </shape>
+ <rotation x="0" y="0" z="90*deg-180*deg/Hcal_inner_symmetry"/>
+ </envelope>
+
+ <type_flags type=" DetType_CALORIMETER + DetType_BARREL + DetType_HADRONIC " />
+
+ <staves material = "Steel235" vis="BlueVis"/>
+
+
+ <!-- select which subsegmentation will be used to fill the DDRec:LayeredCalorimeterData cell dimensions -->
+ <subsegmentation key="slice" value="Hcal_readout_segmentation_slice"/>
+
+ <layer repeat="Hcal_nlayers" vis="SeeThrough">
+ <slice material="FloatGlass" thickness="HcalSD_glass_anode_thickness" vis="Invisible"/>
+ <slice material="RPCGAS2" thickness="HcalSD_sensitive_gas_gap" sensitive="yes" limits="cal_limits" vis="YellowVis"/>
+ <slice material="FloatGlass" thickness="HcalSD_glass_cathode_thickness" vis="Invisible"/>
+ <slice material="G4_POLYSTYRENE" thickness = "Hcal_scintillator_thickness" sensitive = "yes" limits="cal_limits" vis="CyanVis" />
+ <slice material="Air" thickness="Hcal_chamber_thickness - ( HcalSD_glass_anode_thickness + HcalSD_sensitive_gas_gap + HcalSD_glass_cathode_thickness + Hcal_scintillator_thickness)" vis="Invisible" />
+ </layer>
+ </detector>
+ </detectors>
+
+ <readouts>
+ <readout name="HcalBarrelCollection">
+ <segmentation type="MultiSegmentation" key="slice">
+ <segmentation name="RPCgrid" type="CartesianGridXY" key_value="1" grid_size_x="Hcal_cell_size" grid_size_y="Hcal_cell_size" />
+ <segmentation name="Scigrid" type="TiledLayerGridXY" key_value="3" grid_size_x="3" grid_size_y="3.03248"/>
+ </segmentation>
+ <hits_collections>
+ <hits_collection name="HCalBarrelRPCHits" key="slice" key_value="1"/>
+ <hits_collection name="HcalBarrelRegCollection" key="slice" key_value="3"/>
+ </hits_collections>
+ <id>system:5,module:3,stave:4,tower:5,layer:6,slice:4,x:32:-16,y:-16</id>
+ </readout>
+ </readouts>
+
+
+</lccdd>
diff --git a/Detector/DetCRD/compact/CRD_common_v01/SHcalSc04_Endcaps_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/SHcalSc04_Endcaps_v01_01.xml
new file mode 100644
index 0000000000000000000000000000000000000000..6b11b773d72d530ae92aaf640de0b6e8d705d168
--- /dev/null
+++ b/Detector/DetCRD/compact/CRD_common_v01/SHcalSc04_Endcaps_v01_01.xml
@@ -0,0 +1,90 @@
+<lccdd>
+ <define>
+ <constant name="SDHCal_cell_size" value="10*mm"/>
+ <constant name="AHCal_cell_size" value="10*mm"/>
+ <constant name="Hcal_endcap_lateral_structure_thickness" value="5.0*mm"/>
+ <constant name="Hcal_endcap_layer_air_gap" value="2.5*mm"/>
+ <constant name="Hcal_steel_cassette_thickness" value="0.5*mm"/>
+ <constant name="Hcal_services_outer_FR4_thickness" value="2.8*mm"/>
+ <constant name="Hcal_services_outer_Cu_thickness" value="0.4*mm"/>
+ <constant name="Hcal_endcap_services_module_width" value="100.0*mm"/>
+ <constant name="Hcal_endcap_nlayers" value="40"/>
+ <constant name="Hcal_endcap_env_thickness" value="Hcal_endcap_zmax-Hcal_endcap_zmin"/>
+ <constant name="Hcal_x_modul" value="12"/>
+ <constant name="Hcal_x_width" value="Hcal_endcap_inner_radius*2/2"/><!--350-->
+ <constant name="Hcal_y_height" value="Hcal_x_width*Hcal_x_modul/2"/><!--2100-->
+ <constant name="Hcal_hole_height" value="Hcal_endcap_inner_radius"/>
+ <constant name="Hcal_r_max" value="Hcal_y_height/cos(pi/Hcal_endcap_symmetry)"/><!--2174.08-->
+ <constant name="Hcal_x_top" value="Hcal_y_height*tan(pi/Hcal_endcap_symmetry)"/> <!--562.69-->
+ <constant name="Hcal_x_point" value="(Hcal_y_height+Hcal_x_top*tan(2*pi/Hcal_endcap_symmetry))/(1+tan(2*pi/Hcal_endcap_symmetry))"/><!--1537.30-->
+ <constant name="Hcal_y5" value="Hcal_y_height-(Hcal_x_width*2-Hcal_x_top)*tan(2*pi/Hcal_endcap_symmetry)"/>
+ <constant name="Hcal_y4" value="Hcal_y5-Hcal_x_width*tan(2*pi/Hcal_endcap_symmetry)"/>
+ <constant name="Hcal_y3" value="Hcal_y4-Hcal_x_width*tan(2*pi/Hcal_endcap_symmetry)"/>
+ <constant name="Hcal_y2" value="Hcal_x_top+Hcal_x_width/tan(2*pi/Hcal_endcap_symmetry)"/>
+
+ </define>
+ <detectors>
+ <detector id="DetID_HCAL_ENDCAP" name="HcalEndcap" type="SHcalSc04_Endcaps_v01" readout="HcalEndcapsReadout" vis="GreenVis" calorimeterType="HAD_ENDCAP">
+ <comment>Hadron Calorimeter Endcap</comment>
+
+ <envelope vis="ILD_HCALVis">
+ <shape type="BooleanShape" operation="Subtraction" material="Air"><!--2. create center box hole -->
+ <shape type="BooleanShape" operation="Subtraction" material="Air"><!--1. create Endcaps envelope -->
+ <shape type="Tube" rmin="0.0" rmax="Solenoid_inner_radius" dz="Hcal_endcap_zmax + env_safety"/>
+ <shape type="Tube" rmin="0.0" rmax="Solenoid_inner_radius" dz="Hcal_endcap_zmin - env_safety"/>
+ </shape>
+ <shape type="Box" dx="Hcal_endcap_inner_radius - env_safety" dy="Hcal_endcap_inner_radius - env_safety"
+ dz="Hcal_endcap_zmax + 2.0*env_safety"/>
+ </shape>
+ <rotation x="0" y="0" z="0"/>
+ </envelope>
+
+ <type_flags type=" DetType_CALORIMETER + DetType_ENDCAP + DetType_HADRONIC " />
+
+ <material name="Steel235"/><!-- radiator and the thickness has been defined in the main xml file-->
+
+ <dimensions numsides="Hcal_x_modul" >
+ <dimensions id="1" y_offset="Hcal_x_top/2" dim_x="Hcal_hole_height" dim_y="Hcal_x_top" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="2" y_offset="Hcal_y2/2" dim_x="Hcal_hole_height" dim_y="Hcal_y2" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="3" y_offset="Hcal_y3/2" dim_x="Hcal_hole_height" dim_y="Hcal_y3" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="4" y_offset="Hcal_y4/2" dim_x="Hcal_hole_height" dim_y="Hcal_y4" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="5" y_offset="Hcal_y5/2" dim_x="Hcal_hole_height" dim_y="Hcal_y5" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="6" y_offset="Hcal_y_height/2+Hcal_hole_height/2" dim_x="Hcal_hole_height" dim_y="Hcal_y_height-Hcal_hole_height" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="7" y_offset="Hcal_y_height/2+Hcal_hole_height/2" dim_x="Hcal_hole_height" dim_y="Hcal_y_height-Hcal_hole_height" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="8" y_offset="Hcal_y5/2" dim_x="Hcal_hole_height" dim_y="Hcal_y5" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="9" y_offset="Hcal_y4/2" dim_x="Hcal_hole_height" dim_y="Hcal_y4" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="10" y_offset="Hcal_y3/2" dim_x="Hcal_hole_height" dim_y="Hcal_y3" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="11" y_offset="Hcal_y2/2" dim_x="Hcal_hole_height" dim_y="Hcal_y2" dim_z="Hcal_endcap_env_thickness"/>
+ <dimensions id="12" y_offset="Hcal_x_top/2" dim_x="Hcal_hole_height" dim_y="Hcal_x_top" dim_z="Hcal_endcap_env_thickness"/>
+ </dimensions>
+
+ <!-- select which subsegmentation will be used to fill the DDRec:LayeredCalorimeterData cell dimensions -->
+ <subsegmentation key="slice" value="Hcal_readout_segmentation_slice"/>
+
+ <layer repeat="Hcal_endcap_nlayers" vis="SeeThrough">
+ <slice material="FloatGlass" thickness="HcalSD_glass_anode_thickness" vis="Invisible"/>
+ <slice material="RPCGAS2" thickness="HcalSD_sensitive_gas_gap" sensitive="yes" limits="cal_limits" vis="YellowVis"/>
+ <slice material="FloatGlass" thickness="HcalSD_glass_cathode_thickness" vis="Invisible"/>
+ <slice material="G4_POLYSTYRENE" thickness = "Hcal_scintillator_thickness" sensitive = "yes" limits="cal_limits" vis="CyanVis" />
+ <slice material="Air" thickness="Hcal_chamber_thickness - ( HcalSD_glass_anode_thickness + HcalSD_sensitive_gas_gap + HcalSD_glass_cathode_thickness + Hcal_scintillator_thickness)" vis="Invisible" />
+ </layer>
+ </detector>
+ </detectors>
+
+ <readouts>
+ <readout name="HcalEndcapsReadout">
+ <segmentation type="MultiSegmentation" key="slice">
+ <segmentation name="RPCgrid" type="CartesianGridXY" key_value="1" grid_size_x="SDHCal_cell_size" grid_size_y="SDHCal_cell_size" offset_x="SDHCal_cell_size/2.0" offset_y="SDHCal_cell_size/2.0" />
+ <segmentation name="Scigrid" type="CartesianGridXY" key_value="3" grid_size_x="AHCal_cell_size" grid_size_y="AHCal_cell_size" offset_x="AHCal_cell_size/2.0" offset_y="AHCal_cell_size/2.0" />
+ </segmentation>
+ <hits_collections>
+ <hits_collection name="HCalEndcapRPCHits" key="slice" key_value="1"/>
+ <hits_collection name="HcalEndcapsCollection" key="slice" key_value="3"/>
+ </hits_collections>
+ <id>system:5,module:3,stave:3,tower:5,layer:6,slice:4,x:32:-16,y:-16</id>
+ </readout>
+ </readouts>
+
+
+</lccdd>
+
diff --git a/Detector/DetCRD/compact/Standalone/Standalone-SHcalSc04.xml b/Detector/DetCRD/compact/Standalone/Standalone-SHcalSc04.xml
new file mode 100644
index 0000000000000000000000000000000000000000..42edeebec1c972c0eb4cfa5665c02759a332e711
--- /dev/null
+++ b/Detector/DetCRD/compact/Standalone/Standalone-SHcalSc04.xml
@@ -0,0 +1,51 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+ xmlns:xs="http://www.w3.org/2001/XMLSchema"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+ <info name="StandaloneEcalRotCrystal"
+ title="CepC standalone calorimeter with rotated crystal"
+ author="C.D.Fu"
+ url="http://cepc.ihep.ac.cn"
+ status="developing"
+ version="v01">
+ <comment>CepC detector simulation models used for detector study </comment>
+ </info>
+
+ <includes>
+ <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/>
+ <gdmlFile ref="../CRD_common_v01/materials.xml"/>
+ </includes>
+
+ <define>
+ <constant name="world_size" value="6*m"/>
+ <constant name="world_x" value="world_size"/>
+ <constant name="world_y" value="world_size"/>
+ <constant name="world_z" value="world_size"/>
+
+ <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/>
+ </define>
+
+ <include ref="./Dimensions_v01_01.xml"/>
+
+ <!--include ref="../CRD_common_v01/Coil_Simple_v01_01.xml"/-->
+ <include ref="../CRD_common_v01/SHcalSc04_Barrel_v04_01.xml"/>
+ <include ref="../CRD_common_v01/SHcalSc04_Endcaps_v01_01.xml"/>
+
+ <fields>
+ <field name="InnerSolenoid" type="solenoid"
+ inner_field="Field_nominal_value"
+ outer_field="0"
+ zmax="SolenoidCoil_half_length"
+ inner_radius="SolenoidCoil_center_radius"
+ outer_radius="Solenoid_outer_radius">
+ </field>
+ <field name="OuterSolenoid" type="solenoid"
+ inner_field="0"
+ outer_field="Field_outer_nominal_value"
+ zmax="SolenoidCoil_half_length"
+ inner_radius="Solenoid_outer_radius"
+ outer_radius="Yoke_barrel_inner_radius">
+ </field>
+ </fields>
+
+</lccdd>
diff --git a/Detector/DetCRD/scripts/CRD-Sim.py b/Detector/DetCRD/scripts/CRD-Sim.py
index a4ad74c40c73b4b30ba43e51560e7d755a9b6d77..0e734b785f0c6442b320bc75e43c3d10472bc8a5 100644
--- a/Detector/DetCRD/scripts/CRD-Sim.py
+++ b/Detector/DetCRD/scripts/CRD-Sim.py
@@ -89,6 +89,23 @@ detsimalg.AnaElems = [
]
detsimalg.RootDetElem = "WorldDetElemTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+from Configurables import DriftChamberSensDetTool
+dc_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+dc_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+ dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+ dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+ dedx_simtool.material_Z = 2
+ dedx_simtool.material_A = 4
+ dedx_simtool.scale = 10
+ dedx_simtool.resolution = 0.0001
+
# output
from Configurables import PodioOutput
out = PodioOutput("outputalg")
diff --git a/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py b/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py
index 1810a5749abca5cab4561466e90b32ddb05ef4ab..19836013be8fcc3b1dcbcd6993473de64b1a0ba5 100644
--- a/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py
+++ b/Detector/DetCRD/scripts/CRD_o1_v01-SimRec.py
@@ -84,6 +84,23 @@ detsimalg.AnaElems = [
]
detsimalg.RootDetElem = "WorldDetElemTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+from Configurables import DriftChamberSensDetTool
+dc_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+dc_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+ dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+ dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+ dedx_simtool.material_Z = 2
+ dedx_simtool.material_A = 4
+ dedx_simtool.scale = 10
+ dedx_simtool.resolution = 0.0001
+
from Configurables import MarlinEvtSeeder
evtseeder = MarlinEvtSeeder("EventSeeder")
diff --git a/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py b/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py
index acee999de77be93b4b760ebbd2e2e23c3a91698d..253ea21a71125f20f9216094dd2071e2fdd950e4 100644
--- a/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py
+++ b/Detector/DetCRD/scripts/CRD_o1_v02-SimRec.py
@@ -84,6 +84,23 @@ detsimalg.AnaElems = [
]
detsimalg.RootDetElem = "WorldDetElemTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+from Configurables import DriftChamberSensDetTool
+dc_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+dc_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+ dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+ dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+ dedx_simtool.material_Z = 2
+ dedx_simtool.material_A = 4
+ dedx_simtool.scale = 10
+ dedx_simtool.resolution = 0.0001
+
from Configurables import MarlinEvtSeeder
evtseeder = MarlinEvtSeeder("EventSeeder")
diff --git a/Detector/DetInterface/include/DetInterface/IGeomSvc.h b/Detector/DetInterface/include/DetInterface/IGeomSvc.h
index d4f508b99779f44e6b30811f4296ac722d8cf9c8..51c26406a307d3e66724d207a730bbeca823c216 100644
--- a/Detector/DetInterface/include/DetInterface/IGeomSvc.h
+++ b/Detector/DetInterface/include/DetInterface/IGeomSvc.h
@@ -27,11 +27,6 @@ namespace dd4hep {
class StructExtension;
-namespace gear{
- class ZPlanarParametersImpl;
- class GearParametersImpl;
-}
-class TMaterial;
// class G4VUserDetectorConstruction;
class GAUDI_API IGeomSvc : virtual public IService {
@@ -49,16 +44,6 @@ public:
// short cut to retrieve the Decoder according to the Readout name
virtual Decoder* getDecoder(const std::string& readout_name) = 0;
- // obsolete parameter format, will remove once StructExtension<> validated
- virtual const gear::ZPlanarParametersImpl* getVXDParameters() = 0;
-
- virtual const dd4hep::rec::ZPlanarData* getVXDData() = 0;
- virtual const dd4hep::rec::ConicalSupportData* getBeamPipeData() =0;
-
- virtual const std::map<std::string,double>& getDetParameters(std::string s) = 0;
- virtual double getDetParameter(std::string set_name, std::string par_name) = 0;
- virtual TMaterial* getMaterial(std::string s) = 0;
-
virtual ~IGeomSvc() {}
};
diff --git a/Detector/GeomSvc/src/GeomSvc.cpp b/Detector/GeomSvc/src/GeomSvc.cpp
index 54736e6c617c0e4c2591fd14817641751fe67f1e..f7ae4a2bd6c04f6f728e51a8c251e1cd9a31575e 100644
--- a/Detector/GeomSvc/src/GeomSvc.cpp
+++ b/Detector/GeomSvc/src/GeomSvc.cpp
@@ -22,7 +22,7 @@
DECLARE_COMPONENT(GeomSvc)
GeomSvc::GeomSvc(const std::string& name, ISvcLocator* svc)
-: base_class(name, svc), m_dd4hep_geo(nullptr), m_vxdData(nullptr), m_beamPipeData(nullptr){
+: base_class(name, svc), m_dd4hep_geo(nullptr){
}
@@ -38,39 +38,12 @@ GeomSvc::initialize() {
// if failed to load the compact, a runtime error will be thrown.
m_dd4hep_geo->fromCompact(m_dd4hep_xmls.value());
- dd4hep::DetElement world = m_dd4hep_geo->world();
- //info() << world.type() << " " << world.path() << " " << world.placementPath() << endmsg;
- const std::map<std::string, dd4hep::DetElement>& subs = world.children();
- for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){
- dd4hep::DetElement sub = it->second;
- info() << it->first << " " << sub.path() << " " << sub.placementPath() << endmsg;
- if(it->first=="Tube"){
- dd4hep::Volume vol = sub.volume();
- dd4hep::Solid solid = vol.solid();
- //info() << " " << solid.type() << " " << solid.name() << endmsg;
- const std::map<std::string, dd4hep::DetElement>& pipes = sub.children();
- for(std::map<std::string, dd4hep::DetElement>::const_iterator it=pipes.begin();it!=pipes.end();it++){
- dd4hep::DetElement pipe = it->second;
- //info() << " " << it->first << " " << pipe.id() << " " << pipe.path() << " " << pipe.placementPath() << endmsg;
- }
- try{
- m_beamPipeData = sub.extension<dd4hep::rec::ConicalSupportData>();
- }
- catch(std::runtime_error& e){
- info() << e.what() << " " << m_beamPipeData << endmsg;
- }
- }
- if(it->first=="VXD"){
- sc = convertVXD(sub);
- }
- }
return sc;
}
StatusCode
GeomSvc::finalize() {
StatusCode sc;
- if(m_vxdParameters) delete m_vxdParameters;
return sc;
}
@@ -119,484 +92,3 @@ GeomSvc::getDecoder(const std::string& readout_name) {
return decoder;
}
-
-
-const std::map<std::string,double>& GeomSvc::getDetParameters(std::string name){
- if(m_detParameters.find(name)!=m_detParameters.end()) return m_detParameters[name];
- else{
- char message[200];
- sprintf(message,"GeomSvc has not the parameter set named %s", name.c_str());
- throw std::runtime_error(message);
- }
-}
-
-double GeomSvc::getDetParameter(std::string set_name, std::string par_name){
- std::map<std::string, std::map<std::string,double> >::iterator it=m_detParameters.find(set_name);
- if(it!=m_detParameters.end()){
- if(it->second.find(par_name)!=it->second.end()) return it->second[par_name];
- }
- char message[200];
- sprintf(message,"GeomSvc has not the parameter named %s in set %s", par_name.c_str(), set_name.c_str());
- throw std::runtime_error(message);
-}
-
-StatusCode GeomSvc::convertVXD(dd4hep::DetElement& vxd){
- StatusCode sc;
- //fucd: another method to obtain parameters, but not fully for KalDet
- bool extensionDataValid = true;
- try{
- m_vxdData = vxd.extension<dd4hep::rec::ZPlanarData>();
- }
- catch(std::runtime_error& e){
- extensionDataValid = false;
- info() << e.what() << " " << m_vxdData << endmsg;
- }
-
- std::vector<helpLayer> helpSensitives;
- std::vector<helpLayer> helpLadders;
- std::vector<int> helpNumberLadders;
- std::vector<double> helpPhi0;
- int helpCount=0;
- int type=0;
- double shellInnerRadius, shellOuterRadius, shellHalfLength, gap, shellRadLength;
- int nLadders=0;
- double phi0=0;
- helpLayer thisLadder;
- double beryllium_ladder_block_length=0,end_electronics_half_z=0,side_band_electronics_width=0;
- double rAlu, drAlu, rSty, drSty, dzSty, rInner, aluEndcapZ, aluHalfZ, alu_RadLen, Cryostat_dEdx;
- double VXDSupportDensity, VXDSupportZeff, VXDSupportAeff, VXDSupportRadLen;
- dd4hep::Volume vxd_vol = vxd.volume();
- for(int i=0;i<vxd_vol->GetNdaughters();i++){
- TGeoNode* daughter = vxd_vol->GetNode(i);
- std::string nodeName = daughter->GetName();
- //info() << daughter->GetName() << endmsg;
- if(nodeName=="VXD_support_assembly_0"){
- TGeoNode* shell = FindNode(daughter, "SupportShell");
- if(shell){
- const TGeoShape* shape_shell = shell->GetVolume()->GetShape();
- //fucd: IsA() method does not always work for TGeoTube, sometimes, strange?
- //if(shape_shell->IsA()==TGeoTube::Class()){
- if(shape_shell->TestShapeBit(TGeoTube::kGeoTube)){
- const TGeoTube* tube = (const TGeoTube*) shape_shell;
- shellInnerRadius = tube->GetRmin()*CLHEP::cm;
- shellOuterRadius = tube->GetRmax()*CLHEP::cm;
- shellHalfLength = tube->GetDz()*CLHEP::cm;
- }
- else{
- error() << shell->GetName() << " is not a TGeoTube! Shape bits = " << shape_shell->TestShapeBits(0xFFFFFFFF) << endmsg;
- }
- TGeoMaterial* mat = shell->GetMedium()->GetMaterial();
- shellRadLength = mat->GetRadLen()*CLHEP::cm;
- }
- TGeoNode* block = FindNode(daughter, "BerylliumAnnulusBlock");
- if(block){
- const TGeoShape* shape_block = block->GetVolume()->GetShape();
- if(shape_block->IsA()==TGeoBBox::Class()){
- const TGeoBBox* box = (const TGeoBBox*) shape_block;
- beryllium_ladder_block_length = box->GetDY()*CLHEP::cm;
- }
- else{
- error() << block->GetName() << " is not a TGeoTube! Shape bits = " << shape_block->TestShapeBits(0xFFFFFFFF) << endmsg;
- }
- }
- TGeoNode* skin = FindNode(daughter, "CryostatAluSkinBarrel");
- if(skin){
- const TGeoShape* shape_skin = skin->GetVolume()->GetShape();
- if(shape_skin->TestShapeBit(TGeoTube::kGeoTube)){
- const TGeoTube* tube = (const TGeoTube*) shape_skin;
- rAlu = tube->GetRmin()*CLHEP::cm;
- drAlu = tube->GetRmax()*CLHEP::cm - rAlu;
- aluHalfZ = tube->GetDz()*CLHEP::cm;
- //info() << rmin << "," << rmax << "," << zhalf << endmsg;
- }
- else{
- error() << skin->GetName() << " is not a TGeoTube! Shape bits = " << shape_skin->TestShapeBits(0xFFFFFFFF) << endmsg;
- }
- }
- TGeoNode* foam = FindNode(daughter, "CryostatFoamBarrel");
- if(foam){
- const TGeoShape* shape_foam = foam->GetVolume()->GetShape();
- if(shape_foam->TestShapeBit(TGeoTube::kGeoTube)){
- const TGeoTube* tube = (const TGeoTube*) shape_foam;
- rSty = tube->GetRmin()*CLHEP::cm;
- drSty = tube->GetRmax()*CLHEP::cm - rSty;
- dzSty = tube->GetDz()*CLHEP::cm;
- //info() << rmin << "," << rmax << "," << zhalf << endmsg;
- }
- else{
- error() << foam->GetName() << " is not a TGeoTube! Shape bits = " << shape_foam->TestShapeBits(0xFFFFFFFF) << endmsg;
- }
- }
- TGeoNode* skinEnd = FindNode(daughter, "CryostatAluSkinEndPlateInner");
- if(skinEnd){
- const TGeoShape* shape_skinEnd = skinEnd->GetVolume()->GetShape();
- if(shape_skinEnd->TestShapeBit(TGeoTube::kGeoTube)){
- const TGeoTube* tube = (const TGeoTube*) shape_skinEnd;
- rInner = tube->GetRmin()*CLHEP::cm;
- double rmax = tube->GetRmax()*CLHEP::cm;
- drAlu = tube->GetDz()*CLHEP::cm*2;
- //info() << rmin << "," << rmax << "," << zhalf << endmsg;
- }
- else{
- error() << skinEnd->GetName() << " is not a TGeoTube! Shape bits = " << shape_skinEnd->TestShapeBits(0xFFFFFFFF) << endmsg;
- }
- }
- TGeoNode* shellEnd = FindNode(daughter, "EndPlateShell_outer");
- if(shellEnd){
- const TGeoShape* shape_shellEnd = shellEnd->GetVolume()->GetShape();
- if(shape_shellEnd->TestShapeBit(TGeoTube::kGeoTube)){
- const TGeoTube* tube = (const TGeoTube*) shape_shellEnd;
- double rmin = tube->GetRmin()*CLHEP::cm;
- double rmax = tube->GetRmax()*CLHEP::cm;
- double zhalf = tube->GetDz()*CLHEP::cm;
- //info() << rmin << "," << rmax << "," << zhalf << endmsg;
- }
- else{
- error() << shellEnd->GetName() << " is not a TGeoTube! Shape bits = " << shape_shellEnd->TestShapeBits(0xFFFFFFFF) << endmsg;
- }
- }
-
- }
- else if(nodeName=="layer_assembly_0_1"){
- if(TGeoNode* side_band = FindNode(daughter, "ElectronicsBand")){
- const TGeoShape* shape_band = side_band->GetVolume()->GetShape();
- if(shape_band->IsA()==TGeoBBox::Class()){
- const TGeoBBox* box = (const TGeoBBox*) shape_band;
- side_band_electronics_width = box->GetDX()*CLHEP::cm*2;
- //info() << "fucd: "<< box->GetDX() << " " << box->GetDY() << " " << box->GetDZ() <<endmsg;
- }
- else{
- error() << "ElectronicsBand is not a TGeoBBox!!!"<< endmsg;
- }
- }
- if(TGeoNode* end = FindNode(daughter, "ElectronicsEnd")){
- const TGeoShape* shape_end = end->GetVolume()->GetShape();
- if(shape_end->IsA()==TGeoBBox::Class()){
- const TGeoBBox* box = (const TGeoBBox*) shape_end;
- end_electronics_half_z = box->GetDY()*CLHEP::cm;
- //info() << "fucd: " << box->GetDX() << " " << box->GetDY() << " " << box->GetDZ() << endmsg;
- }
- else{
- error() << "ElectronicsEnd is not a TGeoBBox!!!"<< endmsg;
- }
- }
- }
-
- /*
- for(int j=0;j<daughter->GetNdaughters();j++){
- TGeoNode* next = daughter->GetDaughter(j);
- info() << "fucd: " << next->GetName() << endmsg;
- }
- */
- }
-
- const std::map<std::string, dd4hep::DetElement>& components = vxd.children();
- for(std::map<std::string, dd4hep::DetElement>::const_iterator it=components.begin();it!=components.end();it++){
- dd4hep::DetElement component = it->second;
- dd4hep::Volume vol = component.volume();
- dd4hep::PlacedVolume phys = component.placement();
- TGeoMaterial* mat = vol->GetMaterial();
- const TGeoShape* shape = vol->GetShape();
- const dd4hep::PlacedVolumeExtension::VolIDs& ids = phys.volIDs();
- //info() << " " << it->first << " " << vol->GetName() << " " << component.id() << " " << component.path() << " " << component.placementPath() << endmsg;
- //info() << " " << shape->GetName() << " " << vol.solid().name() << endmsg;
- //info() << " " << mat->GetName() << " " << mat->GetRadLen() << endmsg;
- //info() << " " << ids.str() << endmsg;
- //info() << " " << vol->GetNdaughters() << endmsg;
- //double radL = mat->GetRadLen();
- //dd4hep::Solid solid = vol.solid();
- //info() << " " << sh->TestShapeBit(TGeoShape::kGeoBox) << " " << sh->GetName() << " " << phys.material().radLength() << endmsg;
- if(vol.isSensitive()&&shape->IsA()==TGeoBBox::Class()){
- int iLayer = ids.find("layer")->second;
- int iModule = ids.find("module")->second;
- int iSide = ids.find("side")->second;
- //info() << " layer=" << iLayer << " module=" << iModule << mat->GetName() << endmsg;
- if(iModule==0&&iLayer==helpCount+1){
- helpCount++;
- helpSensitives.push_back(thisLadder);
- helpLadders.push_back(thisLadder);
- helpNumberLadders.push_back(nLadders);
- helpPhi0.push_back(phi0);
- nLadders = 0;
- thisLadder.length = 0;
- }
- if(iLayer == helpCount){
- if(iModule == 0){
- const TGeoBBox* box = (const TGeoBBox*) shape;
- double width = box->GetDX()*CLHEP::cm;
- double length = box->GetDY()*CLHEP::cm;
- double thickness = box->GetDZ()*CLHEP::cm;
- TGeoMatrix* matrix = phys->GetMatrix();
- const double* pos = matrix->GetTranslation();
- const double* rot_data = matrix->GetRotationMatrix();
- TGeoRotation rot;
- rot.SetMatrix(rot_data);
- double theta,phi,psi;
- rot.GetAngles(phi,theta,psi);
- phi *= TMath::DegToRad();
- theta *= TMath::DegToRad();
- psi *= TMath::DegToRad();
- phi0 = -TMath::PiOver2()+phi;
- double distance = fabs(cos(phi0)*sin(theta)*pos[0]+sin(phi0)*sin(theta)*pos[1]+cos(theta)*pos[2]);
- double offset = sqrt(pos[0]*pos[0]+pos[1]*pos[1]-distance*distance)*pos[0]/fabs(pos[0])*CLHEP::cm;
- distance *= CLHEP::cm;
- distance -= thickness;
- double radL = mat->GetRadLen()*CLHEP::cm;
- //info() << " -> " << helpCount << ": " << distance << " " << cos(atan2(pos[1],pos[0])-phi)*sqrt(pos[0]*pos[0]+pos[1]*pos[1]) << endmsg;
- thisLadder.distance = distance;
- thisLadder.offset = offset;
- thisLadder.thickness = thickness;
- thisLadder.length += length;
- thisLadder.width = width;
- thisLadder.radLength = radL;
- thisLadder.z = pos[2]*CLHEP::cm;
- }
- if(iModule==nLadders) nLadders++;
- }
- }
- //info() << " " << vol.solid().type() << " " << vol.solid().name() << " " << vol->GetNdaughters() << endmsg;
- else if(it->first=="VXD_support"){
- helpCount++;
- helpSensitives.push_back(thisLadder);
- helpLadders.push_back(thisLadder);
- helpNumberLadders.push_back(nLadders);
- helpPhi0.push_back(phi0);
- nLadders = 0;
- if(vol->GetNdaughters()==0) error() << "!!!!!!!!!" << endmsg;
-
- int nFlexCable = 0, nFoamSpacer=0, nMetalTraces=0;
- int currentLayer = -1;
- double tFlexCable, tFoamSpacer, tMetalTraces;
- double radLFlexCable, radLFoamSpacer, radLMetalTraces;
- double intLFlexCable, intLFoamSpacer, intLMetalTraces;
- double dFlexCable, dFoamSpacer, dMetalTraces;
- double metalZeff, metalZAeff, foamZeff, foamZAeff, flexZeff, flexZAeff;
- for(int i=0;i<vol->GetNdaughters();i++){
- TGeoNode* daughter = vol->GetNode(i);
- TGeoMaterial* matDaughter = daughter->GetMedium()->GetMaterial();
- const TGeoShape* shape_sup = daughter->GetVolume()->GetShape();
- TGeoMatrix* matrix = daughter->GetMatrix();
- const double* pos = matrix->GetTranslation();
- const double* rot_data = matrix->GetRotationMatrix();
- TGeoRotation rot;
- rot.SetMatrix(rot_data);
- double theta,phi,psi;
- rot.GetAngles(phi,theta,psi);
- phi *= TMath::DegToRad();
- theta *= TMath::DegToRad();
- psi *= TMath::DegToRad();
- phi0 = -TMath::PiOver2()+phi;
- std::string phy_name = daughter->GetName();
- if(phy_name.find("FoamSpacer")==-1&&phy_name.find("FlexCable")==-1&&phy_name.find("MetalTraces")==-1){
- //info() << phy_name <<endmsg;
- continue;
- }
- int iLayer = atoi(phy_name.substr(phy_name.find("_")+1,2).c_str());
- if(iLayer!=currentLayer){
- //info() << tFoamSpacer << "," << tFlexCable << "," << tMetalTraces << endmsg;
- helpLadders[currentLayer].thickness = tFoamSpacer+tFlexCable+tMetalTraces;
- helpLadders[currentLayer].radLength = helpLadders[currentLayer].thickness / (tFoamSpacer/radLFoamSpacer+tFlexCable/radLFlexCable+tMetalTraces/radLMetalTraces);
- nFlexCable = 0;
- nFoamSpacer=0;
- nMetalTraces=0;
- currentLayer=iLayer;
- }
- //info() << "ss pos=" << pos[0] << "," << pos[1] << "," << pos[2] << " distance=" << sqrt(pos[0]*pos[0]+pos[1]*pos[1]) << endmsg;
- //info() << "ss rot=" << phi << "," << theta << "," << psi << endmsg;
- //info() << "ss " << daughter->GetName() << " " << daughter->GetVolume()->GetName() << " " << endmsg;
- if(shape_sup->IsA()==TGeoBBox::Class()&&(nFoamSpacer==0||nFlexCable==0||nMetalTraces==0)){
- const TGeoBBox* box = (const TGeoBBox*) shape_sup;
- //info() << phy_name.substr(phy_name.find("_")+1,2) << " " << iLayer << " " << box->GetDX() << "," << box->GetDY() << "," << box->GetDZ() << endmsg;
- //info() << "fucd: pos " << pos[0] << " " << pos[1] << " " << pos[2] << endmsg;
- double distance = fabs(cos(phi0)*sin(theta)*pos[0]+sin(phi0)*sin(theta)*pos[1]+cos(theta)*pos[2]);
- double offset = sqrt(pos[0]*pos[0]+pos[1]*pos[1]-distance*distance)*pos[0]/fabs(pos[0])*CLHEP::cm;
- distance -= box->GetDZ();
- distance *= CLHEP::cm;
- if(helpLadders[iLayer].distance == helpSensitives[iLayer].distance) helpLadders[iLayer].distance = distance;
- else helpLadders[iLayer].distance = TMath::Min(helpLadders[iLayer].distance, distance);
- //info() << phy_name << " " << distance << " " << offset << endmsg;
- if(phy_name.find("FoamSpacer")!=-1&&nFoamSpacer==0){
- helpLadders[iLayer].offset = offset;
- tFoamSpacer = box->GetDZ()*CLHEP::cm;
- radLFoamSpacer = matDaughter->GetRadLen()*CLHEP::cm;
- intLFoamSpacer = matDaughter->GetIntLen()*CLHEP::cm;
- dFoamSpacer = matDaughter->GetDensity()*CLHEP::g/CLHEP::cm3;
- //fucd: A calculated by TGeoMaterial class is not equal to Zeff/ZAeff, Zeff = sum(Zi*Ai/sumA), ZAeff = sum(Zi/Ai*A/totalA)
- // use Zeff and ZAeff to keep same with Mokka case
- //foamZ = matDaughter->GetZ();
- //foamA = matDaughter->GetA();
- double totalA = 0, Zeff = 0, ZAeff = 0;
- for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
- totalA += matDaughter->GetElement(iEle)->A();
- }
- for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
- double A, Z, w;
- // by fucd: w!=A/totalA, strange! to fix
- matDaughter->GetElementProp(A,Z,w,iEle);
- Zeff += Z*w;
- ZAeff += Z/A*w;
- //info() << std::setprecision(16) << Z << " " << A << " " << A/totalA << " " << w << endmsg;
- }
- foamZeff = Zeff;
- foamZAeff = ZAeff;
- nFoamSpacer++;
- }
- if(phy_name.find("FlexCable")!=-1&&nFlexCable==0){
- helpLadders[iLayer].width = box->GetDX()*CLHEP::cm;
- helpLadders[iLayer].length = box->GetDY()*CLHEP::cm-beryllium_ladder_block_length*2-end_electronics_half_z*2;
- tFlexCable = box->GetDZ()*CLHEP::cm;
- radLFlexCable = matDaughter->GetRadLen()*CLHEP::cm;
- intLFlexCable = matDaughter->GetIntLen()*CLHEP::cm;
- dFlexCable = matDaughter->GetDensity()*CLHEP::g/CLHEP::cm3;
- double Zeff = 0, ZAeff = 0;
- for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
- double A, Z, w;
- matDaughter->GetElementProp(A,Z,w,iEle);
- Zeff += Z*w;
- ZAeff += Z/A*w;
- //std::cout << std::setprecision(16) << Z << " " << A << " " << w << std::endl;
- }
- flexZeff = Zeff;
- flexZAeff = ZAeff;
- nFlexCable++;
- }
- if(phy_name.find("MetalTraces")!=-1&&nMetalTraces==0){
- tMetalTraces = box->GetDZ()*CLHEP::cm;
- radLMetalTraces = matDaughter->GetRadLen()*CLHEP::cm;
- intLMetalTraces = matDaughter->GetIntLen()*CLHEP::cm;
- dMetalTraces = matDaughter->GetDensity()*CLHEP::g/CLHEP::cm3;
- double totalA = 0, Zeff = 0, ZAeff = 0;
- for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
- totalA += matDaughter->GetElement(iEle)->A();
- }
- for(int iEle = 0; iEle<matDaughter->GetNelements(); iEle++){
- double A, Z, w;
- matDaughter->GetElementProp(A,Z,w,iEle);
- Zeff += Z*w;
- ZAeff += Z/A*w;
- //info() << Z << " " << A << " " << w << endmsg;
- }
- metalZeff = Zeff;
- metalZAeff = ZAeff;
- nMetalTraces++;
- }
- }
- }
- {
- //info() << tFoamSpacer << "," << tFlexCable << "," << tMetalTraces << endmsg;
- double tSupport = tMetalTraces + tFoamSpacer + tFlexCable;
- helpLadders[currentLayer].thickness = tSupport;
- helpLadders[currentLayer].radLength = helpLadders[currentLayer].thickness / (tFoamSpacer/radLFoamSpacer+tFlexCable/radLFlexCable+tMetalTraces/radLMetalTraces);
- nFlexCable = 0;
- nFoamSpacer=0;
- nMetalTraces=0;
-
- //calculations of thickness fractions of each layer of the support
- double metalTF = tMetalTraces / tSupport;
- double foamTF = tFoamSpacer / tSupport;
- double flexTF = tFlexCable / tSupport;
- //info() << foamTF << "," << flexTF << "," << metalTF << endmsg;
- //info() << dFoamSpacer/(CLHEP::kg/CLHEP::cm3) << "," << dFlexCable/(CLHEP::kg/CLHEP::cm3) << "," << dMetalTraces/(CLHEP::kg/CLHEP::cm3) << endmsg;
- //info() << foamZeff << " " << flexZeff << " " << metalZeff << endmsg;
- //info() << foamZAeff << " " << flexZAeff << " " << metalZAeff << endmsg;
- G4double elemVol = 1*CLHEP::cm3;
- G4double VXDSupportMass = (foamTF*dFoamSpacer + flexTF*dFlexCable + metalTF*dMetalTraces)*elemVol;
- VXDSupportDensity = VXDSupportMass/elemVol;
-
- G4double foamFM = 100. * ((foamTF*(elemVol)*dFoamSpacer) / VXDSupportMass) ;
- G4double kaptonFM = 100. * ((flexTF*(elemVol)*dFlexCable) / VXDSupportMass) ;
- G4double metalFM = 100. * ((metalTF*(elemVol)*dMetalTraces) / VXDSupportMass) ;
-
- VXDSupportRadLen = helpLadders[currentLayer].radLength;
-
- VXDSupportZeff = (metalFM/100.)*metalZeff + (kaptonFM/100.)*flexZeff + (foamFM/100.)*foamZeff;
- G4double VXDSupportZAeff = (metalFM/100.)*metalZAeff + (kaptonFM/100.)*flexZAeff + (foamFM/100.)*foamZAeff;
- VXDSupportAeff = VXDSupportZeff / VXDSupportZAeff;
- G4double VXDSupportIntLength = 1. / ((metalTF/intLMetalTraces) + (flexTF/intLFlexCable) + (foamTF/intLFoamSpacer));
- VXDSupportDensity = 1000*VXDSupportDensity/(CLHEP::g/CLHEP::cm3);
- //info() << "fucd: " << VXDSupportZeff << " " << VXDSupportAeff << " " << VXDSupportRadLen << " " << VXDSupportIntLength << " " << VXDSupportDensity << endmsg;
- //info() << intLMetalTraces << " " << intLFlexCable << " " << intLFoamSpacer <<endmsg;
- }
- }
- //info() << it->first << endmsg;
- }
- if(end_electronics_half_z>0 && side_band_electronics_width==0) type = gear::ZPlanarParametersImpl::CCD ;
- if(side_band_electronics_width>0 && end_electronics_half_z==0 ) type = gear::ZPlanarParametersImpl::CMOS ;
- if(side_band_electronics_width>0 && end_electronics_half_z>0) type = gear::ZPlanarParametersImpl::HYBRID ;
-
- m_vxdParameters = new gear::ZPlanarParametersImpl(type, shellInnerRadius, shellOuterRadius, shellHalfLength, gap, shellRadLength );
- // by fucd: debug info, if validated enough, merge them in future
- info() << "=====================from convertor==============================" << endmsg;
- info() << type << " " << shellInnerRadius << " " << shellOuterRadius << " " << shellHalfLength << " " << gap << " " << shellRadLength << endmsg;
- for(int i=0;i<helpCount;i++){
- m_vxdParameters->addLayer(helpNumberLadders[i] , helpPhi0[i] ,
- helpLadders[i].distance , helpLadders[i].offset, helpLadders[i].thickness*2 ,
- helpLadders[i].length , helpLadders[i].width*2 , helpLadders[i].radLength ,
- helpSensitives[i].distance, helpSensitives[i].offset , helpSensitives[i].thickness*2 ,
- helpSensitives[i].length , helpSensitives[i].width*2 , helpSensitives[i].radLength ) ;
- info() << "fucd " << i << ": " << helpNumberLadders[i] << ", " << helpPhi0[i] << ", "
- << helpLadders[i].distance << ", " << helpLadders[i].offset << ", " << helpLadders[i].thickness*2 << ", " << helpLadders[i].length << ", "
- << helpLadders[i].width*2 << ", " << helpLadders[i].radLength << ", " << helpSensitives[i].distance << ", " << helpSensitives[i].offset << ", "
- << helpSensitives[i].thickness*2 << ", " << helpSensitives[i].length << ", " << helpSensitives[i].width*2 << ", " << helpSensitives[i].radLength << endmsg;
- }
- //m_vxdInfra = new gear::GearParametersImpl;
- //CryostatAlRadius, CryostatAlThickness, CryostatAlInnerR, CryostatAlZEndCap, CryostatAlHalfZ
- //m_vxdInfra->setDoubleVal("CryostatAlRadius",rAlu);
- //m_vxdInfra->setDoubleVal("CryostatAlThickness",drAlu);
- //m_vxdInfra->setDoubleVal("CryostatAlInnerR",rInner);
- //m_vxdInfra->setDoubleVal("CryostatAlZEndCap",aluEndcapZ=dzSty + drSty + drAlu / 2);
- //m_vxdInfra->setDoubleVal("CryostatAlHalfZ",aluHalfZ= dzSty + drSty);
- // change GearParametersImpl to map
- std::map<std::string,double> vxdInfra;
- vxdInfra["CryostatAlRadius"] = rAlu;
- vxdInfra["CryostatAlThickness"] = drAlu;
- vxdInfra["CryostatAlInnerR"] = rInner;
- vxdInfra["CryostatAlZEndCap"] = dzSty+drSty+drAlu/2;
- vxdInfra["CryostatAlHalfZ"] = dzSty+drSty;
- m_detParameters["VXDInfra"] = vxdInfra;
- //effective A different with what in Mokka, fix them as Mokka's
- //m_materials["VXDSupportMaterial"] = new TMaterial("VXDSupportMaterial", "", VXDSupportAeff, VXDSupportZeff, VXDSupportDensity, VXDSupportRadLen, 0.);
- m_materials["VXDSupportMaterial"] = new TMaterial("VXDSupportMaterial", "", 2.075865162e+01, 1.039383117e+01, 2.765900000e+02/1000, 1.014262421e+02, 0.);
-
- info() << "=====================from ZPlanarData==============================" << endmsg;
- if(m_vxdData){
- info() << m_vxdData->rInnerShell << " " << m_vxdData->rOuterShell << " " << m_vxdData->zHalfShell << " " << m_vxdData->gapShell << endmsg;
- const std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& layers = m_vxdData->layers;
- for(int i=0;i<layers.size();i++){
- const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = layers[i];
- info() << i << ": " << thisLayer.ladderNumber << "," << thisLayer.phi0 << "," << thisLayer.distanceSupport << "," << thisLayer.offsetSupport << ","
- << thisLayer.thicknessSupport << "," << thisLayer.zHalfSupport << "," << thisLayer.widthSupport << "," << "NULL,"
- << thisLayer.distanceSensitive << "," << thisLayer.offsetSensitive << "," << thisLayer.thicknessSensitive << "," << thisLayer.zHalfSensitive << ","
- << thisLayer.widthSensitive << ",NULL" << endmsg;
- }
- }
- info() << rAlu << " " << drAlu << " " << rInner << " " << aluEndcapZ << " " << aluHalfZ << endmsg;
- //info() << m_materials["VXDSupportMaterial"] << endmsg;
- return sc;
-}
-
-TGeoNode* GeomSvc::FindNode(TGeoNode* mother, char* name){
- TGeoNode* next = 0;
- if(mother->GetNdaughters()!=0){
- for(int i=0;i<mother->GetNdaughters();i++){
- TGeoNode* daughter = mother->GetDaughter(i);
- std::string s = daughter->GetName();
- //info() << "current: " << s << " search for" << name << endmsg;
- if(s.find(name)!=-1){
- next = daughter;
- break;
- }
- else{
- next = FindNode(daughter, name);
- }
- }
- }
- return next;
-}
-
-TMaterial* GeomSvc::getMaterial(std::string name){
- std::map<std::string, TMaterial*>::const_iterator it = m_materials.find(name);
- if(it!=m_materials.end()) return it->second;
- else return 0;
-
-}
diff --git a/Detector/GeomSvc/src/GeomSvc.h b/Detector/GeomSvc/src/GeomSvc.h
index 4ef39d2f0644c25b666145d4ee4a5caffd8420e1..119e2ea5870db041ad48bf4bf0175af290af51cd 100644
--- a/Detector/GeomSvc/src/GeomSvc.h
+++ b/Detector/GeomSvc/src/GeomSvc.h
@@ -34,45 +34,15 @@ class GeomSvc: public extends<Service, IGeomSvc> {
dd4hep::DetElement getDD4HepGeo() override;
dd4hep::Detector* lcdd() override;
- const gear::ZPlanarParametersImpl* getVXDParameters() override {return m_vxdParameters;};
- const dd4hep::rec::ZPlanarData* getVXDData() override {return m_vxdData;};
- const dd4hep::rec::ConicalSupportData* getBeamPipeData() override {return m_beamPipeData;};
-
- const std::map<std::string,double>& getDetParameters(std::string name) override;
- double getDetParameter(std::string set_name, std::string par_name) override;
- TMaterial* getMaterial(std::string name);
-
private:
- StatusCode convertVXD(dd4hep::DetElement& sub);
-
- Decoder* getDecoder(const std::string& readout_name) override;
-
+ Decoder* getDecoder(const std::string& readout_name) override;
+
private:
- TGeoNode* FindNode(TGeoNode* mother, char* name);
// DD4hep XML compact file path
Gaudi::Property<std::string> m_dd4hep_xmls{this, "compact"};
//
dd4hep::Detector* m_dd4hep_geo;
-
-
- gear::ZPlanarParametersImpl* m_vxdParameters{nullptr};
- dd4hep::rec::ZPlanarData* m_vxdData{nullptr};
- dd4hep::rec::ConicalSupportData* m_beamPipeData{nullptr};
-
- //gear::GearParametersImpl* m_vxdInfra;
- std::map<std::string, std::map<std::string,double> > m_detParameters;
- std::map<std::string, TMaterial*> m_materials;
- struct helpLayer {
- double distance =0;
- double offset =0;
- double thickness =0;
- double length =0;
- double width =0;
- double radLength =0;
- double z =0;
- double foam_spacer_radLength =0;
- };
};
#endif // GeomSvc_h
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h
index 87debdd87da87e2e46ae6305103e421e00ac03f2..d645f2db44812b4cbdd2281f66b7ce7657be6824 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h
@@ -25,7 +25,8 @@
#include <string>
-typedef std::vector<edm4hep::CalorimeterHit *> CalorimeterHitVector;
+//typedef std::vector<edm4hep::CalorimeterHit *> CalorimeterHitVector;
+typedef std::vector<const edm4hep::CalorimeterHit *> CalorimeterHitVector;
namespace gear { class GearMgr; }
@@ -174,25 +175,25 @@ private:
* @brief Get common calo hit properties: position, parent address, input energy and time
*
*/
- void GetCommonCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, PandoraApi::CaloHit::Parameters &caloHitParameters) const;
+ void GetCommonCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, PandoraApi::CaloHit::Parameters &caloHitParameters) const;
/**
* @brief Get end cap specific calo hit properties: cell size, absorber radiation and interaction lengths, normal vector
*
*/
- void GetEndCapCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
+ void GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const;
- void GetEndCapCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+ void GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const;
/**
* @brief Get barrel specific calo hit properties: cell size, absorber radiation and interaction lengths, normal vector
*
*/
- void GetBarrelCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
+ void GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
unsigned int barrelSymmetryOrder, float barrelPhi0, unsigned int staveNumber, PandoraApi::CaloHit::Parameters &caloHitParameters,
float &absorberCorrection) const;
- void GetBarrelCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+ void GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
unsigned int barrelSymmetryOrder, float barrelPhi0, unsigned int staveNumber, PandoraApi::CaloHit::Parameters &caloHitParameters,
float &absorberCorrection) const;
@@ -200,13 +201,13 @@ private:
* @brief Get number of active layers from position of a calo hit to the edge of the detector
*
*/
- int GetNLayersFromEdge(edm4hep::CalorimeterHit *const pCaloHit) const;
+ int GetNLayersFromEdge(const edm4hep::CalorimeterHit *const pCaloHit) const;
/**
* @brief Get the maximum radius of a calo hit in a polygonal detector structure
*
*/
- float GetMaximumRadius(edm4hep::CalorimeterHit *const pCaloHit, const unsigned int symmetryOrder, const float phi0) const;
+ float GetMaximumRadius(const edm4hep::CalorimeterHit *const pCaloHit, const unsigned int symmetryOrder, const float phi0) const;
/**
* @brief Get the layer coding string from the provided cell id encoding string
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/MCParticleCreator.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/MCParticleCreator.h
index ba37ad54820b5970d0da2ee150d42b57d2987ea4..4e446cf9c59ed15ea3f28034f99adedd8f407917 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/MCParticleCreator.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/MCParticleCreator.h
@@ -78,7 +78,7 @@ private:
const Settings m_settings; ///< The mc particle creator settings
const pandora::Pandora *m_pPandora; ///< Address of the pandora object to create the mc particles
const float m_bField; ///< The bfield
- std::map<unsigned int, edm4hep::MCParticle*>* m_id_pMC_map;
+ std::map<unsigned int, const edm4hep::MCParticle*>* m_id_pMC_map;
};
inline void MCParticleCreator::Reset()
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h
index b8b8ec1844545b8ea0456a5f95b284e58e58eaf3..e6059bc7a36bd3d569ef6711d3d9333dc31be42b 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h
@@ -22,7 +22,7 @@ namespace gear { class GearMgr; }
class CollectionMaps;
-typedef std::vector<edm4hep::Track *> TrackVector;
+typedef std::vector<const edm4hep::Track *> TrackVector;
typedef std::set<unsigned int> TrackList;
typedef std::map<edm4hep::Track, int> TrackToPidMap;
/*
@@ -123,7 +123,7 @@ public:
*/
pandora::StatusCode CreateTrackAssociations(const CollectionMaps& collectionMaps);
- edm4hep::Track* GetTrackAddress(const CollectionMaps& collectionMaps, const edm4hep::Track& pTrack );
+ const edm4hep::Track* GetTrackAddress(const CollectionMaps& collectionMaps, const edm4hep::Track& pTrack );
/**
* @brief Create tracks, insert user code here
*
@@ -196,7 +196,7 @@ private:
* @brief Copy track states stored in tracks to pandora track parameters
*
*/
- void GetTrackStates(edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const;
+ void GetTrackStates(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const;
/**
* @brief Copy track state from track state instance to pandora input track state
@@ -208,13 +208,13 @@ private:
* @brief Obtain track time when it reaches ECAL
*
*/
- float CalculateTrackTimeAtCalorimeter(edm4hep::Track *const pTrack) const;
+ float CalculateTrackTimeAtCalorimeter(const edm4hep::Track *const pTrack) const;
/**
* @brief Decide whether track reaches the ecal surface
*
*/
- void TrackReachesECAL(edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const;
+ void TrackReachesECAL(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const;
/**
* @brief Determine whether a track can be used to form a pfo under the following conditions:
@@ -222,7 +222,7 @@ private:
* 2) if the track proves to have no cluster associations
*
*/
- void DefineTrackPfoUsage(edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const;
+ void DefineTrackPfoUsage(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const;
/**
* @brief Whether track passes the quality cuts required in order to be used to form a pfo
@@ -230,19 +230,19 @@ private:
*
* @return boolean
*/
- bool PassesQualityCuts(edm4hep::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const;
+ bool PassesQualityCuts(const edm4hep::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const;
/**
* @brief Get number of hits in TPC of a track
*
*/
- int GetNTpcHits(edm4hep::Track *const pTrack) const;
+ int GetNTpcHits(const edm4hep::Track *const pTrack) const;
/**
* @brief Get number of hits in FTD of a track
*
*/
- int GetNFtdHits(edm4hep::Track *const pTrack) const;
+ int GetNFtdHits(const edm4hep::Track *const pTrack) const;
const Settings m_settings; ///< The track creator settings
const pandora::Pandora *m_pPandora; ///< Address of the pandora object to create tracks and track relationships
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
index d13b5e5ad5c709e0daf3a3eec0fc6bd819d65ff4..29b9adcbfa023d2abbb5280e7a028e4aff224439 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
@@ -240,7 +240,7 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col
try
{
if(m_settings.m_debug) std::cout<<"CaloHitCreator for "<<tmp_col_name<<std::endl;
- auto pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(tmp_col_name))->second;
+ auto & pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(tmp_col_name))->second;
const int nElements(pCaloHitCollection.size());
if (0 == nElements)
@@ -276,8 +276,8 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col
{
try
{
- auto pCaloHit0 = pCaloHitCollection.at(i);
- auto pCaloHit = &(pCaloHit0);
+ //auto pCaloHit = const_cast<edm4hep::CalorimeterHit*>(&(pCaloHitCollection.at(i)));
+ auto pCaloHit = &(pCaloHitCollection.at(i));
if (NULL == pCaloHit)
throw ("CreateECalCaloHits pCaloHit Collection type mismatch");
@@ -395,7 +395,7 @@ pandora::StatusCode CaloHitCreator::CreateHCalCaloHits(const CollectionMaps& col
if(collectionMaps.collectionMap_CaloHit.find(tmp_col_name) == collectionMaps.collectionMap_CaloHit.end()) { if(m_settings.m_debug) std::cout<<"not find "<<tmp_col_name<<std::endl; continue;}
try
{
- auto pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(tmp_col_name))->second;
+ auto & pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(tmp_col_name))->second;
const int nElements(pCaloHitCollection.size());
if (0 == nElements)
@@ -431,8 +431,7 @@ pandora::StatusCode CaloHitCreator::CreateHCalCaloHits(const CollectionMaps& col
{
try
{
- auto pCaloHit0 = pCaloHitCollection.at(i);
- auto pCaloHit = &(pCaloHit0);
+ auto pCaloHit = &(pCaloHitCollection.at(i));
if (NULL == pCaloHit)
throw ("CreateHCalCaloHits Collection type mismatch");
@@ -520,7 +519,7 @@ pandora::StatusCode CaloHitCreator::CreateMuonCaloHits(const CollectionMaps& col
if(collectionMaps.collectionMap_CaloHit.find(tmp_col_name) == collectionMaps.collectionMap_CaloHit.end()) {if(m_settings.m_debug) std::cout<<"not find "<<tmp_col_name<<std::endl; continue;}
try
{
- auto pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(tmp_col_name))->second;
+ auto & pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(tmp_col_name))->second;
const int nElements(pCaloHitCollection.size());
if (0 == nElements)
@@ -558,8 +557,7 @@ pandora::StatusCode CaloHitCreator::CreateMuonCaloHits(const CollectionMaps& col
{
try
{
- auto pCaloHit0 = pCaloHitCollection.at(i);
- auto pCaloHit = &(pCaloHit0);
+ auto pCaloHit = &(pCaloHitCollection.at(i));
if (NULL == pCaloHit)
throw ("Muon Collection type mismatch");
@@ -661,7 +659,7 @@ pandora::StatusCode CaloHitCreator::CreateLCalCaloHits(const CollectionMaps& col
if(collectionMaps.collectionMap_CaloHit.find(*iter) == collectionMaps.collectionMap_CaloHit.end()) {if(m_settings.m_debug)std::cout<<"not find "<<(*iter)<<std::endl; continue;}
try
{
- auto pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(*iter))->second;
+ auto & pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(*iter))->second;
const int nElements(pCaloHitCollection.size());
if (0 == nElements)
@@ -679,8 +677,7 @@ pandora::StatusCode CaloHitCreator::CreateLCalCaloHits(const CollectionMaps& col
{
try
{
- auto pCaloHit0 = pCaloHitCollection.at(i);
- auto pCaloHit = &(pCaloHit0);
+ auto pCaloHit = &(pCaloHitCollection.at(i));
if (NULL == pCaloHit)
throw ("LCal Collection type mismatch");
@@ -736,7 +733,7 @@ pandora::StatusCode CaloHitCreator::CreateLHCalCaloHits(const CollectionMaps& co
if(collectionMaps.collectionMap_CaloHit.find(*iter) == collectionMaps.collectionMap_CaloHit.end()) {if(m_settings.m_debug) std::cout<<"not find "<<(*iter)<<std::endl; continue;}
try
{
- auto pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(*iter))->second;
+ auto & pCaloHitCollection = (collectionMaps.collectionMap_CaloHit.find(*iter))->second;
const int nElements(pCaloHitCollection.size());
if (0 == nElements)
@@ -753,8 +750,7 @@ pandora::StatusCode CaloHitCreator::CreateLHCalCaloHits(const CollectionMaps& co
{
try
{
- auto pCaloHit0 = pCaloHitCollection.at(i);
- auto pCaloHit = &(pCaloHit0);
+ auto pCaloHit = &(pCaloHitCollection.at(i));
if (NULL == pCaloHit)
throw ("LHCal Collection type mismatch");
@@ -802,7 +798,7 @@ pandora::StatusCode CaloHitCreator::CreateLHCalCaloHits(const CollectionMaps& co
//------------------------------------------------------------------------------------------------------------------------------------------
-void CaloHitCreator::GetCommonCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, PandoraApi::CaloHit::Parameters &caloHitParameters) const
+void CaloHitCreator::GetCommonCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, PandoraApi::CaloHit::Parameters &caloHitParameters) const
{
const float pCaloHitPosition[3]={pCaloHit->getPosition()[0], pCaloHit->getPosition()[1], pCaloHit->getPosition()[2]};
const pandora::CartesianVector positionVector(pCaloHitPosition[0], pCaloHitPosition[1], pCaloHitPosition[2]);
@@ -817,7 +813,7 @@ void CaloHitCreator::GetCommonCaloHitProperties(edm4hep::CalorimeterHit *const p
//------------------------------------------------------------------------------------------------------------------------------------------
-void CaloHitCreator::GetEndCapCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
+void CaloHitCreator::GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const
{
caloHitParameters.m_hitRegion = pandora::ENDCAP;
@@ -862,7 +858,7 @@ void CaloHitCreator::GetEndCapCaloHitProperties(edm4hep::CalorimeterHit *const p
}
//------------------------------------------------------------------------------------------------------------------------------------------
-void CaloHitCreator::GetEndCapCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+void CaloHitCreator::GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const
{
caloHitParameters.m_hitRegion = pandora::ENDCAP;
@@ -914,7 +910,7 @@ void CaloHitCreator::GetEndCapCaloHitProperties(edm4hep::CalorimeterHit *const p
}
//------------------------------------------------------------------------------------------------------------------------------------------
-void CaloHitCreator::GetBarrelCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
+void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
unsigned int barrelSymmetryOrder, float barrelPhi0, unsigned int staveNumber, PandoraApi::CaloHit::Parameters &caloHitParameters,
float &absorberCorrection) const
{
@@ -977,7 +973,7 @@ void CaloHitCreator::GetBarrelCaloHitProperties(edm4hep::CalorimeterHit *const p
}
}
-void CaloHitCreator::GetBarrelCaloHitProperties(edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
unsigned int barrelSymmetryOrder, float barrelPhi0, unsigned int staveNumber, PandoraApi::CaloHit::Parameters &caloHitParameters,
float &absorberCorrection) const
{
@@ -1054,7 +1050,7 @@ void CaloHitCreator::GetBarrelCaloHitProperties(edm4hep::CalorimeterHit *const p
//------------------------------------------------------------------------------------------------------------------------------------------
-int CaloHitCreator::GetNLayersFromEdge(edm4hep::CalorimeterHit *const pCaloHit) const
+int CaloHitCreator::GetNLayersFromEdge(const edm4hep::CalorimeterHit *const pCaloHit) const
{
// Calo hit coordinate calculations
const float barrelMaximumRadius(this->GetMaximumRadius(pCaloHit, m_hCalBarrelOuterSymmetry, m_hCalBarrelOuterPhi0));
@@ -1096,7 +1092,7 @@ int CaloHitCreator::GetNLayersFromEdge(edm4hep::CalorimeterHit *const pCaloHit)
//------------------------------------------------------------------------------------------------------------------------------------------
-float CaloHitCreator::GetMaximumRadius(edm4hep::CalorimeterHit *const pCaloHit, const unsigned int symmetryOrder, const float phi0) const
+float CaloHitCreator::GetMaximumRadius(const edm4hep::CalorimeterHit *const pCaloHit, const unsigned int symmetryOrder, const float phi0) const
{
const float pCaloHitPosition[3]={pCaloHit->getPosition()[0], pCaloHit->getPosition()[1], pCaloHit->getPosition()[2]};
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp
index ef457a7b4a0691bd0df5666d59709efa87272844..1912eee92592f0b8071fcabed38259a61c049f87 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/MCParticleCreator.cpp
@@ -22,7 +22,7 @@ MCParticleCreator::MCParticleCreator(const Settings &settings, const pandora::Pa
m_pPandora(pPandora),
m_bField(settings.m_bField)
{
-m_id_pMC_map = new std::map<unsigned int, edm4hep::MCParticle*>;
+m_id_pMC_map = new std::map<unsigned int, const edm4hep::MCParticle*>;
}
//------------------------------------------------------------------------------------------------------------------------------------------
@@ -41,19 +41,20 @@ pandora::StatusCode MCParticleCreator::CreateMCParticles(const CollectionMaps& c
if(collectionMaps.collectionMap_MC.find(*iter) == collectionMaps.collectionMap_MC.end()) continue;
try
{
- auto pMCParticleCollection = (collectionMaps.collectionMap_MC.find(*iter))->second;
+ auto & pMCParticleCollection = (collectionMaps.collectionMap_MC.find(*iter))->second;
std::cout<<"Do CreateMCParticles, collection:"<<(*iter)<<", size="<<pMCParticleCollection.size()<<std::endl;
for (int im = 0; im < pMCParticleCollection.size(); im++)
{
try
{
- auto pMcParticle = pMCParticleCollection.at(im);
+ auto & pMcParticle = pMCParticleCollection.at(im);
PandoraApi::MCParticle::Parameters mcParticleParameters;
mcParticleParameters.m_energy = sqrt(pMcParticle.getMomentum()[0] * pMcParticle.getMomentum()[0] + pMcParticle.getMomentum()[1] * pMcParticle.getMomentum()[1] + pMcParticle.getMomentum()[2] * pMcParticle.getMomentum()[2] + pMcParticle.getMass() * pMcParticle.getMass());
mcParticleParameters.m_particleId = pMcParticle.getPDG();
mcParticleParameters.m_mcParticleType = pandora::MC_3D;
mcParticleParameters.m_pParentAddress = &pMcParticle;
unsigned int p_id = pMcParticle.id();
+ //auto p_mc = const_cast<edm4hep::MCParticle*>(&pMcParticle);
auto p_mc = &pMcParticle;
(*m_id_pMC_map) [p_id] = p_mc;
mcParticleParameters.m_momentum = pandora::CartesianVector(pMcParticle.getMomentum()[0], pMcParticle.getMomentum()[1],
@@ -104,7 +105,7 @@ pandora::StatusCode MCParticleCreator::CreateMCParticles(const CollectionMaps& c
pandora::StatusCode MCParticleCreator::CreateCaloHitToMCParticleRelationships(const CollectionMaps& collectionMaps, const CalorimeterHitVector &calorimeterHitVector) const
{
- typedef std::map<edm4hep::MCParticle *, float> MCParticleToEnergyWeightMap;
+ typedef std::map<const edm4hep::MCParticle *, float> MCParticleToEnergyWeightMap;
MCParticleToEnergyWeightMap mcParticleToEnergyWeightMap;
for (StringVector::const_iterator iter = m_settings.m_CaloHitRelationCollections.begin(), iterEnd = m_settings.m_CaloHitRelationCollections.end();
@@ -173,7 +174,8 @@ pandora::StatusCode MCParticleCreator::CreateTrackToMCParticleRelationships(cons
const pandora::Helix helixFit(pTrack->getTrackStates(0).phi, pTrack->getTrackStates(0).D0, pTrack->getTrackStates(0).Z0, pTrack->getTrackStates(0).omega, pTrack->getTrackStates(0).tanLambda, m_bField);
const float recoMomentum(helixFit.GetMomentum().GetMagnitude());
// Use momentum magnitude to identify best mc particle
- edm4hep::MCParticle *pBestMCParticle = NULL;
+ //edm4hep::MCParticle *pBestMCParticle = NULL;
+ int best_mc_id = -1;
float bestDeltaMomentum(std::numeric_limits<float>::max());
try
{
@@ -193,15 +195,18 @@ pandora::StatusCode MCParticleCreator::CreateTrackToMCParticleRelationships(cons
const float deltaMomentum(std::fabs(recoMomentum - trueMomentum));
if (deltaMomentum < bestDeltaMomentum)
{
- pBestMCParticle =((*m_id_pMC_map)[ipa.id()]);
+ //pBestMCParticle =((*m_id_pMC_map)[ipa.id()]);
+ best_mc_id = ipa.id() ;
bestDeltaMomentum = deltaMomentum;
}
}
}
}
- if (NULL == pBestMCParticle) continue;
- PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackToMCParticleRelationship(*m_pPandora, pTrack, pBestMCParticle));
+ //if (NULL == pBestMCParticle) continue;
+ if (best_mc_id == -1) continue;
+ //PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackToMCParticleRelationship(*m_pPandora, pTrack, pBestMCParticle));
+ PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::SetTrackToMCParticleRelationship(*m_pPandora, pTrack, (*m_id_pMC_map)[best_mc_id] ) );
}
catch (pandora::StatusCodeException &statusCodeException)
{
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
index ae33a07c79a7bd768fb6dcd1dc6e26e5de11d49a..940d116f14d5682f6adfdd530deba7777a0c9d0d 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
@@ -291,20 +291,20 @@ pandora::StatusCode TrackCreator::ExtractKinks(const CollectionMaps& collectionM
if(collectionMaps.collectionMap_Vertex.find(*iter) == collectionMaps.collectionMap_Vertex.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;}
try
{
- auto pKinkCollection = (collectionMaps.collectionMap_Vertex.find(*iter))->second;
+ auto & pKinkCollection = (collectionMaps.collectionMap_Vertex.find(*iter))->second;
for (int i = 0, iMax = pKinkCollection.size(); i < iMax; ++i)
{
try
{
- auto pVertex0 = pKinkCollection.at(i);
+ auto & pVertex0 = pKinkCollection.at(i);
auto pVertex = &(pVertex0);
if (NULL == pVertex) throw ("Collection type mismatch");
//std::cout<<"pVertex getChi2="<<pVertex->getChi2()<<std::endl;
//std::cout<<"pReconstructedParticle en="<<pReconstructedParticle.getEnergy()<<",type="<<pReconstructedParticle.getType()<<std::endl;
- auto pReconstructedParticle = pVertex->getAssociatedParticle();
+ auto & pReconstructedParticle = pVertex->getAssociatedParticle();
if (this->IsConflictingRelationship(pReconstructedParticle))continue;
const int vertexPdgCode(pReconstructedParticle.getType());
@@ -395,17 +395,17 @@ pandora::StatusCode TrackCreator::ExtractProngsAndSplits(const CollectionMaps& c
if(collectionMaps.collectionMap_Vertex.find(*iter) == collectionMaps.collectionMap_Vertex.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;}
try
{
- auto pProngOrSplitCollection = (collectionMaps.collectionMap_Vertex.find(*iter))->second;
+ auto & pProngOrSplitCollection = (collectionMaps.collectionMap_Vertex.find(*iter))->second;
for (int i = 0, iMax = pProngOrSplitCollection.size(); i < iMax; ++i)
{
try
{
- auto pVertex0 = pProngOrSplitCollection.at(i);
+ auto & pVertex0 = pProngOrSplitCollection.at(i);
auto pVertex = &(pVertex0);
if (NULL == pVertex) throw ("Collection type mismatch");
- const edm4hep::ReconstructedParticle pReconstructedParticle = pVertex->getAssociatedParticle();
+ const edm4hep::ReconstructedParticle & pReconstructedParticle = pVertex->getAssociatedParticle();
if (this->IsConflictingRelationship(pReconstructedParticle))continue;
@@ -461,18 +461,18 @@ pandora::StatusCode TrackCreator::ExtractV0s(const CollectionMaps& collectionMap
if(collectionMaps.collectionMap_Vertex.find(*iter) == collectionMaps.collectionMap_Vertex.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;}
try
{
- auto pV0Collection = (collectionMaps.collectionMap_Vertex.find(*iter))->second;
+ auto & pV0Collection = (collectionMaps.collectionMap_Vertex.find(*iter))->second;
for (int i = 0, iMax = pV0Collection.size(); i < iMax; ++i)
{
try
{
- auto pVertex0 = pV0Collection.at(i);
+ auto & pVertex0 = pV0Collection.at(i);
auto pVertex = &(pVertex0);
if (NULL == pVertex) throw ("Collection type mismatch");
- const edm4hep::ReconstructedParticle pReconstructedParticle = pVertex->getAssociatedParticle();
+ const edm4hep::ReconstructedParticle & pReconstructedParticle = pVertex->getAssociatedParticle();
if (this->IsConflictingRelationship(pReconstructedParticle))continue;
@@ -547,15 +547,15 @@ bool TrackCreator::IsConflictingRelationship(const edm4hep::ReconstructedParticl
return false;
}
-edm4hep::Track* TrackCreator::GetTrackAddress(const CollectionMaps& collectionMaps, const edm4hep::Track& pTrack )
+const edm4hep::Track* TrackCreator::GetTrackAddress(const CollectionMaps& collectionMaps, const edm4hep::Track& pTrack )
{
for (StringVector::const_iterator iter = m_settings.m_trackCollections.begin(), iterEnd = m_settings.m_trackCollections.end(); iter != iterEnd; ++iter)
{
if(collectionMaps.collectionMap_Track.find(*iter) == collectionMaps.collectionMap_Track.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;}
- auto pTrackCollection = (collectionMaps.collectionMap_Track.find(*iter))->second;
+ auto & pTrackCollection = (collectionMaps.collectionMap_Track.find(*iter))->second;
for (int i = 0, iMax = pTrackCollection.size(); i < iMax; ++i)
{
- auto pTrack0 = pTrackCollection.at(i);
+ auto & pTrack0 = pTrackCollection.at(i);
if (pTrack.id() == pTrack0.id()) return (&pTrack0);
}
}
@@ -572,15 +572,14 @@ pandora::StatusCode TrackCreator::CreateTracks(const CollectionMaps& collectionM
if(collectionMaps.collectionMap_Track.find(*iter) == collectionMaps.collectionMap_Track.end()) { std::cout<<"not find "<<(*iter)<<std::endl; continue;}
try
{
- auto pTrackCollection = (collectionMaps.collectionMap_Track.find(*iter))->second;
+ auto & pTrackCollection = (collectionMaps.collectionMap_Track.find(*iter))->second;
if(m_settings.m_debug) std::cout<<"TrackSize:"<<pTrackCollection.size()<<std::endl;
for (int i = 0, iMax = pTrackCollection.size(); i < iMax; ++i)
{
try
{
- auto pTrack0 = pTrackCollection.at(i);
- auto pTrack = (&pTrack0);
+ auto pTrack = &(pTrackCollection.at(i));
if (NULL == pTrack) throw ("Collection type mismatch");
if(m_settings.m_debug){
@@ -665,7 +664,7 @@ pandora::StatusCode TrackCreator::CreateTracks(const CollectionMaps& collectionM
//------------------------------------------------------------------------------------------------------------------------------------------
-void TrackCreator::GetTrackStates(edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const
+void TrackCreator::GetTrackStates(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const
{
// for DD4HEP, 0 is IP, 1 is AtFirstHit, 2 is AtLastHit, 3 is AtCalo
edm4hep::TrackState pTrackState = pTrack->getTrackStates(0);
@@ -691,7 +690,7 @@ void TrackCreator::GetTrackStates(edm4hep::Track *const pTrack, PandoraApi::Trac
//------------------------------------------------------------------------------------------------------------------------------------------
-float TrackCreator::CalculateTrackTimeAtCalorimeter(edm4hep::Track *const pTrack) const
+float TrackCreator::CalculateTrackTimeAtCalorimeter(const edm4hep::Track *const pTrack) const
{
const pandora::Helix helix(pTrack->getTrackStates(0).phi, pTrack->getTrackStates(0).D0, pTrack->getTrackStates(0).Z0, pTrack->getTrackStates(0).omega, pTrack->getTrackStates(0).tanLambda, m_bField);
const pandora::CartesianVector &referencePoint(helix.GetReferencePoint());
@@ -764,7 +763,7 @@ void TrackCreator::CopyTrackState(const edm4hep::TrackState & pTrackState, pando
//------------------------------------------------------------------------------------------------------------------------------------------
-void TrackCreator::TrackReachesECAL(edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const
+void TrackCreator::TrackReachesECAL(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const
{
// Calculate hit position information
@@ -847,7 +846,7 @@ void TrackCreator::TrackReachesECAL(edm4hep::Track *const pTrack, PandoraApi::Tr
//------------------------------------------------------------------------------------------------------------------------------------------
-void TrackCreator::DefineTrackPfoUsage(edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const
+void TrackCreator::DefineTrackPfoUsage(const edm4hep::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const
{
bool canFormPfo(false);
bool canFormClusterlessPfo(false);
@@ -932,7 +931,7 @@ void TrackCreator::DefineTrackPfoUsage(edm4hep::Track *const pTrack, PandoraApi:
//------------------------------------------------------------------------------------------------------------------------------------------
-bool TrackCreator::PassesQualityCuts(edm4hep::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const
+bool TrackCreator::PassesQualityCuts(const edm4hep::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const
{
// First simple sanity checks
if (trackParameters.m_trackStateAtCalorimeter.Get().GetPosition().GetMagnitude() < m_settings.m_minTrackECalDistanceFromIp)
@@ -1009,7 +1008,7 @@ bool TrackCreator::PassesQualityCuts(edm4hep::Track *const pTrack, const Pandora
//------------------------------------------------------------------------------------------------------------------------------------------
-int TrackCreator::GetNTpcHits(edm4hep::Track *const pTrack) const
+int TrackCreator::GetNTpcHits(const edm4hep::Track *const pTrack) const
{
// ATTN
//fg: hit numbers are now given in different order wrt LOI:
@@ -1028,7 +1027,7 @@ int TrackCreator::GetNTpcHits(edm4hep::Track *const pTrack) const
//------------------------------------------------------------------------------------------------------------------------------------------
-int TrackCreator::GetNFtdHits(edm4hep::Track *const pTrack) const
+int TrackCreator::GetNFtdHits(const edm4hep::Track *const pTrack) const
{
// ATTN
//fg: hit numbers are now given in different order wrt LOI:
diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp
index e793bdbe8ac98c691ce7a6e085964d41709456f7..81c3611f63eafe0e0eb537b127b4cb44c0045a74 100644
--- a/Service/GearSvc/src/GearSvc.cpp
+++ b/Service/GearSvc/src/GearSvc.cpp
@@ -206,9 +206,9 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){
double phi0=0;
helpLayer thisLadder;
double beryllium_ladder_block_length=0,end_electronics_half_z=0,side_band_electronics_width=0;
- double rAlu=0, drAlu, rSty, drSty, dzSty, rInner, aluEndcapZ, aluHalfZ, alu_RadLen, Cryostat_dEdx;
- double VXDSupportDensity, VXDSupportZeff, VXDSupportAeff, VXDSupportRadLen, VXDSupportIntLen=0;
- double styDensity, styZeff, styAeff, styRadLen, styIntLen;
+ double rAlu=0, drAlu=0, rSty=0, drSty=0, dzSty=0, rInner=0, aluEndcapZ=0, aluHalfZ=0, alu_RadLen=0, Cryostat_dEdx=0;
+ double VXDSupportDensity=0, VXDSupportZeff=0, VXDSupportAeff=0, VXDSupportRadLen=0, VXDSupportIntLen=0;
+ double styDensity=0, styZeff=0, styAeff=0, styRadLen=0, styIntLen=0;
dd4hep::Volume vxd_vol = vxd.volume();
for(int i=0;i<vxd_vol->GetNdaughters();i++){
TGeoNode* daughter = vxd_vol->GetNode(i);
@@ -402,11 +402,11 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){
int nFlexCable = 0, nFoamSpacer=0, nMetalTraces=0;
int currentLayer = -1;
- double tFlexCable, tFoamSpacer, tMetalTraces;
- double radLFlexCable, radLFoamSpacer, radLMetalTraces;
- double intLFlexCable, intLFoamSpacer, intLMetalTraces;
- double dFlexCable, dFoamSpacer, dMetalTraces;
- double metalZeff, metalZAeff, foamZeff, foamZAeff, flexZeff, flexZAeff;
+ double tFlexCable=0, tFoamSpacer=0, tMetalTraces=0;
+ double radLFlexCable=0, radLFoamSpacer=0, radLMetalTraces=0;
+ double intLFlexCable=0, intLFoamSpacer=0, intLMetalTraces=0;
+ double dFlexCable=0, dFoamSpacer=0, dMetalTraces=0;
+ double metalZeff=0, metalZAeff=0, foamZeff=0, foamZAeff=0, flexZeff=0, flexZAeff=0;
for(int i=0;i<vol->GetNdaughters();i++){
TGeoNode* daughter = vol->GetNode(i);
TGeoMaterial* matDaughter = daughter->GetMedium()->GetMaterial();
diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp
index 53d9673eabf4eaf6e8746bd3d3f4c8b8110ef67b..fb2c1126b79293482e3c66c6ebb3ffe9681302d7 100644
--- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp
+++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.cpp
@@ -138,7 +138,16 @@ AnExampleDetElemTool::ConstructSDandField() {
warning() << "TimeProjectionChamberSensDetTool is not found, and default tracker SD will be used" << endmsg;
}
}
-
+ else {
+ m_tracker_sdtool = ToolHandle<ISensDetTool>("GenericTrackerSensDetTool");
+ if (m_tracker_sdtool) {
+ info() << "Find the GenericTrackerSensDetTool" << endmsg;
+ g4sd = m_tracker_sdtool->createSD(nam);
+ }
+ else {
+ warning() << "GenericTrackerSensDetTool is not found. " << endmsg;
+ }
+ }
}
}
diff --git a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h
index 9606d5cc22620f9fe91d33a24943f6f1256d393d..1b3a923f1f317a3c5e35b9d5d54eb30c84a0b097 100644
--- a/Simulation/DetSimGeom/src/AnExampleDetElemTool.h
+++ b/Simulation/DetSimGeom/src/AnExampleDetElemTool.h
@@ -37,6 +37,7 @@ private:
ToolHandle<ISensDetTool> m_calo_sdtool;
ToolHandle<ISensDetTool> m_driftchamber_sdtool;
ToolHandle<ISensDetTool> m_tpc_sdtool;
+ ToolHandle<ISensDetTool> m_tracker_sdtool;
};
#endif
diff --git a/Simulation/DetSimSD/CMakeLists.txt b/Simulation/DetSimSD/CMakeLists.txt
index eed28a788b4893b6a771aa965f15759b3f2addfa..fa8ba9ef3e13d27da31d8c25763278ba795ef47c 100644
--- a/Simulation/DetSimSD/CMakeLists.txt
+++ b/Simulation/DetSimSD/CMakeLists.txt
@@ -13,6 +13,8 @@ gaudi_add_module(DetSimSD
src/TimeProjectionChamberSensDetTool.cpp
src/TimeProjectionChamberSensitiveDetector.cpp
+ src/GenericTrackerSensDetTool.cpp
+ src/GenericTrackerSensitiveDetector.cpp
LINK DetSimInterface
DetInterface
${DD4hep_COMPONENT_LIBRARIES}
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp b/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..480f509aa57c563e90c032aac64abcaca4957a31
--- /dev/null
+++ b/Simulation/DetSimSD/src/GenericTrackerSensDetTool.cpp
@@ -0,0 +1,40 @@
+#include "GenericTrackerSensDetTool.h"
+
+#include "G4VSensitiveDetector.hh"
+
+#include "DD4hep/Detector.h"
+
+#include "GenericTrackerSensitiveDetector.h"
+
+#include "CLHEP/Units/SystemOfUnits.h"
+
+DECLARE_COMPONENT(GenericTrackerSensDetTool)
+
+StatusCode GenericTrackerSensDetTool::initialize() {
+ StatusCode sc;
+ debug() << "initialize() " << endmsg;
+
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if (!m_geosvc) {
+ error() << "Failed to find GeomSvc." << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ return AlgTool::initialize();
+}
+
+StatusCode GenericTrackerSensDetTool::finalize() {
+ StatusCode sc;
+
+ return sc;
+}
+
+G4VSensitiveDetector* GenericTrackerSensDetTool::createSD(const std::string& name) {
+ debug() << "createSD for " << name << endmsg;
+
+ dd4hep::Detector* dd4hep_geo = m_geosvc->lcdd();
+
+ G4VSensitiveDetector* sd = new GenericTrackerSensitiveDetector(name, *dd4hep_geo);
+
+ return sd;
+}
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensDetTool.h b/Simulation/DetSimSD/src/GenericTrackerSensDetTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..9f46d38111a8eeb6f56517b16474e4909d0aef4f
--- /dev/null
+++ b/Simulation/DetSimSD/src/GenericTrackerSensDetTool.h
@@ -0,0 +1,35 @@
+#ifndef GenericTrackerSensDetTool_h
+#define GenericTrackerSensDetTool_h
+
+/*
+ * GenericTrackerSensDetTool is used to create Time Projection Chamber SD.
+ */
+
+#include "GaudiKernel/AlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "DetSimInterface/ISensDetTool.h"
+#include "DetInterface/IGeomSvc.h"
+
+#include "DD4hep/DD4hepUnits.h"
+
+class GenericTrackerSensDetTool: public extends<AlgTool, ISensDetTool> {
+
+ public:
+
+ using extends::extends;
+
+ /// Overriding initialize and finalize
+ StatusCode initialize() override;
+ StatusCode finalize() override;
+
+ /// Override ISensDetTool
+ virtual G4VSensitiveDetector* createSD(const std::string& name) override;
+
+private:
+
+ // in order to initialize SD, we need to get the lcdd()
+ SmartIF<IGeomSvc> m_geosvc;
+
+};
+
+#endif
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d250079ef0a8aa5f15e50fb801b1de90c46d3a17
--- /dev/null
+++ b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.cpp
@@ -0,0 +1,59 @@
+#include "GenericTrackerSensitiveDetector.h"
+
+#include "G4Step.hh"
+#include "G4VProcess.hh"
+//#include "G4SteppingManager.hh"
+#include "G4SDManager.hh"
+//#include "UserTrackInformation.hh"
+#include "DD4hep/DD4hepUnits.h"
+
+GenericTrackerSensitiveDetector::GenericTrackerSensitiveDetector(const std::string& name,
+ dd4hep::Detector& description)
+ : DDG4SensitiveDetector(name, description),
+ m_hc(nullptr){
+
+ G4String CollName=name+"Collection";
+ collectionName.insert(CollName);
+}
+
+void GenericTrackerSensitiveDetector::Initialize(G4HCofThisEvent* HCE){
+ m_hc = new HitCollection(GetName(), collectionName[0]);
+ int HCID = G4SDManager::GetSDMpointer()->GetCollectionID(m_hc);
+ HCE->AddHitsCollection( HCID, m_hc );
+}
+
+G4bool GenericTrackerSensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*){
+
+ G4TouchableHandle touchPost = step->GetPostStepPoint()->GetTouchableHandle();
+ G4TouchableHandle touchPre = step->GetPreStepPoint()->GetTouchableHandle();
+ dd4hep::sim::Geant4StepHandler h(step);
+ dd4hep::Position prePos = h.prePos();
+ dd4hep::Position postPos = h.postPos();
+ dd4hep::Position direction = postPos - prePos;
+ dd4hep::Position position = mean_direction(prePos,postPos);
+ double hit_len = direction.R();
+ if (hit_len > 0) {
+ double new_len = mean_length(h.preMom(),h.postMom())/hit_len;
+ direction *= new_len/hit_len;
+ }
+ dd4hep::sim::Geant4TrackerHit* hit = nullptr;
+ hit = new dd4hep::sim::Geant4TrackerHit(h.track->GetTrackID(),
+ h.track->GetDefinition()->GetPDGEncoding(),
+ step->GetTotalEnergyDeposit(),
+ h.track->GetGlobalTime());
+
+ if ( hit ) {
+ hit->cellID = getCellID( step ) ;
+ hit->energyDeposit = step->GetTotalEnergyDeposit();
+ hit->position = position;
+ hit->momentum = direction;
+ hit->length = hit_len;
+ m_hc->insert(hit);
+ return true;
+ }
+ throw std::runtime_error("new() failed: Cannot allocate hit object");
+ return false;
+}
+
+void GenericTrackerSensitiveDetector::EndOfEvent(G4HCofThisEvent* HCE){
+}
diff --git a/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.h b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.h
new file mode 100644
index 0000000000000000000000000000000000000000..63257295e9cf8f4f580fbdcb246e7734d33c64ee
--- /dev/null
+++ b/Simulation/DetSimSD/src/GenericTrackerSensitiveDetector.h
@@ -0,0 +1,24 @@
+// *********************************************************
+//
+// $Id: GenericTrackerSensitiveDetector.hh,v 1.0 2022/03/27
+
+#ifndef GenericTrackerSensitiveDetector_h
+#define GenericTrackerSensitiveDetector_h
+
+#include "DetSimSD/DDG4SensitiveDetector.h"
+#include "DDG4/Defs.h"
+
+class GenericTrackerSensitiveDetector: public DDG4SensitiveDetector {
+ public:
+ GenericTrackerSensitiveDetector(const std::string& name, dd4hep::Detector& description);
+
+ void Initialize(G4HCofThisEvent* HCE);
+ G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
+ void EndOfEvent(G4HCofThisEvent* HCE);
+
+ protected:
+
+ HitCollection* m_hc = nullptr;
+
+};
+#endif
diff --git a/Utilities/KalDet/src/ild/common/MaterialDataBase.cc b/Utilities/KalDet/src/ild/common/MaterialDataBase.cc
index a6b8f6e63488dd7026ad2a2f958ecf8292374966..f181b8b53de84a990f3b4351d455241e0a344bf8 100644
--- a/Utilities/KalDet/src/ild/common/MaterialDataBase.cc
+++ b/Utilities/KalDet/src/ild/common/MaterialDataBase.cc
@@ -205,9 +205,12 @@ void MaterialDataBase::createMaterials(const gear::GearMgr& gearMgr, IGeomSvc* g
// VXD Support Material
if(geoSvc){
- TMaterial* vxdsupport = geoSvc->getMaterial("VXDSupportMaterial");
- if(vxdsupport) this->addMaterial(vxdsupport, "VXDSupportMaterial");
- else std::cout << "Material VXDSupportMaterial not found" << std::endl;
+ //obselete
+ //TMaterial* vxdsupport = geoSvc->getMaterial("VXDSupportMaterial");
+ //if(vxdsupport) this->addMaterial(vxdsupport, "VXDSupportMaterial");
+ //else std::cout << "Material VXDSupportMaterial not found" << std::endl;
+ std::cout << "geoSvc should be set as NULL before new GeomSvc interface ready!" << std::endl;
+ exit(1);
}
else{
try{
diff --git a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc
index fe92415eb745cdfcbf883be2fd2c9ac3dc878bdb..64dd65463d1cf178c0c85e7d87bd079390591e21 100644
--- a/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc
+++ b/Utilities/KalDet/src/ild/support/ILDSupportKalDetector.cc
@@ -35,7 +35,22 @@ TVKalDetector(10)
std::vector<double> z, rInner, rOuter;
// streamlog_out(DEBUG1) << "ILDSupportKalDetector building beampipe using GEAR " << std::endl ;
if(geoSvc){
- const dd4hep::rec::ConicalSupportData* pBeamPipeData = geoSvc->getBeamPipeData();
+ dd4hep::DetElement world = geoSvc->getDD4HepGeo();
+ dd4hep::DetElement pipe;
+ const std::map<std::string, dd4hep::DetElement>& subs = world.children();
+ for(std::map<std::string, dd4hep::DetElement>::const_iterator it=subs.begin();it!=subs.end();it++){
+ if(it->first!="Tube") continue;
+ pipe = it->second;
+ }
+ dd4hep::rec::ConicalSupportData* pBeamPipeData = nullptr;
+ try{
+ pBeamPipeData = pipe.extension<dd4hep::rec::ConicalSupportData>();
+ }
+ catch(std::runtime_error& e){
+ std::cout << e.what() << " " << pBeamPipeData << std::endl;
+ std::cerr << "Cannot find extension data in Beampipe!" << std::endl;
+ }
+
const std::vector<dd4hep::rec::ConicalSupportData::Section>& sections = pBeamPipeData->sections;
const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
bz = field.z()/dd4hep::tesla;
diff --git a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc
index 5cacd3d76205a515295102ec6c6f865bd24776a4..3043a7f392acd79700a441c826b3cbd3f5f32522 100644
--- a/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc
+++ b/Utilities/KalDet/src/ild/vxd/ILDVXDKalDetector.cc
@@ -382,7 +382,6 @@ void ILDVXDKalDetector::setupGearGeom( const gear::GearMgr& gearMgr ){
void ILDVXDKalDetector::setupGearGeom( IGeomSvc* geoSvc){
- /*
dd4hep::DetElement world = geoSvc->getDD4HepGeo();
dd4hep::DetElement vxd;
const std::map<std::string, dd4hep::DetElement>& subs = world.children();
@@ -398,14 +397,13 @@ void ILDVXDKalDetector::setupGearGeom( IGeomSvc* geoSvc){
std::cout << e.what() << " " << vxdData << std::endl;
throw GaudiException(e.what(), "FATAL", StatusCode::FAILURE);
}
- */
+
const dd4hep::Direction& field = geoSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
_bZ = field.z()/dd4hep::tesla;
+
+ const std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& layers = vxdData->layers;
- const gear::ZPlanarParametersImpl* pVXDDetMain = geoSvc->getVXDParameters();
- const gear::VXDLayerLayout& pVXDLayerLayout = pVXDDetMain->getVXDLayerLayout();
-
- _nLayers = pVXDLayerLayout.getNLayers();
+ _nLayers = layers.size();
_VXDgeo.resize(_nLayers);
//SJA:FIXME: for now the support is taken as the same size the sensitive
@@ -414,16 +412,17 @@ void ILDVXDKalDetector::setupGearGeom( IGeomSvc* geoSvc){
// for a significant distance
for( int layer=0; layer < _nLayers; ++layer){
- _VXDgeo[layer].nLadders = pVXDLayerLayout.getNLadders(layer);
- _VXDgeo[layer].phi0 = pVXDLayerLayout.getPhi0(layer);
+ const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = layers[layer];
+ _VXDgeo[layer].nLadders = thisLayer.ladderNumber;
+ _VXDgeo[layer].phi0 = thisLayer.phi0;
_VXDgeo[layer].dphi = 2*M_PI / _VXDgeo[layer].nLadders;
- _VXDgeo[layer].senRMin = pVXDLayerLayout.getSensitiveDistance(layer);
- _VXDgeo[layer].supRMin = pVXDLayerLayout.getLadderDistance(layer);
- _VXDgeo[layer].length = pVXDLayerLayout.getSensitiveLength(layer) * 2.0 ; // note: gear for historical reasons uses the halflength
- _VXDgeo[layer].width = pVXDLayerLayout.getSensitiveWidth(layer);
- _VXDgeo[layer].offset = pVXDLayerLayout.getSensitiveOffset(layer);
- _VXDgeo[layer].senThickness = pVXDLayerLayout.getSensitiveThickness(layer);
- _VXDgeo[layer].supThickness = pVXDLayerLayout.getLadderThickness(layer);
+ _VXDgeo[layer].senRMin = thisLayer.distanceSensitive;
+ _VXDgeo[layer].supRMin = thisLayer.distanceSupport;
+ _VXDgeo[layer].length = thisLayer.zHalfSensitive * 2.0 ; // note: gear for historical reasons uses the halflength
+ _VXDgeo[layer].width = thisLayer.widthSensitive;
+ _VXDgeo[layer].offset = thisLayer.offsetSensitive;
+ _VXDgeo[layer].senThickness = thisLayer.thicknessSensitive;
+ _VXDgeo[layer].supThickness = thisLayer.thicknessSupport;
//std::cout << layer << ": " << _VXDgeo[layer].nLadders << " " << _VXDgeo[layer].phi0 << " " << _VXDgeo[layer].dphi << " " << _VXDgeo[layer].senRMin
// << " " << _VXDgeo[layer].supRMin << " " << _VXDgeo[layer].length << " " << _VXDgeo[layer].width << " " << _VXDgeo[layer].offset
// << " " << _VXDgeo[layer].senThickness << " " << _VXDgeo[layer].supThickness << std::endl;
@@ -432,26 +431,20 @@ void ILDVXDKalDetector::setupGearGeom( IGeomSvc* geoSvc){
// layer, this can optionally be changed, e.g. in the case of the FPCCD where the
// epitaxial layer is 15 mu thick (in a 50 mu wafer)
_relative_position_of_measurement_surface = 0.5 ;
+
// Cryostat
- try {
- _vxd_Cryostat.alRadius = geoSvc->getDetParameter("VXDInfra","CryostatAlRadius");
- _vxd_Cryostat.alThickness = geoSvc->getDetParameter("VXDInfra","CryostatAlThickness");
- _vxd_Cryostat.alInnerR = geoSvc->getDetParameter("VXDInfra","CryostatAlInnerR");
- _vxd_Cryostat.alZEndCap = geoSvc->getDetParameter("VXDInfra","CryostatAlZEndCap");
- _vxd_Cryostat.alHalfZ = geoSvc->getDetParameter("VXDInfra","CryostatAlHalfZ");
+ // FIXME not support by ZPlanarData, not install in current CEPC
+ //_vxd_Cryostat.shellInnerR = vxdData->rInnerShell;
+ //_vxd_Cryostat.shellThickness = vxdData->rOuterShell;
+ //_vxd_Cryostat.shelllHalfZ = vxdData->zHalfShell;
+ //_vxd_Cryostat.alRadius = geoSvc->getDetParameter("VXDInfra","CryostatAlRadius");
+ //_vxd_Cryostat.alThickness = geoSvc->getDetParameter("VXDInfra","CryostatAlThickness");
+ //_vxd_Cryostat.alInnerR = geoSvc->getDetParameter("VXDInfra","CryostatAlInnerR");
+ //_vxd_Cryostat.alZEndCap = geoSvc->getDetParameter("VXDInfra","CryostatAlZEndCap");
+ //_vxd_Cryostat.alHalfZ = geoSvc->getDetParameter("VXDInfra","CryostatAlHalfZ");
+ _vxd_Cryostat.exists = false;
- _vxd_Cryostat.shellInnerR = pVXDDetMain->getShellInnerRadius();
- _vxd_Cryostat.shellThickness = pVXDDetMain->getShellOuterRadius() - _vxd_Cryostat.shellInnerR;
- _vxd_Cryostat.shelllHalfZ = pVXDDetMain->getShellHalfLength();
-
- _vxd_Cryostat.exists = true;
- //std::cout << "VXDInfra: " << _vxd_Cryostat.alRadius << " " << _vxd_Cryostat.alThickness << " " << _vxd_Cryostat.alInnerR << " " << _vxd_Cryostat.alZEndCap << " "
- // << _vxd_Cryostat.alHalfZ << " " << _vxd_Cryostat.shellInnerR << " " << _vxd_Cryostat.shellThickness << " " << _vxd_Cryostat.shelllHalfZ << std::endl;
- }
- catch (std::runtime_error& e) {
- std::cout << e.what() << std::endl ;
- _vxd_Cryostat.exists = false;
-
- }
+ //std::cout << "VXDInfra: " << _vxd_Cryostat.alRadius << " " << _vxd_Cryostat.alThickness << " " << _vxd_Cryostat.alInnerR << " " << _vxd_Cryostat.alZEndCap << " "
+ // << _vxd_Cryostat.alHalfZ << " " << _vxd_Cryostat.shellInnerR << " " << _vxd_Cryostat.shellThickness << " " << _vxd_Cryostat.shelllHalfZ << std::endl;
}