slicPandora/src
diff -u -r1.15 -r1.16
--- CalorimeterHitProcessor.cpp 24 Nov 2010 16:44:08 -0000 1.15
+++ CalorimeterHitProcessor.cpp 26 Nov 2010 14:52:27 -0000 1.16
@@ -1,4 +1,10 @@
-// $Id: CalorimeterHitProcessor.cpp,v 1.15 2010/11/24 16:44:08 speckmay Exp $
+// $Id: CalorimeterHitProcessor.cpp,v 1.16 2010/11/26 14:52:27 speckmay Exp $
+
+// variable names change from cdr to trunk version. use the new names, and keep the transformation of the names in defines
+#define m_nCellRadiationLengths m_nRadiationLengths
+#define m_nCellInteractionLengths m_nInteractionLengths
+#define m_cellNormalVector m_normalVector
+#define POST_CDR_VERSION 1
#include "CalorimeterHitProcessor.h"
@@ -142,6 +148,10 @@
// Position in mm.
params.m_positionVector = pandora::CartesianVector(pos[0], pos[1], pos[2]);
+ #ifdef POST_CDR_VERSION
+ params.m_expectedDirection = pandora::CartesianVector(pos[0], pos[1], pos[2]).GetUnitVector();
+ #endif
+
// Digital calorimeter setting.
params.m_isDigital = xsubdet->m_isDigital;
@@ -199,10 +209,10 @@
params.m_cellThickness = xlayerParams.m_cellThickness;
// Number of radiation lengths in layer.
- params.m_nRadiationLengths = layerParams.m_nRadiationLengths;
+ params.m_nCellRadiationLengths = layerParams.m_nRadiationLengths;
// Number of interaction lengths in layer.
- params.m_nInteractionLengths = layerParams.m_nInteractionLengths;
+ params.m_nCellInteractionLengths = layerParams.m_nInteractionLengths;
// Type of hit.
params.m_hitType = xsubdet->m_inputHitType;
@@ -213,16 +223,20 @@
// Barrel case uses module.
if (xsubdet->m_inputDetectorRegion.Get() == pandora::BARREL)
{
- params.m_normalVector = *(xsubdet->m_normalVectors.at(module));
+ params.m_cellNormalVector = *(xsubdet->m_normalVectors.at(module));
}
// Endcap case is computed.
else
{
- params.m_normalVector = pandora::CartesianVector(0., 0., pos[2]/fabs(pos[2]));
+ params.m_cellNormalVector = pandora::CartesianVector(0., 0., pos[2]/fabs(pos[2]));
}
// Number of interaction lengths from the IP.
- params.m_nInteractionLengthsFromIp = 0.f;
+ float intLengthsFromIp = 0.f;
+ float radLengthsFromIp = 0.f;
+ GetPathLengths(hit, radLengthsFromIp, intLengthsFromIp);
+ params.m_nInteractionLengthsFromIp = intLengthsFromIp;
+ params.m_nRadiationLengthsFromIp = radLengthsFromIp;
// Return completed calorimeter parameters.
return params;
@@ -237,10 +251,481 @@
std::cout << " time: " << params.m_time.Get() << " [ns]" << std::endl;
std::cout << " cell sizes: " << params.m_cellSizeU.Get() << ", " << params.m_cellSizeV.Get() << " [mm]" << std::endl;
std::cout << " cell thick: " << params.m_cellThickness.Get() << " [mm]" << std::endl;
- std::cout << " radLengths: " << params.m_nRadiationLengths.Get() << std::endl;
- std::cout << " intLengths: " << params.m_nInteractionLengths.Get() << std::endl;
- std::cout << " normalVec: " << params.m_normalVector.Get() << std::endl;
+ std::cout << " radLengths: " << params.m_nCellRadiationLengths.Get() << std::endl;
+ std::cout << " intLengths: " << params.m_nCellInteractionLengths.Get() << std::endl;
+ std::cout << " normalVec: " << params.m_cellNormalVector.Get() << std::endl;
std::cout << " hitType: " << params.m_hitType.Get() << std::endl;
std::cout << " detectorRegion: " << params.m_detectorRegion.Get() << std::endl;
std::cout << " outerLayer: " << params.m_isInOuterSamplingLayer.Get() << std::endl;
}
+
+
+
+void CalorimeterHitProcessor::GetPathLengths(const EVENT::CalorimeterHit *const pCaloHit, float &nRadiationLengthsFromIp,
+ float &nInteractionLengthsFromIp)
+{
+ try
+ {
+ static bool initialized = false;
+
+ // Coordinates of the sub-detectors in one quadrant
+ static float rMinECalBarrel = 0.f;
+ static float rMaxECalBarrel = 0.f;
+ static float zMinECalBarrel = 0.f;
+ static float zMaxECalBarrel = 0.f;
+ static float innerPhiECalBarrel = 0.f;
+ static int innerSymECalBarrel = 0;
+ static float outerPhiECalBarrel = 0.f;
+ static int outerSymECalBarrel = 0;
+
+ static float rMinECalEndCap = 0.f;
+ static float rMaxECalEndCap = 0.f;
+ static float zMinECalEndCap = 0.f;
+ static float zMaxECalEndCap = 0.f;
+ static float innerPhiECalEndCap = 0.f;
+ static int innerSymECalEndCap = 0;
+ static float outerPhiECalEndCap = 0.f;
+ static int outerSymECalEndCap = 0;
+
+ static float rMinHCalBarrel = 0.f;
+ static float rMaxHCalBarrel = 0.f;
+ static float zMinHCalBarrel = 0.f;
+ static float zMaxHCalBarrel = 0.f;
+ static float innerPhiHCalBarrel = 0.f;
+ static int innerSymHCalBarrel = 0;
+ static float outerPhiHCalBarrel = 0.f;
+ static int outerSymHCalBarrel = 0;
+
+ static float rMinHCalEndCap = 0.f;
+ static float rMaxHCalEndCap = 0.f;
+ static float zMinHCalEndCap = 0.f;
+ static float zMaxHCalEndCap = 0.f;
+ static float innerPhiHCalEndCap = 0.f;
+ static int innerSymHCalEndCap = 0;
+ static float outerPhiHCalEndCap = 0.f;
+ static int outerSymHCalEndCap = 0;
+
+ static float rMinCoil = 0.f;
+ static float rMaxCoil = 0.f;
+ static float zMinCoil = 0.f;
+ static float zMaxCoil = 0.f;
+// static float innerPhiCoil = 0.f;
+// static int innerSymCoil = 0;
+// static float outerPhiCoil = 0.f;
+// static int outerSymCoil = 0;
+
+ static float rMinTracker = 0.f;
+ static float rMaxTracker = 0.f;
+ static float zMinTracker = 0.f;
+ static float zMaxTracker = 0.f;
+// static float innerPhiTracker = 0.f;
+// static int innerSymTracker = 0;
+// static float outerPhiTracker = 0.f;
+// static int outerSymTracker = 0;
+
+ static float rMinMuonBarrel = 0.f;
+ static float rMaxMuonBarrel = 0.f;
+ static float zMinMuonBarrel = 0.f;
+ static float zMaxMuonBarrel = 0.f;
+ static float innerPhiMuonBarrel = 0.f;
+ static int innerSymMuonBarrel = 0;
+ static float outerPhiMuonBarrel = 0.f;
+ static int outerSymMuonBarrel = 0;
+
+ static float rMinMuonEndCap = 0.f;
+ static float rMaxMuonEndCap = 0.f;
+ static float zMinMuonEndCap = 0.f;
+ static float zMaxMuonEndCap = 0.f;
+ static float innerPhiMuonEndCap = 0.f;
+ static int innerSymMuonEndCap = 0;
+ static float outerPhiMuonEndCap = 0.f;
+ static int outerSymMuonEndCap = 0;
+
+ static float avgIntLengthECalBarrel = 0.f;
+ static float avgIntLengthECalEndCap = 0.f;
+ static float avgIntLengthHCalBarrel = 0.f;
+ static float avgIntLengthHCalEndCap = 0.f;
+ static float avgIntLengthCoil = 0.f;
+// static float avgIntLengthTracker = 0.f;
+ static float avgIntLengthMuonBarrel = 0.f;
+ static float avgIntLengthMuonEndCap = 0.f;
+
+ static float avgRadLengthECalBarrel = 0.f;
+ static float avgRadLengthECalEndCap = 0.f;
+ static float avgRadLengthHCalBarrel = 0.f;
+ static float avgRadLengthHCalEndCap = 0.f;
+ static float avgRadLengthCoil = 0.f;
+// static float avgRadLengthTracker = 0.f;
+ static float avgRadLengthMuonBarrel = 0.f;
+ static float avgRadLengthMuonEndCap = 0.f;
+
+ if (!initialized)
+ {
+ const PandoraApi::Geometry::Parameters geometryParameters;
+
+ DetectorGeometry* geom = getJobManager()->getDetectorGeometry();
+ JobConfig::CalorimeterTypes calTypes = getJobManager()->getJobConfig()->getCalorimeterTypes();
+
+ // no path length info for tracker computed yet
+ rMinTracker = 0.;
+ rMaxTracker = 0.;
+ zMinTracker = 0;
+ zMaxTracker = 0.;
+
+ // coil
+ rMinCoil = geom->getGeometryParameters()->m_coilInnerRadius.Get();
+ rMaxCoil = geom->getGeometryParameters()->m_coilOuterRadius.Get();
+ zMinCoil = 0.f;
+ zMaxCoil = geom->getGeometryParameters()->m_coilZExtent.Get();
+ avgRadLengthCoil = geom->m_coilRadiationLength;
+ avgIntLengthCoil = geom->m_coilInteractionLength;
+
+
+ // Get the EM barrel SubDetector parameters.
+ std::string calType = "EM_BARREL";
+ PandoraApi::GeometryParameters::SubDetectorParameters* subdet = geom->getSubDetectorFromType(calType);
+ DetectorGeometry::ExtraSubDetectorParameters* xsubdet = geom->getExtraSubDetectorParametersFromType(calType);
+
+ ComputeAverageLengthsFromExtraSubDetectorParameters(xsubdet, avgIntLengthECalBarrel, avgRadLengthECalBarrel);
+
+ rMinECalBarrel = subdet->m_innerRCoordinate.Get();
+ rMaxECalBarrel = subdet->m_outerRCoordinate.Get();
+ zMinECalBarrel = subdet->m_innerZCoordinate.Get();
+ zMaxECalBarrel = subdet->m_outerZCoordinate.Get();
+ innerPhiECalBarrel = subdet->m_innerPhiCoordinate.Get();
+ outerPhiECalBarrel = subdet->m_outerPhiCoordinate.Get();
+ innerSymECalBarrel = subdet->m_innerSymmetryOrder.Get();
+ outerSymECalBarrel = subdet->m_outerSymmetryOrder.Get();
+
+ // Get the HCal barrel SubDetector parameters.
+ calType = "HAD_BARREL";
+ subdet = geom->getSubDetectorFromType(calType);
+ xsubdet = geom->getExtraSubDetectorParametersFromType(calType);
+
+ ComputeAverageLengthsFromExtraSubDetectorParameters(xsubdet, avgIntLengthHCalBarrel, avgRadLengthHCalBarrel);
+
+ rMinHCalBarrel = subdet->m_innerRCoordinate.Get();
+ rMaxHCalBarrel = subdet->m_outerRCoordinate.Get();
+ zMinHCalBarrel = subdet->m_innerZCoordinate.Get();
+ zMaxHCalBarrel = subdet->m_outerZCoordinate.Get();
+ innerPhiHCalBarrel = subdet->m_innerPhiCoordinate.Get();
+ outerPhiHCalBarrel = subdet->m_outerPhiCoordinate.Get();
+ innerSymHCalBarrel = subdet->m_innerSymmetryOrder.Get();
+ outerSymHCalBarrel = subdet->m_outerSymmetryOrder.Get();
+
+
+ // Get the ECal EndCap SubDetector parameters.
+ calType = "EM_ENDCAP";
+ subdet = geom->getSubDetectorFromType(calType);
+ xsubdet = geom->getExtraSubDetectorParametersFromType(calType);
+
+ ComputeAverageLengthsFromExtraSubDetectorParameters(xsubdet, avgIntLengthECalEndCap, avgRadLengthECalEndCap);
+
+ rMinECalEndCap = subdet->m_innerRCoordinate.Get();
+ rMaxECalEndCap = subdet->m_outerRCoordinate.Get();
+ zMinECalEndCap = subdet->m_innerZCoordinate.Get();
+ zMaxECalEndCap = subdet->m_outerZCoordinate.Get();
+ innerPhiECalEndCap = subdet->m_innerPhiCoordinate.Get();
+ outerPhiECalEndCap = subdet->m_outerPhiCoordinate.Get();
+ innerSymECalEndCap = subdet->m_innerSymmetryOrder.Get();
+ outerSymECalEndCap = subdet->m_outerSymmetryOrder.Get();
+
+
+ // Get the HCal EndCap SubDetector parameters.
+ calType = "HAD_ENDCAP";
+ subdet = geom->getSubDetectorFromType(calType);
+ xsubdet = geom->getExtraSubDetectorParametersFromType(calType);
+
+ ComputeAverageLengthsFromExtraSubDetectorParameters(xsubdet, avgIntLengthHCalEndCap, avgRadLengthHCalEndCap);
+
+ rMinHCalEndCap = subdet->m_innerRCoordinate.Get();
+ rMaxHCalEndCap = subdet->m_outerRCoordinate.Get();
+ zMinHCalEndCap = subdet->m_innerZCoordinate.Get();
+ zMaxHCalEndCap = subdet->m_outerZCoordinate.Get();
+ innerPhiHCalEndCap = subdet->m_innerPhiCoordinate.Get();
+ outerPhiHCalEndCap = subdet->m_outerPhiCoordinate.Get();
+ innerSymHCalEndCap = subdet->m_innerSymmetryOrder.Get();
+ outerSymHCalEndCap = subdet->m_outerSymmetryOrder.Get();
+
+ // Get the Muon Barrel SubDetector parameters.
+ calType = "MUON_BARREL";
+ subdet = geom->getSubDetectorFromType(calType);
+ xsubdet = geom->getExtraSubDetectorParametersFromType(calType);
+
+ ComputeAverageLengthsFromExtraSubDetectorParameters(xsubdet, avgIntLengthMuonBarrel, avgRadLengthMuonBarrel);
+
+ rMinMuonBarrel = subdet->m_innerRCoordinate.Get();
+ rMaxMuonBarrel = subdet->m_outerRCoordinate.Get();
+ zMinMuonBarrel = subdet->m_innerZCoordinate.Get();
+ zMaxMuonBarrel = subdet->m_outerZCoordinate.Get();
+ innerPhiMuonBarrel = subdet->m_innerPhiCoordinate.Get();
+ outerPhiMuonBarrel = subdet->m_outerPhiCoordinate.Get();
+ innerSymMuonBarrel = subdet->m_innerSymmetryOrder.Get();
+ outerSymMuonBarrel = subdet->m_outerSymmetryOrder.Get();
+
+ // Get the Muon EndCap SubDetector parameters.
+ calType = "MUON_ENDCAP";
+ subdet = geom->getSubDetectorFromType(calType);
+ xsubdet = geom->getExtraSubDetectorParametersFromType(calType);
+
+ ComputeAverageLengthsFromExtraSubDetectorParameters(xsubdet, avgIntLengthMuonEndCap, avgRadLengthMuonEndCap);
+
+ rMinMuonEndCap = subdet->m_innerRCoordinate.Get();
+ rMaxMuonEndCap = subdet->m_outerRCoordinate.Get();
+ zMinMuonEndCap = subdet->m_innerZCoordinate.Get();
+ zMaxMuonEndCap = subdet->m_outerZCoordinate.Get();
+ innerPhiMuonEndCap = subdet->m_innerPhiCoordinate.Get();
+ outerPhiMuonEndCap = subdet->m_outerPhiCoordinate.Get();
+ innerSymMuonEndCap = subdet->m_innerSymmetryOrder.Get();
+ outerSymMuonEndCap = subdet->m_outerSymmetryOrder.Get();
+
+ initialized = true;
+ }
+
+ const float *pPosition = pCaloHit->getPosition();
+ pandora::CartesianVector positionVector(pPosition[0], pPosition[1], pPosition[2]);
+
+ float radius, phi, z;
+ positionVector.GetCylindricalCoordinates(radius, phi, z);
+ positionVector.SetValues(radius, 0.f, std::fabs(z));
+
+ float corrRMin = 0.f, corrRMax = 0.f;
+
+ ComputePolygonCorrection(innerSymECalBarrel, rMinECalBarrel, true, innerPhiECalBarrel, phi, corrRMin);
+ ComputePolygonCorrection(innerSymECalBarrel, rMaxECalBarrel, true, outerPhiECalBarrel, phi, corrRMax);
+ const float eCalBarrelPathLength(ComputePathLengthFromIPInRectangle(positionVector, corrRMin, zMinECalBarrel, corrRMax, zMaxECalBarrel));
+
+ ComputePolygonCorrection(innerSymHCalBarrel, rMinHCalBarrel, true, innerPhiHCalBarrel, phi, corrRMin);
+ ComputePolygonCorrection(innerSymHCalBarrel, rMaxHCalBarrel, true, outerPhiHCalBarrel, phi, corrRMax);
+ const float hCalBarrelPathLength(ComputePathLengthFromIPInRectangle(positionVector, corrRMin, zMinHCalBarrel, corrRMax, zMaxHCalBarrel));
+
+ ComputePolygonCorrection(innerSymMuonBarrel, rMinMuonBarrel, true, innerPhiMuonBarrel, phi, corrRMin);
+ ComputePolygonCorrection(innerSymMuonBarrel, rMaxMuonBarrel, true, outerPhiMuonBarrel, phi, corrRMax);
+ const float muonBarrelPathLength(ComputePathLengthFromIPInRectangle(positionVector, corrRMin, zMinMuonBarrel, corrRMax, zMaxMuonBarrel));
+
+ ComputePolygonCorrection(innerSymECalEndCap, rMinECalEndCap, true, innerPhiECalEndCap, phi, corrRMin);
+ ComputePolygonCorrection(innerSymECalEndCap, rMaxECalEndCap, true, outerPhiECalEndCap, phi, corrRMax);
+ const float eCalEndCapPathLength(ComputePathLengthFromIPInRectangle(positionVector, corrRMin, zMinECalEndCap, corrRMax, zMaxECalEndCap));
+
+ ComputePolygonCorrection(innerSymHCalEndCap, rMinHCalEndCap, true, innerPhiHCalEndCap, phi, corrRMin);
+ ComputePolygonCorrection(innerSymHCalEndCap, rMaxHCalEndCap, true, outerPhiHCalEndCap, phi, corrRMax);
+ const float hCalEndCapPathLength(ComputePathLengthFromIPInRectangle(positionVector, corrRMin, zMinHCalEndCap, corrRMax, zMaxHCalEndCap));
+
+ ComputePolygonCorrection(innerSymMuonEndCap, rMinMuonEndCap, true, innerPhiMuonEndCap, phi, corrRMin);
+ ComputePolygonCorrection(innerSymMuonEndCap, rMaxMuonEndCap, true, outerPhiMuonEndCap, phi, corrRMax);
+ const float muonEndCapPathLength(ComputePathLengthFromIPInRectangle(positionVector, corrRMin, zMinMuonEndCap, corrRMax, zMaxMuonEndCap));
+
+ const float coilPathLength(ComputePathLengthFromIPInRectangle(positionVector, rMinCoil, zMinCoil, rMaxCoil, zMaxCoil));
+
+ nRadiationLengthsFromIp = (eCalBarrelPathLength * avgRadLengthECalBarrel) +
+ (hCalBarrelPathLength * avgRadLengthHCalBarrel) +
+ (muonBarrelPathLength * avgRadLengthMuonBarrel) +
+ (eCalEndCapPathLength * avgRadLengthECalEndCap) +
+ (hCalEndCapPathLength * avgRadLengthHCalEndCap) +
+ (muonEndCapPathLength * avgRadLengthMuonEndCap) +
+ (coilPathLength * avgRadLengthCoil);
+
+ nInteractionLengthsFromIp = (eCalBarrelPathLength * avgIntLengthECalBarrel) +
+ (hCalBarrelPathLength * avgIntLengthHCalBarrel) +
+ (muonBarrelPathLength * avgIntLengthMuonBarrel) +
+ (eCalEndCapPathLength * avgIntLengthECalEndCap) +
+ (hCalEndCapPathLength * avgIntLengthHCalEndCap) +
+ (muonEndCapPathLength * avgIntLengthMuonEndCap) +
+ (coilPathLength * avgIntLengthCoil);
+ }
+ catch (...)
+ {
+ std::cout << "CalorimeterHitProcessor: calculation of path length failed" << std::endl;
+ throw;
+ }
+}
+
+
+void CalorimeterHitProcessor::ComputeAverageLengthsFromExtraSubDetectorParameters(DetectorGeometry::ExtraSubDetectorParameters* xsubdet, float& avgIntLength, float& avgRadLength)
+{
+ // compute the average interaction and radiation lengths
+ float sumCellThicknesses = 0.f;
+ avgIntLength = 0.f;
+ avgRadLength = 0.f;
+ for (DetectorGeometry::ExtraLayerParametersList::iterator itLay = xsubdet->m_extraLayerParams.begin(), itLayEnd = xsubdet->m_extraLayerParams.end(); itLay != itLayEnd; ++itLay )
+ {
+ DetectorGeometry::ExtraLayerParameters& layParameters = (*itLay);
+ sumCellThicknesses += layParameters.m_cellThickness.Get();
+ avgIntLength += layParameters.m_intLength;
+ avgRadLength += layParameters.m_radLength;
+ }
+ avgIntLength = (sumCellThicknesses>0.f ? avgIntLength / sumCellThicknesses: 0.f);
+ avgRadLength = (sumCellThicknesses>0.f ? avgRadLength / sumCellThicknesses: 0.f);
+}
+
+
+float CalorimeterHitProcessor::ComputePathLengthFromIPInRectangle(const pandora::CartesianVector &position, float rMin, float zMin,
+ float rMax, float zMax)
+{
+ // compute cuts with rectangle borders
+ float phi, radius, z;
+ position.GetCylindricalCoordinates(radius, phi, z);
+
+ bool valid[4];
+ float xInt[4], zInt[4];
+
+ valid[0] = IntersectLines2D(0.f, 0.f, radius, z, rMin, zMin, rMin, zMax, xInt[0], zInt[0]); // first edge of rectangle at rMin
+ valid[1] = IntersectLines2D(0.f, 0.f, radius, z, rMin, zMin, rMax, zMin, xInt[1], zInt[1]); // first edge of rectangle at zMin
+ valid[2] = IntersectLines2D(0.f, 0.f, radius, z, rMax, zMax, rMax, zMin, xInt[2], zInt[2]); // first edge of rectangle at rMax
+ valid[3] = IntersectLines2D(0.f, 0.f, radius, z, rMax, zMax, rMin, zMax, xInt[3], zInt[3]); // first edge of rectangle at zMax
+
+ int indexFirstPoint = -1;
+ int indexSecondPoint = -1;
+
+ for (int i = 0; i < 4; ++i)
+ {
+ if (valid[i])
+ {
+ if (indexFirstPoint == -1)
+ {
+ indexFirstPoint = i;
+ }
+ else if (indexSecondPoint == -1)
+ {
+ indexSecondPoint = i;
+ }
+ else
+ {
+ std::cout << "ERROR calohitcreator problem at computing path length from IP to calohit in rectangle" << std::endl;
+ }
+ }
+ }
+
+ if (indexFirstPoint == -1)
+ return 0.f;
+
+ pandora::CartesianVector intersectionA( xInt[indexFirstPoint], 0.f, zInt[indexFirstPoint] );
+
+ if (indexSecondPoint == -1)
+ {
+ const float length(pandora::CartesianVector(position - intersectionA).GetMagnitude());
+ return length;
+ }
+
+ pandora::CartesianVector intersectionB(xInt[indexSecondPoint], 0.f, zInt[indexSecondPoint]);
+
+ const float length(pandora::CartesianVector(intersectionA - intersectionB).GetMagnitude());
+ return length;
+}
+
+
+
+
+bool CalorimeterHitProcessor::IntersectLines2D(float lineAXStart, float lineAYStart, float lineAXEnd, float lineAYEnd,
+ float lineBXStart, float lineBYStart, float lineBXEnd, float lineBYEnd, float &xIntersect, float &yIntersect)
+{
+ // Slopes of the two lines, take max float value instead of infinity
+ float k0(std::numeric_limits<float>::max()), k1(std::numeric_limits<float>::max());
+
+ bool parallelToY_A = false;
+ bool parallelToY_B = false;
+
+ if (std::fabs(lineAXEnd-lineAXStart) > std::numeric_limits<float>::epsilon())
+ {
+ k0 = (lineAYEnd - lineAYStart) / (lineAXEnd - lineAXStart);
+ }
+ else
+ {
+ parallelToY_A = true;
+ }
+
+ if (std::fabs(lineBXEnd - lineBXStart) > std::numeric_limits<float>::epsilon())
+ {
+ k1 = (lineBYEnd - lineBYStart) / (lineBXEnd - lineBXStart);
+ }
+ else
+ {
+ parallelToY_B = true;
+ }
+
+ if (parallelToY_A && parallelToY_B)
+ {
+ xIntersect = 0.f;
+ yIntersect = 0.f;
+
+ return false;
+ }
+
+ if (parallelToY_A)
+ {
+ xIntersect = lineAXStart;
+ yIntersect = lineBXStart + k1 * (xIntersect - lineBXStart);
+ }
+ else if (parallelToY_B)
+ {
+ xIntersect = lineBXStart;
+ yIntersect = lineAXStart + k0 * (xIntersect - lineAXStart);
+ }
+ else
+ {
+ const float b0 = -1;
+ const float b1 = -1;
+
+ const float c0 = (lineAYStart - k0 * lineAXStart);
+ const float c1 = (lineBYStart - k1 * lineBXStart);
+
+ const float determinant(k0 * b1 - k1 * b0);
+
+ if (0.f == determinant)
+ {
+ std::cout << "ERROR zero determinant in interaction length calculator" << std::endl;
+ return false;
+ }
+
+ // use Kramers rule to compute xi and yi
+ xIntersect=((b0 * c1 - b1 * c0) / determinant);
+ yIntersect=((k1 * c0 - k0 * c1) / determinant);
+ }
+
+ // check if intersections are within end-points of lines
+ // check x coordinate, line A
+ if (!((xIntersect >= lineAXStart && xIntersect <= lineAXEnd) || (xIntersect >= lineAXEnd && xIntersect <= lineAXStart)))
+ return false;
+
+ // check x coordinate, line B
+ if (!((xIntersect >= lineBXStart && xIntersect <= lineBXEnd) || (xIntersect >= lineBXEnd && xIntersect <= lineBXStart)))
+ return false;
+
+ // check y coordinate, line A
+ if (!((yIntersect >= lineAYStart && yIntersect <= lineAYEnd) || (yIntersect >= lineAYEnd && yIntersect <= lineAYStart)))
+ return false;
+
+ // check y coordinate, line B
+ if (!((yIntersect >= lineBYStart && yIntersect <= lineBYEnd) || (yIntersect >= lineBYEnd && yIntersect <= lineBYStart)))
+ return false;
+
+ return true;
+}
+
+
+void CalorimeterHitProcessor::ComputePolygonCorrection(int symmetryOrder, float distanceToIp, bool isClosestDistance, float phi0, float phi, float& correctedDistance)
+{
+ correctedDistance = distanceToIp;
+ if (symmetryOrder > 2)
+ {
+ static const float pi(std::acos(-1.));
+
+ const float symmOrderAngle = (pi/float(symmetryOrder));
+ float rotPhi = phi - trunc(phi/symmOrderAngle ) *symmOrderAngle;
+ float rotPhi0 = phi0- trunc(phi0/symmOrderAngle) *symmOrderAngle;
+
+
+// std::cout << "phi " << phi << " rotphi " << rotPhi << " phi0 " << phi0 << " rotphi0 " << rotPhi0 << " symordang " << symmOrderAngle;
+
+ rotPhi -= rotPhi0;
+// std::cout << " rotphi- " << rotPhi << std::endl;
+
+ if (isClosestDistance)
+ correctedDistance = distanceToIp/cos(rotPhi);
+ else
+ correctedDistance = cos(symmOrderAngle/2.0)*distanceToIp/cos(rotPhi);
+
+// std::cout << "distance " << distanceToIp << " corrected distance " << correctedDistance << std::endl;
+ }
+}
+
slicPandora/src
diff -u -r1.18 -r1.19
--- DetectorGeometry.cpp 5 Nov 2010 14:16:06 -0000 1.18
+++ DetectorGeometry.cpp 26 Nov 2010 14:52:27 -0000 1.19
@@ -148,6 +148,8 @@
layerExtra.m_emSamplingFraction = emSamplingFrac;
layerExtra.m_hadSamplingFraction = hadSamplingFrac;
layerExtra.m_cellThickness = cellThickness;
+ layerExtra.m_intLength = intLen;
+ layerExtra.m_radLength = radLen;
// Register extra layer parameters.
extras.m_extraLayerParams.push_back(layerExtra);
@@ -251,7 +253,7 @@
// Tracking parameters.
TiXmlElement* tracking = root->FirstChildElement("tracking");
float tinnerR, touterR, tz;
- tinnerR = touterR = tz;
+ tinnerR = touterR = tz = 0.f;
tracking->QueryFloatAttribute("innerR", &tinnerR);
tracking->QueryFloatAttribute("outerR", &touterR);
tracking->QueryFloatAttribute("z", &tz);
@@ -270,23 +272,28 @@
// Coil and B-field.
TiXmlElement* coil = root->FirstChildElement("coil");
- float cinnerR, couterR, cz, bfield;
+ float cinnerR, couterR, cz, bfield, coilRadLen = 0.f, coilIntLen = 0.f;
cinnerR = couterR = cz = bfield = 0.;
coil->QueryFloatAttribute("innerR", &cinnerR);
coil->QueryFloatAttribute("outerR", &couterR);
coil->QueryFloatAttribute("z", &cz);
coil->QueryFloatAttribute("bfield", &bfield);
+ coil->QueryFloatAttribute("radLen", &coilRadLen);
+ coil->QueryFloatAttribute("intLen", &coilIntLen);
geom.m_coilInnerRadius = cinnerR;
geom.m_coilOuterRadius = couterR;
geom.m_coilZExtent = cz;
+ m_coilRadiationLength = coilRadLen;
+ m_coilInteractionLength = coilIntLen;
+
SimpleBFieldCalculator::m_innerBField = bfield;
// FIXME: Next four are just dummy parameters for the moment.
- geom.m_nRadLengthsInZGap = 0.;
- geom.m_nIntLengthsInZGap = 0.;
- geom.m_nRadLengthsInRadialGap = 0.;
- geom.m_nIntLengthsInRadialGap = 0.;
+// geom.m_nRadLengthsInZGap = 0.;
+// geom.m_nIntLengthsInZGap = 0.;
+// geom.m_nRadLengthsInRadialGap = 0.;
+// geom.m_nIntLengthsInRadialGap = 0.;
#ifdef DETECTOR_GEOMETRY_DEBUG
// Print coil and field.