CubicFilter.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 <fstream>

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

#ifndef Gw_CubicFilter
#include "CubicFilter.h"
#endif

#ifndef Gw_Link
#include "Link.h"
#endif

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

using namespace Gw;

//
ClassImp(CubicFilter)


//__________________________________________________________
CubicFilter::CubicFilter() 
  : CSFilter("CubicFilter", "Cubic Gating"),
        fGateList(new TList()),
        IsPeak(true)
{ 
  // default constructor
  fLog.GetProcessName() = "CubicFilter";
  fGateList->SetOwner(true);
  for (Int_t i = 0; i < 4; ++i) {
        fGateList->Add(new Gate());
  }
}


//__________________________________________________________
CubicFilter::~CubicFilter()
{
  // default destructor
  delete fGateList;
}

//__________________________________________________________
void CubicFilter::Reset()
{
  CSFilter::Reset();
  ResetGates();
}

//__________________________________________________________
Bool_t CubicFilter::IsListFilled(Int_t idx)
{
   Gate* gate = static_cast<Gate*> (fGateList->At(idx));
   
   if (gate) 
          if (gate->GetNofGates() > 0)  return true;
   
        return  false;
}

//__________________________________________________________
void CubicFilter::ImportGates(const Char_t* name,  Option_t* opt)
{
  
  Gate* gate = static_cast<Gate*> (fGateList->At(0));

  IsPeak = false;
  std::ifstream filein(name);
  Double_t gateInf;
  Double_t gateSup;
  
  while ( filein.good() ) {
        filein >> gateInf >> gateSup;
        gate->AddGate(gateInf, gateSup, opt);
  }
  
  filein.close();
}

//__________________________________________________________
void CubicFilter::AddGate(Double_t gateInf, Double_t gateSup, Int_t axis, Option_t* opt)
{
   
   if (axis < fGateList->GetEntries()) {
          Gate* gate = static_cast<Gate*> (fGateList->At(axis));
          gate->AddGate(gateInf, gateSup, opt);
   } else {
          fLog.SetProcessMethod("AddGate(Double_t, Double_t, Int_t, Option_t* )");
          fLog << error << "axis not defined" << dolog;
   }
}

//__________________________________________________________
void CubicFilter::MakeGates(Option_t* opt)
{
  TString tmp(opt);
  tmp.ToLower();
  
  if (IsPeak) {
        ResetGates();
        if (tmp.Contains("peak"))
          FillGates();
        else 
          FillBkgGates();
        IsPeak = true;
  }
}

//__________________________________________________________
Bool_t CubicFilter::IsCombined()
{
   Int_t nList = 0;
   for (Int_t iGate = 0; iGate < fGateList->GetEntries(); ++iGate) {
          
          Gate* gate = static_cast<Gate*> ( fGateList->At(iGate) );
          if (gate->GetNofGates() > 0) nList++;
   }
   
   return (nList > 1);
}  

//__________________________________________________________
Int_t CubicFilter::IsInside(Double_t* energies, Int_t mult, Short_t* inGate, Int_t iGate)
{   
   Int_t nGate = 0;
   Gate* gate  = static_cast<Gate*> ( fGateList->At(iGate) );
   
   if (gate->GetNofGates() <= 0) {
          fLog.SetProcessMethod("IsInside(Double_t*, Int_t, Short_t*, Int_t )");
          fLog << error << "Gate index not defined" << dolog;
          return 0;
   }
   
   for(Int_t i = 0; i < mult; ++i) {
          inGate[i] = 0;
          Bool_t IsInGate = false;
          
          for(Int_t k = 0; k < gate->GetNofGates(); ++k) {
                 if( (energies[i] >= gate->GetGateInf(k)) && (energies[i] <= gate->GetGateSup(k)) ) {
                        IsInGate = true;
                 }
          }  
          
          if (IsInGate) {
                 inGate[i] = 1;
                 nGate++;
          }
   }
   
   if (nGate < fCondition) {
          for(Int_t i = 0; i < mult; ++i) 
                 inGate[i] = 0;
   }
   
   if (nGate > fCondition+fDimension) {
          for(Int_t i = 0; i < mult; ++i) 
                 inGate[i] = 1;
   }
   
   return nGate;
}

// private methods

//__________________________________________________________
void CubicFilter::ResetGates()
{
  
  for (Int_t i = 0; i < fGateList->GetEntries(); ++i) {
        Gate* gate = static_cast<Gate*> ( fGateList->At(i) );
        gate->Reset();
  }
}

//__________________________________________________________
void CubicFilter::FillGates()
{
  Gate* gate   = 0x0;
  Peak1D* peak = 0x0;
  IsPeak = true;

  TIter next(fPeakList);
  while ( (peak = (Peak1D*)next()) )  {
        Int_t idx = Peak1D::GetGateAxis(peak->GetLineColorPeak());
        if (idx == -1) continue;  // should not happen
        gate = static_cast<Gate*> ( fGateList->At(idx) );
        gate->AddGate(peak->GetPosition()-peak->GetFWHM()/2., peak->GetPosition()+peak->GetFWHM()/2.);
  }
  IsBkg = false;
}

//__________________________________________________________
void CubicFilter::FillBkgGates()
{
  Gate* gate   = 0x0;
  Peak1D* peak = 0x0;
  
  TIter next(fPeakList);
  while ( (peak = (Peak1D*)next()) )  {
        
        Int_t idx = Peak1D::GetGateAxis(peak->GetLineColorPeak());
        if (idx == -1) continue;
        
        gate = static_cast<Gate*> ( fGateList->At(idx) );
  
        TVector2 vecL = peak->GetBkgLeft();
        TVector2 vecR = peak->GetBkgRight();

        gate->AddGate(vecL.X(), vecL.Y());
        gate->AddGate(vecR.X(), vecR.Y());

  }
  IsBkg = true;
}


// Gate
ClassImp(Gate)

//__________________________________________________________
Gate::Gate() 
: TObject(),
  fGateInf(),
  fGateSup()
{ 
  // default constructor
}

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

//__________________________________________________________
void Gate::Reset()
{
  fGateInf.Reset(0.);
  fGateSup.Reset(0.);
  fGateInf.Set(0);
  fGateSup.Set(0);
}

//__________________________________________________________
void Gate::AddGate(Double_t gateInf, Double_t gateSup, Option_t* opt)
{
  TString option(opt);
  
  option.ToLower();
  if (option.Contains("reset")) {
        Reset();
  }
  
  fGateInf.Set(fGateInf.GetSize()+1);
  fGateInf.AddAt(gateInf, fGateInf.GetSize()-1);
  
  fGateSup.Set(fGateSup.GetSize()+1);
  fGateSup.AddAt(gateSup, fGateSup.GetSize()-1);
}