Multi part debug

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;
 

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);
   }
 
 }

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;