GeneralBrokenLineInterfaceHelpers.cc

Go to the documentation of this file.

#ifdef USE_GBL

#include "GeneralBrokenLineInterfaceHelpers.h"

#include <iostream>
#include <string>

// include what you need from lcio
#include <EVENT/LCEvent.h>
#include <IMPL/LCCollectionVec.h>
#include <EVENT/Track.h>
#include <EVENT/TrackerHit.h>
#include <EVENT/TrackerPulse.h>

#include "IMPL/TrackImpl.h"
#include "IMPL/LCFlagImpl.h"
#include <Exceptions.h>
#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>

#include "marlin/Global.h"

// GEAR
#include "GEAR.h"
#include <gear/TPCParameters.h>
#include <gear/TPCModule.h>
//#include <gear/BField.h>
#include <gearimpl/Vector3D.h>

// the gbl includes
#include "GblPoint.h"

// cr: needed for accessing the magnetic field
#include "GlobalFieldProcessor.h"
#include "CLHEP/Vector/ThreeVector.h"

//#include "TMath.h"

using namespace lcio;
using namespace marlin;
using namespace std;
using namespace gear;

namespace marlintpc {


gblHelperHit::gblHelperHit(const EVENT::TrackerHit& aHit, const simpleHelix& refHelix, const bool milleOut, const bool enMID,
                const int milleCalcMethod, const double phiSeed) :
                _isMeas(true), _hitPhi(_calcLocalPhi(aHit, enMID, phiSeed)), _isValid(
                                refHelix.getExpectedPlanePos(aHit.getPosition(), _hitPhi + 0.5 * M_PI, _xyPred, _zPred, _s, _seedPhi, _seedLambda)), _eScaled(
                                _calcEScaled(aHit)), _l2M(_getL2M()), _res(_getRes()), _prec(_getPrec(aHit.getCovMatrix())), _numGlobals(0), _globalLab(), _globalDer() {
        // needs only to be done if _writeMillepedeOut is true
        if (milleOut) {
                TVectorD dir(3);
                dir[0] = cos(_seedPhi) * cos(_seedLambda);
                dir[1] = sin(_seedPhi) * cos(_seedLambda);
                dir[2] = sin(_seedLambda);
                // calculate global derivatives
//              std::cout<< "_seedPhi=" << _seedPhi <<std::endl;
//              std::cout<< "_hitPhi=" << _hitPhi <<std::endl;
//              std::cout<< "_seedPhi-_hitPhi=" << _seedPhi-_hitPhi <<std::endl;

                trackerAlcaCalculator calc;
                switch (milleCalcMethod) {
                case 0:
                        calc.calcModuleAlignment(aHit, dir);
                        break;
                case 1:
                        calc.calcTrackDistortions(aHit);
                        break;
                case 2:
                        calc.calcHitCharge(aHit);
                        break;
                case 3:
                        calc.calcDr(aHit, _seedPhi - _hitPhi);
                        break;
                default:
                        calc.calcModuleAlignment(aHit, dir);
                        break;
                }
                _numGlobals = calc.getNumGlobals();
                _globalLab = calc.getGlobalLabels();
                _globalDer.ResizeTo(2, _numGlobals);
                _globalDer = calc.getGlobalDerivatives();
        }
}


gblHelperHit::gblHelperHit(const double phi) :
                _isMeas(false) {
        _s = 0.;
        _seedPhi = phi;
        _isValid = true;
}


TMatrixD gblHelperHit::_getL2M() {
        TMatrixD projectLocalToMeasurement(2, 2);
        // Inverse (Pu = du/dm) is easy
        projectLocalToMeasurement[0][0] = cos(_hitPhi - _seedPhi);
        projectLocalToMeasurement[0][1] = 0.;
        projectLocalToMeasurement[1][0] = sin(_hitPhi - _seedPhi) * sin(_seedLambda);
        projectLocalToMeasurement[1][1] = cos(_seedLambda);
        // Invert
        projectLocalToMeasurement.Invert();
        //cout << " proL2M " << endl;
        //projectLocalToMeasurement.Print();

        return projectLocalToMeasurement;
}


TMatrixDSym gblHelperHit::_getPrec(FloatVec covMatrix) {
        TMatrixDSym precision(2);
        const double cosPhi = cos(_hitPhi);
        const double sinPhi = sin(_hitPhi);
        const double varRad = cosPhi * cosPhi * covMatrix[0] + 2. * sinPhi * cosPhi * covMatrix[1] + sinPhi * sinPhi * covMatrix[2];
        const double varPhi = sinPhi * sinPhi * covMatrix[0] - 2. * sinPhi * cosPhi * covMatrix[1] + cosPhi * cosPhi * covMatrix[2];
        const double varZ = covMatrix[5];
        const double derXY = tan(_seedPhi - _hitPhi);
        const double derZS = tan(_seedLambda) / cos(_seedPhi - _hitPhi);
        // covariance matrix
        precision[0][0] = varPhi + varRad * derXY * derXY;
        precision[0][1] = (varRad > 1.0E-6 * varPhi) ? varRad * derXY * derZS : 0.; // avoid correlations only from float precision
        precision[1][0] = precision[0][1];
        precision[1][1] = varZ + varRad * derZS * derZS;
        // needs to be inverted
        const double det = precision[0][0] * precision[1][1] - precision[0][1] * precision[1][0];
        if (det)
                precision.InvertFast();
        else {
                precision.Zero();
                streamlog_out(ERROR) << "Singular covariance matrix" << std::endl;
        }
        return precision;
}

TVectorD gblHelperHit::_getRes() {
        TVectorD residuals(2);
        // predictions are relative to the measurements
        residuals[0] = -_xyPred;
        residuals[1] = -_zPred;
        return residuals;
}


double gblHelperHit::_calcLocalPhi(const EVENT::TrackerHit& aHit, const bool encodedModuleID, const double phiTrack) {
        FloatVec covMatrix = aHit.getCovMatrix();
        double phi; // radial direction

        // analyze (XY) covariance matrix
        // std::cout << " cov-xy " << aHit.getCellID0() << ": " << covMatrix[0] << " " << covMatrix[1] << " " << covMatrix[2] << " " << covMatrix[5] << std::endl;
        const double trace = covMatrix[0] + covMatrix[2];
        const double det = covMatrix[0] * covMatrix[2] - covMatrix[1] * covMatrix[1];
        // std::cout << " T,D: " << trace << " " << det << " " << 4.0*det/(trace*trace) << std::endl;
        if (4.0 * det / (trace * trace) < 0.9999) {
                // eigenvalues are significantly different
                // hit direction from covariance matrix (perpendicular to measurement direction)
                double phiCov;
                if (4.0 * det / (trace * trace) < 0.000001) {
                        // rank of XY covariance matrix is 1, hit direction has eigenvalue/variance 0.
                        phiCov = 0.5 * atan2(-2.0 * covMatrix[1], covMatrix[2] - covMatrix[0]);
                } else {
                        // rank of XY covariance matrix is 2, hit direction has larger eigenvalue/variance (than measurement direction)
                        phiCov = 0.5 * atan2(2.0 * covMatrix[1], covMatrix[0] - covMatrix[2]);
                }
                // hit direction from GEAR
                UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
                encoder.setValue(aHit.getCellID0());
                const gear::TPCModule &module =
                                encodedModuleID ?
                                                Global::GEAR->getTPCParameters().getModule(encoder[ILDCellID0::module]) :
                                                Global::GEAR->getTPCParameters().getNearestModule(aHit.getPosition()[0], aHit.getPosition()[1]);

                if (module.getLocalPadLayout().getCoordinateType() == gear::PadRowLayout2D::POLAR) {
                        Vector2D correction = module.getOffset();
                        phi = atan2(aHit.getPosition()[1] - correction[1], aHit.getPosition()[0] - correction[0]);
                        // std::cout << " Polar " << phi << std::endl;
                } else {
                        phi = module.getAngle() - 0.5 * M_PI;
                        // std::cout << " Cartesian " << phi << std::endl;
                }
                if (abs(phi - phiCov) > 0.00001)
                        streamlog_out(WARNING) << " Inconsistent hit direction, check hit errors, use TPCHit_FixErrors_forTracking " << std::endl;
        } else {
                // eigenvalues are similar
                // std::cout << " similar eigenvalues, phiTrack " << phiTrack << std::endl;
                phi = phiTrack; // measurement is perpendicular to track direction
        }
        return phi;
}


double gblHelperHit::_calcEScaled(const EVENT::TrackerHit& aHit) {
        return aHit.getEDep() * fabs(cos(_hitPhi - _seedPhi) * cos(_seedLambda));
}

unsigned int gblHelperHit::getNumGlobals() const {
        return _numGlobals;
}

std::vector<int> gblHelperHit::getGlobalLabels() const {
        return _globalLab;
}

const TMatrixD& gblHelperHit::getGlobalDerivatives() const {
        return _globalDer;
}


trackerAlcaCalculator::trackerAlcaCalculator() :
                _numDer(0), _globalLab(), _globalDer() {
//      switch (calcMethod) {
//      case 0:
//              trackerAlcaCalculator::calcModuleAlignment(aHit, tdir);
//              break;
//      case 1:
//              trackerAlcaCalculator::calcTrackDistortions(aHit);
//              break;
//      case 2:
//              trackerAlcaCalculator::calcHitCharge(aHit);
//              break;
//      default:
//              trackerAlcaCalculator::calcModuleAlignment(aHit, tdir);
//      }
}

unsigned int trackerAlcaCalculator::getNumGlobals() const {
        return _numDer;
}

std::vector<int> trackerAlcaCalculator::getGlobalLabels() const {
        return _globalLab;
}

const TMatrixD& trackerAlcaCalculator::getGlobalDerivatives() const {
        return _globalDer;
}


void trackerAlcaCalculator::calcTrackDistortions(const EVENT::TrackerHit& aHit) {

        // get module and row number
//      UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
//      encoder.setValue(aHit.getCellID0());
        const gear::TPCModule &theModule = Global::GEAR->getTPCParameters().getNearestModule(aHit.getPosition()[0],
                        aHit.getPosition()[1]);
        const int padIndex = theModule.getNearestPad(aHit.getPosition()[0], aHit.getPosition()[1]);
        const int nModules = Global::GEAR->getTPCParameters().getNModules();
        const int moduleNum = theModule.getModuleID();
//        const int nRows = theModule.getNRows();
//      const int layer = encoder[ILDCellID0::layer];
        const int rowNum = theModule.getRowNumber(padIndex);
        const int rowLabel = 100 * moduleNum + rowNum + 1;

        //
        // === row-wise offsets

        _numDer = 2;
        _globalLab.push_back(rowLabel);
        _globalLab.push_back(rowLabel + 100 * nModules);
        _globalDer.ResizeTo(2, _numDer);
        _globalDer[0][0] = 1.;          // offset in XY
        _globalDer[1][1] = 1.;          // offset in Z
}


void trackerAlcaCalculator::calcModuleAlignment(const TrackerHit & aHit, const TVectorD& tdir) {

        UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
        encoder.setValue(aHit.getCellID0());

        // ===  pad plane alignment  ===
        // --- hit position in pad plane system
        //int modNum = encoder[ILDCellID0::module];
        const gear::TPCModule &module = Global::GEAR->getTPCParameters().getNearestModule(aHit.getPosition()[0], aHit.getPosition()[1]);
        int modNum = module.getModuleID();
        Vector2D correction = module.getOffset();
        const double xHit = aHit.getPosition()[0] - correction[0];
        const double yHit = aHit.getPosition()[1] - correction[1];
        const double zHit = aHit.getPosition()[2];
        const double phiHit = atan2(yHit, xHit);
        // measurement direction (-sin(_phiHit), cos(_phiHit))
        const double ex = -sin(phiHit);
        const double ey = cos(phiHit);
        // --- normal vector (to (virtual) measurement plane)
        TVectorD normal(3);
        normal[0] = ey;
        normal[1] = -ex;
        //cout << " tdir " << endl; tdir.Print();
        //cout << " normal " << endl; _normal.Print();
        // --- dr/dm (residual vs measurement)
        TMatrixD drdm(3, 3);
        drdm.UnitMatrix();
        // product tdir * normal (scalar)
        double tn = tdir[0] * normal[0] + tdir[1] * normal[1] + tdir[2] * normal[2];
        // (subtract) product tdir * transposed(normal) (matrix)
        for (unsigned int i = 0; i < 3; i++)
                for (unsigned int j = 0; j < 3; j++)
                        drdm[i][j] -= tdir[i] * normal[j] / tn;
        //cout << " drdm " << endl; _drdm.Print();
        // --- dm/dg (measurement vs 6 rigid body parameters)
        TMatrixD dmdg(3, 6);
        dmdg[0][0] = 1.;
        dmdg[0][4] = zHit;
        dmdg[0][5] = -yHit;
        dmdg[1][1] = 1.;
        dmdg[1][3] = -zHit;
        dmdg[1][5] = xHit;
        dmdg[2][2] = 1.;
        dmdg[2][3] = yHit;
        dmdg[2][4] = -xHit;
        //cout << " dmdg " << endl; _dmdg.Print();
        // --- drl/drg (local vs global residuals)
        TMatrixD drldrg(3, 3);
        drldrg[0][0] = ex;
        drldrg[0][1] = ey; // in pad direction
        drldrg[1][2] = 1.;                        // in drift direction
        drldrg[2][0] = -ey;
        drldrg[2][1] = ex; // in row direction
        //cout << " drldrg " << endl; _drldrg.Print();
        // --- drl/dg (local residuals vs rigid body parameters)
        //TMatrixD _tmp = _drdm * _dmdg;
        //cout << " tmp " << endl; _tmp.Print();
        TMatrixD drldg = drldrg * drdm * dmdg;
        //cout << " drldg " << endl; _drldg.Print();
        _numDer = 6; // +2; // for distorions
        _globalDer.ResizeTo(2, _numDer);
        for (unsigned int i = 0; i < 6; i++) {
                _globalLab.push_back(10 * modNum + i + 1);
                _globalDer[0][i] = drldg[0][i];  // derivatives for pad measurement
                _globalDer[1][i] = drldg[1][i];  // derivatives for drift measurement
        }
        /* higher order distortions along rows
         const double row = rowNum;
         const double relRow = (row-14.5)/14.0;    // relative row number in [-1 .. +1]
         _globalLab.push_back(10*modNum + 7);       // distortion for pad measurement
         _globalDer[0][6] = pow(relRow,2) - 1./3.; // (rescaled) 2nd order Legendre polynomial
         _globalLab.push_back(10*modNum + 8);       // distortion for drift measurement
         _globalDer[1][7] = pow(relRow,2) - 1./3.; // (rescaled) 2nd order Legendre polynomial
         */
}


void trackerAlcaCalculator::calcHitCharge(const TrackerHit & aHit) {

        // get module and row number
//      UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
//      encoder.setValue(aHit.getCellID0());
        const gear::TPCModule &theModule = Global::GEAR->getTPCParameters().getNearestModule(aHit.getPosition()[0],
                        aHit.getPosition()[1]);
        const int padIndex = theModule.getNearestPad(aHit.getPosition()[0], aHit.getPosition()[1]);
        const int nModules = Global::GEAR->getTPCParameters().getNModules();
        const int moduleNum = theModule.getModuleID();
//        const int nRows = theModule.getNRows();
//      const int layer = encoder[ILDCellID0::layer];
        const int rowNum = theModule.getRowNumber(padIndex);
        const int rowLabel = 100 * moduleNum + rowNum + 1;

        //calc charge bin
        double hitMaxPulseCharge = 0;
        for (LCObjectVec::const_iterator pulseIter = aHit.getRawHits().begin(); pulseIter < aHit.getRawHits().end(); ++pulseIter) {
                TrackerPulse const * pulse = dynamic_cast<TrackerPulse const *>(*pulseIter);
                if (pulse->getCharge() > hitMaxPulseCharge) {
                        hitMaxPulseCharge = pulse->getCharge();
                }
        }
        const double lgQ = log10(hitMaxPulseCharge);
        int bin;
        if (lgQ >= 4)
                bin = 40;
        else if (lgQ < 2)
                bin = 1;
        else {
                bin = static_cast<int>((lgQ - 2) / 2 * 40) + 1;
        }
        if (bin < 1) bin = 1;

        _numDer = 4;
        _globalLab.push_back(rowLabel);
        _globalLab.push_back(rowLabel + 100 * nModules);
        _globalLab.push_back(bin + 2 * 100 * nModules);
        _globalLab.push_back(bin + 2 * 100 * nModules + 100);
        _globalDer.ResizeTo(2, _numDer);
        _globalDer[0][0] = 1.;  // offset in XY
        _globalDer[1][1] = 1.;  // offset in Z
        _globalDer[0][2] = 1.;  // offset in XY
        _globalDer[1][3] = 1.;  // offset in Z
}


void trackerAlcaCalculator::calcDr(const TrackerHit & aHit, double phi) {

        // row number
        UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
        encoder.setValue(aHit.getCellID0());
        int rowNum = encoder[ILDCellID0::layer];
        //calc charge bin
        _numDer = 2;
        _globalLab.push_back(rowNum + 1);
        _globalLab.push_back(rowNum + 1 + 100);
        _globalDer.ResizeTo(2, _numDer);
        _globalDer[0][0] = 1.;            // offset in XY
        _globalDer[0][1] = tan(phi);  // offset in XY
}


trackerAlcaSelector::trackerAlcaSelector(const EVENT::Track& aTrack, const int milleMinHits) :
                _selectionFlag(true) {
        // for LPTPC require minimum track length (~ 1.5 modules)
        unsigned int minHits = milleMinHits;
        _selectionFlag = (aTrack.getTrackerHits().size() >= minHits);
}

const bool trackerAlcaSelector::isSelected() const {
        return _selectionFlag;
}

} // close namespace brackets
#endif // conditional compilation if GBL program package is available on the system