Public Member Functions | Private Member Functions | Private Attributes | List of all members
gar::ParamSim Class Reference
Inheritance diagram for gar::ParamSim:
art::EDAnalyzer art::detail::Analyzer art::detail::LegacyModule art::Observer art::ModuleBase

Public Member Functions

 ParamSim (fhicl::ParameterSet const &p)
 
 ParamSim (ParamSim const &)=delete
 
 ParamSim (ParamSim &&)=delete
 
ParamSimoperator= (ParamSim const &)=delete
 
ParamSimoperator= (ParamSim &&)=delete
 
virtual void beginJob () override
 
void analyze (art::Event const &e) override
 
- Public Member Functions inherited from art::EDAnalyzer
 EDAnalyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDAnalyzer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Analyzer
virtual ~Analyzer () noexcept
 
 Analyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 Analyzer (Table< Config > const &config)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Observer
 ~Observer () noexcept
 
 Observer (Observer const &)=delete
 
 Observer (Observer &&)=delete
 
Observeroperator= (Observer const &)=delete
 
Observeroperator= (Observer &&)=delete
 
void registerProducts (ProductDescriptions &, ModuleDescription const &)
 
void fillDescriptions (ModuleDescription const &)
 
fhicl::ParameterSetID selectorConfig () const
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

void ClearVectors ()
 
void SaveGtruthMCtruth (art::Event const &e)
 
void TreatMCParticles (art::Event const &e)
 
void TreatTPCVisible (const float ecaltime)
 
void TreatTPCNotVisible (const float ecaltime)
 
void FillCommonVariables (const int pdg, const TLorentzVector &momentum, const TLorentzVector &position, const TLorentzVector &positionEnd, const int mctrackid, const int mothertrackid, const int motherpdg, const float ptrue, const float angle, const std::string mcp_process, const std::string mcp_endprocess, const float time)
 
void ComputeTrkLength (const simb::MCParticle &mcp)
 
void DoRangeCalculation (float &preco, float &angle_reco)
 
void DoGluckSternCalculation (float &preco, float &angle_reco)
 
void TPCParticleIdentification ()
 
void TreatNeutrons (const float ecaltime)
 
void TreatPhotons (const float ecaltime)
 
void TreatOthers (const float ecaltime)
 
bool CheckVectorSize ()
 

Private Attributes

TTree * fTree
 
std::unordered_map< int, TH2F * > m_pidinterp
 
const geo::GeometryCorefGeo
 pointer to the geometry More...
 
CAFHelperfHelper
 
CLHEP::HepRandomEngine & fEngine
 random engine More...
 
std::string fGeneratorLabel
 
std::string fGeantLabel
 module label for geant4 simulated hits More...
 
bool fCorrect4origin
 
float fOrigin [3]
 
const std::vector< int > neutrinos = {12, 14, 16}
 
const std::vector< int > pdg_neutral = {22, 2112, 111, 130, 310, 311, 2114}
 
const std::vector< int > pdg_charged = {211, 13, 2212, 321, 1000010020, 11}
 
const float gastpc_len = 2.
 
const float gastpc_B = 0.5
 
const float gastpc_padPitch = 0.1
 
const float gastpc_X0 = 1300.
 
const float sigmaP_short = 0.1
 
const float sigma_x = 0.1
 
const float NeutronECAL_detEff [2] = {0.2, 0.4}
 
const float sigmaNeutronECAL_first = 0.11
 
const float ECAL_stock = 0.06
 
const float ECAL_const = 0.02
 
TF1 * fRes
 
const int nLayers = 60
 
const double ECAL_MIP_Res = 0.23
 
const double MIP2GeV_factor = 0.814 / 1000
 
const float fECALTimeResolution = 1.
 
TParticlePDG * neutron = TDatabasePDG::Instance()->GetParticle(2112)
 
const float neutron_mass = neutron->Mass()
 
unsigned int fRun
 
unsigned int fEvent
 
unsigned int fSubRun
 
std::vector< int > mode
 
std::vector< int > ccnc
 
std::vector< int > ntype
 
std::vector< int > gint
 
std::vector< int > weight
 
std::vector< int > tgtpdg
 
std::vector< int > gt_t
 
std::vector< int > intert
 
std::vector< int > detected
 
std::vector< double > q2
 
std::vector< double > w
 
std::vector< double > y
 
std::vector< double > x
 
std::vector< double > theta
 
std::vector< double > t
 
std::vector< double > mctime
 
std::vector< double > mcnupx
 
std::vector< double > mcnupy
 
std::vector< double > mcnupz
 
std::vector< double > vertx
 
std::vector< double > verty
 
std::vector< double > vertz
 
unsigned int _nFSP
 
std::vector< int > mctrkid
 
std::vector< int > motherid
 
std::vector< int > pdgmother
 
std::vector< int > truepdg
 
std::vector< int > _MCPStartX
 
std::vector< int > _MCPStartY
 
std::vector< int > _MCPStartZ
 
std::vector< int > _MCPEndX
 
std::vector< int > _MCPEndY
 
std::vector< int > _MCPEndZ
 
std::vector< std::string_MCProc
 
std::vector< std::string_MCEndProc
 
std::vector< double > trkLen
 
std::vector< double > trkLenPerp
 
std::vector< double > truep
 
std::vector< double > truepx
 
std::vector< double > truepy
 
std::vector< double > truepz
 
std::vector< double > _angle
 
std::vector< int > recopid
 
std::vector< int > recopidecal
 
std::vector< double > prob_arr
 
std::vector< double > anglereco
 
std::vector< double > _preco
 
std::vector< double > erecon
 
std::vector< double > etime
 
std::vector< unsigned int > isFidStart
 
std::vector< unsigned int > isTPCStart
 
std::vector< unsigned int > isCaloStart
 
std::vector< unsigned int > isInBetweenStart
 
std::vector< unsigned int > isThroughCaloStart
 
std::vector< unsigned int > isFidEnd
 
std::vector< unsigned int > isTPCEnd
 
std::vector< unsigned int > isCaloEnd
 
std::vector< unsigned int > isInBetweenEnd
 
std::vector< unsigned int > isThroughCaloEnd
 
std::vector< unsigned int > isBarrelStart
 
std::vector< unsigned int > isEndcapStart
 
std::vector< unsigned int > isBarrelEnd
 
std::vector< unsigned int > isEndcapEnd
 

Additional Inherited Members

- Public Types inherited from art::EDAnalyzer
using WorkerType = WorkerT< EDAnalyzer >
 
using ModuleType = EDAnalyzer
 
- Protected Member Functions inherited from art::Observer
std::string const & processName () const
 
bool wantAllEvents () const noexcept
 
bool wantEvent (ScheduleID id, Event const &e) const
 
Handle< TriggerResultsgetTriggerResults (Event const &e) const
 
 Observer (fhicl::ParameterSet const &config)
 
 Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 250 of file ParamSim_module.cc.

Constructor & Destructor Documentation

gar::ParamSim::ParamSim ( fhicl::ParameterSet const &  p)
explicit

Definition at line 359 of file ParamSim_module.cc.

360  : EDAnalyzer(p),
362  {
363  fGeo = gar::providerFrom<geo::GeometryGAr>();
364 
365  fGeneratorLabel = p.get<std::string>("GeneratorLabel","genie");
366  fGeantLabel = p.get<std::string>("GEANTLabel","geant");
367  fCorrect4origin = p.get<bool>("Correct4Origin", false);
368 
369  consumes<std::vector<simb::MCTruth> >(fGeneratorLabel);
370  consumes<std::vector<simb::GTruth> >(fGeneratorLabel);
371  consumes<std::vector<simb::MCParticle> >(fGeantLabel);
372  }
base_engine_t & createEngine(seed_t seed)
CLHEP::HepRandomEngine & fEngine
random engine
std::string string
Definition: nybbler.cc:12
std::string fGeneratorLabel
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
std::string fGeantLabel
module label for geant4 simulated hits
p
Definition: test.py:223
const geo::GeometryCore * fGeo
pointer to the geometry
gar::ParamSim::ParamSim ( ParamSim const &  )
delete
gar::ParamSim::ParamSim ( ParamSim &&  )
delete

Member Function Documentation

void gar::ParamSim::analyze ( art::Event const &  e)
overridevirtual

Implements art::EDAnalyzer.

Definition at line 485 of file ParamSim_module.cc.

