diff --git a/benchmarks/rich/include/ChargedParticle.h b/benchmarks/rich/include/ChargedParticle.h index 607557871c1f0802bca1095c43631230b92828dd..873b65fc9367071c4f6da70bb25a2f3d58711dc0 100644 --- a/benchmarks/rich/include/ChargedParticle.h +++ b/benchmarks/rich/include/ChargedParticle.h @@ -19,8 +19,9 @@ namespace benchmarks { // accessors float GetMomentum() { return m_momentum; }; - float GetPDG() { return m_pdg; }; - + int GetPDG() { return m_pdg; }; + int GetGeneratorStatus() {return m_generatorstatus;}; + int GetThrown(){return m_thrown;}; // set info from a single MCParticle void SetSingleParticle(const edm4hep::MCParticle& mc_part); // set info from a (thrown) MCParticle of a collection @@ -31,7 +32,8 @@ namespace benchmarks { // members float m_momentum = 0.0; int m_pdg = 0; - + int m_generatorstatus=0; + int m_thrown = 0; // logger std::shared_ptr<spdlog::logger> m_log; diff --git a/benchmarks/rich/src/ChargedParticle.cc b/benchmarks/rich/src/ChargedParticle.cc index 90531d06762fdf1d2a78b1f6727453d654682e63..dda1918d1565f5b5a89c54027a0608501b0ac3cd 100644 --- a/benchmarks/rich/src/ChargedParticle.cc +++ b/benchmarks/rich/src/ChargedParticle.cc @@ -9,6 +9,7 @@ namespace benchmarks { void ChargedParticle::SetSingleParticle(const edm4hep::MCParticle& mc_part) { m_momentum = edm4hep::utils::p(mc_part); m_pdg = mc_part.getPDG(); + m_generatorstatus = mc_part.getGeneratorStatus(); } // set info from a (thrown) MCParticle of a collection @@ -19,9 +20,10 @@ namespace benchmarks { this->SetSingleParticle(mc_part); n_thrown++; } + m_thrown = n_thrown; } if(n_thrown == 0) m_log->warn("ChargedParticle::SetSingleParticle: no particle found"); - if(n_thrown > 1) m_log->warn("ChargedParticle::SetSingleParticle: found {} particles (more than one)", n_thrown); + if(n_thrown > 3) m_log->warn("ChargedParticle::SetSingleParticle: found {} particles (more than one)", n_thrown); } } diff --git a/benchmarks/rich/src/CherenkovPIDAnalysis.cc b/benchmarks/rich/src/CherenkovPIDAnalysis.cc index 84b424c73016aa587b9fc64e66bc5174688db577..8b04a71b57384935eda25ef31ca93e29129b088a 100644 --- a/benchmarks/rich/src/CherenkovPIDAnalysis.cc +++ b/benchmarks/rich/src/CherenkovPIDAnalysis.cc @@ -113,13 +113,14 @@ namespace benchmarks { ) { m_log->trace("-> {} Radiator:", m_rad_name); - + //printf("Rad name: %s\n", m_rad_name.c_str()); // get thrown momentum from a single MCParticle auto part = std::make_unique<ChargedParticle>(m_log); part->SetSingleParticle(mc_parts); auto thrown_momentum = part->GetMomentum(); auto thrown_pdg = part->GetPDG(); - + auto thrown_genstat = part->GetGeneratorStatus(); + auto thrown_part = part->GetThrown(); // loop over `CherenkovParticleID` objects for(const auto& cherenkov_pid : cherenkov_pids) { @@ -130,48 +131,77 @@ namespace benchmarks { } // estimate the charged particle momentum using the momentum of the first TrackPoint at this radiator's entrance - /* + auto charged_particle = cherenkov_pid.getChargedParticle(); if(!charged_particle.isAvailable()) { m_log->warn("Charged particle not available in this radiator"); continue; } if(charged_particle.points_size()==0) { m_log->warn("Charged particle has no TrackPoints in this radiator"); continue; } auto charged_particle_momentum = edm4hep::utils::magnitude( charged_particle.getPoints(0).momentum ); m_log->trace(" Charged Particle p = {} GeV at radiator entrance", charged_particle_momentum); - m_log->trace(" If it is a pion, E = {} GeV", std::hypot(charged_particle_momentum, Tools::GetPDGMass(211))); - */ + m_log->trace(" If it is a kaonon, E = {} GeV", std::hypot(charged_particle_momentum, Tools::GetPDGMass(211))); + // alternatively: use `thrown_momentum` (FIXME: will not work for multi-track events) - auto charged_particle_momentum = thrown_momentum; + //auto charged_particle_momentum = thrown_momentum; auto charged_particle_pdg = thrown_pdg; auto charged_particle_mass = Tools::GetPDGMass(charged_particle_pdg); m_log->trace(" Charged Particle PDG={}, mass={} GeV, p={} GeV from truth", charged_particle_pdg, charged_particle_mass, charged_particle_momentum); // get average Cherenkov angle: `theta_rec` - double theta_rec = 0.0; - for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons()) - theta_rec += theta; - theta_rec /= cherenkov_pid.getNpe(); - + double theta_rec = 0.0; double Npe = 0.0; + TH1D * htmp = new TH1D("htmp","htmp",1500,0,300); + for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons()){ + + if(m_rad_name == "Gas") + if(theta*1e3 > 10 && theta*1e3<=45){ + htmp->Fill(theta*1e3); + //theta_rec += theta; + Npe++; + } + if(m_rad_name == "Aerogel") + if(theta*1e3 > 150 && theta*1e3<=220){ + htmp->Fill(theta*1e3); + //theta_rec += theta; + Npe++; + } + + }// theta,phi + theta_rec = htmp->GetMean(); + htmp->Delete(); + + + //theta_rec /= cherenkov_pid.getNpe(); + //theta_rec /= Npe; //cherenkov_pid.getNpe(); + //Npe =0.0; // calculate expected Cherenkov angle `theta_exp` and residual `theta_resid`, // using refractive index from MC truth - auto mc_rindex = cherenkov_pid.getRefractiveIndex(); // average refractive index for photons used in this `cherenkov_pid` + //auto mc_rindex = cherenkov_pid.getRefractiveIndex(); // average refractive index for photons used in this `cherenkov_pid` + double mc_rindex =0. ; + if(m_rad_name == "Gas") mc_rindex = 1.00076; + else if(m_rad_name == "Aerogel") mc_rindex = 1.019; auto theta_exp = TMath::ACos( TMath::Hypot( charged_particle_momentum, charged_particle_mass ) / ( mc_rindex * charged_particle_momentum ) ); - auto theta_resid = theta_rec - theta_exp; + if(m_rad_name == "Gas") printf("---> %lf %lf\n",theta_exp*1e3,theta_rec); + printf("---> %d %d %d %lf\n",thrown_genstat,thrown_pdg,thrown_part,thrown_momentum); + auto theta_resid = 0.; /*if(charged_particle_pdg == thrown_pdg)*/ + theta_resid = theta_rec - 1e3*theta_exp; + //printf("-->Index Mom: %s %7.2f %7.2f\n",m_rad_name.c_str(),mc_rindex,charged_particle_momentum); + //printf("-->Angle Res: %7.2f %7.2f %7.2f\n",theta_resid*1e3,theta_rec*1e3,theta_exp*1e3); // fill PID histograms - m_npe_dist->Fill(cherenkov_pid.getNpe()); + //m_npe_dist->Fill(cherenkov_pid.getNpe()); + m_npe_dist->Fill(Npe); m_npe_vs_p->Fill(charged_particle_momentum,cherenkov_pid.getNpe()); - m_theta_dist->Fill(theta_rec*1e3); // [rad] -> [mrad] + m_theta_dist->Fill(theta_rec);//*1e3); // [rad] -> [mrad] m_theta_vs_p->Fill(charged_particle_momentum,theta_rec*1e3); // [rad] -> [mrad] - m_thetaResid_dist->Fill(theta_resid*1e3); // [rad] -> [mrad] + m_thetaResid_dist->Fill(theta_resid);//*1e3); // [rad] -> [mrad] m_thetaResid_vs_p->Fill(charged_particle_momentum,theta_resid*1e3); // [rad] -> [mrad] for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons()) m_photonTheta_vs_photonPhi->Fill(phi,theta*1e3); // [rad] -> [mrad] m_mcWavelength_dist->Fill( Tools::HC / cherenkov_pid.getPhotonEnergy() ); // energy [GeV] -> wavelength [nm] m_mcRindex_dist->Fill(mc_rindex); - + Npe=0; // find the PDG hypothesis with the highest weight float max_weight = -1000; int pdg_max_weight = 0;