Class Reference for E1039 Core & Analysis Software
DimuAnaRUS.cc
Go to the documentation of this file.
1 #include <fstream>
2 #include <iomanip>
3 #include <TSystem.h>
4 #include <TFile.h>
5 #include <TTree.h>
6 #include <TChain.h>
7 #include <TH1D.h>
8 #include <TCanvas.h>
9 #include <interface_main/SQRun.h>
10 #include <interface_main/SQHitVector.h>
11 #include <interface_main/SQEvent.h>
12 #include <interface_main/SQMCEvent.h>
13 #include <interface_main/SQRun.h>
14 #include <interface_main/SQTrackVector.h>
15 #include <interface_main/SQDimuonVector.h>
16 #include <interface_main/SQTrackVector_v1.h>
17 #include <ktracker/FastTracklet.h>
18 #include <ktracker/SRecEvent.h>
19 #include <fun4all/Fun4AllReturnCodes.h>
20 #include <phool/PHNodeIterator.h>
21 #include <phool/PHIODataNode.h>
22 #include <phool/getClass.h>
23 #include <geom_svc/GeomSvc.h>
24 #include "DimuAnaRUS.h"
25 using namespace std;
26 
27 DimuAnaRUS::DimuAnaRUS(const std::string& name)
28  : SubsysReco (name),
29  m_file(0),
30  m_tree(0),
31  m_tree_name("tree"),
32  m_file_name("output.root"),
33  m_evt(0),
34  m_sp_map(0),
35  m_hit_vec(0),
36  m_sq_trk_vec(0),
37  m_sq_dim_vec(0),
38  saveDimuonOnly(false),
39  true_mode(false),
40  reco_dimu_mode(false),
41  reco_mode(false),
42  data_trig_mode(false),
43  mc_trig_mode(true)
44 {
45  ;
46 }
47 
48 DimuAnaRUS::~DimuAnaRUS()
49 {
50  ;
51 }
52 
53 int DimuAnaRUS::Init(PHCompositeNode* startNode)
54 {
55  return Fun4AllReturnCodes::EVENT_OK;
56 }
57 
58 int DimuAnaRUS::InitRun(PHCompositeNode* startNode)
59 {
60 
61  m_file = new TFile(m_file_name.c_str(), "RECREATE");
62 
63  if (!m_file || m_file->IsZombie()) {
64  std::cerr << "Error: Could not create file " << m_file_name << std::endl;
65  exit(1);
66  } else {
67  std::cout << "File " << m_file->GetName() << " opened successfully." << std::endl;
68  }
69 
70  m_tree = new TTree(m_tree_name.c_str(), "Tree for storing events");
71  if (!m_tree) {
72  std::cerr << "Error: Could not create tree " << m_tree_name << std::endl;
73  exit(1);
74  } else {
75  std::cout << "Tree " << m_tree->GetName() << " created successfully." << std::endl;
76  }
77 
78  m_tree->Branch("eventID", &eventID, "eventID/I");
79  m_tree->Branch("runID", &runID, "runID/I");
80  m_tree->Branch("spillID", &spillID, "spillID/I");
81  m_tree->Branch("eventID", &eventID, "eventID/I");
82  m_tree->Branch("rfID", &rfID, "rfID/I");
83  m_tree->Branch("turnID", &turnID, "turnID/I");
84  m_tree->Branch("rfIntensity", rfIntensity, "rfIntensity[33]/I");
85  m_tree->Branch("fpgaTrigger", fpgaTrigger, "fpgaTrigger[5]/I");
86  m_tree->Branch("nimTrigger", nimTrigger, "nimTrigger[5]/I");
87 
88  m_tree->Branch("hitID", &hitID);
89  m_tree->Branch("hitTrackID", &hitTrackID);
90  m_tree->Branch("processID", &processID);
91  m_tree->Branch("detectorID", &detectorID);
92  m_tree->Branch("elementID", &elementID);
93  m_tree->Branch("tdcTime", &tdcTime);
94  m_tree->Branch("driftDistance", &driftDistance);
95 
96  if (true_mode) {
97  m_tree->Branch("gCharge", &gCharge);
98  m_tree->Branch("trackID", &trackID);
99  m_tree->Branch("gvx", &gvx);
100  m_tree->Branch("gvy", &gvy);
101  m_tree->Branch("gvz", &gvz);
102  m_tree->Branch("gpx", &gpx);
103  m_tree->Branch("gpy", &gpy);
104  m_tree->Branch("gpz", &gpz);
105 
106  m_tree->Branch("gx_st1", &gx_st1);
107  m_tree->Branch("gy_st1", &gy_st1);
108  m_tree->Branch("gz_st1", &gz_st1);
109  m_tree->Branch("gpx_st1", &gpx_st1);
110  m_tree->Branch("gpy_st1", &gpy_st1);
111  m_tree->Branch("gpz_st1", &gpz_st1);
112 
113  m_tree->Branch("gx_st3", &gx_st3);
114  m_tree->Branch("gy_st3", &gy_st3);
115  m_tree->Branch("gz_st3", &gz_st3);
116  m_tree->Branch("gpx_st3", &gpx_st3);
117  m_tree->Branch("gpy_st3", &gpy_st3);
118  m_tree->Branch("gpz_st3", &gpz_st3);
119 }
120 
121  if (reco_mode==true){
122  m_tree->Branch("rec_track_id", &rec_track_id);
123  m_tree->Branch("rec_track_charge", &rec_track_charge);
124  m_tree->Branch("rec_track_vx", &rec_track_vx);
125  m_tree->Branch("rec_track_vy", &rec_track_vy);
126  m_tree->Branch("rec_track_vz", &rec_track_vz);
127  m_tree->Branch("rec_track_px", &rec_track_px);
128  m_tree->Branch("rec_track_py", &rec_track_py);
129  m_tree->Branch("rec_track_pz", &rec_track_pz);
130  m_tree->Branch("rec_track_x_st1", &rec_track_x_st1);
131  m_tree->Branch("rec_track_y_st1", &rec_track_y_st1);
132  m_tree->Branch("rec_track_z_st1", &rec_track_z_st1);
133  m_tree->Branch("rec_track_px_st1", &rec_track_px_st1);
134  m_tree->Branch("rec_track_py_st1", &rec_track_py_st1);
135  m_tree->Branch("rec_track_pz_st1", &rec_track_pz_st1);
136  m_tree->Branch("rec_track_x_st3", &rec_track_x_st3);
137  m_tree->Branch("rec_track_y_st3", &rec_track_y_st3);
138  m_tree->Branch("rec_track_z_st3", &rec_track_z_st3);
139  m_tree->Branch("rec_track_px_st3", &rec_track_px_st3);
140  m_tree->Branch("rec_track_py_st3", &rec_track_py_st3);
141  m_tree->Branch("rec_track_pz_st3", &rec_track_pz_st3);
142  m_tree->Branch("rec_track_num_hits", &rec_track_num_hits);
143  m_tree->Branch("rec_track_chisq", &rec_track_chisq);
144  m_tree->Branch("rec_track_chisq_tgt", &rec_track_chisq_tgt);
145  m_tree->Branch("rec_track_chisq_dump", &rec_track_chisq_dump);
146  m_tree->Branch("rec_track_chisq_upstream", &rec_track_chisq_upstream);
147  m_tree->Branch("rec_track_x_tgt", &rec_track_x_tgt);
148  m_tree->Branch("rec_track_y_tgt", &rec_track_y_tgt);
149  m_tree->Branch("rec_track_z_tgt", &rec_track_z_tgt);
150  m_tree->Branch("rec_track_px_tgt", &rec_track_px_tgt);
151  m_tree->Branch("rec_track_py_tgt", &rec_track_py_tgt);
152  m_tree->Branch("rec_track_pz_tgt", &rec_track_pz_tgt);
153  m_tree->Branch("rec_track_x_dump", &rec_track_x_dump);
154  m_tree->Branch("rec_track_y_dump", &rec_track_y_dump);
155  m_tree->Branch("rec_track_z_dump", &rec_track_z_dump);
156  m_tree->Branch("rec_track_px_dump", &rec_track_px_dump);
157  m_tree->Branch("rec_track_py_dump", &rec_track_py_dump);
158  m_tree->Branch("rec_track_pz_dump", &rec_track_pz_dump);
159  m_tree->Branch("rec_hit_ids", &rec_hit_ids);
160  m_tree->Branch("rec_track_hit_x", &rec_track_hit_x);
161  m_tree->Branch("rec_track_hit_y", &rec_track_hit_y);
162 
163 
164  if (reco_dimu_mode == true) {
165  m_tree->Branch("rec_dimuon_id", &rec_dimuon_id);
166  m_tree->Branch("rec_dimuon_true_id", &rec_dimuon_true_id);
167  m_tree->Branch("rec_dimuon_track_id_pos", &rec_dimuon_track_id_pos);
168  m_tree->Branch("rec_dimuon_track_id_neg", &rec_dimuon_track_id_neg);
169 
170  m_tree->Branch("rec_dimuon_x", &rec_dimuon_x);
171  m_tree->Branch("rec_dimuon_y", &rec_dimuon_y);
172  m_tree->Branch("rec_dimuon_z", &rec_dimuon_z);
173 
174  m_tree->Branch("rec_dimuon_px_pos", &rec_dimuon_px_pos);
175  m_tree->Branch("rec_dimuon_py_pos", &rec_dimuon_py_pos);
176  m_tree->Branch("rec_dimuon_pz_pos", &rec_dimuon_pz_pos);
177 
178  m_tree->Branch("rec_dimuon_px_neg", &rec_dimuon_px_neg);
179  m_tree->Branch("rec_dimuon_py_neg", &rec_dimuon_py_neg);
180  m_tree->Branch("rec_dimuon_pz_neg", &rec_dimuon_pz_neg);
181 
182  m_tree->Branch("rec_dimuon_px_pos_tgt", &rec_dimuon_px_pos_tgt);
183  m_tree->Branch("rec_dimuon_py_pos_tgt", &rec_dimuon_py_pos_tgt);
184  m_tree->Branch("rec_dimuon_pz_pos_tgt", &rec_dimuon_pz_pos_tgt);
185 
186  m_tree->Branch("rec_dimuon_px_neg_tgt", &rec_dimuon_px_neg_tgt);
187  m_tree->Branch("rec_dimuon_py_neg_tgt", &rec_dimuon_py_neg_tgt);
188  m_tree->Branch("rec_dimuon_pz_neg_tgt", &rec_dimuon_pz_neg_tgt);
189 
190  m_tree->Branch("rec_dimuon_px_pos_dump", &rec_dimuon_px_pos_dump);
191  m_tree->Branch("rec_dimuon_py_pos_dump", &rec_dimuon_py_pos_dump);
192  m_tree->Branch("rec_dimuon_pz_pos_dump", &rec_dimuon_pz_pos_dump);
193 
194  m_tree->Branch("rec_dimuon_px_neg_dump", &rec_dimuon_px_neg_dump);
195  m_tree->Branch("rec_dimuon_py_neg_dump", &rec_dimuon_py_neg_dump);
196  m_tree->Branch("rec_dimuon_pz_neg_dump", &rec_dimuon_pz_neg_dump);
197 
198  }
199  }
200 
201 
202  if (reco_mode) {
203  m_sq_trk_vec = findNode::getClass<SQTrackVector>(startNode, "SQRecTrackVector");
204  if(reco_dimu_mode==true) m_sq_dim_vec = findNode::getClass<SQDimuonVector>(startNode, "SQRecDimuonVector_PM");
205 
206  if (!m_sq_trk_vec) {
207  std::cerr << "Error: m_sq_trk_vec is null!" << std::endl;
208  return Fun4AllReturnCodes::ABORTEVENT;
209  }
210 
211  if(reco_dimu_mode==true){
212  if (!m_sq_dim_vec) {
213  std::cerr << "Error: m_sq_dim_vec is null!" << std::endl;
214  return Fun4AllReturnCodes::ABORTEVENT;
215  }
216  }
217  }
218  m_evt = findNode::getClass<SQEvent>(startNode, "SQEvent");
219 
220  if(true_mode == true){
221  m_hit_vec = findNode::getClass<SQHitVector>(startNode, "SQHitVector");
222  m_vec_trk = findNode::getClass<SQTrackVector>(startNode, "SQTruthTrackVector");
223  if (!m_evt || !m_hit_vec || !m_vec_trk) {
224  return Fun4AllReturnCodes::ABORTEVENT;
225  }
226  }
227  else {
228  m_hit_vec = findNode::getClass<SQHitVector>(startNode, "SQHitVector");
229  if (!m_evt || !m_hit_vec) {
230  return Fun4AllReturnCodes::ABORTEVENT;
231  }
232  }
233 
234  if(reco_mode ==true && data_trig_mode ==true){
235  //cout << "inside the data roadset mode: "<<endl;
236  SQRun* sq_run = findNode::getClass<SQRun>(startNode, "SQRun");
237  if (!sq_run) std::cout << "Error: SQRun is null!" << std::endl;
238  if (!sq_run) return Fun4AllReturnCodes::ABORTEVENT;
239  int LBtop = sq_run->get_v1495_id(2);
240  int LBbot = sq_run->get_v1495_id(3);
241  int ret = m_rs.LoadConfig(LBtop, LBbot);
242  if (ret != 0) {
243  cout << "!!WARNING!! OnlMonTrigEP::InitRunOnlMon(): roadset.LoadConfig returned " << ret << ".\n";
244  }
245  cout <<"Roadset " << m_rs.str(1) << endl;
246  }
247 
248  if(reco_mode ==true && mc_trig_mode ==true){
249  int ret = m_rs.LoadConfig(131);
250  if (ret != 0) {
251  cout << "!!WARNING!! OnlMonTrigEP::InitRunOnlMon(): roadset.LoadConfig returned " << ret << ".\n";
252  }
253  //cout <<"Roadset " << m_rs.str(1) << endl;
254  }
255 
256  return Fun4AllReturnCodes::EVENT_OK;
257 }
258 
259 
260 int DimuAnaRUS::process_event(PHCompositeNode* startNode)
261 {
262 
263  int proc_id=0;
264  int SourceFlag=0;
265  eventID = m_evt->get_event_id();
266 
267  runID = m_evt->get_run_id();
268  spillID = m_evt->get_spill_id();
269  rfID = m_evt->get_qie_rf_id();
270  turnID = m_evt->get_qie_turn_id();
271 
272  fpgaTrigger[0] = m_evt->get_trigger(SQEvent::MATRIX1);
273  fpgaTrigger[1] = m_evt->get_trigger(SQEvent::MATRIX2);
274  fpgaTrigger[2] = m_evt->get_trigger(SQEvent::MATRIX3);
275  fpgaTrigger[3] = m_evt->get_trigger(SQEvent::MATRIX4);
276  fpgaTrigger[4] = m_evt->get_trigger(SQEvent::MATRIX5);
277 
278  nimTrigger[0] = m_evt->get_trigger(SQEvent::NIM1);
279  nimTrigger[1] = m_evt->get_trigger(SQEvent::NIM2);
280  nimTrigger[2] = m_evt->get_trigger(SQEvent::NIM3);
281  nimTrigger[3] = m_evt->get_trigger(SQEvent::NIM4);
282  nimTrigger[4] = m_evt->get_trigger(SQEvent::NIM5);
283  for (int i = -16; i <= 16; ++i) {
284  rfIntensity[i + 16] = m_evt->get_qie_rf_intensity(i);
285  }
286 
287 
288  if (m_hit_vec && true_mode==false) {
289  ResetHitBranches();
290  for (int ihit = 0; ihit < m_hit_vec->size(); ++ihit) {
291  SQHit* hit = m_hit_vec->at(ihit);
292  hitTrackID.push_back(hit->get_track_id());
293  hitID.push_back(hit->get_hit_id());
294  detectorID.push_back(hit->get_detector_id());
295  elementID.push_back(hit->get_element_id());
296  tdcTime.push_back(hit->get_tdc_time());
297  driftDistance.push_back(hit->get_drift_distance());
298  }
299  }
300 
301  const double muon_mass = 0.10566;
302  TLorentzVector mu1;
303  TLorentzVector mu2;
304  if(true_mode == true){
305  ResetTrueBranches();
306  ResetHitBranches();
307  for (unsigned int ii = 0; ii < m_vec_trk->size(); ii++) {
308 
309  SQTrack* trk = m_vec_trk->at(ii);
310  gCharge.push_back(trk->get_charge());
311  trackID.push_back(trk->get_track_id());
312  gvx.push_back(trk->get_pos_vtx().X());
313  gvy.push_back(trk->get_pos_vtx().Y());
314  gvz.push_back(trk->get_pos_vtx().Z());
315  gpx.push_back(trk->get_mom_vtx().Px());
316  gpy.push_back(trk->get_mom_vtx().Py());
317  gpz.push_back(trk->get_mom_vtx().Pz());
318 
319  gx_st1.push_back(trk->get_pos_st1().X());
320  gy_st1.push_back(trk->get_pos_st1().Y());
321  gz_st1.push_back(trk->get_pos_st1().Z());
322  gpx_st1.push_back(trk->get_mom_st1().Px());
323  gpy_st1.push_back(trk->get_mom_st1().Py());
324  gpz_st1.push_back(trk->get_mom_st1().Pz());
325 
326  gx_st3.push_back(trk->get_pos_st3().X());
327  gy_st3.push_back(trk->get_pos_st3().Y());
328  gz_st3.push_back(trk->get_pos_st3().Z());
329  gpx_st3.push_back(trk->get_mom_st3().Px());
330  gpy_st3.push_back(trk->get_mom_st3().Py());
331  gpz_st3.push_back(trk->get_mom_st3().Pz());
332 
333  int SourceFlag=0;
334  double z = trk->get_pos_vtx().Z();
335  if (z <= -296. && z >= -304.) SourceFlag = 1; // target
336  else if (z >= 0. && z < 250) SourceFlag = 2; // dump
337  else if (z > -296. && z < 0.) SourceFlag = 3; // gap
338  else if (z < -304) SourceFlag = 0; // upstream
339 
340  if (m_hit_vec) {
341  for (int ihit = 0; ihit < m_hit_vec->size(); ++ihit) {
342  SQHit* hit = m_hit_vec->at(ihit);
343  //if (!hit->is_in_time()) continue;
344  proc_id=12;
345  if(hit->get_track_id() != trk->get_track_id()) continue;
346  int processID_;
347  if(trk->get_charge() >0) processID_ =proc_id;
348  else processID_ = proc_id +10 ;
349  int sourceFlag_= SourceFlag;
350  unsigned int encodedValue = EncodeProcess(processID_, sourceFlag_);
351  //cout << "encodedValue: "<< encodedValue<< endl;
352  //cout << "DecodeSourceFlag: "<< DecodeSourceFlag(encodedValue) << endl;
353  //cout << "DecodeProcessID: "<< DecodeProcessID(encodedValue) << endl;
354  //cout << "det Name: "<< hit->get_detector_id() << " Hit Id: "<< hit->get_hit_id () << endl;
355  processID.push_back(encodedValue);
356  hitID.push_back(hit->get_hit_id());
357  hitTrackID.push_back(hit->get_track_id());
358  detectorID.push_back(hit->get_detector_id());
359  elementID.push_back(hit->get_element_id());
360  tdcTime.push_back(hit->get_tdc_time());
361  driftDistance.push_back(hit->get_drift_distance());
362  }
363  }
364  }
365  }
366 
367 int index = -1;
368 
369  if(reco_mode == true){
370  ResetRecoBranches();
371 
372  for (auto it = m_sq_trk_vec->begin(); it != m_sq_trk_vec->end(); ++it) {
373  index+=1;
374  SRecTrack* trk = dynamic_cast<SRecTrack*>(*it);
375 
376  //cout << "track id: "<< trk->get_rec_track_id () <<" charge: " << trk->get_charge() << "index: "<<endl;
377 
378  rec_track_id.push_back(index); // where i_trk is the index in the reco track loop
379  // Basic track info
380  rec_track_charge.push_back(trk->get_charge());
381  rec_track_vx.push_back(trk->get_pos_vtx().X());
382  rec_track_vy.push_back(trk->get_pos_vtx().Y());
383  rec_track_vz.push_back(trk->get_pos_vtx().Z());
384  rec_track_px.push_back(trk->get_mom_vtx().Px());
385  rec_track_py.push_back(trk->get_mom_vtx().Py());
386  rec_track_pz.push_back(trk->get_mom_vtx().Pz());
387 
388  // Station 1
389  rec_track_x_st1.push_back(trk->get_pos_st1().X());
390  rec_track_y_st1.push_back(trk->get_pos_st1().Y());
391  rec_track_z_st1.push_back(trk->get_pos_st1().Z());
392  rec_track_px_st1.push_back(trk->get_mom_st1().Px());
393  rec_track_py_st1.push_back(trk->get_mom_st1().Py());
394  rec_track_pz_st1.push_back(trk->get_mom_st1().Pz());
395 
396  // Station 3
397  rec_track_x_st3.push_back(trk->get_pos_st3().X());
398  rec_track_y_st3.push_back(trk->get_pos_st3().Y());
399  rec_track_z_st3.push_back(trk->get_pos_st3().Z());
400  rec_track_px_st3.push_back(trk->get_mom_st3().Px());
401  rec_track_py_st3.push_back(trk->get_mom_st3().Py());
402  rec_track_pz_st3.push_back(trk->get_mom_st3().Pz());
403 
404  // Chi-squared
405  rec_track_chisq.push_back(trk->get_chisq());
406  rec_track_chisq_tgt.push_back(trk->getChisqTarget());
407  rec_track_chisq_dump.push_back(trk->get_chisq_dump());
408  rec_track_chisq_upstream.push_back(trk->get_chisq_upstream());
409 
410  // Number of hits
411  rec_track_num_hits.push_back(trk->get_num_hits());
412 
413  // Target region
414  rec_track_x_tgt.push_back(trk->get_pos_target().X());
415  rec_track_y_tgt.push_back(trk->get_pos_target().Y());
416  rec_track_z_tgt.push_back(trk->get_pos_target().Z());
417  rec_track_px_tgt.push_back(trk->get_mom_target().Px());
418  rec_track_py_tgt.push_back(trk->get_mom_target().Py());
419  rec_track_pz_tgt.push_back(trk->get_mom_target().Pz());
420 
421  // Dump region
422  rec_track_x_dump.push_back(trk->get_pos_dump().X());
423  rec_track_y_dump.push_back(trk->get_pos_dump().Y());
424  rec_track_z_dump.push_back(trk->get_pos_dump().Z());
425  rec_track_px_dump.push_back(trk->get_mom_dump().Px());
426  rec_track_py_dump.push_back(trk->get_mom_dump().Py());
427  rec_track_pz_dump.push_back(trk->get_mom_dump().Pz());
428 
429  // Now fill hit IDs for this track
430  std::vector<int> hit_ids;
431  std::vector<double> hit_ids_pos_x;
432  std::vector<double> hit_ids_pos_y;
433  for (int i_hit = 0; i_hit < trk->get_num_hits(); ++i_hit) {
434  //cout <<"charge: "<< trk->get_charge()<< " hit id of the tracks: " << trk->get_hit_id(i_hit) << endl;
435  hit_ids.push_back(trk->get_hit_id(i_hit));
436  auto [det_id, det_z_pos] = GetDetElemIDFromHitID(trk->get_hit_id(i_hit));
437  Double_t x = 0.0, y = 0.0;
438  trk->getExpPositionFast(static_cast<Double_t>(det_z_pos), x, y);
439  //cout << "x, y pos: (" << x << ", " << y << ")" << endl;
440  hit_ids_pos_x.push_back(x);
441  hit_ids_pos_y.push_back(y);
442  }
443  rec_hit_ids.push_back(hit_ids);
444  rec_track_hit_x.push_back(hit_ids_pos_x);
445  rec_track_hit_y.push_back(hit_ids_pos_y);
446  }
447 
448  if(reco_dimu_mode==true){
449  ResetRecoDimuBranches();
450  for (auto it = m_sq_dim_vec->begin(); it != m_sq_dim_vec->end(); it++) {
451  SRecDimuon& sdim = dynamic_cast<SRecDimuon&>(**it);
452  int trk_id_pos = sdim.get_track_id_pos();
453  int trk_id_neg = sdim.get_track_id_neg();
454  SRecTrack& trk_pos = dynamic_cast<SRecTrack&>(*(m_sq_trk_vec->at(trk_id_pos)));
455  SRecTrack& trk_neg = dynamic_cast<SRecTrack&>(*(m_sq_trk_vec->at(trk_id_neg)));
456 
457 
458  rec_dimuon_id.push_back(sdim.get_dimuon_id());
459  rec_dimuon_true_id.push_back(sdim.get_rec_dimuon_id());
460  rec_dimuon_track_id_pos.push_back(trk_id_pos);
461  rec_dimuon_track_id_neg.push_back(trk_id_neg);
462  rec_dimuon_px_pos.push_back(sdim.get_mom_pos().Px());
463  rec_dimuon_py_pos.push_back(sdim.get_mom_pos().Py());
464  rec_dimuon_pz_pos.push_back(sdim.get_mom_pos().Pz());
465  rec_dimuon_px_neg.push_back(sdim.get_mom_neg().Px());
466  rec_dimuon_py_neg.push_back(sdim.get_mom_neg().Py());
467  rec_dimuon_pz_neg.push_back(sdim.get_mom_neg().Pz());
468  rec_dimuon_x.push_back(sdim.get_pos().X());
469  rec_dimuon_y.push_back(sdim.get_pos().Y());
470  rec_dimuon_z.push_back(sdim.get_pos().Z());
471  // ===== Target hypothesis =====
472  sdim.calcVariables(1); // 1 = target
473  rec_dimuon_px_pos_tgt.push_back(sdim.p_pos_target.Px());
474  rec_dimuon_py_pos_tgt.push_back(sdim.p_pos_target.Py());
475  rec_dimuon_pz_pos_tgt.push_back(sdim.p_pos_target.Pz());
476  rec_dimuon_px_neg_tgt.push_back(sdim.p_neg_target.Px());
477  rec_dimuon_py_neg_tgt.push_back(sdim.p_neg_target.Py());
478  rec_dimuon_pz_neg_tgt.push_back(sdim.p_neg_target.Pz());
479  // ===== Dump hypothesis =====
480  sdim.calcVariables(2); // 2 = dump
481  rec_dimuon_px_pos_dump.push_back(sdim.p_pos_dump.Px());
482  rec_dimuon_py_pos_dump.push_back(sdim.p_pos_dump.Py());
483  rec_dimuon_pz_pos_dump.push_back(sdim.p_pos_dump.Pz());
484  rec_dimuon_px_neg_dump.push_back(sdim.p_neg_dump.Px());
485  rec_dimuon_py_neg_dump.push_back(sdim.p_neg_dump.Py());
486  rec_dimuon_pz_neg_dump.push_back(sdim.p_neg_dump.Pz());
487  }
488  }
489  }
490  m_tree->Fill();
491  return Fun4AllReturnCodes::EVENT_OK;
492 }
493 
494 int DimuAnaRUS::End(PHCompositeNode* startNode)
495 {
496  m_file->cd();
497  m_file->Write();
498  m_file->Close();
499  return Fun4AllReturnCodes::EVENT_OK;
500 }
501 
502 
503 void DimuAnaRUS::ResetHitBranches() {
504  hitID.clear();
505  processID.clear();
506  hitTrackID.clear();
507  detectorID.clear();
508  elementID.clear();
509  tdcTime.clear();
510  driftDistance.clear();
511 }
512 void DimuAnaRUS::ResetTrueBranches() {
513  gCharge.clear(); trackID.clear();
514  gvx.clear(); gvy.clear(); gvz.clear();
515  gpx.clear(); gpy.clear(); gpz.clear();
516  gx_st1.clear(); gy_st1.clear(); gz_st1.clear();
517  gpx_st1.clear(); gpy_st1.clear(); gpz_st1.clear();
518  gx_st3.clear(); gy_st3.clear(); gz_st3.clear();
519  gpx_st3.clear(); gpy_st3.clear(); gpz_st3.clear();
520 }
521 
522 void DimuAnaRUS::ResetRecoBranches() {
523  rec_track_id.clear();
524  rec_track_charge.clear();
525  rec_track_vx.clear(); rec_track_vy.clear(); rec_track_vz.clear();
526  rec_track_px.clear(); rec_track_py.clear(); rec_track_pz.clear();
527  rec_track_x_st1.clear(); rec_track_y_st1.clear(); rec_track_z_st1.clear();
528  rec_track_px_st1.clear(); rec_track_py_st1.clear(); rec_track_pz_st1.clear();
529  rec_track_x_st3.clear(); rec_track_y_st3.clear(); rec_track_z_st3.clear();
530  rec_track_px_st3.clear(); rec_track_py_st3.clear(); rec_track_pz_st3.clear();
531  rec_track_chisq.clear();
532  rec_track_chisq_tgt.clear();
533  rec_track_chisq_dump.clear();
534  rec_track_chisq_upstream.clear();
535  rec_track_num_hits.clear();
536  rec_track_x_tgt.clear(); rec_track_y_tgt.clear(); rec_track_z_tgt.clear();
537  rec_track_px_tgt.clear(); rec_track_py_tgt.clear(); rec_track_pz_tgt.clear();
538  rec_track_x_dump.clear(); rec_track_y_dump.clear(); rec_track_z_dump.clear();
539  rec_track_px_dump.clear(); rec_track_py_dump.clear(); rec_track_pz_dump.clear();
540  rec_hit_ids.clear();
541  rec_track_hit_x.clear();
542  rec_track_hit_y.clear();
543 }
544 
545 void DimuAnaRUS::ResetRecoDimuBranches() {
546  rec_dimuon_id.clear(); rec_dimuon_true_id.clear(); rec_dimuon_track_id_pos.clear(); rec_dimuon_track_id_neg.clear();
547  rec_dimuon_x.clear(); rec_dimuon_y.clear(); rec_dimuon_z.clear();
548  rec_dimuon_px_pos.clear(); rec_dimuon_py_pos.clear(); rec_dimuon_pz_pos.clear();
549  rec_dimuon_px_neg.clear(); rec_dimuon_py_neg.clear(); rec_dimuon_pz_neg.clear();
550  rec_dimuon_px_pos_tgt.clear(); rec_dimuon_py_pos_tgt.clear(); rec_dimuon_pz_pos_tgt.clear();
551  rec_dimuon_px_neg_tgt.clear(); rec_dimuon_py_neg_tgt.clear(); rec_dimuon_pz_neg_tgt.clear();
552  rec_dimuon_px_pos_dump.clear(); rec_dimuon_py_pos_dump.clear(); rec_dimuon_pz_pos_dump.clear();
553  rec_dimuon_px_neg_dump.clear(); rec_dimuon_py_neg_dump.clear(); rec_dimuon_pz_neg_dump.clear();
554 }
555 
556 
557 
558 unsigned int DimuAnaRUS::EncodeProcess(int processID, int sourceFlag) {
559  return ( (sourceFlag & 0x3) << 5 ) | // 2 bits for sourceFlag (1-3)
560  (processID & 0x1F); // 5 bits for processID (11-24)
561  }
562 int DimuAnaRUS::DecodeSourceFlag(unsigned int encoded) {
563  return (encoded >> 5) & 0x3;
564 }
565 int DimuAnaRUS::DecodeProcessID(unsigned int encoded) {
566  return encoded & 0x1F;
567 }
568 
569 std::pair<int, int> DimuAnaRUS::GetDetElemIDFromHitID(int hit_id) const {
570  if (!m_hit_vec) return {-1, -1};
571 
572  for (const auto* hit : *m_hit_vec) {
573  if (hit->get_hit_id() == hit_id) {
574  GeomSvc::UseDbSvc(true);
575  GeomSvc* geom = GeomSvc::instance();
576  return {hit->get_detector_id(), geom->getDetectorZ0(geom->getDetectorName(hit->get_detector_id()))};
577  }
578  }
579 
580  return {-1, -1}; // Not found
581 }
582 
SQTrackVector_v1.h
DimuAnaRUS::DecodeProcessID
static int DecodeProcessID(unsigned int encoded)
Definition: DimuAnaRUS.cc:565
DimuAnaRUS::EncodeProcess
unsigned int EncodeProcess(int processID, int sourceFlag)
Definition: DimuAnaRUS.cc:558
DimuAnaRUS::process_event
int process_event(PHCompositeNode *startNode)
Definition: DimuAnaRUS.cc:260
DimuAnaRUS::InitRun
int InitRun(PHCompositeNode *startNode)
Definition: DimuAnaRUS.cc:58
DimuAnaRUS::reco_dimu_mode
bool reco_dimu_mode
Definition: DimuAnaRUS.h:49
DimuAnaRUS::End
int End(PHCompositeNode *startNode)
Called at the end of all processing.
Definition: DimuAnaRUS.cc:494
DimuAnaRUS::GetDetElemIDFromHitID
std::pair< int, int > GetDetElemIDFromHitID(int hit_id) const
Definition: DimuAnaRUS.cc:569
DimuAnaRUS::DimuAnaRUS
DimuAnaRUS(const std::string &name="DimuAnaRUS")
Definition: DimuAnaRUS.cc:27
DimuAnaRUS::proc_id
int proc_id
Definition: DimuAnaRUS.h:47
DimuAnaRUS::ResetTrueBranches
void ResetTrueBranches()
Definition: DimuAnaRUS.cc:512
DimuAnaRUS::mc_trig_mode
bool mc_trig_mode
Definition: DimuAnaRUS.h:44
DimuAnaRUS::~DimuAnaRUS
virtual ~DimuAnaRUS()
Definition: DimuAnaRUS.cc:48
DimuAnaRUS::data_trig_mode
bool data_trig_mode
Definition: DimuAnaRUS.h:43
DimuAnaRUS::m_rs
UtilTrigger::TrigRoadset m_rs
Definition: DimuAnaRUS.h:40
DimuAnaRUS::ResetRecoBranches
void ResetRecoBranches()
Definition: DimuAnaRUS.cc:522
DimuAnaRUS::ResetHitBranches
void ResetHitBranches()
Definition: DimuAnaRUS.cc:503
DimuAnaRUS::ResetRecoDimuBranches
void ResetRecoDimuBranches()
Definition: DimuAnaRUS.cc:545
DimuAnaRUS::true_mode
bool true_mode
Definition: DimuAnaRUS.h:41
DimuAnaRUS::reco_mode
bool reco_mode
Definition: DimuAnaRUS.h:42
DimuAnaRUS::DecodeSourceFlag
static int DecodeSourceFlag(unsigned int encoded)
Definition: DimuAnaRUS.cc:562
DimuAnaRUS::SourceFlag
int SourceFlag
Definition: DimuAnaRUS.h:48
DimuAnaRUS::Init
int Init(PHCompositeNode *startNode)
Definition: DimuAnaRUS.cc:53
GeomSvc
User interface class about the geometry of detector planes.
Definition: GeomSvc.h:164
GeomSvc::instance
static GeomSvc * instance()
singlton instance
Definition: GeomSvc.cxx:212
GeomSvc::UseDbSvc
static bool UseDbSvc()
Definition: GeomSvc.h:352
GeomSvc::getDetectorZ0
double getDetectorZ0(const std::string detectorName) const
Definition: GeomSvc.h:213
GeomSvc::getDetectorName
std::string getDetectorName(const int &detectorID) const
Definition: GeomSvc.h:223
SQDimuonVector::begin
virtual ConstIter begin() const =0
SQDimuonVector::end
virtual ConstIter end() const =0
SQEvent::get_run_id
virtual int get_run_id() const =0
Return the run ID.
SQEvent::MATRIX2
@ MATRIX2
Definition: SQEvent.h:28
SQEvent::NIM5
@ NIM5
Definition: SQEvent.h:26
SQEvent::MATRIX1
@ MATRIX1
Definition: SQEvent.h:27
SQEvent::NIM4
@ NIM4
Definition: SQEvent.h:25
SQEvent::NIM2
@ NIM2
Definition: SQEvent.h:23
SQEvent::NIM1
@ NIM1
Definition: SQEvent.h:22
SQEvent::MATRIX3
@ MATRIX3
Definition: SQEvent.h:29
SQEvent::NIM3
@ NIM3
Definition: SQEvent.h:24
SQEvent::MATRIX4
@ MATRIX4
Definition: SQEvent.h:30
SQEvent::MATRIX5
@ MATRIX5
Definition: SQEvent.h:31
SQEvent::get_qie_rf_intensity
virtual int get_qie_rf_intensity(const short i) const =0
Return the i-th QIE RF intensity, where i=-16...+16.
SQEvent::get_trigger
virtual bool get_trigger(const SQEvent::TriggerMask i) const =0
Return the trigger bit (fired or not) of the selected trigger channel.
SQEvent::get_qie_turn_id
virtual int get_qie_turn_id() const =0
Return the QIE turn ID.
SQEvent::get_qie_rf_id
virtual int get_qie_rf_id() const =0
Return the QIE RF ID.
SQEvent::get_spill_id
virtual int get_spill_id() const =0
Return the spill ID.
SQEvent::get_event_id
virtual int get_event_id() const =0
Return the event ID, which is unique per run.
SQHitVector::at
virtual const SQHit * at(const size_t idkey) const =0
SQHitVector::size
virtual size_t size() const =0
SQHit
An SQ interface class to hold one detector hit.
Definition: SQHit.h:20
SQHit::get_drift_distance
virtual float get_drift_distance() const
Return the drift distance of this hit. Probably the value is not properly set at present....
Definition: SQHit.h:57
SQHit::get_element_id
virtual short get_element_id() const
Return the element ID of this hit.
Definition: SQHit.h:45
SQHit::get_hit_id
virtual int get_hit_id() const
Return the ID of this hit.
Definition: SQHit.h:39
SQHit::get_tdc_time
virtual float get_tdc_time() const
Return the TDC time (nsec) of this hit.
Definition: SQHit.h:54
SQHit::get_detector_id
virtual short get_detector_id() const
Return the detector ID of this hit.
Definition: SQHit.h:42
SQHit::get_track_id
virtual int get_track_id() const
Return the track ID associated with this hit. Probably the value is not properly set at present.
Definition: SQHit.h:66
SQRun
An SQ interface class to hold the run-level info.
Definition: SQRun.h:18
SQRun::get_v1495_id
virtual int get_v1495_id(const int chan) const
Return the ID of the chan-th V1495 boards.
Definition: SQRun.h:53
SQTrackVector::begin
virtual ConstIter begin() const =0
SQTrackVector::at
virtual const SQTrack * at(const size_t id) const =0
SQTrackVector::end
virtual ConstIter end() const =0
SQTrackVector::size
virtual size_t size() const =0
SQTrack
An SQ interface class to hold one true or reconstructed track.
Definition: SQTrack.h:8
SQTrack::get_track_id
virtual int get_track_id() const =0
Return the track ID, which is unique per event(?).
SQTrack::get_pos_vtx
virtual TVector3 get_pos_vtx() const =0
Return the track position at vertex.
SQTrack::get_charge
virtual int get_charge() const =0
Return the charge, i.e. +1 or -1.
SQTrack::get_mom_vtx
virtual TLorentzVector get_mom_vtx() const =0
Return the track momentum at vertex.
SQTrack::get_mom_st1
virtual TLorentzVector get_mom_st1() const =0
Return the track momentum at Station 1.
SQTrack::get_pos_st3
virtual TVector3 get_pos_st3() const =0
Return the track position at Station 3.
SQTrack::get_pos_st1
virtual TVector3 get_pos_st1() const =0
Return the track position at Station 1.
SQTrack::get_mom_st3
virtual TLorentzVector get_mom_st3() const =0
Return the track momentum at Station 3.
SRecDimuon::get_rec_dimuon_id
virtual int get_rec_dimuon_id() const
Return the dimuon ID of associated reconstructed dimuon. Valid only if this object holds truth dimuon...
Definition: SRecEvent.h:314
SRecDimuon::get_dimuon_id
virtual int get_dimuon_id() const
SQDimuon virtual functions.
Definition: SRecEvent.h:311
SRecDimuon::calcVariables
void calcVariables(int choice=0)
Calculate the kinematic vairables, 0 = vertex, 1 = target, 2 = dump.
Definition: SRecEvent.cxx:633
SRecDimuon::p_neg_target
TLorentzVector p_neg_target
Definition: SRecEvent.h:370
SRecDimuon::get_pos
virtual TVector3 get_pos() const
Return the dimuon position at vertex.
Definition: SRecEvent.h:326
SRecDimuon::p_pos_target
TLorentzVector p_pos_target
Track momentum projections at different location.
Definition: SRecEvent.h:369
SRecDimuon::get_track_id_neg
virtual int get_track_id_neg() const
Return the track ID of the negative track.
Definition: SRecEvent.h:323
SRecDimuon::get_track_id_pos
virtual int get_track_id_pos() const
Return the track ID of the positive track.
Definition: SRecEvent.h:320
SRecDimuon::p_pos_dump
TLorentzVector p_pos_dump
Definition: SRecEvent.h:371
SRecDimuon::get_mom_neg
virtual TLorentzVector get_mom_neg() const
Return the momentum of the negative track at vertex.
Definition: SRecEvent.h:335
SRecDimuon::p_neg_dump
TLorentzVector p_neg_dump
Definition: SRecEvent.h:372
SRecDimuon::get_mom_pos
virtual TLorentzVector get_mom_pos() const
Return the momentum of the positive track at vertex.
Definition: SRecEvent.h:332
SRecTrack::get_pos_st3
virtual TVector3 get_pos_st3() const
Return the track position at Station 3.
Definition: SRecEvent.h:72
SRecTrack::getChisqTarget
Double_t getChisqTarget()
Definition: SRecEvent.h:199
SRecTrack::get_mom_st3
virtual TLorentzVector get_mom_st3() const
Return the track momentum at Station 3.
Definition: SRecEvent.h:81
SRecTrack::get_charge
virtual int get_charge() const
Return the charge, i.e. +1 or -1.
Definition: SRecEvent.h:60
SRecTrack::get_mom_target
virtual TLorentzVector get_mom_target() const
Definition: SRecEvent.h:93
SRecTrack::getExpPositionFast
void getExpPositionFast(Double_t z, Double_t &x, Double_t &y, Int_t iNode=-1)
Definition: SRecEvent.cxx:267
SRecTrack::get_pos_st1
virtual TVector3 get_pos_st1() const
Return the track position at Station 1.
Definition: SRecEvent.h:69
SRecTrack::get_hit_id
virtual int get_hit_id(const int i) const
Definition: SRecEvent.h:96
SRecTrack::get_chisq_upstream
virtual double get_chisq_upstream() const
Definition: SRecEvent.h:88
SRecTrack::get_mom_st1
virtual TLorentzVector get_mom_st1() const
Return the track momentum at Station 1.
Definition: SRecEvent.h:78
SRecTrack::get_pos_dump
virtual TVector3 get_pos_dump() const
Definition: SRecEvent.h:91
SRecTrack::get_mom_dump
virtual TLorentzVector get_mom_dump() const
Definition: SRecEvent.h:94
SRecTrack::get_num_hits
virtual int get_num_hits() const
Return the number of hits associated to this track.
Definition: SRecEvent.h:63
SRecTrack::get_pos_target
virtual TVector3 get_pos_target() const
Definition: SRecEvent.h:90
SRecTrack::get_chisq_dump
virtual double get_chisq_dump() const
Definition: SRecEvent.h:86
SRecTrack::get_mom_vtx
virtual TLorentzVector get_mom_vtx() const
Return the track momentum at vertex.
Definition: SRecEvent.h:75
SRecTrack::get_chisq
virtual double get_chisq() const
Definition: SRecEvent.h:84
SRecTrack::get_pos_vtx
virtual TVector3 get_pos_vtx() const
Return the track position at vertex.
Definition: SRecEvent.h:66
UtilTrigger::TrigRoadset::LoadConfig
int LoadConfig(const std::string dir)
Definition: TrigRoadset.cc:23
UtilTrigger::TrigRoadset::str
std::string str(const int level=0) const
Definition: TrigRoadset.cc:79