diff --git a/Reconstruction/ParticleID/src/TofRecAlg.cpp b/Reconstruction/ParticleID/src/TofRecAlg.cpp
index 3dce6848fa43d070afabbcf9a1bf9e196a66840d..f7509088083a327a10a69006b04a118c6da350c7 100644
--- a/Reconstruction/ParticleID/src/TofRecAlg.cpp
+++ b/Reconstruction/ParticleID/src/TofRecAlg.cpp
@@ -11,70 +11,33 @@
#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
#include "CLHEP/Units/SystemOfUnits.h"
-#include <math.h>
+#include <cmath>
#include "UTIL/ILDConf.h"
#include "DetIdentifier/CEPCConf.h"
#include "UTIL/CellIDEncoder.h"
+using namespace edm4hep;
+
DECLARE_COMPONENT( TofRecAlg )
//------------------------------------------------------------------------------
TofRecAlg::TofRecAlg( const std::string& name, ISvcLocator* pSvcLocator )
: Algorithm( name, pSvcLocator ) {
- declareProperty("MCParticleCollection", _inMCColHdl, "Handle of the Input MCParticle collection");
- declareProperty("TrackCollection", _inTrackColHdl, "Handle of the Input Track collection from CEPCSW");
+ // input
+ declareProperty("CompleteTracks", m_FulTrkCol, "handler of the input track particle association collection");
+ // output
+ declareProperty("RecTofCollection", m_RecToFCol, "handler of the collection of reconstructed tof");
}
//------------------------------------------------------------------------------
StatusCode TofRecAlg::initialize(){
- info() << "Booking Ntuple" << endmsg;
-
- NTuplePtr nt1(ntupleSvc(), "MyTuples/ToF");
- if ( !nt1 ) {
- m_tuple = ntupleSvc()->book("MyTuples/ToF",CLID_ColumnWiseTuple,"ToF result");
- if ( 0 != m_tuple ) {
- m_tuple->addItem ("ntrk", m_nTracks, 0, 500 ).ignore();
- m_tuple->addIndexedItem("ToFt", m_nTracks, m_ToFt).ignore();
- m_tuple->addIndexedItem("ToFterr", m_nTracks, m_ToFterr).ignore();
- m_tuple->addIndexedItem("ToFx", m_nTracks, m_ToFx).ignore();
- m_tuple->addIndexedItem("ToFy", m_nTracks, m_ToFy).ignore();
- m_tuple->addIndexedItem("ToFz", m_nTracks, m_ToFz).ignore();
- m_tuple->addIndexedItem("length", m_nTracks, m_length).ignore();
- m_tuple->addIndexedItem("lengtherr", m_nTracks, m_lengtherr).ignore();
- m_tuple->addIndexedItem("exp_time", m_nTracks, 5, m_exp_time).ignore();
- m_tuple->addIndexedItem("exp_timeerr", m_nTracks, 5, m_exp_timeerr).ignore();//above need to be stored in edm4hep
-
- m_tuple->addIndexedItem("trackState0", m_nTracks, 15, trackState0).ignore();
- m_tuple->addIndexedItem("trackState1", m_nTracks, 15, trackState1).ignore();
- m_tuple->addIndexedItem("trackState2", m_nTracks, 15, trackState2).ignore();
- m_tuple->addIndexedItem("trackState3", m_nTracks, 15, trackState3).ignore();
- m_tuple->addIndexedItem("trackState4", m_nTracks, 15, trackState4).ignore();
-
- m_tuple->addIndexedItem("p_atIP", m_nTracks, m_p_atIP).ignore();
- m_tuple->addIndexedItem("perr_atIP", m_nTracks, m_perr_atIP).ignore();
- m_tuple->addIndexedItem("p_atLastHit", m_nTracks, m_p_atLast).ignore();
- m_tuple->addIndexedItem("perr_atLastHit", m_nTracks, m_perr_atLast).ignore();
- }
- else { // did not manage to book the N tuple....
- fatal() << "Cannot book MyTuples/ToF" <<endmsg;
- return StatusCode::FAILURE;
- }
+ debug() << "Initializing..." << endmsg;
+ if (m_method == "Simple") {
+ m_pid_svc = service("SimplePIDSvc");
}
- else{
- m_tuple = nt1;
+ else {
+ m_pid_svc = nullptr;
}
-
- auto geomSvc = service<IGeomSvc>("GeomSvc");
- if(geomSvc){
- const dd4hep::Direction& field = geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
- m_field = field.z()/dd4hep::tesla;
- info() << "Magnetic field will obtain from GeomSvc = " << m_field << " tesla" << endmsg;
- }
- else{
- info() << "Failed to find GeomSvc ..." << endmsg;
- info() << "Magnetic field will use what input through python option for this algorithm namse as Field, now " << m_field << " tesla" << endmsg;
- }
-
_nEvt = 0;
return StatusCode::SUCCESS;
}
@@ -82,234 +45,116 @@ StatusCode TofRecAlg::initialize(){
//------------------------------------------------------------------------------
StatusCode TofRecAlg::execute(){
- const edm4hep::TrackCollection* trackCols = nullptr;
- auto rectofCol = m_rectofCol.createAndPut();
-
+ const edm4hep::TrackCollection* fultrkcol = nullptr;
+ auto rectofcol = m_RecToFCol.createAndPut();
+
try {
- trackCols = _inTrackColHdl.get();
+ fultrkcol = m_FulTrkCol.get();
}
catch ( GaudiException &e ) {
- debug() << "Collection " << _inTrackColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ debug() << "Complete track particle association collection " << m_FulTrkCol.fullKey() << " is unavailable in event " << _nEvt << endmsg;
+ _nEvt++;
return StatusCode::SUCCESS;
}
-
-
- if(trackCols){
-
- m_nTracks = 0;
-
- for(auto track : *trackCols){
- // only take the first track
- // if(m_nTracks>=1)break;
-
- int NHits = track.trackerHits_size();
- int det_id = 0;
- bool hasSETHit = false;
- bool hasFTDHit = false;
- UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string );
- for ( int hit=0; hit<NHits; hit++ ) {
- edm4hep::TrackerHit SimTHit = track.getTrackerHits(hit);
- encoder.setValue(SimTHit.getCellID());
- det_id = encoder[lcio::ILDCellID0::subdet];
- if ( det_id == lcio::ILDDetID::SET ){
- hasSETHit = true;
- m_ToFt[m_nTracks] = SimTHit.getTime();
- m_ToFx[m_nTracks] = SimTHit.getPosition()[0];
- m_ToFy[m_nTracks] = SimTHit.getPosition()[1];
- m_ToFz[m_nTracks] = SimTHit.getPosition()[2];
- }
- else if( det_id == lcio::ILDDetID::FTD ){
- const int cellID = SimTHit.getCellID();
- encoder.setValue(cellID);
- int layer = encoder[lcio::ILDCellID0::layer];
- if(layer==4){
- m_ToFt[m_nTracks] = SimTHit.getTime();
- m_ToFx[m_nTracks] = SimTHit.getPosition()[0];
- m_ToFy[m_nTracks] = SimTHit.getPosition()[1];
- m_ToFz[m_nTracks] = SimTHit.getPosition()[2];
- hasFTDHit = true;
- }else{
- m_ToFt[m_nTracks] = -999;
- m_ToFx[m_nTracks] = -999;
- m_ToFy[m_nTracks] = -999;
- m_ToFz[m_nTracks] = -999;
- }
- }
- else {
- m_ToFt[m_nTracks] = -999;
- m_ToFx[m_nTracks] = -999;
- m_ToFy[m_nTracks] = -999;
- m_ToFz[m_nTracks] = -999;
- }
- }//end track hits
- if ( !hasSETHit && !hasFTDHit ) continue;
- auto rectof = rectofCol->create();
- std::array<float, 5> tmp_timeExp_list{ -99.9, -99.9, -99.9, -99.9, -99.9 };
+ hasFTDHit = false;
+ hasSETHit = false;
+
+ for(auto track : *fultrkcol){
- // smearing the measured time by gaussian 0.05 for intrinsic and 0.02 for bunchcrossing, giving the measured terr a constant
- m_ToFterr[m_nTracks] = std::sqrt(tof_resolution*tof_resolution + bunchcrossing*bunchcrossing);
- m_ToFt[m_nTracks] += normal_distribution_tof(generator_tof);
- m_ToFt[m_nTracks] += normal_distribution_bunch(generator_bunch);
+ FindToFHits( track, hasFTDHit, hasSETHit, Toft, Tofx, Tofy, Tofz );
- // found ToF hit, then to evaluate trajectory length
- for( std::vector<edm4hep::TrackState>::const_iterator it = track.trackStates_begin(); it != track.trackStates_end(); it++ ){
-
- edm4hep::TrackState trackState = *it;
- int _location = trackState.location;
- (*trackStateMap[_location])[m_nTracks][0] = _location;
- (*trackStateMap[_location])[m_nTracks][1] = trackState.D0;
- (*trackStateMap[_location])[m_nTracks][2] = trackState.phi;
- (*trackStateMap[_location])[m_nTracks][3] = trackState.omega;
- (*trackStateMap[_location])[m_nTracks][4] = trackState.Z0;
- (*trackStateMap[_location])[m_nTracks][5] = trackState.tanLambda;
- (*trackStateMap[_location])[m_nTracks][6] = trackState.time;
- (*trackStateMap[_location])[m_nTracks][7] = trackState.referencePoint[0];
- (*trackStateMap[_location])[m_nTracks][8] = trackState.referencePoint[1];
- (*trackStateMap[_location])[m_nTracks][9] = trackState.referencePoint[2];
- (*trackStateMap[_location])[m_nTracks][10] = std::sqrt(trackState.covMatrix[0]);
- (*trackStateMap[_location])[m_nTracks][11] = std::sqrt(trackState.covMatrix[2]);
- (*trackStateMap[_location])[m_nTracks][12] = std::sqrt(trackState.covMatrix[5]);
- (*trackStateMap[_location])[m_nTracks][13] = std::sqrt(trackState.covMatrix[9]);
- (*trackStateMap[_location])[m_nTracks][14] = std::sqrt(trackState.covMatrix[14]);
-
+ if ( !hasSETHit && !hasFTDHit ) {
+ debug() << "No SET or FTD hit found for this track" << endmsg;
+ continue;
}
- //expected time calculation
- if ( hasSETHit || hasFTDHit ){
-
- // trackState at IP
- float phi0 = (*trackStateMap[1])[m_nTracks][2];
- if (phi0 < 0) phi0 += 2*M_PI;
- float phi0err = (*trackStateMap[1])[m_nTracks][11];
- //float z0 = (*trackStateMap[1])[m_nTracks][4];
- //float z0err = (*trackStateMap[1])[m_nTracks][13];
- float omega0 = (*trackStateMap[1])[m_nTracks][3];
- float omega0err = (*trackStateMap[1])[m_nTracks][12];
- float totz0 = (*trackStateMap[1])[m_nTracks][4] + (*trackStateMap[1])[m_nTracks][9];
- float totz0err = (*trackStateMap[1])[m_nTracks][13];
- float tanL0 = (*trackStateMap[1])[m_nTracks][5];
- float tanL0err = (*trackStateMap[1])[m_nTracks][14];
-
- float radius0 = 1./std::abs(omega0);
- float radius0err = omega0err/(omega0*omega0);
-
- float pT_tmp = radius0 * m_field * c * 10E-13;
- float pTerr_tmp = radius0err * m_field * c * 10E-13;
-
- float pZ_tmp = pT_tmp*tanL0;
- float pZerr_tmp = std::sqrt( pT_tmp*pT_tmp*tanL0err*tanL0err + tanL0*tanL0*pTerr_tmp*pTerr_tmp );
-
- float p_atIP = std::sqrt( pT_tmp*pT_tmp + pZ_tmp*pZ_tmp );
- float perr_atIP = std::sqrt( (pT_tmp*pT_tmp / p_atIP*p_atIP)*(pTerr_tmp*pTerr_tmp) + (pZ_tmp*pZ_tmp / p_atIP*p_atIP)*(pZerr_tmp*pZerr_tmp) );
-
- // trackState at the last hit
- float phi1 = (*trackStateMap[3])[m_nTracks][2];
- if (phi1 < 0) phi1 += 2*M_PI;
- float phi1err = (*trackStateMap[3])[m_nTracks][11];
- //float z1 = (*trackStateMap[3])[m_nTracks][4];
- //float z1err = (*trackStateMap[3])[m_nTracks][13];
- float omega1 = (*trackStateMap[3])[m_nTracks][3];
- float omega1err = (*trackStateMap[3])[m_nTracks][12];
- float totz1 = (*trackStateMap[3])[m_nTracks][4] + (*trackStateMap[3])[m_nTracks][9];
- float totz1err = (*trackStateMap[3])[m_nTracks][13];
- float tanL1 = (*trackStateMap[3])[m_nTracks][5];
- float tanL1err = (*trackStateMap[3])[m_nTracks][14];
- float positionx = (*trackStateMap[3])[m_nTracks][7];
- float positiony = (*trackStateMap[3])[m_nTracks][8];
- float positionz = (*trackStateMap[3])[m_nTracks][9];
-
- float radius1 = 1./std::abs(omega1);
- float radius1err = omega1err/(omega1*omega1);
-
- float pT_tmp1 = radius1 * m_field * c * 10E-13;
- float pTerr_tmp1 = radius1err * m_field * c * 10E-13;
-
- float pZ_tmp1 = pT_tmp1*tanL1;
- float pZerr_tmp1 = std::sqrt( pT_tmp1*pT_tmp1*tanL1err*tanL1err + tanL1*tanL1*pTerr_tmp1*pTerr_tmp1 );
-
- float p_atLast = std::sqrt(pT_tmp1*pT_tmp1 + pZ_tmp1*pZ_tmp1);
- float perr_atLast = std::sqrt( (pT_tmp1*pT_tmp1 / p_atLast*p_atLast)*(pTerr_tmp1*pTerr_tmp1) + (pZ_tmp1*pZ_tmp1 / p_atLast*p_atLast)*(pZerr_tmp1*pZerr_tmp1) );
-
- //assum all tracks end in one circle, for multi cirlce tracks, the result does not make sense
- float deltaphi = phi1 -phi0;
- if ( omega0 * deltaphi > 0.0 ){
- deltaphi = 2. * M_PI - std::abs( deltaphi );
- }
-
- // calculate the length
- float _length = std::sqrt( radius0*radius0 * (deltaphi)*(deltaphi) + (totz1-totz0)*(totz1-totz0) );//in mm
- float _lengtherr = 0.0;
- float _ldr2 = (radius0*radius0) * (deltaphi)*(deltaphi)*(deltaphi)*(deltaphi) / (_length*_length);
- float _ldphi02 = radius0*radius0*radius0*radius0 * (deltaphi)*(deltaphi) / (_length*_length);
- float _ldphi12 = _ldphi02;
- float _ldz02 = (totz1-totz0)*(totz1-totz0) / (_length*_length);
- float _ldz12 = _ldz02;
-
- _lengtherr = std::sqrt( _ldr2*radius0err*radius0err + _ldphi02*phi0err*phi0err + _ldphi12*phi1err*phi1err + _ldz02*totz0err*totz0err + _ldz12*totz1err*totz1err );
- m_length[m_nTracks] = _length;
- m_lengtherr[m_nTracks] = _lengtherr;
- m_p_atIP[m_nTracks] = p_atIP;
- m_perr_atIP[m_nTracks] = perr_atIP;
- m_p_atLast[m_nTracks] = p_atLast;
- m_perr_atLast[m_nTracks] = perr_atLast;
-
- //momenta to kg m/s
- float p_tmp = p_atIP * GeV2kgms;
- float perr_tmp = perr_atIP * GeV2kgms;
-
- for (auto it = masses.begin(); it != masses.end(); ++it){//expected times and errors for e, mu, pi, k, proton
- int idx = it->first;
- double assum = it->second;
- double m_tmp = assum * GeV2kg;
- double E_tmp = std::sqrt(p_tmp * p_tmp * c * c + m_tmp * m_tmp * c * c * c * c);
- double v_tmp = p_tmp * c * c / E_tmp;
- double verr_tmp = perr_tmp * ( c * c * E_tmp - (p_tmp * p_tmp * c * c * c * c/E_tmp))/ (E_tmp * E_tmp);
- double t_tmp = 10e5 * _length / v_tmp;
- double terr_tmp = 10e5 * std::sqrt( (_lengtherr/v_tmp)*(_lengtherr/v_tmp) + ( (_length*_length)/(v_tmp*v_tmp*v_tmp*v_tmp) )*verr_tmp*verr_tmp );
- m_exp_time[m_nTracks][idx] = t_tmp;
- m_exp_timeerr[m_nTracks][idx] = terr_tmp;
- tmp_timeExp_list[idx] = (float) t_tmp;
- debug()<<"idx : "<< idx <<" assum: "<< assum <<" omega: "<<omega0<<" length: "<<_length<<" phi1: "<<phi1<<" phi0: "<<phi0<<" phi1-phi0: "<<deltaphi<<" exp_t: "<<t_tmp<<" obs_t: "<<m_ToFt[m_nTracks]<<endmsg;
- }
- rectof.setTimeExp( tmp_timeExp_list );
- rectof.setTime( m_ToFt[m_nTracks] );
- rectof.setSigma( m_ToFterr[m_nTracks] );
- rectof.setPosition( { positionx, positiony, positionz } );
- rectof.setPathLength( { _length, _length, _length, _length, _length } );
-
- }else{
- m_length[m_nTracks] = -999;
- m_lengtherr[m_nTracks] = -999;
- m_p_atIP[m_nTracks] = -999;
- m_perr_atIP[m_nTracks] = -999;
- m_p_atLast[m_nTracks] = -999;
- m_perr_atLast[m_nTracks] = -999;
- for (auto it = masses.begin(); it != masses.end(); ++it){
- int idx = it->first;
- m_exp_time[m_nTracks][idx] = -999;
- m_exp_timeerr[m_nTracks][idx] = -999;
+ for( std::vector<edm4hep::TrackState>::const_iterator it = track.trackStates_begin(); it != track.trackStates_end(); it++ ){
+ edm4hep::TrackState trackState = *it;
+ if ( trackState.location == TrackState::AtIP ){
+ position0x = trackState.referencePoint[0];
+ position0y = trackState.referencePoint[1];
+ position0z = trackState.referencePoint[2];
+ debug()<<"position0x = "<<position0x<<" position0y = "<<position0y<<" position0z = "<<position0z<<endmsg;
+ first_hit_pos = {position0x, position0y, position0z};//for trajectory length calculation
}
- rectof.setTimeExp( tmp_timeExp_list );
- rectof.setTime( -999 );
- rectof.setSigma( -999 );
- rectof.setPosition( { -999, -999, -999 } );
- rectof.setPathLength( { -999, -999, -999, -999, -999 } );
- }
- m_nTracks++;
-
- }// end track iteration
-
- }//end track collection
+ if ( trackState.location == TrackState::AtLastHit ){
+ position1x = trackState.referencePoint[0];
+ position1y = trackState.referencePoint[1];
+ position1z = trackState.referencePoint[2];
+ debug()<<"position1x = "<<position1x<<" position1y = "<<position1y<<" position1z = "<<position1z<<endmsg;
+ last_hit_pos = {position1x, position1y, position1z};
+ }
+ }// find track states at IP and OTK
- m_tuple->write();
+ // smearing the measured time by gaussian 0.05 for intrinsic and 0.02 for bunchcrossing, giving the measured terr a constant
+ Tofterr = std::sqrt( tof_resolution * tof_resolution + bunchcrossing * bunchcrossing );
+ Toft += normal_distribution_tof( generator_tof );
+ Toft += normal_distribution_bunch( generator_bunch );
+
+ // length of the trajectory using simple pid svc
+ double _length, _p_atIP, _costheta;
+ m_pid_svc->getTrackPars( first_hit_pos, last_hit_pos, track, TrackState::AtIP, _length, _p_atIP, _costheta );
+ debug()<<"_length = "<<_length<<" _p_atIP = "<<_p_atIP<<" _costheta = "<<_costheta<<endmsg;
+
+ std::array<float, 5> exp_times = { -99.9, -99.9, -99.9, -99.9, -99.9 };
+ // evaluate expected times and errors for e, mu, pi, k, proton hypotheses
+ for (auto it = masses.begin(); it != masses.end(); ++it){
+ int idx = it->first;
+ double assum = it->second;//mass in GeV
+ double energy = std::sqrt( _p_atIP * _p_atIP + assum * assum );
+ double velocity = ( _p_atIP / energy ) * c_mm_ns;//mm/ns
+ double t_exp = _length / velocity;//ns
+ exp_times[idx] = t_exp;
+ debug()<<" idx : "<< idx <<" assum: "<< assum <<" length: "<<_length<<" exp_t: "<<t_exp<<" obs_t: "<<Toft<<endmsg;
+ }//end loop over masses
+
+ auto rectof = rectofcol->create();
+ rectof.setTime(Toft);
+ rectof.setTimeExp( exp_times );
+ rectof.setPosition({Tofx, Tofy, Tofz});
+ rectof.setSigma(Tofterr);
+ float _flength = float(_length);//avoid the warning of narrowing conversion
+ rectof.setPathLength({_flength, _flength, _flength, _flength, _flength});//all particles have the same path length for now, becasue the path length is calculated by the simple pid svc with helix.
+ rectof.setTrack(track);
+
+ }//end loop over tracks
_nEvt++;
return StatusCode::SUCCESS;
-}
+}// end execute
//------------------------------------------------------------------------------
StatusCode TofRecAlg::finalize(){
debug() << "Finalizing..." << endmsg;
-
return StatusCode::SUCCESS;
}
+
+void TofRecAlg::FindToFHits(const edm4hep::Track& _track, bool& _hasFTDHit, bool& _hasSETHit, double& _Toft, double& _Tofx, double& _Tofy, double& _Tofz){
+ int NHits = _track.trackerHits_size();
+ int det_id = 0;
+ UTIL::BitField64 encoder( lcio::ILDCellID0::encoder_string );
+ for ( int hit=0; hit<NHits; hit++ ) {
+ edm4hep::TrackerHit SimTHit = _track.getTrackerHits(hit);
+ encoder.setValue(SimTHit.getCellID());
+ det_id = encoder[lcio::ILDCellID0::subdet];
+ if ( det_id == lcio::ILDDetID::SET ){
+ _hasSETHit = true;
+ _Toft = SimTHit.getTime();
+ _Tofx = SimTHit.getPosition()[0];
+ _Tofy = SimTHit.getPosition()[1];
+ _Tofz = SimTHit.getPosition()[2];
+ }
+ else if( det_id == lcio::ILDDetID::FTD ){
+ const int cellID = SimTHit.getCellID();
+ encoder.setValue(cellID);
+ int layer = encoder[lcio::ILDCellID0::layer];
+ if(layer==4){
+ _Toft = SimTHit.getTime();
+ _Tofx = SimTHit.getPosition()[0];
+ _Tofy = SimTHit.getPosition()[1];
+ _Tofz = SimTHit.getPosition()[2];
+ _hasFTDHit = true;
+ }//find ftd hit
+ }//find ftd hit
+ }//end track hits loop
+}
\ No newline at end of file
diff --git a/Reconstruction/ParticleID/src/TofRecAlg.h b/Reconstruction/ParticleID/src/TofRecAlg.h
index 1ddd298f0be5415a1babc57deec3a51393513f76..9d2bce99b1cbcc4d551e7494926f554b4840e3cc 100644
--- a/Reconstruction/ParticleID/src/TofRecAlg.h
+++ b/Reconstruction/ParticleID/src/TofRecAlg.h
@@ -5,12 +5,20 @@
#include "GaudiKernel/Algorithm.h"
#include "edm4hep/MCParticleCollection.h"
-#include "edm4hep/TrackCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"
+
+#include "edm4hep/TrackerHit.h"
+#include "edm4hep/TrackCollection.h"
#include "edm4hep/TrackerHitCollection.h"
-#include <random>
-#include "GaudiKernel/NTuple.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/RecDqdx.h"
+#include "edm4hep/RecDqdxCollection.h"
#include "edm4cepc/RecTofCollection.h"
+#include "edm4hep/TrackState.h"
+#include "SimplePIDSvc/ISimplePIDSvc.h"
+#include "edm4hep/Vector3d.h"
+#include <random>
class TofRecAlg : public Algorithm {
public:
@@ -23,64 +31,17 @@ class TofRecAlg : public Algorithm {
StatusCode finalize() override;
private:
- DataHandle<edm4hep::MCParticleCollection> _inMCColHdl{"MCParticle", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackCollection> _inTrackColHdl{"CompleteTracks", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::RecTofCollection> m_rectofCol{"RecTofCollection", Gaudi::DataHandle::Writer, this};
-
- Gaudi::Property<double> m_field{this, "Field", 3.0};
+ DataHandle<edm4hep::TrackCollection> m_FulTrkCol{"CompleteTracks", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::RecTofCollection> m_RecToFCol{"RecTofCollection", Gaudi::DataHandle::Writer, this};
+ Gaudi::Property<std::string> m_method{this, "Method", "Simple"};
+ SmartIF<ISimplePIDSvc> m_pid_svc;
- NTuple::Tuple* m_tuple;
- NTuple::Item<long> m_nTracks;
- NTuple::Array<float> m_x;
- NTuple::Array<float> m_y;
- NTuple::Array<float> m_z;
- NTuple::Array<float> m_px;
- NTuple::Array<float> m_py;
- NTuple::Array<float> m_pz;
- NTuple::Array<float> m_p;
- NTuple::Array<float> m_d0;
- NTuple::Array<float> m_phi0;
- NTuple::Array<float> m_omega;
- NTuple::Array<float> m_z0;
- NTuple::Array<float> m_tanLambda;
- NTuple::Array<float> m_sigma_d0;
- NTuple::Array<float> m_sigma_phi0;
- NTuple::Array<float> m_sigma_omega;
- NTuple::Array<float> m_sigma_z0;
- NTuple::Array<float> m_sigma_tanLambda;
+ void FindToFHits(const edm4hep::Track& _track, bool& _hasFTDHit, bool& _hasSETHit, double& _Toft, double& _Tofx, double& _Tofy, double& _Tofz);
- NTuple::Array<float> m_ToFt;
- NTuple::Array<float> m_ToFterr;
- NTuple::Array<float> m_ToFx;
- NTuple::Array<float> m_ToFy;
- NTuple::Array<float> m_ToFz;
-
- NTuple::Matrix<float> trackState0;
- NTuple::Matrix<float> trackState1;
- NTuple::Matrix<float> trackState2;
- NTuple::Matrix<float> trackState3;
- NTuple::Matrix<float> trackState4;
- NTuple::Matrix<float> trackState5;
-
- std::map<int, NTuple::Matrix<float>*> trackStateMap = {
- {0, &trackState0},
- {1, &trackState1},
- {2, &trackState2},
- {3, &trackState3},
- {4, &trackState4},
- {5, &trackState5}
- };
-
- NTuple::Array<float> m_length;
- NTuple::Array<float> m_lengtherr;
-
- NTuple::Array<float> m_p_atIP;
- NTuple::Array<float> m_perr_atIP;
- NTuple::Array<float> m_p_atLast;
- NTuple::Array<float> m_perr_atLast;
-
- NTuple::Matrix<double> m_exp_time;
- NTuple::Matrix<double> m_exp_timeerr;
+ double Tofx, Tofy, Tofz, Toft, Tofterr;
+ double position0x, position0y, position0z, position1x, position1y, position1z;
+ edm4hep::Vector3d first_hit_pos;
+ edm4hep::Vector3d last_hit_pos;
int _nEvt;
@@ -92,10 +53,10 @@ class TofRecAlg : public Algorithm {
{4, 0.938272},
};
- float GeV2kgms = 5.344286e-19;//momenta in GeV to kg m/s
- float GeV2kg = 1.78266192e-27;//mass in GeV to kg
- float c = 2.99792458e8;//spead of light in m/s
- float q = 1.60217662e-19;//electron charge in C
+ float c_mm_ns = 2.99792458e2;//spead of light in mm/ns
+
+ bool hasFTDHit = false;
+ bool hasSETHit = false;
float bunchcrossing = 0.02;//ns
float tof_resolution = 0.05;//ns
@@ -103,8 +64,5 @@ class TofRecAlg : public Algorithm {
std::normal_distribution<float> normal_distribution_tof{0, tof_resolution};
std::default_random_engine generator_bunch;
std::normal_distribution<float> normal_distribution_bunch{0, bunchcrossing};
-
-
};
-
-#endif
+#endif
\ No newline at end of file