Commit 0fb6159c authored by Paul Schütze's avatar Paul Schütze
Browse files

Merge branch 'seedcharge' into 'master' 

DetectorHistogrammer: Correctly Calculate Cluster Seed Charge

See merge request allpix-squared/allpix-squared!545
parents 8816b8eb 1328a0de

src/modules/DetectorHistogrammer/Cluster.cpp

+24 −18
Original line number Diff line number Diff line 
/**

 * @file

 * @brief Implementation of object with a cluster of PixelHits

 * @copyright Copyright (c) 2017-2020 CERN and the Allpix Squared authors.

 * @copyright Copyright (c) 2017-2021 CERN and the Allpix Squared authors.

 * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".

 * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an

 * Intergovernmental Organization or submit itself to any jurisdiction.
@@ -10,6 +10,7 @@ 
#include "Cluster.hpp"



#include <algorithm>

#include <cmath>

#include <set>



#include "core/utils/log.h"
@@ -17,32 +18,38 @@ 


using namespace allpix;



Cluster::Cluster(const PixelHit* seedPixelHit) : seedPixelHit_(seedPixelHit) {

    pixelHits_.insert(seedPixelHit);

    clusterCharge_ = seedPixelHit->getSignal();

Cluster::Cluster(const PixelHit* seed_pixel_hit) : seed_pixel_hit_(seed_pixel_hit) {

    pixel_hits_.insert(seed_pixel_hit);

    cluster_charge_ = seed_pixel_hit->getSignal();



    minX_ = seedPixelHit->getPixel().getIndex().x();

    minX_ = seed_pixel_hit->getPixel().getIndex().x();

    maxX_ = minX_;

    minY_ = seedPixelHit->getPixel().getIndex().y();

    minY_ = seed_pixel_hit->getPixel().getIndex().y();

    maxY_ = minY_;



    for(const auto* mc_particle : seedPixelHit->getMCParticles()) {

    for(const auto* mc_particle : seed_pixel_hit->getMCParticles()) {

        mc_particles_.insert(mc_particle);

    }

}



bool Cluster::addPixelHit(const PixelHit* pixelHit) {

    auto ret = pixelHits_.insert(pixelHit);

bool Cluster::addPixelHit(const PixelHit* pixel_hit) {

    auto ret = pixel_hits_.insert(pixel_hit);

    if(ret.second == true) {

        clusterCharge_ += pixelHit->getSignal();

        unsigned int pixX = pixelHit->getPixel().getIndex().x();

        unsigned int pixY = pixelHit->getPixel().getIndex().y();

        cluster_charge_ += pixel_hit->getSignal();

        unsigned int pixX = pixel_hit->getPixel().getIndex().x();

        unsigned int pixY = pixel_hit->getPixel().getIndex().y();

        minX_ = std::min(pixX, minX_);

        maxX_ = std::max(pixX, maxX_);

        minY_ = std::min(pixY, minY_);

        maxY_ = std::max(pixY, maxY_);



        for(const auto* mc_particle : pixelHit->getMCParticles()) {

        // Update seed pixel if new charge is larger:

        if(std::signbit(seed_pixel_hit_->getSignal()) == std::signbit(pixel_hit->getSignal()) &&

           std::abs(seed_pixel_hit_->getSignal()) < std::abs(pixel_hit->getSignal())) {

            seed_pixel_hit_ = pixel_hit;

        }



        for(const auto* mc_particle : pixel_hit->getMCParticles()) {

            mc_particles_.insert(mc_particle);

        }
@@ -51,14 +58,13 @@ bool Cluster::addPixelHit(const PixelHit* pixelHit) { 
    return false;

}



ROOT::Math::XYVectorD Cluster::getPosition() const {

    ROOT::Math::XYVectorD meanPos = ROOT::Math::XYVectorD(0., 0.);

ROOT::Math::XYZPoint Cluster::getPosition() const {

    ROOT::Math::XYZVector meanPos;

    for(const auto& pixel : this->getPixelHits()) {

        meanPos += ROOT::Math::XYVectorD(pixel->getPixel().getIndex().x() * pixel->getSignal(),

                                         pixel->getPixel().getIndex().y() * pixel->getSignal());

        meanPos = pixel->getPixel().getLocalCenter() * pixel->getSignal() + meanPos;

    }

    meanPos /= getCharge();

    return meanPos;

    return static_cast<ROOT::Math::XYZPoint>(meanPos);

}



std::pair<unsigned int, unsigned int> Cluster::getSizeXY() const {

src/modules/DetectorHistogrammer/Cluster.hpp

+19 −19
Original line number Diff line number Diff line 
/**

 * @file

 * @brief Definition of object with a cluster containing several PixelHits

 * @copyright Copyright (c) 2017-2020 CERN and the Allpix Squared authors.

 * @copyright Copyright (c) 2017-2021 CERN and the Allpix Squared authors.

 * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md".

 * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an

 * Intergovernmental Organization or submit itself to any jurisdiction.
@@ -10,7 +10,7 @@ 
#ifndef ALLPIX_DETECTOR_HISTOGRAMMER_CLUSTER_H

#define ALLPIX_DETECTOR_HISTOGRAMMER_CLUSTER_H



#include <Math/DisplacementVector2D.h>

#include <Math/Vector3D.h>



#include <set>
@@ -23,27 +23,27 @@ namespace allpix { 
    public:

        /**

         * @brief Construct a cluster

         * @param seedPixelHit PixelHit to start the cluster with

         * @param seed_pixel_hit PixelHit to start the cluster with

         */

        explicit Cluster(const PixelHit* seedPixelHit);

        explicit Cluster(const PixelHit* seed_pixel_hit);



        /**

         * @brief Get the signal data for the hit

         * @return Digitized signal

         * @brief Get the total accumulated signal of the cluster

         * @return Total cluster charge

         */

        double getCharge() const { return clusterCharge_; }

        double getCharge() const { return cluster_charge_; }



        /**

         * @brief Add PixelHit to the cluster

         * @param pixelHit PixelHit to be added

         * @brief Add a PixelHit to the cluster

         * @param pixel_hit PixelHit to be added

         */

        bool addPixelHit(const PixelHit* pixelHit);

        bool addPixelHit(const PixelHit* pixel_hit);



        /**

         * @brief Get the cluster size

         * @return cluster size

         */

        unsigned long int getSize() const { return pixelHits_.size(); }

        unsigned long int getSize() const { return pixel_hits_.size(); }



        /**

         * @brief Get the cluster sizes in x and y
@@ -52,16 +52,16 @@ namespace allpix { 
        std::pair<unsigned int, unsigned int> getSizeXY() const;



        /**

         * @brief Get the weighted mean cluster position

         * @brief Get the charge-weighted mean cluster position in local coordinates

         * @return weighted mean cluster position

         */

        ROOT::Math::XYVectorD getPosition() const;

        ROOT::Math::XYZPoint getPosition() const;



        /**

         * @brief Get the seed PixelHit

         * @brief Get the seed PixelHit, i.e. the PixelHit with the largest charge

         * @return Pointer to PixelHit

         */

        const PixelHit* getSeedPixelHit() const { return seedPixelHit_; }

        const PixelHit* getSeedPixelHit() const { return seed_pixel_hit_; }



        /**

         * @brief Get the PixelHit at coordinates x and y
@@ -75,7 +75,7 @@ namespace allpix { 
         * @brief Get the PixelHits contained in this cluster

         * @return List of all contained PixelHits

         */

        std::set<const PixelHit*> getPixelHits() const { return pixelHits_; }

        std::set<const PixelHit*> getPixelHits() const { return pixel_hits_; }



        /**

         * @brief Get all MCParticles related to the cluster
@@ -84,12 +84,12 @@ namespace allpix { 
        std::set<const MCParticle*> getMCParticles() const;



    private:

        const PixelHit* seedPixelHit_;

        const PixelHit* seed_pixel_hit_;



        std::set<const PixelHit*> pixelHits_;

        std::set<const PixelHit*> pixel_hits_;

        std::set<const MCParticle*> mc_particles_;



        double clusterCharge_{};

        double cluster_charge_{};



        unsigned int minX_, minY_, maxX_, maxY_;

    };

src/modules/DetectorHistogrammer/DetectorHistogrammerModule.cpp

+21 −19
Original line number Diff line number Diff line
@@ -289,6 +289,10 @@ void DetectorHistogrammerModule::initialize() { 
    cluster_charge = CreateHistogram<TH1D>(

        "cluster_charge", cluster_charge_title.c_str(), 1000, 0., static_cast<double>(max_cluster_charge));



    std::string cluster_seed_charge_title = "Cluster seed charge for " + detector_->getName() + ";seed charge [ke];clusters";

    cluster_seed_charge = CreateHistogram<TH1D>(

        "cluster_seed_charge", cluster_seed_charge_title.c_str(), 1000, 0., static_cast<double>(max_cluster_charge));



    std::string pixel_charge_title = "Pixel charge for " + detector_->getName() + ";pixel charge [ke];pixels";

    pixel_charge =

        CreateHistogram<TH1D>("pixel_charge", pixel_charge_title.c_str(), 1000, 0., static_cast<double>(max_cluster_charge));
@@ -348,8 +352,10 @@ void DetectorHistogrammerModule::run(Event* event) { 
        cluster_size_y->Fill(clusSizesXY.second);



        auto clusterPos = clus.getPosition();

        LOG(DEBUG) << "Cluster at coordinates " << clusterPos << " with charge " << Units::display(clus.getCharge(), "ke");

        cluster_map->Fill(clusterPos.x(), clusterPos.y());

        auto [cluster_x, cluster_y] = detector_->getModel()->getPixelIndex(clusterPos);

        LOG(DEBUG) << "Cluster at indices " << cluster_x << ", " << cluster_y << "(" << clusterPos

                   << " local coordinates) with charge " << Units::display(clus.getCharge(), "ke");

        cluster_map->Fill(cluster_x, cluster_y);

        cluster_charge->Fill(static_cast<double>(Units::convert(clus.getCharge(), "ke")));

        charge_sum += clus.getCharge();
@@ -366,8 +372,6 @@ void DetectorHistogrammerModule::run(Event* event) { 


        LOG(TRACE) << "Matching primaries: " << intersection.size();

        for(const auto& particle : intersection) {

            auto pitch = detector_->getModel()->getPixelSize();



            auto particlePos = particle->getLocalReferencePoint() + track_smearing(track_resolution_);

            LOG(DEBUG) << "MCParticle at " << Units::display(particlePos, {"mm", "um"});
@@ -387,16 +391,15 @@ void DetectorHistogrammerModule::run(Event* event) { 
            cluster_charge_map->Fill(

                inPixel_um_x, inPixel_um_y, static_cast<double>(Units::convert(clus.getCharge(), "ke")));



            // Retrieve the pixel to which this MCParticle points:

            const auto* pixel = clus.getPixelHit(static_cast<unsigned int>(xpixel), static_cast<unsigned int>(ypixel));

            if(pixel != nullptr) {

            // Retrieve the seed pixel:

            const auto* seed_pixel = clus.getSeedPixelHit();

            seed_charge_map->Fill(

                    inPixel_um_x, inPixel_um_y, static_cast<double>(Units::convert(pixel->getSignal(), "ke")));

            }

                inPixel_um_x, inPixel_um_y, static_cast<double>(Units::convert(seed_pixel->getSignal(), "ke")));

            cluster_seed_charge->Fill(static_cast<double>(Units::convert(seed_pixel->getSignal(), "ke")));



            // Calculate residual with cluster position:

            auto residual_um_x = static_cast<double>(Units::convert(particlePos.x() - clusterPos.x() * pitch.x(), "um"));

            auto residual_um_y = static_cast<double>(Units::convert(particlePos.y() - clusterPos.y() * pitch.y(), "um"));

            auto residual_um_x = static_cast<double>(Units::convert(particlePos.x() - clusterPos.x(), "um"));

            auto residual_um_y = static_cast<double>(Units::convert(particlePos.y() - clusterPos.y(), "um"));

            residual_x->Fill(residual_um_x);

            residual_y->Fill(residual_um_y);

            residual_x_vs_x->Fill(inPixel_um_x, std::fabs(residual_um_x));
@@ -419,8 +422,6 @@ void DetectorHistogrammerModule::run(Event* event) { 


    // Calculate efficiency: search for matching clusters for all primary MCParticles

    for(auto& particle : primary_particles) {

        auto pitch = detector_->getModel()->getPixelSize();



        // Calculate 2D local position of particle:

        auto particlePos = particle->getLocalReferencePoint() + track_smearing(track_resolution_);
@@ -432,10 +433,9 @@ void DetectorHistogrammerModule::run(Event* event) { 
        auto inPixel_um_x = static_cast<double>(Units::convert(inPixelPos.x(), "um"));

        auto inPixel_um_y = static_cast<double>(Units::convert(inPixelPos.y(), "um"));



        auto matched_cluster =

            std::find_if(clusters.begin(), clusters.end(), [this, &particlePos, &pitch](const Cluster& clus) {

                return (std::fabs(clus.getPosition().x() * pitch.x() - particlePos.x()) < matching_cut_.x()) &&

                       (std::fabs(clus.getPosition().y() * pitch.y() - particlePos.y()) < matching_cut_.y());

        auto matched_cluster = std::find_if(clusters.begin(), clusters.end(), [this, &particlePos](const Cluster& clus) {

            return (std::fabs(clus.getPosition().x() - particlePos.x()) < matching_cut_.x()) &&

                   (std::fabs(clus.getPosition().y() - particlePos.y()) < matching_cut_.y());

        });



        // Do we have a match?
@@ -489,6 +489,7 @@ void DetectorHistogrammerModule::finalize() { 
    auto efficiency_map_histogram = efficiency_map->Merge();

    auto n_cluster_histogram = n_cluster->Merge();

    auto cluster_charge_histogram = cluster_charge->Merge();

    auto cluster_seed_charge_histogram = cluster_seed_charge->Merge();

    auto cluster_charge_map_histogram = cluster_charge_map->Merge();

    auto seed_charge_map_histogram = seed_charge_map->Merge();

    auto pixel_charge_histogram = pixel_charge->Merge();
@@ -604,6 +605,7 @@ void DetectorHistogrammerModule::finalize() { 
    efficiency_map_histogram->Write();

    n_cluster_histogram->Write();

    cluster_charge_histogram->Write();

    cluster_seed_charge_histogram->Write();

    pixel_charge_histogram->Write();

    cluster_charge_map_histogram->Write();

    seed_charge_map_histogram->Write();

src/modules/DetectorHistogrammer/DetectorHistogrammerModule.hpp

+1 −1
Original line number Diff line number Diff line
@@ -100,7 +100,7 @@ namespace allpix { 
        Histogram<TH1D> event_size;

        Histogram<TH1D> cluster_size, cluster_size_x, cluster_size_y;

        Histogram<TH1D> n_cluster;

        Histogram<TH1D> cluster_charge, pixel_charge, total_charge;

        Histogram<TH1D> cluster_charge, pixel_charge, total_charge, cluster_seed_charge;

    };

} // namespace allpix