MatPlayer.cpp

Go to the documentation of this file.

//
//
// MatPlayer Version 0.1
// MatPlayer is used to analyse a 2D histogram by putting gates interactively.
//
// A MatPlayer object must be connected to a TCanvas and a TH2D histogram. It is done by using SetCanvas & SetMatrix.
// Once a canvas and a matrix has been linked, any key pressed on the defined canvas gives a gate (TBox) that can be
// modify graphically. The gates on the current canvas are use to  
//
//      Olivier Stezowski - 
//
// 
//

//#ifndef __CINT__

/*
#ifndef ROOT_KeySymbols
#include "RTypes.h"
#endif

#ifndef ROOT_TROOT
#include "TROOT.h"
#endif

#ifndef ROOT_TROOT
#include <TFrame.h>
#endif

#ifndef ROOT_TApplication
#include "TApplication.h"
#endif

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

*/

#include "KeySymbols.h"

#ifndef ROOT_TStyle
#include "TStyle.h"
#endif

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

#ifndef ROOT_TSpectrum
#include "TSpectrum.h"
#endif

#ifndef ROOT_TFrame
#include "TFrame.h"
#endif

#ifndef ROOT_TROOT
#include "TROOT.h"
#endif

#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,0)
#include "TMath.h"
#endif

#ifndef GW_MATPLAYER_H
#include "MatPlayer.h"
#endif

using namespace std;
using namespace Gw;

ClassImp(MatPlayer);

Int_t MatPlayer::fLastX = 0, MatPlayer::fLastY = 0; TList *MatPlayer::fgPlayers = NULL;

MatPlayer::MatPlayer() : TNamed("Matplayer","To set gates interactively in a Canvas")
{
// constructor   
        fMat = NULL; fPx = fPy = fTotx = fToty = fBx = fBy = NULL;
   
        fLastGates.SetOwner(kTRUE); 
        // fGatesY.SetOwner(kTRUE); fGatesBX.SetOwner(kTRUE); fGatesBY.SetOwner(kTRUE);
        
        gStyle->SetOptStat(0);
  
        fStep.first = fStep.second = 1.0; 
        
        fWidthGate.first  = new TFormula("MatPlayer_WX","pol0"); fWidthGate.first->SetParameter(0,2.5);
        fWidthGate.second = new TFormula("MatPlayer_WY","pol0"); fWidthGate.second->SetParameter(0,2.5); 
  
        fBoxX.SetFillColor(2);  fBoxX.SetFillStyle(3002);
        fBoxY.SetFillColor(4);  fBoxY.SetFillStyle(3013);
        fBoxBX.SetFillColor(16); fBoxBX.SetFillStyle(3016); 
        fBoxBY.SetFillColor(16); fBoxBY.SetFillStyle(3020);
        
        if ( fgPlayers == NULL ) { fgPlayers = new TList(); fgPlayers->SetOwner(kTRUE); }
        fgPlayers->Add(this); 
}

MatPlayer::~MatPlayer()
{
// destructor
        Clear();
        if ( fWidthGate.first ) delete fWidthGate.first; if (fWidthGate.second ) delete fWidthGate.second;
        
        fgPlayers->Remove(this);
}   

void MatPlayer::Clear(const Option_t */*opt*/)
{
        if ( fTotx ) delete fTotx; fTotx = NULL; if ( fToty ) delete fToty; fToty = NULL;
        
        if ( fPx ) delete fPx; fPx = NULL; if ( fPy ) delete fPy; fPy = NULL;
        if ( fBx ) delete fBx; fBx = NULL; if ( fBy ) delete fBy; fBy = NULL;
}

