kopecna/EWP-BplusToKstMuMu-AngAna
kopecna/EWP-BplusToKstMuMu-AngAna
2
1
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity
0f8ca4f87d
EWP-BplusToKstMuMu-AngAna/Code/Selection/BDTSelection.cpp

2650 lines
128 KiB
C++
Executable File
Raw Blame History

//offline selection of BDT variables for B+->Kst+mumu
//david gerick
#include "GlobalFunctions.hh"
#include "Utils.hpp"
#include "Paths.hpp"
using namespace std;
using namespace RooFit ;
//////////////////////////////////////////////////////////////////////
/// preSelection()
/// loops over all .root-files created by subjobs in ganga and converts into one tree. trigger, kinematic, PID and
/// some more cuts are applied. Q^2-Binning applied, Kshort becomes extra branch for decay inside or outside of
/// VELO. histograms of inv. mass distributions are created and saved in separate .root-file
///
/// runAllYearsAndMagPol()
/// preSelection is executed for 2011 + 2012 as well as MagUp and MagDown of the decay
///
/// runAllSignalData()
/// all data for this decay is evaluated
///
/// runAllSignalMC()
/// runs all MC files for the signal channel
///
/// runAllRefMC()
/// same as above but for reference channel
///
/// runComplete()
/// every available data for all three (sub)decays is evaluated, data and MC for signal and reference channel
///
/// getTCKs()
/// return the exact events numbers for every MagDown and MagUp, Year and TCK configuration
///
/// plotHistos()
/// Selected histograms generated in the preSelection() function are plotted and nicely formatted before
/// saved as .eps file
int copyDataForTMVA(std::string year = "2011", std::string magnet = "down") {
string string_input = GetInputFile(year,magnet,true,false,false,false,smallSample);
TChain* tree = new TChain("DecayTree");
tree->Add(string_input.c_str());
string output_path = string_input;
replace(output_path,".root","_tmp.root");
if(tree->GetEntries() == 0){
//If it complains while running, it's just trying to read also from David's preselected files
coutERROR("No entries found in TTree from tfile " + string_input + ". Exit program!");
coutERROR("Exit program!");
return 0;
}
TFile *output = new TFile(output_path.c_str(),"RECREATE");
tree->SetBranchStatus("*",1);
TTree * NewTree = tree->CloneTree(0);
bool TMedBKGCAT = true;
bool TMed_noPi0 = true;
bool TMed = true;
bool TM_gammas = true;
NewTree->Branch("TMedBKGCAT", &TMedBKGCAT, "TMedBKGCAT/b");
NewTree->Branch("TMed_noPi0", &TMed_noPi0, "TMed_noPi0/b");
NewTree->Branch("TMed", &TMed, "TMed/b");
NewTree->Branch("TM_gammas", &TM_gammas, "TM_gammas/b");
for (int evt = 0; evt < tree->GetEntries(); evt++){
tree->GetEntry(evt);
TMedBKGCAT = true;
TMed_noPi0 = true;
TMed = true;
TM_gammas = true;
NewTree->Fill();
}
output->cd();
NewTree->Write("",TObject::kWriteDelete);
output->Close();
cout << "[INFO]\t\t New tree created." << endl;
//careful, possibly leaky (wouldnẗ life be so much easier if one could just change ROOT tree on the spot...
string command = "mv " + string_input + " ";
replace(string_input,".root","_backup.root");
command += string_input;
system(command.c_str());
replace(string_input,"_backup.root",".root");
command = "mv " + output_path + " " + string_input;
system(command.c_str());
return 1;
}
int copyDataForTMVAAll(int Run = 1){
for (auto& year : yearsData(Run)){
if (copyDataForTMVA(year,"down") ==0) return 0;
if (copyDataForTMVA(year,"up") ==0) return 0;
}
return 1;
}
//the function to loop over all events and merge them into one file; apply cut-based selection while merging
int preSelection(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false, bool B0 = false, bool K1 = false, bool Inc = false) {
TChain* tree = 0;
//Reference channel and PHSP are both MC data
if(ReferenceChannel)MC = true;
if(PHSP)MC = true;
if(K1)MC = true;
if(Inc)MC = true;
if(ReferenceChannel && PHSP){
coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
PHSP = false;
}
//Kst2Kpluspi0Resolved case
if(Kst2Kpluspi0Resolved){
tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
if (B0 || K1 || Inc){
tree->Add(GetInputFileBkg(year,magnet,false,ReferenceChannel,B0, K1, Inc).c_str());
coutDebug("Adding " + GetInputFileBkg(year,magnet,false,ReferenceChannel,B0, K1, Inc));
}
else{
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
coutDebug("Adding " + GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample));
}
}
//Kst2Kpluspi0Merged case
if(Kst2Kpluspi0Merged){
tree=new TChain("b2KstKpi0mumuMergedTuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
cout << "[DEBUG]\t\t Openning " << GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str() << endl;
}
//Kst2Kspiplus case
if(Kst2Kspiplus){
tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
}
int N = tree->GetEntries();
cout << "Old files contain " << N << " events" << endl;
if(N == 0){
coutERROR("Tree has no entries. Exit!");
return 0;
}
//Disable all branches
tree->SetBranchStatus("*",0);
//activate needed branches
tree->SetBranchStatus("B_plus_*Decision*",1) ;
tree->SetBranchStatus("*PID*",1) ;
tree->SetBranchStatus("*ProbNN*",1) ;
//tree->SetBranchStatus("*GhostProb",1) ;
tree->SetBranchStatus("nCandidate",1) ;
tree->SetBranchStatus("*Tracks",1) ;
tree->SetBranchStatus("nVeloClusters", 1) ;
tree->SetBranchStatus("nTTClusters", 1) ;
tree->SetBranchStatus("nITClusters", 1) ;
tree->SetBranchStatus("nOTClusters", 1) ;
tree->SetBranchStatus("nSPDHits", 1) ;
tree->SetBranchStatus("nPV",1) ;
tree->SetBranchStatus("eventNumber",1) ;
tree->SetBranchStatus("runNumber",1) ;
tree->SetBranchStatus("GpsTime",1) ;
tree->SetBranchStatus("EventInSequence",1) ;
tree->SetBranchStatus("totCandidates",1) ;
tree->SetBranchStatus("Polarity",1) ;
//keep the TCKinfo for KS (2012early/late) samples:
if(Kst2Kspiplus)// && year == "2012" && PHSP)
tree->SetBranchStatus("HLT1TCK", 1);
//activate cone variables
if(IsDST(year,MC,ReferenceChannel,PHSP) || Inc || (B0 && ReferenceChannel) || K1 ){
tree->SetBranchStatus("*_cp*",1);
tree->SetBranchStatus("*deltaEta*",1);
tree->SetBranchStatus("*deltaPhi*",1);
tree->SetBranchStatus("*cmult*",1);
tree->SetBranchStatus("*pxasy*",1);
tree->SetBranchStatus("*pyasy*",1);
tree->SetBranchStatus("*pzasy*",1);
tree->SetBranchStatus("*ptasy*",1);
tree->SetBranchStatus("*pasy*",1);
}
tree->SetBranchStatus("*M",1) ;
tree->SetBranchStatus("*MERR",1) ;
tree->SetBranchStatus("*MM",1) ;
tree->SetBranchStatus("*MMERR",1) ;
tree->SetBranchStatus( "*P", 1 );
tree->SetBranchStatus( "*PX", 1 );
tree->SetBranchStatus( "*PY", 1);
tree->SetBranchStatus( "*PZ", 1);
tree->SetBranchStatus( "*PE", 1);
tree->SetBranchStatus( "*PT", 1 );
tree->SetBranchStatus( "*TAU", 1 );
tree->SetBranchStatus( "*ETA", 1 ); //This should be fixed now
tree->SetBranchStatus( "*PHI", 1 );
tree->SetBranchStatus( "*FD*", 1 );
tree->SetBranchStatus( "*IP*", 1 );
tree->SetBranchStatus( "*IPCHI2_OWNPV", 1 );
tree->SetBranchStatus( "*FDCHI2_OWNPV",1 );
tree->SetBranchStatus( "*OWNPV*", 1 );
tree->SetBranchStatus( "*DIRA_OWNPV",1);
tree->SetBranchStatus( "*ENDVERTEX_CHI2",1 );
tree->SetBranchStatus( "*PV_Z*",1 );
tree->SetBranchStatus( "*ENDVERTEX*",1 );
tree->SetBranchStatus( "*DOCA*",1 ); //This should be fixed now
tree->SetBranchStatus( "*TRACK_CHI2NDOF",1) ;
tree->SetBranchStatus( "*TRACK_Type",1) ;
tree->SetBranchStatus( "*chi2",1 );
tree->SetBranchStatus( "*nDOF",1 );
tree->SetBranchStatus( "*status",1 );
tree->SetBranchStatus( "*ctau",1 );
tree->SetBranchStatus( "*isMuon",1 );
tree->SetBranchStatus( "*hasRich",1 );
tree->SetBranchStatus( "*PID*",1 );
tree->SetBranchStatus( "*_ID",1 );
tree->SetBranchStatus( "*Phi",1 );
tree->SetBranchStatus( "*ThetaL",1 );
tree->SetBranchStatus( "*ThetaK",1 );
if(MC){
tree->SetBranchStatus("*TRUE*",1) ;
tree->SetBranchStatus("*BKGCAT",1) ;
tree->SetBranchStatus("*MOTHER_ID",1) ;
tree->SetBranchStatus( "*CosTheta", 1 );
}
if (Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
tree->SetBranchStatus("gamma*ProbNN*",0) ; //Turn off gamma ProbNN
tree->SetBranchStatus("*_CL",1);
tree->SetBranchStatus("gamma*CaloHypo_E",1); //possibly use for tests, doesn't hurt to include it
tree->SetBranchStatus("gamma*CaloHypo_X",1);
tree->SetBranchStatus("gamma*CaloHypo_Y",1);
tree->SetBranchStatus("gamma*CaloHypo_Z",1);
}
//DTF convergence
if(UseDTF){
tree->SetBranchStatus("B_plus_DTF_status",1);
tree->SetBranchStatus("B_plus_DTF_chi2",1);
}
//define variables
ULong64_t nEventNumber;
Int_t nSubSetNumber; //split the subset into three subsets for BDT test&training, obtimization, application
TLorentzVector LorVec_B_plus;
TLorentzVector LorVec_K_star_plus;
TLorentzVector LorVec_K_plus;
TLorentzVector LorVec_pi_zero_merged;
TLorentzVector LorVec_pi_zero_resolved;
TLorentzVector LorVec_K_short;
TLorentzVector LorVec_pi_plus;
TLorentzVector LorVec_KS_pi_plus;
TLorentzVector LorVec_KS_pi_minus;
TLorentzVector LorVec_mu_minus;
TLorentzVector LorVec_mu_plus;
TLorentzVector LorVec_gamma1;
TLorentzVector LorVec_gamma2;
TLorentzVector LorVec_DiGamma;
//potentially misIDed particles:
TLorentzVector LorVec_pi_plus_MisIDed;
TLorentzVector LorVec_mu_plus_MisIDed;
TLorentzVector LorVec_B_plus_From_DoubleMisIDed;
TLorentzVector LorVec_Jpsi_From_DoubleMisIDed;
TLorentzVector LorVec_K_star_plus_From_DoubleMisIDed;
//recombined particles:
TLorentzVector LorVec_DiMuon;
TLorentzVector LorVec_KaonPion;
TLorentzVector LorVec_ReconB;
TLorentzVector LorVec_ReconBPiZero;
TLorentzVector LorVec_KshortMuMu;
TLorentzVector LorVec_KplusMuMu;
///B+
Double_t B_plus_M = 0.;
Double_t B_plus_MM = 0.;
Double_t B_plus_PX = 0.;
Double_t B_plus_PY = 0.;
Double_t B_plus_PZ = 0.;
Double_t B_plus_PE = 0.;
Double_t B_plus_PT = 0.;
Float_t f_B_plus_M_DTF[100];
Float_t f_B_plus_MERR_DTF[100];
Double_t B_plus_M_DTF = 0.;
Double_t B_plus_MERR_DTF = 0.;
Double_t B_plus_P_DTF = 0.;
Double_t B_plus_PX_DTF = 0.;
Double_t B_plus_PY_DTF = 0.;
Double_t B_plus_PZ_DTF = 0.;
Double_t B_plus_PE_DTF = 0.;
Double_t B_plus_PT_DTF = 0.;
Double_t B_plus_ETA_DTF = 0.;
Double_t B_plus_PHI_DTF = 0.;
///K*+
Double_t K_star_plus_M = 0.;
Double_t K_star_plus_MM = 0.;
Double_t K_star_plus_PX = 0.;
Double_t K_star_plus_PY = 0.;
Double_t K_star_plus_PZ = 0.;
Double_t K_star_plus_PE = 0.;
Double_t K_star_plus_PT = 0.;
Float_t f_K_star_plus_M_DTF[100];
Double_t K_star_plus_M_DTF = 0.;
Double_t K_star_plus_P_DTF = 0.;
Double_t K_star_plus_PX_DTF = 0.;
Double_t K_star_plus_PY_DTF = 0.;
Double_t K_star_plus_PZ_DTF = 0.;
Double_t K_star_plus_PE_DTF = 0.;
Double_t K_star_plus_PT_DTF = 0.;
Double_t K_star_plus_ETA_DTF= 0.;
Double_t K_star_plus_PHI_DTF= 0.;
///K+
Double_t K_plus_PX = 0.;
Double_t K_plus_PY = 0.;
Double_t K_plus_PZ = 0.;
Double_t K_plus_PE = 0.;
Double_t K_plus_PT = 0.;
Float_t f_K_plus_PX_DTF[100];
Float_t f_K_plus_PY_DTF[100];
Float_t f_K_plus_PZ_DTF[100];
Float_t f_K_plus_PE_DTF[100];
Float_t f_K_plus_ID_DTF[100];
Double_t K_plus_M_DTF = 0.; //Just pdg mass
Double_t K_plus_P_DTF = 0.;
Double_t K_plus_PX_DTF = 0.;
Double_t K_plus_PY_DTF = 0.;
Double_t K_plus_PZ_DTF = 0.;
Double_t K_plus_PE_DTF = 0.;
Double_t K_plus_PT_DTF = 0.;
Double_t K_plus_ETA_DTF= 0.;
Double_t K_plus_PHI_DTF= 0.;
Int_t K_plus_ID_DTF = 0;
Double_t K_plus_ProbNNk = 0.;
///pi0 (merged)
Double_t pi_zero_merged_M = 0.;
Double_t pi_zero_merged_PX = 0.;
Double_t pi_zero_merged_PY = 0.;
Double_t pi_zero_merged_PZ = 0.;
Double_t pi_zero_merged_PE = 0.;
Double_t pi_zero_merged_PT = 0.;
///pi0 (resolved)
Double_t pi_zero_resolved_M = 0.;
Double_t pi_zero_resolved_PX = 0.;
Double_t pi_zero_resolved_PY = 0.;
Double_t pi_zero_resolved_PZ = 0.;
Double_t pi_zero_resolved_PE = 0.;
Double_t pi_zero_resolved_PT = 0.;
Float_t f_pi_zero_resolved_M_DTF[100];
Double_t pi_zero_resolved_M_DTF = 0.;
Double_t pi_zero_resolved_P_DTF = 0.;
Double_t pi_zero_resolved_PX_DTF = 0.;
Double_t pi_zero_resolved_PY_DTF = 0.;
Double_t pi_zero_resolved_PZ_DTF = 0.;
Double_t pi_zero_resolved_PT_DTF = 0.;
Double_t pi_zero_resolved_PE_DTF = 0.;
Double_t pi_zero_resolved_ETA_DTF = 0.;
Double_t pi_zero_resolved_PHI_DTF = 0.;
///Kshort
Double_t K_short_PX = 0.;
Double_t K_short_PY = 0.;
Double_t K_short_PZ = 0.;
Double_t K_short_PE = 0.;
Double_t K_short_PT = 0.;
Float_t f_K_short_M_DTF[100];
Double_t K_short_M_DTF = 0.;
Double_t K_short_PT_DTF = 0.;
Double_t K_short_ETA_DTF = 0.;
Double_t K_short_PHI_DTF = 0.;
///pi+ (KS)
Double_t Ks_pi_plus_PX = 0.;
Double_t Ks_pi_plus_PY = 0.;
Double_t Ks_pi_plus_PZ = 0.;
Double_t Ks_pi_plus_PE = 0.;
Double_t Ks_pi_plus_PT = 0.;
Float_t f_Ks_pi_plus_PX_DTF[100];
Float_t f_Ks_pi_plus_PY_DTF[100];
Float_t f_Ks_pi_plus_PZ_DTF[100];
Float_t f_Ks_pi_plus_PE_DTF[100];
Float_t f_Ks_pi_plus_ID_DTF[100];
Double_t Ks_pi_plus_PX_DTF = 0.;
Double_t Ks_pi_plus_PY_DTF = 0.;
Double_t Ks_pi_plus_PZ_DTF = 0.;
Double_t Ks_pi_plus_PE_DTF = 0.;
Double_t Ks_pi_plus_PT_DTF = 0.;
Double_t Ks_pi_plus_ETA_DTF = 0.;
Double_t Ks_pi_plus_PHI_DTF = 0.;
Int_t Ks_pi_plus_ID_DTF = 0;
///pi- (KS)
Double_t Ks_pi_minus_PX = 0.;
Double_t Ks_pi_minus_PY = 0.;
Double_t Ks_pi_minus_PZ = 0.;
Double_t Ks_pi_minus_PE = 0.;
Double_t Ks_pi_minus_PT = 0.;
Float_t f_Ks_pi_minus_PX_DTF[100];
Float_t f_Ks_pi_minus_PY_DTF[100];
Float_t f_Ks_pi_minus_PZ_DTF[100];
Float_t f_Ks_pi_minus_PE_DTF[100];
Float_t f_Ks_pi_minus_ID_DTF[100];
Double_t Ks_pi_minus_PX_DTF = 0.;
Double_t Ks_pi_minus_PY_DTF = 0.;
Double_t Ks_pi_minus_PZ_DTF = 0.;
Double_t Ks_pi_minus_PE_DTF = 0.;
Double_t Ks_pi_minus_PT_DTF = 0.;
Double_t Ks_pi_minus_ETA_DTF = 0.;
Double_t Ks_pi_minus_PHI_DTF = 0.;
Int_t Ks_pi_minus_ID_DTF = 0;
///pi+
Double_t pi_plus_PX = 0.;
Double_t pi_plus_PY = 0.;
Double_t pi_plus_PZ = 0.;
Double_t pi_plus_PE = 0.;
Double_t pi_plus_PT = 0.;
Float_t f_pi_plus_PX_DTF[100];
Float_t f_pi_plus_PY_DTF[100];
Float_t f_pi_plus_PZ_DTF[100];
Float_t f_pi_plus_PE_DTF[100];
Float_t f_pi_plus_ID_DTF[100];
Double_t pi_plus_PX_DTF = 0.;
Double_t pi_plus_PY_DTF = 0.;
Double_t pi_plus_PZ_DTF = 0.;
Double_t pi_plus_PE_DTF = 0.;
Double_t pi_plus_PT_DTF = 0.;
Double_t pi_plus_ETA_DTF= 0.;
Double_t pi_plus_PHI_DTF= 0.;
Int_t pi_plus_ID_DTF = 0;
///mu-
Double_t mu_minus_PX = 0.;
Double_t mu_minus_PY = 0.;
Double_t mu_minus_PZ = 0.;
Double_t mu_minus_PE = 0.;
Double_t mu_minus_PT = 0.;
Float_t f_mu_minus_PX_DTF[100];
Float_t f_mu_minus_PY_DTF[100];
Float_t f_mu_minus_PZ_DTF[100];
Float_t f_mu_minus_PE_DTF[100];
Float_t f_mu_minus_ID_DTF[100];
Double_t mu_minus_M_DTF = 0.;
Double_t mu_minus_P_DTF = 0.;
Double_t mu_minus_PX_DTF = 0.;
Double_t mu_minus_PY_DTF = 0.;
Double_t mu_minus_PZ_DTF = 0.;
Double_t mu_minus_PE_DTF = 0.;
Double_t mu_minus_PT_DTF = 0.;
Double_t mu_minus_ETA_DTF= 0.;
Double_t mu_minus_PHI_DTF= 0.;
Int_t mu_minus_ID_DTF = 0;
///mu+
Double_t mu_plus_PX = 0.;
Double_t mu_plus_PY = 0.;
Double_t mu_plus_PZ = 0.;
Double_t mu_plus_PE = 0.;
Double_t mu_plus_PT = 0.;
Float_t f_mu_plus_PX_DTF[100];
Float_t f_mu_plus_PY_DTF[100];
Float_t f_mu_plus_PZ_DTF[100];
Float_t f_mu_plus_PE_DTF[100];
Float_t f_mu_plus_ID_DTF[100];
Double_t mu_plus_M_DTF = 0.;
Double_t mu_plus_P_DTF = 0.;
Double_t mu_plus_PX_DTF = 0.;
Double_t mu_plus_PY_DTF = 0.;
Double_t mu_plus_PZ_DTF = 0.;
Double_t mu_plus_PE_DTF = 0.;
Double_t mu_plus_PT_DTF = 0.;
Double_t mu_plus_ETA_DTF= 0.;
Double_t mu_plus_PHI_DTF= 0.;
Int_t mu_plus_ID_DTF = 0;
///gamma1
Double_t gamma1_PX = 0.;
Double_t gamma1_PY = 0.;
Double_t gamma1_PZ = 0.;
Double_t gamma1_PE = 0.;
Double_t gamma1_PT = 0.;
Float_t f_gamma1_PX_DTF[100];
Float_t f_gamma1_PY_DTF[100];
Float_t f_gamma1_PZ_DTF[100];
Float_t f_gamma1_PE_DTF[100];
Float_t f_gamma1_ID_DTF[100];
Double_t gamma1_M_DTF = 0.;
Double_t gamma1_P_DTF = 0.;
Double_t gamma1_PX_DTF = 0.;
Double_t gamma1_PY_DTF = 0.;
Double_t gamma1_PZ_DTF = 0.;
Double_t gamma1_PE_DTF = 0.;
Double_t gamma1_PT_DTF = 0.;
Double_t gamma1_ETA_DTF= 0.;
Double_t gamma1_PHI_DTF= 0.;
Int_t gamma1_ID_DTF = 0;
///gamma2
Double_t gamma2_PX = 0.;
Double_t gamma2_PY = 0.;
Double_t gamma2_PZ = 0.;
Double_t gamma2_PE = 0.;
Double_t gamma2_PT = 0.;
Float_t f_gamma2_PX_DTF[100];
Float_t f_gamma2_PY_DTF[100];
Float_t f_gamma2_PZ_DTF[100];
Float_t f_gamma2_PE_DTF[100];
Float_t f_gamma2_ID_DTF[100];
Double_t gamma2_M_DTF = 0.;
Double_t gamma2_P_DTF = 0.;
Double_t gamma2_PX_DTF = 0.;
Double_t gamma2_PY_DTF = 0.;
Double_t gamma2_PZ_DTF = 0.;
Double_t gamma2_PE_DTF = 0.;
Double_t gamma2_PT_DTF = 0.;
Double_t gamma2_ETA_DTF= 0.;
Double_t gamma2_PHI_DTF= 0.;
Int_t gamma2_ID_DTF = 0;
//vertex quality:
Double_t B_plus_DIRA_OWNPV = 0.;
Double_t B_plus_FDCHI2_OWNPV = 0.;
Double_t B_plus_IPCHI2_OWNPV = 0.;
Double_t K_star_plus_IPCHI2_OWNPV = 0.;
Double_t mu_plus_IPCHI2_OWNPV = 0.;
Double_t mu_minus_IPCHI2_OWNPV = 0.;
Double_t K_star_plus_FDCHI2_OWNPV = 0.;
//particle identification
Double_t K_plus_PIDK = 0.;
Double_t pi_plus_PIDK = 0.;
Double_t mu_plus_ProbNNmu = 0.;
Double_t mu_minus_ProbNNmu = 0.;
//track type variables
Int_t Ks_pi_plus_TRACK_Type;
Int_t Ks_pi_minus_TRACK_Type;
//confidence level
Double_t gamma1_CL = 0.;
Double_t gamma2_CL = 0.;
//Angles
Double_t ctl = 0.;
Double_t ctk = 0.;
Double_t phi = 0.;
//trigger lines
//L0
Bool_t B_plus_L0MuonDecision_TOS = false;
//HLT1
Bool_t B_plus_Hlt1TrackAllL0Decision_TOS = false;
Bool_t B_plus_Hlt1TrackMuonDecision_TOS = false;
//HLT2
Bool_t B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2Topo2BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2Topo3BodyBBDTDecision_TOS = false;
// Bool_t B_plus_Hlt2SingleMuonDecision_TOS = false;
Bool_t B_plus_Hlt2DiMuonDetachedDecision_TOS = false;
//isMuon variables
Bool_t mu_minus_isMuon = false;
Bool_t mu_plus_isMuon = false;
//hasRich variables
Bool_t pi_plus_hasRich = false;
Bool_t K_plus_hasRich = false;
//DTF status
Float_t B_plus_DTF_status[100];
Float_t B_plus_DTF_chi2[100];
///////////////////////////////////////
/// LINK VARIABLES TO BRANCHES
///////////////////////////////////////
//event numbers
tree -> SetBranchAddress( "eventNumber" , &nEventNumber );
//trigger decisions
tree -> SetBranchAddress( "B_plus_L0MuonDecision_TOS" , &B_plus_L0MuonDecision_TOS );
if(year == "2011" || year == "2012"){
tree -> SetBranchAddress( "B_plus_Hlt1TrackAllL0Decision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );
}
else if(year == "2015" || year == "2016" || year == "2017" || year == "2018"){
tree -> SetBranchAddress( "B_plus_Hlt1TrackMVADecision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );
if(year == "2016" || year == "2017" || year == "2018"){ //found to be FALSE for all 2015 data
tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS );
}
}
else{
coutERROR("Incorrect year given for the assosiation of HLT1 lines with branches! Year given: " + year);
return 0;
}
tree -> SetBranchAddress( "B_plus_Hlt1TrackMuonDecision_TOS" , &B_plus_Hlt1TrackMuonDecision_TOS );
// tree -> SetBranchAddress( "B_plus_Hlt2SingleMuonDecision_TOS" , &B_plus_Hlt2SingleMuonDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2DiMuonDetachedDecision_TOS" , &B_plus_Hlt2DiMuonDetachedDecision_TOS );
//particle identification (PID)
tree -> SetBranchAddress( "mu_minus_isMuon", &mu_minus_isMuon );
tree -> SetBranchAddress( "mu_plus_isMuon", &mu_plus_isMuon );
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
tree -> SetBranchAddress( "K_plus_hasRich", &K_plus_hasRich );
if(UsePIDgen && PHSP){ //TODO add '&& MC' and '&& ReferenceChannel' once the PIDgen tuples are available for these sets
tree -> SetBranchAddress( "K_plus_PIDK_corr", &K_plus_PIDK );
tree -> SetBranchAddress( "K_plus_ProbNNk", &K_plus_ProbNNk );
tree -> SetBranchAddress( "mu_plus_ProbNNmu_corr", &mu_plus_ProbNNmu );
tree -> SetBranchAddress( "mu_minus_ProbNNmu_corr", &mu_minus_ProbNNmu );
}
else{
//tree -> SetBranchAddress( "K_plus_PIDK", &K_plus_PIDK ); replaced by probNN
tree -> SetBranchAddress( "K_plus_ProbNNk", &K_plus_ProbNNk );
tree -> SetBranchAddress( "mu_plus_ProbNNmu", &mu_plus_ProbNNmu );
tree -> SetBranchAddress( "mu_minus_ProbNNmu", &mu_minus_ProbNNmu );
}
}
if(Kst2Kspiplus){
tree -> SetBranchAddress( "pi_plus_hasRich", &pi_plus_hasRich );
if(UsePIDgen && PHSP){ //TODO add '&& MC' and '&& ReferenceChannel' once the PIDgen tuples are available for these sets
tree -> SetBranchAddress( "pi_plus_PIDK_corr", &pi_plus_PIDK );
tree -> SetBranchAddress( "mu_plus_ProbNNmu_corr", &mu_plus_ProbNNmu );
tree -> SetBranchAddress( "mu_minus_ProbNNmu_corr", &mu_minus_ProbNNmu );
}
else{
tree -> SetBranchAddress( "pi_plus_PIDK", &pi_plus_PIDK );
tree -> SetBranchAddress( "mu_plus_ProbNNmu", &mu_plus_ProbNNmu );
tree -> SetBranchAddress( "mu_minus_ProbNNmu", &mu_minus_ProbNNmu );
}
}
//particle masses and momenta
tree -> SetBranchAddress( "B_plus_M" , &B_plus_M );
tree -> SetBranchAddress( "B_plus_MM" , &B_plus_MM );
tree -> SetBranchAddress( "B_plus_PX" , &B_plus_PX );
tree -> SetBranchAddress( "B_plus_PY" , &B_plus_PY );
tree -> SetBranchAddress( "B_plus_PZ" , &B_plus_PZ );
tree -> SetBranchAddress( "B_plus_PE" , &B_plus_PE );
tree -> SetBranchAddress( "B_plus_PT" , &B_plus_PT );
tree -> SetBranchAddress( "B_plus_DTF_M" , &f_B_plus_M_DTF);
tree -> SetBranchAddress( "B_plus_DTF_MERR", &f_B_plus_MERR_DTF);
tree -> SetBranchAddress( "K_star_plus_M" , &K_star_plus_M );
tree -> SetBranchAddress( "K_star_plus_MM" , &K_star_plus_MM );
tree -> SetBranchAddress( "K_star_plus_PX" , &K_star_plus_PX );
tree -> SetBranchAddress( "K_star_plus_PY" , &K_star_plus_PY );
tree -> SetBranchAddress( "K_star_plus_PZ" , &K_star_plus_PZ );
tree -> SetBranchAddress( "K_star_plus_PE" , &K_star_plus_PE );
tree -> SetBranchAddress( "K_star_plus_PT" , &K_star_plus_PT );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_M" , &f_K_star_plus_M_DTF );
if(Kst2Kpluspi0Merged){
tree -> SetBranchAddress( "pi_zero_merged_PX" , &pi_zero_merged_PX );
tree -> SetBranchAddress( "pi_zero_merged_PY" , &pi_zero_merged_PY );
tree -> SetBranchAddress( "pi_zero_merged_PZ" , &pi_zero_merged_PZ );
tree -> SetBranchAddress( "pi_zero_merged_PE" , &pi_zero_merged_PE );
tree -> SetBranchAddress( "pi_zero_merged_PT" , &pi_zero_merged_PT );
tree -> SetBranchAddress( "pi_zero_merged_M" , &pi_zero_merged_M );
}
if(Kst2Kpluspi0Resolved){
tree -> SetBranchAddress( "pi_zero_resolved_PX" , &pi_zero_resolved_PX );
tree -> SetBranchAddress( "pi_zero_resolved_PY" , &pi_zero_resolved_PY );
tree -> SetBranchAddress( "pi_zero_resolved_PZ" , &pi_zero_resolved_PZ );
tree -> SetBranchAddress( "pi_zero_resolved_PE" , &pi_zero_resolved_PE );
tree -> SetBranchAddress( "pi_zero_resolved_PT" , &pi_zero_resolved_PT );
tree -> SetBranchAddress( "pi_zero_resolved_M" , &pi_zero_resolved_M );
tree -> SetBranchAddress( "B_plus_DTF_pi0_M" , &f_pi_zero_resolved_M_DTF );
tree -> SetBranchAddress( "gamma1_PX" , &gamma1_PX );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PX" , &f_gamma1_PX_DTF );
tree -> SetBranchAddress( "gamma1_PY" , &gamma1_PY );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PY" , &f_gamma1_PY_DTF );
tree -> SetBranchAddress( "gamma1_PZ" , &gamma1_PZ );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PZ" , &f_gamma1_PZ_DTF );
tree -> SetBranchAddress( "gamma1_PE" , &gamma1_PE );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_PE" , &f_gamma1_PE_DTF );
tree -> SetBranchAddress( "gamma1_PT" , &gamma1_PT );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_0_ID" , &f_gamma1_ID_DTF );
tree -> SetBranchAddress( "gamma2_PX" , &gamma2_PX );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PX" , &f_gamma2_PX_DTF );
tree -> SetBranchAddress( "gamma2_PY" , &gamma2_PY );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PY" , &f_gamma2_PY_DTF );
tree -> SetBranchAddress( "gamma2_PZ" , &gamma2_PZ );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PZ" , &f_gamma2_PZ_DTF );
tree -> SetBranchAddress( "gamma2_PE" , &gamma2_PE );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_PE" , &f_gamma2_PE_DTF );
tree -> SetBranchAddress( "gamma2_PT" , &gamma2_PT );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_pi0_gamma_ID" , &f_gamma2_ID_DTF );
tree -> SetBranchAddress( "gamma1_CL" , &gamma1_CL );
tree -> SetBranchAddress( "gamma2_CL" , &gamma2_CL );
}
if(Kst2Kpluspi0Merged || Kst2Kpluspi0Resolved){
tree -> SetBranchAddress( "K_plus_PX" , &K_plus_PX );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PX" , &f_K_plus_PX_DTF );
tree -> SetBranchAddress( "K_plus_PY" , &K_plus_PY );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PY" , &f_K_plus_PY_DTF );
tree -> SetBranchAddress( "K_plus_PZ" , &K_plus_PZ );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PZ" , &f_K_plus_PZ_DTF );
tree -> SetBranchAddress( "K_plus_PE" , &K_plus_PE );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_PE" , &f_K_plus_PE_DTF );
tree -> SetBranchAddress( "K_plus_PT" , &K_plus_PT );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_Kplus_ID" , &f_K_plus_ID_DTF );
}
if(Kst2Kspiplus){
tree -> SetBranchAddress( "K_short_PX" , &K_short_PX );
tree -> SetBranchAddress( "K_short_PY" , &K_short_PY );
tree -> SetBranchAddress( "K_short_PZ" , &K_short_PZ );
tree -> SetBranchAddress( "K_short_PE" , &K_short_PE );
tree -> SetBranchAddress( "K_short_PT" , &K_short_PT );
tree -> SetBranchAddress( "B_plus_DTF_KS0_M" , &f_K_short_M_DTF );
tree -> SetBranchAddress( "pi_plus_PX" , &pi_plus_PX );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PX" , &f_pi_plus_PX_DTF );
tree -> SetBranchAddress( "pi_plus_PY" , &pi_plus_PY );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PY" , &f_pi_plus_PY_DTF );
tree -> SetBranchAddress( "pi_plus_PZ" , &pi_plus_PZ );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PZ" , &f_pi_plus_PZ_DTF );
tree -> SetBranchAddress( "pi_plus_PE" , &pi_plus_PE );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_PE" , &f_pi_plus_PE_DTF );
tree -> SetBranchAddress( "pi_plus_PT" , &pi_plus_PT );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_piplus_ID" , &f_pi_plus_ID_DTF );
tree -> SetBranchAddress( "Ks_pi_minus_PX" , &Ks_pi_minus_PX );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PX" , &f_Ks_pi_minus_PX_DTF );
tree -> SetBranchAddress( "Ks_pi_minus_PY" , &Ks_pi_minus_PY );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PY" , &f_Ks_pi_minus_PY_DTF );
tree -> SetBranchAddress( "Ks_pi_minus_PZ" , &Ks_pi_minus_PZ );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PZ" , &f_Ks_pi_minus_PZ_DTF );
tree -> SetBranchAddress( "Ks_pi_minus_PE" , &Ks_pi_minus_PE );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_PE" , &f_Ks_pi_minus_PE_DTF );
tree -> SetBranchAddress( "Ks_pi_minus_PT" , &Ks_pi_minus_PT );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_ID" , &f_Ks_pi_minus_ID_DTF );
tree -> SetBranchAddress( "Ks_pi_plus_PX" , &Ks_pi_plus_PX );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PX" , &f_Ks_pi_plus_PX_DTF );
tree -> SetBranchAddress( "Ks_pi_plus_PY" , &Ks_pi_plus_PY );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PY" , &f_Ks_pi_plus_PY_DTF );
tree -> SetBranchAddress( "Ks_pi_plus_PZ" , &Ks_pi_plus_PZ );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PZ" , &f_Ks_pi_plus_PZ_DTF );
tree -> SetBranchAddress( "Ks_pi_plus_PE" , &Ks_pi_plus_PE );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_PE" , &f_Ks_pi_plus_PE_DTF );
tree -> SetBranchAddress( "Ks_pi_plus_PT" , &Ks_pi_plus_PT );
tree -> SetBranchAddress( "B_plus_DTF_Kst_892_plus_KS0_piplus_0_ID" , &f_Ks_pi_plus_ID_DTF );
}
tree -> SetBranchAddress( "mu_minus_PX" , &mu_minus_PX );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PX" , &f_mu_minus_PX_DTF );
tree -> SetBranchAddress( "mu_minus_PY" , &mu_minus_PY );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PY" , &f_mu_minus_PY_DTF );
tree -> SetBranchAddress( "mu_minus_PZ" , &mu_minus_PZ );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PZ" , &f_mu_minus_PZ_DTF );
tree -> SetBranchAddress( "mu_minus_PE" , &mu_minus_PE );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_PE" , &f_mu_minus_PE_DTF );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_0_ID" , &f_mu_minus_ID_DTF );
tree -> SetBranchAddress( "mu_minus_PT" , &mu_minus_PT );
tree -> SetBranchAddress( "mu_plus_PX" , &mu_plus_PX );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PX" , &f_mu_plus_PX_DTF );
tree -> SetBranchAddress( "mu_plus_PY" , &mu_plus_PY );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PY" , &f_mu_plus_PY_DTF );
tree -> SetBranchAddress( "mu_plus_PZ" , &mu_plus_PZ );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PZ" , &f_mu_plus_PZ_DTF );
tree -> SetBranchAddress( "mu_plus_PE" , &mu_plus_PE );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_PE" , &f_mu_plus_PE_DTF );
tree -> SetBranchAddress( "mu_plus_PT" , &mu_plus_PT );
tree -> SetBranchAddress( "B_plus_DTF_J_psi_1S_muminus_ID" , &f_mu_plus_ID_DTF );
tree -> SetBranchAddress( "B_plus_DTF_status", &B_plus_DTF_status);
tree -> SetBranchAddress( "B_plus_DTF_chi2", &B_plus_DTF_chi2);
tree -> SetBranchAddress( "B_plus_DIRA_OWNPV" , &B_plus_DIRA_OWNPV );
tree -> SetBranchAddress( "B_plus_FDCHI2_OWNPV" , &B_plus_FDCHI2_OWNPV );
tree -> SetBranchAddress( "B_plus_IPCHI2_OWNPV" , &B_plus_IPCHI2_OWNPV );
tree -> SetBranchAddress( "K_star_plus_IPCHI2_OWNPV" , &K_star_plus_IPCHI2_OWNPV );
tree -> SetBranchAddress( "mu_plus_IPCHI2_OWNPV" , &mu_plus_IPCHI2_OWNPV );
tree -> SetBranchAddress( "mu_minus_IPCHI2_OWNPV" , &mu_minus_IPCHI2_OWNPV );
tree -> SetBranchAddress( "K_star_plus_FDCHI2_OWNPV" , &K_star_plus_FDCHI2_OWNPV );
//Angles
tree -> SetBranchAddress( "B_plus_ThetaL", &ctl);
tree -> SetBranchAddress( "B_plus_ThetaK", &ctk);
tree -> SetBranchAddress( "B_plus_Phi", &phi);
//track types
if(Kst2Kspiplus){
tree -> SetBranchAddress( "Ks_pi_plus_TRACK_Type" , &Ks_pi_plus_TRACK_Type );
tree -> SetBranchAddress( "Ks_pi_minus_TRACK_Type" , &Ks_pi_minus_TRACK_Type );
}
//-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//Create output file
TFile* output = 0;
TFile* histos = 0;
// Kst2Kpluspi0Resolved case
if(Kst2Kpluspi0Resolved){
if (B0 || K1 || Inc){
output = new TFile(GetInputFileBkg(year,magnet,true,ReferenceChannel,B0,K1,Inc).c_str(),"RECREATE");
histos = new TFile(GetHistosFileBkg(year,magnet,ReferenceChannel,B0,K1,Inc).c_str(),"RECREATE");
}
else{
output = new TFile(GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
histos = new TFile(GetHistosFile(year,magnet,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
}
}
// Kst2Kpluspi0Merged case
if(Kst2Kpluspi0Merged){
output = new TFile(GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
histos = new TFile(GetHistosFile(year,magnet,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
}
// Kst2Kspiplus case
if(Kst2Kspiplus){
output = new TFile(GetInputFile(year,magnet,true,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
histos = new TFile(GetHistosFile(year,magnet,MC,ReferenceChannel,PHSP,smallSample).c_str(),"RECREATE");
}
output->cd();
TTree * summary_tree = tree->CloneTree(0);
//-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
int nCounterTrigger = 0;
int nCounterKinematic = 0;
int nCounterPrV = 0;
int nCounterPoV = 0;
int nCounterPID = 0;
int nCounterJpsiResonance = 0;
int nCounterKstarMassWindow = 0;
int nCounterB0KsmumuMass = 0;
int nCounterDiMuonMassInQ2Range = 0;
int n_DoubleMisID = 0;
histos->cd();
TH1 * h_JpsiMass = new TH1D("h_DiMuonMass" , "mass spectra of combined muons", 250, 0, 5000);
h_JpsiMass->GetXaxis()->SetTitle("DiMuonMass [MeV]");
h_JpsiMass->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_JpsiMassCut = new TH1D("h_DiMuonMassw/oResonances" , "mass spectra of combined muons after resonance cut", 250, 0, 5000);
h_JpsiMassCut->GetXaxis()->SetTitle("DiMuonMass [MeV]");
h_JpsiMassCut->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_JpsiMassResonance = new TH1D("h_DiMuonMassOnlyJpsiResonance" , "mass spectra of combined muons inside the Jpsi band", 75, 2500, 4000);
h_JpsiMassResonance->GetXaxis()->SetTitle("DiMuonMass [MeV]");
h_JpsiMassResonance->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_Psi2sMassResonance = new TH1D("h_DiMuonMassOnlyPsi2Resonance" , "mass spectra of combined muons inside the Psi2 band", 75, 2500, 4000);
h_Psi2sMassResonance->GetXaxis()->SetTitle("DiMuonMass [MeV]");
h_Psi2sMassResonance->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_CharmoniumResonance = new TH1D("h_DiMuonMassCharmoniumResonance" , "mass spectra of combined muons inside the charmonium bands", 75, 2500, 4000);
h_CharmoniumResonance->GetXaxis()->SetTitle("DiMuonMass [MeV]");
h_CharmoniumResonance->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_BplusMass_B02KsmumuMass = new TH1D("h_BplusMass_B02KsmumuMass" , "mass spectra of reconstructed B0 from KS and di-muon", 175, 3500, 6000);
h_BplusMass_B02KsmumuMass->GetXaxis()->SetTitle("Reconstructed B0 mass [MeV]");
h_BplusMass_B02KsmumuMass->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_BplusMass_DTF = new TH1D("h_BplusMass_DTF" , "Mass spectra from DTF", 125, 4500, 7000);
h_BplusMass_DTF->GetXaxis()->SetTitle("B^{+} Mass from DTF [MeV]");
h_BplusMass_DTF->GetYaxis()->SetTitle("Counts / 20MeV");
TH2 * h_BplusMass_DiMuonMass = new TH2D("h_BplusMass_DiMuonMass" , "2D mass plot after pre-selection", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
h_BplusMass_DiMuonMass->GetXaxis()->SetTitle("B^{+} mass [MeV]");
h_BplusMass_DiMuonMass->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
TH2 * h_BplusMass_DiMuonMassRaw = new TH2D("h_BplusMass_DiMuonMassRaw" , "2D mass plot after stripping", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
h_BplusMass_DiMuonMassRaw->GetXaxis()->SetTitle("B^{+} mass [MeV]");
h_BplusMass_DiMuonMassRaw->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
TH2 * h_BplusMass_DiMuonMassDD = new TH2D("h_BplusMass_DiMuonMassDD" , "2D mass plot after pre-selection (DD tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
h_BplusMass_DiMuonMassDD->GetXaxis()->SetTitle("B^{+} mass [MeV]");
h_BplusMass_DiMuonMassDD->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
TH2 * h_BplusMass_DiMuonMassRawDD = new TH2D("h_BplusMass_DiMuonMassRawDD" , "2D mass plot after stripping (DD tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
h_BplusMass_DiMuonMassRawDD->GetXaxis()->SetTitle("B^{+} mass [MeV]");
h_BplusMass_DiMuonMassRawDD->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
TH2 * h_BplusMass_DiMuonMassLL = new TH2D("h_BplusMass_DiMuonMassLL" , "2D mass plot after pre-selection (LL tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
h_BplusMass_DiMuonMassLL->GetXaxis()->SetTitle("B^{+} mass [MeV]");
h_BplusMass_DiMuonMassLL->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
TH2 * h_BplusMass_DiMuonMassRawLL = new TH2D("h_BplusMass_DiMuonMassRawLL" , "2D mass plot after stripping (LL tracks)", 150 ,get_cut_B_plus_M_low(year), cut_B_plus_M_high, 150, 1000., 4500.);
h_BplusMass_DiMuonMassRawLL->GetXaxis()->SetTitle("B^{+} mass [MeV]");
h_BplusMass_DiMuonMassRawLL->GetYaxis()->SetTitle("#mu^{+}#mu^{-} mass [MeV]");
TH1 * h_ReconstructBmassInCharmoniumResonance = new TH1D("h_K*JpsiMassInCharmoniumResonance" , "mass spectra of K* and Jpsi inside the charmonium bands", 75, 5000, 6500);
h_ReconstructBmassInCharmoniumResonance->GetXaxis()->SetTitle("K*JpsiMass [MeV]");
h_ReconstructBmassInCharmoniumResonance->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_ReconstructBmassOffResonance = new TH1D("h_K*JpsiMassOutsideOffCharmoniumResonance" , "mass spectra of K* and DiMuon without the charmonium bands", 75, 5000, 6500);
h_ReconstructBmassOffResonance->GetXaxis()->SetTitle("K*mumuMass [MeV]");
h_ReconstructBmassOffResonance->GetYaxis()->SetTitle("Counts / 20MeV");
TH1 * h_KstarMass = new TH1D("h_KstarMass" , "mass spectra of combined kaon and pion", 90, 500, 1400);
h_KstarMass->GetXaxis()->SetTitle("KstarMass [MeV]");
h_KstarMass->GetYaxis()->SetTitle("Counts / 10MeV");
TH1 * h_KstarMasscut = new TH1D("h_KstarMassIn100MeV window" , "mass spectra of combined kaon and pion in 100MeV window", 90, 500, 1400);
h_KstarMasscut->GetXaxis()->SetTitle("KstarMass [MeV]");
h_KstarMasscut->GetYaxis()->SetTitle("Counts / 10MeV");
TH2 * h_KstarMass2D = new TH2D("h_KstarMass2D" , "mass spectrum of Kstarplus", 90, 500, 1400, 90, 500, 1400);
h_KstarMass2D->GetXaxis()->SetTitle("mass from tree [MeV]");
h_KstarMass2D->GetYaxis()->SetTitle("mass from K and pi [MeV]");
TH2 * h_diMuonReconstructedBMass2D = new TH2D("h_diMuonReconstructedBMass2D" , "mass spectrum of dimuon vs B+candidates", 150, 5000, 6500, 150, 2500, 4000);
TH2 * h_diMuonReconstructedBMass2DFixPiZero = new TH2D("h_diMuonReconstructedBMass2DFixPiZero" , "mass spectrum of dimuon vs B+candidates with fixed m_{#pi^0}", 150, 5000, 6500, 150, 2500, 4000);
h_diMuonReconstructedBMass2D->GetYaxis()->SetTitle("m_{#mu#mu} [MeV]");
if (Kst2Kspiplus)
h_diMuonReconstructedBMass2D->GetXaxis()->SetTitle("m_{K_{s}#pi^{+}#mu#mu} [MeV]");
else{
if (Kst2Kpluspi0Merged)
h_diMuonReconstructedBMass2D->GetXaxis()->SetTitle("m_{K_{+}#pi^{0}#mu#mu} [MeV]");
else
h_diMuonReconstructedBMass2D->GetXaxis()->SetTitle("m_{K_{+}#gamma#gamma#mu#mu} [MeV]");
h_diMuonReconstructedBMass2DFixPiZero->GetXaxis()->SetTitle("m_{K_{+}#pi^{0}#mu#mu} [MeV]");
}
//effect of PID cut on angular distribution
TH1 * h_BeforePIDcuts_ctl = new TH1D("h_BeforePIDcuts_ctl" , "cos(#theta_{L}) before PID cuts; cos(#theta_{L}); entries", 100, -1., +1.);
TH1 * h_BeforePIDcuts_ctk = new TH1D("h_BeforePIDcuts_ctk" , "cos(#theta_{K}) before PID cuts; cos(#theta_{K}); entries", 100, -1., +1.);
TH1 * h_BeforePIDcuts_phi = new TH1D("h_BeforePIDcuts_phi" , "#phi before PID cuts; #phi (rad.); entries", 100, -TMath::Pi(), +TMath::Pi());
TH1 * h_AfterPIDcuts_ctl = new TH1D("h_AfterPIDcuts_ctl" , "cos(#theta_{L}) after PID cuts; cos(#theta_{L}); entries", 100, -1., +1.);
TH1 * h_AfterPIDcuts_ctk = new TH1D("h_AfterPIDcuts_ctk" , "cos(#theta_{K}) after PID cuts; cos(#theta_{K}); entries", 100, -1., +1.);
TH1 * h_AfterPIDcuts_phi = new TH1D("h_AfterPIDcuts_phi" , "#phi after PID cuts; #phi (rad.); entries", 100, -TMath::Pi(), +TMath::Pi());
/// K_s pi+
TH2 * h_KshortDaughters = new TH2D("h_KshortDaughters" , "track types of the Kshort daughter pions", 10, -0.5, 9.5, 10, -0.5, 9.5);
h_KshortDaughters->GetXaxis()->SetTitle("pi minus");
h_KshortDaughters->GetYaxis()->SetTitle("pi plus");
TH2 * h_BplusM_vs_BplusM_DoubleMisIDed = new TH2D("h_BplusM_vs_BplusM_DoubleMisIDed" , "comparison of B+ mass versus B+ mass (mu+ and pi+ double misIDed", 100, 5000, 6000, 100, 5000, 6000);
h_BplusM_vs_BplusM_DoubleMisIDed->GetXaxis()->SetTitle("B_plus_M [MeV]");
h_BplusM_vs_BplusM_DoubleMisIDed->GetYaxis()->SetTitle("B_plus_M (double misIDed) [MeV]");
TH2 * h_BplusM_vs_JpsiM_DoubleMisIDed = new TH2D("h_BplusM_vs_JpsiM_DoubleMisIDed" , "overlay of J_psi mass and B+ mass (both mu+ and pi+ double misIDed", 50, 5000, 6000, 200, 0, 4000);
h_BplusM_vs_JpsiM_DoubleMisIDed->GetXaxis()->SetTitle("B_plus_M (double misIDed) [MeV]");
h_BplusM_vs_JpsiM_DoubleMisIDed->GetYaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");
TH2 * h_KstarM_vs_JpsiM_DoubleMisIDed = new TH2D("h_KstarM_vs_JpsiM_DoubleMisIDed" , "overlay of J_psi mass and K*+ mass (both mu+ and pi+ double misIDed", 125, 500, 3000, 200, 0, 4000);
h_KstarM_vs_JpsiM_DoubleMisIDed->GetXaxis()->SetTitle("K_star_mass (double misIDed) [MeV]");
h_KstarM_vs_JpsiM_DoubleMisIDed->GetYaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");
TH1 * h_Jpsi_Spectrum_MisIDed_Muon = new TH1D("h_Jpsi_Spectrum_MisIDed_Muon" , "Invariant mass of Jpsi with a pion misIDed as a muon", 200, 0, 4000);
h_Jpsi_Spectrum_MisIDed_Muon->GetXaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");
h_Jpsi_Spectrum_MisIDed_Muon->GetYaxis()->SetTitle("Candidates/20MeV");
TH1 * h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon = new TH1D("h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon" , "Invariant mass of Jpsi with a pion misIDed as a muon (Veto Resonant DiMuons)", 200, 0, 4000);
h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->GetXaxis()->SetTitle("J_psi_M (double misIDed) [MeV]");
h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->GetYaxis()->SetTitle("Candidates/20MeV");
TH1 * h_Q2_Spectrum_MisIDed_Muon = new TH1D("h_Q2_Spectrum_MisIDed_Muon" , "Q2 of non-resonant DiMuons with a pion misIDed as a muon", 100, 0, 20);
h_Q2_Spectrum_MisIDed_Muon->GetXaxis()->SetTitle("Q2 (double misIDed) [GeV^{2}]");
h_Q2_Spectrum_MisIDed_Muon->GetYaxis()->SetTitle("Candidates / 0.2 GeV^{2}");
/// K+ pi0
TH1 * h_pi0_resolved_candidates_M = new TH1D("h_pi0_resolved_candidates_M" , "#pi^{0} resolved mass spectra", 70, 100, 170);
h_pi0_resolved_candidates_M->GetXaxis()->SetTitle("#pi^{0} resolved mass [MeV]");
h_pi0_resolved_candidates_M->GetYaxis()->SetTitle("Candidates/1MeV");
TH1 * h_pi0_merged_candidates_M = new TH1D("h_pi0_merged_candidates_M" , "#pi^{0} merged mass spectra", 70, 100, 170);
h_pi0_merged_candidates_M->GetXaxis()->SetTitle("#pi^{0} merged mass [MeV]");
h_pi0_merged_candidates_M->GetYaxis()->SetTitle("Candidates/1MeV");
TH1 * h_pi0_merged_candidates_pT = new TH1D("h_pi0_merged_candidates_pT" , "p_{T} spectra of merged pion", 100, 0, 20000);
h_pi0_merged_candidates_pT->GetXaxis()->SetTitle("p_{T} (#pi^{0} merged) [MeV]");
h_pi0_merged_candidates_pT->GetYaxis()->SetTitle("Candidates/200MeV");
TH1 * h_pi0_resolved_candidates_pT = new TH1D("h_pi0_resolved_candidates_pT" , "p_{T} spectra of resolved pion", 100, 0, 20000);
h_pi0_resolved_candidates_pT->GetXaxis()->SetTitle("p_{T} (#pi^{0} resolved) [MeV]");
h_pi0_resolved_candidates_pT->GetYaxis()->SetTitle("Candidates/200MeV");
TH1 * h_KGammaGammaMuMu_mass_before_preselection = new TH1D("h_KGammaGammaMuMu_mass_before_preselection" , "m_{K#gamma#gamma#mu#mu before}", 100, get_cut_B_plus_M_low(year), cut_B_plus_M_high);
h_KGammaGammaMuMu_mass_before_preselection->GetXaxis()->SetTitle("m_{K#gamma#gamma#mu#mu} [MeV] before");
h_KGammaGammaMuMu_mass_before_preselection->GetYaxis()->SetTitle(Form("Counts / %.2fMeV", TMath::Abs(cut_B_plus_M_high - get_cut_B_plus_M_low(year))/100));
TH1 * h_KGammaGammaMuMu_mass_after_preselection= new TH1D("h_KGammaGammaMuMu_mass_after_preselection" , "m_{K#gamma#gamma#mu#mu after}", 100, get_cut_B_plus_M_low(year), cut_B_plus_M_high);
h_KGammaGammaMuMu_mass_after_preselection->GetXaxis()->SetTitle("m_{K#gamma#gamma#mu#mu} [MeV] after");
h_KGammaGammaMuMu_mass_after_preselection->GetYaxis()->SetTitle(Form("Counts / %.2fMeV", TMath::Abs(cut_B_plus_M_high - get_cut_B_plus_M_low(year))/100));
TH1 * h_KPiMuMu_KGammaGammaMuMu_mass= new TH1D("h_KPiMuMu_KGammaGammaMuMu_mass" , "m_{K#pi_{0}#mu#mu}", 100, get_cut_B_plus_M_low(year), cut_B_plus_M_high);
h_KPiMuMu_KGammaGammaMuMu_mass->GetXaxis()->SetTitle("m_{K#pi_{0}#mu#mu}[MeV]");
h_KPiMuMu_KGammaGammaMuMu_mass->GetYaxis()->SetTitle(Form("Counts / %.2fMeV", TMath::Abs(cut_B_plus_M_high - get_cut_B_plus_M_low(year))/100));
TH2 * h_pi0Mass2D = new TH2D("h_pi0Mass2D" , "mass spectrum of #pi^{0} ", 70, 100, 170, 70, 100, 170);
h_pi0Mass2D->GetXaxis()->SetTitle("mass from tree [MeV]");
h_pi0Mass2D->GetYaxis()->SetTitle("mass from #gamma #gamma [MeV]");
//-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//Define binning in DiMuon mass^2:
//bin# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
float f_BinBoundaries[19] = {0.1, 0.98, 1.1, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 11.0, 11.75, 12.5, 15.0, 16.0, 17.0, 18.0, 19.0, 20.0, 21.0}; //in (GeV/c2)^2
for(int i = 0; i < 19; i++)f_BinBoundaries[i] *= 1e6; //in (MeV/c2)^2
//phi -resonance: bin #2
//jpsi -resonance: bin #9
//psi2s-resonance: bin #12
unsigned int nDiMuonMassBin = 0;
summary_tree->Branch("nDiMuonMassBin", &nDiMuonMassBin, "nDiMuonMassBin/I");
Double_t dQ2;
summary_tree->Branch("Q2", &dQ2, "dQ2/D");
unsigned int Is_Phi_JPsi_Psi2S = 0;
summary_tree->Branch("Is_Phi_JPsi_Psi2S", &Is_Phi_JPsi_Psi2S, "Is_Phi_JPsi_Psi2S/I");
//generate subsets for BDT training, etc. corresponding to the eventnumber
//modulus(3) of the eventnumber
summary_tree->Branch("nSubSetNumber", &nSubSetNumber, "nSubSetNumber/I");
//check, whether the Kshort Decays within the velo
//this requires both Kshort daughter pions to be tracktype 3
unsigned int nKshortDecayInVeLo = 0;
if(Kst2Kspiplus){
summary_tree->Branch("KshortDecayInVeLo", &nKshortDecayInVeLo, "nKshortDecayInVeLo/I");
}
//save invariant masses of certain combinations into TBranches
Double_t dDiMuonM;
summary_tree->Branch("DiMuonInvMass", &dDiMuonM, "dDiMuonM/D");
Double_t dKshortMuMuM;
if(Kst2Kspiplus)summary_tree->Branch("KshortMuonsInvMass", &dKshortMuMuM, "dKshortMuMuM/D");
Double_t dKaonPionM;
summary_tree->Branch("KaonPionInvMass", &dKaonPionM, "dKaonPionM/D");
Double_t dReconstructedBM;
summary_tree->Branch("BplusInvMass", &dReconstructedBM, "dReconstructedBM/D");
Double_t dJpsiMisIDed_M;
if(Kst2Kspiplus)summary_tree->Branch("JpsiMisIDedM", &dJpsiMisIDed_M, "JpsiMisIDedM/D");
Double_t dReconstructedBMFixPiZero = 0;
if (Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){ //Add B+ mass with fixed pi0 mass
summary_tree->Branch("BplusInvMass_fixPiZero", &dReconstructedBMFixPiZero, "dReconstructedBMFixPiZero/D");
}
//Branch with gamma gamma angle
Double_t open_angle = 0.;
if (Kst2Kpluspi0Resolved){
summary_tree->Branch("open_angle", &open_angle, "open_angle/D");
}
//get same name branches with pTasymmetry in cones
//for the muDST of the data and Run 2 MC, the conept is gotton from the RELINFO tupletools
//while for Run 1 MC (which is full DST) the conept is determined properly by the default tool and created for multiple cone angles
Double_t B_plus_ConeDeltaETA, B_plus_ConeP, B_plus_ConePasym, B_plus_ConePT, B_plus_ConePTasym;
//create new branches for all trees with the 'new' cone variables! this is needed for the TMVA such that they have the same name
summary_tree->Branch("B_plus_NEW_ConeDeltaETA", &B_plus_ConeDeltaETA, "B_plus_ConeDeltaETA/D");
summary_tree->Branch("B_plus_NEW_ConeP", &B_plus_ConeP, "B_plus_ConeP/D");
summary_tree->Branch("B_plus_NEW_ConePasym", &B_plus_ConePasym, "B_plus_ConePasym/D");
summary_tree->Branch("B_plus_NEW_ConePT", &B_plus_ConePT, "B_plus_ConePT/D");
summary_tree->Branch("B_plus_NEW_ConePTasym", &B_plus_ConePTasym, "B_plus_ConePTasym/D");
//Rename DTF branches and save them
summary_tree->Branch("B_plus_M_DTF", &B_plus_M_DTF, "B_plus_M_DTF/D");
summary_tree->Branch("B_plus_MERR_DTF", &B_plus_MERR_DTF, "B_plus_MERR_DTF/D");
summary_tree->Branch("B_plus_P_DTF", &B_plus_P_DTF, "B_plus_P_DTF/D");
summary_tree->Branch("B_plus_PX_DTF", &B_plus_PX_DTF, "B_plus_PX_DTF/D");
summary_tree->Branch("B_plus_PY_DTF", &B_plus_PY_DTF, "B_plus_PY_DTF/D");
summary_tree->Branch("B_plus_PZ_DTF", &B_plus_PZ_DTF, "B_plus_PZ_DTF/D");
summary_tree->Branch("B_plus_PE_DTF", &B_plus_PE_DTF, "B_plus_PE_DTF/D");
summary_tree->Branch("B_plus_PT_DTF", &B_plus_PT_DTF, "B_plus_PT_DTF/D");
summary_tree->Branch("B_plus_ETA_DTF", &B_plus_ETA_DTF, "B_plus_ETA_DTF/D");
summary_tree->Branch("B_plus_PHI_DTF", &B_plus_PHI_DTF, "B_plus_PHI_DTF/D");
summary_tree->Branch("K_star_plus_M_DTF", &K_star_plus_M_DTF, "K_star_plus_M_DTF/D");
summary_tree->Branch("K_star_plus_P_DTF", &K_star_plus_P_DTF, "K_star_plus_P_DTF/D");
summary_tree->Branch("K_star_plus_PX_DTF", &K_star_plus_PX_DTF, "K_star_plus_PX_DTF/D");
summary_tree->Branch("K_star_plus_PY_DTF", &K_star_plus_PY_DTF, "K_star_plus_PY_DTF/D");
summary_tree->Branch("K_star_plus_PZ_DTF", &K_star_plus_PZ_DTF, "K_star_plus_PZ_DTF/D");
summary_tree->Branch("K_star_plus_PE_DTF", &K_star_plus_PE_DTF, "K_star_plus_PE_DTF/D");
summary_tree->Branch("K_star_plus_PT_DTF", &K_star_plus_PT_DTF, "K_star_plus_PT_DTF/D");
summary_tree->Branch("K_star_plus_ETA_DTF", &K_star_plus_ETA_DTF,"K_star_plus_ETA_DTF/D");
summary_tree->Branch("K_star_plus_PHI_DTF", &K_star_plus_PHI_DTF,"K_star_plus_PHI_DTF/D");
summary_tree->Branch("mu_plus_M_DTF", &mu_plus_M_DTF, "mu_plus_M_DTF/D");
summary_tree->Branch("mu_plus_P_DTF", &mu_plus_P_DTF, "mu_plus_P_DTF/D");
summary_tree->Branch("mu_plus_PX_DTF", &mu_plus_PX_DTF, "mu_plus_PX_DTF/D");
summary_tree->Branch("mu_plus_PY_DTF", &mu_plus_PY_DTF, "mu_plus_PY_DTF/D");
summary_tree->Branch("mu_plus_PZ_DTF", &mu_plus_PZ_DTF, "mu_plus_PZ_DTF/D");
summary_tree->Branch("mu_plus_PE_DTF", &mu_plus_PE_DTF, "mu_plus_PE_DTF/D");
summary_tree->Branch("mu_plus_PT_DTF", &mu_plus_PT_DTF, "mu_plus_PT_DTF/D");
summary_tree->Branch("mu_plus_ETA_DTF", &mu_plus_ETA_DTF, "mu_plus_ETA_DTF/D");
summary_tree->Branch("mu_plus_PHI_DTF", &mu_plus_PHI_DTF, "mu_plus_PHI_DTF/D");
summary_tree->Branch("mu_plus_ID_DTF", &mu_plus_ID_DTF, "mu_plus_ID_DTF/I");
summary_tree->Branch("mu_minus_M_DTF", &mu_minus_M_DTF, "mu_minus_M_DTF/D");
summary_tree->Branch("mu_minus_P_DTF", &mu_minus_P_DTF, "mu_minus_P_DTF/D");
summary_tree->Branch("mu_minus_PX_DTF", &mu_minus_PX_DTF, "mu_minus_PX_DTF/D");
summary_tree->Branch("mu_minus_PY_DTF", &mu_minus_PY_DTF, "mu_minus_PY_DTF/D");
summary_tree->Branch("mu_minus_PZ_DTF", &mu_minus_PZ_DTF, "mu_minus_PZ_DTF/D");
summary_tree->Branch("mu_minus_PE_DTF", &mu_minus_PE_DTF, "mu_minus_PE_DTF/D");
summary_tree->Branch("mu_minus_PT_DTF", &mu_minus_PT_DTF, "mu_minus_PT_DTF/D");
summary_tree->Branch("mu_minus_ETA_DTF", &mu_minus_ETA_DTF, "mu_minus_ETA_DTF/D");
summary_tree->Branch("mu_minus_PHI_DTF", &mu_minus_PHI_DTF, "mu_minus_PHI_DTF/D");
summary_tree->Branch("mu_minus_ID_DTF", &mu_minus_ID_DTF, "mu_minus_ID_DTF/I");
if (Kst2Kpluspi0Resolved){
summary_tree->Branch("pi_zero_resolved_M_DTF", &pi_zero_resolved_M_DTF, "pi_zero_resolved_M_DTF/D");
summary_tree->Branch("pi_zero_resolved_P_DTF", &pi_zero_resolved_P_DTF, "pi_zero_resolved_P_DTF/D");
summary_tree->Branch("pi_zero_resolved_PX_DTF", &pi_zero_resolved_PX_DTF, "pi_zero_resolved_PX_DTF/D");
summary_tree->Branch("pi_zero_resolved_PY_DTF", &pi_zero_resolved_PY_DTF, "pi_zero_resolved_PY_DTF/D");
summary_tree->Branch("pi_zero_resolved_PZ_DTF", &pi_zero_resolved_PZ_DTF, "pi_zero_resolved_PZ_DTF/D");
summary_tree->Branch("pi_zero_resolved_PE_DTF", &pi_zero_resolved_PE_DTF, "pi_zero_resolved_PE_DTF/D");
summary_tree->Branch("pi_zero_resolved_PT_DTF", &pi_zero_resolved_PT_DTF, "pi_zero_resolved_PT_DTF/D");
summary_tree->Branch("pi_zero_resolved_ETA_DTF", &pi_zero_resolved_ETA_DTF, "pi_zero_resolved_ETA_DTF/D");
summary_tree->Branch("pi_zero_resolved_PHI_DTF", &pi_zero_resolved_PHI_DTF, "pi_zero_resolved_PHI_DTF/D");
summary_tree->Branch("K_plus_M_DTF", &K_plus_M_DTF, "K_plus_M_DTF/D");
summary_tree->Branch("K_plus_P_DTF", &K_plus_P_DTF, "K_plus_P_DTF/D");
summary_tree->Branch("K_plus_PX_DTF", &K_plus_PX_DTF, "K_plus_PX_DTF/D");
summary_tree->Branch("K_plus_PY_DTF", &K_plus_PY_DTF, "K_plus_PY_DTF/D");
summary_tree->Branch("K_plus_PZ_DTF", &K_plus_PZ_DTF, "K_plus_PZ_DTF/D");
summary_tree->Branch("K_plus_PE_DTF", &K_plus_PE_DTF, "K_plus_PE_DTF/D");
summary_tree->Branch("K_plus_PT_DTF", &K_plus_PT_DTF, "K_plus_PT_DTF/D");
summary_tree->Branch("K_plus_ETA_DTF", &K_plus_ETA_DTF, "K_plus_ETA_DTF/D");
summary_tree->Branch("K_plus_PHI_DTF", &K_plus_PHI_DTF, "K_plus_PHI_DTF/D");
summary_tree->Branch("K_plus_ID_DTF", &K_plus_ID_DTF, "K_plus_ID_DTF/I");
summary_tree->Branch("gamma1_M_DTF", &gamma1_M_DTF, "gamma1_M_DTF/D");
summary_tree->Branch("gamma1_P_DTF", &gamma1_P_DTF, "gamma1_P_DTF/D");
summary_tree->Branch("gamma1_PX_DTF", &gamma1_PX_DTF, "gamma1_PX_DTF/D");
summary_tree->Branch("gamma1_PY_DTF", &gamma1_PY_DTF, "gamma1_PY_DTF/D");
summary_tree->Branch("gamma1_PZ_DTF", &gamma1_PZ_DTF, "gamma1_PZ_DTF/D");
summary_tree->Branch("gamma1_PE_DTF", &gamma1_PE_DTF, "gamma1_PE_DTF/D");
summary_tree->Branch("gamma1_PT_DTF", &gamma1_PT_DTF, "gamma1_PT_DTF/D");
summary_tree->Branch("gamma1_ETA_DTF", &gamma1_ETA_DTF, "gamma1_ETA_DTF/D");
summary_tree->Branch("gamma1_PHI_DTF", &gamma1_PHI_DTF, "gamma1_PHI_DTF/D");
summary_tree->Branch("gamma1_ID_DTF", &gamma1_ID_DTF, "gamma1_ID_DTF/I");
summary_tree->Branch("gamma2_M_DTF", &gamma2_M_DTF, "gamma2_M_DTF/D");
summary_tree->Branch("gamma2_P_DTF", &gamma2_P_DTF, "gamma2_P_DTF/D");
summary_tree->Branch("gamma2_PX_DTF", &gamma2_PX_DTF, "gamma2_PX_DTF/D");
summary_tree->Branch("gamma2_PY_DTF", &gamma2_PY_DTF, "gamma2_PY_DTF/D");
summary_tree->Branch("gamma2_PZ_DTF", &gamma2_PZ_DTF, "gamma2_PZ_DTF/D");
summary_tree->Branch("gamma2_PE_DTF", &gamma2_PE_DTF, "gamma2_PE_DTF/D");
summary_tree->Branch("gamma2_PT_DTF", &gamma2_PT_DTF, "gamma2_PT_DTF/D");
summary_tree->Branch("gamma2_ETA_DTF", &gamma2_ETA_DTF, "gamma2_ETA_DTF/D");
summary_tree->Branch("gamma2_PHI_DTF", &gamma2_PHI_DTF, "gamma2_PHI_DTF/D");
summary_tree->Branch("gamma2_ID_DTF", &gamma2_ID_DTF, "gamma2_ID_DTF/I");
}
if(Kst2Kspiplus){
summary_tree->Branch("K_short_M_DTF", &K_short_M_DTF, "K_short_M_DTF/D");
summary_tree->Branch("K_short_PT_DTF", &K_short_PT_DTF, "K_short_PT_DTF/D");
summary_tree->Branch("K_short_ETA_DTF", &K_short_ETA_DTF, "K_short_ETA_DTF/D");
summary_tree->Branch("K_short_PHI_DTF", &K_short_PHI_DTF, "K_short_PHI_DTF/D");
summary_tree->Branch("pi_plus_PX_DTF", &pi_plus_PX_DTF, "pi_plus_PX_DTF/D");
summary_tree->Branch("pi_plus_PY_DTF", &pi_plus_PY_DTF, "pi_plus_PY_DTF/D");
summary_tree->Branch("pi_plus_PZ_DTF", &pi_plus_PZ_DTF, "pi_plus_PZ_DTF/D");
summary_tree->Branch("pi_plus_PT_DTF", &pi_plus_PT_DTF, "pi_plus_PT_DTF/D");
summary_tree->Branch("pi_plus_ETA_DTF", &pi_plus_ETA_DTF, "pi_plus_ETA_DTF/D");
summary_tree->Branch("pi_plus_PHI_DTF", &pi_plus_PHI_DTF, "pi_plus_PHI_DTF/D");
summary_tree->Branch("pi_plus_PE_DTF", &pi_plus_PE_DTF, "pi_plus_PE_DTF/D");
summary_tree->Branch("pi_plus_ID_DTF", &pi_plus_ID_DTF, "pi_plus_ID_DTF/I");
summary_tree->Branch("Ks_pi_plus_PX_DTF", &Ks_pi_plus_PX_DTF, "Ks_pi_plus_PX_DTF/D");
summary_tree->Branch("Ks_pi_plus_PY_DTF", &Ks_pi_plus_PY_DTF, "Ks_pi_plus_PY_DTF/D");
summary_tree->Branch("Ks_pi_plus_PZ_DTF", &Ks_pi_plus_PZ_DTF, "Ks_pi_plus_PZ_DTF/D");
summary_tree->Branch("Ks_pi_plus_PT_DTF", &Ks_pi_plus_PT_DTF, "Ks_pi_plus_PT_DTF/D");
summary_tree->Branch("Ks_pi_plus_ETA_DTF", &Ks_pi_plus_ETA_DTF, "Ks_pi_plus_ETA_DTF/D");
summary_tree->Branch("Ks_pi_plus_PHI_DTF", &Ks_pi_plus_PHI_DTF, "Ks_pi_plus_PHI_DTF/D");
summary_tree->Branch("Ks_pi_plus_PE_DTF", &Ks_pi_plus_PE_DTF, "Ks_pi_plus_PE_DTF/D");
summary_tree->Branch("Ks_pi_plus_ID_DTF", &Ks_pi_plus_ID_DTF, "Ks_pi_plus_ID_DTF/I");
summary_tree->Branch("Ks_pi_minus_PX_DTF", &Ks_pi_minus_PX_DTF, "Ks_pi_minus_PX_DTF/D");
summary_tree->Branch("Ks_pi_minus_PY_DTF", &Ks_pi_minus_PY_DTF, "Ks_pi_minus_PY_DTF/D");
summary_tree->Branch("Ks_pi_minus_PZ_DTF", &Ks_pi_minus_PZ_DTF, "Ks_pi_minus_PZ_DTF/D");
summary_tree->Branch("Ks_pi_minus_PT_DTF", &Ks_pi_minus_PT_DTF, "Ks_pi_minus_PT_DTF/D");
summary_tree->Branch("Ks_pi_minus_ETA_DTF", &Ks_pi_minus_ETA_DTF, "Ks_pi_minus_ETA_DTF/D");
summary_tree->Branch("Ks_pi_minus_PHI_DTF", &Ks_pi_minus_PHI_DTF, "Ks_pi_minus_PHI_DTF/D");
summary_tree->Branch("Ks_pi_minus_PE_DTF", &Ks_pi_minus_PE_DTF, "Ks_pi_minus_PE_DTF/D");
summary_tree->Branch("Ks_pi_minus_ID_DTF", &Ks_pi_minus_ID_DTF, "Ks_pi_minus_ID_DTF/I");
}
//link variables according to MC or data for all cone variables
if (IsDST(year,MC,ReferenceChannel, PHSP) || (MC && B0 && ReferenceChannel) ){
tree -> SetBranchAddress( "B_plus_deltaEta_1.00" , &B_plus_ConeDeltaETA );
tree -> SetBranchAddress( "B_plus_cp_1.00" , &B_plus_ConeP );
tree -> SetBranchAddress( "B_plus_pasy_1.00" , &B_plus_ConePasym );
tree -> SetBranchAddress( "B_plus_cpt_1.00" , &B_plus_ConePT );
tree -> SetBranchAddress( "B_plus_ptasy_1.00" , &B_plus_ConePTasym );
}
else{
tree -> SetBranchAddress( "B_plus_CONEDELTAETA" , &B_plus_ConeDeltaETA );
tree -> SetBranchAddress( "B_plus_CONEP" , &B_plus_ConeP );
tree -> SetBranchAddress( "B_plus_CONEPASYM" , &B_plus_ConePasym );
tree -> SetBranchAddress( "B_plus_CONEPT" , &B_plus_ConePT );
tree -> SetBranchAddress( "B_plus_CONEPTASYM" , &B_plus_ConePTasym );
}
bool IsBplus = true;
///////////////////////////////////////////////////////////////////////////
///
/// loop over events
///
///////////////////////////////////////////////////////////////////////////
for(int i=0; i< N; i++){
//output the different counters!
if ((0ul == (i % 10000ul) || i + 1 == N) && i != 0) {
cout << "Read event " << i << "/" << N
<< ": tTRG true: " << nCounterTrigger << "/" << i //events passed trigger cut
<< " In K*+ mass: " << nCounterKstarMassWindow << "/" << i //events passed K* mass cut
<< " kinematics: " << nCounterKinematic << "/" << i //events passed kinematics cut
<< " PID: " << nCounterPID << "/" << i //events passed PID cuts
<< " Q2 Range: " << nCounterDiMuonMassInQ2Range << "/" << i //events passed Q^2 cut
<< " not Jpsi: " << nCounterJpsiResonance << "/" << i //events passed trigger cut
<< " B02Ksmm: " << nCounterB0KsmumuMass << "/" << nCounterDiMuonMassInQ2Range //events passed trigger cut
<< endl;
}
tree->GetEntry(i);
//If use DTF, check for convergence and kill the loop at an early stage
if (UseDTF && B_plus_DTF_status[0] > cut_DTF_status) continue;
if (UseDTF && B_plus_DTF_chi2[0] > cut_DTF_chi2) continue; //Sanity check on chi2
//use variable to use consistent pion and muon daugthers for the complex conjugated cases (B-)
//Correlation between B+ and muons/pions:
//B_plus_ID mu_minus_ID mu_plus_ID KS_piplus_ID KS_piminus_ID
//+521 +13 -13 +211 -211
//-521 -13 +13 -211 +211
//in DTF only final state particles are assign an ID, therefore, use the charged daughter of K*+ (i.e. K+ or pi+) for identification
//B+ ID > 0 <==> K+/pi+ ID > 0
if(Kst2Kspiplus)IsBplus = f_pi_plus_ID_DTF[0] > 0;
else IsBplus = f_K_plus_ID_DTF[0] > 0;
//convert Float_t (arrays) to Double_t
//momenta
if(Kst2Kspiplus){
pi_plus_PX_DTF = (Double_t) f_pi_plus_PX_DTF[0];
pi_plus_PY_DTF = (Double_t) f_pi_plus_PY_DTF[0];
pi_plus_PZ_DTF = (Double_t) f_pi_plus_PZ_DTF[0];
pi_plus_PE_DTF = (Double_t) f_pi_plus_PE_DTF[0];
pi_plus_ID_DTF = std::round(f_pi_plus_ID_DTF[0]);
if(IsBplus){
Ks_pi_minus_PX_DTF = (Double_t) f_Ks_pi_plus_PX_DTF[0];
Ks_pi_minus_PY_DTF = (Double_t) f_Ks_pi_plus_PY_DTF[0];
Ks_pi_minus_PZ_DTF = (Double_t) f_Ks_pi_plus_PZ_DTF[0];
Ks_pi_minus_PE_DTF = (Double_t) f_Ks_pi_plus_PE_DTF[0];
Ks_pi_minus_ID_DTF = std::round(f_Ks_pi_plus_ID_DTF[0]);
Ks_pi_plus_PX_DTF = (Double_t) f_Ks_pi_minus_PX_DTF[0];
Ks_pi_plus_PY_DTF = (Double_t) f_Ks_pi_minus_PY_DTF[0];
Ks_pi_plus_PZ_DTF = (Double_t) f_Ks_pi_minus_PZ_DTF[0];
Ks_pi_plus_PE_DTF = (Double_t) f_Ks_pi_minus_PE_DTF[0];
Ks_pi_plus_ID_DTF = std::round(f_Ks_pi_minus_ID_DTF[0]);
}
else{
Ks_pi_plus_PX_DTF = (Double_t) f_Ks_pi_plus_PX_DTF[0];
Ks_pi_plus_PY_DTF = (Double_t) f_Ks_pi_plus_PY_DTF[0];
Ks_pi_plus_PZ_DTF = (Double_t) f_Ks_pi_plus_PZ_DTF[0];
Ks_pi_plus_PE_DTF = (Double_t) f_Ks_pi_plus_PE_DTF[0];
Ks_pi_plus_ID_DTF = std::round(f_Ks_pi_plus_ID_DTF[0]);
Ks_pi_minus_PX_DTF = (Double_t) f_Ks_pi_minus_PX_DTF[0];
Ks_pi_minus_PY_DTF = (Double_t) f_Ks_pi_minus_PY_DTF[0];
Ks_pi_minus_PZ_DTF = (Double_t) f_Ks_pi_minus_PZ_DTF[0];
Ks_pi_minus_PE_DTF = (Double_t) f_Ks_pi_minus_PE_DTF[0];
Ks_pi_minus_ID_DTF = std::round(f_Ks_pi_minus_ID_DTF[0]);
}
}
if(Kst2Kpluspi0Resolved){
K_plus_PX_DTF = (Double_t) f_K_plus_PX_DTF[0];
K_plus_PY_DTF = (Double_t) f_K_plus_PY_DTF[0];
K_plus_PZ_DTF = (Double_t) f_K_plus_PZ_DTF[0];
K_plus_PE_DTF = (Double_t) f_K_plus_PE_DTF[0];
K_plus_ID_DTF = std::round(f_K_plus_ID_DTF[0]);
gamma1_PX_DTF = (Double_t) f_gamma1_PX_DTF[0];
gamma1_PY_DTF = (Double_t) f_gamma1_PY_DTF[0];
gamma1_PZ_DTF = (Double_t) f_gamma1_PZ_DTF[0];
gamma1_PE_DTF = (Double_t) f_gamma1_PE_DTF[0];
gamma1_ID_DTF = std::round(f_gamma1_ID_DTF[0]);
gamma2_PX_DTF = (Double_t) f_gamma2_PX_DTF[0];
gamma2_PY_DTF = (Double_t) f_gamma2_PY_DTF[0];
gamma2_PZ_DTF = (Double_t) f_gamma2_PZ_DTF[0];
gamma2_PE_DTF = (Double_t) f_gamma2_PE_DTF[0];
gamma2_ID_DTF = std::round(f_gamma2_ID_DTF[0]);
}
if(IsBplus){
mu_minus_PX_DTF = (Double_t) f_mu_plus_PX_DTF[0];
mu_minus_PY_DTF = (Double_t) f_mu_plus_PY_DTF[0];
mu_minus_PZ_DTF = (Double_t) f_mu_plus_PZ_DTF[0];
mu_minus_PE_DTF = (Double_t) f_mu_plus_PE_DTF[0];
mu_minus_ID_DTF = std::round(f_mu_plus_ID_DTF[0]);
mu_plus_PX_DTF = (Double_t) f_mu_minus_PX_DTF[0];
mu_plus_PY_DTF = (Double_t) f_mu_minus_PY_DTF[0];
mu_plus_PZ_DTF = (Double_t) f_mu_minus_PZ_DTF[0];
mu_plus_PE_DTF = (Double_t) f_mu_minus_PE_DTF[0];
mu_plus_ID_DTF = std::round(f_mu_minus_ID_DTF[0]);
}
else{
mu_plus_PX_DTF = (Double_t) f_mu_plus_PX_DTF[0];
mu_plus_PY_DTF = (Double_t) f_mu_plus_PY_DTF[0];
mu_plus_PZ_DTF = (Double_t) f_mu_plus_PZ_DTF[0];
mu_plus_PE_DTF = (Double_t) f_mu_plus_PE_DTF[0];
mu_plus_ID_DTF = std::round(f_mu_plus_ID_DTF[0]);
mu_minus_PX_DTF = (Double_t) f_mu_minus_PX_DTF[0];
mu_minus_PY_DTF = (Double_t) f_mu_minus_PY_DTF[0];
mu_minus_PZ_DTF = (Double_t) f_mu_minus_PZ_DTF[0];
mu_minus_PE_DTF = (Double_t) f_mu_minus_PE_DTF[0];
mu_minus_ID_DTF = std::round(f_mu_minus_ID_DTF[0]);
}
//convert masses Float_t to Double_t
B_plus_M_DTF = (Double_t) f_B_plus_M_DTF[0];
B_plus_MERR_DTF = (Double_t) f_B_plus_MERR_DTF[0];
K_star_plus_M_DTF = (Double_t) f_K_star_plus_M_DTF[0];
if(Kst2Kspiplus){
K_short_M_DTF = (Double_t) f_K_short_M_DTF[0];
}
if(Kst2Kpluspi0Resolved){
pi_zero_resolved_M_DTF = (Double_t) f_pi_zero_resolved_M_DTF[0];
}
//define the Lorentz vectors
if(UseDTF){ //start useDTF for LorentzVectors
if(UsePDGmIn4Dvectors){
LorVec_mu_minus.SetXYZM (mu_minus_PX_DTF, mu_minus_PY_DTF, mu_minus_PZ_DTF, PDGMASS.MU);
LorVec_mu_plus.SetXYZM (mu_plus_PX_DTF, mu_plus_PY_DTF, mu_plus_PZ_DTF, PDGMASS.MU);
if(Kst2Kspiplus){
/// Ks, pi+
LorVec_pi_plus.SetXYZM (pi_plus_PX_DTF, pi_plus_PY_DTF, pi_plus_PZ_DTF, PDGMASS.PI_PLUS);
LorVec_KS_pi_plus.SetXYZM (Ks_pi_plus_PX_DTF, Ks_pi_plus_PY_DTF, Ks_pi_plus_PZ_DTF, PDGMASS.PI_PLUS);
LorVec_KS_pi_minus.SetXYZM (Ks_pi_minus_PX_DTF, Ks_pi_minus_PY_DTF, Ks_pi_minus_PZ_DTF, PDGMASS.PI_PLUS);
LorVec_pi_plus_MisIDed.SetXYZM (pi_plus_PX_DTF, pi_plus_PY_DTF, pi_plus_PZ_DTF, PDGMASS.MU);
LorVec_mu_plus_MisIDed.SetXYZM (mu_plus_PX_DTF, mu_plus_PY_DTF, mu_plus_PZ_DTF, PDGMASS.PI_PLUS);
}
else{
/// K+, pi0
LorVec_K_plus.SetXYZM (K_plus_PX_DTF, K_plus_PY_DTF, K_plus_PZ_DTF, PDGMASS.K_PLUS);
LorVec_gamma1.SetXYZM (gamma1_PX_DTF, gamma1_PY_DTF, gamma1_PZ_DTF, PDGMASS.GAMMA);
LorVec_gamma2.SetXYZM (gamma2_PX_DTF, gamma2_PY_DTF, gamma2_PZ_DTF, PDGMASS.GAMMA);
}
}
else{
LorVec_mu_minus.SetPxPyPzE (mu_minus_PX_DTF, mu_minus_PY_DTF, mu_minus_PZ_DTF, mu_minus_PE_DTF);
LorVec_mu_plus.SetPxPyPzE (mu_plus_PX_DTF, mu_plus_PY_DTF, mu_plus_PZ_DTF, mu_plus_PE_DTF);
if(Kst2Kspiplus){
/// Ks, pi+
LorVec_pi_plus.SetPxPyPzE (pi_plus_PX_DTF, pi_plus_PY_DTF, pi_plus_PZ_DTF, pi_plus_PE_DTF);
LorVec_KS_pi_plus.SetPxPyPzE (Ks_pi_plus_PX_DTF, Ks_pi_plus_PY_DTF, Ks_pi_plus_PZ_DTF, Ks_pi_plus_PE_DTF);
LorVec_KS_pi_minus.SetPxPyPzE (Ks_pi_minus_PX_DTF, Ks_pi_minus_PY_DTF, Ks_pi_minus_PZ_DTF, Ks_pi_minus_PE_DTF);
LorVec_pi_plus_MisIDed.SetXYZM (pi_plus_PX_DTF, pi_plus_PY_DTF, pi_plus_PZ_DTF, PDGMASS.MU);
LorVec_mu_plus_MisIDed.SetXYZM (mu_plus_PX_DTF, mu_plus_PY_DTF, mu_plus_PZ_DTF, PDGMASS.PI_PLUS);
}
else{
/// K+, pi0
LorVec_K_plus.SetPxPyPzE (K_plus_PX_DTF, K_plus_PY_DTF, K_plus_PZ_DTF, K_plus_PE_DTF);
LorVec_gamma1.SetPxPyPzE (gamma1_PX_DTF, gamma1_PY_DTF, gamma1_PZ_DTF, gamma1_PE_DTF);
LorVec_gamma2.SetPxPyPzE (gamma2_PX_DTF, gamma2_PY_DTF, gamma2_PZ_DTF, gamma2_PE_DTF);
}
}
//directly recombine the non-stable particles for later use (K*+ not really needed/used for Kst2Kspiplus)
if(Kst2Kspiplus){
LorVec_K_short = LorVec_KS_pi_plus + LorVec_KS_pi_minus;
LorVec_KaonPion = LorVec_K_short + LorVec_pi_plus;
}
if(Kst2Kpluspi0Resolved){
LorVec_DiGamma = LorVec_gamma1 + LorVec_gamma2;
LorVec_pi_zero_resolved = LorVec_DiGamma;
LorVec_KaonPion = LorVec_K_plus + LorVec_pi_zero_resolved;
}
LorVec_K_star_plus = LorVec_KaonPion;
LorVec_DiMuon = LorVec_mu_minus + LorVec_mu_plus;
LorVec_B_plus = LorVec_K_star_plus + LorVec_DiMuon;
//determine transvers momenta PT for all particles:
B_plus_P_DTF = LorVec_B_plus.P();
B_plus_PX_DTF = LorVec_B_plus.Px();
B_plus_PY_DTF = LorVec_B_plus.Py();
B_plus_PZ_DTF = LorVec_B_plus.Pz();
B_plus_PT_DTF = LorVec_B_plus.Pt();
B_plus_PE_DTF = LorVec_B_plus.E();
B_plus_ETA_DTF = LorVec_B_plus.Eta();
B_plus_PHI_DTF = LorVec_B_plus.Phi();
K_star_plus_M_DTF = LorVec_K_star_plus.M();
K_star_plus_P_DTF = LorVec_K_star_plus.P();
K_star_plus_PX_DTF = LorVec_K_star_plus.Px();
K_star_plus_PY_DTF = LorVec_K_star_plus.Py();
K_star_plus_PZ_DTF = LorVec_K_star_plus.Pz();
K_star_plus_PT_DTF = LorVec_K_star_plus.Pt();
K_star_plus_PE_DTF = LorVec_K_star_plus.E();
K_star_plus_ETA_DTF = LorVec_K_star_plus.Eta();
K_star_plus_PHI_DTF = LorVec_K_star_plus.Phi();
mu_plus_P_DTF = LorVec_mu_plus.P();
mu_plus_M_DTF = LorVec_mu_plus.M();
mu_plus_PT_DTF = LorVec_mu_plus.Pt();
mu_plus_ETA_DTF = LorVec_mu_plus.Eta();
mu_plus_PHI_DTF = LorVec_mu_plus.Phi();
mu_minus_P_DTF = LorVec_mu_minus.P();
mu_minus_M_DTF = LorVec_mu_minus.M();
mu_minus_PT_DTF = LorVec_mu_minus.Pt();
mu_minus_ETA_DTF = LorVec_mu_minus.Eta();
mu_minus_PHI_DTF = LorVec_mu_minus.Phi();
if(Kst2Kspiplus){
K_short_PT_DTF = LorVec_K_short.Pt();
K_short_ETA_DTF = LorVec_K_short.Eta();
K_short_PHI_DTF = LorVec_K_short.Phi();
Ks_pi_plus_PT_DTF = LorVec_KS_pi_plus.Pt();
Ks_pi_plus_ETA_DTF = LorVec_KS_pi_plus.Eta();
Ks_pi_plus_PHI_DTF = LorVec_KS_pi_plus.Phi();
Ks_pi_minus_PT_DTF = LorVec_KS_pi_minus.Pt();
Ks_pi_minus_ETA_DTF = LorVec_KS_pi_minus.Eta();
Ks_pi_minus_PHI_DTF = LorVec_KS_pi_minus.Phi();
pi_plus_PT_DTF = LorVec_pi_plus.Pt();
pi_plus_ETA_DTF = LorVec_pi_plus.Eta();
pi_plus_PHI_DTF = LorVec_pi_plus.Phi();
}
if(Kst2Kpluspi0Resolved){
K_plus_P_DTF = LorVec_K_plus.P();
K_plus_M_DTF = LorVec_K_plus.M();
K_plus_PT_DTF = LorVec_K_plus.Pt();
K_plus_ETA_DTF = LorVec_K_plus.Eta();
K_plus_PHI_DTF = LorVec_K_plus.Phi();
pi_zero_resolved_P_DTF = LorVec_pi_zero_resolved.P();
pi_zero_resolved_PX_DTF = LorVec_pi_zero_resolved.Px();
pi_zero_resolved_PY_DTF = LorVec_pi_zero_resolved.Py();
pi_zero_resolved_PZ_DTF = LorVec_pi_zero_resolved.Pz();
pi_zero_resolved_PT_DTF = LorVec_pi_zero_resolved.Pt();
pi_zero_resolved_PE_DTF = LorVec_pi_zero_resolved.E();
pi_zero_resolved_ETA_DTF= LorVec_pi_zero_resolved.Eta();
pi_zero_resolved_PHI_DTF= LorVec_pi_zero_resolved.Phi();
gamma1_M_DTF = LorVec_gamma1.M();
gamma1_P_DTF = LorVec_gamma1.P();
gamma1_PT_DTF = LorVec_gamma1.Pt();
gamma1_ETA_DTF = LorVec_gamma1.Eta();
gamma1_PHI_DTF = LorVec_gamma1.Phi();
gamma2_M_DTF = LorVec_gamma2.M();
gamma2_P_DTF = LorVec_gamma2.P();
gamma2_PT_DTF = LorVec_gamma2.Pt();
gamma2_ETA_DTF = LorVec_gamma2.Eta();
gamma2_PHI_DTF = LorVec_gamma2.Phi();
}
} //end useDTF for LorentzVectors
else{ //NON-DTF variables
if(UsePDGmIn4Dvectors){ //USE PDG masses
LorVec_B_plus.SetXYZM (B_plus_PX, B_plus_PY, B_plus_PZ, PDGMASS.B_PLUS);
LorVec_K_star_plus.SetXYZM (K_star_plus_PX, K_star_plus_PY, K_star_plus_PZ, PDGMASS.K_STAR_PLUS);
LorVec_mu_minus.SetXYZM (mu_minus_PX, mu_minus_PY, mu_minus_PZ, PDGMASS.MU);
LorVec_mu_plus.SetXYZM (mu_plus_PX, mu_plus_PY, mu_plus_PZ, PDGMASS.MU);
if(Kst2Kspiplus){
/// Ks, pi+
LorVec_K_short.SetXYZM (K_short_PX, K_short_PY, K_short_PZ, PDGMASS.K_SHORT);
LorVec_pi_plus.SetXYZM (pi_plus_PX, pi_plus_PY, pi_plus_PZ, PDGMASS.PI_PLUS);
LorVec_KS_pi_plus.SetXYZM (Ks_pi_plus_PX, Ks_pi_plus_PY, Ks_pi_plus_PZ, PDGMASS.PI_PLUS);
LorVec_KS_pi_minus.SetXYZM (Ks_pi_minus_PX, Ks_pi_minus_PY, Ks_pi_minus_PZ, PDGMASS.PI_PLUS);
LorVec_pi_plus_MisIDed.SetXYZM (pi_plus_PX, pi_plus_PY, pi_plus_PZ, PDGMASS.MU);
LorVec_mu_plus_MisIDed.SetXYZM (mu_plus_PX, mu_plus_PY, mu_plus_PZ, PDGMASS.PI_PLUS);
}
else{
/// K+, pi0
LorVec_K_plus.SetXYZM (K_plus_PX, K_plus_PY, K_plus_PZ, PDGMASS.K_PLUS);
LorVec_pi_zero_merged.SetXYZM (pi_zero_merged_PX, pi_zero_merged_PY, pi_zero_merged_PZ, PDGMASS.PI_ZERO);
LorVec_pi_zero_resolved.SetXYZM (pi_zero_resolved_PX, pi_zero_resolved_PY, pi_zero_resolved_PZ, PDGMASS.PI_ZERO);
LorVec_gamma1.SetXYZM (gamma1_PX, gamma1_PY, gamma1_PZ, PDGMASS.GAMMA);
LorVec_gamma2.SetXYZM (gamma2_PX, gamma2_PY, gamma2_PZ, PDGMASS.GAMMA);
}
}
else{
LorVec_B_plus.SetPxPyPzE (B_plus_PX, B_plus_PY, B_plus_PZ, B_plus_PE);
LorVec_K_star_plus.SetPxPyPzE (K_star_plus_PX, K_star_plus_PY, K_star_plus_PZ, K_star_plus_PE);
LorVec_mu_minus.SetPxPyPzE (mu_minus_PX, mu_minus_PY, mu_minus_PZ, mu_minus_PE);
LorVec_mu_plus.SetPxPyPzE (mu_plus_PX, mu_plus_PY, mu_plus_PZ, mu_plus_PE);
if(Kst2Kspiplus){
/// Ks, pi+
LorVec_K_short.SetPxPyPzE (K_short_PX, K_short_PY, K_short_PZ, K_short_PE);
LorVec_pi_plus.SetPxPyPzE (pi_plus_PX, pi_plus_PY, pi_plus_PZ, pi_plus_PE);
LorVec_KS_pi_plus.SetPxPyPzE (Ks_pi_plus_PX, Ks_pi_plus_PY, Ks_pi_plus_PZ, Ks_pi_plus_PE);
LorVec_KS_pi_minus.SetPxPyPzE (Ks_pi_minus_PX, Ks_pi_minus_PY, Ks_pi_minus_PZ, Ks_pi_minus_PE);
LorVec_pi_plus_MisIDed.SetXYZM (pi_plus_PX, pi_plus_PY, pi_plus_PZ, PDGMASS.MU);
LorVec_mu_plus_MisIDed.SetXYZM (mu_plus_PX, mu_plus_PY, mu_plus_PZ, PDGMASS.PI_PLUS);
}
else{
/// K+, pi0
LorVec_K_plus.SetPxPyPzE (K_plus_PX, K_plus_PY, K_plus_PZ, K_plus_PE);
LorVec_pi_zero_merged.SetPxPyPzE (pi_zero_merged_PX, pi_zero_merged_PY, pi_zero_merged_PZ, pi_zero_merged_PE);
LorVec_pi_zero_resolved.SetPxPyPzE (pi_zero_resolved_PX, pi_zero_resolved_PY, pi_zero_resolved_PZ, pi_zero_resolved_PE);
LorVec_gamma1.SetPxPyPzE (gamma1_PX, gamma1_PY, gamma1_PZ, gamma1_PE);
LorVec_gamma2.SetPxPyPzE (gamma2_PX, gamma2_PY, gamma2_PZ, gamma2_PE);
}
}
//Recombination of 'Parent'-Particles
LorVec_DiMuon = LorVec_mu_minus + LorVec_mu_plus;
if(Kst2Kspiplus){
LorVec_KaonPion = LorVec_K_short + LorVec_pi_plus;
}
if(Kst2Kpluspi0Merged){
LorVec_KaonPion = LorVec_K_plus + LorVec_pi_zero_merged;
}
if(Kst2Kpluspi0Resolved){
LorVec_DiGamma = LorVec_gamma1 + LorVec_gamma2;
LorVec_KaonPion = LorVec_K_plus + LorVec_DiGamma;
}
}
//Recombination of Particles
if(Kst2Kspiplus){
LorVec_KshortMuMu = LorVec_K_short + LorVec_DiMuon;
}
if(Kst2Kpluspi0Merged){
LorVec_KplusMuMu = LorVec_K_plus+ LorVec_DiMuon;
}
if(Kst2Kpluspi0Resolved){
LorVec_KplusMuMu = LorVec_K_plus+ LorVec_DiMuon;
LorVec_ReconBPiZero = LorVec_DiMuon + LorVec_K_plus + LorVec_pi_zero_resolved;
dReconstructedBMFixPiZero = LorVec_ReconBPiZero.M();
open_angle =LorVec_gamma1.Angle(LorVec_gamma2.Vect());
}
LorVec_K_star_plus = LorVec_KaonPion;
LorVec_ReconB = LorVec_DiMuon + LorVec_K_star_plus;
//double misIDed
if(Kst2Kspiplus){
LorVec_K_star_plus_From_DoubleMisIDed = LorVec_K_short + LorVec_mu_plus_MisIDed;
LorVec_Jpsi_From_DoubleMisIDed = LorVec_pi_plus_MisIDed + LorVec_mu_minus;
LorVec_B_plus_From_DoubleMisIDed = LorVec_Jpsi_From_DoubleMisIDed + LorVec_K_star_plus_From_DoubleMisIDed;
}
//save masses to doubles!
dDiMuonM = LorVec_DiMuon.M();
dKaonPionM = LorVec_KaonPion.M();
dReconstructedBM = LorVec_ReconB.M();
if(Kst2Kspiplus) dKshortMuMuM = LorVec_KshortMuMu.M();
//LL (Kshort decays inside the velo)
if(Kst2Kspiplus){
h_KshortDaughters->Fill(Ks_pi_minus_TRACK_Type, Ks_pi_plus_TRACK_Type);
if(Ks_pi_minus_TRACK_Type == 3 && Ks_pi_plus_TRACK_Type == 3){
nKshortDecayInVeLo = 1;
}
//DD (Kshort leaves velo)
else nKshortDecayInVeLo = 0;
if(nKshortDecayInVeLo)//cut on pion daughter PT for LL tracks!
if(Ks_pi_plus_PT < 250. || Ks_pi_minus_PT < 250.)
continue;
}
h_BplusMass_DiMuonMassRaw->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
if(Kst2Kspiplus){
if(nKshortDecayInVeLo == 0)h_BplusMass_DiMuonMassRawDD->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
else h_BplusMass_DiMuonMassRawLL->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
}
//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//trigger selection
//L0 stage
if(B_plus_L0MuonDecision_TOS == 0) continue;
//HLT 1 stage
if((B_plus_Hlt1TrackAllL0Decision_TOS == 0) && (B_plus_Hlt1TrackMuonDecision_TOS == 0)) continue;
//HLT2 stage
if((B_plus_Hlt2Topo2BodyBBDTDecision_TOS == 0) &&
(B_plus_Hlt2Topo3BodyBBDTDecision_TOS == 0) &&
(B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS == 0) &&
(B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS == 0) &&
//(B_plus_Hlt2SingleMuonDecision_TOS== 0) && //remove SingleMuon from selection, since its exclusive yield is 1% and pre-scaling messes up data/MC agreement
(B_plus_Hlt2DiMuonDetachedDecision_TOS== 0) &&
(B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS == 0) &&
(B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS == 0)
) continue;
nCounterTrigger++;
//---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//some mass plots
h_KstarMass2D->Fill(K_star_plus_M, LorVec_KaonPion.M());
h_KstarMass->Fill(LorVec_KaonPion.M());
if (Kst2Kpluspi0Merged){
h_KGammaGammaMuMu_mass_before_preselection->Fill(dReconstructedBM);
h_pi0_merged_candidates_pT->Fill(LorVec_pi_zero_merged.Perp());
h_pi0_merged_candidates_M->Fill(LorVec_pi_zero_merged.M());
}
if (Kst2Kpluspi0Resolved){
h_pi0_resolved_candidates_pT->Fill(LorVec_pi_zero_resolved.Perp());
h_KGammaGammaMuMu_mass_before_preselection->Fill(dReconstructedBM);
h_pi0_resolved_candidates_M->Fill(LorVec_DiGamma.M());
h_pi0Mass2D->Fill(pi_zero_resolved_M,LorVec_DiGamma.M());
}
h_KstarMass2D->Fill(K_star_plus_M, LorVec_KaonPion.M());
h_KstarMass->Fill(LorVec_KaonPion.M());
//mass cuts
if(TMath::Abs((UseDTF ? K_star_plus_M_DTF : K_star_plus_M) - PDGMASS.K_STAR_PLUS) > K_star_plus_mass_difference) continue; //remove non-resonant K+pi
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
if((UseDTF ? LorVec_K_star_plus.Pt() : K_star_plus_PT) < cut_K_star_plus_pT)continue;
//if(TMath::Abs(LorVec_KaonPion.M() - PDGMASS.K_STAR_PLUS) > 150) continue; //remove non-resonant K+pi
}
h_KstarMasscut->Fill(LorVec_KaonPion.M());
nCounterKstarMassWindow++;
if((UseDTF ? B_plus_M_DTF : B_plus_M) < get_cut_B_plus_M_low(year) || (UseDTF ? B_plus_M_DTF : B_plus_M) > cut_B_plus_M_high)continue;
//Check for non-DTF mass to allign with the older stripping
if (Kst2Kpluspi0Resolved){
if(B_plus_M < cut_B_plus_M_low_basic) continue;
}
//kinematic cuts
if(LorVec_B_plus.Theta() <= cut_kin_Theta_low || LorVec_B_plus.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_K_star_plus.Theta() <= cut_kin_Theta_low || LorVec_K_star_plus.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_mu_minus.Theta() <= cut_kin_Theta_low || LorVec_mu_minus.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_mu_plus.Theta() <= cut_kin_Theta_low || LorVec_mu_plus.Theta() >= cut_kin_Theta_up)continue;
/// newly added
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
if ((UseDTF ? LorVec_B_plus.Pt() : B_plus_PT) < cut_B_plus_pT) continue;
if (B_plus_DIRA_OWNPV < cut_B_plus_DIRA) continue;
}
if(Kst2Kspiplus){
if(LorVec_K_short.Theta() <= cut_kin_Theta_low || LorVec_K_short.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_pi_plus.Theta() <= cut_kin_Theta_low || LorVec_pi_plus.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_KS_pi_plus.Theta() <= cut_kin_Theta_low || LorVec_KS_pi_plus.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_KS_pi_minus.Theta() <= cut_kin_Theta_low || LorVec_KS_pi_minus.Theta() >= cut_kin_Theta_up)continue;
}
if(Kst2Kpluspi0Merged){
if(LorVec_K_plus.Theta() <= cut_kin_Theta_low || LorVec_K_plus.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_pi_zero_merged.Theta() <= cut_kin_Theta_low || LorVec_pi_zero_merged.Theta() >= cut_kin_Theta_up)continue;
if((UseDTF ? LorVec_K_plus.Pt() : K_plus_PT) < 1350) continue;
}
if(Kst2Kpluspi0Resolved){
if(LorVec_K_plus.Theta() <= cut_kin_Theta_low || LorVec_K_plus.Theta() >= cut_kin_Theta_up)continue;
if(LorVec_pi_zero_resolved.Theta() <= cut_kin_Theta_low || LorVec_pi_zero_resolved.Theta() >= cut_kin_Theta_up)continue;
}
if(LorVec_mu_minus.Angle(LorVec_mu_plus.Vect()) <= cut_muon_angle)continue;
if(Kst2Kspiplus){
// if(LorVec_K_short.Angle(LorVec_pi_plus.Vect()) <= cut_muon_angle)continue;
if(LorVec_pi_plus.Angle(LorVec_KS_pi_minus.Vect()) <= cut_muon_angle)continue;
if(LorVec_pi_plus.Angle(LorVec_KS_pi_plus.Vect()) <= cut_muon_angle)continue;
if(LorVec_KS_pi_plus.Angle(LorVec_KS_pi_minus.Vect()) <= cut_muon_angle)continue;
}
if(Kst2Kpluspi0Merged)if(LorVec_K_plus.Angle(LorVec_pi_zero_merged.Vect()) <= cut_muon_angle)continue;
if(Kst2Kpluspi0Resolved)if(LorVec_K_plus.Angle(LorVec_pi_zero_resolved.Vect()) <= cut_muon_angle)continue;
/// newly added
if(Kst2Kpluspi0Merged)if((UseDTF ? LorVec_pi_zero_merged.Pt() : pi_zero_merged_PT) < cut_pi_zero_PT)continue;
if(Kst2Kpluspi0Resolved)if((UseDTF ? LorVec_pi_zero_resolved.Pt() : pi_zero_resolved_PT) < cut_pi_zero_PT)continue;
/// newly added
//flight distance cuts
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
if (B_plus_FDCHI2_OWNPV < cut_B_plus_FDCHI2) continue;
if ((!Inc && !K1 && !B0 ) && K_star_plus_FDCHI2_OWNPV < cut_K_star_plus_FDCHI2) continue;
}
nCounterKinematic++;
//primary vertex cut
if (B_plus_IPCHI2_OWNPV < cut_IPCHI2_OWNPV_low)continue;
if (K_star_plus_IPCHI2_OWNPV < cut_IPCHI2_OWNPV_low)continue;
if (mu_plus_IPCHI2_OWNPV < cut_muon_IPCHI2_OWNPV_low)continue;
if (mu_minus_IPCHI2_OWNPV < cut_muon_IPCHI2_OWNPV_low)continue;
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged)if (B_plus_IPCHI2_OWNPV > cut_IPCHI2_OWNPV_high)continue;
nCounterPrV++;
/// newly added
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
if (B_plus_ConePTasym < cut_B_plus_ConePTasym)continue;
}
//position vertex cuts?
nCounterPoV++;
h_BeforePIDcuts_ctl->Fill(TMath::Cos(ctl));
h_BeforePIDcuts_ctk->Fill(TMath::Cos(ctk));
h_BeforePIDcuts_phi->Fill(phi);
//muon cuts
if(!mu_plus_isMuon || !mu_minus_isMuon)continue;
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
if (LorVec_DiMuon.M() > cut_DiMuon_M) continue; //Should be in stripping already
}
//particle ID
if(mu_minus_ProbNNmu <= cut_muon_ProbNNmu)continue;
if(mu_plus_ProbNNmu <= cut_muon_ProbNNmu)continue;
if(Kst2Kspiplus){
if(!pi_plus_hasRich)continue;
//if(pi_plus_PIDK >= 25)continue;
}
if(Kst2Kpluspi0Merged){
if(!K_plus_hasRich)continue;
if(K_plus_ProbNNk <= cut_K_plus_ProbNNk)continue;
}
if(Kst2Kpluspi0Resolved){
if(!K_plus_hasRich)continue;
if(K_plus_ProbNNk <= cut_K_plus_ProbNNk)continue;
if(gamma1_CL <= cut_gamma_CL)continue;
if(gamma2_CL <= cut_gamma_CL)continue;
}
h_AfterPIDcuts_ctl->Fill(TMath::Cos(ctl));
h_AfterPIDcuts_ctk->Fill(TMath::Cos(ctk));
h_AfterPIDcuts_phi->Fill(phi);
nCounterPID++;
h_JpsiMass->Fill(LorVec_DiMuon.M());
//only keep events inside q2 region
if(LorVec_DiMuon.M2() < f_BinBoundaries[0] || LorVec_DiMuon.M2() > f_BinBoundaries[18])continue;
nCounterDiMuonMassInQ2Range++;
h_diMuonReconstructedBMass2D->Fill(dReconstructedBM,dDiMuonM);
if(Kst2Kpluspi0Resolved || Kst2Kpluspi0Merged){
h_diMuonReconstructedBMass2DFixPiZero->Fill(dReconstructedBMFixPiZero,dDiMuonM);
}
if (Kst2Kpluspi0Merged){
h_KGammaGammaMuMu_mass_after_preselection->Fill(dReconstructedBM);
}
if (Kst2Kpluspi0Resolved){
h_KGammaGammaMuMu_mass_after_preselection->Fill(dReconstructedBM);
h_KPiMuMu_KGammaGammaMuMu_mass->Fill(dReconstructedBMFixPiZero);
}
h_BplusMass_B02KsmumuMass->Fill(LorVec_KshortMuMu.M());
if(Kst2Kspiplus){ //B0 veto
if(TMath::Abs(LorVec_KshortMuMu.M() - PDGMASS.B_ZERO) < 50){
//h_BplusMass_B02KsmumuMass->Fill(B_plus_M);
nCounterB0KsmumuMass++;
continue;
}
}
///////////////////////////////////////
///////////////////////////////////////
///////////////////////////////////////
nDiMuonMassBin = 0;
Is_Phi_JPsi_Psi2S = 0;
while(LorVec_DiMuon.M2() >= f_BinBoundaries[nDiMuonMassBin]){
nDiMuonMassBin++;
}
if(nDiMuonMassBin == 2){
Is_Phi_JPsi_Psi2S = 1;
}
else if(nDiMuonMassBin == 9){
h_JpsiMassResonance->Fill(LorVec_DiMuon.M());
h_CharmoniumResonance->Fill(LorVec_DiMuon.M());
h_ReconstructBmassInCharmoniumResonance->Fill(LorVec_ReconB.M());
Is_Phi_JPsi_Psi2S = 1;
}
else if(nDiMuonMassBin == 12){
h_Psi2sMassResonance->Fill(LorVec_DiMuon.M());
h_CharmoniumResonance->Fill(LorVec_DiMuon.M());
h_ReconstructBmassInCharmoniumResonance->Fill(LorVec_ReconB.M());
Is_Phi_JPsi_Psi2S = 1;
}
else{
h_JpsiMassCut->Fill(LorVec_DiMuon.M());
h_ReconstructBmassOffResonance->Fill(LorVec_ReconB.M());
nCounterJpsiResonance++;
}
/*if((LorVec_DiMuon.M2() < 8e6 || LorVec_DiMuon.M2() > 11e6) && (LorVec_DiMuon.M2() < 12.5e6 || LorVec_DiMuon.M2() > 15e6) && (LorVec_DiMuon.M2() < 0.98e6 || LorVec_DiMuon.M2() > 1.1e6)){
h_JpsiMassCut->Fill(LorVec_DiMuon.M());
nCounterJpsiResonance++;
}*/
///////////////////////////////////////
if(Kst2Kpluspi0Merged || Kst2Kpluspi0Resolved){
// if (nDiMuonMassBin == 9 || nDiMuonMassBin == 12) continue;
nCounterB0KsmumuMass++;
//h_BplusMass_B02KsmumuMass->Fill(B_plus_M);
}
///very soft cuts to ensure a more even distribution for the BDT training:
if(Kst2Kspiplus){
if((UseDTF ? LorVec_B_plus.Pt() : B_plus_PT) > 55000.)continue;
if((UseDTF ? LorVec_K_short.Pt() : K_short_PT) > (nKshortDecayInVeLo == 1 ? 9000. : 14000.))continue;
if(TMath::Log(B_plus_IPCHI2_OWNPV) < -8.)continue;
}
//get subset number from modulus(3) of the eventnumber
//range 0 to 2
nSubSetNumber = nEventNumber % 3;
//plot the double misIDed reconstructed mass distributions for Kshort Channel:
if(Kst2Kspiplus){
dJpsiMisIDed_M = LorVec_Jpsi_From_DoubleMisIDed.M();
h_Jpsi_Spectrum_MisIDed_Muon->Fill(dJpsiMisIDed_M);
h_BplusM_vs_BplusM_DoubleMisIDed->Fill(LorVec_B_plus_From_DoubleMisIDed.M(), dReconstructedBM);
h_BplusM_vs_JpsiM_DoubleMisIDed->Fill(LorVec_B_plus_From_DoubleMisIDed.M(), dJpsiMisIDed_M);
h_KstarM_vs_JpsiM_DoubleMisIDed->Fill(LorVec_K_star_plus_From_DoubleMisIDed.M(), dJpsiMisIDed_M);
if(nDiMuonMassBin != 2 && nDiMuonMassBin != 9 && nDiMuonMassBin != 12){
h_Q2_Spectrum_MisIDed_Muon->Fill(LorVec_Jpsi_From_DoubleMisIDed.M2() * 1e-6);
h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->Fill(dJpsiMisIDed_M);
}
}
//remove double mis-identified events:
if(Kst2Kspiplus){
if(TMath::Abs(dJpsiMisIDed_M - PDGMASS.J_PSI) < 50 && TMath::Abs(LorVec_K_star_plus_From_DoubleMisIDed.M() - PDGMASS.K_STAR_PLUS) < 30){
n_DoubleMisID++;
continue;
}
}
//plot for 2D histos
h_BplusMass_DiMuonMass->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
if(Kst2Kspiplus){
if(nKshortDecayInVeLo == 1)h_BplusMass_DiMuonMassLL->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
else h_BplusMass_DiMuonMassDD->Fill(UseDTF ? B_plus_M_DTF : B_plus_M, dDiMuonM);
}
dQ2 = LorVec_DiMuon.M2();
//---------------------------------------------------------------------------------
summary_tree->Fill();
}
output->cd();
coutInfo("Remove " + to_string(n_DoubleMisID)+ " double mis-IDed events");
coutInfo("Save " + to_string(summary_tree->GetEntries())+ " events to new file!");
summary_tree->Write("",TObject::kWriteDelete);
output->Close();
coutInfo("Save quality histograms to seperate file!");
histos->cd();
h_JpsiMass->Write("",TObject::kWriteDelete);
h_JpsiMassCut->Write("",TObject::kWriteDelete);
h_JpsiMassResonance->Write("",TObject::kWriteDelete);
h_Psi2sMassResonance->Write("",TObject::kWriteDelete);
h_CharmoniumResonance->Write("",TObject::kWriteDelete);
h_ReconstructBmassInCharmoniumResonance->Write("",TObject::kWriteDelete);
h_ReconstructBmassOffResonance->Write("",TObject::kWriteDelete);
h_KstarMass->Write("",TObject::kWriteDelete);
h_KstarMasscut->Write("",TObject::kWriteDelete);
h_KstarMass2D->Write("",TObject::kWriteDelete);
h_BplusMass_B02KsmumuMass->Write("",TObject::kWriteDelete);
h_KshortDaughters->Write("",TObject::kWriteDelete);
h_Jpsi_Spectrum_MisIDed_Muon->Write("",TObject::kWriteDelete);
h_Jpsi_Spectrum_ResonantVeto_MisIDed_Muon->Write("",TObject::kWriteDelete);
h_Q2_Spectrum_MisIDed_Muon->Write("",TObject::kWriteDelete);
h_BplusM_vs_BplusM_DoubleMisIDed->Write("",TObject::kWriteDelete);
h_BplusM_vs_JpsiM_DoubleMisIDed->Write("",TObject::kWriteDelete);
h_KstarM_vs_JpsiM_DoubleMisIDed->Write("",TObject::kWriteDelete);
h_pi0_merged_candidates_M->Write("",TObject::kWriteDelete);
h_pi0_resolved_candidates_M->Write("",TObject::kWriteDelete);
h_pi0_merged_candidates_pT->Write("",TObject::kWriteDelete);
h_pi0_resolved_candidates_pT->Write("",TObject::kWriteDelete);
h_KGammaGammaMuMu_mass_before_preselection->Write("",TObject::kWriteDelete);
h_KGammaGammaMuMu_mass_after_preselection->Write("",TObject::kWriteDelete);
h_KPiMuMu_KGammaGammaMuMu_mass->Write("",TObject::kWriteDelete);
h_pi0Mass2D->Write("",TObject::kWriteDelete);
h_diMuonReconstructedBMass2D->Write("",TObject::kWriteDelete);
h_diMuonReconstructedBMass2DFixPiZero->Write("",TObject::kWriteDelete);
h_BplusMass_DiMuonMass->Write("",TObject::kWriteDelete);
h_BplusMass_DiMuonMassDD->Write("",TObject::kWriteDelete);
h_BplusMass_DiMuonMassLL->Write("",TObject::kWriteDelete);
h_BplusMass_DiMuonMassRaw->Write("",TObject::kWriteDelete);
h_BplusMass_DiMuonMassRawDD->Write("",TObject::kWriteDelete);
h_BplusMass_DiMuonMassRawLL->Write("",TObject::kWriteDelete);
h_BeforePIDcuts_ctl->Write("",TObject::kWriteDelete);
h_BeforePIDcuts_ctk->Write("",TObject::kWriteDelete);
h_BeforePIDcuts_phi->Write("",TObject::kWriteDelete);
h_AfterPIDcuts_ctl->Write("",TObject::kWriteDelete);
h_AfterPIDcuts_ctk->Write("",TObject::kWriteDelete);
h_AfterPIDcuts_phi->Write("",TObject::kWriteDelete);
histos->Close();
coutInfo("All done for " + year + " " + magnet + string(MC ? (PHSP ? " (PHSP)" : (ReferenceChannel ? " (RefMC) " : " (MC)") ) : "") + "!");
delete output;
delete histos;
//Add TM branches for TMVA running
if(!MC) copyDataForTMVA(year,magnet);
return 1;
}
//run multiple or all functions:
int runAllYearsAndMagPol(bool MC = false, bool ReferenceChannel = false, bool PHSP = false, bool B0 = false, bool K1 = false, bool Inc = false, Int_t Run = 1){
std::vector<string> years = yearsVector(MC,ReferenceChannel,PHSP,Run);
string Magnets[] = {"down", "up"};
for(unsigned int y = 0; y < years.size(); y++){ //loop over all years
for(string & magnet : Magnets){ //loop over both magnet polaritys
if(preSelection(years.at(y), magnet, MC, ReferenceChannel, PHSP, B0, K1, Inc) == 0){
coutERROR("Failed for " + years.at(y) + " " + magnet);
return 0;
}
}
}
return 1;
}
//signal channels
int runAllSignalData(Int_t Run = 1){
if(runAllYearsAndMagPol(false, false, false, false, false, false, Run) == 0){
coutERROR("Failed for signal data sample of decay " + TheDecay + "!");
return 0;
}
coutInfo("All signal data for decay " + TheDecay + " successfully processed!!");
return 1;
}
int runAllSignalMC(Int_t Run = 1){
if(runAllYearsAndMagPol(true, false, false, false, false, false, Run) == 0){
coutERROR("Failed for signal MC sample of decay " + TheDecay + "!");
return 0;
}
coutInfo("All signal MC for decay " + TheDecay + " successfully processed!!");
return 1;
}
//reference MC channels
int runAllRefMC(Int_t Run = 1){
if(runAllYearsAndMagPol(true, true, false, false, false, false, Run) == 0){
coutERROR("Failed for reference MC sample of decay " + TheDecay + "!");
return 0;
}
coutInfo("All reference MC for decay " + TheDecay + " successfully processed!!");
return 1;
}
//reference MC channels
int runAllPHSPMC(Int_t Run = 1){
if(runAllYearsAndMagPol(true, false, true, false, false, false, Run) == 0){
coutERROR("Failed for PSHP MC sample of decay " + TheDecay + "!" );
return 0;
}
coutInfo("All PHSP MC for decay " + TheDecay +" successfully processed!!" );
return 1;
}
//B0toKst MC channels
int runAllB0MC(Int_t Run = 1, bool ReferenceChannel = false){
//BtoXMuMu not available at the moment, set ReferenceChannel to true
if(runAllYearsAndMagPol(true, ReferenceChannel, false, true, false, false, Run) == 0){
coutERROR("Failed for inclusive MC sample of decay " + TheDecay + "!" );
// return 0;
}
coutInfo("All inclusive MC for decay " + TheDecay +" successfully processed!!" );
return 1;
}
//BtoK1mumu MC channels
int runAllK1MC(bool ReferenceChannel = false){
if(preSelection("2012" , "down", true, ReferenceChannel, false, false, true, false) == 0){
coutERROR("Failed for K1 MC sample of decay " + TheDecay + "!" );
return 0; // Easier than specifying only 2012 available
}
if(preSelection("2012" , "up", true, ReferenceChannel, false, false, true, false) == 0){
coutERROR("Failed for K1 MC sample of decay " + TheDecay + "!" );
return 0; // Easier than specifying only 2012 available
}
coutInfo("All K1 MC for decay " + TheDecay +" successfully processed!!" );
return 1;
}
//Inclusive MC channels
int runAllIncMC(Int_t Run = 1){
//BtoXMuMu not available at the moment, set ReferenceChannel to true
bool ReferenceChannel = true;
string Magnets[] = {"down", "up"};
for(auto yr: yearsInc(Run)){ //loop over all years
for(string & magnet : Magnets){ //loop over both magnet polaritys
if(preSelection(yr, magnet, true, ReferenceChannel, false, false, false, true) == 0){
coutERROR("Failed for " + yr + " " + magnet);
return 0;
}
}
}
coutInfo("All inclusive MC for decay " + TheDecay +" successfully processed!!" );
return 1;
}
//All background MC channels
int runAllBkgMC(){
runAllB0MC(1,0);
runAllB0MC(1,1);
runAllB0MC(2,0);
runAllB0MC(2,1);
runAllK1MC(0); //only 2012 available
runAllK1MC(1); //only 2012 available
runAllIncMC(12);
coutInfo("All inclusive MC for decay " + TheDecay +" successfully processed!!" );
return 1;
}
int getTCKs(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false) {
gStyle -> SetOptStat(0);
LHCbStyle();
gROOT->SetBatch(kTRUE);
TChain* tree = nullptr;
if(ReferenceChannel)MC = true;
if(PHSP)MC = true;
if(ReferenceChannel && PHSP){
coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
PHSP = false;
}
//Kst2Kpluspi0Resolved case
if(Kst2Kpluspi0Resolved){
tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
}
//Kst2Kpluspi0Merged case
if(Kst2Kpluspi0Merged){
tree=new TChain("b2KstKpi0mumuMergedTuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
}
//Kst2Kspiplus case
if(Kst2Kspiplus){
tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
}
if(tree == nullptr){
coutERROR("No files found!");
return 0;
}
int N = tree->GetEntries();
if(N == 0)coutERROR("No events found in file!");
else coutInfo("Old files contain " + to_string(N)+ " events");
//Disable all branches
tree->SetBranchStatus("*",0);
//activate needed branch
UInt_t HLT1TCK;
tree->SetBranchStatus("HLT1TCK",1);
tree->SetBranchAddress("HLT1TCK", &HLT1TCK);
coutInfo("Get Min and Max from HLT1TCK branch");
Double_t HLT1TCKMin = tree->GetMinimum("HLT1TCK");
Double_t HLT1TCKMax = tree->GetMaximum("HLT1TCK");
UInt_t runNumber;
tree->SetBranchStatus("runNumber",1);
tree->SetBranchAddress("runNumber", &runNumber);
coutInfo("Get Min and Max from runNumber branch");
Double_t runNumberMin = tree->GetMinimum("runNumber");
Double_t runNumberMax = tree->GetMaximum("runNumber");
TCanvas * cRunNumberTCKs = new TCanvas("RunNumberTCKs", "RunNumberTCKs");
cRunNumberTCKs->cd();
TH2D * hRunNumberTCK = new TH2D(Form("%s_RunNumber_vs_TCKs_used_%s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()),Form("%s: RunNumber vs. TCKs used in %s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()), 100, runNumberMin, runNumberMax, 100, HLT1TCKMin, HLT1TCKMax);
//vector to count total number of each TCK. Index 0 is the sum of all TCKs, so should be equal to number of events
std::vector<Int_t>TCKcounts;
std::vector<UInt_t>FoundTCKs;
coutInfo("Start loop over " + to_string(N) + " events!");
//Fill in the first entry
tree->GetEntry(0);
hRunNumberTCK->Fill(runNumber, HLT1TCK);
FoundTCKs.push_back(HLT1TCK);
TCKcounts.push_back(1);
//Check for different TCKs in a loop
for(int ev = 1; ev < N; ev++){
tree->GetEntry(ev);
hRunNumberTCK->Fill(runNumber, HLT1TCK);
if(ev%100000 == 0)std::cout << Form("[COUNTER]\t\t %i / %i", ev, N) << std::endl;
unsigned int TCKindex = 0;
//If the TCK is different to all previous TCKs, fill it into FoundTCKs
while(FoundTCKs.at(TCKindex) != HLT1TCK){
TCKindex++;
if(FoundTCKs.size() == TCKindex){
FoundTCKs.push_back(HLT1TCK);
TCKcounts.push_back(0);
TCKindex--;
break;
}
}
//Add count to a given TCK
TCKcounts.at(TCKindex)++;
}
if(TCKcounts.size() != FoundTCKs.size()){
coutERROR("Number of found TCKs and number of counters is not equal! Exit program...");
return 0;
}
unsigned int i = 1;//, j = 0;
Int_t TmpTCKcount, TmpFoundTCK;
// order TCKs in order (says David)
while(i < TCKcounts.size()){
unsigned int j = i;
while(FoundTCKs.at(j) < FoundTCKs.at(j - 1)){
TmpFoundTCK = FoundTCKs.at(j);
FoundTCKs.at(j) = FoundTCKs.at(j - 1);
FoundTCKs.at(j - 1) = TmpFoundTCK;
TmpTCKcount = TCKcounts.at(j);
TCKcounts.at(j) = TCKcounts.at(j - 1);
TCKcounts.at(j - 1) = TmpTCKcount;
j--;
if(j ==0)break;
}
i++;
}
bool print_latex = true;
if(print_latex)
std::cout << year << " " << magnet;
else{
coutInfo("Number of events of each used TCKs in " + year + " " + magnet +":");
std::cout << "[INFO]" << std::endl;
coutInfo("TCK\tEvtCounter\tPercentage");
std::cout << "[INFO]" << std::endl;
}
for(unsigned int index = 0; index < TCKcounts.size(); index++){
//std::cout << Form("[INFO]\t\t%d \t\t%d\t\t%.2f", FoundTCKs.at(index), TCKcounts.at(index), 100./N*TCKcounts.at(index)) << std::endl;
if(print_latex)
std::cout << "\t&" << hex << FoundTCKs.at(index) << "\t&" << dec << TCKcounts.at(index) << "\t&" << std::setprecision(2) << std::fixed << 100./N*TCKcounts.at(index) << "\\\\" << std::endl;
else
std::cout << "[INFO]" << hex << FoundTCKs.at(index) << "\t" << dec << TCKcounts.at(index) << "\t" << 100./N*TCKcounts.at(index) << std::endl;
}
std::cout << "[INFO]" << std::endl;
TCanvas * cTCKs = new TCanvas("TCKs", "TCKs");
TH1D * hTCKs = new TH1D(Form("%s_TCKs_used_%s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()),Form("%s: TCKs used in %s%s", TheDecay.c_str(), year.c_str(), magnet.c_str()), TCKcounts.size(), 0, TCKcounts.size());
//Label TCK bins with TCK
for(unsigned bin = 0; bin < TCKcounts.size(); bin++){
for(int counts = 0; counts < TCKcounts.at(bin); counts++){
hTCKs->Fill(bin);
}
hTCKs->GetXaxis()->SetBinLabel(bin + 1, Form("%x", FoundTCKs.at(bin)));
//hTCKs->GetXaxis()->ChangeLabel(bin + 1, 45., -1., -1, -1, -1, Form("%i", FoundTCKs.at(bin)));
}
//Save the TCK histograms
gStyle->SetOptStat(0);
string suffix = "";
if(MC){
if(PHSP)suffix.append("_PHSP");
else if(ReferenceChannel)suffix.append("_Ref");
else suffix.append("_MC");
}
cTCKs->cd();
hTCKs->GetXaxis()->SetTitle("TCKid");
hTCKs->GetYaxis()->SetTitle("Events / TCK");
hTCKs->SetMarkerStyle(kFullCircle);
hTCKs->Draw("PE");
cTCKs->SaveAs(Form("%s/TCKinfo/TCKs%s%s%s%s.eps",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));
cTCKs->SaveAs(Form("%s/TCKinfo/TCKs%s%s%s%s.root",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));
cRunNumberTCKs->cd();
hRunNumberTCK->GetYaxis()->SetTitle("TCKid");
hRunNumberTCK->GetXaxis()->SetTitle("RunNumber");
hRunNumberTCK->Draw("COLZ");
cRunNumberTCKs->SaveAs(Form("%s/TCKinfo/RunNumberTCKs%s%s%s%s.eps",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));
cRunNumberTCKs->SaveAs(Form("%s/TCKinfo/RunNumberTCKs%s%s%s%s.root",thePath.c_str(),TheDecay.c_str(),year.c_str(),magnet.c_str(),suffix.c_str()));
return 1;
}
int ExclusiveTriggerLineSelection(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false) {
gStyle -> SetOptStat(0);
LHCbStyle();
gROOT->SetBatch(kTRUE);
TChain* tree = nullptr;
if(ReferenceChannel)MC = true;
if(PHSP)MC = true;
if(ReferenceChannel && PHSP){
coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
PHSP = false;
}
//Kst2Kpluspi0Resolved case
if(Kst2Kpluspi0Resolved){
tree=new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
}
//Kst2Kpluspi0Merged case
if(Kst2Kpluspi0Merged){
tree=new TChain("b2KstKpi0mumuMergedTuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
}
//Kst2Kspiplus case
if(Kst2Kspiplus){
tree=new TChain("b2KstKs0pimumu_Tuple/DecayTree");
tree->Add(GetInputFile(year,magnet,false,MC,ReferenceChannel,PHSP,smallSample).c_str());
}
if(tree == nullptr){
coutERROR("No files found!");
return 0;
}
int N = tree->GetEntries();
if(N == 0)coutERROR("No events found in file!");
else coutInfo("Old files contain " + to_string(N)+ " events");
//Disable all branches
tree->SetBranchStatus("*",0);
//activate needed branch
//L0
Bool_t B_plus_L0MuonDecision_TOS = false;
//HLT1
Bool_t B_plus_Hlt1TrackAllL0Decision_TOS = false;
Bool_t B_plus_Hlt1TrackMuonDecision_TOS = false;
//HLT2
Bool_t B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2Topo2BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2Topo3BodyBBDTDecision_TOS = false;
Bool_t B_plus_Hlt2SingleMuonDecision_TOS = false;
Bool_t B_plus_Hlt2DiMuonDetachedDecision_TOS = false;
tree -> SetBranchAddress( "B_plus_L0MuonDecision_TOS" , &B_plus_L0MuonDecision_TOS );
if(year == "2011" || year == "2012"){
tree -> SetBranchAddress( "B_plus_Hlt1TrackAllL0Decision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyBBDTDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );
}
else if(year == "2015" || year == "2016" || year == "2017" || year == "2018"){
tree -> SetBranchAddress( "B_plus_Hlt1TrackMVADecision_TOS" , &B_plus_Hlt1TrackAllL0Decision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo2BodyDecision_TOS" , &B_plus_Hlt2Topo2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2Topo3BodyDecision_TOS" , &B_plus_Hlt2Topo3BodyBBDTDecision_TOS );
if(year == "2016" || year == "2017" || year == "2018"){ //found to be FALSE for all 2015 data
tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu2BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2TopoMuMu3BodyDecision_TOS" , &B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS );
}
}
else{
coutERROR("Incorrect year given for the assosiation of HLT1 lines with branches! Year given: " + year);
return 0;
}
tree -> SetBranchAddress( "B_plus_Hlt1TrackMuonDecision_TOS" , &B_plus_Hlt1TrackMuonDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2SingleMuonDecision_TOS" , &B_plus_Hlt2SingleMuonDecision_TOS );
tree -> SetBranchAddress( "B_plus_Hlt2DiMuonDetachedDecision_TOS" , &B_plus_Hlt2DiMuonDetachedDecision_TOS );
//count exclusive trigger decisions:
Int_t ex_Hlt1TrackAllL0 = 0,
ex_Hlt1TrackMuon = 0;
Int_t ex_Hlt2SingleMuon = 0,
ex_Hlt2DiMuonDetached = 0,
ex_Hlt2TopoMu2Body = 0,
ex_Hlt2TopoMu3Body = 0,
ex_Hlt2Topo2Body = 0,
ex_Hlt2Topo3Body = 0,
ex_Hlt2TopoMuMu2Body = 0,
ex_Hlt2TopoMuMu3Body = 0;
coutInfo("Start loop over " +to_string(N)+ " events!");
for(int ev = 0; ev < N; ev++){
tree->GetEntry(ev);
if(ev%100000 == 0)std::cout << Form("[COUNTER]\t\t %i / %i", ev, N) << std::endl;
//HLT1
if(B_plus_Hlt1TrackAllL0Decision_TOS && !B_plus_Hlt1TrackMuonDecision_TOS)
ex_Hlt1TrackAllL0++;
if(!B_plus_Hlt1TrackAllL0Decision_TOS && B_plus_Hlt1TrackMuonDecision_TOS)
ex_Hlt1TrackMuon++;
//HLT2
if(B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2SingleMuon++;
else if(!B_plus_Hlt2SingleMuonDecision_TOS && B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2DiMuonDetached++;
else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2TopoMu2Body++;
else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2TopoMu3Body++;
else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2Topo2Body++;
else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2Topo3Body++;
else if(year == "2016" || year == "2017" || year == "2018"){
if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2TopoMuMu2Body++;
else if(!B_plus_Hlt2SingleMuonDecision_TOS && !B_plus_Hlt2DiMuonDetachedDecision_TOS && !B_plus_Hlt2TopoMu2BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMu3BodyBBDTDecision_TOS && !B_plus_Hlt2Topo2BodyBBDTDecision_TOS && !B_plus_Hlt2Topo3BodyBBDTDecision_TOS && !B_plus_Hlt2TopoMuMu2BodyBBDTDecision_TOS && B_plus_Hlt2TopoMuMu3BodyBBDTDecision_TOS)
ex_Hlt2TopoMuMu3Body++;
}
}
bool print_latex = true;
std::cout << std::endl;
std::cout << year << " " << magnet << ":" << std::endl;
std::cout << "Trigger line\texclusive events\tpercentage" << std::endl;
std::cout << "HLT1" << std::endl;
std::cout << "Hlt1TrackAllL0\t\t" << (print_latex ? " &" : "") << ex_Hlt1TrackAllL0 << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt1TrackAllL0 / N << std::endl;
std::cout << "Hlt1TrackMuon\t\t" << (print_latex ? " &" : "") << ex_Hlt1TrackMuon << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt1TrackMuon / N << std::endl;
std::cout << "HLT2" << std::endl;
std::cout << "Hlt2SingleMuon\t\t" << (print_latex ? " &" : "") << ex_Hlt2SingleMuon << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2SingleMuon / N << std::endl;
std::cout << "Hlt2DiMuonDetached\t" << (print_latex ? " &" : "") << ex_Hlt2DiMuonDetached << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2DiMuonDetached / N << std::endl;
std::cout << "Hlt2TopoMu2Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMu2Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMu2Body / N << std::endl;
std::cout << "Hlt2TopoMu3Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMu3Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMu3Body / N << std::endl;
std::cout << "Hlt2Topo2Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2Topo2Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2Topo2Body / N << std::endl;
std::cout << "Hlt2Topo3Body\t\t" << (print_latex ? " &" : "") << ex_Hlt2Topo3Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2Topo3Body / N << std::endl;
std::cout << "Hlt2TopoMuMu2Body\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMuMu2Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMuMu2Body / N << std::endl;
std::cout << "Hlt2TopoMuMu3Body\t" << (print_latex ? " &" : "") << ex_Hlt2TopoMuMu3Body << "\t" << (print_latex ? " &" : "") << 100.0 * ex_Hlt2TopoMuMu3Body / N << std::endl;
return 1;
}
Int_t plotHistos(string year = "2011", string magnet = "down", bool MC = false, bool ReferenceChannel = false, bool PHSP = false){
//Reference channel and PHSP are both MC data
if(ReferenceChannel)MC = true;
if(PHSP)MC = true;
if(ReferenceChannel && PHSP){
coutWarning("Cannot set boolean of reference channel and phase-space MC at the same time! Process Reference Channel!");
PHSP = false;
}
gStyle -> SetOptStat(0);
LHCbStyle();
TFile * histos = nullptr;
if(Kst2Kpluspi0Resolved){
if(!MC)histos = new TFile(Form("%s/data/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
else{
if(ReferenceChannel)histos = new TFile(Form("%s/data/MC/RefKplusPi0Resolved/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
else if(PHSP)histos = new TFile(Form("%s/data/PHSP/KplusPi0Resolved/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
else histos = new TFile(Form("%s/data/MC/KplusPi0Resolved/%s%s/%s%s_pi0Resolved_histos.root",path_to_output_KplusPizero.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
}
}
if(Kst2Kspiplus){
if(!MC)histos = new TFile(Form("%s/data/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
else{
if(ReferenceChannel)histos = new TFile(Form("%s/data/MC/RefKshortPiplus/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
else if(PHSP)histos = new TFile(Form("%s/data/PHSP/KshortPiplus/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
else histos = new TFile(Form("%s/data/MC/KshortPiplus/%s%s/%s%s_piplus_histos.root",path_to_output_KshortPiplus.c_str(), year.c_str(),magnet.c_str(),year.c_str(),magnet.c_str()),"OPEN");
}
}
if(histos == nullptr){
coutERROR("File with histograms not found or opened! Exit...");
return 0;
}
TCanvas * c1 = new TCanvas("c1", "c1");
c1->cd();
TH2D * h_KstarM_vs_JpsiM_DoubleMisIDed = nullptr;
h_KstarM_vs_JpsiM_DoubleMisIDed = (TH2D * ) histos->Get("h_KstarM_vs_JpsiM_DoubleMisIDed");
if(h_KstarM_vs_JpsiM_DoubleMisIDed == nullptr){
coutERROR("Histogram 'h_KstarM_vs_JpsiM_DoubleMisIDed' not found file! Exit...");
return 0;
}
h_KstarM_vs_JpsiM_DoubleMisIDed->GetXaxis()->SetTitle((UseDTF ? "m([#mu^{+}]_{#pi^{+}}K_{S}^{0})" : "m([#mu^{+}]_{#pi^{+}}#pi^{-}#pi^{+})"));
h_KstarM_vs_JpsiM_DoubleMisIDed->GetYaxis()->SetTitle("m([#pi^{+}]_{#mu^{+}}#mu^{-})");
h_KstarM_vs_JpsiM_DoubleMisIDed->Draw("COLZ");
//insert veto-region as overlay:
Float_t xLow = PDGMASS.K_STAR_PLUS - 30.;
Float_t xHigh = PDGMASS.K_STAR_PLUS + 30.;
Float_t yLow = PDGMASS.J_PSI - 50.;
Float_t yHigh = PDGMASS.J_PSI + 50.;
TGraphErrors * box = new TGraphErrors();
box->SetPoint(0, PDGMASS.K_STAR_PLUS, PDGMASS.J_PSI);
box->SetPointError(0, 30., 50.);
box->SetFillColor(kRed);
box->SetFillStyle(3445);
box->Draw("2SAME");
TLine * boxframe[4];
boxframe[0] = new TLine(xLow, yLow, xHigh, yLow);
boxframe[1] = new TLine(xHigh, yLow, xHigh, yHigh);
boxframe[2] = new TLine(xHigh, yHigh, xLow, yHigh);
boxframe[3] = new TLine(xLow, yHigh, xLow, yLow);
for(int i = 0; i < 4; i++){
boxframe[i]->SetLineWidth(2);
boxframe[i]->SetLineColor(kRed);
boxframe[i]->Draw("SAME");
}
c1->Print(Form("%s/PreSelectionPlots/%s_%s%s%s%s%s_DoubleMisID_2Dplot.eps",
thePath.c_str(),
TheDecay.c_str(),
year.c_str(),
magnet.c_str(),
(MC ? "_MC" : ""),
(ReferenceChannel ? "_Ref" : ""),
(PHSP ? "_PHSP" : "")));
TCanvas * c2 = new TCanvas("c2", "c2");
c2->cd();
TH1D * h_BplusMass_B02KsmumuMass = nullptr;
h_BplusMass_B02KsmumuMass = (TH1D * ) histos->Get("h_BplusMass_B02KsmumuMass");
if(h_BplusMass_B02KsmumuMass == nullptr){
coutERROR("Histogram 'h_BplusMass_B02KsmumuMass' not found file! Exit...");
return 0;
}
h_BplusMass_B02KsmumuMass->GetXaxis()->SetTitle((UseDTF ? "m(K_{S}^{0}#mu^{+}#mu^{-})" : "m(#pi^{-}#pi^{+}#mu^{+}#mu^{-})"));
h_BplusMass_B02KsmumuMass->GetYaxis()->SetTitle("candidates / 20MeV");
Float_t MaxContent = -1.;
for(int b = 1; b <= h_BplusMass_B02KsmumuMass->GetXaxis()->GetNbins(); b++){
if(h_BplusMass_B02KsmumuMass->GetBinContent(b) > MaxContent)MaxContent = h_BplusMass_B02KsmumuMass->GetBinContent(b);
}
h_BplusMass_B02KsmumuMass->GetYaxis()->SetRangeUser(0, MaxContent * 1.2);
h_BplusMass_B02KsmumuMass->Draw();
TGraphErrors * box2 = new TGraphErrors();
box2->SetPoint(0, PDGMASS.B_ZERO, MaxContent * 0.6);
box2->SetPointError(0, 50., MaxContent * 0.6);
box2->SetFillColor(kRed);
box2->SetFillStyle(3445);
box2->Draw("2SAME");
TLine * vetoframe[2];
vetoframe[0] = new TLine(PDGMASS.B_ZERO - 50., 0., PDGMASS.B_ZERO - 50., MaxContent * 1.2);
vetoframe[1] = new TLine(PDGMASS.B_ZERO + 50., 0., PDGMASS.B_ZERO + 50., MaxContent * 1.2);
for(int i = 0; i < 2; i++){
vetoframe[i]->SetLineWidth(2);
vetoframe[i]->SetLineColor(kRed);
vetoframe[i]->Draw("SAME");
}
c2->Print(Form("%s/PreSelectionPlots/%s_%s%s%s%s%s_B0plusRandomPion.eps",
thePath.c_str(),
TheDecay.c_str(),
year.c_str(),
magnet.c_str(),
(MC ? "_MC" : ""),
(ReferenceChannel ? "_Ref" : ""),
(PHSP ? "_PHSP" : "")));
return 1;
}
void compare_PIDgen_Eventnumbers(int year = 2011, string magnet = "down"){
string subDecay = Kst2Kspiplus ? "KshortPiplus" : "KplusPi0Resolved";
TChain * treePIDgen, * tree;
if(Kst2Kspiplus){
treePIDgen = new TChain("b2KstKs0pimumu_Tuple/DecayTree");
tree = new TChain("b2KstKs0pimumu_Tuple/DecayTree");
}
else{
treePIDgen = new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
tree = new TChain("b2KstKpi0mumuResolvedTuple/DecayTree");
}
treePIDgen->Add((path_to_data + "/data/PHSP/" + subDecay + "/" + std::to_string(year) + magnet + "/*B2KstmumuPIDgen_*.root").c_str());
tree ->Add((path_to_data + "/data/PHSP/" + subDecay + "/" + std::to_string(year) + magnet + "/*B2Kstmumu_*.root").c_str());
std::cout << "------------------------------------" << std::endl;
std::cout << "----" << year << " " << magnet << "----" << std::endl;
std::cout << "------------------------------------" << std::endl;
std::cout << "PIDgen tuple entries=" << treePIDgen->GetEntries() << std::endl;
std::cout << "Normal tuple entries=" << tree ->GetEntries() << std::endl;
std::cout << "------------------------------------" << std::endl;
}
Reference in New Issue View Git Blame Copy Permalink