GSPlayer.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004-2006 by Olivier Stezowski & Christian Finck        *
 *   stezow(AT)ipnl.in2p3.fr, cfinck(AT)ires.in2p3.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.             *
 ***************************************************************************/

#include <Riostream.h>
#include <vector>

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

#ifndef Gw_GSPlayer
#include "GSPlayer.h"
#endif

#ifndef ROOT_TCanvas
#include "TCanvas.h"
#endif

#ifndef ROOT_TFormula
#include "TFormula.h"
#endif

#ifndef ROOT_TH1
#include "TH1.h"
#endif

#ifndef ROOT_TH2
#include "TH2.h"
#endif

#ifndef ROOT_KeySymbols
#include "KeySymbols.h"
#endif

#ifndef ROOT_TList
#include "TList.h"
#endif

#ifndef ROOT_TLine
#include "TLine.h"
#endif

#ifndef ROOT_TPolyLine
#include "TPolyLine.h"
#endif

#ifndef ROOT_TSpline
#include "TSpline.h"
#endif

#ifndef ROOT_TString
#include "TString.h"
#endif

#ifndef ROOT_TSystem
#include "TSystem.h"
#endif

#ifndef ROOT_TVector2
#include "TVector2.h"
#endif

#ifndef ROOT_TVirtualPad
#include "TVirtualPad.h"
#endif


using namespace Gw;

ClassImp(GSPlayer)

//__________________________________________________________
GSPlayer::GSPlayer() 
: BaseSpectrumPlayer()
{ 
  fLog.GetProcessName() = "GSPlayer";

  //Default Constructor
}

//__________________________________________________________
GSPlayer::~GSPlayer()
{
  // default destructor
}

//__________________________________________________________
Bool_t GSPlayer::SetParameter(const char* name, Double_t value)
{
  // set parameter
  Bool_t ok = BaseSpectrumPlayer::SetParameter(name,value);
  if ( !ok ) {
        // set search parameter
//      if (strncmp(name, "FWHM", strlen(name)) == 0)
//      fDefaultFwhm = TParameter<Double_t>(name, value);   
  }
  if ( !ok ) {
    fLog.SetProcessMethod("SetParameter(const char*, Double_t )");
    fLog << error << "Unkown parameter: " << name << dolog; 
  }
  return ok;
}

//__________________________________________________________
Int_t GSPlayer::FindPeaks(TH1 *histo, Option_t* o)
{
  // placed peak with a moving on line
  
  Int_t xBinMin  = histo->GetXaxis()->GetFirst();
  Int_t xBinMax  = histo->GetXaxis()->GetLast();
  Double_t xMin  = histo->GetXaxis()->GetBinCenter(xBinMin);
  Double_t xMax  = histo->GetXaxis()->GetBinCenter(xBinMax);
  Double_t xMin0 = xMin;
  Double_t xMax0 = xMax;
  Double_t yMax  = histo->GetMaximum()*1.05;
  Double_t yMin  = histo->GetMinimum()*0.95; 
  
  Double_t range  = xMax - xMin; 
  Double_t width  = GetDefaultPeakWidth()->Eval(xMin)*15;   
  Int_t nSpec     = int(range/width + 0.5);
  Int_t xBin      = 0;
  
  Int_t result    = 0;
  Int_t eventType = 0;; 
  
  // if Contains +, add to the list 
  TString opt = o;
  if ( !opt.Contains("+") ) 
    GetPeakList()->Delete();
  
   Option_t* optPeak;
   if ( opt.Contains("bg") ) 
          optPeak = "bg";
   else 
          optPeak = "";
   

           
  // activated search for 1D spectrum
  if ( histo->GetDimension() == 1 ) {
    // starts search once button1 is pressed
    do {
      gSystem->ProcessEvents();
      eventType = TVirtualPad::Pad()->GetEvent();
    } while(eventType != kButton1Down);
        
    // loop over ranges
    for (Int_t k = 0; k < nSpec; ++k) {
      xMax = xMin+width;
      histo->SetAxisRange(xMin, xMax);
      xMin = xMax;
      TLine line;
      xBinMin = histo->GetXaxis()->GetFirst();
      xBinMax = histo->GetXaxis()->GetLast();
      yMax    = histo->GetMaximum();
      yMin    = histo->GetMinimum(); 

      // loop over bins of the given range
      for (Int_t i = xBinMin; i < xBinMax; ++i) {
        
        xBin = int(histo->GetXaxis()->GetBinCenter(i) + 0.5);
         
        // set line
        line.SetX1(xBin); line.SetY1(yMin);
        line.SetX2(xBin); line.SetY2(yMax);
        line.Draw();
        TVirtualPad::Pad()->Update();
        
        // keep position of line.
        gSystem->ProcessEvents();
        // press q to proceed next range
        if (TVirtualPad::Pad()->GetEvent() == kKeyPress) {
          Int_t keysym = TVirtualPad::Pad()->GetEventY();
                   if ( (EKeySym)keysym == kKey_s) { 
                          fCreator->CreatePeak(histo->FindBin(xBin), optPeak);
                          TVirtualPad::Pad()->Update();
                   }
        }        
        
                 Int_t kk;
                 do {
                        gSystem->ProcessEvents();
                        eventType = TVirtualPad::Pad()->GetEvent();
                        kk = TVirtualPad::Pad()->GetEventY();
                 } while((EKeySym)kk == kKey_s);
                 
        gSystem->Sleep(100);
      }
    }
    histo->SetAxisRange(xMin0, xMax0);
    TVirtualPad::Pad()->Update();
  }
  
  if ( histo->GetDimension() == 2) {
    fLog.SetProcessMethod("FindPeaks(TH1*, Option_t* )");
    fLog << error << "2D search peak not implemented yet" << dolog; 
    return 0;
  }
  
  CollectPeaks();
  result = GetPeakList()->GetEntries();
  
  return result;
}


//__________________________________________________________
TH1* GSPlayer::Background(const TH1* histo, Option_t* /*opt*/)
{
  // generates bkg for polymarker (for the moment)
  TH1* bkg        = 0x0;
  TPolyLine* line = 0x0;
  
  TList* list = TVirtualPad::Pad()->GetListOfPrimitives();
  for (Int_t i = 0; i < list->GetEntries(); ++i) {
    TObject* o       = list->At(i);
    const char* name = o->ClassName();
    TString tmp(name);
    
    // check polyline in pad
    if ( tmp.CompareTo("TPolyLine") == 0 )
      line = static_cast<TPolyLine*>(o);

    // check bkg spectrum in pad
    if ( o->InheritsFrom("TH1")) {
      bkg = static_cast<TH1*>(o);
      tmp = bkg->GetName();
      if (tmp.Contains("Bkg"))
        return bkg;
      else
        bkg = 0x0;
    }
  }
  
  if (!line && !bkg) {
    fLog.SetProcessMethod("Background(const TH1*)");
    fLog << error << "No background line/spectrum defined" << dolog; 
    return bkg;
  }
  
  Double_t* x = line->GetX();
  Double_t* y = line->GetY();
  Int_t     n = line->Size();
  
  bkg = (TH1*)histo->Clone();
  Int_t xBinMin  = bkg->FindBin(x[0]);
  Int_t xBinMax  = bkg->FindBin(x[n-1]);
  bkg->Reset("ICE");
  
  TSpline3 spline("Bkg", x, y, n);

  // building backgroud spectrum
  for (Int_t i = xBinMin; i < xBinMax; ++i) {
    Double_t xBin = histo->GetXaxis()->GetBinCenter(i);
    bkg->SetBinContent(i, spline.Eval(xBin));
  }
  
  return bkg;
}