diff --git a/JugPID/src/components/IRTAlgorithm.cpp b/JugPID/src/components/IRTAlgorithm.cpp index b606d3607c089718682e11d73a7f95e4fee89f8d..9065a2f4a8af175763b999e0e2bd3ea259dcfc01 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,16 +326,16 @@ 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; - for(auto rptr: m_IrtDet->Radiators()) + for(auto rptr: m_IrtDet->Radiators()) particle->StartRadiatorHistory(std::make_pair(rptr.second, new RadiatorHistory())); // Distribute photons over radiators where they were presumably produced, to @@ -344,8 +344,11 @@ 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(); + 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; @@ -353,14 +356,15 @@ StatusCode Jug::PID::IRTAlgorithm::execute( void ) //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); + // 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); - //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; @@ -375,7 +379,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 +391,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 +406,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 +442,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 +490,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 +590,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 +598,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 +645,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); #endif } } //for rctrack (mctrack) diff --git a/JugPID/src/components/IRTAlgorithm.h b/JugPID/src/components/IRTAlgorithm.h index d6300caa5a3a494d92cd38bba33f163471b6971f..e38d187514075d0d5058e18b35ce3d204dcd2136 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 diff --git a/JugPID/tests/options/irt.py b/JugPID/tests/options/irt.py new file mode 100644 index 0000000000000000000000000000000000000000..657060d616fe70dab236fd75f6c3153ba03219f8 --- /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 + ) +