#include <LBNFWormPlotter.hh>

Public Member Functions
	LBNFWormPlotter ()

virtual	~LBNFWormPlotter ()

void	CreatePlot (const G4String &fileName)

G4ThreeVector	GetXAxis () const

G4ThreeVector	GetYAxis () const

G4ThreeVector	GetZAxis () const

G4double	GetXMin () const

G4double	GetXMax () const

G4double	GetYMin () const

G4double	GetYMax () const

G4double	GetZ0 () const

G4String	GetXName () const

G4String	GetYName () const

void	SetXAxis (const G4ThreeVector &vect)

void	SetYAxis (const G4ThreeVector &vect)

void	SetXName (const G4String &word)

void	SetYName (const G4String &word)

void	SetXMin (const G4double &val)

void	SetXMax (const G4double &val)

void	SetYMin (const G4double &val)

void	SetYMax (const G4double &val)

void	SetZ0 (const G4double &val)

Protected Member Functions
std::vector< G4ThreeVector >	findIntersections (const G4ThreeVector &rayPosition, const G4ThreeVector &rayDirection) const

G4bool	insidePlot (const G4ThreeVector &thePoint) const

Protected Attributes
G4ThreeVector	xAxis

G4ThreeVector	yAxis

G4ThreeVector	zAxis

G4double	xMin

G4double	xMax

G4double	yMin

G4double	yMax

G4double	z0

G4bool	gotBounds

G4double	ABSq

G4double	ADSq

G4ThreeVector	AB

G4ThreeVector	AD

G4ThreeVector	PointA

G4String	xName

G4String	yName

G4RayShooter *	theRayShooter

LBNFWormPlotMessenger *	theMessenger

G4EventManager *	theEventManager

Detailed Description

Definition at line 49 of file LBNFWormPlotter.hh.

Constructor & Destructor Documentation

LBNFWormPlotter::LBNFWormPlotter ( )

Definition at line 42 of file LBNFWormPlotter.cc.

                                  :
     xAxis(G4ThreeVector(0.0, 0.0, 1.0)), // along z
     yAxis(G4ThreeVector(0.0, 1.0, 0.0)), // along y
     zAxis(G4ThreeVector(-1.0, 0.0, 0.0)), // cross-product of x and y
     xMin(-1.0*CLHEP::m),
     xMax(1.0*CLHEP::m),
     yMin(-1.0*CLHEP::m),
     yMax(1.0*CLHEP::m),
     z0(0.0*CLHEP::m),
     gotBounds(false),
     ABSq(0.0),
     ADSq(0.0),
     AB(G4ThreeVector(0.0, 0.0, 0.0)),
     AD(G4ThreeVector(0.0, 0.0, 0.0)),
     PointA(G4ThreeVector(0.0, 0.0, 0.0)),
     xName(""),
     yName(""),
     theRayShooter(0),
     theMessenger(0),
     theEventManager(0)
 {
 
     theRayShooter = new G4RayShooter();
     theMessenger = LBNFWormPlotMessenger::GetInstance(this);
     theEventManager = G4EventManager::GetEventManager();
 
 }

LBNFWormPlotter::~LBNFWormPlotter ( )

virtual

Definition at line 70 of file LBNFWormPlotter.cc.

                                   {
 
     delete theRayShooter;
     delete theMessenger;
 
 }

Member Function Documentation

void LBNFWormPlotter::CreatePlot ( const G4String & fileName )

