Line data Source code
1 : #include <assert.h>
2 : #include <math.h>
3 : #include <bitset>
4 : #include <fstream>
5 : #include <iostream>
6 : #include <random>
7 : #include <string>
8 : #include <vector>
9 :
10 : typedef uint16_t adc_t;
11 :
12 : /// Struct to hold information from simulated CAL digis read in from external file
13 : struct calhit
14 : {
15 : int evt; ///< Event Number
16 : int crystalId; ///< Crystal ID
17 : int recoDigiId; ///< Digit ID
18 : int recoDigiT0; ///< TDC Value
19 : int recoDigiSamples; ///< Sample count
20 : std::vector<double> waveform; ///< Reconstructed waveform
21 :
22 : int rocID; ///< ROC ID of CAL digit
23 : int linkID; ///< Link ID of CAL digit
24 : int apdID; ///< APD ID of CAL digit
25 : };
26 :
27 : /// Struct to hold information from simulated TRK digis read in from external file
28 : struct trkhit
29 : {
30 : int evt; ///< Event Number
31 : int strawIdx; ///< Straw Index
32 : int recoDigiT0; ///< TDC Value 0
33 : int recoDigiT1; ///< TDC Value 1
34 : int recoDigiSamples; ///< Sample count
35 : std::vector<double> waveform; ///< Reconstructed waveform
36 :
37 : int rocID; ///< ROC ID of TRK digit
38 : int linkID; ///< Link ID of TRK digit
39 : };
40 :
41 : // Waveform for TRK:
42 : double f(double t, double tau, double sigma, double offset);
43 : // Waveform for CAL
44 : double logn(double x, double eta, double sigma, double Epeak, double norm);
45 : // CAL Digitizer specs:
46 : // 200M samples per second => 5ns steps
47 : // Zero-suppression on board
48 : // 12 bit resolution over 2V dynamic range
49 :
50 : std::vector<adc_t> generateDMABlockHeader(size_t byteCount);
51 :
52 : enum PacketType
53 : {
54 : PacketType_TRK = 1,
55 : PacketType_CAL = 2,
56 : };
57 :
58 0 : unsigned getOptionValue(int* index, char** argv[])
59 : {
60 0 : char* arg = (*argv)[*index];
61 0 : if (arg[2] == '\0')
62 : {
63 0 : (*index)++;
64 0 : return strtoul((*argv)[*index], nullptr, 0);
65 : }
66 0 : int offset = 2;
67 0 : if (arg[2] == '=')
68 : {
69 0 : offset = 3;
70 : }
71 :
72 0 : return strtoul(&(arg[offset]), nullptr, 0);
73 : }
74 :
75 0 : std::string getOptionString(int* index, char** argv[])
76 : {
77 0 : char* arg = (*argv)[*index];
78 0 : if (arg[2] == '\0')
79 : {
80 0 : (*index)++;
81 0 : return std::string((*argv)[*index]);
82 : }
83 0 : int offset = 2;
84 0 : if (arg[2] == '=')
85 : {
86 0 : offset = 3;
87 : }
88 :
89 0 : return std::string(&(arg[offset]));
90 : }
91 :
92 0 : void printHelpMsg()
93 : {
94 0 : std::cout << "Usage: packetGenerator [options]" << std::endl;
95 0 : std::cout << "Options are:" << std::endl
96 0 : << " -h: This message." << std::endl
97 0 : << " -n <number>: Number of timestamps to generate. (Default: 200000)" << std::endl
98 0 : << " -o <number>: Starting Timestamp offest. (Default: 1)." << std::endl
99 0 : << " -v: Verbose mode" << std::endl
100 0 : << " -V: Ridiculously verbose mode" << std::endl
101 0 : << " -f <path>: Filename for output" << std::endl
102 0 : << " -F: Do not output a file" << std::endl
103 0 : << " -S: <number>: PRNG seed value" << std::endl
104 0 : << " -t: Generate Tracker packets (this is the default, default file name is TRK_packets.bin)"
105 0 : << std::endl
106 0 : << " -c: Generate Calorimeter packets (default file name is CAL_packets.bin)" << std::endl;
107 0 : exit(0);
108 : }
109 :
110 0 : int main(int argc, char** argv)
111 : {
112 0 : bool verbose = false;
113 0 : bool veryverbose = false;
114 0 : bool save_adc_values = true;
115 0 : double nevents = 200000; // Number of hits to generate
116 :
117 0 : bool read_cal_digis_from_file = true;
118 0 : std::string inputCalDigiFile = "artdaq_cut_cal.txt";
119 0 : std::ifstream inputCalDigiStream;
120 0 : std::vector<calhit> calHitVector; // Vector of cal hit digi data
121 0 : std::vector<std::vector<calhit> > calEventVector; // Vector of vectors of cal hit digi data for each event
122 : // For now, crystal IDs start at minCrystalID and end at minCrystalID+number_of_links*rocs_per_link-1
123 0 : size_t number_of_crystals_per_roc = 10; // 20 channels, 2 APDs per crystal
124 0 : size_t minCrystalID = 200;
125 :
126 0 : bool read_trk_digis_from_file = true;
127 0 : std::string inputTrkDigiFile = "artdaq_cut_trk.txt";
128 0 : std::ifstream inputTrkDigiStream;
129 0 : std::vector<trkhit> trkHitVector; // Vector of trk hit digi data
130 0 : std::vector<std::vector<trkhit> > trkEventVector; // Vector of vectors of cal hit digi data for each event
131 : // For now, crystal IDs start at minCrystalID and end at minCrystalID+number_of_links*rocs_per_link-1
132 0 : size_t number_of_straws_per_roc = 96; // Each panel in the tracker has 96 straws
133 0 : size_t minStrawID = 960;
134 :
135 0 : size_t SEED = 100; // Random seed
136 :
137 : // Packet type to generate
138 0 : PacketType packetType = PacketType_TRK;
139 :
140 0 : std::string outputFile = "";
141 : // size_t numADCSamples = 8;
142 0 : size_t numADCSamples = 12;
143 :
144 : // The timestamp count will be offset by the following amount:
145 0 : size_t starting_timestamp = 1;
146 :
147 0 : for (int optind = 1; optind < argc; ++optind)
148 : {
149 0 : if (argv[optind][0] == '-')
150 : {
151 0 : switch (argv[optind][1])
152 : {
153 0 : case 'n':
154 0 : nevents = getOptionValue(&optind, &argv);
155 0 : break;
156 0 : case 'o':
157 0 : starting_timestamp = getOptionValue(&optind, &argv);
158 0 : break;
159 0 : case 'v':
160 0 : verbose = true;
161 0 : break;
162 0 : case 'V':
163 0 : verbose = true;
164 0 : veryverbose = true;
165 0 : break;
166 0 : case 'f':
167 0 : save_adc_values = true;
168 0 : outputFile = getOptionString(&optind, &argv);
169 0 : break;
170 0 : case 'F':
171 0 : save_adc_values = false;
172 0 : break;
173 0 : case 't':
174 0 : packetType = PacketType_TRK;
175 0 : break;
176 0 : case 'c':
177 0 : packetType = PacketType_CAL;
178 0 : break;
179 0 : case 'S':
180 0 : SEED = getOptionValue(&optind, &argv);
181 0 : break;
182 0 : default:
183 0 : std::cout << "Unknown option: " << argv[optind] << std::endl;
184 0 : printHelpMsg();
185 0 : break;
186 0 : case 'h':
187 0 : printHelpMsg();
188 0 : break;
189 : }
190 : }
191 : }
192 :
193 0 : size_t max_DMA_block_size = 32000; // Maximum size in bytes of a DMA block
194 : // Normally a DMA block begins when a new timestamp begins, however
195 : // if the size of the DMA block will exceed the limit within the current
196 : // timestamp, a new block is created
197 :
198 0 : if (outputFile == "" && packetType == PacketType_TRK)
199 : {
200 0 : outputFile = "TRK_packets.bin";
201 : }
202 0 : else if (outputFile == "" && packetType == PacketType_CAL)
203 : {
204 0 : outputFile = "CAL_packets.bin";
205 : }
206 0 : std::ofstream binFile;
207 0 : if (save_adc_values)
208 : {
209 0 : binFile.open(outputFile, std::ios::out | std::ios::app | std::ios::binary);
210 : }
211 :
212 0 : int number_of_links = 6;
213 0 : int rocs_per_link = 5;
214 :
215 0 : if (read_cal_digis_from_file && packetType == PacketType_CAL)
216 : {
217 0 : inputCalDigiStream.open(inputCalDigiFile, std::ifstream::in);
218 : int curEvt;
219 0 : while (inputCalDigiStream >> curEvt)
220 : {
221 0 : calhit curHit;
222 0 : curHit.evt = curEvt;
223 0 : inputCalDigiStream >> curHit.crystalId;
224 0 : inputCalDigiStream >> curHit.recoDigiId;
225 0 : inputCalDigiStream >> curHit.recoDigiT0;
226 0 : inputCalDigiStream >> curHit.recoDigiSamples;
227 : double curSample;
228 0 : for (int curSampleNum = 0; curSampleNum < curHit.recoDigiSamples; curSampleNum++)
229 : {
230 0 : inputCalDigiStream >> curSample;
231 0 : curHit.waveform.push_back(curSample);
232 : }
233 :
234 0 : curHit.linkID = int((curHit.crystalId - minCrystalID) / (rocs_per_link * number_of_crystals_per_roc));
235 0 : curHit.rocID =
236 0 : ((curHit.crystalId - minCrystalID) - (rocs_per_link * number_of_crystals_per_roc) * curHit.linkID) /
237 : number_of_crystals_per_roc;
238 0 : curHit.apdID = curHit.recoDigiId % 2; // Even is APD 0, Odd is APD 1
239 :
240 0 : calHitVector.push_back(curHit);
241 0 : }
242 0 : inputCalDigiStream.close();
243 :
244 0 : size_t firstEvtNum = calHitVector[0].evt;
245 0 : size_t lastEvtNum = calHitVector[calHitVector.size() - 1].evt;
246 :
247 0 : std::cout << "Read in " << calHitVector.size() << " calorimeter digi hits" << std::endl;
248 0 : std::cout << "First CAL event num: " << firstEvtNum << std::endl;
249 0 : std::cout << "Last CAL event num: " << lastEvtNum << std::endl;
250 :
251 0 : if (lastEvtNum < nevents)
252 : {
253 : std::cout << "WARNING: nevents is greater than number of events in CAL digi file. "
254 0 : << "Setting nevents to " << lastEvtNum + 1 << std::endl;
255 0 : nevents = lastEvtNum + 1;
256 : }
257 :
258 0 : std::vector<calhit> curEventVector;
259 0 : size_t prevEvent = calHitVector[0].evt;
260 : size_t curEvent;
261 0 : for (size_t curHitVecPos = 0; curHitVecPos < calHitVector.size(); curHitVecPos++)
262 : {
263 0 : curEvent = calHitVector[curHitVecPos].evt;
264 0 : if (curEvent != prevEvent)
265 : {
266 0 : calEventVector.push_back(curEventVector);
267 0 : curEventVector.clear();
268 0 : prevEvent = curEvent;
269 : }
270 0 : curEventVector.push_back(calHitVector[curHitVecPos]);
271 : }
272 0 : calEventVector.push_back(curEventVector);
273 :
274 0 : std::cout << "Length of calEventVector: " << calEventVector.size() << std::endl;
275 0 : }
276 :
277 0 : if (read_trk_digis_from_file && packetType == PacketType_TRK)
278 : {
279 0 : inputTrkDigiStream.open(inputTrkDigiFile, std::ifstream::in);
280 : int curEvt;
281 0 : while (inputTrkDigiStream >> curEvt)
282 : {
283 0 : trkhit curHit;
284 0 : curHit.evt = curEvt;
285 0 : inputTrkDigiStream >> curHit.strawIdx;
286 0 : inputTrkDigiStream >> curHit.recoDigiT0;
287 0 : inputTrkDigiStream >> curHit.recoDigiT1;
288 0 : inputTrkDigiStream >> curHit.recoDigiSamples;
289 :
290 : // The tracker should have a fixed number of digitizer samples (12)
291 0 : assert(static_cast<int>(numADCSamples) == curHit.recoDigiSamples);
292 :
293 : double curSample;
294 0 : for (int curSampleNum = 0; curSampleNum < curHit.recoDigiSamples; curSampleNum++)
295 : {
296 0 : inputTrkDigiStream >> curSample;
297 0 : curHit.waveform.push_back(curSample);
298 : }
299 :
300 0 : curHit.linkID = int((curHit.strawIdx - minStrawID) / (rocs_per_link * number_of_straws_per_roc));
301 0 : curHit.rocID = ((curHit.strawIdx - minStrawID) - (rocs_per_link * number_of_straws_per_roc) * curHit.linkID) /
302 : number_of_straws_per_roc;
303 :
304 0 : trkHitVector.push_back(curHit);
305 0 : }
306 0 : inputTrkDigiStream.close();
307 :
308 0 : size_t firstEvtNum = trkHitVector[0].evt;
309 0 : size_t lastEvtNum = trkHitVector[trkHitVector.size() - 1].evt;
310 :
311 0 : std::cout << "Read in " << trkHitVector.size() << " tracker digi hits" << std::endl;
312 0 : std::cout << "First TRK event num: " << firstEvtNum << std::endl;
313 0 : std::cout << "Last TRK event num: " << lastEvtNum << std::endl;
314 :
315 0 : if (lastEvtNum < nevents)
316 : {
317 : std::cout << "WARNING: nevents is greater than number of events in TRK digi file. "
318 0 : << "Setting nevents to " << lastEvtNum + 1 << std::endl;
319 0 : nevents = lastEvtNum + 1;
320 : }
321 :
322 0 : std::vector<trkhit> curEventVector;
323 0 : size_t prevEvent = trkHitVector[0].evt;
324 : size_t curEvent;
325 0 : for (size_t curHitVecPos = 0; curHitVecPos < trkHitVector.size(); curHitVecPos++)
326 : {
327 0 : curEvent = trkHitVector[curHitVecPos].evt;
328 0 : if (curEvent != prevEvent)
329 : {
330 0 : trkEventVector.push_back(curEventVector);
331 0 : curEventVector.clear();
332 0 : prevEvent = curEvent;
333 : }
334 0 : curEventVector.push_back(trkHitVector[curHitVecPos]);
335 : }
336 0 : trkEventVector.push_back(curEventVector);
337 :
338 0 : std::cout << "Length of trkEventVector: " << trkEventVector.size() << std::endl;
339 0 : }
340 :
341 : double peakVal;
342 : int noise_level;
343 : int sigma_noise;
344 : double threshold;
345 0 : auto nominal_tau = -0.9;
346 : double nominal_sigma;
347 0 : double nominal_offset = 10;
348 0 : auto nominal_eta = 0.03;
349 0 : double nominal_Epeak = 10;
350 0 : double nominal_norm = 1000;
351 0 : double sigma_tau = 2;
352 : double sigma_sigma;
353 0 : double sigma_offset = 2;
354 0 : auto sigma_eta = 0.03;
355 0 : double sigma_Epeak = 10;
356 0 : double sigma_norm = 15;
357 : double minTime;
358 : double maxTime;
359 : double stepSize;
360 :
361 0 : switch (packetType)
362 : {
363 0 : case PacketType_TRK:
364 : // Set maximum waveform value for scaling purposes
365 0 : peakVal = 0.02;
366 :
367 0 : noise_level = 200;
368 0 : sigma_noise = 15;
369 :
370 0 : threshold = 0.0001; // Detection threshold
371 :
372 : // Nominal values for hit generation
373 0 : nominal_tau = 20;
374 0 : nominal_sigma = 20;
375 0 : nominal_offset = 10;
376 : // Uncertainties on nominal values for use in randomizing hits
377 0 : sigma_tau = 2;
378 0 : sigma_sigma = 2;
379 0 : sigma_offset = 2;
380 :
381 0 : minTime = 0;
382 0 : maxTime = 300; // nanoseconds
383 0 : stepSize = 20; // 20ns => 50M samples per second
384 0 : break;
385 0 : case PacketType_CAL:
386 : // Set maximum logn value for scaling purposes
387 0 : peakVal = 32; // Actually about 28.8372
388 0 : noise_level = 140;
389 0 : sigma_noise = 10;
390 :
391 0 : threshold = 0.3; // Detection threshold
392 :
393 : // Nominal values for hit generation
394 0 : nominal_eta = -0.9;
395 0 : nominal_sigma = 55;
396 0 : nominal_Epeak = 200;
397 0 : nominal_norm = 1000;
398 : // Uncertainties on nominal values for use in randomizing hits
399 0 : sigma_eta = 0.03;
400 0 : sigma_sigma = 5;
401 0 : sigma_Epeak = 10;
402 0 : sigma_norm = 15;
403 :
404 0 : minTime = 0;
405 0 : maxTime = 2000; // nanoseconds
406 0 : stepSize = 5; // 5ns => 200M samples per second
407 0 : break;
408 0 : default:
409 0 : std::cout << "ERROR: Invalid packetType: " << packetType << std::endl;
410 0 : return 1;
411 : }
412 :
413 : // PRNG initialization
414 : // Seed with a real random value, if available
415 0 : std::random_device r;
416 0 : std::default_random_engine generator(r());
417 0 : generator.seed(SEED);
418 :
419 0 : std::normal_distribution<double> tau_distribution(nominal_tau, sigma_tau);
420 0 : std::normal_distribution<double> sigma_distribution(nominal_sigma, sigma_sigma);
421 0 : std::normal_distribution<double> offset_distribution(nominal_offset, sigma_offset);
422 :
423 0 : std::uniform_int_distribution<int> strawIndex_distribution(0, 23039);
424 :
425 : // For now, just use a random value between 0 and 50000 for the TDC
426 0 : std::uniform_int_distribution<int> TDC_distribution(0, 50000);
427 :
428 0 : std::normal_distribution<double> eta_distribution(nominal_eta, sigma_eta);
429 0 : std::normal_distribution<double> Epeak_distribution(nominal_Epeak, sigma_Epeak);
430 0 : std::normal_distribution<double> norm_distribution(nominal_norm, sigma_norm);
431 :
432 0 : std::uniform_int_distribution<int> crystalID_distribution(0, 1860);
433 0 : std::uniform_int_distribution<int> apdID_distribution(0, 1);
434 :
435 0 : std::normal_distribution<double> noise_distribution(noise_level, sigma_noise);
436 :
437 : // timeStampVector holds 1 vector per timestamp
438 : // Each of these vectors contains all the datablocks for that timestamp
439 : // Each datablock vector contains the adc_t values for the header packet
440 : // and associated data packets
441 0 : std::vector<std::vector<std::vector<adc_t> > > timeStampVector;
442 :
443 0 : for (size_t eventNum = 0; eventNum < nevents; eventNum++)
444 : {
445 : // Create a vector to hold all Link ID / ROC ID combinations to simulate
446 0 : std::vector<std::pair<int, int> > rocLinkVector;
447 0 : for (auto curLinkID = 0; curLinkID < number_of_links; curLinkID++)
448 : {
449 0 : for (auto curROCID = 0; curROCID < rocs_per_link; curROCID++)
450 : {
451 0 : std::pair<int, int> curPair(curLinkID, curROCID);
452 0 : rocLinkVector.push_back(curPair);
453 : }
454 : }
455 : // // Randomize the order in which the Links and ROCs are received
456 : // std::shuffle(rocLinkVector.begin(),rocLinkVector.end(),generator);
457 :
458 : // Vector to hold all the DataBlocks for this event
459 0 : std::vector<std::vector<adc_t> > curDataBlockVector;
460 :
461 : // Loop over the ROC/link pairs and generate datablocks for each ROC on
462 : // all the links
463 0 : auto targetNumROCs = rocLinkVector.size();
464 0 : for (size_t curPairNum = 0; curPairNum < targetNumROCs; curPairNum++)
465 : {
466 0 : size_t linkID = rocLinkVector[curPairNum].first;
467 0 : size_t rocID = rocLinkVector[curPairNum].second;
468 :
469 0 : if (read_cal_digis_from_file && packetType == PacketType_CAL)
470 : {
471 0 : std::vector<calhit> curHitVector;
472 : // Find all hits for this event coming from the specified Link/ROC
473 0 : for (size_t curHitIdx = 0; curHitIdx < calEventVector[eventNum].size(); curHitIdx++)
474 : {
475 0 : if (calEventVector[eventNum][curHitIdx].rocID == static_cast<int>(rocID) &&
476 0 : calEventVector[eventNum][curHitIdx].linkID == static_cast<int>(linkID))
477 : {
478 0 : curHitVector.push_back(calEventVector[eventNum][curHitIdx]);
479 : }
480 : }
481 :
482 0 : if (curHitVector.size() == 0)
483 : {
484 : // No hits, so just fill a header packet and no data packets
485 0 : std::vector<adc_t> curDataBlock;
486 : // Add the header packet to the DataBlock (leaving including a placeholder for
487 : // the number of packets in the DataBlock);
488 0 : adc_t null_adc = 0;
489 : // First 16 bits of header (reserved values)
490 0 : curDataBlock.push_back(null_adc);
491 : // Second 16 bits of header (ROC ID, packet type, and link ID):
492 0 : std::bitset<16> curROCID = rocID; // 4 bit ROC ID
493 0 : std::bitset<16> headerPacketType = 5; // 4 bit Data packet header type is 5
494 0 : headerPacketType <<= 4; // Shift left by 4
495 0 : std::bitset<16> curLinkID = linkID; // 3 bit link ID
496 0 : curLinkID <<= 8; // Shift left by 8
497 0 : std::bitset<16> secondEntry = (curROCID | headerPacketType | curLinkID);
498 0 : secondEntry[15] = 1; // valid bit
499 0 : curDataBlock.push_back((adc_t)secondEntry.to_ulong());
500 : // Third 16 bits of header (number of data packets is 0)
501 0 : curDataBlock.push_back(null_adc);
502 : // Fourth through sixth 16 bits of header (timestamp)
503 0 : uint64_t timestamp = eventNum + starting_timestamp;
504 0 : curDataBlock.push_back(static_cast<adc_t>(timestamp & 0xFFFF));
505 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 16) & 0xFFFF));
506 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 32) & 0xFFFF));
507 :
508 : // Seventh 16 bits of header (data packet format version and status)
509 0 : adc_t status = 0; // 0 Corresponds to "Timestamp has valid data"
510 0 : adc_t formatVersion = (5 << 8); // Using 5 for now
511 0 : curDataBlock.push_back(formatVersion + status);
512 : // Eighth 16 bits of header (Unassigned)
513 0 : curDataBlock.push_back(null_adc);
514 :
515 : // Fill in the byte count field of the header packet
516 0 : adc_t numBytes = 16; // Just the header packet
517 0 : curDataBlock[0] = numBytes;
518 :
519 0 : curDataBlockVector.push_back(curDataBlock);
520 0 : }
521 : else
522 : {
523 0 : for (size_t curHitIdx = 0; curHitIdx < curHitVector.size(); curHitIdx++)
524 : {
525 : // Generate a DataBlock for every hit on this ROC
526 :
527 0 : auto curHit = curHitVector[curHitIdx];
528 :
529 0 : std::vector<adc_t> curDataBlock;
530 : // Add the header packet to the DataBlock (leaving including a placeholder for
531 : // the number of packets in the DataBlock);
532 0 : adc_t null_adc = 0;
533 : // First 16 bits of header (reserved values)
534 0 : curDataBlock.push_back(null_adc);
535 : // Second 16 bits of header (ROC ID, packet type, and link ID):
536 0 : std::bitset<16> curROCID = rocID; // 4 bit ROC ID
537 0 : std::bitset<16> headerPacketType = 5; // 4 bit Data packet header type is 5
538 0 : headerPacketType <<= 4; // Shift left by 4
539 0 : std::bitset<16> curLinkID = linkID; // 3 bit link ID
540 0 : curLinkID <<= 8; // Shift left by 8
541 0 : auto secondEntry = (curROCID | headerPacketType | curLinkID);
542 0 : secondEntry[15] = 1; // valid bit
543 0 : curDataBlock.push_back(static_cast<adc_t>(secondEntry.to_ulong()));
544 : // Third 16 bits of header (number of data packets is 0)
545 0 : curDataBlock.push_back(null_adc);
546 : // Fourth through sixth 16 bits of header (timestamp)
547 0 : uint64_t timestamp = eventNum + starting_timestamp;
548 0 : curDataBlock.push_back(static_cast<adc_t>(timestamp & 0xFFFF));
549 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 16) & 0xFFFF));
550 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 32) & 0xFFFF));
551 :
552 : // Seventh 16 bits of header (data packet format version and status)
553 0 : adc_t status = 0; // 0 Corresponds to "Timestamp has valid data"
554 0 : adc_t formatVersion = (5 << 8); // Using 5 for now
555 0 : curDataBlock.push_back(formatVersion + status);
556 : // Eighth 16 bits of header (Unassigned)
557 0 : curDataBlock.push_back(null_adc);
558 :
559 : // Create a vector of adc_t values corresponding to
560 : // the content of CAL data packets.
561 0 : std::vector<adc_t> packetVector;
562 :
563 : // Fill the data packets:
564 : // Assume the 0th apd is always read out before the second
565 0 : adc_t crystalID = curHit.crystalId;
566 0 : adc_t apdID = curHit.recoDigiId % 2;
567 0 : adc_t IDNum = ((apdID << 12) | crystalID);
568 :
569 0 : packetVector.push_back(IDNum);
570 0 : packetVector.push_back((adc_t)(curHit.recoDigiT0));
571 0 : packetVector.push_back((adc_t)(curHit.recoDigiSamples));
572 0 : for (auto sampleIdx = 0; sampleIdx < curHit.recoDigiSamples; sampleIdx++)
573 : {
574 0 : adc_t scaledVal = static_cast<adc_t>(curHit.waveform[sampleIdx]);
575 0 : packetVector.push_back(scaledVal);
576 : }
577 : // Pad any empty space in the last packet with 0s
578 0 : size_t padding_slots = 8 - ((curHit.recoDigiSamples - 5) % 8);
579 0 : if (padding_slots < 8)
580 : {
581 0 : for (size_t i = 0; i < padding_slots; i++)
582 : {
583 0 : packetVector.push_back((adc_t)0);
584 : }
585 : }
586 :
587 : // Fill in the number of data packets entry in the header packet
588 0 : adc_t numDataPackets = static_cast<adc_t>(packetVector.size() / 8);
589 0 : curDataBlock[2] = numDataPackets;
590 :
591 : // Fill in the byte count field of the header packet
592 0 : adc_t numBytes = (numDataPackets + 1) * 16;
593 0 : curDataBlock[0] = numBytes;
594 :
595 : // Append the data packets after the header packet in the DataBlock
596 0 : curDataBlock.insert(curDataBlock.end(), packetVector.begin(), packetVector.end());
597 0 : curDataBlockVector.push_back(curDataBlock);
598 0 : }
599 : }
600 0 : }
601 0 : else if (read_trk_digis_from_file && packetType == PacketType_TRK)
602 : {
603 0 : std::vector<trkhit> curHitVector;
604 : // Find all hits for this event coming from the specified Link/ROC
605 0 : for (size_t curHitIdx = 0; curHitIdx < trkEventVector[eventNum].size(); curHitIdx++)
606 : {
607 0 : if (trkEventVector[eventNum][curHitIdx].rocID == (int)rocID &&
608 0 : trkEventVector[eventNum][curHitIdx].linkID == (int)linkID)
609 : {
610 0 : curHitVector.push_back(trkEventVector[eventNum][curHitIdx]);
611 : }
612 : }
613 :
614 0 : if (curHitVector.size() == 0)
615 : {
616 : // No hits, so just fill a header packet and no data packets
617 0 : std::vector<adc_t> curDataBlock;
618 : // Add the header packet to the DataBlock (leaving including a placeholder for
619 : // the number of packets in the DataBlock);
620 0 : adc_t null_adc = 0;
621 : // First 16 bits of header (reserved values)
622 0 : curDataBlock.push_back(null_adc);
623 : // Second 16 bits of header (ROC ID, packet type, and link ID):
624 0 : std::bitset<16> curROCID = rocID; // 4 bit ROC ID
625 0 : std::bitset<16> headerPacketType = 5; // 4 bit Data packet header type is 5
626 0 : headerPacketType <<= 4; // Shift left by 4
627 0 : std::bitset<16> curLinkID = linkID; // 3 bit link ID
628 0 : curLinkID <<= 8; // Shift left by 8
629 0 : std::bitset<16> secondEntry = (curROCID | headerPacketType | curLinkID);
630 0 : secondEntry[15] = 1; // valid bit
631 0 : curDataBlock.push_back((adc_t)secondEntry.to_ulong());
632 : // Third 16 bits of header (number of data packets is 0)
633 0 : curDataBlock.push_back(null_adc);
634 : // Fourth through sixth 16 bits of header (timestamp)
635 0 : uint64_t timestamp = eventNum + starting_timestamp;
636 0 : curDataBlock.push_back(static_cast<adc_t>(timestamp & 0xFFFF));
637 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 16) & 0xFFFF));
638 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 32) & 0xFFFF));
639 :
640 : // Seventh 16 bits of header (data packet format version and status)
641 0 : adc_t status = 0; // 0 Corresponds to "Timestamp has valid data"
642 0 : adc_t formatVersion = (5 << 8); // Using 5 for now
643 0 : curDataBlock.push_back(formatVersion + status);
644 : // Eighth 16 bits of header (Unassigned)
645 0 : curDataBlock.push_back(null_adc);
646 :
647 : // Fill in the byte count field of the header packet
648 0 : adc_t numBytes = 16; // Just the header packet
649 0 : curDataBlock[0] = numBytes;
650 :
651 0 : curDataBlockVector.push_back(curDataBlock);
652 0 : }
653 : else
654 : {
655 0 : for (size_t curHitIdx = 0; curHitIdx < curHitVector.size(); curHitIdx++)
656 : {
657 : // Generate a DataBlock for every hit on this ROC
658 :
659 0 : trkhit curHit = curHitVector[curHitIdx];
660 :
661 0 : std::vector<adc_t> curDataBlock;
662 : // Add the header packet to the DataBlock (leaving including a placeholder for
663 : // the number of packets in the DataBlock);
664 0 : adc_t null_adc = 0;
665 : // First 16 bits of header (reserved values)
666 0 : curDataBlock.push_back(null_adc);
667 : // Second 16 bits of header (ROC ID, packet type, and link ID):
668 0 : std::bitset<16> curROCID = rocID; // 4 bit ROC ID
669 0 : std::bitset<16> headerPacketType = 5; // 4 bit Data packet header type is 5
670 0 : headerPacketType <<= 4; // Shift left by 4
671 0 : std::bitset<16> curLinkID = linkID; // 3 bit link ID
672 0 : curLinkID <<= 8; // Shift left by 8
673 0 : std::bitset<16> secondEntry = (curROCID | headerPacketType | curLinkID);
674 0 : secondEntry[15] = 1; // valid bit
675 0 : curDataBlock.push_back((adc_t)secondEntry.to_ulong());
676 : // Third 16 bits of header (number of data packets is 0)
677 0 : curDataBlock.push_back(null_adc);
678 : // Fourth through sixth 16 bits of header (timestamp)
679 0 : uint64_t timestamp = eventNum + starting_timestamp;
680 0 : curDataBlock.push_back(static_cast<adc_t>(timestamp & 0xFFFF));
681 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 16) & 0xFFFF));
682 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 32) & 0xFFFF));
683 :
684 : // Seventh 16 bits of header (data packet format version and status)
685 0 : adc_t status = 0; // 0 Corresponds to "Timestamp has valid data"
686 0 : adc_t formatVersion = (5 << 8); // Using 5 for now
687 0 : curDataBlock.push_back(formatVersion + status);
688 : // Eighth 16 bits of header (Unassigned)
689 0 : curDataBlock.push_back(null_adc);
690 :
691 : // Create a vector of adc_t values corresponding to
692 : // the content of TRK data packets.
693 0 : std::vector<adc_t> packetVector;
694 :
695 : // Fill the data packets:
696 : // Assume the 0th apd is always read out before the second
697 0 : adc_t strawIndex = curHit.strawIdx;
698 :
699 0 : uint32_t TDC0 = curHit.recoDigiT0;
700 0 : uint32_t TDC1 = curHit.recoDigiT1;
701 :
702 0 : adc_t TDC0_low = TDC0;
703 0 : adc_t TDC0_high = TDC0 >> 16;
704 0 : adc_t TDC1_low = TDC1 << 8;
705 0 : adc_t TDC1_high = TDC1 >> 8;
706 :
707 0 : adc_t TDC0_high_TDC1_low = TDC0_high | TDC1_low;
708 :
709 0 : packetVector.push_back(strawIndex);
710 0 : packetVector.push_back(TDC0_low);
711 0 : packetVector.push_back(TDC0_high_TDC1_low);
712 0 : packetVector.push_back(TDC1_high);
713 :
714 0 : for (int sampleIdx = 0; sampleIdx < curHit.recoDigiSamples; sampleIdx++)
715 : {
716 0 : adc_t scaledVal = static_cast<adc_t>(curHit.waveform[sampleIdx]);
717 0 : packetVector.push_back(scaledVal);
718 : }
719 : // Pad any empty space in the last packet with 0s
720 0 : size_t padding_slots = 8 - ((curHit.recoDigiSamples - 4) % 8);
721 0 : if (padding_slots < 8)
722 : {
723 0 : for (size_t i = 0; i < padding_slots; i++)
724 : {
725 0 : packetVector.push_back((adc_t)0);
726 : }
727 : }
728 :
729 : // Fill in the number of data packets entry in the header packet
730 :
731 0 : adc_t numDataPackets = static_cast<adc_t>(packetVector.size() / 8);
732 0 : curDataBlock[2] = numDataPackets;
733 :
734 : // Fill in the byte count field of the header packet
735 0 : adc_t numBytes = (numDataPackets + 1) * 16;
736 0 : curDataBlock[0] = numBytes;
737 :
738 : // Append the data packets after the header packet in the DataBlock
739 0 : curDataBlock.insert(curDataBlock.end(), packetVector.begin(), packetVector.end());
740 0 : curDataBlockVector.push_back(curDataBlock);
741 0 : }
742 : }
743 0 : }
744 : else
745 : {
746 : // Generate packet content based on approximate waveforms instead of more complete simulation
747 0 : std::vector<adc_t> curDataBlock;
748 : // Add the header packet to the DataBlock (leaving including a placeholder for
749 : // the number of packets in the DataBlock);
750 0 : adc_t null_adc = 0;
751 : // First 16 bits of header (reserved values)
752 0 : curDataBlock.push_back(null_adc);
753 : // Second 16 bits of header (ROC ID, packet type, and link ID):
754 0 : std::bitset<16> curROCID = rocID; // 4 bit ROC ID
755 0 : std::bitset<16> headerPacketType = 5; // 4 bit Data packet header type is 5
756 0 : headerPacketType <<= 4; // Shift left by 4
757 0 : std::bitset<16> curLinkID = linkID; // 3 bit link ID
758 0 : curLinkID <<= 8; // Shift left by 8
759 0 : std::bitset<16> secondEntry = (curROCID | headerPacketType | curLinkID);
760 0 : secondEntry[15] = 1; // valid bit
761 0 : curDataBlock.push_back((adc_t)secondEntry.to_ulong());
762 : // Third 16 bits of header (number of data packets is 0)
763 0 : curDataBlock.push_back(null_adc);
764 : // Fourth through sixth 16 bits of header (timestamp)
765 0 : uint64_t timestamp = eventNum + starting_timestamp;
766 0 : curDataBlock.push_back(static_cast<adc_t>(timestamp & 0xFFFF));
767 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 16) & 0xFFFF));
768 0 : curDataBlock.push_back(static_cast<adc_t>((timestamp >> 32) & 0xFFFF));
769 :
770 : // Seventh 16 bits of header (data packet format version and status)
771 0 : adc_t status = 0; // 0 Corresponds to "Timestamp has valid data"
772 0 : adc_t formatVersion = (5 << 8); // Using 5 for now
773 0 : curDataBlock.push_back(formatVersion + status);
774 : // Eighth 16 bits of header (Unassigned)
775 0 : curDataBlock.push_back(null_adc);
776 :
777 : // Vector to hold raw digitized waveform
778 0 : std::vector<double> digiVector;
779 0 : double hitTime = 0.0;
780 :
781 : // Create a vector of adc_t values corresponding to
782 : // the content of TRK/CAL data packets.
783 0 : std::vector<adc_t> packetVector;
784 :
785 : // Sample distributions
786 0 : double cur_tau = tau_distribution(generator);
787 0 : double cur_sigma = sigma_distribution(generator);
788 0 : double cur_offset = offset_distribution(generator);
789 :
790 0 : double cur_eta = eta_distribution(generator);
791 0 : double cur_Epeak = Epeak_distribution(generator);
792 0 : double cur_norm = norm_distribution(generator);
793 :
794 : // Control flags
795 0 : bool inWindow = false;
796 0 : bool exitWindow = false;
797 0 : double curVal = 0;
798 : size_t padding_slots;
799 :
800 0 : switch (packetType)
801 : {
802 0 : case PacketType_TRK: {
803 : // Generate TRK data packets
804 :
805 0 : for (double curTime = minTime; curTime < maxTime && !exitWindow; curTime += stepSize)
806 : {
807 0 : curVal = f(curTime, cur_tau, cur_sigma, cur_offset);
808 0 : if (!inWindow && curVal > threshold)
809 : {
810 0 : inWindow = true;
811 0 : curTime -= 2 * stepSize; // Save samples starting 2 steps before the threshold is reached
812 0 : hitTime = curTime;
813 0 : curVal = f(curTime, cur_tau, cur_sigma, cur_offset);
814 : }
815 0 : else if (curVal < threshold && inWindow && curTime > hitTime + 2 * stepSize)
816 : {
817 0 : inWindow = false;
818 0 : exitWindow = true;
819 : }
820 :
821 0 : if (inWindow)
822 : {
823 0 : digiVector.push_back(curVal);
824 : }
825 : }
826 :
827 0 : adc_t strawIndex = strawIndex_distribution(generator);
828 0 : uint32_t TDC0 = TDC_distribution(generator);
829 0 : uint32_t TDC1 = TDC_distribution(generator);
830 :
831 0 : adc_t TDC0_low = TDC0;
832 0 : adc_t TDC0_high = TDC0 >> 16;
833 0 : adc_t TDC1_low = TDC1 << 8;
834 0 : adc_t TDC1_high = TDC1 >> 8;
835 :
836 0 : adc_t TDC0_high_TDC1_low = TDC0_high | TDC1_low;
837 :
838 0 : packetVector.push_back(strawIndex);
839 0 : packetVector.push_back(TDC0_low);
840 0 : packetVector.push_back(TDC0_high_TDC1_low);
841 0 : packetVector.push_back(TDC1_high);
842 :
843 0 : for (size_t i = 0; i < numADCSamples; i++)
844 : {
845 0 : adc_t scaledVal = 0;
846 : // Scale the function value relative to a peak
847 : // value of around 0.02 and convert to a 12 bit integer
848 : // stored in a 16 bit adc_t
849 0 : if (i < digiVector.size())
850 : {
851 0 : scaledVal = static_cast<adc_t>(digiVector[i] / (1.0 * peakVal) * (1 << 12));
852 : }
853 : // Add some noise
854 0 : adc_t noise = static_cast<adc_t>(noise_distribution(generator));
855 0 : packetVector.push_back(scaledVal + noise);
856 : }
857 : // Pad any empty space in the last packet with 0s
858 0 : padding_slots = 8 - ((numADCSamples - 4) % 8);
859 0 : if (padding_slots < 8)
860 : {
861 0 : for (size_t i = 0; i < padding_slots; i++)
862 : {
863 0 : packetVector.push_back((adc_t)0);
864 : }
865 : }
866 : }
867 0 : break;
868 0 : case PacketType_CAL: {
869 : {
870 : // Approximate CAL waveforms on the fly instead of reading them in from a file
871 0 : for (double curTime = minTime; curTime < maxTime && !exitWindow; curTime += stepSize)
872 : {
873 0 : curVal = logn(curTime, cur_eta, cur_sigma, cur_Epeak, cur_norm);
874 0 : if (!inWindow && curVal > threshold)
875 : {
876 0 : inWindow = true;
877 0 : curTime -= 5 * stepSize; // Save samples starting 5 steps before the threshold is reached
878 0 : hitTime = curTime;
879 0 : curVal = logn(curTime, cur_eta, cur_sigma, cur_Epeak, cur_norm);
880 : }
881 0 : else if (curVal < threshold && inWindow && curTime > hitTime + 5 * stepSize)
882 : {
883 0 : inWindow = false;
884 0 : exitWindow = true;
885 : }
886 :
887 0 : if (inWindow)
888 : {
889 0 : digiVector.push_back(curVal);
890 : }
891 : }
892 :
893 0 : adc_t apdID = apdID_distribution(generator);
894 0 : adc_t crystalID = crystalID_distribution(generator);
895 0 : adc_t IDNum = ((apdID << 12) | crystalID);
896 0 : packetVector.push_back(IDNum);
897 0 : packetVector.push_back((adc_t)hitTime);
898 0 : packetVector.push_back((adc_t)digiVector.size());
899 0 : for (size_t i = 0; i < digiVector.size(); i++)
900 : {
901 : // Scale the function value relative to a peak possible
902 : // value of around 35 and convert to a 12 bit integer
903 : // stored in a 16 bit adc_t
904 0 : adc_t scaledVal = static_cast<adc_t>(digiVector[i] / (1.0 * peakVal) * (1 << 12));
905 : // Add some noise
906 0 : adc_t noise = static_cast<adc_t>(noise_distribution(generator));
907 0 : packetVector.push_back(scaledVal + noise);
908 : }
909 : // Pad any empty space in the last packet with 0s
910 0 : padding_slots = 8 - ((digiVector.size() - 5) % 8);
911 0 : if (padding_slots < 8)
912 : {
913 0 : for (size_t i = 0; i < padding_slots; i++)
914 : {
915 0 : packetVector.push_back((adc_t)0);
916 : }
917 : }
918 : }
919 : }
920 0 : break;
921 0 : default:
922 0 : break;
923 : }
924 :
925 : // Fill in the number of data packets entry in the header packet
926 0 : auto numDataPackets = static_cast<adc_t>(packetVector.size() / 8);
927 0 : curDataBlock[2] = numDataPackets;
928 :
929 : // Fill in the byte count field of the header packet
930 0 : adc_t numBytes = (numDataPackets + 1) * 16;
931 0 : curDataBlock[0] = numBytes;
932 :
933 : // Append the data packets after the header packet in the DataBlock
934 0 : curDataBlock.insert(curDataBlock.end(), packetVector.begin(), packetVector.end());
935 0 : curDataBlockVector.push_back(curDataBlock);
936 0 : }
937 :
938 : } // Done generating DataBlocks for this timestamp
939 0 : timeStampVector.push_back(curDataBlockVector);
940 :
941 0 : } // Done generating data for all timestamps
942 :
943 0 : std::cout << "Size of timestamp vector: " << timeStampVector.size() << std::endl;
944 :
945 0 : if (verbose)
946 : {
947 0 : std::cout << "\tNumber of DataBlocks for each timestamp vector: {";
948 0 : for (size_t i = 0; i < timeStampVector.size(); i++)
949 : {
950 0 : std::cout << timeStampVector[i].size();
951 0 : if (i < timeStampVector.size() - 1)
952 : {
953 0 : std::cout << ", ";
954 : }
955 : }
956 0 : std::cout << "}" << std::endl;
957 :
958 0 : std::cout << "\tNumber of packets in each data block: {";
959 0 : for (size_t i = 0; i < timeStampVector.size(); i++)
960 : {
961 0 : std::cout << "{";
962 0 : for (size_t j = 0; j < timeStampVector[i].size(); j++)
963 : {
964 0 : std::cout << (timeStampVector[i][j].size()) / 8;
965 0 : if (j < timeStampVector[i].size() - 1)
966 : {
967 0 : std::cout << ", ";
968 : }
969 : }
970 0 : std::cout << "}";
971 0 : if (i < timeStampVector.size() - 1)
972 : {
973 0 : std::cout << ",";
974 : }
975 : }
976 0 : std::cout << "}" << std::endl;
977 : }
978 :
979 0 : std::vector<adc_t> masterVector;
980 0 : for (size_t i = 0; i < timeStampVector.size(); i++)
981 : {
982 : // Determine how to divide DataBlocks between DMABlocks within each timestamp
983 0 : std::vector<size_t> dataBlockPartition;
984 0 : std::vector<size_t> dataBlockPartitionSizes;
985 0 : for (size_t j = 0, curDMABlockSize = 0, numDataBlocksInCurDMABlock = 0; j < timeStampVector[i].size(); j++)
986 : {
987 0 : numDataBlocksInCurDMABlock++; // Increment number of DataBlocks in the current DMA block
988 0 : curDMABlockSize += timeStampVector[i][j].size() * 2; // Size of current data block in 8bit words
989 0 : assert(curDMABlockSize <= max_DMA_block_size - 8);
990 0 : if (j == timeStampVector[i].size() - 1)
991 : {
992 0 : dataBlockPartition.push_back(numDataBlocksInCurDMABlock);
993 0 : dataBlockPartitionSizes.push_back(curDMABlockSize + 8);
994 : }
995 0 : else if (curDMABlockSize + (2 * timeStampVector[i][j + 1].size()) > max_DMA_block_size - 8)
996 : {
997 0 : dataBlockPartition.push_back(numDataBlocksInCurDMABlock);
998 0 : dataBlockPartitionSizes.push_back(curDMABlockSize + 8);
999 0 : curDMABlockSize = 0;
1000 0 : numDataBlocksInCurDMABlock = 0;
1001 : }
1002 : }
1003 :
1004 : // Break the DataBlocks into DMABlocks and add DMABlock headers
1005 0 : for (size_t curDMABlockNum = 0, curDataBlockNum = 0; curDMABlockNum < dataBlockPartition.size(); curDMABlockNum++)
1006 : {
1007 0 : size_t numDataBlocksInCurDMABlock = dataBlockPartition[curDMABlockNum];
1008 0 : size_t curDMABlockSize = dataBlockPartitionSizes[curDMABlockNum];
1009 0 : std::vector<adc_t> header = generateDMABlockHeader(curDMABlockSize);
1010 0 : for (size_t adcNum = 0; adcNum < header.size(); adcNum++)
1011 : {
1012 0 : masterVector.push_back(header[adcNum]);
1013 : }
1014 :
1015 0 : for (size_t j = 0; j < numDataBlocksInCurDMABlock; j++)
1016 : {
1017 0 : std::vector<adc_t> curDataBlock = timeStampVector[i][curDataBlockNum];
1018 0 : for (size_t adcNum = 0; adcNum < curDataBlock.size(); adcNum++)
1019 : {
1020 0 : masterVector.push_back(curDataBlock[adcNum]);
1021 : }
1022 0 : curDataBlockNum++;
1023 0 : }
1024 0 : }
1025 :
1026 0 : if (verbose)
1027 : {
1028 : // Print out number of DataBlocks in each DMABlock
1029 0 : std::cout << "\tTimestamp " << i + starting_timestamp << " DataBlock partition: {";
1030 0 : for (size_t k = 0; k < dataBlockPartition.size(); k++)
1031 : {
1032 0 : std::cout << dataBlockPartition[k];
1033 0 : if (k < dataBlockPartition.size() - 1)
1034 : {
1035 0 : std::cout << ", ";
1036 : }
1037 : }
1038 0 : std::cout << "}" << std::endl;
1039 : }
1040 0 : } // Close loop over timestamps
1041 :
1042 0 : std::cout << std::endl
1043 0 : << "Length of final adc_t array: " << masterVector.size() << std::endl;
1044 :
1045 : // Print contents of final adc_t array:
1046 0 : if (veryverbose)
1047 : {
1048 0 : std::cout << "Contents of final adc_t array: " << std::endl;
1049 0 : for (size_t i = 0; i < masterVector.size(); i++)
1050 : {
1051 0 : std::bitset<16> curEntry = masterVector[i];
1052 0 : std::cout << "\t" << curEntry.to_string() << " " << curEntry.to_ulong() << std::endl;
1053 0 : if (i > 0 && (i + 1) % 4 == 0)
1054 : {
1055 0 : std::cout << std::endl;
1056 : }
1057 : }
1058 : }
1059 :
1060 : // Save contents of final adc_t array to binary file
1061 0 : std::cout << "Writing generated data to file " << outputFile << std::endl;
1062 0 : if (save_adc_values)
1063 : {
1064 0 : for (size_t i = 0; i < masterVector.size(); i++)
1065 : {
1066 0 : binFile.write(reinterpret_cast<const char*>(&(masterVector[i])), sizeof(adc_t));
1067 : }
1068 0 : binFile.close();
1069 : }
1070 :
1071 0 : return 0;
1072 0 : }
1073 :
1074 0 : double f(double t, double tau, double sigma, double offset)
1075 : {
1076 0 : double E = 2.718281828459045;
1077 0 : double retval = (pow(E, (pow(sigma, 2) + 2 * offset * tau - 2 * t * tau) / (2. * pow(tau, 2))) *
1078 0 : (-pow(sigma, 2) + (-offset + t) * tau)) /
1079 0 : pow(tau, 3);
1080 0 : if (retval < 0)
1081 : {
1082 0 : retval = 0;
1083 : }
1084 0 : return retval;
1085 : }
1086 :
1087 0 : double logn(double x, double eta, double sigma, double Epeak, double norm)
1088 : {
1089 : double Aterm;
1090 : double logterms0, s0;
1091 : double logn, logterm;
1092 : double expterm;
1093 0 : double pigreco = 3.14159265;
1094 :
1095 : // double f = 2.35;
1096 0 : double f = 20.0;
1097 :
1098 0 : logterms0 = eta * f / 2 + sqrt(1 + pow((eta * f / 2), 2));
1099 0 : s0 = (2 / f) * log(logterms0);
1100 :
1101 0 : Aterm = eta / (sqrt(2 * pigreco) * sigma * s0);
1102 :
1103 0 : logterm = 1 - (eta / sigma) * (x - Epeak);
1104 :
1105 0 : if (logterm < 0)
1106 : {
1107 0 : logterm = 0.0001;
1108 : }
1109 0 : expterm = log(logterm) / s0;
1110 0 : expterm = -0.5 * pow(expterm, 2);
1111 :
1112 0 : logn = norm * Aterm * exp(expterm);
1113 0 : return logn;
1114 : }
1115 :
1116 0 : std::vector<adc_t> generateDMABlockHeader(size_t theCount)
1117 : {
1118 0 : std::bitset<64> byteCount = theCount;
1119 0 : std::bitset<16> byteCount0 = 0;
1120 0 : std::bitset<16> byteCount1 = 0;
1121 0 : std::bitset<16> byteCount2 = 0;
1122 0 : std::bitset<16> byteCount3 = 0;
1123 0 : for (int i = 0; i < 16; i++)
1124 : {
1125 0 : byteCount0[i] = byteCount[i + 16 * 0];
1126 0 : byteCount1[i] = byteCount[i + 16 * 1];
1127 0 : byteCount2[i] = byteCount[i + 16 * 2];
1128 0 : byteCount3[i] = byteCount[i + 16 * 3];
1129 : }
1130 0 : std::vector<adc_t> header;
1131 0 : header.push_back((adc_t)byteCount0.to_ulong());
1132 0 : header.push_back((adc_t)byteCount1.to_ulong());
1133 0 : header.push_back((adc_t)byteCount2.to_ulong());
1134 0 : header.push_back((adc_t)byteCount3.to_ulong());
1135 :
1136 0 : return header;
1137 0 : }
|
