diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 47d126f..1df546a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -143,6 +143,7 @@ include: - local: 'benchmarks/timing/config.yml' - local: 'benchmarks/b0_tracker/config.yml' - local: 'benchmarks/insert_muon/config.yml' + - local: 'benchmarks/insert_tau/config.yml' - local: 'benchmarks/zdc_sigma/config.yml' - local: 'benchmarks/zdc_lambda/config.yml' - local: 'benchmarks/insert_neutron/config.yml' @@ -170,6 +171,7 @@ deploy_results: - "collect_results:zdc" - "collect_results:zdc_lyso" - "collect_results:insert_muon" + - "collect_results:insert_tau" - "collect_results:zdc_photon" - "collect_results:zdc_pi0" - "collect_results:femc_electron" diff --git a/Snakefile b/Snakefile index 1390e86..9c12548 100644 --- a/Snakefile +++ b/Snakefile @@ -8,13 +8,14 @@ include: "benchmarks/material_scan/Snakefile" include: "benchmarks/tracking_performances/Snakefile" include: "benchmarks/tracking_performances_dis/Snakefile" include: "benchmarks/lfhcal/Snakefile" +include: "benchmarks/zdc_lyso/Snakefile" include: "benchmarks/insert_muon/Snakefile" include: "benchmarks/zdc_lambda/Snakefile" -include: "benchmarks/zdc_lyso/Snakefile" include: "benchmarks/zdc_photon/Snakefile" include: "benchmarks/zdc_pi0/Snakefile" include: "benchmarks/zdc_sigma/Snakefile" include: "benchmarks/insert_neutron/Snakefile" +include: "benchmarks/insert_tau/Snakefile" include: "benchmarks/femc_electron/Snakefile" include: "benchmarks/femc_photon/Snakefile" include: "benchmarks/femc_pi0/Snakefile" diff --git a/benchmarks/insert_tau/Snakefile b/benchmarks/insert_tau/Snakefile new file mode 100644 index 0000000..57e4390 --- /dev/null +++ b/benchmarks/insert_tau/Snakefile @@ -0,0 +1,68 @@ +def get_n_events(wildcards): + energy = float(wildcards.P) + n_events = 1000 + n_events = int(n_events // ((energy / 20) ** 0.5)) + return n_events + +rule insert_tau_generate: + input: + script="benchmarks/insert_tau/analysis/gen_particles.cxx", + params: + N_EVENTS=get_n_events, + th_max=7.0, + th_min=1.7, + output: + GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc", + shell: + """ +root -l -b -q '{input.script}({params.N_EVENTS},"{output.GEN_FILE}", "tau-", {params.th_min}, {params.th_max}, 0., 360., {wildcards.P})' +""" + +rule insert_tau_simulate: + input: + GEN_FILE="sim_output/insert_tau/tau-_{P}GeV.hepmc", + warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root", + params: + N_EVENTS=get_n_events, + PHYSICS_LIST="FTFP_BERT", + output: + SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV_{INDEX}.edm4hep.root" + shell: + """ +# Running simulation +npsim \ + --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \ + --numberOfEvents {params.N_EVENTS} \ + --skipNEvents $(( {params.N_EVENTS} * {wildcards.INDEX} )) \ + --physicsList {params.PHYSICS_LIST} \ + --inputFiles {input.GEN_FILE} \ + --outputFile {output.SIM_FILE} +""" + +rule insert_tau_recon: + input: + SIM_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_sim_tau-_{P}GeV_{INDEX}.edm4hep.root", + params: + N_EVENTS=get_n_events, + output: + REC_FILE="sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV_{INDEX}.edm4eic.root", + shell: + """ +eicrecon {input.SIM_FILE} -Ppodio:output_file={output.REC_FILE} -Pdd4hep:xml_files=$DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml -Ppodio:output_collections=MCParticles,HcalEndcapPInsertRecHits,HcalEndcapPInsertClusters,HcalEndcapPInsertSubcellHits,EcalEndcapPInsertRecHits,EcalEndcapPInsertClusters,EcalEndcapPClusters,LFHCALClusters -Pjana:nevents={params.N_EVENTS} +""" + +rule insert_tau_analysis: + input: + expand("sim_output/insert_tau/{DETECTOR_CONFIG}_rec_tau-_{P}GeV_{INDEX}.edm4eic.root", + P=[20, 30, 40, 50, 60, 80, 100], + DETECTOR_CONFIG=["{DETECTOR_CONFIG}"], + INDEX=range(5), + ), + script="benchmarks/insert_tau/analysis/tau_plots.py", + output: + results_dir=directory("results/{DETECTOR_CONFIG}/insert_tau"), + shell: + """ +mkdir -p {output.results_dir} +python {input.script} {output.results_dir} +""" diff --git a/benchmarks/insert_tau/analysis/gen_particles.cxx b/benchmarks/insert_tau/analysis/gen_particles.cxx new file mode 100644 index 0000000..0877a85 --- /dev/null +++ b/benchmarks/insert_tau/analysis/gen_particles.cxx @@ -0,0 +1,127 @@ +#include "HepMC3/GenEvent.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/WriterAscii.h" +#include "HepMC3/Print.h" + +#include "TRandom3.h" +#include "TVector3.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace HepMC3; + +// Generate single electron as input to the Insert simulation. +// -- +// We generate events with a constant polar angle with respect to +// the proton direction and then rotate so that the events are given +// in normal lab coordinate system +// -- +void gen_particles( + int n_events = 1000, + const char* out_fname = "gen_particles.hepmc", + TString particle_name = "e-", + double th_min = 3., // Minimum polar angle, in degrees + double th_max = 3., // Maximum polar angle, in degrees + double phi_min = 0., // Minimum azimuthal angle, in degrees + double phi_max = 360., // Maximum azimuthal angle, in degrees + double p = 10., // Momentum in GeV/c + int dist = 0, //Momentum distribution: 0=fixed, 1=uniform, 2=Gaussian + int useCrossingAngle=1 // 0= no rotation, 1 = -25 mrad + ) +{ + WriterAscii hepmc_output(out_fname); + int events_parsed = 0; + GenEvent evt(Units::GEV, Units::MM); + + // Random number generator + TRandom3 *r1 = new TRandom3(0); //Use time as random seed + + // Getting generated particle information + TDatabasePDG *pdg = new TDatabasePDG(); + TParticlePDG *particle = pdg->GetParticle(particle_name); + const double mass = particle->Mass(); + const int pdgID = particle->PdgCode(); + + for (events_parsed = 0; events_parsed < n_events; events_parsed++) { + + //Set the event number + evt.set_event_number(events_parsed); + + // FourVector(px,py,pz,e,pdgid,status) + // type 4 is beam + // pdgid 11 - electron + // pdgid 111 - pi0 + // pdgid 2212 - proton + GenParticlePtr p1 = + std::make_shared(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared( + FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4); + + // Define momentum with respect to proton direction + double phi = r1->Uniform(phi_min*TMath::DegToRad(),phi_max*TMath::DegToRad()); + double th = r1->Uniform(th_min*TMath::DegToRad(),th_max*TMath::DegToRad()); + + //Total momentum distribution + double pevent = -1; + if(dist==0){ //fixed + pevent = p; + } + else if(dist==1){ //Uniform: +-50% variation + pevent = p*(1. + r1->Uniform(-0.5,0.5) ); + } + else if(dist==2){ //Gaussian: Sigma = 0.1*mean + while(pevent<0) //Avoid negative values + pevent = r1->Gaus(p,0.1*p); + } + + double px = pevent * std::cos(phi) * std::sin(th); + double py = pevent * std::sin(phi) * std::sin(th); + double pz = pevent * std::cos(th); + TVector3 pvec(px,py,pz); + + //Rotate to lab coordinate system + double cross_angle = -25./1000.*useCrossingAngle; //in Rad + TVector3 pbeam_dir(sin(cross_angle),0,cos(cross_angle)); //proton beam direction + pvec.RotateY(-pbeam_dir.Theta()); // Theta is returned positive, beam in negative X + // type 1 is final state + // pdgid 11 - electron 0.510 MeV/c^2 + GenParticlePtr p3 = std::make_shared( + FourVector( + pvec.X(), pvec.Y(), pvec.Z(), + sqrt(pevent*pevent + (mass * mass))), + pdgID, 1); + + //If wanted, set non-zero vertex + double vx = 0.; + double vy = 0.; + double vz = 0.; + double vt = 0.; + + GenVertexPtr v1 = std::make_shared(); + evt.shift_position_by(FourVector(vx, vy, vz, vt)); + v1->add_particle_in(p1); + v1->add_particle_in(p2); + + v1->add_particle_out(p3); + evt.add_vertex(v1); + + if (events_parsed == 0) { + std::cout << "First event: " << std::endl; + Print::listing(evt); + } + + hepmc_output.write_event(evt); + if (events_parsed % 100 == 0) { + std::cout << "Event: " << events_parsed << std::endl; + } + evt.clear(); + } + hepmc_output.close(); + std::cout << "Events parsed and written: " << events_parsed << std::endl; +} diff --git a/benchmarks/insert_tau/analysis/tau_plots.py b/benchmarks/insert_tau/analysis/tau_plots.py new file mode 100644 index 0000000..77c2dd3 --- /dev/null +++ b/benchmarks/insert_tau/analysis/tau_plots.py @@ -0,0 +1,154 @@ +import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib as mpl, awkward as ak, sys +import mplhep as hep +hep.style.use("CMS") + +plt.rcParams['figure.facecolor']='white' +plt.rcParams['savefig.facecolor']='white' +plt.rcParams['savefig.bbox']='tight' + +plt.rcParams["figure.figsize"] = (7, 7) + +config=sys.argv[1].split("/")[1] #results/{config}/{benchmark_name} +outdir=sys.argv[1]+"/" +try: + import os + os.mkdir(outdir[:-1]) +except: + pass + +import uproot as ur +arrays_sim={} +momenta=20, 30, 40, 50, 60,80,100 +for p in momenta: + arrays_sim[p] = ur.concatenate({ + f'sim_output/insert_tau/{config}_rec_tau-_{p}GeV_{index}.edm4eic.root': 'events' + for index in range(5) + }) + + +for a in arrays_sim.values(): + #recon + Etot=0 + px=0 + py=0 + pz=0 + for det in "HcalEndcapPInsert", "EcalEndcapPInsert", "EcalEndcapP", "LFHCAL": + + E=a[f'{det}Clusters.energy'] + + #todo apply corrections depending on whether this is an electromagnetic or hadronic shower. + + x=a[f'{det}Clusters.position.x'] + y=a[f'{det}Clusters.position.y'] + z=a[f'{det}Clusters.position.z'] + r=np.sqrt(x**2+y**2+z**2) + Etot=Etot+np.sum(E, axis=-1) + px=px+np.sum(E*x/r,axis=-1) + py=py+np.sum(E*y/r,axis=-1) + pz=pz+np.sum(E*z/r,axis=-1) + + a['jet_p_recon']=np.sqrt(px**2+py**2+pz**2) + a['jet_E_recon']=Etot + + a['jet_theta_recon']=np.arctan2(np.hypot(px*np.cos(-.025)-pz*np.sin(-.025),py), + pz*np.cos(-.025)+px*np.sin(-.025)) + + #truth + m=a['MCParticles.mass'][::,2:] + px=a['MCParticles.momentum.x'][::,2:] + py=a['MCParticles.momentum.y'][::,2:] + pz=a['MCParticles.momentum.z'][::,2:] + E=np.sqrt(m**2+px**2+py**2+pz**2) + status=a['MCParticles.simulatorStatus'][::,2:] + PDG=a['MCParticles.PDG'][::,2:] + + #find the hadronic part: initial-state tau - all leptons + selection=1*(PDG==15)-1*(np.abs(PDG)==16) + is_hadronic=1*(np.sum((PDG==-14)+(PDG==-12), axis=-1)==0) + + px_hfs, py_hfs, pz_hfs= np.sum(px*selection,axis=-1)*is_hadronic,np.sum(py*selection,axis=-1)*is_hadronic, np.sum(pz*selection,axis=-1)*is_hadronic + + a['hfs_p_truth']=np.sqrt(px_hfs**2+py_hfs**2+pz_hfs**2) + a['hfs_E_truth']=np.sum(E*selection,axis=-1)*is_hadronic + + + a['hfs_theta_truth']=np.arctan2(np.hypot(px_hfs*np.cos(-.025)-pz_hfs*np.sin(-.025),py_hfs), + pz_hfs*np.cos(-.025)+px_hfs*np.sin(-.025)) + a['hfs_eta_truth']=-np.log(np.tan(a['hfs_theta_truth']/2)) + a['n_mu']=np.sum(np.abs(PDG)==13, axis=-1) + a['n_e']=np.sum(np.abs(PDG)==13, axis=-1) + a['hfs_mass_truth']=np.sqrt(a['hfs_E_truth']**2-a['hfs_p_truth']**2) + +for a in arrays_sim.values(): + selection=(a['hfs_eta_truth']>3.1) & (a['hfs_eta_truth']<3.8)\ + &(a['n_mu']==0)&(a['n_e']==0)&(a['hfs_mass_truth']>.140)&(a['jet_E_recon']>0) + +Etruth=[] +Erecon=[] + +theta_truth=[] +theta_recon=[] + +eta_max=3.7 +eta_min=3.3 +for a in arrays_sim.values(): + selection=(a['hfs_eta_truth']>eta_min) & (a['hfs_eta_truth'].140)&(a['jet_E_recon']>1) + theta_truth=np.concatenate((theta_truth,a['hfs_theta_truth'][selection])) + theta_recon=np.concatenate((theta_recon,a['jet_theta_recon'][selection])) + Etruth=np.concatenate((Etruth,a['hfs_E_truth'][selection])) + Erecon=np.concatenate((Erecon,a['jet_E_recon'][selection])) + +plt.figure() +plt.scatter(theta_truth, theta_recon, 1) +plt.xlabel("$\\theta^{hfs}_{\\rm truth}$ [rad]") +plt.ylabel("$\\theta^{hfs}_{\\rm rec}$ [rad]") +plt.title(f"$E_{{\\tau}}$=20-100 GeV, ${eta_min}<\\eta_{{hfs}}<{eta_max}$") +plt.plot([0.04,0.1], [0.04, 0.1], color='tab:orange') +plt.ylim(0, 0.15) +plt.savefig(outdir+"/theta_scatter.pdf") + +plt.figure() +plt.scatter(Etruth, Erecon, 1) +plt.xlabel("$E^{hfs}_{\\rm truth}$ [GeV]") +plt.ylabel("$E^{hfs}_{\\rm rec}$ [GeV]") +plt.title(f"$E_{{\\tau}}$=20-100 GeV, ${eta_min}<\\eta_{{hfs}}<{eta_max}$") +plt.plot((0,100), (0, 100), color='tab:orange') +plt.savefig(outdir+"/energy_scatter.pdf") + +def gauss(x, A,mu, sigma): + return A * np.exp(-(x-mu)**2/(2*sigma**2)) +from scipy.optimize import curve_fit + +res=[] +dres=[] +emid=[] +ew=[] +partitions=(20,30, 40, 60,80,100) +for emin, emax in zip(partitions[:-1], partitions[1:]): + + y,x = np.histogram((theta_recon-theta_truth)[(Etruth>emin)&(Etruth