ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CopyAndSubtractJets.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CopyAndSubtractJets.cc
1 #include "CopyAndSubtractJets.h"
2 
3 #include "TowerBackground.h"
4 
5 #include <calobase/RawTower.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <calobase/RawTowerGeom.h>
8 #include <calobase/RawTowerGeomContainer.h>
9 
10 #include <g4jets/Jet.h>
11 #include <g4jets/JetMap.h>
12 #include <g4jets/JetMapv1.h>
13 #include <g4jets/Jetv1.h>
14 
16 #include <fun4all/SubsysReco.h>
17 
18 #include <phool/PHCompositeNode.h>
19 #include <phool/PHIODataNode.h>
20 #include <phool/PHNode.h>
21 #include <phool/PHNodeIterator.h>
22 #include <phool/PHObject.h>
23 #include <phool/getClass.h>
24 #include <phool/phool.h>
25 
26 // standard includes
27 #include <cmath>
28 #include <iostream>
29 #include <map>
30 #include <vector>
31 #include <utility>
32 
34  : SubsysReco(name)
35  , _use_flow_modulation(false)
36 {
37 }
38 
40 {
41  CreateNode(topNode);
42 
44 }
45 
47 {
48  if (Verbosity() > 0)
49  std::cout << "CopyAndSubtractJets::process_event: entering, with _use_flow_modulation = " << _use_flow_modulation << std::endl;
50 
51  // pull out needed calo tower info
52  RawTowerContainer *towersEM3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER");
53  RawTowerContainer *towersIH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
54  RawTowerContainer *towersOH3 = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
55 
56  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
57  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
58 
59  // pull out jets and background
60  JetMap *unsub_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_HIRecoSeedsRaw_r02");
61  JetMap *sub_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_HIRecoSeedsSub_r02");
62 
63  TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerBackground_Sub1");
64 
65  std::vector<float> background_UE_0 = background->get_UE(0);
66  std::vector<float> background_UE_1 = background->get_UE(1);
67  std::vector<float> background_UE_2 = background->get_UE(2);
68 
69  float background_v2 = background->get_v2();
70  float background_Psi2 = background->get_Psi2();
71 
72  if (Verbosity() > 0)
73  {
74  std::cout << "CopyAndSubtractJets::process_event: entering with # unsubtracted jets = " << unsub_jets->size() << std::endl;
75  std::cout << "CopyAndSubtractJets::process_event: entering with # subtracted jets = " << sub_jets->size() << std::endl;
76  }
77 
78  // iterate over old jets
79  int ijet = 0;
80  for (JetMap::Iter iter = unsub_jets->begin(); iter != unsub_jets->end(); ++iter)
81  {
82  Jet *this_jet = iter->second;
83 
84  float this_pt = this_jet->get_pt();
85  float this_phi = this_jet->get_phi();
86  float this_eta = this_jet->get_eta();
87 
88  Jet *new_jet = new Jetv1();
89 
90  float new_total_px = 0;
91  float new_total_py = 0;
92  float new_total_pz = 0;
93  float new_total_e = 0;
94 
95  //if (this_jet->get_pt() < 5) continue;
96 
97  if (Verbosity() > 1 && this_jet->get_pt() > 5)
98  std::cout << "CopyAndSubtractJets::process_event: unsubtracted jet with pt / eta / phi = " << this_pt << " / " << this_eta << " / " << this_phi << std::endl;
99 
100  for (Jet::ConstIter comp = this_jet->begin_comp(); comp != this_jet->end_comp(); ++comp)
101  {
102  RawTower *tower = 0;
103  RawTowerGeom *tower_geom = 0;
104 
105  double comp_e = 0;
106  double comp_eta = 0;
107  double comp_phi = 0;
108 
109  int comp_ieta = 0;
110 
111  double comp_background = 0;
112 
113  if ((*comp).first == 5)
114  {
115  tower = towersIH3->getTower((*comp).second);
116  tower_geom = geomIH->get_tower_geometry(tower->get_key());
117 
118  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
119  comp_background = background_UE_1.at(comp_ieta);
120  }
121  else if ((*comp).first == 7)
122  {
123  tower = towersOH3->getTower((*comp).second);
124  tower_geom = geomOH->get_tower_geometry(tower->get_key());
125 
126  comp_ieta = geomOH->get_etabin(tower_geom->get_eta());
127  comp_background = background_UE_2.at(comp_ieta);
128  }
129  else if ((*comp).first == 13)
130  {
131  tower = towersEM3->getTower((*comp).second);
132  tower_geom = geomIH->get_tower_geometry(tower->get_key());
133 
134  comp_ieta = geomIH->get_etabin(tower_geom->get_eta());
135  comp_background = background_UE_0.at(comp_ieta);
136  }
137 
138  if (tower)
139  comp_e = tower->get_energy();
140  if (tower_geom)
141  {
142  comp_eta = tower_geom->get_eta();
143  comp_phi = tower_geom->get_phi();
144  }
145 
146  if (Verbosity() > 4 && this_jet->get_pt() > 5)
147  {
148  std::cout << "CopyAndSubtractJets::process_event: --> constituent in layer " << (*comp).first << ", has unsub E = " << comp_e << ", is at ieta #" << comp_ieta << ", and has UE = " << comp_background << std::endl;
149  }
150 
151  // flow modulate background if turned on
153  {
154  comp_background = comp_background * (1 + 2 * background_v2 * cos(2 * (comp_phi - background_Psi2)));
155  if (Verbosity() > 4 && this_jet->get_pt() > 5)
156  std::cout << "CopyAndSubtractJets::process_event: --> --> flow mod, at phi = " << comp_phi << ", v2 and Psi2 are = " << background_v2 << " , " << background_Psi2 << ", UE after modulation = " << comp_background << std::endl;
157  }
158 
159  // update constituent energy based on the background
160  double comp_sub_e = comp_e - comp_background;
161 
162  // now define new kinematics
163 
164  double comp_px = comp_sub_e / cosh(comp_eta) * cos(comp_phi);
165  double comp_py = comp_sub_e / cosh(comp_eta) * sin(comp_phi);
166  double comp_pz = comp_sub_e * tanh(comp_eta);
167 
168  new_total_px += comp_px;
169  new_total_py += comp_py;
170  new_total_pz += comp_pz;
171  new_total_e += comp_sub_e;
172  }
173 
174  new_jet->set_px(new_total_px);
175  new_jet->set_py(new_total_py);
176  new_jet->set_pz(new_total_pz);
177  new_jet->set_e(new_total_e);
178  new_jet->set_id(ijet);
179 
180  sub_jets->insert(new_jet);
181 
182  if (Verbosity() > 1 && this_pt > 5)
183  {
184  std::cout << "CopyAndSubtractJets::process_event: old jet #" << ijet << ", old px / py / pz / e = " << this_jet->get_px() << " / " << this_jet->get_py() << " / " << this_jet->get_pz() << " / " << this_jet->get_e() << std::endl;
185  std::cout << "CopyAndSubtractJets::process_event: new jet #" << ijet << ", new px / py / pz / e = " << new_jet->get_px() << " / " << new_jet->get_py() << " / " << new_jet->get_pz() << " / " << new_jet->get_e() << std::endl;
186  }
187 
188  ijet++;
189  }
190 
191  if (Verbosity() > 0)
192  {
193  std::cout << "CopyAndSubtractJets::process_event: exiting with # subtracted jets = " << sub_jets->size() << std::endl;
194  }
195 
197 }
198 
200 {
202 }
203 
205 {
206  PHNodeIterator iter(topNode);
207 
208  // Looking for the DST node
209  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
210  if (!dstNode)
211  {
212  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
214  }
215 
216  // Looking for the ANTIKT node
217  PHCompositeNode *antiktNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "ANTIKT"));
218  if (!antiktNode)
219  {
220  std::cout << PHWHERE << "ANTIKT node not found, doing nothing." << std::endl;
221  }
222 
223  // Looking for the TOWER node
224  PHCompositeNode *towerNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "TOWER"));
225  if (!towerNode)
226  {
227  std::cout << PHWHERE << "TOWER node not found, doing nothing." << std::endl;
228  }
229 
230  // store the new jet collection
231  JetMap *test_jets = findNode::getClass<JetMap>(topNode, "AntiKt_Tower_HIRecoSeedsSub_r02");
232  if (!test_jets)
233  {
234  if (Verbosity() > 0) std::cout << "CopyAndSubtractJets::CreateNode : creating AntiKt_Tower_HIRecoSeedsSub_r02 node " << std::endl;
235 
236  JetMap *sub_jets = new JetMapv1();
237  PHIODataNode<PHObject> *subjetNode = new PHIODataNode<PHObject>(sub_jets, "AntiKt_Tower_HIRecoSeedsSub_r02", "PHObject");
238  towerNode->addNode(subjetNode);
239  }
240  else
241  {
242  std::cout << "CopyAndSubtractJets::CreateNode : AntiKt_Tower_HIRecoSeedsSub_r02 already exists! " << std::endl;
243  }
244 
246 }