Part 9

//construction.cc


#include  "construction.hh"

 MyDetectorConstruction::MyDetectorConstruction()
 {}
 
  MyDetectorConstruction::~MyDetectorConstruction()
  {}
   G4VPhysicalVolume *MyDetectorConstruction::Construct()
   {
   
   G4NistManager *nist = G4NistManager::Instance();
   
   G4Material *SiO2 =new G4Material ("SiO2", 2.201*g/cm3,2);
   SiO2->AddElement(nist->FindOrBuildElement("Si"),1);
   SiO2->AddElement(nist->FindOrBuildElement("O"),2);
   
   G4Material *H2O =new G4Material ("H2O", 1.000*g/cm3,2);
   H2O->AddElement(nist->FindOrBuildElement("H"),2);
   H2O->AddElement(nist->FindOrBuildElement("O"),1);
   
   G4Element *C = nist->FindOrBuildElement("C");
   
   G4Material *Aerogel = new G4Material("Aerogel",0.200*g/cm3,3);
   Aerogel -> AddMaterial(SiO2,62.5*perCent);
   Aerogel-> AddMaterial(H2O,37.4*perCent);
   Aerogel-> AddElement(C,0.1*perCent);
   
   G4double energy[2] = {1.239841939*eV/0.2 ,1.239841939*eV/0.9};
   G4double rindexAerogel[2]={1.1,1.1};
    G4double rindexWorld[2]={1.0,1.0};
   
   
   
   G4MaterialPropertiesTable *mptAerogel =new G4MaterialPropertiesTable();
   mptAerogel->AddProperty("RINDEX",energy,rindexAerogel,2);
   
   
   
   
   Aerogel->SetMaterialPropertiesTable(mptAerogel);
   
   G4Material *worldMat= nist->FindOrBuildMaterial("G4_AIR");
   
   G4MaterialPropertiesTable *mptWorld= new G4MaterialPropertiesTable();
   mptWorld->AddProperty("RINDEX",energy,rindexWorld, 2);
   
   
   worldMat->SetMaterialPropertiesTable(mptWorld);
   
   
   
   
   
   
    
   G4Box *solidworld=new G4Box("solidworld",0.5*m,0.5*m,0.5*m);
   
   G4LogicalVolume *logicworld=new G4LogicalVolume(solidworld,worldMat,"logicWorld");
   
   G4VPhysicalVolume *physWorld= new G4PVPlacement(0, G4ThreeVector(0.,0.,0.),logicworld,"physworld",0,false,0,true);
   
   G4Box *solidRadiator =new G4Box ("solidRadiator", 0.4*m ,0.4*m, 0.05*m);
   
   G4LogicalVolume *logicRadiator =new G4LogicalVolume(solidRadiator,Aerogel, "logicalRadiator");
   
   G4VPhysicalVolume *physRadiator= new G4PVPlacement(0, G4ThreeVector(0.,0.,0.),logicRadiator,"physRadiator",logicworld,false,0,true);
   
   
   G4Box *solidDetector =new G4Box("solidDetector",0.005*m,0.005*m,0.005*m); //solid detector define
   
   
   logicDetector =new G4LogicalVolume(solidDetector, worldMat , "logicDetector");
   
   
   for (G4int i=0;i<100;i++)
   {
   for(G4int j = 0; j<100;j++)
   
   {
   G4VPhysicalVolume *phyDetector=new G4PVPlacement(0, G4ThreeVector(-0.5*m+(i+0.5)*m/100,-0.5*m+(j+0.5)*m/100,0.49*m),logicDetector,"physDetector",logicworld,false,j+i*100,true);
   
   }
   
   
   
   }
   
   
   
   
   return physWorld;
   
   }
   
   void MyDetectorConstruction::ConstructSDandField()
   {
   MySensitiveDetector *sensDet=new MySensitiveDetector("SensitiveDetector");
   
   logicDetector->SetSensitiveDetector(sensDet);
   
   }





//generator.cc

#include "generator.hh"

MyPrimaryGenerator::MyPrimaryGenerator()
{
    fParticleGun = new G4ParticleGun(1);
G4ParticleTable *particleTable=G4ParticleTable::GetParticleTable();
  
  G4String particleName="proton";
  G4ParticleDefinition *particle =particleTable->FindParticle("proton");
  
  G4ThreeVector pos(0.,0.,0.);
  
  G4ThreeVector mom(0.,0.,1.);
  
    fParticleGun->SetParticlePosition(pos);
        fParticleGun->SetParticleMomentumDirection(mom);
          fParticleGun->SetParticleMomentum(100. *GeV);
         fParticleGun->SetParticleDefinition(particle);
   // G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
  //  G4ParticleDefinition *particle = particleTable->FindParticle("chargedgeantino");

  //  G4ThreeVector pos(0.,0.,-20.*km);
  //  G4ThreeVector mom(0.,0.,1.);

  //  fParticleGun->SetParticlePosition(pos);
   // fParticleGun->SetParticleMomentumDirection(mom);
  //  fParticleGun->SetParticleMomentum(100.*GeV);
  //  fParticleGun->SetParticleDefinition(particle);
}

