BmnSiliconLayer.cxx

Go to the documentation of this file.

#include "BmnSiliconLayer.h"
 
#include "TMath.h"
 
#include <iostream>
 
using namespace TMath;
using namespace std;
 
Int_t BmnSiliconLayer::fUniqueIdL = 0;
Int_t BmnSiliconLayer::fUniqueIdU = 0;
 
BmnSiliconLayer::BmnSiliconLayer()
{
 
    Verbosity = true;
 
    ZoneID = -1;
    LayerType = LowerStripLayer;
    FirstStripNumber = 0;
 
    XMinLayer = 0.0;
    XMaxLayer = 10.0;
    YMinLayer = 0.0;
    YMaxLayer = 10.0;
 
    Pitch = 0.01;                       // cm
    AngleDeg = 2.5;                     // in degrees from a vertical line where a plus value is clockwise
    AngleRad = AngleDeg * Pi() / 180;   // in radians
 
    StripOrder = LeftToRight;
 
    if (AngleDeg <= 0.0 && AngleDeg >= -90.0) {
        SetStripNumberingBorders(LeftBottom, RightTop);
    } else {
        if (AngleDeg > 0.0 && AngleDeg <= 90.0) {
            SetStripNumberingBorders(LeftTop, RightBottom);
        }
    }
 
    ClusterFindingThreshold = 0.0;
 
    InitializeLayer();
}
 
BmnSiliconLayer::BmnSiliconLayer(Int_t zone_id,
                                 StripLayerType layer_type,
                                 Int_t first_strip_number,
                                 Double_t xsize,
                                 Double_t ysize,
                                 Double_t xorig,
                                 Double_t yorig,
                                 Double_t pitch,
                                 Double_t adeg)
{
 
    Verbosity = true;
 
    ZoneID = zone_id;
    LayerType = layer_type;
    FirstStripNumber = first_strip_number;
 
    XMinLayer = xorig;
    XMaxLayer = xorig + xsize;
    YMinLayer = yorig;
    YMaxLayer = yorig + ysize;
 
    Pitch = pitch;
    AngleDeg = adeg;                    // in degrees from a vertical line where a plus value is clockwise
    AngleRad = AngleDeg * Pi() / 180;   // in radians
 
    StripOrder = LeftToRight;
 
    if (AngleDeg <= 0.0 && AngleDeg >= -90.0) {
        SetStripNumberingBorders(LeftBottom, RightTop);
    } else {
        if (AngleDeg > 0.0 && AngleDeg <= 90.0) {
            SetStripNumberingBorders(LeftTop, RightBottom);
        }
    }
 
    ClusterFindingThreshold = 0.0;
 
    InitializeLayer();
}
 
BmnSiliconLayer::~BmnSiliconLayer() {}
 
void BmnSiliconLayer::InitializeLayer()
{
    Int_t NStrips = GetNStrips();
 
    // strips
    Strips.clear();
    Strips.resize(NStrips, 0.0);
 
    // matches
    ResetStripMatches();
    ResetStripDigitNumberMatches();
 
    // strip hits
    ResetStripHits();
 
    BmnSiliconLayer::fUniqueIdL = 0;
    BmnSiliconLayer::fUniqueIdU = 0;
}
 
void BmnSiliconLayer::ResetLayer()
{
    InitializeLayer();
}
 
Int_t BmnSiliconLayer::GetNStrips()
{
    Double_t xleft_layer = ConvertNormalPointToStripX(XLeftPointOfStripNumbering, YLeftPointOfStripNumbering);
    Double_t xright_layer = ConvertNormalPointToStripX(XRightPointOfStripNumbering, YRightPointOfStripNumbering);
    Double_t n_strips = (xright_layer - xleft_layer) / Pitch;
 
    if ((n_strips - (Int_t)n_strips) < 1E-10) {
        return (Int_t)n_strips;
    } else {
        return (Int_t)(n_strips + 1);
    }
}
 
void BmnSiliconLayer::SetPitch(Double_t pitch)
{
    Pitch = pitch;
 
    InitializeLayer();
}
 
void BmnSiliconLayer::SetLayerSizes(Double_t xsize, Double_t ysize, Double_t xorig, Double_t yorig)
{
    XMinLayer = xorig;
    XMaxLayer = xorig + xsize;
    YMinLayer = yorig;
    YMaxLayer = yorig + ysize;
 
    InitializeLayer();
}
 