486  {
487  ClearVectors();
488  fRun = e.run();
489  fSubRun = e.subRun();
490  fEvent = e.id().event();
491 
494 
495  //Checks
496  if(CheckVectorSize()) {
497  fTree->Fill();
498  } else {
499  std::cerr << "Event " << fEvent << std::endl;
500  std::cerr << "Number of FSP " << _nFSP << std::endl;
501 
502  std::cerr << "Size of pdgmother " << pdgmother.size() << std::endl;
503  std::cerr << "Size of truepdg " << truepdg.size() << std::endl;
504  std::cerr << "Size of mctime " << mctime.size() << std::endl;
505  std::cerr << "Size of mctrkid " << mctrkid.size() << std::endl;
506  std::cerr << "Size of motherid " << motherid.size() << std::endl;
507  std::cerr << "Size of _MCPStartX " << _MCPStartX.size() << std::endl;
508  std::cerr << "Size of _MCPStartY " << _MCPStartY.size() << std::endl;
509  std::cerr << "Size of _MCPStartZ " << _MCPStartZ.size() << std::endl;
510  std::cerr << "Size of _MCPEndX " << _MCPEndX.size() << std::endl;
511  std::cerr << "Size of _MCPEndY " << _MCPEndY.size() << std::endl;
512  std::cerr << "Size of _MCPEndZ " << _MCPEndZ.size() << std::endl;
513  std::cerr << "Size of _MCProc " << _MCProc.size() << std::endl;
514  std::cerr << "Size of _MCEndProc " << _MCEndProc.size() << std::endl;
515  std::cerr << "Size of trkLen " << trkLen.size() << std::endl;
516  std::cerr << "Size of trkLenPerp " << trkLenPerp.size() << std::endl;
517  std::cerr << "Size of truep " << truep.size() << std::endl;
518  std::cerr << "Size of truepx " << truepx.size() << std::endl;
519  std::cerr << "Size of truepy " << truepy.size() << std::endl;
520  std::cerr << "Size of truepz " << truepz.size() << std::endl;
521  std::cerr << "Size of _angle " << _angle.size() << std::endl;
522 
523  //Reco values
524  std::cerr << "Size of detected " << detected.size() << std::endl;
525  std::cerr << "Size of recopid " << recopid.size() << std::endl;
526  std::cerr << "Size of recopidecal " << recopidecal.size() << std::endl;
527  // std::cerr << "Size of prob_arr " << prob_arr.size() << std::endl;
528  std::cerr << "Size of anglereco " << anglereco.size() << std::endl;
529  std::cerr << "Size of _preco " << _preco.size() << std::endl;
530  std::cerr << "Size of erecon " << erecon.size() << std::endl;
531  std::cerr << "Size of etime " << etime.size() << std::endl;
532 
533  //Geometry
534  std::cerr << "Size of isFidStart " << isFidStart.size() << std::endl;
535  std::cerr << "Size of isTPCStart " << isTPCStart.size() << std::endl;
536  std::cerr << "Size of isCaloStart " << isCaloStart.size() << std::endl;
537  std::cerr << "Size of isThroughCaloStart " << isThroughCaloStart.size() << std::endl;
538  std::cerr << "Size of isInBetweenStart " << isInBetweenStart.size() << std::endl;
539  std::cerr << "Size of isBarrelStart " << isBarrelStart.size() << std::endl;
540  std::cerr << "Size of isEndcapStart " << isEndcapStart.size() << std::endl;
541 
542  std::cerr << "Size of isFidEnd " << isFidEnd.size() << std::endl;
543  std::cerr << "Size of isTPCEnd " << isTPCEnd.size() << std::endl;
544  std::cerr << "Size of isCaloEnd " << isCaloEnd.size() << std::endl;
545  std::cerr << "Size of isThroughCaloEnd " << isThroughCaloEnd.size() << std::endl;
546  std::cerr << "Size of isInBetweenEnd " << isInBetweenEnd.size() << std::endl;
547  std::cerr << "Size of isBarrelEnd " << isBarrelEnd.size() << std::endl;
548  std::cerr << "Size of isEndcapEnd " << isEndcapEnd.size() << std::endl;
549 
550  std::cerr << "Event with wrong vector sizes... skipped" << std::endl;
551  }
552 
553  }
std::vector< double > _preco
std::vector< double > truepz
unsigned int fRun
std::vector< unsigned int > isFidEnd
std::vector< std::string > _MCEndProc
std::vector< unsigned int > isBarrelStart
std::vector< int > detected
unsigned int fEvent
std::vector< std::string > _MCProc
std::vector< int > recopidecal
std::vector< double > _angle
std::vector< double > truepx
std::vector< unsigned int > isTPCStart
std::vector< int > motherid
std::vector< int > _MCPEndY
std::vector< unsigned int > isEndcapEnd
std::vector< int > truepdg
const double e
std::vector< unsigned int > isTPCEnd
std::vector< int > mctrkid
std::vector< unsigned int > isFidStart
std::vector< int > _MCPEndZ
std::vector< unsigned int > isThroughCaloEnd
std::vector< double > truepy
std::vector< double > trkLen
std::vector< unsigned int > isEndcapStart
void TreatMCParticles(art::Event const &e)
std::vector< int > recopid
void SaveGtruthMCtruth(art::Event const &e)
std::vector< unsigned int > isCaloStart
std::vector< unsigned int > isThroughCaloStart
std::vector< unsigned int > isCaloEnd
unsigned int fSubRun
std::vector< int > _MCPStartZ
std::vector< double > truep
std::vector< unsigned int > isBarrelEnd
std::vector< int > _MCPStartY
unsigned int _nFSP
std::vector< double > mctime
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< unsigned int > isInBetweenEnd
std::vector< unsigned int > isInBetweenStart
std::vector< int > _MCPStartX
std::vector< int > pdgmother
QTextStream & endl(QTextStream &s)
std::vector< double > trkLenPerp
std::vector< double > etime
void gar::ParamSim::beginJob ( )
overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 375 of file ParamSim_module.cc.

376  {
377  fOrigin[0] = fGeo->GetOriginX();
378  fOrigin[1] = fGeo->GetOriginY();
379  fOrigin[2] = fGeo->GetOriginZ();
380 
381  fRes = new TF1("fRes", "TMath::Sqrt ( [0]*[0]/x + [1]*[1] )", 3);
382  fRes->FixParameter(0, ECAL_stock);
383  fRes->FixParameter(1, ECAL_const);
384  fHelper = new CAFHelper(fGeo, fEngine);
385 
386  // read the PID parametrization ntuple from T. Junk
387  TString filename = "${DUNE_PARDATA_DIR}/MPD/dedxPID/dedxpidmatrices8kevcm.root";
388  TFile pidfile(filename, "READ");
389 
390  m_pidinterp.clear();
391  char str[11];
392  for (int q = 0; q < 501; ++q)
393  {
394  sprintf(str, "%d", q);
395  std::string s = "pidmatrix";
396  s.append(str);
397  // read the 500 histograms one by one; each histogram is a
398  // 6 by 6 matrix of probabilities for a given momentum value
399  m_pidinterp.insert( std::make_pair(q, (TH2F*) pidfile.Get(s.c_str())->Clone("pidinterp")) );
400  }
401 
402  // pidfile.Close();
403 
405  fTree = tfs->make<TTree>("caf","caf tree");
406 
407  //Event number, //Run number, //Sub Run number
408  fTree->Branch("Event", &fEvent);
409  fTree->Branch("Run", &fRun);
410  fTree->Branch("SubRun", &fSubRun);
411 
412  //MC Truth info (Generator)
413  fTree->Branch("mode", &mode);
414  fTree->Branch("q2", &q2);
415  fTree->Branch("w", &w);
416  fTree->Branch("y", &y);
417  fTree->Branch("x", &x);
418  fTree->Branch("theta", &theta);
419  fTree->Branch("t", &t);
420  fTree->Branch("ntype", &ntype);
421  fTree->Branch("ccnc", &ccnc);
422  fTree->Branch("gint", &gint);
423  fTree->Branch("tgtpdg", &tgtpdg);
424  fTree->Branch("weight", &weight);
425  fTree->Branch("gt_t", &gt_t);
426  fTree->Branch("intert", &intert);
427  fTree->Branch("mcnupx", &mcnupx);
428  fTree->Branch("mcnupy", &mcnupy);
429  fTree->Branch("mcnupz", &mcnupz);
430  fTree->Branch("vertx", &vertx);
431  fTree->Branch("verty", &verty);
432  fTree->Branch("vertz", &vertz);
433 
434  //Number of final state particle (primaries)
435  fTree->Branch("nFSP", &_nFSP);
436  //MC Particle info
437  fTree->Branch("detected", &detected);
438  fTree->Branch("pdgmother", &pdgmother);
439  fTree->Branch("MCPTime", &mctime);
440  fTree->Branch("MCPStartX", &_MCPStartX);
441  fTree->Branch("MCPStartY", &_MCPStartY);
442  fTree->Branch("MCPStartZ", &_MCPStartZ);
443  fTree->Branch("motherid", &motherid);
444  fTree->Branch("mctrkid", &mctrkid);
445  fTree->Branch("truepx", &truepx);
446  fTree->Branch("truepy", &truepy);
447  fTree->Branch("truepz", &truepz);
448  fTree->Branch("MCPEndX", &_MCPEndX);
449  fTree->Branch("MCPEndY", &_MCPEndY);
450  fTree->Branch("MCPEndZ", &_MCPEndZ);
451  fTree->Branch("MCProc", &_MCProc);
452  fTree->Branch("MCEndProc", &_MCEndProc);
453  fTree->Branch("angle", &_angle);
454  fTree->Branch("truep", &truep);
455  fTree->Branch("truepdg", &truepdg);
456  //Reco info
457  fTree->Branch("recopid", &recopid);
458  fTree->Branch("recopidecal", &recopidecal);
459  fTree->Branch("trkLen", &trkLen);
460  fTree->Branch("trkLenPerp", &trkLenPerp);
461  fTree->Branch("preco", &_preco);
462  fTree->Branch("anglereco", &anglereco);
463  fTree->Branch("erecon", &erecon);
464  fTree->Branch("etime", &etime);
465  fTree->Branch("prob_arr", &prob_arr);
466  //Geometry
467  fTree->Branch("isFidStart", &isFidStart);
468  fTree->Branch("isTPCStart", &isTPCStart);
469  fTree->Branch("isCaloStart", &isCaloStart);
470  fTree->Branch("isInBetweenStart", &isInBetweenStart);
471  fTree->Branch("isThroughCaloStart", &isThroughCaloStart);
472  fTree->Branch("isBarrelStart", &isBarrelStart);
473  fTree->Branch("isEndcapStart", &isEndcapStart);
474 
475  fTree->Branch("isFidEnd", &isFidEnd);
476  fTree->Branch("isTPCEnd", &isTPCEnd);
477  fTree->Branch("isCaloEnd", &isCaloEnd);
478  fTree->Branch("isThroughCaloEnd", &isThroughCaloEnd);
479  fTree->Branch("isInBetweenEnd", &isInBetweenEnd);
480  fTree->Branch("isBarrelEnd", &isBarrelEnd);
481  fTree->Branch("isEndcapEnd", &isEndcapEnd);
482  }
std::vector< double > _preco
std::vector< int > weight
std::vector< double > truepz
std::vector< int > ccnc
unsigned int fRun
CLHEP::HepRandomEngine & fEngine
random engine
std::vector< unsigned int > isFidEnd
std::vector< double > y
std::vector< double > mcnupy
std::vector< std::string > _MCEndProc
std::string string
Definition: nybbler.cc:12
std::vector< unsigned int > isBarrelStart
std::vector< int > detected
float GetOriginX() const
Definition: GeometryCore.h:548
std::vector< int > ntype
unsigned int fEvent
std::vector< std::string > _MCProc
std::unordered_map< int, TH2F * > m_pidinterp
std::vector< int > recopidecal
string filename
Definition: train.py:213
std::vector< double > _angle
std::vector< double > truepx
std::vector< int > gint
std::vector< unsigned int > isTPCStart
std::vector< double > vertx
float GetOriginZ() const
Definition: GeometryCore.h:552
std::vector< int > motherid
std::vector< double > mcnupx
std::vector< double > vertz
std::vector< int > _MCPEndY
std::vector< unsigned int > isEndcapEnd
std::vector< int > truepdg
std::vector< double > verty
std::vector< unsigned int > isTPCEnd
const float ECAL_const
std::vector< int > mctrkid
std::vector< unsigned int > isFidStart
std::vector< int > gt_t
std::vector< int > _MCPEndZ
std::vector< double > w
std::vector< unsigned int > isThroughCaloEnd
float GetOriginY() const
Definition: GeometryCore.h:550
std::vector< int > mode
std::vector< double > q2
std::vector< double > truepy
std::vector< double > trkLen
std::vector< unsigned int > isEndcapStart
std::vector< double > prob_arr
std::vector< int > recopid
std::vector< unsigned int > isCaloStart
std::vector< double > x
std::vector< unsigned int > isThroughCaloStart
std::vector< unsigned int > isCaloEnd
unsigned int fSubRun
std::vector< int > _MCPStartZ
std::vector< double > truep
std::vector< unsigned int > isBarrelEnd
std::vector< int > _MCPStartY
unsigned int _nFSP
std::vector< double > mcnupz
CAFHelper * fHelper
std::vector< double > theta
std::vector< double > t
std::vector< double > mctime
const geo::GeometryCore * fGeo
pointer to the geometry
std::vector< int > intert
const float ECAL_stock
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< unsigned int > isInBetweenEnd
std::vector< int > tgtpdg
static QCString * s
Definition: config.cpp:1042
std::vector< unsigned int > isInBetweenStart
std::vector< int > _MCPStartX
static QCString str
std::vector< int > pdgmother
QCString & append(const char *s)
Definition: qcstring.cpp:383
std::vector< double > trkLenPerp
std::vector< double > etime
bool gar::ParamSim::CheckVectorSize ( )
private

