Merge branch 'rricesmi/calypso-master' into 'master'

Example job to read Tracker RDO data and associate with Truth

See merge request faser/calypso!54

Control/CalypsoExample/RDOReadExample/CMakeLists.txt

+atlas_subdir( RDOReadExample )
+
+atlas_depends_on_subdirs( PRIVATE
+                    Generators/GeneratorObjects
+                    Control/AthenaBaseComps
+                    Tracker/TrackerSimEvent
+		    Tracker/TrackerRawEvent/TrackerRawData
+		    Tracker/TrackerRawEvent/TrackerSimData
+		    Tracker/TrackerDetDescr/TrackerIdentifier
+        )
+
+atlas_add_component( RDOReadExample
+                    src/RDOReadAlg.cxx
+                    src/components/RDOReadExample_entries.cxx
+                    LINK_LIBRARIES AthenaBaseComps GeneratorObjects TrackerSimEvent TrackerRawData TrackerSimData TrackerIdentifier
+        )
+
+#atlas_install_joboptions( share/*.py )
+atlas_install_python_modules( python/*.py )

Control/CalypsoExample/RDOReadExample/README.md

+Example algorithm to access RDO data & its corresponding simulation data (given generated from simulation)
+
+For each RDO, uses the identifier and TrackerSimDataCollection to get the TrackerSimData associated with it
+
+Then loops through the energy deposits of that TrackerSimData to find the particle barcode of the highest energy deposit
+
+Then find the FaserSiHit corresponding to the particle with the same barcode collision at the same place as the RDO
+
+And plots the RDO GroupSize versus the incident angle the particle made when it was detected
+
+Currently reads from a RDO root file named myRDO.pool.root, which can be generated from the digitization:
+https://gitlab.cern.ch/faser/calypso/tree/master/Tracker/TrackerDigitization/FaserSCT_Digitization
+
+After calypso is installed and compiled & digitization generated
+
+run > python python/RDOReadExample/RDOReadExampleConfig.py
\ No newline at end of file

Control/CalypsoExample/RDOReadExample/python/RDOReadExampleConfig.py

+#!/usr/bin/env python
+import sys
+from AthenaConfiguration.ComponentFactory import CompFactory
+
+def RDOReadExampleCfg(flags, name="RDOReadExampleAlg", **kwargs):
+
+    from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg
+    a = FaserGeometryCfg(flags)
+
+#    from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg
+#    a.merge(MagneticFieldSvcCfg(flags))
+#    fieldSvc = a.getService("FaserFieldSvc")
+
+    from RDOReadExample.RDOReadExampleConf import RDOReadAlg
+    a.addEventAlgo(RDOReadAlg(name, **kwargs))
+#    a.getEventAlgo(name).FieldService = fieldSvc
+
+    thistSvc = CompFactory.THistSvc()
+    thistSvc.Output += ["HIST DATAFILE='rdoReadHist.root' OPT='RECREATE'"]
+    a.addService(thistSvc)
+
+
+    return a
+
+
+if __name__ == "__main__":
+    # from AthenaCommon.Logging import log, logging
+    from AthenaCommon.Constants import VERBOSE, INFO
+    from AthenaCommon.Configurable import Configurable
+
+    from CalypsoConfiguration.AllConfigFlags import ConfigFlags
+
+    Configurable.configurableRun3Behavior = True
+    
+# Flags for this job
+    ConfigFlags.Input.isMC = True                                # Needed to bypass autoconfig
+    ConfigFlags.IOVDb.GlobalTag = "OFLCOND-XXXX-XXX-XX"          # Needed to bypass autoconfig, only the "OFLCOND" matters at the moment
+    ConfigFlags.GeoModel.FaserVersion     = "FASER-00"           # Default FASER geometry
+    ConfigFlags.Detector.SimulateFaserSCT = True
+    ConfigFlags.Input.Files = ["myRDO.pool.root"]
+    ConfigFlags.lock()
+
+# Configure components
+    from AthenaConfiguration.MainServicesConfig import MainServicesSerialCfg
+    from AthenaPoolCnvSvc.PoolReadConfig import PoolReadCfg
+    acc = MainServicesSerialCfg()
+    acc.merge(PoolReadCfg(ConfigFlags))
+
+# Set things up to create a conditions DB with neutral Tracker alignment transforms
+    acc.merge(RDOReadExampleCfg(ConfigFlags))
+    # acc.addService(CompFactory.IOVRegistrationSvc(PayloadTable=False))
+    # acc.getService("IOVRegistrationSvc").OutputLevel = VERBOSE
+    # acc.getService("IOVDbSvc").dbConnection = "sqlite://;schema=ALLP200.db;dbname=OFLP200" 
+
+# Configure verbosity    
+    # ConfigFlags.dump()
+    # logging.getLogger('forcomps').setLevel(VERBOSE)
+    # acc.foreach_component("*").OutputLevel = VERBOSE
+    # acc.foreach_component("*ClassID*").OutputLevel = INFO
+    # log.setLevel(VERBOSE)
+    
+# Execute and finish
+
+
+    sys.exit(int(acc.run(maxEvents=-1).isFailure()))

Control/CalypsoExample/RDOReadExample/python/__init__.py

+#  __author__ = 'Ryan Rice-Smith'

Control/CalypsoExample/RDOReadExample/src/RDOReadAlg.cxx

+#include "RDOReadAlg.h"
+#include "TrackerSimEvent/FaserSiHitIdHelper.h"
+#include "StoreGate/StoreGateSvc.h"
+#include "StoreGate/StoreGate.h"
+#include "TrackerIdentifier/FaserSCT_ID.h"
+
+
+RDOReadAlg::RDOReadAlg(const std::string& name, ISvcLocator* pSvcLocator)
+: AthHistogramAlgorithm(name, pSvcLocator) { m_hist = nullptr; }
+
+RDOReadAlg::~RDOReadAlg() { }
+
+StatusCode RDOReadAlg::initialize()
+{
+    // initialize a histogram 
+    // letter at end of TH1 indicated variable type (D double, F float etc)
+    m_hist = new TH1D("GroupSize", "RDO Group Size", 8, 0, 8); //first string is root object name, second is histogram title
+    m_hprof = new TProfile("IncAngleGroup", "Mean Group Size vs Incident Angle", 10, -.1, .1 ,0,5);
+    m_incAnglHist = new TH1D("IncAngleHist", "Incident Angle Count", 10, -.1, .1);
+    ATH_CHECK(histSvc()->regHist("/HIST/myhist", m_hist));
+    ATH_CHECK(histSvc()->regHist("/HIST/myhistprof", m_hprof));
+    ATH_CHECK(histSvc()->regHist("/HIST/myhistAngl", m_incAnglHist));
+
+    // initialize data handle keys
+    ATH_CHECK( m_mcEventKey.initialize() );
+    ATH_CHECK( m_faserSiHitKey.initialize() );
+    ATH_CHECK( m_faserRdoKey.initialize());
+    ATH_CHECK( m_sctMap.initialize());
+    ATH_MSG_INFO( "Using GenEvent collection with key " << m_mcEventKey.key());
+    ATH_MSG_INFO( "Using Faser SiHit collection with key " << m_faserSiHitKey.key());
+    ATH_MSG_INFO( "Using FaserSCT RDO Container with key " << m_faserRdoKey.key());//works
+    ATH_MSG_INFO( "Using SCT_SDO_Map with key "<< m_sctMap.key());
+    return StatusCode::SUCCESS;
+}
+
+StatusCode RDOReadAlg::execute()
+{
+    // Handles created from handle keys behave like pointers to the corresponding container
+    SG::ReadHandle<McEventCollection> h_mcEvents(m_mcEventKey);
+    ATH_MSG_INFO("Read McEventContainer with " << h_mcEvents->size() << " events");
+    if (h_mcEvents->size() == 0) return StatusCode::FAILURE;
+
+    SG::ReadHandle<FaserSiHitCollection> h_siHits(m_faserSiHitKey);
+    ATH_MSG_INFO("Read FaserSiHitCollection with " << h_siHits->size() << " hits");
+
+    SG::ReadHandle<FaserSCT_RDO_Container> h_sctRDO(m_faserRdoKey);
+    SG::ReadHandle<TrackerSimDataCollection> h_collectionMap(m_sctMap);
+
+    //Looping through RDO's
+    for( const auto& collection : *h_sctRDO)
+    {
+        for(const auto& rawdata : *collection)
+        {
+            //Each RDO has an identifier connecting it to the simulation data that made it
+            //Allows access of the simulation data that made that RDO
+            const auto identifier = rawdata->identify();
+            ATH_MSG_INFO("map size "<<h_collectionMap->size());
+            if( h_collectionMap->count(identifier) == 0)
+            {
+                ATH_MSG_INFO("no map found w/identifier "<<identifier);
+                continue;
+            }
+
+            //Collection map takes identifier and returns simulation data
+            const auto& simdata = h_collectionMap->find(rawdata->identify())->second;
+            const auto& deposits = simdata.getdeposits();
+//            ATH_MSG_INFO("deposits size "<<deposits.size());
+
+//            ATH_MSG_INFO("identifier: "<<rawdata->identify());
+//            ATH_MSG_INFO("group size of "<<rawdata->getGroupSize());
+            m_hist->Fill(rawdata->getGroupSize());
+
+            //loop through deposits to find one w/ highest energy & get barcode
+            float highestDep = 0;
+            int barcode = 0;
+            for( const auto& depositPair : deposits)
+            {
+                if( depositPair.second > highestDep)
+                {
+                    highestDep = depositPair.second;
+                    barcode = depositPair.first->barcode();
+                    depositPair.first->print(std::cout);
+                    ATH_MSG_INFO("pdg id "<<depositPair.first->pdg_id());
+                }
+            }
+            ATH_MSG_INFO("final barcode of: "<<barcode);
+
+
+            //Helper function to get hit location information from RDO identifier
+            const FaserSCT_ID* pix;
+            StoreGateSvc* detStore(nullptr);
+            {
+                detStore = StoreGate::pointer("DetectorStore");
+                if (detStore->retrieve(pix, "FaserSCT_ID").isFailure()) { pix = 0; }
+            }
+            int station = pix->station(identifier);
+            int plane = pix->layer(identifier);
+            int row = pix->phi_module(identifier);
+            int module = pix->eta_module(identifier);
+            int sensor = pix->side(identifier);
+
+            ATH_MSG_INFO("trying to match hit to stat/plane/row/mod/sens: "<<station<<" "<<plane<<" "<<row<<" "<<module<<" "<<sensor);
+            for (const FaserSiHit& hit : *h_siHits)
+            {
+//                ATH_MSG_INFO("hit w/vals "<<hit.getStation()<<" "<<hit.getPlane()<<" "<<hit.getRow()<<" "<<hit.getModule()<<" "<<hit.getSensor()<<" barcode: "<<hit.trackNumber());
+                //set of conditions to confirm looking at same particle in same place for SiHit as RDO
+                if(hit.getStation() == station 
+                    && hit.getPlane() == plane
+                    && hit.getRow() == row
+                    && hit.getModule() == module
+                    && hit.getSensor() == sensor
+                    && hit.trackNumber() == barcode)
+                {
+                    ATH_MSG_INFO("matched particle and plotting w/ barcode "<<barcode);
+                    //here we plot point of angle vs countsize!
+                    float delx = hit.localEndPosition().x() - hit.localStartPosition().x();
+                    float dely = hit.localEndPosition().y() - hit.localStartPosition().y();
+                    float delz = hit.localEndPosition().z() - hit.localStartPosition().z();
+                    float norm = sqrt(delx*delx + dely*dely + delz*delz);
+                    ATH_MSG_INFO("acos val and group val of "<<acos(abs(delx)/norm)<<" "<<rawdata->getGroupSize());
+                    float ang = acos(delx/norm);
+                    if(ang > 1.5)
+                    {
+                        ang = ang - 3.1415;
+                    }
+                    m_hprof->Fill(ang, rawdata->getGroupSize(),1);
+                    m_incAnglHist->Fill(ang);
+                    break;
+                }
+            }
+        }
+    }
+
+    return StatusCode::SUCCESS;
+}
+
+StatusCode RDOReadAlg::finalize()
+{
+    return StatusCode::SUCCESS;
+}
\ No newline at end of file

Control/CalypsoExample/RDOReadExample/src/RDOReadAlg.h

+#include "AthenaBaseComps/AthHistogramAlgorithm.h"
+#include "GeneratorObjects/McEventCollection.h"
+#include "TrackerSimEvent/FaserSiHitCollection.h"
+#include "TrackerRawData/FaserSCT_RDO_Container.h"
+#include "TrackerSimData/TrackerSimDataCollection.h"
+#include <TH1.h>
+#include <math.h>
+#include <TProfile.h>
+
+/* RDORead reading example - Ryan Rice-Smith, UC Irvine */
+
+class RDOReadAlg : public AthHistogramAlgorithm
+{
+    public:
+    RDOReadAlg(const std::string& name, ISvcLocator* pSvcLocator);
+
+    virtual ~RDOReadAlg();
+
+    StatusCode initialize();
+    StatusCode execute();
+    StatusCode finalize();
+
+    private:
+    TH1* m_hist;  // Example histogram
+    TH1* m_incAnglHist;
+    TProfile* m_hprof;
+
+    // Read handle keys for data containers
+    // Any other event data can be accessed identically
+    // Note the key names ("GEN_EVENT" or "SCT_Hits") are Gaudi properties and can be configured at run-time
+    SG::ReadHandleKey<McEventCollection> m_mcEventKey       { this, "McEventCollection", "GEN_EVENT" };
+    SG::ReadHandleKey<FaserSiHitCollection> m_faserSiHitKey { this, "FaserSiHitCollection", "SCT_Hits" };
+    SG::ReadHandleKey<FaserSCT_RDO_Container> m_faserRdoKey { this, "FaserSCT_RDO_Container", "SCT_RDOs"};
+    SG::ReadHandleKey<TrackerSimDataCollection> m_sctMap {this, "TrackerSimDataCollection", "SCT_SDO_Map"};
+};
\ No newline at end of file

Control/CalypsoExample/RDOReadExample/src/components/RDOReadExample_entries.cxx

+#include "../RDOReadAlg.h"
+
+DECLARE_COMPONENT( RDOReadAlg ) //change alg name if update file/class etc
\ No newline at end of file