#include "EVENT/LCIO.h" #include "EVENT/LCRunHeader.h" #include "EVENT/LCCollection.h" #include "EVENT/LCParameters.h" #include "EVENT/ReconstructedParticle.h" #include "EVENT/MCParticle.h" #include "MarlinProcessor.h" #include "TLorentzVector.h" #include "TVector3.h" #include "TFile.h" #include "TH1F.h" #include "TH2F.h" // Marlin stuff #include // the event display using namespace lcio; MarlinProcessor aMarlinProcessor; MarlinProcessor::MarlinProcessor() : marlin::Processor("MarlinProcessor") { _nRun =0; _nEvt =0; //_myMcTree = NULL; registerProcessorParameter( "Printing" , "Print certain messages" , _printing, (int)1 ) ; // the following parameters don't do anything yet registerProcessorParameter( "IntegratedLuminosity" , "Intergrated Luminosity for plots" , _integratedLuminosity, (float)500. ) ; registerProcessorParameter( "ElectronPolarization" , "Scale to this electron polarization" , _electronPolarization, (float)0. ) ; registerProcessorParameter( "PositronPolarization" , "Scale to this positron polarization" , _positronPolarization, (float)0. ) ; std::string rootFileName; rootFileName = "MarlinProcessor.root"; registerOptionalParameter( "RootFileName" , "Name of Root file", _rootFileName, rootFileName); } void MarlinProcessor::init() { std::cout << " MarlinProcessor::init " << std::endl; //_dataStore = new DSTFileProcessDataStore(); //_dataStore->registerAnalysis("ZHmumuX"); //_dataStore->registerAnalysisCut("ZHmumuX","Preselection"); //_dataStore->registerAnalysisCut("ZHmumuX","mumu"); //_dataStore->registerAnalysisCut("ZHmumuX","Z->mumu(tight)"); //_dataStore->registerAnalysisCut("ZHmumuX","ZH->mumull"); //_dataStore->registerAnalysisCut("ZHmumuX","ZH->mumuqq"); //_dataStore->registerAnalysisCut("ZHmumuX","Selected"); //_dataStore->registerAnalysis("ZHeeX"); //_dataStore->registerAnalysisCut("ZHeeX","Preselection"); //_dataStore->registerAnalysisCut("ZHeeX","ee"); //_dataStore->registerAnalysisCut("ZHeeX","Selected"); // Histograms _rootFile = new TFile(_rootFileName.c_str(),"recreate"); _fPFA = new TH1F("fPFA", "total reco energy",1000, 0.,500.0); _fPFAnu = new TH1F("fPFAnu", "total reco energy + enu",1000, 0.,500.0); _fNtracks = new TH1F("fNtracks", "Tracks (charged PFOs)",100, 0.,100.0); _fEvis = new TH1F("fEvis", "Evis",250, 0.,500.0); _fPt = new TH1F("fPt", "Event pT",250, 0.,500.0); _fCosMis = new TH1F("fCosMis", "CosMis (charged PFOs)",100, 0.,1.0); _fNtracksP = new TH1F("fNtracksP", "Tracks (charged PFOs)",100, 0.,100.0); _fEvisP = new TH1F("fEvisP", "Evis",250, 0.,500.0); _fPtP = new TH1F("fPtP", "Event pT",250, 0.,500.0); _fCosMisP = new TH1F("fCosMisP", "CosMis (charged PFOs)",100, 0.,1.0); _fNtracksS = new TH1F("fNtracksS", "Tracks (charged PFOs)",100, 0.,100.0); _fEvisS = new TH1F("fEvisS", "Evis",250, 0.,500.0); _fPtS = new TH1F("fPtS", "Event pT",250, 0.,500.0); _fCosMisS = new TH1F("fCosMisS", "CosMis (charged PFOs)",100, 0.,1.0); _fNtracksSa = new TH1F("fNtracksSa", "Tracks (charged PFOs)",100, 0.,100.0); _fEvisSa = new TH1F("fEvisSa", "Evis",250, 0.,500.0); _fPtSa = new TH1F("fPtSa", "Event pT",250, 0.,500.0); _fCosMisSa = new TH1F("fCosMisSa", "CosMis (charged PFOs)",100, 0.,1.0); _fNtracksSb = new TH1F("fNtracksSb", "Tracks (charged PFOs)",100, 0.,100.0); _fEvisSb = new TH1F("fEvisSb", "Evis",250, 0.,500.0); _fPtSb = new TH1F("fPtSb", "Event pT",250, 0.,500.0); _fCosMisSb = new TH1F("fCosMisSb", "CosMis (charged PFOs)",100, 0.,1.0); _fNtracksSee = new TH1F("fNtracksSee", "Tracks (charged PFOs) ee",100, 0.,100.0); _fEvisSee = new TH1F("fEvisSee", "Evis ee",250, 0.,500.0); _fPtSee = new TH1F("fPtSee", "Event pT ee",250, 0.,500.0); _fCosMisSee = new TH1F("fCosMisSee", "CosMis (charged PFOs) ee",100, 0.,1.0); //NEW HISTOGRAMS _fEresPMC = new TH1F("fEresPMC", "Energy Resolution of PANDORA / MCParticle",500, -0.1, 1.5); _fEresPMC2 = new TH1F("fEresPMC2", "Energy Resolution of PANDORA / MCParticle",500, 0.6, 1.1); _fEresPMCEta = new TH1F("fEresPMCEta", "Energy Resolution of PANDORA / MCParticle as a function of Eta",500, -5., 5.); _fEresPMCEta2 = new TH1F("fEresPMCEta2", "Energy Resolution of PANDORA / MCParticle as a function of Eta",500, -5., 5.); _dEtaEresRange = new TH1D("dEtaEresRange", "Eta distribition weighted by the mean energy resolution",3, 0, 1.4); //Z _fMCZ0M = new TH1F("fMCZ0M", "Reconstructed MC Z mass",200, 0., 250.0); _fMCZ0E = new TH1F("fMCZ0E", "Reconstructed MC Z energy",200, 0., 300.0); _fMCZ0Et = new TH1F("fMCZ0Et", "Reconstructed MC Z transverse energy",200, 0., 300.0); _fMCZ0Pt = new TH1F("fMCZ0Pt", "Reconstructed MC Z transverse momentum",100, 0., 250.0); _fMCZ0Eta = new TH1F("fMCZ0Eta", "Reconstructed MC Z pseudo rapidity",100, -5., 5.0); _fPANZ0M = new TH1F("fPANZ0M", "Reconstructed PANDORA Z mass",200, 0., 250.0); _fPANZ0E = new TH1F("fPANZ0E", "Reconstructed PANDORA Z energy",200, 0., 300.0); _fPANZ0Et = new TH1F("fPANZ0Et", "Reconstructed PANDORA Z transverse energy",200, 0., 300.0); _fPANZ0Pt = new TH1F("fPANZ0Pt", "Reconstructed PANDORA Z transverse momentum",100, 0., 250.0); _fPANZ0Eta = new TH1F("fPANZ0Eta", "Reconstructed PANDORA Z pseudo rapidity",100, -5., 5.0); //HIGGS _fMCHM = new TH1F("fMCHM", "Reconstructed MC H mass",200, 0., 250.0); _fMCHE = new TH1F("fMCHE", "Reconstructed MC H energy",200, 0., 300.0); _fMCHEt = new TH1F("fMCHEt", "Reconstructed MC H transverse energy",200, 0., 300.0); _fMCHPt = new TH1F("fMCHPt", "Reconstructed MC H transverse momentum",100, 0., 250.0); _fMCHEta = new TH1F("fMCHEta", "Reconstructed MC H pseudo rapidity",100, -5., 5.0); _fPANHM = new TH1F("fPANHM", "Reconstructed PANDORA H mass",200, 0., 250.0); _fPANHE = new TH1F("fPANHE", "Reconstructed PANDORA H energy",200, 0., 300.0); _fPANHEt = new TH1F("fPANHEt", "Reconstructed PANDORA H transverse energy",200, 0., 300.0); _fPANHPt = new TH1F("fPANHPt", "Reconstructed PANDORA H transverse momentum",100, 0., 250.0); _fPANHEta = new TH1F("fPANHEta", "Reconstructed PANDORA H pseudo rapidity",100, -5., 5.0); _fEresHb = new TH1F("fEresHb", "Energy Resolution of Truth Tagged b / MCParticle b",500, -0.1, 1.5); _fEresHb->GetXaxis()->SetTitle("E_{b}/GeV"); _fEresHb->GetYaxis()->SetTitle("Events"); _fEresHbbar = new TH1F("fEresHbbar", "Energy Resolution of Truth Tagged bbar / MCParticle bbar",500, -0.1, 1.5); _fEresHbbar->GetXaxis()->SetTitle("E_{bbar}/GeV"); _fEresHbbar->GetYaxis()->SetTitle("Events"); //ERES OF Z _fEresPANMCee = new TH1F("fEresPANMCee", "Energy Resolution of PANDORA e- e+ / MCParticle e- e+",500, -0.1, 1.5); _fEresPANMCeeEta1 = new TH1F("fEresPANMCeeEta1", "Energy Resolution of PANDORA electron / MCParticle electron for |eta|<=0.7",500, -0.1, 1.5); _fEresPANMCeeEta2 = new TH1F("fEresPANMCeeEta2", "Energy Resolution of PANDORA electron / MCParticle electron for 0.7<|eta|<=1.4",500, -0.1, 1.5); _fEresPANMCeeEta3 = new TH1F("fEresPANMCeeEta3", "Energy Resolution of PANDORA electron / MCParticle electron for |eta|>1.4",500, -0.1, 1.5); _fEresPANMCem = new TH1F("fEresPANMCem", "Energy Resolution of PANDORA e- / MCParticle e-",500, -0.1, 1.5); _fEresPANMCep = new TH1F("fEresPANMCep", "Energy Resolution of PANDORA e+ / MCParticle e+",500, -0.1, 1.5); _fEresPANEem = new TH1F("fEresPANEem", "Energy of PANDORA e^{-}",500, 0., 300.); _fEresPANEem->GetXaxis()->SetTitle("E/GeV"); _fEresPANEem->GetYaxis()->SetTitle("Events"); _fEresPANEep = new TH1F("fEresPANEep", "Energy of PANDORA e+",500, 0., 300.); _fEresMCEem = new TH1F("fEresMCEem", "Energy of MC e-",500, 0., 300.); _fEresMCEep = new TH1F("fEresMCEep", "Energy of MC e+",500, 0., 300.); _fEcalEem = new TH1F("fEcalEem", "Reconstructed PANDORA Z e^{-} ECAL Energy",400, 0., 300.0); _fEcalEem->GetXaxis()->SetTitle("E/GeV"); _fEcalEem->GetYaxis()->SetTitle("Events"); _fEMEem = new TH1F("fEMEem", "Reconstructed PANDORA Z electron ECAL+LCAL+BCAL energy",400, 0., 160.0); _fHcalEem = new TH1F("fHcalEem", "Reconstructed PANDORA Z electron HCAL energy",400, 0., 20.0); _fHADEem = new TH1F("fHADEem", "Reconstructed PANDORA Z electron HCAL+LHCAL energy",400, 0., 20.0); _fEcalEep = new TH1F("fEcalEep", "Reconstructed PANDORA Z e^{+} ECAL Energy",200, 0., 300.0); _fEcalEep->GetXaxis()->SetTitle("E/GeV"); _fEcalEep->GetYaxis()->SetTitle("Events"); _fEMEep = new TH1F("fEMEep", "Reconstructed PANDORA Z positron ECAL+LCAL+BCAL energy",200, 0., 300.0); _fHcalEep = new TH1F("fHcalEep", "Reconstructed PANDORA Z positron HCAL energy",200, 0., 300.0); _fHADEep = new TH1F("fHADEep", "Reconstructed PANDORA Z positron HCAL+LHCAL energy",200, 0., 20.0); // _fM4bb = new TH1F("fM4bb", "True b-jet di-jet Invariant Mass",200, 0., 350.0); _fM4bb->GetXaxis()->SetTitle("m_{b bbar}/GeV"); _fM4bb->GetYaxis()->SetTitle("Events"); _fE4dbb = new TH1F("fE4dbb", "E of true b-jet di-jets",200, 0., 500.0); _fEt4dbb = new TH1F("fEt4dbb", "Et of true b-jet di-jets",200, 0., 500.0); _fPt4dbb = new TH1F("fPt4dbb", "Pt of true b-jet di-jets",200, 0., 500.0); _fPR4dbb = new TH1F("fPR4dbb", "Eta of true b dijets",100, -5., 5.0); _fEresHbb = new TH1F("fEresHbb", "Energy Resolution of Truth Tagged Higgs b or bbar / MC Higgs b or bbar",500, -0.1, 2.0); _fEresHbb->GetXaxis()->SetTitle("E_{b bbar}/GeV"); _fEresHbb->GetYaxis()->SetTitle("Events"); _fEresHbb1 = new TH1F("fEresHbb1", "Energy Resolution of Truth Tagged Higgs b or bbar / MC Higgs b or bbar for |eta|<=0.7",500, -0.1, 2.0); _fEresHbb2 = new TH1F("fEresHbb2", "Energy Resolution of Truth Tagged Higgs b or bbar / MC Higgs b or bbar for 0.7<|eta|<=1.4",500, -0.1, 2.0); _fEresHbb3 = new TH1F("fEresHbb3", "Energy Resolution of Truth Tagged Higgs b or bbar / MC Higgs b or bbar for |eta|>1.4",500, -0.1, 2.0); //Jet PFO Ratios _fJetResCEb = new TH1F("fJetResCEb", "b Jet Energy of Charged PFOs / Total PFOs",500, -0.5, 1.5); _fJetResCEbbar = new TH1F("fJetResCEbbar", "bbar Jet Energy of Charged PFOs / Total PFOs",500, -0.5, 1.5); _fJetResCNb = new TH1F("fJetResCNb", "b Jet Number of Charged PFOs / Total PFOs",500, -0.5, 1.5); _fJetResCNbbar = new TH1F("fJetResCNbbar", "bbar Jet Number of Charged PFOs / Total PFOs",500, -0.5, 1.5); _fJetResNEb = new TH1F("fJetResNEb", "b Jet Energy of Neutral PFOs / Total PFOs",500, -0.5, 1.5); _fJetResNEbbar = new TH1F("fJetResNEbbar", "bbar Jet Energy of Neutral PFOs / Total PFOs",500, -0.5, 1.5); _fJetResNNb = new TH1F("fJetResNNb", "b Jet Number of Neutral PFOs / Total PFOs",500, -0.5, 1.5); _fJetResNCNbbar = new TH1F("fJetResNCNbbar", "bbar Jet Number of Neutral PFOs / Total PFOs",500, -0.5, 1.5); //WOLF HISTOGRAMS _fEresWOLFEem = new TH1F("fEresWOLFEem", "Energy of WOLF e^{-}",500, 0., 300.); _fEresWOLFEem->GetXaxis()->SetTitle("E/GeV"); _fEresWOLFEem->GetYaxis()->SetTitle("Events"); _fEresWOLFMCem = new TH1F("fEresWOLFMCem", "Energy Resolution of WOLF e- / MCParticle e-",500, -0.1, 1.5); _fEresWOLFMCem->GetYaxis()->SetTitle("Events"); _fEresWOLFMCep = new TH1F("fEresWOLFMCep", "Energy Resolution of WOLF e- / MCParticle e-",500, -0.1, 1.5); _fEresWOLFMCep->GetYaxis()->SetTitle("Events"); _fWOLFZ0M = new TH1F("fWOLFZ0M", "Reconstructed WOLF Z mass",200, 0., 250.0); _fWOLFZ0E = new TH1F("fWOLFZ0E", "Reconstructed WOLF Z energy",200, 0., 300.0); _fWOLFZ0Et = new TH1F("fWOLFZ0Et", "Reconstructed WOLF Z transverse energy",200, 0., 300.0); _fWOLFZ0Pt = new TH1F("fWOLFZ0Pt", "Reconstructed WOLF Z transverse momentum",100, 0., 250.0); _fWOLFZ0Eta = new TH1F("fWOLFZ0Eta", "Reconstructed WOLF Z pseudo rapidity",100, -5., 5.0); _fEcalWOLFEem = new TH1F("fEcalWOLFEem", "Reconstructed WOLF Z e^{-} ECAL Energy",400, 0., 300.0); _fEcalWOLFEem->GetXaxis()->SetTitle("E/GeV"); _fEcalWOLFEem->GetYaxis()->SetTitle("Events"); _fEMWOLFEem = new TH1F("fEMWOLFEem", "Reconstructed WOLF Z electron ECAL+LCAL+BCAL energy",400, 0., 160.0); _fHcalWOLFEem = new TH1F("fHcalWOLFEem", "Reconstructed WOLF Z electron HCAL energy",400, 0., 20.0); _fHADWOLFEem = new TH1F("fHADWOLFEem", "Reconstructed WOLF Z electron HCAL+LHCAL energy",400, 0., 20.0); _fEcalWOLFEep = new TH1F("fEcalWOLFEep", "Reconstructed WOLF Z e^{+} ECAL Energy",200, 0., 300.0); _fEcalWOLFEep->GetXaxis()->SetTitle("E/GeV"); _fEcalWOLFEep->GetYaxis()->SetTitle("Events"); _fEMWOLFEep = new TH1F("fEMWOLFEep", "Reconstructed WOLF Z positron ECAL+LCAL+BCAL energy",200, 0., 300.0); _fHcalWOLFEep = new TH1F("fHcalWOLFEep", "Reconstructed WOLF Z positron HCAL energy",200, 0., 300.0); _fHADWOLFEep = new TH1F("fHADWOLFEep", "Reconstructed WOLF Z positron HCAL+LHCAL energy",200, 0., 20.0); _fEresWOLFEem10Eta = new TH1F("fEresWOLFEem10Eta", "Pseudo Rapidity of Reconstructed WOLF e^{-} with E<10GeV",100, -5., 5.0); _fEresWOLFEep10Eta = new TH1F("fEresWOLFEep10Eta", "Pseudo Rapidity of Reconstructed WOLF e^{+} with E<10GeV",100, -5., 5.0); //4th Year Histograms _fM12M34 = new TH2F("fM12M34", "jet 12 vs 34 masses ",200, 0., 200.0, 200, 0., 200.0); _fM12 = new TH1F("fM12", "2jet di-jet invariant mass",200, 0., 500.0); _fM34 = new TH1F("fM34", "4jet di-jet invariant mass (combination closest to mass of Higgs)",200, 0., 500.0); _fMbb = new TH1F("fMbb", "b-jet mass",200, 0., 500.0); _fMcc = new TH1F("fMcc", "c-jet mass",200, 0., 500.0); _fMone = new TH1F("fMone", "1jet invariant mass",200, 0., 5.0); _fM2bb = new TH1F("fM2bb", "2jet true b-jet di-jet invariant mass",200, 0., 500.0); _fM2cc = new TH1F("fM2cc", "2jet true c-jet di-jet invariant mass",200, 0., 500.0); _fM3bb = new TH1F("fM3bb", "3jet true b-jet di-jet invariant mass",200, 0., 500.0); _fM3cc = new TH1F("fM3cc", "3jet true c-jet di-jet invariant mass",200, 0., 500.0); _fM4cc = new TH1F("fM4cc", "4jet true c-jet di-jet invariant mass",200, 0., 500.0); _fMFT4dbb = new TH1F("fMFT4dbb", "Correct true flavour tagged b-jet di-jet invariant mass",200, 0., 500.0); _fMFT4dcc = new TH1F("fMFT4dcc", "Correct true flavour tagged c-jet di-jet invariant mass",200, 0., 500.0); _dNFTi4bb = new TH1D("dNFTi4bb", "Total number of each flavour tagged b-jets",5, 0, 5); _dNFT4bb = new TH1D("dNFT4bb", "Total number of each flavour tagged b-jets",5, 0, 5); _dNFTi4cc = new TH1D("dNFTi4cc", "Total number of each flavour tagged c-jets",4, 0, 4); _dNmFTi4bb = new TH1D("dNmFTi4bb", "Total number of incorrect flavour tagged b-jet di-jets",1, 0, 2); _dNmFTi4cc = new TH1D("dNmFTi4cc", "Total number of incorrect flavour tagged c-jet di-jets",1, 0, 2); _fPRmFT4bb = new TH1F("fPRmFT4bb", "Eta of mistagged FT b-jets",100, -5., 5.0); _fPRmFT4cc = new TH1F("fPRmFT4cc", "Eta of mistagged FT c-jets",100, -5., 5.0); _fPRmFT4dbb = new TH1F("fPRmFT4dbb", "Eta of mistagged FT true b-jet di-jets",100, -5., 5.0); _fPRmFT4dcc = new TH1F("fPRmFT4dcc", "Eta of mistagged FT true c-jet di-jets",100, -5., 5.0); _fPRFT4nb = new TH1F("fPRFT4nb", "Eta of FT true non b-jets",100, -5., 5.0); _fPRFT4nc = new TH1F("fPRFT4nc", "Eta of FT true non c-jets",100, -5., 5.0); _fPtmFT4bb = new TH1F("fPtmFT4bb", "Pt of FT mistagged b-jets",100, 0., 250.0); _fPtmFT4cc = new TH1F("fPtmFT4cc", "Pt of FT mistagged c-jets",100, 0., 250.0); _fFTPR4bb = new TH1F("fFTPR4bb", "Eta of correct FT b-jets",100, -5., 5.0); _fFTPR4cc = new TH1F("fFTPR4cc", "Eta of correct FT c-jets",100, -5., 5.0); _fFTPR4dbb = new TH1F("fFTPR4dbb", "Eta of correct FT b-jets di-jets",100, -5., 5.0); _fFTPR4dcc = new TH1F("fFTPR4dcc", "Eta of correct FT c-jets di-jets",100, -5., 5.0); _fPt4dcc = new TH1F("fPt4dcc", "Pt of true c-jets di-jets",200, 0., 500.0); _fPt4dbb2550 = new TH1F("fPt4dbb2550", "Pt of true b-jet di-jets 25-50GeV",200, 0., 500.0); _fPt4dcc2550 = new TH1F("fPt4dcc2550", "Pt of true c-jet di-jets 25-50GeV",200, 0., 500.0); _fPt4dbb5075 = new TH1F("fPt4dbb5075", "Pt of true b-jet di-jets 50-75GeV",200, 0., 500.0); _fPt4dcc5075 = new TH1F("fPt4dcc5075", "Pt of true c-jet di-jets 50-75GeV",200, 0., 500.0); _fPt4dbb75 = new TH1F("fPt4dbb75", "Pt of true b-jet di-jets 75Gev And Above",200, 0., 500.0); _fPt4dcc75 = new TH1F("fPt4dcc75", "Pt of true c-jet di-jets 75Gev And Above",200, 0., 500.0); _fPR4bb = new TH1F("fPR4bb", "Eta of true b-jets",100, -5., 5.0); _fPR4cc = new TH1F("fPR4cc", "Eta of true c-jet di-jets",100, -5., 5.0); _fEta4dbb2550 = new TH1F("fEta4dbb2550", "Eta of true b-jet di-jets 25-50GeV",100, -5.0, 5.0); _fEta4dcc2550 = new TH1F("fEta4dcc2550", "Eta of true c-jet di-jets 25-50GeV",100, -5.0, 5.0); _fEta4dbb5075 = new TH1F("fEta4dbb5075", "Eta of true b-jet di-jets 50-75GeV",100, -5.0, 5.0); _fEta4dcc5075 = new TH1F("fEta4dcc5075", "Eta of true c-jet di-jets 50-75GeV",100, -5.0, 5.0); _fEta4dbb75 = new TH1F("fEta4dbb75", "Eta of true b-jet di-jets 75Gev And Above",100, -5.0, 5.0); _fEta4dcc75 = new TH1F("fEta4dcc75", "Eta of true c-jet di-jets 75Gev And Above",100, -5.0, 5.0); _fEtaFT4dbb2550 = new TH1F("fEtaFT4dbb2550", "Eta of correct FT b-jet di-jets 25-50GeV",100, -5.0, 5.0); _fEtaFT4dcc2550 = new TH1F("fEtaFT4dcc2550", "Eta of correct FT c-jet di-jets 25-50GeV",100, -5.0, 5.0); _fEtaFT4dbb5075 = new TH1F("fEtaFT4dbb5075", "Eta of correct FT b-jet di-jets 50-75GeV",100, -5.0, 5.0); _fEtaFT4dcc5075 = new TH1F("fEtaFT4dcc5075", "Eta of correct FT c-jet di-jets 50-75GeV",100, -5.0, 5.0); _fEtaFT4dbb75 = new TH1F("fEtaFT4dbb75", "Eta of correct FT b-jet di-jets 75Gev And Above",100, -5.0, 5.0); _fEtaFT4dcc75 = new TH1F("fEtaFT4dcc75", "Eta of correct FT c-jet di-jets 75Gev And Above",100, -5.0, 5.0); _fPtFT4dbb = new TH1F("fPtFT4dbb", "Pt of correct FT b-jet di-jets",200, 0.0, 500.0); _fPtFT4dcc = new TH1F("fPtFT4dcc", "Pt of correct FT c-jet di-jets",200, 0.0, 500.0); _fPtFT4dbb2550 = new TH1F("fPtFT4dbb2550", "Pt of correct FT b-jet di-jets 25-50GeV",200, 0.0, 500.0); _fPtFT4dcc2550 = new TH1F("fPtFT4dcc2550", "Pt of correct FT c-jet di-jets 25-50GeV",200, 0.0, 500.0); _fPtFT4dbb5075 = new TH1F("fPtFT4dbb5075", "Pt of correct FT b-jet di-jets 50-75GeV",200, 0.0, 500.0); _fPtFT4dcc5075 = new TH1F("fPtFT4dcc5075", "Pt of correct FT c-jet di-jets 50-75GeV",200, 0.0, 500.0); _fPtFT4dbb75 = new TH1F("fPtFT4dbb75", "Pt of correct FT b-jet di-jets 75Gev And Above",200, 0.0, 500.0); _fPtFT4dcc75 = new TH1F("fPtFT4dcc75", "Pt of correct FT c-jet di-jets 75Gev And Above",200, 0.0, 500.0); _fFTVNNbb = new TH1F("fFTVNNbb", "FT value for Neural Net for true b-jets",100, 0.0, 1.0); //Have changed from 50 bins to 100 to look at X-section _fFTVNN4nb = new TH1F("fFTVNN4nb", "FT value for Neural Net for all non true b-jets",100, 0.0, 1.0); _fFTVNN4all = new TH1F("fFTVNN4all", "FT value for Neural Net for all jets",100, 0.0, 1.0); _fPtNN4bb = new TH1F("fPtNN4bb", "Pt of Neural Net for b-jets",200, 0.0, 500.0); _fPtNN4cc = new TH1F("fPtNN4cc", "Pt of Neural Net for c-jets",200, 0.0, 500.0); _fEt4bb = new TH1F("fEt4bb", "Et of true b-jets",100, 0., 250.0); _fE4bb = new TH1F("fE4bb", "E of true b-jets",100, 0., 250.0); _fEt4cc = new TH1F("fEt4cc", "Et of true c-jets",100, 0., 250.0); _fEtFT4bb = new TH1F("fEtFT4bb", "Et of correct FT b-jets",100, 0., 250.0); _fEtFT4cc = new TH1F("fEtFT4cc", "Et of correct FT c-jets",100, 0., 250.0); _fEtmFT4bb = new TH1F("fEtmFT4bb", "Et of mistagged FT b-jets",100, 0., 250.0); _fEtmFT4cc = new TH1F("fEtmFT4cc", "Et of mistagged FT c-jets",100, 0., 250.0); _dXsecuubb = new TH1D("dXsecuubb", "mumubb Cross section measurement required values",6, 0, 6); _dXsecuucc = new TH1D("dXsecuucc", "mumucc Cross section measurement required values",6, 0, 6); _fEtFT4nb = new TH1F("fEtFT4nb", "Et of FT true non b-jets",100, 0.0, 250.0); _fPtFT4nb = new TH1F("fPtFT4nb", "Pt of FT true non b-jets",100, 0.0, 250.0); _fPtFT4bb = new TH1F("fPtFT4bb", "Pt of correct FT b-jets",100, 0., 250.0); _fEW = new TH1F("fEW", "Event weight for each correct FT b-jet di-jet event",2, 0.0, 2.0); // Monte Carlo Distributions _fmhmcmm = new TH1F("fmhmcmm", "mhmcmm" ,200, 100., 200.0); _fmzmcmm = new TH1F("fmzmcmm", "mzmcmm" ,100, 60., 120.0); _fmzmcee = new TH1F("fmzmcee", "mzmcee" ,100, 60., 120.0); _fmhmcee = new TH1F("fmhmcee", "mhmcee" ,200, 100., 200.0); // Reconstructed Distributions _fmzmhee = new TH2F("fmzmhee", "mzmhee" ,100,60.,100.,200, 100., 200.); _fmzmhmm = new TH2F("fmzmhmm", "mzmhmm" ,100,60.,100.,200, 100., 200.); _fmzmm = new TH1F("fmzmm", "mzmm" ,100, 60., 120.0); _fmzee = new TH1F("fmzee", "mzee" ,100, 60., 120.0); _fmhee = new TH1F("fmhee", "mhee" ,200, 100., 200.0); _fmhmm = new TH1F("fmhmm", "mhmm" ,200, 100., 200.0); _fmzmmS = new TH1F("fmzmmS", "Sel mzmm" ,100, 60., 120.0); _fmzeeS = new TH1F("fmzeeS", "Sel mzee" ,100, 60., 120.0); _fmheeS = new TH1F("fmheeS", "Sel mhee" ,200, 100., 200.0); _fmhmmS = new TH1F("fmhmmS", "Sel mhmm" ,200, 100., 200.0); _fmzmmraw = new TH1F("fmzmmraw", "mzmmraw" ,100, 60., 120.0); _fmhmmraw = new TH1F("fmhmmraw", "mhmmraw" ,200, 100., 200.0); _fmzeeraw = new TH1F("fmzeeraw", "mzeeraw" ,100, 60., 120.0); _fmheeraw = new TH1F("fmheeraw", "mheeraw" ,200, 100., 200.0); // Distributions for good/bad ee events _fmheeg = new TH1F("fmheeg", "mheeg" ,200, 100., 200.0); _fmheeb = new TH1F("fmheeb", "mheeb" ,200, 100., 200.0); _fmzeeg = new TH1F("fmzeeg", "mzeeg" ,100, 60., 120.0); _fmzeeb = new TH1F("fmzeeb", "mzeeb" ,100, 60., 120.0); _fmheegraw = new TH1F("fmheegraw", "mheegraw" ,200, 100., 200.0); _fmheebraw = new TH1F("fmheebraw", "mheebraw" ,200, 100., 200.0); _fmzeegraw = new TH1F("fmzeegraw", "mzeegraw" ,100, 60., 120.0); _fmzeebraw = new TH1F("fmzeebraw", "mzeebraw" ,100, 60., 120.0); _fBremR = new TH1F("fBremR", "BremR" ,2000, 0., 2000.0); _eop = new TH1F("eop", "eop" ,200, 0., 2.0); // determine helicity fractions for incoming electron/positrons as function of polarization _helicityFractionRR = 0.25*(1+_electronPolarization)*(1+_positronPolarization); _helicityFractionRL = 0.25*(1+_electronPolarization)*(1-_positronPolarization); _helicityFractionLR = 0.25*(1-_electronPolarization)*(1+_positronPolarization); _helicityFractionLL = 0.25*(1-_electronPolarization)*(1-_positronPolarization); std::cout << " Electron/Positron Helicity fractions " << std::endl; std::cout << " RR : " << _helicityFractionRR << std::endl; std::cout << " RL : " << _helicityFractionRL << std::endl; std::cout << " LR : " << _helicityFractionLR << std::endl; std::cout << " LL : " << _helicityFractionLL << std::endl; } void MarlinProcessor::processRunHeader( LCRunHeader* run) { _nRun++ ; //std::cout << " MarlinProcessor::processRunHeader " << _nEvt << std::endl; _runNumber = run->getParameters().getIntVal("DataRunNumber"); } void MarlinProcessor::processEvent( LCEvent * evt ) { // get event weight // float weight = evt->getWeight(); _crossSection = evt->getParameters().getFloatVal("CrossSection_fb"); std::string process = evt->getParameters().getStringVal("Process"); _wgt = 1.; // for now! if(_process!=process){ //_dataStore->newProcess(process); _signal = false; _process = process; _helicityWeight = 0.0; if(_process.find("ep+1.0_em+1.0")!=std::string::npos)_helicityWeight = _helicityFractionRR; if(_process.find("ep-1.0_em+1.0")!=std::string::npos)_helicityWeight = _helicityFractionRL; if(_process.find("ep+1.0_em-1.0")!=std::string::npos)_helicityWeight = _helicityFractionLR; if(_process.find("ep-1.0_em-1.0")!=std::string::npos)_helicityWeight = _helicityFractionLL; // std::cout << " Process = " << _process << " helicity weight : " << _helicityWeight << std::endl; // std::cout << " Cross section = " << _crossSection << std::endl; int n = 1; //!!!_dataStore->getExpectedEvents(); if(n!=0)_wgt = _integratedLuminosity*_crossSection*_helicityWeight/n; _wgt = 1.; // for now! // std::cout << " Weight = " << _wgt << std::endl; //_dataStore->setWeight(_wgt); // if(_process.find("e2e2h")!=std::string::npos)_signal=true; } //_dataStore->addEvent(); // if((_nEvt/10000)*10000==_nEvt)std::cout << " MarlinProcessor::processEvent : " << _nEvt << std::endl; _nEvt++; //Count number of events _dXsecuubb->Fill(0.0,_wgt); _dXsecuucc->Fill(0.0,_wgt); //if(!_signal)return; // _eventLuminosity = 1.0/_crossSection/_helicityWeight; // clear all global variables (vectors) _MCPvec.clear(); _pfovec.clear(); _RecoParvec.clear(); _onejetvec.clear(); _twojetvec.clear(); _threejetvec.clear(); _fourjetvec.clear(); _MCele.clear(); _Pele.clear(); _Pbbbar.clear(); _Z0MC.clear(); _HbbMC.clear(); _Z0PANDORA.clear(); _onejetTruthvec.clear(); _twojetTruthvec.clear(); _threejetTruthvec.clear(); _fourjetTruthvec.clear(); _fourjetTruthPartonvec.clear(); _fourjetFTvec.clear(); _Zmmfound.clear(); // reset Monte Carlo infomation //if(_myMcTree!=NULL)delete _myMcTree; //_myMcTree = NULL; // Read the PFOs into a vector this->ReadPFOs(evt); // Analyse MC truth infomation this->AnalyseMC(evt); // Analyse Reconstructed Event this->AnalyseEvent(evt); } void MarlinProcessor::end(){ //_dataStore->close(); //_dataStore->dumpFileList(); //_dataStore->fileSummary(); //_dataStore->analysisSummary(); // End of program - write out root file _rootFile->Write(); _rootFile->Close(); std::cout << " Closed root file " << std::endl; } void MarlinProcessor::AnalyseMC( LCEvent* evt ) { _good = true; _signal = false; typedef const std::vector StringVec ; StringVec* strVec = evt->getCollectionNames() ; for(StringVec::const_iterator name=strVec->begin(); name!=strVec->end(); name++){ LCCollection* col = evt->getCollection(*name); if (col->getTypeName() == LCIO::MCPARTICLE ) { //_myMcTree = new MyMCTree(col); //_myMcTree->PerfectPFO(0); //_pMcPFOs = *(_myMcTree->getMCPFOs()); int nParticles = col->getNumberOfElements() ; // std::cout << " In MCPARTICLE section with number of elements " << nParticles << std::endl; // std::cout << " MCPARTICLE with first element address is " << col->getElementAt(0) << std::endl; // std::cout << " MCPARTICLE with second element address is " << col->getElementAt(1) << std::endl; //INCLUDING TESTING PART HERE ReconstructedParticle* MCparticle = NULL; float MCpenergy = 0.0; double MCpmomentum; double MCpmass; for(int i=0;i(col->getElementAt(i)); // MCpenergy = MCparticle->getEnergy(); // MCpmomentum = MCparticle->getMomentum(); // MCpmass = MCparticle->getMass(); // EVENT::MCParticleVec _daughters = part->getDaughters(); // EVENT::MCParticleVec _parents = part->getParents(); // // int ipdg = part->getPDG(); // std::cout << "In i = " << i << std::endl; //NOT RIGHT std::cout << "MCpmomentum[0]: " << MCparticle->getMomentum()[0] << std::endl; // std::cout << " MCpenergy: " << MCparticle->getEnergy() << ". MCpmomentum: " << MCparticle->getMomentum() << ". MCpmass: " << MCparticle->getMass() << "." << std::endl; } } } // loops over all generated particles expected to give a reconstructed particle std::vector mcMuPlus; std::vector mcMuMinus; std::vector mcEPlus; std::vector mcEMinus; //for(uint imc=0;imc<_pMcPFOs.size();imc++){ // MCParticle* mcpart = _pMcPFOs[imc]->getMCParticle(); // float px = _pMcPFOs[imc]->getMomentum()[0]; // float py = _pMcPFOs[imc]->getMomentum()[1]; // float pz = _pMcPFOs[imc]->getMomentum()[2]; // float e= sqrt(px*px+py*py+pz*pz); // EVENT::MCParticleVec parents = mcpart->getParents(); // bool ok = false; // if(parents.size()==0)ok=true; // if(parents.size()==1 && parents[0]->getPDG()==94 )ok=true; // if(parents.size()==1 && abs(parents[0]->getPDG())==11)ok=true; // if(parents.size()==1 && abs(parents[0]->getPDG())==13)ok=true; // if(ok){ // if(mcpart->getPDG()==-13){ // float e= sqrt(px*px+py*py+pz*pz); // TLorentzVector pmu(px,py,pz,e); // mcMuMinus.push_back(pmu); // } // if(mcpart->getPDG()==13){ // float e= sqrt(px*px+py*py+pz*pz); // TLorentzVector pmu(px,py,pz,e); // mcMuPlus.push_back(pmu); // } // if(mcpart->getPDG()==-11){ // float e= sqrt(px*px+py*py+pz*pz); // TLorentzVector pe(px,py,pz,e); // mcEMinus.push_back(pe); // EVENT::MCParticleVec daughters = mcpart->getDaughters(); // for(unsigned int jmc=0;jmcgetVertex()[0]; // float y = daughters[jmc]->getVertex()[1]; // float z = daughters[jmc]->getVertex()[2]; // float l = sqrt(x*x+y*y+z*z); // float r = sqrt(x*x+y*y); // //std::cout << " " << jmc << " " << daughters[jmc]->getPDG() << " " << daughters[jmc]->getEnergy() << " r = " << r << std::endl; // if( daughters[jmc]->getEnergy()>0.5 && fabs(z)<2240 && r<1850)_good=false; // if( daughters[jmc]->getEnergy()>0.5 && fabs(z)<2240)_fBremR->Fill(r,_wgt); // } // } // // if(mcpart->getPDG()==11){ // float e= sqrt(px*px+py*py+pz*pz); // TLorentzVector pe(px,py,pz,e); // mcEPlus.push_back(pe); // EVENT::MCParticleVec daughters = mcpart->getDaughters(); // for(unsigned int jmc=0;jmcgetVertex()[0]; // float y = daughters[jmc]->getVertex()[1]; // float z = daughters[jmc]->getVertex()[2]; // float r = sqrt(x*x+y*y); // float l = sqrt(x*x+y*y+z*z); // // std::cout << " " << jmc << " " << daughters[jmc]->getPDG() << " " << daughters[jmc]->getEnergy() << " r = " << r << std::endl; // if( daughters[jmc]->getEnergy()>0.5 && fabs(z)<2240 && r<1850)_good=false; // if( daughters[jmc]->getEnergy()>0.5 && fabs(z)<2240)_fBremR->Fill(r,_wgt); // } // } // } // // float cosTheta = fabs(pz)/sqrt(px*px+py*py+pz*pz); // if(cosTheta>0.995){ // if(_pMcPFOs[imc]->isPhoton() && cosTheta>0.9995){ // _mcIsrEnergy+= mcpart->getEnergy(); // } // } // _mcTotalEnergy+=mcpart->getEnergy(); //} float bestdmz = 999.; TLorentzVector pmumu; TLorentzVector pMuPlus; TLorentzVector pMuMinus; for(unsigned int i=0;i0)std::cout << " MC Z->mumu : " << pmumu.M() << " " << precoil.M() << std::endl; _fmzmcmm->Fill(pmumu.M(),_wgt); _fmhmcmm->Fill(precoil.M(),_wgt); } bestdmz = 999.; TLorentzVector pee; TLorentzVector pEPlus; TLorentzVector pEMinus; for(unsigned int i=0;i0)std::cout << " MC Z->ee : " << pee.M() << " " << precoil.M() << " " << costp << " : " << costm << std::endl; if(fabs(pee.M()-91.2)<25.)_signal = true; if(fabs(costp)>0.975||costm>0.975)_signal=false; bool fp = false; bool fm = false; for(unsigned int i=0;i<_pfovec.size();i++){ if(_pfovec[i]->getCharge()!=0){ float px = _pfovec[i]->getMomentum()[0]; float py = _pfovec[i]->getMomentum()[1]; float pz = _pfovec[i]->getMomentum()[2]; float e = _pfovec[i]->getEnergy(); float cosp = (px*pEPlus.Px()+py*pEPlus.Py()+pz*pEPlus.Pz())/e/pEPlus.E(); float cosm = (px*pEMinus.Px()+py*pEMinus.Py()+pz*pEMinus.Pz())/e/pEMinus.E(); if(cosp>0.99)fp = true; if(cosm>0.99)fm = true; } } if(!fp||!fm)_signal=false; _fmzmcee->Fill(pee.M(),_wgt); _fmhmcee->Fill(precoil.M(),_wgt); } //std::vector mcQuarks = *(_myMcTree->getQuarks()); float sumE2 = 0; //for(unsigned int i=0;igetEnergy(); // sumE2 += eq*eq; //} } void MarlinProcessor::ReadPFOs( LCEvent* evt ) { _pfovec.clear(); _RecoParvec.clear(); _MCPvec.clear(); _Pele.clear(); _Pbbbar.clear(); _MCele.clear(); _Z0MC.clear(); _HbbMC.clear(); _Z0PANDORA.clear(); typedef const std::vector StringVec ; StringVec* strVec = evt->getCollectionNames() ; for(StringVec::const_iterator name=strVec->begin(); name!=strVec->end(); name++){ LCCollection* col = evt->getCollection(*name); unsigned int nelem = col->getNumberOfElements(); if (col->getTypeName() == LCIO::RECONSTRUCTEDPARTICLE ) { if(*name=="PandoraPFANewPFOs"){ // std::cout << "In PandoraPFANewPFOs, with nelem = " << nelem << std::endl; for(unsigned int i=0;i(col->getElementAt(i)); _pfovec.push_back(recoPart); } } } } } void MarlinProcessor::AnalyseEvent( LCEvent* evt ) { _MCPvec.clear(); _pfovec.clear(); _RecoParvec.clear(); _Pele.clear(); _Pbbbar.clear(); _MCele.clear(); _Z0MC.clear(); _HbbMC.clear(); _Z0PANDORA.clear(); if(_nEvt==1)streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent" << std::endl; streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent: Event Number " << _nEvt << std::endl; ReconstructedParticle* Zmm = NULL; ReconstructedParticle* Zee = NULL; typedef const std::vector StringVec ; StringVec* strVec = evt->getCollectionNames() ; for(StringVec::const_iterator name=strVec->begin(); name!=strVec->end(); name++){ LCCollection* col = evt->getCollection(*name); unsigned int nelem = col->getNumberOfElements(); // on the first and second event print out available collections //CHANGE FOR EVERY EVENT //if(_nEvt==1 || _nEvt==2){ // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent Input Collections : " << *name << " elements : " << nelem << std::endl; if (col->getTypeName() == LCIO::CLUSTER ) { std::vector sdNames; col->getParameters().getStringVals( "ClusterSubdetectorNames" , sdNames ); for(unsigned i=0;i< sdNames.size();++i){ // if(_printing>0)std::cout << " Index = " << i << " Name = " << sdNames[i] << std::endl; //} } } int ipdg=0; int firstH = 0; EVENT::MCParticleVec _Higgs2parents; if(*name=="MCParticle"){ //if(*name=="MCParticlesSkimmed"){ for(unsigned int i=0;i(col->getElementAt(i)); ipdg = part->getPDG(); // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent MCParticle : Mass = "<< part->getMass() << " Energy = " << part->getEnergy() << " PDG = " << ipdg << std::endl; //Zee Part if(ipdg==25 && firstH==0){//checking for the first Higgs particle // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent MCParticle : Higgs FOUND !!!!!!!!!!!!!! "<< std::endl; firstH++; //parents and daughters of the Higgs EVENT::MCParticleVec _parents = part->getParents(); //beam e+e- MCParticle* beamele = dynamic_cast( _parents[0] ) ; //only need to look at either e+ or e- EVENT::MCParticleVec _daughters = beamele->getDaughters(); //H and e+e- just after scatter // std::cout<<" ndaughters = "<<_daughters.size()<getPDG()==5)_HbbMC.push_back(_daughters[id]); if(_daughters[id]->getPDG()==-5)_HbbMC.push_back(_daughters[id]); } } } } } if (col->getTypeName() == LCIO::RECONSTRUCTEDPARTICLE ) { if(*name=="PandoraPFANewPFOs"){ for(unsigned int i=0;i(col->getElementAt(i)); _pfovec.push_back(recoPart); ipdg = recoPart->getType(); //Zee Part if(ipdg==11 || ipdg==-11){//if e- or e+ // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent PandoraPFANewPFOs : Found e- or e+!!!!!!!!!!!!!! "<< std::endl; _Pele.push_back(recoPart); } if(ipdg==5 || ipdg==-5){//if b or bbar // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent PandoraPFANewPFOs : Found b or bbar!!!!!!!!!!!!!! "<< std::endl; _Pbbbar.push_back(recoPart); } // if(ipdg==25){ //check if Higgs // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent PandoraPFANewPFOs : Higgs FOUND = "<< std::endl; // } } } if (col->getTypeName() == LCIO::RECONSTRUCTEDPARTICLE ) { if(*name=="RecoParticles"){ for(unsigned int i=0;i(col->getElementAt(i)); _RecoParvec.push_back(recoPart); ipdg = recoPart->getType(); //streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent RecoParticles of type: " << ipdg << std::endl; } } } // Look for reconstructed Zs if(nelem==1){ if(*name=="ZmumuPFOs"){ Zmm = dynamic_cast(col->getElementAt(0)); //Count number of Zmm candidates _dXsecuubb->Fill(1.0,_wgt); _dXsecuucc->Fill(1.0,_wgt); _Zmmfound.push_back(1); } if(*name=="ZeePFOs"){ Zee = dynamic_cast(col->getElementAt(0)); // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent ZeePFOs : Mass = "<< Zee->getMass() << " Energy = " << Zee->getEnergy() << " PDG = " << Zee->getType() << std::endl; // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent Found a Zee candidate" <(col->getElementAt(i)); _onejetvec.push_back(recoPart); } } // look at event when reconstructed as two jets if(*name=="Durham_2Jets"){ // creates a vector of 2 pointers to the reconstructed jets for(unsigned int i=0;i(col->getElementAt(i)); _twojetvec.push_back(recoPart); } } // look at event when reconstructed as three jets if(*name=="Durham_3Jets"){ // creates a vector of 3 pointers to the reconstructed jets for(unsigned int i=0;i(col->getElementAt(i)); _threejetvec.push_back(recoPart); } } // look at event when reconstructed as four jets if(*name=="Durham_4Jets"){ // creates a vector of 4 pointers to the reconstructed jets for(unsigned int i=0;i(col->getElementAt(i)); _fourjetvec.push_back(recoPart); } } } // look at truth info if(*name=="TrueJetFlavour_1Jets"){ // creates a vector of 2 pointers to the jet flavour info for(unsigned int i=0;i(col->getElementAt(i)); _onejetTruthvec.push_back(floatVec); } } if(*name=="TrueJetFlavour_2Jets"){ // creates a vector of 2 pointers to the jet flavour info for(unsigned int i=0;i(col->getElementAt(i)); _twojetTruthvec.push_back(floatVec); } } if(*name=="TrueJetFlavour_3Jets"){ // creates a vector of 3 pointers to the jet flavour info for(unsigned int i=0;i(col->getElementAt(i)); _threejetTruthvec.push_back(floatVec); } } if(*name=="TrueJetFlavour_4Jets"){ // creates a vector of 4 pointers to the jet flavour info for(unsigned int i=0;i(col->getElementAt(i)); _fourjetTruthvec.push_back(floatVec); } } if(*name=="FlavourTag"){ //std::cout<<"FlavourTag was found with nelem " << nelem << " for event " << _nEvt << std::endl; for(unsigned int i=0;i(col->getElementAt(i)); _fourjetFTvec.push_back(floatVec); } } } // if(_onejetvec.size()==1)this->AnalyseOneJetEvent(); // if(_twojetvec.size()==2)this->AnalyseTwoJetEvent(); // if(_threejetvec.size()==3)this->AnalyseThreeJetEvent(); if(_HbbMC.size()==2)this->AnalyseFourJetEvent(); //Analyse only when MC b bar from the Higgs was found // if(_Pele.size()==1 && _MCele.size()!=0)this->AnalysePandoraPFANewPFOs(); //Shouldn't use, was made for electron gun and is wrong this->AnalyseEmEp(); //Z reconstruction and e- e+ energy resolutions //WOLF // if(_RecoParvec.size()!=0)this->AnalyseRecoParticlesEvent(); // if(_HbbMC.size()==2 && _fourjetvec.size()==4)this->AnalyseWOLFFourJetEvent(); //Analyse only when MC b bar from the Higgs was found // global variables for cuts int ntracks = 0; float px = 0.; float py=0.;float pz=0.; float e=0.; for(unsigned int i=0;i<_pfovec.size();i++){ if(_pfovec[i]->getCharge()!=0)ntracks++; px += _pfovec[i]->getMomentum()[0]; py += _pfovec[i]->getMomentum()[1]; pz += _pfovec[i]->getMomentum()[2]; e += _pfovec[i]->getEnergy(); } TLorentzVector totalMass(px,py,pz,e); //TEST PART FOR MCPvec // std::cout << " MCPvec size = " << _MCPvec.size() << std::endl; for(unsigned int i=0;i<_MCPvec.size();i++){ px += _MCPvec[i]->getMomentum()[0]; py += _MCPvec[i]->getMomentum()[1]; pz += _MCPvec[i]->getMomentum()[2]; e += _MCPvec[i]->getEnergy(); // std::cout << "for i (of) = " << i << " ("<< i << ") " << "px = " << px << " py = " << py << " pz = " << pz << " e = " << e << std::endl; } TLorentzVector MCPtotalMass(px,py,pz,e); float MCPe = MCPtotalMass.E(); float MCPpt = MCPtotalMass.Pt(); // std::cout << " MCPe = " << MCPe << " MCPpt = " << MCPpt << std::endl; float evis = totalMass.E(); float pt = totalMass.Pt(); float cosMis = 0; if(totalMass.P()>0)cosMis=fabs(totalMass.Pz())/totalMass.P(); _fEvis->Fill(evis,_wgt); _fCosMis->Fill(cosMis,_wgt); _fPt->Fill(pt,_wgt); _fPFA->Fill(50.,1.); //_fPFA->Fill(totalMass.M(),_wgt); _fPFAnu->Fill(evis+_mcNuEnergy,_wgt); _fNtracks->Fill((float)ntracks,_wgt); // apply preselection bool presel = true; if(ntracks<4)presel=false; if(evis<100)presel=false; if(!presel)return; //_dataStore->passedAnalysisCut("ZHmumuX","Preselection"); _fEvisP->Fill(evis,_wgt); _fCosMisP->Fill(cosMis,_wgt); _fPtP->Fill(pt,_wgt); _fNtracksP->Fill((float)ntracks,_wgt); // muon analysis if(Zmm!=NULL){ float dmz = Zmm->getMass() - 91.2; TLorentzVector pmumu(Zmm->getMomentum()[0],Zmm->getMomentum()[1],Zmm->getMomentum()[2],Zmm->getEnergy()); TLorentzVector precoil = TLorentzVector(0.,0.,0., 250.)-pmumu; TLorentzVector gamma(0.,0.,0.,0.); TLorentzVector Zpfos(0.,0.,0.,0.); if(fabs(dmz)<50){ _fmzmm->Fill(pmumu.M(),_wgt); _fmhmm->Fill(precoil.M(),_wgt); _fmzmhmm->Fill(pmumu.M(),precoil.M(),_wgt); //_dataStore->passedAnalysisCut("ZHmumuX","mumu"); // for comparison look at Z without FSR ReconstructedParticleVec zDecay = Zmm->getParticles(); for(unsigned int i=0;igetType()==22)gamma+= TLorentzVector(zDecay[i]->getMomentum()[0], zDecay[i]->getMomentum()[1],zDecay[i]->getMomentum()[2],zDecay[i]->getEnergy() ); Zpfos += TLorentzVector(zDecay[i]->getMomentum()[0], zDecay[i]->getMomentum()[1],zDecay[i]->getMomentum()[2],zDecay[i]->getEnergy() ); } TLorentzVector recoH = totalMass - Zpfos; // if(_printing>0)std::cout << " Reco Z->mumu : " << pmumu.M() << " " << precoil.M() << " " << recoH.M() << std::endl; TLorentzVector pmumug = pmumu -gamma; TLorentzVector precoilg = TLorentzVector(0.,0.,0., 250.)-pmumug; _fmzmmraw->Fill(pmumug.M(),_wgt); _fmhmmraw->Fill(precoilg.M(),_wgt); // selection if(fabs(dmz)<10){ //_dataStore->passedAnalysisCut("ZHmumuX","Z->mumu(tight)"); bool passed = false; if(ntracks>8){ if(cosMis<0.98){ passed = true; if(precoil.M()>110){ _fEvisSa->Fill(evis,_wgt); _fCosMisSa->Fill(cosMis,_wgt); _fPtSa->Fill(pt,_wgt); _fNtracksSa->Fill((float)ntracks,_wgt); //_dataStore->passedAnalysisCut("ZHmumuX","ZH->mumuqq"); } } }else{ if(cosMis<0.98 && evis < 240.){ passed = true; if(precoil.M()>110){ _fEvisSb->Fill(evis,_wgt); _fCosMisSb->Fill(cosMis,_wgt); _fPtSb->Fill(pt,_wgt); _fNtracksSb->Fill((float)ntracks,_wgt); //_dataStore->passedAnalysisCut("ZHmumuX","ZH->mumull"); } } } if(passed){ //_dataStore->passedAnalysisCut("ZHmumuX","Selected"); _fmzmmS->Fill(pmumu.M(),_wgt); _fmhmmS->Fill(precoil.M(),1.); } } } } // eeX analysis if(Zee==NULL && _signal){ // if(_printing>0)std::cout << " ZFinder failed " << std::endl; } if(Zee!=NULL){ float dmz = Zee->getMass() - 91.2; TLorentzVector pee(Zee->getMomentum()[0],Zee->getMomentum()[1],Zee->getMomentum()[2],Zee->getEnergy()); TLorentzVector precoil = TLorentzVector(0.,0.,0., 250.)-pee; TLorentzVector gamma(0.,0.,0.,0.); TLorentzVector Zpfos(0.,0.,0.,0.); // if(_printing>0)std::cout << " Reco Z->ee : " << pee.M() << " " << precoil.M() << " " << Zee->getMass() << std::endl; //if(fabs(dmz)<50){ //something weird here if(1){ //_dataStore->passedAnalysisCut("ZHeeX","ee"); _fmzee->Fill(pee.M(),_wgt); // streamlog_out(MESSAGE) << "MarlinProcessor::AnalyseEvent Fill histo "<Fill(precoil.M(),_wgt); _fmzmhee->Fill(pee.M(),precoil.M(),_wgt); if(_good){ _fmzeeg->Fill(pee.M(),_wgt); _fmheeg->Fill(precoil.M(),_wgt); }else{ _fmzeeb->Fill(pee.M(),_wgt); _fmheeb->Fill(precoil.M(),_wgt); } // for comparison look at Z without FSR ReconstructedParticleVec zDecay = Zee->getParticles(); for(unsigned int i=0;igetType()==22)gamma+= TLorentzVector(zDecay[i]->getMomentum()[0], zDecay[i]->getMomentum()[1],zDecay[i]->getMomentum()[2],zDecay[i]->getEnergy() ); Zpfos += TLorentzVector(zDecay[i]->getMomentum()[0], zDecay[i]->getMomentum()[1],zDecay[i]->getMomentum()[2],zDecay[i]->getEnergy() ); } TLorentzVector recoH = totalMass - Zpfos; TLorentzVector peeg = pee - gamma; TLorentzVector precoilg = TLorentzVector(0.,0.,0., 250.)-peeg; _fmzeeraw->Fill(peeg.M(),_wgt); _fmheeraw->Fill(precoilg.M(),_wgt); if(_good){ _fmzeegraw->Fill(peeg.M(),_wgt); _fmheegraw->Fill(precoilg.M(),_wgt); }else{ _fmzeebraw->Fill(peeg.M(),_wgt); _fmheebraw->Fill(precoilg.M(),_wgt); } // selection if(fabs(dmz)<10){ bool passed = false; if(ntracks>8){ if(cosMis<0.98){ passed = true; } }else{ if(cosMis<0.98 && evis < 240.){ passed = true; } } if(passed){ //_dataStore->passedAnalysisCut("ZHeeX","Selected"); _fmzeeS->Fill(pee.M(),_wgt); _fmheeS->Fill(precoil.M(),_wgt); if(precoil.M()>110){ _fEvisSee->Fill(evis,_wgt); _fCosMisSee->Fill(cosMis,_wgt); _fPtSee->Fill(pt,_wgt); _fNtracksSee->Fill((float)ntracks,_wgt); } } } } } } void MarlinProcessor::AnalyseOneJetEvent() { std::vector fourMomenta; unsigned int i=0; // if(_printing>0){ // std::cout << " Reco Jet " << i << " Energy = " << _onejetvec[i]->getEnergy() << " Mass = " << _onejetvec[i]->getMass() << std::endl; // } // create a root four-vector object from reconstructed jet TLorentzVector jet(_onejetvec[i]->getMomentum()[0] ,_onejetvec[i]->getMomentum()[1] ,_onejetvec[i]->getMomentum()[2], _onejetvec[i]->getEnergy() ); // add it to vector of four-momenta fourMomenta.push_back(jet); //print out the the flavour of jet // std::cout << " The OneJet flavour is " << _onejetTruthvec[0]->at(0) << std::endl; // calculate jet invariant mass TLorentzVector onejet = fourMomenta[0]; // if(_printing>0)std::cout << " m12: " << m12 << std::endl; // std::cout << " The OneJet mass and Energy is " << onejet.M() << " " << onejet.E() << std::endl; // histogram it _fMone->Fill(onejet.M(),_wgt); return; } void MarlinProcessor::AnalyseTwoJetEvent() { std::vector fourMomenta; std::vector bjetn; std::vector cjetn; int bjets=0; int cjets=0; for(unsigned int i=0;i<_twojetvec.size();i++){ // if(_printing>0){ // std::cout << " 2 Jet Collection: Reco Jet " << i << " Energy = " << _twojetvec[i]->getEnergy() << std::endl; // } if (i==0){ bjets=0; cjets=0;} // create a root four-vector object from reconstructed jet TLorentzVector jet(_twojetvec[i]->getMomentum()[0] ,_twojetvec[i]->getMomentum()[1] ,_twojetvec[i]->getMomentum()[2], _twojetvec[i]->getEnergy() ); // add it to vector of four-momenta fourMomenta.push_back(jet); //if b or c jet, store jet number if ( _twojetTruthvec[i]->at(0) == 5.0){ bjetn.push_back(i); bjets++; } if ( _twojetTruthvec[i]->at(0) == 4.0){ cjetn.push_back(i); cjets++; } } // calculate b jets mass TLorentzVector bjetv; float mbb = 0.0; //Check for 2 b jets: if(bjetn.size()==2){ //mass int one=bjetn[0], two=bjetn[1]; bjetv = fourMomenta[one] + fourMomenta[two]; mbb=bjetv.M(); // histogram it _fM2bb->Fill(mbb,_wgt); } // calculate c jets mass TLorentzVector cjetv; float mcc = 0.0; //Check for 2 c jets: if(cjetn.size()==2){ //mass int one=cjetn[0], two=cjetn[1]; cjetv = fourMomenta[one] + fourMomenta[two]; mcc=cjetv.M(); // histogram it _fM2cc->Fill(mcc,_wgt); } // calculate di-jet invariant mass TLorentzVector jet12 = fourMomenta[0] + fourMomenta[1]; float m12 = jet12.M(); // if(_printing>0)std::cout << " m12: " << m12 << std::endl; // histogram it _fM12->Fill(m12,_wgt); return; } void MarlinProcessor::AnalyseThreeJetEvent() { std::vector fourMomenta; std::vector bjetn; std::vector cjetn; int bjets=0; int cjets=0; for(unsigned int i=0;i<_threejetvec.size();i++){ // if(_printing>0){ // std::cout << " 2 Jet Collection: Reco Jet " << i << " Energy = " << _threejetvec[i]->getEnergy() << std::endl; // } if (i==0){ bjets=0; cjets=0;} // create a root four-vector object from reconstructed jet TLorentzVector jet(_threejetvec[i]->getMomentum()[0] ,_threejetvec[i]->getMomentum()[1] ,_threejetvec[i]->getMomentum()[2], _threejetvec[i]->getEnergy() ); // add it to vector of four-momenta fourMomenta.push_back(jet); //if b or c jet, store jet number if ( _threejetTruthvec[i]->at(0) == 5.0){ bjetn.push_back(i); bjets++; } if ( _threejetTruthvec[i]->at(0) == 4.0){ cjetn.push_back(i); cjets++; } } // calculate b jets mass TLorentzVector bjetv; float mbb = 0.0; //Check for 2 b jets: if(bjetn.size()==2){ //mass int one=bjetn[0], two=bjetn[1]; bjetv = fourMomenta[one] + fourMomenta[two]; mbb=bjetv.M(); // histogram it _fM3bb->Fill(mbb,_wgt); } // calculate c jets mass TLorentzVector cjetv; float mcc = 0.0; //Check for 2 c jets: if(cjetn.size()==2){ //mass int one=cjetn[0], two=cjetn[1]; cjetv = fourMomenta[one] + fourMomenta[two]; mcc=cjetv.M(); // histogram it _fM3cc->Fill(mcc,_wgt); } return; } void MarlinProcessor::AnalyseFourJetEvent() { std::vector fourMomenta; std::vector bjetn; std::vector cjetn; std::vector FTbjetn; std::vector FTcjetn; int bjets=0; int cjets=0; int FTbjets=0; int FTcjets=0; //Look at MCParticles of the Higgs decay first std::vector MCHb; std::vector MCHbbar; std::vector MCHdecay; int nHb=0, nHbbar=0; for(unsigned int i=0;i<_HbbMC.size();i++){ // create a root four-vector object from reconstructed jet TLorentzVector jet(_HbbMC[i]->getMomentum()[0] ,_HbbMC[i]->getMomentum()[1] ,_HbbMC[i]->getMomentum()[2], _HbbMC[i]->getEnergy() ); MCHdecay.push_back(jet); } for(unsigned int i=0;i<_HbbMC.size();i++){ if(_HbbMC[i]->getPDG()==5){TLorentzVector jet(_HbbMC[i]->getMomentum()[0] ,_HbbMC[i]->getMomentum()[1] ,_HbbMC[i]->getMomentum()[2], _HbbMC[i]->getEnergy() ); MCHb.push_back(jet); nHb++;} if(_HbbMC[i]->getPDG()==-5){TLorentzVector jet(_HbbMC[i]->getMomentum()[0] ,_HbbMC[i]->getMomentum()[1] ,_HbbMC[i]->getMomentum()[2], _HbbMC[i]->getEnergy() ); MCHbbar.push_back(jet);nHbbar++;} } if(nHb==1 && nHbbar==1){ //Combined Vector of e- e+ TLorentzVector p_Hreco = MCHb[0]+MCHbbar[0]; //Fill histograms for this MC H particle _fMCHM->Fill(p_Hreco.M(), _wgt); _fMCHE->Fill(p_Hreco.E(), _wgt); _fMCHEt->Fill(p_Hreco.Et(), _wgt); _fMCHPt->Fill(p_Hreco.Pt(), _wgt); _fMCHEta->Fill(p_Hreco.Eta(), _wgt); } if( nHb==1 && nHbbar==1){ //LOOK AT THE TRUTH TAGGED JETS int nb=0, nbbar=0, bNum = 0, bbarNum =0; for(unsigned int i=0;i<_fourjetvec.size();i++){ //std::cout << " 4 Jet Collection: Reco Jet " << i << " Energy = " << _fourjetvec[i]->getEnergy() << std::endl; // create a root four-vector object from reconstructed jet TLorentzVector jet(_fourjetvec[i]->getMomentum()[0] ,_fourjetvec[i]->getMomentum()[1] ,_fourjetvec[i]->getMomentum()[2], _fourjetvec[i]->getEnergy() ); // add it to vector of four-momenta fourMomenta.push_back(jet); //LOOK AT B-JETS //if b jt, store jet number if ( _fourjetTruthvec[i]->at(0) == 5.0){ //true b jet event //_fFTVNNbb->Fill(_fourjetFTvec[i]->at(0),_wgt); //_fEt4bb->Fill(jet.Et(),_wgt); //_fPR4bb->Fill(jet.Eta(),_wgt); //_fPt4bb->Fill(jet.Pt(),_wgt); if ( _fourjetTruthvec[i]->at(1) == 1){nb++;bNum=i;} if ( _fourjetTruthvec[i]->at(1) == -1){nbbar++;bbarNum=i;} bjetn.push_back(i); bjets++; } } TLorentzVector jet1 = fourMomenta[0]; TLorentzVector jet2 = fourMomenta[1]; TLorentzVector jet3 = fourMomenta[2]; TLorentzVector jet4 = fourMomenta[3]; // calculate jet invariant mass // std::cout << " Jet 1's mass and energy is " << jet1.M() << " " << jet1.E() << std::endl; // std::cout << " Jet 2's mass and energy is " << jet2.M() << " " << jet2.E() << std::endl; // std::cout << " Jet 3's mass and energy is " << jet3.M() << " " << jet3.E() << std::endl; // std::cout << " Jet 4's mass and energy is " << jet4.M() << " " << jet4.E() << std::endl; //print out the the flavour of each jet // std::cout << " Jet 1, 2, 3 and 4's flavour tag is " << _fourjetTruthvec[0]->at(0) << " " << _fourjetTruthvec[1]->at(0) << " " << _fourjetTruthvec[2]->at(0) << " " << _fourjetTruthvec[3]->at(0) << std::endl; // calculate Higgs mass from truth tagged jets //ENERGY SCALING FACTOR float Escaling = 1.0; TLorentzVector bjetv; float mbb = 0.0; float PRbb = 0.0, Ptbb = 0.0; //Check for 2 b jets with one b and one bbar: if(bjetn.size()==2 && nb==1 && nbbar==1){ int one=bjetn[0], two=bjetn[1]; //mass, eta and Pt bjetv = fourMomenta[one]; PRbb = (bjetv.PseudoRapidity())*Escaling; Ptbb = (bjetv.Pt())*Escaling; _fPR4dbb->Fill(PRbb,_wgt); _fPt4dbb->Fill(Ptbb,_wgt); _fE4dbb->Fill(((bjetv.E())*Escaling),_wgt); _fEt4dbb->Fill(((bjetv.Et())*Escaling),_wgt); //if (bjetv.E() >= 25.0 && bjetv.E() <= 50.0){ _fPt4dbb2550->Fill(Ptbb,_wgt); _fEta4dbb2550->Fill(PRbb,_wgt);} //if (bjetv.E() > 50.0 && bjetv.E() < 75.0){_fPt4dbb5075->Fill(Ptbb,_wgt); _fEta4dbb5075->Fill(PRbb,_wgt);} //if (bjetv.E() > 75.0){_fPt4dbb75->Fill(Ptbb,_wgt); _fEta4dbb75->Fill(PRbb,_wgt);} bjetv = fourMomenta[two]; PRbb = (bjetv.PseudoRapidity())*Escaling; Ptbb = (bjetv.Pt())*Escaling; _fPR4dbb->Fill(PRbb,_wgt); _fPt4dbb->Fill(Ptbb,_wgt); _fE4dbb->Fill(((bjetv.E())*Escaling),_wgt); _fEt4dbb->Fill(((bjetv.Et())*Escaling),_wgt); //if (bjetv.E() >= 25.0 && bjetv.E() <= 50.0){ _fPt4dbb2550->Fill(Ptbb,_wgt); _fEta4dbb2550->Fill(PRbb,_wgt);} //if (bjetv.E() > 50.0 && bjetv.E() < 75.0){_fPt4dbb5075->Fill(Ptbb,_wgt); _fEta4dbb5075->Fill(PRbb,_wgt);} //if (bjetv.E() > 75.0){_fPt4dbb75->Fill(Ptbb,_wgt); _fEta4dbb75->Fill(PRbb,_wgt);} bjetv = fourMomenta[one] + fourMomenta[two]; mbb=(bjetv.M())*Escaling; // histogram it _fM4bb->Fill(mbb,_wgt); _fE4bb->Fill(((bjetv.E())*Escaling),_wgt); _fEt4bb->Fill(((bjetv.Et())*Escaling),_wgt); } //Fill energy resolution histograms //Check for 2 b jets: if(bjetn.size()==2 && nb==1 && nbbar==1){//only consider H->b bbar events which match with recojets of b bbar if(nHb==1 && nb==1){//b quarks bjetv = fourMomenta[bNum]; //std::cout << " B QUARK bjetv.E() MCHb[0].E() " << bjetv.E() << " " << MCHb[0].E() << std::endl; _fEresHbb ->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); //H resolution _fEresHb ->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); if( fabs( MCHb[0].Eta() ) <=0.7 ) _fEresHbb1->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); if( fabs( MCHb[0].Eta() ) >0.7 && fabs( MCHb[0].Eta() ) <=1.4 ) _fEresHbb2->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); if( fabs( MCHb[0].Eta() ) >1.4 ) _fEresHbb3->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); } if(nHbbar==1 && nbbar==1){//b bar quarks bjetv = fourMomenta[bbarNum]; //std::cout << " B QUARK bjetv.E() MCHb[0].E() " << bjetv.E() << " " << MCHbbar[0].E() << std::endl; _fEresHbb ->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt);//H resolution _fEresHbbar ->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); if( fabs( MCHbbar[0].Eta() ) <=0.7 ) _fEresHbb1->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); if( fabs( MCHbbar[0].Eta() ) >0.7 && fabs( MCHbbar[0].Eta() ) <=1.4 ) _fEresHbb2->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); if( fabs( MCHbbar[0].Eta() ) >1.4 ) _fEresHbb3->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); } } } //Analyse the constituent PFOs for jet _fourjetvec[i] : // for(unsigned int i=0;i<_fourjetvec.size();i++){ // ReconstructedParticleVec constituents = _fourjetvec[i]->getParticles(); // std::cout<< "For jet " << i << " constituent particles size is " << constituents.size() <getCharge()<<" "<getEnergy()<getParticles(); // std::cout<< "For jet " << i << " constituent particles size is " << constituents.size() <getCharge()<<" "<getEnergy()<at(0) == 5.0){ if ( _fourjetTruthvec[i]->at(1) == 1){bJNum=i;} if ( _fourjetTruthvec[i]->at(1) == -1){bbarJNum=i;} } } //Analyse the constituent PFOs of the jets float NChargeb=0.0, NNeutb=0.0, NChargebbar=0.0, NNeutbbar=0.0; float EChargeb=0.0, ENeutb=0.0, EChargebbar=0.0, ENeutbbar=0.0; float NTotalPFOsb=0.0, ETotalPFOsb=0.0, NTotalPFOsbbar=0.0, ETotalPFOsbbar=0.0; //b jet ReconstructedParticleVec bconstituents = _fourjetvec[bJNum]->getParticles(); // std::cout<< "For jet " << bJNum << " bconstituents particles size is " << bconstituents.size() <getCharge() != 0.0){NChargeb+= 1.0; EChargeb+=bconstituents[ii]->getEnergy();} if(bconstituents[ii]->getCharge() == 0.0){NNeutb+= 1.0; ENeutb+=bconstituents[ii]->getEnergy();} // std::cout<< bconstituents[ii]->getCharge()<<" "<getEnergy()<getEnergy(); } //bbar jet ReconstructedParticleVec bbarconstituents = _fourjetvec[bbarJNum]->getParticles(); // std::cout<< "For jet " << bbarJNum << " bbarconstituents particles size is " << bbarconstituents.size() <getCharge() != 0.0){NChargebbar+= 1.0; EChargebbar+=bbarconstituents[ii]->getEnergy();} if(bbarconstituents[ii]->getCharge() == 0.0){NNeutbbar+= 1.0; ENeutbbar+=bbarconstituents[ii]->getEnergy();} // std::cout<< bbarconstituents[ii]->getCharge()<<" "<getEnergy()<getEnergy(); } // std::cout<< "Charged PFO & neutral PFO number and energy for b jet " << NChargeb << " " << NNeutb << " " << EChargeb << " " << ENeutb <Fill((EChargeb/ETotalPFOsb), _wgt); _fJetResCEbbar->Fill((EChargebbar/ETotalPFOsbbar), _wgt); _fJetResCNb->Fill((NChargeb/NTotalPFOsb), _wgt); _fJetResCNbbar->Fill((NChargebbar/NTotalPFOsbbar), _wgt); _fJetResNEb->Fill((ENeutb/ETotalPFOsb), _wgt); _fJetResNEbbar->Fill((ENeutbbar/ETotalPFOsbbar), _wgt); _fJetResNNb->Fill((NNeutb/NTotalPFOsb), _wgt); _fJetResNCNbbar->Fill((NNeutbbar/NTotalPFOsbbar), _wgt); return; } void MarlinProcessor::AnalysePandoraPFANewPFOs() { //Was created for electron gun simulations and is wrong!!! std::vector PANDORAe; //PANDORA PART for(unsigned int i=0;i<_Pele.size();i++){ // create a root four-vector object from reconstructed jet TLorentzVector jet(_Pele[i]->getMomentum()[0] ,_Pele[i]->getMomentum()[1] ,_Pele[i]->getMomentum()[2], _Pele[i]->getEnergy() ); PANDORAe.push_back(jet); } // std::cout << " PandoraPFANewPFOs ELECTRON Mass, Energy, Pt and PseudoRapidity " << PANDORAe[0].M() << " " << PANDORAe[0].E() << " " << PANDORAe[0].Pt() << " " << PANDORAe[0].PseudoRapidity() << std::endl; //MCPARTICLE PART // std::cout<<"MCParticles size wasn't zero, size = " << _MCele.size() << std::endl; std::vector MCPARTICLEe; std::vector MCelectron; int pcounter=0; float partElectronEnergy = 100.0; float closeness = 100.0; float partEnergyHighest = 0.0; for(unsigned int i=0;i<_MCele.size();i++){ // create a root four-vector object from reconstructed jet TLorentzVector jet(_MCele[i]->getMomentum()[0] ,_MCele[i]->getMomentum()[1] ,_MCele[i]->getMomentum()[2], _MCele[i]->getEnergy() ); // add it to vector of four-momenta MCelectron.push_back(jet); if(fabs(jet.E()-partElectronEnergy)Fill(0.1,_wgt); if( fabs( PANDORAe[0].Eta() ) >0.7 && fabs( PANDORAe[0].Eta() ) <=1.4 ) _dEtaEresRange->Fill(0.8,_wgt); if( fabs( PANDORAe[0].Eta() ) >1.4 ) _dEtaEresRange->Fill(1.3,_wgt); return; } void MarlinProcessor::AnalyseEmEp() { // std::cout << "In AnalyseEmEp()" << std::endl; //Z reconstruction /////////MC PART std::vector MCem; std::vector MCep; int MCNem=0, MCNep=0, MCemNum = 0, MCepNum =0; for(unsigned int i=0;i<_Z0MC.size();i++){//Get the MC e- and e+ if(_Z0MC[i]->getPDG()==11){ TLorentzVector jet(_Z0MC[i]->getMomentum()[0] ,_Z0MC[i]->getMomentum()[1] ,_Z0MC[i]->getMomentum()[2], _Z0MC[i]->getEnergy() ); MCem.push_back(jet); MCNem++; MCemNum=i; } if(_Z0MC[i]->getPDG()==-11){ TLorentzVector jet(_Z0MC[i]->getMomentum()[0] ,_Z0MC[i]->getMomentum()[1] ,_Z0MC[i]->getMomentum()[2], _Z0MC[i]->getEnergy() ); MCep.push_back(jet); MCNep++; MCepNum=i; } } if(MCNem==1 && MCNep==1){ //Combined Vector of e- e+ TLorentzVector p_Z0reco = MCem[0]+MCep[0]; //Fill histograms for this MC Z particle _fMCZ0M->Fill(p_Z0reco.M(), _wgt); _fMCZ0E->Fill(p_Z0reco.E(), _wgt); _fMCZ0Et->Fill(p_Z0reco.Et(), _wgt); _fMCZ0Pt->Fill(p_Z0reco.Pt(), _wgt); _fMCZ0Eta->Fill(p_Z0reco.Eta(), _wgt); } /////////PANDORA PART std::vector PANDORAem; std::vector PANDORAep; float emcloseness=1000.0, epcloseness=1000.0;//very high value as starting point int PANemNum=0, PANepNum=0; int testem=0, testep=0; if(_Pele.size()>0){//Make sure PANDORA e- or e+ exist for(unsigned int i=0;i<_Pele.size();i++){ //Look at PANDORA e- and e+ to compare with the MC e- and e+ if(_Pele[i]->getType()==11){ if(fabs(_Pele[i]->getEnergy()-MCem[0].E())getEnergy()-MCem[0].E()); PANemNum=i; testem=1; } } if(_Pele[i]->getType()==-11){ if(fabs(_Pele[i]->getEnergy()-MCep[0].E())getEnergy()-MCep[0].E()); PANepNum=i; testep=1; } } } } if(testem==1 && testep==1){//if both e- and e+ were reconstructed TLorentzVector jetem(_Pele[PANemNum]->getMomentum()[0] ,_Pele[PANemNum]->getMomentum()[1] ,_Pele[PANemNum]->getMomentum()[2], _Pele[PANemNum]->getEnergy() ); PANDORAem.push_back(jetem); TLorentzVector jetep(_Pele[PANepNum]->getMomentum()[0] ,_Pele[PANepNum]->getMomentum()[1] ,_Pele[PANepNum]->getMomentum()[2], _Pele[PANepNum]->getEnergy() ); PANDORAep.push_back(jetep); TLorentzVector p_Z0recoPANDORA = PANDORAem[0]+PANDORAep[0]; //Fill histograms for reconstructed Z particle _fPANZ0M->Fill(p_Z0recoPANDORA.M(), _wgt); _fPANZ0E->Fill(p_Z0recoPANDORA.E(), _wgt); _fPANZ0Et->Fill(p_Z0recoPANDORA.Et(), _wgt); _fPANZ0Pt->Fill(p_Z0recoPANDORA.Pt(), _wgt); _fPANZ0Eta->Fill(p_Z0recoPANDORA.Eta(), _wgt); } //ENERGY RESOLUTION HISTOGRAMS - requires both e- and e+ PANDORA particles. //Fill energy resolution histograms if( MCNem==1 && testem==1 && MCNep==1 && testep==1){//electrons and positrons check // std::cout << " PANDORAem[0].E() MCem[0].E() " << PANDORAem[0].E() << " " << MCem[0].E() << std::endl; _fEresPANMCee ->Fill((PANDORAem[0].E()/MCem[0].E()), _wgt); _fEresPANMCem ->Fill((PANDORAem[0].E()/MCem[0].E()), _wgt); _fEresPANEem ->Fill(PANDORAem[0].E(), _wgt); _fEresMCEem ->Fill(MCem[0].E(), _wgt); if( fabs( MCem[0].Eta() ) <=0.7 ) _fEresPANMCeeEta1->Fill((PANDORAem[0].E()/MCem[0].E()),_wgt); if( fabs( MCem[0].Eta() ) >0.7 && fabs( MCem[0].Eta() ) <=1.4 ) _fEresPANMCeeEta2->Fill((PANDORAem[0].E()/MCem[0].E()),_wgt); if( fabs( MCem[0].Eta() ) >1.4 ) _fEresPANMCeeEta3->Fill((PANDORAem[0].E()/MCem[0].E()),_wgt); _fEresPANMCee ->Fill((PANDORAep[0].E()/MCep[0].E()), _wgt); _fEresPANMCep ->Fill((PANDORAep[0].E()/MCep[0].E()), _wgt); _fEresPANEep ->Fill(PANDORAep[0].E(), _wgt); _fEresMCEep ->Fill(MCep[0].E(), _wgt); // std::cout << " PANDORAep[0].E() MCep[0].E() " << PANDORAep[0].E() << " " << MCep[0].E() << std::endl; if( fabs( MCep[0].Eta() ) <=0.7 ) _fEresPANMCeeEta1->Fill((PANDORAep[0].E()/MCep[0].E()),_wgt); if( fabs( MCep[0].Eta() ) >0.7 && fabs( MCep[0].Eta() ) <=1.4 ) _fEresPANMCeeEta2->Fill((PANDORAep[0].E()/MCep[0].E()),_wgt); if( fabs( MCep[0].Eta() ) >1.4 ) _fEresPANMCeeEta3->Fill((PANDORAep[0].E()/MCep[0].E()),_wgt); } else{ // std::cout << " PANDORA DIDN'T FIND BOTH 1 ELECTRON AND 1 POSITRON FOR RECONSTRUCTION SO HAVE SKIPPED THIS EVENT" << std::endl; } //Electron EM and HAD energies e- in the Z reconstruction if(MCNem==1 && testem==1 && MCNep==1 && testep==1){ ClusterVec cPANem = _Pele[PANemNum]->getClusters(); //std::cout<<"CLUSTER VEC SIZE OF PANem : " << cPANem.size() << std::endl; float ecalonly = 0.; float lcalonly = 0.; float bcalonly = 0.; float hcalonly = 0.; float lhcalonly = 0.; // WANT E(ecal+lcal) versus E(hcal+lhcal)? if(cPANem.size()==1){ // parameter ClusterSubdetectorNames[string]: 0=ecal, 1=hcal, 2=yoke, 3=lcal, 4=lhcal, 5=bcal ecalonly = cPANem[0]->getSubdetectorEnergies()[0]; lcalonly = cPANem[0]->getSubdetectorEnergies()[3]; bcalonly = cPANem[0]->getSubdetectorEnergies()[5]; hcalonly = cPANem[0]->getSubdetectorEnergies()[1]; lhcalonly = cPANem[0]->getSubdetectorEnergies()[4]; _fEcalEem->Fill(ecalonly, _wgt); _fEMEem->Fill((ecalonly+lcalonly+bcalonly), _wgt); _fHcalEem->Fill(hcalonly, _wgt); _fHADEem->Fill((hcalonly+lhcalonly), _wgt); // std::cout<<"IN TEST EM CLUSTER!! ecal, lcal, bcal, hcal and lhcal : " << ecalonly << " " << lcalonly << " " << bcalonly << " " << hcalonly << " " << lhcalonly << std::endl; } } //Electron EM and HAD energies e+ in the Z reconstruction if(MCNem==1 && testem==1 && MCNep==1 && testep==1){ ClusterVec cPANep = _Pele[PANepNum]->getClusters(); //std::cout<<"CLUSTER VEC SIZE OF PANem : " << cPANem.size() << std::endl; float ecalonly = 0.; float lcalonly = 0.; float bcalonly = 0.; float hcalonly = 0.; float lhcalonly = 0.; // WANT E(ecal+lcal) versus E(hcal+lhcal)? if(cPANep.size()==1){ // parameter ClusterSubdetectorNames[string]: 0=ecal, 1=hcal, 2=yoke, 3=lcal, 4=lhcal, 5=bcal ecalonly = cPANep[0]->getSubdetectorEnergies()[0]; lcalonly = cPANep[0]->getSubdetectorEnergies()[3]; bcalonly = cPANep[0]->getSubdetectorEnergies()[5]; hcalonly = cPANep[0]->getSubdetectorEnergies()[1]; lhcalonly = cPANep[0]->getSubdetectorEnergies()[4]; _fEcalEep->Fill(ecalonly, _wgt); _fEMEep->Fill((ecalonly+lcalonly+bcalonly), _wgt); _fHcalEep->Fill(hcalonly, _wgt); _fHADEep->Fill((hcalonly+lhcalonly), _wgt); // std::cout<<"IN TEST EM CLUSTER!! ecal, lcal, bcal, hcal and lhcal : " << ecalonly << " " << lcalonly << " " << bcalonly << " " << hcalonly << " " << lhcalonly << std::endl; } } return; } void MarlinProcessor::AnalyseRecoParticlesEvent() { //ENERGY SCALING FACTOR float Escaling = 1.0; // float Escaling = 2.040816327; //MC PART for Z daughters of e- and e+ std::vector MCem; std::vector MCep; int MCNem=0, MCNep=0, MCemNum = 0, MCepNum =0; for(unsigned int i=0;i<_Z0MC.size();i++){//Get the MC e- and e+ if(_Z0MC[i]->getPDG()==11){ TLorentzVector jet(_Z0MC[i]->getMomentum()[0] ,_Z0MC[i]->getMomentum()[1] ,_Z0MC[i]->getMomentum()[2], _Z0MC[i]->getEnergy() ); MCem.push_back(jet); MCNem++; MCemNum=i; } if(_Z0MC[i]->getPDG()==-11){ TLorentzVector jet(_Z0MC[i]->getMomentum()[0] ,_Z0MC[i]->getMomentum()[1] ,_Z0MC[i]->getMomentum()[2], _Z0MC[i]->getEnergy() ); MCep.push_back(jet); MCNep++; MCepNum=i; } } TLorentzVector p_Z0recoMC = MCem[0]+MCep[0]; // std::cout<<"MC Z0 PARTICLE M, E, Eta, Phi are: " << p_Z0recoMC.M() << " " << p_Z0recoMC.E() << " " << p_Z0recoMC.Eta() << " " << p_Z0recoMC.Phi()<< std::endl; std::cout<<"MCem and MCep M, E, Eta, Phi are: " << MCem[0].M() << " " << MCem[0].E() << " " << MCem[0].Eta() << " " << MCem[0].Phi() << " " << MCep[0].M() << " " << MCep[0].E() << " " << MCep[0].Eta() << " " << MCep[0].Phi() << std::endl; //std::cout<<"IN ANALYSE RECO PARTICLES EVENT!!!!" << std::endl; // std::cout<<"_RecoParvec.size() = " << _RecoParvec.size() << std::endl; std::vector RecoPar; std::vector RecoParTest; double emDRcloseness=100.0, epDRcloseness=100.0;//very high value as starting point int WOLFemNum=0, WOLFepNum=0; int testem=0, testep=0; float NumCPFO=0.0, NumNPFO=0.0; double deltaREm=0.0, deltaREp=0.0; float emcloseness=1000.0, epcloseness=1000.0;//very high value as starting point if(_RecoParvec.size()>1){//Try to make sure WOLF e- and e+ exist for(unsigned int i=0;i<_RecoParvec.size();i++){ //Look at WOLF e- and e+ to compare with the MC e- and e+ TLorentzVector jet(_RecoParvec[i]->getMomentum()[0] ,_RecoParvec[i]->getMomentum()[1] ,_RecoParvec[i]->getMomentum()[2], _RecoParvec[i]->getEnergy() ); RecoParTest.push_back(jet); // std::cout<<"WOLF particle " << i << " E, Eta, Phi and type of " << RecoParTest[i].E() << " " << RecoParTest[i].Eta() << " " << RecoParTest[i].Phi() << " " << _RecoParvec[i]->getType() << std::endl; deltaREm=MCem[0].DeltaR(RecoParTest[i]); deltaREp=MCep[0].DeltaR(RecoParTest[i]); // std::cout<<"deltaREm and deltaREp are: " << deltaREm << " " << deltaREp << std::endl; //std::cout<<"jet M and E are: " << RecoParTest[i].M() << " " << RecoParTest[i].E() << std::endl; // if(fabs(_RecoParvec[i]->getEnergy()-MCem[0].E())getEnergy()-MCem[0].E()); // WOLFemNum=i; // testem=1; // } // } // for(unsigned int i=0;i<_RecoParvec.size();i++){ //NOW FIND THE SECOND CLOSEST // if(i!=WOLFemNum){ // if(fabs(_RecoParvec[i]->getEnergy()-MCep[0].E())getEnergy()-MCep[0].E()); // WOLFepNum=i; // testep=1; // } // } // } //OLD DeltaR WAY //DeltaR = ( (Mcem.Eta()-jet.Eta())*(MCem.Eta()-jet.Eta()) )+( (Mcem.Phi()-jet.Phi())*(MCem.Phi()-jet.Phi()) ) // deltaREm=MCem[0].DeltaR(jet[0]); // deltaREp=MCep[0].DeltaR(jet[0]); // std::cout<<"deltaREm and deltaREp are: " << deltaREm << " " << deltaREp << std::endl; //Test for em then ep //if(fabs(deltaREm)getEnergy()-MCem[0].E())getEnergy()-MCem[0].E()); // std::cout<<"IN DELTAREM TEST AGAIN " << std::endl; //emDRcloseness=fabs(deltaREm); WOLFemNum=i; testem=1; } } //if(fabs(deltaREp)getEnergy()-MCep[0].E())getEnergy()-MCep[0].E()); // std::cout<<"IN DELTAREM TEST AGAIN " << std::endl; // std::cout<<"IN DELTAREP TEST AGAIN " << std::endl; epDRcloseness=fabs(deltaREp); WOLFepNum=i; testep=1; } } // if (deltaREm <0.1 || deltaREp <0.1){ // std::cout<<"WOLF particle " << i << " E, Eta, Phi and type of " << RecoParTest[i].E() << " " << RecoParTest[i].Eta() << " " << RecoParTest[i].Phi() << " " << _RecoParvec[i]->getType() << std::endl; // std::cout<<"deltaREm and deltaREp values are " << deltaREm << " " << deltaREp << std::endl; // } } } // std::cout<<"emDRcloseness with particle number is " << emDRcloseness << " " << WOLFemNum << std::endl; // std::cout<<"epDRcloseness with particle number is " << epDRcloseness << " " << WOLFepNum << std::endl; std::vector WOLFem; std::vector WOLFep; // if(testem==1 && testep==1 && WOLFemNum!=WOLFepNum){//if both e- and e+ were reconstructed and not the same WOLF reconstructed particle TLorentzVector jetem(_RecoParvec[WOLFemNum]->getMomentum()[0] ,_RecoParvec[WOLFemNum]->getMomentum()[1] ,_RecoParvec[WOLFemNum]->getMomentum()[2], _RecoParvec[WOLFemNum]->getEnergy() ); WOLFem.push_back(jetem); TLorentzVector jetep(_RecoParvec[WOLFepNum]->getMomentum()[0] ,_RecoParvec[WOLFepNum]->getMomentum()[1] ,_RecoParvec[WOLFepNum]->getMomentum()[2], _RecoParvec[WOLFepNum]->getEnergy() ); WOLFep.push_back(jetep); TLorentzVector p_Z0recoWOLF = WOLFem[0]+WOLFep[0]; //ENERGY RESOLUTION PLOTS _fEresWOLFEem->Fill((Escaling*WOLFem[0].E()), _wgt); _fEresWOLFMCem->Fill(((Escaling*WOLFem[0].E())/(MCem[0].E())), _wgt); _fEresWOLFMCep->Fill(((Escaling*WOLFep[0].E())/(MCep[0].E())), _wgt); if(WOLFem[0].E()<10.0){ _fEresWOLFEem10Eta->Fill((Escaling*WOLFem[0].Eta()), _wgt); // std::cout<<"WOLFem Energy is less than 10!! for particle " << WOLFemNum << std::endl; // std::cout<<"RECO Z0 PARTICLE WOLFem[0] E and Eta() is: " << WOLFem[0].E() << " " << WOLFem[0].Eta() << std::endl; } if(WOLFep[0].E()<10.0){ // std::cout<<"WOLFep Energy is less than 10!! for particle " << WOLFepNum << std::endl; _fEresWOLFEep10Eta->Fill((Escaling*WOLFep[0].Eta()), _wgt); } //Fill histograms for reconstructed Z particle _fWOLFZ0M->Fill((Escaling*p_Z0recoWOLF.M()), _wgt); _fWOLFZ0E->Fill((Escaling*p_Z0recoWOLF.E()), _wgt); _fWOLFZ0Et->Fill((Escaling*p_Z0recoWOLF.Et()), _wgt); _fWOLFZ0Pt->Fill((Escaling*p_Z0recoWOLF.Pt()), _wgt); _fWOLFZ0Eta->Fill((Escaling*p_Z0recoWOLF.Eta()), _wgt); // std::cout<<"RECO Z0 PARTICLE M and E is: " << p_Z0recoWOLF.M() << " " << p_Z0recoWOLF.E() << std::endl; // std::cout<<"FROM em and ep M and E : " << WOLFem[0].M() << " " << WOLFem[0].E() << " " << WOLFep[0].M() << " " << WOLFep[0].E() <getClusters(); //std::cout<<"CLUSTER VEC SIZE OF PANem : " << cPANem.size() << std::endl; float ecalonly = 0.; float lcalonly = 0.; float bcalonly = 0.; float hcalonly = 0.; float lhcalonly = 0.; // WANT E(ecal+lcal) versus E(hcal+lhcal)? if(cPANem.size()==1){ // parameter ClusterSubdetectorNames[string]: 0=ecal, 1=hcal, 2=yoke, 3=lcal, 4=lhcal, 5=bcal ecalonly = cPANem[0]->getSubdetectorEnergies()[0]; lcalonly = cPANem[0]->getSubdetectorEnergies()[3]; bcalonly = cPANem[0]->getSubdetectorEnergies()[5]; hcalonly = cPANem[0]->getSubdetectorEnergies()[1]; lhcalonly = cPANem[0]->getSubdetectorEnergies()[4]; _fEcalWOLFEem->Fill((Escaling*ecalonly), _wgt); _fEMWOLFEem->Fill((Escaling*(ecalonly+lcalonly+bcalonly)), _wgt); _fHcalWOLFEem->Fill((Escaling*hcalonly), _wgt); _fHADWOLFEem->Fill((Escaling*(hcalonly+lhcalonly)), _wgt); // std::cout<<"IN TEST EM CLUSTER!! ecal, lcal, bcal, hcal and lhcal : " << ecalonly << " " << lcalonly << " " << bcalonly << " " << hcalonly << " " << lhcalonly << std::endl; } } //Electron EM and HAD energies e+ in the Z reconstruction if(testem==1 && testep==1){ ClusterVec cPANep = _RecoParvec[WOLFepNum]->getClusters(); //std::cout<<"CLUSTER VEC SIZE OF PANem : " << cPANem.size() << std::endl; float ecalonly = 0.; float lcalonly = 0.; float bcalonly = 0.; float hcalonly = 0.; float lhcalonly = 0.; // WANT E(ecal+lcal) versus E(hcal+lhcal)? if(cPANep.size()==1){ // parameter ClusterSubdetectorNames[string]: 0=ecal, 1=hcal, 2=yoke, 3=lcal, 4=lhcal, 5=bcal ecalonly = cPANep[0]->getSubdetectorEnergies()[0]; lcalonly = cPANep[0]->getSubdetectorEnergies()[3]; bcalonly = cPANep[0]->getSubdetectorEnergies()[5]; hcalonly = cPANep[0]->getSubdetectorEnergies()[1]; lhcalonly = cPANep[0]->getSubdetectorEnergies()[4]; _fEcalWOLFEep->Fill((Escaling*ecalonly), _wgt); _fEMWOLFEep->Fill((Escaling*(ecalonly+lcalonly+bcalonly)), _wgt); _fHcalWOLFEep->Fill((Escaling*hcalonly), _wgt); _fHADWOLFEep->Fill((Escaling*(hcalonly+lhcalonly)), _wgt); // std::cout<<"IN TEST EM CLUSTER!! ecal, lcal, bcal, hcal and lhcal : " << ecalonly << " " << lcalonly << " " << bcalonly << " " << hcalonly << " " << lhcalonly << std::endl; } } return; } void MarlinProcessor::AnalyseWOLFFourJetEvent() { //ENERGY SCALING FACTOR float Escaling = 1.0; // float Escaling = 1.724137931; std::vector fourMomenta; std::vector bjetn; std::vector cjetn; std::vector FTbjetn; std::vector FTcjetn; int bjets=0; int cjets=0; int FTbjets=0; int FTcjets=0; //Look at MCParticles of the Higgs decay first std::vector MCHb; std::vector MCHbbar; std::vector MCHdecay; int nHb=0, nHbbar=0; for(unsigned int i=0;i<_HbbMC.size();i++){ // create a root four-vector object from reconstructed jet TLorentzVector jet(_HbbMC[i]->getMomentum()[0] ,_HbbMC[i]->getMomentum()[1] ,_HbbMC[i]->getMomentum()[2], _HbbMC[i]->getEnergy() ); MCHdecay.push_back(jet); } for(unsigned int i=0;i<_HbbMC.size();i++){ if(_HbbMC[i]->getPDG()==5){TLorentzVector jet(_HbbMC[i]->getMomentum()[0] ,_HbbMC[i]->getMomentum()[1] ,_HbbMC[i]->getMomentum()[2], _HbbMC[i]->getEnergy() ); MCHb.push_back(jet); nHb++;} if(_HbbMC[i]->getPDG()==-5){TLorentzVector jet(_HbbMC[i]->getMomentum()[0] ,_HbbMC[i]->getMomentum()[1] ,_HbbMC[i]->getMomentum()[2], _HbbMC[i]->getEnergy() ); MCHbbar.push_back(jet);nHbbar++;} } if(nHb==1 && nHbbar==1){ //Combined Vector of e- e+ TLorentzVector p_Hreco = MCHb[0]+MCHbbar[0]; //Fill histograms for this MC H particle _fMCHM->Fill(p_Hreco.M(), _wgt); _fMCHE->Fill(p_Hreco.E(), _wgt); _fMCHEt->Fill(p_Hreco.Et(), _wgt); _fMCHPt->Fill(p_Hreco.Pt(), _wgt); _fMCHEta->Fill(p_Hreco.Eta(), _wgt); } if( nHb==1 && nHbbar==1){ //LOOK AT THE TRUTH TAGGED JETS int nb=0, nbbar=0, bNum = 0, bbarNum =0; for(unsigned int i=0;i<_fourjetvec.size();i++){ std::cout << " 4 Jet Collection: Reco Jet " << i << " Energy = " << _fourjetvec[i]->getEnergy() << std::endl; // create a root four-vector object from reconstructed jet TLorentzVector jet(_fourjetvec[i]->getMomentum()[0] ,_fourjetvec[i]->getMomentum()[1] ,_fourjetvec[i]->getMomentum()[2], _fourjetvec[i]->getEnergy() ); // add it to vector of four-momenta fourMomenta.push_back(jet); //LOOK AT B-JETS //if b jt, store jet number if ( _fourjetTruthvec[i]->at(0) == 5.0){ //true b jet event //_fFTVNNbb->Fill(_fourjetFTvec[i]->at(0),_wgt); //_fEt4bb->Fill(jet.Et(),_wgt); //_fPR4bb->Fill(jet.Eta(),_wgt); //_fPt4bb->Fill(jet.Pt(),_wgt); if ( _fourjetTruthvec[i]->at(1) == 1){nb++;bNum=i;} if ( _fourjetTruthvec[i]->at(1) == -1){nbbar++;bbarNum=i;} bjetn.push_back(i); bjets++; } } TLorentzVector jet1 = fourMomenta[0]; TLorentzVector jet2 = fourMomenta[1]; TLorentzVector jet3 = fourMomenta[2]; TLorentzVector jet4 = fourMomenta[3]; // calculate jet invariant mass // std::cout << " Jet 1's mass and energy is " << jet1.M() << " " << jet1.E() << std::endl; // std::cout << " Jet 2's mass and energy is " << jet2.M() << " " << jet2.E() << std::endl; // std::cout << " Jet 3's mass and energy is " << jet3.M() << " " << jet3.E() << std::endl; // std::cout << " Jet 4's mass and energy is " << jet4.M() << " " << jet4.E() << std::endl; //print out the the flavour of each jet // std::cout << " Jet 1, 2, 3 and 4's flavour tag is " << _fourjetTruthvec[0]->at(0) << " " << _fourjetTruthvec[1]->at(0) << " " << _fourjetTruthvec[2]->at(0) << " " << _fourjetTruthvec[3]->at(0) << std::endl; // calculate Higgs mass from truth tagged jets TLorentzVector bjetv; float mbb = 0.0; float PRbb = 0.0, Ptbb = 0.0; //Check for 2 b jets with one b and one bbar: if(bjetn.size()==2 && nb==1 && nbbar==1){ int one=bjetn[0], two=bjetn[1]; //mass, eta and Pt bjetv = fourMomenta[one]; PRbb = (bjetv.PseudoRapidity())*Escaling; Ptbb = (bjetv.Pt())*Escaling; _fPR4dbb->Fill(PRbb,_wgt); _fPt4dbb->Fill(Ptbb,_wgt); _fE4dbb->Fill(((bjetv.E())*Escaling),_wgt); _fEt4dbb->Fill(((bjetv.Et())*Escaling),_wgt); //if (bjetv.E() >= 25.0 && bjetv.E() <= 50.0){ _fPt4dbb2550->Fill(Ptbb,_wgt); _fEta4dbb2550->Fill(PRbb,_wgt);} //if (bjetv.E() > 50.0 && bjetv.E() < 75.0){_fPt4dbb5075->Fill(Ptbb,_wgt); _fEta4dbb5075->Fill(PRbb,_wgt);} //if (bjetv.E() > 75.0){_fPt4dbb75->Fill(Ptbb,_wgt); _fEta4dbb75->Fill(PRbb,_wgt);} bjetv = fourMomenta[two]; PRbb = (bjetv.PseudoRapidity())*Escaling; Ptbb = (bjetv.Pt())*Escaling; _fPR4dbb->Fill(PRbb,_wgt); _fPt4dbb->Fill(Ptbb,_wgt); _fE4dbb->Fill(((bjetv.E())*Escaling),_wgt); _fEt4dbb->Fill(((bjetv.Et())*Escaling),_wgt); //if (bjetv.E() >= 25.0 && bjetv.E() <= 50.0){ _fPt4dbb2550->Fill(Ptbb,_wgt); _fEta4dbb2550->Fill(PRbb,_wgt);} //if (bjetv.E() > 50.0 && bjetv.E() < 75.0){_fPt4dbb5075->Fill(Ptbb,_wgt); _fEta4dbb5075->Fill(PRbb,_wgt);} //if (bjetv.E() > 75.0){_fPt4dbb75->Fill(Ptbb,_wgt); _fEta4dbb75->Fill(PRbb,_wgt);} bjetv = fourMomenta[one] + fourMomenta[two]; mbb=(bjetv.M())*Escaling; std::cout << " mbb is " << mbb << std::endl; // histogram it _fM4bb->Fill(mbb,_wgt); _fE4bb->Fill(((bjetv.E())*Escaling),_wgt); _fEt4bb->Fill(((bjetv.Et())*Escaling),_wgt); } //Fill energy resolution histograms //Check for 2 b jets: if(bjetn.size()==2 && nb==1 && nbbar==1){//only consider H->b bbar events which match with recojets of b bbar if(nHb==1 && nb==1){//b quarks bjetv = fourMomenta[bNum]; //std::cout << " B QUARK bjetv.E() MCHb[0].E() " << bjetv.E() << " " << MCHb[0].E() << std::endl; _fEresHbb ->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); //H resolution _fEresHb ->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); if( fabs( MCHb[0].Eta() ) <=0.7 ) _fEresHbb1->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); if( fabs( MCHb[0].Eta() ) >0.7 && fabs( MCHb[0].Eta() ) <=1.4 ) _fEresHbb2->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); if( fabs( MCHb[0].Eta() ) >1.4 ) _fEresHbb3->Fill((((bjetv.E())*Escaling)/MCHb[0].E()), _wgt); } if(nHbbar==1 && nbbar==1){//b bar quarks bjetv = fourMomenta[bbarNum]; //std::cout << " B QUARK bjetv.E() MCHb[0].E() " << bjetv.E() << " " << MCHbbar[0].E() << std::endl; _fEresHbb ->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt);//H resolution _fEresHbbar ->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); if( fabs( MCHbbar[0].Eta() ) <=0.7 ) _fEresHbb1->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); if( fabs( MCHbbar[0].Eta() ) >0.7 && fabs( MCHbbar[0].Eta() ) <=1.4 ) _fEresHbb2->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); if( fabs( MCHbbar[0].Eta() ) >1.4 ) _fEresHbb3->Fill((((bjetv.E())*Escaling)/MCHbbar[0].E()), _wgt); } } } //Analyse the constituent PFOs for jet _fourjetvec[i] : // for(unsigned int i=0;i<_fourjetvec.size();i++){ // ReconstructedParticleVec constituents = _fourjetvec[i]->getParticles(); // std::cout<< "For jet " << i << " constituent particles size is " << constituents.size() <getCharge()<<" "<getEnergy()<