Skip to content

Commit

Permalink
Split the SD for HGCal for silicon and scintillators
Browse files Browse the repository at this point in the history
  • Loading branch information
Sunanda committed Jul 10, 2018
1 parent 4d160db commit ec638ad
Show file tree
Hide file tree
Showing 10 changed files with 488 additions and 55 deletions.
4 changes: 4 additions & 0 deletions Geometry/HGCalCommonData/interface/HGCalDDDConstants.h
View file
Expand Up @@ -38,6 +38,8 @@ class HGCalDDDConstants {
int lay, bool reco) const;
double cellSizeHex(int type) const;
double cellThickness(int layer, int waferU, int waferV) const;
double distFromEdgeHex(double x, double y, double z) const;
double distFromEdgeTrap(double x, double y, double z) const;
void etaPhiFromPosition(const double x, const double y,
const double z, const int layer,
int& ieta, int& iphi, int& type,
Expand Down Expand Up @@ -129,6 +131,8 @@ class HGCalDDDConstants {
std::pair<int,float> getIndex(int lay, bool reco) const;
bool isValidCell(int layindex, int wafer, int cell) const;
bool waferInLayerTest(int wafer, int lay, bool full) const;
bool waferGlobal2Local(double& xx, double& yy, int& wafer,
int& celltyp) const;

const double k_horizontalShift = 1.0;
const float dPhiMin = 0.02;
Expand Down
121 changes: 113 additions & 8 deletions Geometry/HGCalCommonData/src/HGCalDDDConstants.cc
View file
Expand Up @@ -15,6 +15,7 @@
#include <numeric>

//#define EDM_ML_DEBUG

static const int maxType = 2;
static const int minType = 0;

Expand Down Expand Up @@ -242,10 +243,89 @@ double HGCalDDDConstants::cellSizeHex(int type) const {
return cell;
}

double HGCalDDDConstants::distFromEdgeHex(double x, double y, double z) const {

if (z < 0) x = -x;
double dist(0);
//Input x, y in Geant4 unit and transformed to CMSSW standard
double xx = HGCalParameters::k_ScaleFromDDD*x;
double yy = HGCalParameters::k_ScaleFromDDD*y;
int sizew = (int)(hgpar_->waferPosX_.size());
int wafer = sizew;
for (int k=0; k<sizew; ++k) {
double dx = std::abs(xx-hgpar_->waferPosX_[k]);
double dy = std::abs(yy-hgpar_->waferPosY_[k]);
if (dx <= rmax_ && dy <= hexside_) {
if ((dy <= 0.5*hexside_) || (dx*tan30deg_ <= (hexside_-dy))) {
wafer = k;
xx -= hgpar_->waferPosX_[k];
yy -= hgpar_->waferPosY_[k];
break;
}
}
}
if (wafer < sizew) {
if (std::abs(yy) < 0.5*hexside_) {
dist = rmax_ - std::abs(xx);
} else {
dist = 0.5*((rmax_-std::abs(xx))-sqrt3_*(std::abs(yy)-0.5*hexside_));
}
} else {
dist = 0;
}
dist *= HGCalParameters::k_ScaleToDDD;
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "DistFromEdgeHex: Local " << xx << ":"
<< yy << " wafer " << wafer << " flag "
<< (wafer < sizew) << " Distance " << rmax_
<< ":" << (rmax_-std::abs(xx)) << ":"
<< (std::abs(yy)-0.5*hexside_) << ":"
<< 0.5*hexside_ << ":" << dist;
#endif
return dist;
}

double HGCalDDDConstants::distFromEdgeTrap(double x, double y, double z) const {

int lay = getLayer(z,false);
double xx = (z < 0) ? -x : x;
int indx = layerIndex(lay,false);
double zz = HGCalParameters::k_ScaleToDDD*hgpar_->zLayerHex_[indx];
double r = std::sqrt(x*x+y*y+zz*zz);
double theta = (r == 0. ? 0. : std::acos(std::max(std::min(zz/r,1.0),-1.0)));
double stheta= std::sin(theta);
double phi = (r*stheta == 0. ? 0. : std::atan2(y,xx));
if (phi < 0) phi += (2.0*M_PI);
double eta = (std::abs(stheta) == 1.0 ? 0. : -std::log(std::abs(std::tan(0.5*theta))) );
double cell = hgpar_->dPhiEtaBH_[indx];
int ieta = 1 + (int)((std::abs(eta)-hgpar_->etaMinBH_)/cell);
int iphi = 1 + (int)(phi/cell);
double rr = std::sqrt(x*x+y*y);
double dphi = std::max(0.0,(0.5*cell-std::abs(phi-(iphi-0.5)*cell)));
double dist = ((eta > hgpar_->etaMinBH_+(ieta-0.5)*cell) ?
(rr - zz/std::sinh(hgpar_->etaMinBH_+ieta*cell)) :
(zz/std::sinh(hgpar_->etaMinBH_+(ieta-1)*cell) - rr));
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HGCalGeom") << "DistFromEdgeTrap: Global " << x << ":"
<< y << ":" << z << " Layer " << lay
<< " Index " << indx << " xx|zz " << xx << ":"
<< zz << " Eta " << eta << ":"<< ieta << ":"
<< (hgpar_->etaMinBH_+(ieta-0.5)*cell)
<< " Phi " << phi << ":" << iphi << ":"
<< (iphi-0.5)*cell << " cell " << cell << " R "
<< rr << " Dphi " << dphi << " Dist " << dist
<< ":" << rr*dphi << ":"
<< rr-zz/std::sinh(hgpar_->etaMinBH_+ieta*cell)
<< ":" << zz/std::sinh(hgpar_->etaMinBH_+(ieta-1)*cell)-rr;
#endif
return ((dist > rr*dphi) ? rr*dphi : dist);
}

int HGCalDDDConstants::getLayer(double z, bool reco) const {

unsigned int k = 0;
double zz = std::abs(z);
double zz = (reco ? std::abs(z) :
HGCalParameters::k_ScaleFromDDD*std::abs(z));
const auto& zLayerHex = hgpar_->zLayerHex_;
std::find_if(zLayerHex.begin()+1,zLayerHex.end(),[&k,&zz,&zLayerHex](double zLayer){ ++k; return zz < 0.5*(zLayerHex[k-1]+zLayerHex[k]);});
int lay = k;
Expand Down Expand Up @@ -750,20 +830,21 @@ int HGCalDDDConstants::waferFromCopy(int copy) const {
void HGCalDDDConstants::waferFromPosition(const double x, const double y,
int& wafer, int& icell,
int& celltyp) const {
//Input x, y in Geant4 unit and transformed to conform +z convention
double xx(HGCalParameters::k_ScaleFromDDD*x), yy(HGCalParameters::k_ScaleFromDDD*y);
//Input x, y in Geant4 unit and transformed to CMSSW standard
double xx = HGCalParameters::k_ScaleFromDDD*x;
double yy = HGCalParameters::k_ScaleFromDDD*y;
int size_ = (int)(hgpar_->waferCopy_.size());
wafer = size_;
for (int k=0; k<size_; ++k) {
double dx = std::abs(xx-hgpar_->waferPosX_[k]);
double dy = std::abs(yy-hgpar_->waferPosY_[k]);
if (dx <= rmax_ && dy <= hexside_) {
if ((dy <= 0.5*hexside_) || (dx*tan30deg_ <= (hexside_-dy))) {
wafer = k;
celltyp = hgpar_->waferTypeT_[k];
xx -= hgpar_->waferPosX_[k];
yy -= hgpar_->waferPosY_[k];
break;
wafer = k;
celltyp = hgpar_->waferTypeT_[k];
xx -= hgpar_->waferPosX_[k];
yy -= hgpar_->waferPosY_[k];
break;
}
}
}
Expand Down Expand Up @@ -1067,6 +1148,30 @@ bool HGCalDDDConstants::waferInLayerTest(int wafer, int lay, bool full) const {
return in;
}

bool HGCalDDDConstants::waferGlobal2Local(double& xx, double& yy, int& wafer,
int& celltyp) const {
std::cout << xx <<":" << yy << ":" << wafer << ":" << celltyp << ":" << hgpar_ << "\n";
int sizew = (int)(hgpar_->waferPosX_.size());
wafer = sizew;
for (int k=0; k<sizew; ++k) {
double dx = std::abs(xx-hgpar_->waferPosX_[k]);
double dy = std::abs(yy-hgpar_->waferPosY_[k]);
if (dx <= rmax_ && dy <= hexside_) {
if ((dy <= 0.5*hexside_) || (dx*tan30deg_ <= (hexside_-dy))) {
wafer = k;
celltyp = hgpar_->waferTypeT_[k];
xx -= hgpar_->waferPosX_[k];
yy -= hgpar_->waferPosY_[k];
break;
}
}
}
edm::LogVerbatim("HGCalGeom") << "WaferGlobal2Local: Global " << xx << ":"
<< yy << " wafer " << wafer << " Type "
<< celltyp << " flag " << (wafer < sizew);
return (wafer < sizew);
}

#include "FWCore/Utilities/interface/typelookup.h"

TYPELOOKUP_DATA_REG(HGCalDDDConstants);
22 changes: 14 additions & 8 deletions Geometry/HGCalCommonData/test/HGCalNumberingTester.cc
View file
Expand Up @@ -111,6 +111,7 @@ void HGCalNumberingTester::analyze( const edm::Event& iEvent, const edm::EventSe
std::string flg;
int subsec(0);
int loff = hgdc.firstLayer();
double scl = (reco_ ? 10.0 : 1.0);
for (unsigned int k=0; k<positionX_.size(); ++k) {
float localx(positionX_[k]), localy(positionY_[k]);
for (unsigned int i=0; i<hgdc.layers(reco_); ++i) {
Expand Down Expand Up @@ -145,8 +146,9 @@ void HGCalNumberingTester::analyze( const edm::Event& iEvent, const edm::EventSe
<< ", " << i+loff << "), assignCell o/p (" << kxy[0] << ":"
<< kxy[1] << ":" << kxy[2] << ") locateCell o/p ("
<< xy.first << ", " << xy.second << ")," << " final ("
<< lxy[0] << ":" << lxy[1] << ":" << lxy[2] << ")" << flg
<< std::endl;
<< lxy[0] << ":" << lxy[1] << ":" << lxy[2] << ") Dist "
<< hgdc.distFromEdgeTrap(scl*localx,scl*localy,scl*zpos)
<< " " << flg << std::endl;
kxy = hgdc.assignCellTrap(-localx,-localy,zpos,i+loff,reco_);
xy = hgdc.locateCellTrap(i+loff,kxy[0],kxy[1],reco_);
lxy = hgdc.assignCellTrap(xy.first,xy.second,zpos,i+loff,reco_);
Expand All @@ -155,21 +157,24 @@ void HGCalNumberingTester::analyze( const edm::Event& iEvent, const edm::EventSe
<< ", " << i+loff << "), assignCell o/p (" << kxy[0] << ":"
<< kxy[1] << ":" << kxy[2] << ") locateCell o/p ("
<< xy.first << ", " << xy.second << ")," << " final ("
<< lxy[0] << ":" << lxy[1] << ":" << lxy[2] << ")" << flg
<< std::endl;
<< lxy[0] << ":" << lxy[1] << ":" << lxy[2] << ") Dist "
<< hgdc.distFromEdgeTrap(scl*localx,scl*localy,scl*zpos)
<< " " << flg << std::endl;
} else {
std::array<int,5> kxy, lxy;
kxy = hgdc.assignCellHex(localx,localy,i+loff,reco_);
xy = hgdc.locateCell(i+loff,kxy[0],kxy[1],kxy[3],kxy[4],reco_,true);
lxy = hgdc.assignCellHex(xy.first,xy.second,i+loff,reco_);
flg = (kxy == lxy) ? " " : " ***** Error *****";
double zpos = hgdc.waferZ(i+loff,reco_);
std::cout << "Input: (" << localx << "," << localy << ", " << i+loff
<< "), assignCell o/p (" << kxy[0] << ":" << kxy[1]
<< ":" << kxy[2] << ":" << kxy[3] << ":" << kxy[4]
<< ") locateCell o/p (" << xy.first << ", " << xy.second
<< ")," << " final (" << lxy[0] << ":" << lxy[1] << ":"
<< lxy[2] << ":" << lxy[3] << ":" << lxy[4] << ")" << flg
<< std::endl;
<< lxy[2] << ":" << lxy[3] << ":" << lxy[4] << ") Dist "
<< hgdc.distFromEdgeHex(scl*localx,scl*localy,scl*zpos)
<< " " << flg << std::endl;
kxy = hgdc.assignCellHex(-localx,-localy,i+loff,reco_);
xy = hgdc.locateCell(i+loff,kxy[0],kxy[1],kxy[3],kxy[4],reco_,true);
lxy = hgdc.assignCellHex(xy.first,xy.second,i+loff,reco_);
Expand All @@ -179,8 +184,9 @@ void HGCalNumberingTester::analyze( const edm::Event& iEvent, const edm::EventSe
<< ":" << kxy[2] << ":" << kxy[3] << ":" << kxy[4]
<< ") locateCell o/p (" << xy.first << ", " << xy.second
<< ")," << " final (" << lxy[0] << ":" << lxy[1] << ":"
<< lxy[2] << ":" << lxy[3] << ":" << lxy[4] << ")" << flg
<< std::endl;
<< lxy[2] << ":" << lxy[3] << ":" << lxy[4] << ") Dist "
<< hgdc.distFromEdgeHex(scl*localx,scl*localy,scl*zpos)
<< " " << flg << std::endl;
}
if (k == 0 && i==0 && detType_ == 1) {
std::vector<int> ncells = hgdc.numberCells(i+1,reco_);
Expand Down
2 changes: 1 addition & 1 deletion Geometry/HGCalSimData/data/hgcsensv9.xml
View file
Expand Up @@ -15,7 +15,7 @@
</SpecPar>
<SpecPar name="hgchesci">
<PartSelector path="//HGCalHEScintillatorSensitive.*"/>
<Parameter name="SensitiveDetector" value="HGCalSensitiveDetector" eval="false"/>
<Parameter name="SensitiveDetector" value="HGCScintillatorSensitiveDetector" eval="false"/>
<Parameter name="ReadOutName" value="HGCHitsHEback" eval="false"/>
</SpecPar>
</SpecParSection>
Expand Down
54 changes: 54 additions & 0 deletions SimG4CMS/Calo/interface/HGCScintSD.h
View file
@@ -0,0 +1,54 @@
#ifndef SimG4CMS_HGCScintSD_h
#define SimG4CMS_HGCScintSD_h
///////////////////////////////////////////////////////////////////////////////
// File: HGCScintSD.h
// Description: Stores hits of the High Granularity Calorimeter (HGC) in the
// appropriate container (post TDR version)
///////////////////////////////////////////////////////////////////////////////

#include "SimG4CMS/Calo/interface/CaloSD.h"
#include "SimG4Core/Notification/interface/BeginOfJob.h"
#include "SimG4CMS/Calo/interface/HGCalNumberingScheme.h"

#include <string>

class DDCompactView;
class HGCalDDDConstants;
class G4LogicalVolume;
class G4Step;

class HGCScintSD : public CaloSD, public Observer<const BeginOfJob *> {

public:

HGCScintSD(const std::string& , const DDCompactView &,
const SensitiveDetectorCatalog &, edm::ParameterSet const &,
const SimTrackManager*);
~HGCScintSD() override;

uint32_t setDetUnitId(const G4Step* step) override;

protected:

double getEnergyDeposit(const G4Step*) override;
void update(const BeginOfJob *) override;
void initRun() override;
bool filterHit(CaloG4Hit*, double) override;

private:

uint32_t setDetUnitId(int, int, int, int, G4ThreeVector &);
bool isItinFidVolume (const G4ThreeVector&);

const HGCalDDDConstants* hgcons_;
HGCalNumberingScheme* numberingScheme_;
DetId::Detector mydet_;
std::string nameX_;
HGCalGeometryMode::GeometryMode geom_mode_;
double eminHit_, slopeMin_, distanceFromEdge_;
int levelT1_, levelT2_;
bool storeAllG4Hits_, fiducialCut_, useBirk_;
double birk1_, birk2_, birk3_, weight_;
};

#endif // HGCScintSD_h
12 changes: 7 additions & 5 deletions SimG4CMS/Calo/interface/HGCalSD.h
View file
Expand Up @@ -14,6 +14,7 @@
#include <string>

class DDCompactView;
class HGCalDDDConstants;
class G4LogicalVolume;
class G4Step;

Expand All @@ -37,19 +38,20 @@ class HGCalSD : public CaloSD, public Observer<const BeginOfJob *> {
private:

uint32_t setDetUnitId(int, int, int, int, G4ThreeVector &);
bool isItinFidVolume (const G4ThreeVector&) {return true;}
bool isItinFidVolume (const G4ThreeVector&);

const HGCalDDDConstants* hgcons_;
HGCalNumberingScheme* numberingScheme_;
HGCMouseBite* mouseBite_;
DetId::Detector mydet_;
std::string nameX_;
HGCalGeometryMode::GeometryMode geom_mode_;
double eminHit_, slopeMin_, mouseBiteCut_;
double eminHit_, slopeMin_;
double distanceFromEdge_, mouseBiteCut_, weight_;
int levelT1_, levelT2_;
bool storeAllG4Hits_, rejectMB_, waferRot_;
bool useBirk_, isScint_;
bool storeAllG4Hits_;
bool fiducialCut_, rejectMB_, waferRot_;
const double tan30deg_;
double birk1_, birk2_, birk3_, weight_;
std::vector<double> angles_;
};

Expand Down
3 changes: 3 additions & 0 deletions SimG4CMS/Calo/plugins/module.cc
View file
Expand Up @@ -2,6 +2,7 @@
#include "SimG4CMS/Calo/interface/HCalSD.h"
#include "SimG4CMS/Calo/interface/HGCSD.h"
#include "SimG4CMS/Calo/interface/HGCalSD.h"
#include "SimG4CMS/Calo/interface/HGCScintSD.h"
#include "SimG4CMS/Calo/interface/CaloTrkProcessing.h"
#include "SimG4CMS/Calo/interface/HcalTestAnalysis.h"
#include "SimG4Core/SensitiveDetector/interface/SensitiveDetectorPluginFactory.h"
Expand All @@ -17,6 +18,8 @@ typedef HGCSD HGCSensitiveDetector;
DEFINE_SENSITIVEDETECTOR(HGCSensitiveDetector);
typedef HGCalSD HGCalSensitiveDetector;
DEFINE_SENSITIVEDETECTOR(HGCalSensitiveDetector);
typedef HGCScintSD HGCScintillatorSensitiveDetector;
DEFINE_SENSITIVEDETECTOR(HGCScintillatorSensitiveDetector);

DEFINE_SENSITIVEDETECTOR(CaloTrkProcessing);

Expand Down

0 comments on commit ec638ad

Please sign in to comment.