void MatPlayer::SetCanvas()
{
// to connect the current canvas to the interactive gating 
        TVirtualPad *pad = TVirtualPad::Pad(); TCanvas *c = NULL;
        
        if ( pad ) {
                c = pad->GetCanvas();
                c->Clear(); c->Divide(1,1,0.001,0.001); c->cd(1);
        }
        else {
                if ( (c = dynamic_cast<TCanvas *>(gROOT->GetListOfCanvases()->FindObject("DefaultMatPlayer"))) == NULL ) {
                                c = new TCanvas("DefaultMatPlayer","Canvas to interact with MatPlayer",1);
                                c->Divide(1,1,0.001,0.001); c->cd(1);
                }
        }
                                
        c->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)",
                   "Gw::MatPlayer",this,"XEventAction(Int_t,Int_t,Int_t,TObject*)");    
        
        cout << " - Info - Interactive gating (by pressing any key) is now possible on any pad of the canvas " << c->GetName() << endl;
}

void MatPlayer::SetMatrix(const char *matname)
{
        TObject *obj = gDirectory->Get(matname);
        if ( obj == NULL ) {
                cout << "Cannot find matrix named " << matname << " in the current TDirectory " << endl; 
        } else {        
                if ( obj->InheritsFrom("TH2") ) {
                        SetMatrix(dynamic_cast<TH2 *>(obj)); SetName(obj->GetName());
                }
                else 
                        cout << "Object named " << matname << " is not a matrix !! " << endl;
        }
}

void MatPlayer::SetMatrix(TH2 *mat)
{
        TString hname; TH1D *Px, *Py, *Bx, *By, *Totx, *Toty;
                        
// nothing checked to avoid wrong pointers for mat .. except that this function is private
        if ( fMat == mat ) return;
        
        fMat = mat; Totx = Toty = Px = Py = Bx = By = NULL;
        
// it creates necessary histograms (or checks if they already exist) from mat on the same location (TDirectory) than the matrix
        if ( mat->GetDirectory() ) { // mat belongs to a TDirectory .. so check if the histograms have already been created
                hname = mat->GetName(); hname += "_tx" ;
                Totx  = dynamic_cast<TH1D *>( mat->GetDirectory()->Get(hname.Data()) );
                hname = mat->GetName(); hname += "_ty" ;
                Toty  = dynamic_cast<TH1D *>( mat->GetDirectory()->Get(hname.Data()) );
                hname = mat->GetName(); hname += "_px" ;
                Px    = dynamic_cast<TH1D *>( mat->GetDirectory()->Get(hname.Data()) );
                hname = mat->GetName(); hname += "_py" ;
                Py    = dynamic_cast<TH1D *>( mat->GetDirectory()->Get(hname.Data()) );
                hname = mat->GetName(); hname += "_bx" ;                
                Bx    = dynamic_cast<TH1D *>( mat->GetDirectory()->Get(hname.Data()) );
                hname = mat->GetName(); hname += "_by" ;
                By    = dynamic_cast<TH1D *>( mat->GetDirectory()->Get(hname.Data()) );         
        }
        // histograms to store full projection
        if ( Totx == NULL ) { hname = mat->GetName(); hname += "_tx"; Totx = mat->ProjectionX(hname.Data()); }
        if ( Toty == NULL ) { hname = mat->GetName(); hname += "_ty"; Toty = mat->ProjectionY(hname.Data()); }
                
        // histograms to store gated spectra
        if ( Px == NULL ) { hname = mat->GetName(); hname += "_px"; Px = (TH1D *)Totx->Clone(hname.Data()); Px->Reset("ICE"); }
        if ( Py == NULL ) { hname = mat->GetName(); hname += "_py"; Py = (TH1D *)Toty->Clone(hname.Data()); Py->Reset("ICE"); }
        
        // histograms to store background spectra
        if ( Bx == NULL ) hname = mat->GetName(); hname += "_bx"; Bx = (TH1D *)Totx->Clone(hname.Data()); Bx->Reset("ICE"); 
        if ( By == NULL ) hname = mat->GetName(); hname += "_by"; By = (TH1D *)Toty->Clone(hname.Data()); By->Reset("ICE");     
        
// now keep a local copy protected from delete
        // remove previous projection
        Clear();
        
        // clone local histograms and set directory to NULL     
        fPx = dynamic_cast<TH1D *>( Px->Clone("MatPlayer_px") ); fPx->SetDirectory(0); 
        fPy = dynamic_cast<TH1D *>( Py->Clone("MatPlayer_py") ); fPy->SetDirectory(0);
        fBx = dynamic_cast<TH1D *>( Bx->Clone("MatPlayer_bx") ); fBx->SetDirectory(0); 
        fBy = dynamic_cast<TH1D *>( By->Clone("MatPlayer_by") ); fBy->SetDirectory(0);
        
        fTotx = dynamic_cast<TH1D *>( Totx->Clone("MatPlayer_tx") ); fTotx->SetDirectory(0); 
        fToty = dynamic_cast<TH1D *>( Toty->Clone("MatPlayer_ty") ); fToty->SetDirectory(0);
        
        Totx->Draw();
}

