Input gammas, event

[janmajay]

I want to simulate Cerenkov radiation. I want to shoot gammas to water. There are .txt file which has 100000 gammas of different energies. I wrote the following code

ReadFiles();

G4ThreeVector pos(0., 0., 0.);
//G4ThreeVector mon(0., 0., 1.);
G4ThreeVector mon = G4RandomDirection();
for(auto a =0;a<100000;a++){
G4ParticleDefinition* particle = G4ParticleTable::GetParticleTable()->FindParticle("gamma");
 fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticleEnergy(gammasv[a]*keV);
fParticleGun->SetParticlePosition(pos);
fParticleGun->SetParticleMomentumDirection(mon);

fParticleGun->GeneratePrimaryVertex(anEvent); 

}
This creates only one event. I want 100000 gammas as 100000 events.

[mkelsey]

That looks like the XXXGeneratePrimaryVertex()XXX GeneratePrimaries() code (I’m guessing, because you didn’t show us). The stuff inside the loop should be the only stuff in your XXXGeneratePrimaryVertex()XXX GeneratePrimaries() function.

Read the data into your array at the beginning of the run. Don’t bother with the loop at all. In XXXGeneratePrimaryVertex()XXX GeneratePrimaries(), get the event ID and use that to index into your array.

Do you want each gamma thrown in a random direction, or do you want all of them going the same way?

[janmajay]

Yes, this is GeneratePrimaryVertex() , the .hh file looks like this
class MyPrimaryGenerator : public G4VUserPrimaryGeneratorAction
{
public:
MyPrimaryGenerator();
~MyPrimaryGenerator();

virtual void GeneratePrimaries(G4Event*);
 //G4GeneralParticleSource* GetParticleGun() { return fParticleGun; };//for gps
void ReadFiles();

private:
G4ParticleGun fParticleGun;
//G4GeneralParticleSource fParticleGun; //new gun
std::vector gammasv;
};
. I am sorry to say, I did not understand, how should I loop with eventID, and if I will not take care of loop inside the GeneratePrimaryVertex(), how should I generate all those gammas?

[mkelsey]

First, my mistake typing: “GeneratePrimaryVertex()” is the function we use in our own different generators. You’re correct that the function is GeneratePrimaries(); I’ve edited my initial response above.

Notice that GeneratePrimaries() takes a G4Event pointer as argument. It is called at the beginning of each individual event, and G4Event has a unique event ID. So the general idea is that you’ll open the file at the beginning of the job, and either read all the data into a buffer (what you’re doing now), or read one line at a time within GeneratePrimaries(). Then, you can use the unique event ID in place of your [a] index, to get a different single gamma for each event.

This is code you will need to write yourself. I don’t know if there are examples which deal with reading primaries as inputs from a file.

[janmajay]

Thanks a lot, it seems working. But the only problem is when I am inputting this much huge number of gammas, around 10000000, simulation stops.

[mkelsey]

10 million shouldn’t be a big deal. You can’t do more than 2.14 billion in a job, because G4 event ID is an int not an unsigned. Is it really “stopped,” or is it just taking a long time to finish? Before your /run/beamOn line, you could add a line in your macro file, /run/printProgress 1000. It’ll print a line every 1000 events so you can tell that the job is still moving.

[janmajay]

I did this,
void MyPrimaryGenerator::GeneratePrimaries(G4Event anEvent)
ReadFiles();
G4int eventNumber =anEvent->GetEventID();
G4ThreeVector pos(0., 0., 0.);
//G4ThreeVector mon(0., 0., 1.);
G4ThreeVector mon = G4RandomDirection();
G4int fNEvents = gammasv.size();
//fNEvents = fGammaEnergies.size();
//G4double energy = gammasv[eventNumber];
G4ParticleDefinition particle = G4ParticleTable::GetParticleTable()->FindParticle(“gamma”);
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticleEnergy(gammasv[eventNumber] *keV);
fParticleGun->SetParticlePosition(pos);
fParticleGun->SetParticleMomentumDirection(mon);

fParticleGun->GeneratePrimaryVertex(anEvent);

void MyPrimaryGenerator::ReadFiles()
{

std::vector gammasv;

std::ifstream file("/home/graviton/gamma_cascade/output_file.txt");
 
 std::string line;
 while (std::getline(file, line)) {
    std::istringstream linestream(line);
    double gamma;
    linestream >> gamma;
        gammasv.push_back(gamma);
        
        }

std::ofstream binary_file_stream("gammas.bin", std::ios::binary);
binary_file_stream.write((char *)gammasv.data(), gammasv.size() * sizeof(double));
binary_file_stream.close();

std::ifstream binary_file("gammas.bin", std::ios::binary);


// Get the size of the file
binary_file.seekg(0, binary_file.end);
int size = binary_file.tellg();
binary_file.seekg(0, binary_file.beg);

// Read the vector of gammas from the binary file
//std::vector<double> gammasv(size / sizeof(double));
binary_file.read((char *)gammasv.data(), size);
binary_file.close();

and if I do even /run/beamOn 100. it takes enormous time and then it hanges , and then get killed .

[mkelsey]

Ah, you’re reopening and reading the file on every event; that’s what’s taking the huge amount of time.

Are you running this with multiple threads, or with just a single thread? If you’re using just a single thread, then you know that eventID==0 is the start of the run, so you can write,

G4int eventNumber = anEvent->GetEventID();
if (eventNumber==0) ReadFiles();

using your old, original version of ReadFiles() where you fill the array of gamma energies.

If you’re running multi-threaded, it gets a bit trickier, because you’ll need to create a separate global singleton class to handle the file, and provide a const accessor function so the different threads can ask the singleton for each gamma energy.

[janmajay]

It worked. Thanks a lot. Thanks for dealing with my ignorance with valuable patience.