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 boost::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]) );
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()) {
825 hltCount = hltTriggerResultHandle->size();
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;
1126 int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1127 int class2 = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1154 if( class1==10 || class1 ==20)
1156 if( class1==110 || class1 ==120)
1158 if( class1>=30 && class1 <=34)
1160 if( class1>=130 && class1 <=134)
1171 if( class2==10 || class2 ==20)
1173 if( class2==110 || class2 ==120)
1175 if( class2>=30 && class2 <=34)
1177 if( class2>=130 && class2 <=134)
1189 DetId firstElehottestDetId =
getHottestDetId( zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->seed()->hitsAndFractions() , hits, ehits ).
first;
1190 DetId secondElehottestDetId =
getHottestDetId( zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->seed()->hitsAndFractions() , hits, ehits ).
first;
1192 bool firstElectronIsOnModuleBorder(
false);
1193 bool secondElectronIsOnModuleBorder(
false);
1199 if( firstElehottestDetId.subdetId() ==
EcalBarrel)
1204 if( firstElehottestDetId.subdetId() ==
EcalBarrel && !firstElectronIsOnModuleBorder )
1211 if( secondElehottestDetId.subdetId() ==
EcalBarrel)
1216 if( secondElehottestDetId.subdetId() ==
EcalBarrel && !secondElectronIsOnModuleBorder )
1223 if ( firstElehottestDetId.subdetId() ==
EcalBarrel && firstElectronIsOnModuleBorder ){
1230 if ( secondElehottestDetId.subdetId() ==
EcalBarrel && secondElectronIsOnModuleBorder ){
1240 if(class1<100 && class2<100){
1242 if(class1==0 && class2==0)
BBZN_gg++;
1243 if(class1<21 && class2<21)
BBZN_tt++;
1244 if(class1<21 || class2<21)
BBZN_t0++;
1248 if(class1>=100 && class2>=100){
1250 if(class1==100 && class2==100)
EEZN_gg++;
1251 if(class1<121 && class2<121)
EEZN_tt++;
1252 if(class1<121 || class2<121)
EEZN_t0++;
1256 if( (class1<100 && class2>=100) || (class2<100 && class1>=100)){
1258 if( (class1==0 && class2==100)||(class2==0 && class1==100) )
EBZN_gg++;
1259 if( ( class1<21 && class2<121) ||(class2<21 && class1<121) )
EBZN_tt++;
1260 if( class2<21 || class1<21 || class2<121 || class1<121 )
EBZN_t0++;
1272 bool selectionBool=
false;
1277 float theta1 = 2. * atan(
exp(- zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->eta()) );
1278 bool ET_1 = ( (zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->energy() *
sin( theta1) ) > 20.);
1280 float theta2 = 2. * atan(
exp(- zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->eta()) );
1281 bool ET_2 = ( (zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->energy() *
sin( theta2) ) > 20.);
1284 bool HoE_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->hadronicOverEm() < 0.115);
1285 bool HoE_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->hadronicOverEm() < 0.115);
1287 bool DeltaPhiIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1288 bool DeltaPhiIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1290 bool DeltaEtaIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1291 bool DeltaEtaIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1295 if(! (invMassBool &&
1298 DeltaPhiIn_1 && DeltaPhiIn_2 &&
1299 DeltaEtaIn_1 && DeltaEtaIn_2
1311 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1312 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1313 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 40 &&
1314 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 140);
1319 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==0 ||
1320 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==10 ||
1321 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==20 ||
1322 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==100 ||
1323 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==110 ||
1324 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==120
1326 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() == 0 ||
1327 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 10 ||
1328 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 20 ||
1329 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 100 ||
1330 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 110 ||
1331 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 120
1337 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() == 0 ||
1338 zeeCandidates[myBestZ].first->getRecoElectron()->classification() == 100
1340 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() == 0 ||
1341 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 100
1347 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() >=30 &&
1348 zeeCandidates[myBestZ].first->getRecoElectron()->classification() <=34)
1350 ((zeeCandidates[myBestZ].
first->getRecoElectron()->classification() >=130 &&
1351 zeeCandidates[myBestZ].first->getRecoElectron()->classification() <=134))
1354 ( (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() >=30 &&
1355 zeeCandidates[myBestZ].second->getRecoElectron()->classification() <=34)
1357 ((zeeCandidates[myBestZ].
second->getRecoElectron()->classification() >=130 &&
1358 zeeCandidates[myBestZ].second->getRecoElectron()->classification() <=134))
1369 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==0 ||
1370 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==10 ||
1371 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==20
1372 ) && zeeCandidates[myBestZ].
second->getRecoElectron()->classification()>=100
1373 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()!=140
1379 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() ==0 ||
1380 zeeCandidates[myBestZ].second->getRecoElectron()->classification() ==10 ||
1381 zeeCandidates[myBestZ].second->getRecoElectron()->classification() ==20
1382 ) && zeeCandidates[myBestZ].
first->getRecoElectron()->classification()>=100
1383 && zeeCandidates[myBestZ].first->getRecoElectron()->classification()!=140
1393 zeeCandidates[myBestZ].first->getRecoElectron()->classification()>=100 &&
1394 zeeCandidates[myBestZ].second->getRecoElectron()->classification()>= 100 &&
1395 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 140 &&
1396 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 140);
1401 zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 &&
1402 zeeCandidates[myBestZ].second->getRecoElectron()->classification()< 100 &&
1403 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1404 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 40);
1409 !(zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 &&
1410 zeeCandidates[myBestZ].second->getRecoElectron()->classification()>=100) &&
1411 !(zeeCandidates[myBestZ].
first->getRecoElectron()->classification()>=100 &&
1412 zeeCandidates[myBestZ].second->getRecoElectron()->classification()<100) );
1415 float ele1EnergyCorrection(1.);
1416 float ele2EnergyCorrection(1.);
1425 if (invMassBool && selectionBool)
1438 std::map<HepMC::GenParticle*,const reco::GsfElectron*> myMCmap;
1440 std::vector<const reco::GsfElectron*> dauElectronCollection;
1442 dauElectronCollection.push_back(zeeCandidates[myBestZ].
first->getRecoElectron() );
1443 dauElectronCollection.push_back(zeeCandidates[myBestZ].
second->getRecoElectron() );
1445 fillMCmap(&dauElectronCollection,mcEle,myMCmap);
1457 if(zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()<100 )
1461 if(zeeCandidates[myBestZ].first->getRecoElectron()->classification()>=100 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()>=100 )
1478 std::cout <<
"Added event to algorithm" << std::endl;
1490 std::cout<<
"[ZeeCalibration] Starting loop number " << iLoop<<std::endl;
1499 std::cout<<
"[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1529 std::cout<<
"[ZeeCalibration] Ending loop " << iLoop<<std::endl;
1541 std::cout<<
"Optimized coefficients " << optimizedCoefficients.size() <<std::endl;
1547 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1561 std::cout<<
"size "<<optimizedCoefficients.size()<<std::endl;
1563 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1565 calibCoeff[ieta] *= optimizedCoefficients[ieta];
1571 std::cout<< ieta <<
" " << optimizedCoefficients[ieta] <<std::endl;
1574 std::vector<DetId> ringIds;
1586 for (
unsigned int iid=0; iid<ringIds.size();++iid){
1589 EBDetId myEBDetId(ringIds[iid]);
1595 EEDetId myEEDetId(ringIds[iid]);
1596 if(myEEDetId.
zside() < 0)
1599 if(myEEDetId.
zside() > 0)
1605 ical->setValue( ringIds[iid], *(
ical->getMap().find(ringIds[iid]) ) * optimizedCoefficients[ieta] );
1642 double parResidual[3];
1643 double errparResidual[3];
1655 std::cout<<
"Fit on residuals, sigma is "<<parResidual[2]<<
" +/- "<<errparResidual[2]<<std::endl;
1694 h1_seedOverSC_=
new TH1F(
"h1_seedOverSC",
"h1_seedOverSC", 400, 0., 2.);
1701 h2_fEtaBarrelGood_ =
new TH2F(
"fEtaBarrelGood",
"fEtaBarrelGood",800,-4.,4.,800,0.8,1.2);
1705 h2_fEtaBarrelBad_ =
new TH2F(
"fEtaBarrelBad",
"fEtaBarrelBad",800,-4.,4.,800,0.8,1.2);
1709 h2_fEtaEndcapGood_ =
new TH2F(
"fEtaEndcapGood",
"fEtaEndcapGood",800,-4.,4.,800,0.8,1.2);
1713 h2_fEtaEndcapBad_ =
new TH2F(
"fEtaEndcapBad",
"fEtaEndcapBad",800,-4.,4.,800,0.8,1.2);
1717 for (
int i=0;
i<2;
i++)
1721 sprintf(histoName,
"h_eleEffEta_%d",
i);
1722 h_eleEffEta_[
i] =
new TH1F(histoName,histoName, 150, 0., 2.7);
1725 sprintf(histoName,
"h_eleEffPhi_%d",
i);
1726 h_eleEffPhi_[
i] =
new TH1F(histoName,histoName, 400, -4., 4.);
1729 sprintf(histoName,
"h_eleEffPt_%d",
i);
1730 h_eleEffPt_[
i] =
new TH1F(histoName,histoName, 200, 0., 200.);
1744 for (
int i=0;
i<25;
i++)
1748 sprintf(histoName,
"h_ESCEtrueVsEta_%d",
i);
1754 sprintf(histoName,
"h_ESCEtrue_%d",
i);
1756 h_ESCEtrue_[
i] =
new TH1F(histoName,histoName, 300,0.,1.5);
1758 sprintf(histoName,
"h2_chi2_%d",
i);
1759 h2_chi2_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 5);
1761 sprintf(histoName,
"h2_iterations_%d",
i);
1762 h2_iterations_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 15);
1764 sprintf(histoName,
"h_ESCcorrEtrueVsEta_%d",
i);
1770 sprintf(histoName,
"h_ESCcorrEtrue_%d",
i);
1774 sprintf(histoName,
"h2_xtalRecalibCoeffBarrel_%d",
i);
1780 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapMinus_%d",
i);
1785 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapPlus_%d",
i);
1811 h1_zMassResol_ =
new TH1F(
"zMassResol",
"zMassResol", 200, -50., 50.);
1815 h1_eleEtaResol_ =
new TH1F(
"eleEtaResol",
"eleEtaResol", 100, -0.01, 0.01);
1831 h1_elePhiResol_ =
new TH1F(
"elePhiResol",
"elePhiResol", 100, -0.01, 0.01);
1836 h1_zEtaResol_ =
new TH1F(
"zEtaResol",
"zEtaResol", 200, -1., 1.);
1841 h1_zPhiResol_ =
new TH1F(
"zPhiResol",
"zPhiResol", 200, -1., 1.);
1845 h1_nEleReco_ =
new TH1F(
"nEleReco",
"Number of reco electrons",10,-0.5,10.5);
1858 h1_occupancy_ =
new TH1F(
"occupancy",
"occupancy",1000,0,10000);
1883 h1_ZCandMult_ =
new TH1F(
"ZCandMult",
"Multiplicity of Z candidates in one event",10,-0.5,10.5);
1886 h1_reco_ZMass_ =
new TH1F(
"reco_ZMass",
"Inv. mass of 2 reco Electrons",200,0.,150.);
1890 h1_reco_ZMassCorr_ =
new TH1F(
"reco_ZMassCorr",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1894 h1_reco_ZMassCorrBB_ =
new TH1F(
"reco_ZMassCorrBB",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1898 h1_reco_ZMassCorrEE_ =
new TH1F(
"reco_ZMassCorrEE",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1904 h2_coeffVsEta_=
new TH2F(
"h2_calibCoeffVsEta",
"h2_calibCoeffVsEta",249,-124,125, 200, 0., 2.);
1908 h2_coeffVsEtaGrouped_=
new TH2F(
"h2_calibCoeffVsEtaGrouped",
"h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1912 h2_zMassVsLoop_=
new TH2F(
"h2_zMassVsLoop",
"h2_zMassVsLoop",1000,0,40, 90, 80.,95.);
1914 h2_zMassDiffVsLoop_=
new TH2F(
"h2_zMassDiffVsLoop",
"h2_zMassDiffVsLoop",1000,0,40, 100, -1., 1.);
1918 h2_zWidthVsLoop_=
new TH2F(
"h2_zWidthVsLoop",
"h2_zWidthVsLoop",1000,0,40, 100, 0.,10.);
1920 h2_coeffVsLoop_=
new TH2F(
"h2_coeffVsLoop",
"h2_coeffVsLoop",1000,0,40, 100, 0., 2.);
1922 h2_residualSigma_=
new TH2F(
"h2_residualSigma",
"h2_residualSigma",1000, 0, 40, 100, 0., .5);
1924 h2_miscalRecal_ =
new TH2F(
"h2_miscalRecal",
"h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1928 h2_miscalRecalEB_ =
new TH2F(
"h2_miscalRecalEB",
"h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1932 h2_miscalRecalEE_ =
new TH2F(
"h2_miscalRecalEE",
"h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1936 h1_mc_ =
new TH1F(
"h1_residualMiscalib",
"h1_residualMiscalib", 200, -0.2, 0.2);
1937 h1_mcEB_ =
new TH1F(
"h1_residualMiscalibEB",
"h1_residualMiscalibEB", 200, -0.2, 0.2);
1938 h1_mcEE_ =
new TH1F(
"h1_residualMiscalibEE",
"h1_residualMiscalibEE", 200, -0.2, 0.2);
1940 for (
int i=0;
i<25;
i++)
1960 sprintf(histoName,
"h1_residualMiscalibParz_%d",
i);
1961 h1_mcParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1962 sprintf(histoName,
"h1_residualMiscalibEBParz_%d",
i);
1963 h1_mcEBParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1964 sprintf(histoName,
"h1_residualMiscalibEEParz_%d",
i);
1965 h1_mcEEParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1975 float p0 = 1.00153e+00;
1976 float p1 = 3.29331e-02;
1977 float p2 = 1.21187e-03;
1979 double x = (double) fabs(scEta);
1981 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
1989 float p0 = 1.06819e+00;
1990 float p1 = -1.53189e-02;
1991 float p2 = 4.01707e-04 ;
1993 double x = (double) fabs(scEta);
1995 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2001 float p0 = 1.17382e+00;
2002 float p1 = -6.52319e-02;
2003 float p2 = 6.26108e-03;
2005 double x = (double) fabs(scEta);
2007 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2013 float p0 = 9.99782e-01 ;
2014 float p1 = 1.26983e-02;
2015 float p2 = 2.16344e-03;
2017 double x = (double) fabs(scEta);
2019 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2026 void ZeeCalibration::fillMCmap(
const std::vector<const reco::GsfElectron*>* electronCollection,
const std::vector<HepMC::GenParticle*>& mcEle,std::map<HepMC::GenParticle*,const reco::GsfElectron*>& myMCmap)
2028 for (
unsigned int i=0;
i<mcEle.size();
i++)
2032 for (
unsigned int j=0;
j<electronCollection->size();
j++)
2034 float dr=
EvalDR(mcEle[
i]->momentum().pseudoRapidity(),(*(*electronCollection)[
j]).
eta(),mcEle[
i]->momentum().
phi(),(*(*electronCollection)[j]).
phi());
2037 myMatchEle = (*electronCollection)[
j];
2041 myMCmap.insert(std::pair<HepMC::GenParticle*,const reco::GsfElectron*>(mcEle[
i],myMatchEle));
2049 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2050 float DPhi = Phi - Phi_ref ;
2053 float DEta = Eta - Eta_ref ;
2055 float DR =
sqrt( DEta*DEta + DPhi*DPhi );
2062 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2063 return (Phi - Phi_ref);
2069 for (
unsigned int i=0;
i<mcEle.size();
i++)
2072 h_eleEffEta_[0]->Fill(fabs(mcEle[
i]->momentum().pseudoRapidity()));
2076 std::map<HepMC::GenParticle*,const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[
i]);
2077 if (mIter == associationMap.end() )
2084 h_eleEffEta_[1]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
2119 if(k>=0 && k<=84)index = k - 85;
2121 if(k>=85 && k<=169)index = k - 84;
2123 if(k>=170 && k<=208)index = - k + 84;
2125 if(k>=209 && k<=247)index = k - 123;
2131 if(k>=0 && k<=71)index = k - 72;
2133 if(k>=72 && k<=143)index = k - 71;
2143 double correction(1.);
2176 double maxEnergy = -9999.;
2179 std::pair<DetId, double> myPair (
DetId(0), -9999.);
2182 for( std::vector<std::pair<DetId,float> >::const_iterator idIt=mySCRecHits.begin(); idIt != mySCRecHits.end(); idIt++){
2186 hottestRecHit = & (* ( ebhits->
find((*idIt).first) ) );
2188 if( hottestRecHit == & (*( ebhits->
end())) )
2190 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2194 else if (idIt->first.subdetId() ==
EcalEndcap )
2196 hottestRecHit = & (* ( eehits->
find((*idIt).first) ) );
2197 if( hottestRecHit == & (*( eehits->
end())) )
2199 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2206 if(hottestRecHit && hottestRecHit->
energy() > maxEnergy){
2208 maxEnergy = hottestRecHit->
energy();
2210 myPair.first = hottestRecHit ->
id();
2211 myPair.second = maxEnergy;
2228 short ieta = myEBDetId.
ieta();
2229 short iphi = myEBDetId.
iphi();
2233 myBool = (
abs( ieta ) == 1 ||
abs( ieta ) == 25
2234 ||
abs( ieta ) ==26 ||
abs( ieta ) == 45
2235 ||
abs( ieta ) ==46 ||
abs( ieta ) == 65
2236 ||
abs( ieta ) ==66 ||
abs( ieta ) == 85 );
2238 for(
int i = 0;
i < 19;
i++){
2240 if(iphi == ( 20*
i + 1 ) || iphi == 20*
i )
2256 bool isNearCrack =
false;
2268 dist +=
pow( v1[
i]-v2[
i], 2 );
2324 for (
int i=0;
i<2;
i++)
2382 std::ofstream
fout(
"ZeeStatistics.txt");
2385 std::cout <<
"Cannot open output file.\n";
2388 fout<<
"ZeeStatistics"<<std::endl;
2390 fout<<
"##########################RECO#########################"<<std::endl;
2391 fout<<
"##################Zee with Barrel-Barrel electrons: "<<
BBZN<<std::endl;
2394 fout<<
"##################Zee with Barrel-Endcap electrons: "<<
EBZN<<std::endl;
2396 fout<<
"##################Zee with Endcap-Endcap electrons: "<<
EEZN<<std::endl;
2399 fout<<
"\n"<<std::endl;
2401 fout<<
"##########################GEN#########################"<<std::endl;
2402 fout<<
"##################Zee with Barrel-Barrel electrons: "<<(float)
MCZBB/
NEVT<<std::endl;
2403 fout<<
"##################Zee with Barrel-Endcap electrons: "<<(float)
MCZEB/
NEVT<<std::endl;
2404 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)
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_
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]
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]
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
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_