EnsdfLevelSchemeReader.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.             *
 ***************************************************************************/

#ifndef Gw_EnsdfLevelSchemeReader
#include "EnsdfLevelSchemeReader.h"
#endif

// root include's


#ifndef ROOT_TStyle
#include <TStyle.h>
#endif

#ifndef ROOT_TObjArray
#include <TObjArray.h>
#endif

#ifndef ROOT_TList
#include <TList.h>
#endif

// gw include's
#ifndef GW_ENSDF_H
#include "ENSDF.h"
#endif

#ifndef GW_LEVELSCHEME_H
#include "LevelScheme.h"
#endif

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

#ifndef GW_XGAMMALINK_H
#include "XGammaLink.h"
#endif


using namespace Gw;

//__________________________________________________________
EnsdfLevelSchemeReader::EnsdfLevelSchemeReader(const char* filename, Option_t* opt)
: BaseLevelSchemeReader(filename, opt)
{
    //Default Constructor
    fLog.GetProcessName() = "EnsdfLevelSchemeReader";
}

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


//__________________________________________________________
TString EnsdfLevelSchemeReader::GetDsid()
{
    TString dsid = fOptions.Data();
    
    if ( dsid.First(' ') > 0 ) { // at least two strings
        dsid.Remove(0,dsid.First(' ')+1);
        if ( dsid.Length() < 3 )
            dsid = "ADOPTED LEVELS, GAMMAS";
    } else {
        dsid = "ADOPTED LEVELS, GAMMAS";
    }
    return dsid;
}

//__________________________________________________________
TString EnsdfLevelSchemeReader::GetNuclide()
{
    TString nuclide = fOptions.Data();
    
    if ( nuclide.First(' ') > 0 ) { // at least two strings
        nuclide.Remove(nuclide.First(' '));
    } else {
        if (nuclide.Length() >= 5)
            nuclide.Remove(5);
        else
            nuclide.Remove(nuclide.Length());
    }
    return nuclide;
}


