18 #include "CLHEP/Vector/LorentzVector.h" 63 #include "TGraphErrors.h" 82 std::cout<<
"[ZeeCalibration] Starting the ctor"<<std::endl;
114 outputFile_ = TFile::Open(outputFileName_.c_str(),
"RECREATE");
116 myTree =
new TTree(
"myTree",
"myTree");
118 myTree->Branch(
"zMass",&
mass4tree,
"mass/F");
157 findingRecord<EcalIntercalibConstantsRcd> () ;
159 for(
int i = 0;
i<50;
i++){
169 std::cout<<
"[ZeeCalibration] Done with the ctor"<<std::endl;
183 std::shared_ptr<EcalIntercalibConstants>
186 std::cout <<
"@SUB=ZeeCalibration::produceEcalIntercalibConstants" << std::endl;
209 if(ieta==0)
continue;
243 std::cout<<
"Writing histos..."<<std::endl;
272 for(
unsigned int i =0;
i<25;
i++){
309 Double_t
mean[25] = {0.};
310 Double_t
num[25] = {0.};
311 Double_t meanErr[25] = {0.};
312 Float_t
rms[25] = {0.};
313 Float_t tempRms[10][25];
315 for(
int ia = 0; ia<10; ia++){
316 for(
int ib = 0;
ib<25;
ib++){
318 tempRms[ia][
ib] = 0.;
331 bool isNearCrack =
false;
396 if(
k>=170 &&
k<=204){
436 if(
k>=205 &&
k<=208){
484 for(
int ic = 0; ic< 17; ic++){
486 mean[ic] = mean[ic] / num[ic];
488 meanErr[ic] = 1. /
sqrt(meanErr[ic]);
494 for(
int ic = 0; ic< 25; ic++){
495 for(
int id = 0;
id< 10;
id++){
497 if(tempRms[
id][ic] > 0.){
499 rms[ic] += (tempRms[
id][ic] - mean[j])*(tempRms[
id][ic] - mean[j]);
504 rms[ic] =
sqrt(rms[ic]);
511 Double_t xtalEta[25] = {1.4425, 1.3567,1.2711,1.1855,
513 0.7468,0.6612,0.5756,0.4897,0.3985,0.3117,0.2250,0.1384,0.0487,
514 1.546, 1.651, 1.771, 1.908, 2.071, 2.267, 2.516, 2.8};
516 Double_t zero[25] = {0.026};
518 for(
int j = 0; j <25; j++)
525 px->SetXTitle(
"Eta channel");
526 px->SetYTitle(
"recalibCoeff");
546 double weightSumMeanBarrel = 0.;
547 double weightSumMeanEndcap = 0.;
585 std::cout<<
"Weight sum mean on channels in Barrel is :"<<weightSumMeanBarrel<<std::endl;
586 std::cout<<
"Weight sum mean on channels in Endcap is :"<<weightSumMeanEndcap<<std::endl;
598 TGraphErrors* graph =
new TGraphErrors(25,xtalEta,mean,zero,meanErr);
602 double zero50[50] = { 0. };
605 residualSigmaGraph->SetName(
"residualSigmaGraph");
606 residualSigmaGraph->Draw(
"APL");
607 residualSigmaGraph->Write();
610 coefficientDistanceAtIterationGraph->SetName(
"coefficientDistanceAtIterationGraph");
611 coefficientDistanceAtIterationGraph->Draw(
"APL");
612 coefficientDistanceAtIterationGraph->Write();
614 Float_t noError[250] = {0.};
616 Float_t ringInd[250];
617 for(
int i =0;
i<250;
i++)
621 graphCoeff->SetName(
"graphCoeff");
622 graphCoeff->Draw(
"APL");
661 std::cout<<
"[ZeeCalibration] Entering duringLoop"<<std::endl;
680 std::cout<<
"[ZeeCalibration::beginOfJob] Histograms booked "<<std::endl;
698 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EB miscal file"<<std::endl;
707 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EE miscal file"<<std::endl;
711 std::cout <<
" theAlgorithm_->getNumberOfChannels() " 721 std::vector<DetId> ringIds;
735 for (
unsigned int iid=0; iid<ringIds.size();++iid)
737 float miscalib=* (miscalibMap->
get().
getMap().
find(ringIds[iid]) );
750 ical = std::make_shared<EcalIntercalibConstants>();
756 std::vector<DetId> ringIds;
768 for (
unsigned int iid=0; iid<ringIds.size();++iid){
772 EBDetId myEBDetId(ringIds[iid]);
778 EEDetId myEEDetId(ringIds[iid]);
779 if(myEEDetId.
zside() < 0)
782 if(myEEDetId.
zside() > 0)
787 ical->setValue( ringIds[iid], *(miscalibMap->
get().
getMap().
find(ringIds[iid]) ) );
810 for(
unsigned int iHLT=0; iHLT<200; ++iHLT) {
815 std::cout<<
"[ZeeCalibration::duringLoop] Done with initializing aHLTresults[] "<<std::endl;
821 if(!hltTriggerResultHandle.
isValid()) {
831 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillHLTInfo(hltTriggerResultHandle); "<<std::endl;
849 std::cout<<
"[ZeeCalibration::duringLoop] End HLT section"<<std::endl;
855 std::vector<HepMC::GenParticle*> mcEle;
857 float myGenZMass(-1);
867 const HepMC::GenEvent * myGenEvent = hepProd->
GetEvent();
873 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillZMCInfo( & (*myGenEvent) ); "<<std::endl;
877 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
878 p != myGenEvent->particles_end(); ++
p ) {
880 if ( (*p)->pdg_id() == 23 && (*p)->status()==2){
882 myGenZMass = (*p)->momentum().m();
895 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
896 p != myGenEvent->particles_end(); ++
p ) {
898 if (
abs( (*p)->pdg_id() ) == 11 )
900 mcEle.push_back( (*
p) );
907 if(mcEle.size()==2 && fabs(mcEle[0]->momentum().
eta())<2.4 && fabs(mcEle[1]->momentum().
eta())<2.4 ){
910 if( fabs(mcEle[0]->momentum().
eta())<1.479 && fabs(mcEle[1]->momentum().
eta())<1.479 )
MCZBB++;
912 if( (fabs(mcEle[0]->momentum().
eta())>1.479 && fabs(mcEle[1]->momentum().
eta())<1.479) || (fabs(mcEle[0]->momentum().
eta())<1.479 && fabs(mcEle[1]->momentum().
eta())>1.479) )
MCZEB++;
914 if( fabs(mcEle[0]->momentum().eta())>1.479 && fabs(mcEle[1]->momentum().eta())>1.479 )
MCZEE++;
928 std::cerr <<
"Error! can't get the product EBRecHitCollection " << std::endl;
937 std::cerr <<
"Error! can't get the product EERecHitCollection " << std::endl;
947 std::cerr <<
"Error! can't get the product SuperClusterCollection "<< std::endl;
952 std::cout<<
"scCollection->size()"<<scCollection->size()<<std::endl;
953 for(reco::SuperClusterCollection::const_iterator scIt = scCollection->begin(); scIt != scCollection->end(); scIt++)
964 std::cerr <<
"Error! can't get the product IslandSuperClusterCollection "<< std::endl;
969 std::cout<<
"scCollection->size()"<<scIslandCollection->size()<<std::endl;
972 if( ( scCollection->size()+scIslandCollection->size() ) < 2)
980 std::cerr <<
"Error! can't get the product ElectronCollection "<< std::endl;
993 if(electronCollection->size() < 2)
996 if ( !hits && !ehits){
1001 if (hits->
size() == 0 && ehits->
size() == 0){
1002 std::cout <<
"hits->size() == 0" << std::endl;
1006 if (!electronCollection){
1007 std::cout <<
"!electronCollection" << std::endl;
1011 if (electronCollection->size() == 0){
1012 std::cout <<
"electronCollection->size() == 0" << std::endl;
1029 std::cout <<
" Starting with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1036 std::cout <<
" Done with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1042 std::vector<calib::CalibElectron> calibElectrons;
1048 for(
unsigned int e_it = 0 ; e_it != electronCollection->size() ; e_it++)
1052 std::cout << calibElectrons.back().getRecoElectron()->superCluster()->energy() <<
" " << calibElectrons.back().getRecoElectron()->energy() << std::endl;
1060 std::cout <<
"Filled histos" << std::endl;
1064 std::vector<std::pair<calib::CalibElectron*,calib::CalibElectron*> > zeeCandidates;
1068 double DeltaMinvMin(5000.);
1070 if (calibElectrons.size() < 2)
1073 for(
unsigned int e_it = 0 ; e_it != calibElectrons.size() - 1 ; e_it++){
1074 for(
unsigned int p_it = e_it + 1 ; p_it != calibElectrons.size() ; p_it++)
1077 std::cout << e_it <<
" " << calibElectrons[e_it].getRecoElectron()->charge() <<
" " << p_it <<
" " << calibElectrons[p_it].getRecoElectron()->charge() << std::endl;
1079 if (calibElectrons[e_it].getRecoElectron()->
charge() * calibElectrons[p_it].getRecoElectron()->
charge() != -1)
1088 std::cout <<
"#######################mass "<<
mass << std::endl;
1091 zeeCandidates.push_back(std::pair<calib::CalibElectron*,calib::CalibElectron*>(&(calibElectrons[e_it]),&(calibElectrons[p_it])));
1092 double DeltaMinv = fabs(
mass -
MZ);
1094 if( DeltaMinv < DeltaMinvMin)
1096 DeltaMinvMin = DeltaMinv;
1097 myBestZ=zeeCandidates.size()-1;
1108 if(zeeCandidates.size()==0 || myBestZ==-1 )
1115 std::cout <<
"Found ZCandidates " << myBestZ << std::endl;
1123 int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1124 int class2 = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1155 if( class1==10 || class1 ==20)
1157 if( class1==110 || class1 ==120)
1159 if( class1>=30 && class1 <=34)
1161 if( class1>=130 && class1 <=134)
1172 if( class2==10 || class2 ==20)
1174 if( class2==110 || class2 ==120)
1176 if( class2>=30 && class2 <=34)
1178 if( class2>=130 && class2 <=134)
1193 bool firstElectronIsOnModuleBorder(
false);
1194 bool secondElectronIsOnModuleBorder(
false);
1205 if( firstElehottestDetId.
subdetId() ==
EcalBarrel && !firstElectronIsOnModuleBorder )
1212 if( secondElehottestDetId.subdetId() ==
EcalBarrel)
1217 if( secondElehottestDetId.subdetId() ==
EcalBarrel && !secondElectronIsOnModuleBorder )
1224 if ( firstElehottestDetId.
subdetId() ==
EcalBarrel && firstElectronIsOnModuleBorder ){
1231 if ( secondElehottestDetId.subdetId() ==
EcalBarrel && secondElectronIsOnModuleBorder ){
1241 if(class1<100 && class2<100){
1243 if(class1==0 && class2==0)
BBZN_gg++;
1244 if(class1<21 && class2<21)
BBZN_tt++;
1245 if(class1<21 || class2<21)
BBZN_t0++;
1249 if(class1>=100 && class2>=100){
1251 if(class1==100 && class2==100)
EEZN_gg++;
1252 if(class1<121 && class2<121)
EEZN_tt++;
1253 if(class1<121 || class2<121)
EEZN_t0++;
1257 if( (class1<100 && class2>=100) || (class2<100 && class1>=100)){
1259 if( (class1==0 && class2==100)||(class2==0 && class1==100) )
EBZN_gg++;
1260 if( ( class1<21 && class2<121) ||(class2<21 && class1<121) )
EBZN_tt++;
1261 if( class2<21 || class1<21 || class2<121 || class1<121 )
EBZN_t0++;
1273 bool selectionBool=
false;
1278 float theta1 = 2. * atan(
exp(- zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->eta()) );
1279 bool ET_1 = ( (zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->energy() *
sin( theta1) ) > 20.);
1281 float theta2 = 2. * atan(
exp(- zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->eta()) );
1282 bool ET_2 = ( (zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->energy() *
sin( theta2) ) > 20.);
1285 bool HoE_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->hadronicOverEm() < 0.115);
1286 bool HoE_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->hadronicOverEm() < 0.115);
1288 bool DeltaPhiIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1289 bool DeltaPhiIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1291 bool DeltaEtaIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1292 bool DeltaEtaIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1296 if(! (invMassBool &&
1299 DeltaPhiIn_1 && DeltaPhiIn_2 &&
1300 DeltaEtaIn_1 && DeltaEtaIn_2
1310 int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1311 int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1417 float ele1EnergyCorrection(1.);
1418 float ele2EnergyCorrection(1.);
1427 if (invMassBool && selectionBool)
1440 std::map<HepMC::GenParticle*,const reco::GsfElectron*> myMCmap;
1442 std::vector<const reco::GsfElectron*> dauElectronCollection;
1444 dauElectronCollection.push_back(zeeCandidates[myBestZ].
first->getRecoElectron() );
1445 dauElectronCollection.push_back(zeeCandidates[myBestZ].
second->getRecoElectron() );
1447 fillMCmap(&dauElectronCollection,mcEle,myMCmap);
1459 int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1460 int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1461 if(c1st<100 && c2nd<100 )
1465 if(c1st>=100 && c2nd>=100 )
1482 std::cout <<
"Added event to algorithm" << std::endl;
1494 std::cout<<
"[ZeeCalibration] Starting loop number " << iLoop<<std::endl;
1503 std::cout<<
"[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1533 std::cout<<
"[ZeeCalibration] Ending loop " << iLoop<<std::endl;
1545 std::cout<<
"Optimized coefficients " << optimizedCoefficients.size() <<std::endl;
1551 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1565 std::cout<<
"size "<<optimizedCoefficients.size()<<std::endl;
1567 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1569 calibCoeff[ieta] *= optimizedCoefficients[ieta];
1575 std::cout<< ieta <<
" " << optimizedCoefficients[ieta] <<std::endl;
1578 std::vector<DetId> ringIds;
1590 for (
unsigned int iid=0; iid<ringIds.size();++iid){
1593 EBDetId myEBDetId(ringIds[iid]);
1599 EEDetId myEEDetId(ringIds[iid]);
1600 if(myEEDetId.
zside() < 0)
1603 if(myEEDetId.
zside() > 0)
1609 ical->setValue( ringIds[iid], *(
ical->getMap().find(ringIds[iid]) ) * optimizedCoefficients[ieta] );
1646 double parResidual[3];
1647 double errparResidual[3];
1659 std::cout<<
"Fit on residuals, sigma is "<<parResidual[2]<<
" +/- "<<errparResidual[2]<<std::endl;
1698 h1_seedOverSC_=
new TH1F(
"h1_seedOverSC",
"h1_seedOverSC", 400, 0., 2.);
1705 h2_fEtaBarrelGood_ =
new TH2F(
"fEtaBarrelGood",
"fEtaBarrelGood",800,-4.,4.,800,0.8,1.2);
1709 h2_fEtaBarrelBad_ =
new TH2F(
"fEtaBarrelBad",
"fEtaBarrelBad",800,-4.,4.,800,0.8,1.2);
1713 h2_fEtaEndcapGood_ =
new TH2F(
"fEtaEndcapGood",
"fEtaEndcapGood",800,-4.,4.,800,0.8,1.2);
1717 h2_fEtaEndcapBad_ =
new TH2F(
"fEtaEndcapBad",
"fEtaEndcapBad",800,-4.,4.,800,0.8,1.2);
1721 for (
int i=0;
i<2;
i++)
1725 sprintf(histoName,
"h_eleEffEta_%d",
i);
1726 h_eleEffEta_[
i] =
new TH1F(histoName,histoName, 150, 0., 2.7);
1729 sprintf(histoName,
"h_eleEffPhi_%d",
i);
1730 h_eleEffPhi_[
i] =
new TH1F(histoName,histoName, 400, -4., 4.);
1733 sprintf(histoName,
"h_eleEffPt_%d",
i);
1734 h_eleEffPt_[
i] =
new TH1F(histoName,histoName, 200, 0., 200.);
1748 for (
int i=0;
i<25;
i++)
1752 sprintf(histoName,
"h_ESCEtrueVsEta_%d",
i);
1758 sprintf(histoName,
"h_ESCEtrue_%d",
i);
1760 h_ESCEtrue_[
i] =
new TH1F(histoName,histoName, 300,0.,1.5);
1762 sprintf(histoName,
"h2_chi2_%d",
i);
1763 h2_chi2_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 5);
1765 sprintf(histoName,
"h2_iterations_%d",
i);
1766 h2_iterations_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 15);
1768 sprintf(histoName,
"h_ESCcorrEtrueVsEta_%d",
i);
1774 sprintf(histoName,
"h_ESCcorrEtrue_%d",
i);
1778 sprintf(histoName,
"h2_xtalRecalibCoeffBarrel_%d",
i);
1784 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapMinus_%d",
i);
1789 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapPlus_%d",
i);
1815 h1_zMassResol_ =
new TH1F(
"zMassResol",
"zMassResol", 200, -50., 50.);
1819 h1_eleEtaResol_ =
new TH1F(
"eleEtaResol",
"eleEtaResol", 100, -0.01, 0.01);
1835 h1_elePhiResol_ =
new TH1F(
"elePhiResol",
"elePhiResol", 100, -0.01, 0.01);
1840 h1_zEtaResol_ =
new TH1F(
"zEtaResol",
"zEtaResol", 200, -1., 1.);
1845 h1_zPhiResol_ =
new TH1F(
"zPhiResol",
"zPhiResol", 200, -1., 1.);
1849 h1_nEleReco_ =
new TH1F(
"nEleReco",
"Number of reco electrons",10,-0.5,10.5);
1862 h1_occupancy_ =
new TH1F(
"occupancy",
"occupancy",1000,0,10000);
1887 h1_ZCandMult_ =
new TH1F(
"ZCandMult",
"Multiplicity of Z candidates in one event",10,-0.5,10.5);
1890 h1_reco_ZMass_ =
new TH1F(
"reco_ZMass",
"Inv. mass of 2 reco Electrons",200,0.,150.);
1894 h1_reco_ZMassCorr_ =
new TH1F(
"reco_ZMassCorr",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1898 h1_reco_ZMassCorrBB_ =
new TH1F(
"reco_ZMassCorrBB",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1902 h1_reco_ZMassCorrEE_ =
new TH1F(
"reco_ZMassCorrEE",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1908 h2_coeffVsEta_=
new TH2F(
"h2_calibCoeffVsEta",
"h2_calibCoeffVsEta",249,-124,125, 200, 0., 2.);
1912 h2_coeffVsEtaGrouped_=
new TH2F(
"h2_calibCoeffVsEtaGrouped",
"h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1916 h2_zMassVsLoop_=
new TH2F(
"h2_zMassVsLoop",
"h2_zMassVsLoop",1000,0,40, 90, 80.,95.);
1918 h2_zMassDiffVsLoop_=
new TH2F(
"h2_zMassDiffVsLoop",
"h2_zMassDiffVsLoop",1000,0,40, 100, -1., 1.);
1922 h2_zWidthVsLoop_=
new TH2F(
"h2_zWidthVsLoop",
"h2_zWidthVsLoop",1000,0,40, 100, 0.,10.);
1924 h2_coeffVsLoop_=
new TH2F(
"h2_coeffVsLoop",
"h2_coeffVsLoop",1000,0,40, 100, 0., 2.);
1926 h2_residualSigma_=
new TH2F(
"h2_residualSigma",
"h2_residualSigma",1000, 0, 40, 100, 0., .5);
1928 h2_miscalRecal_ =
new TH2F(
"h2_miscalRecal",
"h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1932 h2_miscalRecalEB_ =
new TH2F(
"h2_miscalRecalEB",
"h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1936 h2_miscalRecalEE_ =
new TH2F(
"h2_miscalRecalEE",
"h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1940 h1_mc_ =
new TH1F(
"h1_residualMiscalib",
"h1_residualMiscalib", 200, -0.2, 0.2);
1941 h1_mcEB_ =
new TH1F(
"h1_residualMiscalibEB",
"h1_residualMiscalibEB", 200, -0.2, 0.2);
1942 h1_mcEE_ =
new TH1F(
"h1_residualMiscalibEE",
"h1_residualMiscalibEE", 200, -0.2, 0.2);
1944 for (
int i=0;
i<25;
i++)
1964 sprintf(histoName,
"h1_residualMiscalibParz_%d",
i);
1965 h1_mcParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1966 sprintf(histoName,
"h1_residualMiscalibEBParz_%d",
i);
1967 h1_mcEBParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1968 sprintf(histoName,
"h1_residualMiscalibEEParz_%d",
i);
1969 h1_mcEEParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1979 float p0 = 1.00153e+00;
1980 float p1 = 3.29331e-02;
1981 float p2 = 1.21187e-03;
1983 double x = (double) fabs(scEta);
1985 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
1993 float p0 = 1.06819e+00;
1994 float p1 = -1.53189e-02;
1995 float p2 = 4.01707e-04 ;
1997 double x = (double) fabs(scEta);
1999 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2005 float p0 = 1.17382e+00;
2006 float p1 = -6.52319e-02;
2007 float p2 = 6.26108e-03;
2009 double x = (double) fabs(scEta);
2011 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2017 float p0 = 9.99782e-01 ;
2018 float p1 = 1.26983e-02;
2019 float p2 = 2.16344e-03;
2021 double x = (double) fabs(scEta);
2023 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2032 for (
unsigned int i=0;
i<mcEle.size();
i++)
2036 for (
unsigned int j=0;j<electronCollection->size();j++)
2038 float dr=
EvalDR(mcEle[
i]->momentum().pseudoRapidity(),(*(*electronCollection)[j]).
eta(),mcEle[
i]->momentum().
phi(),(*(*electronCollection)[j]).
phi());
2041 myMatchEle = (*electronCollection)[j];
2045 myMCmap.insert(std::pair<HepMC::GenParticle*,const reco::GsfElectron*>(mcEle[
i],myMatchEle));
2053 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2054 float DPhi = Phi - Phi_ref ;
2057 float DEta = Eta - Eta_ref ;
2059 float DR =
sqrt( DEta*DEta + DPhi*DPhi );
2066 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2067 return (Phi - Phi_ref);
2073 for (
unsigned int i=0;
i<mcEle.size();
i++)
2076 h_eleEffEta_[0]->Fill(fabs(mcEle[
i]->momentum().pseudoRapidity()));
2080 std::map<HepMC::GenParticle*,const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[
i]);
2081 if (mIter == associationMap.end() )
2088 h_eleEffEta_[1]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
2123 if(k>=0 && k<=84)index = k - 85;
2125 if(k>=85 && k<=169)index = k - 84;
2127 if(k>=170 && k<=208)index = - k + 84;
2129 if(k>=209 && k<=247)index = k - 123;
2135 if(k>=0 && k<=71)index = k - 72;
2137 if(k>=72 && k<=143)index = k - 71;
2147 double correction(1.);
2184 std::pair<DetId, double> myPair (
DetId(0), -9999.);
2187 for( std::vector<std::pair<DetId,float> >::const_iterator idIt=mySCRecHits.begin(); idIt != mySCRecHits.end(); idIt++){
2191 hottestRecHit = & (* ( ebhits->
find((*idIt).first) ) );
2193 if( hottestRecHit == & (*( ebhits->
end())) )
2195 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2199 else if (idIt->first.subdetId() ==
EcalEndcap )
2201 hottestRecHit = & (* ( eehits->
find((*idIt).first) ) );
2202 if( hottestRecHit == & (*( eehits->
end())) )
2204 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2213 maxEnergy = hottestRecHit->
energy();
2215 myPair.first = hottestRecHit ->
id();
2233 short ieta = myEBDetId.
ieta();
2234 short iphi = myEBDetId.
iphi();
2238 myBool = (
abs( ieta ) == 1 ||
abs( ieta ) == 25
2239 ||
abs( ieta ) ==26 ||
abs( ieta ) == 45
2240 ||
abs( ieta ) ==46 ||
abs( ieta ) == 65
2241 ||
abs( ieta ) ==66 ||
abs( ieta ) == 85 );
2243 for(
int i = 0;
i < 19;
i++){
2245 if(iphi == ( 20*
i + 1 ) || iphi == 20*
i )
2261 bool isNearCrack =
false;
2273 dist +=
pow( v1[
i]-v2[
i], 2 );
2329 for (
int i=0;
i<2;
i++)
2387 std::ofstream
fout(
"ZeeStatistics.txt");
2390 std::cout <<
"Cannot open output file.\n";
2393 fout<<
"ZeeStatistics"<<std::endl;
2395 fout<<
"##########################RECO#########################"<<std::endl;
2396 fout<<
"##################Zee with Barrel-Barrel electrons: "<<
BBZN<<std::endl;
2399 fout<<
"##################Zee with Barrel-Endcap electrons: "<<
EBZN<<std::endl;
2401 fout<<
"##################Zee with Endcap-Endcap electrons: "<<
EEZN<<std::endl;
2404 fout<<
"\n"<<std::endl;
2406 fout<<
"##########################GEN#########################"<<std::endl;
2407 fout<<
"##################Zee with Barrel-Barrel electrons: "<<(
float)
MCZBB/
NEVT<<std::endl;
2408 fout<<
"##################Zee with Barrel-Endcap electrons: "<<(
float)
MCZEB/
NEVT<<std::endl;
2409 fout<<
"##################Zee with Endcap-Endcap electrons: "<<(
float)
MCZEE/
NEVT<<std::endl;
TH1 * weightedRescaleFactor[50][250]
virtual void endOfJob()
Called at end of job.
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
const math::XYZPoint & position() const
cluster centroid position
const ZIterativeAlgorithmWithFitPlots * getHistos() const
std::string outputFileName_
static const int MIN_IPHI
double fEtaEndcapGood(double scEta) const
std::string scIslandCollection_
std::string electronCollection_
int BARREL_ELECTRONS_AFTER_BORDER_CUT
TH1F * h1_eventsBeforeEWKSelection_
double fEtaBarrelBad(double scEta) const
float computeCoefficientDistanceAtIteration(float v1[250], float v2[250], int size)
virtual Status endOfLoop(const edm::EventSetup &, unsigned int iLoop)
Called at end of loop.
virtual void beginOfJob()
Called at beginning of job.
virtual double eta() const final
momentum pseudorapidity
const self & getMap() const
void fillEleInfo(const reco::GsfElectronCollection *)
void fillEleMCInfo(const HepMC::GenEvent *)
const std::vector< float > & getOptimizedChiSquare() const
int GOLDEN_ELECTRONS_IN_ENDCAP
void fillZMCInfo(const HepMC::GenEvent *)
TH1F * h1_eventsAfterBorderSelection_
int getNumberOfChannels() const
edm::InputTag hlTriggerResults_
Sin< T >::type sin(const T &t)
bool accept() const
Has at least one path accepted the event?
bool parseXMLMiscalibFile(std::string configFile)
void fillMCmap(const std::vector< const reco::GsfElectron * > *electronCollection, const std::vector< HepMC::GenParticle * > &mcEle, std::map< HepMC::GenParticle *, const reco::GsfElectron * > &myMCmap)
const std::vector< float > & getOptimizedCoefficients() const
TH1F * h1_occupancyEndcap_
bool addEvent(calib::CalibElectron *, calib::CalibElectron *, float)
void writeMCEleHistograms()
virtual Status duringLoop(const edm::Event &, const edm::EventSetup &)
Called at each event.
TH1F * h1_eventsBeforeBorderSelection_
TH1F * h1_electronCosTheta_SC_
void writeMCZHistograms()
double fEtaEndcapBad(double scEta) const
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
std::string rechitCollection_
void setWhatProduced(T *iThis, const es::Label &iLabel=es::Label())
TH1F * h1_occupancyBarrel_
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 iphi() const
get the crystal iphi
uint32_t rawId() const
get the raw id
virtual double phi() const final
momentum azimuthal angle
U second(std::pair< T, U > const &p)
TH1F * h1_electronCosTheta_TK_
~ZeeCalibration()
Destructor.
TH2F * h_ESCcorrEtrueVsEta_[25]
std::string rechitProducer_
const EcalIntercalibConstants & get()
int ringNumberCorrector(int k)
std::string scIslandProducer_
int SILVER_ELECTRONS_IN_BARREL
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
int TOTAL_ELECTRONS_IN_ENDCAP
int BARREL_ELECTRONS_BEFORE_BORDER_CUT
unsigned int size() const
Get number of paths stored.
void fillHLTInfo(edm::Handle< edm::TriggerResults >)
int SHOWER_ELECTRONS_IN_ENDCAP
virtual SuperClusterRef superCluster() const
reference to a SuperCluster
int getNumberOfIterations() const
std::string erechitProducer_
Abs< T >::type abs(const T &t)
double energy() const
cluster energy
ZIterativeAlgorithmWithFit * theAlgorithm_
int ieta() const
get the crystal ieta
unsigned int electronSelection_
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
bool xtalIsOnModuleBorder(EBDetId myEBDetId)
TH2F * h2_xtalRecalibCoeffBarrel_[25]
float EvalDPhi(float Phi, float Phi_ref)
TH1F * h1_eventsAfterEWKSelection_
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
TH1F * h1_preshowerOverSC_
TH1F * h1_electronCosTheta_SC_TK_
double coefficientDistanceAtIteration[50]
const_iterator end() const
TH2F * h2_zMassDiffVsLoop_
std::string scCollection_
const std::vector< float > & getOptimizedCoefficientsError() const
const HepMC::GenEvent * GetEvent() const
static const int MAX_IPHI
DetId id() const
get the id
T const * product() const
void writeEleHistograms()
edm::ParameterSet theParameterSet
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Classification classification() 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
float initCalibCoeff[250]
int SILVER_ELECTRONS_IN_ENDCAP
TH1F * h1_reco_ZMassCorrBB_
void writeLine(EBDetId const &, float)
iterator find(key_type k)
T perp() const
Magnitude of transverse component.
virtual std::shared_ptr< EcalIntercalibConstants > produceEcalIntercalibConstants(const EcalIntercalibConstantsRcd &iRecord)
Produce Ecal interCalibrations.
virtual void startingNewLoop(unsigned int iLoop)
Called at beginning of loop.
TH1F * h1_reco_ZMassCorrEE_
TH2F * h2_xtalMiscalibCoeffBarrel_
void fillZInfo(std::pair< calib::CalibElectron *, calib::CalibElectron * > myZeeCandidate)
void fillEleInfo(std::vector< HepMC::GenParticle * > &a, std::map< HepMC::GenParticle *, const reco::GsfElectron * > &b)
const_iterator find(uint32_t rawId) const
TH1F * h1_weightSumMeanBarrel_
double fEtaBarrelGood(double scEta) const
TH2F * h_ESCEtrueVsEta_[25]
TH2F * h2_iterations_[25]
const CaloClusterPtr & seed() const
seed BasicCluster
ZeeCalibration(const edm::ParameterSet &iConfig)
Constructor.
void Reset(std::vector< TH2F > &depth)
double preshowerEnergy() const
energy deposited in preshower
float calibCoeffError[250]
TH2F * h2_xtalRecalibCoeffEndcapPlus_[25]
int CRACK_ELECTRONS_IN_ENDCAP
int CRACK_ELECTRONS_IN_BARREL
TH2F * h2_fEtaEndcapGood_
TH2F * h2_xtalRecalibCoeffEndcapMinus_[25]
void bookEleMCHistograms()
std::string electronProducer_
double getEtaCorrection(const reco::GsfElectron *)
TH2F * h2_fEtaBarrelGood_
std::string erechitCollection_
Power< A, B >::type pow(const A &a, const B &b)
TH1F * h1_borderElectronClassification_
TH2F * h2_xtalMiscalibCoeffEndcapMinus_
double sigmaErrorArray[50]
int SHOWER_ELECTRONS_IN_BARREL
const std::vector< int > & getOptimizedIterations() const
int TOTAL_ELECTRONS_IN_BARREL
TH1F * h1_weightSumMeanEndcap_
std::shared_ptr< EcalIntercalibConstants > ical