---

shower.cxx

Go to the documentation of this file.

// Simple test of showering
//
//  Try to contain a gamma ray in a very large medium
// $Id: shower.cxx,v 1.1 2000/09/01 20:25:16 burnett Exp $
#include "BaseApp.h"
#include "SlabWorld.h"
#include "SmplStat.h"

#include "gismo/Interactor.h"
#include "geometry/Box.h"
#include "gismo/EGSInteractor.h"
#include "gismo/GheishaInteractor.h"

#  include <iomanip>

#ifdef _MSC_VER
  using namespace std;
#endif

const char*     particleName=   "gamma";
const char*     materialName=   "Pb";
unsigned long   numberOfEvents= 10;  // default number of events
float           pBeam =         2.0f; //2.0;  // initial momentem
float           cut   =         0.00001f; //2.0;  // k.e. cutOff
float           thick =         10000.0f;  // thickness of slab

SampleStatistic energyStat;

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

const char* help=
"Help for the test program shower:\n\n"
"   shower [[[[particleName] materialName] eventCount] pBeam]keCutOff]\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 [10] number of events to generate\n"
"   pBeam      [2.0 GeV] the beam momentum\n"
"   keCutOff   [.00001 GeV] the kinetic energy cutoff\n";

class ShowerDetector : public Detector
{
public:
 ShowerDetector():Detector(){}
 void score(MCParticle*);
 void clear(){e_dedx = 0; e_gamma = 0; e_nucl = 0; e_neutr = 0; e_stop=0;}
 static int   evtNum;   // event counter
 static float e_dedx;   // de/dx energy deposited
 static float e_gamma;  // gamma energy (interacted or stopped)
 static float e_nucl;   // nuclear energy (interacted)
 static float e_neutr;  // neutrino energy (not rejestered)
 static float e_stop;   // stopped particle energy (not rejestered)
private:
 int mult;     // count no of hits in this event

};
int   ShowerDetector::evtNum = 0;
float ShowerDetector::e_dedx=0;
float ShowerDetector::e_gamma=0;
float ShowerDetector::e_nucl=0;
float ShowerDetector::e_neutr=0;
float ShowerDetector::e_stop=0;

#include "gismo/PDT.h"
void ShowerDetector::score( MCParticle *p)
{
    static long gammaId = -1;
    static long eplusId, eminusId, pId, nId, deutrId, tritnId, alphaId, nueId, numId;

    if(gammaId == -1) {
       PDT *thePDT = GParticle::thePDT();
       gammaId = (thePDT->lookUp("gamma"))->idCode();
       eplusId = (thePDT->lookUp("e+"))->idCode();
       eminusId= (thePDT->lookUp("e-"))->idCode();
       nueId   = (thePDT->lookUp("nu_e"))->idCode();
       numId   = (thePDT->lookUp("nu_mu"))->idCode();
       pId =     (thePDT->lookUp("p"))->idCode();
       nId =     (thePDT->lookUp("n"))->idCode();
       deutrId = (thePDT->lookUp("deutron"))->idCode();
       tritnId = (thePDT->lookUp("triton"))->idCode();
       alphaId = (thePDT->lookUp("alpha"))->idCode();
   }
   long trackId = p->idCode();
   float eLoss = 0.;
   int numProd = p->numChildren();

// Get the energy loss due neutrino's which leave
  if(trackId == nueId || trackId == numId) {
      e_neutr += p->e();
      return;
   }

// Get the energy loss due to ionization (dE/dX).
   if( p->charge())  e_dedx += p->energyLoss();

// Get the energy loss due to stopped particles (mostly masses)
   if(p->status() == MCParticle::STOPPED) {
      if(trackId == nId )  e_stop += p->e() - p->mass();
      else if(trackId != pId      && trackId != deutrId &&
              trackId != tritnId  && trackId != alphaId &&
              trackId != eminusId && trackId != eplusId)
              e_stop += p->e();
   }
// Get the energy loss due to gamma's which  stop via photo-electric.
   if(trackId == gammaId) {
        if(p->status() == MCParticle::INTERACTED && numProd == 1) {
           GParticle *kid = p->child(0);
           e_gamma += p->e() - (kid->e() - kid->mass());
      }
   }


// If its a interacting hadron correct for nuclear fragments. The correction is
// done by adding the kinetic energy of the incomming paritcle to the energy loss
// and then subtracting the energy (or kinetic energy in the case of p or n) of all
// daughters recognized by gismo.
// NOTE: remove this code whennuclear fragmenets are tracked by gismo
// Deutrons, tritons, and alphas now tracked... wba 15-june-95

    if(trackId != eplusId && trackId != eminusId && trackId != gammaId
       && p->status() == MCParticle::INTERACTED) { // Must be hadronic

         if(trackId == pId     || trackId == nId     ||
            trackId == deutrId || trackId == tritnId ||
            trackId == alphaId )                 e_nucl += p->e() - p->mass();
         else                                    e_nucl += p->e();

      for(int i=0; i<numProd; i++) {
         GParticle *prod = p->child(i);
         if(!prod) break;

         long prodId = prod->idCode();
         if(prodId == pId     || prodId == nId     ||
            prodId == deutrId || prodId == tritnId ||
            prodId == alphaId )   e_nucl -= prod->e() - prod->mass();
         else                     e_nucl -= prod->e();
       }
     }
   return;
  }
