BmnECALRaw2Digit.cxx

Go to the documentation of this file.

#include "BmnECALRaw2Digit.h"
 
#include "TBox.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TRandom.h"
#include "TSystem.h"
#include "TText.h"
 
#include <TError.h>
#include <TF1.h>
#include <TMinuit.h>
#include <fstream>
 
#define SHIFT 0
 
using namespace std;
 
BmnECALRaw2Digit::BmnECALRaw2Digit(ULong_t runId)
{
    const char* dir = getenv("VMCWORKDIR");
 
    if (runId < 1170) {
        Error(__func__, "No ECAL in run %lu\n", runId);
        return;
    }
 
    Int_t run = 0;
    Int_t pos = 0;
 
    if (runId < 2000) {
        run = 6;
        pos = 1;
    } else if (runId < 5186) {
        run = 7;
        if (runId < 3503)
            pos = 2;   // 2579..3500 - Pos. 2
        else if (runId < 3514)
            pos = 1;   // 3503..3513 - Pos. 1
        else if (runId < 3589)
            pos = 2;   // 3514..3588 - Pos. 2
        else if (runId < 4153)
            pos = 3;   // 3589..4152 - Pos. 3
        else if (runId < 4307)
            pos = 4;   // 4153..4305 - Pos. 4
        else if (runId < 4906)
            pos = 5;   // 4307..4905 - Pos. 5
        else if (runId < 4976)
            pos = 4;   // 4906..4972 - Pos. 4
        else
            pos = 5;   // 4976..5184 - Pos. 5
    }
 
    if (run && pos) {
        TString map = TString::Format("%s/input/ECAL_map_period_%d.txt", dir, run);
        TString cal = TString::Format("%s/parameters/ecal/ECAL_calibration_run%d.txt", dir, run);
        TString geo = TString::Format("%s/geometry/ECAL_v3_run%d_pos%d.root", dir, run, pos);
        LoadAdcMap(map.Data());
        LoadCalibration(cal.Data());
        LoadGeometry(geo.Data());
    } else {
        Error(__func__, "Unable to determine ECAL position for runId %lu\n", runId);
    }
}
 
BmnECALRaw2Digit::BmnECALRaw2Digit(const char* map_file_name, const char* calibr_file_name, const char* geo_file_name)
{
    LoadAdcMap(map_file_name);
    LoadCalibration(calibr_file_name);
    LoadGeometry(geo_file_name);
}
 
void BmnECALRaw2Digit::print()
{
    BmnECALMapElement::PrintTitle();
    for (Int_t chan = 1; chan < MAX_ECAL_CHANNELS + 1; chan++) {
        if (fECALMapArray[chan].GetChan() != -1) {
            fECALMapArray[chan].Print();
        }
    }
    printf("Pedestal: ");
    fPedestal.Print();
    printf("Signal: ");
    fSignal.Print();
}
 
Bool_t BmnECALRaw2Digit::LoadAdcMap(const char* map_file_name)
{
    ifstream in;
    in.open(map_file_name);
    if (!in.is_open()) {
        Fatal(__func__, "Loading ECAL Map from file: %s - file open error!\n", map_file_name);
        return kFALSE;
    }
    Info(__func__, "Loading ECAL Map from file: %s\n", map_file_name);
    TString dummy;
    in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy;
    while (!in.eof()) {
        ULong_t adcId;
        UShort_t adcChan;
        Int_t chan;
        Float_t x, y;
        in >> std::hex >> adcId >> std::dec >> adcChan >> chan >> x >> y;
        if (!in.good())
            break;
        if (adcId == 0 || adcChan == 0 || chan <= 0)
            continue;
        fECALMapArray[chan].SetAdcMap(chan, adcId, adcChan);
    }
    in.close();
    return kTRUE;
}
 
