18 #include "CLHEP/Vector/LorentzVector.h"
64 #include "TGraphErrors.h"
83 std::cout<<
"[ZeeCalibration] Starting the ctor"<<std::endl;
115 outputFile_ = TFile::Open(outputFileName_.c_str(),
"RECREATE");
117 myTree =
new TTree(
"myTree",
"myTree");
119 myTree->Branch(
"zMass",&
mass4tree,
"mass/F");
158 findingRecord<EcalIntercalibConstantsRcd> () ;
160 for(
int i = 0;
i<50;
i++){
170 std::cout<<
"[ZeeCalibration] Done with the ctor"<<std::endl;
184 boost::shared_ptr<EcalIntercalibConstants>
187 std::cout <<
"@SUB=ZeeCalibration::produceEcalIntercalibConstants" << std::endl;
210 if(ieta==0)
continue;
244 std::cout<<
"Writing histos..."<<std::endl;
273 for(
unsigned int i =0;
i<25;
i++){
310 Double_t
mean[25] = {0.};
311 Double_t
num[25] = {0.};
312 Double_t meanErr[25] = {0.};
313 Float_t
rms[25] = {0.};
314 Float_t tempRms[10][25];
316 for(
int ia = 0; ia<10; ia++){
317 for(
int ib = 0;
ib<25;
ib++){
319 tempRms[ia][
ib] = 0.;
332 bool isNearCrack =
false;
397 if(
k>=170 &&
k<=204){
437 if(
k>=205 &&
k<=208){
485 for(
int ic = 0; ic< 17; ic++){
487 mean[ic] = mean[ic] / num[ic];
489 meanErr[ic] = 1. /
sqrt(meanErr[ic]);
495 for(
int ic = 0; ic< 25; ic++){
496 for(
int id = 0;
id< 10;
id++){
498 if(tempRms[
id][ic] > 0.){
500 rms[ic] += (tempRms[id][ic] - mean[
j])*(tempRms[
id][ic] - mean[j]);
505 rms[ic] =
sqrt(rms[ic]);
512 Double_t xtalEta[25] = {1.4425, 1.3567,1.2711,1.1855,
514 0.7468,0.6612,0.5756,0.4897,0.3985,0.3117,0.2250,0.1384,0.0487,
515 1.546, 1.651, 1.771, 1.908, 2.071, 2.267, 2.516, 2.8};
517 Double_t
zero[25] = {0.026};
519 for(
int j = 0; j <25; j++)
526 px->SetXTitle(
"Eta channel");
527 px->SetYTitle(
"recalibCoeff");
547 double weightSumMeanBarrel = 0.;
548 double weightSumMeanEndcap = 0.;
586 std::cout<<
"Weight sum mean on channels in Barrel is :"<<weightSumMeanBarrel<<std::endl;
587 std::cout<<
"Weight sum mean on channels in Endcap is :"<<weightSumMeanEndcap<<std::endl;
599 TGraphErrors*
graph =
new TGraphErrors(25,xtalEta,mean,zero,meanErr);
603 double zero50[50] = { 0. };
606 residualSigmaGraph->SetName(
"residualSigmaGraph");
607 residualSigmaGraph->Draw(
"APL");
608 residualSigmaGraph->Write();
611 coefficientDistanceAtIterationGraph->SetName(
"coefficientDistanceAtIterationGraph");
612 coefficientDistanceAtIterationGraph->Draw(
"APL");
613 coefficientDistanceAtIterationGraph->Write();
615 Float_t noError[250] = {0.};
617 Float_t ringInd[250];
618 for(
int i =0;
i<250;
i++)
622 graphCoeff->SetName(
"graphCoeff");
623 graphCoeff->Draw(
"APL");
662 std::cout<<
"[ZeeCalibration] Entering duringLoop"<<std::endl;
681 std::cout<<
"[ZeeCalibration::beginOfJob] Histograms booked "<<std::endl;
699 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EB miscal file"<<std::endl;
708 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EE miscal file"<<std::endl;
712 std::cout <<
" theAlgorithm_->getNumberOfChannels() "
722 std::vector<DetId> ringIds;
736 for (
unsigned int iid=0; iid<ringIds.size();++iid)
738 float miscalib=* (miscalibMap->
get().
getMap().
find(ringIds[iid]) );
757 std::vector<DetId> ringIds;
769 for (
unsigned int iid=0; iid<ringIds.size();++iid){
773 EBDetId myEBDetId(ringIds[iid]);
779 EEDetId myEEDetId(ringIds[iid]);
780 if(myEEDetId.
zside() < 0)
783 if(myEEDetId.
zside() > 0)
788 ical->setValue( ringIds[iid], *(miscalibMap->
get().
getMap().
find(ringIds[iid]) ) );
811 for(
unsigned int iHLT=0; iHLT<200; ++iHLT) {
816 std::cout<<
"[ZeeCalibration::duringLoop] Done with initializing aHLTresults[] "<<std::endl;
822 if(!hltTriggerResultHandle.isValid()) {
826 hltCount = hltTriggerResultHandle->size();
832 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillHLTInfo(hltTriggerResultHandle); "<<std::endl;
850 std::cout<<
"[ZeeCalibration::duringLoop] End HLT section"<<std::endl;
856 std::vector<HepMC::GenParticle*> mcEle;
858 float myGenZMass(-1);
868 const HepMC::GenEvent * myGenEvent = hepProd->GetEvent();
874 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillZMCInfo( & (*myGenEvent) ); "<<std::endl;
878 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
879 p != myGenEvent->particles_end(); ++
p ) {
881 if ( (*p)->pdg_id() == 23 && (*p)->status()==2){
883 myGenZMass = (*p)->momentum().m();
896 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
897 p != myGenEvent->particles_end(); ++
p ) {
899 if (
abs( (*p)->pdg_id() ) == 11 )
901 mcEle.push_back( (*
p) );
908 if(mcEle.size()==2 && fabs(mcEle[0]->momentum().
eta())<2.4 && fabs(mcEle[1]->momentum().
eta())<2.4 ){
911 if( fabs(mcEle[0]->momentum().
eta())<1.479 && fabs(mcEle[1]->momentum().
eta())<1.479 )
MCZBB++;
913 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++;
915 if( fabs(mcEle[0]->momentum().eta())>1.479 && fabs(mcEle[1]->momentum().eta())>1.479 )
MCZEE++;
929 std::cerr <<
"Error! can't get the product EBRecHitCollection " << std::endl;
938 std::cerr <<
"Error! can't get the product EERecHitCollection " << std::endl;
948 std::cerr <<
"Error! can't get the product SuperClusterCollection "<< std::endl;
953 std::cout<<
"scCollection->size()"<<scCollection->size()<<std::endl;
954 for(reco::SuperClusterCollection::const_iterator scIt = scCollection->begin(); scIt != scCollection->end(); scIt++)
965 std::cerr <<
"Error! can't get the product IslandSuperClusterCollection "<< std::endl;
970 std::cout<<
"scCollection->size()"<<scIslandCollection->size()<<std::endl;
973 if( ( scCollection->size()+scIslandCollection->size() ) < 2)
981 std::cerr <<
"Error! can't get the product ElectronCollection "<< std::endl;
994 if(electronCollection->size() < 2)
997 if ( !hits && !ehits){
1002 if (hits->size() == 0 && ehits->size() == 0){
1003 std::cout <<
"hits->size() == 0" << std::endl;
1007 if (!electronCollection){
1008 std::cout <<
"!electronCollection" << std::endl;
1012 if (electronCollection->size() == 0){
1013 std::cout <<
"electronCollection->size() == 0" << std::endl;
1030 std::cout <<
" Starting with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1037 std::cout <<
" Done with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1043 std::vector<calib::CalibElectron> calibElectrons;
1049 for(
unsigned int e_it = 0 ; e_it != electronCollection->size() ; e_it++)
1053 std::cout << calibElectrons.back().getRecoElectron()->superCluster()->energy() <<
" " << calibElectrons.back().getRecoElectron()->energy() << std::endl;
1061 std::cout <<
"Filled histos" << std::endl;
1065 std::vector<std::pair<calib::CalibElectron*,calib::CalibElectron*> > zeeCandidates;
1069 double DeltaMinvMin(5000.);
1071 if (calibElectrons.size() < 2)
1074 for(
unsigned int e_it = 0 ; e_it != calibElectrons.size() - 1 ; e_it++){
1075 for(
unsigned int p_it = e_it + 1 ; p_it != calibElectrons.size() ; p_it++)
1078 std::cout << e_it <<
" " << calibElectrons[e_it].getRecoElectron()->charge() <<
" " << p_it <<
" " << calibElectrons[p_it].getRecoElectron()->charge() << std::endl;
1080 if (calibElectrons[e_it].getRecoElectron()->
charge() * calibElectrons[p_it].getRecoElectron()->
charge() != -1)
1089 std::cout <<
"#######################mass "<<
mass << std::endl;
1092 zeeCandidates.push_back(std::pair<calib::CalibElectron*,calib::CalibElectron*>(&(calibElectrons[e_it]),&(calibElectrons[p_it])));
1093 double DeltaMinv = fabs(
mass -
MZ);
1095 if( DeltaMinv < DeltaMinvMin)
1097 DeltaMinvMin = DeltaMinv;
1098 myBestZ=zeeCandidates.size()-1;
1109 if(zeeCandidates.size()==0 || myBestZ==-1 )
1116 std::cout <<
"Found ZCandidates " << myBestZ << std::endl;
1127 int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1128 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)
1190 DetId firstElehottestDetId =
getHottestDetId( zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->seed()->hitsAndFractions() , hits, ehits ).
first;
1191 DetId secondElehottestDetId =
getHottestDetId( zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->seed()->hitsAndFractions() , hits, ehits ).
first;
1193 bool firstElectronIsOnModuleBorder(
false);
1194 bool secondElectronIsOnModuleBorder(
false);
1200 if( firstElehottestDetId.subdetId() ==
EcalBarrel)
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
1312 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1313 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1314 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 40 &&
1315 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 140);
1320 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==0 ||
1321 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==10 ||
1322 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==20 ||
1323 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==100 ||
1324 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==110 ||
1325 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==120
1327 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() == 0 ||
1328 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 10 ||
1329 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 20 ||
1330 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 100 ||
1331 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 110 ||
1332 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 120
1338 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() == 0 ||
1339 zeeCandidates[myBestZ].first->getRecoElectron()->classification() == 100
1341 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() == 0 ||
1342 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 100
1348 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() >=30 &&
1349 zeeCandidates[myBestZ].first->getRecoElectron()->classification() <=34)
1351 ((zeeCandidates[myBestZ].
first->getRecoElectron()->classification() >=130 &&
1352 zeeCandidates[myBestZ].first->getRecoElectron()->classification() <=134))
1355 ( (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() >=30 &&
1356 zeeCandidates[myBestZ].second->getRecoElectron()->classification() <=34)
1358 ((zeeCandidates[myBestZ].
second->getRecoElectron()->classification() >=130 &&
1359 zeeCandidates[myBestZ].second->getRecoElectron()->classification() <=134))
1370 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==0 ||
1371 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==10 ||
1372 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==20
1373 ) && zeeCandidates[myBestZ].
second->getRecoElectron()->classification()>=100
1374 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()!=140
1380 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() ==0 ||
1381 zeeCandidates[myBestZ].second->getRecoElectron()->classification() ==10 ||
1382 zeeCandidates[myBestZ].second->getRecoElectron()->classification() ==20
1383 ) && zeeCandidates[myBestZ].
first->getRecoElectron()->classification()>=100
1384 && zeeCandidates[myBestZ].first->getRecoElectron()->classification()!=140
1394 zeeCandidates[myBestZ].first->getRecoElectron()->classification()>=100 &&
1395 zeeCandidates[myBestZ].second->getRecoElectron()->classification()>= 100 &&
1396 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 140 &&
1397 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 140);
1402 zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 &&
1403 zeeCandidates[myBestZ].second->getRecoElectron()->classification()< 100 &&
1404 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1405 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 40);
1410 !(zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 &&
1411 zeeCandidates[myBestZ].second->getRecoElectron()->classification()>=100) &&
1412 !(zeeCandidates[myBestZ].
first->getRecoElectron()->classification()>=100 &&
1413 zeeCandidates[myBestZ].second->getRecoElectron()->classification()<100) );
1416 float ele1EnergyCorrection(1.);
1417 float ele2EnergyCorrection(1.);
1426 if (invMassBool && selectionBool)
1439 std::map<HepMC::GenParticle*,const reco::GsfElectron*> myMCmap;
1441 std::vector<const reco::GsfElectron*> dauElectronCollection;
1443 dauElectronCollection.push_back(zeeCandidates[myBestZ].
first->getRecoElectron() );
1444 dauElectronCollection.push_back(zeeCandidates[myBestZ].
second->getRecoElectron() );
1446 fillMCmap(&dauElectronCollection,mcEle,myMCmap);
1458 if(zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()<100 )
1462 if(zeeCandidates[myBestZ].first->getRecoElectron()->classification()>=100 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()>=100 )
1479 std::cout <<
"Added event to algorithm" << std::endl;
1491 std::cout<<
"[ZeeCalibration] Starting loop number " << iLoop<<std::endl;
1500 std::cout<<
"[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1530 std::cout<<
"[ZeeCalibration] Ending loop " << iLoop<<std::endl;
1542 std::cout<<
"Optimized coefficients " << optimizedCoefficients.size() <<std::endl;
1548 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1562 std::cout<<
"size "<<optimizedCoefficients.size()<<std::endl;
1564 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1566 calibCoeff[ieta] *= optimizedCoefficients[ieta];
1572 std::cout<< ieta <<
" " << optimizedCoefficients[ieta] <<std::endl;
1575 std::vector<DetId> ringIds;
1587 for (
unsigned int iid=0; iid<ringIds.size();++iid){
1590 EBDetId myEBDetId(ringIds[iid]);
1596 EEDetId myEEDetId(ringIds[iid]);
1597 if(myEEDetId.
zside() < 0)
1600 if(myEEDetId.
zside() > 0)
1606 ical->setValue( ringIds[iid], *(
ical->getMap().find(ringIds[iid]) ) * optimizedCoefficients[ieta] );
1643 double parResidual[3];
1644 double errparResidual[3];
1656 std::cout<<
"Fit on residuals, sigma is "<<parResidual[2]<<
" +/- "<<errparResidual[2]<<std::endl;
1695 h1_seedOverSC_=
new TH1F(
"h1_seedOverSC",
"h1_seedOverSC", 400, 0., 2.);
1702 h2_fEtaBarrelGood_ =
new TH2F(
"fEtaBarrelGood",
"fEtaBarrelGood",800,-4.,4.,800,0.8,1.2);
1706 h2_fEtaBarrelBad_ =
new TH2F(
"fEtaBarrelBad",
"fEtaBarrelBad",800,-4.,4.,800,0.8,1.2);
1710 h2_fEtaEndcapGood_ =
new TH2F(
"fEtaEndcapGood",
"fEtaEndcapGood",800,-4.,4.,800,0.8,1.2);
1714 h2_fEtaEndcapBad_ =
new TH2F(
"fEtaEndcapBad",
"fEtaEndcapBad",800,-4.,4.,800,0.8,1.2);
1718 for (
int i=0;
i<2;
i++)
1722 sprintf(histoName,
"h_eleEffEta_%d",
i);
1723 h_eleEffEta_[
i] =
new TH1F(histoName,histoName, 150, 0., 2.7);
1726 sprintf(histoName,
"h_eleEffPhi_%d",
i);
1727 h_eleEffPhi_[
i] =
new TH1F(histoName,histoName, 400, -4., 4.);
1730 sprintf(histoName,
"h_eleEffPt_%d",
i);
1731 h_eleEffPt_[
i] =
new TH1F(histoName,histoName, 200, 0., 200.);
1745 for (
int i=0;
i<25;
i++)
1749 sprintf(histoName,
"h_ESCEtrueVsEta_%d",
i);
1755 sprintf(histoName,
"h_ESCEtrue_%d",
i);
1757 h_ESCEtrue_[
i] =
new TH1F(histoName,histoName, 300,0.,1.5);
1759 sprintf(histoName,
"h2_chi2_%d",
i);
1760 h2_chi2_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 5);
1762 sprintf(histoName,
"h2_iterations_%d",
i);
1763 h2_iterations_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 15);
1765 sprintf(histoName,
"h_ESCcorrEtrueVsEta_%d",
i);
1771 sprintf(histoName,
"h_ESCcorrEtrue_%d",
i);
1775 sprintf(histoName,
"h2_xtalRecalibCoeffBarrel_%d",
i);
1781 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapMinus_%d",
i);
1786 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapPlus_%d",
i);
1812 h1_zMassResol_ =
new TH1F(
"zMassResol",
"zMassResol", 200, -50., 50.);
1816 h1_eleEtaResol_ =
new TH1F(
"eleEtaResol",
"eleEtaResol", 100, -0.01, 0.01);
1832 h1_elePhiResol_ =
new TH1F(
"elePhiResol",
"elePhiResol", 100, -0.01, 0.01);
1837 h1_zEtaResol_ =
new TH1F(
"zEtaResol",
"zEtaResol", 200, -1., 1.);
1842 h1_zPhiResol_ =
new TH1F(
"zPhiResol",
"zPhiResol", 200, -1., 1.);
1846 h1_nEleReco_ =
new TH1F(
"nEleReco",
"Number of reco electrons",10,-0.5,10.5);
1859 h1_occupancy_ =
new TH1F(
"occupancy",
"occupancy",1000,0,10000);
1884 h1_ZCandMult_ =
new TH1F(
"ZCandMult",
"Multiplicity of Z candidates in one event",10,-0.5,10.5);
1887 h1_reco_ZMass_ =
new TH1F(
"reco_ZMass",
"Inv. mass of 2 reco Electrons",200,0.,150.);
1891 h1_reco_ZMassCorr_ =
new TH1F(
"reco_ZMassCorr",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1895 h1_reco_ZMassCorrBB_ =
new TH1F(
"reco_ZMassCorrBB",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1899 h1_reco_ZMassCorrEE_ =
new TH1F(
"reco_ZMassCorrEE",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1905 h2_coeffVsEta_=
new TH2F(
"h2_calibCoeffVsEta",
"h2_calibCoeffVsEta",249,-124,125, 200, 0., 2.);
1909 h2_coeffVsEtaGrouped_=
new TH2F(
"h2_calibCoeffVsEtaGrouped",
"h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1913 h2_zMassVsLoop_=
new TH2F(
"h2_zMassVsLoop",
"h2_zMassVsLoop",1000,0,40, 90, 80.,95.);
1915 h2_zMassDiffVsLoop_=
new TH2F(
"h2_zMassDiffVsLoop",
"h2_zMassDiffVsLoop",1000,0,40, 100, -1., 1.);
1919 h2_zWidthVsLoop_=
new TH2F(
"h2_zWidthVsLoop",
"h2_zWidthVsLoop",1000,0,40, 100, 0.,10.);
1921 h2_coeffVsLoop_=
new TH2F(
"h2_coeffVsLoop",
"h2_coeffVsLoop",1000,0,40, 100, 0., 2.);
1923 h2_residualSigma_=
new TH2F(
"h2_residualSigma",
"h2_residualSigma",1000, 0, 40, 100, 0., .5);
1925 h2_miscalRecal_ =
new TH2F(
"h2_miscalRecal",
"h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1929 h2_miscalRecalEB_ =
new TH2F(
"h2_miscalRecalEB",
"h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1933 h2_miscalRecalEE_ =
new TH2F(
"h2_miscalRecalEE",
"h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1937 h1_mc_ =
new TH1F(
"h1_residualMiscalib",
"h1_residualMiscalib", 200, -0.2, 0.2);
1938 h1_mcEB_ =
new TH1F(
"h1_residualMiscalibEB",
"h1_residualMiscalibEB", 200, -0.2, 0.2);
1939 h1_mcEE_ =
new TH1F(
"h1_residualMiscalibEE",
"h1_residualMiscalibEE", 200, -0.2, 0.2);
1941 for (
int i=0;
i<25;
i++)
1961 sprintf(histoName,
"h1_residualMiscalibParz_%d",
i);
1962 h1_mcParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1963 sprintf(histoName,
"h1_residualMiscalibEBParz_%d",
i);
1964 h1_mcEBParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1965 sprintf(histoName,
"h1_residualMiscalibEEParz_%d",
i);
1966 h1_mcEEParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1976 float p0 = 1.00153e+00;
1977 float p1 = 3.29331e-02;
1978 float p2 = 1.21187e-03;
1980 double x = (double) fabs(scEta);
1982 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
1990 float p0 = 1.06819e+00;
1991 float p1 = -1.53189e-02;
1992 float p2 = 4.01707e-04 ;
1994 double x = (double) fabs(scEta);
1996 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2002 float p0 = 1.17382e+00;
2003 float p1 = -6.52319e-02;
2004 float p2 = 6.26108e-03;
2006 double x = (double) fabs(scEta);
2008 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2014 float p0 = 9.99782e-01 ;
2015 float p1 = 1.26983e-02;
2016 float p2 = 2.16344e-03;
2018 double x = (double) fabs(scEta);
2020 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2027 void ZeeCalibration::fillMCmap(
const std::vector<const reco::GsfElectron*>* electronCollection,
const std::vector<HepMC::GenParticle*>& mcEle,std::map<HepMC::GenParticle*,const reco::GsfElectron*>& myMCmap)
2029 for (
unsigned int i=0;
i<mcEle.size();
i++)
2033 for (
unsigned int j=0;
j<electronCollection->size();
j++)
2035 float dr=
EvalDR(mcEle[
i]->momentum().pseudoRapidity(),(*(*electronCollection)[
j]).
eta(),mcEle[
i]->momentum().
phi(),(*(*electronCollection)[j]).
phi());
2038 myMatchEle = (*electronCollection)[
j];
2042 myMCmap.insert(std::pair<HepMC::GenParticle*,const reco::GsfElectron*>(mcEle[
i],myMatchEle));
2050 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2051 float DPhi = Phi - Phi_ref ;
2054 float DEta = Eta - Eta_ref ;
2056 float DR =
sqrt( DEta*DEta + DPhi*DPhi );
2063 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2064 return (Phi - Phi_ref);
2070 for (
unsigned int i=0;
i<mcEle.size();
i++)
2073 h_eleEffEta_[0]->Fill(fabs(mcEle[
i]->momentum().pseudoRapidity()));
2077 std::map<HepMC::GenParticle*,const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[
i]);
2078 if (mIter == associationMap.end() )
2085 h_eleEffEta_[1]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
2120 if(k>=0 && k<=84)index = k - 85;
2122 if(k>=85 && k<=169)index = k - 84;
2124 if(k>=170 && k<=208)index = - k + 84;
2126 if(k>=209 && k<=247)index = k - 123;
2132 if(k>=0 && k<=71)index = k - 72;
2134 if(k>=72 && k<=143)index = k - 71;
2144 double correction(1.);
2177 double maxEnergy = -9999.;
2180 std::pair<DetId, double> myPair (
DetId(0), -9999.);
2183 for( std::vector<std::pair<DetId,float> >::const_iterator idIt=mySCRecHits.begin(); idIt != mySCRecHits.end(); idIt++){
2187 hottestRecHit = & (* ( ebhits->
find((*idIt).first) ) );
2189 if( hottestRecHit == & (*( ebhits->
end())) )
2191 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2195 else if (idIt->first.subdetId() ==
EcalEndcap )
2197 hottestRecHit = & (* ( eehits->
find((*idIt).first) ) );
2198 if( hottestRecHit == & (*( eehits->
end())) )
2200 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2207 if(hottestRecHit && hottestRecHit->
energy() > maxEnergy){
2209 maxEnergy = hottestRecHit->
energy();
2211 myPair.first = hottestRecHit ->
id();
2212 myPair.second = maxEnergy;
2229 short ieta = myEBDetId.
ieta();
2230 short iphi = myEBDetId.
iphi();
2234 myBool = (
abs( ieta ) == 1 ||
abs( ieta ) == 25
2235 ||
abs( ieta ) ==26 ||
abs( ieta ) == 45
2236 ||
abs( ieta ) ==46 ||
abs( ieta ) == 65
2237 ||
abs( ieta ) ==66 ||
abs( ieta ) == 85 );
2239 for(
int i = 0;
i < 19;
i++){
2241 if(iphi == ( 20*
i + 1 ) || iphi == 20*
i )
2257 bool isNearCrack =
false;
2269 dist +=
pow( v1[
i]-v2[
i], 2 );
2325 for (
int i=0;
i<2;
i++)
2383 ofstream
fout(
"ZeeStatistics.txt");
2386 std::cout <<
"Cannot open output file.\n";
2389 fout<<
"ZeeStatistics"<<std::endl;
2391 fout<<
"##########################RECO#########################"<<std::endl;
2392 fout<<
"##################Zee with Barrel-Barrel electrons: "<<
BBZN<<std::endl;
2395 fout<<
"##################Zee with Barrel-Endcap electrons: "<<
EBZN<<std::endl;
2397 fout<<
"##################Zee with Endcap-Endcap electrons: "<<
EEZN<<std::endl;
2400 fout<<
"\n"<<std::endl;
2402 fout<<
"##########################GEN#########################"<<std::endl;
2403 fout<<
"##################Zee with Barrel-Barrel electrons: "<<(float)
MCZBB/
NEVT<<std::endl;
2404 fout<<
"##################Zee with Barrel-Endcap electrons: "<<(float)
MCZEB/
NEVT<<std::endl;
2405 fout<<
"##################Zee with Endcap-Endcap electrons: "<<(float)
MCZEE/
NEVT<<std::endl;
virtual void endOfJob()
Called at end of job.
T getParameter(std::string const &) const
float NewCalibCoeff[nMaxChannels]
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.
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)
std::pair< DetId, double > getHottestDetId(std::vector< std::pair< DetId, float > > mySCRecHits, const EBRecHitCollection *ebhits, const EERecHitCollection *eehits)
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)
float calibCoeffError[nMaxChannels]
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
U second(std::pair< T, U > const &p)
TH1F * h1_electronCosTheta_TK_
~ZeeCalibration()
Destructor.
TH2F * h_ESCcorrEtrueVsEta_[25]
virtual float phi() const GCC11_FINAL
momentum azimuthal angle
boost::shared_ptr< EcalIntercalibConstants > ical
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
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_
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]
virtual float eta() const GCC11_FINAL
momentum pseudorapidity
const_iterator end() const
TH2F * h2_zMassDiffVsLoop_
std::string scCollection_
const std::vector< float > & getOptimizedCoefficientsError() const
static const int MAX_IPHI
DetId id() const
get the id
void writeEleHistograms()
edm::ParameterSet theParameterSet
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Classification classification() const
TH1F * h_ESCcorrEtrue_[25]
TH2F * h2_xtalMiscalibCoeffEndcapPlus_
static const int MAX_IETA
virtual boost::shared_ptr< EcalIntercalibConstants > produceEcalIntercalibConstants(const EcalIntercalibConstantsRcd &iRecord)
Produce Ecal interCalibrations.
float calibCoeff[nMaxChannels]
int SILVER_ELECTRONS_IN_ENDCAP
EcalIntercalibConstantMap EcalIntercalibConstants
TH1F * h1_reco_ZMassCorrBB_
void writeLine(EBDetId const &, float)
iterator find(key_type k)
T perp() const
Magnitude of transverse component.
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
TH2F * h2_xtalRecalibCoeffEndcapPlus_[25]
int CRACK_ELECTRONS_IN_ENDCAP
int CRACK_ELECTRONS_IN_BARREL
TH1 * weightedRescaleFactor[nMaxIterations][nMaxChannels]
TH2F * h2_fEtaEndcapGood_
TH2F * h2_xtalRecalibCoeffEndcapMinus_[25]
void bookEleMCHistograms()
std::string electronProducer_
double getEtaCorrection(const reco::GsfElectron *)
TH2F * h2_fEtaBarrelGood_
std::string erechitCollection_
tuple size
Write out results.
Power< A, B >::type pow(const A &a, const B &b)
TH1F * h1_borderElectronClassification_
float initCalibCoeff[nMaxChannels]
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_