40 #define TWOPI 6.283185308 74 if ( clusPhiSize_ % 2 == 0 || clusEtaSize_ % 2 == 0)
75 edm::LogError(
"AlCaPi0RecHitsProducerError") <<
"Size of eta/phi for simple clustering should be odd numbers";
239 hMinvPi0EB_ = ibooker.
book1D(
"Pi0InvmassEB",
"Pi0 Invariant Mass in EB",100,0.,0.5);
242 hMinvPi0EE_ = ibooker.
book1D(
"Pi0InvmassEE",
"Pi0 Invariant Mass in EE",100,0.,0.5);
245 hMinvEtaEB_ = ibooker.
book1D(
"EtaInvmassEB",
"Eta Invariant Mass in EB",100,0.,0.85);
248 hMinvEtaEE_ = ibooker.
book1D(
"EtaInvmassEE",
"Eta Invariant Mass in EE",100,0.,0.85);
252 hPt1Pi0EB_ = ibooker.
book1D(
"Pt1Pi0EB",
"Pt 1st most energetic Pi0 photon in EB",100,0.,20.);
255 hPt1Pi0EE_ = ibooker.
book1D(
"Pt1Pi0EE",
"Pt 1st most energetic Pi0 photon in EE",100,0.,20.);
258 hPt1EtaEB_ = ibooker.
book1D(
"Pt1EtaEB",
"Pt 1st most energetic Eta photon in EB",100,0.,20.);
261 hPt1EtaEE_ = ibooker.
book1D(
"Pt1EtaEE",
"Pt 1st most energetic Eta photon in EE",100,0.,20.);
264 hPt2Pi0EB_ = ibooker.
book1D(
"Pt2Pi0EB",
"Pt 2nd most energetic Pi0 photon in EB",100,0.,20.);
267 hPt2Pi0EE_ = ibooker.
book1D(
"Pt2Pi0EE",
"Pt 2nd most energetic Pi0 photon in EE",100,0.,20.);
270 hPt2EtaEB_ = ibooker.
book1D(
"Pt2EtaEB",
"Pt 2nd most energetic Eta photon in EB",100,0.,20.);
273 hPt2EtaEE_ = ibooker.
book1D(
"Pt2EtaEE",
"Pt 2nd most energetic Eta photon in EE",100,0.,20.);
301 hS4S91Pi0EB_ = ibooker.
book1D(
"S4S91Pi0EB",
"S4S9 1st most energetic Pi0 photon in EB",50,0.,1.);
304 hS4S91Pi0EE_ = ibooker.
book1D(
"S4S91Pi0EE",
"S4S9 1st most energetic Pi0 photon in EE",50,0.,1.);
307 hS4S91EtaEB_ = ibooker.
book1D(
"S4S91EtaEB",
"S4S9 1st most energetic Eta photon in EB",50,0.,1.);
310 hS4S91EtaEE_ = ibooker.
book1D(
"S4S91EtaEE",
"S4S9 1st most energetic Eta photon in EE",50,0.,1.);
313 hS4S92Pi0EB_ = ibooker.
book1D(
"S4S92Pi0EB",
"S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
316 hS4S92Pi0EE_ = ibooker.
book1D(
"S4S92Pi0EE",
"S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
319 hS4S92EtaEB_ = ibooker.
book1D(
"S4S92EtaEB",
"S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
322 hS4S92EtaEE_ = ibooker.
book1D(
"S4S92EtaEE",
"S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
353 std::vector<EcalRecHit> seeds;
356 vector<EBDetId> usedXtals;
397 for(itb=rhEBpi0->
begin(); itb!=rhEBpi0->
end(); ++itb){
400 double energy = itb->energy();
415 etot+= itb->energy();
432 vector<float> etClus;
433 vector<float> etaClus;
434 vector<float> thetaClus;
435 vector<float> phiClus;
436 vector<EBDetId> max_hit;
438 vector< vector<EcalRecHit> > RecHitsCluster;
439 vector< vector<EcalRecHit> > RecHitsCluster5x5;
440 vector<float> s4s9Clus;
441 vector<float> s9s25Clus;
450 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
451 EBDetId seed_id = itseed->id();
452 std::vector<EBDetId>::const_iterator usedIds;
454 bool seedAlreadyUsed=
false;
455 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
456 if(*usedIds==seed_id){
457 seedAlreadyUsed=
true;
462 if(seedAlreadyUsed)
continue;
464 std::vector<std::pair<DetId,float> > clus_used;
467 vector<EcalRecHit> RecHitsInWindow;
468 vector<EcalRecHit> RecHitsInWindow5x5;
470 double simple_energy = 0;
472 for (std::vector<DetId >::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
475 bool HitAlreadyUsed=
false;
476 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
482 if(HitAlreadyUsed)
continue;
489 usedXtals.push_back(*det);
490 RecHitsInWindow.push_back(
EBRecHits[nn]);
491 clus_used.push_back(std::make_pair(*det,1));
492 simple_energy = simple_energy +
EBRecHits[
nn].energy();
497 if(simple_energy <= 0)
continue;
504 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
510 float et_s = simple_energy *
sin(theta_s);
517 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
519 int seed_ieta = seed_id.
ieta();
520 int seed_iphi = seed_id.
iphi();
527 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
530 int ieta = det.
ieta();
531 int iphi = det.
iphi();
535 float en = RecHitsInWindow[j].energy();
540 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
541 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
542 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
543 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
551 if(e3x3 <= 0)
continue;
553 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
557 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
558 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
562 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
565 if(itdet0 != usedXtals.end())
continue;
573 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
580 if(e5x5 <= 0)
continue;
582 eClus.push_back(simple_energy);
583 etClus.push_back(et_s);
584 etaClus.push_back(clus_pos.eta());
585 thetaClus.push_back(theta_s);
586 phiClus.push_back(clus_pos.phi());
587 s4s9Clus.push_back(s4s9_max);
588 s9s25Clus.push_back(e3x3/e5x5);
589 RecHitsCluster.push_back(RecHitsInWindow);
590 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
607 for(Int_t
i=0 ;
i<nClus ;
i++){
608 for(Int_t j=
i+1 ; j<nClus ; j++){
613 float p0x = etClus[
i] *
cos(phiClus[
i]);
614 float p1x = etClus[j] *
cos(phiClus[j]);
615 float p0y = etClus[
i] *
sin(phiClus[i]);
616 float p1y = etClus[j] *
sin(phiClus[j]);
617 float p0z = eClus[
i] *
cos(thetaClus[i]);
618 float p1z = eClus[j] *
cos(thetaClus[j]);
621 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
625 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
632 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
633 for(Int_t
k=0 ;
k<nClus ;
k++){
638 if(
k==i ||
k==j)
continue;
639 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
641 float dretacl = fabs(etaClus[k] - pairVect.Eta());
642 float drcl = ClusVect.DeltaR(pairVect);
646 Iso = Iso + etClus[
k];
647 IsoClus.push_back(k);
694 for(itb=rhEBeta->
begin(); itb!=rhEBeta->
end(); ++itb){
697 double energy = itb->energy();
712 etot+= itb->energy();
729 vector<float> etClus;
730 vector<float> etaClus;
731 vector<float> thetaClus;
732 vector<float> phiClus;
733 vector<EBDetId> max_hit;
735 vector< vector<EcalRecHit> > RecHitsCluster;
736 vector< vector<EcalRecHit> > RecHitsCluster5x5;
737 vector<float> s4s9Clus;
738 vector<float> s9s25Clus;
746 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
747 EBDetId seed_id = itseed->id();
748 std::vector<EBDetId>::const_iterator usedIds;
750 bool seedAlreadyUsed=
false;
751 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
752 if(*usedIds==seed_id){
753 seedAlreadyUsed=
true;
758 if(seedAlreadyUsed)
continue;
760 std::vector<std::pair<DetId,float> > clus_used;
763 vector<EcalRecHit> RecHitsInWindow;
764 vector<EcalRecHit> RecHitsInWindow5x5;
766 double simple_energy = 0;
768 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
771 bool HitAlreadyUsed=
false;
772 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
778 if(HitAlreadyUsed)
continue;
785 usedXtals.push_back(*det);
786 RecHitsInWindow.push_back(
EBRecHits[nn]);
787 clus_used.push_back(std::make_pair(*det,1));
788 simple_energy = simple_energy +
EBRecHits[
nn].energy();
793 if(simple_energy <= 0)
continue;
800 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
806 float et_s = simple_energy *
sin(theta_s);
813 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
815 int seed_ieta = seed_id.
ieta();
816 int seed_iphi = seed_id.
iphi();
823 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
826 int ieta = det.
ieta();
827 int iphi = det.
iphi();
831 float en = RecHitsInWindow[j].energy();
836 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
837 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
838 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
839 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
847 if(e3x3 <= 0)
continue;
849 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
853 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
854 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
858 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
861 if(itdet0 != usedXtals.end())
continue;
869 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
876 if(e5x5 <= 0)
continue;
878 eClus.push_back(simple_energy);
879 etClus.push_back(et_s);
880 etaClus.push_back(clus_pos.eta());
881 thetaClus.push_back(theta_s);
882 phiClus.push_back(clus_pos.phi());
883 s4s9Clus.push_back(s4s9_max);
884 s9s25Clus.push_back(e3x3/e5x5);
885 RecHitsCluster.push_back(RecHitsInWindow);
886 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
903 for(Int_t
i=0 ;
i<nClus ;
i++){
904 for(Int_t j=
i+1 ; j<nClus ; j++){
909 float p0x = etClus[
i] *
cos(phiClus[
i]);
910 float p1x = etClus[j] *
cos(phiClus[j]);
911 float p0y = etClus[
i] *
sin(phiClus[i]);
912 float p1y = etClus[j] *
sin(phiClus[j]);
913 float p0z = eClus[
i] *
cos(thetaClus[i]);
914 float p1z = eClus[j] *
cos(thetaClus[j]);
917 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
921 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
928 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
929 for(Int_t
k=0 ;
k<nClus ;
k++){
934 if(
k==i ||
k==j)
continue;
935 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
937 float dretacl = fabs(etaClus[k] - pairVect.Eta());
938 float drcl = ClusVect.DeltaR(pairVect);
942 Iso = Iso + etClus[
k];
943 IsoClus.push_back(k);
1004 std::vector<EcalRecHit> seedsEndCap;
1005 seedsEndCap.clear();
1007 vector<EEDetId> usedXtalsEndCap;
1008 usedXtalsEndCap.clear();
1013 for (ite=rhEEpi0->
begin(); ite!=rhEEpi0->
end(); ite++) {
1014 double energy = ite->energy();
1030 etot+= ite->energy();
1040 vector<float> eClusEndCap;
1041 vector<float> etClusEndCap;
1042 vector<float> etaClusEndCap;
1043 vector<float> thetaClusEndCap;
1044 vector<float> phiClusEndCap;
1045 vector< vector<EcalRecHit> > RecHitsClusterEndCap;
1046 vector< vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1047 vector<float> s4s9ClusEndCap;
1048 vector<float> s9s25ClusEndCap;
1057 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1058 EEDetId seed_id = itseed->id();
1059 std::vector<EEDetId>::const_iterator usedIds;
1061 bool seedAlreadyUsed=
false;
1062 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1063 if(*usedIds==seed_id){
1064 seedAlreadyUsed=
true;
1069 if(seedAlreadyUsed)
continue;
1071 std::vector<std::pair<DetId,float> > clus_used;
1074 vector<EcalRecHit> RecHitsInWindow;
1075 vector<EcalRecHit> RecHitsInWindow5x5;
1077 float simple_energy = 0;
1079 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1082 bool HitAlreadyUsed=
false;
1083 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1085 HitAlreadyUsed=
true;
1090 if(HitAlreadyUsed)
continue;
1097 usedXtalsEndCap.push_back(*det);
1098 RecHitsInWindow.push_back(
EERecHits[nn]);
1099 clus_used.push_back(std::make_pair(*det,1));
1100 simple_energy = simple_energy +
EERecHits[
nn].energy();
1105 if( simple_energy <= 0)
continue;
1110 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1111 float et_s = simple_energy *
sin(theta_s);
1120 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1122 int ixSeed = seed_id.
ix();
1123 int iySeed = seed_id.
iy();
1127 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
1129 int dx = ixSeed - det_this.
ix();
1130 int dy = iySeed - det_this.
iy();
1132 float en = RecHitsInWindow[j].energy();
1134 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1135 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1136 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1137 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1145 if(e3x3 <= 0)
continue;
1147 eClusEndCap.push_back(simple_energy);
1148 etClusEndCap.push_back(et_s);
1149 etaClusEndCap.push_back(clus_pos.eta());
1150 thetaClusEndCap.push_back(theta_s);
1151 phiClusEndCap.push_back(clus_pos.phi());
1152 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1153 s9s25ClusEndCap.push_back(e3x3/e5x5);
1154 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1155 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1173 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1174 for(Int_t j=
i+1 ; j<nClusEndCap ; j++){
1178 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1179 float p1x = etClusEndCap[j] *
cos(phiClusEndCap[j]);
1180 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1181 float p1y = etClusEndCap[j] *
sin(phiClusEndCap[j]);
1182 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1183 float p1z = eClusEndCap[j] *
cos(thetaClusEndCap[j]);
1186 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1188 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1195 vector<int> IsoClus;
1198 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1199 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1204 if(
k==i ||
k==j)
continue;
1207 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1208 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1209 float drcl = clusVect.DeltaR(pairVect);
1212 Iso = Iso + etClusEndCap[
k];
1213 IsoClus.push_back(k);
1258 std::vector<EcalRecHit> seedsEndCap;
1259 seedsEndCap.clear();
1261 vector<EEDetId> usedXtalsEndCap;
1262 usedXtalsEndCap.clear();
1267 for (ite=rhEEeta->
begin(); ite!=rhEEeta->
end(); ite++) {
1268 double energy = ite->energy();
1284 etot+= ite->energy();
1294 vector<float> eClusEndCap;
1295 vector<float> etClusEndCap;
1296 vector<float> etaClusEndCap;
1297 vector<float> thetaClusEndCap;
1298 vector<float> phiClusEndCap;
1299 vector< vector<EcalRecHit> > RecHitsClusterEndCap;
1300 vector< vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1301 vector<float> s4s9ClusEndCap;
1302 vector<float> s9s25ClusEndCap;
1311 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1312 EEDetId seed_id = itseed->id();
1313 std::vector<EEDetId>::const_iterator usedIds;
1315 bool seedAlreadyUsed=
false;
1316 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1317 if(*usedIds==seed_id){
1318 seedAlreadyUsed=
true;
1323 if(seedAlreadyUsed)
continue;
1325 std::vector<std::pair<DetId,float> > clus_used;
1328 vector<EcalRecHit> RecHitsInWindow;
1329 vector<EcalRecHit> RecHitsInWindow5x5;
1331 float simple_energy = 0;
1333 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1336 bool HitAlreadyUsed=
false;
1337 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1339 HitAlreadyUsed=
true;
1344 if(HitAlreadyUsed)
continue;
1351 usedXtalsEndCap.push_back(*det);
1352 RecHitsInWindow.push_back(
EERecHits[nn]);
1353 clus_used.push_back(std::make_pair(*det,1));
1354 simple_energy = simple_energy +
EERecHits[
nn].energy();
1359 if( simple_energy <= 0)
continue;
1364 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1365 float et_s = simple_energy *
sin(theta_s);
1374 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1376 int ixSeed = seed_id.
ix();
1377 int iySeed = seed_id.
iy();
1381 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
1383 int dx = ixSeed - det_this.
ix();
1384 int dy = iySeed - det_this.
iy();
1386 float en = RecHitsInWindow[j].energy();
1388 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1389 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1390 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1391 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1399 if(e3x3 <= 0)
continue;
1401 eClusEndCap.push_back(simple_energy);
1402 etClusEndCap.push_back(et_s);
1403 etaClusEndCap.push_back(clus_pos.eta());
1404 thetaClusEndCap.push_back(theta_s);
1405 phiClusEndCap.push_back(clus_pos.phi());
1406 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1407 s9s25ClusEndCap.push_back(e3x3/e5x5);
1408 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1409 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1424 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1425 for(Int_t j=
i+1 ; j<nClusEndCap ; j++){
1429 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1430 float p1x = etClusEndCap[j] *
cos(phiClusEndCap[j]);
1431 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1432 float p1y = etClusEndCap[j] *
sin(phiClusEndCap[j]);
1433 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1434 float p1z = eClusEndCap[j] *
cos(thetaClusEndCap[j]);
1437 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1439 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1446 vector<int> IsoClus;
1449 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1450 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1455 if(
k==i ||
k==j)
continue;
1458 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1459 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1460 float drcl = clusVect.DeltaR(pairVect);
1463 Iso = Iso + etClusEndCap[
k];
1464 IsoClus.push_back(k);
1531 if(neta > 0) neta -= 1;
1532 if(nphi > 359) nphi=nphi-360;
1538 mdiff=(neta1-neta2);
1550 if(nphi1>nphi2) mdiff=-mdiff;
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEEeta_
double seleEtaBeltDREndCap_
void analyze(const edm::Event &e, const edm::EventSetup &c) override
T getParameter(std::string const &) const
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
MonitorElement * hiYDistrEEeta_
Distribution of rechits in iy EE (eta)
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * hMinvPi0EE_
Pi0 invariant mass in EE.
double seleMinvMinPi0EndCap_
MonitorElement * hiXDistrEEeta_
Distribution of rechits in ix EE (eta)
MonitorElement * hiPhiDistrEBpi0_
Distribution of rechits in iPhi (pi0)
MonitorElement * hPt1Pi0EB_
Pt of the 1st most energetic Pi0 photon in EB.
MonitorElement * hPt2Pi0EE_
Pt of the 2nd most energetic Pi0 photon in EE.
MonitorElement * hiYDistrEEpi0_
Distribution of rechits in iy EE (pi0)
void beginLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &context) override
double ptMinForIsolationEtaEndCap_
MonitorElement * hPtPi0EE_
Pi0 Pt in EE.
MonitorElement * hMinvEtaEB_
Eta invariant mass in EB.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
std::vector< EEDetId > detIdEERecHits
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
PositionCalc posCalculator_
double seleMinvMinEtaEndCap_
MonitorElement * hS4S92EtaEB_
S4S9 of the 2nd most energetic eta photon.
std::vector< EBDetId > detIdEBRecHits
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEBeta_
double selePtGammaEndCap_
for pi0->gg endcap
MonitorElement * hMeanRecHitEnergyEBeta_
Distribution of Mean energy per rechit EB (eta)
Sin< T >::type sin(const T &t)
double seleMinvMaxPi0EndCap_
std::string folderName_
DQM folder name.
std::vector< EcalRecHit >::const_iterator const_iterator
MonitorElement * hMeanRecHitEnergyEEpi0_
Distribution of Mean energy per rechit EE (pi0)
MonitorElement * hNRecHitsEBpi0_
Distribution of number of RecHits EB (pi0)
MonitorElement * hPt2EtaEE_
Pt of the 2nd most energetic Eta photon in EE.
MonitorElement * hRechitEnergyEBpi0_
Energy Distribution of rechits EB (pi0)
std::string fileName_
Output file name if required.
MonitorElement * hRechitEnergyEEpi0_
Energy Distribution of rechits EE (pi0)
double seleMinvMaxEtaEndCap_
std::vector< EcalRecHit > EERecHits
MonitorElement * hPtEtaEB_
Eta Pt in EB.
MonitorElement * hIsoEtaEB_
Eta Iso EB.
bool saveToFile_
Write to file.
double seleS4S9GammaEtaEndCap_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
MonitorElement * hIsoPi0EB_
Pi0 Iso EB.
MonitorElement * hS4S91Pi0EE_
S4S9 of the 1st most energetic pi0 photon EE.
MonitorElement * hS4S91EtaEB_
S4S9 of the 1st most energetic eta photon.
double ptMinForIsolationEta_
int iphi() const
get the crystal iphi
MonitorElement * hMinvPi0EB_
Pi0 invariant mass in EB.
MonitorElement * hS4S92Pi0EE_
S4S9 of the 2nd most energetic pi0 photon EE.
void convxtalid(int &, int &)
double ptMinForIsolationEndCap_
MonitorElement * hiEtaDistrEBpi0_
Distribution of rechits in iEta (pi0)
double selePi0BeltDetaEndCap_
MonitorElement * hS4S91Pi0EB_
S4S9 of the 1st most energetic pi0 photon.
MonitorElement * hiXDistrEEpi0_
Distribution of rechits in ix EE (pi0)
MonitorElement * hiEtaDistrEBeta_
Distribution of rechits in iEta (eta)
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEBpi0_
object to monitor
bool isMonEBpi0_
which subdet will be monitored
double selePtGammaEtaEndCap_
for eta->gg endcap
Cos< T >::type cos(const T &t)
int diff_nphi_s(int, int)
double ptMinForIsolation_
int diff_neta_s(int, int)
MonitorElement * book1D(Args &&...args)
double seleXtalMinEnergy_
Abs< T >::type abs(const T &t)
DQMSourcePi0(const edm::ParameterSet &)
int ieta() const
get the crystal ieta
MonitorElement * hPtEtaEE_
Eta Pt in EE.
double seleXtalMinEnergyEndCap_
double seleS9S25GammaEtaEndCap_
MonitorElement * hEventEnergyEEpi0_
Distribution of total event energy EE (pi0)
const_iterator end() const
MonitorElement * hNRecHitsEEeta_
Distribution of number of RecHits EE (eta)
double selePtGammaEta_
for eta->gg barrel
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEEpi0_
object to monitor
void setCurrentFolder(const std::string &fullpath)
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
double seleEtaBeltDetaEndCap_
T const * product() const
XYZPointD XYZPoint
point in space with cartesian internal representation
MonitorElement * hiPhiDistrEBeta_
Distribution of rechits in iPhi (eta)
MonitorElement * hRechitEnergyEBeta_
Energy Distribution of rechits EB (eta)
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
MonitorElement * hPt1EtaEE_
Pt of the 1st most energetic Eta photon in EE.
MonitorElement * hEventEnergyEBpi0_
Distribution of total event energy EB (pi0)
MonitorElement * hPt1EtaEB_
Pt of the 1st most energetic Eta photon in EB.
MonitorElement * hEventEnergyEBeta_
Distribution of total event energy EB (eta)
MonitorElement * hIsoEtaEE_
Eta Iso EE.
MonitorElement * hS4S92Pi0EB_
S4S9 of the 2nd most energetic pi0 photon.
MonitorElement * hPt2EtaEB_
Pt of the 2nd most energetic Eta photon in EB.
std::vector< EcalRecHit > EBRecHits
MonitorElement * hMinvEtaEE_
Eta invariant mass in EE.
MonitorElement * hPtPi0EB_
Pi0 Pt in EB.
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=0)
MonitorElement * hS4S92EtaEE_
S4S9 of the 2nd most energetic eta photon EE.
unsigned int prescaleFactor_
Monitor every prescaleFactor_ events.
MonitorElement * hMeanRecHitEnergyEEeta_
Distribution of Mean energy per rechit EE (eta)
MonitorElement * hIsoPi0EE_
Pi0 Iso EE.
double seleS4S9GammaEndCap_
MonitorElement * hNRecHitsEBeta_
Distribution of number of RecHits EB (eta)
MonitorElement * hPt1Pi0EE_
Pt of the 1st most energetic Pi0 photon in EE.
MonitorElement * hS4S91EtaEE_
S4S9 of the 1st most energetic eta photon EE.
void endLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c) override
MonitorElement * hRechitEnergyEEeta_
Energy Distribution of rechits EE (eta)
MonitorElement * hMeanRecHitEnergyEBpi0_
Distribution of Mean energy per rechit EB (pi0)
MonitorElement * hEventEnergyEEeta_
Distribution of total event energy EE (eta)
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
MonitorElement * hNRecHitsEEpi0_
Distribution of number of RecHits EE (pi0)
double selePi0BeltDREndCap_
const_iterator begin() const
double seleS9S25GammaEta_
MonitorElement * hPt2Pi0EB_
Pt of the 2nd most energetic Pi0 photon in EB.
void endRun(const edm::Run &r, const edm::EventSetup &c) override
double clusSeedThrEndCap_