BmnEcalDigitizer.cxx

Go to the documentation of this file.

/* 
 * File:   BmnEcalDigitizer.cxx
 * Author: Petr Alekseev
 * 
 * Created on 05.08.2020, 19:31
 */
 
#include "BmnEcalDigitizer.h"
#include "BmnEcalPoint.h"
 
#include "FairLogger.h"
#include "FairRootManager.h"
 
#include "TGeoManager.h"
 
#include <iostream>
 
BmnEcalDigitizer::BmnEcalDigitizer() {
}
 
BmnEcalDigitizer::~BmnEcalDigitizer() {
}
 
InitStatus BmnEcalDigitizer::Init() {
 
    FairRootManager* ioman = FairRootManager::Instance();
    fArrayOfEcalPoints = (TClonesArray*) ioman->GetObject("EcalPoint");
    if (fArrayOfEcalPoints == nullptr)
    {
        LOG(error)<<"BmnEcalDigitizer::Init() branch 'EcalPoint' not found! Task will be deactivated";
        SetActive(kFALSE);
        return kERROR;
    }
 
    fArrayOfEcalDigits = new TClonesArray("BmnECALDigit");
    ioman->Register("EcalDigit", "Ecal", fArrayOfEcalDigits, kTRUE);
    
    if (LoadGeometry() != 0)
    {
        SetActive(kFALSE);
        return kERROR;
    }
    
    Info(__func__,"ECAL digitizer ready");
    return kSUCCESS;
    
}
 
void BmnEcalDigitizer::Exec(Option_t* opt)
{
    if (!IsActive()) return;
    
    // Initialize
    fArrayOfEcalDigits->Delete();
 
    for (Int_t i = 0; i < fCellsSize; i++) {
        fCells[i].SetAmp(0.);
        fCells[i].SetStartTime(0.);
    }
    
    // Collect points
    Int_t N = fArrayOfEcalPoints->GetEntries();
    for (Int_t i = 0; i < N; i++) {
        BmnEcalPoint * p = (BmnEcalPoint *)fArrayOfEcalPoints->At(i);
        
        Int_t ch = p->GetCopyMother();
        if (ch < fCellsSize) {
            if (fCells[ch].GetChannel() == ch) {
                Float_t energyLoss = p->GetEnergyLoss();
                if (energyLoss > 0) {
                    Float_t pointTime = p->GetTime();
                    if (fMaxPointTime > 0 && pointTime > fMaxPointTime) continue;
                    if (fFiberSOL > 0.) {
                        pointTime += (fFiberLength - fLayerThickness * (p->GetCopy()+1))/fFiberSOL; 
                    }
                    fCells[ch].SetAmp(fCells[ch].GetAmp() + energyLoss);
                    fCells[ch].SetStartTime(fCells[ch].GetStartTime() + pointTime * energyLoss);
                }
            } else {
                Error(__func__, "ECAL ch %d was not initialized",ch);
            }   
        } else {
            Error(__func__,"ECAL ch %d ignored",ch);
        }
    }
    
    // Store digits
    for (Int_t i = 0; i < fCellsSize; i++) {
        if (fCells[i].GetChannel() != i) continue;
        
        Float_t amp = fCells[i].GetAmp() * 1000.;
        
        if (amp == 0.) continue;
        if (amp < fThreshold) continue;
        
        Float_t time = fCells[i].GetStartTime()/fCells[i].GetAmp();
        
        BmnECALDigit * p = new((*fArrayOfEcalDigits)[fArrayOfEcalDigits->GetEntriesFast()]) BmnECALDigit();
        *p = fCells[i];
        p->SetAmp(amp);
        p->SetPeakAmp(amp);
        p->SetStartTime(time);
        p->SetPeakTime(time + fPeakTimeDelay);
    }
    
}
 
void BmnEcalDigitizer::Print(Option_t *option) const {
 
    printf("\n" "\e[1;92m");
    printf("BMN ECAL Digitizer\n");
    printf("ECAL geometry fileN: %s\n", fEcalGeometryFileName);
    printf("ECAL interaction depth shift: %f\n", fDepthShift);    
    printf("ECAL cells coords (((\n");
    Float_t x, y, z;
    for (Int_t i = 1; i < fCellsSize; i++) {
        if (fCells[i].GetChannel() == i) {
            fCells[i].GetLabCoords(x,y,z);
            printf("    % 4d [%3.0f,%3.0f] -> [%7.2f,%7.2f,%7.2f]\n", i, fCells[i].GetX(), fCells[i].GetY(), x, y, z);
        } else {
            //printf("    % 04d - Not initialized\n", i);
        }
    }
    printf("))) ECAL cells coords\n"); 
    printf("\e[0m" "\n");
    
}
 
