LBNESteppingAction.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 // LBNESteppingAction.cc
3 // $Id: LBNESteppingAction.cc,v 1.1.1.1.2.16 2013/12/02 20:24:52 lebrun Exp $
4 //----------------------------------------------------------------------
5 
6 //C++
7 #include <string>
8 
9 #include "LBNESteppingAction.hh"
10 
11 #include "G4ProcessManager.hh"
12 #include "G4SteppingManager.hh"
13 #include "G4Track.hh"
14 #include "G4ParticleDefinition.hh"
15 #include "G4ParticleTypes.hh"
16 #include "G4TrajectoryContainer.hh"
17 #include "G4RunManager.hh"
18 #include "G4VPhysicalVolume.hh"
19 #include "G4Event.hh"
20 #include "G4VTouchable.hh"
21 #include "G4TouchableHistory.hh"
22 #include "G4EventManager.hh"
23 #include "G4FieldManager.hh"
24 #include "G4TrackingManager.hh"
25 
26 #include "LBNETrackInformation.hh"
27 #include "LBNETrajectory.hh"
28 #include "LBNETrackInformation.hh"
29 #include "LBNEAnalysis.hh"
30 #include "LBNEEventAction.hh"
31 #include "LBNERunManager.hh"
32 #include "LBNEMagneticField.hh"
34 #include "LBNEVolumePlacements.hh"
35 #include "LBNEQuickPiToNu.hh"
37 #include "LBNFDeDxMap.hh"
38 #include "TFile.h"
39 #include "TTree.h"
40 
41 //----------------------------------------------------------------------
43 {
44  pRunManager=(LBNERunManager*)LBNERunManager::GetRunManager();
45  fKillTrackingThreshold = 0.050*CLHEP::GeV;
46  if (pRunManager->GetVerboseLevel() > 0) {
47  std::cerr << " LBNESteppingAction::LBNESteppingAction called ... " << std::endl;
48  }
50  fEvtIdPrevious = -1;
52  fStudyGeantinoMode=G4String("none");
53  fKeyVolumeForOutput=G4String("blank");
54  fKeyVolumeForOutputTo=G4String("blank");
63  fGenerateVoxelsZScale = 25.0;
64  fGenerateVoxelsZOffset = 5000.0;
65  fGenerateVoxelsZMax = 25000.0;
67  fGenerateVoxelsRMax = 2000.0;
69  fOutMuonSculptedAbsorberFilename="MuonFluxAtSculptedAbsorber";
70  //
71  // April 2017 study for MARS, cmpare meson rate..
72  //
73  bool doSergeiStriganovApr2017 = false;
74  if (doSergeiStriganovApr2017) {
75  std::string dirNameTmp("/dune/data/users/lebrun/LBNFMARSApril2017/");
76  std::string fNameTgt(dirNameTmp); fNameTgt += std::string("DumpOptimizedEngineered_Tgt.txt");
77  std::ofstream *fOutTargetSergei = new std::ofstream(fNameTgt.c_str());
78  fOutPtrsForMarsCmpApr2017.push_back(fOutTargetSergei);
79  //
80  std::string fNameHornA(dirNameTmp); fNameHornA += std::string("DumpOptimizedEngineered_HornA.txt");
81  std::ofstream *fOutHornASergei = new std::ofstream(fNameHornA.c_str());
82  fOutPtrsForMarsCmpApr2017.push_back(fOutHornASergei);
83  //
84  std::string fNameHornB(dirNameTmp); fNameHornB += std::string("DumpOptimizedEngineered_HornB.txt");
85  std::ofstream *fOutHornBSergei = new std::ofstream(fNameHornB.c_str());
86  fOutPtrsForMarsCmpApr2017.push_back(fOutHornBSergei);
87  //
88  std::string fNameHornC(dirNameTmp); fNameHornC += std::string("DumpOptimizedEngineered_HornC.txt");
89  std::ofstream *fOutHornCSergei = new std::ofstream(fNameHornC.c_str());
90  fOutPtrsForMarsCmpApr2017.push_back(fOutHornCSergei);
91  for (size_t k=0; k!= fOutPtrsForMarsCmpApr2017.size(); k++) {
92  (*fOutPtrsForMarsCmpApr2017[k]) << " evt trId pdg x y z px py pz ek " << std::endl;
93  }
94  }
95 }
96 //----------------------------------------------------------------------
98 {
99  if (fOutStudy.is_open()) {
100  if (fStudyGeantinoMode.find("RZVoxels") != std::string::npos) {
101  std::cerr << " Dumping in for RZVoxels, num Voxels = " << fRZVoxelsData.size() << std::endl;
103  itMRZ != fRZVoxelsData.end(); itMRZ++) {
104  int ik = itMRZ->first;
105  int iz = ik/fGenerateVoxelsIRKMax;
106  int ir = ik - iz*fGenerateVoxelsIRKMax;
107  const double z = fGenerateVoxelsZScale*iz - fGenerateVoxelsZOffset; //see in GenerateRZVoxel..
108  const double r = std::exp(static_cast<double>(ir)/fGenerateVoxelsRScale);
109  fOutStudy << " " << ik << " " << r << " " << z << " " << itMRZ->second << std::endl;
110  }
111  }
112  fOutStudy.close();
113  }
114  if (fOutStepStudy.is_open()) fOutStepStudy.close();
115  if(makeSteppingTuple) {
116  steppingTupleFile->cd();
117  steppingTuple->Write();
118  steppingTupleFile->Close();
119  }
120 
121  delete pMessenger;
122  // delete
123  //
124  for (size_t k=0; k!= fOutPtrsForMarsCmpApr2017.size(); k++) {
125  fOutPtrsForMarsCmpApr2017[k]->close();
127  }
128 
129 }
130 
131 //----------------------------------------------------------------------
132 void LBNESteppingAction::UserSteppingAction(const G4Step * theStep)
133 {
134  //
135  // Temporary timing study: kill the track outside the target region,
136  //
137 // {
138 // G4StepPoint *prePtr = theStep->GetPostStepPoint();
139 // G4ThreeVector posT = prePtr->GetPosition();
140 // const double rrSq = posT[0]*posT[0] + posT[1]*posT[1];
141 // G4Track *theTrack = theStep->GetTrack();
142 // if (rrSq > 310.) {
143 // if (posT[2] < 2500.) theTrack->SetTrackStatus(fStopAndKill);
144 // } else {
145 // if (posT[2] > 2500.) theTrack->SetTrackStatus(fStopAndKill);
146 // }
147 // }
148  // Killing tracks that are stuck on a boundary crossing
149  // This should not happen if the geometry has no overlap and is strictly
150  // correct. For the nominal case, except for the hadron absorver,
151  // no overlaps have been detected at construction stage.
152  // However, the tests are "weak", i.e., do not explore all corners.
153  // Let us protect our-self here again infinite loop
154  //
155  if (theStep->GetStepLength() < 1.0e-15) {
157  if (fNConsecutivSmallSteps > 5000) {
158  G4StepPoint *prePtr = theStep->GetPreStepPoint();
159  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
160  std::cerr << " LBNESteppingAction, too man consecutive small steps " << std::endl
161  << ".. from " << volPre->GetName()
162  << " detected at " << prePtr->GetPosition() <<
163  " last step length " << theStep->GetStepLength() << std::endl;
164  G4Track* aTrack= theStep->GetTrack();
165  std::cerr << " ... for track " << aTrack->GetTrackID() << " At " << aTrack->GetPosition()
166  << " P " << aTrack->GetMomentumDirection() << " E " << aTrack->GetTotalEnergy()
167  << " name " << aTrack->GetParticleDefinition()->GetParticleName() << std::endl;
168  aTrack->SetTrackStatus(fStopAndKill);
169  std::cerr << " Track killed .... " << std::endl;
171  }
172 
173  } else fNConsecutivSmallSteps=0;
174  // Last check on memory explosion: No more than 500 mgB usage. Check than our step counter (can't trust G4 anylonger)
175  // gets too large..
176  // G4 10 seems be a bit more demanding. Relax these limits.
178  if (fNumStepsCurrentTrack > 200000) {
179  std::cerr << " LBNESteppingAction... Anomalously large number of steps.. kill this event.. " << std::endl;
180  G4Exception("LBNESteppingAction", " ", EventMustBeAborted, " Too many steps.. ");
181  return;
182 // const LBNERunAction * runAct = dynamic_cast<const LBNERunAction*>(pRunManager->GetUserRunAction());
183 // runAct->CheckMemoryUsage(1000000);
184  }
185  //
186  // June 2017: we were asked to provide dedx in the horns.
187  //
188  const LBNERunAction * runAct = dynamic_cast<const LBNERunAction*>(pRunManager->GetUserRunAction());
189  if (runAct->GetDoComputeEDepInHorns()) {
190  LBNFDeDxMap *aMapDedX = runAct->GetG4MARSDeDxMap();
191  if (aMapDedX != 0) aMapDedX->fill(theStep);
192  G4StepPoint *postPtr = theStep->GetPreStepPoint();
193  // We don't track past the horns, not interested in neutrino spectrum in this slow tedious mode..
194  // Assume extended chase. In fact, I am asked to do it only for horn A
195  if (postPtr->GetPosition()[2] > 3000.) theStep->GetTrack()->SetTrackStatus(fStopAndKill);
196  }
197  //
198  // Code bloat, April 2017..
199  //
200  if (fOutPtrsForMarsCmpApr2017.size() > 0) this->CheckHadronsMarsCmpApr2017(theStep);
201 
202  // Debugging the bad neutrons in 4.9.6.p02
203 // G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
204 
205 // if (evtno == 32135) {
206 // const G4Track* aTrack= theStep->GetTrack();
207 // if (aTrack->GetTrackID() == 79) {
208 // G4StepPoint *prePtr = theStep->GetPreStepPoint();
209 // G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
210 // std::cout << " LBNESteppingAction::StudyPropagation, bad track 79 " << std::endl
211 // << ".. from " << volPre->GetName()
212 // << " detected at " << prePtr->GetPosition() <<
213 // " length " << theStep->GetStepLength() << std::endl;
214 // Found :
215 // LBNESteppingAction::StudyPropagation, bad track 79
216 //.. from Horn1IOTransCont detected at (31.0918,-3.71985,109.441) length 1.77636e-15
217 // LBNESteppingAction::StudyPropagation, bad track 79
218 //.. from Horn1IOTransInnerPart3 detected at (31.0918,-3.71985,109.441) length 0
219 // LBNESteppingAction::StudyPropagation, bad track 79
220 //.. from Horn1IOTransCont detected at (31.0918,-3.71985,109.441) length 1.77636e-15
221 // For ever...
222 //
223 // }
224 // }
225 
226  int verboseLevel =
227  G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager()->GetverboseLevel();
228  if(verboseLevel > 3)
229  {
230  G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
231  std::cout << "Event " << evtno << ": LBNESteppingAction::UserSteppingAction() Called." << std::endl;
232  }
233  if((fStudyGeantinoMode.find("Mag") != std::string::npos)
234  && (pRunManager->GetCurrentEvent()->GetEventID() < 5))
235  this->dumpStepCheckVolumeAndFields(theStep);
236 
238 
240 
244  if (fOutStudy.is_open()) {
246 // std::cerr << " Stepping .... doStudyParticleThroughHorns.. And quit " << std::endl; exit(2);
247  StudyParticleThroughHorns(theStep);
248  }
249  else {
250 
251  if (fStudyGeantinoMode.find("Absorb") != std::string::npos) StudyAbsorption(theStep);
252  if (fStudyGeantinoMode.find("Propa") != std::string::npos) StudyPropagation(theStep);
253  if (fStudyGeantinoMode.find("PropCO") != std::string::npos) StudyCheckOverlap(theStep);
254  if (fStudyGeantinoMode.find("MagTilt") != std::string::npos) StudyCheckMagneticTilts(theStep);
255  if (fStudyGeantinoMode.find("RZVoxels") != std::string::npos) GenerateVolumeCrossingsRZMap(theStep);
256  }
257  }
258  if (makeSteppingTuple)
259  StudyAbsorption(theStep);
261 
262  // Add energy loss if requested..
263  const LBNERunAction *runUserAction = reinterpret_cast<const LBNERunAction *>(G4RunManager::GetRunManager()->GetUserRunAction());
264  if (runUserAction->GetDoComputeEDepInGraphite()) {
265  const std::string tgt("Target");
266  const std::string tg = theStep->GetTrack()->GetMaterial()->GetName();
267  if (tg == tgt) fEnergyDepInGraphite += theStep->GetTotalEnergyDeposit(); // Not true for all targets!!!! Must have the right material name.
268  }
269  if (runUserAction->GetDoComputeEDepInArgonGas()) {
270  const G4StepPoint *prePt = theStep->GetPreStepPoint();
271  if ((prePt != 0) && (prePt->GetPhysicalVolume() != 0)) {
272  const G4LogicalVolume *aVol = prePt->GetPhysicalVolume()->GetLogicalVolume();
273  std::string theMatName(aVol->GetMaterial()->GetName());
274  if (theMatName == std::string("Argon")) {
275  std::string theVolName(aVol->GetName());
276  if (theVolName.find("orn1") != std::string::npos) {
277  fEnergyDepInArgonGasHorn1 += theStep->GetTotalEnergyDeposit();
278 // std::cerr << " fEnergyDepInArgonGasHorn1 now " << fEnergyDepInArgonGasHorn1 << std::endl;
279  } else if (theVolName.find("orn2") != std::string::npos) {
280  fEnergyDepInArgonGasHorn2 += theStep->GetTotalEnergyDeposit();
281 // std::cerr << " fEnergyDepInArgonGasHorn2 now " << fEnergyDepInArgonGasHorn2 << std::endl;
282  } else {
283  std::string msg(" Unknow volume that contains Argon, volName "); msg+=theVolName;
284  G4Exception("LBNESteppingaction", "", FatalException, msg.c_str());
285  }
286  }
287  }
288  }
289 
290  //
291  // Adding checks for Perfect Focsuing..
292  //
293 // const G4Track *theTrack = theStep->GetTrack();
294 // const G4DynamicParticle *pDef = theTrack->GetDynamicParticle();
295 // const int pId = pDef->GetDefinition()->GetPDGEncoding();
296 // if (((std::abs(pId) == 211) || (std::abs(pId) == 321)) && (theStep->GetPreStepPoint()->GetPosition()[2] < 2000.*CLHEP::mm)) {
297 // const LBNETrackInformation *oldTrInfo = reinterpret_cast<const LBNETrackInformation*>(theTrack->GetUserInformation());
298 // std::cerr << " Pion or Kaon, Id " << theTrack->GetTrackID() << " at " << theStep->GetPreStepPoint()->GetPosition()
299 // << " Momentum " << theTrack->GetMomentum() << " Flag " << oldTrInfo->GetPFFlag() << std::endl;
300 // }
301 
303  //---tracking planes
304  TrkPlnLogical =
305  ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
306  ->GetUserDetectorConstruction())->GetMuonDetectorLogical();
307 
308  G4StepPoint *postStep = theStep->GetPostStepPoint();
309  if(postStep->GetPhysicalVolume()){
310  G4LogicalVolume *postStepVolume =
311  postStep->GetPhysicalVolume()->GetLogicalVolume();
312  if(postStepVolume){
313  // std::cout << postStepVolume->GetName() << std::endl;
314  if(postStepVolume == TrkPlnLogical){
316  }
317  }
318  }
319  }
320 
322 
323  CheckInTargetOutput(theStep);
324  }
325  /*
326 if(pRunManager->GetCreateTrkPlaneDPOutput()){
327 
328  CheckInTrackingDetectorDPPlane(theStep);
329 }
330 */
331 
332 // To stepping method to save the particle info that makes through the tracking plane at the end of Horn1 and horn 2. Amit Bashyal
334  //---tracking planes
336  ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
337  ->GetUserDetectorConstruction())->GetHorn1PlaneLogical();
338 
339  G4StepPoint *postStep = theStep->GetPostStepPoint();
340  if(postStep->GetPhysicalVolume()){
341  G4LogicalVolume *postStepVolume =
342  postStep->GetPhysicalVolume()->GetLogicalVolume();
343  if(postStepVolume){
344  // std::cout << postStepVolume->GetName() << std::endl;
345  if(postStepVolume == TrkPlnH1Logical){
346 
348  }
349  }
350  }
351 }
352 
354  //---tracking planes
356  ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
357  ->GetUserDetectorConstruction())->GetHorn2PlaneLogical();
358 
359  G4StepPoint *postStep = theStep->GetPostStepPoint();
360  if(postStep->GetPhysicalVolume()){
361  G4LogicalVolume *postStepVolume =
362  postStep->GetPhysicalVolume()->GetLogicalVolume();
363  if(postStepVolume){
364  // std::cout << postStepVolume->GetName() << std::endl;
365  if(postStepVolume == TrkPlnH2Logical){
367  }
368  }
369  }
370 }
371 //Although the tracking plane at the end (or start) of Decay Pipe has been created and the method to save the particle info is set, the Hadron monitor is likely to
372 // be in this position.
374  //---tracking planes
376  ((LBNEDetectorConstruction*)G4RunManager::GetRunManager()
377  ->GetUserDetectorConstruction())->GetDecayPipePlaneLogical();
378 
379  G4StepPoint *postStep = theStep->GetPostStepPoint();
380  if((postStep != 0) && (postStep->GetPhysicalVolume() != 0)){
381  G4LogicalVolume *postStepVolume =
382  postStep->GetPhysicalVolume()->GetLogicalVolume();
383  if(postStepVolume){
384  // std::cout << postStepVolume->GetName() << std::endl;
385  if(postStepVolume == TrkPlnDPLogical){
387  }
388  }
389  }
390 }
391 
392 //Sculpted absorber alcove tracking plane
395  }
396 
397 }
398 
399 
400 //More PPFX stuff...yea!!! A. Bashyal
401 //----------------------------------------------------------------------------------------------------------------
403 {
404 G4Track * aTrack = theStep->GetTrack();
405 if(theStep->GetPreStepPoint()->GetProcessDefinedStep()==NULL)return;
406  G4ThreeVector mom = G4ThreeVector(0.,0.,0.);
407  //Basically imported ditto from NumiSteppingAction for the sake of dk2nu ntuples..A Bashyal
408  G4int nSecAtRest = fpSteppingManager->GetfN2ndariesAtRestDoIt();
409  G4int nSecAlong = fpSteppingManager->GetfN2ndariesAlongStepDoIt();
410  G4int nSecPost = fpSteppingManager->GetfN2ndariesPostStepDoIt();
411  G4int nSecTotal = nSecAtRest+nSecAlong+nSecPost;
412  G4TrackVector* secVec = fpSteppingManager->GetfSecondary();
413  //std::cout<<(*secVec).size()<<" "<<nSecTotal<<std::endl;
414  //G4bool isShower = false;
415  LBNETrackInformation* trackInformation = (LBNETrackInformation*)(aTrack->GetUserInformation());
416  if(nSecTotal==0)return;
417  //if(nSecTotal>0) {
418 
419  for(size_t lp1=(*secVec).size()-nSecTotal; lp1<(*secVec).size(); ++lp1) {
420  G4Track* secTrack = (*secVec)[lp1];
421  // std::cout<<"trackid for secondary tracks "<<(*secVec)[lp1]->GetTrackID()<<std::endl;
422  G4ParticleDefinition* def = secTrack->GetDefinition();
423  //std::cout<<"LBNESteppingAction... "<<def<<std::endl;
424  // Showering particles (e+/-,gamma) are not interesting for neutrino production,
425  // skip them
426  if (def == G4Electron::Electron() || def == G4Positron::Positron() ||
427  def == G4Gamma::Gamma())continue;
428  mom = theStep->GetPreStepPoint()->GetMomentum();
429 
430  }
431  // if(!isShower) mom = theStep->GetPostStepPoint()->GetMomentum();
432 
433  // std::cout<<"LBNESteppingAction:: "<<mom<<std::endl;
434 
435  // LBNETrackInformation* trackInformation = (LBNETrackInformation*)(aTrack->GetUserInformation());
436  // std::cout<<"I am here SteppingAction"<<std::endl;
437  if (trackInformation != 0) {
438  // std::cout<<" Now inside the if statement LBNESteppingaction"<<std::endl;
439  //std::cout<<"LBNESteppingAction:: "<<mom<<std::endl;
440  // LBNETrackInformation* trackInformation = new LBNETrackInformation;
441  // assert(trackInformation != NULL);
442  trackInformation->SetParentMomentumAtThisProduction(mom);
443  //if(mom != G4ThreeVector(0.,0.,0.)) std::cout<<"checking the function "<<mom<<std::endl;
444  aTrack->SetUserInformation(trackInformation);
445  } else {
446  // std::cout<<"And the else statement.."<<std::endl;
447  LBNETrackInformation* newtrackInformation = new LBNETrackInformation();
448  //= dynamic_cast<LBNETrackInformation*>(trackInformation->GetUserInformation());
449  newtrackInformation->SetParentMomentumAtThisProduction(mom);
450  aTrack->SetUserInformation(newtrackInformation);
451 
452  }
453 
454  //}
455  // }
456 
457 }
458 
459 //----------------------------------------------------------------------
460 //----------------------------------------------------------------------
462 {
463 
464  const LBNERunAction *runUserAction = reinterpret_cast<const LBNERunAction*>(G4RunManager::GetRunManager()->GetUserRunAction());
465  if (runUserAction->GetDoComputeEDepInHorns()) return;
466  if (runUserAction->GetDoComputeEDepInGraphite()) {
467  const std::string tgt("Target");
468  const std::string tg = theStep->GetTrack()->GetMaterial()->GetName();
469  if (tgt == tg) return; // we let them go until they escape the target..
470  }
471  if (runUserAction->GetDoComputeEDepInArgonGas()) {
472  const G4StepPoint *prePt = theStep->GetPreStepPoint();
473  if (prePt->GetPhysicalVolume() !=0 && (prePt->GetPosition()[2] > 9900.)) {
474  std::string vName = theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
475  G4Track * theTrack = theStep->GetTrack();
476  if (vName.find("Tunnel") != std::string::npos) theTrack->SetTrackStatus(fStopAndKill); // We don't produce neutrino,
477  // just compute energy loss in Argon, to compute ionization rate...
478  return;
479  }
480  }
481 
483  if(theStep->GetPreStepPoint()->GetPhysicalVolume() != NULL &&
484  theStep->GetPostStepPoint()->GetPhysicalVolume()!=NULL) {
485 
486  std::string preStepPointName =
487  theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
488  std::string postStepPointName =
489  theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
490  if (preStepPointName=="HadronAbsorberSculpBackWall" || postStepPointName=="HadronAbsorberSculpBackWall" ||
491  preStepPointName=="HadronAbsorberSculpAirBuffer" || postStepPointName=="HadronAbsorberSculpAirBuffer" ||
492  preStepPointName=="MuonTrackingVol" || postStepPointName=="MuonTrackingVol" )
493  return;//Don't kill tracks in back wall, muon alcove or tracking plane
494 
495  }//if step volumes exist
496  }//if saving alcove output
497 
498  G4Track * theTrack = theStep->GetTrack();
499  G4ParticleDefinition * particleDefinition = theTrack->GetDefinition();
500  //
501  //kill the track if it is below the kill tracking threshold and it is NOT a neutrino
502  //
503  if ( (theTrack->GetKineticEnergy() < fKillTrackingThreshold) &&
504  (particleDefinition != G4NeutrinoE::NeutrinoEDefinition())&&
505  (particleDefinition != G4NeutrinoMu::NeutrinoMuDefinition()) &&
506  (particleDefinition != G4NeutrinoTau::NeutrinoTauDefinition()) &&
507  (particleDefinition != G4AntiNeutrinoE::AntiNeutrinoEDefinition()) &&
508  (particleDefinition != G4AntiNeutrinoMu::AntiNeutrinoMuDefinition()) &&
509  (particleDefinition != G4AntiNeutrinoTau::AntiNeutrinoTauDefinition()))
510  {
511  theTrack->SetTrackStatus(fStopAndKill);
512  }
513 
514 }
515 
516 
517 //----------------------------------------------------------------------
518 void LBNESteppingAction::CheckDecay(const G4Step * theStep)
519 {
520  G4Track * theTrack = theStep->GetTrack();
521  G4ParticleDefinition * particleDefinition = theTrack->GetDefinition();
522 
523 
524  // Check if the Pi+, Pi-, K+, K-, K0L, mu+ or mu- decayed and set Ndecay code:
525  // 1 K0L -> nu_e pi- e+
526  // 2 K0L -> anti_nu_e pi+ e-
527  // 3 K0L -> nu_mu pi- mu+
528  // 4 K0L -> anti_nu_mu pi+ mu-
529  // 5 K+ -> nu_mu mu+
530  // 6 K+ -> nu_e pi0 e+
531  // 7 K+ -> nu_mu pi0 mu+
532  // 8 K- -> anti_nu_mu mu-
533  // 9 K- -> anti_nu_e pi0 e-
534  // 10 K- -> anti_nu_mu pi0 mu-
535  // 11 mu+ -> anti_nu_mu nu_e e+
536  // 12 mu- -> nu_mu anti_nu_e e-
537  // 13 pi+ -> nu_mu mu+
538  // 14 pi- -> anti_nu_mu mu-
539 
540  if (theStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL)
541  {
542  G4int decay_code=0;
543  if (theStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName()=="Decay")
544  {
545  G4int nSecAtRest = fpSteppingManager->GetfN2ndariesAtRestDoIt();
546  G4int nSecAlong = fpSteppingManager->GetfN2ndariesAlongStepDoIt();
547  G4int nSecPost = fpSteppingManager->GetfN2ndariesPostStepDoIt();
548  G4int nSecTotal = nSecAtRest+nSecAlong+nSecPost;
549  G4TrackVector* secVec = fpSteppingManager->GetfSecondary();
550 
551  if (particleDefinition==G4PionPlus::PionPlusDefinition())
552  {
553  for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
554  {
555  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
556  decay_code=13;
557  }
558  }
559  if (particleDefinition==G4PionMinus::PionMinusDefinition())
560  {
561  for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
562  {
563  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
564  decay_code=14;
565  }
566  }
567  if (particleDefinition==G4KaonPlus::KaonPlusDefinition())
568  {
569  for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
570  {
571  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
572  {if (nSecTotal==2) decay_code=5;
573  if (nSecTotal==3) decay_code=7;}
574  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_e")
575  decay_code=6;
576  }
577  }
578  if (particleDefinition==G4KaonMinus::KaonMinusDefinition())
579  {
580  for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
581  {
582  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
583  {if (nSecTotal==2) decay_code=8;
584  if (nSecTotal==3) decay_code=10;}
585  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_e")
586  decay_code=9;
587  }
588  }
589  if (particleDefinition==G4KaonZeroLong::KaonZeroLongDefinition())
590  {
591  for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
592  {
593  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_e")
594  decay_code=1;
595  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_e")
596  decay_code=2;
597  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
598  decay_code=3;
599  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
600  decay_code=4;
601  }
602  }
603  if (particleDefinition==G4MuonPlus::MuonPlusDefinition())
604  {
605  for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
606  {
607  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="anti_nu_mu")
608  decay_code=11;
609  }
610  }
611  if (particleDefinition==G4MuonMinus::MuonMinusDefinition())
612  {
613  for (size_t partno=(*secVec).size()-nSecTotal;partno<(*secVec).size();partno++)
614  {
615  if ((*secVec)[partno]->GetDefinition()->GetParticleName()=="nu_mu")
616  decay_code=12;
617  }
618  }
619 
620  LBNETrackInformation* oldinfo=(LBNETrackInformation*)(theTrack->GetUserInformation());
621  if (oldinfo!=0)
622  {
623  oldinfo->SetDecayCode(decay_code);
624  theTrack->SetUserInformation(oldinfo);
625  }
626  else
627  {
629  newinfo->SetDecayCode(decay_code);
630  theTrack->SetUserInformation(newinfo);
631  }
632  }
633  }
634 }
635 
636 //----------------------------------------------------------------------
637 void LBNESteppingAction::CheckInHornEndPlane(const G4Step * theStep)
638 {
639 
640  G4Track * theTrack = theStep->GetTrack();
641 
642  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
643 
644  std::string preStepPointName =
645  theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
646  std::string postStepPointName =
647  theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
648 
649  if (((preStepPointName == G4String("TargetHallAndHorn1")) && (postStepPointName == G4String("Tunnel"))) ||
650  ((preStepPointName == G4String("Horn2Hall")) && (postStepPointName == G4String("Tunnel"))))
651  {
652  theTrack->SetTrackStatus(fStopAndKill);
653  G4VTrajectory* currTrajectory = 0;
654  currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());
655 
656  if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;
657 
658  LBNETrajectory* currLBNETrajectory = 0;
659  currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);
660 
661  if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;
662 
663  LBNETrajectory* newtrajectory = 0;
664  newtrajectory = new LBNETrajectory(*currLBNETrajectory);
665  if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;
666 
667  newtrajectory ->AppendStep(theStep);
668 
670  analysis->FillTrackingNtuple(*theTrack, newtrajectory);
671 
672  delete newtrajectory;
673  }
674 
675 
676 }
677 
678 //----------------------------------------------------------------------
680 {
681 
682  //G4Track * theTrack = theStep->GetTrack();
683 
684  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL || theStep->GetPostStepPoint()->GetPhysicalVolume()==NULL) return;
685 
686  G4int pdg = abs(theStep->GetTrack()->GetDefinition()->GetPDGEncoding());
687  if (pdg == 12 || pdg == 14 || pdg == 16) return;//Don't save neutrino info
688 
689  std::string preStepPointName =
690  theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
691  std::string postStepPointName =
692  theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
693 
694  if ((preStepPointName == G4String("MuonTrackingVol")) || (postStepPointName == G4String("MuonTrackingVol")))
695  {
696 
697  G4VTrajectory* currTrajectory = 0;
698  currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());
699 
700  if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;
701 
702  LBNETrajectory* currLBNETrajectory = 0;
703  currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);
704 
705  if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;
706 
708  analysis->FillAlcoveTrackingPlaneData(*theStep);
709 
710  }
711 
712 
713 }
714 
715 //----------------------------------------------------------------------
716 void LBNESteppingAction::CheckInTgtEndPlane(const G4Step * theStep)
717 {
718 
719  G4Track * theTrack = theStep->GetTrack();
720 
721  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
722 
723  std::string preStepPointName = theStep->GetPreStepPoint()->GetPhysicalVolume()->GetName();
724 
725  std::stringstream endName;
726  endName << "Horn1TopLevelDown";
727 
728  if (preStepPointName.find(endName.str()) != std::string::npos)
729  {
730  theTrack->SetTrackStatus(fStopAndKill);
731  G4VTrajectory* currTrajectory = 0;
732  currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());
733 
734  if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;
735 
736  LBNETrajectory* currLBNETrajectory = 0;
737  currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);
738 
739  if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;
740 
741  LBNETrajectory* newtrajectory = 0;
742  newtrajectory = new LBNETrajectory(*currLBNETrajectory);
743  if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;
744 
745  newtrajectory ->AppendStep(theStep);
746 
748  analysis->FillTrackingNtuple(*theTrack, newtrajectory);
749 
750  delete newtrajectory;
751  }
752 
753 
754 }
755 void LBNESteppingAction::OpenAscii(const char *fname) {
756 
757  std::cerr << " LBNESteppingAction::OpenAscii Start with filename " << std::string(fname) << std::endl;
758  fOutStudy.open(fname);
759  if (!fOutStudy.is_open()) {
760  std::string descr("Can not open output file ");
761  descr += std::string(fname);
762  G4Exception("LBNESteppingAction::OpenAscii", "I/O Problem ", FatalException, descr.c_str());
763  }
764  if (fStudyGeantinoMode.find("Absorb") != std::string::npos) {
765  fOutStudy <<
766  " id x y z ILDecayChan ILHorn1Neck ILHorn2Entr ILNCDecayChan ILNCHorn1Neck ILNCHorn2Entr ILAlHorn2Entr" << std::endl;
767  } else if (fStudyGeantinoMode.find("Propa") != std::string::npos) {
768  fOutStudy << " id x y z xo yo zo zPost step matPre matPost " << std::endl;
769  } else if (fStudyGeantinoMode.find("PropCO") != std::string::npos) {
770  fOutStudy << " id x y z xo yo zo step matPre matPost " << std::endl;
771  } else if (fStudyGeantinoMode.find("MagTilt") != std::string::npos) {
772  fOutStudy << " id x y z xp yp zp " << std::endl;
773  } else if (fStudyGeantinoMode.find("RZVoxels") != std::string::npos) {
774  fOutStudy << " ik r z nCr " << std::endl;
775  } else if (doStudyParticleThroughHorns) {
776  fOutStudy << " evt tr id sNum x y z px py pz w " << std::endl;
777  }
778  fOutStudy.flush();
779  std::cerr << " LBNESteppingAction::OpenAscii " << std::string(fname) << std::endl;
780 }
781 
782 void LBNESteppingAction::OpenNtuple(const char *fname) {
783 
784  std::cout << " LBNESteppingAction::OpenNtuple Start with filename " << std::string(fname) << std::endl;
785  steppingTupleFile = new TFile(fname,"RECREATE");
786 
787  if (steppingTupleFile->IsZombie()) {
788  std::string descr("Can not open output file ");
789  descr += std::string(fname);
790  G4Exception("LBNESteppingAction::OpenNtuple", "I/O Problem ", FatalException, descr.c_str());
791  }
792  makeSteppingTuple = true;
793  steppingTuple = new TTree("steppingTuple","steppingTuple");
794  steppingTuple->SetDirectory(steppingTupleFile);
795  steppingTuple->Branch("x",&steppingTuple_x,"x/D");
796  steppingTuple->Branch("y",&steppingTuple_y,"y/D");
797  steppingTuple->Branch("z",&steppingTuple_z,"z/D");
798  steppingTuple->Branch("density",&steppingTuple_density,"density/D");
799 
800  std::cout << " LBNESteppingAction::OpenNtuple " << std::string(fname) << std::endl;
801 }
802 
803 void LBNESteppingAction::StudyAbsorption(const G4Step * theStep) {
804 //
805 //make sure we are dealing with a geantino, or a mu, to include lengthening of step due to curling in
806 // B Field
807 //
808  G4Track * theTrack = theStep->GetTrack();
809  if (((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) && (
810  ((theTrack->GetParticleDefinition()->GetParticleName()).find("mu+") == std::string::npos ))) return;
811 
812 
813 
814 
815  G4StepPoint* prePtr = theStep->GetPreStepPoint();
816  if (prePtr == 0) return;
817 /*
818  if ( theTrack->GetNextVolume() == 0 ) {
819  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
820  for (size_t k=0; k!=3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
821  fOutStudy << " " << totalAbsDecayChan << " " << totalAbsHorn1Neck
822  << " " << totalAbsHorn2Entr;
823  fOutStudy << " " << waterAbsDecayChan << " " << waterAbsHorn1Neck
824  << " " << waterAbsHorn2Entr << " " << alumAbsHorn2Entr << std::endl;
825  return;
826  }
827 */
828  //
829  // I set the position of the geantino production vertex at Z=0.;
830  //
831 // G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
832  if (fEvtIdPrevious != pRunManager->GetCurrentEvent()->GetEventID() ) {
833 // std::cerr << " Evt id " <<
834 // pRunManager->GetCurrentEvent()->GetEventID() <<
835 // " Starting point, z = " << prePtr->GetPosition()[2] << std::endl;
836  totalAbsDecayChan= 0.;
839  waterAbsDecayChan= 0.;
842  alumAbsHorn2Entr=0.;
843  goneThroughHorn1Neck = false;
844  goneThroughHorn2Entr = false;
845  fEvtIdPrevious = pRunManager->GetCurrentEvent()->GetEventID();
846  return;
847  }
848 
849  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
850  const double ll = theStep->GetStepLength();
851  G4StepPoint* postPtr = theStep->GetPostStepPoint();
852 /*
853  if (postPtr == NULL) {
854  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
855  for (size_t k=0; k!=3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
856  fOutStudy << " " << totalAbsDecayChan << " " << totalAbsHorn1Neck
857  << " " << totalAbsHorn2Entr;
858  fOutStudy << " " << waterAbsDecayChan << " " << waterAbsHorn1Neck
859  << " " << waterAbsHorn2Entr << " " << alumAbsHorn2Entr << std::endl;
860  return;
861  }
862 */
863  if (postPtr == NULL) return;
864  G4VPhysicalVolume *physVol = postPtr->GetPhysicalVolume();
865  std::string vName(physVol->GetName());
866  G4Material *material = theTrack->GetMaterial();
867 
868  if(makeSteppingTuple) {
869  steppingTuple_x = postPtr->GetPosition()[0];
870  steppingTuple_y = postPtr->GetPosition()[1];
871  steppingTuple_z = postPtr->GetPosition()[2];
872  steppingTuple_density = material->GetDensity();
873  steppingTuple->Fill();
874  }
875 
876  if (pRunManager->GetCurrentEvent()->GetEventID() < -3) {
877  const double r = std::sqrt(postPtr->GetPosition()[0]*postPtr->GetPosition()[0] +
878  postPtr->GetPosition()[1]*postPtr->GetPosition()[1]);
879  const double t = r/std::abs(postPtr->GetPosition()[2] + 515.25); // ZOrigin = -515.25
880  // debugging only valid for Zorigin of -515.
881  std::cerr << " r = " << r << " theta " << t << " z = " << postPtr->GetPosition()[2] <<
882  " In " << vName << " material " << material->GetName()
883  << " InterLength " << material->GetNuclearInterLength() << std::endl;
884  }
885  if (postPtr->GetPosition()[2] > 500.) goneThroughHorn1Neck=true; // approximate...
886  if (postPtr->GetPosition()[2] > 6360.) goneThroughHorn2Entr=true; //truly approximate.
887  if (ll < 1.0e-10) return;
888 
889  totalAbsDecayChan += ll/material->GetNuclearInterLength();
890  if (vName.find("WaterLayer") != std::string::npos) waterAbsDecayChan += ll/material->GetNuclearInterLength();
891  if (!goneThroughHorn1Neck) {
892  totalAbsHorn1Neck += ll/material->GetNuclearInterLength();
893  if (vName.find("WaterLayer") != std::string::npos) waterAbsHorn1Neck += ll/material->GetNuclearInterLength();
894  }
895  if (!goneThroughHorn2Entr) {
896  totalAbsHorn2Entr += ll/material->GetNuclearInterLength();
897 // if (theTrack->GetTrackLength() < (6000.0*CLHEP::mm))
898 // std::cerr << " trackLength = " << theTrack->GetTrackLength() << " Z = " << postPtr->GetPosition()[2] <<
899 // " Abs L " << totalAbsHorn2Entr << std::endl;
900 
901  if (vName.find("WaterLayer") != std::string::npos) {
902  waterAbsHorn2Entr += ll/material->GetNuclearInterLength();
903  } else {
904  if ((vName.find("Horn1") != std::string::npos) && (material->GetName().find("Alumin") != std::string::npos))
905  alumAbsHorn2Entr += ll/material->GetNuclearInterLength();
906  }
907  }
908  if (vName.find("DecayPipe") != std::string::npos) {
909  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
910  for (size_t k=0; k!=3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
911  fOutStudy << " " << totalAbsDecayChan << " " << totalAbsHorn1Neck
912  << " " << totalAbsHorn2Entr;
913  fOutStudy << " " << waterAbsDecayChan << " " << waterAbsHorn1Neck
914  << " " << waterAbsHorn2Entr << " " << alumAbsHorn2Entr << std::endl;
915  theTrack->SetTrackStatus(fStopAndKill);
916  return;
917  }
918 
919 }
920 void LBNESteppingAction::StudyPropagation(const G4Step * theStep) {
921 //
922 //make sure we are dealing with a geantino...
923 //
924  G4Track * theTrack = theStep->GetTrack();
925  if ((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) return;
926  G4StepPoint* prePtr = theStep->GetPreStepPoint();
927  if (prePtr == 0) return;
928  G4StepPoint* postPtr = theStep->GetPostStepPoint();
929  if (postPtr == 0) return;
930  if (theStep->GetStepLength() < 0.1*CLHEP::mm) return;
931  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
932  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
933  std::string volNamePost(volPost->GetName());
934  std::string volNamePre(volPre->GetName());
935 
936  if ((volNamePost.find(fKeyVolumeForOutput.c_str()) != std::string::npos) ||
937  (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos)) {
938  std::cout << " LBNESteppingAction::StudyPropagation, critical volume " << std::endl
939  << ".. from " << volNamePre << " to " << volNamePost
940  << " detected at " << prePtr->GetPosition() << " going to " << postPtr->GetPosition() << std::endl;
941  const G4VTouchable *preHist = prePtr->GetTouchable();
942  const G4VTouchable *postHist = postPtr->GetTouchable();
943  const G4int nDepthPre = preHist->GetHistoryDepth();
944  const G4int nDepthPost = postHist->GetHistoryDepth();
945  std::cerr << " .... History depth for pre point " << nDepthPre << " Post " << nDepthPost << std::endl;
946  for (int k=0; k!= std::max(nDepthPre, nDepthPost); k++) {
947  if ((k <nDepthPre) && (k <nDepthPost))
948  std::cerr << " ............. Translation at depth, pre ... " << preHist->GetTranslation(k)
949  << " Post " << postHist->GetTranslation(k) << std::endl;
950  else if (k <nDepthPre)
951  std::cerr << " ..............Translation at depth, pre ... " << preHist->GetTranslation(k) << std::endl;
952  else if (k <nDepthPost)
953  std::cerr << " ..............Translation at depth, post ... " << postHist->GetTranslation(k) << std::endl;
954  }
955  std::cerr << " ------------------------------------------- " << std::endl;
956  }
957  if (volPre->GetMaterial()->GetName() == volPost->GetMaterial()->GetName()) return;
958  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
959  for (int k=0; k != 3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
960  for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
961  fOutStudy << " " << postPtr->GetPosition()[2];
962  fOutStudy << " " << theStep->GetStepLength();
963  fOutStudy << " " << volPre->GetMaterial()->GetName();
964  fOutStudy << " " << volPost->GetMaterial()->GetName();
965  fOutStudy << std::endl;
966  G4String vName(volPre->GetName());
967 // if (vName.find("DecayPipe") != std::string::npos) {
968 // theTrack->SetTrackStatus(fStopAndKill);
969 // }
970 }
971 void LBNESteppingAction::StudyCheckOverlap(const G4Step * theStep) {
972 //
973 //make sure we are dealing with a geantino...
974 //
975  G4Track * theTrack = theStep->GetTrack();
976  // Skip this cut for X-checking the muon-genatinos
977 // if ((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) return;
978  if( theTrack->GetNextVolume() == 0 ) {
979  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
980  std::cout << " Out of world with track length = " << theTrack->GetTrackLength() << std::endl;
981  return;
982  }
983  G4StepPoint* prePtr = theStep->GetPreStepPoint();
984  if (prePtr == 0) return;
985  G4StepPoint* postPtr = theStep->GetPostStepPoint();
986  if (postPtr == 0) {
987  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
988  std::cout << " No Post Point, Last step at Z = " << prePtr->GetPosition()[2] << " r " <<
989  std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
990  prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl;
991 
992  return;
993  }
994  if (postPtr->GetPhysicalVolume() == 0) {
995  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
996  std::cout << " No Post Point Volume Last step at Z = " << prePtr->GetPosition()[2] << " r " <<
997  std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
998  prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl;
999 
1000  return;
1001  }
1002  if (prePtr->GetPhysicalVolume() == 0) return;
1003  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
1004  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
1005  if (volPre == 0) return;
1006  if (volPost == 0) return;
1007  std::string volNamePost(volPost->GetName());
1008  std::string volNamePre(volPre->GetName());
1009  std::string matNamePre = volPre->GetMaterial()->GetName();
1010  std::string matNamePost = volPost->GetMaterial()->GetName();
1011  //
1012  if (pRunManager->GetCurrentEvent()->GetEventID() < 3) {
1013  std::cout << " at Z = " << prePtr->GetPosition()[2] << " r " <<
1014  std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
1015  prePtr->GetPosition()[1]*prePtr->GetPosition()[1])
1016  << ", " << volNamePre << " to " << postPtr->GetPosition() << ", "
1017  << volNamePost << " mat from " << matNamePre << " to " << matNamePost;
1018  if (std::abs(theTrack->GetParticleDefinition()->GetPDGEncoding()) == 13) {
1019  std::cout << " Ek " << theTrack->GetKineticEnergy()/CLHEP::GeV << std::endl;
1020  } else { std::cout << std::endl; }
1021  /* These debugging statement may cause a fatal exception, because A and/or Z not defined for mixtures..
1022  if (matNamePost.find("Targ") != std::string::npos) {
1023  std::cout << " Target material found, Z " << volPost->GetMaterial()->GetZ() << " A "
1024  << volPost->GetMaterial()->GetA() << std::endl;
1025  }
1026  if (matNamePost.find("DecayPipe") != std::string::npos) {
1027  std::cout << " Decay Pipe material found, Z " << volPost->GetMaterial()->GetZ() << " A "
1028  << volPost->GetMaterial()->GetA()/mole << " density " << volPost->GetMaterial()->GetDensity()/(g/cm3) << std::endl;
1029  }
1030  */
1031  }
1032 // if (((volNamePost.find(fKeyVolumeForOutput.c_str()) != std::string::npos) ||
1033 // (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos)) &&
1034 // ( (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos) ||
1035 // (volNamePre.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos))) {
1036  if ( (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos) &&
1037  (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos)) {
1038  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
1039  for (int k=0; k != 3; k++) fOutStudy << " " << prePtr->GetPosition()[k];
1040  for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
1041  fOutStudy << " " << theStep->GetStepLength();
1042  fOutStudy << " " << volPre->GetMaterial()->GetName();
1043  fOutStudy << " " << volPost->GetMaterial()->GetName();
1044  fOutStudy << std::endl;
1045  }
1046 }
1047 
1048 void LBNESteppingAction::StudyCheckMagneticTilts(const G4Step * theStep) {
1049 //
1050 //make sure we are dealing with a geantino...
1051 //
1052  G4Track * theTrack = theStep->GetTrack();
1053  if ((theTrack->GetParticleDefinition()->GetParticleName()).find("mu") == std::string::npos) return;
1054  if( theTrack->GetNextVolume() == 0 ) {
1055  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
1056  std::cout << " Out of world with track length = " << theTrack->GetTrackLength() << std::endl;
1057  return;
1058  }
1059  G4StepPoint* prePtr = theStep->GetPreStepPoint();
1060  if (prePtr == 0) return;
1061  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
1062  std::cout << " at Z = " << prePtr->GetPosition()[2] << " r " <<
1063  std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
1064  prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl;
1065  G4StepPoint* postPtr = theStep->GetPostStepPoint();
1066  if (postPtr == 0) {
1067  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
1068  std::cout << " No Post Point, Last step at Z = " << prePtr->GetPosition()[2] << " r " <<
1069  std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
1070  prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl;
1071 
1072  return;
1073  }
1074  if (postPtr->GetPhysicalVolume() == 0) {
1075  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
1076  std::cout << " No Post Point Volume Last step at Z = " << prePtr->GetPosition()[2] << " r " <<
1077  std::sqrt(prePtr->GetPosition()[0]*prePtr->GetPosition()[0] +
1078  prePtr->GetPosition()[1]*prePtr->GetPosition()[1]) << std::endl;
1079 
1080  return;
1081  }
1082  if (prePtr->GetPhysicalVolume() == 0) return;
1083  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
1084  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
1085  if (volPre == 0) return;
1086  if (volPost == 0) return;
1087  std::string volNamePost(volPost->GetName());
1088  std::string volNamePre(volPre->GetName());
1089 // if (((volNamePost.find(fKeyVolumeForOutput.c_str()) != std::string::npos) ||
1090 // (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos)) &&
1091 // ( (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos) ||
1092 // (volNamePre.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos))) {
1093  if ( (volNamePre.find(fKeyVolumeForOutput.c_str()) != std::string::npos) &&
1094  (volNamePost.find(fKeyVolumeForOutputTo.c_str()) != std::string::npos)) {
1095  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
1096  for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
1097  for (int k=0; k != 3; k++) fOutStudy << " " << theTrack->GetMomentumDirection()[k];
1098  fOutStudy << std::endl;
1099  }
1100 }
1101 
1103 //
1104 //make sure we are dealing with a geantino...
1105 //
1106  G4Track * theTrack = theStep->GetTrack();
1107  if ((theTrack->GetParticleDefinition()->GetParticleName()).find("geantino") == std::string::npos) return;
1108  if( theTrack->GetNextVolume() == 0 ) {
1109  if (pRunManager->GetCurrentEvent()->GetEventID() < 3)
1110  std::cout << " Out of world with track length = " << theTrack->GetTrackLength() << std::endl;
1111  return;
1112  }
1113  G4StepPoint* prePtr = theStep->GetPreStepPoint();
1114  if (prePtr == 0) return;
1115  G4StepPoint* postPtr = theStep->GetPostStepPoint();
1116  if (postPtr == 0) return;
1117  if (postPtr->GetPhysicalVolume() == 0) return;
1118  if (prePtr->GetPhysicalVolume() == 0) return;
1119  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
1120  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
1121  if (volPre == 0) return;
1122  if (volPost == 0) return;
1123  double zPost = postPtr->GetPosition()[2];
1124  double rPost = std::sqrt(postPtr->GetPosition()[0]*postPtr->GetPosition()[0] +
1125  postPtr->GetPosition()[1]*postPtr->GetPosition()[1]);
1126  zPost = std::min(zPost, fGenerateVoxelsZMax);
1127  zPost = std::max(zPost, fGenerateVoxelsZScale);
1128  rPost = std::min(rPost, fGenerateVoxelsRMax );
1129  rPost = std::max(rPost, 1.);
1130  double lnR = std::log(rPost);
1131  int iz = static_cast<int>((zPost + fGenerateVoxelsZOffset)/fGenerateVoxelsZScale);
1132  int ir = static_cast<int>(lnR*fGenerateVoxelsRScale);
1133  int ik = iz*fGenerateVoxelsIRKMax + ir;
1135  if(itMRZ == fRZVoxelsData.end())
1136  fRZVoxelsData.insert(std::pair<int, int>(ik, 1));
1137  else
1138  itMRZ->second += 1;
1139 }
1140 
1142 
1143  int idEvt = pRunManager->GetCurrentEvent()->GetEventID();
1144  if (fEvtIdPrevious != idEvt) {
1145  if (fOutStepStudy.is_open()) fOutStepStudy.close();
1146  std::ostringstream fNameStrStr; fNameStrStr << "./StepSudiesMagn_Evt" << idEvt << ".txt";
1147  std::string fNameStr(fNameStrStr.str());
1148  fOutStepStudy.open(fNameStr.c_str());
1149  fOutStepStudy << " x y z xp yp dpt bphi bz radIC radOC factR vName " << std::endl;
1150  fEvtIdPrevious = idEvt;
1151  }
1152  G4StepPoint* prePtr = theStep->GetPreStepPoint();
1153  if (theStep->GetPreStepPoint()->GetPhysicalVolume() == 0) return;
1154  if (prePtr == 0) return;
1155  G4Track * theTrack = theStep->GetTrack();
1156  if ( theTrack->GetNextVolume() == 0 ) {
1157  if (fOutStepStudy.is_open()) fOutStepStudy.close();
1158  return;
1159  }
1160  G4StepPoint* postPtr = theStep->GetPostStepPoint();
1161  if (postPtr == 0) {
1162  if (fOutStepStudy.is_open()) fOutStepStudy.close();
1163  return;
1164  }
1165  for (size_t k=0; k !=3; k++) fOutStepStudy << " " << postPtr->GetPosition()[k];
1166  G4ThreeVector postVec = postPtr->GetMomentum();
1167  G4ThreeVector preVec = prePtr->GetMomentum();
1168  fOutStepStudy << " " << postVec[0]/postVec[2];
1169  fOutStepStudy << " " << postVec[1]/postVec[2];
1170  const double ptPre = std::sqrt(preVec[0]*preVec[0] + preVec[1]*preVec[1])/CLHEP::GeV;
1171  const double ptPost = std::sqrt(postVec[0]*postVec[0] + postVec[1]*postVec[1])/CLHEP::GeV;
1172  fOutStepStudy << " " << ptPost - ptPre;
1173  // Use our our utility to navigate to find the right field
1174  const G4Field *aField = 0;
1175  const LBNEMagneticFieldHorn *myField=0;
1176  const LBNFMagneticFieldPolyconeHorn *myFieldPH=0;
1177  double effInnerRad = 0; double effOuterRad = 0; double effSkinD = 0.;
1179  if ((aPlacementHandler->GetNumberOfHornsPolycone() == 0) &&
1180  (postPtr->GetPosition()[2] > -500.) && (postPtr->GetPosition()[2] < 5000.)) {
1181  const LBNEVolumePlacementData *pl =
1182  aPlacementHandler->Find("FieldHorn1", "Horn1PolyM1", "DetectorConstruction");
1183  aField = pl->fCurrent->GetFieldManager()->GetDetectorField();
1184  if (aField != 0) {
1185  myField = dynamic_cast<const LBNEMagneticFieldHorn *> (aField);
1186  if (myField != 0) {
1187  effInnerRad = myField->GetEffectiveInnerRad();
1188  effOuterRad = myField->GetEffectiveOuterRad();
1189  effSkinD = myField->GetEffectiveSkinDepthFact();
1190  }
1191  }
1192  } else if ((aPlacementHandler->GetNumberOfHornsPolycone() == 0) && (postPtr->GetPosition()[2] < 13000.)) {
1193  const LBNEVolumePlacementData *pl =
1194  aPlacementHandler->Find("FieldHorn2", "Horn2Hall", "DetectorConstruction");
1195  aField = pl->fCurrent->GetFieldManager()->GetDetectorField();
1196  myField = dynamic_cast<const LBNEMagneticFieldHorn *> (aField);
1197  if (myField != 0) {
1198  effInnerRad =myField->GetEffectiveInnerRad();
1199  effOuterRad = myField->GetEffectiveOuterRad();
1200  effSkinD = myField->GetEffectiveSkinDepthFact();
1201  }
1202  } else if ((aPlacementHandler->GetNumberOfHornsPolycone() != 0) && (postPtr->GetPosition()[2] < 22000.)) {
1203  G4FieldManager *aFieldMgr = postPtr->GetPhysicalVolume()->GetLogicalVolume()->GetFieldManager();
1204  if (aFieldMgr != 0) {
1205  aField = aFieldMgr->GetDetectorField();
1206  myFieldPH = dynamic_cast<const LBNFMagneticFieldPolyconeHorn *> (aField);
1207  if (myFieldPH != 0) {
1208  effInnerRad =myFieldPH->GetEffectiveInnerRad();
1209  effOuterRad = myFieldPH->GetEffectiveOuterRad();
1210  effSkinD = myFieldPH->GetEffectiveSkinDepthFact();
1211  }
1212  }
1213  }
1214 
1215  if (aField == 0) {
1216  fOutStepStudy << " 0. 0. 0. 0. 0. ";
1217  } else {
1218  double pos[3]; for (size_t k=0; k!= 3; ++k) pos[k] = postPtr->GetPosition()[k];
1219  double bf[3];
1220  if (myField != 0) myField->GetFieldValue(pos, &bf[0]);
1221  if (myFieldPH != 0) {
1222  myFieldPH->GetFieldValue(pos, &bf[0]); // We repeat the checks on effective radii..
1223  effInnerRad =myFieldPH->GetEffectiveInnerRad();
1224  effOuterRad = myFieldPH->GetEffectiveOuterRad();
1225  effSkinD = myFieldPH->GetEffectiveSkinDepthFact();
1226  }
1227  if((myField == 0) && (myFieldPH == 0)) aField->GetFieldValue(pos, &bf[0]);
1228  fOutStepStudy << " " << std::sqrt(bf[0]*bf[0] + bf[1]*bf[1])/CLHEP::tesla;
1229  fOutStepStudy << " " << bf[2]/CLHEP::tesla;
1230  fOutStepStudy << " " << effInnerRad << " " << effOuterRad << " " << effSkinD;
1231  }
1232  fOutStepStudy << " " << postPtr->GetPhysicalVolume()->GetLogicalVolume()->GetName() << std::endl;
1233 }
1234 
1236 
1237  G4Track * theTrack = theStep->GetTrack();
1238  if (theTrack == 0) return;
1239  //
1240  // Simplified version for animation with muons...
1241  //
1242  /*
1243  if (theTrack->GetPosition()[2] > 15000.) return;
1244  G4ThreeVector direction = theTrack->GetMomentumDirection();
1245  const double eTot = theTrack->GetTotalEnergy();
1246  const double mass= theTrack->GetParticleDefinition()->GetPDGMass();
1247  const double pTot = std::sqrt(eTot*eTot - mass*mass);
1248  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
1249  fOutStudy << " " << theTrack->GetParticleDefinition()->GetPDGEncoding();
1250  const int stepNum = theTrack->GetCurrentStepNumber();
1251  fOutStudy << " " << stepNum;
1252  if (stepNum < 2) {
1253  for (int k=0; k != 3; k++) fOutStudy << " " << theTrack->GetVertexPosition()[k];
1254  for (int k=0; k != 3; k++) fOutStudy << " " << 0.001*pTot*theTrack->GetVertexMomentumDirection()[k];
1255  fOutStudy << std::endl;
1256  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
1257  fOutStudy << " " << theTrack->GetParticleDefinition()->GetPDGEncoding();
1258  fOutStudy << " " << stepNum;
1259  }
1260  for (int k=0; k != 3; k++) fOutStudy << " " << theTrack->GetPosition()[k];
1261  for (int k=0; k != 3; k++) fOutStudy << " " << 0.001*pTot*direction[k];
1262  fOutStudy << std::endl;
1263  return;
1264  */
1265  G4StepPoint* prePtr = theStep->GetPreStepPoint();
1266  G4StepPoint* postPtr = theStep->GetPostStepPoint();
1267  if (prePtr->GetPhysicalVolume() == 0) return;
1268  if (postPtr == 0) return;
1269  if (postPtr->GetPhysicalVolume() == 0) return;
1270  //
1271  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
1272  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
1273  std::string volNamePost(volPost->GetName());
1274  std::string volNamePre(volPre->GetName());
1275  bool ofInterest = false;
1276  // Entering the corridor
1277  if ((volNamePost == std::string("Tunnel")) &&
1278  (volNamePre.find("Horn1") != std::string::npos )) ofInterest=true;
1279  if ((volNamePre.find("Horn2") != std::string::npos) &&
1280  (volNamePost == std::string("Tunnel"))) ofInterest=true;
1281  if ((volNamePre == std::string("Tunnel")) &&
1282  (volNamePost == std::string("DecayPipeHall"))) ofInterest=true;
1283  if (!ofInterest) return;
1284  const double eTot = theTrack->GetTotalEnergy();
1285  const double mass= theTrack->GetParticleDefinition()->GetPDGMass();
1286  const double pTot = std::sqrt(eTot*eTot - mass*mass);
1287  if (pTot < 0.1*CLHEP::GeV) return;
1288  if (pTot > 50.*CLHEP::GeV) return;
1289  int pdgInfo = theTrack->GetParticleDefinition()->GetPDGEncoding();
1290  // check only the muons for now.
1291  if (std::abs(pdgInfo) != 211) return;
1292  fOutStudy << " " << pRunManager->GetCurrentEvent()->GetEventID();
1293  fOutStudy << " " << theTrack->GetTrackID();
1294  const int stepNum = theTrack->GetCurrentStepNumber();
1295  fOutStudy << " " << pdgInfo << " " << stepNum;
1296  for (int k=0; k != 3; k++) fOutStudy << " " << postPtr->GetPosition()[k];
1297  G4ThreeVector direction = theTrack->GetMomentumDirection();
1298  for (int k=0; k != 3; k++) fOutStudy << " " << 0.001*pTot*direction[k];
1299  LBNETrackInformation* info = reinterpret_cast<LBNETrackInformation*>(theTrack->GetUserInformation());
1300  fOutStudy << " " << info->GetNImpWt(); // Importance weight
1301  fOutStudy << std::endl;
1302 
1303 }
1304 void LBNESteppingAction::StudyPionsThroughHorn2 (const G4Step* theStep) {
1305  const G4Track * theTrack = theStep->GetTrack();
1306  if (theTrack == 0) return;
1307  int pID = theTrack->GetParticleDefinition()->GetPDGEncoding();
1308  if (std::abs(pID) != 211) return;
1309  G4StepPoint* prePtr = theStep->GetPreStepPoint();
1310  if (prePtr == 0) return;
1311  if (prePtr->GetPhysicalVolume() == 0) return;
1312  G4StepPoint* postPtr = theStep->GetPostStepPoint();
1313  if (postPtr == 0) return;
1314  if (postPtr->GetPhysicalVolume() == 0) return;
1315  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
1316  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
1317  std::string volNamePost(volPost->GetName());
1318  std::string volNamePre(volPre->GetName());
1319  if ((volNamePre == std::string("Horn2Part3Ring")) &&
1320  (volNamePost.find("Horn2") != std::string::npos)) {
1322  int sign = (pID > 0) ? 1 : -1;
1323  ptr->RingHangerPion(theTrack->GetTrackID(), sign, theTrack->GetMomentum(), postPtr->GetPosition());
1324  }
1325 }
1326 
1328 {
1329  // std::cout << "I am here" << std::endl;
1330 
1331  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
1332 
1333  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
1334  theStep->GetPostStepPoint()->GetPhysicalVolume()){
1335  G4LogicalVolume *preStepVolume =
1336  theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1337  G4LogicalVolume *postStepVolume =
1338  theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1339  if(postStepVolume == TrkPlnLogical){
1340  if(preStepVolume != TrkPlnLogical){
1341  G4VTrajectory* currTrajectory = 0;
1342  currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());
1343 
1344  if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;
1345 
1346  LBNETrajectory* currLBNETrajectory = 0;
1347  currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);
1348  if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;
1349 
1350  LBNETrajectory* newtrajectory = 0;
1351  newtrajectory = new LBNETrajectory(*currLBNETrajectory);
1352  if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;
1353 
1354  newtrajectory ->AppendStep(theStep);
1355 
1357  //analysis->FillTrackingNtuple(*theTrack, newtrajectory);
1358  analysis->FillTrackingPlaneData(*theStep);
1359 
1360  delete newtrajectory;
1361  }
1362  }
1363  }
1364 }
1365 
1366 
1367 void LBNESteppingAction::CheckInTargetOutput(const G4Step *theStep)
1368 {
1369 
1370  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
1371 std::string preStepPointName =
1372 theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName();
1373 //std::string postStepPointName =
1374 //theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
1375 std::stringstream endName, posvolname;
1376 endName<<"Target";
1377 posvolname<<"TargetHallandHorn1";
1378 if(preStepPointName.find(endName.str())!=std::string::npos)
1379 //if(preStepPointName==postStepPointName) return;
1380 //if(preStepPointName!=postStepPointName)
1381 {
1382 //theTrack->SetTrackStatus(fStopAndKill);
1383 G4VTrajectory* currTrajectory = 0;
1384  currTrajectory = (G4EventManager::GetEventManager()->GetTrackingManager()->GimmeTrajectory());
1385 
1386  if(!currTrajectory) std::cout << "Current trajectory is invalid" << std::endl;
1387 
1388  LBNETrajectory* currLBNETrajectory = 0;
1389  currLBNETrajectory = dynamic_cast<LBNETrajectory*>(currTrajectory);
1390 
1391  if(!currLBNETrajectory) std::cout << "Cast failed"<< std::endl;
1392 
1393  LBNETrajectory* newtrajectory = 0;
1394  newtrajectory = new LBNETrajectory(*currLBNETrajectory);
1395  if(!newtrajectory) std::cout << "new trajectory is invalid" << std::endl;
1396 
1397  newtrajectory ->AppendStep(theStep);
1398 
1400  analysis->FillTargetOutputData(*theStep);
1401 
1402  delete newtrajectory;
1403 }
1404 }
1405 
1406 
1407 // Amit Bashyal
1408 // Does this still work ? Paul Lebrun, May 26 2016.
1410 {
1411  // std::cout << "I am here" << std::endl;
1412 // G4Track * theTrack = theStep->GetTrack();
1413 
1414  if ((theStep->GetPreStepPoint() == 0) || (theStep->GetPostStepPoint() == 0)) return;
1415  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
1416 
1417  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
1418  theStep->GetPostStepPoint()->GetPhysicalVolume()){
1419  G4LogicalVolume *preStepVolume =
1420  theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1421  G4LogicalVolume *postStepVolume =
1422  theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1423  if(postStepVolume == TrkPlnH1Logical){
1424  if(preStepVolume != TrkPlnH1Logical){
1426  // analysis->FillTrackingPlaneH1Data(*theTrack, newtrajectory);
1427  // std::cout<<"Filling H1Trackingplane Ntuples"<<std::endl;
1428  analysis->FillTrackingPlaneH1Data(*theStep);
1429  }
1430  }
1431  }
1432 
1433 }
1434 //Amit Bashyal
1436 {
1437  // std::cout << "I am here" << std::endl;
1438  if ((theStep->GetPreStepPoint() == 0) || (theStep->GetPostStepPoint() == 0)) return;
1439  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
1440 
1441  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
1442  theStep->GetPostStepPoint()->GetPhysicalVolume()){
1443  G4LogicalVolume *preStepVolume =
1444  theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1445  G4LogicalVolume *postStepVolume =
1446  theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1447  if(postStepVolume == TrkPlnH2Logical){
1448  if(preStepVolume != TrkPlnH2Logical){
1450  analysis->FillTrackingPlaneH2Data(*theStep);
1451  }
1452  }
1453  }
1454 }
1455 //Amit Bashyal
1457 {
1458  // std::cout << "I am here" << std::endl;
1459  if(theStep->GetPreStepPoint()->GetPhysicalVolume() == NULL) return;
1460 
1461  if(theStep->GetPreStepPoint()->GetPhysicalVolume() &&
1462  theStep->GetPostStepPoint()->GetPhysicalVolume()){
1463  G4LogicalVolume *preStepVolume =
1464  theStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1465  G4LogicalVolume *postStepVolume =
1466  theStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume();
1467  if(postStepVolume == TrkPlnDPLogical){
1468  if(preStepVolume != TrkPlnDPLogical){
1470  analysis->FillTrackingPlaneDPData(*theStep);
1471  }
1472  }
1473  }
1474 }
1477  // Set a file name that can be unique.. Particularly, tailored for the use on FermiGrid
1478  // But not mandated.
1479  int clusterNum = 0;
1480  int procNum =0;
1481  const char *clusterEnv = getenv("CLUSTER");
1482  if (clusterEnv != 0) clusterNum = atoi(clusterEnv);
1483  const char *procEnv = getenv("PROCESS");
1484  if (procEnv != 0) procNum = atoi(procEnv);
1485  std::string fName("./hadronFluxFromTargetASCII.txt");
1486  if ((clusterNum != 0) || (procNum != 0)) {
1487 // const char *userEnv = getenv("USER"); // assume it always works
1488 // std::string aUserStr(userEnv);
1489  std::ostringstream fNStrStr;
1490  fNStrStr << "./hadronFluxFromTargetASCII_"
1491  << clusterNum << "_" << procNum << ".txt";
1492  fName = fNStrStr.str();
1493  }
1495  fStrHadronFluxFromTargetASCII.open(fName.c_str());
1496  fStrHadronFluxFromTargetASCII << " evt tr id x y z px py pz w " << std::endl;
1497  std::cerr << " LBNESteppingAction::InitiateHadronFluxFromTargetASCII, Opening file " <<
1498  fName << std::endl;
1500  std::cerr << " O.K., ready for HadronFluxFromTargetASCII.... in... " << this << std::endl;
1501 }
1502 void LBNESteppingAction::FillHadronFluxFromTargetASCII(const G4Step *theStep) const {
1503 
1504  G4Track * theTrack = theStep->GetTrack();
1505  if (theTrack->GetPosition()[2] > 3500.) return; // Hopefully, the target is not longer..
1506 // std::cerr << " LBNESteppingAction::FillHadronFluxFromTargetASCII At z " << theTrack->GetPosition()[2] << std::endl;
1507  G4StepPoint* prePtr = theStep->GetPreStepPoint();
1508  if (prePtr == 0) return;
1509  G4StepPoint* postPtr = theStep->GetPostStepPoint();
1510  if (postPtr == 0) return;
1511  if (postPtr->GetPhysicalVolume() == 0) return;
1512  if (prePtr->GetPhysicalVolume() == 0) return;
1513  G4LogicalVolume *volPost = postPtr->GetPhysicalVolume()->GetLogicalVolume();
1514  G4LogicalVolume *volPre = prePtr->GetPhysicalVolume()->GetLogicalVolume();
1515  std::string volNamePost(volPost->GetName());
1516  std::string volNamePre(volPre->GetName());
1517 // std::cerr << " LBNESteppingAction::FillHadronFluxFromTargetASCII volNamePre "
1518 // << volNamePre << " Post " << volNamePost << std::endl;
1519 // if (((volNamePost.find("TargetNoSplitHelium") != std::string::npos) ||
1520 // (volNamePost.find("TargetNoSplitHeCont") != std::string::npos)) &&
1521 // (volNamePre.find("TargetUsptreamSimpleCyl") != std::string::npos)) {
1522  if ((volNamePre.find("TargetNoSplitHelium") != std::string::npos) &&
1523  (volNamePost.find("TargetNoSplitHeContainer") != std::string::npos)) {
1524  double pTotSq = 0; for (int k=0; k != 3; k++) pTotSq += theTrack->GetMomentum()[k]*theTrack->GetMomentum()[k];
1525  if (pTotSq < 2500) return;
1526  G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
1527  fStrHadronFluxFromTargetASCII << " " << evtno << " " << theTrack->GetTrackID() << " "
1528  << theTrack->GetParticleDefinition()->GetPDGEncoding();
1529 
1530  for (int k=0; k != 3; k++) fStrHadronFluxFromTargetASCII << " " << postPtr->GetPosition()[k];
1531  for (int k=0; k != 3; k++) fStrHadronFluxFromTargetASCII << " " << (theTrack->GetMomentum()[k])/CLHEP::GeV;
1532  LBNETrackInformation* info = (LBNETrackInformation*)(theTrack->GetUserInformation());
1533  fStrHadronFluxFromTargetASCII << " " << info->GetNImpWt();
1535  //
1536  // To speed up the detection of missing pion (with respect to FLUKA)
1537  //
1538 // theTrack->SetTrackStatus(fStopAndKill);
1539  }
1540 }
1541 
1542 
1544 // June 16 2015: add the capability to generate
1545  if (!fOutMuonSculptedAbsorber.is_open()){
1546  //fOutMuonSculptedAbsorber.open("./MuonFluxAtSculptedAbsorber.txt");
1547  G4cout <<"Opening muon file!!!!! "<< fOutMuonSculptedAbsorberFilename<<G4endl;
1549  fOutMuonSculptedAbsorber << " evt tr id slice x y z px py pz ek w g2 " << std::endl;
1550  }
1551  const G4Track* aTrack = aStep->GetTrack();
1552  int pdg = aTrack->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
1553  if (std::abs(pdg) != 13) return;
1554  const G4StepPoint *prePt = aStep->GetPreStepPoint();
1555  if (prePt == 0) return;
1556  if (prePt->GetPhysicalVolume() == 0) return;
1557  if (prePt->GetPosition()[2] < 200000.) return; // Rough cut!.
1558  const G4StepPoint *postPt = aStep->GetPostStepPoint();
1559  if (postPt->GetPhysicalVolume() == 0) return;
1560  G4LogicalVolume *volPost = postPt->GetPhysicalVolume()->GetLogicalVolume();
1561  G4LogicalVolume *volPre = prePt->GetPhysicalVolume()->GetLogicalVolume();
1562  std::string volNamePost(volPost->GetName());
1563  std::string volNamePre(volPre->GetName());
1564  //
1565  // debugging ...
1566  //
1567 // if ((volNamePre.find("HadronAbsorberSculpBackWall") == 0) || (volNamePre.find("HadronAbsorberSculpEndIronAC") == 0)) {
1568 // std::cerr << " Got the end of HA, vol name Pre " << volNamePre << " .. And quit " << std::endl; exit(2);
1569 // }
1570  int aSlice=-100;
1571  bool ofInterestLBNE = ((volNamePre == std::string("Tunnel")) && (volNamePost == std::string("HadronAbsorberTop"))) ||
1572  ((volNamePre == std::string("AH_Muon_alkair")) && (volNamePost == std::string("detLogical")));
1573 
1574  bool ofInterestLBNF = ((volNamePre == std::string("Tunnel")) && (volNamePost == std::string("HadronAbsorberSculpTop"))) ||
1575  ((volNamePost == std::string("Tunnel")) && (volNamePre == std::string("HadronAbsorberSculpTop"))) ||
1576  (volNamePre.find("HadronAbsorberSculp") != std::string::npos);
1577  if ((!ofInterestLBNF) && (!ofInterestLBNE)) return;
1578  if (volNamePre == std::string("Tunnel")) {
1579  aSlice = -2;
1581  }
1582  if (ofInterestLBNE) {
1583  if (volNamePre == std::string("AH_Muon_alkair")) aSlice = 30;
1584  } else {
1585  if ((volNamePre == std::string("HadronAbsorberSculpAibufferDiag")) && (volNamePost == std::string("HadronAbsorberSculpTop"))) {
1586  aSlice = -1;
1587  } else if ((volNamePre == std::string("HadronAbsorberSculpSpoilerAir")) && (volNamePost == std::string("HadronAbsorberSculpTop"))) {
1588  aSlice = 0;
1589  } else if (volNamePost == std::string("HadronAbsorberSculpTop")) {
1590  size_t iPosU = volNamePre.find("-");
1591  if (iPosU != std::string::npos) {
1592  iPosU++;
1593  std::string tSliceStr = volNamePre.substr(iPosU, std::string::npos);
1594  aSlice = atoi(tSliceStr.c_str());
1595  }
1596  if (volNamePre == std::string("HadronAbsorberSculpBackWall")) aSlice = 30;
1597  if (volNamePre == std::string("HadronAbsorberSculpEndIronAC")) aSlice = 29;
1598  }
1599  }
1600  if (aSlice < -5) return;
1601  G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
1602  fOutMuonSculptedAbsorber << " " << evtno << " " << aTrack->GetTrackID() << " " << pdg << " " << aSlice;
1603  for (size_t k=0; k!=3; k++) fOutMuonSculptedAbsorber << " " << prePt->GetPosition()[k];
1604  for (size_t k=0; k!=3; k++) fOutMuonSculptedAbsorber << " " << aTrack->GetMomentum()[k]/CLHEP::GeV;
1605  LBNETrackInformation* info = (LBNETrackInformation*)(aTrack->GetUserInformation());
1606  fOutMuonSculptedAbsorber << " " << aTrack->GetDynamicParticle()->GetKineticEnergy()/CLHEP::GeV << " "
1607  << info->GetNImpWt();
1609  else fOutMuonSculptedAbsorber << " 0";
1611 }
1613 
1614  size_t kSelOut = 99;
1615  const G4StepPoint *prePt = aStep->GetPreStepPoint();
1616  if (prePt == 0) return;
1617  if (prePt->GetPhysicalVolume() == 0) return;
1618  const G4StepPoint *postPt = aStep->GetPostStepPoint();
1619  if (postPt == 0) return;
1620  if (postPt->GetPhysicalVolume() == 0) return;
1621  G4LogicalVolume *volPost = postPt->GetPhysicalVolume()->GetLogicalVolume();
1622  G4LogicalVolume *volPre = prePt->GetPhysicalVolume()->GetLogicalVolume();
1623  std::string volNamePost(volPost->GetName());
1624  std::string volNamePre(volPre->GetName());
1625  if ((volNamePre == std::string("TargetNoSplitHeContainer")) &&
1626  (volNamePost == std::string("TargetNoSplitM1"))) kSelOut = 0;
1627  else if ((volNamePre == std::string("Horn1PolyM1")) &&
1628  (volNamePost == std::string("TargetHallAndHorn1"))) kSelOut = 1;
1629  else if ((volNamePre == std::string("LBNFSimpleHorn2Container")) &&
1630  (volNamePost == std::string("Tunnel"))) kSelOut = 2;
1631  else if ((volNamePre == std::string("LBNFSimpleHorn3Container")) &&
1632  (volNamePost == std::string("Tunnel"))) kSelOut = 3;
1633 
1634  if (kSelOut == 99) return;
1635  G4Track* aTrack = aStep->GetTrack();
1636  int pdg = aTrack->GetDynamicParticle()->GetDefinition()->GetPDGEncoding();
1637  bool isMeson = (std::abs(pdg) == 211) || (std::abs(pdg) == 321);
1638  if (!isMeson) return;
1639  G4ThreeVector xyz = postPt->GetPosition();
1640  G4ThreeVector pMom = postPt->GetMomentum();
1641  G4int evtno = pRunManager->GetCurrentEvent()->GetEventID();
1642  (*fOutPtrsForMarsCmpApr2017[kSelOut]) << " " << evtno << " " << aTrack->GetTrackID() << " " << pdg << " "
1643  << xyz[0]/10. << " " << xyz[1]/10. << " " << xyz[2]/10. << " "
1644  << pMom[0]/1000. << " " << pMom[1]/1000. << " " << pMom[2]/1000. << " "
1645  << postPt->GetKineticEnergy()/1000. << std::endl;
1646  if (kSelOut == 3) aTrack->SetTrackStatus(fStopAndKill);
1647 }
void FillTrackingPlaneH1Data(const G4Step &aStep)
intermediate_table::iterator iterator
void fill(const G4Step *aStep)
Definition: LBNFDeDxMap.cc:238
LBNESteppingActionMessenger * pMessenger
void StudyCheckOverlap(const G4Step *)
bool GetDoComputeEDepInGraphite() const
void msg(const char *fmt,...)
Definition: message.cpp:107
bool GetCreateTrkPlaneH1Output() const
const LBNEVolumePlacementData * Find(const G4String &name, const char *motherName, const char *method) const
void dumpStepCheckVolumeAndFields(const G4Step *)
bool GetDoComputeEDepInArgonGas() const
void SetParentMomentumAtThisProduction(G4ThreeVector mom)
std::string string
Definition: nybbler.cc:12
std::vector< std::ofstream * > fOutPtrsForMarsCmpApr2017
static LBNEVolumePlacements * Instance()
bool GetCreateAlcoveTrackingOutput() const
std::map< int, int > fRZVoxelsData
void FillTrackingPlaneH2Data(const G4Step &aStep)
void CheckInTrackingDetectorDPPlane(const G4Step *theStep)
void GenerateVolumeCrossingsRZMap(const G4Step *theStep)
void CheckInTrackingDetectorPlane(const G4Step *theStep)
void CheckInTargetOutput(const G4Step *theStep)
bool GetCreateTrkPlaneDPOutput() const
G4LogicalVolume * TrkPlnH1Logical
virtual void GetFieldValue(const double Point[3], double *Bfield) const
intermediate_table::const_iterator const_iterator
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:45
G4LogicalVolume * TrkPlnH2Logical
bool GetCreateTargetOutput() const
std::ofstream fOutMuonSculptedAbsorber
void StudyMuonSculptedAbsorberFlux(const G4Step *theStep)
void FillTrackingNtuple(const G4Track &track, LBNETrajectory *currTrajectory)
void StudyCheckMagneticTilts(const G4Step *)
void CheckInTgtEndPlane(const G4Step *theStep)
void FillAlcoveTrackingPlaneData(const G4Step &aStep)
double GetEffectiveInnerRad() const
virtual void GetFieldValue(const double Point[3], double *Bfield) const
void InitiateHadronFluxFromTargetASCII() const
static LBNEAnalysis * getInstance()
T abs(T value)
void StudyPropagation(const G4Step *)
LBNFDeDxMap * GetG4MARSDeDxMap() const
bool GetDoComputeEDepInHorns() const
bool GetCreateTrkPlaneOutput() const
void FillTrackingPlaneData(const G4Step &aStep)
std::string getenv(std::string const &name)
Definition: getenv.cc:15
double z
void StudyParticleThroughHorns(const G4Step *)
static LBNEQuickPiToNuVect * Instance()
void SetDecayCode(G4int decaycode)
std::ofstream fStrHadronFluxFromTargetASCII
void StudyPionsThroughHorn2(const G4Step *)
static int max(int a, int b)
G4LogicalVolume * TrkPlnDPLogical
LBNERunManager * pRunManager
void OpenNtuple(const char *fname)
int sign(double val)
Definition: UtilFunc.cxx:103
void CheckInTrackingDetectorH2Plane(const G4Step *theStep)
std::ofstream fOutStepStudy
void KillNonNuThresholdParticles(const G4Step *theStep)
void FillTrackingPlaneDPData(const G4Step &aStep)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
void CheckHadronsMarsCmpApr2017(const G4Step *theStep)
G4double GetNImpWt() const
int GetVerboseLevel() const
G4LogicalVolume * TrkPlnLogical
virtual void AppendStep(const G4Step *aStep)
void OpenAscii(const char *fname)
void CheckInAlcoveTrackingPlane(const G4Step *theStep)
double GetEffectiveOuterRad() const
void CheckInTrackingDetectorH1Plane(const G4Step *theStep)
void CheckDecay(const G4Step *theStep)
double GetEffectiveSkinDepthFact() const
int GetNumberOfHornsPolycone() const
void FillTargetOutputData(const G4Step &aStep)
void SetMomentumInfoForParticle(const G4Step *theStep)
void FillHadronFluxFromTargetASCII(const G4Step *theStep) const
G4String fOutMuonSculptedAbsorberFilename
void CheckInHornEndPlane(const G4Step *theStep)
QTextStream & endl(QTextStream &s)
void StudyAbsorption(const G4Step *)
virtual void UserSteppingAction(const G4Step *)
void RingHangerPion(int trNum1, int sign, G4ThreeVector p, G4ThreeVector xPos)
bool GetCreateTrkPlaneH2Output() const