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

Fix RPCGeometryBuilder to set correct chamber surface #14848

Merged
merged 8 commits into from Jun 17, 2016
Merged

Fix RPCGeometryBuilder to set correct chamber surface #14848

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
b8a8139
code cleanup minor speedup using set
jhgoh Jun 10, 2016
06d1d0b
constify variables
jhgoh Jun 10, 2016
19f6ecf
constify and cleanup
jhgoh Jun 10, 2016
0ef3692
compute RPCChamber bounds by envelope of all rolls
jhgoh Jun 10, 2016
bbd4719
should be half height
jhgoh Jun 10, 2016
58c632c
should be half height
jhgoh Jun 10, 2016
6725764
typo
jhgoh Jun 11, 2016
ec40e57
thickness from Z difference
jhgoh Jun 11, 2016
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
188 changes: 107 additions & 81 deletions Geometry/RPCGeometryBuilder/src/RPCGeometryBuilderFromCondDB.cc
View file
Open in desktop
Expand Up @@ -16,6 +16,7 @@

#include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
#include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
#include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"

#include "DataFormats/GeometryVector/interface/Basic3DVector.h"

Expand All @@ -28,57 +29,45 @@ RPCGeometryBuilderFromCondDB::RPCGeometryBuilderFromCondDB(bool comp11) :
theComp11Flag(comp11)
{ }

RPCGeometryBuilderFromCondDB::~RPCGeometryBuilderFromCondDB()
RPCGeometryBuilderFromCondDB::~RPCGeometryBuilderFromCondDB()
{ }

RPCGeometry* RPCGeometryBuilderFromCondDB::build(const RecoIdealGeometry& rgeo)
{
const std::vector<DetId>& detids(rgeo.detIds());
RPCGeometry* geometry = new RPCGeometry();

std::string name;
std::vector<double>::const_iterator tranStart;
std::vector<double>::const_iterator shapeStart;
std::vector<double>::const_iterator rotStart;
std::vector<std::string>::const_iterator strStart;

for (unsigned int id=0; id<detids.size(); id++){


for (unsigned int id=0; id<detids.size(); ++id){
RPCDetId rpcid(detids[id]);
RPCDetId chid(rpcid.region(),rpcid.ring(),rpcid.station(),
rpcid.sector(),rpcid.layer(),rpcid.subsector(),0);

tranStart=rgeo.tranStart(id);
shapeStart=rgeo.shapeStart(id);
rotStart=rgeo.rotStart(id);
strStart=rgeo.strStart(id);
name=*(strStart);
rpcid.sector(),rpcid.layer(),rpcid.subsector(),0);

const auto tranStart = rgeo.tranStart(id);
const auto shapeStart = rgeo.shapeStart(id);
const auto rotStart = rgeo.rotStart(id);
const std::string& name = *rgeo.strStart(id);

Surface::PositionType pos(*(tranStart)/cm,*(tranStart+1)/cm, *(tranStart+2)/cm);
// CLHEP way
Surface::RotationType rot(*(rotStart+0),*(rotStart+1),*(rotStart+2),
*(rotStart+3),*(rotStart+4),*(rotStart+5),
*(rotStart+6),*(rotStart+7),*(rotStart+8));

std::vector<float> pars;
*(rotStart+3),*(rotStart+4),*(rotStart+5),
*(rotStart+6),*(rotStart+7),*(rotStart+8));

RPCRollSpecs* rollspecs= 0;
Bounds* bounds = 0;

// if (dpar.size()==4){
if ( rgeo.shapeEnd(id) - shapeStart == 4 ) {
float width = *(shapeStart+0)/cm;
float length = *(shapeStart+1)/cm;
float thickness = *(shapeStart+2)/cm;
float nstrip = *(shapeStart+3);
//RectangularPlaneBounds*
bounds =
new RectangularPlaneBounds(width,length,thickness);
pars.push_back(width);
pars.push_back(length);
pars.push_back(nstrip);
const float width = *(shapeStart+0)/cm;
const float length = *(shapeStart+1)/cm;
const float thickness = *(shapeStart+2)/cm;
const float nstrip = *(shapeStart+3);
//RectangularPlaneBounds*
bounds = new RectangularPlaneBounds(width,length,thickness);
const std::vector<float> pars = {width, length, nstrip};

if (!theComp11Flag) {
//Correction of the orientation to get the REAL geometry.
//Correction of the orientation to get the REAL geometry.
//Change of axes for the +z part only.
//Including the 0 whell
if ( *(tranStart+2) > -1500.) { //tran[2] >-1500. ){
Expand All @@ -89,76 +78,113 @@ RPCGeometry* RPCGeometryBuilderFromCondDB::build(const RecoIdealGeometry& rgeo)
}
}


rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCBarrel,name,pars);


}else{
float be = *(shapeStart+0)/cm;
float te = *(shapeStart+1)/cm;
float ap = *(shapeStart+2)/cm;
float ti = *(shapeStart+3)/cm;
float nstrip = *(shapeStart+4);
// TrapezoidalPlaneBounds*
bounds =
new TrapezoidalPlaneBounds(be,te,ap,ti);
pars.push_back(be); //b/2;
pars.push_back(te); //B/2;
pars.push_back(ap); //h/2;
pars.push_back(nstrip);


}
else {
const float be = *(shapeStart+0)/cm;
const float te = *(shapeStart+1)/cm;
const float ap = *(shapeStart+2)/cm;
const float ti = *(shapeStart+3)/cm;
const float nstrip = *(shapeStart+4);
// TrapezoidalPlaneBounds*
bounds = new TrapezoidalPlaneBounds(be,te,ap,ti);
const std::vector<float> pars = {be /*b/2*/, te /*B/2*/, ap /*h/2*/, nstrip};

rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCEndcap,name,pars);

