GEMC  1.8
Geant4 Monte-Carlo Framework
DC_hitprocess.cc
Go to the documentation of this file.
1 // %%%%%%%%%%%%%
2 // gemc headers
3 // %%%%%%%%%%%%%
4 #include "DC_hitprocess.h"
5 #include "CLHEP/Random/RandGaussT.h"
6 
7 // NOTE:
8 // Need to run gemc with -USE_PHYSICSL=gemc
9 // otherwise the step size inside the cell is too small
10 // to calculate DOCA correctly
11 
12 PH_output DC_HitProcess :: ProcessHit(MHit* aHit, gemc_opts Opt)
13 {
14  HCname = "DC Hit Process";
15 
16  double mini_stagger_shift_R2 = Opt.args["DC_MSTAG_R2"].arg*mm; // Mini Stagger shift for Region 2
17  double mini_stagger_shift_R3 = Opt.args["DC_MSTAG_R3"].arg*mm; // Mini Stagger shift for Region 3
18  double mini_stagger_shift;
19 
20  string hd_msg = Opt.args["LOG_MSG"].args + " DC Hit Process " ;
21  double HIT_VERBOSITY = Opt.args["HIT_VERBOSITY"].arg;
22 
23  PH_output out;
24  out.identity = aHit->GetId();
25  int sector = out.identity[0].id;
26  int SL = out.identity[1].id;
27  int Layer = out.identity[2].id;
28  int nwire = out.identity[3].id;
29 
30  mini_stagger_shift = 0;
31  if(SL > 2 && SL <= 4) // shift only for region 2
32  {
33  mini_stagger_shift = mini_stagger_shift_R2;
34  if(Layer == 2 || Layer == 4 || Layer == 6) mini_stagger_shift *= -1;
35  }
36  else if(SL > 4) // shift only for region 3
37  {
38  mini_stagger_shift = mini_stagger_shift_R3;
39  if(Layer == 2 || Layer == 4 || Layer == 6) mini_stagger_shift *= -1;
40  }
41 
42  // nwire position information
43  double NWIRES = 113;
44  double ylength = aHit->GetDetector().dimensions[3];
45  double deltay = 2.0*ylength/NWIRES;
46  double WIRE_Y = nwire*deltay + mini_stagger_shift;
47 
48  // drift velocities
49  double drift_velocity = 0;
50  if(SL == 1 || SL == 2) drift_velocity = 0.053;
51  if(SL == 3 || SL == 4) drift_velocity = 0.026;
52  if(SL == 5 || SL == 6) drift_velocity = 0.036;
53 
54  // %%%%%%%%%%%%%%%%%%%
55  // Raw hit information
56  // %%%%%%%%%%%%%%%%%%%
57  int nsteps = aHit->GetPos().size();
58  // Identifying the fastest - given by time + doca(s) / drift velocity
59  int trackId = -1;
60  double minTime = 10000;
61  double signal_t = 0;
62 
63  vector<int> stepTrackId = aHit->GetTIds();
64  vector<double> stepTime = aHit->GetTime();
65  vector<G4double> Edep = aHit->GetEdep();
66  vector<G4ThreeVector> pos = aHit->GetPos();
67  vector<G4ThreeVector> Lpos = aHit->GetLPos();
68 
69  for(int s=0; s<nsteps; s++)
70  {
71  G4ThreeVector DOCA(0, Lpos[s].y() + ylength - WIRE_Y, Lpos[s].z()); // local cylinder
72  signal_t = stepTime[s]/ns + DOCA.mag()/drift_velocity;
73 
74  if(signal_t < minTime && Edep[s] >= aHit->GetThreshold())
75  {
76  trackId = stepTrackId[s];
77  minTime = signal_t;
78  }
79  // debug mode:
80 /* cout << " step: " << s+1 << " w: " << nwire
81  << " L: " << Layer << " SL: " << SL
82  << " ms: " << mini_stagger_shift << " doca: " << DOCA.mag()
83  << " time: " << stepTime[s]/ns << " tdc: " << signal_t
84  << " step ID: " << stepTrackId[s] << " fast: " << trackId
85  << " Edep (MeV): " << Edep[s]/MeV << " Threshold: " << aHit->GetThreshold() << endl;*/
86  }
87 
88  // If no step pass the threshold, getting the fastest signal - given by time + doca(s) / drift velocity
89  if(trackId == -1)
90  {
91  for(int s=0; s<nsteps; s++)
92  {
93  G4ThreeVector DOCA(0, Lpos[s].y() + ylength - WIRE_Y, Lpos[s].z());
94  signal_t = stepTime[s]/ns + DOCA.mag()/drift_velocity;
95  if(signal_t < minTime)
96  {
97  trackId = stepTrackId[s];
98  minTime = signal_t;
99  }
100  }
101  }
102 
103  // Get Total Energy deposited of the fastest track
104  double Etot = 0;
105  int count_track_steps = 0;
106  for(int s=0; s<nsteps; s++)
107  if(stepTrackId[s] == trackId)
108  {
109  Etot = Etot + Edep[s];
110  count_track_steps++;
111  }
112 
113 
114  // average global, local positions of fastest track
115  double x, y, z;
116  double lx, ly, lz;
117  x = y = z = lx = ly = lz = 0;
118 
119  /* if(Etot>0)*/
120  for(int s=0; s<nsteps; s++)
121  {
122  if(stepTrackId[s] == trackId)
123  {
124  /* x = x + pos[s].x()*Edep[s]/Etot;
125  y = y + pos[s].y()*Edep[s]/Etot;
126  z = z + pos[s].z()*Edep[s]/Etot;
127  lx = lx + Lpos[s].x()*Edep[s]/Etot;
128  ly = ly + Lpos[s].y()*Edep[s]/Etot;
129  lz = lz + Lpos[s].z()*Edep[s]/Etot;*/
130  x = x + pos[s].x()/count_track_steps;
131  y = y + pos[s].y()/count_track_steps;
132  z = z + pos[s].z()/count_track_steps;
133  lx = lx + Lpos[s].x()/count_track_steps;
134  ly = ly + Lpos[s].y()/count_track_steps;
135  lz = lz + Lpos[s].z()/count_track_steps;
136  }
137  }
138  // else
139  // {
140  // x = pos[0].x();
141  // y = pos[0].y();
142  // z = pos[0].z();
143  // lx = Lpos[0].x();
144  // ly = Lpos[0].y();
145  // lz = Lpos[0].z();
146  // }
147 
148  // average time of fastest track
149  double time = 0;
150  vector<G4double> times = aHit->GetTime();
151  for(int s=0; s<nsteps; s++) if(stepTrackId[s] == trackId) time = time + times[s]/count_track_steps;
152 
153  // Energy of the track
154  double Ene = aHit->GetE();
155 
156  out.raws.push_back(Etot);
157  out.raws.push_back(x);
158  out.raws.push_back(y);
159  out.raws.push_back(z);
160  out.raws.push_back(lx);
161  out.raws.push_back(ly);
162  out.raws.push_back(lz);
163  out.raws.push_back(time);
164  out.raws.push_back((double) aHit->GetPID());
165  out.raws.push_back(aHit->GetVert().getX());
166  out.raws.push_back(aHit->GetVert().getY());
167  out.raws.push_back(aHit->GetVert().getZ());
168  out.raws.push_back(Ene);
169  out.raws.push_back((double) aHit->GetmPID());
170  out.raws.push_back(aHit->GetmVert().getX());
171  out.raws.push_back(aHit->GetmVert().getY());
172  out.raws.push_back(aHit->GetmVert().getZ());
173 
174 
175 
176  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
177  // Distance 1: smear by 300 um
178  // Drift Time 1: using 50 um/ns
179  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
180 
181 
182  // Finding DOCA
183  double doca = 10000;
184  double LR = 0;
185 
186  for(int s=0; s<nsteps; s++)
187  {
188  G4ThreeVector DOCA(0, Lpos[s].y() + ylength - WIRE_Y, Lpos[s].z());
189  if(DOCA.mag() <= doca && stepTrackId[s] == trackId )
190  {
191  doca = DOCA.mag();
192  if(DOCA.y() >=0 ) LR = 1;
193  else LR = -1;
194 
195  }
196  // // debug mode:
197  // cout << " wire: " << nwire << " layer: " << iden[2].id << " SL: " << SL
198  // << " doca: " << doca << " Distance: " << DOCA.mag() << " start Time: " << stepTime[s]/ns
199  // << " signal: " << signal_t << " step: " << s+1 << " track ID: " << stepTrackId[s]
200  // << " fast track: " << trackId << " Edep (MeV): " << Edep[s]/MeV << endl;
201 
202  }
203 
204 
205 
206 
207  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
208  // Digitization
209  // DC ID:
210  // sector, SL, Layer, nwire
211  // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
212 
213  double sdoca = -1;
214  sdoca = fabs(CLHEP::RandGauss::shoot(doca, 0.3));
215  double time1 = doca/drift_velocity;
216  double stime1 = sdoca/drift_velocity;
217 
218  out.raws.push_back(LR);
219  out.raws.push_back(doca);
220  out.raws.push_back(sdoca);
221  out.raws.push_back(time1);
222  out.raws.push_back(stime1);
223 
224  out.dgtz.push_back(sector);
225  out.dgtz.push_back(SL);
226  out.dgtz.push_back(Layer);
227  out.dgtz.push_back(nwire);
228 
229  if(HIT_VERBOSITY>4)
230  cout << hd_msg
231  << " PID: " << aHit->GetPID()
232  << " Sector: " << sector
233  << " Superlayer: " << SL
234  << " Layer: " << Layer
235  << " Cell Size: " << deltay
236  << " nwire: " << nwire
237  << " L/R: " << LR
238  << " doca: " << doca
239  << " sdoca: " << sdoca
240  << " time1: " << time1
241  << " stime1: " << stime1 << endl;
242 
243  return out;
244 }
245 
246 vector<identifier> DC_HitProcess :: ProcessID(vector<identifier> id, G4Step* aStep, detector Detector, gemc_opts Opt)
247 {
248  vector<identifier> yid = id;
249  double mini_stagger_shift_R2 = Opt.args["DC_MSTAG_R2"].arg*mm; // Mini Stagger shift for Region 2
250  double mini_stagger_shift_R3 = Opt.args["DC_MSTAG_R3"].arg*mm; // Mini Stagger shift for Region 3
251  double mini_stagger_shift;
252 
253  double NWIRES = 113;
254  int SL = yid[1].id;
255  int Layer = yid[2].id;
256 
257  mini_stagger_shift = 0;
258  if(SL > 2 && SL <= 4) // shift only for region 2
259  {
260  mini_stagger_shift = mini_stagger_shift_R2;
261  if(Layer == 2 || Layer == 4 || Layer == 6) mini_stagger_shift *= -1;
262  }
263  else if(SL > 4) // shift only for region 3
264  {
265  mini_stagger_shift = mini_stagger_shift_R3;
266  if(Layer == 2 || Layer == 4 || Layer == 6) mini_stagger_shift *= -1;
267  }
268 
269  G4StepPoint *prestep = aStep->GetPreStepPoint();
270  G4StepPoint *poststep = aStep->GetPostStepPoint();
271  G4VTouchable* TH = (G4VTouchable*) aStep->GetPreStepPoint()->GetTouchable();
272 
273  string name = TH->GetVolume(0)->GetName();
274  G4ThreeVector xyz = poststep->GetPosition();
275  G4ThreeVector Lxyz = prestep->GetTouchableHandle()->GetHistory()
276  ->GetTopTransform().TransformPoint(xyz);
277 
278 
279  double ylength = Detector.dimensions[3];
280  double deltay = 2.0*ylength/NWIRES;
281  double loc_y = Lxyz.y() + ylength - mini_stagger_shift;
282 
283  int nwire = (int) floor(loc_y/deltay);
284  if(nwire <= 0 ) nwire = 1;
285  if(nwire == 113) nwire = 112;
286 
287  yid[3].id = nwire;
288 
289  if(fabs( (nwire+1)*deltay - loc_y ) < fabs( nwire*deltay - loc_y ) && nwire != 112 )
290  yid[3].id = nwire + 1;
291 
292 /* if(nwire > NWIRES) cout << " SuperLayer: " << SL << " layer: " << Layer
293  << " wire: " << nwire << " length: " << ylength
294  << " ly: " << Lxyz.y() << " deltay: " << deltay
295  << " loc_y: " << loc_y << " nwire*deltay: " << fabs( nwire*deltay - loc_y ) << endl;*/
296 
297  yid[3].id_sharing = 1;
298 
299  return yid;
300 }
301 
302 
303 
304 
305 
PH_output
Definition: MPHBaseClass.h:23
MHit::GetThreshold
double GetThreshold()
Definition: MHit.h:104
DC_HitProcess::ProcessHit
PH_output ProcessHit(MHit *, gemc_opts)
Method to process the hit.
Definition: DC_hitprocess.cc:12
MHit::GetVert
G4ThreeVector GetVert()
Definition: MHit.h:72
PH_output::raws
vector< double > raws
Raw information.
Definition: MPHBaseClass.h:26
MHit::GetPID
int GetPID()
Definition: MHit.h:107
MHit::GetmPID
int GetmPID()
Definition: MHit.h:122
MHit::GetId
vector< identifier > GetId()
Definition: MHit.h:96
MPHBaseClass::HCname
string HCname
Hit Collection name.
Definition: MPHBaseClass.h:41
MHit::GetmVert
G4ThreeVector GetmVert()
Definition: MHit.h:127
PH_output::identity
vector< identifier > identity
Identifier.
Definition: MPHBaseClass.h:28
MHit::GetLPos
vector< G4ThreeVector > GetLPos()
Definition: MHit.h:69
MHit::GetE
double GetE()
Definition: MHit.h:89
MHit::GetPos
vector< G4ThreeVector > GetPos()
Definition: MHit.h:65
MHit
Definition: MHit.h:29
MHit::GetTime
vector< double > GetTime()
Definition: MHit.h:82
gemc_opts::args
map< string, opts > args
Options map.
Definition: usage.h:68
MHit::GetTIds
vector< int > GetTIds()
Definition: MHit.h:94
detector::dimensions
vector< double > dimensions
vector of dimensions. Size, units depends on solid type
Definition: detector.h:63
gemc_opts
Definition: usage.h:59
MHit::GetEdep
vector< double > GetEdep()
Definition: MHit.h:76
MHit::GetDetector
detector GetDetector()
Definition: MHit.h:101
PH_output::dgtz
vector< int > dgtz
Digitized information.
Definition: MPHBaseClass.h:27
DC_HitProcess::ProcessID
vector< identifier > ProcessID(vector< identifier >, G4Step *, detector, gemc_opts)
Method to calculate new identifier.
Definition: DC_hitprocess.cc:246
DC_hitprocess.h
Generated on Wed Jul 6 2016 11:12:55 for GEMC by   doxygen 1.8.11