Fixed refractive index for dRICH residuals

benchmarks/rich/src/CherenkovPIDAnalysis.cc

@@ -231,7 +231,10 @@ namespace benchmarks {
 
       // 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`
+      float mc_rindex=0; 
+      if(m_rad_name=="Aerogel") mc_rindex =1.019; 
+      else if(m_rad_name=="Gas") mc_rindex = 1.00076; 
+      else mc_rindex = cherenkov_pid.getRefractiveIndex(); // average 
       auto theta_exp = TMath::ACos(
           TMath::Hypot( charged_particle_momentum, charged_particle_mass ) /
           ( mc_rindex * charged_particle_momentum )
@@ -239,6 +242,20 @@ namespace benchmarks {
       auto theta_resid = theta_rec - theta_exp;
       auto theta_resid_mrad = theta_resid * 1e3; // [rad] -> [mrad]
 
+      for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons()){
+        m_photonTheta_vs_photonPhi->Fill( phi, theta*1e3);
+        m_speTheta_dist->Fill(theta*1e3);
+        m_speResid_dist->Fill((theta-theta_exp)*1e3);
+        m_spetheta_vs_p->Fill(charged_particle_momentum,theta*1e3);
+        m_spethetaResid_vs_p->Fill(charged_particle_momentum,(theta-theta_exp)*1e3);
+        m_spetheta_vs_eta->Fill(charged_particle_eta,theta*1e3);
+        m_spethetaResid_vs_eta->Fill(charged_particle_eta,(theta-theta_exp)*1e3);
+        m_speTheta_vs_mcWavelength->Fill(Tools::HC/ cherenkov_pid.getPhotonEnergy(),theta*1e3 );
+        m_speResid_vs_mcWavelength->Fill(Tools::HC/ cherenkov_pid.getPhotonEnergy(),(theta-theta_exp)*1e3 );
+      } // [rad] -> [mrad]
+
+      m_speTheta_dist->Fit("gaus");
+
       // fill PID histograms
       m_npe_dist->Fill(   cherenkov_pid.getNpe());
       m_npe_vs_p->Fill(   charged_particle_momentum, cherenkov_pid.getNpe());
@@ -249,6 +266,7 @@ namespace benchmarks {
       m_thetaResid_dist->Fill(   theta_resid_mrad);
       m_thetaResid_vs_p->Fill(   charged_particle_momentum, theta_resid_mrad);
       m_thetaResid_vs_eta->Fill( charged_particle_eta,      theta_resid_mrad);
+      /*
       for(const auto& [theta,phi] : cherenkov_pid.getThetaPhiPhotons()){
         m_photonTheta_vs_photonPhi->Fill( phi, theta*1e3);
         m_speTheta_dist->Fill(theta*1e3);
@@ -257,9 +275,11 @@ namespace benchmarks {
         m_spethetaResid_vs_p->Fill(charged_particle_momentum,(theta-theta_exp)*1e3);
         m_spetheta_vs_eta->Fill(charged_particle_eta,theta*1e3);
         m_spethetaResid_vs_eta->Fill(charged_particle_eta,(theta-theta_exp)*1e3);
-	m_speTheta_vs_mcWavelength->Fill(Tools::HC/ cherenkov_pid.getPhotonEnergy(),theta*1e3 );
-	m_speResid_vs_mcWavelength->Fill(Tools::HC/ cherenkov_pid.getPhotonEnergy(),(theta-theta_exp)*1e3 );
+      	m_speTheta_vs_mcWavelength->Fill(Tools::HC/ cherenkov_pid.getPhotonEnergy(),theta*1e3 );
+	      m_speResid_vs_mcWavelength->Fill(Tools::HC/ cherenkov_pid.getPhotonEnergy(),(theta-theta_exp)*1e3 );
       } // [rad] -> [mrad]
+      */
+      
       m_mcWavelength_dist->Fill(  Tools::HC / cherenkov_pid.getPhotonEnergy() ); // energy [GeV] -> wavelength [nm]
       m_mcRindex_dist->Fill( mc_rindex);     
       // find the PDG hypothesis with the highest weight