18 #include "CLHEP/Vector/LorentzVector.h"
66 #include "TGraphErrors.h"
85 std::cout<<
"[ZeeCalibration] Starting the ctor"<<std::endl;
117 outputFile_ = TFile::Open(outputFileName_.c_str(),
"RECREATE");
119 myTree =
new TTree(
"myTree",
"myTree");
121 myTree->Branch(
"zMass",&
mass4tree,
"mass/F");
160 findingRecord<EcalIntercalibConstantsRcd> () ;
162 for(
int i = 0;
i<50;
i++){
172 std::cout<<
"[ZeeCalibration] Done with the ctor"<<std::endl;
186 boost::shared_ptr<EcalIntercalibConstants>
189 std::cout <<
"@SUB=ZeeCalibration::produceEcalIntercalibConstants" << std::endl;
212 if(ieta==0)
continue;
246 std::cout<<
"Writing histos..."<<std::endl;
275 for(
unsigned int i =0;
i<25;
i++){
312 Double_t
mean[25] = {0.};
313 Double_t
num[25] = {0.};
314 Double_t meanErr[25] = {0.};
315 Float_t
rms[25] = {0.};
316 Float_t tempRms[10][25];
318 for(
int ia = 0; ia<10; ia++){
319 for(
int ib = 0; ib<25; ib++){
321 tempRms[ia][ib] = 0.;
334 bool isNearCrack =
false;
399 if(
k>=170 &&
k<=204){
439 if(
k>=205 &&
k<=208){
487 for(
int ic = 0; ic< 17; ic++){
489 mean[ic] = mean[ic] / num[ic];
491 meanErr[ic] = 1. /
sqrt(meanErr[ic]);
497 for(
int ic = 0; ic< 25; ic++){
498 for(
int id = 0;
id< 10;
id++){
500 if(tempRms[
id][ic] > 0.){
502 rms[ic] += (tempRms[
id][ic] - mean[
j])*(tempRms[
id][ic] - mean[j]);
507 rms[ic] =
sqrt(rms[ic]);
514 Double_t xtalEta[25] = {1.4425, 1.3567,1.2711,1.1855,
516 0.7468,0.6612,0.5756,0.4897,0.3985,0.3117,0.2250,0.1384,0.0487,
517 1.546, 1.651, 1.771, 1.908, 2.071, 2.267, 2.516, 2.8};
519 Double_t
zero[25] = {0.026};
521 for(
int j = 0; j <25; j++)
528 px->SetXTitle(
"Eta channel");
529 px->SetYTitle(
"recalibCoeff");
549 double weightSumMeanBarrel = 0.;
550 double weightSumMeanEndcap = 0.;
588 std::cout<<
"Weight sum mean on channels in Barrel is :"<<weightSumMeanBarrel<<std::endl;
589 std::cout<<
"Weight sum mean on channels in Endcap is :"<<weightSumMeanEndcap<<std::endl;
601 TGraphErrors*
graph =
new TGraphErrors(25,xtalEta,mean,zero,meanErr);
605 double zero50[50] = { 0. };
608 residualSigmaGraph->SetName(
"residualSigmaGraph");
609 residualSigmaGraph->Draw(
"APL");
610 residualSigmaGraph->Write();
613 coefficientDistanceAtIterationGraph->SetName(
"coefficientDistanceAtIterationGraph");
614 coefficientDistanceAtIterationGraph->Draw(
"APL");
615 coefficientDistanceAtIterationGraph->Write();
617 Float_t noError[250] = {0.};
619 Float_t ringInd[250];
620 for(
int i =0;
i<250;
i++)
624 graphCoeff->SetName(
"graphCoeff");
625 graphCoeff->Draw(
"APL");
664 std::cout<<
"[ZeeCalibration] Entering duringLoop"<<std::endl;
683 std::cout<<
"[ZeeCalibration::beginOfJob] Histograms booked "<<std::endl;
701 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EB miscal file"<<std::endl;
710 std::cout<<
"[ZeeCalibration::beginOfJob] Parsed EE miscal file"<<std::endl;
714 std::cout <<
" theAlgorithm_->getNumberOfChannels() "
724 std::vector<DetId> ringIds;
738 for (
unsigned int iid=0; iid<ringIds.size();++iid)
740 float miscalib=* (miscalibMap->
get().
getMap().
find(ringIds[iid]) );
759 std::vector<DetId> ringIds;
771 for (
unsigned int iid=0; iid<ringIds.size();++iid){
775 EBDetId myEBDetId(ringIds[iid]);
781 EEDetId myEEDetId(ringIds[iid]);
782 if(myEEDetId.
zside() < 0)
785 if(myEEDetId.
zside() > 0)
790 ical->setValue( ringIds[iid], *(miscalibMap->
get().
getMap().
find(ringIds[iid]) ) );
813 for(
unsigned int iHLT=0; iHLT<200; ++iHLT) {
818 std::cout<<
"[ZeeCalibration::duringLoop] Done with initializing aHLTresults[] "<<std::endl;
824 if(!hltTriggerResultHandle.isValid()) {
828 hltCount = hltTriggerResultHandle->size();
834 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillHLTInfo(hltTriggerResultHandle); "<<std::endl;
852 std::cout<<
"[ZeeCalibration::duringLoop] End HLT section"<<std::endl;
858 std::vector<HepMC::GenParticle*> mcEle;
860 float myGenZMass(-1);
870 const HepMC::GenEvent * myGenEvent = hepProd->GetEvent();
876 std::cout<<
"[ZeeCalibration::duringLoop] Done with myZeePlots_->fillZMCInfo( & (*myGenEvent) ); "<<std::endl;
882 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
883 p != myGenEvent->particles_end(); ++
p ) {
885 if ( (*p)->pdg_id() == 23 && (*p)->status()==2){
887 myGenZMass = (*p)->momentum().m();
901 for ( HepMC::GenEvent::particle_const_iterator
p = myGenEvent->particles_begin();
902 p != myGenEvent->particles_end(); ++
p ) {
904 if (
abs( (*p)->pdg_id() ) == 11 )
906 mcEle.push_back( (*
p) );
913 if(mcEle.size()==2 && fabs(mcEle[0]->momentum().
eta())<2.4 && fabs(mcEle[1]->momentum().
eta())<2.4 ){
916 if( fabs(mcEle[0]->momentum().
eta())<1.479 && fabs(mcEle[1]->momentum().
eta())<1.479 )
MCZBB++;
918 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++;
920 if( fabs(mcEle[0]->momentum().eta())>1.479 && fabs(mcEle[1]->momentum().eta())>1.479 )
MCZEE++;
934 std::cerr <<
"Error! can't get the product EBRecHitCollection " << std::endl;
943 std::cerr <<
"Error! can't get the product EERecHitCollection " << std::endl;
953 std::cerr <<
"Error! can't get the product SuperClusterCollection "<< std::endl;
970 std::cerr <<
"Error! can't get the product IslandSuperClusterCollection "<< std::endl;
975 std::cout<<
"scCollection->size()"<<scIslandCollection->size()<<std::endl;
978 if( (
scCollection->size()+scIslandCollection->size() ) < 2)
986 std::cerr <<
"Error! can't get the product ElectronCollection "<< std::endl;
1002 if ( !hits && !ehits){
1007 if (hits->size() == 0 && ehits->size() == 0){
1008 std::cout <<
"hits->size() == 0" << std::endl;
1013 std::cout <<
"!electronCollection" << std::endl;
1018 std::cout <<
"electronCollection->size() == 0" << std::endl;
1035 std::cout <<
" Starting with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1042 std::cout <<
" Done with myZeePlots_->fillEleInfo(electronCollection); " << std::endl;
1048 std::vector<calib::CalibElectron> calibElectrons;
1058 std::cout << calibElectrons.back().getRecoElectron()->superCluster()->energy() <<
" " << calibElectrons.back().getRecoElectron()->energy() << std::endl;
1066 std::cout <<
"Filled histos" << std::endl;
1070 std::vector<std::pair<calib::CalibElectron*,calib::CalibElectron*> > zeeCandidates;
1074 double DeltaMinvMin(5000.);
1076 if (calibElectrons.size() < 2)
1079 for(
unsigned int e_it = 0 ; e_it != calibElectrons.size() - 1 ; e_it++){
1080 for(
unsigned int p_it = e_it + 1 ; p_it != calibElectrons.size() ; p_it++)
1083 std::cout << e_it <<
" " << calibElectrons[e_it].getRecoElectron()->charge() <<
" " << p_it <<
" " << calibElectrons[p_it].getRecoElectron()->charge() << std::endl;
1085 if (calibElectrons[e_it].getRecoElectron()->
charge() * calibElectrons[p_it].getRecoElectron()->
charge() != -1)
1094 std::cout <<
"#######################mass "<<
mass << std::endl;
1097 zeeCandidates.push_back(std::pair<calib::CalibElectron*,calib::CalibElectron*>(&(calibElectrons[e_it]),&(calibElectrons[p_it])));
1098 double DeltaMinv = fabs(
mass -
MZ);
1100 if( DeltaMinv < DeltaMinvMin)
1102 DeltaMinvMin = DeltaMinv;
1103 myBestZ=zeeCandidates.size()-1;
1114 if(zeeCandidates.size()==0 || myBestZ==-1 )
1121 std::cout <<
"Found ZCandidates " << myBestZ << std::endl;
1132 int class1 = zeeCandidates[myBestZ].first->getRecoElectron()->classification();
1133 int class2 = zeeCandidates[myBestZ].second->getRecoElectron()->classification();
1160 if( class1==10 || class1 ==20)
1162 if( class1==110 || class1 ==120)
1164 if( class1>=30 && class1 <=34)
1166 if( class1>=130 && class1 <=134)
1177 if( class2==10 || class2 ==20)
1179 if( class2==110 || class2 ==120)
1181 if( class2>=30 && class2 <=34)
1183 if( class2>=130 && class2 <=134)
1195 DetId firstElehottestDetId =
getHottestDetId( zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->seed()->hitsAndFractions() , hits, ehits ).
first;
1196 DetId secondElehottestDetId =
getHottestDetId( zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->seed()->hitsAndFractions() , hits, ehits ).
first;
1198 bool firstElectronIsOnModuleBorder(
false);
1199 bool secondElectronIsOnModuleBorder(
false);
1205 if( firstElehottestDetId.subdetId() ==
EcalBarrel)
1210 if( firstElehottestDetId.subdetId() ==
EcalBarrel && !firstElectronIsOnModuleBorder )
1217 if( secondElehottestDetId.subdetId() ==
EcalBarrel)
1222 if( secondElehottestDetId.subdetId() ==
EcalBarrel && !secondElectronIsOnModuleBorder )
1229 if ( firstElehottestDetId.subdetId() ==
EcalBarrel && firstElectronIsOnModuleBorder ){
1236 if ( secondElehottestDetId.subdetId() ==
EcalBarrel && secondElectronIsOnModuleBorder ){
1246 if(class1<100 && class2<100){
1248 if(class1==0 && class2==0)
BBZN_gg++;
1249 if(class1<21 && class2<21)
BBZN_tt++;
1250 if(class1<21 || class2<21)
BBZN_t0++;
1254 if(class1>=100 && class2>=100){
1256 if(class1==100 && class2==100)
EEZN_gg++;
1257 if(class1<121 && class2<121)
EEZN_tt++;
1258 if(class1<121 || class2<121)
EEZN_t0++;
1262 if( (class1<100 && class2>=100) || (class2<100 && class1>=100)){
1264 if( (class1==0 && class2==100)||(class2==0 && class1==100) )
EBZN_gg++;
1265 if( ( class1<21 && class2<121) ||(class2<21 && class1<121) )
EBZN_tt++;
1266 if( class2<21 || class1<21 || class2<121 || class1<121 )
EBZN_t0++;
1278 bool selectionBool=
false;
1283 float theta1 = 2. * atan(
exp(- zeeCandidates[myBestZ].
first->getRecoElectron()->superCluster()->eta()) );
1284 bool ET_1 = ( (zeeCandidates[myBestZ].first->getRecoElectron()->superCluster()->energy() *
sin( theta1) ) > 20.);
1286 float theta2 = 2. * atan(
exp(- zeeCandidates[myBestZ].
second->getRecoElectron()->superCluster()->eta()) );
1287 bool ET_2 = ( (zeeCandidates[myBestZ].second->getRecoElectron()->superCluster()->energy() *
sin( theta2) ) > 20.);
1290 bool HoE_1 = (zeeCandidates[myBestZ].first->getRecoElectron()->hadronicOverEm() < 0.115);
1291 bool HoE_2 = (zeeCandidates[myBestZ].second->getRecoElectron()->hadronicOverEm() < 0.115);
1293 bool DeltaPhiIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1294 bool DeltaPhiIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaPhiSuperClusterTrackAtVtx() < 0.090);
1296 bool DeltaEtaIn_1 = ( zeeCandidates[myBestZ].first->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1297 bool DeltaEtaIn_2 = ( zeeCandidates[myBestZ].second->getRecoElectron()->deltaEtaSuperClusterTrackAtVtx() < 0.0090);
1301 if(! (invMassBool &&
1304 DeltaPhiIn_1 && DeltaPhiIn_2 &&
1305 DeltaEtaIn_1 && DeltaEtaIn_2
1317 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1318 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1319 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 40 &&
1320 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 140);
1325 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==0 ||
1326 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==10 ||
1327 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==20 ||
1328 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==100 ||
1329 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==110 ||
1330 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==120
1332 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() == 0 ||
1333 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 10 ||
1334 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 20 ||
1335 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 100 ||
1336 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 110 ||
1337 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 120
1343 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() == 0 ||
1344 zeeCandidates[myBestZ].first->getRecoElectron()->classification() == 100
1346 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() == 0 ||
1347 zeeCandidates[myBestZ].second->getRecoElectron()->classification() == 100
1353 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() >=30 &&
1354 zeeCandidates[myBestZ].first->getRecoElectron()->classification() <=34)
1356 ((zeeCandidates[myBestZ].
first->getRecoElectron()->classification() >=130 &&
1357 zeeCandidates[myBestZ].first->getRecoElectron()->classification() <=134))
1360 ( (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() >=30 &&
1361 zeeCandidates[myBestZ].second->getRecoElectron()->classification() <=34)
1363 ((zeeCandidates[myBestZ].
second->getRecoElectron()->classification() >=130 &&
1364 zeeCandidates[myBestZ].second->getRecoElectron()->classification() <=134))
1375 (zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==0 ||
1376 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==10 ||
1377 zeeCandidates[myBestZ].first->getRecoElectron()->classification() ==20
1378 ) && zeeCandidates[myBestZ].
second->getRecoElectron()->classification()>=100
1379 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()!=140
1385 (zeeCandidates[myBestZ].
second->getRecoElectron()->classification() ==0 ||
1386 zeeCandidates[myBestZ].second->getRecoElectron()->classification() ==10 ||
1387 zeeCandidates[myBestZ].second->getRecoElectron()->classification() ==20
1388 ) && zeeCandidates[myBestZ].
first->getRecoElectron()->classification()>=100
1389 && zeeCandidates[myBestZ].first->getRecoElectron()->classification()!=140
1399 zeeCandidates[myBestZ].first->getRecoElectron()->classification()>=100 &&
1400 zeeCandidates[myBestZ].second->getRecoElectron()->classification()>= 100 &&
1401 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 140 &&
1402 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 140);
1407 zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 &&
1408 zeeCandidates[myBestZ].second->getRecoElectron()->classification()< 100 &&
1409 zeeCandidates[myBestZ].first->getRecoElectron()->classification()!= 40 &&
1410 zeeCandidates[myBestZ].second->getRecoElectron()->classification()!= 40);
1415 !(zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 &&
1416 zeeCandidates[myBestZ].second->getRecoElectron()->classification()>=100) &&
1417 !(zeeCandidates[myBestZ].
first->getRecoElectron()->classification()>=100 &&
1418 zeeCandidates[myBestZ].second->getRecoElectron()->classification()<100) );
1421 float ele1EnergyCorrection(1.);
1422 float ele2EnergyCorrection(1.);
1431 if (invMassBool && selectionBool)
1444 std::map<HepMC::GenParticle*,const reco::GsfElectron*> myMCmap;
1446 std::vector<const reco::GsfElectron*> dauElectronCollection;
1448 dauElectronCollection.push_back(zeeCandidates[myBestZ].
first->getRecoElectron() );
1449 dauElectronCollection.push_back(zeeCandidates[myBestZ].
second->getRecoElectron() );
1451 fillMCmap(&dauElectronCollection,mcEle,myMCmap);
1463 if(zeeCandidates[myBestZ].first->getRecoElectron()->classification()<100 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()<100 )
1467 if(zeeCandidates[myBestZ].first->getRecoElectron()->classification()>=100 && zeeCandidates[myBestZ].second->getRecoElectron()->classification()>=100 )
1484 std::cout <<
"Added event to algorithm" << std::endl;
1496 std::cout<<
"[ZeeCalibration] Starting loop number " << iLoop<<std::endl;
1505 std::cout<<
"[ZeeCalibration] exiting from startingNewLoop" << std::endl;
1535 std::cout<<
"[ZeeCalibration] Ending loop " << iLoop<<std::endl;
1547 std::cout<<
"Optimized coefficients " << optimizedCoefficients.size() <<std::endl;
1553 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1567 std::cout<<
"size "<<optimizedCoefficients.size()<<std::endl;
1569 for (
unsigned int ieta=0;ieta<optimizedCoefficients.size();ieta++)
1571 calibCoeff[ieta] *= optimizedCoefficients[ieta];
1577 std::cout<< ieta <<
" " << optimizedCoefficients[ieta] <<std::endl;
1580 std::vector<DetId> ringIds;
1592 for (
unsigned int iid=0; iid<ringIds.size();++iid){
1595 EBDetId myEBDetId(ringIds[iid]);
1601 EEDetId myEEDetId(ringIds[iid]);
1602 if(myEEDetId.
zside() < 0)
1605 if(myEEDetId.
zside() > 0)
1611 ical->setValue( ringIds[iid], *(
ical->getMap().find(ringIds[iid]) ) * optimizedCoefficients[ieta] );
1648 double parResidual[3];
1649 double errparResidual[3];
1661 std::cout<<
"Fit on residuals, sigma is "<<parResidual[2]<<
" +/- "<<errparResidual[2]<<std::endl;
1700 h1_seedOverSC_=
new TH1F(
"h1_seedOverSC",
"h1_seedOverSC", 400, 0., 2.);
1707 h2_fEtaBarrelGood_ =
new TH2F(
"fEtaBarrelGood",
"fEtaBarrelGood",800,-4.,4.,800,0.8,1.2);
1711 h2_fEtaBarrelBad_ =
new TH2F(
"fEtaBarrelBad",
"fEtaBarrelBad",800,-4.,4.,800,0.8,1.2);
1715 h2_fEtaEndcapGood_ =
new TH2F(
"fEtaEndcapGood",
"fEtaEndcapGood",800,-4.,4.,800,0.8,1.2);
1719 h2_fEtaEndcapBad_ =
new TH2F(
"fEtaEndcapBad",
"fEtaEndcapBad",800,-4.,4.,800,0.8,1.2);
1723 for (
int i=0;
i<2;
i++)
1727 sprintf(histoName,
"h_eleEffEta_%d",
i);
1728 h_eleEffEta_[
i] =
new TH1F(histoName,histoName, 150, 0., 2.7);
1731 sprintf(histoName,
"h_eleEffPhi_%d",
i);
1732 h_eleEffPhi_[
i] =
new TH1F(histoName,histoName, 400, -4., 4.);
1735 sprintf(histoName,
"h_eleEffPt_%d",
i);
1736 h_eleEffPt_[
i] =
new TH1F(histoName,histoName, 200, 0., 200.);
1750 for (
int i=0;
i<25;
i++)
1754 sprintf(histoName,
"h_ESCEtrueVsEta_%d",
i);
1760 sprintf(histoName,
"h_ESCEtrue_%d",
i);
1762 h_ESCEtrue_[
i] =
new TH1F(histoName,histoName, 300,0.,1.5);
1764 sprintf(histoName,
"h2_chi2_%d",
i);
1765 h2_chi2_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 5);
1767 sprintf(histoName,
"h2_iterations_%d",
i);
1768 h2_iterations_[
i] =
new TH2F(histoName,histoName, 1000,-150,150, 1000, -1, 15);
1770 sprintf(histoName,
"h_ESCcorrEtrueVsEta_%d",
i);
1776 sprintf(histoName,
"h_ESCcorrEtrue_%d",
i);
1780 sprintf(histoName,
"h2_xtalRecalibCoeffBarrel_%d",
i);
1786 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapMinus_%d",
i);
1791 sprintf(histoName,
"h2_xtalRecalibCoeffEndcapPlus_%d",
i);
1817 h1_zMassResol_ =
new TH1F(
"zMassResol",
"zMassResol", 200, -50., 50.);
1821 h1_eleEtaResol_ =
new TH1F(
"eleEtaResol",
"eleEtaResol", 100, -0.01, 0.01);
1837 h1_elePhiResol_ =
new TH1F(
"elePhiResol",
"elePhiResol", 100, -0.01, 0.01);
1842 h1_zEtaResol_ =
new TH1F(
"zEtaResol",
"zEtaResol", 200, -1., 1.);
1847 h1_zPhiResol_ =
new TH1F(
"zPhiResol",
"zPhiResol", 200, -1., 1.);
1851 h1_nEleReco_ =
new TH1F(
"nEleReco",
"Number of reco electrons",10,-0.5,10.5);
1864 h1_occupancy_ =
new TH1F(
"occupancy",
"occupancy",1000,0,10000);
1889 h1_ZCandMult_ =
new TH1F(
"ZCandMult",
"Multiplicity of Z candidates in one event",10,-0.5,10.5);
1892 h1_reco_ZMass_ =
new TH1F(
"reco_ZMass",
"Inv. mass of 2 reco Electrons",200,0.,150.);
1896 h1_reco_ZMassCorr_ =
new TH1F(
"reco_ZMassCorr",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1900 h1_reco_ZMassCorrBB_ =
new TH1F(
"reco_ZMassCorrBB",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1904 h1_reco_ZMassCorrEE_ =
new TH1F(
"reco_ZMassCorrEE",
"Inv. mass of 2 corrected reco Electrons",200,0.,150.);
1910 h2_coeffVsEta_=
new TH2F(
"h2_calibCoeffVsEta",
"h2_calibCoeffVsEta",249,-124,125, 200, 0., 2.);
1914 h2_coeffVsEtaGrouped_=
new TH2F(
"h2_calibCoeffVsEtaGrouped",
"h2_calibCoeffVsEtaGrouped", 200, 0., 3., 200, 0.6, 1.4);
1918 h2_zMassVsLoop_=
new TH2F(
"h2_zMassVsLoop",
"h2_zMassVsLoop",1000,0,40, 90, 80.,95.);
1920 h2_zMassDiffVsLoop_=
new TH2F(
"h2_zMassDiffVsLoop",
"h2_zMassDiffVsLoop",1000,0,40, 100, -1., 1.);
1924 h2_zWidthVsLoop_=
new TH2F(
"h2_zWidthVsLoop",
"h2_zWidthVsLoop",1000,0,40, 100, 0.,10.);
1926 h2_coeffVsLoop_=
new TH2F(
"h2_coeffVsLoop",
"h2_coeffVsLoop",1000,0,40, 100, 0., 2.);
1928 h2_residualSigma_=
new TH2F(
"h2_residualSigma",
"h2_residualSigma",1000, 0, 40, 100, 0., .5);
1930 h2_miscalRecal_ =
new TH2F(
"h2_miscalRecal",
"h2_miscalRecal", 500, 0., 2., 500, 0., 2.);
1934 h2_miscalRecalEB_ =
new TH2F(
"h2_miscalRecalEB",
"h2_miscalRecalEB", 500, 0., 2., 500, 0., 2.);
1938 h2_miscalRecalEE_ =
new TH2F(
"h2_miscalRecalEE",
"h2_miscalRecalEE", 500, 0., 2., 500, 0., 2.);
1942 h1_mc_ =
new TH1F(
"h1_residualMiscalib",
"h1_residualMiscalib", 200, -0.2, 0.2);
1943 h1_mcEB_ =
new TH1F(
"h1_residualMiscalibEB",
"h1_residualMiscalibEB", 200, -0.2, 0.2);
1944 h1_mcEE_ =
new TH1F(
"h1_residualMiscalibEE",
"h1_residualMiscalibEE", 200, -0.2, 0.2);
1946 for (
int i=0;
i<25;
i++)
1966 sprintf(histoName,
"h1_residualMiscalibParz_%d",
i);
1967 h1_mcParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1968 sprintf(histoName,
"h1_residualMiscalibEBParz_%d",
i);
1969 h1_mcEBParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1970 sprintf(histoName,
"h1_residualMiscalibEEParz_%d",
i);
1971 h1_mcEEParz_[
i] =
new TH1F(histoName,histoName, 200, -0.2, 0.2);
1981 float p0 = 1.00153e+00;
1982 float p1 = 3.29331e-02;
1983 float p2 = 1.21187e-03;
1985 double x = (double) fabs(scEta);
1987 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
1995 float p0 = 1.06819e+00;
1996 float p1 = -1.53189e-02;
1997 float p2 = 4.01707e-04 ;
1999 double x = (double) fabs(scEta);
2001 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2007 float p0 = 1.17382e+00;
2008 float p1 = -6.52319e-02;
2009 float p2 = 6.26108e-03;
2011 double x = (double) fabs(scEta);
2013 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2019 float p0 = 9.99782e-01 ;
2020 float p1 = 1.26983e-02;
2021 float p2 = 2.16344e-03;
2023 double x = (double) fabs(scEta);
2025 return 1. / ( p0 + p1*x*x + p2*x*x*x*
x );
2034 for (
unsigned int i=0;
i<mcEle.size();
i++)
2038 for (
unsigned int j=0;
j<electronCollection->size();
j++)
2040 float dr=
EvalDR(mcEle[
i]->momentum().pseudoRapidity(),(*(*electronCollection)[
j]).
eta(),mcEle[
i]->momentum().
phi(),(*(*electronCollection)[j]).
phi());
2043 myMatchEle = (*electronCollection)[
j];
2047 myMCmap.insert(std::pair<HepMC::GenParticle*,const reco::GsfElectron*>(mcEle[
i],myMatchEle));
2055 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2056 float DPhi = Phi - Phi_ref ;
2059 float DEta = Eta - Eta_ref ;
2061 float DR =
sqrt( DEta*DEta + DPhi*DPhi );
2068 if (Phi_ref<0) Phi_ref = 2*
TMath::Pi() + Phi_ref;
2069 return (Phi - Phi_ref);
2075 for (
unsigned int i=0;
i<mcEle.size();
i++)
2078 h_eleEffEta_[0]->Fill(fabs(mcEle[
i]->momentum().pseudoRapidity()));
2082 std::map<HepMC::GenParticle*,const reco::GsfElectron*>::const_iterator mIter = associationMap.find(mcEle[
i]);
2083 if (mIter == associationMap.end() )
2090 h_eleEffEta_[1]->Fill(fabs(mcEle[i]->momentum().pseudoRapidity()));
2125 if(k>=0 && k<=84)index = k - 85;
2127 if(k>=85 && k<=169)index = k - 84;
2129 if(k>=170 && k<=208)index = - k + 84;
2131 if(k>=209 && k<=247)index = k - 123;
2137 if(k>=0 && k<=71)index = k - 72;
2139 if(k>=72 && k<=143)index = k - 71;
2149 double correction(1.);
2182 double maxEnergy = -9999.;
2185 std::pair<DetId, double> myPair (
DetId(0), -9999.);
2188 for( std::vector<std::pair<DetId,float> >::const_iterator idIt=mySCRecHits.begin(); idIt != mySCRecHits.end(); idIt++){
2192 hottestRecHit = & (* ( ebhits->
find((*idIt).first) ) );
2194 if( hottestRecHit == & (*( ebhits->
end())) )
2196 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2200 else if (idIt->first.subdetId() ==
EcalEndcap )
2202 hottestRecHit = & (* ( eehits->
find((*idIt).first) ) );
2203 if( hottestRecHit == & (*( eehits->
end())) )
2205 std::cout<<
"@@@@@@@@@@@@@@@@@@@@@@@@@@@ NO RECHIT FOUND SHOULD NEVER HAPPEN"<<std::endl;
2212 if(hottestRecHit && hottestRecHit->
energy() > maxEnergy){
2214 maxEnergy = hottestRecHit->
energy();
2216 myPair.first = hottestRecHit ->
id();
2217 myPair.second = maxEnergy;
2234 short ieta = myEBDetId.
ieta();
2235 short iphi = myEBDetId.
iphi();
2239 myBool = (
abs( ieta ) == 1 ||
abs( ieta ) == 25
2240 ||
abs( ieta ) ==26 ||
abs( ieta ) == 45
2241 ||
abs( ieta ) ==46 ||
abs( ieta ) == 65
2242 ||
abs( ieta ) ==66 ||
abs( ieta ) == 85 );
2244 for(
int i = 0;
i < 19;
i++){
2246 if(iphi == ( 20*
i + 1 ) || iphi == 20*
i )
2262 bool isNearCrack =
false;
2274 dist +=
pow( v1[
i]-v2[
i], 2 );
2330 for (
int i=0;
i<2;
i++)
2388 ofstream
fout(
"ZeeStatistics.txt");
2391 std::cout <<
"Cannot open output file.\n";
2394 fout<<
"ZeeStatistics"<<std::endl;
2396 fout<<
"##########################RECO#########################"<<std::endl;
2397 fout<<
"##################Zee with Barrel-Barrel electrons: "<<
BBZN<<std::endl;
2400 fout<<
"##################Zee with Barrel-Endcap electrons: "<<
EBZN<<std::endl;
2402 fout<<
"##################Zee with Endcap-Endcap electrons: "<<
EEZN<<std::endl;
2405 fout<<
"\n"<<std::endl;
2407 fout<<
"##########################GEN#########################"<<std::endl;
2408 fout<<
"##################Zee with Barrel-Barrel electrons: "<<(float)
MCZBB/
NEVT<<std::endl;
2409 fout<<
"##################Zee with Barrel-Endcap electrons: "<<(float)
MCZEB/
NEVT<<std::endl;
2410 fout<<
"##################Zee with Endcap-Endcap electrons: "<<(float)
MCZEE/
NEVT<<std::endl;
static bool validDetId(int i, int j)
check if a valid index combination
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_
T perp() const
Magnitude of transverse component.
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()
Exp< T >::type exp(const T &t)
virtual Status duringLoop(const edm::Event &, const edm::EventSetup &)
Called at each event.
SuperClusterRef superCluster() const
reference to a SuperCluster
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
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
std::string rechitCollection_
virtual double eta() const
momentum pseudorapidity
void setWhatProduced(T *iThis, const es::Label &iLabel=es::Label())
TH1F * h1_occupancyBarrel_
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]
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
int TOTAL_ELECTRONS_IN_ENDCAP
int BARREL_ELECTRONS_BEFORE_BORDER_CUT
void fillHLTInfo(edm::Handle< edm::TriggerResults >)
int SHOWER_ELECTRONS_IN_ENDCAP
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]
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
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)
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.
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 *)
virtual double phi() const
momentum azimuthal angle
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 double par[8 *NPar][4]
const std::vector< int > & getOptimizedIterations() const
int TOTAL_ELECTRONS_IN_BARREL
TH1F * h1_weightSumMeanEndcap_