EvapLMOF.cpp

Go to the documentation of this file.

#ifndef ROOT_TMath
#include "TMath.h"
#endif

#ifndef ROOT_TVector3
#include "TVector3.h"
#endif

#ifndef ROOT_TLeafS
#include "TLeafS.h"
#endif
#ifndef ROOT_TLeafI
#include "TLeafI.h"
#endif
#ifndef ROOT_TLeafF
#include "TLeafF.h"
#endif
#ifndef ROOT_TLeafB
#include "TLeafB.h"
#endif

#ifndef ROOT_TRandom
#include "TRandom.h"
#endif

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

#include "EvapLMOF.h"

TRandom *therand;

// anonymous namespace to put EvapOR import that come usiing f2c for some fortran functions of the package
//
// plm2.f, getang.f
//
namespace {

/* Common Block Declarations */
struct xqang_ { Float_t sum[15300] /* was [300][51] */; } xqang_1;

/* Table of constant values */
static Double_t c_b4 = 1e-15;

/* Subroutine */ Int_t plm2_(void)
{
    /* Initialized data */
    static Float_t bin[3] = { .1667f,.5f,.8333f };

    /* System generated locals */
    Long_t i__1;
    Double_t d__1;

    /* Builtin functions */
    // Double_t sqrt(Double_t), d_sign(Double_t *, Double_t *);

    /* Local variables */
    static Long_t l, m, n;
    static Double_t z__;
    static Long_t ii, ik, jk, ll, mm;
    static Double_t pl[900]     /* was [30][30] */;
    static Long_t iang;
    static Double_t temp, plmm, dsum[15300]     /* was [300][51] */;
    static Long_t iiter;

/*     CALCULATE ASSOCIATED LEGENDRE FUNTIONS AND FORM TABLE OF RUNNING */
/*     DSUMS FOR LATER SELECTION OF ANGLE. */
/*     TO SAVE SPACE THE RUNNING DSUM  DSUM(N,IANG)  IS STORED FOR EACH */
/*     P(L,M) AT  N=L*(L-1)/2+M . THE POLYNOMIAL IS OF ORDER (L-1,M-1). */
/*     IANG IS A COS(ANGLE) INDEX FROM .05 TO .95 IN 50 STEPS. */
/* ******** WARNING. FOR  M = L  AND L ABOVE 20 RESULTS ARE INACCURATE. */
/*     THIS IS A DUPLICATE OF THE SUBROUTINE PLM, IN THE MAIN PROGRAM */
/*     A FINER GRANULARITY IN THE ANGLE RANGE IS REQUIRED HERE. */

/*     THIS FORMS FACTORIALS FOR LATER  C.G.  CALCULATIONS */
    for (ii = 1; ii <= 3; ++ii) {
        bin[ii - 1] *= .02f;
    }
    for (ik = 1; ik <= 300; ++ik) {
        for (jk = 1; jk <= 50; ++jk) {
/* L10: */
            dsum[ik + jk * 300 - 301] = 0.f;
        }
    }
    for (iang = 1; iang <= 50; ++iang) {
        for (iiter = 1; iiter <= 3; ++iiter) {
            z__ = iang * .02f - bin[iiter - 1];
            pl[0] = 1.f;
            pl[1] = z__;
/* Computing 2nd power */
            d__1 = z__;
            temp = 1.f - d__1 * d__1;
            pl[31] = TMath::Sqrt(temp);
            if (TMath::Abs(pl[31]) < 1e-15) {
                pl[31] = TMath::Sign(c_b4, pl[31]);
            }
/* Computing 2nd power */
            d__1 = z__;
            pl[2] = (d__1 * d__1 * 3.f - 1.f) * .5f;
            pl[32] = pl[31] * 3.f * z__;
            pl[62] = temp * 3.f;
            for (l = 4; l <= 24; ++l) {
                ll = l - 1;
                pl[l - 1] = (((ll << 1) - 1) * z__ * pl[l - 2] - (ll - 1) * 
                        pl[l - 3]) / ll;
                pl[l + 29] = (((ll << 1) - 1) * z__ * pl[l + 28] - ll * pl[l 
                        + 27]) / (ll - 1);
                i__1 = l;
                for (m = 3; m <= i__1; ++m) {
                    mm = m - 1;
                    pl[l + m * 30 - 31] = (ll + mm - 1) * pl[l - 1 + (m - 1) *
                             30 - 31] - (ll - mm + 1) * z__ * pl[l + (m - 1) *
                             30 - 31];
                    pl[l + m * 30 - 31] /= pl[31];
/*     INDEX OF (L,M) = N */
/* L20: */
                }
            }
            for (l = 1; l <= 24; ++l) {
                i__1 = l;
                for (m = 1; m <= i__1; ++m) {
                    plmm = pl[l + m * 30 - 31] * 1e-10f;
                    if (l > 11 && m > (l - 1) / 2) {
                        plmm = pl[l + m * 30 - 31] * 1e-20f;
                    }
                    n = l * (l - 1) / 2 + m;
/* L30: */
/* Computing 2nd power */
                    d__1 = plmm;
                    dsum[n + iang * 300 - 301] += d__1 * d__1;
                }
            }
/* L40: */
        }
    }

    for (n = 1; n <= 300; ++n) {
        for (iang = 2; iang <= 50; ++iang) {
/* L50: */
            dsum[n + iang * 300 - 301] += dsum[n + (iang - 1) * 300 - 301];
        }
        if (dsum[n + 14699] == 0.f) {
            dsum[n + 14699] = 1e-20;
        }
        for (iang = 1; iang <= 50; ++iang) {
            dsum[n + iang * 300 - 301] /= dsum[n + 14699];
/* L60: */
            xqang_1.sum[n + iang * 300 - 301] = dsum[n + iang * 300 - 301];
        }
    }

    return 0;
} /* plm2_ */

/* Subroutine */ Int_t getang_(Int_t *l, Int_t *m, Float_t *ct, Float_t *theta, Float_t *phi)
{
    /* Initialized data */
    static Short_t isetup = 0;

    /* Builtin functions */
    // double Double_t acos(Double_t);

    /* Local variables */
    static Int_t n;
    static Float_t x;
    static Int_t im, ict;
//    static Int_t iogr;

/*  GETANG chooses an angle for emission randomly based on the (l,m) */
/*   values from pace. igrid gives the number of cos theta bins (in 0-90 */
/*   range. */
    if (isetup == 0) {
        isetup = 1;
        plm2_();
        therand = new TRandom(0); // replace raninit_rtc__(&c__2) as a random generator 
//      iogr = 50;
    }
    im = TMath::Abs(*m);
    n = (*l) * (*l + 1) / 2 + im + 1;
    x = therand->Rndm(); // ranf_(&idummy);
    for (ict = 1; ict <= 50; ++ict) {
        *ct = (Float_t) ict;
        if (x <= xqang_1.sum[n + ict * 300 - 301]) {
            goto L100;
        }
    }
L100:
//    *ct = (*ct - ranf_(&idummy)) * .02f;
    *ct = (*ct - therand->Rndm()) * .02f;
    x = therand->Rndm(); // ranf_(&idummy);
    if (x < .5f) {
        *ct = -(*ct);
    }
//    *phi = ranf_(&idummy) * 6.2831853999999998f;
    *phi = therand->Rndm() * 6.2831853999999998f;
    *theta = TMath::ACos(*ct);
    return 0;
} /* getang_ */

} // anonymous namespace for EvapOR like functions


