Differences

This shows you the differences between two versions of the page.

--- software:analysis [2020/04/21 06:17] – himuser
+++ software:analysis [2022/01/14 18:44] (current) – external edit 127.0.0.1
@@ Line 1: / Line 1: @@
-====== Analysis ======
+===== Analysis tools =====
-There are already several files which help you to analyse root files from primasoft:
+===Novosibirsk function===
-  - PrimaSimAnalysis (can be found here: /geant4/src/PRIMASimAnalysis.cc, /geant4/include/PRIMASimAnalysis.hh)
+The Novosibirsk function is typically used to fit monoenergetic electromagnetic calorimeter spectra.
-  - Config file for Analysis (examples can be found here: /geant4/sim_ana_utilities/ana_default.conf)
+Explanations on this function are given here {{ :software:novosibirsk.pdf |}}.
-  - Rootmacros (examples can be found here: /geant4/rootmacros)
+A C++ code snippet to calculate the function using ROOT classes is given below. The interface is chosen to be usable with the ROOT fitting routines, e.g. ''TH1::Fit()''. It takes two arrays as input. From the first one only the first element is used as ''x'', From the second array, the values of the fit parameters are taken.
+  #include "TMath.h"
-This is what they do:
+  /*   p[0] = normalisation constant
-== PrimaSimAnalysis ==
+       p[1] = peak position
-  * Sorts hits into clusters
+       p[2] = sigma (energy  resolution)
-  * Single crystal threshold is applied. The value van be set in the config file, default is 2 MeV. This is analogue the to hit detection threshold in a measurement.
+       p[3] = eta (asymmetry parameter)
-  * Different crystal setups can be used:
+       p[4] = constant background
-  *   * all crystals are used (768)
+  */
-  *   * the inner 640 crystals are used symmetrically (640)
+  double novosibirsk(double *x, double *p)
-  *   * the inner 640 crystals are used in an oval shape (640_oval)
+  {
-  *   * the inner 560 crystals are used symmetrically (560)
+     double xi=2.35482004503094933e+00;
+     double Y=(x[0]-p[1])*p[3]/p[2];
+     if(Y>=1) return p[4];
+     double lnX = log( 1. - Y );
+     double s0  = 2/xi*TMath::ASinH(p[3]*xi*0.5);
+     double s02 = s0*s0;
+     double f = p[0]*exp ( -0.5/s02 * lnX*lnX - 0.5*s02) + p[4];
+     return f;
+  }