Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove the St_geant_Maker dependencies in the event generator to starsim bridge #616

Merged
merged 2 commits into from
Nov 2, 2023
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
90 changes: 75 additions & 15 deletions StRoot/StarGenerator/BASE/AgStarReader.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,12 @@ ClassImp(AgStarReader);
#include "TDatabasePDG.h"
#include "TParticlePDG.h"

#include "St_geant_Maker/St_geant_Maker.h"
#include "TLorentzVector.h"
#include <vector>
#include <map>
using namespace std;

#include "TGiant3.h"
#include <cassert>

#include "StarGenerator/UTIL/StarParticleData.h"
#include "StarGenerator/BASE/StarParticleStack.h"
Expand All @@ -21,18 +20,69 @@ AgStarReader *AgStarReader::mInstance = 0;
// TODO: Refactor TDatabasePDG to StarParticleData

//
// Setup the interface with starsim
// Setup the interface with starsim and direct to geant3
//
#define ageventread F77_NAME(ageventread,AGEVENTREAD)
#define agsvert F77_NAME(agsvert, AGSVERT)
#define agskine F77_NAME(agskine, AGSKINE)
#define gskine F77_NAME(gskine, GSKINE)
#define gsvert F77_NAME(gsvert, GSVERT)

// Comis routine for obtaining the address of common blocks
#define gcomad F77_NAME(gcomad,GCOMAD)

extern "C" {
void type_of_call ageventread() { AgStarReader::Instance().ReadEvent(); }
void type_of_call agsvert( float *vertex, int *nb, int *nt, float *ubuf, int *nu, int *nv );
void type_of_call agskine( float *plab, int *ip, int *nv, float *ubuf, int *nb, int *nt );
void type_of_call gsvert( float *, int &, int &, float *, int &, int &);
void type_of_call gskine( float *, int &, int &, float *, int &, int &);
void type_of_call gcomad(DEFCHARD, int*& DEFCHARL);
};

//
//----------GCTRAK ... normally I would rather poke my eye out rather than hardcode
// this... but this structure has existed w/in GEANT3 unchanged for 3 decades or so.
// Think we can get away with this here...
//
#define MAXMEC 30
typedef struct {
float vect[7];
float getot;
float gekin;
int vout[7];
int nmec;
int lmec[MAXMEC];
int namec[MAXMEC];
int nstep;
int maxnst;
float destep;
float destel;
float safety;
float sleng;
float step;
float snext;
float sfield;
float tofg;
float gekrat;
float upwght;
int ignext;
int inwvol;
int istop;
int igauto;
int iekbin;
int ilosl;
int imull;
int ingoto;
int nldown;
int nlevin;
int nlsav;
int istory;
} Gctrak_t;

Gctrak_t* fGctrak = 0;


// ----------------------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------------------
Expand All @@ -45,7 +95,10 @@ AgStarReader::AgStarReader() : TObject(), mStack(0), mParticleData(&StarParticle
// ----------------------------------------------------------------------------------------------------
AgStarReader &AgStarReader::Instance()
{
if ( !mInstance ) mInstance = new AgStarReader();
if ( !mInstance ) {
mInstance = new AgStarReader();
gcomad(PASSCHARD("GCTRAK"),(int*&) fGctrak PASSCHARL("GCTRAK"));
}
return (*mInstance);
}
// ----------------------------------------------------------------------------------------------------
Expand All @@ -54,29 +107,27 @@ AgStarReader &AgStarReader::Instance()
void AgStarReader::ReadEvent()
{

TGiant3 *geant3 = St_geant_Maker::instance()->Geant3();

TParticle *part = 0; // particle
int itrk = -1; // track number


struct Vertex_t {
Double_t x, y, z, t;
double x, y, z, t;
int idx;
};

int idvtx = 0;
map<int, int> start_vtx;

float* ubuf = 0;
int nwbuf=0;

while( (part=mStack->PopNextTrack(itrk)) )
{

// First check that the status of the particle is stable. If not,
// continue on to the next particle.
if ( part->GetStatusCode() != 1 ) continue;

// part->Print();

// Get the parent particle and lookup the vertex id
int parent = part->GetFirstMother();
int myvtx = start_vtx[parent];
Expand All @@ -91,11 +142,17 @@ void AgStarReader::ReadEvent()
float(part->Vz())
};

myvtx = geant3->Gsvert( v, 0, 0 );
assert(myvtx==idvtx);
// Set the vertex of the track
{
int nt = 0;
int nb = 0;
gsvert(v, nt, nb, ubuf, nwbuf, myvtx);
assert(myvtx==idvtx);
}

// Set time of flight
geant3->Gctrak()->tofg = part->T();
assert(fGctrak);
fGctrak->tofg = part->T();

}

Expand All @@ -117,8 +174,11 @@ void AgStarReader::ReadEvent()
}
}


geant3->Gskine( plab, g3id, myvtx );
// And add the particle
{
int nt = 0;
gskine( plab, g3id, myvtx, ubuf, nwbuf, nt );
}

}

Expand Down