LBNESteppingAction.hh

Go to the documentation of this file.

//
// LBNESteppingAction.hh
//

#ifndef LBNESteppingAction_H
#define LBNESteppingAction_H 1
#include <fstream>
#include <iostream>
#include "globals.hh"
#include "G4UserSteppingAction.hh"
#include "G4LogicalVolume.hh"
#include "LBNESteppingActionMessenger.hh"

class G4EventManager;
class LBNEEventAction;
class LBNERunManager;
class TFile;
class TTree;

//
// Investigating why fixing the Ring Hanger diameter gave no change at 3GeV. 
// We will try to correlate the pion that crosses the 

class LBNESteppingAction : public G4UserSteppingAction
{
  
 public:
  LBNESteppingAction();
  virtual ~LBNESteppingAction();
  
  virtual void UserSteppingAction(const G4Step*);

   void KillNonNuThresholdParticles(const G4Step * theStep);
   void CheckDecay(const G4Step * theStep);
   void SetMomentumInfoForParticle(const G4Step * theStep);//A. Bashyal ...ppfx

   void CheckInHornEndPlane(const G4Step * theStep);

   void CheckInTgtEndPlane(const G4Step * theStep);
   
   void OpenAscii(const char *fname); // for Geantino studies. (Absorption studies..) 
   void OpenNtuple(const char *fname); // for Geantino studies. (Absorption studies..) 
                                      // Extended to Particle flux with PoT at exit of Horn1, entrance of Horn2, exit of Horn2 

  
 
private:
   
   LBNERunManager *pRunManager;
   G4EventManager *EvtManager;
   LBNEEventAction *LBNEEvtAct;
   G4LogicalVolume *TrkPlnLogical;   // tracking plane stuff
   G4LogicalVolume *TrkPlnH1Logical; //tracking plane in the end of Horn1
   G4LogicalVolume *TrkPlnH2Logical; //tracking plane in the end of Horn2
   G4LogicalVolume  *DecayPipeHall;
   G4LogicalVolume *TrkPlnDPLogical;   //tracking plane at the end of DecayPipe
//
// This partly "private code" from P. Lebrun, but could be used as template for other studies. 
//  Should not hurt if not activated..
// 
   LBNESteppingActionMessenger *pMessenger;
   std::ofstream fOutStudy;
   std::ofstream fOutStepStudy;
  bool makeSteppingTuple = false;
  TFile *steppingTupleFile;
  TTree *steppingTuple;
  double steppingTuple_x;
  double steppingTuple_y;
  double steppingTuple_z;
  double steppingTuple_density;

   double fKillTrackingThreshold;
   

   // Study Geantino and absorption length patterns in the target/horn . 
   
   double totalAbsDecayChan; // total absorption probability. 
   double totalAbsHorn1Neck; //At various Z locations.. 
   double totalAbsHorn2Entr; 
   double waterAbsDecayChan; 
   double waterAbsHorn1Neck; 
   double waterAbsHorn2Entr; 
   double alumAbsHorn2Entr; // to check inner conductor horm geometry.
   
   int fNConsecutivSmallSteps;
   int fNumTracksKilledAsStuck;
   mutable int fNumStepsCurrentTrack;
   
   bool goneThroughHorn1Neck;
   bool goneThroughHorn2Entr;
   bool doStudyParticleThroughHorns;
   mutable bool doGenerateHadronFluxFromTargetASCII;
   mutable std::ofstream fStrHadronFluxFromTargetASCII;
   mutable bool GenerateMuonSculptedAbsorberFlux;

   mutable bool GenerateMuonLBNEAbsorberFlux;
   mutable std::ofstream fOutMuonSculptedAbsorber;
   mutable G4String fOutMuonSculptedAbsorberFilename;
   G4String fStudyGeantinoMode;
   G4String fKeyVolumeForOutput;
   G4String fKeyVolumeForOutputTo;
   int fEvtIdPrevious; 
   //
   // More code bloat, quick study for Sergei Striganov, April 3 - 2017
   //
   mutable std::vector<std::ofstream*> fOutPtrsForMarsCmpApr2017; 
  
   mutable double fEnergyDepInGraphite; 
   mutable double fEnergyDepInArgonGasHorn1; 
   mutable double fEnergyDepInArgonGasHorn2;
   