void BmnSiliconLayer::SetAngleDeg(Double_t deg)
{
    if (Abs(deg) <= 90.0) {
        AngleDeg = deg;
    } else {
        AngleDeg = 0.0;
        if (Verbosity)
            cerr << "WARNING: BmnSiliconLayer::SetAngleDeg(): the strip angle is incorrect and has been set to zero\n";
    }
 
    AngleRad = AngleDeg * Pi() / 180;
 
    InitializeLayer();
}
 
Bool_t BmnSiliconLayer::SetStripNumberingBorders(Double_t x_left, Double_t y_left, Double_t x_right, Double_t y_right)
{
    if (x_right < x_left)
        return false;
 
    XLeftPointOfStripNumbering = x_left;
    YLeftPointOfStripNumbering = y_left;
    XRightPointOfStripNumbering = x_right;
    YRightPointOfStripNumbering = y_right;
 
    ResetLayer();
 
    return true;
}
 
Bool_t BmnSiliconLayer::SetStripNumberingBorders(StripBorderPoint left, StripBorderPoint right)
{
 
    if (left != LeftTop && left != LeftBottom) {
        if (Verbosity)
            cerr << "WARNING: SetStripNumberingBorders: left strip border point is incorrect\n";
        return false;
    }
    if (right != RightTop && right != RightBottom) {
        if (Verbosity)
            cerr << "WARNING: SetStripNumberingBorders: right strip border point is incorrect\n";
        return false;
    }
 
    XLeftPointOfStripNumbering = XMinLayer;
    XRightPointOfStripNumbering = XMaxLayer;
 
    switch (left) {
        case LeftTop:
            YLeftPointOfStripNumbering = YMaxLayer;
            break;
        case LeftBottom:
            YLeftPointOfStripNumbering = YMinLayer;
            break;
        default:
            return false;
    }
 
    switch (right) {
        case RightTop:
            YRightPointOfStripNumbering = YMaxLayer;
            break;
        case RightBottom:
            YRightPointOfStripNumbering = YMinLayer;
            break;
        default:
            return false;
    }
 
    ResetLayer();
 
    return true;
}
 
Bool_t BmnSiliconLayer::SetStripNumberingOrder(StripNumberingDirection strip_direction)
{
    StripOrder = strip_direction;
 
    ResetLayer();
 
    return true;
}
 
Bool_t BmnSiliconLayer::AddDeadZone(Int_t n_points, Double_t* x_points, Double_t* y_points)
{
    DeadZoneOfStripLayer dead_zone;
    Bool_t status = dead_zone.SetDeadZone(n_points, x_points, y_points);
    if (status == true) {
        DeadZones.push_back(dead_zone);
        return true;
    } else {
        return false;
    }
}
 
Bool_t BmnSiliconLayer::AddDeadZone(DeadZoneOfStripLayer dead_zone)
{
    if (dead_zone.GetNPoints() > 2) {
        DeadZones.push_back(dead_zone);
        return true;
    } else {
        return false;
    }
}
 
Bool_t BmnSiliconLayer::IsPointInsideDeadZones(Double_t x, Double_t y)
{
    for (size_t izone = 0; izone < DeadZones.size(); ++izone)
        if (DeadZones[izone].IsInside(x, y))
            return true;
 
    return false;
}
 
Bool_t BmnSiliconLayer::IsPointInsideStripLayer(Double_t x, Double_t y)
{
    if (x >= XMinLayer && x <= XMaxLayer && y >= YMinLayer && y <= YMaxLayer && !IsPointInsideDeadZones(x, y))
        return true;
    else
        return false;
}
 
Bool_t BmnSiliconLayer::SetStripSignal(Int_t strip_num, Double_t signal)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)Strips.size() + FirstStripNumber)) {
        Strips[strip_num - FirstStripNumber] = signal;
        return true;
    } else
        return false;
}
 
Bool_t BmnSiliconLayer::AddStripSignal(Int_t strip_num, Double_t signal)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)Strips.size() + FirstStripNumber)) {
        Strips[strip_num - FirstStripNumber] += signal;
        return true;
    } else
        return false;
}
 
Double_t BmnSiliconLayer::GetStripSignal(Int_t strip_num)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)Strips.size() + FirstStripNumber))
        return Strips[strip_num - FirstStripNumber];
    else
        return -1;
}
 
Bool_t BmnSiliconLayer::SetStripMatch(Int_t strip_num, BmnMatch mc_match)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)StripMatches.size() + FirstStripNumber)) {
        StripMatches[strip_num - FirstStripNumber] = mc_match;
        return true;
    } else
        return false;
}
 
Bool_t BmnSiliconLayer::AddLinkToStripMatch(Int_t strip_num, Double_t weight, Int_t mc_num)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)StripMatches.size() + FirstStripNumber)) {
        StripMatches[strip_num - FirstStripNumber].AddLink(weight, mc_num);
        return true;
    } else
        return false;
}
 
