2 circular targets

Getting Started
64

Hi @anna thank you. Then, do you confirm that I’m accurately collecting:

  1. The released energy inside the 2 targets
  2. The PDGID of the particle that is releasing the energy
  3. The kinetic energy of the particle that is releasing the energy
  4. The directions of the particle that is releasing the energy?

Anyway, I sent the code to my supervisor by e.mail and he said that in this way I just have a line for each event, but the primary particle could intercat so I need more lines for each event…do you know what I should do?

65

Someone can help for this please?

66

After filling all the fields inside an ntuple you need to also do analysisManager->AddNtupleRow() (so not only in the event action but also inside of stepping action).
However, you are always filling in only part of the ntuple (either accumulated energies which will have always 1 entry per event, or only volume1, or only volume2), so you can consider creating 3 ntuples for each of those cases and call ::AddNtupleRow(<ntupleID>). Otherwise the un-filled columns will yield 0 entries.
You can read more in the user guide and have a look how to create more than one ntuple in the extended analysis/AnaEx01 example.

67

Thank you @anna I’m not sure to understand what you mean…anyway

  1. In the B1RunAction.cc
    a. I wrote
    analysisManager->CreateNtuple("B1Volume1", "Volume 1");
    analysisManager->CreateNtupleDColumn("ParticleID");
    analysisManager->CreateNtupleDColumn("kinEnergy");
    analysisManager->CreateNtupleDColumn("MomentumDirection.x");
    analysisManager->CreateNtupleDColumn("MomentumDirection.y");
    analysisManager->CreateNtupleDColumn("MomentumDirection.z");
    analysisManager->CreateNtupleDColumn("OutputEnergy");
    analysisManager->FinishNtuple();
    //ntupla id=1
    analysisManager->CreateNtuple("B1Volume2", "Volume 2");
    analysisManager->CreateNtupleDColumn("ParticleID2");
    analysisManager->CreateNtupleDColumn("kinEnergy2");
    analysisManager->CreateNtupleDColumn("MomentumDirection.x2");
    analysisManager->CreateNtupleDColumn("MomentumDirection.y2");
    analysisManager->CreateNtupleDColumn("MomentumDirection.z2");
    analysisManager->CreateNtupleDColumn("OutputEnergy2");
    analysisManager->FinishNtuple();

to create 2 Ntuples one for the volume 1 and one for the volume 2.

B1RunAction.cc (8.4 KB)

  1. In 1 Stepping action I wrote
    a. if (volume == fScoringVolume) {
    int pdg = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
    analysisManager->FillNtupleDColumn(0,0, pdg);
    double kinEnergy = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
    analysisManager->FillNtupleDColumn(0,1, kinEnergy);
    double MomDirx = step->GetTrack()->GetMomentumDirection().x();
    analysisManager->FillNtupleDColumn(0,2, MomDirx);
    double MomDiry = step->GetTrack()->GetMomentumDirection().y();
    analysisManager->FillNtupleDColumn(0,3, MomDiry);
    double MomDirz = step->GetTrack()->GetMomentumDirection().z();
    analysisManager->FillNtupleDColumn(0,4, MomDirz);
    analysisManager->AddNtupleRow();

}
if (volume == fScoringVolume2) {
int pdg2 = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
analysisManager->FillNtupleDColumn(1,0, pdg2);
double kinEnergy2 = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
analysisManager->FillNtupleDColumn(1,1, kinEnergy2);
double MomDirx2 = step->GetTrack()->GetMomentumDirection().x();
analysisManager->FillNtupleDColumn(1,2, MomDirx2);
double MomDiry2 = step->GetTrack()->GetMomentumDirection().y();
analysisManager->FillNtupleDColumn(1,3, MomDiry2);
double MomDirz2 = step->GetTrack()->GetMomentumDirection().z();
analysisManager->FillNtupleDColumn(1,4, MomDirz2);
analysisManager->AddNtupleRow();
}
B1SteppingAction.cc (6.6 KB)

  1. In B1EventAction I wrote:
    a. // fill ntupla 0
    analysisManager->FillNtupleDColumn(0,5, fEdep);
    analysisManager->AddNtupleRow();
    // fill ntupla 1
    analysisManager->FillNtupleDColumn(1,5, fEdep2Tar);
    analysisManager->AddNtupleRow();

