HCalCosmicCalib/src/HCalCalibTree.cc

0001 #include "HCalCalibTree.h"
0002
0003 #include <calobase/TowerInfoContainerv2.h>
0004 #include <calobase/TowerInfov2.h>
0005
0006 #include <fun4all/Fun4AllServer.h>
0007 #include <fun4all/Fun4AllHistoManager.h>
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phool/getClass.h>
0011
0012 #include <TFile.h>
0013 #include <TH1F.h>
0014
0015 #include <Event/Event.h>
0016 #include <Event/packet.h>
0017 #include <cassert>
0018 #include <sstream>
0019 #include <string>
0020
0021 using namespace std;
0022
0023 HCalCalibTree::HCalCalibTree(const std::string &name, const std::string &filename, const std::string &prefix)
0024     : SubsysReco(name), detector("HCALIN"), prefix(prefix), outfilename(filename)
0025 {
0026 }
0027
0028 HCalCalibTree::~HCalCalibTree() {
0029   delete hm;
0030 }
0031
0032 int HCalCalibTree::Init(PHCompositeNode *) {
0033   std::cout << std::endl << "HCalCalibTree::Init" << std::endl;
0034   hm = new Fun4AllHistoManager(Name());
0035   outfile = new TFile(outfilename.c_str(), "RECREATE");
0036
0037   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0038     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0039       std::string channel_histname = "h_channel_" + std::to_string(ieta) + "_" + std::to_string(iphi);
0040       h_channel_hist[ieta][iphi] = new TH1F(channel_histname.c_str(), "", 200, 0, 10000);
0041     }
0042   }
0043   h_waveformchi2 = new TH2F("h_waveformchi2", "", 1000, 0, 10000, 1000, 0, 100000);
0044   h_waveformchi2->GetXaxis()->SetTitle("peak (ADC)");
0045   h_waveformchi2->GetYaxis()->SetTitle("chi2");
0046
0047   h_check = new TH1F("h_check", "", 1000, 0, 0.01);
0048   h_eventnumber_record = new TH1F("h_eventnumber_record", "", 2, 0, 2);
0049
0050   if (prefix == "TOWERINFO_SIM_") {
0051     tower_threshold = tower_threshold_sim;
0052     vert_threshold = vert_threshold_sim;
0053     veto_threshold = veto_threshold_sim;
0054   } else {
0055     tower_threshold = tower_threshold_data;
0056     vert_threshold = vert_threshold_data;
0057     veto_threshold = veto_threshold_data;
0058   }
0059   std::cout << "Offline cut applied: " << std::endl;
0060   std::cout << "tower_threshold = " << tower_threshold << std::endl;
0061   std::cout << "vert_threshold = " << vert_threshold << std::endl;
0062   std::cout << "veto_threshold = " << veto_threshold << std::endl;
0063
0064   Fun4AllServer *se = Fun4AllServer::instance();
0065   se -> registerHistoManager(hm);
0066
0067   event = 0;
0068   goodevent = 0;
0069   return 0;
0070 }
0071
0072 int HCalCalibTree::process_event(PHCompositeNode *topNode) {
0073   if (event % 100 == 0) std::cout << "HCalCalibTree::process_event " << event << std::endl;
0074   goodevent_check = 0;
0075   process_towers(topNode);
0076   event++;
0077   if (goodevent_check == 1) goodevent++;
0078   return Fun4AllReturnCodes::EVENT_OK;
0079 }
0080
0081 int HCalCalibTree::process_towers(PHCompositeNode *topNode) {
0082   ostringstream nodenamev2;
0083   nodenamev2.str("");
0084   nodenamev2 << prefix << detector;
0085
0086   TowerInfoContainer *towers = findNode::getClass<TowerInfoContainer>(topNode, nodenamev2.str());
0087   if (!towers ) {
0088      std::cout << std::endl << "Didn't find node " << nodenamev2.str() << std::endl;
0089      return Fun4AllReturnCodes::EVENT_OK;
0090   }
0091
0092   int size = towers->size();
0093   for (int channel = 0; channel < size; channel++) {
0094     TowerInfo *tower = towers->get_tower_at_channel(channel);
0095     float energy = tower->get_energy();
0096     float chi2 = tower->get_chi2();
0097     unsigned int towerkey = towers->encode_key(channel);
0098     int ieta = towers->getTowerEtaBin(towerkey);
0099     int iphi = towers->getTowerPhiBin(towerkey);
0100     m_peak[ieta][iphi] = energy;
0101     m_chi2[ieta][iphi] = chi2;
0102     h_waveformchi2->Fill(m_peak[ieta][iphi], m_chi2[ieta][iphi]);
0103     h_check->Fill(m_peak[ieta][iphi]);
0104     if (m_chi2[ieta][iphi] > 10000) m_peak[ieta][iphi] = 0;
0105   }
0106
0107   // Apply cut
0108   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0109     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0110       if (m_peak[ieta][iphi] < tower_threshold) continue; //tower cut
0111       int up = iphi + 1;
0112       int down = iphi - 1;
0113       if (up > 63) up -= 64;
0114       if (down < 0) down += 64;
0115       if (m_peak[ieta][up] < vert_threshold || m_peak[ieta][down] < vert_threshold) continue; //vert cut
0116       if (ieta != 0 && (m_peak[ieta-1][up] > veto_threshold || m_peak[ieta-1][iphi] > veto_threshold || m_peak[ieta-1][down] > veto_threshold)) continue; // left veto cut
0117       if (ieta != 23 && (m_peak[ieta+1][up] > veto_threshold || m_peak[ieta+1][iphi] > veto_threshold || m_peak[ieta+1][down] > veto_threshold)) continue; //right veto cut
0118       h_channel_hist[ieta][iphi]->Fill(m_peak[ieta][iphi]);
0119       goodevent_check = 1;
0120     }
0121   }
0122   return Fun4AllReturnCodes::EVENT_OK;
0123 }
0124
0125 int HCalCalibTree::ResetEvent(PHCompositeNode *topNode) {
0126   return Fun4AllReturnCodes::EVENT_OK;
0127 }
0128
0129 int HCalCalibTree::End(PHCompositeNode * /*topNode*/) {
0130   std::cout << "HCalCalibTree::End" << std::endl;
0131   std::cout << "Number of events: " << event << std::endl;
0132   std::cout << "Number of good events: " << goodevent << std::endl;
0133   h_eventnumber_record->SetBinContent(1, event);
0134   h_eventnumber_record->SetBinContent(2, goodevent);
0135   outfile->cd();
0136   for (int ieta = 0; ieta < n_etabin; ++ieta) {
0137     for (int iphi = 0; iphi < n_phibin; ++iphi) {
0138       h_channel_hist[ieta][iphi]->Write();
0139       delete h_channel_hist[ieta][iphi];
0140     }
0141   }
0142   h_waveformchi2->Write();
0143   h_check->Write();
0144   h_eventnumber_record->Write();
0145   outfile->Close();
0146   delete outfile;
0147   hm->dumpHistos(outfilename, "UPDATE");
0148   return 0;
0149 }