Back to home page

EIC code displayed by LXR

 
 

    


Warning, /tutorial-jana2/episodes/04-factory.md is written in an unsupported language. File is not indexed.

0001 ---
0002 title: "Creating or modifying a JANA factory in order to implement a reconstruction algorithm"
0003 teaching: 15
0004 exercises: 20
0005 ---
0006 
0007 ::::::::::::::::::::::::::::::::::::::::::::: questions
0008 
0009 - How to write a reconstruction algorithm in EICrecon?
0010 
0011 :::::::::::::::::::::::::::::::::::::::::::::
0012 
0013 ::::::::::::::::::::::::::::::::::::::::::::: objectives
0014 
0015 - Learn how to create a new factory in EICrecon that supplies a reconstruction algorithm for all to use.
0016 - Understand the directory structure for where the factory should be placed in the source tree.
0017 - Understand how to use a generic algorithm in a JANA factory.
0018 
0019 :::::::::::::::::::::::::::::::::::::::::::::
0020 
0021 ::::::::::::::::::::::::::::::::::::::::::::: callout
0022 
0023 Note: The following episode presents a somewhat outdated view, and some commands may not function.
0024 If you are only interested in analyzing already-reconstructed data, then there is no requirement
0025 to use a plugin as described below; the output ROOT file works too.
0026 
0027 :::::::::::::::::::::::::::::::::::::::::::::
0028 
0029 ## Introduction
0030 
0031 Now that you've learned about JANA plugins and JEventProcessors, let's talk about JFactories. JFactories are another essential JANA component just like JEventProcessors and JEventSources. While JEventProcessors are used for _aggregating_ results from each event into a structured output such as a histogram or a file, JFactories are used for computing those results in an organized way.
0032 
0033 ### When do I use a JFactory?
0034 
0035 - If you have an input file and need to read data model objects from it, use a JEventSource.
0036 - If you have an output file (or histogram) and wish to write data model objects to it, use a JEventProcessor.
0037 - If you have some data model objects and wish to produce a new data model object, use a JFactory.
0038 
0039 ### Why should I prefer writing a JFactory?
0040 
0041 1. They make your code reusable. Different people can use your results later without having to understand the specifics of what you did.
0042 
0043 2. If you are consuming some data which doesn't look right to you, JFactories make it extremely easy to pinpoint exactly which code produced this data.
0044 
0045 3. EICrecon needs to run multithreaded, and using JFactories can help steer you away from introducing thorny parallelism bugs.
0046 
0047 4. You can simply ask for the results you need and the JFactory will provide it. If nobody needs the results from the JFactory, it won't be run. If the results were already in the input file, it won't be run. If there are multiple consumers, the results are only computed once and then cached. If the JFactory relies on results from other JFactories, it will call them transparently and recursively.
0048 
0049 ### When do I create my own plugin?
0050 
0051 - If you are doing a one-off prototype, it's fine to just use a ROOT macro.
0052 - If you are writing code you'll probably return to, we recommend putting the code in a standalone (i.e. outside of the EICrecon source tree) plugin.
0053 - If you are writing code other people will probably want to run, we recommend adding your plugin to the EICrecon source tree.
0054 - If you are writing a JFactory, we recommend adding it to the EICrecon source tree, either to an existing plugin or to a new one.
0055 
0056 ## Algorithms vs Factories
0057 
0058 In general, a Factory is a programming pattern for constructing objects in an abstract way. Oftentimes, the Factory is calling an algorithm under the hood.
0059 This algorithm may be very generic. For instance, we may have a Factory that produces Cluster objects for a barrel calorimeter, and it calls a clustering algorithm
0060 that doesn't care at all about barrel calorimeters, just the position and energy of energy of each CalorimeterHit object. Perhaps multiple factories for creating clusters
0061 for completely different detectors are all using the same algorithm.
0062 
0063 Note that Gaudi provides an abstraction called "Algorithm" which is essentially its own version of a JFactory. In EICrecon, we have been separating out _generic algorithms_ from the old Gaudi and new JANA code so that these can be developed and tested independently. To see an example of how a generic algorithm is being implemented, look at these examples:
0064 
0065 ```
0066 src/detectors/EEMC/RawCalorimeterHit_factory_EcalEndcapNRawHits.h
0067 src/algorithms/calorimetry/CalorimeterHitDigi.h
0068 src/algorithms/calorimetry/CalorimeterHitDigi.cc
0069 ```
0070 
0071 Using generic algorithms makes things slightly more complex. However, the generic algorithms can be recycled for use in multiple detector systems which adds some simplification.
0072 
0073 ## Parallelism considerations
0074 
0075 JEventProcessors observe the entire _event stream_, and require a _critical section_ where only one thread is allowed to modify a shared resource (such as a histogram) at any time.
0076 JFactories, on the other hand, only observe a single event at a time, and work on each event independently. Each worker thread is given an independent event with its own set of factories. This means that for a given JFactory instance, there will be only one thread working on one event at any time. You get the benefits of multithreading _without_ having to make each JFactory thread-safe.
0077 
0078 You can write JFactories in an almost-functional style, but you can also cache some data on the JFactory that will stick around from event-to-event. This is useful for things like conditions and geometry data, where for performance reasons you don't want to be doing a deep lookup on every event. Instead, you can write callbacks such as `BeginRun()`, where you can update your cached values when the run number changes.
0079 
0080 Note that just because the JFactory _can_ be called in parallel doesn't mean it always will. If you call event->Get() from inside `JEventProcessor::ProcessSequential`, in particular, the factory will run single-threaded and slow everything down. However, if you call it using `Prefetch` instead, it will run in parallel and you may get a speed boost.
0081 
0082 ## How do I use an existing JFactory?
0083 
0084 Using an existing JFactory is extremely easy! Any time you are someplace where you have access to a `JEvent` object, do this:
0085 
0086 ```c++
0087 auto clusters = event->Get<edm4eic::Cluster>("EcalEndcapNIslandClusters");
0088 
0089 for (auto c : clusters) {
0090   // ... do something with a cluster
0091 }
0092 ```
0093 
0094 As you can see, it doesn't matter whether the `Cluster` objects were calculated from some simpler objects, or were simply loaded from a file. This is a very powerful concept.
0095 
0096 One thing we might want to do is to swap one factory for another, possibly even at runtime. This is easy to do if you just make the factory tag be a parameter:
0097 
0098 ```c++
0099 std::string my_cluster_source = "EcalEndcapNIslandClusters";  // Make this be a parameter
0100 app->SetDefaultParameter("MyPlugin:MyAnalysis:my_cluster_source", my_cluster_source, "Cluster source for MyAnalysis");
0101 auto clusters = event->Get<edm4eic::Cluster>(my_cluster_source);
0102 ```
0103 
0104 ## How do I create a new JFactory?
0105 
0106 We are going to add a new JFactory inside EICrecon.
0107 
0108 `src/detectors/EEMC/Cluster_factory_EcalEndcapNIslandClusters.h`:
0109 
0110 ```c++
0111 #pragma once
0112 
0113 #include <edm4eic/Cluster.h>
0114 #include <JANA/JFactoryT.h>
0115 #include <services/log/Log_service.h>
0116 
0117 class Cluster_factory_EcalEndcapNIslandClusters : public JFactoryT<edm4eic::Cluster> {
0118 public:
0119 
0120     Cluster_factory_EcalEndcapNIslandClusters(); // Constructor
0121 
0122     void Init() override;
0123     // Gets called exactly once at the beginning of the JFactory's life
0124 
0125     void ChangeRun(const std::shared_ptr<const JEvent> &event) override {};
0126     // Gets called on events where the run number has changed (before Process())
0127 
0128     void Process(const std::shared_ptr<const JEvent> &event) override;
0129     // Gets called on every event
0130 
0131     void Finish() override {};
0132     // Gets called exactly once at the end of the JFactory's life
0133 
0134 private:
0135     float m_scaleFactor;
0136     std::shared_ptr<spdlog::logger> m_log;
0137 
0138 };
0139 ```
0140 
0141 `src/detectors/EEMC/Cluster_factory_EcalEndcapNIslandClusters.cc`:
0142 
0143 ```c++
0144 #include "Cluster_factory_EcalEndcapNIslandClusters.h"
0145 
0146 #include <edm4eic/ProtoCluster.h>
0147 #include <JANA/JEvent.h>
0148 
0149 
0150 Cluster_factory_EcalEndcapNIslandClusters::Cluster_factory_EcalEndcapNIslandClusters() {
0151 
0152     SetTag("EcalEndcapNIslandClusters");
0153 }
0154 
0155 
0156 void Cluster_factory_EcalEndcapNIslandClusters::Init() {
0157     auto app = GetApplication();
0158 
0159     // This is an example of how to declare a configuration parameter that
0160     // can be set at run time. e.g. with -PEEMC:EcalEndcapNIslandClusters:scaleFactor=0.97
0161     m_scaleFactor =0.98;
0162     app->SetDefaultParameter("EEMC:EcalEndcapNIslandClusters:scaleFactor", m_scaleFactor, "Energy scale factor");
0163 
0164     // This is how you access shared resources using the JService interface
0165     m_log = app->GetService<Log_service>()->logger("EcalEndcapNIslandClusters");
0166 }
0167 
0168 
0169 void Cluster_factory_EcalEndcapNIslandClusters::Process(const std::shared_ptr<const JEvent> &event) {
0170 
0171     m_log->info("Processing event {}", event->GetEventNumber());
0172 
0173     // Grab inputs
0174     auto protoclusters = event->Get<edm4eic::ProtoCluster>("EcalEndcapNIslandProtoClusters");
0175 
0176     // Loop over protoclusters and turn each into a cluster
0177     std::vector<edm4eic::Cluster*> outputClusters;
0178     for( auto proto : protoclusters ) {
0179 
0180         // ======================
0181         // Algorithm goes here!
0182         // ======================
0183 
0184         auto cluster = new edm4eic::Cluster(
0185             0, // type
0186             energy * m_scaleFactor,
0187             sqrt(energyError_squared),
0188             time,
0189             timeError,
0190             proto->hits_size(),
0191             position,
0192             edm4eic::Cov3f(), // positionError,
0193             0.0,              // intrinsicTheta,
0194             0.0,              // intrinsicPhi,
0195             edm4eic::Cov2f()  // intrinsicDirectionError
0196             );
0197 
0198         outputClusters.push_back( cluster );
0199     }
0200 
0201     // Hand ownership of algorithm objects over to JANA
0202     Set(outputClusters);
0203 }
0204 ```
0205 
0206 We can now fill in the algorithm with anything we like!
0207 
0208 ```c++
0209         // Grab inputs
0210         auto protoclusters = event->Get<edm4eic::ProtoCluster>("EcalEndcapNIslandProtoClusters");
0211 
0212         // Loop over protoclusters and turn each into a cluster
0213         std::vector<edm4eic::Cluster*> outputClusters;
0214         for( auto proto : protoclusters ) {
0215 
0216             // Fill cumulative values by looping over all hits in proto cluster
0217             float energy = 0;
0218             double energyError_squared = 0.0;
0219             float time = 1.0E8;
0220             float timeError;
0221             edm4hep::Vector3f position;
0222             double sum_weights = 0.0;
0223             for( uint32_t ihit=0; ihit<proto->hits_size() ; ihit++){
0224                 auto const &hit = proto->getHits(ihit);
0225                 auto weight = proto->getWeights(ihit);
0226                 energy += hit.getEnergy();
0227                 energyError_squared += std::pow(hit.getEnergyError(), 2.0);
0228                 if( hit.getTime() < time ){
0229                     time = hit.getTime();            // use earliest time
0230                     timeError = hit.getTimeError();  // use error of earliest time
0231                 }
0232                 auto &p = hit.getPosition();
0233                 position.x += p.x*weight;
0234                 position.y += p.y*weight;
0235                 position.z += p.z*weight;
0236                 sum_weights += weight;
0237             }
0238 
0239             // Normalize position
0240             position.x /= sum_weights;
0241             position.y /= sum_weights;
0242             position.z /= sum_weights;
0243 
0244             // Create a cluster object from values accumulated from hits above
0245             auto cluster = new edm4eic::Cluster(
0246                 0, // type (?))
0247                 energy * m_scaleFactor,
0248                 sqrt(energyError_squared),
0249                 time,
0250                 timeError,
0251                 proto->hits_size(),
0252                 position,
0253 
0254                 // Not sure how to calculate these last few
0255                 edm4eic::Cov3f(), // positionError,
0256                 0.0, // intrinsicTheta,
0257                 0.0, // intrinsicPhi,
0258                 edm4eic::Cov2f() // intrinsicDirectionError
0259                 );
0260 
0261             outputClusters.push_back( cluster );
0262         }
0263 
0264         // Hand ownership of algorithm objects over to JANA
0265         Set(outputClusters);
0266 ```
0267 
0268 You can't pass JANA a JFactory directly (because it needs to create an arbitrary number of them on the fly). Instead you register a `JFactoryGenerator` object:
0269 
0270 `src/detectors/EEMC/EEMC.cc`
0271 
0272 ```c++
0273 // In your plugin's init
0274 
0275 #include <JANA/JFactoryGenerator.h>
0276 // ...
0277 #include "Cluster_factory_EcalEndcapNIslandClusters.h"
0278 
0279 extern "C" {
0280     void InitPlugin(JApplication *app) {
0281         InitJANAPlugin(app);
0282         // ...
0283 
0284         app->Add(new JFactoryGeneratorT<Cluster_factory_EcalEndcapNIslandClusters>());
0285      }
0286 ```
0287 
0288 Finally, we go ahead and trigger the factory (remember, factories won't do anything unless activated by a JEventProcessor).
0289 
0290 ```bash
0291 eicrecon in.root -Ppodio:output_file=out.root -Ppodio:output_collections=EcalEndcapNIslandClusters -Pjana:nevents=10
0292 ```
0293 
0294 ::::::::::::::::::::::::::::::::::::::::::::: challenge
0295 
0296 ## Exercise
0297 
0298 Your exercise is to get this JFactory working! You can tweak the algorithm, add log messages, add additional config parameters, etc.
0299 
0300 ::::::::::::::: solution
0301 
0302 Place the header and source files under `src/detectors/EEMC/`, register the factory generator in
0303 `EEMC.cc`, rebuild EICrecon, and run the `eicrecon` command above. A successful run writes
0304 `out.root` containing the `EcalEndcapNIslandClusters` collection and prints the `m_log->info`
0305 message for each processed event. Experiment by changing `scaleFactor` on the command line with
0306 `-PEEMC:EcalEndcapNIslandClusters:scaleFactor=0.97`.
0307 
0308 :::::::::::::::
0309 
0310 :::::::::::::::::::::::::::::::::::::::::::::
0311 
0312 ::::::::::::::::::::::::::::::::::::::::::::: keypoints
0313 
0314 - Create a factory for reconstructing single subdetector data or for global reconstruction.
0315 
0316 :::::::::::::::::::::::::::::::::::::::::::::