using namespace std;
using namespace Gw;

Short_t EvapLMOF::VERBOSE = 0;  
const Short_t EvapLMOF::MAXCASCADE; const Short_t EvapLMOF::MAXPARTICLES; 
        
EvapLMOF::EvapLMOF() : fHeader(NULL), fGlobalRecord(NULL), fIdRecord(NULL), fDataRecord(NULL), fIdOnly(false), fInfile(NULL), fCurrent(-1), fMaxFiles(0), fStatus(kEmpty)
{
        SetPath();
}

EvapLMOF::~EvapLMOF()
{
        if ( EvapLMOF::VERBOSE ) cout << " IN EvapLMOF::~EvapLMOF()" << endl;
        
        if ( fHeader != NULL ) { 
                if ( fHeader->IsStatus(Gws::Buffer::kCorrupt) && EvapLMOF::VERBOSE ) cout << " Header buffer corrupted " << endl;
                delete fHeader; fHeader = NULL;
        }
        if ( fGlobalRecord != NULL ) { 
                if ( fGlobalRecord->IsStatus(Gws::Buffer::kCorrupt) && EvapLMOF::VERBOSE ) cout << " Global buffer corrupted " << endl; 
                delete fGlobalRecord; fGlobalRecord = NULL; 
        }
        if ( fIdRecord != NULL ) { 
                if ( fIdRecord->IsStatus(Gws::Buffer::kCorrupt) && EvapLMOF::VERBOSE ) cout << " Id buffer corrupted " << endl;         
                delete fIdRecord;  fIdRecord= NULL; 
        }
        if ( fDataRecord != NULL ) { 
                if ( fDataRecord->IsStatus(Gws::Buffer::kCorrupt) && EvapLMOF::VERBOSE ) cout << " Data buffer corrupted " << endl;                     
                delete fDataRecord; fDataRecord = NULL; 
        }       
        
        if ( fInfile ) { ::fclose(fInfile); fInfile = NULL; } 
        
        if ( EvapLMOF::VERBOSE ) cout << " OUT EvapLMOF::~EvapLMOF()" << endl;
}
        