Definition at line 77 of file LBNFWormPlotter.cc.

                                                          {
 
     zAxis = xAxis.cross(yAxis);
 
     G4cout<<"Creating worm plot: "<<fileName<<G4endl;
     // Print plot dimension information
     G4cout<<"xAxis unit vector = "<<xAxis<<G4endl;
     G4cout<<"yAxis unit vector = "<<yAxis<<G4endl;
     G4cout<<"zAxis unit vector = "<<zAxis<<G4endl;
     G4cout<<"x limits = ("<<xMin/CLHEP::cm<<", "<<xMax/CLHEP::cm<<") cm"<<G4endl;
     G4cout<<"y limits = ("<<yMin/CLHEP::cm<<", "<<yMax/CLHEP::cm<<") cm"<<G4endl;
     G4cout<<"z0 position = "<<z0/CLHEP::cm<<" cm"<<G4endl;
 
     // Check if we have got the squared lengths of the plot axes
     if (gotBounds == false) {    
         
         AB = (xMax - xMin)*xAxis;
         ABSq = AB.mag2();
         AD = (yMax - yMin)*yAxis;
         ADSq = AD.mag2();
         
         PointA = xMin*xAxis + yMin*yAxis;
         
         gotBounds = true;
         
         //G4cout<<"*** ABSq = "<<ABSq<<", ADSq = "<<ADSq<<", PointA = "<<PointA<<G4endl;
         
     }
 
     // Use 1,000 lines per axis
     G4int nX(1000), nY(1000);
     G4double dx = (xMax - xMin)/(nX*1.0);
     G4double dy = (yMax - yMin)/(nY*1.0);
     G4cout<<"dx = "<<dx/CLHEP::cm<<", dy = "<<dy/CLHEP::cm<<" cm"<<G4endl;
     
     // Confirm process(es) of Geantino is initialized
     G4VPhysicalVolume* pWorld = G4TransportationManager::GetTransportationManager()->
         GetNavigatorForTracking()->GetWorldVolume();
     G4RegionStore::GetInstance()->UpdateMaterialList(pWorld);
     
     G4ProductionCutsTable::GetProductionCutsTable()->UpdateCoupleTable(pWorld);
     
     G4ProcessVector* pVector = G4Geantino::GeantinoDefinition()->GetProcessManager()->GetProcessList();
     for (G4int j=0; j < pVector->size(); ++j) {
         (*pVector)[j]->BuildPhysicsTable(*(G4Geantino::GeantinoDefinition()));
     }
     
     // Close geometry and set the application state
     G4GeometryManager* geomManager = G4GeometryManager::GetInstance();
     geomManager->OpenGeometry();
     geomManager->CloseGeometry(1,0);
     
     G4ThreeVector center(0,0,0);
     G4Navigator* navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
     navigator->LocateGlobalPointAndSetup(center,0,false);
     
     G4StateManager* theStateMan = G4StateManager::GetStateManager();
     theStateMan->SetNewState(G4State_GeomClosed); 
     
     TGraph* theGraph = new TGraph();
     theGraph->SetMarkerStyle(21);
     theGraph->SetMarkerSize(0.02);
     
     // Find intersections along a horizontal grid of lines, starting from the bottom and going up
     G4ThreeVector x0 = xMin*xAxis + z0*zAxis;
     for (G4int iy = 0; iy < nY; iy++) {
         
         // Starting position for the ray
         G4ThreeVector yVect = (iy*dy + yMin)*yAxis;
         G4ThreeVector rayStart = yVect + x0;
         //G4cout<<"Start ray position = "<<rayStart<<G4endl;
         
         std::vector<G4ThreeVector> hIntersections = this->findIntersections(rayStart, xAxis);
         
         G4int nI = hIntersections.size();
         //G4cout<<"List of "<<nI<<" horizontal intersections:"<<G4endl;
         
         for (G4int j = 0; j < nI; j++) {
             
             G4ThreeVector hP = hIntersections[j];
             
             //G4cout<<hP<<G4endl;
             
             // Store points in the graph
             theGraph->SetPoint(theGraph->GetN(), hP.dot(xAxis)/CLHEP::cm, hP.dot(yAxis)/CLHEP::cm);
             
         }
         
     }
     
     // Find intersections along a vertical grid of lines
     G4ThreeVector y0 = yMin*yAxis + z0*zAxis;
     for (G4int ix = 0; ix < nX; ix++) {
         
         // Starting position for the ray
         G4ThreeVector xVect = (ix*dx + xMin)*xAxis;
         G4ThreeVector rayStart = xVect + y0;
         //G4cout<<"Start ray position = "<<rayStart<<G4endl;
         
         std::vector<G4ThreeVector> vIntersections = this->findIntersections(rayStart, yAxis);
         
         G4int nI = vIntersections.size();
         //G4cout<<"List of "<<nI<<" vertical intersections:"<<G4endl;
         for (G4int j = 0; j < nI; j++) {
             
             G4ThreeVector vP = vIntersections[j];
             //G4cout<<vP<<G4endl;
             
             // Store points in the graph
             theGraph->SetPoint(theGraph->GetN(), vP.dot(xAxis)/CLHEP::cm, vP.dot(yAxis)/CLHEP::cm);
             
         }
         
     }
     
     TCanvas theCanvas("theCanvas", "", 900, 700);
     gROOT->SetStyle("Plain");
     gStyle->SetOptStat(0);
     theCanvas.UseCurrentStyle();
     
     theCanvas.SetGrid();
     
     // Create a 2d histogram to specify the drawing axes, then overlay the graph
     TH2D* theHist = new TH2D("axes", "", 2, xMin/CLHEP::cm, xMax/CLHEP::cm, 2, yMin/CLHEP::cm, yMax/CLHEP::cm);
     theHist->SetDirectory(0);
     TString xLabel(xName), yLabel(yName);
     xLabel.ReplaceAll("\"", "");
     yLabel.ReplaceAll("\"", "");
     theHist->GetXaxis()->SetTitle(xLabel);
     theHist->GetYaxis()->SetTitle(yLabel);
     theHist->Draw();
     
     theGraph->Draw("psame");
     theCanvas.Update();
     
     theCanvas.Print(fileName);
     
     delete theGraph;
     delete theHist;
     
     // Restore state to idle
     theStateMan->SetNewState(G4State_Idle);   
 
 }

