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

#include "BaseGEM.h"

// Required for the wigner symbols / legendre in MathMore
#include "RConfigure.h"
#include "TFile.h"
#include "TMath.h"
#include "TLorentzVector.h"
#include "TF1.h"
#include "TROOT.h"

#include "Spin.h"
#include "GammaLink.h"
#include "NuclearLevel.h"

using namespace Gw;

#ifdef R__HAS_MATHMORE
#include "Math/SpecFuncMathMore.h"
//
namespace { // anonymous namespace for specifics to Angular Distribution. Hidden from the user point of view because of the root dependancy to wigner symbol
    
    Gw::Spin s_zero;

    // Population of one sub-state m with gaussian center at 0 and with a width sigma
    Double_t P(const Gw::Spin &Ii, const Gw::Spin &m, Double_t sigma)
    {
        Double_t temp = 0.0;
        
        for (Gw::Spin i = -Ii ; i <= Ii ; i+=1) {
            temp += TMath::Exp(-(i.Get()*i.Get())/(2*sigma*sigma));
            if ( gDebug > 0 )
                printf("temp = %f , sigma = %f , exp = %f\n",temp,sigma,TMath::Exp(-(i.Get()*i.Get())/(2*sigma*sigma)));
        }
        
        return TMath::Exp(-(m.Get()*m.Get())/(2*sigma*sigma))/temp;
    }
    Double_t ClebschGordan(const Gw::Spin &j1, const Gw::Spin &m1, const Gw::Spin &j2, const Gw::Spin &m2, const Gw::Spin &J, const Gw::Spin &M)
    {
        Double_t CG = TMath::Power(-1,j1.Get()-j2.Get()+M.Get())
                    * TMath::Sqrt(2*J.Get()+1)
                    * ROOT::Math::wigner_3j(j1.GetInTwoHalfUnit(),j2.GetInTwoHalfUnit(),J.GetInTwoHalfUnit(),m1.GetInTwoHalfUnit(),m2.GetInTwoHalfUnit(),-M.GetInTwoHalfUnit());
        //Double_t CG1 = TMath::Power(-1,j1.Get()-j2.Get()+M.Get());
        //Double_t CG2 = TMath::Sqrt(2*J.Get()+1);
        //Double_t CG3 = ROOT::Math::wigner_3j(j1.GetInTwoHalfUnit(),j2.GetInTwoHalfUnit(),J.GetInTwoHalfUnit(),m1.GetInTwoHalfUnit(),m2.GetInTwoHalfUnit(),-M.GetInTwoHalfUnit());
        
        //printf("Calculated Parameters CG1 %f CG2 %f CG3 %f\n",CG1,CG2,CG3);
        
        return CG;
    }
    Double_t B(const Gw::Spin &Ii, Int_t lambda)
    {
        Double_t result = 0, CG = 0;
        Gw::Spin s_lambda(lambda,1);
        
        if( !Ii.IsHalfInt() ) // I1 integer
        {
            CG = ClebschGordan(Ii,s_zero,Ii,s_zero,s_lambda,s_zero);
            result = TMath::Sqrt(2*Ii.Get()+1)*TMath::Power(-1,Ii.Get())*CG;
        }
        else  // Ii half integer
        {
            Gw::Spin s_1_2(1,2), m_s_1_2(-1,2);
            
            CG = ClebschGordan(Ii,m_s_1_2,Ii,s_1_2,s_lambda,s_zero);
            result = TMath::Sqrt(2*Ii.Get()+1)*(-1)*TMath::Power(-1,Ii.Get()+0.5)*CG;
        }
        //printf("Calculated Parameters CG %f result %f\n",CG,result);
        
        return result;
    }
    Double_t B(const Gw::Spin &Ii, Double_t lambda, Double_t sigma)
    {
        
        Double_t result = 0.0, CG = 0.0; Gw::Spin s_lambda(lambda,1);
        
        for (Gw::Spin m = -Ii ; m <= Ii ; m+=1)
        {
            Gw::Spin m_m = -m;
            CG = ClebschGordan(Ii,m_m,Ii,m,s_lambda,s_zero);
            
            result += TMath::Power(-1,Ii.Get()+m.Get())*CG*P(Ii,m,sigma);
        }
        result = TMath::Sqrt(2*Ii.Get()+1)*result;
        
        return result;
    }
    // F-coefficients ordinaires dans le cas d'une transition avec melange
    Double_t F(const Gw::Spin &Ii, const Gw::Spin &If, Int_t L1, Int_t L2, Double_t lambda)
    {
        Double_t threeJ, sixJ, result;
        
        threeJ = ROOT::Math::wigner_3j(2*L1,2*L2,2*lambda,2,-2,0);
        sixJ   = ROOT::Math::wigner_6j(2*L1,2*L2,2*lambda,Ii.GetInTwoHalfUnit(),Ii.GetInTwoHalfUnit(),If.GetInTwoHalfUnit());
        
        result =  TMath::Power(-1,1+Ii.Get()+If.Get())*TMath::Sqrt((2*lambda+1)*(2*L1+1)*(2*L2+1)*(2*Ii.Get()+1))*threeJ*sixJ;
        
        //printf("Calculated Parameters three %f six %f result %f\n",threeJ,sixJ,result);
        
        return result;
    }
    // ROOT suitable function to draw angular distribution up to order 4.
    Double_t Wth(Double_t *x, Double_t *par)
    {
        Double_t p2,p4;
        p2 = ROOT::Math::legendre(2,TMath::Cos(x[0]));
        p4 = ROOT::Math::legendre(4,TMath::Cos(x[0]));
        
        return par[0]+par[1]*p2+par[2]*p4;
    }
    Double_t WthDeg(Double_t *x, Double_t *par)
    {
        Double_t p2,p4;
        p2 = ROOT::Math::legendre(2,TMath::Cos(x[0]*TMath::Pi()/180.));
        p4 = ROOT::Math::legendre(4,TMath::Cos(x[0]*TMath::Pi()/180.));
        
        return par[0]+par[1]*p2+par[2]*p4;
    }
    // ROOT suitable function to draw angular distribution up to order 4
    Double_t WthCos(Double_t *x, Double_t *par)
    {
        Double_t p2,p4;
        p2 = ROOT::Math::legendre(2,x[0]);
        p4 = ROOT::Math::legendre(4,x[0]);
        
        return par[0]+par[1]*p2+par[2]*p4;
    }
}
#endif