MyPrimaryGenerator::~MyPrimaryGenerator()
{
    delete fParticleGun;
}

void MyPrimaryGenerator::GeneratePrimaries(G4Event *anEvent)
{
  
         
           fParticleGun->GeneratePrimaryVertex(anEvent);
       
         
        
  
  
   // G4ParticleDefinition* particle = fParticleGun->GetParticleDefinition();

    //if(particle == G4ChargedGeantino::ChargedGeantino())
    //{
      //  G4int Z = 27;
      //  G4int A = 60;

      //  G4double charge   = 0.*eplus;
      //  G4double energy = 0.*keV;

      //  G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon(Z,A,energy);
      //  fParticleGun->SetParticleDefinition(ion);
      //  fParticleGun->SetParticleCharge(charge);
   // }

   // fParticleGun->GeneratePrimaryVertex(anEvent);
}





//generator.hh

#ifndef GENERATOR_HH
#define GENERATOR_HH

#include "G4VUserPrimaryGeneratorAction.hh"

#include "G4ParticleGun.hh"
#include "G4SystemOfUnits.hh"
#include "G4ParticleTable.hh"
//#include "G4IonTable.hh"
//#include "G4ChargedGeantino.hh"
//#include "Randomize.hh"

class MyPrimaryGenerator : public G4VUserPrimaryGeneratorAction
{
public:
    MyPrimaryGenerator();
    ~MyPrimaryGenerator();

    virtual void GeneratePrimaries(G4Event*);

private:
    G4ParticleGun *fParticleGun;
};

#endif



//detector.hh

#ifndef DETECTOR_HH
#define DETECTOR_HH

#include "G4VSensitiveDetector.hh"
#include "G4RunManager.hh"
#include "G4SystemOfUnits.hh"

#include "g4root.hh"

class MySensitiveDetector : public G4VSensitiveDetector
{
public:
    MySensitiveDetector(G4String);
    ~MySensitiveDetector();
    
private:
    virtual G4bool ProcessHits(G4Step *, G4TouchableHistory *);
    
   // G4PhysicsOrderedFreeVector *quEff;
};

#endif






//detector.cc


#include "detector.hh"
#include "g4root.hh"
MySensitiveDetector::MySensitiveDetector(G4String name) : G4VSensitiveDetector(name)
{
    /*quEff = new G4PhysicsOrderedFreeVector();

    std::ifstream datafile;
    datafile.open("eff.dat");

    while(1)
    {
        G4double wlen, queff;

        datafile >> wlen >> queff;

        if(datafile.eof())
            break;

        G4cout << wlen << " " << queff << std::endl;

        quEff->InsertValues(wlen, queff/100.);*/
    }

//    datafile.close();

//    quEff->SetSpline(false);}

MySensitiveDetector::~MySensitiveDetector()
{}

G4bool MySensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist)


{
   G4Track *track = aStep->GetTrack();

    track->SetTrackStatus(fStopAndKill); //when photon hits the detector it kills photon's trajectory

    G4StepPoint *preStepPoint = aStep->GetPreStepPoint();   //photon's start point
    G4StepPoint *postStepPoint = aStep->GetPostStepPoint(); //photon's end point

    G4ThreeVector posPhoton = preStepPoint->GetPosition();
   /* G4ThreeVector momPhoton = preStepPoint->GetMomentum();

    G4double time = preStepPoint->GetGlobalTime();

    G4double wlen = (1.239841939*eV/momPhoton.mag())*1E+03;*/

   // G4cout << "Photon position: " << posPhoton << G4endl;
    
    
    
    
    

   const G4VTouchable *touchable = aStep->GetPreStepPoint()->GetTouchable();

    G4int copyNo = touchable->GetCopyNumber();

    //G4cout << "Copy number: " << copyNo << G4endl;

    G4VPhysicalVolume *physVol = touchable->GetVolume();
    G4ThreeVector posDetector = physVol->GetTranslation();
      G4cout << "Detector position: " << posDetector << G4endl;

   // #ifndef G4MULTITHREADED
      
        //G4cout << "Photon wavelength: " << wlen << G4endl;
    
    
    /*#endif*/

    G4int evt = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();

    G4AnalysisManager *man = G4AnalysisManager::Instance();

   // man->FillNtupleIColumn(0, 0, evt);
    //man->FillNtupleDColumn(0, 1, posPhoton[0]);
    //man->FillNtupleDColumn(0, 2, posPhoton[1]);
    //man->FillNtupleDColumn(0, 3, posPhoton[2]);
   // man->FillNtupleDColumn(0, 4, time);
    //man->FillNtupleDColumn(0, 5, wcslen);
   // man->AddNtupleRow(0);

   /* if(G4UniformRand() < quEff->Value(wlen))*/
    {
        man->FillNtupleIColumn(0, evt);
        man->FillNtupleDColumn(1,  posDetector[0]);
        man->FillNtupleDColumn(2,  posDetector[1]);
        man->FillNtupleDColumn(3, posDetector[2]);
        man->AddNtupleRow(0);
    }

    return true;
}








