Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:44

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer
0004 // Class:      HiFJRhoFlowModulationProducer
0005 
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/Common/interface/View.h"
0008 #include "DataFormats/JetReco/interface/Jet.h"
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0011 #include "DataFormats/Provenance/interface/EventID.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
0018 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
0019 #include "DataFormats/JetReco/interface/JetCollection.h"
0020 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0021 #include "FWCore/Framework/interface/stream/EDProducer.h"
0022 #include "FWCore/Utilities/interface/StreamID.h"
0023 #include "FWCore/Utilities/interface/EDGetToken.h"
0024 #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h"
0025 
0026 #include "TF1.h"
0027 #include "TH1.h"
0028 #include "TMath.h"
0029 #include "TMinuitMinimizer.h"
0030 
0031 #include <algorithm>
0032 #include <cmath>
0033 #include <memory>
0034 #include <string>
0035 #include <utility>
0036 #include <vector>
0037 
0038 namespace {
0039   double lineFunction(double* x, double* par) { return par[0]; }
0040   double flowFunction(double* x, double* par) {
0041     return par[0] * (1. + 2. * (par[1] * std::cos(2. * (x[0] - par[2])) + par[3] * std::cos(3. * (x[0] - par[4]))));
0042   }
0043 };  // namespace
0044 
0045 class HiFJRhoFlowModulationProducer : public edm::stream::EDProducer<> {
0046 public:
0047   explicit HiFJRhoFlowModulationProducer(const edm::ParameterSet&);
0048   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0049 
0050 private:
0051   void produce(edm::Event&, const edm::EventSetup&) override;
0052 
0053   const bool doEvtPlane_;
0054   const bool doFreePlaneFit_;
0055   const bool doJettyExclusion_;
0056   const int evtPlaneLevel_;
0057   const edm::EDGetTokenT<reco::JetView> jetToken_;
0058   const edm::EDGetTokenT<reco::PFCandidateCollection> pfCandsToken_;
0059   const edm::EDGetTokenT<reco::EvtPlaneCollection> evtPlaneToken_;
0060   std::unique_ptr<TF1> lineFit_p_;
0061   std::unique_ptr<TF1> flowFit_p_;
0062 };
0063 HiFJRhoFlowModulationProducer::HiFJRhoFlowModulationProducer(const edm::ParameterSet& iConfig)
0064     : doEvtPlane_(iConfig.getParameter<bool>("doEvtPlane")),
0065       doFreePlaneFit_(iConfig.getParameter<bool>("doFreePlaneFit")),
0066       doJettyExclusion_(iConfig.getParameter<bool>("doJettyExclusion")),
0067       evtPlaneLevel_(iConfig.getParameter<int>("evtPlaneLevel")),
0068       jetToken_(doJettyExclusion_ ? consumes<reco::JetView>(iConfig.getParameter<edm::InputTag>("jetTag"))
0069                                   : edm::EDGetTokenT<reco::JetView>()),
0070       pfCandsToken_(consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"))),
0071       evtPlaneToken_(consumes<reco::EvtPlaneCollection>(iConfig.getParameter<edm::InputTag>("EvtPlane"))) {
0072   produces<std::vector<double>>("rhoFlowFitParams");
0073   TMinuitMinimizer::UseStaticMinuit(false);
0074   lineFit_p_ = std::make_unique<TF1>("lineFit", lineFunction, -TMath::Pi(), TMath::Pi());
0075   flowFit_p_ = std::make_unique<TF1>("flowFit", flowFunction, -TMath::Pi(), TMath::Pi());
0076 }
0077 
0078 // ------------ method called to produce the data  ------------
0079 void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0080   int nPhiBins = 10;
0081   auto const& pfCands = iEvent.get(pfCandsToken_);
0082   static constexpr int kMaxEvtPlane = 29;
0083   std::array<float, kMaxEvtPlane> hiEvtPlane;
0084 
0085   if (doEvtPlane_) {
0086     auto const& evtPlanes = iEvent.get(evtPlaneToken_);
0087     assert(evtPlanes.size() == hi::NumEPNames);
0088     std::transform(evtPlanes.begin(), evtPlanes.end(), hiEvtPlane.begin(), [this](auto const& ePlane) -> float {
0089       return ePlane.angle(evtPlaneLevel_);
0090     });
0091   }
0092 
0093   edm::Handle<reco::JetView> jets;
0094   if (doJettyExclusion_)
0095     iEvent.getByToken(jetToken_, jets);
0096 
0097   int nFill = 0;
0098 
0099   double eventPlane2 = -100;
0100   double eventPlane3 = -100;
0101 
0102   constexpr int nParamVals = 9;
0103   auto rhoFlowFitParamsOut = std::make_unique<std::vector<double>>(nParamVals, 1e-6);
0104 
0105   rhoFlowFitParamsOut->at(0) = 0;
0106   rhoFlowFitParamsOut->at(1) = 0;
0107   rhoFlowFitParamsOut->at(2) = 0;
0108   rhoFlowFitParamsOut->at(3) = 0;
0109   rhoFlowFitParamsOut->at(4) = 0;
0110 
0111   double eventPlane2Cos = 0;
0112   double eventPlane2Sin = 0;
0113 
0114   double eventPlane3Cos = 0;
0115   double eventPlane3Sin = 0;
0116 
0117   std::vector<bool> pfcuts(pfCands.size(), false);
0118   int iCand = -1;
0119   for (auto const& pfCandidate : pfCands) {
0120     iCand++;
0121     if (pfCandidate.particleId() != 1 || std::abs(pfCandidate.eta()) > 1.0 || pfCandidate.pt() < .3 ||
0122         pfCandidate.pt() > 3.) {
0123       continue;
0124     }
0125 
0126     if (doJettyExclusion_) {
0127       bool isGood = true;
0128       for (auto const& jet : *jets) {
0129         if (deltaR2(jet, pfCandidate) < .16) {
0130           isGood = false;
0131           break;
0132         }
0133       }
0134       if (!isGood) {
0135         continue;
0136       }
0137     }
0138 
0139     nFill++;
0140     pfcuts[iCand] = true;
0141 
0142     if (!doEvtPlane_) {
0143       eventPlane2Cos += std::cos(2 * pfCandidate.phi());
0144       eventPlane2Sin += std::sin(2 * pfCandidate.phi());
0145 
0146       eventPlane3Cos += std::cos(3 * pfCandidate.phi());
0147       eventPlane3Sin += std::sin(3 * pfCandidate.phi());
0148     }
0149   }
0150 
0151   if (!doEvtPlane_) {
0152     eventPlane2 = std::atan2(eventPlane2Sin, eventPlane2Cos) / 2.;
0153     eventPlane3 = std::atan2(eventPlane3Sin, eventPlane3Cos) / 3.;
0154   } else {
0155     eventPlane2 = hiEvtPlane[hi::HF2];
0156     eventPlane3 = hiEvtPlane[hi::HF3];
0157   }
0158   int pfcuts_count = 0;
0159   if (nFill >= 100 && eventPlane2 > -99) {
0160     nPhiBins = std::max(10, nFill / 30);
0161 
0162     std::string name = "phiTestIEta4_" + std::to_string(iEvent.id().event()) + "_h";
0163     std::string nameFlat = "phiTestIEta4_Flat_" + std::to_string(iEvent.id().event()) + "_h";
0164     std::unique_ptr<TH1F> phi_h = std::make_unique<TH1F>(name.data(), "", nPhiBins, -TMath::Pi(), TMath::Pi());
0165     phi_h->SetDirectory(nullptr);
0166     for (auto const& pfCandidate : pfCands) {
0167       if (pfcuts.at(pfcuts_count))
0168         phi_h->Fill(pfCandidate.phi());
0169       pfcuts_count++;
0170     }
0171     flowFit_p_->SetParameter(0, 10);
0172     flowFit_p_->SetParameter(1, 0);
0173     flowFit_p_->SetParameter(2, eventPlane2);
0174     flowFit_p_->SetParameter(3, 0);
0175     flowFit_p_->SetParameter(4, eventPlane3);
0176     if (!doFreePlaneFit_) {
0177       flowFit_p_->FixParameter(2, eventPlane2);
0178       flowFit_p_->FixParameter(4, eventPlane3);
0179     }
0180 
0181     lineFit_p_->SetParameter(0, 10);
0182 
0183     phi_h->Fit(flowFit_p_.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi());
0184     phi_h->Fit(lineFit_p_.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi());
0185     rhoFlowFitParamsOut->at(0) = flowFit_p_->GetParameter(0);
0186     rhoFlowFitParamsOut->at(1) = flowFit_p_->GetParameter(1);
0187     rhoFlowFitParamsOut->at(2) = flowFit_p_->GetParameter(2);
0188     rhoFlowFitParamsOut->at(3) = flowFit_p_->GetParameter(3);
0189     rhoFlowFitParamsOut->at(4) = flowFit_p_->GetParameter(4);
0190 
0191     rhoFlowFitParamsOut->at(5) = flowFit_p_->GetChisquare();
0192     rhoFlowFitParamsOut->at(6) = flowFit_p_->GetNDF();
0193 
0194     rhoFlowFitParamsOut->at(7) = lineFit_p_->GetChisquare();
0195     rhoFlowFitParamsOut->at(8) = lineFit_p_->GetNDF();
0196 
0197     phi_h.reset();
0198     pfcuts.clear();
0199   }
0200 
0201   iEvent.put(std::move(rhoFlowFitParamsOut), "rhoFlowFitParams");
0202 }
0203 
0204 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0205 void HiFJRhoFlowModulationProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0206   edm::ParameterSetDescription desc;
0207   desc.add<bool>("doEvtPlane", false);
0208   desc.add<edm::InputTag>("EvtPlane", edm::InputTag("hiEvtPlane"));
0209   desc.add<bool>("doJettyExclusion", false);
0210   desc.add<bool>("doFreePlaneFit", false);
0211   desc.add<bool>("doFlatTest", false);
0212   desc.add<edm::InputTag>("jetTag", edm::InputTag("ak4PFJets"));
0213   desc.add<edm::InputTag>("pfCandSource", edm::InputTag("particleFlow"));
0214   desc.add<int>("evtPlaneLevel", 0);
0215   descriptions.add("hiFJRhoFlowModulationProducer", desc);
0216 }
0217 
0218 DEFINE_FWK_MODULE(HiFJRhoFlowModulationProducer);