EcalCalibAlgos/src/Pi0FixedMassWindowCalibration.cc

0001 #include "Calibration/EcalCalibAlgos/interface/Pi0FixedMassWindowCalibration.h"
0002
0003 // System include files
0004
0005 // Framework
0006
0007 // Conditions database
0008
0009 #include "Calibration/Tools/interface/Pi0CalibXMLwriter.h"
0010
0011 // Reconstruction Classes
0012 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0013 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0014 // Geometry
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 // EgammaCoreTools
0021 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0022 #include "DataFormats/EgammaReco/interface/ClusterShapeFwd.h"
0023
0024 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0025
0026 //const double Pi0Calibration::PDGPi0Mass =  0.1349766;
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   // The verbosity level
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   // The names of the produced cluster collections
0052   barrelClusterCollection_ = iConfig.getParameter<std::string>("barrelClusterCollection");
0053
0054   // Island algorithm parameters
0055   double barrelSeedThreshold = iConfig.getParameter<double>("IslandBarrelSeedThr");
0056   double endcapSeedThreshold = iConfig.getParameter<double>("IslandEndcapSeedThr");
0057
0058   // Selection algorithm parameters
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   // Parameters for the position calculation:
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   //AssociationMap
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 // Close files, etc.
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 // Initialize algorithm
0172
0173 void Pi0FixedMassWindowCalibration::beginOfJob() {
0174   //std::cout << "[Pi0FixedMassWindowCalibration] beginOfJob "<<std::endl;
0175
0176   isfirstcall_ = true;
0177 }
0178
0179 void Pi0FixedMassWindowCalibration::endOfJob() {
0180   std::cout << "[Pi0FixedMassWindowCalibration] endOfJob" << endl;
0181
0182   // Write new calibration constants
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 // Called at beginning of loop
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 // Called at end of loop
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   // replace old calibration constants with new one
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 // Called at each event
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   // this chunk used to belong to beginJob(isetup). Moved here
0273   // with the beginJob without arguments migration
0274
0275   // get the ecal geometry:
0276   const auto& geometry = setup.getData(geometryToken_);
0277
0278   if (isfirstcall_) {
0279     // initialize arrays
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     // get initial constants out of DB
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     // loop over all barrel crystals
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       // get the initial calibration constants
0305       EcalIntercalibConstantMap::const_iterator itcalib = imap.find(eb.rawId());
0306       if (itcalib == imap.end()) {
0307         // FIXME -- throw error
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     //cout << " EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
0339     //cout << " EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
0340     //cout << " EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;
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   //  cout<<" Recalib size: "<<recalibEcalRecHitCollection->size()<<endl;
0350   for (EcalRecHitCollection::const_iterator aRecHitEB = recalibEcalRecHitCollection->begin();
0351        aRecHitEB != recalibEcalRecHitCollection->end();
0352        aRecHitEB++) {
0353     //    EBDetId ebrhdetid = aRecHitEB->detid();
0354     //cout << " [recalibrated] EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
0355     //cout << " [recalibrated] EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
0356     //cout << " [recalibrated] EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;
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   reco::BasicClusterCollection clusters;
0373   clusters = island_p->makeClusters(ecalRecHitBarrelCollection,geometry_p,&topology,geometryES_p,IslandClusterAlgo::barrel);
0374
0375   //Create associated ClusterShape objects.
0376   std::vector <reco::ClusterShape> ClusVec;
0377   for (int erg=0;erg<int(clusters.size());++erg){
0378     reco::ClusterShape TestShape = shapeAlgo_.Calculate(clusters[erg],ecalRecHitBarrelCollection,geometry_p,&topology);
0379     ClusVec.push_back(TestShape);
0380   }
0381
0382   //Put clustershapes in event, but retain a Handle on them.
0383   std::unique_ptr<reco::ClusterShapeCollection> clustersshapes_p(new reco::ClusterShapeCollection);
0384   clustersshapes_p->assign(ClusVec.begin(), ClusVec.end());
0385
0386   cout<<"[Pi0Calibration] Basic Cluster collection size: "<<clusters.size()<<endl;
0387   cout<<"[Pi0Calibration] Basic Cluster Shape Collection size: "<<clustersshapes_p->size()<<endl;
0388
0389   int iClus=0;
0390   for(reco::BasicClusterCollection::const_iterator aClus = clusters.begin(); aClus != clusters.end(); aClus++) {
0391     cout<<" CLUSTER : #,NHits,e,et,eta,phi,e2x2,e3x3,e5x5: "<<iClus<<" "<<aClus->getHitsByDetId().size()<<" "<<aClus->energy()<<" "<<aClus->energy()*sin(aClus->position().theta())<<" "<<aClus->position().eta()<<" "<<aClus->position().phi()<<" "<<(*clustersshapes_p)[iClus].e2x2()<<" "<<(*clustersshapes_p)[iClus].e3x3()<<" "<<(*clustersshapes_p)[iClus].e5x5()<<endl;
0392     iClus++;
0393   }
0394   */
0395
0396   // recalibrated clusters
0397   reco::BasicClusterCollection clusters_recalib;
0398   clusters_recalib = island_p->makeClusters(
0399       recalibEcalRecHitCollection, geometry_p, &topology, geometryES_p, IslandClusterAlgo::barrel);
0400
0401   //Create associated ClusterShape objects.
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   //Put clustershapes in event, but retain a Handle on them.
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   // pizero selection
0419
0420   // Get ECAL Barrel Island Basic Clusters collection
0421   // ECAL Barrel Island Basic Clusters
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   // indexes to the RecHits assiciated with
0433   // ECAL Barrel Island Basic Clusters
0434   int nIslandBCEBRecHits[MAXBCEB];
0435   //  int indexIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
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       //       indexIslandBCEBRecHits[i][j] = 0;
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       // need to get hit by DetID in order to get energy
0488       aHit = recHitsEB_map->find((*hit).first);
0489       //cout << "       RecHit #: "<<irhcount<<" from Basic Cluster with Energy: "<<aHit->second.energy()<<endl;
0490
0491       EBDetId sel_rh = aHit->second.detid();
0492       //cout << "       RecHit: z,ieta,iphi "<<sel_rh.zside()<<" "<<sel_rh.ieta()<<" "<<sel_rh.iphi()<<endl;
0493       //cout << "       RecHit: tower_ieta,tower_iphi "<<sel_rh.tower_ieta()<<" "<<sel_rh.tower_iphi()<<endl;
0494       //cout << "       RecHit: iSM, ic "<<sel_rh.ism()<<" "<<sel_rh.ic()<<endl;
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   // Selection, based on ECAL Barrel Basic Clusters
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           //cout << " pi0_theta, pi0_eta, pi0_phi "<<pi0_theta<<" "<<pi0_eta<<" "<<pi0_phi<<endl;
0535
0536           // belt isolation
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           //float dr_pi0 = sqrt ( (etaIslandBCEB[i]-etaIslandBCEB[j])*(etaIslandBCEB[i]-etaIslandBCEB[j]) + (phiIslandBCEB[i]-phiIslandBCEB[j])*(phiIslandBCEB[i]-phiIslandBCEB[j]) );
0551
0552           //cout <<" pi0 pt,dr:  "<<pt_pi0<<" "<<dr_pi0<<endl;
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             //float s9Esc_1 = e3x3IslandBCEB[i]/eIslandBCEB[i];
0565             //float s9Esc_2 = e3x3IslandBCEB[j]/eIslandBCEB[j];
0566
0567             if (s4s9_1 > selePi0S4S9GammaOneMin_ && s4s9_2 > selePi0S4S9GammaTwoMin_ &&
0568                 s9s25_1 > selePi0S9S25GammaOneMin_ && s9s25_2 > selePi0S9S25GammaTwoMin_ &&
0569                 (et_belt / pt_pi0 < selePi0EtBeltIsoRatioMax_)) {
0570               //good pizero candidate
0571               if (m_inv > (selePi0MinvMeanFixed_ - 2 * selePi0MinvSigmaFixed_) &&
0572                   m_inv < (selePi0MinvMeanFixed_ + 2 * selePi0MinvSigmaFixed_)) {
0573                 //fill wxtals and mwxtals weights
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       }  // End of the "j" loop over BCEB
0612     }    // End of the "i" loop over BCEB
0613
0614   } else {
0615     cout << " Not enough ECAL Barrel Basic Clusters: " << nIslandBCEB << endl;
0616   }
0617
0618   return kContinue;
0619 }
0620
0621 // ----------------------------------------------------------------------------