Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-05-10 03:54:21

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 nBz = 0;
0108   int nBzb = 0;
0109 
0110   for (HepMC::GenEvent::particle_const_iterator p = Evt->particles_begin(); p != Evt->particles_end(); ++p) {
0111     // General
0112     TLorentzVector thePart4m(0., 0., 0., 0.);
0113     HepMC::GenVertex* endvert = (*p)->end_vertex();
0114     HepMC::GenVertex* prodvert = (*p)->production_vertex();
0115     float gamma = (*p)->momentum().e() / (*p)->generated_mass();
0116     float dectime = 0.0;
0117     int mixed = -1;  // mixed is: -1 = unmixed
0118                      //           0 = mixed (before mixing)
0119                      //           1 = mixed (after mixing)
0120     if (endvert && prodvert) {
0121       dectime = (endvert->position().t() - prodvert->position().t()) * mmcToPs / gamma;
0122 
0123       // Mixed particle ?
0124       for (HepMC::GenVertex::particles_in_const_iterator p2 = prodvert->particles_in_const_begin();
0125            p2 != prodvert->particles_in_const_end();
0126            ++p2) {
0127         if ((*p)->pdg_id() + (*p2)->pdg_id() == 0) {
0128           mixed = 1;
0129           gamma = (*p2)->momentum().e() / (*p2)->generated_mass();
0130           HepMC::GenVertex* mixvert = (*p2)->production_vertex();
0131           dectime = (prodvert->position().t() - mixvert->position().t()) * mmcToPs / gamma;
0132         }
0133       }
0134       for (HepMC::GenVertex::particles_out_const_iterator ap = endvert->particles_out_const_begin();
0135            ap != endvert->particles_out_const_end();
0136            ++ap) {
0137         if ((*p)->pdg_id() + (*ap)->pdg_id() == 0)
0138           mixed = 0;
0139       }
0140     }
0141 
0142     hGeneralId->Fill((*p)->pdg_id());
0143 
0144     // --------------------------------------------------------------
0145     if (std::abs((*p)->pdg_id()) == 521)  // B+/-
0146     // || abs((*p)->pdg_id()/100) == 4 || abs((*p)->pdg_id()/100) == 3)
0147     {
0148       htbPlus->Fill(dectime);
0149     }
0150     // if ((*p)->pdg_id() == 443) *undecayed << (*p)->pdg_id() << std::endl;
0151     // --------------------------------------------------------------
0152 
0153     if ((*p)->pdg_id() == 531 /* && endvert */) {  // B_s
0154       // nbs++;
0155       hPtbs->Fill((*p)->momentum().perp());
0156       hPbs->Fill(sqrt(pow((*p)->momentum().px(), 2) + pow((*p)->momentum().py(), 2) + pow((*p)->momentum().pz(), 2)));
0157       hPhibs->Fill((*p)->momentum().phi());
0158       hEtabs->Fill((*p)->momentum().pseudoRapidity());
0159 
0160       for (HepMC::GenVertex::particles_out_const_iterator ap = endvert->particles_out_const_begin();
0161            ap != endvert->particles_out_const_end();
0162            ++ap) {
0163         hIdBsDaugs->Fill((*ap)->pdg_id());
0164       }
0165 
0166       if (mixed == 1) {
0167         htbsMix->Fill(dectime);
0168       } else if (mixed == -1) {
0169         htbsUnmix->Fill(dectime);
0170       }
0171     }
0172     // --------------------------------------------------------------
0173     if (std::abs((*p)->pdg_id()) == 511 /* && endvert */) {  // B0
0174       if (mixed != 0) {
0175         nB++;
0176         if ((*p)->pdg_id() > 0) {
0177           nBz++;
0178         } else {
0179           nBzb++;
0180         }
0181       }
0182       int isJpsiKs = 0;
0183       int isSemilept = 0;
0184       for (HepMC::GenVertex::particles_out_const_iterator bp = endvert->particles_out_const_begin();
0185            bp != endvert->particles_out_const_end();
0186            ++bp) {
0187         // Check invariant mass consistency ...
0188         TLorentzVector theDaug4m(
0189             (*bp)->momentum().px(), (*bp)->momentum().py(), (*bp)->momentum().pz(), (*bp)->momentum().e());
0190         thePart4m += theDaug4m;
0191         hIdBDaugs->Fill((*bp)->pdg_id());
0192         if ((*bp)->pdg_id() == 443 || (*bp)->pdg_id() == 310)
0193           isJpsiKs++;
0194         if ((*bp)->pdg_id() == 22) {
0195           hPtRadPho->Fill((*bp)->momentum().perp());
0196           hPhiRadPho->Fill((*bp)->momentum().phi());
0197           hEtaRadPho->Fill((*bp)->momentum().pseudoRapidity());
0198         }
0199         if (std::abs((*bp)->pdg_id()) == 11 || std::abs((*bp)->pdg_id()) == 13 || std::abs((*bp)->pdg_id()) == 15)
0200           isSemilept++;
0201       }
0202 
0203       hMinvb->Fill(sqrt(thePart4m.M2()));
0204       /* if (fabs(sqrt(thePart4m.M2())-5.28) > 0.05) {
0205      *undecayed << sqrt(thePart4m.M2()) << "  " << (*p)->pdg_id() << " --> " ;
0206          for ( GenVertex::particles_out_const_iterator bp = endvert->particles_out_const_begin(); bp != endvert->particles_out_const_end(); ++bp ) {
0207        *undecayed << (*bp)->pdg_id() << " " ;
0208          }
0209          *undecayed << std::endl;
0210        } */
0211 
0212       if (isSemilept) {
0213         if (mixed == 1) {
0214           htbMix->Fill(dectime);
0215           if ((*p)->pdg_id() > 0) {
0216             htbMixPlus->Fill(dectime);
0217           } else {
0218             htbMixMinus->Fill(dectime);
0219           }
0220         } else if (mixed == -1) {
0221           htbUnmix->Fill(dectime);
0222         }
0223       }
0224 
0225       if (isJpsiKs == 2) {
0226         if ((*p)->pdg_id() * mixed < 0) {
0227           htbbarJpsiKs->Fill(dectime);
0228         } else {
0229           htbJpsiKs->Fill(dectime);
0230         }
0231       }
0232     }
0233     // --------------------------------------------------------------
0234     if ((*p)->pdg_id() == 443) {  //Jpsi
0235       nJpsi++;
0236       for (HepMC::GenVertex::particles_in_const_iterator ap = prodvert->particles_in_const_begin();
0237            ap != prodvert->particles_in_const_end();
0238            ++ap) {
0239         hIdJpsiMot->Fill((*ap)->pdg_id());
0240       }
0241     }
0242     // --------------------------------------------------------------
0243     if ((*p)->pdg_id() == 333) {  // phi
0244       if (endvert) {
0245         for (HepMC::GenVertex::particles_out_const_iterator cp = endvert->particles_out_const_begin();
0246              cp != endvert->particles_out_const_end();
0247              ++cp) {
0248           hIdPhiDaugs->Fill((*cp)->pdg_id());
0249         }
0250       }
0251     }
0252     // --------------------------------------------------------------
0253     if ((*p)->pdg_id() == 13) {  // mu+
0254       for (HepMC::GenVertex::particles_in_const_iterator p2 = prodvert->particles_in_const_begin();
0255            p2 != prodvert->particles_in_const_end();
0256            ++p2) {
0257         if (std::abs((*p2)->pdg_id()) == 511) {  // B0
0258           hPtmu->Fill((*p)->momentum().perp());
0259           hPmu->Fill(
0260               sqrt(pow((*p)->momentum().px(), 2) + pow((*p)->momentum().py(), 2) + pow((*p)->momentum().pz(), 2)));
0261           hPhimu->Fill((*p)->momentum().phi());
0262           hEtamu->Fill((*p)->momentum().pseudoRapidity());
0263           for (HepMC::GenVertex::particles_out_const_iterator p3 = prodvert->particles_out_const_begin();
0264                p3 != prodvert->particles_out_const_end();
0265                ++p3) {
0266             if ((*p3)->pdg_id() == -13) {  // mu-
0267               TLorentzVector pmu1(
0268                   (*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0269               TLorentzVector pmu2(
0270                   (*p3)->momentum().px(), (*p3)->momentum().py(), (*p3)->momentum().pz(), (*p3)->momentum().e());
0271               TLorentzVector pb0(
0272                   (*p2)->momentum().px(), (*p2)->momentum().py(), (*p2)->momentum().pz(), (*p2)->momentum().e());
0273               TLorentzVector ptot = pmu1 + pmu2;
0274               TVector3 booster = ptot.BoostVector();
0275               TLorentzVector leptdir = (((*p2)->pdg_id() > 0) ? pmu1 : pmu2);
0276 
0277               leptdir.Boost(-booster);
0278               pb0.Boost(-booster);
0279               hmumuMassSqr->Fill(ptot.M2());
0280               if (cos(leptdir.Vect().Angle(pb0.Vect())) > 0) {
0281                 hmumuMassSqrPlus->Fill(ptot.M2());
0282               } else {
0283                 hmumuMassSqrMinus->Fill(ptot.M2());
0284               }
0285             }
0286           }
0287         }
0288       }
0289     }
0290     // --------------------------------------------------------------
0291     // Calculate helicity angles to test polarization
0292     // (from Hrivnac et al., J. Phys. G21, 629)
0293 
0294     if ((*p)->pdg_id() == 5122 /* && endvert */) {  // lambdaB
0295 
0296       TLorentzVector pMuP;
0297       TLorentzVector pProt;
0298       TLorentzVector pLambda0;
0299       TLorentzVector pJpsi;
0300       TLorentzVector pLambdaB;
0301       TVector3 enne;
0302 
0303       nbs++;
0304       if (!endvert) {
0305         *undecayed << (*p)->pdg_id() << std::endl;
0306       } else {
0307         *decayed << (*p)->pdg_id() << " --> ";
0308         for (HepMC::GenVertex::particles_out_const_iterator bp = endvert->particles_out_const_begin();
0309              bp != endvert->particles_out_const_end();
0310              ++bp) {
0311           *decayed << (*bp)->pdg_id() << " ";
0312         }
0313       }
0314       *decayed << std::endl;
0315 
0316       pLambdaB.SetPxPyPzE((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz(), (*p)->momentum().e());
0317       enne = -(pLambdaB.Vect().Cross(TVector3(0., 0., 1.))).Unit();
0318 
0319       if (endvert) {
0320         for (HepMC::GenVertex::particles_out_const_iterator p2 = endvert->particles_out_const_begin();
0321              p2 != endvert->particles_out_const_end();
0322              ++p2) {
0323           if ((*p2)->pdg_id() == 443) {  // J/psi
0324             pJpsi.SetPxPyPzE(
0325                 (*p2)->momentum().px(), (*p2)->momentum().py(), (*p2)->momentum().pz(), (*p2)->momentum().e());
0326             HepMC::GenVertex* psivert = (*p2)->end_vertex();
0327             if (psivert) {
0328               for (HepMC::GenVertex::particles_out_const_iterator p3 = psivert->particles_out_const_begin();
0329                    p3 != psivert->particles_out_const_end();
0330                    ++p3) {
0331                 if ((*p3)->pdg_id() == -13) {  // mu+
0332                   pMuP.SetPxPyPzE(
0333                       (*p3)->momentum().px(), (*p3)->momentum().py(), (*p3)->momentum().pz(), (*p3)->momentum().e());
0334                 }
0335               }
0336             }
0337           }
0338           if ((*p2)->pdg_id() == 3122) {  // Lambda0
0339             pLambda0.SetPxPyPzE(
0340                 (*p2)->momentum().px(), (*p2)->momentum().py(), (*p2)->momentum().pz(), (*p2)->momentum().e());
0341             HepMC::GenVertex* lamvert = (*p2)->end_vertex();
0342             if (lamvert) {
0343               for (HepMC::GenVertex::particles_out_const_iterator p3 = lamvert->particles_out_const_begin();
0344                    p3 != lamvert->particles_out_const_end();
0345                    ++p3) {
0346                 if (std::abs((*p3)->pdg_id()) == 2212) {  // p
0347                   pProt.SetPxPyPzE(
0348                       (*p3)->momentum().px(), (*p3)->momentum().py(), (*p3)->momentum().pz(), (*p3)->momentum().e());
0349                 }
0350               }
0351             }
0352           }
0353         }
0354       }
0355 
0356       TVector3 booster1 = pLambdaB.BoostVector();
0357       TVector3 booster2 = pLambda0.BoostVector();
0358       TVector3 booster3 = pJpsi.BoostVector();
0359 
0360       pLambda0.Boost(-booster1);
0361       pJpsi.Boost(-booster1);
0362       hCosThetaLambda->Fill(cos(pLambda0.Vect().Angle(enne)));
0363 
0364       pProt.Boost(-booster2);
0365       hCosTheta1->Fill(cos(pProt.Vect().Angle(pLambda0.Vect())));
0366       TVector3 tempY = (pLambda0.Vect().Cross(enne)).Unit();
0367       TVector3 xyProj = (enne.Dot(pProt.Vect())) * enne + (tempY.Dot(pProt.Vect())) * tempY;
0368       // find the sign of phi
0369       TVector3 crossProd = xyProj.Cross(enne);
0370       float tempPhi = (crossProd.Dot(pLambda0.Vect()) > 0 ? xyProj.Angle(enne) : -xyProj.Angle(enne));
0371       hPhi1->Fill(tempPhi);
0372 
0373       pMuP.Boost(-booster3);
0374       hCosTheta2->Fill(cos(pMuP.Vect().Angle(pJpsi.Vect())));
0375       tempY = (pJpsi.Vect().Cross(enne)).Unit();
0376       xyProj = (enne.Dot(pMuP.Vect())) * enne + (tempY.Dot(pMuP.Vect())) * tempY;
0377       // find the sign of phi
0378       crossProd = xyProj.Cross(enne);
0379       tempPhi = (crossProd.Dot(pJpsi.Vect()) > 0 ? xyProj.Angle(enne) : -xyProj.Angle(enne));
0380       hPhi2->Fill(tempPhi);
0381     }
0382   }
0383   // ---------------------------------------------------------
0384 
0385   hnB->Fill(nB);
0386   hnJpsi->Fill(nJpsi);
0387   hnBz->Fill(nBz);
0388   hnBzb->Fill(nBzb);
0389   // *undecayed << "-------------------------------" << std::endl;
0390   // *decayed << "-------------------------------" << std::endl;
0391 
0392   return;
0393 }
0394 
0395 void EvtGenTestAnalyzer::endJob() {
0396   /*
0397   TObjArray Hlist(0);
0398   Hlist.Add(hGeneralId);
0399   Hlist.Add(hIdPhiDaugs);
0400   Hlist.Add(hIdJpsiMot);
0401   Hlist.Add(hnB);
0402   Hlist.Add(hnBz);
0403   Hlist.Add(hnBzb);
0404   Hlist.Add(hnJpsi);
0405   Hlist.Add(hMinvb);
0406   Hlist.Add(hPtbs);
0407   Hlist.Add(hPbs);
0408   Hlist.Add(hPhibs);
0409   Hlist.Add(hEtabs);
0410   Hlist.Add(hPtmu);
0411   Hlist.Add(hPmu);
0412   Hlist.Add(hPhimu);
0413   Hlist.Add(hEtamu);
0414   Hlist.Add(hPtRadPho);
0415   Hlist.Add(hPhiRadPho);
0416   Hlist.Add(hEtaRadPho);
0417   Hlist.Add(htbJpsiKs);
0418   Hlist.Add(htbbarJpsiKs);
0419   Hlist.Add(htbPlus);
0420   Hlist.Add(htbsUnmix);
0421   Hlist.Add(htbsMix);
0422   Hlist.Add(htbUnmix);
0423   Hlist.Add(htbMix);
0424   Hlist.Add(htbMixPlus);
0425   Hlist.Add(htbMixMinus);
0426   Hlist.Add(hmumuMassSqr);
0427   Hlist.Add(hmumuMassSqrPlus);
0428   Hlist.Add(hmumuMassSqrMinus);
0429   Hlist.Add(hIdBsDaugs);
0430   Hlist.Add(hIdBDaugs);
0431   Hlist.Add(hCosTheta1);
0432   Hlist.Add(hCosTheta2);
0433   Hlist.Add(hPhi1);
0434   Hlist.Add(hPhi2);
0435   Hlist.Add(hCosThetaLambda);
0436   Hlist.Write();
0437   fOutputFile->Close();
0438   */
0439   std::cout << "N_events = " << nevent << "\n";
0440   std::cout << "N_LambdaB = " << nbs << "\n";
0441   return;
0442 }
0443 
0444 DEFINE_FWK_MODULE(EvtGenTestAnalyzer);