AntiSigma/src/SigmaTimingNtuple.cc

0001 #include "SigmaTimingNtuple.h"
0002
0003 #include <g4main/PHG4HitContainer.h>
0004 #include <g4main/PHG4Hit.h>
0005
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4VtxPoint.h>
0009
0010 #include <fun4all/Fun4AllHistoManager.h>
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012
0013 #include <phool/getClass.h>
0014 #include <phool/PHRandomSeed.h>
0015
0016 #include <TFile.h>
0017 #include <TH1.h>
0018 #include <TH2.h>
0019 #include <TNtuple.h>
0020
0021 #include <gsl/gsl_const.h>
0022 #include <gsl/gsl_randist.h>
0023
0024 #include <boost/foreach.hpp>
0025
0026 #include<sstream>
0027
0028 using namespace std;
0029
0030 // speed of light in cm/ns
0031 const double C_light = GSL_CONST_MKS_SPEED_OF_LIGHT / 1e7;
0032 const double MASSNEUTRON = 0.939565330;
0033 const double MASSPION = 0.13957018;
0034
0035 SigmaTimingNtuple::SigmaTimingNtuple(const std::string &name, const std::string &filename):
0036   SubsysReco( name ),
0037   nblocks(0),
0038   hm(NULL),
0039   _filename(filename),
0040   outfile(NULL)
0041 {
0042   RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
0043   unsigned int seed = PHRandomSeed(); // fixed seed is handled in this funtcion
0044   cout << Name() << " random seed: " << seed << endl;
0045   gsl_rng_set(RandomGenerator,seed);
0046   return;
0047 }
0048
0049 SigmaTimingNtuple::~SigmaTimingNtuple()
0050 {
0051   //  delete ntup;
0052   delete hm;
0053 }
0054
0055
0056 int
0057 SigmaTimingNtuple::Init( PHCompositeNode* )
0058 {
0059   ostringstream hname, htit;
0060   hm = new  Fun4AllHistoManager(Name());
0061   outfile = new TFile(_filename.c_str(), "RECREATE");
0062
0063   ntupprim = new TNtuple("prim", "Primaries", "pid:phi:theta:eta:e:p");
0064   ntupsec = new TNtuple("sec", "Secondaries", "pid:phi:theta:eta:e:p");
0065   ntupt = new TNtuple("tnt", "G4Hits", "ID:t:dtotal");
0066   ntupsigma = new TNtuple("sigma","Sigmas","mass:inv1:inv2:inv3:t:tres");
0067 //  ntupavet = new TNtuple("time", "G4Hits", "ID:avt");
0068   //  ntup->SetDirectory(0);
0069   TH1 *h1 = new TH1F("edep1GeV","edep 0-1GeV",1000,0,1);
0070   eloss.push_back(h1);
0071   h1 = new TH1F("edep100GeV","edep 0-100GeV",1000,0,100);
0072   eloss.push_back(h1);
0073   return 0;
0074 }
0075
0076 int
0077 SigmaTimingNtuple::process_event( PHCompositeNode* topNode )
0078 {
0079   ostringstream nodename;
0080
0081   map<int, PHG4Particle*>::const_iterator particle_iter;
0082   PHG4TruthInfoContainer *_truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0083   if (_truth_container->size() > 3)
0084   {
0085     return Fun4AllReturnCodes::ABORTEVENT;
0086   }
0087   double sigmass = 0;
0088
0089   PHG4TruthInfoContainer::ConstRange primary_range =
0090     _truth_container->GetPrimaryParticleRange();
0091   float ntvars[6] = {0};
0092   for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
0093        particle_iter != primary_range.second; ++particle_iter)
0094   {
0095     PHG4Particle *particle = particle_iter->second;
0096
0097     ntvars[0] = particle->get_pid();
0098     ntvars[1] = atan2(particle->get_py(), particle->get_px());
0099     ntvars[2] = atan(sqrt(particle->get_py()*particle->get_py() +
0100               particle->get_px()*particle->get_px()) /
0101              particle->get_pz());
0102     if (ntvars[2] < 0)
0103     {
0104       ntvars[2]+=M_PI;
0105     }
0106     ntvars[3] = 0.5*log((particle->get_e()+particle->get_pz())/
0107             (particle->get_e()-particle->get_pz()));
0108     ntvars[4] = particle->get_e();
0109     ntvars[5] = sqrt(particle->get_py()*particle->get_py() +
0110              particle->get_px()*particle->get_px() +
0111              particle->get_pz()*particle->get_pz());
0112     sigmass = sqrt(ntvars[4]*ntvars[4] - ntvars[5]*ntvars[5]);
0113     // cout << " primary particle " << particle->get_name()
0114     //   << ", pid: " << particle->get_pid()
0115     //   << ", mass: " << sigmass << endl;
0116     ntupprim->Fill(ntvars);
0117   }
0118   bool anti_neutron_ok = false;
0119   bool piminus_ok = false;
0120   int anti_neutron_trk_id = 0;
0121 //  double piminus_mom = NAN;
0122   double anti_neutron_mom = NAN;
0123
0124   double E_an = NAN;
0125   double px_an = NAN;
0126   double py_an = NAN;
0127   double pz_an = NAN;
0128
0129   double E_pi = NAN;
0130   double px_pi = NAN;
0131   double py_pi = NAN;
0132   double pz_pi = NAN;
0133
0134   double p_scale = NAN;
0135   double e_an_tof = NAN;
0136   double e_an_tof_tres = NAN;
0137   double p_scale_tres = NAN;
0138
0139   double t_ave = NAN;
0140   double t_res = NAN;
0141
0142   PHG4TruthInfoContainer::ConstRange secondary_range =
0143     _truth_container->GetSecondaryParticleRange();
0144   for (PHG4TruthInfoContainer::ConstIterator particle_iter = secondary_range.first;
0145        particle_iter != secondary_range.second; ++particle_iter)
0146   {
0147     PHG4Particle *particle = particle_iter->second;
0148     ntvars[0] = particle->get_pid();
0149     ntvars[1] = atan2(particle->get_py(), particle->get_px());
0150     ntvars[2] = atan(sqrt(particle->get_py()*particle->get_py() +
0151               particle->get_px()*particle->get_px()) /
0152              particle->get_pz());
0153     if (ntvars[2] < 0)
0154     {
0155       ntvars[2]+=M_PI;
0156     }
0157     ntvars[3] = 0.5*log((particle->get_e()+particle->get_pz())/
0158             (particle->get_e()-particle->get_pz()));
0159     ntvars[4] = particle->get_e();
0160     ntvars[5] = sqrt(particle->get_py()*particle->get_py() +
0161              particle->get_px()*particle->get_px() +
0162              particle->get_pz()*particle->get_pz());
0163     ntupsec->Fill(ntvars);
0164     if (particle->get_pid() == -211)
0165     {
0166       piminus_ok = true;
0167 //      piminus_mom = ntvars[5];
0168       E_pi = ntvars[4];
0169       px_pi = particle->get_px();
0170       py_pi = particle->get_py();
0171       pz_pi = particle->get_pz();
0172     }
0173     else if (particle->get_pid() == -2112)
0174     {
0175       anti_neutron_ok = true;
0176       anti_neutron_trk_id = particle->get_track_id();
0177       anti_neutron_mom = ntvars[5];
0178       E_an = ntvars[4];
0179       px_an = particle->get_px();
0180       py_an = particle->get_py();
0181       pz_an = particle->get_pz();
0182     }
0183     // cout << " secondary particle " << particle->get_name()
0184     //   << ", pid: " << particle->get_pid()
0185     //   << ", G4 track id: " << particle->get_track_id() << endl;
0186   }
0187   if (! (piminus_ok && anti_neutron_ok))
0188   {
0189     return Fun4AllReturnCodes::ABORTEVENT;
0190   }
0191   set<string>::const_iterator iter;
0192   vector<TH1 *>::const_iterator eiter;
0193
0194   PHG4HitContainer *hits = nullptr;
0195
0196   hits = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_CALO");
0197
0198   if (hits)
0199   {
0200     //          double numhits = hits->size();
0201     //          nhits[i]->Fill(numhits);
0202     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0203     for ( PHG4HitContainer::ConstIterator hit_iter = hit_range.first ; hit_iter !=  hit_range.second; hit_iter++ )
0204
0205     {
0206       if (hit_iter->second->get_trkid() == anti_neutron_trk_id)
0207       {
0208 //  cout << "G4 trkid: " << hit_iter->second->get_trkid() << endl;
0209     t_ave = hit_iter->second->get_avg_t();
0210     double dtotal = get_dtotal(hit_iter->second, ntvars[0],ntvars[1]);
0211     ntupt->Fill(hit_iter->second->get_layer(),t_ave,dtotal);
0212     double dist = sqrt(hit_iter->second->get_avg_x()*hit_iter->second->get_avg_x()+
0213                hit_iter->second->get_avg_y()*hit_iter->second->get_avg_y()+
0214                hit_iter->second->get_avg_z()*hit_iter->second->get_avg_z());
0215     double beta = dist/t_ave/C_light;
0216     double neutMom = MASSNEUTRON * beta / sqrt(1.0 - beta*beta);
0217     e_an_tof = sqrt(neutMom*neutMom+MASSNEUTRON*MASSNEUTRON);
0218 //  cout << "neut mom: " << neutMom << ", orig: " << anti_neutron_mom << endl;
0219     p_scale = neutMom/anti_neutron_mom;
0220     t_res = t_ave+gsl_ran_gaussian(RandomGenerator, 1.);
0221     double beta_tres = dist/t_res/C_light;
0222 double neutMom_tres = MASSNEUTRON * beta_tres / sqrt(1.0 - beta_tres*beta_tres);
0223     e_an_tof_tres = sqrt(neutMom_tres*neutMom_tres+MASSNEUTRON*MASSNEUTRON);
0224     p_scale_tres = neutMom_tres/anti_neutron_mom;
0225
0226       }
0227     }
0228     double im2 = MASSNEUTRON*MASSNEUTRON +MASSPION*MASSPION +2*(E_an * E_pi - (px_an*px_pi+py_an*py_pi+pz_an*pz_pi));
0229     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));
0230     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));
0231 //    cout << "im2: " << im2 << ", inv mass: " << sqrt(im2) << endl;
0232     ntupsigma->Fill(sigmass,sqrt(im2),sqrt(im2a),sqrt(im2b),t_ave,t_res);
0233   }
0234
0235   return 0;
0236 }
0237
0238 int
0239 SigmaTimingNtuple::End(PHCompositeNode * topNode)
0240 {
0241   outfile->cd();
0242   //  ntup->Write();
0243   ntupt->Write();
0244   ntupprim->Write();
0245   ntupsec->Write();
0246   ntupsigma->Write();
0247   outfile->Write();
0248   outfile->Close();
0249   delete outfile;
0250   hm->dumpHistos(_filename, "UPDATE");
0251   return 0;
0252 }
0253
0254 void
0255 SigmaTimingNtuple::AddNode(const std::string &name, const int detid)
0256 {
0257   _node_postfix.insert(name);
0258   _detid[name] = detid;
0259   return;
0260 }
0261
0262 double
0263 SigmaTimingNtuple::get_dtotal(const PHG4Hit *hit, const double phi, const double theta)
0264 {
0265   double phi_hit =  atan2(hit->get_avg_y(),hit->get_avg_x());
0266   double theta_hit = atan(sqrt(hit->get_avg_x()*hit->get_avg_x()+
0267                    hit->get_avg_y()*hit->get_avg_y()) /
0268               hit->get_avg_z());
0269   double diffphi =  phi_hit - phi;
0270   if (diffphi > M_PI)
0271   {
0272     diffphi -= 2*M_PI;
0273   }
0274   else if (diffphi < - M_PI)
0275   {
0276     diffphi += 2*M_PI;
0277   }
0278   if (theta_hit < 0)
0279   {
0280     theta_hit += M_PI;
0281   }
0282   double difftheta = theta_hit-theta;
0283   double deltasqrt = sqrt(diffphi*diffphi+difftheta*difftheta);
0284   return deltasqrt;
0285 }