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;
404 for(itb=rhEBpi0->
begin(); itb!=rhEBpi0->
end(); ++itb){
407 double energy = itb->energy();
422 etot+= itb->energy();
439 vector<float> etClus;
440 vector<float> etaClus;
441 vector<float> thetaClus;
442 vector<float> phiClus;
443 vector<EBDetId> max_hit;
445 vector< vector<EcalRecHit> > RecHitsCluster;
446 vector< vector<EcalRecHit> > RecHitsCluster5x5;
447 vector<float> s4s9Clus;
448 vector<float> s9s25Clus;
457 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
458 EBDetId seed_id = itseed->id();
459 std::vector<EBDetId>::const_iterator usedIds;
461 bool seedAlreadyUsed=
false;
462 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
463 if(*usedIds==seed_id){
464 seedAlreadyUsed=
true;
469 if(seedAlreadyUsed)
continue;
471 std::vector<std::pair<DetId,float> > clus_used;
474 vector<EcalRecHit> RecHitsInWindow;
475 vector<EcalRecHit> RecHitsInWindow5x5;
477 double simple_energy = 0;
479 for (std::vector<DetId >::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
482 bool HitAlreadyUsed=
false;
483 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
489 if(HitAlreadyUsed)
continue;
496 usedXtals.push_back(*det);
497 RecHitsInWindow.push_back(
EBRecHits[nn]);
498 clus_used.push_back(std::make_pair(*det,1));
499 simple_energy = simple_energy +
EBRecHits[
nn].energy();
504 if(simple_energy <= 0)
continue;
511 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
517 float et_s = simple_energy *
sin(theta_s);
524 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
526 int seed_ieta = seed_id.
ieta();
527 int seed_iphi = seed_id.
iphi();
534 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
537 int ieta = det.
ieta();
538 int iphi = det.
iphi();
542 float en = RecHitsInWindow[j].energy();
547 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
548 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
549 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
550 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
558 if(e3x3 <= 0)
continue;
560 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
564 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
565 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
569 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
572 if(itdet0 != usedXtals.end())
continue;
580 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
587 if(e5x5 <= 0)
continue;
589 eClus.push_back(simple_energy);
590 etClus.push_back(et_s);
591 etaClus.push_back(clus_pos.eta());
592 thetaClus.push_back(theta_s);
593 phiClus.push_back(clus_pos.phi());
594 s4s9Clus.push_back(s4s9_max);
595 s9s25Clus.push_back(e3x3/e5x5);
596 RecHitsCluster.push_back(RecHitsInWindow);
597 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
614 for(Int_t
i=0 ;
i<nClus ;
i++){
615 for(Int_t j=
i+1 ; j<nClus ; j++){
620 float p0x = etClus[
i] *
cos(phiClus[
i]);
621 float p1x = etClus[j] *
cos(phiClus[j]);
622 float p0y = etClus[
i] *
sin(phiClus[i]);
623 float p1y = etClus[j] *
sin(phiClus[j]);
624 float p0z = eClus[
i] *
cos(thetaClus[i]);
625 float p1z = eClus[j] *
cos(thetaClus[j]);
628 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
632 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
639 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
640 for(Int_t
k=0 ;
k<nClus ;
k++){
645 if(
k==i ||
k==j)
continue;
646 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
648 float dretacl = fabs(etaClus[k] - pairVect.Eta());
649 float drcl = ClusVect.DeltaR(pairVect);
653 Iso = Iso + etClus[
k];
654 IsoClus.push_back(k);
701 for(itb=rhEBeta->
begin(); itb!=rhEBeta->
end(); ++itb){
704 double energy = itb->energy();
719 etot+= itb->energy();
736 vector<float> etClus;
737 vector<float> etaClus;
738 vector<float> thetaClus;
739 vector<float> phiClus;
740 vector<EBDetId> max_hit;
742 vector< vector<EcalRecHit> > RecHitsCluster;
743 vector< vector<EcalRecHit> > RecHitsCluster5x5;
744 vector<float> s4s9Clus;
745 vector<float> s9s25Clus;
753 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
754 EBDetId seed_id = itseed->id();
755 std::vector<EBDetId>::const_iterator usedIds;
757 bool seedAlreadyUsed=
false;
758 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
759 if(*usedIds==seed_id){
760 seedAlreadyUsed=
true;
765 if(seedAlreadyUsed)
continue;
767 std::vector<std::pair<DetId,float> > clus_used;
770 vector<EcalRecHit> RecHitsInWindow;
771 vector<EcalRecHit> RecHitsInWindow5x5;
773 double simple_energy = 0;
775 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
778 bool HitAlreadyUsed=
false;
779 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
785 if(HitAlreadyUsed)
continue;
792 usedXtals.push_back(*det);
793 RecHitsInWindow.push_back(
EBRecHits[nn]);
794 clus_used.push_back(std::make_pair(*det,1));
795 simple_energy = simple_energy +
EBRecHits[
nn].energy();
800 if(simple_energy <= 0)
continue;
807 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
813 float et_s = simple_energy *
sin(theta_s);
820 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
822 int seed_ieta = seed_id.
ieta();
823 int seed_iphi = seed_id.
iphi();
830 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
833 int ieta = det.
ieta();
834 int iphi = det.
iphi();
838 float en = RecHitsInWindow[j].energy();
843 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
844 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
845 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
846 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
854 if(e3x3 <= 0)
continue;
856 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
860 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
861 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
865 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
868 if(itdet0 != usedXtals.end())
continue;
876 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
883 if(e5x5 <= 0)
continue;
885 eClus.push_back(simple_energy);
886 etClus.push_back(et_s);
887 etaClus.push_back(clus_pos.eta());
888 thetaClus.push_back(theta_s);
889 phiClus.push_back(clus_pos.phi());
890 s4s9Clus.push_back(s4s9_max);
891 s9s25Clus.push_back(e3x3/e5x5);
892 RecHitsCluster.push_back(RecHitsInWindow);
893 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
910 for(Int_t
i=0 ;
i<nClus ;
i++){
911 for(Int_t j=
i+1 ; j<nClus ; j++){
916 float p0x = etClus[
i] *
cos(phiClus[
i]);
917 float p1x = etClus[j] *
cos(phiClus[j]);
918 float p0y = etClus[
i] *
sin(phiClus[i]);
919 float p1y = etClus[j] *
sin(phiClus[j]);
920 float p0z = eClus[
i] *
cos(thetaClus[i]);
921 float p1z = eClus[j] *
cos(thetaClus[j]);
924 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
928 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
935 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
936 for(Int_t
k=0 ;
k<nClus ;
k++){
941 if(
k==i ||
k==j)
continue;
942 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
944 float dretacl = fabs(etaClus[k] - pairVect.Eta());
945 float drcl = ClusVect.DeltaR(pairVect);
949 Iso = Iso + etClus[
k];
950 IsoClus.push_back(k);
1011 std::vector<EcalRecHit> seedsEndCap;
1012 seedsEndCap.clear();
1014 vector<EEDetId> usedXtalsEndCap;
1015 usedXtalsEndCap.clear();
1020 for (ite=rhEEpi0->
begin(); ite!=rhEEpi0->
end(); ite++) {
1021 double energy = ite->energy();
1037 etot+= ite->energy();
1047 vector<float> eClusEndCap;
1048 vector<float> etClusEndCap;
1049 vector<float> etaClusEndCap;
1050 vector<float> thetaClusEndCap;
1051 vector<float> phiClusEndCap;
1052 vector< vector<EcalRecHit> > RecHitsClusterEndCap;
1053 vector< vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1054 vector<float> s4s9ClusEndCap;
1055 vector<float> s9s25ClusEndCap;
1064 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1065 EEDetId seed_id = itseed->id();
1066 std::vector<EEDetId>::const_iterator usedIds;
1068 bool seedAlreadyUsed=
false;
1069 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1070 if(*usedIds==seed_id){
1071 seedAlreadyUsed=
true;
1076 if(seedAlreadyUsed)
continue;
1078 std::vector<std::pair<DetId,float> > clus_used;
1081 vector<EcalRecHit> RecHitsInWindow;
1082 vector<EcalRecHit> RecHitsInWindow5x5;
1084 float simple_energy = 0;
1086 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1089 bool HitAlreadyUsed=
false;
1090 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1092 HitAlreadyUsed=
true;
1097 if(HitAlreadyUsed)
continue;
1104 usedXtalsEndCap.push_back(*det);
1105 RecHitsInWindow.push_back(
EERecHits[nn]);
1106 clus_used.push_back(std::make_pair(*det,1));
1107 simple_energy = simple_energy +
EERecHits[
nn].energy();
1112 if( simple_energy <= 0)
continue;
1117 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1118 float et_s = simple_energy *
sin(theta_s);
1127 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1129 int ixSeed = seed_id.
ix();
1130 int iySeed = seed_id.
iy();
1134 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
1136 int dx = ixSeed - det_this.
ix();
1137 int dy = iySeed - det_this.
iy();
1139 float en = RecHitsInWindow[j].energy();
1141 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1142 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1143 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1144 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1152 if(e3x3 <= 0)
continue;
1154 eClusEndCap.push_back(simple_energy);
1155 etClusEndCap.push_back(et_s);
1156 etaClusEndCap.push_back(clus_pos.eta());
1157 thetaClusEndCap.push_back(theta_s);
1158 phiClusEndCap.push_back(clus_pos.phi());
1159 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1160 s9s25ClusEndCap.push_back(e3x3/e5x5);
1161 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1162 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1180 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1181 for(Int_t j=
i+1 ; j<nClusEndCap ; j++){
1185 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1186 float p1x = etClusEndCap[j] *
cos(phiClusEndCap[j]);
1187 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1188 float p1y = etClusEndCap[j] *
sin(phiClusEndCap[j]);
1189 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1190 float p1z = eClusEndCap[j] *
cos(thetaClusEndCap[j]);
1193 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1195 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1202 vector<int> IsoClus;
1205 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1206 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1211 if(
k==i ||
k==j)
continue;
1214 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1215 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1216 float drcl = clusVect.DeltaR(pairVect);
1219 Iso = Iso + etClusEndCap[
k];
1220 IsoClus.push_back(k);
1265 std::vector<EcalRecHit> seedsEndCap;
1266 seedsEndCap.clear();
1268 vector<EEDetId> usedXtalsEndCap;
1269 usedXtalsEndCap.clear();
1274 for (ite=rhEEeta->
begin(); ite!=rhEEeta->
end(); ite++) {
1275 double energy = ite->energy();
1291 etot+= ite->energy();
1301 vector<float> eClusEndCap;
1302 vector<float> etClusEndCap;
1303 vector<float> etaClusEndCap;
1304 vector<float> thetaClusEndCap;
1305 vector<float> phiClusEndCap;
1306 vector< vector<EcalRecHit> > RecHitsClusterEndCap;
1307 vector< vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1308 vector<float> s4s9ClusEndCap;
1309 vector<float> s9s25ClusEndCap;
1318 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1319 EEDetId seed_id = itseed->id();
1320 std::vector<EEDetId>::const_iterator usedIds;
1322 bool seedAlreadyUsed=
false;
1323 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1324 if(*usedIds==seed_id){
1325 seedAlreadyUsed=
true;
1330 if(seedAlreadyUsed)
continue;
1332 std::vector<std::pair<DetId,float> > clus_used;
1335 vector<EcalRecHit> RecHitsInWindow;
1336 vector<EcalRecHit> RecHitsInWindow5x5;
1338 float simple_energy = 0;
1340 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1343 bool HitAlreadyUsed=
false;
1344 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1346 HitAlreadyUsed=
true;
1351 if(HitAlreadyUsed)
continue;
1358 usedXtalsEndCap.push_back(*det);
1359 RecHitsInWindow.push_back(
EERecHits[nn]);
1360 clus_used.push_back(std::make_pair(*det,1));
1361 simple_energy = simple_energy +
EERecHits[
nn].energy();
1366 if( simple_energy <= 0)
continue;
1371 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1372 float et_s = simple_energy *
sin(theta_s);
1381 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1383 int ixSeed = seed_id.
ix();
1384 int iySeed = seed_id.
iy();
1388 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
1390 int dx = ixSeed - det_this.
ix();
1391 int dy = iySeed - det_this.
iy();
1393 float en = RecHitsInWindow[j].energy();
1395 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1396 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1397 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1398 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1406 if(e3x3 <= 0)
continue;
1408 eClusEndCap.push_back(simple_energy);
1409 etClusEndCap.push_back(et_s);
1410 etaClusEndCap.push_back(clus_pos.eta());
1411 thetaClusEndCap.push_back(theta_s);
1412 phiClusEndCap.push_back(clus_pos.phi());
1413 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1414 s9s25ClusEndCap.push_back(e3x3/e5x5);
1415 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1416 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1431 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1432 for(Int_t j=
i+1 ; j<nClusEndCap ; j++){
1436 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1437 float p1x = etClusEndCap[j] *
cos(phiClusEndCap[j]);
1438 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1439 float p1y = etClusEndCap[j] *
sin(phiClusEndCap[j]);
1440 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1441 float p1z = eClusEndCap[j] *
cos(thetaClusEndCap[j]);
1444 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1446 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1453 vector<int> IsoClus;
1456 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1457 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1462 if(
k==i ||
k==j)
continue;
1465 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1466 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1467 float drcl = clusVect.DeltaR(pairVect);
1470 Iso = Iso + etClusEndCap[
k];
1471 IsoClus.push_back(k);
1538 if(neta > 0) neta -= 1;
1539 if(nphi > 359) nphi=nphi-360;
1545 mdiff=(neta1-neta2);
1557 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
virtual 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_