Double_t MatPlayer::GetWidthGate(Double_t pos, Option_t *o)
{
        TString opt = o;
        if ( opt.Contains("Gx") ) return fWidthGate.first->Eval(pos);  
        else return fWidthGate.second->Eval(pos); 
}

void MatPlayer::SetGateX(Option_t *o)
{
        Double_t posx, posy;    
        bool isTH2; TString opt = o; Double_t w, x1, x2, y1, y2; TH1 *hist; TObject *obj; TFrame *frame; 
        
        // check for a valid TPad and the first histogram
        TVirtualPad *pad = TVirtualPad::Pad(); if ( pad == NULL ) return;
        posx = pad->PadtoX(pad->AbsPixeltoX(fLastX));
        posy = pad->PadtoY(pad->AbsPixeltoY(fLastY));

        TIter next(pad->GetListOfPrimitives()); // look for an histogram
        frame = pad->GetFrame(); 
        isTH2 = false;
        hist = NULL; 
        while ( (obj = next()) ) { 
                if ( obj->InheritsFrom("TH1") ) { 
                                hist = dynamic_cast<TH1 *>(obj); if ( obj->InheritsFrom("TH2") ) isTH2 = true;  break; }
        }
        if ( hist == NULL || frame == NULL ) { cout << "Cannot add a gate, no spectrum selected " << endl; return; } 
        
        // default is a horizontal gate on X
        TBox *box = &fBoxX; 
        w = GetWidthGate(posx,"Gx");
        x1 = posx-w/2.0; y1 = frame->GetY1(); x2 = posx+w/2.0; y2 = frame->GetY2();
        
        // gate on Y depends if the histogram is a 1D or 2D histogram
        if ( opt.Contains("Gy") ) { 
                box = &fBoxY;
                if ( isTH2 ) { // horizontal gate
                        w = GetWidthGate(posy,"Gy"); 
                        x1 = frame->GetX1(); y1 = posy-w/2.0; x2 = frame->GetX2(); y2 =posy+w/2.0;
                }
                else { // vertical gate
                        w = GetWidthGate(posx,"Gy");            
                        x1 = posx-w/2.0; y1 = frame->GetY1(); x2 = posx+w/2.0; y2 = frame->GetY2();
                }
        }
        if ( opt.Contains("Bx") ) { 
                box = &fBoxBX; w = GetWidthGate(posx,"Gx"); 
                x1 = posx-w/2.0; y1 = frame->GetY1(); x2 = posx+w/2.0; y2 = frame->GetY2();
        }
        if ( opt.Contains("By") ) { 
                box = &fBoxBY; 
                if ( isTH2 ) { // horizontal gate
                        w = GetWidthGate(posy,"Gy"); 
                        x1 = frame->GetX1(); y1 = posy-w/2.0; x2 = frame->GetX2(); y2 =posy+w/2.0;
                }
                else { // vertical gate
                        w = GetWidthGate(posx,"Gy");
                        x1 = posx-w/2.0; y1 = frame->GetY1(); x2 = posx+w/2.0; y2 = frame->GetY2();
                }               
        }
        
        // Add one gates
        box->DrawBox(x1,y1,x2,y2); pad->Update();
}

bool MatPlayer::IsBox(const TBox *box, const TBox *boxref) const
{
        return ( box->GetFillStyle() == boxref->GetFillStyle()  && box->GetFillColor() == boxref->GetFillColor() );
}