   mutable bool fThisParticleGotToHAFront; // for computing the probability that a muon detected 
                                           // at the back of HA was created in the HA.  
   std::map<int,int> fRZVoxelsData; // Yet an other tool to debug geometry problem. 
                                   // contains the multiplicity on volume crossings for Geantinos

   double fGenerateVoxelsZScale ;
   double fGenerateVoxelsZOffset;
   double fGenerateVoxelsZMax;
   double fGenerateVoxelsRScale;
   double fGenerateVoxelsRMax;
   int fGenerateVoxelsIRKMax;
      
   void StudyAbsorption(const G4Step*); 
   void StudyPropagation(const G4Step*); 
   void StudyCheckOverlap(const G4Step*);
   void GenerateVolumeCrossingsRZMap(const G4Step * theStep); 
   // The above allow for a semi-blind test of the geometry,
    // by filling a set of 2D (R-Z plane) of voxels.
   void StudyCheckMagneticTilts(const G4Step*); 
   void dumpStepCheckVolumeAndFields(const G4Step*);
   void StudyParticleThroughHorns(const G4Step*);
   void StudyPionsThroughHorn2(const G4Step*);
  
   
public:
 void CheckInTrackingDetectorH1Plane(const G4Step *theStep);
  void CheckInTrackingDetectorH2Plane(const G4Step *theStep);
  void CheckInTrackingDetectorDPPlane(const G4Step *theStep);
  void CheckInTargetOutput(const G4Step *theStep);
  void CheckInTrackingDetectorPlane(const G4Step *theStep);
  void CheckInAlcoveTrackingPlane(const G4Step *theStep);
   void  SetKillTrackingThreshold(double t) {fKillTrackingThreshold = t;} 
   double GetKillTrackingThreshold() const {return fKillTrackingThreshold;} 
   G4String GetStudyGeantinoMode() const { return fStudyGeantinoMode; } 
   void SetStudyGeantinoMode(G4String v) {fStudyGeantinoMode = v; }
   void SetKeyVolumeForOutput(G4String v) {fKeyVolumeForOutput = v; }
   void SetKeyVolumeForOutputTo(G4String v) {fKeyVolumeForOutputTo = v; }
   void SetStudyParticleThroughHorns(bool t) {doStudyParticleThroughHorns = t;}
   void InitiateHadronFluxFromTargetASCII() const; // Pseudo Const.. 
   void FillHadronFluxFromTargetASCII(const G4Step *theStep) const;// Pseudo Const.. 
   inline int GetNumTracksKilledAsStuck() const { return fNumTracksKilledAsStuck;}
   inline void ResetNumSteps() const {fNumStepsCurrentTrack = 0;} // Not const on mutable. Screw you.. 
   inline void ResetEnergyDepInGraphite() const { fEnergyDepInGraphite = 0.;} // Same 
   inline void ResetEnergyDepInArgonGas() const {  
                fEnergyDepInArgonGasHorn1= 0.;  fEnergyDepInArgonGasHorn2= 0.; } // Same 
   inline double GetEnergyDepInGraphite() const { return fEnergyDepInGraphite; }
   inline double GetEnergyDepInArgonGasH1() const { return fEnergyDepInArgonGasHorn1; }
   inline double GetEnergyDepInArgonGasH2() const { return fEnergyDepInArgonGasHorn2; }
   void SetGenerateMuonSculptedAbsorberFlux(bool t){GenerateMuonSculptedAbsorberFlux = t;}
   void SetMuonSculptedAbsorberFluxFilename(G4String str = "./MuonFluxAtSculptedAbsorber.txt")
     {fOutMuonSculptedAbsorberFilename = str;  }
   void StudyMuonSculptedAbsorberFlux(const G4Step *theStep); // Also for LBNE Absorber, saving a method. 
   void StudySculptedAbsorberTrackingFlux(const G4Step *theStep); // Also for LBNE Absorber, saving a method. 
   void CloseSculptedAbsorberTrackingFlux();
   inline void ResetFlagParticleGotToHAFront () const {fThisParticleGotToHAFront=false; } // pseudo const, for HA muon 
   // contamination studies.
   
 private:
   void CheckHadronsMarsCmpApr2017(const G4Step *theStep);
    
};

#endif