//sim.cc


#include <iostream>

#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4VisManager.hh"
#include "G4VisExecutive.hh"
#include "G4UIExecutive.hh"
#include "construction.cc"
#include "physics.cc"
#include "action.cc"
 #include "generator.cc"
 #include "detector.cc"
 #include "run.cc"

int  main(int argc, char** argv)

{ G4RunManager *runManager = new G4RunManager();
runManager->SetUserInitialization(new MyDetectorConstruction());
runManager->SetUserInitialization(new MyPhysicsList());
runManager->SetUserInitialization(new MyActionInitialization());


runManager->Initialize();

 G4UIExecutive *ui = 0;
 
 if(argc==1){
 ui= new G4UIExecutive(argc, argv);
}
G4VisManager *vismanager =new G4VisExecutive ();
vismanager->Initialize();
G4UImanager *UImanager = G4UImanager::GetUIpointer();
if(ui)
{

 UImanager->ApplyCommand("/control/execute vis.mac");
ui->SessionStart();
}
else
{
G4String command="/control/execute ";
G4String fileName=argv[1];
UImanager->ApplyCommand(command+fileName);
}
return 0;

}






//vis.mac


/run/initialize
/vis/open OGL
/vis/viewer/set/viewpointVector 1 1 1
/vis/drawVolume
/vis/viewer/set/autoRefresh true
/vis/scene/add/trajectories smooth
/vis/scene/endOfEventAction accumulate
/vis/scene/add/scale 10 cm
/vis/scene/add/axes
/vis/scene/add/eventID
/*/gun/particle proton
/vis/ogl/set/displayListLimit 5000000*/





//CMakeLists.txt

cmake_minimum_required(VERSION 2.6 FATAL_ERROR)

project(Simulation)

find_package(Geant4 REQUIRED  ui_all vis_all)

include(${Geant4_USE_FILE})

file(GLOB sources ${PROJECT_SOURCE_DIR}/*.cc)
file(GLOB headers ${PROJECT_SOURCE_DIR}/*.hh)

file(GLOB MACRO_FILES

"*.mac"

)

file(COPY ${MACRO_FILES} DESTINATION ${PROJECT_BINARY_DIR}) 

add_executable(sim sim.cc ${source} ${headers})
target_link_libraries(sim ${Geant4_LIBRARIES})

add_custom_target(Simulation DEPENDS sim)




//run.mac


/run/initialize
/gun/momentumAmp 0.5 GeV
/run/beamOn 100
/gun/momentumAmp 1 GeV
/run/beamOn 100
/gun/momentumAmp 5 GeV
/run/beamOn 100




//run.cc

#include "run.hh"

MyRunAction::MyRunAction()
{ G4AnalysisManager *man = G4AnalysisManager::Instance();
   man->CreateNtuple("Hits", "Hits");
    man->CreateNtupleIColumn("fEvent");
    man->CreateNtupleDColumn("fX");
    man->CreateNtupleDColumn("fY");
    man->CreateNtupleDColumn("fZ");
    man->FinishNtuple(0);
}

MyRunAction::~MyRunAction()
{}

void MyRunAction::BeginOfRunAction(const G4Run* run)
{ 
  G4AnalysisManager *man = G4AnalysisManager::Instance();
  
     G4int runID = run->GetRunID();

    std::stringstream strRunID;
    strRunID << runID;

    man->OpenFile("output"+strRunID.str()+".root");

 
  //man->CreateNtuple("Photons", "Photons");
    //man->CreateNtupleIColumn("fEvent");  //N means integer
    //man->CreateNtupleDColumn("fX");
    //man->CreateNtupleDColumn("fY");
    //man->CreateNtupleDColumn("fZ");
   // man->CreateNtupleDColumn("fT");
   // man->CreateNtupleDColumn("fWlen");
//    man->FinishNtuple(0);

   

    //man->CreateNtuple("Scoring", "Scoring");
    //man->CreateNtupleDColumn("fEdep");
    //man->FinishNtuple(2);




/*
    G4AnalysisManager *man = G4AnalysisManager::Instance();

    G4int runID = run->GetRunID();

    std::stringstream strRunID;
    strRunID << runID;

    man->OpenFile("output"+strRunID.str()+".root");
*/
}

void MyRunAction::EndOfRunAction(const G4Run*)
{
    G4AnalysisManager *man = G4AnalysisManager::Instance();

    man->Write();
    man->CloseFile(); 
}




//run.hh





#ifndef RUN_HH  //header file
#define RUN_HH

#include "G4UserRunAction.hh"
#include "G4Run.hh"

#include "g4root.hh"

class MyRunAction : public G4UserRunAction
{
public:
    MyRunAction();
    ~MyRunAction();
    
    virtual void BeginOfRunAction(const G4Run*);
    virtual void EndOfRunAction(const G4Run*);
};

#endif