GeantLMOF.cpp

Go to the documentation of this file.

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

#ifndef ROOT_TObjArray
#include "TObjArray.h"
#endif

#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,0)
#include "TList.h"
#include "TMath.h" 
#endif

#include "GeantLMOF.h"

#ifndef Gw_AgataGeometryTransformer
#include "AgataGeometryTransformer.h"
#endif

#define MAXEMITTED  200
#define MAXIMPACTS 1000

using namespace std;
using namespace Gw;


GeantLMOF::GeantLMOF()
{
        SetPath(); fCurrent = -1; fMaxFiles = 0; fStatus = kEmpty; fDistFactor = fEnerFactor = 1.0;

        fEnergy.reserve(MAXEMITTED); fDx.reserve(MAXEMITTED); fDy.reserve(MAXEMITTED); fDz.reserve(MAXEMITTED);
        fType.reserve(MAXEMITTED);

        fNdet.reserve(MAXIMPACTS); 
        fEvNum.reserve(MAXIMPACTS); 
        fNInt.reserve(MAXIMPACTS); 
        fPrimary.reserve(MAXIMPACTS);

        fEdep.reserve(MAXIMPACTS); 
        fPx.reserve(MAXIMPACTS); fPy.reserve(MAXIMPACTS); fPz.reserve(MAXIMPACTS);
        fNseg.reserve(MAXIMPACTS);

        fGamma = NULL; fFastOut = false; SetRecoil();
}

GeantLMOF::~GeantLMOF()
{
        if ( fGamma ) delete fGamma; fGamma = NULL;
}
        
TGeoVolume *GeantLMOF::ImportAgataG4(const char *asolid, const char *aclust, const char *aeuler, const char */*top*/)
{
  return AgataGeometryTransformer::ImportAgata(asolid, aclust, aeuler);
}

void GeantLMOF::ShowCurrentConditions()
{
        cout << "--> Energy Factor " << fEnerFactor << " and Distance Factor " << fDistFactor << endl; 
        cout << "--> Path to data:\n\t" << fPath << endl;
        cout << "--> Chain of files: " << endl;
        for (int i = 0; i < fMaxFiles; i++ ) {
                if ( i == fCurrent ) cout << " [c] ";
                else cout << " [ ] ";
                cout << fEventsRead[i] << " events read in " << fFileName[i] << endl;
        }
        cout << "--> Recoil: " << fRecoilBeta << " " << fRecoilSigma << " " << fRecoilDx << " " << fRecoilDy << " " << fRecoilDz << " " << fRecoilAngle << endl;
}

void GeantLMOF::ShowCurrentEvent()
{
        UInt_t i;

        cout << "list of emmited particles " << endl;
        for ( i = 0; i < fNbEmitted; i++ ) {
                cout << " " << fType[i] 
                << " " << fEnergy[i] 
                << " " << fDx[i]
                << " " << fDy[i]
                << " " << fDz[i]
                << " " << fEvNum[i] << endl;
        }
        cout << "list of impacts " << endl;
        for ( i = 0; i < fNbImpacts; i++ ) {
                cout << " " << fNdet[i] 
                << " " << fEdep[i] 
                << " " << fPx[i]
                << " " << fPy[i]
                << " " << fPz[i]
                << " " << fNseg[i] 
                << " " << fPrimary[i] << endl;
        }
}


void GeantLMOF::ReadHeader()
{
        fStatus = kGood;
        while( fInfile.good() ) { 

                // to get a new line
                for (int i = 0; i < MAXCHARACTER; i++) { buf[i] = 0; }
                fInfile.getline(buf,MAXCHARACTER);      

                TString id(buf,MAXCHARACTER);
                if ( id.Contains("RECOIL") ) {
                        istringstream theline(buf+6); 
                        theline.clear();
                        theline >> fRecoilBeta >> fRecoilSigma >> fRecoilDx >> fRecoilDy >> fRecoilDz >> fRecoilAngle;
                }
                if ( id.Contains("DISTFACTOR") ) {
                        istringstream theline(buf+10); 
                        theline.clear();
                        theline >> fDistFactor;
                }
                if ( id.Contains("ENERFACTOR") ) {
                        istringstream theline(buf+10); 
                        theline.clear();
                        theline >> fEnerFactor;
                }               
                // look for beginning of events 
                if ( buf[0] == '$' ) { return; }
        }
        cout << " Cannot find beginning of events $ !! " << endl;
        fInfile.close(); fInfile.clear(); fStatus = kFail; // nothing has been found, so status is kFail
}

