BOSS 6.6.4.p03
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:
5//Author: Matthias Ullrich
6//Created: May, 2012
7//Modified:
8//Comment:
9//---------------------------------------------------------------------------//
10//$Id: BesTofDigitizerEcV4.cc
11
12#include "BesTofDigitizerEcV4.hh"
13#include "BesTofDigi.hh"
14#include "BesTofHit.hh"
15#include "G4DigiManager.hh"
16#include "BesTofGeoParameter.hh"
17#include "Randomize.hh"
18#include "TMath.h"
19#include <math.h>
20#include "GaudiKernel/ISvcLocator.h"
21#include "GaudiKernel/Bootstrap.h"
22#include "GaudiKernel/IDataProviderSvc.h"
23#include "G4Svc/IG4Svc.h"
24#include "G4Svc/G4Svc.h"
25#include "TH1F.h"
26#include "TNtuple.h"
27#include "TFile.h"
28
29
30BesTofDigitizerEcV4::BesTofDigitizerEcV4()
31{
32
33
34 m_storedata=false;
35
36 if(m_storedata)
37 {
38 FileName = "output_mrpc_digitizerecv4.root";
39 m_TreeFile=new TFile(TString(FileName),"RECREATE","TreeFile");
40 m_ecv4tree = new TTree("ecv4tree", ""); //Create Tree
41 m_ecv4tree->Branch("m_partId",&m_partID,"m_partId/D");
42 m_ecv4tree->Branch("m_stripidentifier",&m_stripidentifier,"m_stripidentifier/D");
43 m_ecv4tree->Branch("m_trackIndex",&m_trackIndex,"m_trackIndex/D");
44 m_ecv4tree->Branch("m_signal_pc",&m_signal_pc,"m_signal_pc/D");
45 m_ecv4tree->Branch("m_time_threshold",&m_time_threshold,"m_time_threshold/D");
46 m_ecv4tree->Branch("m_time_1sthit",&m_time_1sthit,"m_time_1sthit/D");
47 m_ecv4tree->Branch("m_time_1",&m_time_1,"m_time_1/D");
48 m_ecv4tree->Branch("m_time_2",&m_time_2,"m_time_2/D");
49 m_ecv4tree->Branch("m_trans_time",&m_trans_time,"m_trans_time/D");
50 m_ecv4tree->Branch("m_firedstrip",&m_firedstrip,"m_firedstrip/D");
51 m_ecv4tree->Branch("m_numberions",&m_numberions,"m_numberions/D");
52 m_ecv4tree->Branch("m_numberofgaps_with_avalanches",&m_numberofgaps_with_avalanches,"m_numberofgaps_with_avalanches/D");
53 m_ecv4tree->Branch("m_numberofstripsfired",&m_numberofstripsfired,"m_numberofstripsfired/D");
54 m_ecv4tree->Branch("m_multihit",&m_multihit,"m_multihit/D");
55 m_ecv4tree->Branch("m_inefficiency",&m_inefficiency,"m_inefficiency/D");
56 m_ecv4tree->Branch("m_particle_true",&m_particle_true,"m_particle_true/D");
57 m_ecv4tree->Branch("m_particle_has_signal",&m_particle_has_signal,"m_particle_has_signal/D");
58
59 }
60
61}
62
63
64BesTofDigitizerEcV4::~BesTofDigitizerEcV4()
65{
66
67 if(m_storedata)
68 {
69 m_TreeFile->Write();
70 m_TreeFile->Close();
71 }
72
73}
74
75void BesTofDigitizerEcV4::Digitize(ScintSingle* single_module, BesTofDigitsCollection* DC)
76{
77
78 trackIndex=0;
79 partId=-999;
80 module_mrpc=-999;
81 time=-999.;
82
83
84 m_besTofDigitsCollection = DC;
85 G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
86 G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
87 m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
88 if (m_THC)
89 {
90
91
92
93 G4int n_hit=single_module->GetHitIndexes_mrpc()->size();
94 //G4cout << "Vor der for Schleife im TofDigitizerMRPC. NoHits " << n_hit << G4endl;
95 //G4cout << "time_1 before for loop: " << time << G4endl;
96
97
98 vector<adc_info> avalanche_info;//Will store the basic information for signal creation etc.
99 avalanche_info.clear();
100
101 partId = single_module->GetPartId();
102 module_mrpc = single_module->GetModule_mrpc();
103
104
105 //G4double edep_module = single_module->GetEdep()/keV;
106
107
108
109
110
111 //Process the hits
112 int deadarea=0;//DeadArea Counter: How often can't we associate an readoutstrip?
113
114
115 for(G4int i=0;i<n_hit;i++)
116 {
117
118
119 BesTofHit* hit = (*m_THC)[( *(single_module->GetHitIndexes_mrpc()) )[i]];
120
121
122
123
124
125 G4ThreeVector help = hit->GetPos();
126 G4int readoutstripnumber = Calculate_Readoutstrip_number(help.x()/mm,help.y()/mm,partId,module_mrpc);
127
128 if( readoutstripnumber==0 ) {deadarea++; continue;}
129
130 //G4cout << "BesTofDigitizerEcV4 Readoutstripnumber: " << partId<< "_" << module_mrpc<< " " << readoutstripnumber << G4endl;
131
132 adc_info struct_adc_info;
133 struct_adc_info.gap =0; struct_adc_info.readoutstrip =0; struct_adc_info.ions =0; struct_adc_info.zdistance =0.; struct_adc_info.trackindex=0;struct_adc_info.time=0.;
134 struct_adc_info.pdg_code=0;//Intialize all values
135
136 struct_adc_info.gap=Get_gapnumber(partId, help.z());
137 struct_adc_info.readoutstrip=readoutstripnumber;
138 struct_adc_info.ions=hit->GetIons();
139 struct_adc_info.zdistance=Distance_for_avalanche(struct_adc_info.gap,help.z(), partId);
140 struct_adc_info.time = hit->GetTime();
141 struct_adc_info.trackindex=hit->GetG4Index();
142 struct_adc_info.pdg_code=hit->GetPDGcode();
143 struct_adc_info.transtime=Calculate_strip_transition_time(help.x()/mm,help.y()/mm,partId,module_mrpc);
144 struct_adc_info.x=help.x()/mm;
145 struct_adc_info.y=help.y()/mm;
146 if(struct_adc_info.ions>0 && struct_adc_info.readoutstrip!=0 && struct_adc_info.time!=0.) avalanche_info.push_back(struct_adc_info);
147
148 /*
149 if(struct_adc_info.pdg_code==2212)
150 {
151 G4ThreeVector help2 = hit->GetMomentum();
152 cout << "#Matthias: Proton Impuls = "<<struct_adc_info.gap
153 << " " << sqrt(help2.x()*help2.x()+help2.y()*help2.y()+help2.z()*help2.z())/GeV << endl;
154
155 //cout << "#Matthias: Proton Impuls x= " << help2.x() << endl;
156 //cout << "#Matthias: Proton Impuls y= " << help2.y() << endl;
157 //cout << "#Matthias: Proton Impuls z= " << help2.z() << endl;
158
159 }
160 */
161 // G4cout <<"Gap " << struct_adc_info.gap << " Readoutstrup "<< struct_adc_info.readoutstrip << " Ions "<<struct_adc_info.ions << " "<<struct_adc_info.zdistance << G4endl;
162 //G4cout << "Particle PDG code "<< hit->GetPDGcode() << G4endl;
163
164 }//Close for-loop over hits
165
166
167
168
169 G4int avalanche_info_size = avalanche_info.size();
170 //G4cout << "avalanche_info.size() = " << avalanche_info_size << G4endl;
171
172
173
174 //Determine which readoutstrips has been fired
175 vector <int> fired_readoutstrips;
176
177 fired_readoutstrips.clear();
178 for(G4int i=0; i<avalanche_info_size; i++)
179 {
180 //if(i==0) fired_readoutstrips.push_back(avalanche_info[i].readoutstrip);
181
182 bool dlh=true;
183 for(G4int j=0;j<fired_readoutstrips.size();j++)
184 {
185 if(fired_readoutstrips[j]==avalanche_info[i].readoutstrip){dlh=false; break;}
186 //if(dlh==true) G4cout << "fired_readoutstrips[ " << j<<"] "<< fired_readoutstrips[j] << G4endl;
187 }
188 if(dlh==true)
189 {
190 fired_readoutstrips.push_back(avalanche_info[i].readoutstrip);
191 }
192 }
193
194
195
196 //Determine the adcinfo for each fired strips with the smallest time -> These information will be stored within the digis
197 vector <int> mybest_adcinfo_inreadoutstrips(fired_readoutstrips.size());
198 for(G4int i=0;i< fired_readoutstrips.size();i++)
199 {
200 double time_help=10000.;
201 for(G4int j=0; j<avalanche_info_size; j++)
202 {
203 if((time_help>avalanche_info[j].time) && avalanche_info[j].readoutstrip==fired_readoutstrips[i])
204 {
205 time_help = avalanche_info[j].time;
206 mybest_adcinfo_inreadoutstrips[i]=j;
207 }//close if
208 }//close for j
209 }//close for i
210
211
212
213
214 //Calculate the signal strength and determine the time information and store these information into the digi collection
215 double stripidentifier=0;
216 for(G4int i=0; i<fired_readoutstrips.size();i++)
217 {
218
219 //m_particle_true=0;
220 //m_particle_has_signal=0;
221
222 //Determine the total number of ions!
223 G4int sumions=0;
224 for(G4int j=0; j<avalanche_info_size; j++)
225 {
226 if(avalanche_info[j].readoutstrip==fired_readoutstrips[i])
227 {
228 sumions = sumions+ avalanche_info[j].ions;
229 //if(avalanche_info[j].pdg_code==2212) m_particle_true=1;
230 }
231
232 }
233 //if(sumions <= 10) continue;
234
235 //Determine the total number of gaps with avalanches
236 G4int numberofgaps_with_avalanches=0;
237 for(G4int j=1;j<=12;j++)
238 {
239 for(G4int k=0; k<avalanche_info_size; k++) { if(avalanche_info[k].gap ==j && avalanche_info[k].ions >0){numberofgaps_with_avalanches++; break;} }
240 }
241 //if(numberofgaps_with_avalanches < 4) continue;
242
243
244
245 G4double signal_pc=0; G4double time_of_threshold_exceedance=0;
246 simulate_avalanche(fired_readoutstrips[i], &avalanche_info,&time_of_threshold_exceedance,&signal_pc);
247 signal_pc=signal_pc;//1e-12*A*s;
248 time_of_threshold_exceedance=time_of_threshold_exceedance*s;//Stores the time when the threshold is reached.
249
250
251 if(signal_pc<0.015) continue; //We only store signals wich are over the threshold
252
253 //if(m_particle_true==1) m_particle_has_signal=1;
254
255 //G4cout << "singal " << setprecision(10) << signal_pc << " pC" << endl;
256 //G4cout << "test_time " << time_of_threshold_exceedance / ns << endl;
257
258
259 trackIndex = avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].trackindex;
260
261
262
263 // We add the threshold crossing time for the avalanche. This will be corrected later within the reconstruction.
264 time = avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time + time_of_threshold_exceedance + avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime;
265
266 //std::cout << "Digitizer time " << avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time << std::endl;
267 //std::cout<< "Digitizer threshold time " << time_of_threshold_exceedance<<std::endl;
268 //std::cout << "Digitizer time+ threshold " << time << std::endl;
269
270
271
272 //time smear #Correct
273 //Uncertainities which are not considered in the digitization process. Following values are assumed:
274 //Avalanche Process : Already considered with time of threshold!
275 //FEE (NINO) : Depending on the charge, will be calculated in the next step! = NINO_RES
276 //Hit position : 10 ps
277 //Readout : 27 ps
278 //Electronics/cable : 20 ps
279 //-----------------------------------------
280 // Sum =sqrt(10**2 + 27**2 +20**2 NINO_RES**2 ) time smear
281
282 double femtocolomb= 1000.;
283
284 double NINO_RES=0;
285 NINO_RES = Get_NINO_leadingedge_resolution(signal_pc*femtocolomb);
286 double smear_time = sqrt(NINO_RES*NINO_RES + 27.0*27.0 + 10.0*10.0 + 20.0*20.0); //This is in ps
287 time = Smear_gaussian(time,0,smear_time*0.001*ns);
288
289
290
291 stripidentifier= Produce_unique_identifier(module_mrpc,avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].readoutstrip);
292
293 //Resolution of Pulse length:
294 // 27 ps TDC
295 // Timeresolution of pulse length
296 double NINO_pulselength_resolution = Get_NINO_pulselength_resolution( signal_pc*femtocolomb );
297 double smear_time_2 = sqrt(NINO_pulselength_resolution*NINO_pulselength_resolution + 27.0*27.0 );
298 double time2 =chargetotime(signal_pc*femtocolomb)+time; //The function returns the time in ns
299 time2 = Smear_gaussian(time2,0,smear_time_2*0.001*ns);
300
301
302 if(m_storedata)
303 {
304 m_partID=partId;
305 m_stripidentifier=stripidentifier;
306 m_trackIndex=trackIndex;
307 m_signal_pc=signal_pc;
308 m_time_threshold= time_of_threshold_exceedance;
309 m_time_1sthit=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].time;
310 m_time_1=time;
311 m_time_2=time2;
312 m_trans_time=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].transtime;
313 m_firedstrip=fired_readoutstrips[i];
314 m_numberions=sumions;
315 m_numberofgaps_with_avalanches=numberofgaps_with_avalanches;
316 m_numberofstripsfired = fired_readoutstrips.size();
317 m_x=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].x;
318 m_y=avalanche_info[(mybest_adcinfo_inreadoutstrips[i])].y;
319 m_ecv4tree->Fill();
320 }
321
322
323 //std::cout << "Digitizer time "<<time << std::endl;
324 //std::cout << "Digitizer time2 "<<time2 <<std::endl;
325
326
327 BesTofDigi *digi = new BesTofDigi;
328 digi->SetTrackIndex(trackIndex);
329 digi->SetPartId(partId);
330 digi->SetModule_mrpc(stripidentifier);
331 digi->SetTime1_mrpc(time);
332 digi->SetTime2_mrpc(time2);
333 m_besTofDigitsCollection->insert(digi);
334
335
336 }//close i
337
338
339
340
341 }//close if(m_THC)
342
343} //Digitize()
344
345
346
347
348
349//////////////////////////////////////////////////////////
350//////////////////////////////////////////////////////////
351//////
352////// Function used within the Digitization process ////
353//////
354/////////////////////////////////////////////////////////
355/////////////////////////////////////////////////////////
356
357
358
359G4double BesTofDigitizerEcV4::Smear_gaussian(G4double variable, G4double gauss_mean,G4double gauss_sigma)
360{
361
362 G4double gaussrand = G4RandGauss::shoot(gauss_mean,gauss_sigma);//returns x-value
363 return (variable+gaussrand);
364
365}
366
367
368
369//This function caluclates depending on the energy deposition which strips cause a signal!!
370
371G4int BesTofDigitizerEcV4::Calculate_Readoutstrip_number(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
372
373{
374
375//Within this functions all numbers are in mm or in degree!!!!
376
377 G4int readoutstripnumber=0;
378
379 G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm);
380 G4double phi= atan2(y_mm,x_mm)*180./pi;
381
382
383 if(-180<phi && phi<90){phi=phi+270;} // Gehe in ein Koordinatensystem, was oben 0 hat und dann nach rechts/links herum anwaechst! Kommt auf die Seite an!
384 else {phi=phi-90;}
385
386
387 G4double result_phi=0;
388 if(partId_f!=3 && partId_f!=6) result_phi=0;
389 else result_phi=9.97312; //This is the angle to which the first and second layer are shifted in degree!!!
390
391
392
393 if(partId_f==3 || partId_f==4) result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
394 else result_phi = (module_mrpc_f-1)*20 - phi + result_phi; // If > 0 Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
395
396 //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
397 if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
398 if(partId_f==5 && module_mrpc_f==1 & result_phi<-350) result_phi=result_phi+360;
399
400
401
402 result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5 - 3.82; //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transformation of the readout strips...
403 // 3.82 is the ks_trafo
404
405 double result_x=fabs((result_radius/cos(result_phi*pi/180.0)+634.56-159.5+3.82)*sin(result_phi*pi/180)); //Calculate the resulting position x!
406
407
408 if(result_phi>=0 && result_x>2.)
409 {
410 if( -12.5 <result_radius && result_radius < 12.5 && result_x <= 44 ) readoutstripnumber=2;
411 if( (-12.5 +(25.0+4.0)*1 ) <result_radius && result_radius < (12.5+(25.0+4.0)*1) && result_x <=46) readoutstripnumber=4;
412 if( (-12.5 +(25.0+4.0)*2 ) <result_radius && result_radius < (12.5+(25.0+4.0)*2) && result_x <=49) readoutstripnumber=6;
413 if( (-12.5 +(25.0+4.0)*3 ) <result_radius && result_radius < (12.5+(25.0+4.0)*3) && result_x <=51) readoutstripnumber=8;
414 if( (-12.5 +(25.0+4.0)*4 ) <result_radius && result_radius < (12.5+(25.0+4.0)*4) && result_x <=54) readoutstripnumber=10;
415 if( (-12.5 +(25.0+4.0)*5 ) <result_radius && result_radius < (12.5+(25.0+4.0)*5) && result_x <=56) readoutstripnumber=12;
416 if( (-12.5 +(25.0+4.0)*6 ) <result_radius && result_radius < (12.5+(25.0+4.0)*6) && result_x <=59) readoutstripnumber=14;
417 if( (-12.5 +(25.0+4.0)*7 ) <result_radius && result_radius < (12.5+(25.0+4.0)*7) && result_x <=61) readoutstripnumber=16;
418 if( (-12.5 +(25.0+4.0)*8 ) <result_radius && result_radius < (12.5+(25.0+4.0)*8) && result_x <=64) readoutstripnumber=18;
419 if( (-12.5 +(25.0+4.0)*9 ) <result_radius && result_radius < (12.5+(25.0+4.0)*9) && result_x <=66) readoutstripnumber=20;
420 if( (-12.5 +(25.0+4.0)*10 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10) && result_x <=69) readoutstripnumber=22;
421 if( (-12.5 +(25.0+4.0)*11 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11) && result_x <=71) readoutstripnumber=24;
422 }
423 else if(result_phi<0 && result_x>2.)
424 {
425 if( -12.5 <result_radius && result_radius < 12.5 && result_x <= 44 ) readoutstripnumber=1;
426 if( (-12.5 +(25.0+4.0)*1 ) <result_radius && result_radius < (12.5+(25.0+4.0)*1) && result_x<=46 ) readoutstripnumber=3;
427 if( (-12.5 +(25.0+4.0)*2 ) <result_radius && result_radius < (12.5+(25.0+4.0)*2) && result_x<=49) readoutstripnumber=5;
428 if( (-12.5 +(25.0+4.0)*3 ) <result_radius && result_radius < (12.5+(25.0+4.0)*3) && result_x<=51) readoutstripnumber=7;
429 if( (-12.5 +(25.0+4.0)*4 ) <result_radius && result_radius < (12.5+(25.0+4.0)*4) && result_x<=54) readoutstripnumber=9;
430 if( (-12.5 +(25.0+4.0)*5 ) <result_radius && result_radius < (12.5+(25.0+4.0)*5) && result_x<=56) readoutstripnumber=11;
431 if( (-12.5 +(25.0+4.0)*6 ) <result_radius && result_radius < (12.5+(25.0+4.0)*6) && result_x<=59) readoutstripnumber=13;
432 if( (-12.5 +(25.0+4.0)*7 ) <result_radius && result_radius < (12.5+(25.0+4.0)*7) && result_x<=61) readoutstripnumber=15;
433 if( (-12.5 +(25.0+4.0)*8 ) <result_radius && result_radius < (12.5+(25.0+4.0)*8) && result_x<=64) readoutstripnumber=17;
434 if( (-12.5 +(25.0+4.0)*9 ) <result_radius && result_radius < (12.5+(25.0+4.0)*9) && result_x<=66) readoutstripnumber=19;
435 if( (-12.5 +(25.0+4.0)*10 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10) && result_x<=69) readoutstripnumber=21;
436 if( (-12.5 +(25.0+4.0)*11 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11) && result_x<=71) readoutstripnumber=23;
437 }
438
439 //G4cout << " result_radius " << result_radius << G4endl;
440 //G4cout << " result_x " << result_x << G4endl;
441 // G4cout << G4endl<< "phi "<< phi << G4endl;
442 // G4cout << G4endl<< "result_phi "<< result_phi << G4endl;
443 //G4cout << G4endl<< "module_mrpc_f "<< module_mrpc_f*20<< G4endl;
444 //G4cout << G4endl<< "Number: "<< readoutstripnumber<< G4endl;
445
446 return readoutstripnumber;
447
448}
449
450
451
452
453//This function is almost the same as above, but it does not consider the "dead" readout area! This is important for the reconstruction!
454//So it is possible to determine the readoutstrip, even if no strip is hit and one can do the extrapolation for the MDC Tracks
455G4int BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
456{
457
458 //Within this functions all numbers are in mm or in degree!!!!
459 //G4cout << "x | y | part | module " << x_mm << " " << y_mm << " " << partId_f << " " << module_mrpc_f << G4endl;
460
461 G4int readoutstripnumber=0;
462
463 G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm);
464 G4double phi= atan2(y_mm,x_mm)*180./pi;
465
466
467 if(-180<phi && phi<90){phi=phi+270;} // Gehe in ein Koordinatensystem, was oben 0 hat und dann nach rechts/links herum anwächst! Kommt auf die Seite an!
468 else {phi=phi-90;}
469
470
471 G4double result_phi=0;
472 if(partId_f!=3 && partId_f!=6) result_phi=0;
473 else result_phi=9.97312; //This is the angle to which the first and second layer are shifted in degree!!!
474 if(partId_f==3 || partId_f==4) result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
475 else result_phi = (module_mrpc_f-1)*20 - phi + result_phi; // If > 0 Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
476
477 //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
478 if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
479 if(partId_f==5 && module_mrpc_f==1 & result_phi<-350) result_phi=result_phi+360;
480
481
482
483 result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5-3.82; //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transformation of the readout strips...
484
485 //Die +-12.0 sind beliebing gewaehlt. In diesesm Bereich werden Treffer noch als mrpc treffer gewertet.
486 if(result_phi>=0)
487 {
488 if( (-12.5 -12.0) <=result_radius && result_radius < 12.5+2.0 ) readoutstripnumber=2;
489 if( (-12.5 +(25.0+4.0)*1-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*1)+2.0 ) readoutstripnumber=4;
490 if( (-12.5 +(25.0+4.0)*2-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*2) +2.0) readoutstripnumber=6;
491 if( (-12.5 +(25.0+4.0)*3-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*3) +2.0) readoutstripnumber=8;
492 if( (-12.5 +(25.0+4.0)*4-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*4) +2.0) readoutstripnumber=10;
493 if( (-12.5 +(25.0+4.0)*5-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*5) +2.0) readoutstripnumber=12;
494 if( (-12.5 +(25.0+4.0)*6-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*6) +2.0) readoutstripnumber=14;
495 if( (-12.5 +(25.0+4.0)*7-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*7) +2.0) readoutstripnumber=16;
496 if( (-12.5 +(25.0+4.0)*8-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*8) +2.0) readoutstripnumber=18;
497 if( (-12.5 +(25.0+4.0)*9-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*9) +2.0) readoutstripnumber=20;
498 if( (-12.5 +(25.0+4.0)*10-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*10) +2.0) readoutstripnumber=22;
499 if( (-12.5 +(25.0+4.0)*11-2.0 ) <=result_radius && result_radius < (12.5+(25.0+4.0)*11) +12.0) readoutstripnumber=24;
500 }
501 else if(result_phi<0)
502 {
503 if( (-12.5-12.0) <=result_radius && result_radius < (12.5+2.0) ) readoutstripnumber=1;
504 if( (-12.5 +(25.0+4.0)*1 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*1) +2.0) readoutstripnumber=3;
505 if( (-12.5 +(25.0+4.0)*2 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*2) +2.0) readoutstripnumber=5;
506 if( (-12.5 +(25.0+4.0)*3 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*3) +2.0) readoutstripnumber=7;
507 if( (-12.5 +(25.0+4.0)*4 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*4) +2.0) readoutstripnumber=9;
508 if( (-12.5 +(25.0+4.0)*5 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*5) +2.0) readoutstripnumber=11;
509 if( (-12.5 +(25.0+4.0)*6 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*6) +2.0) readoutstripnumber=13;
510 if( (-12.5 +(25.0+4.0)*7 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*7) +2.0) readoutstripnumber=15;
511 if( (-12.5 +(25.0+4.0)*8 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*8) +2.0) readoutstripnumber=17;
512 if( (-12.5 +(25.0+4.0)*9 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*9) +2.0) readoutstripnumber=19;
513 if( (-12.5 +(25.0+4.0)*10 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*10) +2.0) readoutstripnumber=21;
514 if( (-12.5 +(25.0+4.0)*11 -2.0) <=result_radius && result_radius < (12.0+(25.0+4.0)*11) +12.0) readoutstripnumber=23;
515 }
516
517 // G4cout << " result_radius " << result_radius << G4endl;
518 // G4cout << G4endl<< "phi "<< phi << G4endl;
519 // G4cout << G4endl<< "result_phi "<< result_phi << G4endl;
520 // G4cout << G4endl<< "module_mrpc_f "<< module_mrpc_f*20<< G4endl;
521 // G4cout << G4endl<< "Number: "<< readoutstripnumber<< G4endl;
522
523 return readoutstripnumber;
524
525}
526
527
528
529
530
531
532
533
534G4double BesTofDigitizerEcV4::Calculate_strip_transition_time(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
535{
536 G4double signal_velocity_copper_mrpc_strip= 0.8*0.299792458; // unit: percentage * c_0 * mm/ps
537
538
539 //Within this functions all numbers are in mm or in degree!!!!
540 G4double length_to_border_of_the_readoutstrip=0;
541
542 G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm); //Radius where the particle did hit the detector
543 G4double phi= atan2(y_mm,x_mm)*180./pi; //Phi angle of the particle hitting the detector
544
545 /*
546 cout << "module "<< module_mrpc_f<< endl;
547 cout << "partID "<< partId_f<< endl;
548 cout << "r "<< result_radius<<endl;
549 cout << "phi "<< phi/180.*pi<<endl;
550 cout << "x "<<x_mm/10. << endl;
551 cout << "y "<<y_mm/10. << endl;
552 */
553
554
555 if(-180<phi && phi<90){phi=phi+270;} // Transformation into a coordinate-system, whose coordinates do increase to the right/left side. Left or right does depend on the endcap!
556 else {phi=phi-90;}
557
558 G4double result_phi=0;
559 if(partId_f!=3 && partId_f!=6) result_phi=0;
560 else result_phi=9.97312; //This is the angle to which the first and second layer are shifted in degree!!!
561
562
563 if(partId_f==3 || partId_f==4) result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
564 else result_phi = (module_mrpc_f-1)*20 - phi + result_phi; // If > 0 Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
565 //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
566
567 if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
568 if(partId_f==5 && module_mrpc_f==1 & result_phi<-350) result_phi=result_phi+360;
569
570 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
571 //The Phi is now the is the phi angle for each module and not for the whole assembly
572 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
573
574
575
576 result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5 - 3.82; //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transformation of the readout strips...
577 // 3.82 is the ks_trafo
578
579 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
580 //The radius is now the the radius within one module. The point of origin for the
581 //radius is the middle of the lowest readout cell
582 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
583
584
585 double result_x=fabs((result_radius/cos(result_phi*pi/180.0)+634.56-159.5+3.82)*sin(result_phi*pi/180)); //Calculate the resulting position x!
586
587 //The two main if cases are not necessary, but I did only copy the fucnction from the one BesTofDigitizerEcV4::Calculate_Readoutstrip_number(....)
588 if(result_phi>=0 && result_x>2.)
589 {
590 if( -12.5 <result_radius && result_radius < 12.5 && result_x <= 44 )
591 {
592 length_to_border_of_the_readoutstrip=sqrt( (44-result_x)*(44-result_x) + fabs((result_radius-0.)*(result_radius-0.)));
593 //cout << "y_pos " << fabs((result_radius-0.)*(result_radius-0.))<<endl;
594 }
595 if( (-12.5 +(25.0+4.0)*1 ) <result_radius && result_radius < (12.5+(25.0+4.0)*1) && result_x <=46)
596 {
597 length_to_border_of_the_readoutstrip=sqrt( (46-result_x)*(46-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
598 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*1)*(result_radius-(25.0+4.0)*1))<<endl;
599 }
600 if( (-12.5 +(25.0+4.0)*2 ) <result_radius && result_radius < (12.5+(25.0+4.0)*2) && result_x <=49)
601 {
602 length_to_border_of_the_readoutstrip=sqrt( (49-result_x)*(49-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
603 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*2)*(result_radius-(25.0+4.0)*2))<<endl;
604 }
605 if( (-12.5 +(25.0+4.0)*3 ) <result_radius && result_radius < (12.5+(25.0+4.0)*3) && result_x <=51)
606 {
607 length_to_border_of_the_readoutstrip=sqrt( (51-result_x)*(51-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
608 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*3)*(result_radius-(25.0+4.0)*3))<<endl;
609 }
610 if( (-12.5 +(25.0+4.0)*4 ) <result_radius && result_radius < (12.5+(25.0+4.0)*4) && result_x <=54)
611 {
612 length_to_border_of_the_readoutstrip=sqrt( (54-result_x)*(54-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
613 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*4)*(result_radius-(25.0+4.0)*4))<<endl;
614 }
615 if( (-12.5 +(25.0+4.0)*5 ) <result_radius && result_radius < (12.5+(25.0+4.0)*5) && result_x <=56)
616 {
617 length_to_border_of_the_readoutstrip=sqrt( (56-result_x)*(56-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
618 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*5)*(result_radius-(25.0+4.0)*5))<<endl;
619 }
620 if( (-12.5 +(25.0+4.0)*6 ) <result_radius && result_radius < (12.5+(25.0+4.0)*6) && result_x <=59)
621 {
622 length_to_border_of_the_readoutstrip=sqrt( (59-result_x)*(59-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
623 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*6)*(result_radius-(25.0+4.0)*6))<<endl;
624 }
625 if( (-12.5 +(25.0+4.0)*7 ) <result_radius && result_radius < (12.5+(25.0+4.0)*7) && result_x <=61)
626 {
627 length_to_border_of_the_readoutstrip=sqrt( (61-result_x)*(61-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
628 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*7)*(result_radius-(25.0+4.0)*7))<<endl;
629 }
630 if( (-12.5 +(25.0+4.0)*8 ) <result_radius && result_radius < (12.5+(25.0+4.0)*8) && result_x <=64)
631 {
632 length_to_border_of_the_readoutstrip=sqrt( (64-result_x)*(64-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
633 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*8)*(result_radius-(25.0+4.0)*8))<<endl;
634 }
635 if( (-12.5 +(25.0+4.0)*9 ) <result_radius && result_radius < (12.5+(25.0+4.0)*9) && result_x <=66)
636 {
637 length_to_border_of_the_readoutstrip=sqrt( (66-result_x)*(66-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
638 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*9)*(result_radius-(25.0+4.0)*9))<<endl;
639 }
640 if( (-12.5 +(25.0+4.0)*10 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10) && result_x <=69)
641 {
642 length_to_border_of_the_readoutstrip=sqrt( (69-result_x)*(69-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
643 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*10)*(result_radius-(25.0+4.0)*10))<<endl;
644 }
645 if( (-12.5 +(25.0+4.0)*11 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11) && result_x <=71)
646 {
647 length_to_border_of_the_readoutstrip=sqrt( (71-result_x)*(71-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
648 //cout << "y_pos " << fabs((result_radius-(25.0+4.0)*11)*(result_radius-(25.0+4.0)*11))<<endl;
649 }
650 }
651 else if(result_phi<0 && result_x>2.)
652 {
653 if( -12.5 <result_radius && result_radius < 12.5 && result_x <= 44 )
654 length_to_border_of_the_readoutstrip=sqrt( (44-result_x)*(44-result_x) + fabs((result_radius-0.)*(result_radius-0.)));
655 if( (-12.5 +(25.0+4.0)*1 ) <result_radius && result_radius < (12.5+(25.0+4.0)*1) && result_x <=46)
656 length_to_border_of_the_readoutstrip=sqrt( (46-result_x)*(46-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
657 if( (-12.5 +(25.0+4.0)*2 ) <result_radius && result_radius < (12.5+(25.0+4.0)*2) && result_x <=49)
658 length_to_border_of_the_readoutstrip=sqrt( (49-result_x)*(49-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
659 if( (-12.5 +(25.0+4.0)*3 ) <result_radius && result_radius < (12.5+(25.0+4.0)*3) && result_x <=51)
660 length_to_border_of_the_readoutstrip=sqrt( (51-result_x)*(51-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
661 if( (-12.5 +(25.0+4.0)*4 ) <result_radius && result_radius < (12.5+(25.0+4.0)*4) && result_x <=54)
662 length_to_border_of_the_readoutstrip=sqrt( (54-result_x)*(54-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
663 if( (-12.5 +(25.0+4.0)*5 ) <result_radius && result_radius < (12.5+(25.0+4.0)*5) && result_x <=56)
664 length_to_border_of_the_readoutstrip=sqrt( (56-result_x)*(56-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
665 if( (-12.5 +(25.0+4.0)*6 ) <result_radius && result_radius < (12.5+(25.0+4.0)*6) && result_x <=59)
666 length_to_border_of_the_readoutstrip=sqrt( (59-result_x)*(59-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
667 if( (-12.5 +(25.0+4.0)*7 ) <result_radius && result_radius < (12.5+(25.0+4.0)*7) && result_x <=61)
668 length_to_border_of_the_readoutstrip=sqrt( (61-result_x)*(61-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
669 if( (-12.5 +(25.0+4.0)*8 ) <result_radius && result_radius < (12.5+(25.0+4.0)*8) && result_x <=64)
670 length_to_border_of_the_readoutstrip=sqrt( (64-result_x)*(64-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
671 if( (-12.5 +(25.0+4.0)*9 ) <result_radius && result_radius < (12.5+(25.0+4.0)*9) && result_x <=66)
672 length_to_border_of_the_readoutstrip=sqrt( (66-result_x)*(66-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
673 if( (-12.5 +(25.0+4.0)*10 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10) && result_x <=69)
674 length_to_border_of_the_readoutstrip=sqrt( (69-result_x)*(69-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
675 if( (-12.5 +(25.0+4.0)*11 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11) && result_x <=71)
676 length_to_border_of_the_readoutstrip=sqrt( (71-result_x)*(71-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
677 }
678
679
680 double transition_time = length_to_border_of_the_readoutstrip/signal_velocity_copper_mrpc_strip/1000.; //length is in mm , 3*10^8 m/s = 0.3 mm/ps --> Unit of transitionstime is ns (due to division by 1000.)
681 return transition_time;
682
683}//close calculate_strip_transition_time
684
685
686
687
688
689//This function is the same as above, but the readoutsttrips are assumed to be a little bit
690//larger (continues). The reason is that during the reconstruction the extrapolated
691//is sometimes in between some chip or at the border. At the border the pads are extended by 2mm
692G4double BesTofDigitizerEcV4::Calculate_strip_transition_time_cont(G4double x_mm,G4double y_mm, G4int partId_f, G4int module_mrpc_f)
693{
694 G4double signal_velocity_copper_mrpc_strip= 0.8*0.299792458; // unit: percentage * c_0 * mm/ps
695
696
697 //Within this functions all numbers are in mm or in degree!!!!
698 G4double length_to_border_of_the_readoutstrip=0;
699
700 G4double result_radius= sqrt(x_mm*x_mm+y_mm*y_mm); //Radius where the particle did hit the detector
701 G4double phi= atan2(y_mm,x_mm)*180./pi; //Phi angle of the particle hitting the detector
702
703 /*
704 cout << "module "<< module_mrpc_f<< endl;
705 cout << "partID "<< partId_f<< endl;
706 cout << "r "<< result_radius<<endl;
707 cout << "phi "<< phi/180.*pi<<endl;
708 cout << "x "<<x_mm/10. << endl;
709 cout << "y "<<y_mm/10. << endl;
710 */
711
712
713 if(-180<phi && phi<90){phi=phi+270;} // Transformation into a coordinate-system, whose coordinates do increase to the right/left side. Left or right does depend on the endcap!
714 else {phi=phi-90;}
715
716 G4double result_phi=0;
717 if(partId_f!=3 && partId_f!=6) result_phi=0;
718 else result_phi=9.97312; //This is the angle to which the first and second layer are shifted in degree!!!
719
720
721 if(partId_f==3 || partId_f==4) result_phi = -((19-module_mrpc_f-1)*20 - phi + result_phi);
722 else result_phi = (module_mrpc_f-1)*20 - phi + result_phi; // If > 0 Readoutstripnumber -> 2,4,6,8,10,12,14.16,18,20,22,24
723 //The 18th and the 1st module are located between 350 and 10 degree. This has to be taken into account
724
725 if(partId_f==4 && module_mrpc_f==18 && result_phi>350) result_phi=result_phi-360;
726 if(partId_f==5 && module_mrpc_f==1 & result_phi<-350) result_phi=result_phi+360;
727
728 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
729 //The Phi is now the is the phi angle for each module and not for the whole assembly
730 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
731
732
733
734 result_radius= result_radius*cos(result_phi*pi/180.0) - 634.56 + 159.5 - 3.82; //634.56 is the origin of coordinates from the module, 159.5 is the distance from the transform\
735ation of the readout strips...
736 // 3.82 is the ks_trafo
737
738 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
739 //The radius is now the the radius within one module. The point of origin for the
740 //radius is the middle of the lowest readout cell
741 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
742
743
744 double result_x=fabs((result_radius/cos(result_phi*pi/180.0)+634.56-159.5+3.82)*sin(result_phi*pi/180)); //Calculate the resulting position x!
745
746//The two main if cases are not necessary, but I did only copy the fucnction from the one BesTofDigitizerEcV4::Calculate_Readoutstrip_number(....)
747 if(result_phi>=0 && result_x>0.)
748 {
749 if( (-12.5-2.0) <result_radius && result_radius < (12.5+2.0) && result_x <= (44+2.0) )
750 length_to_border_of_the_readoutstrip=sqrt( (44-result_x)*(44-result_x) + fabs((result_radius-0.)*(result_radius-0.)));
751 if( (-12.5 +(25.0+4.0)*1-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*1+2.0) && result_x <=(46+2.0))
752 length_to_border_of_the_readoutstrip=sqrt( (46-result_x)*(46-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
753 if( (-12.5 +(25.0+4.0)*2-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*2+2.0) && result_x <=(49+2.0))
754 length_to_border_of_the_readoutstrip=sqrt( (49-result_x)*(49-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
755 if( (-12.5 +(25.0+4.0)*3-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*3+2.0) && result_x <=(51+2.0))
756 length_to_border_of_the_readoutstrip=sqrt( (51-result_x)*(51-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
757 if( (-12.5 +(25.0+4.0)*4-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*4+2.0) && result_x <=(54+2.0))
758 length_to_border_of_the_readoutstrip=sqrt( (54-result_x)*(54-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
759 if( (-12.5 +(25.0+4.0)*5-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*5+2.0) && result_x <=(56+2.0))
760 length_to_border_of_the_readoutstrip=sqrt( (56-result_x)*(56-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
761 if( (-12.5 +(25.0+4.0)*6-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*6+2.0) && result_x <=(59+2.0))
762 length_to_border_of_the_readoutstrip=sqrt( (59-result_x)*(59-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
763 if( (-12.5 +(25.0+4.0)*7-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*7+2.0) && result_x <=(61+2.0))
764 length_to_border_of_the_readoutstrip=sqrt( (61-result_x)*(61-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
765 if( (-12.5 +(25.0+4.0)*8-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*8+2.0) && result_x <=(64+2.0))
766 length_to_border_of_the_readoutstrip=sqrt( (64-result_x)*(64-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
767 if( (-12.5 +(25.0+4.0)*9-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*9+2.0) && result_x <=(66+2.0))
768 length_to_border_of_the_readoutstrip=sqrt( (66-result_x)*(66-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
769 if( (-12.5 +(25.0+4.0)*10-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10+2.0) && result_x <=(69+2.0))
770 length_to_border_of_the_readoutstrip=sqrt( (69-result_x)*(69-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
771 if( (-12.5 +(25.0+4.0)*11-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11+2.0) && result_x <=(71+2.0))
772 length_to_border_of_the_readoutstrip=sqrt( (71-result_x)*(71-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
773 }
774 else if(result_phi<0 && result_x>0.)
775 {
776 if( (-12.5-2.0) <result_radius && result_radius < (12.5+2.0) && result_x <= (44+2.0) )
777 length_to_border_of_the_readoutstrip=sqrt( (44-result_x)*(44-result_x) + fabs((result_radius-0.)*(result_radius-0.)));
778 if( (-12.5 +(25.0+4.0)*1-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*1+2.0) && result_x <=(46+2.0))
779 length_to_border_of_the_readoutstrip=sqrt( (46-result_x)*(46-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*1)));
780 if( (-12.5 +(25.0+4.0)*2-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*2+2.0) && result_x <=(49+2.0))
781 length_to_border_of_the_readoutstrip=sqrt( (49-result_x)*(49-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*2)));
782 if( (-12.5 +(25.0+4.0)*3-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*3+2.0) && result_x <=(51+2.0))
783 length_to_border_of_the_readoutstrip=sqrt( (51-result_x)*(51-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*3)));
784 if( (-12.5 +(25.0+4.0)*4-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*4+2.0) && result_x <=(54+2.0))
785 length_to_border_of_the_readoutstrip=sqrt( (54-result_x)*(54-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*4)));
786 if( (-12.5 +(25.0+4.0)*5-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*5+2.0) && result_x <=(56+2.0))
787 length_to_border_of_the_readoutstrip=sqrt( (56-result_x)*(56-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*5)));
788 if( (-12.5 +(25.0+4.0)*6-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*6+2.0) && result_x <=(59+2.0))
789 length_to_border_of_the_readoutstrip=sqrt( (59-result_x)*(59-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*6)));
790 if( (-12.5 +(25.0+4.0)*7-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*7+2.0) && result_x <=(61+2.0))
791 length_to_border_of_the_readoutstrip=sqrt( (61-result_x)*(61-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*7)));
792 if( (-12.5 +(25.0+4.0)*8-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*8+2.0) && result_x <=(64+2.0))
793 length_to_border_of_the_readoutstrip=sqrt( (64-result_x)*(64-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*8)));
794 if( (-12.5 +(25.0+4.0)*9-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*9+2.0) && result_x <=(66+2.0))
795 length_to_border_of_the_readoutstrip=sqrt( (66-result_x)*(66-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*9)));
796 if( (-12.5 +(25.0+4.0)*10-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*10+2.0) && result_x <=(69+2.0))
797 length_to_border_of_the_readoutstrip=sqrt( (69-result_x)*(69-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*10)));
798 if( (-12.5 +(25.0+4.0)*11-2.0 ) <result_radius && result_radius < (12.5+(25.0+4.0)*11+2.0) && result_x <=(71+2.0))
799 length_to_border_of_the_readoutstrip=sqrt( (71-result_x)*(71-result_x) + fabs((result_radius-(25.0+4.0)*1.)*(result_radius-(25.0+4.0)*11)));
800 }
801
802 double transition_time = length_to_border_of_the_readoutstrip/signal_velocity_copper_mrpc_strip/1000.; //length is in mm , 3*10^8 m/s = 0.3 mm/ps --> Unit of transitionstime is ns (due to division by 1000.)
803 return transition_time;
804
805}//close calculate_strip_transition_time
806
807
808//This is the transitiontime from the middle of the module
809//This value will be used if the reconstruction did not work proberly!
810G4double BesTofDigitizerEcV4::Average_transition_time(G4int my_module)
811{
812 G4double signal_velocity_copper_mrpc_strip= 0.8*0.299792458; // unit: percentage * c_0 * mm/ps
813 G4double length_to_border_of_the_readoutstrip=0.;
814 if(my_module==1 || my_module==2) length_to_border_of_the_readoutstrip=sqrt( (44.-22.)*(44.-22.));
815 if(my_module==3 || my_module==4) length_to_border_of_the_readoutstrip=sqrt( (46.-23.)*(46.-23.));
816 if(my_module==5 || my_module==6) length_to_border_of_the_readoutstrip=sqrt( (49.-24.5)*(49.-24.5));
817 if(my_module==7 || my_module==8) length_to_border_of_the_readoutstrip=sqrt( (51.-25.5)*(51.-25.5));
818 if(my_module==9 || my_module==10) length_to_border_of_the_readoutstrip=sqrt( (54.-27.)*(54.-27.));
819 if(my_module==11 || my_module==12) length_to_border_of_the_readoutstrip=sqrt( (56.-28.)*(56.-28.));
820 if(my_module==13 || my_module==14) length_to_border_of_the_readoutstrip=sqrt( (59.-29.5)*(59.-29.5));
821 if(my_module==15 || my_module==16) length_to_border_of_the_readoutstrip=sqrt( (61.-30.5)*(61.-30.5));
822 if(my_module==17 || my_module==18) length_to_border_of_the_readoutstrip=sqrt( (64.-32.)*(64.-32.));
823 if(my_module==19 || my_module==20) length_to_border_of_the_readoutstrip=sqrt( (66.-33.)*(66.-33.));
824 if(my_module==21 || my_module==22) length_to_border_of_the_readoutstrip=sqrt( (69.-34.5)*(69.-34.5));
825 if(my_module==23 || my_module==24) length_to_border_of_the_readoutstrip=sqrt( (71.-35.5)*(71.-35.5));
826
827 double transition_time = length_to_border_of_the_readoutstrip/signal_velocity_copper_mrpc_strip/1000.;
828 return transition_time;
829}
830
831
832
833
834
835
836
837
838
839
840
841G4int BesTofDigitizerEcV4::Produce_unique_identifier(G4int module_mrpc_f, G4int readoutstripnumber_f)
842{
843 G4int returnvalue=0;
844
845
846 returnvalue= module_mrpc_f*25 + readoutstripnumber_f;
847
848 return returnvalue;
849}
850
851
852
853
854G4int BesTofDigitizerEcV4::Get_module_mrpc_from_unique_identifier(G4int unique_identifier_f)
855{
856
857 G4int returnvalue=0;
858
859 returnvalue= unique_identifier_f/25;
860
861
862 return returnvalue;
863}
864
865
866G4int BesTofDigitizerEcV4::Get_stripnumber_from_unique_identifier(G4int unique_identifier_f)
867{
868
869 G4int returnvalue=0;
870
871 returnvalue= unique_identifier_f%25;
872
873 return returnvalue;
874}
875
876
877
878
879G4int BesTofDigitizerEcV4::Get_gapnumber(G4int partID, G4double z_position)
880{
881
882
883
884 G4double pos_west_1=-1343.5*mm; // = partID 5
885 G4double pos_west_2=-1368.5*mm; // = partID 6
886 G4double pos_east_1=+1343.5*mm; // = partID 4
887 G4double pos_east_2=+1368.5*mm; // = partID 3
888
889
890 //We transform our z-coordinate into the system of one MRPC module. PoV= Interactionpoint:
891 //The gaps 1-6 will have negative values, the gaps 7-12 will have positive values!
892 G4double z_zero=0;
893 if(partID==3)
894 {z_zero= z_position-pos_east_2;}
895 else if(partID==4)
896 {z_zero= z_position-pos_east_1;}
897 else if(partID==5)
898 {z_zero = pos_west_1 - z_position ;}//z_position is negative in this case!!!!
899 else if(partID==6)
900 {z_zero = pos_west_2 - z_position ;}//z_position is negative in this case!!!!
901 else
902 {G4cout << "TofDigitzer: Error in part ID" <<G4endl; return 0;}
903
904 //std::cout << "z_pos " << z_position << std::endl;
905 //std::cout << "z_zero " << std::setprecision(100) <<z_zero << std::endl;
906
907
908 //There is a porblem with the numbers. The correct ones are the round up and down ones. However in this case to often 0 will be returned!
909 if((-4.5999*mm) >= z_zero && z_zero>(-4.8201*mm) ) {return 1;}
910 else if((-3.9799*mm) >= z_zero && z_zero>(-4.2001*mm) ){return 2;}
911 else if((-3.3599*mm) >= z_zero && z_zero>(-3.5801*mm) ){return 3;}
912 else if((-2.7399*mm) >= z_zero && z_zero>(-2.9601*mm) ){return 4;}
913 else if((-2.1199*mm) >= z_zero && z_zero>(-2.3401*mm) ){return 5;}
914 else if((-1.4999*mm) >= z_zero && z_zero>(-1.7201*mm)) {return 6;}
915 else if((1.4999*mm) <= z_zero && z_zero<=(1.7201*mm) ) {return 7;}
916 else if((2.1199*mm) <= z_zero && z_zero<=(2.3401*mm) ) {return 8;}
917 else if((2.7399*mm) <= z_zero && z_zero<=(2.9601*mm) ) {return 9;}
918 else if((3.3599*mm) <= z_zero && z_zero<=(3.5801*mm) ) {return 10;}
919 else if((3.9799*mm) <= z_zero && z_zero<=(4.2001*mm) ) {return 11;}
920 else if((4.5999*mm) <= z_zero && z_zero<=(4.8201*mm) ){return 12;}
921 else {G4cout << "TofDigitzer: Error in Hit-Position! " <<G4endl; return 0; }
922
923}
924
925
926
927G4double BesTofDigitizerEcV4::Distance_for_avalanche(G4int gap, G4double z_position, G4int partID)
928{
929
930 //All the values are in mm
931 G4double pos_west_1=-1343.5*mm; // = partID 5
932 G4double pos_west_2=-1368.5*mm; // = partID 6
933 G4double pos_east_1=+1343.5*mm; // = partID 4
934 G4double pos_east_2=+1368.5*mm; // = partID 3
935
936
937 //We transform our z-coordinate into the system of one MRPC module. PoV= Interactionpoint:
938 //The gaps 1-6 will have negative values, the gaps 7-12 will have positive values!
939 G4double z_zero=0;
940 if(partID==3)
941 {z_zero= z_position-pos_east_2;}
942 else if(partID==4)
943 {z_zero= z_position-pos_east_1;}
944 else if(partID==5)
945 {z_zero = pos_west_1 - z_position ;}//z_position is negative in this case!!!!
946 else if(partID==6)
947 {z_zero = pos_west_2 - z_position ;}//z_position is negative in this case!!!!
948 else
949 {G4cout << "TofDigitzer: Error in part ID" <<G4endl; return 0;}
950
951 //cout << "z_zero "<< setprecision(100)<< z_zero << G4endl;
952
953 /*
954 if(gap==1) G4cout << (-z_zero-4.60*mm) << " "<< gap << G4endl ;
955 if(gap==2) G4cout <<(-z_zero-3.98*mm)<< " "<< gap << G4endl;
956 if(gap==3) G4cout<< (-z_zero-3.36*mm)<< " "<< gap << G4endl;
957 if(gap==4) G4cout<< (-z_zero-2.74*mm)<< " "<< gap << G4endl;
958 if(gap==5) G4cout<< (-z_zero-2.12*mm)<< " "<< gap << G4endl;
959 if(gap==6) G4cout<< (-z_zero-1.5*mm)<< " "<< gap << G4endl;
960 if(gap==7) G4cout<< (1.72-z_zero*mm)<< " "<< gap << G4endl;
961 if(gap==8) G4cout<< (2.34-z_zero*mm)<< " "<< gap << G4endl;
962 if(gap==9) G4cout<< (2.96-z_zero*mm)<< " "<< gap << G4endl;
963 if(gap==10)G4cout<< (3.58-z_zero*mm)<< " "<< gap << G4endl;
964 if(gap==11)G4cout<< (4.20-z_zero*mm)<< " "<< gap << G4endl;
965 if(gap==12)G4cout<< (4.82-z_zero*mm)<< " "<< gap << G4endl;
966 */
967
968
969
970
971
972 if(gap==1) return (-z_zero-4.60*mm);
973 if(gap==2) return (-z_zero-3.98*mm);
974 if(gap==3) return (-z_zero-3.36*mm);
975 if(gap==4) return (-z_zero-2.74*mm);
976 if(gap==5) return (-z_zero-2.12*mm);
977 if(gap==6) return (-z_zero-1.5*mm);
978 if(gap==7) return (1.72-z_zero*mm);
979 if(gap==8) return (2.34-z_zero*mm);
980 if(gap==9) return (2.96-z_zero*mm);
981 if(gap==10)return (3.58-z_zero*mm);
982 if(gap==11)return (4.20-z_zero*mm);
983 if(gap==12)return (4.82-z_zero*mm);
984
985
986}
987
988//////////////////////////////////////////////////////////
989//////////////////////////////////////////////////////////
990//////
991////// Function used for avalanche simulation ////
992//////
993/////////////////////////////////////////////////////////
994/////////////////////////////////////////////////////////
995
996
997
998void BesTofDigitizerEcV4::simulate_avalanche(G4int actual_strip, vector<adc_info> * avalanche_info, double * avalanche_threshold_time, double *induced_charge)
999{
1000
1001 //Originally this function was a stand alone program. That's why the style of programming may be a little strange:
1002 //It depends on a method described in the Dissertation of Christian Lippmann and is known in there as 1-D Model.
1003 //Simulation with this model have been done for a chamber of the same type as implemented: Simulation study on the operation of a multi-gap resistive plate chamber - Meas. Sci. Technol. 17 (2006) 123-127
1004 //
1005 // This function use SI values as input. However it returns charge in picoColomb. The time is returned in seconds!
1006 //
1007 //General sequence of the programm:
1008 //1. Initialize all the values: Distribut the avalanche along the gap and set the number of primary electrons.
1009 //2. Invoke stepping action to simulate the avalanche evolution through the gap.
1010 //3. Construct the signal
1011
1012
1013
1014 int number_of_gaps=12;
1015 double alpha=144800.; // Townsend coefficient: /m
1016 double eta=5013; /// Attachment coefficient: /m
1017
1018 double thick_gap = 0.00022; // m //gap thickness
1019 const int steps =200;//Into how many steps is the gap divided?
1020 double user_stepwidth=thick_gap/(double)steps;
1021 double v_drift = 210.9e3; // m/s
1022
1023 double E_weight = 1./(6.*0.22e-3+(5.*0.4e-3+2*0.55e-3)/(8.)+2.*0.07e-3/(3.1)); //epsilon_glass= 8 / epsilon_mylar=3.1 /d_mylar=0.07e-3 , d_g= 0.22e-3 ,
1024 double physics_e= 1.602176462e-19;
1025 double threshold=93622.7;//Threshold in "number of electrons" when the signal will be detected: value_in_fc / physics_e
1026
1027 //double alpha=76750; double eta=8152; double v_drift=163.1e3; //85 kV/cm
1028 //double alpha=82060; double eta=7728; double v_drift=167.3e3; //87.5 kV/cm
1029 //double alpha=87560; double eta=7358; double v_drift=171.1e3;//90 kv/cm
1030 //double alpha=99020; double eta=6725; double v_drift=179.7e3;//95 kv/cm
1031 //double alpha=115200; double eta=6033; double v_drift=191.2e3;//102 kv/cm
1032 //double alpha=124700; double eta=5683; double v_drift=197.8e3;//106 kv/cm
1033 //double alpha=134600; double eta=5374; double v_drift=204.3e3;//110 kv/cm
1034 //double alpha=144500; double eta=5028; double v_drift=210.9e3;//114kv/cm
1035 //double alpha=154700; double eta=4823; double v_drift=217.5e3;//118kv/cm
1036
1037
1038
1039
1040//Intialize
1041 vector<int> primary_clusters_in_gap(number_of_gaps);
1042 for(int i=0;i<number_of_gaps;i++)primary_clusters_in_gap[i]=0;
1043
1044 vector< vector< vector<double> > > step_avalanche(number_of_gaps);
1045 vector< vector< vector<long int> > > elec_avalche(number_of_gaps);
1046
1047 //How many cluster do we have for the different steps:
1048 for(int i=0;i<number_of_gaps;i++)
1049 {
1050 int counter_avalanches=0;
1051 for(int j=0;j<avalanche_info->size();j++)
1052 {
1053 if(actual_strip==(*avalanche_info)[j].readoutstrip && ((i+1)==(*avalanche_info)[j].gap) )
1054 {counter_avalanches++;}
1055 }
1056
1057 primary_clusters_in_gap[i]=counter_avalanches;
1058 step_avalanche[i].resize(primary_clusters_in_gap[i]);
1059 elec_avalche[i].resize(primary_clusters_in_gap[i]);
1060 }
1061
1062
1063 for(int j=0;j<number_of_gaps;j++)
1064 {
1065 for(int i=0;i<primary_clusters_in_gap[j];i++)
1066 {
1067 step_avalanche[j][i].resize(steps+2);
1068 elec_avalche[j][i].resize(steps+2);
1069 }
1070 }
1071
1072 for(int h=0;h<number_of_gaps;h++) {for(int i=0;i<primary_clusters_in_gap[h];i++){for(int j=0;j<steps+2;j++){step_avalanche[h][i][j]=0.;elec_avalche[h][i][j]=0;}}} //Set all arrays equal to zero
1073
1074 vector<int> avalanches_in_gap(number_of_gaps);
1075
1076
1077
1078
1079//Distribute the cluster along the gap and set the number of ions
1080 for(int h=0;h<number_of_gaps;h++)
1081 {
1082
1083
1084 int i=0;
1085 for(int j=0;j<avalanche_info->size();j++)
1086 {
1087 if(actual_strip==(*avalanche_info)[j].readoutstrip && ((h+1)==(*avalanche_info)[j].gap) )
1088 {
1089 step_avalanche[h][i][0]=thick_gap-((*avalanche_info)[j].zdistance/m);
1090 elec_avalche[h][i][0]=(*avalanche_info)[j].ions;
1091 i++;
1092 }
1093 }
1094
1095 avalanches_in_gap[h]=primary_clusters_in_gap[h];
1096
1097 }//close for h (number_of_gaps)
1098
1099 /*
1100 for(int h=0;h<number_of_gaps;h++)
1101 {
1102 for(int j=0;j<primary_clusters_in_gap[h];j++)
1103 {
1104 G4cout << "Gap " << h+1 << " Avalanche " << j+1 << " Pos in gap " << step_avalanche[h][j][0] << G4endl;
1105 }
1106 }
1107 */
1108
1109
1110 for(int h=0;h<number_of_gaps;h++)
1111 {
1112 // G4cout << "Avalacnhes in gap " << h+1 << " are " << avalanches_in_gap[h] << G4endl;
1113 long int tot_elec_in_gap=0;
1114 if(avalanches_in_gap[h]!=0)
1115 {
1116
1117 int i=0;
1118
1119
1120 vector<bool> is_saturated(avalanches_in_gap[h]);
1121 for(int e=0;e<avalanches_in_gap[h];e++) is_saturated[e]=false;
1122
1123
1124
1125 do //Steppingaction: Moves the avalanche through the gap!
1126 {
1127
1128
1129 tot_elec_in_gap=0;
1130 for(int j=0;j<primary_clusters_in_gap[h];j++)
1131 {
1132
1133
1134
1135 elec_avalche[h][j][i+1]=adcsignal_simulate_step(elec_avalche[h][j][i],user_stepwidth,alpha,eta,is_saturated[j]);
1136 step_avalanche[h][j][i+1]= step_avalanche[h][j][i] + user_stepwidth;
1137 tot_elec_in_gap=tot_elec_in_gap + elec_avalche[h][j][i+1];
1138 //if((h+1)==12) G4cout<< "Gap "<<h+1<<" cluster "<<j+1<<" = "<< elec_avalche[h][j][i+1]<<" Pos "<<step_avalanche[h][j][i+1]<< " tot " <<tot_elec_in_gap<<G4endl;
1139
1140 if(elec_avalche[h][j][i+1]>1.5e7){is_saturated[j]=true;}else{is_saturated[j]=false;}
1141 }
1142
1143
1144 i++;
1145
1146
1147 for(int j=0;j<primary_clusters_in_gap[h];j++)
1148 {
1149
1150 if((step_avalanche[h][j][i])>=thick_gap)
1151 {
1152 primary_clusters_in_gap[h]=(primary_clusters_in_gap[h]-1);
1153 }
1154 }
1155
1156
1157 }
1158 while(step_avalanche[h][0][i]<thick_gap); //If the first avalanche leaves the gap, quit the stepping action
1159
1160
1161 }//Close if(avalanches_in_gap!=0)
1162 }//close for over gaps
1163
1164
1165
1166
1167//Produce signal
1168 double signal=0;//stores the total sum of the induced signal!
1169 bool over_threshold = false;
1170
1171 for(int i=0;i<steps;i++)
1172 {
1173 for(int h=0;h<number_of_gaps;h++)
1174 {
1175 for(int j=0;j<avalanches_in_gap[h];j++)
1176 {
1177 signal= signal + elec_avalche[h][j][i];
1178 }
1179 }
1180
1181 if(over_threshold == false && signal>threshold) {over_threshold=true; *avalanche_threshold_time=(thick_gap/steps/v_drift*(i+1));}
1182
1183 }//close the for loop over steps.
1184
1185
1186
1187
1188 //Return the value in pC:
1189
1190 //G4cout <<" *avalanche_threshold_time " << *avalanche_threshold_time << G4endl;
1191 //G4cout << "Induced Charge "<< signal*E_weight*v_drift*physics_e*1e12 << G4endl;
1192 *induced_charge=signal*E_weight*v_drift*physics_e*1e12 * thick_gap/steps/v_drift;
1193
1194
1195
1196}//close simulate_avalanche
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208//This fuction return the number of electrons for the next step
1209long int BesTofDigitizerEcV4::adcsignal_simulate_step(int n_elec,double stepwidth,double alpha,double eta, bool saturated)
1210{
1211
1212 if(n_elec<=150 && (saturated == false))
1213 {
1214
1215 long int n_elec_step =0;
1216
1217 for(int j=0;j<n_elec;j++)
1218 {
1219 double s=G4UniformRand();
1220
1221 n_elec_step = n_elec_step + adcsignal_get_n_electron(s, stepwidth, alpha, eta);
1222
1223 }
1224 return n_elec_step;
1225 }
1226
1227 if(n_elec>150 && saturated == false)
1228 {
1229 double sigma=adcsignal_get_sigma(alpha,eta,stepwidth);
1230 double nbar = exp((alpha-eta)*stepwidth);
1231 double my_sigma = G4RandGauss::shoot(0,(sqrt(n_elec)*sigma));
1232
1233 return (long int)(n_elec*nbar+ my_sigma);
1234 }
1235
1236 if(saturated==true)
1237 {
1238
1239 return n_elec;
1240 }//close if(is_saturated==true)
1241
1242}//closesimulate step
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252//Return the number of electrons within the next step
1253long int BesTofDigitizerEcV4::adcsignal_get_n_electron(double s, double stepwidth, double alpha, double eta)
1254{
1255
1256 double nbar = exp((alpha-eta)*stepwidth);
1257 double s_border= eta/alpha*(nbar-1.)/(nbar-eta/alpha);
1258 double k = eta/alpha;
1259
1260
1261 if(s<s_border)
1262 {
1263 return 0;
1264
1265 }
1266 else
1267 {
1268 return (long int) (1. + (1/log(((nbar-1)/(nbar-k))) * log(((nbar-k)*(s-1)/(k-1)/nbar))));
1269 }
1270
1271
1272}//close get_n_electron
1273
1274
1275double BesTofDigitizerEcV4::adcsignal_get_sigma(double alpha,double eta, double stepwidth)
1276{
1277
1278 double k=eta/alpha;
1279 double nbar = exp((alpha-eta)*stepwidth);
1280
1281 return (sqrt( (1.+k)/(1.-k) * nbar * (nbar-1.) )) ;
1282
1283}//close get_sigma
1284
1285
1286
1287
1288
1289
1290double BesTofDigitizerEcV4::chargetotime(double charge)
1291{
1292
1293
1294 double time=0.;
1295 if(charge>=13.) time=(-4.50565/log(charge*0.0812208)+16.6653)*ns;
1296
1297 return time;
1298
1299}//close chargetotime
1300
1301
1302
1303
1304
1305
1306
1307double BesTofDigitizerEcV4::Get_NINO_pulselength_resolution(double induced_charge_fc)
1308{
1309
1310 double result = (64.3326*exp(-induced_charge_fc/25.4638 + 0.944184)+19.4846 );
1311
1312 return result;
1313
1314}
1315
1316
1317
1318double BesTofDigitizerEcV4::Get_NINO_leadingedge_resolution(double induced_charge_fc)
1319{
1320
1321 double result =0;
1322
1323 if(induced_charge_fc > 50.)
1324 {
1325 result = 67.6737*exp(-induced_charge_fc/50.9995-0.27755)+9.06223;
1326 }
1327 else
1328 {
1329 result = 176.322-2.98345*induced_charge_fc;
1330 }
1331
1332 return result;
1333
1334}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
BesTofDigitsCollection
G4TDigiCollection< BesTofDigi > BesTofDigitsCollection
Definition: BesTofDigi.hh:83
BesTofHitsCollection
G4THitsCollection< BesTofHit > BesTofHitsCollection
Definition: BesTofHit.hh:108
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
time
Double_t time
Definition: DataBase/tau_mode.c:7
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
alpha
const double alpha
Definition: FastVertexFit.cxx:4
help
std::string help()
Definition: HelloServer.cpp:35
s
XmlRpcServer s
Definition: HelloServer.cpp:11
BesTofDigi::SetTime1_mrpc
void SetTime1_mrpc(G4double time_mrpc)
Definition: BesTofDigi.hh:47
BesTofDigi::SetPartId
void SetPartId(G4int partId)
Definition: BesTofDigi.hh:39
BesTofDigi::SetTrackIndex
void SetTrackIndex(G4int index)
Definition: BesTofDigi.hh:38
BesTofDigi::SetTime2_mrpc
void SetTime2_mrpc(G4double time_mrpc2)
Definition: BesTofDigi.hh:48
BesTofDigi::SetModule_mrpc
void SetModule_mrpc(G4int module_mrpc)
Definition: BesTofDigi.hh:46
BesTofDigitizerEcV4::adcsignal_get_n_electron
long int adcsignal_get_n_electron(double s, double stepwidth, double alpha, double eta)
Definition: BesTofDigitizerEcV4.cc:1253
BesTofDigitizerEcV4::Average_transition_time
static G4double Average_transition_time(G4int my_module)
Definition: BesTofDigitizerEcV4.cc:810
BesTofDigitizerEcV4::Calculate_strip_transition_time
static G4double Calculate_strip_transition_time(G4double x_mm, G4double y_mm, G4int partId_f, G4int module_mrpc_f)
Definition: BesTofDigitizerEcV4.cc:534
BesTofDigitizerEcV4::Distance_for_avalanche
G4double Distance_for_avalanche(G4int gap, G4double z_position, G4int partID)
Definition: BesTofDigitizerEcV4.cc:927
BesTofDigitizerEcV4::Digitize
virtual void Digitize(ScintSingle *, BesTofDigitsCollection *)
Definition: BesTofDigitizerEcV4.cc:75
BesTofDigitizerEcV4::Produce_unique_identifier
static G4int Produce_unique_identifier(G4int module_mrpc_f, G4int readoutstripnumber_f)
Definition: BesTofDigitizerEcV4.cc:841
BesTofDigitizerEcV4::Get_NINO_pulselength_resolution
double Get_NINO_pulselength_resolution(double induced_charge_fc)
Definition: BesTofDigitizerEcV4.cc:1307
BesTofDigitizerEcV4::simulate_avalanche
void simulate_avalanche(G4int actual_strip, vector< adc_info > *avalanche_info, double *avalanche_threshold_time, double *induced_charge)
Definition: BesTofDigitizerEcV4.cc:998
BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum
static G4int Calculate_Readoutstrip_number_continuum(G4double x_mm, G4double y_mm, G4int partId_f, G4int module_mrpc_f)
Definition: BesTofDigitizerEcV4.cc:455
BesTofDigitizerEcV4::Get_module_mrpc_from_unique_identifier
static G4int Get_module_mrpc_from_unique_identifier(G4int unique_identifier_f)
Definition: BesTofDigitizerEcV4.cc:854
BesTofDigitizerEcV4::Calculate_strip_transition_time_cont
static G4double Calculate_strip_transition_time_cont(G4double x_mm, G4double y_mm, G4int partId_f, G4int module_mrpc_f)
Definition: BesTofDigitizerEcV4.cc:692
BesTofDigitizerEcV4::Calculate_Readoutstrip_number
static G4int Calculate_Readoutstrip_number(G4double x, G4double y, G4int partId_f, G4int module_mrpc_f)
Definition: BesTofDigitizerEcV4.cc:371
BesTofDigitizerEcV4::Get_NINO_leadingedge_resolution
double Get_NINO_leadingedge_resolution(double induced_charge_fc)
Definition: BesTofDigitizerEcV4.cc:1318
BesTofDigitizerEcV4::adcsignal_get_sigma
double adcsignal_get_sigma(double alpha, double eta, double stepwidth)
Definition: BesTofDigitizerEcV4.cc:1275
BesTofDigitizerEcV4::Smear_gaussian
G4double Smear_gaussian(G4double, G4double, G4double)
Definition: BesTofDigitizerEcV4.cc:359
BesTofDigitizerEcV4::Get_gapnumber
G4int Get_gapnumber(G4int partID, G4double z_position)
Definition: BesTofDigitizerEcV4.cc:879
BesTofDigitizerEcV4::Get_stripnumber_from_unique_identifier
static G4int Get_stripnumber_from_unique_identifier(G4int unique_identifier_f)
Definition: BesTofDigitizerEcV4.cc:866
BesTofDigitizerEcV4::adcsignal_simulate_step
long int adcsignal_simulate_step(int n_elec, double stepwidth, double alpha, double eta, bool saturated)
Definition: BesTofDigitizerEcV4.cc:1209
BesTofDigitizerV::m_besTofDigitsCollection
BesTofDigitsCollection * m_besTofDigitsCollection
Definition: BesTofDigitizerV.hh:38
BesTofHit::GetIons
G4int GetIons()
Definition: BesTofHit.hh:80
BesTofHit::GetPDGcode
G4int GetPDGcode()
Definition: BesTofHit.hh:73
BesTofHit::GetTime
G4double GetTime()
Definition: BesTofHit.hh:66
BesTofHit::GetPos
G4ThreeVector GetPos()
Definition: BesTofHit.hh:65
BesTofHit::GetG4Index
G4int GetG4Index()
Definition: BesTofHit.hh:59
ScintSingle::GetHitIndexes_mrpc
vector< G4int > * GetHitIndexes_mrpc()
Definition: ScintSingle.hh:48
ScintSingle::GetModule_mrpc
G4int GetModule_mrpc()
Definition: ScintSingle.hh:49
ScintSingle::GetPartId
G4int GetPartId()
Definition: ScintSingle.hh:44
BesTofDigitizerEcV4::adc_info
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::trackindex
G4int trackindex
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::time
G4double time
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::y
G4double y
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::readoutstrip
G4int readoutstrip
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::gap
G4int gap
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::x
G4double x
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::zdistance
G4double zdistance
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::transtime
G4double transtime
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::pdg_code
G4int pdg_code
Definition: BesTofDigitizerEcV4.hh:51
BesTofDigitizerEcV4::adc_info::ions
G4int ions
Definition: BesTofDigitizerEcV4.hh:51
pi
const float pi
Definition: vector3.h:133
ns
#define ns(x)
Definition: xmltok.c:1504