RandObj.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004 by Olivier Stezowski                               *
 *   stezow(AT)ipnl.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.             *
 ***************************************************************************/

#ifndef  ROOT_TRandom
#include <TRandom.h>
#endif
#ifndef  ROOT_TH1D
#include <TH1D.h>
#endif

#ifndef GW_RANDOM_H
#include "Random.h"
#endif

#ifndef GW_RANDOBJ_H
#include "RandObj.h"
#endif

#include <iostream>

// namespace to be used
using namespace std;
using namespace Gw;  

ClassImp(RandObj);
 
RandObj::RandObj()  { fRand = Random::Instance(); fObj.SetOwner(kFALSE); }
RandObj::~RandObj() { fObj.Clear();}

void RandObj::Add(TObject *obj, Float_t weight)
{
        Int_t which; 
        
        // check if obj already in the collection
        which = fObj.IndexOf(obj);
        
        if (which == -1 ) { // new object
                if ( weight < 0 ) 
                        {fObj.Add(obj); fWeight.push_back(-weight);}
                else 
                        {fObj.Add(obj); fWeight.push_back(weight); }
        }
        else { // obj already in the collection
                if ( (fWeight[which] + weight) > 0 ) fWeight[which] = fWeight[which] + weight; 
                else {
                  fObj.Remove(obj); fObj.Compress(); fWeight.erase(fWeight.begin()+which);
                }
        }
        ComputeCumu();
}

void RandObj::ComputeCumu()
{
        unsigned int i;
        
        fCumu.clear(); if ( fObj.GetEntries() == 0 ) return; // in case of empty collection
        
        // new size
        fCumu.resize(fObj.GetEntries());
        
        // now compute the cumulative distribution
        fCumu[0] = fWeight[0];
        for( i = 1; i < fCumu.size(); i++ ){ fCumu[i] = fCumu[i-1] + fWeight[i]; }
        
        Double_t sum = fCumu[fCumu.size()-1];
        for( i = 0; i < fCumu.size() ; i++ ){ if ( sum != 0.0 ) fCumu[i] = fCumu[i] / sum; }
}

void RandObj::Clear(Option_t *options)
{
        fObj.Clear(options); fWeight.clear(); fCumu.clear();
}
    
TObject * RandObj::Rand()
{
        Int_t sup, inf, middle; inf = 0; sup = fCumu.size(); 
        
        if ( sup == 1 ) return fObj[0]; 
                
        Double_t r = fRand->Current()->Rndm();
        
        while( (sup-inf) > 1) {
                middle = (sup+inf)/2;
                if (r == fCumu[middle-1]) return fObj[middle-1];
                if (r  < fCumu[middle-1]) 
                        sup = middle;
                else 
                        inf = middle;
        }
        return fObj[inf];
}


TObject * RandObj::Rand(Int_t &slot)
{
        Int_t sup, inf, middle; inf = 0; sup = fCumu.size(); 
        
        if ( sup == 1 ) { slot = 0; return fObj[0]; }
        
        Double_t r = fRand->Current()->Rndm();

        while( (sup-inf) > 1) {
                middle = (sup+inf)/2;
                if (r == fCumu[middle-1]) { slot = middle-1; return fObj[middle-1]; }
                if (r  < fCumu[middle-1]) 
                          sup = middle;
                else 
                          inf = middle;  
        }
        slot = inf; return fObj[slot];
}

void RandObj::ls(Option_t* /*opt*/) const
{
 cout << " # of objects " << fObj.GetEntries() 
      << " dim cumulative " << fCumu.size()
      << " dim cumulative " << fWeight.size() << endl;
}

void RandObj::FillRandom(TH1 *h, Int_t ntimes)
{
        Int_t slot; 
         
        if ( fObj.GetEntries() > 0 && ntimes > 0 ) {
                h = new TH1D("RandObj","Random distribution from RandObj::FillRandom",fObj.GetEntries(),0,fObj.GetEntries());
                for (Int_t i = 0; i < ntimes;  i++ ) 
                        { RandObj::Rand(slot); h->Fill(slot); }
        }
        else h = NULL;
}