BmnHodo.cxx

Go to the documentation of this file.

#include "BmnHodo.h"
 
#include "BmnHodoGeo.h"
#include "BmnHodoPoint.h"
#include "CbmStack.h"
#include "FairGeoInterface.h"
#include "FairGeoLoader.h"
#include "FairGeoNode.h"
#include "FairGeoRootBuilder.h"
#include "FairMCPoint.h"
#include "FairRootManager.h"
#include "FairVolume.h"
#include "TClonesArray.h"
#include "TGeoArb8.h"
#include "TGeoMCGeometry.h"
#include "TGeoManager.h"
#include "TLorentzVector.h"
#include "TParticle.h"
#include "TVector3.h"
#include "TVirtualMC.h"
 
#include <iostream>
// add on for debug
// #include "FairGeoG3Builder.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "TObjArray.h"
#include "TParticlePDG.h"
 
// -----   Default constructor   -------------------------------------------
BmnHodo::BmnHodo()
{
    fHodoCollection = new TClonesArray("BmnHodoPoint");
    volDetector = 0;
    fEventID = -1;
    fHodoStickSensVolId = 0;
}
 
// -----   Standard constructor   ------------------------------------------
BmnHodo::BmnHodo(const char* name, Bool_t active)
    : FairDetector(name, active)
{
    fHodoCollection = new TClonesArray("BmnHodoPoint");
    volDetector = 0;
    fEventID = -1;
    fHodoStickSensVolId = 0;
    fVerboseLevel = 1;
}
 
// -----   Destructor   ----------------------------------------------------
BmnHodo::~BmnHodo()
{
    if (fHodoCollection) {
        fHodoCollection->Delete();
        delete fHodoCollection;
    }
}
 
// -----   Public method Intialize   ---------------------------------------
void BmnHodo::Initialize()
{
    // Init function
 
    FairDetector::Initialize();
    // FairRun* sim = FairRun::Instance();
    // FairRuntimeDb* rtdb = sim->GetRuntimeDb();
 
    fHodoStickSensVolId = gMC->VolId("hodo_stick_sens");
}
 
// -------------------------------------------------------------------------
BmnHodoPoint* BmnHodo::GetHit(Int_t i) const
{
    // Returns the hit for the specified layer.
    // ---
 
    return (BmnHodoPoint*)fHodoCollection->At(i);
}
 
//_____________________________________________________________________________
BmnHodoPoint* BmnHodo::GetHit(Int_t sticksens, Int_t stick) const
{
    // Returns the hit for the specified vsc and module.
    // ---
 
    BmnHodoPoint* hit;
    Int_t nofHits = fHodoCollection->GetEntriesFast();
    for (Int_t i = 0; i < nofHits; i++) {
        hit = GetHit(i);
        if (hit->GetCopy() == sticksens && hit->GetCopyMother() == stick)
            return hit;
    }
 
    return 0;
}
 
