gbl_example.cxx

Go to the documentation of this file.

/*
 * this is an example of the usage
 * it's based on the example1.cpp in GBL by Claus Kleinwort
 */
 
#define AIDATT_USE_DD4HEP 1
#define USE_LCIO 1
#define USE_GBL 1
 
// c++
#include <iostream>
#include <map>
#include <string>
 
// lcio
#include "lcio.h"
#include "IO/LCWriter.h"
#include "EVENT/LCEvent.h"
#include "EVENT/LCCollection.h"
#include "EVENT/SimTrackerHit.h"
#include "UTIL/LCTrackerConf.h"
 
#include <IMPL/LCCollectionVec.h>
#include "IMPL/TrackImpl.h"
 
// DD4hep
#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
#include "DDRec/SurfaceHelper.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DDRec/CellIDPositionConverter.h"
 
 
 
// aidaTT
#include "aidaTT/DD4hepGeometry.hh"
#include "aidaTT/AidaTT.hh"
#include "aidaTT/ConstantSolenoidBField.hh"
#include "aidaTT/analyticalPropagation.hh"
#include "aidaTT/simplifiedPropagation.hh"
#include "aidaTT/GBLInterface.hh"
#include "aidaTT/fitResults.hh"
 
using namespace std;
using namespace lcio;
 
 
int main(int argc, char** argv)
{
    if(argc < 2)
        {
            std::cout << " usage: ./gbl_example ILDEx.xml ILDExSimu.slcio" << std::endl ;
            return 1;
        }
 
    std::string inFile =  argv[1] ;
 
    dd4hep::Detector& theDet = dd4hep::Detector::getInstance();
    theDet.fromCompact(inFile);
 
    dd4hep::DetElement world = theDet.world() ;
 
    aidaTT::DD4hepGeometry geom(theDet);
 
    const auto& surfaces = geom.getSurfaces() ;
 
    dd4hep::rec::CellIDPositionConverter cellid_converter(theDet);
 
    dd4hep::rec::SurfaceManager& surfMan = *theDet.extension<dd4hep::rec::SurfaceManager>() ;
    //auto surfMap = surfMan.map( "world" ) ;
 
    // create map of surfaces
    std::map< long64, const aidaTT::ISurface* > surfMap ;
    for(const auto& surf : surfaces)
    {
      surfMap[surf->id() ] = (surf) ;
      //surfMap[surf->id() ]->setProperty(dd4hep::rec::SurfaceType::OrthogonalToZ,true);
    }
 
 
 
    std::string lcioFileName = argv[2] ;
 
    LCReader* rdr = LCFactory::getInstance()->createLCReader() ;
    rdr->open(lcioFileName) ;
 
    std::string outFileName = "trackTest.slcio";
 
    LCWriter* wrt = LCFactory::getInstance()->createLCWriter() ;
    wrt->open(outFileName) ;
 
    LCEvent* evt = 0 ;
 
    std::vector< std::string > colNames ;
    colNames.push_back("SiTrackerBarrelHits") ;
    //colNames.push_back("SiVertexBarrelHits") ;
    //colNames.push_back("SiTrackerEndcapHits") ;
    //colNames.push_back("SiVertexEndcapHits") ;
 
    UTIL::BitField64 idDecoder(LCTrackerCellID::encoding_string()) ;
 
 
    // the aidaTT stuff
 
    // with B = 3.5T => omega = 3* 10^(-4) * 3.5 / sqrt(6)
    aidaTT::trackParameters bogusTP;
 
    //~ const double initialOmega = 3.5*1e-5;
    //~ const double initialTanLambda = 0.;
    //~ const double initialPhi = 5*M_PI_4;
 
    const double initialOmega = 1.5 * 1e-4;
    const double initialTanLambda = 1.;
    const double initialPhi = M_PI_2;
 
    bogusTP.setTrackParameters(aidaTT::Vector5(initialOmega, initialTanLambda, initialPhi, 0., 0.));
 
    // create the different objects needed for fitting
    // first a constant field parallel to z, 1T
    aidaTT::ConstantSolenoidBField*  bfield = new aidaTT::ConstantSolenoidBField(3.5);
 
    // create the propagation object
    aidaTT::analyticalPropagation* propagation = new aidaTT::analyticalPropagation();
    //aidaTT::simplifiedPropagation* propagation = new aidaTT::simplifiedPropagation();
 
    // create the fitter object
    aidaTT::GBLInterface* fitter = new aidaTT::GBLInterface();
 
 
    while((evt = rdr->readNextEvent()) != 0)
        {
 
            IMPL::LCCollectionVec* outputTrackCollection = new IMPL::LCCollectionVec(EVENT::LCIO::TRACK);
            std::string outName = "trytracks";
 
            std::cout << " reading an event " << std::endl;
            // now create a trajectory object to be fitted
            aidaTT::trajectory theMaster(bogusTP, fitter,  propagation, &geom);
 
            for(unsigned icol = 0, ncol = colNames.size() ; icol < ncol ; ++icol)
                {
                    LCCollection* col = NULL;
                    try
                    {
                      col = evt->getCollection(colNames[ icol ]) ;
                    }
                    catch(DataNotAvailableException &e)
                    {
                      continue;
                    }
 
                    int nHit = col->getNumberOfElements() ;
 
                    for(int i = 0 ; i < nHit ; ++i)
                        {
                            SimTrackerHit* sHit = (SimTrackerHit*) col->getElementAt(i) ;
                            long64 id = sHit->getCellID0() ;
                            //idDecoder.setValue(id) ;
 
                            //const auto si = surfMap->find( cellid_converter.findContext(id)->identifier );
                            //    // Get the Surface (http://test-dd4hep.web.cern.ch/test-dd4hep/doxygen/html/classdd4hep_1_1rec_1_1_i_surface.html)
                            //std::cout << surfMap->count( cellid_converter.findContext(id)->identifier ) << "\n";
                            //dd4hep::rec::ISurface* surf = (si != surfMap->end() ?  si->second  : nullptr);
 
                            const aidaTT::ISurface* surf = surfMap[ id ] ;
 
                            dd4hep::rec::Vector3D  pos1;
                            pos1 = surf->origin();
                            double recalcPos[3] = {pos1.x(), pos1.y(), pos1.z()}; // could it be more ugly?
                            if(surf->type().isSensitive()) {
 
                              std::vector<double> precisionDummy;
                              precisionDummy.push_back(1. / 0.0001);
                              precisionDummy.push_back(1. / 0.00012);
                              //for(unsigned int i = 0; i < 3; ++i)
                              //  recalcPos[i] = sHit->getPosition()[i] * dd4hep::mm;
 
                              theMaster.addMeasurement(recalcPos, precisionDummy, *surf, sHit);
                            }
                        }
 
                    theMaster.prepareForFitting();
                    theMaster.fit();
 
                    const aidaTT::fitResults& result = *(theMaster.getFitResults());
 
                    TrackImpl* theTrack = new TrackImpl();
                    theTrack->setOmega(result.estimatedParameters()(0)) ;
                    theTrack->setTanLambda(result.estimatedParameters()(1));
                    theTrack->setPhi(result.estimatedParameters()(2));
                    theTrack->setD0(result.estimatedParameters()(3));
                    theTrack->setZ0(result.estimatedParameters()(4));
 
                    TrackImpl* theInitialTrack = new TrackImpl();
                    theInitialTrack->setOmega(initialOmega) ;
                    theInitialTrack->setTanLambda(initialTanLambda);
                    theInitialTrack->setPhi(initialPhi);
                    theInitialTrack->setD0(0.);
                    theInitialTrack->setZ0(0.);
 
 
                    outputTrackCollection->addElement(theTrack);
                    outputTrackCollection->addElement(theInitialTrack);
 
                }
            evt->addCollection(outputTrackCollection , outName);
 
            wrt->writeEvent(evt);
 
        }
 
 
    return 0;
}