Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Phase2-hgx124 Split the SD for HGCal for silicon and scintillators #23772

Merged
merged 2 commits into from Jul 12, 2018
Merged

Phase2-hgx124 Split the SD for HGCal for silicon and scintillators #23772

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ec638ad
Split the SD for HGCal for silicon and scintillators
Jul 10, 2018
512425c
Take care of Kevin's suggestions
Jul 10, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 2 additions & 0 deletions Geometry/HGCalCommonData/interface/HGCalDDDConstants.h
View file
Open in desktop
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
96 changes: 88 additions & 8 deletions Geometry/HGCalCommonData/src/HGCalDDDConstants.cc
View file
Open in desktop
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,88 @@ double HGCalDDDConstants::cellSizeHex(int type) const {
return cell;
}

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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bsunanda - please, add a comment what is calculated here. Is the point(x,y) inside of the polygon?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@bsunanda it takes longer for me to check the code without comment on what it suppose to do.


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_ &&
((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 std::min(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 +829,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
22 changes: 14 additions & 8 deletions Geometry/HGCalCommonData/test/HGCalNumberingTester.cc
View file
Open in desktop
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
Open in desktop
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
28 changes: 14 additions & 14 deletions SimG4CMS/Calo/interface/HGCSD.h
View file
Open in desktop
Expand Up @@ -24,9 +24,10 @@ class HGCSD : public CaloSD, public Observer<const BeginOfJob *> {

public:

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

uint32_t setDetUnitId(const G4Step* step) override;

Expand All @@ -43,18 +44,17 @@ class HGCSD : public CaloSD, public Observer<const BeginOfJob *> {
int, int, G4ThreeVector &);
bool isItinFidVolume (const G4ThreeVector&) {return true;}

std::string nameX_;

HGCalGeometryMode::GeometryMode geom_mode_;
HGCNumberingScheme* numberingScheme_;
HGCMouseBite* mouseBite_;
double eminHit_;
ForwardSubdetector myFwdSubdet_;
double slopeMin_;
int levelT_;
bool storeAllG4Hits_, rejectMB_, waferRot_;
double mouseBiteCut_;
std::vector<double> angles_;
std::string nameX_;
HGCalGeometryMode::GeometryMode geom_mode_;
std::unique_ptr<HGCNumberingScheme> numberingScheme_;
std::unique_ptr<HGCMouseBite> mouseBite_;
double eminHit_;
ForwardSubdetector myFwdSubdet_;
double slopeMin_;
int levelT_;
bool storeAllG4Hits_, rejectMB_, waferRot_;
double mouseBiteCut_;
std::vector<double> angles_;
};

#endif // HGCSD_h
55 changes: 55 additions & 0 deletions SimG4CMS/Calo/interface/HGCScintSD.h
View file
Open in desktop
@@ -0,0 +1,55 @@
#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 = default;

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_;
std::unique_ptr<HGCalNumberingScheme> numberingScheme_;
DetId::Detector mydet_;
std::string nameX_;
HGCalGeometryMode::GeometryMode geom_mode_;
double eminHit_, slopeMin_, distanceFromEdge_;
int levelT1_, levelT2_;
bool storeAllG4Hits_, fiducialCut_;
bool useBirk_;
double birk1_, birk2_, birk3_, weight_;
};

#endif // HGCScintSD_h
35 changes: 19 additions & 16 deletions SimG4CMS/Calo/interface/HGCalSD.h
View file
Open in desktop
Expand Up @@ -14,16 +14,18 @@
#include <string>

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

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

public:

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

uint32_t setDetUnitId(const G4Step* step) override;

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

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

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

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

#endif // HGCalSD_h
3 changes: 3 additions & 0 deletions SimG4CMS/Calo/plugins/module.cc
View file
Open in desktop
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