Skip to content

Commit

Permalink
Merge branch 'WmassNanoProd_10_6_26' of https://github.com/WMass/cmssw …
Browse files Browse the repository at this point in the history 
…into WmassNanoProd_10_6_26
  • Loading branch information
@dbruschi
dbruschi committed Jan 2, 2023
2 parents 2c6a6e6 + a79cdfc commit 2cd9b66
Show file tree
Hide file tree
Showing 3 changed files with 138 additions and 109 deletions.
156 changes: 60 additions & 96 deletions Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@

#include "Alignment/CommonAlignment/interface/Utilities.h"


// #include "DataFormats/Math/interface/Point3D.h"


// #include "../interface/OffsetMagneticField.h"
Expand Down Expand Up @@ -377,8 +377,8 @@ ResidualGlobalCorrectionMakerBase::beginRun(edm::Run const& run, edm::EventSetup
// // const bool align2d = ispixel || isglued || isendcap;
// const DetId parmdetid = isglued ? DetId(gluedid) : det->geographicalId();

// const bool align2d = ispixel || isendcap;
const bool align2d = ispixel || isendcap || (alignGlued_ && isglued);
const bool align2d = ispixel || isendcap;
// const bool align2d = ispixel || isendcap || (alignGlued_ && isglued);
// const bool align2d = true;
// const bool align2d = ispixel;
// const bool align2d = false;
Expand All @@ -398,15 +398,15 @@ ResidualGlobalCorrectionMakerBase::beginRun(edm::Run const& run, edm::EventSetup
// parmset.emplace(1, det->geographicalId());
// }

parmset.emplace(0, aligndetid);
parmset.emplace(0, det->geographicalId());
parmset.emplace(2, aligndetid);
parmset.emplace(3, aligndetid);
parmset.emplace(4, aligndetid);
parmset.emplace(5, aligndetid);
parmset.emplace(5, det->geographicalId());

if (align2d) {
//local y alignment parameters only for pixels and wedge modules
parmset.emplace(1, aligndetid);
parmset.emplace(1, det->geographicalId());
}

// bfield and material parameters are associated to glued detids where applicable
Expand Down Expand Up @@ -558,6 +558,8 @@ ResidualGlobalCorrectionMakerBase::beginRun(edm::Run const& run, edm::EventSetup

const Vector3DBase<double, GlobalTag> dpos = surfaceDIdeal.position() - surfaceD.position();

// std::cout << "dpos: " << dpos << std::endl;

dx = dpos.x();
dy = dpos.y();
dz = dpos.z();
Expand Down Expand Up @@ -800,119 +802,81 @@ ResidualGlobalCorrectionMakerBase::beginRun(edm::Run const& run, edm::EventSetup
const DetId parmdetid = isglued ? DetId(gluedid) : det->geographicalId();
const GeomDet* parmDet = isglued ? globalGeometry->idToDet(parmdetid) : det;

const DetId aligndetid = alignGlued_ ? parmdetid : det->geographicalId();

const bool align2d = detidparms.count(std::make_pair(1, aligndetid));
const bool align2d = detidparms.count(std::make_pair(1, det->geographicalId()));

const Surface &surface = det->surface();

GloballyPositioned<double> surfaceD = surfaceToDouble(surface);

Matrix<double, 2, 2> Rglued = Matrix<double, 2, 2>::Identity();


//TODO restore alignment application functionality

if (isglued) {
GloballyPositioned<double> surfaceGlued = surfaceToDouble(parmDet->surface());

applyAlignment(surfaceGlued, parmdetid);

if (false) {
auto const dpos = surfaceD.position() - surfaceGlued.position();
const double gluedot = surfaceD.rotation().z()*surfaceGlued.rotation().z();
const double gluecross = surfaceD.rotation().z().cross(surfaceGlued.rotation().z()).mag();

std::cout << "surface pos:\n" << surfaceD.position() << std::endl;
std::cout << "surfaceglued pos:\n" << surfaceGlued.position() << std::endl;
std::cout << "gluedot = " << gluedot <<" gluecross = " << gluecross << " dpos:\n" << dpos << std::endl;
}

surfacemapD_[parmdetid] = surfaceGlued;

// recreate plane using relative position and orientation from ideal geometry, enforcing that the plane is parallel to the glued one (but preserving the relative orientation of the local z axis in case they are flipped)
const Surface &surfaceGluedIdealPre = globalGeometryIdeal->idToDet(parmDet->geographicalId())->surface();
const Surface &surfaceIdealPre = globalGeometryIdeal->idToDet(det->geographicalId())->surface();

const double zdot = surfaceGluedIdealPre.rotation().z()*surfaceIdealPre.rotation().z();
const double zsign = std::copysign(1.0, zdot);

const GloballyPositioned<double> surfaceGluedIdeal = surfaceToDouble(surfaceGluedIdealPre);
const GloballyPositioned<double> surfaceIdeal = surfaceToDouble(surfaceIdealPre, zsign*surfaceGluedIdeal.rotation().z());

if (false) {
// const GloballyPositioned<double> surfaceIdealOrig = surfaceToDouble(surfaceIdealPre);

std::cout << "zdot = " << zdot << std::endl;
}

const Vector3DBase<double, GlobalTag> uxpre(surfaceIdeal.rotation().x());
const Vector3DBase<double, GlobalTag> uypre(surfaceIdeal.rotation().y());
if (alignGlued_) {

//TODO apply partial alignment to surfaceGlued here

// recreate plane using relative position and orientation from ideal geometry, enforcing that the plane is parallel to the glued one (but preserving the relative orientation of the local z axis in case they are flipped)
// only the out-of-plane DOF's are preserved (ie the spacing of the planes), whereas the in-plane DOF's (position and orientation of local x and y axes) are left as-is from the nominal geometry
const Surface &surfaceGluedIdealPre = globalGeometryIdeal->idToDet(parmDet->geographicalId())->surface();
const Surface &surfaceIdealPre = globalGeometryIdeal->idToDet(det->geographicalId())->surface();

const Point3DBase<double, LocalTag> poslocal = surfaceGluedIdeal.toLocal(surfaceIdeal.position());
const double zdot = surfaceGluedIdealPre.rotation().z()*surfaceIdealPre.rotation().z();
const double zsign = std::copysign(1.0, zdot);

const Vector3DBase<double, LocalTag> uxlocal = surfaceGluedIdeal.toLocal(uxpre);
const Vector3DBase<double, LocalTag> uylocal = surfaceGluedIdeal.toLocal(uypre);
const GloballyPositioned<double> surfaceGluedIdeal = surfaceToDouble(surfaceGluedIdealPre);
const GloballyPositioned<double> surfaceIdeal = surfaceToDouble(surfaceIdealPre, zsign*surfaceGluedIdeal.rotation().z());

auto const poslocalideal = surfaceGluedIdeal.toLocal(surfaceIdeal.position());
auto const poslocalnom = surfaceGlued.toLocal(surfaceD.position());

const Point3DBase<double, LocalTag> poslocal(poslocalnom.x(), poslocalnom.y(), poslocalideal.z());

auto const posglobal = surfaceGlued.toGlobal(poslocal);

auto const uzglobal = zsign*surfaceGlued.rotation().z();

auto const &gy = surfaceD.rotation().y();
auto const &gz = uzglobal;

const Point3DBase<double, GlobalTag> posglobal = surfaceGlued.toGlobal(poslocal);
const Vector3DBase<double, GlobalTag> uxglobal = surfaceGlued.toGlobal(uxlocal);
const Vector3DBase<double, GlobalTag> uyglobal = surfaceGlued.toGlobal(uylocal);
auto const &uzglobal = zsign*surfaceGlued.rotation().z();
const Matrix<double, 3, 1> vy(gy.x(), gy.y(), gy.z());
const Matrix<double, 3, 1> vz(gz.x(), gz.y(), gz.z());

const TkRotation<double> tkrot(uxglobal.x(), uxglobal.y(), uxglobal.z(),
uyglobal.x(), uyglobal.y(), uyglobal.z(),
uzglobal.x(), uzglobal.y(), uzglobal.z());
const Matrix<double, 3, 1> &uz = vz;
const Matrix<double, 3, 1> ux = vy.cross(vz).normalized();
const Matrix<double, 3, 1> uy = uz.cross(ux);

surfaceD = GloballyPositioned<double>(posglobal, tkrot);
const TkRotation<double> tkrot(ux[0], ux[1], ux[2],
uy[0], uy[1], uy[2],
uz[0], uz[1], uz[2]);

if (false) {
const Point3DBase<double, LocalTag> posrel = surfaceGlued.toLocal(surfaceD.position());

const Vector3DBase<double, LocalTag> uxrel = surfaceGlued.toLocal(Vector3DBase<double, GlobalTag>(surfaceD.rotation().x()));
const Vector3DBase<double, LocalTag> uyrel = surfaceGlued.toLocal(Vector3DBase<double, GlobalTag>(surfaceD.rotation().y()));
const Vector3DBase<double, LocalTag> uzrel = surfaceGlued.toLocal(Vector3DBase<double, GlobalTag>(surfaceD.rotation().z()));

const Point3DBase<double, LocalTag> posrelideal = surfaceGluedIdeal.toLocal(surfaceIdeal.position());

const Vector3DBase<double, LocalTag> uxrelideal = surfaceGluedIdeal.toLocal(Vector3DBase<double, GlobalTag>(surfaceIdeal.rotation().x()));
const Vector3DBase<double, LocalTag> uyrelideal = surfaceGluedIdeal.toLocal(Vector3DBase<double, GlobalTag>(surfaceIdeal.rotation().y()));
const Vector3DBase<double, LocalTag> uzrelideal = surfaceGluedIdeal.toLocal(Vector3DBase<double, GlobalTag>(surfaceIdeal.rotation().z()));
surfaceD = GloballyPositioned<double>(posglobal, tkrot);

std::cout << "posrel: " << posrel << std::endl;
std::cout << "uxrel: " << uxrel << std::endl;
std::cout << "uyrel: " << uyrel << std::endl;
std::cout << "uzrel: " << uzrel << std::endl;

std::cout << "posrelideal: " << posrelideal << std::endl;
std::cout << "uxrelideal: " << uxrelideal << std::endl;
std::cout << "uyrelideal: " << uyrelideal << std::endl;
std::cout << "uzrelideal: " << uzrelideal << std::endl;
//TODO apply partial alignment to surfaceD here

const Vector3DBase<double, GlobalTag> uxglued(surfaceGlued.rotation().x());
const Vector3DBase<double, GlobalTag> uyglued(surfaceGlued.rotation().y());

const Vector3DBase<double, LocalTag> lxalt = surfaceD.toLocal(uxglued);
const Vector3DBase<double, LocalTag> lyalt = surfaceD.toLocal(uyglued);

// rotation matrix (jacobian dlocal/dglued)
Rglued(0, 0) = lxalt.x();
Rglued(0, 1) = lyalt.x();
Rglued(1, 0) = lxalt.y();
Rglued(1, 1) = lyalt.y();

}

// const GloballyPositioned<double> surfacemod(posglobal, tkrot);
// //
// std::cout << "surfaceD position:\n" << surfaceD.position() << "\nrotation:\n" << surfaceD.rotation() << std::endl;
// //
// std::cout << "surfacemod position:\n" << surfacemod.position() << "\nrotation:\n" << surfacemod.rotation() << std::endl;

const Vector3DBase<double, GlobalTag> uxglued(surfaceGlued.rotation().x());
const Vector3DBase<double, GlobalTag> uyglued(surfaceGlued.rotation().y());

const Vector3DBase<double, LocalTag> lxalt = surfaceD.toLocal(uxglued);
const Vector3DBase<double, LocalTag> lyalt = surfaceD.toLocal(uyglued);

// rotation matrix (jacobian dlocal/dglued)
Rglued(0, 0) = lxalt.x();
Rglued(0, 1) = lyalt.x();
Rglued(1, 0) = lxalt.y();
Rglued(1, 1) = lyalt.y();
}
else {
applyAlignment(surfaceD, det->geographicalId());
surfacemapD_[parmdetid] = surfaceGlued;
}

surfacemapD_[det->geographicalId()] = surfaceD;
rgluemap_[det->geographicalId()] = Rglued;


}

}
Expand Down
25 changes: 16 additions & 9 deletions Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerG4e.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -679,13 +679,13 @@ void ResidualGlobalCorrectionMakerG4e::produce(edm::Event &iEvent, const edm::Ev

const uint32_t gluedid = trackerTopology->glued(hit->geographicalId());
const bool isglued = gluedid != 0;
const DetId parmdetid = isglued ? DetId(gluedid) : hit->geographicalId();
// const DetId parmdetid = isglued ? DetId(gluedid) : hit->geographicalId();

const DetId aligndetid = alignGlued_ ? parmdetid : hit->geographicalId();
// const DetId aligndetid = alignGlued_ ? parmdetid : hit->geographicalId();


// const bool align2d = detidparms.count(std::make_pair(1, hit->geographicalId()));
const bool align2d = detidparms.count(std::make_pair(1, aligndetid));
const bool align2d = detidparms.count(std::make_pair(1, hit->geographicalId()));
// const bool align2d = detidparms.count(std::make_pair(1, aligndetid));
//
if (align2d) {
nvalidalign2d += 1;
Expand Down Expand Up @@ -1591,7 +1591,8 @@ void ResidualGlobalCorrectionMakerG4e::produce(edm::Event &iEvent, const edm::Ev
break;
}

const bool align2d = detidparms.count(std::make_pair(1, aligndetid));
const bool align2d = detidparms.count(std::make_pair(1, preciseHit->geographicalId()));
// const bool align2d = detidparms.count(std::make_pair(1, aligndetid));

const Matrix<double, 2, 2> &Rglued = rgluemap_.at(preciseHit->geographicalId());
const GloballyPositioned<double> &surfaceglued = surfacemapD_.at(parmdetid);
Expand Down Expand Up @@ -1783,6 +1784,9 @@ void ResidualGlobalCorrectionMakerG4e::produce(edm::Event &iEvent, const edm::Ev
const double localdydzval = localparmsalign[2];
const double localxval = localparmsalign[3];
const double localyval = localparmsalign[4];

const double localxvalorig = localparmsalignprop[3];
const double localyvalorig = localparmsalignprop[4];

//standard case

Expand All @@ -1803,15 +1807,16 @@ void ResidualGlobalCorrectionMakerG4e::produce(edm::Event &iEvent, const edm::Ev
// dy/dtheta_y
Aval(1,4) = -localxval*localdydzval;
// dx/dtheta_z
Aval(0,5) = localyval;
Aval(0,5) = localyvalorig;
// dy/dtheta_z
Aval(1,5) = -localxval;
Aval(1,5) = -localxvalorig;

// const Matrix<double, 2, 6> &A = alignGlued_ ? Rglued*Aval : Aval;

Matrix<double, 2, 6> A = R*Aval;
if (alignGlued_) {
A = R*Rglued*Aval;
// glued alignment dofs only for out-of-plance
A.middleCols<3>(2) = R*Rglued*Aval.middleCols<3>(2);
}

const unsigned int xresglobalidx = dores ? detidparms.at(std::make_pair(8, hit->geographicalId())) : 0;
Expand Down Expand Up @@ -1904,7 +1909,9 @@ void ResidualGlobalCorrectionMakerG4e::produce(edm::Event &iEvent, const edm::Ev
}

for (unsigned int idim=0; idim<nlocalalignment; ++idim) {
const unsigned int xglobalidx = detidparms.at(std::make_pair(alphaidxs[idim], aligndetid));
const unsigned int iidx = alphaidxs[idim];
const DetId ialigndetid = iidx > 1 && iidx < 5 ? aligndetid : preciseHit->geographicalId();
const unsigned int xglobalidx = detidparms.at(std::make_pair(iidx, ialigndetid));
globalidxv[iparm] = xglobalidx;
iparm++;
if (alphaidxs[idim]==0) {
Expand Down
Loading

0 comments on commit 2cd9b66

Please sign in to comment.