new file: scripts/compact_geo_test2.cxx

scripts/compact_geo_test2.cxx

+#include "UTIL/LCTOOLS.h"
+#include "UTIL/Operators.h"
+#include <UTIL/CellIDDecoder.h>
+//#include "DD4hep/LCDD.h"
+//#include "DD4hep/DD4hepUnits.h"
+//#include "DD4hep/Objects.h"
+// check directory
+#include <sys/stat.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include "TMath.h"
+#include"time.h"
+#include <dirent.h>
+#include <string>
+#include <vector>
+#include <map>
+struct stat sb;
+
+std::vector<std::string> get_files_in_directory(std::string path = "data/.") {
+  std::vector<std::string> files;
+  if (stat(path.c_str(), &sb) == 0 && S_ISDIR(sb.st_mode)) { 
+    DIR*    dir;
+    dirent* pdir;
+    dir = opendir(path.c_str());
+    while( pdir = readdir(dir) ){
+        string rfile=path + "/" + pdir->d_name;
+        if (rfile.find(".slcio") != std::string::npos) 
+             files.push_back(rfile); 
+    }
+    } else {
+        cout << "Directory " << path << " does not exist!!" << endl;
+        std::exit(0); 
+        }
+    return files;
+}
+//______________________________________________________________________________
+
+
+void compact_geo_test(
+    const char* DIRNAME            = "data/electrons",
+    const char* result_file_prefix = "SID_1T",
+    const char* compact_file       = "detectors/SiD/slic/sieic3/sieic3_compact.xml")
+    //const char* compact_file       = "detectors/SiD/compact/sid_working/sidloi3_v00.xml")
+{
+
+  using namespace EVENT ;
+  using namespace UTIL;
+  using namespace IMPL ;
+
+  DD4hep::Geometry::LCDD& lcdd = DD4hep::Geometry::LCDD::getInstance();
+  lcdd.fromCompact(compact_file);
+  const double position[3]={0,0,0}; // position to calculate magnetic field at (the origin in this case)
+  double bField[3]={0,0,0}; 
+  lcdd.field().magneticField(position,bField); 
+  _Bz = bField[2]/dd4hep::tesla; 
+
+  std::cout << " Magnetic Field Bz = " << _Bz << std::endl;
+
+  TH1F * hR0   = new TH1F("hR0","R;R [mm];",100,0,1000);
+  TH1F * hZ0   = new TH1F("hZ0","Z;Z [mm];",100,-1000,1000);
+  TH2F * hZR0   = new TH2F("hZR0","ZR;Z [mm];R [mm];",100,-1000,1000,100,0,1000);
+
+  auto readout = lcdd.readout("SiTrackerBarrelHits");
+  auto seg     = readout.segmentation();
+  std::cout << " num col " << readout.numCollections() << std::endl;
+  auto id      = readout.idSpec();
+  std::cout << " id spec " << id.toString() << std::endl;
+  std::cout << " field " << id.fieldDescription() << std::endl;
+  std::cout << " decoder " << id.decoder() << std::endl;
+  auto bf = id.decoder();
+  std::cout << " barrel "   << bf->index("barrel")<<std::endl;
+  std::cout << " barrel "   << (*bf)["barrel"].name() <<std::endl;
+  std::cout << " barrel "   << (*bf)["barrel"].offset() <<std::endl;
+  std::cout << " barrel "   << (*bf)["barrel"].width() <<std::endl;
+  std::cout << " barrel "   << (*bf)["barrel"].value() <<std::endl;
+
+  int  max_files    = 2;
+  int  nEvents      = 0  ;
+  int  n_files_read = 0;
+  auto files        = get_files_in_directory(DIRNAME);
+
+  IO::LCReader* lcReader = IOIMPL::LCFactory::getInstance()->createLCReader() ;
+  EVENT::LCEvent* evt = 0 ;
+  std::string mcpType("std::vector<IMPL::MCParticleImpl*>") ;
+  std::string mcpName("MCParticle") ;
+
+  std::sort(files.begin(),files.end());
+
+  for(auto aFile : files) {
+
+    std::cout << "Opening " << aFile << '\n';
+    lcReader->open( aFile.c_str() ) ;
+    n_files_read++;
+
+    //----------- the event loop -----------
+    while( (evt = lcReader->readNextEvent()) != 0 ) {
+
+      if(nEvents%1000==0) {
+        if(nEvents==0)UTIL::LCTOOLS::dumpEvent( evt ) ;
+        std::cout << "Event " << nEvents << '\n';
+      }
+
+      IMPL::LCCollectionVec* thrown_col = (IMPL::LCCollectionVec*) evt->getCollection( "MCParticle"  ) ;
+      IMPL::LCCollectionVec* tracks_col = (IMPL::LCCollectionVec*) evt->getCollection( "Tracks"  ) ;
+      IMPL::LCCollectionVec* simtrackhits_col = (IMPL::LCCollectionVec*) evt->getCollection("SiTrackerBarrelHits") ;
+
+      //UTIL::LCTOOLS::printMCParticles( col ) ;
+
+      int nMCP    = thrown_col->getNumberOfElements() ;
+      int nTracks = tracks_col->getNumberOfElements() ;
+      bool hasOneTrack  = (nTracks==1);
+      bool hasTwoTracks = (nTracks==2);
+      int nSimTrackHits = simtrackhits_col->getNumberOfElements() ;
+
+      for(int i_track=0 ; i_track<nSimTrackHits ; ++i_track) {
+        auto  ahit = (EVENT::SimTrackerHit*) simtrackhits_col->getElementAt(i_track);
+        double hit_time  = ahit->getTime();
+        double EDep  = ahit->getEDep();
+        //std::cout << "id " << ahit->id() << std::endl;
+        //if(hit_time != 0.0) {
+        //  std::cout << " time " << hit_time << '\n';
+        //}
+        TVector3 x0(ahit->getPosition());
+        TVector3 x1(ahit->getPosition());
+        x1.Print();
+        auto cellid     = seg.cellID(
+            DD4hep::Geometry::Position(x0.x(),x0.y(),x0.z()),
+            DD4hep::Geometry::Position(x1.x(),x1.y(),x1.z()),
+            5);
+        // std::cout << " cellid " << cellid << std::endl;
+        auto pos = seg.position(cellid);
+        //std::cout << " x=" << pos.x() << ", y=" << pos.y() << ", z=" << pos.z() << std::endl;
+        //std::cout << " R   = " << pos.R() << std::endl;
+        //std::cout << " rho = " << pos.Rho() << std::endl;
+        hR0->Fill(pos.R());
+        hZR0->Fill(pos.Z(),pos.R());
+
+        //if(nTracks>1) {
+        //  std::cout <<  (*track) << std::endl;
+        //}
+      }
+      nEvents++;
+    }
+    // -------- end of event loop -----------
+
+    lcReader->close() ;
+
+    if(n_files_read >= max_files) break;
+  }
+
+      //auto bf      = seg.decoder();
+      //std::cout << bf->valueString() << std::endl;
+
+      //auto pos     = seg.position(1000);
+      //std::cout << " x=" << pos.x() << ", y=" << pos.y() << ", z=" << pos.z() << std::endl;
+
+  TCanvas * c = new TCanvas();
+  hR0->Draw();
+
+  c = new TCanvas();
+  hZR0->Draw("cols");
+  return;
+
+}