Bool_t BmnECALRaw2Digit::LoadCalibration(const char* calibr_file_name)
{
    ifstream in;
    in.open(calibr_file_name);
    if (!in.is_open()) {
        Fatal(__func__, "Loading ECAL Calibration from file: %s - file open error!\n", calibr_file_name);
        return kFALSE;
    }
    Info(__func__, "Loading ECAL Calibration from file: %s\n", calibr_file_name);
    TString dummy;
    in >> dummy >> dummy >> dummy >> dummy;
    while (!in.eof()) {
        int chan;
        float coeff, err;
        in >> chan >> coeff >> err;
        if (!in.good())
            break;
        if (chan <= 0)
            continue;
        fECALMapArray[chan].SetCoeff(chan, coeff);
    }
    in.close();
    return kTRUE;
}
 
Bool_t BmnECALRaw2Digit::LoadGeometry(const char* geo_file_name)
{
    TGeoVolume* top = TGeoVolume::Import(geo_file_name, "TOP");
    if (!top) {
        Fatal(__func__, "Volume TOP not found in %s\n", geo_file_name);
        return kFALSE;
    }
    TGeoNode* ecal = top->GetNode(0);
    if (!ecal) {
        Fatal(__func__, "Unexpected geometry %s\n", geo_file_name);
        return kFALSE;
    }
 
    TGeoNode* ecal1 = ecal->GetDaughter(0);
    if (ecal1) {
        Int_t n = ecal1->GetNdaughters();
        if (n > 504) {
            Fatal(__func__, "Expected ecal node 1 with 504 daughters or less in %s\n", geo_file_name);
            return kFALSE;
        }
        for (Int_t i = 0; i < n; i++) {
            Int_t ch = ecal1->GetDaughter(i)->GetNumber();
            if (ch < 1 || ch > 504) {
                Fatal(__func__, "Unexpected chan=%d at ecal node 1 in %s\n", ch, geo_file_name);
                return kFALSE;
            }
            fECALMapArray[ch].SetCoords(ecal1->GetDaughter(i), ecal1);
        }
    }
 
    TGeoNode* ecal2 = ecal->GetDaughter(1);
    if (ecal2) {
        Int_t n = ecal2->GetNdaughters();
        if (n > 504) {
            Fatal(__func__, "Expected ecal node 2 with 504 daughters or less in %s\n", geo_file_name);
            return kFALSE;
        }
        for (Int_t i = 0; i < n; i++) {
            Int_t ch = ecal2->GetDaughter(i)->GetNumber();
            if (ch < 505 || ch > 1008) {
                Fatal(__func__, "Unexpected chan=%d at ecal node 2 in %s\n", ch, geo_file_name);
                return kFALSE;
            }
            fECALMapArray[ch].SetCoords(ecal2->GetDaughter(i), ecal2);
        }
    }
 
    return kTRUE;
}
 
void BmnECALRaw2Digit::FillEvent(TClonesArray* data, TClonesArray* ecaldigit)
{
    Int_t N = data->GetEntriesFast();
    for (Int_t i = 0; i < N; i++) {
        BmnADCDigit<Short_t>* digit = (BmnADCDigit<Short_t>*)data->At(i);
        Int_t chan;
        for (chan = 1; chan < MAX_ECAL_CHANNELS + 1; chan++) {
            if (fECALMapArray[chan].GetChan() != chan)
                continue;
            if (fECALMapArray[chan].IsAdcMatch(digit->GetSerial(), digit->GetChannel()))
                break;
        }
        if (chan > MAX_ECAL_CHANNELS)
            continue;
        UseNovosibirskFit(digit, chan, ecaldigit);
    }
}
 
