Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:22 

0001 #include <iostream>
0002 
0003 #include "DQMOffline/EGamma/plugins/PiZeroAnalyzer.h"
0004 
0005 //#define TWOPI 6.283185308
0006 //
0007 
0008 /** \class PiZeroAnalyzer
0009  **
0010  **
0011  **  $Id: PiZeroAnalyzer
0012  **  authors:
0013  **   Nancy Marinelli, U. of Notre Dame, US
0014  **   Jamie Antonelli, U. of Notre Dame, US
0015  **
0016  ***/
0017 
0018 using namespace std;
0019 
0020 PiZeroAnalyzer::PiZeroAnalyzer(const edm::ParameterSet& pset)
0021     : caloGeometryToken_{esConsumes()}, caloTopologyToken_{esConsumes()} {
0022   fName_ = pset.getUntrackedParameter<std::string>("Name");
0023   prescaleFactor_ = pset.getUntrackedParameter<int>("prescaleFactor", 1);
0024 
0025   barrelEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit> > >(
0026       pset.getParameter<edm::InputTag>("barrelEcalHits"));
0027   endcapEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit> > >(
0028       pset.getParameter<edm::InputTag>("endcapEcalHits"));
0029 
0030   nEvt_ = 0;
0031 
0032   // parameters for Pizero finding
0033   seleXtalMinEnergy_ = pset.getParameter<double>("seleXtalMinEnergy");
0034   clusSeedThr_ = pset.getParameter<double>("clusSeedThr");
0035   clusEtaSize_ = pset.getParameter<int>("clusEtaSize");
0036   clusPhiSize_ = pset.getParameter<int>("clusPhiSize");
0037   ParameterLogWeighted_ = pset.getParameter<bool>("ParameterLogWeighted");
0038   ParameterX0_ = pset.getParameter<double>("ParameterX0");
0039   ParameterT0_barl_ = pset.getParameter<double>("ParameterT0_barl");
0040   ParameterW0_ = pset.getParameter<double>("ParameterW0");
0041 
0042   selePtGammaOne_ = pset.getParameter<double>("selePtGammaOne");
0043   selePtGammaTwo_ = pset.getParameter<double>("selePtGammaTwo");
0044   seleS4S9GammaOne_ = pset.getParameter<double>("seleS4S9GammaOne");
0045   seleS4S9GammaTwo_ = pset.getParameter<double>("seleS4S9GammaTwo");
0046   selePtPi0_ = pset.getParameter<double>("selePtPi0");
0047   selePi0Iso_ = pset.getParameter<double>("selePi0Iso");
0048   selePi0BeltDR_ = pset.getParameter<double>("selePi0BeltDR");
0049   selePi0BeltDeta_ = pset.getParameter<double>("selePi0BeltDeta");
0050   seleMinvMaxPi0_ = pset.getParameter<double>("seleMinvMaxPi0");
0051   seleMinvMinPi0_ = pset.getParameter<double>("seleMinvMinPi0");
0052 
0053   posCalcParameters_ = pset.getParameter<edm::ParameterSet>("posCalcParameters");
0054 }
0055 
0056 PiZeroAnalyzer::~PiZeroAnalyzer() {}
0057 
0058 void PiZeroAnalyzer::bookHistograms(DQMStore::IBooker& ibooker,
0059                                     edm::Run const& /* iRun */,
0060                                     edm::EventSetup const& /* iSetup */) {
0061   currentFolder_.str("");
0062   currentFolder_ << "Egamma/PiZeroAnalyzer/";
0063   ibooker.setCurrentFolder(currentFolder_.str());
0064 
0065   hMinvPi0EB_ = ibooker.book1D("Pi0InvmassEB", "Pi0 Invariant Mass in EB", 100, 0., 0.5);
0066   hMinvPi0EB_->setAxisTitle("Inv Mass [GeV] ", 1);
0067 
0068   hPt1Pi0EB_ = ibooker.book1D("Pt1Pi0EB", "Pt 1st most energetic Pi0 photon in EB", 100, 0., 20.);
0069   hPt1Pi0EB_->setAxisTitle("1st photon Pt [GeV] ", 1);
0070 
0071   hPt2Pi0EB_ = ibooker.book1D("Pt2Pi0EB", "Pt 2nd most energetic Pi0 photon in EB", 100, 0., 20.);
0072   hPt2Pi0EB_->setAxisTitle("2nd photon Pt [GeV] ", 1);
0073 
0074   hPtPi0EB_ = ibooker.book1D("PtPi0EB", "Pi0 Pt in EB", 100, 0., 20.);
0075   hPtPi0EB_->setAxisTitle("Pi0 Pt [GeV] ", 1);
0076 
0077   hIsoPi0EB_ = ibooker.book1D("IsoPi0EB", "Pi0 Iso in EB", 50, 0., 1.);
0078   hIsoPi0EB_->setAxisTitle("Pi0 Iso", 1);
0079 }
0080 
0081 void PiZeroAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& esup) {
0082   using namespace edm;
0083 
0084   if (nEvt_ % prescaleFactor_)
0085     return;
0086   nEvt_++;
0087   LogInfo("PiZeroAnalyzer") << "PiZeroAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_
0088                             << "\n";
0089 
0090   // Get EcalRecHits
0091   bool validEcalRecHits = true;
0092   Handle<EcalRecHitCollection> barrelHitHandle;
0093   EcalRecHitCollection barrelRecHits;
0094   e.getByToken(barrelEcalHits_token_, barrelHitHandle);
0095   if (!barrelHitHandle.isValid()) {
0096     edm::LogError("PhotonProducer") << "Error! Can't get the product: barrelEcalHits_token_";
0097     validEcalRecHits = false;
0098   }
0099 
0100   Handle<EcalRecHitCollection> endcapHitHandle;
0101   e.getByToken(endcapEcalHits_token_, endcapHitHandle);
0102   EcalRecHitCollection endcapRecHits;
0103   if (!endcapHitHandle.isValid()) {
0104     edm::LogError("PhotonProducer") << "Error! Can't get the product: endcapEcalHits_token_";
0105     validEcalRecHits = false;
0106   }
0107 
0108   if (validEcalRecHits)
0109     makePizero(esup, barrelHitHandle, endcapHitHandle);
0110 }
0111 
0112 void PiZeroAnalyzer::makePizero(const edm::EventSetup& es,
0113                                 const edm::Handle<EcalRecHitCollection> rhEB,
0114                                 const edm::Handle<EcalRecHitCollection> rhEE) {
0115   const EcalRecHitCollection* hitCollection_p = rhEB.product();
0116 
0117   //edm::ESHandle<CaloGeometry> geoHandle;
0118   //es.get<CaloGeometryRecord>().get(geoHandle);
0119   auto geoHandle = es.getHandle(caloGeometryToken_);
0120 
0121   //edm::ESHandle<CaloTopology> theCaloTopology;
0122   //es.get<CaloTopologyRecord>().get(theCaloTopology);
0123   auto theCaloTopology = es.getHandle(caloTopologyToken_);
0124 
0125   const CaloSubdetectorTopology* topology_p;
0126   const CaloSubdetectorGeometry* geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0127   const CaloSubdetectorGeometry* geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0128 
0129   // Parameters for the position calculation:
0130   PositionCalc posCalculator_ = PositionCalc(posCalcParameters_);
0131   //
0132   std::map<DetId, EcalRecHit> recHitsEB_map;
0133   //
0134   std::vector<EcalRecHit> seeds;
0135 
0136   seeds.clear();
0137   //
0138   vector<EBDetId> usedXtals;
0139   usedXtals.clear();
0140   //
0141   EcalRecHitCollection::const_iterator itb;
0142   //
0143   static const int MAXCLUS = 2000;
0144   int nClus = 0;
0145   vector<float> eClus;
0146   vector<float> etClus;
0147   vector<float> etaClus;
0148   vector<float> phiClus;
0149   vector<EBDetId> max_hit;
0150   vector<vector<EcalRecHit> > RecHitsCluster;
0151   vector<float> s4s9Clus;
0152 
0153   // find cluster seeds in EB
0154   for (itb = rhEB->begin(); itb != rhEB->end(); ++itb) {
0155     EBDetId id(itb->id());
0156     double energy = itb->energy();
0157     if (energy > seleXtalMinEnergy_) {
0158       std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
0159       recHitsEB_map.insert(map_entry);
0160     }
0161     if (energy > clusSeedThr_)
0162       seeds.push_back(*itb);
0163   }  // Eb rechits
0164 
0165   sort(seeds.begin(), seeds.end(), [](auto& x, auto& y) { return (x.energy() > y.energy()); });
0166   for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
0167     EBDetId seed_id = itseed->id();
0168     std::vector<EBDetId>::const_iterator usedIds;
0169 
0170     bool seedAlreadyUsed = false;
0171     for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
0172       if (*usedIds == seed_id) {
0173         seedAlreadyUsed = true;
0174         //cout<< " Seed with energy "<<itseed->energy()<<" was used !"<<endl;
0175         break;
0176       }
0177     }
0178     if (seedAlreadyUsed)
0179       continue;
0180     topology_p = theCaloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel);
0181     std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
0182     //std::vector<DetId> clus_used;
0183     std::vector<std::pair<DetId, float> > clus_used;
0184 
0185     vector<EcalRecHit> RecHitsInWindow;
0186 
0187     double simple_energy = 0;
0188 
0189     for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
0190       // EBDetId EBdet = *det;
0191       //      cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi "<<EBdet.iphi()<<endl;
0192       bool HitAlreadyUsed = false;
0193       for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
0194         if (*usedIds == *det) {
0195           HitAlreadyUsed = true;
0196           break;
0197         }
0198       }
0199       if (HitAlreadyUsed)
0200         continue;
0201       if (recHitsEB_map.find(*det) != recHitsEB_map.end()) {
0202         //      cout<<" Used det "<< EBdet<<endl;
0203         std::map<DetId, EcalRecHit>::iterator aHit;
0204         aHit = recHitsEB_map.find(*det);
0205         usedXtals.push_back(*det);
0206         RecHitsInWindow.push_back(aHit->second);
0207         clus_used.push_back(std::pair<DetId, float>(*det, 1.));
0208         simple_energy = simple_energy + aHit->second.energy();
0209       }
0210     }
0211 
0212     math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
0213     float theta_s = 2. * atan(exp(-clus_pos.eta()));
0214     float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
0215     float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
0216     //      float p0z_s = simple_energy * cos(theta_s);
0217     float et_s = sqrt(p0x_s * p0x_s + p0y_s * p0y_s);
0218 
0219     //cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<" "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
0220 
0221     eClus.push_back(simple_energy);
0222     etClus.push_back(et_s);
0223     etaClus.push_back(clus_pos.eta());
0224     phiClus.push_back(clus_pos.phi());
0225     max_hit.push_back(seed_id);
0226     RecHitsCluster.push_back(RecHitsInWindow);
0227     //Compute S4/S9 variable
0228     //We are not sure to have 9 RecHits so need to check eta and phi:
0229     float s4s9_[4];
0230     for (int i = 0; i < 4; i++)
0231       s4s9_[i] = itseed->energy();
0232     for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
0233       //cout << " Simple cluster rh, ieta, iphi : "<<((EBDetId)RecHitsInWindow[j].id()).ieta()<<" "<<((EBDetId)RecHitsInWindow[j].id()).iphi()<<endl;
0234       if ((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() - 1 && seed_id.ieta() != 1) ||
0235           (seed_id.ieta() == 1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() - 2))) {
0236         if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
0237             ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
0238           s4s9_[0] += RecHitsInWindow[j].energy();
0239         } else {
0240           if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()) {
0241             s4s9_[0] += RecHitsInWindow[j].energy();
0242             s4s9_[1] += RecHitsInWindow[j].energy();
0243           } else {
0244             if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
0245                 ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
0246               s4s9_[1] += RecHitsInWindow[j].energy();
0247             }
0248           }
0249         }
0250       } else {
0251         if (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()) {
0252           if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
0253               ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
0254             s4s9_[0] += RecHitsInWindow[j].energy();
0255             s4s9_[3] += RecHitsInWindow[j].energy();
0256           } else {
0257             if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
0258                 ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
0259               s4s9_[1] += RecHitsInWindow[j].energy();
0260               s4s9_[2] += RecHitsInWindow[j].energy();
0261             }
0262           }
0263         } else {
0264           if ((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() + 1 && seed_id.ieta() != -1) ||
0265               (seed_id.ieta() == -1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() + 2))) {
0266             if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
0267                 ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
0268               s4s9_[3] += RecHitsInWindow[j].energy();
0269             } else {
0270               if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()) {
0271                 s4s9_[2] += RecHitsInWindow[j].energy();
0272                 s4s9_[3] += RecHitsInWindow[j].energy();
0273               } else {
0274                 if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
0275                     ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
0276                   s4s9_[2] += RecHitsInWindow[j].energy();
0277                 }
0278               }
0279             }
0280           } else {
0281             cout << " (EBDetId)RecHitsInWindow[j].id()).ieta() " << ((EBDetId)RecHitsInWindow[j].id()).ieta()
0282                  << " seed_id.ieta() " << seed_id.ieta() << endl;
0283             cout << " Problem with S4 calculation " << endl;
0284             return;
0285           }
0286         }
0287       }
0288     }
0289     s4s9Clus.push_back(*max_element(s4s9_, s4s9_ + 4) / simple_energy);
0290     //    cout<<" s4s9Clus[0] "<<s4s9_[0]/simple_energy<<" s4s9Clus[1] "<<s4s9_[1]/simple_energy<<" s4s9Clus[2] "<<s4s9_[2]/simple_energy<<" s4s9Clus[3] "<<s4s9_[3]/simple_energy<<endl;
0291     //    cout<<" Max "<<*max_element( s4s9_,s4s9_+4)/simple_energy<<endl;
0292     nClus++;
0293     if (nClus == MAXCLUS)
0294       return;
0295   }  //  End loop over seed clusters
0296 
0297   // cout<< " Pi0 clusters: "<<nClus<<endl;
0298 
0299   // Selection, based on Simple clustering
0300   //pi0 candidates
0301   static const int MAXPI0S = 200;
0302   int npi0_s = 0;
0303 
0304   vector<EBDetId> scXtals;
0305   scXtals.clear();
0306 
0307   if (nClus <= 1)
0308     return;
0309   for (Int_t i = 0; i < nClus; i++) {
0310     for (Int_t j = i + 1; j < nClus; j++) {
0311       //      cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"  etClus[j] "<<etClus[j]<<endl;
0312       if (etClus[i] > selePtGammaOne_ && etClus[j] > selePtGammaTwo_ && s4s9Clus[i] > seleS4S9GammaOne_ &&
0313           s4s9Clus[j] > seleS4S9GammaTwo_) {
0314         float theta_0 = 2. * atan(exp(-etaClus[i]));
0315         float theta_1 = 2. * atan(exp(-etaClus[j]));
0316 
0317         float p0x = eClus[i] * sin(theta_0) * cos(phiClus[i]);
0318         float p1x = eClus[j] * sin(theta_1) * cos(phiClus[j]);
0319         float p0y = eClus[i] * sin(theta_0) * sin(phiClus[i]);
0320         float p1y = eClus[j] * sin(theta_1) * sin(phiClus[j]);
0321         float p0z = eClus[i] * cos(theta_0);
0322         float p1z = eClus[j] * cos(theta_1);
0323 
0324         float pt_pi0 = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
0325         //      cout<<" pt_pi0 "<<pt_pi0<<endl;
0326         if (pt_pi0 < selePtPi0_)
0327           continue;
0328         float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
0329                            (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
0330         if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
0331           //New Loop on cluster to measure isolation:
0332           vector<int> IsoClus;
0333           IsoClus.clear();
0334           float Iso = 0;
0335           TVector3 pi0vect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
0336           for (Int_t k = 0; k < nClus; k++) {
0337             if (k == i || k == j)
0338               continue;
0339             TVector3 Clusvect = TVector3(eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * cos(phiClus[k]),
0340                                          eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * sin(phiClus[k]),
0341                                          eClus[k] * cos(2. * atan(exp(-etaClus[k]))));
0342             float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
0343             float drclpi0 = Clusvect.DeltaR(pi0vect);
0344 
0345             if ((drclpi0 < selePi0BeltDR_) && (dretaclpi0 < selePi0BeltDeta_)) {
0346               Iso = Iso + etClus[k];
0347               IsoClus.push_back(k);
0348             }
0349           }
0350 
0351           if (Iso / pt_pi0 < selePi0Iso_) {
0352             hMinvPi0EB_->Fill(m_inv);
0353             hPt1Pi0EB_->Fill(etClus[i]);
0354             hPt2Pi0EB_->Fill(etClus[j]);
0355             hPtPi0EB_->Fill(pt_pi0);
0356             hIsoPi0EB_->Fill(Iso / pt_pi0);
0357 
0358             npi0_s++;
0359           }
0360 
0361           if (npi0_s == MAXPI0S)
0362             return;
0363         }
0364       }
0365     }
0366   }
0367 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!