File indexing completed on 2023-05-10 03:54:21
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <iostream>
0015 #include <fstream>
0016
0017 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0018
0019
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
0050
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
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;
0118
0119
0120 if (endvert && prodvert) {
0121 dectime = (endvert->position().t() - prodvert->position().t()) * mmcToPs / gamma;
0122
0123
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)
0146
0147 {
0148 htbPlus->Fill(dectime);
0149 }
0150
0151
0152
0153 if ((*p)->pdg_id() == 531 ) {
0154
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 ) {
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
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
0205
0206
0207
0208
0209
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) {
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) {
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) {
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) {
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) {
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
0292
0293
0294 if ((*p)->pdg_id() == 5122 ) {
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) {
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) {
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) {
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) {
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
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
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
0390
0391
0392 return;
0393 }
0394
0395 void EvtGenTestAnalyzer::endJob() {
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439 std::cout << "N_events = " << nevent << "\n";
0440 std::cout << "N_LambdaB = " << nbs << "\n";
0441 return;
0442 }
0443
0444 DEFINE_FWK_MODULE(EvtGenTestAnalyzer);