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   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     //cout << " ECAL Barrel RecHit #,E,time,det,subdetid: "<<nRecHitsEB<<" "<<aRecHitEB->energy()<<" "<<aRecHitEB->time()<<" "<<aRecHitEB->detid().det()<<" "<<aRecHitEB->detid().subdetId()<<endl;
0339
0340     EBDetId ebrhdetid = aRecHitEB->detid();
0341     //cout << " EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
0342     //cout << " EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
0343     //cout << " EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;
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   //  cout<<" Recalib size: "<<recalibEcalRecHitCollection->size()<<endl;
0355   int irecalib = 0;
0356   for (EcalRecHitCollection::const_iterator aRecHitEB = recalibEcalRecHitCollection->begin();
0357        aRecHitEB != recalibEcalRecHitCollection->end();
0358        aRecHitEB++) {
0359     //cout << " [recalibrated] ECAL Barrel RecHit #,E,time,det,subdetid: "<<irecalib<<" "<<aRecHitEB->energy()<<" "<<aRecHitEB->time()<<" "<<aRecHitEB->detid().det()<<" "<<aRecHitEB->detid().subdetId()<<endl;
0360
0361     //    EBDetId ebrhdetid = aRecHitEB->detid();
0362     //cout << " [recalibrated] EBDETID: z,ieta,iphi "<<ebrhdetid.zside()<<" "<<ebrhdetid.ieta()<<" "<<ebrhdetid.iphi()<<endl;
0363     //cout << " [recalibrated] EBDETID: tower_ieta,tower_iphi "<<ebrhdetid.tower_ieta()<<" "<<ebrhdetid.tower_iphi()<<endl;
0364     //cout << " [recalibrated] EBDETID: iSM, ic "<<ebrhdetid.ism()<<" "<<ebrhdetid.ic()<<endl;
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   reco::BasicClusterCollection clusters;
0383   clusters = island_p->makeClusters(ecalRecHitBarrelCollection,geometry_p,&topology,geometryES_p,IslandClusterAlgo::barrel);
0384
0385   //Create associated ClusterShape objects.
0386   std::vector <reco::ClusterShape> ClusVec;
0387   for (int erg=0;erg<int(clusters.size());++erg){
0388     reco::ClusterShape TestShape = shapeAlgo_.Calculate(clusters[erg],ecalRecHitBarrelCollection,geometry_p,&topology);
0389     ClusVec.push_back(TestShape);
0390   }
0391
0392   //Put clustershapes in event, but retain a Handle on them.
0393   std::unique_ptr<reco::ClusterShapeCollection> clustersshapes_p(new reco::ClusterShapeCollection);
0394   clustersshapes_p->assign(ClusVec.begin(), ClusVec.end());
0395
0396   cout<<"[Pi0Calibration] Basic Cluster collection size: "<<clusters.size()<<endl;
0397   cout<<"[Pi0Calibration] Basic Cluster Shape Collection size: "<<clustersshapes_p->size()<<endl;
0398
0399   int iClus=0;
0400   for(reco::BasicClusterCollection::const_iterator aClus = clusters.begin(); aClus != clusters.end(); aClus++) {
0401     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;
0402     iClus++;
0403   }
0404   */
0405
0406   // recalibrated clusters
0407   reco::BasicClusterCollection clusters_recalib;
0408   clusters_recalib = island_p->makeClusters(
0409       recalibEcalRecHitCollection, geometry_p, &topology, geometryES_p, IslandClusterAlgo::barrel);
0410
0411   //Create associated ClusterShape objects.
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   //Put clustershapes in event, but retain a Handle on them.
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   // pizero selection
0429
0430   // Get ECAL Barrel Island Basic Clusters collection
0431   // ECAL Barrel Island Basic Clusters
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   // indexes to the RecHits assiciated with
0443   // ECAL Barrel Island Basic Clusters
0444   int nIslandBCEBRecHits[MAXBCEB];
0445   //  int indexIslandBCEBRecHits[MAXBCEB][MAXBCEBRH];
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       //       indexIslandBCEBRecHits[i][j] = 0;
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       // need to get hit by DetID in order to get energy
0498       aHit = recHitsEB_map->find((*hit).first);
0499       //cout << "       RecHit #: "<<irhcount<<" from Basic Cluster with Energy: "<<aHit->second.energy()<<endl;
0500
0501       EBDetId sel_rh = aHit->second.detid();
0502       //cout << "       RecHit: z,ieta,iphi "<<sel_rh.zside()<<" "<<sel_rh.ieta()<<" "<<sel_rh.iphi()<<endl;
0503       //cout << "       RecHit: tower_ieta,tower_iphi "<<sel_rh.tower_ieta()<<" "<<sel_rh.tower_iphi()<<endl;
0504       //cout << "       RecHit: iSM, ic "<<sel_rh.ism()<<" "<<sel_rh.ic()<<endl;
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   // Selection, based on ECAL Barrel Basic Clusters
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           //cout << " pi0_theta, pi0_eta, pi0_phi "<<pi0_theta<<" "<<pi0_eta<<" "<<pi0_phi<<endl;
0545
0546           // belt isolation
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           //float dr_pi0 = sqrt ( (etaIslandBCEB[i]-etaIslandBCEB[j])*(etaIslandBCEB[i]-etaIslandBCEB[j]) + (phiIslandBCEB[i]-phiIslandBCEB[j])*(phiIslandBCEB[i]-phiIslandBCEB[j]) );
0561
0562           //cout <<" pi0 pt,dr:  "<<pt_pi0<<" "<<dr_pi0<<endl;
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             //float s9Esc_1 = e3x3IslandBCEB[i]/eIslandBCEB[i];
0575             //float s9Esc_2 = e3x3IslandBCEB[j]/eIslandBCEB[j];
0576
0577             if (s4s9_1 > selePi0S4S9GammaOneMin_ && s4s9_2 > selePi0S4S9GammaTwoMin_ &&
0578                 s9s25_1 > selePi0S9S25GammaOneMin_ && s9s25_2 > selePi0S9S25GammaTwoMin_ &&
0579                 (et_belt / pt_pi0 < selePi0EtBeltIsoRatioMax_)) {
0580               //good pizero candidate
0581               if (m_inv > (selePi0MinvMeanFixed_ - 2 * selePi0MinvSigmaFixed_) &&
0582                   m_inv < (selePi0MinvMeanFixed_ + 2 * selePi0MinvSigmaFixed_)) {
0583                 //fill wxtals and mwxtals weights
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       }  // End of the "j" loop over BCEB
0622     }    // End of the "i" loop over BCEB
0623
0624   } else {
0625     cout << " Not enough ECAL Barrel Basic Clusters: " << nIslandBCEB << endl;
0626   }
0627
0628   return kContinue;
0629 }
0630
0631 // ----------------------------------------------------------------------------