ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CaloCalibEmc_Pi0.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CaloCalibEmc_Pi0.cc
1 #include "CaloCalibEmc_Pi0.h"
2 
3 #include <calobase/RawCluster.h>
4 #include <calobase/RawClusterContainer.h>
5 #include <calobase/RawClusterUtility.h>
6 #include <calobase/RawTower.h>
7 #include <calobase/RawTowerContainer.h>
8 #include <calobase/RawTowerGeom.h>
9 #include <calobase/RawTowerGeomContainer.h>
10 
12 #include <fun4all/SubsysReco.h>
13 
14 #include <phool/PHCompositeNode.h>
15 #include <phool/getClass.h>
16 #include <phool/phool.h>
17 
18 #include <TFile.h>
19 #include <TNtuple.h>
20 #include <TStyle.h>
21 #include <TF1.h>
22 #include <TGraph.h>
23 
25 
27 #include <g4vertex/GlobalVertex.h>
28 
29 #include <TLorentzVector.h>
30 
31 #include <TMath.h>
32 #include <TTree.h>
33 #include <TFile.h>
34 #include <TH1.h>
35 #include <TH2.h>
36 #include <TH3.h>
37 #include <TStyle.h>
38 #include <TSystem.h>
39 
40 #include <cstdlib>
41 #include <iostream>
42 #include <memory>
43 
44 using namespace std;
45 
46 //____________________________________________________________________________..
47 CaloCalibEmc_Pi0::CaloCalibEmc_Pi0(const std::string& name, const std::string& filename):
48  SubsysReco(name)
49  , _ievent(0)
50  , cal_output(0)
51  , _caloname("CEMC")
52  , _filename(filename)
53  , mass_eta(0)
54  , mass_eta_phi(0)
55  , pairInvMassTotal(0)
56  ,_eventTree(0)
57  ,_eventNumber(-1)
58  ,_nClusters(-1)
59  ,maxTowerEta(-1)
60  ,maxTowerPhi(-1)
61  ,alphaCut(-1.0)
62 {
63 
64 
65  for (int nj=0; nj < 96; nj++)
66  {
67  eta_hist[nj] = NULL;
68  for (int nk=0; nk < 256; nk++)
69  {
70  cemc_hist_eta_phi[nj][nk] = NULL;
71  }
72  }
73 
74 }
75 
76 
77 //____________________________________________________________________________..
79 {
80  cout << "LiteCaloEval::Init(PHCompositeNode *topNode) Initializing" << endl;
81 
82  _ievent = 0;
83 
84 
85  cal_output = new TFile(_filename.c_str(),"RECREATE");
86 
87 
88  pairInvMassTotal = new TH1F("pairInvMassTotal", "Pair Mass Histo", 70, 0.0, 0.7);
89  mass_eta = new TH2F("mass_eta", "2d Pair Mass Histo", 70, 0.0, 0.7, 400,-1.5,1.5);
90  mass_eta_phi = new TH3F("mass_eta_phi", "3d Pair Mass Histo", 70, 0.0, 0.7, 150,-1.5,1.5, 256, -3.142, 3.142);
91 
92 
93 
94  // histo to record every tower by tower locations
95  for (int i = 0; i<96; i++) // eta rows
96  {
97  for (int j = 0; j<258; j++) // phi columns
98  {
99  TString i1;
100  TString j1;
101  i1.Form("%d",i);
102  j1.Form("%d",j);
103  TString hist_name = "emc_ieta" + i1 + "_phi"+ j1;
104 
105  cemc_hist_eta_phi[i][j] = new TH1F(hist_name.Data(),"Hist_ieta_phi_",70,0.0,0.7);
106 
107  }
108  }
109 
110  // histo to record each eta locations (with all phis included in each)
111  for (int i = 0; i < 96; i++)
112  {
113  gStyle->SetOptFit(1);
114  TString a;
115  a.Form("%d",i);
116  TString b = "eta_" + a;
117 
118  eta_hist[i] = new TH1F (b.Data(),"eta and all phi's",70,0.0,0.7);
119 
120  }
121 
122  if (topNode != 0)
123  {
124  // TTree declare
125  _eventTree = new TTree("_eventTree", "An event level info Tree");
126  // TTree branches
127  _eventTree->Branch("_eventNumber", &_eventNumber, "_eventNumber/I");
128  _eventTree->Branch("_nClusters", &_nClusters, "_nClusters/I");
129  _eventTree->Branch("_clusterIDs", _clusterIDs, "_clusterIDs[_nClusters]/I");
130  _eventTree->Branch("_clusterEnergies", _clusterEnergies, "_clusterEnergies[_nClusters]/F");
131  _eventTree->Branch("_clusterPts", _clusterPts, "_clusterPts[_nClusters]/F");
132  _eventTree->Branch("_clusterEtas", _clusterEtas, "_clusterEtas[_nClusters]/I");
133  _eventTree->Branch("_clusterPhis", _clusterPhis, "_clusterPhis[_nClusters]/I");
134  _eventTree->Branch("_maxTowerEtas", _maxTowerEtas, "_maxTowerEtas[_nClusters]/I");
135  _eventTree->Branch("_maxTowerPhis", _maxTowerPhis, "_maxTowerPhis[_nClusters]/I");
136  }
137 
139 }
140 
141 //____________________________________________________________________________..
143 {
144  if (_ievent % 50 == 0)
145  {
146  cout << endl;
147  cout << "Beginning of the event " << _ievent << endl;
148  cout << "====================================" << endl;
149  }
150 
151  _eventNumber = _ievent++;
152 
153  // create a cluster object
154  RawClusterContainer *recal_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_POS_COR_CEMC");
155 
156  // create a tower object
157  string towernode = "TOWER_CALIB_" + _caloname;
158  RawTowerContainer* _towers = findNode::getClass<RawTowerContainer>(topNode, towernode.c_str());
159  if (!_towers)
160  {
161  cerr << PHWHERE << " ERROR: Can't find " << towernode << endl;
162  }
163 
164  // create a tower geometry object
165  string towergeomnode = "TOWERGEOM_" + _caloname;
166  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnode.c_str());
167  if (!towergeom)
168  {
169  cout << PHWHERE << ": Could not find node " << towergeomnode << endl;
171  }
172 
173 
174 
175  // Get Vertex
176  float vx = 0;
177  float vy = 0;
178  float vz = 0;
179  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
180  if (vertexmap)
181  {
182  if (!vertexmap->empty())
183  {
184  GlobalVertex *vtx = (vertexmap->begin()->second);
185  vx = vtx->get_x();
186  vy = vtx->get_y();
187  vz = vtx->get_z();
188  }
189  }
191 
192  // cout << "vtx " << vx << " " << vy << endl;
193  // CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
194  CLHEP::Hep3Vector vertex(vx, vy, vz);
195 
196  // ------------------------------
197 
198 
199  // loop over the clusters
200  RawClusterContainer::ConstRange t_rbegin_end = recal_clusters->getClusters();
202 
203  RawCluster * savCs[10000]; // savingClusters that has 1 GeV or more
204  int iCs = 0;
205  int inCs = 0;
206 
207  // saving the clusters
208  for (t_rclusiter = t_rbegin_end.first; t_rclusiter != t_rbegin_end.second; ++t_rclusiter)
209  {
210  RawCluster *t_recalcluster = t_rclusiter->second;
211 
212  float cluse= t_recalcluster->get_energy();
213  if (cluse > 0.1) inCs++;
214 
215  if (cluse > 1.0) savCs[iCs++] = t_recalcluster;
216 
217  }
218 
219  _nClusters = iCs;
220  if (_nClusters > 60)
222 
223 
224  // looping on the saved clusters savCs[]
225  // outer loop (we want to do pair of the loops)
226  for( int jCs=0; jCs<iCs; jCs++)
227  {
228 
229  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*savCs[jCs], vertex);
230 
231  // CLHEP::Hep3Vector E_vec_cluster(PxCluster, PyCluster, ...);
232 
233 
234 
235  _clusterIDs[jCs] = savCs[jCs]->get_id();
236  _clusterEnergies[jCs] = savCs[jCs]->get_energy();
237 
238  float tt_clus_energy = E_vec_cluster.mag();
239  if (tt_clus_energy > 1.0 ) // it is greater than 1.0 no need to do if here but lets keep tho
240  {
241  float tt_clus_eta = E_vec_cluster.pseudoRapidity();
242  float tt_clus_pt = E_vec_cluster.perp();
243  float tt_clus_phi = E_vec_cluster.getPhi();
244  _clusterPts[jCs] = tt_clus_pt;
245  _clusterEtas[jCs] = tt_clus_eta;
246  _clusterPhis[jCs] = tt_clus_phi;
247 
248 
249 
250 
251  // another loop to go into the saved cluster
252  for(int kCs=0; kCs<iCs; kCs++)
253  {
254  if(jCs==kCs) continue;
255 
256  CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetEVec(*savCs[kCs], vertex);
257 
258  float tt2_clus_energy = E_vec_cluster2.mag();
259  if (tt2_clus_energy > 1.0) // again this will be greater than 1.0, but lets keep
260  {
261  // lets do alpha cut here: this is needed tho
262  alphaCut = abs(tt_clus_energy - tt2_clus_energy)/(tt_clus_energy + tt_clus_energy);
263  if (alphaCut <= 0.8)
264  {
265  float tt2_clus_eta = E_vec_cluster2.pseudoRapidity();
266  float tt2_clus_pt = E_vec_cluster2.perp();
267  float tt2_clus_phi = E_vec_cluster2.getPhi();
268 
269  TLorentzVector pho1, pho2, pi0lv;
270 
271  pho1.SetPtEtaPhiE(tt_clus_pt, tt_clus_eta,tt_clus_phi,tt_clus_energy);
272  pho2.SetPtEtaPhiE(tt2_clus_pt, tt2_clus_eta,tt2_clus_phi,tt2_clus_energy);
273 
274  if (pho1.DeltaR(pho2) > 0.8) continue;
275 
276  pi0lv=pho1+pho2;
277  float pairInvMass=pi0lv.M();
278 
279 
280 
281 
282  //vector to hold all the towers etas, phis, and energy in this cluster
283  std::vector<int> toweretas;
284  std::vector<int> towerphis;
285  std::vector<float> towerenergies;
286 
287  // loop over the towers from the outer loop cluster
288  // and find the max tower location and save the
289  // histogram on that max tower location with this
290  // invariant mass
291 
292  RawCluster::TowerConstRange towers = savCs[jCs]->get_towers();
294 
295  for (toweriter = towers.first; toweriter != towers.second; ++toweriter)
296  {
297  RawTower *tower = _towers->getTower(toweriter->first);
298 
299  int towereta = tower->get_bineta();
300  int towerphi = tower->get_binphi();
301  double towerenergy = tower->get_energy();
302 
303  // put the eta, phi, energy into corresponding vectors
304  toweretas.push_back(towereta);
305  towerphis.push_back(towerphi);
306  towerenergies.push_back(towerenergy);
307  }
308 
309  // cout << endl;
310  // cout << "Cluster energy: " << tt_clus_energy << endl;
311  // cout << "Total number of towers (getNTowers()): " << savCs[jCs]->getNTowers() << endl;
312  // cout << "Total number of towers size(toweretas): " << toweretas.size() << endl;
313  // float maxTowerEnergy = *max_element(towerenergies.begin(), towerenergies.end());
314  // cout << "The maxTowerEnergy: " << maxTowerEnergy << endl;
315  int maxTowerIndex = max_element(towerenergies.begin(), towerenergies.end()) - towerenergies.begin();
316  maxTowerEta = toweretas[maxTowerIndex];
317  maxTowerPhi = towerphis[maxTowerIndex];
318 
319  _maxTowerEtas[jCs] = maxTowerEta;
320  _maxTowerPhis[jCs] = maxTowerPhi;
321 
322 
323 
324  if (tt_clus_energy > 2.5 && tt2_clus_energy > 1.5)
325  {
326  pairInvMassTotal-> Fill(pairInvMass);
327  mass_eta->Fill(pairInvMass,tt_clus_eta); mass_eta_phi->Fill(pairInvMass,tt_clus_eta, tt_clus_phi);
328 
329  }
330 
331 
332  // fill the tower by tower histograms with invariant mass
333  cemc_hist_eta_phi[maxTowerEta][maxTowerPhi]->Fill(pairInvMass);
334  eta_hist[maxTowerEta]->Fill(pairInvMass);
335  }
336 
337  }
338 
339  }
340  }
341  }
342  _eventTree->Fill();
343 
344  // _ievent++;
345 
346  // pause few seconds
347 // std::chrono::seconds dura(2);
348 // std::this_thread::sleep_for(dura);
349 // std::cout << "Waited 2 sec\n";
350 
352 }
353 
354 //____________________________________________________________________________..
356 {
357 
358  cal_output->cd();
359 // _eventTree->Write();
360  cal_output->Write();
361  cal_output->Close();
362  delete cal_output;
363 
365 }
366 
367 
368 
369 //______________________________________________________________________________..
370 void CaloCalibEmc_Pi0::Loop(TString _filename, int nevts)
371 {
372  TFile *f = new TFile(_filename);
373 
374  TTree *t1 = (TTree*)f->Get("_eventTree");
375 
376  // Declaration of leaf types
379  Int_t _clusterIDs[10000]; //[_nClusters]
380  Float_t _clusterEnergies[10000]; //[_nClusters]
381  Float_t _clusterPts[10000];
382  Int_t _clusterEtas[10000]; //[_nClusters]
383  Int_t _clusterPhis[10000]; //[_nClusters]
384  Int_t _maxTowerEtas[10000];
385  Int_t _maxTowerPhis[10000];
386 
387  // Set Branches
388  t1->SetBranchAddress("_eventNumber", &_eventNumber);
389  t1->SetBranchAddress("_nClusters", &_nClusters);
390  t1->SetBranchAddress("_clusterIDs", _clusterIDs);
391  t1->SetBranchAddress("_clusterEnergies", _clusterEnergies);
392  t1->SetBranchAddress("_clusterPts", _clusterPts);
393  t1->SetBranchAddress("_clusterEtas", _clusterEtas);
394  t1->SetBranchAddress("_clusterPhis", _clusterPhis);
395  t1->SetBranchAddress("_maxTowerEtas", _maxTowerEtas);
396  t1->SetBranchAddress("_maxTowerPhis", _maxTowerPhis);
397 
398 
399 
400 
401  // pre-loop to save all the clusters LorentzVector
402 
403  TLorentzVector *savClusLV[10000];
404 
405  int nEntries = (int)t1->GetEntriesFast();
406  int nevts2 = nevts;
407 
408  if (nevts < 0 || nEntries < nevts)
409  nevts2 = nEntries;
410 
411  for (int i=0; i < nevts2; i++)
412  {
413  // load the ith instance of the TTree
414  t1->GetEntry(i);
415 
416  // calibration correction will be applied here
417  //
418  int nClusters = _nClusters;
419  for (int j=0; j < nClusters; j++)
420  {
421  // float px, py, pz;
422  float pt, eta, phi, E;
423  pt = _clusterPts[j];
424  eta = _clusterEtas[j];
425  phi = _clusterPhis[j];
426  E = _clusterEnergies[j];
427  // px = pt * cos(phi);
428  // py = pt * sin(phi);
429  // pz = pt * sinh(eta);
430  savClusLV[j] = new TLorentzVector();
431  savClusLV[j]->SetPtEtaPhiE(pt, eta, phi, E);
432  }
433 
434  // Next we need to have the clusterCut to only take
435  // certain clusters of certain energy
436  // But lets keep it simple for now
437 
438  TLorentzVector *pho1, *pho2;
439  int iCs = nClusters;
440  for( int jCs=0; jCs<iCs; jCs++)
441  {
442  pho1 = savClusLV[jCs];
443 
444  // another loop to go into the saved cluster
445  for(int kCs=0; kCs<iCs; kCs++)
446  {
447  if(jCs==kCs) continue;
448 
449  pho2 = savClusLV[kCs];
450 
451  TLorentzVector pi0lv;
452 
453  if (pho1->DeltaR(*pho2) > 0.8) continue;
454 
455  pi0lv=*pho1+*pho2;
456  float pairInvMass=pi0lv.M();
457  //cout << pairInvMass << endl;
458 
459 
460 
461  // fill the tower by tower histograms with invariant mass
462  cemc_hist_eta_phi[_maxTowerEtas[jCs]][_maxTowerPhis[jCs]]->Fill(pairInvMass);
463  eta_hist[_maxTowerEtas[jCs]]->Fill(pairInvMass);
464 
465  }
466 
467  }
468 
469  }
470 }
471 
472 
473 
474 // _______________________________________________________________..
476 {
477  TFile *parFile = new TFile("parFile", "RECREATE");
478 
479  TTree *parTree = new TTree("parTree", "Tree with fit mean saved");
480 
481  int ieta, iphi;
482  float corrval[100000];
483 
484  parTree->Branch("ieta", &ieta, "ieta/I");
485  parTree->Branch("iphi", &iphi, "iphi/I");
486  parTree->Branch("corrval", corrval, "corrval[ieta*1000+iphi]");
487 
488  TF1 *fit_result;
489 
490  for (int i=0; i<96; i++)
491  {
492  for (int j=0; j<=258; j++)
493  {
494  if (j!=258)
495  {
496  cemc_hist_eta_phi[i][j]->Fit("gaus", "", "", 0.095, 0.175);
497  fit_result = cemc_hist_eta_phi[i][j]->GetFunction("gaus");
498  ieta = i;
499  iphi = j;
500  corrval[i*1000+j] = fit_result->GetParameter(1);
501  }
502  else
503  {
504  eta_hist[i]->Fit("gaus", "", "", 0.095, 0.175);
505  fit_result = eta_hist[i]->GetFunction("gaus");
506  ieta = i;
507  iphi = j;
508  corrval[i*1000+j] = fit_result->GetParameter(1);
509  }
510  }
511  }
512  parTree->Fill();
513  parTree->Write();
514  parFile->Close();
515 }
516 
517