Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
sPHAnalysis_calo.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file sPHAnalysis_calo.cc
1 #include "sPHAnalysis_calo.h"
2 
3 #include <cmath>
4 #include <cstdlib>
5 #include <cstdio>
6 #include <iostream>
7 #include <iomanip>
8 #include <fstream>
9 
10 #include "TFile.h"
11 #include "TNtuple.h"
12 #include "TH1F.h"
13 #include "TH1D.h"
14 #include "TF1.h"
15 #include "TLorentzVector.h"
16 #include "TRandom2.h"
17 
18 #include <globalvertex/GlobalVertexMap.h>
19 #include <globalvertex/GlobalVertex.h>
20 
21 #include <calobase/RawClusterContainer.h>
22 #include <calobase/RawCluster.h>
23 #include <calobase/RawClusterv1.h>
24 #include <calobase/RawTowerv2.h>
25 #include <calobase/RawTowerGeom.h>
26 #include <calobase/RawTowerContainer.h>
27 #include <calobase/TowerInfov1.h>
28 #include <calobase/TowerInfoContainerv1.h>
29 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
30 #include <calobase/RawTowerGeomContainer.h>
31 
32 #include <phool/getClass.h>
33 #include <phool/recoConsts.h>
34 #include <phool/PHCompositeNode.h>
35 #include <phool/PHIODataNode.h>
36 #include <phool/PHRandomSeed.h>
37 #include <fun4all/Fun4AllReturnCodes.h>
38 
39 #include <gsl/gsl_rng.h>
40 
41 using namespace std;
42 
43 //==============================================================
44 
45 sPHAnalysis_calo::sPHAnalysis_calo(const std::string &name, const std::string &filename) : SubsysReco(name)
46 {
47  OutputNtupleFile=nullptr;
48  OutputFileName=filename;
49  EventNumber=0;
50  ntp_notracking=nullptr;
51  h_mult=nullptr;
52  h_ecore= nullptr;
53  _whattodo = 0;
54  _rng = nullptr;
55 }
56 
57 //==============================================================
58 
59 int sPHAnalysis_calo::Init(PHCompositeNode *topNode)
60 {
61  std::cout << "sPHAnalysis_calo::Init started..." << endl;
62  OutputNtupleFile = new TFile(OutputFileName.c_str(),"RECREATE");
63  std::cout << "sPHAnalysis_calo::Init: output file " << OutputFileName.c_str() << " opened." << endl;
64 
65  ntp_notracking = new TNtuple("ntp_notracking","","mass:pt:e1:e2:hoe1:hoe2");
66  h_mult = new TH1D("h_mult","",1000,0.,10000.);
67  h_ecore = new TH1D("h_ecore","",200,0.,20.);
68 
69  char hname[99];
70  for(int i=0; i<256*96; i++) {
71  sprintf(hname,"h_pedestal_%d",i);
72  h_pedestal[i] = new TH1F(hname,hname,10000,-0.5,9999.5);
73  }
74 
75  _rng = new TRandom2();
76  _rng->SetSeed(0);
77 
78  cout << "sPHAnalysis_calo::Init ended." << endl;
79  return Fun4AllReturnCodes::EVENT_OK;
80 }
81 
82 //==============================================================
83 
84 int sPHAnalysis_calo::InitRun(PHCompositeNode *topNode)
85 {
86  cout << "sPHAnalysis_calo::InitRun started..." << endl;
87  cout << "sPHAnalysis_calo::InitRun ended." << endl;
88  return Fun4AllReturnCodes::EVENT_OK;
89 }
90 
91 //==============================================================
92 
93 int sPHAnalysis_calo::process_event(PHCompositeNode *topNode)
94 {
95  if(_whattodo==0) {
96  return process_event_test(topNode); // simulation
97  } else if(_whattodo==1) {
98  return process_event_data(topNode); // data
99  } else {
100  cerr << "ERROR: wrong choice of what to do." << endl; return Fun4AllReturnCodes::ABORTRUN;
101  }
102 }
103 
104 //======================================================================
105 
106 int sPHAnalysis_calo::process_event_data(PHCompositeNode *topNode) {
107  EventNumber++;
108  float tmp1[99];
109 
110  cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
111  if(EventNumber==1) topNode->print();
112 
113  TowerInfoContainer* _towersCEMC = findNode::getClass<TowerInfoContainerv1>(topNode, "TOWERS_CEMC");
114  if (!_towersCEMC) { std::cerr<<"No TOWERS_CEMC Node"<<std::endl; return Fun4AllReturnCodes::ABORTEVENT; }
115  unsigned int ntowers = _towersCEMC->size();
116  cout << " # of towers = " << ntowers << endl;
117 
118  for (unsigned int channel = 0; channel < ntowers; channel++)
119  {
120  TowerInfo *cemcinfo_raw = _towersCEMC->get_tower_at_channel(channel);
121  float raw_amplitude = cemcinfo_raw->get_energy();
122  float raw_time = cemcinfo_raw->get_time();
123  unsigned int key = _towersCEMC->encode_key(channel);
124  unsigned int iphi = _towersCEMC->getTowerPhiBin(key);
125  unsigned int ieta = _towersCEMC->getTowerEtaBin(key);
126  h_pedestal[channel]->Fill(raw_amplitude);
127  if(channel>0 && channel<20) cout << channel << " " << iphi << " " << ieta << " " << raw_amplitude << " " << raw_time << endl;
128  }
129 
130  cout << "sPHAnalysis_calo::procedss_event_data() ended." << endl;
131  return Fun4AllReturnCodes::EVENT_OK;
132 }
133 
134 //======================================================================
135 
136 int sPHAnalysis_calo::process_event_test(PHCompositeNode *topNode) {
137  EventNumber++;
138  float tmp1[99];
139 
140  cout<<"--------------------------- EventNumber = " << EventNumber-1 << endl;
141  if(EventNumber==1) topNode->print();
142 
143 // Find event vertex
144  GlobalVertexMap *global_vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
145  //cout << "Number of GlobalVertexMap entries = " << global_vtxmap->size() << endl;
146  double Zvtx = 0.;
147  for (GlobalVertexMap::Iter iter = global_vtxmap->begin(); iter != global_vtxmap->end(); ++iter)
148  {
149  GlobalVertex *vtx = iter->second;
150  if(vtx->get_id()==1) Zvtx = vtx->get_z(); // BBC vertex (?)
151  //cout << "Global vertex: " << vtx->get_id() << " " << vtx->get_z() << " " << vtx->get_t() << endl;
152  }
153  cout << "Global vertex Z = " << Zvtx << endl;
154 
155 // Find node with CEMC clusters
156  RawClusterContainer* cemc_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
157  if(!cemc_clusters) {
158  cerr << PHWHERE << " ERROR: Can not find CLUSTER_CEMC node." << endl;
159  return Fun4AllReturnCodes::ABORTEVENT;
160  }
161  cout << " Number of CEMC clusters = " << cemc_clusters->size() << endl;
162  h_mult->Fill(cemc_clusters->size());
163 
164 // Fetch CEMC geometry
165  RawTowerGeomContainer* _geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
166  if(!_geomEM) std::cerr<<"No TOWERGEOM_CEMC"<<std::endl;
167  RawTowerGeomContainer* _geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
168  if(!_geomIH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
169  RawTowerGeomContainer* _geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
170  if(!_geomIH) std::cerr<<"No TOWERGEOM_HCALIN"<<std::endl;
171 
172  RawTowerContainer* _towersRawEM = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
173  if (!_towersRawEM) std::cerr<<"No TOWER_CALIB_CEMC Node"<<std::endl;
174  RawTowerContainer* _towersRawIH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
175  if (!_towersRawIH) std::cerr<<"No TOWER_CALIB_HCALIN Node"<<std::endl;
176  RawTowerContainer* _towersRawOH = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
177  if (!_towersRawOH) std::cerr<<"No TOWER_CALIB_HCALOUT Node"<<std::endl;
178 
179  vector<TLorentzVector> electrons;
180  vector<double> vhoe;
181 
182  // loop over CEMC clusters with E > 2 GeV
183  RawClusterContainer::Range begin_end = cemc_clusters->getClusters();
184  RawClusterContainer::Iterator iter;
185  for (iter = begin_end.first; iter != begin_end.second; ++iter)
186  {
187  RawCluster* cluster = iter->second;
188  if(!cluster) { cout << "ERROR: bad cluster pointer = " << cluster << endl; continue; }
189  else {
190  double cemc_ecore = cluster->get_ecore();
191  h_ecore->Fill(cemc_ecore);
192  if(cemc_ecore<2.0) continue;
193  double cemc_x = cluster->get_x();
194  double cemc_y = cluster->get_y();
195  double cemc_z = cluster->get_z() - Zvtx; // correct for event vertex position
196  //double cemc_r = cluster->get_r();
197  double cemc_r = sqrt(pow(cemc_x,2)+pow(cemc_y,2));
198  double cemc_rr = sqrt(pow(cemc_r,2)+pow(cemc_z,2));
199  double cemc_eta = asinh( cemc_z / cemc_r );
200  double cemc_phi = atan2( cemc_y, cemc_x );
201  double cemc_px = cemc_ecore * (cemc_x/cemc_rr);
202  double cemc_py = cemc_ecore * (cemc_y/cemc_rr);
203  double cemc_pz = cemc_ecore * (cemc_z/cemc_rr);
204  //cout << "CEMC cluster: " << cemc_ecore << " " << cemc_eta << " " << cemc_phi << endl;
205 
206  // find closest CEMC tower in eta-phi space
207  double distem = 99999.;
208  RawTower* thetowerem = nullptr;
209  RawTowerContainer::ConstRange begin_end_rawEM = _towersRawEM->getTowers();
210  for (RawTowerContainer::ConstIterator rtiter = begin_end_rawEM.first; rtiter != begin_end_rawEM.second; ++rtiter) {
211  RawTower *tower = rtiter->second;
212  RawTowerGeom *tower_geom = _geomEM->get_tower_geometry(tower->get_key());
213  double cemc_tower_phi = tower_geom->get_phi();
214  double cemc_tower_x = tower_geom->get_center_x();
215  double cemc_tower_y = tower_geom->get_center_y();
216  double cemc_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
217  double cemc_tower_r = sqrt(pow(cemc_tower_x,2)+pow(cemc_tower_y,2));
218  double cemc_tower_eta = asinh( cemc_tower_z / cemc_tower_r );
219  double tmpdist = sqrt(pow(cemc_eta-cemc_tower_eta,2)+pow(cemc_phi-cemc_tower_phi,2));
220  if(tmpdist<distem) { distem = tmpdist; thetowerem = tower; }
221  }
222  RawTowerGeom *thetower_geom_em = _geomEM->get_tower_geometry(thetowerem->get_key());
223  unsigned int ietaem = thetower_geom_em->get_bineta();
224  unsigned int jphiem = thetower_geom_em->get_binphi();
225 
226  // find closest HCALIN tower in eta-phi space
227  double distin = 99999.;
228  RawTower* thetowerin = nullptr;
229  RawTowerContainer::ConstRange begin_end_rawIN = _towersRawIH->getTowers();
230  for (RawTowerContainer::ConstIterator rtiter = begin_end_rawIN.first; rtiter != begin_end_rawIN.second; ++rtiter) {
231  RawTower *tower = rtiter->second;
232  RawTowerGeom *tower_geom = _geomIH->get_tower_geometry(tower->get_key());
233  double hcalin_tower_phi = tower_geom->get_phi();
234  double hcalin_tower_x = tower_geom->get_center_x();
235  double hcalin_tower_y = tower_geom->get_center_y();
236  double hcalin_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
237  double hcalin_tower_r = sqrt(pow(hcalin_tower_x,2)+pow(hcalin_tower_y,2));
238  double hcalin_tower_eta = asinh( hcalin_tower_z / hcalin_tower_r );
239  double tmpdist = sqrt(pow(cemc_eta-hcalin_tower_eta,2)+pow(cemc_phi-hcalin_tower_phi,2));
240  if(tmpdist<distin) { distin = tmpdist; thetowerin = tower; }
241  }
242  RawTowerGeom *thetower_geom = _geomIH->get_tower_geometry(thetowerin->get_key());
243  unsigned int ieta = thetower_geom->get_bineta();
244  unsigned int jphi = thetower_geom->get_binphi();
245 
246  // Calcuate 3x3 energy deposit in HCALIN. Inner HCAL has 24 bins in eta and 64 bins in phi
247  double e3x3in = 0.;
248  if(ieta<1 || ieta>22) continue; // ignore clusters in edge towers
249  for(unsigned int i=0; i<=2; i++) {
250  for(unsigned int j=0; j<=2; j++) {
251  unsigned int itmp = ieta-1+i;
252  unsigned int jtmp = 0;
253  if(jphi==0 && j==0) { jtmp = 63; } // wrap around
254  else if(jphi==63 && j==2) { jtmp = 0; } // wrap around
255  else { jtmp = jphi-1+j; }
256  RawTower* tmptower = _towersRawIH->getTower(itmp,jtmp);
257  if(tmptower) { e3x3in += tmptower->get_energy(); }
258  }
259  }
260  //h_notracking_eoh->Fill(e3x3in/cemc_ecore);
261 
262  // find closest HCALOUT tower in eta-phi space
263  double distout = 99999.;
264  RawTower* thetowerout = nullptr;
265  RawTowerContainer::ConstRange begin_end_rawON = _towersRawOH->getTowers();
266  for (RawTowerContainer::ConstIterator rtiter = begin_end_rawON.first; rtiter != begin_end_rawON.second; ++rtiter) {
267  RawTower *tower = rtiter->second;
268  RawTowerGeom *tower_geom = _geomOH->get_tower_geometry(tower->get_key());
269  double hcalout_tower_phi = tower_geom->get_phi();
270  double hcalout_tower_x = tower_geom->get_center_x();
271  double hcalout_tower_y = tower_geom->get_center_y();
272  double hcalout_tower_z = tower_geom->get_center_z() - Zvtx; // correct for event vertex
273  double hcalout_tower_r = sqrt(pow(hcalout_tower_x,2)+pow(hcalout_tower_y,2));
274  double hcalout_tower_eta = asinh( hcalout_tower_z / hcalout_tower_r );
275  double tmpdist = sqrt(pow(cemc_eta-hcalout_tower_eta,2)+pow(cemc_phi-hcalout_tower_phi,2));
276  if(tmpdist<distout) { distout = tmpdist; thetowerout = tower; }
277  }
278  RawTowerGeom *thetower_geomout = _geomOH->get_tower_geometry(thetowerout->get_key());
279  unsigned int ietaout = thetower_geomout->get_bineta();
280  unsigned int jphiout = thetower_geomout->get_binphi();
281 
282  // Calcuate 3x3 energy deposit in HCALOUT. Outer HCAL has 24 bins in eta and 64 bins in phi
283  double e3x3out = 0.;
284  if(ietaout<1 || ietaout>22) continue; // ignore clusters in edge towers
285  for(unsigned int i=0; i<=2; i++) {
286  for(unsigned int j=0; j<=2; j++) {
287  unsigned int itmp = ietaout-1+i;
288  unsigned int jtmp = 0;
289  if(jphiout==0 && j==0) { jtmp = 63; } // wrap around
290  else if(jphiout==63 && j==2) { jtmp = 0; } // wrap around
291  else { jtmp = jphiout-1+j; }
292  RawTower* tmptower = _towersRawOH->getTower(itmp,jtmp);
293  if(tmptower) { e3x3out += tmptower->get_energy(); }
294  }
295  }
296 
297  if(e3x3in/cemc_ecore>0.06) continue; // reject hadrons, 90% eID efficiency is with 0.058 cut
298 
299  TLorentzVector tmp = TLorentzVector(cemc_px,cemc_py,cemc_pz,cemc_ecore);
300  electrons.push_back(tmp);
301  vhoe.push_back(e3x3in/cemc_ecore);
302 
303  } // valid CEMC cluster
304  } // end loop over CEMC clusters
305  cout << "number of selected electrons = " << electrons.size() << " " << vhoe.size() << endl;
306 
307 // Make Upsilons
308 
309  if(electrons.size()>=2) {
310  for(long unsigned int i=0; i<electrons.size()-1; i++) {
311  for(long unsigned int j=i+1; j<electrons.size(); j++) {
312  TLorentzVector pair = electrons[i]+electrons[j];
313  cout << i << " " << j << endl;
314  tmp1[0] = pair.M();
315  tmp1[1] = pair.Pt();
316  cout << pair.M() << " " << pair.Pt() << endl;
317  tmp1[2] = (electrons[i]).Pt();
318  tmp1[3] = (electrons[j]).Pt();
319  cout << vhoe[i] << " " << vhoe[j] << endl;
320  tmp1[4] = vhoe[i];
321  tmp1[5] = vhoe[j];
322  cout << "filling..." << endl;
323  ntp_notracking->Fill(tmp1);
324  cout << "done." << endl;
325  }
326  }
327  }
328 
329  cout << "sPHAnalysis_calo::procedss_event_test() ended." << endl;
330  return Fun4AllReturnCodes::EVENT_OK;
331 }
332 
333 //======================================================================
334 
335 int sPHAnalysis_calo::End(PHCompositeNode *topNode)
336 {
337  std::cout << "sPHAnalysis_calo::End() started, Number of processed events = " << EventNumber << endl;
338  OutputNtupleFile->Write();
339  OutputNtupleFile->Close();
340  return Fun4AllReturnCodes::EVENT_OK;
341 }
342 
343