// -----   Public method ProcessHits  --------------------------------------
Bool_t BmnHodo::ProcessHits(FairVolume* vol)
{
    /** This method is called from the MC stepping */
 
    Int_t copyNoStickSens;
    Int_t copyNoSTICKSENSCom, copyNoSTICKCom;
 
    Int_t ivol;
    TLorentzVector tPos1, tMom1;
    TLorentzVector tPos, tMom;
 
    Int_t sticksens, stick;
 
    Double_t time = 0;
    Double_t length = 0;
 
    // TParticle* part;
    // Double_t charge;
 
    // Double_t  QCF=1; //quenching for Birk
    // Double_t  BirkConst = 12.6; //0.126 mm/MeV for polystyrene
    // 0.126 *(0.1/0.001) = 12.6 cm/GeV
    //(0.126 mm/MeV - from Wikipedia, 0.07943mm/MeV є– in Geant4)
 
    if (gMC->CurrentVolID(copyNoStickSens) != fHodoStickSensVolId)
        return kFALSE;
 
    ivol = vol->getMCid();
 
    if (gMC->CurrentVolID(copyNoStickSens) == fHodoStickSensVolId) {
        gMC->CurrentVolOffID(0, sticksens);
        gMC->CurrentVolOffID(1, stick);
        copyNoSTICKSENSCom = sticksens;
        copyNoSTICKCom = stick;
    }
 
    // cout <<"2_1_0 " << ivol <<" " <<gMC->CurrentVolOffName(2) << " " <<gMC->CurrentVolOffName(1) << " " <<
    // gMC->CurrentVolOffName(0) <<std::endl;
 
    if (gMC->IsTrackEntering()) {
 
        ResetParameters();
 
        fELoss = 0.;
        time = 0.;
        length = 0.;
 
        TLorentzVector PosIn;
        gMC->TrackPosition(PosIn);
        fPos.SetXYZ(PosIn.X(), PosIn.Y(), PosIn.Z());
 
        TLorentzVector MomIn;
        gMC->TrackMomentum(MomIn);
        fMom.SetXYZ(MomIn.Px(), MomIn.Py(), MomIn.Pz());
 
        fTime = gMC->TrackTime() * 1.0e09;
        fLength = gMC->TrackLength();
    }   // if (gMC->IsTrackEntering())
 
    if (gMC->IsTrackInside()) {
 
        gMC->TrackPosition(tPos);
        gMC->TrackMomentum(tMom);
        length += gMC->TrackStep();
 
        fELoss += gMC->Edep();
        /*
//Birk corrections
    if(gMC->TrackStep()>0) QCF = 1.+(BirkConst/gMC->TrackStep())*gMC->Edep();
    else QCF = 1;
    fELoss +=(gMC->Edep())/QCF;
    */
        time += gMC->TrackTime() * 1.0e09;   // nsec
 
        if (gMC->IsTrackStop() || gMC->IsTrackDisappeared()) {
 
            // part    = gMC->GetStack()->GetCurrentTrack();
            // charge = part->GetPDG()->Charge() / 3. ;
 
            // Create BmnZdcPoint
            fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
            // time    = gMC->TrackTime() * 1.0e09;
            // length  = gMC->TrackLength();
            // gMC->TrackPosition(tPos);
            // gMC->TrackMomentum(tMom);
            // gMC->CurrentVolOffID(2, module);
            // gMC->CurrentVolOffID(1, slice);
 
            if (fELoss > 0) {
                if (copyNoSTICKSENSCom == sticksens && copyNoSTICKCom == stick) {
                    if (!GetHit(sticksens, stick)) {
                        AddHit(fTrackID, ivol, sticksens, stick, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
                               TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time, length, fELoss);
                    } else {
                        GetHit(sticksens, stick)
                            ->AddSTICK(fTrackID, ivol, sticksens, stick, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
                                       TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time, length, fELoss);
                    }
                }   // if(copyNoSLICECom==slice && copyNoCELLCom==cell)
 
            }   // if(fELoss>0)
        }   // if ( gMC->IsTrackStop() || gMC->IsTrackDisappeared() )
    }   // if ( gMC->IsTrackInside())
 
    if (gMC->IsTrackExiting()) {
        // part  = gMC->GetStack()->GetCurrentTrack();
 
        // Create BmnZdcPoint
        fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
        time += gMC->TrackTime() * 1.0e09;
        length += gMC->TrackLength();
 
        fELoss += gMC->Edep();
        /*
//Birk corrections
    if(gMC->TrackStep()>0) QCF = 1.+(BirkConst/gMC->TrackStep())*gMC->Edep();
    else QCF = 1;
    fELoss +=(gMC->Edep())/QCF;
    */
 
        gMC->TrackPosition(tPos);
        gMC->TrackMomentum(tMom);
 
        if (fELoss > 0) {
            if (copyNoSTICKSENSCom == sticksens && copyNoSTICKCom == stick) {
                if (!GetHit(sticksens, stick)) {
                    AddHit(fTrackID, ivol, sticksens, stick, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
                           TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time, length, fELoss);
                } else {
                    GetHit(sticksens, stick)
                        ->AddSTICK(fTrackID, ivol, sticksens, stick, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
                                   TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time, length, fELoss);
                }
            }   // if(copyNoSLICECom==slice && copyNoCELLCom==cell)
 
        }   // if(fELoss>0)
    }   // if ( gMC->IsTrackExiting()) {
 
    Int_t points = gMC->GetStack()->GetCurrentTrack()->GetMother(1);
    points = (points & (~(1 << 30))) | (1 << 30);
 
    gMC->GetStack()->GetCurrentTrack()->SetMother(1, points);
 
    ((CbmStack*)gMC->GetStack())->AddPoint(kHODO);
 
    /*//S.Merts
    Int_t ivol = vol->getMCid();
    if (gMC->IsTrackEntering()) {
        ResetParameters();
 
        fELoss = 0.;
 
        TLorentzVector PosIn;
        gMC->TrackPosition(PosIn);
        fPos.SetXYZ(PosIn.X(), PosIn.Y(), PosIn.Z());
 
        TLorentzVector MomIn;
        gMC->TrackMomentum(MomIn);
        fMom.SetXYZ(MomIn.Px(), MomIn.Py(), MomIn.Pz());
 
        fTime = gMC->TrackTime() * 1.0e09;
        fLength = gMC->TrackLength();
    }
 
    fELoss += gMC->Edep();
    if (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) {
        fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
        fVolumeID = vol->getMCid();
        if (fELoss == 0.) {
            return kFALSE;
        }
        AddHit(fTrackID, fVolumeID, fPos, fMom, fTime, fLength, fELoss);
 
        ((CbmStack*)gMC->GetStack())->AddPoint(kHODO);
    }
*/
    return kTRUE;
}
 
