| 1 |
// $Id: FillerSuperClusters.cc,v 1.14 2011/10/09 23:28:48 bendavid Exp $
|
| 2 |
|
| 3 |
#include "MitProd/TreeFiller/interface/FillerSuperClusters.h"
|
| 4 |
#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
|
| 5 |
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
|
| 6 |
#include "MitAna/DataTree/interface/BasicCluster.h"
|
| 7 |
#include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
|
| 8 |
#include "DataFormats/EgammaReco/interface/PreshowerClusterFwd.h"
|
| 9 |
#include "MitAna/DataTree/interface/SuperClusterCol.h"
|
| 10 |
#include "MitAna/DataTree/interface/Names.h"
|
| 11 |
#include "MitProd/ObjectService/interface/ObjectService.h"
|
| 12 |
#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
|
| 13 |
//#include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
|
| 14 |
#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
|
| 15 |
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
|
| 16 |
#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
|
| 17 |
#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
|
| 18 |
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
|
| 19 |
#include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
|
| 20 |
#include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
|
| 21 |
#include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
|
| 22 |
#include "Geometry/Records/interface/CaloGeometryRecord.h"
|
| 23 |
|
| 24 |
using namespace std;
|
| 25 |
using namespace edm;
|
| 26 |
using namespace mithep;
|
| 27 |
|
| 28 |
//--------------------------------------------------------------------------------------------------
|
| 29 |
FillerSuperClusters::FillerSuperClusters(const ParameterSet &cfg, const char *name, bool active) :
|
| 30 |
BaseFiller (cfg,name,active),
|
| 31 |
edmName_ (Conf().getUntrackedParameter<string>("edmName","hybridSuperClusters")),
|
| 32 |
mitName_ (Conf().getUntrackedParameter<string>("mitName","SuperClusters")),
|
| 33 |
basicClusterMapName_ (Conf().getUntrackedParameter<string>("basicClusterMapName",
|
| 34 |
"BasicClusterMap")),
|
| 35 |
psClusterMapName_ (Conf().getUntrackedParameter<string>("psClusterMapName",
|
| 36 |
"")),
|
| 37 |
caloTowerDetIdMapName_(Conf().getUntrackedParameter<string>("caloTowerDetIdMapName",
|
| 38 |
"CaloTowerDetIdMap")),
|
| 39 |
superClusterMapName_ (Conf().getUntrackedParameter<string>("superClusterMapName",
|
| 40 |
"SuperClusterMap")),
|
| 41 |
superClusterIdMapName_(Conf().getUntrackedParameter<string>("superClusterIdMapName",
|
| 42 |
"SuperClusterIdMap")),
|
| 43 |
caloTowerName_ (Conf().getUntrackedParameter<string>("caloTowerName","towerMaker")),
|
| 44 |
basicClusterMap_ (0),
|
| 45 |
psClusterMap_ (0),
|
| 46 |
caloTowerDetIdMap_ (0),
|
| 47 |
superClusters_ (new mithep::SuperClusterArr(25)),
|
| 48 |
superClusterMap_ (new mithep::SuperClusterMap),
|
| 49 |
superClusterIdMap_ (new mithep::SuperClusterIdMap)
|
| 50 |
{
|
| 51 |
// Constructor.
|
| 52 |
}
|
| 53 |
|
| 54 |
//--------------------------------------------------------------------------------------------------
|
| 55 |
FillerSuperClusters::~FillerSuperClusters()
|
| 56 |
{
|
| 57 |
// Destructor.
|
| 58 |
|
| 59 |
delete superClusters_;
|
| 60 |
delete superClusterMap_;
|
| 61 |
delete superClusterIdMap_;
|
| 62 |
}
|
| 63 |
|
| 64 |
//--------------------------------------------------------------------------------------------------
|
| 65 |
void FillerSuperClusters::BookDataBlock(TreeWriter &tws)
|
| 66 |
{
|
| 67 |
// Add super cluster branch to tree and get pointers to maps.
|
| 68 |
|
| 69 |
tws.AddBranch(mitName_,&superClusters_);
|
| 70 |
OS()->add<SuperClusterArr>(superClusters_,mitName_);
|
| 71 |
|
| 72 |
if (!basicClusterMapName_.empty()) {
|
| 73 |
basicClusterMap_ = OS()->get<BasicClusterMap>(basicClusterMapName_);
|
| 74 |
if (basicClusterMap_)
|
| 75 |
AddBranchDep(mitName_,basicClusterMap_->GetBrName());
|
| 76 |
}
|
| 77 |
|
| 78 |
if (!psClusterMapName_.empty()) {
|
| 79 |
psClusterMap_ = OS()->get<PsClusterMap>(psClusterMapName_);
|
| 80 |
if (psClusterMap_)
|
| 81 |
AddBranchDep(mitName_,psClusterMap_->GetBrName());
|
| 82 |
}
|
| 83 |
|
| 84 |
if (!caloTowerDetIdMapName_.empty()) {
|
| 85 |
caloTowerDetIdMap_ = OS()->get<CaloTowerDetIdMap>(caloTowerDetIdMapName_);
|
| 86 |
if (caloTowerDetIdMap_)
|
| 87 |
AddBranchDep(mitName_,caloTowerDetIdMap_->GetBrName());
|
| 88 |
}
|
| 89 |
|
| 90 |
if (!superClusterMapName_.empty()) {
|
| 91 |
superClusterMap_->SetBrName(mitName_);
|
| 92 |
OS()->add<SuperClusterMap>(superClusterMap_,superClusterMapName_);
|
| 93 |
}
|
| 94 |
if (!superClusterIdMapName_.empty()) {
|
| 95 |
superClusterIdMap_->SetBrName(mitName_);
|
| 96 |
OS()->add<SuperClusterIdMap>(superClusterIdMap_,superClusterIdMapName_);
|
| 97 |
}
|
| 98 |
}
|
| 99 |
|
| 100 |
//--------------------------------------------------------------------------------------------------
|
| 101 |
void FillerSuperClusters::FillDataBlock(const edm::Event &event,
|
| 102 |
const edm::EventSetup &setup)
|
| 103 |
{
|
| 104 |
// Fill the collection.
|
| 105 |
|
| 106 |
superClusters_->Delete();
|
| 107 |
superClusterMap_->Reset();
|
| 108 |
superClusterIdMap_->Reset();
|
| 109 |
|
| 110 |
Handle<reco::SuperClusterCollection> hSuperClusterProduct;
|
| 111 |
GetProduct(edmName_, hSuperClusterProduct, event);
|
| 112 |
|
| 113 |
Handle<CaloTowerCollection> hCaloTowerProduct;
|
| 114 |
GetProduct(caloTowerName_, hCaloTowerProduct, event);
|
| 115 |
|
| 116 |
superClusterMap_->SetEdmProductId(hSuperClusterProduct.id().id());
|
| 117 |
const reco::SuperClusterCollection inSuperClusters = *(hSuperClusterProduct.product());
|
| 118 |
|
| 119 |
//lazy tools for supercluster ecal timing
|
| 120 |
EcalClusterLazyTools lazyTools(event, setup,
|
| 121 |
edm::InputTag("reducedEcalRecHitsEB"),
|
| 122 |
edm::InputTag("reducedEcalRecHitsEE"));
|
| 123 |
|
| 124 |
////mustache id
|
| 125 |
//reco::Mustache mustache;
|
| 126 |
|
| 127 |
//es shape variables
|
| 128 |
edm::Handle<EcalRecHitCollection> hESRecHits;
|
| 129 |
event.getByLabel("reducedEcalRecHitsES" , hESRecHits);
|
| 130 |
|
| 131 |
edm::ESHandle<CaloGeometry> pGeometry;
|
| 132 |
setup.get<CaloGeometryRecord>().get(pGeometry);
|
| 133 |
|
| 134 |
edm::ESHandle<CaloTopology> pTopology;
|
| 135 |
setup.get<CaloTopologyRecord>().get(pTopology);
|
| 136 |
|
| 137 |
EcalPreshowerTopology topology_p(pGeometry);
|
| 138 |
|
| 139 |
//map of preshower rechits for shape calculations
|
| 140 |
std::map<DetId, EcalRecHit> esmap;
|
| 141 |
if (hESRecHits.isValid()) {
|
| 142 |
EcalRecHitCollection::const_iterator it;
|
| 143 |
for (it = hESRecHits->begin(); it != hESRecHits->end(); ++it) {
|
| 144 |
// remove bad ES rechits
|
| 145 |
if (it->recoFlag()==1 || it->recoFlag()==14 || (it->recoFlag()<=10 && it->recoFlag()>=5)) continue;
|
| 146 |
//Make the map of DetID, EcalRecHit pairs
|
| 147 |
esmap.insert(std::make_pair(it->id(), *it));
|
| 148 |
}
|
| 149 |
}
|
| 150 |
|
| 151 |
// loop through all super clusters
|
| 152 |
for (reco::SuperClusterCollection::const_iterator inSC = inSuperClusters.begin();
|
| 153 |
inSC != inSuperClusters.end(); ++inSC) {
|
| 154 |
|
| 155 |
mithep::SuperCluster *outSC = superClusters_->Allocate();
|
| 156 |
new (outSC) mithep::SuperCluster();
|
| 157 |
|
| 158 |
outSC->SetXYZ(inSC->x(),inSC->y(),inSC->z());
|
| 159 |
outSC->SetEnergy(inSC->energy());
|
| 160 |
outSC->SetRawEnergy(inSC->rawEnergy());
|
| 161 |
outSC->SetPreshowerEnergy(inSC->preshowerEnergy());
|
| 162 |
outSC->SetPhiWidth(inSC->phiWidth());
|
| 163 |
outSC->SetEtaWidth(inSC->etaWidth());
|
| 164 |
|
| 165 |
//Compute Hadronic Energy behind the supercluster (within DR < 0.15)
|
| 166 |
EgammaTowerIsolation towerIsoDepth1(0.15,0.,0.,1,hCaloTowerProduct.product()) ;
|
| 167 |
EgammaTowerIsolation towerIsoDepth2(0.15,0.,0.,2,hCaloTowerProduct.product()) ;
|
| 168 |
outSC->SetHcalDepth1Energy(towerIsoDepth1.getTowerESum(&(*inSC)));
|
| 169 |
outSC->SetHcalDepth2Energy(towerIsoDepth2.getTowerESum(&(*inSC)));
|
| 170 |
|
| 171 |
////ecal timing information
|
| 172 |
//outSC->SetTime (lazyTools.SuperClusterTime(*inSC, event));
|
| 173 |
//outSC->SetSeedTime(lazyTools.SuperClusterSeedTime(*inSC));
|
| 174 |
outSC->SetTime (-9999);
|
| 175 |
outSC->SetSeedTime(-9999);
|
| 176 |
|
| 177 |
//preshower shape variables
|
| 178 |
if (inSC->seed()->hitsAndFractions().at(0).first.subdetId()==EcalEndcap) {
|
| 179 |
std::vector<float> phoESShape = getESShape(getESHits(inSC->x(), inSC->y(), inSC->z(), esmap, *pGeometry.product(), &topology_p, 0));
|
| 180 |
|
| 181 |
if (phoESShape.size()>=2) {
|
| 182 |
outSC->SetPsEffWidthSigmaXX(phoESShape[0]);
|
| 183 |
outSC->SetPsEffWidthSigmaYY(phoESShape[1]);
|
| 184 |
}
|
| 185 |
}
|
| 186 |
|
| 187 |
// set the seed
|
| 188 |
if (basicClusterMap_ && inSC->seed().isNonnull())
|
| 189 |
outSC->SetSeed(basicClusterMap_->GetMit(inSC->seed()));
|
| 190 |
|
| 191 |
// add basic clusters that belong to this super cluster and tag with mustache id
|
| 192 |
if (basicClusterMap_) {
|
| 193 |
std::vector<reco::CaloCluster> bcs;
|
| 194 |
std::vector<unsigned int> insidebcs;
|
| 195 |
std::vector<unsigned int> outsidebcs;
|
| 196 |
for (reco::CaloCluster_iterator bc = inSC->clustersBegin(); bc != inSC->clustersEnd(); ++bc) {
|
| 197 |
if (bc->isNonnull()) {
|
| 198 |
outSC->AddCluster(basicClusterMap_->GetMit(*bc));
|
| 199 |
bcs.push_back(**bc);
|
| 200 |
}
|
| 201 |
}
|
| 202 |
|
| 203 |
//mustache.MustacheClust(bcs,insidebcs,outsidebcs);
|
| 204 |
//for (unsigned int iclus=0; iclus<insidebcs.size(); ++iclus)
|
| 205 |
// const_cast<mithep::BasicCluster*>(outSC->Cluster(insidebcs.at(iclus)))->SetInsideMustache(kTRUE);
|
| 206 |
//for (unsigned int iclus=0; iclus<outsidebcs.size(); ++iclus)
|
| 207 |
// const_cast<mithep::BasicCluster*>(outSC->Cluster(outsidebcs.at(iclus)))->SetInsideMustache(kFALSE);
|
| 208 |
}
|
| 209 |
|
| 210 |
// add preshower clusters that belong to this super cluster and tag with mustache id
|
| 211 |
if (psClusterMap_) {
|
| 212 |
std::vector<reco::CaloCluster> bcs;
|
| 213 |
std::vector<unsigned int> insidebcs;
|
| 214 |
std::vector<unsigned int> outsidebcs;
|
| 215 |
for(reco::CaloCluster_iterator bc = inSC->preshowerClustersBegin(); bc != inSC->preshowerClustersEnd(); ++bc) {
|
| 216 |
if (bc->isNonnull()) {
|
| 217 |
outSC->AddPsClust(psClusterMap_->GetMit(*bc));
|
| 218 |
bcs.push_back(**bc);
|
| 219 |
}
|
| 220 |
}
|
| 221 |
|
| 222 |
//mustache.MustacheClust(bcs,insidebcs,outsidebcs);
|
| 223 |
//for (unsigned int iclus=0; iclus<insidebcs.size(); ++iclus)
|
| 224 |
// const_cast<mithep::PsCluster*>(outSC->PsClust(insidebcs.at(iclus)))->SetInsideMustache(kTRUE);
|
| 225 |
//for (unsigned int iclus=0; iclus<outsidebcs.size(); ++iclus)
|
| 226 |
// const_cast<mithep::PsCluster*>(outSC->PsClust(outsidebcs.at(iclus)))->SetInsideMustache(kFALSE);
|
| 227 |
}
|
| 228 |
|
| 229 |
//compute preshower energy sums by plane
|
| 230 |
double psenergyplane1 = 0.;
|
| 231 |
double psenergyplane2 = 0.;
|
| 232 |
for(reco::CaloCluster_iterator bc = inSC->preshowerClustersBegin(); bc != inSC->preshowerClustersEnd(); ++bc) {
|
| 233 |
const reco::PreshowerCluster *pscluster = dynamic_cast<const reco::PreshowerCluster*>(bc->get());
|
| 234 |
if (!pscluster) continue;
|
| 235 |
|
| 236 |
if (pscluster->plane()==1) psenergyplane1 += pscluster->energy();
|
| 237 |
else if (pscluster->plane()==2) psenergyplane2 += pscluster->energy();
|
| 238 |
}
|
| 239 |
outSC->SetPreshowerEnergyPlane1(psenergyplane1);
|
| 240 |
outSC->SetPreshowerEnergyPlane2(psenergyplane2);
|
| 241 |
|
| 242 |
//add super cluster det ids to the id map and also fill supercluster-calotower associations
|
| 243 |
const std::vector< std::pair<DetId, float> > &pairs = inSC->hitsAndFractions();
|
| 244 |
for (std::vector< std::pair<DetId, float> >::const_iterator ipair = pairs.begin();
|
| 245 |
ipair < pairs.end(); ++ipair) {
|
| 246 |
|
| 247 |
const DetId &ihit = ipair->first;
|
| 248 |
|
| 249 |
if (caloTowerDetIdMap_) {
|
| 250 |
if (caloTowerDetIdMap_->HasMit(ihit)) {
|
| 251 |
const mithep::CaloTower *matchedTower = caloTowerDetIdMap_->GetMit(ihit);
|
| 252 |
if (!outSC->HasTower(matchedTower)) {
|
| 253 |
outSC->AddTower(matchedTower);
|
| 254 |
}
|
| 255 |
}
|
| 256 |
}
|
| 257 |
superClusterIdMap_->Add(ihit,outSC);
|
| 258 |
}
|
| 259 |
|
| 260 |
// add super cluster to the map
|
| 261 |
reco::SuperClusterRef theRef(hSuperClusterProduct, inSC-inSuperClusters.begin());
|
| 262 |
superClusterMap_->Add(theRef, outSC);
|
| 263 |
|
| 264 |
}
|
| 265 |
superClusters_->Trim();
|
| 266 |
}
|
| 267 |
|
| 268 |
//horrible code below copied from globe for preshower cluster shape calculations
|
| 269 |
std::vector<float> FillerSuperClusters::getESHits(double X, double Y, double Z, std::map<DetId, EcalRecHit> rechits_map, const CaloGeometry& geometry, CaloSubdetectorTopology *topology_p, int row) {
|
| 270 |
std::vector<float> esHits;
|
| 271 |
|
| 272 |
const GlobalPoint point(X,Y,Z);
|
| 273 |
|
| 274 |
const CaloSubdetectorGeometry *geometry_p ;
|
| 275 |
geometry_p = geometry.getSubdetectorGeometry (DetId::Ecal,EcalPreshower) ;
|
| 276 |
|
| 277 |
DetId esId1 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 1);
|
| 278 |
DetId esId2 = (dynamic_cast<const EcalPreshowerGeometry*>(geometry_p))->getClosestCellInPlane(point, 2);
|
| 279 |
ESDetId esDetId1 = (esId1 == DetId(0)) ? ESDetId(0) : ESDetId(esId1);
|
| 280 |
ESDetId esDetId2 = (esId2 == DetId(0)) ? ESDetId(0) : ESDetId(esId2);
|
| 281 |
|
| 282 |
std::map<DetId, EcalRecHit>::iterator it;
|
| 283 |
ESDetId next;
|
| 284 |
ESDetId strip1;
|
| 285 |
ESDetId strip2;
|
| 286 |
|
| 287 |
strip1 = esDetId1;
|
| 288 |
strip2 = esDetId2;
|
| 289 |
|
| 290 |
EcalPreshowerNavigator theESNav1(strip1, topology_p);
|
| 291 |
theESNav1.setHome(strip1);
|
| 292 |
|
| 293 |
EcalPreshowerNavigator theESNav2(strip2, topology_p);
|
| 294 |
theESNav2.setHome(strip2);
|
| 295 |
|
| 296 |
if (row == 1) {
|
| 297 |
if (strip1 != ESDetId(0)) strip1 = theESNav1.north();
|
| 298 |
if (strip2 != ESDetId(0)) strip2 = theESNav2.east();
|
| 299 |
} else if (row == -1) {
|
| 300 |
if (strip1 != ESDetId(0)) strip1 = theESNav1.south();
|
| 301 |
if (strip2 != ESDetId(0)) strip2 = theESNav2.west();
|
| 302 |
}
|
| 303 |
|
| 304 |
// Plane 1
|
| 305 |
if (strip1 == ESDetId(0)) {
|
| 306 |
for (int i=0; i<31; ++i) esHits.push_back(0);
|
| 307 |
} else {
|
| 308 |
|
| 309 |
it = rechits_map.find(strip1);
|
| 310 |
if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
|
| 311 |
else esHits.push_back(0);
|
| 312 |
//cout<<"center : "<<strip1<<" "<<it->second.energy()<<endl;
|
| 313 |
|
| 314 |
// east road
|
| 315 |
for (int i=0; i<15; ++i) {
|
| 316 |
next = theESNav1.east();
|
| 317 |
if (next != ESDetId(0)) {
|
| 318 |
it = rechits_map.find(next);
|
| 319 |
if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
|
| 320 |
else esHits.push_back(0);
|
| 321 |
//cout<<"east "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
|
| 322 |
} else {
|
| 323 |
for (int j=i; j<15; j++) esHits.push_back(0);
|
| 324 |
break;
|
| 325 |
//cout<<"east "<<i<<" : "<<next<<" "<<0<<endl;
|
| 326 |
}
|
| 327 |
}
|
| 328 |
|
| 329 |
// west road
|
| 330 |
theESNav1.setHome(strip1);
|
| 331 |
theESNav1.home();
|
| 332 |
for (int i=0; i<15; ++i) {
|
| 333 |
next = theESNav1.west();
|
| 334 |
if (next != ESDetId(0)) {
|
| 335 |
it = rechits_map.find(next);
|
| 336 |
if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
|
| 337 |
else esHits.push_back(0);
|
| 338 |
//cout<<"west "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
|
| 339 |
} else {
|
| 340 |
for (int j=i; j<15; j++) esHits.push_back(0);
|
| 341 |
break;
|
| 342 |
//cout<<"west "<<i<<" : "<<next<<" "<<0<<endl;
|
| 343 |
}
|
| 344 |
}
|
| 345 |
}
|
| 346 |
|
| 347 |
if (strip2 == ESDetId(0)) {
|
| 348 |
for (int i=0; i<31; ++i) esHits.push_back(0);
|
| 349 |
} else {
|
| 350 |
|
| 351 |
it = rechits_map.find(strip2);
|
| 352 |
if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
|
| 353 |
else esHits.push_back(0);
|
| 354 |
//cout<<"center : "<<strip2<<" "<<it->second.energy()<<endl;
|
| 355 |
|
| 356 |
// north road
|
| 357 |
for (int i=0; i<15; ++i) {
|
| 358 |
next = theESNav2.north();
|
| 359 |
if (next != ESDetId(0)) {
|
| 360 |
it = rechits_map.find(next);
|
| 361 |
if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
|
| 362 |
else esHits.push_back(0);
|
| 363 |
//cout<<"north "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
|
| 364 |
} else {
|
| 365 |
for (int j=i; j<15; j++) esHits.push_back(0);
|
| 366 |
break;
|
| 367 |
//cout<<"north "<<i<<" : "<<next<<" "<<0<<endl;
|
| 368 |
}
|
| 369 |
}
|
| 370 |
|
| 371 |
// south road
|
| 372 |
theESNav2.setHome(strip2);
|
| 373 |
theESNav2.home();
|
| 374 |
for (int i=0; i<15; ++i) {
|
| 375 |
next = theESNav2.south();
|
| 376 |
if (next != ESDetId(0)) {
|
| 377 |
it = rechits_map.find(next);
|
| 378 |
if (it->second.energy() > 1.0e-10 && it != rechits_map.end()) esHits.push_back(it->second.energy());
|
| 379 |
else esHits.push_back(0);
|
| 380 |
//cout<<"south "<<i<<" : "<<next<<" "<<it->second.energy()<<endl;
|
| 381 |
} else {
|
| 382 |
for (int j=i; j<15; j++) esHits.push_back(0);
|
| 383 |
break;
|
| 384 |
//cout<<"south "<<i<<" : "<<next<<" "<<0<<endl;
|
| 385 |
}
|
| 386 |
}
|
| 387 |
}
|
| 388 |
|
| 389 |
return esHits;
|
| 390 |
}
|
| 391 |
|
| 392 |
std::vector<float> FillerSuperClusters::getESShape(std::vector<float> ESHits0)
|
| 393 |
{
|
| 394 |
std::vector<float> esShape;
|
| 395 |
|
| 396 |
const int nBIN = 21;
|
| 397 |
float esRH_F[nBIN];
|
| 398 |
float esRH_R[nBIN];
|
| 399 |
for (int idx=0; idx<nBIN; idx++) {
|
| 400 |
esRH_F[idx] = 0.;
|
| 401 |
esRH_R[idx] = 0.;
|
| 402 |
}
|
| 403 |
|
| 404 |
for(int ibin=0; ibin<((nBIN+1)/2); ibin++) {
|
| 405 |
if (ibin==0) {
|
| 406 |
esRH_F[(nBIN-1)/2] = ESHits0[ibin];
|
| 407 |
esRH_R[(nBIN-1)/2] = ESHits0[ibin+31];
|
| 408 |
} else {
|
| 409 |
esRH_F[(nBIN-1)/2+ibin] = ESHits0[ibin];
|
| 410 |
esRH_F[(nBIN-1)/2-ibin] = ESHits0[ibin+15];
|
| 411 |
esRH_R[(nBIN-1)/2+ibin] = ESHits0[ibin+31];
|
| 412 |
esRH_R[(nBIN-1)/2-ibin] = ESHits0[ibin+31+15];
|
| 413 |
}
|
| 414 |
}
|
| 415 |
|
| 416 |
// ---- Effective Energy Deposit Width ---- //
|
| 417 |
double EffWidthSigmaXX = 0.;
|
| 418 |
double EffWidthSigmaYY = 0.;
|
| 419 |
double totalEnergyXX = 0.;
|
| 420 |
double totalEnergyYY = 0.;
|
| 421 |
double EffStatsXX = 0.;
|
| 422 |
double EffStatsYY = 0.;
|
| 423 |
for (int id_X=0; id_X<21; id_X++) {
|
| 424 |
totalEnergyXX += esRH_F[id_X];
|
| 425 |
EffStatsXX += esRH_F[id_X]*(id_X-10)*(id_X-10);
|
| 426 |
totalEnergyYY += esRH_R[id_X];
|
| 427 |
EffStatsYY += esRH_R[id_X]*(id_X-10)*(id_X-10);
|
| 428 |
}
|
| 429 |
EffWidthSigmaXX = (totalEnergyXX>0.) ? sqrt(fabs(EffStatsXX / totalEnergyXX)) : 0.;
|
| 430 |
EffWidthSigmaYY = (totalEnergyYY>0.) ? sqrt(fabs(EffStatsYY / totalEnergyYY)) : 0.;
|
| 431 |
|
| 432 |
esShape.push_back(EffWidthSigmaXX);
|
| 433 |
esShape.push_back(EffWidthSigmaYY);
|
| 434 |
|
| 435 |
return esShape;
|
| 436 |
}
|
| 437 |
|
