File indexing completed on 2021-09-03 22:28:00
0001 #include "Calibration/EcalCalibAlgos/interface/Pi0FixedMassWindowCalibration.h"
0002
0003
0004
0005
0006
0007
0008
0009 #include "Calibration/Tools/interface/Pi0CalibXMLwriter.h"
0010
0011
0012 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0013 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0014
0015 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0016 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0017 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0018 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0019
0020
0021 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0022 #include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"
0023
0024 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0025
0026
0027
0028 using namespace std;
0029
0030
0031
0032 Pi0FixedMassWindowCalibration::Pi0FixedMassWindowCalibration(const edm::ParameterSet& iConfig)
0033 : theMaxLoops(iConfig.getUntrackedParameter<unsigned int>("maxLoops", 0)),
0034 ecalHitsProducer_(iConfig.getParameter<std::string>("ecalRecHitsProducer")),
0035 barrelHits_(iConfig.getParameter<std::string>("barrelHitCollection")),
0036 recHitToken_(consumes<EcalRecHitCollection>(edm::InputTag(ecalHitsProducer_, barrelHits_))),
0037 intercalibConstantsToken_(esConsumes()),
0038 geometryToken_(esConsumes()) {
0039 std::cout << "[Pi0FixedMassWindowCalibration] Constructor " << std::endl;
0040
0041 std::string verbosityString = iConfig.getParameter<std::string>("VerbosityLevel");
0042 if (verbosityString == "DEBUG")
0043 verbosity = IslandClusterAlgo::pDEBUG;
0044 else if (verbosityString == "WARNING")
0045 verbosity = IslandClusterAlgo::pWARNING;
0046 else if (verbosityString == "INFO")
0047 verbosity = IslandClusterAlgo::pINFO;
0048 else
0049 verbosity = IslandClusterAlgo::pERROR;
0050
0051
0052 barrelClusterCollection_ = iConfig.getParameter<std::string>("barrelClusterCollection");
0053
0054
0055 double barrelSeedThreshold = iConfig.getParameter<double>("IslandBarrelSeedThr");
0056 double endcapSeedThreshold = iConfig.getParameter<double>("IslandEndcapSeedThr");
0057
0058
0059 selePi0PtGammaOneMin_ = iConfig.getParameter<double>("selePi0PtGammaOneMin");
0060 selePi0PtGammaTwoMin_ = iConfig.getParameter<double>("selePi0PtGammaTwoMin");
0061
0062 selePi0DRBelt_ = iConfig.getParameter<double>("selePi0DRBelt");
0063 selePi0DetaBelt_ = iConfig.getParameter<double>("selePi0DetaBelt");
0064
0065 selePi0PtPi0Min_ = iConfig.getParameter<double>("selePi0PtPi0Min");
0066
0067 selePi0S4S9GammaOneMin_ = iConfig.getParameter<double>("selePi0S4S9GammaOneMin");
0068 selePi0S4S9GammaTwoMin_ = iConfig.getParameter<double>("selePi0S4S9GammaTwoMin");
0069 selePi0S9S25GammaOneMin_ = iConfig.getParameter<double>("selePi0S9S25GammaOneMin");
0070 selePi0S9S25GammaTwoMin_ = iConfig.getParameter<double>("selePi0S9S25GammaTwoMin");
0071
0072 selePi0EtBeltIsoRatioMax_ = iConfig.getParameter<double>("selePi0EtBeltIsoRatioMax");
0073
0074 selePi0MinvMeanFixed_ = iConfig.getParameter<double>("selePi0MinvMeanFixed");
0075 selePi0MinvSigmaFixed_ = iConfig.getParameter<double>("selePi0MinvSigmaFixed");
0076
0077
0078 edm::ParameterSet posCalcParameters = iConfig.getParameter<edm::ParameterSet>("posCalcParameters");
0079 posCalculator_ = PositionCalc(posCalcParameters);
0080 shapeAlgo_ = ClusterShapeAlgo(posCalcParameters);
0081 clustershapecollectionEB_ = iConfig.getParameter<std::string>("clustershapecollectionEB");
0082
0083
0084 barrelClusterShapeAssociation_ = iConfig.getParameter<std::string>("barrelShapeAssociation");
0085
0086 const std::vector<std::string> seedflagnamesEB =
0087 iConfig.getParameter<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEB");
0088 const std::vector<int> seedflagsexclEB = StringToEnumValue<EcalRecHit::Flags>(seedflagnamesEB);
0089
0090 const std::vector<std::string> seedflagnamesEE =
0091 iConfig.getParameter<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEE");
0092 const std::vector<int> seedflagsexclEE = StringToEnumValue<EcalRecHit::Flags>(seedflagnamesEE);
0093
0094 const std::vector<std::string> flagnamesEB =
0095 iConfig.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEB");
0096 const std::vector<int> flagsexclEB = StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
0097
0098 const std::vector<std::string> flagnamesEE =
0099 iConfig.getParameter<std::vector<std::string>>("RecHitFlagToBeExcludedEE");
0100 const std::vector<int> flagsexclEE = StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
0101
0102 island_p = new IslandClusterAlgo(barrelSeedThreshold,
0103 endcapSeedThreshold,
0104 posCalculator_,
0105 seedflagsexclEB,
0106 seedflagsexclEE,
0107 flagsexclEB,
0108 flagsexclEE,
0109 verbosity);
0110
0111 theParameterSet = iConfig;
0112 }
0113
0114
0115
0116
0117 Pi0FixedMassWindowCalibration::~Pi0FixedMassWindowCalibration() { delete island_p; }
0118
0119 void Pi0FixedMassWindowCalibration::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0120 edm::ParameterSetDescription desc;
0121
0122 desc.add<unsigned int>("maxLoops", 0);
0123 desc.add<std::string>("ecalRecHitsProducer", "");
0124 desc.add<std::string>("barrelHitCollection", "");
0125
0126 desc.add<std::string>("VerbosityLevel", "");
0127 desc.add<std::string>("barrelClusterCollection", "");
0128
0129 desc.add<double>("IslandBarrelSeedThr", 0);
0130 desc.add<double>("IslandEndcapSeedThr", 0);
0131
0132 desc.add<double>("selePi0PtGammaOneMin", 0);
0133 desc.add<double>("selePi0PtGammaTwoMin", 0);
0134
0135 desc.add<double>("selePi0DRBelt", 0);
0136 desc.add<double>("selePi0DetaBelt", 0);
0137
0138 desc.add<double>("selePi0PtPi0Min", 0);
0139
0140 desc.add<double>("selePi0S4S9GammaOneMin", 0);
0141 desc.add<double>("selePi0S4S9GammaTwoMin", 0);
0142 desc.add<double>("selePi0S9S25GammaOneMin", 0);
0143 desc.add<double>("selePi0S9S25GammaTwoMin", 0);
0144
0145 desc.add<double>("selePi0EtBeltIsoRatioMax", 0);
0146
0147 desc.add<double>("selePi0MinvMeanFixed", 0);
0148 desc.add<double>("selePi0MinvSigmaFixed", 0);
0149
0150 edm::ParameterSetDescription posCalcParameters;
0151 posCalcParameters.add<bool>("LogWeighted", true);
0152 posCalcParameters.add<double>("T0_barl", 7.4);
0153 posCalcParameters.add<double>("T0_endc", 3.1);
0154 posCalcParameters.add<double>("T0_endcPresh", 1.2);
0155 posCalcParameters.add<double>("W0", 4.2);
0156 posCalcParameters.add<double>("X0", 0.89);
0157 desc.add<edm::ParameterSetDescription>("posCalcParameters", posCalcParameters);
0158
0159 desc.add<std::string>("clustershapecollectionEB", "islandBarrelShape");
0160 desc.add<std::string>("barrelShapeAssociation", "islandBarrelShapeAssoc");
0161
0162 desc.add<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEB", {});
0163 desc.add<std::vector<std::string>>("SeedRecHitFlagToBeExcludedEE", {});
0164 desc.add<std::vector<std::string>>("RecHitFlagToBeExcludedEB", {});
0165 desc.add<std::vector<std::string>>("RecHitFlagToBeExcludedEE", {});
0166
0167 descriptions.add("Pi0FixedMassWindowCalibration", desc);
0168 }
0169
0170
0171
0172
0173 void Pi0FixedMassWindowCalibration::beginOfJob() {
0174
0175
0176 isfirstcall_ = true;
0177 }
0178
0179 void Pi0FixedMassWindowCalibration::endOfJob() {
0180 std::cout << "[Pi0FixedMassWindowCalibration] endOfJob" << endl;
0181
0182
0183
0184 Pi0CalibXMLwriter barrelWriter(EcalBarrel, 99);
0185
0186 std::vector<DetId>::const_iterator barrelIt = barrelCells.begin();
0187 for (; barrelIt != barrelCells.end(); barrelIt++) {
0188 EBDetId eb(*barrelIt);
0189 int ieta = eb.ieta();
0190 int iphi = eb.iphi();
0191 int sign = eb.zside() > 0 ? 1 : 0;
0192 barrelWriter.writeLine(eb, newCalibs_barl[abs(ieta) - 1][iphi - 1][sign]);
0193 }
0194 }
0195
0196
0197
0198 void Pi0FixedMassWindowCalibration::startingNewLoop(unsigned int iLoop) {
0199 for (int sign = 0; sign < 2; sign++) {
0200 for (int ieta = 0; ieta < 85; ieta++) {
0201 for (int iphi = 0; iphi < 360; iphi++) {
0202 wxtals[ieta][iphi][sign] = 0.;
0203 mwxtals[ieta][iphi][sign] = 0.;
0204 }
0205 }
0206 }
0207 std::cout << "[Pi0FixedMassWindowCalibration] Starting loop number " << iLoop << std::endl;
0208 }
0209
0210
0211
0212
0213 edm::EDLooper::Status Pi0FixedMassWindowCalibration::endOfLoop(const edm::EventSetup& iSetup, unsigned int iLoop) {
0214 std::cout << "[Pi0FixedMassWindowCalibration] Ending loop " << iLoop << std::endl;
0215
0216 for (int sign = 0; sign < 2; sign++) {
0217 for (int ieta = 0; ieta < 85; ieta++) {
0218 for (int iphi = 0; iphi < 360; iphi++) {
0219 if (wxtals[ieta][iphi][sign] == 0) {
0220 newCalibs_barl[ieta][iphi][sign] = oldCalibs_barl[ieta][iphi][sign];
0221 } else {
0222 newCalibs_barl[ieta][iphi][sign] =
0223 oldCalibs_barl[ieta][iphi][sign] * (mwxtals[ieta][iphi][sign] / wxtals[ieta][iphi][sign]);
0224 }
0225 cout << " New calibration constant: ieta iphi sign - old,mwxtals,wxtals,new: " << ieta << " " << iphi << " "
0226 << sign << " - " << oldCalibs_barl[ieta][iphi][sign] << " " << mwxtals[ieta][iphi][sign] << " "
0227 << wxtals[ieta][iphi][sign] << " " << newCalibs_barl[ieta][iphi][sign] << endl;
0228 }
0229 }
0230 }
0231
0232 Pi0CalibXMLwriter barrelWriter(EcalBarrel, iLoop + 1);
0233
0234 std::vector<DetId>::const_iterator barrelIt = barrelCells.begin();
0235 for (; barrelIt != barrelCells.end(); barrelIt++) {
0236 EBDetId eb(*barrelIt);
0237 int ieta = eb.ieta();
0238 int iphi = eb.iphi();
0239 int sign = eb.zside() > 0 ? 1 : 0;
0240 barrelWriter.writeLine(eb, newCalibs_barl[abs(ieta) - 1][iphi - 1][sign]);
0241 if (iphi == 1) {
0242 std::cout << "Calib constant for barrel crystal "
0243 << " (" << ieta << "," << iphi << ") changed from " << oldCalibs_barl[abs(ieta) - 1][iphi - 1][sign]
0244 << " to " << newCalibs_barl[abs(ieta) - 1][iphi - 1][sign] << std::endl;
0245 }
0246 }
0247
0248
0249
0250 for (int sign = 0; sign < 2; sign++) {
0251 for (int ieta = 0; ieta < 85; ieta++) {
0252 for (int iphi = 0; iphi < 360; iphi++) {
0253 oldCalibs_barl[ieta][iphi][sign] = newCalibs_barl[ieta][iphi][sign];
0254 newCalibs_barl[ieta][iphi][sign] = 0;
0255 }
0256 }
0257 }
0258
0259 if (iLoop == theMaxLoops - 1 || iLoop >= theMaxLoops)
0260 return kStop;
0261 else
0262 return kContinue;
0263 }
0264
0265
0266
0267
0268 edm::EDLooper::Status Pi0FixedMassWindowCalibration::duringLoop(const edm::Event& event, const edm::EventSetup& setup) {
0269 using namespace edm;
0270 using namespace std;
0271
0272
0273
0274
0275
0276 const auto& geometry = setup.getData(geometryToken_);
0277
0278 if (isfirstcall_) {
0279
0280
0281 for (int sign = 0; sign < 2; sign++) {
0282 for (int ieta = 0; ieta < 85; ieta++) {
0283 for (int iphi = 0; iphi < 360; iphi++) {
0284 oldCalibs_barl[ieta][iphi][sign] = 0.;
0285 newCalibs_barl[ieta][iphi][sign] = 0.;
0286 wxtals[ieta][iphi][sign] = 0.;
0287 mwxtals[ieta][iphi][sign] = 0.;
0288 }
0289 }
0290 }
0291
0292
0293
0294 const auto& pIcal = setup.getData(intercalibConstantsToken_);
0295 const auto& imap = pIcal.getMap();
0296 std::cout << "imap.size() = " << imap.size() << std::endl;
0297
0298
0299 barrelCells = geometry.getValidDetIds(DetId::Ecal, EcalBarrel);
0300 std::vector<DetId>::const_iterator barrelIt;
0301 for (barrelIt = barrelCells.begin(); barrelIt != barrelCells.end(); barrelIt++) {
0302 EBDetId eb(*barrelIt);
0303
0304
0305 EcalIntercalibConstantMap::const_iterator itcalib = imap.find(eb.rawId());
0306 if (itcalib == imap.end()) {
0307
0308 }
0309 EcalIntercalibConstant calib = (*itcalib);
0310 int sign = eb.zside() > 0 ? 1 : 0;
0311 oldCalibs_barl[abs(eb.ieta()) - 1][eb.iphi() - 1][sign] = calib;
0312 if (eb.iphi() == 1)
0313 std::cout << "Read old constant for crystal "
0314 << " (" << eb.ieta() << "," << eb.iphi() << ") : " << calib << std::endl;
0315 }
0316 isfirstcall_ = false;
0317 }
0318
0319 nevent++;
0320
0321 if ((nevent < 100 && nevent % 10 == 0) || (nevent < 1000 && nevent % 100 == 0) ||
0322 (nevent < 10000 && nevent % 100 == 0) || (nevent < 100000 && nevent % 1000 == 0) ||
0323 (nevent < 10000000 && nevent % 1000 == 0))
0324 std::cout << "[Pi0FixedMassWindowCalibration] Events processed: " << nevent << std::endl;
0325
0326 recHitsEB_map = new std::map<DetId, EcalRecHit>();
0327
0328 EcalRecHitCollection* recalibEcalRecHitCollection(new EcalRecHitCollection);
0329
0330 int nRecHitsEB = 0;
0331 Handle<EcalRecHitCollection> pEcalRecHitBarrelCollection;
0332 event.getByToken(recHitToken_, pEcalRecHitBarrelCollection);
0333 const EcalRecHitCollection* ecalRecHitBarrelCollection = pEcalRecHitBarrelCollection.product();
0334 cout << " ECAL Barrel RecHits # " << ecalRecHitBarrelCollection->size() << endl;
0335 for (EcalRecHitCollection::const_iterator aRecHitEB = ecalRecHitBarrelCollection->begin();
0336 aRecHitEB != ecalRecHitBarrelCollection->end();
0337 aRecHitEB++) {
0338
0339
0340 EBDetId ebrhdetid = aRecHitEB->detid();
0341
0342
0343
0344
0345 int sign = ebrhdetid.zside() > 0 ? 1 : 0;
0346 EcalRecHit aHit(aRecHitEB->id(),
0347 aRecHitEB->energy() * oldCalibs_barl[abs(ebrhdetid.ieta()) - 1][ebrhdetid.iphi() - 1][sign],
0348 aRecHitEB->time());
0349 recalibEcalRecHitCollection->push_back(aHit);
0350
0351 nRecHitsEB++;
0352 }
0353
0354
0355 int irecalib = 0;
0356 for (EcalRecHitCollection::const_iterator aRecHitEB = recalibEcalRecHitCollection->begin();
0357 aRecHitEB != recalibEcalRecHitCollection->end();
0358 aRecHitEB++) {
0359
0360
0361
0362
0363
0364
0365
0366 std::pair<DetId, EcalRecHit> map_entry(aRecHitEB->id(), *aRecHitEB);
0367 recHitsEB_map->insert(map_entry);
0368
0369 irecalib++;
0370 }
0371
0372 const CaloSubdetectorGeometry* geometry_p;
0373
0374 std::string clustershapetag;
0375 geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0376 EcalBarrelTopology topology{geometry};
0377
0378 const CaloSubdetectorGeometry* geometryES_p;
0379 geometryES_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407 reco::BasicClusterCollection clusters_recalib;
0408 clusters_recalib = island_p->makeClusters(
0409 recalibEcalRecHitCollection, geometry_p, &topology, geometryES_p, IslandClusterAlgo::barrel);
0410
0411
0412 std::vector<reco::ClusterShape> ClusVec_recalib;
0413 for (int erg = 0; erg < int(clusters_recalib.size()); ++erg) {
0414 reco::ClusterShape TestShape_recalib =
0415 shapeAlgo_.Calculate(clusters_recalib[erg], recalibEcalRecHitCollection, geometry_p, &topology);
0416 ClusVec_recalib.push_back(TestShape_recalib);
0417 }
0418
0419
0420 std::unique_ptr<reco::ClusterShapeCollection> clustersshapes_p_recalib(new reco::ClusterShapeCollection);
0421 clustersshapes_p_recalib->assign(ClusVec_recalib.begin(), ClusVec_recalib.end());
0422
0423 cout << "[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster collection size: " << clusters_recalib.size()
0424 << endl;
0425 cout << "[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster Shape Collection size: "
0426 << clustersshapes_p_recalib->size() << endl;
0427
0428
0429
0430
0431
0432 static const int MAXBCEB = 200;
0433 static const int MAXBCEBRH = 200;
0434 int nIslandBCEB;
0435 float eIslandBCEB[MAXBCEB];
0436 float etIslandBCEB[MAXBCEB];
0437 float etaIslandBCEB[MAXBCEB];
0438 float phiIslandBCEB[MAXBCEB];
0439 float e2x2IslandBCEB[MAXBCEB];
0440 float e3x3IslandBCEB[MAXBCEB];
0441 float e5x5IslandBCEB[MAXBCEB];
0442
0443
0444 int nIslandBCEBRecHits[MAXBCEB];
0445
0446 int ietaIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0447 int iphiIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0448 int zsideIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0449 float eIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0450
0451 nIslandBCEB = 0;
0452 for (int i = 0; i < MAXBCEB; i++) {
0453 eIslandBCEB[i] = 0;
0454 etIslandBCEB[i] = 0;
0455 etaIslandBCEB[i] = 0;
0456 phiIslandBCEB[i] = 0;
0457 e2x2IslandBCEB[i] = 0;
0458 e3x3IslandBCEB[i] = 0;
0459 e5x5IslandBCEB[i] = 0;
0460 nIslandBCEBRecHits[i] = 0;
0461 for (int j = 0; j < MAXBCEBRH; j++) {
0462
0463 ietaIslandBCEBRecHits[i][j] = 0;
0464 iphiIslandBCEBRecHits[i][j] = 0;
0465 zsideIslandBCEBRecHits[i][j] = 0;
0466 eIslandBCEBRecHits[i][j] = 0;
0467 }
0468 }
0469
0470 int iClus_recalib = 0;
0471 for (reco::BasicClusterCollection::const_iterator aClus = clusters_recalib.begin(); aClus != clusters_recalib.end();
0472 aClus++) {
0473 cout << " CLUSTER [recalibration] : #,NHits,e,et,eta,phi,e2x2,e3x3,e5x5: " << iClus_recalib << " " << aClus->size()
0474 << " " << aClus->energy() << " " << aClus->energy() * sin(aClus->position().theta()) << " "
0475 << aClus->position().eta() << " " << aClus->position().phi() << " "
0476 << (*clustersshapes_p_recalib)[iClus_recalib].e2x2() << " "
0477 << (*clustersshapes_p_recalib)[iClus_recalib].e3x3() << " "
0478 << (*clustersshapes_p_recalib)[iClus_recalib].e5x5() << endl;
0479
0480 eIslandBCEB[nIslandBCEB] = aClus->energy();
0481 etIslandBCEB[nIslandBCEB] = aClus->energy() * sin(aClus->position().theta());
0482 etaIslandBCEB[nIslandBCEB] = aClus->position().eta();
0483 phiIslandBCEB[nIslandBCEB] = aClus->position().phi();
0484
0485 e2x2IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e2x2();
0486 e3x3IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e3x3();
0487 e5x5IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e5x5();
0488
0489 nIslandBCEBRecHits[nIslandBCEB] = aClus->size();
0490
0491 std::vector<std::pair<DetId, float>> hits = aClus->hitsAndFractions();
0492 std::vector<std::pair<DetId, float>>::iterator hit;
0493 std::map<DetId, EcalRecHit>::iterator aHit;
0494
0495 int irhcount = 0;
0496 for (hit = hits.begin(); hit != hits.end(); hit++) {
0497
0498 aHit = recHitsEB_map->find((*hit).first);
0499
0500
0501 EBDetId sel_rh = aHit->second.detid();
0502
0503
0504
0505
0506 ietaIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.ieta();
0507 iphiIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.iphi();
0508 zsideIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.zside();
0509 eIslandBCEBRecHits[nIslandBCEB][irhcount] = aHit->second.energy();
0510
0511 irhcount++;
0512 }
0513 nIslandBCEB++;
0514 iClus_recalib++;
0515 }
0516
0517
0518
0519 if (nIslandBCEB > 1) {
0520 for (int i = 0; i < nIslandBCEB; i++) {
0521 for (int j = i + 1; j < nIslandBCEB; j++) {
0522 if (etIslandBCEB[i] > selePi0PtGammaOneMin_ && etIslandBCEB[j] > selePi0PtGammaOneMin_) {
0523 float theta_0 = 2. * atan(exp(-etaIslandBCEB[i]));
0524 float theta_1 = 2. * atan(exp(-etaIslandBCEB[j]));
0525
0526 float p0x = eIslandBCEB[i] * sin(theta_0) * cos(phiIslandBCEB[i]);
0527 float p1x = eIslandBCEB[j] * sin(theta_1) * cos(phiIslandBCEB[j]);
0528
0529 float p0y = eIslandBCEB[i] * sin(theta_0) * sin(phiIslandBCEB[i]);
0530 float p1y = eIslandBCEB[j] * sin(theta_1) * sin(phiIslandBCEB[j]);
0531
0532 float p0z = eIslandBCEB[i] * cos(theta_0);
0533 float p1z = eIslandBCEB[j] * cos(theta_1);
0534
0535 float pi0_px = p0x + p1x;
0536 float pi0_py = p0y + p1y;
0537 float pi0_pz = p0z + p1z;
0538
0539 float pi0_ptot = sqrt(pi0_px * pi0_px + pi0_py * pi0_py + pi0_pz * pi0_pz);
0540
0541 float pi0_theta = acos(pi0_pz / pi0_ptot);
0542 float pi0_eta = -log(tan(pi0_theta / 2));
0543 float pi0_phi = atan(pi0_py / pi0_px);
0544
0545
0546
0547
0548 float et_belt = 0;
0549 for (Int_t k = 0; k < nIslandBCEB; k++) {
0550 if ((k != i) && (k != j)) {
0551 float dr_pi0_k = sqrt((etaIslandBCEB[k] - pi0_eta) * (etaIslandBCEB[k] - pi0_eta) +
0552 (phiIslandBCEB[k] - pi0_phi) * (phiIslandBCEB[k] - pi0_phi));
0553 float deta_pi0_k = fabs(etaIslandBCEB[k] - pi0_eta);
0554 if ((dr_pi0_k < selePi0DRBelt_) && (deta_pi0_k < selePi0DetaBelt_))
0555 et_belt = et_belt + etIslandBCEB[k];
0556 }
0557 }
0558
0559 float pt_pi0 = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
0560
0561
0562
0563 if (pt_pi0 > selePi0PtPi0Min_) {
0564 float m_inv = sqrt((eIslandBCEB[i] + eIslandBCEB[j]) * (eIslandBCEB[i] + eIslandBCEB[j]) -
0565 (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
0566 cout << " pi0 (pt>2.5 GeV) m_inv = " << m_inv << endl;
0567
0568 float s4s9_1 = e2x2IslandBCEB[i] / e3x3IslandBCEB[i];
0569 float s4s9_2 = e2x2IslandBCEB[j] / e3x3IslandBCEB[j];
0570
0571 float s9s25_1 = e3x3IslandBCEB[i] / e5x5IslandBCEB[i];
0572 float s9s25_2 = e3x3IslandBCEB[j] / e5x5IslandBCEB[j];
0573
0574
0575
0576
0577 if (s4s9_1 > selePi0S4S9GammaOneMin_ && s4s9_2 > selePi0S4S9GammaTwoMin_ &&
0578 s9s25_1 > selePi0S9S25GammaOneMin_ && s9s25_2 > selePi0S9S25GammaTwoMin_ &&
0579 (et_belt / pt_pi0 < selePi0EtBeltIsoRatioMax_)) {
0580
0581 if (m_inv > (selePi0MinvMeanFixed_ - 2 * selePi0MinvSigmaFixed_) &&
0582 m_inv < (selePi0MinvMeanFixed_ + 2 * selePi0MinvSigmaFixed_)) {
0583
0584 cout << " Pi0 Good candidate : minv = " << m_inv << endl;
0585 for (int kk = 0; kk < nIslandBCEBRecHits[i]; kk++) {
0586 int ieta_xtal = ietaIslandBCEBRecHits[i][kk];
0587 int iphi_xtal = iphiIslandBCEBRecHits[i][kk];
0588 int sign = zsideIslandBCEBRecHits[i][kk] > 0 ? 1 : 0;
0589 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0590 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] + eIslandBCEBRecHits[i][kk] / e3x3IslandBCEB[i];
0591 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0592 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] +
0593 (selePi0MinvMeanFixed_ / m_inv) * (selePi0MinvMeanFixed_ / m_inv) *
0594 (eIslandBCEBRecHits[i][kk] / e3x3IslandBCEB[i]);
0595 cout << "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: " << ieta_xtal
0596 << " " << iphi_xtal << " " << sign << " " << eIslandBCEBRecHits[i][kk] << " "
0597 << e3x3IslandBCEB[i] << " " << wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << " "
0598 << mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << endl;
0599 }
0600
0601 for (int kk = 0; kk < nIslandBCEBRecHits[j]; kk++) {
0602 int ieta_xtal = ietaIslandBCEBRecHits[j][kk];
0603 int iphi_xtal = iphiIslandBCEBRecHits[j][kk];
0604 int sign = zsideIslandBCEBRecHits[j][kk] > 0 ? 1 : 0;
0605 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0606 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] + eIslandBCEBRecHits[j][kk] / e3x3IslandBCEB[j];
0607 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0608 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] +
0609 (selePi0MinvMeanFixed_ / m_inv) * (selePi0MinvMeanFixed_ / m_inv) *
0610 (eIslandBCEBRecHits[j][kk] / e3x3IslandBCEB[j]);
0611 cout << "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: " << ieta_xtal
0612 << " " << iphi_xtal << " " << sign << " " << eIslandBCEBRecHits[j][kk] << " "
0613 << e3x3IslandBCEB[j] << " " << wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << " "
0614 << mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << endl;
0615 }
0616 }
0617 }
0618 }
0619 }
0620
0621 }
0622 }
0623
0624 } else {
0625 cout << " Not enough ECAL Barrel Basic Clusters: " << nIslandBCEB << endl;
0626 }
0627
0628 return kContinue;
0629 }
0630
0631