LuminosityCounterGoodRuns/src/triggercountmodule.cc

0001 #include "triggercountmodule.h"
0002 #include <globalvertex/GlobalVertexMap.h>
0003 #include <globalvertex/GlobalVertexMapv1.h>
0004 #include <globalvertex/GlobalVertex.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/getClass.h>
0008 #include <ffaobjects/EventHeaderv1.h>
0009 #include <ffarawobjects/Gl1Packetv2.h>
0010 using namespace std;
0011
0012
0013 //____________________________________________________________________________..
0014 triggercountmodule::triggercountmodule(const std::string &filename, int rn, int segn, int maxseg, int debug, const std::string &name):
0015   SubsysReco(name)//).c_str())
0016 {
0017   _rn = rn;
0018   _filename = filename;
0019   _name = name;
0020
0021   _tree = new TTree("_tree","a tree to count events");
0022   _debug = debug;
0023   _nseg = segn;
0024   _lastseg = maxseg;
0025 }
0026
0027 //____________________________________________________________________________..
0028 triggercountmodule::~triggercountmodule()
0029 {
0030
0031 }
0032
0033 //____________________________________________________________________________..
0034 int triggercountmodule::Init(PHCompositeNode *topNode)
0035 {
0036
0037   if(_debug > 1) cout << "Begin init: " << endl;
0038
0039   _tree->Branch("badFlag",&_badFlag,"badFlag/I");
0040   _tree->Branch("startLive",_startLive,"startLive[64]/l");
0041   _tree->Branch("endLive",_endLive,"endLive[64]/l");
0042   _tree->Branch("startScal",_startScal,"startScal[64]/l");
0043   _tree->Branch("endScal",_endScal,"endScal[64]/l");
0044   _tree->Branch("nevt",&_evtn,"nevt/l");
0045   _tree->Branch("avgPS",_avgPS,"avgPS[64]/D");
0046   _tree->Branch("trigCounts",_trigCounts,"trigCounts[6][64]/l");
0047   _tree->Branch("eMBDlive",_eMBDlive,"eMBDlive[5]/D");
0048   _tree->Branch("startBCO",&_startBCO,"startBCO/l");
0049   _tree->Branch("endBCO",&_endBCO,"endBCO/l");
0050   _tree->Branch("nBunch",&_nBunch,"nBunch/I");
0051   _tree->Branch("startRaw",_startRaw,"startRaw[64]/l");
0052   _tree->Branch("endRaw",_endRaw,"endRaw[64]/l");
0053   _mbzhist = new TH1D("mbzhist","",600,-300,300);
0054   _bunchHist = new TH1D("bunchHist","",121,0,121);
0055   if(_debug > 0)
0056     {
0057       cout << "Printing tree: " << endl;
0058       _tree->Print();
0059     }
0060   return Fun4AllReturnCodes::EVENT_OK;
0061 }
0062
0063 //____________________________________________________________________________..
0064 int triggercountmodule::InitRun(PHCompositeNode *topNode)
0065 {
0066   if(_debug > 1) cout << "Initializing!" << endl;
0067   return Fun4AllReturnCodes::EVENT_OK;
0068 }
0069
0070 //____________________________________________________________________________..
0071
0072
0073 int triggercountmodule::process_event(PHCompositeNode *topNode)
0074 {
0075
0076   if(_debug > 1) cout << endl << endl << endl << "triggercountmodule: Beginning event processing" << endl;
0077
0078     Gl1Packetv2* gl1 = findNode::getClass<Gl1Packetv2>(topNode, "GL1Packet");
0079   if(!gl1)
0080     {
0081       cout << "No trigger info!" << endl;
0082       return Fun4AllReturnCodes::ABORTRUN;
0083     }
0084
0085   _bunchHist->Fill(gl1->lValue(0,"BunchNumber"));
0086
0087   if(_evtn == 0)
0088     {
0089       _startBCO = gl1->lValue(0,"BCO");
0090       for(int i=0; i<64; ++i)
0091     {
0092       _startLive[i] = gl1->lValue(i,1);
0093       if(_debug > 2) cout << gl1->lValue(i,1) << endl;
0094       _startScal[i] = gl1->lValue(i,2);
0095       _startRaw[i] = gl1->lValue(i,0);
0096     }
0097     }
0098   //cout << _evtn << endl;
0099   ++_evtn;
0100   _endBCO = gl1->lValue(0,"BCO");
0101   for(int i=0; i<64; ++i)
0102     {
0103       _endLive[i] = gl1->lValue(i,1);
0104       if(_debug > 2) cout << "end" << gl1->lValue(i,1) << endl;
0105       _endScal[i] = gl1->lValue(i,2);
0106       _endRaw[i] = gl1->lValue(i,0);
0107     }
0108
0109   if(_debug > 1) cout << "Getting gl1 trigger vector from: " << gl1 << endl;
0110   long long unsigned int _trigvec = gl1->getScaledVector();
0111
0112   int ismbtrig = (_trigvec >> 10) & 1;
0113
0114   GlobalVertexMap* gvtxmap = findNode::getClass<GlobalVertexMapv1>(topNode, "GlobalVertexMap");
0115
0116   float zvtx = NAN;
0117
0118   if(gvtxmap)
0119     {
0120       GlobalVertex* gvtx = gvtxmap->begin()->second;
0121       if(gvtx)
0122         {
0123           auto startIter = gvtx->find_vertexes(_vtxtype);
0124           auto endIter = gvtx->end_vertexes();
0125           for(auto iter = startIter; iter != endIter; ++iter)
0126             {
0127               const auto &[type, vertexVec] = *iter;
0128               if(type != _vtxtype) continue;
0129               for(const auto *vertex : vertexVec)
0130                 {
0131                   if(!vertex) continue;
0132                   zvtx = vertex->get_z();
0133                 }
0134             }
0135         }
0136     }
0137
0138   int zthresh[5] = {30,60,200,10,1000};
0139
0140
0141
0142   for(int i=0; i<5; ++i)
0143     {
0144       if(abs(zvtx) < zthresh[i] && !std::isnan(zvtx))
0145     {
0146       for(int j=0; j<64; ++j)
0147         {
0148           if((_trigvec >> j) & 1) _trigCounts[i][j]++;
0149         }
0150     }
0151     }
0152
0153   for(int j=0; j<64; ++j)
0154     {
0155       if((_trigvec >> j) & 1) _trigCounts[5][j]++;
0156     }
0157
0158   if(!ismbtrig)
0159     {
0160       if(_debug > 2) cout << "non MB event, not counting" << endl;
0161       return Fun4AllReturnCodes::ABORTEVENT;
0162     }
0163
0164   if(std::isnan(zvtx))
0165     {
0166       if(_debug > 0) cout << "zvtx is NAN after attempting to grab it. ABORT EVENT!" << endl;
0167       return Fun4AllReturnCodes::ABORTEVENT;
0168     }
0169   if(abs(zvtx) > 300)
0170     {
0171       if(_debug > 2) cout << "zvtx > 300. zvtx: " << zvtx << endl;
0172       return Fun4AllReturnCodes::ABORTEVENT;
0173     }
0174
0175   _mbzhist->Fill(zvtx);
0176
0177   return Fun4AllReturnCodes::EVENT_OK;
0178
0179 }
0180 //____________________________________________________________________________..
0181 int triggercountmodule::ResetEvent(PHCompositeNode *topNode)
0182 {
0183   if (Verbosity() > 0)
0184     {
0185       std::cout << "triggercountmodule::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0186     }
0187   return Fun4AllReturnCodes::EVENT_OK;
0188 }
0189
0190 //____________________________________________________________________________..
0191 int triggercountmodule::EndRun(const int runnumber)
0192 {
0193   if (Verbosity() > 0)
0194     {
0195       std::cout << "triggercountmodule::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0196     }
0197   return Fun4AllReturnCodes::EVENT_OK;
0198 }
0199
0200 //____________________________________________________________________________..
0201 int triggercountmodule::End(PHCompositeNode *topNode)
0202 {
0203   if (Verbosity() > 0)
0204     {
0205       std::cout << "triggercountmodule::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0206     }
0207   //cout << _evtn << endl;
0208   if(_nseg != _lastseg && _evtn < 49500)
0209     {
0210       _badFlag = 1;
0211     }
0212
0213   for(int i=0; i<64; ++i)
0214     {
0215       if(_endScal[i] == _startScal[i]) _avgPS[i] = 0;
0216       else _avgPS[i] = ((float)(_endLive[i]-_startLive[i]))/(_endScal[i]-_startScal[i]);
0217       if(_debug > 2) cout << "realend" << _endLive[i] << endl;
0218     }
0219   for(int i=0; i<5; ++i)
0220     {
0221       _eMBDlive[i] = _avgPS[10]*_trigCounts[i][10];
0222     }
0223
0224   int avgBunchCounts = _bunchHist->Integral() / 121.0;
0225
0226   for(int i=1; i<122; ++i)
0227     {
0228       if(_bunchHist->GetBinContent(i) > avgBunchCounts/5)
0229     {
0230       _nBunch++;
0231     }
0232     }
0233   _outfile = new TFile(_filename.c_str(), "RECREATE");
0234   _outfile->cd();
0235   _bunchHist->Write();
0236   _mbzhist->Write();
0237   if(_debug > 0) _tree->Print();
0238   _tree->Fill();
0239   _tree->Write();
0240   _outfile->Write();
0241   _outfile->Close();
0242
0243
0244   return Fun4AllReturnCodes::EVENT_OK;
0245 }
0246
0247 //____________________________________________________________________________..
0248 int triggercountmodule::Reset(PHCompositeNode *topNode)
0249 {
0250   if (Verbosity() > 0)
0251     {
0252       std::cout << "triggercountmodule::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0253     }
0254   return Fun4AllReturnCodes::EVENT_OK;
0255 }
0256
0257 //____________________________________________________________________________..
0258 void triggercountmodule::Print(const std::string &what) const
0259 {
0260   std::cout << "triggercountmodule::Print(const std::string &what) const Printing info for " << what << std::endl;
0261 }