Inheritance diagram for evgen::ldm::MeVPrtlTestRayTrace:

Public Member Functions
	MeVPrtlTestRayTrace (fhicl::ParameterSet const &p)

	MeVPrtlTestRayTrace (MeVPrtlTestRayTrace const &)=delete

	MeVPrtlTestRayTrace (MeVPrtlTestRayTrace &&)=delete

MeVPrtlTestRayTrace &	operator= (MeVPrtlTestRayTrace const &)=delete

MeVPrtlTestRayTrace &	operator= (MeVPrtlTestRayTrace &&)=delete

void	analyze (const art::Event &e) override

	~MeVPrtlTestRayTrace () noexcept

Private Attributes
std::unique_ptr < evgen::ldm::IMesonGen >	fGenTool

std::unique_ptr < evgen::ldm::IMeVPrtlFlux >	fFluxTool

std::vector< std::unique_ptr < evgen::ldm::IRayTrace > >	fRayTools

unsigned	fNCall

unsigned	fNEvt

TTree *	fTree

std::vector< double >	fBranchWeights

MeVPrtlTruth *	fMeVPrtl

Detailed Description

Definition at line 60 of file MeVPrtlTestRayTrace_module.cc.

Constructor & Destructor Documentation

evgen::ldm::MeVPrtlTestRayTrace::MeVPrtlTestRayTrace ( fhicl::ParameterSet const & p )

explicit

Definition at line 90 of file MeVPrtlTestRayTrace_module.cc.

   : EDAnalyzer{p}
 {
   // bring in the tools
   fGenTool = art::make_tool<IMesonGen>(p.get<fhicl::ParameterSet>("MesonGen"));
   fFluxTool = art::make_tool<IMeVPrtlFlux>(p.get<fhicl::ParameterSet>("Flux"));
 
   for (auto const &rayconfig: p.get<std::vector<fhicl::ParameterSet>>("RayTraces")) {
     fRayTools.push_back(art::make_tool<IRayTrace>(rayconfig));
     fBranchWeights.push_back(0);
   }
 
   fNCall = p.get<unsigned>("NCall", 10000);
   fNEvt = 0;
 
   art::ServiceHandle<art::TFileService> tfs;
   fTree = tfs->make<TTree>("testraytrace", "testraytrace");
   fMeVPrtl = new evgen::ldm::MeVPrtlTruth();
   fTree->Branch("mevprtl", &fMeVPrtl);
   for (unsigned iray = 0; iray < fRayTools.size(); iray++) {
     fTree->Branch((std::string(fRayTools[iray]->Name()) + "_wgt").c_str(), &fBranchWeights[iray]);
   }
 
 }

evgen::ldm::MeVPrtlTestRayTrace::MeVPrtlTestRayTrace ( MeVPrtlTestRayTrace const & )

delete

evgen::ldm::MeVPrtlTestRayTrace::MeVPrtlTestRayTrace ( MeVPrtlTestRayTrace && )

delete

evgen::ldm::MeVPrtlTestRayTrace::~MeVPrtlTestRayTrace ( )

inlinenoexcept

Definition at line 75 of file MeVPrtlTestRayTrace_module.cc.

75 {}

Member Function Documentation

void evgen::ldm::MeVPrtlTestRayTrace::analyze ( const art::Event & e )

override

Definition at line 115 of file MeVPrtlTestRayTrace_module.cc.

 {
   // get the next MeVPrtl Truth 
   while (1) {
     simb::MCFlux kaon = fGenTool->GetNext();
 
     evgen::ldm::KaonParent kaonp(kaon);
     bool is_kaon = kaonp.kaon_pdg != 0;
 
     // (void) is_kaon;
     if (is_kaon) { 
      std::cout << "Flux is kaon (" << is_kaon << "). Weight: " << kaonp.weight << ". Produced with energy: " << kaonp.mom.E() 
              << " M=" << kaonp.mom.M() << " P=(" << kaonp.mom.Px() << ", " << kaonp.mom.Py() << ", " << kaonp.mom.Pz() << ") At: ("
              << kaonp.pos.X() << ", " << kaonp.pos.Y() << ", " << kaonp.pos.Z() << ")" << std::endl;
     }
 
     bool success;
 
     evgen::ldm::MeVPrtlFlux flux;
     double flux_weight;
 
     success = fFluxTool->MakeFlux(kaon, flux, flux_weight);
     if (!success) continue;
 
     std::cout << "New flux. E=" << flux.mom.E() << " At: (" << flux.pos.X() << ", " << flux.pos.Y() << ", " << flux.pos.Z() << ")" << std::endl;
     std::cout << "P=(" << flux.mom.Px() << ", " << flux.mom.Py() << ", " << flux.mom.Pz() << ")" << std::endl;
     std::cout << "Flux weight: " << flux_weight << std::endl;
 
 
     // See if an intersection is possible
     double costh_crit = minKinematicCosTheta(flux.kmom.M(), flux.sec.M(), flux.mom.M(), flux.kmom.E());
     TVector3 det(0., 0., 0.); // detector should be near origin
     double costh = flux.kmom.Vect().Unit().Dot((det - flux.pos.Vect()).Unit());
     std::cout << "COSTH CRIT: " << costh_crit << " DETECTOR COSTH: " << costh << std::endl;
     if (costh < costh_crit) continue;
 
     std::cout << "CALLING RAY TOOLS\n";
     fNEvt ++;
     unsigned iray = 0;
     std::vector<double> weightsum(fRayTools.size(), 0.);
     std::vector<double> time(fRayTools.size(), 0.);
     for (auto const &r: fRayTools) {
       std::array<TVector3, 2> intersection;
       double ray_weight;
       std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
 
       for (unsigned icall = 0; icall < fNCall; icall++) {
         success = r->IntersectDetector(flux, intersection, ray_weight);
         if (!success) ray_weight = 0.;
         weightsum[iray] += ray_weight;
       }
 
       std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
       std::chrono::duration<double, std::milli> duration = t2 - t1;
       time[iray] = duration.count();
 
       iray ++;
     }
 
     std::cout << "DONE CALLING RAY TOOLS\n";
     for (unsigned iray = 0; iray < fRayTools.size(); iray++) {
       std::cout << fRayTools[iray]->Name() << " " << (time[iray] / fNCall) << " " << (weightsum[iray] / fNCall) << std::endl;
     }
 
     evgen::ldm::MeVPrtlDecay decay;
     std::array<TVector3, 2> intersection;
     *fMeVPrtl = evgen::ldm::MeVPrtlTruth(flux, decay, intersection, flux_weight, 1., 1., 0.);
     for (unsigned iray = 0; iray < fRayTools.size(); iray++) {
       fBranchWeights[iray] = weightsum[iray] / fNCall;
     }
 
     fTree->Fill();
 
     break;
   }
 }