int BmnEcalDigitizer::LoadGeometry()
{
    Bool_t loadFromFile = fEcalGeometryFileName != nullptr;
    TGeoNode *ecal1 = 0, *ecal2 = 0;
    if (gGeoManager) {
        TGeoVolume * ecal = gGeoManager->FindVolumeFast("ecal");
        if (ecal) {
            if (ecal->GetNdaughters() > 0) ecal1 = ecal->GetNode(0);
            if (ecal->GetNdaughters() > 1) ecal2 = ecal->GetNode(1);
            loadFromFile = kFALSE;
        } else {
            Info(__func__, "Ecal geometry not found by TGeoManager\n");
        }
    } 
    
    if (loadFromFile) {
        Info(__func__, "Loading coordinates of ECAL cells from %s\n", fEcalGeometryFileName);
        TGeoVolume * top = TGeoVolume::Import(fEcalGeometryFileName,"TOP");
 
        if (!top) {
            LOG(error)<<"<BmnEcalDigitizer::LoadGeometry>: Volume TOP not found in "<<(fEcalGeometryFileName == nullptr? "(null)" : fEcalGeometryFileName);
            //Fatal(__func__, "Volume TOP not found in %s\n", fEcalGeometryFileName);
            return -1;
        }
        
        TGeoNode * ecal = top->GetNode(0);
    
        if (!ecal /*|| ecal->GetNdaughters() < 2*/) {
            LOG(error)<<"<BmnEcalDigitizer::LoadGeometry>: Unexpected geometry structure "<<(fEcalGeometryFileName == nullptr? "(null)" : fEcalGeometryFileName);
            //Fatal(__func__, "Unexpected geometry structure %s\n",fEcalGeometryFileName);
            return -1;
        }        
        ecal1 = ecal->GetDaughter(0);
        ecal2 = ecal->GetDaughter(1);
    }
    
    if (ecal1 == 0 && ecal2 == 0) {
        Info(__func__, "BmnEcalDigitizer was set in inactive state, because ECAL geometry was not found");
        return -1;
    }
    
    for (Int_t i = 0; i < fCellsSize; i++) {
        fCells[i].SetChannel(0);
    }
 
    Double_t coords[] = {0.,0.,0.};
    Double_t ecalCoords[3];
    Double_t labCoords[3];
    
    coords[2] = fDepthShift;
  
    if (ecal1) {
        Int_t n = ecal1->GetNdaughters();
        if (n > 504) {
            LOG(error)<<"<BmnEcalDigitizer::LoadGeometry>: Expected ecal node 1 with 504 daughters or less in "<<(fEcalGeometryFileName == nullptr? "(null)" : fEcalGeometryFileName);
            return -1;
        }
        for (Int_t i = 0; i < n; i++) {
            Int_t ch = ecal1->GetDaughter(i)->GetNumber();
            if (ch < 1 || ch > 504) {
                LOG(error)<<"<BmnEcalDigitizer::LoadGeometry>: Unexpected chan "<<ch<<" at ecal node 1 in "<<(fEcalGeometryFileName == nullptr? "(null)" : fEcalGeometryFileName);
                return -1;
            }
            ecal1->GetDaughter(i)->LocalToMaster(coords, ecalCoords);
            ecal1->LocalToMaster(ecalCoords,labCoords);
            fCells[ch].SetChannel(ch);
            fCells[ch].SetX(ecalCoords[0]);
            fCells[ch].SetY(ecalCoords[1]);
            fCells[ch].SetLabCoords(labCoords[0],labCoords[1],labCoords[2]);
        }
    }
 
    if (ecal2) {
        Int_t n = ecal2->GetNdaughters();
        if (n > 504) {
            LOG(error)<<"<BmnEcalDigitizer::LoadGeometry>: Expected ecal node 2 with 504 daughters or less in "<<(fEcalGeometryFileName == nullptr? "(null)" : fEcalGeometryFileName);
            return -1;
        }
        for (Int_t i = 0; i < n; i++) {
            Int_t ch = ecal2->GetDaughter(i)->GetNumber();
            if (ch < 505 || ch > 1008) {
                LOG(error)<<"<BmnEcalDigitizer::LoadGeometry>: Unexpected chan="<<ch<<" at ecal node 2 in "<<(fEcalGeometryFileName == nullptr? "(null)" : fEcalGeometryFileName);
                return -1;
            }
            ecal2->GetDaughter(i)->LocalToMaster(coords, ecalCoords);
            ecal2->LocalToMaster(ecalCoords,labCoords);
            fCells[ch].SetChannel(ch);
            fCells[ch].SetX(ecalCoords[0]);
            fCells[ch].SetY(ecalCoords[1]);
            fCells[ch].SetLabCoords(labCoords[0],labCoords[1],labCoords[2]);
        }
    }
 
    return 0;
}