---

mscat.cxx

Go to the documentation of this file.

/*
        test multiple scattering and  SlabWorld
 $Id: mscat.cxx,v 1.2 2001/10/06 14:18:53 burnett Exp $

*/

#include "BaseApp.h"
#include "SlabWorld.h"
#include "SmplStat.h"
#include "gismo/Units.h"
#include "geometry/CoordTransform.h"
#include "gismo/EGSInteractor.h"
#include "gismo/GheishaInteractor.h"


#include <iomanip>
SampleStatistic thetaSquared;
SampleStatistic energyLoss;
SampleStatistic multiplicity;
SampleStatistic energyDeposited;

unsigned long numberOfEvents=1000;  // default number of events
unsigned long number=0;
float   pBeam = 2.0;  // initial momentem
float   thickness = 1.0;  // thickness of slab
const char* materialName = "Al";
const char* particleName = "mu+";
float   incidentAngle=0.;

Vector pinit(0,0,pBeam);

float   sigmaMS;            // expected projected MS sigma
float   expectedEloss; //
float   initialEnergy;
//ofstream fout("/annih.out");  //temp

using namespace std;

template<class T> inline T sqr(T x){return x*x;}

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 mscat:\n\n"
"   mscat [[[[particleName] materialName] eventCount] pBeam ] materialsize]\n\n"
"where:\n"
"   particleName [mu+] 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"
"   thickness    [1.0 cm] the amount of material\n";



//-------------------------------------------------------
// Detector class that measures the particles emerging from the target
class Spectrometer : public Detector
{
public:
    Spectrometer():Detector(){}
    void score(MCParticle*);
    void clear(){mult=0; m_edep=0;}
    void generateResponse(){};
private:
    int mult;     // count no of hits in this event
    float m_edep;   // cumulative deposited energy
};

void Spectrometer::score( MCParticle *p)
{
  if( p->charge()==0 || p->pTot()==0 ) return;
  multiplicity += ++mult;
  float deltaE = initialEnergy - p->e();
  m_edep += deltaE;

  if( mult > 1 ) return;       // only use first particle

  double theta = pinit.dot(*p)/p->momentum().mag()/pBeam;
  theta = theta>=1? 0 : acos(theta);
  thetaSquared += sqr(theta);

  energyLoss += deltaE;

}

//-------------------------------------------------------
// Detector class to instrument the target
class Calorimeter: public Detector // this to accumulate energy in the slab
{
public:
    Calorimeter():Detector(){}
    void score(MCParticle* p){
        m_e += p->energyLoss(); 
        ++m_n;
    }
    void clear(){ m_n=0; m_e=0;}
    void generateResponse(){ 
        energyDeposited += m_e;
    }
private:
    float m_e;
    int m_n;
};


//-------------------------------------------------------


class TestApp :  public BaseApp
{
  public:
   void endEvent()
   { if( eventNumber%10==0 )
        //              fout << eventNumber << ' ';
                number++;  // so global is accessible
     //fout << *generator;
   }
   void init(ostream& =cout);
   void last(ostream& = cout);
};

void TestApp::init(ostream& cout)
{

        // define the "beam"
        pinit.rotateY(incidentAngle);
        setBeam(particleName, pinit, Point(0,0,-1));

        // define the geometry
        SlabWorld* w = new SlabWorld;
        setWorld(w);
        w->addSlab(thickness,materialName,new Calorimeter );
        w->addSlab(0.1f,"vacuum",new Spectrometer() );


        // calculate expected multiple scattering sigma and eloss
        GParticle p(particleName);
        p.setMomentum(pinit);
        initialEnergy = p.e();
        double ptotal = p.pTot();
        double beta   = ptotal/initialEnergy;
        double eta      = ptotal/p.mass();

        Material& mat = *Material::hatch(materialName);
        double X0 = mat.radiationLength();
        double step = thickness/p.cosTheta();
        double x = step/X0;

        expectedEloss = step * mat.dedx(beta,eta) * sqr(p.charge());

        sigmaMS = Units::fromMeV(13.6f)
                * sqr(p.charge()) * sqrt(x) * (1.+ 0.038*log(x))
                  /(beta* ptotal);

        cout << "\n\tSimple gismo scattering test";
        cout << "\n\t============================\n";
        static int wd=25;
        cout <<        std::setw(wd) << "particle: " << particleName 
                << '\n' << std::setw(wd) << "momentum: " << ptotal << " GeV"
                << '\n' << std::setw(wd) << "KE: " << (p.e()-p.mass()) << " GeV"
                << '\n' << std::setw(wd) << "incident angle: " << incidentAngle 
                << '\n' << std::setw(wd) << "material: " << materialName
                << '\n' << std::setw(wd) << "thickness: " <<thickness << " cm"
                << '\n' << std::setw(wd) << "Radiation length: " << X0 
                << '\n' << std::setw(wd) << "expected m.s. sigma: " << sigmaMS 
                << '\n' << std::setw(wd) << "expected E loss: " << expectedEloss*1e3 <<" MeV"
//              << '\n' << std::setw(wd) << "Interactions performed by: "<< p.interactor()->className()

                << ", which is "<< expectedEloss*1e3/step/mat.density() <<" MeV/g\n";
        cout << "\nbegin run for "<< numberOfEvents<<" events\n";
}
void TestApp::last(ostream& cout)
{
   const int width=45;
   if( thetaSquared.samples()>0)
   {
        cout << std::setw(width) << "Mean square angular deviation / expected: "
        << thetaSquared.mean() / (2*sqr(sigmaMS))
        << " +/- " << thetaSquared.stdDev()
            /sqrt((double)thetaSquared.samples())
                         /(2*sqr(sigmaMS))
        << '\n';
        cout <<  std::setw(width) <<"Average  energy loss / expected: "
        << energyLoss.mean() / expectedEloss
        << " +/- " << energyLoss.stdDev()
            / sqrt((double)energyLoss.samples())/expectedEloss
        << '\n';
        cout <<  std::setw(width) << "energy loss fraction fluctuation: " 
             <<  energyLoss.stdDev()/energyLoss.mean() <<'\n';
    }
        cout << std::setw(width) <<"Average Multiplicity emerging: " << multiplicity.mean() << '\n';
        cout << std::setw(width) <<"Mean Energy Deposited: " << energyDeposited.mean()*1e3 
             << " MeV, or "<< energyDeposited.mean()*1e6/thickness*1e-4 << " keV/um"<<'\n';
        cout << std::setw(width) << "Energy dep. fractional rms: "
             <<  energyDeposited.stdDev()/energyDeposited.mean() <<'\n';
}

//-------------------------------------------------------
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]); pinit=Vector(0,0,pBeam);}
  if( argc>5 ) {thickness = atof(argv[5]);}
  
  TestApp tester;

  tester.addInteractor(new EGSInteractor("electromagnetic") );
  tester.addInteractor(new GheishaInteractor("hadronic"));

  tester.init();
  tester.run(numberOfEvents);
  tester.last();
  return 0;
}

---