Normal incidence

src/test_pixel_gap_cuts.cpp

@@ -5,15 +5,27 @@
 #include <TApplication.h>
 #include <TH2D.h>
 #include <TCanvas.h>
-#include <TStyle.h>
+#include "TSystem.h"
+#include "TStyle.h"
+#include "TRegexp.h"
+#include "TCanvas.h"
+#include "TApplication.h"
+#include "TBox.h"
+#include "ROOT/RDataFrame.hxx"
+
 
 #include <podio/ROOTFrameReader.h>
 #include <podio/Frame.h>
 #include "DDRec/CellIDPositionConverter.h"
 
 #include <edm4eic/MCRecoTrackerHitAssociationCollection.h>
-
+#include <edm4hep/MCParticleCollection.h>
 #include <services/geometry/richgeo/ReadoutGeo.h>
+#include <services/geometry/richgeo/IrtGeoDRICH.h>
+
+using namespace ROOT;
+using namespace ROOT::VecOps;
+using namespace edm4hep;
 
 int main(int argc, char** argv) {
 
@@ -50,28 +62,51 @@ int main(int argc, char** argv) {
 
   // ReadoutGeo
   richgeo::ReadoutGeo geo("DRICH", &det, &cellid_converter, logger);
+  richgeo::IrtGeoDRICH drichGeo(&det, &cellid_converter, logger);
 
   // open input file
   auto reader = podio::ROOTFrameReader();
   reader.openFile(root_file_name);
   const std::string tree_name = "events";
 
   // local hits histogram
+  double pi = TMath::Pi();  
   auto h = new TH2D("h","local MC SiPM hits",10000,-15,15,10000,-15,15);
-
+  auto h1 = new TH1D("h1","Photon Incidence angle (filtered ring); degrees",180,0,90);
+  auto h2 = new TH2D("h2","Photon Incidence angle Vs Lambda(all);#lambda(nm);angle(degrees)",2000,0.,1000.,180,0.,90.);
+  auto h3 = new TH2D("h3","Photon Incidence angle Vs Lambda(X filtered);#lambda(nm);angle(degrees)",2000,0.,1000.,180,0.,90.);
+  auto h4 = new TH2D("h4","sim hits; X;Y",4000,-2000,2000,4000,-2000,2000);
   // event loop
   for(unsigned e=0; e<reader.getEntries(tree_name); e++) {
     logger->trace("EVENT {}", e);
     auto frame = podio::Frame(reader.readNextEntry(tree_name));
     const auto& hit_assocs  = frame.get<edm4eic::MCRecoTrackerHitAssociationCollection>("DRICHRawHitsAssociations");
+    auto isThrown = [](RVec<MCParticleData> parts){
+      return Filter(parts, [](auto p){ return p.generatorStatus==1; } );
+    };
+    //if(!isThrown) continue;
     if(!hit_assocs.isValid())
       throw std::runtime_error("cannot find hit associations");
 
     for(const auto& hit_assoc : hit_assocs) {
       for(const auto& sim_hit : hit_assoc.getSimHits()) {
-
         auto cellID = sim_hit.getCellID();
         auto pos    = sim_hit.getPosition();
+	h4->Fill(pos.x,pos.y);
+        auto mom    = sim_hit.getMomentum();
+        TVector3 p; p.SetX(mom.x); p.SetY(mom.y); p.SetZ(mom.z);
+	double Lambda = (1239.8/(p.Mag()*1.0e+9));
+        auto normZ  = drichGeo.GetSensorSurface(cellID);
+	double cosAng = (normZ.Unit()).Dot(p.Unit());
+	double angle = pi- acos(cosAng);
+	h2->Fill(Lambda,angle*(180/pi));
+	if (pos.x> 1180 && pos.x<1260){
+	  if(pos.y> -45 && pos.y<50){
+	    h1->Fill(angle*(180/pi));
+	  }
+	}
+	if(pos.x>1500 || pos.x<1000)
+	  h3->Fill(Lambda,angle*(180/pi));
         dd4hep::Position pos_global(pos.x*dd4hep::mm, pos.y*dd4hep::mm, pos.z*dd4hep::mm);
         auto pos_local = geo.GetSensorLocalPosition(cellID, pos_global);
         h->Fill(pos_local.y()/dd4hep::mm, pos_local.x()/dd4hep::mm);
@@ -92,13 +127,17 @@ int main(int argc, char** argv) {
   }
 
   gStyle->SetOptStat(0);
-  auto c = new TCanvas();
-  h->Draw();
+  auto c = new TCanvas(); c->Divide(2,2);
+  c->cd(1);h4->Draw();
+  c->cd(2);gPad->SetLogy();h1->Draw();
+  c->cd(3);h2->Draw("colz");
+  c->cd(4);h3->Draw("colz");
   fmt::print("NUMBER OF DIGITIZED PHOTONS: {}\n", h->GetEntries());
   if(interactiveOn) {
     fmt::print("\n\npress ^C to exit.\n\n");
     app->Run();
   } else {
     c->SaveAs("out/pixel_gaps.png");
+    c->SaveAs("out/pixel_gaps.root");
   }
 }