Bool_t BmnSiliconLayer::SetStripDigitNumberMatch(Int_t strip_num, BmnMatch digit_num_match)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)StripDigitNumberMatches.size() + FirstStripNumber)) {
        StripDigitNumberMatches[strip_num - FirstStripNumber] = digit_num_match;
        return true;
    } else
        return false;
}
 
Bool_t BmnSiliconLayer::AddLinkToStripDigitNumberMatch(Int_t strip_num, Double_t weight, Int_t digit_num)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)StripDigitNumberMatches.size() + FirstStripNumber)) {
        StripDigitNumberMatches[strip_num - FirstStripNumber].AddLink(weight, digit_num);
        return true;
    } else
        return false;
}
 
BmnMatch BmnSiliconLayer::GetStripMatch(Int_t strip_num)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)StripMatches.size() + FirstStripNumber))
        return StripMatches[strip_num - FirstStripNumber];
    else
        return BmnMatch();   // return an empty match
}
 
BmnMatch BmnSiliconLayer::GetStripDigitNumberMatch(Int_t strip_num)
{
    if (strip_num >= FirstStripNumber && strip_num < ((Int_t)StripDigitNumberMatches.size() + FirstStripNumber))
        return StripDigitNumberMatches[strip_num - FirstStripNumber];
    else
        return BmnMatch();   // return an empty match
}
 
void BmnSiliconLayer::ResetStripMatches()
{
    Int_t NStrips = GetNStrips();
 
    StripMatches.clear();
    StripMatches.resize(NStrips, BmnMatch());
}
 
void BmnSiliconLayer::ResetStripDigitNumberMatches()
{
    Int_t NStrips = GetNStrips();
 
    StripDigitNumberMatches.clear();
    StripDigitNumberMatches.resize(NStrips, BmnMatch());
}
 
Double_t BmnSiliconLayer::GetStripHitPos(Int_t num)
{
    if (num >= 0 && num < (Int_t)StripHits.size())
        return StripHits[num];
 
    return -1.0;
}
 
Double_t BmnSiliconLayer::GetStripHitTotalSignal(Int_t num)
{
    if (num >= 0 && num < (Int_t)StripHitsTotalSignal.size())
        return StripHitsTotalSignal[num];
 
    return -1.0;
}
 
Double_t BmnSiliconLayer::GetStripHitError(Int_t num)
{
    if (num >= 0 && num < (Int_t)StripHitsErrors.size())
        return StripHitsErrors[num];
 
    return -1.0;
}
 
Int_t BmnSiliconLayer::GetStripHitClusterSize(Int_t num)
{
    if (num >= 0 && num < (Int_t)StripHitsClusterSize.size())
        return StripHitsClusterSize[num];
 
    return -1.0;
}
 
void BmnSiliconLayer::ResetStripHits()
{
    StripHits.clear();
    StripHitsTotalSignal.clear();
    StripHitsErrors.clear();
    StripHitsClusterSize.clear();
    StripClusters.clear();
}
 
Double_t BmnSiliconLayer::ConvertNormalPointToStripX(Double_t x, Double_t y)
{
    Double_t XRotationCenter;
    Double_t YRotationCenter;
 
    if (StripOrder == LeftToRight) {
        XRotationCenter = XLeftPointOfStripNumbering;
        YRotationCenter = YLeftPointOfStripNumbering;
    } else {
        if (StripOrder == RightToLeft) {
            XRotationCenter = XRightPointOfStripNumbering;
            YRotationCenter = YRightPointOfStripNumbering;
        }
    }
    Double_t StripX =
        (x - XRotationCenter) * Cos(-AngleRad) + (y - YRotationCenter) * Sin(-AngleRad) + XRotationCenter;   // see
    return StripX;
}
 
Double_t BmnSiliconLayer::ConvertNormalPointToStripY(Double_t x, Double_t y)
{
    Double_t XRotationCenter;
    Double_t YRotationCenter;
 
    if (StripOrder == LeftToRight) {
        XRotationCenter = XLeftPointOfStripNumbering;
        YRotationCenter = YLeftPointOfStripNumbering;
    } else {
        if (StripOrder == RightToLeft) {
            XRotationCenter = XRightPointOfStripNumbering;
            YRotationCenter = YRightPointOfStripNumbering;
        }
    }
    Double_t StripY =
        -(x - XRotationCenter) * Sin(-AngleRad) + (y - YRotationCenter) * Cos(-AngleRad) + YRotationCenter;   // see
    return StripY;
}
 
