diff --git a/Digitisers/SimpleDigi/CMakeLists.txt b/Digitisers/SimpleDigi/CMakeLists.txt
index c3edd2d784a9272328d4afb075bebfff05925ef5..39f60902538c88005ccca1d64f9e486e9db2c0f2 100644
--- a/Digitisers/SimpleDigi/CMakeLists.txt
+++ b/Digitisers/SimpleDigi/CMakeLists.txt
@@ -12,6 +12,7 @@ gaudi_depends_on_subdirs(
Service/GearSvc
Service/EventSeeder
Service/TrackSystemSvc
+ Utilities/DataHelper
)
set(SimpleDigi_srcs src/*.cpp)
@@ -19,5 +20,5 @@ set(SimpleDigi_srcs src/*.cpp)
# Modules
gaudi_add_module(SimpleDigi ${SimpleDigi_srcs}
INCLUDE_DIRS K4FWCore GaudiKernel GaudiAlgLib ${CLHEP_INCLUDE_DIR} gear ${GSL_INCLUDE_DIRS} ${LCIO_INCLUDE_DIRS}
- LINK_LIBRARIES K4FWCore GaudiKernel GaudiAlgLib ${CLHEP_LIBRARIES} ${GEAR_LIBRARIES} ${GSL_LIBRARIES} ${LCIO_LIBRARIES} EDM4HEP::edm4hep EDM4HEP::edm4hepDict
+ LINK_LIBRARIES K4FWCore GaudiKernel GaudiAlgLib ${CLHEP_LIBRARIES} ${GEAR_LIBRARIES} ${GSL_LIBRARIES} ${LCIO_LIBRARIES} EDM4HEP::edm4hep EDM4HEP::edm4hepDict DataHelperLib
)
diff --git a/Digitisers/SimpleDigi/src/TPCDigiAlg.cpp b/Digitisers/SimpleDigi/src/TPCDigiAlg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ebb1df487d017ec86f5216ec958e222114e5f505
--- /dev/null
+++ b/Digitisers/SimpleDigi/src/TPCDigiAlg.cpp
@@ -0,0 +1,1635 @@
+/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+
+#include "TPCDigiAlg.h"
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <map>
+#include <cmath>
+#include <algorithm>
+#include <array>
+
+#include <gsl/gsl_randist.h>
+
+#include "DataHelper/Circle.h"
+#include "DataHelper/SimpleHelix.h"
+#include "DataHelper/LCCylinder.h"
+
+#include "EventSeeder/IEventSeeder.h"
+//#include <IMPL/LCFlagImpl.h>
+//#include <IMPL/LCRelationImpl.h>
+//#include <IMPL/TrackerHitImpl.h>
+
+//stl exception handler
+#include <stdexcept>
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "voxel.h"
+
+#include "GearSvc/IGearSvc.h"
+// STUFF needed for GEAR
+#include <gear/GEAR.h>
+#include <gear/TPCParameters.h>
+#include <gear/PadRowLayout2D.h>
+#include <gear/BField.h>
+//
+#include "UTIL/ILDConf.h"
+
+#define TRKHITNCOVMATRIX 6
+
+//using namespace lcio ;
+using namespace std ;
+
+DECLARE_COMPONENT(TPCDigiAlg)
+
+bool compare_phi( Voxel_tpc* a, Voxel_tpc* b) {
+ return ( a->getPhiIndex() < b->getPhiIndex() ) ;
+}
+
+bool compare_z( Voxel_tpc* a, Voxel_tpc* b) {
+ return ( a->getZIndex() < b->getZIndex() ) ;
+}
+
+
+TPCDigiAlg::TPCDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+ : GaudiAlgorithm(name, svcLoc),
+ _nEvt(0)
+{
+
+ declareProperty("HeaderCol", _headerCol);
+ // modify processor description
+ //_description = "Produces TPC TrackerHit collection from SimTrackerHit collection, smeared in RPhi and Z. A search is made for adjacent hits on a pad row, if they are closer in z and r-phi than the steering parameters _doubleHitResRPhi (default value 2.0 mm) and _doubleHitResZ (default value 5.0 mm) they are considered to overlap. Clusters of hits smaller than _maxMerge (default value 3) are merged into a single tracker hit, with the position given as the average poision of the hits in phi and in z. Clusters which have _maxMerge hits or more are determined to be identifiable as multiple hits, and are not added to the tracker hit collection. This of course means that good hits caught up in a cluster of background hits will be lossed." ;
+
+ // register steering parameters: name, description, class-variable, default value
+ declareProperty("TPCCollection",_padRowHitColHdl,
+ "The default pad-row based SimTrackerHit collection");
+ //registerInputCollection( LCIO::SIMTRACKERHIT,
+ // "TPCPadRowHitCollectionName" ,
+ // "Name of the default pad-row based SimTrackerHit collection" ,
+ // _padRowHitColName ,
+ // std::string("TPCCollection") ) ;
+
+ declareProperty("TPCSpacePointCollection",_spacePointColHdl,
+ "Additional space point collection which provides additional guide hits between pad row centers.");
+ //registerInputCollection( LCIO::SIMTRACKERHIT,
+ // "TPCSpacePointCollectionName" ,
+ // "Name of the additional space point collection which provides additional guide hits between pad row centers." ,
+ // _spacePointColName ,
+ // std::string("TPCSpacePointCollection") ) ;
+
+ declareProperty("TPCLowPtCollection",_lowPtHitsColHdl,
+ "The LowPt SimTrackerHit collection Produced by Mokka TPC Driver TPC0X");
+ //registerInputCollection( LCIO::SIMTRACKERHIT,
+ // "TPCLowPtCollectionName" ,
+ // "Name of the LowPt SimTrackerHit collection Produced by Mokka TPC Driver TPC0X" ,
+ // _lowPtHitscolName ,
+ // std::string("TPCLowPtCollection") ) ;
+ //declareProperty("MCParticle", _mcColHdl, "The intput MCParticle collection");
+
+ declareProperty("TPCTrackerHitsCol",_TPCTrackerHitsColHdl,
+ "The output TrackerHit collection");
+ //registerOutputCollection( LCIO::TRACKERHIT,
+ // "TPCTrackerHitsCol" ,
+ // "Name of the Output TrackerHit collection" ,
+ // _TPCTrackerHitsCol ,
+ // std::string("TPCTrackerHits") ) ;
+
+ declareProperty("TPCTrackerHitAssCol", _TPCAssColHdl, "Handle of TrackerHit SimTrackHit relation collection");
+ //declareProperty("TPCTrackerHitRelations",_outRelCol,
+ // "The TrackerHit SimTrackHit relation collection");
+ //registerOutputCollection(LCIO::LCRELATION,
+ // "SimTrkHitRelCollection",
+ // "Name of TrackerHit SimTrackHit relation collection",
+ // _outRelColName,
+ // std::string("TPCTrackerHitRelations"));
+
+ declareProperty("UseRawHitsToStoreSimhitPointer",
+ _use_raw_hits_to_store_simhit_pointer=bool(false),
+ "Store the pointer to the SimTrackerHits in RawHits (deprecated) ");
+ //registerProcessorParameter("UseRawHitsToStoreSimhitPointer",
+ // "Store the pointer to the SimTrackerHits in RawHits (deprecated) ",
+ // _use_raw_hits_to_store_simhit_pointer,
+ // bool(false));
+
+ declareProperty("PointResolutionPadPhi",_pointResoPadPhi=0.900,
+ "Pad Phi Resolution constant in TPC");
+ //registerProcessorParameter( "PointResolutionPadPhi" ,
+ // "Pad Phi Resolution constant in TPC" ,
+ // _pointResoPadPhi ,
+ // (float)0.900) ;
+
+ declareProperty("RejectCellID0",_rejectCellID0=1,
+ "whether or not to use hits without proper cell ID (pad row)");
+ //registerProcessorParameter( "RejectCellID0" ,
+ // "whether or not to use hits without proper cell ID (pad row)" ,
+ // _rejectCellID0 ,
+ // (int)1) ;
+
+ declareProperty("PointResolutionRPhi",_pointResoRPhi0=0.050,
+ "R-Phi Resolution constant in TPC");
+ //registerProcessorParameter( "PointResolutionRPhi" ,
+ // "R-Phi Resolution constant in TPC" ,
+ // _pointResoRPhi0 ,
+ // (float)0.050) ;
+
+ declareProperty("DiffusionCoeffRPhi",_diffRPhi=0.025,
+ "R-Phi Diffusion Coefficent in TPC");
+ //registerProcessorParameter( "DiffusionCoeffRPhi" ,
+ // "R-Phi Diffusion Coefficent in TPC" ,
+ // _diffRPhi ,
+ // (float)0.025) ;
+
+ declareProperty("N_eff",_nEff=22,
+ "Number of Effective electrons per pad in TPC");
+ //registerProcessorParameter( "N_eff" ,
+ // "Number of Effective electrons per pad in TPC" ,
+ // _nEff ,
+ // (int)22) ;
+
+ declareProperty("PointResolutionZ",_pointResoZ0=0.4,
+ "TPC Z Resolution Coefficent independent of diffusion");
+ //registerProcessorParameter( "PointResolutionZ" ,
+ // "TPC Z Resolution Coefficent independent of diffusion" ,
+ // _pointResoZ0 ,
+ // (float)0.4) ;
+
+ declareProperty("DiffusionCoeffZ",_diffZ=0.08,
+ "Z Diffusion Coefficent in TPC");
+ //registerProcessorParameter( "DiffusionCoeffZ" ,
+ // "Z Diffusion Coefficent in TPC" ,
+ // _diffZ ,
+ // (float)0.08) ;
+
+ declareProperty("HitSortingBinningZ",_binningZ=5.0,
+ "Defines spatial slice in Z");
+ //registerProcessorParameter( "HitSortingBinningZ" ,
+ // "Defines spatial slice in Z" ,
+ // _binningZ ,
+ // (float)5.0) ;
+
+ declareProperty("HitSortingBinningRPhi",_binningRPhi=2.0,
+ "Defines spatial slice in RP");
+ //registerProcessorParameter( "HitSortingBinningRPhi" ,
+ // "Defines spatial slice in RP" ,
+ // _binningRPhi ,
+ // (float)2.0) ;
+
+
+ declareProperty("DoubleHitResolutionZ",_doubleHitResZ=5.0,
+ "Defines the minimum distance for two seperable hits in Z");
+ //registerProcessorParameter( "DoubleHitResolutionZ" ,
+ // "Defines the minimum distance for two seperable hits in Z" ,
+ // _doubleHitResZ ,
+ // (float)5.0) ;
+
+ declareProperty("DoubleHitResolutionRPhi",_doubleHitResRPhi=2.0,
+ "Defines the minimum distance for two seperable hits in RPhi");
+ //registerProcessorParameter( "DoubleHitResolutionRPhi" ,
+ // "Defines the minimum distance for two seperable hits in RPhi" ,
+ // _doubleHitResRPhi ,
+ // (float)2.0) ;
+
+ declareProperty("MaxClusterSizeForMerge",_maxMerge=3,
+ "Defines the maximum number of adjacent hits which can be merged");
+ //registerProcessorParameter( "MaxClusterSizeForMerge" ,
+ // "Defines the maximum number of adjacent hits which can be merged" ,
+ // _maxMerge ,
+ // (int)3) ;
+}
+
+
+StatusCode TPCDigiAlg::initialize()
+{
+ debug() << "in TPCDigiAlg initialize()" <<endmsg;
+
+// // From GNU documentation:
+// // A replacement for the standard terminate_handler which prints
+// // more information about the terminating exception (if any) on stderr. Call ...
+// std::set_terminate (__gnu_cxx::__verbose_terminate_handler);
+//
+//#ifdef DIGIPLOTS
+// /// Hook an AIDA implementation -----------------------------------------------
+//
+// // First create a pointer to the "IAnalysisFactory" of a specific AIDA
+// // implementation. This factory can then be used to produce all other
+// // factories.
+// _AF = AIDA_createAnalysisFactory();
+//
+// // Create a ITreeFactory. -----------------------------------------------------
+// // A ITree can be used to store AIDA objects in memory or on disk.
+//
+// _TRF = _AF->createTreeFactory();
+//
+// /// Create a ITree object which is bound to a file. ---------------------------
+// // You must always create a "ITree" object to create any other factory.
+// /*
+// * Creates a new tree and associates it with a store.
+// * The store is assumed to be read/write.
+// * The store will be created if it does not exist.
+// * @param storeName The name of the store, if empty (""), the tree is
+// * created in memory and therefore will not be associated
+// * with a file.
+// * @param storeType Implementation specific string, may control store type
+// * @param readOnly If true the store is opened readonly, an exception if it
+// * does not exist
+// * @param createNew If false the file must exist, if true the file will be
+// * created
+// * @param options Other options, currently are not specified
+// */
+// // ITree * ITreeFactory::create(const std::string & storeName,
+// // const std::string & storeType = "",
+// // bool readOnly = false,
+// // bool createNew = false,
+// // const std::string & options = "") ;
+//
+// _TREE = _TRF->create("TPCDigi.root",
+// "root",
+// false,
+// true);
+//
+// /// Create an IHistogramFactory which is bound to the tree "*_TREE". -----------
+//
+// /*
+// * Create an IHistogramFactory.
+// * @param tree The ITree which created histograms will be associated to.
+// * @return The IHistogramFactory.
+// */
+// // IHistogramFactory * IAnalysisFactory::createHistogramFactory(ITree & tree);
+//
+// _HF = _AF->createHistogramFactory(*_TREE);
+//
+// _TREE->mkdir("Histograms");
+//
+// /*
+// * Create a IHistogram1D.
+// * @param path The path of the created IHistogram. The path can either
+// * be a relative or full path.
+// * ("/folder1/folder2/dataName" and
+// * "../folder/dataName" are valid paths).
+// * All the directories in the path must exist. The
+// * characther `/` cannot be used in names; it is only
+// * used to delimit directories within paths.
+// * @param title The title of the IHistogram1D.
+// * @param nBins The number of bins of the x axis.
+// * @param lowerEdge The lower edge of the x axis.
+// * @param upperEdge The upper edge of the x axis.
+// * @param options The options for the IHistogram1D. The default is "".
+// * "type=efficiency" for an efficiency IHistogram1D.
+// * @return The newly created IHistogram1D.
+// */
+//
+//
+//
+// _phiDiffHisto = _HF->createHistogram1D("Histograms/phi_diff",
+// "Calculated Phi - Track Phi",
+// 201, -0.05, 0.05);
+//
+// _thetaDiffHisto = _HF->createHistogram1D("Histograms/theta_diff",
+// "Calculated Theta - Track Theta",
+// 201, -0.05, 0.05);
+//
+// _phiRelHisto = _HF->createHistogram1D("Histograms/padPhi",
+// "Phi Relative to the Pad",
+// 201, 0.0, 6.3);
+//
+// _thetaRelHisto = _HF->createHistogram1D("Histograms/padtheta",
+// "Theta Relative to the pad",
+// 201, 0.0, 6.3);
+//
+// _rPhiDiffHisto = _HF->createHistogram1D("Histograms/rPhiDiff",
+// "rPhi_rec - rPhi_sim",
+// 201, -1.0, 1.0);
+//
+// _zDiffHisto = _HF->createHistogram1D("Histograms/zDiff",
+// "Z_rec - Z_sim",
+// 201, -1.0, 1.0);
+//
+// _zPullHisto = _HF->createHistogram1D("Histograms/zPull",
+// "(z_rec - z_sim) / Sigma_z",
+// 201, -10.0, 10.0);
+//
+// _phiDistHisto = _HF->createHistogram1D("Histograms/phiDist",
+// "phi_rec - Phi_sim",
+// 201, -1.0, 1.0);
+//
+// _rPhiPullHisto = _HF->createHistogram1D("Histograms/rPhiPull",
+// "(rPhi_rec - rPhi_sim) / Sigma_rPhi",
+// 201, -10.0, 10.0);
+//
+// _zSigmaVsZHisto = _HF->createHistogram2D("Histograms/zSigmaVsZ",
+// "z Sigma vs Z ",
+// 3000, 0.0, 3000.0,
+// 201, -0.20, 5.20);
+//
+// _zSigmaHisto = _HF->createHistogram1D("Histograms/zSigma",
+// "z Sigma ",
+// 201, -0.20, 5.20);
+//
+// _rPhiSigmaHisto = _HF->createHistogram1D("Histograms/rPhiSigma",
+// "rPhi Sigma",
+// 201, -0.20, 0.20);
+//
+// _radiusCheckHisto = _HF->createHistogram1D("Histograms/radiusCheck",
+// "R_hit - TPC Rmin - ((RowIndex + 0.5 )* padheight)",
+// 201, -0.20, 0.20);
+//
+// _ResidualsRPhiHisto = _HF->createHistogram1D("Histograms/ResidualsRPhi",
+// "MC Track Phi - Hit Phi",
+// 50, -0.001, 0.001);
+//
+// _NSimTPCHitsHisto = _HF->createHistogram1D("Histograms/SimTPCHits",
+// "Number of SimTPC Hits",
+// 100, 0.0, 1000000.0);
+//
+// _NBackgroundSimTPCHitsHisto = _HF->createHistogram1D("Histograms/NBackgroundSimTPCHits",
+// "Number of Background SimTPC Hits",
+// 100, 0.0, 1000000.0);
+//
+// _NPhysicsSimTPCHitsHisto = _HF->createHistogram1D("Histograms/NPhysicsSimTPCHits",
+// "Number of PhysicsSimTPC Hits",
+// 100, 0.0, 100000.0);
+//
+// _NPhysicsAbove02GeVSimTPCHitsHisto = _HF->createHistogram1D("Histograms/NPhysicsAbove02GeVTPCHits",
+// "Number of PhysicsSimTPC Hits above 0.2GeV pt",
+// 100, 0.0, 100000.0);
+//
+// _NPhysicsAbove1GeVSimTPCHitsHisto = _HF->createHistogram1D("Histograms/NPhysicsAbove1GeVPtTPCHits",
+// "Number of PhysicsSimTPC Hits above 1.0 GeV pt",
+// 100, 0.0, 100000.0);
+//
+// _NRecTPCHitsHisto = _HF->createHistogram1D("Histograms/NRecTPCHits",
+// "Number of Rec TPC Hits",
+// 50, 0.0, 100000.0);
+//
+// _NLostPhysicsTPCHitsHisto = _HF->createHistogram1D("Histograms/NLostPhysicsTPCHits",
+// "Number of PhysicsSimTPC Hits Lost",
+// 100, 0.0, 5000.0);
+//
+// _NLostPhysicsAbove02GeVPtTPCHitsHisto = _HF->createHistogram1D("Histograms/NLostPhysicsAbove02GeVPtTPCHits",
+// "Number of PhysicsSimTPC Hits Lost above 0.2 GeV pt",
+// 100, 0.0, 5000.0);
+//
+// _NLostPhysicsAbove1GeVPtTPCHitsHisto = _HF->createHistogram1D("Histograms/NLostPhysicsAbove1GeVPtTPCHits",
+// "Number of PhysicsSimTPC Hits Lost above 1.0 GeV pt",
+// 100, 0.0, 1000.0);
+//
+// _NRevomedHitsHisto = _HF->createHistogram1D("Histograms/NRevomedHits",
+// "Number of Removed TPC hits",
+// 100, 0.0, 1000000.0);
+//
+//
+// _NKeptPhysicsTPCHitsHistoPercent = _HF->createHistogram1D("Histograms/NKeptPhysicsTPCHitsPercent",
+// "Number of PhysicsSimTPC Hits Kept",
+// 303, 0.0, 1.01);
+//
+// _NKeptPhysicsAbove02GeVPtTPCHitsHistoPercent = _HF->createHistogram1D("Histograms/NKeptPhysicsAbove02GeVPtTPCHitsPercent",
+// "Number of PhysicsSimTPC Hits Kept above 0.2 GeV pt",
+// 303, 0.0, 1.01);
+//
+// _NKeptPhysicsAbove1GeVPtTPCHitsHistoPercent = _HF->createHistogram1D("Histograms/NKeptPhysicsAbove1GeVPtTPCHitsPercent",
+// "Number of PhysicsSimTPC Hits Kept above 1.0 GeV pt",
+// 303, 0.0, 1.01);
+//
+//#endif
+//
+//
+// printParameters() ;
+
+ // get the GEAR manager
+ auto _gear = service<IGearSvc>("GearSvc");
+ if ( !_gear ) {
+ error() << "Failed to find GearSvc ..." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ _GEAR = _gear->getGearMgr();
+
+
+ // initialize gsl random generator
+ _random = gsl_rng_alloc(gsl_rng_ranlxs2);
+ _SEEDER = service<IEventSeeder>("EventSeeder");
+ _SEEDER->registerAlg(this);
+
+ _cellid_encoder = 0 ;
+ _nRun = 0 ;
+ _nEvt = 0 ;
+
+ debug()<<" __LINE__"<<endmsg;
+ //_cellid_encoder = new BitField64( lcio::ILDCellID0::encoder_string ) ;
+ debug()<<" __LINE__"<<endmsg;
+ return GaudiAlgorithm::initialize();
+}
+
+//TPCDigiAlg::processRunHeader( LCRunHeader* run)
+//{
+// _nRun++ ;
+//}
+
+StatusCode TPCDigiAlg::execute()
+{
+ debug() << "in TPCDigiAlg execute()" <<endmsg;
+
+ //auto header = _headerCol.get()->at(0);
+ //int evtNo = header.getEventNumber();
+ //int runNo = header.getRunNumber();
+
+ unsigned int thisSeed = _SEEDER->getSeed(this, _nEvt, 0);
+ gsl_rng_set( _random, thisSeed);
+
+ //debug() << "seed set to " << thisSeed << " for event number "<< evtNo << endmsg;
+
+ int numberOfVoxelsCreated(0);
+
+ _NSimTPCHits = 0;
+ _NBackgroundSimTPCHits = 0;
+ _NPhysicsSimTPCHits = 0;
+ _NPhysicsAbove02GeVSimTPCHits = 0;
+ _NPhysicsAbove1GeVSimTPCHits = 0;
+ _NRecTPCHits = 0;
+
+ _NLostPhysicsTPCHits = 0;
+ _NLostPhysicsAbove02GeVPtTPCHits = 0;
+ _NLostPhysicsAbove1GeVPtTPCHits = 0;
+ _NRevomedHits = 0;
+
+ static bool firstEvent = true;
+ _tpcHitMap.clear();
+ _tpcRowHits.clear();
+
+ // fg: make sure we have one message in the log file with run and event number for the DBD production ....
+ info() << " ========= processing event "
+ << std::setw(9) << _nEvt/*evtNo*/ << " run "
+ << std::setw(9) << 0/*runNo*/
+ << " ========= " << endmsg;
+
+
+ if(firstEvent==true) {
+ if (! _use_raw_hits_to_store_simhit_pointer ) {
+
+ debug() << "The relations to SimTrackerHits are now stored in relation collection TPCTrackerHitRelations\n SimTrackerHits are no longer stored in RawTrackerHits. Enable this deprecated feature by setting UseRawHitsToStoreSimhitPointer to true in steering file." << endmsg;
+
+ }
+ else{
+
+ debug() << "SimTrackerHits will be stored in RawTrackerHits. This is a deprecated please use the relations stored in TPCTrackerHitRelations" << endmsg;
+
+ }
+
+
+ }
+
+ firstEvent = false ;
+
+ const gear::TPCParameters& gearTPC = _GEAR->getTPCParameters() ;
+ const gear::PadRowLayout2D& padLayout = gearTPC.getPadLayout() ;
+ // this gets the center of the first pad in the pad layout
+ const gear::Vector2D padCoord = padLayout.getPadCenter(1) ;
+
+ // this assumes that the pad width in r*phi is the same for all pads
+ _padWidth = padLayout.getPadWidth(0)*padCoord[0];
+ // set size of row_hits to hold (n_rows) vectors
+ _tpcRowHits.resize(padLayout.getNRows());
+
+ //// created the collection which will be written out
+ _trkhitVec = _TPCTrackerHitsColHdl.createAndPut();
+ _relCol = _TPCAssColHdl.createAndPut();//LCRELATION
+
+ //zhang TODO
+ // to store the weights
+ //LCFlagImpl lcFlag(0) ;
+ //lcFlag.setBit( LCIO::LCREL_WEIGHTED ) ;
+ //_relCol->setFlag( lcFlag.getFlag() ) ;
+ //
+
+ //auto mcCol = _mcColHdl.get();
+ //auto mcp = mcCol->at(0);
+ //debug() << "First MCParticle " << mcp << endmsg;
+
+ // first deal with the pad-row based hits from Mokka
+ const edm4hep::SimTrackerHitCollection* STHcol = nullptr;
+ try {
+ STHcol = _padRowHitColHdl.get();
+ }
+ catch ( GaudiException &e ) {
+ debug() << "Collection " << _padRowHitColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ return StatusCode::SUCCESS;
+ }
+
+ _cellid_encoder = new UTIL::BitField64( lcio::ILDCellID0::encoder_string ) ;
+ //int det_id = 0 ;
+ //if ( (STHcol != nullptr) && (STHcol->size()>0) ) {
+ // auto SimTHit = STHcol->at( 0 ) ;
+ // _cellid_encoder->setValue(SimTHit.getCellID()) ;
+ // if ( (*_cellid_encoder)[lcio::ILDCellID0::subdet]!= ILDDetID::TPC ){
+ // //fatal() << "unsupported detector ID NOT TPC. det_id = " << det_id << endmsg;
+ // return StatusCode::FAILURE;
+ // }
+ //}else{
+ // debug() << "Collection " << _padRowHitCol.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ // return StatusCode::SUCCESS;
+ //}
+
+ if( STHcol != nullptr ){
+
+ int n_sim_hits = STHcol->size() ;
+
+ //TODO
+ //LCFlagImpl colFlag( STHcol->getFlag() ) ;
+
+ _NSimTPCHits = n_sim_hits;
+
+ debug() << "number of Pad-Row based SimHits = " << n_sim_hits << endmsg;
+
+ edm4hep::ConstMCParticle nMinus2MCP;
+ edm4hep::ConstMCParticle previousMCP;
+ edm4hep::ConstSimTrackerHit nMinus2SimHit;
+ edm4hep::ConstSimTrackerHit previousSimTHit;
+
+ debug() << "processing nhit=" << n_sim_hits << endmsg;
+ // loop over all the pad row based sim hits
+ for(unsigned int i = 0; i<n_sim_hits; i++){
+ auto SimTHit = STHcol->at(i);
+
+ // this will used for nominaml smearing for very low pt rubish, so set it to zero initially
+ double ptSqrdMC = 0;
+
+ float edep;
+ double padPhi(0.0);
+ double padTheta (0.0);
+
+ auto& pos = SimTHit.getPosition();
+ debug() << "processing hit id = " << SimTHit.id() << endmsg;
+ debug() << " SimTHit= "
+ << " x = " << pos[0]
+ << " y = " << pos[1]
+ << " z = " << pos[2]
+ << endmsg;
+
+ CLHEP::Hep3Vector thisPoint(pos[0],pos[1],pos[2]);
+ double padheight = padLayout.getPadHeight(padLayout.getNearestPad(thisPoint.perp(),thisPoint.phi()));
+
+ const double bField = _GEAR->getBField().at( gear::Vector3D( 0., 0., 0.) ).z() ;
+ // conversion constant. r = pt / (FCT*bField)
+ const double FCT = 2.99792458E-4;
+ bool found_mc = false;
+ edm4hep::ConstMCParticle mcp;
+ try{ // protect crash while MCParticle unavailable
+ mcp = SimTHit.getMCParticle() ;
+ }
+ catch(...){
+ debug() << "catch throw MCParticle not available" << endmsg;
+ }
+ // increase the counters for the different classification of simhits
+ //yzhang FIXME mcp!=NULL
+ if(mcp.isAvailable()){
+
+ // get the pt of the MCParticle, this will used later to uses nominal smearing for low momentum rubish
+ const edm4hep::Vector3f momentumMC= mcp.getMomentum() ;
+ ptSqrdMC = momentumMC[0]*momentumMC[0]+momentumMC[1]*momentumMC[1] ;
+
+ debug() << " mcp id = " << mcp.id()
+ << " px = " << momentumMC[0]
+ << " py = " << momentumMC[1]
+ << " pz = " << momentumMC[2]
+ << endmsg;
+
+ // SJA:FIXME: the fact that it is a physics hit relies on the fact that for overlay
+ // the pointer to the mcp is set to NULL. This distinction may not always be true ...
+ ++_NPhysicsSimTPCHits ;
+ if( ptSqrdMC > (0.2*0.2) ) ++_NPhysicsAbove02GeVSimTPCHits ;
+ if( ptSqrdMC > 1.0 ) ++_NPhysicsAbove1GeVSimTPCHits ;
+
+#ifdef DIGIPLOTS
+// if(mcp) plotHelixHitResidual(mcp, &thisPoint);
+#endif
+ } else {
+ debug() << "MCParticle not available" << endmsg;
+ ++_NBackgroundSimTPCHits;
+ }
+
+ // if the hits contain the momentum of the particle use this to calculate the angles relative to the pad
+ //yzhang TODO
+ //colFlag.bitSet(LCIO::THBIT_MOMENTUM) == true
+ //if(colFlag.bitSet(LCIO::THBIT_MOMENTUM))
+
+ bool colFlag_bitSet = true;
+ if(colFlag_bitSet){
+ // FIXME: fucd, since momentum from SimTrackerHit available, change MCParticle's to SimTrackerHit's, which better?
+ //const edm4hep::Vector3f mcpMomentum = mcp.getMomentum() ;
+ const edm4hep::Vector3f mcpMomentum = SimTHit.getMomentum();
+
+ CLHEP::Hep3Vector mom(mcpMomentum[0],mcpMomentum[1],mcpMomentum[2]);
+ // const double pt = mom.perp();
+ // const double radius = pt / (FCT*bField);
+
+ // const double tanLambda = mom.z()/pt;
+ padPhi = fabs(thisPoint.deltaPhi(mom));
+ padTheta = mom.theta();
+ }
+ else { // LCIO::THBIT_MOMENTUM not set
+
+ // as the momentum vector is not available from the hits use triplets of
+ // hits to fit a circle and calculate theta and phi relative to the pad
+
+ if ((!mcp.isAvailable()) || (sqrt(ptSqrdMC) / (FCT*bField)) < ( padheight / (0.1 * CLHEP::twopi))) {
+ // if the hit has no record of it MCParticle then there is no way to know if this hit has consecutive hits from the same MCParticle
+ // so just set nominal values theta=phi=90
+ // here make a cut for particles which will suffer more than a 10 percent change in phi over the distance of the pad
+ // R > padheight/(0.1*2PI)
+ // in both cases set the angles to 90 degrees
+ padTheta = CLHEP::twopi/4.0 ;
+ padPhi = CLHEP::twopi/4.0 ;
+ }
+ else{
+ edm4hep::ConstSimTrackerHit nextSimTHit;
+ edm4hep::ConstSimTrackerHit nPlus2SimHit;
+ edm4hep::ConstMCParticle nextMCP;
+ edm4hep::ConstMCParticle nPlus2MCP;
+ // if there is at least one more hit after this one, set the pointer to the MCParticle for the next hit
+ if (i < (n_sim_hits-1) ) {
+ nextSimTHit = STHcol->at( i+1 ) ;
+ nextMCP = nextSimTHit.getMCParticle() ;
+ }
+ else{ // set make sure that the pointers are set back to NULL so that the comparisons later hold
+ //nextSimTHit = edm4hep::ConstSimTrackerHit;
+ //nextMCP = edm4hep::ConstMCParticle;
+ }
+ // if there is at least two more hits after this one, set the pointer to the MCParticle for the next but one hit
+ if (i < (n_sim_hits-2) ) {
+ nPlus2SimHit = STHcol->at( i+2 );
+ nPlus2MCP = nPlus2SimHit.getMCParticle() ;
+ }
+ else{ // set make sure that the pointers are set back to NULL so that the comparisons later hold
+ //_nPlus2SimHit = edm4hep::ConstSimTrackerHit;
+ //_nPlus2MCP = edm4hep::ConstMCParticle;
+ }
+
+ if ( mcp==previousMCP && mcp==nextMCP ) { // middle hit of 3 from the same MCParticle
+
+ CLHEP::Hep3Vector precedingPoint(previousSimTHit.getPosition()[0],previousSimTHit.getPosition()[1],previousSimTHit.getPosition()[2]) ;
+ CLHEP::Hep3Vector followingPoint(nextSimTHit.getPosition()[0],nextSimTHit.getPosition()[1],nextSimTHit.getPosition()[2]) ;
+
+ debug() << "address of _previousSimTHit = " << previousSimTHit
+ << " x = " << previousSimTHit.getPosition()[0]
+ << " y = " << previousSimTHit.getPosition()[1]
+ << " z = " << previousSimTHit.getPosition()[2]
+ << std::endl ;
+
+ debug() << "address of _nextSimTHit = " << nextSimTHit
+ << " x = " << nextSimTHit.getPosition()[0]
+ << " y = " << nextSimTHit.getPosition()[1]
+ << " z = " << nextSimTHit.getPosition()[2]
+ << std::endl ;
+
+ // get phi and theta using functions defined below
+ padPhi = getPadPhi( &thisPoint, &precedingPoint, &thisPoint, &followingPoint);
+ padTheta = getPadTheta(&precedingPoint, &thisPoint, &followingPoint);
+
+ }
+ else if ( mcp==nextMCP && mcp==nPlus2MCP ) { // first hit of 3 from the same MCParticle
+
+ CLHEP::Hep3Vector followingPoint(nextSimTHit.getPosition()[0],nextSimTHit.getPosition()[1],nextSimTHit.getPosition()[2]) ;
+ CLHEP::Hep3Vector nPlus2Point(nPlus2SimHit.getPosition()[0],nPlus2SimHit.getPosition()[1],nPlus2SimHit.getPosition()[2]) ;
+
+ // get phi and theta using functions defined below
+ padPhi = getPadPhi( &thisPoint, &thisPoint, &followingPoint, &nPlus2Point);
+ padTheta = getPadTheta(&thisPoint, &followingPoint, &nPlus2Point);
+
+ }
+ else if ( mcp==previousMCP && mcp==nMinus2MCP ) { // last hit of 3 from the same MCParticle
+ auto& posMinus = nMinus2SimHit.getPosition();
+ auto& posPrevious = previousSimTHit.getPosition();
+
+ CLHEP::Hep3Vector nMinus2Point(posMinus[0],posMinus[1],posMinus[2]);
+ CLHEP::Hep3Vector precedingPoint(posPrevious[0],posPrevious[1],posPrevious[2]);
+
+ // get phi and theta using functions defined below
+ padPhi = getPadPhi( &thisPoint, &nMinus2Point, &precedingPoint, &thisPoint);
+ padTheta = getPadTheta(&nMinus2Point, &precedingPoint, &thisPoint);
+
+ }
+ else{ // the hit is isolated as either a single hit, or a pair of hits, from a single MCParticle
+ padTheta = CLHEP::twopi/4.0 ;
+ padPhi = CLHEP::twopi/4.0 ;
+ }
+ }
+
+//#ifdef DIGIPLOTS
+// if(colFlag.bitSet(LCIO::THBIT_MOMENTUM)) {
+//
+// const float * mcpMomentum = SimTHit.getMomentum() ;
+//
+// CLHEP::Hep3Vector mom(mcpMomentum[0],mcpMomentum[1],mcpMomentum[2]);
+//
+// double trackPhi = mom.phi();
+//
+// if(trackPhi<0.0) trackPhi=trackPhi+twopi;
+// if(trackPhi>twopi) trackPhi=trackPhi-twopi;
+// if(trackPhi>twopi/2.0) trackPhi = trackPhi - twopi/2.0 ;
+//
+// double localPhi = thisPoint.phi() - padPhi;
+//
+// _phiRelHisto->fill(padPhi);
+// _phiDiffHisto->fill((fabs(localPhi - trackPhi))/trackPhi);
+// _thetaRelHisto->fill(padTheta);
+// _thetaDiffHisto->fill( (sin(padTheta) - sin(mom.theta()))/sin(mom.theta()) );
+//
+// debug() << "track Phi = " << trackPhi * (360.0 / twopi) << endmsg;
+// debug() << "localPhi = " << localPhi * (360.0 / twopi) << endmsg;
+// debug() << "pad Phi = " << padPhi * (360.0 / twopi) << endmsg;
+// debug() << "pad Phi from track mom = " << ( thisPoint.phi() - trackPhi ) * (360.0 / twopi) << endmsg;
+// debug() << "padTheta = " << padTheta * (360.0 / twopi) << endmsg;
+// debug() << "padTheta from track mom = " << mom.theta() * (360.0 / twopi) << endmsg;
+//
+// }
+//#endif
+
+ }
+ debug() << "padPhi = " << padPhi << " padTheta = " << padTheta << endmsg;
+ // int pad = padLayout.getNearestPad(thisPoint.perp(),thisPoint.phi());
+ int layerNumber = SimTHit.getCellID();
+
+ if(_rejectCellID0 && (layerNumber<1)) {
+ continue;
+ }
+
+ edep = SimTHit.getEDep();
+
+ // Calculate Point Resolutions according to Ron's Formula
+
+ // sigma_{RPhi}^2 = sigma_0^2 + Cd^2/N_{eff} * L_{drift}
+
+ // sigma_0^2 = (50micron)^2 + (900micron*sin(phi))^2
+ // Cd^2/N_{eff}} = 25^2/(22/sin(theta)*h/6mm)
+ // Cd = 25 ( microns / cm^(1/2) )
+ // (this is for B=4T, h is the pad height = pad-row pitch in mm,
+ // theta is the polar angle)
+
+ // sigma_{z}^2 = (400microns)^2 + L_{drift}cm * (80micron/sqrt(cm))^2
+
+ double aReso =_pointResoRPhi0*_pointResoRPhi0 + (_pointResoPadPhi*_pointResoPadPhi * sin(padPhi)*sin(padPhi)) ;
+ double driftLength = gearTPC.getMaxDriftLength() - (fabs(thisPoint.z()));
+
+ if (driftLength <0) {
+ debug() << " TPCDigiAlg : Warning! driftLength < 0 " << driftLength << " --> Check out your GEAR file!!!!" << endmsg;
+ debug() << "Setting driftLength to 0.1" << endmsg;
+ debug() << "gearTPC.getMaxDriftLength() = " << gearTPC.getMaxDriftLength() << endmsg;
+ driftLength = 0.10;
+ }
+
+ padheight = padLayout.getPadHeight(padLayout.getNearestPad(thisPoint.perp(),thisPoint.phi()));
+
+ double bReso = ( (_diffRPhi * _diffRPhi) / _nEff ) * sin(padTheta) * ( 6.0 / (padheight) ) * ( 4.0 / bField ) ;
+
+ double tpcRPhiRes = sqrt( aReso + bReso * (driftLength / 10.0) ); // driftLength in cm
+
+ double tpcZRes = sqrt(( _pointResoZ0 * _pointResoZ0 )
+ +
+ ( _diffZ * _diffZ ) * (driftLength / 10.0) ); // driftLength in cm
+
+ int padIndex = padLayout.getNearestPad(thisPoint.perp(),thisPoint.phi());
+
+ const gear::DoubleVec & planeExt = padLayout.getPlaneExtent() ;
+ double TPCPadPlaneRMin = planeExt[0] ;
+ double TPCPadPlaneRMax = planeExt[1] ;
+
+ int iRowHit = padLayout.getRowNumber(padIndex);
+ int iPhiHit = padLayout.getPadNumber(padIndex);
+ int NBinsZ = (int) ((2.0 * gearTPC.getMaxDriftLength()) / _binningZ);
+ int iZHit = (int) ( (float) NBinsZ * ( gearTPC.getMaxDriftLength() + thisPoint.z() ) / ( 2.0 * gearTPC.getMaxDriftLength() ) ) ;
+
+ if(iZHit<0) iZHit=0;
+ if(iZHit>NBinsZ) iZHit=NBinsZ;
+
+ // make sure that the hit lies at the middle of the pad ring
+ thisPoint.setPerp(padLayout.getPadCenter(padIndex)[0]);
+
+ if( (thisPoint.perp() < TPCPadPlaneRMin) || (thisPoint.perp() > TPCPadPlaneRMax) ) {
+ debug() << "Hit R not in TPC " << endmsg;
+ debug() << "R = " << thisPoint.perp() << endmsg;
+ debug() << "the tpc InnerRadius = " << TPCPadPlaneRMin << endmsg;
+ debug() << "the tpc OuterRadius = " << TPCPadPlaneRMax << endmsg;
+ debug() << "Hit Dropped " << endmsg;
+ continue;
+ }
+
+ if( (fabs(thisPoint.z()) > gearTPC.getMaxDriftLength()) ) {
+ debug() << "Hit Z not in TPC " << endmsg;
+ debug() << "Z = " << thisPoint.z() << endmsg;
+ debug() << "the tpc Max Z = " << gearTPC.getMaxDriftLength() << endmsg;
+ debug() << "Hit Dropped " << endmsg;
+ continue;
+ }
+
+ // create a tpc voxel hit and store it for this row
+ Voxel_tpc * atpcVoxel = new Voxel_tpc(iRowHit,iPhiHit,iZHit, thisPoint, edep, tpcRPhiRes, tpcZRes);
+
+ _tpcRowHits.at(iRowHit).push_back(atpcVoxel);
+ ++numberOfVoxelsCreated;
+
+ // store the simhit pointer for this tpcvoxel hit in the hit map
+ _tpcHitMap[atpcVoxel] = SimTHit;
+
+ // move the pointers on
+ nMinus2MCP = previousMCP;
+ previousMCP = mcp ;
+ nMinus2SimHit = previousSimTHit;
+ previousSimTHit = SimTHit;
+
+ }
+ }
+
+ // now process the LowPt collection
+ try {
+ auto STHcolLowPt = _lowPtHitsColHdl.get();
+
+ if(nullptr!=STHcolLowPt){
+
+ int n_sim_hitsLowPt = STHcolLowPt->size() ;
+
+ _NBackgroundSimTPCHits += n_sim_hitsLowPt;
+ _NSimTPCHits += n_sim_hitsLowPt;
+
+ _padWidth = padLayout.getPadWidth(0)*padCoord[0];
+
+ debug() << "number of LowPt hits:" << n_sim_hitsLowPt << endmsg;
+
+ // loop over the LowPt hit collection
+ for(auto SimTHit : *STHcolLowPt){
+
+ CLHEP::Hep3Vector thisPoint(SimTHit.getPosition()[0],SimTHit.getPosition()[1],SimTHit.getPosition()[2]);
+
+ const gear::DoubleVec & planeExt = padLayout.getPlaneExtent() ;
+ double TPCPadPlaneRMin = planeExt[0] ;
+ double TPCPadPlaneRMax = planeExt[1] ;
+
+ int NBinsZ = (int) ((2.0 * gearTPC.getMaxDriftLength()) / _binningZ);
+
+ if( (thisPoint.perp() < TPCPadPlaneRMin) || (thisPoint.perp() > TPCPadPlaneRMax) ) {
+ debug() << "Hit R not in TPC " << endmsg;
+ debug() << "R = " << thisPoint.perp() << endmsg;
+ debug() << "the tpc InnerRadius = " << TPCPadPlaneRMin << endmsg;
+ debug() << "the tpc OuterRadius = " << TPCPadPlaneRMax << endmsg;
+ debug() << "Hit Dropped " << endmsg;
+ continue;
+ }
+
+ if( (fabs(thisPoint.z()) > gearTPC.getMaxDriftLength()) ) {
+ debug() << "Hit Z not in TPC " << endmsg;
+ debug() << "Z = " << thisPoint.z() << endmsg;
+ debug() << "the tpc Max Z = " << gearTPC.getMaxDriftLength() << endmsg;
+ debug() << "Hit Dropped " << endmsg;
+ continue;
+ }
+
+ int padIndex = padLayout.getNearestPad(thisPoint.perp(),thisPoint.phi());
+ //double padheight = padLayout.getPadHeight(padIndex);
+
+ int iRowHit = padLayout.getRowNumber(padIndex);
+ int iPhiHit = padLayout.getPadNumber(padIndex);
+ int iZHit = (int) ( (float) NBinsZ *
+ ( gearTPC.getMaxDriftLength() + thisPoint.z() ) / ( 2.0 * gearTPC.getMaxDriftLength() ) ) ;
+
+ // shift the hit in r-phi to the nearest pad-row centre
+ thisPoint.setPerp(padLayout.getPadCenter(padIndex)[0]);
+
+ // set the resolutions to the pads to digital like values
+ double tpcRPhiRes = _padWidth;
+ double tpcZRes = _binningZ;
+
+ // create a tpc voxel hit for this simhit and store it for this tpc pad row
+ Voxel_tpc * atpcVoxel = new Voxel_tpc(iRowHit,iPhiHit,iZHit, thisPoint, SimTHit.getEDep(), tpcRPhiRes, tpcZRes);
+
+ _tpcRowHits.at(iRowHit).push_back(atpcVoxel);
+ ++numberOfVoxelsCreated;
+
+ // store the simhit pointer for this voxel hit in a map
+ _tpcHitMap[atpcVoxel] = SimTHit;
+
+ }
+ }
+ }catch(GaudiException e){
+ warning() << "Catch exception when read LowPtHitsCol" << endmsg;
+ }
+
+ int number_of_adjacent_hits(0);
+
+ debug() << "finished looping over simhits, number of voxels = " << numberOfVoxelsCreated << endmsg;
+
+ int numberOfhitsTreated(0);
+
+ vector <Voxel_tpc *> row_hits;
+
+ // loop over the tpc rows containing hits and check for merged hits
+ for (unsigned int i = 0; i<_tpcRowHits.size(); ++i){
+
+ row_hits = _tpcRowHits.at(i);
+ std::sort(row_hits.begin(), row_hits.end(), compare_phi );
+
+ // double loop over the hits in this row
+ for (unsigned int j = 0; j<row_hits.size(); ++j){
+
+ ++numberOfhitsTreated;
+
+ for (unsigned int k = j+1; k<row_hits.size(); ++k){
+
+ if(row_hits[k]->getPhiIndex() > (row_hits[j]->getPhiIndex())+2){ // SJA:FIXME: here we need an OR to catch the wrap around
+ break; // only compare hits in adjacent phi bins
+ }
+
+ // look to see if the two hit occupy the same pad in phi or if not whether they are within the r-phi double hit resolution
+ else if( row_hits[k]->getPhiIndex()==row_hits[j]->getPhiIndex()
+ ||
+ ( (fabs(row_hits[k]->getHep3Vector().deltaPhi(row_hits[j]->getHep3Vector()))) * row_hits[j]->getR()) < _doubleHitResRPhi ) {
+
+ // if neighboring in phi then compare z
+ map <Voxel_tpc*,edm4hep::SimTrackerHit> ::iterator it;
+
+ edm4hep::SimTrackerHit Hit1;
+ edm4hep::SimTrackerHit Hit2;
+
+ // search of the simhit pointers in the tpchit map
+ it=_tpcHitMap.find(row_hits[j]);
+ if(it!= _tpcHitMap.end()) {
+ Hit1 = it->second ; // hit found
+ }
+
+ it=_tpcHitMap.find(row_hits[k]);
+ if(it!= _tpcHitMap.end()) {
+ Hit2 = it->second ; // hit found
+ }
+
+ double pathlengthZ1(0.0);
+ double pathlengthZ2(0.0);
+
+ if( Hit1.isAvailable() && Hit2.isAvailable() ){ // if both sim hits were found
+
+ // check if the track momentum has been stored for the hits
+ bool momentum_set = true;
+
+ if( STHcol != NULL ){
+ //yzhang TODO
+ //LCFlagImpl colFlag( STHcol->getFlag() ) ;
+ //momentum_set = momentum_set && colFlag.bitSet(LCIO::THBIT_MOMENTUM) ;
+ }
+
+ //if( STHcolLowPt != NULL ){
+ //yzhang TODO
+ //LCFlagImpl colFlag( STHcolLowPt->getFlag() ) ;
+ //momentum_set = momentum_set && colFlag.bitSet(LCIO::THBIT_MOMENTUM) ;
+ //}
+
+ if( momentum_set ){
+
+ const edm4hep::Vector3f Momentum1 = Hit1.getMomentum() ;
+ const edm4hep::Vector3f Momentum2 = Hit1.getMomentum() ;
+
+ CLHEP::Hep3Vector mom1(Momentum1[0],Momentum1[1],Momentum1[2]);
+ CLHEP::Hep3Vector mom2(Momentum2[0],Momentum2[1],Momentum2[2]);
+
+ pathlengthZ1 = fabs( Hit1.getPathLength() * mom1.cosTheta() );
+ pathlengthZ2 = fabs( Hit2.getPathLength() * mom2.cosTheta() );
+ }
+ else {
+ pathlengthZ1 = _doubleHitResZ ; // assume the worst i.e. that the track is moving in z
+ pathlengthZ2 = _doubleHitResZ ; // assume the worst i.e. that the track is moving in z
+ }
+
+ double dZ = fabs(row_hits[j]->getZ() - row_hits[k]->getZ());
+
+ double spacial_coverage = 0.5*(pathlengthZ1 + pathlengthZ2) + _binningZ;
+
+ if( (dZ - spacial_coverage) < _doubleHitResZ ){
+
+ row_hits[j]->setAdjacent(row_hits[k]);
+ row_hits[k]->setAdjacent(row_hits[j]);
+ ++number_of_adjacent_hits;
+
+ }
+ } else {
+ debug() << "Hit1=" << Hit1 << "Hit2=" << Hit2 << endmsg;
+ }
+ }
+ }
+ }
+
+
+ // now all hits have been checked for adjacent hits, go throught and write out the hits or merge
+
+ for (unsigned int j = 0; j<row_hits.size(); ++j){
+
+ Voxel_tpc* seed_hit = row_hits[j];
+
+ if(seed_hit->IsMerged() || seed_hit->IsClusterHit()) {
+ continue;
+ }
+
+ if(seed_hit->getNumberOfAdjacent()==0){ // no adjacent hits so smear and write to hit collection
+ writeVoxelToHit(seed_hit);
+ } else if(seed_hit->getNumberOfAdjacent() < (_maxMerge)){ // potential 3-hit cluster, can use simple average merge.
+
+ vector <Voxel_tpc*>* hitsToMerge = new vector <Voxel_tpc*>;
+
+ int clusterSize = seed_hit->clusterFind(hitsToMerge);
+
+ if( clusterSize <= _maxMerge ){ // merge cluster
+ seed_hit->setIsMerged();
+ writeMergedVoxelsToHit(hitsToMerge);
+ }
+ delete hitsToMerge;
+ }
+ }
+ }
+
+ int numberOfHits(0);
+ // count up the number of hits merged or lost
+ for (unsigned int i = 0; i<_tpcRowHits.size(); ++i){
+ row_hits = _tpcRowHits.at(i);
+ for (unsigned int j = 0; j<row_hits.size(); ++j){
+ numberOfHits++;
+ Voxel_tpc* seed_hit = row_hits[j];
+ if(seed_hit->IsMerged() || seed_hit->IsClusterHit() || seed_hit->getNumberOfAdjacent() > _maxMerge ) {
+ ++_NRevomedHits;
+ auto mcp = (_tpcHitMap[ seed_hit ]).getMCParticle() ;
+ if(mcp.isAvailable()) {
+ ++_NLostPhysicsTPCHits;
+ const edm4hep::Vector3f mom= mcp.getMomentum() ;
+ double ptSQRD = mom[0]*mom[0]+mom[1]*mom[1] ;
+ if( ptSQRD > (0.2*0.2) ) ++_NLostPhysicsAbove02GeVPtTPCHits ;
+ if( ptSQRD > 1.0 ) ++_NLostPhysicsAbove1GeVPtTPCHits ;
+ }
+ }
+ }
+ }
+
+ debug() << "the number of adjacent hits is " << number_of_adjacent_hits << " _doubleHitResZ " << _doubleHitResZ << endmsg;
+ debug() << "number of rec_hits = " << _NRecTPCHits << endmsg ;
+ debug() << "finished row hits " << numberOfHits << " " << numberOfhitsTreated << endmsg;
+
+ // set the parameters to decode the type information in the collection
+ // for the time being this has to be done manually
+ // in the future we should provide a more convenient mechanism to
+ // decode this sort of meta information
+
+ //StringVec typeNames ;
+ //IntVec typeValues ;
+ //typeNames.push_back( LCIO::TPCHIT ) ;
+ //typeValues.push_back( 1 ) ;
+ //yzhang TODO
+ //_trkhitVec->parameters().setValues("TrackerHitTypeNames" , typeNames ) ;
+ //_trkhitVec->parameters().setValues("TrackerHitTypeValues" , typeValues ) ;
+
+ //// add the collection to the event
+ //evt->addCollection( _trkhitVec , _TPCTrackerHitsCol ) ;
+ //evt->addCollection( _relCol , _outRelColName ) ;
+
+ // delete voxels
+ for (unsigned int i = 0; i<_tpcRowHits.size(); ++i){
+ vector <Voxel_tpc *>* current_row = &_tpcRowHits.at(i);
+ for (unsigned int j = 0; j<current_row->size(); ++j){
+ delete current_row->at(j);
+ }
+ }
+
+//#ifdef DIGIPLOTS
+ //_NSimTPCHitsHisto->fill(_NSimTPCHits);
+ //_NBackgroundSimTPCHitsHisto->fill(_NBackgroundSimTPCHits);
+ //_NPhysicsSimTPCHitsHisto->fill(_NPhysicsSimTPCHits);
+ //_NPhysicsAbove02GeVSimTPCHitsHisto->fill(_NPhysicsAbove02GeVSimTPCHits);
+ //_NPhysicsAbove1GeVSimTPCHitsHisto->fill(_NPhysicsAbove1GeVSimTPCHits);
+ //_NRecTPCHitsHisto->fill(_NRecTPCHits);
+
+ //_NLostPhysicsTPCHitsHisto->fill(_NLostPhysicsTPCHits);
+ //_NLostPhysicsAbove02GeVPtTPCHitsHisto->fill(_NLostPhysicsAbove02GeVPtTPCHits);
+ //_NLostPhysicsAbove1GeVPtTPCHitsHisto->fill(_NLostPhysicsAbove1GeVPtTPCHits);
+ //_NRevomedHitsHisto->fill(_NRevomedHits);
+
+ //_NKeptPhysicsTPCHitsHistoPercent->fill( (float)(_NPhysicsSimTPCHits-_NLostPhysicsTPCHits) / (float)_NPhysicsSimTPCHits );
+ //_NKeptPhysicsAbove02GeVPtTPCHitsHistoPercent->fill( (float)(_NPhysicsAbove02GeVSimTPCHits-_NLostPhysicsAbove02GeVPtTPCHits) / (float)_NPhysicsAbove02GeVSimTPCHits );
+ //_NKeptPhysicsAbove1GeVPtTPCHitsHistoPercent->fill( (float)(_NPhysicsAbove1GeVSimTPCHits-_NLostPhysicsAbove1GeVPtTPCHits) / (float)_NPhysicsAbove1GeVSimTPCHits );
+//#endif
+
+ debug() << "_NSimTPCHits = " << _NSimTPCHits << endmsg;
+ debug() << "_NBackgroundSimTPCHits = " << _NBackgroundSimTPCHits << endmsg;
+ debug() << "_NPhysicsSimTPCHits = " << _NPhysicsSimTPCHits << endmsg;
+ debug() << "_NPhysicsAbove02GeVSimTPCHits = " << _NPhysicsAbove02GeVSimTPCHits << endmsg;
+ debug() << "_NPhysicsAbove1GeVSimTPCHits = " << _NPhysicsAbove1GeVSimTPCHits << endmsg;
+ debug() << "_NRecTPCHits = " << _NRecTPCHits<< endmsg;
+ debug() << "_NLostPhysicsTPCHits = " << _NLostPhysicsTPCHits << endmsg;
+ debug() << "_NLostPhysicsAbove02GeVPtTPCHits = " << _NLostPhysicsAbove02GeVPtTPCHits << endmsg;
+ debug() << "_NLostPhysicsAbove1GeVPtTPCHits = " << _NLostPhysicsAbove1GeVPtTPCHits << endmsg;
+ debug() << "_NRevomedHits = " << _NRevomedHits << endmsg;
+
+ _nEvt++;
+ //Clear the maps and the end of the event.
+ _tpcHitMap.clear();
+ _tpcRowHits.clear();
+
+ delete _cellid_encoder ;
+ return StatusCode::SUCCESS;
+}
+
+
+
+StatusCode TPCDigiAlg::finalize()
+{
+
+#ifdef DIGIPLOTS
+ //_TREE->commit();
+ //_TREE->cd("/Histograms");
+ //_TREE->ls("..");
+
+ //_TREE->close();
+ //info() << "DIGICHECKPLOTS Finished" << endmsg;
+#endif
+
+ gsl_rng_free(_random);
+ info() << "TPCDigiAlg::end() " << name()
+ << " processed " << _nEvt << " events in " << _nRun << " runs "
+ << endl ;
+
+ return StatusCode::SUCCESS;
+}
+
+void TPCDigiAlg::writeVoxelToHit( Voxel_tpc* aVoxel){
+
+ const gear::TPCParameters& gearTPC = _GEAR->getTPCParameters() ;
+ const gear::PadRowLayout2D& padLayout = gearTPC.getPadLayout() ;
+ const gear::Vector2D padCoord = padLayout.getPadCenter(1) ;
+
+ Voxel_tpc* seed_hit = aVoxel;
+
+ // if( seed_hit->getRowIndex() > 5 ) return ;
+ debug() << "==============" << endmsg;
+ //store hit variables
+ edm4hep::TrackerHit trkHit;// = _trkhitVec->create();
+ //now the hit pos has to be smeared
+
+ double tpcRPhiRes = seed_hit->getRPhiRes();
+ double tpcZRes = seed_hit->getZRes();
+
+ CLHEP::Hep3Vector point(seed_hit->getX(),seed_hit->getY(),seed_hit->getZ());
+
+ double unsmearedPhi = point.phi();
+
+ double randrp = gsl_ran_gaussian(_random,tpcRPhiRes);
+ double randz = gsl_ran_gaussian(_random,tpcZRes);
+
+ point.setPhi( point.phi() + randrp/ point.perp() );
+ point.setZ( point.z() + randz );
+
+ // make sure the hit is not smeared beyond the TPC Max DriftLength
+ if( fabs(point.z()) > gearTPC.getMaxDriftLength() ) point.setZ( (fabs(point.z()) / point.z() ) * gearTPC.getMaxDriftLength() );
+ debug() << "==============" << endmsg;
+ edm4hep::Vector3d pos(point.x(),point.y(),point.z());
+ trkHit.setPosition(pos);
+ trkHit.setEDep(seed_hit->getEDep());
+ // trkHit->setType( 500 );
+
+ // int side = lcio::ILDDetID::barrel ;
+ //
+ // if( pos[2] < 0.0 ) side = 1 ;
+ //change to Marlin's, fucd
+ //map<Voxel_tpc*,edm4hep::SimTrackerHit>::iterator it=_tpcHitMap.find(seed_hit);
+ //assert(_tpcHitMap.end() != it);
+
+ //const int celId = it->second.getCellID();
+ //_cellid_encoder->setValue( celId );
+ //trkHit.setCellID( _cellid_encoder->lowWord() );
+
+ //int side = (*_cellid_encoder)[lcio::ILDCellID0::side];
+ //int layer = (*_cellid_encoder)[lcio::ILDCellID0::layer];
+ //int module = (*_cellid_encoder)[lcio::ILDCellID0::module];
+ //int sensor = (*_cellid_encoder)[lcio::ILDCellID0::sensor];
+
+ //debug() << "Hit = " << " has celId " << celId << endmsg;
+ //debug() << "side = " << side << endmsg;
+ //debug() << "layerNumber = " << layer << endmsg;
+ //debug() << "moduleNumber = " << module << endmsg;
+ //debug() << "sensorNumber = " << sensor << endmsg;
+
+
+ (*_cellid_encoder)[ lcio::ILDCellID0::subdet ] = lcio::ILDDetID::TPC ;
+ (*_cellid_encoder)[ lcio::ILDCellID0::layer ] = seed_hit->getRowIndex() ;
+ (*_cellid_encoder)[ lcio::ILDCellID0::module ] = 0 ;
+
+ //SJA:FIXME: for now don't use side
+ // (*_cellid_encoder)[ lcio::ILDCellID0::side ] = side ;
+ (*_cellid_encoder)[ lcio::ILDCellID0::side ] = lcio::ILDDetID::barrel ;
+ trkHit.setCellID(_cellid_encoder->lowWord());
+ //_cellid_encoder->setCellID( &trkHit ) ;
+
+
+ // check values for inf and nan
+ if( std::isnan(unsmearedPhi) || std::isinf(unsmearedPhi) || std::isnan(tpcRPhiRes) || std::isinf(tpcRPhiRes) ) {
+ std::stringstream errorMsg;
+ errorMsg << "\nProcessor: TPCDigiAlg \n"
+ << "unsmearedPhi = "
+ << unsmearedPhi
+ << " tpcRPhiRes = "
+ << tpcRPhiRes
+ << "\n" ;
+ throw errorMsg.str();
+ }
+ debug() << "==============" << endmsg;
+ // For no error in R
+ std::array<float,TRKHITNCOVMATRIX> covMat={sin(unsmearedPhi)*sin(unsmearedPhi)*tpcRPhiRes*tpcRPhiRes,
+ -cos(unsmearedPhi)*sin(unsmearedPhi)*tpcRPhiRes*tpcRPhiRes,
+ cos(unsmearedPhi)*cos(unsmearedPhi)*tpcRPhiRes*tpcRPhiRes,
+ 0.,
+ 0.,
+ float(tpcZRes*tpcZRes)};
+
+ trkHit.setCovMatrix(covMat);
+
+ if( !_tpcHitMap[seed_hit].isAvailable()){
+ std::stringstream errorMsg;
+ errorMsg << "\nProcessor: TPCDigiAlg \n"
+ << "SimTracker Pointer is NULL throwing exception\n"
+ << "\n" ;
+ throw errorMsg.str();
+ }
+ debug() << "==============" << endmsg;
+ if(pos[0]*pos[0]+pos[1]*pos[1]>0.0){
+ // push back the SimTHit for this TrackerHit
+
+ if (_use_raw_hits_to_store_simhit_pointer) {
+ trkHit.addToRawHits(_tpcHitMap[seed_hit].getObjectID());
+ }
+
+ auto rel = _relCol->create();
+ rel.setRec (trkHit);
+ rel.setSim (_tpcHitMap[seed_hit]);
+ rel.setWeight( 1.0 );
+
+ _trkhitVec->push_back( trkHit );
+ _NRecTPCHits++;
+ }
+
+ debug() << "==============" << endmsg;
+#ifdef DIGIPLOTS
+// edm4hep::SimTrackerHit* theSimHit = _tpcHitMap[seed_hit];
+// double rSimSqrd = theSimHit->getPosition()[0]*theSimHit->getPosition()[0] + theSimHit->getPosition()[1]*theSimHit->getPosition()[1];
+//
+// double phiSim = atan2(theSimHit->getPosition()[1],theSimHit->getPosition()[0]);
+//
+// double rPhiDiff = (point.phi() - phiSim)*sqrt(rSimSqrd);
+// double rPhiPull = ((point.phi() - phiSim)*sqrt(rSimSqrd))/(sqrt((covMat[2])/(cos(point.phi())*cos(point.phi()))));
+//
+// double zDiff = point.getZ() - theSimHit->getPosition()[2];
+// double zPull = zDiff/sqrt(covMat[5]);
+//
+//
+ //_rPhiDiffHisto->fill(rPhiDiff);
+ //_rPhiPullHisto->fill(rPhiPull);
+ //_phiDistHisto->fill(point.phi() - phiSim);
+ //_zDiffHisto->fill(zDiff);
+ //_zPullHisto->fill(zPull);
+
+ //_zSigmaVsZHisto->fill(seed_hit->getZ(),sqrt(covMat[5]));
+ //_rPhiSigmaHisto->fill(sqrt((covMat[2])/(cos(point.phi())*cos(point.phi()))));
+ //_zSigmaHisto->fill(sqrt(covMat[5]));
+#endif
+}
+
+void TPCDigiAlg::writeMergedVoxelsToHit( vector <Voxel_tpc*>* hitsToMerge){
+
+ const gear::TPCParameters& gearTPC = _GEAR->getTPCParameters() ;
+ const gear::PadRowLayout2D& padLayout = gearTPC.getPadLayout() ;
+ const gear::Vector2D padCoord = padLayout.getPadCenter(1) ;
+
+ edm4hep::TrackerHit trkHit;// = _trkhitVec->create();
+
+ double sumZ = 0;
+ double sumPhi = 0;
+ double sumEDep = 0;
+ // double R = 0;
+ double lastR = 0;
+
+ unsigned number_of_hits_to_merge = hitsToMerge->size();
+
+
+ for(unsigned int ihitCluster = 0; ihitCluster < number_of_hits_to_merge; ++ihitCluster){
+
+ sumZ += hitsToMerge->at(ihitCluster)->getZ();
+ sumPhi += hitsToMerge->at(ihitCluster)->getPhi();
+ sumEDep += hitsToMerge->at(ihitCluster)->getEDep();
+ hitsToMerge->at(ihitCluster)->setIsMerged();
+ lastR = hitsToMerge->at(ihitCluster)->getR();
+
+ if (_use_raw_hits_to_store_simhit_pointer) {
+ trkHit.addToRawHits(_tpcHitMap[hitsToMerge->at(ihitCluster)].getObjectID());
+ }
+
+ auto rel = _relCol->create();
+ rel.setRec (trkHit);
+ rel.setSim (_tpcHitMap[ hitsToMerge->at(ihitCluster) ]);
+ rel.setWeight( float(1.0/number_of_hits_to_merge) );
+
+ }
+
+ double avgZ = sumZ/(hitsToMerge->size());
+ double avgPhi = sumPhi/(hitsToMerge->size());
+
+ CLHEP::Hep3Vector* mergedPoint = new CLHEP::Hep3Vector(1.0,1.0,1.0);
+ mergedPoint->setPerp(lastR);
+ mergedPoint->setPhi(avgPhi);
+ mergedPoint->setZ(avgZ);
+
+ //store hit variables
+
+ // first the hit pos has to be smeared------------------------------------------------
+
+ //FIXME: which errors should we use for smearing the merged hits ?
+ // this might be a bit large ....
+ double tpcRPhiRes = _padWidth;
+ double tpcZRes = _binningZ;
+
+ CLHEP::Hep3Vector point( mergedPoint->getX(), mergedPoint->getY(), mergedPoint->getZ() ) ;
+
+// double unsmearedPhi = point.phi();
+
+ double randrp = gsl_ran_gaussian(_random,tpcRPhiRes);
+ double randz = gsl_ran_gaussian(_random,tpcZRes);
+
+ point.setPhi( point.phi() + randrp/ point.perp() );
+ point.setZ( point.z() + randz );
+
+ // make sure the hit is not smeared beyond the TPC Max DriftLength
+ if( fabs(point.z()) > gearTPC.getMaxDriftLength() ) point.setZ( (fabs(point.z()) / point.z() ) * gearTPC.getMaxDriftLength() );
+
+ double pos[3] = {point.x(),point.y(),point.z()};
+
+ //---------------------------------------------------------------------------------
+ trkHit.setPosition(pos);
+ trkHit.setEDep(sumEDep);
+ // trkHit->setType( 500 );
+
+ // SJA:FIXME: here you can use the value 2 but not 3 which is odd as the width of the field is 1, only 0 and 1 should be allowed?
+ int side = 1 ;
+ int padIndex = padLayout.getNearestPad(mergedPoint->perp(),mergedPoint->phi());
+ int row = padLayout.getRowNumber(padIndex);
+
+ if( pos[2] < 0.0 ) side = 1 ;
+
+ (*_cellid_encoder)[ lcio::ILDCellID0::subdet ] = lcio::ILDDetID::TPC ;
+ (*_cellid_encoder)[ lcio::ILDCellID0::layer ] = row ;
+ (*_cellid_encoder)[ lcio::ILDCellID0::module ] = 0 ;
+ ////SJA:FIXME: for now don't use side
+ //// (*_cellid_encoder)[ lcio::ILDCellID0::side ] = side ;
+ (*_cellid_encoder)[ lcio::ILDCellID0::side ] = 0 ;
+
+ //yzhang FIXME?
+ //_cellid_encoder->setCellID( &trkHit ) ;
+ trkHit.setCellID( _cellid_encoder->lowWord() );
+
+ double phi = mergedPoint->getPhi();
+
+ // check values for inf and nan
+ if( std::isnan(phi) || std::isinf(phi) || std::isnan(tpcRPhiRes) || std::isinf(tpcRPhiRes) ) {
+ std::stringstream errorMsg;
+ errorMsg << "\nProcessor: TPCDigiAlg \n"
+ << "phi = "
+ << phi
+ << " tpcRPhiRes = "
+ << tpcRPhiRes
+ << "\n" ;
+ throw errorMsg.str();
+ }
+
+ // For no error in R
+ std::array<float,TRKHITNCOVMATRIX> covMat={sin(phi)*sin(phi)*tpcRPhiRes*tpcRPhiRes,
+ -cos(phi)*sin(phi)*tpcRPhiRes*tpcRPhiRes,
+ cos(phi)*cos(phi)*tpcRPhiRes*tpcRPhiRes,
+ 0.,
+ 0.,
+ float(tpcZRes*tpcZRes)};
+
+ trkHit.setCovMatrix(covMat);
+
+ if(pos[0]*pos[0]+pos[1]*pos[1]>0.0){
+ //yzhang TODO
+ _trkhitVec->push_back( trkHit );
+ ++_nRechits;
+ } else {
+ //???yzhang TODO
+ //delete trkHit;
+ }
+
+ delete mergedPoint;
+
+}
+
+
+#ifdef DIGIPLOTS
+//void TPCDigiAlg::plotHelixHitResidual( edm4hep::MCParticle *mcp, CLHEP::Hep3Vector *thisPoint){
+//
+// const double bField = _GEAR->getBField().at( gear::Vector3D( 0., 0., 0.) ).z() ;
+// const double FCT = 2.99792458E-4;
+// double charge = mcp->getCharge();
+// const double *mom = mcp->getMomentum();
+// double pt = sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+// double radius = pt / (FCT*bField);
+// double tanLambda = mom[2]/pt;
+// double omega = charge / radius;
+//
+// if(pt>1.0) {
+//
+// //FIXME SJA: this is only valid for tracks from the IP and should be done correctly for non prompt tracks
+// double Z0 = 0.;
+// double D0 = 0.;
+//
+// LCVector3D refPoint(0.,0.,0);
+//
+// SimpleHelix* helix = new SimpleHelix(D0,
+// atan2(mom[1],mom[0]),
+// omega,
+// Z0,
+// tanLambda,
+// refPoint);
+//
+//
+// // an almost "infinite" cylinder in z
+// LCVector3D startCylinder(0.,0.,-1000000.0);
+// LCVector3D endCylinder(0.,0.,1000000.0);
+// bool endplane=true;
+//
+// LCCylinder cylinder(startCylinder,endCylinder,thisPoint->perp(),endplane);
+//
+// bool pointExists = false;
+//
+// double pathlength = helix->getIntersectionWithCylinder( cylinder, pointExists);
+//
+// LCErrorMatrix* errors = new LCErrorMatrix();
+//
+// if(pointExists){
+//
+// LCVector3D intersection = helix->getPosition(pathlength, errors);
+//
+// double intersectionPhi = atan2(intersection[1],intersection[0]);
+// double residualRPhi = ((intersectionPhi-thisPoint->phi())) ;
+// _ResidualsRPhiHisto->fill(residualRPhi);
+//
+// }
+//
+// delete errors;
+// delete helix;
+//
+// const gear::TPCParameters& gearTPC = _GEAR->getTPCParameters() ;
+// const gear::PadRowLayout2D& padLayout = gearTPC.getPadLayout() ;
+//
+// int row = padLayout.getRowNumber(padLayout.getNearestPad(thisPoint->perp(),thisPoint->phi()));
+// int pad = padLayout.getNearestPad(thisPoint->perp(),thisPoint->phi());
+//
+// double rHit_diff = thisPoint->perp()
+// - padLayout.getPlaneExtent()[0]
+// - (( row + 0.5 )
+// * padLayout.getPadHeight(pad)) ;
+//
+// _radiusCheckHisto->fill(rHit_diff);
+//
+// // info() << "$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$#$" << endmsg;
+// // info() << "thisPoint->perp() = " << thisPoint->perp() << endmsg;
+// // info() << "TPC Sensitive rMin = " << padLayout.getPlaneExtent()[0] << endmsg;
+// // info() << "Row number + 0.5 = " << row + 0.5 << endmsg;
+// // info() << "Pad Height = " << padLayout.getPadHeight(pad) << endmsg;
+// // info() << "Row Height = " << padLayout.getRowHeight(row) << endmsg;
+// // info() << "R_hit - TPC Rmin - ((RowIndex + 0.5 )* padheight) = " << rHit_diff << endmsg;
+//
+// }
+// return;
+//}
+#endif
+
+
+double TPCDigiAlg::getPadPhi( CLHEP::Hep3Vector *thisPoint, CLHEP::Hep3Vector* firstPoint, CLHEP::Hep3Vector* middlePoint, CLHEP::Hep3Vector* lastPoint){
+
+ CLHEP::Hep2Vector firstPointRPhi(firstPoint->x(),firstPoint->y());
+ CLHEP::Hep2Vector middlePointRPhi(middlePoint->x(),middlePoint->y());
+ CLHEP::Hep2Vector lastPointRPhi(lastPoint->x(),lastPoint->y());
+
+ // check that the points are not the same, at least at the level of a tenth of a micron
+ if( (fabs( firstPointRPhi.x() - middlePointRPhi.x() ) < 1.e-05 && fabs( firstPointRPhi.y() - middlePointRPhi.y() ) < 1.e-05)
+ ||
+ (fabs( middlePointRPhi.x() - lastPointRPhi.x() ) < 1.e-05 && fabs( middlePointRPhi.y() - lastPointRPhi.y() ) < 1.e-05)
+ ||
+ (fabs( firstPointRPhi.x() - lastPointRPhi.x() ) < 1.e-05 && fabs( firstPointRPhi.y() - lastPointRPhi.y() ) < 1.e-05)
+ ) {
+
+ warning() << " TPCDigiAlg::getPadPhi "
+ << "2 of the 3 SimTracker hits passed to Circle Fit are the same hit taking pad phi as PI/2\n"
+ << " firstPoint->x() " << firstPoint->x()
+ << " firstPoint->y() " << firstPoint->y()
+ << " firstPoint->z() " << firstPoint->z()
+ << " middlePoint->x() " << middlePoint->x()
+ << " middlePoint->y() " << middlePoint->y()
+ << " middlePoint->z() " << middlePoint->z()
+ << " lastPoint->x() " << lastPoint->x()
+ << " lastPoint->y() " << lastPoint->y()
+ << " lastPoint.z() " << lastPoint->z()
+ << std::endl ;
+
+ return CLHEP::twopi/4.0 ;
+
+ }
+
+
+
+ Circle theCircle(&firstPointRPhi, &middlePointRPhi, &lastPointRPhi);
+
+ double localPhi = atan2((thisPoint->y() - theCircle.GetCenter()->y()) ,(thisPoint->x() - theCircle.GetCenter()->x())) + (CLHEP::twopi/4.0) ;
+
+ if(localPhi>CLHEP::twopi) localPhi=localPhi - CLHEP::twopi;
+ if(localPhi<0.0) localPhi=localPhi + CLHEP::twopi;
+ if(localPhi>CLHEP::pi) localPhi = localPhi - CLHEP::pi ;
+
+ double pointPhi = thisPoint->phi();
+
+ if(pointPhi>CLHEP::twopi) pointPhi=pointPhi - CLHEP::twopi;
+ if(pointPhi<0.0) pointPhi=pointPhi + CLHEP::twopi;
+ if(pointPhi>CLHEP::pi) pointPhi = pointPhi - CLHEP::pi;
+
+ double padPhi = fabs(pointPhi - localPhi);
+
+ // check that the value returned is reasonable
+ if( std::isnan(padPhi) || std::isinf(padPhi) ) {
+ std::stringstream errorMsg;
+ errorMsg << "\nProcessor: TPCDigiAlg \n"
+ << "padPhi = "
+ << padPhi
+ << "\n" ;
+ throw errorMsg.str();
+ }
+
+ return padPhi;
+
+}
+
+double TPCDigiAlg::getPadTheta(CLHEP::Hep3Vector* firstPoint, CLHEP::Hep3Vector* middlePoint, CLHEP::Hep3Vector* lastPoint){
+
+ // Calculate thetaPad for current hit
+ CLHEP::Hep2Vector firstPointRPhi(firstPoint->x(),firstPoint->y()) ;
+ CLHEP::Hep2Vector middlePointRPhi(middlePoint->x(),middlePoint->y());
+ CLHEP::Hep2Vector lastPointRPhi(lastPoint->x(),lastPoint->y());
+
+ // check that the points are not the same, at least at the level of a tenth of a micron
+ if( (fabs( firstPointRPhi.x() - middlePointRPhi.x() ) < 1.e-05 && fabs( firstPointRPhi.y() - middlePointRPhi.y() ) < 1.e-05)
+ ||
+ (fabs( middlePointRPhi.x() - lastPointRPhi.x() ) < 1.e-05 && fabs( middlePointRPhi.y() - lastPointRPhi.y() ) < 1.e-05)
+ ||
+ (fabs( firstPointRPhi.x() - lastPointRPhi.x() ) < 1.e-05 && fabs( firstPointRPhi.y() - lastPointRPhi.y() ) < 1.e-05)
+ ) {
+
+ warning() << " TPCDigiAlg::getPadTheta "
+ << "2 of the 3 SimTracker hits passed to Circle Fit are the same hit taking pad phi as PI/2\n"
+ << " firstPoint->x() " << firstPoint->x()
+ << " firstPoint->y() " << firstPoint->y()
+ << " firstPoint->z() " << firstPoint->z()
+ << " middlePoint->x() " << middlePoint->x()
+ << " middlePoint->y() " << middlePoint->y()
+ << " middlePoint->z() " << middlePoint->z()
+ << " lastPoint->x() " << lastPoint->x()
+ << " lastPoint->y() " << lastPoint->y()
+ << " lastPoint.z() " << lastPoint->z()
+ << endmsg;
+
+ return CLHEP::twopi/4.0 ;
+
+ }
+
+
+ Circle theCircle(&firstPointRPhi, &middlePointRPhi, &lastPointRPhi);
+
+ double deltaPhi = firstPoint->deltaPhi(*lastPoint);
+
+ double pathlength = fabs(deltaPhi) * theCircle.GetRadius();
+
+ double padTheta = atan ( pathlength / fabs(lastPoint->z() - firstPoint->z()) ) ;
+
+ double pathlength1 = 2.0 * asin( ( sqrt (
+ (middlePointRPhi.x() - firstPointRPhi.x()) * (middlePointRPhi.x()-firstPointRPhi.x())
+ +
+ (middlePointRPhi.y()-firstPointRPhi.y()) * (middlePointRPhi.y()-firstPointRPhi.y())
+ ) / 2.0 ) / theCircle.GetRadius() ) * theCircle.GetRadius() ;
+
+
+ double pathlength2 = 2.0 * asin( ( sqrt (
+ (lastPointRPhi.x()-middlePointRPhi.x()) * (lastPointRPhi.x()-middlePointRPhi.x())
+ +
+ (lastPointRPhi.y()-middlePointRPhi.y()) * (lastPointRPhi.y()-middlePointRPhi.y())
+ ) / 2.0 ) / theCircle.GetRadius() ) * theCircle.GetRadius() ;
+
+
+ padTheta = atan ((fabs(pathlength1 + pathlength2)) / (fabs(lastPoint->z() - firstPoint->z())) ) ;
+
+ // check that the value returned is reasonable
+ if( std::isnan(padTheta) || std::isinf(padTheta) ) {
+ std::stringstream errorMsg;
+ errorMsg << "\nProcessor: TPCDigiAlg \n"
+ << "padTheta = "
+ << padTheta
+ << "\n" ;
+ throw errorMsg.str();
+ }
+
+ return padTheta;
+
+}
diff --git a/Digitisers/SimpleDigi/src/TPCDigiAlg.h b/Digitisers/SimpleDigi/src/TPCDigiAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..a1dc6019602753a3b197a692413efe7b82f374f1
--- /dev/null
+++ b/Digitisers/SimpleDigi/src/TPCDigiAlg.h
@@ -0,0 +1,271 @@
+/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+
+/*
+Evolved version of TPCDigi that provides additional functionality to deal with background. Couple to the Mokka Sensitive Detector Driver TPCSD03.cc
+
+SJA:FIXME: Still needs to be tidied up for production release.
+
+Three cases can be consider in the treatment of SimTrackerHits
+i) A clean isolated hit; this will be smeared according to the parametric point resolution and converted to a TrackerHit
+ii) Two or Three hits which are considered to be closer than the double hit resolution and which therefore cannot be viewed as seperable hits. These will be merged and be assigned a large associated measurement error.
+iii) A continuous set of hits within one pad row which cannot be resolved as single hits, these are condidered to be charaterisable as background hits created by extremely low pt charged particles (pt < 10MeV) and therefore are removed from the hit collection.
+
+The Driver has been modified to take an additional collection of SimTrackerHits which are produced by the Mokka TPC Sensitive Driver TPCSD03.cc. These hits are produced for particles which have very low pt and often do not move outside of the dimensions of a single pad row. These hits need to be treated differently as they do not cross any geometric boundaries in a Padrow based TPC Geometry. This negates the need to voxalise the TPC in Geant4 which has proved in the past to be prohibitive in terms of processing time due to the vastly increased number of geometric volumes.
+
+Steve Aplin 26 June 2009 (DESY)
+
+*/
+
+#ifndef TPCDigiAlg_h
+#define TPCDigiAlg_h 1
+
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "FWCore/DataHandle.h"
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+
+#include <gsl/gsl_rng.h>
+
+//#ifdef MARLIN_USE_AIDA
+
+//#include <marlin/AIDAProcessor.h>
+//#include <AIDA/IHistogramFactory.h>
+//#include <AIDA/ICloud1D.h>
+//#include <AIDA/IHistogram1D.h>
+
+
+
+//#define DIGIPLOTS
+
+
+
+//#ifdef DIGIPLOTS
+// includes all AIDA header files
+//#include <AIDA/AIDA.h>
+//#endif
+
+//#endif
+
+#include <string>
+#include <vector>
+#include <map>
+
+//#include <EVENT/LCCollection.h>
+//#include <EVENT/SimTrackerHit.h>
+//#include <EVENT/MCParticle.h>
+//#include <IMPL/LCCollectionVec.h>
+//#include <IMPL/TrackerHitImpl.h>
+#include <UTIL/BitField64.h>
+#include <edm4hep/MCParticle.h>
+#include "edm4hep/SimTrackerHit.h"
+
+#include "CLHEP/Vector/TwoVector.h"
+class Voxel_tpc;
+
+
+
+
+//using namespace lcio ;
+//using namespace marlin ;
+//#ifdef MARLIN_USE_AIDA
+//using namespace AIDA ;
+//#endif
+
+
+namespace gear { class GearMgr; }
+class IEventSeeder;
+
+/** ====== TPCDigiAlg ====== <br>
+ *
+ * This algorithm depends on Circle.h from MarlinUtil
+ *
+ * Caution: This digitiser presently does not process space-point like SimTrackerHits which have been flagged with CellIDs set to the negetive row number. This must be implemented in future.
+ *Produces TPC TrackerHit collection from SimTrackerHit collection, smeared in r-phi and z.
+ * Double hits are identified but are currently not added to the collection. This may be change
+ * at a later date when criteria for their seperation is defined. The resolutions are defined in
+ * the GEAR stearing file.
+ *
+ * Resolution in r-phi is calculated according to the formular <br>
+ * sigma_{point}^2 = sigma_0^2 + Cd^2/N_{eff} * L_{drift}
+ * Cd^2/N_{eff}} = 25^2/(22/sin(theta)*h/6mm)
+ * Cd = 25 ( microns / cm^(1/2) )
+ * (this is for B=4T, h is the pad height = pad-row pitch in mm,
+ * theta is the polar angle)
+ *
+ * At the moment resolution in z assumed to be independent of drift length. <br>
+ *
+ * The type of TPC TrackerHit is set to 500 via method TrackerHitImpl::setType(int type) <br>
+ * <h4>Input collections and prerequisites</h4>
+ * Algorithm requires collections of SimTrackerHits in TPC <br>
+ * <h4>Output</h4>
+ * Processor produces collection of digitized TrackerHits in TPC <br>
+ * @param CollectionName The name of input SimTrackerHit collection <br>
+ * (default name STpc01_TPC)
+ * @param RejectCellID0 Whether or not to reject SimTrackerHits with Cell ID 0. Mokka drivers
+ * TPC00-TPC03 encoded the pad row number in the cell ID, which should always be non-zero anyway.
+ * Drivers TPC04 and TPC05 do not simulate pad rows and thus have the cell ID set to zero for all hits.
+ * You will need to set RejectCellID0 to 0 in order to use this processor with these drivers, but note
+ * that the implications for track reconstruction are not strictly defined. Mokka driver TPC06 uses
+ * a mixed approach with one hit per pad row having non-zero cell ID, extra hits having 0 cell ID.
+ * Typically, unless you use TPC04 or TPC05, you should not touch this parameter. <br>
+ * (default value 1)
+ * @param TPCTrackerHitsCol The name of output collection of TrackerHits <br>
+ * (default name TPCTrackerHits) <br>
+ * <br>
+ * @authors S. Aplin, DESY and A.Raspereza, MPI
+ *
+ * Changed 7/9/07 so that the const and diffusion resolution terms are taken as processor parameters rather than the gear file.
+ * The parameters _pixZ and pixRP were also changed from gear parameters to processor parameters
+ * clare.lynch@bristol.ac.uk
+ *
+ */
+class TPCDigiAlg: public GaudiAlgorithm{
+
+public:
+
+ TPCDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+
+ /** Called at the begin of the job before anything is read.
+ * Use to initialize the processor, e.g. book histograms.
+ */
+ virtual StatusCode initialize();
+
+ /** Called for every event - the working horse.
+ */
+ virtual StatusCode execute();
+
+ /** Called after data processing for clean up.
+ */
+ virtual StatusCode finalize();
+
+ void writeVoxelToHit( Voxel_tpc* aVoxel) ;
+ void writeMergedVoxelsToHit( std::vector <Voxel_tpc*>* hitList ) ;
+ void plotHelixHitResidual(edm4hep::MCParticle *mcp, CLHEP::Hep3Vector *thisPointRPhi);
+ double getPadPhi( CLHEP::Hep3Vector* thisPointRPhi, CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi);
+ double getPadTheta( CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi );
+
+protected:
+
+ DataHandle<edm4hep::EventHeaderCollection> _headerCol{"EventHeader", Gaudi::DataHandle::Reader, this};
+ /** Input collection name.
+ */
+ //std::string _padRowHitColName ;
+ //std::string _spacePointColName ;
+ //std::string _lowPtHitscolName ;
+ DataHandle<edm4hep::SimTrackerHitCollection> _padRowHitColHdl{"TPCCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> _spacePointColHdl{"TPCSpacePointCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> _lowPtHitsColHdl{"TPCLowPtCollection", Gaudi::DataHandle::Reader, this};
+ //DataHandle<edm4hep::MCParticleCollection> _mcColHdl{"MCParticle", Gaudi::DataHandle::Reader, this};
+
+ /** Output collection name.
+ */
+ DataHandle<edm4hep::TrackerHitCollection> _TPCTrackerHitsColHdl{"TPCTrackerHits", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection> _TPCAssColHdl{"TPCTrackerHitAss", Gaudi::DataHandle::Writer, this};
+ edm4hep::TrackerHitCollection* _trkhitVec;
+ edm4hep::MCRecoTrackerAssociationCollection* _relCol;
+
+ bool _use_raw_hits_to_store_simhit_pointer;
+
+ int _rejectCellID0;
+ float _padWidth;
+
+ int _nRun ;
+ int _nEvt ;
+
+ //edm4hep::ConstMCParticle _mcp;
+ //edm4hep::ConstMCParticle _previousMCP;
+ //edm4hep::ConstMCParticle _nextMCP;
+ //edm4hep::ConstMCParticle _nMinus2MCP;
+ //edm4hep::ConstMCParticle _nPlus2MCP;
+
+ //edm4hep::SimTrackerHit _SimTHit;
+ //edm4hep::SimTrackerHit _previousSimTHit;
+ //edm4hep::SimTrackerHit _nextSimTHit;
+ //edm4hep::SimTrackerHit _nPlus2SimHit;
+ //edm4hep::SimTrackerHit _nMinus2SimHit;
+
+ gear::GearMgr* _GEAR;
+ IEventSeeder* _SEEDER;
+ // gsl random number generator
+ gsl_rng * _random ;
+
+ float _pointResoRPhi0; // Coefficient for RPhi point res independant of drift length
+ float _pointResoPadPhi; // Coefficient for the point res dependance on relative phi angle to the pad verticle
+ float _diffRPhi; // Coefficient for the rphi point res dependance on diffusion
+ int _nEff; // number of effective electrons
+
+
+ float _pointResoZ0; // Coefficient Z point res independant of drift length
+ float _diffZ; // Coefficient for the Z point res dependance on diffusion
+
+ float _binningZ;
+ float _binningRPhi;
+ float _doubleHitResZ;
+ float _doubleHitResRPhi;
+ int _maxMerge;
+
+ int _nRechits;
+
+ std::vector< std::vector <Voxel_tpc *> > _tpcRowHits;
+ std::map< Voxel_tpc *,edm4hep::SimTrackerHit > _tpcHitMap;
+
+ UTIL::BitField64* _cellid_encoder;
+
+ int _NSimTPCHits;
+ int _NBackgroundSimTPCHits;
+ int _NPhysicsSimTPCHits;
+ int _NPhysicsAbove02GeVSimTPCHits;
+ int _NPhysicsAbove1GeVSimTPCHits;
+ int _NRecTPCHits;
+
+ int _NLostPhysicsTPCHits;
+ int _NLostPhysicsAbove02GeVPtTPCHits;
+ int _NLostPhysicsAbove1GeVPtTPCHits;
+ int _NRevomedHits;
+
+
+//#ifdef DIGIPLOTS
+// IAnalysisFactory * _AF;
+// ITreeFactory * _TRF;
+// ITree * _TREE;
+// IHistogramFactory * _HF;
+// IHistogram1D * _phiDiffHisto;
+// IHistogram1D * _thetaDiffHisto;
+// IHistogram1D * _phiRelHisto;
+// IHistogram1D * _thetaRelHisto;
+//
+// IHistogram1D * _phiDistHisto;
+// IHistogram1D * _rPhiPullHisto;
+// IHistogram1D * _rPhiDiffHisto;
+// IHistogram1D * _zDiffHisto;
+// IHistogram1D * _zPullHisto;
+// IHistogram2D * _zSigmaVsZHisto;
+// IHistogram1D * _zSigmaHisto;
+// IHistogram1D * _rPhiSigmaHisto;
+// IHistogram1D * _radiusCheckHisto;
+// IHistogram1D * _ResidualsRPhiHisto;
+//
+// IHistogram1D * _NSimTPCHitsHisto;
+// IHistogram1D * _NBackgroundSimTPCHitsHisto;
+// IHistogram1D * _NPhysicsSimTPCHitsHisto;
+// IHistogram1D * _NPhysicsAbove02GeVSimTPCHitsHisto;
+// IHistogram1D * _NPhysicsAbove1GeVSimTPCHitsHisto;
+// IHistogram1D * _NRecTPCHitsHisto;
+//
+// IHistogram1D * _NLostPhysicsTPCHitsHisto;
+// IHistogram1D * _NLostPhysicsAbove02GeVPtTPCHitsHisto;
+// IHistogram1D * _NLostPhysicsAbove1GeVPtTPCHitsHisto;
+// IHistogram1D * _NRevomedHitsHisto;
+//
+// IHistogram1D * _NKeptPhysicsTPCHitsHistoPercent;
+// IHistogram1D * _NKeptPhysicsAbove02GeVPtTPCHitsHistoPercent;
+// IHistogram1D * _NKeptPhysicsAbove1GeVPtTPCHitsHistoPercent;
+//
+//#endif
+
+
+} ;
+#endif
diff --git a/Digitisers/SimpleDigi/src/voxel.cpp b/Digitisers/SimpleDigi/src/voxel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bfa801ac16dce7e5199f7b14fdba56e060b1807b
--- /dev/null
+++ b/Digitisers/SimpleDigi/src/voxel.cpp
@@ -0,0 +1,60 @@
+
+
+#include <iostream>
+#include "voxel.h"
+
+using namespace std;
+
+Voxel_tpc::Voxel_tpc(){
+}
+Voxel_tpc::Voxel_tpc(int row, int phi, int z, double pos[3], double posRPhi[2], double edep, double RPhiRes, double ZRes)
+{
+ _row_index = row;
+ _phi_index = phi;
+ _z_index = z;
+ _coord.setX(pos[0]);
+ _coord.setY(pos[1]);
+ _coord.setZ(pos[2]);
+ _edep = edep;
+ _rPhiRes = RPhiRes;
+ _zRes = ZRes;
+ _isMerged = false;
+ _isClusterHit = false;
+}
+
+Voxel_tpc::Voxel_tpc(int row, int phi, int z, CLHEP::Hep3Vector coord, double edep, double RPhiRes, double ZRes)
+{
+ _row_index = row;
+ _phi_index = phi;
+ _z_index = z;
+ _coord=coord;
+ _edep = edep;
+ _rPhiRes = RPhiRes;
+ _zRes = ZRes;
+ _isMerged = false;
+ _isClusterHit = false;
+}
+
+Voxel_tpc::~Voxel_tpc()
+{
+}
+
+//bool Voxel_tpc::compare_phi( Voxel_tpc * & a, Voxel_tpc * & b){
+
+// return ( a->getPhiIndex() < b->getPhiIndex());
+
+//}
+
+int Voxel_tpc::clusterFind(vector <Voxel_tpc*>* hitList){
+
+ if(!this->IsClusterHit()){
+ hitList->push_back(this);
+ this->setIsClusterHit();
+ for(int i=0; i<this->getNumberOfAdjacent();++i){
+ getAdjacent(i)->clusterFind(hitList);
+ }
+ }
+
+ return hitList->size();
+}
+
diff --git a/Digitisers/SimpleDigi/src/voxel.h b/Digitisers/SimpleDigi/src/voxel.h
new file mode 100644
index 0000000000000000000000000000000000000000..fe3b2ab68bf2bcdc27ba1b70a73900d4cc814dfc
--- /dev/null
+++ b/Digitisers/SimpleDigi/src/voxel.h
@@ -0,0 +1,61 @@
+#ifndef _voxel_included_
+#define _voxel_included_
+
+// A header file which defines a voxel class for the TPC
+#include <vector>
+//#include "ThreeVector.h"
+#include <CLHEP/Vector/ThreeVector.h>
+
+using namespace std;
+
+class Voxel_tpc{
+
+ public:
+ Voxel_tpc();
+ // the intialation in the constructor here would be preferable though I don't know how to intialise
+ // the array xyz[3] here with pos[3], for the mean time the constructor will be put in the .cc file
+ // Voxel_tpc(int row, int phi, int z, double pos[3]) : row_index(row), phi_index(phi), z_index(z){}
+ Voxel_tpc(int row, int phi, int z, double pos[3], double posRPhi[2], double edep, double rPhiRes, double zRes);
+ Voxel_tpc(int row, int phi, int z, CLHEP::Hep3Vector coord, double edep, double rPhiRes, double zRes);
+ ~Voxel_tpc();
+
+ void setAdjacent(Voxel_tpc * p_voxel) { _adjacent_voxels.push_back(p_voxel);};
+ void setIsClusterHit() { _isClusterHit = true;};
+ void setIsMerged() { _isMerged = true;};
+ bool IsClusterHit() { return _isClusterHit;};
+ bool IsMerged() { return _isMerged;};
+ int clusterFind(vector <Voxel_tpc*>* hitList);
+
+
+ int getRowIndex() {return _row_index;};
+ int getPhiIndex() {return _phi_index;};
+ int getZIndex() {return _z_index;};
+ Voxel_tpc * getFirstAdjacent() {return *(_adjacent_voxels.begin());};
+ Voxel_tpc * getAdjacent(int i) {return _adjacent_voxels[i];};
+ int getNumberOfAdjacent() {return _adjacent_voxels.size();};
+ double getX() {return _coord.x();};
+ double getY() {return _coord.y();};
+ double getZ() {return _coord.z();};
+ double getR() {return _coord.perp();};
+ double getPhi() {return _coord.phi();};
+ double getEDep() {return _edep;};
+ double getRPhiRes() {return _rPhiRes;};
+ double getZRes() {return _zRes;};
+ const CLHEP::Hep3Vector getHep3Vector() {return _coord;};
+ // bool compare_phi( Voxel_tpc * & a, Voxel_tpc * & b);
+
+
+
+ private:
+ int _row_index;
+ int _phi_index;
+ int _z_index;
+ vector <Voxel_tpc *> _adjacent_voxels;
+ CLHEP::Hep3Vector _coord;
+ double _edep;
+ double _rPhiRes;
+ double _zRes;
+ bool _isMerged;
+ bool _isClusterHit;
+};
+#endif
diff --git a/Utilities/DataHelper/CMakeLists.txt b/Utilities/DataHelper/CMakeLists.txt
index a0b346f3580fbbf1567c01dda9710193118d43a0..a93dc10a07a2ab5cd239c4dea1cd1ef23565cecf 100644
--- a/Utilities/DataHelper/CMakeLists.txt
+++ b/Utilities/DataHelper/CMakeLists.txt
@@ -1,13 +1,12 @@
gaudi_subdir(DataHelper v0r0)
+find_package(CLHEP REQUIRED;CONFIG)
find_package(EDM4HEP REQUIRED)
find_package(GSL REQUIRED )
message("GSL: ${GSL_LIBRARIES} ")
message("GSL INCLUDE_DIRS: ${GSL_INCLUDE_DIRS} ")
-INCLUDE_DIRECTORIES(${GSL_INCLUDE_DIRS})
-
-gaudi_depends_on_subdirs()
+# gaudi_depends_on_subdirs()
set(DataHelperLib_srcs src/*.cc src/*.cpp)
@@ -15,5 +14,6 @@ set(DataHelperLib_srcs src/*.cc src/*.cpp)
gaudi_add_library(DataHelperLib ${DataHelperLib_srcs}
PUBLIC_HEADERS DataHelper
- LINK_LIBRARIES EDM4HEP::edm4hep EDM4HEP::edm4hepDict ${GSL_LIBRARIES}
+ INCLUDE_DIRS ${CLHEP_INCLUDE_DIR} ${GSL_INCLUDE_DIRS}
+ LINK_LIBRARIES EDM4HEP::edm4hep EDM4HEP::edm4hepDict ${GSL_LIBRARIES} ${CLHEP_LIBRARIES}
)
diff --git a/Utilities/DataHelper/DataHelper/Circle.h b/Utilities/DataHelper/DataHelper/Circle.h
new file mode 100644
index 0000000000000000000000000000000000000000..f380976cd05a76b485dc0f584242e7f46d29851a
--- /dev/null
+++ b/Utilities/DataHelper/DataHelper/Circle.h
@@ -0,0 +1,31 @@
+// Circle.h: interface for the Circle class.
+// Circle class.
+// Purpose : Represent the circle object
+// Input : 3 different points
+// Process : Calcuate the radius and center
+// Output : Circle
+//
+// This class originally designed for representation of discretized curvature information
+// of sequential pointlist
+// KJIST CAD/CAM Ryu, Jae Hun ( ryu@geguri.kjist.ac.kr)
+// Last update : 1999. 7. 4
+
+
+#include "CLHEP/Vector/TwoVector.h"
+
+class Circle
+{
+public:
+ double GetRadius();
+ CLHEP::Hep2Vector* GetCenter();
+ Circle(CLHEP::Hep2Vector *p1, CLHEP::Hep2Vector *p2, CLHEP::Hep2Vector *p3); // p1, p2, p3 are co-planar
+ Circle();
+ virtual ~Circle();
+
+private:
+ double CalcCircle(CLHEP::Hep2Vector *pt1, CLHEP::Hep2Vector *pt2, CLHEP::Hep2Vector *pt3);
+ bool IsPerpendicular(CLHEP::Hep2Vector *pt1, CLHEP::Hep2Vector *pt2, CLHEP::Hep2Vector *pt3);
+ double m_dRadius;
+ CLHEP::Hep2Vector m_Center;
+};
+
diff --git a/Utilities/DataHelper/DataHelper/LCCylinder.h b/Utilities/DataHelper/DataHelper/LCCylinder.h
new file mode 100644
index 0000000000000000000000000000000000000000..924a55922004b9fe7373231a003366fa78d41ede
--- /dev/null
+++ b/Utilities/DataHelper/DataHelper/LCCylinder.h
@@ -0,0 +1,113 @@
+#ifndef LCCylinder_H
+#define LCCylinder_H 1
+
+// #include "CLHEP/Vector/ThreeVector.h"
+#include <DataHelper/LCGeometryTypes.h>
+
+/** Definition of a LCCylinder describing a geometrical cylinder in 3D space.
+ * @author T.Kraemer, DESY
+ * @version $Id: LCCylinder.h,v 1.1 2006-10-16 15:23:32 tkraemer Exp $
+ */
+
+class LCCylinder {
+
+public:
+
+ /**
+ * Constructor from two points and a radius.
+ * @param point1 point1 is the start point of the cylinder axis
+ * @param point2 point2 is the end point of the cylinder axis
+ * @param radius radius is the radius of the xylinder
+ * @param endPlane endPlane switches if cylinder is open or not
+ */
+ LCCylinder(LCVector3D point1,
+ LCVector3D point2,
+ double radius,
+ bool endPlane = false) ;
+
+ /**
+ * Constructor from one point, th Axis of the cylinder and the radius
+ * @param radius radius of cylinder.
+ * @param point point in the middle of the axis of the xylinder
+ * @param axis axis is the orientation of the xylinder axis.
+ * the length of axis is the half length of the cylinder
+ * @param endPlane endPlane switches if cylinder is open or not
+ */
+ LCCylinder(double radius, LCVector3D point, LCVector3D axis, bool endPlane) ;
+
+ /** Copy constructor.
+ * @param cylinder cylinder is an other LCCylinder.
+ */
+ LCCylinder(const LCCylinder & cylinder) ;
+
+ /**
+ * Destructor. */
+ ~LCCylinder() {}
+
+ /**
+ * Assignment. */
+ LCCylinder & operator=(const LCCylinder & rhs) ;
+
+ /**
+ * startpoint of cylinder axis */
+ LCVector3D startPoint() const ;
+
+ /**
+ * end point of cylinder axis */
+ LCVector3D endPoint() const ;
+
+ /**
+ * orientation of cylinder axis. the return vector is normalised */
+ LCVector3D axisDirection() const ;
+
+ /**
+ * length of cylinder axis */
+ double length() const ;
+
+ /**
+ * Radius of cylinder. */
+ double radius() const ;
+
+ /**
+ * Distance of a point to the cylinder.
+ * @param point point is a point in space
+ */
+ double distance(const LCVector3D & point) const ;
+
+ /**
+ * Projection of a point on to the surface of the cylinder.
+ * @param point point is a point in space.
+ * @param code code gives an integer code for the type of area the
+ * point gets projected to:
+ * 0 : No projection possible (point is not normal above the surface)
+ * 1 : prjection hits plane at start point ( startPoint() )
+ * 2 : prjection hits plane at end point ( endPoint() )
+ * 3 : projection hits the finite tube of the cylinder
+ */
+ LCVector3D projectPoint(const LCVector3D & point, int & code) const ;
+
+ /**
+ * Checks if a given point is inside the clyinder.
+ * @param point point is a point in space.
+ */
+ bool isInside(const LCVector3D & point) const ;
+
+ /**
+ * Test for equality. */
+ bool operator==(const LCCylinder & rhs) const ;
+
+ /**
+ * Test for inequality. */
+ bool operator!=(const LCCylinder & rhs) const ;
+
+protected:
+
+ double _radius;
+ bool _endPlane;
+
+ LCVector3D _axisSstartPoint;
+ LCVector3D _axisEndPoint;
+
+};
+
+#endif /* ifndef LCCylinder_H */
diff --git a/Utilities/DataHelper/DataHelper/LCGeometryTypes.h b/Utilities/DataHelper/DataHelper/LCGeometryTypes.h
new file mode 100644
index 0000000000000000000000000000000000000000..22a39b7c159b2574961264a817fa0658edfd100d
--- /dev/null
+++ b/Utilities/DataHelper/DataHelper/LCGeometryTypes.h
@@ -0,0 +1,38 @@
+#ifndef LCGeometryTypes_H
+#define LCGeometryTypes_H 1
+
+// #include "CLHEP/Geometry/Point3D.h"
+// #include "CLHEP/Geometry/Vector3D.h"
+#include "CLHEP/Geometry/Plane3D.h"
+#include "CLHEP/Matrix/SymMatrix.h"
+
+#include "CLHEP/Vector/ThreeVector.h"
+
+
+/** @file Definition of geometry types used in ILC software - currently use CHEP.
+ * @author gaede
+ * @version $Id: LCGeometryTypes.h,v 1.3 2006-10-11 16:03:24 tkraemer Exp $
+ */
+
+//using namespace CLHEP ;
+
+
+// typedef HepGeom::Point3D<double> LCPoint3D ;
+// typedef HepGeom::Vector3D<double> LCVector3D ;
+
+// typedef CLHEP::Hep3Vector LCPoint3D ;
+typedef CLHEP::Hep3Vector LCVector3D ;
+
+
+typedef CLHEP::HepSymMatrix LCErrorMatrix ;
+
+// typedef HepGeom::Plane3D<double> LCPlane3D ;
+
+typedef CLHEP::HepLorentzVector LCLorentzVector ;
+
+
+
+
+
+
+#endif /* ifndef LCGeometryTypes_H */
diff --git a/Utilities/DataHelper/DataHelper/LCLine3D.h b/Utilities/DataHelper/DataHelper/LCLine3D.h
new file mode 100644
index 0000000000000000000000000000000000000000..8ef1afee10c42b980eab1b98cdc8b7ebb88b2b33
--- /dev/null
+++ b/Utilities/DataHelper/DataHelper/LCLine3D.h
@@ -0,0 +1,138 @@
+#ifndef LCLine3D_H
+#define LCLine3D_H 1
+
+// #include "CLHEP/Vector/ThreeVector.h"
+#include <DataHelper/LCGeometryTypes.h>
+#include <DataHelper/LCPlane3D.h>
+
+/** Definition of a LCLine3D describing a geometrical line in 3D space.
+ * @author T.Kraemer, DESY
+ * @version $Id: LCLine3D.h,v 1.8 2006-11-03 16:43:13 tkraemer Exp $
+ */
+
+class LCLine3D {
+
+public:
+
+ /** Standard constructor:
+ * Initializes a line along the x-axis.
+ */
+ LCLine3D();
+
+ /**
+ * Constructor from a point and a direction.
+ * @param point Point is a point of the line
+ * @param direction Direction is the directional vector of the line.
+ */
+ LCLine3D(const LCVector3D & point, const LCVector3D & direction) ;
+
+ /**
+ * Constructor from a point and a direction.
+ * @param point Point is a point of the line
+ * @param direction Direction is the directional vector of the line.
+ * @param reference reference point of the line.
+ */
+ LCLine3D(const LCVector3D & point,
+ const LCVector3D & direction,
+ const LCVector3D & reference) ;
+
+ /**
+ * Constructor using the canonical parameterization.
+ * @param d0 d0 is the point of closest approach in the xy plane.
+ * @param phi0 phi0 is the angle in the xy plane.
+ * @param z0 z0 is the z coordinate of the point of closest approach.
+ * @param tanLambda tanLambda is the angle of with respect to the xy plane.
+ */
+ LCLine3D(double d0, double phi0, double z0, double tanLambda) ;
+
+ /**
+ * Constructor using the canonical parameterization.
+ * @param d0 d0 is the point of closest approach in the xy plane.
+ * @param phi0 phi0 is the angle in the xy plane.
+ * @param z0 z0 is the z coordinate of the point of closest approach.
+ * @param tanLambda tanLambda is the angle of with respect to the xy plane.
+ * @param reference reference point of the line.
+ */
+ LCLine3D(double d0, double phi0, double z0, double tanLambda,
+ const LCVector3D & reference) ;
+
+ /** Copy constructor.
+ * @param line line is an other LCLine3D.
+ */
+ LCLine3D(const LCLine3D & line) ;
+
+ /**
+ * Destructor. */
+ ~LCLine3D() {}
+
+ /**
+ * set the Parameters for a line using a point and a direction.
+ * @param point Point is a point of the line
+ * @param direction Direction is the directional vector of the line.
+ * @param reference reference point of the line.
+ */
+ bool set(const LCVector3D & point,
+ const LCVector3D & direction,
+ const LCVector3D & reference) ;
+
+ /**
+ * Set the Parameters of a line using the canonical parameterization.
+ * @param d0 d0 is the point of closest approach in the xy plane.
+ * @param phi0 phi0 is the angle in the xy plane.
+ * @param z0 z0 is the z coordinate of the point of closest approach.
+ * @param tanLambda tanLambda is the angle of with respect to the xy plane.
+ * @param reference reference point of the line.
+ */
+ bool set(double d0, double phi0, double z0, double tanLambda,
+ const LCVector3D & reference) ;
+
+ /**
+ * Assignment. */
+ LCLine3D & operator=(const LCLine3D & rhs) ;
+
+ /**
+ * Position is the point of the line after a distance s.
+ * Is is given with respect to the point of closes approach to the origen of
+ * the coordinate system.
+ * @param s s is the path length along the line */
+ LCVector3D position(const double s = 0) const ;
+
+ /** Direction of the line
+ */
+ LCVector3D direction() const ;
+
+ /**
+ * Distance of a point to the line.
+ * @param point point is a point in space
+ */
+ double distance(const LCVector3D & point) const ;
+
+ /**
+ * Projection of a point on to the line.
+ * @param point point is a point in space.
+ */
+ double projectPoint(const LCVector3D & point) const ;
+
+ /**
+ * Test for equality. */
+ bool operator==(const LCLine3D & rhs) const ;
+
+ /**
+ * Test for inequality. */
+ bool operator!=(const LCLine3D & rhs) const ;
+
+ /** Pathlength at closest intersection point with plane - undefined
+ * if pointExists==false.
+ */
+ double intersectionWithPlane(const LCPlane3D plane, bool& pointExists) const ;
+
+protected:
+
+ LCVector3D _point;
+ LCVector3D _direction;
+ LCVector3D _reference;
+};
+
+std::ostream & operator << (std::ostream &os, const LCLine3D &l) ;
+
+#endif /* ifndef LCLine3D_H */
diff --git a/Utilities/DataHelper/DataHelper/LCPlane3D.h b/Utilities/DataHelper/DataHelper/LCPlane3D.h
new file mode 100644
index 0000000000000000000000000000000000000000..8f9bdc18e85fbbc3b80460bfca95bdfa5722a507
--- /dev/null
+++ b/Utilities/DataHelper/DataHelper/LCPlane3D.h
@@ -0,0 +1,120 @@
+#ifndef LCPlane3D_H
+#define LCPlane3D_H 1
+
+// #include "CLHEP/Vector/ThreeVector.h"
+#include <DataHelper/LCGeometryTypes.h>
+
+/** Definition of a LCPlane3D describing a geometrical plane in 3D space.
+ * @author T.Kraemer, DESY
+ * @version $Id: LCPlane3D.h,v 1.3 2006-10-19 15:59:26 tkraemer Exp $
+ */
+
+class LCPlane3D {
+
+public:
+
+ /**
+ * Constructor from four numbers - creates plane a*x+b*y+c*z+d=0.
+ * @param a
+ * @param b
+ * @param c
+ * @param d
+ */
+ LCPlane3D(double a = 0, double b = 0, double c = 1, double d = 0) ;
+
+ /**
+ * Constructor from normal and point.
+ * @param normal vector pointing in the direction of the normal.
+ * This vector does not have to be normalised.
+ * @param point Point on the plane.
+ */
+ LCPlane3D(LCVector3D normal, LCVector3D point) ;
+
+ /**
+ * Constructor from three different points.
+ * @param point1 Point on the plane.
+ * @param point2 Point on the plane.
+ * @param point3 Point on the plane.
+ */
+ LCPlane3D(LCVector3D point1, LCVector3D point2, LCVector3D point3) ;
+
+ /** Constructor for a plane using a normal and the Distance between Origen
+ * and the plane.
+ * @param normal vector pointing in the direction of the normal.
+ * This vector does not have to be normalised.
+ * @param distance distance is the distance from the origen to the plane.
+ */
+ LCPlane3D(LCVector3D normal, double distance) ;
+
+ /** Copy constructor.
+ * @param plane plane is an other LCPlane3D.
+ */
+ LCPlane3D(const LCPlane3D & plane) ;
+
+ /**
+ * Destructor. */
+ ~LCPlane3D() {}
+
+ /**
+ * Assignment. */
+
+ LCPlane3D & operator=(const LCPlane3D & rhs) ;
+
+ /**
+ * Returns the a-coefficient in the plane equation: a*x+b*y+c*z+d=0. */
+ double a() const ;
+
+ /**
+ * Returns the b-coefficient in the plane equation: a*x+b*y+c*z+d=0. */
+ double b() const ;
+
+ /**
+ * Returns the c-coefficient in the plane equation: a*x+b*y+c*z+d=0. */
+ double c() const ;
+
+ /**
+ * Returns the free member of the plane equation: a*x+b*y+c*z+d=0. */
+ double d() const ;
+
+ /**
+ * Returns normal. */
+ LCVector3D normal() const ;
+
+ /**
+ * Normalization. */
+ LCPlane3D & normalize() ;
+
+ /**
+ * Distance of a point to the plane.
+ * The value of the distance is
+ * - negative if the point and the origen are on the same side of the plane
+ * - positive if the Point and the origen are on opposite sides of the
+ * plane.
+ * @param point point is a point in space
+ */
+ double distance(const LCVector3D & point) const ;
+
+ /**
+ * Projection of a point on to the plane.
+ * @param point point is a point in space.
+ */
+ LCVector3D projectPoint(const LCVector3D & point) const ;
+
+ /**
+ * Projection of the origin onto the plane. */
+ LCVector3D projectPoint() const ;
+
+ /**
+ * Test for equality. */
+ bool operator==(const LCPlane3D & plane) const ;
+
+ /**
+ * Test for inequality. */
+ bool operator!=(const LCPlane3D & plane) const ;
+
+protected:
+ double _a, _b, _c, _d;
+};
+
+std::ostream & operator << (std::ostream &os, const LCPlane3D &p) ;
+#endif /* ifndef LCPlane3D_H */
diff --git a/Utilities/DataHelper/DataHelper/SimpleHelix.h b/Utilities/DataHelper/DataHelper/SimpleHelix.h
new file mode 100644
index 0000000000000000000000000000000000000000..3c289ebabcd5137f61526791b53470349f948f08
--- /dev/null
+++ b/Utilities/DataHelper/DataHelper/SimpleHelix.h
@@ -0,0 +1,127 @@
+#ifndef SimpleHelix_H
+#define SimpleHelix_H 1
+
+#include "DataHelper/Trajectory.h"
+
+/** Simple helix trajectory.
+ * @author T.Kraemer, DESY
+ * @version $Id: SimpleHelix.h,v 1.7 2007-06-20 18:47:25 samson Exp $
+ */
+
+class SimpleHelix : public Trajectory {
+
+public:
+
+ virtual ~SimpleHelix() {}
+
+ /** Construct Helix from canonical parameters.
+ */
+ SimpleHelix( double d0, double phi0, double omega,
+ double z0, double tanLambda,
+ LCVector3D referencePoint, LCErrorMatrix* errors=0) ;
+
+ /** Position at path length s - s==0 corresponds to P.C.A to the origin.
+ * @param s path length
+ * @param errors return argument - not computed if NULL
+ */
+ virtual LCVector3D getPosition(double s, LCErrorMatrix* errors=0) const ;
+
+ /** Direction at path length s, i.e. (dx/ds,dy/ds,dz/ds)
+ * @param s path length
+ * @param errors return argument - not computed if NULL
+ */
+ virtual LCVector3D getDirection(double s, LCErrorMatrix* errors=0) const ;
+
+ /** Full covariance Matrix of x,y,z,px,py,pz
+ * @param s path length
+ */
+ virtual LCErrorMatrix getCovarianceMatrix( double s) const ;
+
+ /** Pathlength at point on trajectory closest to given position.
+ * In order to get the distance use for example: <br>
+ * LCVector3D pt = t.getPosition( t.getPathAtClosestPoint( p ) ) ; <br>
+ * double d = LCVector3D( pt - p ).mag() ; <br>
+ */
+ virtual double getPathAt(const LCVector3D position ) const ;
+
+ /** get the radius of the helix.
+ */
+ virtual double getRadius() const ;
+
+ /*----------------------------------------------------------------------*/
+
+ /** Pathlength at closest intersection point with plane - undefined
+ * if pointExists==false.
+ */
+ virtual double getIntersectionWithPlane( LCPlane3D p, bool& pointExists) const ;
+
+ /** Pathlength at closest intersection point with cylinder - undefined
+ * if pointExists==false.
+ * @param cylinder cylinder object to intersect with
+ */
+ virtual double getIntersectionWithCylinder(const LCCylinder & cylinder,
+ bool & pointExists) const ;
+
+ /** Pathlength at the start and end point of the trajectory.
+ */
+ virtual double getStart() const ;
+ virtual double getEnd() const ;
+
+ /** Set pathlength at the start and end point of the trajectory.
+ * @param start pathlength of start point
+ * @param end pathlength of end point
+ */
+ virtual bool setStart(double s);
+ virtual bool setEnd(double s);
+ virtual bool setStartEnd(double start, double end);
+
+ virtual void printProperties();
+
+protected:
+
+ SimpleHelix() {}
+
+ virtual void init() ;
+
+ virtual double getCentreX() const ;
+ virtual double getCentreY() const ;
+ virtual double getWindingLength() const ;
+ virtual double getPitch();
+ double _d0;
+ double _phi0;
+ double _omega;
+ double _z0;
+ double _tanLambda;
+
+ // double _Bz;
+
+ double _helixStart;
+ double _helixEnd;
+
+ static const char* _names[];
+
+ static const double _a; // = 2.99792458E-4;
+ static const double _pi; // = 3.14159265358979323846;
+
+ LCVector3D _reference;
+ LCErrorMatrix* _errors;
+}; // class
+
+
+
+// /** Physical trajectory describing a (charged) particle's path in a B
+// * field and material.
+// * @author F.Gaede, DESY
+// * @version $Id: SimpleHelix.h,v 1.7 2007-06-20 18:47:25 samson Exp $
+// */
+
+// class PhysicalSimpleLine : public SimpleLine{
+
+// /** Particle's momentum at path length s. Implementations will have to have knowledge
+// * about the particle type, B-field and material.
+// */
+// virtual LCLorentzVector get4Momentum( double s ) const ;
+// }
+
+
+#endif /* ifndef SimpleLine_H */
diff --git a/Utilities/DataHelper/DataHelper/Trajectory.h b/Utilities/DataHelper/DataHelper/Trajectory.h
new file mode 100644
index 0000000000000000000000000000000000000000..3bf388cc465cfade770103dc54b218a06cdf68a1
--- /dev/null
+++ b/Utilities/DataHelper/DataHelper/Trajectory.h
@@ -0,0 +1,81 @@
+#ifndef Trajectory_H
+#define Trajectory_H 1
+
+#include <DataHelper/LCGeometryTypes.h>
+#include <DataHelper/LCPlane3D.h>
+#include <DataHelper/LCCylinder.h>
+
+/** Abstract trajectory interface describing a geometrical path in 3D space.
+ * @author F.Gaede, DESY
+ * @version $Id: Trajectory.h,v 1.4 2006-10-24 08:54:22 tkraemer Exp $
+ */
+
+class Trajectory {
+
+public:
+
+ /** Position at path length s - s==0 corresponds to P.C.A to the origin.
+ * @param s path length
+ * @param errors 3x3 matrix, return argument - not computed if NULL
+ */
+ virtual LCVector3D getPosition(double s, LCErrorMatrix* errors=0) const = 0;
+
+ /** Direction at path length s, i.e. (dx/ds,dy/ds,dz/ds)
+ * @param s path length
+ * @param errors 3x3 matrix, return argument - not computed if NULL
+ */
+ virtual LCVector3D getDirection(double s, LCErrorMatrix* errors=0) const = 0;
+
+ /** Full covariance Matrix of x,y,z,px,py,pz
+ * @param s path length
+ */
+ virtual LCErrorMatrix getCovarianceMatrix( double s) const = 0;
+
+
+ /** Pathlength at point on trajectory closest to given position.
+ * In order to get the distance use for example: <br>
+ * LCVector3D pt = t.getPosition( t.getPathAtClosestPoint( p ) ) ; <br>
+ * double d = LCVector3D( pt - p ).mag() ; <br>
+ */
+
+ virtual double getPathAt(const LCVector3D position ) const = 0;
+
+
+ /*----------------------------------------------------------------------*/
+
+ /** Pathlength at closest intersection point with plane - undefined
+ * if pointExists==false.
+ */
+ virtual double getIntersectionWithPlane( LCPlane3D p, bool& pointExists) const = 0 ;
+
+
+ /** Pathlength at closest intersection point with cylinder - undefined
+ * if pointExists==false.
+ * @param cylinder cylinder object to intersect with
+ */
+ virtual double getIntersectionWithCylinder(const LCCylinder & cylinder,
+ bool & pointExists) const = 0;
+
+ virtual ~Trajectory(){ ; } ;
+}; // class
+
+
+
+/** Physical trajectory describing a (charged) particle's path in a B
+ * field and material.
+ * @author F.Gaede, DESY
+ * @version $Id: Trajectory.h,v 1.4 2006-10-24 08:54:22 tkraemer Exp $
+ */
+
+class PhysicalTrajectory : public Trajectory{
+
+ /** Particle's momentum at path length s. Implementations will have to have knowledge
+ * about the particle type, B-field and material.
+ */
+ virtual LCLorentzVector get4Momentum( double s ) const = 0;
+
+ virtual ~PhysicalTrajectory(){ ; } ;
+};
+
+
+#endif /* ifndef Trajectory_H */
diff --git a/Utilities/DataHelper/src/Circle.cc b/Utilities/DataHelper/src/Circle.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c88c1a418349f34a8b8019a51d53919bca086c62
--- /dev/null
+++ b/Utilities/DataHelper/src/Circle.cc
@@ -0,0 +1,132 @@
+// Circle.cpp: implementation of the Circle class.
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "DataHelper/Circle.h"
+
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+Circle::Circle()
+{
+ this->m_dRadius=-1; // error checking
+}
+
+Circle::~Circle()
+{
+
+}
+
+Circle::Circle(CLHEP::Hep2Vector *V1, CLHEP::Hep2Vector *V2, CLHEP::Hep2Vector *V3)
+{
+ this->m_dRadius=-1; // error checking
+
+ CLHEP::Hep2Vector *pt1=new CLHEP::Hep2Vector; CLHEP::Hep2Vector *pt2=new CLHEP::Hep2Vector; CLHEP::Hep2Vector *pt3=new CLHEP::Hep2Vector;
+
+ *pt1=*V1; *pt2=*V2; *pt3=*V3;
+
+ if (!this->IsPerpendicular(pt1, pt2, pt3) ) this->CalcCircle(pt1, pt2, pt3);
+ else if (!this->IsPerpendicular(pt1, pt3, pt2) ) this->CalcCircle(pt1, pt3, pt2);
+ else if (!this->IsPerpendicular(pt2, pt1, pt3) ) this->CalcCircle(pt2, pt1, pt3);
+ else if (!this->IsPerpendicular(pt2, pt3, pt1) ) this->CalcCircle(pt2, pt3, pt1);
+ else if (!this->IsPerpendicular(pt3, pt2, pt1) ) this->CalcCircle(pt3, pt2, pt1);
+ else if (!this->IsPerpendicular(pt3, pt1, pt2) ) this->CalcCircle(pt3, pt1, pt2);
+ else {
+ std::cout << "The three pts are perpendicular to axis" << std::endl;
+ // pt1->trace(); pt2->trace(); pt3->trace();
+ delete pt1; delete pt2; delete pt3;
+ this->m_dRadius=-1;
+ return ;
+ }
+ delete pt1; delete pt2; delete pt3;
+
+}
+
+bool Circle::IsPerpendicular(CLHEP::Hep2Vector *pt1, CLHEP::Hep2Vector *pt2, CLHEP::Hep2Vector *pt3)
+// Check the given point are perpendicular to x or y axis
+{
+ double yDelta_a= pt2->y() - pt1->y();
+ double xDelta_a= pt2->x() - pt1->x();
+ double yDelta_b= pt3->y() - pt2->y();
+ double xDelta_b= pt3->x() - pt2->x();
+
+
+// TRACE(" yDelta_a: %f xDelta_a: %f \n",yDelta_a,xDelta_a);
+// TRACE(" yDelta_b: %f xDelta_b: %f \n",yDelta_b,xDelta_b);
+
+ // checking whether the line of the two pts are vertical
+ if (fabs(xDelta_a) <= 0.000000001 && fabs(yDelta_b) <= 0.000000001){
+ std::cout << "The points are pependicular and parallel to x-y axis" << std::endl;
+ return false;
+ }
+
+ if (fabs(yDelta_a) <= 0.0000001){
+// TRACE(" A line of two point are perpendicular to x-axis 1\n");
+ return true;
+ }
+ else if (fabs(yDelta_b) <= 0.0000001){
+// TRACE(" A line of two point are perpendicular to x-axis 2\n");
+ return true;
+ }
+ else if (fabs(xDelta_a)<= 0.000000001){
+// TRACE(" A line of two point are perpendicular to y-axis 1\n");
+ return true;
+ }
+ else if (fabs(xDelta_b)<= 0.000000001){
+// TRACE(" A line of two point are perpendicular to y-axis 2\n");
+ return true;
+ }
+ else return false ;
+}
+
+double Circle::CalcCircle(CLHEP::Hep2Vector *pt1, CLHEP::Hep2Vector *pt2, CLHEP::Hep2Vector *pt3)
+{
+ double yDelta_a= pt2->y() - pt1->y();
+ double xDelta_a= pt2->x() - pt1->x();
+ double yDelta_b= pt3->y() - pt2->y();
+ double xDelta_b= pt3->x() - pt2->x();
+
+ if (fabs(xDelta_a) <= 0.000000001 && fabs(yDelta_b) <= 0.000000001){
+ // TRACE("Calc cirlce \n");
+ this->m_Center.setX(0.5*(pt2->x() + pt3->x()));
+ this->m_Center.setY(0.5*(pt1->y() + pt2->y()));
+ // this->m_dRadius = m_Center.howNear(*pt1); // calc. radius
+ this->m_dRadius = sqrt((pt1->x()-m_Center.x())*(pt1->x()-m_Center.x()) + (pt1->y()-m_Center.y())*(pt1->y()-m_Center.y()));
+ // TRACE(" Center: %f %f %f\n", m_Center.x(), m_Center.y(), m_Center.z());
+ // TRACE(" radius: %f %f %f\n", length(&m_Center,pt1), length(&m_Center,pt2),length(&m_Center,pt3));
+
+ return this->m_dRadius;
+ }
+
+ // IsPerpendicular() assure that xDelta(s) are not zero
+ double aSlope=yDelta_a/xDelta_a; //
+ double bSlope=yDelta_b/xDelta_b;
+ if (fabs(aSlope-bSlope) <= 0.000000001){ // checking whether the given points are colinear.
+ std::cout << "The three pts are colinear" << std::endl;
+ return -1;
+ }
+
+ // calc center
+ this->m_Center.setX((aSlope*bSlope*(pt1->y() - pt3->y()) + bSlope*(pt1->x() + pt2->x())
+ - aSlope*(pt2->x()+pt3->x()) )/(2* (bSlope-aSlope) ));
+ this->m_Center.setY(-1*(m_Center.x() - (pt1->x()+pt2->x())/2)/aSlope + (pt1->y()+pt2->y())/2);
+
+ //this->m_dRadius= m_Center.howNear(*pt1); // calc. radius
+ this->m_dRadius = sqrt((pt1->x()-m_Center.x())*(pt1->x()-m_Center.x()) + (pt1->y()-m_Center.y())*(pt1->y()-m_Center.y()));
+ // TRACE(" Center: %f %f %f\n", m_Center.x(), m_Center.y(), m_Center.z());
+ // TRACE(" radius: %f %f %f\n", length(&m_Center,pt1), length(&m_Center,pt2),length(&m_Center,pt3));
+ return this->m_dRadius;
+}
+
+CLHEP::Hep2Vector* Circle::GetCenter()
+{
+ return &this->m_Center;
+
+}
+
+double Circle::GetRadius()
+{
+ return this->m_dRadius;
+}
diff --git a/Utilities/DataHelper/src/LCCylinder.cc b/Utilities/DataHelper/src/LCCylinder.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d2f0aeac93760c27b604ca6cc62028ff1bc7eb68
--- /dev/null
+++ b/Utilities/DataHelper/src/LCCylinder.cc
@@ -0,0 +1,221 @@
+#include <DataHelper/LCCylinder.h>
+#include <DataHelper/LCPlane3D.h>
+#include <DataHelper/LCLine3D.h>
+
+#include <iostream>
+
+#include <cmath>
+#include <float.h>
+#include <exception>
+
+LCCylinder::LCCylinder(LCVector3D point1,
+ LCVector3D point2,
+ double radius,
+ bool endPlane)
+{
+ _axisSstartPoint = point1;
+ _axisEndPoint = point2;
+
+ _radius = fabs(radius);
+ _endPlane = endPlane;
+}
+
+LCCylinder::LCCylinder(double radius,
+ LCVector3D point,
+ LCVector3D axis,
+ bool endPlane)
+{
+ _axisSstartPoint = point - axis ;
+ _axisEndPoint = point + axis ;
+
+ _radius = fabs(radius);
+ _endPlane = endPlane;
+}
+
+LCCylinder::LCCylinder(const LCCylinder & cylinder)
+{
+ _axisSstartPoint = cylinder._axisSstartPoint ;
+ _axisEndPoint = cylinder._axisEndPoint ;
+
+ _radius = cylinder._radius ;
+ _endPlane = cylinder._endPlane;
+}
+
+LCCylinder & LCCylinder::operator=(const LCCylinder & rhs)
+{
+ _axisSstartPoint = rhs._axisSstartPoint ;
+ _axisEndPoint = rhs._axisEndPoint ;
+
+ _radius = rhs._radius ;
+ _endPlane = rhs._endPlane;
+
+ return *this;
+}
+
+LCVector3D LCCylinder::startPoint() const
+{
+ return _axisSstartPoint;
+}
+
+LCVector3D LCCylinder::endPoint() const
+{
+ return _axisEndPoint ;
+}
+
+LCVector3D LCCylinder::axisDirection() const
+{
+ return (_axisEndPoint - _axisSstartPoint).unit();
+}
+
+double LCCylinder::length() const
+{
+ return (_axisEndPoint - _axisSstartPoint).mag();
+}
+
+double LCCylinder::radius() const
+{
+ return _radius;
+}
+
+double LCCylinder::distance(const LCVector3D & point) const
+{
+ int dummy ;
+ return (point - projectPoint( point, dummy ) ).mag() ;
+}
+
+LCVector3D LCCylinder::projectPoint(const LCVector3D & point, int & code) const
+{
+ // std::cout << "problem!" << std::endl;
+ // std::cout << "asp: " << _axisSstartPoint << " aep: " << _axisEndPoint << " point: " << point << std::endl;
+ LCLine3D a( _axisSstartPoint , axisDirection() );
+ // std::cout << a << std::endl;
+ double s = a.projectPoint( _axisSstartPoint ) ;
+ double e = a.projectPoint( _axisEndPoint ) ;
+ double p = a.projectPoint( point ) ;
+ double d = a.distance( point ) ;
+ // std::cout << "s: " << s << " e: " << e << " p: " << p << " d: " << d << std::endl;
+
+ double drp = fabs( d - radius() ) ;
+ double dsp = fabs( s - p ) ;
+ double dep = fabs( e - p ) ;
+
+
+ LCVector3D projection;
+
+ // classify in which region the point is located:
+ // point between the two planes at both ends
+ if ( p >= s && p <= e )
+ {
+ if (_endPlane && (d <= radius()) )
+ {
+ if ( (drp <= dsp) && (drp <= dep) )
+ {
+ projection = ( point - a.position(p) ).unit() ;
+ if (projection.mag() < 0.00001) projection = axisDirection().orthogonal().unit() ;
+ projection *= radius() ;
+ projection = a.position(p) + projection;
+ code = 3;
+ return projection;
+ }
+ else if (dsp <= dep)
+ {
+ LCPlane3D sPlane(-axisDirection(),_axisSstartPoint);
+ projection = sPlane.projectPoint( point );
+ code = 1;
+ return projection;
+ }
+ else // if ( dep < dsp )
+ {
+ LCPlane3D ePlane(axisDirection(),_axisEndPoint);
+ projection = ePlane.projectPoint( point );
+ code = 2;
+ return projection;
+ }
+ }
+ else
+ {
+ projection = ( point - a.position(p) ).unit() ;
+ if (projection.mag() < 0.00001) projection = axisDirection().orthogonal().unit() ;
+ projection *= radius() ;
+ projection = a.position(p) + projection;
+ code = 3;
+ return projection;
+ }
+ }
+ else // outside the two planes at the end
+ {
+ if (_endPlane && (d <= radius()) )
+ {
+ if ( p < s)
+ {
+ LCPlane3D sPlaneo(-axisDirection(),_axisSstartPoint);
+ projection = sPlaneo.projectPoint( point );
+ code = 1;
+ return projection;
+ }
+ else // if ( p > e )
+ {
+ LCPlane3D ePlaneo(axisDirection(),_axisEndPoint);
+ projection = ePlaneo.projectPoint( point );
+ code = 2;
+ return projection;
+ }
+ }
+ else
+ {
+ if ( p < s)
+ {
+ projection = ( point - a.position(p) ).unit() ;
+ if (projection.mag() < 0.00001) projection = axisDirection().orthogonal().unit() ;
+ projection *= radius() ;
+ projection = a.position(s) + projection;
+ code = 0;
+ return projection;
+ }
+ else // if ( p > e )
+ {
+ projection = ( point - a.position(p) ).unit() ;
+ if (projection.mag() < 0.00001) projection = axisDirection().orthogonal().unit() ;
+ projection *= radius() ;
+ projection = a.position(e) + projection;
+ code = 0;
+ return projection;
+ }
+ }
+ }
+
+ std::cout << "Classification faild!!!" << std::endl;
+
+ return projection;
+}
+
+bool LCCylinder::isInside(const LCVector3D & point) const
+{
+ LCLine3D a( _axisSstartPoint , axisDirection() );
+
+ if ( radius() < a.distance( point ) ) return false ;
+
+ double s = a.projectPoint( _axisSstartPoint ) ;
+ double e = a.projectPoint( _axisEndPoint ) ;
+ double p = a.projectPoint( point ) ;
+ if (p < s || p > e) return false;
+
+ return true ;
+}
+
+bool LCCylinder::operator==(const LCCylinder & rhs) const
+{
+ return (_radius == rhs._radius &&
+ _axisSstartPoint == rhs._axisSstartPoint &&
+ _axisEndPoint == rhs._axisEndPoint &&
+ _endPlane == rhs._endPlane ) ;
+}
+
+bool LCCylinder::operator!=(const LCCylinder & rhs) const
+{
+ return (_radius != rhs._radius ||
+ _axisSstartPoint != rhs._axisSstartPoint ||
+ _axisEndPoint != rhs._axisEndPoint ||
+ _endPlane != rhs._endPlane) ;
+}
+
diff --git a/Utilities/DataHelper/src/LCLine3D.cc b/Utilities/DataHelper/src/LCLine3D.cc
new file mode 100644
index 0000000000000000000000000000000000000000..58d1e587a13aa73326b13cfa5e689a9cd995ec68
--- /dev/null
+++ b/Utilities/DataHelper/src/LCLine3D.cc
@@ -0,0 +1,175 @@
+#include <DataHelper/LCLine3D.h>
+#include <DataHelper/LCPlane3D.h>
+#include <iostream>
+
+#include <cmath>
+#include <float.h>
+#include <exception>
+
+LCLine3D::LCLine3D()
+{
+ _reference.set(0.,0.,0.);
+ _point.set(0.,0.,0.);
+ _direction.set(1.,0.,0.);
+}
+
+LCLine3D::LCLine3D(const LCVector3D & point, const LCVector3D & direction)
+{
+ set( point, direction, LCVector3D(0.,0.,0.) );
+}
+
+LCLine3D::LCLine3D(const LCVector3D & point,
+ const LCVector3D & direction,
+ const LCVector3D & reference)
+{
+ set( point, direction, reference );
+}
+
+LCLine3D::LCLine3D(double d0, double phi0, double z0, double tanLambda)
+{
+ set( d0, phi0, z0, tanLambda, LCVector3D(0.,0.,0.) );
+}
+
+LCLine3D::LCLine3D(double d0, double phi0, double z0, double tanLambda,
+ const LCVector3D & reference)
+{
+ set( d0, phi0, z0, tanLambda, reference );
+}
+
+LCLine3D::LCLine3D(const LCLine3D & line)
+{
+ _point = line._point;
+ _direction = line._direction;
+ _reference = line._reference;
+}
+
+bool LCLine3D::set(const LCVector3D & point,
+ const LCVector3D & direction,
+ const LCVector3D & reference)
+{
+ // std::cout << "LCLine " << _reference << " " << point << " " << direction << std::endl;
+
+ _reference = reference;
+ _direction = direction.unit();
+ if (_direction.mag2() == 0)
+ {
+ return false;
+ }
+
+// calculate _point to be the PCA to the reference point according to the
+// definition given in LC-LC-DET-2006-004:
+// the x,y compnents have to beh teh PCA, the z compnente is calculated
+// after that.
+
+ LCVector3D p = point, d = _direction;
+ p.setZ(0.);
+ d.setZ(0.);
+
+ if (d.mag() != 0.)
+ {
+ double sFaktor = 1./d.mag();
+ d = d.unit();
+ double s = ( - p*d ) / d.mag2() ;
+ // x,y componentes:
+ // _point = p + s * d ;
+ // z component:
+ _point = ( (point + s*sFaktor*_direction) );
+ }
+ else
+ {
+ _point = point;
+ _point.setZ(0.);
+ }
+ // std::cout << "LCLine: " << _reference << " " << _point << " " << _direction << std::endl;
+ return true;
+}
+
+bool LCLine3D::set(double d0, double phi0, double z0, double tanLambda,
+ const LCVector3D & reference)
+{
+ _reference = reference;
+ _direction.set( cos(phi0), sin(phi0), tanLambda );
+ _direction = _direction.unit();
+ if (d0 == 0.)
+ {
+ _point.set(0.,0.,z0);
+ }
+ else
+ {
+ _point.set( ( d0*sin(phi0) ), ( d0*cos(phi0) ), z0 );
+ }
+ return true;
+}
+
+LCLine3D & LCLine3D::operator=(const LCLine3D & rhs)
+{
+ _point = rhs._point;
+ _direction = rhs._direction;
+ _reference = rhs._reference;
+
+ return *this;
+}
+
+LCVector3D LCLine3D::position(const double s) const
+{
+ return (_reference+_point + s*_direction) ;
+}
+
+LCVector3D LCLine3D::direction() const
+{
+ return _direction;
+}
+
+double LCLine3D::distance(const LCVector3D & point) const
+{
+ return ( point - position( projectPoint( point ) ) ).mag() ;
+}
+
+double LCLine3D::projectPoint(const LCVector3D & point) const
+{
+ // the last therm : (...) / _direction.mag2() is not there becaus
+ // the _direction vector is normalised.
+ // return ( 2*point*_direction - _point*_direction ) / _direction.mag2() ;
+ // return ( point*_direction - (_reference+_point)*_direction ) / _direction.mag2() ;
+ double x = ( point*_direction - (_reference+_point)*_direction ) / _direction.mag2() ;
+ // std::cout << "x: " << x << std::endl;
+ // std::cout << "point: " << point
+ // << " direction: " << _direction << " _point: " << _point
+ // << " d.mag: " << _direction.mag2() << std::endl;
+
+
+ return x;
+}
+
+bool LCLine3D::operator==(const LCLine3D & rhs) const
+{
+ return (_point == rhs._point &&
+ _direction == rhs._direction &&
+ _reference == rhs._reference);
+}
+
+bool LCLine3D::operator!=(const LCLine3D & rhs) const
+{
+ return (_point != rhs._point ||
+ _direction != rhs._direction ||
+ _reference != rhs._reference) ;
+}
+
+double LCLine3D::intersectionWithPlane(const LCPlane3D plane, bool& pointExists) const
+{
+ double c = direction() * plane.normal() ;
+
+ if (c == 0)
+ { // no interaction
+ pointExists = false;
+ return DBL_MAX;
+ }
+
+ pointExists = true;
+ return - ( position() * plane.normal() + plane.d() ) / c ;
+}
+
+std::ostream & operator << (std::ostream &os, const LCLine3D &l)
+{
+ return os << l.position() << "+s*" << l.direction() ;
+}
diff --git a/Utilities/DataHelper/src/LCPlane3D.cc b/Utilities/DataHelper/src/LCPlane3D.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d7e298eeda015f7588f4116cd740a6e3a8c17bd6
--- /dev/null
+++ b/Utilities/DataHelper/src/LCPlane3D.cc
@@ -0,0 +1,181 @@
+#include <DataHelper/LCPlane3D.h>
+
+#include <string>
+#include <sstream>
+#include <iostream>
+#include <cmath>
+#include <float.h>
+#include <exception>
+// #include <CLHEP/Vector/ThreeVector.h>
+
+LCPlane3D::LCPlane3D(double a, double b, double c, double d)
+{
+ _a = a;
+ _b = b;
+ _c = c;
+ _d = d;
+ normalize();
+
+ // _d has to be negative to give the distance with the normal vector pointing
+ // away from the Origen !!!
+ if (_d > 0)
+ {
+ _a = -_a;
+ _b = -_b;
+ _c = -_c;
+ _d = -_d;
+ }
+}
+
+LCPlane3D::LCPlane3D(LCVector3D normal, LCVector3D point)
+{
+ LCVector3D n = normal.unit();
+
+ _a = n.x();
+ _b = n.y();
+ _c = n.z();
+ _d = - n*point;
+}
+
+LCPlane3D::LCPlane3D(LCVector3D point1, LCVector3D point2, LCVector3D point3)
+{
+ LCVector3D n = ( (point2-point1).cross(point3-point1) ).unit();
+ _a = n.x() ;
+ _b = n.y() ;
+ _c = n.z() ;
+
+ _d = - n*point1;
+}
+
+LCPlane3D::LCPlane3D(LCVector3D normal, double distance)
+{
+ LCVector3D n = normal.unit() ;
+ _a = n.x();
+ _b = n.y();
+ _c = n.z();
+ _d = -distance;
+}
+
+LCPlane3D::LCPlane3D(const LCPlane3D & plane)
+ : _a(plane._a),
+ _b(plane._b),
+ _c(plane._c),
+ _d(plane._d)
+{}
+
+LCPlane3D & LCPlane3D::operator=(const LCPlane3D & rhs)
+{
+ _a = rhs._a;
+ _b = rhs._b;
+ _c = rhs._c;
+ _d = rhs._d;
+
+ return *this;
+}
+
+double LCPlane3D::a() const
+{
+ return _a ;
+}
+
+double LCPlane3D::b() const
+{
+ return _b ;
+}
+
+double LCPlane3D::c() const
+{
+ return _c ;
+}
+
+double LCPlane3D::d() const
+{
+ return _d ;
+}
+
+LCVector3D LCPlane3D::normal() const
+{
+ LCVector3D n(_a,_b,_c);
+ return n.unit();
+}
+
+LCPlane3D & LCPlane3D::normalize()
+{
+ double norm = sqrt( _a*_a + _b*_b + _c*_c );
+
+ if (norm > 0.)
+ {
+ _a /= norm;
+ _b /= norm;
+ _c /= norm;
+ _d /= norm;
+ }
+
+ return *this;
+}
+
+double LCPlane3D::distance(const LCVector3D & point) const
+{
+ return _a*point.x() + _b*point.y() + _c*point.z() + _d ;
+}
+
+LCVector3D LCPlane3D::projectPoint(const LCVector3D & point) const
+{
+ double k = distance(point) / ( _a*_a + _b*_b + _c*_c );
+ return LCVector3D( point.x()-_a*k, point.y()-_b*k, point.z()-_c*k);
+}
+
+LCVector3D LCPlane3D::projectPoint() const
+{
+ double k = -_d / ( _a*_a + _b*_b + _c*_c );
+ return LCVector3D( _a*k, _b*k, _c*k);
+}
+
+bool LCPlane3D::operator==(const LCPlane3D & plane) const
+{
+ return ( _a == plane._a &&
+ _b == plane._b &&
+ _c == plane._c &&
+ _d == plane._d );
+}
+
+bool LCPlane3D::operator!=(const LCPlane3D & plane) const
+{
+ return ( _a != plane._a ||
+ _b != plane._b ||
+ _c != plane._c ||
+ _d != plane._d );
+}
+
+std::ostream & operator << (std::ostream &os, const LCPlane3D &p)
+{
+ std::stringstream returnString ;
+ bool isFirst = true;
+
+ returnString << "(" ;
+ if (p.a() != 0)
+ {
+ returnString << p.a() << "*x" ;
+ isFirst = false;
+ }
+ if (p.b() != 0)
+ {
+ if (!isFirst && (p.b() > 0.) ) returnString << "+";
+ returnString << p.b() << "*y" ;
+ isFirst = false;
+ }
+ if (p.c() != 0)
+ {
+ if (!isFirst && (p.c() > 0) ) returnString << "+";
+ returnString << p.c() << "*z" ;
+ isFirst = false;
+ }
+ if (p.d() != 0)
+ {
+ if (!isFirst && (p.d() > 0) ) returnString << "+";
+ returnString << p.d() ;
+ }
+ returnString << "=0)" ;
+
+ return os << returnString.str() ;
+}
diff --git a/Utilities/DataHelper/src/SimpleHelix.cc b/Utilities/DataHelper/src/SimpleHelix.cc
new file mode 100644
index 0000000000000000000000000000000000000000..17011f5723be174b174c0ae1a52734730a443851
--- /dev/null
+++ b/Utilities/DataHelper/src/SimpleHelix.cc
@@ -0,0 +1,432 @@
+#include "DataHelper/SimpleHelix.h"
+#include "CLHEP/Matrix/Matrix.h"
+#include "CLHEP/Matrix/Vector.h"
+
+#include <iostream>
+#include <iomanip>
+#include <DataHelper/LCLine3D.h>
+#include <cmath>
+#include <float.h>
+#include <exception>
+
+const double SimpleHelix::_a = 2.99792458E-4;
+const double SimpleHelix::_pi = M_PI;
+
+SimpleHelix::SimpleHelix( double d0, double phi0, double omega,
+ double z0, double tanLambda,
+ LCVector3D referencePoint, LCErrorMatrix* errors)
+{
+ init();
+
+ _d0 = d0;
+ _phi0 = phi0;
+ _omega = omega;
+ _z0 = z0;
+ _tanLambda = tanLambda;
+
+ _reference = referencePoint;
+
+ // _Bz = Bz;
+
+ if ( errors == 0 )
+ {
+ _errors = new LCErrorMatrix(5,0);
+ }
+ else
+ {
+ if ( errors->num_row() != 5 )
+ throw("The error matrix has to be a 5x5 symmetric matrix in SimpleHelix Constructor");
+ *_errors = *errors;
+ }
+}
+
+void SimpleHelix::init()
+{
+ _d0 = 0;
+ _phi0 = 0;
+ _omega = 1;
+ _z0 = 0;
+ _tanLambda = 0;
+ // _Bz = 4;
+
+ _reference.set(0.,0.,0.);
+
+ _helixStart = -DBL_MAX;
+ _helixEnd = DBL_MAX;
+
+}
+
+double SimpleHelix::getCentreX() const
+{
+ return ( _reference.x() + ((1/_omega) - _d0) * sin(_phi0) );
+}
+
+double SimpleHelix::getCentreY() const
+{
+ return ( _reference.y() - ((1/_omega) - _d0) * cos(_phi0) );
+}
+
+double SimpleHelix::getWindingLength() const
+{
+ return 2*_pi*sqrt(1+_tanLambda*_tanLambda)/fabs(_omega);
+}
+
+LCVector3D SimpleHelix::getPosition(double s, LCErrorMatrix* errors) const
+{
+ LCVector3D x;
+
+ double xc = getCentreX();
+ double yc = getCentreY();
+
+ double varphi0 = _phi0 + ((_omega * _pi) / (2*fabs(_omega)));
+ double w = _omega / sqrt(1 + _tanLambda*_tanLambda);
+
+ x.setX(xc + fabs(1/_omega) * cos( w*s - varphi0 ) );
+ x.setY(yc + fabs(1/_omega) * sin( (-w*s) + varphi0) );
+ x.setZ(_reference.z() + _z0 + s*_tanLambda/sqrt(1 + _tanLambda*_tanLambda) );
+
+ return x;
+}
+
+LCVector3D SimpleHelix::getDirection(double s, LCErrorMatrix* errors) const
+{
+ LCVector3D t;
+
+ double varphi0 = _phi0 + ((_omega * _pi) / (2*fabs(_omega)));
+ double w = _omega / sqrt(1 + _tanLambda*_tanLambda);
+
+ t.setX(-w*fabs(1/_omega) * sin( (w*s) - varphi0) );
+ t.setY(-w*fabs(1/_omega) * cos( (-w*s) + varphi0) );
+ t.setZ(_tanLambda/sqrt(1 + _tanLambda*_tanLambda) );
+
+ return t.unit();
+}
+
+LCErrorMatrix SimpleHelix::getCovarianceMatrix( double s) const
+{
+ return LCErrorMatrix( 6 , 0 ) ;
+}
+
+double SimpleHelix::getPathAt(const LCVector3D position ) const
+{
+ double sStart = (position.z() - _reference.z() - _z0)
+ * sqrt(1 + _tanLambda*_tanLambda) / _tanLambda;
+
+ if (sStart <= _helixStart) sStart = _helixStart;
+ if (sStart >= _helixEnd) sStart = _helixEnd;
+
+ double startRange = sStart - getWindingLength();
+ double endRange = sStart + getWindingLength();
+
+ if (startRange <= _helixStart) startRange = _helixStart;
+ if (endRange >= _helixEnd) endRange = _helixEnd;
+
+ ///### int nSteps = 100;
+ int nSteps = 10;
+ double stepWidth = (endRange - startRange)/((double)nSteps);
+
+ double sOfMin = 0;
+ double distMinSQ = DBL_MAX;
+ double s = startRange;
+ LCVector3D x;
+ for (int i = 0; i <= nSteps; i++)
+ {
+ x = getPosition(s);
+ double distsq = (x-position).mag2() ;
+// ###### vector
+ if (distsq < distMinSQ)
+ {
+ distMinSQ = distsq;
+ sOfMin = s;
+ }
+ s += stepWidth;
+ }
+
+ double minStepWidth = 0.000001;
+ while (stepWidth > minStepWidth)
+ {
+ stepWidth /= 10;
+
+ x = getPosition(sOfMin + stepWidth);
+ double lp = (x-position).mag2();
+ x = getPosition(sOfMin - stepWidth);
+ double lm = (x-position).mag2();
+
+ double upOrDown = 0;
+ double lastDistance = 0;
+ nSteps = 10 ;
+ if (distMinSQ <= lp && distMinSQ <= lm)
+ {
+ nSteps = 0 ;
+ lastDistance = distMinSQ;
+ }
+ else if (lm<lp)
+ {
+ upOrDown = -1;
+ lastDistance = lm;
+ }
+ else if (lp<lm)
+ {
+ upOrDown = 1;
+ lastDistance = lp;
+ }
+ else
+ {
+ // cout << "Da laeuft was falsch!!! "
+ // << " distMinSQ " << distMinSQ << " lm " << lm << " lp " << lp << endl;
+ }
+
+ double step = 0, l = 0, lastStep = sOfMin;
+ for (int i = 1; i<=nSteps;i++)
+ {
+ step = sOfMin + (upOrDown*stepWidth*(double)i);
+ x = getPosition(step);
+ l = (x-position).mag2();
+ if (l <= lastDistance)
+ {
+ lastStep = step;
+ lastDistance = l;
+ }
+ }
+ sOfMin = lastStep ;
+ distMinSQ = lastDistance;
+ }
+ return sOfMin;
+}
+
+double SimpleHelix::getIntersectionWithPlane( LCPlane3D p,
+ bool& pointExists) const
+{
+
+ // Calculation of the region of s (parameter of Helix), where the helix
+ // interacts with the plane. This is done by two lines parallel to the
+ // axis of the helix (in this parametrisation the z-axis) in the distance
+ // of the helix radius. Both lines on opposite sides of the helix.
+ // One on the side nearest to the plane, the other exactly on the
+ // other side, away from the plane with respect to the Origen.
+
+ LCVector3D helixCentre( getCentreX(), getCentreY(), 0.);
+ LCVector3D normalXY = p.normal();
+ normalXY.setZ(0.);
+ // If plane is perpendicular to helix axis, set normalXY to a value
+ // perpendicular to the helix axis to enable intersection calculation.
+ if (normalXY.x() == 0. && normalXY.y() == 0. ) normalXY.set(1.,0.,0.);
+ normalXY = ( normalXY.unit() )/_omega;
+ LCVector3D lineDirection(0.,0.,1.);
+
+ LCLine3D frontLine( (helixCentre + normalXY) , lineDirection);
+ LCLine3D backLine( (helixCentre - normalXY) , lineDirection);
+ bool parallelBack, parallelFront;
+ double zStart = ( frontLine.position( frontLine.intersectionWithPlane(p,parallelFront) ) ).z();
+ double zEnd = ( backLine.position( backLine.intersectionWithPlane(p,parallelBack) ) ).z();
+ if ( !(parallelBack && parallelFront) )
+ { // plane is parallel to helix axis
+ if ( fabs( p.distance(helixCentre) ) > (1/_omega) )
+ { // Helix never hits plane!
+ pointExists = false ;
+ return 0;
+ }
+ else
+ {
+ pointExists = true ;
+ zStart = -DBL_MAX ;
+ zEnd = DBL_MAX ;
+ }
+ }
+
+ double sStart = ( zStart - _reference.z() - _z0)
+ * sqrt(1 + _tanLambda*_tanLambda) / _tanLambda;
+ double sEnd = ( zEnd - _reference.z() - _z0)
+ * sqrt(1 + _tanLambda*_tanLambda) / _tanLambda;
+
+ if (sStart > sEnd)
+ { // rong order for Start and Endpoint
+ double temp = sEnd;
+ sEnd = sStart;
+ sStart = temp;
+ }
+
+ // produce an artificial gap between sStart and sEnd
+ if (sStart > -DBL_MAX && sStart < DBL_MAX ) sStart -= 0.00001;
+ if (sEnd > -DBL_MAX && sEnd < DBL_MAX ) sEnd += 0.00001;
+
+ if ( (sStart < 0.) && (sEnd < 0.) )
+ { // Intersection is in backwards direction
+ pointExists = false ;
+ return 0;
+ }
+ else if ( (sStart < 0.) && (sEnd > 0.) )
+ { // intersection region starts in backwards direction
+ sStart = 0;
+ }
+
+ double s = sStart, sOld = -DBL_MAX;
+ double epsilon = 0.0000001;
+
+ while ( (s-sOld) > epsilon )
+ {
+ double d = fabs( p.distance( getPosition(s) ) ) ;
+ sOld = s;
+ s += d;
+ if (s > sEnd)
+ { // Problem! no intersection !
+ std::cout << "ERROR: No intersection found!!!" << std::endl;
+ pointExists = false ;
+ return 0;
+ }
+ }
+
+ pointExists = true ;
+ return s ;
+}
+
+double SimpleHelix::getIntersectionWithCylinder(const LCCylinder & cylinder,
+ bool & pointExists) const
+{
+ LCVector3D helixCentre( getCentreX(), getCentreY(), 0.);
+ LCVector3D axisDirection(0.,0.,1.);
+ LCLine3D axisLine(helixCentre , axisDirection);
+
+ LCVector3D middlePoint =
+ (cylinder.startPoint() + cylinder.endPoint())/2.;
+ double minDistance = sqrt( 0.25*cylinder.length()*cylinder.length()
+ + cylinder.radius()*cylinder.radius() );
+
+
+ if ( axisLine.distance(middlePoint) > (minDistance+getRadius()) )
+ {
+ pointExists = false ;
+ return 0 ;
+ }
+
+ double sProject = axisLine.projectPoint(middlePoint);
+ double sStart = sProject - minDistance ;
+ double sEnd = sProject + minDistance ;
+
+ double sHelixStart = ( (axisLine.position(sStart)).z() - _reference.z() - _z0)
+ * sqrt(1 + _tanLambda*_tanLambda) / _tanLambda;
+ double sHelixEnd = ( (axisLine.position(sEnd)).z() - _reference.z() - _z0)
+ * sqrt(1 + _tanLambda*_tanLambda) / _tanLambda;
+
+ if (sHelixStart > sHelixEnd)
+ {
+ double temp = sHelixStart;
+ sHelixStart = sHelixEnd;
+ sHelixEnd = temp;
+ }
+
+ if ( (sHelixStart < 0.) && (sHelixEnd < 0.) )
+ { // Intersection is in backwards direction
+ pointExists = false ;
+ return 0;
+ }
+ else if ( (sHelixStart < 0.) && (sHelixEnd > 0.) )
+ { // intersection region starts in backwards direction
+ sHelixStart = 0;
+ }
+
+ double s = sHelixStart, sOld = -DBL_MAX;
+ double epsilon = 0.0000001;
+
+ while ( (s-sOld) > epsilon )
+ {
+ double d = fabs( cylinder.distance( getPosition(s) ) ) ;
+ sOld = s;
+ if (s > sHelixEnd)
+ { // Problem! no intersection !
+ pointExists = false ;
+ return 0;
+ }
+ s += d;
+ }
+
+ pointExists = true ;
+ return s ;
+}
+
+double SimpleHelix::getStart() const
+{
+ return _helixStart;
+}
+
+double SimpleHelix::getEnd() const
+{
+ return _helixEnd;
+}
+
+bool SimpleHelix::setStart(double s)
+{
+ if (s <= _helixEnd)
+ {
+ _helixStart = s;
+ return true;
+ }
+
+ return false;
+}
+
+bool SimpleHelix::setEnd(double s)
+{
+ if (s >= _helixStart)
+ {
+ _helixEnd = s;
+ return true;
+ }
+
+ return false;
+}
+
+bool SimpleHelix::setStartEnd(double start, double end)
+{
+ if (start <= end)
+ {
+ _helixStart = start;
+ _helixEnd = end;
+ return true;
+ }
+
+ return false;
+}
+
+void SimpleHelix::printProperties()
+{
+ using namespace std;
+ int colWidth = 10;
+ cout << "*******************************************************************************" << endl;
+ cout << "*" << endl;
+ cout << "* Helix Parameters:" << endl;
+ cout << "* -----------------" << endl;
+ cout << "*" << endl;
+ cout << "* d_0 " << setw(colWidth) << _d0 << endl;
+ cout << "* phi_0 " << setw(colWidth) << _phi0 << endl;
+ cout << "* Omega " << setw(colWidth) << _omega
+ << " (Radius: " << (1/_omega) << ")" << endl;
+ cout << "* z_0 " << setw(colWidth) << _z0 << endl;
+ cout << "* tan lambda " << setw(colWidth) << _tanLambda << endl;
+ cout << "*" << endl;
+ cout << "* Reference point "
+ << setw(colWidth) << _reference.x() << " "
+ << setw(colWidth) << _reference.y() << " "
+ << setw(colWidth) << _reference.z() << endl;
+ cout << "* pitch " << setw(colWidth) << getPitch() << endl;
+ cout << "* Winding length " << setw(colWidth) << getWindingLength() << endl;
+ cout << "* Centre of arc "
+ << setw(colWidth) << getCentreX() << " "
+ << setw(colWidth) << getCentreY() << endl;
+ cout << "* Begin and end "
+ << setw(colWidth) << getStart()
+ << setw(colWidth) << getEnd() << endl;
+ cout << "*" << endl;
+ cout << "*******************************************************************************" << endl;
+}
+
+double SimpleHelix::getPitch()
+{
+ return 2*_pi*_tanLambda/fabs(_omega);
+}
+
+double SimpleHelix::getRadius() const
+{
+ return 1/_omega;
+}
+