B1EventAction.cc (3.8 KB)
So I compiled and runned, it looks like without errors

log.txt (98.0 KB)

Then I opened the ROOT file it looks like i’ve two ntuple
but you see…I get 24 events in the B1Volume1 and 0 event in the B1Volume2

image
image732×1007 29.7 KB

Moreover, you wrote to create ntuples for each case…what else ntupla should I create?

68

AddNtupleRow is missing an argument, you only fill the first ntuple.
By third ntuple I meant e.g. overall event statistics, filled only once per event (the one from event action).

69

Thank you @anna. I forgot that in the second AddNtupleRow I had to add the ntupla id!
So

  1. In the B1EventAction
    a. I sobstitued
    // fill ntupla 1
    analysisManager->FillNtupleDColumn(1,5, fEdep2Tar);
    analysisManager->AddNtupleRow();

by
// fill ntupla 1
analysisManager->FillNtupleDColumn(1,5, fEdep2Tar);
analysisManager->AddNtupleRow(1);

B1EventAction.cc (3.8 KB)

  1. In B1SteppingAction.
    a. I substituted
    if (volume == fScoringVolume2) {
    int pdg2 = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
    analysisManager->FillNtupleDColumn(1,0, pdg2);
    double kinEnergy2 = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
    analysisManager->FillNtupleDColumn(1,1, kinEnergy2);
    double MomDirx2 = step->GetTrack()->GetMomentumDirection().x();
    analysisManager->FillNtupleDColumn(1,2, MomDirx2);
    double MomDiry2 = step->GetTrack()->GetMomentumDirection().y();
    analysisManager->FillNtupleDColumn(1,3, MomDiry2);
    double MomDirz2 = step->GetTrack()->GetMomentumDirection().z();
    analysisManager->FillNtupleDColumn(1,4, MomDirz2);
    analysisManager->AddNtupleRow();
    }

by

if (volume == fScoringVolume2) {
int pdg2 = step->GetTrack()->GetParticleDefinition()->GetPDGEncoding();
analysisManager->FillNtupleDColumn(1,0, pdg2);
double kinEnergy2 = step->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
analysisManager->FillNtupleDColumn(1,1, kinEnergy2);
double MomDirx2 = step->GetTrack()->GetMomentumDirection().x();
analysisManager->FillNtupleDColumn(1,2, MomDirx2);
double MomDiry2 = step->GetTrack()->GetMomentumDirection().y();
analysisManager->FillNtupleDColumn(1,3, MomDiry2);
double MomDirz2 = step->GetTrack()->GetMomentumDirection().z();
analysisManager->FillNtupleDColumn(1,4, MomDirz2);
analysisManager->AddNtupleRow(1);
}

B1SteppingAction.cc (6.6 KB)

The I compiled and runned without errors
log.txt (97.9 KB)

and I opened the ROOT file and now I get 2 TTRee having 12 events each one.

It looks like I get:

  1. Same particle ID in the 2 targets

    image
    image860×617 28.4 KB

    image
    image865×617 28.7 KB

  2. Same kinetic energy in the 2 targets

    image
    image973×686 38.6 KB

    image
    image1048×751 49.7 KB

  3. Different released energy in the 2 targets

    image
    image855×620 32.8 KB

    image
    image858×605 32 KB

  4. Different momentum direction along x in the 2 targets

    image
    image856×622 33 KB

    image
    image850×617 31.7 KB

  5. Different momentum direction along y in the 2 targets

image
image865×625 39.8 KB

image
image867×628 39.2 KB

  1. Different momentum direction along z in the 2 targets
    image
    image863×614 33.1 KB

    image
    image857×617 32.9 KB

Do you think that are the results right to get
a. Particle ID
b. Kinetic energy
c. Direction
d. Released energy in the two targets
about the particle? And are the results also about secondary particle or just about primary particles?

I simulated a e+ beam, and by the graph of particleid in the two targets I see there are 6 e+, 5 photons and 1 electron…so I guess that the photons and the electron are the secondary particles…
do you think that results are right to get information about primary and secondary particle in the targets?? I’m not sure, because you wrote to create ntuples for each case…what else ntupla should I create?