Int_t MatPlayer::InitGating(Option_t *o)
{
// Init gating - It looks in the current pad for gates (boxes) with the good characteristics (color and fill style).
// It returns the number of gates found
        
        TObject *obj; TString opt = o;
        
        if ( TVirtualPad::Pad() == NULL ) return 0;
        
        TIter next(TVirtualPad::Pad()->GetListOfPrimitives()); // look for an histogram
        bool isTH2 = false, overlap;
        while ( (obj = next()) ) { 
                if ( obj->InheritsFrom("TH1") ) { 
                         if ( obj->InheritsFrom("TH2") ) isTH2 = true;  break; }
        } // while
        
// look for gates and cuts and prepare them for gating
        TBox *boxnew, *boxadd ; fLastGates.Delete(); // remove previous gates

        // re-start iteration and look for the right gates
        next.Reset();
        while ( (boxnew = (TBox *)next()) ) { // look for a box and keep a copy that can be i=used for gating
                if ( boxnew->InheritsFrom("TBox") && opt.Contains("x") ) {
                        // Gate on X
                        if ( IsBox(boxnew,&fBoxX) ) {
                                // look for an existing gate that overlaps this one
                                TIter existing(&fLastGates);
                                overlap = false;
                                while ( (boxadd = (TBox *)existing()) ) {
                                        if ( (boxnew->GetX1() >= boxadd->GetX1() && boxnew->GetX1() <= boxadd->GetX2()) ||
                                                 (boxnew->GetX2() >= boxadd->GetX1() && boxnew->GetX2() <= boxadd->GetX2()) ) { overlap = true; break; }
                                }
                                if ( overlap == false ) { fLastGates.Add(boxnew->Clone()); } // new gate or
                                else { // change current box settings
                                        boxadd->SetX1( TMath::Min(boxadd->GetX1(),boxnew->GetX1()) );
                                        boxadd->SetX2( TMath::Max(boxadd->GetX2(),boxnew->GetX2()) );
                                }
                        }
                        // Background gate on X
                        if ( IsBox(boxnew,&fBoxBX) ) {
                                // look for an existing gate that overlaps this one
                                TIter existing(&fLastGates);
                                overlap = false;                                
                                while ( (boxadd = (TBox *)existing()) ) {
                                        if ( (boxnew->GetX1() >= boxadd->GetX1() && boxnew->GetX1() <= boxadd->GetX2()) ||
                                                 (boxnew->GetX2() >= boxadd->GetX1() && boxnew->GetX2() <= boxadd->GetX2()) )  { overlap = kTRUE; break; }
                                }
                                if ( overlap == false ) { fLastGates.Add(boxnew->Clone()); } // new gate or
                                else { // change current box settings
                                        boxadd->SetX1( TMath::Min(boxadd->GetX1(),boxnew->GetX1()) );
                                        boxadd->SetX2( TMath::Max(boxadd->GetX2(),boxnew->GetX2()) );
                                }
                         }                      
                } // x
                if ( boxnew->InheritsFrom("TBox") && opt.Contains("y") ) {
                        // Gate on Y ... Y1 and Y2 are used for real gating ... so in case gates are set on a 1D histogram
                        // one should use X1 and X2
                        Double_t realY1, realY2;
                        if ( isTH2 ) { // 
                                realY1 = boxnew->GetY1();
                                realY2 = boxnew->GetY2();
                        }
                        else {
                                realY1 = boxnew->GetX1();
                                realY2 = boxnew->GetX2();       
                        }
                        if ( IsBox(boxnew,&fBoxY) ) {
                                // look for an existing gate that overlaps this one
                                TIter existing(&fLastGates);
                                overlap = false;                                
                                while ( (boxadd = (TBox *)existing()) ) {
                                        if ( (realY1 >= boxadd->GetY1() && realY1 <= boxadd->GetY2()) ||
                                             (realY2 >= boxadd->GetY1() && realY2 <= boxadd->GetY2()) ) { overlap = true; break; }
                                }
                                if ( overlap == false ) { 
                                        TBox *box = dynamic_cast<TBox *>(boxnew->Clone());
                                        box->SetY1(realY1);
                                        box->SetY2(realY2);     
                                        fLastGates.Add(box);
                                } // new gate 
                                else { // change current box settings
                                        boxadd->SetY1( TMath::Min(boxadd->GetY1(),realY1) );
                                        boxadd->SetY2( TMath::Max(boxadd->GetY2(),realY2) );
                                }
                        }
                        // Background gate on Y
                        if ( IsBox(boxnew,&fBoxBY) ) {
                                // look for an existing gate that overlaps this one
                                TIter existing(&fLastGates);
                                overlap = false;                                
                                while ( (boxadd = (TBox *)existing()) ) {
                                        if ( (realY1 >= boxadd->GetY1() && realY1 <= boxadd->GetY2()) ||
                                             (realY2 >= boxadd->GetY1() && realY2 <= boxadd->GetY2()) ) { overlap = true; break; }
                                }
                                if ( overlap == false ) { 
                                        TBox *box = dynamic_cast<TBox *>(boxnew->Clone());
                                        box->SetY1(realY1);
                                        box->SetY2(realY2);     
                                        fLastGates.Add(box);
                                } // new gate 
                                else { // change current box settings
                                        boxadd->SetX1( TMath::Min(boxadd->GetX1(),realY1) );
                                        boxadd->SetX2( TMath::Max(boxadd->GetX2(),realY2) );
                                }
                         }      
                } // y  
        } // while   
 
        // now merge if overlaps
        cout << " Nb Gates found in the current pad : " << fLastGates.GetSize() << endl; return fLastGates.GetSize(); 
}