//Change of axes for the forward
Basic3DVector<float> newX(1.,0.,0.);
Basic3DVector<float> newY(0.,0.,1.);
// if (tran[2] > 0. )
newY *= -1;
Basic3DVector<float> newZ(0.,1.,0.);
rot.rotateAxes (newX, newY,newZ);
rot.rotateAxes (newX, newY,newZ);
}

BoundPlane* bp = new BoundPlane(pos,rot,bounds);
ReferenceCountingPointer<BoundPlane> surf(bp);
RPCRoll* r=new RPCRoll(rpcid,surf,rollspecs);
RPCRoll* r = new RPCRoll(rpcid,surf,rollspecs);
geometry->add(r);

std::list<RPCRoll *> rls;
if (chids.find(chid)!=chids.end()){
rls = chids[chid];
}
rls.push_back(r);
chids[chid]=rls;

auto rls = chids.find(chid);
if ( rls == chids.end() ) rls = chids.insert(std::make_pair(chid, std::set<RPCRoll*>())).first;
rls->second.insert(r);
}
// Create the RPCChambers and store them on the Geometry
for( std::map<RPCDetId, std::list<RPCRoll *> >::iterator ich=chids.begin();
ich != chids.end(); ich++){
RPCDetId chid = ich->first;
std::list<RPCRoll * > rls = ich->second;

// compute the overall boundplane. At the moment we use just the last
// surface

// Create the RPCChambers and store them on the Geometry
for ( auto ich=chids.begin(); ich != chids.end(); ++ich ) {
const RPCDetId& chid = ich->first;
const auto& rls = ich->second;

// compute the overall boundplane.
BoundPlane* bp=0;
for(std::list<RPCRoll *>::iterator rl=rls.begin();
rl!=rls.end(); rl++){
const BoundPlane& bps = (*rl)->surface();
bp = const_cast<BoundPlane *>(&bps);
if ( !rls.empty() ) {
// First set the baseline plane to calculate relative poisions
const auto& refSurf = (*rls.begin())->surface();
if ( chid.region() == 0 ) {
float corners[6] = {0,0,0,0,0,0};
for ( auto rl : rls ) {
const double h2 = rl->surface().bounds().length()/2;
const double w2 = rl->surface().bounds().width()/2;
const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-w2,-h2,0)));
const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+w2,+h2,0)));
corners[0] = std::min(corners[0], x1y1AtRef.x());
corners[1] = std::min(corners[1], x1y1AtRef.y());
corners[2] = std::max(corners[2], x2y2AtRef.x());
corners[3] = std::max(corners[3], x2y2AtRef.y());

