firstdemo.C

Go to the documentation of this file.

#ifndef __CINT__
#include <Link.h>
#include <XGammaLink.h>

#include <TStopwatch.h>
#include <TH1.h>
#include <TTree.h>

#include <iostream>
using namespace std;
#endif

using namespace Gw;

void CheckLS(const Char_t *name = "197Pb.ags")
{
        LevelScheme *lev = new LevelScheme;     
        if ( lev->InitAGS(name) != 0 ) lev->Draw();
}

void CheckX(const Char_t *name = "simu2.ags")
{
        GammaLink *link;
        
        LevelScheme lev(name,"ags");
        if ( lev.GetLinks().GetSize() == 0 )  { printf("No links in %s \n",name); }
        
        TIter iter(&lev.GetLinks()); // start an iterator on the list of links
        
        Int_t count = 0;
        while ( link = (GammaLink *)iter() ) { // next link
                if ( link->InheritsFrom("XGammaLink") ) count++; 
        }
        printf("%d x-links found in %s \n",count,name);
}

void CGate(TH1 *spgated, Cascade *cas, Int_t g1, Int_t g2)
{
        Int_t i, iend, nb_in, flag[200];
        
        iend = cas->GetSize(); nb_in = 0;
        for (i = 0; i < iend; i++ ) { // count how many gammas in the gate [g1,g2]

                flag[i] = 0;

                gam = (GammaLink *)cas->At(i);
                if ( gam->GetEnergy().Get() > g1 && gam->GetEnergy().Get() < g2 ) { nb_in++; flag[i] = 1; }
        } // i
        if ( nb_in == 0 ) return;       // no gammas in the gate
        if ( nb_in == 1 ) {     // keep only gammas outside the gate
          for (i = 0; i < iend; i++ ) {
                gam = (GammaLink *)cas->At(i);
                if ( flag[i] == 0 ) spgated->Fill(gam->GetEnergy().Get());
          } // i
        } else {                        // keep all gammas
          for (i = 0; i < iend; i++ ) {
                gam = (GammaLink *)cas->At(i); spgated->Fill(gam->GetEnergy().Get());
          } // i
        }
}

void simu1(Int_t nbcas = 100)
{
        // name of the gls file and definition of the gate
        const Char_t nameAGS[] = "simu2.ags";
        Float_t gl1 = 572, gr1 = 574, gl2 = 656, gr2 = 657;

        // load AGS file
        BaseGEM random;
        if ( random.InitAGS(nameAGS) == 0 ) {
                cout << " GLS level Scheme not correctly loaded " << endl; return;
        }

        // some spectra to be filled
        TH1F *K = new TH1F("K","K",100,0,100);
        TH1F *h = new TH1F("ProjTot","ProjTot",8192,0,4096);
        TH1F *hg1 = new TH1F("Gated1","Gated1",8192,0,4096); TH1F *hg2 = new TH1F("Gated2","Gated2",8192,0,4096);

        // start the Monte-Carlo
        Cascade cas; GammaLink *gam;
        for (Int_t i = 0; i < nbcas; i++ ) {

                random.DoCascade(cas,""); // to get the next cascade

                K->Fill(cas.GetSize()) ;
                for (Int_t j = 0; j < cas.GetSize(); j++ ) { // to set the total projection
                        gam = (GammaLink *)cas.At(j); h->Fill(gam->GetEnergy().Get());
                }
                CGate(hg1,&cas,gl1,gr1); CGate(hg2,&cas,gl2,gr2);
        }
        hg1->Draw();
}

void oneLS(const Char_t *nameAGS, Int_t nbcas = 100)
{
        // load the AGS file
        BaseGEM random;
        if ( random.InitAGS(nameAGS) == 0 ) {
                cout << " AGS level Scheme not correctly loaded " << endl; return;
        }

        // some spectra to be filled
        TH1F *K = new TH1F("K","K",100,0,100);
        TH1F *h = new TH1F("ProjTot","ProjTot",8192,0,4096);

        // start the Monte-Carlo
        Cascade cas; GammaLink *gam;
        for (Int_t i = 0; i < nbcas; i++ ) {

                random.DoCascade(&cas,"");

                K->Fill(cas.GetSize()) ;
                for (Int_t j = 0; j < cas.GetSize(); j++ ) {
                        gam = (GammaLink *)cas.At(j); h->Fill(gam->GetEnergy().Get());
                } // j
        } //i
        h->Draw();
}

