From ad0b46d6b9dffa362505b778346773a285a64aa4 Mon Sep 17 00:00:00 2001 From: Juliette Alimena Date: Fri, 14 Oct 2022 10:12:12 +0200 Subject: [PATCH 01/25] runFinal improvements --- config/FCCAnalysisRun.py | 120 +++++++++++++++++- .../bsm/LLPs/DisplacedHNL/analysis_final.py | 20 ++- 2 files changed, 135 insertions(+), 5 deletions(-) diff --git a/config/FCCAnalysisRun.py b/config/FCCAnalysisRun.py index d6dbe8f8d3..1e75ca7319 100644 --- a/config/FCCAnalysisRun.py +++ b/config/FCCAnalysisRun.py @@ -123,6 +123,30 @@ def getElement(rdfModule, element, isFinal=False): return {} else: print('The option <{}> is not available in presel analysis'.format(element)) + elif element=='doScale': + if isFinal: + print('The variable <{}> is optional in your analysis_final.py file return empty dictionary'.format(element)) + return {} + else: print('The option <{}> is not available in presel analysis'.format(element)) + + elif element=='intLumi': + if isFinal: + print('The variable <{}> is optional in your analysis_final.py file, return empty dictionary. However, if you set doScale, then it is mandatory!'.format(element)) + return {} + else: print('The option <{}> is not available in presel analysis'.format(element)) + + elif element=='saveTabular': + if isFinal: + print('The variable <{}> is optional in your analysis_final.py file return empty dictionary'.format(element)) + return {} + else: print('The option <{}> is not available in presel analysis'.format(element)) + + elif element=='cutLabels': + if isFinal: + print('The variable <{}> is optional in your analysis_final.py file return empty dictionary'.format(element)) + return {} + else: print('The option <{}> is not available in presel analysis'.format(element)) + elif element=='geometryFile': print('The variable <{}> is optional in your analysys.py file, return default value empty string'.format(element)) if isFinal: print('The option <{}> is not available in final analysis'.format(element)) @@ -711,6 +735,8 @@ def runFinal(rdfModule): processEvents={} eventsTTree={} processList={} + saveTab=[] + efficiencyList=[] inputDir = getElement(rdfModule,"inputDir", True) if inputDir!="": @@ -723,6 +749,23 @@ def runFinal(rdfModule): if not os.path.exists(outputDir) and outputDir!='': os.system("mkdir -p {}".format(outputDir)) + cutList = getElement(rdfModule,"cutList", True) + length_cuts_names = max([len(cut) for cut in cutList]) + cutLabels = getElement(rdfModule,"cutLabels", True) + + # save a table in a separate tex file + saveTabular = getElement(rdfModule,"saveTabular", True) + if saveTabular: + # option to rewrite the cuts in a better way for the table. otherwise, take them from the cutList + if cutLabels: + cutNames = list(cutLabels.values()) + else: + cutNames = [cut for cut in cutList] + + cutNames.insert(0,' ') + saveTab.append(cutNames) + efficiencyList.append(cutNames) + for pr in getElement(rdfModule,"processList", True): processEvents[pr]=0 eventsTTree[pr]=0 @@ -774,10 +817,9 @@ def runFinal(rdfModule): print('processed events ',processEvents) print('events in ttree ',eventsTTree) - cutList = getElement(rdfModule,"cutList", True) - length_cuts_names = max([len(cut) for cut in cutList]) - histoList = getElement(rdfModule,"histoList", True) + doScale = getElement(rdfModule,"doScale", True) + intLumi = getElement(rdfModule,"intLumi", True) doTree = getElement(rdfModule,"doTree", True) for pr in getElement(rdfModule,"processList", True): @@ -795,6 +837,10 @@ def runFinal(rdfModule): histos_list = [] tdf_list = [] count_list = [] + cuts_list = [] + cuts_list.append(pr) + eff_list=[] + eff_list.append(pr) # Define all histos, snapshots, etc... print ('----> Defining snapshots and histograms') @@ -825,11 +871,45 @@ def runFinal(rdfModule): print ('----> Done') nevents_real += all_events + uncertainty = ROOT.Math.sqrt(all_events) + + if doScale: + all_events = all_events*1.*procDict[pr]["crossSection"]*procDict[pr]["kfactor"]*procDict[pr]["matchingEfficiency"]*intLumi/processEvents[pr] + uncertainty = ROOT.Math.sqrt(all_events)*procDict[pr]["crossSection"]*procDict[pr]["kfactor"]*procDict[pr]["matchingEfficiency"]*intLumi/processEvents[pr] + print(' Printing scaled number of events!!! ') print ('----> Cutflow') print (' {cutname:{width}} : {nevents}'.format(cutname='All events', width=16+length_cuts_names, nevents=all_events)) + + if saveTabular: + cuts_list.append('{nevents:.2e} $\\pm$ {uncertainty:.2e}'.format(nevents=all_events,uncertainty=uncertainty)) # scientific notation - recomended for backgrounds + # cuts_list.append('{nevents:.3f} $\\pm$ {uncertainty:.3f}'.format(nevents=all_events,uncertainty=uncertainty)) # float notation - recomended for signals with few events + eff_list.append(1.) #start with 100% efficiency + for i, cut in enumerate(cutList): - print (' After selection {cutname:{width}} : {nevents}'.format(cutname=cut, width=length_cuts_names, nevents=count_list[i].GetValue())) + neventsThisCut = count_list[i].GetValue() + neventsThisCut_raw = neventsThisCut + uncertainty = ROOT.Math.sqrt(neventsThisCut_raw) + if doScale: + neventsThisCut = neventsThisCut*1.*procDict[pr]["crossSection"]*procDict[pr]["kfactor"]*procDict[pr]["matchingEfficiency"]*intLumi/processEvents[pr] + uncertainty = ROOT.Math.sqrt(neventsThisCut_raw)*procDict[pr]["crossSection"]*procDict[pr]["kfactor"]*procDict[pr]["matchingEfficiency"]*intLumi/processEvents[pr] + print (' After selection {cutname:{width}} : {nevents}'.format(cutname=cut, width=length_cuts_names, nevents=neventsThisCut)) + + # Saving the number of events, uncertainty and efficiency for the output-file + if saveTabular and cut != 'selNone': + if neventsThisCut != 0: + cuts_list.append('{nevents:.2e} $\\pm$ {uncertainty:.2e}'.format(nevents=neventsThisCut,uncertainty=uncertainty)) # scientific notation - recomended for backgrounds + # cuts_list.append('{nevents:.3f} $\\pm$ {uncertainty:.3f}'.format(nevents=neventsThisCut,uncertainty=uncertainty)) # # float notation - recomended for signals with few events + prevNevents = cuts_list[-2].split() + eff_list.append('{eff:.3f}'.format(eff=1.*neventsThisCut/all_events)) + # if number of events is zero, the previous uncertainty is saved instead: + elif '$\\pm$' in cuts_list[-1]: + cut = (cuts_list[-1]).split() + cuts_list.append('$\\leq$ {uncertainty}'.format(uncertainty=cut[2])) + eff_list.append('0.') + else: + cuts_list.append(cuts_list[-1]) + eff_list.append('0.') # And save everything print ('----> Saving outputs') @@ -850,6 +930,38 @@ def runFinal(rdfModule): validfile = testfile(fout_list[i]) if not validfile: continue + if saveTabular and cut != 'selNone': + saveTab.append(cuts_list) + efficiencyList.append(eff_list) + + if saveTabular: + f = open(outputDir+"outputTabular.txt","w") + # Printing the number of events in format of a LaTeX table + print('\\begin{table}[H] \n \\centering \n \\resizebox{\\textwidth}{!}{ \n \\begin{tabular}{|l||',end='',file=f) + print('c|' * (len(cuts_list)-1),end='',file=f) + print('} \hline',file=f) + for i, row in enumerate(saveTab): + print(' ', end='', file=f) + print(*row, sep = ' & ', end='', file=f) + print(' \\\\ ', file=f) + if (i == 0): + print(' \\hline',file=f) + print(' \\hline \n \\end{tabular}} \n \\caption{Caption} \n \\label{tab:my_label} \n\\end{table}', file=f) + + # Efficiency: + print('\n\nEfficiency: ', file=f) + print('\\begin{table}[H] \n \\centering \n \\resizebox{\\textwidth}{!}{ \n \\begin{tabular}{|l||',end='',file=f) + print('c|' * (len(cuts_list)-1),end='',file=f) + print('} \hline',file=f) + for i in range(len(eff_list)): + print(' ', end='', file=f) + v = [row[i] for row in efficiencyList] + print(*v, sep = ' & ', end='', file=f) + print(' \\\\ ', file=f) + if (i == 0): + print(' \\hline',file=f) + print(' \\hline \n \\end{tabular}} \n \\caption{Caption} \n \\label{tab:my_label} \n\\end{table}', file=f) + f.close() elapsed_time = time.time() - start_time print ('==============================SUMMARY==============================') diff --git a/examples/FCCee/bsm/LLPs/DisplacedHNL/analysis_final.py b/examples/FCCee/bsm/LLPs/DisplacedHNL/analysis_final.py index 5aa46d97de..0b11d73245 100644 --- a/examples/FCCee/bsm/LLPs/DisplacedHNL/analysis_final.py +++ b/examples/FCCee/bsm/LLPs/DisplacedHNL/analysis_final.py @@ -3,9 +3,18 @@ #inputDir = "/eos/user/j/jalimena/FCCeeLLP/" #inputDir = "output_stage1/" -#Input directory where the files produced at the final-selection level are +#Output directory where the files produced at the final-selection level are outputDir = "output_finalSel/" +#Integrated luminosity for scaling number of events (required only if setting doScale to true) +#intLumi = 150e6 #pb^-1 + +#Scale event yields by intLumi and cross section (optional) +#doScale = True + +#Save event yields in a table (optional) +#saveTabular = True + processList = { #run over the full statistics from stage1 'p8_ee_Zee_ecm91':{}, @@ -46,6 +55,15 @@ # "sel2RecoEle_vetoes_MissingEnergyGt10_chi2Gt1_LxyzGt5": "n_RecoElectrons==2 && n_RecoMuons==0 && n_RecoPhotons==0 && n_RecoJets==0 && n_RecoPhotons==0 && RecoMissingEnergy_p[0]>10 && RecoDecayVertex.chi2>1 && Reco_Lxyz>5", #displaced vertex } +###Dictionary for prettier names of cuts (optional) +cutLabels = { + "selNone": "Before selection", + "sel2RecoEle": "Exactly 2 electrons", + "sel2RecoEle_vetoes": "Veto photons, muons, and jets", + "sel2RecoEle_vetoes_MissingEnergyGt10": "$\\not\\! p >$ 10 GeV", + "sel2RecoEle_vetoes_MissingEnergyGt10_absD0Gt0p5": "Electron $|d_0| >$ 0.5 mm", +} + ###Dictionary for the ouput variable/hitograms. The key is the name of the variable in the output files. "name" is the name of the variable in the input file, "title" is the x-axis label of the histogram, "bin" the number of bins of the histogram, "xmin" the minimum x-axis value and "xmax" the maximum x-axis value. histoList = { From 864c3228b09991c53c97523b9dfda5cda1043631 Mon Sep 17 00:00:00 2001 From: clementhelsens Date: Mon, 17 Oct 2022 20:02:13 +0200 Subject: [PATCH 02/25] bla --- .../vertexing/Exercises/MyAnalysis.cc | 174 --------------- .../vertexing/Exercises/MyAnalysis.h | 201 ++++++++++++++---- 2 files changed, 156 insertions(+), 219 deletions(-) delete mode 100644 examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.cc diff --git a/examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.cc b/examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.cc deleted file mode 100644 index 4ec759060b..0000000000 --- a/examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.cc +++ /dev/null @@ -1,174 +0,0 @@ -#include "FCCAnalyses/MyAnalysis.h" -#include - -namespace FCCAnalyses{ - -namespace MyAnalysis { - - - double sum_momentum_tracks( const VertexingUtils::FCCAnalysesVertex& vertex) { - double sum = 0; - ROOT::VecOps::RVec< TVector3 > momenta = vertex.updated_track_momentum_at_vertex ; - int n = momenta.size(); - for (int i=0; i < n; i++) { - TVector3 p = momenta[i]; - double px = p[0]; - double py = p[1]; - double pt = sqrt(pow(px,2)+pow(py,2)) ; - sum += pt; - } - return sum; - } - - -// ------------------------------------------------------------------------------------------------- - - double tau3mu_vertex_mass( const VertexingUtils::FCCAnalysesVertex& vertex ) { - double muon_mass = 0.1056; - TLorentzVector tau; - ROOT::VecOps::RVec< TVector3 > momenta = vertex.updated_track_momentum_at_vertex ; - int n = momenta.size(); - for (int ileg=0; ileg < n; ileg++) { - TVector3 track_momentum = momenta[ ileg ]; - TLorentzVector leg; - leg.SetXYZM( track_momentum[0], track_momentum[1], track_momentum[2], muon_mass ) ; - tau += leg; - } - return tau.M(); - } - - - double tau3mu_raw_mass( const ROOT::VecOps::RVec& legs ) { - double muon_mass = 0.1056; - TLorentzVector tau; - int n = legs.size(); - for (int ileg=0; ileg < n; ileg++) { - TLorentzVector leg; - leg.SetXYZM(legs[ileg].momentum.x, legs[ileg].momentum.y, legs[ileg].momentum.z, muon_mass ); - tau += leg; - } - return tau.M(); - } - - -// ------------------------------------------------------------------------------------------------- - -ROOT::VecOps::RVec< ROOT::VecOps::RVec > build_triplets( const ROOT::VecOps::RVec& in , float total_charge ) { - - ROOT::VecOps::RVec< ROOT::VecOps::RVec< edm4hep::ReconstructedParticleData> > results; - float charge =0; - int n = in.size(); - if ( n < 3 ) return results; - - for (int i=0; i < n; i++) { - edm4hep::ReconstructedParticleData pi = in[i]; - float charge_i = pi.charge ; - - for (int j=i+1; j < n; j++) { - edm4hep::ReconstructedParticleData pj = in[j]; - float charge_j = pj.charge ; - - for (int k=j+1; k < n; k++) { - edm4hep::ReconstructedParticleData pk = in[k]; - float charge_k = pk.charge ; - float charge_tot = charge_i + charge_j + charge_k; - if ( charge_tot == total_charge ) { - ROOT::VecOps::RVec a_triplet; - a_triplet.push_back( pi ); - a_triplet.push_back( pj ); - a_triplet.push_back( pk ); - results.push_back( a_triplet ); - } - - } - - } - - } - return results; -} - - -ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex > build_AllTauVertexObject( - const ROOT::VecOps::RVec< ROOT::VecOps::RVec >& triplets, - const ROOT::VecOps::RVec& allTracks ) { - ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex > results; - int ntriplets = triplets.size(); - for (int i=0; i < ntriplets; i++) { - ROOT::VecOps::RVec legs = triplets[i]; - - ROOT::VecOps::RVec the_tracks = ReconstructedParticle2Track::getRP2TRK( legs, allTracks ); - VertexingUtils::FCCAnalysesVertex vertex = VertexFitterSimple::VertexFitter_Tk( 2, the_tracks ); - results.push_back( vertex ); - } - return results; -} - - -ROOT::VecOps::RVec< double > build_AllTauMasses( const ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex>& vertices ) { - ROOT::VecOps::RVec< double > results; - for ( auto& v: vertices) { - double mass = tau3mu_vertex_mass( v ); - results.push_back( mass ); - } - return results; -} - -// ------------------------------------------------------------------------------------------------- - - - -selRP_Fakes::selRP_Fakes( float arg_fakeRate, float arg_mass ) : m_fakeRate(arg_fakeRate), m_mass( arg_mass) { - unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); - std::default_random_engine generator (seed); - m_generator = generator; - std::uniform_real_distribution flatdis(0.,1.); - m_flat.param( flatdis.param() ); -}; - -std::vector selRP_Fakes::operator() (ROOT::VecOps::RVec in) { - std::vector result; - for (size_t i = 0; i < in.size(); ++i) { - auto & p = in[i]; - float arandom = m_flat (m_generator ); - if ( arandom <= m_fakeRate) { - edm4hep::ReconstructedParticleData reso = p; - //reso.momentum.x = p.momentum.x ; - //reso.momentum.y = p.momentum.y ; - //reso.momentum.z = p.momentum.z ; - reso.mass = m_mass; - //reso.charge = p.charge; - result.push_back( reso ); - } - } - return result; -} - - - -// --- for tests... - -float get_p(const edm4hep::MCParticleData& p) { - TLorentzVector tlv; - tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass); - return tlv.P(); -} - -float get_e(const edm4hep::MCParticleData& p) { - TLorentzVector tlv; - tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass); - return tlv.E(); -} - -float get_theta(const edm4hep::MCParticleData& p) { - TLorentzVector tlv; - tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass); - return tlv.Theta(); -} - - - - -} //end NS -} // end NS - diff --git a/examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.h b/examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.h index d90e0bfeef..4ec759060b 100644 --- a/examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.h +++ b/examples/FCCee/tutorials/vertexing/Exercises/MyAnalysis.h @@ -1,63 +1,174 @@ - -#ifndef MYANALYSIS_ANALYZERS_H -#define MYANALYSIS_ANALYZERS_H - -#include -#include - -#include "TLorentzVector.h" -#include "ROOT/RVec.hxx" -#include "edm4hep/ReconstructedParticleData.h" -#include "edm4hep/ParticleIDData.h" -#include "edm4hep/TrackState.h" - -#include "VertexingUtils.h" -#include "ReconstructedParticle2Track.h" -#include "VertexFitterSimple.h" - -#include -#include - +#include "FCCAnalyses/MyAnalysis.h" +#include namespace FCCAnalyses{ namespace MyAnalysis { - double sum_momentum_tracks( const VertexingUtils::FCCAnalysesVertex& vertex ); - - double tau3mu_vertex_mass( const VertexingUtils::FCCAnalysesVertex& vertex ); - double tau3mu_raw_mass( const ROOT::VecOps::RVec& legs ) ; - - ROOT::VecOps::RVec< ROOT::VecOps::RVec > build_triplets( - const ROOT::VecOps::RVec& in , float total_charge) ; + double sum_momentum_tracks( const VertexingUtils::FCCAnalysesVertex& vertex) { + double sum = 0; + ROOT::VecOps::RVec< TVector3 > momenta = vertex.updated_track_momentum_at_vertex ; + int n = momenta.size(); + for (int i=0; i < n; i++) { + TVector3 p = momenta[i]; + double px = p[0]; + double py = p[1]; + double pt = sqrt(pow(px,2)+pow(py,2)) ; + sum += pt; + } + return sum; + } + + +// ------------------------------------------------------------------------------------------------- + + double tau3mu_vertex_mass( const VertexingUtils::FCCAnalysesVertex& vertex ) { + double muon_mass = 0.1056; + TLorentzVector tau; + ROOT::VecOps::RVec< TVector3 > momenta = vertex.updated_track_momentum_at_vertex ; + int n = momenta.size(); + for (int ileg=0; ileg < n; ileg++) { + TVector3 track_momentum = momenta[ ileg ]; + TLorentzVector leg; + leg.SetXYZM( track_momentum[0], track_momentum[1], track_momentum[2], muon_mass ) ; + tau += leg; + } + return tau.M(); + } + + + double tau3mu_raw_mass( const ROOT::VecOps::RVec& legs ) { + double muon_mass = 0.1056; + TLorentzVector tau; + int n = legs.size(); + for (int ileg=0; ileg < n; ileg++) { + TLorentzVector leg; + leg.SetXYZM(legs[ileg].momentum.x, legs[ileg].momentum.y, legs[ileg].momentum.z, muon_mass ); + tau += leg; + } + return tau.M(); + } + + +// ------------------------------------------------------------------------------------------------- + +ROOT::VecOps::RVec< ROOT::VecOps::RVec > build_triplets( const ROOT::VecOps::RVec& in , float total_charge ) { + + ROOT::VecOps::RVec< ROOT::VecOps::RVec< edm4hep::ReconstructedParticleData> > results; + float charge =0; + int n = in.size(); + if ( n < 3 ) return results; + + for (int i=0; i < n; i++) { + edm4hep::ReconstructedParticleData pi = in[i]; + float charge_i = pi.charge ; + + for (int j=i+1; j < n; j++) { + edm4hep::ReconstructedParticleData pj = in[j]; + float charge_j = pj.charge ; + + for (int k=j+1; k < n; k++) { + edm4hep::ReconstructedParticleData pk = in[k]; + float charge_k = pk.charge ; + float charge_tot = charge_i + charge_j + charge_k; + if ( charge_tot == total_charge ) { + ROOT::VecOps::RVec a_triplet; + a_triplet.push_back( pi ); + a_triplet.push_back( pj ); + a_triplet.push_back( pk ); + results.push_back( a_triplet ); + } + + } + + } + + } + return results; +} + + +ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex > build_AllTauVertexObject( + const ROOT::VecOps::RVec< ROOT::VecOps::RVec >& triplets, + const ROOT::VecOps::RVec& allTracks ) { + ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex > results; + int ntriplets = triplets.size(); + for (int i=0; i < ntriplets; i++) { + ROOT::VecOps::RVec legs = triplets[i]; + + ROOT::VecOps::RVec the_tracks = ReconstructedParticle2Track::getRP2TRK( legs, allTracks ); + VertexingUtils::FCCAnalysesVertex vertex = VertexFitterSimple::VertexFitter_Tk( 2, the_tracks ); + results.push_back( vertex ); + } + return results; +} + + +ROOT::VecOps::RVec< double > build_AllTauMasses( const ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex>& vertices ) { + ROOT::VecOps::RVec< double > results; + for ( auto& v: vertices) { + double mass = tau3mu_vertex_mass( v ); + results.push_back( mass ); + } + return results; +} + +// ------------------------------------------------------------------------------------------------- + + + +selRP_Fakes::selRP_Fakes( float arg_fakeRate, float arg_mass ) : m_fakeRate(arg_fakeRate), m_mass( arg_mass) { + unsigned seed = std::chrono::system_clock::now().time_since_epoch().count(); + std::default_random_engine generator (seed); + m_generator = generator; + std::uniform_real_distribution flatdis(0.,1.); + m_flat.param( flatdis.param() ); +}; - ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex > build_AllTauVertexObject( - const ROOT::VecOps::RVec< ROOT::VecOps::RVec >& triplets, - const ROOT::VecOps::RVec& allTracks ) ; +std::vector selRP_Fakes::operator() (ROOT::VecOps::RVec in) { + std::vector result; + for (size_t i = 0; i < in.size(); ++i) { + auto & p = in[i]; + float arandom = m_flat (m_generator ); + if ( arandom <= m_fakeRate) { + edm4hep::ReconstructedParticleData reso = p; + //reso.momentum.x = p.momentum.x ; + //reso.momentum.y = p.momentum.y ; + //reso.momentum.z = p.momentum.z ; + reso.mass = m_mass; + //reso.charge = p.charge; + result.push_back( reso ); + } + } + return result; +} - ROOT::VecOps::RVec< double > build_AllTauMasses( const ROOT::VecOps::RVec< VertexingUtils::FCCAnalysesVertex>& vertices ) ; -struct selRP_Fakes { - selRP_Fakes( float arg_fakeRate, float arg_mass ); - float m_fakeRate = 1e-3; - float m_mass = 0.106; // muon mass - std::default_random_engine m_generator; - std::uniform_real_distribution m_flat; - std::vector operator() (ROOT::VecOps::RVec in); -}; +// --- for tests... +float get_p(const edm4hep::MCParticleData& p) { + TLorentzVector tlv; + tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass); + return tlv.P(); +} +float get_e(const edm4hep::MCParticleData& p) { + TLorentzVector tlv; + tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass); + return tlv.E(); +} +float get_theta(const edm4hep::MCParticleData& p) { + TLorentzVector tlv; + tlv.SetXYZM(p.momentum.x, p.momentum.y, p.momentum.z, p.mass); + return tlv.Theta(); +} -float get_p(const edm4hep::MCParticleData& p) ; -float get_e(const edm4hep::MCParticleData& p) ; -float get_theta(const edm4hep::MCParticleData& in) ; -}//end NS MyAnalysis -}//end NS FCCAnalyses -#endif +} //end NS +} // end NS From ceeedbdf6744f202106c01c95056b2f7122d75c6 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Tue, 18 Oct 2022 15:46:40 +0200 Subject: [PATCH 03/25] add Weaver interface --- analyzers/dataframe/FCCAnalyses/WeaverUtils.h | 27 +++++++++ analyzers/dataframe/src/WeaverUtils.cc | 55 +++++++++++++++++++ 2 files changed, 82 insertions(+) create mode 100644 analyzers/dataframe/FCCAnalyses/WeaverUtils.h create mode 100644 analyzers/dataframe/src/WeaverUtils.cc diff --git a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h new file mode 100644 index 0000000000..0887c303d7 --- /dev/null +++ b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h @@ -0,0 +1,27 @@ +#ifndef FCCAnalyses_WEAVERUTILS_h +#define FCCAnalyses_WEAVERUTILS_h + +#include + +namespace FCCAnalyses { + namespace WeaverUtils { + + + /// Compute all weights given a collection of input variables + /// \note This helper should not be used directly in RDataFrame examples + ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>&); + + /// Setup the ONNXRuntime instance using Weaver-provided parameters + void setup_weaver(const std::string&, const std::string&, const rv::RVec&); + + /// Compute all weights given an unspecified collection of input variables + template + ROOT::VecOps::RVec > get_weights(Args&&... args) { + return compute_weights(std::vector>{std::forward(args)...}); + } + /// Get one specific weight previously computed + rv::RVec get_weight(const ROOT::VecOps::RVec >&, int); + } // namespace JetFlavourUtils +} // namespace FCCAnalyses + +#endif diff --git a/analyzers/dataframe/src/WeaverUtils.cc b/analyzers/dataframe/src/WeaverUtils.cc new file mode 100644 index 0000000000..29215dfb72 --- /dev/null +++ b/analyzers/dataframe/src/WeaverUtils.cc @@ -0,0 +1,55 @@ +#include "FCCAnalyses/WeaverUtils.h" +#include "ONNXRuntime/WeaverInterface.h" + +#include + +namespace FCCAnalyses { + std::unique_ptr gWeaver; + + namespace WeaverUtils { + void setup_weaver(const std::string& onnx_filename, + const std::string& json_filename, + const rv::RVec& vars) { + gWeaver = std::make_unique(onnx_filename, json_filename, vars); + } + + ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>& vars) { + if (!gWeaver) + throw std::runtime_error("Weaver interface is not initialised!"); + ROOT::VecOps::RVec > out; + if (vars.empty()) // no variables registered + return out; + size_t num_obj = vars.at(0).size(); + if (num_obj == 0) // no jets to categorise + return out; + // transform a collection of {var1 -> {jet1 -> {constit1, constit2, ...}, jet2 -> {...}, ...}, var2 -> {...}} + // into a collection of {jet -> {var1 -> {constit1, constit2, ...}, var2 -> {...}, ...}} + + + for (size_t i = 0; i < num_obj; ++i) { + ROOT::VecOps::RVec> obj_sc_vars; + size_t num_constits = vars.at(0).at(i).size(); + for (size_t k = 0; k < vars.size(); ++k) { + ROOT::VecOps::RVec constit_vars; + for (size_t j = 0; j < num_constits; ++j) + constit_vars.push_back((float)vars.at(k).at(i).at(j)); + obj_sc_vars.push_back(constit_vars); + } + out.emplace_back(gWeaver->run(obj_sc_vars)); + } + return out; + } + + ROOT::VecOps::RVec get_weight(const ROOT::VecOps::RVec >& jets_weights, int weight) { + if (weight < 0) + throw std::runtime_error("Invalid index requested for jet flavour weight."); + ROOT::VecOps::RVec out; + for (const auto& jet_weights : jets_weights) { + if (weight >= jet_weights.size()) + throw std::runtime_error("Flavour weight index exceeds the number of weights registered."); + out.emplace_back(jet_weights.at(weight)); + } + return out; + } + } // namespace JetFlavourUtils +} // namespace FCCAnalyses From 6a238a9470bae13362b8cd3afb3f2cc9890637b4 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Tue, 18 Oct 2022 15:51:46 +0200 Subject: [PATCH 04/25] add Weaver interface --- analyzers/dataframe/FCCAnalyses/WeaverUtils.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h index 0887c303d7..6c51157e53 100644 --- a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h +++ b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h @@ -12,7 +12,7 @@ namespace FCCAnalyses { ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>&); /// Setup the ONNXRuntime instance using Weaver-provided parameters - void setup_weaver(const std::string&, const std::string&, const rv::RVec&); + void setup_weaver(const std::string&, const std::string&, const ROOT::VecOps::RVec&); /// Compute all weights given an unspecified collection of input variables template @@ -20,7 +20,7 @@ namespace FCCAnalyses { return compute_weights(std::vector>{std::forward(args)...}); } /// Get one specific weight previously computed - rv::RVec get_weight(const ROOT::VecOps::RVec >&, int); + ROOT::VecOps::RVec get_weight(const ROOT::VecOps::RVec >&, int); } // namespace JetFlavourUtils } // namespace FCCAnalyses From cdc2804e08decdaa87664e52e55065ed2cd44fb7 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Tue, 18 Oct 2022 15:56:17 +0200 Subject: [PATCH 05/25] add Weaver interface --- analyzers/dataframe/FCCAnalyses/WeaverUtils.h | 2 +- analyzers/dataframe/src/WeaverUtils.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h index 6c51157e53..4d2266c204 100644 --- a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h +++ b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h @@ -9,7 +9,7 @@ namespace FCCAnalyses { /// Compute all weights given a collection of input variables /// \note This helper should not be used directly in RDataFrame examples - ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>&); + ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>>&); /// Setup the ONNXRuntime instance using Weaver-provided parameters void setup_weaver(const std::string&, const std::string&, const ROOT::VecOps::RVec&); diff --git a/analyzers/dataframe/src/WeaverUtils.cc b/analyzers/dataframe/src/WeaverUtils.cc index 29215dfb72..5c9293e56e 100644 --- a/analyzers/dataframe/src/WeaverUtils.cc +++ b/analyzers/dataframe/src/WeaverUtils.cc @@ -13,7 +13,7 @@ namespace FCCAnalyses { gWeaver = std::make_unique(onnx_filename, json_filename, vars); } - ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>& vars) { + ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>>& vars) { if (!gWeaver) throw std::runtime_error("Weaver interface is not initialised!"); ROOT::VecOps::RVec > out; From bcac230d0f5e6b20945991c73cdbf038bef5abde Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Tue, 18 Oct 2022 15:57:50 +0200 Subject: [PATCH 06/25] add Weaver interface --- analyzers/dataframe/FCCAnalyses/WeaverUtils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h index 4d2266c204..942f12bcc8 100644 --- a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h +++ b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h @@ -17,7 +17,7 @@ namespace FCCAnalyses { /// Compute all weights given an unspecified collection of input variables template ROOT::VecOps::RVec > get_weights(Args&&... args) { - return compute_weights(std::vector>{std::forward(args)...}); + return compute_weights(std::vector>>{std::forward(args)...}); } /// Get one specific weight previously computed ROOT::VecOps::RVec get_weight(const ROOT::VecOps::RVec >&, int); From 816cff52e5ecf1a28b79d0ca75ea198f54488443 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Tue, 18 Oct 2022 15:59:51 +0200 Subject: [PATCH 07/25] add Weaver interface --- analyzers/dataframe/src/WeaverUtils.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/analyzers/dataframe/src/WeaverUtils.cc b/analyzers/dataframe/src/WeaverUtils.cc index 5c9293e56e..f075b4bebd 100644 --- a/analyzers/dataframe/src/WeaverUtils.cc +++ b/analyzers/dataframe/src/WeaverUtils.cc @@ -4,17 +4,17 @@ #include namespace FCCAnalyses { - std::unique_ptr gWeaver; + std::unique_ptr gWeaver2; namespace WeaverUtils { void setup_weaver(const std::string& onnx_filename, const std::string& json_filename, const rv::RVec& vars) { - gWeaver = std::make_unique(onnx_filename, json_filename, vars); + gWeaver2 = std::make_unique(onnx_filename, json_filename, vars); } ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>>& vars) { - if (!gWeaver) + if (!gWeaver2) throw std::runtime_error("Weaver interface is not initialised!"); ROOT::VecOps::RVec > out; if (vars.empty()) // no variables registered @@ -35,7 +35,7 @@ namespace FCCAnalyses { constit_vars.push_back((float)vars.at(k).at(i).at(j)); obj_sc_vars.push_back(constit_vars); } - out.emplace_back(gWeaver->run(obj_sc_vars)); + out.emplace_back(gWeaver2->run(obj_sc_vars)); } return out; } From 4f1aa354dfadb1aaef4726ab6af296d075aff5cc Mon Sep 17 00:00:00 2001 From: Laurent Forthomme Date: Wed, 19 Oct 2022 03:04:35 +0200 Subject: [PATCH 08/25] WeaverUtils objects remnaming --- analyzers/dataframe/FCCAnalyses/WeaverUtils.h | 10 +-- analyzers/dataframe/src/WeaverUtils.cc | 35 +++++----- .../tutorials/analysis_CaloClustering.py | 69 +++++++++++++++++++ 3 files changed, 92 insertions(+), 22 deletions(-) create mode 100644 examples/FCCee/tutorials/analysis_CaloClustering.py diff --git a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h index 942f12bcc8..39bd043594 100644 --- a/analyzers/dataframe/FCCAnalyses/WeaverUtils.h +++ b/analyzers/dataframe/FCCAnalyses/WeaverUtils.h @@ -6,22 +6,22 @@ namespace FCCAnalyses { namespace WeaverUtils { - /// Compute all weights given a collection of input variables /// \note This helper should not be used directly in RDataFrame examples - ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>>&); + ROOT::VecOps::RVec> compute_weights( + const ROOT::VecOps::RVec>>&); /// Setup the ONNXRuntime instance using Weaver-provided parameters void setup_weaver(const std::string&, const std::string&, const ROOT::VecOps::RVec&); /// Compute all weights given an unspecified collection of input variables template - ROOT::VecOps::RVec > get_weights(Args&&... args) { + ROOT::VecOps::RVec> get_weights(Args&&... args) { return compute_weights(std::vector>>{std::forward(args)...}); } /// Get one specific weight previously computed - ROOT::VecOps::RVec get_weight(const ROOT::VecOps::RVec >&, int); - } // namespace JetFlavourUtils + ROOT::VecOps::RVec get_weight(const ROOT::VecOps::RVec>&, int); + } // namespace WeaverUtils } // namespace FCCAnalyses #endif diff --git a/analyzers/dataframe/src/WeaverUtils.cc b/analyzers/dataframe/src/WeaverUtils.cc index f075b4bebd..b7b6d0b270 100644 --- a/analyzers/dataframe/src/WeaverUtils.cc +++ b/analyzers/dataframe/src/WeaverUtils.cc @@ -13,43 +13,44 @@ namespace FCCAnalyses { gWeaver2 = std::make_unique(onnx_filename, json_filename, vars); } - ROOT::VecOps::RVec > compute_weights(const ROOT::VecOps::RVec>>& vars) { + ROOT::VecOps::RVec> compute_weights( + const ROOT::VecOps::RVec>>& vars) { if (!gWeaver2) throw std::runtime_error("Weaver interface is not initialised!"); - ROOT::VecOps::RVec > out; + ROOT::VecOps::RVec> out; if (vars.empty()) // no variables registered return out; size_t num_obj = vars.at(0).size(); - if (num_obj == 0) // no jets to categorise + if (num_obj == 0) // no objects to categorise return out; - // transform a collection of {var1 -> {jet1 -> {constit1, constit2, ...}, jet2 -> {...}, ...}, var2 -> {...}} - // into a collection of {jet -> {var1 -> {constit1, constit2, ...}, var2 -> {...}, ...}} - + // transform a collection of {var1 -> {object1 -> {input1, input2, ...}, object2 -> {...}, ...}, var2 -> {...}} + // into a collection of {object -> {var1 -> {input1, input2, ...}, var2 -> {...}, ...}} for (size_t i = 0; i < num_obj; ++i) { ROOT::VecOps::RVec> obj_sc_vars; - size_t num_constits = vars.at(0).at(i).size(); + size_t num_inputs = vars.at(0).at(i).size(); for (size_t k = 0; k < vars.size(); ++k) { - ROOT::VecOps::RVec constit_vars; - for (size_t j = 0; j < num_constits; ++j) - constit_vars.push_back((float)vars.at(k).at(i).at(j)); - obj_sc_vars.push_back(constit_vars); + ROOT::VecOps::RVec input_vars; + for (size_t j = 0; j < num_inputs; ++j) + input_vars.push_back((float)vars.at(k).at(i).at(j)); + obj_sc_vars.push_back(input_vars); } out.emplace_back(gWeaver2->run(obj_sc_vars)); } return out; } - ROOT::VecOps::RVec get_weight(const ROOT::VecOps::RVec >& jets_weights, int weight) { + ROOT::VecOps::RVec get_weight(const ROOT::VecOps::RVec>& objects_weights, + int weight) { if (weight < 0) - throw std::runtime_error("Invalid index requested for jet flavour weight."); + throw std::runtime_error("Invalid index requested for object weight " + std::to_string(weight) + "."); ROOT::VecOps::RVec out; - for (const auto& jet_weights : jets_weights) { - if (weight >= jet_weights.size()) + for (const auto& object_weights : objects_weights) { + if (weight >= object_weights.size()) throw std::runtime_error("Flavour weight index exceeds the number of weights registered."); - out.emplace_back(jet_weights.at(weight)); + out.emplace_back(object_weights.at(weight)); } return out; } - } // namespace JetFlavourUtils + } // namespace WeaverUtils } // namespace FCCAnalyses diff --git a/examples/FCCee/tutorials/analysis_CaloClustering.py b/examples/FCCee/tutorials/analysis_CaloClustering.py new file mode 100644 index 0000000000..716ee73da2 --- /dev/null +++ b/examples/FCCee/tutorials/analysis_CaloClustering.py @@ -0,0 +1,69 @@ + +#Mandatory: List of processes +processList = { + #'p8_ee_ZH_ecm240':{'fraction':0.2, 'chunks':2, 'output':'p8_ee_ZH_ecm240_out'} + 'p8_noBES_ee_H_Hbb_ecm125':{'fraction':0.01, 'chunks':1, 'output':'test_out'} +} + +#Mandatory: Production tag when running over EDM4Hep centrally produced events, this points to the yaml files for getting sample statistics +prodTag = "FCCee/spring2021/IDEA/" + +#Optional: output directory, default is local running directory +outputDir = "." + +#Optional +nCPUS = 8 +runBatch = False +#batchQueue = "longlunch" +#compGroup = "group_u_FCC.local_gen" + +#Optional test file +testFile = "/eos/experiment/fcc/ee/tmp/4Lolo/photons_MVA1.root" + +#Mandatory: RDFanalysis class where the use defines the operations on the TTree +class RDFanalysis(): + #__________________________________________________________ + #Mandatory: analysers funtion to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 + def analysers(df): + from ROOT import WeaverUtils + from ROOT import gInterpreter + from os import getenv + test_inputs_path = getenv('TEST_INPUT_DATA_DIR', '/eos/experiment/fcc/ee/tutorial/PNet_pi0Gamma/v1') + + gInterpreter.Declare(""" + template + ROOT::VecOps::RVec > get_vector(const ROOT::VecOps::RVec& qty) { return ROOT::VecOps::RVec >(1, qty); } + """) + weaver = WeaverUtils.setup_weaver(test_inputs_path + '/fccee_pi_vs_gamma_v1.onnx', + test_inputs_path + '/preprocess_fccee_pi_vs_gamma_v1.json', + ('recocells_e', 'recocells_theta', 'recocells_phi', 'recocells_radius', 'recocells_layer')) + + df2 = (df + #.Define("cells_e", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_energy); return v;") + #.Define("cells_theta", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_theta); return v;") + #.Define("cells_phi", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_phi); return v;") + #.Define("cells_radius", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_radius); return v;") + #.Define("cells_layer", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_layer); return v;") + .Define("cells_e", "get_vector(maxEnergyCluster_cells_energy)") + .Define("cells_theta", "get_vector(maxEnergyCluster_cells_theta)") + .Define("cells_phi", "get_vector(maxEnergyCluster_cells_phi)") + .Define("cells_radius", "get_vector(maxEnergyCluster_cells_radius)") + .Define("cells_layer", "get_vector(maxEnergyCluster_cells_layer)") + + # run inference + .Define("MVAVec", "WeaverUtils::get_weights(cells_e, cells_theta, cells_phi, cells_radius, cells_layer)") + #.Define("MVAVec", "WeaverUtils::get_weights(maxEnergyCluster_cells_energy, maxEnergyCluster_cells_theta, maxEnergyCluster_cells_phi, maxEnergyCluster_cells_radius, maxEnergyCluster_cells_layer)") + + # recast output + .Define("Cluster_isPhoton", "WeaverUtils::get_weight(MVAVec, 0)") + .Define("Cluster_isPi0", "WeaverUtils::get_weight(MVAVec, 1)") + ) + return df2 + + #__________________________________________________________ + #Mandatory: output function, please make sure you return the branchlist as a python list + def output(): + branchList = [ + 'Cluster_isPhoton', 'Cluster_isPi0', + ] + return branchList From 7e98868f363d41d15026fbe8ce05b44f81ddcb37 Mon Sep 17 00:00:00 2001 From: Laurent Forthomme Date: Wed, 19 Oct 2022 03:13:04 +0200 Subject: [PATCH 09/25] Cleanup --- .../tutorials/analysis_CaloClustering.py | 22 ++++++++----------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/examples/FCCee/tutorials/analysis_CaloClustering.py b/examples/FCCee/tutorials/analysis_CaloClustering.py index 716ee73da2..c4ba4e7713 100644 --- a/examples/FCCee/tutorials/analysis_CaloClustering.py +++ b/examples/FCCee/tutorials/analysis_CaloClustering.py @@ -32,27 +32,23 @@ def analysers(df): gInterpreter.Declare(""" template - ROOT::VecOps::RVec > get_vector(const ROOT::VecOps::RVec& qty) { return ROOT::VecOps::RVec >(1, qty); } - """) + ROOT::VecOps::RVec > as_vector(const ROOT::VecOps::RVec& qty) { + return ROOT::VecOps::RVec >(1, qty); + }""") weaver = WeaverUtils.setup_weaver(test_inputs_path + '/fccee_pi_vs_gamma_v1.onnx', test_inputs_path + '/preprocess_fccee_pi_vs_gamma_v1.json', ('recocells_e', 'recocells_theta', 'recocells_phi', 'recocells_radius', 'recocells_layer')) df2 = (df - #.Define("cells_e", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_energy); return v;") - #.Define("cells_theta", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_theta); return v;") - #.Define("cells_phi", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_phi); return v;") - #.Define("cells_radius", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_radius); return v;") - #.Define("cells_layer", "ROOT::VecOps::RVec> v; v.push_back(maxEnergyCluster_cells_layer); return v;") - .Define("cells_e", "get_vector(maxEnergyCluster_cells_energy)") - .Define("cells_theta", "get_vector(maxEnergyCluster_cells_theta)") - .Define("cells_phi", "get_vector(maxEnergyCluster_cells_phi)") - .Define("cells_radius", "get_vector(maxEnergyCluster_cells_radius)") - .Define("cells_layer", "get_vector(maxEnergyCluster_cells_layer)") + # define input variables + .Define("cells_e", "as_vector(maxEnergyCluster_cells_energy)") + .Define("cells_theta", "as_vector(maxEnergyCluster_cells_theta)") + .Define("cells_phi", "as_vector(maxEnergyCluster_cells_phi)") + .Define("cells_radius", "as_vector(maxEnergyCluster_cells_radius)") + .Define("cells_layer", "as_vector(maxEnergyCluster_cells_layer)") # run inference .Define("MVAVec", "WeaverUtils::get_weights(cells_e, cells_theta, cells_phi, cells_radius, cells_layer)") - #.Define("MVAVec", "WeaverUtils::get_weights(maxEnergyCluster_cells_energy, maxEnergyCluster_cells_theta, maxEnergyCluster_cells_phi, maxEnergyCluster_cells_radius, maxEnergyCluster_cells_layer)") # recast output .Define("Cluster_isPhoton", "WeaverUtils::get_weight(MVAVec, 0)") From 9e585df716a0ca44bfe764fa6d8601ade501c4b6 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Wed, 19 Oct 2022 08:06:17 +0200 Subject: [PATCH 10/25] add analysis generator --- .../FCCee/tutorials/analysis_generator.py | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 examples/FCCee/tutorials/analysis_generator.py diff --git a/examples/FCCee/tutorials/analysis_generator.py b/examples/FCCee/tutorials/analysis_generator.py new file mode 100644 index 0000000000..14a4abdd4c --- /dev/null +++ b/examples/FCCee/tutorials/analysis_generator.py @@ -0,0 +1,26 @@ +testFile='http://fccsw.web.cern.ch/tutorials/october2020/tutorial1/kk_tautau_10000.e4h.root' + +# Mandatory: RDFanalysis class where the user defines the operations on the TTree +class RDFanalysis(): + #__________________________________________________________ + # Mandatory: analysers function to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 + def analysers(df): + df2 = (df + #.Alias("GenParticles0", "GenParticles#0.index") + .Define("MC_p", "FCCAnalyses::MCParticle::get_p(GenParticles)") + .Define("MC_theta","FCCAnalyses::MCParticle::get_theta(GenParticles)") + .Define("MC_charge","FCCAnalyses::MCParticle::get_charge(GenParticles)") + .Define("MC_phi","FCCAnalyses::MCParticle::get_phi(GenParticles)") + ) + return df2 + + #__________________________________________________________ + # Mandatory: output function, please make sure you return the branchlist as a python list + def output(): + branchList = [ + "MC_p", + "MC_theta", + "MC_charge", + "MC_phi", + ] + return branchList From ff0af40ff68bba2734eeaa789c7b94485a5b998e Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Wed, 19 Oct 2022 08:08:33 +0200 Subject: [PATCH 11/25] clean up --- examples/FCCee/tutorials/analysis_generator.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/FCCee/tutorials/analysis_generator.py b/examples/FCCee/tutorials/analysis_generator.py index 14a4abdd4c..7c6e06d119 100644 --- a/examples/FCCee/tutorials/analysis_generator.py +++ b/examples/FCCee/tutorials/analysis_generator.py @@ -6,12 +6,12 @@ class RDFanalysis(): # Mandatory: analysers function to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 def analysers(df): df2 = (df - #.Alias("GenParticles0", "GenParticles#0.index") - .Define("MC_p", "FCCAnalyses::MCParticle::get_p(GenParticles)") - .Define("MC_theta","FCCAnalyses::MCParticle::get_theta(GenParticles)") - .Define("MC_charge","FCCAnalyses::MCParticle::get_charge(GenParticles)") - .Define("MC_phi","FCCAnalyses::MCParticle::get_phi(GenParticles)") - ) + #namespace FCCAnalyses loaded, could be removed in case + .Define("MC_p", "FCCAnalyses::MCParticle::get_p(GenParticles)") + .Define("MC_theta", "FCCAnalyses::MCParticle::get_theta(GenParticles)") + .Define("MC_charge", "FCCAnalyses::MCParticle::get_charge(GenParticles)") + .Define("MC_phi", "FCCAnalyses::MCParticle::get_phi(GenParticles)") + ) return df2 #__________________________________________________________ From 94c7fcf1dff8cae4436e8558132208a3e8f8d204 Mon Sep 17 00:00:00 2001 From: clementhelsens Date: Wed, 19 Oct 2022 11:15:10 +0200 Subject: [PATCH 12/25] more on gen --- .../FCCee/tutorials/analysis_generator.py | 28 +++++++++++++++---- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/examples/FCCee/tutorials/analysis_generator.py b/examples/FCCee/tutorials/analysis_generator.py index 7c6e06d119..c990ad5ae3 100644 --- a/examples/FCCee/tutorials/analysis_generator.py +++ b/examples/FCCee/tutorials/analysis_generator.py @@ -5,12 +5,28 @@ class RDFanalysis(): #__________________________________________________________ # Mandatory: analysers function to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 def analysers(df): + import ROOT + ROOT.gInterpreter.Declare(""" + + float scalarProduct(const edm4hep::Vector3f& in1, const edm4hep::Vector3f& in2 ){ + + return ((in1.x*in2.x + in1.y*in2.y + in1.z*in2.z)/sqrt(pow(in1.x,2)+pow(in1.y,2)+pow(in1.z,2))/sqrt(pow(in2.x,2)+pow(in2.y,2)+pow(in2.z,2) )); +}; +""") + df2 = (df - #namespace FCCAnalyses loaded, could be removed in case - .Define("MC_p", "FCCAnalyses::MCParticle::get_p(GenParticles)") - .Define("MC_theta", "FCCAnalyses::MCParticle::get_theta(GenParticles)") - .Define("MC_charge", "FCCAnalyses::MCParticle::get_charge(GenParticles)") - .Define("MC_phi", "FCCAnalyses::MCParticle::get_phi(GenParticles)") + + #namespace FCCAnalyses loaded, could be removed in case + .Define("MC_p", "FCCAnalyses::MCParticle::get_p(GenParticles)") + .Define("MC_theta", "FCCAnalyses::MCParticle::get_theta(GenParticles)") + .Define("MC_charge", "FCCAnalyses::MCParticle::get_charge(GenParticles)") + .Define("MC_phi", "FCCAnalyses::MCParticle::get_phi(GenParticles)") + .Define("tauplusvec", ROOT.MCParticle.sel_pdgID(-15, 0),["GenParticles"]) + .Define("tauminusvec", ROOT.MCParticle.sel_pdgID(15, 0),["GenParticles"]) + .Filter("tauminusvec.size()==1") + .Define("tauplus", "tauplusvec[0]") + .Define("scalarProd", "scalarProduct(tauminusvec[0].momentum, tauplusvec[0].momentum)") + .Define("tauplus_px", "tauplus.momentum.x") ) return df2 @@ -22,5 +38,7 @@ def output(): "MC_theta", "MC_charge", "MC_phi", + "tauplus_px", + "scalarProd" ] return branchList From 954797bd0e609c7427173492e667546b9e92f131 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Wed, 19 Oct 2022 11:37:20 +0200 Subject: [PATCH 13/25] change to cvmfs --- config/common_defaults.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/config/common_defaults.py b/config/common_defaults.py index b6c1668105..f9453ac34e 100644 --- a/config/common_defaults.py +++ b/config/common_defaults.py @@ -1 +1,2 @@ -deffccdicts = "/afs/cern.ch/work/f/fccsw/public/FCCDicts/" +#deffccdicts = "/afs/cern.ch/work/f/fccsw/public/FCCDicts/" +deffccdicts = "/cvmfs/fcc.cern.ch/FCCDicts/" From a0ca3c875e2bb8754ec6331e0d69402482fbd952 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Wed, 19 Oct 2022 14:47:36 +0200 Subject: [PATCH 14/25] template vec --- analyzers/dataframe/FCCAnalyses/Utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analyzers/dataframe/FCCAnalyses/Utils.h b/analyzers/dataframe/FCCAnalyses/Utils.h index b4f7301259..78c5616d0d 100644 --- a/analyzers/dataframe/FCCAnalyses/Utils.h +++ b/analyzers/dataframe/FCCAnalyses/Utils.h @@ -7,7 +7,7 @@ namespace FCCAnalyses { namespace Utils { template inline auto getsize( T& vec){ return vec.size();}; - + template inline ROOT::VecOps::RVec > as_vector(const ROOT::VecOps::RVec& in){return ROOT::VecOps::RVec >(1, in);}; } } From 2fdc73015c078f9d73573155a7e16d0221129f42 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Wed, 19 Oct 2022 15:46:58 +0200 Subject: [PATCH 15/25] fix example --- .../tutorials/analysis_CaloClustering.py | 34 +++---------------- 1 file changed, 5 insertions(+), 29 deletions(-) diff --git a/examples/FCCee/tutorials/analysis_CaloClustering.py b/examples/FCCee/tutorials/analysis_CaloClustering.py index c4ba4e7713..f4faa5d7f9 100644 --- a/examples/FCCee/tutorials/analysis_CaloClustering.py +++ b/examples/FCCee/tutorials/analysis_CaloClustering.py @@ -1,22 +1,3 @@ - -#Mandatory: List of processes -processList = { - #'p8_ee_ZH_ecm240':{'fraction':0.2, 'chunks':2, 'output':'p8_ee_ZH_ecm240_out'} - 'p8_noBES_ee_H_Hbb_ecm125':{'fraction':0.01, 'chunks':1, 'output':'test_out'} -} - -#Mandatory: Production tag when running over EDM4Hep centrally produced events, this points to the yaml files for getting sample statistics -prodTag = "FCCee/spring2021/IDEA/" - -#Optional: output directory, default is local running directory -outputDir = "." - -#Optional -nCPUS = 8 -runBatch = False -#batchQueue = "longlunch" -#compGroup = "group_u_FCC.local_gen" - #Optional test file testFile = "/eos/experiment/fcc/ee/tmp/4Lolo/photons_MVA1.root" @@ -30,22 +11,17 @@ def analysers(df): from os import getenv test_inputs_path = getenv('TEST_INPUT_DATA_DIR', '/eos/experiment/fcc/ee/tutorial/PNet_pi0Gamma/v1') - gInterpreter.Declare(""" - template - ROOT::VecOps::RVec > as_vector(const ROOT::VecOps::RVec& qty) { - return ROOT::VecOps::RVec >(1, qty); - }""") weaver = WeaverUtils.setup_weaver(test_inputs_path + '/fccee_pi_vs_gamma_v1.onnx', test_inputs_path + '/preprocess_fccee_pi_vs_gamma_v1.json', ('recocells_e', 'recocells_theta', 'recocells_phi', 'recocells_radius', 'recocells_layer')) df2 = (df # define input variables - .Define("cells_e", "as_vector(maxEnergyCluster_cells_energy)") - .Define("cells_theta", "as_vector(maxEnergyCluster_cells_theta)") - .Define("cells_phi", "as_vector(maxEnergyCluster_cells_phi)") - .Define("cells_radius", "as_vector(maxEnergyCluster_cells_radius)") - .Define("cells_layer", "as_vector(maxEnergyCluster_cells_layer)") + .Define("cells_e", "Utils::as_vector(maxEnergyCluster_cells_energy)") + .Define("cells_theta", "Utils::as_vector(maxEnergyCluster_cells_theta)") + .Define("cells_phi", "Utils::as_vector(maxEnergyCluster_cells_phi)") + .Define("cells_radius", "Utils::as_vector(maxEnergyCluster_cells_radius)") + .Define("cells_layer", "Utils::as_vector(maxEnergyCluster_cells_layer)") # run inference .Define("MVAVec", "WeaverUtils::get_weights(cells_e, cells_theta, cells_phi, cells_radius, cells_layer)") From cc0e4fabe408538b26882d106f64b5462fc50ac5 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Wed, 19 Oct 2022 17:22:39 +0200 Subject: [PATCH 16/25] test meta --- config/FCCAnalysisRun.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/config/FCCAnalysisRun.py b/config/FCCAnalysisRun.py index 48a2058908..d670363182 100644 --- a/config/FCCAnalysisRun.py +++ b/config/FCCAnalysisRun.py @@ -516,6 +516,19 @@ def runLocal(rdfModule, fileList, args): if nevents_meta>nevents_local:n[0]=nevents_meta p = ROOT.TParameter(int)( "eventsProcessed", n[0]) p.Write() + + if args.test: + outf2 = ROOT.TFile(fileListRoot[0]) + outt2_1 = outf2.Get("metadata") + outt2_2 = outf2.Get("run_metadata") + outt2_3 = outf2.Get("evt_metadata") + outt2_4 = outf2.Get("col_metadata") + outf.cd() + outt2_1.Write() + outt2_2.Write() + outt2_3.Write() + outt2_4.Write() + outf.Write() outf.Close() From 57a4e1c04a5197746678df6e4064abb26a514567 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 10:15:59 +0200 Subject: [PATCH 17/25] test bug fix --- addons/ONNXRuntime/src/WeaverInterface.cc | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/addons/ONNXRuntime/src/WeaverInterface.cc b/addons/ONNXRuntime/src/WeaverInterface.cc index d698f40eec..e278730bd2 100644 --- a/addons/ONNXRuntime/src/WeaverInterface.cc +++ b/addons/ONNXRuntime/src/WeaverInterface.cc @@ -102,8 +102,13 @@ rv::RVec WeaverInterface::run(const rv::RVec& constituen ConstituentVars jc; for (size_t j = 0; j < params.var_names.size(); ++j) { // transform and add the proper amount of padding const auto& var_name = params.var_names.at(j); - if (std::find(variables_names_.begin(), variables_names_.end(), "pfcand_mask") == variables_names_.end()) - jc = ConstituentVars(constituents.at(0).size(), 1.f); + //if (std::find(variables_names_.begin(), variables_names_.end(), "pfcand_mask") == variables_names_.end()) + // jc = ConstituentVars(constituents.at(0).size(), 1.f); + + if (var_name.find("_mask") != std::string::npos) { + jc = ConstituentVars(constituents.at(0).size(), 1.f); + } + else jc = constituents.at(variablePos(var_name)); const auto& var_info = params.info(var_name); From 5c9f7dfb1ebfe00b639bf55374aafd3d689dea9e Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 11:40:40 +0200 Subject: [PATCH 18/25] comment meta --- config/FCCAnalysisRun.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/config/FCCAnalysisRun.py b/config/FCCAnalysisRun.py index d670363182..e72ac53493 100644 --- a/config/FCCAnalysisRun.py +++ b/config/FCCAnalysisRun.py @@ -517,17 +517,17 @@ def runLocal(rdfModule, fileList, args): p = ROOT.TParameter(int)( "eventsProcessed", n[0]) p.Write() - if args.test: - outf2 = ROOT.TFile(fileListRoot[0]) - outt2_1 = outf2.Get("metadata") - outt2_2 = outf2.Get("run_metadata") - outt2_3 = outf2.Get("evt_metadata") - outt2_4 = outf2.Get("col_metadata") - outf.cd() - outt2_1.Write() - outt2_2.Write() - outt2_3.Write() - outt2_4.Write() +# if args.test: +# outf2 = ROOT.TFile(fileListRoot[0]) +# outt2_1 = outf2.Get("metadata") +# outt2_2 = outf2.Get("run_metadata") +# outt2_3 = outf2.Get("evt_metadata") +# outt2_4 = outf2.Get("col_metadata") +# outf.cd() +# outt2_1.Write() +# outt2_2.Write() +# outt2_3.Write() +# outt2_4.Write() outf.Write() outf.Close() From 30c654d3dba03675e35adaaa35eff44f086e3e17 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 13:08:53 +0200 Subject: [PATCH 19/25] add filter --- .../dataframe/FCCAnalyses/CaloNtupleizer.h | 12 +++ analyzers/dataframe/src/CaloNtupleizer.cc | 14 +++- .../tutorials/analysis_CaloClustering.py | 75 +++++++++++++++---- 3 files changed, 85 insertions(+), 16 deletions(-) diff --git a/analyzers/dataframe/FCCAnalyses/CaloNtupleizer.h b/analyzers/dataframe/FCCAnalyses/CaloNtupleizer.h index 8306a927ed..1b9b9a8291 100644 --- a/analyzers/dataframe/FCCAnalyses/CaloNtupleizer.h +++ b/analyzers/dataframe/FCCAnalyses/CaloNtupleizer.h @@ -19,6 +19,18 @@ namespace CaloNtupleizer{ void loadGeometry(std::string xmlGeometryPath, std::string readoutName); +/// select layers +struct sel_layers { +public: + sel_layers(int arg_min=0, int arg_max=10); + ROOT::VecOps::RVec operator()(const ROOT::VecOps::RVec& in); + +private: + int _min;//> min layer + int _max;//> max layer +}; + + // calo hits (single cells) ROOT::VecOps::RVec getCaloHit_x (const ROOT::VecOps::RVec& in); ROOT::VecOps::RVec getCaloHit_y (const ROOT::VecOps::RVec& in); diff --git a/analyzers/dataframe/src/CaloNtupleizer.cc b/analyzers/dataframe/src/CaloNtupleizer.cc index 65050269af..daee5fdade 100644 --- a/analyzers/dataframe/src/CaloNtupleizer.cc +++ b/analyzers/dataframe/src/CaloNtupleizer.cc @@ -11,7 +11,7 @@ namespace FCCAnalyses{ namespace CaloNtupleizer{ -dd4hep::DDSegmentation::BitFieldCoder* m_decoder; +std::unique_ptr m_decoder; void loadGeometry(std::string xmlGeometryPath, std::string readoutName){ dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); @@ -22,6 +22,18 @@ void loadGeometry(std::string xmlGeometryPath, std::string readoutName){ } +sel_layers::sel_layers(int arg_min, int arg_max) : _min(arg_min), _max(arg_max) {}; +ROOT::VecOps::RVec sel_layers::operator() (const ROOT::VecOps::RVec& in) { + + ROOT::VecOps::RVec res; + for (auto & p: in){ + dd4hep::DDSegmentation::CellID cellId = p.cellID; + int layer = m_decoder->get(cellId, "layer"); + if (layer>_min || layer<_max)res.emplace_back(p); + } + return res; +} + // calo hit ROOT::VecOps::RVec getCaloHit_x (const ROOT::VecOps::RVec& in){ ROOT::VecOps::RVec result; diff --git a/examples/FCCee/tutorials/analysis_CaloClustering.py b/examples/FCCee/tutorials/analysis_CaloClustering.py index f4faa5d7f9..b02fa15087 100644 --- a/examples/FCCee/tutorials/analysis_CaloClustering.py +++ b/examples/FCCee/tutorials/analysis_CaloClustering.py @@ -1,34 +1,76 @@ -#Optional test file -testFile = "/eos/experiment/fcc/ee/tmp/4Lolo/photons_MVA1.root" +from os import getenv +FCCDETECTORS=getenv('FCCDETECTORS') +geometryFile = FCCDETECTORS+"/Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.xml" +readoutName = "ECalBarrelPhiEta" + +import ROOT +ROOT.gInterpreter.Declare(""" +template +ROOT::VecOps::RVec myRange(ROOT::VecOps::RVec& vec, std::size_t begin, std::size_t end) +{ + ROOT::VecOps::RVec ret; + ret.reserve(end - begin); + for (auto i = begin; i < end; ++i) + ret.push_back(vec[i]); + return ret; +} +""") + + #Mandatory: RDFanalysis class where the use defines the operations on the TTree class RDFanalysis(): #__________________________________________________________ #Mandatory: analysers funtion to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 def analysers(df): + from ROOT import WeaverUtils - from ROOT import gInterpreter from os import getenv test_inputs_path = getenv('TEST_INPUT_DATA_DIR', '/eos/experiment/fcc/ee/tutorial/PNet_pi0Gamma/v1') - weaver = WeaverUtils.setup_weaver(test_inputs_path + '/fccee_pi_vs_gamma_v1.onnx', test_inputs_path + '/preprocess_fccee_pi_vs_gamma_v1.json', ('recocells_e', 'recocells_theta', 'recocells_phi', 'recocells_radius', 'recocells_layer')) + df2 = (df - # define input variables - .Define("cells_e", "Utils::as_vector(maxEnergyCluster_cells_energy)") - .Define("cells_theta", "Utils::as_vector(maxEnergyCluster_cells_theta)") - .Define("cells_phi", "Utils::as_vector(maxEnergyCluster_cells_phi)") - .Define("cells_radius", "Utils::as_vector(maxEnergyCluster_cells_radius)") - .Define("cells_layer", "Utils::as_vector(maxEnergyCluster_cells_layer)") - - # run inference + + .Define("maxEnergyCluster_index","std::distance(CaloClusters.energy.begin(),std::max_element(CaloClusters.energy.begin(), CaloClusters.energy.end()))") + .Define("minEnergyCluster_index","std::distance(CaloClusters.energy.begin(),std::min_element(CaloClusters.energy.begin(), CaloClusters.energy.end()))") + .Define("clusters_n","CaloClusters.energy.size()") + .Define("clusters_energy","CaloClusters.energy") + + .Define("maxEnergyCluster_firstCell_index","CaloClusters[maxEnergyCluster_index].hits_begin") + .Define("maxEnergyCluster_lastCell_index" ,"CaloClusters[maxEnergyCluster_index].hits_end") + .Define("maxEnergyCluster_cells", "myRange(PositionedCaloClusterCells, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") + + .Define("maxEnergyCluster_cells_energy", "CaloNtupleizer::getCaloHit_energy(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_phi", "CaloNtupleizer::getCaloHit_phi(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_theta", "CaloNtupleizer::getCaloHit_theta(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_layer" , "CaloNtupleizer::getCaloHit_layer(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_n" , "maxEnergyCluster_cells.size()" ) + + .Define("maxEnergyCluster_cells_x", "myRange(PositionedCaloClusterCells.position.x, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") + .Define("maxEnergyCluster_cells_y", "myRange(PositionedCaloClusterCells.position.y, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") + .Define("maxEnergyCluster_cells_radius", "sqrt(pow(maxEnergyCluster_cells_x,2)+pow(maxEnergyCluster_cells_y,2))") + + .Define("maxEnergyCluster_energy", "CaloClusters.energy[maxEnergyCluster_index]") + .Define("maxEnergyCluster_x", "CaloClusters.position.x[maxEnergyCluster_index]") + .Define("maxEnergyCluster_y", "CaloClusters.position.y[maxEnergyCluster_index]") + .Define("maxEnergyCluster_z", "CaloClusters.position.z[maxEnergyCluster_index]") + .Define("maxEnergyCluster_phi", "CaloNtupleizer::getCaloCluster_phi(CaloClusters)[maxEnergyCluster_index]") + .Define("maxEnergyCluster_theta", "CaloNtupleizer::getCaloCluster_theta(CaloClusters)[maxEnergyCluster_index]") + + .Define("cells_e", "Utils::as_vector(maxEnergyCluster_cells_energy)") + .Define("cells_theta", "Utils::as_vector(maxEnergyCluster_cells_theta)") + .Define("cells_phi", "Utils::as_vector(maxEnergyCluster_cells_phi)") + .Define("cells_radius", "Utils::as_vector(maxEnergyCluster_cells_radius)") + .Define("cells_layer", "Utils::as_vector(maxEnergyCluster_cells_layer)") + .Define("MVAVec", "WeaverUtils::get_weights(cells_e, cells_theta, cells_phi, cells_radius, cells_layer)") - # recast output .Define("Cluster_isPhoton", "WeaverUtils::get_weight(MVAVec, 0)") - .Define("Cluster_isPi0", "WeaverUtils::get_weight(MVAVec, 1)") + .Define("Cluster_isPi0", "WeaverUtils::get_weight(MVAVec, 1)") + ) return df2 @@ -36,6 +78,9 @@ def analysers(df): #Mandatory: output function, please make sure you return the branchlist as a python list def output(): branchList = [ - 'Cluster_isPhoton', 'Cluster_isPi0', + "maxEnergyCluster_index","minEnergyCluster_index","clusters_n","clusters_energy", + "maxEnergyCluster_cells_energy","maxEnergyCluster_cells_phi","maxEnergyCluster_cells_theta","maxEnergyCluster_cells_layer","maxEnergyCluster_cells_n","maxEnergyCluster_cells_radius", + "maxEnergyCluster_energy", "maxEnergyCluster_x", "maxEnergyCluster_y", "maxEnergyCluster_z", "maxEnergyCluster_phi", "maxEnergyCluster_theta", + "Cluster_isPhoton","Cluster_isPi0","cells_e" ] return branchList From b43ba489aa4573c2ab797947d67a162e897269e4 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 14:59:15 +0200 Subject: [PATCH 20/25] bug fix --- analyzers/dataframe/src/CaloNtupleizer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analyzers/dataframe/src/CaloNtupleizer.cc b/analyzers/dataframe/src/CaloNtupleizer.cc index daee5fdade..8ff0d021b7 100644 --- a/analyzers/dataframe/src/CaloNtupleizer.cc +++ b/analyzers/dataframe/src/CaloNtupleizer.cc @@ -11,7 +11,7 @@ namespace FCCAnalyses{ namespace CaloNtupleizer{ -std::unique_ptr m_decoder; +dd4hep::DDSegmentation::BitFieldCoder* m_decoder; void loadGeometry(std::string xmlGeometryPath, std::string readoutName){ dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance()); From 38168cadb7b5b6ef38c2064a79da7c7532b5e5c2 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 15:05:12 +0200 Subject: [PATCH 21/25] bug fix --- analyzers/dataframe/src/CaloNtupleizer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/analyzers/dataframe/src/CaloNtupleizer.cc b/analyzers/dataframe/src/CaloNtupleizer.cc index 8ff0d021b7..c162cf112b 100644 --- a/analyzers/dataframe/src/CaloNtupleizer.cc +++ b/analyzers/dataframe/src/CaloNtupleizer.cc @@ -29,7 +29,7 @@ ROOT::VecOps::RVec sel_layers::operator() (const R for (auto & p: in){ dd4hep::DDSegmentation::CellID cellId = p.cellID; int layer = m_decoder->get(cellId, "layer"); - if (layer>_min || layer<_max)res.emplace_back(p); + if (layer>_min && layer<_max)res.emplace_back(p); } return res; } From b1b766e1be4e32a030f2801644a2abc58523b8d7 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 15:38:46 +0200 Subject: [PATCH 22/25] add mva tuto --- .../tutorials/analysis_CaloClustering.py | 86 ---------------- .../FCCee/tutorials/analysis_tutorial_mva.py | 98 +++++++++++++++++++ 2 files changed, 98 insertions(+), 86 deletions(-) delete mode 100644 examples/FCCee/tutorials/analysis_CaloClustering.py create mode 100644 examples/FCCee/tutorials/analysis_tutorial_mva.py diff --git a/examples/FCCee/tutorials/analysis_CaloClustering.py b/examples/FCCee/tutorials/analysis_CaloClustering.py deleted file mode 100644 index b02fa15087..0000000000 --- a/examples/FCCee/tutorials/analysis_CaloClustering.py +++ /dev/null @@ -1,86 +0,0 @@ -from os import getenv -FCCDETECTORS=getenv('FCCDETECTORS') -geometryFile = FCCDETECTORS+"/Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.xml" -readoutName = "ECalBarrelPhiEta" - -import ROOT -ROOT.gInterpreter.Declare(""" -template -ROOT::VecOps::RVec myRange(ROOT::VecOps::RVec& vec, std::size_t begin, std::size_t end) -{ - ROOT::VecOps::RVec ret; - ret.reserve(end - begin); - for (auto i = begin; i < end; ++i) - ret.push_back(vec[i]); - return ret; -} -""") - - - -#Mandatory: RDFanalysis class where the use defines the operations on the TTree -class RDFanalysis(): - #__________________________________________________________ - #Mandatory: analysers funtion to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 - def analysers(df): - - from ROOT import WeaverUtils - from os import getenv - test_inputs_path = getenv('TEST_INPUT_DATA_DIR', '/eos/experiment/fcc/ee/tutorial/PNet_pi0Gamma/v1') - weaver = WeaverUtils.setup_weaver(test_inputs_path + '/fccee_pi_vs_gamma_v1.onnx', - test_inputs_path + '/preprocess_fccee_pi_vs_gamma_v1.json', - ('recocells_e', 'recocells_theta', 'recocells_phi', 'recocells_radius', 'recocells_layer')) - - - df2 = (df - - .Define("maxEnergyCluster_index","std::distance(CaloClusters.energy.begin(),std::max_element(CaloClusters.energy.begin(), CaloClusters.energy.end()))") - .Define("minEnergyCluster_index","std::distance(CaloClusters.energy.begin(),std::min_element(CaloClusters.energy.begin(), CaloClusters.energy.end()))") - .Define("clusters_n","CaloClusters.energy.size()") - .Define("clusters_energy","CaloClusters.energy") - - .Define("maxEnergyCluster_firstCell_index","CaloClusters[maxEnergyCluster_index].hits_begin") - .Define("maxEnergyCluster_lastCell_index" ,"CaloClusters[maxEnergyCluster_index].hits_end") - .Define("maxEnergyCluster_cells", "myRange(PositionedCaloClusterCells, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") - - .Define("maxEnergyCluster_cells_energy", "CaloNtupleizer::getCaloHit_energy(maxEnergyCluster_cells)" ) - .Define("maxEnergyCluster_cells_phi", "CaloNtupleizer::getCaloHit_phi(maxEnergyCluster_cells)" ) - .Define("maxEnergyCluster_cells_theta", "CaloNtupleizer::getCaloHit_theta(maxEnergyCluster_cells)" ) - .Define("maxEnergyCluster_cells_layer" , "CaloNtupleizer::getCaloHit_layer(maxEnergyCluster_cells)" ) - .Define("maxEnergyCluster_cells_n" , "maxEnergyCluster_cells.size()" ) - - .Define("maxEnergyCluster_cells_x", "myRange(PositionedCaloClusterCells.position.x, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") - .Define("maxEnergyCluster_cells_y", "myRange(PositionedCaloClusterCells.position.y, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") - .Define("maxEnergyCluster_cells_radius", "sqrt(pow(maxEnergyCluster_cells_x,2)+pow(maxEnergyCluster_cells_y,2))") - - .Define("maxEnergyCluster_energy", "CaloClusters.energy[maxEnergyCluster_index]") - .Define("maxEnergyCluster_x", "CaloClusters.position.x[maxEnergyCluster_index]") - .Define("maxEnergyCluster_y", "CaloClusters.position.y[maxEnergyCluster_index]") - .Define("maxEnergyCluster_z", "CaloClusters.position.z[maxEnergyCluster_index]") - .Define("maxEnergyCluster_phi", "CaloNtupleizer::getCaloCluster_phi(CaloClusters)[maxEnergyCluster_index]") - .Define("maxEnergyCluster_theta", "CaloNtupleizer::getCaloCluster_theta(CaloClusters)[maxEnergyCluster_index]") - - .Define("cells_e", "Utils::as_vector(maxEnergyCluster_cells_energy)") - .Define("cells_theta", "Utils::as_vector(maxEnergyCluster_cells_theta)") - .Define("cells_phi", "Utils::as_vector(maxEnergyCluster_cells_phi)") - .Define("cells_radius", "Utils::as_vector(maxEnergyCluster_cells_radius)") - .Define("cells_layer", "Utils::as_vector(maxEnergyCluster_cells_layer)") - - .Define("MVAVec", "WeaverUtils::get_weights(cells_e, cells_theta, cells_phi, cells_radius, cells_layer)") - - .Define("Cluster_isPhoton", "WeaverUtils::get_weight(MVAVec, 0)") - .Define("Cluster_isPi0", "WeaverUtils::get_weight(MVAVec, 1)") - - ) - return df2 - - #__________________________________________________________ - #Mandatory: output function, please make sure you return the branchlist as a python list - def output(): - branchList = [ - "maxEnergyCluster_index","minEnergyCluster_index","clusters_n","clusters_energy", - "maxEnergyCluster_cells_energy","maxEnergyCluster_cells_phi","maxEnergyCluster_cells_theta","maxEnergyCluster_cells_layer","maxEnergyCluster_cells_n","maxEnergyCluster_cells_radius", - "maxEnergyCluster_energy", "maxEnergyCluster_x", "maxEnergyCluster_y", "maxEnergyCluster_z", "maxEnergyCluster_phi", "maxEnergyCluster_theta", - "Cluster_isPhoton","Cluster_isPi0","cells_e" - ] - return branchList diff --git a/examples/FCCee/tutorials/analysis_tutorial_mva.py b/examples/FCCee/tutorials/analysis_tutorial_mva.py new file mode 100644 index 0000000000..3fd9bf1a9e --- /dev/null +++ b/examples/FCCee/tutorials/analysis_tutorial_mva.py @@ -0,0 +1,98 @@ +from os import getenv +FCCDETECTORS=getenv('FCCDETECTORS') +geometryFile = FCCDETECTORS+"/Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.xml" +readoutName = "ECalBarrelPhiEta" + +#Mandatory: RDFanalysis class where the use defines the operations on the TTree +class RDFanalysis(): + #__________________________________________________________ + #Mandatory: analysers funtion to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 + def analysers(df): + + #util function to get a sub-set of a vector + ROOT.gInterpreter.Declare(""" + template + ROOT::VecOps::RVec myRange(ROOT::VecOps::RVec& vec, std::size_t begin, std::size_t end){ + ROOT::VecOps::RVec ret; + ret.reserve(end - begin); + for (auto i = begin; i < end; ++i) + ret.push_back(vec[i]); + return ret; + } +""") + + #Get the weaver Deep learning + from ROOT import WeaverUtils + from os import getenv + test_inputs_path = getenv('TEST_INPUT_DATA_DIR', '/eos/experiment/fcc/ee/tutorial/PNet_pi0Gamma/v1') + weaver = WeaverUtils.setup_weaver(test_inputs_path + '/fccee_pi_vs_gamma_v1.onnx', + test_inputs_path + '/preprocess_fccee_pi_vs_gamma_v1.json', + ('recocells_e', 'recocells_theta', 'recocells_phi', 'recocells_radius', 'recocells_layer')) + + df2 = (df + #Define the index of the highest and lowest energetic cluster in the CaloClusters collection + .Define("maxEnergyCluster_index", "std::distance(CaloClusters.energy.begin(),std::max_element(CaloClusters.energy.begin(), CaloClusters.energy.end()))") + .Define("minEnergyCluster_index", "std::distance(CaloClusters.energy.begin(),std::min_element(CaloClusters.energy.begin(), CaloClusters.energy.end()))") + + #Define the number of clusters and their energies + .Define("clusters_n", "CaloClusters.energy.size()") + .Define("clusters_energy", "CaloClusters.energy") + + #For the max energetic cluster, define the index of the first and last cells in the PositionedCaloClusterCells collection + .Define("maxEnergyCluster_firstCell_index", "CaloClusters[maxEnergyCluster_index].hits_begin") + .Define("maxEnergyCluster_lastCell_index" , "CaloClusters[maxEnergyCluster_index].hits_end") + + #Using the first and last cells indices, build a sub-collection of cells. Now we have a collection of Cells from the highest energetic cluster + .Define("maxEnergyCluster_cells", "myRange(PositionedCaloClusterCells, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") + + #uncomment to filter cells that would not belong to the last two layers of the calorimeter (need comment the previous line) + #.Define("maxEnergyCluster_cellsFull", "myRange(PositionedCaloClusterCells, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") + #.Define("maxEnergyCluster_cells", ROOT.CaloNtupleizer.sel_layers(0, 10),["maxEnergyCluster_cellsFull"]) + + #Define energy, phi, theta, layer and n from the sub cell collection + .Define("maxEnergyCluster_cells_energy", "CaloNtupleizer::getCaloHit_energy(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_phi", "CaloNtupleizer::getCaloHit_phi(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_theta", "CaloNtupleizer::getCaloHit_theta(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_layer" , "CaloNtupleizer::getCaloHit_layer(maxEnergyCluster_cells)" ) + .Define("maxEnergyCluster_cells_n" , "maxEnergyCluster_cells.size()" ) + + #Define the x and y position to compute the radius + .Define("maxEnergyCluster_cells_x", "myRange(PositionedCaloClusterCells.position.x, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") + .Define("maxEnergyCluster_cells_y", "myRange(PositionedCaloClusterCells.position.y, maxEnergyCluster_firstCell_index, maxEnergyCluster_lastCell_index)") + .Define("maxEnergyCluster_cells_radius", "sqrt(pow(maxEnergyCluster_cells_x,2)+pow(maxEnergyCluster_cells_y,2))") + + #The Deep learning interface needs to use vectors (as we might want to evaluate several clusters in the same event) + .Define("cells_e", "Utils::as_vector(maxEnergyCluster_cells_energy)") + .Define("cells_theta", "Utils::as_vector(maxEnergyCluster_cells_theta)") + .Define("cells_phi", "Utils::as_vector(maxEnergyCluster_cells_phi)") + .Define("cells_radius", "Utils::as_vector(maxEnergyCluster_cells_radius)") + .Define("cells_layer", "Utils::as_vector(maxEnergyCluster_cells_layer)") + + #Call the inference + .Define("MVAVec", "WeaverUtils::get_weights(cells_e, cells_theta, cells_phi, cells_radius, cells_layer)") + + #The result is a vector (type photon or pi0) of vector (number of cluster) + .Define("Cluster_isPhoton", "WeaverUtils::get_weight(MVAVec, 0)") + .Define("Cluster_isPi0", "WeaverUtils::get_weight(MVAVec, 1)") + + ) + return df2 + + #__________________________________________________________ + #Mandatory: output function, please make sure you return the branchlist as a python list + def output(): + branchList = [ + "maxEnergyCluster_index", + "minEnergyCluster_index", + "clusters_n", + "clusters_energy", + "maxEnergyCluster_cells_energy", + "maxEnergyCluster_cells_phi", + "maxEnergyCluster_cells_theta", + "maxEnergyCluster_cells_layer", + "maxEnergyCluster_cells_n", + "maxEnergyCluster_cells_radius", + "Cluster_isPhoton", + "Cluster_isPi0" + ] + return branchList From bcacea91593fa0b0e6d403ec8f338eeab75d9047 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 15:47:46 +0200 Subject: [PATCH 23/25] add test files --- examples/FCCee/tutorials/analysis_tutorial_mva.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/examples/FCCee/tutorials/analysis_tutorial_mva.py b/examples/FCCee/tutorials/analysis_tutorial_mva.py index 3fd9bf1a9e..89df534274 100644 --- a/examples/FCCee/tutorials/analysis_tutorial_mva.py +++ b/examples/FCCee/tutorials/analysis_tutorial_mva.py @@ -1,3 +1,8 @@ +#test files of single photons and pi0 at 10GeV +#testFile='/eos/experiment/fcc/ee/tutorial/pi0GammaLAr2022/edm4hepFormat_smallSampleNotUsedForTraining/output_caloFullSim_10GeV_pdgId_22_noiseFalse.root' +#testFile='/eos/experiment/fcc/ee/tutorial/pi0GammaLAr2022/edm4hepFormat_smallSampleNotUsedForTraining/output_caloFullSim_10GeV_pdgId_111_noiseFalse.root' + +#Get the detector geometry from os import getenv FCCDETECTORS=getenv('FCCDETECTORS') geometryFile = FCCDETECTORS+"/Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.xml" From 316a8b6a52eb5d10018e4f87eee4640c116c3b76 Mon Sep 17 00:00:00 2001 From: Clement Helsens Date: Thu, 20 Oct 2022 16:54:51 +0200 Subject: [PATCH 24/25] add root import --- examples/FCCee/tutorials/analysis_tutorial_mva.py | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/FCCee/tutorials/analysis_tutorial_mva.py b/examples/FCCee/tutorials/analysis_tutorial_mva.py index 89df534274..326e816636 100644 --- a/examples/FCCee/tutorials/analysis_tutorial_mva.py +++ b/examples/FCCee/tutorials/analysis_tutorial_mva.py @@ -15,6 +15,7 @@ class RDFanalysis(): def analysers(df): #util function to get a sub-set of a vector + import ROOT ROOT.gInterpreter.Declare(""" template ROOT::VecOps::RVec myRange(ROOT::VecOps::RVec& vec, std::size_t begin, std::size_t end){ From 91025ade92f54408e749635b6e93d26858613b2e Mon Sep 17 00:00:00 2001 From: Michele Selvaggi Date: Mon, 31 Oct 2022 12:43:12 +0100 Subject: [PATCH 25/25] Jet Flavour data pre-processing and inference (#224) * added methods for computing tagging observables * add stage1/2 macros for producing jet ntuples * first commit README * added all stage script * commit inference example * polished readme * Remove commented code * Apply suggestions from code review Co-authored-by: Valentin Volkl --- .../FCCAnalyses/JetConstituentsUtils.h | 240 ++++ .../FCCAnalyses/ReconstructedParticle2Track.h | 37 + .../dataframe/src/JetConstituentsUtils.cc | 1027 ++++++++++++++++- .../src/ReconstructedParticle2Track.cc | 260 ++++- examples/FCCee/weaver/README.md | 323 ++++++ examples/FCCee/weaver/all_stages.py | 145 +++ examples/FCCee/weaver/analysis_inference.py | 263 +++++ examples/FCCee/weaver/stage1.py | 193 ++++ examples/FCCee/weaver/stage2.cpp | 535 +++++++++ 9 files changed, 2999 insertions(+), 24 deletions(-) create mode 100644 examples/FCCee/weaver/README.md create mode 100644 examples/FCCee/weaver/all_stages.py create mode 100644 examples/FCCee/weaver/analysis_inference.py create mode 100644 examples/FCCee/weaver/stage1.py create mode 100644 examples/FCCee/weaver/stage2.cpp diff --git a/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h b/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h index 7399f008f8..22909c8dfe 100644 --- a/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h +++ b/analyzers/dataframe/FCCAnalyses/JetConstituentsUtils.h @@ -3,6 +3,13 @@ #include "ROOT/RVec.hxx" #include "edm4hep/ReconstructedParticle.h" +#include "edm4hep/MCParticle.h" +#include "fastjet/JetDefinition.hh" + +#include "TMath.h" +#include "TVector3.h" +#include "TRotation.h" +#include "TLorentzVector.h" namespace FCCAnalyses { namespace JetConstituentsUtils { @@ -14,18 +21,251 @@ namespace FCCAnalyses { rv::RVec build_constituents(const rv::RVec&, const rv::RVec&); + rv::RVec build_constituents_cluster(const rv::RVec& rps, + const std::vector>& indices); + /// Retrieve the constituents of an indexed jet in event FCCAnalysesJetConstituents get_jet_constituents(const rv::RVec&, int); /// Retrieve the constituents of a collection of indexed jets in event rv::RVec get_constituents(const rv::RVec&, const rv::RVec&); + + //sorting jets + ROOT::VecOps::RVec jets_sorting_on_nconst(const rv::RVec&); + ROOT::VecOps::RVec jets_sorting_on_energy(const rv::RVec&); + + + rv::RVec get_Bz(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + rv::RVec get_pt(const rv::RVec&); rv::RVec get_e(const rv::RVec&); rv::RVec get_theta(const rv::RVec&); rv::RVec get_phi(const rv::RVec&); rv::RVec get_type(const rv::RVec&); rv::RVec get_charge(const rv::RVec&); + + //displacement + rv::RVec get_d0(const rv::RVec&, + const ROOT::VecOps::RVec&); + + rv::RVec get_z0(const rv::RVec& , + const ROOT::VecOps::RVec&); + + rv::RVec get_phi0(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_omega(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_tanLambda(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + + rv::RVec XPtoPar_dxy(const rv::RVec&, + const ROOT::VecOps::RVec&, + const TVector3&, + const float&); + rv::RVec XPtoPar_dz(const rv::RVec&, + const ROOT::VecOps::RVec&, + const TVector3&, + const float&); + rv::RVec XPtoPar_phi(const rv::RVec&, + const ROOT::VecOps::RVec&, + const TVector3&, + const float&); + rv::RVec XPtoPar_C(const rv::RVec&, + const ROOT::VecOps::RVec&, + const TVector3&, + const float&); + rv::RVec XPtoPar_ct(const rv::RVec&, + const ROOT::VecOps::RVec&, + const TVector3&, + const float&); + + //covariance matrix + //diagonal + rv::RVec get_omega_cov(const rv::RVec&, + const ROOT::VecOps::RVec&); + + rv::RVec get_d0_cov(const rv::RVec&, + const ROOT::VecOps::RVec& ); + + rv::RVec get_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_phi0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_tanlambda_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + //off-diag + rv::RVec get_d0_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_phi0_d0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_phi0_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_tanlambda_phi0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_tanlambda_d0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_tanlambda_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_omega_tanlambda_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_omega_phi0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_omega_d0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_omega_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + + rv::RVec get_dndx(const rv::RVec& jcs, + const rv::RVec& dNdx, + const rv::RVec& trackdata, + const rv::RVec JetsConstituents_isChargedHad); + + rv::RVec get_Sip2dVal(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_Sip2dVal_cluster(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + + rv::RVec get_Sip2dVal_clusterV(const rv::RVec& jets, + const rv::RVec& D0, + const rv::RVec& phi0, + const TVector3& V, + const float Bz); + + + rv::RVec get_Sip2dSig(const rv::RVec& Sip2dVals, + const rv::RVec& err2_D0); + + rv::RVec get_Sip3dVal(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + + rv::RVec get_Sip3dVal_cluster(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_Sip3dVal_clusterV(const rv::RVec& jets, + const rv::RVec& D0, + const rv::RVec& Z0, + const rv::RVec& phi0, + const TVector3& V, + const float Bz); + + rv::RVec get_Sip3dSig(const rv::RVec& Sip3dVals, + const rv::RVec& err2_D0, + const rv::RVec& err2_Z0); + + rv::RVec get_JetDistVal(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_JetDistVal_cluster(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks); + + rv::RVec get_JetDistVal_clusterV(const rv::RVec& jets, + const rv::RVec& jcs, + const rv::RVec& D0, + const rv::RVec& Z0, + const rv::RVec& phi0, + const TVector3& V, + const float Bz); + + rv::RVec get_JetDistSig(const rv::RVec& JetDistVal, + const rv::RVec& err2_D0, + const rv::RVec& err2_Z0); + + rv::RVec get_mtof(const rv::RVec& jcs, + const rv::RVec& track_L, + const rv::RVec& trackdata, + const rv::RVec& trackerhits, + const rv::RVec JetsConstituents_isChargedHad); + + + rv::RVec get_PIDs(const ROOT::VecOps::RVec< int > recin, + const ROOT::VecOps::RVec< int > mcin, + const rv::RVec& RecPart, + const rv::RVec& Particle, + const rv::RVec& Jets); + + rv::RVec get_PIDs_cluster(const ROOT::VecOps::RVec< int > recin, + const ROOT::VecOps::RVec< int > mcin, + const rv::RVec& RecPart, + const rv::RVec& Particle, + const std::vector>& indices); + + rv::RVec get_isMu(const rv::RVec&); + rv::RVec get_isEl(const rv::RVec&); + rv::RVec get_isChargedHad(const rv::RVec& PIDs, + const rv::RVec& jcs); + rv::RVec get_isGamma(const rv::RVec&); + rv::RVec get_isNeutralHad(const rv::RVec& PIDs, + const rv::RVec& jcs); + + //countings + int count_jets(rv::RVec jets); + rv::RVec count_consts(rv::RVec jets); + rv::RVec count_type(const rv::RVec& isType); + + + + rv::RVec get_erel(const rv::RVec& jets, + const rv::RVec& jcs); + rv::RVec get_erel_cluster(const rv::RVec& jets, + const rv::RVec& jcs); + + rv::RVec get_erel_log(const rv::RVec& jets, + const rv::RVec& jcs); + rv::RVec get_erel_log_cluster(const rv::RVec& jets, + const rv::RVec& jcs); + + rv::RVec get_thetarel(const rv::RVec& jets, + const rv::RVec& jcs); + rv::RVec get_thetarel_cluster(const rv::RVec& jets, + const rv::RVec& jcs); + + rv::RVec get_phirel(const rv::RVec& jets, + const rv::RVec& jcs); + rv::RVec get_phirel_cluster(const rv::RVec& jets, + const rv::RVec& jcs); + + //residues + rv::RVec compute_tlv_jets(const rv::RVec& jets); + rv::RVec sum_tlv_constituents(const rv::RVec& jets); + float InvariantMass(const TLorentzVector& tlv1, const TLorentzVector& tlv2); + rv::RVec compute_residue_energy(const rv::RVec& tlv_jet, + const rv::RVec& sum_tlv_jcs); + rv::RVec compute_residue_pt(const rv::RVec& tlv_jet, + const rv::RVec& sum_tlv_jcs); + rv::RVec compute_residue_phi(const rv::RVec& tlv_jet, + const rv::RVec& sum_tlv_jcs); + rv::RVec compute_residue_theta(const rv::RVec& tlv_jet, + const rv::RVec& sum_tlv_jcs); + rv::RVec compute_residue_px(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs); + rv::RVec compute_residue_py(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs); + rv::RVec compute_residue_pz(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs); + } // namespace JetConstituentsUtils } // namespace FCCAnalyses diff --git a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h index 91d6aff624..507276f137 100644 --- a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h +++ b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle2Track.h @@ -8,11 +8,48 @@ #include "ROOT/RVec.hxx" #include "edm4hep/ReconstructedParticleData.h" #include "edm4hep/TrackState.h" +#include +#include +#include +#include namespace FCCAnalyses{ namespace ReconstructedParticle2Track{ + //compute the magnetic field Bz + ROOT::VecOps::RVec getRP2TRK_Bz(const ROOT::VecOps::RVec& rps, + const ROOT::VecOps::RVec& tracks); //here computed for all particles passed + + float Bz(const ROOT::VecOps::RVec& rps, + const ROOT::VecOps::RVec& tracks); //here only computed for the first charged particle encountered + + + ROOT::VecOps::RVec XPtoPar_dxy(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& x, + const float& Bz); + + ROOT::VecOps::RVec XPtoPar_dz(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz); + + ROOT::VecOps::RVec XPtoPar_phi(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz); + + ROOT::VecOps::RVec XPtoPar_C(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz); + + ROOT::VecOps::RVec XPtoPar_ct(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz); + /// Return the D0 of a track to a reconstructed particle ROOT::VecOps::RVec getRP2TRK_D0 (ROOT::VecOps::RVec in, ROOT::VecOps::RVec tracks); diff --git a/analyzers/dataframe/src/JetConstituentsUtils.cc b/analyzers/dataframe/src/JetConstituentsUtils.cc index 13ea04580e..ab0e8d131c 100644 --- a/analyzers/dataframe/src/JetConstituentsUtils.cc +++ b/analyzers/dataframe/src/JetConstituentsUtils.cc @@ -1,5 +1,27 @@ #include "FCCAnalyses/JetConstituentsUtils.h" #include "FCCAnalyses/ReconstructedParticle.h" +#include "FCCAnalyses/ReconstructedParticle2Track.h" +#include "FCCAnalyses/ReconstructedParticle2MC.h" +#include "edm4hep/MCParticleData.h" +#include "edm4hep/Track.h" +#include "edm4hep/TrackerHitData.h" +#include "edm4hep/TrackData.h" +#include "edm4hep/ReconstructedParticleData.h" + + + +#include "FCCAnalyses/JetClusteringUtils.h" +//#include "FCCAnalyses/ExternalRecombiner.h" +#include "fastjet/JetDefinition.hh" +#include "fastjet/PseudoJet.hh" +#include "fastjet/Selector.hh" + +/* ************************* +//COMMENTS +1. Neutral particles (Clusters??) +2. units of measurement? + +************************ */ namespace FCCAnalyses { namespace JetConstituentsUtils { @@ -8,12 +30,30 @@ namespace FCCAnalyses { rv::RVec jcs; for (const auto& jet : jets) { auto& jc = jcs.emplace_back(); - for (auto it = jet.particles_begin; it < jet.particles_end; ++it) + float energy_jet = jet.energy; + float energy_const = 0; + for (auto it = jet.particles_begin; it < jet.particles_end; ++it) { jc.emplace_back(rps.at(it)); + energy_const += rps.at(it).energy; + } + } + return jcs; + } + + rv::RVec build_constituents_cluster(const rv::RVec& rps, + const std::vector>& indices) { + rv::RVec jcs; + for (const auto& jet_index : indices) { + FCCAnalysesJetConstituents jc; + for(const auto& const_index : jet_index) { + jc.push_back(rps.at(const_index)); + } + jcs.push_back(jc); } return jcs; } + FCCAnalysesJetConstituents get_jet_constituents(const rv::RVec& csts, int jet) { if (jet < 0) return FCCAnalysesJetConstituents(); @@ -39,14 +79,39 @@ namespace FCCAnalyses { return out; }; + + /// This function simply applies the 2 args functions per vector of Rec Particles to a vector of vectors of Rec Particles + auto cast_constituent_2 = [](const auto& jcs, const auto& coll, auto&& meth) { + rv::RVec out; + for (const auto& jc : jcs) { + out.emplace_back(meth(jc, coll)); + } + return out; + }; + + auto cast_constituent_4 = [](const auto& jcs, const auto& coll1, const auto& coll2, const auto& coll3,auto&& meth) { + rv::RVec out; + for (const auto& jc : jcs) { + out.emplace_back(meth(jc, coll1, coll2, coll3)); + } + return out; + }; + + + rv::RVec get_Bz(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Bz); + } + + rv::RVec get_pt(const rv::RVec& jcs) { return cast_constituent(jcs, ReconstructedParticle::get_pt); } - + rv::RVec get_e(const rv::RVec& jcs) { return cast_constituent(jcs, ReconstructedParticle::get_e); } - + rv::RVec get_theta(const rv::RVec& jcs) { return cast_constituent(jcs, ReconstructedParticle::get_theta); } @@ -62,5 +127,961 @@ namespace FCCAnalyses { rv::RVec get_charge(const rv::RVec& jcs) { return cast_constituent(jcs, ReconstructedParticle::get_charge); } + + + //sorting + ROOT::VecOps::RVec jets_sorting_on_nconst(const rv::RVec& jets){ + ROOT::VecOps::RVec nconst; + ROOT::VecOps::RVec out; + for (const auto& jet : jets) { + nconst.push_back(jet.particles_end - jet.particles_begin); + } + auto indices = ROOT::VecOps::Argsort(nconst); + for (int index = 0; index < jets.size(); ++index) { + out.push_back( jets.at( indices.at(indices.size()-1-index) ) ); + } + return out; + } + + ROOT::VecOps::RVec jets_sorting_on_energy(const rv::RVec& jets){ + ROOT::VecOps::RVec energy; + ROOT::VecOps::RVec out; + for (const auto& jet : jets) { + energy.push_back(jet.energy); + } + auto indices = ROOT::VecOps::Argsort(energy); + for (int index = 0; index < jets.size(); ++index) { + out.push_back( jets.at( indices.at(indices.size()-1-index) ) ); + } + return out; + } + + //displacement (wrt (0,0,0)) + rv::RVec get_d0(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0); + } + + rv::RVec get_z0(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Z0); + } + + rv::RVec get_phi0(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi); + } + + rv::RVec get_omega(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_omega); + } + + rv::RVec get_tanLambda(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_tanLambda); + } + + + rv::RVec XPtoPar_dxy(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_dxy); + } + + rv::RVec XPtoPar_dz(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_dz); + } + + rv::RVec XPtoPar_phi(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_phi); + } + + rv::RVec XPtoPar_C(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_C); + } + + rv::RVec XPtoPar_ct(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + return cast_constituent_4(jcs, tracks, V, Bz, ReconstructedParticle2Track::XPtoPar_ct); + } + + + //Covariance matrix elements of tracks parameters + //diagonal + rv::RVec get_omega_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_omega_cov); + } + + rv::RVec get_d0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0_cov); + } + + rv::RVec get_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Z0_cov); + } + + rv::RVec get_phi0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi_cov); + } + + rv::RVec get_tanlambda_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_tanLambda_cov); + } + //off-diagonal + rv::RVec get_d0_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_z0_cov); + } + + rv::RVec get_phi0_d0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_phi0_cov); + } + + rv::RVec get_phi0_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi0_z0_cov); + } + + rv::RVec get_tanlambda_phi0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi0_tanlambda_cov); + } + + rv::RVec get_tanlambda_d0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_tanlambda_cov); + } + + rv::RVec get_tanlambda_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_z0_tanlambda_cov); + } + + rv::RVec get_omega_tanlambda_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_omega_tanlambda_cov); + } + + rv::RVec get_omega_phi0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi0_omega_cov); + } + + rv::RVec get_omega_d0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_d0_omega_cov); + } + + rv::RVec get_omega_z0_cov(const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + return cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_omega_z0_cov); + } + + + //neutrals are set to 0; muons and electrons are also set to 0; + // only charged hads are considered (mtof used to disctriminate charged kaons and pions) + rv::RVec get_dndx(const rv::RVec& jcs, + const rv::RVec& dNdx, //ETrackFlow_2 + const rv::RVec& trackdata, //Eflowtrack + const rv::RVec JetsConstituents_isChargedHad) { + rv::RVec out; + for(int i = 0; i < jcs.size(); ++i){ + FCCAnalysesJetConstituents ct = jcs.at(i); + FCCAnalysesJetConstituentsData isChargedHad = JetsConstituents_isChargedHad.at(i); + FCCAnalysesJetConstituentsData tmp; + for(int j = 0; j < ct.size(); ++j) { + if (ct.at(j).tracks_begin < trackdata.size() && (int)isChargedHad.at(j) == 1) { + tmp.push_back( dNdx.at( trackdata.at(ct.at(j).tracks_begin).dxQuantities_begin).value ); + } else { + tmp.push_back(0.); + } + } + out.push_back(tmp); + } + return out; + } + + rv::RVec get_Sip2dVal(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + rv::RVec out; + rv::RVec D0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0); + rv::RVec phi0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi); + + + for(int i = 0; i < jets.size(); ++i){ + TVector2 p(jets[i].momentum.x, jets[i].momentum.y); + FCCAnalysesJetConstituentsData cprojs; + for (int j = 0; j < jcs[i].size(); ++j) { + if (D0.at(i).at(j) != -9) { + TVector2 d0( - D0.at(i).at(j) * TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j) * TMath::Cos(phi0.at(i).at(j)) ); + cprojs.push_back( TMath::Sign(1, d0*p) * fabs(D0.at(i).at(j)) ); + } else { + cprojs.push_back(-9.); + } + } + out.push_back(cprojs); + } + return out; + } + + rv::RVec get_Sip2dVal_cluster(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + rv::RVec out; + rv::RVec D0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0); + rv::RVec phi0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi); + + for(int i = 0; i < jets.size(); ++i){ + TVector2 p(jets[i].px(), jets[i].py()); + FCCAnalysesJetConstituentsData cprojs; + for (int j = 0; j < jcs[i].size(); ++j) { + if (D0.at(i).at(j) != -9) { + TVector2 d0( - D0.at(i).at(j) * TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j) * TMath::Cos(phi0.at(i).at(j)) ); + cprojs.push_back( TMath::Sign(1, d0*p) * fabs(D0.at(i).at(j)) ); + } else { + cprojs.push_back(-9.); + } + } + out.push_back(cprojs); + } + return out; + } + + + rv::RVec get_Sip2dVal_clusterV(const rv::RVec& jets, + const rv::RVec& D0, + const rv::RVec& phi0, + const TVector3& V, + const float Bz) { + rv::RVec out; + + for(int i = 0; i < jets.size(); ++i){ + TVector2 p(jets[i].px(), jets[i].py()); + FCCAnalysesJetConstituentsData cprojs; + for (int j = 0; j < D0[i].size(); ++j) { + if (D0.at(i).at(j) != -9) { + TVector2 d0( - D0.at(i).at(j) * TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j) * TMath::Cos(phi0.at(i).at(j)) ); + cprojs.push_back( TMath::Sign(1, d0*p) * fabs(D0.at(i).at(j)) ); + } else { + cprojs.push_back(-9.); + } + } + out.push_back(cprojs); + } + return out; + } + + + /// The functions get_Sip2dSig and get_Sip2dVal can be made independent; + /// I passed to the former the result of the latter, avoiding the recomputation + rv::RVec get_Sip2dSig(const rv::RVec& Sip2dVals, + const rv::RVec& err2_D0) { + rv::RVec out; + for(int i = 0; i < Sip2dVals.size(); ++i) { + FCCAnalysesJetConstituentsData s; + for(int j = 0; j < Sip2dVals.at(i).size(); ++j) { + if(err2_D0.at(i).at(j) > 0) { + s.push_back( Sip2dVals.at(i).at(j)/std::sqrt(err2_D0.at(i).at(j)) ); + } else { + s.push_back(-9); + } + } + out.push_back(s); + } + return out; + } + + + rv::RVec get_Sip3dVal(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + rv::RVec out; + rv::RVec D0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0); + rv::RVec Z0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Z0); + rv::RVec phi0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi); + + for(int i = 0; i < jets.size(); ++i){ + TVector3 p(jets[i].momentum.x, jets[i].momentum.y, jets[i].momentum.z); + FCCAnalysesJetConstituentsData cprojs; + for (int j = 0; j < jcs[i].size(); ++j){ + if(D0.at(i).at(j) != -9) { + TVector3 d( - D0.at(i).at(j) * TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j)*TMath::Cos(phi0.at(i).at(j)), Z0.at(i).at(j) ); + cprojs.push_back( TMath::Sign(1, d*p) * fabs( sqrt( D0.at(i).at(j)*D0.at(i).at(j) + Z0.at(i).at(j)*Z0.at(i).at(j) ) ) ); + } else { + cprojs.push_back(-9); + } + } + out.push_back(cprojs); + } + return out; + } + + rv::RVec get_Sip3dVal_cluster(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + rv::RVec out; + rv::RVec D0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0); + rv::RVec Z0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Z0); + rv::RVec phi0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi); + + for(int i = 0; i < jets.size(); ++i){ + TVector3 p(jets[i].px(), jets[i].py(), jets[i].pz()); + FCCAnalysesJetConstituentsData cprojs; + for (int j = 0; j < jcs[i].size(); ++j){ + if(D0.at(i).at(j) != -9) { + TVector3 d( - D0.at(i).at(j) * TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j)*TMath::Cos(phi0.at(i).at(j)), Z0.at(i).at(j) ); + cprojs.push_back( TMath::Sign(1, d*p) * fabs( sqrt( D0.at(i).at(j)*D0.at(i).at(j) + Z0.at(i).at(j)*Z0.at(i).at(j) ) ) ); + } else { + cprojs.push_back(-9); + } + } + out.push_back(cprojs); + } + return out; + } + + + rv::RVec get_Sip3dVal_clusterV(const rv::RVec& jets, + const rv::RVec& D0, + const rv::RVec& Z0, + const rv::RVec& phi0, + const TVector3& V, + const float Bz) { + rv::RVec out; + + for(int i = 0; i < jets.size(); ++i){ + TVector3 p(jets[i].px(), jets[i].py(), jets[i].pz()); + FCCAnalysesJetConstituentsData cprojs; + for (int j = 0; j < D0[i].size(); ++j){ + if(D0.at(i).at(j) != -9) { + TVector3 d( - D0.at(i).at(j) * TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j)*TMath::Cos(phi0.at(i).at(j)), Z0.at(i).at(j) ); + cprojs.push_back( TMath::Sign(1, d*p) * fabs( sqrt( D0.at(i).at(j)*D0.at(i).at(j) + Z0.at(i).at(j)*Z0.at(i).at(j) ) ) ); + } else { + cprojs.push_back(-9); + } + } + out.push_back(cprojs); + } + return out; + } + + + rv::RVec get_Sip3dSig(const rv::RVec& Sip3dVals, + const rv::RVec& err2_D0, + const rv::RVec& err2_Z0) { + rv::RVec out; + for(int i = 0; i < Sip3dVals.size(); ++i) { + FCCAnalysesJetConstituentsData s; + for(int j = 0; j < Sip3dVals.at(i).size(); ++j) { + if (err2_D0.at(i).at(j) > 0.) { + s.push_back( Sip3dVals.at(i).at(j)/sqrt( err2_D0.at(i).at(j) + err2_Z0.at(i).at(j) ) ); + } else { + s.push_back(-9); + } + } + out.push_back(s); + } + return out; + } + + + rv::RVec get_JetDistVal(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + rv::RVec out; + rv::RVec D0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0); + rv::RVec Z0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Z0); + rv::RVec phi0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi); + for(int i = 0; i < jets.size(); ++i){ + FCCAnalysesJetConstituentsData tmp; + TVector3 p_jet(jets[i].momentum.x, jets[i].momentum.y, jets[i].momentum.z); + FCCAnalysesJetConstituents ct = jcs.at(i); + for(int j = 0; j < ct.size(); ++j) { + if (D0.at(i).at(j) != -9) { + TVector3 d( - D0.at(i).at(j)* TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j)*TMath::Cos(phi0.at(i).at(j)), Z0.at(i).at(j) ); + TVector3 p_ct(ct[j].momentum.x, ct[j].momentum.y, ct[j].momentum.z); + TVector3 r_jet(0.0, 0.0, 0.0); + TVector3 n = p_ct.Cross(p_jet).Unit(); //What if they are parallel? + tmp.push_back( n.Dot(d - r_jet) ); + } else { + tmp.push_back(-9); + } + } + out.push_back(tmp); + } + return out; + } + + + rv::RVec get_JetDistVal_cluster(const rv::RVec& jets, + const rv::RVec& jcs, + const ROOT::VecOps::RVec& tracks) { + rv::RVec out; + rv::RVec D0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_D0); + rv::RVec Z0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_Z0); + rv::RVec phi0 = cast_constituent_2(jcs, tracks, ReconstructedParticle2Track::getRP2TRK_phi); + for(int i = 0; i < jets.size(); ++i){ + FCCAnalysesJetConstituentsData tmp; + TVector3 p_jet(jets[i].px(), jets[i].py(), jets[i].pz()); + FCCAnalysesJetConstituents ct = jcs.at(i); + for(int j = 0; j < ct.size(); ++j) { + if (D0.at(i).at(j) != -9) { + TVector3 d( - D0.at(i).at(j)* TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j)*TMath::Cos(phi0.at(i).at(j)), Z0.at(i).at(j) ); + TVector3 p_ct(ct[j].momentum.x, ct[j].momentum.y, ct[j].momentum.z); + TVector3 r_jet(0.0, 0.0, 0.0); + TVector3 n = p_ct.Cross(p_jet).Unit(); //What if they are parallel? + tmp.push_back( n.Dot(d - r_jet) ); + } else { + tmp.push_back(-9); + } + } + out.push_back(tmp); + } + return out; + } + + rv::RVec get_JetDistVal_clusterV(const rv::RVec& jets, + const rv::RVec& jcs, + const rv::RVec& D0, + const rv::RVec& Z0, + const rv::RVec& phi0, + const TVector3& V, + const float Bz) { + rv::RVec out; + + for(int i = 0; i < jets.size(); ++i){ + FCCAnalysesJetConstituentsData tmp; + TVector3 p_jet(jets[i].px(), jets[i].py(), jets[i].pz()); + FCCAnalysesJetConstituents ct = jcs.at(i); + for(int j = 0; j < ct.size(); ++j) { + if (D0.at(i).at(j) != -9) { + TVector3 d( - D0.at(i).at(j)* TMath::Sin(phi0.at(i).at(j)) , D0.at(i).at(j)*TMath::Cos(phi0.at(i).at(j)), Z0.at(i).at(j) ); + TVector3 p_ct(ct[j].momentum.x, ct[j].momentum.y, ct[j].momentum.z); + TVector3 r_jet(0.0, 0.0, 0.0); + TVector3 n = p_ct.Cross(p_jet).Unit(); //What if they are parallel? + tmp.push_back( n.Dot(d - r_jet) ); + } else { + tmp.push_back(-9); + } + } + out.push_back(tmp); + } + return out; + } + + + rv::RVec get_JetDistSig(const rv::RVec& JetDistVal, + const rv::RVec& err2_D0, + const rv::RVec& err2_Z0) { + rv::RVec out; + for(int i = 0; i < JetDistVal.size(); ++i) { + FCCAnalysesJetConstituentsData tmp; + for(int j = 0; j < JetDistVal.at(i).size(); ++j) { + if (err2_D0.at(i).at(j) > 0) { + float err3d = std::sqrt(err2_D0.at(i).at(j) + err2_Z0.at(i).at(j)); + float jetdistsig = JetDistVal.at(i).at(j) / err3d; + tmp.push_back(jetdistsig); + } else { + tmp.push_back(-9.); + } + } + out.push_back(tmp); + } + return out; + } + + + //we measure L, tof; mtof in GeV + //neutrals are set to 0; muons and electrons are set to their mass; + // only charged hads are considered (mtof used to disctriminate charged kaons and pions) + rv::RVec get_mtof(const rv::RVec& jcs, + const rv::RVec& track_L, + const rv::RVec& trackdata, + const rv::RVec& trackerhits, + const rv::RVec JetsConstituents_Pids) { + rv::RVec out; + for(int i = 0; i < jcs.size(); ++i){ + FCCAnalysesJetConstituents ct = jcs.at(i); + FCCAnalysesJetConstituentsData pids = JetsConstituents_Pids.at(i); + FCCAnalysesJetConstituentsData tmp; + for(int j = 0; j < ct.size(); ++j) { + //if(ct.at(j).tracks_begin > 0 && ct.at(j).tracks_begin < trackdata.size()) { //??????????? CHECK!!! + if (ct.at(j).tracks_begin < trackdata.size()) { + if( abs(pids.at(j)) == 11) { + tmp.push_back(0.00051099895); + } else if (abs(pids.at(j)) == 13) { + tmp.push_back(105.65837); + } else { + float Tin = trackerhits.at(trackdata.at(ct.at(j).tracks_begin).trackerHits_begin).time; + float Tout = trackerhits.at(trackdata.at(ct.at(j).tracks_begin).trackerHits_end-1).time; //one track and two hits per recon. particle are assumed + float tof = (Tout - Tin); + float L = track_L.at(ct.at(j).tracks_begin) * 0.001; + //std::cout << "tof: " << tof << " - L: " << L << << std::endl; + float beta = L/(tof * 2.99792458e+8) ; + float p = std::sqrt( ct.at(j).momentum.x*ct.at(j).momentum.x + ct.at(j).momentum.y*ct.at(j).momentum.y + ct.at(j).momentum.z * ct.at(j).momentum.z ); + //std::cout << "tof: " << tof << " - L: " << L << " - beta: " << beta << " - momentum: " << p << " - mtof: " << p * std::sqrt(1/(beta*beta)-1) << std::endl; + if (beta < 1. && beta > 0.) { + tmp.push_back( p * std::sqrt(1/(beta*beta)-1) ); + } else { + tmp.push_back(0.13957039); + } + } + } else { + //float E = ct.at(j).energy; + //tmp.pushback( E * std::sqrt(1-beta*beta) ); + tmp.push_back(0.); + } + } + out.push_back(tmp); + } + return out; + } + + + //kinematics const/jet + rv::RVec get_erel_log(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + float e_jet = jets.at(i).energy; + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + float val = (e_jet > 0.) ? jc.energy / e_jet : 1.; + float erel_log = float(std::log10(val)); + jet_csts.emplace_back(erel_log); + } + } + return out; + } + + rv::RVec get_erel_log_cluster(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + float e_jet = jets.at(i).E(); + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + float val = (e_jet > 0.) ? jc.energy / e_jet : 1.; + float erel_log = float(std::log10(val)); + jet_csts.emplace_back(erel_log); + } + } + return out; + } + + + rv::RVec get_erel(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + double e_jet = jets.at(i).energy; + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + float val = (e_jet > 0.) ? jc.energy / e_jet : 1.; + float erel = val; + jet_csts.emplace_back(erel); + } + } + return out; + } + + rv::RVec get_erel_cluster(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + double e_jet = jets.at(i).E(); + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + float val = (e_jet > 0.) ? jc.energy / e_jet : 1.; + float erel = val; + jet_csts.emplace_back(erel); + } + } + return out; + } + + rv::RVec get_thetarel(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + TLorentzVector tlv_jet; + tlv_jet.SetXYZM(jets.at(i).momentum.x, jets.at(i).momentum.y, jets.at(i).momentum.z, jets.at(i).mass); + float theta_jet = tlv_jet.Theta(); + float phi_jet = tlv_jet.Phi(); + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + TLorentzVector tlv_const; + tlv_const.SetXYZM(jc.momentum.x, jc.momentum.y, jc.momentum.z, jc.mass); + TVector3 v_const = tlv_const.Vect(); + v_const.RotateZ(-phi_jet); + v_const.RotateY(-theta_jet); + float theta_rel = v_const.Theta(); + jet_csts.emplace_back(theta_rel); + } + } + return out; + } + + rv::RVec get_thetarel_cluster(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + TLorentzVector tlv_jet; + tlv_jet.SetXYZM(jets.at(i).px(), jets.at(i).py(), jets.at(i).pz(), jets.at(i).E()); + float theta_jet = tlv_jet.Theta(); + float phi_jet = tlv_jet.Phi(); + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + TLorentzVector tlv_const; + tlv_const.SetXYZM(jc.momentum.x, jc.momentum.y, jc.momentum.z, jc.energy); + TVector3 v_const = tlv_const.Vect(); + v_const.RotateZ(-phi_jet); + v_const.RotateY(-theta_jet); + float theta_rel = v_const.Theta(); + jet_csts.emplace_back(theta_rel); + } + } + return out; + } + + + rv::RVec get_phirel(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + TLorentzVector tlv_jet; + tlv_jet.SetXYZM(jets.at(i).momentum.x, jets.at(i).momentum.y, jets.at(i).momentum.z, jets.at(i).mass); + float theta_jet = tlv_jet.Theta(); + float phi_jet = tlv_jet.Phi(); + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + TLorentzVector tlv_const; + tlv_const.SetXYZM(jc.momentum.x, jc.momentum.y, jc.momentum.z, jc.mass); + TVector3 v_const = tlv_const.Vect(); + v_const.RotateZ(-phi_jet); + v_const.RotateY(-theta_jet); + float phi_rel = v_const.Phi(); + jet_csts.emplace_back(phi_rel); + } + } + return out; + } + + rv::RVec get_phirel_cluster(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { + auto& jet_csts = out.emplace_back(); + TLorentzVector tlv_jet; + tlv_jet.SetXYZM(jets.at(i).px(), jets.at(i).py(), jets.at(i).pz(), jets.at(i).E()); + float theta_jet = tlv_jet.Theta(); + float phi_jet = tlv_jet.Phi(); + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { + TLorentzVector tlv_const; + tlv_const.SetXYZM(jc.momentum.x, jc.momentum.y, jc.momentum.z, jc.mass); + TVector3 v_const = tlv_const.Vect(); + v_const.RotateZ(-phi_jet); + v_const.RotateY(-theta_jet); + float phi_rel = v_const.Phi(); + jet_csts.emplace_back(phi_rel); + } + } + return out; + } + + + + //Identification + + rv::RVec get_PIDs(const ROOT::VecOps::RVec< int > recin, + const ROOT::VecOps::RVec< int > mcin, + const rv::RVec& RecPart, + const rv::RVec& Particle, + const rv::RVec& jets) { + rv::RVec out; + FCCAnalysesJetConstituentsData PIDs = FCCAnalyses::ReconstructedParticle2MC::getRP2MC_pdg(recin, mcin, RecPart, Particle); + + for (const auto& jet : jets) { + FCCAnalysesJetConstituentsData tmp; + for (auto it = jet.particles_begin; it < jet.particles_end; ++it) { + tmp.push_back(PIDs.at(it)); + } + out.push_back(tmp); + } + return out; + } + + rv::RVec get_PIDs_cluster(const ROOT::VecOps::RVec< int > recin, + const ROOT::VecOps::RVec< int > mcin, + //const rv::RVec& jcs, + const rv::RVec& RecPart, + const rv::RVec& Particle, + const std::vector>& indices) { + rv::RVec out; + FCCAnalysesJetConstituentsData PIDs = FCCAnalyses::ReconstructedParticle2MC::getRP2MC_pdg(recin, mcin, RecPart, Particle); + + for (const auto& jet_index : indices) { + FCCAnalysesJetConstituentsData tmp; + for (const auto& const_index : jet_index) { + tmp.push_back(PIDs.at(const_index)); + } + out.push_back(tmp); + } + return out; + } + + + rv::RVec get_isMu(const rv::RVec& PIDs) { + rv::RVec out; + for(int i = 0; i < PIDs.size(); ++i) { + FCCAnalysesJetConstituentsData is_Mu; + for (int j = 0; j < PIDs.at(i).size(); ++j) { + if ( abs(PIDs.at(i).at(j)) == 13) { + is_Mu.push_back(1.); + } else { + is_Mu.push_back(0.); + } + } + out.push_back(is_Mu); + } + return out; + } + + + rv::RVec get_isEl(const rv::RVec& PIDs) { + rv::RVec out; + for(int i = 0; i < PIDs.size(); ++i) { + FCCAnalysesJetConstituentsData is_El; + FCCAnalysesJetConstituentsData pids = PIDs.at(i); + for (int j = 0; j < pids.size(); ++j) { + if ( abs(pids.at(j)) == 11) { + is_El.push_back(1.); + } else { + is_El.push_back(0.); + } + } + out.push_back(is_El); + } + return out; + } + + rv::RVec get_isChargedHad(const rv::RVec& PIDs, + const rv::RVec& jcs) { + rv::RVec out; + for(int i = 0; i < PIDs.size(); ++i) { + FCCAnalysesJetConstituentsData is_ChargedHad; + FCCAnalysesJetConstituents ct = jcs.at(i); + FCCAnalysesJetConstituentsData pids = PIDs.at(i); + for (int j = 0; j < pids.size(); ++j) { + if (ct.at(j).charge != 0 && abs(pids.at(j)) != 11 && abs(pids.at(j)) != 13) { + is_ChargedHad.push_back(1.); + } else { + is_ChargedHad.push_back(0.); + } + } + out.push_back(is_ChargedHad); + } + return out; + } + + rv::RVec get_isGamma(const rv::RVec& PIDs) { + rv::RVec out; + for(int i = 0; i < PIDs.size(); ++i) { + FCCAnalysesJetConstituentsData is_Gamma; + FCCAnalysesJetConstituentsData pids = PIDs.at(i); + for (int j = 0; j < pids.size(); ++j) { + if ( abs(pids.at(j)) == 22) { + is_Gamma.push_back(1.); + } else { + is_Gamma.push_back(0.); + } + } + out.push_back(is_Gamma); + } + return out; + } + + rv::RVec get_isNeutralHad(const rv::RVec& PIDs, + const rv::RVec& jcs) { + rv::RVec out; + for(int i = 0; i < PIDs.size(); ++i) { + FCCAnalysesJetConstituentsData is_NeutralHad; + FCCAnalysesJetConstituents ct = jcs.at(i); + FCCAnalysesJetConstituentsData pids = PIDs.at(i); + for (int j = 0; j < pids.size(); ++j) { + if (ct.at(j).charge == 0 && abs(pids.at(j)) != 22 ) + is_NeutralHad.push_back(1.); + else + is_NeutralHad.push_back(0.); + } + out.push_back(is_NeutralHad); + } + return out; + } + + + //countings + int count_jets(rv::RVec jets) { + return jets.size(); + } + + rv::RVec count_consts(rv::RVec jets) { + rv::RVec out; + for(int i = 0; i < jets.size(); ++i) { + out.push_back(jets.at(i).size()); + } + return out; + } + + rv::RVec count_type(const rv::RVec& isType) { + rv::RVec out; + for(int i = 0; i < isType.size(); ++i){ + int count = 0; + rv::RVec istype = isType.at(i); + for(int j = 0; j < istype.size(); ++j){ + if( (int)(istype.at(j)) == 1) count++; + } + out.push_back(count); + } + return out; + } + + //compute residues + rv::RVec compute_tlv_jets(const rv::RVec& jets) { + rv::RVec out; + for(const auto& jet : jets) { + TLorentzVector tlv_jet; + tlv_jet.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E()); + out.push_back(tlv_jet); + } + return out; + } + + rv::RVec sum_tlv_constituents(const rv::RVec& jets) { + rv::RVec out; + for (int i = 0; i < jets.size(); ++i) { + TLorentzVector sum_tlv; // initialized by (0., 0., 0., 0.) + FCCAnalysesJetConstituents jcs = jets.at(i); + for (const auto& jc : jcs) { + TLorentzVector tlv; + tlv.SetPxPyPzE(jc.momentum.x, jc.momentum.y, jc.momentum.z, jc.energy); + sum_tlv += tlv; + } + out.push_back(sum_tlv); + } + return out; + } + + float InvariantMass(const TLorentzVector& tlv1, const TLorentzVector& tlv2) { + float E = tlv1.E() + tlv2.E(); + float px = tlv1.Px() + tlv2.Px(); + float py = tlv1.Py() + tlv2.Py(); + float pz = tlv1.Pz() + tlv2.Pz(); + return std::sqrt(E*E - px*px - py*py - pz*pz); + } + + rv::RVec compute_residue_energy(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs) { + rv::RVec out; + for(int i = 0; i < tlv_jet.size(); ++i) { + float de = (sum_tlv_jcs.at(i).E() - tlv_jet.at(i).E())/tlv_jet.at(i).E(); + out.push_back(de); + } + return out; + } + + rv::RVec compute_residue_px(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs) { + rv::RVec out; + for(int i = 0; i < tlv_jet.size(); ++i) { + float dpx = (sum_tlv_jcs.at(i).Px() - tlv_jet.at(i).Px())/tlv_jet.at(i).Px(); + out.push_back(dpx); + } + return out; + } + + rv::RVec compute_residue_py(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs) { + rv::RVec out; + for(int i = 0; i < tlv_jet.size(); ++i) { + float dpy = (sum_tlv_jcs.at(i).Py() - tlv_jet.at(i).Py())/tlv_jet.at(i).Py(); + out.push_back(dpy); + } + return out; + } + + rv::RVec compute_residue_pz(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs) { + rv::RVec out; + for(int i = 0; i < tlv_jet.size(); ++i) { + float dpz = (sum_tlv_jcs.at(i).Pz() - tlv_jet.at(i).Pz())/tlv_jet.at(i).Pz(); + out.push_back(dpz); + } + return out; + } + + rv::RVec compute_residue_pt(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs) { + rv::RVec out; + for(int i = 0; i < tlv_jet.size(); ++i) { + double pt_jet = std::sqrt( tlv_jet.at(i).Px()*tlv_jet.at(i).Px() + tlv_jet.at(i).Py()*tlv_jet.at(i).Py() ); + double pt_jcs = std::sqrt( sum_tlv_jcs.at(i).Px()*sum_tlv_jcs.at(i).Px() + sum_tlv_jcs.at(i).Py()*sum_tlv_jcs.at(i).Py() ); + double dpt = ( pt_jcs - pt_jet)/pt_jet; + out.push_back(dpt); + } + return out; + } + + rv::RVec compute_residue_phi(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs) { + rv::RVec out; + for(int i = 0; i < tlv_jet.size(); ++i) { + double phi_jet = tlv_jet.at(i).Phi(); + double phi_jcs = sum_tlv_jcs.at(i).Phi(); + double dphi = (phi_jcs - phi_jet)/phi_jet; + out.push_back(dphi); + } + return out; + } + + rv::RVec compute_residue_theta(const rv::RVec& tlv_jet, const rv::RVec& sum_tlv_jcs) { + rv::RVec out; + for(int i = 0; i < tlv_jet.size(); ++i) { + double theta_jet = tlv_jet.at(i).Theta(); + double theta_jcs = sum_tlv_jcs.at(i).Theta(); + double dtheta = (theta_jcs - theta_jet)/theta_jet; + out.push_back(dtheta); + } + return out; + } + + } // namespace JetConstituentsUtils } // namespace FCCAnalyses diff --git a/analyzers/dataframe/src/ReconstructedParticle2Track.cc b/analyzers/dataframe/src/ReconstructedParticle2Track.cc index fc8a7b5a64..38de350179 100644 --- a/analyzers/dataframe/src/ReconstructedParticle2Track.cc +++ b/analyzers/dataframe/src/ReconstructedParticle2Track.cc @@ -4,6 +4,224 @@ namespace FCCAnalyses{ namespace ReconstructedParticle2Track{ + ROOT::VecOps::RVec getRP2TRK_Bz(const ROOT::VecOps::RVec& rps, const ROOT::VecOps::RVec& tracks) { + const double c_light = 2.99792458e8; + const double a = c_light * 1e3 * 1e-15; //[omega] = 1/mm + ROOT::VecOps::RVec out; + + for(auto & p: rps) { + if(p.tracks_begin < tracks.size()) { + double pt= sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y); + double Bz= tracks.at(p.tracks_begin).omega / a * pt * std::copysign(1.0, p.charge); + out.push_back(Bz); + } else { + out.push_back(-9.); + } + } + return out; + } + + float Bz(const ROOT::VecOps::RVec& rps, const ROOT::VecOps::RVec& tracks) { + const double c_light = 2.99792458e8;// speed of light m/sec; + const double a = c_light * 1e3 * 1e-15; //[omega] = 1/mm + + double Bz = -9; + + for(auto & p: rps) { + if(p.tracks_begin < tracks.size()) { + double pt= sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y); + Bz= tracks.at(p.tracks_begin).omega / a * pt * std::copysign(1.0, p.charge); + } + } + return Bz; + } + + ROOT::VecOps::RVec XPtoPar_dxy(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + const double cSpeed = 2.99792458e8 * 1.0e-9; + + ROOT::VecOps::RVec out; + + for (const auto & rp: in) { + + if( rp.tracks_begin < tracks.size()) { + + float D0_wrt0 = tracks.at(rp.tracks_begin).D0; + float Z0_wrt0 = tracks.at(rp.tracks_begin).Z0; + float phi0_wrt0 = tracks.at(rp.tracks_begin).phi; + + TVector3 X( - D0_wrt0 * TMath::Sin(phi0_wrt0) , D0_wrt0 * TMath::Cos(phi0_wrt0) , Z0_wrt0); + TVector3 x = X - V; + + TVector3 p(rp.momentum.x, rp.momentum.y, rp.momentum.z); + + double a = - rp.charge * Bz * cSpeed; + double pt = p.Pt(); + double r2 = x(0) * x(0) + x(1) * x(1); + double cross = x(0) * p(1) - x(1) * p(0); + double D=-9; + if (pt * pt - 2 * a * cross + a * a * r2 > 0) { + double T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2); + if (pt < 10.0) D = (T - pt) / a; + else D = (-2 * cross + a * r2) / (T + pt); + } + out.push_back(D); + + } else { + out.push_back(-9.); + } + } + return out; + } + + + + ROOT::VecOps::RVec XPtoPar_dz(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + const double cSpeed = 2.99792458e8 * 1.0e-9; //Reduced speed of light ??? + + ROOT::VecOps::RVec out; + + for (const auto & rp: in) { + + if( rp.tracks_begin < tracks.size()) { + + float D0_wrt0 = tracks.at(rp.tracks_begin).D0; + float Z0_wrt0 = tracks.at(rp.tracks_begin).Z0; + float phi0_wrt0 = tracks.at(rp.tracks_begin).phi; + + TVector3 X( - D0_wrt0 * TMath::Sin(phi0_wrt0) , D0_wrt0 * TMath::Cos(phi0_wrt0) , Z0_wrt0); + TVector3 x = X - V; + + TVector3 p(rp.momentum.x, rp.momentum.y, rp.momentum.z); + + double a = - rp.charge * Bz * cSpeed; + double pt = p.Pt(); + double C = a/(2 * pt); + double r2 = x(0) * x(0) + x(1) * x(1); + double cross = x(0) * p(1) - x(1) * p(0); + double T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2); + double D; + if (pt < 10.0) D = (T - pt) / a; + else D = (-2 * cross + a * r2) / (T + pt); + double B = C * TMath::Sqrt(TMath::Max(r2 - D * D, 0.0) / (1 + 2 * C * D)); + double st = TMath::ASin(B) / C; + double ct = p(2) / pt; + double z0; + double dot = x(0) * p(0) + x(1) * p(1); + if (dot > 0.0) z0 = x(2) - ct * st; + else z0 = x(2) + ct * st; + + out.push_back(z0); + } else { + out.push_back(-9.); + } + } + return out; + } + + ROOT::VecOps::RVec XPtoPar_phi(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + const double cSpeed = 2.99792458e8 * 1.0e-9; //Reduced speed of light ??? + + ROOT::VecOps::RVec out; + + for (const auto & rp: in) { + + if( rp.tracks_begin < tracks.size()) { + + float D0_wrt0 = tracks.at(rp.tracks_begin).D0; + float Z0_wrt0 = tracks.at(rp.tracks_begin).Z0; + float phi0_wrt0 = tracks.at(rp.tracks_begin).phi; + + TVector3 X( - D0_wrt0 * TMath::Sin(phi0_wrt0) , D0_wrt0 * TMath::Cos(phi0_wrt0) , Z0_wrt0); + TVector3 x = X - V; + + TVector3 p(rp.momentum.x, rp.momentum.y, rp.momentum.z); + + double a = - rp.charge * Bz * cSpeed; + double pt = p.Pt(); + double r2 = x(0) * x(0) + x(1) * x(1); + double cross = x(0) * p(1) - x(1) * p(0); + double T = TMath::Sqrt(pt * pt - 2 * a * cross + a * a * r2); + double phi0 = TMath::ATan2((p(1) - a * x(0)) / T, (p(0) + a * x(1)) / T); + + out.push_back(phi0); + + } else { + out.push_back(-9.); + } + } + return out; + } + + ROOT::VecOps::RVec XPtoPar_C(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + const double cSpeed = 2.99792458e8 * 1.0e3 * 1.0e-15; + ROOT::VecOps::RVec out; + + for (const auto & rp: in) { + + if( rp.tracks_begin < tracks.size()) { + + TVector3 p(rp.momentum.x, rp.momentum.y, rp.momentum.z); + + double a = std::copysign(1.0, rp.charge) * Bz * cSpeed; + double pt = p.Pt(); + double C = a/(2 * pt); + + out.push_back(C); + } else { + out.push_back(-9.); + } + } + return out; + } + + ROOT::VecOps::RVec XPtoPar_ct(const ROOT::VecOps::RVec& in, + const ROOT::VecOps::RVec& tracks, + const TVector3& V, + const float& Bz) { + + const double cSpeed = 2.99792458e8 * 1.0e-9; + ROOT::VecOps::RVec out; + + for (const auto & rp: in) { + + if( rp.tracks_begin < tracks.size()) { + + TVector3 p(rp.momentum.x, rp.momentum.y, rp.momentum.z); + double pt = p.Pt(); + + double ct = p(2) / pt; + + out.push_back(ct); + + } else { + out.push_back(-9.); + } + } + return out; + } + + + + + + + ROOT::VecOps::RVec getRP2TRK_D0(ROOT::VecOps::RVec in, ROOT::VecOps::RVec tracks) { @@ -11,7 +229,7 @@ getRP2TRK_D0(ROOT::VecOps::RVec in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin i for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin for (auto & p: in) { if (p.tracks_begin in, for (auto & p: in) { if (p.tracks_begin in for (auto & p: in) { if (p.tracks_begin`. +The line `.Define("RP_px", "ReconstructedParticle::get_px(ReconstructedParticles)")` takes this branch and for all entries computes px of each particle; the output of this call is a branch called _RP_px_ containing an `RVec` per each event. + +The jet clustering is performed using the 4-momenta of the reconstructed particles. This operation returns two outputs: + - `jets_ee_genkt` : `RVec< fastjet::Pseudojet >` ; `Pseudojet` methods and attributes allow to access the overall jet properties (implemented in `JetClusteringUtils.cc`). + - `jetconstituents_ee_genkt` : `RVec< RVec < int > >` , which is a vector of vectors of integer labels which were assigned to the particles during the clustering by the function `JetClusteringUtils::set_pseudoJets` : +``` +std::vector JetClusteringUtils::set_pseudoJets(ROOT::VecOps::RVec px, + ROOT::VecOps::RVec py, + ROOT::VecOps::RVec pz, + ROOT::VecOps::RVec e) { + std::vector result; + unsigned index = 0; + for (size_t i = 0; i < px.size(); ++i) { + result.emplace_back(px.at(i), py.at(i), pz.at(i), e.at(i)); + result.back().set_user_index(index); + ++index; + } + return result; +} +``` +and that are now used to associate the particles to the belonging jet. In fact, by calling +``` +.Define("JetsConstituents", "JetConstituentsUtils::build_constituents_cluster(ReconstructedParticles, jetconstituents_ee_genkt)") #build jet constituents +``` +we create an RVec::< RVec < ReconstructedParticleData > > in which the first index _i_ runs over the jets identified in the event and the second _j_ over the particles belonging to the _i_-th jet : +``` +rv::RVec build_constituents_cluster(const rv::RVec& rps, + const std::vector>& indices) { + //std::cout << "... Building jets-constituents after clustering " << std::endl; + rv::RVec jcs; + for (const auto& jet_index : indices) { + FCCAnalysesJetConstituents jc; + for(const auto& const_index : jet_index) { + jc.push_back(rps.at(const_index)); + } + jcs.push_back(jc); + } + return jcs; + } +``` +IMPORTANT: this function associates properly jet-constituents because the labels are set to the particles following the order of the particles as they present in the initial entry; so the index coincides with the particle position in the input branch. + +The goodness of the association particles-jets was proven by plotting the histograms of the residuals $P^\mu_{jet}-\sum{P^{\mu}_{constituents}})/P^\mu_{jet}$ with $P^\mu = (E, \vec{p})$. Actually the residuals were computed on other 4 quantities which are in one to one relation with 4-momenta: $E, p_t, \phi, \theta$. + + PLOTS!!! + +The structure RVec::< RVec < type > > is the key data structure for treating jets constituents in an event and will be mantained also when computing the features of the constituents for this stage. +Let's see an example of a function computing a feature of the constituents of the jets in one event (the same structure is mantained for other functions): +``` +rv::RVec get_erel_log_cluster(const rv::RVec& jets, + const rv::RVec& jcs) { + rv::RVec out; + for (size_t i = 0; i < jets.size(); ++i) { // external loop: first index (i) running over the jets in the event + auto& jet_csts = out.emplace_back(); //in the ext. loop we fill the external vector with vectors + //(one per jet) + float e_jet = jets.at(i).E(); + auto csts = get_jet_constituents(jcs, i); + for (const auto& jc : csts) { //second index (j) running over the constituents + //of the i-th jet + float val = (e_jet > 0.) ? jc.energy / e_jet : 1.; + float erel_log = float(std::log10(val)); + jet_csts.emplace_back(erel_log); //in the int. loop we fill the internal vectors (jets) + //with constituents data + } + } + return out; + } +``` +Notice: the functions are written considering one event, then the processing of all the events in the tree is performed by `fccanalysis run` command (see `produceTrainingTrees_mp.py`). + +The computed features were also validated by comparing them with the ones obtained in a previous study by Michele Selvaggi and Loukas Gouskos (samples spring2021). We obtained the following plots. + + PLOTS!!! + +At the end of this stage we have a tree in which each entry is an event; the features of the jets are saved as `RVec < float > ` and the features of the constituents of the jets are saved as `RVec < RVec < float > >` ; the return type is actually a pointer to these structures. We still need to rearrange the structure from a per-event tree of pointers to RVec to a per-jet tree of arrays (an ntuple). + +#### What if implicit clustering ? + ... + +### Stage_ntuple : `stage2.cpp` +The main goal of this stage is to rearrange the tree obtained in _stage1_ to a per-jet format, but other tasks are accomplished: +* setting the flags of the class which the jets belong to; +* checking the number of events actually considered is the wanted one; +* there is a $\sim 30\%$ cases in which the clustering returns more than 2 jets, and $\sim$ few per million cases in which less than 2 jets are returned; so in the first case just the two higher energy jets are considered, while in the second case no jet is considered; a count of this events is printed to stdout. + +`stage2.cpp` takes 4 arguments: `USAGE: ./stage [root_inFileName] [root_outFileName] N_i N_f` +1. [root_inFileName] : path to input file in the form `path_to_stage1file/stage1_infilename` , +2. [root_outFileName] : path to output file in the form `path_to_outputdir/outfilename` , +3. N_i : index of the event from where start reading the tree , +4. N_f : index of the event ehrtr to stop reading the tree . + +Our choices are implemented in the app `all_stages.py`, so will be explained in the next section. +Now, let's go through the code. + +###### Loop structure + +The general structure of the loop is: +``` +loop : events { + ... + loop : jets { + ... + loop : constituents { + ... + } + ntuple.Fill() + } +} +``` + +The position of ntuple.Fill() inside the loop structure determines the per-jet structure. +We modified this basic structure: +* we insert where to start looping by `N_i` and how many events to consider by `Nevents_Max = N_f - N_i`; +* we loop over the events from N_i to the end of the tree, but introduce an external counter `saved_events_counts` which grows only if the event has been saved; this is the reliable counter! In fact, if the event is "strange" we don't save any jet, we skip to the next event. The loop breaks when `saved_events_counts == Nevents_Max` if the loop has not reached the end of the input file. Let's study different cases: + - If `nentries - N_i < Nevents_Max` $\implies$ in the file there are less events than required; no error produced, but in stdout will be printed the actual number of saved events; + - If `nentries - N_i > Nevents_Max` $\implies$ in the file there are exactly $Nevents_{Max}$; + - If `nentries - N_i = Nevents_Max` but there are strange events I will have less saved events, no error, but I know from stdout. + +Are considered "strange" events those in which the clustering returns `njets < 2`; in that case we skip the loop. If `njets >= 2` the loop is performed on the first two jets only (the ones having more ENERGY, they're ordered in stage1; the successives not expected: leak in clustering). + +``` +int N_i = atoi(argv[3]); //where to start reading the tree +int N_f = atoi(argv[4]); +int Nevents_Max = N_f - N_i; //number of events to be saved +int saved_events_counts = 0; //counts the number of events actually saved +int nentries = ev->GetEntries(); //total number of events in the tree + +for(int i = N_i+1; i < nentries; ++i) { // Loop over the events + ev->GetEntry(i); + njet = nJets; + + if (njet < 2) { //exclude the events with less than two jets + continue ; + } + + for(int j=0; j < 2; ++j) { //Loop over the jets inside the i-th event + //we only take the first two jets (the ones having more ENERGY, they're ordered in stage1), + //the third not expected: leak in clustering + + ... + + nconst = (count_Const->at(j)); + + for(int k = 0; k < nconst; ++k){ //Loop over the constituents of the j-th jet in the i-th event + ... + } + ntuple->Fill(); + } + saved_events_counts += 1; //we count the num of events saved + if (saved_events_counts == Nevents_Max) { //interrupt the loop if Nevents_max events have already been saved + break; + } +} +``` + + +###### Translating +As an example, we take one jet feature (stored as `RVec < float> *` in the per-event tree) and one constituent feature (`RVec < RVec < float> > *`) and follow them through the code. + +``` +//Setting variables for reading +int nJets; //number of jets in the event +int nconst = 0; //number of constituents of the jets +ROOT::VecOps::RVec *Jets_e=0; +ROOT::VecOps::RVec > *JetsConstituents_e = 0; + +ev->SetBranchAddress("njet", &nJets); +ev->SetBranchAddress("nconst", &count_Const); +ev->SetBranchAddress("Jets_e", &Jets_e); +ev->SetBranchAddress("JetsConstituents_e", &JetsConstituents_e); + +int njet = 0; +int nconst = 0; +double recojet_e; +float pfcand_e[1000] = {0.}; //here we initialize wit a large size + +ntuple->Branch("nconst", &nconst, "nconst/I"); +ntuple->Branch("recojet_e", &recojet_e); +ntuple->Branch("pfcand_e", pfcand_e, "pfcand_e[nconst]/F"); + +loop : i over events + loop : j over jets + + recojet_e = (*Jets_e)[j]; //Pointer usage + nconst = (count_Const->at(j)); + + loop : k over constituents + pfcand_e[k] = (JetsConstituents_e->at(j))[k]; //k-th element of the j-th vector + //pointed by JetsConstituents_e +``` + +In the ntuple, a jet overall property will be just a `double`. +Since we don't know a priori how many constituents the jet has, we initialize a constituent feature as an array with a number of elements larger than possible number of constituents, to be sure that we succeed in reading all of them from the input file (this array works as a temporary "home" before being sent to the output file); we initialize all the values to 0. When we create the branch of the ntuple (`ntuple->Branch("pfcand_e", pfcand_e, "pfcand_e[nconst]/F")`), we pass as the size of the entry the actual number of constituents `nconst` which has already been read into another branch (during the loop), in order not to save all the padding zeros contained in the local array. + +###### Setting the flags +We read the name and take the character in the name corresponding to the class (in this case last letter before .root); we fix the flag in the beginning and never change it anymore; since it is pointing to the ntuple branch, everytime I call `ntuple.Fill()` the same value will be added to the branch. +``` +std::string infileName(argv[1]); +char flavour = infileName[infileName.length()-6]; + +float is_q = 0.; +float is_b = 0.; +float is_c = 0.; +float is_s = 0.; +float is_g = 0.; + +if (flavour == 'q') {is_q = 1.;} +if (flavour == 'b') {is_b = 1.;} +if (flavour == 'c') {is_c = 1.;} +if (flavour == 's') {is_s = 1.;} +if (flavour == 'g') {is_g = 1.;} +if (flavour == 't') {is_t = 1.;} + +ntuple->Branch("pfcand_isMu", pfcand_isMu, "pfcand_isMu[nconst]/F"); +ntuple->Branch("pfcand_isEl", pfcand_isEl, "pfcand_isEl[nconst]/F"); +ntuple->Branch("pfcand_isChargedHad", pfcand_isChargedHad, "pfcand_isChargedHad[nconst]/F"); +ntuple->Branch("pfcand_isGamma", pfcand_isGamma, "pfcand_isGamma[nconst]/F"); +ntuple->Branch("pfcand_isNeutralHad", pfcand_isNeutralHad, "pfcand_isNeutralHad[nconst]/F"); + +loop : events { + ... + loop: jets { + ... + loop: constituents { + ... + } + ntuple.Fill() + } +} +``` + diff --git a/examples/FCCee/weaver/all_stages.py b/examples/FCCee/weaver/all_stages.py new file mode 100644 index 0000000000..cbcfd4b679 --- /dev/null +++ b/examples/FCCee/weaver/all_stages.py @@ -0,0 +1,145 @@ +import sys +import os +import argparse +import multiprocessing as mp +import subprocess +from subprocess import Popen, PIPE +from datetime import date +import time + + +# ________________________________________________________________________________ +def main(): + parser = argparse.ArgumentParser() + + parser.add_argument( + "--indir", + help="path input directory", + default="/eos/experiment/fcc/ee/generation/DelphesEvents/pre_fall2022_training/IDEA/", + ) + parser.add_argument( + "--outdir", + help="path output directory", + default="/eos/experiment/fcc/ee/jet_flavour_tagging/pre_fall2022_training/IDEA/", + ) + parser.add_argument("--nev", help="number of events", type=int, default=1000000) + parser.add_argument("--fsplit", help="fraction train/test", type=float, default=0.9) + parser.add_argument("--dry", help="dry run ", action='store_true') + + args = parser.parse_args() + inDIR = args.indir + outDIR = args.outdir + + # create necessary subdirectories + actualDIR = subprocess.check_output(["bash", "-c", "pwd"], universal_newlines=True)[:-1] + username = subprocess.check_output(["bash", "-c", "echo $USER"], universal_newlines=True)[:-1] + today = date.today().strftime("%Y%b%d") + subdir = username + "_" + today + subprocess.call(["bash", "-c", "mkdir -p {}".format(subdir)], cwd=outDIR) + OUTDIR = outDIR + username + "_" + today + "/" + print(OUTDIR) + + # set total number of events + N = args.nev + frac_split = args.fsplit + N_split = int(frac_split * N) + + print("") + print("=================================================") + print("") + print(" INDIR = {}".format(inDIR)) + print(" OUTDIR = {}".format(OUTDIR)) + print(" NEVENTS = {}".format(N)) + print(" FRAC SPLIT = {}".format(frac_split)) + print("") + + # compile stage2 c++ code + cmd_compile = "g++ -o stage2 stage2.cpp `root-config --cflags --libs` -Wall" + subprocess.check_call(cmd_compile, shell=True, stdout=None, stderr=None) + + # create commands + #cmdbase_stage1 = "fccanalysis run stage1.py --nevents {}".format(N) + cmdbase_stage1 = "fccanalysis run stage1.py".format(N) + opt1_out = " --output {}stage1_ee_ZH_vvCLASS.root ".format(OUTDIR) + opt1_in = " --files-list {}p8_ee_ZH_Znunu_HCLASS_ecm240/*.root ".format(inDIR) + wait = " ; sleep 30;" + cmd_stage1 = cmdbase_stage1 + opt1_out + opt1_in + + cmdbase_stagentuple = "./stage2 DIRstage1_ee_ZH_vvCLASS.root DIRntuple_MOD_ee_ZH_vvCLASS.root " + cmd_stage2 = cmdbase_stagentuple.replace("DIR", OUTDIR) + + samples = ["bb", "cc", "ss", "gg", "qq"] + mods = ["train", "test"] + + # create files storing stdout and stderr + list_stdout = [open(OUTDIR + "{}_stdout.txt".format(sample), "w") for sample in samples] + list_stderr = [open(OUTDIR + "{}_stderr.txt".format(sample), "w") for sample in samples] + + ###=== RUN STAGE 1 + for i, sample in enumerate(samples): + + if sample == "qq": + cmd_stage1_f = ( + cmdbase_stage1 + + opt1_out.replace("CLASS", "qq") + + " --files-list {}p8_ee_ZH_Znunu_Huu_ecm240/*.root {}p8_ee_ZH_Znunu_Hdd_ecm240/*.root".format( + inDIR, inDIR + ) + ) + else: + cmd_stage1_f = cmd_stage1.replace("CLASS", sample) + + print(cmd_stage1_f) + # run stage1 + start1_time = time.time() + if not args.dry: + subprocess.check_call( + cmd_stage1_f, shell=True, stdout=list_stdout[i], stderr=list_stderr[i] + ) + end1_time = time.time() + list_stdout[i].write("Stage1 time: {} \n".format(end1_time - start1_time)) + + ###=== RUN STAGE 2 + threads = [] + for i, sample in enumerate(samples): + + cmd_stage2_train = cmd_stage2.replace("CLASS", sample).replace( + "MOD", mods[0] + ) + " {} {} ".format(0, N_split) + cmd_stage2_test = cmd_stage2.replace("CLASS", sample).replace( + "MOD", mods[1] + ) + " {} {} ".format(N_split, N) + print(cmd_stage2_train) + print(cmd_stage2_test) + + if not args.dry: + thread = mp.Process( + target=ntuplizer, + args=(cmd_stage2_train, cmd_stage2_test, list_stdout[i], list_stderr[i]), + ) + thread.start() + threads.append(thread) + + for proc in threads: + proc.join() + + + for i in range(len(list_stdout)): + list_stdout[i].close() + list_stderr[i].close() + + +# ________________________________________________________________________________ +def ntuplizer(cmd_stage2_train, cmd_stage2_test, f_stdout, f_stderr): + + start2_time = time.time() + subprocess.check_call(cmd_stage2_train, shell=True, stdout=f_stdout, stderr=f_stderr) + subprocess.check_call(cmd_stage2_test, shell=True, stdout=f_stdout, stderr=f_stderr) + end2_time = time.time() + f_stdout.write("stage2 time (run only): {} \n".format(end2_time - start2_time)) + + +# _______________________________________________________________________________________ +if __name__ == "__main__": + main() + diff --git a/examples/FCCee/weaver/analysis_inference.py b/examples/FCCee/weaver/analysis_inference.py new file mode 100644 index 0000000000..b449ffad7a --- /dev/null +++ b/examples/FCCee/weaver/analysis_inference.py @@ -0,0 +1,263 @@ +#Mandatory: List of processes +processList = { + #prefall2022 samples (generated centrally) + 'p8_ee_ZH_Znunu_Hbb_ecm240':{}, #1030000 events + 'p8_ee_ZH_Znunu_Hcc_ecm240':{}, #1060000 + 'p8_ee_ZH_Znunu_Hss_ecm240':{}, #1060000 + 'p8_ee_ZH_Znunu_Hgg_ecm240':{'fraction':0.5}, #2000000 + 'p8_ee_ZH_Znunu_Huu_ecm240':{'fraction':0.5}, #we take only half sample for uu,dd because they will go into qq label which contains both + 'p8_ee_ZH_Znunu_Hdd_ecm240':{'fraction':0.5}, #and we want for qq same number of jets as other classes; the two files 2080000 events in total, 1040000 each? +} + +#Mandatory: Production tag when running over EDM4Hep centrally produced events, this points to the yaml files for getting sample statistics +#prodTag = "FCCee/spring2021/IDEA/" +prodTag = "FCCee/pre_fall2022_training/IDEA/" #for prefall2022 samples + +#Optional: output directory, default is local running directory +#outputDir = "/eos/home-a/adelvecc/FCCevaluate/" + +#Optional +nCPUS = 8 +runBatch = False +#batchQueue = "longlunch" +#compGroup = "group_u_FCC.local_gen" + + +#USER DEFINED CODE + +#Mandatory: RDFanalysis class where the use defines the operations on the TTree +class RDFanalysis(): + #__________________________________________________________ + #Mandatory: analysers funtion to define the analysers to process, please make sure you return the last dataframe, in this example it is df2 + def analysers(df): + + from ROOT import JetFlavourUtils + weaver = JetFlavourUtils.setup_weaver( + "/eos/experiment/fcc/ee/jet_flavour_tagging/pre_fall2022_training/IDEA/selvaggi_2022Oct30/fccee_flavtagging_edm4hep_v2.onnx", #name of the trained model exported + "/eos/experiment/fcc/ee/jet_flavour_tagging/pre_fall2022_training/IDEA/selvaggi_2022Oct30/preprocess_fccee_flavtagging_edm4hep_v2.json", #.json file produced by weaver during training + ( + "pfcand_erel_log", #list of the training variables, + "pfcand_thetarel", #will be used for predictions as well + "pfcand_phirel", + "pfcand_dxy", + "pfcand_dz", + "pfcand_dptdpt", + "pfcand_detadeta", + "pfcand_dphidphi", + "pfcand_dxydxy", + "pfcand_dzdz", + "pfcand_dxydz", + "pfcand_dphidxy", + "pfcand_dlambdadz", + "pfcand_dxyc", + "pfcand_dxyctgtheta", + "pfcand_phic", + "pfcand_phidz", + "pfcand_phictgtheta", + "pfcand_cdz", + "pfcand_cctgtheta", + "pfcand_mtof", + "pfcand_dndx", + "pfcand_charge", + "pfcand_isMu", + "pfcand_isEl", + "pfcand_isChargedHad", + "pfcand_isGamma", + "pfcand_isNeutralHad", + "pfcand_btagSip2dVal", + "pfcand_btagSip2dSig", + "pfcand_btagSip3dVal", + "pfcand_btagSip3dSig", + "pfcand_btagJetDistVal", + "pfcand_btagJetDistSig", + ), + ) + + df2 = ( + df + ### COMPUTE THE VARIABLES FOR INFERENCE OF THE TRAINING MODEL + ### This section should be equal to the one used for training + + ### MC primary vertex + .Define("MC_PrimaryVertex", "FCCAnalyses::MCParticle::get_EventPrimaryVertex(21)( Particle )" ) + + # CLUSTERING + #define the RP px, py, pz and e + .Define("RP_px", "ReconstructedParticle::get_px(ReconstructedParticles)") + .Define("RP_py", "ReconstructedParticle::get_py(ReconstructedParticles)") + .Define("RP_pz", "ReconstructedParticle::get_pz(ReconstructedParticles)") + .Define("RP_e", "ReconstructedParticle::get_e(ReconstructedParticles)") + .Define("RP_m", "ReconstructedParticle::get_mass(ReconstructedParticles)") + .Define("RP_q", "ReconstructedParticle::get_charge(ReconstructedParticles)") + + #build pseudo jets with the RP + .Define("pseudo_jets", "JetClusteringUtils::set_pseudoJets(RP_px, RP_py, RP_pz, RP_e)") + #run jet clustering with all reconstructed particles. ee_genkt_algorithm, R=1.5, inclusive clustering, E-scheme + .Define("FCCAnalysesJets_ee_genkt", "JetClustering::clustering_ee_genkt(1.5, 0, 0, 0, 0, -1)(pseudo_jets)") + #get the jets out of the struct + .Define("jets_ee_genkt", "JetClusteringUtils::get_pseudoJets(FCCAnalysesJets_ee_genkt)") + #get the jets constituents out of the struct + .Define("jetconstituents_ee_genkt","JetClusteringUtils::get_constituents(FCCAnalysesJets_ee_genkt)") + + #===== COMPUTE TRAINING FEATURES + + .Define("JetsConstituents", "JetConstituentsUtils::build_constituents_cluster(ReconstructedParticles, jetconstituents_ee_genkt)") #build jet constituents lists + + ### Types of particles + .Alias("MCRecoAssociations0", "MCRecoAssociations#0.index") + .Alias("MCRecoAssociations1", "MCRecoAssociations#1.index") + .Define("JetsConstituents_Pids", "JetConstituentsUtils::get_PIDs_cluster(MCRecoAssociations0, MCRecoAssociations1, ReconstructedParticles, Particle, jetconstituents_ee_genkt)") + + .Define("pfcand_isMu", "JetConstituentsUtils::get_isMu(JetsConstituents_Pids)") + .Define("pfcand_isEl", "JetConstituentsUtils::get_isEl(JetsConstituents_Pids)") + .Define("pfcand_isChargedHad", "JetConstituentsUtils::get_isChargedHad(JetsConstituents_Pids, JetsConstituents)") + .Define("pfcand_isGamma", "JetConstituentsUtils::get_isGamma(JetsConstituents_Pids)") + .Define("pfcand_isNeutralHad", "JetConstituentsUtils::get_isNeutralHad(JetsConstituents_Pids, JetsConstituents)") + + ### Kinematics, displacement, PID + + .Define("pfcand_erel", "JetConstituentsUtils::get_erel_cluster(jets_ee_genkt, JetsConstituents)") + .Define("pfcand_erel_log", "JetConstituentsUtils::get_erel_log_cluster(jets_ee_genkt, JetsConstituents)") + .Define("pfcand_thetarel", "JetConstituentsUtils::get_thetarel_cluster(jets_ee_genkt, JetsConstituents)") + .Define("pfcand_phirel", "JetConstituentsUtils::get_phirel_cluster(jets_ee_genkt, JetsConstituents)") + + .Define("pfcand_charge", "JetConstituentsUtils::get_charge(JetsConstituents)") + .Define("pfcand_dndx", "JetConstituentsUtils::get_dndx(JetsConstituents, EFlowTrack_2, EFlowTrack, pfcand_isChargedHad)") + .Define("pfcand_mtof", "JetConstituentsUtils::get_mtof(JetsConstituents, EFlowTrack_L, EFlowTrack, TrackerHits, JetsConstituents_Pids)") + + .Define("Bz", "ReconstructedParticle2Track::Bz(ReconstructedParticles, EFlowTrack_1)") + + .Define("pfcand_dxy", "JetConstituentsUtils::XPtoPar_dxy(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("pfcand_dz", "JetConstituentsUtils::XPtoPar_dz(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("pfcand_phi0", "JetConstituentsUtils::XPtoPar_phi(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("pfcand_C", "JetConstituentsUtils::XPtoPar_C(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("pfcand_ct", "JetConstituentsUtils::XPtoPar_ct(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + + .Define("pfcand_dptdpt", "JetConstituentsUtils::get_omega_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dxydxy", "JetConstituentsUtils::get_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dzdz", "JetConstituentsUtils::get_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dphidphi", "JetConstituentsUtils::get_phi0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_detadeta", "JetConstituentsUtils::get_tanlambda_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dxydz", "JetConstituentsUtils::get_d0_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dphidxy", "JetConstituentsUtils::get_phi0_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_phidz", "JetConstituentsUtils::get_phi0_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_phictgtheta", "JetConstituentsUtils::get_tanlambda_phi0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dxyctgtheta", "JetConstituentsUtils::get_tanlambda_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dlambdadz", "JetConstituentsUtils::get_tanlambda_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_cctgtheta", "JetConstituentsUtils::get_omega_tanlambda_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_phic", "JetConstituentsUtils::get_omega_phi0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_dxyc", "JetConstituentsUtils::get_omega_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("pfcand_cdz", "JetConstituentsUtils::get_omega_z0_cov(JetsConstituents, EFlowTrack_1)") + + .Define("pfcand_btagSip2dVal", "JetConstituentsUtils::get_Sip2dVal_clusterV(jets_ee_genkt, pfcand_dxy, pfcand_phi0, MC_PrimaryVertex, Bz)") + .Define("pfcand_btagSip2dSig", "JetConstituentsUtils::get_Sip2dSig(pfcand_btagSip2dVal, pfcand_dxydxy)") + .Define("pfcand_btagSip3dVal", "JetConstituentsUtils::get_Sip3dVal_clusterV(jets_ee_genkt, pfcand_dxy, pfcand_dz, pfcand_phi0, MC_PrimaryVertex, Bz)") + .Define("pfcand_btagSip3dSig", "JetConstituentsUtils::get_Sip3dSig(pfcand_btagSip3dVal, pfcand_dxydxy, pfcand_dzdz)") + .Define("pfcand_btagJetDistVal", "JetConstituentsUtils::get_JetDistVal_clusterV(jets_ee_genkt, JetsConstituents, pfcand_dxy, pfcand_dz, pfcand_phi0, MC_PrimaryVertex, Bz)") + .Define("pfcand_btagJetDistSig", "JetConstituentsUtils::get_JetDistSig(pfcand_btagJetDistVal, pfcand_dxydxy, pfcand_dzdz)") + + + ##### RUN INFERENCE (fixed by the previous section) + + .Define( + "MVAVec", + "JetFlavourUtils::get_weights(\ + pfcand_erel_log,\ + pfcand_thetarel,\ + pfcand_phirel,\ + pfcand_dxy,\ + pfcand_dz,\ + pfcand_dptdpt,\ + pfcand_dphidphi,\ + pfcand_detadeta,\ + pfcand_dxydxy,\ + pfcand_dzdz,\ + pfcand_dxydz,\ + pfcand_dphidxy,\ + pfcand_dlambdadz,\ + pfcand_dxyc,\ + pfcand_dxyctgtheta,\ + pfcand_phic,\ + pfcand_phidz,\ + pfcand_phictgtheta,\ + pfcand_cdz,\ + pfcand_cctgtheta,\ + pfcand_mtof,\ + pfcand_dndx,\ + pfcand_charge,\ + pfcand_isMu,\ + pfcand_isEl,\ + pfcand_isChargedHad,\ + pfcand_isGamma,\ + pfcand_isNeutralHad,\ + pfcand_btagSip2dVal,\ + pfcand_btagSip2dSig,\ + pfcand_btagSip3dVal,\ + pfcand_btagSip3dSig,\ + pfcand_btagJetDistVal,\ + pfcand_btagJetDistSig\ + )", + ) + + ##### RECAST OUTPUT (get predictions per each sample) + .Define("recojet_isG", "JetFlavourUtils::get_weight(MVAVec, 0)") + .Define("recojet_isQ", "JetFlavourUtils::get_weight(MVAVec, 1)") + .Define("recojet_isS", "JetFlavourUtils::get_weight(MVAVec, 2)") + .Define("recojet_isC", "JetFlavourUtils::get_weight(MVAVec, 3)") + .Define("recojet_isB", "JetFlavourUtils::get_weight(MVAVec, 4)") + + ##### COMPUTE OBSERVABLES FOR ANALYSIS + #if not changing training etc... but only interested in the analysis using a trained model (fixed classes), you should only operate in this section. + #if you're interested in saving variables used for training don't need to compute them again, just + #add them to the list in at the end of the code + + #EXAMPLE + + #EVENT LEVEL + .Define("njet", "JetConstituentsUtils::count_jets(JetsConstituents)") + + #JET LEVEL + #jet kinematics + .Define("recojet_pt", "JetClusteringUtils::get_pt(jets_ee_genkt)") + .Define("recojet_e", "JetClusteringUtils::get_e(jets_ee_genkt)") + .Define("recojet_mass", "JetClusteringUtils::get_m(jets_ee_genkt)") + .Define("recojet_phi", "JetClusteringUtils::get_phi(jets_ee_genkt)") + .Define("recojet_theta", "JetClusteringUtils::get_theta(jets_ee_genkt)") + + .Define("tlv_jets", "JetConstituentsUtils::compute_tlv_jets(jets_ee_genkt)") + .Define("invariant_mass", "JetConstituentsUtils::InvariantMass(tlv_jets[0], tlv_jets[1])") + + #counting types of particles composing the jet + .Define("nconst", "JetConstituentsUtils::count_consts(JetsConstituents)") + .Define("nmu", "JetConstituentsUtils::count_type(pfcand_isMu)") + .Define("nel", "JetConstituentsUtils::count_type(pfcand_isEl)") + .Define("nchargedhad", "JetConstituentsUtils::count_type(pfcand_isChargedHad)") + .Define("nphoton", "JetConstituentsUtils::count_type(pfcand_isGamma)") + .Define("nneutralhad", "JetConstituentsUtils::count_type(pfcand_isNeutralHad)") + + #CONSTITUENTS LEVEL + .Define("pfcand_e", "JetConstituentsUtils::get_e(JetsConstituents)") + .Define("pfcand_pt", "JetConstituentsUtils::get_pt(JetsConstituents)") + .Define("pfcand_theta", "JetConstituentsUtils::get_theta(JetsConstituents)") + .Define("pfcand_phi", "JetConstituentsUtils::get_phi(JetsConstituents)") + + ) + + return df2 + + #__________________________________________________________ + #SAVE PREDICTIONS & OBSERVABLES FOR ANALYSIS + #Mandatory: output function, please make sure you return the branchlist as a python list + def output(): + branchList = [ + #predictions + 'recojet_isG', 'recojet_isQ', 'recojet_isS', 'recojet_isC', 'recojet_isB', + #observables + 'recojet_mass', 'recojet_e', 'recojet_pt', + 'invariant_mass', + 'nconst', 'nchargedhad', + 'pfcand_e', 'pfcand_pt', 'pfcand_phi', + 'pfcand_erel', + 'pfcand_erel_log', + ] + return branchList diff --git a/examples/FCCee/weaver/stage1.py b/examples/FCCee/weaver/stage1.py new file mode 100644 index 0000000000..aeb35a483e --- /dev/null +++ b/examples/FCCee/weaver/stage1.py @@ -0,0 +1,193 @@ +processList = { + #prefall2022 samples (generated centrally) + 'p8_ee_ZH_Znunu_Hbb_ecm240':{}, #1030000 events + 'p8_ee_ZH_Znunu_Hcc_ecm240':{}, #1060000 + 'p8_ee_ZH_Znunu_Hss_ecm240':{}, #1060000 + 'p8_ee_ZH_Znunu_Hgg_ecm240':{'fraction':0.5}, #2000000 + 'p8_ee_ZH_Znunu_Huu_ecm240':{'fraction':0.5}, #we take only half sample for uu,dd because they will go into qq label which contains both + 'p8_ee_ZH_Znunu_Hdd_ecm240':{'fraction':0.5}, #and we want for qq same number of jets as other classes; the two files 2080000 events in total, 1040000 each? +} + + +#prodTag = "FCCee/spring2021/IDEA/" #for spring2022 samples +prodTag = "FCCee/pre_fall2022_training/IDEA/" #for prefall2022 samples + +#outputDir = "" + +#Optional: ncpus, default is 4 +nCPUS = 8 + + +class RDFanalysis(): + + def analysers(df): + + df2 = ( + df + + #===== VERTEX + + #MC primary vertex + .Define("MC_PrimaryVertex", "FCCAnalyses::MCParticle::get_EventPrimaryVertex(21)( Particle )" ) + + #===== CLUSTERING + #define the RP px, py, pz and e + .Define("RP_px", "ReconstructedParticle::get_px(ReconstructedParticles)") + .Define("RP_py", "ReconstructedParticle::get_py(ReconstructedParticles)") + .Define("RP_pz", "ReconstructedParticle::get_pz(ReconstructedParticles)") + .Define("RP_e", "ReconstructedParticle::get_e(ReconstructedParticles)") + .Define("RP_m", "ReconstructedParticle::get_mass(ReconstructedParticles)") + .Define("RP_q", "ReconstructedParticle::get_charge(ReconstructedParticles)") + + #build pseudo jets with the RP, using the interface that takes px,py,pz,E + #.Define("pseudo_jets", "JetClusteringUtils::set_pseudoJets_xyzm(RP_px, RP_py, RP_pz, RP_m)") + .Define("pseudo_jets", "JetClusteringUtils::set_pseudoJets(RP_px, RP_py, RP_pz, RP_e)") + #run jet clustering with all reconstructed particles. ee_genkt_algorithm, R=1.5, inclusive clustering, E-scheme + .Define("FCCAnalysesJets_ee_genkt", "JetClustering::clustering_ee_genkt(1.5, 0, 0, 0, 0, -1)(pseudo_jets)") + #get the jets out of the struct + .Define("jets_ee_genkt", "JetClusteringUtils::get_pseudoJets(FCCAnalysesJets_ee_genkt)") + #get the jets constituents out of the struct + .Define("jetconstituents_ee_genkt","JetClusteringUtils::get_constituents(FCCAnalysesJets_ee_genkt)") + + + #===== OBSERVABLES + #JET LEVEL + ###.Define("Jets_px", "JetClusteringUtils::get_px(jets_ee_genkt)") #jets_ee_genkt_px + ###.Define("Jets_py", "JetClusteringUtils::get_py(jets_ee_genkt)") + ###.Define("Jets_pz", "JetClusteringUtils::get_pz(jets_ee_genkt)") + + .Define("Jets_pt", "JetClusteringUtils::get_pt(jets_ee_genkt)") + .Define("Jets_e", "JetClusteringUtils::get_e(jets_ee_genkt)") + .Define("Jets_mass", "JetClusteringUtils::get_m(jets_ee_genkt)") + .Define("Jets_phi", "JetClusteringUtils::get_phi(jets_ee_genkt)") + .Define("Jets_theta", "JetClusteringUtils::get_theta(jets_ee_genkt)") + + #CONSTITUENT LEVEL + .Define("JetsConstituents", "JetConstituentsUtils::build_constituents_cluster(ReconstructedParticles, jetconstituents_ee_genkt)") #build jet constituents lists + + #getting the types of particles + .Alias("MCRecoAssociations0", "MCRecoAssociations#0.index") + .Alias("MCRecoAssociations1", "MCRecoAssociations#1.index") + .Define("JetsConstituents_Pids", "JetConstituentsUtils::get_PIDs_cluster(MCRecoAssociations0, MCRecoAssociations1, ReconstructedParticles, Particle, jetconstituents_ee_genkt)") + .Define("JetsConstituents_isMu", "JetConstituentsUtils::get_isMu(JetsConstituents_Pids)") + .Define("JetsConstituents_isEl", "JetConstituentsUtils::get_isEl(JetsConstituents_Pids)") + .Define("JetsConstituents_isChargedHad", "JetConstituentsUtils::get_isChargedHad(JetsConstituents_Pids, JetsConstituents)") + .Define("JetsConstituents_isGamma", "JetConstituentsUtils::get_isGamma(JetsConstituents_Pids)") + .Define("JetsConstituents_isNeutralHad", "JetConstituentsUtils::get_isNeutralHad(JetsConstituents_Pids, JetsConstituents)") + + #kinematics, displacement, PID + .Define("JetsConstituents_e", "JetConstituentsUtils::get_e(JetsConstituents)") + .Define("JetsConstituents_pt", "JetConstituentsUtils::get_pt(JetsConstituents)") + .Define("JetsConstituents_theta", "JetConstituentsUtils::get_theta(JetsConstituents)") + .Define("JetsConstituents_phi", "JetConstituentsUtils::get_phi(JetsConstituents)") + .Define("JetsConstituents_charge", "JetConstituentsUtils::get_charge(JetsConstituents)") + + .Define("JetsConstituents_erel", "JetConstituentsUtils::get_erel_cluster(jets_ee_genkt, JetsConstituents)") + .Define("JetsConstituents_erel_log", "JetConstituentsUtils::get_erel_log_cluster(jets_ee_genkt, JetsConstituents)") + .Define("JetsConstituents_thetarel", "JetConstituentsUtils::get_thetarel_cluster(jets_ee_genkt, JetsConstituents)") + .Define("JetsConstituents_phirel", "JetConstituentsUtils::get_phirel_cluster(jets_ee_genkt, JetsConstituents)") + + .Define("JetsConstituents_dndx", "JetConstituentsUtils::get_dndx(JetsConstituents, EFlowTrack_2, EFlowTrack, JetsConstituents_isChargedHad)") + .Define("JetsConstituents_mtof", "JetConstituentsUtils::get_mtof(JetsConstituents, EFlowTrack_L, EFlowTrack, TrackerHits, JetsConstituents_Pids)") + + .Define("JetsConstituents_d0_wrt0", "JetConstituentsUtils::get_d0(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_z0_wrt0", "JetConstituentsUtils::get_z0(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_phi0_wrt0", "JetConstituentsUtils::get_phi0(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_omega_wrt0", "JetConstituentsUtils::get_omega(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_tanlambda_wrt0", "JetConstituentsUtils::get_tanLambda(JetsConstituents, EFlowTrack_1)") + + .Define("JetsConstituents_Bz", "JetConstituentsUtils::get_Bz(JetsConstituents, EFlowTrack_1)") + .Define("Bz", "ReconstructedParticle2Track::Bz(ReconstructedParticles, EFlowTrack_1)") + + #.Define("JetsConstituents_Par", "ReconstructedParticle2Track::XPtoPar(ReconstructedParticles, EFlowTrack_1, MC_PrimaryVertex, Bz)") + #.Define("JetsConstituents_dxy", "JetConstituentsUtils::XPtoPar_dxy(JetsConstituents_Par)") + #.Define("JetsConstituents_dxy", "JetConstituentsUtils::XPtoPar_dxy(JetsConstituents_Par)") + #.Define("JetsConstituents_dz", "JetConstituentsUtils::XPtoPar_dz(JetsConstituents_Par)") + #.Define("JetsConstituents_phi0", "JetConstituentsUtils::XPtoPar_phi0(JetsConstituents_Par)") + #.Define("JetsConstituents_C", "JetConstituentsUtils::XPtoPar_C(JetsConstituents_Par)") + #.Define("JetsConstituents_ct", "JetConstituentsUtils::XPtoPar_ct(JetsConstituents_Par)") + + .Define("JetsConstituents_dxy", "JetConstituentsUtils::XPtoPar_dxy(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("JetsConstituents_dz", "JetConstituentsUtils::XPtoPar_dz(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("JetsConstituents_phi0", "JetConstituentsUtils::XPtoPar_phi(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("JetsConstituents_C", "JetConstituentsUtils::XPtoPar_C(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + .Define("JetsConstituents_ct", "JetConstituentsUtils::XPtoPar_ct(JetsConstituents, EFlowTrack_1, MC_PrimaryVertex, Bz)") + + .Define("JetsConstituents_omega_cov", "JetConstituentsUtils::get_omega_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_d0_cov", "JetConstituentsUtils::get_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_z0_cov", "JetConstituentsUtils::get_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_phi0_cov", "JetConstituentsUtils::get_phi0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_tanlambda_cov", "JetConstituentsUtils::get_tanlambda_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_d0_z0_cov", "JetConstituentsUtils::get_d0_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_phi0_d0_cov", "JetConstituentsUtils::get_phi0_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_phi0_z0_cov", "JetConstituentsUtils::get_phi0_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_tanlambda_phi0_cov", "JetConstituentsUtils::get_tanlambda_phi0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_tanlambda_d0_cov", "JetConstituentsUtils::get_tanlambda_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_tanlambda_z0_cov", "JetConstituentsUtils::get_tanlambda_z0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_omega_tanlambda_cov", "JetConstituentsUtils::get_omega_tanlambda_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_omega_phi0_cov", "JetConstituentsUtils::get_omega_phi0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_omega_d0_cov", "JetConstituentsUtils::get_omega_d0_cov(JetsConstituents, EFlowTrack_1)") + .Define("JetsConstituents_omega_z0_cov", "JetConstituentsUtils::get_omega_z0_cov(JetsConstituents, EFlowTrack_1)") + + .Define("JetsConstituents_Sip2dVal", "JetConstituentsUtils::get_Sip2dVal_clusterV(jets_ee_genkt, JetsConstituents_dxy, JetsConstituents_phi0, MC_PrimaryVertex, Bz)") + .Define("JetsConstituents_Sip2dSig", "JetConstituentsUtils::get_Sip2dSig(JetsConstituents_Sip2dVal, JetsConstituents_d0_cov)") + .Define("JetsConstituents_Sip3dVal", "JetConstituentsUtils::get_Sip3dVal_clusterV(jets_ee_genkt, JetsConstituents_dxy, JetsConstituents_dz, JetsConstituents_phi0, MC_PrimaryVertex, Bz)") + .Define("JetsConstituents_Sip3dSig", "JetConstituentsUtils::get_Sip3dSig(JetsConstituents_Sip3dVal, JetsConstituents_d0_cov, JetsConstituents_z0_cov)") + .Define("JetsConstituents_JetDistVal", "JetConstituentsUtils::get_JetDistVal_clusterV(jets_ee_genkt, JetsConstituents, JetsConstituents_dxy, JetsConstituents_dz, JetsConstituents_phi0, MC_PrimaryVertex, Bz)") + .Define("JetsConstituents_JetDistSig", "JetConstituentsUtils::get_JetDistSig(JetsConstituents_JetDistVal, JetsConstituents_d0_cov, JetsConstituents_z0_cov)") + + #counting the types of particles per jet + .Define("njet", "JetConstituentsUtils::count_jets(JetsConstituents)") + .Define("nconst", "JetConstituentsUtils::count_consts(JetsConstituents)") + .Define("nmu", "JetConstituentsUtils::count_type(JetsConstituents_isMu)") + .Define("nel", "JetConstituentsUtils::count_type(JetsConstituents_isEl)") + .Define("nchargedhad", "JetConstituentsUtils::count_type(JetsConstituents_isChargedHad)") + .Define("nphoton", "JetConstituentsUtils::count_type(JetsConstituents_isGamma)") + .Define("nneutralhad", "JetConstituentsUtils::count_type(JetsConstituents_isNeutralHad)") + + #compute the residues jet-constituents on significant kinematic variables as a check + .Define("tlv_jets", "JetConstituentsUtils::compute_tlv_jets(jets_ee_genkt)") + .Define("sum_tlv_jcs", "JetConstituentsUtils::sum_tlv_constituents(JetsConstituents)") + .Define("de", "JetConstituentsUtils::compute_residue_energy(tlv_jets, sum_tlv_jcs)") + .Define("dpt", "JetConstituentsUtils::compute_residue_pt(tlv_jets, sum_tlv_jcs)") + .Define("dphi", "JetConstituentsUtils::compute_residue_phi(tlv_jets, sum_tlv_jcs)") + .Define("dtheta", "JetConstituentsUtils::compute_residue_theta(tlv_jets, sum_tlv_jcs)") + + #.Define("Invariant_mass", "ROOT::VecOps::InvariantMasses(jets_ee_genkt[0].pt(), jets_ee_genkt[0].rap(), jets_ee_genkt[0].phi(), jets_ee_genkt[0].m(), jets_ee_genkt[1].pt(), jets_ee_genkt[1].rap(), jets_ee_genkt[1].phi(), jets_ee_genkt[1].m())); + .Define("invariant_mass", "JetConstituentsUtils::InvariantMass(tlv_jets[0], tlv_jets[1])") + ) + return df2 + + def output(): + branchList = [ + #'RP_px', 'RP_py','RP_pz','RP_e', 'RP_m', 'RP_q', + #'Jets_px', 'Jets_py', 'Jets_pz', + 'Jets_e', 'Jets_mass', 'Jets_pt', 'Jets_phi', 'Jets_theta', + 'JetsConstituents_e', 'JetsConstituents_pt', 'JetsConstituents_theta', 'JetsConstituents_phi', 'JetsConstituents_charge', + 'JetsConstituents_erel', 'JetsConstituents_erel_log', 'JetsConstituents_thetarel', 'JetsConstituents_phirel', + 'JetsConstituents_dndx', 'JetsConstituents_mtof', + + 'JetsConstituents_d0_wrt0', 'JetsConstituents_z0_wrt0', 'JetsConstituents_phi0_wrt0', 'JetsConstituents_omega_wrt0', 'JetsConstituents_tanlambda_wrt0', + 'Bz', 'JetsConstituents_Bz', + #'JetsConstituents_Par', + 'JetsConstituents_dxy', 'JetsConstituents_dz', 'JetsConstituents_phi0', 'JetsConstituents_C', 'JetsConstituents_ct', + + 'JetsConstituents_omega_cov', 'JetsConstituents_d0_cov', 'JetsConstituents_z0_cov', 'JetsConstituents_phi0_cov', 'JetsConstituents_tanlambda_cov', + 'JetsConstituents_d0_z0_cov', 'JetsConstituents_phi0_d0_cov', 'JetsConstituents_phi0_z0_cov', + 'JetsConstituents_tanlambda_phi0_cov', 'JetsConstituents_tanlambda_d0_cov', 'JetsConstituents_tanlambda_z0_cov', + 'JetsConstituents_omega_tanlambda_cov', 'JetsConstituents_omega_phi0_cov', 'JetsConstituents_omega_d0_cov', 'JetsConstituents_omega_z0_cov', + 'JetsConstituents_Sip2dVal', 'JetsConstituents_Sip2dSig', + 'JetsConstituents_Sip3dVal', 'JetsConstituents_Sip3dSig', + 'JetsConstituents_JetDistVal', 'JetsConstituents_JetDistSig', + #'JC_Jet0_Pids', + 'JetsConstituents_isMu', + 'JetsConstituents_isEl', + 'JetsConstituents_isChargedHad', + 'JetsConstituents_isGamma', + 'JetsConstituents_isNeutralHad', + 'njet', 'nconst', + 'nmu', 'nel', 'nchargedhad', 'nphoton', 'nneutralhad', + 'de', 'dpt', 'dphi', 'dtheta', + 'invariant_mass' + ] + return branchList diff --git a/examples/FCCee/weaver/stage2.cpp b/examples/FCCee/weaver/stage2.cpp new file mode 100644 index 0000000000..17377a2a53 --- /dev/null +++ b/examples/FCCee/weaver/stage2.cpp @@ -0,0 +1,535 @@ +//standard library header files +#include +#include +#include +#include +#include + + + +//ROOT header files +#include "TH1F.h" +#include "TCanvas.h" +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TStyle.h" +#include +#include "ROOT/RVec.hxx" + + + + +int main(int argc, char* argv[]) { + + //usage + if( argc!= 5 ) { + std::cerr << "USAGE: ./to_jetntuple [root_inFileName] [root_outFileName] N_i N_f" << std::endl; + exit(1); + } + + //std::string inDir = ""; + std::string infileName(argv[1]); + //Opening the input file containing the tree (output of fcc analyses_jets_stage1.py) + TFile* infile = TFile::Open(infileName.c_str()); + if( !infile->IsOpen() ){ + std::cerr << "Problems opening root file. Exiting." << std::endl; + exit(-1); + } + std::cout << "-> Opened file " << infileName.c_str() << std::endl; + //Get pointer to tree object + + //std::cout << "infileName lenght: " << infileName.length() << std::endl; + char flavour = infileName[infileName.length()-6]; + std::cout << "flavour: " << flavour << std::endl; + + TTree* ev = (TTree*)infile->Get("events"); + if(!ev) { + std::cerr << "null pointer for TTree! Exiting." << std::endl; + exit(-2); + } + std::cout << "-> Opened tree " << "events" << std::endl; + //variables to be read from the tree + //event properties + ROOT::VecOps::RVec *Jets_e=0; + ROOT::VecOps::RVec *Jets_mass=0; + ROOT::VecOps::RVec *Jets_pt = 0; + ROOT::VecOps::RVec *Jets_phi = 0; + ROOT::VecOps::RVec *Jets_theta = 0; + + int nJets; + //float Bz_i; + + //properties of constituents + + ROOT::VecOps::RVec > *JetsConstituents_e = 0; + ROOT::VecOps::RVec > *JetsConstituents_pt = 0; + ROOT::VecOps::RVec > *JetsConstituents_theta = 0; + ROOT::VecOps::RVec > *JetsConstituents_phi = 0; + ROOT::VecOps::RVec > *JetsConstituents_charge = 0; + + ROOT::VecOps::RVec > *JetsConstituents_erel = 0; + ROOT::VecOps::RVec > *JetsConstituents_erel_log = 0; + ROOT::VecOps::RVec > *JetsConstituents_thetarel = 0; + ROOT::VecOps::RVec > *JetsConstituents_phirel = 0; + + ROOT::VecOps::RVec > *JetsConstituents_dndx = 0; + ROOT::VecOps::RVec > *JetsConstituents_mtof = 0; + + //ROOT::VecOps::RVec > *JetsConstituents_d0_wrt0 = 0; + //ROOT::VecOps::RVec > *JetsConstituents_z0_wrt0 = 0; + //ROOT::VecOps::RVec > *JetsConstituents_phi0_wrt0 = 0; + //ROOT::VecOps::RVec > *JetsConstituents_omega_wrt0 = 0; + //ROOT::VecOps::RVec > *JetsConstituents_tanlambda_wrt0 = 0; + + ROOT::VecOps::RVec > *JetsConstituents_dxy = 0; + ROOT::VecOps::RVec > *JetsConstituents_dz = 0; + //ROOT::VecOps::RVec > *JetsConstituents_phi0 = 0; + //ROOT::VecOps::RVec > *JetsConstituents_C = 0; + //ROOT::VecOps::RVec > *JetsConstituents_ct = 0; + + ROOT::VecOps::RVec > *JetsConstituents_omega_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_d0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_z0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_phi0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_tanlambda_cov = 0; + + ROOT::VecOps::RVec > *JetsConstituents_d0_z0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_phi0_d0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_phi0_z0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_tanlambda_phi0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_tanlambda_d0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_tanlambda_z0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_omega_tanlambda_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_omega_phi0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_omega_d0_cov = 0; + ROOT::VecOps::RVec > *JetsConstituents_omega_z0_cov = 0; + + ROOT::VecOps::RVec > *JetsConstituents_Sip2dVal = 0; + ROOT::VecOps::RVec > *JetsConstituents_Sip2dSig = 0; + ROOT::VecOps::RVec > *JetsConstituents_Sip3dVal = 0; + ROOT::VecOps::RVec > *JetsConstituents_Sip3dSig = 0; + ROOT::VecOps::RVec > *JetsConstituents_JetDistVal = 0; + ROOT::VecOps::RVec > *JetsConstituents_JetDistSig = 0; + + ROOT::VecOps::RVec > *JetsConstituents_isMu = 0; + ROOT::VecOps::RVec > *JetsConstituents_isEl = 0; + ROOT::VecOps::RVec > *JetsConstituents_isChargedHad = 0; + ROOT::VecOps::RVec > *JetsConstituents_isGamma = 0; + ROOT::VecOps::RVec > *JetsConstituents_isNeutralHad = 0; + + ROOT::VecOps::RVec* count_Const = 0; + ROOT::VecOps::RVec* count_Mu = 0; + ROOT::VecOps::RVec* count_El = 0; + ROOT::VecOps::RVec* count_ChargedHad = 0; + ROOT::VecOps::RVec* count_Photon = 0; + ROOT::VecOps::RVec* count_NeutralHad = 0; + + //Set the info for each branch of the tree to correspond to our data + + ev->SetBranchAddress("Jets_e", &Jets_e); + ev->SetBranchAddress("Jets_mass", &Jets_mass); + ev->SetBranchAddress("Jets_pt", &Jets_pt); + ev->SetBranchAddress("Jets_phi", &Jets_phi); + ev->SetBranchAddress("Jets_theta", &Jets_theta); + + ev->SetBranchAddress("JetsConstituents_e", &JetsConstituents_e); + ev->SetBranchAddress("JetsConstituents_pt", &JetsConstituents_pt); + ev->SetBranchAddress("JetsConstituents_theta", &JetsConstituents_theta); + ev->SetBranchAddress("JetsConstituents_phi", &JetsConstituents_phi); + ev->SetBranchAddress("JetsConstituents_charge", &JetsConstituents_charge); + + ev->SetBranchAddress("JetsConstituents_erel", &JetsConstituents_erel); + ev->SetBranchAddress("JetsConstituents_erel_log", &JetsConstituents_erel_log); + ev->SetBranchAddress("JetsConstituents_thetarel", &JetsConstituents_thetarel); + ev->SetBranchAddress("JetsConstituents_phirel", &JetsConstituents_phirel); + + ev->SetBranchAddress("JetsConstituents_dndx", &JetsConstituents_dndx); + ev->SetBranchAddress("JetsConstituents_mtof", &JetsConstituents_mtof); + + //ev->SetBranchAddress("JetsConstituents_d0_wrt0", &JetsConstituents_d0_wrt0); + //ev->SetBranchAddress("JetsConstituents_z0_wrt0", &JetsConstituents_z0_wrt0); + //ev->SetBranchAddress("JetsConstituents_phi0_wrt0", &JetsConstituents_phi0_wrt0); + //ev->SetBranchAddress("JetsConstituents_omega_wrt0", &JetsConstituents_omega_wrt0); + //ev->SetBranchAddress("JetsConstituents_tanlambda_wrt0", &JetsConstituents_tanlambda_wrt0); + + ev->SetBranchAddress("JetsConstituents_dxy", &JetsConstituents_dxy); + ev->SetBranchAddress("JetsConstituents_dz", &JetsConstituents_dz); + //ev->SetBranchAddress("JetsConstituents_phi0", &JetsConstituents_phi0); + //ev->SetBranchAddress("JetsConstituents_C", &JetsConstituents_C); + //ev->SetBranchAddress("JetsConstituents_ct", &JetsConstituents_ct); + + ev->SetBranchAddress("JetsConstituents_omega_cov", &JetsConstituents_omega_cov); + ev->SetBranchAddress("JetsConstituents_d0_cov", &JetsConstituents_d0_cov); + ev->SetBranchAddress("JetsConstituents_z0_cov", &JetsConstituents_z0_cov); + ev->SetBranchAddress("JetsConstituents_phi0_cov", &JetsConstituents_phi0_cov); + ev->SetBranchAddress("JetsConstituents_tanlambda_cov", &JetsConstituents_tanlambda_cov); + + ev->SetBranchAddress("JetsConstituents_d0_z0_cov", &JetsConstituents_d0_z0_cov); + ev->SetBranchAddress("JetsConstituents_phi0_d0_cov", &JetsConstituents_phi0_d0_cov); + ev->SetBranchAddress("JetsConstituents_phi0_z0_cov", &JetsConstituents_phi0_z0_cov); + ev->SetBranchAddress("JetsConstituents_tanlambda_phi0_cov", &JetsConstituents_tanlambda_phi0_cov); + ev->SetBranchAddress("JetsConstituents_tanlambda_d0_cov", &JetsConstituents_tanlambda_d0_cov); + ev->SetBranchAddress("JetsConstituents_tanlambda_z0_cov", &JetsConstituents_tanlambda_z0_cov); + ev->SetBranchAddress("JetsConstituents_omega_phi0_cov", &JetsConstituents_omega_phi0_cov); + ev->SetBranchAddress("JetsConstituents_omega_d0_cov", &JetsConstituents_omega_d0_cov); + ev->SetBranchAddress("JetsConstituents_omega_z0_cov", &JetsConstituents_omega_z0_cov); + ev->SetBranchAddress("JetsConstituents_omega_tanlambda_cov", &JetsConstituents_omega_tanlambda_cov); + + ev->SetBranchAddress("JetsConstituents_Sip2dVal", &JetsConstituents_Sip2dVal); + ev->SetBranchAddress("JetsConstituents_Sip2dSig", &JetsConstituents_Sip2dSig); + ev->SetBranchAddress("JetsConstituents_Sip3dVal", &JetsConstituents_Sip3dVal); + ev->SetBranchAddress("JetsConstituents_Sip3dSig", &JetsConstituents_Sip3dSig); + ev->SetBranchAddress("JetsConstituents_JetDistVal", &JetsConstituents_JetDistVal); + ev->SetBranchAddress("JetsConstituents_JetDistSig", &JetsConstituents_JetDistSig); + + ev->SetBranchAddress("JetsConstituents_isMu", &JetsConstituents_isMu); + ev->SetBranchAddress("JetsConstituents_isEl", &JetsConstituents_isEl); + ev->SetBranchAddress("JetsConstituents_isChargedHad", &JetsConstituents_isChargedHad); + ev->SetBranchAddress("JetsConstituents_isGamma", &JetsConstituents_isGamma); + ev->SetBranchAddress("JetsConstituents_isNeutralHad", &JetsConstituents_isNeutralHad); + + //ev->SetBranchAddress("Bz", &Bz_i); + + ev->SetBranchAddress("njet", &nJets); + ev->SetBranchAddress("nconst", &count_Const); + ev->SetBranchAddress("nmu", &count_Mu); + ev->SetBranchAddress("nel", &count_El); + ev->SetBranchAddress("nchargedhad", &count_ChargedHad); + ev->SetBranchAddress("nphoton", &count_Photon); + ev->SetBranchAddress("nneutralhad", &count_NeutralHad); + + + //we defined how we read the tree. Now we + //need to define how we write the ntuple. + + //std::string outDir(""); + std::string outfileName(argv[2]); + TFile* outfile = new TFile(outfileName.c_str(),"recreate"); + std::cout << "-> Opened outfile " << std::endl; + TTree* ntuple = new TTree("tree", "jets_Ntuple"); + std::cout << "-> Opened ntuple " << std::endl; + + //variables to write + //jet + double recojet_e, recojet_mass, recojet_pt, recojet_phi, recojet_theta; + + //constituents + float pfcand_e[1000] = {0.}; + float pfcand_pt[1000] = {0.}; + float pfcand_charge[1000] = {0.}; + float pfcand_theta[1000] = {0.}; + float pfcand_phi[1000] = {0.}; + + float pfcand_erel[1000] = {0.}; + float pfcand_erel_log[1000] = {0.}; + float pfcand_thetarel[1000] = {0.}; + float pfcand_phirel[1000] = {0.}; + + float pfcand_dndx[1000] = {0.}; + float pfcand_mtof[1000] = {0.}; + + float pfcand_dxy[1000] = {0.}; + float pfcand_dz[1000] = {0.}; + + float pfcand_dptdpt[1000] = {0.}; + float pfcand_dxydxy[1000] = {0.}; + float pfcand_dzdz[1000] = {0.}; + float pfcand_dphidphi[1000] = {0.}; + float pfcand_detadeta[1000] = {0.}; + + float pfcand_dxydz[1000] = {0.}; + float pfcand_dphidxy[1000] = {0.}; + float pfcand_phidz[1000] = {0.}; + float pfcand_phictgtheta[1000] = {0.}; + float pfcand_dxyctgtheta[1000] = {0.}; + float pfcand_dlambdadz[1000] = {0.}; + float pfcand_cctgtheta[1000] = {0.}; + float pfcand_phic[1000] = {0.}; + float pfcand_dxyc[1000] = {0.}; + float pfcand_cdz[1000] = {0.}; + + float pfcand_btagSip2dVal[1000] = {0.}; + float pfcand_btagSip2dSig[1000] = {0.}; + float pfcand_btagSip3dVal[1000] = {0.}; + float pfcand_btagSip3dSig[1000] = {0.}; + float pfcand_btagJetDistVal[1000] = {0.}; + float pfcand_btagJetDistSig[1000] = {0.}; + + float pfcand_isMu[1000] = {0.}; + float pfcand_isEl[1000] ={0.}; + float pfcand_isChargedHad[1000] ={0.}; + float pfcand_isGamma[1000] ={0.}; + float pfcand_isNeutralHad[1000] ={0.}; + + //counting species + int njet = 0; + int nconst = 0; //number of constituents of the jets + int anomaly_njets_counts_less = 0; + int anomaly_njets_counts_more = 0; + int anomaly_njets_counts = 0; + int saved_events_counts = 0; //this variable will be used to count the number of events actually saved + //int Nevents_Max = 1000000; // maximum number of events to be saved + int nphotons = 0; + int ncharged = 0; + int nchargedhad = 0; + int nneutralhad = 0; + int nel = 0; + int nmu = 0; + + + //set flags + float is_q = 0.; + float is_b = 0.; + float is_c = 0.; + float is_s = 0.; + float is_g = 0.; + float is_t = 0.; + + if (flavour == 'q') {is_q = 1.;} + if (flavour == 'b') {is_b = 1.;} + if (flavour == 'c') {is_c = 1.;} + if (flavour == 's') {is_s = 1.;} + if (flavour == 'g') {is_g = 1.;} + if (flavour == 't') {is_t = 1.;} + + std::cout << "is_q: " << is_q << std::endl; + + // In an n-tuple, we assign each variable to its own branch. + ntuple->Branch("recojet_e", &recojet_e); + ntuple->Branch("recojet_mass", &recojet_mass); + ntuple->Branch("recojet_pt", &recojet_pt); + ntuple->Branch("recojet_phi", &recojet_phi); + ntuple->Branch("recojet_theta", &recojet_theta); + + + ntuple->Branch("recojet_isQ", &is_q); + ntuple->Branch("recojet_isB", &is_b); + ntuple->Branch("recojet_isC", &is_c); + ntuple->Branch("recojet_isS", &is_s); + ntuple->Branch("recojet_isG", &is_g); + ntuple->Branch("recojet_isT", &is_t); + + ntuple->Branch("nconst", &nconst, "nconst/I"); + ntuple->Branch("nphotons", &nphotons, "nphotons/I"); + ntuple->Branch("ncharged", &ncharged, "ncharged/I"); + ntuple->Branch("nneutralhad", &nneutralhad, "nneutralhad/I"); + ntuple->Branch("nchargedhad", &nchargedhad, "nchargedhad/I"); + ntuple->Branch("nel", &nel, "nel/I"); + ntuple->Branch("nmu", &nmu, "nmu/I"); + + ntuple->Branch("pfcand_e", pfcand_e, "pfcand_e[nconst]/F"); + ntuple->Branch("pfcand_pt", pfcand_pt, "pfcand_pt[nconst]/F"); + ntuple->Branch("pfcand_charge", pfcand_charge, "pfcand_charge[nconst]/F"); + ntuple->Branch("pfcand_theta", pfcand_theta, "pfcand_theta[nconst]/F"); + ntuple->Branch("pfcand_phi", pfcand_phi, "pfcand_phi[nconst]/F"); + + ntuple->Branch("pfcand_erel", pfcand_erel, "pfcand_erel[nconst]/F"); + ntuple->Branch("pfcand_erel_log", pfcand_erel_log, "pfcand_erel_log[nconst]/F"); + ntuple->Branch("pfcand_thetarel", pfcand_thetarel, "pfcand_thetarel[nconst]/F"); + ntuple->Branch("pfcand_phirel", pfcand_phirel, "pfcand_phirel[nconst]/F"); + + ntuple->Branch("pfcand_dndx", pfcand_dndx, "pfcand_dndx[nconst]/F"); + ntuple->Branch("pfcand_mtof", pfcand_mtof, "pfcand_mtof[nconst]/F"); + + ntuple->Branch("pfcand_dxy", pfcand_dxy, "pfcand_dxy[nconst]/F"); + ntuple->Branch("pfcand_dz", pfcand_dz, "pfcand_dz[nconst]/F"); + + ntuple->Branch("pfcand_dptdpt", pfcand_dptdpt, "pfcand_dptdpt[nconst]/F"); + ntuple->Branch("pfcand_dxydxy", pfcand_dxydxy, "pfcand_dxydxy[nconst]/F"); + ntuple->Branch("pfcand_dzdz", pfcand_dzdz, "pfcand_dzdz[nconst]/F"); + ntuple->Branch("pfcand_dphidphi", pfcand_dphidphi, "pfcand_dphidphi[nconst]/F"); + ntuple->Branch("pfcand_detadeta", pfcand_detadeta, "pfcand_detadeta[nconst]/F"); + + ntuple->Branch("pfcand_dxydz", pfcand_dxydz, "pfcand_dxydz[nconst]/F"); + ntuple->Branch("pfcand_dphidxy", pfcand_dphidxy, "pfcand_dphidxy[nconst]/F"); + ntuple->Branch("pfcand_phidz", pfcand_phidz, "pfcand_phidz[nconst]/F"); + ntuple->Branch("pfcand_phictgtheta", pfcand_phictgtheta, "pfcand_phictgtheta[nconst]/F"); + ntuple->Branch("pfcand_dxyctgtheta", pfcand_dxyctgtheta, "pfcand_dxyctgtheta[nconst]/F"); + ntuple->Branch("pfcand_dlambdadz", pfcand_dlambdadz, "pfcand_dlambdadz[nconst]/F"); + ntuple->Branch("pfcand_cctgtheta", pfcand_cctgtheta, "pfcand_cctgtheta[nconst]/F"); + ntuple->Branch("pfcand_phic", pfcand_phic, "pfcand_phic[nconst]/F"); + ntuple->Branch("pfcand_dxyc", pfcand_dxyc, "pfcand_dxyc[nconst]/F"); + ntuple->Branch("pfcand_cdz", pfcand_cdz, "pfcand_cdz[nconst]/F"); + + ntuple->Branch("pfcand_btagSip2dVal", pfcand_btagSip2dVal, "pfcand_btagSip2dVal[nconst]/F"); + ntuple->Branch("pfcand_btagSip2dSig", pfcand_btagSip2dSig, "pfcand_btagSip2dSig[nconst]/F"); + ntuple->Branch("pfcand_btagSip3dVal", pfcand_btagSip3dVal, "pfcand_btagSip3dVal[nconst]/F"); + ntuple->Branch("pfcand_btagSip3dSig", pfcand_btagSip3dSig, "pfcand_btagSip3dSig[nconst]/F"); + ntuple->Branch("pfcand_btagJetDistVal", pfcand_btagJetDistVal, "pfcand_btagJetDistVal[nconst]/F"); + ntuple->Branch("pfcand_btagJetDistSig", pfcand_btagJetDistSig, "pfcand_btagJetDistSig[nconst]/F"); + + ntuple->Branch("pfcand_isMu", pfcand_isMu, "pfcand_isMu[nconst]/F"); + ntuple->Branch("pfcand_isEl", pfcand_isEl, "pfcand_isEl[nconst]/F"); + ntuple->Branch("pfcand_isChargedHad", pfcand_isChargedHad, "pfcand_isChargedHad[nconst]/F"); + ntuple->Branch("pfcand_isGamma", pfcand_isGamma, "pfcand_isGamma[nconst]/F"); + ntuple->Branch("pfcand_isNeutralHad", pfcand_isNeutralHad, "pfcand_isNeutralHad[nconst]/F"); + + int N_i = atoi(argv[3]); + int N_f = atoi(argv[4]); + int Nevents_Max = N_f - N_i; // maximum number of events to be saved + + //Run over each entry in the tree: + int nentries = ev->GetEntries(); //total number of events in the file + + std::cout<< " " << std::endl; + std::cout << "-> number of events contained in the tree: " << nentries <GetEntry(i); + + //njet = (*Jets_e).size(); + njet = nJets; + + if(i % 10000 == 0) { + std::cout<< "-----" << std::endl; + std::cout << "-> event: " << i << " - " << "#jets: " << njet << " - (*Jets_e).size(): " << (*Jets_e).size() << std::endl; + std::cout << "-----" << std::endl; + } + + if(njet != 2) { + anomaly_njets_counts += 1; + if (njet > 2) { + anomaly_njets_counts_more += 1; + } else { + anomaly_njets_counts_less += 1; + } + } + + if (njet < 2) { //exclude the events with less than two jets + continue ; + } + + //run over the jets in the event + for(int j=0; j < 2; ++j) { //we only take the two leadingjets + + recojet_e = (*Jets_e)[j]; + //jet_e = jets_e->at(j); + recojet_mass = (*Jets_mass)[j]; + recojet_pt = (*Jets_pt)[j]; + recojet_phi = (*Jets_phi)[j]; + recojet_theta = (*Jets_theta)[j]; + //n_constituents = (*jets_ends)[j] - (*jets_begins)[j]; + //nconst = (JetsConstituents_e->at(j)).size(); + nconst = (count_Const->at(j)); + nel = (count_El->at(j)); + nmu = (count_Mu->at(j)); + nchargedhad = (count_ChargedHad->at(j)); + nphotons = (count_Photon->at(j)); + nneutralhad = (count_NeutralHad->at(j)); + + + if(i % 10000 == 0) { + std::cout << "-> jet: " << j << " - " << "nconst: " << nconst << "-> (JetsConstituents_e->at(j)).size(): " << (JetsConstituents_e->at(j)).size() << std::endl; + std::cout << "-> jet_e: " << recojet_e << " - jet_pt: " << recojet_pt << std::endl; + std::cout<< "-----" << std::endl; + } + + for(int k = 0; k < nconst; ++k){ + pfcand_e[k] = (JetsConstituents_e->at(j))[k]; + //std::cout << k << ' ' << JC_e[k] << std::endl; + pfcand_pt[k] = (JetsConstituents_pt->at(j))[k]; + pfcand_theta[k] = (JetsConstituents_theta->at(j))[k]; + pfcand_phi[k] = (JetsConstituents_phi->at(j))[k]; + pfcand_charge[k] = (JetsConstituents_charge->at(j))[k]; + + pfcand_erel[k] = (JetsConstituents_erel->at(j))[k]; + pfcand_erel_log[k] = (JetsConstituents_erel_log->at(j))[k]; + pfcand_thetarel[k] = (JetsConstituents_thetarel->at(j))[k]; + pfcand_phirel[k] = (JetsConstituents_phirel->at(j))[k]; + + pfcand_dndx[k] = (JetsConstituents_dndx->at(j))[k]/1000.; //transformed in mm + pfcand_mtof[k] = (JetsConstituents_mtof->at(j))[k]; + + pfcand_dxy[k] = (JetsConstituents_dxy->at(j))[k]; + //pfcand_dz[k] = (JetsConstituents_dz->at(j))[k]; + //std::cout<at(j))[k]<at(j))[k])) pfcand_dz[k] = -9; + else pfcand_dz[k] = (JetsConstituents_dz->at(j))[k]; + pfcand_dptdpt[k] = (JetsConstituents_omega_cov->at(j))[k]; + pfcand_dxydxy[k] = (JetsConstituents_d0_cov->at(j))[k]; + pfcand_dzdz[k] = (JetsConstituents_z0_cov->at(j))[k]; + pfcand_dphidphi[k] = (JetsConstituents_phi0_cov->at(j))[k]; + pfcand_detadeta[k] = (JetsConstituents_tanlambda_cov->at(j))[k]; + + pfcand_dxydz[k] = (JetsConstituents_d0_z0_cov->at(j))[k]; + pfcand_dphidxy[k] = (JetsConstituents_phi0_d0_cov->at(j))[k]; //***** + pfcand_phidz[k] = (JetsConstituents_phi0_z0_cov->at(j))[k]; + + pfcand_phictgtheta[k] = (JetsConstituents_tanlambda_phi0_cov->at(j))[k]; + pfcand_dxyctgtheta[k] = (JetsConstituents_tanlambda_d0_cov->at(j))[k]; + pfcand_dlambdadz[k] = (JetsConstituents_tanlambda_z0_cov->at(j))[k]; + + pfcand_cctgtheta[k] = (JetsConstituents_omega_tanlambda_cov->at(j))[k]; + pfcand_phic[k] = (JetsConstituents_omega_phi0_cov->at(j))[k]; + pfcand_dxyc[k] = (JetsConstituents_omega_d0_cov->at(j))[k]; + pfcand_cdz[k] = (JetsConstituents_omega_z0_cov->at(j))[k]; + + pfcand_btagSip2dVal[k] = (JetsConstituents_Sip2dVal->at(j))[k]; + pfcand_btagSip2dSig[k] = (JetsConstituents_Sip2dSig->at(j))[k]; + //pfcand_btagSip3dVal[k] = (JetsConstituents_Sip3dVal->at(j))[k]; + //pfcand_btagSip3dSig[k] = (JetsConstituents_Sip3dSig->at(j))[k]; + //pfcand_btagJetDistVal[k] = (JetsConstituents_JetDistVal->at(j))[k]; + //pfcand_btagJetDistSig[k] = (JetsConstituents_JetDistSig->at(j))[k]; + + if (isnan((JetsConstituents_Sip3dVal->at(j))[k]))pfcand_btagSip3dVal[k] = -9; + else pfcand_btagSip3dVal[k] = JetsConstituents_Sip3dVal->at(j)[k]; + if (isnan((JetsConstituents_Sip3dSig->at(j))[k]))pfcand_btagSip3dSig[k] = -9; + else pfcand_btagSip3dSig[k] = JetsConstituents_Sip3dSig->at(j)[k]; + if (isnan((JetsConstituents_JetDistVal->at(j))[k]))pfcand_btagJetDistVal[k] = -9; + else pfcand_btagJetDistVal[k] = JetsConstituents_JetDistVal->at(j)[k]; + if (isnan((JetsConstituents_JetDistSig->at(j))[k]))pfcand_btagJetDistSig[k] = -9; + else pfcand_btagJetDistSig[k] = JetsConstituents_JetDistSig->at(j)[k]; + + pfcand_isMu[k] = (JetsConstituents_isMu->at(j))[k]; + pfcand_isEl[k] = (JetsConstituents_isEl->at(j))[k]; + pfcand_isChargedHad[k] = (JetsConstituents_isChargedHad->at(j))[k]; + pfcand_isGamma[k] = (JetsConstituents_isGamma->at(j))[k]; + pfcand_isNeutralHad[k] = (JetsConstituents_isNeutralHad->at(j))[k]; + + //std::cout << "k: "<< k << std::endl; + //counting species + /*if ( (JetsConstituents_isMu->at(j))[k] == 1 ) { + nmu += 1; + ncharged += 1; + } else if ((JetsConstituents_isEl->at(j))[k] == 1) { + nel += 1; + ncharged += 1; + } else if ( (JetsConstituents_isGamma->at(j))[k] == 1 ) { + nphotons += 1; + } else if ((JetsConstituents_isChargedHad->at(j))[k] == 1) { + nchargedhad += 1; + ncharged += 1; + } else if ((JetsConstituents_isNeutralHad->at(j))[k] == 1) { + nneutralhad += 1; + }*/ + } + ntuple->Fill(); + } + saved_events_counts += 1; //we count the num of events saved + if (saved_events_counts == Nevents_Max) { //interrupt the loop if Nevents_max events have already been saved + break; + } + } + + std::cout << "-> number of entries run: " << nentries < number of entries considered: " << saved_events_counts < number of events with njets != 2: " << anomaly_njets_counts << std::endl; + std::cout << "--> number of events with njets > 2: " << anomaly_njets_counts_more << std::endl; + std::cout << "--> number of events with njets < 2: " << anomaly_njets_counts_less << std::endl; + + + //outfile->cd(); + //outfile->mkdir("deepntuplizer/"); + //outfile->cd("deepntuplizer/"); + //ntuple->SetDirectory(gDirectory); + ntuple->Write(); + infile->Close(); + outfile->Close(); + std::cout << "-> Closed files "<