RadLevelSchemeReader.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_RadLevelSchemeReader
#include "RadLevelSchemeReader.h"
#endif

// root include's 
#ifndef ROOT_TString
#include <TString.h>
#endif

#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;

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

}

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


//__________________________________________________________
Int_t RadLevelSchemeReader::Check()
{
  
  static const UInt_t kMaxLengthLine = 120;
  
  // logger message  
  fLog.SetProcessMethod("Check()");
  
  // open file to read data.
  FILE *file = fopen(fFileName.Data(),"r");
  if ( file == NULL ) {
    fLog << error << Form(" Cannot open file %s in LevelScheme::InitAGS \n", fFileName.Data()) << dolog;
    return 0;
  }
  
  // Initialisations from Radware
  Int_t j, j1; Char_t line[kMaxLengthLine];
  
  Int_t z, nblevel, nbgamma, nbband, nbtext;
  Float_t sizex, sizey, aspect_ratio, maxtan, maxdist, 
  cascadewidth, cascadesep, arroww, arrowl;
  
  Int_t band, flag, sflag, elflag;
  Float_t e, de, dxl, dxr, sldx, sldy, eldx, eldy;
  
  Char_t eol, jpi[16], kay[16];
  
  Int_t  li, lf, flg, elflg;
  Float_t intensity, dintensity, branching, conv, mix, dconv, dmix, dbranching, x1, x2;//,lamba;
  char em, clamba;
  
  // temporary collection to keep the numeric label associated in the AGS file 
  // with the gammas, the levels, the bands and the comments.
  // if everything goes right, the four collections will be emptied in the corresponding TList
  // in which the slot number is not necessary
        
//  TObjArray fBand, fLevel, fLink, fText;
//  fBand.SetOwner(kTRUE); fLevel.SetOwner(kTRUE); 
//  fLink.SetOwner(kTRUE); fText.SetOwner(kTRUE);
  
  TString lab; 
  NuclearLevel *level; Cascade *cascade; GammaLink *gamma;
  
  // read some comments (first five lines of file)
  //    printf("Level Scheme File header:\n");
  for (int j = 0; j < 5; ++j) {
    fgets(line, kMaxLengthLine, file);  // printf("%s", line);
  }
  fLog << info << Form(" ---> Reading file %s",fFileName.Data());
  fgets(line, kMaxLengthLine, file);
  if (strncmp(line, " ASCII GLS file format version 1.0", 34)) {
    fLog << error << "Cannot read ASCII gls file" << dolog; 
        fclose(file);
    fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
    return 0;
  }
  fgets(line, kMaxLengthLine, file);
  
  // read some global data      
  fgets(line, kMaxLengthLine, file);
  if (sscanf(line, "%i%i%i%i%i%f%f%f",
             &z,
             &nblevel, &nbgamma, &nbband, &nbtext,
             &sizex, &sizey, &aspect_ratio) != 8) {
    fLog << error << "Cannot read ASCII gls file" << dolog; 
        fclose(file);
    fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
    return 0;
  }

  fLog << info << Form("\t %d levels, %d gammas and %d bands supposed to be in this level scheme \n",
                       nblevel,nbgamma,nbband) << nline;;
  fLevel.Expand(nblevel);
        fLink.Expand(nbgamma); 
  fBand.Expand(nbband); 
        fText.Expand(nbtext);
  
  //    fName.SetTextSize(sizey); // default text size in LevelScheme's name
  
  fgets(line, kMaxLengthLine, file);
  fgets(line, kMaxLengthLine, file);
  if (sscanf(line, "%f%f%f%f%f%f",
             &maxtan, &maxdist,
             &cascadewidth, &cascadesep,
             &arroww, &arrowl) != 6) {
    fLog << error << "Cannot read ASCII gls file" << dolog; 
        fclose(file);
    fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
    return 0;
  }
  
  //    fLevelL.SetLineWidth();
  
  fgets(line, kMaxLengthLine, file);
  fgets(line, kMaxLengthLine, file);
  
  Float_t arrowbreak, levelbreak, lsl_csx, lsl_csy, lel_csx, lel_csy, gel_csx, gel_csy;
  if (sscanf(line, "%f%f%f%f%f%f%f%f",
             &arrowbreak, &levelbreak,
             &lsl_csx, &lsl_csy,
             &lel_csx, &lel_csy,
             &gel_csx, &gel_csy) != 8) {
    fLog << error << "Cannot read ASCII gls file" << dolog; 
        fclose(file);
    fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
    return 0;
  }
  fgets(line, kMaxLengthLine, file);
  fgets(line, kMaxLengthLine, file);
  
  // Reading each individual level
  
  fLog << debug << Form(" Reading level") << nline;
  for (j = 0; j < nblevel; j++) { // loop over nblevel
    
    fLog << debug << Form(" Reading level %d \n ",j) << nline;
    
    fgets(line, kMaxLengthLine, file);
    if (sscanf(line, "%*6c%f%f%8c%6c%i%i%i%i %c",
               &e, &de, jpi, kay,
               &band, &flag,
               &sflag, &elflag, &eol) != 9) {
      fLog << error << "Cannot read ASCII gls file" << dolog; 
          fclose(file);
      fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
      return 0;
    }
    
    band--;
    // printf(" jpi %s \n ",jpi);
    jpi[8] = '\0';
    
    if (eol == '&') {
      fgets(line, kMaxLengthLine, file);
      if (sscanf(line, "%*2c%f%f%f%f%f%f",
                 &sldx, &sldy,
                 &eldx, &eldy,
                 &dxl, &dxr) != 6) {
        fLog << error << "Cannot read ASCII gls file" << dolog; 
                fclose(file);
        fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
        return 0;
      }
    } else {
      sldx = sldy = 0.f;
      eldx = eldy = 0.f;
      dxl  = dxr  = 0.f;
    }
    
    // Init this new level and add it to the list of level
    fLog << debug << Form(" Initialisation level %d \n ",j) << nline;
    level = new NuclearLevel(); fLevel.AddAt(level,j);
    
    // In Radware, a level belongs to a band. In this package, a Cascade
    // is a given decay path. So a cascade is defined by a sequence of links.
    // Here, the tempband collection is used to store temporary the levels.
    // The final cascades will be composed a links one they are read.
    
    if ( fBand.At(band) == NULL ) { // doesn't exist 
      cascade = new Cascade(); fBand.AddAt(cascade,band);       
    } else { 
      cascade = (Cascade *)fBand.At(band); 
    }
    
    cascade->Add(level); // add this level to the band # band
    
    fLog << debug << Form(" Set Level properties \n ") << nline;                
    // flags for level
    if ( flag == 1 ) 
      level->SetStyle(Level::kTentative); // tentative
    if ( flag == 2 ) level->SetLineWidth(2*gStyle->GetLineWidth()); // isomer
    
    sscanf(jpi, "%i", &j1); 
    if (!strncmp(jpi + 5, "/2", 2)) {
      //                        printf(" spin %d \n ",j1);
      level->GetSpin().Set(j1,2);
      // spin = (float) j1 / 2.f;
    } else {
      //                        printf(" spin %d \n ",j1);
      level->GetSpin().Set(j1); 
      // spin = (float) j1;
    }
    level->SetColor(Level::kUnknown);
    if (jpi[7] == '-') {
      level->GetParity().Set(Parity::kMinus);
      level->SetColor(Level::kMinus);   
      // parity = -1.f;
    }
    else {
      level->GetParity().Set(Parity::kPlus);
      level->SetColor(Level::kPlus);    
      // parity = 1.f;
    }
    
    // flags for spin and parity
    fLog << debug << Form(" \t Spin & parity \n ") << nline;    
    if ( sflag == 0 ) { // both are known
      level->GetSpin().SetInfo(InfoData::kKnown);
      level->GetParity().SetInfo(InfoData::kKnown);
      level->SetLabel(jpi,0);   
    }
    if ( sflag == 1 ) { // spin is tentative, parity known
      level->GetSpin().SetInfo(InfoData::kTentative);   
      level->GetParity().SetInfo(InfoData::kKnown);     
      lab = jpi; lab.Insert(6,")"); lab.Insert(0,"("); level->SetLabel(lab.Data(),0);           
    }           
    if ( sflag == 2 ) { // spin is tentative, parity known
      level->GetSpin().SetInfo(InfoData::kKnown);       
      level->GetParity().SetInfo(InfoData::kTentative);                 
      lab = jpi; lab.Insert(8,")"); lab.Insert(6,"("); level->SetLabel(lab.Data(),0);           
    }           
    if ( sflag == 3 ) { // both tentative 
      level->GetSpin().SetInfo(InfoData::kTentative);   
      level->GetParity().SetInfo(InfoData::kTentative);                 
      lab = jpi; lab.Insert(8,")"); lab.Insert(0,"("); level->SetLabel(lab.Data(),0);           
    }           
    if ( sflag == 4 ) { // spin known, parity unknown
      level->GetSpin().SetInfo(InfoData::kKnown);       
      level->GetParity().SetInfo(InfoData::kUnknown);                   
      lab = jpi; lab.Remove(7,1); level->SetLabel(lab.Data(),0);                
    }           
    if ( sflag == 5 ) { // spin tentative, parity unknown
      level->GetSpin().SetInfo(InfoData::kTentative);   
      level->GetParity().SetInfo(InfoData::kUnknown);                   
      lab = jpi; lab.Remove(7,1); lab.Insert(6,")"); lab.Insert(0,"("); level->SetLabel(lab.Data(),0);          
    }           
    if ( sflag == 6 ) { // both are unknown
      level->GetSpin().SetInfo(InfoData::kUnknown);
      level->GetParity().SetInfo(InfoData::kUnknown);   
      lab = "-?"; level->SetLabel(lab.Data(),0);        
    }
    //K quantum number - not yet implemented
    //    sscanf(kay, "%i", &j1);
    //     if (!strncmp(kay + 4, "/2", 2)) {
    //     k = (float) j1 / 2.f;
    //     } else {
    //     k = (float) j1;
    //     }
     
    fLog << debug << Form(" Set labels \n ") << nline;  
    
    // spin label
    //level->GetLabel(0).SetTextSize(lsl_csy);
    level->GetLabel(0).SetX(sldx); level->GetLabel(0).SetY(sldy);
    //level->GetLabel(0).SetTextSize(sizey);level->GetLabel(0).SetX(sldx);level->GetLabel(0).SetY(sldy);
    
    // for isomers
    //level->GetLabel(2).SetTextSize(sizey); 
    
    // energy label
    //level->GetLabel(0).SetTextSize(lel_csy);
    // latex->SetX(eldx);latex->SetY(eldy);
    
    // additional label
    //latex->SetTextSize(sizey);        
    
    fLog << debug << Form(" Set level offset \n ") << nline ;   
    // set additional offset from the default position.
    fLog << debug << Form(" offset left side %f \n ",dxl) << nline ;    
    level->SetX1(dxl);
    fLog << debug << Form(" offset right side %f \n ",dxr) << nline ;   
    level->SetX2(dxr);          
    fLog << debug << Form(" vertical offset %f \n ",e) << nline ;       
    level->SetY1(e); 
    fLog << debug << Form(" vertical  offset %f \n ",e) << nline ;      
    level->SetY2(e);    
    
    level->GetEnergy().Set(e,de);                       
  } // j < nblevels
  fLog << info << Form("\t %d levels correctly read in %d bands \n",
                       fLevel.GetEntries(),fBand.GetEntries()) << nline ;
  
  fgets(line, kMaxLengthLine, file);
  
  // Start reading each band
  char name[8]; Float_t x0, nx;
  for (j = 0; j < nbband; ++j) {
    
    fgets(line, kMaxLengthLine, file);
    if (sscanf(line, "%*6c%8c%f%f%f%f%f%f",
               name, &x0, &nx,
               &sldx, &sldy,
               &eldx, &eldy) != 7) {
      fLog << error << "Cannot read ASCII gls file" << dolog; 
          fclose(file);
      fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
      return 0;
    }
    
    // in case a band exist without levels in it.
    if ( fBand.At(j) == 0 ) fBand.AddAt(new Cascade(),j);
    cascade = (Cascade *)fBand.At(j);
    
    // set cascade's name and position
    TLatex &lat = cascade->GetLatex("name"); 
    TString tmp(name);
    while (tmp[0] == ' ') tmp.Remove(0,1); // remove blanks in name
    Float_t textSize = sizey/(nx*5);
    lat.SetText(x0, sizey,tmp.Data()); lat.SetTextSize(textSize);
    Level::SetDefaultLabelSize(textSize);
    Link::SetDefaultLabelSize(textSize);

    // now for this band, compute the level position and width
    for( j1 = 0; j1 < cascade->GetSize() ; j1++){
      
      level = (NuclearLevel *)cascade->At(j1);
      //printf("x0 = %f ; nx/2 = %f ; level->GetX1()=%f ; SetX1 = %f \n",x0,nx/2,level->GetX1(),x0 - nx/2 + level->GetX1());
      //printf("x0 = %f ; nx/2 = %f ; level->GetX2()=%f ; SetX2 = %f \n",x0,nx/2,level->GetX2(),x0 + nx/2 + level->GetX2());
      level->SetX1( x0 - nx/2 + level->GetX1() );
      level->SetX2( x0 + nx/2 + level->GetX2() );
      level->SetDefaultLabels();
    } // j1
  } // nbbands
  fLog << info << Form("\t %d bands correctly read \n",fBand.GetSize()) << nline;
  
  // Start reading each individual gamma
  
  fgets(line, kMaxLengthLine, file);
  fgets(line, kMaxLengthLine, file);
  fgets(line, kMaxLengthLine, file);
  
  for (j = 0; j < nbgamma; ++j) {
    
    // read informations from the files
    fgets(line, kMaxLengthLine, file);
    eol = '&';
    intensity = dintensity = 0.0f;
    if (sscanf(line, "%*6c%f%f%*3c%c%*c%c%i%i%f%f %c",
               &e, &de, &em, &clamba,
               &li, &lf,
               &intensity, &dintensity, &eol) < 6) {
      fLog << error << "Cannot read ASCII gls file" << dolog; 
          fclose(file);
      fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
      return 0;
    }
    
    if (clamba == ' ') clamba = '0';
//    lamba = (float) (clamba - '0');
    li--;
    lf--;
    if (eol == '&') {
      fgets(line, kMaxLengthLine, file);
      if (sscanf(line, "%*2c%f%f%f%f%f%f %c",
                 &conv, &dconv,
                 &branching, &dbranching,
                 &mix, &dmix, &eol) != 7) {
        fLog << error << "Cannot read ASCII gls file" << dolog; 
                fclose(file);
        fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
        return 0;
      }
    } else {
      conv  = dconv  = 0.0f;
      branching = dbranching = 0.0f;
      mix  = dmix  = 0.0f;
      //                        calc_alpha(j);
    }
    if (eol == '&') {
      fgets(line, kMaxLengthLine, file);
      if (sscanf(line, "%*2c%f%f%f%f%i%i",
                 &x1, &x2,
                 &eldx, &eldy,
                 &flg, &elflg) != 6) {
        fLog << error << "Cannot read ASCII gls file" << dolog; 
                fclose(file);
        fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
        return 0;
      }
    } else {
      x1 = x2 = 0.f;
      eldx = eldy = 0.f;
      flg = elflg = 0;
    }           
    if ( li >= fLevel.GetEntries() ) {
      fLog << error << Form("- ! - gamma-ray %d associated with a non-existing initial level %d \n",j,li) << dolog;
      fclose(file);
      fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
      return 0; 
    }
    if ( lf >= fLevel.GetEntries() ) {
      fLog << error << Form(" - ! - gamma-ray %d associated with a non-existing final level %d \n",j,lf) << dolog;
      fclose(file);
      fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
    }
    
    // now creates a new gamma link
    if( conv > 0 ){ // a converted gamma-ray
      //                        gamma = new XGammaLink();               
      gamma = new GammaLink();          
    } else {
      gamma = new GammaLink();
    }
    if ( gamma == NULL ) 
                fLog << error << " Pointer is NULL" << nline ;
    
    fLink.AddAt(gamma,j);
    // set characteristics
    gamma->GetEnergy().Set(e,de);
    if( conv > 0 ) { 
      // conv = 0.0;
      //gamma->GetStrength().Set(intensity*(1+conv),
      //                         intensity*(1+conv)*((dintensity/intensity)+(dconv/conv)));             
      gamma->GetStrength().Set(intensity,dintensity);                                   
    } else {
      gamma->GetStrength().Set(intensity,dintensity);
    }
    gamma->GetMixing().Set(mix,dmix);
    gamma->GetConversion().Set(conv,dconv);
    
    lab = em + clamba;
    gamma->SetEM(lab.Data());
    gamma->SetIL((Level *)fLevel.At(li));
    gamma->SetFL((Level *)fLevel.At(lf));       
    gamma->SetDefaultLabels();
    // to draw the link
    if ( gamma->GetIL() != NULL || gamma->GetFL() != NULL ) {
      Double_t xi = (gamma->GetIL()->GetX2()+gamma->GetIL()->GetX1())/2.0;
      Double_t xf = (gamma->GetFL()->GetX2()+gamma->GetFL()->GetX1())/2.0;
      if ( x1 == 0.00 && x2 != 0.00 ) {
        xi = (x2 - TMath::Abs(gamma->GetIL()->GetX2()-gamma->GetIL()->GetX1()))/2.0; 
        xf = (x2 - TMath::Abs(gamma->GetIL()->GetX2()-gamma->GetIL()->GetX1()))/2.0;        
      }
      if ( x1 != 0.00 && x2 != 0.00 ) { 
        xi = x1 - TMath::Abs(gamma->GetIL()->GetX2()-gamma->GetIL()->GetX1()) / 2.0; 
        xf = x2 - TMath::Abs(gamma->GetIL()->GetX2()-gamma->GetIL()->GetX1()) / 2.0;    
      }
      gamma->SetPoints(xi,xf); 
      if ( xi > gamma->GetIL()->GetX2() ) gamma->GetIL()->SetExtraLineX2(xi);
      if ( xi < gamma->GetIL()->GetX1() ) gamma->GetIL()->SetExtraLineX1(xi);
      
      if ( xf > gamma->GetFL()->GetX2() ) gamma->GetFL()->SetExtraLineX2(xf);
      if ( xf < gamma->GetFL()->GetX1() ) gamma->GetFL()->SetExtraLineX1(xf);     
    }
    else { 
      fLog << warning << Form("WARNING one gamma %f is not linked \n", gamma->GetEnergy().Get()) << nline ; 
    } 
  } // j - gamma
  fLog << info << Form("\t %d gammas correctly read \n",fLink.GetEntries()) << nline;
  
  // now read additional text to be written on the level scheme
  Char_t text[40]; Int_t dumi; Float_t size, x, y;
  if (nbtext > 0) {
    fgets(line, kMaxLengthLine, file);
    fgets(line, kMaxLengthLine, file);
    for (j = 0; j < nbtext; ++j) {
      fgets(line, kMaxLengthLine, file);
      if (sscanf(line, "%*7c%40c%i", text, &dumi) != 2) {
        fLog << error << "Cannot read ASCII gls file" << dolog ; 
                fclose(file);
        fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
        return 0;
      }
      
      fgets(line, kMaxLengthLine, file);
      if (sscanf(line, "%*2c%f%f%f%f",
                 &size, &size,&x, &y) != 4) {
        fLog << error << "Cannot read ASCII gls file" << dolog; 
                fclose(file);
        fBand.Clear(); fLevel.Clear(); fLink.Clear(); fText.Clear();
        return 0;
      }
      
      TLatex *latex = new TLatex(x,y,text); latex->SetTextSize(size); fText.AddAt(latex,j);
    } // j - labels
    fLog << info << Form("\t %d extra comments correctly read \n",fText.GetEntries()) << nline;
  } // if nbtext        
  fclose(file);
  fLog << info << Form("  ... Reading file %s DONE \n", fFileName.Data()) << nline;
  fLog << dolog; 

  fIsChecked = true;
  return 1;
}

 
//__________________________________________________________
void RadLevelSchemeReader::FillLabels( LevelScheme& levelScheme)
{  
  // set name and reference
  levelScheme.SetName(fOptions.Data()); 
  TString lab = "From Radware AGS file "; lab += fFileName.Data(); 
  levelScheme.SetReference(lab.Data());
  
}