diff --git a/examples/TOC b/examples/TOC new file mode 100644 index 0000000000000000000000000000000000000000..74cee7020a0eeac9f75f180b4666af46aa1cd482 --- /dev/null +++ b/examples/TOC @@ -0,0 +1,14 @@ +- [NPdet Examples](#npdet-examples) + * [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) \ No newline at end of file diff --git a/src/ConceptDetectors/topside/.gitignore b/src/ConceptDetectors/topside/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..74ed96fff717bf86d816a5e0ec581b70dca70666 --- /dev/null +++ b/src/ConceptDetectors/topside/.gitignore @@ -0,0 +1 @@ +!*.png diff --git a/src/GenericDetectors/calorimeters/src/PolyhedraEndcapCalorimeter2_geo.cpp b/src/GenericDetectors/calorimeters/src/PolyhedraEndcapCalorimeter2_geo.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f9782d49425513f086a63fa70ea9c104e3328a2e --- /dev/null +++ b/src/GenericDetectors/calorimeters/src/PolyhedraEndcapCalorimeter2_geo.cpp @@ -0,0 +1,134 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== +// +// Specialized generic detector constructor +// +//========================================================================== +#include "DD4hep/DetFactoryHelper.h" +#include "XML/Layering.h" + +using namespace std; +using namespace dd4hep; +using namespace dd4hep::detail; + +static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) { + xml_det_t x_det = e; + xml_dim_t dim = x_det.dimensions(); + int det_id = x_det.id(); + bool reflect = x_det.reflect(true); + string det_name = x_det.nameStr(); + Material air = description.air(); + int numsides = dim.numsides(); + double rmin = dim.rmin(); + double rmax = dim.rmax()*std::cos(M_PI/numsides); + double zmin = dim.zmin(); + Layering layering(x_det); + double totalThickness = layering.totalThickness(); + Volume endcapVol("endcap",PolyhedraRegular(numsides,rmin,rmax,totalThickness),air); + DetElement endcap("endcap",det_id); + + int l_num = 1; + int layerType = 0; + double layerZ = -totalThickness/2; + + endcapVol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),x_det.visStr()); + + for(xml_coll_t xc(x_det,_U(layer)); xc; ++xc) { + xml_comp_t x_layer = xc; + double l_thick = layering.layer(l_num-1)->thickness(); + string l_name = _toString(layerType,"layer%d"); + int l_repeat = x_layer.repeat(); + Volume l_vol(l_name,PolyhedraRegular(numsides,rmin,rmax,l_thick),air); + vector<PlacedVolume> sensitives; + + int s_num = 1; + double sliceZ = -l_thick/2; + for(xml_coll_t xs(x_layer,_U(slice)); xs; ++xs) { + xml_comp_t x_slice = xs; + string s_name = _toString(s_num,"slice%d"); + double s_thick = x_slice.thickness(); + Material s_mat = description.material(x_slice.materialStr()); + Volume s_vol(s_name,PolyhedraRegular(numsides,rmin,rmax,s_thick),s_mat); + + s_vol.setVisAttributes(description.visAttributes(x_slice.visStr())); + sliceZ += s_thick/2; + PlacedVolume s_phv = l_vol.placeVolume(s_vol,Position(0,0,sliceZ)); + s_phv.addPhysVolID("slice",s_num); + if ( x_slice.isSensitive() ) { + sens.setType("calorimeter"); + s_vol.setSensitiveDetector(sens); + sensitives.push_back(s_phv); + } + sliceZ += s_thick/2; + s_num++; + } + l_vol.setVisAttributes(description.visAttributes(x_layer.visStr())); + if ( l_repeat <= 0 ) throw std::runtime_error(x_det.nameStr()+"> Invalid repeat value"); + for(int j=0; j<l_repeat; ++j) { + string phys_lay = _toString(l_num,"layer%d"); + layerZ += l_thick/2; + DetElement layer_elt(endcap, phys_lay, l_num); + PlacedVolume pv = endcapVol.placeVolume(l_vol,Position(0,0,layerZ)); + pv.addPhysVolID("layer", l_num); + layer_elt.setPlacement(pv); + for(size_t ic=0; ic<sensitives.size(); ++ic) { + PlacedVolume sens_pv = sensitives[ic]; + DetElement comp_elt(layer_elt,sens_pv.volume().name(),l_num); + comp_elt.setPlacement(sens_pv); + } + layerZ += l_thick/2; + ++l_num; + } + ++layerType; + } + + double z_pos = zmin+totalThickness/2; + PlacedVolume pv; + // Reflect it. + if ( reflect ) { + Assembly assembly(det_name); + DetElement both_endcaps(det_name,det_id); + Volume motherVol = description.pickMotherVolume(both_endcaps); + DetElement sdetA = endcap; + Ref_t(sdetA)->SetName((det_name+"_A").c_str()); + DetElement sdetB = endcap.clone(det_name+"_B",x_det.id()); + + pv = assembly.placeVolume(endcapVol,Transform3D(RotationZYX(M_PI/numsides,0,0), + Position(0,0,z_pos))); + pv.addPhysVolID("barrel", 1); + sdetA.setPlacement(pv); + + pv = assembly.placeVolume(endcapVol,Transform3D(RotationZYX(M_PI/numsides,M_PI,0), + Position(0,0,-z_pos))); + pv.addPhysVolID("barrel", 2); + sdetB.setPlacement(pv); + + pv = motherVol.placeVolume(assembly); + pv.addPhysVolID("system", det_id); + both_endcaps.setPlacement(pv); + both_endcaps.add(sdetA); + both_endcaps.add(sdetB); + return both_endcaps; + } + Volume motherVol = description.pickMotherVolume(endcap); + pv = motherVol.placeVolume(endcapVol,Transform3D(RotationZYX(M_PI/numsides,0,0), + Position(0,0,z_pos))); + pv.addPhysVolID("system", det_id); + pv.addPhysVolID("barrel", 1); + endcap.setPlacement(pv); + Ref_t(endcap)->SetName(det_name.c_str()); + return endcap; +} + +DECLARE_DETELEMENT(DD4hep_PolyhedraEndcapCalorimeter2,create_detector) +DECLARE_DEPRECATED_DETELEMENT(PolyhedraEndcapCalorimeter2,create_detector) diff --git a/src/GenericDetectors/compact/N2_index_of_refraction.json b/src/GenericDetectors/compact/N2_index_of_refraction.json new file mode 100644 index 0000000000000000000000000000000000000000..5ee3707f42e4d836b99fa6698584941aaddec20d --- /dev/null +++ b/src/GenericDetectors/compact/N2_index_of_refraction.json @@ -0,0 +1,35 @@ +{ + "material": "N2 gas", + "property": "index of refraction", + "data": { + "wavelengths": [ + 150.0, 160.0, 170.0, 180.0, 190.0, 200.0, 210.0, 220.0, 230.0, 240.0, 250.0, + 260.0, 270.0, 280.0, 290.0, 300.0, 310.0, 320.0, 330.0, 340.0, 350.0, 360.0, + 370.0, 380.0, 390.0, 400.0, 410.0, 420.0, 430.0, 440.0, 450.0, 460.0, 470.0, + 480.0, 490.0, 500.0, 510.0, 520.0, 530.0, 540.0, 550.0, 560.0, 570.0, 580.0, + 590.0, 600.0, 610.0, 620.0, 630.0, 640.0, 650.0, 660.0, 670.0, 680.0, 690.0, + 700.0, 710.0, 720.0, 730.0, 740.0, 750.0, 760.0, 770.0, 780.0, 790.0, + 800.0], + "energies": [8.26667, 7.75, 7.29412, 6.88889, 6.52632, 6.2, 5.90476, + 5.63636, 5.3913, 5.16667, 4.96, 4.76923, 4.59259, 4.42857, 4.27586, + 4.13333, 4.0, 3.875, 3.75758, 3.64706, 3.54286, 3.44444, 3.35135, + 3.26316, 3.17949, 3.1, 3.02439, 2.95238, 2.88372, 2.81818, 2.75556, + 2.69565, 2.6383, 2.58333, 2.53061, 2.48, 2.43137, 2.38462, 2.33962, + 2.2963, 2.25455, 2.21429, 2.17544, 2.13793, 2.10169, 2.06667, + 2.03279, 2.0, 1.96825, 1.9375, 1.90769, 1.87879, 1.85075, 1.82353, + 1.7971, 1.77143, 1.74648, 1.72222, 1.69863, 1.67568, 1.65333, + 1.63158, 1.61039, 1.58974, 1.56962, 1.55], + "index of refraction": [1.0004, 1.00038, + 1.00037, 1.00036, 1.00035, 1.00034, 1.00034, 1.00033, 1.00033, + 1.00032, 1.00032, 1.00032, 1.00032, 1.00032, 1.00031, 1.00031, + 1.00031, 1.00031, 1.00031, 1.00031, 1.00031, 1.00031, 1.00031, + 1.00031, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, + 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, + 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, + 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, + 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, 1.0003, + 1.0003, 1.0003, 1.0003] + }, + "comments": "Standard conditions: 0 °C; 101 325 Pa.", + "reference": "https://refractiveindex.info/?shelf=main&book=N2&page=Griesmann" +} diff --git a/src/GenericDetectors/compact/forward_rich.xml b/src/GenericDetectors/compact/forward_rich.xml new file mode 100644 index 0000000000000000000000000000000000000000..729cee7740488398aa3c9c7c55bfdef039286ff4 --- /dev/null +++ b/src/GenericDetectors/compact/forward_rich.xml @@ -0,0 +1,154 @@ +<lccdd> + <info name="beam_pipe" title="Beam pipe test" + author="Whitney Armstrong" + url="https://eicweb.phy.anl.gov/EIC/NPDet" + status="development" + version=""> + <comment>EIC Beam Pipe</comment> + </info> + + <includes> + <gdmlFile ref="elements.xml"/> + <gdmlFile ref="materials.xml"/> + <gdmlFile ref="material_properties.xml"/> + </includes> + + <material_properties> + + <material_properties_table name="N2derp"> + <data file="N2_index_of_refraction.json" /> + </material_properties_table> + + </material_properties> + + <derps> + </derps> + + <define> + <constant name="world_side" value="10*m"/> + <constant name="world_x" value="world_side"/> + <constant name="world_y" value="world_side"/> + <constant name="world_z" value="10*world_side"/> + + <constant name="CrossingAngle" value="0.020*rad"/> + + <constant name="CentralBeamPipe_length" value="50.0*cm"/> + <constant name="CentralBeamPipe_thickness" value="0.1*mm"/> + <constant name="CentralBeamPipe_radius" value="3.5*cm"/> + <constant name="CentralBeamPipe_z" value="0.0*cm"/> + + <constant name="UpStreamBeamPipe_length" value="150.0*cm"/> + <constant name="UpStreamBeamPipe_thickness" value="0.1*mm"/> + <constant name="UpStreamBeamPipe_radius" value="5*cm"/> + <constant name="UpStreamBeamPipe_z" value="-1.0*(CentralBeamPipe_length+UpStreamBeamPipe_length)/2.0"/> + + <constant name="DownStreamBeamPipe_length" value="250.0*cm"/> + <constant name="DownStreamBeamPipe_thickness" value="0.1*mm"/> + <constant name="DownStreamBeamPipe_radius" value="5*cm"/> + <constant name="DownStreamBeamPipe_z" value="1.0*(CentralBeamPipe_length+DownStreamBeamPipe_length)/2.0"/> + + <constant name="Place_Center" value="0*cm"/> + + <constant name="ForwardTrackerPlane_z0" value="400*cm"/> + + <constant name="tracking_region_radius" value="1.0*cm"/> + <constant name="tracking_region_zmax" value="2.0*m"/> + <constant name="tracker_region_rmax" value="tracking_region_radius"/> + <constant name="tracker_region_zmax" value="tracking_region_zmax"/> + + + </define> + + <limits> + <limitset name="cal_limits"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + </limitset> + <limitset name="SiTrackerBarrelRegionLimitSet"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + <limit name="track_length_max" particles="*" value="5.0" unit="mm" /> + <limit name="time_max" particles="*" value="5.0" unit="ns" /> + <limit name="ekin_min" particles="*" value="0.01" unit="MeV" /> + <limit name="range_min" particles="*" value="5.0" unit="mm" /> + </limitset> + </limits> + <regions> + <region name="SiTrackerBarrelRegion" eunit="MeV" lunit="mm" cut="0.001" threshold="0.001"> + <limitsetref name="SiTrackerBarrelRegionLimitSet"/> + </region> + </regions> + + <limits> + <limitset name="Tracker_limits"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + </limitset> + </limits> + + <comment>Common Generic visualization attributes</comment> + <display> + <vis name="InvisibleNoDaughters" showDaughters="false" visible="false"/> + <vis name="InvisibleWithDaughters" showDaughters="true" visible="false"/> + <vis name="GreenVis" alpha="0.5" r= "0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="RedVis" alpha="0.5" r= "1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="BlueVis" alpha="0.5" r= "0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="OrangeVis" alpha="0.5" r= "1.0" g="0.45" b="0.0" showDaughters="true" visible="true"/> + <vis name="RedGreenVis" alpha="0.5" r= "1.0" g="1.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="BlueGreenVis" alpha="0.5" r= "0.0" g="1.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="PurpleVis" alpha="0.5" r= "1.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="DoubleRedG" alpha="0.5" r= "2.0" g=".10" b="0.0" showDaughters="true" visible="true"/> + <vis name="RBG015" alpha="0.5" r= "0.0" g=".2" b="1.0" showDaughters="true" visible="true"/> + <vis name="RBG510" alpha="0.5" r= "1.0" g=".2" b="0.0" showDaughters="true" visible="true"/> + <vis name="RBG" alpha="0.5" r= "1.0" g="1.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="GrayVis" alpha="1.0" r= "0.75" g="0.75" b="0.75" showDaughters="true" visible="true"/> + </display> + +<detectors> + + <detector id="5" name="ForwardRich" type="ForwardRICH" vis="BlueGreenVis" > + </detector> + + <!-- + <detector + name="TrackerPlanes" + type="ForwardPlaneTracker" + vis="PurpleVis" + id="1" + limits="Tracker_limits" + readout="TPCollection" + insideTrackingVolume="true"> + + <layer z_offset="ForwardTrackerPlane_z0+10*cm" id="0"> + <module name="SinglePlane" thickness="1*cm" x_half="100*cm" y_half="100*cm" material="CarbonFiber" vis="BlueVis" /> + <sensitive thickness="0.1*mm" material="Air" vis="GreenVis" /> + </layer> + <layer z_offset="ForwardTrackerPlane_z0+20*cm" id="1"> + <module name="SinglePlane1" thickness="1*cm" x_half="100*cm" y_half="100*cm" material="CarbonFiber" vis="BlueVis" /> + <sensitive thickness="0.1*mm" material="Air" vis="GreenVis" /> + </layer> + <layer z_offset="ForwardTrackerPlane_z0+40*cm" id="2"> + <module name="SinglePlane2" thickness="1*cm" x_half="100*cm" y_half="100*cm" material="CarbonFiber" vis="BlueVis" /> + <sensitive thickness="0.1*mm" material="Air" vis="GreenVis" /> + </layer> + <layer z_offset="ForwardTrackerPlane_z0+70*cm" id="3"> + <module name="SinglePlane3" thickness="1*cm" x_half="100*cm" y_half="100*cm" material="CarbonFiber" vis="BlueVis" /> + <sensitive thickness="0.1*mm" material="Air" vis="GreenVis" /> + </layer> + <layer z_offset="ForwardTrackerPlane_z0+100*cm" id="4"> + <module name="SinglePlane4" thickness="1*cm" x_half="100*cm" y_half="100*cm" material="CarbonFiber" vis="BlueVis" /> + <sensitive thickness="0.1*mm" material="Air" vis="GreenVis" /> + </layer> + </detector> +--> + + </detectors> + + <!-- + <readouts> + <readout name="TPCollection"> + <segmentation type="CartesianGridXY" grid_size_x="10.0*cm" grid_size_y="10.0*cm" /> + <id>system:5,layer:9,module:8,x:32:-16,y:-16</id> + </readout> + </readouts> + +--> + +</lccdd> diff --git a/src/GenericDetectors/compact/gps.mac b/src/GenericDetectors/compact/gps.mac new file mode 100644 index 0000000000000000000000000000000000000000..a590d9270d50c53926d5ad601904287ea41329a2 --- /dev/null +++ b/src/GenericDetectors/compact/gps.mac @@ -0,0 +1,23 @@ +/control/verbose 2 +/run/initialize + +/gps/verbose 2 +/gps/particle e- +/gps/number 3 + +/gps/ene/type Gauss +/gps/ene/mono 9.0 GeV +/gps/ene/sigma 3.0 GeV + +/gps/pos/type Volume +/gps/pos/shape Cylinder +/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 0 + +#/gps/direction 1 0 0 +/gps/ang/type iso +/gps/ang/mintheta 178 deg + +/run/beamOn 100 diff --git a/src/GenericDetectors/compact/material_properties.xml b/src/GenericDetectors/compact/material_properties.xml new file mode 100644 index 0000000000000000000000000000000000000000..1415e40d36f2dda06c60bde46043fcaa49a624a7 --- /dev/null +++ b/src/GenericDetectors/compact/material_properties.xml @@ -0,0 +1,29 @@ +<?xml version="1.0" encoding="UTF-8"?> +<material_properties> + + <material_properties_table name="N2"> + <data file="N2_index_of_refraction.json" /> + </material_properties_table> + <!-- + <material name="PolystyreneFoam"> + <D type="density" value="0.0056" unit="g/cm3"/> + <fraction n="1.0" ref="Polystyrene"/> + </material> + + <material name="Kapton"> + <D value="1.43" unit="g/cm3" /> + <composite n="22" ref="C"/> + <composite n="10" ref="H" /> + <composite n="2" ref="N" /> + <composite n="5" ref="O" /> + </material> + + <material name="PEEK"> + <D value="1.37" unit="g/cm3" /> + <composite n="19" ref="C"/> + <composite n="12" ref="H" /> + <composite n="3" ref="O" /> + </material> + --> + +</material_properties> diff --git a/src/GenericDetectors/compact/run_example b/src/GenericDetectors/compact/run_example new file mode 100755 index 0000000000000000000000000000000000000000..dc29f6723fb84546ee1c3876fb11f38065404ae0 --- /dev/null +++ b/src/GenericDetectors/compact/run_example @@ -0,0 +1,32 @@ +#!/bin/bash + +rm -f forward_rich.slcio +rm -f forward_rich.root + +#ddsim --runType run --enableG4GPS \ +# --macroFile gps.mac \ +# --compactFile detectors/SiD/slic/sieic3/sieic3_compact.xml \ +# --outputFile slic_ddsim_out.slcio + +#ddsim --runType run --enableG4GPS \ +# --macroFile gps.mac \ +# --compactFile detectors/SiD/compact/sid_working/sidloi3_v00.xml \ +# --outputFile example1_ddsim_out.slcio + + +#ddsim --runType run --enableG4GPS \ +# --macroFile gps.mac \ +# --compactFile ./gem_tracker_disc.xml \ +# --outputFile test_tracker_disc.root + +ddsim --runType run --enableG4GPS \ + --macroFile gps.mac \ + --compactFile ./roman_pots.xml \ + --outputFile roman_pot_test.root + +#anajob example1_ddsim_out.slcio | tail -n 30 > collections.txt + +#rm outputfile.slcio +#Marlin steering_files/basicsteering.xml +# +#anajob outputfile.slcio | tail -n 30 > collections2.txt diff --git a/src/GenericDetectors/include/SiDigiHitInfo.h b/src/GenericDetectors/include/SiDigiHitInfo.h new file mode 100644 index 0000000000000000000000000000000000000000..7aed147fc6cc687e615e96992d88f1bf700c987f --- /dev/null +++ b/src/GenericDetectors/include/SiDigiHitInfo.h @@ -0,0 +1,22 @@ +#ifndef SiDigiHitInfo_HH +#define SiDigiHitInfo_HH + + + +namespace npdet { + namespace digi { + + struct SiDigiHitInfo { + int64_t volumeID; + int64_t cellID; + double pathLength; + double EDep; + ROOT::Math::XYZTVector position; + ROOT::Math::XYZTVector momentum; + }; + + } +} + +#endif + diff --git a/src/GenericDetectors/trackers/src/RPC_sim_7.f b/src/GenericDetectors/trackers/src/RPC_sim_7.f new file mode 100644 index 0000000000000000000000000000000000000000..354cfd4814622e738981499b7e7516235b87f296 --- /dev/null +++ b/src/GenericDetectors/trackers/src/RPC_sim_7.f @@ -0,0 +1,342 @@ + program rpc_sim + +* Author J.Repond. Version of April 21, 2010 +* added transfer of interaction layer information to atuple + + implicit none + + integer i1,i2,i3,i4,i5,i6 + integer i_event,n_event,n_point + integer i_pad_x,i_pad_y + integer hit(0:350,0:350,0:150) + integer z_store(0:10000),z_tmp(10000),f_tmp(10000),n_tmp + integer i_int,n_int,inter_layer,inter_layer_new + + real r1,r2,r3,r4,r5 + real e,x,y,z,Q + real T,Q0 + real par(3) + real x_store(0:10000),y_store(0:10000) + real x_tmp(10000),y_tmp(10000),dist_cut,e_cut + real x_int,y_int,radius_max + real q_int,q_tot,sigma1,sigma2,ratio + real q_integral1,q_integral2,q_integral,q_check + real q_pad(0:600,0:600,0:600) + real x_deposit,y_deposit + real x_crack(2),x_scaler,y_crack(4),y_scaler + real dist(10000) + + character*12 file0 + character*3 file1 + character*7 file2 + character*4 file3 + character*25 file4 + +* data par /0.00057,1.735,0.0012/ ! 6.2 kV +* data par /0.00007,1.951,0.0010/ ! 6.3 kV data par /0.00001,2.167,0.0008/ ! 6.4 kV + data par /0.00001,2.167,0.0008/ ! 6.4 kV + + +* units are cm + + data n_int /1000/ ! Number of points in distributed charge + data radius_max /4.0/ ! Maximum extent of charge, default = 4.0 + data y_crack /0.005,32.525,65.065,97.595/ ! Cracks in y + +* the famous 6 parameters + + data sigma1 /0.078/ ! Width of narrow Gaussian + data sigma2 /0.28/ ! Width of wide Gaussian + data ratio /0.310/ ! Ratio of contribution from slope1 and slope2; default = 0.345 + data T /0.27/ ! Threshold (in pC); default = 0.3645 + data dist_cut /0.092/ ! Minimum distance between hits in x,y plane; default = 0.092 + data e_cut /0.15/ ! Distance at which charge not attenuated anymore, value needs to be > dist_cut + data Q0 /+0.125/ ! offset in charge; default = +0.201 + data file2 /'.atuple'/ + +* MC output files + +* data file0 /'../mc1/Muon.'/ + data file0 /'../mc1/Posi.'/ +* data file0 /'../mc1/Pion.'/ + data file3 /'.txt'/ +* data file4 /'/HEPdata/ILC/repond/Posi.'/ + data file4 /'/HEPdata/ILC/repond/Muon.'/ +* data file4 /'/HEPdata/ILC/repond/Test.'/ + +* check distance cuts + + if(dist_cut.gt.e_cut) then + write(6,*) 'Problem with distance cut values' + stop + endif + +* setup + + i_event=0 + n_event=0 + + write(6,*) 'Type run number' + read(5,*) file1 +* open(unit=61,status='OLD',file=file0//file1//file3) + open(unit=61,status='OLD',file=file4//file1//file3) ! for Posi.007.txt and beyond + open(unit=62,status='UNKNOWN',file=file4//file1//file2) ! the output atuple + + read(61,*) z_store(0),x_store(0),y_store(0) + inter_layer=x_store(0) + +* calculate integral over charge distribution + + q_integral1=sigma1**2 + q_integral2=sigma2**2 + +** get x,y,z of energy deposit + + 1 continue + i_event=i_event+1 + n_event=n_event+1 + if(n_event.ge.100) then + write(6,*) i_event + n_event=0 + endif + if(i_event.gt.60000) stop + +* read in from file (generated by MC) + + n_point=1 + 5 continue + + read(61,*,err=4,end=4) z_store(n_point), + > x_store(n_point),y_store(n_point) + + if(z_store(n_point).ge.0) then + x_store(n_point)=x_store(n_point)/10. + y_store(n_point)=y_store(n_point)/10. + n_point=n_point+1 + goto 5 + else + x_store(0)=x_store(n_point) + y_store(0)=y_store(n_point) + z_store(0)=z_store(n_point) + n_point=n_point-1 + endif + inter_layer_new=x_store(0) + +* clean up hits (eliminate close by hits) + + do i1=1,n_point + f_tmp(i1)=1 ! 0 = discard point, 1 = keep point + enddo + + do i1=1,n_point-1 + if(f_tmp(i1).eq.1) then + do i2=i1+1,n_point + if(z_store(i1).eq.z_store(i2).and. + > f_tmp(i2).eq.1) then + r1=sqrt((x_store(i1)-x_store(i2))**2+ + > (y_store(i1)-y_store(i2))**2) + if(r1.lt.dist_cut) f_tmp(i2)=0 + endif + enddo + endif + enddo + + n_tmp=0 + do i1=1,n_point + if(f_tmp(i1).eq.1) then + n_tmp=n_tmp+1 + x_tmp(n_tmp)=x_store(i1) + y_tmp(n_tmp)=y_store(i1) + z_tmp(n_tmp)=z_store(i1) + endif + enddo + + do i1=1,n_tmp + x_store(i1)=x_tmp(i1) + y_store(i1)=y_tmp(i1) + z_store(i1)=z_tmp(i1) + enddo + n_point=n_tmp + +* calculate distance of each point to closest point + + do i1=1,n_point + dist(i1)=999. + do i2=1,n_point + if(i2.ne.i1) then + if(z_store(i1).eq.z_store(i2)) then + r1=sqrt((x_store(i1)-x_store(i2))**2+(y_store(i1)-y_store(i2))**2) + if(r1.lt.dist(i1)) then + dist(i1)=r1 + endif + endif + endif + enddo + enddo + +* reset pad array and pad charges + + do i1=0,150 + do i2=0,150 + do i3=0,150 + hit(i3,i2,i1)=0 + q_pad(i3,i2,i1)=0. + enddo + enddo + enddo + +* loop over all points in event + + do i1=1,n_point + +* get x,y + + x=x_store(i1) + y=y_store(i1) + +* determine id of pad hit + + i_pad_x=int(x) + i_pad_y=int(y) + +* determine minimum distance to boundaries in x and y + + x_crack(1)=x+48.0 + x_crack(2)=96.0-(x+48.0) + x_scaler=amin1(x_crack(1),x_crack(2)) + + y_scaler=999. + do i2=1,4 + r1=abs(y-y_crack(i2)+48.8) + if(r1.le.y_scaler) y_scaler=r1 + enddo + +* generate charge: Distributed like RABBIT data + + 2 r1=rand(0)*12000. + r2=par(1)*r1**par(2)*exp(-r1*par(3)) + r3=rand(0)*40. + if(r2.gt.r3) then + Q=r1/1100. + else + goto 2 + endif + +* Introduce offset + + Q=Q-Q0 + +* Adjust charge if close to other points + + if(dist(i1).gt.dist_cut.and.dist(i1).le.e_cut) then + r1=0.5*(1.+(dist(i1)-dist_cut)/(e_cut-dist_cut)) + Q=r1*Q + endif + +* Adjust charge close to boundary + + if(x_scaler.lt.10.0) Q=Q*((x_scaler/10.0)**0.7) + if(y_scaler.lt.10.0) Q=Q*((y_scaler/10.0)**0.25) + if(y_scaler.lt.1.0) Q=Q/6. + +* Distribute charge according to STAR measurements + + +* integrate with n_int points + + i_int=0 + q_check=0. + 3 continue + x_int=2.*radius_max*(rand(0)-.5) + y_int=2.*radius_max*(rand(0)-.5) + r1=sqrt(x_int**2+y_int**2) + if(r1.gt.radius_max) goto 3 + i_int=i_int+1 + +* calculate charge + + q_int=Q/real(n_int)*8.* + > ((1.-ratio)/q_integral1*exp(-0.5*((r1/sigma1)**2))+ + > ratio/q_integral2*exp(-0.5*((r1/sigma2)**2))) + + q_check=q_check+q_int +* write(63,*) Q,q_int,r1 + +* deposit charge on corresponding pad + + x_deposit=x+x_int+48.0 ! x of deposited charge in cm, centered onto plane + y_deposit=y+y_int+48.8 ! y of deposited charge in cm, centered onto plane + + if(x_deposit.lt.0.0.or.x_deposit.gt.96.0) q_int=0. ! set charge to zero if off the boards in x + + i2=int(x_deposit) ! pad index in x + + if(y_deposit.lt.0.27) then + q_int=0. ! set charge to zero if off the boards in y + elseif(y_deposit.gt.32.27.and.y_deposit.lt.32.8) then + q_int=0. + elseif(y_deposit.gt.64.8.and.y_deposit.lt.65.33) then + q_int=0. + elseif(y_deposit.gt.97.33) then + q_int=0. + endif + + if(y_deposit.le.32.27) then ! pad index in y + i3=int(y_deposit-0.27) + elseif(y_deposit.le.64.80) then + i3=int(y_deposit-0.8) + elseif(y_deposit.le.97.33) then + i3=int(y_deposit-1.33) + else + i3=0 + endif + +* write(70,*) 'x,y,x_int,y_int',x,y,x_int,y_int +* write(70,*) 'x_deposit,y_deposit',x_deposit,y_deposit +* write(70,*) 'padx,pady',i2,i3 + + i4=z_store(i1) ! pad index in z + + q_pad(i2,i3,i4)=q_pad(i2,i3,i4)+q_int + if(i_int.lt.n_int) goto 3 + + enddo + +* determine pads hit + + do i1=0,150 + do i2=0,150 + do i3=0,150 + if(q_pad(i1,i2,i3).gt.T) hit(i1,i2,i3)=1 + enddo + enddo + enddo + +* write out hits + + write(62,*) i_event,inter_layer,-1,-1 + do i3=0,150 + do i2=0,200 + do i1=0,200 + if(hit(i1,i2,i3).eq.1) then + write(62,*) i_event,i1,i2,i3 + endif + enddo + enddo + enddo + +* check total charge + + write(60,*) q,q_check + +* get next event + + inter_layer=inter_layer_new + goto 1 + +* finish up + + 4 continue + + stop + end diff --git a/src/tracking/CMakeLists.txt b/src/tracking/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..894af1c6d2ab3caedcf647ba6094d4a72bf23243 --- /dev/null +++ b/src/tracking/CMakeLists.txt @@ -0,0 +1,64 @@ +cmake_minimum_required(VERSION 3.3 FATAL_ERROR) + +# see https://rix0r.nl/blog/2015/08/13/cmake-guide/ +# Must use GNUInstallDirs to install libraries into correct +# locations on all platforms. +include(GNUInstallDirs) + +find_package( LCIO REQUIRED ) + +find_package( DD4hep REQUIRED COMPONENTS DDCore DDG4 DDRec ) + +find_package(ROOT REQUIRED COMPONENTS Geom GenVector MathCore) +include(${ROOT_USE_FILE}) + +FIND_PACKAGE( GBL ) +FIND_PACKAGE( aidaTT ) + +IF( GBL_FOUND ) + ADD_DEFINITIONS( "-DUSE_GBL" ) + ADD_DEFINITIONS( "-DGBL_EIGEN_SUPPORT_ROOT" ) +ELSE() + MESSAGE( STATUS "GBL not found, track fitting with GBL will not be available." ) +ENDIF() + +IF( DD4hep_FOUND ) + ADD_DEFINITIONS( "-DAIDATT_USE_DD4HEP" ) +ELSE() + MESSAGE( STATUS "DD4hep not found, the geometry decription from DD4hep will not be available." ) +ENDIF() + +IF( LCIO_FOUND ) + ADD_DEFINITIONS( "-DUSE_LCIO" ) +ELSE() + MESSAGE( STATUS "LCIO not found, LCIO persistency will not be available." ) +ENDIF() + +IF( streamlog_FOUND ) + ADD_DEFINITIONS( "-DAIDATT_USE_STREAMLOG" ) +ENDIF() + +set(exe_name gbl_example) +add_executable(${exe_name} src/${exe_name}.cxx) +target_include_directories(${exe_name} PUBLIC + PRIVATE include ) +target_compile_features(${exe_name} + PUBLIC cxx_auto_type + PUBLIC cxx_trailing_return_types + PRIVATE cxx_variadic_templates ) + +target_link_libraries( ${exe_name} + PUBLIC DDCore DDG4 DDRec ROOT::Core ROOT::Hist ROOT::GenVector ROOT::Gpad + ${LCIO_LIBRARIES} ${aidaTT_LIBRARIES} ${GBL_LIBRARIES} ) + +message("aidaTT_LIBRARIES : ${aidaTT_LIBRARIES}") + +INCLUDE_DIRECTORIES( SYSTEM ${aidaTT_INCLUDE_DIRS} ${GBL_INCLUDE_DIRS} ) +#$TARGET_LINK_LIBRARIES( ${exe_name} ${GBL_LIBRARIES} ) +ADD_DEFINITIONS( ${${GBL_LIBRARIES}_DEFINITIONS} ) + +install(TARGETS ${exe_name} + EXPORT NPDetTargets + RUNTIME DESTINATION bin ) + + diff --git a/src/tracking/src/gbl_example.cxx b/src/tracking/src/gbl_example.cxx new file mode 100644 index 0000000000000000000000000000000000000000..580ca2dd39f15f3b24cbadf40027def298e26ba7 --- /dev/null +++ b/src/tracking/src/gbl_example.cxx @@ -0,0 +1,222 @@ +/* + * this is an example of the usage + * it's based on the example1.cpp in GBL by Claus Kleinwort + */ + +#define AIDATT_USE_DD4HEP 1 +#define USE_LCIO 1 +#define USE_GBL 1 + +// c++ +#include <iostream> +#include <map> +#include <string> + +// lcio +#include "lcio.h" +#include "IO/LCWriter.h" +#include "EVENT/LCEvent.h" +#include "EVENT/LCCollection.h" +#include "EVENT/SimTrackerHit.h" +#include "UTIL/LCTrackerConf.h" + +#include <IMPL/LCCollectionVec.h> +#include "IMPL/TrackImpl.h" + +// DD4hep +#include "DD4hep/Detector.h" +#include "DD4hep/DD4hepUnits.h" +#include "DDRec/SurfaceHelper.h" +#include "DDRec/SurfaceManager.h" +#include "DDRec/Surface.h" +#include "DDRec/CellIDPositionConverter.h" + + + +// aidaTT +#include "aidaTT/DD4hepGeometry.hh" +#include "aidaTT/AidaTT.hh" +#include "aidaTT/ConstantSolenoidBField.hh" +#include "aidaTT/analyticalPropagation.hh" +#include "aidaTT/simplifiedPropagation.hh" +#include "aidaTT/GBLInterface.hh" +#include "aidaTT/fitResults.hh" + +using namespace std; +using namespace lcio; + + +int main(int argc, char** argv) +{ + if(argc < 2) + { + std::cout << " usage: ./gbl_example ILDEx.xml ILDExSimu.slcio" << std::endl ; + return 1; + } + + /// dd4hep stuff + std::string inFile = argv[1] ; + + /// preamble: load the geo info, get all surfaces => entry point for intersection calculation + dd4hep::Detector& theDet = dd4hep::Detector::getInstance(); + theDet.fromCompact(inFile); + + dd4hep::DetElement world = theDet.world() ; + + aidaTT::DD4hepGeometry geom(theDet); + + const auto& surfaces = geom.getSurfaces() ; + + dd4hep::rec::CellIDPositionConverter cellid_converter(theDet); + + dd4hep::rec::SurfaceManager& surfMan = *theDet.extension<dd4hep::rec::SurfaceManager>() ; + //auto surfMap = surfMan.map( "world" ) ; + + // create map of surfaces + std::map< long64, const aidaTT::ISurface* > surfMap ; + for(const auto& surf : surfaces) + { + surfMap[surf->id() ] = (surf) ; + //surfMap[surf->id() ]->setProperty(dd4hep::rec::SurfaceType::OrthogonalToZ,true); + } + + + + /// lcio stuff + std::string lcioFileName = argv[2] ; + + LCReader* rdr = LCFactory::getInstance()->createLCReader() ; + rdr->open(lcioFileName) ; + + std::string outFileName = "trackTest.slcio"; + + LCWriter* wrt = LCFactory::getInstance()->createLCWriter() ; + wrt->open(outFileName) ; + + LCEvent* evt = 0 ; + + std::vector< std::string > colNames ; + colNames.push_back("SiTrackerBarrelHits") ; + //colNames.push_back("SiVertexBarrelHits") ; + //colNames.push_back("SiTrackerEndcapHits") ; + //colNames.push_back("SiVertexEndcapHits") ; + + UTIL::BitField64 idDecoder(LCTrackerCellID::encoding_string()) ; + + + // the aidaTT stuff + + /// create some bogus track parameters: + ///~ test data is p = 3./sqrt(2) * (1 1 .1)^T + /// => d0 = 0; tanLambda = 1/sqrt(2), phi0 = pi/4, z0 = 0 + // with B = 3.5T => omega = 3* 10^(-4) * 3.5 / sqrt(6) + aidaTT::trackParameters bogusTP; + + //~ const double initialOmega = 3.5*1e-5; + //~ const double initialTanLambda = 0.; + //~ const double initialPhi = 5*M_PI_4; + + const double initialOmega = 1.5 * 1e-4; + const double initialTanLambda = 1.; + const double initialPhi = M_PI_2; + + bogusTP.setTrackParameters(aidaTT::Vector5(initialOmega, initialTanLambda, initialPhi, 0., 0.)); + + // create the different objects needed for fitting + // first a constant field parallel to z, 1T + aidaTT::ConstantSolenoidBField* bfield = new aidaTT::ConstantSolenoidBField(3.5); + + // create the propagation object + aidaTT::analyticalPropagation* propagation = new aidaTT::analyticalPropagation(); + //aidaTT::simplifiedPropagation* propagation = new aidaTT::simplifiedPropagation(); + + // create the fitter object + aidaTT::GBLInterface* fitter = new aidaTT::GBLInterface(); + + + while((evt = rdr->readNextEvent()) != 0) + { + + IMPL::LCCollectionVec* outputTrackCollection = new IMPL::LCCollectionVec(EVENT::LCIO::TRACK); + std::string outName = "trytracks"; + + std::cout << " reading an event " << std::endl; + // now create a trajectory object to be fitted + aidaTT::trajectory theMaster(bogusTP, fitter, propagation, &geom); + + for(unsigned icol = 0, ncol = colNames.size() ; icol < ncol ; ++icol) + { + LCCollection* col = NULL; + try + { + col = evt->getCollection(colNames[ icol ]) ; + } + catch(DataNotAvailableException &e) + { + continue; + } + + int nHit = col->getNumberOfElements() ; + + for(int i = 0 ; i < nHit ; ++i) + { + SimTrackerHit* sHit = (SimTrackerHit*) col->getElementAt(i) ; + long64 id = sHit->getCellID0() ; + //idDecoder.setValue(id) ; + + //const auto si = surfMap->find( cellid_converter.findContext(id)->identifier ); + // // Get the Surface (http://test-dd4hep.web.cern.ch/test-dd4hep/doxygen/html/classdd4hep_1_1rec_1_1_i_surface.html) + //std::cout << surfMap->count( cellid_converter.findContext(id)->identifier ) << "\n"; + //dd4hep::rec::ISurface* surf = (si != surfMap->end() ? si->second : nullptr); + + const aidaTT::ISurface* surf = surfMap[ id ] ; + + dd4hep::rec::Vector3D pos1; + pos1 = surf->origin(); + double recalcPos[3] = {pos1.x(), pos1.y(), pos1.z()}; // could it be more ugly? + if(surf->type().isSensitive()) { + + std::vector<double> precisionDummy; + precisionDummy.push_back(1. / 0.0001); + precisionDummy.push_back(1. / 0.00012); + //for(unsigned int i = 0; i < 3; ++i) + // recalcPos[i] = sHit->getPosition()[i] * dd4hep::mm; + + theMaster.addMeasurement(recalcPos, precisionDummy, *surf, sHit); + } + } + + theMaster.prepareForFitting(); + theMaster.fit(); + + const aidaTT::fitResults& result = *(theMaster.getFitResults()); + + TrackImpl* theTrack = new TrackImpl(); + theTrack->setOmega(result.estimatedParameters()(0)) ; + theTrack->setTanLambda(result.estimatedParameters()(1)); + theTrack->setPhi(result.estimatedParameters()(2)); + theTrack->setD0(result.estimatedParameters()(3)); + theTrack->setZ0(result.estimatedParameters()(4)); + + TrackImpl* theInitialTrack = new TrackImpl(); + theInitialTrack->setOmega(initialOmega) ; + theInitialTrack->setTanLambda(initialTanLambda); + theInitialTrack->setPhi(initialPhi); + theInitialTrack->setD0(0.); + theInitialTrack->setZ0(0.); + + + outputTrackCollection->addElement(theTrack); + outputTrackCollection->addElement(theInitialTrack); + + } + evt->addCollection(outputTrackCollection , outName); + + wrt->writeEvent(evt); + + } + + + return 0; +} +