Definition at line 1556 of file ParamSim_module.cc.

1557  {
1558  bool isOK = true;
1559 
1560  if(_nFSP != pdgmother.size()) isOK = false;
1561  if(_nFSP != truepdg.size()) isOK = false;
1562  if(_nFSP != mctime.size()) isOK = false;
1563  if(_nFSP != mctrkid.size()) isOK = false;
1564  if(_nFSP != motherid.size()) isOK = false;
1565  if(_nFSP != _MCPStartX.size()) isOK = false;
1566  if(_nFSP != _MCPStartY.size()) isOK = false;
1567  if(_nFSP != _MCPStartZ.size()) isOK = false;
1568  if(_nFSP != _MCPEndX.size()) isOK = false;
1569  if(_nFSP != _MCPEndY.size()) isOK = false;
1570  if(_nFSP != _MCPEndZ.size()) isOK = false;
1571  if(_nFSP != _MCProc.size()) isOK = false;
1572  if(_nFSP != _MCEndProc.size()) isOK = false;
1573  if(_nFSP != trkLen.size()) isOK = false;
1574  if(_nFSP != trkLenPerp.size()) isOK = false;
1575  if(_nFSP != truep.size()) isOK = false;
1576  if(_nFSP != truepx.size()) isOK = false;
1577  if(_nFSP != truepy.size()) isOK = false;
1578  if(_nFSP != truepz.size()) isOK = false;
1579  if(_nFSP != _angle.size()) isOK = false;
1580 
1581  //Reco values
1582  if(_nFSP != recopid.size()) isOK = false;
1583  if(_nFSP != detected.size()) isOK = false;
1584  if(_nFSP != recopidecal.size()) isOK = false;
1585  // if(_nFSP != prob_arr.size()) isOK = false;
1586  if(_nFSP != anglereco.size()) isOK = false;
1587  if(_nFSP != _preco.size()) isOK = false;
1588  if(_nFSP != erecon.size()) isOK = false;
1589  if(_nFSP != etime.size()) isOK = false;
1590 
1591  //Geometry
1592  if(_nFSP != isFidStart.size()) isOK = false;
1593  if(_nFSP != isTPCStart.size()) isOK = false;
1594  if(_nFSP != isCaloStart.size()) isOK = false;
1595  if(_nFSP != isThroughCaloStart.size()) isOK = false;
1596  if(_nFSP != isInBetweenStart.size()) isOK = false;
1597  if(_nFSP != isBarrelStart.size()) isOK = false;
1598  if(_nFSP != isEndcapStart.size()) isOK = false;
1599 
1600  if(_nFSP != isFidEnd.size()) isOK = false;
1601  if(_nFSP != isTPCEnd.size()) isOK = false;
1602  if(_nFSP != isCaloEnd.size()) isOK = false;
1603  if(_nFSP != isThroughCaloEnd.size()) isOK = false;
1604  if(_nFSP != isInBetweenEnd.size()) isOK = false;
1605  if(_nFSP != isBarrelEnd.size()) isOK = false;
1606  if(_nFSP != isEndcapEnd.size()) isOK = false;
1607 
1608  return isOK;
1609  }
std::vector< double > _preco
std::vector< double > truepz
std::vector< unsigned int > isFidEnd
std::vector< std::string > _MCEndProc
std::vector< unsigned int > isBarrelStart
std::vector< int > detected
std::vector< std::string > _MCProc
std::vector< int > recopidecal
std::vector< double > _angle
std::vector< double > truepx
std::vector< unsigned int > isTPCStart
std::vector< int > motherid
std::vector< int > _MCPEndY
std::vector< unsigned int > isEndcapEnd
std::vector< int > truepdg
std::vector< unsigned int > isTPCEnd
std::vector< int > mctrkid
std::vector< unsigned int > isFidStart
std::vector< int > _MCPEndZ
std::vector< unsigned int > isThroughCaloEnd
std::vector< double > truepy
std::vector< double > trkLen
std::vector< unsigned int > isEndcapStart
std::vector< int > recopid
std::vector< unsigned int > isCaloStart
std::vector< unsigned int > isThroughCaloStart
std::vector< unsigned int > isCaloEnd
std::vector< int > _MCPStartZ
std::vector< double > truep
std::vector< unsigned int > isBarrelEnd
std::vector< int > _MCPStartY
unsigned int _nFSP
std::vector< double > mctime
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< unsigned int > isInBetweenEnd
std::vector< unsigned int > isInBetweenStart
std::vector< int > _MCPStartX
std::vector< int > pdgmother
std::vector< double > trkLenPerp
std::vector< double > etime
void gar::ParamSim::ClearVectors ( )
private

Definition at line 1480 of file ParamSim_module.cc.

1481  {
1482  //Generator values
1483  mode.clear();
1484  ccnc.clear();
1485  ntype.clear();
1486  gint.clear();
1487  weight.clear();
1488  tgtpdg.clear();
1489  gt_t.clear();
1490  intert.clear();
1491  q2.clear();
1492  w.clear();
1493  y.clear();
1494  x.clear();
1495  theta.clear();
1496  t.clear();
1497  mcnupx.clear();
1498  mcnupy.clear();
1499  mcnupz.clear();
1500  vertx.clear();
1501  verty.clear();
1502  vertz.clear();
1503 
1504  //MC Particle values
1505  _nFSP = 0;
1506  detected.clear();
1507  pdgmother.clear();
1508  truepdg.clear();
1509  mctime.clear();
1510  mctrkid.clear();
1511  motherid.clear();
1512  _MCPStartX.clear();
1513  _MCPStartY.clear();
1514  _MCPStartZ.clear();
1515  _MCPEndX.clear();
1516  _MCPEndY.clear();
1517  _MCPEndZ.clear();
1518  _MCProc.clear();
1519  _MCEndProc.clear();
1520  trkLen.clear();
1521  trkLenPerp.clear();
1522  truep.clear();
1523  truepx.clear();
1524  truepy.clear();
1525  truepz.clear();
1526  _angle.clear();
1527 
1528  //Reco values
1529  recopid.clear();
1530  recopidecal.clear();
1531  prob_arr.clear();
1532  anglereco.clear();
1533  _preco.clear();
1534  erecon.clear();
1535  etime.clear();
1536 
1537  //Geometry
1538  isFidStart.clear();
1539  isTPCStart.clear();
1540  isCaloStart.clear();
1541  isThroughCaloStart.clear();
1542  isInBetweenStart.clear();
1543  isBarrelStart.clear();
1544  isEndcapStart.clear();
1545 
1546  isFidEnd.clear();
1547  isTPCEnd.clear();
1548  isCaloEnd.clear();
1549  isThroughCaloEnd.clear();
1550  isInBetweenEnd.clear();
1551  isBarrelEnd.clear();
1552  isEndcapEnd.clear();
1553  }
std::vector< double > _preco
std::vector< int > weight
std::vector< double > truepz
std::vector< int > ccnc
std::vector< unsigned int > isFidEnd
std::vector< double > y
std::vector< double > mcnupy
std::vector< std::string > _MCEndProc
std::vector< unsigned int > isBarrelStart
std::vector< int > detected
std::vector< int > ntype
std::vector< std::string > _MCProc
std::vector< int > recopidecal
std::vector< double > _angle
std::vector< double > truepx
std::vector< int > gint
std::vector< unsigned int > isTPCStart
std::vector< double > vertx
std::vector< int > motherid
std::vector< double > mcnupx
std::vector< double > vertz
std::vector< int > _MCPEndY
std::vector< unsigned int > isEndcapEnd
std::vector< int > truepdg
std::vector< double > verty
std::vector< unsigned int > isTPCEnd
std::vector< int > mctrkid
std::vector< unsigned int > isFidStart
std::vector< int > gt_t
std::vector< int > _MCPEndZ
std::vector< double > w
std::vector< unsigned int > isThroughCaloEnd
std::vector< int > mode
std::vector< double > q2
std::vector< double > truepy
std::vector< double > trkLen
std::vector< unsigned int > isEndcapStart
std::vector< double > prob_arr
std::vector< int > recopid
std::vector< unsigned int > isCaloStart
std::vector< double > x
std::vector< unsigned int > isThroughCaloStart
std::vector< unsigned int > isCaloEnd
std::vector< int > _MCPStartZ
std::vector< double > truep
std::vector< unsigned int > isBarrelEnd
std::vector< int > _MCPStartY
unsigned int _nFSP
std::vector< double > mcnupz
std::vector< double > theta
std::vector< double > t
std::vector< double > mctime
std::vector< int > intert
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< unsigned int > isInBetweenEnd
std::vector< int > tgtpdg
std::vector< unsigned int > isInBetweenStart
std::vector< int > _MCPStartX
std::vector< int > pdgmother
std::vector< double > trkLenPerp
std::vector< double > etime
void gar::ParamSim::ComputeTrkLength ( const simb::MCParticle mcp)
private

Definition at line 739 of file ParamSim_module.cc.