void MatPlayer::ShowGates(Option_t *o)
{
//
// look for an histogram in the current pad. If one is found, it displays the boxes saved in fLastGates
//
        
        TObject *obj; TBox *box; TString opt = o;
        
        if ( TVirtualPad::Pad() == NULL ) return ;
        
        TIter next(TVirtualPad::Pad()->GetListOfPrimitives()); // look for an histogram
//      bool isTH2 = false;
        while ( (obj = next()) ) { 
                if ( obj->InheritsFrom("TH1") ) { 
             if ( obj->InheritsFrom("TH2") ) /*isTH2 = true;*/  break; }
        } // while
        
        TFrame *frame = TVirtualPad::Pad()->GetFrame(); 

        TIter existing(&fLastGates);
        while ( (box = (TBox *) existing()) ) { // 
                // Gate on X
                if ( IsBox(box,&fBoxX) ) {
                        box->DrawBox(box->GetX1(),frame->GetY1(),box->GetX2(),frame->GetY2());
                }
                if ( IsBox(box,&fBoxBX) ) {
                        box->DrawBox(box->GetX1(),frame->GetY1(),box->GetX2(),frame->GetY2());
                }
        }
}

void MatPlayer::Reset()
{
// to reset the conditions of gating
        TObject *obj; TList toremove;

        if ( TVirtualPad::Pad() == NULL ) return;

        TIter next(TVirtualPad::Pad()->GetListOfPrimitives());
        while ( (obj = next()) ) { // look for a box and keep a reference
                if ( obj->InheritsFrom("TBox") ) { toremove.Add(obj); }
        }
        TIter remove(&toremove);
        while ( (obj = remove()) ) { // remove lines
                        TVirtualPad::Pad()->GetListOfPrimitives()->Remove(obj);
        }
        toremove.Delete();

        TVirtualPad::Pad()->Update();
}