//__________________________________________________________
Int_t EnsdfLevelSchemeReader::Check()
{
    Int_t ok = 1;
    TString nuclide = GetNuclide();
    TString dsid    = GetDsid();
    
    fLog.SetProcessMethod("Check()");
    
    // open the file and check out if the dataset id exists
    ENSDF ens; unsigned int idset = 0;
    if ( ! ens.Open(fFileName.Data()) ) {
        fLog << error << "Cannot open file " << fFileName.Data() << nline ;
        ok = 0;
    } else {
        if ( (idset = ens.IsDataSet(nuclide.Data(),dsid.Data())) == 0 ) {
            fLog << error << "Cannot find dataset " << fOptions.Data()  << nline ;
            ok = 0;
        }
    }
    
    // if ok, it reads the data set and build the level scheme.
    Level *level; Cascade *cascade; Link *link;
    
    if ( ok == 1 ) { // a data set has been found, now the level scheme is built
        
        fLog << info << Form(" ---> Reading %s in dataset %s \n",nuclide.Data(),dsid.Data()) << nline ;
        
        // loop over the dataset and fill temporary collections.
        char *record = ens.BlankRecord(); bool is_record = ens.FirstRecord(idset,record); Int_t line = -1;
        while( is_record ) {
            
            // to count lines
            line++;
            
            TString tmp;
            tmp = "";
            //
            if ( (level = ens.GetLevel(record,tmp)) ) { // one new level
                fLevel.Add( level );
                // cout << " add level " << endl;
                // level->ls();
                // band treatment
                if ( tmp.Length() > 0 && tmp[0] != ' ' ) {
                    // look for a cascade with a name Energy[tmp[0]]
                    TIter iter(&fBand); TObject *obj; TString cname;
                    
                    cname = "Energy"; cname += tmp[0]; cascade = 0x0;
                    while ( (obj = iter() ) ) {
                        if ( cname == obj->GetName() ) {
                            cascade = (Cascade *)obj; break;
                        }
                    }
                    if ( cascade == 0x0 ) { // doesn't exist yet
                        cascade = new Cascade();
                        cascade->SetName(cname.Data()); fBand.Add(cascade);
                    }
                    // modify the level graphical position
                    Int_t whichband = fBand.IndexOf(cascade);
                    
                    // levels belonging to this cascade are moved at another x position
                    Double_t width;
                    width = TMath::Abs(level->GetX2() - level->GetX1());
                    
                    if ( whichband % 2 == 0 ) { // on the left
                        level->SetX1(level->GetX1()-1.1*width*(1+whichband/2));
                        level->SetX2(level->GetX2()-1.1*width*(1+whichband/2));
                    }
                    else { // on the right
                        level->SetX1(level->GetX1()+1.1*width*(1+whichband/2));
                        level->SetX2(level->GetX2()+1.1*width*(1+whichband/2));
                    }
                    
                    cascade->Add( level );
                    // printf("Energy offset %s \n ",tmp.Data());
                } // there is an unknown energy offset
                if ( tmp.Length() > 1 && tmp[1] != ' ' ) {
                    // look for a cascade with a name Spin[tmp[0]]
                    TIter iter(&fBand); TObject *obj; TString cname;
                    
                    cname = "Spin"; cname += tmp[1]; cascade = 0x0;
                    while ( (obj = iter() ) ) {
                        if ( cname == obj->GetName() ) {
                            cascade = (Cascade *)obj; break;
                        }
                    }
                    if ( cascade == 0x0 ) { // doesn't exist yet
                        cascade = new Cascade();
                        cascade->SetName(cname.Data()); fBand.Add(cascade);
                    }
                    // modify the level graphical position
                    // NOT FOR SPIN OFFSET, ONLY FOR ENERGY OFFSET
                    //Int_t whichband = fBand.IndexOf(cascade);
                    
                    // levels belonging to this cascade are moved at another x position
                    //Double_t width;
                    //width = TMath::Abs(level->GetX2() - level->GetX1());
                    
                    //if ( whichband % 2 == 0 ) { // on the left
                    //  level->SetX1(level->GetX1()-1.1*width*(1+whichband/2));
                    //  level->SetX2(level->GetX2()-1.1*width*(1+whichband/2));
                    //}
                    //else { // on the right
                    //  level->SetX1(level->GetX1()+1.1*width*(1+whichband/2));
                    //  level->SetX2(level->GetX2()+1.1*width*(1+whichband/2));
                    //}
                    cascade->Add( level );
                    // printf("Spin offset %s \n ",tmp.Data());
                } // there is an unknown spin offset
            }
            
            if ( (link = ens.GetLink(record,tmp)) ) { // one new link
                
                // look for the initial & final level attached with this link
                if ( fLevel.GetEntries() == 0 ) {
                    fLog << warning
                            << Form("Gamma with no initial level \n") << "\n from record at line " << line << nline ;
                    delete link;
                }
                else {
                    // last level is the initial level
                    NuclearLevel *initial, *final;
                    initial = (NuclearLevel*) fLevel.Last() ;
                    
                    link->SetIL( initial );
                    //cout << " add link " << endl;
                    //link->ls();
                    // Now look for the final level which gives an energy difference
                    // close to the gamma-ray energy
                    Float_t emin, eini, efin, egamma;
                    
                    egamma = ((GammaLink*)link)->GetEnergy().Get();
                    eini = initial->GetEnergy().Get();
                    emin = 2.0*egamma;
                    
                    //printf("gamma %f from level %f \n",egamma,eini);
                    
                    final = 0x0;
                    for (int i = 0; i < fLevel.GetEntries(); i++ ) { //
                        
                        final = (NuclearLevel*)fLevel.At(i);
                        if ( final == initial )
                            continue;
                        
                        efin = final->GetEnergy().GetValue();
                        if  ( TMath::Abs(eini-efin-egamma) < emin ) {
                            link->SetFL( final );
                            link->SetPoints();
                            emin = TMath::Abs(eini-efin-egamma) ;
                        }
                    } // i
                    if ( link->GetFL() == 0x0 ) {
                        fLog << warning
                                << Form("Gamma [record at line %d] from level %f, does not have a final level",line,eini) << nline ;
                        delete link;
                    } else {
                        // fLog << warning << Form("Gamma ok %f not linked with final level",egamma) << nline ;
                        fLink.Add( link );
                    }
                }
            }
            is_record = ens.NextRecord(idset,record);
        } // ok 
        
        delete [] record;
        
    } // if ok == 1
    
    fIsChecked = true;
        
    fLog << dolog;

    return ok;
}

//__________________________________________________________
void EnsdfLevelSchemeReader::FillLabels( LevelScheme& levelScheme)
{
    TString nuclide = GetNuclide();
    TString dsid    = GetDsid();  
    
    fLog << info << Form("  ... Reading dataset %s %s DONE", nuclide.Data(), dsid.Data()) << dolog;             
    // set name and reference
    levelScheme.SetName(nuclide.Data()); 
    dsid.Prepend("From ENSDF dataset "); 
    levelScheme.SetReference(dsid.Data());
}