ClassImp(BaseGEM);

BaseGEM::BaseGEM():
    TObject(),
    fLevelScheme(0x0), fFeeding(), fRandFeeding(), fRandDown(), fAD(), fRand0(), fTmp(), fEoC(),
    fDirection(kDown),
    fMinEnergyFactor(10.0),
    fADType(0),
    fCutLifeTime(0.0000000001), /* 0.1*ns from Krane 1971 */
    fSigma(0.0)
{
        fTmp.SetOwner(kTRUE);
    fRandDown.SetOwner(kTRUE);
    fFeeding.SetOwner(kTRUE);
    fRandFeeding.SetOwner(kTRUE);
    fAD.SetOwner(kTRUE);
}

BaseGEM::~BaseGEM()
{
        BaseGEM::Clear(""); // clear everything
    
    if ( fLevelScheme ) {
        delete fLevelScheme;
        fLevelScheme = 0x0;
    }
}

void BaseGEM::Clear(Option_t *o)
{
        TString opt = o;
        if ( opt.Contains("all") ) {
        delete fLevelScheme; fLevelScheme = 0x0;
    }
        
        fRand0.Clear(); fFeeding.Clear(); fRandFeeding.Clear(); fRandDown.Clear(); fTmp.Clear(); fAD.Clear();
    fMinEnergyFactor = 10.0;
    fDirection = kDown;
}

