Reference Manual

Go to the documentation of this file.
00001 #include <functional>
00002 #include <algorithm>
00003 #include <iostream>
00004 #include <sstream>
00005 #include <utility>
00006 #include <string>
00007 #include <vector>
00008 #include <memory>
00009 #include <cmath>
00010 
00011 #include <boost/bind.hpp>
00012 #include <boost/ptr_container/ptr_vector.hpp>
00013 
00014 #include <TH1.h>
00015 
00016 #include <HepMC/GenEvent.h>
00017 #include <HepMC/SimpleVector.h>
00018 
00019 #include "FWCore/Framework/interface/EDAnalyzer.h"
00020 #include "FWCore/Framework/interface/MakerMacros.h"
00021 #include "FWCore/Framework/interface/Event.h"
00022 #include "FWCore/Framework/interface/EventSetup.h"
00023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00024 #include "FWCore/Utilities/interface/InputTag.h"
00025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00026 #include "FWCore/ServiceRegistry/interface/Service.h"
00027 
00028 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00029 
00030 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00031 
00032 #include "GeneratorInterface/LHEInterface/interface/JetInput.h"
00033 #include "GeneratorInterface/LHEInterface/interface/JetClustering.h"
00034 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
00035 
00036 using namespace lhef;
00037 
00038 class LHEAnalyzer : public edm::EDAnalyzer {
00039     public:
00040         explicit LHEAnalyzer(const edm::ParameterSet &params);
00041         virtual ~LHEAnalyzer();
00042 
00043     protected:
00044         virtual void analyze(const edm::Event &event,
00045                              const edm::EventSetup &es);
00046 
00047 
00048     private:
00049         void fillDJRSched(unsigned int min, unsigned int max);
00050 
00051         edm::InputTag                           sourceLabel;
00052         JetInput                                jetInput;
00053 
00054         double                                  defaultDeltaCut;
00055         double                                  defaultPtCut;
00056         double                                  maxEta;
00057         bool                                    useEt;
00058         double                                  ptFraction;
00059 
00060         unsigned int                            binsDelta;
00061         double                                  minDelta;
00062         double                                  maxDelta;
00063         unsigned int                            binsPt;
00064         double                                  minPt;
00065         double                                  maxPt;
00066         unsigned int                            minDJR;
00067         unsigned int                            maxDJR;
00068 
00069         boost::ptr_vector<JetClustering>        deltaClustering;
00070         std::auto_ptr<JetClustering>            ptClustering;
00071         std::vector<unsigned int>               djrSched;
00072 
00073         std::vector<TH1*>                       histoDelta;
00074         std::vector<TH1*>                       histoPt;
00075 };
00076 
00077 LHEAnalyzer::LHEAnalyzer(const edm::ParameterSet &params) :
00078         sourceLabel(params.getParameter<edm::InputTag>("src")),
00079         jetInput(params.getParameter<edm::ParameterSet>("jetInput")),
00080         defaultDeltaCut(params.getParameter<double>("defaultDeltaCut")),
00081         defaultPtCut(params.getParameter<double>("defaultPtCut")),
00082         maxEta(params.getParameter<double>("maxEta")),
00083         useEt(params.getParameter<bool>("useEt")),
00084         ptFraction(params.getUntrackedParameter<double>("ptFraction", 0.75)),
00085         binsDelta(params.getParameter<unsigned int>("binsDelta")),
00086         minDelta(params.getParameter<double>("minDelta")),
00087         maxDelta(params.getParameter<double>("maxDelta")),
00088         binsPt(params.getParameter<unsigned int>("binsPt")),
00089         minPt(params.getParameter<double>("minPt")),
00090         maxPt(params.getParameter<double>("maxPt")),
00091         minDJR(params.getParameter<unsigned int>("minDJR")),
00092         maxDJR(params.getParameter<unsigned int>("maxDJR"))
00093 {
00094         edm::ParameterSet jetClusPSet =
00095                 params.getParameter<edm::ParameterSet>("jetClustering");
00096 
00097         for(unsigned int i = 0; i <= binsDelta; i++) {
00098                 double deltaCut =
00099                         minDelta + (maxDelta - minDelta) * i / binsDelta;
00100                 jetClusPSet.addParameter("coneRadius", deltaCut);
00101                 edm::ParameterSet tmp;
00102                 tmp.addParameter("algorithm", jetClusPSet);
00103                 deltaClustering.push_back(
00104                         new JetClustering(tmp, defaultPtCut * ptFraction));
00105         }
00106 
00107         jetClusPSet.addParameter("coneRadius", defaultDeltaCut);
00108         edm::ParameterSet tmp;
00109         tmp.addParameter("algorithm", jetClusPSet);
00110         ptClustering.reset(new JetClustering(tmp, minPt * ptFraction));
00111 
00112         fillDJRSched(minDJR <= 0 ? 1 : minDJR, maxDJR - 1);
00113 
00114         edm::Service<TFileService> fs;
00115         for(unsigned int i = minDJR; i < maxDJR; i++) {
00116                 std::ostringstream ss, ss2;
00117                 ss << (i + 1) << "#rightarrow" << i << " jets";
00118                 ss2 << i;
00119                 TH1 *h = fs->make<TH1D>(("delta" + ss2.str()).c_str(),
00120                                         ("DJR " + ss.str()).c_str(),
00121                                         binsDelta, minDelta, maxDelta);
00122                 h->SetXTitle("p_{T} [GeV/c^2]");
00123                 h->SetYTitle("#delta#sigma [mb]");
00124 
00125                 if (i == 0) {
00126                         h->Delete();
00127                         h = 0;
00128                 }
00129                 histoDelta.push_back(h);
00130 
00131                 std::string what = useEt ? "E" : "p";
00132                 h = fs->make<TH1D>(("pt" + ss2.str()).c_str(),
00133                                    ("DJR " + ss.str()).c_str(), binsPt,
00134                                    std::log10(minPt), std::log10(maxPt));
00135                 h->SetXTitle(("log_{10}(" + what +
00136                                         "_{T} [GeV/c^{2}])").c_str());
00137                 h->SetYTitle("#delta#sigma [mb]");
00138 
00139                 histoPt.push_back(h);
00140         }
00141 }
00142 
00143 LHEAnalyzer::~LHEAnalyzer()
00144 {
00145 }
00146 
00147 void LHEAnalyzer::fillDJRSched(unsigned int min, unsigned int max)
00148 {
00149         unsigned int middle = (min + max) / 2;
00150 
00151         djrSched.push_back(middle);
00152 
00153         if (min < middle)
00154                 fillDJRSched(min, middle - 1);
00155         if (middle < max)
00156                 fillDJRSched(middle + 1, max);
00157 }
00158 
00159 void LHEAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es)
00160 {
00161         using boost::bind;
00162         typedef JetClustering::Jet Jet;
00163 
00164         edm::Handle<edm::HepMCProduct> hepmc;
00165         event.getByLabel(sourceLabel, hepmc);
00166 
00167         std::auto_ptr<HepMC::GenEvent> clonedEvent;
00168         const HepMC::GenEvent *genEvent = hepmc->GetEvent();
00169         if (!genEvent->signal_process_vertex()) {
00170                 clonedEvent.reset(new HepMC::GenEvent(*genEvent));
00171                 const HepMC::GenVertex *signalVertex =
00172                         LHEEvent::findSignalVertex(clonedEvent.get());
00173                 clonedEvent->set_signal_process_vertex(
00174                         const_cast<HepMC::GenVertex*>(signalVertex));
00175                 genEvent = clonedEvent.get();
00176         }
00177 
00178         JetInput::ParticleVector particles = jetInput(genEvent);
00179 
00180         std::vector<Jet> ptJets = (*ptClustering)(particles);
00181         std::sort(ptJets.begin(), ptJets.end(),
00182                   bind(std::greater<double>(),
00183                        bind(useEt ? &Jet::et : &Jet::pt, _1),
00184                        bind(useEt ? &Jet::et : &Jet::pt, _2)));
00185 
00186         typedef std::pair<int, int> Pair;
00187         std::vector<Pair> deltaJets(maxDJR - minDJR + 1,
00188                                     Pair(-1, binsDelta + 1));
00189 
00190         for(std::vector<unsigned int>::const_iterator djr = djrSched.begin();
00191             djr != djrSched.end(); ++djr) {
00192 //std::cout << "DJR schedule " << (*djr + 1) << " -> " << *djr << std::endl;
00193                 int result = -1;
00194                 for(;;) {
00195 //for(int i = minDJR; i <= maxDJR; i++)
00196 //std::cout << "+++ " << i << ": (" << deltaJets[i - minDJR].first << ", " << deltaJets[i - minDJR].second << ")" << std::endl;
00197                         int upper = binsDelta + 1;
00198                         for(int i = *djr; i >= (int)minDJR; i--) {
00199                                 if (deltaJets[i - minDJR].second <=
00200                                                         (int)binsDelta) {
00201                                         upper = deltaJets[i - minDJR].second;
00202                                         break;
00203                                 }
00204                         }
00205                         int lower = -1;
00206                         for(int i = *djr + 1; i <= (int)maxDJR; i++) {
00207                                 if (deltaJets[i - minDJR].first >= 0) {
00208                                         lower = deltaJets[i - minDJR].first;
00209                                         break;
00210                                 }
00211                         }
00212 //std::cout << "\t" << lower << " < " << upper << std::endl;
00213 
00214                         result = (lower + upper + 2) / 2 - 1;
00215                         if (result == lower)
00216                                 break;
00217                         else if (result < lower) {
00218                                 result = -1;
00219                                 break;
00220                         }
00221 
00222                         std::vector<Jet> jets =
00223                                 deltaClustering[result](particles);
00224                         unsigned int nJets = 0;
00225                         for(std::vector<Jet>::const_iterator iter =
00226                                 jets.begin(); iter != jets.end(); ++iter)
00227                                 if ((useEt ? iter->et() : iter->pt())
00228                                                         > defaultPtCut)
00229                                         nJets++;
00230 
00231 //std::cout << "\t---(" << *djr << ")--> bin " << result << ": " << nJets << " jets" << std::endl;
00232 
00233                         if (nJets < minDJR)
00234                                 nJets = minDJR;
00235                         else if (nJets > maxDJR)
00236                                 nJets = maxDJR;
00237 
00238                         for(int j = nJets; j >= (int)minDJR; j--) {
00239                                 if (deltaJets[j - minDJR].first < 0 ||
00240                                     result > deltaJets[j - minDJR].first)
00241                                         deltaJets[j - minDJR].first = result;
00242                         }
00243                         for(int j = nJets; j <= (int)maxDJR; j++) {
00244                                 if (deltaJets[j - minDJR].second <
00245                                                         (int)binsDelta ||
00246                                     result < deltaJets[j - minDJR].second)
00247                                         deltaJets[j - minDJR].second = result;
00248                         }
00249                 }
00250 
00251 //std::cout << "final " << *djr << ": " << result << std::endl;
00252                 TH1 *h = histoDelta[*djr - minDJR];
00253                 h->Fill(h->GetBinCenter(result + 1));
00254 
00255                 h = histoPt[*djr - minDJR];
00256                 if (*djr >= ptJets.size())
00257                         h->Fill(-999.0);
00258                 else if (useEt)
00259                         h->Fill(std::log10(ptJets[*djr].et()));
00260                 else
00261                         h->Fill(std::log10(ptJets[*djr].pt()));
00262         }
00263 
00264         if (minDJR <= 0) {
00265                 TH1 *h = histoPt[0];
00266                 if (minDJR >= ptJets.size())
00267                         h->Fill(-999.0);
00268                 else if (useEt)
00269                         h->Fill(std::log10(ptJets[minDJR].et()));
00270                 else
00271                         h->Fill(std::log10(ptJets[minDJR].pt()));
00272         }
00273 }
00274 
00275 DEFINE_FWK_MODULE(LHEAnalyzer);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/GeneratorInterface/LHEInterface/plugins/LHEAnalyzer.cc