void MatPlayer::ProjectionY(Option_t *option)
{
// Gating on a gamma-gamma matrix
        TString opt = option;
  
        if ( fMat == NULL ) {
                cout << "Cannot Project in MatPlayer::ProjectionY, no matrix selected ! " << endl; return;
        }
        if ( InitGating("Gx") == 0 ) return; // init gating from boxes drawn in the current pad
  
// start gating by reseting first the histograms
        fPy->Reset("ICE"); fBy->Reset("ICE");
  
        bool isbg; Float_t w, wb; TBox *box; TH1D *h, *hb;
        
        isbg = false; w = wb = 0.0; h = hb = NULL;
        TIter next(&fLastGates);
        while ( (box = (TBox *)next()) ) { // make successive projections
                if ( IsBox(box,&fBoxX) )  { // projection
                        w += box->GetX2() - box->GetX1();
                        fPy->Add(  h = fMat->ProjectionY("MatPlayer_tmp_py",fPx->FindBin(box->GetX1()),fPx->FindBin(box->GetX2()) ) );
                }
                if ( IsBox(box,&fBoxBX) ) { // background               
                        isbg = true; wb += box->GetX2() - box->GetX1();
                        fBy->Add( hb = fMat->ProjectionY("MatPlayer_tmp_by",fPx->FindBin(box->GetX1()),fPx->FindBin(box->GetX2()) ) );
                }
        }  
        if ( h ) delete h; if ( hb ) delete hb;
        
        // do background substraction
    if ( isbg ) {
        fPy->Add(fBy,-w/wb);
        hb = dynamic_cast<TH1D *>( fMat->GetDirectory()->Get( Form("%s_by",fMat->GetName()) ) );
        if ( hb ) {
            hb->Reset("ICE");
            hb->Add(fBy);
            hb->Scale(w/wb);
        }
    }
    // move content of gated in local in global
    h = dynamic_cast<TH1D *>( fMat->GetDirectory()->Get( Form("%s_py",fMat->GetName()) ) );
    if ( h ) {
        h->Reset("ICE");
        h->Add(fPy);
    }
        
        fPy->Draw(); cout << " Projection on Y DONE \n ";
}


void MatPlayer::XEventAction(Int_t event, Int_t px, Int_t py, TObject */*sel*/)
{
// Do action depending on X actions
        
        TVirtualPad *pad = gROOT->GetSelectedPad(); if ( pad == NULL ) return;

        switch (event) {
                case kMouseMotion:
                        // to keep track of the last mouse position
                        fLastX = px; fLastY = py;
                        break;

                case kKeyPress:
                {
                        Int_t key = px, keysym = py;
                        switch ((EKeySym)keysym) {
                                case kKey_Left:
                                        cout<<"left arrow"<<endl;
                                        break;
                                case kKey_Right:
                                        cout<<"rightarrow"<<endl;
                                        break;
                                case kKey_Down:
                                        cout<<"down arrow"<<endl;
                                        break;
                                case kKey_Up:
                                        cout<<"up arrow"<<endl;
                                        break;
                                case kKey_x:
                                        // set a gate on the x axis
                                        SetGateX("Gx"); 
                                        break;
                                case kKey_X:
                                        // 
                                        break;                                  
                                case kKey_y:
                                        // set a gate at y position                             
                                        SetGateX("Gy");
                                        break;
                                case kKey_Y:
                                        ProjectionY(""); pad->Update();
                                        break;                                  
                                case kKey_b:
                                        SetGateX("Bx"); 
                                        break;
                                case kKey_B:
                                        SetGateX("By");                                 
                                        break;  
                                case kKey_h:
                                        cout << "Keys used to interact with MatPlayer: " << endl;
                                        cout << "\t x: add a gate on the x-axis " << endl;
                                        cout << "\t y: add a gate on the y-axis " << endl;
                                        cout << "\t b: background gate on the x-axis " << endl;
                                        cout << "\t B: background gate on the y-axis " << endl;                                 
                                        cout << "\t X: Processed X-Gating " << endl;
                                        cout << "\t Y: Processed Y-Gating " << endl;
                                        cout << "\t r: Remove all gates in the current pad " << endl;
                                        cout << "\t s: Show last gates in the current pad " << endl;
                                        break;  
                                case kKey_r:
                                        Reset(); pad->Update();
                                case kKey_s:
                                        ShowGates(""); pad->Update();
                                        break;  
                                default :
                                break;                                          
                        }       
                        if ( gDebug ) cout << "key " << key << "key " << keysym << endl;
                }
                break;
                
                default :
                break;
        } // switch
}

