PiZeroAnalyzer.cc

Go to the documentation of this file.

#include <iostream>

#include "DQMOffline/EGamma/plugins/PiZeroAnalyzer.h"

//#define TWOPI 6.283185308
//

using namespace std;

PiZeroAnalyzer::PiZeroAnalyzer(const edm::ParameterSet& pset)
    : caloGeometryToken_{esConsumes()}, caloTopologyToken_{esConsumes()} {
  fName_ = pset.getUntrackedParameter<std::string>("Name");
  prescaleFactor_ = pset.getUntrackedParameter<int>("prescaleFactor", 1);

  barrelEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit> > >(
      pset.getParameter<edm::InputTag>("barrelEcalHits"));
  endcapEcalHits_token_ = consumes<edm::SortedCollection<EcalRecHit, edm::StrictWeakOrdering<EcalRecHit> > >(
      pset.getParameter<edm::InputTag>("endcapEcalHits"));

  nEvt_ = 0;

  // parameters for Pizero finding
  seleXtalMinEnergy_ = pset.getParameter<double>("seleXtalMinEnergy");
  clusSeedThr_ = pset.getParameter<double>("clusSeedThr");
  clusEtaSize_ = pset.getParameter<int>("clusEtaSize");
  clusPhiSize_ = pset.getParameter<int>("clusPhiSize");
  ParameterLogWeighted_ = pset.getParameter<bool>("ParameterLogWeighted");
  ParameterX0_ = pset.getParameter<double>("ParameterX0");
  ParameterT0_barl_ = pset.getParameter<double>("ParameterT0_barl");
  ParameterW0_ = pset.getParameter<double>("ParameterW0");

  selePtGammaOne_ = pset.getParameter<double>("selePtGammaOne");
  selePtGammaTwo_ = pset.getParameter<double>("selePtGammaTwo");
  seleS4S9GammaOne_ = pset.getParameter<double>("seleS4S9GammaOne");
  seleS4S9GammaTwo_ = pset.getParameter<double>("seleS4S9GammaTwo");
  selePtPi0_ = pset.getParameter<double>("selePtPi0");
  selePi0Iso_ = pset.getParameter<double>("selePi0Iso");
  selePi0BeltDR_ = pset.getParameter<double>("selePi0BeltDR");
  selePi0BeltDeta_ = pset.getParameter<double>("selePi0BeltDeta");
  seleMinvMaxPi0_ = pset.getParameter<double>("seleMinvMaxPi0");
  seleMinvMinPi0_ = pset.getParameter<double>("seleMinvMinPi0");

  posCalcParameters_ = pset.getParameter<edm::ParameterSet>("posCalcParameters");
}

PiZeroAnalyzer::~PiZeroAnalyzer() {}

void PiZeroAnalyzer::bookHistograms(DQMStore::IBooker& ibooker,
                                    edm::Run const& /* iRun */,
                                    edm::EventSetup const& /* iSetup */) {
  currentFolder_.str("");
  currentFolder_ << "Egamma/PiZeroAnalyzer/";
  ibooker.setCurrentFolder(currentFolder_.str());

  hMinvPi0EB_ = ibooker.book1D("Pi0InvmassEB", "Pi0 Invariant Mass in EB", 100, 0., 0.5);
  hMinvPi0EB_->setAxisTitle("Inv Mass [GeV] ", 1);

  hPt1Pi0EB_ = ibooker.book1D("Pt1Pi0EB", "Pt 1st most energetic Pi0 photon in EB", 100, 0., 20.);
  hPt1Pi0EB_->setAxisTitle("1st photon Pt [GeV] ", 1);

  hPt2Pi0EB_ = ibooker.book1D("Pt2Pi0EB", "Pt 2nd most energetic Pi0 photon in EB", 100, 0., 20.);
  hPt2Pi0EB_->setAxisTitle("2nd photon Pt [GeV] ", 1);

  hPtPi0EB_ = ibooker.book1D("PtPi0EB", "Pi0 Pt in EB", 100, 0., 20.);
  hPtPi0EB_->setAxisTitle("Pi0 Pt [GeV] ", 1);

  hIsoPi0EB_ = ibooker.book1D("IsoPi0EB", "Pi0 Iso in EB", 50, 0., 1.);
  hIsoPi0EB_->setAxisTitle("Pi0 Iso", 1);
}

void PiZeroAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& esup) {
  using namespace edm;

  if (nEvt_ % prescaleFactor_)
    return;
  nEvt_++;
  LogInfo("PiZeroAnalyzer") << "PiZeroAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_
                            << "\n";

  // Get EcalRecHits
  bool validEcalRecHits = true;
  Handle<EcalRecHitCollection> barrelHitHandle;
  EcalRecHitCollection barrelRecHits;
  e.getByToken(barrelEcalHits_token_, barrelHitHandle);
  if (!barrelHitHandle.isValid()) {
    edm::LogError("PhotonProducer") << "Error! Can't get the product: barrelEcalHits_token_";
    validEcalRecHits = false;
  }

  Handle<EcalRecHitCollection> endcapHitHandle;
  e.getByToken(endcapEcalHits_token_, endcapHitHandle);
  EcalRecHitCollection endcapRecHits;
  if (!endcapHitHandle.isValid()) {
    edm::LogError("PhotonProducer") << "Error! Can't get the product: endcapEcalHits_token_";
    validEcalRecHits = false;
  }

  if (validEcalRecHits)
    makePizero(esup, barrelHitHandle, endcapHitHandle);
}

void PiZeroAnalyzer::makePizero(const edm::EventSetup& es,
                                const edm::Handle<EcalRecHitCollection> rhEB,
                                const edm::Handle<EcalRecHitCollection> rhEE) {
  const EcalRecHitCollection* hitCollection_p = rhEB.product();

  //edm::ESHandle<CaloGeometry> geoHandle;
  //es.get<CaloGeometryRecord>().get(geoHandle);
  auto geoHandle = es.getHandle(caloGeometryToken_);

  //edm::ESHandle<CaloTopology> theCaloTopology;
  //es.get<CaloTopologyRecord>().get(theCaloTopology);
  auto theCaloTopology = es.getHandle(caloTopologyToken_);

  const CaloSubdetectorTopology* topology_p;
  const CaloSubdetectorGeometry* geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
  const CaloSubdetectorGeometry* geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);

  // Parameters for the position calculation:
  PositionCalc posCalculator_ = PositionCalc(posCalcParameters_);
  //
  std::map<DetId, EcalRecHit> recHitsEB_map;
  //
  std::vector<EcalRecHit> seeds;

  seeds.clear();
  //
  vector<EBDetId> usedXtals;
  usedXtals.clear();
  //
  EcalRecHitCollection::const_iterator itb;
  //
  static const int MAXCLUS = 2000;
  int nClus = 0;
  vector<float> eClus;
  vector<float> etClus;
  vector<float> etaClus;
  vector<float> phiClus;
  vector<EBDetId> max_hit;
  vector<vector<EcalRecHit> > RecHitsCluster;
  vector<float> s4s9Clus;

  // find cluster seeds in EB
  for (itb = rhEB->begin(); itb != rhEB->end(); ++itb) {
    EBDetId id(itb->id());
    double energy = itb->energy();
    if (energy > seleXtalMinEnergy_) {
      std::pair<DetId, EcalRecHit> map_entry(itb->id(), *itb);
      recHitsEB_map.insert(map_entry);
    }
    if (energy > clusSeedThr_)
      seeds.push_back(*itb);
  }  // Eb rechits

  sort(seeds.begin(), seeds.end(), [](auto& x, auto& y) { return (x.energy() > y.energy()); });
  for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
    EBDetId seed_id = itseed->id();
    std::vector<EBDetId>::const_iterator usedIds;

    bool seedAlreadyUsed = false;
    for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
      if (*usedIds == seed_id) {
        seedAlreadyUsed = true;
        //cout<< " Seed with energy "<<itseed->energy()<<" was used !"<<endl;
        break;
      }
    }
    if (seedAlreadyUsed)
      continue;
    topology_p = theCaloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel);
    std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
    //std::vector<DetId> clus_used;
    std::vector<std::pair<DetId, float> > clus_used;

    vector<EcalRecHit> RecHitsInWindow;

    double simple_energy = 0;

    for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
      // EBDetId EBdet = *det;
      //      cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi "<<EBdet.iphi()<<endl;
      bool HitAlreadyUsed = false;
      for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
        if (*usedIds == *det) {
          HitAlreadyUsed = true;
          break;
        }
      }
      if (HitAlreadyUsed)
        continue;
      if (recHitsEB_map.find(*det) != recHitsEB_map.end()) {
        //      cout<<" Used det "<< EBdet<<endl;
        std::map<DetId, EcalRecHit>::iterator aHit;
        aHit = recHitsEB_map.find(*det);
        usedXtals.push_back(*det);
        RecHitsInWindow.push_back(aHit->second);
        clus_used.push_back(std::pair<DetId, float>(*det, 1.));
        simple_energy = simple_energy + aHit->second.energy();
      }
    }

    math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
    float theta_s = 2. * atan(exp(-clus_pos.eta()));
    float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
    float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
    //      float p0z_s = simple_energy * cos(theta_s);
    float et_s = sqrt(p0x_s * p0x_s + p0y_s * p0y_s);

    //cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<" "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;

    eClus.push_back(simple_energy);
    etClus.push_back(et_s);
    etaClus.push_back(clus_pos.eta());
    phiClus.push_back(clus_pos.phi());
    max_hit.push_back(seed_id);
    RecHitsCluster.push_back(RecHitsInWindow);
    //Compute S4/S9 variable
    //We are not sure to have 9 RecHits so need to check eta and phi:
    float s4s9_[4];
    for (int i = 0; i < 4; i++)
      s4s9_[i] = itseed->energy();
    for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
      //cout << " Simple cluster rh, ieta, iphi : "<<((EBDetId)RecHitsInWindow[j].id()).ieta()<<" "<<((EBDetId)RecHitsInWindow[j].id()).iphi()<<endl;
      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))) {
        if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
            ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
          s4s9_[0] += RecHitsInWindow[j].energy();
        } else {
          if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()) {
            s4s9_[0] += RecHitsInWindow[j].energy();
            s4s9_[1] += RecHitsInWindow[j].energy();
          } else {
            if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
                ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
              s4s9_[1] += RecHitsInWindow[j].energy();
            }
          }
        }
      } else {
        if (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()) {
          if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
              ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
            s4s9_[0] += RecHitsInWindow[j].energy();
            s4s9_[3] += RecHitsInWindow[j].energy();
          } else {
            if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
                ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
              s4s9_[1] += RecHitsInWindow[j].energy();
              s4s9_[2] += RecHitsInWindow[j].energy();
            }
          }
        } else {
          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))) {
            if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
                ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
              s4s9_[3] += RecHitsInWindow[j].energy();
            } else {
              if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()) {
                s4s9_[2] += RecHitsInWindow[j].energy();
                s4s9_[3] += RecHitsInWindow[j].energy();
              } else {
                if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
                    ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
                  s4s9_[2] += RecHitsInWindow[j].energy();
                }
              }
            }
          } else {
            cout << " (EBDetId)RecHitsInWindow[j].id()).ieta() " << ((EBDetId)RecHitsInWindow[j].id()).ieta()
                 << " seed_id.ieta() " << seed_id.ieta() << endl;
            cout << " Problem with S4 calculation " << endl;
            return;
          }
        }
      }
    }
    s4s9Clus.push_back(*max_element(s4s9_, s4s9_ + 4) / simple_energy);
    //    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;
    //    cout<<" Max "<<*max_element( s4s9_,s4s9_+4)/simple_energy<<endl;
    nClus++;
    if (nClus == MAXCLUS)
      return;
  }  //  End loop over seed clusters

  // cout<< " Pi0 clusters: "<<nClus<<endl;

  // Selection, based on Simple clustering
  //pi0 candidates
  static const int MAXPI0S = 200;
  int npi0_s = 0;

  vector<EBDetId> scXtals;
  scXtals.clear();

  if (nClus <= 1)
    return;
  for (Int_t i = 0; i < nClus; i++) {
    for (Int_t j = i + 1; j < nClus; j++) {
      //      cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"  etClus[j] "<<etClus[j]<<endl;
      if (etClus[i] > selePtGammaOne_ && etClus[j] > selePtGammaTwo_ && s4s9Clus[i] > seleS4S9GammaOne_ &&
          s4s9Clus[j] > seleS4S9GammaTwo_) {
        float theta_0 = 2. * atan(exp(-etaClus[i]));
        float theta_1 = 2. * atan(exp(-etaClus[j]));

        float p0x = eClus[i] * sin(theta_0) * cos(phiClus[i]);
        float p1x = eClus[j] * sin(theta_1) * cos(phiClus[j]);
        float p0y = eClus[i] * sin(theta_0) * sin(phiClus[i]);
        float p1y = eClus[j] * sin(theta_1) * sin(phiClus[j]);
        float p0z = eClus[i] * cos(theta_0);
        float p1z = eClus[j] * cos(theta_1);

        float pt_pi0 = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
        //      cout<<" pt_pi0 "<<pt_pi0<<endl;
        if (pt_pi0 < selePtPi0_)
          continue;
        float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
                           (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
        if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
          //New Loop on cluster to measure isolation:
          vector<int> IsoClus;
          IsoClus.clear();
          float Iso = 0;
          TVector3 pi0vect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
          for (Int_t k = 0; k < nClus; k++) {
            if (k == i || k == j)
              continue;
            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]))));
            float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
            float drclpi0 = Clusvect.DeltaR(pi0vect);

            if ((drclpi0 < selePi0BeltDR_) && (dretaclpi0 < selePi0BeltDeta_)) {
              Iso = Iso + etClus[k];
              IsoClus.push_back(k);
            }
          }

          if (Iso / pt_pi0 < selePi0Iso_) {
            hMinvPi0EB_->Fill(m_inv);
            hPt1Pi0EB_->Fill(etClus[i]);
            hPt2Pi0EB_->Fill(etClus[j]);
            hPtPi0EB_->Fill(pt_pi0);
            hIsoPi0EB_->Fill(Iso / pt_pi0);

            npi0_s++;
          }

          if (npi0_s == MAXPI0S)
            return;
        }
      }
    }
  }
}