70

As I wrote above, each case = either target1, or target 2, or event information (sum of energy in targets, etc). This would mean you need another ntuple for the information you fill in the event action, as there will be only one entry per event (with sum of energy in target 1 and sum of energy in target 2).

71

To see what happens in your simulation, I suggest to study the output after adding /tracking/verbose 1 to your .mac macro. Initial particle will have parentID=0 and trackID=1. Next, all the particles will have an incremented trackID and parentID pointing to the parent. You will see what and where is produced.

72

Thank you @anna

To store the sum of energy

  1. In B1SteppingAction.cc
    a. I wrote

    if (volume == fScoringVolume) {
    G4double edepStep = step->GetTotalEnergyDeposit();
    fEventAction->AddEdep(edepStep);
    return (edepStep);
    }
    if (volume == fScoringVolume2) {
    G4double edepStep2 = step->GetTotalEnergyDeposit();
    fEventAction->AddEdep2(edepStep2);
    return (edepStep2);
    }
    if (volume == fScoringVolume || volume == fScoringVolume2) {
    G4double edepStep2tartot = edepStep+edepStep2;
    fEventAction->AddEdeptot(edepStep2tartot);
    }

B1SteppingAction.cc (6.8 KB)

  1. In B1EventAction.cc
    a. I added fEdep2TarSum = 0.;
    b. I added fRunAction->AddEdeptot(fEdep);
    c. I added analysisManager->FillNtupleDColumn(2,0, fEdep2TarSum);
    analysisManager->AddNtupleRow(2);

B1EventAction.cc (4.3 KB)

  1. In B1RunAction.cc
    a. I added
    analysisManager->CreateNtuple("B1VolumeTot", "Volume Tot");
    analysisManager->CreateNtupleDColumn("OutputEnergyTot");
    analysisManager->FinishNtuple();
    b. I added
    void B1RunAction::AddEdeptot(G4double edep)
    {
    fEdep += edep;
    }

B1RunAction.cc (8.8 KB)

  1. In B1EventAction.hh
    a. I added void void AddEdeptot(G4double edep) { fEdep2TarSum += edep; }
    b. I added G4double fEdep2TarSum;

B1EventAction.hh (3.0 KB)

  1. In B1RunAction.hh
    a. I added void AddEdeptot (G4double edep);

B1RunAction.hh (2.9 KB)

But compiling I get these errors

image
image1105×650 28.6 KB

I understand they are due to

return (edepStep); and return (edepStep2);…

To watch what is happening

In my macro file I had the /tracking/verbose 1 line yet…what should I do to see what is happening?

run1.mac.txt (559 Bytes)

73

I am not sure why you introduced those changes. Do you need a variable that returns the sum of first and second target? What I meant is that the ntuple you fill in event action could be declared as a third ntuple, with Id=2 - which means the columns " OutputEnergy" and " OutputEnergy2" will be created within analysisManager->CreateNtuple("EventVariables", "Event");.
And yes, the errors mean you try to return a value in a method that is not supposed to return anything. Also, return means exit from the method, so your third “if” will never happen as one of the first 2 will.

Regarding the tracking verbosity, you need to analyse the output that you obtain. It should start with something like that:

G4WT0 > *********************************************************************************************************
G4WT0 > * G4Track Information:   Particle = e+,   Track ID = 1,   Parent ID = 0
G4WT0 > *********************************************************************************************************
G4WT0 > 
G4WT0 > Step#    X(mm)    Y(mm)    Z(mm) KinE(MeV)  dE(MeV) StepLeng TrackLeng  NextVolume ProcName
G4WT0 >     0        0        0   -1e+03   4.5e+04        0        0         0       World initStep
G4WT0 >     1        0        0     -609   4.5e+04   0.0805      391       391       World eIoni
74

Thank you @anna. No I don’t think to need the sum of the released energies in the two targets (my supervisor didn’t ask me it).Reading yuor previous message

As I wrote above, each case = either target1, or target 2, or event information (sum of energy in targets, etc).

I understood that you were suggesting me to do it!

