FirstTracking.C

Go to the documentation of this file.

#ifndef __CINT__

#ifndef ROOT_TH1
#include "TH1.h"
#endif
#ifndef ROOT_TH2
#include "TH2.h"
#endif
#ifndef ROOT_TGeoManager
#include "TGeoManager.h"
#endif
#ifndef ROOT_TGeoTrack
#include "TGeoTrack.h"
#endif

#endif

#include "GeantLMOF.h"
#include "GammaTracker.h"
#include "XtermProgressBar.h"

using namespace Gw;

const char TRACKERNAME[] = "OFT"; // Tracker you would like to use for in these functions

void EPTracking(const char *input, UInt_t mult = 1, UInt_t MaxEvents = 1000000)
{
        const Int_t window = 2; // window to integrate peaks in the output spectrum
        
// load the tacker you would like to use        
        GammaTracker *tracker = GammaTracker::GetTracker(TRACKERNAME);
        if ( tracker == NULL ) {
                cout << " Cannot run, unknown GammaTracker " << TRACKERNAME << endl; return;
        }
        
// input file(s) decoded by GeantLMOF reader
        GeantLMOF reader; reader << fastout; reader.AddFile(input); 

// spectra definitions
        const int BINRANGE = 2000; const Stat_t MINX = 0.; const Stat_t MAXX = 2000.;
        
        TH1F *hinput  = new TH1F("Input","From Geant",BINRANGE,MINX,MAXX);
        TH1F *houtput = new TH1F("Output","After Tracking",BINRANGE,MINX,MAXX);
        
        TH1F *ep = new TH1F("PhotopeakEfficiency","PhotopeakEfficiency",BINRANGE,MINX,MAXX);

// start reading
        XtermProgressBar progress(0,MaxEvents);
        
        while( reader.NextEvent(mult) ) { // loop on events
                
                progress++; 

                // fill the tracker with the current event
                reader >> tracker; 
                for (UInt_t i = 0; i < reader.NbEmitted(); i++ ) { // fill input
                        hinput->Fill(reader.GetEnerFactor()*reader.GetEnergy()[i]);
                }
                
                // do the tracking
                tracker->DoTracking();
                for (UInt_t i = 0; i < tracker->GetOutputN(); i++ ) { // fill output
                        houtput->Fill(tracker->GetEnerFactor()*tracker->GetOutputE()[i]);
                }               
                
                if ( progress.IsFinished() ) break;

        } // while reader
        cout << endl; reader.ShowCurrentConditions(); 
        
// now compute photopeak efficiency
        for(Int_t i = 0 ; i < BINRANGE ; i++ ){
                if ( hinput->GetBinContent(i+1) > 1 ) {
                        Stat_t peak = houtput->Integral(i+1-window,i+1+window) - 
                                        0.5*(houtput->Integral(i-2*window,i-window)+houtput->Integral(i+window,i+2*window));
                        
                        ep->SetBinContent(i+1,peak/hinput->GetBinContent(i+1));
                }
        }
        
        delete tracker;
}

void GGTracking(const char *input, UInt_t mult = 1, UInt_t MaxEvents = 1000000)
{       
// load the tacker you would like to use        
        GammaTracker *tracker = GammaTracker::GetTracker(TRACKERNAME);
        if ( tracker == NULL ) {
                cout << " Cannot run, unknown GammaTracker " <<  TRACKERNAME << endl; return;
        }
        
// input file(s) decoded by GeantLMOF reader
// be carefull, fastout is called, remove this call is your GEANT file contains other particles that gammas
        GeantLMOF reader; reader << fastout; reader.AddFile(input); 

// matrix definition
        const int BINRANGE = 2000; const Stat_t MINX = 0.; const Stat_t MAXX = 2000.;
        TH2F *matrix = new TH2F("GxG","Gamma-Gamma Matrix",BINRANGE, MINX, MAXX, BINRANGE, MINX, MAXX);

// start reading
        XtermProgressBar progress(0,MaxEvents);
        
        while( reader.NextEvent(mult) ) { // loop on events
                
                progress++; 
                
                // fill the tracker with the current event
                reader >> tracker; 
                
                // do the tracking
                tracker->DoTracking(); if ( tracker->GetOutputN() < 2 ) continue;

                for (UInt_t i = 0; i < tracker->GetOutputN(); i++ ) {
                        Float_t ecos1 = TMath::Cos(tracker->GetOutputTheta1()[i]);
                        Float_t edop1 = (tracker->GetEnerFactor()*tracker->GetOutputE()[i]) 
                                        / ( 1 + reader.GetRecoilBeta() * ecos1 ) ;
                                
                        for (UInt_t j = i+1; j < tracker->GetOutputN(); j++ ) {
                                
                                Float_t ecos2 = TMath::Cos(tracker->GetOutputTheta1()[j]);
                                Float_t edop2 = (tracker->GetEnerFactor()*tracker->GetOutputE()[j]) 
                                        / ( 1 + reader.GetRecoilBeta() * ecos2 ) ;
                                
                                matrix->Fill(edop1,edop2); matrix->Fill(edop2,edop1); 
                        } // j 
                } // i
                
                if ( progress.IsFinished() ) break;
        } // while reader
        reader.ShowCurrentConditions(); 
        
        delete tracker;
}

void TestImportAgata()
{
        // the world is created first
        TGeoManager  *geom = new TGeoManager("Agata","AGATA");
        TGeoMaterial *mat  = new TGeoMaterial("Vacuum",0,0,0); TGeoMedium *med = new TGeoMedium("Vacuum",1,mat);
        
        TGeoVolume *top = gGeoManager->MakeBox("World",med,100,100,100); gGeoManager->SetTopVolume(top);
                
        // agata is added to the world
        Gw::GeantLMOF::ImportAgataG4();
        
        // the geometry MUST be closed once defined !
        geom->CloseGeometry();
}

void VisuTrack(const char *input, UInt_t mult = 1, UInt_t MaxEvents = 10)
{
// the world is created first
        TGeoManager  *geom = new TGeoManager("Agata","AGATA");
        TGeoMaterial *mat  = new TGeoMaterial("Vacuum",0,0,0); TGeoMedium *med = new TGeoMedium("Vacuum",1,mat);
        
        TGeoVolume *top = gGeoManager->MakeBox("World",med,100,100,100); gGeoManager->SetTopVolume(top);
                
        // agata is added to the World
        Gw::GeantLMOF::ImportAgataG4();
        
        // the geometry MUST be closed once defined !
        geom->CloseGeometry();
        
// input file(s) decoded by GeantLMOF reader
        GeantLMOF reader; reader << fastout; reader.AddFile(input); 
        
// start reading
        XtermProgressBar progress(0,MaxEvents); 
        
        while( reader.NextEvent(mult) ) { // loop on events
        
                progress++; if ( reader.NbImpacts() < 1 ) continue;

                Int_t trackid = geom->AddTrack(Int_t(progress.GetCurrentPosition()),22); TVirtualGeoTrack *primary = geom->GetTrack(trackid);
                
                primary->AddPoint(0,0,0,0);
                for( UInt_t i = 0; i < reader.NbImpacts() ; i++ ){
                        primary->AddPoint(reader.GetPx()[i]/10,reader.GetPy()[i]/10,reader.GetPz()[i]/10,0);
                        
                }
                if ( progress.IsFinished() ) break;     
        }       
        geom->GetTopVolume()->Draw("ogl");
}