std::vector< G4ThreeVector > LBNFWormPlotter::findIntersections	(	const G4ThreeVector &	rayPosition,
		const G4ThreeVector &	rayDirection
	)		const

protected

Definition at line 222 of file LBNFWormPlotter.cc.

                                                                                                        {
 
     // Vector of intersection points
     std::vector<G4ThreeVector> theIntersections;
     
     G4int iEvent(1);
     G4Event* anEvent = new G4Event(iEvent);
     
     // Use geantino "ray" to find the intersections with the geometry along the given direction.
     // This will find intersections that may extend beyond the "plot boundary", so only store
     // the points we need
     theRayShooter->Shoot(anEvent, rayStart, rayDirection.unit());
     theEventManager->ProcessOneEvent(anEvent);
     
     // Fractional distance cut (|P1-P0|/|P0|) to find steps across boundaries
     G4double fracCut(1e-6);
     
     G4TrajectoryContainer* trajectoryContainer = anEvent->GetTrajectoryContainer();
     if (trajectoryContainer) {
         
         //G4int nTraj = trajectoryContainer->entries();
         //G4cout<<"nTraj = "<<nTraj<<G4endl;
         
         LBNETrajectory* trajectory = dynamic_cast<LBNETrajectory*>( (*trajectoryContainer)[0] );
         
         if (trajectory) {
             
             G4int nPoints = trajectory->GetPointEntries();
             //G4cout<<"Number of points along trajectory = "<<nPoints<<G4endl;
             
             G4int iP(0);
             G4int pInt(-1);
             
             for (iP = 0; iP < nPoints; iP++) {
                 
                 G4VTrajectoryPoint* thePoint = trajectory->GetPoint(iP);
                 
                 //G4cout<<"Trajectory point "<<iP<<" pre-step volume is "<<trajectory->GetPreStepVolumeName(iP)<<G4endl;
                 
                 if (thePoint) {
                     
                     G4ThreeVector thePos = thePoint->GetPosition();
                     // First, check to see if the point is inside the "plot square". If not, ignore it.
                     G4bool inside = this->insidePlot(thePos);
                     
                     if (inside) {
                         
                         //G4cout<<"Inside? Yes"<<G4endl;
                         pInt++;
                         //G4cout<<"The point "<<pInt<<" = "<<thePos<<G4endl;      
                         
                         // Usually we get two very nearby step points on either side of a volume region change.
                         // Check the previous "intersection" point and if its very close to the current point, average them
                         if (pInt == 0) {
                             
                             //G4cout<<"Adding first point to theIntersections"<<G4endl;
                             
                             theIntersections.push_back(thePos);
                             
                         } else {
                             
                             G4int pInt2 = pInt - 1;
                             G4ThreeVector prevPos = theIntersections[pInt2];
                             G4ThreeVector diffPos = thePos - prevPos;
                             G4double d1 = prevPos.mag();
                             G4double fracD(1.0);
                             if (d1 > 0.0) {fracD = diffPos.mag()/d1;}
                             
                             //G4cout<<"pInt2 = "<<pInt2<<"; prevPos = "<<prevPos<<"; fracD = "<<fracD<<G4endl;
                             
                             if (fracD < fracCut) {
                                 
                                 // The current point is very close to the previous one. Average them
                                 // and update the previous point stored in the vector
                                 G4ThreeVector meanPos = 0.5*(prevPos + thePos);
                                 //G4cout<<"Point - currentPoint fracD = "<<fracD<<"; average = "<<meanPos<<G4endl;
                                 //G4cout<<"Intersection "<<pInt2<<" = "<<theIntersections[pInt2]<<G4endl;
                                 G4ThreeVector& vectPos = theIntersections[pInt2];
                                 vectPos.setX(meanPos.x());
                                 vectPos.setY(meanPos.y());
                                 vectPos.setZ(meanPos.z());
                                 //G4cout<<"New Intersection "<<pInt2<<" = "<<theIntersections[pInt2]<<G4endl;
                                 
                                 // Decrement pInt, since we have merged the two points
                                 pInt--;
                                 
                             } else {
                                 
                                 // The current point is not close to the previous one. Add it to the vector
                                 //G4cout<<"Adding "<<thePos<<G4endl;
                                 theIntersections.push_back(thePos);
                                 
                             } // d1 > 0.0
                             
                         } // Start comparison with previous point
                         
                     } // inside
                     
                 } // thePoint
                 
             } // iP loop
             
         } // trajectory
         
     } // trajectoryContainer
     
     delete anEvent;
     
     return theIntersections;
     
 }

