QwGeant4
QweakSimSteppingAction.cc
Go to the documentation of this file.
1 //=============================================================================
2 //
3 // ---------------------------
4 // | Doxygen File Information |
5 // ---------------------------
6 //
7 /**
8 
9  \file QweakSimSteppingAction.cc
10 
11  $Revision: 1.4 $
12  $Date: 2006/05/05 21:46:05 $
13 
14  \author Klaus Hans Grimm
15 
16 */
17 
18 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
19 
21 
22 // geant4 includes
23 #include "G4TrackVector.hh"
24 #include "G4ParticleDefinition.hh"
25 #include "G4Gamma.hh"
26 #include "G4Electron.hh"
27 #include "G4Positron.hh"
28 #include "G4OpticalPhoton.hh"
29 #include "G4PionMinus.hh"
30 #include "G4PionPlus.hh"
31 #include "G4OpBoundaryProcess.hh"
32 #include "G4SDManager.hh"
33 #include "G4UnitsTable.hh"
34 
35 // user includes
39 #include "QweakSimEPEvent.hh"
41 #include "QweakSimPMTOnly_PMTSD.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 
46 {
47 
48  G4cout << "###### Calling QweakSimSteppingAction::QweakSimSteppingAction() " << G4endl;
49 
50  myEventCounter = 0;
51  fSecondary = NULL;
52  myUserInfo = myUInfo;
53 
54  myEvent = myEPEvent;
55 
56  evtGenStatus = 0;
58  //RandomPositionZ = myEvent->GetVertexZ();
59 
60  //std::ofstream EventDataFile("Event.dat", std::ios::out);
61 
62  G4cout << "###### Leaving QweakSimSteppingAction::QweakSimSteppingAction() " << G4endl;
63 
64 }
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 void QweakSimSteppingAction::UserSteppingAction(const G4Step* theStep){
68 
69  fSecondary = fpSteppingManager->GetfSecondary();
70 
71  // ** all the track info is postStep **
72  G4Track* theTrack = theStep->GetTrack();
73  G4ParticleDefinition* particleType = theTrack->GetDefinition();
74  G4StepPoint* thePrePoint = theStep->GetPreStepPoint();
75  G4VPhysicalVolume* thePrePV = thePrePoint->GetPhysicalVolume();
76  G4StepPoint* thePostPoint = theStep->GetPostStepPoint();
77  G4VPhysicalVolume* thePostPV = thePostPoint->GetPhysicalVolume();
78  G4TouchableHistory* theTouchable = (G4TouchableHistory*)(thePrePoint->GetTouchable());
79  // G4int ReplicaNo = 0;
80  G4String particleName = theTrack->GetDefinition()->GetParticleName();
81  G4ProcessManager* pm = particleType->GetProcessManager();
82 
83 //get material
84  G4Material* theMaterial = theTrack->GetMaterial();
85 
86  // G4int nprocesses = pm->GetProcessListLength();
87  // G4ProcessVector* pv = pm->GetProcessList();
88  // G4VSteppingVerbose* theVerbStep = G4VSteppingVerbose::GetInstance();
89  // G4double charge = particleType->GetPDGCharge();
90 
91  // G4int nSecAtRest = GetNumOfAtRestSecondaries();
92  // G4int nSecAlong = GetNumOfAlongStepSecondaries();
93  // G4int nSecPost = GetNumOfPostStepSecondaries();
94  // G4int nSecTotal = GetTotalNumOfSecondaries();
95 
96  //RandomPositionZ = myEvent->GetVertexZ();
98 
99  //jpan@nuclear.uwinnipeg.ca Thu Apr 16 01:33:14 CDT 2009
100  // check if it is primary
101 
102  //to get the kryptonie to work
103  if(theMaterial->GetName()=="Kryptonite")
104  {
105  theTrack->SetTrackStatus(fKillTrackAndSecondaries);
106  }
107 
108  G4int parentID = theTrack->GetParentID();
109  if( (particleType==G4Electron::ElectronDefinition()||particleType==G4PionMinus::PionMinusDefinition()||particleType==G4PionPlus::PionPlusDefinition()) && parentID==0 ){
110 
111  //jpan: to account for the energy loss before the event generation,
112  // force to change primary momentum direction here
113 
114  //scattering only occur inside reaction region of the target, only occur once
115  G4ThreeVector thePosition = theTrack->GetPosition();
116  G4double theX = thePosition.getX();
117  G4double theY = thePosition.getY();
118  G4double theZ = thePosition.getZ();
119  // G4cout << "Track pos:: " << theX/mm << "\t" << theY/mm << "\t" << theZ/mm << G4endl;
120  // G4cout << "RandomPositionZ :: " << RandomPositionZ << G4endl;
121 
122  G4String procName = thePostPoint->GetProcessDefinedStep()->GetProcessName();
123  // dEE in MeV -> all Elosses are in MeV
124  G4double dEE = thePrePoint->GetKineticEnergy()/MeV-thePostPoint->GetKineticEnergy()/MeV;
125 
126  if(myUserInfo->GetPrimaryEventNumber() %2!=0){
127 
128  // if( theZ > targetCenterPositionZ+35*cm*0.5+5*(2.54*cm*0.001) || sqrt(theX*theX+theY*theY)>2.54*cm)
129  // {
130  // evtGenStatus = 0;
131  // }
132 
133  if( myUserInfo->EvtGenStatus == 0){ // true for primary event # = 0,2,4,...
134  // see QweakSimPrimaryGeneratorAction.cc
135 
136  // // various energy losses at the target
137  // // IMP:: all the Elosses upto the vertex, theZ, is stored here
138  if ( thePrePV && thePostPV) { // Added check to prevent a seg fault
139  if(((thePrePV->GetName()).contains("QweakTarget") ||
140  (thePostPV->GetName()).contains("QweakTarget"))) {
141  if (procName.compare("msc")==0)
142  myUserInfo->AddTodEMscIn(dEE);
143  else if(procName.compare("eIoni")==0)
144  myUserInfo->AddTodEIonIn(dEE);
145  else if(procName.compare("eBrem")==0)
147  }
148  }
149 
150  G4double theStepLength = theStep->GetStepLength();
151 
152  // if the z-position of the particle, theZ, is within
153  // the stepLength of RandomPositionZ,
154  // then kill the particle, and reset the origin vertex Z as theZ
155  if( fabs(theZ - RandomPositionZ)<=theStepLength && sqrt(theX*theX+theY*theY)<2.54*cm){
156 
157  std::vector< G4double > CrossSection;
158  for (int i = 0; i<16; i++) { CrossSection.push_back(0.0); }
159 
160  G4double WeightN, Q2, E_out, theta, phi;
161  G4double Asymmetry;
162  G4ThreeVector MomentumDirection = theTrack->GetMomentumDirection();
163  G4double E_in = theTrack->GetTotalEnergy()/MeV; //Event generator needs units of MeV
164 
165  // evaluate the MomentumDirection, E_out ..
166  myEvent->GetanEvent(E_in, CrossSection, WeightN, Q2, E_out, MomentumDirection, theta, phi, Asymmetry);
167 
168  //theTrack->SetKineticEnergy(E_out*MeV);
169  //theTrack->SetMomentumDirection(MomentumDirection);
170 
172 
173  // set track info
174  //theTrack->SetVertexPosition(thePosition);
175  //theTrack->SetVertexMomentumDirection(MomentumDirection);
176  //theTrack->SetVertexKineticEnergy(E_out);
177 
178  //fill user info
179  myUserInfo->StoreTrackID(theTrack->GetTrackID());
180  myUserInfo->StoreGlobalTime(theTrack->GetGlobalTime());
183  // myUserInfo->StoreOriginVertexPositionZ(theZ);
185 
186  myUserInfo->StoreOriginVertexMomentumDirectionX(MomentumDirection.getX());
187  myUserInfo->StoreOriginVertexMomentumDirectionY(MomentumDirection.getY());
188  myUserInfo->StoreOriginVertexMomentumDirectionZ(MomentumDirection.getZ());
189 
192  //myUserInfo->StoreOriginVertexKineticEnergy(theTrack->GetKineticEnergy());
193  //myUserInfo->StoreOriginVertexTotalEnergy(theTrack->GetTotalEnergy());
194 
195  myUserInfo->StorePreScatteringKineticEnergy(E_in - 0.511*MeV);
196  if(particleType==G4Electron::ElectronDefinition()){
197  myUserInfo->StoreOriginVertexKineticEnergy(E_out - 0.511*MeV);
198  // G4cout << "Stepping Action E_out:: " << (E_out - 0.511*MeV)/MeV << G4endl;
199  }
200  if(particleType==G4PionMinus::PionMinusDefinition()||particleType==G4PionPlus::PionPlusDefinition())
201  {
202  myUserInfo->StoreOriginVertexKineticEnergy(E_out - 139.57*MeV);
203  }
205 
206  myUserInfo->StorePrimaryQ2(Q2*0.000001); //in units of GeV^2
207  myUserInfo->StoreCrossSection(CrossSection[0]);
209  myUserInfo->StoreCrossSectionBornTotal (CrossSection[1]);
210  myUserInfo->StoreCrossSectionBornInelastic(CrossSection[2]);
211  myUserInfo->StoreCrossSectionBornQE (CrossSection[3]);
212  myUserInfo->StoreCrossSectionRadTotal (CrossSection[4]);
213  myUserInfo->StoreCrossSectionRadElastic (CrossSection[5]);
214  myUserInfo->StoreCrossSectionRadQE (CrossSection[6]);
215  myUserInfo->StoreCrossSectionRadDIS (CrossSection[7]);
216  myUserInfo->StoreCrossSectionRadTotalIntOnly (CrossSection[8]);
218  myUserInfo->StoreCrossSectionRadQEIntOnly (CrossSection[10]);
219  myUserInfo->StoreCrossSectionRadDISIntOnly (CrossSection[11]);
220  myUserInfo->StoreCrossSectionRadElasticPeak (CrossSection[12]);
221  myUserInfo->StoreAsymmetry ( Asymmetry );
222  //myUserInfo->StorePrimaryEventNumber(myEventCounter);
224  myUserInfo->StorePDGcode(theTrack->GetDefinition()->GetPDGEncoding());
225 
226  // print the stored values
227  // G4cout << "*********** myEventCounter = " << myEventCounter << G4endl;
228  // G4cout << "Pre scat KE:: "<<E_in<< G4endl;
229  // G4cout << "Org ver scat KE:: "<<E_out<< G4endl;
230 
231  // myUserInfo->Print();
232  theTrack->SetTrackStatus(fStopAndKill);
233  } // end of if( fabs(theZ - RandomPositionZ)<=theStepLength && sqrt(theX*theX+theY*theY)<2.54*cm)
234  } // end of if( myUserInfo->EvtGenStatus == 0){ // true for primary event # = 0,2,4,...
235  } // end of if(myUserInfo->GetPrimaryEventNumber() %2!=0)
236  else{
237 
238  // need this here, else the track might get out of physical volume,
239  // and crash the program
240  if(!thePostPV || !thePrePV) return;
241 
242  // various energy losses at the target
243  if(((thePrePV->GetName()).contains("QweakTarget") ||
244  (thePostPV->GetName()).contains("QweakTarget"))) {
245  if (procName.compare("msc")==0)
247  else if(procName.compare("eIoni")==0)
249  else if(procName.compare("eBrem")==0)
250  myUserInfo->AddTodEBremOut(dEE);
251 
252  }
253  }
254  }
255 
256 //now this is handled in the TrackingAction with the control of the tracking flag
257 // else //secondary, umcomment to disregard all secondaries to speed up the primary particle simulation
258 // {
259 // theTrack->SetTrackStatus(fStopAndKill);
260 // return;
261 // }
262 
263 //jpan@nuclear.uwinnipeg
264 //kill a track if it is in collimators or shielding wall
265 // if(thePrePV->GetName()=="CollimatorHousing" || thePrePV->GetName()=="ShieldingWallHousing"){
266 // theTrack->SetTrackStatus(fStopAndKill); return;
267 // }
268 
269 // Peiqing : storing secondary information makes the simulation very slow
270 // I commented them out. Uncomment them only in case we really
271 // want to spend time to study the details of secondaries.
272 //
273 // QweakSimTrackInformation* info = (QweakSimTrackInformation*)(theTrack->GetUserInformation());
274 //
275 // for(G4int i = GetTrackVectorSize()-nSecTotal; i < GetTrackVectorSize(); i++)
276 // {
277 // if((*fSecondary)[i]->GetUserInformation()==0)
278 // {
279 // QweakSimTrackInformation* infoNew = new QweakSimTrackInformation(info);
280 //
281 // infoNew->StoreParticleDefinition(GetSecondaryParticleDefinition(i));
282 // infoNew->StoreParentEnergy(theTrack->GetTotalEnergy());
283 // infoNew->StorePrimaryKineticEnergy(GetSecondaryParticleKineticEnergy(i));
284 // infoNew->StoreCerenkovHitEnergy(-1,-1.0*MeV);
285 // infoNew->StoreCreatorProcess(GetSecondaryCreatorProcessName(i));
286 // infoNew->StoreOriginVertex(GetSecondaryParticleOrigin(i));
287 // (*fSecondary)[i]->SetUserInformation(infoNew);
288 // }
289 //
290 // if(particleType==G4Electron::ElectronDefinition() && theTrack->GetParentID() == 0 &&
291 // // !strcmp(thePrePV->GetName(),"CerenkovDetector_Physical")){
292 // (!strcmp(thePrePV->GetName(),"QuartzBar_PhysicalRight") || !strcmp(thePrePV->GetName(),"QuartzBar_PhysicalLeft") ))
293 // {
294 // if(GetSecondaryParticleDefinition(i) == G4OpticalPhoton::OpticalPhotonDefinition() &&
295 // GetSecondaryParticleTotalEnergy(i)/eV <= 4.9594)
296 // {
297 // myUserInfo->IncrementCerenkovOpticalPhotonCount();
298 // myUserInfo->StoreCerenkovPhotonEnergy(GetSecondaryParticleTotalEnergy(i));
299 // }
300 // }
301 // }
302 
303 
304 //jpan@nuclear.uwinnipeg.ca
305 // commented out the cout statements for speeding up
306 
307 // G4cout << "Particle Name = " << particleType->GetParticleName() << G4endl;
308 
309 // if(!strcmp(thePrePV->GetName(),"CerenkovDetector_Physical")){
310 //!strcmp(thePrePV->GetName(),"LightGuide_PhysicalRight") || !strcmp(thePrePV->GetName(),"LightGuide_PhysicalLeft")
311 
312 // Peiqing: commented out the followings for speeding up
313 //
314 // if(!strcmp(thePrePV->GetName(),"QuartzBar_PhysicalRight") || !strcmp(thePrePV->GetName(),"QuartzBar_PhysicalLeft"))
315 // {
316 //
317 // myUserInfo->AddCerenkovEnergyDeposit(theStep->GetTotalEnergyDeposit());
318 //
319 // if(theTrack->GetParentID() > 0 && (particleType==G4Electron::ElectronDefinition() ||
320 // particleType==G4Positron::PositronDefinition() ||
321 // particleType==G4Gamma::GammaDefinition())){
322 //
323 // // if(!strcmp(myUserInfo->GetStoredStepVolumeName(),"CerenkovContainer_Physical") &&
324 // // !strcmp(thePrePV->GetName(),"CerenkovDetector_Physical")
325 // // ){
326 //
327 // if(!strcmp(myUserInfo->GetStoredStepVolumeName(),"ActiveArea_Physical")){
328 //
329 // myUserInfo->StoreCerenkovSecondaryParticleInfo(theTrack->GetVertexPosition(),
330 // theTrack->GetMomentum(),
331 // theTrack->GetTotalEnergy(),
332 // charge);
333 // }
334 // }
335 // }
336 
337  if(particleType==G4Electron::ElectronDefinition())
338  {
339  G4ThreeVector worldPos = thePrePoint->GetPosition();
340  G4ThreeVector localPos = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPos);
341 
342  if((!strcmp(myUserInfo->GetStoredStepVolumeName(),"QuartzBar_PhysicalRight") ||
343  !strcmp(myUserInfo->GetStoredStepVolumeName(),"QuartzBar_PhysicalLeft")) &&
344  (!strcmp(thePrePV->GetName(),"ActiveArea_Physical") ||
345  !strcmp(thePrePV->GetName(),"CerenkovMasterContainer_Physical")))
346  {
347  // if(!strcmp(myUserInfo->GetStoredStepVolumeName(),"CerenkovDetector_Physical") &&
348  // !strcmp(thePrePV->GetName(),"CerenkovContainer_Physical")){
350  }
351  }
352 
353 // QweakSimTrackInformation *TrackInfo = (QweakSimTrackInformation*)theTrack->GetUserInformation();
354 
355 //pqwang: optical photon process
356 
357  G4OpBoundaryProcessStatus boundaryStatus=Undefined;
358  static G4OpBoundaryProcess* boundary=NULL;
359 
360  if(!boundary){
361 // G4ProcessManager* pm
362 // = theStep->GetTrack()->GetDefinition()->GetProcessManager();
363  G4int nprocesses = pm->GetProcessListLength();
364  G4ProcessVector* pv = pm->GetProcessList();
365  for(G4int i=0;i<nprocesses;i++){
366  if((*pv)[i]->GetProcessName()=="OpBoundary"){
367  boundary = (G4OpBoundaryProcess*)(*pv)[i];
368  break;
369  }
370  }
371  }
372 
373  if(!thePostPV){
374  return;
375  }
376 
377  if(particleType==G4OpticalPhoton::OpticalPhotonDefinition()){
378  boundaryStatus=boundary->GetStatus();
379  if(thePostPoint->GetStepStatus()==fGeomBoundary){
380  switch(boundaryStatus){
381  case Absorption:
382  {
383  //std::cout<<"Absorption"<<std::endl;
384  break;
385  }
386  case Detection:
387  {
388  //std::cout<<"Detected a photon."<<std::endl;
389  G4SDManager* SDman = G4SDManager::GetSDMpointer();
390  if(strcmp(thePostPV->GetName(),"PMTOnlyCathode_Physical") == 0) {
391  //std::cout << "Optical Photon hit pmtonl cathode" << std::endl;
392  G4String pmtonlSDName="PMTOnly_PMTSD";
393  QweakSimPMTOnly_PMTSD* pmtonlSD = (QweakSimPMTOnly_PMTSD*)SDman->FindSensitiveDetector(pmtonlSDName);
394  if(pmtonlSD) {
395  //myUserInfo->GetCurrentPMTHit()->SetHitValid(True);
396  pmtonlSD->ProcessHits_constStep(theStep,NULL);
397  theTrack->SetTrackStatus(fStopAndKill);
398  }
399  break;
400  } else {
401  G4String CerenkovSDName="/CerenkovPMTSD";
402  QweakSimCerenkov_PMTSD* cerenkovpmtSD = (QweakSimCerenkov_PMTSD*)SDman->FindSensitiveDetector(CerenkovSDName);
403  if(cerenkovpmtSD) {
404  //myUserInfo->GetCurrentPMTHit()->SetHitValid(True);
405  cerenkovpmtSD->ProcessHits_constStep(theStep,NULL);
406  theTrack->SetTrackStatus(fStopAndKill);
407  }
408  break;
409  }
410  }
411  case FresnelReflection:
412  {
413  //std::cout<<"FresnelReflection"<<std::endl;
414  break;
415  }
416  case TotalInternalReflection:
417  {
418  //std::cout<<"TotalInternalReflection"<<std::endl;
419  break;
420  }
421  case SpikeReflection:
422  {
423  //std::cout<<"SpikeReflection"<<std::endl;
424  break;
425  }
426  default:
427  {
428  //std::cout<<"Undefined"<<std::endl;
429  break;
430  }
431  }
432  }
433  } // end of optical photon process
434 
435 
436 // if(particleType==G4OpticalPhoton::OpticalPhotonDefinition()){
437 //
438 // if(!strcmp(myUserInfo->GetStoredStepVolumeName(),"PMTEntranceWindow_Physical")){
439 // if(!strcmp(thePrePV->GetName(),"Cathode_Physical")){
440 //
441 // myUserInfo->GetCurrentPMTHit()->SetHitValid(True);
442 // // theTrack->SetTrackStatus(fStopAndKill);
443 // }
444 // }
445 // }
446 
447  myUserInfo->StoreStepVolumeName(thePrePV->GetName());
448 
449 //======================================================================
450 // Stolen from GATE code:
451 //
452 // In a few random cases, a particle gets 'stuck' in an
453 // an infinite loop in the geometry. It then oscillates until GATE
454 // crashes on some out-of-memory error.
455 // To prevent this from happening, I've added below a quick fix where
456 // particles get killed when their step number gets absurdely high
457 
458 //jpan@nuclear.uwinnipeg.ca
459 //The following codes cause the electron tracks killed at the edge of magnetic field
460 //This must be avoid by increasing the 10000 steps to 100000.
461 
462  if ( theStep->GetTrack()->GetCurrentStepNumber() > 100000 )
463  theStep->GetTrack()->SetTrackStatus(fStopAndKill);
464  return;
465 } // end of QweakSimSteppingAction::UserSteppingAction
466 
467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
468 
470 {
471  if(!fSecondary) return -1;
472 
473  return (*fSecondary).size() - GetTotalNumOfSecondaries();
474 }
475 
476 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
478 {
479  if(!fSecondary) return 0;
480  return (*fSecondary).size();
481 }
482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
484 {
485  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return NULL;
486 
487  return (*fSecondary)[idx]->GetDefinition();
488 }
489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
491 {
492  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return "undefined";
493 
494  return (*fSecondary)[idx]->GetDefinition()->GetParticleName();
495 }
496 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
498 {
499  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return -1;
500 
501  return (*fSecondary)[idx]->GetTotalEnergy();
502 }
503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
505 {
506  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return -1;
507 
508  return (*fSecondary)[idx]->GetKineticEnergy();
509 }
510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
512 {
513  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return 1e6;
514 
515  return (*fSecondary)[idx]->GetPosition().x();
516 }
517 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
519 {
520  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return 1e6;
521 
522  return (*fSecondary)[idx]->GetPosition().y();
523 }
524 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
526 {
527  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return 1e6;
528 
529  return (*fSecondary)[idx]->GetPosition().z();
530 }
531 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
533 {
534  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return G4ThreeVector(1e6,1e6,1e6);
535 
536  return (*fSecondary)[idx]->GetPosition();
537 }
538 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
540 {
541  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return G4ThreeVector(1e6,1e6,1e6);
542 
543  return (*fSecondary)[idx]->GetMomentumDirection();
544 }
545 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
547 {
548  if(!fSecondary || idx >= GetTrackVectorSize() || idx < GetTrackVectorStartIndex()) return "undefined";
549  return (*fSecondary)[idx]->GetCreatorProcess()->GetProcessName();
550 }
551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
552 
553 
554 
555  //----------------------------------------------------------------------------------------------------------------------------------------------------------
556  // here we are only interested if the particle is crossing a boundary
557  // If you want to identify the first step in a volume: pick fGeomBoundary status in preStepPoint
558  // If you want to identify a step out of a volume : pick fGeomBoundary status in postStepPoint
559 
560 // if( (theStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) && (strcmp( theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName(),"CerenkovDetector_Physical") == 0) ) //OK
561 // //if( (theStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) && (strcmp( theStep->GetPostStepPoint()->GetPhysicalVolume()->GetName(),"CerenkovDetector_Physical") == 0) ) // crash
562 // {
563 // G4cout << "------------------------" << G4endl;
564 // //G4cout << "PreStepVolumeName = " << theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << G4endl;
565 // //G4cout << "PostStepVolumeName = " << theStep->GetPostStepPoint()->GetPhysicalVolume()->GetName() << G4endl;
566 // G4cout << "------------------------" << G4endl;
567 
568 // // QweakSimTrackInformation* trackInfo;
569 // // trackInfo = (QweakSimTrackInformation*)(theTrack->GetUserInformation());
570 // // trackInfo->SetTrackingStatus(1);
571 // // trackInfo->SetImpactTrackInformationForCerenkov(theTrack);
572 
573 // }
574 
575  //----------------------------------------------------------------------------------------------------------------------------------------------------------
576 
577 // if(particleType==G4Electron::ElectronDefinition())
578 // {
579 // G4ThreeVector worldPos = thePrePoint->GetPosition();
580 // G4ThreeVector localPos = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPos);
581 
582 // if( (myUserInfo->GetStoredStepVolumeName().compare("CerenkovDetector_Physical")==0) &&
583 // (thePrePV->GetName().compare("CerenkovContainer_Physical")==0) )
584 // {
585 // myUserInfo->StoreLocalCerenkovExitPosition(localPos);
586 // }
587 // }
588  //----------------------------------------------------------------------------------------------------------------------------------------------------------
589 
590 // if(particleType==G4OpticalPhoton::OpticalPhotonDefinition()){
591 // if(!strcmp(myUserInfo->GetStoredStepVolumeName(),"PMTEntranceWindow_Physical")){
592 //
593 //
594 // if(!strcmp(thePrePV->GetName(),"Cathode_Physical")){
595 //
596 // QweakSimTrackInformation *TrackInfo = (QweakSimTrackInformation*)theTrack->GetUserInformation();
597 // G4cout << "Particle History For This Photon" << G4endl;
598 // G4cout << "Event ID = " << myUserInfo->GetPrimaryEventNumber() << G4endl;
599 // G4cout << "Hit ID = " << myUserInfo->GetCurrentPMTHit()->GetHitID() << G4endl;
600 // for(int i = 0; i < TrackInfo->GetParticleHistoryLength(); i++){
601 // G4cout << "Particle "<< i <<" = " << TrackInfo->GetParticleDefinitionAtIndex(i)->GetParticleName();
602 // if(i == TrackInfo->GetParticleHistoryLength()-1)
603 // G4cout << " Eng = " << theTrack->GetTotalEnergy()/MeV;
604 // else
605 // G4cout << " Eng = " << TrackInfo->GetParentEnergyAtIndex(i+1)/MeV;
606 //
607 // G4cout << " Cerenkov Hit Energy = " << TrackInfo->GetCerenkovHitEnergyAtIndex(i) << G4endl;
608 // }
609 // // G4int index = 0;
610 // // if(TrackInfo->GetParticleHistoryLength() < 3)
611 // // // myUserInfo->SetPhotonFromPrimary(TrackInfo->GetParticleDefinitionAtIndex(index));
612 // // else{
613 // // index = TrackInfo->GetParticleHistoryLength()-3;
614 // // // myUserInfo->SetPhotonFromParticle(TrackInfo->GetParticleDefinitionAtIndex(index));
615 // // }
616 //
617 // myUserInfo->GetCurrentPMTHit()->SetHitValid(True);
618 // // theTrack->SetTrackStatus(fStopAndKill);
619 // }
620 // }
621 // }
622 
623 
624 
625 
626 
627 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
628 
629 
void GetanEvent(G4double E_in, std::vector< G4double > &CrossSection, G4double &weight_n, G4double &Q2, G4double &E_out, G4ThreeVector &MomentumDirection, G4double &theta, G4double &phi, G4double &Asymmetry)
G4String GetSecondaryCreatorProcessName(G4int idx)
void StoreCrossSectionRadTotal(G4double cs)
void UserSteppingAction(const G4Step *)
G4ThreeVector GetSecondaryParticleMomentum(G4int idx)
G4double GetOriginVertexPositionZ() const
void StoreOriginVertexPhiAngle(G4double phi)
void StoreOriginVertexPositionX(G4double vx)
void StoreCrossSectionRadDISIntOnly(G4double cs)
void StoreCrossSectionRadElasticPeak(G4double cs)
void StoreCrossSectionBornTotal(G4double cs)
void StoreOriginVertexPositionY(G4double vy)
void StoreStepVolumeName(G4String name)
void StoreGlobalTime(G4double gtime)
void StoreAsymmetry(G4double asym)
void StoreOriginVertexTotalEnergy(G4double etot)
void StoreCrossSectionBornQE(G4double cs)
G4double GetSecondaryParticleTotalEnergy(G4int idx)
G4String GetSecondaryParticleName(G4int idx)
void StoreCrossSectionRadQEIntOnly(G4double cs)
void StoreOriginVertexPositionZ(G4double vz)
void StoreCrossSectionRadTotalIntOnly(G4double cs)
void StoreOriginVertexMomentumDirectionZ(G4double vz)
void StoreOriginVertexMomentumDirectionX(G4double vx)
void StoreCrossSectionRadDIS(G4double cs)
void StoreCrossSectionBornInelastic(G4double cs)
void StoreCrossSectionRadElastic(G4double cs)
void StorePreScatteringKineticEnergy(G4double ekin)
void StoreOriginVertexKineticEnergy(G4double ekin)
QweakSimSteppingAction(QweakSimUserInformation *myUInfo, QweakSimEPEvent *myEPEvent)
void StoreCrossSectionWeight(G4double csw)
G4ParticleDefinition * GetSecondaryParticleDefinition(G4int idx)
void StoreCrossSectionRadQE(G4double cs)
G4double GetSecondaryParticleXOrigin(G4int idx)
void StoreCrossSectionRadElasticIntOnly(G4double cs)
G4bool ProcessHits_constStep(const G4Step *aStep, G4TouchableHistory *ROhist)
G4double GetSecondaryParticleYOrigin(G4int idx)
QweakSimUserInformation * myUserInfo
G4double GetSecondaryParticleKineticEnergy(G4int idx)
G4ThreeVector GetSecondaryParticleOrigin(G4int idx)
G4double GetSecondaryParticleZOrigin(G4int idx)
void StoreOriginVertexThetaAngle(G4double theta)
G4bool ProcessHits_constStep(const G4Step *, G4TouchableHistory *)
void StoreOriginVertexMomentumDirectionY(G4double vy)
void StoreLocalCerenkovExitPosition(G4ThreeVector ep)