HepMCCandAlgos/test/HZZ4muAnalysisExample.C

0001
0002 {
0003
0004 // Attention :
0005 // This example will ONLY work if you execute GenH190... cfg
0006 // with the option of writing edm::Event out (via PoolOutputModule)
0007
0008
0009 // No need for this here anymore,
0010 // since it already in the rootlogon.C
0011 //
0012 //   gSystem->Load("libFWCoreFWLite") ;
0013 //   FWLiteEnabler::enable() ;
0014
0015    TFile* f = new TFile("../../../Configuration/Examples/data/pythiaH190ZZ4mu.root") ;
0016
0017    TTree* tevt = f->Get("Events") ;
0018
0019    int nevt = tevt->GetEntries() ;
0020    std::cout << " Number of Events = " << nevt << std::endl ;
0021
0022    cout << " We will print out some information for the 1st event" << endl ;
0023
0024
0025    edm::HepMCProduct EvtProd ;
0026    TBranch* bhepmc =
0027       tevt->GetBranch( "edmHepMCProduct_source__Gen.obj") ;
0028
0029    bhepmc->SetAddress( & EvtProd ) ;
0030
0031    HepMC::GenEvent*    Evt = 0 ;
0032    HepMC::GenVertex*   Vtx = 0 ;
0033    HepMC::GenParticle* Part = 0 ;
0034    int                 NVtx = 0 ;
0035    int                 NPrt = 0 ;
0036
0037    HepMC::GenVertex*   HiggsDecVtx = 0 ;
0038
0039    HepMC::GenEvent::particle_iterator pit ;
0040    HepMC::GenEvent::vertex_iterator vit ;
0041    HepMC::GenVertex::particle_iterator vpit ;
0042
0043    // Somehow I'm having problems when trying to store (pointers to) selected muons
0044    // in a containers std::vector<HepMC::GenParticle*> (despite existing dictionary...),
0045    // thus I'm doing a dirty trick of storing momenta in one container (std::vector)
0046    // and the PDG's (for checking charge) in another (TArrayI)
0047    //
0048    // I'm still using vector of CLHEP 4-vectors because there's NO dictionary for
0049    // std::vector<HepMC::FourVector>
0050    // I'll fix it later
0051    //
0052    std::vector<CLHEP::HepLorentzVector> StableMuMom ;
0053
0054    //
0055    // I'm reserving it for 20 int's... I wish it was automatic
0056    //
0057    TArrayI StablePDG(20) ;
0058
0059    TH1D* Hist2muMass = new TH1D( "Hist2muMass", "test 2-mu mass", 100,  60., 120. ) ;
0060    TH1D* Hist4muMass = new TH1D( "Hist4muMass", "test 4-mu mass", 100, 170., 210. ) ;
0061
0062
0063    for ( int iev=0; iev<nevt; iev++ )
0064    {
0065
0066       bhepmc->GetEntry(iev) ;
0067
0068       Evt = EvtProd->GetEvent() ;
0069
0070       if ( iev == 0 ) Evt->print() ;
0071
0072       HiggsDecVtx = 0 ;
0073       StableMuMom.clear() ;
0074
0075       for ( pit=Evt->particles_begin(); pit!=Evt->particles_end(); pit++)
0076       {
0077      Part = (*pit) ;
0078      if ( Part->pdg_id() == 25 )
0079      {
0080         HiggsDecVtx = Part->end_vertex() ;
0081         if ( HiggsDecVtx != 0 )
0082         {
0083            if ( iev == 0 ) cout << " Found Higgs with valid decay vertex : "
0084                                 << HiggsDecVtx->barcode()
0085                             << " in event " << iev << endl ;
0086            break ;
0087         }
0088      }
0089       }
0090
0091       if ( HiggsDecVtx == 0 )
0092       {
0093          cout << "There is NO Higgs in the event " << iev << endl ;
0094       }
0095       else
0096       {
0097          for ( pit=Evt->particles_begin(); pit!=Evt->particles_end(); pit++)
0098      {
0099         Part = (*pit) ;
0100         if ( abs(Part->pdg_id()) == 13 && Part->status() == 1 )
0101         {
0102            // here's the "dirty trick" itself, storing info
0103            // in 2 differnt containers...
0104            //
0105            StableMuMom.push_back( CLHEP::HepLorentzVector(Part->momentum().px(),
0106                                                           Part->momentum().py(),
0107                                   Part->momentum().pz(),
0108                                   Part->momentum().e() ) ) ;
0109            StablePDG[StableMuMom.size()-1] = Part->pdg_id() ;
0110         }
0111      }
0112       }
0113
0114       if ( iev == 0 ) cout << " Found " << StableMuMom.size() << " stable muons" << endl ;
0115
0116       HepMC::FourVector Mom2mu ;
0117
0118       for ( unsigned int i=0; i<StableMuMom.size(); i++ )
0119       {
0120      for ( unsigned int j=i+1; j<StableMuMom.size(); j++ )
0121      {
0122         if ( StablePDG.At(i)*StablePDG.At(j) > 0 ) continue ; // skip same charge pairs
0123         if ( iev == 0 )
0124         {
0125            cout << " Stable particles id-s: " << StablePDG.At(i) << " "
0126                                               << StablePDG.At(j) << endl ;
0127             }
0128
0129         double XMass2mu = HepMC::FourVector( (StableMuMom[i].px()+StableMuMom[j].px()),
0130                                     (StableMuMom[i].py()+StableMuMom[j].py()),
0131                                 (StableMuMom[i].pz()+StableMuMom[j].pz()),
0132                             (StableMuMom[i].e() +StableMuMom[j].e()) ).m() ;
0133
0134         if ( iev == 0 ) cout << " 2-mu inv.mass = " << XMass2mu << endl ;
0135         Hist2muMass->Fill( XMass2mu ) ;
0136      }
0137       }
0138
0139       if ( StableMuMom.size() == 4 )
0140       {
0141          double px, py, pz, e ;
0142      px=py=pz=e=0. ;
0143      for ( unsigned int i=0; i<4; i++ )
0144      {
0145         //Mom4mu += StableMuMom[i] ;
0146         px += StableMuMom[i].px() ;
0147         py += StableMuMom[i].py() ;
0148         pz += StableMuMom[i].pz() ;
0149         e  += StableMuMom[i].e() ;
0150      }
0151      double XMass4mu = HepMC::FourVector( px, py, pz, e ).m() ;
0152      Hist4muMass->Fill( XMass4mu ) ;
0153          if ( iev == 0 ) cout << " 4-mu inv.mass = " << XMass4mu << endl ;
0154       }
0155    }
0156
0157    TCanvas* cnv = new TCanvas("cnv") ;
0158
0159    cnv->Divide(2,1) ;
0160
0161    cnv->cd(1) ;
0162    Hist2muMass->Draw() ;
0163    cnv->cd(2) ;
0164    Hist4muMass->Draw() ;
0165
0166
0167 }