void EvapLMOF::ShowCurrentConditions(std::ostream &out)
{
        out << "--> Path to data:\n\t" << fPath << endl;
        out << "--> Chain of files: " << endl;
        for (int i = 0; i < fMaxFiles; i++ ) {
                if ( i == fCurrent ) cout << " [c] ";
                else out << " [ ] ";
                out << fEventsRead[i] << " events read in " << fFileName[i] << endl;
        }
        if ( fWarning != "" )   out << " Warning: " << fWarning.Data() << endl;
        if ( fStatus == kFail ) out << " Last Error: " << fError.Data() << endl;
        if ( fStatus == kEoE )  out << " All events have been read " << endl;   
}

void EvapLMOF::ShowCurrentEvent(std::ostream &out)
{
        Short_t i;
        
        out << " -- GLOBAL RECORD " << endl;
        out << "AC: " << fGlobal.AC << ", ZC: " << fGlobal.ZC << ", NCASC: " << fGlobal.NCASC << ", MDIR: " << fGlobal.MDIR
             << ", JCAS: " << fGlobal.JCASC << ", ENER: " << fGlobal.ENER << ", SIGMA: " << fGlobal.SIGMA 
             << ", VRC: "  << fGlobal.VRC << ", MAX_P: "  << fGlobal.MAX_P << endl; 
        out << " -- ID RECORD " << endl;        
        out << "MG: " << fId.MG << ", EE: " << fId.EE << ", JE: " <<  fId.JE << ", L: " << fId.L
             << ", Z: " << fId.Z << ", N: " << fId.N << ", NUM: " << fId.NUM 
             << ", WT: "  << fId.WT << ", FIS: "  <<  fId.FIS << endl;  
             
        if ( fIdOnly ) return; // no need to show the full cascade in idonly mode
        
        out << " -- DATA RECORD " << endl;
        for ( i = 0; i < fId.NUM; i++ ) { // cascade
                out << "Cascade " << i << endl << "MODE: " << fData.MODE[i] << ", EI: " << fData.EI[i] << ", JI: " << fData.JI[i] << ", EP: " <<  fData.EP[i] << ", LP: " << fData.LP[i] << ", MP: " << fData.MP[i] << endl;
                out << "VE0: " << fKin.VE[i][0] << ", VE1: " << fKin.VE[i][1] << ", VE2: " << fKin.VE[i][2] << ", THETA: " << fKin.THETA[i] << ", PHI: " << fKin.PHI[i] << ", EEJ: " << fKin.EEJ[i] << ", JF: " << fKin.JF[i] << endl;  
        }         
        for (  i = 1; i < MAXPARTICLES; i++ ) { // statitics per cascade
                if ( fData.MUL[i] > 0 ) out << " MODE " << i << " MULT: " << fData.MUL[i] << ", ES: " << fData.ES[i] << ", DJ: " << fData.DJ[i] << endl;
        }  
        out << "ZB: " << fKin.ZB << ", AB: " << fKin.AB << ", ER: " << fKin.ER  << ", VO0: " << fKin.V0[0] << ", VO1: " << fKin.V0[1] << ", VO2: " << fKin.V0[2] << endl;
}

void EvapLMOF::SetRecord(Gws::Memory::EEndian e)
{
        if ( EvapLMOF::VERBOSE ) cout << " IN EvapLMOF::SetRecord(Gws::Memory::EEndian e)" << endl;
        
        if ( fHeader != NULL ) { // already exist
                if ( fHeader->IsBytes(e) ) { // with the good endian type
                        fHeader->Reset();
                }
                else { delete fHeader;  fHeader = NULL; } 
        }
        if ( fHeader == NULL ) { fHeader = Gws::Buffer::New(e,4); }
        
        if ( fGlobalRecord != NULL ) { // already exist
                if ( fGlobalRecord->IsBytes(e) ) { // with the good endian type
                        fGlobalRecord->Reset();
                }
                else { delete fGlobalRecord; fGlobalRecord = NULL; } 
        }
        if ( fGlobalRecord == NULL ) { fGlobalRecord = Gws::Buffer::New(e,40+4); } 
        // 40+4 because of the additionnal int at the end that corresponds to the length of the buffer
                
        if ( fIdRecord != NULL ) { // already exist
                if ( fIdRecord->IsBytes(e) ) { // with the good endian type
                        fIdRecord->Reset();
                }
                else { delete fIdRecord; fIdRecord = NULL; } 
        }
        if ( fIdRecord == NULL ) { fIdRecord = Gws::Buffer::New(e,14+4); }      
        // 14+4 because of the additionnal int at the end that corresponds to the length of the buffer
                
        if (  fDataRecord != NULL ) { // already exist
                if ( fDataRecord->IsBytes(e) ) { // with the good endian type
                        fDataRecord->Reset();
                }
                else { delete fDataRecord; fDataRecord = NULL; } 
        }
        if ( fDataRecord == NULL ) { fDataRecord = Gws::Buffer::New(e,8*MAXCASCADE); }  
        
        if ( EvapLMOF::VERBOSE ) cout << " OUT EvapLMOF::SetRecord(Gws::Memory::EEndian e)" << endl;
}

