Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:09 

0001 #include "CaloAna.h"
0002 
0003 // G4Hits includes
0004 #include <TLorentzVector.h>
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 
0008 // G4Cells includes
0009 #include <g4detectors/PHG4Cell.h>
0010 #include <g4detectors/PHG4CellContainer.h>
0011 
0012 // Tower includes
0013 #include <calobase/RawCluster.h>
0014 #include <calobase/RawClusterContainer.h>
0015 #include <calobase/RawClusterUtility.h>
0016 #include <calobase/RawTower.h>
0017 #include <calobase/RawTowerContainer.h>
0018 #include <calobase/RawTowerGeom.h>
0019 #include <calobase/RawTowerGeomContainer.h>
0020 #include <calobase/TowerInfo.h>
0021 #include <calobase/TowerInfoContainer.h>
0022 #include <calobase/TowerInfoDefs.h>
0023 
0024 // Cluster includes
0025 #include <calobase/RawCluster.h>
0026 #include <calobase/RawClusterContainer.h>
0027 
0028 #include <fun4all/Fun4AllHistoManager.h>
0029 #include <fun4all/Fun4AllReturnCodes.h>
0030 
0031 #include <phool/getClass.h>
0032 
0033 #include <globalvertex/GlobalVertex.h>
0034 #include <globalvertex/GlobalVertexMap.h>
0035 
0036 // MBD
0037 #include <mbd/BbcGeom.h>
0038 #include <mbd/MbdPmtContainerV1.h>
0039 #include <mbd/MbdPmtHit.h>
0040 
0041 #include <TFile.h>
0042 #include <TH1.h>
0043 #include <TH2.h>
0044 #include <TNtuple.h>
0045 #include <TTree.h>
0046 #include <TProfile2D.h>
0047 
0048 #include <Event/Event.h>
0049 #include <Event/packet.h>
0050 #include <cassert>
0051 #include <sstream>
0052 #include <string>
0053 
0054 TH2F* LogYHist2D(const char* name, const char* title, int xbins_in, double_t xmin, double_t xmax, int ybins_in, double_t ymin, double_t ymax);
0055 
0056 using namespace std;
0057 
0058 CaloAna::CaloAna(const std::string& name, const std::string& filename)
0059   : SubsysReco(name)
0060   , detector("HCALIN")
0061   , outfilename(filename)
0062 {
0063   _eventcounter = 0;
0064 }
0065 
0066 CaloAna::~CaloAna()
0067 {
0068   delete hm;
0069   delete g4hitntuple;
0070   delete g4cellntuple;
0071   delete towerntuple;
0072   delete clusterntuple;
0073 }
0074 
0075 int CaloAna::Init(PHCompositeNode*)
0076 {
0077   hm = new Fun4AllHistoManager(Name());
0078   // create and register your histos (all types) here
0079 
0080   outfile = new TFile(outfilename.c_str(), "RECREATE");
0081   // correlation plots
0082 
0083   h_emcal_mbd_correlation = new TH2F("h_emcal_mbd_correlation", ";emcal;mbd", 100, 0, 1, 100, 0, 1);
0084   h_ohcal_mbd_correlation = new TH2F("h_ohcal_mbd_correlation", ";ohcal;mbd", 100, 0, 1, 100, 0, 1);
0085   h_ihcal_mbd_correlation = new TH2F("h_ihcal_mbd_correlation", ";ihcal;mbd", 100, 0, 1, 100, 0, 1);
0086   h_emcal_hcal_correlation = new TH2F("h_emcal_hcal_correlation", ";emcal;hcal", 100, 0, 1, 100, 0, 1);
0087   h_cemc_etaphi = new TH2F("h_cemc_etaphi", ";eta;phi", 96, 0, 96, 256, 0, 256);
0088   h_hcalin_etaphi = new TH2F("h_ihcal_etaphi", ";eta;phi", 24, 0, 24, 64, 0, 64);
0089   h_hcalout_etaphi = new TH2F("h_ohcal_etaphi", ";eta;phi", 24, 0, 24, 64, 0, 64);
0090 
0091   h_cemc_etaphi_wQA = new TH2F("h_cemc_etaphi_wQA", ";eta;phi", 96, 0, 96, 256, 0, 256);
0092   h_hcalin_etaphi_wQA = new TH2F("h_ihcal_etaphi_wQA", ";eta;phi", 24, 0, 24, 64, 0, 64);
0093   h_hcalout_etaphi_wQA = new TH2F("h_ohcal_etaphi_wQA", ";eta;phi", 24, 0, 24, 64, 0, 64);
0094 
0095   h_cemc_e_chi2 =  LogYHist2D("h_cemc_e_chi2", "", 500,-2,30, 1000, 0.5, 5e6);
0096   h_ihcal_e_chi2 =  LogYHist2D("h_ihcal_e_chi2", "", 500,-2,30, 1000, 0.5, 5e6);
0097   h_ohcal_e_chi2 =  LogYHist2D("h_ohcal_e_chi2", "", 500,-2,30, 1000, 0.5, 5e6);
0098 
0099   h_cemc_etaphi_time = new TProfile2D("h_cemc_etaphi_time", ";eta;phi", 96, 0, 96, 256, 0, 256,-10,10);
0100   h_hcalin_etaphi_time = new TProfile2D("h_ihcal_etaphi_time", ";eta;phi", 24, 0, 24, 64, 0, 64 ,-10,10);
0101   h_hcalout_etaphi_time = new TProfile2D("h_ohcal_etaphi_time", ";eta;phi", 24, 0, 24, 64, 0, 64 ,-10,10);
0102 
0103   h_cemc_etaphi_badChi2 = new TProfile2D("h_cemc_etaphi_badChi2", ";eta;phi", 96, 0, 96, 256, 0, 256,-10,10);
0104   h_hcalin_etaphi_badChi2 = new TProfile2D("h_ihcal_etaphi_badChi2", ";eta;phi", 24, 0, 24, 64, 0, 64 ,-10,10);
0105   h_hcalout_etaphi_badChi2 = new TProfile2D("h_ohcal_etaphi_badChi2", ";eta;phi", 24, 0, 24, 64, 0, 64 ,-10,10);
0106 
0107 
0108   // 1D distributions
0109   h_InvMass = new TH1F("h_InvMass", "Invariant Mass", 120, 0, 1.2);
0110 
0111   // ZDC QA plots
0112   hzdcSouthraw = new TH1D("hzdcSouthraw", "hzdcSouthraw", 1500, 0, 15000);
0113   hzdcNorthraw = new TH1D("hzdcNorthraw", "hzdcNorthraw", 1500, 0, 15000);
0114   hzdcSouthcalib = new TH1D("hzdcSouthcalib", "hzdcSouthcalib", 1500, 0, 15000);
0115   hzdcNorthcalib = new TH1D("hzdcNorthcalib", "hzdcNorthcalib", 1500, 0, 15000);
0116   h_totalzdc_e = new TH1D("h_totalzdc_e", "", 200, 0, 2e4);
0117   h_emcal_zdc_correlation = new TH2F("h_zdc_emcal_correlation", ";emcal;zdc", 100, 0, 1, 100, 0, 1);
0118 
0119   // vertex distributions
0120   hvtx_z_raw = new TH1D("hvtx_z_raw", "hvtx_z_raw", 201, -100.5, 100.5);
0121   hvtx_z_cut = new TH1D("hvtx_z_cut", "hvtx_z_cut", 201, -100.5, 100.5);
0122 
0123   // raw timing information
0124   hzdctime_cut = new TH1D("hzdctime_cut", "hzdctime_cut", 50, -17.5, 32.5);
0125   hemcaltime_cut = new TH1D("hemcaltime_cut", "hemcaltime_cut", 50, -17.5, 32.5);
0126   hihcaltime_cut = new TH1D("hihcaltime_cut", "hihcaltime_cut", 50, -17.5, 32.5);
0127   hohcaltime_cut = new TH1D("hohcaltime_cut", "hohcaltime_cut", 50, -17.5, 32.5);
0128 
0129   // extracted timing information
0130   hzdctime = new TH1D("hzdctime", "hzdctime", 50, -17.5, 32.5);
0131   hemcaltime = new TH1D("hemcaltime", "hemcaltime", 50, -17.5, 32.5);
0132   hihcaltime = new TH1D("hihcaltime", "hihcaltime", 50, -17.5, 32.5);
0133   hohcaltime = new TH1D("hohcaltime", "hohcaltime", 50, -17.5, 32.5);
0134 
0135   // cluster QA
0136   h_etaphi_clus = new TH2F("h_etaphi_clus", "", 140, -1.2, 1.2, 64, -1 * TMath::Pi(), TMath::Pi());
0137   h_clusE = new TH1F("h_clusE", "", 100, 0, 10);
0138 
0139   return 0;
0140 }
0141 
0142 int CaloAna::process_event(PHCompositeNode* topNode)
0143 {
0144   _eventcounter++;
0145 
0146   process_towers(topNode);
0147 
0148   return Fun4AllReturnCodes::EVENT_OK;
0149 }
0150 
0151 int CaloAna::process_towers(PHCompositeNode* topNode)
0152 {
0153   std::cout << _eventcounter << std::endl;
0154 
0155   float totalcemc = 0.;
0156   float totalihcal = 0.;
0157   float totalohcal = 0.;
0158   float totalmbd = 0.;
0159   float totalzdc = 0.;
0160   float totalzdcsouthraw = 0.;
0161   float totalzdcnorthraw = 0.;
0162   float totalzdcsouthcalib = 0.;
0163   float totalzdcnorthcalib = 0.;
0164 
0165   float emcaldownscale = 1000000 / 800;
0166   float ihcaldownscale = 40000 / 300;
0167   float ohcaldownscale = 250000 / 600;
0168   float mbddownscale = 2000.0;
0169   float zdcdownscale = 1e4;
0170 
0171   float emcal_hit_threshold = 0.5; // GeV
0172   float ohcal_hit_threshold = 0.5;
0173   float ihcal_hit_threshold = 0.25;
0174 
0175   int max_zdc_t = -1;
0176   int max_emcal_t = -1;
0177   int max_ihcal_t = -1;
0178   int max_ohcal_t = -1;
0179 
0180   // get time estimate
0181   max_zdc_t = Getpeaktime(hzdctime_cut);
0182   max_emcal_t = Getpeaktime(hemcaltime_cut);
0183   max_ihcal_t = Getpeaktime(hihcaltime_cut);
0184   max_ohcal_t = Getpeaktime(hohcaltime_cut);
0185 
0186   //----------------------------------get vertex------------------------------------------------------//
0187 
0188   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0189   if (!vertexmap)
0190   {
0191     // std::cout << PHWHERE << " Fatal Error - GlobalVertexMap node is missing"<< std::endl;
0192     std::cout << "CaloAna GlobalVertexMap node is missing" << std::endl;
0193     // return Fun4AllReturnCodes::ABORTRUN;
0194   }
0195   float vtx_z = NAN;
0196   if (vertexmap && !vertexmap->empty())
0197   {
0198     GlobalVertex* vtx = vertexmap->begin()->second;
0199     if (vtx)
0200     {
0201       vtx_z = vtx->get_z();
0202     }
0203     hvtx_z_raw->Fill(vtx_z);
0204   }
0205 
0206   vector<float> ht_eta;
0207   vector<float> ht_phi;
0208 
0209   if (!m_vtxCut || abs(vtx_z) < _vz)
0210   {
0211     hvtx_z_cut->Fill(vtx_z);
0212 
0213     //----------------------------------------------tower energies -----------------------------------------------//
0214 
0215     {
0216       TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0217       if (towers)
0218       {
0219         int size = towers->size();  // online towers should be the same!
0220         for (int channel = 0; channel < size; channel++)
0221         {
0222           TowerInfo* tower = towers->get_tower_at_channel(channel);
0223           float offlineenergy = tower->get_energy();
0224           unsigned int towerkey = towers->encode_key(channel);
0225           int ieta = towers->getTowerEtaBin(towerkey);
0226           int iphi = towers->getTowerPhiBin(towerkey);
0227           int _time = tower->get_time();
0228           h_cemc_e_chi2->Fill( offlineenergy,tower->get_chi2());
0229           float _timef = tower->get_time_float();
0230           hemcaltime_cut->Fill(_time);
0231           bool isGood = ! (tower->get_isBadChi2());
0232           if (!isGood && offlineenergy > 0.2) 
0233           {
0234             ht_eta.push_back(ieta);
0235             ht_phi.push_back(iphi);
0236           }
0237 
0238           if (_time > (max_emcal_t - _range) && _time < (max_emcal_t + _range))
0239           {
0240             totalcemc += offlineenergy;
0241             hemcaltime->Fill(_time);
0242             if (offlineenergy > emcal_hit_threshold)
0243             {
0244               h_cemc_etaphi_time->Fill(ieta, iphi, _timef);
0245               h_cemc_etaphi->Fill(ieta, iphi);
0246               if (isGood) h_cemc_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
0247               if (tower->get_isBadChi2()){
0248                 h_cemc_etaphi_badChi2->Fill(ieta,iphi,1);
0249               }
0250               else{
0251                 h_cemc_etaphi_badChi2->Fill(ieta,iphi,0);
0252               }
0253             }
0254           }
0255         }
0256       }
0257     }
0258 
0259     {
0260       TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0261       if (towers)
0262       {
0263         int size = towers->size();  // online towers should be the same!
0264         for (int channel = 0; channel < size; channel++)
0265         {
0266           TowerInfo* tower = towers->get_tower_at_channel(channel);
0267           float offlineenergy = tower->get_energy();
0268           unsigned int towerkey = towers->encode_key(channel);
0269           int ieta = towers->getTowerEtaBin(towerkey);
0270           int iphi = towers->getTowerPhiBin(towerkey);
0271           int _time = tower->get_time();
0272           float _timef = tower->get_time_float();
0273           hihcaltime_cut->Fill(_time);
0274           h_ihcal_e_chi2->Fill( offlineenergy,tower->get_chi2());
0275           bool isGood = !(tower->get_isBadChi2());
0276           if (_time > (max_ihcal_t - _range) && _time < (max_ihcal_t + _range))
0277           {
0278             totalihcal += offlineenergy;
0279             hihcaltime->Fill(_time);
0280 
0281             if (offlineenergy > ihcal_hit_threshold)
0282             {
0283               h_hcalin_etaphi->Fill(ieta, iphi);
0284               h_hcalin_etaphi_time->Fill(ieta, iphi, _timef);
0285               if (isGood) h_hcalin_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
0286               if (tower->get_isBadChi2()){
0287                 h_hcalin_etaphi_badChi2->Fill(ieta,iphi,1);
0288               }
0289               else{
0290                 h_hcalin_etaphi_badChi2->Fill(ieta,iphi,0);
0291               }
0292             }
0293           }
0294         }
0295       }
0296     }
0297     {
0298       TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0299       if (towers)
0300       {
0301         int size = towers->size();  // online towers should be the same!
0302         for (int channel = 0; channel < size; channel++)
0303         {
0304           TowerInfo* tower = towers->get_tower_at_channel(channel);
0305           float offlineenergy = tower->get_energy();
0306           unsigned int towerkey = towers->encode_key(channel);
0307           int ieta = towers->getTowerEtaBin(towerkey);
0308           int iphi = towers->getTowerPhiBin(towerkey);
0309           int _time = tower->get_time();
0310           float _timef = tower->get_time_float();
0311           hihcaltime_cut->Fill(_time);
0312           hohcaltime_cut->Fill(_time);
0313           h_ohcal_e_chi2->Fill( offlineenergy,tower->get_chi2());
0314           bool isGood = !(tower->get_isBadChi2());
0315 
0316           if (_time > (max_ohcal_t - _range) && _time < (max_ohcal_t + _range))
0317           {
0318             totalohcal += offlineenergy;
0319             hohcaltime->Fill(_time);
0320 
0321             if (offlineenergy > ohcal_hit_threshold)
0322             {
0323               h_hcalout_etaphi->Fill(ieta, iphi);
0324               h_hcalout_etaphi_time->Fill(ieta, iphi, _timef);
0325               if (isGood) h_hcalout_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
0326               if (tower->get_isBadChi2()){
0327                 h_hcalout_etaphi_badChi2->Fill(ieta,iphi,1);
0328               }
0329               else{
0330                 h_hcalout_etaphi_badChi2->Fill(ieta,iphi,0);
0331               }
0332             }
0333           }
0334         }
0335       }
0336     }
0337 
0338     {
0339       TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_ZDC");
0340       if (towers)
0341       {
0342         int size = towers->size();  // online towers should be the same!
0343         for (int channel = 0; channel < size; channel++)
0344         {
0345           TowerInfo* tower = towers->get_tower_at_channel(channel);
0346           float offlineenergy = tower->get_energy();
0347           int _time = towers->get_tower_at_channel(channel)->get_time();
0348           hzdctime_cut->Fill(_time);
0349           if (channel == 0 || channel == 2 || channel == 4)
0350           {
0351             totalzdcsouthcalib += offlineenergy;
0352           }
0353           if (channel == 8 || channel == 10 || channel == 12)
0354           {
0355             totalzdcnorthcalib += offlineenergy;
0356           }
0357           if (channel == 0 || channel == 2 || channel == 4 || channel == 8 || channel == 10 || channel == 12)
0358           {
0359             if (_time > (max_zdc_t - _range) && _time < (max_zdc_t + _range))
0360             {
0361               totalzdc += offlineenergy;
0362               hzdctime->Fill(_time);
0363             }
0364           }
0365         }
0366       }
0367     }
0368 
0369     {
0370       TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_ZDC");
0371       if (towers)
0372       {
0373         int size = towers->size();  // online towers should be the same!
0374         for (int channel = 0; channel < size; channel++)
0375         {
0376           TowerInfo* tower = towers->get_tower_at_channel(channel);
0377           float offlineenergy = tower->get_energy();
0378           if (channel == 0 || channel == 2 || channel == 4)
0379           {
0380             totalzdcsouthraw += offlineenergy;
0381           }
0382           if (channel == 8 || channel == 10 || channel == 12)
0383           {
0384             totalzdcnorthraw += offlineenergy;
0385           }
0386         }
0387       }
0388     }
0389     //--------------------------------------- MBD ---------------------------------------------------//
0390     MbdPmtContainer* bbcpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0391     if (!bbcpmts)
0392     {
0393       std::cout << "makeMBDTrees::process_event: Could not find MbdPmtContainer, aborting" << std::endl;
0394       return Fun4AllReturnCodes::ABORTEVENT;
0395     }
0396 
0397     int nPMTs = bbcpmts->get_npmt();
0398     for (int i = 0; i < nPMTs; i++)
0399     {
0400       MbdPmtHit* mbdpmt = bbcpmts->get_pmt(i);
0401       float pmtadc = mbdpmt->get_q();
0402       totalmbd += pmtadc;
0403     }
0404 
0405     h_emcal_mbd_correlation->Fill(totalcemc / emcaldownscale, totalmbd / mbddownscale);
0406     h_ihcal_mbd_correlation->Fill(totalihcal / ihcaldownscale, totalmbd / mbddownscale);
0407     h_ohcal_mbd_correlation->Fill(totalohcal / ohcaldownscale, totalmbd / mbddownscale);
0408     h_emcal_hcal_correlation->Fill(totalcemc / emcaldownscale, totalohcal / ohcaldownscale);
0409     h_emcal_zdc_correlation->Fill(totalcemc / emcaldownscale, totalzdc / zdcdownscale);
0410     h_totalzdc_e->Fill(totalzdc);
0411 
0412     hzdcSouthraw->Fill(totalzdcsouthraw);
0413     hzdcNorthraw->Fill(totalzdcnorthraw);
0414     hzdcSouthcalib->Fill(totalzdcsouthcalib);
0415     hzdcNorthcalib->Fill(totalzdcnorthcalib);
0416 
0417   }  // vertex cut
0418   //------------------------------------------------- pi 0 --------------------------------------------------------//
0419 
0420   RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_POS_COR_CEMC");
0421   if (!clusterContainer)
0422   {
0423     std::cout << PHWHERE << "funkyCaloStuff::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0424     return 0;
0425   }
0426 
0427   //////////////////////////////////////////
0428   // geometry for hot tower/cluster masking
0429   std::string towergeomnodename = "TOWERGEOM_CEMC";
0430   RawTowerGeomContainer* m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
0431   if (!m_geometry)
0432   {
0433     std::cout << Name() << "::" 
0434               << "CreateNodeTree"
0435               << ": Could not find node " << towergeomnodename << std::endl;
0436     throw std::runtime_error("failed to find TOWERGEOM node in RawClusterDeadHotMask::CreateNodeTree");
0437   }
0438 
0439   // cuts
0440   float emcMinClusE1 = 1.5;  // 0.5;
0441   float emcMinClusE2 = 0.8;  // 0.5;
0442   float emcMaxClusE = 32;
0443   float minDr = 0.08;
0444   float maxAlpha = 0.8;
0445 
0446   // if (totalmbd < 0.1 * mbddownscale  )
0447   if (totalcemc < 0.2 * emcaldownscale)
0448   {
0449     RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0450     RawClusterContainer::ConstIterator clusterIter;
0451     RawClusterContainer::ConstIterator clusterIter2;
0452 
0453     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0454     {
0455       RawCluster* recoCluster = clusterIter->second;
0456 
0457       CLHEP::Hep3Vector vertex(0, 0, 0);
0458       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0459 
0460       float clusE = E_vec_cluster.mag();
0461       float clus_eta = E_vec_cluster.pseudoRapidity();
0462       float clus_phi = E_vec_cluster.phi();
0463       float clus_pt = E_vec_cluster.perp();
0464       float clus_chisq = recoCluster->get_chi2();
0465 
0466       h_clusE->Fill(clusE);
0467 
0468       if (clusE < emcMinClusE1 || clusE > emcMaxClusE)
0469       {
0470         continue;
0471       }
0472       if (abs(clus_eta) > 0.7)
0473       {
0474         continue;
0475       }
0476       if (clus_chisq > 4)
0477       {
0478         continue;
0479       }
0480 
0481       if (dynMaskClus)
0482       {
0483         // loop over the towers in the cluster
0484         RawCluster::TowerConstRange towers = recoCluster->get_towers();
0485         RawCluster::TowerConstIterator toweriter;
0486         bool hotClus = false;
0487         for (toweriter = towers.first; toweriter != towers.second; ++toweriter)
0488         {
0489           int towereta = m_geometry->get_tower_geometry(toweriter->first)->get_bineta();
0490           int towerphi = m_geometry->get_tower_geometry(toweriter->first)->get_binphi();
0491           for (size_t i = 0; i < ht_eta.size(); i++)
0492           {
0493             if (towerphi == ht_phi[i] && towereta == ht_eta[i]) hotClus = true;
0494           }
0495         }
0496         if (hotClus == true) continue;
0497       }
0498 
0499       h_etaphi_clus->Fill(clus_eta, clus_phi);
0500 
0501       TLorentzVector photon1;
0502       photon1.SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clusE);
0503 
0504       for (clusterIter2 = clusterEnd.first; clusterIter2 != clusterEnd.second; clusterIter2++)
0505       {
0506         if (clusterIter == clusterIter2)
0507         {
0508           continue;
0509         }
0510         RawCluster* recoCluster2 = clusterIter2->second;
0511 
0512         CLHEP::Hep3Vector vertex2(0, 0, 0);
0513         CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*recoCluster2, vertex2);
0514 
0515         float clus2E = E_vec_cluster2.mag();
0516         float clus2_eta = E_vec_cluster2.pseudoRapidity();
0517         float clus2_phi = E_vec_cluster2.phi();
0518         float clus2_pt = E_vec_cluster2.perp();
0519         float clus2_chisq = recoCluster2->get_chi2();
0520 
0521         if (clus2E < emcMinClusE2 || clus2E > emcMaxClusE)
0522         {
0523           continue;
0524         }
0525         if (abs(clus2_eta) > 0.7)
0526         {
0527           continue;
0528         }
0529         if (clus2_chisq > 4)
0530         {
0531           continue;
0532         }
0533 
0534         if (dynMaskClus)
0535         {
0536           // loop over the towers in the cluster
0537           RawCluster::TowerConstRange towers2 = recoCluster2->get_towers();
0538           RawCluster::TowerConstIterator toweriter2;
0539           bool hotClus = false;
0540           for (toweriter2 = towers2.first; toweriter2 != towers2.second; ++toweriter2)
0541           {
0542             int towereta = m_geometry->get_tower_geometry(toweriter2->first)->get_bineta();
0543             int towerphi = m_geometry->get_tower_geometry(toweriter2->first)->get_binphi();
0544 
0545             for (size_t i = 0; i < ht_eta.size(); i++)
0546             {
0547               if (towerphi == ht_phi[i] && towereta == ht_phi[i]) hotClus = true;
0548             }
0549           }
0550           if (hotClus == true) continue;
0551         }
0552 
0553         TLorentzVector photon2;
0554         photon2.SetPtEtaPhiE(clus2_pt, clus2_eta, clus2_phi, clus2E);
0555 
0556         if (sqrt(pow(clusE - clus2E, 2)) / (clusE + clus2E) > maxAlpha) continue;
0557 
0558         if (photon1.DeltaR(photon2) < minDr) continue;
0559         TLorentzVector pi0 = photon1 + photon2;
0560         h_InvMass->Fill(pi0.M());
0561       }
0562     }
0563   }
0564 
0565   return Fun4AllReturnCodes::EVENT_OK;
0566 }
0567 
0568 int CaloAna::End(PHCompositeNode* /*topNode*/)
0569 {
0570   outfile->cd();
0571 
0572   outfile->Write();
0573   outfile->Close();
0574   delete outfile;
0575   hm->dumpHistos(outfilename, "UPDATE");
0576   return 0;
0577 }
0578 
0579 int CaloAna::Getpeaktime(TH1* h)
0580 {
0581   int getmaxtime, tcut = -1;
0582 
0583   for (int bin = 1; bin < h->GetNbinsX() + 1; bin++)
0584   {
0585     double c = h->GetBinContent(bin);
0586     double max = h->GetMaximum();
0587     int bincenter = h->GetBinCenter(bin);
0588     if (max == c)
0589     {
0590       getmaxtime = bincenter;
0591       if (getmaxtime != -1) tcut = getmaxtime;
0592     }
0593   }
0594 
0595   return tcut;
0596 }
0597 
0598 
0599 
0600 TH2F* LogYHist2D(const char* name, const char* title, int xbins_in, double_t xmin, double_t xmax, int ybins_in, double_t ymin, double_t ymax)    {
0601   Double_t logymin = TMath::Log10(ymin);
0602   Double_t logymax = TMath::Log10(ymax);
0603   Double_t binwidth = (logymax - logymin) / ybins_in;
0604   Double_t ybins[ybins_in + 1];
0605 
0606   for (Int_t i = 0; i <= ybins_in + 1; i++)
0607     ybins[i] = TMath::Power(10, logymin + i * binwidth);
0608 
0609   TH2F* h = new TH2F(name, title, xbins_in, xmin, xmax, ybins_in, ybins);
0610 
0611   return h;
0612 }
0613
[ 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!