7 #include "GaudiKernel/Bootstrap.h"
8 #include "GaudiKernel/ISvcLocator.h"
9 #include "GaudiKernel/IMessageSvc.h"
10 #include "GaudiKernel/IDataProviderSvc.h"
11 #include "GaudiKernel/INTupleSvc.h"
12 #include "GaudiKernel/NTuple.h"
13 #include "GaudiKernel/SmartDataPtr.h"
31 for (
int i : pdgids ){
38 NTupleFilePtr file1(
ntupleSvc(),
"/NTUPLES/FILE1");
40 SmartDataPtr<NTuple::Directory>
41 ntdir(
ntupleSvc(),
"/NTUPLES/FILE1/StepNtuple");
43 if ( !ntdir ) ntdir =
ntupleSvc()->createDirectory(file1,
"StepNtuple");
45 NTuplePtr
nt(
ntupleSvc(),
"/NTUPLES/FILE1/StepNtuple/10");
47 nt =
ntupleSvc()->book (ntdir.ptr(), 10,CLID_ColumnWiseTuple,
"GEANT4 Step NTuple");
52 if(
nt->addItem (
"nstep",
m_nsteps,0 ,50000)!=StatusCode::SUCCESS
77 ||
nt->addItem (
"evtid",
m_evtid)!=StatusCode::SUCCESS)
101 if(!
ntupleSvc()->writeRecord(
"/NTUPLES/FILE1/StepNtuple/10").isSuccess()) {
111 int pdg = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
119 std::cout<<
"ACH139: SG_StepNtuple: other code thinks "<<pdg<<
" is an Rhadron!"<<std::endl;
127 G4Material *
mat = aStep->GetTrack()->GetMaterial();
128 double minA=1500000.;
129 for (
unsigned int i=0;
i<
mat->GetNumberOfElements();++
i){
130 if (
mat->GetElement(
i) && minA>
mat->GetElement(
i)->GetN()){
131 minA=
mat->GetElement(
i)->GetN();
141 if (firstslow || aStep->GetTrack()->GetCurrentStepNumber()<=1 || aStep->GetPostStepPoint()->GetKineticEnergy()==0.){
145 int id = std::abs(aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetPDGEncoding());
154 if (aStep->GetPreStepPoint()->GetGlobalTime()==0)
m_rhadronIndex++;
157 m_pdg[
m_nsteps]=aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
160 m_mass[
m_nsteps]=aStep->GetTrack()->GetDefinition()->GetPDGMass();
162 m_t1[
m_nsteps]=aStep->GetPreStepPoint()->GetGlobalTime();
163 m_t2[
m_nsteps]=aStep->GetPostStepPoint()->GetGlobalTime();
164 G4ThreeVector pos1=aStep->GetPreStepPoint()->GetPosition();
168 G4ThreeVector pos2=aStep->GetPostStepPoint()->GetPosition();
172 m_ke1[
m_nsteps]=aStep->GetPreStepPoint()->GetKineticEnergy();
173 m_ke2[
m_nsteps]=aStep->GetPostStepPoint()->GetKineticEnergy();
174 m_pt1[
m_nsteps]=aStep->GetPreStepPoint()->GetMomentum().perp();
175 m_pt2[
m_nsteps]=aStep->GetPostStepPoint()->GetMomentum().perp();
179 m_v2[
m_nsteps]=aStep->GetPostStepPoint()->GetVelocity();
192 if (std::abs(pdg)>1000000 && std::abs(pdg)<10000000){
193 std::cout<<
"ACH129: SG_StepNtuple: Killing non-rh track with pdg "<<pdg<<std::endl;
195 aStep->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
196 const G4TrackVector *tv = aStep->GetSecondary();
198 for (
unsigned int i=0;
i<tv->size();
i++){
199 G4Track *
t = (*tv)[
i];
200 t->SetTrackStatus(fKillTrackAndSecondaries);
210 if (
id==1000021 ||
id==1000005 ||
id==1000006 ||
id==1000512 ||
id==1000522 ||
id==1000991 ||
id==1000993 ||
211 id==1000612 ||
id==1000622 ||
id==1000632 ||
id==1000642 ||
id==1000652 ||
id==1005211 ||
212 id==1006113 ||
id==1006211 ||
id==1006213 ||
id==1006223 ||
id==1006311 ||
213 id==1006313 ||
id==1006321 ||
id==1006323 ||
id==1006333 ||
214 id==1009111 ||
id==1009113 ||
id==1009211 ||
id==1009213 ||
id==1009311 ||
215 id==1009313 ||
id==1009321 ||
id==1009323 ||
id==1009223 ||
id==1009333 ||
216 id==1092112 ||
id==1091114 ||
id==1092114 ||
id==1092212 ||
id==1092214 ||
id==1092224 ||
217 id==1093114 ||
id==1093122 ||
id==1093214 ||
id==1093224 ||
id==1093314 ||
id==1093324 ||
id==1093334)