128 std::map<DetId, EcalRecHit> recHitsEB_map;
130 std::vector<EcalRecHit> seeds;
134 vector<EBDetId> usedXtals;
139 static const int MAXCLUS = 2000;
142 vector<float> etClus;
143 vector<float> etaClus;
144 vector<float> phiClus;
145 vector<EBDetId> max_hit;
146 vector< vector<EcalRecHit> > RecHitsCluster;
147 vector<float> s4s9Clus;
150 for(itb=rhEB->begin(); itb!=rhEB->end(); ++itb){
152 double energy = itb->energy();
154 std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
155 recHitsEB_map.insert(map_entry);
161 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
162 EBDetId seed_id = itseed->id();
163 std::vector<EBDetId>::const_iterator usedIds;
165 bool seedAlreadyUsed=
false;
166 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
167 if(*usedIds==seed_id){
168 seedAlreadyUsed=
true;
173 if(seedAlreadyUsed)
continue;
177 std::vector<std::pair<DetId, float> > clus_used;
179 vector<EcalRecHit> RecHitsInWindow;
181 double simple_energy = 0;
183 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
186 bool HitAlreadyUsed=
false;
187 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
193 if(HitAlreadyUsed)
continue;
194 if (recHitsEB_map.find(*det) != recHitsEB_map.end()){
196 std::map<DetId, EcalRecHit>::iterator aHit;
197 aHit = recHitsEB_map.find(*det);
198 usedXtals.push_back(*det);
199 RecHitsInWindow.push_back(aHit->second);
200 clus_used.push_back( std::pair<DetId, float>(*det, 1.) );
201 simple_energy = simple_energy + aHit->second.energy();
206 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
207 float p0x_s = simple_energy *
sin(theta_s) *
cos(clus_pos.phi());
208 float p0y_s = simple_energy *
sin(theta_s) *
sin(clus_pos.phi());
210 float et_s =
sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
214 eClus.push_back(simple_energy);
215 etClus.push_back(et_s);
216 etaClus.push_back(clus_pos.eta());
217 phiClus.push_back(clus_pos.phi());
218 max_hit.push_back(seed_id);
219 RecHitsCluster.push_back(RecHitsInWindow);
223 for(
int i=0;
i<4;
i++)s4s9_[
i]= itseed->energy();
224 for(
unsigned int j=0; j<RecHitsInWindow.size();j++){
226 if((((
EBDetId)RecHitsInWindow[j].
id()).ieta() == seed_id.
ieta()-1 && seed_id.
ieta()!=1 ) || ( seed_id.
ieta()==1 && (((
EBDetId)RecHitsInWindow[j].
id()).ieta() == seed_id.
ieta()-2))){
227 if(((
EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.
iphi()-1 ||((
EBDetId)RecHitsInWindow[j].
id()).iphi()-360 == seed_id.
iphi()-1 ){
228 s4s9_[0]+=RecHitsInWindow[j].energy();
230 if(((
EBDetId)RecHitsInWindow[j].
id()).iphi() == seed_id.
iphi()){
231 s4s9_[0]+=RecHitsInWindow[j].energy();
232 s4s9_[1]+=RecHitsInWindow[j].energy();
234 if(((
EBDetId)RecHitsInWindow[j].
id()).iphi() == seed_id.
iphi()+1 ||((
EBDetId)RecHitsInWindow[j].
id()).iphi()-360 == seed_id.
iphi()+1 ){
235 s4s9_[1]+=RecHitsInWindow[j].energy();
240 if(((
EBDetId)RecHitsInWindow[j].
id()).ieta() == seed_id.
ieta()){
241 if(((
EBDetId)RecHitsInWindow[j].
id()).iphi() == seed_id.
iphi()-1 ||((
EBDetId)RecHitsInWindow[j].
id()).iphi()-360 == seed_id.
iphi()-1 ){
242 s4s9_[0]+=RecHitsInWindow[j].energy();
243 s4s9_[3]+=RecHitsInWindow[j].energy();
245 if(((
EBDetId)RecHitsInWindow[j].
id()).iphi() == seed_id.
iphi()+1 ||((
EBDetId)RecHitsInWindow[j].
id()).iphi()-360 == seed_id.
iphi()+1 ){
246 s4s9_[1]+=RecHitsInWindow[j].energy();
247 s4s9_[2]+=RecHitsInWindow[j].energy();
251 if((((
EBDetId)RecHitsInWindow[j].
id()).ieta() == seed_id.
ieta()+1 && seed_id.
ieta()!=-1 ) || ( seed_id.
ieta()==-1 && (((
EBDetId)RecHitsInWindow[j].
id()).ieta() == seed_id.
ieta()+2))){
252 if(((
EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.
iphi()-1 ||((
EBDetId)RecHitsInWindow[j].
id()).iphi()-360 == seed_id.
iphi()-1 ){
253 s4s9_[3]+=RecHitsInWindow[j].energy();
255 if(((
EBDetId)RecHitsInWindow[j].
id()).iphi() == seed_id.
iphi()){
256 s4s9_[2]+=RecHitsInWindow[j].energy();
257 s4s9_[3]+=RecHitsInWindow[j].energy();
259 if(((
EBDetId)RecHitsInWindow[j].
id()).iphi() == seed_id.
iphi()+1 ||((
EBDetId)RecHitsInWindow[j].
id()).iphi()-360 == seed_id.
iphi()+1 ){
260 s4s9_[2]+=RecHitsInWindow[j].energy();
265 cout<<
" (EBDetId)RecHitsInWindow[j].id()).ieta() "<<((
EBDetId)RecHitsInWindow[j].
id()).ieta()<<
" seed_id.ieta() "<<seed_id.
ieta()<<endl;
266 cout<<
" Problem with S4 calculation "<<endl;
return;
271 s4s9Clus.push_back(*max_element( s4s9_,s4s9_+4)/simple_energy);
275 if (nClus == MAXCLUS)
return;
282 static const int MAXPI0S = 200;
285 vector<EBDetId> scXtals;
288 if (nClus <= 1)
return;
289 for(Int_t
i=0 ;
i<nClus ;
i++){
290 for(Int_t j=
i+1 ; j<nClus ; j++){
293 float theta_0 = 2. * atan(
exp(-etaClus[
i]));
294 float theta_1 = 2. * atan(
exp(-etaClus[j]));
296 float p0x = eClus[
i] *
sin(theta_0) *
cos(phiClus[
i]);
297 float p1x = eClus[j] *
sin(theta_1) *
cos(phiClus[j]);
298 float p0y = eClus[
i] *
sin(theta_0) *
sin(phiClus[
i]);
299 float p1y = eClus[j] *
sin(theta_1) *
sin(phiClus[j]);
300 float p0z = eClus[
i] *
cos(theta_0);
301 float p1z = eClus[j] *
cos(theta_1);
303 float pt_pi0 =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
306 float m_inv =
sqrt ( (eClus[
i] + eClus[j])*(eClus[
i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
313 TVector3 pi0vect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
314 for(Int_t
k=0 ;
k<nClus ;
k++){
315 if(
k==
i ||
k==j)
continue;
316 TVector3 Clusvect = TVector3(eClus[
k] *
sin(2. * atan(
exp(-etaClus[
k]))) *
cos(phiClus[k]), eClus[k] *
sin(2. * atan(
exp(-etaClus[k]))) *
sin(phiClus[k]) , eClus[k] *
cos(2. * atan(
exp(-etaClus[k]))));
317 float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
318 float drclpi0 = Clusvect.DeltaR(pi0vect);
322 Iso = Iso + etClus[
k];
323 IsoClus.push_back(k);
340 if(npi0_s == MAXPI0S)
return;
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
double clusSeedThr_
parameters needed for pizero finding
MonitorElement * hIsoPi0EB_
Sin< T >::type sin(const T &t)
MonitorElement * hPtPi0EB_
MonitorElement * hMinvPi0EB_
std::vector< EcalRecHit >::const_iterator const_iterator
edm::ParameterSet posCalcParameters_
int iphi() const
get the crystal iphi
Cos< T >::type cos(const T &t)
int ieta() const
get the crystal ieta
virtual std::vector< DetId > getWindow(const DetId &id, const int &northSouthSize, const int &eastWestSize) const
XYZPointD XYZPoint
point in space with cartesian internal representation
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=0)
MonitorElement * hPt2Pi0EB_
MonitorElement * hPt1Pi0EB_
double seleXtalMinEnergy_
Public Types inherited from edm::stream::EDAnalyzer< edm::RunSummaryCache< dqmDetails::NoCache >, edm::LuminosityBlockSummaryCache< dqmDetails::NoCache > >
1.8.11
