sPHENIX_sims/Calib/Calib.cc

0001 #include "Calib.h"
0002 #include <fun4all/Fun4AllReturnCodes.h>
0003 #include <fun4all/PHTFileServer.h>
0004 #include <fun4all/Fun4AllServer.h>
0005
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4Hit.h>
0008 #include <g4main/PHG4TruthInfoContainer.h>
0009 #include <g4main/PHG4VtxPoint.h>
0010 #include <g4main/PHG4Particle.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>
0014 #include <phool/PHNode.h>
0015 #include <phool/PHNodeIterator.h>
0016 #include <phool/PHObject.h>
0017 #include <phool/getClass.h>
0018 #include <phool/phool.h>
0019
0020
0021 #include <calobase/RawTowerDefs.h>
0022 #include <calobase/RawTowerContainer.h>
0023 #include <calobase/RawTower.h>
0024 #include <calobase/RawTowerGeomContainer.h>
0025 #include <calobase/RawTowerGeom.h>
0026 #include <calobase/RawClusterContainer.h>
0027 #include <calobase/RawCluster.h>
0028 #include <calobase/RawTowerv1.h>
0029
0030 #include <phgeom/PHGeomUtility.h>
0031 #include <TTree.h>
0032 #include <TH2D.h>
0033 #include <TProfile2D.h>
0034 #include <TVector3.h>
0035 #include <TRandom3.h>
0036 #include <TMath.h>
0037
0038 #include <fstream>
0039
0040
0041
0042 using namespace std;
0043
0044 //----------------------------------------------------------------------------//
0045 //-- Constructor:
0046 //--  simple initialization
0047 //----------------------------------------------------------------------------//
0048
0049 Calib::Calib(const string &name) :
0050   SubsysReco(name), _event_tree(NULL)
0051 {
0052     //initialize
0053     _event = 0;
0054     _outfile_name = "Tower_test.root";
0055
0056 }
0057
0058 //----------------------------------------------------------------------------//
0059 //-- Init():
0060 //--   Intialize all histograms, trees, and ntuples
0061 //----------------------------------------------------------------------------//
0062 int Calib::Init(PHCompositeNode *topNode) {
0063
0064     cout << PHWHERE << " Opening file " << _outfile_name << endl;
0065
0066     PHTFileServer::get().open(_outfile_name, "RECREATE");
0067     PHTFileServer::get().cd(_outfile_name);
0068
0069     Double_t xbinlimits[25] =
0070      {
0071          -1.10001,-1.00833,-0.916665,-0.824998,-0.733332,-0.641665,-0.549998,-0.458332,-0.366665,-0.274998,-0.183332,-0.091665,0.,0.0916683,0.183335,0.275002,0.366668,0.458335,0.550002,0.641668,0.733335,0.825002,0.916668,1.00833,1.10001
0072    };
0073
0074    Double_t ybinlimits[65] =
0075     {
0076         -3.16565,-3.01936,-2.92119,-2.82301,-2.72484,-2.62666,-2.52849,-2.43032,-2.33214,-2.23397,-2.13579,-2.03762,-1.93944,-1.84127,-1.74309,-1.64492,-1.54674,-1.44857,-1.35039,-1.25222,-1.15404,-1.05587,-0.95769,-0.859519,-0.761344,-0.663169,-0.564994,-0.46682,-0.368645,-0.27047,-0.172295,-0.074121,0.024068900,0.122229,0.220404,0.318578,0.416753,0.514928,0.613103,0.711278,0.809452,0.907627,1.0058,1.10398,1.20215,1.30033,1.3985,1.49668,1.59485,1.69303,1.7912,1.88937,1.98755,2.08572,2.1839,2.28207,2.38025,2.47842,2.5766,2.67477,2.77295,2.87112,2.9693,3.06747,3.16565
0077     };
0078
0079     Long_t xbinnumb = sizeof(xbinlimits)/sizeof(Double_t) - 1;
0080     Long_t ybinnumb = sizeof(ybinlimits)/sizeof(Double_t) - 1;
0081
0082
0083     hprof2d  = new TProfile2D("hprof2d","Profile of tower energy versus eta and phi", xbinnumb,xbinlimits, ybinnumb,ybinlimits);
0084     hprof2d->SetXTitle("#eta");
0085     hprof2d->SetYTitle("#phi");
0086
0087     _event_tree = new TTree("event", "Calib => event info");
0088
0089     _event_tree->Branch("event", &_event, "_event/I");
0090     _event_tree->Branch("oCalib_ohcal", &_oCalib_ohcal, "_oCalib_ohcal/F");
0091
0092
0093     return Fun4AllReturnCodes::EVENT_OK;
0094 }
0095
0096 int Calib::InitRun(PHCompositeNode *topNode)
0097 {
0098     return Fun4AllReturnCodes::EVENT_OK;
0099 }
0100
0101 //----------------------------------------------------------------------------//
0102 //-- process_event():
0103 //--   Call user instructions for every event.
0104 //--   This function contains the analysis structure.
0105 //----------------------------------------------------------------------------//
0106
0107 int Calib::process_event(PHCompositeNode *topNode)
0108 {
0109     _event++;
0110
0111     GetNodes(topNode);
0112
0113     fill_tree(topNode);
0114
0115     return Fun4AllReturnCodes::EVENT_OK;
0116 }
0117
0118 //----------------------------------------------------------------------------//
0119 //-- End():
0120 //--   End method, wrap everything up
0121 //----------------------------------------------------------------------------//
0122
0123 int Calib::EndRun(PHCompositeNode *topNode)
0124 {
0125
0126     return Fun4AllReturnCodes::EVENT_OK;
0127
0128 }
0129
0130 int Calib::End(PHCompositeNode *topNode)
0131 {
0132
0133     PHTFileServer::get().cd(_outfile_name);
0134     PHTFileServer::get().write(_outfile_name);
0135
0136     return Fun4AllReturnCodes::EVENT_OK;
0137 }
0138
0139 //----------------------------------------------------------------------------//
0140 //-- fill_tree():
0141 //--   Fill the various trees...
0142 //----------------------------------------------------------------------------//
0143
0144 void Calib::fill_tree(PHCompositeNode *topNode)
0145 {
0146
0147   cout << _event << endl;
0148
0149
0150
0151     //ohcal  towers
0152        int ietabin, iphibin = -1;
0153       float oh_ocal = 0.0;
0154       float OHCAL_o[24][64] = {{0.0}};
0155       float OHCAL_tower_eta[24] = {0.0};
0156       float OHCAL_tower_phi[64] = {0.0};
0157
0158     //loop over ohcal towers
0159     RawTowerContainer::ConstRange begin_end_oo = _ohcal_towers_o->getTowers();
0160       RawTowerContainer::ConstIterator oh_o = begin_end_oo.first;
0161       for (; oh_o != begin_end_oo.second; ++oh_o) {
0162         RawTowerDefs::keytype towerid_o_o = oh_o->first;
0163         RawTower *rawtower_o_o = _ohcal_towers_o->getTower(towerid_o_o);
0164         if(rawtower_o_o) {
0165         RawTowerGeom *tgeo_o_o = _ohcal_towergeom->get_tower_geometry(towerid_o_o);
0166         // get towers energy
0167         oh_ocal  = rawtower_o_o->get_energy();
0168
0169         // binning in eta and phi
0170         ietabin = tgeo_o_o->get_column();
0171         iphibin = tgeo_o_o->get_row();
0172         if((iphibin >= 0)&&(ietabin >= 0))
0173         {
0174             OHCAL_o[ietabin][iphibin] += oh_ocal;
0175             OHCAL_tower_eta[ietabin] = tgeo_o_o->get_eta();
0176             OHCAL_tower_phi[iphibin] = tgeo_o_o->get_phi();
0177         }
0178
0179        }
0180       }
0181
0182
0183
0184     //---- Store energy ------//
0185     _oCalib_ohcal = 0.0;
0186
0187     for(int i=0; i<24; i++){
0188           for(int j=0; j<64; j++){
0189               _oCalib_ohcal += OHCAL_o[i][j];
0190               hprof2d->Fill(OHCAL_tower_eta[i], OHCAL_tower_phi[j],OHCAL_o[i][j]);
0191           }
0192     }
0193
0194
0195 _event_tree->Fill();
0196
0197
0198 return;
0199
0200 }
0201
0202 //----------------------------------------------------------------------------//
0203 //-- GetNodes():
0204 //--   Get all the all the required nodes off the node tree
0205 //----------------------------------------------------------------------------//
0206 int Calib::GetNodes(PHCompositeNode * topNode) {
0207
0208     //oHCAL
0209
0210     _ohcal_towers_o = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
0211     if (!_ohcal_towers_o)
0212       {
0213     std::cout << PHWHERE << ": Could not find node TOWER_CALIB_HCALOUT" << std::endl;
0214       return Fun4AllReturnCodes::ABORTEVENT;
0215       }
0216
0217
0218       //towergeom
0219     _ohcal_towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0220     if (! _ohcal_towergeom)
0221       {
0222     cout << PHWHERE << ": Could not find node TOWERGEOM_HCALOUT" << endl;
0223       return Fun4AllReturnCodes::ABORTEVENT;
0224       }
0225
0226      return Fun4AllReturnCodes::EVENT_OK;
0227 }
0228