void MatPlayer::SubstractBG(const char *matname, Option_t *option)
{
        TString opt = option;
        
        // Test Option 
        if ( (!opt.Contains("RAD")) && (!opt.Contains("SNIP")) ) {
                cout << "Bad options - Choose 'RAD' or 'SNIP' " << endl; return;
        }
        
        TObject *obj = gROOT->FindObject(matname);
        if ( obj == NULL ) {
                cout << "Cannot find matrix named " << matname << endl;
        } else {        
                if ( obj->InheritsFrom("TH2") ) {
                        if ( opt.Contains("RAD") ) SubstractBG_RAD((TH2 *)obj);
                        if ( opt.Contains("SNIP")) SubstractBG_SNIP((TH2 *)obj);
                }
                else cout << "Object named "<< matname << " is not a matrix !! ";
        }
}

void MatPlayer::SubstractBG_RAD(TH2 *mat)
{
        Int_t i, j; TSpectrum worker; // use to substract a background on the x and y projection
        
        TH1D *proji, *projj;   
        Float_t *ptoti, *ptotj, *bgi, *bgj; Stat_t T;
        
// to get the projection on X and Y
        proji = mat->ProjectionX("tmp_rad_x"); projj = mat->ProjectionY("tmp_rad_y");
        
//Smooth projection
        proji->Smooth(2); projj->Smooth(2); 
        
// to normalize with the total number of counts in the matrix.
        T = mat->Integral();
        cout << " Number of counts in the gamma-gamma matrix: " << T << endl;
        
        ptoti = new Float_t[mat->GetNbinsX()]; bgi = new Float_t[mat->GetNbinsX()];
        ptotj = new Float_t[mat->GetNbinsY()]; bgj = new Float_t[mat->GetNbinsY()];
        for ( i = 0; i < mat->GetNbinsX(); i++ ) { ptoti[i] = bgi[i] = proji->GetBinContent(i+1); } 
        for ( j = 0; j < mat->GetNbinsY(); j++ ) { ptotj[j] = bgj[j] = projj->GetBinContent(j+1); }
        
        delete proji; delete projj;

        // worker.Background(bgi,mat->GetNbinsX(),40);worker.Background(bgj,mat->GetNbinsY(),40);                       
        // worker.Background1(bgi,mat->GetNbinsX(),40);worker.Background1(bgj,mat->GetNbinsY(),40);     
        //worker.Background1General(bgi,mat->GetNbinsX(),20,BACK1_INCREASING_WINDOW,BACK1_ORDER2,BACK1_INCLUDE_COMPTON);
        //worker.Background1General(bgj,mat->GetNbinsY(),20,BACK1_INCREASING_WINDOW,BACK1_ORDER2,BACK1_INCLUDE_COMPTON);
        
// now do the background substraction
        cout << " Background substraction on " << mat->GetName() << " by using the RADFORD procedure ..." << endl;
        
        Stat_t c, cmin = 0;
        for ( i = 0; i < mat->GetNbinsX(); i++ ) {
                for ( j = 0; j < mat->GetNbinsY(); j++ ) {
                
                        cout << " [" << i << "," << j << "] \r" ;
                        // printf(" [%d,%d] \r",i,j); 
                        
                        c = mat->GetBinContent(i+1,j+1);
                        c = c  - (ptoti[i]*ptotj[j] - (ptoti[i]-bgi[i])*(ptotj[j]-bgj[j]))/T;
                        if ( c < 0 && c < cmin ) { cmin = c; }
                                                          
                        mat->SetBinContent(i+1,j+1,c);
                } // j
        } // i  
        cout << " ... DONE. Lowest value in the background substracted matrix " << cmin << endl;
        
        delete [] ptoti; delete [] ptotj; delete [] bgi; delete [] bgj;
}

