added documentation to the .C files

Signed-off-by: acolijn <auke.pieter.colijn@gmail.com>

Author: acolijn

calibration/analyzer.C

@@ -1,4 +1,21 @@
 #define analyzer_cxx
+/*---------------------------------------------------------------------------------------------------*/
+//
+// analyzer.C Routine to analyze spectra and calculate the rate of the sources
+//
+// Usage:
+//  prompt> #include "analyzer.C"
+//  prompt> analyzer ana(<directory_with_energy_calibrated_rootfiles>,<analysis_output_rootfile>)
+//  prompt> ana.Loop()
+//
+// To inspect the fit results it is possible to plot the fit results. In analyzer.h  set:
+// #define PLOT_ON_SCREEN 1
+// 
+// To set the time interval in which a spectrum is calculated, in analyzer.h set:
+// #define TIME_INTERVAL <interval_in_seconds>
+//
+// A.P.
+/*---------------------------------------------------------------------------------------------------*/
 #include "analyzer.h"
 // RooFit include files
 #include <RooRealVar.h>
@@ -59,6 +76,10 @@ const float base_max_val = 2000;
 
 
 /*----------------------------------------------------------------------------------------------------*/
+
+//
+// Gaussian function + 2nd order polynomial for simple rate fitting
+//
 Double_t fitf(Double_t *v, Double_t *par)
 {
     Double_t arg = 0;
@@ -74,7 +95,6 @@ void analyzer::fit_spectrum(int ichannel, double *fit_range){
     //
     // RooFit based spectrum fitter
     //
-    //
     cout <<"analyzer::fit_spectrum  channel = "<<ichannel<<endl;
     //
     // identify the peaks in our specified energy range and store their peak ID.
@@ -106,11 +126,14 @@ void analyzer::fit_spectrum(int ichannel, double *fit_range){
     //
     string mc_file="";
     if       (ichannel == 2 || ichannel == 3){
-        mc_file = "/user/z37/Modulation/analysis/calibration/MC_ti44_modulation.root";
+//        mc_file = "/user/z37/Modulation/analysis/calibration/MC_ti44_modulation.root";
+        mc_file = "MC_ti44_modulation.root";
     } else if(ichannel == 4 || ichannel == 5){
-        mc_file = "/user/z37/Modulation/analysis/calibration/MC_co60_modulation.root";
+//        mc_file = "/user/z37/Modulation/analysis/calibration/MC_co60_modulation.root";
+        mc_file = "MC_co60_modulation.root";
     } else if(ichannel == 6 || ichannel == 7){
-        mc_file = "/user/z37/Modulation/analysis/calibration/MC_cs137_modulation.root";
+//        mc_file = "/user/z37/Modulation/analysis/calibration/MC_cs137_modulation.root";
+        mc_file = "MC_cs137_modulation.root";
     }
 
     TFile *f_mc = new TFile(mc_file.c_str(),"READONLY");
@@ -267,7 +290,7 @@ void analyzer::fit_spectrum(int ichannel){
     //
     // spectrum is a function of the energy
     //
-    RooRealVar E("E (keV)","E (keV)",emin,emax);
+    RooRealVar E("E","E (keV)",emin,emax);
 
     //
     // the background template for each of the sources obtained from a GEANT4 simulation

calibration/ecal.C

@@ -11,15 +11,29 @@
 #include <vector>
 #include <stdio.h>
 
+/*---------------------------------------------------------------------------------------------------*/
+//
+// ecal.C Routine to do the energy calibration of a modulation run
+//       
+// Usage:
+//  prompt> #include "ecal.C"
+//  prompt> ecal e(<directory_with_uncalibrated_rootfiles>,<calibration_output_rootfile>)
+//  prompt> e.Loop()
+//
+// Important parameter in ecal.h: CALIBRATION_MODE 0 : one calibration for full run (not adviced)
+//                                CALIBRATION_MODE 1 : re-calibrate every TIME_INTERVAL seconds (see ecal.h)
+//
+// A.P.
 /*---------------------------------------------------------------------------------------------------*/
 // (1) You need to identify where the first peak is in the histograms of the integrals.
 // (2) Set range_low and range_high to find the range in which the peak with cal_energy should be
 const double range_low[NUMBER_OF_CHANNELS] ={0.16e-6,0.14e-6,0.05e-6,0.05e-6,0.,0.,0.,0.};
 const double range_high[NUMBER_OF_CHANNELS]={1e-6,1e-6,1e-6,1e-6,1e-6,1e-6,1e-6,1e-6};
 
+/*---------------------------------------------------------------------------------------------------*/
 float source_energy[NUMBER_OF_CHANNELS][MAX_PEAKS] =
 //
-// the energy peaks you wish to select should be in this list
+// the energy peaks you wish to select for the calibration should be in this list
 // NOTE: the first peak should be the highest in the spectrum (sub-optimal, but handy for finding)
 //
 {

calibration/ecal.h

@@ -11,22 +11,17 @@
 
 /*----------------------------------------------------------------------------*/
 
-#define CHANNEL0 0
-#define CHANNEL1 1
-#define CHANNEL2 2
-#define CHANNEL3 3
-#define CHANNEL4 4
-#define CHANNEL5 5
-#define CHANNEL6 6
-#define CHANNEL7 7
-
 #define NUMBER_OF_CHANNELS 8
 #define MAX_PEAKS 5
 #define MAX_PARAMETERS 3
 
 #define BEFORE_CALIBRATION 0
 #define AFTER_CALIBRATION 1
 
+//
+// CALIBRATION_MODE 0 = one calibration for full run
+// CALIBRATION_MODE 1 = one calibration every TIME_INTERVAL seconds
+//
 #define CALIBRATION_MODE 1
 #define TIME_INTERVAL 1800 
 

monitor/fitspectrum.C

@@ -31,6 +31,24 @@ float source_energy[NUMBER_OF_CHANNELS][MAX_PEAKS] =
     {662.,-1,-1,-1,-1}  // channel7: 137Cs
 };
 
+/*---------------------------------------------------------------------------------------------*/
+//
+// fitspectrum.C - RooFit based spectrum fitter. This routine can fit the spectrum for a full
+//                 modulation experiment run. The same algorithm is used in the analysis/calibration/analyzer.C
+//                 routine to fit the spectra after the chosen time interval.
+//
+//                 In this routine it is easy to develop new fitting ideas.
+//
+// Arguments
+//            data_file    : name of the data file to fit
+//            mc_file      : Geant4 prediction of the background spectrum. The spectra are generated
+//                           with the modexp/modusim package from GitHub
+//            ichannel     : channel number [0..7]
+//            e0           : lowest energy for fit (keV)
+//            e1           : highest energy for fit (keV)
+//
+// A.P. 2015-07-02
+//
 /*---------------------------------------------------------------------------------------------*/
 void fitspectrum(string data_file, string mc_file, int ichannel, double e0, double e1){
     double fit_range[2];
@@ -113,11 +131,7 @@ void fitspectrum(string data_file, string mc_file, int ichannel, double e0, doub
     for (int id=0; id<nselect; id++){
         sprintf(vname,"gaus%i",id);
         cout <<"vname >>"<<vname<<"<< pk  = "<< (pk_mean[id])->getValV()<<endl;
-        
-//        RooGaussian *gg = new RooGaussian(vname,vname,E,*pk_mean[id],*pk_sigma[id]);
         pk_gaus.push_back(new RooGaussian(vname,vname,E,*pk_mean.at(id),*pk_sigma.at(id)));
-        cout <<"dodo"<<endl;
-//        delete gg;
     }
     cout <<"analyzer::fit_spectrum Define photo peak Gaussians ---- DONE"<<endl;
     //
@@ -128,17 +142,10 @@ void fitspectrum(string data_file, string mc_file, int ichannel, double e0, doub
     //
     // compose the joined pdf for background + signal
     //
-    
     cout <<"analyzer::fit_spectrum Compose the combined pdf"<<endl;
 
     RooAddPdf *sum;
     if (nselect==1)  sum = new RooAddPdf("sum","g1+bg",RooArgList(*pk_gaus[0],bg),RooArgList(*pk_frac[0]));
-    //      if(peak_id[0] == 2) {
-    //         sum = new RooAddPdf("sum","g1+g2+bg",RooArgList(*pk_gaus[0],*pk_gaus_tail[0],bg),RooArgList(*pk_frac[0],*pk_frac_tail[0]));
-    //      } else {
-    //         sum = new RooAddPdf("sum","g1+bg",RooArgList(*pk_gaus[0],bg),RooArgList(*pk_frac[0]));
-    //      }
-    //    }
     if (nselect==2)  sum = new RooAddPdf("sum","g1+g2+bg",RooArgList(*pk_gaus[0],*pk_gaus[1],bg),RooArgList(*pk_frac[0],*pk_frac[1]));
     if (nselect==3)  sum = new RooAddPdf("sum","g1+g2+g3+bg",RooArgList(*pk_gaus[0],*pk_gaus[1],*pk_gaus[2],bg),RooArgList(*pk_frac[0],*pk_frac[1],*pk_frac[2]));
 

monitor/lifetime.C

@@ -12,6 +12,25 @@
 
 const int nmax = 100000;
 
+/*----------------------------------------------------------------------------------------------------*/
+//
+// Fit the half-life to the tree data produced by the analyzer.C code. 
+//
+// Usage from the ROOT6 command line
+//
+// prompt> #include<lifetime.C>
+// prompt> lifetime l(<directory_with_ANA_files>)
+// prompt> l.Life(<channel>, <peak>, <type>, <save>)
+//
+// <channel> : [0..7]
+// <peak>    : [0..2] peak identifier in the spectrum. at the moment [0] for BG, [0..2] for 44Ti & 60Co, [0] for 137Cs
+// <type>    : what to plot : "life" = rate vs time wit exponential fit
+//                            "pull" = (data-fit)/fit_error
+// <save>    : to file or not?
+//
+// A.P.
+//
+/*----------------------------------------------------------------------------------------------------*/
 void lifetime::Life(int channel_sel, int peak_sel, string type, bool save)
 {
    TCanvas *c1 = new TCanvas("c1","c1",600,400);
@@ -59,6 +78,9 @@ void lifetime::Life(int channel_sel, int peak_sel, string type, bool save)
    char cmd[128];
    gStyle->SetOptFit(111);
    if(type == "life"){
+     //
+     // plot rate vs time and fitted exponential
+     //
      g1->SetMarkerStyle(24);
      g1->Draw("AP");
 
@@ -68,6 +90,9 @@ void lifetime::Life(int channel_sel, int peak_sel, string type, bool save)
      g1->GetYaxis()->SetTitle("rate (Hz)");
 
    } else if (type == "pull"){
+     //
+     // plot pull distribution
+     //
      double res;
      for(int i=0; i<n; i++){
        res = (R[i] - f1->Eval(t[i]) )/dR[i];

monitor/lifetime.h

@@ -61,7 +61,7 @@ lifetime::lifetime(string fname) : fChain(0)
 {
     char cmd[256];
     TChain *tree = new TChain("ana");
-    sprintf(cmd,"%s*.root",fname.c_str());
+    sprintf(cmd,"%s/*.root",fname.c_str());
     tree->Add(cmd);
    Init(tree);
 }