diff --git a/scripts/compact_geo_test2.cxx b/scripts/compact_geo_test2.cxx new file mode 100644 index 0000000000000000000000000000000000000000..683be9a14f09557adb8bb078139fb76a66568e77 --- /dev/null +++ b/scripts/compact_geo_test2.cxx @@ -0,0 +1,165 @@ +#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; + +}