File indexing completed on 2024-09-07 04:34:58
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 Handle<EcalRecHitCollection> pEcalRecHitBarrelCollection;
0331 event.getByToken(recHitToken_, pEcalRecHitBarrelCollection);
0332 const EcalRecHitCollection* ecalRecHitBarrelCollection = pEcalRecHitBarrelCollection.product();
0333 cout << " ECAL Barrel RecHits # " << ecalRecHitBarrelCollection->size() << endl;
0334 for (EcalRecHitCollection::const_iterator aRecHitEB = ecalRecHitBarrelCollection->begin();
0335 aRecHitEB != ecalRecHitBarrelCollection->end();
0336 aRecHitEB++) {
0337 EBDetId ebrhdetid = aRecHitEB->detid();
0338
0339
0340
0341
0342 int sign = ebrhdetid.zside() > 0 ? 1 : 0;
0343 EcalRecHit aHit(aRecHitEB->id(),
0344 aRecHitEB->energy() * oldCalibs_barl[abs(ebrhdetid.ieta()) - 1][ebrhdetid.iphi() - 1][sign],
0345 aRecHitEB->time());
0346 recalibEcalRecHitCollection->push_back(aHit);
0347 }
0348
0349
0350 for (EcalRecHitCollection::const_iterator aRecHitEB = recalibEcalRecHitCollection->begin();
0351 aRecHitEB != recalibEcalRecHitCollection->end();
0352 aRecHitEB++) {
0353
0354
0355
0356
0357
0358 std::pair<DetId, EcalRecHit> map_entry(aRecHitEB->id(), *aRecHitEB);
0359 recHitsEB_map->insert(map_entry);
0360 }
0361
0362 const CaloSubdetectorGeometry* geometry_p;
0363
0364 std::string clustershapetag;
0365 geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0366 EcalBarrelTopology topology{geometry};
0367
0368 const CaloSubdetectorGeometry* geometryES_p;
0369 geometryES_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397 reco::BasicClusterCollection clusters_recalib;
0398 clusters_recalib = island_p->makeClusters(
0399 recalibEcalRecHitCollection, geometry_p, &topology, geometryES_p, IslandClusterAlgo::barrel);
0400
0401
0402 std::vector<reco::ClusterShape> ClusVec_recalib;
0403 for (int erg = 0; erg < int(clusters_recalib.size()); ++erg) {
0404 reco::ClusterShape TestShape_recalib =
0405 shapeAlgo_.Calculate(clusters_recalib[erg], recalibEcalRecHitCollection, geometry_p, &topology);
0406 ClusVec_recalib.push_back(TestShape_recalib);
0407 }
0408
0409
0410 std::unique_ptr<reco::ClusterShapeCollection> clustersshapes_p_recalib(new reco::ClusterShapeCollection);
0411 clustersshapes_p_recalib->assign(ClusVec_recalib.begin(), ClusVec_recalib.end());
0412
0413 cout << "[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster collection size: " << clusters_recalib.size()
0414 << endl;
0415 cout << "[Pi0FixedMassWindowCalibration][recalibration] Basic Cluster Shape Collection size: "
0416 << clustersshapes_p_recalib->size() << endl;
0417
0418
0419
0420
0421
0422 static const int MAXBCEB = 200;
0423 static const int MAXBCEBRH = 200;
0424 int nIslandBCEB;
0425 float eIslandBCEB[MAXBCEB];
0426 float etIslandBCEB[MAXBCEB];
0427 float etaIslandBCEB[MAXBCEB];
0428 float phiIslandBCEB[MAXBCEB];
0429 float e2x2IslandBCEB[MAXBCEB];
0430 float e3x3IslandBCEB[MAXBCEB];
0431 float e5x5IslandBCEB[MAXBCEB];
0432
0433
0434 int nIslandBCEBRecHits[MAXBCEB];
0435
0436 int ietaIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0437 int iphiIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0438 int zsideIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0439 float eIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
0440
0441 nIslandBCEB = 0;
0442 for (int i = 0; i < MAXBCEB; i++) {
0443 eIslandBCEB[i] = 0;
0444 etIslandBCEB[i] = 0;
0445 etaIslandBCEB[i] = 0;
0446 phiIslandBCEB[i] = 0;
0447 e2x2IslandBCEB[i] = 0;
0448 e3x3IslandBCEB[i] = 0;
0449 e5x5IslandBCEB[i] = 0;
0450 nIslandBCEBRecHits[i] = 0;
0451 for (int j = 0; j < MAXBCEBRH; j++) {
0452
0453 ietaIslandBCEBRecHits[i][j] = 0;
0454 iphiIslandBCEBRecHits[i][j] = 0;
0455 zsideIslandBCEBRecHits[i][j] = 0;
0456 eIslandBCEBRecHits[i][j] = 0;
0457 }
0458 }
0459
0460 int iClus_recalib = 0;
0461 for (reco::BasicClusterCollection::const_iterator aClus = clusters_recalib.begin(); aClus != clusters_recalib.end();
0462 aClus++) {
0463 cout << " CLUSTER [recalibration] : #,NHits,e,et,eta,phi,e2x2,e3x3,e5x5: " << iClus_recalib << " " << aClus->size()
0464 << " " << aClus->energy() << " " << aClus->energy() * sin(aClus->position().theta()) << " "
0465 << aClus->position().eta() << " " << aClus->position().phi() << " "
0466 << (*clustersshapes_p_recalib)[iClus_recalib].e2x2() << " "
0467 << (*clustersshapes_p_recalib)[iClus_recalib].e3x3() << " "
0468 << (*clustersshapes_p_recalib)[iClus_recalib].e5x5() << endl;
0469
0470 eIslandBCEB[nIslandBCEB] = aClus->energy();
0471 etIslandBCEB[nIslandBCEB] = aClus->energy() * sin(aClus->position().theta());
0472 etaIslandBCEB[nIslandBCEB] = aClus->position().eta();
0473 phiIslandBCEB[nIslandBCEB] = aClus->position().phi();
0474
0475 e2x2IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e2x2();
0476 e3x3IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e3x3();
0477 e5x5IslandBCEB[nIslandBCEB] = (*clustersshapes_p_recalib)[nIslandBCEB].e5x5();
0478
0479 nIslandBCEBRecHits[nIslandBCEB] = aClus->size();
0480
0481 std::vector<std::pair<DetId, float>> hits = aClus->hitsAndFractions();
0482 std::vector<std::pair<DetId, float>>::iterator hit;
0483 std::map<DetId, EcalRecHit>::iterator aHit;
0484
0485 int irhcount = 0;
0486 for (hit = hits.begin(); hit != hits.end(); hit++) {
0487
0488 aHit = recHitsEB_map->find((*hit).first);
0489
0490
0491 EBDetId sel_rh = aHit->second.detid();
0492
0493
0494
0495
0496 ietaIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.ieta();
0497 iphiIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.iphi();
0498 zsideIslandBCEBRecHits[nIslandBCEB][irhcount] = sel_rh.zside();
0499 eIslandBCEBRecHits[nIslandBCEB][irhcount] = aHit->second.energy();
0500
0501 irhcount++;
0502 }
0503 nIslandBCEB++;
0504 iClus_recalib++;
0505 }
0506
0507
0508
0509 if (nIslandBCEB > 1) {
0510 for (int i = 0; i < nIslandBCEB; i++) {
0511 for (int j = i + 1; j < nIslandBCEB; j++) {
0512 if (etIslandBCEB[i] > selePi0PtGammaOneMin_ && etIslandBCEB[j] > selePi0PtGammaOneMin_) {
0513 float theta_0 = 2. * atan(exp(-etaIslandBCEB[i]));
0514 float theta_1 = 2. * atan(exp(-etaIslandBCEB[j]));
0515
0516 float p0x = eIslandBCEB[i] * sin(theta_0) * cos(phiIslandBCEB[i]);
0517 float p1x = eIslandBCEB[j] * sin(theta_1) * cos(phiIslandBCEB[j]);
0518
0519 float p0y = eIslandBCEB[i] * sin(theta_0) * sin(phiIslandBCEB[i]);
0520 float p1y = eIslandBCEB[j] * sin(theta_1) * sin(phiIslandBCEB[j]);
0521
0522 float p0z = eIslandBCEB[i] * cos(theta_0);
0523 float p1z = eIslandBCEB[j] * cos(theta_1);
0524
0525 float pi0_px = p0x + p1x;
0526 float pi0_py = p0y + p1y;
0527 float pi0_pz = p0z + p1z;
0528
0529 float pi0_ptot = sqrt(pi0_px * pi0_px + pi0_py * pi0_py + pi0_pz * pi0_pz);
0530
0531 float pi0_theta = acos(pi0_pz / pi0_ptot);
0532 float pi0_eta = -log(tan(pi0_theta / 2));
0533 float pi0_phi = atan(pi0_py / pi0_px);
0534
0535
0536
0537
0538 float et_belt = 0;
0539 for (Int_t k = 0; k < nIslandBCEB; k++) {
0540 if ((k != i) && (k != j)) {
0541 float dr_pi0_k = sqrt((etaIslandBCEB[k] - pi0_eta) * (etaIslandBCEB[k] - pi0_eta) +
0542 (phiIslandBCEB[k] - pi0_phi) * (phiIslandBCEB[k] - pi0_phi));
0543 float deta_pi0_k = fabs(etaIslandBCEB[k] - pi0_eta);
0544 if ((dr_pi0_k < selePi0DRBelt_) && (deta_pi0_k < selePi0DetaBelt_))
0545 et_belt = et_belt + etIslandBCEB[k];
0546 }
0547 }
0548
0549 float pt_pi0 = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
0550
0551
0552
0553 if (pt_pi0 > selePi0PtPi0Min_) {
0554 float m_inv = sqrt((eIslandBCEB[i] + eIslandBCEB[j]) * (eIslandBCEB[i] + eIslandBCEB[j]) -
0555 (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
0556 cout << " pi0 (pt>2.5 GeV) m_inv = " << m_inv << endl;
0557
0558 float s4s9_1 = e2x2IslandBCEB[i] / e3x3IslandBCEB[i];
0559 float s4s9_2 = e2x2IslandBCEB[j] / e3x3IslandBCEB[j];
0560
0561 float s9s25_1 = e3x3IslandBCEB[i] / e5x5IslandBCEB[i];
0562 float s9s25_2 = e3x3IslandBCEB[j] / e5x5IslandBCEB[j];
0563
0564
0565
0566
0567 if (s4s9_1 > selePi0S4S9GammaOneMin_ && s4s9_2 > selePi0S4S9GammaTwoMin_ &&
0568 s9s25_1 > selePi0S9S25GammaOneMin_ && s9s25_2 > selePi0S9S25GammaTwoMin_ &&
0569 (et_belt / pt_pi0 < selePi0EtBeltIsoRatioMax_)) {
0570
0571 if (m_inv > (selePi0MinvMeanFixed_ - 2 * selePi0MinvSigmaFixed_) &&
0572 m_inv < (selePi0MinvMeanFixed_ + 2 * selePi0MinvSigmaFixed_)) {
0573
0574 cout << " Pi0 Good candidate : minv = " << m_inv << endl;
0575 for (int kk = 0; kk < nIslandBCEBRecHits[i]; kk++) {
0576 int ieta_xtal = ietaIslandBCEBRecHits[i][kk];
0577 int iphi_xtal = iphiIslandBCEBRecHits[i][kk];
0578 int sign = zsideIslandBCEBRecHits[i][kk] > 0 ? 1 : 0;
0579 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0580 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] + eIslandBCEBRecHits[i][kk] / e3x3IslandBCEB[i];
0581 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0582 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] +
0583 (selePi0MinvMeanFixed_ / m_inv) * (selePi0MinvMeanFixed_ / m_inv) *
0584 (eIslandBCEBRecHits[i][kk] / e3x3IslandBCEB[i]);
0585 cout << "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: " << ieta_xtal
0586 << " " << iphi_xtal << " " << sign << " " << eIslandBCEBRecHits[i][kk] << " "
0587 << e3x3IslandBCEB[i] << " " << wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << " "
0588 << mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << endl;
0589 }
0590
0591 for (int kk = 0; kk < nIslandBCEBRecHits[j]; kk++) {
0592 int ieta_xtal = ietaIslandBCEBRecHits[j][kk];
0593 int iphi_xtal = iphiIslandBCEBRecHits[j][kk];
0594 int sign = zsideIslandBCEBRecHits[j][kk] > 0 ? 1 : 0;
0595 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0596 wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] + eIslandBCEBRecHits[j][kk] / e3x3IslandBCEB[j];
0597 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] =
0598 mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] +
0599 (selePi0MinvMeanFixed_ / m_inv) * (selePi0MinvMeanFixed_ / m_inv) *
0600 (eIslandBCEBRecHits[j][kk] / e3x3IslandBCEB[j]);
0601 cout << "[Pi0FixedMassWindowCalibration] eta, phi, sign, e, e3x3, wxtals and mwxtals: " << ieta_xtal
0602 << " " << iphi_xtal << " " << sign << " " << eIslandBCEBRecHits[j][kk] << " "
0603 << e3x3IslandBCEB[j] << " " << wxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << " "
0604 << mwxtals[abs(ieta_xtal) - 1][iphi_xtal - 1][sign] << endl;
0605 }
0606 }
0607 }
0608 }
0609 }
0610
0611 }
0612 }
0613
0614 } else {
0615 cout << " Not enough ECAL Barrel Basic Clusters: " << nIslandBCEB << endl;
0616 }
0617
0618 return kContinue;
0619 }
0620
0621