LBNFWormPlotter.cc

Go to the documentation of this file.

// Code that implements the "worm" plotter using G4 intersections and ROOT graphs.
// John Back, March 2016

#include "LBNFWormPlotter.hh"

#include "G4EventManager.hh"

#include "G4UnitsTable.hh"
#include "G4VFigureFileMaker.hh"

#include "G4GeometryManager.hh"
#include "G4SDManager.hh"
#include "G4StateManager.hh"
#include "G4Event.hh"
#include "G4TrajectoryContainer.hh"

#include "G4VisAttributes.hh"
#include "G4UImanager.hh"
#include "G4TransportationManager.hh"
#include "G4RegionStore.hh"
#include "G4ProductionCutsTable.hh"
#include "G4VVisManager.hh"

#include "G4ProcessManager.hh"
#include "G4ProcessVector.hh"
#include "G4Geantino.hh"
#include "G4RayShooter.hh"

#include "G4UserLimits.hh"

#include "LBNFWormPlotMessenger.hh"
#include "LBNETrajectory.hh"
#include "G4VTrajectoryPoint.hh"

#include "TGraph.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TH2.h"

LBNFWormPlotter::LBNFWormPlotter() :
    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() {

    delete theRayShooter;
    delete theMessenger;

}

void LBNFWormPlotter::CreatePlot(const G4String& fileName) {

    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& rayStart,
                                                              const G4ThreeVector& rayDirection) const {

    // 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;
    
}

G4bool LBNFWormPlotter::insidePlot(const G4ThreeVector& thePoint) const {
    
    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;
    
}