Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:13:26 

0001 #include "LargeRLENC_LEDPedestalScan.h"
0002 // Just run a version of the LargeRL enc that only calls the processing and no cuts
0003 LEDPedestalScan::LEDPedestalScan(const int n_run, const int n_segment, const bool doE2C, const bool act_loc, const std::string& name)
0004 {
0005     if(n_run > 1000) 
0006         encCalc=new LargeRLENC(n_run, n_segment, 1.0, false, true);
0007     else 
0008         encCalc=new LargeRLENC(n_run, n_segment, 1.0, false, false);        
0009     this->run=n_run;
0010     this->segement=n_segment;
0011     this->local=act_loc;
0012     this->run_e2c=doE2C;
0013     jet_loc=new TH2F("jet_loc", "Fake Jet locations all jets; #eta; #varphi; #N_{hits}", 24, -1.1, 1.1, 64, 0, 2*PI); 
0014     jet_loc_center=new TH2F("jet_loc_center", "Fake Jet center all jets; #eta; #varphi; #N_{hits}", 24, -1.1, 1.1, 64, 0, 2*PI); 
0015     lead_jet_loc=new TH2F("lead_jet_loc", "Fake Jet locations leading jets; #eta; #varphi; #N_{hits}", 24, -1.1, 1.1, 64, 0, 2*PI); 
0016     lead_jet_loc_center=new TH2F("lead_jet_loc_center", "Fake Jet center leading jets; #eta; #varphi; #N_{hits}", 24, -1.1, 1.1, 64, 0, 2*PI); 
0017     this->energy_thresholds.push_back( encCalc->Thresholds);
0018     for(int i= 0; i<(int) energy_thresholds.size(); i++){   
0019         ohcal_energy_i.push_back(new TH2F(Form("ohcal_energy_%d_MeV",(int)(energy_thresholds[i][3]*1000)), Form("OHCAL energy, E_{tow} #geq %f MeV; #eta; #varphi; E [GeV]",(float)(energy_thresholds[i][3]*1000)), 24, -1.1, 1.1, 64, 0, 2*PI));
0020         ihcal_energy_i.push_back(new TH2F(Form("ihcal_energy_%d_MeV",(int)(energy_thresholds[i][2]*1000)), Form("IHCAL energy, E_{tow} #geq %f MeV; #eta; #varphi; E [GeV]",(float)(energy_thresholds[i][2]*1000)), 24, -1.1, 1.1, 64, 0, 2*PI));
0021         emcal_energy_i.push_back(new TH2F(Form("emcal_energy_%d_MeV",(int)(energy_thresholds[i][1]*1000)), Form("EMCAL energy, E_{tow} #geq %f MeV; #eta; #varphi; E [GeV]",(float)(energy_thresholds[i][1]*1000)), 24, -1.1, 1.1, 64, 0, 2*PI));
0022         ohcal_energy_node.push_back(new TH2F(Form("ohcal_energy_node_%d_MeV",(int)(energy_thresholds[i][3]*1000)), Form("OHCAL energy, E_{tow} #geq %f MeV; #eta; #varphi; E [GeV]",(float)(energy_thresholds[i][3]*1000.)), 24, -1.1, 1.1, 64, 0, 2*PI));
0023         ihcal_energy_node.push_back(new TH2F(Form("ihcal_energy_node_%d_MeV",(int)(energy_thresholds[i][2]*1000)), Form("IHCAL energy, E_{tow} #geq %f MeV; #eta; #varphi; E [GeV]", (float)(energy_thresholds[i][2]*1000.)),24, -1.1, 1.1, 64, 0, 2*PI));
0024         emcal_energy_node.push_back(new TH2F(Form("emcal_energy_node_%d_MeV",(int)( energy_thresholds[i][1]*1000)), Form("EMCAL energy, E_{tow} #geq %f MeV; #eta; #varphi; E [GeV]", (float)(energy_thresholds[i][1]*1000.)), 24, -1.1, 1.1, 64, 0, 2*PI));
0025     }
0026     ohcal_hit_plot=new TH1F("ohc_hit", "Outer HCAL total hits above threshold; Tower Threshold [MeV]; Hits (E_{tow} > E_{thresh}", 5000, -0.5, 4999.5);
0027     ihcal_hit_plot=new TH1F("ihc_hit", "Inner HCAL total hits above threshold; Tower Threshold [MeV]; Hits (E_{tow} > E_{thresh}", 5000, -0.5, 4999.5);
0028     emcal_hit_plot=new TH1F("emc_hit", "Electormagnetic Calorimeter total hits above threshold; Tower Threshold [MeV]; Hits (E_{tow} > E_{thresh}", 5000, -0.5, 4999.5);
0029     allcal_hit_plot=new TH1F("amc_hit", "All Calorimeter total hits above threshold (Hcal Tower binning); Tower Threshold [MeV]; Hits (E_{tow} > E_{thresh}", 5000, -0.5, 4999.5);
0030     ohcal_energy=new TH1F("ohcal_energy", "Outer HCal Tower Energy; E_{tower} [MeV]; N_{towers}", 2000, -1000.5, 999.5);
0031     ihcal_energy=new TH1F("ihcal_energy", "Inner HCal Tower Energy; E_{tower} [MeV]; N_{towers}", 2000, -1000.5, 999.5);
0032     emcal_energy=new TH1F("emcal_energy", "Electormagnetic Calorimeter Tower Energy; E_{tower} [MeV]; N_{towers}", 1000, -1000.5, 999.5);
0033     allcal_energy=new TH1F("allcal_energy", "All Calorimeter Tower Energy (hcal binning); E_{tower} [MeV]; N_{towers}", 5000, -0.5, 4999.5);
0034     ohcal_energy_total=new TH1F("ohcal_total_E", "Outer HCal Total energy; E_{Calo} [GeV]; N_{event}", 100, -0.5, 99.5);
0035     ihcal_energy_total=new TH1F("ihcal_total_E", "Inner HCal Total Energy; E_{Calo} [GeV]; N_{event}", 100, -0.5, 99.5);
0036     emcal_energy_total=new TH1F("emcal_total_E", "Electormagnetic Calorimeter Total Energy; E_{Cal} [GeV]; N_{event}", 1000, -0.5, 999.5);
0037     allcal_energy_total=new TH1F("allcal_total_E", "All Calorimeter Total Energy; E_{Cal} [GeV]; N_{event}", 1000, -0.5, 999.5);
0038     
0039     //ohcal_energy->SetCanExtend(TH1::kAllAxes);
0040 //  ihcal_energy->SetCanExtend(TH1::kAllAxes);
0041     //emcal_energy->SetCanExtend(TH1::kAllAxes);
0042 //  ohcal_energy_total->SetCanExtend(TH1::kAllAxes);
0043 //  ihcal_energy_total->SetCanExtend(TH1::kAllAxes);
0044 //  emcal_energy_total->SetCanExtend(TH1::kAllAxes);
0045     n_evts=0;
0046 }
0047 
0048 int LEDPedestalScan::process_event(PHCompositeNode* topNode)
0049 {
0050     n_evts++;
0051     float emcal_e=0, ohcal_e=0, ihcal_e=0;
0052     std::cout<<"Running Event " <<n_evts<<std::endl;
0053     std::array<float, 3> vertex {0., 0., 0.};
0054 
0055     if(run > 1000 && !local){
0056         bool did_trig=encCalc->triggerCut(false, topNode);
0057         auto jets = findNode::getClass<JetContainerv1>(topNode, "AntiKt_unsubtracted_r04");
0058         float max_jet_pt=0;
0059         for(auto j:*jets) if(j->get_pt() >= max_jet_pt) max_jet_pt=j->get_pt(); 
0060         if(!did_trig || max_jet_pt > 30. || max_jet_pt < 10.  ) return Fun4AllReturnCodes::ABORTEVENT;
0061         else std::cout<<"Found a Jet 10 trigger and jet with pt " <<max_jet_pt <<std::endl;
0062     }
0063     std::map<std::array<float, 3>, float> emcal_towers, ihcal_towers, ohcal_towers;
0064     std::string ohcal_energy_towers="TOWERINFO_CALIB_HCALOUT", ihcal_energy_towers="TOWERINFO_CALIB_HCALIN", emcal_energy_towers="TOWERINFO_CALIB_CEMC";
0065 //  std::string ohcal_energy_towers="TOWERS_HCALOUT", ihcal_energy_towers="TOWERS_HCALIN", emcal_energy_towers="TOWERS_CEMC";
0066     auto emcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC"   );
0067     auto emcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, emcal_energy_towers      );
0068     auto ihcal_geom= findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0069     auto ihcal_tower_energy= findNode::getClass<TowerInfoContainer>(topNode, ihcal_energy_towers     );
0070     auto ohcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0071     auto ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, ohcal_energy_towers      );
0072     if(n_evts==1) encCalc->MakeEMCALRetowerMap(emcal_geom,emcal_tower_energy, ohcal_geom, ihcal_tower_energy);
0073     std::map<std::pair<float, float>, float> allcal_tows;
0074     for(int n=0; n<(int) ohcal_tower_energy->size(); n++){
0075         ohcal_geom->set_calorimeter_id(RawTowerDefs::HCALOUT);
0076         auto key=ohcal_tower_energy->encode_key(n);
0077         int phibin=ohcal_tower_energy->getTowerPhiBin(key);
0078         int etabin=ohcal_tower_energy->getTowerEtaBin(key);
0079         float phicenter=ohcal_geom->get_phicenter(phibin);
0080         float etacenter=ohcal_geom->get_etacenter(etabin);
0081         allcal_tows[std::make_pair(etacenter, phicenter)]=0.;
0082     }
0083     for(int n=0; n<(int) emcal_tower_energy->size(); n++){
0084         emcal_geom->set_calorimeter_id(RawTowerDefs::CEMC);
0085         auto key=emcal_tower_energy->encode_key(n);
0086         auto tower=emcal_tower_energy->get_tower_at_channel(n);
0087         int phibin=emcal_tower_energy->getTowerPhiBin(key);
0088         int etabin=emcal_tower_energy->getTowerEtaBin(key);
0089         float phicenter=emcal_geom->get_phicenter(phibin);
0090         float etacenter=emcal_geom->get_etacenter(etabin);
0091         float te=tower->get_energy();
0092         if(this->local){
0093             te=te-tower->get_pedestal();
0094         }
0095         for(int i=0; i<(int)energy_thresholds.size(); i++){
0096             if(te > energy_thresholds[i][1]) 
0097                 emcal_energy_node[i]->Fill(etacenter, phicenter, te);
0098             }
0099         if(tower->get_energy() > energy_thresholds[0][1])
0100             encCalc->addTower(n, emcal_tower_energy, emcal_geom, &emcal_towers, RawTowerDefs::CEMC);
0101         std::pair<float, float> hcal_tow=encCalc->emcal_lookup_table[n];
0102         if(allcal_tows.find(hcal_tow) != allcal_tows.end()) allcal_tows[hcal_tow]+=te;
0103         else allcal_tows[hcal_tow]=te;
0104         float tea=te*1000.;
0105         emcal_energy->Fill(tea);
0106         for(int i=0; i<(int)tea; i++)emcal_hit_plot->Fill(i); 
0107         emcal_e+=te;
0108     }
0109     for(int n=0; n<(int) ihcal_tower_energy->size(); n++){
0110         ihcal_geom->set_calorimeter_id(RawTowerDefs::HCALIN);
0111         auto key=ihcal_tower_energy->encode_key(n);
0112         auto tower=ihcal_tower_energy->get_tower_at_channel(n);
0113         int phibin=ihcal_tower_energy->getTowerPhiBin(key);
0114         int etabin=ihcal_tower_energy->getTowerEtaBin(key);
0115         float phicenter=ihcal_geom->get_phicenter(phibin);
0116         float etacenter=ihcal_geom->get_etacenter(etabin);
0117         float te=tower->get_energy();
0118         if(this->local){
0119             te=te-tower->get_pedestal();
0120         }
0121         for(int i=0; i<(int)ihcal_energy_node.size(); i++)
0122             if(tower->get_energy() > energy_thresholds[i][2]) 
0123                 ihcal_energy_node[i]->Fill(etacenter, phicenter, tower->get_energy());
0124         if(tower->get_energy() > energy_thresholds[0][2])
0125              encCalc->addTower(n, ihcal_tower_energy, ihcal_geom, &ihcal_towers, RawTowerDefs::HCALIN);
0126         std::pair<float, float> hcal_tow {etacenter, phicenter};
0127         if(allcal_tows.find(hcal_tow) != allcal_tows.end()) allcal_tows[hcal_tow]+=te;
0128         else allcal_tows[hcal_tow]=te;
0129         float tea=te*1000.;
0130         ihcal_energy->Fill(tea);
0131         for(int i=0; i<(int)tea; i++)ihcal_hit_plot->Fill(i); 
0132         ihcal_e+=te;
0133     }
0134     for(int n=0; n<(int) ohcal_tower_energy->size(); n++){
0135         ohcal_geom->set_calorimeter_id(RawTowerDefs::HCALOUT);
0136         auto key=ohcal_tower_energy->encode_key(n);
0137         auto tower=ohcal_tower_energy->get_tower_at_channel(n);
0138         int phibin=ohcal_tower_energy->getTowerPhiBin(key);
0139         int etabin=ohcal_tower_energy->getTowerEtaBin(key);
0140         float phicenter=ohcal_geom->get_phicenter(phibin);
0141         float etacenter=ohcal_geom->get_etacenter(etabin);
0142         float te=tower->get_energy();
0143         if(this->local){
0144             te=te-tower->get_pedestal();
0145         }
0146         for(int i=0; i<(int)ohcal_energy_node.size(); i++)
0147             if(tower->get_energy() > energy_thresholds[i][3]) 
0148                 ohcal_energy_node[i]->Fill(etacenter, phicenter, tower->get_energy());
0149         if(tower->get_energy() > energy_thresholds[0][3])
0150              encCalc->addTower(n, ohcal_tower_energy, ohcal_geom, &ohcal_towers, RawTowerDefs::HCALOUT);    
0151         std::pair<float, float> hcal_tow {etacenter, phicenter};
0152         if(allcal_tows.find(hcal_tow) != allcal_tows.end()) allcal_tows[hcal_tow]+=te;
0153         else allcal_tows[hcal_tow]=te;
0154         float tea=te*1000.;
0155         ohcal_energy->Fill(tea);
0156         for(int i=0; i<(int)tea; i++)ohcal_hit_plot->Fill(i); 
0157         ohcal_e+=te;
0158     }
0159     float allcal_e=0.;
0160     for(auto a:allcal_tows){
0161         float em=a.second * 0.001;
0162         allcal_energy->Fill(em);
0163         allcal_e+=a.second;
0164         for(int i=0; i<(int)em; i++) allcal_hit_plot->Fill(i);
0165     }
0166     //try to build a jet like object
0167     std::map<std::array<float, 3>, float> total_array;
0168     for(auto e: emcal_towers)
0169     {
0170         for(int j=0; j<(int)emcal_energy_i.size(); j++)
0171             if(e.second > energy_thresholds[j][1]*16.) 
0172                 emcal_energy_i[j]->Fill(e.first[0], e.first[1], e.second);
0173         total_array[e.first]=e.second; //retower allready handeled
0174     }
0175     for(auto i: ihcal_towers)
0176     {
0177         for(int j=0; j<(int)ihcal_energy_i.size(); j++)
0178             if(i.second > energy_thresholds[j][2]) 
0179                 ihcal_energy_i[j]->Fill(i.first[0], i.first[1], i.second);
0180         if(total_array.find(i.first) != total_array.end()) total_array[i.first]+=i.second;
0181         else total_array[i.first]=i.second;
0182     }
0183     for(auto o: ohcal_towers)
0184     {
0185         for(int j=0; j<(int)ohcal_energy_i.size(); j++)
0186             if(o.second > energy_thresholds[j][3]) 
0187                 ohcal_energy_i[j]->Fill(o.first[0], o.first[1], o.second);
0188         if(total_array.find(o.first) != total_array.end()) total_array[o.first]+=o.second;
0189         else total_array[o.first]=o.second;
0190     }
0191     if(emcal_e > 0 ) emcal_energy_total->Fill(emcal_e);
0192     if(ohcal_e > 0 ) ohcal_energy_total->Fill(ohcal_e);
0193     if(ihcal_e > 0 ) ihcal_energy_total->Fill(ihcal_e);
0194     if(allcal_e > 0 ) allcal_energy_total->Fill(allcal_e);
0195     if(run_e2c){    
0196     JetContainerv1* jets=getJets(0.4, vertex, total_array);
0197     float lead_phi=0., max_e=0.;
0198     if(jets->size() > 0)
0199     {
0200         std::vector<float> fake_ratio;
0201         for(int i=0; i<(int)jets->size(); i++)fake_ratio.push_back(0.5);
0202         encCalc->eventCut->passesTheCut(jets, fake_ratio, fake_ratio, 0.5, vertex, fake_ratio);
0203         lead_phi = encCalc->eventCut->getLeadPhi();
0204     }
0205     else{
0206         for(auto t: total_array)
0207         {
0208             if(t.second > max_e)
0209             {
0210                 max_e = t.second;
0211                 lead_phi = t.first.at(1);
0212             }
0213             else continue;
0214         }
0215     }   
0216     encCalc->CaloRegion(emcal_towers, ihcal_towers, ohcal_towers, 0.05, "E",  vertex, lead_phi);
0217     }
0218     return Fun4AllReturnCodes::EVENT_OK;
0219 }
0220 JetContainerv1* LEDPedestalScan::getJets(float radius, std::array<float, 3> vertex, std::map<std::array<float, 3>, float> towers)
0221 {
0222     JetContainerv1* fastjetCont=new JetContainerv1();
0223     std::vector<fastjet::PseudoJet> jet_objs;
0224     fastjet::JetDefinition fjd(fastjet::antikt_algorithm, radius);
0225     for(auto t:towers){
0226         float E=t.second;
0227         if(E < 0.03) continue; //min 10 MeV per tower per detector
0228         float pz=E*tanh(t.first.at(0));
0229         float pt=E/(float)cosh(t.first.at(0));
0230         float px=pt*cos(t.first.at(1));
0231         float py=pt*sin(t.first.at(1));
0232         jet_objs.push_back(fastjet::PseudoJet(px, py, pz, E));
0233     }
0234     fastjet::ClusterSequence cs(jet_objs, fjd);
0235     auto js=cs.inclusive_jets();
0236     bool is_max=true;
0237     for(auto j:js){
0238         auto jet=fastjetCont->add_jet();
0239         jet->set_px(j.px());
0240         jet->set_py(j.py());
0241         jet->set_pz(j.pz());
0242         jet->set_e(j.e());
0243         for(auto cmp:j.constituents())
0244         {
0245             jet->insert_comp(Jet::SRC::CEMC_TOWER, cmp.user_index());
0246             jet_loc->Fill(cmp.eta(), cmp.phi());
0247             if(is_max) lead_jet_loc->Fill(cmp.eta(), cmp.phi());
0248         }
0249         jet_loc_center->Fill(j.eta(), j.phi());
0250         if(is_max){
0251             lead_jet_loc_center->Fill(j.eta(), j.phi());
0252             is_max=false;
0253         }
0254     }
0255     return fastjetCont;
0256 }
0257 void LEDPedestalScan::Print(const std::string &what) const
0258 {
0259     if(run_e2c) encCalc->Print(what);
0260     TFile* f1=new TFile(Form("LEDPedestal_scan_run_%d_segement_%d.root", run, segement), "RECREATE");
0261     jet_loc->Write();
0262     jet_loc_center->Write();
0263     lead_jet_loc->Write();
0264     lead_jet_loc_center->Write();
0265     for(int i= 0; i<(int)emcal_energy_i.size(); i++){
0266         emcal_energy_i[i]->Write();
0267         ihcal_energy_i[i]->Write();
0268         ohcal_energy_i[i]->Write();
0269         emcal_energy_node[i]->Write();
0270         ihcal_energy_node[i]->Write();
0271         ohcal_energy_node[i]->Write();
0272     }
0273     ohcal_energy->Write();
0274     ohcal_hit_plot->Write(); 
0275     ihcal_energy->Write();
0276     ihcal_hit_plot->Write(); 
0277     emcal_energy->Write();
0278     emcal_hit_plot->Write(); 
0279     allcal_energy->Write();
0280     allcal_hit_plot->Write(); 
0281     ohcal_energy_total->Write();
0282     ihcal_energy_total->Write();
0283     emcal_energy_total->Write();    
0284     allcal_energy_total->Write();
0285     f1->Close();
0286     return;
0287 }
0288
[ 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!