From 1b1471ea7c8b34f84954fe2367a2ed533a31f521 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Fri, 12 Aug 2022 18:37:12 -0400 Subject: [PATCH 1/6] fix: update `IRTAlgorithm` w.r.t. data model changes --- JugPID/src/components/IRTAlgorithm.cpp | 48 ++++++++++++-------------- JugPID/src/components/IRTAlgorithm.h | 8 ++--- 2 files changed, 27 insertions(+), 29 deletions(-) diff --git a/JugPID/src/components/IRTAlgorithm.cpp b/JugPID/src/components/IRTAlgorithm.cpp index b606d360..9a4d5e3e 100644 --- a/JugPID/src/components/IRTAlgorithm.cpp +++ b/JugPID/src/components/IRTAlgorithm.cpp @@ -161,7 +161,7 @@ StatusCode Jug::PID::IRTAlgorithm::initialize( void ) // Input hit collection; m_inputHitCollection = - std::make_unique>((dname + "Hits").c_str(), + std::make_unique>((dname + "Hits").c_str(), Gaudi::DataHandle::Reader, this); // Output PID info collection; m_outputCherenkovPID = @@ -309,15 +309,15 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // At this point have a reference to a 'mctrack', either this or that way; // Now just follow the logic of TrackParamTruthInit.cpp; - // genStatus = 1 means thrown G4Primary, but dd4gun uses genStatus == 0 - if (mctrack.genStatus() > 1 ) { + // generatorStatus = 1 means thrown G4Primary, but dd4gun uses generatorStatus == 0 + if (mctrack.getGeneratorStatus() > 1 ) { if (msgLevel(MSG::DEBUG)) { - debug() << "ignoring particle with genStatus = " << mctrack.genStatus() << endmsg; + debug() << "ignoring particle with generatorStatus = " << mctrack.getGeneratorStatus() << endmsg; } continue; } //if - const double charge = m_pidSvc->particle(mctrack.pdgID()).charge; + const double charge = m_pidSvc->particle(mctrack.getPDG()).charge; if (abs(charge) < std::numeric_limits::epsilon()) { if (msgLevel(MSG::DEBUG)) { debug() << "ignoring neutral particle" << endmsg; @@ -326,12 +326,12 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) } //if // FIXME: consider only primaries for the time being; - if (mctrack.g4Parent()) continue; - //printf("@@@ %d %3d %d\n", mctrack.ID(), mctrack.pdgID(), mctrack.g4Parent()); + if (mctrack.parents_size()>0) continue; + //printf("@@@ %d %3d %d\n", mctrack.ID(), mctrack.getPDG(), mctrack.parents_size()); // FIXME: a hack; remove muons; - if (mctrack.pdgID() == 13) continue; + if (mctrack.getPDG() == 13) continue; - auto particle = new ChargedParticle(mctrack.pdgID()); + auto particle = new ChargedParticle(mctrack.getPDG()); event->AddChargedParticle(particle); // Start radiator history for all known radiators; @@ -344,16 +344,14 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) //printf("%3d\n", mctracks.size()); for(const auto &hit: hits) { // FIXME: range checks; - const auto &phtrack = mctracks[hit.truth().trackID]; - //const auto &phtrack = mctracks[hit.truth().trackID-1]; + const auto &phtrack = hit.getMCParticle(); // FIXME: yes, use MC truth here; not really needed I guess; // FIXME: range checks; - //printf("Here: %d %d %3d!\n", phtrack.parents()[0], mctrack.ID(), hit.truth().trackID); //if (phtrack.parents()[0] != mctrack.ID()) continue; // Vertex where photon was created; - const auto &vs = phtrack.vs(); + const auto &vs = phtrack.getVertex(); TVector3 vtx(vs.x, vs.y, vs.z); //printf("@@@ %f %f %f\n", vs.x, vs.y, vs.z); { @@ -375,7 +373,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // Simulate QE & geometric sensor efficiency; FIXME: hit.energy() is numerically // in GeV units, but Gaudi::Units::GeV = 1000; prefer to convert photon energies // to [eV] in all places by hand; - double eVenergy = 1E9*hit.energyDeposit(); + double eVenergy = 1E9*hit.getEDep(); //printf("%f\n", eVenergy); if (!QE_pass(eVenergy, m_rngUni()) || m_rngUni() > /*m_GeometricEfficiency.value()**/m_SafetyFactor.value()) { @@ -387,7 +385,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // Add the point to the radiator history buffer; FIXME: all this is not // needed if ACTS trajectory parameterization is used; - particle->FindRadiatorHistory(radiator)->AddStepBufferPoint(phtrack.time(), vtx); + particle->FindRadiatorHistory(radiator)->AddStepBufferPoint(phtrack.getTime(), vtx); //printf("Here-2A\n"); continue; @@ -402,13 +400,13 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // FIXME: '0' - for now assume a single photon detector type for the whole ERICH (DRICH), // which must be a reasonable assumption; eventually may want to generalize; const auto pd = m_IrtDet->m_PhotonDetectors[0]; - uint64_t vcopy = hit.cellID() & m_ReadoutCellMask; + uint64_t vcopy = hit.getCellID() & m_ReadoutCellMask; { const auto irt = pd->GetIRT(vcopy); auto sensor = dynamic_cast(irt->tail()->GetSurface()); double pitch = m_SensorPixelPitch.value(), sens = m_SensorPixelSize.value();//pitch - gap; - const auto &x = hit.position(); + const auto &x = hit.getPosition(); TVector3 lpt(x.x, x.y, x.z); double lxy[2] = {sensor->GetLocalX(lpt), sensor->GetLocalY(lpt)}, lxyPixels[2] = {0.0}; //printf("%7.2f %7.2f\n", lxy[0], lxy[1]); @@ -438,7 +436,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // Start vertex and momentum; photon->SetVertexPosition(vtx); { - const auto &p = phtrack.ps(); + const auto &p = phtrack.getMomentum(); photon->SetVertexMomentum(TVector3(p.x, p.y, p.z)); } @@ -486,7 +484,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // FIXME: need it not at vertex, but in the radiator; as coded here, this can // hardly work once the magnetic field is turned on; but this is a best guess // if nothing else is available; - const auto &p = mctrack.ps(); + const auto &p = mctrack.getMomentum(); p0 = TVector3(p.x, p.y, p.z); @@ -586,7 +584,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) eicd::CherenkovPdgHypothesis hypothesis; // Flip sign if needed; FIXME: use reconstructed one?; - if (m_pidSvc->particle(mctrack.pdgID()).charge < 0.0) pdg_table[ip] *= -1; + if (m_pidSvc->particle(mctrack.getPDG()).charge < 0.0) pdg_table[ip] *= -1; // Storage model is not exactly efficient, but it is simple for users; hypothesis.radiator = ir; @@ -594,7 +592,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) hypothesis.npe = hypo->GetNpe (radiator); hypothesis.weight = hypo->GetWeight(radiator); - cbuffer.addoptions(hypothesis); + cbuffer.addToOptions(hypothesis); } //for ip @@ -641,18 +639,18 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) rdata.rindex = npe ? ri / npe : 0.0; rdata.wavelength = npe ? wl / npe : 0.0; //printf("@@@ %7.2f\n", 1000*rdata.theta); - cbuffer.addangles(rdata); + cbuffer.addToAngles(rdata); } } //for ir // FIXME: well, and what does go here instead of 0?; - cbuffer.ID({0, algorithmID()}); + // cbuffer.ID({0, algorithmID()}); // Reference to either MC track or a reconstructed track; #ifdef _USE_RECONSTRUCTED_TRACKS_ - cbuffer.recID(rctrack.ID()); + cbuffer.recID(rctrack.ID()); // do not use, out of date #else - cbuffer.recID(mctrack.ID()); + //cbuffer.setAssociatedParticle(mctrack); // TODO: wait for https://eicweb.phy.anl.gov/EIC/eicd/-/merge_requests/86 merge #endif } } //for rctrack (mctrack) diff --git a/JugPID/src/components/IRTAlgorithm.h b/JugPID/src/components/IRTAlgorithm.h index d6300caa..e38d1875 100644 --- a/JugPID/src/components/IRTAlgorithm.h +++ b/JugPID/src/components/IRTAlgorithm.h @@ -11,8 +11,8 @@ #include "JugBase/IGeoSvc.h" #include "JugBase/IParticleSvc.h" -#include "dd4pod/Geant4ParticleCollection.h" -#include "dd4pod/TrackerHitCollection.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" #include "eicd/TrajectoryCollection.h" #include "eicd/ReconstructedParticleCollection.h" @@ -62,8 +62,8 @@ namespace Jug::PID { StatusCode finalize( void ) override; // Input collections; - std::unique_ptr> m_inputHitCollection; - DataHandle m_inputMCParticles { + std::unique_ptr> m_inputHitCollection; + DataHandle m_inputMCParticles { "inputMCParticles", Gaudi::DataHandle::Reader, this -- GitLab From 16ea0826615678d4a2fb7ec504432f474e0c84f5 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Fri, 12 Aug 2022 19:25:29 -0400 Subject: [PATCH 2/6] relation `eicd::CherenkovParticleID` to `edm4hep::MCParticle` is now possible --- JugPID/src/components/IRTAlgorithm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/JugPID/src/components/IRTAlgorithm.cpp b/JugPID/src/components/IRTAlgorithm.cpp index 9a4d5e3e..ebf0dc72 100644 --- a/JugPID/src/components/IRTAlgorithm.cpp +++ b/JugPID/src/components/IRTAlgorithm.cpp @@ -650,7 +650,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) #ifdef _USE_RECONSTRUCTED_TRACKS_ cbuffer.recID(rctrack.ID()); // do not use, out of date #else - //cbuffer.setAssociatedParticle(mctrack); // TODO: wait for https://eicweb.phy.anl.gov/EIC/eicd/-/merge_requests/86 merge + cbuffer.setAssociatedParticle(mctrack); #endif } } //for rctrack (mctrack) -- GitLab From 5af163c1fe14957a328589dbdc46acf5fd626ad3 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Sat, 13 Aug 2022 01:35:30 -0400 Subject: [PATCH 3/6] debugging statements for `IRTAlgorithm` --- JugPID/src/components/IRTAlgorithm.cpp | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/JugPID/src/components/IRTAlgorithm.cpp b/JugPID/src/components/IRTAlgorithm.cpp index ebf0dc72..1f4ab1e6 100644 --- a/JugPID/src/components/IRTAlgorithm.cpp +++ b/JugPID/src/components/IRTAlgorithm.cpp @@ -270,6 +270,7 @@ StatusCode Jug::PID::IRTAlgorithm::initialize( void ) StatusCode Jug::PID::IRTAlgorithm::execute( void ) { + printf("####################################### new event\n"); // Input collection(s); const auto &hits = *m_inputHitCollection->get(); const auto &mctracks = *m_inputMCParticles.get(); @@ -335,8 +336,10 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) event->AddChargedParticle(particle); // Start radiator history for all known radiators; - for(auto rptr: m_IrtDet->Radiators()) + for(auto rptr: m_IrtDet->Radiators()) { particle->StartRadiatorHistory(std::make_pair(rptr.second, new RadiatorHistory())); + printf("found radiator %s\n",rptr.second->GetAlternativeMaterialName()); + } // Distribute photons over radiators where they were presumably produced, to // create a structure which would mimic a standalone G4 stepper behavior; @@ -345,6 +348,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) for(const auto &hit: hits) { // FIXME: range checks; const auto &phtrack = hit.getMCParticle(); + printf("here0\n"); // FIXME: yes, use MC truth here; not really needed I guess; // FIXME: range checks; @@ -353,11 +357,13 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // Vertex where photon was created; const auto &vs = phtrack.getVertex(); TVector3 vtx(vs.x, vs.y, vs.z); - //printf("@@@ %f %f %f\n", vs.x, vs.y, vs.z); + printf("@@@ photon vtx: %f %f %f\n", vs.x, vs.y, vs.z); { - //const auto &ps = phtrack.ps(); - //double phi = TVector3(ps.x, ps.y, ps.z).Phi(); - //printf("%7.2f\n", phi); + const auto &ps = phtrack.getMomentum(); + double phi = TVector3(ps.x, ps.y, ps.z).Phi(); + double theta = TVector3(ps.x, ps.y, ps.z).Theta(); + printf("@@@ theta, phi = %7.2f %7.2f\n",theta,phi); + if(phtrack.getPDG()!=-22) fprintf(stderr,"ERROR: this was not an optical photon\n"); //if (fabs(phi - M_PI/2) > M_PI/4 && fabs(phi + M_PI/2) > M_PI/4) continue; } //if (vs.z < 2500 || vs.z > 2700) continue; @@ -365,16 +371,17 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // FIXME: this is in general correct, but needs refinement; TVector3 ip(0,0,0); // FIXME: unify with the code below; + printf("guessing radiator..\n"); auto radiator = m_IrtDet->GuessRadiator(vtx, (vtx - ip).Unit()); // FIXME: do it better later; if (!radiator) continue; - //printf("Here-1\n"); + printf("Here-1\n"); // Simulate QE & geometric sensor efficiency; FIXME: hit.energy() is numerically // in GeV units, but Gaudi::Units::GeV = 1000; prefer to convert photon energies // to [eV] in all places by hand; double eVenergy = 1E9*hit.getEDep(); - //printf("%f\n", eVenergy); + printf("%f\n", eVenergy); if (!QE_pass(eVenergy, m_rngUni()) || m_rngUni() > /*m_GeometricEfficiency.value()**/m_SafetyFactor.value()) { @@ -387,10 +394,10 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // needed if ACTS trajectory parameterization is used; particle->FindRadiatorHistory(radiator)->AddStepBufferPoint(phtrack.getTime(), vtx); - //printf("Here-2A\n"); + printf("Here-2A\n"); continue; } //if - //printf("Here-2B\n"); + printf("Here-2B\n"); //if (vs.z < 2500 || vs.z > 2700) continue; //if (vs.z < 2200 || vs.z > 2400) continue; -- GitLab From 0b44052f8b94ebb3c5117389351101dadad586c5 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Sat, 13 Aug 2022 21:39:00 -0400 Subject: [PATCH 4/6] modify some debugging printouts --- JugPID/src/components/IRTAlgorithm.cpp | 31 ++++++++++++++------------ 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/JugPID/src/components/IRTAlgorithm.cpp b/JugPID/src/components/IRTAlgorithm.cpp index 1f4ab1e6..33048b1b 100644 --- a/JugPID/src/components/IRTAlgorithm.cpp +++ b/JugPID/src/components/IRTAlgorithm.cpp @@ -338,7 +338,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // Start radiator history for all known radiators; for(auto rptr: m_IrtDet->Radiators()) { particle->StartRadiatorHistory(std::make_pair(rptr.second, new RadiatorHistory())); - printf("found radiator %s\n",rptr.second->GetAlternativeMaterialName()); + // printf("found radiator %s\n",rptr.second->GetAlternativeMaterialName()); } // Distribute photons over radiators where they were presumably produced, to @@ -348,7 +348,11 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) for(const auto &hit: hits) { // FIXME: range checks; const auto &phtrack = hit.getMCParticle(); - printf("here0\n"); + // printf("here0\n"); + if(phtrack.getPDG()!=-22) { + fprintf(stderr,"ERROR: this hit was not from an optical photon: PDG = %d\n",phtrack.getPDG()); + continue; + } // FIXME: yes, use MC truth here; not really needed I guess; // FIXME: range checks; @@ -357,31 +361,30 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // Vertex where photon was created; const auto &vs = phtrack.getVertex(); TVector3 vtx(vs.x, vs.y, vs.z); - printf("@@@ photon vtx: %f %f %f\n", vs.x, vs.y, vs.z); + // printf("@@@ photon vtx: %f %f %f\n", vs.x, vs.y, vs.z); { - const auto &ps = phtrack.getMomentum(); - double phi = TVector3(ps.x, ps.y, ps.z).Phi(); - double theta = TVector3(ps.x, ps.y, ps.z).Theta(); - printf("@@@ theta, phi = %7.2f %7.2f\n",theta,phi); - if(phtrack.getPDG()!=-22) fprintf(stderr,"ERROR: this was not an optical photon\n"); - //if (fabs(phi - M_PI/2) > M_PI/4 && fabs(phi + M_PI/2) > M_PI/4) continue; + // const auto &ps = phtrack.getMomentum(); + // double phi = TVector3(ps.x, ps.y, ps.z).Phi(); + // double theta = TVector3(ps.x, ps.y, ps.z).Theta(); + // printf("@@@ theta, phi = %7.2f %7.2f\n",theta,phi); + //if (fabs(phi - M_PI/2) > M_PI/4 && fabs(phi + M_PI/2) > M_PI/4) continue; } //if (vs.z < 2500 || vs.z > 2700) continue; // FIXME: this is in general correct, but needs refinement; TVector3 ip(0,0,0); // FIXME: unify with the code below; - printf("guessing radiator..\n"); + // printf("guessing radiator..\n"); auto radiator = m_IrtDet->GuessRadiator(vtx, (vtx - ip).Unit()); // FIXME: do it better later; if (!radiator) continue; - printf("Here-1\n"); + // printf("Here-1\n"); // Simulate QE & geometric sensor efficiency; FIXME: hit.energy() is numerically // in GeV units, but Gaudi::Units::GeV = 1000; prefer to convert photon energies // to [eV] in all places by hand; double eVenergy = 1E9*hit.getEDep(); - printf("%f\n", eVenergy); + // printf("%f\n", eVenergy); if (!QE_pass(eVenergy, m_rngUni()) || m_rngUni() > /*m_GeometricEfficiency.value()**/m_SafetyFactor.value()) { @@ -394,10 +397,10 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // needed if ACTS trajectory parameterization is used; particle->FindRadiatorHistory(radiator)->AddStepBufferPoint(phtrack.getTime(), vtx); - printf("Here-2A\n"); + // printf("Here-2A\n"); continue; } //if - printf("Here-2B\n"); + // printf("Here-2B\n"); //if (vs.z < 2500 || vs.z > 2700) continue; //if (vs.z < 2200 || vs.z > 2400) continue; -- GitLab From 720aa827cdb8d166ac95ed8f9df6cc8f4b4b0a40 Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Mon, 15 Aug 2022 19:32:30 -0400 Subject: [PATCH 5/6] new file: JugPID/tests/options/irt.py --- JugPID/tests/options/irt.py | 84 +++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) create mode 100644 JugPID/tests/options/irt.py diff --git a/JugPID/tests/options/irt.py b/JugPID/tests/options/irt.py new file mode 100644 index 00000000..657060d6 --- /dev/null +++ b/JugPID/tests/options/irt.py @@ -0,0 +1,84 @@ +import os +from Gaudi.Configuration import * +from GaudiKernel import SystemOfUnits as units + +from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc +from Configurables import PodioInput +from Configurables import Jug__PID__IRTAlgorithm as IRTAlgorithm + +### SETTINGS #################### +compact_file = os.environ['DETECTOR_PATH'] + '/' + os.environ['DETECTOR'] +'_drich_only.xml' +sim_file = os.environ['JUGGLER_SIM_FILE'] +irt_auxfile = os.environ['JUGGLER_IRT_AUXFILE'] +out_file = os.environ['JUGGLER_REC_FILE'] +################################# + +# Well, for now only need DRICH geometry for this example; +geo_service = GeoSvc("GeoSvc", detectors=[compact_file]) + +# Input file after 'npsim' pass; +podioevent = EICDataSvc("EventDataSvc", inputs=[sim_file]) + +# S13660-3050AE-08 SiPM quantum efficiency [(wavelength [nm], q.e.)] +# Note: this is consistent with S13361-3050AE-08 (for DRICH) +qe_data = [ + (325, 0.04), + (340, 0.10), + (350, 0.20), + (370, 0.30), + (400, 0.35), + (450, 0.40), + (500, 0.38), + (550, 0.35), + (600, 0.27), + (650, 0.20), + (700, 0.15), + (750, 0.12), + (800, 0.08), + (850, 0.06), + (900, 0.04) +] + +radiators = [ + "Aerogel zbins=5 smearing=gaussian 2mrad rindex=1.0190 attenuation[mm]=48.0", + "GasVolume zbins=10 smearing=gaussian 5mrad rindex=1.00076" +] + +podioinput = PodioInput( + "PodioReader", + # Input collections: MC truth tracks and DRICH raw hits (photons); + collections=["MCParticles", "DRICHHits"], + OutputLevel=DEBUG + ) + +irtrec = IRTAlgorithm( + # Input collections: MC truth tracks and DRICH raw hits (photons); + inputMCParticles="MCParticles", + + # DRICH optics configuration produced by DRICH_geo.cpp code along with the dd4hep XML file; + ConfigFile=irt_auxfile, + Radiators=[ (r) for r in radiators ], + + # SiPM PDE; FIXME: units.eV coefficient gives extra x1000 (?); + QEcurve=[ ((1239.84/a), b) for a, b in qe_data ], + # Rebin the QE in that many equidistant bins internally; + QEbins="100", + # SiPM geometric fill factor and "safety factor" for the photon count estimates; + #GeometricEfficiency="1.00", + SafetyFactor="0.70", + ) + +# Output ROOT file; keep the input collections as well, append DRICH PID tables; +out = PodioOutput("out", filename=out_file) +out.outputCommands = ["keep *"] + +ApplicationMgr( + TopAlg = [podioinput, irtrec, out], + EvtSel = 'NONE', + # Process that many events; + EvtMax = 50000, + ExtSvc = [podioevent], + OutputLevel = DEBUG, + PluginDebugLevel = 2 + ) + -- GitLab From 542884eede49fdbf609a44db5e4744f687a3a53a Mon Sep 17 00:00:00 2001 From: Christopher Dilks Date: Mon, 15 Aug 2022 20:03:05 -0400 Subject: [PATCH 6/6] modified: JugPID/src/components/IRTAlgorithm.cpp --- JugPID/src/components/IRTAlgorithm.cpp | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/JugPID/src/components/IRTAlgorithm.cpp b/JugPID/src/components/IRTAlgorithm.cpp index 33048b1b..9065a2f4 100644 --- a/JugPID/src/components/IRTAlgorithm.cpp +++ b/JugPID/src/components/IRTAlgorithm.cpp @@ -270,7 +270,6 @@ StatusCode Jug::PID::IRTAlgorithm::initialize( void ) StatusCode Jug::PID::IRTAlgorithm::execute( void ) { - printf("####################################### new event\n"); // Input collection(s); const auto &hits = *m_inputHitCollection->get(); const auto &mctracks = *m_inputMCParticles.get(); @@ -336,10 +335,8 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) event->AddChargedParticle(particle); // Start radiator history for all known radiators; - for(auto rptr: m_IrtDet->Radiators()) { + for(auto rptr: m_IrtDet->Radiators()) particle->StartRadiatorHistory(std::make_pair(rptr.second, new RadiatorHistory())); - // printf("found radiator %s\n",rptr.second->GetAlternativeMaterialName()); - } // Distribute photons over radiators where they were presumably produced, to // create a structure which would mimic a standalone G4 stepper behavior; @@ -348,7 +345,6 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) for(const auto &hit: hits) { // FIXME: range checks; const auto &phtrack = hit.getMCParticle(); - // printf("here0\n"); if(phtrack.getPDG()!=-22) { fprintf(stderr,"ERROR: this hit was not from an optical photon: PDG = %d\n",phtrack.getPDG()); continue; @@ -356,6 +352,7 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // FIXME: yes, use MC truth here; not really needed I guess; // FIXME: range checks; + //printf("Here: %d %d %3d!\n", phtrack.parents()[0], mctrack.ID(), hit.truth().trackID); //if (phtrack.parents()[0] != mctrack.ID()) continue; // Vertex where photon was created; @@ -374,17 +371,16 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // FIXME: this is in general correct, but needs refinement; TVector3 ip(0,0,0); // FIXME: unify with the code below; - // printf("guessing radiator..\n"); auto radiator = m_IrtDet->GuessRadiator(vtx, (vtx - ip).Unit()); // FIXME: do it better later; if (!radiator) continue; - // printf("Here-1\n"); + //printf("Here-1\n"); // Simulate QE & geometric sensor efficiency; FIXME: hit.energy() is numerically // in GeV units, but Gaudi::Units::GeV = 1000; prefer to convert photon energies // to [eV] in all places by hand; double eVenergy = 1E9*hit.getEDep(); - // printf("%f\n", eVenergy); + //printf("%f\n", eVenergy); if (!QE_pass(eVenergy, m_rngUni()) || m_rngUni() > /*m_GeometricEfficiency.value()**/m_SafetyFactor.value()) { @@ -397,10 +393,10 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) // needed if ACTS trajectory parameterization is used; particle->FindRadiatorHistory(radiator)->AddStepBufferPoint(phtrack.getTime(), vtx); - // printf("Here-2A\n"); + //printf("Here-2A\n"); continue; } //if - // printf("Here-2B\n"); + //printf("Here-2B\n"); //if (vs.z < 2500 || vs.z > 2700) continue; //if (vs.z < 2200 || vs.z > 2400) continue; -- GitLab