Bool_t BmnECALRaw2Digit::UseNovosibirskFit(BmnADCDigit<Short_t>* digit, Int_t chan, TClonesArray* ecaldigit)
{
    UInt_t N = digit->GetNSamples();
    //    Short_t* s = (Short_t*)digit->GetValue();
    auto& s = digit->GetValue();
 
    Float_t pedestal = 0.;
    for (Int_t i = fPedestal.sample0; i < fPedestal.sample0 + fPedestal.samples; i++) {
        if (i >= (Int_t)N)
            return 0;
        pedestal += (Float_t)s[i];
    }
 
    pedestal /= fPedestal.samples;
    if (pedestal < fPedestal.min || pedestal > fPedestal.max)
        return 0;
 
    TGraph graph(N);
    Int_t min_n = 0;
    Int_t min_cnt = 0;
    Short_t min_v = 0;
    Float_t scale = fECALMapArray[chan].GetCoeff();
    for (Int_t n = 0; n < (Int_t)N; n++) {
        Short_t& v = s[n];
        if (min_v > v) {
            min_v = v;
            min_n = n;
            min_cnt = 1;
        } else if (min_v == v) {
            min_n = n;
            min_cnt++;
        }
 
        Float_t y = v;
        y = (y - pedestal) * scale;
        graph.SetPoint(n, n, y);
    }
 
    if (min_v > fSignal.max || min_n - min_cnt < fSignal.sample0 || min_n > fSignal.sample0 + fSignal.samples) {
        return 0;
    }
 
    // Initialize fitting function
    TF1 pFit("Novosibirsk", Novosibirsk, 0., N, 4);
    pFit.SetParameters(min_v, 0.01, 5., min_n);
 
    // Remove points at the flat peak if any
    while (min_cnt-- > 0) {
        graph.RemovePoint(min_n--);
    }
 
    graph.Fit("Novosibirsk", "QR0", "");
 
    if (!gMinuit->fCstatu.EqualTo("CONVERGED "))
        return kFALSE;
    TF1* pRes = graph.GetFunction("Novosibirsk");
 
    Float_t Tp = pRes->GetMinimumX();
    if (Tp < fSignal.sample0)
        return kFALSE;
    if (Tp > fSignal.sample0 + fSignal.samples)
        return kFALSE;
 
    Float_t Vp = pRes->GetMinimum();
    // if (Vp > fSignal.max || Vp < fSignal.min) return kFALSE;
 
    Float_t To = pRes->GetX(0.5 * Vp, 0, Tp);
 
    Float_t dVdT = pRes->Derivative(To);
    if (dVdT >= 0)
        return kFALSE;
 
    To = To - 0.5 * Vp / dVdT;
 
    if (To < fPedestal.sample0)
        return kFALSE;
    if (To > fSignal.sample0 + fSignal.samples)
        return kFALSE;
 
    Float_t V = pRes->Integral(To, To + fSignal.samples) / fSignal.samples;
 
    BmnECALDigit* ecal = (BmnECALDigit*)ecaldigit->ConstructedAt(ecaldigit->GetEntriesFast());
    ecal->Set(fECALMapArray + chan, -V, -Vp, To, Tp);
 
    return kTRUE;
}
 
Double_t BmnECALRaw2Digit::Novosibirsk(Double_t* x, Double_t* p)
{
    // As Defined in RooNovosibirsk.cxx
    // If tail=eta=0 the Belle distribution becomes gaussia
    double tail = p[1];
    double width = p[2];
    double peak = p[3];
    if (TMath::Abs(tail) < 1.e-7)
        return p[0] * TMath::Exp(-0.5 * TMath::Power(((x[0] - peak) / width), 2));
    Double_t arg = 1.0 - (x[0] - peak) * tail / width;
    if (arg < 1.e-6)
        return 0.0;   // Argument of logaritem negative. Real continuation -> function equals zero
 
    Double_t log = TMath::Log(arg);
    static const Double_t xi = 2.3548200450309494;   // 2 Sqrt( Ln(4) )
    Double_t width_zero = (2.0 / xi) * TMath::ASinH(tail * xi * 0.5);
    Double_t width_zero2 = width_zero * width_zero;
    Double_t exponent = (-0.5 / (width_zero2)*log * log) - (width_zero2 * 0.5);
 
    return p[0] * TMath::Exp(exponent);
}