Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1352) +++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1352) @@ -0,0 +1,19 @@ + -*-*- END -*-*- + 2002/06/10: Thomas Bretz + + * MElectron.[h,cc], MPhoton.[h,cc], MParticle.[h,cc]: + - implemented SetNewDirection, SetNewPosition + - implemented DoInvCompton + - =operator + + * MPairProduction.cc: + - changed to use SetNewDirection Method + + * phys.C: + - changed to new style functions + + * MGeinIRPhoton.cc, MGenPrimaryParticle.[h,cc]: + - make compile (not used) + - changed Process to get angle + + Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1352) @@ -34,4 +34,6 @@ #include +#include "MPhoton.h" + ClassImp(MElectron); @@ -238,5 +240,5 @@ -Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k = NULL) +Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep) { Double_t E = x[0]; @@ -255,4 +257,7 @@ Double_t e = pow(10, fP.GetRandom()); + if (ep) + *ep = e; + Double_t omega = e*E/E0/E0; Double_t Gamma = 4.*omega; @@ -269,11 +274,55 @@ } -Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z) +Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z, Double_t *ep) { return EnergyLoss(&E, &z); } -Double_t MElectron::GetEnergyLoss() const -{ - return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ); -} +Double_t MElectron::GetEnergyLoss(Double_t *ep) const +{ + return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ, ep); +} + +#include + +MPhoton *MElectron::DoInvCompton() +{ + Double_t E0 = 511e-6; //[GeV] + + Double_t epsilon; + Double_t e = GetEnergyLoss(&epsilon); + /* + Double_t gamma = fEnergy/E0; + Double_t beta = sqrt(1-1/(gamma*gamma)); + Double_t bega = sqrt(gamma*gamma-1); + + Double_t theta = TMath::Pi()/4; + + // er: photon energy before interaction rest frame + // e: photon energy after interaction lab + Double_t er = gamma*epsilon*(1-beta*cos(theta)); // photon energy rest frame + Double_t r1 = fEnergy/e; + Double_t r2 = E0/er; + Double_t ctheta = (r1-r2-1)/(beta*r1-1); + Double_t tg = sqrt(1-ctheta*ctheta)/gamma/(beta+ctheta); + */ + fEnergy -= e; + + MPhoton &p = *new MPhoton(e, fZ); + p = *this; + + /* + cout << "(" << atan(tg)*180/TMath::Pi() << ")" << flush; + + static TRandom rand(0); + p->SetNewDirection(atan(tg), rand.Uniform(TMath::Pi()*2)); + */ + + // MISSING: Electron angle + // + // E1 = fEnergy (after!) + // + // sin(t) = (epsilon sin(theta) - e sin(atan(tg)))/sqrt(E1*E1-E0*E0) + + return &p; +} Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.h (revision 1352) @@ -5,4 +5,6 @@ #include "MParticle.h" #endif + +class MPhoton; class MElectron : public MParticle @@ -14,4 +16,6 @@ fZ = z; } + + void operator=(MParticle &p) { MParticle::operator=(p); } // ---------------------------------------------------------------- @@ -31,8 +35,10 @@ 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); + static Double_t EnergyLoss(Double_t *E, Double_t *z=NULL, Double_t *ep=NULL); + static Double_t GetEnergyLoss(Double_t E, Double_t z=0, Double_t *ep=NULL); - Double_t GetEnergyLoss() const; + Double_t GetEnergyLoss(Double_t *ep) const; + + MPhoton *DoInvCompton(); ClassDef(MElectron, 1) Index: /trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MGenIRPhoton.cc (revision 1352) @@ -34,5 +34,5 @@ #include "MParList.h" -#include "MParticle.h" +#include "MPhoton.h" ClassImp(MGenIRPhoton); @@ -104,8 +104,11 @@ { const Double_t pi2 = TMath::Pi() * 2; - MParticle *p = new MParticle(MParticle::kEGamma); - p->SetAngle(fRand.Uniform(pi2)); - p->SetEnergy(fSrc->GetRandom()*1e-9); + MPhoton *p = new MPhoton; + /* + MParticle *p = new MParticle(MParticle::kEGamma); + p->SetAngle(fRand.Uniform(pi2)); + p->SetEnergy(fSrc->GetRandom()*1e-9); + */ return p; Index: /trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MGenPrimaryParticle.cc (revision 1352) @@ -71,6 +71,6 @@ Bool_t MGenPrimaryParticle::Process() { - MParticle *p = new MParticle(MParticle::kEGamma); - p->SetEnergy(fSrc->GetRandom()); +// MParticle *p = new MParticle(MParticle::kEGamma); +// p->SetEnergy(fSrc->GetRandom()); return kTRUE; } Index: /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc (revision 1352) @@ -32,4 +32,5 @@ #include #include +#include #include "TList.h" @@ -70,5 +71,5 @@ // -------------------------------------------------------------------------- -Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, TList *list) +Bool_t MPairProduction::Process(MParticle *gamma, MParticle *phot, const Double_t theta, TList *list) { // @@ -78,5 +79,5 @@ const Double_t E0 = 511e-6; // [GeV] - const Double_t theta = phot->GetAngle(); // [2pi] + //const Double_t theta = phot->GetAngle(); // [2pi] const Double_t Ep = phot->GetEnergy(); // [GeV] const Double_t Eg = gamma->GetEnergy(); // [GeV] @@ -120,7 +121,8 @@ 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()); + MElectron &p0 = *new MElectron(E+dE, gamma->GetZ()); + MElectron &p1 = *new MElectron(E-dE, gamma->GetZ()); + p0 = *gamma; // Set Position and direction + p1 = *gamma; // Set Position and direction const Double_t E1 = E0/(E+dE); @@ -130,6 +132,6 @@ const Double_t beta2 = sqrt(1.-E2*E2); - p[0]->SetBeta(beta1); - p[1]->SetBeta(beta2); + p0.SetBeta(beta1); + p1.SetBeta(beta2); const Double_t Bscp = Bsind*cphi; @@ -140,9 +142,11 @@ const Double_t tan2 = (spg*(Bcosd-1) + Bscp)/((cpg*(Bcosd-1) - Bscp)); - p[0]->SetAngle(atan(tan1)); - p[1]->SetAngle(atan(tan2)); + static TRandom rand(0); + Double_t rnd = rand.Uniform(TMath::Pi() * 2); + p0.SetNewDirection(atan(tan1), rnd); + p1.SetNewDirection(atan(tan2), rnd+TMath::Pi()); - list->Add(p[0]); - list->Add(p[1]); + list->Add(&p0); + list->Add(&p1); return kTRUE; Index: /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.h (revision 1352) @@ -19,5 +19,5 @@ virtual ~MPairProduction(); - Bool_t Process(MParticle *g, MParticle *p, TList *l); + Bool_t Process(MParticle *g, MParticle *p, const Double_t theta, TList *l); ClassDef(MPairProduction, 0) // Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (revision 1352) @@ -7,8 +7,53 @@ #include "MParticle.h" +#include +#include + ClassImp(MParticle); +/*********************** + * + * 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) +{ + const Double_t c = 299792458; // [m/s] + const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + const Double_t r = x[0] /pc*ly*1e3; // [m] + + const 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; + + const Double_t c = 299792458; // [m/s] + const Double_t H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1] + + const Double_t ly = 3600.*24.*365.*c; // [m/ly] + const Double_t pc = 1./3.258; // [pc/ly] + + const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m] + + return R * pc/ly/1e3; // [kpc] +} + MParticle::MParticle(ParticleType_t t, const char *name, const char *title) - : fPType(t) + : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0) { // @@ -25,2 +70,87 @@ } +void MParticle::SetNewDirection(Double_t theta, Double_t phi) +{ + TMatrixD B(3, 3); + + B(0, 0) = cos(fTheta)*cos(fPsi); + B(1, 0) = cos(fTheta)*sin(fPsi); + B(2, 0) = -sin(fTheta); + + B(0, 1) = -sin(fPsi); + B(1, 1) = cos(fPsi); + B(2, 1) = 0; + + B(0, 2) = sin(fTheta)*cos(fPsi); + B(1, 2) = sin(fTheta)*sin(fPsi); + B(2, 2) = cos(fTheta); + + // ------------------------------ + + TVectorD r(3); + + r(0) = sin(theta)*cos(phi); + r(1) = sin(theta)*sin(phi); + r(2) = cos(theta); + + // ------------------------------ + + r *= B; + + fTheta = sqrt(r(0)*r(0)+r(1)*r(1)); // Numerically bad: acos(r(2)); + fPsi = atan2(r(1), r(0)); + + if (fTheta*2 > TMath::Pi()) + fTheta = fabs(fTheta-TMath::Pi()); +} + +#include +#include + +Bool_t MParticle::SetNewPosition(Double_t dr) +{ + Bool_t rc=kTRUE; + + TVectorD x(3); + + x(0) = sin(fTheta)*cos(fPsi); + x(1) = sin(fTheta)*sin(fPsi); + x(2) = cos(fTheta); + + x *= dr; + + // ------------------------------ + + const Double_t R = RofZ(&fZ); + + if (x(2) > R*cos(fTheta)) + { + x *= R/dr; + rc = kFALSE; + } + + // ------------------------------ + + TVectorD r(3); + + r(0) = fR*cos(fPhi); + r(1) = fR*sin(fPhi); + r(2) = R; + + // ------------------------------ + + r -= x; + + fR = sqrt(r(0)*r(0)+r(1)*r(1)); + fPhi = atan2(r(1), r(0)); + fZ = ZofR(&r(2)); + + return rc; +} + +Bool_t MParticle::SetNewPosition() +{ + static TRandom rand(0); + Double_t r = rand.Exp(GetInteractionLength()); + return SetNewPosition(r); +} Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1352) @@ -17,6 +17,11 @@ Double_t fEnergy; // [GeV] Energy Double_t fBeta; // [1] beta - Double_t fAngle; // [rad] Angle + Double_t fZ; // [1] Red shift + Double_t fR; // [kpc] Radius from line of view + Double_t fPhi; // [rad] Phi@z from line of view + + Double_t fTheta; // [rad] Direction of momentum, angle || line of view + Double_t fPsi; // [rad] Direction of momentum, az. angle public: @@ -24,16 +29,41 @@ // ~MParticle(); + void operator=(MParticle &p) + { + fZ = p.fZ; + fR = p.fR; + fPhi = p.fPhi; + fTheta = p.fTheta; + fPsi = p.fPsi; + } + // void Fill(const MParContainer *inp); // void Draw(Option_t *opt=NULL); + virtual Double_t GetInteractionLength() const = 0; 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; } + /* + void SetPhi(Double_t p) { fPhi = p; } + void SetR(Double_t r) { fR = r; } + void SetTheta(Double_t t) { fTheta = t; } + void SetPsi(Double_t p) { fPsi = p; } + */ + + void SetNewDirection(Double_t theta, Double_t phi); + Bool_t SetNewPosition(Double_t dr); + Bool_t SetNewPosition(); Double_t GetEnergy() const { return fEnergy; } - Double_t GetAngle() const { return fAngle; } + Double_t GetZ() const { return fZ; } + Double_t GetPhi() const { return fPhi; } + Double_t GetR() const { return fR; } + + Double_t GetTheta() const { return fTheta; } + Double_t GetPsi() const { return fPsi; } ClassDef(MParticle, 1) // Container which holds hostograms for the Hillas parameters Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h (revision 1352) @@ -15,4 +15,6 @@ } + void operator=(MParticle &p) { MParticle::operator=(p); } + static Double_t Planck(Double_t *x, Double_t *k=NULL); static Double_t Sigma_gg(Double_t *x, Double_t *k=NULL); Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1351) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1352) @@ -6,4 +6,5 @@ #include #include +#include #include #include @@ -689,5 +690,5 @@ void DoIt() { - Double_t rcygnusx3 = 100; //30/3.258; // [~10kpc] + Double_t rcygnusx3 = 5000; //100; //30/3.258; // [~10kpc] Double_t startz = ZofR(&rcygnusx3); @@ -697,5 +698,5 @@ //Double_t nphot = 50; - Double_t runtime = 18*60*60; // [s] + Double_t runtime = 15*60; ///*18*60*/60; // [s] Double_t lo = 2e4; @@ -717,4 +718,23 @@ gStyle->SetOptStat(10); + + TH2D position; + TH2D angle; + TH1D histpos; + TH1D hista; + + MBinning binspolx; + MBinning binspoly; + binspolx.SetEdges(32, -180, 180); + binspoly.SetEdgesLog(20, 1e-9, 1e-5); + MH::SetBinning(&position, &binspolx, &binspoly); + MH::SetBinning(&angle, &binspolx, &binspoly); + MH::SetBinning(&histpos, &binspoly); + MH::SetBinning(&hista, &binspoly); + + hista.SetTitle("Particle Angle"); + angle.SetTitle("Particle Angle"); + histpos.SetTitle("Particle Position"); + position.SetTitle("Particle Position"); TH1D hist; @@ -735,9 +755,4 @@ MH::SetBinning(&histsrc, &bins); - TH1D hista; - MBinning binsa; - binsa.SetEdgesLog(16, 1e-12, 1e-8); - MH::SetBinning(&hista, &binsa); - MH::MakeDefCanvas(); @@ -791,29 +806,19 @@ 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 (rSetNewPosition()) { - cout << (l==0?">":"!") << flush; + cout << "!" << flush; + hist.Fill(p->GetEnergy(), p->GetEnergy()); + position.Fill(p->GetPhi()*kRad2Deg, p->GetR()); + angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg); + histpos.Fill(p->GetR()); + hista.Fill(p->GetTheta()*kRad2Deg); + fcas << p->GetEnergy() << endl; + listg.Remove(p); continue; } - - // - // Set new z value - // - r -= l; - z = ZofR(&r); - p->SetZ(z); // @@ -822,9 +827,9 @@ MPhoton photon; - phot.SetParameter(0, z); + phot.SetParameter(0, p->GetZ()); photon.SetEnergy(pow(10, phot.GetRandom())); - photon.SetAngle(rand.Uniform(TMath::Pi() * 2)); - - if (!pair.Process(p, &photon, &liste)) + Double_t theta = rand.Uniform(TMath::Pi() * 2); + + if (!pair.Process(p, &photon, theta, &liste)) { cout << "0" << flush;; @@ -832,7 +837,6 @@ } + listg.Remove(p); cout << "." << flush; - - listg.Remove(p); } @@ -847,15 +851,7 @@ { 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) + if (!e->SetNewPosition()) { cout << "!" << flush; @@ -863,43 +859,20 @@ } - z = ZofR(&r); - e->SetZ(z); - - Double_t ep = e->GetEnergyLoss(); - - if (E-epDoInvCompton(); + if (e->GetEnergy()SetEnergy(E-ep); - - if (epGetEnergy()Modified(); gPad->Update(); @@ -943,5 +917,24 @@ MH::MakeDefCanvas(); - hista.SetXTitle("[rad]"); + position.SetXTitle("Pos: \\Phi [\\circ]"); + position.SetYTitle("Pos: R [kpc]"); + position.DrawCopy("surf1pol"); + gPad->SetLogy(); + + MH::MakeDefCanvas(); + angle.SetXTitle("Angle: \\Psi [\\circ]"); + angle.SetYTitle("Angle: \\Theta [\\circ]"); + angle.DrawCopy("surf1pol"); + gPad->SetLogy(); + + MH::MakeDefCanvas(); + histpos.SetXTitle("R [kpc]"); + histpos.SetYTitle("Counts"); + histpos.DrawCopy(); + gPad->SetLogx(); + + MH::MakeDefCanvas(); + hista.SetXTitle("\\Theta [\\circ]"); + hista.SetYTitle("Counts"); hista.DrawCopy(); gPad->SetLogx();