void EvapLMOF::ReadHeader()
{       
        if ( EvapLMOF::VERBOSE ) cout << " IN EvapLMOF::ReadHeader()" << endl;
        
        // read magic numbers to determine if it is a .pax file and set the endian type.
        Int_t magic1 = 0, magic2 = 0;

        fStatus = kGood; 

        ::fread(&magic1,sizeof(Int_t),1,fInfile); ::fread(&magic2,sizeof(Int_t),1,fInfile);
        if ( ::ferror(fInfile) == 0 ) {
                ::rewind(fInfile); // rewind to start at the beginning of the file
                if ( magic1 == 0x00000028 && magic2 == 0x77777777 ) { /* cout << "LITTLE" << endl; */ SetRecord(Gws::Memory::kLittle); return; }
                if ( magic1 == 0x28000000 && magic2 == 0x77777777 ) { /* cout << "BIG" << endl; */ SetRecord(Gws::Memory::kBig); return; }      
        }       
        ::fclose(fInfile); fInfile = NULL; fStatus = kFail; // nothing good has been found, so status is kFail  

        if ( EvapLMOF::VERBOSE ) cout << " OUT EvapLMOF::ReadHeader()" << endl;
}

void EvapLMOF::NextFile()
{
        if ( EvapLMOF::VERBOSE ) cout << " IN EvapLMOF::NextFile()" << endl;
        
        if ( fMaxFiles == 0 ) { fStatus = kEmpty; return; }

        // in case the user wants to switch in a middle of a file
        if ( fInfile ) { if ( ::ferror(fInfile) ) cout << "WARNING: ferror detected in EvapLMOF::NextFile" << endl; ::fclose(fInfile); fInfile = NULL; } 
        
        if ( fCurrent < fMaxFiles ) {
                Int_t tmp = fCurrent;
                for (Int_t i = fCurrent+1; i < fMaxFiles; i++ ) { // look for a next good file

                        std::string fullname; 
                        fullname = fPath; fullname += fFileName[i]; fInfile = ::fopen(fullname.data(),"r");
                        
                        if ( fInfile ) { ReadHeader(); if (fStatus==kGood ) { tmp = i; break; } }
                } // i
                if ( tmp == fCurrent ) { fStatus = kEoE; } // fCurrent is the last file
                else { fCurrent = tmp; }
        }
        else { fStatus = kEoE; }
        
        if ( EvapLMOF::VERBOSE ) cout << " OUT EvapLMOF::NextFile()" << endl;
}

Bool_t EvapLMOF::ReadGlobal()
{
        Int_t id; 

        ::fread(fGlobalRecord->GetBuffer(),fGlobalRecord->Size(),1,fInfile); if ( ::ferror(fInfile) != 0 ) return kFALSE;
        
        // start decoding
        fGlobalRecord->SetOffset();  
        (*fGlobalRecord) >> id; if ( id != 0x77777777 ) { fError = " Wrong GLOBAL record id "; return kFALSE; } // check if it is a global record 
                
        (*fGlobalRecord) >> fGlobal.ZC;
        (*fGlobalRecord) >> fGlobal.AC;
        (*fGlobalRecord) >> fGlobal.NCASC;
        (*fGlobalRecord) >> fGlobal.MDIR;
        (*fGlobalRecord) >> fGlobal.JCASC;
        (*fGlobalRecord) >> fGlobal.ENER;
        (*fGlobalRecord) >> fGlobal.SIGMA;
        (*fGlobalRecord) >> fGlobal.VRC;
        (*fGlobalRecord) >> fGlobal.MAX_P;
        
        return ( !fGlobalRecord->IsStatus(Gws::Buffer::kFail) );
}