corners[4] = std::min(corners[4], x1y1AtRef.z());
corners[5] = std::max(corners[5], x1y1AtRef.z());
}
const LocalPoint lpOfCentre((corners[0]+corners[2])/2, (corners[1]+corners[3])/2, 0);
const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
auto bounds = new RectangularPlaneBounds((corners[2]-corners[0])/2, (corners[3]-corners[1])/2, (corners[5]-corners[4])+0.5);
bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
}
else {
float cornersLo[3] = {0,}, cornersHi[3] = {0,};
float cornersZ[2] = {0,};
Copy link
Contributor

Choose a reason for hiding this comment

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

are the higher (above 0) indices intentionally uninitialized?
You will likely get values from the previous calls in cornersLo[1, 2] etc

Copy link
Contributor Author

Choose a reason for hiding this comment

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

they are initialized to be zero by default.

Copy link
Contributor

Choose a reason for hiding this comment

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

where are they initialized to zero? float cornersZ[2] = {0,} only sets the cornersZ[0]

for ( auto rl : rls ) {
const double h2 = rl->surface().bounds().length()/2;
const double w2 = rl->surface().bounds().width()/2;
const auto& topo = dynamic_cast<const TrapezoidalStripTopology&>(rl->specificTopology());
const double r = topo.radius();
const double wAtLo = w2/r*(r-h2); // tan(theta/2) = (w/2)/r = x/(r-h/2)
const double wAtHi = w2/r*(r+h2);

const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtLo, -h2, 0)));
const auto x2y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtLo, -h2, 0)));
const auto x1y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtHi, +h2, 0)));
const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtHi, +h2, 0)));

cornersLo[0] = std::min(cornersLo[0], x1y1AtRef.x());
cornersLo[1] = std::max(cornersLo[1], x2y1AtRef.x());
cornersLo[2] = std::min(cornersLo[2], x1y1AtRef.y());

cornersHi[0] = std::min(cornersHi[0], x1y2AtRef.x());
cornersHi[1] = std::max(cornersHi[1], x2y2AtRef.x());
cornersHi[2] = std::max(cornersHi[2], x1y2AtRef.y());

cornersZ[0] = std::min(cornersZ[0], x1y1AtRef.z());
cornersZ[1] = std::max(cornersZ[1], x1y1AtRef.z());
}
const LocalPoint lpOfCentre((cornersHi[0]+cornersHi[1])/2, (cornersLo[2]+cornersHi[2])/2, 0);
const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
auto bounds = new TrapezoidalPlaneBounds((cornersLo[1]-cornersLo[0])/2, (cornersHi[1]-cornersHi[0])/2, (cornersHi[2]-cornersLo[2])/2, (cornersZ[1]-cornersZ[0])+0.5);
bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
}
}

ReferenceCountingPointer<BoundPlane> surf(bp);
// Create the chamber
RPCChamber* ch = new RPCChamber (chid, surf);
// Create the chamber
RPCChamber* ch = new RPCChamber (chid, surf);
// Add the rolls to rhe chamber
for(std::list<RPCRoll *>::iterator rl=rls.begin();
rl!=rls.end(); rl++){
ch->add(*rl);
}
for(auto rl : rls ) ch->add(rl);
// Add the chamber to the geometry
geometry->add(ch);

}
return geometry;
}



4 changes: 2 additions & 2 deletions Geometry/RPCGeometryBuilder/src/RPCGeometryBuilderFromCondDB.h
View file
Open in desktop
Expand Up @@ -11,7 +11,7 @@
#include <CondFormats/GeometryObjects/interface/RecoIdealGeometry.h>
#include <string>
#include <map>
#include <list>
#include <set>


class RPCGeometry;
Expand All @@ -30,7 +30,7 @@ class RPCGeometryBuilderFromCondDB


private:
std::map<RPCDetId,std::list<RPCRoll *> > chids;
std::map<RPCDetId,std::set<RPCRoll *> > chids;
Copy link
Contributor

Choose a reason for hiding this comment

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

randomizes the order of RPCRoll between jobs or IOVs

bool theComp11Flag;

};
Expand Down