Int_t BaseGEM::Import(const Char_t *filename, Option_t *opt)
{
    Int_t ok = 0; Gw::LevelScheme *newlev = 0x0;
    
        // init level scheme and test if initialisation has succeded
    TString lname(filename), lopt(opt);
    if ( lname.Contains(".root")) { // different way but could be integrated in lev->Import through the ReadTObject method
        TFile lroot(filename);
        if ( lroot.IsOpen() ) {
            newlev = (Gw::LevelScheme *)lroot.Get(opt);
            if ( newlev )
                ok = 1;
            else
                printf("Impossible to find levelscheme %s in %s \n",filename,opt);
        }
        else
            return ok;
    }
    else {
        newlev = new Gw::LevelScheme();
        ok = newlev->Import(filename,opt);
    }
    
    if ( newlev && ok ) {
        if ( fLevelScheme ) {
            delete fLevelScheme;
        }
        fLevelScheme = newlev;
        InitRandom();
    }
    return ok;
}

Bool_t BaseGEM::SetParameter(const char *name, Int_t value)
{
    TString tname(name);
    if ( tname == "ADType" ) {
        fADType = value;
        return true;
    }
    return false;
}

Bool_t BaseGEM::SetParameter(const char *name, Double_t value)
{
    TString tname(name);
    if ( tname == "MinEnergyFactor" ) {
        fMinEnergyFactor = value;
        return true;
    }
    if ( tname == "CutLifeTime" ) {
        fCutLifeTime = value;
        return true;
    }
    return false;
}

