14 #include <TGraphErrors.h> 21 #include <CLHEP/Vector/LorentzVector.h> 71 std::cout<<
"[ZeeCalibration] Starting the ctor"<<std::endl;
103 outputFile_ = TFile::Open(outputFileName_.c_str(),
"RECREATE");
105 myTree =
new TTree(
"myTree",
"myTree");
107 myTree->Branch(
"zMass",&
mass4tree,
"mass/F");
146 findingRecord<EcalIntercalibConstantsRcd> () ;
148 for(
int i = 0;
i<50;
i++){
158 std::cout<<
"[ZeeCalibration] Done with the ctor"<<std::endl;
172 std::shared_ptr<EcalIntercalibConstants>
175 std::cout <<
"@SUB=ZeeCalibration::produceEcalIntercalibConstants" << std::endl;
198 if(ieta==0)
continue;
232 std::cout<<
"Writing histos..."<<std::endl;
261 for(
unsigned int i =0;
i<25;
i++){
298 Double_t
mean[25] = {0.};
299 Double_t
num[25] = {0.};
300 Double_t meanErr[25] = {0.};
301 Float_t
rms[25] = {0.};
302 Float_t tempRms[10][25];
304 for(
int ia = 0; ia<10; ia++){
305 for(
int ib = 0;
ib<25;
ib++){
307 tempRms[ia][
ib] = 0.;
320 bool isNearCrack =
false;
385 if(
k>=170 &&
k<=204){
425 if(
k>=205 &&
k<=208){
473 for(
int ic = 0; ic< 17; ic++){
475 mean[ic] = mean[ic] / num[ic];
477 meanErr[ic] = 1. /
sqrt(meanErr[ic]);
483 for(
int ic = 0; ic< 25; ic++){
484 for(
int id = 0;
id< 10;
id++){
486 if(tempRms[
id][ic] > 0.){
488 rms[ic] += (tempRms[
id][ic] - mean[j])*(tempRms[
id][ic] - mean[j]);
493 rms[ic] =
sqrt(rms[ic]);
500 Double_t xtalEta[25] = {1.4425, 1.3567,1.2711,1.1855,
502 0.7468,0.6612,0.5756,0.4897,0.3985,0.3117,0.2250,0.1384,0.0487,
503 1.546, 1.651, 1.771, 1.908, 2.071, 2.267, 2.516, 2.8};
505 Double_t zero[25] = {0.026};
507 for(
int j = 0; j <25; j++)
514 px->SetXTitle(
"Eta channel");
515 px->SetYTitle(
"recalibCoeff");
535 double weightSumMeanBarrel = 0.;
536 double weightSumMeanEndcap = 0.;
574 std::cout<<
"Weight sum mean on channels in Barrel is :"<<weightSumMeanBarrel<<std::endl;
575 std::cout<<
"Weight sum mean on channels in Endcap is :"<<weightSumMeanEndcap<<std::endl;
587 TGraphErrors* graph =
new TGraphErrors(25,xtalEta,mean,zero,meanErr);
591 double zero50[50] = { 0. };
594 residualSigmaGraph->SetName(
"residualSigmaGraph");
595 residualSigmaGraph->Draw(
"APL");
596 residualSigmaGraph->Write();
599 coefficientDistanceAtIterationGraph->SetName(
"coefficientDistanceAtIterationGraph");
600 coefficientDistanceAtIterationGraph->Draw(
"APL");
601 coefficientDistanceAtIterationGraph->Write();
603 Float_t noError[250] = {0.};
605 Float_t ringInd[250];
606 for(
int i =0;
i<250;
i++)
610 graphCoeff->SetName(
"graphCoeff");
611 graphCoeff->Draw(
"APL");
650 std::cout<<
"[ZeeCalibration] Entering duringLoop"<<std::endl;
669 std::cout<<
"[ZeeCalibration::beginOfJob] Histograms booked "<<std::endl;
687 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EB miscal file"<<std::endl;
696 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EE miscal file"<<std::endl;
700 std::cout <<
" theAlgorithm_->getNumberOfChannels() " 710 std::vector<DetId> ringIds;
724 for (
unsigned int iid=0; iid<ringIds.size();++iid)
726 float miscalib=* (miscalibMap->
get().
getMap().
find(ringIds[iid]) );
739 ical = std::make_shared<EcalIntercalibConstants>();
745 std::vector<DetId> ringIds;
757 for (
unsigned int iid=0; iid<ringIds.size();++iid){
761 EBDetId myEBDetId(ringIds[iid]);
767 EEDetId myEEDetId(ringIds[iid]);
768 if(myEEDetId.
zside() < 0)
771 if(myEEDetId.
zside() > 0)
776 ical->setValue( ringIds[iid], *(miscalibMap->
get().
getMap().
find(ringIds[iid]) ) );
799 for(
unsigned int iHLT=0; iHLT<200; ++iHLT) {
804 std::cout<<
"[ZeeCalibration::duringLoop] Done with initializing aHLTresults[] "<<std::endl;
810 if(!hltTriggerResultHandle.
isValid()) {
820 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillHLTInfo(hltTriggerResultHandle); "<<std::endl;
838 std::cout<<
"[ZeeCalibration::duringLoop] End HLT section"<<std::endl;
844 std::vector<HepMC::GenParticle*> mcEle;
846 float myGenZMass(-1);
862 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillZMCInfo( & (*myGenEvent) ); "<<std::endl;
866 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
867 p != myGenEvent->particles_end(); ++
p ) {
869 if ( (*p)->pdg_id() == 23 && (*p)->status()==2){
871 myGenZMass = (*p)->momentum().m();
884 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
885 p != myGenEvent->particles_end(); ++
p ) {
887 if (
abs( (*p)->pdg_id() ) == 11 )
889 mcEle.push_back( (*
p) );
896 if(mcEle.size()==2 && fabs(mcEle[0]->momentum().
eta())<2.4 && fabs(mcEle[1]->momentum().
eta())<2.4 ){
899 if( fabs(mcEle[0]->momentum().
eta())<1.479 && fabs(mcEle[1]->momentum().
eta())<1.479 )
MCZBB++;
901 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++;
903 if( fabs(mcEle[0]->momentum().eta())>1.479 && fabs(mcEle[1]->momentum().eta())>1.479 )
MCZEE++;
917 std::cerr <<
"Error! can't get the product EBRecHitCollection " << std::endl;
926 std::cerr <<
"Error! can't get the product EERecHitCollection " << std::endl;
936 std::cerr <<
"Error! can't get the product SuperClusterCollection "<< std::endl;
941 std::cout<<
"scCollection->size()"<<scCollection->size()<<std::endl;
942 for(reco::SuperClusterCollection::const_iterator scIt = scCollection->begin(); scIt != scCollection->end(); scIt++)
953 std::cerr <<
"Error! can't get the product IslandSuperClusterCollection "<< std::endl;
958 std::cout<<
"scCollection->size()"<<scIslandCollection->size()<<std::endl;
961 if( ( scCollection->size()+scIslandCollection->size() ) < 2)
969 std::cerr <<
"Error! can't get the product ElectronCollection "<< std::endl;
982 if(electronCollection->size() < 2)
985 if ( !hits && !ehits){
991 std::cout <<
"hits->size() == 0" << std::endl;
995 if (!electronCollection){
996 std::cout <<
"!electronCollection" << std::endl;
1000 if (electronCollection->empty()){
1001 std::cout <<
"electronCollection->size() == 0" << std::endl;
1018 std::cout <<
" Starting with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1025 std::cout <<
" Done with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1031 std::vector<calib::CalibElectron> calibElectrons;
1037 for(
unsigned int e_it = 0 ; e_it != electronCollection->size() ; e_it++)
1041 std::cout << calibElectrons.back().getRecoElectron()->superCluster()->energy() <<
" " << calibElectrons.back().getRecoElectron()->energy() << std::endl;
1049 std::cout <<
"Filled histos" << std::endl;
1053 std::vector<std::pair<calib::CalibElectron*,calib::CalibElectron*> > zeeCandidates;
1057 double DeltaMinvMin(5000.);
1059 if (calibElectrons.size() < 2)
1062 for(
unsigned int e_it = 0 ; e_it != calibElectrons.size() - 1 ; e_it++){
1063 for(
unsigned int p_it = e_it + 1 ; p_it != calibElectrons.size() ; p_it++)
1066 std::cout << e_it <<
" " << calibElectrons[e_it].getRecoElectron()->charge() <<
" " << p_it <<
" " << calibElectrons[p_it].getRecoElectron()->charge() << std::endl;
1068 if (calibElectrons[e_it].getRecoElectron()->
charge() * calibElectrons[p_it].getRecoElectron()->
charge() != -1)
1077 std::cout <<
"#######################mass "<<
mass << std::endl;
1080 zeeCandidates.push_back(std::pair<calib::CalibElectron*,calib::CalibElectron*>(&(calibElectrons[e_it]),&(calibElectrons[p_it])));
1081 double DeltaMinv = fabs(
mass -
MZ);
1083 if( DeltaMinv < DeltaMinvMin)
1085 DeltaMinvMin = DeltaMinv;
1086 myBestZ=zeeCandidates.size()-1;
1097 if(zeeCandidates.empty() || myBestZ==-1 )
1104 std::cout <<
"Found ZCandidates " << myBestZ << std::endl;
1112 int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1113 int class2 = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1144 if( class1==10 || class1 ==20)
1146 if( class1==110 || class1 ==120)
1148 if( class1>=30 && class1 <=34)
1150 if( class1>=130 && class1 <=134)
1161 if( class2==10 || class2 ==20)
1163 if( class2==110 || class2 ==120)
1165 if( class2>=30 && class2 <=34)
1167 if( class2>=130 && class2 <=134)
1182 bool firstElectronIsOnModuleBorder(
false);
1183 bool secondElectronIsOnModuleBorder(
false);
1194 if( firstElehottestDetId.
subdetId() ==
EcalBarrel && !firstElectronIsOnModuleBorder )
1201 if( secondElehottestDetId.subdetId() ==
EcalBarrel)
1206 if( secondElehottestDetId.subdetId() ==
EcalBarrel && !secondElectronIsOnModuleBorder )
1213 if ( firstElehottestDetId.
subdetId() ==
EcalBarrel && firstElectronIsOnModuleBorder ){
1220 if ( secondElehottestDetId.subdetId() ==
EcalBarrel && secondElectronIsOnModuleBorder ){
1230 if(class1<100 && class2<100){
1232 if(class1==0 && class2==0)
BBZN_gg++;
1233 if(class1<21 && class2<21)
BBZN_tt++;
1234 if(class1<21 || class2<21)
BBZN_t0++;
1238 if(class1>=100 && class2>=100){
1240 if(class1==100 && class2==100)
EEZN_gg++;
1241 if(class1<121 && class2<121)
EEZN_tt++;
1242 if(class1<121 || class2<121)
EEZN_t0++;
1246 if( (class1<100 && class2>=100) || (class2<100 && class1>=100)){
1248 if( (class1==0 && class2==100)||(class2==0 && class1==100) )
EBZN_gg++;
1249 if( ( class1<21 && class2<121) ||(class2<21 && class1<121) )
EBZN_tt++;
1250 if( class2<21 || class1<21 || class2<121 || class1<121 )
EBZN_t0++;
1262 bool selectionBool=
false;
1267 float theta1 = 2. * atan(
exp(- zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->eta()) );
1268 bool ET_1 = ( (zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->energy() *
sin( theta1) ) > 20.);
1270 float theta2 = 2. * atan(
exp(- zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->eta()) );
1271 bool ET_2 = ( (zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->energy() *
sin( theta2) ) > 20.);
1274 bool HoE_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->hadronicOverEm() < 0.115);
1275 bool HoE_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->hadronicOverEm() < 0.115);
1277 bool DeltaPhiIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1278 bool DeltaPhiIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1280 bool DeltaEtaIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1281 bool DeltaEtaIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1285 if(! (invMassBool &&
1288 DeltaPhiIn_1 && DeltaPhiIn_2 &&
1289 DeltaEtaIn_1 && DeltaEtaIn_2
1299 int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1300 int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1406 float ele1EnergyCorrection(1.);
1407 float ele2EnergyCorrection(1.);
1416 if (invMassBool && selectionBool)
1429 std::map<HepMC::GenParticle*,const reco::GsfElectron*> myMCmap;
1431 std::vector<const reco::GsfElectron*> dauElectronCollection;
1433 dauElectronCollection.push_back(zeeCandidates[myBestZ].
first->getRecoElectron() );
1434 dauElectronCollection.push_back(zeeCandidates[myBestZ].
second->getRecoElectron() );
1436 fillMCmap(&dauElectronCollection,mcEle,myMCmap);
1448 int c1st = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1449 int c2nd = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1450 if(c1st<100 && c2nd<100 )
1454 if(c1st>=100 && c2nd>=100 )
1471 std::cout <<
"Added event to algorithm" << std::endl;
1483 std::cout<<
"[ZeeCalibration] Starting loop number " << iLoop<<std::endl;
1492 std::cout<<
"[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1522 std::cout<<
"[ZeeCalibration] Ending loop " << iLoop<<std::endl;
1534 std::cout<<
"Optimized coefficients " << optimizedCoefficients.size() <<std::endl;
1540 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1554 std::cout<<
"size "<<optimizedCoefficients.size()<<std::endl;
1556 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1558 calibCoeff[ieta] *= optimizedCoefficients[ieta];
1564 std::cout<< ieta <<
" " << optimizedCoefficients[ieta] <<std::endl;
1567 std::vector<DetId> ringIds;
1579 for (
unsigned int iid=0; iid<ringIds.size();++iid){
1582 EBDetId myEBDetId(ringIds[iid]);
1588 EEDetId myEEDetId(ringIds[iid]);
1589 if(myEEDetId.
zside() < 0)
1592 if(myEEDetId.
zside() > 0)
1598 ical->setValue( ringIds[iid], *(
ical->getMap().find(ringIds[iid]) ) * optimizedCoefficients[ieta] );
1635 double parResidual[3];
1636 double errparResidual[3];
1648 std::cout<<
"Fit on residuals, sigma is "<<parResidual[2]<<
" +/- "<<errparResidual[2]<<std::endl;
1687 h1_seedOverSC_=
new TH1F(
"h1_seedOverSC",
"h1_seedOverSC", 400, 0., 2.);
1694 h2_fEtaBarrelGood_ =
new TH2F(
"fEtaBarrelGood",
"fEtaBarrelGood",800,-4.,4.,800,0.8,1.2);
1698 h2_fEtaBarrelBad_ =
new TH2F(
"fEtaBarrelBad",
"fEtaBarrelBad",800,-4.,4.,800,0.8,1.2);
1702 h2_fEtaEndcapGood_ =
new TH2F(
"fEtaEndcapGood",
"fEtaEndcapGood",800,-4.,4.,800,0.8,1.2);
1706 h2_fEtaEndcapBad_ =
new TH2F(
"fEtaEndcapBad",
"fEtaEndcapBad",800,-4.,4.,800,0.8,1.2);
1710 for (
int i=0;
i<2;
i++)
1714 sprintf(histoName,
"h_eleEffEta_%d",
i);
1715 h_eleEffEta_[
i] =
new TH1F(histoName,histoName, 150, 0., 2.7);
1718 sprintf(histoName,
"h_eleEffPhi_%d",
i);
1719 h_eleEffPhi_[
i] =
new TH1F(histoName,histoName, 400, -4., 4.);
1722 sprintf(histoName,
"h_eleEffPt_%d",
i);
1723 h_eleEffPt_[
i] =
new TH1F(histoName,histoName, 200, 0., 200.);
1737 for (
int i=0;
i<25;
i++)
1741 sprintf(histoName,
"h_ESCEtrueVsEta_%d",
i);
1747 sprintf(histoName,
"h_ESCEtrue_%d",
i);
1749 h_ESCEtrue_[
i] =
new TH1F(histoName,histoName, 300,0.,1.5);
1751 sprintf(histoName,
"h2_chi2_%d",
i);
1752 h2_chi2_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 5);
1754 sprintf(histoName,
"h2_iterations_%d",
i);
1755 h2_iterations_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 15);
1757 sprintf(histoName,
"h_ESCcorrEtrueVsEta_%d",
i);
1763 sprintf(histoName,
"h_ESCcorrEtrue_%d",
i);
1767 sprintf(histoName,
"h2_xtalRecalibCoeffBarrel_%d",
i);
1773 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapMinus_%d",
i);
1778 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapPlus_%d",
i);
1804 h1_zMassResol_ =
new TH1F(
"zMassResol",
"zMassResol", 200, -50., 50.);
1808 h1_eleEtaResol_ =
new TH1F(
"eleEtaResol",
"eleEtaResol", 100, -0.01, 0.01);
1824 h1_elePhiResol_ =
new TH1F(
"elePhiResol",
"elePhiResol", 100, -0.01, 0.01);
1829 h1_zEtaResol_ =
new TH1F(
"zEtaResol",
"zEtaResol", 200, -1., 1.);
1834 h1_zPhiResol_ =
new TH1F(
"zPhiResol",
"zPhiResol", 200, -1., 1.);
1838 h1_nEleReco_ =
new TH1F(
"nEleReco",
"Number of reco electrons",10,-0.5,10.5);
1851 h1_occupancy_ =
new TH1F(
"occupancy",
"occupancy",1000,0,10000);
1876 h1_ZCandMult_ =
new TH1F(
"ZCandMult",
"Multiplicity of Z candidates in one event",10,-0.5,10.5);
1879 h1_reco_ZMass_ =
new TH1F(
"reco_ZMass",
"Inv. mass of 2 reco Electrons",200,0.,150.);
1883 h1_reco_ZMassCorr_ =
new TH1F(
"reco_ZMassCorr",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1887 h1_reco_ZMassCorrBB_ =
new TH1F(
"reco_ZMassCorrBB",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1891 h1_reco_ZMassCorrEE_ =
new TH1F(
"reco_ZMassCorrEE",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1897 h2_coeffVsEta_=
new TH2F(
"h2_calibCoeffVsEta",
"h2_calibCoeffVsEta",249,-124,125, 200, 0., 2.);
1901 h2_coeffVsEtaGrouped_=
new TH2F(
"h2_calibCoeffVsEtaGrouped",
"h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1905 h2_zMassVsLoop_=
new TH2F(
"h2_zMassVsLoop",
"h2_zMassVsLoop",1000,0,40, 90, 80.,95.);
1907 h2_zMassDiffVsLoop_=
new TH2F(
"h2_zMassDiffVsLoop",
"h2_zMassDiffVsLoop",1000,0,40, 100, -1., 1.);
1911 h2_zWidthVsLoop_=
new TH2F(
"h2_zWidthVsLoop",
"h2_zWidthVsLoop",1000,0,40, 100, 0.,10.);
1913 h2_coeffVsLoop_=
new TH2F(
"h2_coeffVsLoop",
"h2_coeffVsLoop",1000,0,40, 100, 0., 2.);
1915 h2_residualSigma_=
new TH2F(
"h2_residualSigma",
"h2_residualSigma",1000, 0, 40, 100, 0., .5);
1917 h2_miscalRecal_ =
new TH2F(
"h2_miscalRecal",
"h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1921 h2_miscalRecalEB_ =
new TH2F(
"h2_miscalRecalEB",
"h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1925 h2_miscalRecalEE_ =
new TH2F(
"h2_miscalRecalEE",
"h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1929 h1_mc_ =
new TH1F(
"h1_residualMiscalib",
"h1_residualMiscalib", 200, -0.2, 0.2);
1930 h1_mcEB_ =
new TH1F(
"h1_residualMiscalibEB",
"h1_residualMiscalibEB", 200, -0.2, 0.2);
1931 h1_mcEE_ =
new TH1F(
"h1_residualMiscalibEE",
"h1_residualMiscalibEE", 200, -0.2, 0.2);
1933 for (
int i=0;
i<25;
i++)
1953 sprintf(histoName,
"h1_residualMiscalibParz_%d",
i);
1954 h1_mcParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1955 sprintf(histoName,
"h1_residualMiscalibEBParz_%d",
i);
1956 h1_mcEBParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1957 sprintf(histoName,
"h1_residualMiscalibEEParz_%d",
i);
1958 h1_mcEEParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1968 float p0 = 1.00153e+00;
1969 float p1 = 3.29331e-02;
1970 float p2 = 1.21187e-03;
1972 double x = (double) fabs(scEta);
1974 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
1982 float p0 = 1.06819e+00;
1983 float p1 = -1.53189e-02;
1984 float p2 = 4.01707e-04 ;
1986 double x = (double) fabs(scEta);
1988 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
1994 float p0 = 1.17382e+00;
1995 float p1 = -6.52319e-02;
1996 float p2 = 6.26108e-03;
1998 double x = (double) fabs(scEta);
2000 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2006 float p0 = 9.99782e-01 ;
2007 float p1 = 1.26983e-02;
2008 float p2 = 2.16344e-03;
2010 double x = (double) fabs(scEta);
2012 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2021 for (
unsigned int i=0;
i<mcEle.size();
i++)
2025 for (
unsigned int j=0;j<electronCollection->size();j++)
2027 float dr=
EvalDR(mcEle[
i]->momentum().pseudoRapidity(),(*(*electronCollection)[j]).
eta(),mcEle[
i]->momentum().
phi(),(*(*electronCollection)[j]).
phi());
2030 myMatchEle = (*electronCollection)[j];
2034 myMCmap.insert(std::pair<HepMC::GenParticle*,const reco::GsfElectron*>(mcEle[
i],myMatchEle));
2042 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2043 float DPhi = Phi - Phi_ref ;
2046 float DEta = Eta - Eta_ref ;
2048 float DR =
sqrt( DEta*DEta + DPhi*DPhi );
2055 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2056 return (Phi - Phi_ref);
2062 for (
unsigned int i=0;
i<mcEle.size();
i++)
2065 h_eleEffEta_[0]->Fill(fabs(mcEle[
i]->momentum().pseudoRapidity()));
2069 std::map<HepMC::GenParticle*,const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[
i]);
2070 if (mIter == associationMap.end() )
2077 h_eleEffEta_[1]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
2112 if(k>=0 && k<=84)index = k - 85;
2114 if(k>=85 && k<=169)index = k - 84;
2116 if(k>=170 && k<=208)index = - k + 84;
2118 if(k>=209 && k<=247)index = k - 123;
2124 if(k>=0 && k<=71)index = k - 72;
2126 if(k>=72 && k<=143)index = k - 71;
2136 double correction(1.);
2173 std::pair<DetId, double> myPair (
DetId(0), -9999.);
2176 for( std::vector<std::pair<DetId,float> >::const_iterator idIt=mySCRecHits.begin(); idIt != mySCRecHits.end(); idIt++){
2180 hottestRecHit = & (* ( ebhits->
find((*idIt).first) ) );
2182 if( hottestRecHit == & (*( ebhits->
end())) )
2184 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2188 else if (idIt->first.subdetId() ==
EcalEndcap )
2190 hottestRecHit = & (* ( eehits->
find((*idIt).first) ) );
2191 if( hottestRecHit == & (*( eehits->
end())) )
2193 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2202 maxEnergy = hottestRecHit->
energy();
2204 myPair.first = hottestRecHit ->
id();
2222 short ieta = myEBDetId.
ieta();
2223 short iphi = myEBDetId.
iphi();
2227 myBool = (
abs( ieta ) == 1 ||
abs( ieta ) == 25
2228 ||
abs( ieta ) ==26 ||
abs( ieta ) == 45
2229 ||
abs( ieta ) ==46 ||
abs( ieta ) == 65
2230 ||
abs( ieta ) ==66 ||
abs( ieta ) == 85 );
2232 for(
int i = 0;
i < 19;
i++){
2234 if(iphi == ( 20*
i + 1 ) || iphi == 20*
i )
2250 bool isNearCrack =
false;
2262 dist +=
pow( v1[
i]-v2[
i], 2 );
2318 for (
int i=0;
i<2;
i++)
2376 std::ofstream
fout(
"ZeeStatistics.txt");
2379 std::cout <<
"Cannot open output file.\n";
2382 fout<<
"ZeeStatistics"<<std::endl;
2384 fout<<
"##########################RECO#########################"<<std::endl;
2385 fout<<
"##################Zee with Barrel-Barrel electrons: "<<
BBZN<<std::endl;
2388 fout<<
"##################Zee with Barrel-Endcap electrons: "<<
EBZN<<std::endl;
2390 fout<<
"##################Zee with Endcap-Endcap electrons: "<<
EEZN<<std::endl;
2393 fout<<
"\n"<<std::endl;
2395 fout<<
"##########################GEN#########################"<<std::endl;
2396 fout<<
"##################Zee with Barrel-Barrel electrons: "<<(
float)
MCZBB/
NEVT<<std::endl;
2397 fout<<
"##################Zee with Barrel-Endcap electrons: "<<(
float)
MCZEB/
NEVT<<std::endl;
2398 fout<<
"##################Zee with Endcap-Endcap electrons: "<<(
float)
MCZEE/
NEVT<<std::endl;
TH1 * weightedRescaleFactor[50][250]
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
double eta() const final
momentum pseudorapidity
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)
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?
constexpr uint32_t rawId() const
get the raw id
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()
Status duringLoop(const edm::Event &, const edm::EventSetup &) override
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_
~ZeeCalibration() override
Destructor.
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
U second(std::pair< T, U > const &p)
TH1F * h1_electronCosTheta_TK_
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 beginOfJob() override
Called at beginning of job.
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
int getNumberOfIterations() const
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
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_
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_
void startingNewLoop(unsigned int iLoop) override
Called at beginning of loop.
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
Status endOfLoop(const edm::EventSetup &, unsigned int iLoop) override
Called at end of loop.
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.
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)
SuperClusterRef superCluster() const override
reference to a SuperCluster
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 *)
double phi() const final
momentum azimuthal angle
TH2F * h2_fEtaBarrelGood_
std::string erechitCollection_
Power< A, B >::type pow(const A &a, const B &b)
TH1F * h1_borderElectronClassification_
TH2F * h2_xtalMiscalibCoeffEndcapMinus_
void endOfJob() override
Called at end of job.
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