forked from sPHENIX-Collaboration/coresoftware
/
PHG4MvtxDigitizer.cc
211 lines (170 loc) · 6.71 KB
/
PHG4MvtxDigitizer.cc
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
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
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
// This is the new trackbase container version
#include "PHG4MvtxDigitizer.h"
// Move to new storage containers
#include <trackbase/TrkrDefs.h>
#include <trackbase/TrkrHitv2.h> // for TrkrHit
#include <trackbase/TrkrHitSet.h>
#include <trackbase/TrkrHitSetContainer.h>
#include <trackbase/TrkrHitTruthAssoc.h>
#include <g4detectors/PHG4CylinderGeom.h>
#include <g4detectors/PHG4CylinderGeomContainer.h>
#include <fun4all/Fun4AllReturnCodes.h>
#include <fun4all/SubsysReco.h> // for SubsysReco
#include <phool/PHCompositeNode.h>
#include <phool/PHNode.h> // for PHNode
#include <phool/PHNodeIterator.h>
#include <phool/PHRandomSeed.h>
#include <phool/getClass.h>
#include <phool/phool.h> // for PHWHERE
#include <gsl/gsl_rng.h> // for gsl_rng_alloc
#include <cstdlib> // for exit
#include <iostream>
#include <set>
using namespace std;
PHG4MvtxDigitizer::PHG4MvtxDigitizer(const string &name)
: SubsysReco(name)
, _energy_threshold(0.95e-6)
{
unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
cout << Name() << " random seed: " << seed << endl;
RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
gsl_rng_set(RandomGenerator, seed);
if (Verbosity() > 0)
cout << "Creating PHG4MvtxDigitizer with name = " << name << endl;
}
PHG4MvtxDigitizer::~PHG4MvtxDigitizer()
{
gsl_rng_free(RandomGenerator);
}
int PHG4MvtxDigitizer::InitRun(PHCompositeNode *topNode)
{
//-------------
// Add Hit Node
//-------------
PHNodeIterator iter(topNode);
// Looking for the DST node
PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
if (!dstNode)
{
cout << PHWHERE << "DST Node missing, doing nothing." << endl;
return Fun4AllReturnCodes::ABORTRUN;
}
CalculateMvtxLadderCellADCScale(topNode);
//----------------
// Report Settings
//----------------
if (Verbosity() > 0)
{
cout << "====================== PHG4MvtxDigitizer::InitRun() =====================" << endl;
for (std::map<int, unsigned short>::iterator iter = _max_adc.begin();
iter != _max_adc.end();
++iter)
{
cout << " Max ADC in Layer #" << iter->first << " = " << iter->second << endl;
}
for (std::map<int, float>::iterator iter = _energy_scale.begin();
iter != _energy_scale.end();
++iter)
{
cout << " Energy per ADC in Layer #" << iter->first << " = " << 1.0e6 * iter->second << " keV" << endl;
}
cout << "===========================================================================" << endl;
}
return Fun4AllReturnCodes::EVENT_OK;
}
int PHG4MvtxDigitizer::process_event(PHCompositeNode *topNode)
{
// This code now only does the Mvtx
DigitizeMvtxLadderCells(topNode);
return Fun4AllReturnCodes::EVENT_OK;
}
void PHG4MvtxDigitizer::CalculateMvtxLadderCellADCScale(PHCompositeNode *topNode)
{
// defaults to 8-bit ADC, short-axis MIP placed at 1/4 dynamic range
PHG4CylinderGeomContainer *geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
if (!geom_container) return;
if (Verbosity())
cout << "Found CYLINDERGEOM_MVTX node" << endl;
PHG4CylinderGeomContainer::ConstRange layerrange = geom_container->get_begin_end();
for (PHG4CylinderGeomContainer::ConstIterator layeriter = layerrange.first;
layeriter != layerrange.second;
++layeriter)
{
int layer = layeriter->second->get_layer();
float thickness = (layeriter->second)->get_pixel_thickness();
float pitch = (layeriter->second)->get_pixel_x();
float length = (layeriter->second)->get_pixel_z();
float minpath = pitch;
if (length < minpath) minpath = length;
if (thickness < minpath) minpath = thickness;
float mip_e = 0.003876 * minpath;
if (Verbosity())
cout << "mip_e = " << mip_e << endl;
if (_max_adc.find(layer) == _max_adc.end())
{
_max_adc[layer] = 255;
_energy_scale[layer] = mip_e / 64;
}
}
return;
}
void PHG4MvtxDigitizer::DigitizeMvtxLadderCells(PHCompositeNode *topNode)
{
//----------
// Get Nodes
//----------
// new containers
//=============
// Get the TrkrHitSetContainer node
TrkrHitSetContainer *trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
if (!trkrhitsetcontainer)
{
cout << "Could not locate TRKR_HITSET node, quit! " << endl;
exit(1);
}
// Get the TrkrHitTruthAssoc node
auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
// Digitization
// We want all hitsets for the Mvtx
TrkrHitSetContainer::ConstRange hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::mvtxId);
for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first;
hitset_iter != hitset_range.second;
++hitset_iter)
{
// we have an itrator to one TrkrHitSet for the mvtx from the trkrHitSetContainer
// get the hitset key so we can find the layer
TrkrDefs::hitsetkey hitsetkey = hitset_iter->first;
int layer = TrkrDefs::getLayer(hitsetkey);
if (Verbosity() > 1) cout << "PHG4MvtxDigitizer: found hitset with key: " << hitsetkey << " in layer " << layer << endl;
// get all of the hits from this hitset
TrkrHitSet *hitset = hitset_iter->second;
TrkrHitSet::ConstRange hit_range = hitset->getHits();
std::set<TrkrDefs::hitkey> hits_rm;
for (TrkrHitSet::ConstIterator hit_iter = hit_range.first;
hit_iter != hit_range.second;
++hit_iter)
{
TrkrHit *hit = hit_iter->second;
// Convert the signal value to an ADC value and write that to the hit
//unsigned int adc = hit->getEnergy() / (TrkrDefs::MvtxEnergyScaleup *_energy_scale[layer]);
if (Verbosity() > 0)
cout << " PHG4MvtxDigitizer: found hit with key: " << hit_iter->first << " and signal " << hit->getEnergy() / TrkrDefs::MvtxEnergyScaleup << " in layer " << layer << std::endl;
// Remove the hits with energy under threshold
bool rm_hit = false;
if ( (hit->getEnergy() / TrkrDefs::MvtxEnergyScaleup) < _energy_threshold) rm_hit = true;
unsigned short adc = (unsigned short) (hit->getEnergy() / (TrkrDefs::MvtxEnergyScaleup *_energy_scale[layer]));
if (adc > _max_adc[layer]) adc = _max_adc[layer];
hit->setAdc(adc);
if (rm_hit) hits_rm.insert(hit_iter->first);
}
for( const auto& key:hits_rm )
{
if (Verbosity() > 0) cout << " PHG4MvtxDigitizer: remove hit with key: " << key << endl;
hitset->removeHit(key);
if( hittruthassoc ) hittruthassoc->removeAssoc(hitsetkey, key);
}
}
// end new containers
//===============
return;
}