740  {
741  //start track length
742  //***************************************************************************************************************/
743 
744  // calculate the total and the transverse track lengths and restrict the
745  // tracklength to be above the gas TPC track length threshold
746  double tracklen = 0.;
747  double tracklen_perp = 0.;
748 
749  //CAREFUL No offset for the trajectory points (origin for them is the TPC?)??????
750  //TODO check if the mcp point is within the TPC volume! Skip for mcp in the ECAL (showers)
751  //TODO Link showers to original mcp?
752  for(size_t itraj = 1; itraj < mcp.Trajectory().size(); itraj++)
753  {
754  float xTraj = mcp.Trajectory().X(itraj);
755  float yTraj = mcp.Trajectory().Y(itraj);
756  float zTraj = mcp.Trajectory().Z(itraj);
757 
758  //Traj point+1
759  TVector3 point(xTraj - fOrigin[0], yTraj - fOrigin[1], zTraj - fOrigin[2]);
760 
761  //point is not in the TPC anymore - stop traj loop
762  if(not fHelper->PointInTPC(point))
763  continue;
764 
765  // find the length of the track by getting the distance between each hit
766  TVector3 diff(xTraj - mcp.Trajectory().X(itraj - 1), yTraj - mcp.Trajectory().Y(itraj - 1), zTraj - mcp.Trajectory().Z(itraj - 1));
767  // perp length
768  TVector2 tracklen_perp_vec(zTraj - mcp.Trajectory().Z(itraj - 1), yTraj - mcp.Trajectory().Y(itraj - 1));
769  // Summing up
770  tracklen += diff.Mag();
771  tracklen_perp += tracklen_perp_vec.Mod();
772  }
773 
774  trkLen.push_back(tracklen);
775  trkLenPerp.push_back(tracklen_perp);
776  }
double Z(const size_type i) const
Definition: MCTrajectory.h:151
double X(const size_type i) const
Definition: MCTrajectory.h:149
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
bool PointInTPC(const TVector3 &point)
double Y(const size_type i) const
Definition: MCTrajectory.h:150
std::vector< double > trkLen
size_type size() const
Definition: MCTrajectory.h:166
CAFHelper * fHelper
std::vector< double > trkLenPerp
void gar::ParamSim::DoGluckSternCalculation ( float &  preco,
float &  angle_reco 
)
private

Definition at line 797 of file ParamSim_module.cc.

798  {
799  //Case where the endpoint is not in the TPC, should be able to use the Gluckstern formula
800  // calculate number of trackpoints
801  float nHits = round (trkLen.at(_nFSP) / gastpc_padPitch);
802  // measurement term in Gluckstern formula
803  float fracSig_meas = sqrt(720./(nHits+4)) * ((0.01*gastpc_padPitch*truep.at(_nFSP)) / (0.3 * gastpc_B * 0.0001 *trkLenPerp.at(_nFSP)*trkLenPerp.at(_nFSP)));
804  // multiple Coulomb scattering term in Gluckstern formula
805  float fracSig_MCS = (0.052*sqrt(1.43)) / (gastpc_B * sqrt(gastpc_X0*trkLenPerp.at(_nFSP)*0.0001));
806  // momentum resoltion from the two terms above
807  float sigmaP = truep.at(_nFSP) * sqrt( fracSig_meas*fracSig_meas + fracSig_MCS*fracSig_MCS );
808  // now Gaussian smear the true momentum using the momentum resolution
809  preco = fHelper->GaussianSmearing( truep.at(_nFSP), sigmaP );
810 
811  // measurement term in the Gluckstern formula for calculating the
812  // angular resolution
813  float sigma_angle_1 = ((sigma_x * sigma_x * 0.0001) / trkLen.at(_nFSP)*trkLen.at(_nFSP)*0.0001) * (12*(nHits-1))/(nHits*(nHits+1));
814  // scattering term in Gluckstern formula
815  float sigma_angle_2 = (0.015*0.015 / (3. * truep.at(_nFSP) * truep.at(_nFSP))) * (trkLen.at(_nFSP)/gastpc_X0);
816  // angular resolution from the two terms above
817  float sigma_angle = sqrt(sigma_angle_1 + sigma_angle_2);
818  // now Gaussian smear the true angle using the angular resolution
819  angle_reco = fHelper->GaussianSmearing(_angle.at(_nFSP), sigma_angle);
820  }
float GaussianSmearing(const float mean, const float sigma)
std::vector< double > _angle
const float sigma_x
const float gastpc_padPitch
std::vector< double > trkLen
std::vector< double > truep
unsigned int _nFSP
const float gastpc_X0
CAFHelper * fHelper
const float gastpc_B
std::vector< double > trkLenPerp
void gar::ParamSim::DoRangeCalculation ( float &  preco,
float &  angle_reco 
)
private

Definition at line 779 of file ParamSim_module.cc.

780  {
781  // calculate number of trackpoints
782  float nHits = round (trkLen.at(_nFSP) / gastpc_padPitch);
783  // angular resolution first term
784  float sigma_angle_1 = ((sigma_x * sigma_x * 0.0001) / trkLen.at(_nFSP)*trkLen.at(_nFSP)*0.0001) * (12*(nHits-1))/(nHits*(nHits+1));
785  // scattering term in Gluckstern formula
786  float sigma_angle_2 = (0.015*0.015 / (3. * truep.at(_nFSP) * truep.at(_nFSP))) * (trkLen.at(_nFSP)/gastpc_X0);
787  // angular resolution from the two terms above
788  float sigma_angle_short = sqrt(sigma_angle_1 + sigma_angle_2);
789 
790  //reconstructed angle
791  angle_reco = fHelper->GaussianSmearing(_angle.at(_nFSP), sigma_angle_short);
792  //reconstructed momentum
794  }
const float sigmaP_short
float GaussianSmearing(const float mean, const float sigma)
std::vector< double > _angle
const float sigma_x
const float gastpc_padPitch
std::vector< double > trkLen
std::vector< double > truep
unsigned int _nFSP
const float gastpc_X0
CAFHelper * fHelper
void gar::ParamSim::FillCommonVariables ( const int  pdg,
const TLorentzVector &  momentum,
const TLorentzVector &  position,
const TLorentzVector &  positionEnd,
const int  mctrackid,
const int  mothertrackid,
const int  motherpdg,
const float  ptrue,
const float  angle,
const std::string  mcp_process,
const std::string  mcp_endprocess,
const float  time 
)
private

Definition at line 683 of file ParamSim_module.cc.

684  {
685  const TVector3 spoint(position.X() - fOrigin[0], position.Y() - fOrigin[1], position.Z() - fOrigin[2]);
686  const TVector3 epoint(positionEnd.X() - fOrigin[0], positionEnd.Y() - fOrigin[1], positionEnd.Z() - fOrigin[2]);
687 
688  truepdg.push_back(pdg);
689  truepx.push_back(momentum.X());
690  truepy.push_back(momentum.Y());
691  truepz.push_back(momentum.Z());
692 
693  if(fCorrect4origin){
694  _MCPStartX.push_back(position.X() - fOrigin[0]);
695  _MCPStartY.push_back(position.Y() - fOrigin[1]);
696  _MCPStartZ.push_back(position.Z() - fOrigin[2]);
697  _MCPEndX.push_back(positionEnd.X() - fOrigin[0]);
698  _MCPEndY.push_back(positionEnd.Y() - fOrigin[1]);
699  _MCPEndZ.push_back(positionEnd.Z() - fOrigin[2]);
700  } else {
701  _MCPStartX.push_back(position.X());
702  _MCPStartY.push_back(position.Y());
703  _MCPStartZ.push_back(position.Z());
704  _MCPEndX.push_back(positionEnd.X());
705  _MCPEndY.push_back(positionEnd.Y());
706  _MCPEndZ.push_back(positionEnd.Z());
707  }
708 
709  // save the true momentum
710  truep.push_back(ptrue);
711  // save the true angle
712  _angle.push_back(angle);
713  //Save MC process
714  _MCProc.push_back(mcp_process);
715  _MCEndProc.push_back(mcp_endprocess);
716  mctime.push_back(time);
717  mctrkid.push_back(mctrackid);
718  motherid.push_back(mothertrackid);
719  pdgmother.push_back(motherpdg);
720 
721  isFidStart.push_back(fHelper->PointInFiducial(spoint));
722  isTPCStart.push_back(fHelper->PointInTPC(spoint));
723  isCaloStart.push_back(fHelper->PointInCalo(spoint));
724  isThroughCaloStart.push_back(fHelper->isThroughCalo(spoint));
725  isInBetweenStart.push_back(fHelper->PointStopBetween(spoint));
726  isBarrelStart.push_back(fHelper->isBarrel(spoint));
727  isEndcapStart.push_back(fHelper->isEndcap(spoint));
728  //Check where endpoint of mcp is
729  isFidEnd.push_back(fHelper->PointInFiducial(epoint));
730  isTPCEnd.push_back(fHelper->PointInTPC(epoint));
731  isCaloEnd.push_back(fHelper->PointInCalo(epoint));
732  isThroughCaloEnd.push_back(fHelper->isThroughCalo(epoint));
733  isInBetweenEnd.push_back(fHelper->PointStopBetween(epoint));
734  isBarrelEnd.push_back(fHelper->isBarrel(epoint));
735  isEndcapEnd.push_back(fHelper->isEndcap(epoint));
736  }
bool PointInFiducial(const TVector3 &point)
std::vector< double > truepz
std::vector< unsigned int > isFidEnd
std::vector< std::string > _MCEndProc
std::vector< unsigned int > isBarrelStart
std::vector< std::string > _MCProc
bool PointInTPC(const TVector3 &point)
std::vector< double > _angle
std::vector< double > truepx
std::vector< unsigned int > isTPCStart
bool isBarrel(const TVector3 &point)
std::vector< int > motherid
std::vector< int > _MCPEndY
std::vector< unsigned int > isEndcapEnd
std::vector< int > truepdg
std::vector< unsigned int > isTPCEnd
std::vector< int > mctrkid
std::vector< unsigned int > isFidStart
bool PointStopBetween(const TVector3 &point)
std::vector< int > _MCPEndZ
std::vector< unsigned int > isThroughCaloEnd
std::vector< double > truepy
bool isThroughCalo(const TVector3 &point)
bool isEndcap(const TVector3 &point)
std::vector< unsigned int > isEndcapStart
std::vector< unsigned int > isCaloStart
std::vector< unsigned int > isThroughCaloStart
std::vector< unsigned int > isCaloEnd
std::vector< int > _MCPStartZ
std::vector< double > truep
std::vector< unsigned int > isBarrelEnd
std::vector< int > _MCPStartY
CAFHelper * fHelper
std::vector< double > mctime
def momentum(x1, x2, x3, scale=1.)
bool PointInCalo(const TVector3 &point)
std::vector< int > _MCPEndX
std::vector< unsigned int > isInBetweenEnd
std::vector< unsigned int > isInBetweenStart
std::vector< int > _MCPStartX
std::vector< int > pdgmother
ParamSim& gar::ParamSim::operator= ( ParamSim const &  )
delete
ParamSim& gar::ParamSim::operator= ( ParamSim &&  )
delete
void gar::ParamSim::SaveGtruthMCtruth ( art::Event const &  e)
private

Definition at line 556 of file ParamSim_module.cc.

557  {
558  auto MCTHandle = e.getHandle< std::vector<simb::MCTruth> >(fGeneratorLabel);
559  if (!MCTHandle) {
560  throw cet::exception("ParamSim") << " No simb::MCTruth branch."
561  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
562  }
563 
564  auto GTHandle = e.getHandle< std::vector<simb::GTruth> >(fGeneratorLabel);
565  if (!GTHandle) {
566  throw cet::exception("ParamSim") << " No simb::GTruth branch."
567  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
568  }
569 
570  // save MCTruth info
571  for ( auto const& mct : (*MCTHandle) ) {
572  if (mct.NeutrinoSet()) {
573  simb::MCNeutrino nuw = mct.GetNeutrino();
574  ccnc.push_back(nuw.CCNC());
575  ntype.push_back(nuw.Nu().PdgCode());
576  q2.push_back(nuw.QSqr());
577  w.push_back(nuw.W());
578  y.push_back(nuw.Y());
579  x.push_back(nuw.X());
580  theta.push_back(nuw.Theta());
581  mode.push_back(nuw.Mode());
582  intert.push_back(nuw.InteractionType());
583  if(fCorrect4origin){
584  vertx.push_back(nuw.Nu().EndX() - fOrigin[0]);
585  verty.push_back(nuw.Nu().EndY() - fOrigin[1]);
586  vertz.push_back(nuw.Nu().EndZ() - fOrigin[2]);
587  } else {
588  vertx.push_back(nuw.Nu().EndX());
589  verty.push_back(nuw.Nu().EndY());
590  vertz.push_back(nuw.Nu().EndZ());
591  }
592  mcnupx.push_back(nuw.Nu().Px());
593  mcnupy.push_back(nuw.Nu().Py());
594  mcnupz.push_back(nuw.Nu().Pz());
595  } // end MC info from MCTruth
596  }
597 
598  // save GTruth info
599  for ( auto const& gt : (*GTHandle) ) {
600  gint.push_back(gt.fGint);
601  tgtpdg.push_back(gt.ftgtPDG);
602  weight.push_back(gt.fweight);
603  gt_t.push_back(gt.fgT);
604  }
605  }
int PdgCode() const
Definition: MCParticle.h:212
int CCNC() const
Definition: MCNeutrino.h:148
std::vector< int > weight
double QSqr() const
Definition: MCNeutrino.h:157
std::vector< int > ccnc
double Theta() const
angle between incoming and outgoing leptons, in radians
Definition: MCNeutrino.cxx:63
double Py(const int i=0) const
Definition: MCParticle.h:231
double EndZ() const
Definition: MCParticle.h:228
std::vector< double > y
std::vector< double > mcnupy
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
std::vector< int > ntype
double Px(const int i=0) const
Definition: MCParticle.h:230
std::string fGeneratorLabel
double EndY() const
Definition: MCParticle.h:227
std::vector< int > gint
std::vector< double > vertx
std::vector< double > mcnupx
std::vector< double > vertz
int InteractionType() const
Definition: MCNeutrino.h:150
const double e
std::vector< double > verty
double W() const
Definition: MCNeutrino.h:154
double Y() const
Definition: MCNeutrino.h:156
std::vector< int > gt_t
std::vector< double > w
double X() const
Definition: MCNeutrino.h:155
std::vector< int > mode
std::vector< double > q2
std::vector< double > x
std::vector< double > mcnupz
double Pz(const int i=0) const
Definition: MCParticle.h:232
std::vector< double > theta
std::vector< int > intert
double EndX() const
Definition: MCParticle.h:226
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149
std::vector< int > tgtpdg
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
void gar::ParamSim::TPCParticleIdentification ( )
private

Definition at line 1054 of file ParamSim_module.cc.

1055  {
1056  //start pid
1057  //***************************************************************************************************************/
1058 
1059  int pdg = truepdg.at(_nFSP);
1060 
1061  for (int pidm = 0; pidm < 6; ++pidm)
1062  {
1063  if ( abs(pdg) == pdg_charged.at(pidm) )
1064  {
1065  float p = _preco.at(_nFSP);
1066  std::vector<double> vec;
1067  std::vector<std::string> pnamelist = {"#pi", "#mu", "p", "K", "d", "e"};
1068  std::vector<std::string> recopnamelist = {"#pi", "#mu", "p", "K", "d", "e"};
1069 
1070  int qclosest = 0;
1071  float dist = 100000000.;
1072 
1073  for (int q = 0; q < 501; ++q)
1074  {
1075  //Check the title and the reco momentum take only the one that fits
1076  std::string fulltitle = m_pidinterp[q]->GetTitle();
1077  unsigned first = fulltitle.find("=");
1078  unsigned last = fulltitle.find("GeV");
1079  std::string substr = fulltitle.substr(first+1, last - first-1);
1080  float pidinterp_mom = std::atof(substr.c_str());
1081  //calculate the distance between the bin and mom, store the q the closest
1082  float disttemp = std::abs(pidinterp_mom - p);
1083 
1084  if( disttemp < dist ) {
1085  dist = disttemp;
1086  qclosest = q;
1087  }
1088  } // closes the "pidmatrix" loop
1089 
1090  //loop over the columns (true pid)
1091  std::vector< std::pair<float, std::string> > v_prob;
1092  //get true particle name
1093  std::string trueparticlename = m_pidinterp[qclosest]->GetXaxis()->GetBinLabel(pidm+1);
1094  // std::cout << trueparticlename << std::endl;
1095  if ( trueparticlename == pnamelist[pidm] )
1096  {
1097  //loop over the rows (reco pid)
1098  for (int pidr = 0; pidr < 6; ++pidr)
1099  {
1100  std::string recoparticlename = m_pidinterp[qclosest]->GetYaxis()->GetBinLabel(pidr+1);
1101  // std::cout << recoparticlename << std::endl;
1102  if (recoparticlename == recopnamelist[pidr])
1103  {
1104  float prob = m_pidinterp[qclosest]->GetBinContent(pidm+1,pidr+1);
1105  prob_arr.push_back(prob);
1106  //Need to check random number value and prob value then associate the recopdg to the reco prob
1107  v_prob.push_back( std::make_pair(prob, recoparticlename) );
1108  }
1109  }
1110 
1111  int pid = -1;
1112  if(v_prob.size() > 1)
1113  {
1114  //Order the vector of prob
1115  std::sort(v_prob.begin(), v_prob.end());
1116  //Throw a random number between 0 and 1
1117  float random_number = fHelper->GetRamdomNumber();
1118  //Make cumulative sum to get the range
1119  std::partial_sum(v_prob.begin(), v_prob.end(), v_prob.begin(), [](const std::pair<float, std::string>& _x, const std::pair<float, std::string>& _y){return std::pair<float, std::string>(_x.first + _y.first, _y.second);});
1120  for(size_t ivec = 0; ivec < v_prob.size()-1; ivec++)
1121  {
1122  if( random_number < v_prob.at(ivec+1).first && random_number >= v_prob.at(ivec).first )
1123  {
1124  pid = pdg_charged.at( std::distance( recopnamelist.begin(), std::find(recopnamelist.begin(), recopnamelist.end(), v_prob.at(ivec+1).second) ) );
1125  }
1126  }
1127  }
1128  else
1129  {
1130  pid = pdg_charged.at( std::distance( recopnamelist.begin(), std::find(recopnamelist.begin(), recopnamelist.end(), v_prob.at(0).second) ) );
1131  }
1132 
1133  recopid.push_back( pid );
1134  } // closes the if statement
1135  } // end if charged in pdg list
1136  } // end loop pidm
1137 
1138  //end pid
1139  //***************************************************************************************************************/
1140  }
std::vector< double > _preco
std::string string
Definition: nybbler.cc:12
std::unordered_map< int, TH2F * > m_pidinterp
T abs(T value)
std::vector< int > truepdg
p
Definition: test.py:223
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
std::vector< double > prob_arr
std::vector< int > recopid
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
unsigned int _nFSP
CAFHelper * fHelper
const std::vector< int > pdg_charged
void gar::ParamSim::TreatMCParticles ( art::Event const &  e)
private

Definition at line 608 of file ParamSim_module.cc.

609  {
610  auto MCPHandle = e.getHandle< std::vector<simb::MCParticle> >(fGeantLabel);
611  if (!MCPHandle) {
612  throw cet::exception("anatree") << " No simb::MCParticle branch."
613  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
614  }
615 
616  std::unordered_map<int, int> TrackIdToIndex;
617  int index = 0;
618  for ( auto const& mcp : (*MCPHandle) ) {
619  int TrackId = mcp.TrackId();
620  TrackIdToIndex[TrackId] = index++;
621  }
622 
623  //Loop over the mcp
624  for ( auto const& mcp : (*MCPHandle) ) {
625 
626  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
627  const TParticlePDG* definition = databasePDG->GetParticle( mcp.PdgCode() );
628  if (definition == nullptr) continue;
629 
630  const std::string mcp_process = mcp.Process();
631  const std::string mcp_endprocess = mcp.EndProcess();
632  const int mctrackid = mcp.TrackId();
633  const int mothertrackid = mcp.Mother();
634  const int pdg = mcp.PdgCode();
635  const TLorentzVector& position = mcp.Position(0);
636  const TVector3 spoint(position.X() - fOrigin[0], position.Y() - fOrigin[1], position.Z() - fOrigin[2]);
637  const TLorentzVector& positionEnd = mcp.EndPosition();
638  const TVector3 epoint(positionEnd.X() - fOrigin[0], positionEnd.Y() - fOrigin[1], positionEnd.Z() - fOrigin[2]);
639  const TLorentzVector& momentum = mcp.Momentum(0);
640  const TVector3 mom(momentum.X(), momentum.Y(), momentum.Z());
641  const float ptrue = mom.Mag();
642  const float angle = atan(mom.X() / mom.Z());
643  const float time = mcp.T();
644  int motherpdg = 0;
645 
646  //Find out mother pdg
647  if(TrackIdToIndex.find(mothertrackid) != TrackIdToIndex.end()) {
648  motherpdg = (*MCPHandle).at(TrackIdToIndex[mothertrackid]).PdgCode();
649  }
650 
651  //need to ignore neutrals for this - put the value to 0
652  auto result = std::find(pdg_neutral.begin(), pdg_neutral.end(), abs(pdg));
653  bool isNeutral = (result != pdg_neutral.end()) ? true : false;
654 
655  if( isNeutral )
656  {
657  trkLen.push_back(-1);
658  trkLenPerp.push_back(-1);
659  }
660  else {
661  ComputeTrkLength(mcp);
662  }
663 
664  float ecaltime = fHelper->GaussianSmearing(time, fECALTimeResolution);
665 
666  //Store common mcp properties
667  FillCommonVariables(pdg, momentum, position, positionEnd, mctrackid, mothertrackid, motherpdg, ptrue, angle, mcp_process, mcp_endprocess, time);
668 
669  //Treat visible particles in the TPC
670  if( trkLen.at(_nFSP) > gastpc_len ) {
671  TreatTPCVisible(ecaltime);
672  }
673  else
674  {
675  TreatTPCNotVisible(ecaltime);
676  }// end trkLen.at(_nFSP) < gastpc_len
677 
678  _nFSP++;
679  }//end loop mcp
680  }
void TreatTPCVisible(const float ecaltime)
static QCString result
std::string string
Definition: nybbler.cc:12
float GaussianSmearing(const float mean, const float sigma)
void ComputeTrkLength(const simb::MCParticle &mcp)
std::string fGeantLabel
module label for geant4 simulated hits
T abs(T value)
const double e
const std::vector< int > pdg_neutral
std::vector< double > trkLen
void FillCommonVariables(const int pdg, const TLorentzVector &momentum, const TLorentzVector &position, const TLorentzVector &positionEnd, const int mctrackid, const int mothertrackid, const int motherpdg, const float ptrue, const float angle, const std::string mcp_process, const std::string mcp_endprocess, const float time)
unsigned int _nFSP
CAFHelper * fHelper
const float fECALTimeResolution
const float gastpc_len
def momentum(x1, x2, x3, scale=1.)
void TreatTPCNotVisible(const float ecaltime)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
std::vector< double > trkLenPerp
void gar::ParamSim::TreatNeutrons ( const float  ecaltime)
private