Well… the if now I understand what you mean, I’ve to create a third ntupla whitin I’ve a column for the released energy in target 1 and a second column for the released energy in target 2…true?
Then

  1. In B1RunAction.cc
    a. I added
    analysisManager->CreateNtuple("B1EventVariables", "Events");
    analysisManager->CreateNtupleDColumn("ReleasedEnergy");
    analysisManager->CreateNtupleDColumn("ReleasedEnergy2");
    analysisManager->FinishNtuple();

B1RunAction.cc (8.6 KB)

  1. In B1EventAction.cc
    a. I added
    analysisManager->FillNtupleDColumn(2,0, fEdep);
    analysisManager->FillNtupleDColumn(2,1, fEdep2Tar);
    analysisManager->AddNtupleRow(2);

B1EventAction.cc (4.0 KB)

Then I compiled without errors and I see there are 3 Ntuples, the first one for the Volume1, the second one for the Volume 2 and the last one that I just created to print both relased energy in first and second target. But I think there is a problem, because, you see in the first 2 Ntuples I get 12 events, instead in the third one I get just 5 events! So I’m not filling all the events ( it’s the same problem that I had before to create a second Ntupla!)

image
image715×596 18.8 KB

image
image714×599 18.3 KB

image
image718×292 8.04 KB

Regarding the tracking verbosity, yes when I run the simulation, I get the table that you showed me. You can see it in this file
log.txt (60.7 KB)

I read the are Events having Track ID >1 and Parent ID>0, but I’m seeing that in this output I’ve 68 tables, instead in the ROOT file I’ve just 12 events in the Volume 1 and 12 event in the Volume 2…so maybe I’m not storing all the events in the Root file?

Thank you

75

You save only particles that traverse your detector. The last but one column says which volume the step enters. Check this output to see how many particles (and which) create entry in your ntuple. The row with your volume will be exactly the step that you store in the ntuple (momentum direction, energy, PDG…)

76

ah ok @anna, then I’ve just 12 events, just because the others don’t pass go through the 2 targets, but out of them, true?

Instead, what should I do to fix the problem about the third ntupla?

77

If that’s how many you counted then yes.

It’s not a problem with the 3rd ntuple, but you may create it (in run action), and create OutputEnergy and OutputEnergy columns inside it, and fill this ntuple (with probably ID=2) in the event action. Just as you did the other two ntuples.

78

@anna, could be possible that I’m storing more events than how many events I’ve?

My two targets are defined

G4Tubs* solidEnv =
new G4Tubs("Envelope", //its name
pRMin, pRMax, pDz, pSPhi, pDPhi); //its size*/

    and

G4Tubs* solidEnv2 =
new G4Tubs("Envelope2", //its name
pRMin, pRMax, pDz, pSPhi, pDPhi); //its size*/

then, their names are " Envelope" and “Envelope2” in the log file I read the rows, for the target 1

  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

15 -0.0131 0.00543 -1.5 4.5e+04 0.0312 146 999 Envelope Transportation

  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

11 -0.00536 -0.00191 -1.5 4.5e+04 0.000176 0.639 999 Envelope Transportation

  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0
8 -0.00152   -0.004     -1.5   4.5e+04   0.0152     73.6       999    Envelope Transportation
  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

11 -0.00126 0.00137 -1.5 4.5e+04 0.0264 140 999 Envelope Transportation
12 -0.00128 0.0014 1.28 4.5e+04 0.595 2.78 1e+03 Envelope eIoni

  • G4Track Information: Particle = e-, Track ID = 12, Parent ID = 1
0 -0.00128   0.0014     1.28      1.65        0        0         0    Envelope initStep
  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

12 -0.00435 0.00198 -1.5 4.5e+04 0.00186 15.2 999 Envelope Transportation

and, for the target 2

  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

17 -0.0139 0.00571 21.5 4.5e+04 0.00374 20 1.02e+03 Envelope2 Transportation

  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

14 -0.00571 -0.00213 21.5 4.5e+04 0.00198 18.7 1.02e+03 Envelope2 Transportation
15 -0.00573 -0.00214 22.3 4.5e+04 0.159 0.774 1.02e+03 Envelope2 eIoni

  • G4Track Information: Particle = e-, Track ID = 13, Parent ID = 1
0 -0.00573 -0.00214     22.3       5.9        0        0         0   Envelope2 initStep
  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

