Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 //////////////////////////////////////////////////////////////////////////////////////////////////
0002 //                                              //
0003 //                                              //
0004 //      Calorimeter Tower Jets N Point-Energy Correlator Study              //
0005 //      Author: Skadi Grossberndt                           //
0006 //      Date of First Commit: 12 July 2024                      //
0007 //      Date of Last Update:  15 Oct 2024                       //
0008 //      version: v2.0                                   //
0009 //                                              //
0010 //                                              //
0011 //      This project is a study on energy correlators using particle jets and       //
0012 //  jets  definied over calorimeter towers. Feasibility study to evaluate pp prospects  //
0013 //  for sPHENIX in 2024 with minimal tracking and momentum resolution           //
0014 //                                              //
0015 //////////////////////////////////////////////////////////////////////////////////////////////////
0016 
0017 
0018 #include "CalorimeterTowerENC.h"
0019  
0020 CalorimeterTowerENC::CalorimeterTowerENC(int n_run, int n_segment, float jet_min, const std::string &name){
0021     //this is the constructor
0022     n_evts=0;
0023     this->jet_cutoff=jet_min;
0024     std::stringstream run, seg;
0025     run.width(10);
0026     run.fill('0');
0027     run <<n_run;
0028     seg.width(5);
0029     seg.fill('0');
0030     seg << n_segment;
0031     outfilename="Calorimeter_tower_ENC-"+run.str()+"-"+seg.str()+".root";
0032     jethits=new TH2F("jethits", "Location of particle energy deposition from truth jets; #eta; #varphi; N_{hits}", 200, -1.15, 1.05, 200, -3.1416, 3.1415);
0033     calojethits=new TH2F("calojethits", "Location of energy deposition from calo only anti-k_{T} jets; #eta; #varphi; N_{hits}", 200, -1.15, 1.05, 200, -3.1416, 3.1415);
0034     comptotows=new TH2F("comp_to_towers", "Particle energy deposition versus EMCal tower energy per EMCAL tower; E_{particle} [GeV]; E_{emcal} [GeV]; N", 50, -0.5, 49.5, 50, -0.05, 4.95 ); 
0035     number_of_jets=new TH1F("nJ", "Number of Jets per event; N_{jet}; events", 100, -0.5, 99.5);
0036     //cont_pos=new TH1F("Continuous_approx", "Monte Carlo approximation to ideal resolution; R_L; < #sum  pairlikes >", 60, -0.05, 0.85);
0037     EM_energy=new TH1F("emcal_energy", "Total energy in emcal; E [GeV]; events", 200, -0.25, 99.75);
0038     IH_energy=new TH1F("ihcal_energy", "Total energy in ihcal; E [GeV]; events", 200, -0.25, 99.75);
0039     OH_energy=new TH1F("ohcal_energy", "Total energy in ohcal; E [GeV]; events", 200, -0.25, 99.75);
0040     
0041     Particles=new MethodHistograms("Particles");
0042     EMCal=new MethodHistograms("EMCAL");
0043     IHCal=new MethodHistograms("IHCAL");
0044     OHCal=new MethodHistograms("OHCAL");
0045     AllCal=new MethodHistograms("ALLCAL");
0046     EMCalKT=new MethodHistograms("EMCAL_kt");
0047     IHCalKT=new MethodHistograms("IHCAL_kt");
0048     OHCalKT=new MethodHistograms("OHCAL_kt");
0049     AllCalKT=new MethodHistograms("ALLCAL_kt");
0050     histogram_map["parts"]=Particles;
0051     histogram_map["emcal"]=EMCal;
0052     histogram_map["ihcal"]=IHCal;
0053     histogram_map["ohcal"]=OHCal;
0054     histogram_map["emcalkt"]=EMCalKT;
0055     histogram_map["ihcalkt"]=IHCalKT;
0056     histogram_map["ohcalkt"]=OHCalKT;
0057     histogram_map["allcal"]=AllCal;
0058     histogram_map["allcalkt"]=AllCalKT;
0059     for(auto h:this->histogram_map){
0060         std::string typelabel = h.first;
0061         //auto hists=h.second;
0062         std::cout<<"Making the histograms for " <<h.first <<std::endl;
0063 //      for(auto h1:hists->histsvector) std::cout<<"Have a histogram with that has name " <<h1->GetName()  <<std::endl;
0064     }
0065 }
0066 
0067 int CalorimeterTowerENC::Init(PHCompositeNode *topNode)
0068 {   
0069     //book histograms and TTrees
0070     histogram_map["parts"]->E->GetName();   
0071     return 0; 
0072 }
0073 
0074 float CalorimeterTowerENC::getPt(PHG4Particle* p)
0075 {
0076     float pt=0;
0077     pt=pow(p->get_px(), 2) + pow(p->get_py(), 2);
0078     pt=sqrt(pt);
0079     return pt;
0080 }
0081 
0082 float CalorimeterTowerENC::getR(std::pair<float, float> p1, std::pair<float,float> p2){
0083     float eta_dist=p1.first - p2.first;
0084     float phi_dist=p1.second - p2.second;
0085     phi_dist=abs(phi_dist);
0086     if(phi_dist > PI) phi_dist=2*PI-phi_dist;
0087     float R=sqrt(pow(eta_dist,2) + pow(phi_dist, 2));
0088     return R;
0089 }
0090 void CalorimeterTowerENC::GetENCCalo(std::map<int, float> allcal, RawTowerGeomContainer_Cylinderv1* geom, TowerInfoContainer* data, MethodHistograms* histograms, float jete_tot)
0091 {
0092     geom->set_calorimeter_id(RawTowerDefs::HCALOUT);
0093     int ntow=0; 
0094     float jete=0.0;
0095     for(auto e:allcal) jete+=e.second;
0096     for(auto e:allcal)
0097     {
0098         if(e.second >= 0.01) ntow++; //10MeV threshold
0099         if(!geom) continue;
0100         auto key_1=data->encode_key(e.first);
0101 //      auto tower_2 = data->get_tower_at_channel(e.first);
0102         int phibin_1=data->getTowerPhiBin(key_1);
0103         int etabin_1=data->getTowerEtaBin(key_1);
0104         if(etabin_1 < 0 || phibin_1 < 0 ||  phibin_1 >= geom->get_phibins() || etabin_1 >= geom->get_etabins() ) continue;
0105         float phi_center_1 = geom->get_phicenter(phibin_1);
0106         float eta_center_1 = geom->get_etacenter(etabin_1);
0107         for(auto f:allcal)
0108         {
0109             if(e.first >= f.first) continue;
0110             auto key_2 = data->encode_key(f.first);
0111 //          auto tower_2 = data->get_tower_at_channel(f.first);
0112             int phibin_2 = data->getTowerPhiBin(key_2);
0113             int etabin_2 = data->getTowerEtaBin(key_2);
0114             float phi_center_2 = geom->get_phicenter(phibin_2);
0115             float eta_center_2 = geom->get_etacenter(etabin_2);
0116             if(etabin_2 < 0 || phibin_2 < 0 ||  phibin_2 >= geom->get_phibins() || etabin_2 >= geom->get_etabins() ) continue;
0117             std::pair<float, float> tower_center_1 {eta_center_1, phi_center_1}, tower_center_2 {eta_center_2, phi_center_2};
0118             float R_12=getR(tower_center_1, tower_center_2);
0119             histograms->R_geom->Fill(R_12);
0120             float jet2=pow(jete,2);
0121             float jet2t=pow(jete_tot,2);
0122             if (n_evts < 2) std::cout<<"The jet energy divsor is " <<jet2<<std::endl;
0123             float e2c=e.second*f.second / jet2;
0124             float e2c_all=e.second*f.second /jet2t; 
0125             if(abs(e.second) <= 0.01 )continue; //10 MeV threshold for consideration
0126             if(abs(f.second) <= 0.01) continue; //10 MeV threshold for consideration
0127             if(e2c != 0 ) histograms->R->Fill(R_12);
0128             e2c=e2c/((float) Nj);
0129             e2c_all=e2c_all/((float) Nj);
0130             std::cout<<"The 2 point energy correlator is " <<e2c <<std::endl;
0131             if(!std::isnormal(e2c) ) continue;
0132             histograms->E2C->Fill(R_12, e2c/*/(float) Nj*/);
0133             for(auto g:allcal)
0134             {
0135                 if( f.first >= g.first || e.first >= g.first) continue;
0136                 if(g.first < 0 /*|| k >= (int) data->size()*/ ) continue;
0137                 auto key_3 = data->encode_key(g.first);
0138                 int phibin_3 = data->getTowerPhiBin(key_3);
0139                 int etabin_3 = data->getTowerEtaBin(key_3);
0140                 if(etabin_3 < 0 || phibin_3 < 0 ||  phibin_3 >= geom->get_phibins() || etabin_3 >= geom->get_etabins() ) continue;
0141                 float phi_center_3 = geom->get_phicenter(phibin_3);
0142                 float eta_center_3 = geom->get_etacenter(etabin_3);
0143                 std::pair<float, float> tower_center_3 {eta_center_3, phi_center_3};
0144                 float R_13=getR(tower_center_1, tower_center_3);
0145                 float R_23=getR(tower_center_2, tower_center_3);
0146                 float e3c=e.second*f.second*g.second;
0147                 float jete3=pow(jete, 3);
0148                 float jete3t=pow(jete_tot, 3);
0149                 float e3c_all=e3c/jete3t;
0150                 e3c=e3c/jete3;
0151                 float R_L = std::max(R_12, R_13);
0152                 R_L=std::max(R_L, R_23);
0153                 e3c=e3c/((float) Nj);
0154                 e3c_all=e3c_all/((float) Nj);
0155                 if(!std::isnormal(e3c)) continue;
0156                 histograms->E3C->Fill(R_L, e3c/*/(float) Nj*/);
0157                 if(std::isnormal(e3c_all)) histograms->E3C_pt->Fill(R_L, e3c_all);
0158             }
0159         }
0160     }
0161     histograms->N->Fill(ntow);
0162     histograms->E->Fill(jete);
0163 }
0164         
0165 void CalorimeterTowerENC::GetENCCalo(PHCompositeNode* topNode, std::unordered_set<int> tower_set, TowerInfoContainer* data, RawTowerGeomContainer_Cylinderv1* geom, RawTowerDefs::CalorimeterId calo, float jete, std::string caloh, int npt, float jete_tot, std::map<int, float>* allcal_e)
0166 {
0167     geom->set_calorimeter_id(calo);
0168     MethodHistograms* histograms=histogram_map.at(caloh);
0169     histograms->E->Fill(jete);
0170     int ntow=0;
0171     std::cout<<"taking data on a jet size " <<tower_set.size() <<std::endl;
0172     for(auto i:tower_set){
0173         if(!geom) continue;
0174         if(i < 0/* || i >= (int) data->size()*/ ) continue;
0175         auto key_1 = data->encode_key(i);
0176         auto tower_1 = data->get_tower_at_channel(i);
0177         int phibin_1 = data->getTowerPhiBin(key_1);
0178         int etabin_1 = data->getTowerEtaBin(key_1);
0179         if(!tower_1) continue;
0180     //  std::cout<<"Phibin: " <<phibin_1 <<" Etabin: " <<etabin_1 <<std::endl;
0181         if(etabin_1 < 0 || phibin_1 < 0 ||  phibin_1 >= geom->get_phibins() || etabin_1 >= geom->get_etabins() ) continue;
0182         float phi_center_1 = geom->get_phicenter(phibin_1);
0183         float eta_center_1 = geom->get_etacenter(etabin_1);
0184         float energy_1= tower_1->get_energy(); 
0185         int ih=i;
0186         if(caloh.find("emcal") != std::string::npos && emcal_lookup.find(i) !=emcal_lookup.end()) ih=emcal_lookup[i];
0187         if(allcal_e->find(ih) == allcal_e->end()) (*allcal_e)[ih]=energy_1;
0188         else allcal_e->at(i)+=energy_1;
0189          ntow++;
0190         //there is probably a smart way to collect the group of sizes n, but I can really just do it by hand for rn 
0191         for(auto j:tower_set){
0192             if (i >= j) continue;
0193             if(j < 0 /*|| j >= (int) data->size()*/ ) continue;
0194             auto key_2 = data->encode_key(j);
0195             auto tower_2 = data->get_tower_at_channel(j);
0196             int phibin_2 = data->getTowerPhiBin(key_2);
0197             int etabin_2 = data->getTowerEtaBin(key_2);
0198             float phi_center_2 = geom->get_phicenter(phibin_2);
0199             float eta_center_2 = geom->get_etacenter(etabin_2);
0200             if(etabin_2 < 0 || phibin_2 < 0 ||  phibin_2 >= geom->get_phibins() || etabin_2 >= geom->get_etabins() ) continue;
0201             float energy_2 = tower_2->get_energy();
0202             std::pair<float, float> tower_center_1 {eta_center_1, phi_center_1}, tower_center_2 {eta_center_2, phi_center_2};
0203             float R_12=getR(tower_center_1, tower_center_2);
0204             histograms->R_geom->Fill(R_12);
0205             if(n_evts < 2 && energy_1 >0 && energy_2 >0 )   std::cout<<"Seperation is R="<<R_12<<std::endl;
0206             float jet2=pow(jete,2);
0207             float jet2t=pow(jete_tot,2);
0208             if (n_evts < 2&& energy_1 >0 && energy_2 >0 ) std::cout<<"The jet energy divsor is " <<jet2 <<" Total is " <<jet2t <<std::endl;
0209             float e2c=energy_1*energy_2 / jet2;
0210             float e2c_all=energy_1*energy_2 /jet2t; 
0211             if(abs(energy_1) <= 0.001 )continue; //1 MeV threshold for consideration
0212             if(abs(energy_2) <= 0.001) continue; //1 MeV threshold for consideration
0213             if(e2c != 0 ) histograms->R->Fill(R_12);
0214             e2c=e2c/((float) Nj);
0215             e2c_all=e2c_all/((float) Nj);
0216             if(!std::isnormal(e2c) ) continue;
0217             histograms->E2C->Fill(R_12, e2c/*/(float) Nj*/);
0218             if(std::isnormal(e2c_all)) histograms->E2C_pt->Fill(R_12, e2c_all);
0219             if( npt >= 3){
0220                 for( auto k:tower_set){
0221                     if( j >= k || i >= k) continue;
0222                     if(k < 0 /*|| k >= (int) data->size()*/ ) continue;
0223                     auto key_3 = data->encode_key(k);
0224                     auto tower_3 = data->get_tower_at_channel(k);
0225                     int phibin_3 = data->getTowerPhiBin(key_3);
0226                     int etabin_3 = data->getTowerEtaBin(key_3);
0227                     if(etabin_3 < 0 || phibin_3 < 0 ||  phibin_3 >= geom->get_phibins() || etabin_3 >= geom->get_etabins() ) continue;
0228                     float phi_center_3 = geom->get_phicenter(phibin_3);
0229                     float eta_center_3 = geom->get_etacenter(etabin_3);
0230                     float energy_3 = tower_3->get_energy();
0231                     std::pair<float, float> tower_center_3 {eta_center_3, phi_center_3};
0232                     float R_13=getR(tower_center_1, tower_center_3);
0233                     float R_23=getR(tower_center_2, tower_center_3);
0234                     if(abs(energy_3) <= 0.001)continue;
0235                     float e3c=energy_3*energy_2*energy_1;
0236                     float jete3=pow(jete, 3);
0237                     float jete3t=pow(jete_tot, 3);
0238                     float e3c_all=e3c/jete3t;
0239                     e3c=e3c/jete3;
0240                     float R_L = std::max(R_12, R_13);
0241                     R_L=std::max(R_L, R_23);
0242                     e3c=e3c/((float) Nj);
0243                     e3c_all=e3c_all/((float) Nj);
0244                     if(!std::isnormal(e3c)) continue;
0245                     histograms->E3C->Fill(R_L, e3c/*/(float) Nj*/);
0246                     if(std::isnormal(e3c_all)) histograms->E3C_pt->Fill(R_L, e3c_all);
0247                 }
0248             }
0249         }
0250     }
0251     histograms->N->Fill(ntow);
0252 }
0253 void CalorimeterTowerENC::GetE2C(PHCompositeNode* topNode, std::unordered_set<int> em, std::unordered_set<int> ih, std::unordered_set<int> oh, float jete)
0254 {
0255     //this is the processor that allows for read out of the tower energies given the set of towers in the phi-eta plane
0256     auto emcal_geom =  findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC");
0257     auto emcal_tower_energy =  findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0258     auto ihcal_geom= findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0259     auto ihcal_tower_energy= findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0260     auto ohcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0261     auto ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0262     float jete_em = 0, jete_ih=0, jete_oh=0;
0263     for(auto i:em)
0264     {
0265         try{
0266             auto tower_e=emcal_tower_energy->get_tower_at_channel(i);
0267             jete_em+=tower_e->get_energy();
0268         }
0269         catch(std::exception& e) { continue; }
0270     }
0271     for(auto i:ih)
0272     {
0273         try{
0274             auto tower_e=ihcal_tower_energy->get_tower_at_channel(i);
0275             jete_ih+=tower_e->get_energy();
0276         }
0277         catch(std::exception& e) { continue; }
0278     }
0279     for(auto i:oh)
0280     {
0281         try{
0282             auto tower_e=ohcal_tower_energy->get_tower_at_channel(i);
0283             jete_oh+=tower_e->get_energy();
0284         }
0285         catch(std::exception& e) { continue; }
0286     }
0287     std::map<int, float> allcal_e;
0288     histogram_map["emcal"]->E->Fill(jete_em);
0289     histogram_map["ihcal"]->E->Fill(jete_ih);
0290     histogram_map["ohcal"]->E->Fill(jete_oh);
0291     histogram_map["emcal"]->N->Fill(em.size());
0292     histogram_map["ihcal"]->N->Fill(ih.size());
0293     histogram_map["ohcal"]->N->Fill(oh.size());
0294     GetENCCalo(topNode,em, emcal_tower_energy, emcal_geom, RawTowerDefs::CEMC,    jete_em, "emcal", 2, jete, &allcal_e);
0295     GetENCCalo(topNode,ih, ihcal_tower_energy, ihcal_geom, RawTowerDefs::HCALIN,  jete_ih, "ihcal", 2, jete, &allcal_e);
0296     GetENCCalo(topNode,oh, ohcal_tower_energy, ohcal_geom, RawTowerDefs::HCALOUT, jete_oh, "ohcal", 2, jete, &allcal_e);
0297     return;
0298 }
0299 void CalorimeterTowerENC::GetE3C(PHCompositeNode* topNode, std::unordered_set<int> em, std::unordered_set<int> ih, std::unordered_set<int> oh, std::map<int, float>* emtowere, bool is_kt, float jete)
0300 {
0301     //this is the processor that allows for read out of the tower energies given the set of towers in the phi-eta plane
0302     auto emcal_geom =  findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC");
0303     auto emcal_tower_energy =  findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0304     auto ihcal_geom= findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0305     auto ihcal_tower_energy= findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0306     auto ohcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0307     auto ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0308     float jete_em = 0, jete_ih=0, jete_oh=0;
0309     for(int n=0; n<(int) ohcal_tower_energy->size(); n++){
0310          jete_oh+=ohcal_tower_energy->get_tower_at_channel(n)->get_energy();
0311     //  std::cout<<"For tower at channel " <<n <<" the energy is : " <<ohcal_tower_energy->get_tower_at_channel(n)->get_energy() <<std::endl;
0312     }
0313 //  if (n_evts < 10 ) std::cout<<"Total energy found in the ohcal is " <<jete_oh <<std::endl;
0314     jete_oh=0;
0315     for(auto i:em)
0316     {
0317         try{
0318             if( i <= 0 ) continue;
0319 //          if(n_evts <10 ) std::cout<<"The tower number in use is " <<i <<std::endl;
0320             auto tower_e=emcal_tower_energy->get_tower_at_channel(i);
0321             if(tower_e){
0322                 auto et=tower_e->get_energy();
0323                 /*if(et == 0 && i > 1) et+=emcal_tower_energy->get_tower_at_channel(i-1)->get_energy();
0324                 if(et == 0 && i > 1) et+=emcal_tower_energy->get_tower_at_channel(i+1)->get_energy();*/
0325                 jete_em+=et;
0326 //              if(n_evts < 10 ) std::cout<<"The added energy is " <<et <<" total energy " <<jete_em <<std::endl;; 
0327                 (*emtowere)[i]+=tower_e->get_energy();
0328             }
0329         }
0330         catch(std::exception& e) { continue; }
0331     }
0332     for(auto i:ih)
0333     {
0334         try{
0335             if( i <= 0 ) continue;
0336             auto tower_e=ihcal_tower_energy->get_tower_at_channel(i);
0337             if(tower_e) jete_ih+=tower_e->get_energy();
0338         }
0339         catch(std::exception& e) { continue; }
0340     }
0341     for(auto i:oh)
0342     {
0343         try{
0344             if( i <= 0 ) continue;
0345             auto tower_e=ohcal_tower_energy->get_tower_at_channel(i);
0346     //      std::cout<<"Jet Tower: We get the energy for the tower at channel " <<i <<" as " <<tower_e->get_energy() <<std::endl;
0347             if(tower_e) jete_oh+=tower_e->get_energy();
0348         }
0349         catch(std::exception& e) { continue; }
0350     }
0351     //std::cout<<"Checked " <<em.size() <<" towers in the outer hcal \n Got an energy in the ohcal part of the jet of " <<jete_oh <<std::endl; 
0352 //  histogram_map["ohcal"]->E->Fill(jete_oh);
0353     std::vector<std::thread> calo_thread; 
0354     std::string emcal_label="emcal", ihcal_label="ihcal", ohcal_label="ohcal", all="allcal";
0355     if(!is_kt){
0356         emcal_label+="kt";
0357         ihcal_label+="kt";
0358         ohcal_label+="kt";
0359         all+="kt";
0360     }
0361     std::map<int, float> allcal, allcal_em, allcal_ih, allcal_oh;
0362     /*calo_thread.push_back(std::thread(&CalorimeterTowerENC::*/GetENCCalo(/*, this,*/topNode,em, emcal_tower_energy, emcal_geom, RawTowerDefs::CEMC, jete_em, emcal_label, 3, jete, &allcal_em)/*)*/;
0363 /*  calo_thread.push_back(std::thread(&CalorimeterTowerENC::*/GetENCCalo(/*, this,*/ topNode,ih, ihcal_tower_energy, ihcal_geom, RawTowerDefs::HCALIN,  jete_ih, ihcal_label, 3, jete, &allcal_ih)/*)*/;
0364 /*  calo_thread.push_back(std::thread(&CalorimeterTowerENC::*/GetENCCalo(/*, this, */topNode,oh, ohcal_tower_energy, ohcal_geom, RawTowerDefs::HCALOUT, jete_oh, ohcal_label, 3, jete, &allcal_oh/*)*/);
0365 //  for(int t=0; t<(int)calo_thread.size(); t++) calo_thread.at(t).join();
0366     for(auto e:allcal_em)
0367     {
0368         if(allcal.find(e.first)==allcal.end()) allcal[e.first]=e.second;
0369         else allcal[e.first]+=e.second;
0370     }
0371     for(auto e:allcal_ih)
0372     {
0373         if(allcal.find(e.first)==allcal.end()) allcal[e.first]=e.second;
0374         else allcal[e.first]+=e.second;
0375     }
0376     for(auto e:allcal_oh)
0377     {
0378         if(allcal.find(e.first)==allcal.end()) allcal[e.first]=e.second;
0379         else allcal[e.first]+=e.second;
0380     }
0381     //MethodHistograms* hs=histogram_map[all];
0382     //GetENCCalo(allcal, ohcal_geom, ohcal_tower_energy, hs, jete);
0383     return;
0384 }
0385 
0386 void CalorimeterTowerENC::GetE2C(PHCompositeNode* topNode, std::map<PHG4Particle*, std::pair<float, float>> jet)
0387 {
0388     //this is the processor that allows for particle read of jets phi-eta plane
0389     float jet_total_e=0;
0390     for(auto p1:jet) jet_total_e+=p1.first->get_e();
0391     histogram_map["parts"]->E->Fill(jet_total_e);
0392     histogram_map["parts"]->N->Fill(jet.size());
0393     for(auto p1it=jet.begin(); p1it != jet.end(); ++p1it )
0394     {
0395         auto p1=(*p1it);
0396         for(auto p2it=p1it; p2it != jet.end(); ++p2it){
0397             auto p2=(*p2it);
0398             if(p1it == p2it) continue;
0399             double pe1=p1.first->get_e();
0400             double pe2=p2.first->get_e();
0401             float R12=getR(p1.second, p2.second);
0402             histogram_map["parts"]->R->Fill(R12);
0403             float e2c=pe1*pe2/jet_total_e;
0404             e2c=e2c/(float)Nj;
0405             histogram_map["parts"]->E2C->Fill(R12, e2c);
0406         }
0407     }
0408     return;
0409             
0410 }
0411 
0412 void CalorimeterTowerENC::GetE3C(PHCompositeNode* topNode, std::map<PHG4Particle*, std::pair<float, float>> jet)
0413 {
0414     //this is the processor that allows for particle read of jet in the phi-eta plane
0415     float jet_total_e=0;
0416     for(auto p1:jet) jet_total_e+=p1.first->get_e();
0417 //  std::cout<<"The histogram map has size " <<histogram_map.size() <<std::endl;
0418 //  for(auto h:this->histogram_map) for(auto h1:h.second->histsvector) std::cout<<"The histogram map for " <<h.first <<" has a histogram with name " <<h1->GetName() <<std::endl;
0419     this->histogram_map["parts"]->E->Fill(jet_total_e);
0420     this->histogram_map["parts"]->N->Fill(jet.size());
0421     for(auto p1it=jet.begin(); p1it != jet.end(); ++p1it )
0422     {
0423         auto p1=(*p1it);
0424         for(auto p2it=p1it; p2it != jet.end(); ++p2it){
0425             auto p2=(*p2it);
0426             if(p1it == p2it) continue;
0427             double pe1=p1.first->get_e();
0428             double pe2=p2.first->get_e();
0429             float R12=getR(p1.second, p2.second);
0430             histogram_map["parts"]->R->Fill(R12);
0431             float e2c=pe1*pe2/(float) pow(jet_total_e, 2);
0432             e2c=e2c/(float)Nj;
0433             histogram_map["parts"]->E2C->Fill(R12, e2c);
0434             for(auto p3it=p2it; p3it != jet.end() ; ++p3it){
0435                 auto p3=(*p3it);
0436                 if(p1it == p3it || p2it == p3it) continue;
0437                 double pe3=p3.first->get_e();
0438                 float R13=getR(p1.second, p3.second);
0439                 float R23=getR(p2.second, p3.second);
0440                 float e3c=e2c*pe3/(float)jet_total_e;
0441                 float R_L=std::max(R12, R23);
0442                 R_L=std::max(R_L, R13);
0443                 histogram_map["parts"]->E3C->Fill(R_L, e3c);
0444             }
0445         }
0446     }
0447     return;
0448 }
0449     
0450 std::pair<std::map<float, std::map<float, int>>, std::pair<float, float>> CalorimeterTowerENC::GetTowerMaps(RawTowerGeomContainer_Cylinderv1* geom, RawTowerDefs::CalorimeterId caloid, TowerInfoContainer* calokey)
0451 {
0452     //ge thte geometry map for the towers in the first event to may it easy to serarc
0453     //the idea is put the menan in and the associate a key with it 
0454     geom->set_calorimeter_id(caloid);
0455     std::pair<float, float> deltas{0,0};
0456     std::map<float, std::map<float, int>> bins;
0457     std::map<int, std::pair<float, float>> gs;
0458     geom->set_calorimeter_id(caloid);
0459     for(int i = 0; i<(int) calokey->size(); i++){
0460         try{
0461     //      std::cout<<"The channel number is " <<i <<std::endl;
0462             int key=calokey->encode_key(i);
0463             //int tower=calokey->get_tower_at_channel(key);
0464             int phibin=calokey->getTowerPhiBin(key);
0465             int etabin=calokey->getTowerEtaBin(key);
0466             float eta=geom->get_etacenter(etabin);
0467             float phi=geom->get_phicenter(phibin);
0468             std::pair g { eta, phi};
0469             if( deltas.first == 0){
0470                 std::pair <double, double> etabounds=geom->get_etabounds(etabin);
0471                 deltas.first=abs(etabounds.first - etabounds.second)/2.;
0472                 std::pair <double, double> phibounds=geom->get_phibounds(etabin);
0473                 deltas.second=abs(phibounds.first - phibounds.second)/2.;
0474             }
0475             gs[i]=g;
0476         }
0477         catch(std::exception& e){ 
0478             std::pair<float, float> g {-5.0, -20.0 };
0479             gs[i]=g;
0480         } 
0481     }
0482     int i=0;
0483     for(auto g1: gs){
0484         auto g=g1.second;
0485         if(g.first == -5  && i > 1){
0486             if((int)(gs.at(i-1).first*4.4/deltas.first) % 8 != 7 && (int)((3.1416/deltas.second)*gs.at(i-1).second) % 2 == 1  ){
0487                 g.first=gs.at(i-1).first+2*deltas.first; 
0488                 g.second=gs.at(i-1).second;
0489             }
0490             else if ((int)((3.1416/deltas.second)*gs.at(i-1).second) % 2 != 1 ){
0491                 g.first=gs.at(i-1).first;
0492                 g.second=gs.at(i-1).second + 2*deltas.second;
0493             }
0494             else if ((int)((4.4/deltas.first)*gs.at(i-1).first) % 8 == 7  && (int)((3.1416/deltas.second)*gs.at(i-1).second) % 2 == 1 && (4.4/deltas.first)*gs.at(i-1).first != 1.1 - deltas.first){
0495                 g.first=gs.at(i-1).first+2*deltas.first;
0496                 g.second=gs.at(i-1).second-16*deltas.second;
0497             }
0498             else if( (int) ((4.4/deltas.first)*gs.at(i-1).first) == 1.1 - deltas.first && (int)((3.1416/deltas.second)*gs.at(i-1).second) % 2 == 1 ) {
0499                 gs.at(i-1).first = -1.1 + deltas.first; 
0500                 gs.at(i-1).second = gs.at(i-1).second + 2*deltas.second;
0501             }
0502         }
0503         bins[g.first][g.second]=g1.first;
0504         i++;
0505     }
0506     return std::make_pair(bins, deltas);
0507 }
0508 int CalorimeterTowerENC::GetTowerNumber(std::pair<float, float> part, std::map<float, std::map< float, int>> Calomap, std::pair<float, float> delta)
0509 {
0510     int keynumber=-1;
0511     part.second+=PI;
0512 //  std::cout<<"The particle is centered around eta = " <<part.first <<" phi = " <<part.second <<std::endl;
0513     if (abs(part.first) > 1.2 ) return keynumber;
0514     for(auto t:Calomap)
0515     {
0516         if( part.first < t.first + delta.first && part.first >= t.first - delta.first){
0517             //checking the eta bounds first
0518             for (auto tp:t.second){
0519                 if(part.second < tp.first + delta.second && part.second >= tp.first - delta.second){
0520     //              if(n_evts < 10 ) std::cout<<"Found a good bin centered in eta at " <<t.first <<" in phi at "<<tp.first <<" with bin number " <<tp.second <<std::endl;
0521                     keynumber=tp.second;
0522                     break;
0523                 }
0524             }
0525         }
0526     }
0527 //  if(n_evts < 10 ) std::cout<<" The particle is at eta= " <<part.first <<" and phi = " <<part.second <<std::endl;
0528     return keynumber;
0529 }
0530 int CalorimeterTowerENC::RecordHits(PHCompositeNode* topNode, Jet* truth_jet, std::vector<std::unordered_set<int>> kt_towers){
0531     //record where the particles are located and which towers are seeing hits
0532     std::set<std::pair<float, float>> particle_coords; //This will store the location of each particle which then I can use tho capture the relevant towe
0533     std::unordered_set<PHG4Particle*> jetparticles;
0534     std::set<std::pair<float, float>> circle_jet;
0535     std::map<int, float> jettowenergy, emtowerenergy, emtowerenergykt;
0536     auto _truthinfo=findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0537 //  auto _truthhits=findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BH_1");
0538     auto particles=_truthinfo->GetMap();    
0539     //std::cout<<"Indicies for the truth map are ";
0540     //for(auto p:particles) std::cout<<p.first <<std::endl; 
0541     for( auto& iter:truth_jet->get_comp_vec()){
0542         Jet::SRC source = iter.first;
0543         unsigned int idx = iter.second;
0544         if( source != Jet::PARTICLE){
0545             std::cout<<"This jet has data that is not a particle" <<std::endl;
0546             return -1;
0547         }
0548             if (particles.find(idx) == particles.end()) std::cout<<"Index is not in map for index " <<idx <<std::endl;  
0549             else { 
0550                 PHG4Particle* truth_part = _truthinfo->GetParticle(idx);
0551             //  PHG4Shower* truth_shower = _truthinfo->GetShower(idx);
0552                 jetparticles.insert(truth_part);
0553             //  auto showerhits=truth_shower->g4Hit_ids()
0554             }
0555         //jethits.insert(truth_hit);
0556     }
0557     std::map<PHG4Particle*, std::pair<float, float>> jet_particle_map;
0558     float axis_phi=truth_jet->get_phi(), axis_eta=truth_jet->get_eta();
0559     float edge_R=0;
0560     float jet_energy=truth_jet->get_e();
0561     int rdivbins=Particles->R_geom->GetNbinsX();    
0562     float rbinsides=Particles->R_geom->GetBinWidth(5);
0563     for(auto p:jetparticles){
0564         std::pair<float, float> particle_coord; 
0565         particle_coord.first = atanh(p->get_pz()/sqrt(pow(p->get_py(),2) + pow(p->get_px(),2) + pow(p->get_pz(), 2)));
0566         particle_coord.second = atan2(p->get_py(), p->get_px());
0567         float R=sqrt(pow(particle_coord.second-axis_phi, 2) + pow(particle_coord.first - axis_eta, 2));
0568         if(edge_R < R) edge_R=R;
0569     }
0570     for(auto p:jetparticles){
0571            //take in the identified particles that make up the subjet and then we overap that in the towers
0572         std::pair<float, float> particle_coord; 
0573         particle_coord.first = atanh(p->get_pz()/sqrt(pow(p->get_py(),2) + pow(p->get_px(),2) + pow(p->get_pz(), 2)));
0574         particle_coord.second = atan2(p->get_py(), p->get_px());
0575         //float angle_off=acos((particle_coord.second-axis_phi)/(particle_coord.first-axis_eta)); //this is the anglular offset of the particle circle
0576         particle_coords.insert(particle_coord); 
0577         jet_particle_map[p]=particle_coord;
0578         jethits->Fill(particle_coord.first, particle_coord.second);
0579 //      if ( n_evts < 10 ) std::cout<<"The particle we are looking for is located at eta = " <<particle_coord.first 
0580 //              <<" phi = " << particle_coord.second <<std::endl;
0581         int towernumberemcal=GetTowerNumber(particle_coord, EMCALMAP.first, EMCALMAP.second);
0582         if(jettowenergy.find(towernumberemcal) == jettowenergy.end())jettowenergy[towernumberemcal]=0; 
0583         jettowenergy[towernumberemcal]+=p->get_e();
0584         float R=sqrt(pow(particle_coord.second-axis_phi, 2) + pow(particle_coord.first - axis_eta, 2));
0585         float ptr=sqrt(pow(p->get_px(),2)+ pow(p->get_py(), 2))/(float) truth_jet->get_pt();
0586         histogram_map["parts"]->R_pt->Fill(R/(float)edge_R, ptr);
0587         float intersectionbase=(pow(edge_R, 2) + pow(R, 2))/(2*edge_R*R);
0588         for(int i=0; i < rdivbins; i++)
0589         {
0590             //this fills the "continous" case approximation
0591             float R_L=i*rbinsides;
0592             float intersection=intersectionbase+pow(R_L,2)/(2*edge_R*R);
0593             float circumfrence_included=0.0;
0594             if(intersection >= 1) {
0595                 //this is the case of all of the ring within the jet
0596                 circumfrence_included=2*3.14159;
0597             }
0598             else if(intersection <= -1){
0599                 circumfrence_included=0.;
0600             }
0601             else{
0602                 try{
0603                     circumfrence_included=2*acos(intersection);
0604                 }
0605                 catch(std::exception& e) {};
0606                 if( R_L+R > edge_R && circumfrence_included > 3.14159) circumfrence_included = 2*3.14159-circumfrence_included;
0607                 else if(R_L +R < edge_R && circumfrence_included < 3.14159 ) circumfrence_included = 2*3.14159 - circumfrence_included;
0608             }
0609             circumfrence_included=R_L*circumfrence_included;
0610             if(!std::isnormal(circumfrence_included)) continue;
0611             Particles->R_geom->Fill(R_L, circumfrence_included);
0612         }
0613     }
0614     //find the center of the jet and go to R=0.4 to try to correct for particles not being in the right place
0615     //this is a fast way to get around the G4Hit version
0616     std::pair<float, float> jet_center {0,0};
0617     for(auto p:jet_particle_map){
0618         jet_center.first+=p.second.first;
0619         jet_center.second+=p.second.second;
0620     }
0621     jet_center.first=jet_center.first/(float)jet_particle_map.size();
0622     jet_center.second=jet_center.second/(float)jet_particle_map.size();
0623     int netabin=1+0.4/(float) EMCALMAP.second.first, nphibin=1+0.4/(float) EMCALMAP.second.second;
0624     for(int etabin=0; etabin< netabin; etabin++)
0625     {
0626         for(int phibin=0; phibin<nphibin; phibin++)
0627         {
0628             float deta=EMCALMAP.second.first*etabin;
0629             float dphi=EMCALMAP.second.second*phibin;
0630             float R=sqrt(pow(deta,2)+pow(dphi,2));
0631             if(R > 0.4 ) continue;
0632             std::pair<float, float> tow_center, neg_tow_center;
0633             tow_center.first=jet_center.first+deta;
0634             tow_center.second=jet_center.first+dphi;
0635             circle_jet.insert(tow_center);
0636             tow_center.first=jet_center.first+deta;
0637             tow_center.second=jet_center.first-dphi;
0638             circle_jet.insert(tow_center);
0639             tow_center.first=jet_center.first-deta;
0640             tow_center.second=jet_center.first-dphi;
0641             circle_jet.insert(tow_center);
0642             tow_center.first=jet_center.first-deta;
0643             tow_center.second=jet_center.first+dphi;
0644             circle_jet.insert(tow_center);
0645             }
0646     }//this gets the jet centers using emcal binning, using the undorded set allows for double writes to be ok
0647     //std::cout<<"The circluar jet should encompass " <<circle_jet.size() <<" bins in the emcal" <<std::endl;   
0648 //  getE2C(jet_particle_map, topNode); 
0649     std::unordered_set<int> jet_towers_ihcal, jet_towers_ohcal, jet_towers_emcal; 
0650     for(auto pc:/*particle_coords*/ circle_jet)
0651     {
0652         if(abs(pc.first) > 1.2) continue;
0653 //      if ( n_evts < 10 ) std::cout<<"The particle we are looking for is located at eta = " <<pc.first 
0654 //              <<" phi = " << pc.second <<std::endl;
0655         int tne=GetTowerNumber(pc, EMCALMAP.first, EMCALMAP.second);
0656         if(tne != -1 ){ 
0657             jet_towers_emcal.insert(tne);
0658 //          if(n_evts < 10 ) std::cout<<"The identified tower number for the emcal is " <<tne <<std::endl;
0659         }
0660         emtowerenergy[tne]=0;
0661         emtowerenergykt[tne]=0;
0662         int tni=GetTowerNumber(pc, IHCALMAP.first, IHCALMAP.second);
0663         if(tni != -1) jet_towers_ihcal.insert(tni);
0664 //      if(n_evts < 10) std::cout<<"Particle added to the ihcal in bin " <<tni <<std::endl;
0665         int tno=GetTowerNumber(pc, OHCALMAP.first, OHCALMAP.second);
0666         if (tno != -1) jet_towers_ohcal.insert(tno);
0667 //      if (n_evts < 10 ) std::cout<<"Particle added to the ohcal in bin " <<tno <<std::endl;
0668     }
0669 //  if(n_evts < 10) std::cout<<"We have loaded in particles " <<particle_coords.size() <<std::endl;
0670 //  getE2C(topNode, jet_towers_ihcal, jet_towers_ohcal, jet_towers_emcal);
0671     bool kt_tows_good=false;
0672     if(jet_towers_emcal.size() == 0 || jet_towers_ohcal.size() == 0 ) return 1;
0673     try{
0674         std::cout<<" The kt map has size " <<kt_towers.size() <<std::endl;
0675     }
0676     catch(std::exception& e) {std::cout<<"Couldn't access the map to get size" <<std::endl;}
0677     if(kt_towers.size() > 0){
0678     try{
0679         std::cout<<"The tows emcal bit has size " <<kt_towers[0].size() <<std::endl;
0680     }
0681     catch(std::exception& e){std::cout<<"Nothing in the 0th element " <<std::endl;}
0682     try{
0683         std::cout<<"The tows ihcal bit has size " <<kt_towers[1].size() <<std::endl;
0684     }
0685     catch(std::exception& e){std::cout<<"Nothing in the 1st element " <<std::endl;}
0686     try{
0687         std::cout<<"The tows ohcal bit has size " <<kt_towers[2].size() <<std::endl;
0688         kt_tows_good=true;
0689     }
0690     catch(std::exception& e){std::cout<<"Nothing in the 2nd element " <<std::endl;}
0691     }
0692 try{    GetE3C(topNode, jet_particle_map); }//get both in one call
0693     catch(std::exception& e) {std::cout<<"unable to run the particle map, \n Exception " <<e.what() <<std::endl;}
0694     try{GetE3C(topNode, jet_towers_emcal, jet_towers_ihcal, jet_towers_ohcal, &emtowerenergy, false, jet_energy); }//get both values filled in one 
0695     catch(std::exception& e) {std::cout<<"unable to run the tower map, \n Exception " <<e.what() <<std::endl;}
0696     try{if(kt_towers.size() > 0 && kt_tows_good) GetE3C(topNode, kt_towers[0], kt_towers[1], kt_towers[2], &emtowerenergykt, true, jet_energy); }
0697     catch(std::exception& e) {std::cout<<"unable to run the kt map, \n Exception " <<e.what() <<std::endl;}
0698     for(auto m:jettowenergy) comptotows->Fill(m.second, emtowerenergy[m.first]);
0699     return 1;
0700 }
0701 std::map<std::pair<float, float>, std::vector<std::unordered_set<int>>> CalorimeterTowerENC::FindAntiKTTowers(PHCompositeNode* topNode)
0702 {
0703         auto emcal_geom =  findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC");
0704         auto emcal_tower_energy =  findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0705         auto ihcal_geom= findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0706         auto ihcal_tower_energy= findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0707         auto ohcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0708         auto ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0709         //This is a fake anti-kt process, really using E_T instead to allow for jet ID without tracking 
0710         //For this I will just mao all the energy onto a "single" calo, with the Get Tower number to squash 
0711         //emcal down to hcal 
0712         std::map<int, float> emcal_to_hcal_energy;
0713         std::map<int, std::pair<float, float>> hcal_tower_pos;
0714         std::map<std::pair<float, float>, std::vector<std::unordered_set<int>>> tower_jets;     
0715         for(int i=0; i<(int)emcal_tower_energy->size(); i++)
0716         {
0717             auto tower =emcal_tower_energy->get_tower_at_channel(i);
0718             float energy=tower->get_energy();
0719             if(energy==0) continue;
0720             auto key=emcal_tower_energy->encode_key(i);
0721             int phibin = emcal_tower_energy->getTowerPhiBin(key);
0722             int etabin = emcal_tower_energy->getTowerEtaBin(key);
0723             auto phi=emcal_geom->get_phicenter(phibin);
0724             auto eta=emcal_geom->get_etacenter(etabin); 
0725             energy=energy/cosh(eta);
0726             std::pair<float, float> coords {eta, phi};
0727             int hcal_tower_n=-1;
0728             if(emcal_lookup.find(i) == emcal_lookup.end()){
0729                 int t=GetTowerNumber(coords, OHCALMAP.first, OHCALMAP.second);
0730                 emcal_lookup[i]=t;
0731                 hcal_lookup[t].push_back(i);
0732             }
0733             hcal_tower_n=emcal_lookup[i]; 
0734             if(emcal_to_hcal_energy.find(hcal_tower_n) == emcal_to_hcal_energy.end())
0735                 emcal_to_hcal_energy[hcal_tower_n]=energy;
0736             else emcal_to_hcal_energy[hcal_tower_n]+=energy;
0737         }
0738         for(int i=0; i<(int)ihcal_tower_energy->size(); i++)
0739         {
0740             auto tower =ihcal_tower_energy->get_tower_at_channel(i);
0741             float energy=tower->get_energy();
0742             if(energy==0) continue;
0743             auto key=ihcal_tower_energy->encode_key(i);
0744             int etabin = ihcal_tower_energy->getTowerEtaBin(key);
0745             auto eta=ihcal_geom->get_etacenter(etabin); 
0746             energy=energy/cosh(eta);
0747             if(emcal_to_hcal_energy.find(i) == emcal_to_hcal_energy.end())
0748                 emcal_to_hcal_energy[i]=energy;
0749             else emcal_to_hcal_energy[i]+=energy;
0750         }   
0751         for(int i=0; i<(int)ohcal_tower_energy->size(); i++)
0752         {
0753             auto tower =ohcal_tower_energy->get_tower_at_channel(i);
0754             float energy=tower->get_energy();
0755             if(energy==0) continue;
0756             auto key=ohcal_tower_energy->encode_key(i);
0757             int phibin = ohcal_tower_energy->getTowerPhiBin(key);
0758             int etabin = ohcal_tower_energy->getTowerEtaBin(key);
0759             auto phi=ohcal_geom->get_phicenter(phibin);
0760             auto eta=ohcal_geom->get_etacenter(etabin); 
0761             energy=energy/cosh(eta);
0762             std::pair<float, float> coords {eta, phi};
0763             if(emcal_to_hcal_energy.find(i) == emcal_to_hcal_energy.end())
0764                 emcal_to_hcal_energy[i]=energy;
0765             else emcal_to_hcal_energy[i]+=energy;
0766             hcal_tower_pos[i]= coords;
0767 
0768         }   
0769         float etmin=1.0; //10 MeV threshold chosen because
0770         float max_found=0;
0771         int max_bin=0;
0772         std::cout<<"The energy output is non zero in how many towers?: " <<emcal_to_hcal_energy.size() <<std::endl;
0773         for(auto e:emcal_to_hcal_energy){
0774             if(e.second > max_found){
0775                 max_found=e.second;
0776                 max_bin=e.first;
0777             }
0778         }
0779         std::cout<<"The maximum energy is " <<max_found <<std::endl;
0780         while( max_found > etmin) 
0781         {
0782             //now I need to find the jets that look correct 
0783             //auto tower=ohcal_tower_energy->get_tower_at_channel(max_bin);
0784             //auto key=ohcal_tower_energy->encode_key(max_bin);
0785             std::map<int, float> candidates;
0786 //          if(n_evts < 5) std::cout<<"Looking for jet number " <<tower_jets.size() <<std::endl;
0787             std::unordered_set<int> jettows;
0788             auto central=hcal_tower_pos[max_bin];
0789             float jet_pt=emcal_to_hcal_energy[max_bin];
0790             for(auto c:hcal_tower_pos){
0791                 float R=getR(central, c.second);
0792         //      std::cout<<" The seperation is " <<R <<" for a particle at " <<c.second.second <<std::endl;
0793                 if(R < 0.5 && emcal_to_hcal_energy[c.first] > 0.001 ){ candidates[c.first]=R;
0794                 jet_pt+=emcal_to_hcal_energy[c.first];
0795                 }
0796             }
0797             candidates[max_bin]=0.;
0798         //  if(candidates.size() > 2) std::cout<<"The number of candidate towers is " <<candidates.size() <<" with total E_T " <<jet_pt  <<std::endl;
0799             if(jet_pt < 0.005*jet_cutoff){
0800                 for(auto c:candidates){
0801         //          std::cout<<"The Number of bins still in play is now " <<emcal_to_hcal_energy.size() <<std::endl;
0802         //          std::cout<<"The tower needing to leave has bin number=" <<c.first <<" and R=" <<c.second <<std::endl;
0803                     if(emcal_to_hcal_energy.find(c.first) != emcal_to_hcal_energy.end()) emcal_to_hcal_energy.erase(c.first);
0804         //          std::cout<<"Now the energy bins should have decreased by 1? " <<emcal_to_hcal_energy.size() <<std::endl;
0805                 }
0806         //      std::cout<<"Jet Energy was below half a percent of the cutoff" <<std::endl; 
0807                 max_found=0.0;
0808                 for(auto e:emcal_to_hcal_energy){
0809                 if(e.second > max_found){
0810                     max_found=e.second;
0811                     max_bin=e.first;
0812                 }
0813             }
0814                 
0815                 continue;
0816             } 
0817             bool found_all_parts=false;
0818             jet_pt=0.0;
0819             while(!found_all_parts)
0820             {
0821                 //to recombine, need to make an assumption that all the particles are massless
0822                 //I think this is close enough for right now 
0823                 if(candidates.size() < 2){
0824                     for(auto c:candidates){
0825                         if(emcal_to_hcal_energy.find(c.first) !=emcal_to_hcal_energy.end()){
0826                              emcal_to_hcal_energy.erase(c.first);
0827                         }
0828                     }
0829                      break;
0830                 }
0831                 std::map<int, float> di;
0832                 for( auto c:candidates) di[c.first] = pow(emcal_to_hcal_energy[c.first], -2);
0833                 std::map<std::pair<int, int>, float> dij;
0834                 for(auto c:candidates){
0835                     for(auto d:candidates){
0836                         if(c.first >= d.first)continue; //avoid double counting
0837                         std::pair <int, int> couple {c.first, d.first};
0838                         float aktc=pow(emcal_to_hcal_energy[c.first], -2);
0839                         float aktd=pow(emcal_to_hcal_energy[d.first], -2);
0840                         float delta=getR(hcal_tower_pos[c.first], hcal_tower_pos[d.first]);
0841                         delta=pow(delta/0.4, 2);
0842                         dij[couple]=(std::min(aktc, aktd) * delta);
0843                     }
0844                 }
0845                 std::pair<int,int> mergecoords {-1,-1};
0846                 float minmerge=10000000.;
0847                 for(auto d:di){
0848                      if(d.second < minmerge){
0849                         minmerge=d.second; 
0850                         mergecoords.first=d.first;
0851                     }
0852                 }
0853                 for(auto d:dij){
0854                     if(d.second < minmerge){
0855                         minmerge=d.second;
0856                         mergecoords=d.first;
0857                     }
0858                 }
0859                 if(mergecoords.second != -1){
0860                     //this means that we have to merge particles 
0861                     int c = mergecoords.first, d = mergecoords.second;
0862                     int newparticleindex=((--candidates.end())->first) + mergecoords.first*mergecoords.second;
0863                     float E=0., px=0., py=0.,pz=0.;
0864                     E=emcal_to_hcal_energy[c]*cosh(hcal_tower_pos[c].first);
0865                     E+=emcal_to_hcal_energy[d]*cosh(hcal_tower_pos[d].first);
0866                     px=emcal_to_hcal_energy[c]*cos(hcal_tower_pos[c].second);
0867                     px+=emcal_to_hcal_energy[d]*cos(hcal_tower_pos[d].second);
0868                     py=emcal_to_hcal_energy[c]*sin(hcal_tower_pos[c].second);
0869                     py+=emcal_to_hcal_energy[d]*sin(hcal_tower_pos[d].second);
0870                     E=emcal_to_hcal_energy[c]*sinh(hcal_tower_pos[c].first);
0871                     E+=emcal_to_hcal_energy[d]*sinh(hcal_tower_pos[d].first);
0872                     float eta=0., phi=0.;
0873                     eta=atanh(pz/E);
0874                     phi=atan(py/px);
0875                     std::pair<float, float> coords {eta, phi};
0876                     candidates[newparticleindex]=getR(central, coords);
0877                     candidates.erase(c);
0878                     candidates.erase(d);
0879                     //only record the tower indexs
0880                     if(c < (int)ohcal_tower_energy->size()) jettows.insert(c);
0881                     if(d < (int)ohcal_tower_energy->size()) jettows.insert(d);
0882                     emcal_to_hcal_energy[newparticleindex]=sqrt(abs(pow(E,2)-pow(pz, 2)));
0883                     emcal_to_hcal_energy.erase(c);
0884                     emcal_to_hcal_energy.erase(d);
0885                     hcal_tower_pos[newparticleindex]=coords;
0886                     jet_pt=emcal_to_hcal_energy[newparticleindex];
0887 //                  std::cout<<"Merged particles have E_{T} = " <<jet_pt <<std::endl;
0888                 }
0889                 else{
0890                     found_all_parts=true; 
0891                     for(auto c:candidates){
0892                         if(emcal_to_hcal_energy.find(c.first) != emcal_to_hcal_energy.end()) 
0893                             emcal_to_hcal_energy.erase(c.first);
0894                         }
0895                     continue;
0896                 }
0897                 if(candidates.size() <= 1){
0898                     found_all_parts=true;
0899                     continue;
0900                 }
0901                 continue;
0902             }
0903             for(auto e:emcal_to_hcal_energy){
0904                 if(e.second > max_found){
0905                     max_found=e.second;
0906                     max_bin=e.first;
0907                 }
0908             }
0909             std::unordered_set<int> emcaltowers;
0910             for(auto c:jettows)
0911             {
0912                 for(auto b:hcal_lookup[c]) emcaltowers.insert(b);
0913             }
0914             std::vector<std::unordered_set<int>> tows {emcaltowers, jettows, jettows};
0915             if(jet_pt > 0 ) std::cout<<"Identified a Jet with a E_t of " <<jet_pt <<std::endl;
0916             if(jet_pt > 0.01*jet_cutoff){ //trying a much lower cutoff, and use a qucik cluster trick above to decrease the number of trys it takes
0917                 tower_jets[central]=tows;
0918                 calojethits->Fill(central.first, central.second);
0919             }
0920                 
0921             continue;
0922         }
0923 
0924         std::cout<<"Identified jets: " <<tower_jets.size() <<std::endl; 
0925         return tower_jets;
0926     
0927 }
0928 std::map<Jet*, std::pair<float,float>> CalorimeterTowerENC::MatchKtTowers(std::map<std::pair<float, float>, std::vector<std::unordered_set<int>>> jettow, JetContainer* jets){
0929     std::map<Jet*, std::pair<float, float>> coords;
0930     for(auto jet:*jets){
0931         float jeteta=jet->get_eta();
0932         float jetphi=jet->get_phi();
0933         std::pair<float, float> jetcoord {jeteta, jetphi};
0934         float rmin=10.0;
0935         std::pair<float, float> closest {0.0, 0.0};
0936         for(auto t:jettow){
0937             float R=getR(jetcoord, t.first);
0938             if(R < rmin){
0939                 rmin=R;
0940                 closest=t.first;
0941             }
0942         }
0943         coords[jet]=closest;
0944     }
0945     std::cout<<"Matched the jets in the calo to the nearest jet axis" <<std::endl;
0946     return coords;
0947 }
0948 int CalorimeterTowerENC::process_event(PHCompositeNode *topNode){
0949     //need to process the events, this is really going to be a minor thing, need to pull my existing ENC code to do it better
0950     //for the first event build the  tower number map that is needed
0951     n_evts++;
0952     std::cout<<"Running on event " <<n_evts<<std::endl;
0953     try{
0954         Nj=0;
0955         auto emcal_geom =  findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_CEMC");
0956         auto emcal_tower_energy =  findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0957         auto ihcal_geom= findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALIN");
0958         auto ihcal_tower_energy= findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0959         auto ohcal_geom=findNode::getClass<RawTowerGeomContainer_Cylinderv1>(topNode, "TOWERGEOM_HCALOUT");
0960         auto ohcal_tower_energy=findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0961         if(EMCALMAP.first.size() == 0 )  EMCALMAP=GetTowerMaps(emcal_geom, RawTowerDefs::CEMC, emcal_tower_energy);
0962         if(IHCALMAP.first.size() == 0 )  IHCALMAP=GetTowerMaps(ihcal_geom, RawTowerDefs::HCALIN, ihcal_tower_energy);
0963         if(OHCALMAP.first.size() == 0 )  OHCALMAP=GetTowerMaps(ohcal_geom, RawTowerDefs::HCALOUT, ohcal_tower_energy);
0964         auto jets = findNode::getClass<JetContainer> (topNode, "AntiKt_Truth_r04");
0965         float emcal_energy=0, ihcal_energy=0, ohcal_energy=0;   
0966         try{
0967         for(int i=0; i<(int)emcal_tower_energy->size(); i++ ){
0968             auto tower = emcal_tower_energy->get_tower_at_channel(i);
0969             float energy=tower->get_energy();
0970             emcal_energy+=energy;
0971         }
0972         for(int i=0; i<(int)ihcal_tower_energy->size(); i++ ){
0973             auto tower = ihcal_tower_energy->get_tower_at_channel(i);
0974             float energy=tower->get_energy();
0975             ihcal_energy+=energy;
0976         }
0977         for(int i=0; i<(int)ohcal_tower_energy->size(); i++ ){
0978             auto tower = ohcal_tower_energy->get_tower_at_channel(i);
0979             float energy=tower->get_energy();
0980             ohcal_energy+=energy;
0981         }
0982         EM_energy->Fill(emcal_energy);
0983         OH_energy->Fill(ohcal_energy);
0984         IH_energy->Fill(ihcal_energy);
0985         try{Nj=jets->size();
0986         auto towerktjets=FindAntiKTTowers(topNode);
0987         try{auto coordinatedjets=MatchKtTowers(towerktjets, jets);
0988         std::vector<std::thread> jetThreads;
0989         try{for(auto j:*jets){
0990             if(j->get_pt() < jet_cutoff){ //put in a 10 GeV cut on the jets
0991                 Nj--;
0992                 continue;
0993             }
0994             this->histogram_map["parts"]->pt->Fill(j->get_pt());
0995             /*jetThreads.push_back(std::thread(&*/CalorimeterTowerENC::RecordHits(/*, this,*/ topNode, j, towerktjets[coordinatedjets[j]]);
0996         }
0997     //  for(int t=0; t<(int)jetThreads.size(); t++) jetThreads.at(t).join();
0998         number_of_jets->Fill(Nj);
0999         }
1000         catch(std::exception& e4) {std::cout<<"The threads objects are the problems \n Exception: "<<e4.what() <<std::endl;}
1001         }
1002         catch(std::exception& e3) {std::cout<<"The matches are the problems \n Exception: "<<e3.what() <<std::endl;}
1003         }
1004         catch(std::exception& e2) {std::cout<<"The anti-kt objects are the problems \n Exception: "<<e2.what() <<std::endl;}
1005         }
1006         catch(std::exception& e1) {std::cout<<"The energy objects are the problems \n Exception: "<<e1.what() <<std::endl;}
1007     }
1008     catch(std::exception& e ) {
1009         std::cout<<"Failed to run on event " <<n_evts <<"\n Skipping Event" <<" \n Exception was " <<e.what() <<std::endl;  
1010         n_evts--;
1011     }
1012     return Fun4AllReturnCodes::EVENT_OK;
1013     
1014     
1015 }
1016 int CalorimeterTowerENC::End(PHCompositeNode *topNode){
1017     return 1;
1018 }
1019 void CalorimeterTowerENC::Print(const std::string &what) const
1020 {
1021     std::cout<<"Printing the files out now to a file named: " <<outfilename <<std::endl;
1022     TFile* f=new TFile(outfilename.c_str(), "RECREATE");
1023     //printing the stuff out 
1024     for(auto hs:histogram_map){
1025         //hs.second->R->Scale(1/(float)hs.second->N->GetEntries()); //average over number of jets
1026         hs.second->E2C->Scale(1/(float)hs.second->E2C->GetBinWidth(5)); //correct for binning
1027         hs.second->E3C->Scale(1/(float)hs.second->E3C->GetBinWidth(5)); //correct for binning
1028         //hs.second->E2C->Scale(1/(float)n_evts); //average over events
1029         //hs.second->E3C->Scale(1/(float)n_evts); //average over events
1030         for(auto h:hs.second->histsvector){
1031             h->Write();
1032         }
1033         hs.second->R_pt->Write();
1034     }
1035     jethits->Write();
1036     comptotows->Write();
1037     number_of_jets->Write();
1038     EM_energy->Write();
1039     IH_energy->Write();
1040     OH_energy->Write();
1041     f->Write();
1042     f->Close();
1043     return;
1044 }
[ 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!