Bool_t EvapLMOF::ReadId()
{
        Short_t id, idok = 0xFFFF, idfission = 0xFFEE, idnodata = 0xFFFD; Char_t c; UChar_t uc;
                 
        ::fread(fIdRecord->GetBuffer(),fIdRecord->Size(),1,fInfile); if ( ::ferror(fInfile) != 0 ) return kFALSE;       
        
        // start decoding
        fIdRecord->SetOffset();  
        (*fIdRecord) >> id; 
        if ( !(id == idok || id == idfission || id == idnodata) ) { fError = " Wrong ID record "; return kFALSE; } // check if it is a valid id record
                
        if ( id == idfission ) fId.FIS = kTRUE; // the cascade ends up by fission
        else fId.FIS = kFALSE;
        
        fId.MG = 0; 
        
        (*fIdRecord) >> id; fId.EE = id;
        (*fIdRecord) >> uc; fId.JE = Short_t(uc);
        (*fIdRecord) >> uc; fId.L  = Short_t(uc); 
        (*fIdRecord) >> uc; fId.N  = Short_t(uc); 
        (*fIdRecord) >> c;  fId.Z =  Short_t(c);  if ( c < 0 ) { fId.Z = - Short_t(c); fId.FIS = kTRUE; } // -Z if fission (from evapor sources)
        (*fIdRecord) >> fId.WT; 
        (*fIdRecord) >> fId.NUM;
        
        if ( id == idnodata ) fId.NUM = 0; // no data record exists
        
        return ( !fIdRecord->IsStatus(Gws::Buffer::kFail) );
}