void twoLS(Int_t nbcas = 100)
{
        // load the two AGS file
        BaseGEM pb197, pb198;
        if ( pb197.InitAGS("197Pb.ags") == 0 ) {
                cout << " 197Pb.ags level Scheme not correctly loaded " << endl; return;
        }
        if ( pb198.InitAGS("198Pb.ags") == 0 ) {
                cout << " 198Pb.ags level Scheme not correctly loaded " << endl; return;
        }
        // create the RandObj to select between 197Pb 90 % & 198Pb 10 %
        RandObj channel;
        channel.Add(&pb197,0.9);channel.Add(&pb198,0.1); 

        // some spectra to be filled
        TH1F *hpb197 = new TH1F("Pb197","Pb197",8192,0,4096); hpb197->SetLineColor(2);
        TH1F *hpb198 = new TH1F("Pb198","Pb198",8192,0,4096); hpb198->SetLineColor(4);

        // start the Monte-Carlo
        Cascade cas; BaseGEM *gem; GammaLink *gam;
        for (Int_t i = 0; i < nbcas; i++ ) {

                gem = (BaseGEM *)channel.Rand();        // select one nucleus
                gem->DoCascade(&cas);                   // determine a cascade on this nucleus

                if ( gem == &pb197 ) {
                        for (Int_t j = 0; j < cas.GetSize(); j++ ) {
                                gam = (GammaLink *)cas.At(j); hpb197->Fill(gam->GetEnergy().Get());
                        } // j
                } else {
                        for (Int_t j = 0; j < cas.GetSize(); j++ ) {
                                gam = (GammaLink *)cas.At(j); hpb198->Fill(gam->GetEnergy().Get());
                        } // j
                } // gem
        } //i
        hpb197->Draw(); hpb198->Draw("same");
}

void showIE(Int_t nbcas = 100)
{
        // name of the gls file and definition of the gate
        const Char_t nameAGS[] = "197Pb.ags";
        Float_t spin_min = 0, spin_max = 50, e_min = 0, e_max = 10000;

        // load AGS file
        BaseGEM random;
        if ( random.InitAGS(nameAGS) == 0 ) {
                cout << nameAGS << " level Scheme not correctly loaded " << endl; return;
        }

        // to display the entry point 
        TH2F *bidim = new TH2F("I*E","I*E",50,spin_min,spin_max,100,e_min,e_max);

        // start the Monte-Carlo
        Cascade li_gamma; 
        GammaLink *gam; NuclearLevel *levi;

        Float_t max_e, spin; 
        for (Int_t i = 0; i < nbcas; i++ ) {

                random.DoCascade(&li_gamma);

                max_e = 0;
                for (Int_t j = 0; j < li_gamma.GetSize(); j++ ) { // loop on links to fin the highest level
                        gam = (GammaLink *)li_gamma.At(j);
                        levi = (NuclearLevel *)gam->GetIL(); 
                        if ( levi->GetEnergy().Get() > max_e && levi->GetSpin().Get() > 0 ) { 
                                max_e = levi->GetEnergy().Get();
                                spin  = levi->GetSpin().Get();
                        }
                }
                bidim->Fill(spin,max_e);
        }
        bidim->Draw();
}

void simu1BG(Int_t nbcas = 100000)
{
        // name of the gls file and definition of the gate
        const Char_t nameAGS[] = "simu1.ags";

        // load AGS file
        BaseGEM random;
        if ( random.InitAGS(nameAGS) == 0 ) {
                cout << " GLS level Scheme not correctly loaded " << endl; return;
        }
        // create the function for the background
        TF1 bg("background","100-100*x/4096",0,4096);

        // create the RandObj to select between sim1.ags and bg
        RandObj channel;
        channel.Add(&random,0.001);channel.Add(&bg,0.999); 

        // some spectra to be filled
        TH1F *h = new TH1F("ProjTot","ProjTot",8192,0,4096);

        // start the Monte-Carlo
        TObject *obj; BaseGEM *gem; GammaLink *gam;

        Cascade cas; Int_t which_channel;
        for (Int_t i = 0; i < nbcas; i++ ) {

                obj = channel.Rand(which_channel); // channel selection

                if ( which_channel == 0 ) { // level scheme

                        random.DoCascade(&cas); // to get the next cascade

                        for (Int_t j = 0; j < cas.GetSize(); j++ ) { // to set the total projection
                                gam = (GammaLink *)cas.At(j); h->Fill(gam->GetEnergy().Get());
                        } // j
                } 
                else { // background
                        h->Fill(bg.GetRandom(1,4096));
                }
        }
        h->Draw();
}

void simu1X(Int_t nbcas = 100000)
{
        TObject *obj; XGammaLink *xlink;
        
        // name of the gls file and definition of the gate
        const Char_t nameAGS[] = "simu2.ags";
        const Int_t nbx = 2; 
        Float_t xe[nbx] = { 80, 82 }; Float_t xi[nbx] = { 80, 20 };
        
        // load the AGS file
        BaseGEM random;
        if ( random.InitAGS(nameAGS) == 0 ) {
                cout << " GLS level Scheme not correctly loaded " << endl; return;
        }
        
        // start an iterator on the list of links to find out a XGammaLink
        TIter iter(&random.GetLS().GetLinks()); 
        while ( obj = iter() ) { 
                if ( obj->InheritsFrom("XGammaLink") ) {
                        xlink = (XGammaLink *)obj; xlink->SetX(nbx,xe,xi);
                }
        }

        // some spectra to be filled
        TH1F *h = new TH1F("ProjTot","ProjTot",8192,0,4096);
        
        // start the Monte-Carlo
        Cascade cas; GammaLink *gam;
        for (Int_t i = 0; i < nbcas; i++ ) {

                random.DoCascade(&cas); // to get the next cascade

                for (Int_t j = 0; j < cas.GetSize(); j++ ) { // to set the total projection
                        gam = (GammaLink *)cas.At(j); h->Fill(gam->GetEnergy().Get());
                }
        }
        h->Draw();
}