35 barrelEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >(pset.
getParameter<
edm::InputTag>(
"barrelEcalHits"));
36 endcapEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > >(pset.
getParameter<
edm::InputTag>(
"endcapEcalHits"));
46 seleXtalMinEnergy_ = pset.
getParameter<
double> (
"seleXtalMinEnergy");
47 clusSeedThr_ = pset.
getParameter<
double> (
"clusSeedThr");
50 ParameterLogWeighted_ = pset.
getParameter<
bool> (
"ParameterLogWeighted");
51 ParameterX0_ = pset.
getParameter<
double> (
"ParameterX0");
52 ParameterT0_barl_ = pset.
getParameter<
double> (
"ParameterT0_barl");
53 ParameterW0_ = pset.
getParameter<
double> (
"ParameterW0");
55 selePtGammaOne_ = pset.
getParameter<
double> (
"selePtGammaOne");
56 selePtGammaTwo_ = pset.
getParameter<
double> (
"selePtGammaTwo");
57 seleS4S9GammaOne_ = pset.
getParameter<
double> (
"seleS4S9GammaOne");
58 seleS4S9GammaTwo_ = pset.
getParameter<
double> (
"seleS4S9GammaTwo");
61 selePi0BeltDR_ = pset.
getParameter<
double> (
"selePi0BeltDR");
62 selePi0BeltDeta_ = pset.
getParameter<
double> (
"selePi0BeltDeta");
63 seleMinvMaxPi0_ = pset.
getParameter<
double> (
"seleMinvMaxPi0");
64 seleMinvMinPi0_ = pset.
getParameter<
double> (
"seleMinvMinPi0");
94 if (verbosity_ > 0 ) {
111 currentFolder_.str(
"");
112 currentFolder_ <<
"Egamma/PiZeroAnalyzer/";
118 hMinvPi0EB_ =
dbe_->
book1D(
"Pi0InvmassEB",
"Pi0 Invariant Mass in EB",100,0.,0.5);
121 hPt1Pi0EB_ =
dbe_->
book1D(
"Pt1Pi0EB",
"Pt 1st most energetic Pi0 photon in EB",100,0.,20.);
124 hPt2Pi0EB_ =
dbe_->
book1D(
"Pt2Pi0EB",
"Pt 2nd most energetic Pi0 photon in EB",100,0.,20.);
127 hPtPi0EB_ =
dbe_->
book1D(
"PtPi0EB",
"Pi0 Pt in EB",100,0.,20.);
130 hIsoPi0EB_ =
dbe_->
book1D(
"IsoPi0EB",
"Pi0 Iso in EB",50,0.,1.);
148 if (nEvt_% prescaleFactor_ )
return;
150 LogInfo(
"PiZeroAnalyzer") <<
"PiZeroAnalyzer Analyzing event number: " << e.
id() <<
" Global Counter " << nEvt_ <<
"\n";
154 bool validEcalRecHits=
true;
157 e.
getByToken(barrelEcalHits_token_, barrelHitHandle);
158 if (!barrelHitHandle.isValid()) {
159 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product: barrelEcalHits_token_";
160 validEcalRecHits=
false;
164 e.
getByToken(endcapEcalHits_token_, endcapHitHandle);
166 if (!endcapHitHandle.isValid()) {
167 edm::LogError(
"PhotonProducer") <<
"Error! Can't get the product: endcapEcalHits_token_";
168 validEcalRecHits=
false;
171 if (validEcalRecHits) makePizero(esup, barrelHitHandle, endcapHitHandle);
199 std::map<DetId, EcalRecHit> recHitsEB_map;
201 std::vector<EcalRecHit> seeds;
205 vector<EBDetId> usedXtals;
210 static const int MAXCLUS = 2000;
213 vector<float> etClus;
214 vector<float> etaClus;
215 vector<float> phiClus;
216 vector<EBDetId> max_hit;
217 vector< vector<EcalRecHit> > RecHitsCluster;
218 vector<float> s4s9Clus;
221 for(itb=rhEB->begin(); itb!=rhEB->end(); ++itb){
223 double energy = itb->energy();
224 if (energy > seleXtalMinEnergy_) {
225 std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
226 recHitsEB_map.insert(map_entry);
228 if (energy > clusSeedThr_) seeds.push_back(*itb);
232 for (std::vector<EcalRecHit>::iterator itseed=seeds.begin(); itseed!=seeds.end(); itseed++) {
233 EBDetId seed_id = itseed->id();
234 std::vector<EBDetId>::const_iterator usedIds;
236 bool seedAlreadyUsed=
false;
237 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
238 if(*usedIds==seed_id){
239 seedAlreadyUsed=
true;
244 if(seedAlreadyUsed)
continue;
246 std::vector<DetId> clus_v = topology_p->
getWindow(seed_id,clusEtaSize_,clusPhiSize_);
248 std::vector<std::pair<DetId, float> > clus_used;
250 vector<EcalRecHit> RecHitsInWindow;
252 double simple_energy = 0;
254 for (std::vector<DetId>::iterator det=clus_v.begin(); det!=clus_v.end(); det++) {
257 bool HitAlreadyUsed=
false;
258 for(usedIds=usedXtals.begin(); usedIds!=usedXtals.end(); usedIds++){
264 if(HitAlreadyUsed)
continue;
265 if (recHitsEB_map.find(*det) != recHitsEB_map.end()){
267 std::map<DetId, EcalRecHit>::iterator aHit;
268 aHit = recHitsEB_map.find(*det);
269 usedXtals.push_back(*det);
270 RecHitsInWindow.push_back(aHit->second);
271 clus_used.push_back( std::pair<DetId, float>(*det, 1.) );
272 simple_energy = simple_energy + aHit->second.energy();
276 math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used,hitCollection_p,geometry_p,geometryES_p);
277 float theta_s = 2. * atan(
exp(-clus_pos.eta()));
278 float p0x_s = simple_energy *
sin(theta_s) *
cos(clus_pos.phi());
279 float p0y_s = simple_energy *
sin(theta_s) *
sin(clus_pos.phi());
281 float et_s =
sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
285 eClus.push_back(simple_energy);
286 etClus.push_back(et_s);
287 etaClus.push_back(clus_pos.eta());
288 phiClus.push_back(clus_pos.phi());
289 max_hit.push_back(seed_id);
290 RecHitsCluster.push_back(RecHitsInWindow);
294 for(
int i=0;
i<4;
i++)s4s9_[
i]= itseed->energy();
295 for(
unsigned int j=0;
j<RecHitsInWindow.size();
j++){
297 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))){
298 if(((
EBDetId)RecHitsInWindow[
j].id()).iphi() == seed_id.
iphi()-1 ||((
EBDetId)RecHitsInWindow[
j].
id()).iphi()-360 == seed_id.
iphi()-1 ){
299 s4s9_[0]+=RecHitsInWindow[
j].energy();
301 if(((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()){
302 s4s9_[0]+=RecHitsInWindow[
j].energy();
303 s4s9_[1]+=RecHitsInWindow[
j].energy();
305 if(((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()+1 ||((
EBDetId)RecHitsInWindow[
j].
id()).iphi()-360 == seed_id.
iphi()+1 ){
306 s4s9_[1]+=RecHitsInWindow[
j].energy();
311 if(((
EBDetId)RecHitsInWindow[
j].
id()).ieta() == seed_id.
ieta()){
312 if(((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()-1 ||((
EBDetId)RecHitsInWindow[
j].
id()).iphi()-360 == seed_id.
iphi()-1 ){
313 s4s9_[0]+=RecHitsInWindow[
j].energy();
314 s4s9_[3]+=RecHitsInWindow[
j].energy();
316 if(((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()+1 ||((
EBDetId)RecHitsInWindow[
j].
id()).iphi()-360 == seed_id.
iphi()+1 ){
317 s4s9_[1]+=RecHitsInWindow[
j].energy();
318 s4s9_[2]+=RecHitsInWindow[
j].energy();
322 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))){
323 if(((
EBDetId)RecHitsInWindow[
j].id()).iphi() == seed_id.
iphi()-1 ||((
EBDetId)RecHitsInWindow[
j].
id()).iphi()-360 == seed_id.
iphi()-1 ){
324 s4s9_[3]+=RecHitsInWindow[
j].energy();
326 if(((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()){
327 s4s9_[2]+=RecHitsInWindow[
j].energy();
328 s4s9_[3]+=RecHitsInWindow[
j].energy();
330 if(((
EBDetId)RecHitsInWindow[
j].
id()).iphi() == seed_id.
iphi()+1 ||((
EBDetId)RecHitsInWindow[
j].
id()).iphi()-360 == seed_id.
iphi()+1 ){
331 s4s9_[2]+=RecHitsInWindow[
j].energy();
336 cout<<
" (EBDetId)RecHitsInWindow[j].id()).ieta() "<<((
EBDetId)RecHitsInWindow[
j].
id()).ieta()<<
" seed_id.ieta() "<<seed_id.
ieta()<<endl;
337 cout<<
" Problem with S4 calculation "<<endl;
return;
342 s4s9Clus.push_back(*max_element( s4s9_,s4s9_+4)/simple_energy);
346 if (nClus == MAXCLUS)
return;
353 static const int MAXPI0S = 200;
356 vector<EBDetId> scXtals;
359 if (nClus <= 1)
return;
360 for(Int_t
i=0 ;
i<nClus ;
i++){
361 for(Int_t
j=
i+1 ;
j<nClus ;
j++){
363 if( etClus[
i]>selePtGammaOne_ && etClus[
j]>selePtGammaTwo_ && s4s9Clus[
i]>seleS4S9GammaOne_ && s4s9Clus[
j]>seleS4S9GammaTwo_){
364 float theta_0 = 2. * atan(
exp(-etaClus[
i]));
365 float theta_1 = 2. * atan(
exp(-etaClus[
j]));
367 float p0x = eClus[
i] *
sin(theta_0) *
cos(phiClus[i]);
368 float p1x = eClus[
j] *
sin(theta_1) *
cos(phiClus[j]);
369 float p0y = eClus[
i] *
sin(theta_0) *
sin(phiClus[i]);
370 float p1y = eClus[
j] *
sin(theta_1) *
sin(phiClus[j]);
371 float p0z = eClus[
i] *
cos(theta_0);
372 float p1z = eClus[
j] *
cos(theta_1);
374 float pt_pi0 =
sqrt( (p0x+p1x)*(p0x+p1x) + (p0y+p1y)*(p0y+p1y));
376 if (pt_pi0 < selePtPi0_)
continue;
377 float m_inv =
sqrt ( (eClus[i] + eClus[j])*(eClus[i] + eClus[j]) - (p0x+p1x)*(p0x+p1x) - (p0y+p1y)*(p0y+p1y) - (p0z+p1z)*(p0z+p1z) );
378 if ( (m_inv<seleMinvMaxPi0_) && (m_inv>seleMinvMinPi0_) ){
384 TVector3 pi0vect = TVector3((p0x+p1x), (p0y+p1y), (p0z+p1z));
385 for(Int_t
k=0 ;
k<nClus ;
k++){
386 if(
k==i ||
k==j)
continue;
387 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]))));
388 float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
389 float drclpi0 = Clusvect.DeltaR(pi0vect);
391 if((drclpi0<selePi0BeltDR_) && (dretaclpi0<selePi0BeltDeta_) ){
393 Iso = Iso + etClus[
k];
394 IsoClus.push_back(k);
399 if(Iso/pt_pi0<selePi0Iso_){
401 hMinvPi0EB_->Fill(m_inv);
402 hPt1Pi0EB_->Fill(etClus[i]);
403 hPt2Pi0EB_->Fill(etClus[j]);
404 hPtPi0EB_->Fill(pt_pi0);
405 hIsoPi0EB_->Fill(Iso/pt_pi0);
411 if(npi0_s == MAXPI0S)
return;
426 bool outputMEsInRootFile = parameters_.getParameter<
bool>(
"OutputMEsInRootFile");
428 if(outputMEsInRootFile){
432 edm::LogInfo(
"PiZeroAnalyzer") <<
"Analyzed " << nEvt_ <<
"\n";
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
PiZeroAnalyzer(const edm::ParameterSet &)
void makePizero(const edm::EventSetup &es, const edm::Handle< EcalRecHitCollection > eb, const edm::Handle< EcalRecHitCollection > ee)
Sin< T >::type sin(const T &t)
std::vector< EcalRecHit >::const_iterator const_iterator
int iphi() const
get the crystal iphi
Cos< T >::type cos(const T &t)
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
void setVerbose(unsigned level)
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
virtual void analyze(const edm::Event &, const edm::EventSetup &)
T const * product() const
return(e1-e2)*(e1-e2)+dp *dp
virtual ~PiZeroAnalyzer()
void showDirStructure(void) const
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
void setCurrentFolder(const std::string &fullpath)