Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-06 00:35:31

0001 // -*- C++ -*-
0002 //
0003 // Package:    GeneratorInterface
0004 // Class:      EvtGenTestAnalyzer
0005 //
0006 //
0007 // Description: Module to analyze Pythia-EvtGen HepMCproducts
0008 //
0009 //
0010 // Original Author:  Roberto Covarelli
0011 //         Created:  April 26, 2007
0012 //
0013 
0014 #include <iostream>
0015 #include <fstream>
0016 
0017 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0018 
0019 // essentials !!!
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "DataFormats/Common/interface/Handle.h"
0022 
0023 #include "TFile.h"
0024 #include "TH1.h"
0025 #include "TF1.h"
0026 #include "TLorentzVector.h"
0027 #include "TVector3.h"
0028 #include "TObjArray.h"
0029 
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031 #include "FWCore/ServiceRegistry/interface/Service.h"
0032 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0033 
0034 #include "GeneratorInterface/ExternalDecays/test/EvtGenTestAnalyzer.h"
0035 
0036 EvtGenTestAnalyzer::EvtGenTestAnalyzer(const edm::ParameterSet& pset)
0037     : fOutputFileName(pset.getUntrackedParameter<std::string>("HistOutFile", std::string("TestBs.root"))),
0038       tokenHepMC_(consumes<edm::HepMCProduct>(
0039           edm::InputTag(pset.getUntrackedParameter("moduleLabel", std::string("generator")), "unsmeared"))),
0040       fOutputFile(0) {
0041   usesResource(TFileService::kSharedResource);
0042 }
0043 
0044 void EvtGenTestAnalyzer::beginJob() {
0045   nevent = 0;
0046   nbs = 0;
0047 
0048   edm::Service<TFileService> fs;
0049   // fOutputFile = new TFile(fOutputFileName.c_str(), "RECREATE");
0050   // fHist2muMass  = new TH1D(  "Hist2muMass", "2-mu inv. mass", 100,  60., 120. ) ;
0051   hGeneralId = fs->make<TH1D>("hGeneralId", "LundIDs of all particles", 100, -1000., 1000.);
0052   hnB = fs->make<TH1D>("hnB", "N(B)", 10, 0., 10.);
0053   hnJpsi = fs->make<TH1D>("hnJpsi", "N(Jpsi)", 10, 0., 10.);
0054   hnBz = fs->make<TH1D>("hnBz", "N(B0)", 10, 0., 10.);
0055   hnBzb = fs->make<TH1D>("hnBzb", "N(B0bar)", 10, 0., 10.);
0056   hMinvb = fs->make<TH1D>("hMinvb", "B invariant mass", 100, 5.0, 6.0);
0057   hPtbs = fs->make<TH1D>("hPtbs", "Pt Bs", 100, 0., 50.);
0058   hPbs = fs->make<TH1D>("hPbs", "P Bs", 100, 0., 200.);
0059   hPhibs = fs->make<TH1D>("hPhibs", "Phi Bs", 100, -3.14, 3.14);
0060   hEtabs = fs->make<TH1D>("hEtabs", "Eta Bs", 100, -7.0, 7.0);
0061   hPtmu = fs->make<TH1D>("hPtmu", "Pt Mu", 100, 0., 50.);
0062   hPmu = fs->make<TH1D>("hPmu", "P Mu", 100, 0., 200.);
0063   hPhimu = fs->make<TH1D>("hPhimu", "Phi Mu", 100, -3.14, 3.14);
0064   hEtamu = fs->make<TH1D>("hEtamu", "Eta Mu", 100, -7.0, 7.0);
0065   hPtRadPho = fs->make<TH1D>("hPtRadPho", "Pt radiated photon", 100, 0., 200.);
0066   hPhiRadPho = fs->make<TH1D>("hPhiRadPho", "Phi radiated photon", 100, -3.14, 3.14);
0067   hEtaRadPho = fs->make<TH1D>("hEtaRadPho", "Eta radiated photon", 100, -7.0, 7.0);
0068   htbPlus = fs->make<TH1D>("htbPlus", "B+ proper decay time", 50, 0., 12.);
0069   htbUnmix = fs->make<TH1D>("htbUnmix", "B0 proper decay time (unmixed)", 50, 0., 12.);
0070   htbMix = fs->make<TH1D>("htbMix", "B0 proper decay time (mixed)", 50, 0., 12.);
0071   htbMixPlus = fs->make<TH1D>("htbMixPlus", "B0 proper decay time (mixed positive)", 50, 0., 12.);
0072   htbMixMinus = fs->make<TH1D>("htbMixMinus", "B0 proper decay time (mixed negative)", 50, 0., 12.);
0073   htbsUnmix = fs->make<TH1D>("htbsUnmix", "Bs proper decay time (unmixed)", 50, 0., 12.);
0074   htbsMix = fs->make<TH1D>("htbsMix", "Bs proper decay time (mixed)", 50, 0., 12.);
0075   htbJpsiKs = fs->make<TH1D>("htbJpsiKs", "B0 -> J/#psiK_{s} decay time (B0 tags)", 50, 0., 12.);
0076   htbbarJpsiKs = fs->make<TH1D>("htbbarJpsiKs", "B0 -> J/#psiK_{s} decay time (B0bar tags)", 50, 0., 12.);
0077   hmumuMassSqr = fs->make<TH1D>("hmumuMassSqr", "#mu^{+}#mu^{-} invariant mass squared", 100, -1.0, 25.0);
0078   hmumuMassSqrPlus =
0079       fs->make<TH1D>("hmumuMassSqrPlus", "#mu^{+}#mu^{-} invariant mass squared (cos#theta > 0)", 100, -1.0, 25.0);
0080   hmumuMassSqrMinus =
0081       fs->make<TH1D>("hmumuMassSqrMinus", "#mu^{+}#mu^{-} invariant mass squared (cos#theta < 0)", 100, -1.0, 25.0);
0082   hIdBsDaugs = fs->make<TH1D>("hIdBsDaugs", "LundIDs of the Bs's daughters", 100, -1000., 1000.);
0083   hIdPhiDaugs = fs->make<TH1D>("hIdPhiDaugs", "LundIDs of the phi's daughters", 100, -500., 500.);
0084   hIdJpsiMot = fs->make<TH1D>("hIdJpsiMot", "LundIDs of the J/psi's mother", 100, -500., 500.);
0085   hIdBDaugs = fs->make<TH1D>("hIdBDaugs", "LundIDs of the B's daughters", 100, -1000., 1000.);
0086   hCosTheta1 = fs->make<TH1D>("hCosTheta1", "cos#theta_{1}", 50, -1., 1.);
0087   hCosTheta2 = fs->make<TH1D>("hCosTheta2", "cos#theta_{2}", 50, -1., 1.);
0088   hPhi1 = fs->make<TH1D>("hPhi1", "#phi_{1}", 50, -3.14, 3.14);
0089   hPhi2 = fs->make<TH1D>("hPhi2", "#phi_{2}", 50, -3.14, 3.14);
0090   hCosThetaLambda = fs->make<TH1D>("hCosThetaLambda", "cos#theta_{#Lambda}", 50, -1., 1.);
0091 
0092   decayed = new std::ofstream("decayed.txt");
0093   undecayed = new std::ofstream("undecayed.txt");
0094   return;
0095 }
0096 
0097 void EvtGenTestAnalyzer::analyze(const edm::Event& e, const edm::EventSetup&) {
0098   const edm::Handle<edm::HepMCProduct>& EvtHandle = e.getHandle(tokenHepMC_);
0099 
0100   const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0101   if (Evt)
0102     nevent++;
0103 
0104   const float mmcToPs = 3.3355;
0105   int nB = 0;
0106   int nJpsi = 0;
0107   int nBp = 0;
0108   int nBz = 0;
0109   int nBzb = 0;
0110   int nBzmix = 0;
0111   int nBzunmix = 0;
0112   int nBzKmumu = 0;
0113   int nBJpsiKs = 0;
0114   int nBJpsiKstar = 0;
0115 
0116   for (HepMC::GenEvent::particle_const_iterator p = Evt->particles_begin(); p != Evt->particles_end(); ++p) {
0117     // General
0118     TLorentzVector thePart4m(0., 0., 0., 0.);
0119     HepMC::GenVertex* endvert = (*p)->end_vertex();
0120     HepMC::GenVertex* prodvert = (*p)->production_vertex();
0121     float gamma = (*p)->momentum().e() / (*p)->generated_mass();
0122     float dectime = 0.0;
0123     int mixed = -1;  // mixed is: -1 = unmixed
0124                      //           0 = mixed (before mixing)
0125                      //           1 = mixed (after mixing)
0126     if (endvert && prodvert) {
0127       dectime = (endvert->position().t() - prodvert->position().t()) * mmcToPs / gamma;
0128 
0129       // Mixed particle ?
0130       for (HepMC::GenVertex::particles_in_const_iterator p2 = prodvert->particles_in_const_begin();
0131            p2 != prodvert->particles_in_const_end();
0132            ++p2) {
0133         if ((*p)->pdg_id() + (*p2)->pdg_id() == 0) {
0134           mixed = 1;
0135           gamma = (*p2)->momentum().e() / (*p2)->generated_mass();
0136           HepMC::GenVertex* mixvert = (*p2)->production_vertex();
0137           dectime = (prodvert->position().t() - mixvert->position().t()) * mmcToPs / gamma;
0138         }
0139       }
0140       for (HepMC::GenVertex::particles_out_const_iterator ap = endvert->particles_out_const_begin();
0141            ap != endvert->particles_out_const_end();
0142            ++ap) {
0143         if ((*p)->pdg_id() + (*ap)->pdg_id() == 0)
0144           mixed = 0;
0145       }
0146     }
0147 
0148     hGeneralId->Fill((*p)->pdg_id());
0149 
0150     // --------------------------------------------------------------
0151     if (std::abs((*p)->pdg_id()) == 521)  // B+/-
0152     // || abs((*p)->pdg_id()/100) == 4 || abs((*p)->pdg_id()/100) == 3)
0153     {
0154       nBp++;
0155       htbPlus->Fill(dectime);
0156       int isJpsiKstar = 0;
0157       for (HepMC::GenVertex::particles_out_const_iterator bp = endvert->particles_out_const_begin();
0158            bp != endvert->particles_out_const_end();
0159            ++bp) {
0160         if ((*bp)->pdg_id() == 443 || std::abs((*bp)->pdg_id()) == 323)
0161           isJpsiKstar++;
0162       }
0163       if (isJpsiKstar == 2)
0164         nBJpsiKstar++;
0165     }
0166     // if ((*p)->pdg_id() == 443) *undecayed << (*p)->pdg_id() << std::endl;
0167     // --------------------------------------------------------------
0168 
0169     if ((*p)->pdg_id() == 531 /* && endvert */) {  // B_s
0170       // nbs++;
0171       hPtbs->Fill((*p)->momentum().perp());
0172       hPbs->Fill(sqrt(pow((*p)->momentum().px(), 2) + pow((*p)->momentum().py(), 2) + pow((*p)->momentum().pz(), 2)));
0173       hPhibs->Fill((*p)->momentum().phi());
0174       hEtabs->Fill((*p)->momentum().pseudoRapidity());
0175 
0176       for (HepMC::GenVertex::particles_out_const_iterator ap = endvert->particles_out_const_begin();
0177            ap != endvert->particles_out_const_end();
0178            ++ap) {
0179         hIdBsDaugs->Fill((*ap)->pdg_id());
0180       }
0181 
0182       if (mixed == 1) {
0183         htbsMix->Fill(dectime);
0184       } else if (mixed == -1) {
0185         htbsUnmix->Fill(dectime);
0186       }
0187     }
0188     // --------------------------------------------------------------
0189     if (std::abs((*p)->pdg_id()) == 511 /* && endvert */) {  // B0
0190       if (mixed != 0) {
0191         nB++;
0192         if ((*p)->pdg_id() > 0) {
0193           nBz++;
0194           if (mixed == 1) {
0195             nBzmix++;
0196           } else {
0197             nBzunmix++;
0198           }
0199         } else {
0200           nBzb++;
0201         }
0202       }
0203       int isJpsiKs = 0;
0204       int isKmumu = 0;
0205       int isSemilept = 0;
0206       for (HepMC::GenVertex::particles_out_const_iterator bp = endvert->particles_out_const_begin();
0207            bp != endvert->particles_out_const_end();
0208            ++bp) {
0209         // Check invariant mass consistency ...
0210         TLorentzVector theDaug4m(
0211             (*bp)->momentum().px(), (*bp)->momentum().py(), (*bp)->momentum().pz(), (*bp)->momentum().e());
0212         thePart4m += theDaug4m;
0213         hIdBDaugs->Fill((*bp)->pdg_id());
0214         if ((*bp)->pdg_id() == 443 || (*bp)->pdg_id() == 310)
0215           isJpsiKs++;
0216         if ((*bp)->pdg_id() == 22) {
0217           hPtRadPho->Fill((*bp)->momentum().perp());
0218           hPhiRadPho->Fill((*bp)->momentum().phi());
0219           hEtaRadPho->Fill((*bp)->momentum().pseudoRapidity());
0220         }
0221         if ((*p)->pdg_id() > 0 && (std::abs((*bp)->pdg_id()) == 313 || std::abs((*bp)->pdg_id()) == 13))
0222           isKmumu++;
0223         if (std::abs((*bp)->pdg_id()) == 11 || std::abs((*bp)->pdg_id()) == 13 || std::abs((*bp)->pdg_id()) == 15)
0224           isSemilept++;
0225       }
0226 
0227       hMinvb->Fill(sqrt(thePart4m.M2()));
0228       /* if (fabs(sqrt(thePart4m.M2())-5.28) > 0.05) {
0229      *undecayed << sqrt(thePart4m.M2()) << "  " << (*p)->pdg_id() << " --> " ;
0230          for ( GenVertex::particles_out_const_iterator bp = endvert->particles_out_const_begin(); bp != endvert->particles_out_const_end(); ++bp ) {
0231        *undecayed << (*bp)->pdg_id() << " " ;
0232          }
0233          *undecayed << std::endl;
0234        } */
0235 
0236       if (isSemilept) {
0237         if (mixed == 1) {
0238           htbMix->Fill(dectime);
0239           if ((*p)->pdg_id() > 0) {
0240             htbMixPlus->Fill(dectime);
0241           } else {
0242             htbMixMinus->Fill(dectime);
0243           }
0244         } else if (mixed == -1) {
0245           htbUnmix->Fill(dectime);
0246         }
0247       }
0248 
0249       if (isJpsiKs == 2) {
0250         nBJpsiKs++;
0251         if ((*p)->pdg_id() * mixed < 0) {
0252           htbbarJpsiKs->Fill(dectime);
0253         } else {
0254           htbJpsiKs->Fill(dectime);
0255         }
0256       }
0257 
0258       if (isKmumu == 3)
0259         nBzKmumu++;
0260     }
0261     // --------------------------------------------------------------
0262     if ((*p)->pdg_id() == 443) {  //Jpsi
0263       nJpsi++;
0264       for (HepMC::GenVertex::particles_in_const_iterator ap = prodvert->particles_in_const_begin();
0265            ap != prodvert->particles_in_const_end();
0266            ++ap) {
0267         hIdJpsiMot->Fill((*ap)->pdg_id());
0268       }
0269     }
0270     // --------------------------------------------------------------
0271     if ((*p)->pdg_id() == 333) {  // phi
0272       if (endvert) {
0273         for (HepMC::GenVertex::particles_out_const_iterator cp = endvert->particles_out_const_begin();
0274              cp != endvert->particles_out_const_end();
0275              ++cp) {
0276           hIdPhiDaugs->Fill((*cp)->pdg_id());
0277         }
0278       }
0279     }
0280     // --------------------------------------------------------------
0281     if ((*p)->pdg_id() == 13) {  // mu+
0282       for (HepMC::GenVertex::particles_in_const_iterator p2 = prodvert->particles_in_const_begin();
0283            p2 != prodvert->particles_in_const_end();
0284            ++p2) {
0285         if (std::abs((*p2)->pdg_id()) == 511) {  // B0
0286           hPtmu->Fill((*p)->momentum().perp());
0287           hPmu->Fill(
0288               sqrt(pow((*p)->momentum().px(), 2) + pow((*p)->momentum().py(), 2) + pow((*p)->momentum().pz(), 2)));
0289           hPhimu->Fill((*p)->momentum().phi());
0290           hEtamu->Fill((*p)->momentum().pseudoRapidity());
0291           for (HepMC::GenVertex::particles_out_const_iterator p3 = prodvert->particles_out_const_begin();
0292                p3 != prodvert->particles_out_const_end();
0293                ++p3) {
0294             if ((*p3)->pdg_id() == -13) {  // mu-
0295               TLorentzVector pmu1(
0296                   (*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0297               TLorentzVector pmu2(
0298                   (*p3)->momentum().px(), (*p3)->momentum().py(), (*p3)->momentum().pz(), (*p3)->momentum().e());
0299               TLorentzVector pb0(
0300                   (*p2)->momentum().px(), (*p2)->momentum().py(), (*p2)->momentum().pz(), (*p2)->momentum().e());
0301               TLorentzVector ptot = pmu1 + pmu2;
0302               TVector3 booster = ptot.BoostVector();
0303               TLorentzVector leptdir = (((*p2)->pdg_id() > 0) ? pmu1 : pmu2);
0304 
0305               leptdir.Boost(-booster);
0306               pb0.Boost(-booster);
0307               hmumuMassSqr->Fill(ptot.M2());
0308               if (cos(leptdir.Vect().Angle(pb0.Vect())) > 0) {
0309                 hmumuMassSqrPlus->Fill(ptot.M2());
0310               } else {
0311                 hmumuMassSqrMinus->Fill(ptot.M2());
0312               }
0313             }
0314           }
0315         }
0316       }
0317     }
0318     // --------------------------------------------------------------
0319     // Calculate helicity angles to test polarization
0320     // (from Hrivnac et al., J. Phys. G21, 629)
0321 
0322     if ((*p)->pdg_id() == 5122 /* && endvert */) {  // lambdaB
0323 
0324       TLorentzVector pMuP;
0325       TLorentzVector pProt;
0326       TLorentzVector pLambda0;
0327       TLorentzVector pJpsi;
0328       TLorentzVector pLambdaB;
0329       TVector3 enne;
0330 
0331       nbs++;
0332       if (!endvert) {
0333         *undecayed << (*p)->pdg_id() << std::endl;
0334       } else {
0335         *decayed << (*p)->pdg_id() << " --> ";
0336         for (HepMC::GenVertex::particles_out_const_iterator bp = endvert->particles_out_const_begin();
0337              bp != endvert->particles_out_const_end();
0338              ++bp) {
0339           *decayed << (*bp)->pdg_id() << " ";
0340         }
0341       }
0342       *decayed << std::endl;
0343 
0344       pLambdaB.SetPxPyPzE((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0345       enne = -(pLambdaB.Vect().Cross(TVector3(0., 0., 1.))).Unit();
0346 
0347       if (endvert) {
0348         for (HepMC::GenVertex::particles_out_const_iterator p2 = endvert->particles_out_const_begin();
0349              p2 != endvert->particles_out_const_end();
0350              ++p2) {
0351           if ((*p2)->pdg_id() == 443) {  // J/psi
0352             pJpsi.SetPxPyPzE(
0353                 (*p2)->momentum().px(), (*p2)->momentum().py(), (*p2)->momentum().pz(), (*p2)->momentum().e());
0354             HepMC::GenVertex* psivert = (*p2)->end_vertex();
0355             if (psivert) {
0356               for (HepMC::GenVertex::particles_out_const_iterator p3 = psivert->particles_out_const_begin();
0357                    p3 != psivert->particles_out_const_end();
0358                    ++p3) {
0359                 if ((*p3)->pdg_id() == -13) {  // mu+
0360                   pMuP.SetPxPyPzE(
0361                       (*p3)->momentum().px(), (*p3)->momentum().py(), (*p3)->momentum().pz(), (*p3)->momentum().e());
0362                 }
0363               }
0364             }
0365           }
0366           if ((*p2)->pdg_id() == 3122) {  // Lambda0
0367             pLambda0.SetPxPyPzE(
0368                 (*p2)->momentum().px(), (*p2)->momentum().py(), (*p2)->momentum().pz(), (*p2)->momentum().e());
0369             HepMC::GenVertex* lamvert = (*p2)->end_vertex();
0370             if (lamvert) {
0371               for (HepMC::GenVertex::particles_out_const_iterator p3 = lamvert->particles_out_const_begin();
0372                    p3 != lamvert->particles_out_const_end();
0373                    ++p3) {
0374                 if (std::abs((*p3)->pdg_id()) == 2212) {  // p
0375                   pProt.SetPxPyPzE(
0376                       (*p3)->momentum().px(), (*p3)->momentum().py(), (*p3)->momentum().pz(), (*p3)->momentum().e());
0377                 }
0378               }
0379             }
0380           }
0381         }
0382       }
0383 
0384       TVector3 booster1 = pLambdaB.BoostVector();
0385       TVector3 booster2 = pLambda0.BoostVector();
0386       TVector3 booster3 = pJpsi.BoostVector();
0387 
0388       pLambda0.Boost(-booster1);
0389       pJpsi.Boost(-booster1);
0390       hCosThetaLambda->Fill(cos(pLambda0.Vect().Angle(enne)));
0391 
0392       pProt.Boost(-booster2);
0393       hCosTheta1->Fill(cos(pProt.Vect().Angle(pLambda0.Vect())));
0394       TVector3 tempY = (pLambda0.Vect().Cross(enne)).Unit();
0395       TVector3 xyProj = (enne.Dot(pProt.Vect())) * enne + (tempY.Dot(pProt.Vect())) * tempY;
0396       // find the sign of phi
0397       TVector3 crossProd = xyProj.Cross(enne);
0398       float tempPhi = (crossProd.Dot(pLambda0.Vect()) > 0 ? xyProj.Angle(enne) : -xyProj.Angle(enne));
0399       hPhi1->Fill(tempPhi);
0400 
0401       pMuP.Boost(-booster3);
0402       hCosTheta2->Fill(cos(pMuP.Vect().Angle(pJpsi.Vect())));
0403       tempY = (pJpsi.Vect().Cross(enne)).Unit();
0404       xyProj = (enne.Dot(pMuP.Vect())) * enne + (tempY.Dot(pMuP.Vect())) * tempY;
0405       // find the sign of phi
0406       crossProd = xyProj.Cross(enne);
0407       tempPhi = (crossProd.Dot(pJpsi.Vect()) > 0 ? xyProj.Angle(enne) : -xyProj.Angle(enne));
0408       hPhi2->Fill(tempPhi);
0409     }
0410   }
0411   // ---------------------------------------------------------
0412 
0413   hnB->Fill(nB);
0414   hnJpsi->Fill(nJpsi);
0415   hnBz->Fill(nBz);
0416   hnBzb->Fill(nBzb);
0417   // *undecayed << "-------------------------------" << std::endl;
0418   // *decayed << "-------------------------------" << std::endl;
0419 
0420   // if (nBz > 0) std::cout << "nBz = " << nBz << " nBz (K*mu+mu-) = " << nBzKmumu << " nMix = " << nBzmix << std::endl;
0421   // if (nBz > 0 && nBzKmumu == 0) Evt->print();
0422   // if (nB > 0) std::cout << "nB = " << nB << " nBz (JPsi Ks) = " << nBJpsiKs << std::endl;
0423   // if (nBp > 0) std::cout << "nB = " << nBp << " nBz (JPsi Kstar) = " << nBJpsiKstar << std::endl;
0424   // if (nBp > 0 && nBJpsiKstar == 0) Evt->print();
0425 
0426   return;
0427 }
0428 
0429 void EvtGenTestAnalyzer::endJob() {
0430   /*
0431   TObjArray Hlist(0);
0432   Hlist.Add(hGeneralId);
0433   Hlist.Add(hIdPhiDaugs);
0434   Hlist.Add(hIdJpsiMot);
0435   Hlist.Add(hnB);
0436   Hlist.Add(hnBz);
0437   Hlist.Add(hnBzb);
0438   Hlist.Add(hnJpsi);
0439   Hlist.Add(hMinvb);
0440   Hlist.Add(hPtbs);
0441   Hlist.Add(hPbs);
0442   Hlist.Add(hPhibs);
0443   Hlist.Add(hEtabs);
0444   Hlist.Add(hPtmu);
0445   Hlist.Add(hPmu);
0446   Hlist.Add(hPhimu);
0447   Hlist.Add(hEtamu);
0448   Hlist.Add(hPtRadPho);
0449   Hlist.Add(hPhiRadPho);
0450   Hlist.Add(hEtaRadPho);
0451   Hlist.Add(htbJpsiKs);
0452   Hlist.Add(htbbarJpsiKs);
0453   Hlist.Add(htbPlus);
0454   Hlist.Add(htbsUnmix);
0455   Hlist.Add(htbsMix);
0456   Hlist.Add(htbUnmix);
0457   Hlist.Add(htbMix);
0458   Hlist.Add(htbMixPlus);
0459   Hlist.Add(htbMixMinus);
0460   Hlist.Add(hmumuMassSqr);
0461   Hlist.Add(hmumuMassSqrPlus);
0462   Hlist.Add(hmumuMassSqrMinus);
0463   Hlist.Add(hIdBsDaugs);
0464   Hlist.Add(hIdBDaugs);
0465   Hlist.Add(hCosTheta1);
0466   Hlist.Add(hCosTheta2);
0467   Hlist.Add(hPhi1);
0468   Hlist.Add(hPhi2);
0469   Hlist.Add(hCosThetaLambda);
0470   Hlist.Write();
0471   fOutputFile->Close();
0472   */
0473   std::cout << "N_events = " << nevent << "\n";
0474   std::cout << "N_LambdaB = " << nbs << "\n";
0475   return;
0476 }
0477 
0478 DEFINE_FWK_MODULE(EvtGenTestAnalyzer);