Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 //////////////////////////////////////////////////////////////////////////////////////////////////
0002 //                                              //
0003 //          Large R_L Calorimeter ENC                       //
0004 //                                              //
0005 //          Author: Skaydi Grossberndt                      //
0006 //          Depends on: Calorimeter Tower ENC                   //
0007 //          First Commit date: 18 Oct 2024                      //
0008 //          Most recent Commit: 17 Apr 2025                     //
0009 //          version: v5 single threshold value as derived from the values       //
0010 
0011 //                                              //
0012 //////////////////////////////////////////////////////////////////////////////////////////////////
0013 
0014 #include "LargeRLENC.h"
0015 
0016 
0017 LargeRLENC::LargeRLENC(const int n_run/*=0*/, const int n_segment/*=0*/, const float jet_min_pT/*=0*/, const bool data/*=false*/, const bool pedestal, std::fstream* ofs/*=nullptr*/, const std::string vari/*="E"*/, const std::string& name/* = "LargeRLENC"*/) 
0018 {
0019 
0020     n_evts=0;
0021     this->pedestalData=pedestal;
0022     if(pedestal){
0023         ohcal_min=0.014;
0024         emcal_min=0.062;
0025         ihcal_min=0.004;
0026         all_min=0.1;
0027     }
0028     thresh_mins[0]=all_min;
0029     thresh_mins[1]=emcal_min;
0030     thresh_mins[2]=ihcal_min;
0031     thresh_mins[3]=ohcal_min;
0032     thresh_mins[4]=0.1;
0033     MethodHistograms* fc, *fe, *fi, *fo, *tc, *te, *ti, *to, *ac, *ae, *ai, *ao, *trc, *tre, *tri, *tro;
0034 //set bin widths to tower size
0035     float allcal_thresh=1000*all_min;
0036     float emcal_thresh=1000*emcal_min;
0037     float ohcal_thresh= 1000*ohcal_min;
0038     float ihcal_thresh=1000*ihcal_min; //has much smaller gaps
0039     this->Thresholds=thresh_mins;
0040     std::string ihcal_thresh_s="_"+std::to_string((int)ihcal_thresh)+"_MeV_threshold";
0041     std::string ohcal_thresh_s="_"+std::to_string((int)ohcal_thresh)+"_MeV_threshold";
0042     std::string emcal_thresh_s="_"+std::to_string((int)emcal_thresh)+"_MeV_threshold";
0043     std::string allcal_thresh_s="_"+std::to_string((int)allcal_thresh)+"_MeV_threshold";
0044     
0045     fc=new MethodHistograms("Full_CAL"+allcal_thresh_s, 4*PI, 0.1); 
0046     fe=new MethodHistograms("Full_EMCAL"+emcal_thresh_s, 4*PI, 0.1);
0047     fi=new MethodHistograms("Full_IHCAL"+ihcal_thresh_s, 4*PI, 0.1);
0048     fo=new MethodHistograms("Full_OHCAL"+ohcal_thresh_s, 4*PI, 0.1); //Full OHCAL has R ~ 3.83 (delta eta ~ 2.2, delta phi ~pi)
0049 
0050     tc=new MethodHistograms("Towards_Region_CAL"+allcal_thresh_s, 4*PI, 0.1);
0051     te=new MethodHistograms("Towards_Region_EMCAL"+emcal_thresh_s, 4*PI, 0.1);
0052     ti=new MethodHistograms("Towards_Region_IHCAL"+ihcal_thresh_s, 4*PI, 0.1);
0053     to=new MethodHistograms("Towards_Region_OHCAL"+ohcal_thresh_s, 4*PI, 0.1); //Towards_Region OHCAL has R ~ 3.83 (delta eta ~ 2.2, delta phi ~2pi/3)
0054 
0055     ac=new MethodHistograms("Away_Region_CAL"+allcal_thresh_s, 4*PI, 0.1);
0056     ae=new MethodHistograms("Away_Region_EMCAL"+emcal_thresh_s, 4*PI, 0.1);
0057     ai=new MethodHistograms("Away_Region_IHCAL"+ihcal_thresh_s, 4*PI, 0.1);
0058     ao=new MethodHistograms("Away_Region_OHCAL"+ohcal_thresh_s, 4*PI, 0.1); //Away_Region OHCAL has R ~ 3.83 (delta eta ~ 2.2, delta phi ~2pi/3)
0059 
0060     trc=new MethodHistograms("Transverse_Region_CAL"+allcal_thresh_s, 4*PI, 0.1);
0061     tre=new MethodHistograms("Transverse_Region_EMCAL"+emcal_thresh_s, 4*PI, 0.1);
0062     tri=new MethodHistograms("Transverse_Region_IHCAL"+ihcal_thresh_s, 4*PI, 0.1);
0063     tro=new MethodHistograms("Transverse_Region_OHCAL"+ohcal_thresh_s, 4*PI, 0.1); //Transverse_Region OHCAL has R ~ 3.83 (delta eta ~ 2.2, delta phi ~pi)
0064     if(!pedestalData && !data){
0065         MethodHistograms* fc1=new MethodHistograms("Truth_Full_CAL"+allcal_thresh_s, 4*PI, 0.1); 
0066     
0067         MethodHistograms* tc1=new MethodHistograms("Truth_Towards_Region_CAL"+allcal_thresh_s, 4*PI, 0.1);
0068 
0069         MethodHistograms* ac1=new MethodHistograms("Truth_Away_Region_CAL"+allcal_thresh_s, 4*PI , 0.1);
0070     
0071         MethodHistograms* trc1=new MethodHistograms("Truth_Transverse_Region_CAL"+allcal_thresh_s, 4*PI, 0.1);
0072         std::vector<std::vector<MethodHistograms*>> Region_truth {std::vector<MethodHistograms*>{fc1}, std::vector<MethodHistograms*>{tc1}, std::vector<MethodHistograms*>{ac1}, std::vector<MethodHistograms*>{trc1}};
0073         this->Truth_Region_vector=Region_truth; 
0074     }   
0075     std::vector<MethodHistograms*> hf {fc, fe, fi, fo}, ht {tc, te, ti, to}, ha{ac, ae, ai, ao}, htr{trc, tre, tri, tro};
0076     std::vector<MethodHistograms*> FullCal=hf;
0077     std::vector<MethodHistograms*> TowardRegion=ht;
0078     std::vector<MethodHistograms*> AwayRegion=ha;
0079     std::vector<MethodHistograms*> TransverseRegion=htr;
0080     std::vector<std::vector<MethodHistograms*>> Region_vector_1;
0081     Region_vector_1.push_back(FullCal);
0082     Region_vector_1.push_back(TowardRegion);
0083     Region_vector_1.push_back(AwayRegion);
0084     Region_vector_1.push_back(TransverseRegion);
0085     if(pedestal || !data )
0086     {
0087         for(auto a:Region_vector_1)
0088             for(auto b:a){
0089                 b->E->SetCanExtend(TH1::kAllAxes);
0090                 b->N->SetCanExtend(TH1::kAllAxes);
0091                 b->pt->SetCanExtend(TH1::kAllAxes);
0092                 b->R_pt->SetCanExtend(TH1::kAllAxes);
0093             }
0094     } //led just does kinda odd things with the energy
0095     this->Region_vector=Region_vector_1;
0096     
0097 
0098 
0099     //while the Towards and away region arent't so extensive, haveing Rmax ~ 3, but close enough, transverse region is going to have some discontinuities
0100     this->isRealData=data; //check if the data is real data, if not run the truth as a cross check 
0101     this->nRun=n_run;
0102     this->nSegment=n_segment;
0103     this->jetMinpT=jet_min_pT;
0104     this->algo="Tower";
0105     this->radius="r04";
0106     this->eventCut=new DijetEventCuts();
0107     this->which_variable=vari;
0108     this->output_file_name="Large_RL_ENC_def_"+algo+"_dijets_anti_kT_"+radius+"-"+std::to_string(nRun)+"-"+std::to_string(nSegment)+".root";
0109     DijetQA=new TTree("DijetQA", "Dijet Event QA");
0110     DijetQA->Branch("N_jets", &m_Njets);
0111     DijetQA->Branch("x_j_L",  &m_xjl);
0112     DijetQA->Branch("A_jj_L", &m_Ajjl);
0113 
0114     DijetQA->Branch("Dijet_sets", &m_dijets, "dijet[3]/F"); 
0115     EEC=new TTree("EEC", "Energy Correlator");
0116     EEC->Branch("E_Total",   &m_etotal);
0117     EEC->Branch("E_CEMC",    &m_eemcal);
0118     EEC->Branch("E_IHCAL",   &m_eihcal);
0119     EEC->Branch("E_OHCAL",   &m_eohcal);
0120 
0121     EEC->Branch("Region",    &m_region);
0122     EEC->Branch("Calo",      &m_calo);
0123     EEC->Branch("R_L",   &m_rl);
0124     EEC->Branch("R_min",     &m_rm);
0125     EEC->Branch("R_i",   &m_ri);
0126 
0127     EEC->Branch("3_pt",      &m_e3c /*, "region/C:calo/C:r_l/F:e3c/F"*/);
0128 
0129     emcal_occup=new TH1F("emcal_occup", "Occupancy in the emcal in individual runs; Percent of Towers; N_{evts}", 100, -0.05, 99.5);
0130     ihcal_occup=new TH1F("ihcal_occup", "Occupancy in the ihcal in individual runs; Percent of Towers; N_{evts}", 100, -0.05, 99.5);
0131     ohcal_occup=new TH1F("ohcal_occup", "Occupancy in the ohcal in individual runs; Percent of Towers; N_{evts}", 100, -0.05, 99.5);
0132     ohcal_rat_h=new TH1F("ohcal_rat", "Ratio of energy in ohcal compared to total calorimeter; Ratio of OHCAL Energy; N_{evts}", 200, -0.5, 1.5);
0133     ohcal_rat_occup=new TH2F("ohcal_rat_occup", "Ratio of energy in ohcal and occupancy in the ohcal in individual runs; Energy deposited in OHCAL; Percent of Towers; N_{evts}", 150, 0.045, 0.995, 100, -0.05, 99.5);
0134     ohcal_bad_hits=new TH2F("ohcal_bad_hits", "#eta-#varphi energy deposition of \" Bad Hit\" events; #eta; #varphi; E [GeV]", 24, -1.1, 1.1, 64, -0.0001, 2*PI);
0135     bad_occ_em_oh=new TH2F("bad_occ_em_oh", "EMCAL to OHCAL tower deposition of \" Bad Hit\" events; Percent EMCAL towers; Percent OHCAL Towers; N_{evts}", 100, -0.5, 99.5, 100, -0.5, 99.5);
0136     bad_occ_em_h=new TH2F("em_allh_bad_hits", "Emcal_occ to Average hcal energy deposition of \" Bad Hit\" events;Percent EMCAL Towers; Average Percent HCAL towers; N_{evts}", 100, -0.5, 99.5, 100, -0.5, 99.5);
0137 
0138     IE=new TH1F("IE", "Radial Energy distribution of jets; #sum_{Constit} E_{i} R^{2} /E_{j}; N_{jets}", 1000, -0.05, 1.45); 
0139     badIE=new TH1F("badIE", "Radial Energy distribution of jets in events that fail cuts; #sum_{Constit} E_{i} R^{2} /E_{j}; N_{jets}", 1000, -0.05, 1.45); 
0140     goodIE=new TH1F("goodIE", "Radial Energy distribution of jets that pass dijet cuts; #sum_{Constit} E_{i} R^{2} /E_{j}; N_{jets}", 1000, -0.05, 1.45); 
0141     
0142     E_IE=new TH2F("E_IE", "Radial Energy Distribution of jets as a function of jet energy; E [GeV]; #sum E_{i} R^{2}/E_{j}; N_{jets}", 100, -0.5, 99.5, 100, -0.05, 1.45);
0143     badE_IE=new TH2F("badE_IE", "Radial Energy Distribution of jets as a function of jet energy failing dijet cuts; E [GeV]; #sum E_{i} R^{2}/E_{j}; N_{jets}", 100, -0.5, 99.5, 100, -0.05, 1.45);
0144     goodE_IE=new TH2F("goodE_IE", "Radial Energy Distribution of jets as a function of jet energy passing dijet cuts; E [GeV]; #sum E_{i} R^{2}/E_{j}; N_{jets}", 100, -0.5, 99.5, 100, -0.05, 1.45);
0145     MinpTComp=0.01; //10 MeV cut on tower/components
0146     for(int ci=0; ci < (int) Et_miss_hists.size(); ci++){
0147         std::string Calo_name;
0148         auto EC=&Et_miss_hists[ci];
0149         if(ci == 0){
0150             Calo_name="EMCAL";
0151         }
0152         else if(ci == 1){
0153             Calo_name = "IHCAL";
0154         }
0155         else if(ci == 2) {
0156             Calo_name = "OHCAL";
0157         }
0158         else{
0159             Calo_name= "ERR_CAL_NOT_FOUND";
0160         }
0161         for(int ri=0; ri<(int)EC->size(); ri++){
0162             float Rval=0.;
0163             if(ri==0) Rval=0.1;
0164             if(ri==1) Rval=0.4;
0165             if(ri==2) Rval=1.0;
0166             int r_no_dec=ri+1;
0167             EC->at(ri)=new TH1F(Form("h_%s_ETmiss_r_%d", Calo_name.c_str(), r_no_dec), Form("#tilde{E}_{T} of Jet Event Observable for %s at R=%f; #tilde{E}_{T} [GeV]; N_{event}", Calo_name.c_str(), Rval), 100, -0.5, 99.5);
0168         }
0169     }
0170     std::cout<<"Setup complete"<<std::endl;
0171 }
0172 bool LargeRLENC::triggerCut(bool use_em, PHCompositeNode* topNode)
0173 {
0174     bool good_trigger=false;
0175     TriggerAnalyzer* ta=new TriggerAnalyzer();
0176     ta->UseEmulator(use_em);
0177     ta->decodeTriggers(topNode);
0178     good_trigger = ta->didTriggerFire("Jet 10 GeV");
0179     return good_trigger;
0180 }
0181         
0182 JetContainerv1* LargeRLENC::getJets(std::string input, std::string radius, std::array<float, 3> vertex, float ohcal_rat, PHCompositeNode* topNode)
0183 {
0184     //This is just running the Fastjet reco 
0185     JetContainerv1* fastjetCont=new JetContainerv1();
0186     std::vector<fastjet::PseudoJet> jet_objs;
0187     float radius_float=0.;
0188     std::string rs=""; //striped down string to convert to float
0189     for(auto c:radius)
0190         if(isdigit(c))
0191             rs+=c;
0192     radius_float=stof(rs);
0193     radius_float=radius_float*0.1; //put the value to the correct range
0194     fastjet::JetDefinition fjd (fastjet::antikt_algorithm,  radius_float);
0195 
0196     if(input=="towers" || input.find("ower") != std::string::npos){
0197         TowerJetInput* ti_ohcal=new TowerJetInput(Jet::HCALOUT_TOWERINFO);
0198         TowerJetInput* ti_ihcal=new TowerJetInput(Jet::HCALIN_TOWERINFO);
0199         TowerJetInput* ti_emcal=new TowerJetInput(Jet::CEMC_TOWERINFO);
0200         auto psjets=ti_ohcal->get_input(topNode);
0201         auto psjets_ih=ti_ihcal->get_input(topNode);
0202         auto psjets_em=ti_emcal->get_input(topNode);
0203         for(auto p:psjets_ih) psjets.push_back(p);
0204         for(auto p:psjets_em) psjets.push_back(p);
0205         for(auto p:psjets) jet_objs.push_back(fastjet::PseudoJet(p->get_px(), p->get_py(), p->get_pz(), p->get_e()));
0206 
0207     }
0208     else if(input.find("luster") != std::string::npos){
0209         ClusterJetInput* ci=new ClusterJetInput(Jet::HCALOUT_CLUSTER);
0210         auto psjets=ci->get_input(topNode);
0211         std::cout<<"The tower jet input has n many objects " <<psjets.size() <<std::endl;
0212         for(auto p:psjets) jet_objs.push_back(fastjet::PseudoJet(p->get_px(), p->get_py(), p->get_pz(), p->get_e()));
0213     }
0214     else{
0215         std::cout<<"Could not recognize input method" <<std::endl;
0216     }
0217     fastjet::ClusterSequence cs(jet_objs, fjd); 
0218     auto js=cs.inclusive_jets();
0219     std::cout<<"Fastjet Clusterizer found : " <<js.size() <<std::endl;
0220     for(auto j:js)
0221     {
0222         auto jet=fastjetCont->add_jet();
0223         jet->set_px(j.px());
0224         jet->set_py(j.py());
0225         jet->set_pz(j.pz());
0226         jet->set_e( j.e() );
0227 
0228         for(auto cmp:j.constituents()){
0229             jet->insert_comp(Jet::SRC::HCALOUT_CLUSTER, cmp.user_index());
0230         }
0231     }
0232     return fastjetCont;
0233 }   
0234 std::array<float, 3> LargeRLENC::HadronicEnergyBalence(Jet* jet, float ohcal_ihcal_ratio, PHCompositeNode* topNode)
0235 {
0236     //For the case where we don't already have the source pieces
0237     float hadronic_energy=0., electromagnetic_energy=0.;
0238     float ohcal_conversion=ohcal_ihcal_ratio/(float)(ohcal_ihcal_ratio+1.);
0239     float jet_phi=jet->get_phi(), jet_eta=jet->get_eta();
0240     float i_e=0.;
0241     try{
0242         findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0243     }
0244     catch(std::exception& e){ return {0.,0., 0.};}
0245     auto truth_particles_p=findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0246     try{
0247         truth_particles_p->GetMap();
0248     }
0249     catch(std::exception& e){std::cout<<"Could not find particle map" <<std::endl;
0250         return {0.,0.,0.};
0251     }
0252     std::map<int, PHG4Particle*> truth_particles;
0253     for(const auto& a:truth_particles_p->GetMap()) truth_particles[a.first]=a.second;
0254     for(auto& iter:jet->get_comp_vec()){
0255         Jet::SRC source=iter.first;
0256         if(source != Jet::SRC::PARTICLE && source != Jet::SRC::CHARGED_PARTICLE && source != Jet::SRC::HEPMC_IMPORT){
0257             //don't have a source particle so skip it
0258             continue;
0259         }
0260         else{
0261             unsigned int id=iter.second;
0262             if(truth_particles.find(id) == truth_particles.end()){
0263                 continue;
0264             }
0265             else{
0266                 PHG4Particle* particle = truth_particles.at(id);
0267                 int pid=particle->get_pid();
0268                 if(abs(pid) == 11 || pid== 22){ 
0269                     //electrons, positrons and photons get put in the emcal
0270                     electromagnetic_energy+=particle->get_e();
0271                     float particle_phi=std::atan2(particle->get_py(), particle->get_px());
0272                     float particle_eta=std::atanh(particle->get_pz()/particle->get_e());
0273                     i_e+=particle->get_e()*std::pow(getR(particle_eta, particle_phi, jet_eta, jet_phi),2);
0274                 }
0275                 else if(abs(pid) > 11 && abs(pid) <= 18){
0276                     //don't count neutrinos, muons, tau
0277                     hadronic_energy+=particle->get_e();
0278                     float particle_phi=std::atan2(particle->get_py(), particle->get_px());
0279                     float particle_eta=std::atanh(particle->get_pz()/particle->get_e());
0280                     i_e+=particle->get_e()*std::pow(getR(particle_eta, particle_phi, jet_eta, jet_phi),2);
0281                 }
0282             }
0283         }
0284     }
0285     //assume that the hcal energy split is consistent with the whole detector energy split 
0286     float ohcal_energy=hadronic_energy*ohcal_conversion;
0287     float ohcal_ratio=ohcal_energy/(hadronic_energy+electromagnetic_energy);
0288     float emcal_ratio=electromagnetic_energy/(hadronic_energy+electromagnetic_energy);
0289     i_e=i_e/(hadronic_energy+electromagnetic_energy);
0290     return {ohcal_ratio, emcal_ratio, i_e};
0291 }
0292 std::vector<std::array<float,3>> LargeRLENC::getJetEnergyRatios(JetContainerv1* jets, float ohcal_ihcal_ratio, PHCompositeNode* topNode)
0293 {
0294     //get the energy ballence in each calorimeter for each jet
0295     std::vector<std::array<float, 3>> ohcal_ratio;
0296     auto emcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode,  emcal_energy_towers );
0297     auto ihcal_tower_energy= findNode::getClass<TowerInfoContainer>(topNode, ihcal_energy_towers );
0298     auto ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode,  ohcal_energy_towers );
0299     auto ohcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0300     auto emcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC"   );
0301     auto ihcal_geom= findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0302     ohcal_geom->set_calorimeter_id(RawTowerDefs::HCALOUT);
0303     ihcal_geom->set_calorimeter_id(RawTowerDefs::HCALIN);
0304     emcal_geom->set_calorimeter_id(RawTowerDefs::CEMC);
0305     for(auto j: *jets)
0306     {
0307         if(!j) continue;
0308         bool has_particle=false;
0309         float ohcal_energy=0., allcal_energy=0., emcal_energy=0.; 
0310         float i_e=0;
0311         float jet_phi=j->get_phi(), jet_eta=j->get_eta();
0312         auto cmp_vec=j->get_comp_vec(); 
0313         for(auto iter:cmp_vec){
0314             Jet::SRC source=iter.first; //get source of object
0315             if(source == Jet::SRC::PARTICLE || source == Jet::SRC::CHARGED_PARTICLE || source == Jet::SRC::HEPMC_IMPORT){
0316                 has_particle=true; 
0317                 //if any particle source is found, we can use that. otherwise have to not use this cut
0318             }   
0319             if(  source==Jet::SRC::HEPMC_IMPORT || source > Jet::SRC::HCALOUT_TOWERINFO_SIM || source==Jet::SRC::HCALOUT_CLUSTER || source==Jet::SRC::HCALIN_CLUSTER)continue;
0320             else{
0321                 if(source==Jet::SRC::HCALOUT_TOWER|| source == Jet::SRC::HCALOUT_CLUSTER){
0322                     unsigned int tower_id=iter.second;
0323                     float e=0.;
0324                     try{
0325                         ohcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0326                     }
0327                     catch(std::exception& e){std::cout<<"Bad tower id found for source " <<Jet::SRC::HCALOUT_TOWER <<" despite actual source being " <<source <<" and tower id is " <<tower_id<<std::endl;}
0328                     ohcal_energy+= e;//functionally the same as below, I just want to make things readable
0329                     allcal_energy+= e;
0330                     int phibin=ohcal_tower_energy->getTowerPhiBin(tower_id);
0331                     int etabin=ohcal_tower_energy->getTowerEtaBin(tower_id);
0332                     float phicenter=ohcal_geom->get_phicenter(phibin);
0333                     float etacenter=ohcal_geom->get_etacenter(etabin);  
0334                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0335                 }
0336                 else if( source == Jet::SRC::HCALOUT_TOWER_SUB1  || source == Jet::SRC::HCALOUT_TOWERINFO_SUB1){
0337                     unsigned int tower_id=iter.second;
0338                     ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0339                     float e=0.;
0340                     auto key_id=ohcal_tower_energy->encode_key(tower_id);
0341                     e=ohcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0342                     ohcal_energy+= e;//functionally the same as below, I just want to make things readable
0343                     allcal_energy+= e;
0344                     int phibin=ohcal_tower_energy->getTowerPhiBin(key_id);
0345                     int etabin=ohcal_tower_energy->getTowerEtaBin(key_id);
0346                     float phicenter=ohcal_geom->get_phicenter(phibin);
0347                     float etacenter=ohcal_geom->get_etacenter(etabin);  
0348                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0349                 }
0350                 else if(source== Jet::SRC::HCALOUT_TOWERINFO || source == Jet::SRC::HCALOUT_TOWER_SUB1CS ){
0351                     unsigned int tower_id=iter.second;
0352                     float e=ohcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0353                     ohcal_energy+= e;//functionally the same as below, I just want to make things readable
0354                     allcal_energy+= e;
0355                     int phibin=ohcal_tower_energy->getTowerPhiBin(tower_id);
0356                     int etabin=ohcal_tower_energy->getTowerEtaBin(tower_id);
0357                     float phicenter=ohcal_geom->get_phicenter(phibin);
0358                     float etacenter=ohcal_geom->get_etacenter(etabin);  
0359                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0360                 }
0361                 else if(source== Jet::SRC::HCALOUT_TOWERINFO_SIM  || source==Jet::SRC::HCALOUT_TOWERINFO_EMBED ){
0362                     unsigned int tower_id=iter.second;
0363                     float e=ohcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0364                     ohcal_energy+= e;//functionally the same as below, I just want to make things readable
0365                     allcal_energy+= e;
0366                     int phibin=ohcal_tower_energy->getTowerPhiBin(tower_id);
0367                     int etabin=ohcal_tower_energy->getTowerEtaBin(tower_id);
0368                     float phicenter=ohcal_geom->get_phicenter(phibin);
0369                     float etacenter=ohcal_geom->get_etacenter(etabin);  
0370                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0371                 }
0372                 else if(source== Jet::SRC::HCALIN_TOWER || source == Jet::SRC::HCALIN_CLUSTER  ){
0373                     unsigned int tower_id=iter.second;
0374                     float e=ihcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0375                     allcal_energy+= e;
0376                     int phibin=ihcal_tower_energy->getTowerPhiBin(tower_id);
0377                     int etabin=ihcal_tower_energy->getTowerEtaBin(tower_id);
0378                     float phicenter=ihcal_geom->get_phicenter(phibin);
0379                     float etacenter=ihcal_geom->get_etacenter(etabin);  
0380                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0381                 }
0382                 else if( source == Jet::SRC::HCALIN_TOWER_SUB1 || source == Jet::SRC::HCALIN_TOWERINFO_SUB1){
0383                     unsigned int tower_id=iter.second;
0384                     ihcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0385                     float e=0.;
0386                     e=ihcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0387                     allcal_energy+= e;
0388                     auto key=ihcal_tower_energy->encode_key(tower_id);
0389                     int phibin=ihcal_tower_energy->getTowerPhiBin(key);
0390                     int etabin=ihcal_tower_energy->getTowerEtaBin(key);
0391                     float phicenter=ihcal_geom->get_phicenter(phibin);
0392                     float etacenter=ihcal_geom->get_etacenter(etabin);  
0393                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0394                 }
0395                 else if(source== Jet::SRC::HCALIN_TOWERINFO || source == Jet::SRC::HCALIN_TOWER_SUB1CS ){
0396                     unsigned int tower_id=iter.second;
0397                     float e=ihcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0398                     int phibin=ihcal_tower_energy->getTowerPhiBin(tower_id);
0399                     int etabin=ihcal_tower_energy->getTowerEtaBin(tower_id);
0400                     float phicenter=ihcal_geom->get_phicenter(phibin);
0401                     float etacenter=ihcal_geom->get_etacenter(etabin);  
0402                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0403                     allcal_energy+= e;
0404                 }
0405                 else if(source== Jet::SRC::HCALIN_TOWERINFO_SIM ||  source==Jet::SRC::HCALIN_TOWERINFO_EMBED ){
0406                     unsigned int tower_id=iter.second;
0407                     float e=ihcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0408                     int phibin=ihcal_tower_energy->getTowerPhiBin(tower_id);
0409                     int etabin=ihcal_tower_energy->getTowerEtaBin(tower_id);
0410                     float phicenter=ihcal_geom->get_phicenter(phibin);
0411                     float etacenter=ihcal_geom->get_etacenter(etabin);  
0412                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0413                     allcal_energy+= e;
0414                 }
0415                 
0416                 else if(source== Jet::SRC::CEMC_TOWER || source == Jet::SRC::CEMC_CLUSTER ){
0417                     unsigned int tower_id=iter.second;
0418                     float e=emcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0419                     int phibin=emcal_tower_energy->getTowerPhiBin(tower_id);
0420                     int etabin=emcal_tower_energy->getTowerEtaBin(tower_id);
0421                     float phicenter=emcal_geom->get_phicenter(phibin);
0422                     float etacenter=emcal_geom->get_etacenter(etabin);  
0423                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0424                     allcal_energy+= e;
0425                     emcal_energy+=e;
0426                 }
0427                 else if( source == Jet::SRC::CEMC_TOWER_SUB1 || source== Jet::SRC::CEMC_TOWERINFO_SUB1 ){
0428                     unsigned int tower_id=iter.second;
0429                     emcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0430                     float e=0.;
0431                     e=emcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0432                     allcal_energy+= e;
0433                     emcal_energy+=e;
0434                     auto key=emcal_tower_energy->encode_key(tower_id);
0435                     int phibin=emcal_tower_energy->getTowerPhiBin(key);
0436                     int etabin=emcal_tower_energy->getTowerEtaBin(key);
0437                     float phicenter=emcal_geom->get_phicenter(phibin);
0438                     float etacenter=emcal_geom->get_etacenter(etabin);  
0439                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0440                 }
0441                 else if(source== Jet::SRC::CEMC_TOWERINFO || source == Jet::SRC::CEMC_TOWER_SUB1CS ){
0442                     unsigned int tower_id=iter.second;
0443                     float e=emcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0444                     int phibin=emcal_tower_energy->getTowerPhiBin(tower_id);
0445                     int etabin=emcal_tower_energy->getTowerEtaBin(tower_id);
0446                     float phicenter=emcal_geom->get_phicenter(phibin);
0447                     float etacenter=emcal_geom->get_etacenter(etabin);  
0448                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0449                     allcal_energy+= e;
0450                     emcal_energy+=e;
0451                 }
0452                 else if(source== Jet::SRC::CEMC_TOWERINFO_SIM ||  source==Jet::SRC::CEMC_TOWERINFO_EMBED ){
0453                     unsigned int tower_id=iter.second;
0454                     float e=emcal_tower_energy->get_tower_at_channel(tower_id)->get_energy();
0455                     int phibin=emcal_tower_energy->getTowerPhiBin(tower_id);
0456                     int etabin=emcal_tower_energy->getTowerEtaBin(tower_id);
0457                     float phicenter=emcal_geom->get_phicenter(phibin);
0458                     float etacenter=emcal_geom->get_etacenter(etabin);  
0459                     i_e+=e*std::pow(getR( etacenter, phicenter, jet_eta, jet_phi),2);
0460                     allcal_energy+= e;
0461                     emcal_energy+=e;
0462                 }
0463             }
0464         }
0465         if(ohcal_energy == 0. || allcal_energy == 0. || emcal_energy == 0.)
0466         {
0467     //      std::cout<<"There is an issue getting an energy value, please check status?"<<ohcal_energy <<"   " <<allcal_energy <<std::endl;
0468             if(has_particle) ohcal_ratio.push_back(HadronicEnergyBalence(j, ohcal_ihcal_ratio, topNode)); //if there is no tower source present, use the particles  
0469             else ohcal_ratio.push_back({0.,0.,0.});
0470         }
0471         else{
0472             i_e=(float)i_e/(float)allcal_energy;
0473     //      std::cout<<"The jet moment of inertia is " <<i_e <<std::endl;
0474             ohcal_energy=ohcal_energy/allcal_energy;
0475             emcal_energy=emcal_energy/allcal_energy;
0476             ohcal_ratio.push_back({ohcal_energy, emcal_energy, i_e});
0477         }
0478     }
0479     return ohcal_ratio;
0480 }
0481 
0482 void LargeRLENC::addTower(int n, TowerInfoContainer* energies, RawTowerGeomContainer_Cylinderv1* geom, std::map<std::array<float, 3>, float>* towers, RawTowerDefs::CalorimeterId td)
0483 {
0484     if(!geom) return;
0485     geom->set_calorimeter_id(td);
0486 //  std::cout<<"The map has size " <<towers->size() <<std::endl;
0487     auto key=energies->encode_key(n);
0488     auto tower=energies->get_tower_at_channel(n);
0489     int phibin=energies->getTowerPhiBin(key);
0490     int etabin=energies->getTowerEtaBin(key);
0491     float phicenter=geom->get_phicenter(phibin);
0492     float etacenter=geom->get_etacenter(etabin);
0493     float r=geom->get_radius();
0494 //  std::cout<<"energy is " <<tower->get_energy()<<std::endl;
0495     std::array<float, 3> center {etacenter, phicenter, r};
0496 
0497     if(td != RawTowerDefs::CEMC && (td != RawTowerDefs::HCALOUT && this->pedestalData)) towers->insert(std::make_pair(center, tower->get_energy()));
0498     else if (td == RawTowerDefs::HCALOUT && this->pedestalData) towers->insert(std::make_pair(center, tower->get_energy()/100.));
0499 
0500     else if(td==RawTowerDefs::CEMC){
0501         //retowering it by hand for right now to improve running speed 
0502     /*  for(int j=0; j<24; j++){
0503             float hcalbinval=((j+1)*1.1/12.)-1.1;
0504             if(center[0] < hcalbinval){
0505                 center[0]=hcalbinval;
0506             }
0507             else break;
0508         }
0509         for(auto j=0; j<64; j++)
0510         {
0511             float hcalbinval=(j+1)*2*PI/64.;
0512             if(center[1] < hcalbinval)  center[1]=hcalbinval;
0513             else break;
0514         }*/
0515         center[0]=emcal_lookup_table[n].first;
0516         center[1]=emcal_lookup_table[n].second;
0517         if(towers->find(center) != towers->end()) towers->at(center)+=tower->get_energy();
0518         else towers->insert(std::make_pair(center, tower->get_energy()));
0519     }
0520     return;
0521 }
0522 void LargeRLENC::MakeEMCALRetowerMap(RawTowerGeomContainer_Cylinderv1* em_geom, TowerInfoContainer* emcal, RawTowerGeomContainer_Cylinderv1* h_geom, TowerInfoContainer* hcal )
0523 {
0524     em_geom->set_calorimeter_id(RawTowerDefs::CEMC);
0525     h_geom->set_calorimeter_id(RawTowerDefs::HCALOUT);
0526 
0527     for(int n=0; n<(int) emcal->size(); n++){
0528         auto key=emcal->encode_key(n);
0529         int phibin=emcal->getTowerPhiBin(key);
0530         int etabin=emcal->getTowerEtaBin(key);
0531         float phicenter=em_geom->get_phicenter(phibin);
0532         float etacenter=em_geom->get_etacenter(etabin);
0533         for(int j=0; j<(int) hcal->size(); j++)
0534         {
0535             bool goodPhi=false, goodEta=false;
0536             auto key_h=hcal->encode_key(j);
0537             int phibin_h=hcal->getTowerPhiBin(key_h);
0538             int etabin_h=hcal->getTowerEtaBin(key_h);
0539             float phicenter_h=h_geom->get_phicenter(phibin_h);
0540             float etacenter_h=h_geom->get_etacenter(etabin_h);  
0541             std::pair<double, double> phi_bounds=h_geom->get_phibounds(phibin_h);
0542             std::pair<double, double> eta_bounds=h_geom->get_etabounds(etabin_h);
0543             if(phicenter >= phi_bounds.first && phicenter < phi_bounds.second) goodPhi=true; 
0544             if(etacenter >= eta_bounds.first && etacenter < eta_bounds.second) goodEta=true;
0545             if(goodPhi && goodEta){
0546                 this->emcal_lookup_table[n]=std::make_pair(etacenter_h, phicenter_h);
0547                 break;
0548             }
0549             else continue;
0550         }
0551     }
0552 
0553     return;
0554 }
0555 int LargeRLENC::process_event(PHCompositeNode* topNode)
0556 {
0557     //This is where I can be a bit clever with allowing for the parrelization for easier crosschecks
0558     n_evts++;
0559     std::cout<<"Running on event : " <<n_evts<<std::endl;
0560 
0561     m_emcal.clear();
0562     m_ihcal.clear();
0563     m_ohcal.clear();
0564     std::cout<<"Cleared the vectors" <<std::endl;
0565     auto emcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC"   );
0566 
0567     auto emcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, emcal_energy_towers      );
0568     auto ihcal_geom= findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0569     auto ihcal_tower_energy= findNode::getClass<TowerInfoContainer>(topNode, ihcal_energy_towers     );
0570     auto ohcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0571     auto ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, ohcal_energy_towers      );
0572     //look for the jet objects 
0573     JetContainerv1* jets=NULL;
0574     bool /*foundJetConts=false, */isDijet=false;
0575     if(this->emcal_lookup_table.size() == 0  || n_evts==1)
0576     {
0577         MakeEMCALRetowerMap(emcal_geom, emcal_tower_energy, ohcal_geom, ohcal_tower_energy);
0578         std::cout<<"Lookup table has size: " <<this->emcal_lookup_table.size() <<std::endl;
0579     }
0580     try{
0581         if(!isRealData) jets = findNode::getClass<JetContainerv1>(topNode, "AntiKt_Truth_r04"); //look for the r_04 truth jets
0582         else{
0583             std::string recoJetName="AntiKt_"+algo+radius+"_Sub1";
0584             jets = findNode::getClass<JetContainerv1>(topNode, recoJetName); //check for already reconstructed jets
0585             if(!jets || jets->size() == 0 ){
0586                 jets = findNode::getClass<JetContainerv1>(topNode, "AntiKt_Tower_r04_Sub1"); //explicit backup check
0587                 std::cout<<"Jet container for tower sub1 has size: " <<jets->size() <<std::endl;
0588                 if(!jets || jets->size() == 0){
0589                     jets = findNode::getClass<JetContainerv1>(topNode, "AntiKt_TowerInfo_r04");
0590                     if(jets) std::cout<<"Jet container for towerinfo has size: " <<jets->size() <<std::endl;
0591                 }
0592                 if(!jets || jets->size() == 0){
0593                     jets = findNode::getClass<JetContainerv1>(topNode, "AntiKt_unsubtracted_r04");
0594                     if(jets) std::cout<<"Jet container for unsub has size: " <<jets->size() <<std::endl;
0595                 }
0596                 //try every possible style
0597             }
0598     /*      else {
0599                 foundJetConts=true;
0600                 }*/
0601         }
0602     }
0603     catch(std::exception& e){ std::cout<<"Did not find a jet container object, will attempt to reconstruct the jets" <<std::endl;}
0604 //  if(jets==NULL || !foundJetConts) jets=getJets("Tower", radius, vertex, ohcal_rat, topNode); //if needed will go back to this
0605     if(!jets){
0606         std::cout<<"Didn't find any jets, skipping event" <<std::endl;
0607         return Fun4AllReturnCodes::EVENT_OK;
0608     }
0609     else if(jets->size() == 0){
0610         std::cout<<"Jet Container empty, skipping event" <<std::endl;
0611         return Fun4AllReturnCodes::EVENT_OK;
0612     }
0613     std::array<float, 3> vertex={0.,0.,0.}; //set the initial vertex to origin 
0614     try{
0615         GlobalVertexMap* vertexmap=findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0616 
0617         if(vertexmap){
0618             if(vertexmap->empty())
0619                 std::cout<<"Empty Vertex Map. \n Setting vertex to origin" <<std::endl; 
0620             else{
0621 
0622                 GlobalVertex* gl_vtx=nullptr;
0623                 for(auto vertex_iter:*vertexmap){
0624                     if(vertex_iter.first == GlobalVertex::VTXTYPE::MBD || vertex_iter.first == GlobalVertex::VTXTYPE::SVTX_MBD ) 
0625                     { 
0626                         gl_vtx=vertex_iter.second;
0627                     }
0628                 }
0629                 if(gl_vtx){
0630                     vertex[0]=gl_vtx->get_x();
0631                     m_vtx=vertex[0];
0632                     vertex[1]=gl_vtx->get_y();
0633                     m_vty=vertex[1];
0634                     vertex[2]=gl_vtx->get_z();
0635 
0636                     m_vtz=vertex[2];
0637                 }
0638             }
0639         }
0640     }
0641     catch(std::exception& e){std::cout<<"Could not find the vertex. \n Setting to origin" <<std::endl;}
0642 
0643     int emcal_oc=0; //allow for occupancy to be calculated seperate from the other bits
0644     float emcal_energy=0., ihcal_energy=0., ohcal_energy=0., total_energy=0., ohcal_rat=0.;
0645     std::map<std::array<float, 3>, float> ihcal_towers, emcal_towers, ohcal_towers, truth_pts; //these are just to collect the non-zero towers to decrease the processing time 
0646     for(int n=0; n<(int) emcal_tower_energy->size(); n++){
0647         if(! emcal_tower_energy->get_tower_at_channel(n)->get_isGood()) continue;
0648         float energy=emcal_tower_energy->get_tower_at_channel(n)->get_energy();
0649         if(energy > emcal_min )//put zero supression into effect
0650             addTower(n, emcal_tower_energy, emcal_geom, &emcal_towers, RawTowerDefs::CEMC   );
0651         if(energy > emcal_min) emcal_oc++;
0652         emcal_energy+=energy;
0653     }
0654     for(int n=0; n<(int) ihcal_tower_energy->size(); n++){
0655         if(n >= (int)ihcal_tower_energy->size()) continue;
0656         if(! ihcal_tower_energy->get_tower_at_channel(n)->get_isGood()) continue;
0657         float energy=ihcal_tower_energy->get_tower_at_channel(n)->get_energy();
0658         if(energy > ihcal_min) //zero suppression is at 10 MeV in IHCAL
0659             addTower(n, ihcal_tower_energy, ihcal_geom, &ihcal_towers, RawTowerDefs::HCALIN  );
0660         ihcal_energy+=energy;
0661     }
0662     for(int n=0; n<(int) ohcal_tower_energy->size(); n++){
0663         if(! ohcal_tower_energy->get_tower_at_channel(n)->get_isGood()) continue;
0664         float energy=ohcal_tower_energy->get_tower_at_channel(n)->get_energy();
0665         if(energy > ohcal_min ) //10 MeV cutoff 
0666             addTower(n, ohcal_tower_energy, ohcal_geom, &ohcal_towers, RawTowerDefs::HCALOUT );
0667         ohcal_energy+=energy;
0668     }
0669     
0670     if(!isRealData && !pedestalData){
0671         auto hepmc_gen_event= findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0672         if(hepmc_gen_event){
0673             for( PHHepMCGenEventMap::ConstIter evtIter=hepmc_gen_event->begin(); evtIter != hepmc_gen_event->end(); ++evtIter)
0674             {
0675                 PHHepMCGenEvent* hpev=evtIter->second;
0676                 if(hpev){
0677                     HepMC::GenEvent* ev=hpev->getEvent();   
0678                     if(ev)
0679                     {
0680                         for(HepMC::GenEvent::particle_const_iterator iter=ev->particles_begin(); iter !=ev->particles_end(); ++iter){
0681                             if((*iter))
0682                             {
0683                                 if(!(*iter)->end_vertex() && (*iter)->status() == 1){
0684                                     float px=(*iter)->momentum().px();
0685                                     float py=(*iter)->momentum().py();
0686                                     float pz=(*iter)->momentum().pz();
0687                                     float E=(*iter)->momentum().e();
0688                                     float phi=atan2(py, px);
0689                                     float eta=atanh(px/(std::sqrt(px*px+py*py+pz*pz)));
0690                                     float r=1.;
0691                                     std::array<float, 3> loc {eta, phi, r};
0692                                     truth_pts[loc]=E;
0693 
0694                                 }
0695                             }
0696                         }
0697                     }
0698                 }
0699             }
0700         }
0701     }   
0702 
0703     total_energy=emcal_energy+ihcal_energy+ohcal_energy;
0704     ohcal_rat=ohcal_energy/(float)total_energy; //take the ratio at the whole calo as that is the region of interest
0705     float emcal_occupancy=emcal_oc/((float)emcal_tower_energy->size())*100.;
0706     float ihcal_occupancy=ihcal_towers.size()/((float)ihcal_tower_energy->size())*100.;
0707     float ohcal_occupancy=ohcal_towers.size()/((float)ohcal_tower_energy->size())*100.;
0708     emcal_occup->Fill(emcal_occupancy); //filling in the occupancy plots for easy QA while going 
0709     ihcal_occup->Fill(ihcal_occupancy);
0710     ohcal_occup->Fill(ohcal_occupancy);
0711     ohcal_rat_h->Fill(ohcal_rat);
0712 
0713     std::vector<float> ohcal_jet_rat, emcal_jet_rat, jet_ie;
0714     std::vector<std::array<float, 3>> ohcal_jet_rat_and_ie=getJetEnergyRatios(jets, ohcal_energy/(float)ihcal_energy, topNode);
0715     for(auto h:ohcal_jet_rat_and_ie){
0716         ohcal_jet_rat.push_back(h[0]);
0717         emcal_jet_rat.push_back(h[1]);
0718         jet_ie.push_back(h[2]);
0719     }
0720     int i1=0;
0721     for(auto j: *jets){
0722         if(jet_ie.at(i1) > 0 )IE->Fill(jet_ie.at(i1));
0723         if(jet_ie.at(i1) > 0 )E_IE->Fill(j->get_e(), jet_ie.at(i1));
0724         i1++;
0725     }
0726     bool triggered_event = true;
0727     if(isRealData){
0728         try{
0729             auto gl1packet = findNode::getClass<Gl1Packet>(topNode, "Gl1Packet");
0730             if(gl1packet){
0731                 triggered_event = triggerCut(false, topNode);
0732             }
0733         }
0734         catch(std::exception& e){std::cout<<"Could not find the trigger object" <<std::endl;}
0735     }
0736         
0737     isDijet=eventCut->passesTheCut(jets, ohcal_jet_rat, emcal_jet_rat, ohcal_rat, vertex, jet_ie);  
0738     if(!isDijet || !triggered_event){ //stores some data about the bad cuts to look for any arrising structure
0739         ohcal_rat_occup->Fill(ohcal_rat, ohcal_occupancy);
0740         if(ohcal_rat > 0.9){
0741                 for(auto p:ohcal_towers)ohcal_bad_hits->Fill(p.first[0], p.first[1], p.second);
0742         }
0743         bad_occ_em_oh->Fill(emcal_occupancy, ohcal_occupancy);
0744         bad_occ_em_h->Fill(emcal_occupancy, (ohcal_occupancy+ihcal_occupancy)/2.0);
0745         i1=0;
0746         for(auto j:*jets){
0747             if(jet_ie.at(i1) > 0 ) badIE->Fill(jet_ie.at(i1));
0748             if(jet_ie.at(i1) > 0 )badE_IE->Fill(j->get_e(), jet_ie.at(i1));
0749             i1++;
0750         }
0751             return Fun4AllReturnCodes::EVENT_OK;
0752     }
0753     else{
0754         i1=0;
0755         for(auto j:*jets){
0756             if(jet_ie.at(i1) > 0 )goodIE->Fill(jet_ie.at(i1));
0757             if(jet_ie.at(i1) > 0 )goodE_IE->Fill(j->get_e(), jet_ie.at(i1));
0758             i1++;
0759         }
0760         float lead_phi=eventCut->getLeadPhi();
0761         eventCut->getDijets(jets, &m_dijets);
0762 
0763         m_vertex=vertex;
0764         m_etotal=total_energy;
0765         m_eemcal=emcal_energy;
0766         m_eihcal=ihcal_energy;
0767         m_eohcal=ohcal_energy;
0768         m_philead=eventCut->getLeadPhi();
0769         m_etalead=eventCut->getLeadEta();
0770 
0771         int i=0;
0772         for(auto j:*jets){
0773             if(!j) continue;
0774             Region_vector[0][0]->pt->Fill(j->get_pt());
0775             Region_vector[0][1]->pt->Fill(j->get_pt()*(emcal_jet_rat.at(i)));
0776             Region_vector[0][2]->pt->Fill(j->get_pt()*(1-emcal_jet_rat.at(i)-ohcal_jet_rat.at(i)));
0777             Region_vector[0][3]->pt->Fill(j->get_pt()*(ohcal_jet_rat.at(i)));
0778             i++;
0779         }
0780         if(isDijet) n_with_jets++;
0781         //this is running the actual analysis
0782         LargeRLENC::CaloRegion(emcal_towers, ihcal_towers, ohcal_towers, jetMinpT, which_variable, vertex, lead_phi);
0783         if(!isRealData && !pedestalData) LargeRLENC::TruthRegion(truth_pts, jetMinpT, which_variable, vertex, lead_phi); 
0784         DijetQA->Fill();
0785     }
0786     return Fun4AllReturnCodes::EVENT_OK;
0787 }
0788 void LargeRLENC::TruthRegion(std::map<std::array<float, 3>, float> truth, float jetMinpT, std::string variable_of_interest, std::array<float, 3> vertex, float lead_jet_center_phi)
0789 {
0790     bool transverse=false;
0791     bool energy=true;
0792     float phi_min=0., phi_max=2*PI;
0793     int v_o_i_code=0;
0794     std::map<int, std::pair<float, float>> region_ints;
0795     if(variable_of_interest.find("T") != std::string::npos) v_o_i_code+=2;
0796     if(variable_of_interest.find("p") != std::string::npos) v_o_i_code+=1; 
0797     switch(v_o_i_code)
0798     {
0799         case 0:
0800             transverse=false;
0801             energy=true;
0802             break;
0803         case 1:
0804             transverse=false;
0805             energy=false;
0806             break;
0807         case 2:
0808             transverse=true;
0809             energy=true;
0810             break;
0811         case 3:
0812             transverse=true;
0813             energy=false;
0814             break;
0815     }
0816     for(int region=0; region< 5; region++)
0817     {
0818         switch(region)
0819         {
0820             case LargeRLENC::Regions::Full: //go over the full calorimeter
0821                 phi_min=0;
0822                 phi_max=2*PI;
0823                 break;
0824             case LargeRLENC::Regions::Towards: //toward region of the detector
0825                 phi_min=lead_jet_center_phi - PI/3.;
0826                 phi_max=lead_jet_center_phi + PI/3.;
0827                 break;
0828             case LargeRLENC::Regions::Away: //away region of the detector
0829                 phi_min=lead_jet_center_phi + 2*PI/3.;
0830                 phi_max=lead_jet_center_phi + 4*PI/3.;
0831                 break;
0832             case LargeRLENC::Regions::TransverseMax: //transverse region1
0833                 phi_min=lead_jet_center_phi + PI/3.;
0834                 phi_max=lead_jet_center_phi + 2*PI/3.; 
0835                 break;
0836             case LargeRLENC::Regions::TransverseMin: //transverse region 2 
0837                 phi_min=lead_jet_center_phi - PI/3.;
0838                 phi_max=lead_jet_center_phi - 2*PI/3.; 
0839         }
0840     }
0841     //make sure the bounds are within the region of 0->2*PI
0842         if(phi_min < 0) phi_min+=2*PI;
0843         if(phi_max < 0) phi_max+=2*PI;
0844         if(phi_min > 2*PI) phi_min+=-2*PI;
0845         if(phi_max > 2*PI) phi_max+=-2*PI;
0846         float holding=std::max(phi_min, phi_max);
0847         float holding_2=std::min(phi_min, phi_max);
0848         phi_min=holding_2;
0849         phi_max=holding;        
0850         std::pair<float, float> phi_lim{phi_min, phi_max}; 
0851         SingleCaloENC(truth, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::TRUTH);
0852         truth.clear();
0853         return;
0854 }
0855 void LargeRLENC::CaloRegion(std::map<std::array<float, 3>, float> emcal, std::map<std::array<float, 3>, float> ihcal, std::map<std::array<float, 3>, float> ohcal, float jetMinpT, std::string variable_of_interest, std::array<float, 3> vertex, float lead_jet_center_phi)
0856 { 
0857     std::map<std::array<float, 3>, float>allcal;
0858     for(auto t:emcal) allcal[t.first]=t.second;
0859     for(auto t:ihcal){
0860         if(allcal.find(t.first) != allcal.end()) allcal[t.first]+=t.second;
0861         else allcal[t.first]=t.second;
0862     }
0863     for(auto t:ohcal){
0864         if(allcal.find(t.first) != allcal.end()) allcal[t.first]+=t.second;
0865         else allcal[t.first]=t.second;
0866     }
0867 
0868     bool transverse=false;
0869     bool energy=true;
0870     float phi_min=0., phi_max=2*PI;
0871     int v_o_i_code=0;
0872     std::map<int, std::pair<float, float>> region_ints;
0873 
0874     if(variable_of_interest.find("T") != std::string::npos) v_o_i_code+=2;
0875     if(variable_of_interest.find("p") != std::string::npos) v_o_i_code+=1; 
0876     switch(v_o_i_code)
0877     {
0878         case 0:
0879             transverse=false;
0880             energy=true;
0881             break;
0882         case 1:
0883             transverse=false;
0884             energy=false;
0885             break;
0886         case 2:
0887             transverse=true;
0888             energy=true;
0889             break;
0890         case 3:
0891             transverse=true;
0892             energy=false;
0893             break;
0894     }
0895     for(int region=0; region< 5; region++)
0896     {
0897         switch(region)
0898         {
0899             case LargeRLENC::Regions::Full: //go over the full calorimeter
0900                 phi_min=0;
0901                 phi_max=2*PI;
0902                 break;
0903             case LargeRLENC::Regions::Towards: //toward region of the detector
0904                 phi_min=lead_jet_center_phi - PI/3.;
0905                 phi_max=lead_jet_center_phi + PI/3.;
0906                 break;
0907             case LargeRLENC::Regions::Away: //away region of the detector
0908                 phi_min=lead_jet_center_phi + 2*PI/3.;
0909                 phi_max=lead_jet_center_phi + 4*PI/3.;
0910                 break;
0911             case LargeRLENC::Regions::TransverseMax: //transverse region1
0912                 phi_min=lead_jet_center_phi + PI/3.;
0913                 phi_max=lead_jet_center_phi + 2*PI/3.; 
0914                 break;
0915             case LargeRLENC::Regions::TransverseMin: //transverse region 2 
0916                 phi_min=lead_jet_center_phi - PI/3.;
0917                 phi_max=lead_jet_center_phi - 2*PI/3.; 
0918         }
0919     //make sure the bounds are within the region of 0->2*PI
0920         if(phi_min < 0) phi_min+=2*PI;
0921         if(phi_max < 0) phi_max+=2*PI;
0922         if(phi_min > 2*PI) phi_min+=-2*PI;
0923         if(phi_max > 2*PI) phi_max+=-2*PI;
0924         float holding=std::max(phi_min, phi_max);
0925         float holding_2=std::min(phi_min, phi_max);
0926         phi_min=holding_2;
0927         phi_max=holding;        
0928         std::pair<float, float> phi_lim{phi_min, phi_max}; 
0929         region_ints[region]=phi_lim;
0930     }
0931 
0932 //  std::vector<std::thread> CaloThreads;
0933 //  CaloThreads.push_back(std::thread(&LargeRLENC::SingleCaloENC, this, emcal, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::EMCAL));
0934     //CaloThreads.push_back(std::thread(&LargeRLENC::SingleCaloENC, this, ihcal, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::IHCAL));
0935 //  CaloThreads.push_back(std::thread(&LargeRLENC::SingleCaloENC, this, ohcal, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::OHCAL));
0936 //  if(! this->pedestalData) 
0937     //  CaloThreads.push_back(std::thread(&LargeRLENC::c
0938     if(!pedestalData) SingleCaloENC(allcal, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::All);
0939 //  for(int ct=0; ct<(int)CaloThreads.size(); ct++) CaloThreads[ct].join();
0940     allcal.clear();
0941         SingleCaloENC(emcal, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::EMCAL);
0942     emcal.clear();
0943         SingleCaloENC(ihcal, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::IHCAL);
0944     ihcal.clear();
0945         SingleCaloENC(ohcal, jetMinpT, vertex, transverse, energy, region_ints, LargeRLENC::Calorimeter::OHCAL);
0946     ohcal.clear();
0947     
0948     return;
0949 }
0950 void LargeRLENC::SingleCaloENC(std::map<std::array<float, 3>, float> cal, float jetMinpT, std::array<float, 3> vertex, bool transverse, bool energy, std::map<int, std::pair<float, float>> phi_edges, LargeRLENC::Calorimeter which_calo)
0951 {
0952     //MethodHistograms* ha=NULL, *hc=NULL;
0953     float base_thresh=thresh_mins[which_calo];
0954     auto i=cal.begin();
0955     std::vector<StrippedDownTower*> strippedCalo;
0956     if(which_calo == LargeRLENC::Calorimeter::EMCAL) base_thresh=base_thresh*16.; //correct for the fact that the emcal has ~16 actual towers going to a single tower object
0957     std::array<float, 4> e_region { 0., 0., 0., 0. };
0958     
0959     
0960     //Processs the calorimeter data into my analysis class
0961     //Do the vertex shift and add it into the output 
0962     std::cout<<"Incoming prefiltered towers : " <<cal.size() <<std::endl;
0963     if(cal.size() == 0 ) return; 
0964     while(i !=cal.end()){
0965         auto j=i->first;
0966         auto j_val=i->second;
0967         if(j_val < base_thresh){
0968             ++i;
0969     //      std::cout<<"Failed with tower energy " <<j_val <<std::endl;
0970             continue;
0971         }
0972         float eta_shift=i->first[0], r=i->first[2];
0973         float x=r*cos(i->first[1]), y=r*sin(i->first[1]);
0974         if(vertex[0]!=0) x+=vertex[0];
0975         if(vertex[1]!=0) y+=vertex[1];
0976         r=std::sqrt(std::pow(x,2) + std::pow(y,2));
0977 
0978         float z=r*sinh(eta_shift);
0979         if(vertex[2]!=0)z+=vertex[2];
0980         eta_shift=asinh(z/r); //did the vertex shift
0981         if(eta_shift != eta_shift) std::cout<<"The eta shift returns nan, r="<<r <<" z=" <<z <<" ratio " <<z/r <<std::endl;
0982         j[0]=eta_shift;
0983         if(vertex[0] != 0 && vertex[1] !=0) j[1]=atan2(y,x);
0984 
0985         j[2]=r; //make sure the positions are properly aligned for all caluclations
0986         StrippedDownTower* t =new StrippedDownTower(base_thresh, (StrippedDownTower::Calorimeter) which_calo); 
0987         for(auto x: phi_edges){
0988             float phi_low=x.second.first;
0989             float phi_up=x.second.second;
0990             if(j[1] >= phi_low && j[1] < phi_up){
0991                 t->tag=(StrippedDownTower::Regions)x.first;
0992                 break;
0993             }
0994             else continue;
0995         }
0996         t->r=r;
0997         t->phi=j[1];
0998         t->eta=j[0];
0999         t->E=j_val;
1000         strippedCalo.push_back(t);
1001         ++i;
1002         //std::cout<<"Towers is now : " <<(int)strippedCalo.size() <<std::endl; 
1003     }
1004     //Now get the energy regions
1005     int n_valid_towers=0;
1006     for(auto* t:strippedCalo)
1007     {
1008         StrippedDownTower::Regions region_tag=t->tag;
1009         if(t->E < t->RegionOutput->threshold) break; //because it is ordered
1010         else{
1011             e_region.at(region_tag)+=t->E;
1012             e_region.at(StrippedDownTower::Regions::Full)+=t->E;
1013             n_valid_towers++;   
1014         }
1015     }
1016     std::cout<<"The thresholds and their number of towers are \n Threshold     Towers" <<std::endl;
1017     if(n_valid_towers <= 0 ) return; 
1018     std::cout<<strippedCalo.at(0)->RegionOutput->threshold <<":    " <<n_valid_towers <<std::endl;
1019     //now I need to actualy run the dataset 
1020     std::vector<std::thread> CalculatingThreads;
1021     //std::cout<<"Calo object has size " <<strippedCalo.size() <<std::endl;
1022     for(int i=0; i<(int)strippedCalo.size()-1; i++)
1023     {
1024         if(i%200 == 0 ) std::cout<<"Prepping tower for threading : " <<i <<std::endl;
1025         StrippedDownTower* t = strippedCalo.at(i);
1026         std::vector<StrippedDownTower> CaloCopy;
1027         for(int j=i+1; j<(int)strippedCalo.size(); j++) CaloCopy.push_back(*(strippedCalo[j])); //avoid double counting
1028 //      CalculateENC(t, CaloCopy, transverse, energy);
1029         CalculatingThreads.push_back(std::thread(&LargeRLENC::CalculateENC, this, t, CaloCopy, transverse, energy));
1030     }
1031     std::cout<<"Performing caluclation" <<std::endl;
1032     for(int k=0; k<(int)CalculatingThreads.size(); k++) CalculatingThreads.at(k).join();
1033     //now need to sum over the relevant towers
1034     CalculatingThreads.clear(); //just trying to find memmory leaks
1035 /*  std::array<std::vector<TowerOutput*>, 5> AllTowOutput;
1036     for(int i=0; i< 5; i++) {
1037         std::vector<TowerOutput*> a;
1038         for(auto t:thresholds){
1039             a.push_back(new TowerOutput(t));
1040         }
1041         AllTowOutput[i]=a;
1042     } //just get it the right range */
1043     std::array<int, 5> n_entries_arr;
1044 //  std::array<std::vector<std::vector<TowerOutput*>>,5> ToMerge; 
1045 //  std::array<std::vector<std::set<float>>,5> ToMergeRl; 
1046 //  std::array<std::vector<std::set<std::array<float,3>>>,5> ToMergeRlRmRs; 
1047     for(int n=0; n<5; n++)
1048         n_entries_arr[n]=0;
1049 //  for(int n=0; n<(int)ToMerge.size(); n++)
1050     //  for(int i=0; i < (int) arr_size; i++) 
1051         //  ToMerge[n].push_back(std::vector<TowerOutput*>());
1052 /*  for(int n=0; n<(int)ToMergeRl.size(); n++)
1053         for(int i=0; i < (int) arr_size; i++) 
1054             ToMergeRl[n].push_back(std::set<float>());
1055     for(int n=0; n< (int)ToMergeRlRmRs.size(); n++)
1056         for(int i=0; i < (int) arr_size; i++) */
1057     //      ToMergeRlRmRs[n].push_back(std::set<std::array<float, 3>>());
1058     std::cout<<"Max elements: " <<strippedCalo.size()<<std::endl;
1059     for(int i=0; i<(int)strippedCalo.size(); i++)
1060     {
1061             if(strippedCalo.at(i)->FullOutput->RL.size() > 0 )
1062                  n_entries_arr[0]++;
1063             //ToMerge[0][j].push_back(strippedCalo.at(i)->FullOutput->at(j));
1064         //  for(auto r:strippedCalo.at(i)->FullOutput->at(j)->RL)ToMergeRl[0][j].insert(r);
1065         //  for(auto r:strippedCalo.at(i)->FullOutput->at(j)->RL_RM_RS)ToMergeRlRmRs[0][j].insert(r);
1066                     
1067         
1068         if(strippedCalo.at(i)->RegionOutput->RL.size() > 0 )
1069              n_entries_arr[(int)strippedCalo.at(i)->tag]++;
1070         //  int tagn=strippedCalo.at(i)->tag;
1071         //  ToMerge[tagn][j].push_back(strippedCalo.at(i)->RegionOutput->at(j));
1072         //  for(auto r:strippedCalo.at(i)->RegionOutput->at(j)->RL)ToMergeRl[tagn][j].insert(r);
1073         //  for(auto r:strippedCalo.at(i)->RegionOutput->at(j)->RL_RM_RS)ToMergeRlRmRs[tagn][j].insert(r);
1074             
1075         
1076     }
1077     std::cout<<"Plotting Energy " <<std::endl;
1078     if(which_calo == LargeRLENC::Calorimeter::TRUTH)
1079     {
1080         for(int region=0; region < 4; region++)
1081         {
1082             Truth_Region_vector[region][0]->E->Fill(e_region[region]);
1083             if(region < 3 )
1084                 Truth_Region_vector[region][0]->N->Fill(n_entries_arr[region]);
1085             else if(region == 3)
1086                 Truth_Region_vector[region][0]->N->Fill(n_entries_arr[region]+n_entries_arr[region+1]);
1087         }
1088         for(auto Tow:strippedCalo)
1089         {
1090             int region=Tow->tag;
1091             //std::cout<<"Threshold index " <<threshold_index<<std::endl;
1092             TowerOutput* TF=Tow->FullOutput;
1093             TowerOutput* TR=Tow->RegionOutput;
1094             int region_shift=region;
1095             if(region_shift > 3 ) region_shift=3;
1096             TF->Normalize(e_region[region_shift]);
1097             TR->Normalize(e_region[region_shift]);
1098             auto Hists=Truth_Region_vector[region_shift][0];
1099             auto FullHists=Truth_Region_vector[0][0];
1100             if((int)TF->RL.size() > 0 ){
1101                 for(int i=0; i<(int)TF->RL.size(); i++){
1102                     FullHists->R->Fill(TF->RL.at(i));
1103                     FullHists->E2C->Fill(TF->RL.at(i), TF->E2C.at(i));
1104                     if((int) TF->E3C.size() <= i ) continue; 
1105                         else FullHists->E3C->Fill(TF->RL.at(i), TF->E3C.at(i));
1106                 }
1107             }
1108             if((int)TR->RL.size() > 0 ){
1109                 for(int i=0; i<(int)TR->RL.size(); i++){
1110                     Hists->R->Fill(TR->RL.at(i));
1111                     Hists->E2C->Fill(TR->RL.at(i), TR->E2C.at(i));
1112                         if((int) TR->E3C.size() <= i ) continue; 
1113                         else Hists->E3C->Fill(TR->RL.at(i), TR->E3C.at(i));
1114                 }
1115             }   
1116         }
1117         strippedCalo.clear();
1118         return;
1119     }
1120     if(which_calo == LargeRLENC::Calorimeter::TRUTH) return;
1121     for(int region=0; region < 4; region++)
1122     {
1123         Region_vector[region][which_calo]->E->Fill(e_region[region]);
1124         if(region < 3 )
1125              Region_vector[region][which_calo]->N->Fill(n_entries_arr[region]);
1126         else if(region == 3)
1127             Region_vector[region][which_calo]->N->Fill(n_entries_arr[region]+n_entries_arr[region+1]);
1128     }
1129     std::cout<<"Plotting all the values" <<std::endl;
1130     
1131     for(auto Tow:strippedCalo)
1132     {
1133         int region=Tow->tag;
1134         TowerOutput* F=Tow->FullOutput;
1135         TowerOutput* R=Tow->RegionOutput;
1136         int region_shift=region;
1137         if(region_shift > 3 ) region_shift=3;
1138         F->Normalize(e_region[region_shift]);
1139         R->Normalize(e_region[region_shift]);
1140         auto Hists=Region_vector[region_shift][which_calo];
1141         auto FullHists=Region_vector[0][which_calo];
1142         if((int)F->RL.size() > 0 ){
1143             for(int i=0; i<(int)F->RL.size(); i++){
1144                 FullHists->R->Fill(F->RL.at(i));
1145                 FullHists->E2C->Fill(F->RL.at(i), F->E2C.at(i));
1146                 if((int) F->E3C.size() <= i ) continue; 
1147                 else FullHists->E3C->Fill(F->RL.at(i), F->E3C.at(i));
1148                 }
1149         }
1150         if((int)R->RL.size() > 0 ){
1151             for(int i=0; i<(int)R->RL.size(); i++){
1152                 Hists->R->Fill(R->RL.at(i));
1153                 Hists->E2C->Fill(R->RL.at(i), R->E2C.at(i));
1154                 if((int) R->E3C.size() <= i ) continue; 
1155                 else Hists->E3C->Fill(R->RL.at(i), R->E3C.at(i));
1156             }
1157         }   
1158     }
1159     std::cout<<"Plotted all the values" <<std::endl;
1160     strippedCalo.clear();
1161 /*  std::vector<std::thread> Merger_thread; //not merging, it was too slow, trying to just fill direct now
1162     for(int i=0; i<5; i++)
1163     {
1164         for(int j=0; j<(int)arr_size; j++)
1165         {
1166             Merger_thread.push_back(std::thread(&LargeRLENC::Merger, this, AllTowOutput[i][j], ToMerge[i][j], ToMergeRl[i][j], ToMergeRlRmRs[i][j]));
1167         }
1168     }
1169     std::cout<<"Merging output" <<std::endl;    
1170     //for(int k=0; k<(int)Merger_thread.size(); k++) Merger_thread.at(k).join(); 
1171     std::cout<<"Saving to plots" <<std::endl;   
1172     //now push the data out to the plots
1173     for(int region_index=0; region_index<(int)AllTowOutput.size(); region_index++)
1174     {
1175         for(int threshold_index=0; threshold_index<(int)AllTowOutput[region_index].size(); threshold_index++)
1176         {
1177             if(threshold_index >= (int)Region_vector.size()) continue;
1178             int region_shift=region_index;
1179         //  std::cout<<"Threshold index is : " <<threshold_index <<" Region: " << region_shift <<std::endl;
1180             if(region_shift >= 3 ) region_shift=3;
1181             AllTowOutput[region_index]->CalculateFlatE3C();
1182             AllTowOutput[region_index]->Normalize(e_region[region_shift]);
1183             auto Hists=Region_vector[region_shift][which_calo];
1184             auto Output=AllTowOutput[region_index];
1185             if(region_index < 3) Hists->N->Fill(n_entries_arr[region_index]);
1186             else if (region_index==3)Hists->N->Fill(n_entries_arr[region_index]+n_entries_arr[region_index + 1]);
1187             else n_entries_arr[region_index]=0;
1188             int tower_number=(int)Output->RL.size();
1189 //          std::cout<<"number of entries to fill is " <<(int)Output->E2C.size() <<std::endl;
1190             if(tower_number > 0){
1191                 for(int rl_index=0; rl_index < tower_number; rl_index++)
1192                 {
1193                     Hists->R->Fill(Output->RL.at(rl_index));
1194                     Hists->E2C->Fill(Output->RL.at(rl_index), Output->E2C.at(rl_index));
1195                     if((int)Output->E3C.size() > rl_index) Hists->E3C->Fill(Output->RL.at(rl_index), Output->E3C.at(rl_index));
1196                 }
1197             }
1198             Hists->E->Fill(e_region[region_shift]);
1199         }
1200     }*/
1201      
1202     return;
1203     
1204 }
1205 void LargeRLENC::Merger(TowerOutput* output, std::vector<TowerOutput*> inputs, std::set<float> RLs, std::set<std::array<float, 3>> Rlmss)
1206 {
1207     std::map<float, float> e2c;
1208     std::map<std::array<float, 3>, float> e3c;
1209     for(auto r:RLs) e2c[r]=0.;
1210     for(auto r:Rlmss) e3c[r]=0.;
1211     for(auto t:inputs)
1212     {
1213         for(int i=0; i <(int) t->RL.size(); i++)
1214         {
1215             e2c[t->RL.at(i)]=t->E2C.at(i);
1216         }
1217         for(int i=0; i< (int) t->RL_RM_RS.size(); i++)
1218         {
1219             e3c[t->RL_RM_RS.at(i)]=t->E3C_full_shape.at(i);
1220         }
1221     }
1222     for(auto er:e2c) output->AddE2CValues(er.second, er.first);
1223     for(auto er:e3c) output->AddE3CValues(er.second, er.first);
1224     
1225     return;
1226 }
1227 void LargeRLENC::CalculateENC(StrippedDownTower* tower1, std::vector<StrippedDownTower> towerSet, bool transverse, bool energy)
1228 {
1229     //this is the threaded object, I could try nesting the threads, but for now I am going to try to run without the nested threads for the 3pt
1230     float thresh=tower1->RegionOutput->threshold;
1231     for(int i = 0; i<(int) towerSet.size(); i++){
1232         //run over all other towers
1233         StrippedDownTower* tower2=&towerSet.at(i); 
1234         float R_L=getR(tower1->eta, tower1->phi, tower2->eta, tower2->phi);
1235         float e2c=0.;
1236         if(energy) e2c=(tower1->E) * (tower2->E);
1237         else e2c=(tower1->E * ptoE) * (tower2->E * ptoE); //use a basic swap of momenta ratio for each calo
1238         if(transverse) e2c=e2c/(cosh(tower1->eta) * cosh(tower2->eta));
1239         float energy1=tower1->E, energy2=tower2->E;
1240         int region1=tower1->tag, region2=tower2->tag;
1241         bool sameRegion = ( region1 == region2);
1242         std::vector<std::pair<std::pair<std::array<float, 3>, std::pair<float, float>>, bool>> e3c_relevant;
1243         if(i != (int) towerSet.size() - 1){
1244             for(int j=i+1; j<(int) towerSet.size(); j++){
1245                 StrippedDownTower* tower3=&towerSet.at(j);
1246                 float R_13=getR(tower1->eta, tower1->phi, tower3->eta, tower3->phi);
1247                 float R_23=getR(tower2->eta, tower2->phi, tower3->eta, tower3->phi);
1248                 if(R_13 > R_L || R_23 > R_L ) continue;
1249                 else{
1250                     float Rs = std::min(R_13, R_23);
1251                     float Rm = std::max(R_13, R_23);
1252                     std::array<float, 3> R {R_L, Rm, Rs}; 
1253                     float energy3=tower3->E;
1254                     float e3c=1;
1255                     if(energy) e3c=e2c*energy3;
1256                     else e3c=e2c* (energy3* ptoE);
1257                     if(transverse) e3c=e3c/cosh(tower3->eta);
1258                     std::pair<float, float> energy_e3 {energy3, e3c};
1259                     std::pair<std::array<float, 3>, std::pair<float, float>> e3c_sing {R, energy_e3};
1260                     bool sameRegion3 = ( region1 == tower3->tag);
1261                     e3c_relevant.push_back(std::make_pair(e3c_sing, sameRegion3));
1262                 }
1263             }
1264         } //caluclate and store the e3c to easily iterate over
1265         if(energy1 > thresh && energy2 > thresh){
1266             tower1->FullOutput->AddE2CValues(e2c, R_L);
1267             //std::cout<<"Adding pair " <<R_L <<", "<<e2c <<std::endl;
1268             if(sameRegion)tower1->RegionOutput->AddE2CValues(e2c, R_L); //this is the two point done noew 
1269             for(auto t3:e3c_relevant)
1270             {
1271                 if(t3.first.second.first > thresh){
1272                     tower1->FullOutput->AddE3CValues(t3.first.second.second, t3.first.first);
1273                     if(t3.second) tower1->RegionOutput->AddE3CValues(t3.first.second.second, t3.first.first);
1274                     }
1275                 }   
1276             }   
1277         
1278     }//caluclate the flattend e3c
1279     tower1->FullOutput->CalculateFlatE3C();
1280     tower1->RegionOutput->CalculateFlatE3C();
1281 
1282     return;
1283 }
1284 void LargeRLENC::JetEventObservablesBuilding(std::array<float, 3> central_tower, std::map<std::array<float, 3>, float> calorimeter_input, std::map<float, float>* Etir_output )
1285 {
1286     //This is the Jet Event observable. Ideally I would seperate this out 
1287     //Input is of the form <tower r, tower eta, tower phi>, tower E (vertex corrected)
1288     //Output is <tower r, tower eta, tower phi>, <R=0.1-1.0, ETir> 
1289 <
1290     for(int r=0; r<3; r++)
1291     {
1292         float Rmax=0.1;
1293         if(r==1) Rmax=0.4;
1294         if(r==2) Rmax=1.0;
1295 
1296         (*Etir_output)[Rmax]=0.;
1297     }
1298     (*Etir_output)[-1.0]=0; //full calorimeter holding 
1299     for(auto tower_energy:calorimeter_input){
1300         if(tower_energy.first==central_tower) continue; //don't double count
1301         float RVal=getR(tower_energy.first[1], tower_energy.first[2], central_tower[1], central_tower[2]);
1302         for(auto et=Etir_output->begin(); et != Etir_output->end(); ++et)
1303         {
1304             if(et->first > RVal)et->second+=tower_energy.second/((float) cosh(tower_energy.first[1]));
1305         }
1306     }
1307     return; 
1308 }
1309 float LargeRLENC::getR(float eta1, float phi1, float eta2, float phi2, bool print)
1310 {
1311     float deta=eta1-eta2; 
1312 
1313     float dphi=std::abs(phi1-phi2);
1314     if(std::abs(dphi) > PI ) dphi+=-PI;
1315     float rsq=std::pow(deta, 2)+ std::pow(dphi, 2);
1316     if(print)std::cout<<" The value for R squared is square of "<<dphi <<" + square of " <<deta <<" = "<<rsq <<std::endl;
1317     return std::sqrt(rsq);
1318 
1319 }
1320 void LargeRLENC::Print(const std::string &what) const
1321 {
1322     TFile* f1=new TFile(output_file_name.c_str(), "RECREATE");
1323 
1324     float pct=(float) n_with_jets / (float) n_evts;
1325     pct=pct*100.; 
1326     TH1F* h_percent=new TH1F("h_pct", "Percentage of events with a non-empty jet container; Percent; N_{files}", 100, -0.5, 99.5);
1327     h_percent->Fill(pct);
1328 //  h_percent->Write();
1329     TDirectory* dir_evt=f1->mkdir("event_categorization");
1330     dir_evt->cd();
1331     std::cout<<"Percentage of events where the jets were present: " <<pct <<std::endl; 
1332 
1333     emcal_occup->Write();
1334     ihcal_occup->Write();
1335     ohcal_occup->Write();
1336     ohcal_rat_h->Write();
1337     ohcal_rat_occup->Write();
1338     ohcal_bad_hits->Write();
1339     bad_occ_em_oh->Write();
1340     bad_occ_em_h->Write();
1341     IE->Write();
1342     badIE->Write();
1343     goodIE->Write();
1344     E_IE->Write();
1345     badE_IE->Write();
1346     goodE_IE->Write();
1347     eventCut->JetCuts->Write();
1348     std::vector<TDirectory*> region_dirs;
1349     std::cout<<__LINE__<<std::endl;
1350     auto FullCal=Region_vector.at(0);
1351     auto TowardRegion=Region_vector.at(1);
1352     auto AwayRegion=Region_vector.at(2);
1353     auto TransverseRegion=Region_vector.at(3);
1354     TDirectory* dir_full=f1->mkdir("Full_Calorimeter");
1355     TDirectory* dir_towrd=f1->mkdir("Towards_Region");
1356     TDirectory* dir_away=f1->mkdir("Away_Region");
1357     TDirectory* dir_trans=f1->mkdir("Transverse_Region");
1358     std::cout<<__LINE__<<std::endl;
1359     std::vector<TDirectory*> tds {dir_full, dir_towrd, dir_away, dir_trans};
1360     if(!isRealData && !pedestalData)
1361     {
1362         for(int i=0; i<(int)Truth_Region_vector.size(); i++)
1363         {
1364             tds.at(i)->cd();
1365             auto hr=Truth_Region_vector.at(i);
1366             for(auto hv:hr){
1367                 hv->R_pt->Write();
1368                 hv->E2C->Write();
1369                 hv->E3C->Write();
1370                 hv->R->Write();
1371                 hv->N->Write();
1372                 hv->E->Write();
1373                 hv->pt->Write();
1374             }
1375         }
1376     }   
1377         dir_full->cd();
1378         std::cout<<__LINE__<<std::endl;
1379         for(auto hv:FullCal){
1380             for( auto h:hv->histsvector){
1381                 try{
1382                      if(h) h->Write();
1383                 }
1384                 catch(std::exception& e){continue;}
1385             }
1386             hv->R_pt->Write();
1387             hv->E2C->Write();
1388             hv->E3C->Write();
1389             hv->R->Write();
1390             hv->N->Write();
1391             hv->E->Write();
1392             hv->pt->Write();
1393         }
1394         dir_trans->cd();
1395         std::cout<<__LINE__<<std::endl;
1396         for(auto hv:TransverseRegion){
1397             for( auto h:hv->histsvector)if(h) h->Write();   
1398             hv->R_pt->Write();
1399             hv->E2C->Write();
1400             hv->E3C->Write();
1401             hv->R->Write();
1402             hv->N->Write();
1403             hv->E->Write();
1404             hv->pt->Write();
1405         }
1406         dir_away->cd();
1407         std::cout<<__LINE__<<std::endl;
1408         for(auto hv:AwayRegion){
1409             for( auto h:hv->histsvector)if(h) h->Write();
1410             hv->R_pt->Write();
1411             hv->E2C->Write();
1412             hv->E3C->Write();
1413             hv->R->Write();
1414             hv->N->Write();
1415             hv->E->Write();
1416             hv->pt->Write();
1417         }
1418         dir_towrd->cd();
1419         std::cout<<__LINE__<<std::endl;
1420         for(auto hv:TowardRegion){  
1421     //      for( auto h:hv->histsvector)if(h) h->Write(); //this stopped working for some reason
1422             hv->R_pt->Write();
1423             hv->E2C->Write();
1424             hv->E3C->Write();
1425             hv->R->Write();
1426             hv->N->Write();
1427             hv->E->Write();
1428             hv->pt->Write();
1429             }
1430         
1431         
1432     std::cout<<__LINE__<<std::endl;
1433     f1->cd();
1434     f1->Write();
1435     std::cout<<__LINE__<<std::endl;
1436 
1437     f1->Close();
1438     return; 
1439 }
[ 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!