Changeset 1428 for trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

Timestamp:

07/24/02 09:17:52 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (modified) (1 diff)

trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

@@ -72 +72 @@
 }
 
-Double_t MParticle::Planck(Double_t *x, Double_t *k)
+#include <fstream.h>
+#include <TH2.h>
+#include "../mhist/MBinning.h"
+#include "../mhist/MH.h"
+
+TH2D *hist2;
+
+Double_t MParticle::Planck(Double_t *x, Double_t *k=NULL)
+{
+    static Bool_t isloaded = kFALSE;
+
+    if (!isloaded)
+    {
+        Double_t c   = 299792458;             // [m/s]
+        Double_t e   = 1.602176462e-19;       // [C]
+        Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
+        Double_t hc  = h*c;                   // [GeVm]
+        Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
+
+        ifstream fin("background.txt");
+
+        hist2 = new TH2D;
+
+        MBinning binsz;
+        MBinning binse;
+        binsz.SetEdgesLog(100, 1e-6,  1);    // --> 101 Edges / 100 bins
+        binse.SetEdgesLog(100, 7e-15, 3e-8); // --> 101 Edges / 100 bins
+
+        MH::SetBinning(hist2, &binsz, &binse);
+
+        for (int y=0; y<101; y++)
+        {
+            Double_t val;
+            fin >> val;
+            for (int x=0; x<101; x++)
+            {
+                fin >> val;
+
+                val += 9;
+
+                Double_t z = binsz.GetEdges()[x];
+                Double_t f = z+1;
+
+                Double_t newval = pow(10, val)/konst;
+                hist2->SetBinContent(x, y, newval*f*f*f);
+
+            }
+        }
+        isloaded = kTRUE;
+    }
+
+    //
+    // y = (y1-y0)/(x1-x0) * (x-x0) + y0
+    //
+    Double_t z = k ? k[0] : 0;
+    Double_t E = x[0];
+
+    Int_t i = hist2->GetXaxis()->FindFixBin(z);
+    Int_t j = hist2->GetYaxis()->FindFixBin(E);
+
+    Double_t z1 = hist2->GetXaxis()->GetBinLowEdge(i+1);
+    Double_t z0 = hist2->GetXaxis()->GetBinLowEdge(i);
+
+    Double_t E1 = hist2->GetYaxis()->GetBinLowEdge(j+1);
+    Double_t E0 = hist2->GetYaxis()->GetBinLowEdge(j);
+
+    Double_t n00 = hist2->GetBinContent(i,   j);
+    Double_t n01 = hist2->GetBinContent(i+1, j);
+    Double_t n10 = hist2->GetBinContent(i,   j+1);
+    Double_t n11 = hist2->GetBinContent(i+1, j+1);
+
+    Double_t dz = (z-z0)/(z1-z0);
+    Double_t dE = (E-E0)/(E1-E0);
+
+    Double_t n0 = dz*(n01-n00)+n00;
+    Double_t n1 = dz*(n11-n10)+n10;
+
+    Double_t n = dE*(n1-n0)+n0;
+
+    return n;
+    /*
+     //
+     // TANJA2
+     //
+     Double_t E1 = x[0];
+     Double_t E2 = x[0]/8;
+     return (MParticle::Planck0(&E1, k)+MParticle::Planck0(&E2, k)/40e3)*0.7/0.4;
+     */
+    /*
+     //
+     // TANJA
+     //
+     Double_t E1 = x[0];
+     Double_t E2 = x[0]/8;
+     return Planck0(&E1, k)+Planck0(&E2, k)/5e3;
+     */
+}
+
+Double_t MParticle::Planck0(Double_t *x, Double_t *k)
 {
     //