Back to home page

sPhenix code displayed by LXR
 
 

    

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

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