GEMC  2.3
Geant4 Monte-Carlo Framework
veto_hitprocess.cc
Go to the documentation of this file.
1 
2 // G4 headers
3 #include "G4Poisson.hh"
4 #include "Randomize.hh"
5 
6 // gemc headers
7 #include "veto_hitprocess.h"
8 
9 // CLHEP units
10 #include "CLHEP/Units/PhysicalConstants.h"
11 using namespace CLHEP;
12 
13 map<string, double> veto_HitProcess :: integrateDgt(MHit* aHit, int hitn)
14 {
15  map<string, double> dgtz;
16  vector<identifier> identity = aHit->GetId();
17 
18  int sector = identity[0].id;
19  int veto_id = identity[1].id;
20  int channel = identity[2].id;
21 
22  // Digitization Parameters
23 
24 
25 
26 // double adc_conv=10; // conversion factor from pC to ADC (typical sensitivy of CAEN VME QDC is of 0.1 pC/ch)
27 // double adc_ped=0; // ADC Pedestal
28 // double tdc_conv=1*ns; // TDC conversion factor
29 
30 
31  // initialize ADC and TDC
32  double etotL = 0;
33  double etotR = 0;
34  double timeL = 0;
35  double timeR = 0;
36  int ADC1 = 0;
37  int ADC2 = 0;
38  int ADC3 = 0;
39  int ADC4 = 0;
40  int TDC1 = 4096;
41  int TDC2 = 4096;
42  int TDC3 = 4096;
43  int TDC4 = 4096;
44 
45  double length = 0;
46  // From measurement
47  // spe = 0.36pC
48  // Cosmic = 84pC (~235pe)
49  // Attenuation at 80cm: ~0.85 -> effective att lenght ~350cm
50  //Plastic
51  double paddle_surface = 0; // paddle surface
52  double light_yield = 0; // number of optical photons pruced in the scintillator per MeV of deposited energy
53  double att_length = 0; // light at tenuation length
54  double veff = 0; // light velocity in scintillator
55  //PMT
56  double sensor_surface = 0; // area of photo sensor
57  double sensor_effective_area = 0; // considering only a fraction of the photocathod
58  double sensor_qe = 0; // photo sensor quantum efficiency
59  double sensor_gain = 0; // pmt gain x electron charge in pC (2.2x10^6)x(1.6x10^-7) -> ~0.36pC or 1 to have pe
60  double light_coll = 0; // ratio of photo_sensor area over paddle section ~ light collection efficiency
61  double light_guide_att = 0;
62  double tL = 0;
63  double tR = 0;
64  double peL=0.;
65  double peR = 0;
66  double etot_g4=0.;
67 
68  // Proposal
69  // sector: run over channels, from 0 to N
70  // Oveto == 5
71  // channel: run over the position (1=T 2=B 3=R 4=L 5=D 6=U)
72  if(veto_id==5)
73  {
74  {
75  double optical_coupling[13]= {0., 0.94,0.57, 0.35, 0.7, 0.094, 0.177, 0.52, 0.75, 0.52, 0.52, 0.38, 1.0 };
76  for (int s=0; s<13; s++) optical_coupling[s] = optical_coupling[s]*0.68;
77  light_yield=9200/MeV;
78  veff=13*cm/ns ;
79  sensor_effective_area=0.9;
80  sensor_qe=0.25;
81  sensor_gain=1.;
82 
83  // Upper/lower
84  if(channel==1 || channel==2)
85  {
86  // Get the paddle length: in veto paddles are along z
87  length = aHit->GetDetector().dimensions[2];
88  double s1=aHit->GetDetector().dimensions[0];
89  double s2=aHit->GetDetector().dimensions[1];
90  paddle_surface = 2*s1*2*s2;
91  sensor_surface=pow(2.5*cm,2)*pi; // 2" pmt R-> (2*2.5/2 TBC)
92  att_length=350*cm;
93  light_guide_att=1.0;
94  // cout << " lenght: " << length << " optical-coupled surface: " << paddle_surface <<endl;
95  }
96  if(channel==5 || channel==6)
97  {
98  double s1=aHit->GetDetector().dimensions[0];
99  double s2=aHit->GetDetector().dimensions[3];
100  paddle_surface = 2*s1*2*s2;// surface perpendicular to the pmt position Surf=XxZ
101  sensor_surface=pow((2.5/2)*cm,2)*pi; //
102  att_length=400*cm; // longer att lenght to take into account the perpendicular readout
103  light_guide_att=0.19; // no light guide
104  // cout << " lenght: " << length << " optical-coupled surface: " << paddle_surface <<endl;
105  }
106  if(channel==3 || channel==4)
107  {
108  // Get the paddle length: in veto paddles are along y
109  length = aHit->GetDetector().dimensions[1];
110  double s1=aHit->GetDetector().dimensions[0];
111  double s2=aHit->GetDetector().dimensions[2];
112  sensor_surface=pow(2.5*cm,2)*pi; // 2" pmt R-> (2*2.5/2 TBC)
113  paddle_surface = 2*s1*2*s2;
114  att_length=350*cm;
115  light_guide_att=1.0;
116  // cout << " lenght: " << length << " optical-coupled surface: " << paddle_surface <<endl;
117  }
118 
119  light_coll=sensor_surface/paddle_surface;
120  if (sensor_surface>paddle_surface) light_coll=1.; // no more than the PMT size
121  light_coll=optical_coupling[1]*light_coll*sensor_effective_area*light_guide_att; // coupling [1] identical for all
122  //cout << " light collo " << light_coll <<endl;
123 
124  // Get info about detector material to eveluate Birks effect
125  double birks_constant=aHit->GetDetector().GetLogical()->GetMaterial()->GetIonisation()->GetBirksConstant();
126  // cout << " Birks constant is: " << birks_constant << endl;
127  // cout << aHit->GetDetector().GetLogical()->GetMaterial()->GetName() << endl;
128 
129 
130  double time_min[4] = {0,0,0,0};
131 
132  vector<G4ThreeVector> Lpos = aHit->GetLPos();
133  vector<G4double> Edep = aHit->GetEdep();
134  vector<G4double> Dx = aHit->GetDx();
135  // Charge for each step
136  vector<int> charge = aHit->GetCharges();
137  vector<G4double> times = aHit->GetTime();
138 
139  unsigned int nsteps = Edep.size();
140  double Etot = 0;
141 
142  for(unsigned int s=0; s<nsteps; s++) Etot = Etot + Edep[s];
143 
144  if(Etot>0)
145  {
146  for(unsigned int s=0; s<nsteps; s++)
147  {
148  double dLeft =-10000.;
149  double dRight =-10000.;
150 
151  // Distances from left, right for upper/lower (along z)
152  if(channel==1 || channel==2)
153  {
154  dLeft = length - Lpos[s].z();
155  dRight = length + Lpos[s].z();
156  }
157  // Distances from left, right for other OV (along y)
158 
159  if(channel==3 || channel==4)
160  {
161  dLeft = length + Lpos[s].y();
162  dRight = length - Lpos[s].y();
163  }
164  if(channel==5 || channel==6)
165  {
166  dLeft = Lpos[s].x();
167  dRight = Lpos[s].y();
168  double dCent=sqrt(dLeft*dLeft+dRight*dRight);
169  dRight =dCent;
170 
171  }
172 
173 
174 
175  // cout << "\n Distances: " << endl;
176  // cout << "\t dLeft, dRight " << dLeft << ", " << dRight << endl;
177 
178  // apply Birks effect
179  // double stepl = 0.;
180 
181  // if (s == 0){
182  // stepl = sqrt(pow((Lpos[s+1].x() - Lpos[s].x()),2) + pow((Lpos[s+1].y() - Lpos[s].y()),2) + pow((Lpos[s+1].z() - Lpos[s].z()),2));
183  // }
184  // else {
185  // stepl = sqrt(pow((Lpos[s].x() - Lpos[s-1].x()),2) + pow((Lpos[s].y() - Lpos[s-1].y()),2) + pow((Lpos[s].z() - Lpos[s-1].z()),2));
186  // }
187 
188  double Edep_B = BirksAttenuation(Edep[s],Dx[s],charge[s],birks_constant);
189  Edep_B=Edep[s];
190  etot_g4=etot_g4+Edep_B;
191  // cout << "\t Birks Effect: " << " Edep=" << Edep[s] << " StepL=" << stepl
192  // << " PID =" << pids[s] << " charge =" << charge[s] << " Edep_B=" << Edep_B << endl;
193 
194  //if (light_coll > 1) light_coll = 1.; // To make sure you don't miraculously get more energy than you started with
195 
196  etotL = etotL + Edep_B/2 * exp(-dLeft/att_length) * light_coll;
197  etotR = etotR + Edep_B/2 * exp(-dRight/att_length) * light_coll;
198 
199  // cout << "step: " << s << " etotL, etotR " << etotL << ", " << etotR << endl;
200 
201  timeL= timeL + (times[s] + dLeft/veff) / nsteps;
202  timeR= timeR + (times[s] + dRight/veff) / nsteps;
203 
204  if(etotL > 0.) {
205  if(s==0 || (time_min[0]>(times[s]+dLeft/veff))) time_min[0]=times[s]+dLeft/veff;
206  }
207  // cout << "min " << time_min[0] << "min " << times[s]+dLeft/veff << endl;
208  if(etotR > 0.) {
209  if(s==0 || (time_min[1]>(times[s]+dRight/veff))) time_min[1]=times[s]+ dRight/veff;
210  }
211  }
212  //cout << " etotR " << etotR << " ; etotL" << etotL <<endl;
213  peL=G4Poisson(etotL*light_yield*sensor_qe);
214  peR=G4Poisson(etotR*light_yield*sensor_qe);
215  //cout << " per " << peR << " ; pel" << peL <<endl;
216  //peL=(etotL*light_yield*sensor_qe);
217  //peR=(etotR*light_yield*sensor_qe);
218  //cout << " per " << peR << " ; pel" << peL <<endl;
219 
220  double sigmaTL=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peL+1.));
221  double sigmaTR=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peR+1.));
222  //sigmaTL=0;
223  //sigmaTR=0;
224  tL=(time_min[0]+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
225  tR=(time_min[1]+G4RandGauss::shoot(0.,sigmaTR))*1000.;// time in ps
226  // Digitization for ADC and QDC not used
227  //TDC1=(int) (tL * tdc_conv);
228  //TDC2=(int) (tR * tdc_conv);
229  //if(TDC1<0) TDC1=0;
230  //if(TDC2<0) TDC2=0;
231  //ADC1=(int) (peL*sensor_gain*adc_conv + adc_ped);
232  //ADC2=(int) (peR*sensor_gain*adc_conv + adc_ped);
233 
234  // cout << "ADC1: " << ADC1 << " " << peL << " " << sensor_gain << " " << adc_conv << endl;
235 
236 
237  //cout << "energy right: " << ADC2 / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E left: " << ADC1 / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
238  //cout << "energy forw: " << ADCF / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E back: " << ADCB / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
239 
240  //cout << " Light collection: " << light_coll << endl;
241 
242  }
243  // closes (Etot > 0) loop
244 
245 
246 
247  // ch 1,3 dRight
248  // ch 2,4 dLeft
249  // ch 7,8,9,10,11,12 dRight
250  // ch 5,6 dRight
251  // ignore the other side
252 
253  if(channel==1 || channel==2 )
254  { if(sector==0)
255  {ADC1=peR;
256  TDC1=tR;}
257  else if(sector==1)
258  {ADC1=peL;
259  TDC1=tL;}
260  }
261  if(channel==3 || channel ==4)
262  {
263  ADC1=peR;
264  TDC1=tR;
265  }
266  if(channel==5 || channel==6)
267  {
268  ADC1=peR;
269  TDC1=tR;
270  }
271 
272  // cout << " ADC1: " << ADC1 << " ; TDC1: " << TDC1 << " ; ADC2: "<< etot_g4*1000. << endl;
273  if(verbosity>4)
274  {
275  cout << log_msg << " veto: " << veto_id << ", channel: " << channel << ", sector: " << sector ;
276  cout << log_msg << " Etot=" << Etot/MeV << endl;
277  cout << log_msg << " TDC1=" << TDC1 << " TDC2=" << TDC2 << " ADC1=" << ADC1 << " ADC2=" << ADC2 << endl;
278  }
279 
280 
281  }
282  }
283  else if(veto_id==4)
284  // proposal IV
285  {
286  double veff=13*cm/ns ;// TO BE CHECKED
287  // scintillator sizes
288  double sx=aHit->GetDetector().dimensions[0];
289  double sy=aHit->GetDetector().dimensions[1];
290  double sz=aHit->GetDetector().dimensions[2];
291 
292 // double time_min[4] = {0,0,0,0};
293 
294  vector<G4ThreeVector> Lpos = aHit->GetLPos();
295  vector<G4double> Edep = aHit->GetEdep();
296  vector<G4double> Dx = aHit->GetDx();
297  // Charge for each step
298  vector<int> charge = aHit->GetCharges();
299  vector<G4double> times = aHit->GetTime();
300  unsigned int nsteps = Edep.size();
301  double Etot = 0;
302  double X_hit_ave=0.;
303  double Y_hit_ave=0.;
304  double Z_hit_ave=0.;
305  double T_hit_ave=0.;
306  double dLeft =-10000.;
307 
308  for(unsigned int s=0; s<nsteps; s++) Etot = Etot + Edep[s];
309  if(Etot>0)
310  {
311  for(unsigned int s=0; s<nsteps; s++)
312  {
313  double Edep_B=Edep[s];
314  etot_g4=etot_g4+Edep_B;
315  // average hit position XYZ
316  X_hit_ave=X_hit_ave+Lpos[s].x();
317  Y_hit_ave=Y_hit_ave+Lpos[s].y();
318  Z_hit_ave=Z_hit_ave+Lpos[s].z();
319  // average hit time
320  T_hit_ave=T_hit_ave+times[s];
321 
322  //cout << "X " << Lpos[s].x() << " " << "Y " << Lpos[s].y() << " " << "Z " << Lpos[s].z() << " "<< "T " <<times[s] << " " << endl;
323 
324  }
325  X_hit_ave=X_hit_ave/nsteps;
326  Y_hit_ave=Y_hit_ave/nsteps;
327  Z_hit_ave=Z_hit_ave/nsteps;
328  T_hit_ave=T_hit_ave/nsteps;
329  dLeft =sz-Z_hit_ave;
330  timeL= dLeft/veff+T_hit_ave;
331 
332 
333  double *pe_sipm;// response for a mip (2..05 MeV energy released in 1cm thick)
334 
335  pe_sipm=IVresponseProposal(channel, X_hit_ave, Y_hit_ave, Z_hit_ave,sx,sy,sz);
336 
337  ADC1=G4Poisson(pe_sipm[0]*etot_g4/2.05) ; // Scaling for more/less energy release)
338  ADC2=G4Poisson(pe_sipm[1]*etot_g4/2.05) ; // Scaling for more/less energy release)
339  ADC3=G4Poisson(pe_sipm[2]*etot_g4/2.05) ; // Scaling for more/less energy release)
340  ADC4=G4Poisson(pe_sipm[3]*etot_g4/2.05) ; // Scaling for more/less energy release)
341  double sigmaTL=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peL+1.));
342  sigmaTL=0.;
343  TDC1=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
344  TDC2=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
345  TDC3=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
346  TDC4=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
347 
348  //cout << log_msg << " veto: " << veto_id << ", channel: " << channel << ", sector: " << sector << endl;
349  //cout << "X " << X_hit_ave << " " << "Y " << Y_hit_ave << " " << "Z " << Z_hit_ave << " "<< "T " <<T_hit_ave << " " << endl;
350  //cout << "sipm1 " << pe_sipm[0] << " " << "sipm2 " << pe_sipm[1] << " " << "sipm3 " << pe_sipm[2] << " " << "sipm4 " << pe_sipm[3] << " " << endl;
351  // cout << "dLeft " << dLeft << " " << "timeL " << timeL << " " << endl;
352 
353  //cout << "energy right: " << ADC2 / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E left: " << ADC1 / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
354  //cout << "energy forw: " << ADCF / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E back: " << ADCB / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
355 
356  //cout << " Light collection: " << light_coll << endl;
357 
358  }
359  // closes (Etot > 0) loop
360 
361  }
362 
363  // End Proposal
364 
365 
366 
367 
368  // Outer veto
369  if(veto_id==2)
370  {
371  double optical_coupling[13]= {0., 0.94,0.57, 0.35, 0.7, 0.094, 0.177, 0.52, 0.75, 0.52, 0.52, 0.38, 1.0 };
372  for (int s=0; s<13; s++) optical_coupling[s] = optical_coupling[s]*0.68;
373  light_yield=9200/MeV;
374  veff=13*cm/ns ;
375  sensor_effective_area=0.9;
376  sensor_qe=0.25;
377  sensor_gain=1.;
378 
379  // Upper/lower
380  if(channel==1 || channel==2 || channel==3 || channel ==4)
381  {
382  // Get the paddle length: in veto paddles are along z
383  length = aHit->GetDetector().dimensions[2];
384  double s1=aHit->GetDetector().dimensions[0];
385  double s2=aHit->GetDetector().dimensions[1];
386  paddle_surface = 2*s1*2*s2;
387  sensor_surface=pow(2.5*cm,2)*pi; // 2" pmt R-> (2*2.5/2 TBC)
388  att_length=350*cm;
389  light_guide_att=1.0;
390  // cout << " lenght: " << length << " optical-coupled surface: " << paddle_surface <<endl;
391  }
392  if(channel==5 || channel==6)
393  {
394  double s1=aHit->GetDetector().dimensions[0];
395  double s2=aHit->GetDetector().dimensions[3];
396  paddle_surface = 2*s1*2*s2;// surface perpendicular to the pmt position Surf=XxZ
397  sensor_surface=pow((2.5/2)*cm,2)*pi; //
398  att_length=400*cm; // longer att lenght to take into account the perpendicular readout
399  light_guide_att=0.19; // no light guide
400  // cout << " lenght: " << length << " optical-coupled surface: " << paddle_surface <<endl;
401  }
402  if(channel==7 || channel==8 || channel==9 || channel ==10 || channel ==11 || channel ==12 )
403  {
404  // Get the paddle length: in veto paddles are along y
405  length = aHit->GetDetector().dimensions[1];
406  double s1=aHit->GetDetector().dimensions[0];
407  double s2=aHit->GetDetector().dimensions[2];
408  sensor_surface=pow(2.5*cm,2)*pi; // 2" pmt R-> (2*2.5/2 TBC)
409  paddle_surface = 2*s1*2*s2;
410  att_length=350*cm;
411  light_guide_att=1.0;
412  // cout << " lenght: " << length << " optical-coupled surface: " << paddle_surface <<endl;
413  }
414 
415  light_coll=sensor_surface/paddle_surface;
416  if (sensor_surface>paddle_surface) light_coll=1.; // no more than the PMT size
417  light_coll=optical_coupling[channel]*light_coll*sensor_effective_area*light_guide_att; // Including the coupling efficiency and the pc effective area
418  //cout << " light collo " << light_coll <<endl;
419 
420  // Get info about detector material to eveluate Birks effect
421  double birks_constant=aHit->GetDetector().GetLogical()->GetMaterial()->GetIonisation()->GetBirksConstant();
422  // cout << " Birks constant is: " << birks_constant << endl;
423  // cout << aHit->GetDetector().GetLogical()->GetMaterial()->GetName() << endl;
424 
425 
426  double time_min[4] = {0,0,0,0};
427 
428  vector<G4ThreeVector> Lpos = aHit->GetLPos();
429  vector<G4double> Edep = aHit->GetEdep();
430  vector<G4double> Dx = aHit->GetDx();
431  // Charge for each step
432  vector<int> charge = aHit->GetCharges();
433  vector<G4double> times = aHit->GetTime();
434 
435  unsigned int nsteps = Edep.size();
436  double Etot = 0;
437 
438  for(unsigned int s=0; s<nsteps; s++) Etot = Etot + Edep[s];
439 
440  if(Etot>0)
441  {
442  for(unsigned int s=0; s<nsteps; s++)
443  {
444  double dLeft =-10000.;
445  double dRight =-10000.;
446 
447  // Distances from left, right for upper/lower (along z)
448  if(channel==1 || channel==2 || channel==3 || channel ==4)
449  {
450  dLeft = length - Lpos[s].z();
451  dRight = length + Lpos[s].z();
452  }
453  // Distances from left, right for other OV (along y)
454 
455  if(channel==7 || channel==8 || channel==9 || channel ==10 || channel ==11 || channel ==12 )
456  {
457  dLeft = length + Lpos[s].y();
458  dRight = length - Lpos[s].y();
459  }
460  if(channel==5 || channel==6)
461  {
462  dLeft = Lpos[s].x();
463  dRight = Lpos[s].y();
464  double dCent=sqrt(dLeft*dLeft+dRight*dRight);
465  dRight =dCent;
466 
467  }
468 
469 
470 
471  // cout << "\n Distances: " << endl;
472  // cout << "\t dLeft, dRight " << dLeft << ", " << dRight << endl;
473 
474  // apply Birks effect
475 // double stepl = 0.;
476 
477 // if (s == 0){
478 // stepl = sqrt(pow((Lpos[s+1].x() - Lpos[s].x()),2) + pow((Lpos[s+1].y() - Lpos[s].y()),2) + pow((Lpos[s+1].z() - Lpos[s].z()),2));
479 // }
480 // else {
481 // stepl = sqrt(pow((Lpos[s].x() - Lpos[s-1].x()),2) + pow((Lpos[s].y() - Lpos[s-1].y()),2) + pow((Lpos[s].z() - Lpos[s-1].z()),2));
482 // }
483 
484  double Edep_B = BirksAttenuation(Edep[s],Dx[s],charge[s],birks_constant);
485  Edep_B=Edep[s];
486  etot_g4=etot_g4+Edep_B;
487  // cout << "\t Birks Effect: " << " Edep=" << Edep[s] << " StepL=" << stepl
488  // << " PID =" << pids[s] << " charge =" << charge[s] << " Edep_B=" << Edep_B << endl;
489 
490  //if (light_coll > 1) light_coll = 1.; // To make sure you don't miraculously get more energy than you started with
491 
492  etotL = etotL + Edep_B/2 * exp(-dLeft/att_length) * light_coll;
493  etotR = etotR + Edep_B/2 * exp(-dRight/att_length) * light_coll;
494 
495 // cout << "step: " << s << " etotL, etotR " << etotL << ", " << etotR << endl;
496 
497  timeL= timeL + (times[s] + dLeft/veff) / nsteps;
498  timeR= timeR + (times[s] + dRight/veff) / nsteps;
499 
500  if(etotL > 0.) {
501  if(s==0 || (time_min[0]>(times[s]+dLeft/veff))) time_min[0]=times[s]+dLeft/veff;
502  }
503  // cout << "min " << time_min[0] << "min " << times[s]+dLeft/veff << endl;
504  if(etotR > 0.) {
505  if(s==0 || (time_min[1]>(times[s]+dRight/veff))) time_min[1]=times[s]+ dRight/veff;
506  }
507  }
508  //cout << " etotR " << etotR << " ; etotL" << etotL <<endl;
509  peL=G4Poisson(etotL*light_yield*sensor_qe);
510  peR=G4Poisson(etotR*light_yield*sensor_qe);
511  //cout << " per " << peR << " ; pel" << peL <<endl;
512  //peL=(etotL*light_yield*sensor_qe);
513  //peR=(etotR*light_yield*sensor_qe);
514  //cout << " per " << peR << " ; pel" << peL <<endl;
515 
516  double sigmaTL=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peL+1.));
517  double sigmaTR=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peR+1.));
518  //sigmaTL=0;
519  //sigmaTR=0;
520  tL=(time_min[0]+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
521  tR=(time_min[1]+G4RandGauss::shoot(0.,sigmaTR))*1000.;// time in ps
522  // Digitization for ADC and QDC not used
523  //TDC1=(int) (tL * tdc_conv);
524  //TDC2=(int) (tR * tdc_conv);
525  //if(TDC1<0) TDC1=0;
526  //if(TDC2<0) TDC2=0;
527  //ADC1=(int) (peL*sensor_gain*adc_conv + adc_ped);
528  //ADC2=(int) (peR*sensor_gain*adc_conv + adc_ped);
529 
530  // cout << "ADC1: " << ADC1 << " " << peL << " " << sensor_gain << " " << adc_conv << endl;
531 
532 
533  //cout << "energy right: " << ADC2 / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E left: " << ADC1 / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
534  //cout << "energy forw: " << ADCF / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E back: " << ADCB / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
535 
536  //cout << " Light collection: " << light_coll << endl;
537 
538  }
539  // closes (Etot > 0) loop
540 
541 
542  if(verbosity>4)
543  {
544  cout << log_msg << " veto: " << veto_id << ", channel: " << channel ;
545  cout << log_msg << " Etot=" << Etot/MeV << endl;
546  cout << log_msg << " TDC1=" << TDC1 << " TDC2=" << TDC2 << " ADC1=" << ADC1 << " ADC2=" << ADC2 << endl;
547  }
548 
549 
550  // ch 1,3 dRight
551  // ch 2,4 dLeft
552  // ch 7,8,9,10,11,12 dRight
553  // ch 5,6 dRight
554  // ignore the other side
555 
556  if(channel==1 || channel==3 )
557  {
558  ADC1=peR;
559  TDC1=tR;
560  }
561  if(channel==2 || channel ==4)
562  {
563  ADC1=peL;
564  TDC1=tL;
565  }
566  if(channel==5 || channel==6)
567  {
568  ADC1=peR;
569  TDC1=tR;
570  }
571  if(channel==7 || channel==8 || channel==9 || channel ==10 || channel ==11 || channel ==12 )
572  {
573  ADC1=peR;
574  TDC1=tR;
575 
576  }
577  // cout << " ADC1: " << ADC1 << " ; TDC1: " << TDC1 << " ; ADC2: "<< etot_g4*1000. << endl;
578 
579  } // end of OV
580 
581 
582 
583 
584 
585 
586 
587 // INNER VETO
588  if(veto_id==1)
589  {
590  double veff=13*cm/ns ;// TO BE CHECKED
591  // scintillator sizes
592 // double sx=aHit->GetDetector().dimensions[0];
593 // double sy=aHit->GetDetector().dimensions[1];
594  double sz=aHit->GetDetector().dimensions[2];
595 
596 // double time_min[4] = {0,0,0,0};
597 
598  vector<G4ThreeVector> Lpos = aHit->GetLPos();
599  vector<G4double> Edep = aHit->GetEdep();
600  vector<G4double> Dx = aHit->GetDx();
601  // Charge for each step
602  vector<int> charge = aHit->GetCharges();
603  vector<G4double> times = aHit->GetTime();
604  unsigned int nsteps = Edep.size();
605  double Etot = 0;
606  double X_hit_ave=0.;
607  double Y_hit_ave=0.;
608  double Z_hit_ave=0.;
609  double T_hit_ave=0.;
610  double dLeft =-10000.;
611 
612  for(unsigned int s=0; s<nsteps; s++) Etot = Etot + Edep[s];
613  if(Etot>0)
614  {
615  for(unsigned int s=0; s<nsteps; s++)
616  {
617  double Edep_B=Edep[s];
618  etot_g4=etot_g4+Edep_B;
619  // average hit position XYZ
620  X_hit_ave=X_hit_ave+Lpos[s].x();
621  Y_hit_ave=Y_hit_ave+Lpos[s].y();
622  Z_hit_ave=Z_hit_ave+Lpos[s].z();
623  // average hit time
624  T_hit_ave=T_hit_ave+times[s];
625 
626  //cout << "X " << Lpos[s].x() << " " << "Y " << Lpos[s].y() << " " << "Z " << Lpos[s].z() << " "<< "T " <<times[s] << " " << endl;
627 
628  }
629  X_hit_ave=X_hit_ave/nsteps;
630  Y_hit_ave=Y_hit_ave/nsteps;
631  Z_hit_ave=Z_hit_ave/nsteps;
632  T_hit_ave=T_hit_ave/nsteps;
633  dLeft =sz-Z_hit_ave;
634  timeL= dLeft/veff+T_hit_ave;
635 
636 
637  double *pe_sipm;// response for a mip (2..05 MeV energy released in 1cm thick)
638  pe_sipm=IVresponse(channel, X_hit_ave, Y_hit_ave, Z_hit_ave);
639 
640  ADC1=G4Poisson(pe_sipm[0]*etot_g4/2.05) ; // Scaling for more/less energy release)
641  ADC2=G4Poisson(pe_sipm[1]*etot_g4/2.05) ; // Scaling for more/less energy release)
642  ADC3=G4Poisson(pe_sipm[2]*etot_g4/2.05) ; // Scaling for more/less energy release)
643  ADC4=G4Poisson(pe_sipm[3]*etot_g4/2.05) ; // Scaling for more/less energy release)
644  double sigmaTL=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peL+1.));
645  sigmaTL=0.;
646  TDC1=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
647  TDC2=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
648  TDC3=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
649  TDC4=(timeL+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
650 
651 
652  // cout << "channel " << channel << endl;
653  // cout << "X " << X_hit_ave << " " << "Y " << Y_hit_ave << " " << "Z " << Z_hit_ave << " "<< "T " <<T_hit_ave << " " << endl;
654  // cout << "sipm1 " << pe_sipm[0] << " " << "sipm2 " << pe_sipm[1] << " " << "sipm3 " << pe_sipm[2] << " " << "sipm4 " << pe_sipm[3] << " " << endl;
655  // cout << "dLeft " << dLeft << " " << "timeL" << timeL << " " << endl;
656 
657  //cout << "energy right: " << ADC2 / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E left: " << ADC1 / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
658  //cout << "energy forw: " << ADCF / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E back: " << ADCB / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
659 
660  //cout << " Light collection: " << light_coll << endl;
661 
662  }
663  // closes (Etot > 0) loop
664 
665  }// end of IV
666 
667 
668  //starting paddles
669  if(veto_id==3)
670  {
671  double optical_coupling[3]= {0., 1.,0.37 };
672  for (int s=0; s<3; s++) optical_coupling[s] = optical_coupling[s]*0.34;
673 
674  light_yield=9200/MeV;
675  veff=13*cm/ns ;
676  sensor_surface=pow(1.27*cm,2)*pi; // 1" pmt R-> (2.5/2 TBC)
677  sensor_effective_area=0.9;
678 
679  sensor_qe=0.25;
680  sensor_gain=1.;
681 
682 
683  // Get the paddle length: in veto paddles are along z
684  length = aHit->GetDetector().dimensions[2];
685  double s1=aHit->GetDetector().dimensions[0];
686  double s2=aHit->GetDetector().dimensions[1];
687  paddle_surface = 2*s1*2*s2;
688  att_length=350*cm;
689  light_guide_att=1.;
690  //cout << " lenght: " << length << " optical-coupled surface: " << paddle_surface <<endl;
691 
692  light_coll=sensor_surface/paddle_surface;
693  if (sensor_surface>paddle_surface) light_coll=1.; // no more than the PMT size
694  light_coll=optical_coupling[channel]*light_coll*sensor_effective_area*light_guide_att; // Including the coupling efficiency and the pc effective area
695  // cout << " channel,veto: " << channel << " " << veto_id << " optical-couping: " << optical_coupling[channel] << " light coll: " << light_coll <<endl;
696  //cout << " light collo " << light_coll <<endl;
697 
698  // Get info about detector material to eveluate Birks effect
699  double birks_constant=aHit->GetDetector().GetLogical()->GetMaterial()->GetIonisation()->GetBirksConstant();
700  // cout << " Birks constant is: " << birks_constant << endl;
701  // cout << aHit->GetDetector().GetLogical()->GetMaterial()->GetName() << endl;
702 
703 
704  double time_min[4] = {0,0,0,0};
705 
706  vector<G4ThreeVector> Lpos = aHit->GetLPos();
707  vector<G4double> Edep = aHit->GetEdep();
708  vector<G4double> Dx = aHit->GetDx();
709  // Charge for each step
710  vector<int> charge = aHit->GetCharges();
711  vector<G4double> times = aHit->GetTime();
712 
713  unsigned int nsteps = Edep.size();
714  double Etot = 0;
715 
716  for(unsigned int s=0; s<nsteps; s++) Etot = Etot + Edep[s];
717  //for(unsigned int s=0; s<nsteps; s++) cout << "Energy = " << Edep[s]*1000.*1000. << endl;
718 
719  if(Etot>0)
720  {
721  for(unsigned int s=0; s<nsteps; s++)
722  {
723  double dLeft =-10000.;
724  double dRight =-10000.;
725 
726  // Distances from left, right for upper/lower (along z)
727  dLeft = length - Lpos[s].z();
728  dRight = length + Lpos[s].z();
729 
730 
731  //cout << "\n Distances: " << endl;
732  // cout << "\t dLeft, dRight " << dLeft << ", " << dRight << endl;
733 
734  // apply Birks effect
735  // double stepl = 0.;
736 
737  // if (s == 0){
738  // stepl = sqrt(pow((Lpos[s+1].x() - Lpos[s].x()),2) + pow((Lpos[s+1].y() - Lpos[s].y()),2) + pow((Lpos[s+1].z() - Lpos[s].z()),2));
739  // }
740  // else {
741  // stepl = sqrt(pow((Lpos[s].x() - Lpos[s-1].x()),2) + pow((Lpos[s].y() - Lpos[s-1].y()),2) + pow((Lpos[s].z() - Lpos[s-1].z()),2));
742  // }
743 
744  double Edep_B = BirksAttenuation(Edep[s],Dx[s],charge[s],birks_constant);
745  Edep_B=Edep[s];
746  etot_g4=etot_g4+Edep_B;
747  // cout << "\t Birks Effect: " << " Edep=" << Edep[s] << " StepL=" << stepl
748  // << " PID =" << pids[s] << " charge =" << charge[s] << " Edep_B=" << Edep_B << endl;
749 
750  //if (light_coll > 1) light_coll = 1.; // To make sure you don't miraculously get more energy than you started with
751 
752  etotL = etotL + Edep_B/2 * exp(-dLeft/att_length) * light_coll;
753  etotR = etotR + Edep_B/2 * exp(-dRight/att_length) * light_coll;
754 
755  // cout << "step: " << s << " etotL, etotR " << etotL << ", " << etotR << endl;
756 
757  timeL= timeL + (times[s] + dLeft/veff) / nsteps;
758  timeR= timeR + (times[s] + dRight/veff) / nsteps;
759 
760  if(etotL > 0.) {
761  if(s==0 || (time_min[0]>(times[s]+dLeft/veff))) time_min[0]=times[s]+dLeft/veff;
762  }
763  // cout << "min " << time_min[0] << "min " << times[s]+dLeft/veff << endl;
764  if(etotR > 0.) {
765  if(s==0 || (time_min[1]>(times[s]+dRight/veff))) time_min[1]=times[s]+ dRight/veff;
766  }
767  }
768  //cout << " etotR " << etotR << " ; etotL" << etotL <<endl;
769  peL=G4Poisson(etotL*light_yield*sensor_qe);
770  peR=G4Poisson(etotR*light_yield*sensor_qe);
771  //cout << " per " << peR << " ; pel" << peL <<endl;
772  //peL=(etotL*light_yield*sensor_qe);
773  //peR=(etotR*light_yield*sensor_qe);
774  //cout << " per " << peR << " ; pel" << peL <<endl;
775 
776  double sigmaTL=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peL+1.));
777  double sigmaTR=sqrt(pow(0.2*nanosecond,2.)+pow(1.*nanosecond,2.)/(peR+1.));
778 // sigmaTL=0;
779 // sigmaTR=0;
780  tL=(time_min[0]+G4RandGauss::shoot(0.,sigmaTL))*1000.;//time in ps
781  tR=(time_min[1]+G4RandGauss::shoot(0.,sigmaTR))*1000.;// time in ps
782  // cout << " tL " << tL << " ; timeL" << timeL <<endl;
783  // Digitization for ADC and QDC not used
784  //TDC1=(int) (tL * tdc_conv);
785  //TDC2=(int) (tR * tdc_conv);
786  //if(TDC1<0) TDC1=0;
787  //if(TDC2<0) TDC2=0;
788  //ADC1=(int) (peL*sensor_gain*adc_conv + adc_ped);
789  //ADC2=(int) (peR*sensor_gain*adc_conv + adc_ped);
790 
791  // cout << "ADC1: " << ADC1 << " " << peL << " " << sensor_gain << " " << adc_conv << endl;
792 
793 
794  //cout << "energy right: " << ADC2 / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E left: " << ADC1 / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
795  //cout << "energy forw: " << ADCF / (adc_conv*sensor_gain*sensor_qe*light_yield) << " E back: " << ADCB / (adc_conv*sensor_gain*sensor_qe*light_yield) << endl;
796 
797  //cout << " Light collection: " << light_coll << endl;
798 
799  }
800  // closes (Etot > 0) loop
801 
802 
803  if(verbosity>4)
804  {
805  cout << log_msg << " veto: " << veto_id << ", channel: " << channel ;
806  cout << log_msg << " Etot=" << Etot/MeV << endl;
807  cout << log_msg << " TDC1=" << TDC1 << " TDC2=" << TDC2 << " ADC1=" << ADC1 << " ADC2=" << ADC2 << endl;
808  }
809 
810 
811  // ch 1,3 dRight
812  // ch 2,4 dLeft
813  // ch 7,8,9,10,11,12 dRight
814  // ch 5,6 dRight
815  // ignore the other side
816 
817 
818  ADC1=peR;
819  TDC1=tR;
820  // cout << " ADC1: " << ADC1 << " ; TDC1: " << TDC1 << " ; ADC2: "<< etot_g4*1000. << endl;
821 
822  }// end of paddles
823 
824 
825  dgtz["hitn"] = hitn;
826  dgtz["sector"] = sector;
827  dgtz["veto"] = veto_id;
828  dgtz["channel"] = channel;
829  dgtz["adc1"] = ADC1;// output in pe
830  dgtz["adc2"] = ADC2;//deposited energy in keV
831  dgtz["adc3"] = ADC3;// ignore
832  dgtz["adc4"] = ADC4;// ignore
833  dgtz["tdc1"] = TDC1;// output in ps
834  dgtz["tdc2"] = TDC2;// ignore
835  dgtz["tdc3"] = TDC3;// ignore
836  dgtz["tdc4"] = TDC4;// ignore
837 
838  return dgtz;
839 }
840 
841 
842 vector<identifier> veto_HitProcess :: processID(vector<identifier> id, G4Step *step, detector Detector)
843 {
844  id[id.size()-1].id_sharing = 1;
845  return id;
846 }
847 
848 
849 double* veto_HitProcess::IVresponse(int channel, double xx, double yy,double zz)
850 {
851 
852  // Response of the different IV plastic paddles
853  // ch
854  //
855  //
856  static double response[4];
857 
858  for(unsigned int s=0; s<4; s++)response[s] = 0.;
859 
860  if (channel==1)//top
861  {
862  double x=xx/10.;
863  double y=(1058/2.-zz)/10.;
864 
865  double parm[4][8]={
866  {1.99627e+01, 1.64910e-01, -5.83528e-01, -7.34483e-03, -1.25062e-03, 4.43805e-03, 5.63766e-05, 1.40682e-05},
867  {1.86162e+01, 4.36475e-02, -6.78752e-02, -5.47887e-03, -1.60512e-04, -2.33958e-02, 5.55285e-05, -5.94424e-05},
868  {1.85966e+01, 1.96301e-01, 1.34868e-01, -7.66131e-04, -1.61720e-03, -1.91598e-02, -1.76198e-06, -4.72970e-05},
869  {9.73394e+00, 1.56111e-01, 3.27558e-01, 2.45041e-03, -1.31615e-03, 5.82688e-03, -1.48528e-05, 2.35177e-05}
870 
871  };
872 
873  for(unsigned int s=0; s<4; s++)
874  response[s] = parm[s][7]*x*x*y + parm[s][6]*x*y*y + parm[s][5]*x*x + parm[s][4]*y*y + parm[s][3]*x*y + parm[s][2]*x + parm[s][1]*y + parm[s][0];
875 
876  // cout << " x: " << x << " y: " << y << " zz: " << zz << endl ;
877  //cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
878 
879  }
880  else if (channel==2)//bottom
881  {// Assuming an overall size of 42.8 cm with 4 bars of
882  double x=-(xx-428/2)/10;
883  double y=(1058/2.-zz)/10.;
884 
885 
886  for(unsigned int s=0; s<4; s++) response[s] =0.;
887  if (x<10) response[0]= (-0.000303034)*y*y + (0.00658939)*y + 32.4847; //D1
888  if (x>10 && x <20 ) response[1]= (0.00301674)*y*y + (-0.446544)*y + 27.6374; //D4
889  if (x>20 && x <32.8 ) response[2]= (-0.000275694)*y*y + (0.00124251)*y + 18.8999; //D3
890  if (x>32.8 && x <42.8 ) response[3]= (-0.00139525)*y*y + (0.104993)*y + 18.1047; //D2
891 
892  // cout << " x: " << x << " y: " << y << " zz: " << zz << endl ;
893  // cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
894 
895  }
896  else if (channel==3)// Side Upstream
897  {// Assuming an overall size of 42.8 cm with 4 bars of
898  double x=-xx/10;
899  double y=(yy+346./2)/10.;
900 
901  double parm[4]={-0.04, -0.05, 1.4, 85.};
902 
903 
904  for(unsigned int s=0; s<4; s++) response[s] =0.;
905  response[0] = parm[0]*x*x + parm[1]*y*y + parm[2]*y + parm[3];
906 
907  // cout << " x: " << x << " y: " << y << " yy: " << yy << endl ;
908  // cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
909 
910  }
911  else if (channel==4)// Side Downstream
912  {// Assuming an overall size of 42.8 cm with 4 bars of
913  double x=xx/10;
914  double y=(yy+346./2)/10.;
915 
916  double parm[4]={-0.04, -0.05, 1.4, 75.};
917 
918 
919  for(unsigned int s=0; s<4; s++) response[s] =0.;
920  response[0] = parm[0]*x*x + parm[1]*y*y + parm[2]*y + parm[3];
921 
922  //cout << " x: " << x << " y: " << y << " yy: " << yy << endl ;
923  //cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
924 
925  }
926 
927  else if (channel==5)//Right
928  {
929  double x=-yy/10.;
930  double y=(1058/2.-zz)/10.;
931 
932  double parm[4][8]={
933  {2.34524e+01, 4.28317e-02, -5.91894e-01, -5.13309e-03, -2.47905e-04, -3.44887e-03, 4.25481e-05, -1.03817e-05},
934  {1.68313e+01, 5.36853e-02, -2.14037e-01, -4.80535e-03, -4.65364e-04, -1.66572e-02, 4.89028e-05, -3.33380e-05},
935  {2.50310e+01, -3.10007e-02, 3.57657e-01, -1.39833e-02, 2.99406e-04, -3.23669e-02, 1.27237e-04, -1.13100e-05},
936  {1.74834e+01, 1.83925e-01, 5.36737e-01, 7.09769e-04, -1.64490e-03, 7.48199e-03, 3.43011e-08, 2.11894e-05}
937 
938  };
939 
940  for(unsigned int s=0; s<4; s++)
941  response[s] = parm[s][7]*x*x*y + parm[s][6]*x*y*y + parm[s][5]*x*x + parm[s][4]*y*y + parm[s][3]*x*y + parm[s][2]*x + parm[s][1]*y + parm[s][0];
942 
943  //cout << " x: " << x << " y: " << y << " zz: " << zz << endl ;
944  //cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
945 
946  }
947  else if (channel==6)//Left
948  {
949  double x=-yy/10.;
950  double y=(1058/2.-zz)/10.;
951 
952  double parm[4][8]={
953  {8.12418e+00, 6.61315e-02, -2.99641e-01, -9.10408e-04, -6.79474e-04, 2.00648e-03, 1.24963e-05, -1.73809e-05},
954  {1.19501e+01, 4.76291e-02, -1.77047e-01, 9.27111e-05, -4.63061e-04, -1.40014e-02, 4.39766e-06, -2.93896e-05},
955  {1.68607e+01, -4.15476e-02, 2.54857e-01, -6.87363e-03, 3.26876e-04, -2.65178e-02, 5.62748e-05, -3.56067e-06},
956  {9.73394e+00, 1.56111e-01, 3.27558e-01, 2.45041e-03, -1.31615e-03, 5.82688e-03, -1.48528e-05, 2.35177e-05}
957 
958  };
959 
960  for(unsigned int s=0; s<4; s++)
961  response[s] = parm[s][7]*x*x*y + parm[s][6]*x*y*y + parm[s][5]*x*x + parm[s][4]*y*y + parm[s][3]*x*y + parm[s][2]*x + parm[s][1]*y + parm[s][0];
962 
963  }
964 
965  // cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
966  // cout << " res[0]: " << response[0] << " channel " << channel<<endl;
967  return response;
968 }
969 
970 double* veto_HitProcess::IVresponseProposal(int channel, double xx, double yy,double zz, double sx, double sy, double sz)
971 {
972  // Response of the different IV plastic paddles
973  // ch
974  //
975  //
976  static double response[4];
977 
978  for(unsigned int s=0; s<4; s++)response[s] = 0.;
979 
980  if (channel==1)//top
981  {
982  double x=xx/10.;
983  double y=(sz-zz)/10.;
984 
985  double parm[4][8]={
986  {1.99627e+01, 1.64910e-01, -5.83528e-01, -7.34483e-03, -1.25062e-03, 4.43805e-03, 5.63766e-05, 1.40682e-05},
987  {1.86162e+01, 4.36475e-02, -6.78752e-02, -5.47887e-03, -1.60512e-04, -2.33958e-02, 5.55285e-05, -5.94424e-05},
988  {1.85966e+01, 1.96301e-01, 1.34868e-01, -7.66131e-04, -1.61720e-03, -1.91598e-02, -1.76198e-06, -4.72970e-05},
989  {9.73394e+00, 1.56111e-01, 3.27558e-01, 2.45041e-03, -1.31615e-03, 5.82688e-03, -1.48528e-05, 2.35177e-05}
990 
991  };
992 
993  for(unsigned int s=0; s<4; s++)
994  response[s] = parm[s][7]*x*x*y + parm[s][6]*x*y*y + parm[s][5]*x*x + parm[s][4]*y*y + parm[s][3]*x*y + parm[s][2]*x + parm[s][1]*y + parm[s][0];
995 
996  // cout << " x: " << x << " y: " << y << " zz: " << zz << endl ;
997  //cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
998 
999  }
1000  else if (channel==2)//bottom
1001  {// Assuming an overall size of 42.8 cm with 4 bars of
1002  double x=-(xx-sx)/10;
1003  double y=(sz-zz)/10.;
1004 
1005 
1006  for(unsigned int s=0; s<4; s++) response[s] =0.;
1007  if (x<10) response[0]= (-0.000303034)*y*y + (0.00658939)*y + 32.4847; //D1
1008  if (x>10 && x <20 ) response[1]= (0.00301674)*y*y + (-0.446544)*y + 27.6374; //D4
1009  if (x>20 && x <32.8 ) response[2]= (-0.000275694)*y*y + (0.00124251)*y + 18.8999; //D3
1010  if (x>32.8 && x <42.8 ) response[3]= (-0.00139525)*y*y + (0.104993)*y + 18.1047; //D2
1011 
1012  // cout << " x: " << x << " y: " << y << " zz: " << zz << endl ;
1013  // cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
1014 
1015  }
1016  else if (channel==6)// Side Upstream
1017  {// Assuming an overall size of 42.8 cm with 4 bars of
1018  double x=-xx/10;
1019  double y=(yy+sy)/10.;
1020 
1021  double parm[4]={-0.04, -0.05, 1.4, 85.};
1022 
1023 
1024  for(unsigned int s=0; s<4; s++) response[s] =0.;
1025  response[0] = parm[0]*x*x + parm[1]*y*y + parm[2]*y + parm[3];
1026 
1027  // cout << " x: " << x << " y: " << y << " yy: " << yy << endl ;
1028  // cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
1029 
1030  }
1031  else if (channel==5)// Side Downstream
1032  {// Assuming an overall size of 42.8 cm with 4 bars of
1033  double x=xx/10;
1034  double y=(yy+sy)/10.;
1035 
1036  double parm[4]={-0.04, -0.05, 1.4, 75.};
1037 
1038 
1039  for(unsigned int s=0; s<4; s++) response[s] =0.;
1040  response[0] = parm[0]*x*x + parm[1]*y*y + parm[2]*y + parm[3];
1041 
1042  //cout << " x: " << x << " y: " << y << " yy: " << yy << endl ;
1043  //cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
1044 
1045  }
1046 
1047  else if (channel==3)//Right
1048  {
1049  double x=-yy/10.;
1050  double y=(sz-zz)/10.;
1051 
1052  double parm[4][8]={
1053  {2.34524e+01, 4.28317e-02, -5.91894e-01, -5.13309e-03, -2.47905e-04, -3.44887e-03, 4.25481e-05, -1.03817e-05},
1054  {1.68313e+01, 5.36853e-02, -2.14037e-01, -4.80535e-03, -4.65364e-04, -1.66572e-02, 4.89028e-05, -3.33380e-05},
1055  {2.50310e+01, -3.10007e-02, 3.57657e-01, -1.39833e-02, 2.99406e-04, -3.23669e-02, 1.27237e-04, -1.13100e-05},
1056  {1.74834e+01, 1.83925e-01, 5.36737e-01, 7.09769e-04, -1.64490e-03, 7.48199e-03, 3.43011e-08, 2.11894e-05}
1057 
1058  };
1059 
1060  for(unsigned int s=0; s<4; s++)
1061  response[s] = parm[s][7]*x*x*y + parm[s][6]*x*y*y + parm[s][5]*x*x + parm[s][4]*y*y + parm[s][3]*x*y + parm[s][2]*x + parm[s][1]*y + parm[s][0];
1062 
1063  //cout << " x: " << x << " y: " << y << " zz: " << zz << endl ;
1064  //cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
1065 
1066  }
1067  else if (channel==4)//Left
1068  {
1069  double x=-yy/10.;
1070  double y=(sz-zz)/10.;
1071 
1072  double parm[4][8]={
1073  {8.12418e+00, 6.61315e-02, -2.99641e-01, -9.10408e-04, -6.79474e-04, 2.00648e-03, 1.24963e-05, -1.73809e-05},
1074  {1.19501e+01, 4.76291e-02, -1.77047e-01, 9.27111e-05, -4.63061e-04, -1.40014e-02, 4.39766e-06, -2.93896e-05},
1075  {1.68607e+01, -4.15476e-02, 2.54857e-01, -6.87363e-03, 3.26876e-04, -2.65178e-02, 5.62748e-05, -3.56067e-06},
1076  {9.73394e+00, 1.56111e-01, 3.27558e-01, 2.45041e-03, -1.31615e-03, 5.82688e-03, -1.48528e-05, 2.35177e-05}
1077 
1078  };
1079 
1080  for(unsigned int s=0; s<4; s++)
1081  response[s] = parm[s][7]*x*x*y + parm[s][6]*x*y*y + parm[s][5]*x*x + parm[s][4]*y*y + parm[s][3]*x*y + parm[s][2]*x + parm[s][1]*y + parm[s][0];
1082 
1083  }
1084 
1085  // cout << " res[0]: " << response[0] << " res[1]: " << response[1]<< " res[2]: " << response[2]<< " res[3]: " << response[3] << endl ;
1086  // cout << " res[0]: " << response[0] << " channel " << channel<<endl;
1087  return response;
1088 }
1089 
1090 
1091 double veto_HitProcess::BirksAttenuation(double destep, double stepl, int charge, double birks)
1092 {
1093  //Example of Birk attenuation law in organic scintillators.
1094  //adapted from Geant3 PHYS337. See MIN 80 (1970) 239-244
1095  //
1096  // Taken from GEANT4 examples advanced/amsEcal and extended/electromagnetic/TestEm3
1097  //
1098  double response = destep;
1099  if (birks*destep*stepl*charge != 0.)
1100  {
1101  response = destep/(1. + birks*destep/stepl);
1102  }
1103  return response;
1104 }
1105 
1106 
1107 double veto_HitProcess::BirksAttenuation2(double destep,double stepl,int charge,double birks)
1108 {
1109  //Extension of Birk attenuation law proposed by Chou
1110  // see G.V. O'Rielly et al. Nucl. Instr and Meth A368(1996)745
1111  //
1112  //
1113  double C=9.59*1E-4*mm*mm/MeV/MeV;
1114  double response = destep;
1115  if (birks*destep*stepl*charge != 0.)
1116  {
1117  response = destep/(1. + birks*destep/stepl + C*pow(destep/stepl,2.));
1118  }
1119  return response;
1120 }
1121 
1122 
1123 map< string, vector <int> > veto_HitProcess :: multiDgt(MHit* aHit, int hitn)
1124 {
1125  map< string, vector <int> > MH;
1126 
1127  return MH;
1128 }
1129 
1130 
1131 
1132 
1133 
veto_hitprocess.h
veto_HitProcess::integrateDgt
map< string, double > integrateDgt(MHit *, int)
Definition: veto_hitprocess.cc:13
detector::GetLogical
G4LogicalVolume * GetLogical()
Returns Logical Volume pointer.
Definition: detector.h:113
veto_HitProcess::IVresponse
double * IVresponse(int, double, double, double)
Definition: veto_hitprocess.cc:849
verbosity
double verbosity
Definition: gdml_det_factory.cc:6
MHit::GetId
vector< identifier > GetId()
Definition: Hit.h:103
veto_HitProcess::BirksAttenuation2
double BirksAttenuation2(double, double, int, double)
Definition: veto_hitprocess.cc:1107
MHit::GetLPos
vector< G4ThreeVector > GetLPos()
Definition: Hit.h:76
veto_HitProcess::multiDgt
map< string, vector< int > > multiDgt(MHit *, int)
Definition: veto_hitprocess.cc:1123
MHit
Definition: Hit.h:22
MHit::GetTime
vector< double > GetTime()
Definition: Hit.h:89
veto_HitProcess::BirksAttenuation
double BirksAttenuation(double, double, int, double)
Definition: veto_hitprocess.cc:1091
detector::dimensions
vector< double > dimensions
vector of dimensions. Size, units depends on solid type
Definition: detector.h:77
veto_HitProcess::processID
vector< identifier > processID(vector< identifier >, G4Step *, detector)
Definition: veto_hitprocess.cc:842
MHit::GetEdep
vector< double > GetEdep()
Definition: Hit.h:83
veto_HitProcess::IVresponseProposal
double * IVresponseProposal(int, double, double, double, double, double, double)
Definition: veto_hitprocess.cc:970
MHit::GetDetector
detector GetDetector()
Definition: Hit.h:108
MHit::GetCharges
vector< int > GetCharges()
Definition: Hit.h:116
MHit::GetDx
vector< double > GetDx()
Definition: Hit.h:86
Generated on Wed Jul 6 2016 11:13:42 for GEMC by   doxygen 1.8.11