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

Last change on this file since 1282 was 807, checked in by garnier, 16 years ago

update
File size: 6.7 KB
Line 
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.