42 #define TWOPI 6.283185308
77 if ( clusPhiSize_ % 2 == 0 || clusEtaSize_ % 2 == 0)
78 edm::LogError(
"AlCaPi0RecHitsProducerError") <<
"Size of eta/phi for simple clustering should be odd numbers";
222 hNRecHitsEBpi0_ =
dbe_->book1D(
"nRechitsEBpi0",
"#rechits in pi0 collection EB",100,0.,250.);
225 hNRecHitsEEpi0_ =
dbe_->book1D(
"nRechitsEEpi0",
"#rechits in pi0 collection EE",100,0.,250.);
228 hNRecHitsEBeta_ =
dbe_->book1D(
"nRechitsEBeta",
"#rechits in eta collection EB",100,0.,250.);
231 hNRecHitsEEeta_ =
dbe_->book1D(
"nRechitsEEeta",
"#rechits in eta collection EE",100,0.,250.);
246 hMinvPi0EB_ =
dbe_->book1D(
"Pi0InvmassEB",
"Pi0 Invariant Mass in EB",100,0.,0.5);
249 hMinvPi0EE_ =
dbe_->book1D(
"Pi0InvmassEE",
"Pi0 Invariant Mass in EE",100,0.,0.5);
252 hMinvEtaEB_ =
dbe_->book1D(
"EtaInvmassEB",
"Eta Invariant Mass in EB",100,0.,0.85);
255 hMinvEtaEE_ =
dbe_->book1D(
"EtaInvmassEE",
"Eta Invariant Mass in EE",100,0.,0.85);
259 hPt1Pi0EB_ =
dbe_->book1D(
"Pt1Pi0EB",
"Pt 1st most energetic Pi0 photon in EB",100,0.,20.);
262 hPt1Pi0EE_ =
dbe_->book1D(
"Pt1Pi0EE",
"Pt 1st most energetic Pi0 photon in EE",100,0.,20.);
265 hPt1EtaEB_ =
dbe_->book1D(
"Pt1EtaEB",
"Pt 1st most energetic Eta photon in EB",100,0.,20.);
268 hPt1EtaEE_ =
dbe_->book1D(
"Pt1EtaEE",
"Pt 1st most energetic Eta photon in EE",100,0.,20.);
271 hPt2Pi0EB_ =
dbe_->book1D(
"Pt2Pi0EB",
"Pt 2nd most energetic Pi0 photon in EB",100,0.,20.);
274 hPt2Pi0EE_ =
dbe_->book1D(
"Pt2Pi0EE",
"Pt 2nd most energetic Pi0 photon in EE",100,0.,20.);
277 hPt2EtaEB_ =
dbe_->book1D(
"Pt2EtaEB",
"Pt 2nd most energetic Eta photon in EB",100,0.,20.);
280 hPt2EtaEE_ =
dbe_->book1D(
"Pt2EtaEE",
"Pt 2nd most energetic Eta photon in EE",100,0.,20.);
284 hPtPi0EB_ =
dbe_->book1D(
"PtPi0EB",
"Pi0 Pt in EB",100,0.,20.);
287 hPtPi0EE_ =
dbe_->book1D(
"PtPi0EE",
"Pi0 Pt in EE",100,0.,20.);
290 hPtEtaEB_ =
dbe_->book1D(
"PtEtaEB",
"Eta Pt in EB",100,0.,20.);
293 hPtEtaEE_ =
dbe_->book1D(
"PtEtaEE",
"Eta Pt in EE",100,0.,20.);
308 hS4S91Pi0EB_ =
dbe_->book1D(
"S4S91Pi0EB",
"S4S9 1st most energetic Pi0 photon in EB",50,0.,1.);
311 hS4S91Pi0EE_ =
dbe_->book1D(
"S4S91Pi0EE",
"S4S9 1st most energetic Pi0 photon in EE",50,0.,1.);
314 hS4S91EtaEB_ =
dbe_->book1D(
"S4S91EtaEB",
"S4S9 1st most energetic Eta photon in EB",50,0.,1.);
317 hS4S91EtaEE_ =
dbe_->book1D(
"S4S91EtaEE",
"S4S9 1st most energetic Eta photon in EE",50,0.,1.);
320 hS4S92Pi0EB_ =
dbe_->book1D(
"S4S92Pi0EB",
"S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
323 hS4S92Pi0EE_ =
dbe_->book1D(
"S4S92Pi0EE",
"S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
326 hS4S92EtaEB_ =
dbe_->book1D(
"S4S92EtaEB",
"S4S9 2nd most energetic Pi0 photon in EB",50,0.,1.);
329 hS4S92EtaEE_ =
dbe_->book1D(
"S4S92EtaEE",
"S4S9 2nd most energetic Pi0 photon in EE",50,0.,1.);
360 std::vector<EcalRecHit> seeds;
363 std::vector<EBDetId> usedXtals;
379 bool GetRecHitsCollectionEBpi0 =
true;
385 GetRecHitsCollectionEBpi0 =
false;
389 bool GetRecHitsCollectionEBeta =
true;
397 GetRecHitsCollectionEBeta =
false;
401 bool GetRecHitsCollectionEEpi0 =
true;
407 GetRecHitsCollectionEEpi0 =
false;
411 bool GetRecHitsCollectionEEeta =
true;
417 GetRecHitsCollectionEEeta =
false;
448 if (rhEBpi0.
isValid() && (rhEBpi0->size() > 0) && GetRecHitsCollectionEBpi0){
453 for(itb=rhEBpi0->begin(); itb!=rhEBpi0->end(); ++itb){
456 double energy = itb->energy();
471 etot+= itb->energy();
487 std::vector<float> eClus;
488 std::vector<float> etClus;
489 std::vector<float> etaClus;
490 std::vector<float> thetaClus;
491 std::vector<float> phiClus;
492 std::vector<EBDetId> max_hit;
494 std::vector< std::vector<EcalRecHit> > RecHitsCluster;
495 std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5;
496 std::vector<float> s4s9Clus;
497 std::vector<float> s9s25Clus;
506 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
507 EBDetId seed_id = itseed->id();
508 std::vector<EBDetId>::const_iterator usedIds;
510 bool seedAlreadyUsed=
false;
511 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
512 if(*usedIds==seed_id){
513 seedAlreadyUsed=
true;
518 if(seedAlreadyUsed)
continue;
520 std::vector< std::pair<DetId, float> > clus_used;
523 vector<EcalRecHit> RecHitsInWindow;
524 vector<EcalRecHit> RecHitsInWindow5x5;
526 double simple_energy = 0;
528 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
531 bool HitAlreadyUsed=
false;
532 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
538 if(HitAlreadyUsed)
continue;
545 usedXtals.push_back(*det);
546 RecHitsInWindow.push_back(
EBRecHits[nn]);
547 clus_used.push_back(std::pair<DetId, float>(*det, 1) );
548 simple_energy = simple_energy +
EBRecHits[nn].energy();
553 if(simple_energy <= 0)
continue;
560 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
566 float et_s = simple_energy *
sin(theta_s);
573 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
575 int seed_ieta = seed_id.
ieta();
576 int seed_iphi = seed_id.
iphi();
583 for(
unsigned int j=0;
j<RecHitsInWindow.size();
j++){
586 int ieta = det.
ieta();
587 int iphi = det.
iphi();
591 float en = RecHitsInWindow[
j].energy();
596 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
597 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
598 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
599 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
601 if(
abs(dx)<=1 &&
abs(dy)<=1) e3x3 += en;
602 if(
abs(dx)<=2 &&
abs(dy)<=2) e5x5 += en;
607 if(e3x3 <= 0)
continue;
609 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
613 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
614 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
618 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
621 if(itdet0 != usedXtals.end())
continue;
629 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
636 if(e5x5 <= 0)
continue;
638 eClus.push_back(simple_energy);
639 etClus.push_back(et_s);
640 etaClus.push_back(clus_pos.eta());
641 thetaClus.push_back(theta_s);
642 phiClus.push_back(clus_pos.phi());
643 s4s9Clus.push_back(s4s9_max);
644 s9s25Clus.push_back(e3x3/e5x5);
645 RecHitsCluster.push_back(RecHitsInWindow);
646 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
663 for(Int_t
i=0 ;
i<nClus ;
i++){
664 for(Int_t
j=
i+1 ;
j<nClus ;
j++){
669 float p0x = etClus[
i] *
cos(phiClus[
i]);
670 float p1x = etClus[
j] *
cos(phiClus[
j]);
671 float p0y = etClus[
i] *
sin(phiClus[i]);
672 float p1y = etClus[
j] *
sin(phiClus[j]);
673 float p0z = eClus[
i] *
cos(thetaClus[i]);
674 float p1z = eClus[
j] *
cos(thetaClus[j]);
677 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
681 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
685 std::vector<int> IsoClus;
688 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
689 for(Int_t
k=0 ;
k<nClus ;
k++){
694 if(
k==i ||
k==j)
continue;
695 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
697 float dretacl = fabs(etaClus[k] - pairVect.Eta());
698 float drcl = ClusVect.DeltaR(pairVect);
702 Iso = Iso + etClus[
k];
703 IsoClus.push_back(k);
745 if (rhEBeta.
isValid() && (rhEBeta->size() > 0) && GetRecHitsCollectionEBeta){
750 for(itb=rhEBeta->begin(); itb!=rhEBeta->end(); ++itb){
753 double energy = itb->energy();
768 etot+= itb->energy();
784 std::vector<float> eClus;
785 std::vector<float> etClus;
786 std::vector<float> etaClus;
787 std::vector<float> thetaClus;
788 std::vector<float> phiClus;
789 std::vector<EBDetId> max_hit;
791 std::vector< std::vector<EcalRecHit> > RecHitsCluster;
792 std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5;
793 std::vector<float> s4s9Clus;
794 std::vector<float> s9s25Clus;
802 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
803 EBDetId seed_id = itseed->id();
804 std::vector<EBDetId>::const_iterator usedIds;
806 bool seedAlreadyUsed=
false;
807 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
808 if(*usedIds==seed_id){
809 seedAlreadyUsed=
true;
814 if(seedAlreadyUsed)
continue;
816 std::vector< std::pair<DetId, float> > clus_used;
819 vector<EcalRecHit> RecHitsInWindow;
820 vector<EcalRecHit> RecHitsInWindow5x5;
822 double simple_energy = 0;
824 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
827 bool HitAlreadyUsed=
false;
828 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
834 if(HitAlreadyUsed)
continue;
841 usedXtals.push_back(*det);
842 RecHitsInWindow.push_back(
EBRecHits[nn]);
843 clus_used.push_back(std::pair<DetId, float>(*det, 1));
844 simple_energy = simple_energy +
EBRecHits[nn].energy();
849 if(simple_energy <= 0)
continue;
856 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
862 float et_s = simple_energy *
sin(theta_s);
869 for(
int i=0;
i<4;
i++)s4s9_tmp[
i]= 0;
871 int seed_ieta = seed_id.
ieta();
872 int seed_iphi = seed_id.
iphi();
879 for(
unsigned int j=0;
j<RecHitsInWindow.size();
j++){
882 int ieta = det.
ieta();
883 int iphi = det.
iphi();
887 float en = RecHitsInWindow[
j].energy();
892 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
893 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
894 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
895 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
897 if(
abs(dx)<=1 &&
abs(dy)<=1) e3x3 += en;
898 if(
abs(dx)<=2 &&
abs(dy)<=2) e5x5 += en;
903 if(e3x3 <= 0)
continue;
905 float s4s9_max = *max_element( s4s9_tmp,s4s9_tmp+4)/e3x3;
909 std::vector<DetId> clus_v5x5 = topology_p->
getWindow(seed_id,5,5);
910 for( std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++){
914 std::vector<EBDetId>::iterator itdet0 =
find(usedXtals.begin(),usedXtals.end(),det);
917 if(itdet0 != usedXtals.end())
continue;
925 RecHitsInWindow5x5.push_back(
EBRecHits[nn]);
932 if(e5x5 <= 0)
continue;
934 eClus.push_back(simple_energy);
935 etClus.push_back(et_s);
936 etaClus.push_back(clus_pos.eta());
937 thetaClus.push_back(theta_s);
938 phiClus.push_back(clus_pos.phi());
939 s4s9Clus.push_back(s4s9_max);
940 s9s25Clus.push_back(e3x3/e5x5);
941 RecHitsCluster.push_back(RecHitsInWindow);
942 RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
959 for(Int_t
i=0 ;
i<nClus ;
i++){
960 for(Int_t
j=
i+1 ;
j<nClus ;
j++){
965 float p0x = etClus[
i] *
cos(phiClus[
i]);
966 float p1x = etClus[
j] *
cos(phiClus[
j]);
967 float p0y = etClus[
i] *
sin(phiClus[i]);
968 float p1y = etClus[
j] *
sin(phiClus[j]);
969 float p0z = eClus[
i] *
cos(thetaClus[i]);
970 float p1z = eClus[
j] *
cos(thetaClus[j]);
973 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
977 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
981 std::vector<int> IsoClus;
984 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
985 for(Int_t
k=0 ;
k<nClus ;
k++){
990 if(
k==i ||
k==j)
continue;
991 TVector3 ClusVect = TVector3(etClus[
k] *
cos(phiClus[
k]), etClus[k] *
sin(phiClus[k]) , eClus[k] *
cos(thetaClus[k]));
993 float dretacl = fabs(etaClus[k] - pairVect.Eta());
994 float drcl = ClusVect.DeltaR(pairVect);
998 Iso = Iso + etClus[
k];
999 IsoClus.push_back(k);
1050 if (rhEEpi0.
isValid() && (rhEEpi0->size() > 0) && GetRecHitsCollectionEEpi0){
1060 std::vector<EcalRecHit> seedsEndCap;
1061 seedsEndCap.clear();
1063 std::vector<EEDetId> usedXtalsEndCap;
1064 usedXtalsEndCap.clear();
1069 for (ite=rhEEpi0->begin(); ite!=rhEEpi0->end(); ite++) {
1070 double energy = ite->energy();
1086 etot+= ite->energy();
1096 std::vector<float> eClusEndCap;
1097 std::vector<float> etClusEndCap;
1098 std::vector<float> etaClusEndCap;
1099 std::vector<float> thetaClusEndCap;
1100 std::vector<float> phiClusEndCap;
1101 std::vector< std::vector<EcalRecHit> > RecHitsClusterEndCap;
1102 std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1103 std::vector<float> s4s9ClusEndCap;
1104 std::vector<float> s9s25ClusEndCap;
1113 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1114 EEDetId seed_id = itseed->id();
1115 std::vector<EEDetId>::const_iterator usedIds;
1117 bool seedAlreadyUsed=
false;
1118 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1119 if(*usedIds==seed_id){
1120 seedAlreadyUsed=
true;
1125 if(seedAlreadyUsed)
continue;
1127 std::vector< std::pair<DetId, float> > clus_used;
1129 vector<EcalRecHit> RecHitsInWindow;
1130 vector<EcalRecHit> RecHitsInWindow5x5;
1132 float simple_energy = 0;
1134 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1137 bool HitAlreadyUsed=
false;
1138 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1140 HitAlreadyUsed=
true;
1145 if(HitAlreadyUsed)
continue;
1152 usedXtalsEndCap.push_back(*det);
1153 RecHitsInWindow.push_back(
EERecHits[nn]);
1154 clus_used.push_back(std::pair<DetId, float>(*det, 1));
1155 simple_energy = simple_energy +
EERecHits[nn].energy();
1160 if( simple_energy <= 0)
continue;
1165 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1166 float et_s = simple_energy *
sin(theta_s);
1175 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1177 int ixSeed = seed_id.
ix();
1178 int iySeed = seed_id.
iy();
1182 for(
unsigned int j=0;
j<RecHitsInWindow.size();
j++){
1184 int dx = ixSeed - det_this.
ix();
1185 int dy = iySeed - det_this.
iy();
1187 float en = RecHitsInWindow[
j].energy();
1189 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1190 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1191 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1192 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1194 if(
abs(dx)<=1 &&
abs(dy)<=1) e3x3 += en;
1195 if(
abs(dx)<=2 &&
abs(dy)<=2) e5x5 += en;
1200 if(e3x3 <= 0)
continue;
1202 eClusEndCap.push_back(simple_energy);
1203 etClusEndCap.push_back(et_s);
1204 etaClusEndCap.push_back(clus_pos.eta());
1205 thetaClusEndCap.push_back(theta_s);
1206 phiClusEndCap.push_back(clus_pos.phi());
1207 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1208 s9s25ClusEndCap.push_back(e3x3/e5x5);
1209 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1210 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1228 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1229 for(Int_t
j=
i+1 ;
j<nClusEndCap ;
j++){
1233 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1234 float p1x = etClusEndCap[
j] *
cos(phiClusEndCap[
j]);
1235 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1236 float p1y = etClusEndCap[
j] *
sin(phiClusEndCap[j]);
1237 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1238 float p1z = eClusEndCap[
j] *
cos(thetaClusEndCap[j]);
1241 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1243 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1250 std::vector<int> IsoClus;
1253 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1254 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1259 if(
k==i ||
k==j)
continue;
1262 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1263 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1264 float drcl = clusVect.DeltaR(pairVect);
1267 Iso = Iso + etClusEndCap[
k];
1268 IsoClus.push_back(k);
1304 if (rhEEeta.
isValid() && (rhEEeta->size() > 0) && GetRecHitsCollectionEEeta){
1313 std::vector<EcalRecHit> seedsEndCap;
1314 seedsEndCap.clear();
1316 std::vector<EEDetId> usedXtalsEndCap;
1317 usedXtalsEndCap.clear();
1322 for (ite=rhEEeta->begin(); ite!=rhEEeta->end(); ite++) {
1323 double energy = ite->energy();
1339 etot+= ite->energy();
1349 std::vector<float> eClusEndCap;
1350 std::vector<float> etClusEndCap;
1351 std::vector<float> etaClusEndCap;
1352 std::vector<float> thetaClusEndCap;
1353 std::vector<float> phiClusEndCap;
1354 std::vector< std::vector<EcalRecHit> > RecHitsClusterEndCap;
1355 std::vector< std::vector<EcalRecHit> > RecHitsCluster5x5EndCap;
1356 std::vector<float> s4s9ClusEndCap;
1357 std::vector<float> s9s25ClusEndCap;
1366 for (std::vector<EcalRecHit>::iterator itseed=seedsEndCap.begin(); itseed!=seedsEndCap.end(); itseed++) {
1367 EEDetId seed_id = itseed->id();
1368 std::vector<EEDetId>::const_iterator usedIds;
1370 bool seedAlreadyUsed=
false;
1371 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1372 if(*usedIds==seed_id){
1373 seedAlreadyUsed=
true;
1378 if(seedAlreadyUsed)
continue;
1380 std::vector< std::pair<DetId, float> > clus_used;
1382 vector<EcalRecHit> RecHitsInWindow;
1383 vector<EcalRecHit> RecHitsInWindow5x5;
1385 float simple_energy = 0;
1387 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
1390 bool HitAlreadyUsed=
false;
1391 for(usedIds=usedXtalsEndCap.begin(); usedIds!=usedXtalsEndCap.end(); usedIds++){
1393 HitAlreadyUsed=
true;
1398 if(HitAlreadyUsed)
continue;
1405 usedXtalsEndCap.push_back(*det);
1406 RecHitsInWindow.push_back(
EERecHits[nn]);
1407 clus_used.push_back(std::pair<DetId, float>(*det, 1));
1408 simple_energy = simple_energy +
EERecHits[nn].energy();
1413 if( simple_energy <= 0)
continue;
1418 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
1419 float et_s = simple_energy *
sin(theta_s);
1428 for(
int i=0;
i<4;
i++) s4s9_tmp[
i]= 0;
1430 int ixSeed = seed_id.
ix();
1431 int iySeed = seed_id.
iy();
1435 for(
unsigned int j=0;
j<RecHitsInWindow.size();
j++){
1437 int dx = ixSeed - det_this.
ix();
1438 int dy = iySeed - det_this.
iy();
1440 float en = RecHitsInWindow[
j].energy();
1442 if(dx <= 0 && dy <=0) s4s9_tmp[0] += en;
1443 if(dx >= 0 && dy <=0) s4s9_tmp[1] += en;
1444 if(dx <= 0 && dy >=0) s4s9_tmp[2] += en;
1445 if(dx >= 0 && dy >=0) s4s9_tmp[3] += en;
1447 if(
abs(dx)<=1 &&
abs(dy)<=1) e3x3 += en;
1448 if(
abs(dx)<=2 &&
abs(dy)<=2) e5x5 += en;
1453 if(e3x3 <= 0)
continue;
1455 eClusEndCap.push_back(simple_energy);
1456 etClusEndCap.push_back(et_s);
1457 etaClusEndCap.push_back(clus_pos.eta());
1458 thetaClusEndCap.push_back(theta_s);
1459 phiClusEndCap.push_back(clus_pos.phi());
1460 s4s9ClusEndCap.push_back(*max_element( s4s9_tmp,s4s9_tmp+4)/e3x3);
1461 s9s25ClusEndCap.push_back(e3x3/e5x5);
1462 RecHitsClusterEndCap.push_back(RecHitsInWindow);
1463 RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1478 for(Int_t
i=0 ;
i<nClusEndCap ;
i++){
1479 for(Int_t
j=
i+1 ;
j<nClusEndCap ;
j++){
1483 float p0x = etClusEndCap[
i] *
cos(phiClusEndCap[
i]);
1484 float p1x = etClusEndCap[
j] *
cos(phiClusEndCap[
j]);
1485 float p0y = etClusEndCap[
i] *
sin(phiClusEndCap[i]);
1486 float p1y = etClusEndCap[
j] *
sin(phiClusEndCap[j]);
1487 float p0z = eClusEndCap[
i] *
cos(thetaClusEndCap[i]);
1488 float p1z = eClusEndCap[
j] *
cos(thetaClusEndCap[j]);
1491 float pt_pair =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
1493 float m_inv =
sqrt ( (eClusEndCap[i] + eClusEndCap[j])*(eClusEndCap[i] + eClusEndCap[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
1500 std::vector<int> IsoClus;
1503 TVector3 pairVect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
1504 for(Int_t
k=0 ;
k<nClusEndCap ;
k++){
1509 if(
k==i ||
k==j)
continue;
1512 TVector3 clusVect = TVector3(etClusEndCap[
k] *
cos(phiClusEndCap[
k]), etClusEndCap[k] *
sin(phiClusEndCap[k]) , eClusEndCap[k] *
cos(thetaClusEndCap[k]) ) ;
1513 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1514 float drcl = clusVect.DeltaR(pairVect);
1517 Iso = Iso + etClusEndCap[
k];
1518 IsoClus.push_back(k);
1585 if(neta > 0) neta -= 1;
1586 if(nphi > 359) nphi=nphi-360;
1589 if(nphi >359 || nphi <0 || neta< -85 || neta > 84)
1591 std::cout <<
" unexpected fatal error in HLTPi0RecHitsFilter::convxtalid "<< nphi <<
" " << neta <<
" "
1599 mdiff=(neta1-neta2);
1606 if(
abs(nphi1-nphi2) < (360-
abs(nphi1-nphi2))) {
1610 mdiff=360-
abs(nphi1-nphi2);
1611 if(nphi1>nphi2) mdiff=-mdiff;
double selePi0BeltDREndCap_
T getParameter(std::string const &) const
MonitorElement * hMinvPi0EE_
Pi0 invariant mass in EE.
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * hMinvEtaEE_
Eta invariant mass in EE.
double selePtGammaEndCap_
for pi0->gg endcap
MonitorElement * hEventEnergyEEpi0_
Distribution of total event energy EE (pi0)
MonitorElement * hPt1EtaEE_
Pt of the 1st most energetic Eta photon in EE.
void analyze(const edm::Event &e, const edm::EventSetup &c)
double ptMinForIsolation_
MonitorElement * hRechitEnergyEBpi0_
Energy Distribution of rechits EB (pi0)
MonitorElement * hiXDistrEEeta_
Distribution of rechits in ix EE (eta)
MonitorElement * hEventEnergyEBeta_
Distribution of total event energy EB (eta)
std::vector< EEDetId > detIdEERecHits
bool getByToken(EDGetToken token, Handle< PROD > &result) const
MonitorElement * hIsoPi0EB_
Pi0 Iso EB.
Sin< T >::type sin(const T &t)
MonitorElement * hS4S92EtaEE_
S4S9 of the 2nd most energetic eta photon EE.
double seleMinvMinPi0EndCap_
MonitorElement * hMeanRecHitEnergyEEpi0_
Distribution of Mean energy per rechit EE (pi0)
double seleS9S25GammaEta_
MonitorElement * hPt1Pi0EE_
Pt of the 1st most energetic Pi0 photon in EE.
std::vector< EcalRecHit >::const_iterator const_iterator
double selePi0BeltDetaEndCap_
HLTAlCaMonPi0(const edm::ParameterSet &)
MonitorElement * hNRecHitsEEpi0_
Distribution of number of RecHits EE (pi0)
MonitorElement * hNRecHitsEBpi0_
Distribution of number of RecHits EB (pi0)
MonitorElement * hIsoEtaEB_
Eta Iso EB.
MonitorElement * hPt2EtaEE_
Pt of the 2nd most energetic Eta photon in EE.
double clusSeedThrEndCap_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
double seleS4S9GammaEndCap_
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEEpi0Token_
void endRun(const edm::Run &r, const edm::EventSetup &c)
MonitorElement * hS4S92Pi0EE_
S4S9 of the 2nd most energetic pi0 photon EE.
double seleEtaBeltDREndCap_
double seleMinvMaxPi0EndCap_
MonitorElement * hS4S91Pi0EB_
S4S9 of the 1st most energetic pi0 photon.
MonitorElement * hRechitEnergyEEpi0_
Energy Distribution of rechits EE (pi0)
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEEetaToken_
MonitorElement * hPt1EtaEB_
Pt of the 1st most energetic Eta photon in EB.
int iphi() const
get the crystal iphi
MonitorElement * hPt2Pi0EE_
Pt of the 2nd most energetic Pi0 photon in EE.
MonitorElement * hEventEnergyEBpi0_
Distribution of total event energy EB (pi0)
MonitorElement * hMinvEtaEB_
Eta invariant mass in EB.
double ptMinForIsolationEndCap_
MonitorElement * hPt1Pi0EB_
Pt of the 1st most energetic Pi0 photon in EB.
MonitorElement * hPtPi0EE_
Pi0 Pt in EE.
MonitorElement * hiPhiDistrEBpi0_
Distribution of rechits in iPhi (pi0)
int diff_neta_s(int, int)
bool saveToFile_
Write to file.
edm::InputTag productMonitoredEBeta_
double seleEtaBeltDetaEndCap_
Cos< T >::type cos(const T &t)
MonitorElement * hS4S91Pi0EE_
S4S9 of the 1st most energetic pi0 photon EE.
std::vector< EcalRecHit > EERecHits
MonitorElement * hMeanRecHitEnergyEBeta_
Distribution of Mean energy per rechit EB (eta)
double seleMinvMaxEtaEndCap_
double selePtGammaEtaEndCap_
for eta->gg endcap
MonitorElement * hRechitEnergyEEeta_
Energy Distribution of rechits EE (eta)
MonitorElement * hS4S92EtaEB_
S4S9 of the 2nd most energetic eta photon.
Abs< T >::type abs(const T &t)
double seleS4S9GammaEtaEndCap_
MonitorElement * hiPhiDistrEBeta_
Distribution of rechits in iPhi (eta)
PositionCalc posCalculator_
edm::InputTag productMonitoredEBpi0_
object to monitor
std::string folderName_
DQM folder name.
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEBpi0Token_
unsigned int prescaleFactor_
Monitor every prescaleFactor_ events.
int ieta() const
get the crystal ieta
MonitorElement * hMeanRecHitEnergyEBpi0_
Distribution of Mean energy per rechit EB (pi0)
MonitorElement * hPt2Pi0EB_
Pt of the 2nd most energetic Pi0 photon in EB.
EventAuxiliary const & eventAuxiliary() const
std::string fileName_
Output file name if required.
MonitorElement * hNRecHitsEBeta_
Distribution of number of RecHits EB (eta)
MonitorElement * hiXDistrEEpi0_
Distribution of rechits in ix EE (pi0)
MonitorElement * hiYDistrEEpi0_
Distribution of rechits in iy EE (pi0)
MonitorElement * hS4S91EtaEB_
S4S9 of the 1st most energetic eta photon.
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
void convxtalid(int &, int &)
MonitorElement * hPtPi0EB_
Pi0 Pt in EB.
T const * product() const
void beginLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &context)
XYZPointD XYZPoint
point in space with cartesian internal representation
MonitorElement * hiYDistrEEeta_
Distribution of rechits in iy EE (eta)
double seleXtalMinEnergy_
MonitorElement * hRechitEnergyEBeta_
Energy Distribution of rechits EB (eta)
edm::EDGetTokenT< EcalRecHitCollection > productMonitoredEBetaToken_
MonitorElement * hPt2EtaEB_
Pt of the 2nd most energetic Eta photon in EB.
MonitorElement * hPtEtaEB_
Eta Pt in EB.
MonitorElement * hIsoPi0EE_
Pi0 Iso EE.
bool isMonEBpi0_
which subdet will be monitored
MonitorElement * hEventEnergyEEeta_
Distribution of total event energy EE (eta)
MonitorElement * hiEtaDistrEBpi0_
Distribution of rechits in iEta (pi0)
MonitorElement * hMinvPi0EB_
Pi0 invariant mass in EB.
double ptMinForIsolationEtaEndCap_
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=0)
double seleXtalMinEnergyEndCap_
int diff_nphi_s(int, int)
edm::InputTag productMonitoredEEeta_
MonitorElement * hNRecHitsEEeta_
Distribution of number of RecHits EE (eta)
double ptMinForIsolationEta_
std::vector< EBDetId > detIdEBRecHits
MonitorElement * hS4S92Pi0EB_
S4S9 of the 2nd most energetic pi0 photon.
double seleMinvMinEtaEndCap_
double selePtGammaEta_
for eta->gg barrel
void endLuminosityBlock(const edm::LuminosityBlock &lumiSeg, const edm::EventSetup &c)
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
MonitorElement * hIsoEtaEE_
Eta Iso EE.
edm::InputTag productMonitoredEEpi0_
object to monitor
EventNumber_t event() const
double seleS9S25GammaEtaEndCap_
std::vector< EcalRecHit > EBRecHits
MonitorElement * hPtEtaEE_
Eta Pt in EE.
MonitorElement * hMeanRecHitEnergyEEeta_
Distribution of Mean energy per rechit EE (eta)
MonitorElement * hiEtaDistrEBeta_
Distribution of rechits in iEta (eta)
void beginRun(const edm::Run &r, const edm::EventSetup &c)
MonitorElement * hS4S91EtaEE_
S4S9 of the 1st most energetic eta photon EE.