void BaseGEM::InitRandom()
{
        if ( fLevelScheme->GetLinks().GetSize() == 0 ) return; // nothing to be done
        
        Int_t i, iend, j, jend; Link *gam1 = NULL, *gam2 = NULL;
        
        // clear the previous definition BUT NOT the level scheme which has just been read.
        BaseGEM::Clear("");
    
        // set each collection with the right number.
        fRandDown.Expand(fLevelScheme->GetLinks().GetSize());
        fRandFeeding.Expand(fLevelScheme->GetLinks().GetSize()); fFeeding.Expand(fLevelScheme->GetLinks().GetSize()); fAD.Expand(fLevelScheme->GetLinks().GetSize());
    
        // First it loops to check if all intensities are > 0 and it calculates the lowest positive one MIN .
        // If at least one intensity is found <=0,
        // all the gamma-rays with an intensity <= 0 are modified with an intensity = MIN / 100
        Float_t min = 1000000; Int_t nbwrong = 0;
        iend = jend = fLevelScheme->GetLinks().GetSize();
        for (i = 0; i < iend; i++) {
                gam1 = (Link *)fLevelScheme->GetLinks().At(i);
                if ( gam1->GetStrength().Get() <= 0 ) {
                        nbwrong++;
                }
                else { // more than 0
                        if ( gam1->GetStrength().Get() < min )
                                min = gam1->GetStrength().Get() ;
                }
        } // iend
        if ( nbwrong > 0 ) {
                printf("\n WARNING in BaseGEM::Import \n");
                printf("  %d links have been found with an intensity = 0 in %s \n",nbwrong,fLevelScheme->GetName());
                printf("  The intensity of those links are set to the weakest intensity %f / %f \n \n",min,fMinEnergyFactor);
                for (i = 0; i < iend; i++) {
                        gam1 = (Link *)fLevelScheme->GetLinks().At(i);
                        if ( gam1->GetStrength().Get() <= 0 ) {
                                gam1->GetStrength().Set(min/fMinEnergyFactor);
                        } // if
                } // i
        } // nbwrong
        
        // now for each gamma, a generator is created to go down in the path and to set fRand0 which determine the entry point
        RandObj *down, *downfeeding; TObject *obj;
        Float_t sum_int_up_in, sum_int_up_out, sum_int_down_in, sum_int_down_out;
        
        for (i = 0; i < iend; i++) { // loop on all gammas to build the cascade generator
        
                // current link with its associated generators
                gam1 = (Link *)fLevelScheme->GetLinks().At(i); Gw::GammaLink *gamma = (Gw::GammaLink *)gam1;
        //
        down = new RandObj();
        downfeeding = new RandObj();
                
                // feeding link for gam1
                Link *feeding = new Link();
                feeding->SetIL(gam1->GetIL());
        feeding->SetFL(gam1->GetIL());
                
                obj = new TObject(); // gam1 is one of the link directly in coincidence with gam1
                fTmp.Add(obj);
                obj->SetUniqueID(i);
                downfeeding->Add(obj,gam1->GetStrength().Get());
        
                // look for all gammas up and down in coincidence with gam1
                sum_int_up_in = 0;
                sum_int_up_out  = gam1->GetStrength().Get();
                sum_int_down_in = gam1->GetStrength().Get();
                sum_int_down_out = 0;
        
                for (j = 0; j < jend; j++) { // look for gammas in coincidence with gamma # i
            
                        if ( i == j ) continue;
            
                        gam2 = (Link *)fLevelScheme->GetLinks().At(j);
            
                        if ( gam1->GetIL() == gam2->GetFL() ) { // gammas filling the initial level of gam1
                                sum_int_up_in += gam2->GetStrength().Get();
                        }
                        if ( gam1->GetIL() == gam2->GetIL() ) { // gammas desexciting the initial level of gam1
                                sum_int_up_out += gam2->GetStrength().Get();
                                
                                obj = new TObject();
                                fTmp.Add(obj);
                                obj->SetUniqueID(j);
                                downfeeding->Add(obj,gam2->GetStrength().Get());
                        }
                        if ( gam1->GetFL() == gam2->GetIL() ) { // gammas desexciting the final level of gam1
                                sum_int_down_out += gam2->GetStrength().Get();
                                
                                obj = new TObject();
                                fTmp.Add(obj);
                                obj->SetUniqueID(j);
                                down->Add(obj,gam2->GetStrength().Get());
                        }
            
                        if ( gam1->GetFL() == gam2->GetFL() ) { // gammas filling the final level of gam1
                                sum_int_down_in += gam2->GetStrength().Get();
                        }
                } // jend
                
                // lateral feeding treatment
        // 01-2015 because in the procedurethe way it is done is ti select randomly directly one gamma of the level scheme
        // but through the feeding part, one should give a fraction of the feeinf to all gamma de exciting the level
                if ( sum_int_up_out > sum_int_up_in ) { // the feeding for gam1 has a intensity greater than 0
                        feeding->GetStrength().Set( (sum_int_up_out - sum_int_up_in ) * gam1->GetStrength().Get() / (sum_int_up_out) );
                }
                else { feeding->GetStrength().Set(0.0); }
        
                // isomeric state not connected and end of cascade treatment
                if ( sum_int_down_in > sum_int_down_out ) { down->Add(&fEoC,sum_int_down_in - sum_int_down_out); }
        
                fRand0.Add(feeding,feeding->GetStrength().Get()); fFeeding.AddAt(feeding,i);
                fRandFeeding.AddAt(downfeeding,i);
                fRandDown.AddAt(down,i);
        
        // ADType 0 means no angular distribution i.e. use gRandom->Sphere
        fAD.AddAt(0x0,i);
        if ( fADType == 1 && gam1->InheritsFrom("Gw::GammaLink") && gamma->GetLambda() > 0 ) {

#ifdef R__HAS_MATHMORE
            Int_t lambda_max = 0, L, L1, L2; Double_t a0,a2,a4;  // name and parameters of the function to be used ... limited to order 4
            a0 = 0.0;
            a2 = 0.0;
            a4 = 0.0;
            Gw::Spin Ii, If;
            
            if ( fSigma <= 0.0 ) { // pure orientation of sub-magnetic states
                //
                Ii = ((Gw::NuclearLevel *)gamma->GetIL())->GetSpin();
                If = ((Gw::NuclearLevel *)gamma->GetFL())->GetSpin();
                Double_t delta = gamma->GetMixing().Get();
                
                if ( delta == 0.0 ) { // pure transition
                    L = gamma->GetLambda(); lambda_max = 2*gamma->GetLambda();
                    a0 = B(Ii,0)*F(Ii,If,L,L,0);
                    if ( lambda_max >=2  )
                        a2 = B(Ii,2)*F(Ii,If,L,L,2);
                    if ( lambda_max >=4  )
                        a4 = B(Ii,4)*F(Ii,If,L,L,4);
                }
                else { // mixed transition L, L+1
                    L = gamma->GetLambda(); lambda_max = 2*gamma->GetLambda() + 1, L1 = L, L2 = L+1;
                    // cmompute parameters
                    //if ( L1 == L2 ) { // stay null in case L1 != L2 because |L1-L2| <= labda <= L1+L2
                    //    a0 = B(Ii,0)*F(Ii,If,L1,L2,0);
                    //}
                    if (lambda_max >=2 )
                        a2 = B(Ii,2)*( F(Ii,If,L1,L1,2) + 2*delta*F(Ii,If,L1,L2,2) + delta*delta*F(Ii,If,L2,L2,2)  ) / (1+delta*delta);
                    if (lambda_max >=4 )
                        a4 = B(Ii,4)*( F(Ii,If,L1,L1,4) + 2*delta*F(Ii,If,L1,L2,4) + delta*delta*F(Ii,If,L2,L2,4)  ) / (1+delta*delta);
                }
            }
            else {
                //
                Ii = ((Gw::NuclearLevel *)gamma->GetIL())->GetSpin();
                If = ((Gw::NuclearLevel *)gamma->GetFL())->GetSpin();
                Double_t delta = gamma->GetMixing().Get();
                
                if ( delta == 0.0 ) { // pure transition
                    L = gamma->GetLambda(); lambda_max = 2*gamma->GetLambda();
                    //
                    a0 = B(Ii,0,fSigma)*F(Ii,If,L,L,0);
                    if (lambda_max >=2 )
                        a2 = B(Ii,2,fSigma)*F(Ii,If,L,L,2);
                    if (lambda_max >=4 )
                        a4 = B(Ii,4,fSigma)*F(Ii,If,L,L,4);
                }
                else { // mixed transition L, L+1
                    L = gamma->GetLambda(); lambda_max = 2*gamma->GetLambda() + 1, L1 = L, L2 = L+1;
                    
                    if ( L1 == L2 ) { // stay null in case L1 != L2 because |L1-L2| <= labda <= L1+L2
                        //     a0 = B(Ii,0.fSigma)*F(Ii,If,L1,L2,0);
                    }
                    if (lambda_max >=2 )
                        a2 = B(Ii,2,fSigma)*( F(Ii,If,L1,L1,2) + 2*delta*F(Ii,If,L1,L2,2) + delta*delta*F(Ii,If,L2,L2,2)  ) / (1+delta*delta);
                    if (lambda_max >=4 )
                        a4 = B(Ii,4,fSigma)*( F(Ii,If,L1,L1,4) + 2*delta*F(Ii,If,L1,L2,4) + delta*delta*F(Ii,If,L2,L2,4)  ) / (1+delta*delta);
                    
                }
            }
            
            if ( gamma->GetTau().Get() < fCutLifeTime ) { // only if live time of initial state not long i.e. orientation is not washed out by electric/magnetic fields around
                TString name_function = Form("%s_link_%.1f_%.1f_%.1f",fLevelScheme->GetName(),gamma->GetEnergy().Get(),Ii.Get(),If.Get());
                TF1 *W = new TF1(name_function.Data(),Wth,0,TMath::Pi(),3);
                W->SetParameter(0,1);
                W->SetParameter(1,a2/a0);
                W->SetParameter(2,a4/a0);
                W->SetNpx(1000);
                fAD.AddAt(W,i);
                
                // remove from list of function but add a clone that is not used by this. Note that it does not work with the root Clone method ? ...
                // thus the first function is just removed from list of function while a similar one is kept with just a similar name
                gROOT->GetListOfFunctions()->Remove(W);
                
                W = new TF1(name_function.Data(),Wth,0,TMath::Pi(),3);
                W->SetParameter(0,1);
                W->SetParameter(1,a2/a0);
                W->SetParameter(2,a4/a0);
                W->SetNpx(1000);
            }
#else
        printf(" *** Error, cannot set Angular Distribution in BaseGEM, ROOT MathMore lib probably not installed \n");
#endif
        }
                if ( gDebug > 0 ) {
            printf(" There are %d gammas feeding the gamma # %d \n",downfeeding->GetSize(),i);
                        printf(" There are %d gammas de-exciting the gamma # %d \n",down->GetSize(),i);
        }
        } // iend
}

