FWLite/interface/Scanner.h

0001 #ifndef PhysicsTools_FWLite_Scanner_h
0002 #define PhysicsTools_FWLite_Scanner_h
0003 // these includes are FWLite-safe
0004 #include "DataFormats/FWLite/interface/Handle.h"
0005 #include "DataFormats/FWLite/interface/Event.h"
0006 // these are from ROOT, so they're safe too
0007 #include <TString.h>
0008 #include <TRegexp.h>
0009 #include <TObjString.h>
0010 #include <TObjArray.h>
0011 #include <TDirectory.h>
0012 #include <TEnv.h>
0013 #include <TClass.h>
0014
0015 #if !defined(__CINT__) && !defined(__MAKECINT__)
0016 #include <RooFit.h>
0017 #include <RooArgList.h>
0018 #include <RooDataSet.h>
0019 #include <RooRealVar.h>
0020 #include <RooCategory.h>
0021
0022 #include "PhysicsTools/FWLite/interface/ScannerHelpers.h"
0023 #endif
0024
0025 #include "PhysicsTools/FWLite/interface/EventSelectors.h"
0026
0027 namespace fwlite {
0028
0029   /** \brief fwlite::Scanner<C>, a way to inspect or plots elements of a collection C  by using the StringParser.
0030      *
0031      *  fwlite::Scanner<C>, a way to inspect or plots elements of a collection C  by using the StringParser.
0032      *
0033      *  The collection can be something as easy as std::vector<T>, but also some other fancy EDM collections like
0034      *  RefVector, RefToBaseVector and OwnVector (and probably PtrVector, but it was not tested)
0035      *
0036      *  If you're using something other than std::vector, you must provide the full typename, including all
0037      *  optional template parameters; e.g. you can't have C = edm::RefVector<reco::MuonCollection>, but you need
0038      *  C = edm::RefVector<vector<reco::Muon>,reco::Muon,edm::refhelper::FindUsingAdvance<vector<reco::Muon>,reco::Muon> >
0039      *  In order to figure out what is the correct full name for a collection in an event, open it in ROOT/FWLite,
0040      *  get the branch name including the trailing ".obj" (hint: Events->GetAlias("label")) usually works),
0041      *  and then do Events->GetBranch("xxx.obj")->GetClassName() to get something like edm::Wrapper<X>.
0042      *  then X is what you want to use to create the fwlite::Scanner. Don't use typedefs, they don't work.
0043      *
0044      **/
0045   template <typename Collection>
0046   class Scanner {
0047   public:
0048     /// The type of the Handle to read the Ts from the event. Needed to resolve its Type
0049     typedef fwlite::Handle<Collection> HandleT;
0050
0051     /** Create a Scanner, passing a fwlite Event and the labels (just like you would in 'getByLabel') */
0052     Scanner(fwlite::EventBase *ev, const char *label, const char *instance = "", const char *process = "")
0053         : event_(ev),
0054           label_(label),
0055           instance_(instance),
0056           printFullEventId_(ev->isRealData()),
0057           ignoreExceptions_(false),
0058           exprSep_(":"),
0059           maxEvents_(-1),
0060           maxLinesToPrint_(50) {
0061       objType = helper::Parser::elementType(edm::TypeWithDict(HandleT::TempWrapT::typeInfo()));
0062     }
0063
0064     //------------------------------------------------------------------------------------------------------------------------------------
0065     /** Scan the first nmax entries of the event and print out the values of some expressions.
0066              *
0067              *  The cut is applied to the individual entries. To set Event-wide cuts, use addEventSelector().
0068              *
0069              *  The different expressions are separated by ":", unless changed using setExpressionSeparator.
0070              *  The title of each column is the text of the expression, unless one specifies it differently
0071              *  by using the notation "@label=expression"
0072              *
0073              *  Each row is prefixed by the event id (Run/LS/Event on Data, entry number within the file for MC)
0074              *  and by the index of the object within the collection. The behaviour can be changed through the
0075              *  setPrintFullEventId() method.
0076              *
0077              *  The printing will pause by default every 50 lines (see setMaxLinesToPrint() to change this)
0078              *  Scanning will stop after nmax events.
0079              */
0080     void scan(const char *exprs, const char *cut = "", int nmax = -1) {
0081       helper::ScannerBase scanner(objType);
0082       scanner.setIgnoreExceptions(ignoreExceptions_);
0083
0084       TObjArray *exprArray = TString(exprs).Tokenize(exprSep_);
0085       int rowline = 0;
0086       if (printFullEventId_) {
0087         printf(" : %9s : %4s : %9s : %3s", "RUN", "LUMI", "EVENT", "#IT");
0088         rowline += 3 * 4 + 9 + 4 + 9 + 3 - 1;  // -1 as first char remain blank
0089       } else {
0090         printf(" : %5s : %3s", "EVENT", "#IT");
0091         rowline += 3 + 6 + 3 + 3 - 1;  // -1 as first char remain blank
0092       }
0093       for (int i = 0; i < exprArray->GetEntries(); ++i) {
0094         TString str = ((TObjString *)(*exprArray)[i])->GetString();
0095         std::string lb = str.Data();
0096         std::string ex = str.Data();
0097         if ((ex[0] == '@') && (ex.find('=') != std::string::npos)) {
0098           lb = lb.substr(1, ex.find('=') - 1);
0099           ex = ex.substr(ex.find('=') + 1);
0100         }
0101         scanner.addExpression(ex.c_str());
0102         printf(" : %8s",
0103                (lb.size() > 8 ? lb.substr(lb.size() - 8) : lb)
0104                    .c_str());  // the rightmost part is usually the more interesting one
0105         rowline += 3 + 8;
0106       }
0107       std::cout << " :" << std::endl;
0108       rowline += 2;
0109       delete exprArray;
0110
0111       TString rule('-', rowline);
0112       std::cout << " " << rule << " " << std::endl;
0113
0114       if (strlen(cut))
0115         scanner.setCut(cut);
0116
0117       int iev = 0, line = 0;
0118       for (event_->toBegin(); (iev != nmax) && !event_->atEnd(); ++iev, ++(*event_)) {
0119         if (!selectEvent(*event_))
0120           continue;
0121         handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0122         if (handle_.failedToGet()) {
0123           if (ignoreExceptions_)
0124             continue;
0125         }
0126         const Collection &vals = *handle_;
0127         for (size_t j = 0, n = vals.size(); j < n; ++j) {
0128           if (!scanner.test(&vals[j]))
0129             continue;
0130           if (printFullEventId_) {
0131             const edm::EventAuxiliary &id = event_->eventAuxiliary();
0132             printf(" : %9u : %4u : %9llu : %3lu", id.run(), id.luminosityBlock(), id.event(), (unsigned long)j);
0133           } else {
0134             printf(" : %5d : %3lu", iev, (unsigned long)j);
0135           }
0136           scanner.print(&vals[j]);
0137           std::cout << " :" << std::endl;
0138           if (++line == maxLinesToPrint_) {
0139             line = 0;
0140             if (!wantMore()) {
0141               iev = nmax - 1;  // this is to exit the outer loop
0142               break;           // and this to exit the inner one
0143             }
0144           }
0145         }
0146       }
0147       std::cout << std::endl;
0148     }
0149
0150     //------------------------------------------------------------------------------------------------------------------------------------
0151     /** Count the number of entries that pass a given cut.
0152              *  See setMaxEvents() to specify how many events to loop on when counting.
0153              *  Events can be further selected by using addEventSelector(). */
0154     size_t count(const char *cut) {
0155       helper::ScannerBase scanner(objType);
0156       scanner.setIgnoreExceptions(ignoreExceptions_);
0157
0158       scanner.setCut(cut);
0159
0160       size_t npass = 0;
0161       int iev = 0;
0162       for (event_->toBegin(); !event_->atEnd(); ++(*event_), ++iev) {
0163         if (maxEvents_ > -1 && iev > maxEvents_)
0164           break;
0165         if (!selectEvent(*event_))
0166           continue;
0167         handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0168         const Collection &vals = *handle_;
0169         for (size_t j = 0, n = vals.size(); j < n; ++j) {
0170           if (scanner.test(&vals[j]))
0171             npass++;
0172         }
0173       }
0174       return npass;
0175     }
0176
0177     /** Count the number of events, taking into account setMaxEvents() and the event selectors */
0178     size_t countEvents() {
0179       size_t npass = 0;
0180       int iev = 0;
0181       for (event_->toBegin(); !event_->atEnd(); ++(*event_), ++iev) {
0182         if (maxEvents_ > -1 && iev > maxEvents_)
0183           break;
0184         if (selectEvent(*event_))
0185           npass++;
0186       }
0187       return npass;
0188     }
0189
0190     //------------------------------------------------------------------------------------------------------------------------------------
0191     /** Plot the expression expr for events passing 'cut, into histogram hist.
0192              *  hist is *not* reset before filling it, so it will add to the existing content.
0193              *
0194              *  If "NORM" is specified in the draw options, the output histogram is normalized
0195              *  If "GOFF" is specified in the draw options, the output histogram is not drawn
0196              *
0197              *  See setMaxEvents() to specify how many events to loop on when plotting.
0198              *  Events can be further selected by using addEventSelector(). */
0199     TH1 *draw(const char *expr, const char *cut, TString drawopt, TH1 *hist) {
0200       // prep the machinery
0201       helper::ScannerBase scanner(objType);
0202       scanner.setIgnoreExceptions(ignoreExceptions_);
0203       if (!scanner.addExpression(expr))
0204         return nullptr;
0205       if (strlen(cut))
0206         scanner.setCut(cut);
0207
0208       // check histo
0209       if (hist == nullptr) {
0210         std::cerr << "Method draw(expr, cut, drawopt, hist) cannot be called with null 'hist'. Use the other draw "
0211                      "methods instead."
0212                   << std::endl;
0213         return nullptr;
0214       }
0215
0216       // fill histogram
0217       int iev = 0;
0218       for (event_->toBegin(); !event_->atEnd(); ++(*event_), ++iev) {
0219         if (maxEvents_ > -1 && iev > maxEvents_)
0220           break;
0221         if (!selectEvent(*event_))
0222           continue;
0223         handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0224         const Collection &vals = *handle_;
0225         for (size_t j = 0, n = vals.size(); j < n; ++j) {
0226           scanner.fill1D(&vals[j], hist);
0227         }
0228       }
0229
0230       if (drawopt.Contains("NORM", TString::kIgnoreCase) && (hist->Integral() != 0)) {
0231         hist->Sumw2();
0232         hist->Scale(1.0 / hist->Integral());
0233         // remove the "NORM" because THistPainter doesn't understand it
0234         drawopt(TRegexp("[Nn][Oo][Rr][Mm]")) = "";
0235       }
0236
0237       if (!drawopt.Contains("GOFF", TString::kIgnoreCase))
0238         hist->Draw(drawopt);
0239       return hist;
0240     }
0241
0242     /** Plot the expression 'expr' for events passing 'cut', in a histogram named 'hname'
0243              *  - If htemplate is provided it will be cloned,
0244              *  - otherwise, if "SAME" is among the draw options, it will clone "htemp" (if it's available)
0245              *  - otherwise an automatically binned histogram will be used.
0246              *    in the last case, gEnv->GetValue("Hist.Binning.1D.x", 100) is used for the number of bins
0247              *  See draw(const char *expr, const char *cut, TString drawopt, TH1 *histo) for further documentation */
0248     TH1 *draw(const char *expr,
0249               const char *cut = "",
0250               TString drawopt = "",
0251               const char *hname = "htemp",
0252               const TH1 *htemplate = nullptr) {
0253       TH1 *hist = nullptr;
0254       if (htemplate != nullptr) {
0255         if ((strcmp(hname, "htemp") == 0) && (strcmp(hname, htemplate->GetName()) != 0))
0256           htempDelete();
0257         hist = (TH1 *)htemplate->Clone(hname);
0258       } else if (drawopt.Contains("SAME", TString::kIgnoreCase)) {
0259         hist = getSameH1(hname);
0260       }
0261
0262       // if in the end we found no way to make "hist"
0263       if (hist == nullptr) {
0264         if (strcmp(hname, "htemp") == 0)
0265           htempDelete();
0266         hist = new TH1F(hname, "", gEnv->GetValue("Hist.Binning.1D.x", 100), 0, 0);
0267         hist->SetCanExtend(TH1::kAllAxes);
0268       }
0269       hist->SetTitle((strlen(cut) ? TString(expr) + "{" + cut + "}" : TString(expr)));
0270       hist->GetXaxis()->SetTitle(expr);
0271       return draw(expr, cut, drawopt, hist);
0272     }
0273
0274     /// Make a histogram named hname with nbins from xlow to xhigh, and then call draw().
0275     /// If "SAME" is passed in the draw options, complain and ignore the binning.
0276     TH1 *draw(const char *expr,
0277               int nbins,
0278               double xlow,
0279               double xhigh,
0280               const char *cut = "",
0281               const char *drawopt = "",
0282               const char *hname = "htemp") {
0283       if (TString(drawopt).Contains("SAME", TString::kIgnoreCase)) {
0284         std::cerr << "Binning is ignored when 'SAME' is specified." << std::endl;
0285         TH1 *hsame = getSameH1(hname);
0286         return draw(expr, cut, drawopt, hsame);
0287       }
0288       if (strcmp(hname, "htemp") == 0)
0289         htempDelete();
0290       TH1 *htemp = new TH1F(hname, expr, nbins, xlow, xhigh);
0291       if (strlen(cut))
0292         htemp->SetTitle(TString(expr) + "{" + cut + "}");
0293       htemp->GetXaxis()->SetTitle(expr);
0294       return draw(expr, cut, drawopt, htemp);
0295     }
0296     /// Make a histogram named hname with nbins with boundaries xbins, and then call draw().
0297     /// If "SAME" is passed in the draw options, complain and ignore the binning.
0298     TH1 *draw(const char *expr,
0299               int nbins,
0300               double *xbins,
0301               const char *cut = "",
0302               const char *drawopt = "",
0303               const char *hname = "htemp") {
0304       if (TString(drawopt).Contains("SAME", TString::kIgnoreCase)) {
0305         std::cerr << "Binning is ignored when 'SAME' is specified." << std::endl;
0306         TH1 *hsame = getSameH1(hname);
0307         return draw(expr, cut, drawopt, hsame);
0308       }
0309       if (strcmp(hname, "htemp") == 0)
0310         htempDelete();
0311       TH1 *htemp = new TH1F(hname, expr, nbins, xbins);
0312       if (strlen(cut))
0313         htemp->SetTitle(TString(expr) + "{" + cut + "}");
0314       htemp->GetXaxis()->SetTitle(expr);
0315       return draw(expr, cut, drawopt, htemp);
0316     }
0317
0318     //------------------------------------------------------------------------------------------------------------------------------------
0319     /// Just like draw() except that it uses TProfile. Note that the order is (x,y) while in ROOT it's usually (y,x)!
0320     TProfile *drawProf(TString xexpr, TString yexpr, const char *cut, TString drawopt, TProfile *hist) {
0321       // prep the machinery
0322       helper::ScannerBase scanner(objType);
0323       scanner.setIgnoreExceptions(ignoreExceptions_);
0324       if (!scanner.addExpression(xexpr.Data()))
0325         return nullptr;
0326       if (!scanner.addExpression(yexpr.Data()))
0327         return nullptr;
0328       if (strlen(cut))
0329         scanner.setCut(cut);
0330
0331       // check histo
0332       if (hist == nullptr) {
0333         std::cerr << "Method drawProf(xexpr, yexpr, cut, drawopt, hist) cannot be called with null 'hist'. Use the "
0334                      "other draw methods instead."
0335                   << std::endl;
0336         return nullptr;
0337       }
0338
0339       // fill histogram
0340       int iev = 0;
0341       for (event_->toBegin(); !event_->atEnd(); ++(*event_), ++iev) {
0342         if (maxEvents_ > -1 && iev > maxEvents_)
0343           break;
0344         if (!selectEvent(*event_))
0345           continue;
0346         handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0347         const Collection &vals = *handle_;
0348         for (size_t j = 0, n = vals.size(); j < n; ++j) {
0349           scanner.fillProf(&vals[j], hist);
0350         }
0351       }
0352
0353       if (!strlen(hist->GetTitle()))
0354         hist->SetTitle((strlen(cut) ? yexpr + ":" + xexpr + "{" + cut + "}" : yexpr + ":" + xexpr));
0355       if (!strlen(hist->GetXaxis()->GetTitle()))
0356         hist->GetXaxis()->SetTitle(xexpr);
0357       if (!strlen(hist->GetYaxis()->GetTitle()))
0358         hist->GetYaxis()->SetTitle(yexpr);
0359       if (!TString(drawopt).Contains("GOFF", TString::kIgnoreCase))
0360         hist->Draw(drawopt);
0361       return hist;
0362     }
0363     /// Just like draw() except that it uses TProfile. Note that the order is (x,y) while in ROOT it's usually (y,x)!
0364     TProfile *drawProf(TString xexpr,
0365                        TString yexpr,
0366                        const char *cut = "",
0367                        TString drawopt = "",
0368                        const char *hname = "htemp",
0369                        TProfile *htemplate = nullptr) {
0370       TProfile *hist = nullptr;
0371       if (htemplate != nullptr) {
0372         if ((strcmp(hname, "htemp") == 0) && (strcmp(hname, htemplate->GetName()) != 0))
0373           htempDelete();
0374         hist = (TProfile *)htemplate->Clone(hname);
0375       } else if (drawopt.Contains("SAME", TString::kIgnoreCase)) {
0376         hist = getSameProf(hname);
0377       }
0378
0379       // if in the end we found no way to make "hist"
0380       if (hist == nullptr) {
0381         if (strcmp(hname, "htemp") == 0)
0382           htempDelete();
0383         hist = new TProfile(hname, "", gEnv->GetValue("Hist.Binning.1D.x", 100), 0., 0.);
0384         hist->SetCanExtend(TH1::kAllAxes);
0385       }
0386       return drawProf(xexpr, yexpr, cut, drawopt, hist);
0387     }
0388
0389     /// Just like draw() except that it uses TProfile. Note that the order is (x,y) while in ROOT it's usually (y,x)!
0390     TProfile *drawProf(TString xexpr,
0391                        int bins,
0392                        double xlow,
0393                        double xhigh,
0394                        TString yexpr,
0395                        const char *cut = "",
0396                        const char *drawopt = "",
0397                        const char *hname = "htemp") {
0398       if (TString(drawopt).Contains("SAME", TString::kIgnoreCase)) {
0399         std::cerr << "Binning is ignored when 'SAME' is specified." << std::endl;
0400         TProfile *hsame = getSameProf(hname);
0401         return drawProf(xexpr, yexpr, cut, drawopt, hsame);
0402       }
0403       if (strcmp(hname, "htemp") == 0)
0404         htempDelete();
0405       TProfile *htemp = new TProfile(hname, "", bins, xlow, xhigh);
0406       return drawProf(xexpr, yexpr, cut, drawopt, htemp);
0407     }
0408
0409     //------------------------------------------------------------------------------------------------------------------------------------
0410     /// Just like draw() except that it uses TH2. Note that the order is (x,y) while in ROOT it's usually (y,x)!
0411     TH2 *draw2D(TString xexpr, TString yexpr, const char *cut, TString drawopt, TH2 *hist) {
0412       // prep the machinery
0413       helper::ScannerBase scanner(objType);
0414       scanner.setIgnoreExceptions(ignoreExceptions_);
0415       if (!scanner.addExpression((const char *)xexpr))
0416         return nullptr;
0417       if (!scanner.addExpression((const char *)yexpr))
0418         return nullptr;
0419       if (strlen(cut))
0420         scanner.setCut(cut);
0421
0422       // check histo
0423       if (hist == nullptr) {
0424         std::cerr << "Method draw2D(xexpr, yexpr, cut, drawopt, hist) cannot be called with null 'hist'. Use the other "
0425                      "draw methods instead."
0426                   << std::endl;
0427         return nullptr;
0428       }
0429
0430       // fill histogram
0431       int iev = 0;
0432       for (event_->toBegin(), iev = 0; !event_->atEnd(); ++(*event_), ++iev) {
0433         if (maxEvents_ > -1 && iev > maxEvents_)
0434           break;
0435         if (!selectEvent(*event_))
0436           continue;
0437         handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0438         const Collection &vals = *handle_;
0439         for (size_t j = 0, n = vals.size(); j < n; ++j) {
0440           scanner.fill2D(&vals[j], hist);
0441         }
0442       }
0443
0444       if (!strlen(hist->GetTitle()))
0445         hist->SetTitle((strlen(cut) ? yexpr + ":" + xexpr + "{" + cut + "}" : yexpr + ":" + xexpr));
0446       if (!strlen(hist->GetXaxis()->GetTitle()))
0447         hist->GetXaxis()->SetTitle(xexpr);
0448       if (!strlen(hist->GetYaxis()->GetTitle()))
0449         hist->GetYaxis()->SetTitle(yexpr);
0450       if (!TString(drawopt).Contains("GOFF", TString::kIgnoreCase))
0451         hist->Draw(drawopt);
0452       return hist;
0453     }
0454     /// Just like draw() except that it uses TH2. Note that the order is (x,y) while in ROOT it's usually (y,x)!
0455     /// Note that automatical binning for TH2s is more expensive, as it requires to loop on the events twice!
0456     TH2 *draw2D(TString xexpr,
0457                 TString yexpr,
0458                 const char *cut = "",
0459                 TString drawopt = "",
0460                 const char *hname = "htemp",
0461                 TH2 *htemplate = nullptr) {
0462       TH2 *hist = nullptr;
0463       if (htemplate != nullptr) {
0464         if ((strcmp(hname, "htemp") == 0) && (strcmp(hname, htemplate->GetName()) != 0))
0465           htempDelete();
0466         hist = (TH2 *)htemplate->Clone(hname);
0467       } else if (drawopt.Contains("SAME", TString::kIgnoreCase)) {
0468         hist = getSameH2(hname);
0469       }
0470
0471       // if in the end we found no way to make "hist"
0472       if (hist == nullptr) {
0473         // prep the machinery
0474         helper::ScannerBase scanner(objType);
0475         scanner.setIgnoreExceptions(ignoreExceptions_);
0476         if (!scanner.addExpression((const char *)xexpr))
0477           return nullptr;
0478         if (!scanner.addExpression((const char *)yexpr))
0479           return nullptr;
0480         if (strlen(cut))
0481           scanner.setCut(cut);
0482
0483         if (strcmp(hname, "htemp") == 0)
0484           htempDelete();
0485         // ok this is much more a hack than for the 1D case
0486         double xmin = 0, xmax = -1, ymin = 0, ymax = -1;
0487         int iev;
0488         for (event_->toBegin(), iev = 0; !event_->atEnd(); ++(*event_), ++iev) {
0489           if (maxEvents_ > -1 && iev > maxEvents_)
0490             break;
0491           if (!selectEvent(*event_))
0492             continue;
0493           handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0494           const Collection &vals = *handle_;
0495           for (size_t j = 0, n = vals.size(); j < n; ++j) {
0496             if (!scanner.test(&vals[j]))
0497               continue;
0498             double x = scanner.eval(&vals[j], 0);
0499             double y = scanner.eval(&vals[j], 1);
0500             if ((xmax == -1) || (x >= xmax))
0501               xmax = x;
0502             if ((xmin == 0) || (x <= xmin))
0503               xmin = x;
0504             if ((ymax == -1) || (y >= ymax))
0505               ymax = y;
0506             if ((ymin == 0) || (y <= ymin))
0507               ymin = y;
0508           }
0509         }
0510         hist = new TH2F(hname,
0511                         "",
0512                         gEnv->GetValue("Hist.Binning.2D.x", 20),
0513                         xmin,
0514                         xmax,
0515                         gEnv->GetValue("Hist.Binning.2D.y", 20),
0516                         ymin,
0517                         ymax);
0518       }
0519       return draw2D(xexpr, yexpr, cut, drawopt, hist);
0520     }
0521
0522     /// Just like draw() except that it uses TH2. Note that the order is (x,y) while in ROOT it's usually (y,x)!
0523     TH2 *draw2D(TString xexpr,
0524                 int xbins,
0525                 double xlow,
0526                 double xhigh,
0527                 TString yexpr,
0528                 int ybins,
0529                 double ylow,
0530                 double yhigh,
0531                 const char *cut = "",
0532                 const char *drawopt = "",
0533                 const char *hname = "htemp") {
0534       if (TString(drawopt).Contains("SAME", TString::kIgnoreCase)) {
0535         std::cerr << "Binning is ignored when 'SAME' is specified." << std::endl;
0536         TH2 *hsame = getSameH2(hname);
0537         return draw2D(xexpr, yexpr, cut, drawopt, hsame);
0538       }
0539       if (strcmp(hname, "htemp") == 0)
0540         htempDelete();
0541       TH2 *htemp = new TH2F("htemp", "", xbins, xlow, xhigh, ybins, ylow, yhigh);
0542       return draw2D(xexpr, yexpr, cut, drawopt, htemp);
0543     }
0544
0545     /** Draw a scatter plot of x vs y for events passing the cut. */
0546     TGraph *drawGraph(TString xexpr, TString yexpr, const char *cut, TString drawopt, TGraph *graph) {
0547       // prep the machinery
0548       helper::ScannerBase scanner(objType);
0549       scanner.setIgnoreExceptions(ignoreExceptions_);
0550       if (!scanner.addExpression((const char *)xexpr))
0551         return nullptr;
0552       if (!scanner.addExpression((const char *)yexpr))
0553         return nullptr;
0554       if (strlen(cut))
0555         scanner.setCut(cut);
0556
0557       // make graph, if needed
0558       if (graph == nullptr) {
0559         graph = new TGraph();
0560         graph->SetNameTitle("htemp", (strlen(cut) ? yexpr + ":" + xexpr + "{" + cut + "}" : yexpr + ":" + xexpr));
0561       }
0562
0563       // fill graph
0564       int iev = 0;
0565       for (event_->toBegin(); !event_->atEnd(); ++(*event_), ++iev) {
0566         if (maxEvents_ > -1 && iev > maxEvents_)
0567           break;
0568         if (!selectEvent(*event_))
0569           continue;
0570         handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0571         const Collection &vals = *handle_;
0572         for (size_t j = 0, n = vals.size(); j < n; ++j) {
0573           scanner.fillGraph(&vals[j], graph);
0574         }
0575       }
0576
0577       if (!strlen(graph->GetTitle()))
0578         graph->SetTitle((strlen(cut) ? yexpr + ":" + xexpr + "{" + cut + "}" : yexpr + ":" + xexpr));
0579       if (!strlen(graph->GetXaxis()->GetTitle()))
0580         graph->GetXaxis()->SetTitle(xexpr);
0581       if (!strlen(graph->GetYaxis()->GetTitle()))
0582         graph->GetYaxis()->SetTitle(yexpr);
0583       if (!TString(drawopt).Contains("GOFF", TString::kIgnoreCase))
0584         graph->Draw(drawopt);
0585       return graph;
0586     }
0587
0588     /** Draw a scatter plot of x vs y for events passing the cut. */
0589     TGraph *drawGraph(
0590         TString xexpr, TString yexpr, const char *cut = "", TString drawopt = "AP", const char *gname = "htemp") {
0591       if (strcmp(gname, "htemp") == 0)
0592         htempDelete();
0593       TGraph *graph = new TGraph();
0594       graph->SetNameTitle(gname, (strlen(cut) ? yexpr + ":" + xexpr + "{" + cut + "}" : yexpr + ":" + xexpr));
0595       return drawGraph(xexpr, yexpr, cut, drawopt, graph);
0596     }
0597
0598     //------------------------------------------------------------------------------------------------------------------------------------
0599     /** Fill a RooDataSet.
0600              *  - Real variables are defined just like in the scan() command; a list separated by ":" (see also setExpressionSeparator());
0601              *  - Boolean variables are defined just like cuts, and are created as RooCategory with two states: pass(1) and fail(0).
0602              *  For each variable, the name is taken from the expression itself, or can be specified manuall by using the notation "@name=expr"
0603              *  Note: the dataset contains one entry per item, irrespectively of how entries are distributed among events.
0604              */
0605     RooDataSet *fillDataSet(const char *realvars,
0606                             const char *boolvars,
0607                             const char *cut = "",
0608                             const char *name = "data") {
0609       helper::ScannerBase scanner(objType);
0610       scanner.setIgnoreExceptions(ignoreExceptions_);
0611
0612       RooArgList vars;
0613       TObjArray *exprArray = TString(realvars).Tokenize(exprSep_);
0614       TObjArray *catArray = TString(boolvars).Tokenize(exprSep_);
0615       int nreals = exprArray->GetEntries();
0616       int nbools = catArray->GetEntries();
0617       for (int i = 0; i < nreals; ++i) {
0618         TString str = ((TObjString *)(*exprArray)[i])->GetString();
0619         std::string lb = str.Data();
0620         std::string ex = str.Data();
0621         if ((ex[0] == '@') && (ex.find('=') != std::string::npos)) {
0622           lb = lb.substr(1, ex.find('=') - 1);
0623           ex = ex.substr(ex.find('=') + 1);
0624         }
0625         if (!scanner.addExpression(ex.c_str())) {
0626           std::cerr << "Filed to define real variable '" << lb << "', expr = '" << ex << "'" << std::endl;
0627           return nullptr;
0628         }
0629         // FIXME: I have to leave it dangling on the HEAP otherwise ROOT segfaults...
0630         RooRealVar *var = new RooRealVar(lb.c_str(), lb.c_str(), 0.0);
0631         vars.add(*var);
0632       }
0633       for (int i = 0; i < nbools; ++i) {
0634         TString str = ((TObjString *)(*catArray)[i])->GetString();
0635         std::string lb = str.Data();
0636         std::string ex = str.Data();
0637         if ((ex[0] == '@') && (ex.find('=') != std::string::npos)) {
0638           lb = lb.substr(1, ex.find('=') - 1);
0639           ex = ex.substr(ex.find('=') + 1);
0640         }
0641         if (!scanner.addExtraCut(ex.c_str())) {
0642           std::cerr << "Filed to define bool variable '" << lb << "', cut = '" << ex << "'" << std::endl;
0643           return nullptr;
0644         }
0645         RooCategory *cat = new RooCategory(lb.c_str(), lb.c_str());
0646         cat->defineType("fail", 0);
0647         cat->defineType("pass", 1);
0648         vars.add(*cat);
0649       }
0650
0651       RooDataSet *ds = new RooDataSet(name, name, vars);
0652
0653       if (strlen(cut))
0654         scanner.setCut(cut);
0655       int iev = 0;
0656       for (event_->toBegin(); !event_->atEnd(); ++iev, ++(*event_)) {
0657         if (maxEvents_ > -1 && iev > maxEvents_)
0658           break;
0659         if (!selectEvent(*event_))
0660           continue;
0661         handle_.getByLabel(*event_, label_.c_str(), instance_.c_str(), process_.c_str());
0662         if (handle_.failedToGet()) {
0663           if (ignoreExceptions_)
0664             continue;
0665         }
0666         const Collection &vals = *handle_;
0667         for (size_t j = 0, n = vals.size(); j < n; ++j) {
0668           if (!scanner.test(&vals[j]))
0669             continue;
0670           for (int i = 0; i < nreals; ++i) {
0671             RooRealVar *var = (RooRealVar *)vars.at(i);
0672             var->setVal(scanner.eval(&vals[j], i));
0673           }
0674           for (int i = 0; i < nbools; ++i) {
0675             RooCategory *cat = (RooCategory *)vars.at(i + nreals);
0676             cat->setIndex(int(scanner.test(&vals[j], i + 1)));  // 0 is the event selection cut
0677           }
0678           ds->add(vars);
0679         }
0680       }
0681
0682       delete exprArray;
0683       delete catArray;
0684
0685       return ds;
0686     }
0687
0688     void setPrintFullEventId(bool printIt = true) { printFullEventId_ = printIt; }
0689     void setExpressionSeparator(TString separator) { exprSep_ = separator; }
0690     void setIgnoreExceptions(bool ignoreThem) { ignoreExceptions_ = ignoreThem; }
0691     void setMaxLinesToPrint(int lines) { maxLinesToPrint_ = (lines > 0 ? lines : 2147483647); }
0692
0693     void addEventSelector(fwlite::EventSelector *selector) { eventSelectors_.Add(selector); }
0694     void clearEventSelector() { eventSelectors_.Clear(); }
0695     TObjArray &eventSelectors() { return eventSelectors_; }
0696     bool selectEvent(const fwlite::EventBase &ev) const {
0697       for (int i = 0, n = eventSelectors_.GetEntries(); i < n; ++i) {
0698         if (!((fwlite::EventSelector *)(eventSelectors_[i]))->accept(ev))
0699           return false;
0700       }
0701       return true;
0702     }
0703
0704     void setMaxEvents(int max) { maxEvents_ = max; }
0705
0706   private:
0707     fwlite::EventBase *event_;
0708     std::string label_, instance_, process_;
0709     bool printFullEventId_;
0710     bool ignoreExceptions_;
0711     TString exprSep_;
0712     HandleT handle_;
0713     edm::TypeWithDict objType;
0714
0715     TObjArray eventSelectors_;
0716
0717     int maxEvents_;
0718
0719     int maxLinesToPrint_;
0720     bool wantMore() const {
0721       // ask if user wants more
0722       fprintf(stderr, "Type <CR> to continue or q to quit ==> ");
0723       // read first char
0724       int readch = getchar(), answer = readch;
0725       // poll out remaining chars from buffer
0726       while (readch != '\n' && readch != EOF)
0727         readch = getchar();
0728       // check first char
0729       return !(answer == 'q' || answer == 'Q');
0730     }
0731
0732     void htempDelete() {
0733       if (gDirectory) {
0734         TObject *obj = gDirectory->Get("htemp");
0735         if (obj)
0736           obj->Delete();
0737       }
0738     }
0739
0740     /// Get whatever histogram makes sense for a plot passing "SAME" in drawOpt, and call it hname
0741     /// Currently it won't work if the histogram of which we want to be "SAME" is not called "htemp"
0742     TH1 *getSameH1(const char *hname) {
0743       if (gDirectory && gDirectory->Get("htemp") != nullptr &&
0744           gDirectory->Get("htemp")->IsA()->InheritsFrom(TH1::Class())) {
0745         TH1 *hist = (TH1 *)((TH1 *)gDirectory->Get("htemp"))->Clone(hname);
0746         hist->Reset();
0747         hist->SetLineColor(kBlack);
0748         hist->SetMarkerColor(kBlack);
0749         return hist;
0750       } else {
0751         std::cerr << "There is no 'htemp' histogram from which to 'SAME'." << std::endl;
0752         return nullptr;
0753       }
0754     }
0755
0756     /// Get whatever histogram makes sense for a plot passing "SAME" in drawOpt, and call it hname
0757     /// Currently it won't work if the histogram of which we want to be "SAME" is not called "htemp"
0758     TH2 *getSameH2(const char *hname) {
0759       if (gDirectory && gDirectory->Get("htemp") != nullptr &&
0760           gDirectory->Get("htemp")->IsA()->InheritsFrom(TH2::Class())) {
0761         TH2 *hist = (TH2 *)((TH2 *)gDirectory->Get("htemp"))->Clone(hname);
0762         hist->Reset();
0763         hist->SetLineColor(kBlack);
0764         hist->SetMarkerColor(kBlack);
0765         return hist;
0766       } else {
0767         std::cerr << "There is no 'htemp' histogram from which to 'SAME'." << std::endl;
0768         return nullptr;
0769       }
0770     }
0771
0772     /// Get whatever histogram makes sense for a plot passing "SAME" in drawOpt, and call it hname
0773     /// Currently it won't work if the histogram of which we want to be "SAME" is not called "htemp"
0774     TProfile *getSameProf(const char *hname) {
0775       if (gDirectory && gDirectory->Get("htemp") != nullptr &&
0776           gDirectory->Get("htemp")->IsA()->InheritsFrom(TProfile::Class())) {
0777         TProfile *hist = (TProfile *)((TProfile *)gDirectory->Get("htemp"))->Clone(hname);
0778         hist->Reset();
0779         hist->SetLineColor(kBlack);
0780         hist->SetMarkerColor(kBlack);
0781         return hist;
0782       } else {
0783         std::cerr << "There is no 'htemp' histogram from which to 'SAME'." << std::endl;
0784         return nullptr;
0785       }
0786     }
0787   };
0788 }  // namespace fwlite
0789 #endif  // PhysicsTools_FWLite_Scanner_h