---

interact.cxx

Go to the documentation of this file.

// Test of interactions in Gismo
// $Header: /nfs/slac/g/glast/ground/cvs/heavy/test/interact.cxx,v 1.1 2000/04/05 22:16:29 burnett Exp $

#include "gismo/Medium.h"
#include "gismo/Material.h"
#include "gismo/MCParticle.h"
#include "gismo/Interactor.h"
#include "gismo/Units.h"
#include "heavy/HeavyInteractor.h"
#include "gismo/GheishaInteractor.h"


  const char* particleName="He";
  const char* materialName="Al";
  float pBeam =10.0;
  unsigned eventCount=1;

using namespace std;

Vector pinit(0,0,pBeam);

void WARNING(const char * text){cerr << "WARNING: " << text << '\n';}
void FATAL(const char * text){
    cerr << "FATAL error: " <<text << '\n'; 
    cerr << "type help for instructions\n";
    exit(-1);}

const char* help=
"Help for the heavy test program interact:\n\n"
"   heavy_test [[[[particleName] materialName] eventCount] pBeam]\n\n"
"where:\n"
"   particleName [default He] is the name of a particle to interact:\n"
"                 (\"decay list\" for a list)\n"
"   materialName [Al] the name of a PEGS4 material file, e.g. Al.mat\n"
"   eventCount [1] number of events to generate\n"
"   pBeam      [.0 GeV] the beam momentum\n";


int main(int argc, char* argv[])
{
    int n=0;
    if( argc>++n ) particleName = argv[n];
    if( strcmp(particleName, "help")==0 ) {cout << help; return 0;}
    if( argc>++n ) materialName = argv[n];
    if( argc>++n ) eventCount = atoi(argv[n]);
    if( argc>++n ) pBeam = atof(argv[n]);

    // this sets up interactors for the default PDT
    MCParticle::addInteractor( new HeavyInteractor );
    MCParticle::addInteractor( new GheishaInteractor("hadronic") );


    // set up medium (really only for the material) and particle
    Medium* med = new Medium((Medium*)0,(Shape*)0,materialName);
    MCParticle p(particleName);
    p.setMomentum(Vector(0,0,pBeam));

    cout << "\nInteraction test of "
      << pBeam << " GeV " << particleName
         << " in " << materialName << '\n';
    cout << "Interactions performed by: "<< p.interactor()->className()
          << ", Interaction length: "
          << p.interactionLength(med) << " cm\n";

 

    for( unsigned i=0; i<eventCount; i++) {
                MCParticle q(p);
            if( q.interact(med))
               q.setStatus(MCParticle::INTERACTED);
               cout << q << '\n';
        }
    cout.flush();
    return 0;
}
                                  

---