Double_t BmnSiliconLayer::ConvertPointToStripPosition(Double_t x, Double_t y)
{
    // This function returns real(double) value,
    // where integer part - number of a strip, fractional part - position on this strip (as ratio from begin)
 
    if (StripOrder == LeftToRight) {
        return (ConvertNormalPointToStripX(x, y) - XLeftPointOfStripNumbering) / Pitch + FirstStripNumber;
    } else {
        if (StripOrder == RightToLeft) {
            return (XRightPointOfStripNumbering - ConvertNormalPointToStripX(x, y)) / Pitch + FirstStripNumber;
        }
    }
 
    return -1;
}
 
Double_t BmnSiliconLayer::CalculateStripEquationB(Double_t strip_pos)
{
 
    Double_t b;
 
    if (Abs(AngleDeg) != 90.0) {   // case: y=a*x+b
        Double_t x_strip_shift = ((strip_pos - FirstStripNumber) * Pitch) / Cos(Abs(AngleRad));
        Double_t xcoord;
 
        if (StripOrder == LeftToRight) {
            xcoord = XLeftPointOfStripNumbering + x_strip_shift;
            b = YLeftPointOfStripNumbering - Tan(PiOver2() - AngleRad) * xcoord;
        } else {
            if (StripOrder == RightToLeft) {
                xcoord = XRightPointOfStripNumbering - x_strip_shift;
                b = YRightPointOfStripNumbering - Tan(PiOver2() - AngleRad) * xcoord;
            }
        }
    } else {   // case: y=b
        Double_t y_strip_shift = (strip_pos - FirstStripNumber) * Pitch * Sin(AngleRad);
 
        if (StripOrder == LeftToRight) {
            b = YLeftPointOfStripNumbering - y_strip_shift;
        } else {
            if (StripOrder == RightToLeft) {
                b = YRightPointOfStripNumbering + y_strip_shift;
            }
        }
    }
 
    return b;
}
 
void BmnSiliconLayer::FindClustersAndStripHits()
{
 
    ResetStripHits();
 
    Double_t threshold = ClusterFindingThreshold;
    // Double_t threshold = 500.0;
 
    StripCluster cluster;
 
    Bool_t ascent = false;
    Bool_t descent = false;
 
    // Processing strips
    vector<Double_t> AnalyzableStrips = Strips;
 
    // Smooth strip signal ------------------------------------------------------
    // if(Pitch > 0.079) SmoothStripSignal(AnalyzableStrips, 1, 1, 1.0);
    // else SmoothStripSignal(AnalyzableStrips, 2, 1, 1.0);
    //--------------------------------------------------------------------------
 
    for (Int_t is = 0; is < (Int_t)AnalyzableStrips.size(); is++) {
        if (AnalyzableStrips.at(is) <= threshold) {
            if (descent || ascent) {
                descent = false;
                ascent = false;
                // make strip hit
                MakeStripHit(cluster, AnalyzableStrips, is);
            }
            continue;
        }
 
        if (cluster.GetClusterSize() > 0) {
            if (AnalyzableStrips.at(is) >= (cluster.Signals.at(cluster.GetClusterSize() - 1))) {
                if (descent) {
                    ascent = false;
                    descent = false;
                    // make strip hit
                    MakeStripHit(cluster, AnalyzableStrips, is);
                    // continue;
                }
                ascent = true;
                descent = false;
            } else {
                ascent = false;
                descent = true;
            }
        } else {
            ascent = true;
            descent = false;
        }
 
        cluster.AddStrip(is, AnalyzableStrips.at(is));
    }
 
    if (cluster.GetClusterSize() != 0) {
        // make strip hit
        Int_t lastnum = AnalyzableStrips.size() - 1;
        MakeStripHit(cluster, AnalyzableStrips, lastnum);
    }
}
 
