Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SigmaTimingNtuple.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SigmaTimingNtuple.cc
1 #include "SigmaTimingNtuple.h"
2 
3 #include <g4main/PHG4HitContainer.h>
4 #include <g4main/PHG4Hit.h>
5 
6 #include <g4main/PHG4TruthInfoContainer.h>
7 #include <g4main/PHG4Particle.h>
8 #include <g4main/PHG4VtxPoint.h>
9 
10 #include <fun4all/Fun4AllHistoManager.h>
11 #include <fun4all/Fun4AllReturnCodes.h>
12 
13 #include <phool/getClass.h>
14 #include <phool/PHRandomSeed.h>
15 
16 #include <TFile.h>
17 #include <TH1.h>
18 #include <TH2.h>
19 #include <TNtuple.h>
20 
21 #include <gsl/gsl_const.h>
22 #include <gsl/gsl_randist.h>
23 
24 #include <boost/foreach.hpp>
25 
26 #include<sstream>
27 
28 using namespace std;
29 
30 // speed of light in cm/ns
31 const double C_light = GSL_CONST_MKS_SPEED_OF_LIGHT / 1e7;
32 const double MASSNEUTRON = 0.939565330;
33 const double MASSPION = 0.13957018;
34 
35 SigmaTimingNtuple::SigmaTimingNtuple(const std::string &name, const std::string &filename):
36  SubsysReco( name ),
37  nblocks(0),
38  hm(NULL),
39  _filename(filename),
40  outfile(NULL)
41 {
42  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
43  unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
44  cout << Name() << " random seed: " << seed << endl;
45  gsl_rng_set(RandomGenerator,seed);
46  return;
47 }
48 
49 SigmaTimingNtuple::~SigmaTimingNtuple()
50 {
51  // delete ntup;
52  delete hm;
53 }
54 
55 
56 int
57 SigmaTimingNtuple::Init( PHCompositeNode* )
58 {
59  ostringstream hname, htit;
60  hm = new Fun4AllHistoManager(Name());
61  outfile = new TFile(_filename.c_str(), "RECREATE");
62 
63  ntupprim = new TNtuple("prim", "Primaries", "pid:phi:theta:eta:e:p");
64  ntupsec = new TNtuple("sec", "Secondaries", "pid:phi:theta:eta:e:p");
65  ntupt = new TNtuple("tnt", "G4Hits", "ID:t:dtotal");
66  ntupsigma = new TNtuple("sigma","Sigmas","mass:inv1:inv2:inv3:t:tres");
67 // ntupavet = new TNtuple("time", "G4Hits", "ID:avt");
68  // ntup->SetDirectory(0);
69  TH1 *h1 = new TH1F("edep1GeV","edep 0-1GeV",1000,0,1);
70  eloss.push_back(h1);
71  h1 = new TH1F("edep100GeV","edep 0-100GeV",1000,0,100);
72  eloss.push_back(h1);
73  return 0;
74 }
75 
76 int
77 SigmaTimingNtuple::process_event( PHCompositeNode* topNode )
78 {
79  ostringstream nodename;
80 
81  map<int, PHG4Particle*>::const_iterator particle_iter;
82  PHG4TruthInfoContainer *_truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
83  if (_truth_container->size() > 3)
84  {
85  return Fun4AllReturnCodes::ABORTEVENT;
86  }
87  double sigmass = 0;
88 
89  PHG4TruthInfoContainer::ConstRange primary_range =
90  _truth_container->GetPrimaryParticleRange();
91  float ntvars[6] = {0};
92  for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
93  particle_iter != primary_range.second; ++particle_iter)
94  {
95  PHG4Particle *particle = particle_iter->second;
96 
97  ntvars[0] = particle->get_pid();
98  ntvars[1] = atan2(particle->get_py(), particle->get_px());
99  ntvars[2] = atan(sqrt(particle->get_py()*particle->get_py() +
100  particle->get_px()*particle->get_px()) /
101  particle->get_pz());
102  if (ntvars[2] < 0)
103  {
104  ntvars[2]+=M_PI;
105  }
106  ntvars[3] = 0.5*log((particle->get_e()+particle->get_pz())/
107  (particle->get_e()-particle->get_pz()));
108  ntvars[4] = particle->get_e();
109  ntvars[5] = sqrt(particle->get_py()*particle->get_py() +
110  particle->get_px()*particle->get_px() +
111  particle->get_pz()*particle->get_pz());
112  sigmass = sqrt(ntvars[4]*ntvars[4] - ntvars[5]*ntvars[5]);
113  // cout << " primary particle " << particle->get_name()
114  // << ", pid: " << particle->get_pid()
115  // << ", mass: " << sigmass << endl;
116  ntupprim->Fill(ntvars);
117  }
118  bool anti_neutron_ok = false;
119  bool piminus_ok = false;
120  int anti_neutron_trk_id = 0;
121 // double piminus_mom = NAN;
122  double anti_neutron_mom = NAN;
123 
124  double E_an = NAN;
125  double px_an = NAN;
126  double py_an = NAN;
127  double pz_an = NAN;
128 
129  double E_pi = NAN;
130  double px_pi = NAN;
131  double py_pi = NAN;
132  double pz_pi = NAN;
133 
134  double p_scale = NAN;
135  double e_an_tof = NAN;
136  double e_an_tof_tres = NAN;
137  double p_scale_tres = NAN;
138 
139  double t_ave = NAN;
140  double t_res = NAN;
141 
142  PHG4TruthInfoContainer::ConstRange secondary_range =
143  _truth_container->GetSecondaryParticleRange();
144  for (PHG4TruthInfoContainer::ConstIterator particle_iter = secondary_range.first;
145  particle_iter != secondary_range.second; ++particle_iter)
146  {
147  PHG4Particle *particle = particle_iter->second;
148  ntvars[0] = particle->get_pid();
149  ntvars[1] = atan2(particle->get_py(), particle->get_px());
150  ntvars[2] = atan(sqrt(particle->get_py()*particle->get_py() +
151  particle->get_px()*particle->get_px()) /
152  particle->get_pz());
153  if (ntvars[2] < 0)
154  {
155  ntvars[2]+=M_PI;
156  }
157  ntvars[3] = 0.5*log((particle->get_e()+particle->get_pz())/
158  (particle->get_e()-particle->get_pz()));
159  ntvars[4] = particle->get_e();
160  ntvars[5] = sqrt(particle->get_py()*particle->get_py() +
161  particle->get_px()*particle->get_px() +
162  particle->get_pz()*particle->get_pz());
163  ntupsec->Fill(ntvars);
164  if (particle->get_pid() == -211)
165  {
166  piminus_ok = true;
167 // piminus_mom = ntvars[5];
168  E_pi = ntvars[4];
169  px_pi = particle->get_px();
170  py_pi = particle->get_py();
171  pz_pi = particle->get_pz();
172  }
173  else if (particle->get_pid() == -2112)
174  {
175  anti_neutron_ok = true;
176  anti_neutron_trk_id = particle->get_track_id();
177  anti_neutron_mom = ntvars[5];
178  E_an = ntvars[4];
179  px_an = particle->get_px();
180  py_an = particle->get_py();
181  pz_an = particle->get_pz();
182  }
183  // cout << " secondary particle " << particle->get_name()
184  // << ", pid: " << particle->get_pid()
185  // << ", G4 track id: " << particle->get_track_id() << endl;
186  }
187  if (! (piminus_ok && anti_neutron_ok))
188  {
189  return Fun4AllReturnCodes::ABORTEVENT;
190  }
191  set<string>::const_iterator iter;
192  vector<TH1 *>::const_iterator eiter;
193 
194  PHG4HitContainer *hits = nullptr;
195 
196  hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_CALO");
197 
198  if (hits)
199  {
200  // double numhits = hits->size();
201  // nhits[i]->Fill(numhits);
202  PHG4HitContainer::ConstRange hit_range = hits->getHits();
203  for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter != hit_range.second; hit_iter++ )
204 
205  {
206  if (hit_iter->second->get_trkid() == anti_neutron_trk_id)
207  {
208 // cout << "G4 trkid: " << hit_iter->second->get_trkid() << endl;
209  t_ave = hit_iter->second->get_avg_t();
210  double dtotal = get_dtotal(hit_iter->second, ntvars[0],ntvars[1]);
211  ntupt->Fill(hit_iter->second->get_layer(),t_ave,dtotal);
212  double dist = sqrt(hit_iter->second->get_avg_x()*hit_iter->second->get_avg_x()+
213  hit_iter->second->get_avg_y()*hit_iter->second->get_avg_y()+
214  hit_iter->second->get_avg_z()*hit_iter->second->get_avg_z());
215  double beta = dist/t_ave/C_light;
216  double neutMom = MASSNEUTRON * beta / sqrt(1.0 - beta*beta);
217  e_an_tof = sqrt(neutMom*neutMom+MASSNEUTRON*MASSNEUTRON);
218 // cout << "neut mom: " << neutMom << ", orig: " << anti_neutron_mom << endl;
219  p_scale = neutMom/anti_neutron_mom;
220  t_res = t_ave+gsl_ran_gaussian(RandomGenerator, 1.);
221  double beta_tres = dist/t_res/C_light;
222 double neutMom_tres = MASSNEUTRON * beta_tres / sqrt(1.0 - beta_tres*beta_tres);
223  e_an_tof_tres = sqrt(neutMom_tres*neutMom_tres+MASSNEUTRON*MASSNEUTRON);
224  p_scale_tres = neutMom_tres/anti_neutron_mom;
225 
226  }
227  }
228  double im2 = MASSNEUTRON*MASSNEUTRON +MASSPION*MASSPION +2*(E_an * E_pi - (px_an*px_pi+py_an*py_pi+pz_an*pz_pi));
229  double im2a = MASSNEUTRON*MASSNEUTRON +MASSPION*MASSPION +2*(e_an_tof * E_pi - ((px_an*p_scale)*px_pi+(py_an*p_scale)*py_pi+(pz_an*p_scale)*pz_pi));
230  double im2b = MASSNEUTRON*MASSNEUTRON +MASSPION*MASSPION +2*(e_an_tof_tres * E_pi - ((px_an*p_scale_tres)*px_pi+(py_an*p_scale_tres)*py_pi+(pz_an*p_scale_tres)*pz_pi));
231 // cout << "im2: " << im2 << ", inv mass: " << sqrt(im2) << endl;
232  ntupsigma->Fill(sigmass,sqrt(im2),sqrt(im2a),sqrt(im2b),t_ave,t_res);
233  }
234 
235  return 0;
236 }
237 
238 int
239 SigmaTimingNtuple::End(PHCompositeNode * topNode)
240 {
241  outfile->cd();
242  // ntup->Write();
243  ntupt->Write();
244  ntupprim->Write();
245  ntupsec->Write();
246  ntupsigma->Write();
247  outfile->Write();
248  outfile->Close();
249  delete outfile;
250  hm->dumpHistos(_filename, "UPDATE");
251  return 0;
252 }
253 
254 void
255 SigmaTimingNtuple::AddNode(const std::string &name, const int detid)
256 {
257  _node_postfix.insert(name);
258  _detid[name] = detid;
259  return;
260 }
261 
262 double
263 SigmaTimingNtuple::get_dtotal(const PHG4Hit *hit, const double phi, const double theta)
264 {
265  double phi_hit = atan2(hit->get_avg_y(),hit->get_avg_x());
266  double theta_hit = atan(sqrt(hit->get_avg_x()*hit->get_avg_x()+
267  hit->get_avg_y()*hit->get_avg_y()) /
268  hit->get_avg_z());
269  double diffphi = phi_hit - phi;
270  if (diffphi > M_PI)
271  {
272  diffphi -= 2*M_PI;
273  }
274  else if (diffphi < - M_PI)
275  {
276  diffphi += 2*M_PI;
277  }
278  if (theta_hit < 0)
279  {
280  theta_hit += M_PI;
281  }
282  double difftheta = theta_hit-theta;
283  double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
284  return deltasqrt;
285 }