OneLStoSToGS.C

Go to the documentation of this file.

// Macro to produce from a level scheme, random gamma-ray cascades.
// The Monte- Carlo could be customize through parameters.
// In particular one can ask the random generator (ADTy[e == 1) to produce randomly directions of emissions
// based on the initial, final level and the gamma-ray linking them

//  In this case the root mathmore library is mandatory to have access to wigner_3j and 6j
//   thus root should be build with this option and then gammaware built on top
//   in case mathmore does not exists, all gamma-rays are emitted isotropically
//
// The output is a root tree containing the Pr [Primary] branch
// Such a file could be used as input in SToGS
//
// usage:
// root -l Load.C
// .L OneLStoSToGS.C++
// OneLSWithAngularDistributions("89Y.root","89Y")

// O. Stézowski

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

#include "TStopwatch.h"
#include "TFile.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include "TDatabasePDG.h"
#include "TH1F.h"

#include <SToGS_BaseROOTEvents.h>

void OneLSWithAngularDistributions(const char *lsname, const char *lstitle, Int_t nbcas = 100000, const char *treefilename = "PrimariesFromGW.root")
{
    // Random engine is TRandom3 based on clock ... just remove clock to have identical random sequences
    Gw::Random::Instance()->SetCurrent("large_clock");
    
        TStopwatch watch;
        watch.Start();
    
    // build gem, customize it through paramters and then init with a level scheme
    Gw::BaseGEM gem;
    gem.SetParameter("ADType",1);           // it means do angular distribution for all g-rays. 0 means isotropic for all
    //gem.SetParameter("CutLifeTime",100.0);  // below cut, do angular. Up isotropic
    gem.Import(lsname,lstitle);
    
    // ROOT file for output tree
    TFile rf(treefilename,"RECREATE"); TTree *tree_primary = new TTree("GWLS_to_SToGS","Primaries from GW");
    SBRPEvent pr_event;
    //
    tree_primary->Branch("Pr.",&pr_event);
    
    Gw::Cascade cas; TObjArray directions; Gw::GammaLink *gam; Int_t gPDG = TDatabasePDG::Instance()->GetParticle("gamma")->PdgCode();
    
        for (Int_t i = 0; i < nbcas; i++ ) {
        
        // GAMMA-rays from LS with angular distributions
                gem.DoCascade(cas,directions);
        
        Double_t ESUM = 0.0, MULTTOT = 0; pr_event.Clear("");
        
        MULTTOT = cas.GetSize();
        for (Int_t j = 0; j < MULTTOT; j++ ) {
            gam = (Gw::GammaLink *)cas.At(j);
            
            SBRPHit *stogshit = pr_event.AddHit();
            
            // now copy gam content into stogshit
            stogshit->fE = gam->GetEnergy().Get();
            stogshit->fX = 0.;  // X position
            stogshit->fY = 0.;  // Y position
            stogshit->fZ = 0.;  // Z position
            
            TLorentzVector *lv = (TLorentzVector *)directions.At(j);
            stogshit->fPX = lv->X();
            stogshit->fPY = lv->Y();
            stogshit->fPZ = lv->Z();
            stogshit->fT = lv->T();     //
            
            stogshit->fPDG  = gPDG;
            stogshit->fFlag = j;
            stogshit->fUID  = 0;
            
            ESUM += stogshit->fE;
        
                } // j
        
        pr_event.SetEMult(ESUM,MULTTOT);
        tree_primary->Fill();
    }
    
    tree_primary->Print();
    tree_primary->Write();
    
    watch.Stop();
    cout << "Monte-Carlo time " << endl; watch.Print();
}

/* Initial function used to test BaseGEM with angular distributions
void TestADBaseGEM(Int_t nbcas = 100, const char *rfilename = "PrimariesFromGW.root")
{
    TH1F *testE2 = new TH1F("TestE2","Test",1000,0,180);
    TH1F *testM1 = new TH1F("TestM1","Test",1000,0,180);
    TH1F *testISO = new TH1F("TestISO","Test",1000,-1,1);
    //TH1F *testISO = new TH1F("TestISO","Test",1000,-180,180);
    
    // Random engine is TRandom3 based on clock ... just remove clock to have identical random sequences
    Gw::Random::Instance()->SetCurrent("large_clock");
    
        TStopwatch watch;
        watch.Start();
    
    // build gem, set level scheme
    Gw::BaseGEM gem;
    gem.SetParameter("ADType",1);           // it means do angular distribution for all g-rays. 0 means isotropic for all
    //gem.SetParameter("CutLifeTime",100.0);  // below cut, do angular. Up isotropic
    gem.Import("89Y.root","89Y");
    
    // ROOT file for output tree
    TFile rf(rfilename,"RECREATE"); TTree *tree_primary = new TTree("GWLS_to_SToGS","Primaries from GW");
    SBRPEvent pr_event;
    //
    tree_primary->Branch("Pr.",&pr_event);

    Gw::Cascade cas; TObjArray directions; Gw::GammaLink *gam; Int_t gPDG = TDatabasePDG::Instance()->GetParticle("gamma")->PdgCode();

        for (Int_t i = 0; i < nbcas; i++ ) {
        
        // GAMMA-rays from LS with angular distributions
                gem.DoCascade(cas,directions);
        
        Double_t ESUM = 0.0, MULTTOT = 0; pr_event.Clear("");
        
        MULTTOT = cas.GetSize();
        for (Int_t j = 0; j < MULTTOT; j++ ) {
            gam = (Gw::GammaLink *)cas.At(j);
            
            SBRPHit *stogshit = pr_event.AddHit();
            
            // now copy gam content into stogshit
            stogshit->fE = gam->GetEnergy().Get();
            stogshit->fX = 0.;  // X position
            stogshit->fY = 0.;  // Y position
            stogshit->fZ = 0.;  // Z position
            
            TLorentzVector *lv = (TLorentzVector *)directions.At(j);
            stogshit->fPX = lv->X();
            stogshit->fPY = lv->Y();
            stogshit->fPZ = lv->Z();
            stogshit->fT = lv->T();     //
            
            stogshit->fPDG  = gPDG;
            stogshit->fFlag = j;
            stogshit->fUID  = 0;
            
            // hegamma->Fill(stogshit->fE);
            ESUM += stogshit->fE;
            
            if (stogshit->fE >1743 && stogshit->fE < 1745 ) {
                testE2->Fill(lv->Theta()*180./TMath::Pi());
            }
            if (stogshit->fE >1506 && stogshit->fE < 1508 ) {
                testM1->Fill(lv->Theta()*180./TMath::Pi());
            }
            //if (stogshit->fE >4014 && stogshit->fE < 4016 ) {
            if (stogshit->fE >908 && stogshit->fE < 910 ) {
                //testISO->Fill(lv->Theta()*180./TMath::Pi());
                testISO->Fill(TMath::Cos(lv->Theta()) );
                //testISO->Fill(TMath::Cos(lv->Theta()) );

                //testISO->Fill(lv->Phi()*180./TMath::Pi());
            }
            
                } // j
        
        pr_event.SetEMult(ESUM,MULTTOT);
        tree_primary->Fill();
    }
    
    tree_primary->Print();
    tree_primary->Write();
    
    watch.Stop();
    cout << "Monte-Carlo time " << endl; watch.Print();
    
}
*/