---

interact.cxx

Go to the documentation of this file.

// Test of interactions in Gismo
// $Header: /nfs/slac/g/glast/ground/cvs/gismo/src/test/interact.cxx,v 1.1 2000/09/01 20:25:15 burnett Exp $

#include "gismo/Medium.h"
#include "gismo/Material.h"
#include "gismo/MCParticle.h"
#include "gismo/Interactor.h"
#include "gismo/Units.h"
#include "gismo/EGSInteractor.h"
#include "egs/PEGSData.h"
#include "gismo/GheishaInteractor.h"

#ifdef TESTPOL
#  include "egs/EGSpolInteractor.h"
#  include "gismo/PDT.h"
#endif
#ifdef MAKETUPLE
# include "test/Tuple.h"
#endif
#include <stdlib.h>

#ifdef TESTPOL
  char* particleName="polgamma";
  char* materialName="Pb";
  float pBeam =0.04;
  unsigned eventCount=10;
#else
  const char* particleName="gamma";
  const char* materialName="Al";
  float pBeam =2.0;
  unsigned eventCount=1;

#endif
#ifdef _MSC_VER
using namespace std;
#endif
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 gismo test program interact:\n\n"
"   interact [[[[particleName] materialName] eventCount] pBeam]\n\n"
"where:\n"
"   particleName [default gamma] 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      [2.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 EGSInteractor("electromagnetic") );
    MCParticle::addInteractor( new GheishaInteractor("hadronic") );

    Material::addPath("../PEGS4:./");

#ifdef TESTPOL
    MCParticle::addInteractor( new EGSpolInteractor("polarizedEM") );
    MCParticle::thePDT()->setInteraction("polarizedEM");
    MCParticle::thePDT()->addParticle("polgamma",    33, 2, 0,  0.) ;
#endif
    // 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";

    // special stuff to look at e, gamma interaction probabilites
    const PEGSData* pegs = (PEGSData*)med->material().getMatData(0);
    if( string(particleName)== "e+") {
        float ke = Units::toMeV(p.e()-p.mass());
        float prob1= pegs->pbr1(ke),
            prob2 = pegs->pbr2(ke);
        cout << "Positron probabilities to brem, bhabha, annihilate: " 
            << prob1 <<", "<< prob2-prob1<< ", " << 1.0-prob2 << endl;
    }else if( string(particleName)=="gamma") {
        float ke = Units::toMeV(p.e());
        float prob1 = pegs->gbr1(ke), prob2 = pegs->gbr2(ke);
        cout << "photon probabilities to pair, compton, photo: "
            << prob1 <<", "<< prob2-prob1<< ", " << 1.0-prob2 << endl;
    }else if( string(particleName)=="e-"){
        float prob1 = pegs->ebr1(Units::toMeV(p.e()-p.mass()));
        cout << "electron probabilites to brems, moller: "
            << prob1 << ", " << 1.0-prob1 << endl;
    }


#ifdef MAKETUPLE
    Tuple tuple("Polarized gamma study");
    TupleItem e1("e1"),e2("e2"),
        theta1("theta1"),theta2("theta2"),
        phi1("phi1"),phi2("phi2"),
        psi("psi"), open("open_angle");
    tuple.writeHeader(cout);
#endif
    for( unsigned i=0; i<eventCount; i++) {
        MCParticle q(p);
        if( q.interact(med))
            q.setStatus(MCParticle::INTERACTED);
#ifdef MAKETUPLE
        if( q.numChildren()==2) {
            GParticle& a= *q.child(0),
                b= *q.child(1); // secondaries

            e1 = a.e();
            e2 = b.e();
            theta1=a.theta();
            theta2=b.theta();
            phi1 = a.phi();
            phi2 = b.phi();
            psi = phi2-phi1; if( psi<0 ) psi=psi+2*M_PI; 
            open = acos(a.direction()*b.direction());
            cout << tuple;
        }
#else
        cout << q << '\n';
#endif
    }
    cout.flush();
    return 0;
}
                                  

---