StepWiseHelixGeneralBrokenLineInterfaceProcessor.cc

Go to the documentation of this file.

#ifdef USE_GBL

#include "StepWiseHelixGeneralBrokenLineInterfaceProcessor.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 "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"

//#include "TMath.h"

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

namespace marlintpc {
static bool compare(indexArcPair one, indexArcPair two) {
        return ((one.second) < (two.second));
}

using namespace gbl;

StepWiseHelixGeneralBrokenLineInterfaceProcessor aStepWiseHelixGeneralBrokenLineInterfaceProcessor;

StepWiseHelixGeneralBrokenLineInterfaceProcessor::StepWiseHelixGeneralBrokenLineInterfaceProcessor() :
                Processor("StepWiseHelixGeneralBrokenLineInterfaceProcessor") {
        // the processor description
        _description =
                        "StepWiseHelixGeneralBrokenLineInterfaceProcessor is the interface between MarlinTPC and the General Broken Lines track (re-)fitting package.";

        registerInputCollection(LCIO::TRACK, "InputSeedTracks",
                        "The name of the input collection of track candidates (default: TPCSeedTracks)", _inputCollectionName,
                        std::string("TPCSeedTracks"));

        registerOutputCollection(LCIO::TRACK, "OutputTracks",
                        "The name of the output collection with the fitted tracks(default: TPCTracks)", _outputTrackCollectionName,
                        std::string("TPCTracks"));

        registerProcessorParameter("WriteOutputToStorage", "If true the output collection is written to file", _outputIsPersistent, true);

        registerOptionalParameter("BFieldScaleFactor", "Scales magnetic field (map), use 1.0 (default) for field ON or 0.0 for field OFF",
                        _bfieldScaleFactor, double(1.0));

        registerOptionalParameter("FitOptions",
                        "Set wether multiple passes with h huber, t tukey or c cauchy downweighting should be performed", _fitOptions,
                        std::string(""));

        registerOptionalParameter("WriteMillepedeBinary", "Set to true if you want a Millepede binary file (default:false)",
                        _writeMillepedeOut, false);

        registerOptionalParameter("MillePedeCalcMethod",
                        "Defines which parameters should be calculates for Millepede: 0: module alignment, 1: track distortions, 2: hit charge dependence (has yet to be implemented) (default:0)",
                        _milleCalcMethod, 0);

        registerOptionalParameter("MillePedeMinHits", "Minimum number of hits on track for output to Millepede binary file (default: 42)",
                        _milleMinHits, 42);

        registerOptionalParameter("MillepedeFilename", "Set the name of the Millepede binary file ", _fileNameMillepedeFile,
                        std::string(""));
        registerOptionalParameter("EncodedModuleID", "Module ID encoded in CellID0", _encodedModuleID, bool(true));

        registerOptionalParameter("MaxStep", "Maximal step size (with constant magnetic field)", _maxStep, double(10.0));

        registerProcessorParameter("X0PerUnitLength", "X0 per unit length (TPC gas as homogeneous scatterer)", _x0PerUnitLength,
                        double(0.0));

        registerProcessorParameter("ThickScatterer", "Selects thick (otherwise thin) scatterers for multiple scattering", _thickScatterer,
                        bool(false));

        registerProcessorParameter("DefaultMomentum", "Default momentum for multiple scattering (for Bfield off)", _defaultMomentum,
                        double(10.0));
}

void StepWiseHelixGeneralBrokenLineInterfaceProcessor::init() {
        printParameters();
        if (_writeMillepedeOut) {
                _milleBinary = new MilleBinary(_fileNameMillepedeFile.c_str());
        }

        const gear::TPCParameters& tpcParameters = Global::GEAR->getTPCParameters();
        // TPC center
        _Xcenter = tpcParameters.getDoubleVal("TPC_cylinder_centre_c0");
        _Ycenter = tpcParameters.getDoubleVal("TPC_cylinder_centre_c1");
}

void StepWiseHelixGeneralBrokenLineInterfaceProcessor::processRunHeader(LCRunHeader* run) {
}

void StepWiseHelixGeneralBrokenLineInterfaceProcessor::processEvent(LCEvent* evt) {
        // use an LCCollectionVec instead of LCCollection so we can use iterators later
        LCCollectionVec* seedTrackCollectionVector = NULL;

        // do a try / catch in case there is an empty event
        try {
                seedTrackCollectionVector = dynamic_cast<LCCollectionVec *>(evt->getCollection(_inputCollectionName));
        } catch (DataNotAvailableException &e) {
                // no input data, just return and process the next event
                return;
        }

        streamlog_out(DEBUG) << "in event " << evt->getEventNumber() << " there are " << seedTrackCollectionVector->size() << " tracks. "
                        << endl;

        Vector3D theTPCCenter(_Xcenter, _Ycenter, 0.);
        Vector3D bfield = marlin::Global::GEAR->getBField().at(theTPCCenter);

        const bool bOff(bfield.z() * _bfieldScaleFactor == 0.);
        if (bOff) {
                streamlog_out(DEBUG) << "BField(GEAR) is OFF" << std::endl;
        } else {
                streamlog_out(DEBUG) << "BField(GEAR) is ON" << std::endl;
        }
        // magnetic field (in z direction) at TPC center
        double Bz = bfield.z() * _bfieldScaleFactor; // rho component is accessible with field.rho()

        // output collection for the fitted tracks
        IMPL::LCCollectionVec* outputTrackCollection = new IMPL::LCCollectionVec(EVENT::LCIO::TRACK);

        // hits are to be stored with the track -- set the according flag
        LCFlagImpl trkFlag(0);
        trkFlag.setBit(LCIO::TRBIT_HITS);
        outputTrackCollection->setFlag(trkFlag.getFlag());

        // loop over the seed tracks
        for (LCCollectionVec::iterator inputIter = seedTrackCollectionVector->begin(), endIter = seedTrackCollectionVector->end();
                        inputIter < endIter; ++inputIter) {
                // dynamic cast to lcio::track
                Track* seedTrack = dynamic_cast<Track *>(*inputIter);

                double seedD0 = seedTrack->getD0();
                double seedPhi = seedTrack->getPhi();
                double seedOmega = seedTrack->getOmega();
                double seedZ0 = seedTrack->getZ0();
                double seedTanLambda = seedTrack->getTanLambda();
                vector<TrackerHit*> hitList = seedTrack->getTrackerHits();
                //CHK track parameters are defined at (distance of closest approach to) reference point !?
                const float* refPoint = seedTrack->getReferencePoint();
                const double refPos[] = { refPoint[0], refPoint[1], refPoint[2] };
                //CHK arc-length is measured from reference point
                float rInner = seedTrack->getRadiusOfInnermostHit();

                streamlog_out(DEBUG) << " ++ track number " << endIter - inputIter << " with " << seedTrack->getTrackerHits().size()
                                << " hits on the track has the parameters: [omega, tanLambda, phi0, d0, z0] = [" << seedOmega << "," << seedTanLambda << ","
                                << seedPhi << "," << seedD0 << "," << seedZ0 << "]." << endl;

                if (bOff) {
                        _curvature = false;
                        // current Clupatra fits curved tracks independent of magnetic field
                        // calculate  new seed from straight line defined by (radially) first and last hit
                        if (seedOmega != 0.) calcLineSeed(hitList, refPos, seedOmega, seedPhi, seedD0, seedTanLambda, seedZ0);
                } else {
                        _curvature = true;
                }

                const double seedLambda = atan(seedTanLambda);

                //CHK reference helix
                double seedPar[] = { -seedOmega, seedPhi, -seedD0, seedTanLambda, seedZ0 };
                // use local helix
                localHelix lHelix(seedPar, refPos, _bfieldScaleFactor);
                simpleHelix hlx = lHelix.getSimpleHelix();
                // lHelix.dump();

                vector<indexArcPair> hitWithArcList;
                for (unsigned int i = 0; i < hitList.size(); ++i) {
                        double sArc = hlx.getArcLengthXY(hitList[i]->getPosition());
                        hitWithArcList.push_back(make_pair(i, sArc));
                }
                // find smallest arc-length
                indexArcPair closest(hitWithArcList.front()), first(closest), last(closest);
                for (vector<indexArcPair>::const_iterator thePair = hitWithArcList.begin() + 1; thePair < hitWithArcList.end(); ++thePair) {
                        if (abs(thePair->second) < abs(closest.second)) closest = *thePair;
                        if (thePair->second < first.second) first = *thePair;
                        if (thePair->second > last.second) last = *thePair;
                }
                // add reference point if not at a hit (or close by)
                const double averageStep = (last.second - first.second) / (hitWithArcList.size() - 1);
                if (abs(closest.second / averageStep) > 1.0E-5) {
                        closest = make_pair(-1, 0.);
                        hitWithArcList.push_back(closest);
                }
                sort(hitWithArcList.begin(), hitWithArcList.end(), compare);

                // list of scatterers between measurements (relative arc-length and variance)
                vector<arcVarPair> arcWithVarList;
                if (_x0PerUnitLength > 0.) {
                        double ptot = _defaultMomentum;
                        if (!bOff) ptot = Bz * 0.0002998 / (seedOmega * cos(seedLambda));
                        double XbyX0 = _x0PerUnitLength * (last.second - first.second) / cos(seedLambda); // total material on track
                        // Highland formula, variance per unit (arc-)length in XY
                        double MultScatVar = pow(0.0136 / ptot * (1. + 0.038 * std::log(XbyX0)), 2) * _x0PerUnitLength / cos(seedLambda);
                        if (_thickScatterer) { // use two equivalent thin scatterers
                                arcWithVarList.push_back(make_pair(0.21132, 0.5 * MultScatVar)); // at 0.5-1/sqrt(12.)
                                arcWithVarList.push_back(make_pair(0.78868, 0.5 * MultScatVar)); // at 0.5+1/sqrt(12.)
                        } else
                                arcWithVarList.push_back(make_pair(0.5, MultScatVar)); // at 0.5

                }

                vector<GblPoint> theListOfPoints;
                // number of valid hits (with measurement)
                unsigned int numValid = 0;
                // point 2 point jacobian in local curvilinear parameters
                TMatrixD point2pointJacobian(5, 5);
                point2pointJacobian.UnitMatrix();
                bool makeScatterers = false; // in active volume, between first and last hit?
                // Energy loss estimators.
                vector<double> eLossEstimators;
                // number of hits lost (due to threshold)
                unsigned int numLow = 0;
                for (vector<indexArcPair>::const_iterator thePair = hitWithArcList.begin(); thePair < hitWithArcList.end(); ++thePair) {
                        // cout << " sArc " << thePair->first << ", " << thePair->second << endl;
                        int index = thePair->first;
                        double phiTrack = seedPhi - seedOmega * thePair->second;
                        // propagate to position
                        const double* position = index >= 0 ? hitList[index]->getPosition() : refPos;

                        // insert scatterers
                        if (makeScatterers) {
                                // simple helix (at last point)
                                simpleHelix hlx = lHelix.getSimpleHelix();
                                // distance (to this point)
                                double step = hlx.getArcLengthXY(position);
                                double scatPos[3];
                                for (vector<arcVarPair>::const_iterator theScat = arcWithVarList.begin(); theScat < arcWithVarList.end(); ++theScat) {
                                        // get position of scatterer
                                        hlx.getPosAtArcLength(theScat->first * step, scatPos);
                                        point2pointJacobian = lHelix.propagateTo(scatPos, _maxStep) * point2pointJacobian;
                                        // create a point from the jacobian
                                        GblPoint point(point2pointJacobian);
                                        // reset (integrated) jacobian
                                        point2pointJacobian.UnitMatrix();
                                        // add scatterer
                                        TVectorD scat(2);
                                        scat.Zero();
                                        TVectorD scatPrec(2);
                                        scatPrec[0] = scatPrec[1] = 1.0 / (theScat->second * step);
                                        point.addScatterer(scat, scatPrec);
                                        // store the point in the list that will be handed to the trajectory
                                        theListOfPoints.push_back(point);
                                }
                        }
                        if (_x0PerUnitLength > 0.) {
                                if (index == first.first) makeScatterers = true; // reached first hit
                                if (index == last.first) makeScatterers = false; // reached last hit
                        }

                        point2pointJacobian = lHelix.propagateTo(position, _maxStep) * point2pointJacobian;

                        gblHelperHit aHit =
                                        index >= 0 ?
                                                        gblHelperHit(*hitList[index], lHelix.getSimpleHelix(), _writeMillepedeOut, _encodedModuleID, _milleCalcMethod,
                                                                        phiTrack) :
                                                        gblHelperHit(seedPhi);
                        if (aHit.isValid()) {
                                numValid++;
                                // create a point from the jacobian
                                GblPoint point(point2pointJacobian);
                                // reset (integrated) jacobian
                                point2pointJacobian.UnitMatrix();
                                // now add the measurement point
                                if (aHit.isMeasurement()) {
                                        point.addMeasurement(aHit.getLocalToMeasurementProjection(), aHit.getResiduals(), aHit.getPrecision());
                                        if (aHit.getNumGlobals() > 0) point.addGlobals(aHit.getGlobalLabels(), aHit.getGlobalDerivatives());
                                        if (aHit.getEScaled() > 0.)
                                                eLossEstimators.push_back(1. / sqrt(aHit.getEScaled()));
                                        else
                                                numLow++;
                                }
                                // store the point in the list that will be handed to the trajectory
                                theListOfPoints.push_back(point);
                                // is reference point?
                                if (index == closest.first) _refPointIndex = theListOfPoints.size();
                        } else
                                streamlog_out(DEBUG) << " invalid hit: " << thePair->first << ", " << thePair->second << endl;
                }
                if (numValid < 3) {
                        streamlog_out(DEBUG) << " ++ track number " << endIter - inputIter << " skipped - too few valid hits " << endl;
                        continue;
                }
                _trajectory = new GblTrajectory(theListOfPoints, _curvature);

                if (!_trajectory->fit(_chisquare, _ndf, _lostweight, _fitOptions)) {
                        if (_writeMillepedeOut) {
                                trackerAlcaSelector sel(*seedTrack, _milleMinHits);
                                if (sel.isSelected()) _trajectory->milleOut(*_milleBinary);
                        }
                        TVectorD trackparameters(5);
                        TMatrixDSym trackcovariance(5);
                        if (_ndf > 0) // check if the fit was successful
                                        {
                                // get the results at a given label in local cl track parameters
                                // the track parameters are corrections to the curvilinear track parameters
                                _trajectory->getResults(_refPointIndex, trackparameters, trackcovariance);
                                // the results need to be transformed to the LCIO perigee representation ( d0, phi0, Omega, z0, tanLambda)
                                TMatrixD jacCL2LCIO = hlx.perigeeToLCIOJacobian() * hlx.curvilinearToPerigeeJacobian();
                                TVectorD correctionVec = jacCL2LCIO * trackparameters;
                                TMatrixDSym covarianceMatrix = trackcovariance.Similarity(jacCL2LCIO);
                                /* prepare LCIO output; unknown parts first:
                                 *
                                 * theTrack->setTypeBit (int index, bool val=true) // ?
                                 *
                                 */

                                TrackImpl* theTrack = new TrackImpl();
                                theTrack->setOmega(seedOmega + correctionVec[2]);
                                theTrack->setTanLambda(seedTanLambda + correctionVec[4]);
                                theTrack->setPhi(seedPhi + correctionVec[1]);
                                theTrack->setD0(seedD0 + correctionVec[0]);
                                theTrack->setZ0(seedZ0 + correctionVec[3]);
                                theTrack->setReferencePoint(refPoint);
                                theTrack->setRadiusOfInnermostHit(rInner);
                                // compressed covariance matrix
                                float compressedCovariance[15];
                                int ind = 0;
                                for (int i = 0; i < 5; ++i)
                                        for (int j = 0; j <= i; ++j)
                                                compressedCovariance[ind++] = covarianceMatrix[i][j];
                                theTrack->setCovMatrix(compressedCovariance);

                                theTrack->setChi2(_chisquare);
                                theTrack->setNdf(_ndf);

                                // simple dEdx from median(1/sqrt(EScaled))
                                const unsigned int numEstimators = eLossEstimators.size();
                                if (numEstimators > numLow) {
                                        sort(eLossEstimators.begin(), eLossEstimators.end());
                                        const int nTot = numEstimators + numLow;
                                        const int indMed = nTot / 2;
                                        const double eMed =
                                                        (nTot % 2 == 0) ? 0.5 * (eLossEstimators[indMed - 1] + eLossEstimators[indMed]) : eLossEstimators[indMed]; // median
                                        vector<double> absDev;
                                        for (unsigned int i = 0; i < numEstimators; ++i)
                                                absDev.push_back(abs(eLossEstimators[i] - eMed));
                                        sort(absDev.begin(), absDev.end());
                                        const double eMAD = (nTot % 2 == 0) ? 0.5 * (absDev[indMed - 1] + absDev[indMed]) : absDev[indMed]; // median of absolute deviations, sigma(gauss) = 1.4826 * MAD
                                        const float ndEdx = nTot;
                                        const float dEdx = 1.0 / (eMed * eMed);
                                        const float dEdxError = 2. * 1.4826 * eMAD * dEdx / eMed / sqrt(ndEdx);
                                        streamlog_out(DEBUG) << " dEdx " << ndEdx << " " << dEdx << " " << dEdxError << endl;
                                        theTrack->setdEdx(dEdx);
                                        theTrack->setdEdxError(dEdxError);
                                }
                                // and finally add the hits to the track
                                for (vector<TrackerHit*>::const_iterator theHit = seedTrack->getTrackerHits().begin(), hitEnd =
                                                seedTrack->getTrackerHits().end(); theHit < hitEnd; ++theHit)
                                        theTrack->addHit(*theHit);

                                streamlog_out(DEBUG) << " +++ the fitted track has a chi^2 value of " << _chisquare << " for ndf = " << _ndf
                                                << " and the parameters [omega, tanLambda, phi0, d0, z0] = [" << theTrack->getOmega() << "," << theTrack->getTanLambda()
                                                << "," << theTrack->getPhi() << "," << theTrack->getD0() << "," << theTrack->getZ0() << "]." << endl;
                                outputTrackCollection->addElement(theTrack);
                        }
                }
                // cleanup
                if (_trajectory) delete _trajectory;
        }
        // final storage to the event
        if (_outputIsPersistent) evt->addCollection(outputTrackCollection, _outputTrackCollectionName);
}

void StepWiseHelixGeneralBrokenLineInterfaceProcessor::check(LCEvent* evt) {
}

void StepWiseHelixGeneralBrokenLineInterfaceProcessor::end() {
        // close Millepede-II binary file
        if (_milleBinary) delete _milleBinary;
}


void StepWiseHelixGeneralBrokenLineInterfaceProcessor::calcLineSeed(const std::vector<TrackerHit*> hitList, const double* refPos,
                double& omega, double& phi, double& d0, double& tanLambda, double& z0) const {
        //
        TMatrixD xyzr(2, 4);
        for (vector<TrackerHit*>::const_iterator theHit = hitList.begin(), hitEnd = hitList.end(); theHit < hitEnd; ++theHit) {
                const EVENT::TrackerHit& aHit = **theHit;
                const double x = aHit.getPosition()[0];
                const double y = aHit.getPosition()[1];
                const double z = aHit.getPosition()[2];
                const double r = sqrt(pow(x, 2) + pow(y, 2));
                if (xyzr[0][3] == 0.)                                   // for first hit
                                {
                        xyzr[0][0] = x;
                        xyzr[0][1] = y;
                        xyzr[0][2] = z;
                        xyzr[0][3] = r;
                        xyzr[1][0] = x;
                        xyzr[1][1] = y;
                        xyzr[1][2] = z;
                        xyzr[1][3] = r;
                }
                if (r < xyzr[0][3]) // new innermost hit
                                {
                        xyzr[0][0] = x;
                        xyzr[0][1] = y;
                        xyzr[0][2] = z;
                        xyzr[0][3] = r;
                }
                if (r > xyzr[1][3]) // new outermost hit
                                {
                        xyzr[1][0] = x;
                        xyzr[1][1] = y;
                        xyzr[1][2] = z;
                        xyzr[1][3] = r;
                }
        }
        const double dx = xyzr[1][0] - xyzr[0][0];
        const double dy = xyzr[1][1] - xyzr[0][1];
        const double dz = xyzr[1][2] - xyzr[0][2];
        const double dr = sqrt(pow(dx, 2) + pow(dy, 2));
        const double x1 = xyzr[0][0] - refPos[0];
        const double y1 = xyzr[0][1] - refPos[1];
        const double z1 = xyzr[0][2] - refPos[2];
        // new track parameters
        omega = 0.;
        phi = atan2(dy, dx);
        d0 = y1 * cos(phi) - x1 * sin(phi);
        tanLambda = dz / dr;
        z0 = z1 - tanLambda * (x1 * cos(phi) + y1 * sin(phi));
        //cout << " new seed " << phi << ", " << d0 << ", " << tanLambda << ", " << z0 << endl;
}

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