TestFeeding.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 GW_TESTFEEDING_H
#include "TestFeeding.h"
#endif

#ifndef GW_NUCLEARLEVEL_H
#include "NuclearLevel.h"
#endif

#ifndef GW_GAMMALINK_H
#include "GammaLink.h"
#endif

#define MAX_E1 100
#define MAX_E2 100

using namespace Gw;

TestFeeding::TestFeeding()
{
        Link::SetIL(&fTop);
        
        fNbE1 = fNbE2 = 0; 
        fE1 = new TObjArray(MAX_E1); fE1->SetOwner(kTRUE); fE2 = new TObjArray(MAX_E2); fE2->SetOwner(kTRUE);
        
        // set E1 and E2 properties of the gammas
        GammaLink *gamma;
        
        for (Int_t i = 0; i < MAX_E1; i++ ) {
                gamma = new GammaLink(); gamma->SetEM("E1"); fE1->AddAt(gamma,i);       
        }
        for (Int_t i = 0; i < MAX_E2; i++ ) {
                gamma = new GammaLink(); gamma->SetEM("E2"); fE2->AddAt(gamma,i);       
        }
        // to get the single Rand instance
        fRand = Random::Instance();     
}

TestFeeding::~TestFeeding()
{
        // delete the tempory collection used 
        fNbE1 = fNbE2 = 0;  
        delete fE1; fE1 = NULL; delete fE2; fE2 = NULL; 
}
        
Int_t TestFeeding::DoCascade(TSeqCollection &col, Option_t */*opt*/)
{
        GammaLink *gamma;
        
        // determine what is the energy difference between the initial and the final level
        Float_t e, e_top, e_feeding, delta;
        
        e_top = fTop.GetEnergy().GetValue(); 
        e_feeding = ((NuclearLevel *)fFinal)->GetEnergy().GetValue();
        delta = e_top - e_feeding;
        
        if ( delta < 0.0 )
        return 0;
        
        // now starts the random process up which consist in adding alternatively an E2 and
        // an E2 transition to fill the energy gap delta.
        // E1 are choosen randomly between 0 and 8 MeV. E2 between 0 and 1 MeV
        bool is_E1;
        fNbE1 = 0; fNbE2 = 0;
        
        is_E1 = true; // true means an E1 has just been produced .. so it is E2 to be produced
        do {
                
                if ( fNbE1 > MAX_E1-1 || fNbE2 > MAX_E2-1 ) break; // TODO throw error
                e = fRand->Current()->Rndm();
                
                if ( is_E1 ) { // E2 to be produced

                        // set the proper energy
                        e = 8000 * e; 
                        if ( e > delta ) e = delta; // check out if the random process ends
                        
                        // set the gamma-ray properties and add it to the cascade
                        gamma = static_cast<GammaLink *>(fE2->At(fNbE2++));
                        gamma->GetEnergy().SetValue(e); col.Add(gamma);
                        
                        // next is E1
                        is_E1 = false;
                }
                else { // E1 to be produced

                        // set the proper energy
                        e = 8000 * e; 
                        if ( e > delta ) e = delta; // check out if the random process ends
                        
                        // set the gamma-ray properties and add it to the cascade
                        gamma = static_cast<GammaLink *>(fE1->At(fNbE1++));
                        gamma->GetEnergy().SetValue(e); col.Add(gamma);
                        
                        // next is E2
                        is_E1 = true;   
                }
                delta = delta - e;
        } while (delta > 0) ;
    
    return 1;
}