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

Run2-hcx192 Changes to DDHCalAngular by removing explicit reference to CLHEP #26271

Merged
merged 4 commits into from Apr 12, 2019
Merged

Run2-hcx192 Changes to DDHCalAngular by removing explicit reference to CLHEP #26271

Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b0eb2bf
Changes to DDHCalAngular by removing explicit reference to CLHEP
Mar 27, 2019
3c37ae7
Convert to degree only to find the name of rot matrix
Apr 4, 2019
28783fe
Use a function in DDTypes
Apr 4, 2019
ca03edf
Correct typing error
Apr 4, 2019
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
87 changes: 52 additions & 35 deletions Geometry/HcalAlgo/plugins/DDHCalAngular.cc
View file
Open in desktop
Expand Up @@ -7,13 +7,19 @@
#include <algorithm>

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Math/interface/GeantUnits.h"
#include "DetectorDescription/Core/interface/DDLogicalPart.h"
#include "DetectorDescription/Core/interface/DDCurrentNamespace.h"
#include "Geometry/HcalAlgo/plugins/DDHCalAngular.h"
#include "CLHEP/Units/GlobalSystemOfUnits.h"

//#define EDM_ML_DEBUG
using namespace geant_units;
using namespace geant_units::operators;

DDHCalAngular::DDHCalAngular() {
LogDebug("HCalGeom") << "DDHCalAngular test: Creating an instance";
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HCalGeom") << "DDHCalAngular: Creating an instance";
#endif
}

DDHCalAngular::~DDHCalAngular() {}
Expand All @@ -32,37 +38,44 @@ void DDHCalAngular::initialize(const DDNumericArguments & nArgs,
n = int (nArgs["n"]);
startCopyNo = int (nArgs["startCopyNo"]);
incrCopyNo = int (nArgs["incrCopyNo"]);
LogDebug("HCalGeom") << "DDHCalAngular debug: Parameters for positioning-- "
<< n << " copies in " << rangeAngle/CLHEP::deg
<< " from " << startAngle/CLHEP::deg << "\tShifts "
<< shiftX << ", " << shiftY
<< " along x, y axes; \tZoffest " << zoffset
<< "\tStart and inremental copy nos " << startCopyNo
<< ", " << incrCopyNo;

#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HCalGeom") << "DDHCalAngular: Parameters for positioning "
<< "-- " << n << " copies in "
<< convertRadToDeg(rangeAngle) << " from "
<< convertRadToDeg(startAngle) << "\tShifts "
<< shiftX << ", " << shiftY
<< " along x, y axes; \tZoffest " << zoffset
<< "\tStart and inremental copy nos "
<< startCopyNo << ", " << incrCopyNo;
#endif
rotns = sArgs["RotNameSpace"];
idNameSpace = DDCurrentNamespace::ns();
childName = sArgs["ChildName"];
#ifdef EDM_ML_DEBUG
DDName parentName = parent().name();
LogDebug("HCalGeom") << "DDHCalAngular debug: Parent " << parentName
<< "\tChild " << childName << "\tNameSpace "
<< idNameSpace << "\tRotation Namespace " << rotns;
edm::LogVerbatim("HCalGeom") << "DDHCalAngular: Parent " << parentName
<< "\tChild " << childName << "\tNameSpace "
<< idNameSpace << "\tRotation Namespace "
<< rotns;
#endif
}