Definition at line 1187 of file ParamSim_module.cc.

1188  {
1189  TVector3 epoint(_MCPEndX.at(_nFSP), _MCPEndY.at(_nFSP), _MCPEndZ.at(_nFSP));
1190  //start neutrons
1191  //***************************************************************************************************************/
1192 
1193  if(fHelper->PointInCalo(epoint)) //needs to stop in the ECAL
1194  {
1195  //check if it can be detected by the ECAL
1196  //Assumes 40% efficiency to detect
1197  float random_number = fHelper->GetRamdomNumber();
1198  float true_KE = std::sqrt(truep.at(_nFSP)*truep.at(_nFSP) + neutron_mass*neutron_mass) - neutron_mass;
1199  // float true_KE = ptrue*ptrue / (2*neutron_mass); // in GeV
1200  int index = (true_KE >= 0.05) ? 1 : 0;
1201 
1202  if(random_number > (1 - NeutronECAL_detEff[index]) && true_KE > 0.003)//Threshold of 3 MeV
1203  {
1204  //TODO random is first interaction or rescatter and smear accordingly to Chris's study
1205  //Detected in the ECAL
1206  recopid.push_back(-1); //reco pid set to 0?
1207  detected.push_back(1);
1208  float eres = sigmaNeutronECAL_first * true_KE;
1209  float ereco = fHelper->GaussianSmearing( true_KE, eres );
1210  erecon.push_back(ereco > 0 ? ereco : 0.);
1211  etime.push_back(ecaltime);
1212  _preco.push_back(-1);
1213  anglereco.push_back(-1);
1214  recopidecal.push_back(2112);
1215  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1216  }
1217  else
1218  {
1219  //neutron not detected
1220  detected.push_back(0);
1221  recopid.push_back(-1);
1222  erecon.push_back(-1);
1223  etime.push_back(-1);
1224  _preco.push_back(-1);
1225  anglereco.push_back(-1);
1226  recopidecal.push_back(-1);
1227  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1228  }
1229  } //endpoint is in ECAL
1230  else
1231  {
1232  //Endpoint is not in calo (TPC/isInBetween or outside Calo)
1233  detected.push_back(0);
1234  recopid.push_back(-1);
1235  erecon.push_back(-1);
1236  etime.push_back(-1);
1237  _preco.push_back(-1);
1238  anglereco.push_back(-1);
1239  recopidecal.push_back(-1);
1240  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1241  }
1242 
1243  //End neutrons
1244  //***************************************************************************************************************/
1245  }
std::vector< double > _preco
const float sigmaNeutronECAL_first
std::vector< int > detected
float GaussianSmearing(const float mean, const float sigma)
std::vector< int > recopidecal
std::vector< int > _MCPEndY
const float NeutronECAL_detEff[2]
std::vector< int > _MCPEndZ
const float neutron_mass
std::vector< double > prob_arr
std::vector< int > recopid
std::vector< double > truep
unsigned int _nFSP
CAFHelper * fHelper
bool PointInCalo(const TVector3 &point)
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< double > etime
void gar::ParamSim::TreatOthers ( const float  ecaltime)
private

Definition at line 1323 of file ParamSim_module.cc.

1324  {
1325  int pdg = truepdg.at(_nFSP);
1326  TVector3 epoint(_MCPEndX.at(_nFSP), _MCPEndY.at(_nFSP), _MCPEndZ.at(_nFSP));
1327  //Case for particles that stop or go through ECAL (problematic particles with no track length????)
1328  //Not visible in the TPC and not neutron or gamma or pi0 (otherwise it has been already done above)
1329  if(fHelper->PointInCalo(epoint))
1330  {
1331  detected.push_back(1);
1332  etime.push_back(ecaltime);
1333 
1334  TParticlePDG *part = TDatabasePDG::Instance()->GetParticle(abs(pdg));
1335  float mass = 0.;
1336  if(nullptr != part)
1337  mass = part->Mass();//in GeV
1338 
1339  float etrue = std::sqrt(truep.at(_nFSP)*truep.at(_nFSP) + mass*mass) - mass;
1340  float ECAL_resolution = fRes->Eval(etrue)*etrue;
1341  float ereco = fHelper->GaussianSmearing(etrue, ECAL_resolution);
1342  erecon.push_back((ereco > 0) ? ereco : 0.);
1343 
1344  //Electron
1345  if( abs(pdg) == 11 ){
1346  recopidecal.push_back(11);
1347  }
1348  else if( abs(pdg) == 13 || abs(pdg) == 211 )
1349  {
1350  //Muons and Pions
1351  //ptrue < 480 MeV/c 100% separation
1352  //80% from 480 to 750
1353  //90% up to 750 to 900
1354  //95% over 900
1355  float random_number = fHelper->GetRamdomNumber();
1356 
1357  if(truep.at(_nFSP) < 0.48) {
1358  recopidecal.push_back(abs(pdg));//100% efficiency by range
1359  }
1360  else if(truep.at(_nFSP) >= 0.48 && truep.at(_nFSP) < 0.75)
1361  {
1362  //case muon
1363  if(abs(pdg) == 13)
1364  {
1365  if(random_number > (1 - 0.8)) {
1366  recopidecal.push_back(13);
1367  }
1368  else{
1369  recopidecal.push_back(211);
1370  }
1371  }
1372 
1373  //case pion
1374  if(abs(pdg) == 211)
1375  {
1376  if(random_number > (1 - 0.8)) {
1377  recopidecal.push_back(211);
1378  }
1379  else{
1380  recopidecal.push_back(13);
1381  }
1382  }
1383  }
1384  else if(truep.at(_nFSP) >= 0.75 && truep.at(_nFSP) < 0.9)
1385  {
1386  //case muon
1387  if(abs(pdg) == 13){
1388  if(random_number > (1 - 0.9)) {
1389  recopidecal.push_back(13);
1390  }
1391  else{
1392  recopidecal.push_back(211);
1393  }
1394  }
1395  //case pion
1396  if(abs(pdg) == 211) {
1397  if(random_number > (1 - 0.9)) {
1398  recopidecal.push_back(211);
1399  }
1400  else{
1401  recopidecal.push_back(13);
1402  }
1403  }
1404  }
1405  else
1406  {
1407  //case muon
1408  if(abs(pdg) == 13){
1409  if(random_number > (1 - 0.95)) {
1410  recopidecal.push_back(13);
1411  }
1412  else{
1413  recopidecal.push_back(211);
1414  }
1415  }
1416  //case pion
1417  if(abs(pdg) == 211){
1418  if(random_number > (1 - 0.95)) {
1419  recopidecal.push_back(211);
1420  }
1421  else{
1422  recopidecal.push_back(13);
1423  }
1424  }
1425  }
1426  }
1427  else if( abs(pdg) == 2212 )
1428  {
1429  recopidecal.push_back(2212);//TODO for p/pi separation
1430  }
1431  else {
1432  recopidecal.push_back(-1);
1433  }
1434 
1435  _preco.push_back(-1);
1436  anglereco.push_back(-1);
1437  recopid.push_back(-1);
1438  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1439  }
1440  else if (fHelper->isThroughCalo(epoint))
1441  {
1442  //Case the endpoint is outside the CALO -> it went through the ECAL (mu/pi/p possible)
1443  //the ECAL will see 60 MIPs on average
1444  double Evis = (double)nLayers; //in MIP
1445  //Smearing to account for Gaussian detector noise (Landau negligible)
1446  Evis = fHelper->GaussianSmearing(Evis, ECAL_MIP_Res);
1447  //1 MIP = 0.814 MeV
1448  double Erec = Evis * MIP2GeV_factor;
1449  erecon.push_back((Erec > 0) ? Erec : 0.);
1450  etime.push_back(ecaltime);
1451  detected.push_back(1);
1452 
1453  //Muon/Pions/Protons are reco as Muons (without MuID detector)
1454  if( abs(pdg) == 13 || abs(pdg) == 211 || abs(pdg) == 2212 ) {
1455  recopidecal.push_back(13);
1456  }
1457  else {
1458  recopidecal.push_back(-1);
1459  }
1460 
1461  _preco.push_back(-1);
1462  anglereco.push_back(-1);
1463  recopid.push_back(-1);
1464  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1465  }
1466  else if(fHelper->PointInTPC(epoint) || fHelper->PointStopBetween(epoint))
1467  {
1468  detected.push_back(0);
1469  etime.push_back(-1);
1470  erecon.push_back(-1);
1471  _preco.push_back(-1);
1472  anglereco.push_back(-1);
1473  recopid.push_back(-1);
1474  recopidecal.push_back(-1);
1475  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1476  }
1477  }
std::vector< double > _preco
std::vector< int > detected
const double ECAL_MIP_Res
const double MIP2GeV_factor
float GaussianSmearing(const float mean, const float sigma)
bool PointInTPC(const TVector3 &point)
std::vector< int > recopidecal
std::vector< int > _MCPEndY
T abs(T value)
std::vector< int > truepdg
bool PointStopBetween(const TVector3 &point)
std::vector< int > _MCPEndZ
bool isThroughCalo(const TVector3 &point)
std::vector< double > prob_arr
std::vector< int > recopid
std::vector< double > truep
unsigned int _nFSP
CAFHelper * fHelper
bool PointInCalo(const TVector3 &point)
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< double > etime
void gar::ParamSim::TreatPhotons ( const float  ecaltime)
private

