BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
BesTofDigitizerEcV4.cc
Go to the documentation of this file.
1//---------------------------------------------------------------------------//
2// BOOST --- BESIII Object_Oriented Simulation Tool //
3//---------------------------------------------------------------------------//
4// Description: This is the Digitizer for the MRPC with doubled sided readout
5// Author: An Fenfen
6// Created: Nov, 2015
7
8//---------------------------------------------------------------------------//
9//$Id: BesTofDigitizerEcV4.cc
10
11#include "Gaudi/Interfaces/IOptionsSvc.h"
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/IDataProviderSvc.h"
14#include "GaudiKernel/ISvcLocator.h"
15#include "GaudiKernel/MsgStream.h"
16#include "GaudiKernel/Parsers.h"
17#include "GaudiKernel/PropertyMgr.h"
18
19#include "G4DigiManager.hh"
20#include "G4PhysicalConstants.hh"
21#include "G4SystemOfUnits.hh"
22#include "Randomize.hh"
23#include "TFile.h"
24#include "TH1F.h"
25#include "TMath.h"
26#include "TNtuple.h"
27#include "TROOT.h"
28#include "TSpline.h"
29#include "TofSim/BesTofDigi.hh"
30#include "TofSim/BesTofDigitizerEcV4.hh"
31#include "TofSim/BesTofHit.hh"
32#include <math.h>
33
34using Gaudi::Interfaces::IOptionsSvc;
35
36BesTofDigitizerEcV4::BesTofDigitizerEcV4()
37 : tdcRes_const( "tdcRes_const", -1 )
38 , adcRes_const( "adcRes_const", -1 )
39 , m_rootFlag( "RootFlag", false )
40 , m_fileName( "FileName", "mrpc.root" )
41 , m_V( "E", 7000 )
42 , m_threshold( "Threshold", 5.5e+08 )
43 , m_nstep( "nstep", -1 )
44 , m_E_weight( "E_weight", -1 )
45 , m_saturationFlag( "saturationFlag", true )
46 , m_calTdcRes_charge_flag( "calTdcRes_charge_flag", 0 )
47 , m_charge2Time_flag( "charge2Time_flag", 0 )
48 , m_calAdcRes_charge_flag( "calAdcRes_charge_flag", 0 ) {
49
50 // declare properties
51 auto opt_svc = Gaudi::svcLocator()->service<IOptionsSvc>( "JobOptionsSvc" );
52 opt_svc->bind( "BesTofDigitizerEcV4", &tdcRes_const );
53 opt_svc->bind( "BesTofDigitizerEcV4", &adcRes_const );
54 opt_svc->bind( "BesTofDigitizerEcV4", &m_rootFlag );
55 opt_svc->bind( "BesTofDigitizerEcV4", &m_fileName );
56 opt_svc->bind( "BesTofDigitizerEcV4", &m_V );
57 opt_svc->bind( "BesTofDigitizerEcV4", &m_threshold );
58 opt_svc->bind( "BesTofDigitizerEcV4", &m_nstep );
59 opt_svc->bind( "BesTofDigitizerEcV4", &m_E_weight );
60 opt_svc->bind( "BesTofDigitizerEcV4", &m_saturationFlag );
61 opt_svc->bind( "BesTofDigitizerEcV4", &m_calTdcRes_charge_flag );
62 opt_svc->bind( "BesTofDigitizerEcV4", &m_charge2Time_flag );
63 opt_svc->bind( "BesTofDigitizerEcV4", &m_calAdcRes_charge_flag );
64
65 initial();
66
67 if ( m_rootFlag )
68 {
69 m_file = new TFile( m_fileName.value().c_str(), "RECREATE" );
70 m_tree = new TTree( "mrpc", "mrpc" );
71
72 m_tree->Branch( "event", &m_event, "event/D" );
73 m_tree->Branch( "partId", &m_partId, "partId/D" );
74 m_tree->Branch( "module", &m_module, "module/D" );
75 m_tree->Branch( "time_leading_sphi", &m_time_leading_sphi, "time_leading_sphi/D" );
76 m_tree->Branch( "time_leading_xphi", &m_time_leading_xphi, "time_leading_xphi/D" );
77 m_tree->Branch( "time_trailing_sphi", &m_time_trailing_sphi, "time_trailing_sphi/D" );
78 m_tree->Branch( "time_trailing_xphi", &m_time_trailing_xphi, "time_trailing_xphi/D" );
79 m_tree->Branch( "tdcRes", &m_tdcRes, "tdcRes/D" );
80 m_tree->Branch( "tdcRes_charge", &m_tdcRes_charge, "tdcRes_charge/D" );
81 m_tree->Branch( "adc", &m_adc, "adc/D" );
82 m_tree->Branch( "adcRes", &m_adcRes, "adcRes/D" );
83 m_tree->Branch( "adcRes_charge", &m_adcRes_charge, "adcRes_charge/D" );
84 m_tree->Branch( "strip", &m_strip, "strip/D" );
85 m_tree->Branch( "trkIndex", &m_trkIndex, "trkIndex/D" );
86 m_tree->Branch( "tStart", &m_tStart, "tStart/D" );
87 m_tree->Branch( "tPropagate_sphi", &m_tPropagate_sphi, "tPropagate_sphi/D" );
88 m_tree->Branch( "tPropagate_xphi", &m_tPropagate_xphi, "tPropagate_xphi/D" );
89 m_tree->Branch( "tThreshold", &m_tThreshold, "tThreshold/D" );
90 m_tree->Branch( "charge", &m_charge, "charge/D" );
91 m_tree->Branch( "nhit", &m_nhit, "nhit/I" );
92 m_tree->Branch( "ions_hit", m_ions_hit, "ions_hit[nhit]/D" );
93 m_tree->Branch( "trkIndex_hit", m_trkIndex_hit, "trkIndex_hit[nhit]/D" );
94 m_tree->Branch( "pdgCode_hit", m_pdgCode_hit, "pdgCode_hit[nhit]/D" );
95 m_tree->Branch( "gap_hit", m_gap_hit, "gap_hit[nhit]/D" );
96 m_tree->Branch( "underStrip_hit", m_underStrip_hit, "underStrip_hit[nhit]/D" );
97 m_tree->Branch( "locx_hit", m_locx_hit, "locx_hit[nhit]/D" );
98 m_tree->Branch( "locy_hit", m_locy_hit, "locy_hit[nhit]/D" );
99 m_tree->Branch( "locz_hit", m_locz_hit, "locz_hit[nhit]/D" );
100 m_tree->Branch( "x_hit", m_x_hit, "x_hit[nhit]/D" );
101 m_tree->Branch( "y_hit", m_y_hit, "y_hit[nhit]/D" );
102 m_tree->Branch( "z_hit", m_z_hit, "z_hit[nhit]/D" );
103 m_tree->Branch( "px_hit", m_px_hit, "px_hit[nhit]/D" );
104 m_tree->Branch( "py_hit", m_py_hit, "py_hit[nhit]/D" );
105 m_tree->Branch( "pz_hit", m_pz_hit, "pz_hit[nhit]/D" );
106 }
107}
108
109BesTofDigitizerEcV4::~BesTofDigitizerEcV4() {
110 if ( m_rootFlag )
111 {
112 m_file->Write();
113 m_file->Close();
114 }
115}
116
117void BesTofDigitizerEcV4::initial() {
118 m_param = Param();
119 m_param.setPar( m_nstep, m_E_weight, m_V );
120 m_param.print();
121 m_event = -1;
122 partId = -999;
123 module = -999;
124 tdc_sphi = -999;
125 tdc_xphi = -999;
126 if ( tdcRes_const < 0 ) tdcRes_const = 38;
127 // 45; //sqrt(27*27.+30.*30+20.*20); //ps, 27:TDC; 30:gapNb; 20:cables..
128 tdcRes = -999;
129
130 adc = -999;
131 if ( adcRes_const < 0 ) adcRes_const = 27; // ps TDC
132 adcRes = -999;
133
134 time_leading_sphi = -999;
135 time_leading_xphi = -999;
136 time_trailing_sphi = -999;
137 time_trailing_xphi = -999;
138
139 cout << "Property:" << endl
140 << " FileName= " << m_fileName << " E= " << m_V << " Threshold= " << m_threshold
141 << " nstep= " << m_nstep << " E_weight= " << m_E_weight
142 << " saturationFlag= " << m_saturationFlag << " tdcRes_const= " << tdcRes_const
143 << " adcRes_const= " << adcRes_const
144 << " calTdcRes_charge_flag= " << m_calTdcRes_charge_flag
145 << " charge2Time_flag= " << m_charge2Time_flag
146 << " calAdcRes_charge_flag= " << m_calAdcRes_charge_flag << endl;
147}
148
149void BesTofDigitizerEcV4::Digitize( ScintSingle* single_module, BesTofDigitsCollection* DC ) {
150 m_besTofDigitsCollection = DC;
151 G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
152 G4int THCID = digiManager->GetHitsCollectionID( "BesTofHitsCollection" );
153 m_THC = (BesTofHitsCollection*)( digiManager->GetHitsCollection( THCID ) );
154 if ( !m_THC ) return;
155
156 partId = single_module->GetPartId();
157 module = single_module->GetModule_mrpc();
158 // cout<<"partId= "<<partId<<" module= "<<module<<endl;
159
160 // Process the hits
161 int nstrip = m_param.nstrip;
162 StripStruct stripStruct[12];
163 for ( int i = 0; i < nstrip; i++ )
164 {
165 stripStruct[i].m_param = m_param;
166 stripStruct[i].setPar( m_param.alpha, m_param.eta, m_param.v_drift, m_threshold,
167 m_saturationFlag );
168 }
169
170 BesTofHit* hit;
171 for ( unsigned int i = 0; i < single_module->GetHitIndexes_mrpc()->size(); i++ )
172 {
173 hit = ( *m_THC )[( *( single_module->GetHitIndexes_mrpc() ) )[i]];
174 m_event = hit->GetEvent();
175
176 HitStruct hitStruct;
177 hitStruct.m_param = m_param;
178 hitStruct.trkIndex = hit->GetG4Index();
179 hitStruct.pdgCode = hit->GetPDGcode();
180 hitStruct.ions = hit->GetIons();
181 hitStruct.strip = calStrip( hit->GetLocPos().z() / mm );
182 hitStruct.underStrip = underStrip( hit->GetLocPos().x() / mm, hit->GetLocPos().z() / mm );
183 hitStruct.gap = hit->GetGapNb();
184 hitStruct.glbTime = hit->GetTime() / ns;
185 hitStruct.locx = hit->GetLocPos().x() / mm;
186 hitStruct.locy = hit->GetLocPos().y() / mm;
187 hitStruct.locz = hit->GetLocPos().z() / mm;
188 hitStruct.x = hit->GetPos().x() / mm;
189 hitStruct.y = hit->GetPos().y() / mm;
190 hitStruct.z = hit->GetPos().z() / mm;
191 hitStruct.px = hit->GetMomentum().x() / ( GeV / ( 3e+08 * m / s ) );
192 hitStruct.py = hit->GetMomentum().y() / ( GeV / ( 3e+08 * m / s ) );
193 hitStruct.pz = hit->GetMomentum().z() / ( GeV / ( 3e+08 * m / s ) );
194 // hitStruct.print();
195
196 if ( hitStruct.ions > 0 && hitStruct.glbTime > 0 )
197 { stripStruct[hitStruct.strip].hitStructCol.push_back( hitStruct ); }
198 }
199
200 // test multihit, study the lower peak in charge
201 int hitSize = 0;
202
203 for ( int i = 0; i < nstrip; i++ )
204 {
205 if ( stripStruct[i].hitStructCol.size() == 0 ) continue;
206 stripStruct[i].strip = i;
207 stripStruct[i].calFirstHit();
208 stripStruct[i].avalanche();
209 // stripStruct[i].print();
210
211 if ( stripStruct[i].tThreshold <= 0 ) continue;
212
213 tdc_sphi =
214 stripStruct[i].tStart + stripStruct[i].tThreshold + stripStruct[i].tPropagate_sphi;
215 tdc_xphi =
216 stripStruct[i].tStart + stripStruct[i].tThreshold + stripStruct[i].tPropagate_xphi;
217
218 double tdcRes_charge = 0;
219 if ( m_calTdcRes_charge_flag == 0 )
220 {
221 tdcRes_charge = calTdcRes_charge( stripStruct[i].charge * 1000 ); // ps, charge in fC
222 }
223 else if ( m_calTdcRes_charge_flag == 1 )
224 { tdcRes_charge = calTdcRes_charge1( stripStruct[i].charge * 1000 ); }
225 else if ( m_calTdcRes_charge_flag == 2 ) { tdcRes_charge = 0; }
226
227 tdcRes = sqrt( tdcRes_charge * tdcRes_charge + tdcRes_const * tdcRes_const ) / 1000; // ns
228
229 tdc_sphi = G4RandGauss::shoot( tdc_sphi, tdcRes );
230 tdc_xphi = G4RandGauss::shoot( tdc_xphi, tdcRes );
231
232 if ( m_charge2Time_flag == 0 )
233 {
234 adc = charge2Time( stripStruct[i].charge * 1000 ); // ns, charge in fC
235 }
236 else if ( m_charge2Time_flag == 1 ) { adc = charge2Time1( stripStruct[i].charge * 1000 ); }
237
238 double adcRes_charge = 0;
239 if ( m_calAdcRes_charge_flag == 0 )
240 {
241 adcRes_charge = calAdcRes_charge( stripStruct[i].charge * 1000 ); // ps, charge in fC
242 }
243 else if ( m_calAdcRes_charge_flag == 1 )
244 { adcRes_charge = calAdcRes_charge1( stripStruct[i].charge * 1000 ); }
245 else if ( m_calAdcRes_charge_flag == 2 ) { adcRes_charge = 0; }
246
247 adcRes = sqrt( adcRes_charge * adcRes_charge + adcRes_const * adcRes_const ) / 1000;
248 adc = G4RandGauss::shoot( adc, adcRes );
249 if ( adc < 0 ) adc = 0;
250
251 time_leading_sphi = tdc_sphi;
252 time_leading_xphi = tdc_xphi;
253 time_trailing_sphi = tdc_sphi + adc;
254 time_trailing_xphi = tdc_xphi + adc;
255
256 // save digi information
257 BesTofDigi* digi = new BesTofDigi;
258 digi->SetTrackIndex( stripStruct[i].trkIndex );
259 digi->SetPartId( partId );
260 digi->SetModule( module );
261 digi->SetStrip( stripStruct[i].strip );
262 int mo = ( partId - 3 ) * 36 + module;
263 int st = stripStruct[i].strip;
264 if ( m_param.deadChannel[mo][st] == 0 || m_param.deadChannel[mo][st] == 2 )
265 {
266 // Set dead channel
267 digi->SetForwT1( -999 );
268 digi->SetForwT2( -999 );
269 }
270 else
271 {
272 digi->SetForwT1( time_leading_sphi );
273 digi->SetForwT2( time_trailing_sphi );
274 }
275 if ( m_param.deadChannel[mo][st] == 1 || m_param.deadChannel[mo][st] == 2 )
276 {
277 digi->SetBackT1( -999 );
278 digi->SetBackT2( -999 );
279 }
280 else
281 {
282 digi->SetBackT1( time_leading_xphi );
283 digi->SetBackT2( time_trailing_xphi );
284 }
285 m_besTofDigitsCollection->insert( digi );
286 // cout<<"Print digi info: "
287 // <<" partId= "<<partId
288 // <<" module= "<<module
289 // <<" strip= "<<stripStruct[i].strip
290 // <<" time_leading_sphi= "<<time_leading_sphi
291 // <<" time_leading_xphi= "<<time_leading_xphi
292 // <<" time_trailing_sphi= "<<time_trailing_sphi
293 // <<" time_trailing_xphi= "<<time_trailing_xphi
294 // <<endl;
295
296 // save digi information
297 if ( m_rootFlag )
298 {
299 m_partId = partId;
300 m_module = module;
301 m_time_leading_sphi = time_leading_sphi;
302 m_time_leading_xphi = time_leading_xphi;
303 m_time_trailing_sphi = time_trailing_sphi;
304 m_time_trailing_xphi = time_trailing_xphi;
305 m_tdcRes = tdcRes;
306 m_tdcRes_charge = tdcRes_charge;
307 m_adc = adc;
308 m_adcRes = adcRes;
309 m_adcRes_charge = adcRes_charge;
310
311 m_strip = stripStruct[i].strip;
312 m_trkIndex = stripStruct[i].trkIndex;
313 m_tStart = stripStruct[i].tStart;
314 m_tPropagate_sphi = stripStruct[i].tPropagate_sphi;
315 m_tPropagate_xphi = stripStruct[i].tPropagate_xphi;
316 m_tThreshold = stripStruct[i].tThreshold;
317 m_charge = stripStruct[i].charge;
318
319 m_nhit = stripStruct[i].hitStructCol.size();
320 // cout<<"m_nhit= "<<m_nhit<<endl;
321 for ( int j = 0; j < m_nhit; j++ )
322 {
323 m_ions_hit[j] = stripStruct[i].hitStructCol[j].ions;
324 m_trkIndex_hit[j] = stripStruct[i].hitStructCol[j].trkIndex;
325 m_pdgCode_hit[j] = stripStruct[i].hitStructCol[j].pdgCode;
326 m_gap_hit[j] = stripStruct[i].hitStructCol[j].gap;
327 m_underStrip_hit[j] = stripStruct[i].hitStructCol[j].underStrip;
328 m_locx_hit[j] = stripStruct[i].hitStructCol[j].locx;
329 m_locy_hit[j] = stripStruct[i].hitStructCol[j].locy;
330 m_locz_hit[j] = stripStruct[i].hitStructCol[j].locz;
331 m_x_hit[j] = stripStruct[i].hitStructCol[j].x;
332 m_y_hit[j] = stripStruct[i].hitStructCol[j].y;
333 m_z_hit[j] = stripStruct[i].hitStructCol[j].z;
334 m_px_hit[j] = stripStruct[i].hitStructCol[j].px;
335 m_py_hit[j] = stripStruct[i].hitStructCol[j].py;
336 m_pz_hit[j] = stripStruct[i].hitStructCol[j].pz;
337 }
338 m_tree->Fill();
339 }
340 }
341}
342
343int BesTofDigitizerEcV4::calStrip( double locZ ) {
344 int strip = -1;
345 double stripWidth = m_param.strip_z + m_param.strip_gap; // Strip spread: (24+3)mm
346 int nstrip = m_param.nstrip;
347 // the offset between strip coordinate and gas SD: 0.5
348 double length = locZ + stripWidth * nstrip / 2 - 0.5;
349 if ( length <= 0 ) { strip = 0; }
350 else if ( length < stripWidth * nstrip )
351 {
352 for ( int i = 0; i < nstrip; i++ )
353 {
354 if ( length > i * stripWidth && length < ( i + 1 ) * stripWidth )
355 {
356 strip = i;
357 break;
358 }
359 }
360 }
361 else { strip = nstrip - 1; }
362 if ( strip < 0 ) strip = 0;
363 if ( strip > nstrip - 1 ) strip = nstrip - 1;
364
365 return strip;
366}
367
368bool BesTofDigitizerEcV4::underStrip( double locX, double locZ ) {
369 bool flag = 0;
370 int strip = -1;
371 double stripWidth = m_param.strip_z + m_param.strip_gap; // Strip spread: (24+3)mm
372 int nstrip = m_param.nstrip;
373 // the offset between strip coordinate and gas SD: 0.5
374 double length = locZ + stripWidth * nstrip / 2 - 0.5;
375 if ( length < stripWidth * nstrip )
376 {
377 for ( int i = 0; i < nstrip; i++ )
378 {
379 if ( length > i * stripWidth && length < ( i + 1 ) * stripWidth )
380 {
381 strip = i;
382 if ( length > i * stripWidth + m_param.strip_gap / 2 &&
383 length < ( i + 1 ) * stripWidth - m_param.strip_gap / 2 &&
384 locX > -m_param.strip_x[strip] / 2 && locX < m_param.strip_x[strip] / 2 )
385 flag = 1;
386 break;
387 }
388 }
389 }
390
391 return flag;
392}
393
394void BesTofDigitizerEcV4::StripStruct::calFirstHit() {
395 for ( unsigned int i = 0; i < hitStructCol.size(); i++ )
396 {
397 if ( hitStructCol[i].glbTime < tStart )
398 {
399 tStart = hitStructCol[i].glbTime;
400 trkIndex = hitStructCol[i].trkIndex;
401 hitStructCol[i].calTPropagate();
402 tPropagate_sphi = hitStructCol[i].tPropagate_sphi;
403 tPropagate_xphi = hitStructCol[i].tPropagate_xphi;
404 }
405 }
406}
407
408void BesTofDigitizerEcV4::HitStruct::calTPropagate() {
409 if ( strip < 0 || strip > m_param.nstrip - 1 )
410 {
411 cout << "!! BesTofDigitizerEcV4::HitStruct::calTPropagate Wrong Strip !!!" << endl;
412 return;
413 }
414
415 // It can be minus, consistent with calibration
416 double length_sphi = m_param.strip_x[strip] / 2 - locx; // mm
417 tPropagate_sphi = abs( length_sphi ) / v_propagate;
418
419 double length_xphi = m_param.strip_x[strip] / 2 + locx; // mm
420 tPropagate_xphi = abs( length_xphi ) / v_propagate;
421}
422
423double BesTofDigitizerEcV4::HitStruct::calAvaLength() {
424 // This calculation depends on the arangements of the gasLayer order and the turnover of
425 // gasContainer. all modules have the same local y trends: y larger, 11->0 In units of mm
426 double length = 0;
427 if ( gap >= 0 && gap < m_param.ngap / 2 ) length = m_param.gapWidth / 2 + locy;
428 else if ( gap < m_param.ngap ) length = m_param.gapWidth / 2 - locy;
429 else
430 {
431 cout << "BesTofDigitizerEcV4::StripStruct::calAvaLength Wrong gap calculation !!!"
432 << endl;
433 return -999.0;
434 }
435
436 return length;
437}
438
439void BesTofDigitizerEcV4::StripStruct::avalanche() {
440 // process each hit
441 for ( unsigned int i = 0; i < hitStructCol.size(); i++ )
442 {
443 hitStructCol[i].ava_pos.clear();
444 hitStructCol[i].ava_num.clear();
445
446 hitStructCol[i].ava_pos[0] = hitStructCol[i].calAvaLength();
447 hitStructCol[i].ava_num[0] = hitStructCol[i].ions;
448 // cout<<"i= "<<i<<" gap= "<<hitStructCol[i].gap<<" initial pos=
449 // "<<hitStructCol[i].ava_pos[0]<<endl;
450 for ( int j = 1; j < m_param.nstep; j++ )
451 {
452 hitStructCol[i].ava_pos[j] = hitStructCol[i].ava_pos[j - 1] + m_param.stepWidth;
453 if ( saturationFlag == true && hitStructCol[i].ava_num[j - 1] >
454 1.5e+07 ) // saturation e+07~e+08, ~2pC, Reather limit
455 { hitStructCol[i].ava_num[j] = hitStructCol[i].ava_num[j - 1]; }
456 else { hitStructCol[i].ava_num[j] = calNextN( hitStructCol[i].ava_num[j - 1] ); }
457 if ( hitStructCol[i].ava_pos[j] > m_param.gapWidth ) break;
458 }
459 }
460
461 // decide threshold and charge
462 bool over_threshold = false;
463 long int sum = 0;
464 for ( int i = 0; i < m_param.nstep; i++ )
465 {
466 for ( unsigned int j = 0; j < hitStructCol.size(); j++ )
467 {
468 if ( i < hitStructCol[j].ava_pos.size() &&
469 hitStructCol[j].ava_pos[i] < m_param.gapWidth )
470 { sum += hitStructCol[j].ava_num[i]; }
471 }
472 // cout<<"sum= "<<sum<<" avaSize= "<<hitStructCol.size()<<endl;
473
474 if ( over_threshold == false )
475 {
476 if ( sum > threshold )
477 {
478 over_threshold = true;
479 tThreshold = ( m_param.gapWidth ) / ( m_param.nstep ) / v_drift * ( i + 1 );
480 }
481 }
482 }
483
484 charge = sum * ( m_param.E_weight ) * ( v_drift ) * ( m_param.eCharge ) *
485 ( m_param.gapWidth ) / ( m_param.nstep ) / v_drift;
486}
487
488long int BesTofDigitizerEcV4::StripStruct::calNextN( int num ) {
489 if ( num < 150 )
490 {
491 long int nextN = 0;
492 double rdm;
493 for ( int i = 0; i < num; i++ )
494 {
495 rdm = G4UniformRand();
496 nextN += multiply( rdm );
497 }
498 return nextN;
499 }
500 else
501 {
502 double nbar = exp( ( alpha - eta ) * m_param.stepWidth );
503 double sigma = calSigma();
504 double mean = num * nbar;
505 double resolution = G4RandGauss::shoot( 0, ( sqrt( num ) * sigma ) );
506 long int nextN = mean + resolution;
507 return nextN;
508 }
509}
510
511long int BesTofDigitizerEcV4::StripStruct::multiply( double rdm ) {
512 double nbar = exp( ( alpha - eta ) * m_param.stepWidth );
513 double k = eta / alpha;
514 double rdm_border = k * ( nbar - 1 ) / ( nbar - k );
515 if ( rdm < rdm_border ) { return 0; }
516 else
517 {
518 long int number = 1. + 1. / log( ( nbar - 1. ) / ( nbar - k ) ) *
519 log( ( nbar - k ) * ( rdm - 1 ) / ( k - 1 ) / nbar );
520 return number;
521 }
522}
523
524double BesTofDigitizerEcV4::StripStruct::calSigma() {
525 double nbar = exp( ( alpha - eta ) * m_param.stepWidth );
526 double k = eta / alpha;
527 double sigma = sqrt( ( 1 + k ) / ( 1 - k ) * nbar * ( nbar - 1 ) );
528 return sigma;
529}
530
531double BesTofDigitizerEcV4::calTdcRes_charge( double charge_fC ) {
532 double time = 0;
533 if ( charge_fC < 250 )
534 { time = 100.764 * exp( -charge_fC * 0.0413966 + 0.377154 ) + 13.814; }
535 else { time = 12.8562 + 0.000507299 * charge_fC; }
536 if ( time < 0 ) time = 0;
537 return time; // ps
538}
539
540double BesTofDigitizerEcV4::charge2Time( double charge_fC ) {
541 double time = 0;
542 time = -120.808 / log( charge_fC * 30.1306 ) + 26.6024; // ns
543 if ( time < 0 ) time = 0;
544 return time;
545}
546
547double BesTofDigitizerEcV4::calAdcRes_charge( double charge_fC ) {
548 double time = 0;
549 if ( charge_fC < 250 )
550 { time = 72.6005 * exp( -charge_fC * 0.0302626 + 1.49059 ) + 40.8627; }
551 else { time = 32.6233 + 0.00404149 * charge_fC; }
552 if ( time < 0 ) time = 0;
553 return time; // ps
554}
555
556double BesTofDigitizerEcV4::calTdcRes_charge1( double charge_fC ) {
557 double result = 0;
558 if ( charge_fC > 50. )
559 { result = 67.6737 * exp( -charge_fC / 50.9995 - 0.27755 ) + 9.06223; }
560 else { result = 176.322 - 2.98345 * charge_fC; }
561 if ( result < 0 ) result = 0;
562 return result; // ps
563}
564
565double BesTofDigitizerEcV4::charge2Time1( double charge_fC ) {
566 double time = 0;
567 time = -4.50565 / log( charge_fC * 0.0812208 ) + 16.6653; // ns
568 if ( time < 0 ) time = 0;
569 return time;
570}
571
572double BesTofDigitizerEcV4::calAdcRes_charge1( double charge_fC ) {
573 double time = 64.3326 * exp( -charge_fC / 25.4638 + 0.944184 ) + 19.4846;
574 if ( time < 0 ) time = 0;
575 return time; // ps
576}
577
578BesTofDigitizerEcV4::HitStruct::HitStruct() { initial(); }
579
580// in units of mm or ns
581void BesTofDigitizerEcV4::HitStruct::initial() {
582 trkIndex = -999.0;
583 pdgCode = -999.0;
584 ions = -999.0;
585 strip = -999.0;
586 gap = -999.0;
587 glbTime = -999.0;
588 locx = -999.0;
589 locy = -999.0;
590 locz = -999.0;
591 x = -999.0;
592 y = -999.0;
593 z = -999.0;
594 px = -999.0;
595 py = -999.0;
596 pz = -999.0;
597 v_propagate = 0.5 * 0.299792458e+3; // mm/ns
598 tPropagate_sphi = -999.0;
599 tPropagate_xphi = -999.0;
600}
601
602BesTofDigitizerEcV4::StripStruct::StripStruct() { initial(); }
603
604void BesTofDigitizerEcV4::StripStruct::initial() {
605 // properties to get
606 strip = -999.0;
607 trkIndex = -999.0;
608 tStart = 99999.0;
609 tPropagate_sphi = -999.0;
610 tPropagate_xphi = -999.0;
611 tThreshold = -999.0;
612 charge = -999.0;
613
614 // parameters to tune
615 E = 106;
616 alpha = 144800. / 1000; //-999.0; /mm^-1
617 eta = 5013. / 1000; //-999.0; /mm^-1
618 threshold = 1.5e+08; // Correspond to induced charge of 15 fC
619 v_drift = 210.9e-3; // mm/ns
620
621 hitStructCol.clear();
622}
623
624void BesTofDigitizerEcV4::StripStruct::setPar( double alpha_n, double eta_n, double drift_n,
625 double threshold_n, bool saturationFlag_n ) {
626 alpha = alpha_n;
627 eta = eta_n;
628 v_drift = drift_n;
629 threshold = threshold_n;
630
631 saturationFlag = saturationFlag_n;
632}
633
634void BesTofDigitizerEcV4::Param::setPar( int nstep_n, double E_weight_n, double E_n ) {
635 if ( nstep_n > 0 ) nstep = nstep_n;
636 if ( E_weight_n > 0 ) E_weight = E_weight_n;
637
638 E = E_n;
639 double E_eff = E / 1000 * 2 / 6 / ( gapWidth / 10 ); // kV/cm
640 alpha = getAlpha( E_eff ); // mm^-1
641 eta = getEta( E_eff ); // mm^-1
642 v_drift = getV( E_eff ); // mm/ns
643}
644
645BesTofDigitizerEcV4::Param::Param() {
646 // parameters fixed
647 tofPara = BesTofGeoParameter::GetInstance();
648 nstrip = 12;
649 nmodule = 72;
650 std::stringstream ss;
651 for ( int i = 0; i < nstrip; i++ )
652 {
653 ss.str( "" );
654 ss << "strip_x[" << i << "]";
655 strip_x[i] = tofPara->Get( ss.str().c_str() ); // mm
656 }
657 strip_z = tofPara->Get( "strip_z" );
658 strip_gap = tofPara->Get( "strip_gap" );
659
660 ngap = 12;
661 gapWidth = 0.22; // mm
662 nstep = 200;
663 stepWidth = gapWidth / nstep;
664 E_weight = 1. / ( 6. * 0.22 + ( 5. * 0.4 + 2 * 0.55 ) / 3.7 ); // V/mm
665 eCharge = 1.60217733e-7; // pC
666 tofPara->Get_deadChannel( deadChannel );
667
668 // print();
669}
670
671double BesTofDigitizerEcV4::Param::getAlpha( double E ) {
672 // electric field: kV/cm; alpha: mm-1
673 // kV/cm
674 double e[22] = { 65, 70, 75, 80, 85, 90, 95, 100, 102, 104, 106,
675 108, 110, 112, 114, 116, 118, 120, 125, 130, 135, 140 };
676
677 // mm-1
678 double alpha[22] = {
679 383.5 / 10, 471 / 10, 564.5 / 10, 663.6 / 10, 777.1 / 10, 877 / 10,
680 990.8 / 10, 1106 / 10, 1154 / 10, 1199 / 10, 1253 / 10, 1296 / 10,
681 1344 / 10, 1396 / 10, 1448 / 10, 1502 / 10, 1545 / 10, 1597 / 10,
682 1726 / 10, 1858 / 10, 1992 / 10, 2124 / 10,
683 };
684
685 TSpline3* sp_alpha = new TSpline3( "sp_alpha", e, alpha, 22 );
686 double alphaVal = sp_alpha->Eval( E );
687 delete sp_alpha;
688 return alphaVal;
689}
690
691double BesTofDigitizerEcV4::Param::getEta( double E ) {
692 // electric field: kV/cm; eta: mm-1
693 // kV/cm
694 double e[22] = { 65, 70, 75, 80, 85, 90, 95, 100, 102, 104, 106,
695 108, 110, 112, 114, 116, 118, 120, 125, 130, 135, 140 };
696
697 // mm-1
698 double eta[22] = { 132.6 / 10, 117.2 / 10, 102.6 / 10, 88.26 / 10, 79.81 / 10, 74.0 / 10,
699 66.7 / 10, 62.7 / 10, 61.4 / 10, 57.4 / 10, 55.45 / 10, 54.35 / 10,
700 52.48 / 10, 51.3 / 10, 50.1 / 10, 48.3 / 10, 48.28 / 10, 46.00 / 10,
701 44.08 / 10, 41.67 / 10, 39.97 / 10, 38.04 / 10 };
702
703 TSpline3* sp_eta = new TSpline3( "sp_eta", e, eta, 22 );
704 double etaVal = sp_eta->Eval( E );
705 delete sp_eta;
706 return etaVal;
707}
708
709double BesTofDigitizerEcV4::Param::getV( double E ) {
710 // electric field: kV/cm; velocity: mm/ns
711 // kV/cm
712 double e[22] = { 65, 70, 75, 80, 85, 90, 95, 100, 102, 104, 106,
713 108, 110, 112, 114, 116, 118, 120, 125, 130, 135, 140 };
714
715 // mm/ns
716 double v[22] = { 130.2 / 1000, 138.5 / 1000, 146.7 / 1000, 155.0 / 1000, 163.3 / 1000,
717 171.4 / 1000, 179.7 / 1000, 187.7 / 1000, 191.2 / 1000, 194.5 / 1000,
718 197.9 / 1000, 201.2 / 1000, 204.5 / 1000, 207.6 / 1000, 210.9 / 1000,
719 214.4 / 1000, 217.5 / 1000, 220.9 / 1000, 228.8 / 1000, 237.0 / 1000,
720 244.7 / 1000, 252.9 / 1000 };
721
722 TSpline3* sp_v = new TSpline3( "sp_v", e, v, 22 );
723 double vVal = sp_v->Eval( E );
724 delete sp_v;
725 return vVal;
726}
727
728void BesTofDigitizerEcV4::Param::print() {
729 cout << "Fixed parameters: " << endl;
730 for ( int i = 0; i < nstrip; i++ ) { cout << " strip_x[" << i << "]= " << strip_x[i]; }
731 for ( int i = 0; i < nmodule; i++ )
732 {
733 for ( int j = 0; j < nstrip; j++ )
734 {
735 if ( deadChannel[i][j] != -999 )
736 { cout << " deadChannel[" << i << "][" << j << "]= " << deadChannel[i][j]; }
737 }
738 }
739
740 cout << " strip_z= " << strip_z << " strip_gap= " << strip_gap << " ngap= " << ngap
741 << " gapWidth= " << gapWidth << " nstep= " << nstep << " stepWidth= " << stepWidth
742 << " E_weight= " << E_weight << " eCharge= " << eCharge << " E= " << E
743 << " alpha= " << alpha << " eta= " << eta << " v_drift= " << v_drift << endl;
744}
745
746void BesTofDigitizerEcV4::HitStruct::print() {
747 cout << "Hit information: " << endl;
748 cout << " trkIndex= " << trkIndex << " pdgCode= " << pdgCode << " ions= " << pdgCode
749 << " strip= " << strip << " gap= " << gap << " glbTime= " << glbTime
750 << " locx= " << locx << " locy= " << locy << " locz= " << locz << " x= " << x
751 << " y= " << y << " z= " << z << " px= " << px << " py= " << py << " pz= " << pz
752 << " v_propagate= " << v_propagate << " tPropagate_sphi= " << tPropagate_sphi
753 << " tPropagate_xphi= " << tPropagate_xphi << endl;
754}
755
756void BesTofDigitizerEcV4::StripStruct::print() {
757 cout << "Strip information: " << endl;
758 cout << " strip= " << strip << " trkIndex= " << trkIndex << " tStart= " << tStart
759 << " tPropagate_sphi= " << tPropagate_sphi << " tPropagate_xphi= " << tPropagate_xphi
760 << " tThreshold " << tThreshold << " charge= " << charge << " alpha= " << alpha
761 << " eta= " << eta << " threshold= " << threshold << " v_drift= " << v_drift << endl;
762}
time
Double_t time
Definition DataBase/tau_mode.c:7
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
s
XmlRpcServer s
Definition HelloServer.cpp:11
BesTofDigitsCollection
G4TDigiCollection< BesTofDigi > BesTofDigitsCollection
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TofSim/BesTofDigi.hh:78
BesTofHitsCollection
G4THitsCollection< BesTofHit > BesTofHitsCollection
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TofSim/BesTofHit.hh:113
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition KarLud.h:35
BesTofDigitizerEcV4::calTdcRes_charge1
double calTdcRes_charge1(double charge_fC)
Definition BesTofDigitizerEcV4.cc:556
BesTofDigitizerEcV4::calAdcRes_charge
double calAdcRes_charge(double charge_fC)
Definition BesTofDigitizerEcV4.cc:547
BesTofDigitizerEcV4::calAdcRes_charge1
double calAdcRes_charge1(double charge_fC)
Definition BesTofDigitizerEcV4.cc:572
BesTofDigitizerEcV4::underStrip
bool underStrip(double locX, double locZ)
Definition BesTofDigitizerEcV4.cc:368
BesTofDigitizerEcV4::Digitize
virtual void Digitize(ScintSingle *, BesTofDigitsCollection *)
Definition BesTofDigitizerEcV4.cc:149
BesTofDigitizerEcV4::charge2Time1
double charge2Time1(double charge_fC)
Definition BesTofDigitizerEcV4.cc:565
BesTofDigitizerEcV4::charge2Time
double charge2Time(double charge_fC)
Definition BesTofDigitizerEcV4.cc:540
BesTofDigitizerEcV4::calTdcRes_charge
double calTdcRes_charge(double charge_fC)
Definition BesTofDigitizerEcV4.cc:531
BesTofGeoParameter::GetInstance
static BesTofGeoParameter * GetInstance()
Definition BesTofGeoParameter.cc:22
BesTofDigitizerEcV4::Param::setPar
void setPar(int nstep, double E_weight, double E)
Definition BesTofDigitizerEcV4.cc:634
BesTofDigitizerEcV4::StripStruct::setPar
void setPar(double alpha_n, double eta_n, double drift_n, double threshold, bool saturationFlag=true)
Definition BesTofDigitizerEcV4.cc:624
ns
#define ns(x)
Definition xmltok.c:1355