Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#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_test2(
const char* DIRNAME = "data/tests",
const char* result_file_prefix = "tests",
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);
using namespace DD4hep;
DD4hep::Geometry::DetElement det = lcdd.detector( "SiVertexBarrel" ) ;
//DD4hep::DDRec::SurfaceManager& surfMan = *lcdd.extension<DD4hep::DDRec::SurfaceManager>() ;
//_map = surfMan.map( det.name() ) ;
//if( ! _map ) {
// std::stringstream err ; err << " Could not find surface map for detector: "
// << _subDetName << " in SurfaceManager " ;
// std::cout << err << std::endl;;
//}
auto readout = lcdd.readout("SiTrackerBarrelHits");
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
auto det_el = lcdd.detector("SiVertexBarrel");
auto vol_man = lcdd.volumeManager();
std::cout << " detector type : " << det_el.type() << std::endl;
std::cout << " detector key : " << det_el.key() << std::endl;
std::cout << " detector level : " << det_el.level() << std::endl;
std::cout << " detector id : " << det_el.id() << std::endl;
std::cout << " detector ppath : " << det_el.placementPath() << std::endl;
std::cout << " children : " << det_el.children().size() << std::endl;
auto children = det_el.children();
for(auto& child : children) {
//std::cout << " child : " << child.first << std::endl;
//std::cout << " child type : " << child.second.type() << std::endl;
//std::cout << " child key : " << child.second.key() << std::endl;
//std::cout << " child level : " << child.second.level() << std::endl;
//std::cout << " child id : " << child.second.id() << std::endl;
std::cout << " grandchildren: " << child.second.children().size() << std::endl;
for(auto& gc : child.second.children()) {
//std::cout << " grandchild : " << gc.first << std::endl;
//std::cout << " grandchild name : " << gc.second.name() << std::endl;
//std::cout << " grandchild type : " << gc.second.type() << std::endl;
//std::cout << " grandchild key : " << gc.second.key() << std::endl;
//std::cout << " grandchild level : " << gc.second.level() << std::endl;
std::cout << " grandgrandchildren: " << gc.second.children().size() << std::endl;
for(auto& g2c : gc.second.children()) {
std::cout << " g2child : " << g2c.first << std::endl;
std::cout << " g2child name : " << g2c.second.name() << std::endl;
std::cout << " g2child type : " << g2c.second.type() << std::endl;
std::cout << " g2child key : " << g2c.second.key() << std::endl;
std::cout << " g2child level : " << g2c.second.level() << std::endl;
std::cout << " g3children : " << g2c.second.children().size() << std::endl;
}
}
}
auto seg = readout.segmentation();
std::cout << " num col " << readout.numCollections() << std::endl;
auto id = readout.idSpec();
//std::cout << "id encode " << id.encode({1,0,1,3,1,0,0}) << std::endl;
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);
std::vector<std::string> files = {"data/tests/slic_ddsim_out.slcio"};
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());
//auto id_decoder = DD4hep::DDRec::IDDecoder::getInstance()
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* simtrackhits_col = (IMPL::LCCollectionVec*) evt->getCollection("SiVertexBarrelHits") ;
UTIL::LCTOOLS::printSimTrackerHits( simtrackhits_col );
//UTIL::LCTOOLS::printMCParticles( col ) ;
int nMCP = thrown_col->getNumberOfElements() ;
int nSimTrackHits = simtrackhits_col->getNumberOfElements() ;
CellIDDecoder<SimTrackerHit> lcio_id(simtrackhits_col);
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();
auto cell_id = ahit->getCellID0();
//std::cout << "id " << ahit->getCellID0() << 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()),
// 1);
// std::cout << " cellid " << cellid << std::endl;
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
//if(col->getParameters().getStringVal(LCIO::CellIDEncoding) != ""){
std::cout << " id-fields: ("<< lcio_id(const_cast<EVENT::SimTrackerHit*>(ahit)).valueString() << ")" << std::endl ;
//auto det_el_mod = det_el.children(lcio_bf["layer" ].value()
//auto lcio_bf = lcio_id(const_cast<EVENT::SimTrackerHit*>(ahit));
//auto bf = id.decoder();
//bf->setValue( lcio_bf.getValue() );
//DD4hep::Geometry::IDDescriptor an_id(lcio_id(const_cast<EVENT::SimTrackerHit*>(ahit)).valueString());
//
//auto bf2 = (*bf) ;
//(bf2)["system"] = lcio_bf["system"].value();
//(bf2)["barrel"] = lcio_bf["barrel"].value();
//(bf2)["layer" ] = lcio_bf["layer" ].value();
//(bf2)["module"] = lcio_bf["module"].value();
//(bf2)["sensor"] = 0;//lcio_bf["sensor"].value();
//(bf2)["side" ] = 0;//lcio_bf["side" ].value();
//(bf2)["strip" ] = 0;//lcio_bf["strip" ].value();
//std::cout << " DD4hep bf value " << bf2.getValue() << std::endl;
//std::cout << " Lcio bf value " << lcio_bf.getValue() << std::endl;
//auto pv = vol_man.lookupPlacement(bf2.getValue());
//std::cout << pv.toString() << std::endl;
//std::vector<double> v = {0,0,0};
//pv->LocalToMaster(pv->GetMatrix()->GetTranslation(),&v[0] );
//TVector3 aPos(&v[0] );
//aPos.Print();
////pv.name()
//std::cout << "system" << lcio_bf["system"].value() << std::endl;
//std::cout << "barrel" << lcio_bf["barrel"].value() << std::endl;
//std::cout << "layer " << lcio_bf["layer" ].value() << std::endl;
//std::cout << "module" << lcio_bf["module"].value() << std::endl;
//std::cout << "sensor" << lcio_bf["sensor"].value() << std::endl;
//std::cout << "side " << lcio_bf["side" ].value() << std::endl;
//std::cout << "strip " << lcio_bf["strip" ].value() << std::endl;
//seg.setDecoder( bf );
//auto cell_id2 = lcio_bf.getValue();
//auto pos = seg.position(cell_id2);
//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());
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
//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;
}