Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:14:33 

0001 #include <cassert>
0002 #include <iostream>
0003 
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/PHTFileServer.h>
0006 
0007 #include <phool/getClass.h>
0008 #include <phool/PHCompositeNode.h>
0009 
0010 //#include <calotrigger/TriggerAnalyzer.h>
0011 #include <ffarawobjects/Gl1Packet.h>
0012 #include <globalvertex/GlobalVertexMapv1.h>
0013 
0014 #include <jetbase/JetContainerv1.h>
0015 #include <calobase/TowerInfoContainer.h>
0016 #include <calobase/TowerInfo.h>
0017 
0018 #include <TSQLServer.h>
0019 #include <TSQLResult.h>
0020 #include <TSQLRow.h>
0021 
0022 #include <jetbase/FastJetAlgo.h>
0023 
0024 #include "TH1D.h"
0025 #include "TH3F.h"
0026 #include "TH2D.h"
0027 
0028 //module components
0029 #include "Ana_PPG09_Mod.h"
0030 
0031 
0032 //____________________________________________________________________________..
0033 Ana_PPG09_Mod::Ana_PPG09_Mod(const std::string &recojetname, const std::string& outputfilename):
0034  SubsysReco("PPG09_Mod_"+recojetname)
0035  , m_recoJetName(recojetname)
0036  , m_outputFileName(outputfilename)
0037 {
0038   if(Verbosity() > 5){
0039      std::cout << "Ana_PPG09_Mod::Ana_PPG09_Mod(const std::string &name) Calling ctor" << std::endl;
0040   }
0041 }
0042 
0043 //____________________________________________________________________________..
0044 Ana_PPG09_Mod::~Ana_PPG09_Mod()
0045 {
0046   if(Verbosity() > 5){
0047      std::cout << "Ana_PPG09_Mod::~Ana_PPG09_Mod() Calling dtor" << std::endl;
0048   }
0049 }
0050 
0051 //____________________________________________________________________________..
0052 int Ana_PPG09_Mod::Init(PHCompositeNode *topNode)
0053 {
0054   std::cout << "Ana_PPG09_Mod::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0055   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0056   std::cout << "Ana_PPG09_Mod::Init - Output to " << m_outputFileName << std::endl;
0057   std::cout << "Chosen Jet Cuts: Lead Pt >= " << Lead_RPt_Cut << " GeV, All Pt >= " << All_RPt_Cut << " GeV, |Z-Vertex| <= " << ZVtx_Cut << ", Num. of Constits. > " << NComp_Cut << std::endl;
0058 
0059   char hname[99];
0060 
0061 
0062   //Constructing Histograms
0063   for(int i = 0; i < nTrigs; i++){
0064      //Jet Plots
0065      sprintf(hname,"h_Eta_Phi_Pt_%d",i);
0066      h_Eta_Phi_Pt_[i] = new TH3F(hname,hname,22,-1.1,1.1,64,-3.2,3.2,99,1,100);
0067      
0068      //OHCal Tower Plots
0069      sprintf(hname,"h_oHCal_CS_Eta_Phi_E_%d",i);
0070      h_oHCal_CS_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0071      sprintf(hname,"h_oHCal_C_Eta_Phi_E_%d",i);
0072      h_oHCal_C_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0073      sprintf(hname,"h_oHCal_Raw_Eta_Phi_E_%d",i);
0074      h_oHCal_Raw_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0075      sprintf(hname,"h_oHCal_Jet_Eta_Phi_E_%d",i);
0076      h_oHCal_Jet_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0077      sprintf(hname,"h_oHCal_TE_Sub_Eta_Phi_E_%d",i);
0078      h_oHCal_TE_Sub_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,12000,-600,600);
0079 
0080      //IHCal Tower Plots
0081      sprintf(hname,"h_iHCal_CS_Eta_Phi_E_%d",i);
0082      h_iHCal_CS_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0083      sprintf(hname,"h_iHCal_C_Eta_Phi_E_%d",i);
0084      h_iHCal_C_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0085      sprintf(hname,"h_iHCal_Raw_Eta_Phi_E_%d",i);
0086      h_iHCal_Raw_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0087      sprintf(hname,"h_iHCal_Jet_Eta_Phi_E_%d",i);
0088      h_iHCal_Jet_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0089      sprintf(hname,"h_iHCal_TE_Sub_Eta_Phi_E_%d",i);
0090      h_iHCal_TE_Sub_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,12000,-600,600);
0091 
0092      //EMCal Tower/Retower Plots
0093      sprintf(hname,"h_EMCal_CS_Eta_Phi_E_%d",i);
0094      h_EMCal_CS_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0095      sprintf(hname,"h_EMCal_CR_Eta_Phi_E_%d",i);
0096      h_EMCal_CR_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0097      sprintf(hname,"h_EMCal_C_Eta_Phi_E_%d",i);
0098      h_EMCal_C_Eta_Phi_E_[i] = new TH3F(hname,hname,96,-0.5,95.5,256,-0.5,255.5,200,-100,100);
0099      sprintf(hname,"h_EMCal_Raw_Eta_Phi_E_%d",i);
0100      h_EMCal_Raw_Eta_Phi_E_[i] = new TH3F(hname,hname,96,-0.5,95.5,256,-0.5,255.5,200,-100,100);
0101      sprintf(hname,"h_EMCal_Jet_Eta_Phi_E_%d",i);
0102      h_EMCal_Jet_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,200,-100,100);
0103      sprintf(hname,"h_EMCal_TE_Sub_Eta_Phi_E_%d",i);
0104      h_EMCal_TE_Sub_Eta_Phi_E_[i] = new TH3F(hname,hname,24,-0.5,23.5,64,-0.5,63.5,12000,-600,600);
0105 
0106      //Event Plots
0107      sprintf(hname,"h_ZVtx_%d",i);
0108      h_ZVtx_[i] = new TH1D(hname,hname,200,-100,100);
0109   }
0110   
0111   h_nJetsAboveThresh = new TH2D("h_nJetsAboveThresh","h_nJetsAboveThresh",50,0.,50.,2000,4000.,6000.);
0112   h_EventCount = new TH1D("h_EventCount","Events",5,-2.5,2.5);
0113   h_theJetSpectrum = new TH1F("h_JetSpectrum","h_JetSpectrum",nJetBins, jetBins);
0114   
0115   //grab all trigger information
0116   initializePrescaleInformationFromDB(m_runnumber);
0117   
0118   std::cout << "Ana_PPG09_Mod::init - Initialization completed" << std::endl;
0119   std::cout << "scaleDowns vector size: "<< scaleDowns.size() << std::endl;
0120   for(int i = 0; i < (int)scaleDowns.size(); i++)
0121     {
0122       std::cout << "scaleDowns[" << i << "]: " << scaleDowns[i] << std::endl;
0123     }
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126 
0127 //____________________________________________________________________________..
0128 int Ana_PPG09_Mod::process_event(PHCompositeNode *topNode)
0129 {
0130   
0131   m_event++;
0132   if((m_event % 1000) == 0)  std::cout << "Ana_PPG09_Mod::process_event - Event number = " << m_event << std::endl;
0133 
0134   h_EventCount->Fill(-1);
0135 
0136   GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0137   if(!global_vtxmap || global_vtxmap->empty()){
0138      if(Verbosity() > 0){std::cout << "Aborted Event number = " << m_event << ", no global node" << std::endl;}  
0139      return Fun4AllReturnCodes::ABORTEVENT;
0140   }
0141 
0142   GlobalVertex *vtx = global_vtxmap->begin()->second;
0143   float vtx_z = NAN;
0144   if(!vtx){
0145      if(Verbosity() > 0){std::cout << "Aborted Event number = " << m_event << ", no z-vertex" << std::endl;}
0146      return Fun4AllReturnCodes::ABORTEVENT;
0147   }
0148   else{
0149      vtx_z = vtx->get_z();
0150      if(vtx_z < -ZVtx_Cut || vtx_z > ZVtx_Cut){
0151         if(Verbosity() > 0){std::cout << "Aborted Event number = " << m_event << ", z-vertex out of range" << std::endl;}
0152         return Fun4AllReturnCodes::ABORTEVENT;
0153      }
0154   }
0155 
0156   ///Recording Event Triggers
0157   //triggeranalyzer->decodeTriggers(topNode);
0158   std::vector<int> m_triggers;
0159   Gl1Packet *gl1Packet = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0160   if(gl1Packet){
0161      auto scaled_vector = gl1Packet->getScaledVector();
0162      for(int i = 0; i < 32; i++){
0163         if((scaled_vector & (int)std::pow(2,i)) != 0){
0164            m_triggers.push_back(i);
0165         }
0166      }
0167   }
0168   else{
0169      std::cout << "gl1Packet not found" << std::endl;
0170   }
0171   
0172 
0173   JetContainer* jets_Cont = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0174   if(!jets_Cont){std::cout << "Jets are Missing !" << std::endl;}
0175 
0176   //Checking 
0177   double Lead_Check = 0;
0178   int jetAboveThresh = 0;
0179   for(unsigned int ijet = 0; ijet < jets_Cont->size(); ++ijet){
0180      Jet* jet = jets_Cont->get_jet(ijet);
0181      if(jet->get_pt() > Lead_Check){Lead_Check = jet->get_pt();}
0182      if(jet->get_pt() > 8){jetAboveThresh++;}
0183   }
0184 
0185   h_nJetsAboveThresh -> Fill(m_runnumber,jetAboveThresh);
0186   if(Lead_Check < Lead_RPt_Cut){
0187      if(Verbosity() > 0){std::cout << "Aborted Event number = " << m_event << ", no jet above lead pt threshold" << std::endl;}
0188      return Fun4AllReturnCodes::ABORTEVENT;
0189   }
0190 
0191   int eventPrescale = getEventPrescale(Lead_Check, m_triggers);
0192   
0193   TowerInfoContainer* towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_CEMC");
0194   TowerInfoContainer* towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALIN");
0195   TowerInfoContainer* towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALOUT");
0196    if(!towersEM3 || !towersIH3 || !towersOH3){std::cout << "Raw Towers are Missing !" << std::endl;}
0197 
0198   /*
0199   TowerInfoContainer* towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0200   TowerInfoContainer* towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0201   TowerInfoContainer* towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT"); 
0202   */
0203 
0204   TowerInfoContainer* CtowersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0205   TowerInfoContainer* CtowersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0206   TowerInfoContainer* CtowersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0207   TowerInfoContainer* CRtowersEM3 = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_CEMC_RETOWER");
0208   if(!CtowersEM3 || !CtowersIH3 || !CtowersOH3 || !CRtowersEM3){std::cout << "Calib Towers are Missing !" << std::endl;}
0209 
0210   TowerInfoContainer* CStowersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0211   TowerInfoContainer* CStowersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0212   TowerInfoContainer* CStowersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0213   if(!CStowersEM3 || !CStowersIH3 || !CStowersOH3){std::cout << "Calib Sub. Towers are Missing !" << std::endl;}
0214 
0215   double Index_Plots[nTrigs] = {0};
0216 
0217   h_EventCount->Fill(1);
0218 
0219   
0220   
0221   for(unsigned int ijet = 0; ijet < jets_Cont->size(); ++ijet){
0222      Jet* jet = jets_Cont->get_jet(ijet);
0223 
0224      if(ijet == 0){
0225         for(int k = 0; k < (int)m_triggers.size(); k++){
0226            if(m_triggers.at(k) == 10){ 
0227               Index_Plots[0] = 1;
0228               h_ZVtx_[0]->Fill(vtx_z);
0229            }
0230            if(m_triggers.at(k) == 12){ 
0231               Index_Plots[1] = 1;
0232               h_ZVtx_[1]->Fill(vtx_z);
0233            }
0234            if(m_triggers.at(k) == 17){ 
0235               Index_Plots[2] = 1;
0236               h_ZVtx_[2]->Fill(vtx_z);
0237            }
0238            if(m_triggers.at(k) == 18){
0239               Index_Plots[3] = 1;
0240               h_ZVtx_[3]->Fill(vtx_z);
0241            } 
0242            if(m_triggers.at(k) == 19){ 
0243               Index_Plots[4] = 1;
0244               h_ZVtx_[4]->Fill(vtx_z);
0245            }
0246            if(m_triggers.at(k) == 21){ 
0247               Index_Plots[5] = 1;
0248               h_ZVtx_[5]->Fill(vtx_z);
0249            }
0250            if(m_triggers.at(k) == 22){ 
0251               Index_Plots[6] = 1;
0252               h_ZVtx_[6]->Fill(vtx_z);
0253        }
0254            if(m_triggers.at(k) == 23){
0255               Index_Plots[7] = 1;
0256               h_ZVtx_[7]->Fill(vtx_z);
0257            }
0258         }
0259 
0260         for(int k = 0; k < nTrigs; k++){
0261            if(Index_Plots[k] != 1) continue;
0262        //float scaleDown = scaleDowns[k];
0263        for(int Det = 0; Det < 3; Det++){
0264               TowerInfoContainer* towersCal;
0265               TowerInfoContainer* towersSCal;
0266               TowerInfoContainer* towersRaw;
0267               TowerInfoContainer* towersRCal;
0268           TH3F* RawTowers;
0269               TH3F* CalibTowers;
0270               TH3F* CalibSubTowers;
0271               TH3F* RCalibTowers;
0272           TH3F* SubETowers; 
0273 
0274               if(Det == 0){
0275                  towersCal = CtowersEM3;
0276          towersRCal = CRtowersEM3;
0277                  towersSCal = CStowersEM3;
0278                  towersRaw = towersEM3; 
0279                  RawTowers = h_EMCal_Raw_Eta_Phi_E_[k];
0280                  CalibTowers = h_EMCal_C_Eta_Phi_E_[k];
0281                  CalibSubTowers = h_EMCal_CS_Eta_Phi_E_[k];      
0282                  RCalibTowers = h_EMCal_CR_Eta_Phi_E_[k];
0283          SubETowers = h_EMCal_TE_Sub_Eta_Phi_E_[k];
0284           }          
0285               if(Det == 1){
0286                  towersCal = CtowersIH3;
0287                  towersSCal = CStowersIH3;
0288                  towersRaw = towersIH3;
0289                  RawTowers = h_iHCal_Raw_Eta_Phi_E_[k];
0290                  CalibTowers = h_iHCal_C_Eta_Phi_E_[k];
0291                  CalibSubTowers = h_iHCal_CS_Eta_Phi_E_[k];
0292                  SubETowers = h_iHCal_TE_Sub_Eta_Phi_E_[k];
0293           }
0294               if(Det == 2){
0295                  towersCal = CtowersOH3;
0296                  towersSCal = CStowersOH3;
0297                  towersRaw = towersOH3;
0298                  RawTowers = h_oHCal_Raw_Eta_Phi_E_[k];
0299                  CalibTowers = h_oHCal_C_Eta_Phi_E_[k];
0300                  CalibSubTowers = h_oHCal_CS_Eta_Phi_E_[k];
0301          SubETowers = h_oHCal_TE_Sub_Eta_Phi_E_[k];
0302               }              
0303  
0304               int Chan_Num = -1;
0305               if(Det == 0){Chan_Num = 24576;}
0306               else{Chan_Num = 1536;}
0307               
0308               for (int channel = 0; channel < Chan_Num; channel++){
0309                  //Raw Towers
0310                  float CTowerE = -99;
0311          float CSTowerE = -99;
0312               
0313          TowerInfo *towerRaw = towersRaw->get_tower_at_channel(channel);
0314                  unsigned int channelkeyR = towersRaw->encode_key(channel);
0315                  int ietaR = towersRaw->getTowerEtaBin(channelkeyR);
0316                  int iphiR = towersRaw->getTowerPhiBin(channelkeyR);
0317 
0318                  //Calib Towers
0319                  TowerInfo *towerCal = towersCal->get_tower_at_channel(channel);
0320                  unsigned int channelkeyC = towersCal->encode_key(channel);
0321                  int ietaC = towersCal->getTowerEtaBin(channelkeyC);
0322                  int iphiC = towersCal->getTowerPhiBin(channelkeyC);
0323                  if(!Det == 0){CTowerE = towerCal->get_energy();}
0324 
0325                  //Calib Subtracted Towers
0326          if(channel < 1536){
0327                     TowerInfo *towerSCal = towersSCal->get_tower_at_channel(channel);
0328                     unsigned int channelkeySC = towersSCal->encode_key(channel);
0329                     int ietaSC = towersSCal->getTowerEtaBin(channelkeySC);
0330                     int iphiSC = towersSCal->getTowerPhiBin(channelkeySC); 
0331                     CalibSubTowers->Fill(ietaSC,iphiSC,towerSCal->get_energy(),scaleDowns[k]);
0332                     CSTowerE = towerSCal->get_energy();
0333 
0334                     //Calib. Retowers
0335             if(Det == 0){
0336                    TowerInfo *towerRCal = towersRCal->get_tower_at_channel(channel);
0337                        unsigned int channelkeyRC = towersRCal->encode_key(channel);
0338                        int ietaRC = towersRCal->getTowerEtaBin(channelkeyRC);
0339                        int iphiRC = towersRCal->getTowerPhiBin(channelkeyRC);
0340                RCalibTowers->Fill(ietaRC,iphiRC,towerRCal->get_energy(),scaleDowns[k]); 
0341                CTowerE = towerRCal->get_energy();
0342             }
0343          }
0344     
0345              double SubE = CTowerE - CSTowerE;   
0346          SubETowers->Fill(ietaC,iphiC,SubE,scaleDowns[k]);
0347          CalibTowers->Fill(ietaC,iphiC,towerCal->get_energy(),scaleDowns[k]);
0348                  RawTowers->Fill(ietaR,iphiR,towerRaw->get_energy(),scaleDowns[k]);
0349               }
0350            }
0351         }
0352      }
0353      
0354      bool Eta_Check = check_bad_jet_eta(jet->get_eta(), vtx_z, 0.4);
0355 
0356      if (jet->get_pt() < All_RPt_Cut || jet->size_comp() <= NComp_Cut || Eta_Check) continue;
0357      h_theJetSpectrum -> Fill(jet->get_pt(), eventPrescale);
0358      
0359      for(int k = 0; k < nTrigs; k++){
0360        //std::cout << "on trig index: " << k << "; Index_plots[k] = " << Index_Plots[k] << std::endl;
0361        
0362         if(Index_Plots[k] != 1) continue;
0363     //std::cout << "Scaledown value: " << scaleDowns[k] << std::endl;
0364         if(scaleDowns[k])h_Eta_Phi_Pt_[k]->Fill(jet->get_eta(),jet->get_phi(),jet->get_pt(),liveCounts[k]/scaleDowns[k]);
0365         
0366         auto comp_vec = jet->get_comp_vec();
0367         std::vector<fastjet::PseudoJet> pseudojets;
0368         for(auto comp = comp_vec.begin(); comp != comp_vec.end(); ++comp){
0369            TowerInfo *tower;
0370            unsigned int channel = comp->second;
0371            unsigned int calo = comp->first;
0372           
0373            if (calo == 15 || calo == 30 || calo == 26){
0374               tower = CStowersIH3->get_tower_at_channel(channel);
0375               unsigned int calokey = CStowersIH3->encode_key(channel);
0376               int ieta = CStowersIH3->getTowerEtaBin(calokey);
0377               int iphi = CStowersIH3->getTowerPhiBin(calokey);
0378               h_iHCal_Jet_Eta_Phi_E_[k]->Fill(ieta,iphi,tower->get_energy(),scaleDowns[k]);
0379            }
0380            if (calo == 16 || calo == 31 || calo == 27){
0381               tower = CStowersOH3->get_tower_at_channel(channel);
0382               unsigned int calokey = CStowersOH3->encode_key(channel);
0383               int ieta = CStowersOH3->getTowerEtaBin(calokey);
0384               int iphi = CStowersOH3->getTowerPhiBin(calokey);
0385               h_oHCal_Jet_Eta_Phi_E_[k]->Fill(ieta,iphi,tower->get_energy(),scaleDowns[k]);
0386            }
0387            if (calo == 14 || calo == 29 || calo == 28){
0388               tower = CStowersEM3->get_tower_at_channel(channel);
0389               unsigned int calokey = CStowersEM3->encode_key(channel);
0390               int ieta = CStowersEM3->getTowerEtaBin(calokey);
0391               int iphi = CStowersEM3->getTowerPhiBin(calokey);
0392               h_EMCal_Jet_Eta_Phi_E_[k]->Fill(ieta,iphi,tower->get_energy(),scaleDowns[k]);
0393            }
0394         }
0395      }
0396   }
0397  
0398   return Fun4AllReturnCodes::EVENT_OK;
0399 }
0400 
0401 //____________________________________________________________________________..
0402 int Ana_PPG09_Mod::ResetEvent(PHCompositeNode *topNode)
0403 {
0404   if(Verbosity() > 5){
0405      std::cout << "Ana_PPG09_Mod::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0406   }
0407   return Fun4AllReturnCodes::EVENT_OK;
0408 }
0409 
0410 //____________________________________________________________________________..
0411 int Ana_PPG09_Mod::End(PHCompositeNode *topNode)
0412 {
0413   PHTFileServer::get().cd(m_outputFileName);
0414   std::cout << "Saving histograms" << std::endl;
0415   h_nJetsAboveThresh->Write();
0416   h_EventCount->Write();
0417   h_theJetSpectrum  -> Write();
0418   for(int k = 0; k < nTrigs; k++){
0419      h_Eta_Phi_Pt_[k]->Write();
0420      h_EMCal_Raw_Eta_Phi_E_[k]->Write();
0421      h_iHCal_Raw_Eta_Phi_E_[k]->Write();
0422      h_oHCal_Raw_Eta_Phi_E_[k]->Write();
0423      h_EMCal_CS_Eta_Phi_E_[k]->Write();
0424      h_iHCal_CS_Eta_Phi_E_[k]->Write();
0425      h_oHCal_CS_Eta_Phi_E_[k]->Write();
0426      h_EMCal_C_Eta_Phi_E_[k]->Write();
0427      h_EMCal_CR_Eta_Phi_E_[k]->Write();
0428      h_iHCal_C_Eta_Phi_E_[k]->Write();
0429      h_oHCal_C_Eta_Phi_E_[k]->Write();
0430      h_EMCal_Jet_Eta_Phi_E_[k]->Write();
0431      h_iHCal_Jet_Eta_Phi_E_[k]->Write();
0432      h_oHCal_Jet_Eta_Phi_E_[k]->Write();
0433      h_EMCal_TE_Sub_Eta_Phi_E_[k]->Write();
0434      h_iHCal_TE_Sub_Eta_Phi_E_[k]->Write();
0435      h_oHCal_TE_Sub_Eta_Phi_E_[k]->Write();
0436      h_ZVtx_[k]->Write();
0437   }
0438 
0439   scaleDowns.clear();
0440   liveCounts.clear();
0441 
0442   std::cout << "Ana_PPG09_Mod::End - Output to " << m_outputFileName << std::endl;
0443   if(Verbosity() > 5){
0444      std::cout << "Ana_PPG09_Mod::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0445   }
0446   return Fun4AllReturnCodes::EVENT_OK;
0447 }
0448 
0449 //____________________________________________________________________________..
0450 int Ana_PPG09_Mod::Reset(PHCompositeNode *topNode)
0451 {
0452   if(Verbosity() > 5){
0453      std::cout << "Ana_PPG09_Mod::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0454   }
0455   return Fun4AllReturnCodes::EVENT_OK;
0456 }
0457 //____________________________________________________________________________..
0458 void Ana_PPG09_Mod::initializePrescaleInformationFromDB(int runnumber)
0459 {
0460   std::cout << "Starting trig DB lookup for run "<< runnumber << std::endl;
0461   TSQLServer* db = TSQLServer::Connect("pgsql://sphnxdaqdbreplica:5432/daq","phnxro","");
0462   if (!db || db->IsZombie())
0463     {
0464       std::cerr <<  "[ERROR] DB connection failed for run "
0465         << runnumber << std::endl;
0466       if (db) delete db;
0467       return;
0468     }
0469 
0470   char query[512];
0471   snprintf(query, sizeof(query),
0472        "SELECT s.index, t.triggername, s.live, s.scaled "
0473        "FROM gl1_scalers s "
0474        "JOIN gl1_triggernames t ON (s.index = t.index AND s.runnumber BETWEEN t.runnumber AND t.runnumber_last) "
0475        "WHERE s.runnumber=%d ORDER BY s.index;", runnumber);
0476 
0477   auto *res = db->Query(query);
0478   if (!res)
0479     {
0480       std::cerr << "[ERROR] Query failed for run "
0481         << runnumber  << std::endl;
0482       delete db;
0483       return;
0484     }
0485     
0486   //grab the prescales for each of our triggers of interest
0487   //for(int i = 0; i < (int)sizeof(trigNames)/sizeof(trigNames[0]); ++i)
0488   std::cout << "Looking for " << trigNames.size() << " triggers" << std::endl;
0489   for(int i = 0; i < (int)trigNames.size(); i++)
0490     {
0491       bool trigFound = false;
0492       std::cout << "Looking for trig " << trigNames[i] << std::endl;
0493       while (auto row = res->Next())
0494     {
0495       const char* dbTrig = row->GetField(1);
0496       const char* liveStr = row->GetField(2);
0497       const char* scaledStr = row->GetField(3);
0498 
0499       if (!dbTrig || !liveStr || !scaledStr) { delete row; continue; }
0500 
0501       std::string trigName(dbTrig);
0502       double live = std::atof(liveStr);
0503       double scaled = std::atof(scaledStr);
0504       // check map
0505         
0506       //the trigger string matches one we're looking for
0507       //NB if for some reason there are multiple triggers of the same
0508       //name this would break
0509       if (!strcmp(dbTrig,trigNames[i].c_str()))
0510         {
0511           scaleDowns.push_back(scaled);
0512           liveCounts.push_back(live);
0513           std::cout << "Scale for trig " << trigNames[i] << " found with value " << scaled << std::endl;
0514           std::cout << "Live for trig " << trigNames[i] << " found with value " << live << std::endl;
0515           trigFound = true;
0516           delete row;
0517           break;
0518         }
0519     }
0520       
0521       if(!trigFound)
0522     {
0523       std::cout << "No trigger information for trigger " << trigNames[i] << " in run " << m_runnumber << std::endl;
0524       scaleDowns.push_back(0);
0525       liveCounts.push_back(0);
0526     }
0527     }
0528   delete res; delete db;
0529 }
0530 //____________________________________________________________________________..
0531 int Ana_PPG09_Mod::getEventPrescale(float leadJetPt, std::vector<int> triggerVector)
0532 {
0533   int prescale = 0;  
0534   //need to return pre-scaled from least pre-scaled, efficient trigger
0535   int lowest_prescale = std::numeric_limits<int>::max();
0536   bool trigIsGood = false;
0537   int live = 0;
0538   for(int i = 0; i < (int)triggerVector.size(); ++i)
0539     {
0540       if(std::find(trigIndices.begin(), trigIndices.end(), triggerVector[i]) == trigIndices.end())
0541     {//this trigger is not in our menu
0542       continue;
0543     }
0544       if(scaleDowns[trigToVecIndex[triggerVector[i]]] < lowest_prescale && isTrigEfficient(trigToVecIndex[triggerVector[i]],leadJetPt) && scaleDowns[trigToVecIndex[triggerVector[i]]] > 0)
0545     {
0546       lowest_prescale = scaleDowns[trigToVecIndex[triggerVector[i]]];
0547       live = liveCounts[trigToVecIndex[triggerVector[i]]];
0548       trigIsGood = true;
0549     }
0550     }
0551 
0552   //a fully efficient trigger that fired during this event has been found, and its  pre-scale is the lowest, use its pre-scale
0553   if(lowest_prescale != std::numeric_limits<int>::max() && trigIsGood && lowest_prescale != 0) prescale = live/lowest_prescale;
0554   //else we should default to the MB pre-scale if the MBD trigger fired
0555   else if(std::find(triggerVector.begin(), triggerVector.end(), 10) != triggerVector.end() && scaleDowns[0] != 0) prescale = liveCounts[0]/scaleDowns[0];
0556   //else this is an event where some rare-probe trigger fired without coincidence with the MBD, so we toss it (for now)
0557   else prescale = 0;
0558   return prescale;
0559 }
0560 //____________________________________________________________________________..
0561 bool Ana_PPG09_Mod::isTrigEfficient(int trigIndex, float leadJetPt)
0562 {
0563   
0564   if(leadJetPt > trigCutOffs[trigIndex]) return true;
0565   return true;
0566 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!