37 #define TWOPI 6.283185308 71 if ( clusPhiSize_ % 2 == 0 || clusEtaSize_ % 2 == 0)
72 edm::LogError(
"AlCaPi0RecHitsProducerError") <<
"Size of eta/phi for simple clustering should be odd numbers";
236 hMinvPi0EB_ = ibooker.
book1D(
"Pi0InvmassEB",
"Pi0 Invariant Mass in EB",100,0.,0.5);
239 hMinvPi0EE_ = ibooker.
book1D(
"Pi0InvmassEE",
"Pi0 Invariant Mass in EE",100,0.,0.5);
242 hMinvEtaEB_ = ibooker.
book1D(
"EtaInvmassEB",
"Eta Invariant Mass in EB",100,0.,0.85);
245 hMinvEtaEE_ = ibooker.
book1D(
"EtaInvmassEE",
"Eta Invariant Mass in EE",100,0.,0.85);
249 hPt1Pi0EB_ = ibooker.
book1D(
"Pt1Pi0EB",
"Pt 1st most energetic Pi0 photon in EB",100,0.,20.);
252 hPt1Pi0EE_ = ibooker.
book1D(
"Pt1Pi0EE",
"Pt 1st most energetic Pi0 photon in EE",100,0.,20.);
255 hPt1EtaEB_ = ibooker.
book1D(
"Pt1EtaEB",
"Pt 1st most energetic Eta photon in EB",100,0.,20.);
258 hPt1EtaEE_ = ibooker.
book1D(
"Pt1EtaEE",
"Pt 1st most energetic Eta photon in EE",100,0.,20.);
261 hPt2Pi0EB_ = ibooker.
book1D(
"Pt2Pi0EB",
"Pt 2nd most energetic Pi0 photon in EB",100,0.,20.);
264 hPt2Pi0EE_ = ibooker.
book1D(
"Pt2Pi0EE",
"Pt 2nd most energetic Pi0 photon in EE",100,0.,20.);
267 hPt2EtaEB_ = ibooker.
book1D(
"Pt2EtaEB",
"Pt 2nd most energetic Eta photon in EB",100,0.,20.);
270 hPt2EtaEE_ = ibooker.
book1D(
"Pt2EtaEE",
"Pt 2nd most energetic Eta photon in EE",100,0.,20.);
298 hS4S91Pi0EB_ = ibooker.
book1D(
"S4S91Pi0EB",
"S4S9 1st most energetic Pi0 photon in EB",50,0.,1.);
301 hS4S91Pi0EE_ = ibooker.
book1D(
"S4S91Pi0EE",
"S4S9 1st most energetic Pi0 photon in EE",50,0.,1.);
304 hS4S91EtaEB_ = ibooker.
book1D(
"S4S91EtaEB",
"S4S9 1st most energetic Eta photon in EB",50,0.,1.);
307 hS4S91EtaEE_ = ibooker.
book1D(
"S4S91EtaEE",
"S4S9 1st most energetic Eta photon in EE",50,0.,1.);
310 hS4S92Pi0EB_ = ibooker.
book1D(
"S4S92Pi0EB",
"S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
313 hS4S92Pi0EE_ = ibooker.
book1D(
"S4S92Pi0EE",
"S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
316 hS4S92EtaEB_ = ibooker.
book1D(
"S4S92EtaEB",
"S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
319 hS4S92EtaEE_ = ibooker.
book1D(
"S4S92EtaEE",
"S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
334 std::vector<EcalRecHit> seeds;
337 vector<EBDetId> usedXtals;
378 for(itb=rhEBpi0->
begin(); itb!=rhEBpi0->
end(); ++itb){
381 double energy = itb->energy();
396 etot+= itb->energy();
413 vector<float> etClus;
414 vector<float> etaClus;
415 vector<float> thetaClus;
416 vector<float> phiClus;
417 vector<EBDetId> max_hit;
419 vector< vector<EcalRecHit> > RecHitsCluster;
420 vector< vector<EcalRecHit> > RecHitsCluster5x5;
421 vector<float> s4s9Clus;
422 vector<float> s9s25Clus;
431 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
432 EBDetId seed_id = itseed->id();
433 std::vector<EBDetId>::const_iterator usedIds;
435 bool seedAlreadyUsed=
false;
436 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
437 if(*usedIds==seed_id){
438 seedAlreadyUsed=
true;
443 if(seedAlreadyUsed)
continue;
445 std::vector<std::pair<DetId,float> > clus_used;
448 vector<EcalRecHit> RecHitsInWindow;
449 vector<EcalRecHit> RecHitsInWindow5x5;
451 double simple_energy = 0;
453 for (std::vector<DetId >::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
456 bool HitAlreadyUsed=
false;
457 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
463 if(HitAlreadyUsed)
continue;
470 usedXtals.push_back(*det);
471 RecHitsInWindow.push_back(
EBRecHits[nn]);
472 clus_used.push_back(std::make_pair(*det,1));
473 simple_energy = simple_energy +
EBRecHits[
nn].energy();
478 if(simple_energy <= 0)
continue;
485 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
491 float et_s = simple_energy *
sin(theta_s);
498 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
500 int seed_ieta = seed_id.
ieta();
501 int seed_iphi = seed_id.
iphi();
508 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
511 int ieta = det.
ieta();
512 int iphi = det.
iphi();
516 float en = RecHitsInWindow[j].energy();
521 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
522 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
523 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
524 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
532 if(e3x3 <= 0)
continue;
534 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
538 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
539 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
543 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
546 if(itdet0 != usedXtals.end())
continue;
554 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
561 if(e5x5 <= 0)
continue;
563 eClus.push_back(simple_energy);
564 etClus.push_back(et_s);
565 etaClus.push_back(clus_pos.eta());
566 thetaClus.push_back(theta_s);
567 phiClus.push_back(clus_pos.phi());
568 s4s9Clus.push_back(s4s9_max);
569 s9s25Clus.push_back(e3x3/e5x5);
570 RecHitsCluster.push_back(RecHitsInWindow);
571 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
588 for(Int_t
i=0 ;
i<nClus ;
i++){
589 for(Int_t j=
i+1 ; j<nClus ; j++){
594 float p0x = etClus[
i] *
cos(phiClus[
i]);
595 float p1x = etClus[j] *
cos(phiClus[j]);
596 float p0y = etClus[
i] *
sin(phiClus[i]);
597 float p1y = etClus[j] *
sin(phiClus[j]);
598 float p0z = eClus[
i] *
cos(thetaClus[i]);
599 float p1z = eClus[j] *
cos(thetaClus[j]);
602 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
606 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
613 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
614 for(Int_t
k=0 ;
k<nClus ;
k++){
619 if(
k==i ||
k==j)
continue;
620 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
622 float dretacl = fabs(etaClus[k] - pairVect.Eta());
623 float drcl = ClusVect.DeltaR(pairVect);
627 Iso = Iso + etClus[
k];
628 IsoClus.push_back(k);
675 for(itb=rhEBeta->
begin(); itb!=rhEBeta->
end(); ++itb){
678 double energy = itb->energy();
693 etot+= itb->energy();
710 vector<float> etClus;
711 vector<float> etaClus;
712 vector<float> thetaClus;
713 vector<float> phiClus;
714 vector<EBDetId> max_hit;
716 vector< vector<EcalRecHit> > RecHitsCluster;
717 vector< vector<EcalRecHit> > RecHitsCluster5x5;
718 vector<float> s4s9Clus;
719 vector<float> s9s25Clus;
727 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
728 EBDetId seed_id = itseed->id();
729 std::vector<EBDetId>::const_iterator usedIds;
731 bool seedAlreadyUsed=
false;
732 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
733 if(*usedIds==seed_id){
734 seedAlreadyUsed=
true;
739 if(seedAlreadyUsed)
continue;
741 std::vector<std::pair<DetId,float> > clus_used;
744 vector<EcalRecHit> RecHitsInWindow;
745 vector<EcalRecHit> RecHitsInWindow5x5;
747 double simple_energy = 0;
749 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
752 bool HitAlreadyUsed=
false;
753 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
759 if(HitAlreadyUsed)
continue;
766 usedXtals.push_back(*det);
767 RecHitsInWindow.push_back(
EBRecHits[nn]);
768 clus_used.push_back(std::make_pair(*det,1));
769 simple_energy = simple_energy +
EBRecHits[
nn].energy();
774 if(simple_energy <= 0)
continue;
781 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
787 float et_s = simple_energy *
sin(theta_s);
794 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
796 int seed_ieta = seed_id.
ieta();
797 int seed_iphi = seed_id.
iphi();
804 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
807 int ieta = det.
ieta();
808 int iphi = det.
iphi();
812 float en = RecHitsInWindow[j].energy();
817 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
818 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
819 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
820 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
828 if(e3x3 <= 0)
continue;
830 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
834 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
835 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
839 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
842 if(itdet0 != usedXtals.end())
continue;
850 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
857 if(e5x5 <= 0)
continue;
859 eClus.push_back(simple_energy);
860 etClus.push_back(et_s);
861 etaClus.push_back(clus_pos.eta());
862 thetaClus.push_back(theta_s);
863 phiClus.push_back(clus_pos.phi());
864 s4s9Clus.push_back(s4s9_max);
865 s9s25Clus.push_back(e3x3/e5x5);
866 RecHitsCluster.push_back(RecHitsInWindow);
867 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
884 for(Int_t
i=0 ;
i<nClus ;
i++){
885 for(Int_t j=
i+1 ; j<nClus ; j++){
890 float p0x = etClus[
i] *
cos(phiClus[
i]);
891 float p1x = etClus[j] *
cos(phiClus[j]);
892 float p0y = etClus[
i] *
sin(phiClus[i]);
893 float p1y = etClus[j] *
sin(phiClus[j]);
894 float p0z = eClus[
i] *
cos(thetaClus[i]);
895 float p1z = eClus[j] *
cos(thetaClus[j]);
898 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
902 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
909 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
910 for(Int_t
k=0 ;
k<nClus ;
k++){
915 if(
k==i ||
k==j)
continue;
916 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
918 float dretacl = fabs(etaClus[k] - pairVect.Eta());
919 float drcl = ClusVect.DeltaR(pairVect);
923 Iso = Iso + etClus[
k];
924 IsoClus.push_back(k);
985 std::vector<EcalRecHit> seedsEndCap;
988 vector<EEDetId> usedXtalsEndCap;
989 usedXtalsEndCap.clear();
994 for (ite=rhEEpi0->
begin(); ite!=rhEEpi0->
end(); ite++) {
995 double energy = ite->energy();
1011 etot+= ite->energy();
1021 vector<float> eClusEndCap;
1022 vector<float> etClusEndCap;
1023 vector<float> etaClusEndCap;
1024 vector<float> thetaClusEndCap;
1025 vector<float> phiClusEndCap;
1026 vector< vector<EcalRecHit> > RecHitsClusterEndCap;
1027 vector< vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1028 vector<float> s4s9ClusEndCap;
1029 vector<float> s9s25ClusEndCap;
1038 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1039 EEDetId seed_id = itseed->id();
1040 std::vector<EEDetId>::const_iterator usedIds;
1042 bool seedAlreadyUsed=
false;
1043 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1044 if(*usedIds==seed_id){
1045 seedAlreadyUsed=
true;
1050 if(seedAlreadyUsed)
continue;
1052 std::vector<std::pair<DetId,float> > clus_used;
1055 vector<EcalRecHit> RecHitsInWindow;
1056 vector<EcalRecHit> RecHitsInWindow5x5;
1058 float simple_energy = 0;
1060 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1063 bool HitAlreadyUsed=
false;
1064 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1066 HitAlreadyUsed=
true;
1071 if(HitAlreadyUsed)
continue;
1078 usedXtalsEndCap.push_back(*det);
1079 RecHitsInWindow.push_back(
EERecHits[nn]);
1080 clus_used.push_back(std::make_pair(*det,1));
1081 simple_energy = simple_energy +
EERecHits[
nn].energy();
1086 if( simple_energy <= 0)
continue;
1091 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1092 float et_s = simple_energy *
sin(theta_s);
1101 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1103 int ixSeed = seed_id.
ix();
1104 int iySeed = seed_id.
iy();
1108 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
1110 int dx = ixSeed - det_this.
ix();
1111 int dy = iySeed - det_this.
iy();
1113 float en = RecHitsInWindow[j].energy();
1115 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1116 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1117 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1118 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1126 if(e3x3 <= 0)
continue;
1128 eClusEndCap.push_back(simple_energy);
1129 etClusEndCap.push_back(et_s);
1130 etaClusEndCap.push_back(clus_pos.eta());
1131 thetaClusEndCap.push_back(theta_s);
1132 phiClusEndCap.push_back(clus_pos.phi());
1133 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1134 s9s25ClusEndCap.push_back(e3x3/e5x5);
1135 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1136 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1154 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1155 for(Int_t j=
i+1 ; j<nClusEndCap ; j++){
1159 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1160 float p1x = etClusEndCap[j] *
cos(phiClusEndCap[j]);
1161 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1162 float p1y = etClusEndCap[j] *
sin(phiClusEndCap[j]);
1163 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1164 float p1z = eClusEndCap[j] *
cos(thetaClusEndCap[j]);
1167 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1169 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1176 vector<int> IsoClus;
1179 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1180 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1185 if(
k==i ||
k==j)
continue;
1188 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1189 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1190 float drcl = clusVect.DeltaR(pairVect);
1193 Iso = Iso + etClusEndCap[
k];
1194 IsoClus.push_back(k);
1239 std::vector<EcalRecHit> seedsEndCap;
1240 seedsEndCap.clear();
1242 vector<EEDetId> usedXtalsEndCap;
1243 usedXtalsEndCap.clear();
1248 for (ite=rhEEeta->
begin(); ite!=rhEEeta->
end(); ite++) {
1249 double energy = ite->energy();
1265 etot+= ite->energy();
1275 vector<float> eClusEndCap;
1276 vector<float> etClusEndCap;
1277 vector<float> etaClusEndCap;
1278 vector<float> thetaClusEndCap;
1279 vector<float> phiClusEndCap;
1280 vector< vector<EcalRecHit> > RecHitsClusterEndCap;
1281 vector< vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1282 vector<float> s4s9ClusEndCap;
1283 vector<float> s9s25ClusEndCap;
1292 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1293 EEDetId seed_id = itseed->id();
1294 std::vector<EEDetId>::const_iterator usedIds;
1296 bool seedAlreadyUsed=
false;
1297 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1298 if(*usedIds==seed_id){
1299 seedAlreadyUsed=
true;
1304 if(seedAlreadyUsed)
continue;
1306 std::vector<std::pair<DetId,float> > clus_used;
1309 vector<EcalRecHit> RecHitsInWindow;
1310 vector<EcalRecHit> RecHitsInWindow5x5;
1312 float simple_energy = 0;
1314 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1317 bool HitAlreadyUsed=
false;
1318 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1320 HitAlreadyUsed=
true;
1325 if(HitAlreadyUsed)
continue;
1332 usedXtalsEndCap.push_back(*det);
1333 RecHitsInWindow.push_back(
EERecHits[nn]);
1334 clus_used.push_back(std::make_pair(*det,1));
1335 simple_energy = simple_energy +
EERecHits[
nn].energy();
1340 if( simple_energy <= 0)
continue;
1345 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1346 float et_s = simple_energy *
sin(theta_s);
1355 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1357 int ixSeed = seed_id.
ix();
1358 int iySeed = seed_id.
iy();
1362 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
1364 int dx = ixSeed - det_this.
ix();
1365 int dy = iySeed - det_this.
iy();
1367 float en = RecHitsInWindow[j].energy();
1369 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1370 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1371 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1372 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1380 if(e3x3 <= 0)
continue;
1382 eClusEndCap.push_back(simple_energy);
1383 etClusEndCap.push_back(et_s);
1384 etaClusEndCap.push_back(clus_pos.eta());
1385 thetaClusEndCap.push_back(theta_s);
1386 phiClusEndCap.push_back(clus_pos.phi());
1387 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1388 s9s25ClusEndCap.push_back(e3x3/e5x5);
1389 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1390 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1405 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1406 for(Int_t j=
i+1 ; j<nClusEndCap ; j++){
1410 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1411 float p1x = etClusEndCap[j] *
cos(phiClusEndCap[j]);
1412 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1413 float p1y = etClusEndCap[j] *
sin(phiClusEndCap[j]);
1414 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1415 float p1z = eClusEndCap[j] *
cos(thetaClusEndCap[j]);
1418 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1420 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1427 vector<int> IsoClus;
1430 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1431 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1436 if(
k==i ||
k==j)
continue;
1439 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1440 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1441 float drcl = clusVect.DeltaR(pairVect);
1444 Iso = Iso + etClusEndCap[
k];
1445 IsoClus.push_back(k);
1490 if(neta > 0) neta -= 1;
1491 if(nphi > 359) nphi=nphi-360;
1497 mdiff=(neta1-neta2);
1509 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)
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.
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.
double clusSeedThrEndCap_