From 81e9fad06ffc14370a9063a77b88f27b50d08b63 Mon Sep 17 00:00:00 2001
From: Andreas Suter <andreas.suter@psi.ch>
Date: Tue, 5 Nov 2013 20:35:47 +0000
Subject: [PATCH] slight improvement in the Poisson noise, still not perfect.

---
 .../analyticFakeData/PAddPoissonNoise.cpp     | 12 +++----
 src/tests/analyticFakeData/PAddPoissonNoise.h |  4 +--
 src/tests/analyticFakeData/analyticFakeData.C | 36 ++++++++++++++++---
 3 files changed, 39 insertions(+), 13 deletions(-)

diff --git a/src/tests/analyticFakeData/PAddPoissonNoise.cpp b/src/tests/analyticFakeData/PAddPoissonNoise.cpp
index bc8f2cbe..1e3733e6 100644
--- a/src/tests/analyticFakeData/PAddPoissonNoise.cpp
+++ b/src/tests/analyticFakeData/PAddPoissonNoise.cpp
@@ -50,7 +50,7 @@ PAddPoissonNoise::PAddPoissonNoise(UInt_t seed)
 {
   fValid = true;
 
-  fRandom = new TRandom2(seed);
+  fRandom = new TRandom3(seed);
   if (fRandom == 0) {
     fValid = false;
   }
@@ -144,7 +144,7 @@ Double_t PAddPoissonNoise::PoiDev(const Double_t &mean)
   if (mean < 12.0) {
     if (mean != fOldMean) {
       fOldMean = mean;
-      fG = TMath::Exp(-mean);
+      fG = exp(-mean);
     }
     em = -1.0;
     t  =  1.0;
@@ -155,17 +155,17 @@ Double_t PAddPoissonNoise::PoiDev(const Double_t &mean)
   } else {
     if (mean != fOldMean) {
       fOldMean = mean;
-      fSquareRoot = TMath::Sqrt(2.0*mean);
-      fAlxm = TMath::Log(mean);
+      fSquareRoot = sqrt(2.0*mean);
+      fAlxm = log(mean);
       fG = mean*fAlxm-TMath::LnGamma(mean+1.0);
     }
     do {
       do {
-        y = TMath::Tan(TMath::Pi()*fRandom->Rndm());
+        y = tan(TMath::Pi()*fRandom->Rndm());
         em = fSquareRoot*y+mean;
       } while (em < 0.0);
       em = TMath::Floor(em);
-      t = 0.9*(1.0+y*y)*TMath::Exp(em*fAlxm-TMath::LnGamma(em+1.0)-fG);
+      t = 0.9*(1.0+y*y)*exp(em*fAlxm-TMath::LnGamma(em+1.0)-fG);
     } while (fRandom->Rndm() > t);
   }
 
diff --git a/src/tests/analyticFakeData/PAddPoissonNoise.h b/src/tests/analyticFakeData/PAddPoissonNoise.h
index 2e3ddca4..4a8d1515 100644
--- a/src/tests/analyticFakeData/PAddPoissonNoise.h
+++ b/src/tests/analyticFakeData/PAddPoissonNoise.h
@@ -34,7 +34,7 @@
 
 #include <TObject.h>
 #include <TH1F.h>
-#include <TRandom2.h>
+#include <TRandom3.h>
 
 class PAddPoissonNoise : public TObject
 {
@@ -49,7 +49,7 @@ class PAddPoissonNoise : public TObject
 
   private:
     Bool_t fValid;
-    TRandom2 *fRandom;
+    TRandom3 *fRandom;
 
     Double_t fSquareRoot;
     Double_t fAlxm;
diff --git a/src/tests/analyticFakeData/analyticFakeData.C b/src/tests/analyticFakeData/analyticFakeData.C
index ec8b29c1..8870486f 100644
--- a/src/tests/analyticFakeData/analyticFakeData.C
+++ b/src/tests/analyticFakeData/analyticFakeData.C
@@ -47,7 +47,7 @@ void analyticFakeData(const TString type)
   gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
   decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
 
-  UInt_t runNo = 9001;
+  UInt_t runNo = 9042;
   TString fileName, tstr, label, tstr1;
   fileName.Form("0%d.root", (Int_t)runNo);
 
@@ -164,13 +164,14 @@ void analyticFakeData(const TString type)
   }
 
   // create histos
-  UInt_t n0[8] = {200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0};
-  UInt_t bkg[8] = {1.3, 1.5, 1.0, 1.3, 1.2, 1.1, 1.0, 1.4};
+  Double_t n0[8] = {200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0, 200.0};
+  Double_t bkg[8] = {1.3, 1.5, 1.0, 1.3, 1.2, 1.1, 1.0, 1.4};
   const Double_t tau = 2197.019; // ns
+
   // asymmetry related stuff
   const Double_t timeResolution = 0.1953125; // ns
 //  Double_t a0[8] = {0.12, 0.12, 0.12, 0.12, 0.12, 0.12, 0.12, 0.12};
-  Double_t a0[8] = {0.2201, 0.22, 0.2202, 0.2198, 0.22, 0.2199, 0.22, 0.2203};
+  Double_t a0[8] = {0.21, 0.22, 0.22, 0.24, 0.22, 0.21, 0.22, 0.23};
 //  Double_t a1[8] = {0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05};
   Double_t a1[8] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
   Double_t phase[8] = {5.0*TMath::Pi()/180.0, 50.0*TMath::Pi()/180.0, 95.0*TMath::Pi()/180.0, 140.0*TMath::Pi()/180.0,
@@ -178,9 +179,34 @@ void analyticFakeData(const TString type)
   const Double_t gamma = 0.0000135538817; // gamma/(2pi)
   Double_t bb0 = 200.0; // field in Gauss
   Double_t bb1 = 400.0; // field in Gauss
-  Double_t rate0 = 0.5/1000.0; // in 1/ns
+  Double_t rate0 = 0.10/1000.0; // in 1/ns
   Double_t rate1 = 0.15/1000.0; // in 1/ns
 
+  // fake function parameters header info: only for test purposes
+  if (type.CompareTo("TLemRunHeader")) {
+    TDoubleVector dvec;
+    header->Set("FakeFct/Def", "N0 exp(-t/tau_mu) [1 + A exp(-t lambda) cos(gamma_mu B t + phi)] + bkg");
+    for (UInt_t i=0; i<8; i++)
+      dvec.push_back(n0[i]);
+    header->Set("FakeFct/N0", dvec);
+    dvec.clear();
+    for (UInt_t i=0; i<8; i++)
+      dvec.push_back(bkg[i]);
+    header->Set("FakeFct/bkg", dvec);
+    dvec.clear();
+    for (UInt_t i=0; i<8; i++)
+      dvec.push_back(a0[i]);
+    header->Set("FakeFct/A", dvec);
+    dvec.clear();
+    for (UInt_t i=0; i<8; i++)
+      dvec.push_back(phase[i]);
+    header->Set("FakeFct/phase", dvec);
+    prop.Set("B", bb0, "G");
+    header->Set("FakeFct/B", prop);
+    prop.Set("lambda", rate0, "1/usec");
+    header->Set("FakeFct/lambda", prop);
+  }
+
   TH1F *histo[16];
   char str[128];
   for (UInt_t i=0; i<8; i++) {