// -----   Public method EndOfEvent   -----------------------------------------
void BmnHodo::EndOfEvent()
{
    if (fVerboseLevel)
        Print("");
    Reset();
}
 
// -----   Public method Register   -------------------------------------------
void BmnHodo::Register()
{
    FairRootManager::Instance()->Register("HodoPoint", "Hodo", fHodoCollection, kTRUE);
}
 
// -----   Public method GetCollection   --------------------------------------
TClonesArray* BmnHodo::GetCollection(Int_t iColl) const
{
    if (iColl == 0)
        return fHodoCollection;
    return NULL;
}
 
// -----   Public method Print   ----------------------------------------------
void BmnHodo::Print(Option_t*) const
{
    Int_t nHits = fHodoCollection->GetEntriesFast();
    cout << "-I- BmnHodo: " << nHits << " points registered in this event." << endl;
 
    if (fVerboseLevel > 1)
        for (Int_t i = 0; i < nHits; i++)
            (*fHodoCollection)[i]->Print();
}
 
// -----   Public method Reset   ----------------------------------------------
void BmnHodo::Reset()
{
    fHodoCollection->Delete();
}
 
//-----------------------------------------------------------------------------
void BmnHodo::ConstructGeometry()
{
    TString fileName = GetGeometryFileName();
    if (fileName.EndsWith(".root")) {
        LOG(info) << "Constructing Hodo geometry from ROOT file " << fileName.Data();
        ConstructRootGeometry();
    } else
        LOG(fatal) << "Geometry format of Hodo file " << fileName.Data() << " not supported.";
}
 
// -----   Private method AddHit   --------------------------------------------
BmnHodoPoint* BmnHodo::AddHit(Int_t trackID,
                              Int_t detID,
                              Int_t copyNo,
                              Int_t copyNoMother,
                              TVector3 pos,
                              TVector3 mom,
                              Double_t time,
                              Double_t length,
                              Double_t eLoss)
{
    TClonesArray& clref = *fHodoCollection;
    Int_t size = clref.GetEntriesFast();
    return new (clref[size]) BmnHodoPoint(trackID, detID, copyNo, copyNoMother, pos, mom, time, length, eLoss);
}