Changeset 1373

Timestamp:

06/20/02 16:24:20 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

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

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -8 +8 @@
    * MParticle.cc:
      - made TRandom static (InitRandom)
+
+   * mphys.C:
+     - added support for Emis and Esum
 
 

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -118 +118 @@
 void DoIt()
 {
-    Double_t R      = 500; //MParticle::RofZ(&startz); // [kpc]
+    Double_t R      = 100; //MParticle::RofZ(&startz); // [kpc]
     Double_t startz = MParticle::ZofR(&R);
 
-    const char *filename = "cascade_500kpc_21d_B1e-18_03.root";
-
-    const Double_t B = 1e-18;
+    const char *filename = "data/cascade_delme.root";
+
+    const Double_t B = 0;
 
     cout << "R = " << R << "kpc" << endl;
@@ -131 +131 @@
     MPairProduction pair;
 
-    Double_t runtime = 1; // [s]
+    Double_t runtime = 60; // [s]
 
     Double_t lo = 1e4;
@@ -246 +246 @@
         cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush;
 
+        Double_t Esum=0;
+        Double_t Emis=0;
+
         while (1)
         {
@@ -270 +273 @@
                     p->SetIsPrimary(kFALSE);
                     T1->Fill();
+
+                    Esum += Eg;
 
                     delete listg.Remove(p);
@@ -335 +340 @@
                     {
                         cout << "!" << flush;
+                        Esum += e->GetEnergy();
+                        Emis += e->GetEnergy();
                         break;
                     }
 
                     // WRONG!
-                    Double_t theta = rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
+                    Double_t theta = 0; //rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
                     MPhoton *p = e->DoInvCompton(theta);
 
@@ -357 +364 @@
                             cout << "t" << flush; // ignored
                         delete p;
+                        Esum += p->GetEnergy();
+                        Emis += p->GetEnergy();
                     }
 
-                    if (fabs(e->GetTheta()*3437)<60 &&  // < 60min
-                        e->GetEnergy()>5*lo)
-                        continue;
-
-                    cout << "0" << flush;
-                    break;
+                    if (fabs(e->GetTheta()*3437)>60)  // < 60min
+                    {
+                        cout << "T" << flush;
+                        Esum += e->GetEnergy();
+                        Emis += e->GetEnergy();
+                        break;
+                    }
+
+                    if (e->GetEnergy()<lo)//*5)
+                    {
+                        cout << "E" << flush;
+                        Esum += e->GetEnergy();
+                        Emis += e->GetEnergy();
+                        break;
+                    }
                 }
             }
             liste.Delete();
         }
-        cout << endl;
+        cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl;
 
         Int_t now = (int)(TStopwatch::GetRealTime()- starttime)/5;