Changeset 1429 for trunk/WuerzburgSoft/Thomas

Timestamp:

07/24/02 09:32:37 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (modified) (1 diff)
• phys.C (modified) (10 diffs)

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -1 +1 @@
                                                                   -*-*- END -*-*-
+ 2002/06/21: Thomas Bretz
+
+   * MElectron.cc:
+     - changed upper limit in InteractionLength from 1e-11 to 2e-11
+     - changed upper limit of p_e in EnergyLoss from -10.8 to -10.6
+
+   * MParticle.[h,cc]:
+     - replaced Planck by tanjas background
+     
+   * background.txt:
+     - added
+
+   * MPhoton.[h,cc]:
+     - changed lower limit of Int2-Integration for Interaction Length
+       in case Eg<5e4 to fit better its needs
+     - changed max/min of drawing interaction length
+     - changed lower limit of interaction length plot
+
+   * analp.C:
+     - added plotting 1h1426 data
+     - changed alpha/plot from -2/2 to -1/1 (which means I do the well
+       known E^2 plot of an E^-2 spectrum)
+     - added a cutoff to simulate cutoff spectra
+     - some small corrections for the layout of the plots
+
+   * phys.C:
+     - added a maximum number of inv. Compton processes
+     - don't throw away photons which were reemitted lower than the lower
+       limit of the plot (so they are still written to the file)
+     - changed the lower limit of the PhotonSpect simulation function
+       to a physically more realistic value. No photon which don't have
+       an corresponding angle are produced anymore.
+     - stop the electron only if its energy is less than a thousand
+       of its primary energy
+
+
+
  2002/06/21: Thomas Bretz
 

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -34 +34 @@
 {
     Double_t Ep  = pow(10, x[0]);
+
     Double_t res = MPhoton::Int2(&Ep, k);
-
     return res*1e55; //65/k[0];
-    // return MPhoton::Planck(&Ep, &k[1]);
+
+    //return MPhoton::Planck(&Ep, &k[1]);
 }
 
