Index: trunk/CVSROOT/modules =================================================================== --- trunk/CVSROOT/modules (revision 1348) +++ trunk/CVSROOT/modules (revision 1349) @@ -96,2 +96,7 @@ Cosy MagicSoft/Cosy & tcpip SubsysIO MagicSoft/Control/SubsystemIO +# +# Wuerzburg private modules +# +mphys -dmphys WuerzburgSoft/Thomas/mphys +Thomas MagicSoft/Mars & mmc & mphys Index: trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1349) @@ -0,0 +1,279 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MElectron.h" + +#include + +#include + +ClassImp(MElectron); + +Double_t MElectron::Planck(Double_t *x, Double_t *k=NULL) +{ + // + // Planck, per unit volume, per unit energy + // + // constants moved out of function + // + Double_t E = x[0]; // [GeV] + Double_t z = k ? k[0] : 0; + + Double_t T = 2.96*(z+1); // [K] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + + Double_t EkT = E/kB/T; + + /* + //Double_t c = 299792458; // [m/s] + //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + //Double_t hc = h*c; // [GeVm] + Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ] + */ + + return E*E / (exp(EkT)-1.); // [GeV^2] +} + +Double_t MElectron::Li(Double_t *x, Double_t *k) +{ + Double_t t = x[0]; + + return log(1.-t)/t; +} + +Double_t MElectron::Li2(Double_t *x, Double_t *k=NULL) +{ + // + // Dilog, Li2 + // + Double_t z = x[0]; + + TF1 IntLi("Li", Li, 0, z, 0); + Double_t integ = IntLi.Integral(0, z); + + /* + if (fabs(z)<1) + { + Double_t disum = DiSum(&z); + cout << "Disum (" << z << ") " << disum << "=" << -integ << "\t" << disum-integ << endl; + return disum; + } + */ + + /* + Integral(0, 1) = konst; + Double_t konst = 1./6*TMath::Pi()*TMath::Pi(); + */ + + return -integ; +} + +Double_t MElectron::Flim(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam) mit b-->1 (Maple) +{ + Double_t w = x[0]; + + Double_t w4 = w*4; + Double_t wsqr = w*w; + + Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1); + Double_t u2 = -w4*(wsqr*2 + w*9 + 2); + Double_t d = w4*(w4 + 1); + + Double_t s = -w*2*(1+1); // -2*omega*(1+beta) + Double_t li2 = Li2(&s); + + Double_t res = (u1+u2)/d + li2; + + return res; //<1e-10? 0 : res; +} + +Double_t MElectron::Compton(Double_t *x, Double_t *k) +{ + Double_t E0 = 511e-6; //[GeV] + Double_t E02 = E0*E0; + + Double_t epsilon = x[0]; + Double_t E = k[0]; + // Double_t beta = k[1]; + Double_t z = k[2]; + + Double_t omega = epsilon*E/E02; + + Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1] + return Flim(&omega)*n; +} + + +Double_t MElectron::InteractionLength(Double_t *E, Double_t *k=NULL) +{ + // E = electron energy, ~ TeV(?) 1e12 + // e = photon energy, ~ meV(?) 1e-3 + // mc^2 = electron rest mass energy ~.5keV(?) .5e3 + // + // x^-1 = int( n(epsilon)/2beta * ((mc^2)^2/eE)^2 * int ( omega*sigma(omega), omega=o-..o+), epsilon=0..inf) + // + // o+/- = omage_0 (1 +- beta) + // + // omega_0 = eE/(mc^2)^2 ~1e12*1e-3/.25e6=4e3 + // + // --> x^-1 = (alpha*hc)^2/4pibetaE^2 * int(n(epsilon)/epsilon^2 *( F(o+)-F(o-)), epsilon=0..inf) + // + // F(o) = -o/4 + (9/4 + 1/o + o/2) * ln(1+2o) + 1/8(1+2o) - 3/8 + Li2(-2o) + // + // Li2(x) = int(ln(1-t)/t, t=0..x) + // + // F(o+)~F(2o) = -o/2 + (9/4 + 1/2o + o) * ln(1+4o) + 1/8(1+4o) - 3/8 + Li2(-4o) + // F(o-)~F(0) = 14/8 = 1.75 + + Double_t E0 = 511e-6; // [GeV] + Double_t E02 = E0*E0; // [GeV^2] + Double_t c = 299792458; // [m/s] + Double_t e = 1.602176462e-19; // [C] + Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + Double_t hc = h*c; // [GeVm] + Double_t alpha = 1./137.; // [1] + + Double_t z = k ? k[0] : 0; + + // Double_t beta = sqrt(1-E0/E*E0/E); + Double_t beta = 1; //0.999999999999999999999999999; + + Double_t val[3] = { E[0], beta, z }; // E[GeV] + + Double_t from = 1e-17; + Double_t to = 1e-11; + + /* -------------- old ---------------- + Double_t from = 1e-15; + Double_t to = 1e-11; + eps = [default]; + ----------------------------------- + */ + TF1 func("Compton", Compton, from, to, 3); // [0, inf] + + Double_t integ = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf] + + Double_t aE = alpha/E[0]; // [1/GeV] + + Double_t konst = 2.*E02/hc/beta; // [1 / GeV m] + Double_t ret = konst * (aE*aE) * integ; // [1 / m] + + Double_t ly = 3600.*24.*365.*c; // [m/ly] + Double_t pc = 1./3.258; // [pc/ly] + + return (1./ret)/ly*pc/1000; // [kpc] +} + +Double_t MElectron::GetInteractionLength(Double_t energy, Double_t z) +{ + return InteractionLength(&energy, &z); +} + +Double_t MElectron::GetInteractionLength() const +{ + return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ); +} + +// -------------------------------------------------------------------------- + +Double_t MElectron::p_e(Double_t *x, Double_t *k) +{ + Double_t e = pow(10, x[0]); + Double_t E = k[0]; + Double_t z = k[1]; + + Double_t E0 = 511e-6; //[GeV] + Double_t E02 = E0*E0; + + Double_t omega = e*E/E02; + + Double_t n = Planck(&e, &z); + + Double_t F = Flim(&omega)/omega/omega; + + return n*F*1e26; +} + +Double_t MElectron::G_q(Double_t *x, Double_t *k) +{ + Double_t q = x[0]; + Double_t Gamma = k[0]; + + Double_t Gq = Gamma*q; + + Double_t s1 = 2.*q*log(q); + Double_t s2 = (1.+2.*q); + Double_t s3 = (Gq*Gq)/(1.+Gq)/2.; + + return s1+(s2+s3)*(1.-q); +} + + +Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k = NULL) +{ + Double_t E = x[0]; + Double_t z = k ? k[0] : 0; + + Double_t E0 = 511e-6; //[GeV] + + Double_t lolim = -log10(E)/7*4-13; + + TF1 fP("p_e", p_e, lolim, -10.8, 2); + TF1 fQ("G", G_q, 0, 1., 1); + + fP.SetParameter(0, E); + fP.SetParameter(1, z); + + Double_t e = pow(10, fP.GetRandom()); + + Double_t omega = e*E/E0/E0; + Double_t Gamma = 4.*omega; + + // --old-- fQ.SetRange(1e-6, 1./(1.+ 2.*Gamma)); + fQ.SetParameter(0, Gamma); + + Double_t q = fQ.GetRandom(); + Double_t Gq = Gamma*q; + + Double_t e1 = Gq*E/(1.+Gq); + + return e1; +} + +Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z) +{ + return EnergyLoss(&E, &z); +} + +Double_t MElectron::GetEnergyLoss() const +{ + return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ); +} Index: trunk/WuerzburgSoft/Thomas/mphys/MElectron.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1349) @@ -0,0 +1,41 @@ +#ifndef MARS_MElectron +#define MARS_MElectron + +#ifndef MARS_MParticle +#include "MParticle.h" +#endif + +class MElectron : public MParticle +{ +public: + MElectron(Double_t e=0, Double_t z=0) : MParticle(MParticle::kEGamma) + { + fEnergy = e; + fZ = z; + } + + // ---------------------------------------------------------------- + + static Double_t Planck(Double_t *x, Double_t *k=NULL); + static Double_t Li(Double_t *x, Double_t *k); + static Double_t Li2(Double_t *x, Double_t *k=NULL); + static Double_t Flim(Double_t *x, Double_t *k=NULL); + static Double_t Compton(Double_t *x, Double_t *k); + static Double_t InteractionLength(Double_t *E, Double_t *k=NULL); + static Double_t GetInteractionLength(Double_t E, Double_t z=0); + + Double_t GetInteractionLength() const; + + // ---------------------------------------------------------------- + + static Double_t p_e(Double_t *x, Double_t *k); + static Double_t G_q(Double_t *x, Double_t *k); + static Double_t EnergyLoss(Double_t *E, Double_t *z=NULL); + static Double_t GetEnergyLoss(Double_t E, Double_t z=0); + + Double_t GetEnergyLoss() const; + + ClassDef(MElectron, 1) +}; + +#endif Index: trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc (revision 1349) @@ -0,0 +1,113 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MGenIRPhoton.h" + +#include +#include + +#include "MParList.h" +#include "MParticle.h" + +ClassImp(MGenIRPhoton); + +/* +Double_t Planck(Double_t *x, Double_t *k) +{ + const Double_t E = x[0]; //[eV] + + const Double_t T = 2.87; //[K] + + const Double_t e = 1.602176462e-19; //[C] + const Double_t kB = 1.3806503e-23/e; //[eV/K] + + const Double_t E3 = E*E*E; + const Double_t EkT = E/kB/T; + + //const Double_t h = 1e-9/e*6.62606876e-34; //[GeVs] + //const Double_t h2 = h*h; + //const Double_t c = 299792458; //[m/s] + //const Double_t c3 = c*c*c; + //const Double_t konst = 4/c3/h2; + + return E3/(exp(EkT)-1); +} +*/ + +Double_t MGenIRPhoton::Planck(Double_t *x, Double_t *k=NULL) +{ + // + // Planck, per unit volume, per unit energy + // + // constants moved out of function + // + Double_t E = x[0]; // [GeV] + Double_t z = k ? k[0] : 0; + + Double_t T = 2.96*(z+1); // [K] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + + Double_t EkT = E/kB/T; + + /* + //Double_t c = 299792458; // [m/s] + //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + //Double_t hc = h*c; // [GeVm] + Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ] + */ + + return E*E / (exp(EkT)-1.); // [GeV^2] +} + + +// -------------------------------------------------------------------------- +MGenIRPhoton::MGenIRPhoton() +{ + fSrc = new TF1("Planck", Planck, 0, .5e-2, 0); +} + +MGenIRPhoton::~MGenIRPhoton() +{ + delete fSrc; +} + +// -------------------------------------------------------------------------- +MParticle *MGenIRPhoton::GetRandom() +{ + const Double_t pi2 = TMath::Pi() * 2; + MParticle *p = new MParticle(MParticle::kEGamma); + + p->SetAngle(fRand.Uniform(pi2)); + p->SetEnergy(fSrc->GetRandom()*1e-9); + + return p; +} + Index: trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.h (revision 1349) @@ -0,0 +1,34 @@ +#ifndef MARS_MGenIRPhoton +#define MARS_MGenIRPhoton + +#ifndef ROOT_TRandom +#include +#endif + +class TF1; +class MParticle; + +class MGenIRPhoton +{ +private: + TF1 *fSrc; + + TRandom fRand; + + Double_t fZ; + + static Double_t Planck(Double_t *x, Double_t *k=NULL); + +public: + MGenIRPhoton(); + virtual ~MGenIRPhoton(); + + void SetZ(Double_t z) { fZ = z; } + + MParticle *GetRandom(); + + ClassDef(MGenIRPhoton, 0) // +}; + +#endif + Index: trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc (revision 1349) @@ -0,0 +1,81 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MGenPrimaryParticle.h" + +#include +#include + +#include "MParList.h" +#include "MParticle.h" + +ClassImp(MGenPrimaryParticle); + +Double_t PowerLaw(Double_t *x, Double_t *k) +{ + return 1/x[0]; +} + +// -------------------------------------------------------------------------- +MGenPrimaryParticle::MGenPrimaryParticle(const char *name, const char *title) + : fPi2(TMath::Pi()*2) +{ + fName = name ? name : "MGenPrimaryParticle"; + fTitle = title ? title : ""; + + fSrc = new TF1("Power-Law", PowerLaw, 10, 1e8, 0); // 10GeV-100PeV +} + +MGenPrimaryParticle::~MGenPrimaryParticle() +{ + delete fSrc; +} + +// -------------------------------------------------------------------------- +Bool_t MGenPrimaryParticle::PreProcess(MParList *pList) +{ + fList = (MParList*)pList->FindCreateObj("MParticleList"); + if (!fList) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +Bool_t MGenPrimaryParticle::Process() +{ + MParticle *p = new MParticle(MParticle::kEGamma); + p->SetEnergy(fSrc->GetRandom()); + return kTRUE; +} + +Bool_t MGenPrimaryParticle::PostProcess() +{ + return kTRUE; +} Index: trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.h (revision 1349) @@ -0,0 +1,39 @@ +#ifndef MARS_MGenPrimaryParticle +#define MARS_MGenPrimaryParticle + +#ifndef ROOT_TRandom +#include +#endif + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class TF1; +class MParList; +class MSinglePairInput; + +class MGenPrimaryParticle : public MTask +{ +private: + MParList *fList; + + TF1 *fSrc; + + TRandom fRand; + + Double_t fPi2; + +public: + MGenPrimaryParticle(const char *name=NULL, const char *title=NULL); + ~MGenPrimaryParticle(); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + ClassDef(MGenPrimaryParticle, 0) // +}; + +#endif + Index: trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1349) @@ -0,0 +1,150 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MPairProduction.h" + +#include +#include + +#include "TList.h" +#include "MElectron.h" + +ClassImp(MPairProduction); + +Double_t AngleDistrib(Double_t *x, Double_t *k) +{ + + const Double_t c = x[0]; // cos(alpha) + const Double_t b = k[0]; // sqrt(1-1/s) + + const Double_t b2 = b*b; + const Double_t b4 = b2*b2; + + const Double_t c2 = c*c; + const Double_t c4 = c2*c2; + + const Double_t u = 1 - b4*c4 +2*b2*(1-b2)*(1-c2); + const Double_t d = 1-b2*c2; + + return u/(d*d); +} + +// -------------------------------------------------------------------------- +MPairProduction::MPairProduction() +{ + fAngle = new TF1("AngleDistrib", AngleDistrib, -1, 1, 1); +} + +MPairProduction::~MPairProduction() +{ + delete fAngle; +} + +#include + +// -------------------------------------------------------------------------- +Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, TList *list) +{ + // + // gamma: primary particle from source + // phot: infrared photon from background. (angle = interaction angle) + // + const Double_t E0 = 511e-6; // [GeV] + + const Double_t theta = phot->GetAngle(); // [2pi] + const Double_t Ep = phot->GetEnergy(); // [GeV] + const Double_t Eg = gamma->GetEnergy(); // [GeV] + + const Double_t ctheta = cos(theta); + + const Double_t s = (Ep*Eg/(2*E0*E0))*(1-ctheta); //[1] + if (s<1) + return kFALSE; + + const Double_t Epg = Ep/Eg; // 1+Epg ~ 1 + const Double_t Egp = Eg/Ep; // 1+Egp ~ Egp + + const Double_t phi = atan(sin(theta)/(Egp+ctheta)); + + const Double_t sphi = sin(phi); + const Double_t cphi = cos(phi); + + const Double_t alpha = theta-phi; + const Double_t calpha = cos(alpha); + const Double_t salpha = sin(alpha); + + const Double_t beta = (Eg*cphi+Ep*calpha)/(Ep+Eg); + + // + // gamma1 = 1/gamma = sqrt(1-beta*beta) + // + const Double_t gamma1 = sqrt((sphi*phi+Epg*salpha*Epg*salpha+2*Epg*(1-cphi*calpha))); + + const Double_t Beta = sqrt(1-1/s); //[1] + + fAngle->SetParameter(0, Beta); + + const Double_t psi = atan(gamma1*sphi/(cphi-beta)); + const Double_t delta = acos(fAngle->GetRandom()) - psi; + + const Double_t Bcosd = Beta*cos(delta); + const Double_t Bsind = Beta*sin(delta); + + const Double_t E = sqrt(s)*E0/gamma1; + const Double_t dE = E*Bcosd; + + MElectron *p[2]; + p[0] = new MElectron(E+dE, gamma->GetZ()); + p[1] = new MElectron(E-dE, gamma->GetZ()); + + const Double_t E1 = E0/(E+dE); + const Double_t E2 = E0/(E-dE); + + const Double_t beta1 = sqrt(1.-E1*E1); + const Double_t beta2 = sqrt(1.-E2*E2); + + p[0]->SetBeta(beta1); + p[1]->SetBeta(beta2); + + const Double_t Bscp = Bsind*cphi; + const Double_t spg = sphi/gamma1; + const Double_t cpg = cphi/gamma1; + + const Double_t tan1 = (spg*(Bcosd+1) + Bscp)/((cpg*(Bcosd+1) - Bscp)); + const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp)); + + p[0]->SetAngle(atan(tan1)); + p[1]->SetAngle(atan(tan2)); + + list->Add(p[0]); + list->Add(p[1]); + + return kTRUE; +} + Index: trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 1349) @@ -0,0 +1,27 @@ +#ifndef MARS_MPairProduction +#define MARS_MPairProduction + +#ifndef ROOT_TObject +#include +#endif + +class TF1; +class MParticle; +class TList; + +class MPairProduction +{ +private: + TF1 *fAngle; + +public: + MPairProduction(); + virtual ~MPairProduction(); + + Bool_t Process(MParticle *g, MParticle *p, TList *l); + + ClassDef(MPairProduction, 0) // +}; + +#endif + Index: trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1349) @@ -0,0 +1,26 @@ +/////////////////////////////////////////////////////////////////////// +// +// MParticle +// +/////////////////////////////////////////////////////////////////////// + +#include "MParticle.h" + +ClassImp(MParticle); + +MParticle::MParticle(ParticleType_t t, const char *name, const char *title) + : fPType(t) +{ + // + // default constructor + // creates an a list of histograms for all pixels and both gain channels + // + + // + // set the name and title of this object + // + + fName = name ? name : "MParticle"; + fTitle = title ? title : "Container for a particle"; +} + Index: trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1349) @@ -0,0 +1,43 @@ +#ifndef MARS_MParticle +#define MARS_MParticle + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MParticle : public MParContainer +{ +public: + typedef enum { kEGamma, kEElectron, kEPositron } ParticleType_t; + +private: + const ParticleType_t fPType; + +protected: + Double_t fEnergy; // [GeV] Energy + Double_t fBeta; // [1] beta + Double_t fAngle; // [rad] Angle + Double_t fZ; // [1] Red shift + +public: + MParticle(ParticleType_t t, const char *name=NULL, const char *title=NULL); + // ~MParticle(); + + // void Fill(const MParContainer *inp); + + // void Draw(Option_t *opt=NULL); + + void SetEnergy(Double_t e) { fEnergy = e; } + void SetAngle(Double_t a) { fAngle = a; } + void SetBeta(Double_t b) { fBeta = b; } + void SetZ(Double_t z) { fZ = z; } + + Double_t GetEnergy() const { return fEnergy; } + Double_t GetAngle() const { return fAngle; } + Double_t GetZ() const { return fZ; } + + ClassDef(MParticle, 1) // Container which holds hostograms for the Hillas parameters +}; + +#endif + Index: trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1349) @@ -0,0 +1,193 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) +! Author(s): Thomas Bretz 12/2000 (tbretz@uni-sw.gwdg.de) +! +! Copyright: MAGIC Software Development, 2000-2001 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// // +////////////////////////////////////////////////////////////////////////////// +#include "MPhoton.h" + +#include + +#include + +ClassImp(MPhoton); + +Double_t MPhoton::Planck(Double_t *x, Double_t *k=NULL) +{ + // + // Planck, per unit volume, per unit energy + // + // constants moved out of function, see below + // + Double_t E = x[0]; // [GeV] + Double_t z = k ? k[0] : 0; + + Double_t T = 2.96*(z+1); // [K] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + + Double_t EkT = E/kB/T; + + /* + //Double_t c = 299792458; // [m/s] + //Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + //Double_t hc = h*c; // [GeVm] + + Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ] + */ + + return E*E / (exp(EkT)-1.); // [GeV^2] +} + +Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k=NULL) +{ + Double_t s = x[0]; // omega: CM mass + + Double_t E0 = 511e-6; // [GeV] + Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2 + + Double_t m = E0/s; + + Double_t m2 = m*m; + Double_t beta = sqrt(1.-m2); + Double_t beta2 = 1.-m2; + + Double_t p1 = r0*r0*TMath::Pi()/2; + + // ----- Extreme Relativistic ----- + // return p1*2 * m*m*m* (log(2./m)-1); + // -------------------------------- + + Double_t p2 = m2; + Double_t p3 = 3.-beta2*beta2; + Double_t p4 = log((1.+beta)/(1.-beta)); + Double_t p5 = beta*2*(1.+m2); + + Double_t sigma = p1*p2*(p3*p4-p5); // [m^2] + + return sigma; +} + +Double_t MPhoton::Int1(Double_t *x, Double_t *k=NULL) +{ + Double_t costheta = x[0]; + + Double_t Eg = k[0]; + Double_t Ep = k[1]; + + Double_t E0 = 511e-6; // [GeV] + + Double_t s = Eg*Ep/E0*(1.-costheta)*2; + + if (s Int2==0 <---" << endl; + return 0; + } + + /* Planck constants: konst */ + Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc); + int2 *= konst; // [1 / m] + + Double_t res = 1./ int2; // [m] + res *= pc/ly * 1e-3; // [kpc] + + if (res > 1e50) return 1e50; + if (res < 0) return 1e35; + + return res; //[kpc] +} + +Double_t MPhoton::GetInteractionLength(Double_t energy, Double_t z) +{ + return InteractionLength(&energy, &z); +} + +Double_t MPhoton::GetInteractionLength() const +{ + return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ); +} + Index: trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1349) @@ -0,0 +1,29 @@ +#ifndef MARS_MPhoton +#define MARS_MPhoton + +#ifndef MARS_MParticle +#include "MParticle.h" +#endif + +class MPhoton : public MParticle +{ +public: + MPhoton(Double_t e=0, Double_t z=0) : MParticle(MParticle::kEGamma) + { + fEnergy = e; + fZ = z; + } + + static Double_t Planck(Double_t *x, Double_t *k=NULL); + static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL); + static Double_t Int1(Double_t *x, Double_t *k=NULL); + static Double_t Int2(Double_t *x, Double_t *k); + static Double_t InteractionLength(Double_t *x, Double_t *k=NULL); + static Double_t GetInteractionLength(Double_t energy, Double_t z=0); + + Double_t GetInteractionLength() const; + + ClassDef(MPhoton, 1) +}; + +#endif Index: trunk/WuerzburgSoft/Thomas/mphys/Makefile =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/Makefile (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/Makefile (revision 1349) @@ -0,0 +1,52 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +# +# Handling name of the Root Dictionary Files +# +CINT = Phys + +# +# Library name to creatre +# +LIB = mphys.a + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase + +#------------------------------------------------------------------------------ + +.SUFFIXES: .c .cc .cxx .h .hxx .o + +SRCFILES = MParticle.cc \ + MPhoton.cc \ + MElectron.cc \ + MGenIRPhoton.cc \ + MPairProduction.cc \ + MGenPrimaryParticle.cc + +SRCS = $(SRCFILES) +HEADERS = $(SRCFILES:.cc=.h) +OBJS = $(SRCFILES:.cc=.o) + +############################################################ + +all: $(LIB) + +include ../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib + +mrproper: clean rmbak + +# @endcode + Index: trunk/WuerzburgSoft/Thomas/mphys/PhysIncl.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/PhysIncl.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/PhysIncl.h (revision 1349) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h (revision 1349) @@ -0,0 +1,14 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MParticle+; +#pragma link C++ class MElectron+; +#pragma link C++ class MPhoton+; +#pragma link C++ class MGenPrimaryParticle+; +#pragma link C++ class MGenIRPhoton+; +#pragma link C++ class MPairProduction+; + +#endif Index: trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1349) +++ trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1349) @@ -0,0 +1,1066 @@ +//#ifndef __CINT__ +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +//#endif +#include "mbase/MParContainer.h" +#include "mphys/MPhoton.h" +#include "mphys/MElectron.h" +#include "mphys/MPairProduction.h" +#include "mhist/MBinning.h" +#include "mhist/MH.h" +// 2.96 +// 2.87 +// 2.73 + +/*********************** + * + * calc r from z: + * + * R = 2*c/H0 *(1+z - sqrt(1+z)) + * + * z = -0.5 * (3 + R' +/- sqrt(R'+1)) R' = R*H0/c + * + * H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + * + ************************ + */ + +Double_t ZofR(Double_t *x, Double_t *k=NULL) +{ + Double_t c = 299792458; // [m/s] + Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + Double_t ly = 3600.*24.*365.*c; // [m/ly] + Double_t pc = 1./3.258; // [pc/ly] + Double_t r = x[0] /pc*ly*1e3; // [m] + + Double_t R = r*H0/c; // [1] + + return (R+1+sqrt(R*2+1))/2 - 1; +} + +Double_t RofZ(Double_t *x, Double_t *k=NULL) +{ + Double_t z1 = x[0] + 1; + + Double_t c = 299792458; // [m/s] + Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + Double_t ly = 3600.*24.*365.*c; // [m/ly] + Double_t pc = 1./3.258; // [pc/ly] + + Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m] + + return R * pc/ly/1e3; // [kpc] +} + +Double_t AngleDistrib(Double_t *x, Double_t *k) +{ + + Double_t c = x[0]; // cos(alpha) + Double_t b = k[0]; // sqrt(1-1/s) + + Double_t b2 = b*b; + Double_t b4 = b2*b2; + + Double_t c2 = c*c; + Double_t c4 = c2*c2; + + Double_t u = 1. - b4*c4 +2.*b2*(1.-b2)*(1.-c2); + Double_t d = 1.-b2*c2; + + return u/(d*d); +} + + +// ================================================================ +Double_t DiSum(Double_t *x, Double_t *k=NULL) +{ + Double_t t = x[0]; + + Double_t disum = t; + Double_t add = 0; + + Double_t eps = fabs(t*1e-1); + + Int_t n = 2; + Double_t pow = t*t; // t^2 + + do + { + add = pow/n/n; + + pow *= t; // pow = t^n + n++; + + disum += add; + + } while (fabs(add)>eps); + + return disum; +} + +Double_t F(Double_t *x, Double_t *k=NULL) +{ + Double_t o = x[0]; + Double_t s = -2.*o; + +// if (o<1e-10) +// return 2.125; //-3./8.; + + return -o/4. + (9./4. + 1./o + o/2.) * log(1.+2.*o) + 1./8.*(1.+2.*o) + MElectron::Li2(&s); //- 3./8. +} + +Double_t FlimLi(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam) mit b-->1 (Maple) +{ + Double_t w = x[0]; + + Double_t s = -w*2*(1+1); // -2*omega*(1+beta) + Double_t li2 = MElectron::Li2(&s); + + return fabs(li2); +} + +#include +void rkck(Double_t y[], Double_t dydx[], int n, Double_t x, Double_t h, Double_t yout[], + Double_t yerr[], void (*derivs)(Double_t, Double_t[], Double_t[])) +{ + /* + * --------------------------------------------------------- + * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method + * --------------------------------------------------------- + */ + const Double_t a2 = 0.2; + const Double_t a3 = 0.3; + const Double_t a4 = 0.6; + const Double_t a5 = 1.0; + const Double_t a6 = 0.875; + const Double_t b21 = 0.2; + const Double_t b31 = 3.0/40.0; + const Double_t b32 = 9.0/40.0; + const Double_t b41 = 0.3; + const Double_t b42 = -0.9; + const Double_t b43 = 1.2; + const Double_t b51 = -11.0/54.0; + const Double_t b52 = 2.5; + const Double_t b53 = -70.0/27.0; + const Double_t b54 = 35.0/27.0; + const Double_t b61 = 1631.0/55296.0; + const Double_t b62 = 175.0/512.0; + const Double_t b63 = 575.0/13824.0; + const Double_t b64 = 44275.0/110592.0; + const Double_t b65 = 253.0/4096.0; + const Double_t c1 = 37.0/378.0; + const Double_t c3 = 250.0/621.0; + const Double_t c4 = 125.0/594.0; + const Double_t c6 = 512.0/1771.0; + const Double_t dc5 = -277.00/14336.0; + + const Double_t dc1 = c1-2825.0/27648.0; + const Double_t dc3 = c3-18575.0/48384.0; + const Double_t dc4 = c4-13525.0/55296.0; + const Double_t dc6 = c6-0.25; + + Double_t ak2[n]; + Double_t ak3[n]; + Double_t ak4[n]; + Double_t ak5[n]; + Double_t ak6[n]; + Double_t ytemp[n]; + + for (int i=0; ib ? a : b; +} + +void SolvEq(Double_t y[], Double_t dydx[], int n, Double_t *x, Double_t htry, Double_t eps, + Double_t yscal[], Double_t *hdid, Double_t *hnext, + void (*derivs)(Double_t, Double_t[], Double_t[])) // rkqs +{ + /* + * --------------------------------------------------------- + * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method + * --------------------------------------------------------- + */ + const Double_t SAFETY = 0.9; + const Double_t PGROW = -0.2; + const Double_t PSHRNK = -0.25; + const Double_t ERRCON = 1.89e-4; + + Double_t yerr[n]; + Double_t ytemp[n]; + + Double_t h = htry; + Double_t errmax; + while (1) + { + rkck(y, dydx, n, *x, h, ytemp, yerr, derivs); + + errmax=0.0; + + for (int i=0; i= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h)); + + Double_t xnew= (*x) + h; + + if (xnew != *x) + continue; + + cout << "stepsize underflow in rkqs" << endl; + break; + } + + *hnext = errmax>ERRCON ? SAFETY*h*pow(errmax,PGROW) : 5.0*h; + + *x += (*hdid=h); + + for (int i=0; i=0 ? ZofR(&R) : 0; + + Double_t e = 1.602176462e-19; // [C] + Double_t T = 2.96*(z+1); // [K] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + Double_t kBT = kB*T; // [GeV] + Double_t alpha = 1./137; // [|] + Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + Double_t E0 = 511e-6; // [GeV] + + + Double_t k = TMath::Pi() * alpha * kBT; + + Double_t ln = 4.*kBT/E0/E0; + + return -1./3/h* k*k * (log(ln*E)-1.9805); +} + +void dEdt(Double_t t, Double_t E[], Double_t dedt[]) +{ + dedt[0] = dEdt(E[0], t); + //cout << t << "\t" << E[0] << "\t" << dedt[0] << endl; +} + +Int_t GetSeed() +{ + return (Int_t)fmod(TStopwatch::GetRealTime()*10, 65535); +} + +void DrawDevelopmentHiLim(Double_t E0, Double_t z, Option_t *opt="") +{ + Double_t t = 0; + Double_t E[1] = { E0 }; + Double_t yscal[1] = { 1 }; + + Double_t dedt[1]; + + Double_t eps;// = 1e5; + Double_t step = 5e6; + + Double_t hdid; + Double_t hnext; + + cout << "Start: " << dedt[0] << endl; + + Int_t n = 15000; + Double_t tres[n]; + Double_t Eres[n]; + int i; + for (i=0; i1.5e14 || E[0]<5e6) + break; +// cout << tres[i] << "\t" << Eres[i] << "\t(" << step << ")" << endl; + } + + cout << i << endl; + + TGraph grp(iDraw(opt); + +} + +Double_t EnergyLossRateLoLim(Double_t *x, Double_t *k=NULL) +{ + Double_t t = x[0]; + Double_t E = k[0]; + Double_t t0 = k[1]; + + Double_t c = 299792458; // [m/s] + Double_t ly = 3600.*24.*365.*c; // [m/ly] + Double_t pc = 1./3.258; // [pc/ly] + + Double_t r = t * c / ly * pc / 1000; // [kpc] + Double_t R = RofZ(&k[2]) - r; + Double_t z = k[2]>=0 ? ZofR(&R) : 0; + + Double_t T = 2.96*(z+1); // [K] + // Double_t alpha = 1./137; // [1] + Double_t sigma = 6.653e-29; // [m^2] + Double_t E0 = 511e-6; // [GeV] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + Double_t h = 1e-9/e*6.62606876e-34; // [GeVs] + Double_t hc = h*c; // [GeVm] + Double_t pi = TMath::Pi(); // [1] + Double_t khc = pi*kB/hc; // [1 / K m] + Double_t a = 8./15 * pi * khc*khc*khc*khc * hc; // [Gev / K^4 / m^3] + Double_t konst = 4./3 * sigma * a * T*T*T*T * c / (E0* E0); // [1 / GeV s] + + Double_t ret = 1./(konst*(t-t0) + 1./E); + + return ret; +} + +void DrawDevelopmentLoLim(Double_t t0, Double_t E0, Double_t z=-1, Option_t *opt="") +{ + // 8.7 + + Double_t val[] = { E0, t0, z }; + Double_t t = 1.5e14; + while (EnergyLossRateLoLim(&t, val)<1e4) + t -= 0.01e14; + + TF1 *f1=new TF1("LoLim", EnergyLossRateLoLim, t, 1.5e14, 2); + f1->SetParameter(0, E0); + f1->SetParameter(1, t0); + + f1->Draw(opt); + f1->SetBit(kCanDelete); +} + +// +// (3) Energy loss rate of electrons and 'high energy particle' +// +Double_t DrawDevelopment(Double_t E, Double_t z, Option_t *opt="", TH1 *hist=NULL) +{ + Double_t ly = 3600.*24.*365.; // [s/ly] + Double_t pc = 1./3.258; // [pc/ly] + + TGraph *graph = new TGraph; + graph->SetPoint(0, 0, E); + graph->SetMaximum(E*3); // *MENU* + + cout << "------ " << endl; + + static TRandom rand; + + Double_t x = 0; + for (int i=1; i<10; i++) + { + Double_t l = rand.Exp(MElectron::InteractionLength(&E, &z)); + + if (z>=0) + { + Double_t r = RofZ(&z) - l; + + cout << " " << i << ". R=" << RofZ(&z) << " l=" << l << " z=" << ZofR(&r) << endl; + + z = r>=0 ? ZofR(&r) : 0; + + if (z==0) + cout << "z<0" << endl; + } + + x += l; + + Double_t t = x/pc*ly*1000; + graph->SetPoint(i*2-1, t, E); + + Double_t e1 = MElectron::GetEnergyLoss(E, z<0?0:z); + + E -= e1; + + if (hist) + hist->Fill(e1); + + cout << " Ep=" << e1 << flush; + + graph->SetPoint(i*2, t, E); + } + + graph->SetMinimum(E/3); // *MENU* + graph->Draw(opt); + graph->SetBit(kCanDelete); + + //if (E<31500) + cout << "t=" << x*ly/pc*1000 << "\tE=" << E << "\tz=" << z << endl; + + return E; +} + +void EnergyLossRate() +{ + if (gPad) + gPad->Clear(); + + Double_t E = 1.5e9; // [GeV] + Double_t z = 0.03; + + MBinning bins; + bins.SetEdgesLog(18, 0.1, 1e9); + + TH1D hist; + hist.SetName("Phot"); + hist.SetTitle("Photons from inverse Compton"); + MH::SetBinning(&hist, &bins); + + cout << "Working..." << flush; + + for (int i=0; i<50; i++) + { + cout << i << "." << flush; + DrawDevelopment(E, z, i?"L":"AL", &hist); + } + + //DrawDevelopmentLoLim(2e14, 1.64e2, "Lsame"); // seen + DrawDevelopmentLoLim(1.78e14, 280, z, "Lsame"); // seen + DrawDevelopmentHiLim(E, z, "L"); + gPad->SetLogy(); + + new TCanvas("Photons", "Photons created in inverse Compton"); + hist.DrawCopy(); + gPad->SetLogx(); + + cout << "...done." << endl; +} + +void Scatter() +{ + Double_t E = 1.5e9; // [GeV] + + Int_t nbins = 100; + Double_t lo = E/1e6; + Double_t up = E*10; + + Double_t binsize = log10(up/lo)/nbins; + Double_t *scale=new Double_t[nbins+1]; + for (int i=0; iDraw(); + + TStopwatch clock; + clock.Start(); + + for (int i=0; i<1000; i++) + { + Double_t e1 = epsilon.GetRandom(); + + hist->Fill(e1); + } + + clock.Stop(); + clock.Print(); + + epsilon.DrawCopy("same"); + gPad->SetLogx(); + gPad->SetLogy(); + */ +} + +void DrawILPhoton(Double_t z=0) +{ + if (!gPad) + new TCanvas("IL_Photon", "Mean Interaction Length Photon"); + else + gPad->GetVirtCanvas()->cd(4); + + TF1 f1("length", MPhoton::InteractionLength, 1e4, 1e11, 1); + f1.SetParameter(0, z); + + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + f1.SetMaximum(1e5); + f1.SetLineWidth(1); + + TH1 *h=f1.DrawCopy()->GetHistogram(); + + h->SetTitle("Mean Interaction Length (Photon)"); + h->SetXTitle("E [GeV]"); + h->SetYTitle("x [kpc]"); + + gPad->Modified(); + gPad->Update(); +} + +void DrawILElectron(Double_t z=0) +{ + if (!gPad) + new TCanvas("IL_Electron", "Mean Interaction Length Electron"); + else + gPad->GetVirtCanvas()->cd(3); + + TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0); + f2.SetParameter(0, z); + + gPad->SetLogx(); + gPad->SetLogy(); + gPad->SetGrid(); + f2.SetLineWidth(1); + + TH1 *h=f2.DrawCopy()->GetHistogram(); + + h->SetTitle("Mean Interaction Length (Electron)"); + h->SetXTitle("E [GeV]"); + h->SetYTitle("x [kpc]"); + + gPad->Modified(); + gPad->Update(); +} + +void DrawRZ() +{ + new TCanvas("RZ", "r and z"); + + TF1 f1("ZofR", ZofR, 0, 4.5e5, 0); + TF1 f2("RofZ", RofZ, 0, 5, 0); + + gPad->Divide(2,2); + + gPad->GetVirtCanvas()->cd(1); + TH1 *h = f1.DrawCopy()->GetHistogram(); + + h->SetTitle("z(r)"); + h->SetXTitle("r [kpc]"); + h->SetYTitle("z"); + + gPad->Modified(); + gPad->Update(); + + gPad->GetVirtCanvas()->cd(2); + h = f2.DrawCopy()->GetHistogram(); + + h->SetTitle("r(z)"); + h->SetXTitle("z"); + h->SetYTitle("r [kpc]"); + + gPad->Modified(); + gPad->Update(); +} + +Double_t PrimSpect(Double_t *x, Double_t *k) +{ + return pow(pow(10, x[0]), k[0]); +} + +/* +Double_t Planck(Double_t *x, Double_t *k=NULL) +{ + Double_t E = x[0]; // [GeV] + Double_t z = k ? k[0] : 0; + + Double_t T = 2.96*(z+1); // [K] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + + Double_t EkT = E/kB/T; + + return E*E*E / (exp(EkT)-1.); // [GeV^2] +} +*/ + +Double_t Planck(Double_t *x, Double_t *k=NULL) +{ + // + // Planck, per unit volume, per unit energy + // + // constants moved out of function, see below + // + + // + // This is to get a correct random value in a TF1! + // + Double_t E = pow(10, x[0]); // [GeV] + Double_t z = k ? k[0] : 0; + + Double_t T = 2.96*(z+1); // [K] + Double_t e = 1.602176462e-19; // [C] + Double_t kB = 1e-9/e*1.3806503e-23; // [GeV/K] + + Double_t EkT = E/kB/T; + + return E*E / (exp(EkT)-1.); // [GeV^2] +} + +void DoIt() +{ + Double_t rcygnusx3 = 100; //30/3.258; // [~10kpc] + + Double_t startz = ZofR(&rcygnusx3); + + static TRandom rand(0); + MPairProduction pair; + + //Double_t nphot = 50; + Double_t runtime = 18*60*60; // [s] + + Double_t lo = 2e4; + Double_t hi = 2e10; + Double_t alpha = -2; + + Double_t nbins = 4*log10(hi/lo); + + TF1 phot("PhotonSpectrum", Planck, -15, -10.5, 1); + TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV + + src.SetParameter(0, alpha); + + TList listg; + TList liste; + + listg.SetOwner(); + liste.SetOwner(); + + gStyle->SetOptStat(10); + + TH1D hist; + TH1D histsrc; + + // + // Don't change the order!!! + // + hist.SetFillStyle(0); + hist.SetMarkerStyle(2); + histsrc.SetFillStyle(0); + histsrc.SetMarkerStyle(2); + + MBinning bins; + bins.SetEdgesLog(nbins, lo, hi); + + MH::SetBinning(&hist, &bins); + MH::SetBinning(&histsrc, &bins); + + TH1D hista; + MBinning binsa; + binsa.SetEdgesLog(16, 1e-12, 1e-8); + MH::SetBinning(&hista, &binsa); + + MH::MakeDefCanvas(); + + gPad->SetGrid(); + gPad->SetLogx(); + gPad->SetLogy(); + + histsrc.Draw("P"); + hist.Draw("Psame"); + histsrc.Draw("Csame"); + hist.Draw("Csame"); + + histsrc.SetName("Spectrum"); + histsrc.SetXTitle("E [GeV]"); + histsrc.SetYTitle("E*Counts"); + histsrc.GetXaxis()->SetLabelOffset(-0.015); + histsrc.GetXaxis()->SetTitleOffset(1.1); + + ofstream fsrc("source.dat"); + ofstream fcas("cascade.dat"); + + // ------------------------------ + + TStopwatch clock; + clock.Start(); + + Int_t timer = 0; + + Int_t n=0; + Double_t starttime = TStopwatch::GetRealTime(); // s + while (TStopwatch::GetRealTime() " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush; + + MPhoton *gamma=new MPhoton(E, startz); + listg.Add(gamma); + + while (1) + { + cout << " |P:" << flush; + + TIter NextP(&listg); + MPhoton *p = NULL; + while ((p=(MPhoton*)NextP())) + { + // + // Calculate way until interaction takes place + // + Double_t z = p->GetZ(); + Double_t l = rand.Exp(p->GetInteractionLength()); + Double_t r = RofZ(&z); + + // + // If photon went along us... (l==0: infinite) + // + if (r":"!") << flush; + hist.Fill(p->GetEnergy(), p->GetEnergy()); + fcas << p->GetEnergy() << endl; + listg.Remove(p); + continue; + } + + // + // Set new z value + // + r -= l; + z = ZofR(&r); + p->SetZ(z); + + // + // Sample phtoton from background and interaction angle + // + MPhoton photon; + + phot.SetParameter(0, z); + photon.SetEnergy(pow(10, phot.GetRandom())); + photon.SetAngle(rand.Uniform(TMath::Pi() * 2)); + + if (!pair.Process(p, &photon, &liste)) + { + cout << "0" << flush;; + continue; + } + + cout << "." << flush; + + listg.Remove(p); + } + + if (liste.GetSize()==0 && listg.GetSize()==0) + break; + + cout << " |E" << flush; + + TIter Next(&liste); + MElectron *e = NULL; + while ((e=(MElectron*)Next())) + { + cout << ":" << flush; + + hista.Fill(fabs(e->GetAngle())); + + while(1) + { + Double_t E = e->GetEnergy(); + Double_t z = e->GetZ(); + Double_t l = rand.Exp(e->GetInteractionLength()); + Double_t r = RofZ(&z) - l; + + if (r<0) + { + cout << "!" << flush; + break; + } + + z = ZofR(&r); + e->SetZ(z); + + Double_t ep = e->GetEnergyLoss(); + + if (E-epSetEnergy(E-ep); + + if (epModified(); + gPad->Update(); + timer = now; + } + + } + cout << endl; + + clock.Stop(); + clock.Print(); + + cout << "Created " << n << " gammas from source with E^" << alpha << endl; + cout << "Processed " << Form("%.1f", n/(TStopwatch::GetRealTime()-starttime)) << " gammas/sec." << endl; + + // ------------------------------ + + gPad->Clear(); + + histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n)); + + hist.DrawCopy("P"); + histsrc.DrawCopy("Psame"); + hist.DrawCopy("Csame"); + histsrc.DrawCopy("Csame"); + + gPad->SetGrid(); + gPad->SetLogx(); + gPad->SetLogy(); + + MH::MakeDefCanvas(); + hista.SetXTitle("[rad]"); + hista.DrawCopy(); + gPad->SetLogx(); + +} + +// ------------------------------------------------------------------- +void phys() +{ + /* + Double_t Eg = 1e5; + + Double_t E0 = 511e-6; // [GeV] + Double_t lolim = E0*E0/Eg; + Double_t inf = 4e-12; //E0*E0/Eg * sqrt(Eg); + + // 1e5: 5e-13, 4e-12 + // 1e6: 5e-14, 2e-12 + // 1e8: 5e-16, 1e-12 + // 1e10: 5e-18, 1e-12 + + // 1e5: 2.6e-12, 4e-12 + // 1e6: 2.6e-13, 2e-12 + // 1e8: 2.6e-15, 1e-12 + // 1e10: 1.6e-17, 1e-12 + + TF1 f("int2", MPhoton::Int2, lolim, inf, 2); + f.SetParameter(0, Eg); + f.SetParameter(1, 0); + + MH::MakeDefCanvas(); + gPad->SetLogx(); + f.DrawCopy(); + return; + */ +// MH::MakeDefCanvas(); +// DrawILPhoton(); + +// return ; + + DoIt(); +return; + EnergyLossRate(); + + MH::MakeDefCanvas(); + DrawILElectron(); + + DrawRZ(); + + return; + + /* + cout << "Starting..." << endl; + + MEvtLoop magic; + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + MHInput *input = new MHInput; + MHOutput *output = new MHOutput; + + plist.AddToList(input); + plist.AddToList(output); + + MGenerateInput geninp; + tlist.AddToList(&geninp); + + MFillH finput("MSinglePairInput", input); + tlist.AddToList(&finput); + + MSinglePairProduction prod; + tlist.AddToList(&prod); + + MFillH foutput("MSinglePairOutput", output); + tlist.AddToList(&foutput); + + magic.SetParList(&plist); + if (!magic.Eventloop(10000)) + return; + + input->Draw(); + output->Draw(); + + return; +*/ + + /* + MEvtLoop magic; + MParList plist; + MTaskList tlist; + plist.AddToList(&tlist); + + MHInput *input = new MHInput; + MHOutput *output = new MHOutput; + + plist.AddToList(input); + plist.AddToList(output); + + MGenerateInput geninp; + tlist.AddToList(&geninp); + + MFillH finput("MSinglePairInput", input); + tlist.AddToList(&finput); + + MSinglePairProduction prod; + tlist.AddToList(&prod); + + MFillH foutput("MSinglePairOutput", output); + tlist.AddToList(&foutput); + + magic.SetParList(&plist); + if (!magic.Eventloop(10000)) + return; + + input->Draw(); + output->Draw(); + + return; + */ +} +