Definition at line 1248 of file ParamSim_module.cc.

1249  {
1250  TVector3 epoint(_MCPEndX.at(_nFSP), _MCPEndY.at(_nFSP), _MCPEndZ.at(_nFSP));
1251  //start gammas
1252  //***************************************************************************************************************/
1253 
1254  if( pdgmother.at(_nFSP) != 111 )
1255  {
1256  //Endpoint is in the ECAL
1257  if(fHelper->PointInCalo(epoint))
1258  {
1259  //if they hit the ECAL and smear their energy
1260  float ECAL_resolution = fRes->Eval(truep.at(_nFSP))*truep.at(_nFSP);
1261  float ereco = fHelper->GaussianSmearing(truep.at(_nFSP), ECAL_resolution);
1262  erecon.push_back( (ereco > 0) ? ereco : 0. );
1263  recopid.push_back(-1);
1264  detected.push_back(1);
1265  etime.push_back(ecaltime);
1266  _preco.push_back(-1);
1267  anglereco.push_back(-1);
1268  //reach the ECAL, should be tagged as gamma
1269  recopidecal.push_back(22);
1270  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1271  }
1272  else if(fHelper->PointInTPC(epoint) || fHelper->PointStopBetween(epoint) || fHelper->isThroughCalo(epoint))
1273  {
1274  erecon.push_back(-1);
1275  recopid.push_back(-1);
1276  detected.push_back(0);
1277  etime.push_back(-1);
1278  _preco.push_back(-1);
1279  anglereco.push_back(-1);
1280  //converted so not seen in ECAL
1281  recopidecal.push_back(-1);
1282  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1283  }
1284  }
1285  else
1286  {
1287  //From a pi0
1288  if(fHelper->PointInCalo(epoint))
1289  {
1290  //if they hit the ECAL and smear their energy
1291  float ECAL_resolution = fRes->Eval(truep.at(_nFSP))*truep.at(_nFSP);
1292  float ereco = fHelper->GaussianSmearing(truep.at(_nFSP), ECAL_resolution);
1293  erecon.push_back((ereco > 0) ? ereco : 0.);
1294  recopid.push_back(-1);
1295  detected.push_back(1);
1296  etime.push_back(ecaltime);
1297  _preco.push_back(-1);
1298  anglereco.push_back(-1);
1299  //reaches the ecal
1300  recopidecal.push_back(22);
1301  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1302  }
1303  else if(fHelper->PointInTPC(epoint) || fHelper->PointStopBetween(epoint) || fHelper->isThroughCalo(epoint))
1304  {
1305  //from pi0 and converted in TPC or stopped between TPC/ECAL
1306  erecon.push_back(-1);
1307  recopid.push_back(-1);
1308  detected.push_back(0);
1309  etime.push_back(-1);
1310  _preco.push_back(-1);
1311  anglereco.push_back(-1);
1312  //converted not seen by ecal
1313  recopidecal.push_back(-1);
1314  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1315  }
1316  }
1317 
1318  //end gammas
1319  //***************************************************************************************************************/
1320  }
std::vector< double > _preco
std::vector< int > detected
float GaussianSmearing(const float mean, const float sigma)
bool PointInTPC(const TVector3 &point)
std::vector< int > recopidecal
std::vector< int > _MCPEndY
bool PointStopBetween(const TVector3 &point)
std::vector< int > _MCPEndZ
bool isThroughCalo(const TVector3 &point)
std::vector< double > prob_arr
std::vector< int > recopid
std::vector< double > truep
unsigned int _nFSP
CAFHelper * fHelper
bool PointInCalo(const TVector3 &point)
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< int > pdgmother
std::vector< double > etime
void gar::ParamSim::TreatTPCNotVisible ( const float  ecaltime)
private

Definition at line 1143 of file ParamSim_module.cc.

1144  {
1145  int pdg = truepdg.at(_nFSP);
1146 
1147  //Case neutrinos
1148  if(std::find(neutrinos.begin(), neutrinos.end(), std::abs(pdg)) != neutrinos.end())
1149  {
1150  detected.push_back(0);
1151  etime.push_back(0.);
1152  erecon.push_back(0.);
1153  recopidecal.push_back(-1);
1154  _preco.push_back(-1);
1155  anglereco.push_back(-1);
1156  recopid.push_back(-1);
1157  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1158  }
1159  else if( std::abs(pdg) == 2112 )
1160  {
1161  TreatNeutrons(ecaltime);
1162  }
1163  else if(std::abs(pdg) == 111)
1164  {
1165  //Pi0 case
1166  erecon.push_back(-1);
1167  recopid.push_back(-1);
1168  detected.push_back(0);
1169  recopidecal.push_back(-1);
1170  etime.push_back(-1);
1171  _preco.push_back(-1);
1172  anglereco.push_back(-1);
1173 
1174  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1175  }
1176  else if(std::abs(pdg) == 22)
1177  {
1178  TreatPhotons(ecaltime);
1179  }
1180  else
1181  {
1182  TreatOthers(ecaltime);
1183  }//end case charged but not visible in TPC
1184  }
std::vector< double > _preco
void TreatPhotons(const float ecaltime)
std::vector< int > detected
std::vector< int > recopidecal
void TreatNeutrons(const float ecaltime)
void TreatOthers(const float ecaltime)
T abs(T value)
std::vector< int > truepdg
std::vector< double > prob_arr
std::vector< int > recopid
const std::vector< int > neutrinos
unsigned int _nFSP
std::vector< double > anglereco
std::vector< double > erecon
std::vector< double > etime
void gar::ParamSim::TreatTPCVisible ( const float  ecaltime)
private

Definition at line 823 of file ParamSim_module.cc.

824  {
825  //start tpc
826  //***************************************************************************************************************/
827  int pdg = truepdg.at(_nFSP);
828 
829  if ( std::find(pdg_charged.begin(), pdg_charged.end(), std::abs(pdg)) != pdg_charged.end() )
830  {
831  TVector3 epoint(_MCPEndX.at(_nFSP), _MCPEndY.at(_nFSP), _MCPEndZ.at(_nFSP));
832  float ptrue = truep.at(_nFSP);
833 
834  //Use range instead of Gluckstern for stopping tracks
835  //TODO is that correct? What if it is a scatter in the TPC? Need to check if daughter is same particle
836  float preco = 0.;
837  float angle_reco = 0.;
838 
839  //Case for range, the end point of the mcp is in the TPC, does not reach the ecal
840  if( fHelper->PointInTPC(epoint) )
841  {
842  DoRangeCalculation(preco, angle_reco);
843 
844  if(preco > 0)
845  _preco.push_back(preco);
846  else
847  _preco.push_back(-1);
848  anglereco.push_back(angle_reco);
849  erecon.push_back(-1);
850  recopidecal.push_back(-1);
851  etime.push_back(-1);
852  detected.push_back(-1);
853  }
854  else
855  {
856  DoGluckSternCalculation(preco, angle_reco);
857 
858  // save reconstructed momentum and angle to cafanatree
859  if(preco > 0)
860  _preco.push_back(preco);
861  else
862  _preco.push_back(-1);
863  anglereco.push_back(angle_reco);
864 
865  //Reaches the ECAL and stops there
866  if( fHelper->PointInCalo(epoint) )
867  {
868  //Need energy measurement in ecal
869  TParticlePDG *part = TDatabasePDG::Instance()->GetParticle(abs(pdg));
870  if(nullptr == part)
871  {
872  //deuteron
873  if( pdg == 1000010020 )
874  {
875  float mass = 1.8756;//in GeV mass deuteron
876  float etrue = std::sqrt(ptrue*ptrue + mass*mass) - mass;
877  float ECAL_resolution = fRes->Eval(etrue)*etrue;
878  float ereco = fHelper->GaussianSmearing(etrue, ECAL_resolution);
879  erecon.push_back((ereco > 0) ? ereco : 0.);
880  recopidecal.push_back(-1);
881  detected.push_back(1);
882  etime.push_back(ecaltime);
883  }
884  else
885  {
886  erecon.push_back(-1);
887  recopidecal.push_back(-1);
888  detected.push_back(0);
889  etime.push_back(-1);
890  }
891  }
892  else
893  {
894  //by default should be tagged as an electron as it has a track,
895  //otherwise tag as gamma if not track -> need mis association rate, and use dE/dX in Scintillator?
896  //separation between e and mu/pi should be around 100%
897  //separation mu/pi -> based on Chris study with only the ECAL (no Muon ID detector)
898  //separation with p and mu/pi/e ?? high energy -> confusion with mu/pi, low energy confusion with e
899  //using E/p to ID?
900  float mass = part->Mass();//in GeV
901  float etrue = std::sqrt(ptrue*ptrue + mass*mass) - mass;
902  float ECAL_resolution = fRes->Eval(etrue)*etrue;
903  float ereco = fHelper->GaussianSmearing(etrue, ECAL_resolution);
904  erecon.push_back((ereco > 0) ? ereco : 0.);
905  detected.push_back(1);
906  etime.push_back(ecaltime);
907 
908  //Electron
909  if( abs(pdg) == 11 ){
910  recopidecal.push_back(11);
911  }
912  else if( abs(pdg) == 13 || abs(pdg) == 211 )
913  {
914  //Muons and Pions
915  //ptrue < 480 MeV/c 100% separation
916  //80% from 480 to 750
917  //90% up to 750 to 900
918  //95% over 900
919  float random_number = fHelper->GetRamdomNumber();
920 
921  if(ptrue < 0.48) {
922  recopidecal.push_back(abs(pdg));//100% efficiency by range
923  }
924  else if(ptrue >= 0.48 && ptrue < 0.75)
925  {
926  //case muon
927  if(abs(pdg) == 13)
928  {
929  if(random_number > (1 - 0.8)) {
930  recopidecal.push_back(13);
931  }
932  else{
933  recopidecal.push_back(211);
934  }
935  }
936  //case pion
937  if(abs(pdg) == 211)
938  {
939  if(random_number > (1 - 0.8)) {
940  recopidecal.push_back(211);
941  }
942  else{
943  recopidecal.push_back(13);
944  }
945  }
946  }
947  else if(ptrue >= 0.75 && ptrue < 0.9)
948  {
949  //case muon
950  if(abs(pdg) == 13){
951  if(random_number > (1 - 0.9)) {
952  recopidecal.push_back(13);
953  }
954  else{
955  recopidecal.push_back(211);
956  }
957  }
958  //case pion
959  if(abs(pdg) == 211) {
960  if(random_number > (1 - 0.9)) {
961  recopidecal.push_back(211);
962  }
963  else{
964  recopidecal.push_back(13);
965  }
966  }
967  }
968  else
969  {
970  //case muon
971  if(abs(pdg) == 13){
972  if(random_number > (1 - 0.95)) {
973  recopidecal.push_back(13);
974  }
975  else{
976  recopidecal.push_back(211);
977  }
978  }
979  //case pion
980  if(abs(pdg) == 211){
981  if(random_number > (1 - 0.95)) {
982  recopidecal.push_back(211);
983  }
984  else{
985  recopidecal.push_back(13);
986  }
987  }
988  }
989  }
990  else if( abs(pdg) == 2212 )
991  {
992  recopidecal.push_back(2212);//TODO for p/pi separation
993  }
994  else {
995  recopidecal.push_back(-1);
996  }
997  }
998  }
999  else if( fHelper->isThroughCalo(epoint) )
1000  {
1001  //Case the endpoint is outside the CALO -> it went through the ECAL (mu/pi/p possible)
1002  //the ECAL will see 60 MIPs on average
1003  double Evis = (double)nLayers; //in MIP
1004  //Smearing to account for Gaussian detector noise (Landau negligible)
1005  Evis = fHelper->GaussianSmearing(Evis, ECAL_MIP_Res);
1006  //1 MIP = 0.814 MeV
1007  double Erec = Evis * MIP2GeV_factor;
1008  erecon.push_back((Erec > 0) ? Erec : 0.);
1009  etime.push_back(ecaltime);
1010  detected.push_back(1);
1011 
1012  //Muon/Pions/Protons are reco as Muons (without MuID detector)
1013  if( abs(pdg) == 13 || abs(pdg) == 211 || abs(pdg) == 2212 ) {
1014  recopidecal.push_back(13);
1015  }
1016  else{
1017  recopidecal.push_back(-1);
1018  }
1019  }
1020  else
1021  {
1022  //Does not reach the ECAL???
1023  erecon.push_back(-1);
1024  recopidecal.push_back(-1);
1025  etime.push_back(-1);
1026  detected.push_back(0);
1027  }
1028  } //end endpoint is not in TPC
1029 
1031 
1032  } // end is charged mcp
1033  else
1034  {
1035  //not in the pdglist of particles but visible in TPC?
1036  auto found = std::find(pdg_charged.begin(), pdg_charged.end(), abs(pdg));
1037  if(found == pdg_charged.end())
1038  {
1039  detected.push_back(0);
1040  etime.push_back(-1);
1041  erecon.push_back(-1);
1042  _preco.push_back(-1);
1043  anglereco.push_back(-1);
1044  recopid.push_back(-1);
1045  recopidecal.push_back(-1);
1046  for (int pidr = 0; pidr < 6; ++pidr) prob_arr.push_back(-1);
1047  }
1048  }
1049 
1050  //end tpc
1051  //***************************************************************************************************************/
1052  }
std::vector< double > _preco
void DoGluckSternCalculation(float &preco, float &angle_reco)
std::vector< int > detected
const double ECAL_MIP_Res
const double MIP2GeV_factor
float GaussianSmearing(const float mean, const float sigma)
bool PointInTPC(const TVector3 &point)
std::vector< int > recopidecal
void DoRangeCalculation(float &preco, float &angle_reco)
std::vector< int > _MCPEndY
T abs(T value)
std::vector< int > truepdg
std::vector< int > _MCPEndZ
bool isThroughCalo(const TVector3 &point)
std::vector< double > prob_arr
std::vector< int > recopid
std::vector< double > truep
unsigned int _nFSP
CAFHelper * fHelper
void TPCParticleIdentification()
const std::vector< int > pdg_charged
bool PointInCalo(const TVector3 &point)
std::vector< double > anglereco
std::vector< int > _MCPEndX
std::vector< double > erecon
std::vector< double > etime