10 -0.00128 -0.00456 21.5 4.5e+04 0.00537 20 1.02e+03 Envelope2 Transportation

  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

14 -0.00204 0.00215 21.5 4.5e+04 0.00683 20 1.02e+03 Envelope2 Transportation

  • G4Track Information: Particle = e+, Track ID = 1, Parent ID = 0

15 -0.0039 0.00171 21.5 4.5e+04 0.00222 17.9 1.02e+03 Envelope2 Transportation

Then, I read 7 steps in the target 1 and 7 steps in the target 2…but in the Volume 1 and Volume 2 of the ROOT file I’ve 12 events.…

Regarding of the third Ntupla, I create it in RunAction by the code

analysisManager->CreateNtuple("B1EventVariables", "Events");
analysisManager->CreateNtupleDColumn("ReleasedEnergy");
analysisManager->CreateNtupleDColumn("ReleasedEnergy2");
`analysisManager->FinishNtuple();

and I fill it in EventAction by

analysisManager->FillNtupleDColumn(2,0, fEdep);
analysisManager->FillNtupleDColumn(2,1, fEdep2Tar);
alysisManager->AddNtupleRow(2);

but it doesn’t fill 12 events (as in Volume1 and Volume2), but just 5 events…

79

5 events is probably correct. Isn’t that what you wrote in your macro (after /run/beamOn) ?
What you store in the other ntuples is the number of steps that particles make inside the (sensitive) volumes.

80

Yes @anna in the run1.mac I wrote /run/beamOn 5 , so do you think that it’s correct to get 12 events in the first and second ntupla and to get 5 events in the third one?

And what about the difference between the 7 events counting the ling of the tracking verbosity and the 12 events stored in the ROOT file? Why do I have 7 events in the tracking verbosity and 12 in the Volume 1 and Volume 2 ntupla of the ROOT file?
Thank you

81

5 events in the third ntuple is coming from the initial positrons (which we call 5 events).
If you say you see 7 particles within volume, it fits the scenario where you call AddNtupleRow in the stepping action (7 times = once per particle step within volume) AND in the event action (5 times = once per event). 5+7=12. You need to fill those two first ntuples only in the stepping action, not also in the event action.

82

thank you @anna then

  1. In B1Eventaction.cc
    a. I deleted the lines
    analysisManager->FillNtupleDColumn(0,5, fEdep);
    analysisManager->AddNtupleRow();
    analysisManager->FillNtupleDColumn(1,5, fEdep2Tar);
    analysisManager->AddNtupleRow(1);

B1EventAction.hh (2.7 KB)

  1. In B1SteppingAction.cc
    a. In the if (volume == fScoringVolume) I added the line
    analysisManager->FillNtupleDColumn(0,5, fEdep);
    b. In the if (volume == fScoringVolume2) I added the line
    analysisManager->FillNtupleDColumn(1,5, fEdep2Tar);

B1SteppingAction.hh (2.6 KB)

  1. In B1SteppingAction.hh
    a. I added the lines
    G4double fEdep;
    G4double fEdep2Tar;

the I compiled and runned without errors.

log.txt (97.9 KB)
In the tracking verbosity I read there are 7 steps in the Volume 1 and 7 steps in the Volume 2. Then I opened the ROOT file but you see…

image
image713×989 29.5 KB

I still get 12 events in the first two ntuples, instead I’ve 7 steps in the tracking verbosity…

Then in B1EventAction.cc

  1. I also deleted the lines
    a. fRunAction->AddEdep(fEdep);
    fRunAction->AddEdep2(fEdep);

B1EventAction.hh (2.7 KB)

I compiled and runned again without erros

log.txt (97.9 KB)

and, in the tracking verbosity, I still get 7 steps in Volume 1 and 7 steps in Volume 2. Then I opened the ROOT file and now I get just 7 events in the ROOT file too!

image
image759×849 41.1 KB

But….given that to get 7 events in the ROOT file I had to delete the lines

fRunAction->AddEdep(fEdep);
fRunAction->AddEdep2(fEdep);

is it a problem for the simulation? Will I miss something when I will start the simulations? Near this line there was this comment

accumulate statistics in run action

Thank you

83

You didn’t use the run action, it was a leftover from the example.