FWLite/examples/zPeak.cc

0001 // -*- C++ -*-
0002
0003 // CMS includes
0004 #include "FWCore/Utilities/interface/InputTag.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006 #include "DataFormats/PatCandidates/interface/Muon.h"
0007
0008 #include "PhysicsTools/FWLite/interface/EventContainer.h"
0009 #include "PhysicsTools/FWLite/interface/CommandLineParser.h"
0010
0011 // Root includes
0012 #include "TROOT.h"
0013
0014 using namespace std;
0015
0016 ///////////////////////////
0017 // ///////////////////// //
0018 // // Main Subroutine // //
0019 // ///////////////////// //
0020 ///////////////////////////
0021
0022 int main (int argc, char* argv[])
0023 {
0024    ////////////////////////////////
0025    // ////////////////////////// //
0026    // // Command Line Options // //
0027    // ////////////////////////// //
0028    ////////////////////////////////
0029
0030
0031    // Tell people what this analysis code does and setup default options.
0032    optutl::CommandLineParser parser ("Plot mass of Z candidates");
0033
0034    ////////////////////////////////////////////////
0035    // Change any defaults or add any new command //
0036    //      line options you would like here.     //
0037    ////////////////////////////////////////////////
0038    parser.stringValue ("outputFile") = "zpeak1.root";
0039
0040    // Parse the command line arguments
0041    parser.parseArguments (argc, argv);
0042
0043    //////////////////////////////////
0044    // //////////////////////////// //
0045    // // Create Event Container // //
0046    // //////////////////////////// //
0047    //////////////////////////////////
0048
0049    // This object 'event' is used both to get all information from the
0050    // event as well as to store histograms, etc.
0051    fwlite::EventContainer eventCont (parser);
0052
0053    ////////////////////////////////////////
0054    // ////////////////////////////////// //
0055    // //         Begin Run            // //
0056    // // (e.g., book histograms, etc) // //
0057    // ////////////////////////////////// //
0058    ////////////////////////////////////////
0059
0060    // Setup a style
0061    gROOT->SetStyle ("Plain");
0062
0063    // Book those histograms!
0064    eventCont.add( new TH1F ("Zmass", "Candidate Z mass", 50, 20, 220) );
0065
0066
0067    //////////////////////
0068    // //////////////// //
0069    // // Event Loop // //
0070    // //////////////// //
0071    //////////////////////
0072
0073    // create labels
0074    edm::InputTag muonLabel ("selectedLayer1Muons");
0075
0076    for (eventCont.toBegin(); ! eventCont.atEnd(); ++eventCont)
0077    {
0078       //////////////////////////////////
0079       // Take What We Need From Event //
0080       //////////////////////////////////
0081       edm::Handle<std::vector<pat::Muon> > muonHandle;
0082       eventCont.getByLabel (muonLabel, muonHandle);
0083       assert (muonHandle.isValid());
0084       vector< pat::Muon > const & muonVec = *muonHandle;
0085
0086       // make sure we have enough muons so that something bad doesn't
0087       // happen.
0088       if (muonVec.begin() == muonVec.end())
0089       {
0090          // If this statement is true, then we don't have at least one
0091          // muon, so don't bother.  The reason we need to do this is
0092          // that the logic below assumes that we have at least one and
0093          // we'll end up in an infinite loop without it.
0094          continue;
0095       }
0096
0097       // O.k.  Let's loop through our muons and see what we can find.
0098       const vector< pat::Muon >::const_iterator kEndIter       = muonVec.end();
0099       const vector< pat::Muon >::const_iterator kAlmostEndIter = kEndIter - 1;
0100       for (vector< pat::Muon >::const_iterator outerIter = muonVec.begin();
0101            kAlmostEndIter != outerIter;
0102            ++outerIter)
0103       {
0104          for (vector< pat::Muon >::const_iterator innerIter = outerIter + 1;
0105               kEndIter != innerIter;
0106               ++innerIter)
0107          {
0108             // make sure that we have muons of opposite charge
0109             if (outerIter->charge() * innerIter->charge() >= 0) continue;
0110
0111             // if we're here then we have one positively charged muon
0112             // and one negatively charged muon.
0113             eventCont.hist("Zmass")->Fill( (outerIter->p4() + innerIter->p4()).M() );
0114          } // for innerIter
0115       } // for outerIter
0116    } // for eventCont
0117
0118
0119    ////////////////////////
0120    // ////////////////// //
0121    // // Clean Up Job // //
0122    // ////////////////// //
0123    ////////////////////////
0124
0125    // Histograms will be automatically written to the root file
0126    // specificed by command line options.
0127
0128    // All done!  Bye bye.
0129    return 0;
0130 }