Member Data Documentation

std::vector<double> gar::ParamSim::_angle
private

Definition at line 348 of file ParamSim_module.cc.

std::vector<std::string> gar::ParamSim::_MCEndProc
private

Definition at line 347 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::_MCPEndX
private

Definition at line 346 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::_MCPEndY
private

Definition at line 346 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::_MCPEndZ
private

Definition at line 346 of file ParamSim_module.cc.

std::vector<std::string> gar::ParamSim::_MCProc
private

Definition at line 347 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::_MCPStartX
private

Definition at line 346 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::_MCPStartY
private

Definition at line 346 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::_MCPStartZ
private

Definition at line 346 of file ParamSim_module.cc.

unsigned int gar::ParamSim::_nFSP
private

Definition at line 345 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::_preco
private

Definition at line 351 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::anglereco
private

Definition at line 351 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::ccnc
private

Definition at line 342 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::detected
private

Definition at line 342 of file ParamSim_module.cc.

const float gar::ParamSim::ECAL_const = 0.02
private

Definition at line 323 of file ParamSim_module.cc.

const double gar::ParamSim::ECAL_MIP_Res = 0.23
private

Definition at line 330 of file ParamSim_module.cc.

const float gar::ParamSim::ECAL_stock = 0.06
private

Definition at line 322 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::erecon
private

Definition at line 351 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::etime
private

Definition at line 351 of file ParamSim_module.cc.

bool gar::ParamSim::fCorrect4origin
private

Definition at line 297 of file ParamSim_module.cc.

const float gar::ParamSim::fECALTimeResolution = 1.
private

Definition at line 334 of file ParamSim_module.cc.

CLHEP::HepRandomEngine& gar::ParamSim::fEngine
private

random engine

Definition at line 292 of file ParamSim_module.cc.

unsigned int gar::ParamSim::fEvent
private

Definition at line 340 of file ParamSim_module.cc.

std::string gar::ParamSim::fGeantLabel
private

module label for geant4 simulated hits

Definition at line 296 of file ParamSim_module.cc.

std::string gar::ParamSim::fGeneratorLabel
private

Definition at line 295 of file ParamSim_module.cc.

const geo::GeometryCore* gar::ParamSim::fGeo
private

pointer to the geometry

Definition at line 290 of file ParamSim_module.cc.

CAFHelper* gar::ParamSim::fHelper
private

Definition at line 291 of file ParamSim_module.cc.

float gar::ParamSim::fOrigin[3]
private

Definition at line 300 of file ParamSim_module.cc.

TF1* gar::ParamSim::fRes
private

Definition at line 324 of file ParamSim_module.cc.

unsigned int gar::ParamSim::fRun
private

Definition at line 340 of file ParamSim_module.cc.

unsigned int gar::ParamSim::fSubRun
private

Definition at line 340 of file ParamSim_module.cc.

TTree* gar::ParamSim::fTree
private

Definition at line 286 of file ParamSim_module.cc.

const float gar::ParamSim::gastpc_B = 0.5
private

Definition at line 307 of file ParamSim_module.cc.

const float gar::ParamSim::gastpc_len = 2.
private

Definition at line 306 of file ParamSim_module.cc.

const float gar::ParamSim::gastpc_padPitch = 0.1
private

Definition at line 308 of file ParamSim_module.cc.

const float gar::ParamSim::gastpc_X0 = 1300.
private

Definition at line 309 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::gint
private

Definition at line 342 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::gt_t
private

Definition at line 342 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::intert
private

Definition at line 342 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isBarrelEnd
private

Definition at line 355 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isBarrelStart
private

Definition at line 355 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isCaloEnd
private

Definition at line 354 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isCaloStart
private

Definition at line 353 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isEndcapEnd
private

Definition at line 355 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isEndcapStart
private

Definition at line 355 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isFidEnd
private

Definition at line 354 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isFidStart
private

Definition at line 353 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isInBetweenEnd
private

Definition at line 354 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isInBetweenStart
private

Definition at line 353 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isThroughCaloEnd
private

Definition at line 354 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isThroughCaloStart
private

Definition at line 353 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isTPCEnd
private

Definition at line 354 of file ParamSim_module.cc.

std::vector<unsigned int> gar::ParamSim::isTPCStart
private

Definition at line 353 of file ParamSim_module.cc.

std::unordered_map<int, TH2F*> gar::ParamSim::m_pidinterp
private

Definition at line 287 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::mcnupx
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::mcnupy
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::mcnupz
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::mctime
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::mctrkid
private

Definition at line 346 of file ParamSim_module.cc.

const double gar::ParamSim::MIP2GeV_factor = 0.814 / 1000
private

Definition at line 332 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::mode
private

Definition at line 342 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::motherid
private

Definition at line 346 of file ParamSim_module.cc.

const std::vector<int> gar::ParamSim::neutrinos = {12, 14, 16}
private

Definition at line 301 of file ParamSim_module.cc.

TParticlePDG* gar::ParamSim::neutron = TDatabasePDG::Instance()->GetParticle(2112)
private

Definition at line 335 of file ParamSim_module.cc.

const float gar::ParamSim::neutron_mass = neutron->Mass()
private

Definition at line 336 of file ParamSim_module.cc.

const float gar::ParamSim::NeutronECAL_detEff[2] = {0.2, 0.4}
private

Definition at line 317 of file ParamSim_module.cc.

const int gar::ParamSim::nLayers = 60
private

Definition at line 328 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::ntype
private

Definition at line 342 of file ParamSim_module.cc.

const std::vector<int> gar::ParamSim::pdg_charged = {211, 13, 2212, 321, 1000010020, 11}
private

Definition at line 305 of file ParamSim_module.cc.

const std::vector<int> gar::ParamSim::pdg_neutral = {22, 2112, 111, 130, 310, 311, 2114}
private

Definition at line 303 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::pdgmother
private

Definition at line 346 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::prob_arr
private

Definition at line 351 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::q2
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::recopid
private

Definition at line 350 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::recopidecal
private

Definition at line 350 of file ParamSim_module.cc.

const float gar::ParamSim::sigma_x = 0.1
private

Definition at line 313 of file ParamSim_module.cc.

const float gar::ParamSim::sigmaNeutronECAL_first = 0.11
private

Definition at line 318 of file ParamSim_module.cc.

const float gar::ParamSim::sigmaP_short = 0.1
private

Definition at line 311 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::t
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::tgtpdg
private

Definition at line 342 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::theta
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::trkLen
private

Definition at line 348 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::trkLenPerp
private

Definition at line 348 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::truep
private

Definition at line 348 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::truepdg
private

Definition at line 346 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::truepx
private

Definition at line 348 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::truepy
private

Definition at line 348 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::truepz
private

Definition at line 348 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::vertx
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::verty
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::vertz
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::w
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<int> gar::ParamSim::weight
private

Definition at line 342 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::x
private

Definition at line 343 of file ParamSim_module.cc.

std::vector<double> gar::ParamSim::y
private

Definition at line 343 of file ParamSim_module.cc.


The documentation for this class was generated from the following file: