Changes for pandora PFA

FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml

@@ -48,7 +48,7 @@
   <include ref="HCalBarrel_TileCal_v03.xml"/>
   <include ref="ECalEndcaps_Turbine.xml"/>
   <include ref="HCalEndcaps_ThreeParts_TileCal_v03.xml"/>
-  <include ref="MuonTagger.xml"/>
+  <include ref="MuonTaggerPhiTheta.xml"/>
 
   <fields>
     <!-- FIXME this is a place holder to get a reasonable magnetic field, it should be a real MAP obtained from FEM tools taking into account HCAL return yoke -->

FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/MuonTaggerPhiTheta.xml

@@ -0,0 +1,66 @@
+<?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="FCCeeMuonTagger"
+        title="Muon Tagger for ALLEGRO, with phi-theta segmentation"
+        author="Created by A. Duglishvili from MuonTagger.xml, moving from eta to theta segmentation"
+        url="no"
+        status="development"
+        version="1.0">
+    <comment>
+      Simple muon tagger with phi-theta segmentation - barrel and endcaps
+    </comment>
+  </info>
+
+  <define>
+    <constant name="MuonTaggerBarrelLayers" value="2"/>
+    <constant name="MuonTaggerEndcapLayers" value="2"/>
+  </define>
+
+  <display>
+    <vis name="muon_vis" r="1.0" g="0.1" b="0.1" showDaughters="true" visible="true" />
+  </display> 
+
+  <readouts>
+    <readout name="MuonTaggerBarrelPhiTheta">
+      <segmentation type="FCCSWGridPhiTheta_k4geo" phi_bins="704" offset_phi="-pi+(pi/704.)" grid_size_theta="0.010013373" offset_theta="0.10884388"/>
+      <id>system:4,subsystem:1,layer:5,theta:10,phi:10</id>
+    </readout>
+    <readout name="MuonTaggerEndcapPhiTheta">
+      <segmentation type="FCCSWGridPhiTheta_k4geo" phi_bins="704" offset_phi="-pi+(pi/704.)" grid_size_theta="0.010013373" offset_theta="0.10884388"/>
+      <id>system:4,subsystem:1,layer:5,theta:10,phi:10</id>
+    </readout>
+  </readouts>
+
+  <detectors>
+    <!-- barrel -->
+    <detector id="DetID_Muon_Barrel" name="MuonTaggerBarrel" type="SimpleCylinder_o1_v03" sensitive="true" vis="muon_vis" readout="MuonTaggerBarrelPhiTheta">
+      <!-- added for Pandora -->
+      <type_flags type="DetType_CALORIMETER + DetType_MUON + DetType_BARREL"/>
+      <!-- end -->
+      <sensitive type="SimpleTrackerSD"/>
+      <dimensions rmin="MuonTagger_inner_radius" rmax="MuonTagger_outer_radius"
+                  dz="MuonTagger_half_length" z_offset = "0" material="Polystyrene" phi0="0" deltaphi="360*deg" vis="muon_vis"/>
+      <layers>
+        <layer thickness="(MuonTagger_outer_radius - MuonTagger_inner_radius)/MuonTaggerBarrelLayers" repeat="MuonTaggerBarrelLayers"/>
+      </layers>
+    </detector>
+
+    <!-- endcaps -->
+    <detector id="DetID_Muon_Endcap_1" name="MuonTaggerEndcap" type="SimpleCylinder_o1_v03" sensitive="true" vis="muon_vis" readout="MuonTaggerEndcapPhiTheta">
+      <!-- added for Pandora -->
+      <type_flags type="DetType_CALORIMETER + DetType_MUON + DetType_ENDCAP"/>
+      <!-- end -->
+      <sensitive type="SimpleTrackerSD"/>
+      <dimensions rmin="MuonTaggerEndcap_inner_radius" rmax="MuonTaggerEndcap_outer_radius"
+                  dz="(MuonTaggerEndcap_max_z - MuonTaggerEndcap_min_z)*0.5" z_offset = "MuonTaggerEndcap_min_z + (MuonTaggerEndcap_max_z - MuonTaggerEndcap_min_z)*0.5" material="Polystyrene" phi0="0" deltaphi="360*deg" vis="muon_vis"/>
+      <layers>
+        <layer thickness="(MuonTaggerEndcap_max_z - MuonTaggerEndcap_min_z)/MuonTaggerEndcapLayers" repeat="MuonTaggerEndcapLayers"/>
+      </layers>
+    </detector>
+
+  </detectors>
+
+</lccdd>

FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/ALLEGRO_o1_v04.xml

@@ -48,7 +48,7 @@
   <include ref="HCalBarrel_TileCal_v03.xml"/>
   <include ref="ECalEndcaps_Turbine.xml"/>
   <include ref="HCalEndcaps_ThreeParts_TileCal_v03.xml"/>
-  <include ref="MuonTagger.xml"/>
+  <include ref="MuonTaggerPhiTheta.xml"/>
 
   <fields>
     <!-- FIXME this is a place holder to get a reasonable magnetic field, it should be a real MAP obtained from FEM tools taking into account HCAL return yoke -->

FCCee/ALLEGRO/compact/ALLEGRO_o1_v04/MuonTaggerPhiTheta.xml

@@ -0,0 +1 @@
+../ALLEGRO_o1_v03/MuonTaggerPhiTheta.xml

detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v03_geo.cpp

@@ -683,14 +683,15 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
   envelopePhysVol.addPhysVolID("system", xmlDetElem.id());
   caloDetElem.setPlacement(envelopePhysVol);
 
-  // Create caloData object
+  // Create caloData object for the reconstruction
   auto caloData = new dd4hep::rec::LayeredCalorimeterData;
   caloData->layoutType = dd4hep::rec::LayeredCalorimeterData::BarrelLayout;
   caloDetElem.addExtension<dd4hep::rec::LayeredCalorimeterData>(caloData);
-
+  // Extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm
+  // (for barrel detectors zmin is 0)
   caloData->extent[0] = Rmin;
-  caloData->extent[1] = Rmax; // or r_max ?
-  caloData->extent[2] = 0.;      // NN: for barrel detectors this is 0
+  caloData->extent[1] = Rmax;
+  caloData->extent[2] = 0.;
   caloData->extent[3] = caloDim.dz();
 
 

detector/calorimeter/ECalBarrel_NobleLiquid_InclinedTrapezoids_o1_v04_geo.cpp

@@ -753,21 +753,36 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
   envelopePhysVol.addPhysVolID("system", xmlDetElem.id());
   caloDetElem.setPlacement(envelopePhysVol);
 
-  // Create caloData object
+  // Create caloData object for the reconstruction
   auto caloData = new dd4hep::rec::LayeredCalorimeterData;
   caloData->layoutType = dd4hep::rec::LayeredCalorimeterData::BarrelLayout;
   caloDetElem.addExtension<dd4hep::rec::LayeredCalorimeterData>(caloData);
+  dd4hep::xml::setDetectorTypeFlag(xmlDetElem, caloDetElem);
+
+  // Fill caloData information
 
+  // Extent of the calorimeter in the r-z-plane [ rmin, rmax, zmin, zmax ] in mm
+  // (for barrel detectors zmin is 0)
   caloData->extent[0] = Rmin;
-  caloData->extent[1] = Rmax; // or r_max ?
-  caloData->extent[2] = 0.;      // NN: for barrel detectors this is 0
+  caloData->extent[1] = Rmax;
+  caloData->extent[2] = 0.;
   caloData->extent[3] = caloDim.dz();
 
-  // Set type flags
-  dd4hep::xml::setDetectorTypeFlag(xmlDetElem, caloDetElem);
+  // Information about each layer
+  // double distance : distance from Origin (or the z-axis) to the inner-most face of the layer
+  // double phi0 : phi0 of layer: potential rotation around normal to absorber plane, e.g. if layers are 'staggered' in phi in fwd. calos
+  // double absorberThickness : thickness of the absorber part of the layer. Consider using inner/outer_nRadiationLengths and inner/outer_nInteractionLengths
+  // double inner_nRadiationLengths : Absorber material in front of sensitive element in the layer, units of radiation lengths
+  // double inner_nInteractionLengths : Absorber material in front of sensitive element in the layer, units of radiation lengths
+  // double outer_nRadiationLengths : Absorber material in behind of sensitive element in the layer, units of radiation lengths
+  // double outer_nInteractionLengths : Absorber material in behind of sensitive element in the layer, units of radiation lengths
+  // double inner_thickness : Distance between the innermost face of the layer (closest to IP) and the center of the sensitive element
+  // double outer_thickness : Distance between the center of the sensitive element and the outermost face of the layer
+  // double sensitive_thickness : Thickness of the sensitive element (e.g. scintillator)
+  // double cellSize0 : cell size along the first axis where first is either along the beam (BarrelLayout) or up (EndcapLayout) or the direction closest to that
+  // double cellSize1 : second cell size, perpendicular to the first direction cellSize0 and the depth of the layers
   dd4hep::rec::MaterialManager matMgr(envelopeVol);
   dd4hep::rec::LayeredCalorimeterData::Layer caloLayer;
