//Renata Kopecna #include #include #include #include #include #include #include using namespace fcnc; ///calculate decay angles according to theory convention http://arxiv.org/abs/0805.2525 void bu2kstarmumu_loader::theory_angles(const TLorentzVector& muplus, const TLorentzVector& muminus, const TLorentzVector& kaon, const TLorentzVector& pion, double& costhetal, double& costhetak, double& phi) { //set up boost vectors TLorentzVector b = muplus + muminus + kaon + pion; TLorentzVector mumu = muplus + muminus; TLorentzVector kpi = kaon + pion; TVector3 mumuboost(-mumu.BoostVector()); TVector3 kpiboost(-kpi.BoostVector()); TVector3 bboost(-b.BoostVector()); //determine costhetal TLorentzVector muminusd(muminus); muminusd.Boost(mumuboost); TLorentzVector bd(b); bd.Boost(mumuboost); costhetal = cos(muminusd.Vect().Angle(bd.Vect())); //determine costhetak TLorentzVector kaondd(kaon); kaondd.Boost(kpiboost); TLorentzVector bdd(b); bdd.Boost(kpiboost); costhetak = cos(kaondd.Vect().Angle(bdd.Vect())); //determine phi TLorentzVector kaonddd(kaon); kaonddd.Boost(bboost); TLorentzVector pionddd(pion); pionddd.Boost(bboost); TLorentzVector muminusddd(muminus); muminusddd.Boost(bboost); TLorentzVector muplusddd(muplus); muplusddd.Boost(bboost); TVector3 normalkpi = kaonddd.Vect().Cross(pionddd.Vect()); TVector3 normalmumu = muminusddd.Vect().Cross(muplusddd.Vect()); phi = normalkpi.Angle(normalmumu); TLorentzVector kpiddd(kpi); kpiddd.Boost(bboost); if (normalkpi.Cross(normalmumu).Dot(kpiddd.Vect()) < 0.0) phi = -phi; } ///calculate decay angles according to lhcb convention http://arxiv.org/abs/1304.6325, bzero is true for B0 and false for B0bar void bu2kstarmumu_loader::lhcb_angles(bool bplus, const TLorentzVector& muplus, const TLorentzVector& muminus, const TLorentzVector& kaon, const TLorentzVector& pion, double& costhetal, double& costhetak, double& phi){ //set up boost vectors TLorentzVector b = muplus + muminus + kaon + pion; TLorentzVector mumu = muplus + muminus; TLorentzVector kpi = kaon + pion; TVector3 mumuboost(-mumu.BoostVector()); TVector3 kpiboost(-kpi.BoostVector()); TVector3 bboost(-b.BoostVector()); //determine costhetal TLorentzVector muminusd(muminus); muminusd.Boost(mumuboost); TLorentzVector muplusd(muplus); muplusd.Boost(mumuboost); TLorentzVector bd(b); bd.Boost(mumuboost); if (bplus) costhetal = cos(muplusd.Vect().Angle(-bd.Vect())); else costhetal = cos(muminusd.Vect().Angle(-bd.Vect())); //determine costhetak TLorentzVector kaondd(kaon); kaondd.Boost(kpiboost); TLorentzVector bdd(b); bdd.Boost(kpiboost); costhetak = cos(kaondd.Vect().Angle(-bdd.Vect())); //determine phi TLorentzVector kaonddd(kaon); kaonddd.Boost(bboost); TLorentzVector pionddd(pion); pionddd.Boost(bboost); TLorentzVector muminusddd(muminus); muminusddd.Boost(bboost); TLorentzVector muplusddd(muplus); muplusddd.Boost(bboost); TVector3 normalkpi = kaonddd.Vect().Cross(pionddd.Vect()); TVector3 normalmumu = muplusddd.Vect().Cross(muminusddd.Vect()); TLorentzVector kpiddd(kpi); kpiddd.Boost(bboost); if (bplus) { phi = normalkpi.Angle(normalmumu); if ((normalmumu.Cross(normalkpi)).Dot(kpiddd.Vect()) < 0.0) phi = -phi; } else { phi = normalkpi.Angle(-normalmumu); if ((normalmumu.Cross(normalkpi)).Dot(kpiddd.Vect()) < 0.0) phi = -phi; } spdlog::trace("Calculated angles: \t ctl= {0:0.2f},\t ctk = {1:0.2f}, \t phi = {1:0.2f}", costhetal, costhetak ,phi); return; } //calculate angles (Heidelberg2018). Should be 100% identical to LHCb convention //definition as in LHCb-ANA-2013-097 page 23 void bu2kstarmumu_loader::hd_angles(bool bplus, const TLorentzVector& muplus, const TLorentzVector& muminus, const TLorentzVector& kaon, const TLorentzVector& pion, double& costhetal_lhs, double& costhetal_rhs, double& costhetak_lhs, double& costhetak_rhs, double& phi){ //4D and boost vector of B system TLorentzVector B = (muplus + muminus + kaon + pion); TVector3 boost_B = B.BoostVector(); //ctk TLorentzVector kaon_in_B = kaon; kaon_in_B.Boost(-boost_B); TLorentzVector pion_in_B = pion; pion_in_B.Boost(-boost_B); TLorentzVector kstar_in_B = kaon_in_B + pion_in_B; TLorentzVector kaon_in_kstar = kaon_in_B; kaon_in_kstar.Boost(-kstar_in_B.BoostVector()); costhetak_lhs = kaon_in_kstar.Vect().Unit().Dot(kstar_in_B.Vect().Unit()); TLorentzVector B_in_kstar = B; B_in_kstar.Boost(-boost_B); B_in_kstar.Boost(-kstar_in_B.BoostVector()); costhetak_rhs = kaon_in_kstar.Vect().Unit().Dot(-B_in_kstar.Vect().Unit()); //ctl TLorentzVector muplus_in_B = bplus ? muplus : muminus; muplus_in_B.Boost(-boost_B); TLorentzVector muminus_in_B = bplus ? muminus : muplus; muminus_in_B.Boost(-boost_B); TLorentzVector dimuon_in_B = muplus_in_B + muminus_in_B; TLorentzVector muon_in_dimuon = muplus_in_B; muon_in_dimuon.Boost(-dimuon_in_B.BoostVector()); costhetal_lhs = muon_in_dimuon.Vect().Unit().Dot(dimuon_in_B.Vect().Unit()); TLorentzVector B_in_dimuon = B; B_in_dimuon.Boost(-boost_B); B_in_dimuon.Boost(-dimuon_in_B.BoostVector()); costhetal_rhs = muon_in_dimuon.Vect().Unit().Dot(-B_in_dimuon.Vect().Unit()); //phi TVector3 normalkpi = (kaon_in_B.Vect().Unit().Cross(pion_in_B.Vect().Unit())).Unit(); TVector3 normalmumu = (muplus_in_B.Vect().Unit().Cross(muminus_in_B.Vect().Unit())).Unit(); phi = normalkpi.Angle(normalmumu); if ((bplus ? 1. : -1.)*(normalmumu.Cross(normalkpi)).Dot(kstar_in_B.Vect().Unit()) < 0.0) phi = -phi; } ///B+ case: no further cuts in the tuple. Read (up to) all events without restrictions std::vector bu2kstarmumu_loader::read_full_tuple(int year, std::string filename, std::string treename, bool isMC, bool mctruth, bool genLevelMC, int nevents) { //Set a year for gen lvl MC by hand if(genLevelMC) year = 2017; check_year(year); //open the root file TFile* root_file = new TFile(filename.c_str(), "READ"); assert(root_file->IsOpen()); //Read the tree TTree* tree = static_cast(root_file->Get(treename.c_str())); if (!tree){ //Check if the tree exists, if not, return an empty vector spdlog::warn("Could not get tree {0:s} from file {1:s}", treename, filename); std::vector blank_vector; return blank_vector; //assert(0); } std::vector events; event meas; Double_t mbu = PDGMASS_B; Double_t sigma_mbu = 0.0; //TODO: mbu possibly from the fit? Double_t mkstar = PDGMASS_K_STAR; Double_t muplus_px, muplus_py, muplus_pz, muplus_pe, muminus_px, muminus_py, muminus_pz, muminus_pe, muplus_px_TRUE, muplus_py_TRUE, muplus_pz_TRUE, muplus_pe_TRUE, muminus_px_TRUE, muminus_py_TRUE, muminus_pz_TRUE, muminus_pe_TRUE, kaon_px, kaon_py, kaon_pz, kaon_pe, pion_px, pion_py, pion_pz, pion_pe, kaon_pt, pion_pt, muplus_pt, muminus_pt, Ks_pi_plus_px, Ks_pi_plus_py, Ks_pi_plus_pz, Ks_pi_plus_pt, Ks_pi_plus_pe, Ks_pi_minus_px, Ks_pi_minus_py, Ks_pi_minus_pz, Ks_pi_minus_pt, Ks_pi_minus_pe, gamma1_px, gamma1_py, gamma1_pz, gamma1_pt, gamma1_pe, gamma2_px, gamma2_py, gamma2_pz, gamma2_pt, gamma2_pe, totalweight, delta_weight; Double_t ctl_TT, ctk_TT, phi_TT; Int_t b_id = 521, muplus_id = -13, muminus_id = +13;//, kaon_id, pion_id; Int_t DDLL; UInt_t runno = 0; ULong64_t gpstime = 0; Short_t polarity = 0; if(!genLevelMC){ tree->SetBranchAddress(opts->DTF ? "B_plus_M_DTF" : "B_plus_MM",&mbu); tree->SetBranchAddress(opts->DTF ? "B_plus_MERR_DTF" : "B_plus_MMERR",&sigma_mbu); tree->SetBranchAddress(opts->DTF ? "K_star_plus_M_scaled" : "K_star_plus_MM", &mkstar); tree->SetBranchAddress("Polarity", &polarity); tree->SetBranchAddress("runNumber", &runno); tree->SetBranchAddress("B_plus_ThetaL", &ctl_TT); tree->SetBranchAddress("B_plus_ThetaK", &ctk_TT); tree->SetBranchAddress("B_plus_Phi", &phi_TT); } if (!mctruth && !genLevelMC){ tree->SetBranchAddress("GpsTime", &gpstime); tree->SetBranchAddress((opts->DTF ? "mu_plus_PX_DTF" : "mu_plus_PX"), &muplus_px); tree->SetBranchAddress((opts->DTF ? "mu_plus_PY_DTF" : "mu_plus_PY"), &muplus_py); tree->SetBranchAddress((opts->DTF ? "mu_plus_PZ_DTF" : "mu_plus_PZ"), &muplus_pz); tree->SetBranchAddress((opts->DTF ? "mu_plus_PE_DTF" : "mu_plus_PE"), &muplus_pe); tree->SetBranchAddress((opts->DTF ? "mu_plus_PT_DTF" : "mu_plus_PT"), &muplus_pt); tree->SetBranchAddress((opts->DTF ? "mu_minus_PX_DTF" : "mu_minus_PX"), &muminus_px); tree->SetBranchAddress((opts->DTF ? "mu_minus_PY_DTF" : "mu_minus_PY"), &muminus_py); tree->SetBranchAddress((opts->DTF ? "mu_minus_PZ_DTF" : "mu_minus_PZ"), &muminus_pz); tree->SetBranchAddress((opts->DTF ? "mu_minus_PE_DTF" : "mu_minus_PE"), &muminus_pe); tree->SetBranchAddress((opts->DTF ? "mu_minus_PT_DTF" : "mu_minus_PT"), &muminus_pt); if(opts->KS){ //K_S^0 pi^+ channel if(opts->DTF){ tree->SetBranchAddress("Ks_pi_plus_PX_DTF", &Ks_pi_plus_px); tree->SetBranchAddress("Ks_pi_plus_PY_DTF", &Ks_pi_plus_py); tree->SetBranchAddress("Ks_pi_plus_PZ_DTF", &Ks_pi_plus_pz); tree->SetBranchAddress("Ks_pi_plus_PT_DTF", &Ks_pi_plus_pt); tree->SetBranchAddress("Ks_pi_plus_PE_DTF", &Ks_pi_plus_pe); tree->SetBranchAddress("Ks_pi_minus_PX_DTF", &Ks_pi_minus_px); tree->SetBranchAddress("Ks_pi_minus_PY_DTF", &Ks_pi_minus_py); tree->SetBranchAddress("Ks_pi_minus_PZ_DTF", &Ks_pi_minus_pz); tree->SetBranchAddress("Ks_pi_minus_PT_DTF", &Ks_pi_minus_pt); tree->SetBranchAddress("Ks_pi_minus_PE_DTF", &Ks_pi_minus_pe); } else{ tree->SetBranchAddress("K_short_PX", &kaon_px); tree->SetBranchAddress("K_short_PY", &kaon_py); tree->SetBranchAddress("K_short_PZ", &kaon_pz); tree->SetBranchAddress("K_short_PE", &kaon_pe); tree->SetBranchAddress("K_short_PT", &kaon_pt); } tree->SetBranchAddress((opts->DTF ? "pi_plus_PX_DTF" : "pi_plus_PX"), &pion_px); tree->SetBranchAddress((opts->DTF ? "pi_plus_PY_DTF" : "pi_plus_PY"), &pion_py); tree->SetBranchAddress((opts->DTF ? "pi_plus_PZ_DTF" : "pi_plus_PZ"), &pion_pz); tree->SetBranchAddress((opts->DTF ? "pi_plus_PE_DTF" : "pi_plus_PE"), &pion_pe); tree->SetBranchAddress((opts->DTF ? "pi_plus_PT_DTF" : "pi_plus_PT"), &pion_pt); } else{ //K^+pi^0 channel tree->SetBranchAddress((opts->DTF ? "K_plus_PX_DTF" : "K_short_PX"), &kaon_px); tree->SetBranchAddress((opts->DTF ? "K_plus_PY_DTF" : "K_short_PY"), &kaon_py); tree->SetBranchAddress((opts->DTF ? "K_plus_PZ_DTF" : "K_short_PZ"), &kaon_pz); tree->SetBranchAddress((opts->DTF ? "K_plus_PE_DTF" : "K_short_PE"), &kaon_pe); tree->SetBranchAddress((opts->DTF ? "K_plus_PT_DTF" : "K_short_PT"), &kaon_pt); if(opts->DTF){ tree->SetBranchAddress("gamma1_PX_DTF", &gamma1_px); tree->SetBranchAddress("gamma1_PY_DTF", &gamma1_py); tree->SetBranchAddress("gamma1_PZ_DTF", &gamma1_pz); tree->SetBranchAddress("gamma1_PE_DTF", &gamma1_pe); tree->SetBranchAddress("gamma1_PT_DTF", &gamma1_pt); tree->SetBranchAddress("gamma2_PX_DTF", &gamma2_px); tree->SetBranchAddress("gamma2_PY_DTF", &gamma2_py); tree->SetBranchAddress("gamma2_PZ_DTF", &gamma2_pz); tree->SetBranchAddress("gamma2_PE_DTF", &gamma2_pe); tree->SetBranchAddress("gamma2_PT_DTF", &gamma2_pt); } else{ tree->SetBranchAddress("pi_zero_resolved_PX", &pion_px); tree->SetBranchAddress("pi_zero_resolved_PY", &pion_py); tree->SetBranchAddress("pi_zero_resolved_PZ", &pion_pz); tree->SetBranchAddress("pi_zero_resolved_PE", &pion_pe); tree->SetBranchAddress("pi_zero_resolved_PT", &pion_pt); } } } else if(genLevelMC){ if(opts->KS){ tree->SetBranchAddress("KS0_TRUEP_X", &kaon_px); tree->SetBranchAddress("KS0_TRUEP_Y", &kaon_py); tree->SetBranchAddress("KS0_TRUEP_Z", &kaon_pz); tree->SetBranchAddress("KS0_TRUEP_E", &kaon_pe); tree->SetBranchAddress("piplus0_TRUEP_X", &pion_px); tree->SetBranchAddress("piplus0_TRUEP_Y", &pion_py); tree->SetBranchAddress("piplus0_TRUEP_Z", &pion_pz); tree->SetBranchAddress("piplus0_TRUEP_E", &pion_pe); tree->SetBranchAddress("piplus_TRUEP_X", &Ks_pi_plus_px); tree->SetBranchAddress("piplus_TRUEP_Y", &Ks_pi_plus_py); tree->SetBranchAddress("piplus_TRUEP_Z", &Ks_pi_plus_pz); tree->SetBranchAddress("piplus_TRUEP_E", &Ks_pi_plus_pe); tree->SetBranchAddress("piminus_TRUEP_X", &Ks_pi_minus_px); tree->SetBranchAddress("piminus_TRUEP_Y", &Ks_pi_minus_py); tree->SetBranchAddress("piminus_TRUEP_Z", &Ks_pi_minus_pz); tree->SetBranchAddress("piminus_TRUEP_E", &Ks_pi_minus_pe); } else{ tree->SetBranchAddress("Kplus_TRUEP_X", &kaon_px); tree->SetBranchAddress("Kplus_TRUEP_Y", &kaon_py); tree->SetBranchAddress("Kplus_TRUEP_Z", &kaon_pz); tree->SetBranchAddress("Kplus_TRUEP_E", &kaon_pe); tree->SetBranchAddress("pi0_TRUEP_X", &pion_px); tree->SetBranchAddress("pi0_TRUEP_Y", &pion_py); tree->SetBranchAddress("pi0_TRUEP_Z", &pion_pz); tree->SetBranchAddress("pi0_TRUEP_E", &pion_pe); tree->SetBranchAddress("muplus_TRUEP_X", &muplus_px_TRUE); tree->SetBranchAddress("muplus_TRUEP_Y", &muplus_py_TRUE); tree->SetBranchAddress("muplus_TRUEP_Z", &muplus_pz_TRUE); tree->SetBranchAddress("muplus_TRUEP_E", &muplus_pe_TRUE); tree->SetBranchAddress("muminus_TRUEP_X", &muminus_px_TRUE); tree->SetBranchAddress("muminus_TRUEP_Y", &muminus_py_TRUE); tree->SetBranchAddress("muminus_TRUEP_Z", &muminus_pz_TRUE); tree->SetBranchAddress("muminus_TRUEP_E", &muminus_pe_TRUE); } } else{ //use MC true information if(opts->KS){ tree->SetBranchAddress("K_short_TRUEP_X", &kaon_px); tree->SetBranchAddress("K_short_TRUEP_Y", &kaon_py); tree->SetBranchAddress("K_short_TRUEP_Z", &kaon_pz); tree->SetBranchAddress("K_short_TRUEP_E", &kaon_pe); tree->SetBranchAddress("pi_plus_TRUEP_X", &pion_px); tree->SetBranchAddress("pi_plus_TRUEP_Y", &pion_py); tree->SetBranchAddress("pi_plus_TRUEP_Z", &pion_pz); tree->SetBranchAddress("pi_plus_TRUEP_E", &pion_pe); } else{ tree->SetBranchAddress("K_plus_TRUEP_X", &kaon_px); tree->SetBranchAddress("K_plus_TRUEP_Y", &kaon_py); tree->SetBranchAddress("K_plus_TRUEP_Z", &kaon_pz); tree->SetBranchAddress("K_plus_TRUEP_E", &kaon_pe); tree->SetBranchAddress("pi_zero_resolved_TRUEP_X", &pion_px); tree->SetBranchAddress("pi_zero_resolved_TRUEP_Y", &pion_py); tree->SetBranchAddress("pi_zero_resolved_TRUEP_Z", &pion_pz); tree->SetBranchAddress("pi_zero_resolved_TRUEP_E", &pion_pe); } } //For any MC, always link the TRUE muon as well: if(isMC || mctruth){ tree->SetBranchAddress("mu_plus_TRUEP_X", &muplus_px_TRUE); tree->SetBranchAddress("mu_plus_TRUEP_Y", &muplus_py_TRUE); tree->SetBranchAddress("mu_plus_TRUEP_Z", &muplus_pz_TRUE); tree->SetBranchAddress("mu_plus_TRUEP_E", &muplus_pe_TRUE); tree->SetBranchAddress("mu_minus_TRUEP_X", &muminus_px_TRUE); tree->SetBranchAddress("mu_minus_TRUEP_Y", &muminus_py_TRUE); tree->SetBranchAddress("mu_minus_TRUEP_Z", &muminus_pz_TRUE); tree->SetBranchAddress("mu_minus_TRUEP_E", &muminus_pe_TRUE); } //particle IDs: if(genLevelMC){ tree->SetBranchAddress("Kst_892_plus_MC_MOTHER_ID", &b_id); } else{ tree->SetBranchAddress("B_plus_ID", &b_id); tree->SetBranchAddress("mu_minus_ID", &muminus_id); tree->SetBranchAddress("mu_plus_ID", &muplus_id); } //for K_S^0 channel, get DD or LL track information if(opts->KS && !genLevelMC)tree->SetBranchAddress("KshortDecayInVeLo", &DDLL); //load either sWeights (data) or MC reweights (signal, reference and phase-space MC) if(!genLevelMC){ if(isMC){ if (opts->KS){ tree->SetBranchAddress("weight2D_nLongTracks", &totalweight); tree->SetBranchAddress("delta_weight2D_nLongTracks", &delta_weight); } else{ tree->SetBranchAddress("weight2D_nLongTracks_TM_rndGamma", &totalweight); tree->SetBranchAddress("delta_weight2D_nLongTracks_TM_rndGamma", &delta_weight); } } else tree->SetBranchAddress("N_Bplus_sw", &totalweight); } //determine number of events to load: unsigned int nentries = tree->GetEntries(); if (nevents > 0){ unsigned int n = TMath::Abs(nevents); if(n < nentries)nentries = n; } spdlog::info("Reading {0:d} events from file {1:s}", nentries ,filename); bool DDtracks = strcmp(opts->DDLL.c_str(), "DD") == 0; std::vectorvCorruptedEvents; vCorruptedEvents.clear(); //loop over input tree: for(unsigned int j=0; j < nentries; j++){ tree->GetEntry(j); if (j % (tree->GetEntries()/10) == 0) std::cout << "\r[" << round(double(j*100)/tree->GetEntries()) << "%]" << std::flush; //Is it a B+ or a B-? bool b_plus = b_id > 0; //for genLevelMonteCarlo: define muonID accordingly to B_id: if(genLevelMC){ muminus_id = b_plus ? +13 : -13; muplus_id = b_plus ? -13 : +13; } //skip events in K_S^0 sample that have incorrect tracktype if(opts->KS && !genLevelMC){ if( DDtracks && DDLL)continue; if(!DDtracks && !DDLL)continue; } if(genLevelMC || mctruth){ //MC TRUE samples and GenLevel samples use the MC true information: muminus_px = muminus_px_TRUE; muminus_py = muminus_py_TRUE; muminus_pz = muminus_pz_TRUE; muminus_pe = muminus_pe_TRUE; muplus_px = muplus_px_TRUE; muplus_py = muplus_py_TRUE; muplus_pz = muplus_pz_TRUE; muplus_pe = muplus_pe_TRUE; } //initialize 4D vectors for angel calculations TLorentzVector muplus, muminus, kaon, pion, muplus_TRUE, muminus_TRUE; if(muminus_id == muplus_id){ spdlog::critical("Cannot have mu+ and mu- to have the same ID: {0:d}", muplus_id); spdlog::critical("Event ID = {0:d}"); assert(0); } if (muminus_id == +13){ muminus = TLorentzVector(muminus_px, muminus_py, muminus_pz, muminus_pe); muminus_TRUE = TLorentzVector(muminus_px_TRUE, muminus_py_TRUE, muminus_pz_TRUE, muminus_pe_TRUE); } else if(muminus_id == -13){ muplus = TLorentzVector(muminus_px, muminus_py, muminus_pz, muminus_pe); muplus_TRUE = TLorentzVector(muminus_px_TRUE, muminus_py_TRUE, muminus_pz_TRUE, muminus_pe_TRUE); } else{ spdlog::critical("Found invalid muon ID in mu-: {0:d}", muminus_id); assert(0); } if (muplus_id == -13){ muplus = TLorentzVector(muplus_px, muplus_py, muplus_pz, muplus_pe); muplus_TRUE = TLorentzVector(muplus_px_TRUE, muplus_py_TRUE, muplus_pz_TRUE, muplus_pe_TRUE); } else if(muplus_id == +13){ muminus = TLorentzVector(muplus_px, muplus_py, muplus_pz, muplus_pe); muminus_TRUE = TLorentzVector(muplus_px_TRUE, muplus_py_TRUE, muplus_pz_TRUE, muplus_pe_TRUE); } else{ spdlog::critical("Found invalid muon ID in mu-: {0:d}", muplus_id); assert(0); } if(opts->DTF && !mctruth && !genLevelMC){ //If no true information, re-reconstruct the pion+kaon pair if(opts->KS){ pion = TLorentzVector(pion_px, pion_py, pion_pz, pion_pe); kaon = TLorentzVector(Ks_pi_plus_px, Ks_pi_plus_py, Ks_pi_plus_pz, Ks_pi_plus_pe) + TLorentzVector(Ks_pi_minus_px, Ks_pi_minus_py, Ks_pi_minus_pz, Ks_pi_minus_pe); } else{ pion = TLorentzVector(gamma1_px, gamma1_py, gamma1_pz, gamma1_pe) + TLorentzVector(gamma2_px, gamma2_py, gamma2_pz, gamma2_pe); kaon = TLorentzVector(kaon_px, kaon_py, kaon_pz, kaon_pe); } } else{ //in case of true info, read kaon and pion directly pion = TLorentzVector(pion_px, pion_py, pion_pz, pion_pe); kaon = TLorentzVector(kaon_px, kaon_py, kaon_pz, kaon_pe); } //In case of genLvl MC recalculate B mass if(genLevelMC) mbu = (kaon+pion+muplus+muminus).M(); //de-bugging if (opts->KS){ spdlog::trace("KS_piplus_PX: {0:f}", Ks_pi_plus_px ); spdlog::trace("KS_piplus_PY: {0:f}", Ks_pi_plus_py ); spdlog::trace("KS_piplus_PZ: {0:f}", Ks_pi_plus_pz ); spdlog::trace("KS_piplus_PE: {0:f}", Ks_pi_plus_pe ); spdlog::trace("KS_piminus_PX: {0:f}", Ks_pi_minus_px); spdlog::trace("KS_piminus_PY: {0:f}", Ks_pi_minus_py); spdlog::trace("KS_piminus_PZ: {0:f}", Ks_pi_minus_pz); spdlog::trace("KS_piminus_PE: {0:f}", Ks_pi_minus_pe); } //print vectors if(spdlog::default_logger_raw()->level()==spdlog::level::trace){ spdlog::debug(b_plus ? "[B+]" : "[B-]"); muplus.Print(); muminus.Print(); pion.Print(); kaon.Print(); } //make sure that the TRUE_P's of muons in MC are not identical if(isMC && opts->reject_identical_muon_TRUEID){ if(muplus_TRUE.X() == muminus_TRUE.X() && muplus_TRUE.Y() == muminus_TRUE.Y() && muplus_TRUE.Z() == muminus_TRUE.Z() && muplus_TRUE.T() == muminus_TRUE.T()){ vCorruptedEvents.push_back(j); continue; } } //make sure that mup != mun if(muplus.X() == muminus.X() && muplus.Y() == muminus.Y() && muplus.Z() == muminus.Z() && muplus.T() == muminus.T()){ vCorruptedEvents.push_back(j); continue; } //get pseudo_q2 from opening angle of leptons: Double_t pseudo_q2 = muplus.Angle(muminus.Vect()); // Double_t pseudo_q2 = std::max(muminus.Pt(), muplus.Pt()); //calculate angles Double_t ctl_theory, ctk_theory, phi_theory; Double_t ctl_lhcb, ctk_lhcb, phi_lhcb; Double_t ctl_converted, ctk_converted, phi_converted; Double_t ctl_hd_lhs, ctk_hd_lhs, phi_hd; Double_t ctl_hd_rhs, ctk_hd_rhs; theory_angles( muplus, muminus, kaon, pion, ctl_theory, ctk_theory, phi_theory); lhcb_angles (b_plus, muplus, muminus, kaon, pion, ctl_lhcb, ctk_lhcb, phi_lhcb); hd_angles (b_plus, muplus, muminus, kaon, pion, ctl_hd_lhs, ctl_hd_rhs, ctk_hd_rhs, ctk_hd_lhs, phi_hd); //check and convert theory convention to LHCb convention ctl_converted = b_plus ? ctl_theory : -ctl_theory; ctk_converted = -ctk_theory; if (b_plus) phi_converted = phi_theory > 0 ? TMath::Pi() - phi_theory : -TMath::Pi() - phi_theory; else phi_converted = phi_theory; //convert tupletoolP2VV angles to cos() ctl_TT = cos(ctl_TT); ctk_TT = cos(ctk_TT); spdlog::trace("Event #{0:d}", j); spdlog::trace("B" + std::string(b_plus ? "+" : "-") + "\t\t: Theory \tLHCb \t\tConverted \tHD(lhs) \tHD(rhs)\tTupleToolP2VV"); spdlog::trace("cos(theta_L)\t:{0:0.4f}\t {0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t", ctl_theory, ctl_lhcb, ctl_converted, ctl_hd_lhs, ctl_hd_rhs, ctl_TT); spdlog::trace("cos(theta_k)\t:{0:0.4f}\t {0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t", ctk_theory, ctk_lhcb , ctk_converted, ctk_hd_lhs, ctk_hd_rhs, ctk_TT); spdlog::trace("phit:{0:0.4f}\t {0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t{0:0.4f}\t\t\n", phi_theory, phi_lhcb, phi_converted, phi_hd, phi_TT); //test all angles to be identical: assert(fabs(ctl_lhcb - ctl_converted) < 1.0e-8); assert(fabs(ctk_lhcb - ctk_converted) < 1.0e-8); assert(fabs(phi_lhcb - phi_converted) < 1.0e-8); assert(fabs(ctl_lhcb - ctl_hd_lhs) < 1.0e-8); assert(fabs(ctl_lhcb - ctl_hd_rhs) < 1.0e-8); assert(fabs(ctk_lhcb - ctk_hd_lhs) < 1.0e-8); assert(fabs(ctk_lhcb - ctk_hd_rhs) < 1.0e-8); assert(fabs(phi_lhcb - phi_hd) < 1.0e-8); meas.year = year; meas.m = mbu; meas.sigma_m = sigma_mbu; meas.costhetal = ctl_lhcb; meas.costhetak = ctk_lhcb; meas.phi = phi_lhcb; meas.costhetal_fa = ctl_lhcb; meas.costhetak_fa = ctk_lhcb; meas.phi_fa = phi_lhcb; meas.q2 = (muplus+muminus).Mag2()/1.0e+6;//in gev2 meas.p2 = (kaon+pion).Mag2()/1.0e+6;//in gev2 meas.weight = totalweight; meas.delta_weight = delta_weight; meas.mkpi = mkstar; meas.cp_conjugate = b_plus;//this is B+ meas.max_mu_pt = TMath::Max(muplus.Pt(), muminus.Pt()); meas.kaon_pt = kaon.Pt(); meas.pseudo_q2 = pseudo_q2; meas.magnet = polarity; bool ApplyGenLevelTightCuts = false; if(genLevelMC && ApplyGenLevelTightCuts){ if(muplus.Theta() < 0.005 || muplus.Theta() > 0.4)continue; if(muminus.Theta() < 0.005 || muminus.Theta() > 0.4)continue; if(pion.Theta() < 0.005 || pion.Theta() > 0.4)continue; TLorentzVector piplus(Ks_pi_plus_px, Ks_pi_plus_py, Ks_pi_plus_pz, Ks_pi_plus_pe); TLorentzVector piminus(Ks_pi_minus_px, Ks_pi_minus_py, Ks_pi_minus_pz, Ks_pi_minus_pe); if(piplus.P() <= 1800)continue; if(piminus.P() <= 1800)continue; if(piplus.Theta() < 0.005 || piplus.Theta() > 0.4)continue; if(piminus.Theta() < 0.005 || piminus.Theta() > 0.4)continue; } events.push_back(meas); } if(vCorruptedEvents.size() > 0){ spdlog::warn("Found {0:d} corrupted events with idendical muons: ",vCorruptedEvents.size()); for(UInt_t c = 0; c < vCorruptedEvents.size(); c++) spdlog::warn("Events at: {0:d}",vCorruptedEvents.at(c)); } root_file->Close(); delete root_file; spdlog::info("DONE reading {0:d} events from file " + filename,events.size()); return events; }