void MatPlayer::SubstractBG_SNIP(TH2 *mat)
{
// 2D Background substraction
// Based on the implementation given by M.Morhac NIM A,401 (1997) 113-132
// Based on SNIP procedure (Sensitive Nonlinear Iterative Peak)
        
        Int_t param = 3;
        
        cout << "*********************************************\n";
        cout << "* Background Substraction with SNIP Procedure\n";
        cout << "* Number of iterations : " << param << endl;
        cout << "*********************************************\n";
        cout << "* TH2 name : " << mat->GetName() << endl; 
        
        Int_t dimX, dimY;
        Float_t **v, **w, **y;
        
        Double_t P1, P2, P3, P4, S1, S2, S3, S4, a1, a2;
        
        
        dimX = mat->GetNbinsX();  dimY = mat->GetNbinsY();
        if (dimX <= 0) {cout << "Error, dimX <= 0\n"; return;} 
        if (dimY <= 0) {cout << "Error, dimY <= 0\n"; return;}  
         
        v = new Float_t*[dimX]; w = new Float_t*[dimX], y = new Float_t*[dimX];
        for (Int_t i=0 ; i<dimX ; i++){
                v[i] = new Float_t[dimY];
                w[i] = new Float_t[dimY];
                y[i] = new Float_t[dimY];
        }
        
        // Fill v
        for (Int_t i=0 ; i<dimX ; i++){
                for (Int_t j=0 ; j<dimY ; j++){
                        y[i][j] = mat->GetBinContent(i+1,j+1);
                        v[i][j] = mat->GetBinContent(i+1,j+1);
                }
        }
        
        
        // LLS operator
        for (Int_t i=0 ; i<dimX ; i++){
                for (Int_t j=0 ; j<dimY ; j++){
                        v[i][j] = TMath::Log( TMath::Log( TMath::Sqrt(y[i][j] + 1) + 1) + 1);
                        
                }
        }
        
        //SNIP Algorithm
        for (Int_t p=1 ; p<=param ; p++){
                for (Int_t i=p ; i<dimX-p ; i++){
                        for (Int_t j=p ; j<dimY-p ; j++){
                                P1 = v[i+p][j+p];
                                P2 = v[i-p][j+p];
                                P3 = v[i+p][j-p];
                                P4 = v[i-p][j-p];
                                
                                S1 = v[i+p][j];
                                S2 = v[i][j+p];
                                S3 = v[i][j-p];
                                S4 = v[i-p][j-p];
                                
                                S1 = TMath::Max(S1,(P1+P3)/2) - (P1+P3)/2;
                                S2 = TMath::Max(S2,(P1+P2)/2) - (P1+P2)/2;
                                S3 = TMath::Max(S3,(P3+P4)/2) - (P3+P4)/2;
                                S4 = TMath::Max(S4,(P2+P4)/2) - (P2+P4)/2;
                                
                                a1 = v[i][j];
                                a2 = (S1+S4)/2 + (S2+S3)/2 + (P1+P2+P3+P4)/4;
                                w[i][j] = TMath::Min(a1,a2);
                        } //j
                } //i
                for (Int_t i=p ; i<dimX-p ; i++){
                        for (Int_t j=p ; j<dimY-p ; j++){
                                v[i][j] = w[i][j];      
                        }
                }
                // printf ("[ %i ]\n",p);
                
                //TEST
                /*for (Int_t i=0 ; i<dimX ; i++){
                        for (Int_t j=0 ; j<dimY ; j++){
                                mat->SetBinContent(i+1,j+1,y[i][j]-v[i][j]);
                        } 
                }*/
                //ENDOFTEST
                
                //if (p==1) mat->ProjectionX("p=1")->Draw();
                //else mat->ProjectionX("p=1")->DrawClone("SAME");
                
                
        } //p
        
        
        // Inverse LLS operator
        for (Int_t i=0 ; i<dimX ; i++){
                for (Int_t j=0 ; j<dimY ; j++){
                        v[i][j] = TMath::Exp( TMath::Exp( v[i][j] - 1) -1 ) * TMath::Exp( TMath::Exp( v[i][j]-1) - 1) -1;
                        
                }
        }
        
        //fill mat
        for (Int_t i=0 ; i<dimX ; i++){
                for (Int_t j=0 ; j<dimY ; j++){
                        mat->SetBinContent(i+1,j+1,y[i][j]-v[i][j]);
                } 
        }  
        
        //Delete useless pointers
        for (Int_t i=0 ; i<dimX ; i++){
                delete [] v[i];
                delete [] w[i];
                delete [] y[i];
        }
        delete [] v;
        delete [] w;
        delete [] y;
        
}