/* * hddmInput - functions to handle Monte Carlo generator input to HDGeant * through the standard hddm i/o mechanism. * * Interface: * openInput(filename) - open input stream to file * skipInput(count) - skip next events on open input file * nextInput() - advance to next event on open input stream * loadInput() - push current input event to Geant kine structures * storeInput() - pop current input event from Geant kine structures * closeInput() - close currently open input stream * * Richard Jones * University of Connecticut * July 13, 2001 * * Usage Notes: * 1) Most Monte Carlo generators do not care where the vertex is placed * inside the target, and specify only the final-state particles' * momenta. In this case the vertex position has to be randomized by * the simulation within the beam/target overlap volume. If the vertex * position from the generator is (0,0,0) then the simulation vertex is * generated uniformly inside the cylinder specified by TARGET_LENGTH, * BEAM_DIAMETER, and TARGET_CENTER defined below. * * Revision history: * > Apr 10, 2006 - David Lawrence * Added comments to explain a little what each of these routines is doing. * No functional changes. * * > Dec 15, 2004 - Richard Jones * Changed former behaviour of simulation to overwrite the vertex * coordinates from the input record, if the simulation decides to * override the input values. At present this happens whenever the * input record specifies 0,0,0 for the vertex, but in the future it * may be decided to let the simulator determine the vertex position * in other cases. Since it is not part of the simulation proper, the * decision was made to store this information in the reaction tag. */ #define TARGET_LENGTH 30 #define BEAM_DIAMETER 0.5 #define TARGET_CENTER 65 #include #include #include #include #include s_iostream_t* thisInputStream = 0; s_HDDM_t* thisInputEvent = 0; /*------------------------- /* openInput /*-------------------------*/ int openInput (char* filename) { /* Open HDDM file for reading in "thrown" particle kinematics */ thisInputStream = open_s_HDDM(filename); return (thisInputStream == 0); } /*------------------------- /* skipInput /*-------------------------*/ int skipInput (int count) { /* Skip over the specified number of events in the input file */ /* NOTE: This is NOT safe for cross-platform use!! */ int* buff = (int*) malloc(1000000); while (thisInputStream && (count > 0)) { int ret = fread(buff, sizeof(int), 1, thisInputStream->fd); if (ret == 1) { int nw = *buff; ret = fread(buff, sizeof(int), nw, thisInputStream->fd); --count; } else { break; } } free(buff); return count; } /*------------------------- /* nextInput /*-------------------------*/ int nextInput () { /* Read in the next HDDM event. This only reads it into the HDDM buffer "thisInputEvent" and does not yet define the particles to GEANT. See loadInput for that. */ if (thisInputStream == 0) { return 9; /* input stream was never opened */ } else if (thisInputEvent) { flush_s_HDDM(thisInputEvent, 0); } thisInputEvent = read_s_HDDM(thisInputStream); return (thisInputEvent == 0); } /*------------------------- /* loadInput /*-------------------------*/ int loadInput () { /* Extracts the "thrown" particle 4-vectors and types from the current HDDM buffer "thisInputEvent" and creates a vertex for them(gsvert) and defines the GEANT (gskine) for tracking. */ s_Reactions_t* reacts; int reactCount, ir; reacts = thisInputEvent ->physicsEvents->in[0].reactions; reactCount = reacts->mult; for (ir = 0; ir < reactCount; ir++) { s_Vertices_t* verts; int vertCount, iv; s_Reaction_t* react = &reacts->in[ir]; verts = react->vertices; vertCount = verts->mult; for (iv = 0; iv < vertCount; iv++) { float v[3]; int ntbeam = 0; int nttarg = 0; int nubuf = 0; float ubuf; int nvtx; s_Products_t* prods; int prodCount, ip; s_Vertex_t* vert = &verts->in[iv]; v[0] = vert->origin->vx; v[1] = vert->origin->vy; v[2] = vert->origin->vz; if ((v[0] == 0) && (v[1] == 0) && (v[2] == 0)) { v[0] = 1; v[1] = 1; v[2] = TARGET_CENTER; while (v[0]*v[0] + v[1]*v[1] > 0.25) { int len = 3; grndm_(v,&len); v[0] -= 0.5; v[1] -= 0.5; v[2] -= 0.5; } v[0] *= BEAM_DIAMETER; v[1] *= BEAM_DIAMETER; v[2] *= TARGET_LENGTH; v[2] += TARGET_CENTER; vert->origin->vx = v[0]; vert->origin->vy = v[1]; vert->origin->vz = v[2]; } gsvert_(v, &ntbeam, &nttarg, &ubuf, &nubuf, &nvtx); prods = vert->products; prodCount = prods->mult; for (ip = 0; ip < prodCount; ip++) { int ntrk; float p[3]; Particle_t kind; s_Product_t* prod = &prods->in[ip]; kind = prod->type; p[0] = prod->momentum->px; p[1] = prod->momentum->py; p[2] = prod->momentum->pz; if (prod->decayVertex == 0) { gskine_(p, &kind, &nvtx, &ubuf, &nubuf, &ntrk); } } } } return 0; } /*------------------------- /* storeInput /*-------------------------*/ int storeInput (int runNo, int eventNo, int ntracks) { /* This is called by the built-in generators (coherent brem. and single track) in order to store the "thrown" particle parameters in the output HDDM file. What this actually does is free the input buffer "thisInputEvent" if it exists and creates a new one. When an external generator is used, the thisInputEvent buffer is kept unmodified and this routine is never called.*/ s_PhysicsEvents_t* pes; s_Reactions_t* rs; s_Vertices_t* vs; s_Origin_t* or; s_Products_t* ps; int nvtx, ntbeam, nttarg, itra, nubuf; float vert[3], plab[3], tofg, ubuf[10]; Particle_t kind; if (thisInputEvent) { flush_s_HDDM(thisInputEvent, 0); } thisInputEvent = make_s_HDDM(); thisInputEvent->physicsEvents = pes = make_s_PhysicsEvents(1); pes->in[0].reactions = rs = make_s_Reactions(1); pes->mult = 1; rs->mult = 1; rs->in[0].vertices = vs = make_s_Vertices(99); vs->mult = 0; for (itra = 1; itra <= ntracks; itra++) { gfkine_(&itra,vert,plab,&kind,&nvtx,ubuf,&nubuf); if (nvtx < 1) { return 1; } else { vs->mult = (nvtx < vs->mult)? vs->mult : nvtx; } gfvert_(&nvtx,vert,&ntbeam,&nttarg,&tofg,ubuf,&nubuf); or = vs->in[nvtx-1].origin; ps = vs->in[nvtx-1].products; if (or == HDDM_NULL) { or = make_s_Origin(); vs->in[nvtx-1].origin = or; or->vx = vert[0]; or->vy = vert[1]; or->vz = vert[2]; or->t = tofg; } if (ps == HDDM_NULL) { ps = make_s_Products(ntracks); vs->in[nvtx-1].products = ps; ps->mult = 0; } ps->in[ps->mult].type = kind; ps->in[ps->mult].momentum = make_s_Momentum(); ps->in[ps->mult].momentum->px = plab[0]; ps->in[ps->mult].momentum->py = plab[1]; ps->in[ps->mult].momentum->pz = plab[2]; ps->in[ps->mult].momentum->E = plab[3]; ps->mult++; } pes->in[0].runNo = runNo; pes->in[0].eventNo = eventNo; return 0; } /*------------------------- /* closeInput /*-------------------------*/ int closeInput () { /* Close the HDDM input file */ if (thisInputStream) { close_s_HDDM(thisInputStream); thisInputStream = 0; } return 0; } /* entry points from Fortran */ int openinput_ (char* filename) { int retcode = openInput(strtok(filename," ")); return retcode; } int skipinput_ (int* count) { return skipInput(*count); } int nextinput_ () { return nextInput(); } int loadinput_ () { return loadInput(); } int storeinput_ (int* runNo, int* eventNo, int* ntracks) { return storeInput(*runNo,*eventNo,*ntracks); } int closeinput_ () { return closeInput(); }