/*
 void BaseGEM::ls(Option_t *o) const
 {
 Link *gam1 = NULL; RandObj *robj;  Int_t iend = fRand0.GetSize();
 
 fRand0.ls(o);
 for (Int_t i = 0; i < iend; i++ ) {
 gam1 = (GammaLink *)fRand0.At(i);
 if ( gam1 ) {
 printf("%f keV [%d] is stored with weight %f \n",gam1->GetEnergy().Get(),i,fRand0.WeightAt(i)) ;
 
 robj = (RandObj *)fRandUp.At(i); printf(" * %d gammas feed it \n",robj->GetSize()) ;
 for (Int_t j = 0; j < robj->GetSize(); j++) {
 if ( robj->At(j) == &fEoC ) printf("    Lateral feeding with weight %f \n",robj->WeightAt(j));
 else {
 printf("    Gamma #%d, feeding with weight %f \n",robj->At(j)->GetUniqueID(),robj->WeightAt(j)) ;
 }
 } // j
 robj = (RandObj *)fRandDown.At(i); printf(" * %d gammas desexcites it \n",robj->GetSize()) ;
 for (Int_t j = 0; j < robj->GetSize(); j++) {
 if ( robj->At(j) == &fEoC ) printf("    EoC with weight %f \n",robj->WeightAt(j));
 else {
 printf("    Gamma #%d, desexcitation with weight %f \n",robj->At(j)->GetUniqueID(),robj->WeightAt(j)) ;
 }
 } // j
 } // if gam1
 } // iend
 }
 */

