DDG4/reco/Geant4SurfaceTest.cpp

0001 //  AIDA Detector description implementation
0002 //--------------------------------------------------------------------------
0003 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0004 // All rights reserved.
0005 //
0006 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0007 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0008 //
0009 //  \author Markus Frank
0010 //  \date   2015-10-14
0011 //
0012 //==========================================================================
0013 #ifndef DD4HEP_DDG4_GEANT4SURFACETEST_H
0014 #define DD4HEP_DDG4_GEANT4SURFACETEST_H
0015
0016 // Framework include files
0017 //#include "Geant4SurfaceTest.h"
0018 #include <DDG4/Geant4EventAction.h>
0019
0020 /// Namespace for the AIDA detector description toolkit
0021 namespace dd4hep {
0022
0023   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0024   namespace sim {
0025
0026     /** @class Geant4SurfaceTest Geant4SurfaceTest.h reco/Geant4SurfaceTest.h
0027      */
0028     class Geant4SurfaceTest : public Geant4EventAction     {
0029     public:
0030       /// Standard constructor
0031       Geant4SurfaceTest(Geant4Context* context, const std::string& nam);
0032       /// Destructor
0033       virtual ~Geant4SurfaceTest( );
0034
0035       /// Begin-of-event callback
0036       virtual void begin(const G4Event* event);
0037       /// End-of-event callback
0038       virtual void end(const G4Event* event);
0039
0040     protected:
0041
0042     private:
0043
0044     };
0045   }    // End namespace sim
0046 }      // End namespace dd4hep
0047
0048 #endif // DD4HEP_DDG4_GEANT4SURFACETEST_H
0049
0050 // Framework include files
0051 #include <DD4hep/Detector.h>
0052 #include <DD4hep/DD4hepUnits.h>
0053 #include <DD4hep/InstanceCount.h>
0054
0055 #include <DDG4/Factories.h>
0056 #include <DDG4/Geant4Data.h>
0057 #include <DDG4/Geant4Context.h>
0058 #include <DDG4/Geant4SensDetAction.h>
0059 #include <DDG4/Geant4HitCollection.h>
0060
0061 #include <DDRec/Surface.h>
0062 #include <DDRec/DetectorSurfaces.h>
0063 #include <DDRec/SurfaceManager.h>
0064 #include <DDRec/SurfaceHelper.h>
0065
0066 // Geant 4 includes
0067 #include <G4VHitsCollection.hh>
0068 #include <G4HCofThisEvent.hh>
0069 #include <G4Event.hh>
0070
0071 // C/C++ include files
0072 #include <stdexcept>
0073 #include <map>
0074 #include <sstream>
0075
0076 using namespace dd4hep::DDRec;
0077 using namespace dd4hep::detail;
0078 using namespace dd4hep::sim;
0079
0080 DECLARE_GEANT4ACTION(Geant4SurfaceTest)
0081
0082 /// Standard constructor, initializes variables
0083 Geant4SurfaceTest::Geant4SurfaceTest(Geant4Context* ctxt, const std::string& nam)
0084 : Geant4EventAction(ctxt, nam)
0085 {
0086   InstanceCount::increment(this);
0087 }
0088
0089 /// Default Destructor
0090 Geant4SurfaceTest::~Geant4SurfaceTest() {
0091   InstanceCount::decrement(this);
0092 }
0093
0094 /// Begin-of-event callback
0095 void Geant4SurfaceTest::begin(const G4Event* event)  {
0096   if ( event )  {
0097   }
0098 }
0099
0100 /// End-of-event callback
0101 void Geant4SurfaceTest::end(const G4Event* evt)  {
0102   stringstream sst;
0103   Detector& description = context()->detectorDescription();
0104   SurfaceManager& surfMan = *description.extension< SurfaceManager >() ;
0105   const SurfaceMap& surfMap = *surfMan.map( "world" ) ;
0106   G4HCofThisEvent*  hce = evt->GetHCofThisEvent();
0107
0108   if ( !hce )  {
0109     error("Event:%d: No Geant4 event record (HCE) found.",evt->GetEventID());
0110     return;
0111   }
0112
0113   for(int j=0, nCol=hce->GetNumberOfCollections(); j<nCol; ++j)   {
0114     G4VHitsCollection* hc = hce->GetHC(j);
0115     Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(hc);
0116
0117     if ( !(coll && coll->type().type == typeid(Geant4Tracker::Hit)) )   {
0118       continue;
0119     }
0120
0121     SensitiveDetector sd = coll->sensitive()->sensitiveDetector();
0122     if ( !sd.isValid() )  {
0123       error("Event:%d: Collection:%s Invalid sensitive detector.",
0124             evt->GetEventID(),hc->GetName().c_str());
0125       continue;
0126     }
0127     IDDescriptor idDesc    = sd.readout().idSpec();
0128     const string idStr     = idDesc.fieldDescription();
0129     BitField64   idDecoder(idStr);
0130     info("Event: %d [%s: %d hits] Desc:%s",
0131          evt->GetEventID(),hc->GetName().c_str(),coll->GetSize(),idStr.c_str());
0132     for(size_t i=0, nHit=coll->GetSize(); i<nHit; ++i)    {
0133       Geant4Tracker::Hit* h = coll->hit(i);
0134       const Position& pos = h->position;
0135       int trackID = h->truth.trackID;
0136       idDecoder.setValue( h->cellID ) ;
0137       SurfaceMap::const_iterator si = surfMap.find(h->cellID);
0138       ISurface* surf = (si != surfMap.end() ?  si->second  : 0);
0139
0140       // ===== test that we have a surface with the correct ID for every hit ======================
0141       if ( !surf )   {
0142         error("FAILED: Hit:%ld Track:%d No surface found cell id: %016llX",i,trackID,h->cellID);
0143         continue;
0144       }
0145       info("PASSED: Hit:%ld Track:%d Surface found for cell: %016llX",i,trackID,h->cellID);
0146
0147       Vector3D hit_point(pos.x()*dd4hep::mm,pos.y()*dd4hep::mm,pos.z()*dd4hep::mm);
0148       bool isInside = surf->insideBounds(hit_point);
0149       string flag = "";
0150       if ( h->flag & Geant4Tracker::Hit::HIT_SECONDARY_TRACK ) flag += " TRK:SECONDARY";
0151       if ( h->flag & Geant4Tracker::Hit::HIT_KILLED_TRACK    ) flag += " TRK:KILLED";
0152       if ( h->flag & Geant4Tracker::Hit::HIT_NEW_TRACK       ) flag += " NEW_TRACK";
0153       if ( h->flag & Geant4Tracker::Hit::HIT_STARTED_OUTSIDE ) flag += " START:OUTSIDE";
0154       if ( h->flag & Geant4Tracker::Hit::HIT_STARTED_INSIDE  ) flag += " START:INSIDE";
0155       if ( h->flag & Geant4Tracker::Hit::HIT_STARTED_SURFACE ) flag += " START:SURFACE";
0156       if ( h->flag & Geant4Tracker::Hit::HIT_ENDED_OUTSIDE   ) flag += " END:OUTSIDE";
0157       if ( h->flag & Geant4Tracker::Hit::HIT_ENDED_INSIDE    ) flag += " END:INSIDE";
0158       if ( h->flag & Geant4Tracker::Hit::HIT_ENDED_SURFACE   ) flag += " END:SURFACE";
0159
0160       sst.str("");
0161       if ( isInside )  {
0162         sst << hc->GetName() << " " << "PASSED: Track:" << trackID << " Hit:" << i
0163           //<< " Point " << hit_point
0164             << " is ON SURFACE. "
0165             << "(FLAG:" << h->flag << flag << ").";
0166         print(sst.str().c_str());
0167       }
0168       else   {
0169         double dist = surf->distance(hit_point)/dd4hep::mm;
0170         sst.str("");
0171         sst << hc->GetName() << " FAILED: Track:" << trackID << " Hit:" << i
0172             << " Point: " << hit_point;
0173         error(sst.str().c_str());
0174         //sst.str("");
0175         //sst <<   -> Surface:" <<  *surf;
0176         //error(sst.str().c_str());
0177         sst.str("");
0178         sst << "   -> Hit is NOT ON SURFACE (Flag:" << h->flag << flag << ")"
0179             << " Distance to surface:" << dist << " mm.";
0180         error(sst.str().c_str());
0181         // No need for further tests, if this one failed....
0182         continue;
0183       }
0184       // ====== Test to show that moved hit points are outside their surface
0185       hit_point = hit_point + 1e-3 * surf->normal() ;
0186       sst.str("") ;
0187       sst << "Track:" << trackID << " Hit:" << i << " Moved:" << hit_point;
0188       if ( (isInside=surf->insideBounds(hit_point)) )  {
0189         error("FAILED: %s is ON SURFACE (SHOULD NOT BE)",sst.str().c_str());
0190         continue;
0191       }
0192       info("PASSED: %s is NOT ON SURFACE (SHOULD NOT BE)",sst.str().c_str());
0193     }
0194   }
0195 }