Angular analysis of B+->K*+(K+pi0)mumu
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 

836 lines
37 KiB

#include "../GlobalFunctions.hh"
//////////////////////////////////////////////////////////////////////
/// OptimizePreSelection()
/// loops over all .root-files created by Preselection() defined in BDTSelection.cc.
/// More kinematic, PID and some more cuts are applied and saved in a new file.
///
THStack *stack(vector<TH1F*> histos){
THStack *tmp = new THStack("IncMC","IncMC");
for(unsigned int h = 0; h < histos.size(); h++){
histos.at(h)->SetFillColor(2+h);
tmp->Add(histos.at(h));
}
return tmp;
}
int plotIncSample(std::string year = "2012", std::string magnet = "down") {
bool ReferenceChannel = true;
TFile* input = 0;
TFile* histos = 0;
TTree *tree;
// Kst2Kpluspi0Resolved case
if(Kst2Kpluspi0Resolved){
input = TFile::Open(Form("%s/data/MC/BtoXJpsi/%s%s/%s%s_pi0Resolved%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),(smallSample ? "_small":"") ));
histos = TFile::Open(Form("%s/data/MC/BtoXJpsi/%s%s/%s%s_pi0Resolved_histos%s.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),(smallSample ? "_small":"") ));
}
// Kst2Kpluspi0Merged case
if(Kst2Kpluspi0Merged){
//TODO
}
// Kst2Kspiplus case
if(Kst2Kspiplus){
//TODO
}
input->cd();
cout << "In input" << endl;
tree = (TTree*)input->Get("DecayTree");
cout << "Tree has " << tree->GetEntries() << " entries" << endl;
tree->SetBranchStatus("*",1);
Double_t B_plus_M_DTF;
//assign variables to TRUEID
Int_t B_plus_TRUEID;
Int_t K_star_plus_TRUEID;
Int_t K_plus_TRUEID;
Int_t K_short_TRUEID;
Int_t Ks_pi_plus_TRUEID;
Int_t Ks_pi_minus_TRUEID;
Int_t pi_plus_TRUEID;
Int_t pi_zero_TRUEID;
Int_t mu_plus_TRUEID;
Int_t mu_minus_TRUEID;
Int_t J_psi_TRUEID;
Int_t gamma1_TRUEID;
Int_t gamma2_TRUEID;
//assign variables to TRUE_MOTHER_ID
Int_t B_plus_MOTHER_ID;
Int_t K_star_plus_MOTHER_ID;
Int_t K_plus_MOTHER_ID;
Int_t K_short_MOTHER_ID;
Int_t Ks_pi_plus_MOTHER_ID;
Int_t Ks_pi_minus_MOTHER_ID;
Int_t pi_plus_MOTHER_ID;
Int_t pi_zero_MOTHER_ID;
Int_t mu_plus_MOTHER_ID;
Int_t mu_minus_MOTHER_ID;
Int_t J_psi_MOTHER_ID;
Int_t gamma1_MOTHER_ID;
Int_t gamma2_MOTHER_ID;
//assign variables to TRUE_GD_MOTHER_ID
Int_t K_plus_GD_MOTHER_ID;
Int_t K_plus_GD_GD_MOTHER_ID;
Int_t K_short_GD_MOTHER_ID;
Int_t Ks_pi_plus_GD_MOTHER_ID;
Int_t Ks_pi_minus_GD_MOTHER_ID;
Int_t pi_plus_GD_MOTHER_ID;
Int_t pi_zero_GD_MOTHER_ID;
Int_t mu_plus_GD_MOTHER_ID;
Int_t mu_minus_GD_MOTHER_ID;
Int_t gamma1_GD_MOTHER_ID;
Int_t gamma2_GD_MOTHER_ID;
//assign variables to TRUE_GD_MOTHER_ID
Int_t Ks_pi_plus_GD_GD_MOTHER_ID;
Int_t Ks_pi_minus_GD_GD_MOTHER_ID;
Int_t gamma1_GD_GD_MOTHER_ID;
Int_t gamma2_GD_GD_MOTHER_ID;
Int_t pi_zero_merged_MC_GD_GD_MOTHER_ID;
//link variables to branches
tree -> SetBranchAddress( "B_plus_M_DTF" , &B_plus_M_DTF );
tree -> SetBranchAddress( "B_plus_TRUEID" , &B_plus_TRUEID );
tree -> SetBranchAddress( "K_star_plus_TRUEID" , &K_star_plus_TRUEID );
tree -> SetBranchAddress( "mu_plus_TRUEID" , &mu_plus_TRUEID );
tree -> SetBranchAddress( "mu_minus_TRUEID" , &mu_minus_TRUEID );
tree -> SetBranchAddress( "B_plus_MC_MOTHER_ID" , &B_plus_MOTHER_ID );
tree -> SetBranchAddress( "K_star_plus_MC_MOTHER_ID" , &K_star_plus_MOTHER_ID );
tree -> SetBranchAddress( "mu_plus_MC_MOTHER_ID" , &mu_plus_MOTHER_ID );
tree -> SetBranchAddress( "mu_minus_MC_MOTHER_ID" , &mu_minus_MOTHER_ID );
tree -> SetBranchAddress( "mu_plus_MC_GD_MOTHER_ID" , &mu_plus_GD_MOTHER_ID );
tree -> SetBranchAddress( "mu_minus_MC_GD_MOTHER_ID" , &mu_minus_GD_MOTHER_ID );
if(Kst2Kpluspi0Resolved){
tree -> SetBranchAddress( "K_plus_TRUEID" , &K_plus_TRUEID );
tree -> SetBranchAddress( "pi_zero_resolved_TRUEID" , &pi_zero_TRUEID );
tree -> SetBranchAddress( "gamma1_TRUEID" , &gamma1_TRUEID );
tree -> SetBranchAddress( "gamma2_TRUEID" , &gamma2_TRUEID );
tree -> SetBranchAddress( "K_plus_MC_MOTHER_ID" , &K_plus_MOTHER_ID );
tree -> SetBranchAddress( "pi_zero_resolved_MC_MOTHER_ID" , &pi_zero_MOTHER_ID );
tree -> SetBranchAddress( "gamma1_MC_MOTHER_ID" , &gamma1_MOTHER_ID );
tree -> SetBranchAddress( "gamma2_MC_MOTHER_ID" , &gamma2_MOTHER_ID );
tree -> SetBranchAddress( "K_plus_MC_GD_MOTHER_ID" , &K_plus_GD_MOTHER_ID );
tree -> SetBranchAddress( "K_plus_MC_GD_GD_MOTHER_ID" , &K_plus_GD_GD_MOTHER_ID );
tree -> SetBranchAddress( "pi_zero_resolved_MC_GD_MOTHER_ID" , &pi_zero_GD_MOTHER_ID );
tree -> SetBranchAddress( "gamma1_MC_GD_MOTHER_ID" , &gamma1_GD_MOTHER_ID );
tree -> SetBranchAddress( "gamma2_MC_GD_MOTHER_ID" , &gamma2_GD_MOTHER_ID );
tree -> SetBranchAddress( "gamma1_MC_GD_GD_MOTHER_ID" , &gamma1_GD_GD_MOTHER_ID );
tree -> SetBranchAddress( "gamma2_MC_GD_GD_MOTHER_ID" , &gamma2_GD_GD_MOTHER_ID );
}
if(Kst2Kpluspi0Merged){
tree -> SetBranchAddress( "K_plus_TRUEID" , &K_plus_TRUEID );
tree -> SetBranchAddress( "pi_zero_merged_TRUEID" , &pi_zero_TRUEID );
tree -> SetBranchAddress( "K_plus_MC_MOTHER_ID" , &K_plus_MOTHER_ID );
tree -> SetBranchAddress( "pi_zero_merged_MC_MOTHER_ID" , &pi_zero_MOTHER_ID );
tree -> SetBranchAddress( "K_plus_MC_GD_MOTHER_ID" , &K_plus_GD_MOTHER_ID );
tree -> SetBranchAddress( "pi_zero_merged_MC_GD_MOTHER_ID" , &pi_zero_GD_MOTHER_ID );
tree -> SetBranchAddress( "pi_zero_merged_MC_GD_GD_MOTHER_ID" , &pi_zero_merged_MC_GD_GD_MOTHER_ID );
}
if(Kst2Kspiplus){
tree -> SetBranchAddress( "K_short_TRUEID" , &K_short_TRUEID );
tree -> SetBranchAddress( "pi_plus_TRUEID" , &pi_plus_TRUEID );
tree -> SetBranchAddress( "Ks_pi_plus_TRUEID" , &Ks_pi_plus_TRUEID );
tree -> SetBranchAddress( "Ks_pi_minus_TRUEID" , &Ks_pi_minus_TRUEID );
tree -> SetBranchAddress( "K_short_MC_MOTHER_ID" , &K_short_MOTHER_ID );
tree -> SetBranchAddress( "pi_plus_MC_MOTHER_ID" , &pi_plus_MOTHER_ID );
tree -> SetBranchAddress( "Ks_pi_plus_MC_MOTHER_ID" , &Ks_pi_plus_MOTHER_ID );
tree -> SetBranchAddress( "Ks_pi_minus_MC_MOTHER_ID" , &Ks_pi_minus_MOTHER_ID );
tree -> SetBranchAddress( "K_short_MC_GD_MOTHER_ID" , &K_short_GD_MOTHER_ID );
tree -> SetBranchAddress( "pi_plus_MC_GD_MOTHER_ID" , &pi_plus_GD_MOTHER_ID );
tree -> SetBranchAddress( "Ks_pi_plus_MC_GD_MOTHER_ID" , &Ks_pi_plus_GD_MOTHER_ID );
tree -> SetBranchAddress( "Ks_pi_minus_MC_GD_MOTHER_ID" , &Ks_pi_minus_GD_MOTHER_ID );
tree -> SetBranchAddress( "Ks_pi_plus_MC_GD_GD_MOTHER_ID" , &Ks_pi_plus_GD_GD_MOTHER_ID );
tree -> SetBranchAddress( "Ks_pi_minus_MC_GD_GD_MOTHER_ID" , &Ks_pi_minus_GD_GD_MOTHER_ID );
}
tree -> SetBranchAddress( "J_psi_TRUEID" , &J_psi_TRUEID );
tree -> SetBranchAddress( "J_psi_MC_MOTHER_ID" , &J_psi_MOTHER_ID );
//counter
Int_t B_plus_VETO = 0;
Int_t K_star_plus_VETO = 0;
Int_t K_plus_VETO = 0;
Int_t K_short_VETO = 0;
Int_t Ks_pi_plus_VETO = 0;
Int_t Ks_pi_minus_VETO = 0;
Int_t pi_plus_VETO = 0;
Int_t pi_zero_VETO = 0;
Int_t mu_plus_VETO = 0;
Int_t mu_minus_VETO = 0;
Int_t J_psi_VETO = 0;
Int_t gamma1_VETO = 0;
Int_t gamma2_VETO = 0;
Int_t MOTHER_ID_VETO = 0;
Int_t GD_MOTHER_ID_VETO = 0;
Int_t GD_GD_MOTHER_ID_VETO = 0;
Int_t TRUE_events = 0;
//count the numbers of Kshort daughters, which are pi_zeros
Int_t Kshort_pi_zero_daughters = 0;
//count muons, mis-IDed as pions
Int_t ReconstructedPionsWhichAreMuons = 0;
//======================================
//Get histos of all possible backgrounds
//======================================
//control plots for pi0
Double_t pi_zero_resolved_M = 0.0;
if (Kst2Kpluspi0Resolved)
tree -> SetBranchAddress( "pi_zero_resolved_M" , &pi_zero_resolved_M );
TH1 * h_pi0_resolved_mass= new TH1D("h_pi0_resolved_mass" , "m_{#pi_{0}} resolved", 100, 75, 175);
h_pi0_resolved_mass->GetXaxis()->SetTitle("m_{pi_{0}} resolved [MeV]");
h_pi0_resolved_mass->GetYaxis()->SetTitle("Counts / 1MeV");
TH1 * h_pi0_mismatched_mass= new TH1D("h_pi0_mismatched_mass" , "m_{#pi_{0}} IDed as gamma", 100, 75, 175);
h_pi0_mismatched_mass->GetXaxis()->SetTitle("m_{pi_{0}}[MeV] MisIDed");
h_pi0_mismatched_mass->GetYaxis()->SetTitle("Counts / 1MeV");
TH1 * h_B_plus_TRUEID = new TH1D("h_B_plus_TRUEID", "B^{+} TRUEID", 100000, -1, 99999);
TH1 * h_K_star_plus_TRUEID = new TH1D("h_K_star_plus_TRUEID", "K^{*+} TRUEID", 100000, -1, 99999);
TH1 * h_K_plus_TRUEID = new TH1D("h_K_plus_TRUEID", "K^{+} TRUEID", 100000, -1, 99999);
TH1 * h_pi_zero_TRUEID = new TH1D("h_pi_zero_TRUEID", "#pi^{0} TRUEID", 100000, -1, 99999);
TH1 * h_B_plus_TRUEID_rm = new TH1D("h_B_plus_TRUEID_rm", "B^{+} TRUEID", 100000, -1, 99999);
TH1 * h_K_star_plus_TRUEID_rm= new TH1D("h_K_star_plus_TRUEID_rm", "K^{*+} TRUEID", 100000, -1, 99999);
TH1 * h_K_plus_TRUEID_rm = new TH1D("h_K_plus_TRUEID_rm", "K^{+} TRUEID", 100000, -1, 99999);
TH1 * h_pi_zero_TRUEID_rm = new TH1D("h_pi_zero_TRUEID_rm", "#pi^{0} TRUEID", 100000, -1, 99999);
TH1 * h_B_plus_K1 = new TH1D("h_B_plus_K1", "h_B_plus_K1", 100, cut_B_plus_M_low, cut_B_plus_M_high);
TH1 * h_B_plus_Kst = new TH1D("h_B_plus_Kst", "h_B_plus_Kst", 100, cut_B_plus_M_low, cut_B_plus_M_high);
TH1 * h_B_plus_Kpl = new TH1D("h_B_plus_Kpl", "h_B_plus_Kpl", 100, cut_B_plus_M_low, cut_B_plus_M_high);
vector<std::string> Cuts;
string tmp_string = "";
int nBins = 75;
TH1 * h_BplusMass_tot= new TH1D("h_BplusMass_tot" , "h_BplusMass_tot", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_tot->GetXaxis()->SetTitle("B^{+} mass [MeV]");
//======================================
// Define trees for all channels
//======================================
TTree *tree_correct = tree->CloneTree(0);
tree_correct->SetName("tree_correct");
TH1 * h_BplusMass_correct = new TH1D("h_BplusMass_correct" , "h_BplusMass_correct", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_correct->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_K1_1270 = tree->CloneTree(0);
tree_K1_1270->SetName("tree_K1_1270");
TH1 * h_BplusMass_K1_1270 = new TH1D("h_BplusMass_K1_1270" , "h_BplusMass_K1_1270", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_K1_1270->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_K1_1410 = tree->CloneTree(0);
tree_K1_1410->SetName("tree_K1_1410");
TH1 * h_BplusMass_K1_1410 = new TH1D("h_BplusMass_K1_1410" , "h_BplusMass_K1_1410", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_K1_1410->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_K1_1400 = tree->CloneTree(0);
tree_K1_1400->SetName("tree_K1_1400");
TH1 * h_BplusMass_K1_1400 = new TH1D("h_BplusMass_K1_1400" , "h_BplusMass_K1_1400", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_K1_1400->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_K_plus= tree->CloneTree(0);
tree_K_plus->SetName("tree_K_plus");
TH1 * h_BplusMass_K_plus = new TH1D("h_BplusMass_K_plus" , "h_BplusMass_K_plus", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_K_plus->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_K_plus_misID= tree->CloneTree(0);
tree_K_plus_misID->SetName("tree_K_plus_misID");
TH1 * h_BplusMass_K_plus_misID = new TH1D("h_BplusMass_K_plus_misID" , "h_BplusMass_K_plus_misID", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_K_plus_misID->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_B_star= tree->CloneTree(0);
tree_B_star->SetName("tree_B_star");
TH1 * h_BplusMass_B_star = new TH1D("h_BplusMass_B_star" , "h_BplusMass_B_star", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_B_star->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_Kstar_misID= tree->CloneTree(0);
tree_Kstar_misID->SetName("tree_Kstar_misID");
TH1 * h_BplusMass_Kstar_misID = new TH1D("h_BplusMass_Kstar_misID" , "h_BplusMass_Kstar_misID", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_Kstar_misID->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_bquark= tree->CloneTree(0);
tree_bquark->SetName("tree_bquark");
TH1 * h_BplusMass_bquark = new TH1D("h_BplusMass_bquark" , "h_BplusMass_bquark", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_bquark->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_combinatorics= tree->CloneTree(0);
tree_combinatorics->SetName("tree_combinatorics");
TH1 * h_BplusMass_combinatorics = new TH1D("h_BplusMass_combinatorics" , "h_BplusMass_combinatorics", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_combinatorics->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TTree *tree_others = tree->CloneTree(0);
tree_others->SetName("tree_others");
TH1 * h_BplusMass_others = new TH1D("h_BplusMass_others" , "h_BplusMass_others", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_others->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TH1 * h_BplusMass_JpsiVETO = new TH1D("h_BplusMass_JpsiVETO" , "h_BplusMass_JpsiVETO", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_JpsiVETO->GetXaxis()->SetTitle("B^{+} mass [MeV]");
TH1 * h_BplusMass_GDVETO = new TH1D("h_BplusMass_GDVETO" , "h_BplusMass_GDVETO", nBins ,cut_B_plus_M_low, cut_B_plus_M_high);
h_BplusMass_GDVETO->GetXaxis()->SetTitle("B^{+} mass [MeV]");
//======================================
//check all particles for correct TRUEID
//======================================
//MC loop
cout << "Tree has " << tree->GetEntries() << " entries" << endl;
int nEvents = tree->GetEntries();
std::cout << "Truth matching of " << nEvents << " INC MC events for " << year << magnet << "!" << std::endl;
for(int i = 0; i < nEvents; i++){
if(i%2000==0 && i != 0)std::cout << "Loading MC event " << i << " / " << nEvents << std::endl;
bool flag_other = true;
//others
tree->GetEntry(i);
h_BplusMass_tot->Fill(B_plus_M_DTF);
//remove nonsense from B and Jpsi
if(TMath::Abs(B_plus_TRUEID) != TRUEID.B_PLUS && TMath::Abs(B_plus_TRUEID) != 0 && TMath::Abs(B_plus_TRUEID) != 5){
B_plus_VETO++;
h_B_plus_TRUEID_rm->Fill(B_plus_TRUEID);
continue;
}
if(ReferenceChannel){
if(TMath::Abs(J_psi_TRUEID) != TRUEID.J_PSI){
J_psi_VETO++;
h_BplusMass_JpsiVETO->Fill(B_plus_M_DTF);
continue;
}
if(TMath::Abs(mu_plus_MOTHER_ID) != TRUEID.J_PSI){
MOTHER_ID_VETO++;
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
if(TMath::Abs(mu_minus_MOTHER_ID) != TRUEID.J_PSI){
MOTHER_ID_VETO++;
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
if(TMath::Abs(mu_plus_GD_MOTHER_ID) != TRUEID.B_PLUS){
GD_MOTHER_ID_VETO++;
h_BplusMass_GDVETO->Fill(B_plus_M_DTF);
continue;
}
if(TMath::Abs(mu_minus_GD_MOTHER_ID) != TRUEID.B_PLUS){
GD_MOTHER_ID_VETO++;
h_BplusMass_GDVETO->Fill(B_plus_M_DTF);
continue;
}
}
else{
if(TMath::Abs(J_psi_TRUEID) == TRUEID.J_PSI){
J_psi_VETO++;
h_BplusMass_JpsiVETO->Fill(B_plus_M_DTF);
continue;
}
if(TMath::Abs(mu_plus_MOTHER_ID) != TRUEID.B_PLUS){
MOTHER_ID_VETO++;
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
if(TMath::Abs(mu_minus_MOTHER_ID) != TRUEID.B_PLUS){
MOTHER_ID_VETO++;
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
}
//correct tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_STAR_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.B_PLUS)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.B_PLUS){//Grandmother IDs
tree_correct->Fill();
h_BplusMass_correct->Fill(B_plus_M_DTF);
continue;
}
}
}
//K1 tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_ONE_PLUS || TMath::Abs(K_star_plus_TRUEID) == 0 || TMath::Abs(K_star_plus_TRUEID) == TRUEID.B_PLUS )
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_STAR_PLUS || TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS || TMath::Abs(K_star_plus_MOTHER_ID) == 0)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS || TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.K_ONE_PLUS){//Grandmother IDs
tree_K1_1270->Fill();
h_BplusMass_K1_1270->Fill(B_plus_M_DTF);
continue;
}
}
}
//K1 1400 tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_ONE_PLUS_1400 || TMath::Abs(K_star_plus_TRUEID) == 0|| TMath::Abs(K_star_plus_TRUEID) == TRUEID.B_PLUS )
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_STAR_PLUS || TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle ID
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS_1400)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS || TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS_1400){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.K_ONE_PLUS_1400){//Grandmother IDs
tree_K1_1400->Fill();
h_BplusMass_K1_1410->Fill(B_plus_M_DTF);
tree_others->Fill();
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
}
}
//K1 1410 tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_ONE_PLUS_1410 || TMath::Abs(K_star_plus_TRUEID) == 0 )
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_star_plus_MOTHER_ID) == 0 || TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS_1410)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS || TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS_1410){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == 0 || TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.K_ONE_PLUS_1410){//Grandmother IDs
tree_K1_1410->Fill();
h_BplusMass_K1_1410->Fill(B_plus_M_DTF);
tree_others->Fill();
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
}
}
//K1 1270 tree 2.0
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_STAR_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.K_ONE_PLUS){//Grandmother IDs
tree_K1_1270->Fill();
h_BplusMass_K1_1270->Fill(B_plus_M_DTF);
continue;
}
}
}
//tree_K_plus_misID tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_ONE_PLUS || TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_STAR_PLUS|| TMath::Abs(K_star_plus_TRUEID) == 0 )
if(TMath::Abs(K_plus_TRUEID) == TRUEID.PI_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS || TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_ONE_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.K_ONE_PLUS){//Grandmother IDs
tree_K_plus_misID->Fill();
h_BplusMass_K_plus_misID->Fill(B_plus_M_DTF);
tree_others->Fill();
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
}
}
//tree_K_plus_misID tree 2.0 (K+->pi+)
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_STAR_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.PI_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.B_PLUS)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.K_STAR_PLUS)
if (TMath::Abs(K_plus_GD_GD_MOTHER_ID) == TRUEID.B_PLUS){//Grandmother IDs
tree_K_plus_misID->Fill();
h_BplusMass_K_plus_misID->Fill(B_plus_M_DTF);
tree_others->Fill();
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
}
}
//Kplus tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_PLUS || TMath::Abs(K_star_plus_TRUEID) == 0)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.B_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == 0){//Grandmother IDs
tree_K_plus->Fill();
h_BplusMass_K_plus->Fill(B_plus_M_DTF);
tree_others->Fill();
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
}
}
//K star misID tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.B_PLUS){//Grandmother IDs
tree_Kstar_misID->Fill();
h_BplusMass_Kstar_misID->Fill(B_plus_M_DTF);
tree_others->Fill();
h_BplusMass_others->Fill(B_plus_M_DTF);
continue;
}
}
}
//B+ <-> b quark
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.B_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == 5)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.B_PLUS){ //Mother IDs
tree_bquark->Fill();
h_BplusMass_bquark->Fill(B_plus_M_DTF);
continue;
}
}
//B+ <-> b quark 2.0
if (TMath::Abs(K_star_plus_TRUEID) ==5)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == 0)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.B_PLUS || TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS ){ //Mother IDs
tree_bquark->Fill();
h_BplusMass_bquark->Fill(B_plus_M_DTF);
continue;
}
}
//Bstar tree
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.B_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == 523)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.B_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == 523){//Grandmother IDs
tree_B_star->Fill();
h_BplusMass_B_star->Fill(B_plus_M_DTF);
continue;
}
}
}
//Bstar tree 2.0
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == 523)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.B_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == 523){//Grandmother IDs
tree_B_star->Fill();
h_BplusMass_B_star->Fill(B_plus_M_DTF);
continue;
}
}
}
//Bstar tree 3.0 (K* misID B+)
if (TMath::Abs(K_star_plus_TRUEID) == TRUEID.K_PLUS)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == TRUEID.B_PLUS)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.B_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == 523){//Grandmother IDs
tree_B_star->Fill();
h_BplusMass_B_star->Fill(B_plus_M_DTF);
continue;
}
}
}
//Bstar tree 4.0 (K* misID B+)
if (TMath::Abs(K_star_plus_TRUEID) == 0)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == 0)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.B_PLUS){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == 523 ||TMath::Abs(K_plus_GD_MOTHER_ID) == 10521){//Grandmother IDs
tree_B_star->Fill();
h_BplusMass_B_star->Fill(B_plus_M_DTF);
continue;
}
}
}
//combinatorics tree
if (TMath::Abs(K_star_plus_TRUEID) == 0)
if(TMath::Abs(K_plus_TRUEID) == TRUEID.K_PLUS ||TMath::Abs(K_plus_TRUEID) == 0 ){ //Paticle IDs
if(TMath::Abs(K_star_plus_MOTHER_ID) == 0)
if(TMath::Abs(K_plus_MOTHER_ID) == TRUEID.K_STAR_PLUS ||TMath::Abs(K_plus_TRUEID) == 0 ){ //Mother IDs
if(TMath::Abs(K_plus_GD_MOTHER_ID) == TRUEID.B_PLUS){//Grandmother IDs
tree_combinatorics->Fill();
h_BplusMass_combinatorics->Fill(B_plus_M_DTF);
continue;
}
}
}
//the rest of stuff
tree_others->Fill();
h_BplusMass_others->Fill(B_plus_M_DTF);
}
//output truthmatching results
std::cout << "Resulting TruthMatching data: " << TRUE_events << "/" << nEvents << " are found to be correctly reconstructed!" << std::endl;
std::cout << "The incorrect reconstructed particles are the following:" << std::endl;
std::cout << "(only the first found-incorrect particle is counted)" << std::endl;
std::cout << std::endl;
//std::cout << (B0 ? "B_zero:\t\t" : "B_plus:\t\t") << B_plus_VETO << std::endl;
//std::cout << (B0 ? "K_star_zero:\t" : ( K1 ? "K_one_plus:\t" : "K_star_plus:\t")) << K_star_plus_VETO << std::endl;
std::cout << "mu_plus:\t" << mu_plus_VETO << std::endl;
std::cout << "mu_minus:\t" << mu_minus_VETO << std::endl;
if(Kst2Kpluspi0Resolved){
std::cout << "K_plus:\t\t" << K_plus_VETO << std::endl;
std::cout << "pi_zero:\t" << pi_zero_VETO << std::endl;
std::cout << "gamma1:\t\t" << gamma1_VETO << std::endl;
std::cout << "gamma2:\t\t" << gamma2_VETO << std::endl;
}
if(Kst2Kpluspi0Merged){
std::cout << "K_plus:\t\t" << K_plus_VETO << std::endl;
std::cout << "pi_zero:\t" << pi_zero_VETO << std::endl;
}
if(Kst2Kspiplus){
std::cout << "K_short:\t" << K_short_VETO << std::endl;
std::cout << "Ks_pi_plus:\t" << Ks_pi_plus_VETO << std::endl;
std::cout << "Ks_pi_minus:\t" << Ks_pi_minus_VETO << std::endl;
std::cout << "pi_plus:\t" << pi_plus_VETO << std::endl;
}
std::cout << "J_psi:\t\t" <<J_psi_VETO << std::endl;
if(Kst2Kspiplus){
std::cout << "Pi0 daughers of the Kshort:\t" << Kshort_pi_zero_daughters << std::endl;
std::cout << "Muons mis-IDed as pions:\t" << ReconstructedPionsWhichAreMuons << "/" << Ks_pi_plus_VETO + Ks_pi_minus_VETO + pi_plus_VETO << std::endl;
}
std::cout << "Total incorrect family members of all particles: " << MOTHER_ID_VETO + GD_MOTHER_ID_VETO + GD_GD_MOTHER_ID_VETO << std::endl;
std::cout << "MotherID:\t" << MOTHER_ID_VETO << std::endl;
std::cout << "GDmotherID:\t" << GD_MOTHER_ID_VETO << std::endl;
std::cout << "GDGDmotherID:\t" << GD_GD_MOTHER_ID_VETO << std::endl;
//save to file
TFile * outputFile = nullptr;
// Kst2Kpluspi0Resolved case
if(Kst2Kpluspi0Resolved){
outputFile = new TFile(Form("%s/data/MC/BtoXJpsi/%s%s/%s%s_pi0Resolved_plots.root",path_to_output_KplusPizero.c_str(),year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()), "RECREATE");
outputFile->cd();
tree_correct->Write();
cout << "Correctly id'd entries " << tree_correct->GetEntries() << endl;
tree_K1_1270->Write();
cout << "B+->K1(1270) entries " << tree_K1_1270->GetEntries() << endl;
tree_K1_1400->Write();
cout << "B+->K1(1400) entries " << tree_K1_1400->GetEntries() << endl;
tree_K1_1410->Write();
cout << "B+->K1(1410) entries " << tree_K1_1410->GetEntries() << endl;
tree_K_plus->Write();
cout << "B+->K+pi0 entries " << tree_K_plus->GetEntries() << endl;
tree_K_plus_misID->Write();
cout << "B+->K*->pi+pi0 " << tree_K_plus_misID->GetEntries() << endl;
tree_B_star->Write();
cout << "B*-> B+->K+ entries " << tree_B_star->GetEntries() << endl;
tree_Kstar_misID->Write();
cout << "K* misID entries " << tree_Kstar_misID->GetEntries() << endl;
tree_combinatorics->Write();
cout << "b-quark entries " << tree_bquark->GetEntries() << endl;
tree_bquark->Write();
cout << "Combinatorics " << tree_combinatorics->GetEntries() << endl;
tree_others->Write();
cout << "Other entries " << tree_others->GetEntries() << endl;
THStack *tmp = new THStack("IncMC","IncMC;B^{+} mass [MeV];Entries");
//tmp->GetXaxis()->SetTitle("B^{+} mass [MeV]");
//tmp->GetYaxis()->SetTitle("Entries");
h_BplusMass_others->SetFillColor(kGreen);
h_BplusMass_others->SetLineColor(kGreen);
tmp->Add(h_BplusMass_others);
h_BplusMass_bquark->SetFillColor(kGreen);
h_BplusMass_bquark->SetLineColor(kGreen);
tmp->Add(h_BplusMass_bquark);
h_BplusMass_K1_1270->SetFillColor(kBlue);
h_BplusMass_K1_1270->SetLineColor(kBlue);
tmp->Add(h_BplusMass_K1_1270);
// h_BplusMass_K1_1400->SetFillColor(0);
// h_BplusMass_K1_1400->SetLineColor(0);
// tmp->Add(h_BplusMass_K1_1400);
// h_BplusMass_K1_1410->SetFillColor(0);
// h_BplusMass_K1_1410->SetLineColor(0);
// tmp->Add(h_BplusMass_K1_1410);
// h_BplusMass_K_plus->SetFillColor(0);
// h_BplusMass_K_plus->SetLineColor(0);
// tmp->Add(h_BplusMass_K_plus);
// h_BplusMass_K_plus_misID->SetFillColor(0);
// h_BplusMass_K_plus_misID->SetLineColor(0);
// tmp->Add(h_BplusMass_K_plus_misID);
h_BplusMass_B_star->SetFillColor(kYellow);
h_BplusMass_B_star->SetLineColor(kYellow);
tmp->Add(h_BplusMass_B_star);
// h_BplusMass_Kstar_misID->SetFillColor(0);
// h_BplusMass_Kstar_misID->SetLineColor(0);
// tmp->Add(h_BplusMass_Kstar_misID);
//
// h_BplusMass_combinatorics->SetFillColor(kBlack);
// h_BplusMass_combinatorics->SetLineColor(kBlack);
// tmp->Add(h_BplusMass_combinatorics);
h_BplusMass_JpsiVETO->SetFillColor(kMagenta-10);
h_BplusMass_JpsiVETO->SetLineColor(kMagenta-10);
tmp->Add(h_BplusMass_JpsiVETO);
h_BplusMass_GDVETO->SetFillColor(kOrange-3);
h_BplusMass_GDVETO->SetLineColor(kOrange-3);
tmp->Add(h_BplusMass_GDVETO);
h_BplusMass_correct->SetFillColor(kRed);
h_BplusMass_correct->SetLineColor(kRed);
tmp->Add(h_BplusMass_correct);
tmp->Write();
TCanvas *c1 = new TCanvas("c1","all",900,700);
c1->cd();
gStyle->SetOptStat(0);
// h_BplusMass_correct->SetFillColor(0);
h_BplusMass_correct->Draw();
// h_BplusMass_others->SetFillColor(0);
h_BplusMass_others->Draw("SAME");
// h_BplusMass_K1_1270->SetFillColor(0);
h_BplusMass_K1_1270->Draw("SAME");
//h_BplusMass_K1_1400->Draw("SAME");
//h_BplusMass_K1_1410->Draw("SAME");
//h_BplusMass_K_plus->Draw("SAME");
//h_BplusMass_K_plus_misID->Draw("SAME");
// h_BplusMass_B_star->SetFillColor(0);
h_BplusMass_B_star->Draw("SAME");
//h_BplusMass_Kstar_misID->Draw("SAME");
//h_BplusMass_combinatorics->Draw("SAME");
// h_BplusMass_bquark->SetFillColor(0);
h_BplusMass_bquark->Draw("SAME");
// h_BplusMass_JpsiVETO->SetFillColor(0);
h_BplusMass_JpsiVETO->Draw("SAME");
// h_BplusMass_GDVETO->SetFillColor(0);
h_BplusMass_GDVETO->Draw("SAME");
c1->Write();
TCanvas *c2 = new TCanvas("c2","al2",900,700);
c2->cd();
tmp->Draw();
h_BplusMass_tot->SetLineWidth(3);
h_BplusMass_tot->SetLineColor(kBlack);
h_BplusMass_tot->Draw("SAMEE");
c2->Write();
c2->SaveAs(Form("%s/data/MC/BtoXJpsi/%s%s/%s%s_pi0Resolved%s_BKG.eps",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str(),(smallSample ? "_small":"") ),"eps");
}
// Kst2Kpluspi0Merged case
else if(Kst2Kpluspi0Merged){
//TODO
}
// Kst2Kspiplus case
else if(Kst2Kspiplus){
//TODO
}
else return 0;
if(!outputFile->IsOpen()){
std::cout << "[ERROR]\t\tFile was not opened succesfully!" << std::endl;
return 0;
}
outputFile->cd();
//treeTruthMatched->Write("",TObject::kWriteDelete);
if (Kst2Kpluspi0Resolved){
h_pi0_mismatched_mass->SetLineColor(kRed);
h_pi0_mismatched_mass->SetLineWidth(3);
h_pi0_resolved_mass->SetLineColor(kBlue);
h_pi0_resolved_mass->SetLineWidth(3);
h_pi0_mismatched_mass->Write("",TObject::kWriteDelete);
h_pi0_resolved_mass->Write("",TObject::kWriteDelete);
h_B_plus_TRUEID ->Write("",TObject::kWriteDelete);
h_K_star_plus_TRUEID ->Write("",TObject::kWriteDelete);
h_K_plus_TRUEID ->Write("",TObject::kWriteDelete);
h_pi_zero_TRUEID ->Write("",TObject::kWriteDelete);
h_B_plus_TRUEID_rm ->Write("",TObject::kWriteDelete);
h_K_star_plus_TRUEID_rm ->Write("",TObject::kWriteDelete);
h_K_plus_TRUEID_rm ->Write("",TObject::kWriteDelete);
h_pi_zero_TRUEID_rm ->Write("",TObject::kWriteDelete);
h_B_plus_K1 ->Write("",TObject::kWriteDelete);
h_B_plus_Kst ->Write("",TObject::kWriteDelete);
h_B_plus_Kpl ->Write("",TObject::kWriteDelete);
}
delete outputFile;//->Close();
delete h_pi0_mismatched_mass;
delete h_pi0_resolved_mass;
delete tree;
std::cout << "[DONE]\t\tFinished Truth-Matching for INC sample: " << year << " " << magnet << std::endl;
return 1;
}
int plotIncSampleAll() {
if ( plotIncSample("2011","down") == 0) return 0;
if ( plotIncSample("2011","up") == 0) return 0;
if ( plotIncSample("2012","down") == 0) return 0;
if ( plotIncSample("2012","up") == 0) return 0;
if ( plotIncSample("2015","down") == 0) return 0;
if ( plotIncSample("2015","up") == 0) return 0;
if ( plotIncSample("2016","down") == 0) return 0;
if ( plotIncSample("2016","up") == 0) return 0;
return 1;
}