Bool_t EvapLMOF::ReadData()
{
        Bool_t decode;
        
        // z,n for the different particles
        const Short_t zp[MAXPARTICLES] = { 0,0,1,2,1,1,2,3,8,0,0,0  };
        const Short_t ap[MAXPARTICLES] = { 0,1,1,4,2,3,3,6,16,0,0,0 };  
        
        Int_t lengthstart, lengthend; Short_t id, dz, da; Char_t c;
        
        decode = kTRUE;
        // read the real buffer size 
        ::fread(fHeader->GetBuffer(),fHeader->Size(),1,fInfile); fHeader->SetOffset(); (*fHeader) >> lengthstart; 
        
        if ( lengthstart > (Int_t)fDataRecord->Size() || fId.NUM > MAXCASCADE ) 
                { fWarning = " Large cascade reached and skipped "; ::fseek(fInfile,lengthstart,SEEK_CUR); decode = kFALSE; }
        else 
                ::fread(fDataRecord->GetBuffer(),lengthstart,1,fInfile); // read all cascades in one operation

        // check if equal to the header part
        ::fread(fHeader->GetBuffer(),fHeader->Size(),1,fInfile); fHeader->SetOffset(); (*fHeader) >> lengthend;
        
        if ( ::ferror(fInfile) != 0 )    { fError = " PB to read data record from file "; return kFALSE; }
        if ( lengthstart != lengthend  ) { fError = " Start and end size of data record don't match "; return kFALSE; }
        
        if ( fIdOnly || decode == kFALSE ) return kTRUE; // in case the detail of the cascade is not needed
        
        // rewind buffer to start getting data
        fDataRecord->SetOffset(); 
        ::memset(&fData,0,sizeof(fData)); ::memset(&fKin,0,sizeof(fKin)); fKin.ZB = fId.Z; fKin.AB = fId.Z + fId.N; // starts from residual nucleus 
        
        // first loop to fill the cascade and compute some values
        for (Short_t i = 0; i < fId.NUM; i++ ) {
                (*fDataRecord) >> id; fData.EP[i] = id / 10.;
                (*fDataRecord) >> id; fData.EI[i] = id;
                
                (*fDataRecord) >> c; fData.MODE[i] = Short_t(c);
                (*fDataRecord) >> c;   fData.JI[i] = Short_t(c);
                (*fDataRecord) >> c;   fData.MP[i] = Short_t(c);
                (*fDataRecord) >> c;   fData.LP[i] = Short_t(c);
                
//              if ( fId.MODE[i] < 1 || fId.MODE[i] > MAXPARTICLES-1 )  return kFALSE; // unknown type
                if ( fData.MODE[i] <= fGlobal.MAX_P ) { // this is a charged particle
                        dz = zp[fData.MODE[i]];
                        da = ap[fData.MODE[i]];
                }
                else { // this is a gamma
                        dz = 0;
                        da = 0;         
                }
                
                fKin.ZB += dz; 
                fKin.AB += da; 
                
                if ( fData.MODE[i] == fGlobal.MAX_P + 4 ) { // E2 yrast cascade
                        fData.MUL[fGlobal.MAX_P+4] = fData.LP[i];
                        fData.LP[i] = 0 ;
                        fData.ES[fGlobal.MAX_P+4]  = fData.EI[i];
                }
                else {
                        fData.MUL[fData.MODE[i]] += 1 ;
                        fData.ES[fData.MODE[i]]  += fData.EP[i]; 
                }
                if ( fData.MODE[i] == fGlobal.MAX_P + 3 ) fId.FIS = kTRUE; // ends up by fission ..     
        }       
        // to set the gamma-rays multiplicity
        fId.MG = fData.MUL[fGlobal.MAX_P+1] + fData.MUL[fGlobal.MAX_P+2] + fData.MUL[fGlobal.MAX_P+4];
        
        // residual nucleus
        fKin.V0[0] = fGlobal.VRC * fGlobal.MDIR;
        fKin.V0[1] = 0.0;       
        fKin.V0[2] = fGlobal.VRC * (1.0 - fGlobal.MDIR);
        
        Int_t ac_run = fGlobal.AC, ar = ac_run, ml; Float_t vl[3], ct, avl, fug; 
        
        const Float_t massu = 931.5;
        
        for (Short_t i = 0; i < fId.NUM; i++ ) {
        
                Short_t jmode = fData.MODE[i]; // particle type
                
                if ( i < fId.NUM-1 ) {
                        fData.DJ[jmode] = fData.DJ[jmode] + fData.JI[i] - fData.JI[i+1];
                        fKin.JF[i] = fData.JI[i+1];
                        ml = fData.MP[i] - fData.MP[i+1];
                }
                else { // last data record
                        if( jmode <= fGlobal.MAX_P ) { // charged particle
                                fKin.JF[i] = fData.JI[i]; ml = 0;
                        }
                        else { 
                                if( jmode == (fGlobal.MAX_P + 1) ) { fKin.JF[i] = fData.JI[i]; ml = 0; }
                                else {
                                        if( jmode == (fGlobal.MAX_P + 2) ) {
                                                ml = 0;
                                                fKin.JF[i] = fData.JI[i]-2; fData.DJ[jmode] += 2;
                                        }
                                }
                        }
                }       
                Int_t ll = fData.LP[i], al; // Short_t mdr = mdir;
                al = 0;
                if ( jmode <= fGlobal.MAX_P ) al = ap[jmode] ;
                if( al > 0.0 ) {
                          getang_(&ll,&ml,&ct,&fKin.THETA[i],&fKin.PHI[i]);
                          
                          Double_t st = TMath::Sin(fKin.THETA[i]);
                          Double_t cp = TMath::Cos(fKin.PHI[i]);
                          Double_t sp = TMath::Sin(fKin.PHI[i]);
                          
                          ar = ac_run - al; avl = TMath::Sqrt(2.*fData.EP[i]/(massu*al)); fug= al / ar ;
                          
                          vl[2] =avl*ct; 
                          vl[1] =avl*st*sp; 
                          vl[0] =avl*st*cp;
                          
                          for (Int_t j = 0; j < 3; j++ ) {
                              fKin.VE[i][j] = fKin.V0[j] + vl[j];
                              fKin.V0[j]    = fKin.V0[j] - fug*vl[j];
                          }
                          ac_run = ar ;
                }
                else {
                          for (Int_t j = 0; j < 3; j++ ) {
                             fKin.VE[i][j] = 0.0;
                          }     
                } 
                fKin.EEJ[i] = 0.5 * al * massu * ( fKin.VE[i][0]*fKin.VE[i][0] + fKin.VE[i][1]*fKin.VE[i][1] + fKin.VE[i][2]*fKin.VE[i][2] );
        } // i 
        fKin.ER = 0.5 * ar * massu * ( fKin.V0[0]*fKin.V0[0] + fKin.V0[1]*fKin.V0[1] + fKin.V0[2]*fKin.V0[2] );
        
        return kTRUE;
}
        
void EvapLMOF::Rewind()
{
        if ( EvapLMOF::VERBOSE ) cout << " IN EvapLMOF::Rewind()" << endl;
        
        if ( fInfile ) { ::fclose(fInfile);  fInfile = NULL; }
        
        // change current to 0, reset statistics and load the next file
        fStatus = kEmpty;
        fWarning = "";
        fError = "";
        fCurrent = -1; for ( Int_t i = 0; i < fMaxFiles; i++ ) { fEventsRead[i] = 0.0; }
        
        NextFile();
        
        if ( EvapLMOF::VERBOSE ) cout << " OUT EvapLMOF::Rewind()" << endl;
}

