CalibCo60.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004-2006 by Olivier Stezowski & Christophe Theisen     *
 *   stezow(AT)ipnl.in2p3.fr, christophe.theisen(AT)cea.fr                       *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
 
#ifndef GW_CALIBCO60_H
#include "CalibCo60.h"
#endif

#ifndef ROOT_TROOT
#include "TROOT.h"
#endif
#ifndef ROOT_TPad
#include "TPad.h"
#endif
#ifndef ROOT_TSpectrum
#include "TSpectrum.h"
#endif
#ifndef ROOT_TMath
#include "TMath.h"
#endif
#ifndef ROOT_TCanvas
#include "TCanvas.h"
#endif

#ifndef Gw_Peak1D
#include "Peak1D.h"
#endif

#include <iostream>
#include <fstream>

using namespace std;
using namespace Gw;
 
// ROOT dictionnary
ClassImp(CalibCo60);

//__________________________________________________________________________________________________
CalibCo60::CalibCo60()
{

}

//__________________________________________________________________________________________________
CalibCo60::~CalibCo60()
{
        
}

//__________________________________________________________________________________________________
void CalibCo60::Calibrate(const TH1 */*histo*/, TF1 *function, Double_t channel1, Double_t channel2, Option_t */*AxisOption*/)
{  
  // The user can give the 2 channels corresponding to the 2 peaks
        channel1 = TMath::Abs(channel1);
        channel2 = TMath::Abs(channel2);
        if (channel1 == channel2) {
          printf("info in <Calibco60::DoCalib> parameters are identical !!\n");
          return;
        }
  
        if (channel1 > channel2) { //just invert the 2 given values
                Double_t tmp = channel2;
                channel2 = channel1;
                channel1 = tmp;
        }
  
        //channel1 corresponds to peak1 = 1173.228
        //channel2 corresponds to peak2 = 1332.492

        fRaw.Set(2);
        fRaw[0] = channel1; fRaw[1] = channel2;
        fTabulated.Set(2);
        fTabulated[0] = 1173.228 ; fTabulated[1] = 1332.492;

        BaseCalib::Calibrate(function);
}

//__________________________________________________________________________________________________
void CalibCo60::Calibrate(const TH1 *histo, TF1 *function, Option_t* /*AxisOption*/, Double_t sigma, 
                                                  Option_t* SearchOption, Double_t threshold)
{
        // The assumption is that the 2 most intense peaks are 1.17 and 1.33 MeV respectively
    // The method TSpectrum::Search is used in the range of the given histogram
        // The user have to change the range of the histogram to search peaks
        // The other parameters (sigma,option, threshold) are those of TSpectrum::Search (see the ROOT documenatation for details)
        

        // TSpectrum part
        TSpectrum spectrum;
        Int_t NbPeaks = spectrum.Search(histo,sigma,SearchOption,threshold);
        if (NbPeaks <= 1) {
          printf("info in <CalibCo60::DoCalib> ERROR !! Number of Peaks found < 2 !!\n"); 
          return;
        }
  
        //Double_t *tab_X = spectrum.GetPositionX();
        // Double_t *tab_Y = spectrum.GetPositionY();
        Int_t *index = new Int_t[NbPeaks];
        TMath::Sort(NbPeaks,spectrum.GetPositionY(),index,1);
        Peak1D p1, p2;
        p1.SetPosition(spectrum.GetPositionX()[index[0]]); p1.SetIntensity(spectrum.GetPositionY()[index[0]]);
        p2.SetPosition(spectrum.GetPositionX()[index[1]]); p2.SetIntensity(spectrum.GetPositionY()[index[1]]);

        if (gDebug > 0) {
                printf("Nb Peaks found = %i\n",NbPeaks);
                for (Int_t i=0 ; i<NbPeaks ; i++) {
                        printf("  [%i] Position : %8.3f       Intensity : %8.3f\n",i,spectrum.GetPositionX()[i],spectrum.GetPositionY()[i]);
                }
        }

        fRaw.Set(2);
        fTabulated.Set(2);
        fTabulated[0] = 1173.228 ; fTabulated[1] = 1332.492;
        
        // perform the calibration
        if (p1.GetPosition() == p2.GetPosition()) {
          printf("info in <Calibco60::DoCalib> The 2 most intense peaks have got the same positions !! Abort calibration !!\n"); 
          return;
        }
        if (p1.GetPosition() < p2.GetPosition()) {
                fRaw[0] = p1.GetPosition(); fRaw[1] = p2.GetPosition();
        } else {
                fRaw[0] = p2.GetPosition(); fRaw[1] = p1.GetPosition();
        }
        
        BaseCalib::Calibrate(function);

        //delete pointers
        if (index) delete [] index;
        
}