diff --git a/README.md b/README.md index 1702448f78465a1e48ca432fdc35a76c96e88bb4..ba4f67e32defcc401f1b4089635ac08e936348e6 100644 --- a/README.md +++ b/README.md @@ -1,379 +1,4 @@ -Tutorial Part 1 -=============== +EIC Software tutorial +===================== - * [Introduction](#introduction) - * [How to build a detector from scratch](#how-to-build-a-detector-from-scratch) - + [Compiling a new detector](#compiling-a-new-detector) - + [Building the geometry](#building-the-geometry) - - [Compact detector description entry element](#compact-detector-description-entry-element) - - [Geometry Construction](#geometry-construction) - * [XML Parsing Tip : Provide good default values](#xml-parsing-tip--provide-good-default-values) - - [Critical parts of build_detector](#critical-parts-of-build_detector) - * [The Readout and Bit Fields](#the-readout-and-bit-fields) - * [Simple Reconstruction Overview of scripts](#simple-reconstruction-overview-of-scripts) - + [Dependencies](#dependencies) - + [Running the scripts](#running-the-scripts) - + [To Do](#to-do) - -Introduction ------------- - -To a newcomer it might not be immediately clear why a *generic detector -library* based on `dd4hep` is a good solution to the 'simulation and -reconstruction geometry problem'. However, with the following examples we -intend to show how generic detector libraries can be used to access the -geometry information in simulation and reconstruction. - -We try to emphasize in the following examples that dd4hep (plus a data model) -allows the software components to be only loosely decoupled. This allows for a -'no framework' framework approach where the emphasis is on the underlying -algorithms used (or being designed). Furthermore, using ROOT's TDataframe, the -details of each task are made manifest. - -Furthermore, the 'no framework' framework allows for the each step to be -executed using any tool available, however, the IO and data model should be -fixed. - - -How to build a detector from scratch ------------------------------------- - -Building a new (generic) detector using dd4hep is rather straight forward if we -make the following assumptions. - -1. We will use the built-in sensitive detectors -2. We will use the built-in data model (hits) associated with these detectors - -These items can be customized by using the DD4hep plugin mechanism. This will -be covered in another tutorial. - -### Compiling a new detector - -For this tutorial we will build a simplified Roman Pot detector. -We will discuss the detector built in the source file -`src/GenericDetectors/src/SimpleRomanPot_geo.cpp`. -To compile this detector into the GenericDetectors library the detector needs -to be added to the list of sources in the cmake file -`src/GenericDetectors/CMakeLists.txt`. - -``` -dd4hep_add_plugin(${a_lib_name} SOURCES - src/BeamPipe_geo.cpp - ... - src/SimpleRomanPot_geo.cpp # add this line - ) -``` - -### Building the geometry - -The work of defining the detector is done in a function (here called -`build_detector`) that is registered using the DD4hep plugin macro -`DECLARE_DETELEMENT`. - -``` -static Ref_t build_detector(Detector& dtor, xml_h e, SensitiveDetector sens) -{ - xml_det_t x_det = e; - Material air = dtor.air(); -... -} -DECLARE_DETELEMENT(SimpleRomanPot, build_detector) -``` - -The argument signature of the `build_detector` is: -- `Detector& dtor`: This handle provides the main hook to the detector tree (`dd4hep::Detector`). -- `xml_h e`: Handle to the XML `<detector>` tag associated with the detector in - the "compact" detector description file (more on this later). This provides - the mechanism for feeding in the run-time construction parameters. -- `SensitiveDetector sens`: The sensitive detector to be assigned to the - sensitive volumes/elements of the detector. - -The DD4hep plugin macro `DECLARE_DETELEMENT(SimpleRomanPot, build_detector)` -stamps out the necessary boiler plate code to register a new detector called -`SimpleRomanPot` which is build by calling `build_detector`. - -#### Compact detector description entry element - -The `<detector>` tag defines a new instance of a detector and requires the -attributes "id", "name", and "type". For example: - -``` -<detector id="1" name="MyRomanPot" type="SimpleRomanPot" - vis="RedVis" readout="RomanPotHits" zoffset="1.0*m"> -</detector> -```` - -This defines an instance of the detector named "MyRomanPot" of type -"SimpleRomanPot" (i.e. the type-name given in the first argument of the DD4hep -`DECLARE_DETELEMENT` macro) and with id=1. The additional attributes (vis, -readout, zoffset) will are optional. - -The detector tag is provided as the second argument in the `build_detector` -function. It can be parsed *how ever you want* in order to extract detector -information. The allowed attributes are listed in -[`UnicodeValues.h`](http://test-dd4hep.web.cern.ch/test-dd4hep/doxygen/html/_unicode_values_8h_source.html) -where it is clear how to add new attributes. - - -#### Geometry Construction - -If you are familiar with Geant4 or TGeo geometry construction then this will be -easy. DD4hep has TGeo under hood but there are a few naming tweaks. The -following table will help orient the user. - -| *DD4hep* | *Geant4* | *TGeo | -| ------ | -------- | ---- | -| Solid | G4VSolid | TGeoShape | -| Volume | G4LogicalVolume | TGeoVolume | -| PlacedVolume | G4PVPlacement | TGeoNode | -| Element | G4Element | TGeoElement | -| Material | G4Material | TGeoMaterial/TGeoMedium| - -##### XML Parsing Tip : Provide good default values - -If you have a detector parameter which we later will tweak (while optimizing -the design) try to get the value from the xml element but provide a good -default value. For example: - -``` -double radius = ( x_det.hasAttr(_Unicode(radius)) ) ? x_det.attr<double>(_Unicode(radius)) : 5.0*dd4hep::cm; -``` - -This provides a default radius of 5 cm when `x_det` does not have a "radius" -attribute defined. We will return to this later. - -#### Critical parts of build_detector - - -We will now look at parts of the source file `src/GenericDetectors/src/SimpleRomanPot_geo.cpp`. - -``` -static Ref_t build_detector(Detector& dtor, xml_h e, SensitiveDetector sens) -{ - xml_det_t x_det = e; - Material air = dtor.air(); - Material carbon = dtor.material("CarbonFiber"); - Material silicon = dtor.material("SiliconOxide"); - Material aluminum = dtor.material("Aluminum"); - Material vacuum = dtor.material("Vacuum"); - Material supp_mat = carbon; - Material sens_mat = silicon; - int det_id = x_det.id(); // id=1 - string det_name = x_det.nameStr(); // "MyRomanPot" -``` - -Here we are grabbing the materials that are assumed to be already defined. Also -we are getting the detector id and name defined in the `detector` tag. - -Next we define an [Assembly -volume](http://test-dd4hep.web.cern.ch/test-dd4hep/doxygen/html/classdd4hep_1_1_assembly.html). -Here we also stumble upon the important class -[`dd4hep::DetElement`](http://test-dd4hep.web.cern.ch/test-dd4hep/doxygen/html/classdd4hep_1_1_det_element.html#details). -It is a means of providing the detector hierarchy/tree, but doesn't necessarily -have to map exactly to detector geometry. However, it typically will typically -parallel the geometry (and probably should). - -``` - string module_name = "RomanPot"; - Assembly assembly(det_name + "_assembly"); - DetElement sdet( det_name, det_id); - sens.setType("tracker"); -``` -The last line sets the `SensitiveDetector sens` argument to be the tracker type -(which is a DD4hep built-in). We will soon assign this to sensitive volumes. -`sdet` is associated with the mother detector element by the constructor which -looks up the detector name (here "MyRomanPot"). - -``` - double z_offset = (x_det.hasAttr(_Unicode(zoffset))) ? x_det.attr<double>(_Unicode(zoffset)) : 0.0; - double thickness = (x_det.hasAttr(_Unicode(thickness))) ? x_det.attr<double>(_Unicode(thickness)) : 0.01*dd4hep::cm; -``` -Here we grab attributes and provide default values. We continue with default -values that could also be define through attributes, however, we will want to -add child elements of the detector tag (so the attributes does not grow too -long). - -``` - double rp_chamber_thickness = 5.0*dd4hep::mm; - double rp_chamber_radius = 5.0*dd4hep::cm; - double rp_chamber_length = 50.0*dd4hep::cm; - Tube rp_beam_pipe_tube(rp_chamber_radius, rp_chamber_radius+rp_chamber_thickness, rp_chamber_length/2.0); - Tube rp_beam_vacuum_tube(0.0, rp_chamber_radius+rp_chamber_thickness, rp_chamber_length/2.0); - Tube rp_beam_vacuum_tube2(0.0, rp_chamber_radius, rp_chamber_length/2.0); - - double rp_detector_tube_radius = 2.5*dd4hep::cm; - double rp_detector_tube_length = 20.0*dd4hep::cm; - Tube rp_detector_tube(rp_detector_tube_radius, rp_detector_tube_radius+rp_chamber_thickness, rp_detector_tube_length/2.0); - Tube rp_detector_vacuum_tube(0.0, rp_detector_tube_radius+rp_chamber_thickness, rp_detector_tube_length/2.0); - Tube rp_detector_vacuum_tube2(0.0, rp_detector_tube_radius, rp_detector_tube_length/2.0); - - ROOT::Math::Rotation3D rot_X( ROOT::Math::RotationX(M_PI/2.0) ); - ROOT::Math::Rotation3D rot_Y( ROOT::Math::RotationY(M_PI/2.0) ); - - UnionSolid rp_chamber_tee1(rp_beam_vacuum_tube, rp_detector_vacuum_tube, rot_X); - UnionSolid rp_chamber_tee12(rp_chamber_tee1, rp_detector_vacuum_tube, rot_Y); - - UnionSolid rp_chamber_tee2(rp_beam_vacuum_tube2, rp_detector_vacuum_tube2, rot_X); - UnionSolid rp_chamber_tee22(rp_chamber_tee2, rp_detector_vacuum_tube2, rot_Y); - - SubtractionSolid sub1(rp_chamber_tee12,rp_chamber_tee22); - Volume rp_chamber_vol("rp_chamber_walls_vol", sub1, aluminum); - Volume rp_vacuum_vol("rp_chamber_vacuum_vol", rp_chamber_tee22, vacuum); -``` -The above code builds the up the two solids associated with 3 tubes -intersecting. One volume is the aluminum vacuum chamber walls and the other is -the vacuum contained within this envelope. - -Next we must place these two volumes in the assembly volume (which is just an -empty container-like volume. The PlacedVolume is then associated with a -BitFieldValue in the readout's BitField64 readout string. In this case the -"layer" BitFieldValue. The BitField64 is used to construct unique VolumeIDs and -CellIDs for PlacedVolumes and Segmentations respectively. - -``` - PlacedVolume pv; - pv = assembly.placeVolume( rp_chamber_vol ); - pv = assembly.placeVolume( rp_vacuum_vol ); - pv.addPhysVolID( "layer", 2 ); -``` - -Set the PlacedVolume BitFieldValue ID. "2" in this case. - -``` - double supp_x_half = 1.0*dd4hep::cm; - double supp_y_half = 1.0*dd4hep::cm; - double supp_thickness = 1.0*dd4hep::mm; - double supp_gap_half_width = 1.0*dd4hep::mm; - double sens_thickness = 0.1*dd4hep::mm; - - Box supp_box( supp_x_half, supp_y_half, supp_thickness/2.0 ); - Box sens_box( supp_x_half-supp_gap_half_width, supp_y_half-supp_gap_half_width, sens_thickness/2.0 ); - -``` - -Next we define vectors which are used to define a "surface" (which will later -generate simulation tracker hits). - -``` - // create a measurement plane for the tracking surface attched to the sensitive volume - Vector3D u( 1. , 0. , 0. ) ; - Vector3D v( 0. , 1. , 0. ) ; - Vector3D n( 0. , 0. , 1. ) ; - //Vector3D o( 0. , 0. , 0. ) ; - - // compute the inner and outer thicknesses that need to be assigned to the tracking surface - // depending on wether the support is above or below the sensor - // The tracking surface is used in reconstruction. It provides material thickness - // and radation lengths needed for various algorithms, routines, etc. - double inner_thickness = supp_thickness/2.0; - double outer_thickness = supp_thickness/2.0; - double z_shift = 5.0*dd4hep::mm; - double xy_shift = 15.0*dd4hep::mm; - - SurfaceType type( SurfaceType::Sensitive ) ; -``` - -We now define a simple rectangular pixel sensor. This will be the first of -four: two will come in along the x axis and two along the y axis. - -``` - // ------------- x1 - Volume support1_vol( "xsenseor_supp", supp_box, supp_mat ); - Volume sensor1_vol( "xsenseor_sens", sens_box, sens_mat ); - VolPlane surf1( sensor1_vol, type, inner_thickness , outer_thickness, u,v,n); - sensor1_vol.setSensitiveDetector(sens); -``` -The code above builds two volumes one which will contain the sensitive volume. -The sensitive volume is assigned to be a sensitive detector. - -``` - DetElement layer1_DE( sdet, "layer1_DE", 1 ); - pv = rp_vacuum_vol.placeVolume( support1_vol, Position(xy_shift,0, -z_shift) ); - pv.addPhysVolID("layer", 1 ); - layer1_DE.setPlacement( pv ) ; -``` -The above creates a new DetElement, places the support volume in the vacuum, -sets this PlacedVolume to associate with the layer bitfield, and adds the PV to -the DetElement. Note the DetElement is constructed with the parent element -(sdet) being the first argument. In this way it is clear we are building, -(semi-)parallel to the geometry, a detector element hierarchy. - -``` - DetElement mod1( layer1_DE , "module_1", 1 ); - pv = support1_vol.placeVolume(sensor1_vol, Position(0,0,0)); - pv.addPhysVolID("module", 1 ); - mod1.setPlacement( pv ); -``` -This creates the module DetElement and places the sensitive volume in the -support volume (already placed). It also assocates the pv with the module -bitfield and then adds the PV to the DetElement. - -The above is repeated for the three remaining sensitive elements. - -Finally we get the top level volume to place the assemble volume. Note we are -using the zoffset. This PV is then associated with the top level "system" -bitfieldvalue. - -``` - pv = dtor.pickMotherVolume(sdet).placeVolume(assembly, Position(0,0,z_offset)); - pv.addPhysVolID("system", det_id); // Set the subdetector system ID. - sdet.setPlacement(pv); - - assembly->GetShape()->ComputeBBox() ; - return sdet; -} - -``` - -The Readout and Bit Fields --------------------------- - - -Simple Reconstruction Overview of scripts ---------------------------------------- - -In the examples here we use ROOT and the LCIO data model. -The following scripts should be executed in order: - -1. `run_example` : bash script that runs `ddsim` with the Geant4 macro file - `gps.mac` and the input geometry defined in `gem_tracker_disc.xml`. -2. `scripts/example_hit_position.cxx` : ROOT script showing how to use both the - "segmentation" and "surface" geometry to get the hit position information. -3. `scripts/example_cell_size.cxx` : ROOT script showing how information about - the geometry beyond the position. -4. `scripts/example_digi.cxx` : ROOT script performing a crude digitization of - gem tracker hits. -5. `scripts/example_hit_recon.cxx` : ROOT script that shows how to to use - dd4hep to lookup a 'segment' (cell) position based on the unique cellID - attached to the hit object. -6. `scripts/example_tracking.cxx` : ROOT script showing how to use the - `GenFind` finding library along with `GenFit` to produce tracks. - -### Dependencies - -There are some library dependencies: - -- DD4hep -- lcgeo -- GenFind -- GenFit -- NPdet - - -### Running the scripts - - -``` -./run_example -root scripts/example_digi.cxx++ -root scripts/example_hit_position.cxx++ # no output -root scripts/example_cell_size.cxx++ # no output -root scripts/example_hit_recon.cxx++ -root scripts/example_tracking.cxx++ -``` - -### To Do - -- Show how to build a detector -- Finish track fitting script. -- Create example for basic PFA usage +https://eicweb.phy.anl.gov/Argonne_EIC/tutorial/tutorial_part1 diff --git a/gem_tracker.xml b/gem_tracker.xml index cfc4938d9a58c953160aa927729c9c58ea354214..5fd593a13a3d1bb7e3f348a94a2cb1c03cb66238 100644 --- a/gem_tracker.xml +++ b/gem_tracker.xml @@ -127,7 +127,7 @@ <layer id="5" z="315 *cm" inner_r="109.0*cm" outer_r="237.0*cm" phi0_offset="-0.5*deg" /> </detector> --> - <detector id="2" name="GEMTracker_SIDIS" vis="RedVis" type="GEMTrackerDisc" readout="GEMTrackerHits" > + <detector id="2" name="GEMTracker_SIDIS" vis="RedVis" type="GEMTrackerDiscSOLID" readout="GEMTrackerHits" > <layer id="1" z="-175 *cm" inner_r="36*cm" outer_r="87.0*cm" phi0_offset="0.0*deg" /> <layer id="2" z="-150 *cm" inner_r="21*cm" outer_r="98.0*cm" phi0_offset="0.0*deg" /> <layer id="3" z="-119 *cm" inner_r="25*cm" outer_r="112.0*cm" phi0_offset="0.0*deg" /> @@ -183,7 +183,7 @@ <fields> <field name="GlobalSolenoid" type="solenoid" - inner_field="10.0*tesla" + inner_field="5.0*tesla" outer_field="-0.5*tesla" zmax="4.0*m" outer_radius="2.0*m"> diff --git a/gps.mac b/gps.mac index 6f6e0ceb43187c31691637dc05ba580b6c36aa16..34a060adb34323d971800a0fdb6573c1dd113540 100644 --- a/gps.mac +++ b/gps.mac @@ -2,11 +2,11 @@ /run/initialize /gps/verbose 2 -/gps/particle e- +/gps/particle pi- /gps/number 1 /gps/ene/type Gauss -/gps/ene/mono 9.0 GeV +/gps/ene/mono 3.5 GeV /gps/ene/sigma 3.0 GeV /gps/pos/type Volume @@ -14,9 +14,11 @@ /gps/pos/centre 0.0 0.0 0.0 cm /gps/pos/radius 0.01 cm /gps/pos/halfz 0.01 cm -/gps/position 0 0 -6.0 m +/gps/position 0 0 0 cm -/gps/direction 0 0.1 1.0 -#/gps/ang/type iso +#/gps/direction 0 0.1 1.0 +/gps/ang/type iso +/gps/ang/mintheta 30 degree +/gps/ang/maxtheta 160 degree /run/beamOn 100 diff --git a/scripts/example_cell_size.cxx b/scripts/example_cell_size.cxx index 2bdc1b2f8a09466a93ab86190b64755a85516a1b..40020a3c3475af7253219bf7bc8f4eb10a285d65 100644 --- a/scripts/example_cell_size.cxx +++ b/scripts/example_cell_size.cxx @@ -1,15 +1,7 @@ -// -//#include <vector> -//#include <tuple> -//#include <algorithm> -//#include <iterator> -// -//// DD4hep -//// ----- -//// In .rootlogon.C -//// gSystem->Load("libDDDetectors"); -//// gSystem->Load("libDDG4IO"); -//// gInterpreter->AddIncludePath("/opt/software/local/include"); +R__LOAD_LIBRARY(libfmt.so) +#include "fmt/core.h" +R__LOAD_LIBRARY(libDDG4IO.so) + #include "DD4hep/Detector.h" #include "DDG4/Geant4Data.h" #include "DDRec/CellIDPositionConverter.h" @@ -17,23 +9,10 @@ #include "DDRec/Surface.h" #include "ROOT/RDataFrame.hxx" -//#include "lcio2/MCParticleData.h" -//#include "lcio2/ReconstructedParticleData.h" - -//#include "Math/Vector3D.h" -//#include "Math/Vector4D.h" -//#include "Math/VectorUtil.h" #include "TCanvas.h" -//#include "TLegend.h" -//#include "TMath.h" -//#include "TRandom3.h" -//#include "TFile.h" -//#include "TH1F.h" -//#include "TH1D.h" -//#include "TTree.h" #include "TChain.h" -//#include "TF1.h" -/** Cell size example. + + /** Cell size example. * * Makes use of ROOT's TDataframe * (https://root.cern.ch/doc/master/classROOT_1_1Experimental_1_1TDataFrame.html) @@ -54,7 +33,7 @@ * This example shows how to get both kinds of hit positions. * */ -void example_cell_size(const char* fname = "test_tracker_disc.root"){ +void example_cell_size(const char* fname = "gem_tracker_sim.root"){ using namespace ROOT::Math; @@ -65,14 +44,14 @@ void example_cell_size(const char* fname = "test_tracker_disc.root"){ ROOT::RDataFrame d0(*t, {"GEMTrackerHits","MCParticles"}); //How to get the type of the initial branch: (vector<dd4hep::sim::Geant4Tracker::Hit*>) - //std::cout << t->GetBranch("GEMTrackerHits")->GetClassName() << std::endl; + std::cout << t->GetBranch("MCParticles")->GetClassName() << std::endl; // ------------------------- // Get the DD4hep instance // Load the compact XML file // Initialize the position converter tool dd4hep::Detector& detector = dd4hep::Detector::getInstance(); - detector.fromCompact("gem_tracker_disc.xml"); + detector.fromCompact("gem_tracker.xml"); dd4hep::rec::CellIDPositionConverter cellid_converter(detector); // ------------------------- @@ -132,10 +111,12 @@ void example_cell_size(const char* fname = "test_tracker_disc.root"){ dd4hep::rec::ISurface* surf = (si != surfMap->end() ? si->second : nullptr); dd4hep::rec::Vector3D pos2; // Get the origin of the surface volume - if(!surf) pos2 = surf->origin(); std::cout << " Hit Position : " << pos0 << std::endl; std::cout << "Segmentation-Cell Position : " << pos1 << std::endl; - std::cout << " Surface Origin Position : " << pos2 << std::endl; + if(surf) { + pos2 = surf->origin(); + std::cout << " Surface Origin Position : " << pos2 << std::endl; + } result.push_back( VectorUtil::RotateZ((pos1-pos0), -1.0*pos1.Phi()+TMath::Pi()/2.0).x() ); } return result; @@ -155,9 +136,7 @@ void example_cell_size(const char* fname = "test_tracker_disc.root"){ dd4hep::rec::Vector3D pos2; // Get the origin of the surface volume if(!surf) pos2 = surf->origin(); - //std::cout << " Hit Position : " << pos0 << std::endl; - //std::cout << "Segmentation-Cell Position : " << pos1 << std::endl; - //std::cout << " Surface Origin Position : " << pos2 << std::endl; + //rotate to cell coordinates result.push_back( VectorUtil::RotateZ((pos1-pos0), -1.0*pos1.Phi()+TMath::Pi()/2.0).y() ); } return result; @@ -175,22 +154,11 @@ void example_cell_size(const char* fname = "test_tracker_disc.root"){ return false;}, {"GEMTrackerHits"}) .Define("deltax", deltax, {"GEMTrackerHits"}) .Define("deltay", deltay, {"GEMTrackerHits"}); -// .Define("derp", -// [](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits){ -// std::vector<double> res; -// std::transform(hits.begin(), -// hits.end(), -// std::back_inserter(res), -// [&](auto h) { return h->position.x()/100.0 - 1.0; } ); -// return res; -// }, {"GEMTrackerHits"}) - - - //auto h0 = d1.Histo1D(TH1D("h0", "nhits; ", 100, -10.0,10.0), "deltay"); + auto h0 = d1.Histo2D( TH2D("h0", "nhits; ", 100, -2.0, 2.0, 100, -2.0, 2.0), "deltax", "deltay"); TCanvas* c = new TCanvas(); h0->DrawClone("colz"); + c->SaveAs("scripts/example_cell_size.png"); - //d1.Snapshot("EVENT","derp.root"); } diff --git a/scripts/example_digi.cxx b/scripts/example_digi.cxx index d5f8284783000003e6f53a17b30477a8dca5c1bb..4307be34b7bf94975ba0966ef53e84fd031a2e31 100644 --- a/scripts/example_digi.cxx +++ b/scripts/example_digi.cxx @@ -96,10 +96,11 @@ void example_digi(const char* fname = "gem_tracker_sim.root"){ .Define("RawTrackerHits", digitize_gem_hits, {"GEMTrackerHits"}) ; - auto h0 = d1.Histo1D(TH1D("h0", "nhits; ", 20, 0,20), "nhits"); + auto h0 = d1.Histo1D(TH1D("h0", "Number of Hits in GEM Tracker; N hits ", 20, 0,20), "nhits"); TCanvas* c = new TCanvas(); d1.Snapshot("digitized_EVENT","gem_tracker_digi.root"); h0->DrawClone(); + c->SaveAs("scripts/example_digi.png"); }