G4ThreeVector LBNFWormPlotter::GetXAxis ( ) const

inline

Definition at line 62 of file LBNFWormPlotter.hh.

62 {return xAxis;}

LBNFWormPlotter::xAxis

G4ThreeVector xAxis

Definition: LBNFWormPlotter.hh:89

G4double LBNFWormPlotter::GetXMax ( ) const

inline

Definition at line 67 of file LBNFWormPlotter.hh.

67 {return xMax;}

LBNFWormPlotter::xMax

G4double xMax

Definition: LBNFWormPlotter.hh:93

G4double LBNFWormPlotter::GetXMin ( ) const

inline

Definition at line 66 of file LBNFWormPlotter.hh.

66 {return xMin;}

LBNFWormPlotter::xMin

G4double xMin

Definition: LBNFWormPlotter.hh:92

G4String LBNFWormPlotter::GetXName ( ) const

inline

Definition at line 72 of file LBNFWormPlotter.hh.

72 {return xName;}

LBNFWormPlotter::xName

G4String xName

Definition: LBNFWormPlotter.hh:104

G4ThreeVector LBNFWormPlotter::GetYAxis ( ) const

inline

Definition at line 63 of file LBNFWormPlotter.hh.

63 {return yAxis;}

LBNFWormPlotter::yAxis

G4ThreeVector yAxis

Definition: LBNFWormPlotter.hh:90

G4double LBNFWormPlotter::GetYMax ( ) const

inline

Definition at line 69 of file LBNFWormPlotter.hh.

69 {return yMax;}

LBNFWormPlotter::yMax

G4double yMax

Definition: LBNFWormPlotter.hh:95

G4double LBNFWormPlotter::GetYMin ( ) const

inline

Definition at line 68 of file LBNFWormPlotter.hh.

68 {return yMin;}

LBNFWormPlotter::yMin

G4double yMin

Definition: LBNFWormPlotter.hh:94

G4String LBNFWormPlotter::GetYName ( ) const

inline

Definition at line 73 of file LBNFWormPlotter.hh.

73 {return yName;}

LBNFWormPlotter::yName

G4String yName

Definition: LBNFWormPlotter.hh:105

G4double LBNFWormPlotter::GetZ0 ( ) const

inline

Definition at line 70 of file LBNFWormPlotter.hh.

70 {return z0;}

LBNFWormPlotter::z0

G4double z0

Definition: LBNFWormPlotter.hh:96

G4ThreeVector LBNFWormPlotter::GetZAxis ( ) const

inline

Definition at line 64 of file LBNFWormPlotter.hh.

64 {return zAxis;}

LBNFWormPlotter::zAxis

G4ThreeVector zAxis

Definition: LBNFWormPlotter.hh:91

G4bool LBNFWormPlotter::insidePlot ( const G4ThreeVector & thePoint ) const

protected

Definition at line 335 of file LBNFWormPlotter.cc.

                                                                       {
     
     G4bool inside(false);
     
     G4ThreeVector AP = thePoint - PointA;
     G4double prod1 = AP.dot(AB);
     
     //G4cout<<"prod1 = "<<prod1<<"; ABSq = "<<ABSq<<G4endl;
     
     if (prod1 > 0.0 && prod1 < ABSq) {
         
         G4double prod2 = AP.dot(AD);
         
         //G4cout<<"prod2 = "<<prod2<<"; ADSq = "<<ADSq<<G4endl;
         
         if (prod2 > 0.0 && prod2 < ADSq) {
             
             inside = true;
             
             //G4cout<<"Point "<<thePoint<<" is inside the plot boundary"<<G4endl;
             
         }
         
     }
     
     //G4cout<<"insidePlot:: inside = "<<(G4int)inside<<G4endl;
     
     return inside;
     
 }