Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 1447) +++ /trunk/MagicSoft/Mars/Changelog (revision 1448) @@ -6,4 +6,6 @@ - Fix: removed default arguments of ScaleAxis(...) (did not compile on alphas). + + 2002/07/25: Abelardo Moralejo @@ -16,4 +18,6 @@ * macros/MagicHillas.C: - changed accordingly + + 2002/07/25: Wolfgang Wittek, Thomas Bretz Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1447) +++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog (revision 1448) @@ -1,3 +1,24 @@ -*-*- END -*-*- + 2002/06/26: Thomas Bretz + + * MParticle.[h,cc]: + - added fX data member + - corrected a small bug in the SetNewPosition calculation + (dr*cos(fTheta)>R*cos(fTheta) --> x(2)>R*cos(fTheta)) + - added usage of fX to the =-operator and SetNewPosition + - added fSrcR (unsaved) + + * MElectron.cc: + - added usage of fX to the =-operator and SetNewPosition + - changed to use the more accurate SetNewPosition in case B==0 + + * MPhoton.cc: + - return 1e50 as InteractionLength in case Eg<100GeV + + * phys.C: + - added support for a B-Field-Bubble + + + 2002/06/21: Thomas Bretz Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1447) +++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (revision 1448) @@ -371,4 +371,7 @@ Bool_t MElectron::SetNewPositionB(Double_t B) { + if (B==0) + return SetNewPosition(); + static TRandom rand(0); Double_t x = rand.Exp(GetInteractionLength()); @@ -456,5 +459,5 @@ p(0) = beta_p/beta*x; p(1) = GetPType()==kEElectron?r*sin(rho):-r*sin(rho); - p(2) = r*(1-cos(rho)); + p(2) = r*(1.-cos(rho)); } else @@ -463,4 +466,5 @@ p(1) = beta_o/beta*x; p(2) = 0; + cout << "------------- HEY! --------------" << endl; } @@ -474,4 +478,11 @@ p *= N2; p *= M2; + + if (p(2) R*cos(fTheta)) + const Double_t dx = R - dr*x(2); + + if (dx < 0) { - x *= R/dr; + dr = R/x(2); // R>0 --> x(2)>0 rc = kFALSE; } + else + fX += dr*(1.-x(2)); + + x *= dr; // ------------------------------ @@ -283,7 +287,9 @@ r -= x; - fR = sqrt(r(0)*r(0)+r(1)*r(1)); - fPhi = atan2(r(1), r(0)); - fZ = ZofR(&r(2)); + if (fR!=0) + fPhi = atan2(r(1), r(0)); + + fR = sqrt(r(0)*r(0)+r(1)*r(1)); + fZ = ZofR(&r(2)); return rc; Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1447) +++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h (revision 1448) @@ -22,4 +22,6 @@ protected: + Double_t fSrcR; //! [kpc] Distance from observer to source + Double_t fEnergy; // [GeV] Energy @@ -31,5 +33,5 @@ Double_t fPsi; // [rad] Direction of momentum, az. angle - Double_t fX; // [kpc] real traveled distance (added in version 2) + Double_t fX; // [kpc] real traveled distance minus observers distance (added in version 2) public: @@ -43,4 +45,5 @@ void operator=(MParticle &p) { + fSrcR = p.fSrcR; fZ = p.fZ; fR = p.fR; @@ -60,4 +63,6 @@ void SetIsPrimary(Bool_t is=kTRUE) { is ? SetBit(kEIsPrimary) : ResetBit(kEIsPrimary); } Bool_t IsPrimary() const { return TestBit(kEIsPrimary); } + + void SetSrcR(Double_t r) { fSrcR = r; } // ---------------------------------------------------------------- Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1447) +++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc (revision 1448) @@ -121,4 +121,7 @@ Double_t z = k ? k[0] : 0; + if (Eg<100) + return 1e50; + Double_t val[2] = { Eg, z }; Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C =================================================================== --- /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1447) +++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C (revision 1448) @@ -33,5 +33,5 @@ Double_t PhotonSpect(Double_t *x, Double_t *k=NULL) { - Double_t Ep = pow(10, x[0]); + Double_t Ep = pow(10, x[0]); Double_t res = MPhoton::Int2(&Ep, k); @@ -137,19 +137,23 @@ Double_t R = MParticle::RofZ(&startz); // [kpc] - const char *filename = "cascade_0.03_24_1e2_1e6_256_1.root"; - - const Double_t B = 0; // [T] mean magnetic field - - Double_t runtime = 3.5*60*60; // [s] maximum time to run the simulation - - Int_t nbins = 24; // number of bins produced in energy spectrum - - Double_t lo = 1e2; // lower limit of displayed spectrum - Double_t hi = 1e6; // upper limit of spectrum (cutoff) - - Int_t counter = 256; // maximum number of inv. Compton (-1 means infinite) - - Double_t alpha = -1; // -1 means a E^2 plot - Double_t plot = 1; // 1 means a E^-2 spectrum + const char *filename = "cascade_0.03_18_1e2_1e5_B0_256_1.root"; + //const char *filename = "test.root"; + + const Double_t B = 0;//1e-6*1e-4; // [T=1e4G] mean magnetic field + + Double_t runtime = 45*60; // [s] maximum time to run the simulation + + Int_t nbins = 18; // number of bins produced in energy spectrum + + Double_t lo = 1e2; // lower limit of displayed spectrum + Double_t hi = 1e5; // upper limit of spectrum (cutoff) + + Int_t counter = 256; // maximum number of inv. Compton (-1 means infinite) + + Double_t alpha = -1; // -1 means a E^2 plot + Double_t plot = 1; // 1 means a E^-2 spectrum + + Double_t bubbler = R-1e3*10; // [kpc] + Double_t bubblez = MParticle::ZofR(&bubbler); // -------------------------------------------------------------------- @@ -250,4 +254,5 @@ MPhoton *gamma=new MPhoton(E, startz); gamma->SetIsPrimary(); + gamma->SetSrcR(R); gamma->InitRandom(); listg.Add(gamma); @@ -304,7 +309,6 @@ // phot.SetParameter(1, p->GetZ()); - - Double_t Ep = pow(10, phot.GetRandom()); - if (Ep==pow(10,0)) + Double_t pe = phot.GetRandom(); + if (pe==0) { cout << "z" << flush; @@ -312,4 +316,5 @@ } + Double_t Ep = pow(10, pe); Double_t theta = RandomTheta(Eg, Ep); if (theta==0) @@ -355,6 +360,5 @@ while (test<0 ? true : (test++SetNewPositionB(B)) + if (!e->SetNewPositionB(e->GetZ()>bubblez ? B : 0)) { cout << "!" << flush; @@ -423,4 +427,5 @@ file.Write(); + cout << "Wrote: " << filename << endl; cout << "Created " << n << " gammas from source with E^" << alpha << endl; cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma." << endl; @@ -441,4 +446,25 @@ void phys() { + /* + Double_t Eg = 10; + Double_t E0 = 511e-6; + Double_t z = 0.03; + Double_t lolim = E0*E0/Eg; + Double_t inf = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1)); + if (Eg<5e4) + inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2); + TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2); + phot.SetParameter(0, Eg); + phot.SetParameter(1, z); + + Double_t val[2] = {Eg,z}; + cout << phot.Integral(log10(lolim), log10(inf), val, 1e-2) << endl; + + cout << lolim << " " << inf << endl; + + phot.DrawCopy(); + return; +*/ + DoIt();