source: trunk/examples/novice/gemc/src/DC_hitprocess.cc@ 1226

Last change on this file since 1226 was 807, checked in by garnier, 17 years ago

update
File size: 6.7 KB
RevLine 
[807]1// %%%%%%%%%%%%%
2// gemc headers
3// %%%%%%%%%%%%%
4#include "DC_hitprocess.h"
5#include "CLHEP/Random/RandGaussT.h"
6
7PH_output DC_HitProcess :: ProcessHit(MHit* aHit, gemc_opts Opt)
8{
9 string hd_msg = Opt.args["LOG_MSG"].args + " DC Hit Process " ;
10 double HIT_VERBOSITY = Opt.args["HIT_VERBOSITY"].arg;
11
12 PH_output out;
13 out.identity = aHit->GetId();
14
15 HCname = "DC Hit Process";
16 // %%%%%%%%%%%%%%%%%%%
17 // Raw hit information
18 // %%%%%%%%%%%%%%%%%%%
19 int nsteps = aHit->GetPos().size();
20
21
22
23 // Identifying the fastest track in the hit
24 int trackId = -1;
25 double minTime = 10000;
26 vector<int> stepTrackId = aHit->GetTrackId();
27 vector<double> stepTime = aHit->GetTime();
28 vector<G4double> Edep = aHit->GetEdep();
29
30 for(int s=0; s<nsteps; s++)
31 {
32 if(stepTrackId[s] != trackId && stepTime[s]/ns < minTime && Edep[s] >= aHit->GetThreshold())
33 trackId = stepTrackId[s];
34 }
35
36 // If no step pass the threshold, getting the fastest track
37 if(trackId == -1)
38 {
39 for(int s=0; s<nsteps; s++)
40 {
41 if(stepTrackId[s] != trackId && stepTime[s]/ns < minTime)
42 trackId = stepTrackId[s];
43 }
44 }
45
46 // Get Total Energy deposited of the fastest track
47 double Etot = 0;
48 int count_track_steps = 0;
49 for(int s=0; s<nsteps; s++)
50 if(stepTrackId[s] == trackId)
51 {
52 Etot = Etot + Edep[s];
53 count_track_steps++;
54 }
55
56
57 // average global, local positions of fastest track
58 double x, y, z;
59 double lx, ly, lz;
60 x = y = z = lx = ly = lz = 0;
61 vector<G4ThreeVector> pos = aHit->GetPos();
62 vector<G4ThreeVector> Lpos = aHit->GetLPos();
63
64/* if(Etot>0)*/
65 for(int s=0; s<nsteps; s++)
66 {
67 if(stepTrackId[s] == trackId)
68 {
69/* x = x + pos[s].x()*Edep[s]/Etot;
70 y = y + pos[s].y()*Edep[s]/Etot;
71 z = z + pos[s].z()*Edep[s]/Etot;
72 lx = lx + Lpos[s].x()*Edep[s]/Etot;
73 ly = ly + Lpos[s].y()*Edep[s]/Etot;
74 lz = lz + Lpos[s].z()*Edep[s]/Etot;*/
75 x = x + pos[s].x()/count_track_steps;
76 y = y + pos[s].y()/count_track_steps;
77 z = z + pos[s].z()/count_track_steps;
78 lx = lx + Lpos[s].x()/count_track_steps;
79 ly = ly + Lpos[s].y()/count_track_steps;
80 lz = lz + Lpos[s].z()/count_track_steps;
81 }
82 }
83// else
84// {
85// x = pos[0].x();
86// y = pos[0].y();
87// z = pos[0].z();
88// lx = Lpos[0].x();
89// ly = Lpos[0].y();
90// lz = Lpos[0].z();
91// }
92
93 // average time of fastest track
94 double time = 0;
95 vector<G4double> times = aHit->GetTime();
96 for(int s=0; s<nsteps; s++) if(stepTrackId[s] == trackId) time = time + times[s]/nsteps;
97
98 // Energy of the track
99 double Ene = aHit->GetE();
100
101 out.raws.push_back(Etot);
102 out.raws.push_back(x);
103 out.raws.push_back(y);
104 out.raws.push_back(z);
105 out.raws.push_back(lx);
106 out.raws.push_back(ly);
107 out.raws.push_back(lz);
108 out.raws.push_back(time);
109 out.raws.push_back((double) aHit->GetPID());
110 out.raws.push_back(aHit->GetVert().getX());
111 out.raws.push_back(aHit->GetVert().getY());
112 out.raws.push_back(aHit->GetVert().getZ());
113 out.raws.push_back(Ene);
114 out.raws.push_back((double) aHit->GetmPID());
115 out.raws.push_back(aHit->GetmVert().getX());
116 out.raws.push_back(aHit->GetmVert().getY());
117 out.raws.push_back(aHit->GetmVert().getZ());
118
119
120
121 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
122 // Distance 1: smear by 300 um
123 // Drift Time 1: using 50 um/ns
124 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
125
126 // Finding Wire position
127 double NWIRES = 112;
128 double ylength = aHit->GetDetector().dimensions[3];
129 double deltay = 2.0*ylength/NWIRES;
130 vector<identifier> iden = aHit->GetId();
131 int nwire = iden[3].id;
132 double WIRE_Y = nwire*deltay;
133
134 // Finding DOCA
135 double doca = 10000;
136 double LR = 0;
137
138 for(int s=0; s<nsteps; s++)
139 {
140 G4ThreeVector DOCA(0, Lpos[s].y() + ylength - WIRE_Y, Lpos[s].z());
141 if(DOCA.mag() <= doca && stepTrackId[s] == trackId )
142 {
143 doca = DOCA.mag();
144 if(DOCA.y() >=0 ) LR = 1;
145 else LR = -1;
146
147 }
148 }
149// if(aHit->GetPID() == 2112)
150// cout << " wire=" << nwire << " layer=" << iden[2].id
151// << " doca=" << doca << " nsteps=" << count_track_steps << endl;
152
153
154
155
156 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
157 // Digitization
158 // DC ID:
159 // sector, SL, Layer, Wire
160 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%
161 int sector = out.identity[0].id;
162 int SL = out.identity[1].id;
163 int Layer = out.identity[2].id;
164 int Wire = out.identity[3].id;
165
166 double drift_velocity = 0;
167 if(SL == 1 || SL == 2) drift_velocity = 0.053;
168 if(SL == 3 || SL == 4) drift_velocity = 0.026;
169 if(SL == 5 || SL == 6) drift_velocity = 0.036;
170
171 double sdoca = -1;
172 while(sdoca < 0) sdoca = CLHEP::RandGauss::shoot(doca, 0.3);
173 double time1 = doca/drift_velocity;
174 double stime1 = sdoca/drift_velocity;
175
176 out.raws.push_back(LR);
177 out.raws.push_back(doca);
178 out.raws.push_back(sdoca);
179 out.raws.push_back(time1);
180 out.raws.push_back(stime1);
181
182
183 out.dgtz.push_back(sector);
184 out.dgtz.push_back(SL);
185 out.dgtz.push_back(Layer);
186 out.dgtz.push_back(Wire);
187
188 if(HIT_VERBOSITY>4)
189 cout << hd_msg << " Sector: " << sector
190 << " Superlayer: " << SL
191 << " Layer: " << Layer
192 << " Wire: " << Wire
193 << " L/R: " << LR
194 << " doca: " << doca
195 << " sdoca: " << sdoca
196 << " time1" << time1
197 << " stime1=" << stime1 << endl;
198
199 return out;
200}
201
202vector<identifier> DC_HitProcess :: ProcessID(vector<identifier> id, G4Step* aStep, detector Detector)
203{
204 double NWIRES = 112;
205
206 G4StepPoint *prestep = aStep->GetPreStepPoint();
207 G4StepPoint *poststep = aStep->GetPostStepPoint();
208 G4VTouchable* TH = (G4VTouchable*) aStep->GetPreStepPoint()->GetTouchable();
209
210 string name = TH->GetVolume(0)->GetName();
211 G4ThreeVector xyz = poststep->GetPosition();
212 G4ThreeVector Lxyz = prestep->GetTouchableHandle()->GetHistory()
213 ->GetTopTransform().TransformPoint(xyz);
214
215 double ylength = Detector.dimensions[3];
216 double deltay = 2.0*ylength/NWIRES;
217 double loc_y = Lxyz.y() + ylength;
218
219 int nwire = (int) floor(loc_y/deltay);
220 vector<identifier> yid = id;
221
222 if(nwire > NWIRES) cout << " SuperLayer: " << yid[1].id << " layer: " << yid[2].id
223 << " wire: " << nwire << " length: " << ylength
224 << " ly: " << Lxyz.y() << " deltay: " << deltay
225 << " loc_y: " << loc_y << " nwire*deltay: " << fabs( nwire*deltay - loc_y ) << endl;
226
227
228 yid[3].id = nwire;
229
230 if(fabs( (nwire+1)*deltay - loc_y ) < fabs( nwire*deltay - loc_y ) )
231 yid[3].id = nwire + 1;
232
233 return yid;
234}
235
236
237
238
239
Note: See TracBrowser for help on using the repository browser.