void BmnSiliconLayer::MakeStripHit(StripCluster& cluster, vector<Double_t>& AnalyzableStrips, Int_t& curcnt)
{
 
    Double_t mean_strip_position = 0.0;
    Double_t total_cluster_signal = 0.0;
    Double_t cluster_rms = 0.0;
 
    Int_t NStripsInCluster = cluster.GetClusterSize();
 
    // Share the common strip signal of two adjacent clusters in half
    if (true) {
        if (AnalyzableStrips.at(curcnt) >= AnalyzableStrips.at(curcnt - 1)) {
            cluster.Signals.at(cluster.Signals.size() - 1) *= 0.5;
            AnalyzableStrips.at(curcnt - 1) *= 0.5;
        }
    }
 
    // mean position and cluster signal -----------------------------------------
    for (Int_t i = 0; i < NStripsInCluster; ++i) {
        Double_t strip_num = cluster.Strips.at(i) + FirstStripNumber;
        Double_t signal = cluster.Signals.at(i);
        total_cluster_signal += signal;
        mean_strip_position += (strip_num + 0.5) * signal;
    }
    mean_strip_position /= total_cluster_signal;
 
    if (mean_strip_position < FirstStripNumber)
        mean_strip_position = FirstStripNumber;
    if (mean_strip_position >= AnalyzableStrips.size() + FirstStripNumber)
        mean_strip_position = AnalyzableStrips.size() - 0.001;
    //--------------------------------------------------------------------------
 
    // cluster standard deviation (sigma): RMS ----------------------------------
    //    if(NStripsInCluster > 1) {
    //        for(Int_t i = 0; i < NStripsInCluster; ++i) {
    //            Double_t strip_num = cluster.Strips.at(i) + FirstStripNumber;
    //            Double_t signal = cluster.Signals.at(i);
    //            Double_t residual = (strip_num+0.5) - mean_strip_position;
    //            cluster_rms += residual*residual*signal;
    //        }
    //        cluster_rms /= total_cluster_signal;
    //        cluster_rms = TMath::Sqrt(cluster_rms);
    //    }
    //    else {
    //        cluster_rms = 1.0/TMath::Sqrt(12.0);
    //    }
    //
 
    // AZ, STS, method2
    Double_t sumW = total_cluster_signal;
    Double_t sumWX = 0.;
    Double_t sumWX2 = 0.;
    if (NStripsInCluster > 1) {
        for (Int_t i = 0; i < NStripsInCluster; ++i) {
            Double_t strip_num = cluster.Strips.at(i) + FirstStripNumber;
            Double_t signal = cluster.Signals.at(i);
            sumWX += strip_num * signal;
            sumWX2 += strip_num * strip_num * signal;
        }
        cluster_rms = (sumWX2 - sumWX * sumWX / sumW) / sumW;
    }
 
    else
    {
        cluster_rms = 1.0 / TMath::Sqrt(12.0);
    }
 
    //--------------------------------------------------------------------------
 
    StripHits.push_back(mean_strip_position);
    StripHitsTotalSignal.push_back(total_cluster_signal);
    StripHitsErrors.push_back(cluster_rms);
    StripHitsClusterSize.push_back(NStripsInCluster);
 
    cluster.MeanPosition = mean_strip_position;
    cluster.TotalSignal = total_cluster_signal;
    cluster.SetWidth(NStripsInCluster);
    cluster.Error = cluster_rms;
    if (LayerType == LowerStripLayer) {
        cluster.SetType(0);
        cluster.SetUniqueID(fUniqueIdL++);
    } else {
        cluster.SetType(1);
        cluster.SetUniqueID(fUniqueIdU++);
    }
 
    /* Update strip numbers in the local cluster according to
     * the FirstStripNumber shift before we can use it in the outer code
     */
    for (Int_t i = 0; i < NStripsInCluster; ++i) {
        cluster.Strips.at(i) += FirstStripNumber;
    }
 
    StripClusters.push_back(cluster);
 
    // return to a previous strip
    curcnt--;
    if (curcnt < 0)
        curcnt = 0;
 
    cluster.Clear();
}
 
void BmnSiliconLayer::SmoothStripSignal(vector<Double_t>& AnalyzableStrips,
                                        Int_t NIterations,
                                        Int_t SmoothWindow,
                                        Double_t Weight)
{
    // It's the Simple Moving Average (SMA) method
    // AnalyzableStrips - analyzable strip layer (ref)
    // NIterations - number of smooth iterations (usually 1)
    // SmoothWindow - number of strips on the left-right of the current strip (usually 1)
    // Weight - weight of the current strip (usually 1.0 - for simplicity, greater - for weighted value))
 
    SmoothStrips.clear();
    Int_t NStrips = AnalyzableStrips.size();
 
    for (Int_t iteration = 0; iteration < NIterations; ++iteration) {
        SmoothStrips.clear();
        for (Int_t istrip = 0; istrip < NStrips; ++istrip) {
            Double_t mean_value = 0.0;
            for (Int_t iw = istrip - SmoothWindow; iw <= istrip + SmoothWindow; ++iw) {
                if (iw >= 0 && iw < NStrips) {
                    if (iw == istrip)
                        mean_value += AnalyzableStrips[iw] * Weight;
                    else
                        mean_value += AnalyzableStrips[iw];
                }
            }
            mean_value /= 2.0 * SmoothWindow + Weight;
            SmoothStrips.push_back(mean_value);
        }
        AnalyzableStrips = SmoothStrips;
    }
 
    return;
}