-
   double rad_first = Rmin;
   double rad_last = 0;
   double scale_fact = dR / (-Rmin * cos(angle) + sqrt(pow(Rmax, 2) - pow(Rmin * sin(angle), 2)));
@@ -786,19 +801,19 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
     auto mat = matMgr.createAveragedMaterial(materials);                             // creating average of all the material between two points to calculate X0 and lambda of averaged material
     const double nRadiationLengths = mat.radiationLength();
     const double nInteractionLengths = mat.interactionLength();
-    const double difference_bet_r1r2 = (ivr1 - ivr2).r();
-    const double value_of_x0 = layerHeight[il] / nRadiationLengths;
-    const double value_of_lambda = layerHeight[il] / nInteractionLengths;
+    const double difference_bet_r1r2 = (ivr1 - ivr2).r(); // should be identical to layerHeight[il] * scale_fact
+    const double value_of_x0 = layerHeight[il] / nRadiationLengths; // GM : shouldn't this be multiplied by scale_fact?
+    const double value_of_lambda = layerHeight[il] / nInteractionLengths; // GM : shouldn't this be multiplied by scale_fact?
     std::string str1("LAr");
 
     for (size_t imat = 0; imat < materials.size(); imat++) {
 
       std::string str2(materials.at(imat).first.name());
       if (str1.compare(str2) == 0){
-	  thickness_sen += materials.at(imat).second;
+	      thickness_sen += materials.at(imat).second;
       }
       else {
-	absorberThickness += materials.at(imat).second;
+	      absorberThickness += materials.at(imat).second;
       }
     }
     rad_first = rad_last;
@@ -807,7 +822,9 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
     std::cout << "The radiation length is " << value_of_x0 << " and the interaction length is " << value_of_lambda << std::endl;
 
     caloLayer.distance = rad_first;
+    caloLayer.absorberThickness         = absorberThickness;
     caloLayer.sensitive_thickness       = thickness_sen;
+
     caloLayer.inner_nRadiationLengths   = value_of_x0 / 2.0;
     caloLayer.inner_nInteractionLengths = value_of_lambda / 2.0;
     caloLayer.inner_thickness           = difference_bet_r1r2 / 2.0;
@@ -816,9 +833,9 @@ static dd4hep::detail::Ref_t createECalBarrelInclined(dd4hep::Detector& aLcdd,
     caloLayer.outer_nInteractionLengths = value_of_lambda / 2.0;
     caloLayer.outer_thickness           = difference_bet_r1r2 / 2;
 
-    caloLayer.absorberThickness         = absorberThickness;
-    caloLayer.cellSize0 = 20 * dd4hep::mm;
-    caloLayer.cellSize1 = 20 * dd4hep::mm;
+    // GM: rather retrieve cellDimensions vector from segmentation class
+    caloLayer.cellSize0 = 20 * dd4hep::mm; // GM: rather put delta_eta here and use pandora::POINTING cell type (should actually modify pandora to accept theta grid)
+    caloLayer.cellSize1 = 20 * dd4hep::mm; // GM: rather delta_phi here + pandora::POINTING cell type
 
     caloData->layers.push_back(caloLayer);
   }

detector/calorimeter/ECalEndcap_Turbine_o1_v01_geo.cpp

@@ -1,6 +1,8 @@
 #include "DD4hep/DetFactoryHelper.h"
 #include "DD4hep/Printout.h"
 #include "TMatrixT.h"
+#include "XML/Utilities.h"
+#include <DDRec/DetectorData.h>
 
 // todo: remove gaudi logging and properly capture output
 #define endmsg std::endl
@@ -659,6 +661,22 @@ createECalEndcapTurbine(dd4hep::Detector& aLcdd, dd4hep::xml::Handle_t aXmlEleme
   dd4hep::PlacedVolume envelopePhysVol = motherVol.placeVolume(endcapsAssembly);
   caloDetElem.setPlacement(envelopePhysVol);
   envelopePhysVol.addPhysVolID("system", idDet);
+
+  // Create dummy caloData object for PandoraPFA
+  // FIXME: fill calo and layer data information
+  auto caloData = new dd4hep::rec::LayeredCalorimeterData;
+  caloData->layoutType = dd4hep::rec::LayeredCalorimeterData::EndcapLayout;
+  caloDetElem.addExtension<dd4hep::rec::LayeredCalorimeterData>(caloData);
+
+  // save extent information
+  caloData->extent[0] = dim.rmin1();
+  caloData->extent[1] = dim.rmax1();
+  caloData->extent[2] = dim.z_offset()-dim.dz()/2.;
+  caloData->extent[3] = dim.z_offset()+dim.dz()/2.;
+
+  // Set type flags
+  dd4hep::xml::setDetectorTypeFlag(xmlDetElem, caloDetElem);
+
   return caloDetElem;
 }
 }  // namespace det