void GeantLMOF::NextFile()
{
        if ( fMaxFiles == 0 ) { fStatus = kEmpty; return; }

        if ( fInfile.is_open() ) { fInfile.close(); fInfile.clear(); } // in case the user wants to switch in a middle of a file
        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.open(fullname.data());
                        
                        if ( fInfile.is_open() ) { ReadHeader(); if (fStatus==kGood ) { tmp = i; break; } }
                } // i
                if ( tmp == fCurrent ) { fStatus = kEoE; } // fCurrent is the last file
                else { fCurrent = tmp; }
        }
        else { fStatus = kEoE; return; }
}

UInt_t GeantLMOF::OnlyType(const GeantLMOF &input, Int_t typemin, Int_t typemax)
{       
        fNbImpacts = fNbEmitted = 0; 
        for (UInt_t i = 0; i < input.fNbImpacts; i++) { // to remove particles that are not of selected Types
                if ( input.fType[input.fPrimary[i]] > typemin && input.fType[input.fPrimary[i]] < typemax + 1 ) {
                        fNdet[fNbImpacts] = input.fNdet[i];
                        fEdep[fNbImpacts] = input.fEdep[i];
                        fPx[fNbImpacts]   = input.fPx[i];
                        fPy[fNbImpacts]   = input.fPy[i];
                        fPz[fNbImpacts]   = input.fPz[i];
                        fNseg[fNbImpacts]  = input.fNseg[i];
                        fNInt[fNbImpacts] = input.fNInt[i]; 
                        
                        bool anew = true;
                        for(UInt_t j = 0; j < fNbEmitted; j++ ) { // keep a copy of the emitted if not already done
                                if ( input.fEvNum[i] == fEvNum[j] ) { anew = false; fPrimary[fNbImpacts] = j; break; } 
                        }
                        if ( anew ) {
                                fType[fNbEmitted]       = input.fType[input.fPrimary[i]]; 
                                fEnergy[fNbEmitted]     = input.fEnergy[input.fPrimary[i]];
                                fDx[fNbEmitted]         = input.fDx[input.fPrimary[i]];
                                fDy[fNbEmitted]         = input.fDy[input.fPrimary[i]];
                                fDz[fNbEmitted]         = input.fDz[input.fPrimary[i]];
                                fEvNum[fNbEmitted]      = input.fEvNum[input.fPrimary[i]];      
                                
                                fPrimary[fNbImpacts]    = fNbEmitted++;
                        }
                        fNbImpacts++;
                }
        } // i          

        return fNbImpacts;
}

void GeantLMOF::operator>> (GammaTracker *tracker) 
{
        if ( tracker == NULL ) { return; }
                
        // to convert GeantLMOF units in tracker units
        Double_t efactor = GetEnerFactor() / tracker->GetEnerFactor();
        Double_t pfactor = GetDistFactor() / tracker->GetDistFactor();
        
        if ( fFastOut == true ) { // supposed to be only gammas
                tracker->SetEvent(fNbImpacts,&fEdep[0],&fPx[0],&fPy[0],&fPz[0],efactor,pfactor);
        }
        else  { // to be sure only impacts coming from gammas are there
                if ( fGamma == NULL ) {
                        fGamma = new GeantLMOF(); 
                }
                fGamma->OnlyType((*this)); 
                tracker->SetEvent(fGamma->NbImpacts(),&fGamma->GetEdep()[0],
                                  &fGamma->GetPx()[0],&fGamma->GetPy()[0],&fGamma->GetPz()[0],efactor,pfactor); 
        }
}


