Skip to content
Snippets Groups Projects
Commit 62568bec authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

new file: examples/TOC 

	new file:   src/ConceptDetectors/topside/.gitignore
	new file:   src/GenericDetectors/calorimeters/src/PolyhedraEndcapCalorimeter2_geo.cpp
parent bb59cd38
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1072 additions and 0 deletions
examples/TOC 0 → 100644
+ 14
− 0
View file @ 62568bec
- [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
src/ConceptDetectors/topside/.gitignore 0 → 100644
+ 1
− 0
View file @ 62568bec
!*.png
src/GenericDetectors/calorimeters/src/PolyhedraEndcapCalorimeter2_geo.cpp 0 → 100644
+ 134
− 0
View file @ 62568bec
//==========================================================================
// 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)
src/GenericDetectors/compact/N2_index_of_refraction.json 0 → 100644
+ 35
− 0
View file @ 62568bec
{
"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"
}
src/GenericDetectors/compact/forward_rich.xml 0 → 100644
+ 154
− 0
View file @ 62568bec
<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>
src/GenericDetectors/compact/gps.mac 0 → 100644
+ 23
− 0
View file @ 62568bec
/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
src/GenericDetectors/compact/material_properties.xml 0 → 100644
+ 29
− 0
View file @ 62568bec
<?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>
src/GenericDetectors/compact/run_example 0 → 100755
+ 32
− 0
View file @ 62568bec
#!/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
src/GenericDetectors/include/SiDigiHitInfo.h 0 → 100644
+ 22
− 0
View file @ 62568bec
#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
src/GenericDetectors/trackers/src/RPC_sim_7.f 0 → 100644
+ 342
− 0
View file @ 62568bec
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
src/tracking/CMakeLists.txt 0 → 100644
+ 64
− 0
View file @ 62568bec
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 )
src/tracking/src/gbl_example.cxx 0 → 100644
+ 222
− 0
View file @ 62568bec
/*
* 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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment