plugins/hpDIRCrawHits/hpDIRCrawHitsProcessor.cc

0001
0002 //
0003 // Template for this file generated with eicmkplugin.py
0004 // Author: Nilanga Wickramaarachchi
0005 //
0006
0007 #include "hpDIRCrawHitsProcessor.h"
0008 #include "services/rootfile/RootFile_service.h"
0009
0010 // Include appropriate class headers. e.g.
0011 #include <edm4hep/SimCalorimeterHitCollection.h>
0012
0013 // The following just makes this a JANA plugin
0014 extern "C" {
0015     void InitPlugin(JApplication *app) {
0016         InitJANAPlugin(app);
0017         app->Add(new hpDIRCrawHitsProcessor);
0018     }
0019 }
0020
0021 //-------------------------------------------
0022 // InitWithGlobalRootLock
0023 //-------------------------------------------
0024 void hpDIRCrawHitsProcessor::InitWithGlobalRootLock(){
0025     auto rootfile_svc = GetApplication()->GetService<RootFile_service>();
0026     auto rootfile = rootfile_svc->GetHistFile();
0027     rootfile->mkdir("hpDIRCrawHits")->cd();
0028
0029     // Create histograms here. e.g.
0030     // hEraw  = new TH1D("Eraw",  "BEMC hit energy (raw)",  100, 0, 0.075);
0031     // hEdigi = new TH2D("Edigi", "BEMC hit energy (digi) vs. raw", 200, 0, 2000.0, 100, 0, 0.075);
0032     hHitposition = new TH2D("Hitposition", "hit position Y vs. hit position X; hit position x (mm); hit position y (mm)", 110, 683.5, 1047.875, 110, -175, 189.375);
0033     hHitposition_after = new TH2D("Hitposition_after", "hit position Y vs. hit position X; hit position x (mm); hit position y (mm)", 110, 683.5, 1047.875, 110, -175, 189.375);
0034
0035
0036     hnhits = new TH1D("nhits","; no. of hits per event", 500, 0, 1000);
0037     hnhits_qe = new TH1D("nhits_qe","; no. of hits per event", 500, 0, 1000);
0038     hnhits_qe_primary = new TH1D("nhits_qe_primary","; no. of hits per event", 500, 0, 1000);
0039
0040     hist_num_sim_hits = new TH1I("num_sim_hits","; no. of sim hits for raw hit", 10, 0, 10);
0041     hist_parent_pdg = new TH1I("parent_pdg","; parent PDG", 10000, -5000, 5000);
0042
0043     dirctree = new TTree("dirctree", "Tree for hpDIRC");
0044
0045     /// Event level
0046     dirctree->Branch("nhits", &nhits, "nhits/I");
0047
0048     /// Hit level
0049     dirctree->Branch("hit_pos_x", hit_pos_x, "hit_pos_x[nhits]/D");
0050     dirctree->Branch("hit_pos_y", hit_pos_y, "hit_pos_y[nhits]/D");
0051     dirctree->Branch("hit_pos_z", hit_pos_z, "hit_pos_z[nhits]/D");
0052     dirctree->Branch("mcp_id", mcp_id, "mcp_id[nhits]/I");
0053     dirctree->Branch("pixel_id", pixel_id, "pixel_id[nhits]/I");
0054     dirctree->Branch("hit_time", hit_time, "hit_time[nhits]/D");
0055 }
0056
0057 //-------------------------------------------
0058 // ProcessSequential
0059 //-------------------------------------------
0060 void hpDIRCrawHitsProcessor::ProcessSequential(const std::shared_ptr<const JEvent>& event) {
0061
0062     // Fill histograms here. e.g.
0063     // for (auto hit : rawhits) hEraw->Fill(hit.getEnergy());
0064   Int_t count = 0;
0065   Int_t count_qe = 0;
0066   Int_t count_qe_primary = 0;
0067
0068   for(auto hit : rawhits())
0069     {
0070       count_qe = count_qe+1;
0071
0072       hit_time[count] = hit->getTimeStamp() * timeResolution;
0073       edm4hep::Vector3d hit_position;
0074       int parent_pdg;
0075       bool is_from_primary_track;
0076
0077       for(const auto hit_assoc : in_hit_assocs()) {
0078             if(hit_assoc->getRawHit().isAvailable()) {
0079               if(hit_assoc->getRawHit().id() == hit->id()) {
0080                 if(hit_assoc->simHits_size() > 0) {
0081                   hit_position = hit_assoc->getSimHits(0).getPosition();
0082           hist_num_sim_hits->Fill(hit_assoc->simHits_size());
0083
0084           parent_pdg = hit_assoc->getSimHits(0).getMCParticle().getParents(0).getPDG();
0085           is_from_primary_track = hit_assoc->getSimHits(0).getMCParticle().getParents(0).getGeneratorStatus();
0086         }
0087           }
0088         }
0089       }
0090
0091       if(!is_from_primary_track) continue;
0092       count_qe_primary = count_qe_primary + 1;
0093
0094       hHitposition->Fill(hit_position[0], hit_position[1]);
0095       hist_parent_pdg->Fill(parent_pdg);
0096
0097       hit_pos_x[count] = hit_position[0];
0098       hit_pos_y[count] = hit_position[1];
0099       hit_pos_z[count] = hit_position[2];
0100
0101       if(hit_pos_x[count] > (prism_x_edge + prism_x_dim)) continue;
0102
0103       int mx = (hit_pos_x[count] - (prism_x_edge + gap_x)) / (MCP_total_dim + gap_x);
0104       int my = ((prism_y_edge - gap_y) - hit_pos_y[count]) / (MCP_total_dim + gap_y);
0105       mcp_id[count] = (6 * mx) + my;
0106
0107       float hit_mcp_x_edge = mcp_active_x_edge + (MCP_active_dim + gap_x + 4)*mx;
0108       float hit_mcp_y_edge = mcp_active_y_edge - (MCP_active_dim + gap_y + 4)*my;
0109
0110       if((hit_pos_x[count] - hit_mcp_x_edge) < 0) continue;
0111       if((hit_pos_y[count] - hit_mcp_y_edge) > 0) continue;
0112
0113       int pixel_x = (hit_pos_x[count] - hit_mcp_x_edge) / pixel_dim;
0114       int pixel_y = (hit_mcp_y_edge - hit_pos_y[count]) / pixel_dim;
0115
0116       if(pixel_x > 15 || pixel_y > 15) continue;
0117
0118       if(((hit_pos_x[count] - hit_mcp_x_edge) > 0) && (fmod(hit_pos_x[count] - hit_mcp_x_edge, pixel_dim) == 0)) pixel_x -= 1;
0119       if(((hit_mcp_y_edge - hit_pos_y[count]) > 0) && (fmod(hit_mcp_y_edge - hit_pos_y[count], pixel_dim) == 0)) pixel_y -= 1;
0120
0121       pixel_id[count] = (16 * pixel_x) + pixel_y;
0122
0123       hHitposition_after->Fill(hit_position[0], hit_position[1]);
0124       count = count+1;
0125     }
0126
0127   hnhits->Fill(count);
0128   hnhits_qe->Fill(count_qe);
0129   hnhits_qe_primary->Fill(count_qe_primary);
0130
0131   nhits = count;
0132   dirctree->Fill();
0133 }
0134
0135 //-------------------------------------------
0136 // FinishWithGlobalRootLock
0137 //-------------------------------------------
0138 void hpDIRCrawHitsProcessor::FinishWithGlobalRootLock() {
0139
0140     // Do any final calculations here.
0141 }
0142