Bool_t GeantLMOF::NextEvent(UInt_t asked_mult)
{
        Int_t id; 

        // 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.is_open() || !fInfile.good() || fStatus != kGood ) {  NextFile(); if ( fStatus != kGood ) return kFALSE; }

        // Reset the current event
        UInt_t mult = asked_mult;
        fNbEmitted = 0;
        fNbImpacts = 0; 
        while( fInfile.good() ) { 

                Long_t fp = fInfile.tellg();

                // to get the next line 
                for (int i = 0; i < MAXCHARACTER; i++) { buf[i] = 0; }
                fInfile.getline(buf,MAXCHARACTER); istringstream theline(buf); theline.clear(); 

                theline >> id; 
                if ( id == -100 && fNbEmitted > 0 ) { // next event .. so stop reading this one
                        fEventsRead[fCurrent]++; 
                        break;
                }
                // just to be sure that one are not trying to build events on an already formated file 
                if ( id == -100  ) { mult = kMaxUInt; continue; } 

                if ( id == -101 ) { // source velocity
                        if ( gDebug == 1 ) cout << "New source velocity (-101) ... not yet treated " << endl; continue;
                }
                if ( id == -102 ) { // source position
                        if ( gDebug == 1 ) cout << "New source position (-102) ... not yet treated " << endl; continue;
                }
                if ( id == -103 ) { // fTime emission
                        if ( gDebug ) cout << "New fTime emission (-103) ... not yet treated " << endl; continue;
                }
                if ( id >= 0 ) { // id is the detector number
                        if ( gDebug == 1 ) cout << " Start new impact " << fNbEmitted << endl;

                        Double_t e, x, y, z; 
                        Int_t n;
                        theline >> e; theline >> x >> y >> z; theline >> n;

                        fNdet[fNbImpacts]       = id;
                        fEdep[fNbImpacts]       = e;
                        fPx[fNbImpacts]         = x;    
                        fPy[fNbImpacts]         = y;
                        fPz[fNbImpacts]         = z;
                        fNseg[fNbImpacts]       = n;
//                      fPrimary[fNbImpacts]    = fEvNum[fNbEmitted-1];
                        fPrimary[fNbImpacts]    = fNbEmitted-1;

                        fNbImpacts++;

                        if ( gDebug == 1 ) cout << "End new impact " <<  endl;
                }
                if ( id < 0 ) { // a new particle has been emitted
                        if ( gDebug == 1 ) cout << "Start new particle " << endl; 

                        Double_t e, x, y, z; 
                        Int_t n;
                        theline >> e; theline >> x >> y >> z; theline >> n;
                        
                        fType[fNbEmitted]       = -id;
                        fEnergy[fNbEmitted]     = e;
                        fDx[fNbEmitted]         = x;    
                        fDy[fNbEmitted]         = y;
                        fDz[fNbEmitted]         = z;
                        fEvNum[fNbEmitted]      = n;

                        fNbEmitted++;

                        if ( gDebug == 1 ) cout << "End new particle " << id << " " << fNbEmitted << endl; 
                }
                // to read files without beginning of events flags ... rewind the last line
                if ( fNbEmitted > mult ) { fInfile.seekg(fp,ios::beg); fNbEmitted--; fEventsRead[fCurrent]+=1.0; break; }

                // to avoid out_of_range in vectors
                if ( fNbEmitted > MAXEMITTED - 1 || fNbImpacts > MAXIMPACTS - 1 ) { 
                        cout << " WARNING - Too much emitted particles " << fNbEmitted << " or impacts " << fNbImpacts << endl; 
                        fStatus = kFail; 
                        return kFALSE; 
                }
        } // while in_good
        if ( fNbEmitted > 0 ) return kTRUE;
        return kFALSE;
}