diff --git a/include/JetBranches.h b/include/JetBranches.h index 5af2d35..eaf1e1a 100644 --- a/include/JetBranches.h +++ b/include/JetBranches.h @@ -158,6 +158,8 @@ class JetBranches : public CollectionBranches float _daughters_trackOmega[ LCT_JET_MAX ][ LCT_JET_PARTICLES_MAX ] {} ; /* PFOs track Omega */ float _daughters_trackZ0[ LCT_JET_MAX ][ LCT_JET_PARTICLES_MAX ] {} ; /* PFOs track Z0 */ float _daughters_trackTanLambda[ LCT_JET_MAX ][ LCT_JET_PARTICLES_MAX ] {} ; /* PFOs track TanLambda */ + float _daughters_trackSigmaD0[ LCT_JET_MAX ][ LCT_JET_PARTICLES_MAX ] {} ; // error on d0 from cov matrix + float _daughters_trackSigmaZ0[ LCT_JET_MAX ][ LCT_JET_PARTICLES_MAX ] {} ; // error on Z0 from cov matrix }; /* ----- end of class JetBranches ----- */ diff --git a/include/LCTuple.h b/include/LCTuple.h index 6f0bf80..a7e4efd 100644 --- a/include/LCTuple.h +++ b/include/LCTuple.h @@ -8,7 +8,9 @@ #include "MCParticleFromRelationBranches.h" #include "JetBranches.h" +#include "TrackBranches.h" #include "PIDBranches.h" +#include "RecoParticleBranches.h" using namespace lcio ; @@ -104,6 +106,8 @@ class LCTuple : public Processor { bool _jetColWriteParameters {}; bool _isolepColWriteParameters {}; bool _trkColWriteParameters {}; + bool _trkColStatesParameters {}; /* Enables writing extra trk states parameters */ + bool _trkColHitsParameters {}; /* Enables writing extra trk hits parameters */ bool _cluColWriteParameters {}; bool _sthColWriteParameters {}; bool _trhColWriteParameters {}; @@ -126,11 +130,11 @@ class LCTuple : public Processor { CWBranchesSet* _evtBranches {}; CollectionBranches* _mcpBranches {}; CollectionBranches* _mcpremoveoverlayBranches {}; - CollectionBranches* _recBranches {}; + RecoParticleBranches* _recBranches {}; // CollectionBranches* _jetBranches {}; JetBranches* _jetBranches {}; CollectionBranches* _isolepBranches {}; - CollectionBranches* _trkBranches {}; + TrackBranches* _trkBranches {}; CollectionBranches* _cluBranches {}; CollectionBranches* _sthBranches {}; CollectionBranches* _trhBranches {}; @@ -148,4 +152,4 @@ class LCTuple : public Processor { std::vector _pidBranchDefs {} ; } ; -#endif +#endif \ No newline at end of file diff --git a/include/LCTupleConf.h b/include/LCTupleConf.h index b944fef..9747dae 100644 --- a/include/LCTupleConf.h +++ b/include/LCTupleConf.h @@ -25,15 +25,15 @@ struct CollID : public lcrtrel::LCIntExtension {} ; #define LCT_COLLENTRIES_MAX 1000000 #define LCT_MCPARTICLE_MAX 1000000 #define LCT_RECOPARTICLE_MAX 500000 -#define LCT_TRACK_MAX 300000 -#define LCT_TRACKSTATE_MAX 1000000 +#define LCT_TRACK_MAX 1500000 +#define LCT_TRACKSTATE_MAX 100000 #define LCT_CLUSTER_MAX 500000 #define LCT_RELATION_MAX 1000000 #define LCT_SIMTRACKERHIT_MAX 2000000 #define LCT_TRACKERHIT_MAX 3000000 #define LCT_TRACKERRAWHIT_MAX 5000000 -#define LCT_SIMCALORIMETERHIT_MAX 1000000 -#define LCT_CALORIMETERHIT_MAX 1000000 +#define LCT_SIMCALORIMETERHIT_MAX 3000000 +#define LCT_CALORIMETERHIT_MAX 2000000 #define LCT_PARTICLEID_MAX 1000000 #define LCT_VERTEX_MAX 1000 #define LCT_JET_MAX 200 diff --git a/include/RecoParticleBranches.h b/include/RecoParticleBranches.h index 856c09f..56720d8 100644 --- a/include/RecoParticleBranches.h +++ b/include/RecoParticleBranches.h @@ -27,7 +27,7 @@ class RecoParticleBranches : public CollectionBranches { virtual void initBranches( TTree* tree, const std::string& prefix="" ) ; //const char* prefix=0) ; - virtual void fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt ) ; + virtual void fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt, const EVENT::LCCollection* colTracks, const EVENT::LCCollection* colCluster) ; virtual ~RecoParticleBranches() {} ; @@ -71,6 +71,8 @@ class RecoParticleBranches : public CollectionBranches { // EVENT::ClusterVec _clusters ; // EVENT::TrackVec _tracks ; + int _rcclid[ LCT_RECOPARTICLE_MAX ][5] {} ; + int _rctrid[ LCT_RECOPARTICLE_MAX ][5] {} ; } ; diff --git a/include/TrackBranches.h b/include/TrackBranches.h index 07633e4..fad6374 100644 --- a/include/TrackBranches.h +++ b/include/TrackBranches.h @@ -2,10 +2,8 @@ #define TrackBranches_h 1 #include "LCTupleConf.h" - #include "CollectionBranches.h" - class TTree ; namespace EVENT{ @@ -24,55 +22,80 @@ class TrackBranches : public CollectionBranches { public: TrackBranches() {} ; + + // Following function is used to access parameters set in a steering file + void writeTrkStatesParameters(bool param){ _writeTrkStatesParameters = param; }; + void writeTrkHitsParameters(bool param){ _writeTrkHitsParameters = param; }; - virtual void initBranches( TTree* tree, const std::string& prefix="" ) ; //const char* prefix=0) ; + virtual void initBranches( TTree* tree, const std::string& prefix="" ) ; virtual void fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt ) ; virtual ~TrackBranches() {} ; - private: - int _ntrk {} ; - - int _trori[ LCT_TRACK_MAX ] {} ; - - int _trtyp[ LCT_TRACK_MAX ] {} ; - float _trch2[ LCT_TRACK_MAX ] {} ; - int _trndf[ LCT_TRACK_MAX ] {} ; - float _tredx[ LCT_TRACK_MAX ] {} ; - float _trede[ LCT_TRACK_MAX ] {} ; - float _trrih[ LCT_TRACK_MAX ] {} ; - int _trthn[ LCT_TRACK_MAX ] {} ; // total number of hits - int _trthi[ LCT_TRACK_MAX ][50] {} ; // track hit indices - int _trthd[ LCT_TRACK_MAX ][50] {} ; // track hit subdetector - int _trshn[ LCT_TRACK_MAX ][12] {} ; - int _trnts[ LCT_TRACK_MAX ] {} ; - int _trfts[ LCT_TRACK_MAX ] {} ; - int _trsip[ LCT_TRACK_MAX ] {} ; // track stat atIP - int _trsfh[ LCT_TRACK_MAX ] {} ; // track stat atFirstHit - int _trslh[ LCT_TRACK_MAX ] {} ; // track stat atLastHit - int _trsca[ LCT_TRACK_MAX ] {} ; // track stat atCalorimeter - - int _ntrst {} ; - int _tsloc[ LCT_TRACKSTATE_MAX ] {} ; - float _tsdze[ LCT_TRACKSTATE_MAX ] {} ; - float _tsphi[ LCT_TRACKSTATE_MAX ] {} ; - float _tsome[ LCT_TRACKSTATE_MAX ] {} ; - float _tszze[ LCT_TRACKSTATE_MAX ] {} ; - float _tstnl[ LCT_TRACKSTATE_MAX ] {} ; - float _tscov[ LCT_TRACKSTATE_MAX ] [15] {} ; - float _tsrpx[ LCT_TRACKSTATE_MAX ] {} ; - float _tsrpy[ LCT_TRACKSTATE_MAX ] {} ; - float _tsrpz[ LCT_TRACKSTATE_MAX ] {} ; + bool _writeTrkStatesParameters {} ; /* Whether to write track states parameters */ + bool _writeTrkHitsParameters {} ; /* Whether to write track hits parameters */ + + int _ntrk {} ; // number of tracks + + int _trori[ LCT_TRACK_MAX ] {} ; // track index: ext + int _trtyp[ LCT_TRACK_MAX ] {} ; // type + float _trch2[ LCT_TRACK_MAX ] {} ; // chi2 + int _trndf[ LCT_TRACK_MAX ] {} ; // ndf + float _tredx[ LCT_TRACK_MAX ] {} ; // dEdx + float _trede[ LCT_TRACK_MAX ] {} ; // dEdxError( + float _trrih[ LCT_TRACK_MAX ] {} ; // RadiusOfInnermostHit + int _trnts[ LCT_TRACK_MAX ] {} ; // number of states + int _trthn[ LCT_TRACK_MAX ] {} ; // number of hits + int _trtvhn[ LCT_TRACK_MAX ] {} ; // number of vertex hits + int _trtihn[ LCT_TRACK_MAX ] {} ; // number of inner hits + int _trtohn[ LCT_TRACK_MAX ] {} ; // number of outer hits + int _trtnh[ LCT_TRACK_MAX ] {} ; // number of holes + int _trtout[ LCT_TRACK_MAX ] {} ; // number of outliers (#hit-ndf/2) + float _trome[ LCT_TRACK_MAX ] {} ; // omega at IP state + float _trtnl[ LCT_TRACK_MAX ] {} ; // tan lambda at IP state + float _trthe[ LCT_TRACK_MAX ] {} ; // theta at IP state + float _trdze[ LCT_TRACK_MAX ] {} ; // d0 at IP state + float _trzze[ LCT_TRACK_MAX ] {} ; // z0 at IP state + float _trphi[ LCT_TRACK_MAX ] {} ; // phi at IP state + //float _trcov[ LCT_TRACK_MAX ] [15] {} ; // covariance matrix at IP state (i.e. cov) + float _trk_sigmal0[ LCT_TRACK_MAX ] {} ; // sigma L0: cov[0][0] + float _trk_sigmal1[ LCT_TRACK_MAX ] {} ; // sigma L1: cov[1][1] + float _trk_sigmaphi[ LCT_TRACK_MAX ] {} ; // sigma phi: cov[2][2] + float _trk_sigmatheta[ LCT_TRACK_MAX ] {} ; // sigma theta: cov[3][3] + float _trk_sigmaqoverp[ LCT_TRACK_MAX ] {} ; // sigma qOverp: cov[4][4] + + // track states parameters + int _ntrst {} ; // total number of track states + int _trfts[ LCT_TRACK_MAX ] {} ; // track state index: ext() + int _trsip[ LCT_TRACK_MAX ] {} ; // track state index atIP + int _trsfh[ LCT_TRACK_MAX ] {} ; // track state index atFirstHit + int _trslh[ LCT_TRACK_MAX ] {} ; // track state index atLastHit + int _trsca[ LCT_TRACK_MAX ] {} ; // track state index atCalorimeter + int _tsloc[ LCT_TRACKSTATE_MAX ] {} ; // location + float _tsdze[ LCT_TRACKSTATE_MAX ] {} ; // D0 + float _tsphi[ LCT_TRACKSTATE_MAX ] {} ; // phi + float _tsome[ LCT_TRACKSTATE_MAX ] {} ; // omega + float _tszze[ LCT_TRACKSTATE_MAX ] {} ; // Z0 + float _tstnl[ LCT_TRACKSTATE_MAX ] {} ; // TanLambda + float _tsrpx[ LCT_TRACKSTATE_MAX ] {} ; // x + float _tsrpy[ LCT_TRACKSTATE_MAX ] {} ; // y + float _tsrpz[ LCT_TRACKSTATE_MAX ] {} ; // z + float _tscov[ LCT_TRACKSTATE_MAX ][15] {} ; // covariance matrix + + // track hits parameters + int _trthi[ LCT_TRACK_MAX ][50] {} ; // track hit indices + int _trthd[ LCT_TRACK_MAX ][50] {} ; // track hit subdetector + int _trshn[ LCT_TRACK_MAX ][12] {} ; // track hit per subdetector + float _trthx[ LCT_TRACK_MAX ][50] {} ; // track hit x coord + float _trthy[ LCT_TRACK_MAX ][50] {} ; // track hit y coord + float _trthz[ LCT_TRACK_MAX ][50] {} ; // track hit z coord // EVENT::TrackVec _tracks ; // EVENT::TrackerHitVec _hits ; } ; -#endif - - - +#endif \ No newline at end of file diff --git a/include/VertexBranches.h b/include/VertexBranches.h index 027357b..5ad1712 100644 --- a/include/VertexBranches.h +++ b/include/VertexBranches.h @@ -42,6 +42,7 @@ class VertexBranches : public CollectionBranches { float _vtcov[ LCT_VERTEX_MAX ][6] {} ; float _vtpar[ LCT_VERTEX_MAX ][6] {} ; //arbitrary value -- has to be checked if enough fields + float _vttrchi[LCT_VERTEX_MAX][30] {} ; } ; #endif diff --git a/src/JetBranches.cc b/src/JetBranches.cc index e1a2d5e..3d0c674 100644 --- a/src/JetBranches.cc +++ b/src/JetBranches.cc @@ -58,7 +58,6 @@ void JetBranches::initBranches( TTree* tree, const std::string& pre){ tree->Branch( (pre+"njet").c_str() , &_njet , (pre+"njet/I").c_str() ) ; - // ------------ Default Jet parameters ------------------// tree->Branch( (pre+"jmox").c_str() , _jmox , (pre+"jmox["+pre+"njet]/F").c_str() ) ; tree->Branch( (pre+"jmoy").c_str() , _jmoy , (pre+"jmoy["+pre+"njet]/F").c_str() ) ; @@ -87,9 +86,9 @@ void JetBranches::initBranches( TTree* tree, const std::string& pre){ tree->Branch( (pre+"jnpid" ).c_str() , &_jnpid , (pre+"jnpid/I").c_str() ) ; tree->Branch( (pre+"npfojet").c_str(), &_njetpfo , (pre+"npfojet["+pre+"njet]/I").c_str() ) ; - tree->Branch( (pre+"rcidx").c_str(), &_jetpfoori , (pre+"rcidx["+pre+"njet][LCT_JET_PARTICLES_MAX]/I").c_str() ) ; + tree->Branch( (pre+"rcidx").c_str(), &_jetpfoori , (pre+"rcidx["+pre+"njet][200]/I").c_str() ) ; } // end if - + // 200 //PFO branches if(_writeDaughtersParameters) { @@ -97,17 +96,19 @@ void JetBranches::initBranches( TTree* tree, const std::string& pre){ tree->Branch( (pre+"ndaughters").c_str(), _ndaughters , (pre+"ndaughters["+pre+"njet]/I").c_str() ) ; tree->Branch( (pre+"ntracks").c_str(), _ntracks , (pre+"ntracks["+pre+"njet]/I").c_str() ) ; tree->Branch( (pre+"nclusters").c_str(), _nclusters , (pre+"nclusters["+pre+"njet]/I").c_str() ) ; - tree->Branch( (pre+"daughters_PX").c_str(), _daughters_PX , (pre+"daughters_PX["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_PY").c_str(), _daughters_PY , (pre+"daughters_PY["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_PZ").c_str(), _daughters_PZ , (pre+"daughters_PZ["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_E").c_str(), _daughters_E , (pre+"daughters_E["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_M").c_str(), _daughters_M , (pre+"daughters_M["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_Q").c_str(), _daughters_Q , (pre+"daughters_Q["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_trackD0").c_str(), _daughters_trackD0 , (pre+"daughters_trackD0["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_trackPhi").c_str(), _daughters_trackPhi , (pre+"daughters_trackPhi["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_trackOmega").c_str(), _daughters_trackOmega , (pre+"daughters_trackOmega["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_trackZ0").c_str(), _daughters_trackZ0 , (pre+"daughters_trackZ0["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; - tree->Branch( (pre+"daughters_trackTanLambda").c_str(), _daughters_trackTanLambda , (pre+"daughters_trackTanLambda["+pre+"njet][LCT_JET_PARTICLES_MAX]/F").c_str() ) ; + tree->Branch( (pre+"daughters_PX").c_str(), _daughters_PX , (pre+"daughters_PX["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_PY").c_str(), _daughters_PY , (pre+"daughters_PY["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_PZ").c_str(), _daughters_PZ , (pre+"daughters_PZ["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_E").c_str(), _daughters_E , (pre+"daughters_E["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_M").c_str(), _daughters_M , (pre+"daughters_M["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_Q").c_str(), _daughters_Q , (pre+"daughters_Q["+pre+"njet][200]/I").c_str() ) ; + tree->Branch( (pre+"daughters_trackD0").c_str(), _daughters_trackD0 , (pre+"daughters_trackD0["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_trackPhi").c_str(), _daughters_trackPhi , (pre+"daughters_trackPhi["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_trackOmega").c_str(), _daughters_trackOmega , (pre+"daughters_trackOmega["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_trackZ0").c_str(), _daughters_trackZ0 , (pre+"daughters_trackZ0["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_trackTanLambda").c_str(), _daughters_trackTanLambda , (pre+"daughters_trackTanLambda["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_trackSigmaD0").c_str(), _daughters_trackSigmaD0 , (pre+"daughters_trackSigmaD0["+pre+"njet][200]/F").c_str() ) ; + tree->Branch( (pre+"daughters_trackSigmaZ0").c_str(), _daughters_trackSigmaZ0 , (pre+"daughters_trackSigmaZ0["+pre+"njet][200]/F").c_str() ) ; } @@ -200,17 +201,17 @@ void JetBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt ) if(_writeDaughtersParameters) { for ( size_t j = 0; j < LCT_JET_PARTICLES_MAX ; ++j ) { - _daughters_PX[ i ][ j ] = 0 ; - _daughters_PY[ i ][ j ] = 0 ; - _daughters_PZ[ i ][ j ] = 0 ; - _daughters_E[ i ][ j ] = 0 ; - _daughters_M[ i ][ j ] = 0 ; - _daughters_Q[ i ][ j ] = 0 ; - _daughters_trackD0[ i ][ j ] = 0; - _daughters_trackPhi[ i ][ j ] = 0; - _daughters_trackOmega[ i ][ j ] = 0; - _daughters_trackZ0[ i ][ j ] = 0; - _daughters_trackTanLambda[ i ][ j ] = 0; + _daughters_PX[ i ][ j ] = 0. ; + _daughters_PY[ i ][ j ] = 0. ; + _daughters_PZ[ i ][ j ] = 0. ; + _daughters_E[ i ][ j ] = 0. ; + _daughters_M[ i ][ j ] = 0. ; + _daughters_Q[ i ][ j ] = 0. ; + _daughters_trackD0[ i ][ j ] = 0.; + _daughters_trackPhi[ i ][ j ] = 0.; + _daughters_trackOmega[ i ][ j ] = 0.; + _daughters_trackZ0[ i ][ j ] = 0.; + _daughters_trackTanLambda[ i ][ j ] = 0.; } } } @@ -310,14 +311,15 @@ void JetBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt ) if (abs(_daughters_Q[ i ][ partid ])==0) nclusters++; auto tracks = particles[partid]->getTracks(); - //std::cout << "ntracks = " << tracks.size() << " charge = " << _daughters_Q[ i ][ partid ] << std::endl; - if (tracks.size()>0) { + if (tracks.size()>0 && tracks[0] != NULL) { _daughters_trackD0[ i ][ partid ] = tracks[0]->getD0(); _daughters_trackPhi[ i ][ partid ] = tracks[0]->getPhi(); _daughters_trackOmega[ i ][ partid ] = tracks[0]->getOmega(); _daughters_trackZ0[ i ][ partid ] = tracks[0]->getZ0(); _daughters_trackTanLambda[ i ][ partid ] = tracks[0]->getTanLambda(); + _daughters_trackSigmaD0[ i ][ partid ] = tracks[0]->getCovMatrix()[0]; + _daughters_trackSigmaZ0[ i ][ partid ] = tracks[0]->getCovMatrix()[2]; } } diff --git a/src/LCTuple.cc b/src/LCTuple.cc index 6eec968..5fbb8c5 100644 --- a/src/LCTuple.cc +++ b/src/LCTuple.cc @@ -184,6 +184,18 @@ registerProcessorParameter( "WriteIsoLepCollectionParameters" , false ); + registerProcessorParameter( "TrackCollectionStatesParameters" , + "Switch to write out states parameters for tracks", + _trkColStatesParameters , + false + ); + + registerProcessorParameter( "TrackCollectionHitsParameters" , + "Switch to write out hits parameters for tracks", + _trkColHitsParameters , + false + ); + registerInputCollection( LCIO::CLUSTER, "ClusterCollection" , "Name of the Cluster collection" , @@ -387,6 +399,8 @@ void LCTuple::init() { if( _trkColName.size() ) { _trkBranches = new TrackBranches ; _trkBranches->writeParameters(_trkColWriteParameters); + _trkBranches->writeTrkStatesParameters(_trkColStatesParameters); /* pass the value to TrackBranches */ + _trkBranches->writeTrkHitsParameters(_trkColHitsParameters); /* pass the value to TrackBranches */ _trkBranches->initBranches( _tree ) ; } @@ -551,7 +565,7 @@ void LCTuple::processEvent( LCEvent * evt ) { if( mcpRemoveOverlayCol ) _mcpremoveoverlayBranches->fill( mcpRemoveOverlayCol , evt ) ; if( recCol ) { - _recBranches->fill( recCol , evt ) ; + _recBranches->fill( recCol , evt, trkCol , cluCol ) ; for( auto pidb : _pidBranchesVec ) pidb->fill( recCol , evt ) ; } @@ -687,4 +701,4 @@ void LCTuple::decodePIDBranchDefinitions(){ } streamlog_out( DEBUG5 ) << std::endl ; } -} +} \ No newline at end of file diff --git a/src/MCParticleBranches.cc b/src/MCParticleBranches.cc index f7088c5..97a565a 100644 --- a/src/MCParticleBranches.cc +++ b/src/MCParticleBranches.cc @@ -66,11 +66,14 @@ void MCParticleBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* ev if (_writeparameters) CollectionBranches::fill(col, evt); _nmc = col->getNumberOfElements() ; - + int nmax = 500; + if(_nmc > nmax) _nmc = nmax; + for(int i=0 ; i < _nmc ; ++i){ lcio::MCParticle* mcp = static_cast( col->getElementAt(i) ) ; + _mcori[i] = mcp->ext(); _mcpdg[ i ] = mcp->getPDG() ; @@ -111,5 +114,7 @@ void MCParticleBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* ev _mcda3[ i ] = ( p1.size() > 3 ? p1[3]->ext() - 1 : -1 ) ; _mcda4[ i ] = ( p1.size() > 4 ? p1[4]->ext() - 1 : -1 ) ; + } + } diff --git a/src/RecoParticleBranches.cc b/src/RecoParticleBranches.cc index 0ce6b95..d681d19 100644 --- a/src/RecoParticleBranches.cc +++ b/src/RecoParticleBranches.cc @@ -7,6 +7,7 @@ #include "EVENT/Vertex.h" #include "TTree.h" +#include void RecoParticleBranches::initBranches( TTree* tree, const std::string& pre){ @@ -54,10 +55,13 @@ void RecoParticleBranches::initBranches( TTree* tree, const std::string& pre){ tree->Branch( (pre+"pillh").c_str() , _pillh , (pre+"pillh["+pre+"npid]/F").c_str() ) ; tree->Branch( (pre+"pialg").c_str() , _pialg , (pre+"pialg["+pre+"npid]/I").c_str() ) ; + //tree->Branch( (pre+"rcclid").c_str() , _rcclid , (pre+"rcclid["+pre+"nrec]["+pre+"rcncl]/I").c_str() ) ; + tree->Branch( (pre+"rcclid").c_str() , _rcclid , (pre+"rcclid["+pre+"nrec][5]/I").c_str() ) ; + tree->Branch( (pre+"rctrid").c_str() , _rctrid , (pre+"rctrid["+pre+"nrec][5]/I").c_str() ) ; } -void RecoParticleBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt ){ +void RecoParticleBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt, const EVENT::LCCollection* colTracks, const EVENT::LCCollection* colCluster){ if( !col ) return ; @@ -104,7 +108,8 @@ void RecoParticleBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* _pialg[ i ] = pid->getAlgorithmType() ; } - + std::vector usedClusters; + std::vector usedTracks; //------ fill the Reconstructed particle ---------------------------- for(int i=0 ; i < _nrec ; ++i){ @@ -146,7 +151,55 @@ void RecoParticleBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* _rcnrp[ i ] = rec->getParticles().size(); _rcftr[ i ] = ( rec->getTracks().size() > 0 ? rec->getTracks()[0]->ext() - 1 : -1 ) ; + + for(int s = 0; s < 5; s++){ + _rcclid[i][s] = -1; + _rctrid[i][s] = -1; + } + if(!colCluster) continue; + lcio::ClusterVec clusters = rec->getClusters(); + lcio::Cluster* temp_clus = NULL; + + for(unsigned int c = 0; c < clusters.size(); c++){ + + for(int ccoll = 0; ccoll < colCluster->getNumberOfElements(); ccoll++){ + + if(std::find(usedClusters.begin(), usedClusters.end(), ccoll) != usedClusters.end() ) + continue; + + temp_clus = static_cast( colCluster->getElementAt(ccoll) ); + + if(temp_clus->id() != clusters[c]->id()) continue; + + _rcclid[i][c] = ccoll; + usedClusters.push_back(ccoll); + break; + } + } + + if(!colTracks) continue; + lcio::TrackVec tracks = rec->getTracks(); + lcio::Track* temp_trk = NULL; + + for(unsigned int t = 0; t < tracks.size(); t++){ + + for(int tcoll = 0; tcoll < colTracks->getNumberOfElements(); tcoll++){ + + if(std::find(usedTracks.begin(), usedTracks.end(), tcoll) != usedTracks.end() ) + continue; + + temp_trk = static_cast( colTracks->getElementAt(tcoll) ); + + if(temp_trk->id() != tracks[t]->id()) continue; + + _rctrid[i][t] = tcoll; + usedTracks.push_back(tcoll); + break; + } + } + + usedClusters.clear(); } } diff --git a/src/TrackBranches.cc b/src/TrackBranches.cc index 5167c7d..3a2de96 100644 --- a/src/TrackBranches.cc +++ b/src/TrackBranches.cc @@ -7,6 +7,7 @@ #include "marlin/VerbosityLevels.h" #include "TTree.h" +#include "TMath.h" void TrackBranches::initBranches( TTree* tree, const std::string& pre){ @@ -19,7 +20,6 @@ void TrackBranches::initBranches( TTree* tree, const std::string& pre){ if (_writeparameters) CollectionBranches::initBranches(tree, (pre+"tr").c_str()); tree->Branch( (pre+"ntrk").c_str() , &_ntrk , (pre+"ntrk/I").c_str() ) ; - tree->Branch( (pre+"trori").c_str() , _trori , (pre+"trori["+pre+"ntrk]/I").c_str() ) ; tree->Branch( (pre+"trtyp").c_str() , _trtyp , (pre+"trtyp["+pre+"ntrk]/I").c_str() ) ; @@ -29,29 +29,55 @@ void TrackBranches::initBranches( TTree* tree, const std::string& pre){ tree->Branch( (pre+"trede").c_str() , _trede , (pre+"trede["+pre+"ntrk]/F").c_str() ) ; tree->Branch( (pre+"trrih").c_str() , _trrih , (pre+"trrih["+pre+"ntrk]/F").c_str() ) ; tree->Branch( (pre+"trthn").c_str() , _trthn , (pre+"trthn["+pre+"ntrk]/I").c_str() ) ; - tree->Branch( (pre+"trthi").c_str() , _trthi , (pre+"trthi["+pre+"ntrk][50]/I").c_str() ) ; - tree->Branch( (pre+"trshn").c_str() , _trshn , (pre+"trshn["+pre+"ntrk][12]/I").c_str() ) ; - tree->Branch( (pre+"trthd").c_str() , _trthd , (pre+"trthd["+pre+"ntrk][50]/I").c_str() ) ; tree->Branch( (pre+"trnts").c_str() , _trnts , (pre+"trnts["+pre+"ntrk]/I").c_str() ) ; - tree->Branch( (pre+"trfts").c_str() , _trfts , (pre+"trfts["+pre+"ntrk]/I").c_str() ) ; - tree->Branch( (pre+"trsip").c_str() , _trsip , (pre+"trsip["+pre+"ntrk]/I").c_str() ) ; - tree->Branch( (pre+"trsfh").c_str() , _trsfh , (pre+"trsfh["+pre+"ntrk]/I").c_str() ) ; - tree->Branch( (pre+"trslh").c_str() , _trslh , (pre+"trslh["+pre+"ntrk]/I").c_str() ) ; - tree->Branch( (pre+"trsca").c_str() , _trsca , (pre+"trsca["+pre+"ntrk]/I").c_str() ) ; - - tree->Branch( (pre+"ntrst").c_str() , &_ntrst , (pre+"ntrst/I").c_str() ) ; - - tree->Branch( (pre+"tsloc").c_str() , _tsloc , (pre+"tsloc["+pre+"ntrst]/I").c_str() ) ; - tree->Branch( (pre+"tsdze").c_str() , _tsdze , (pre+"tsdze["+pre+"ntrst]/F").c_str() ) ; - tree->Branch( (pre+"tsphi").c_str() , _tsphi , (pre+"tsphi["+pre+"ntrst]/F").c_str() ) ; - tree->Branch( (pre+"tsome").c_str() , _tsome , (pre+"tsome["+pre+"ntrst]/F").c_str() ) ; - tree->Branch( (pre+"tszze").c_str() , _tszze , (pre+"tszze["+pre+"ntrst]/F").c_str() ) ; - tree->Branch( (pre+"tstnl").c_str() , _tstnl , (pre+"tstnl["+pre+"ntrst]/F").c_str() ) ; - tree->Branch( (pre+"tscov").c_str() , _tscov , (pre+"tscov["+pre+"ntrst][15]/F").c_str() ) ; - tree->Branch( (pre+"tsrpx").c_str() , _tsrpx , (pre+"tsrpx["+pre+"ntrst]/F").c_str() ) ; - tree->Branch( (pre+"tsrpy").c_str() , _tsrpy , (pre+"tsrpy["+pre+"ntrst]/F").c_str() ) ; - tree->Branch( (pre+"tsrpz").c_str() , _tsrpz , (pre+"tsrpz["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"trtvhn").c_str() , _trtvhn , (pre+"trtvhn["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trtihn").c_str() , _trtihn , (pre+"trtihn["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trtohn").c_str() , _trtohn , (pre+"trtohn["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trtnh").c_str() , _trtnh , (pre+"trtnh["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trome").c_str() , _trome , (pre+"trome["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trtnl").c_str() , _trtnl , (pre+"trtnl["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trthe").c_str() , _trthe , (pre+"trthe["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trdze").c_str() , _trdze , (pre+"trdze["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trzze").c_str() , _trzze , (pre+"trzze["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trphi").c_str() , _trphi , (pre+"trphi["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trtout").c_str() , _trtout , (pre+"trtout["+pre+"ntrk]/I").c_str() ) ; + + // tree->Branch( (pre+"trcov").c_str() , _trcov , (pre+"trcov["+pre+"ntrk][15]/F").c_str() ) ; + tree->Branch( (pre+"trk_sigmal0").c_str() , _trk_sigmal0 , (pre+"trk_sigmal0["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trk_sigmal1").c_str() , _trk_sigmal1 , (pre+"trk_sigmal1["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trk_sigmaphi").c_str() , _trk_sigmaphi , (pre+"trk_sigmaphi["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trk_sigmatheta").c_str() , _trk_sigmatheta , (pre+"trk_sigmatheta["+pre+"ntrk]/F").c_str() ) ; + tree->Branch( (pre+"trk_sigmaqoverp").c_str() , _trk_sigmaqoverp , (pre+"trk_sigmaqoverp["+pre+"ntrk]/F").c_str() ) ; + + if (_writeTrkStatesParameters) { + tree->Branch( (pre+"trfts").c_str() , _trfts , (pre+"trfts["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trsip").c_str() , _trsip , (pre+"trsip["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trsfh").c_str() , _trsfh , (pre+"trsfh["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trslh").c_str() , _trslh , (pre+"trslh["+pre+"ntrk]/I").c_str() ) ; + tree->Branch( (pre+"trsca").c_str() , _trsca , (pre+"trsca["+pre+"ntrk]/I").c_str() ) ; + + tree->Branch( (pre+"ntrst").c_str() , &_ntrst , (pre+"ntrst/I").c_str() ) ; + tree->Branch( (pre+"tsloc").c_str() , _tsloc , (pre+"tsloc["+pre+"ntrst]/I").c_str() ) ; + tree->Branch( (pre+"tsdze").c_str() , _tsdze , (pre+"tsdze["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tsphi").c_str() , _tsphi , (pre+"tsphi["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tsome").c_str() , _tsome , (pre+"tsome["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tszze").c_str() , _tszze , (pre+"tszze["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tstnl").c_str() , _tstnl , (pre+"tstnl["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tsrpx").c_str() , _tsrpx , (pre+"tsrpx["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tsrpy").c_str() , _tsrpy , (pre+"tsrpy["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tsrpz").c_str() , _tsrpz , (pre+"tsrpz["+pre+"ntrst]/F").c_str() ) ; + tree->Branch( (pre+"tscov").c_str() , _tscov , (pre+"tscov["+pre+"ntrst][15]/F").c_str() ) ; + } + + if(_writeTrkHitsParameters) { + tree->Branch( (pre+"trthi").c_str() , _trthi , (pre+"trthi["+pre+"ntrk][50]/I").c_str() ) ; + tree->Branch( (pre+"trshn").c_str() , _trshn , (pre+"trshn["+pre+"ntrk][12]/I").c_str() ) ; + tree->Branch( (pre+"trthd").c_str() , _trthd , (pre+"trthd["+pre+"ntrk][50]/I").c_str() ) ; + tree->Branch( (pre+"trthx").c_str() , _trthx , (pre+"trthx["+pre+"ntrk][50]/F").c_str() ) ; + tree->Branch( (pre+"trthy").c_str() , _trthy , (pre+"trthy["+pre+"ntrk][50]/F").c_str() ) ; + tree->Branch( (pre+"trthz").c_str() , _trthz , (pre+"trthz["+pre+"ntrk][50]/F").c_str() ) ; + } } @@ -69,102 +95,137 @@ void TrackBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt ){ if (_writeparameters) CollectionBranches::fill(col, evt); _ntrk = col->getNumberOfElements() ; + + if (_writeTrkStatesParameters) { + //--------- create a helper vector with track states first ------------------------------- + std::vector tsV ; + tsV.reserve( col->getNumberOfElements() * 15 ) ; - //--------- create a helper vector with track states first ------------------------------- - std::vector tsV ; - tsV.reserve( col->getNumberOfElements() * 4 ) ; - - for(int i=0, nrc = col->getNumberOfElements() ; i < nrc ; ++i ){ - - lcio::Track* trk = static_cast( col->getElementAt( i) ) ; - - const EVENT::TrackStateVec & tss = trk->getTrackStates() ; - - for(int j=0, nts = tss.size() ; jext() = tsV.size() ; + for (int i=0, nrc = col->getNumberOfElements() ; i < nrc ; ++i ) { + lcio::Track* trk = static_cast( col->getElementAt( i) ) ; + const EVENT::TrackStateVec & tss = trk->getTrackStates() ; + for (int j=0, nts = tss.size() ; jgetLocation() ; - _tsdze[ i ] = ts->getD0() ; - _tsphi[ i ] = ts->getPhi() ; - _tsome[ i ] = ts->getOmega() ; - _tszze[ i ] = ts->getZ0() ; - _tstnl[ i ] = ts->getTanLambda() ; - _tsrpx[ i ] = ts->getReferencePoint()[0] ; - _tsrpy[ i ] = ts->getReferencePoint()[1] ; - _tsrpz[ i ] = ts->getReferencePoint()[2] ; - - for(int j=0;j<15;++j){ - _tscov[ i ][ j ] = ts->getCovMatrix()[j] ; + _ntrst = tsV.size() ; + streamlog_out( DEBUG ) << " total number of track states : " << _ntrst << std::endl ; + + //---------- fill the track states --------------------------------------- + for (int i=0 ; i < _ntrst ; ++i) { + + lcio::TrackState* ts = tsV[i] ; + _tsloc[ i ] = ts->getLocation() ; + _tsdze[ i ] = ts->getD0() ; + _tsphi[ i ] = ts->getPhi() ; + _tsome[ i ] = ts->getOmega() ; + _tszze[ i ] = ts->getZ0() ; + _tstnl[ i ] = ts->getTanLambda() ; + _tsrpx[ i ] = ts->getReferencePoint()[0] ; + _tsrpy[ i ] = ts->getReferencePoint()[1] ; + _tsrpz[ i ] = ts->getReferencePoint()[2] ; + + for (int j=0;j<15;++j) { + _tscov[ i ][ j ] = ts->getCovMatrix()[j] ; + } } - } - - + //------ fill the Tracks particle ---------------------------- - for(int i=0 ; i < _ntrk ; ++i){ + for (int i=0 ; i < _ntrk ; ++i) { lcio::Track* trk = static_cast( col->getElementAt(i) ) ; - _trori[i] = trk->ext(); - _trtyp[ i ] = trk->getType() ; _trch2[ i ] = trk->getChi2() ; _trndf[ i ] = trk->getNdf() ; _tredx[ i ] = trk->getdEdx() ; _trede[ i ] = trk->getdEdxError() ; _trrih[ i ] = trk->getRadiusOfInnermostHit() ; - _trnts[ i ] = trk->getTrackStates().size() ; - _trfts[ i ] = ( _trnts[i]>0 ? trk->getTrackStates()[0]->ext() -1 : -1 ) ; const lcio::TrackState* ts ; - ts = trk->getTrackState( lcio::TrackState::AtIP ) ; - _trsip[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; - - ts = trk->getTrackState( lcio::TrackState::AtFirstHit ) ; - _trsfh[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; - - ts = trk->getTrackState( lcio::TrackState::AtLastHit ) ; - _trslh[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; - - ts = trk->getTrackState( lcio::TrackState::AtCalorimeter ) ; - _trsca[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; - - - int nshn = ( trk->getSubdetectorHitNumbers().size() < 12 ? - trk->getSubdetectorHitNumbers().size() : 12 ) ; - - for( int j=0; jgetSubdetectorHitNumbers()[j] ; + // save parameters AtIP state + if ( ts ) { + _trome[ i ] = ts->getOmega(); + _trtnl[ i ] = ts->getTanLambda(); + _trthe[ i ] = TMath::PiOver2()-atan(ts->getTanLambda()); + _trdze[ i ] = ts->getD0() ; + _trzze[ i ] = ts->getZ0() ; + _trphi[ i ] = ts->getPhi() ; + } else { + _trome[ i ] = -1; + _trtnl[ i ] = -1; + _trthe[ i ] = -1; + _trdze[ i ] = -1; + _trzze[ i ] = -1; + _trphi[ i ] = -1; + } + + // save sigma values from the covariance matrix + /* + for(int j=0;j<15;++j){ + _trcov[ i ][ j ] = ts->getCovMatrix()[j] ; + } + */ + // save only diagonal values + _trk_sigmal0[i] = ts->getCovMatrix()[0]; + _trk_sigmal1[i] = ts->getCovMatrix()[2]; + _trk_sigmaphi[i] = ts->getCovMatrix()[5]; + _trk_sigmatheta[i] = ts->getCovMatrix()[9]; + _trk_sigmaqoverp[i] = ts->getCovMatrix()[14]; + + if (_writeTrkStatesParameters) { + _trfts[ i ] = ( _trnts[i]>0 ? trk->getTrackStates()[0]->ext() -1 : -1 ) ; + _trsip[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; + + ts = trk->getTrackState( lcio::TrackState::AtFirstHit ) ; + _trsfh[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; + + ts = trk->getTrackState( lcio::TrackState::AtLastHit ) ; + _trslh[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; + + ts = trk->getTrackState( lcio::TrackState::AtCalorimeter ) ; + _trsca[ i ] = ( ts ? ts->ext() - 1 : -1 ) ; + } + + const EVENT::IntVec& subdetectorHits = trk->getSubdetectorHitNumbers(); + if ( trk->getSubdetectorHitNumbers().size() > 2 ) + _trtvhn[ i ] = subdetectorHits[1]+subdetectorHits[2]; // vertex + if ( trk->getSubdetectorHitNumbers().size() > 4 ) + _trtihn[ i ] = subdetectorHits[3]+subdetectorHits[4]; // inner + if ( trk->getSubdetectorHitNumbers().size() > 6 ) + _trtohn[ i ] = subdetectorHits[5]+subdetectorHits[6]; // outer - _trthn[ i ] = trk->getTrackerHits().size(); + _trtnh[ i ] = trk->getNholes(); + _trthn[ i ] = trk->getTrackerHits().size(); + _trtout[ i ] = trk->getTrackerHits().size() - 0.5 * trk->getNdf(); + + if(_writeTrkHitsParameters) { - for( unsigned int ihit=0; ihitgetTrackerHits().size() ; ++ihit ){ - int hit_index = ( trk->getTrackerHits().at(ihit) ? - trk->getTrackerHits().at(ihit)->ext() - 1 : - -1 ); - unsigned det = ( trk->getTrackerHits().at(ihit) ? + int nshn = ( subdetectorHits.size() < 12 ? subdetectorHits.size() : 12); + + for( int j=0; jgetTrackerHits().size() ; ++ihit ){ + int hit_index = ( trk->getTrackerHits().at(ihit) ? + trk->getTrackerHits().at(ihit)->ext() - 1 : -1 ); + unsigned det = ( trk->getTrackerHits().at(ihit) ? (unsigned) (trk->getTrackerHits().at(ihit)->getCellID0() & 0x1f) : 0 ); - _trthd[ i ][ ihit ] = det ; - _trthi[ i ][ ihit ] = hit_index; - + _trthd[ i ][ ihit ] = det ; + _trthi[ i ][ ihit ] = hit_index; + if ( trk->getTrackerHits().at(ihit) ) { + lcio::TrackerHit* hit = static_cast( trk->getTrackerHits().at(ihit) ); + std::cout << " POS " << hit->getPosition( )[0] << std::endl; + _trthx[ i ][ ihit ] = hit->getPosition()[0]; + _trthy[ i ][ ihit ] = hit->getPosition()[1]; + _trthz[ i ][ ihit ] = hit->getPosition()[2]; + } + } } - } -} +} \ No newline at end of file diff --git a/src/VertexBranches.cc b/src/VertexBranches.cc index 26e9c4b..76924ef 100644 --- a/src/VertexBranches.cc +++ b/src/VertexBranches.cc @@ -31,6 +31,7 @@ void VertexBranches::initBranches( TTree* tree, const std::string& pre){ tree->Branch( (pre+"vtprb").c_str() , _vtprb , (pre+"vtprb["+pre+"nvt]/F").c_str() ) ; tree->Branch( (pre+"vtcov").c_str() , _vtcov , (pre+"vtcov["+pre+"nvt][6]/F").c_str() ) ; tree->Branch( (pre+"vtpar").c_str() , _vtpar , (pre+"vtpar["+pre+"nvt][6]/F").c_str() ) ; + tree->Branch( (pre+"vttrchi").c_str() , _vttrchi , (pre+"vttrchi["+pre+"nvt][30]/F").c_str() ) ; } @@ -60,6 +61,12 @@ void VertexBranches::fill(const EVENT::LCCollection* col, EVENT::LCEvent* evt ){ _vtprb[i] = vtx->getProbability(); for(int j=0; j<6 ; ++j) _vtcov[ i ][ j ] = vtx->getCovMatrix()[j] ; // for(int j=0; j<6 ; ++j) _vtpar[ i ][ j ] = vtx->getParameters()[j] ; + + // additional parameters: chi2 of each track + EVENT::FloatVec tracksChi = vtx->getParameters(); + for(unsigned int j = 0; j < tracksChi.size(); j++){ + _vttrchi[i][j] = tracksChi[j]; + } } }