Changeset 1448 for trunk/WuerzburgSoft/Thomas/mphys/phys.C

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/phys.C (modified) (8 diffs)

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();