Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-04-03 08:12:35 

0001 #include "wholeEventEECs.h"
0002 
0003 #include <format>
0004 
0005 
0006 wholeEventEECs::wholeEventEECs(const std::string &name)
0007   : SubsysReco(name)
0008 {
0009 }
0010 
0011 bool wholeEventEECs::isGoodTrigger(PHCompositeNode *topNode)
0012 {
0013     bool goodTrigger = false;
0014     TriggerAnalyzer *ta = new TriggerAnalyzer();
0015     ta->UseEmulator(false);
0016     ta->decodeTriggers(topNode);
0017     goodTrigger = ta->didTriggerFire("Jet 10 GeV");
0018     return goodTrigger;
0019 }
0020 
0021 bool wholeEventEECs::isGoodDijet()
0022 {
0023     if(vtx)
0024     {
0025         m_vtx_z = vtx->get_z();
0026         if(std::abs(m_vtx_z) > 60.0)
0027         {
0028             std::cout << "   vertex not in range" << std::endl;
0029             return false;
0030         };
0031     }
0032     else
0033     {
0034         m_vtx_z = 0.0;
0035     }
0036 
0037     if(jets->size() < 2)
0038     {
0039         std::cout << "   fewer than 2 jets" << std::endl;
0040         return false;
0041     }
0042 
0043     Jet *leadJet = nullptr;
0044     Jet *subleadJet = nullptr;
0045     float lead_pT = 0.0;
0046     float sublead_pT = 0.0;
0047     
0048     for(int i=0; i<(int)jets->size(); i++)
0049     {
0050         Jet *j = jets->get_jet(i);
0051         float pT = j->get_pt();
0052         if(pT > lead_pT)
0053         {
0054             if(leadJet)
0055             {
0056                 subleadJet = leadJet;
0057                 sublead_pT = lead_pT;
0058             }
0059             leadJet = j;
0060             lead_pT = pT;
0061         }
0062     }
0063     
0064     if(!leadJet || !subleadJet)
0065     {
0066         std::cout << "   bad lead or sublead jet" << std::endl;
0067         return false;
0068     }
0069 
0070     if(lead_pT < sublead_pT)
0071     {
0072         std::cout << "   lead pT < sub pT" << std::endl;
0073         return false;
0074     }
0075 
0076     if(lead_pT < 20.9 || sublead_pT < 9.4)
0077     {
0078         std::cout << "   lead or sub pT too low: lead < 20.9? " << lead_pT << "   sub < 9.4? " << sublead_pT << std::endl;
0079         return false;
0080     }
0081 
0082     if(std::abs(leadJet->get_eta()) > 0.7 || std::abs(subleadJet->get_eta()) > 0.7)
0083     {
0084         std::cout << "   lead or sub not in right eta range" << std::endl;
0085         return false;
0086     }
0087 
0088     float dPhi = subleadJet->get_phi() - leadJet->get_phi();
0089     if(dPhi > M_PI) dPhi -= 2*M_PI;
0090     if(dPhi < -M_PI) dPhi += 2*M_PI;
0091     if(std::abs(dPhi) < 0.75*M_PI)
0092     {
0093         std::cout << "   dPhi between sub and lead too small (should be greater than " << 0.75*M_PI << "): " << std::abs(dPhi) << std::endl;
0094         return false;
0095     }
0096 
0097     return true;
0098 }
0099 
0100 std::array<float, 3> wholeEventEECs::correct_for_vertex(std::array<float, 3> center)
0101 {
0102     std::array<float, 3> correctedTower;
0103 
0104     float z = center[2]*sinh(center[0]);
0105     if(m_vtx_z != 0.0)
0106     {
0107         z += m_vtx_z;
0108         float newEta = asinh(z/center[2]);
0109         if(!std::isnan(newEta) && !std::isinf(newEta))
0110         {
0111             correctedTower = {newEta, center[1], center[2]};
0112         }
0113     }
0114     else
0115     {
0116         correctedTower = {center[0], center[1], center[2]};
0117     }
0118 
0119     return correctedTower;
0120 
0121 }
0122 
0123 std::array<float, 5> wholeEventEECs::calc_observables(std::map<std::array<float, 3>, float>::iterator &iterA, std::map<std::array<float, 3>, float>::iterator &iterB)
0124 {
0125     float dEta = std::abs(iterA->first[0] - iterB->first[0]);
0126     float dPhi = iterA->first[1] - iterB->first[1];
0127     if(dPhi > M_PI) dPhi -= 2*M_PI;
0128     if(dPhi < -M_PI) dPhi += 2*M_PI;
0129     dPhi = std::abs(dPhi);
0130 
0131     float dR = sqrt(dEta*dEta + dPhi*dPhi);
0132 
0133 
0134     float pTA = iterA->second/cosh(iterA->first[0]);
0135     float pxA = pTA*cos(iterA->first[1]);
0136     float pyA = pTA*sin(iterA->first[1]);
0137     float pzA = pTA*sinh(iterA->first[0]);
0138 
0139     float pTB = iterB->second/cosh(iterB->first[0]);
0140     float pxB = pTB*cos(iterB->first[1]);
0141     float pyB = pTB*sin(iterB->first[1]);
0142     float pzB = pTB*sinh(iterB->first[0]);
0143 
0144     float cosTheta = (pxA*pxB + pyA*pyB + pzA*pzB)/(iterA->second * iterB->second);
0145 
0146     return std::array<float, 5> {dR, dEta, dPhi, std::acos(cosTheta), (((float)1.0)-cosTheta)/((float)2.0)};
0147 }
0148 
0149 
0150 int wholeEventEECs::InitRun(PHCompositeNode *topNode)
0151 {
0152 
0153 
0154 
0155     towerInfoContainers[0] = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_CEMC_RETOWER");
0156     towerInfoContainers[1] = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_HCALIN");
0157     towerInfoContainers[2] = findNode::getClass<TowerInfoContainer>(topNode,"TOWERINFO_CALIB_HCALOUT");
0158 
0159     if(!towerInfoContainers[0] || !towerInfoContainers[1] || !towerInfoContainers[2])
0160     {
0161         std::cout << "One or more TowerInfoContainers missing. Exiting" << std::endl;
0162         return Fun4AllReturnCodes::ABORTRUN;
0163     }
0164 
0165     emcal_geom = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC");
0166     geoms[0] = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0167     geoms[1] = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0168     geoms[2] = findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0169     if(!geoms[0] || !geoms[1] || !geoms[2])
0170     {
0171         std::cout << "One or more Tower Geometry Containers missing. Exiting" << std::endl;
0172         return Fun4AllReturnCodes::ABORTRUN;
0173     }
0174 
0175     emcal_geom->set_calorimeter_id(RawTowerDefs::CEMC);
0176     geoms[0]->set_calorimeter_id(RawTowerDefs::HCALIN);
0177     geoms[1]->set_calorimeter_id(RawTowerDefs::HCALIN);
0178     geoms[2]->set_calorimeter_id(RawTowerDefs::HCALOUT);
0179     
0180 
0181     vtxMap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0182     if(!vtxMap)
0183     {
0184         vtxMap = findNode::getClass<GlobalVertexMap>(topNode, "MbdVertexMap");
0185         if(!vtxMap)
0186         {
0187             std::cout << "No global vertex map found. Exiting" << std::endl;
0188             return Fun4AllReturnCodes::ABORTRUN;
0189         }
0190     }
0191 
0192     //jets = findNode::getClass<JetContainer>(topNode, "AntiKt_TowerInfo_r04");
0193     jets = findNode::getClass<JetContainer>(topNode, "AntiKt_unsubtracted_r04");
0194     if(!jets)
0195     {
0196         std::cout << "AntiKt_TowerInfo_r04 jet node missing. Exiting" << std::endl;
0197         return Fun4AllReturnCodes::ABORTRUN;
0198     }
0199 
0200     m_outputFile = new TFile(m_outputName.c_str(),"RECREATE");
0201 
0202     nEvents = new TH1D("nEvents","Number of dijet events",1,0.5,1.5);
0203 
0204     for(int i=0; i<5; i++)
0205     {
0206         for(int j=0; j<4; j++)
0207         {
0208             EEC[i][j] = new TH1D(std::format("EEC_{}_{}", obsNames[i], caloNames[j]).c_str(), std::format("Whole Event EEC {};{}", caloNames[j], obsSymbols[i]).c_str(),50,0.0,obsMax[i]);
0209             EEC_corr[i][j] = new TH1D(std::format("EEC_corr_{}_{}", obsNames[i], caloNames[j]).c_str(), std::format("Vertex Corrected Whole Event EEC {};{}", caloNames[j], obsSymbols[i]).c_str(),50,0.0,obsMax[i]);
0210         }
0211 
0212         /*
0213         EEC_RL[i] = new TH1D(std::format("EEC_RL_{}", caloNames[i]).c_str(), std::format("Whole Event EEC {};#DeltaR", caloNames[i]).c_str(),50,0.0,sqrt(2.2*2.2 + M_PI*M_PI));
0214         EEC_eta[i] = new TH1D(std::format("EEC_eta_{}", caloNames[i]).c_str(), std::format("Whole Event EEC {};#Delta#eta", caloNames[i]).c_str(),50,0.0,2.2);
0215         EEC_phi[i] = new TH1D(std::format("EEC_phi_{}", caloNames[i]).c_str(), std::format("Whole Event EEC {};#Delta#phi", caloNames[i]).c_str(),50,0.0,M_PI);
0216         EEC_theta[i] = new TH1D(std::format("EEC_theta_{}", caloNames[i]).c_str(), std::format("Whole Event EEC {};#theta", caloNames[i]).c_str(),50,0.0,M_PI);
0217         EEC_z[i] = new TH1D(std::format("EEC_z_{}", caloNames[i]).c_str(), std::format("Whole Event EEC {};z=(1-cos#theta)/2", caloNames[i]).c_str(),50,0.0,1.0);
0218     
0219         EEC_RL_corr[i] = new TH1D(std::format("EEC_RL_corr_{}", caloNames[i]).c_str(), std::format("Vertex Corrected Whole Event EEC {};#DeltaR", caloNames[i]).c_str(),50,0.0,sqrt(2.2*2.2 + M_PI*M_PI));
0220         EEC_eta_corr[i] = new TH1D(std::format("EEC_eta_corr_{}", caloNames[i]).c_str(), std::format("Vertex Corrected Whole Event EEC {};#Delta#eta", caloNames[i]).c_str(),50,0.0,2.2);
0221         EEC_phi_corr[i] = new TH1D(std::format("EEC_phi_corr_{}", caloNames[i]).c_str(), std::format("Vertex Corrected Whole Event EEC {};#Delta#phi", caloNames[i]).c_str(),50,0.0,M_PI);
0222         EEC_theta_corr[i] = new TH1D(std::format("EEC_theta_corr_{}", caloNames[i]).c_str(), std::format("Vertex Corrected Whole Event EEC {};#theta", caloNames[i]).c_str(),50,0.0,M_PI);
0223         EEC_z_corr[i] = new TH1D(std::format("EEC_z_corr_{}", caloNames[i]).c_str(), std::format("Vertex Corrected Whole Event EEC {};z=(1-cos#theta)/2", caloNames[i]).c_str(),50,0.0,1.0);
0224         */
0225     }
0226 
0227 
0228     return Fun4AllReturnCodes::EVENT_OK;
0229 
0230 }
0231 
0232 int wholeEventEECs::process_event(PHCompositeNode *topNode)
0233 {
0234     for(int i=0; i<4; i++)
0235     {
0236         towers[i].clear();
0237         towersCorrected[i].clear();
0238     }
0239 
0240     /*
0241     EMCalTowers.clear();
0242     IHCalTowers.clear();
0243     OHCalTowers.clear();
0244     AllTowers.clear();
0245 
0246     EMCalCorrectedTowers.clear();
0247     IHCalCorrectedTowers.clear();
0248     OHCalCorrectedTowers.clear();
0249     AllCorrectedTowers.clear();
0250     */
0251 
0252     m_eventIndex++;
0253 
0254 
0255     std::cout << "working on event " << m_eventIndex << ". There have been " << m_nGoodEvents << " good events so far." << std::endl;
0256 
0257     bool goodTrigger = isGoodTrigger(topNode);
0258     if(!goodTrigger)
0259     {
0260         std::cout << "   trigger not found" << std::endl;
0261         return Fun4AllReturnCodes::EVENT_OK;
0262     }
0263 
0264     if(vtxMap->empty())
0265     {
0266         std::cout << "      empty vertex map, setting vz to origin" << std::endl;
0267         vtx = nullptr;
0268     }
0269     else
0270     {
0271         for(auto vItr : *vtxMap)
0272         {
0273             if(vItr.first == 0 || vItr.first == GlobalVertex::VTXTYPE::MBD || vItr.first == GlobalVertex::VTXTYPE::SVTX || vItr.first == GlobalVertex::VTXTYPE::SVTX_MBD)
0274             {
0275                 vtx = vItr.second;
0276             }
0277         }
0278     }
0279 
0280     bool goodDijet = isGoodDijet();
0281     if(!goodDijet)
0282     {
0283         std::cout << "   good dijet not found" << std::endl;
0284         return Fun4AllReturnCodes::EVENT_OK;
0285     }
0286 
0287     m_nGoodEvents++;
0288     nEvents->Fill(1);
0289 
0290     for(int calo = 0; calo < 3; calo++)
0291     {
0292         for(int i = 0; i < (int) towerInfoContainers[calo]->size(); i++)
0293         {
0294             if(!towerInfoContainers[calo]->get_tower_at_channel(i)->get_isGood()) continue;
0295             auto key = towerInfoContainers[calo]->encode_key(i);
0296             auto tower = towerInfoContainers[calo]->get_tower_at_channel(i);
0297             int phiBin = towerInfoContainers[calo]->getTowerPhiBin(key);
0298             int etaBin = towerInfoContainers[calo]->getTowerEtaBin(key);
0299             float phiCenter = geoms[calo]->get_phicenter(phiBin);
0300             float etaCenter = geoms[calo]->get_etacenter(etaBin);
0301             float r = geoms[calo]->get_radius();
0302             if(calo == 0) r = emcal_geom->get_radius();
0303             std::array<float, 3> center {etaCenter, phiCenter, r};
0304             towers[calo].insert(std::make_pair(center, tower->get_energy()));
0305             if(towers[3].find(center) != towers[3].end()) towers[3].at(center) += tower->get_energy();
0306             else towers[3][center] = tower->get_energy();
0307 
0308             std::array<float, 3> newCenter = correct_for_vertex(center);
0309             towersCorrected[calo].insert(std::make_pair(newCenter, tower->get_energy()));
0310             if(towersCorrected[3].find(newCenter) != towersCorrected[3].end()) towersCorrected[3].at(newCenter) += tower->get_energy();
0311             else towersCorrected[3][newCenter] = tower->get_energy();
0312         }
0313     }
0314 
0315     /*
0316     for(int i=0; i<(int) emcalTowerInfoContainer->size(); i++)
0317     {
0318         if(!emcalTowerInfoContainer->get_tower_at_channel(i)->get_isGood()) continue;
0319         
0320         auto key = emcalTowerInfoContainer->encode_key(i);
0321         auto tower = emcalTowerInfoContainer->get_tower_at_channel(i);
0322         int phiBin = emcalTowerInfoContainer->getTowerPhiBin(key);
0323         int etaBin = emcalTowerInfoContainer->getTowerEtaBin(key);
0324         float phiCenter = ihcal_geom->get_phicenter(phiBin);
0325         float etaCenter = ihcal_geom->get_etacenter(etaBin);
0326         float r = emcal_geom->get_radius();
0327         std::array<float, 3> center {etaCenter, phiCenter, r};
0328         EMCalTowers.insert(std::make_pair(center, tower->get_energy()));
0329         if(AllTowers.find(center) != AllTowers.end()) AllTowers.at(center) += tower->get_energy();
0330         else AllTowers[center] = tower->get_energy();
0331 
0332         std::array<float, 3> newCenter = correct_for_vertex(center);
0333         EMCalCorrectedTowers.insert(std::make_pair(newCenter, tower->get_energy()));
0334         if(AllCorrectedTowers.find(newCenter) != AllCorrectedTowers.end()) AllCorrectedTowers.at(newCenter) += tower->get_energy();
0335         else AllCorrectedTowers[newCenter] = tower->get_energy();
0336     }
0337 
0338     for(int i=0; i<(int) ihcalTowerInfoContainer->size(); i++)
0339     {
0340         if(!ihcalTowerInfoContainer->get_tower_at_channel(i)->get_isGood()) continue;
0341         ihcal_geom->set_calorimeter_id(RawTowerDefs::HCALIN);
0342         auto key = ihcalTowerInfoContainer->encode_key(i);
0343         auto tower = ihcalTowerInfoContainer->get_tower_at_channel(i);
0344         int phiBin = ihcalTowerInfoContainer->getTowerPhiBin(key);
0345         int etaBin = ihcalTowerInfoContainer->getTowerEtaBin(key);
0346         float phiCenter = ihcal_geom->get_phicenter(phiBin);
0347         float etaCenter = ihcal_geom->get_etacenter(etaBin);
0348         float r = ihcal_geom->get_radius();
0349         std::array<float, 3> center {etaCenter, phiCenter, r};
0350         IHCalTowers.insert(std::make_pair(center, tower->get_energy()));
0351         if(AllTowers.find(center) != AllTowers.end()) AllTowers.at(center) += tower->get_energy();
0352         else AllTowers[center] = tower->get_energy();
0353 
0354         std::array<float, 3> newCenter = correct_for_vertex(center);
0355         IHCalCorrectedTowers.insert(std::make_pair(newCenter, tower->get_energy()));
0356         if(AllCorrectedTowers.find(newCenter) != AllCorrectedTowers.end()) AllCorrectedTowers.at(newCenter) += tower->get_energy();
0357         else AllCorrectedTowers[newCenter] = tower->get_energy();
0358     }
0359 
0360     for(int i=0; i<(int) ohcalTowerInfoContainer->size(); i++)
0361     {
0362         if(!ohcalTowerInfoContainer->get_tower_at_channel(i)->get_isGood()) continue;
0363         ohcal_geom->set_calorimeter_id(RawTowerDefs::HCALOUT);
0364         auto key = ohcalTowerInfoContainer->encode_key(i);
0365         auto tower = ohcalTowerInfoContainer->get_tower_at_channel(i);
0366         int phiBin = ohcalTowerInfoContainer->getTowerPhiBin(key);
0367         int etaBin = ohcalTowerInfoContainer->getTowerEtaBin(key);
0368         float phiCenter = ohcal_geom->get_phicenter(phiBin);
0369         float etaCenter = ohcal_geom->get_etacenter(etaBin);
0370         float r = ohcal_geom->get_radius();
0371         std::array<float, 3> center {etaCenter, phiCenter, r};
0372         OHCalTowers.insert(std::make_pair(center, tower->get_energy()));
0373         if(AllTowers.find(center) != AllTowers.end()) AllTowers.at(center) += tower->get_energy();
0374         else AllTowers[center] = tower->get_energy();
0375 
0376         std::array<float, 3> newCenter = correct_for_vertex(center);
0377         OHCalCorrectedTowers.insert(std::make_pair(newCenter, tower->get_energy()));
0378         if(AllCorrectedTowers.find(newCenter) != AllCorrectedTowers.end()) AllCorrectedTowers.at(newCenter) += tower->get_energy();
0379         else AllCorrectedTowers[newCenter] = tower->get_energy();
0380     }
0381     */
0382 
0383     for(int calo = 0; calo < 4; calo++)
0384     {
0385         //original calo bins
0386         for(auto it1 = towers[calo].begin(); it1 != towers[calo].end(); ++it1)
0387         {
0388             float pT1 = it1->second / cosh(it1->first[0]);
0389             if(pT1 < 0.3) continue;
0390             for(auto it2 = std::next(it1); it2 != towers[calo].end(); ++it2)
0391             {
0392                 float pT2 = it2->second / cosh(it2->first[0]);
0393                 if(pT2 < 0.3) continue;
0394 
0395                 std::array<float, 5> observables = calc_observables(it1, it2);
0396                 for(int i=0; i<5; i++)
0397                 {
0398                     EEC[i][calo]->Fill(observables[i], pT1 * pT2);
0399                 }
0400             }
0401         }
0402 
0403         //shiftex vertex
0404         for(auto it1 = towersCorrected[calo].begin(); it1 != towersCorrected[calo].end(); ++it1)
0405         {
0406             float pT1 = it1->second / cosh(it1->first[0]);
0407             if(pT1 < 0.3) continue;
0408             for(auto it2 = std::next(it1); it2 != towersCorrected[calo].end(); ++it2)
0409             {
0410                 float pT2 = it2->second / cosh(it2->first[0]);
0411                 if(pT2 < 0.3) continue;
0412 
0413                 std::array<float, 5> observables = calc_observables(it1, it2);
0414                 for(int i=0; i<5; i++)
0415                 {
0416                     EEC_corr[i][calo]->Fill(observables[i], pT1 * pT2);
0417                 }
0418             }
0419         }
0420     }
0421 
0422     
0423 
0424     /*
0425     //no correcting for vertex position
0426     for(auto it1 = EMCalTowers.begin(); it1 != EMCalTowers.end(); ++it1)
0427     {
0428         float pT1 = it1->second / cosh(it1->first[0]);
0429         if(pT1 < 0.3) continue;
0430         for(auto it2 = std::next(it1); it2 != EMCalTowers.end(); ++it2)
0431         {
0432             float pT2 = it2->second / cosh(it2->first[0]);
0433             if(pT2 < 0.3) continue;
0434 
0435 
0436             std::array<float, 5> observables = calc_observables(it1, it2);
0437             for(int i=0; i<5; i++)
0438             {
0439                 EEC[i][0]->Fill(observables[i], pT1 * pT2);
0440             }
0441         }
0442     }
0443 
0444     for(auto it1 = IHCalTowers.begin(); it1 != IHCalTowers.end(); ++it1)
0445     {
0446         float pT1 = it1->second / cosh(it1->first[0]);
0447         if(pT1 < 0.3) continue; 
0448         for(auto it2 = std::next(it1); it2 != IHCalTowers.end(); ++it2)
0449         {
0450             float pT2 = it2->second / cosh(it2->first[0]);
0451             if(pT2 < 0.3) continue;
0452 
0453 
0454             std::array<float, 5> observables = calc_observables(it1, it2);
0455             for(int i=0; i<5; i++)
0456             {
0457                 EEC[i][1]->Fill(observables[i], pT1 * pT2);
0458             }
0459         }
0460     }
0461 
0462     for(auto it1 = OHCalTowers.begin(); it1 != OHCalTowers.end(); ++it1)
0463     {
0464         float pT1 = it1->second / cosh(it1->first[0]);
0465         if(pT1 < 0.3) continue;
0466         for(auto it2 = std::next(it1); it2 != OHCalTowers.end(); ++it2)
0467         {
0468             float pT2 = it2->second / cosh(it2->first[0]);
0469             if(pT2 < 0.3) continue;
0470 
0471 
0472             std::array<float, 5> observables = calc_observables(it1, it2);
0473             for(int i=0; i<5; i++)
0474             {
0475                 EEC[i][2]->Fill(observables[i], pT1 * pT2);
0476             }
0477         }
0478     }
0479 
0480     for(auto it1 = AllTowers.begin(); it1 != AllTowers.end(); ++it1)
0481     {
0482         float pT1 = it1->second / cosh(it1->first[0]);
0483         if(pT1 < 0.3) continue;
0484         for(auto it2 = std::next(it1); it2 != AllTowers.end(); ++it2)
0485         {
0486             float pT2 = it2->second / cosh(it2->first[0]);
0487             if(pT2 < 0.3) continue;
0488 
0489 
0490             std::array<float, 5> observables = calc_observables(it1, it2);
0491             for(int i=0; i<5; i++)
0492             {
0493                 EEC[i][3]->Fill(observables[i], pT1 * pT2);
0494             }
0495         }
0496     }
0497 
0498     //correcting for vertex position
0499     for(auto it1 = EMCalCorrectedTowers.begin(); it1 != EMCalCorrectedTowers.end(); ++it1)
0500     {
0501         float pT1 = it1->second / cosh(it1->first[0]);
0502         if(pT1 < 0.3) continue;
0503         for(auto it2 = std::next(it1); it2 != EMCalCorrectedTowers.end(); ++it2)
0504         {
0505             float pT2 = it2->second / cosh(it2->first[0]);
0506             if(pT2 < 0.3) continue;
0507 
0508 
0509             std::array<float, 5> observables = calc_observables(it1, it2);
0510             for(int i=0; i<5; i++)
0511             {
0512                 EEC_corr[i][0]->Fill(observables[i], pT1 * pT2);
0513             }
0514         }
0515     }
0516 
0517     for(auto it1 = IHCalCorrectedTowers.begin(); it1 != IHCalCorrectedTowers.end(); ++it1)
0518     {
0519         float pT1 = it1->second / cosh(it1->first[0]);
0520         if(pT1 < 0.3) continue;
0521         for(auto it2 = std::next(it1); it2 != IHCalCorrectedTowers.end(); ++it2)
0522         {
0523             float pT2 = it2->second / cosh(it2->first[0]);
0524             if(pT2 < 0.3) continue;
0525 
0526 
0527             std::array<float, 5> observables = calc_observables(it1, it2);
0528             for(int i=0; i<5; i++)
0529             {
0530                 EEC_corr[i][1]->Fill(observables[i], pT1 * pT2);
0531             }
0532         }
0533     }
0534 
0535     for(auto it1 = OHCalCorrectedTowers.begin(); it1 != OHCalCorrectedTowers.end(); ++it1)
0536     {
0537         float pT1 = it1->second / cosh(it1->first[0]);
0538         if(pT1 < 0.3) continue;
0539         for(auto it2 = std::next(it1); it2 != OHCalCorrectedTowers.end(); ++it2)
0540         {
0541             float pT2 = it2->second / cosh(it2->first[0]);
0542             if(pT2 < 0.3) continue;
0543 
0544 
0545             std::array<float, 5> observables = calc_observables(it1, it2);
0546             for(int i=0; i<5; i++)
0547             {
0548                 EEC_corr[i][2]->Fill(observables[i], pT1 * pT2);
0549             }
0550         }
0551     }
0552 
0553     for(auto it1 = AllCorrectedTowers.begin(); it1 != AllCorrectedTowers.end(); ++it1)
0554     {
0555         float pT1 = it1->second / cosh(it1->first[0]);
0556         if(pT1 < 0.3) continue;
0557         for(auto it2 = std::next(it1); it2 != AllCorrectedTowers.end(); ++it2)
0558         {
0559             float pT2 = it2->second / cosh(it2->first[0]);
0560             if(pT2 < 0.3) continue;
0561 
0562 
0563             std::array<float, 5> observables = calc_observables(it1, it2);
0564             for(int i=0; i<5; i++)
0565             {
0566                 EEC_corr[i][3]->Fill(observables[i], pT1 * pT2);
0567             }
0568         }
0569     }
0570     */
0571 
0572     return Fun4AllReturnCodes::EVENT_OK;
0573 }
0574 
0575 int wholeEventEECs::End(PHCompositeNode */*topNode*/)
0576 {
0577     m_outputFile->cd();
0578     nEvents->Write();
0579     for(int i=0; i<5; i++)
0580     {
0581         for(int j=0; j<4; j++)
0582         {
0583             EEC[i][j]->Write();
0584             EEC_corr[i][j]->Write();
0585         }
0586     }
0587     m_outputFile->Close();
0588 
0589     return Fun4AllReturnCodes::EVENT_OK;
0590 }
[ 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!