void DDHCalAngular::execute(DDCompactView& cpv) {

double dphi = rangeAngle/n;
double phi = startAngle;
double phix = startAngle;
int copyNo = startCopyNo;
double theta = 90._deg;

for (int ii=0; ii<n; ii++) {

double phideg = phi/CLHEP::deg;
int iphi;
if (phideg > 0) iphi = int(phideg+0.1);
else iphi = int(phideg-0.1);
if (iphi >= 360) iphi -= 360;
phideg = iphi;
double phideg = convertRadToDeg(phix);
int iphi = std::lround(phideg);
if (iphi >= 360) {
iphi -= 360;
phideg = iphi;
phix = convertDegToRad(phideg);
bsunanda marked this conversation as resolved.
Show resolved Hide resolved
}
double phiy = phix + 90._deg;
DDRotation rotation;
std::string rotstr("NULL");

Expand All @@ -72,26 +85,30 @@ void DDHCalAngular::execute(DDCompactView& cpv) {
rotstr = rotstr + std::to_string(phideg);
rotation = DDRotation(DDName(rotstr, rotns));
if (!rotation) {
LogDebug("HCalGeom") << "DDHCalAngular test: Creating a new rotation "
<< DDName(rotstr, idNameSpace) << "\t90, "
<< phideg << ", 90, " << (phideg+90) << ", 0, 0";
rotation = DDrot(DDName(rotstr, rotns), 90*CLHEP::deg,
phideg*CLHEP::deg, 90*CLHEP::deg,
(90+phideg)*CLHEP::deg, 0*CLHEP::deg, 0*CLHEP::deg);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HCalGeom") << "DDHCalAngular: Creating a rotation "
<< DDName(rotstr, idNameSpace) << "\t90, "
<< phideg << ", 90, " << (phideg+90)
<< ", 0, 0";
#endif
rotation = DDrot(DDName(rotstr, rotns), theta, phix, theta, phiy, 0,0);
}
}

double xpos = shiftX*cos(phi) - shiftY*sin(phi);
double ypos = shiftX*sin(phi) + shiftY*cos(phi);
double xpos = shiftX*cos(phix) - shiftY*sin(phix);
Copy link
Contributor

Choose a reason for hiding this comment

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

otherwise, this is only equivalent to the old implementation for integer values of convertRadToDeg(startAngle)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I have made some tests with new vs old implementation and then test new==old. You see this test fail in many cases. Theta is 90 degree. pi and twopi match. All iphi's are the same but phideg or phi (phix) do not in many cases.

Theta 1.5708:1.5708:1
Pi 3.14159:3.14159:1 TwoPi 6.28319:6.28319:1
phideg 270:270:1 iphi 270:270 phi 4.71239:4.71239:1
phideg 300:300:1 iphi 300:300 phi 5.23599:5.23599:1
phideg 330:330:0 iphi 330:330 phi 5.75959:5.75959:0
phideg 0:0:1 iphi 0:0 phi 0:0:1
phideg 30:30:0 iphi 30:30 phi 0.523599:0.523599:1
phideg 60:60:0 iphi 60:60 phi 1.0472:1.0472:1
phideg 90:90:1 iphi 90:90 phi 1.5708:1.5708:1
phideg 120:120:0 iphi 120:120 phi 2.0944:2.0944:1
phideg 150:150:0 iphi 150:150 phi 2.61799:2.61799:0
phideg 180:180:0 iphi 180:180 phi 3.14159:3.14159:0
phideg 210:210:0 iphi 210:210 phi 3.66519:3.66519:0
phideg 240:240:0 iphi 240:240 phi 4.18879:4.18879:1

So my conclusion is that these are all rounding effects and once there is random # sequence goes crazy.

Copy link
Contributor

Choose a reason for hiding this comment

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

@bsunanda: Could you please explain the print-out above with phideg, iphi, etc.? I don't see any differences.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The first # is the new value, the second is the old value and the third is (old == new). I did not put (old==new) for int32 like iphi. But in many cases old==new fails and shows 0 which is due to some rounding issues. This is the culprit for change is random # sequence. This always causes difference for small samples but statistically they balance out

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Frankly I am now not worried about this difference

Copy link
Contributor

Choose a reason for hiding this comment

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

It would be helpful to see (old - new) when they differ to understand the scale of the differences.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The difference is at the level of e-14. Here they are. For each, there are the new value, old value, difference and test of new==old results
Theta 1.5708:1.5708:0:1
Pi 3.14159:3.14159:0:1 TwoPi 6.28319:6.28319:0:1
phideg 270:270:0:1 iphi 270:270:0:1 phi 4.71239:4.71239:0:1
phideg 300:300:0:1 iphi 300:300:0:1 phi 5.23599:5.23599:0:1
phideg 330:330:5.68434e-14:0 iphi 330:330:0:1 phi 5.75959:5.75959:8.88178e-16:0
phideg 0:0:0:1 iphi 0:0:0:1 phi 0:0:0:1
phideg 30:30:-3.55271e-15:0 iphi 30:30:0:1 phi 0.523599:0.523599:0:1
phideg 60:60:-7.10543e-15:0 iphi 60:60:0:1 phi 1.0472:1.0472:0:1
phideg 90:90:0:1 iphi 90:90:0:1 phi 1.5708:1.5708:0:1
phideg 120:120:-1.42109e-14:0 iphi 120:120:0:1 phi 2.0944:2.0944:0:1
phideg 150:150:-2.84217e-14:0 iphi 150:150:0:1 phi 2.61799:2.61799:-4.44089e-16:0
phideg 180:180:-2.84217e-14:0 iphi 180:180:0:1 phi 3.14159:3.14159:-4.44089e-16:0
phideg 210:210:-2.84217e-14:0 iphi 210:210:0:1 phi 3.66519:3.66519:-8.88178e-16:0
phideg 240:240:-2.84217e-14:0 iphi 240:240:0:1 phi 4.18879:4.18879:0:1

double ypos = shiftX*sin(phix) + shiftY*cos(phix);
DDTranslation tran(xpos, ypos, zoffset);

DDName parentName = parent().name();
cpv.position(DDName(childName,idNameSpace), parentName, copyNo, tran, rotation);
LogDebug("HCalGeom") << "DDHCalAngular test: "
<< DDName(childName, idNameSpace) << " number "
<< copyNo << " positioned in " << parentName << " at "
<< tran << " with " << rotation;
phi += dphi;
cpv.position(DDName(childName,idNameSpace),parentName,copyNo,tran,rotation);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HCalGeom") << "DDHCalAngular: "
<< DDName(childName, idNameSpace)
<< " number " << copyNo << " positioned in "
<< parentName << " at " << tran << " with "
<< rotation;
#endif
phix += dphi;
copyNo += incrCopyNo;
}
}