//-------------------------------------------------------


class ShowerTest :  public BaseApp
{
  public:
   void endEvent();  // override virtual base class
   void init(ostream& =cout);
   void last(ostream& = cout);

};

void ShowerTest::init(ostream& cout)
{
        Material::addPath("../PEGS4:.");


        // start the particle in the middle of the layer
        setBeam(particleName, Vector(0,0,pBeam), Point(0,0,thick/2));

        // define the geometry
        SlabWorld* w = new SlabWorld;
        setWorld(w);
        w->addSlab(thick,materialName,new ShowerDetector());


        // create particle and medium just to print out parameters
        MCParticle p(particleName);
        p.setMomentum(Vector(0,0,pBeam));
        Medium med(0,new Box(1,1,1),materialName);
        Material& mat = med.material();


        cout << "\tGismo Showering test\n";
        cout << "particle: " << particleName
             << "\n     momentum\t" << p.pTot() << " GeV "
             << "\n            E\t" << p.e() << " GeV"
             << "\n           KE\t" << p.e()-p.mass() << " GeV"
             ;
        Material::printAll(cout);
        cout << "\nmaterial: " << materialName
             << "\n   Rad length\t" << mat.radiationLength() << '\n'
             ;
        Interactor& itr = *p.interactor();
        float ecutMin = itr.ecut(&p,&med) - p.mass();
        if(med.kECutOff() > ecutMin) ecutMin = med.kECutOff();
        cout << "Interactions performed by the class \"" << itr.className() <<'"'
             << "\n   Int length\t" << itr.interactionLength(&p,&med) << " cm"
             << "\n      maxStep\t" << itr.maxStepSize(&p,&med) << " cm"
             << "\n   k.e.cutOff\t" << ecutMin << " GeV"
             ;

        cout << "\n\nbegin run for "<< numberOfEvents<< " events "<< '\n';
        cout << "\n\n\t\t Energy Audit Table \n";
        cout << "Evt.No.  \tdE/dX   \tGamma   \tNeutrino   \tNuclear   \tStopped  \tTotal \n";
}
void ShowerTest::endEvent()
{
    ShowerDetector::evtNum++;

    float e_tot = ShowerDetector::e_dedx  + ShowerDetector::e_gamma +
                  ShowerDetector::e_neutr + ShowerDetector::e_nucl  +
                  ShowerDetector::e_stop;

    energyStat += e_tot;
    int old_flags = cout.flags(ios::fixed);
    int old_precision = cout.precision(3);
    cout << setw(5) << ShowerDetector::evtNum
         << setw(16) << ShowerDetector::e_dedx
         << setw(16) << ShowerDetector::e_gamma
         << setw(16) << ShowerDetector::e_neutr
         << setw(16) << ShowerDetector::e_nucl
         << setw(16) << ShowerDetector::e_stop
         << setw(16) << e_tot << endl;
    cout.flush();
    cout.flags(old_flags);
    cout.precision(old_precision);

}
void ShowerTest::last(ostream& cout)
{
   if( energyStat.samples()>0)
   {
        cout << "\nenergy: \tmean = "<< energyStat.mean() << '\n'
             <<       "\t\trms  = "<< energyStat.stdDev()
             <<'\n';
        cout.flush();
   }
}


int main(int argc, char* argv[])
{

        if( argc>1 ) particleName = argv[1];
        if( strcmp(particleName, "help")==0 ) {cout << help; return 0;}
        if( argc>2 ) materialName = argv[2];
        if( argc>3 ) numberOfEvents = atoi(argv[3]);
        if( argc>4 ) pBeam = atof(argv[4]);
        if( argc>5 ) cut = atof(argv[5]);

  ShowerTest tester;
  tester.addInteractor(new EGSInteractor("electromagnetic") );
#ifndef NOGHEISHA
  tester.addInteractor(new GheishaInteractor("hadronic") );
#endif
  tester.init();
  tester.world().setKECutOff(cut);

  tester.run(numberOfEvents);
  tester.last();


  return 0;
}

---