detector/calorimeter/ECalEndcap_Turbine_o1_v02_geo.cpp

@@ -1,6 +1,8 @@
 #include "DD4hep/DetFactoryHelper.h"
 #include "DD4hep/Printout.h"
 #include "TMatrixT.h"
+#include "XML/Utilities.h"
+#include <DDRec/DetectorData.h>
 
 // todo: remove gaudi logging and properly capture output
 #define endmsg std::endl
@@ -663,6 +665,22 @@ namespace det {
       dd4hep::PlacedVolume envelopePhysVol = motherVol.placeVolume(endcapsAssembly);
       caloDetElem.setPlacement(envelopePhysVol);
       envelopePhysVol.addPhysVolID("system", idDet);
+
+      // Create dummy caloData object for PandoraPFA
+      // FIXME: fill calo and layer data information
+      auto caloData = new dd4hep::rec::LayeredCalorimeterData;
+      caloData->layoutType = dd4hep::rec::LayeredCalorimeterData::EndcapLayout;
+      caloDetElem.addExtension<dd4hep::rec::LayeredCalorimeterData>(caloData);
+
+      // save extent information
+      caloData->extent[0] = dim.rmin1();
+      caloData->extent[1] = dim.rmax1();
+      caloData->extent[2] = dim.z_offset()-dim.dz()/2.;
+      caloData->extent[3] = dim.z_offset()+dim.dz()/2.;
+
+      // Set type flags
+      dd4hep::xml::setDetectorTypeFlag(xmlDetElem, caloDetElem);
+
       return caloDetElem;
     }
   }

detector/calorimeter/HCalThreePartsEndcap_o1_v02_geo.cpp

@@ -2,6 +2,7 @@
 #include "DD4hep/DetFactoryHelper.h"
 #include <DDRec/DetectorData.h>
 #include "DD4hep/Printout.h"
+#include "XML/Utilities.h"
 
 using dd4hep::Volume;
 using dd4hep::DetElement;
@@ -451,10 +452,10 @@ static dd4hep::Ref_t createHCalEC(dd4hep::Detector& lcdd, xml_h xmlElement, dd4h
     caloData->layoutType = dd4hep::rec::LayeredCalorimeterData::EndcapLayout;
     caloDetElem.addExtension<dd4hep::rec::LayeredCalorimeterData>(caloData);
 
-    caloData->extent[0] = sensitiveBarrel3Rmin;
-    caloData->extent[1] = sensitiveBarrel3Rmin + moduleDepth3; // 
-    caloData->extent[2] = 0.;      // NN: for barrel detectors this is 0
-    caloData->extent[3] = dzDetector1 + dzDetector2 + dzDetector3; 
+    caloData->extent[0] = sensitiveBarrel3Rmin; // innerRCoordinate
+    caloData->extent[1] = sensitiveBarrel3Rmin + moduleDepth3; // outerRCoordinate
+    caloData->extent[2] = extBarrelOffset1 - dzDetector1; // innerZCoordinate (start of the first part of the Endcap)
+    caloData->extent[3] = extBarrelOffset3 + dzDetector3; // outerZCoordinate (end of the third part of the Endcap)
 
     dd4hep::rec::LayeredCalorimeterData::Layer caloLayer;
 
@@ -489,6 +490,10 @@ static dd4hep::Ref_t createHCalEC(dd4hep::Detector& lcdd, xml_h xmlElement, dd4h
 
         caloData->layers.push_back(caloLayer);
     }
+
+    // Set type flags
+    dd4hep::xml::setDetectorTypeFlag(xmlDet, caloDetElem);
+
     return caloDetElem; 
 }  
 }// namespace det

detector/calorimeter/HCalTileBarrel_o1_v02_geo.cpp

@@ -2,6 +2,9 @@
 #include "DD4hep/DetFactoryHelper.h"
 #include <DDRec/DetectorData.h>
 #include "DD4hep/Printout.h"
+#include "XML/Utilities.h"
+#include "DDRec/MaterialManager.h"
+#include "DDRec/Vector3D.h"
 
 using dd4hep::Volume;
 using dd4hep::DetElement;
