FWLite/examples/secvtxMassTemplates.cc

0001 // -*- C++ -*-
0002
0003 // CMS includes
0004 #include "DataFormats/FWLite/interface/Handle.h"
0005 #include "DataFormats/PatCandidates/interface/Jet.h"
0006 #include "DataFormats/PatCandidates/interface/Muon.h"
0007 #include "DataFormats/HepMCCandidate/interface/FlavorHistory.h"
0008 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0009 #include "Math/GenVector/PxPyPzM4D.h"
0010
0011 #include "PhysicsTools/FWLite/interface/EventContainer.h"
0012 #include "PhysicsTools/FWLite/interface/CommandLineParser.h"
0013
0014 // Root includes
0015 #include "TROOT.h"
0016 #include "TH2F.h"
0017
0018 using namespace std;
0019 using optutl::CommandLineParser;
0020
0021 enum
0022 {
0023    kNormalMode,
0024    kVqqMode,
0025    kLFMode,
0026    kWcMode
0027 };
0028
0029 //////////////////////////
0030 // Forward Declarations //
0031 //////////////////////////
0032
0033 // This subroutine, written by you (below), uses the command line
0034 // arguments and creates an output tag (if any).  This subroutine must
0035 // exist.
0036 void outputNameTagFunc (string &tag);
0037
0038 // Book all histograms to be filled this job.  If wanted, you can skip
0039 // this subroutine and book all histograms in the main subroutine.
0040 void bookHistograms (fwlite::EventContainer &eventCont);
0041
0042 // Calculate the name that should be used for this event based on the
0043 // mode, the HF word, and (if necessary), whether or not it's a W or
0044 // Z.  Returns false if the event should not be processed.
0045 bool calcSampleName (fwlite::EventContainer &eventCont, string &sampleName);
0046
0047 ///////////////////////////
0048 // ///////////////////// //
0049 // // Main Subroutine // //
0050 // ///////////////////// //
0051 ///////////////////////////
0052
0053 int main (int argc, char* argv[])
0054 {
0055    ////////////////////////////////
0056    // ////////////////////////// //
0057    // // Command Line Options // //
0058    // ////////////////////////// //
0059    ////////////////////////////////
0060
0061    // Tell people what this analysis code does and setup default options.
0062    CommandLineParser parser ("Creates SecVtx Mass templates");
0063
0064    //////////////////////////////////////////////////////
0065    // Add any command line options you would like here //
0066    //////////////////////////////////////////////////////
0067    parser.addOption ("mode",         CommandLineParser::kInteger,
0068                      "Normal(0), VQQ(1), LF(2), Wc(3)",
0069                      0);
0070    parser.addOption ("sampleName",   CommandLineParser::kString,
0071                      "Sample name (e.g., top, Wqq, etc.)");
0072
0073    ////////////////////////////////////////////////
0074    parser.stringValue ("outputFile") = "jetPt"; // .root added automatically
0075
0076    // Parse the command line arguments
0077    parser.parseArguments (argc, argv);
0078
0079    //////////////////////////////////
0080    // //////////////////////////// //
0081    // // Create Event Container // //
0082    // //////////////////////////// //
0083    //////////////////////////////////
0084
0085    // This object 'eventCont' is used both to get all information from the
0086    // eventCont as well as to store histograms, etc.
0087    // Parse the command line arguments
0088    parser.parseArguments (argc, argv);
0089
0090    //////////////////////////////////
0091    // //////////////////////////// //
0092    // // Create Event Container // //
0093    // //////////////////////////// //
0094    //////////////////////////////////
0095
0096    fwlite::EventContainer eventCont (parser);
0097
0098    ////////////////////////////////////////
0099    // ////////////////////////////////// //
0100    // //         Begin Run            // //
0101    // // (e.g., book histograms, etc) // //
0102    // ////////////////////////////////// //
0103    ////////////////////////////////////////
0104
0105    // Setup a style
0106    gROOT->SetStyle ("Plain");
0107
0108    // Book those histograms!
0109    bookHistograms (eventCont);
0110
0111    //////////////////////
0112    // //////////////// //
0113    // // Event Loop // //
0114    // //////////////// //
0115    //////////////////////
0116
0117    for (eventCont.toBegin(); !eventCont.atEnd(); ++eventCont)
0118    {
0119       //////////////////////////////////
0120       // Take What We Need From Event //
0121       //////////////////////////////////
0122       fwlite::Handle< vector< pat::Jet > > jetCollection;
0123       jetCollection.getByLabel (eventCont, "selectedLayer1Jets");
0124       assert ( jetCollection.isValid() );
0125
0126       fwlite::Handle< vector< pat::Muon > > goodMuonCollection;
0127       goodMuonCollection.getByLabel (eventCont, "goodLeptons");
0128       assert ( goodMuonCollection.isValid() );
0129
0130       // If we don't have any good leptons, don't bother
0131       if (! goodMuonCollection->size() )
0132       {
0133          continue;
0134       }
0135
0136       // get the sample name for this event
0137       string sampleName;
0138       if ( ! calcSampleName (eventCont, sampleName) )
0139       {
0140          // We don't want this one.
0141          continue;
0142       }
0143
0144       //////////////////////////////////////
0145       // Tagged Jets and Flavor Separator //
0146       //////////////////////////////////////
0147       int numBottom = 0, numCharm = 0, numLight = 0;
0148       int numTags = 0;
0149       double sumVertexMass = 0.;
0150       // Loop over the jets and find out which are tagged
0151       const vector< pat::Jet >::const_iterator kJetEnd = jetCollection->end();
0152       for (vector< pat::Jet >::const_iterator jetIter = jetCollection->begin();
0153            kJetEnd != jetIter;
0154            ++jetIter)
0155       {
0156          // Is this jet tagged and does it have a good secondary vertex
0157          if( jetIter->bDiscriminator("simpleSecondaryVertexBJetTags") < 2.05 )
0158          {
0159             // This jet isn't tagged
0160             continue;
0161          }
0162          reco::SecondaryVertexTagInfo const * svTagInfos
0163             = jetIter->tagInfoSecondaryVertex("secondaryVertex");
0164          if ( svTagInfos->nVertices() <= 0 )
0165          {
0166             // Given that we are using simple secondary vertex
0167             // tagging, I don't think this should ever happen.
0168             // Maybe we should put a counter here just to check.
0169             continue;
0170          } // if we have no secondary verticies
0171
0172          // count it
0173          ++numTags;
0174
0175          // What is the flavor of this jet
0176          int jetFlavor = std::abs( jetIter->partonFlavour() );
0177          if (5 == jetFlavor)
0178          {
0179             ++numBottom;
0180          } // if bottom
0181          else if (4 == jetFlavor)
0182          {
0183             ++numCharm;
0184          } // if charm
0185          else
0186          {
0187             ++numLight;
0188          } // if light flavor
0189
0190          ///////////////////////////
0191          // Calculate SecVtx Mass //
0192          ///////////////////////////
0193          ROOT::Math::LorentzVector< ROOT::Math::PxPyPzM4D<double> > sumVec;
0194          reco::CompositeCandidate vertexCand;
0195          reco::Vertex::trackRef_iterator
0196             kEndTracks = svTagInfos->secondaryVertex(0).tracks_end();
0197          for (reco::Vertex::trackRef_iterator trackIter =
0198                  svTagInfos->secondaryVertex(0).tracks_begin();
0199               trackIter != kEndTracks;
0200               ++trackIter )
0201          {
0202             const double kPionMass = 0.13957018;
0203             ROOT::Math::LorentzVector< ROOT::Math::PxPyPzM4D<double> > vec;
0204             vec.SetPx( (*trackIter)->px() );
0205             vec.SetPy( (*trackIter)->py() );
0206             vec.SetPz( (*trackIter)->pz() );
0207             vec.SetM (kPionMass);
0208             sumVec += vec;
0209          } // for trackIter
0210          sumVertexMass += sumVec.M();
0211          if (2 == numTags)
0212          {
0213             // We've got enough.  Stop.
0214             break;
0215          } // if we have enough tags
0216       } // for jet
0217
0218       ////////////////////////
0219       // General Accounting //
0220       ////////////////////////
0221       int numJets = std::min( (int) jetCollection->size(), 5 );
0222       eventCont.hist( sampleName + "_jettag")->Fill (numJets, numTags);
0223
0224       // If we don't have any tags, don't bother going on
0225       if ( ! numTags)
0226       {
0227          continue;
0228       }
0229
0230       ///////////////////////////////////////
0231       // Calculate average SecVtx mass and //
0232       // fill appropriate histograms.      //
0233       ///////////////////////////////////////
0234       sumVertexMass /= numTags;
0235       string whichtag = "";
0236       if (1 == numTags)
0237       {
0238          // single tag
0239          if      (numBottom)              whichtag = "_b";
0240          else if (numCharm)               whichtag = "_c";
0241          else if (numLight)               whichtag = "_q";
0242          else                             whichtag = "_X";
0243       } else {
0244          // double tags
0245          if      (2 == numBottom)         whichtag = "_bb";
0246          else if (2 == numCharm)          whichtag = "_cc";
0247          else if (2 == numLight)          whichtag = "_qq";
0248          else if (numBottom && numCharm)  whichtag = "_bc";
0249          else if (numBottom && numLight)  whichtag = "_bq";
0250          else if (numCharm  && numLight)  whichtag = "_bq";
0251          else                             whichtag = "_XX";
0252       } // if two tags
0253       string massName = sampleName
0254          + Form("_secvtxMass_%dj_%dt", numJets, numTags);
0255       eventCont.hist(massName)->Fill (sumVertexMass);
0256       eventCont.hist(massName + whichtag)->Fill (sumVertexMass);
0257    } // for eventCont
0258
0259    ////////////////////////
0260    // ////////////////// //
0261    // // Clean Up Job // //
0262    // ////////////////// //
0263    ////////////////////////
0264
0265    // Histograms will be automatically written to the root file
0266    // specificed by command line options.
0267
0268    // All done!  Bye bye.
0269    return 0;
0270 }
0271
0272
0273 //////////////  //////////////////////////////////  //////////////
0274 //////////////  // //////////////////////////// //  //////////////
0275 //////////////  // // Supporting Subroutines // //  //////////////
0276 //////////////  // //////////////////////////// //  //////////////
0277 //////////////  //////////////////////////////////  //////////////
0278
0279
0280 void outputNameTagFunc (string &tag)
0281 {
0282    // If you do not want to give you output filename any "tag" based
0283    // on the command line options, simply do nothing here.  This
0284    // function is designed to be called by fwlite::EventContainer constructor.
0285
0286    // if ( boolValue ("someCondition") )
0287    // {
0288    //    tag += "_someCond";
0289    // }
0290 }
0291
0292
0293 void bookHistograms (fwlite::EventContainer &eventCont)
0294 {
0295    /////////////////////////////////////////////
0296    // First, come up with all possible base   //
0297    // names (E.g., Wbb, Wb2, etc.).           //
0298    /////////////////////////////////////////////
0299    CommandLineParser &parser = eventCont.parser();
0300    CommandLineParser::SVec baseNameVec;
0301    CommandLineParser::SVec beginningVec, endingVec;
0302    switch ( parser.integerValue ("mode") )
0303    {
0304       case kVqqMode:
0305          // We want Wbb, Wb2, .., Zbb, ..  In this case, we completely
0306          // ignore the sampleName that was passed in.
0307          // Starts with
0308          beginningVec.push_back ("X");
0309          beginningVec.push_back ("W");
0310          beginningVec.push_back ("Z");
0311          // Ends with
0312          endingVec.push_back( "bb" );
0313          endingVec.push_back( "b2" );
0314          endingVec.push_back( "cc" );
0315          endingVec.push_back( "c2" );
0316          for (CommandLineParser::SVecConstIter outerIter = beginningVec.begin();
0317               beginningVec.end() != outerIter;
0318               ++outerIter)
0319          {
0320             for (CommandLineParser::SVecConstIter innerIter = endingVec.begin();
0321                  endingVec.end() != innerIter;
0322                  ++innerIter)
0323             {
0324                baseNameVec.push_back( *outerIter + *innerIter);
0325             } // for innerIter
0326          } // for outerIter
0327          break;
0328       case kLFMode:
0329          // just like the default case, except that we do have some
0330          // heavy flavor pieces here, too.
0331          baseNameVec.push_back(parser.stringValue ("sampleName") + "b3");
0332          baseNameVec.push_back(parser.stringValue ("sampleName") + "c3");
0333          // no break because to add just the name as well
0334       default:
0335          // We just want to use the sample name as it was given to us.
0336          baseNameVec.push_back(parser.stringValue ("sampleName"));
0337    } // for switch
0338
0339    ////////////////////////////////////////
0340    // Now the different tagging endings. //
0341    ////////////////////////////////////////
0342    CommandLineParser::SVec singleTagEndingVec, doubleTagEndingVec;
0343    singleTagEndingVec.push_back ("_b");
0344    singleTagEndingVec.push_back ("_c");
0345    singleTagEndingVec.push_back ("_q");
0346    doubleTagEndingVec.push_back ("_bb");
0347    doubleTagEndingVec.push_back ("_cc");
0348    doubleTagEndingVec.push_back ("_qq");
0349    doubleTagEndingVec.push_back ("_bc");
0350    doubleTagEndingVec.push_back ("_bq");
0351    doubleTagEndingVec.push_back ("_cq");
0352
0353    /////////////////////////////////////////
0354    // Finally, let's put it all together. //
0355    /////////////////////////////////////////
0356    for (CommandLineParser::SVecConstIter baseIter = baseNameVec.begin();
0357         baseNameVec.end() != baseIter;
0358         ++baseIter)
0359    {
0360       //////////////////////////////////////////////////////
0361       // For each flavor, one jet/tag counting histogram. //
0362       //////////////////////////////////////////////////////
0363       TString histName = *baseIter + "_jettag";
0364       eventCont.add( new TH2F( histName, histName,
0365                                6, -0.5, 5.5,
0366                                3, -0.5, 2.5) );
0367       for (int jet = 1; jet <= 5; ++jet)
0368       {
0369          for (int tag = 1; tag <= 2; ++tag)
0370          {
0371             ////////////////////////////////////////////
0372             // For each jet/tag, a single secvtx mass //
0373             ////////////////////////////////////////////
0374             if (tag > jet) continue;
0375             histName = *baseIter + Form ("_secvtxMass_%dj_%dt", jet, tag);
0376             eventCont.add( new TH1F( histName, histName, 40, 0, 10) );
0377             CommandLineParser::SVec *vecPtr = &singleTagEndingVec;
0378             if (2 == tag)
0379             {
0380                vecPtr = &doubleTagEndingVec;
0381             } // double tag
0382             for (CommandLineParser::SVecConstIter tagIter = vecPtr->begin();
0383                  vecPtr->end() != tagIter;
0384                  ++tagIter)
0385             {
0386                ////////////////////////////////////////////////////
0387                // And different secvtx mass for each tag ending. //
0388                ////////////////////////////////////////////////////
0389                histName = *baseIter + Form ("_secvtxMass_%dj_%dt", jet, tag)
0390                   + *tagIter;
0391                eventCont.add( new TH1F( histName, histName, 40, 0, 10) );
0392             } // for tagIter
0393          } // for tag
0394       } // for jet
0395    } // for baseIter
0396 }
0397
0398
0399 bool calcSampleName (fwlite::EventContainer &eventCont, string &sampleName)
0400 {
0401    // calculate sample name
0402    CommandLineParser &parser = eventCont.parser();
0403    sampleName = parser.stringValue  ("sampleName");
0404    int mode   = parser.integerValue ("mode");
0405
0406    /////////////////
0407    // Normal Mode //
0408    //// /////////////
0409    if (kNormalMode == mode)
0410    {
0411       // all we want is the sample name, so in this case we're done.
0412       return true;
0413    }
0414    // Get the heavy flavor category
0415    fwlite::Handle< unsigned int > heavyFlavorCategory;
0416    heavyFlavorCategory.getByLabel (eventCont, "flavorHistoryFilter");
0417    assert ( heavyFlavorCategory.isValid() );
0418    int HFcat = (*heavyFlavorCategory);
0419
0420    ///////////////////////
0421    // Light Flavor Mode //
0422    ///////////////////////
0423    if (kLFMode == mode)
0424    {
0425       // Wqq
0426       if (5 == HFcat)
0427       {
0428          sampleName += "b3";
0429       } else if (6 == HFcat)
0430       {
0431          sampleName += "c3";
0432       } else if (11 != HFcat)
0433       {
0434          // skip this event
0435          return false;
0436       } // else if ! 11
0437       return true;
0438    }
0439
0440    /////////////
0441    // Wc Mode //
0442    /////////////
0443    if (kWcMode == mode)
0444    {
0445       // Wc
0446       if (4 != HFcat)
0447       {
0448          // skip this event
0449          return false;
0450       } // if not Wc
0451       return true;
0452    } // else if Wc
0453
0454    //////////////
0455    // Vqq Mode //
0456    //////////////
0457    // MadGraph (at least as CMS has implemented it) has this _lovely_
0458    // feature that if the W or Z is far enough off-shell, it erases
0459    // the W or Z from the event record.  This means that in some
0460    // number of cases, we won't be able to tell whether this is a W or
0461    // Z event by looking for a W or Z in the GenParticle collection.
0462    // (We'll eventually have to be more clever).
0463    sampleName = "X";
0464    fwlite::Handle< vector< reco::GenParticle > > genParticleCollection;
0465    genParticleCollection.getByLabel (eventCont, "genParticles");
0466    assert ( genParticleCollection.isValid() );
0467    // We don't know if it is a W, a Z, or neither
0468    // Iterate over genParticles
0469    const vector< reco::GenParticle>::const_iterator
0470       kGenPartEnd = genParticleCollection->end();
0471    for (vector< reco::GenParticle>::const_iterator gpIter =
0472            genParticleCollection->begin();
0473         gpIter != kGenPartEnd; ++gpIter )
0474    {
0475       if (gpIter->status() == 3 && std::abs(gpIter->pdgId()) == 23)
0476       {
0477          sampleName = "Z";
0478          break;
0479       }
0480       if (gpIter->status() == 3 && std::abs(gpIter->pdgId()) == 24)
0481       {
0482          sampleName = "W";
0483          break;
0484       }
0485    } // for  gpIter
0486    switch (HFcat)
0487    {
0488       // from:
0489       // https://twiki.cern.ch/twiki/bin/viewauth/CMS/SWGuideFlavorHistory
0490       //  1. W+bb with >= 2 jets from the ME (dr > 0.5)
0491       //  2. W+b or W+bb with 1 jet from the ME
0492       //  3. W+cc from the ME (dr > 0.5)
0493       //  4. W+c or W+cc with 1 jet from the ME
0494       //  5. W+bb with 1 jet from the parton shower (dr == 0.0)
0495       //  6. W+cc with 1 jet from the parton shower (dr == 0.0)
0496       //  7. W+bb with >= 2 partons but 1 jet from the ME (dr == 0.0)
0497       //  8. W+cc with >= 2 partons but 1 jet from the ME (dr == 0.0)
0498       //  9. W+bb with >= 2 partons but 2 jets from the PS (dr > 0.5)
0499       // 10. W+cc with >= 2 partons but 2 jets from the PS (dr > 0.5)
0500       // 11. Veto of all the previous (W+ light jets)
0501       case 1:
0502          sampleName += "bb";
0503          break;
0504       case 2:
0505          // Sometimes this is referred to as 'b' (e.g., 'Wb'), but I
0506          // am using the suffix '2' to keep this case clear for when
0507          // we have charm (see below).
0508          sampleName += "b2";
0509          break;
0510       case 3:
0511          sampleName += "cc";
0512          break;
0513       case 4:
0514          // We want to keep this case clear from real W + single charm
0515          // produced (as opposed to two charm quarks produced and one
0516          // goes down the beampipe), so we use 'c2' instead of 'c'.
0517          sampleName += "c2";
0518          break;
0519       default:
0520          // we don't want the rest of the cases.  Return an empty
0521          // string so we know.
0522          return false;
0523    } // switch HFcat
0524    return true;
0525 }