Changeset 1367 for trunk/WuerzburgSoft/Thomas/mphys

Timestamp:

06/14/02 10:41:23 (23 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (modified) (2 diffs)
• MElectron.cc (modified) (4 diffs)
• MPairProduction.cc (modified) (1 diff)
• MPhoton.cc (modified) (2 diffs)
• phys.C (modified) (7 diffs)

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -8 +8 @@
      - changed the eps of all integrals to 1e-2
      - changed the p_e function to use the Compton function
+     - made all function which are used for integration only static
+       and changed the function definition from=0, to=0
 
    * MPhoton.[h,cc]:
@@ -25 +27 @@
    * phys.C:
      - changed the integral eps to 1e-2
+     - small changed to the output
+     - removed integral==0, replaced by Ep==pow(10,0)
      
 

trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

@@ -90 +90 @@
         return DiSum(x);
 
-    TF1 IntLi("Li", Li, 0, z, 0);
+    // TF1 IntLi("Li", Li, 0, z, 0);
+    static TF1 IntLi("Li", Li, 0, 0, 0);
     const Double_t integ = IntLi.Integral(0, z, (Double_t*)NULL, 1e-2);
     return -integ;
@@ -160 +161 @@
     const Double_t z      = k ? k[0] : 0;
 
-    Double_t val[3] = { E[0], z };                         // E[GeV]
-
-    const Double_t from   = 1e-17;
-    const Double_t to     = 1e-11;
-
     /* -------------- old ----------------
        Double_t from   = 1e-15;
@@ -171 +167 @@
        -----------------------------------
     */
-    TF1 func("Compton", Compton, from, to, 2);            // [0, inf]
+    static TF1 func("Compton", Compton, 0, 0, 2);         // [0, inf]
+
+    const Double_t from   = 1e-17;
+    const Double_t to     = 1e-11;
+
+    Double_t val[3] = { E[0], z };                        // E[GeV]
 
     Double_t integ  = func.Integral(from, to, val, 1e-2); // [Gev] [0, inf]
@@ -260 +261 @@
     const Double_t Gamma = 4.*omega;
 
-    // --old-- fQ.SetRange(1e-6, 1./(1.+ 2.*Gamma));
     fQ.SetParameter(0, Gamma);
 

trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
-!   Author(s): Thomas Bretz  12/2000 (tbretz@uni-sw.gwdg.de)
+!   Author(s): Thomas Bretz  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2001
+!   Copyright: MAGIC Software Development, 2000-2002
 !
 !

trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc

@@ -92 +92 @@
     Double_t val[2] = { Eg, Ep };
 
-    const Double_t from = -1.0;
-    const Double_t to   = 1.-E0*E0/(2.*Eg*Ep); // Originally Was: 1.
+    static TF1 f("int1", Int1, 0, 0, 2);
 
-    TF1 f("int1", Int1, from, to, 2);
-
+    const Double_t from   = -1.0;
+    const Double_t to     = 1.-E0*E0/(2.*Eg*Ep);              // Originally Was: 1.
     const Double_t int1   = f.Integral(from, to, val, 1e-2);  // [m^2]
-    const Double_t planck = MParticle::Planck(&Ep, &z); // [GeV^2]
+    const Double_t planck = MParticle::Planck(&Ep, &z);       // [GeV^2]
 
     const Double_t res = planck * int1;
@@ -124 +123 @@
     Double_t val[2] = { Eg, z };
 
+    static TF1 f("int2", Int2, 0, 0, 2);
+
     Double_t lolim = E0*E0/Eg;
     Double_t inf   = Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1);
-
-    TF1 f("int2", Int2, lolim, inf, 2);
-
-    Double_t int2 = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2]
+    Double_t int2  = f.Integral(lolim, inf, val, 1e-2); //[GeV^3 m^2]
 
     if (int2==0)

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -61 +61 @@
     TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0);
 
+    func.SetNpx(50);
+
     Double_t sbar  = pow(10, func.GetRandom());
     Double_t theta = acos(1.-sbar*2/f);
@@ -116 +118 @@
 void DoIt()
 {
-    Double_t startz = 0.003;                   //MParticle::ZofR(&R);
-    Double_t R      = MParticle::RofZ(&startz); // [kpc]
+    Double_t R      = 100; //MParticle::RofZ(&startz); // [kpc]
+    Double_t startz = MParticle::ZofR(&R);
 
     cout << "R = " << R << "kpc" << endl;
@@ -127 +129 @@
     Double_t runtime = 15*60; // [s]
 
-    Double_t lo = 1e4; 
+    Double_t lo = 1e4;
     Double_t hi = 1e11;
     Double_t alpha = -2;
@@ -274 +276 @@
                 phot.SetParameter(0, Eg);
                 phot.SetParameter(1, p->GetZ());
+                /*
                 if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0)
                 {
@@ -279 +282 @@
                     continue;
                 }
-
+                */
                 Double_t Ep    = pow(10, phot.GetRandom());
+                if (Ep==pow(10,0))
+                {
+                    cout << "z" << flush;
+                    continue;
+                }
                 Double_t theta = RandomTheta(Eg, Ep);
                 if (theta==0)
@@ -339 +347 @@
                     else
                     {
-                        cout << "i" << flush; // ignored
+                        if (p->GetEnergy()<=lo)
+                            cout << "e" << flush;
+                        else
+                            cout << "t" << flush; // ignored
                         delete p;
                     }
@@ -374 +385 @@
 
     cout << "Created " << n << " gammas from source with E^" << alpha << endl;
-    cout << "Processed " << Form("%.1f", n/(TStopwatch::GetRealTime()-starttime)) << " gammas/sec." << endl;
+    cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma." << endl;
 
     // ------------------------------