@@ -119 +120 @@
 void DoIt()
 {
-    Double_t R      = 100; //MParticle::RofZ(&startz); // [kpc]
-    Double_t startz = MParticle::ZofR(&R);
-
-    const char *filename = "cascade_delme.root";
-
-    const Double_t B = 0;
+    // a -->  5
+    // b --> 10
+    // c -->  2
+
+    // delme, delmeB starting integration at lolim
+    // delme2, delme2B starting integration at -log10(Eg)-8
+    // delme3,  lolim, 24, 1e2-1e6,   2x inv.Compton
+    // delme3B, lolim, 24, 1e2-1e6,   4x inv.Compton
+    // delme3C, lolim, 24, 1e2-1e6,   8x inv.Compton
+    // delme3D, lolim, 24, 1e2-1e6,  16x inv.Compton
+    // delme3E, lolim, 24, 1e2-1e6,  32x inv.Compton
+
+    // delme3H, lolim, 24, 1e2-1e6, 256x inv.Compton 8*60
+
+    Double_t startz = 0.03; //MParticle::ZofR(&R);
+    Double_t R      = MParticle::RofZ(&startz); // [kpc]
+
+    const char *filename = "delme3H.root";
+
+    const Double_t B = 0; // [T] mean magnetic field
+
+    Double_t runtime = 8*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
+
+    // --------------------------------------------------------------------
 
     cout << "R = " << R << "kpc" << endl;
@@ -131 +160 @@
     MPairProduction pair;
 
-    Double_t runtime = 5*60; // [s]
-
-    Double_t lo = 1e4;
-    Double_t hi = 1e11;
-    Double_t alpha = -2;
-    Double_t plot  =  2;
-
-    Int_t nbins = 21; //4*log10(hi/lo);
-
     TH1D hist;
     TH1D histsrc;
 
-    hist.SetMinimum(pow(lo, plot+alpha)/10);
-    histsrc.SetMinimum(pow(lo, plot+alpha)/10);
+    hist.SetMinimum(pow(lo, plot+alpha)/100);
+    histsrc.SetMinimum(pow(lo, plot+alpha)/100);
 
     TList listg;
@@ -230 +250 @@
         {
             if (listg.GetSize()!=0)
-                cout << " |P:" << flush;
+                cout << " |P" << flush;
 
             TIter NextP(&listg);
@@ -237 +257 @@
             while ((p=(MPhoton*)NextP()))
             {
+                cout << ":" << flush;
                 Double_t Eg = p->GetEnergy();
 
-                if (!p->SetNewPosition())
+                Double_t E0 = 511e-6;
+                Double_t z  = p->GetZ();
+                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(Eg)-8, -10.5, 2);
+                //Double_t E0 = 511e-6; // [GeV]
+                TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2);
+                phot.SetParameter(0, Eg);
+                while (1)
                 {
-                    cout << "!" << flush;
-
-                    hist.Fill(Eg, pow(Eg, plot)*weight);
-
-                    p->SetIsPrimary(kFALSE);
-                    T1->Fill();
-
-                    Esum += Eg;
-
-                    delete listg.Remove(p);
-                    continue;
+                    if (!p->SetNewPosition())
+                    {
+                        cout << "!" << flush;
+
+                        hist.Fill(Eg, pow(Eg, plot)*weight);
+
+                        p->SetIsPrimary(kFALSE);
+                        T1->Fill();
+
+                        Esum += Eg;
+
+                        break;
+                    }
+
+                    //
+                    // Sample phtoton from background and interaction angle
+                    //
+                    phot.SetParameter(1, p->GetZ());
+
+                    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)
+                    {
+                        cout << "t" << flush;
+                        continue;
+                    }
+
+                    if (!pair.Process(p, Ep, theta, &liste))
+                    {
+                        cout << "0" << flush;
+                        continue;
+                    }
+
+                    cout << "." << flush;
+                    break;
                 }
-
-                //
-                // Sample phtoton from background and interaction angle
-                //
-                TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
-                phot.SetParameter(0, Eg);
-                phot.SetParameter(1, p->GetZ());
-
-                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)
-                {
-                    cout << "t" << flush;
-                    continue;
-                }
-
-                if (!pair.Process(p, Ep, theta, &liste))
-                {
-                    cout << "0" << flush;
-                    continue;
-                }
-
                 delete listg.Remove(p);
-                cout << "." << flush;
             }
 
@@ -299 +333 @@
                 e->SetIsPrimary(kFALSE);
 
-                if (e->GetEnergy()<lo)
+                Double_t Ee = e->GetEnergy();
+
+                if (Ee<lo)
                     continue;
 
                 cout << ":" << flush;
-                while(1)
+                int test = counter<0 ? -1 : counter;
+                while (test<0 ? true : (test++<256))
                 {
+
                     if (!e->SetNewPositionB(B))
                     {
@@ -316 +354 @@
                     T2->Fill();
 
+                    cout << "." << flush;
+                    listg.Add(p);
+                    /*
                     if (fabs(p->GetTheta()*3437)<60 &&  // < 60min
                         p->GetEnergy()>lo)
@@ -324 +365 @@
                     else
                     {
-                        if (p->GetEnergy()<=lo)
+                        if (p->GetEnergy()<E*pow(10, alpha)/5 || p->GetEnergy()<=lo)
                             cout << "e" << flush;
                         else
@@ -332 +373 @@
                         Emis += p->GetEnergy();
                     }
-
+                    */
                     if (fabs(e->GetTheta()*3437)>60)  // < 60min
                     {
@@ -339 +380 @@
                     }
 
-                    if (e->GetEnergy()<lo)//*5)
+                    if (e->GetEnergy()<Ee*1e-3) // <2e3
                     {
                         cout << "E" << flush;