147{
148 G4Track* gTrack = aStep->GetTrack() ;
149
150 G4double globalT=gTrack->GetGlobalTime();
151 if(isnan(globalT)){
152 G4cout<<"MdcSD:error, globalT is nan "<<G4endl;
153 return false;
154 }
155 if(globalT > 2000)return false;
156
157
158 G4int charge = gTrack->GetDefinition()->GetPDGCharge();
159 if (charge == 0) return false;
160
161
162 G4double stepLength=aStep->GetStepLength();
163 if(stepLength==0){
164
165 return false;
166 }
167
168 G4double edep = aStep->GetTotalEnergyDeposit() / stepLength;
169 if(edep==0.) return false;
170
171
172 G4StepPoint* prePoint = aStep->GetPreStepPoint() ;
173 G4StepPoint* postPoint = aStep->GetPostStepPoint() ;
174
175
176 G4ThreeVector pointIn = prePoint->GetPosition();
177 G4ThreeVector pointOut = postPoint->GetPosition();
178
179
180 const G4VTouchable *touchable = prePoint->GetTouchable();
181 G4VPhysicalVolume *volume = touchable->GetVolume(0);
182
183 G4double driftD = 0;
184 G4double driftT = 0;
185 G4double edepTemp = 0;
186 G4double lengthTemp = 0;
187 G4int cellId=0;
188 G4int layerId = touchable->GetVolume(1)->GetCopyNo();
189 if(volume->IsReplicated()){
190 cellId = touchable->GetReplicaNumber();
191 }else{
192 cellId=touchable->GetVolume(0)->GetCopyNo();
193 }
194 if(layerId==18&&(cellId==27||cellId==42))return false;
195
197
198
199 if(layerId >= 36) layerId = 36 + (layerId-36)/2;
200 }
201
202 G4double halfLayerLength=mdcGeomSvc->
Layer(layerId)->
Length()/2.;
203 if(((fabs(pointIn.z())-halfLayerLength)>-7.)
204 &&((fabs(pointOut.z())-halfLayerLength)>-7.))return false;
205
206 G4int trackID= gTrack->GetTrackID();
207 G4int truthID, g4TrackID;
209
210 G4double theta=gTrack->GetMomentumDirection().theta();
211
212 G4ThreeVector hitPosition=0;
213 G4double transferT=0;
214 driftD =
Distance(layerId,cellId,pointIn,pointOut,hitPosition,transferT);
215
216 G4double posPhi, wirePhi;
217 posPhi=hitPosition.phi();
218 if(posPhi<0)posPhi += 2*
pi;
219 wirePhi=mdcGeoPointer->
SignalWire(layerId,cellId).
Phi(hitPosition.z());
220
221 G4int posFlag;
222 if(posPhi<=wirePhi){
223 posFlag = 0;
224 }else{
225 posFlag = 1;
226 }
227
228 if(posPhi < 1. && wirePhi > 5.)posFlag = 1;
229 if(posPhi > 5. && wirePhi < 1.)posFlag = 0;
230
231 G4ThreeVector hitLine=pointOut-pointIn;
232 G4double enterAngle=hitLine.phi()-wirePhi;
233 while(enterAngle<-
pi/2.)enterAngle+=
pi;
234 while(enterAngle>
pi/2.)enterAngle-=
pi;
235
236
237 G4int pdg_code = gTrack->GetDefinition()->GetPDGEncoding();
238 G4String process_name="Generator";
239 if(NULL != gTrack->GetCreatorProcess()) process_name = gTrack->GetCreatorProcess()->GetProcessName();
240 G4double trkLen = gTrack->GetTrackLength();
241 G4int isSecondary = (trackID==g4TrackID)? 0:1;
242 G4ThreeVector p3_pre = prePoint->GetMomentum();
243 G4ThreeVector p3_post = postPoint->GetMomentum();
244 G4ThreeVector p3_hit = 0.5*(p3_pre+p3_post);
245
246
248 G4double betaGamma=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma();
249 if(betaGamma<0.01)return false;
250
251
252 G4double eCount=dedxSample(betaGamma, charge, theta);
253 edep=eCount;
254 }
255
258
262 newHit->
SetPos(hitPosition);
273
274
276
277 driftT=mdcCalPointer->
D2T(driftD);
278 globalT+=transferT;
279 driftT+=globalT;
280
283
284 if (hitPointer[layerId][cellId] == -1) {
285 hitsCollection->insert(newHit);
286 G4int NbHits = hitsCollection->entries();
287 hitPointer[layerId][cellId]=NbHits-1;
288 }
289 else
290 {
291 G4int pointer=hitPointer[layerId][cellId];
292 if (g4TrackID == trackID) {
293 G4double preDriftT=(*hitsCollection)[pointer]->GetDriftT();
294 }
295
296 G4double preDriftT = (*hitsCollection)[pointer]->GetDriftT();
297 if (driftT < preDriftT) {
298
299
300
301
302
303
304
305
306
307
308
309 *((*hitsCollection)[pointer])=*newHit;
310 }
311
312 delete newHit;
313 }
314
315
316 if(truthCollection){
317 if(g4TrackID==trackID){
318 G4int pointer=truthPointer[layerId][cellId];
319 if(pointer==-1){
325 truthHit->
SetPos (hitPosition);
332
333 truthCollection->insert(truthHit);
334 G4int NbHits = truthCollection->entries();
335 truthPointer[layerId][cellId]=NbHits-1;
336 }
337 else {
338 if(truthID == (*truthCollection)[pointer]->GetTrackID()){
339 G4double preDriftT=(*truthCollection)[pointer]->GetDriftT();
340 if(driftT<preDriftT){
341 (*truthCollection)[pointer]->SetDriftD(driftD);
342 (*truthCollection)[pointer]->SetDriftT(driftT);
343 (*truthCollection)[pointer]->SetPos(hitPosition);
344 (*truthCollection)[pointer]->SetGlobalT(globalT);
345 (*truthCollection)[pointer]->SetTheta(theta);
346 (*truthCollection)[pointer]->SetPosFlag(posFlag);
347 (*truthCollection)[pointer]->SetEnterAngle(enterAngle);
348 }
349 edepTemp=(*truthCollection)[pointer]->GetEdep();
350 (*truthCollection)[pointer]->SetEdep(edepTemp+edep);
351 } else {
357 truthHit->
SetPos(hitPosition);
364
365 truthCollection->insert(truthHit);
366 G4int NbHits = truthCollection->entries();
367 truthPointer[layerId][cellId]=NbHits-1;
368 }
369 }
370 }
371 }
372
373
374
375
376 return true;
377}
double D2T(double driftDNew)
void SetHitPointer(BesMdcHit *hit)
BesMdcWire SignalWire(int, int)
void SetEdep(G4double de)
void SetPDGCode(G4int code)
void SetDriftT(G4double time)
void SetEnterAngle(G4double angle)
void SetIsSecondary(G4int isSec)
void SetCellNo(G4int cell)
void SetCreatorProcess(G4String proName)
void SetPos(G4ThreeVector xyz)
void SetTrackID(G4int track)
void SetLayerNo(G4int layer)
void SetTheta(G4double angle)
void SetDriftD(G4double distance)
void SetMomentum(G4ThreeVector xyz)
void SetGlobalT(G4double time)
void SetPosFlag(G4int flag)
void SetFlightLength(G4double len)
G4double Distance(G4int, G4int, G4ThreeVector, G4ThreeVector, G4ThreeVector &, G4double &)
void GetCurrentTrackIndex(G4int &trackIndex, G4int &g4TrackId) const