Bool_t EvapLMOF::NextEvent(UInt_t /*i*/)
{
        Int_t bsize; Bool_t isok; 

        // does not read until status is good !
        // if status good, the header part has already been read and the first event is to be read
        if ( fInfile == NULL || ::ferror(fInfile) != 0 || fStatus != kGood ) {  NextFile(); if ( fStatus != kGood ) return kFALSE; }
                
        isok = kTRUE; ::fread(fHeader->GetBuffer(),fHeader->Size(),1,fInfile);
        
        if ( ::feof(fInfile) != 0 ) { // apparently eof is detected only when trying to read after the end
                NextFile(); 
                if ( fStatus != kGood ) return kFALSE; 
                else ::fread(fHeader->GetBuffer(),fHeader->Size(),1,fInfile);
        }        
        if ( ::ferror(fInfile) == 0 ) { 

                fHeader->SetOffset(); (*fHeader) >> bsize ;

//              cout << " Size " << bsize << endl;      
                if ( bsize == 40 ) { // this is a global record
//                      cout << " Block Record " << endl;               
                        if ( ReadGlobal() == kFALSE ) 
                                { /* cout << " GLOBAL record wrong " << endl; */ fStatus = kFail; return kFALSE; }
                        // a new run has been read, so next record must be identified
                        else { ::fread(fHeader->GetBuffer(),fHeader->Size(),1,fInfile); fHeader->SetOffset(); (*fHeader) >> bsize ; }
                }
                if ( bsize == 14  ) { // this is an id record
                        if ( ReadId() == kFALSE )   { /* cout << " ID record wrong " << endl; */ fStatus = kFail; return kFALSE; }
                        else {
                                if ( ReadData() == kFALSE ) { 
                                        /* cout << " Data record wrong " << endl; */ fStatus = kFail; return kFALSE; 
                                }
                                else fEventsRead[fCurrent] += 1.0; // a good event has just been read
                        }
                }
                else { fError = "Wrong size at position "; fError += (::ftell(fInfile)); fStatus = kFail; isok = kFALSE; }
        }
        else { fStatus = kFail; isok = kFALSE; }
        
        return isok;
}

void EvapLMOF::Scan(UInt_t start, UInt_t end)
{
        UInt_t stop;
        
        // to start form the beginning
        Rewind();
        
        // start looping on events
        stop = 0;
        for (UInt_t i = 1; i < end+1; i++ ) { 
                if ( NextEvent() == kFALSE ) { stop = i; break; }
                if ( i >= start+1 ) { cout << " ***** Event # " << i << endl; ShowCurrentEvent(); }
        }
        if ( stop != 0 ) {  
                cout << " !!!!! Breaking at Event # " << stop << endl; cout << fError.Data() <<  endl;
        }
}

