Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "pi0EtaByEta.h"
0002 
0003 #include <g4main/PHG4TruthInfoContainer.h>
0004 #include <g4main/PHG4VtxPoint.h>
0005 
0006 #include <globalvertex/GlobalVertex.h>
0007 #include <globalvertex/GlobalVertexMap.h>
0008 #include <globalvertex/Vertex.h>
0009 
0010 // Tower includes
0011 #include <calobase/RawCluster.h>
0012 #include <calobase/RawClusterContainer.h>
0013 #include <calobase/RawClusterUtility.h>
0014 #include <calobase/RawTowerGeomContainer.h>
0015 #include <calobase/TowerInfo.h>
0016 #include <calobase/TowerInfoContainer.h>
0017 #include <calobase/TowerInfoDefs.h>
0018 
0019 #include <calotrigger/TriggerAnalyzer.h>
0020 #include <ffarawobjects/Gl1Packet.h>
0021 
0022 #include <cdbobjects/CDBTTree.h>  // for CDBTTree
0023 
0024 #include <fun4all/Fun4AllHistoManager.h>
0025 #include <fun4all/Fun4AllReturnCodes.h>
0026 #include <fun4all/SubsysReco.h>  // for SubsysReco
0027 
0028 #include <phool/getClass.h>
0029 #include <phool/phool.h>
0030 
0031 #include <TF1.h>
0032 #include <TFile.h>
0033 #include <TH1.h>
0034 #include <TH2.h>
0035 #include <TH3.h>
0036 #include <TLorentzVector.h>
0037 #include <TNtuple.h>
0038 #include <TStyle.h>
0039 #include <TSystem.h>
0040 #include <TTree.h>
0041 
0042 #include <CLHEP/Vector/ThreeVector.h>  // for Hep3Vector
0043 
0044 #include <cmath>    // for fabs, isnan, M_PI
0045 #include <cstdint>  // for exit
0046 #include <cstdlib>  // for exit
0047 #include <filesystem>
0048 #include <iostream>
0049 #include <map>        // for operator!=, _Rb_tree_con...
0050 #include <stdexcept>  // for runtime_error
0051 #include <string>
0052 #include <utility>
0053 
0054 pi0EtaByEta::pi0EtaByEta(const std::string& name, const std::string& filename)
0055   : SubsysReco(name)
0056   , detector("HCALIN")
0057   , outfilename(filename)
0058 {
0059   h_mass_eta_lt.fill(nullptr);
0060 
0061   if (runTowByTow)
0062   {
0063     for (auto& row : h_mass_tbt_lt)
0064     {
0065       row.fill(nullptr);
0066     }
0067   }
0068 
0069   clusMix = new std::vector<std::vector<std::vector<CLHEP::Hep3Vector>>>();
0070 }
0071 
0072 pi0EtaByEta::~pi0EtaByEta()
0073 {
0074   delete hm;
0075   delete g4hitntuple;
0076   delete g4cellntuple;
0077   delete towerntuple;
0078   delete clusterntuple;
0079 }
0080 
0081 int pi0EtaByEta::Init(PHCompositeNode* /*unused*/)
0082 {
0083   hm = new Fun4AllHistoManager(Name());
0084   // create and register your histos (all types) here
0085 
0086   outfile = new TFile(outfilename.c_str(), "RECREATE");
0087 
0088   // correlation plots
0089   for (int i = 0; i < 96; i++)
0090   {
0091     std::string histoname = "h_mass_eta_lt" + std::to_string(i);
0092     h_mass_eta_lt[i] = new TH1F(histoname.c_str(), "", 50, 0, 0.5);
0093 
0094     if (runTowByTow)
0095     {
0096       for (int j = 0; j < 256; j++)
0097       {
0098         std::string histoname_tbt = "h_mass_tbt_lt_" + std::to_string(i) + "_" + std::to_string(j);
0099         h_mass_tbt_lt[i][j] = new TH1F(histoname_tbt.c_str(), "", 50, 0, 0.5);
0100       }
0101     }
0102   }
0103 
0104   // 3D hist to save inv mass for all towers
0105   if (runTBTCompactMode)
0106   {
0107     h_ieta_iphi_invmass = new TH3F("h_ieta_iphi_invmass", "", 96, 0, 96, 256, 0, 256, 50, 0.0, 0.5);
0108   }
0109 
0110   h_cemc_etaphi = new TH2F("h_cemc_etaphi", "", 96, 0, 96, 256, 0, 256);
0111   h_cemc_etaphi_noCalib = new TH2F("h_cemc_etaphi_noCalib", "", 96, 0, 96, 256, 0, 256);
0112 
0113   // 1D distributions
0114   h_InvMass = new TH1F("h_InvMass", "Invariant Mass", 240, 0, 1.2);
0115   h_InvMassMix = new TH1F("h_InvMassMix", "Invariant Mass", 120, 0, 1.2);
0116 
0117   // cluster QA
0118   h_etaphi_clus = new TH2F("h_etaphi_clus", "", 140, -1.2, 1.2, 64, -1 * M_PI, M_PI);
0119   h_clus_pt = new TH1F("h_clus_pt", "", 100, 0, 10);
0120 
0121   h_emcal_e_eta = new TH1F("h_emcal_e_eta", "", 96, 0, 96);
0122 
0123   h_pt1 = new TH1F("h_pt1", "", 100, 0, 5);
0124   h_pt2 = new TH1F("h_pt2", "", 100, 0, 5);
0125 
0126   h_nclusters = new TH1F("h_nclusters", "", 1000, 0, 1000);
0127 
0128   h_event = new TH1F("h_event", "", 1, 0, 1);
0129 
0130   std::vector<std::vector<CLHEP::Hep3Vector>> temp2 = std::vector<std::vector<CLHEP::Hep3Vector>>();
0131   std::vector<CLHEP::Hep3Vector> temp = std::vector<CLHEP::Hep3Vector>();
0132   for (int i = 0; i < NBinsVtx; i++)
0133   {
0134     temp2.push_back(temp);
0135   }
0136   for (int i = 0; i < NBinsClus; i++)
0137   {
0138     clusMix->push_back(temp2);
0139   }
0140 
0141   trigAna = new TriggerAnalyzer();
0142 
0143   return 0;
0144 }
0145 
0146 int pi0EtaByEta::process_event(PHCompositeNode* topNode)
0147 {
0148   _eventcounter++;
0149 
0150   process_towers(topNode);
0151 
0152   return Fun4AllReturnCodes::EVENT_OK;
0153 }
0154 
0155 int pi0EtaByEta::process_towers(PHCompositeNode* topNode)
0156 {
0157   if ((_eventcounter % 1000) == 0)
0158   {
0159     std::cout << _eventcounter << std::endl;
0160   }
0161 
0162   // float emcaldownscale = 1000000 / 800;
0163 
0164   float emcal_hit_threshold = 0.3;  // GeV
0165 
0166   // cuts
0167   float maxDr = 1.1;
0168   float maxAlpha = 0.6;
0169   float clus_chisq_cut = 0.05;
0170   float nClus_ptCut = 0.5;
0171   int max_nClusCount = 300;
0172 
0173   //--------------------------- trigger and GL1-------------------------------//
0174 
0175   if (reqTrig)
0176   {
0177     trigAna->decodeTriggers(topNode);
0178     bool fireTrig = false;
0179     for (auto bit : triggerList)
0180     {
0181       if (trigAna->didTriggerFire(bit) == true)
0182       {
0183         fireTrig = true;
0184       }
0185     }
0186     if (!fireTrig)
0187     {
0188       return Fun4AllReturnCodes::EVENT_OK;
0189     }
0190   }
0191 
0192   //----------------------------------get vertex------------------------------------------------------//
0193   GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0194   if (!vertexmap)
0195   {
0196     // std::cout << PHWHERE << " Fatal Error - GlobalVertexMap node is missing"<< std::endl;
0197     std::cout << "pi0EtaByEta GlobalVertexMap node is missing" << std::endl;
0198     // return Fun4AllReturnCodes::ABORTRUN;
0199   }
0200 
0201   float vtx_z = 0;
0202   bool found_vertex = false;
0203   if (vertexmap && !vertexmap->empty())
0204   {
0205     GlobalVertex* vtx = vertexmap->begin()->second;
0206     if (vtx)
0207     {
0208       if (m_use_vertextype)
0209       {
0210         auto typeStartIter = vtx->find_vertexes(m_vertex_type);
0211         auto typeEndIter = vtx->end_vertexes();
0212         for (auto iter = typeStartIter; iter != typeEndIter; ++iter)
0213         {
0214           const auto& [type, vertexVec] = *iter;
0215           if (type != m_vertex_type)
0216           {
0217             continue;
0218           }
0219           for (const auto* vertex : vertexVec)
0220           {
0221             if (!vertex)
0222             {
0223               continue;
0224             }
0225             vtx_z = vertex->get_z();
0226             found_vertex = true;
0227           }
0228         }
0229       }
0230       else
0231       {
0232         vtx_z = vtx->get_z();
0233         found_vertex = true;
0234       }
0235     }
0236   }
0237 
0238   //////////////////////////////                                                                  // truth vertex
0239   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0240   float vtx_z_tr = -999;
0241   if (truthinfo)
0242   {
0243     PHG4TruthInfoContainer::VtxRange vtxrange = truthinfo->GetVtxRange();
0244     for (PHG4TruthInfoContainer::ConstVtxIterator iter = vtxrange.first; iter != vtxrange.second; ++iter)
0245     {
0246       PHG4VtxPoint* vtx_tr = iter->second;
0247       if (vtx_tr->get_id() == 1)
0248       {
0249         vtx_z_tr = vtx_tr->get_z();
0250       }
0251     }
0252     if (vtx_z_tr != -999)
0253     {
0254       vtx_z = vtx_z_tr;
0255     }
0256   }
0257 
0258   if (useVertexTruth == true)
0259   {
0260     vtx_z = vtx_z_tr;
0261   }
0262 
0263   if (!found_vertex && reqVertex)
0264   {
0265     return Fun4AllReturnCodes::EVENT_OK;
0266   }
0267 
0268   TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0269   if (towers)
0270   {
0271     int size = towers->size();
0272     for (int channel = 0; channel < size; channel++)
0273     {
0274       TowerInfo* tower = towers->get_tower_at_channel(channel);
0275       float offlineenergy = tower->get_energy();
0276       unsigned int towerkey = towers->encode_key(channel);
0277       int ieta = towers->getTowerEtaBin(towerkey);
0278       int iphi = towers->getTowerPhiBin(towerkey);
0279       bool isGood = tower->get_isGood();
0280       if (offlineenergy > emcal_hit_threshold && isGood)
0281       {
0282         h_cemc_etaphi->Fill(ieta, iphi);
0283       }
0284       if (tower->get_isNoCalib())
0285       {
0286         h_cemc_etaphi_noCalib->Fill(ieta, iphi);
0287       }
0288     }
0289   }
0290 
0291   std::string cluster_node_name = "CLUSTERINFO_CEMC";
0292 
0293   if (use_pdc)
0294   {
0295     cluster_node_name = "CLUSTERINFO_POS_COR_CEMC";
0296   }
0297   RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, cluster_node_name);
0298   if (!clusterContainer)
0299   {
0300     std::cout << PHWHERE << "funkyCaloStuff::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0301     return 0;
0302   }
0303 
0304   //////////////////////////////////////////
0305   // geometry for hot tower/cluster masking
0306   std::string towergeomnodename = "TOWERGEOM_CEMC";
0307   RawTowerGeomContainer* m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
0308   if (!m_geometry)
0309   {
0310     std::cout << Name() << "::"
0311               << "CreateNodeTree"
0312               << ": Could not find node " << towergeomnodename << std::endl;
0313     throw std::runtime_error("failed to find TOWERGEOM node in RawClusterDeadHotMask::CreateNodeTree");
0314   }
0315 
0316   /////////////////////////////
0317   // clusters
0318   RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0319   RawClusterContainer::ConstIterator clusterIter;
0320   RawClusterContainer::ConstIterator clusterIter2;
0321   int nClusCount = 0;
0322   for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0323   {
0324     RawCluster* recoCluster = clusterIter->second;
0325 
0326     CLHEP::Hep3Vector vertex(0, 0, vtx_z);
0327     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0328 
0329     float clus_pt = E_vec_cluster.perp();
0330     float clus_chisq = recoCluster->get_prob();
0331 
0332     if (clus_pt < nClus_ptCut)
0333     {
0334       continue;
0335     }
0336     if (clus_chisq < clus_chisq_cut)
0337     {
0338       continue;
0339     }
0340     nClusCount++;
0341   }
0342 
0343   h_nclusters->Fill(nClusCount);
0344 
0345   if (nClusCount > max_nClusCount)
0346   {
0347     return Fun4AllReturnCodes::EVENT_OK;
0348   }
0349 
0350   if (std::fabs(vtx_z) > vtx_z_cut && doVtxCut)
0351   {
0352     return Fun4AllReturnCodes::EVENT_OK;
0353   }
0354 
0355   h_event->Fill(0);
0356 
0357   float ptClusMax = 4;
0358   float pt1ClusCut = pt1BaseClusCut;
0359   float pt2ClusCut = pt2BaseClusCut;
0360 
0361   if (nClusCount > 30)
0362   {
0363     pt1ClusCut += NclusDeptFac * (nClusCount - 29) / 200.0;
0364     pt2ClusCut += NclusDeptFac * (nClusCount - 29) / 200.0;
0365   }
0366 
0367   float pi0ptcut = 1.22 * (pt1ClusCut + pt2ClusCut);
0368 
0369   for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0370   {
0371     RawCluster* recoCluster = clusterIter->second;
0372 
0373     CLHEP::Hep3Vector vertex(0, 0, vtx_z);
0374     CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0375 
0376     float clusE = E_vec_cluster.mag();
0377     float clus_eta = E_vec_cluster.pseudoRapidity();
0378     float clus_phi = E_vec_cluster.phi();
0379     float clus_pt = E_vec_cluster.perp();
0380     float clus_chisq = recoCluster->get_prob();
0381 
0382     if (clus_chisq < clus_chisq_cut)
0383     {
0384       continue;
0385     }
0386     h_clus_pt->Fill(clus_pt);
0387 
0388     unsigned int lt_eta = recoCluster->get_lead_tower().first;
0389     unsigned int lt_phi = recoCluster->get_lead_tower().second;
0390 
0391     h_etaphi_clus->Fill(clus_eta, clus_phi);
0392 
0393     TLorentzVector photon1;
0394     photon1.SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clusE);
0395 
0396     if (clus_pt < pt1ClusCut || clus_pt > ptClusMax)
0397     {
0398       continue;
0399     }
0400 
0401     for (clusterIter2 = clusterEnd.first; clusterIter2 != clusterIter; clusterIter2++)
0402     {
0403       RawCluster* recoCluster2 = clusterIter2->second;
0404 
0405       CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*recoCluster2, vertex);
0406 
0407       float clus2E = E_vec_cluster2.mag();
0408       float clus2_eta = E_vec_cluster2.pseudoRapidity();
0409       float clus2_phi = E_vec_cluster2.phi();
0410       float clus2_pt = E_vec_cluster2.perp();
0411       float clus2_chisq = recoCluster2->get_prob();
0412 
0413       if (clus2_pt < pt2ClusCut || clus_pt > ptClusMax)
0414       {
0415         continue;
0416       }
0417       if (clus2_chisq < clus_chisq_cut)
0418       {
0419         continue;
0420       }
0421 
0422       TLorentzVector photon2;
0423       photon2.SetPtEtaPhiE(clus2_pt, clus2_eta, clus2_phi, clus2E);
0424 
0425       if (std::fabs(clusE - clus2E) / (clusE + clus2E) > maxAlpha)
0426       {
0427         continue;
0428       }
0429 
0430       if (photon1.DeltaR(photon2) > maxDr)
0431       {
0432         continue;
0433       }
0434 
0435       TLorentzVector pi0 = photon1 + photon2;
0436       if (pi0.Pt() < pi0ptcut)
0437       {
0438         continue;
0439       }
0440 
0441       h_pt1->Fill(photon1.Pt());
0442       h_pt2->Fill(photon2.Pt());
0443 
0444       h_InvMass->Fill(pi0.M());
0445       if (clus2_pt < pt1ClusCut)
0446       {
0447         h_InvMass->Fill(pi0.M());
0448       }
0449 
0450       if (lt_eta > 95)
0451       {
0452         continue;
0453       }
0454       h_mass_eta_lt[lt_eta]->Fill(pi0.M());
0455 
0456       if (runTBTCompactMode)
0457       {
0458         h_ieta_iphi_invmass->Fill(lt_eta, lt_phi, pi0.M());
0459       }  // fill 3D hist for inv mass
0460 
0461       if (runTowByTow)
0462       {
0463         h_mass_tbt_lt[lt_eta][lt_phi]->Fill(pi0.M());
0464       }  // fill 1D inv mass hist for all towers
0465     }
0466   }  // clus1 loop
0467 
0468   return Fun4AllReturnCodes::EVENT_OK;
0469 }
0470 
0471 int pi0EtaByEta::End(PHCompositeNode* /*topNode*/)
0472 {
0473   outfile->cd();
0474 
0475   outfile->Write();
0476   outfile->Close();
0477   delete outfile;
0478   hm->dumpHistos(outfilename, "UPDATE");
0479 
0480   return 0;
0481 }
0482 
0483 TF1* pi0EtaByEta::fitHistogram(TH1* h)
0484 {
0485   TF1* f_sig_initial = new TF1("f_sig_initial", "[0]/[2]/2.5*exp(-0.5*((x-[1])/[2])^2)", 0.05, 0.25);
0486 
0487   f_sig_initial->SetParLimits(0, 0, 10);
0488   f_sig_initial->SetParLimits(1, 0.11, 0.18);
0489   f_sig_initial->SetParLimits(2, 0.007, 0.04);
0490   f_sig_initial->SetParameter(0, 0.1);
0491   f_sig_initial->SetParameter(1, 0.15);
0492   f_sig_initial->SetParameter(2, 0.015);
0493 
0494   h->Fit(f_sig_initial, "QN");
0495 
0496   TF1* fitFunc = new TF1("fitFunc", "[0]/[2]/2.5*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2 + [6]*x^3", h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
0497 
0498   fitFunc->SetParameter(0, f_sig_initial->GetParameter(0));
0499   fitFunc->SetParameter(1, f_sig_initial->GetParameter(1));
0500   fitFunc->SetParameter(2, f_sig_initial->GetParameter(2));
0501   fitFunc->SetParameter(3, 0.0);
0502   fitFunc->SetParameter(4, 0.0);
0503   fitFunc->SetParameter(5, 0.0);
0504   fitFunc->SetParameter(6, 0);
0505 
0506   fitFunc->SetParLimits(0, 0, 10);
0507   fitFunc->SetParLimits(1, 0.11, 0.2);
0508   fitFunc->SetParLimits(2, 0.007, 0.05);
0509   fitFunc->SetParLimits(3, -2, 1);
0510   fitFunc->SetParLimits(4, -1, 40);
0511   fitFunc->SetParLimits(5, -150, 50);
0512   // fitFunc->SetParLimits(6, -1,200 );
0513   fitFunc->SetParLimits(6, -1, 1);
0514 
0515   fitFunc->SetRange(0.06, 0.7);
0516 
0517   // Perform the fit
0518   h->Fit("fitFunc", "QN");
0519 
0520   return fitFunc;
0521 }
0522 
0523 void pi0EtaByEta::fitEtaSlices(const std::string& infile, const std::string& fitOutFile, const std::string& cdbFile)
0524 {
0525   TFile* fin = new TFile(infile.c_str());
0526   std::cout << "getting hists" << std::endl;
0527   TH1* h_peak_eta = new TH1F("h_peak_eta", "", 96, 0, 96);
0528   TH1* h_sigma_eta = new TH1F("h_sigma_eta", "", 96, 0, 96);
0529   TH1* h_p3_eta = new TH1F("h_p3_eta", "", 96, 0, 96);
0530   TH1* h_p4_eta = new TH1F("h_p4_eta", "", 96, 0, 96);
0531   TH1* h_p5_eta = new TH1F("h_p5_eta", "", 96, 0, 96);
0532   TH1* h_p6_eta = new TH1F("h_p6_eta", "", 96, 0, 96);
0533   TH1* h_p0_eta = new TH1F("h_p0_eta", "", 96, 0, 96);
0534   if (!fin)
0535   {
0536     std::cout << "pi0EtaByEta::fitEtaSlices null fin" << std::endl;
0537     exit(1);
0538   }
0539   TH1* h_M_eta[96];
0540   for (int i = 0; i < 96; i++)
0541   {
0542     std::string histoname = "h_mass_eta_lt" + std::to_string(i);
0543     fin->GetObject(histoname.c_str(), h_M_eta[i]);
0544     h_M_eta[i]->Scale(1. / h_M_eta[i]->Integral(), "width");
0545   }
0546 
0547   TF1* fitFunOut[96];
0548   for (int i = 0; i < 96; i++)
0549   {
0550     if (!h_M_eta[i])
0551     {
0552       std::cout << "pi0EtaByEta::fitEtaSlices null hist" << std::endl;
0553     }
0554     if (h_M_eta[i]->GetEntries() == 0)
0555     {
0556       continue;
0557     }
0558 
0559     fitFunOut[i] = fitHistogram(h_M_eta[i]);
0560     std::string funcname = "f_pi0_eta" + std::to_string(i);
0561     fitFunOut[i]->SetName(funcname.c_str());
0562     float mass_val_out = fitFunOut[i]->GetParameter(1);
0563     float mass_err_out = fitFunOut[i]->GetParError(1);
0564     h_peak_eta->SetBinContent(i + 1, mass_val_out);
0565     if (std::isnan(h_M_eta[i]->GetEntries()))
0566     {
0567       h_peak_eta->SetBinError(i + 1, 0);
0568       continue;
0569     }
0570     h_peak_eta->SetBinError(i + 1, mass_err_out);
0571     h_sigma_eta->SetBinContent(i + 1, fitFunOut[i]->GetParameter(2));
0572     h_sigma_eta->SetBinError(i + 1, fitFunOut[i]->GetParError(2));
0573     h_p3_eta->SetBinContent(i + 1, fitFunOut[i]->GetParameter(3));
0574     h_p3_eta->SetBinError(i + 1, fitFunOut[i]->GetParError(3));
0575     h_p4_eta->SetBinContent(i + 1, fitFunOut[i]->GetParameter(4));
0576     h_p4_eta->SetBinError(i + 1, fitFunOut[i]->GetParError(4));
0577     h_p5_eta->SetBinContent(i + 1, fitFunOut[i]->GetParameter(5));
0578     h_p5_eta->SetBinError(i + 1, fitFunOut[i]->GetParError(5));
0579     h_p6_eta->SetBinContent(i + 1, fitFunOut[i]->GetParameter(6));
0580     h_p6_eta->SetBinError(i + 1, fitFunOut[i]->GetParError(6));
0581     h_p0_eta->SetBinContent(i + 1, fitFunOut[i]->GetParameter(0));
0582     h_p0_eta->SetBinError(i + 1, fitFunOut[i]->GetParError(0));
0583   }
0584 
0585   CDBTTree* cdbttree1 = new CDBTTree(cdbFile);
0586   CDBTTree* cdbttree2 = new CDBTTree(cdbFile);
0587 
0588 
0589   float final_mass_target = target_pi0_mass;
0590 
0591   for (int i = 0; i < 96; i++)
0592   {
0593     for (int j = 0; j < 256; j++)
0594     {
0595       if (use_h_target_mass)
0596       {
0597         final_mass_target = h_target_mass->GetBinContent(i + 1);
0598       }
0599       float correction = final_mass_target / h_peak_eta->GetBinContent(i + 1);
0600       if (h_peak_eta->GetBinContent(i + 1) == 0)
0601       {
0602         correction = 0;
0603       }
0604       unsigned int key = TowerInfoDefs::encode_emcal(i, j);
0605       float val1 = cdbttree1->GetFloatValue(key, m_fieldname);
0606       cdbttree2->SetFloatValue(key, m_fieldname, val1 * correction);
0607     }
0608   }
0609 
0610   cdbttree2->Commit();
0611   cdbttree2->WriteCDBTTree();
0612   delete cdbttree2;
0613   delete cdbttree1;
0614 
0615   TFile* fit_out = new TFile(fitOutFile.c_str(), "recreate");
0616   fit_out->cd();
0617   for (auto& i : h_M_eta)
0618   {
0619     i->Write();
0620     delete i;
0621   }
0622   for (auto& i : fitFunOut)
0623   {
0624     i->Write();
0625     delete i;
0626   }
0627 
0628   h_p3_eta->Write();
0629   h_p4_eta->Write();
0630   h_p5_eta->Write();
0631   h_p6_eta->Write();
0632   h_p0_eta->Write();
0633   h_sigma_eta->Write();
0634   h_peak_eta->Write();
0635   fin->Close();
0636 
0637   std::cout << "finish fitting suc" << std::endl;
0638 
0639   return;
0640 }
0641 
0642 // for pi0 tbt fit
0643 void pi0EtaByEta::fitEtaPhiTowers(const std::string& infile, const std::string& fitOutFile, const std::string& cdbFile)
0644 {
0645   TFile* fin = new TFile(infile.c_str());
0646 
0647   std::cout << "getting hists" << std::endl;
0648 
0649   TH2* h_peak_tbt = new TH2F("h_peak_tbt", "", 96, 0, 96, 256, 0, 256);
0650   TH2* h_sigma_tbt = new TH2F("h_sigma_tbt", "", 96, 0, 96, 256, 0, 256);
0651   TH2* h_p3_tbt = new TH2F("h_p3_tbt", "", 96, 0, 96, 256, 0, 256);
0652   TH2* h_p4_tbt = new TH2F("h_p4_tbt", "", 96, 0, 96, 256, 0, 256);
0653   TH2* h_p5_tbt = new TH2F("h_p5_tbt", "", 96, 0, 96, 256, 0, 256);
0654   TH2* h_p6_tbt = new TH2F("h_p6_tbt", "", 96, 0, 96, 256, 0, 256);
0655   TH2* h_p0_tbt = new TH2F("h_p0_tbt", "", 96, 0, 96, 256, 0, 256);
0656   TH1* h_peak_tbt_1D = new TH1F("h_peak_tbt_1D", "", 24576, 0, 24576);
0657 
0658   if (!fin)
0659   {
0660     std::cout << "pi0EtaByEta::fitEtaPhiTowers null fin" << std::endl;
0661     exit(1);
0662   }
0663 
0664   TH1* h_M_tbt[96][256];
0665   TF1* fitFunOut[96][256];
0666 
0667   // creating output file
0668   TFile* fit_out = new TFile(fitOutFile.c_str(), "recreate");
0669 
0670   // we will load, fit and update one eta-bin a time and get into next eta-bin
0671   for (int ieta = 0; ieta < 96; ieta++)
0672   {
0673     std::cout << "fitting in etabin = " << ieta << std::endl;
0674     for (int j = 0; j < 256; j++)
0675     {
0676       std::string histoname = "h_mass_tbt_lt_" + std::to_string(ieta) + "_" + std::to_string(j);
0677       fin->GetObject(histoname.c_str(), h_M_tbt[ieta][j]);
0678 
0679       if (!h_M_tbt[ieta][j])
0680       {
0681         std::cout << "pi0EtaByEta::fitEtaPhiTowers null hist" << std::endl;
0682         gSystem->Exit(1);
0683         exit(1);
0684       }
0685  
0686       if (h_M_tbt[ieta][j]->GetEntries() == 0)
0687       {
0688         continue;
0689       }
0690 
0691       h_M_tbt[ieta][j]->Scale(1. / h_M_tbt[ieta][j]->Integral(), "width");
0692 
0693       fitFunOut[ieta][j] = fitHistogram(h_M_tbt[ieta][j]);
0694       std::string funcname = "f_pi0_tbt_" + std::to_string(ieta) + "_" + std::to_string(j);
0695       fitFunOut[ieta][j]->SetName(funcname.c_str());
0696 
0697       float mass_val_out = fitFunOut[ieta][j]->GetParameter(1);
0698       float mass_err_out = fitFunOut[ieta][j]->GetParError(1);
0699 
0700       if (std::isnan(h_M_tbt[ieta][j]->GetEntries()))
0701       {
0702         h_peak_tbt->SetBinContent(ieta + 1, j + 1, 0);
0703         h_peak_tbt->SetBinError(ieta + 1, j + 1, 0);
0704         continue;
0705       }
0706 
0707       h_peak_tbt->SetBinContent(ieta + 1, j + 1, mass_val_out);
0708       h_peak_tbt->SetBinError(ieta + 1, j + 1, mass_err_out);
0709 
0710       int towerID = TowerInfoDefs::decode_emcal(TowerInfoDefs::encode_emcal(ieta, j));
0711       h_peak_tbt_1D->SetBinContent(towerID + 1, mass_val_out);
0712       h_peak_tbt_1D->SetBinError(towerID + 1, mass_err_out);
0713 
0714       h_sigma_tbt->SetBinContent(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParameter(2));
0715       h_sigma_tbt->SetBinError(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParError(2));
0716 
0717       h_p3_tbt->SetBinContent(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParameter(3));
0718       h_p3_tbt->SetBinError(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParError(3));
0719 
0720       h_p4_tbt->SetBinContent(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParameter(4));
0721       h_p4_tbt->SetBinError(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParError(4));
0722 
0723       h_p5_tbt->SetBinContent(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParameter(5));
0724       h_p5_tbt->SetBinError(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParError(5));
0725 
0726       h_p6_tbt->SetBinContent(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParameter(6));
0727       h_p6_tbt->SetBinError(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParError(6));
0728 
0729       h_p0_tbt->SetBinContent(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParameter(0));
0730       h_p0_tbt->SetBinError(ieta + 1, j + 1, fitFunOut[ieta][j]->GetParError(0));
0731 
0732     }  // fitting for all phi-towers in eta-bin
0733 
0734     fit_out->cd();
0735 
0736     // writing all towers and fit function in given etabin
0737     for (auto& _hist : h_M_tbt[ieta])
0738     {
0739       _hist->Write();
0740       delete _hist;
0741     }
0742 
0743     for (auto& _hist : fitFunOut[ieta])
0744     {
0745       _hist->Write();
0746       delete _hist;
0747     }
0748 
0749   }  // closing eta-bin loop
0750 
0751   // CDBTtree
0752   CDBTTree* cdbttree1 = new CDBTTree(cdbFile);
0753   CDBTTree* cdbttree2 = new CDBTTree(cdbFile);
0754 
0755   float final_mass_target = target_pi0_mass;
0756 
0757   for (int i = 0; i < 96; i++)
0758   {
0759     for (int j = 0; j < 256; j++)
0760     {
0761       if (use_h_target_mass)
0762       {
0763         final_mass_target = h_target_mass->GetBinContent(i + 1, j + 1);
0764       }
0765       float correction = final_mass_target / h_peak_tbt->GetBinContent(i + 1, j + 1);
0766       if (h_peak_tbt->GetBinContent(i + 1, j + 1) == 0)
0767       {
0768         correction = 0;
0769       }
0770       if (h_sigma_tbt->GetBinContent(i + 1, j + 1) > 0.040)
0771       {
0772         correction = 1;
0773       }
0774       unsigned int key = TowerInfoDefs::encode_emcal(i, j);
0775       float val1 = cdbttree1->GetFloatValue(key, m_fieldname);
0776       cdbttree2->SetFloatValue(key, m_fieldname, val1 * correction);
0777     }
0778   }
0779 
0780   // writing out CDBTtree
0781   cdbttree2->Commit();
0782   cdbttree2->WriteCDBTTree();
0783 
0784   delete cdbttree2;
0785   delete cdbttree1;
0786 
0787   // wrriting other histograms in output file
0788   h_p3_tbt->Write();
0789   h_p4_tbt->Write();
0790   h_p5_tbt->Write();
0791   h_p6_tbt->Write();
0792   h_p0_tbt->Write();
0793   h_sigma_tbt->Write();
0794   h_peak_tbt->Write();
0795   h_peak_tbt_1D->Write();
0796 
0797   // closing input and output files
0798   fin->Close();
0799   fit_out->Close();
0800 
0801   std::cout << "finish fitting suc" << std::endl;
0802 
0803   return;
0804 }
0805 
0806 bool pi0EtaByEta::checkOutput(const std::string& file)
0807 {
0808   TFile* fin = new TFile(file.c_str());
0809   TH1* h_peak_eta{nullptr};
0810   fin->GetObject("h_peak_eta", h_peak_eta);
0811 
0812   float final_mass_target = target_pi0_mass;
0813 
0814   int numNotConv = 0;
0815 
0816   for (int i = 0; i < 96; i++)
0817   {
0818     if (use_h_target_mass)
0819     {
0820       final_mass_target = h_target_mass->GetBinContent(i + 1);
0821     }
0822     float corr = 1.0 - (final_mass_target / h_peak_eta->GetBinContent(i + 1));
0823     if (h_peak_eta->GetBinContent(i + 1) == 0)
0824     {
0825       corr = 0;
0826     }
0827     std::cout << "err " << corr << std::endl;
0828     if (std::fabs(corr) > convLev)
0829     {
0830       numNotConv++;
0831     }
0832   }
0833 
0834   delete fin;
0835 
0836   if (numNotConv <= 1)
0837   {
0838     return true;
0839   }
0840   return false;
0841 }
0842 
0843 void pi0EtaByEta::set_massTargetHistFile(const std::string& file)
0844 {
0845   TFile* fin = new TFile(file.c_str());
0846   if (fin)
0847   {
0848     fin->GetObject("h_massTargetHist", h_target_mass);
0849   }
0850   else
0851   {
0852     std::cout << "pi0EtaByEta::set_massTargetHistFile  No mass file found" << std::endl;
0853     return;
0854   }
0855   if (h_target_mass)
0856   {
0857     use_h_target_mass = true;
0858     std::cout << "using target mass histogram" << std::endl;
0859   }
0860   else
0861   {
0862     std::cout << "pi0EtaByEta::set_massTargetHistFile  No mass hist found" << std::endl;
0863   }
0864   return;
0865 }
0866 
0867 // this will split one 3D Hist into 24576 1D hist that contain inv mass distribution for all towers
0868 void pi0EtaByEta::Split3DHist(const std::string& infile, const std::string& out_file)
0869 {
0870   // First we will copy all the content of input file to output file
0871   std::error_code ec;
0872   std::filesystem::copy(infile, out_file, ec);
0873   if (ec)
0874   {
0875     std::cout << "Error copying file from " << infile << " to " << out_file << std::endl;
0876     return;
0877   }
0878 
0879   // Then, we will open the output file after update
0880   TFile* ofile = new TFile(out_file.c_str(), "UPDATE");
0881   if (!ofile)
0882   {
0883     std::cout << "Error opening file: " << out_file << std::endl;
0884     return;
0885   }
0886 
0887   // We will extract 3D hist from updated output file
0888   ofile->GetObject("h_ieta_iphi_invmass", h_ieta_iphi_invmass);
0889   if (!h_ieta_iphi_invmass)
0890   {
0891     std::cout << "Error: 3D Histogram not found in file: " << out_file << std::endl;
0892     ofile->Close();
0893     return;
0894   }
0895 
0896   // Loop over ieta and iphi ranges (x and y axis)
0897   for (int bineta = 1; bineta <= 96; ++bineta)
0898   {
0899     for (int binphi = 1; binphi <= 256; ++binphi)
0900     {
0901       // Loop over third axis in 3D Hist and then fill it into 1D hist
0902       for (int binz = 1; binz <= h_ieta_iphi_invmass->GetNbinsZ(); ++binz)
0903       {
0904         float content = h_ieta_iphi_invmass->GetBinContent(bineta, binphi, binz);
0905         h_mass_tbt_lt[bineta - 1][binphi - 1]->SetBinContent(binz, content);
0906       }
0907     }
0908   }
0909 
0910   // Now we will update all the histogram in the output file
0911   ofile->cd();
0912 
0913   for (auto& i : h_mass_tbt_lt)
0914   {
0915     for (auto& j : i)
0916     {
0917       j->Write();
0918       delete j;
0919     }
0920   }
0921 
0922   ofile->Close();
0923 
0924   std::cout << "Splitting 3D Hist Completed." << std::endl;
0925 
0926   return;
0927 }
[ 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!