Pandora
Pandora source code navigator
Loading...
Searching...
No Matches
LArPandoraInput.cxx
Go to the documentation of this file.
1
8
9#include "larcore/Geometry/Geometry.h"
10#include "larcorealg/Geometry/PlaneGeo.h"
11#include "larcorealg/Geometry/TPCGeo.h"
12#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
13#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h"
14
15#include "lardataobj/RecoBase/Hit.h"
16
17#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
18#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
19
20#include "nusimdata/SimulationBase/MCTruth.h"
21
22#include "larsim/MCCheater/ParticleInventoryService.h"
23
24#include "lardata/DetectorInfoServices/DetectorClocksService.h"
25#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
26#include "lardata/DetectorInfoServices/LArPropertiesService.h"
27
28#include "Api/PandoraApi.h"
31
33
37
38#include "messagefacility/MessageLogger/MessageLogger.h"
39
40#include <limits>
41#include <utility>
42
43namespace lar_pandora {
44
45 void LArPandoraInput::CreatePandoraHits2D(const art::Event& e,
46 const Settings& settings,
47 const LArDriftVolumeMap& driftVolumeMap,
48 const HitVector& hitVector,
49 IdToHitMap& idToHitMap)
50 {
51 mf::LogDebug("LArPandora") << " *** LArPandoraInput::CreatePandoraHits2D(...) *** "
52 << std::endl;
53
54 if (!settings.m_pPrimaryPandora)
55 throw cet::exception("LArPandora")
56 << "CreatePandoraHits2D - primary Pandora instance does not exist ";
57
58 const pandora::Pandora* pPandora(settings.m_pPrimaryPandora);
59
60 art::ServiceHandle<geo::Geometry const> theGeometry;
61 auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e);
63
64 // Loop over ART hits
65 int hitCounter(settings.m_hitCounterOffset);
66
67 lar_content::LArCaloHitFactory caloHitFactory;
68
69 for (HitVector::const_iterator iter = hitVector.begin(), iterEnd = hitVector.end();
70 iter != iterEnd;
71 ++iter) {
72 const art::Ptr<recob::Hit> hit = *iter;
73 const geo::WireID hit_WireID(hit->WireID());
74
75 // Get basic hit properties (view, time, charge)
76 const geo::View_t hit_View(hit->View());
77 const double hit_Charge(hit->Integral());
78 const double hit_Time(hit->PeakTime());
79 const double hit_TimeStart(hit->PeakTimeMinusRMS());
80 const double hit_TimeEnd(hit->PeakTimePlusRMS());
81
82 // Get hit X coordinate and, if using a single global drift volume, remove any out-of-time hits here
83 const double xpos_cm(
84 detProp.ConvertTicksToX(hit_Time, hit_WireID.Plane, hit_WireID.TPC, hit_WireID.Cryostat));
85 const double dxpos_cm(
86 std::fabs(detProp.ConvertTicksToX(
87 hit_TimeEnd, hit_WireID.Plane, hit_WireID.TPC, hit_WireID.Cryostat) -
88 detProp.ConvertTicksToX(
89 hit_TimeStart, hit_WireID.Plane, hit_WireID.TPC, hit_WireID.Cryostat)));
90
91 // Get hit Y and Z coordinates, based on central position of wire
92 auto const xyz = theGeometry->Wire(hit_WireID).GetCenter();
93 const double y0_cm(xyz.Y());
94 const double z0_cm(xyz.Z());
95
96 // Get other hit properties here
97 const double wire_pitch_cm(theGeometry->WirePitch(hit_View)); // cm
98 const double mips(LArPandoraInput::GetMips(detProp, settings, hit_Charge, hit_View));
99
100 // Create Pandora CaloHit
101 lar_content::LArCaloHitParameters caloHitParameters;
102
103 try {
104 caloHitParameters.m_expectedDirection = pandora::CartesianVector(0., 0., 1.);
105 caloHitParameters.m_cellNormalVector = pandora::CartesianVector(0., 0., 1.);
106 caloHitParameters.m_cellSize0 = settings.m_dx_cm;
107 caloHitParameters.m_cellSize1 = (settings.m_useHitWidths ? dxpos_cm : settings.m_dx_cm);
108 caloHitParameters.m_cellThickness = wire_pitch_cm;
109 caloHitParameters.m_cellGeometry = pandora::RECTANGULAR;
110 caloHitParameters.m_time = 0.;
111 caloHitParameters.m_nCellRadiationLengths = settings.m_dx_cm / settings.m_rad_cm;
112 caloHitParameters.m_nCellInteractionLengths = settings.m_dx_cm / settings.m_int_cm;
113 caloHitParameters.m_isDigital = false;
114 caloHitParameters.m_hitRegion = pandora::SINGLE_REGION;
115 caloHitParameters.m_layer = 0;
116 caloHitParameters.m_isInOuterSamplingLayer = false;
117 caloHitParameters.m_inputEnergy = hit_Charge;
118 caloHitParameters.m_mipEquivalentEnergy = mips;
119 caloHitParameters.m_electromagneticEnergy = mips * settings.m_mips_to_gev;
120 caloHitParameters.m_hadronicEnergy = mips * settings.m_mips_to_gev;
121 caloHitParameters.m_pParentAddress = (void*)((intptr_t)(++hitCounter));
122 caloHitParameters.m_larTPCVolumeId =
123 LArPandoraGeometry::GetVolumeID(driftVolumeMap, hit_WireID.Cryostat, hit_WireID.TPC);
125 driftVolumeMap, hit_WireID.Cryostat, hit_WireID.TPC);
126
127 if (hit_View == detType->TargetViewW(hit_WireID.TPC, hit_WireID.Cryostat)) {
128 caloHitParameters.m_hitType = pandora::TPC_VIEW_W;
129 const double wpos_cm(
130 pPandora->GetPlugins()->GetLArTransformationPlugin()->YZtoW(y0_cm, z0_cm));
131 caloHitParameters.m_positionVector = pandora::CartesianVector(xpos_cm, 0., wpos_cm);
132 }
133 else if (hit_View == detType->TargetViewU(hit_WireID.TPC, hit_WireID.Cryostat)) {
134 caloHitParameters.m_hitType = pandora::TPC_VIEW_U;
135 const double upos_cm(
136 pPandora->GetPlugins()->GetLArTransformationPlugin()->YZtoU(y0_cm, z0_cm));
137 caloHitParameters.m_positionVector = pandora::CartesianVector(xpos_cm, 0., upos_cm);
138 }
139 else if (hit_View == detType->TargetViewV(hit_WireID.TPC, hit_WireID.Cryostat)) {
140 caloHitParameters.m_hitType = pandora::TPC_VIEW_V;
141 const double vpos_cm(
142 pPandora->GetPlugins()->GetLArTransformationPlugin()->YZtoV(y0_cm, z0_cm));
143 caloHitParameters.m_positionVector = pandora::CartesianVector(xpos_cm, 0., vpos_cm);
144 }
145 else {
146 throw cet::exception("LArPandora")
147 << "CreatePandoraHits2D - this wire view not recognised (View=" << hit_View << ") ";
148 }
149 }
150 catch (const pandora::StatusCodeException&) {
151 mf::LogWarning("LArPandora")
152 << "CreatePandoraHits2D - invalid calo hit parameter provided, all assigned values must "
153 "be finite, calo hit omitted "
154 << std::endl;
155 continue;
156 }
157
158 // Store the hit address
159 if (hitCounter >= settings.m_uidOffset)
160 throw cet::exception("LArPandora")
161 << "CreatePandoraHits2D - detected an excessive number of hits (" << hitCounter << ") ";
162
163 idToHitMap[hitCounter] = hit;
164
165 // Create the Pandora hit
166 try {
168 pandora::STATUS_CODE_SUCCESS,
169 !=,
170 PandoraApi::CaloHit::Create(*pPandora, caloHitParameters, caloHitFactory));
171 }
172 catch (const pandora::StatusCodeException&) {
173 mf::LogWarning("LArPandora") << "CreatePandoraHits2D - unable to create calo hit, "
174 "insufficient or invalid information supplied "
175 << std::endl;
176 continue;
177 }
178 }
179 }
180
181 //------------------------------------------------------------------------------------------------------------------------------------------
182
184 const LArDriftVolumeList& driftVolumeList)
185 {
186 mf::LogDebug("LArPandora") << " *** LArPandoraInput::CreatePandoraLArTPCs(...) *** "
187 << std::endl;
188
189 if (!settings.m_pPrimaryPandora)
190 throw cet::exception("LArPandora")
191 << "CreatePandoraDetectorGaps - primary Pandora instance does not exist ";
192
193 const pandora::Pandora* pPandora(settings.m_pPrimaryPandora);
194
195 for (const LArDriftVolume& driftVolume : driftVolumeList) {
197
198 try {
199 parameters.m_larTPCVolumeId = driftVolume.GetVolumeID();
200 parameters.m_centerX = driftVolume.GetCenterX();
201 parameters.m_centerY = driftVolume.GetCenterY();
202 parameters.m_centerZ = driftVolume.GetCenterZ();
203 parameters.m_widthX = driftVolume.GetWidthX();
204 parameters.m_widthY = driftVolume.GetWidthY();
205 parameters.m_widthZ = driftVolume.GetWidthZ();
206 parameters.m_wirePitchU = driftVolume.GetWirePitchU();
207 parameters.m_wirePitchV = driftVolume.GetWirePitchV();
208 parameters.m_wirePitchW = driftVolume.GetWirePitchW();
209 parameters.m_wireAngleU = driftVolume.GetWireAngleU();
210 parameters.m_wireAngleV = driftVolume.GetWireAngleV();
211 parameters.m_wireAngleW = driftVolume.GetWireAngleW();
212 parameters.m_sigmaUVW = driftVolume.GetSigmaUVZ();
213 parameters.m_isDriftInPositiveX = driftVolume.IsPositiveDrift();
214 }
215 catch (const pandora::StatusCodeException&) {
216 mf::LogWarning("LArPandora") << "CreatePandoraLArTPCs - invalid tpc parameter provided, "
217 "all assigned values must be finite, tpc omitted "
218 << std::endl;
219 continue;
220 }
221
222 try {
223 PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS,
224 !=,
225 PandoraApi::Geometry::LArTPC::Create(*pPandora, parameters));
226 }
227 catch (const pandora::StatusCodeException&) {
228 mf::LogWarning("LArPandora") << "CreatePandoraLArTPCs - unable to create tpc, insufficient "
229 "or invalid information supplied "
230 << std::endl;
231 continue;
232 }
233 }
234 }
235
236 //------------------------------------------------------------------------------------------------------------------------------------------
237
239 const LArDriftVolumeList& driftVolumeList,
240 const LArDetectorGapList& listOfGaps)
241 {
242 //ATTN - Unlike SP, DP detector gaps are not in the drift direction
243 art::ServiceHandle<geo::Geometry const> theGeometry;
245
246 mf::LogDebug("LArPandora") << " *** LArPandoraInput::CreatePandoraDetectorGaps(...) *** "
247 << std::endl;
248
249 if (!settings.m_pPrimaryPandora)
250 throw cet::exception("LArPandora")
251 << "CreatePandoraDetectorGaps - primary Pandora instance does not exist ";
252
253 const pandora::Pandora* pPandora(settings.m_pPrimaryPandora);
254
255 for (const LArDetectorGap& gap : listOfGaps) {
257
258 try {
259 parameters = detType->CreateLineGapParametersFromDetectorGaps(gap);
260 }
261 catch (const pandora::StatusCodeException&) {
262 mf::LogWarning("LArPandora")
263 << "CreatePandoraDetectorGaps - invalid line gap parameter provided, all assigned values "
264 "must be finite, line gap omitted "
265 << std::endl;
266 continue;
267 }
268 try {
269 PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS,
270 !=,
271 PandoraApi::Geometry::LineGap::Create(*pPandora, parameters));
272 }
273 catch (const pandora::StatusCodeException&) {
274 mf::LogWarning("LArPandora") << "CreatePandoraDetectorGaps - unable to create line gap, "
275 "insufficient or invalid information supplied "
276 << std::endl;
277 continue;
278 }
279 }
280 }
281
282 //------------------------------------------------------------------------------------------------------------------------------------------
283
285 const LArDriftVolumeMap& driftVolumeMap)
286 {
287 mf::LogDebug("LArPandora") << " *** LArPandoraInput::CreatePandoraReadoutGaps(...) *** "
288 << std::endl;
289
290 if (!settings.m_pPrimaryPandora)
291 throw cet::exception("LArPandora")
292 << "CreatePandoraReadoutGaps - primary Pandora instance does not exist ";
293
294 const pandora::Pandora* pPandora(settings.m_pPrimaryPandora);
295
296 art::ServiceHandle<geo::Geometry const> theGeometry;
297 const lariov::ChannelStatusProvider& channelStatus(
298 art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider());
299
301
302 for (auto const& plane : theGeometry->Iterate<geo::PlaneGeo>()) {
303 const float halfWirePitch(0.5f * theGeometry->WirePitch(plane.View()));
304 const unsigned int nWires(theGeometry->Nwires(plane.ID()));
305
306 int firstBadWire(-1), lastBadWire(-1);
307
308 for (unsigned int iwire = 0; iwire < nWires; ++iwire) {
309 const raw::ChannelID_t channel(
310 theGeometry->PlaneWireToChannel(geo::WireID{plane.ID(), iwire}));
311 const bool isBadChannel(channelStatus.IsBad(channel));
312 const bool isLastWire(nWires == (iwire + 1));
313
314 if (isBadChannel && (firstBadWire < 0)) firstBadWire = iwire;
315
316 if (isBadChannel || isLastWire) lastBadWire = iwire;
317
318 if (isBadChannel && !isLastWire) continue;
319
320 if ((firstBadWire < 0) || (lastBadWire < 0)) continue;
321
322 auto const firstXYZ = plane.Wire(firstBadWire).GetCenter();
323 auto const lastXYZ = plane.Wire(lastBadWire).GetCenter();
324
325 firstBadWire = -1;
326 lastBadWire = -1;
327
329
330 try {
331 float xFirst(-std::numeric_limits<float>::max());
332 float xLast(std::numeric_limits<float>::max());
333
334 auto const [icstat, itpc] = std::make_pair(plane.ID().Cryostat, plane.ID().TPC);
335 const unsigned int volumeId(
336 LArPandoraGeometry::GetVolumeID(driftVolumeMap, icstat, itpc));
337 LArDriftVolumeMap::const_iterator volumeIter(driftVolumeMap.find(volumeId));
338
339 if (driftVolumeMap.end() != volumeIter) {
340 xFirst = volumeIter->second.GetCenterX() - 0.5f * volumeIter->second.GetWidthX();
341 xLast = volumeIter->second.GetCenterX() + 0.5f * volumeIter->second.GetWidthX();
342 }
343
344 const geo::View_t iview = plane.View();
345 parameters = detType->CreateLineGapParametersFromReadoutGaps(
346 iview, itpc, icstat, firstXYZ, lastXYZ, halfWirePitch, xFirst, xLast, pPandora);
347 }
348 catch (const pandora::StatusCodeException&) {
349 mf::LogWarning("LArPandora")
350 << "CreatePandoraReadoutGaps - invalid line gap parameter provided, all assigned "
351 "values must be finite, line gap omitted "
352 << std::endl;
353 continue;
354 }
355
356 try {
357 PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS,
358 !=,
359 PandoraApi::Geometry::LineGap::Create(*pPandora, parameters));
360 }
361 catch (const pandora::StatusCodeException&) {
362 mf::LogWarning("LArPandora") << "CreatePandoraReadoutGaps - unable to create line "
363 "gap, insufficient or invalid information supplied "
364 << std::endl;
365 continue;
366 }
367 }
368 }
369 }
370
371 //------------------------------------------------------------------------------------------------------------------------------------------
372
374 const Settings& settings,
375 const MCTruthToMCParticles& truthToParticleMap,
376 const MCParticlesToMCTruth& particleToTruthMap,
377 const RawMCParticleVector& generatorMCParticleVector)
378 {
379 mf::LogDebug("LArPandora") << " *** LArPandoraInput::CreatePandoraMCParticles(...) *** "
380 << std::endl;
381 art::ServiceHandle<cheat::ParticleInventoryService const> particleInventoryService;
382
383 if (!settings.m_pPrimaryPandora)
384 throw cet::exception("LArPandora")
385 << "CreatePandoraMCParticles - primary Pandora instance does not exist ";
386
387 const pandora::Pandora* pPandora(settings.m_pPrimaryPandora);
388
389 // Make indexed list of MC particles
390 MCParticleMap particleMap;
391
392 for (MCParticlesToMCTruth::const_iterator iter = particleToTruthMap.begin(),
393 iterEnd = particleToTruthMap.end();
394 iter != iterEnd;
395 ++iter) {
396 const art::Ptr<simb::MCParticle> particle = iter->first;
397 particleMap[particle->TrackId()] = particle;
398 }
399
400 // Loop over MC truth objects
401 int neutrinoCounter(0);
402
403 lar_content::LArMCParticleFactory mcParticleFactory;
404
405 for (MCTruthToMCParticles::const_iterator iter1 = truthToParticleMap.begin(),
406 iterEnd1 = truthToParticleMap.end();
407 iter1 != iterEnd1;
408 ++iter1) {
409 const art::Ptr<simb::MCTruth> truth = iter1->first;
410
411 if (truth->NeutrinoSet()) {
412 const simb::MCNeutrino neutrino(truth->GetNeutrino());
413 ++neutrinoCounter;
414
415 if (neutrinoCounter >= settings.m_uidOffset)
416 throw cet::exception("LArPandora")
417 << "CreatePandoraMCParticles - detected an excessive number of mc neutrinos ("
418 << neutrinoCounter << ")";
419
420 const int neutrinoID(neutrinoCounter + 9 * settings.m_uidOffset);
421
422 // Create Pandora 3D MC Particle
423 lar_content::LArMCParticleParameters mcParticleParameters;
424
425 try {
426 mcParticleParameters.m_nuanceCode = neutrino.InteractionType();
427 mcParticleParameters.m_process = lar_content::MC_PROC_INCIDENT_NU;
428 mcParticleParameters.m_energy = neutrino.Nu().E();
429 mcParticleParameters.m_momentum =
430 pandora::CartesianVector(neutrino.Nu().Px(), neutrino.Nu().Py(), neutrino.Nu().Pz());
431 mcParticleParameters.m_vertex =
432 pandora::CartesianVector(neutrino.Nu().Vx(), neutrino.Nu().Vy(), neutrino.Nu().Vz());
433 mcParticleParameters.m_endpoint =
434 pandora::CartesianVector(neutrino.Nu().Vx(), neutrino.Nu().Vy(), neutrino.Nu().Vz());
435 mcParticleParameters.m_particleId = neutrino.Nu().PdgCode();
436 mcParticleParameters.m_mcParticleType = pandora::MC_3D;
437 mcParticleParameters.m_pParentAddress = (void*)((intptr_t)neutrinoID);
438 }
439 catch (const pandora::StatusCodeException&) {
440 mf::LogWarning("LArPandora")
441 << "CreatePandoraMCParticles - invalid mc neutrino parameter provided, all assigned "
442 "values must be finite, mc neutrino omitted "
443 << std::endl;
444 continue;
445 }
446
447 try {
449 pandora::STATUS_CODE_SUCCESS,
450 !=,
451 PandoraApi::MCParticle::Create(*pPandora, mcParticleParameters, mcParticleFactory));
452 }
453 catch (const pandora::StatusCodeException&) {
454 mf::LogWarning("LArPandora") << "CreatePandoraMCParticles - unable to create mc "
455 "neutrino, insufficient or invalid information supplied "
456 << std::endl;
457 continue;
458 }
459
460 // Loop over associated particles
461 const MCParticleVector& particleVector = iter1->second;
462
463 for (MCParticleVector::const_iterator iter2 = particleVector.begin(),
464 iterEnd2 = particleVector.end();
465 iter2 != iterEnd2;
466 ++iter2) {
467 const art::Ptr<simb::MCParticle> particle = *iter2;
468 const int trackID(particle->TrackId());
469
470 // Mother/Daughter Links
471 if (particle->Mother() == 0) {
472 try {
474 pandora::STATUS_CODE_SUCCESS,
475 !=,
477 *pPandora, (void*)((intptr_t)neutrinoID), (void*)((intptr_t)trackID)));
478 }
479 catch (const pandora::StatusCodeException&) {
480 mf::LogWarning("LArPandora") << "CreatePandoraMCParticles - unable to create mc "
481 "particle relationship, invalid information supplied "
482 << std::endl;
483 continue;
484 }
485 }
486 }
487 }
488 }
489
490 mf::LogDebug("LArPandora") << " Number of Pandora neutrinos: " << neutrinoCounter
491 << std::endl;
492
493 // Loop over G4 particles
494 int particleCounter(0);
495
496 // Find Primary Generator Particles
497 std::map<const simb::MCParticle, bool> primaryGeneratorMCParticleMap;
498 LArPandoraInput::FindPrimaryParticles(generatorMCParticleVector, primaryGeneratorMCParticleMap);
499
500 for (MCParticleMap::const_iterator iterI = particleMap.begin(), iterEndI = particleMap.end();
501 iterI != iterEndI;
502 ++iterI) {
503 const art::Ptr<simb::MCParticle> particle = iterI->second;
504
505 if (particle->TrackId() != iterI->first)
506 throw cet::exception("LArPandora") << "CreatePandoraMCParticles - mc truth information "
507 "appears to be scrambled in this event";
508
509 if (particle->TrackId() >= settings.m_uidOffset)
510 throw cet::exception("LArPandora")
511 << "CreatePandoraMCParticles - detected an excessive number of MC particles ("
512 << particle->TrackId() << ")";
513
514 ++particleCounter;
515
516 // Find start and end trajectory points
517 int firstT(-1), lastT(-1);
518 LArPandoraInput::GetTrueStartAndEndPoints(settings, particle, firstT, lastT);
519
520 if (firstT < 0 && lastT < 0) {
521 firstT = 0;
522 lastT = 0;
523 }
524
525 // Lookup position and kinematics at start and end points
526 const float vtxX(particle->Vx(firstT));
527 const float vtxY(particle->Vy(firstT));
528 const float vtxZ(particle->Vz(firstT));
529
530 const float endX(particle->Vx(lastT));
531 const float endY(particle->Vy(lastT));
532 const float endZ(particle->Vz(lastT));
533
534 const float pX(particle->Px(firstT));
535 const float pY(particle->Py(firstT));
536 const float pZ(particle->Pz(firstT));
537 const float E(particle->E(firstT));
538
539 // Find the source of the mc particle
540 int nuanceCode(0);
541 const int trackID(particle->TrackId());
542 const simb::Origin_t origin(particleInventoryService->TrackIdToMCTruth(trackID).Origin());
543
544 if (LArPandoraInput::IsPrimaryMCParticle(particle, primaryGeneratorMCParticleMap)) {
545 nuanceCode = 2001;
546 }
547 else if (simb::kCosmicRay == origin) {
548 nuanceCode = 3000;
549 }
550 else if (simb::kSingleParticle == origin) {
551 nuanceCode = 2000;
552 }
553
554 // Create 3D Pandora MC Particle
555 lar_content::LArMCParticleParameters mcParticleParameters;
556
557 try {
558 MCProcessMap processMap;
559 FillMCProcessMap(processMap);
560 mcParticleParameters.m_nuanceCode = nuanceCode;
561 if (processMap.find(particle->Process()) != processMap.end()) {
562 mcParticleParameters.m_process = processMap[particle->Process()];
563 }
564 else {
565 mcParticleParameters.m_process = lar_content::MC_PROC_UNKNOWN;
566 mf::LogWarning("LArPandora")
567 << "CreatePandoraMCParticles - found an unknown process" << std::endl;
568 }
569 mcParticleParameters.m_energy = E;
570 mcParticleParameters.m_particleId = particle->PdgCode();
571 mcParticleParameters.m_momentum = pandora::CartesianVector(pX, pY, pZ);
572 mcParticleParameters.m_vertex = pandora::CartesianVector(vtxX, vtxY, vtxZ);
573 mcParticleParameters.m_endpoint = pandora::CartesianVector(endX, endY, endZ);
574 mcParticleParameters.m_mcParticleType = pandora::MC_3D;
575 mcParticleParameters.m_pParentAddress = (void*)((intptr_t)particle->TrackId());
576 }
577 catch (const pandora::StatusCodeException&) {
578 mf::LogWarning("LArPandora")
579 << "CreatePandoraMCParticles - invalid mc particle parameter provided, all assigned "
580 "values must be finite, mc particle omitted "
581 << std::endl;
582 continue;
583 }
584
585 try {
587 pandora::STATUS_CODE_SUCCESS,
588 !=,
589 PandoraApi::MCParticle::Create(*pPandora, mcParticleParameters, mcParticleFactory));
590 }
591 catch (const pandora::StatusCodeException&) {
592 mf::LogWarning("LArPandora") << "CreatePandoraMCParticles - unable to create mc particle, "
593 "insufficient or invalid information supplied "
594 << std::endl;
595 continue;
596 }
597
598 // Create Mother/Daughter Links between 3D MC Particles
599 const int id_mother(particle->Mother());
600 MCParticleMap::const_iterator iterJ = particleMap.find(id_mother);
601
602 if (iterJ != particleMap.end()) {
603 try {
605 pandora::STATUS_CODE_SUCCESS,
606 !=,
608 *pPandora, (void*)((intptr_t)id_mother), (void*)((intptr_t)particle->TrackId())));
609 }
610 catch (const pandora::StatusCodeException&) {
611 mf::LogWarning("LArPandora") << "CreatePandoraMCParticles - Unable to create mc particle "
612 "relationship, invalid information supplied "
613 << std::endl;
614 continue;
615 }
616 }
617 }
618
619 mf::LogDebug("LArPandora") << "Number of mc particles: " << particleCounter << std::endl;
620 }
621
622 //------------------------------------------------------------------------------------------------------------------------------------------
623
625 const RawMCParticleVector& mcParticleVector,
626 std::map<const simb::MCParticle, bool>& primaryMCParticleMap)
627 {
628 for (const simb::MCParticle& mcParticle : mcParticleVector) {
629 if ("primary" == mcParticle.Process()) {
630 primaryMCParticleMap.emplace(std::make_pair(mcParticle, false));
631 }
632 }
633 }
634
635 //------------------------------------------------------------------------------------------------------------------------------------------
636
638 const art::Ptr<simb::MCParticle>& mcParticle,
639 std::map<const simb::MCParticle, bool>& primaryMCParticleMap)
640 {
641 for (auto& mcParticleIter : primaryMCParticleMap) {
642 if (!mcParticleIter.second) {
643 const simb::MCParticle primaryMCParticle(mcParticleIter.first);
644
645 if (std::fabs(primaryMCParticle.Px() - mcParticle->Px()) <
646 std::numeric_limits<double>::epsilon() &&
647 std::fabs(primaryMCParticle.Py() - mcParticle->Py()) <
648 std::numeric_limits<double>::epsilon() &&
649 std::fabs(primaryMCParticle.Pz() - mcParticle->Pz()) <
650 std::numeric_limits<double>::epsilon()) {
651 mcParticleIter.second = true;
652 return true;
653 }
654 }
655 }
656 return false;
657 }
658
659 //------------------------------------------------------------------------------------------------------------------------------------------
660
662 const IdToHitMap& idToHitMap,
663 const HitsToTrackIDEs& hitToParticleMap)
664 {
665 mf::LogDebug("LArPandora") << " *** LArPandoraInput::CreatePandoraMCLinks(...) *** "
666 << std::endl;
667
668 if (!settings.m_pPrimaryPandora)
669 throw cet::exception("LArPandora")
670 << "CreatePandoraMCLinks2D - primary Pandora instance does not exist ";
671
672 const pandora::Pandora* pPandora(settings.m_pPrimaryPandora);
673
674 for (IdToHitMap::const_iterator iterI = idToHitMap.begin(), iterEndI = idToHitMap.end();
675 iterI != iterEndI;
676 ++iterI) {
677 const int hitID(iterI->first);
678 const art::Ptr<recob::Hit> hit(iterI->second);
679 // const geo::WireID hit_WireID(hit->WireID());
680
681 // Get list of associated MC particles
682 HitsToTrackIDEs::const_iterator iterJ = hitToParticleMap.find(hit);
683
684 if (hitToParticleMap.end() == iterJ) continue;
685
686 const TrackIDEVector& trackCollection = iterJ->second;
687
688 if (trackCollection.size() == 0)
689 throw cet::exception("LArPandora")
690 << "CreatePandoraMCLinks2D - found a hit without any associated MC truth information ";
691
692 // Create links between hits and MC particles
693 for (unsigned int k = 0; k < trackCollection.size(); ++k) {
694 const sim::TrackIDE trackIDE(trackCollection.at(k));
695 const int trackID(std::abs(trackIDE.trackID)); // TODO: Find out why std::abs is needed
696 const float energyFrac(trackIDE.energyFrac);
697
698 try {
700 pandora::STATUS_CODE_SUCCESS,
701 !=,
703 *pPandora, (void*)((intptr_t)hitID), (void*)((intptr_t)trackID), energyFrac));
704 }
705 catch (const pandora::StatusCodeException&) {
706 mf::LogWarning("LArPandora") << "CreatePandoraMCLinks2D - unable to create calo hit to "
707 "mc particle relationship, invalid information supplied "
708 << std::endl;
709 continue;
710 }
711 }
712 }
713 }
714
715 //------------------------------------------------------------------------------------------------------------------------------------------
716
718 const art::Ptr<simb::MCParticle>& particle,
719 int& firstT,
720 int& lastT)
721 {
722 art::ServiceHandle<geo::Geometry const> theGeometry;
723 firstT = -1;
724 lastT = -1;
725
726 for (auto const& tpcid : theGeometry->Iterate<geo::TPCID>()) {
727 int thisfirstT(-1), thislastT(-1);
728 LArPandoraInput::GetTrueStartAndEndPoints(tpcid, particle, thisfirstT, thislastT);
729
730 if (thisfirstT < 0) continue;
731
732 if (firstT < 0 || thisfirstT < firstT) firstT = thisfirstT;
733
734 if (lastT < 0 || thislastT > lastT) lastT = thislastT;
735 }
736 }
737
738 //------------------------------------------------------------------------------------------------------------------------------------------
739
740 void LArPandoraInput::GetTrueStartAndEndPoints(const geo::TPCID& ref_tpcid,
741 const art::Ptr<simb::MCParticle>& particle,
742 int& startT,
743 int& endT)
744 {
745 art::ServiceHandle<geo::Geometry const> theGeometry;
746
747 bool foundStartPosition(false);
748 const int numTrajectoryPoints(static_cast<int>(particle->NumberTrajectoryPoints()));
749
750 for (int nt = 0; nt < numTrajectoryPoints; ++nt) {
751 const geo::Point_t pos{particle->Vx(nt), particle->Vy(nt), particle->Vz(nt)};
752 geo::TPCID tpcID = theGeometry->FindTPCAtPosition(pos);
753
754 if (!tpcID.isValid) continue;
755
756 if (tpcID != ref_tpcid) continue;
757
758 endT = nt;
759
760 if (!foundStartPosition) {
761 startT = endT;
762 foundStartPosition = true;
763 }
764 }
765 }
766
767 //------------------------------------------------------------------------------------------------------------------------------------------
768
769 float LArPandoraInput::GetTrueX0(const art::Event& e,
770 const art::Ptr<simb::MCParticle>& particle,
771 const int nt)
772 {
773 art::ServiceHandle<geo::Geometry const> theGeometry;
774 auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(e);
775 auto const det_prop =
776 art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e, clock_data);
777
778 geo::Point_t const pos{particle->Vx(nt), particle->Vy(nt), particle->Vz(nt)};
779 auto const tpcID = theGeometry->PositionToTPCID(pos);
780
781 const float vtxT(particle->T(nt));
782 const float vtxTDC(clock_data.TPCG4Time2Tick(vtxT));
783 const float vtxTDC0(trigger_offset(clock_data));
784
785 const geo::TPCGeo& theTpc = theGeometry->TPC(tpcID);
786 const float driftDir((theTpc.DriftDirection() == geo::kNegX) ? +1.0 : -1.0);
787 return (driftDir * (vtxTDC - vtxTDC0) * det_prop.GetXTicksCoefficient());
788 }
789
790 //------------------------------------------------------------------------------------------------------------------------------------------
791
792 double LArPandoraInput::GetMips(detinfo::DetectorPropertiesData const& detProp,
793 const Settings& settings,
794 const double hit_Charge,
795 const geo::View_t hit_View)
796 {
797 art::ServiceHandle<geo::Geometry const> theGeometry;
798
799 // TODO: Unite this procedure with other calorimetry procedures under development
800 const double dQdX(hit_Charge / (theGeometry->WirePitch(hit_View))); // ADC/cm
801 const double dQdX_e(dQdX /
802 (detProp.ElectronsToADC() * settings.m_recombination_factor)); // e/cm
803 const double dEdX(settings.m_useBirksCorrection ?
804 detProp.BirksCorrection(dQdX_e) :
805 dQdX_e * 1000. / util::kGeVToElectrons); // MeV/cm
806 double mips(dEdX / settings.m_dEdX_mip);
807
808 if (mips < 0.) mips = settings.m_mips_if_negative;
809
810 if (mips > settings.m_mips_max) mips = settings.m_mips_max;
811
812 return mips;
813 }
814
815 //------------------------------------------------------------------------------------------------------------------------------------------
816
818 {
819 // QGSP_BERT and EM standard physics list mappings
820 processMap["unknown"] = lar_content::MC_PROC_UNKNOWN;
821 processMap["primary"] = lar_content::MC_PROC_PRIMARY;
822 processMap["compt"] = lar_content::MC_PROC_COMPT;
823 processMap["phot"] = lar_content::MC_PROC_PHOT;
824 processMap["annihil"] = lar_content::MC_PROC_ANNIHIL;
825 processMap["eIoni"] = lar_content::MC_PROC_E_IONI;
826 processMap["eBrem"] = lar_content::MC_PROC_E_BREM;
827 processMap["conv"] = lar_content::MC_PROC_CONV;
828 processMap["muIoni"] = lar_content::MC_PROC_MU_IONI;
829 processMap["muMinusCaptureAtRest"] = lar_content::MC_PROC_MU_MINUS_CAPTURE_AT_REST;
830 processMap["neutronInelastic"] = lar_content::MC_PROC_NEUTRON_INELASTIC;
831 processMap["nCapture"] = lar_content::MC_PROC_N_CAPTURE;
832 processMap["hadElastic"] = lar_content::MC_PROC_HAD_ELASTIC;
833 processMap["Decay"] = lar_content::MC_PROC_DECAY;
834 processMap["CoulombScat"] = lar_content::MC_PROC_COULOMB_SCAT;
835 processMap["muBrems"] = lar_content::MC_PROC_MU_BREM;
836 processMap["muPairProd"] = lar_content::MC_PROC_MU_PAIR_PROD;
837 processMap["PhotonInelastic"] = lar_content::MC_PROC_PHOTON_INELASTIC;
838 processMap["hIoni"] = lar_content::MC_PROC_HAD_IONI;
839 processMap["protonInelastic"] = lar_content::MC_PROC_PROTON_INELASTIC;
840 processMap["pi+Inelastic"] = lar_content::MC_PROC_PI_PLUS_INELASTIC;
841 processMap["CHIPSNuclearCaptureAtRest"] = lar_content::MC_PROC_CHIPS_NUCLEAR_CAPTURE_AT_REST;
842 processMap["pi-Inelastic"] = lar_content::MC_PROC_PI_MINUS_INELASTIC;
843 processMap["Transportation"] = lar_content::MC_PROC_TRANSPORTATION;
844 processMap["Rayl"] = lar_content::MC_PROC_RAYLEIGH;
845 processMap["hBrems"] = lar_content::MC_PROC_HAD_BREM;
846 processMap["hPairProd"] = lar_content::MC_PROC_HAD_PAIR_PROD;
847 processMap["ionIoni"] = lar_content::MC_PROC_ION_IONI;
848 processMap["nKiller"] = lar_content::MC_PROC_NEUTRON_KILLER;
849 processMap["ionInelastic"] = lar_content::MC_PROC_ION_INELASTIC;
850 processMap["He3Inelastic"] = lar_content::MC_PROC_HE3_INELASTIC;
851 processMap["alphaInelastic"] = lar_content::MC_PROC_ALPHA_INELASTIC;
852 processMap["anti_He3Inelastic"] = lar_content::MC_PROC_ANTI_HE3_INELASTIC;
853 processMap["anti_alphaInelastic"] = lar_content::MC_PROC_ANTI_ALPHA_INELASTIC;
854 processMap["hFritiofCaptureAtRest"] = lar_content::MC_PROC_HAD_FRITIOF_CAPTURE_AT_REST;
855 processMap["anti_deuteronInelastic"] = lar_content::MC_PROC_ANTI_DEUTERON_INELASTIC;
856 processMap["anti_neutronInelastic"] = lar_content::MC_PROC_ANTI_NEUTRON_INELASTIC;
857 processMap["anti_protonInelastic"] = lar_content::MC_PROC_ANTI_PROTON_INELASTIC;
858 processMap["anti_tritonInelastic"] = lar_content::MC_PROC_ANTI_TRITON_INELASTIC;
859 processMap["dInelastic"] = lar_content::MC_PROC_DEUTERON_INELASTIC;
860 processMap["electronNuclear"] = lar_content::MC_PROC_ELECTRON_NUCLEAR;
861 processMap["photonNuclear"] = lar_content::MC_PROC_PHOTON_NUCLEAR;
862 processMap["kaon+Inelastic"] = lar_content::MC_PROC_KAON_PLUS_INELASTIC;
863 processMap["kaon-Inelastic"] = lar_content::MC_PROC_KAON_MINUS_INELASTIC;
864 processMap["hBertiniCaptureAtRest"] = lar_content::MC_PROC_HAD_BERTINI_CAPTURE_AT_REST;
865 processMap["lambdaInelastic"] = lar_content::MC_PROC_LAMBDA_INELASTIC;
866 processMap["muonNuclear"] = lar_content::MC_PROC_MU_NUCLEAR;
867 processMap["tInelastic"] = lar_content::MC_PROC_TRITON_INELASTIC;
868 processMap["primaryBackground"] = lar_content::MC_PROC_PRIMARY_BACKGROUND;
869 }
870
871 //------------------------------------------------------------------------------------------------------------------------------------------
872 //------------------------------------------------------------------------------------------------------------------------------------------
873
875 : m_pPrimaryPandora(nullptr)
876 , m_useHitWidths(true)
877 , m_useBirksCorrection(false)
878 , m_useActiveBoundingBox(false)
879 , m_uidOffset(100000000)
880 , m_hitCounterOffset(0)
881 , m_dx_cm(0.5)
882 , m_int_cm(84.)
883 , m_rad_cm(14.)
884 , m_dEdX_mip(2.)
885 , m_mips_max(50.)
886 , m_mips_if_negative(0.)
887 , m_mips_to_gev(3.5e-4)
888 , m_recombination_factor(0.63)
889 {}
890
891} // namespace lar_pandora
Interface class for LArPandora producer modules, which reconstruct recob::PFParticles from recob::Hit...
Header file for the lar calo hit class.
Helper functions for extracting detector geometry for use in reconsruction.
Helper functions for providing inputs to pandora.
Header file for the lar transformation plugin interface class.
Header file for the pandora api class.
Header file for the pandora plugin manager class.
#define PANDORA_THROW_RESULT_IF(StatusCode1, Operator, Command)
Definition StatusCodes.h:43
static pandora::StatusCode SetMCParentDaughterRelationship(const pandora::Pandora &pandora, const void *const pParentAddress, const void *const pDaughterAddress)
Set parent-daughter mc particle relationship.
Definition PandoraApi.cc:44
static pandora::StatusCode SetCaloHitToMCParticleRelationship(const pandora::Pandora &pandora, const void *const pCaloHitParentAddress, const void *const pMCParticleParentAddress, const float mcParticleWeight=1)
Set calo hit to mc particle relationship.
Definition PandoraApi.cc:68
LArCaloHitFactory responsible for object creation.
Definition LArCaloHit.h:111
LAr calo hit parameters.
Definition LArCaloHit.h:28
pandora::InputUInt m_daughterVolumeId
The daughter volume id.
Definition LArCaloHit.h:31
pandora::InputUInt m_larTPCVolumeId
The lar tpc volume id.
Definition LArCaloHit.h:30
LArMCParticleFactory responsible for object creation.
LAr mc particle parameters.
pandora::InputInt m_nuanceCode
The nuance code.
pandora::InputInt m_process
The process creating the particle.
drift volume class to hold properties of drift volume
drift volume class to hold properties of drift volume
Empty interface to map pandora to specifics in the LArSoft geometry.
virtual PandoraApi::Geometry::LineGap::Parameters CreateLineGapParametersFromDetectorGaps(const LArDetectorGap &gap) const =0
Create the line gap parameters to give to the pandora API.
virtual PandoraApi::Geometry::LineGap::Parameters CreateLineGapParametersFromReadoutGaps(const geo::View_t view, const geo::TPCID::TPCID_t tpc, const geo::CryostatID::CryostatID_t cstat, const geo::Point_t &firstXYZ, const geo::Point_t &lastXYZ, const float halfWirePitch, const float xFirst, const float xLast, const pandora::Pandora *pPandora) const =0
Create the line gap parameters to give to the pandora API.
virtual geo::View_t TargetViewU(const geo::TPCID::TPCID_t tpc, const geo::CryostatID::CryostatID_t cstat) const =0
Map a LArSoft view to Pandora's U view.
virtual geo::View_t TargetViewV(const geo::TPCID::TPCID_t tpc, const geo::CryostatID::CryostatID_t cstat) const =0
Map a LArSoft view to Pandora's V view.
virtual geo::View_t TargetViewW(const geo::TPCID::TPCID_t tpc, const geo::CryostatID::CryostatID_t cstat) const =0
Map a LArSoft view to Pandora's W view.
static unsigned int GetDaughterVolumeID(const LArDriftVolumeMap &driftVolumeMap, const unsigned int cstat, const unsigned int tpc)
Get daughter volume ID from a specified cryostat/tpc pair.
static unsigned int GetVolumeID(const LArDriftVolumeMap &driftVolumeMap, const unsigned int cstat, const unsigned int tpc)
Get drift volume ID from a specified cryostat/tpc pair.
const pandora::Pandora * m_pPrimaryPandora
static void CreatePandoraReadoutGaps(const Settings &settings, const LArDriftVolumeMap &driftVolumeMap)
Create pandora line gaps to cover any (continuous regions of) bad channels.
static bool IsPrimaryMCParticle(const art::Ptr< simb::MCParticle > &mcParticle, std::map< const simb::MCParticle, bool > &primaryMCParticleMap)
Check whether an MCParticle can be found in a given map.
static void CreatePandoraDetectorGaps(const Settings &settings, const LArDriftVolumeList &driftVolumeList, const LArDetectorGapList &listOfGaps)
Create pandora line gaps to cover dead regions between TPCs in a global drift volume approach.
static float GetTrueX0(const art::Event &evt, const art::Ptr< simb::MCParticle > &particle, const int nT)
Use detector and time services to get a true X offset for a given trajectory point.
static void CreatePandoraHits2D(const art::Event &evt, const Settings &settings, const LArDriftVolumeMap &driftVolumeMap, const HitVector &hitVector, IdToHitMap &idToHitMap)
Create the Pandora 2D hits from the ART hits.
std::map< std::string, lar_content::MCProcess > MCProcessMap
static double GetMips(const detinfo::DetectorPropertiesData &detProp, const Settings &settings, const double hit_Charge, const geo::View_t hit_View)
Convert charge in ADCs to approximate MIPs.
static void CreatePandoraMCLinks2D(const Settings &settings, const HitMap &hitMap, const HitsToTrackIDEs &hitToParticleMap)
Create links between the 2D hits and Pandora MC particles.
static void FillMCProcessMap(MCProcessMap &processMap)
Populate a map from MC process string to enumeration.
static void FindPrimaryParticles(const RawMCParticleVector &mcParticleVector, std::map< const simb::MCParticle, bool > &primaryMCParticleMap)
Find all primary MCParticles in a given vector of MCParticles.
static void CreatePandoraMCParticles(const Settings &settings, const MCTruthToMCParticles &truthToParticles, const MCParticlesToMCTruth &particlesToTruth, const RawMCParticleVector &generatorMCParticleVector)
Create the Pandora MC particles from the MC particles.
static void GetTrueStartAndEndPoints(const Settings &settings, const art::Ptr< simb::MCParticle > &particle, int &startT, int &endT)
Loop over MC trajectory points and identify start and end points within the detector.
static void CreatePandoraLArTPCs(const Settings &settings, const LArDriftVolumeList &driftVolumeList)
Create pandora LArTPCs to represent the different drift volumes in use.
static pandora::StatusCode Create(const pandora::Pandora &pandora, const Parameters &parameters, const pandora::ObjectFactory< Parameters, Object > &factory=pandora::PandoraObjectFactory< Parameters, Object >())
Create a new object from a user factory.
CartesianVector class.
virtual double YZtoU(const double y, const double z) const =0
Transform from (Y,Z) to U position.
virtual double YZtoW(const double y, const double z) const =0
Transform from (Y,Z) to W position.
virtual double YZtoV(const double y, const double z) const =0
Transform from (Y,Z) to V position.
Pandora class.
Definition Pandora.h:40
const PluginManager * GetPlugins() const
Get the pandora plugin instance, providing access to user registered functions and calculators.
Definition Pandora.cc:196
const LArTransformationPlugin * GetLArTransformationPlugin() const
Get the address of the lar transformation plugin.
StatusCodeException class.
@ MC_PROC_HAD_BERTINI_CAPTURE_AT_REST
@ MC_PROC_NEUTRON_INELASTIC
@ MC_PROC_ALPHA_INELASTIC
@ MC_PROC_DEUTERON_INELASTIC
@ MC_PROC_TRITON_INELASTIC
@ MC_PROC_ANTI_DEUTERON_INELASTIC
@ MC_PROC_PI_PLUS_INELASTIC
@ MC_PROC_ANTI_NEUTRON_INELASTIC
@ MC_PROC_ANTI_ALPHA_INELASTIC
@ MC_PROC_CHIPS_NUCLEAR_CAPTURE_AT_REST
@ MC_PROC_PRIMARY_BACKGROUND
@ MC_PROC_ANTI_TRITON_INELASTIC
@ MC_PROC_ELECTRON_NUCLEAR
@ MC_PROC_ANTI_PROTON_INELASTIC
@ MC_PROC_PI_MINUS_INELASTIC
@ MC_PROC_PROTON_INELASTIC
@ MC_PROC_PHOTON_INELASTIC
@ MC_PROC_MU_MINUS_CAPTURE_AT_REST
@ MC_PROC_HAD_FRITIOF_CAPTURE_AT_REST
@ MC_PROC_KAON_PLUS_INELASTIC
@ MC_PROC_LAMBDA_INELASTIC
@ MC_PROC_ANTI_HE3_INELASTIC
@ MC_PROC_KAON_MINUS_INELASTIC
LArPandoraDetectorType * GetDetectorType()
Factory class that returns the correct detector type interface.
std::map< int, art::Ptr< simb::MCParticle > > MCParticleMap
std::map< art::Ptr< simb::MCParticle >, art::Ptr< simb::MCTruth > > MCParticlesToMCTruth
std::vector< sim::TrackIDE > TrackIDEVector
std::map< art::Ptr< recob::Hit >, TrackIDEVector > HitsToTrackIDEs
std::vector< simb::MCParticle > RawMCParticleVector
std::map< art::Ptr< simb::MCTruth >, MCParticleVector > MCTruthToMCParticles
std::vector< art::Ptr< recob::Hit > > HitVector
std::map< int, art::Ptr< recob::Hit > > IdToHitMap
Definition ILArPandora.h:24
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
std::vector< LArDetectorGap > LArDetectorGapList
std::vector< LArDriftVolume > LArDriftVolumeList
std::map< unsigned int, LArDriftVolume > LArDriftVolumeMap