// void BaseGEM::FillRandom(TH1 *h, Int_t ntimes) { fRand0.FillRandom(h,ntimes); }

Int_t BaseGEM::DoCascade(TSeqCollection &cas, Int_t which_first, Option_t *o)
{
        TObject *obj; Link *link; Int_t result = 0;
    
    if ( fRandFeeding.At( which_first ) == 0x0 )
        return result;
    
    TString opt(o);
    if ( !opt.Contains("add",TString::kIgnoreCase) )
        cas.Clear("nodelete");
    
        // Ask the random generator for that feeding for the next link in the level scheme and keep on going with it
        obj = ((RandObj *)fRandFeeding.At( which_first ))->Rand();
        while ( obj != &fEoC ) {
                link = (Link *)fLevelScheme->GetLinks().At(obj->GetUniqueID());
                result += link->DoCascade(cas,o);
                obj = ((RandObj *)fRandDown.At( obj->GetUniqueID() ))->Rand();
        }
        
        if ( fDirection == kUp ) { // if kUp, inverse the cascade
                TObjArray tmp(100);
                tmp.SetOwner(kFALSE); tmp.Clear("nodelete"); tmp.AddAll(&cas); cas.Clear("nodelete");
                
                for (Int_t i = tmp.GetEntries()-1; i >=0; i-- )
            cas.Add(tmp.At(i));
        }
    
    return result;
}

Int_t BaseGEM::DoCascade(TSeqCollection &cas, Option_t *o)
{
    Int_t result = 0;
    
        // select the 'entry' point then call a random cascade from it
        Int_t which_first = 0;
    fRand0.Rand(which_first);
    //
    result += DoCascade(cas,which_first,o);
    
    return result;
}

