11 #include <CLHEP/Vector/LorentzVector.h> 41 : hlTriggerResults_(iConfig.getParameter<
edm::
InputTag>(
"HLTriggerResults")),
42 mcProducer_(iConfig.getUntrackedParameter<
std::
string>(
"mcProducer",
"")),
43 rechitProducer_(iConfig.getParameter<
std::
string>(
"rechitProducer")),
44 rechitCollection_(iConfig.getParameter<
std::
string>(
"rechitCollection")),
45 erechitProducer_(iConfig.getParameter<
std::
string>(
"erechitProducer")),
46 erechitCollection_(iConfig.getParameter<
std::
string>(
"erechitCollection")),
47 scProducer_(iConfig.getParameter<
std::
string>(
"scProducer")),
48 scCollection_(iConfig.getParameter<
std::
string>(
"scCollection")),
49 scIslandProducer_(iConfig.getParameter<
std::
string>(
"scIslandProducer")),
50 scIslandCollection_(iConfig.getParameter<
std::
string>(
"scIslandCollection")),
51 electronProducer_(iConfig.getParameter<
std::
string>(
"electronProducer")),
52 electronCollection_(iConfig.getParameter<
std::
string>(
"electronCollection")),
62 std::cout <<
"[ZeeCalibration] Starting the ctor" << std::endl;
78 myTree =
new TTree(
"myTree",
"myTree");
116 findingRecord<EcalIntercalibConstantsRcd>();
118 for (
int i = 0;
i < 50;
i++) {
126 std::cout <<
"[ZeeCalibration] Done with the ctor" << std::endl;
139 std::cout <<
"@SUB=ZeeCalibration::produceEcalIntercalibConstants" << std::endl;
183 std::cout <<
"Writing histos..." << std::endl;
212 for (
unsigned int i = 0;
i < 25;
i++) {
247 Double_t
mean[25] = {0.};
248 Double_t
num[25] = {0.};
249 Double_t meanErr[25] = {0.};
250 Float_t
rms[25] = {0.};
251 Float_t tempRms[10][25];
253 for (
int ia = 0; ia < 10; ia++) {
254 for (
int ib = 0;
ib < 25;
ib++) {
255 tempRms[ia][
ib] = 0.;
263 bool isNearCrack =
false;
274 if ((
k + 1) % 5 != 0) {
318 if (
k >= 170 &&
k <= 204) {
335 else if (
flag == 4) {
354 if (
k >= 205 &&
k <= 208) {
396 for (
int ic = 0; ic < 17; ic++) {
399 meanErr[ic] = 1. /
sqrt(meanErr[ic]);
403 for (
int ic = 0; ic < 25; ic++) {
404 for (
int id = 0;
id < 10;
id++) {
405 if (tempRms[
id][ic] > 0.) {
415 Double_t xtalEta[25] = {1.4425, 1.3567, 1.2711, 1.1855, 1.10, 1.01, 0.92, 0.83, 0.7468,
416 0.6612, 0.5756, 0.4897, 0.3985, 0.3117, 0.2250, 0.1384, 0.0487, 1.546,
417 1.651, 1.771, 1.908, 2.071, 2.267, 2.516, 2.8};
419 Double_t
zero[25] = {0.026};
421 for (
int j = 0;
j < 25;
j++)
428 px->SetXTitle(
"Eta channel");
429 px->SetYTitle(
"recalibCoeff");
449 double weightSumMeanBarrel = 0.;
450 double weightSumMeanEndcap = 0.;
484 std::cout <<
"Weight sum mean on channels in Barrel is :" << weightSumMeanBarrel << std::endl;
485 std::cout <<
"Weight sum mean on channels in Endcap is :" << weightSumMeanEndcap << std::endl;
497 TGraphErrors* graph =
new TGraphErrors(25, xtalEta,
mean,
zero, meanErr);
501 double zero50[50] = {0.};
504 residualSigmaGraph->SetName(
"residualSigmaGraph");
505 residualSigmaGraph->Draw(
"APL");
506 residualSigmaGraph->Write();
508 TGraphErrors* coefficientDistanceAtIterationGraph =
510 coefficientDistanceAtIterationGraph->SetName(
"coefficientDistanceAtIterationGraph");
511 coefficientDistanceAtIterationGraph->Draw(
"APL");
512 coefficientDistanceAtIterationGraph->Write();
514 Float_t noError[250] = {0.};
516 Float_t ringInd[250];
517 for (
int i = 0;
i < 250;
i++)
520 TGraphErrors* graphCoeff =
522 graphCoeff->SetName(
"graphCoeff");
523 graphCoeff->Draw(
"APL");
556 std::cout <<
"[ZeeCalibration] Entering duringLoop" << std::endl;
572 std::cout <<
"[ZeeCalibration::beginOfJob] Histograms booked " << std::endl;
587 std::cout <<
"[ZeeCalibration::beginOfJob] Parsed EB miscal file" << std::endl;
595 std::cout <<
"[ZeeCalibration::beginOfJob] Parsed EE miscal file" << std::endl;
606 std::vector<DetId> ringIds;
619 for (
unsigned int iid = 0; iid < ringIds.size(); ++iid) {
620 float miscalib = *(miscalibMap->
get().
getMap().
find(ringIds[iid]));
631 ical = std::make_shared<EcalIntercalibConstants>();
636 std::vector<DetId> ringIds;
647 for (
unsigned int iid = 0; iid < ringIds.size(); ++iid) {
651 EBDetId myEBDetId(ringIds[iid]);
653 myEBDetId.
ieta() + 86,
659 EEDetId myEEDetId(ringIds[iid]);
660 if (myEEDetId.
zside() < 0)
666 if (myEEDetId.
zside() > 0)
684 for (
unsigned int iHLT = 0; iHLT < 200; ++iHLT) {
689 std::cout <<
"[ZeeCalibration::duringLoop] Done with initializing aHLTresults[] " << std::endl;
695 if (!hltTriggerResultHandle.
isValid()) {
704 std::cout <<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillHLTInfo(hltTriggerResultHandle); " 721 std::cout <<
"[ZeeCalibration::duringLoop] End HLT section" << std::endl;
726 std::vector<HepMC::GenParticle*> mcEle;
728 float myGenZMass(-1);
741 std::cout <<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillZMCInfo( & (*myGenEvent) ); " << std::endl;
744 for (HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
p != myGenEvent->particles_end();
747 if ((*p)->pdg_id() == 23 && (*p)->status() == 2) {
748 myGenZMass = (*p)->momentum().m();
759 for (HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
p != myGenEvent->particles_end();
761 if (
abs((*p)->pdg_id()) == 11) {
762 mcEle.push_back((*
p));
767 if (mcEle.size() == 2 && fabs(mcEle[0]->momentum().
eta()) < 2.4 && fabs(mcEle[1]->momentum().
eta()) < 2.4) {
770 if (fabs(mcEle[0]->momentum().
eta()) < 1.479 && fabs(mcEle[1]->momentum().
eta()) < 1.479)
773 if ((fabs(mcEle[0]->momentum().
eta()) > 1.479 && fabs(mcEle[1]->momentum().
eta()) < 1.479) ||
774 (fabs(mcEle[0]->momentum().
eta()) < 1.479 && fabs(mcEle[1]->momentum().
eta()) > 1.479))
777 if (fabs(mcEle[0]->momentum().eta()) > 1.479 && fabs(mcEle[1]->momentum().eta()) > 1.479)
785 if (!
phits.isValid()) {
786 std::cerr <<
"Error! can't get the product EBRecHitCollection " << std::endl;
794 std::cerr <<
"Error! can't get the product EERecHitCollection " << std::endl;
801 if (!pSuperClusters.
isValid()) {
802 std::cerr <<
"Error! can't get the product SuperClusterCollection " << std::endl;
808 for (reco::SuperClusterCollection::const_iterator scIt =
scCollection->begin(); scIt !=
scCollection->end(); scIt++) {
809 std::cout << scIt->energy() << std::endl;
816 if (!pIslandSuperClusters.
isValid()) {
817 std::cerr <<
"Error! can't get the product IslandSuperClusterCollection " << std::endl;
822 std::cout <<
"scCollection->size()" << scIslandCollection->size() << std::endl;
825 if ((
scCollection->size() + scIslandCollection->size()) < 2)
832 std::cerr <<
"Error! can't get the product ElectronCollection " << std::endl;
848 if (!
hits && !ehits) {
854 std::cout <<
"hits->size() == 0" << std::endl;
859 std::cout <<
"!electronCollection" << std::endl;
864 std::cout <<
"electronCollection->size() == 0" << std::endl;
879 std::cout <<
" Starting with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
886 std::cout <<
" Done with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
892 std::vector<calib::CalibElectron> calibElectrons;
899 std::cout << calibElectrons.back().getRecoElectron()->superCluster()->energy() <<
" " 900 << calibElectrons.back().getRecoElectron()->energy() << std::endl;
908 std::cout <<
"Filled histos" << std::endl;
912 std::vector<std::pair<calib::CalibElectron*, calib::CalibElectron*> > zeeCandidates;
916 double DeltaMinvMin(5000.);
918 if (calibElectrons.size() < 2)
921 for (
unsigned int e_it = 0; e_it != calibElectrons.size() - 1; e_it++) {
922 for (
unsigned int p_it = e_it + 1; p_it != calibElectrons.size(); p_it++) {
924 std::cout << e_it <<
" " << calibElectrons[e_it].getRecoElectron()->charge() <<
" " << p_it <<
" " 925 << calibElectrons[p_it].getRecoElectron()->charge() << std::endl;
927 if (calibElectrons[e_it].getRecoElectron()->
charge() * calibElectrons[p_it].getRecoElectron()->
charge() != -1)
931 std::pair<calib::CalibElectron*, calib::CalibElectron*>(&(calibElectrons[e_it]), &(calibElectrons[p_it])));
937 std::cout <<
"#######################mass " <<
mass << std::endl;
940 zeeCandidates.push_back(
941 std::pair<calib::CalibElectron*, calib::CalibElectron*>(&(calibElectrons[e_it]), &(calibElectrons[p_it])));
942 double DeltaMinv = fabs(
mass -
MZ);
944 if (DeltaMinv < DeltaMinvMin) {
945 DeltaMinvMin = DeltaMinv;
946 myBestZ = zeeCandidates.size() - 1;
957 if (zeeCandidates.empty() || myBestZ == -1)
964 std::cout <<
"Found ZCandidates " << myBestZ << std::endl;
971 int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
972 int class2 = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1001 if (class1 == 10 || class1 == 20)
1003 if (class1 == 110 || class1 == 120)
1005 if (class1 >= 30 && class1 <= 34)
1007 if (class1 >= 130 && class1 <= 134)
1018 if (class2 == 10 || class2 == 20)
1020 if (class2 == 110 || class2 == 120)
1022 if (class2 >= 30 && class2 <= 34)
1024 if (class2 >= 130 && class2 <= 134)
1035 DetId firstElehottestDetId =
1037 zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->seed()->hitsAndFractions(),
hits, ehits)
1039 DetId secondElehottestDetId =
1041 zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->seed()->hitsAndFractions(),
hits, ehits)
1044 bool firstElectronIsOnModuleBorder(
false);
1045 bool secondElectronIsOnModuleBorder(
false);
1060 if (secondElehottestDetId.subdetId() ==
EcalBarrel)
1065 if (secondElehottestDetId.subdetId() ==
EcalBarrel && !secondElectronIsOnModuleBorder)
1070 if (firstElehottestDetId.
subdetId() ==
EcalBarrel && firstElectronIsOnModuleBorder) {
1077 if (secondElehottestDetId.subdetId() ==
EcalBarrel && secondElectronIsOnModuleBorder) {
1086 if (class1 < 100 && class2 < 100) {
1088 if (class1 == 0 && class2 == 0)
1090 if (class1 < 21 && class2 < 21)
1092 if (class1 < 21 || class2 < 21)
1096 if (class1 >= 100 && class2 >= 100) {
1098 if (class1 == 100 && class2 == 100)
1100 if (class1 < 121 && class2 < 121)
1102 if (class1 < 121 || class2 < 121)
1106 if ((class1 < 100 && class2 >= 100) || (class2 < 100 && class1 >= 100)) {
1108 if ((class1 == 0 && class2 == 100) || (class2 == 0 && class1 == 100))
1110 if ((class1 < 21 && class2 < 121) || (class2 < 21 && class1 < 121))
1112 if (class2 < 21 || class1 < 21 || class2 < 121 || class1 < 121)
1123 bool selectionBool =
false;
1128 float theta1 = 2. * atan(
exp(-zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->eta()));
1129 bool ET_1 = ((zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->energy() *
sin(theta1)) > 20.);
1131 float theta2 = 2. * atan(
exp(-zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->eta()));
1132 bool ET_2 = ((zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->energy() *
sin(theta2)) > 20.);
1134 bool HoE_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->hadronicOverEm() < 0.115);
1135 bool HoE_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->hadronicOverEm() < 0.115);
1137 bool DeltaPhiIn_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1138 bool DeltaPhiIn_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1140 bool DeltaEtaIn_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1141 bool DeltaEtaIn_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1145 if (!(invMassBool && ET_1 && ET_2 && HoE_1 && HoE_2 && DeltaPhiIn_1 && DeltaPhiIn_2 && DeltaEtaIn_1 && DeltaEtaIn_2))
1153 int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1154 int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1156 selectionBool = (myBestZ != -1 && c1st != 40 && c1st != 40 && c2nd != 40 && c2nd != 140);
1162 (myBestZ != -1 && (c1st == 0 || c1st == 10 || c1st == 20 || c1st == 100 || c1st == 110 || c1st == 120) &&
1163 (c2nd == 0 || c2nd == 10 || c2nd == 20 || c2nd == 100 || c2nd == 110 || c2nd == 120));
1167 selectionBool = (myBestZ != -1 && (c1st == 0 || c1st == 100) && (c2nd == 0 || c2nd == 100));
1170 selectionBool = (myBestZ != -1 && ((c1st >= 30 && c1st <= 34) || ((c1st >= 130 && c1st <= 134))) &&
1171 ((c2nd >= 30 && c2nd <= 34) || ((c2nd >= 130 && c2nd <= 134)))
1178 selectionBool = (myBestZ != -1 && (
1180 ((c1st == 0 || c1st == 10 || c1st == 20) && c2nd >= 100 && c2nd != 140)
1184 ((c2nd == 0 || c2nd == 10 || c2nd == 20) && c1st >= 100 && c1st != 140)
1191 selectionBool = (myBestZ != -1 && c1st >= 100 && c2nd >= 100 && c1st != 140 && c2nd != 140);
1196 selectionBool = (myBestZ != -1 && c1st < 100 && c2nd < 100 && c1st != 40 && c2nd != 40);
1201 selectionBool = (myBestZ != -1 && !(c1st < 100 && c2nd >= 100) && !(c1st >= 100 && c2nd < 100));
1203 float ele1EnergyCorrection(1.);
1204 float ele2EnergyCorrection(1.);
1211 if (invMassBool && selectionBool) {
1219 zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) -
1227 std::map<HepMC::GenParticle*, const reco::GsfElectron*> myMCmap;
1229 std::vector<const reco::GsfElectron*> dauElectronCollection;
1231 dauElectronCollection.push_back(zeeCandidates[myBestZ].
first->getRecoElectron());
1232 dauElectronCollection.push_back(zeeCandidates[myBestZ].
second->getRecoElectron());
1234 fillMCmap(&dauElectronCollection, mcEle, myMCmap);
1241 zeeCandidates[myBestZ].
second,
1242 MZ *
sqrt(ele1EnergyCorrection * ele2EnergyCorrection));
1247 zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1249 int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1250 int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1251 if (c1st < 100 && c2nd < 100)
1253 zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1255 if (c1st >= 100 && c2nd >= 100)
1257 zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection));
1260 zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection);
1263 zeeCandidates[myBestZ], ele1EnergyCorrection, ele2EnergyCorrection) -
1272 std::cout <<
"Added event to algorithm" << std::endl;
1281 std::cout <<
"[ZeeCalibration] Starting loop number " << iLoop << std::endl;
1290 std::cout <<
"[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1311 std::cout <<
"[ZeeCalibration] Ending loop " << iLoop << std::endl;
1322 std::cout <<
"Optimized coefficients " << optimizedCoefficients.size() << std::endl;
1328 for (
unsigned int ieta = 0;
ieta < optimizedCoefficients.size();
ieta++) {
1341 std::cout <<
"size " << optimizedCoefficients.size() << std::endl;
1343 for (
unsigned int ieta = 0;
ieta < optimizedCoefficients.size();
ieta++) {
1354 std::vector<DetId> ringIds;
1365 for (
unsigned int iid = 0; iid < ringIds.size(); ++iid) {
1367 EBDetId myEBDetId(ringIds[iid]);
1369 myEBDetId.
ieta() + 86,
1375 EEDetId myEEDetId(ringIds[iid]);
1376 if (myEEDetId.
zside() < 0)
1382 if (myEEDetId.
zside() > 0)
1389 ical->setValue(ringIds[iid], *(
ical->getMap().find(ringIds[iid])) * optimizedCoefficients[
ieta]);
1419 double parResidual[3];
1420 double errparResidual[3];
1432 std::cout <<
"Fit on residuals, sigma is " << parResidual[2] <<
" +/- " << errparResidual[2] << std::endl;
1468 new TH1F(
"h1_eventsBeforeBorderSelection",
"h1_eventsBeforeBorderSelection", 5, 0, 5);
1471 h1_seedOverSC_ =
new TH1F(
"h1_seedOverSC",
"h1_seedOverSC", 400, 0., 2.);
1476 new TH1F(
"h1_borderElectronClassification",
"h1_borderElectronClassification", 55, -5, 50);
1479 h2_fEtaBarrelGood_ =
new TH2F(
"fEtaBarrelGood",
"fEtaBarrelGood", 800, -4., 4., 800, 0.8, 1.2);
1483 h2_fEtaBarrelBad_ =
new TH2F(
"fEtaBarrelBad",
"fEtaBarrelBad", 800, -4., 4., 800, 0.8, 1.2);
1487 h2_fEtaEndcapGood_ =
new TH2F(
"fEtaEndcapGood",
"fEtaEndcapGood", 800, -4., 4., 800, 0.8, 1.2);
1491 h2_fEtaEndcapBad_ =
new TH2F(
"fEtaEndcapBad",
"fEtaEndcapBad", 800, -4., 4., 800, 0.8, 1.2);
1495 for (
int i = 0;
i < 2;
i++) {
1512 new TH2F(
"h2_xtalMiscalibCoeffBarrel",
"h2_xtalMiscalibCoeffBarrel", 171, -85, 85, 360, 0, 360);
1514 new TH2F(
"h2_xtalMiscalibCoeffEndcapMinus",
"h2_xtalMiscalibCoeffEndcapMinus", 100, 0, 100, 100, 0, 100);
1516 new TH2F(
"h2_xtalMiscalibCoeffEndcapPlus",
"h2_xtalMiscalibCoeffEndcapPlus", 100, 0, 100, 100, 0, 100);
1524 for (
int i = 0;
i < 25;
i++) {
1542 sprintf(
histoName,
"h_ESCcorrEtrueVsEta_%d",
i);
1552 sprintf(
histoName,
"h2_xtalRecalibCoeffBarrel_%d",
i);
1558 sprintf(
histoName,
"h2_xtalRecalibCoeffEndcapMinus_%d",
i);
1563 sprintf(
histoName,
"h2_xtalRecalibCoeffEndcapPlus_%d",
i);
1588 h1_zMassResol_ =
new TH1F(
"zMassResol",
"zMassResol", 200, -50., 50.);
1592 h1_eleEtaResol_ =
new TH1F(
"eleEtaResol",
"eleEtaResol", 100, -0.01, 0.01);
1608 h1_elePhiResol_ =
new TH1F(
"elePhiResol",
"elePhiResol", 100, -0.01, 0.01);
1612 h1_zEtaResol_ =
new TH1F(
"zEtaResol",
"zEtaResol", 200, -1., 1.);
1616 h1_zPhiResol_ =
new TH1F(
"zPhiResol",
"zPhiResol", 200, -1., 1.);
1620 h1_nEleReco_ =
new TH1F(
"nEleReco",
"Number of reco electrons", 10, -0.5, 10.5);
1633 h1_occupancy_ =
new TH1F(
"occupancy",
"occupancy", 1000, 0, 10000);
1642 h1_eleClasses_ =
new TH1F(
"eleClasses",
"eleClasses", 301, -1, 300);
1654 h1_ZCandMult_ =
new TH1F(
"ZCandMult",
"Multiplicity of Z candidates in one event", 10, -0.5, 10.5);
1657 h1_reco_ZMass_ =
new TH1F(
"reco_ZMass",
"Inv. mass of 2 reco Electrons", 200, 0., 150.);
1661 h1_reco_ZMassCorr_ =
new TH1F(
"reco_ZMassCorr",
"Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1665 h1_reco_ZMassCorrBB_ =
new TH1F(
"reco_ZMassCorrBB",
"Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1669 h1_reco_ZMassCorrEE_ =
new TH1F(
"reco_ZMassCorrEE",
"Inv. mass of 2 corrected reco Electrons", 200, 0., 150.);
1675 h2_coeffVsEta_ =
new TH2F(
"h2_calibCoeffVsEta",
"h2_calibCoeffVsEta", 249, -124, 125, 200, 0., 2.);
1680 new TH2F(
"h2_calibCoeffVsEtaGrouped",
"h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1684 h2_zMassVsLoop_ =
new TH2F(
"h2_zMassVsLoop",
"h2_zMassVsLoop", 1000, 0, 40, 90, 80., 95.);
1686 h2_zMassDiffVsLoop_ =
new TH2F(
"h2_zMassDiffVsLoop",
"h2_zMassDiffVsLoop", 1000, 0, 40, 100, -1., 1.);
1690 h2_zWidthVsLoop_ =
new TH2F(
"h2_zWidthVsLoop",
"h2_zWidthVsLoop", 1000, 0, 40, 100, 0., 10.);
1692 h2_coeffVsLoop_ =
new TH2F(
"h2_coeffVsLoop",
"h2_coeffVsLoop", 1000, 0, 40, 100, 0., 2.);
1694 h2_residualSigma_ =
new TH2F(
"h2_residualSigma",
"h2_residualSigma", 1000, 0, 40, 100, 0., .5);
1696 h2_miscalRecal_ =
new TH2F(
"h2_miscalRecal",
"h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1700 h2_miscalRecalEB_ =
new TH2F(
"h2_miscalRecalEB",
"h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1704 h2_miscalRecalEE_ =
new TH2F(
"h2_miscalRecalEE",
"h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1708 h1_mc_ =
new TH1F(
"h1_residualMiscalib",
"h1_residualMiscalib", 200, -0.2, 0.2);
1709 h1_mcEB_ =
new TH1F(
"h1_residualMiscalibEB",
"h1_residualMiscalibEB", 200, -0.2, 0.2);
1710 h1_mcEE_ =
new TH1F(
"h1_residualMiscalibEE",
"h1_residualMiscalibEE", 200, -0.2, 0.2);
1712 for (
int i = 0;
i < 25;
i++) {
1731 sprintf(
histoName,
"h1_residualMiscalibParz_%d",
i);
1733 sprintf(
histoName,
"h1_residualMiscalibEBParz_%d",
i);
1735 sprintf(
histoName,
"h1_residualMiscalibEEParz_%d",
i);
1741 float p0 = 1.00153e+00;
1742 float p1 = 3.29331e-02;
1743 float p2 = 1.21187e-03;
1745 double x = (double)fabs(scEta);
1747 return 1. / (p0 +
p1 *
x *
x +
p2 *
x *
x *
x *
x);
1753 float p0 = 1.06819e+00;
1754 float p1 = -1.53189e-02;
1755 float p2 = 4.01707e-04;
1757 double x = (double)fabs(scEta);
1759 return 1. / (p0 +
p1 *
x *
x +
p2 *
x *
x *
x *
x);
1763 float p0 = 1.17382e+00;
1764 float p1 = -6.52319e-02;
1765 float p2 = 6.26108e-03;
1767 double x = (double)fabs(scEta);
1769 return 1. / (p0 +
p1 *
x *
x +
p2 *
x *
x *
x *
x);
1773 float p0 = 9.99782e-01;
1774 float p1 = 1.26983e-02;
1775 float p2 = 2.16344e-03;
1777 double x = (double)fabs(scEta);
1779 return 1. / (p0 +
p1 *
x *
x +
p2 *
x *
x *
x *
x);
1785 const std::vector<HepMC::GenParticle*>& mcEle,
1786 std::map<HepMC::GenParticle*, const reco::GsfElectron*>& myMCmap) {
1787 for (
unsigned int i = 0;
i < mcEle.size();
i++) {
1791 float dr =
EvalDR(mcEle[
i]->momentum().pseudoRapidity(),
1793 mcEle[
i]->momentum().
phi(),
1796 myMatchEle = (*electronCollection)[
j];
1800 myMCmap.insert(std::pair<HepMC::GenParticle*, const reco::GsfElectron*>(mcEle[
i], myMatchEle));
1813 float DEta = Eta - Eta_ref;
1824 return (
Phi - Phi_ref);
1828 std::map<HepMC::GenParticle*, const reco::GsfElectron*>& associationMap) {
1829 for (
unsigned int i = 0;
i < mcEle.size();
i++) {
1830 h_eleEffEta_[0]->Fill(fabs(mcEle[
i]->momentum().pseudoRapidity()));
1834 std::map<HepMC::GenParticle*, const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[
i]);
1835 if (mIter == associationMap.end())
1838 if ((*mIter).second) {
1841 h_eleEffEta_[1]->Fill(fabs(mcEle[
i]->momentum().pseudoRapidity()));
1871 if (
k >= 0 &&
k <= 84)
1874 if (
k >= 85 &&
k <= 169)
1877 if (
k >= 170 &&
k <= 208)
1880 if (
k >= 209 &&
k <= 247)
1886 if (
k >= 0 &&
k <= 71)
1889 if (
k >= 72 &&
k <= 143)
1899 if (
c == 0 ||
c == 10 ||
c == 20)
1902 if (
c == 100 ||
c == 110 ||
c == 120)
1905 if (
c == 30 ||
c == 31 ||
c == 32 ||
c == 33 ||
c == 34)
1908 if (
c == 130 ||
c == 131 ||
c == 132 ||
c == 133 ||
c == 134)
1920 std::pair<DetId, double> myPair(
DetId(0), -9999.);
1922 for (
std::vector<std::pair<DetId, float> >::const_iterator idIt = mySCRecHits.begin(); idIt != mySCRecHits.end();
1925 hottestRecHit = &(*(ebhits->
find((*idIt).first)));
1927 if (hottestRecHit == &(*(ebhits->
end()))) {
1928 std::cout <<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN" << std::endl;
1931 }
else if (idIt->first.subdetId() ==
EcalEndcap) {
1932 hottestRecHit = &(*(eehits->
find((*idIt).first)));
1933 if (hottestRecHit == &(*(eehits->
end()))) {
1934 std::cout <<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN" << std::endl;
1944 myPair.first = hottestRecHit->
id();
1966 for (
int i = 0;
i < 19;
i++) {
1977 for (
int i = 0;
i <
size;
i++) {
1980 bool isNearCrack =
false;
1992 dist +=
pow(v1[
i] - v2[
i], 2);
2041 for (
int i = 0;
i < 2;
i++) {
2079 std::cout <<
"[ CHECK ON BARREL ELECTRON NUMBER ]" 2082 std::cout <<
"[ EFFICIENCY OF THE BORDER SELECTION ]" 2085 std::cout <<
"[ EFFICIENCY OF THE GOLDEN SELECTION ] BARREL: " 2089 std::cout <<
"[ EFFICIENCY OF THE SILVER SELECTION ] BARREL: " 2093 std::cout <<
"[ EFFICIENCY OF THE SHOWER SELECTION ] BARREL: " 2097 std::cout <<
"[ EFFICIENCY OF THE CRACK SELECTION ] BARREL: " 2101 std::ofstream
fout(
"ZeeStatistics.txt");
2104 std::cout <<
"Cannot open output file.\n";
2107 fout <<
"ZeeStatistics" << std::endl;
2109 fout <<
"##########################RECO#########################" << std::endl;
2110 fout <<
"##################Zee with Barrel-Barrel electrons: " <<
BBZN << std::endl;
2114 fout <<
"##################Zee with Barrel-Endcap electrons: " <<
EBZN << std::endl;
2117 fout <<
"##################Zee with Endcap-Endcap electrons: " <<
EEZN << std::endl;
2121 fout <<
"\n" << std::endl;
2123 fout <<
"##########################GEN#########################" << std::endl;
2124 fout <<
"##################Zee with Barrel-Barrel electrons: " << (
float)
MCZBB /
NEVT << std::endl;
2125 fout <<
"##################Zee with Barrel-Endcap electrons: " << (
float)
MCZEB /
NEVT << std::endl;
2126 fout <<
"##################Zee with Endcap-Endcap electrons: " << (
float)
MCZEE /
NEVT << std::endl;
TH1 * weightedRescaleFactor[50][250]
const math::XYZPoint & position() const
cluster centroid position
bool accept() const
Has at least one path accepted the event?
auto setWhatProduced(T *iThis, const es::Label &iLabel={})
const edm::EDGetTokenT< edm::HepMCProduct > hepMCToken_
std::string outputFileName_
T getParameter(std::string const &) const
static const int MIN_IPHI
const std::vector< float > & getOptimizedCoefficientsError() const
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > geometryToken_
int BARREL_ELECTRONS_AFTER_BORDER_CUT
void fillEleInfo(std::vector< HepMC::GenParticle *> &a, std::map< HepMC::GenParticle *, const reco::GsfElectron *> &b)
TH1F * h1_eventsBeforeEWKSelection_
float computeCoefficientDistanceAtIteration(float v1[250], float v2[250], int size)
void fillEleInfo(const reco::GsfElectronCollection *)
void fillEleMCInfo(const HepMC::GenEvent *)
const edm::EDGetTokenT< EERecHitCollection > eeRecHitToken_
int GOLDEN_ELECTRONS_IN_ENDCAP
void fillZMCInfo(const HepMC::GenEvent *)
int iphi() const
get the crystal iphi
TH1F * h1_eventsAfterBorderSelection_
Sin< T >::type sin(const T &t)
T const * product() const
bool parseXMLMiscalibFile(std::string configFile)
void fillZInfo(std::pair< calib::CalibElectron *, calib::CalibElectron *> myZeeCandidate)
TH1F * h1_occupancyEndcap_
bool addEvent(calib::CalibElectron *, calib::CalibElectron *, float)
void writeMCEleHistograms()
Status duringLoop(const edm::Event &, const edm::EventSetup &) override
Called at each event.
TH1F * h1_eventsBeforeBorderSelection_
TH1F * h1_electronCosTheta_SC_
void writeMCZHistograms()
TH1F * h1_reco_ZMassCorr_
static void gausfit(TH1F *histoou, double *par, double *errpar, float nsigmalow, float nsigmaup, double *mychi2, int *iterations)
int GOLDEN_ELECTRONS_IN_BARREL
TH2F * h2_coeffVsEtaGrouped_
TH1F * h1_occupancyVsEta_
std::vector< GsfElectron > GsfElectronCollection
collection of GsfElectron objects
TH1F * h1_occupancyBarrel_
~ZeeCalibration() override
Destructor.
const edm::EDGetTokenT< edm::TriggerResults > trigResultsToken_
Classification classification() const
static bool validDetId(int i, int j)
check if a valid index combination
float EvalDR(float Eta, float Eta_ref, float Phi, float Phi_ref)
int ieta() const
get the crystal ieta
T getUntrackedParameter(std::string const &, T const &) const
U second(std::pair< T, U > const &p)
TH1F * h1_electronCosTheta_TK_
TH2F * h_ESCcorrEtrueVsEta_[25]
unsigned int size() const
Get number of paths stored.
const EcalIntercalibConstants & get()
int ringNumberCorrector(int k)
TupleMultiplicity< TrackerTraits > const *__restrict__ TrackingRecHitSoAConstView< TrackerTraits > TrackerTraits::tindex_type *__restrict__ double *__restrict__ phits
int SILVER_ELECTRONS_IN_BARREL
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
void beginOfJob() override
Called at beginning of job.
int TOTAL_ELECTRONS_IN_ENDCAP
int BARREL_ELECTRONS_BEFORE_BORDER_CUT
void fillHLTInfo(edm::Handle< edm::TriggerResults >)
int SHOWER_ELECTRONS_IN_ENDCAP
Abs< T >::type abs(const T &t)
ZIterativeAlgorithmWithFit * theAlgorithm_
const edm::EDGetTokenT< reco::SuperClusterCollection > islandSCToken_
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
unsigned int electronSelection_
bool xtalIsOnModuleBorder(EBDetId myEBDetId)
TH2F * h2_xtalRecalibCoeffBarrel_[25]
float EvalDPhi(float Phi, float Phi_ref)
TH1F * h1_eventsAfterEWKSelection_
T perp() const
Magnitude of transverse component.
TH1F * h1_preshowerOverSC_
const HepMC::GenEvent * GetEvent() const
const_iterator find(uint32_t rawId) const
TH1F * h1_electronCosTheta_SC_TK_
double coefficientDistanceAtIteration[50]
double fEtaBarrelBad(double scEta) const
void fillMCmap(const std::vector< const reco::GsfElectron *> *electronCollection, const std::vector< HepMC::GenParticle *> &mcEle, std::map< HepMC::GenParticle *, const reco::GsfElectron *> &myMCmap)
const_iterator end() const
TH2F * h2_zMassDiffVsLoop_
double energy() const
cluster energy
void startingNewLoop(unsigned int iLoop) override
Called at beginning of loop.
const edm::EDGetTokenT< reco::SuperClusterCollection > scToken_
int getNumberOfIterations() const
static const int MAX_IPHI
double fEtaEndcapGood(double scEta) const
void writeEleHistograms()
edm::ParameterSet theParameterSet
const ZIterativeAlgorithmWithFitPlots * getHistos() const
double fEtaBarrelGood(double scEta) const
constexpr uint32_t rawId() const
get the raw id
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
const std::vector< int > & getOptimizedIterations() const
TH1F * h_ESCcorrEtrue_[25]
std::pair< DetId, double > getHottestDetId(const std::vector< std::pair< DetId, float > > &mySCRecHits, const EBRecHitCollection *ebhits, const EERecHitCollection *eehits)
TH2F * h2_xtalMiscalibCoeffEndcapPlus_
static const int MAX_IETA
int getNumberOfChannels() const
const edm::EDGetTokenT< EBRecHitCollection > ebRecHitToken_
Status endOfLoop(const edm::EventSetup &, unsigned int iLoop) override
Called at end of loop.
const std::string mcProducer_
DetId id() const
get the id
float initCalibCoeff[250]
const std::vector< float > & getOptimizedCoefficients() const
int SILVER_ELECTRONS_IN_ENDCAP
TH1F * h1_reco_ZMassCorrBB_
const edm::EDGetTokenT< reco::GsfElectronCollection > gsfElectronToken_
const CaloClusterPtr & seed() const
seed BasicCluster
void writeLine(EBDetId const &, float)
iterator find(key_type k)
virtual std::shared_ptr< EcalIntercalibConstants > produceEcalIntercalibConstants(const EcalIntercalibConstantsRcd &iRecord)
Produce Ecal interCalibrations.
TH1F * h1_reco_ZMassCorrEE_
TH2F * h2_xtalMiscalibCoeffBarrel_
const std::vector< float > & getOptimizedChiSquare() const
TH1F * h1_weightSumMeanBarrel_
TH2F * h_ESCEtrueVsEta_[25]
float DPhi(float phi1, float phi2)
const self & getMap() const
TH2F * h2_iterations_[25]
ZeeCalibration(const edm::ParameterSet &iConfig)
Constructor.
float calibCoeffError[250]
TH2F * h2_xtalRecalibCoeffEndcapPlus_[25]
int CRACK_ELECTRONS_IN_ENDCAP
int CRACK_ELECTRONS_IN_BARREL
double fEtaEndcapBad(double scEta) const
TH2F * h2_fEtaEndcapGood_
TH2F * h2_xtalRecalibCoeffEndcapMinus_[25]
void bookEleMCHistograms()
double phi() const final
momentum azimuthal angle
double getEtaCorrection(const reco::GsfElectron *)
double preshowerEnergy() const
energy deposited in preshower
TH2F * h2_fEtaBarrelGood_
Power< A, B >::type pow(const A &a, const B &b)
TH1F * h1_borderElectronClassification_
SuperClusterRef superCluster() const override
reference to a SuperCluster
TH2F * h2_xtalMiscalibCoeffEndcapMinus_
void endOfJob() override
Called at end of job.
double sigmaErrorArray[50]
int SHOWER_ELECTRONS_IN_BARREL
int TOTAL_ELECTRONS_IN_BARREL
TH1F * h1_weightSumMeanEndcap_
double eta() const final
momentum pseudorapidity
std::shared_ptr< EcalIntercalibConstants > ical