@@ -238,22 +241,59 @@ static dd4hep::Ref_t createHCal(dd4hep::Detector& lcdd, xml_det_t xmlDet, dd4hep
   caloData->extent[2] = 0.;      // NN: for barrel detectors this is 0
   caloData->extent[3] = dzDetector; 
 
+  dd4hep::rec::MaterialManager matMgr(envelopeVolume);
   dd4hep::rec::LayeredCalorimeterData::Layer caloLayer;
 
   for (unsigned int idxLayer = 0; idxLayer < layerDepths.size(); ++idxLayer) {
-        const double difference_bet_r1r2 = layerDepths.at(idxLayer); 
-
-        caloLayer.distance                  = layerInnerRadii.at(idxLayer); // radius of the current layer   
-        caloLayer.sensitive_thickness       = difference_bet_r1r2;  // radial dimension of the current layer 
-        caloLayer.inner_thickness           = difference_bet_r1r2 / 2.0;
-        caloLayer.outer_thickness           = difference_bet_r1r2 / 2.0;
-
-        caloData->layers.push_back(caloLayer);
+    const double difference_bet_r1r2 = layerDepths.at(idxLayer);
+    double thickness_sen = 0.;
+    double absorberThickness = 0.;
+
+    // AD: average material radiation length in a given layer depends on the polar angle... not sure how it is used by Pandora... lets calculate it at the angle of 60 degrees...
+    const double angle = 60.*M_PI/180.;
+    dd4hep::rec::Vector3D ivr1 = dd4hep::rec::Vector3D(0., layerInnerRadii.at(idxLayer), 0); // defining starting vector points of the given layer
+    dd4hep::rec::Vector3D ivr2 = dd4hep::rec::Vector3D(0., layerInnerRadii.at(idxLayer) + layerDepths.at(idxLayer), layerDepths.at(idxLayer)*cos(angle)/sin(angle));  // defining end vector points of the given layer
+
+    const dd4hep::rec::MaterialVec &materials = matMgr.materialsBetween(ivr1, ivr2); // calling material manager to get material info between two points
+    auto mat = matMgr.createAveragedMaterial(materials);                             // creating average of all the material between two points to calculate X0 and lambda of averaged material
+    const double nRadiationLengths = layerDepths.at(idxLayer) / mat.radiationLength();
+    const double nInteractionLengths = layerDepths.at(idxLayer) / mat.interactionLength();
+
+    std::string str1("Polystyrene"); // sensitive material
+    for (size_t imat = 0; imat < materials.size(); imat++) {
+      std::string str2(materials.at(imat).first.name());
+      if (str1.compare(str2) == 0){
+          thickness_sen += materials.at(imat).second;
+      }
+      else if(str2!="Air") {
+        absorberThickness += materials.at(imat).second;
+      }
+    }
+    std::cout << "The sensitive thickness is " << thickness_sen << std::endl;
+    std::cout << "The absorber thickness is " << absorberThickness << std::endl;
+    std::cout << "The number of radiation length is " << nRadiationLengths << " and the number of interaction length is " << nInteractionLengths << std::endl;
+
+    caloLayer.distance                  = layerInnerRadii.at(idxLayer); // radius of the current layer   
+    caloLayer.sensitive_thickness	= difference_bet_r1r2;  // radial dimension of the current layer 
+    //caloLayer.sensitive_thickness	= thickness_sen;
+    caloLayer.absorberThickness         = absorberThickness;
+
+    caloLayer.inner_thickness           = difference_bet_r1r2 / 2.0;
+    caloLayer.inner_nRadiationLengths   = nRadiationLengths / 2.0;
+    caloLayer.inner_nInteractionLengths = nInteractionLengths / 2.0;
+    caloLayer.outer_nRadiationLengths   = nRadiationLengths / 2.0;
+    caloLayer.outer_nInteractionLengths = nInteractionLengths / 2.0;
+    caloLayer.outer_thickness           = difference_bet_r1r2 / 2;
+
+    caloLayer.cellSize0 = 20 * dd4hep::mm; // should be updated from DDGeometryCreatorALLEGRO
+    caloLayer.cellSize1 = 20 * dd4hep::mm; // should be updated from DDGeometryCreatorALLEGRO
+    caloData->layers.push_back(caloLayer);
   }
 
+  // Set type flags
+  dd4hep::xml::setDetectorTypeFlag(xmlDet, caloDetElem);
 
   return caloDetElem;
-
 }
 }  // namespace hcal