TLeaf *EvapLMOF::GetLeaf(const char *leafname)
{
        TLeaf *leaf = NULL; TString tmp = leafname;
        
        // leaf corresponding to global record  
        if ( tmp == "AC" )     { leaf = new TLeafI(); leaf->SetAddress(&fGlobal.AC); }
        if ( tmp == "ZC" )     { leaf = new TLeafI(); leaf->SetAddress(&fGlobal.ZC); }
        if ( tmp == "NCASC" )  { leaf = new TLeafI(); leaf->SetAddress(&fGlobal.NCASC); }
        if ( tmp == "MDIR" )  { leaf = new TLeafI(); leaf->SetAddress(&fGlobal.MDIR);  }
        if ( tmp == "JCASC"  ) { leaf = new TLeafI(); leaf->SetAddress(&fGlobal.JCASC);  }
        if ( tmp == "ENER"  )  { leaf = new TLeafF(); leaf->SetAddress(&fGlobal.ENER);  }
        if ( tmp == "SIGMA" )  { leaf = new TLeafF(); leaf->SetAddress(&fGlobal.SIGMA);}
        if ( tmp == "VRC" )    { leaf = new TLeafF(); leaf->SetAddress(&fGlobal.VRC);}
        if ( tmp == "MAX_P" )  { leaf = new TLeafI(); leaf->SetAddress(&fGlobal.MAX_P);}
                        
        // leaf corresponding to ID record
        if ( tmp == "MG" )   { leaf = new TLeafS(); leaf->SetAddress(&fId.MG); }
        if ( tmp == "EE" )   { leaf = new TLeafF(); leaf->SetAddress(&fId.EE); }
        if ( tmp == "JE" )   { leaf = new TLeafF(); leaf->SetAddress(&fId.JE); }
        if ( tmp == "L"  )   { leaf = new TLeafF(); leaf->SetAddress(&fId.L); }
        if ( tmp == "Z" )   { leaf = new TLeafS(); leaf->SetAddress(&fId.Z); }
        if ( tmp == "N"  )   { leaf = new TLeafS(); leaf->SetAddress(&fId.N); }
        if ( tmp == "NUM" )  { leaf = new TLeafS(); leaf->SetAddress(&fId.NUM); }
        if ( tmp == "WT" )   { leaf = new TLeafF(); leaf->SetAddress(&fId.WT); }        
        if ( tmp == "FIS" )  { leaf = new TLeafB(); leaf->SetAddress(&fId.FIS);}
        
        // leaf corresponding to data record
        if ( tmp == "MODE" ) { leaf = new TLeafS(); leaf->SetAddress(fData.MODE); }
        if ( tmp == "EI" )   { leaf = new TLeafF(); leaf->SetAddress(fData.EI); }
        if ( tmp == "EP" )   { leaf = new TLeafF(); leaf->SetAddress(fData.EP); }
        if ( tmp == "JI"  )  { leaf = new TLeafS(); leaf->SetAddress(fData.JI); }
        if ( tmp == "LP"  )  { leaf = new TLeafS(); leaf->SetAddress(fData.LP); }
        if ( tmp == "MP"  )  { leaf = new TLeafS(); leaf->SetAddress(fData.MP); }
        if ( tmp == "MUL" )  { leaf = new TLeafS(); leaf->SetAddress(fData.MUL); }
        if ( tmp == "ES" )   { leaf = new TLeafF(); leaf->SetAddress(fData.ES); }       
        if ( tmp == "DJ" )   { leaf = new TLeafF(); leaf->SetAddress(fData.DJ);}        
        
        // leaf corresponding to kinematic record
        if ( tmp == "V0" )   { leaf = new TLeafF(); leaf->SetAddress(fKin.V0); }
        if ( tmp == "VE" )   { leaf = new TLeafF(); leaf->SetAddress(fKin.VE); }
        if ( tmp == "THETA" ){ leaf = new TLeafF(); leaf->SetAddress(fKin.THETA); }
        if ( tmp == "PHI"  ) { leaf = new TLeafF(); leaf->SetAddress(fKin.PHI); }
        if ( tmp == "ZB"  )  { leaf = new TLeafS(); leaf->SetAddress(&fKin.ZB); }
        if ( tmp == "AB"  )  { leaf = new TLeafS(); leaf->SetAddress(&fKin.AB); }
        if ( tmp == "ER" )   { leaf = new TLeafF(); leaf->SetAddress(&fKin.ER); }
        if ( tmp == "EEJ" )  { leaf = new TLeafF(); leaf->SetAddress(fKin.EEJ); }       
        if ( tmp == "JF" )   { leaf = new TLeafF(); leaf->SetAddress(fKin.JF);} 
                
        return leaf;
}

/*
void EvapLMOF::Map()
{
        Int_t bsize, nbglobal, nbid, nbdata, bytesread;
        
        // to start form the beginning
        nbglobal = nbid = nbdata = bytesread = 0; Rewind();
        
        // does not read until status is good !
        // if status good, the header part has already been read and the first event is to be read
        if ( fInfile == NULL || ::ferror(fInfile) != 0 || fStatus != kGood ) {  NextFile(); if ( fStatus != kGood ) { cout << " Nothing to be read " << endl; return;} }

        // to get the next id record
        bytesread += ::fread(fHeader->GetBuffer(),1,4,fInfile);  
        
        while ( ::ferror(fInfile) == 0 || ::feof(fInfile) == 0 ) { 

                fHeader->SetOffset(); (*fHeader) >> bsize ;

                if ( bsize == 40 ) { // this is a global record
                        if ( ReadGlobal() == kFALSE ) 
                                { cout << " GLOBAL record wrong " << endl; fStatus = kFail; break; }
                        else { bytesread += 44; nbglobal++; }
                        
                        // a new run has been read, so next record must be identified
                        bytesread += ::fread(fHeader->GetBuffer(),1,4,fInfile); fHeader->SetOffset(); (*fHeader) >> bsize ;     
                }
                if ( bsize == 14  ) { // this is an id record
                        if ( ReadId() == kFALSE )   { cout << " ID record wrong " << endl; fStatus = kFail; break; }
                        else { bytesread += 18; nbid++; }
                        
                        bytesread += ::fread(fDataRecord->GetBuffer(),1,fId.NUM*8+8,fInfile);
                        fEventsRead[fCurrent]+=1.0; // a good event has just been read
                }
                else { cout << "Wrong size at position " << (::ftell(fInfile)) << " 14 expected instead of " << bsize << endl; fStatus = kFail; break;}
                
                bytesread += ::fread(fHeader->GetBuffer(),1,4,fInfile); // next id
        } // 

        cout << " # global record " << nbglobal << ", # id record " << nbid << ", bytes read " << bytesread << endl; 
        ShowCurrentConditions();        
        
        // to re-start form the beginning
        Rewind();       
}
*/