/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoEcal/EgammaClusterAlgos/src/Multi5x5BremRecoveryClusterAlgo.cc

Go to the documentation of this file.

#include "RecoEcal/EgammaClusterAlgos/interface/Multi5x5BremRecoveryClusterAlgo.h"
#include "RecoEcal/EgammaCoreTools/interface/BremRecoveryPhiRoadAlgo.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

reco::SuperClusterCollection Multi5x5BremRecoveryClusterAlgo::makeSuperClusters(reco::CaloClusterPtrVector & clustersCollection)
{

        const float etaBorder = 1.479;
        superclusters_v.clear();

        // create vectors of references to clusters of a specific origin...
        reco::CaloClusterPtrVector islandClustersBarrel_v;
        reco::CaloClusterPtrVector islandClustersEndCap_v;

        // ...and populate them:
        for (reco::CaloCluster_iterator it = clustersCollection.begin(); it != clustersCollection.end(); it++)
        {
                reco::CaloClusterPtr cluster_p = *it;
                if (cluster_p->algo() == reco::CaloCluster::multi5x5) 
                {
                        if (fabs(cluster_p->position().eta()) < etaBorder)
                        {
                                islandClustersBarrel_v.push_back(cluster_p);
                        }
                        else
                        {
                                islandClustersEndCap_v.push_back(cluster_p);
                        }
                }
        }

        // make the superclusters from the Barrel clusters - Island
        makeIslandSuperClusters(islandClustersBarrel_v, eb_rdeta_, eb_rdphi_);
        // make the superclusters from the EndCap clusters - Island
        makeIslandSuperClusters(islandClustersEndCap_v, ec_rdeta_, ec_rdphi_);

        return superclusters_v;
}

#include "DataFormats/Math/interface/Vector3D.h"

void Multi5x5BremRecoveryClusterAlgo::makeIslandSuperClusters(reco::CaloClusterPtrVector &clusters_v, 
                double etaRoad, double phiRoad)
{

        std::vector<DetId> usedSeedDetIds;
        usedSeedDetIds.clear();

        for (reco::CaloCluster_iterator currentSeed = clusters_v.begin(); currentSeed != clusters_v.end(); ++currentSeed)
        {

                // check this seed was not already used
                if (std::find(usedSeedDetIds.begin(), usedSeedDetIds.end(), (*currentSeed)->seed()) 
                        != usedSeedDetIds.end()) continue;

                // Does our highest energy cluster have high enough energy?
                // changed this to continue from break (to be robust against the order of sorting of the seed clusters)
                if ((*currentSeed)->energy() * sin((*currentSeed)->position().theta()) < seedTransverseEnergyThreshold) continue;

                // if yes, make it a seed for a new SuperCluster:
                double energy = (*currentSeed)->energy();
                math::XYZVector position_((*currentSeed)->position().X(), 
                                (*currentSeed)->position().Y(), 
                                (*currentSeed)->position().Z());
                position_ *= energy;
                usedSeedDetIds.push_back((*currentSeed)->seed());

LogTrace("EcalClusters") << "*****************************";
LogTrace("EcalClusters") << "******NEW SUPERCLUSTER*******";
LogTrace("EcalClusters") << "Seed R = " << (*currentSeed)->position().Rho();

                reco::CaloClusterPtrVector constituentClusters;
                constituentClusters.push_back(*currentSeed);
                reco::CaloCluster_iterator currentCluster = currentSeed + 1;

                while (currentCluster != clusters_v.end())
                {

                        // if dynamic phi road is enabled then compute the phi road for a cluster
                        // of energy of existing clusters + the candidate cluster.
                        if (dynamicPhiRoad_)
                                phiRoad = phiRoadAlgo_->endcapPhiRoad(energy + (*currentCluster)->energy());

                        // does the cluster match the phi road for this candidate supercluster
                        if (match(*currentSeed, *currentCluster, etaRoad, phiRoad) &&
                                std::find(usedSeedDetIds.begin(), usedSeedDetIds.end(), (*currentCluster)->seed()) == usedSeedDetIds.end())
                        {

                                // add basic cluster to supercluster constituents
                                constituentClusters.push_back(*currentCluster);
                                energy   += (*currentCluster)->energy();
                                position_ += (*currentCluster)->energy() * math::XYZVector((*currentCluster)->position().X(),
                                                (*currentCluster)->position().Y(), 
                                                (*currentCluster)->position().Z());

                                // remove cluster from vector of available clusters
                                usedSeedDetIds.push_back((*currentCluster)->seed());
LogTrace("EcalClusters") << "Cluster R = " << (*currentCluster)->position().Rho();
                        }
                        ++currentCluster;

                }

                position_ /= energy;

LogTrace("EcalClusters") <<"Final SuperCluster R = " << position_.Rho();

                reco::SuperCluster newSuperCluster(energy, 
                                math::XYZPoint(position_.X(), position_.Y(), position_.Z()), 
                                (*currentSeed), 
                                constituentClusters);

                superclusters_v.push_back(newSuperCluster);
LogTrace("EcalClusters") << "created a new supercluster of: ";
LogTrace("EcalClusters") << "Energy = " << newSuperCluster.energy();
LogTrace("EcalClusters") << "Position in (R, phi, theta) = ("
                                << newSuperCluster.position().Rho() << ", "
                                << newSuperCluster.position().phi() << ", "
                                << newSuperCluster.position().theta() << ")" ;


        }

        clusters_v.clear();

}


bool Multi5x5BremRecoveryClusterAlgo::match(reco::CaloClusterPtr seed_p, 
                reco::CaloClusterPtr cluster_p,
                double dEtaMax, double dPhiMax)
{
        math::XYZPoint clusterPosition = cluster_p->position();
        math::XYZPoint seedPosition = seed_p->position();

        double dPhi = acos(cos(seedPosition.phi() - clusterPosition.phi()));

        double dEta = fabs(seedPosition.eta() - clusterPosition.eta());

        if (dEta > dEtaMax) return false;
        if (dPhi > dPhiMax) return false;

        return true;
}