Energy loss calculations: Bethe-Bloch formula

bb.cpp

// Energy loss calculations
// Bethe-Bloch formula according to C.Grupen & B.Shwartz: Particle Detectors, Cambridge
// Author: Dennis Angemi
// Sources: Material provided by prof.ssa Paola La Rocca, DFA, UNICT, Catania

// constants definition
#define K         0.30707   // 4π N_A r^2_e m_ec^2
#define MEC2      0.511     // m_e c^2
#define N_STEPS   10000     // (10k)
#define Z         2.0       // Z projectile (alpha)
#define M         4.0       // M projectile in a.m.u. (alpha)
#define ZT        13        // Z target (13Al)
#define MT        26.98     // M target in a.m.u. (13Al)
#define E_INITIAL 5         // MeV
#define DENSITY   2.7       // g/cm^3
#define THICKNESS 0.0016    // cm
#define MC2       M*931.494 // mc^2
#define DX        THICKNESS/N_STEPS

typedef TGraph*  pTGraph;
typedef TCanvas* pTCanvas;

void bb()
    { 
        Int_t ie; // iterator
        Double_t ipot, e_tot, p, beta, gamma, arg, phi, e_final, de, stopping_power;
        pTGraph gr;
        pTCanvas c1;
        
        gr = new TGraph;
        
        // ionization potential
        ipot = 16.0*TMath::Power(ZT,0.9);
        if(ZT==1) ipot = 21.8;
        ipot = ipot*1.0E-6;
        
        e_final = E_INITIAL;
        
        // header
        cout << "energy,distance" << endl;
        
        // loop on energy
        for(ie=0; ie<N_STEPS; ie++)
            {   
                de = 0;
                e_tot = e_final + MC2;
                p = TMath::Sqrt( (e_final*e_final) + (2.0*e_final*MC2) );
                beta = p/e_tot;
                gamma = 1.0 / TMath::Sqrt( 1.0 - (beta*beta) );  
                
                arg = 2.0*MEC2*gamma*gamma*beta*beta/ipot;
                phi = log(arg) - (beta*beta);
                stopping_power = K*(ZT/MT)*Z*Z*phi*DENSITY/(beta*beta);
                
                e_final -= stopping_power * DX;
                de += stopping_power * DX;
                if(stopping_power<0.1) break;
                
                cout<<e_final<<","<<stopping_power <<endl; 
                
                // add point to char
                gr -> AddPoint(DX*ie, stopping_power);
                // gr -> AddPoint(beta*gamma, -stopping_power[ie]);
            }
            
        cout << "Spessore percorso: " << DX*ie << endl;
        cout << "Energia: " << e_final << endl;
        
        // plotting
        c1 = new TCanvas();
        gr -> SetMarkerColor(2);
        gr -> Draw("AP");
        gr -> SetMarkerStyle(20);
        gr -> SetMarkerSize(0.2);
        gr->GetXaxis()->SetTitle("Distanza percorsa nel bersaglio (cm)");    
	    gr->GetYaxis()->SetTitle("Perdita di energia (MeV/cm)");
         
    }