void BaseGEM::DoAngularDistribution(Int_t which_gamma, TLorentzVector &vector, Bool_t forceiso)
{
    TObject *generator = fAD.At(which_gamma);
    if ( generator == 0x0 || forceiso ) { // means not a function thus isotropic through gRandom
        Double_t x,y,z;
        Gw::Random::Instance()->Current()->Sphere(x,y,z,1.0);
        vector.SetXYZT(x,y,z,0);
        
        // printf("Isotrope for gamma # %d \n",which_gamma);
    }
    else {
        // in two steps. First in 2D through the AD get the theta angle then phi uniform distribution
        TVector3 v; Double_t theta = 0, phi = 0;
        
        TF1 *angdistri2d = (TF1 *)generator;
        theta = angdistri2d->GetRandom();
        phi = Gw::Random::Instance()->Current()->Uniform(0.,TMath::TwoPi());
        v.SetMagThetaPhi(1,theta,phi);
        
        vector.SetVect(v); vector.SetT(0.0);
        
        //printf("Angular Distribution for gamma # %d \n",which_gamma);
    }
}

Int_t BaseGEM::DoCascade(TSeqCollection &cas, TSeqCollection &directions, Int_t which_first, Option_t *o)
{
        TObject *obj; Link *link; Int_t result = 0, result_one_link = 0;

    if ( fRandFeeding.At( which_first ) == 0x0 )
        return result;
    
    TString opt(o);
    if ( !opt.Contains("add",TString::kIgnoreCase) ) {
        cas.Clear("nodelete");
        directions.Clear("nodelete");
    }
    
        // Ask the random generator for that feeding for the next link in the level scheme and keep on going with it
        obj = ((RandObj *)fRandFeeding.At( which_first ))->Rand();
        while ( obj != &fEoC ) {
                link = (Link *)fLevelScheme->GetLinks().At(obj->GetUniqueID());
        result_one_link += link->DoCascade(cas,o);
        
        for (Int_t i = result; i < result + result_one_link; i++) { // make sur the direction collection is synchronize with cas
            if ( directions.At(i) == 0 ) {
                directions.AddAt(new TLorentzVector(),i);
            }
        }
                if ( result_one_link == 1 ) {
            TLorentzVector *lv = (TLorentzVector *)directions.At(result);
            if ( cas.At(result) == link ) { // it means the link itself has been added to the cascade --> AD could be done otherwise X conversion ==> isotropic
                DoAngularDistribution(obj->GetUniqueID(),*lv);
            }
            else
                DoAngularDistribution(obj->GetUniqueID(),*lv,true);
        }
        else {
            for (Int_t i = 0; i < result_one_link; i++) {
                TLorentzVector *lv = (TLorentzVector *)directions.At(result+i);
                DoAngularDistribution(obj->GetUniqueID(),*lv,true);
            }
            printf("*** TO BE IMPLEMENTED in BaseGEM::DoCascade TSeqCollection &cas, TSeqCollection &directions, Int_t which_first, Option_t *o *** /n");
            printf("*** Isotropic has been forced for %d links produced by link # %d *** /n",result_one_link,obj->GetUniqueID());
       }
        result += result_one_link;
                obj = ((RandObj *)fRandDown.At( obj->GetUniqueID() ))->Rand();
        }
        
        if ( fDirection == kUp ) { // if kUp, inverse the cascade
                TObjArray tmpg(30), tmpd(30);
                tmpg.SetOwner(kFALSE); tmpg.Clear("nodelete"); tmpg.AddAll(&cas); cas.Clear("nodelete");
        tmpd.SetOwner(kFALSE); tmpd.Clear("nodelete"); tmpd.AddAll(&directions); directions.Clear("nodelete");

                for (Int_t i = tmpg.GetEntries()-1; i >=0; i-- ) {
            cas.Add(tmpg.At(i));
            directions.Add(tmpd.At(i));
        }
        }
    
    return result;
}

Int_t BaseGEM::DoCascade(TSeqCollection &cas, TSeqCollection &directions, Option_t *o)
{
    TString opt = o;
    
        // Reset the previous cascade.
        if ( !opt.Contains("add",TString::kIgnoreCase) ) {
        cas.Clear("nodelete");
        directions.Clear("nodelete");
    }
    
        // select the 'entry' point then call a random cascade from it
        Int_t which_first = 0;
    fRand0.Rand(which_first);
    //
    return DoCascade(cas,directions,which_first,o);
}