Rad dam01

FastSimulation/CaloHitMakers/interface/CaloHitMaker.h

@@ -14,9 +14,11 @@
 //STL headers
 #include <string>
 #include <map>
+#include <vector>
 
 class CaloGeometryHelper;
 class CalorimeterProperties;
+class HCALProperties;
 
 class CaloHitMaker
 {
@@ -46,13 +48,17 @@ class CaloHitMaker
 
   const CaloGeometryHelper * myCalorimeter;    
   const CalorimeterProperties * theCaloProperties;
+  const HCALProperties * theHCALProperties;
   double moliereRadius;
   double interactionLength;
   double spotEnergy;
 
   bool EMSHOWER;
   bool HADSHOWER;
   bool MIP;
+  bool numbering;
+  std::vector <double> hbLayers_;
+  std::vector <double> heLayers_;
 
  private:
   DetId::Detector base_;

FastSimulation/CaloHitMakers/interface/HcalHitMaker.h

@@ -5,6 +5,7 @@
 #include "FastSimulation/CaloHitMakers/interface/CaloHitMaker.h"
 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
+#include "SimDataFormats/CaloTest/interface/HcalTestNumbering.h"
 
 //#include "Math/GenVector/Transform3D.h"
 #include "FastSimulation/CaloGeometryTools/interface/Transform3DPJ.h"
@@ -39,6 +40,9 @@ class HcalHitMaker : public CaloHitMaker
   /// set the depth in X0 or Lambda0 units depending on showerType
   bool setDepth(double,bool inCm=false); 
 
+  int getLayerNumHE(double);
+  int getLayerNumHB(int,double,double);
+
  private:
     EcalHitMaker& myGrid;
 

FastSimulation/CaloHitMakers/src/CaloHitMaker.cc

@@ -8,7 +8,7 @@
 typedef ROOT::Math::Plane3D::Point Point;
 
 CaloHitMaker::CaloHitMaker(const CaloGeometryHelper * theCalo,DetId::Detector basedet,int subdetn,int cal,unsigned sht)
-  :myCalorimeter(theCalo),theCaloProperties(NULL),base_(basedet),subdetn_(subdetn),onCal_(cal),showerType_(sht)
+  :myCalorimeter(theCalo),theCaloProperties(NULL),theHCALProperties(NULL),base_(basedet),subdetn_(subdetn),onCal_(cal),showerType_(sht)
 {
   //  std::cout << " FamosCalorimeter " << basedet << " " << cal << std::endl;
   EMSHOWER=(sht==0);
@@ -20,6 +20,18 @@ CaloHitMaker::CaloHitMaker(const CaloGeometryHelper * theCalo,DetId::Detector ba
   if(base_==DetId::Ecal&&subdetn_==EcalPreshower&&onCal_)
     theCaloProperties = (PreshowerProperties*)myCalorimeter->layer1Properties(onCal_);
   if(base_==DetId::Hcal&&cal) theCaloProperties = myCalorimeter->hcalProperties(onCal_);
+  if(base_==DetId::Hcal&&cal) theHCALProperties = myCalorimeter->hcalProperties(onCal_);
+
+  if(theHCALProperties && (onCal_==1 || onCal_==2))
+    {
+      hbLayers_ = theHCALProperties->getMapLayersHB();
+      heLayers_ = theHCALProperties->getMapLayersHE();
+      numbering = theHCALProperties->getHcalTestNumbering();
+    } 
+  else 
+    {
+      numbering = false;
+    } 
 
   if(theCaloProperties)
     {

FastSimulation/CaloHitMakers/src/HcalHitMaker.cc

@@ -95,8 +95,25 @@ bool HcalHitMaker::addHit(const XYZPoint& point, unsigned layer)
 
   if(!thecellID.null()  && myDetId.depth()>0)
     {	
-      CaloHitID current_id(thecellID.rawId(),tof,0); //no track yet
-
+      uint32_t NewId;
+      if( numbering )
+      {
+        int lay = 1;
+        if(myDetId.subdet()==1) 
+          lay = getLayerNumHB(myDetId.iphi(),point.X(),point.Y())+1;
+        else if(myDetId.subdet()==2) lay = getLayerNumHE(point.Z())+1;
+        else if(myDetId.subdet()==3) lay = 1;
+        else lay = myDetId.depth();
+
+        int nZside = (myDetId.zside()==1) ? 1 : 0;
+        NewId = HcalTestNumbering::packHcalIndex(myDetId.subdet(),nZside,
+                    myDetId.depth(),myDetId.ietaAbs(),myDetId.iphi(),lay);
+      }
+      else 
+        NewId = thecellID.rawId();
+
+      CaloHitID current_id(NewId,tof,0); //no track yet
+
       //      std::cout << " FamosHcalHitMaker::addHit - the cell num " << cell
       //      		<< std::endl;
 
@@ -180,3 +197,26 @@ HcalHitMaker::setDepth(double depth,bool inCm)
 			   (Point)(origin+planeVec1));
   return true;
 }
+
+int HcalHitMaker::getLayerNumHE(double z)
+{
+  int lay = 0;
+  for (unsigned int i=0; i<heLayers_.size()-1; ++i) {
+    if(abs(z) <=heLayers_[0]) break; 
+    else if (abs(z) > heLayers_[i]) lay = i;
+  }
+  return lay;
+}
+
+int HcalHitMaker::getLayerNumHB(int phi, double x, double y)
+{
+  int lay = 0;
+  double alfa  = (2.*M_PI/174.) + (2.*M_PI/72.0)*(phi-1);  
+  double xPrim = x*cos(alfa) + y*sin(alfa);
+  for (unsigned int i=0; i<hbLayers_.size()-1; ++i) {
+    if(abs(xPrim) <=hbLayers_[0]) break;
+    else if (abs(xPrim) > hbLayers_[i]) lay = i;
+  }
+  return lay;
+}
+

FastSimulation/CalorimeterProperties/interface/HCALProperties.h

@@ -70,6 +70,11 @@ class HCALProperties : public CalorimeterProperties
 
   int eta2ieta(double eta) const; 
 
+  /// functions for rad. damage in HE
+  const std::vector<double> &getMapLayersHB() const {return hbLayers_;}
+  const std::vector<double> &getMapLayersHE() const {return heLayers_;}
+  bool getHcalTestNumbering() const {return HcalNumbering_;}
+
  private:
   double hOPi;
   double spotFrac;
@@ -85,7 +90,9 @@ class HCALProperties : public CalorimeterProperties
   double HCALinteractionLength_;
   std::vector <double> etatow_;// HCAL towers eta edges
   std::vector <double> hcalDepthLam_;// HCAL depth for each tower ieta
-
+  bool HcalNumbering_;
+  std::vector <double> hbLayers_;    // HB layers r-position
+  std::vector <double> heLayers_;    // HE layers z-position
 };
 
 #endif

FastSimulation/CalorimeterProperties/src/HCALProperties.cc

@@ -23,7 +23,9 @@ HCALProperties::HCALProperties(const edm::ParameterSet& fastDet) : CalorimeterPr
   HCALinteractionLength_= fastDetHCAL.getParameter<double>("HCALinteractionLength");
   etatow_=fastDetHCAL.getParameter<std::vector<double>>("HCALetatow");
   hcalDepthLam_=fastDetHCAL.getParameter<std::vector<double>>("HCALDepthLam");
-
+  hbLayers_=fastDetHCAL.getParameter<std::vector<double>>("layersHB");
+  heLayers_=fastDetHCAL.getParameter<std::vector<double>>("layersHE");
+  HcalNumbering_=fastDetHCAL.getUntrackedParameter<bool>("HcalTestNumbering",false);
   // in principle this splitting into 42 bins may change with future detectors, but let's add a protection to make sure that differences are not typos in the configuration file:
   if (etatow_.size() != 42) std::cout << " HCALProperties::eta2ieta - WARNING: here we expect 42 entries instead of " << etatow_.size() << "; is the change intentional?" << std::endl;
   // splitting of  28-th tower is taken into account (2.65-2.853-3.0)

FastSimulation/Calorimetry/interface/CalorimetryManager.h

@@ -13,6 +13,7 @@
 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
 #include "FastSimulation/CaloHitMakers/interface/PreshowerHitMaker.h"
+#include "SimDataFormats/CaloTest/interface/HcalTestNumbering.h"
 
 // For the uint32_t
 //#include <boost/cstdint.hpp>
@@ -118,6 +119,7 @@ class CalorimetryManager{
   //Digitizer
   bool EcalDigitizer_;
   bool HcalDigitizer_;
+  bool HcalTestNumbering_;
   std::vector<double> samplingHBHE_;
   std::vector<double> samplingHF_;
   std::vector<double> samplingHO_;

FastSimulation/Calorimetry/python/Calorimetry_cff.py

@@ -69,7 +69,11 @@
                   HCALinteractionLength= cms.double(15.05),
 
                   HCALetatow=cms.vdouble( 0.000, 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870,    0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830,    1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853, 3.000, 3.139, 3.314, 3.489,    3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191),
-                  HCALDepthLam=cms.vdouble( 8.930, 9.001, 9.132, 8.912, 8.104, 8.571, 8.852, 9.230, 9.732, 10.29,          10.95, 11.68, 12.49, 12.57, 12.63,  6.449, 5.806, 8.973, 8.934,  8.823,          8.727, 8.641, 8.565, 8.496, 8.436, 8.383, 8.346, 8.307, 8.298,  8.281,          9.442, 9.437, 9.432, 9.429, 9.432, 9.433, 9.430, 9.437, 9.442, 9.446, 9.435)
+                  HCALDepthLam=cms.vdouble( 8.930, 9.001, 9.132, 8.912, 8.104, 8.571, 8.852, 9.230, 9.732, 10.29,          10.95, 11.68, 12.49, 12.57, 12.63,  6.449, 5.806, 8.973, 8.934,  8.823,          8.727, 8.641, 8.565, 8.496, 8.436, 8.383, 8.346, 8.307, 8.298,  8.281,          9.442, 9.437, 9.432, 9.429, 9.432, 9.433, 9.430, 9.437, 9.442, 9.446, 9.435),
+
+                  HcalTestNumbering = cms.untracked.bool(False),
+                  layersHB = cms.vdouble(181.1, 188.7, 194.7, 200.7, 206.7, 212.7, 218.7, 224.7, 230.7, 236.7, 242.7, 249.3, 255.9, 262.5, 269.1, 275.7, 282.3, 288.8),
+                  layersHE = cms.vdouble(400.458, 408.718, 416.978, 425.248, 433.508, 441.768, 450.038, 458.298, 466.558, 474.828, 483.088, 491.348, 499.618, 507.878, 516.138, 524.398, 532.668, 540.928, 549.268)
                   ),
 
             BarrelCalorimeterProperties = cms.PSet(

FastSimulation/Calorimetry/src/CalorimetryManager.cc

@@ -550,7 +550,17 @@ void CalorimetryManager::reconstructHCAL(const FSimTrack& myTrack)
   if(emeas > 0.) {  
     DetId cell = myCalorimeter_->getClosestCell(trackPosition.Vect(),false,false);
 	double tof = (myCalorimeter_->getHcalGeometry()->getGeometry(cell)->getPosition().mag())/29.98;//speed of light
-	CaloHitID current_id(cell.rawId(),tof,myTrack.id());
+        uint32_t newId;
+        if( HcalTestNumbering_ )
+        {
+          HcalDetId myDetId(cell);
+          int nZside = (myDetId.zside()==1) ? 1 : 0;
+          newId = HcalTestNumbering::packHcalIndex(myDetId.subdet(),nZside,
+                        myDetId.depth(),myDetId.ietaAbs(),myDetId.iphi(),1);
+        } 
+        else newId = cell.rawId();
+
+	CaloHitID current_id(newId,tof,myTrack.id());
 	std::map<CaloHitID,float> hitMap;
 	hitMap[current_id] = emeas;
 	updateHCAL(hitMap,myTrack.id());
@@ -827,7 +837,17 @@ void CalorimetryManager::HDShowerSimulation(const FSimTrack& myTrack){//,
 	{
 	  DetId cell = myCalorimeter_->getClosestCell(trackPosition.Vect(),false,false);
 	  double tof = (myCalorimeter_->getHcalGeometry()->getGeometry(cell)->getPosition().mag())/29.98;//speed of light
-      CaloHitID current_id(cell.rawId(),tof,myTrack.id());
+          uint32_t newId;
+          if( HcalTestNumbering_ )
+          {
+            HcalDetId myDetId(cell);
+            int nZside = (myDetId.zside()==1) ? 1 : 0;
+            newId = HcalTestNumbering::packHcalIndex(myDetId.subdet(),nZside,
+                    myDetId.depth(),myDetId.ietaAbs(),myDetId.iphi(),myDetId.depth());
+          } 
+          else newId = cell.rawId();
+
+          CaloHitID current_id(newId,tof,myTrack.id());
 	  std::map<CaloHitID,float> hitMap;
 	  hitMap[current_id] = emeas;
 	  updateHCAL(hitMap,myTrack.id());
@@ -1195,6 +1215,10 @@ void CalorimetryManager::readParameters(const edm::ParameterSet& fastCalo) {
   timeShiftHF_  = HCALparameters.getParameter< std::vector<double> >("timeShiftHF");
   timeShiftHO_  = HCALparameters.getParameter< std::vector<double> >("timeShiftHO");
 
+  // VA
+  edm::ParameterSet s0 = fastCalo.getParameter<edm::ParameterSet>("CalorimeterProperties");
+  edm::ParameterSet s1 = s0.getParameter<edm::ParameterSet>("HadronicCalorimeterProperties");
+  HcalTestNumbering_   = s1.getUntrackedParameter<bool>("HcalTestNumbering",false);
 }
 
 void CalorimetryManager::respCorr(double p) {
@@ -1293,23 +1317,30 @@ void CalorimetryManager::updateHCAL(const std::map<CaloHitID,float>& hitMap, int
 	float time = mapitr->first.timeSlice();
 	//put energy into uncalibrated state for digitizer && correct timing
 	if(HcalDigitizer_){
-	  HcalDetId hdetid = HcalDetId(mapitr->first.unitID());
+          HcalDetId hdetid = HcalDetId(mapitr->first.unitID());
+          if( HcalTestNumbering_ ) 
+          {
+            int det, z, depth, eta, phi, lay;
+            HcalTestNumbering::unpackHcalIndex(hdetid,det,z,depth,eta,phi,lay);
+            HcalDetId newCell((HcalSubdetector)det,eta,phi,depth);
+            hdetid = newCell;
+          }
 	  if (hdetid.subdetId()== HcalBarrel){
-        energy /= samplingHBHE_[hdetid.ietaAbs()-1]; //re-convert to GeV
-		time = timeShiftHB_[hdetid.ietaAbs()-ietaShiftHB_];
-      }
+            energy /= samplingHBHE_[hdetid.ietaAbs()-1]; //re-convert to GeV
+	    time = timeShiftHB_[hdetid.ietaAbs()-ietaShiftHB_];
+          }
 	  else if (hdetid.subdetId()== HcalEndcap){
 	    energy /= samplingHBHE_[hdetid.ietaAbs()-1]; //re-convert to GeV
-		time = timeShiftHE_[hdetid.ietaAbs()-ietaShiftHE_];
+            time = timeShiftHE_[hdetid.ietaAbs()-ietaShiftHE_];
 	  }
 	  else if (hdetid.subdetId()== HcalForward){
 	    if(hdetid.depth()== 1) energy *= samplingHF_[0];
 	    if(hdetid.depth()== 2) energy *= samplingHF_[1];
-		time = timeShiftHF_[hdetid.ietaAbs()-ietaShiftHF_];
+	    time = timeShiftHF_[hdetid.ietaAbs()-ietaShiftHF_];
 	  }
 	  else if (hdetid.subdetId()== HcalOuter){
-        energy /= samplingHO_[hdetid.ietaAbs()-1];
-		time = timeShiftHO_[hdetid.ietaAbs()-ietaShiftHO_];
+            energy /= samplingHO_[hdetid.ietaAbs()-1];
+            time = timeShiftHO_[hdetid.ietaAbs()-ietaShiftHO_];
 	  }
 	}