Changeset 1448 for trunk/WuerzburgSoft

Timestamp:

07/26/02 14:49:58 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (modified) (1 diff)
• MElectron.cc (modified) (5 diffs)
• MParticle.cc (modified) (3 diffs)
• MParticle.h (modified) (4 diffs)
• MPhoton.cc (modified) (1 diff)
• phys.C (modified) (8 diffs)

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -1 +1 @@
                                                                   -*-*- 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
 

trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

@@ -371 +371 @@
 Bool_t MElectron::SetNewPositionB(Double_t B)
 {
+    if (B==0)
+        return SetNewPosition();
+
     static TRandom rand(0);
     Double_t x = rand.Exp(GetInteractionLength());
@@ -456 +459 @@
         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 +466 @@
         p(1) = beta_o/beta*x;
         p(2) = 0;
+        cout << "------------- HEY! --------------" << endl;
     }
 
@@ -474 +478 @@
     p *= N2;
     p *= M2;
+
+    if (p(2)<x) // happens sometimes in case B==0
+    {
+        cout << "----- HA: " << B << " " << x << " " << p(2) << " " << x-p(2) << endl;
+        p(2)=x;
+    }
+
     s -= p;
 
@@ -479 +490 @@
     fPhi = atan2(s(1), s(0));
     fZ   = ZofR(&s(2));
+    fX  += x-p(2);
 
     // -----------------------------

trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

@@ -198 +198 @@
 
 MParticle::MParticle(ParticleType_t t, const char *name, const char *title)
-    : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0)
+    : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0), fX(0)
 {
     //
@@ -259 +259 @@
     x(2) = cos(fTheta);
 
-    x *= dr;
-
     // ------------------------------
 
     const Double_t R = RofZ(&fZ);
 
-    if (x(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 +287 @@
     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;

trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

@@ -22 +22 @@
 
 protected:
+    Double_t fSrcR;   //! [kpc] Distance from observer to source
+
     Double_t fEnergy; // [GeV] Energy
 
@@ -31 +33 @@
     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 +45 @@
     void operator=(MParticle &p)
     {
+        fSrcR = p.fSrcR;
         fZ = p.fZ;
         fR = p.fR;
@@ -60 +63 @@
     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; }
 
     // ----------------------------------------------------------------

trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc

@@ -121 +121 @@
     Double_t z  = k ? k[0] : 0;
 
+    if (Eg<100)
+        return 1e50;
+
     Double_t val[2] = { Eg, z };
 

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -33 +33 @@
 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 +137 @@
     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 +254 @@
         MPhoton *gamma=new MPhoton(E, startz);
         gamma->SetIsPrimary();
+        gamma->SetSrcR(R);
         gamma->InitRandom();
         listg.Add(gamma);
@@ -304 +309 @@
                     //
                     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 +316 @@
                     }
 
+                    Double_t Ep = pow(10, pe);
                     Double_t theta = RandomTheta(Eg, Ep);
                     if (theta==0)
@@ -355 +360 @@
                 while (test<0 ? true : (test++<counter))
                 {
-
-                    if (!e->SetNewPositionB(B))
+                    if (!e->SetNewPositionB(e->GetZ()>bubblez ? B : 0))
                     {
                         cout << "!" << flush;
@@ -423 +427 @@
     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 +446 @@
 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();