All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
sbnd::NuAnaAlg Class Reference

#include <NuAnaAlg.h>

Public Member Functions

 NuAnaAlg ()
 
void configureGeometry (art::ServiceHandle< geo::Geometry >)
 
void configureReWeight (const std::vector< reweight > &, const std::vector< std::vector< float >> &)
 
void calcWeight (art::Ptr< simb::MCTruth >, art::Ptr< simb::GTruth >, std::vector< std::vector< float >> &)
 
void packFluxWeight (art::Ptr< simb::MCFlux > flux, std::vector< std::vector< float >> &)
 
unsigned int prepareSigmas (int, unsigned int, std::vector< std::vector< float > > &)
 
void parseWeights (const std::vector< std::string > &, std::vector< reweight > &)
 
void packNeutrinoInfo (simb::MCNeutrino *neutrino, int &nuchan, int &inno, double &enugen, int &isCC, int &mode, double &thetaLep, double &phiLep, double &Elep, std::vector< float > &neutMom, std::vector< float > &vertex)
 
void packFluxInfo (art::Ptr< simb::MCFlux > flux, int &ptype, int &tptype, int &ndecay, std::vector< float > &neutVertexInWindow, std::vector< float > &ParentVertex, std::vector< float > &nuParentMomAtDecay, std::vector< float > &nuParentMomAtProd, std::vector< float > &nuParentMomTargetExit)
 
void packGenieInfo (art::Ptr< simb::MCTruth > truth, std::vector< int > &GeniePDG, std::vector< std::vector< float >> &GenieMomentum, std::vector< std::string > &GenieProc, int &NPi0FinalState, int &NGamma, int &NChargedPions)
 
void packLarg4Info (art::Handle< std::vector< simb::MCParticle > > mclarg4, int, int, int, int, std::vector< std::vector< float > > &leptonPos, std::vector< std::vector< float > > &leptonMom, std::vector< std::vector< float > > &p1PhotonConversionPos, std::vector< std::vector< float > > &p1PhotonConversionMom, std::vector< std::vector< float > > &p2PhotonConversionPos, std::vector< std::vector< float > > &p2PhotonConversionMom, std::vector< std::vector< float > > &miscPhotonConversionPos, std::vector< std::vector< float > > &miscPhotonConversionMom, std::vector< std::vector< float > > &pionPos, std::vector< std::vector< float > > &pionMom, std::vector< std::vector< std::vector< float > > > &chargedPionPos, std::vector< std::vector< std::vector< float > > > &chargedPionMom, std::vector< int > &chargePionSign)
 

Private Member Functions

art::Ptr< simb::MCParticle > getParticleByID (art::Handle< std::vector< simb::MCParticle > > &mclistLARG4, int) const
 
art::Ptr< simb::MCParticle > getParticleByPDG (art::Handle< std::vector< simb::MCParticle > > &mclistLARG4, int) const
 
void pack4Vector (const TLorentzVector &input, std::vector< float > &output) const
 
void pack3Vector (const TVector3 &input, std::vector< float > &output) const
 
bool isInTPC (const TVector3 &) const
 
void GetPhotonConversionInfo (art::Ptr< simb::MCParticle > photon, TLorentzVector &ConversionPos, TLorentzVector &ConversionMom)
 

Private Attributes

std::vector< std::vector
< rwgt::NuReweight * > > 		reweightVector
 
double xlow
 
double xhigh
 
double ylow
 
double yhigh
 
double zlow
 
double zhigh
 

Detailed Description

Definition at line 30 of file NuAnaAlg.h.

Constructor & Destructor Documentation

sbnd::NuAnaAlg::NuAnaAlg ( )

Definition at line 6 of file NuAnaAlg.cxx.

6  {
7  }

Member Function Documentation

void sbnd::NuAnaAlg::calcWeight ( art::Ptr< simb::MCTruth >  mctruth,
art::Ptr< simb::GTruth >  gtruth,
std::vector< std::vector< float >> &  weights 
)

Definition at line 257 of file NuAnaAlg.cxx.

259  {
260  // return reweight.CalcWeight(*mctruth,*gtruth);
261  // if (weights.size() == 0) weights.resize(1);
262  // weights.front().push_back( reweight -> CalcWeight(*mctruth,*gtruth) );
263 
264  // weights needs to be the size of the reweighting vector
265  if (weights.size() != reweightVector.size())
266  weights.resize(reweightVector.size());
267 
268  for (unsigned int i_weight = 0; i_weight < reweightVector.size(); i_weight ++){
269  if (weights[i_weight].size() != reweightVector[i_weight].size()){
270  weights[i_weight].resize(reweightVector[i_weight].size());
271  }
272  for (unsigned int i_reweightingKnob = 0;
273  i_reweightingKnob < reweightVector[i_weight].size();
274  i_reweightingKnob ++)
275  {
276  weights[i_weight][i_reweightingKnob]
277  = reweightVector[i_weight][i_reweightingKnob]
278  -> CalcWeight(*mctruth,*gtruth);
279  }
280  }
281 
282  return;
283  }
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
sbnd::NuAnaAlg::reweightVector
std::vector< std::vector< rwgt::NuReweight * > > reweightVector
Definition: NuAnaAlg.h:120
void sbnd::NuAnaAlg::configureGeometry ( art::ServiceHandle< geo::Geometry >  geom)

Definition at line 14 of file NuAnaAlg.cxx.

14  {
15 
16  xlow = - geom -> DetHalfWidth();
17  xhigh = geom -> DetHalfWidth();
18  ylow = - geom -> DetHalfHeight();
19  yhigh = geom -> DetHalfHeight();
20  zlow = 0.0;
21  zhigh = geom -> DetLength();
22 
23  if (geom -> DetectorName() == "uboone_basic"){
24  xlow = 0.0;
25  xhigh = 2*xhigh;
26  }
27 
28  std::cout << "The dimensions of this detector are read to be:\n"
29  << " x: " << xlow << " to " << xhigh << "\n"
30  << " y: " << ylow << " to " << yhigh << "\n"
31  << " z: " << zlow << " to " << zhigh << "\n";
32 
33  return;
34  }
sbnd::NuAnaAlg::xhigh
double xhigh
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::xlow
double xlow
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::zhigh
double zhigh
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::ylow
double ylow
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::yhigh
double yhigh
Definition: NuAnaAlg.h:123
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
sbnd::NuAnaAlg::zlow
double zlow
Definition: NuAnaAlg.h:123
void sbnd::NuAnaAlg::configureReWeight ( const std::vector< reweight > &  weights,
const std::vector< std::vector< float >> &  reweightingSigmas 
)

This function configures the reweighting machinery. I am guess that, due to the overhead in configuring the genie reweighting, it's faster to make one reweight object for each weight needed. This will also make a "total" reweight object that will set all the switches on for all weighting parameters.

Definition at line 44 of file NuAnaAlg.cxx.

45  {
46 
47  // Expand the vector to the correct number of physical knobs plus 1
48  reweightVector.resize(weights.size()+1);
49 
50  // if (weights.size()+1 != reweightingSigmas.size()){
51  // std::cerr << "Error configuring the reweights, the number of weights must be"
52  // << " equal to the number of boundaries (range).\n";
53  // exit(-1);
54  // }
55 
56  // don't forget the last vector with all of the weights
57  reweightVector.back().resize(reweightingSigmas.front().size());
58  for (auto & ptr : reweightVector.back()) ptr = new rwgt::NuReweight;
59 
60  for (unsigned int i_weight = 0; i_weight < reweightingSigmas.front().size(); ++i_weight)
61  {
62  // reweightVector.back().at(i_weight)
63  // -> ReweightNCEL(reweightingSigmas[kNCEL][i_weight]);
64  reweightVector.back().at(i_weight)
65  -> ReweightQEMA(reweightingSigmas[kQEMA][i_weight]);
66  // reweightVector.back().at(i_weight)
67  // -> ReweightQEVec(reweightingSigmas[kQEVec][i_weight]);
68  reweightVector.back().at(i_weight)
69  -> ReweightResGanged(reweightingSigmas[kResGanged][i_weight]);
70  reweightVector.back().at(i_weight)
71  -> ReweightCCRes(reweightingSigmas[kCCRes][i_weight]);
72  reweightVector.back().at(i_weight)
73  -> ReweightNCRes(reweightingSigmas[kNCRes][i_weight]);
74  // reweightVector.back().at(i_weight)
75  // -> ReweightCoh(reweightingSigmas[kCoh][i_weight]);
76  reweightVector.back().at(i_weight)
77  -> ReweightNonResRvp1pi(reweightingSigmas[kNonResRvp1pi][i_weight]);
78  reweightVector.back().at(i_weight)
79  -> ReweightNonResRvbarp1pi(reweightingSigmas[kNonResRvbarp1pi][i_weight]);
80  reweightVector.back().at(i_weight)
81  -> ReweightNonResRvp2pi(reweightingSigmas[kNonResRvp2pi][i_weight]);
82  reweightVector.back().at(i_weight)
83  -> ReweightNonResRvbarp2pi(reweightingSigmas[kNonResRvbarp2pi][i_weight]);
84  // reweightVector.back().at(i_weight)
85  // -> ReweightResDecay(reweightingSigmas[kResDecay][i_weight]);
86  reweightVector.back().at(i_weight)
87  -> ReweightNC(reweightingSigmas[kNC][i_weight]);
88  // reweightVector.back().at(i_weight)
89  // -> ReweightDIS(reweightingSigmas[kDIS][i_weight]);
90  reweightVector.back().at(i_weight)
91  -> ReweightDISnucl(reweightingSigmas[kDISnucl][i_weight]);
92  // reweightVector.back().at(i_weight)
93  // -> ReweightAGKY(reweightingSigmas[kAGKY][i_weight]);
94  }
95 
96 
97  // loop over the physical knobs and expand to the correct number of weights
98  for (unsigned int i_reweightingKnob = 0;
99  i_reweightingKnob < reweightVector.size()-1;
100  i_reweightingKnob++)
101  {
102 
103  // resize this row to accomodate all of the weight points
104  reweightVector[i_reweightingKnob].resize(
105  reweightingSigmas[i_reweightingKnob].size());
106 
107  for (unsigned int weight_point = 0;
108  weight_point < reweightingSigmas[i_reweightingKnob].size();
109  weight_point++){
110 
111  // Figure out what is the value going in to this reweight
112  // double stepSize = (reweightingSigmas[i_reweightingKnob]
113  // - rangeLow[i_reweightingKnob])/(nWeights-1);
114  // double reweightingValue = rangeLow[i_reweightingKnob]
115  // + weight_point*stepSize;
116 
117  reweightVector[i_reweightingKnob][weight_point] = new rwgt::NuReweight;
118 
119  switch (weights[i_reweightingKnob]){
120  case kNCEL:
121  // reweightVector[i_reweightingKnob][weight_point]
122  // -> ReweightNCEL(reweightingSigmas[kNCEL][weight_point]);
123  break;
124  case kQEMA:
125  reweightVector[i_reweightingKnob][weight_point]
126  -> ReweightQEMA(reweightingSigmas[kQEMA][weight_point]);
127  break;
128  case kQEVec:
129  reweightVector[i_reweightingKnob][weight_point]
130  -> ReweightQEVec(reweightingSigmas[kQEVec][weight_point]);
131  break;
132  case kResGanged:
133  reweightVector[i_reweightingKnob][weight_point]
134  -> ReweightResGanged(reweightingSigmas[kResGanged][weight_point]);
135  break;
136  case kCCRes:
137  reweightVector[i_reweightingKnob][weight_point]
138  -> ReweightCCRes(reweightingSigmas[kCCRes][weight_point]);
139  break;
140  case kNCRes:
141  reweightVector[i_reweightingKnob][weight_point]
142  -> ReweightNCRes(reweightingSigmas[kNCRes][weight_point]);
143  break;
144  case kCoh:
145  // reweightVector[i_reweightingKnob][weight_point]
146  // -> ReweightCoh(reweightingSigmas[kCoh][weight_point]);
147  break;
148  case kNonResRvp1pi:
149  reweightVector[i_reweightingKnob][weight_point]
150  -> ReweightNonResRvp1pi(reweightingSigmas[kNonResRvp1pi][weight_point]);
151  break;
152  case kNonResRvbarp1pi:
153  reweightVector[i_reweightingKnob][weight_point]
154  -> ReweightNonResRvbarp1pi(reweightingSigmas[kNonResRvbarp1pi][weight_point]);
155  break;
156  case kNonResRvp2pi:
157  reweightVector[i_reweightingKnob][weight_point]
158  -> ReweightNonResRvp2pi(reweightingSigmas[kNonResRvp2pi][weight_point]);
159  break;
160  case kNonResRvbarp2pi:
161  reweightVector[i_reweightingKnob][weight_point]
162  -> ReweightNonResRvbarp2pi(reweightingSigmas[kNonResRvbarp2pi][weight_point]);
163  break;
164  case kResDecay:
165  // reweightVector[i_reweightingKnob][weight_point]
166  // -> ReweightResDecay(reweightingSigmas[kResDecay][weight_point]);
167  break;
168  case kNC:
169  reweightVector[i_reweightingKnob][weight_point]
170  -> ReweightNC(reweightingSigmas[kNC][weight_point]);
171  break;
172  case kDIS:
173  // reweightVector[i_reweightingKnob][weight_point]
174  // -> ReweightDIS(reweightingSigmas[kDIS][weight_point]);
175  break;
176  case kDISnucl:
177  reweightVector[i_reweightingKnob][weight_point]
178  -> ReweightDISnucl(reweightingSigmas[kDISnucl][weight_point]);
179  break;
180  case kAGKY:
181  // reweightVector[i_reweightingKnob][weight_point]
182  // -> ReweightAGKY(reweightingSigmas[kAGKY][weight_point]);
183  break;
184  case kNReWeights:
185  break;
186  }
187 
188  } //loop over nWeights
189  } // loop over physical knobs
190 
191 
192 
193 
194 
195 
196 
197  // std::cout << "\n\n\nsetup finished, running"
198  // << " configure on each weight......\n\n\n";
199 
200  // // Tell all of the reweight drivers to configure themselves:
201  // for(auto & vec : reweightVector){
202  // for (auto & driver : vec){
203  // driver -> Configure();
204  // }
205  // }
206 
207  return;
208 
209  }
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
sbnd::NuAnaAlg::reweightVector
std::vector< std::vector< rwgt::NuReweight * > > reweightVector
Definition: NuAnaAlg.h:120
art::Ptr< simb::MCParticle > sbnd::NuAnaAlg::getParticleByID ( art::Handle< std::vector< simb::MCParticle > > &  mclistLARG4,
int  TrackId 
) const
private

Definition at line 396 of file NuAnaAlg.cxx.

399  {
400  for(unsigned int i = 0; i < mclistLARG4 -> size(); i ++){
401  art::Ptr<simb::MCParticle> particle(mclistLARG4,i);
402  if ( particle -> TrackId() == TrackId){
403  return particle;
404  }
405  }
406  art::Ptr<simb::MCParticle> part;
407  return part;
408  }
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
art::Ptr< simb::MCParticle > sbnd::NuAnaAlg::getParticleByPDG ( art::Handle< std::vector< simb::MCParticle > > &  mclistLARG4,
int  PDG 
) const
private

Definition at line 409 of file NuAnaAlg.cxx.

412  {
413 
414  for(unsigned int i = 0; i < mclistLARG4 -> size(); i ++){
415  art::Ptr<simb::MCParticle> particle(mclistLARG4,i);
416  if ( particle -> PdgCode() == PDG){
417  return particle;
418  }
419  }
420  art::Ptr<simb::MCParticle> part;
421  return part;
422  }
PDG
process_name stream1 can override from command line with o or output services DetectorPropertiesService services DetectorPropertiesService services DetectorPropertiesService services DetectorPropertiesService physics analyzers pmtresponse NeutronTrackingCut services LArG4Parameters gaussian physics producers generator PDG
Definition: prodsingle_libraryana.fcl:81
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
void sbnd::NuAnaAlg::GetPhotonConversionInfo ( art::Ptr< simb::MCParticle >  photon,
TLorentzVector &  ConversionPos,
TLorentzVector &  ConversionMom 
)
private

Definition at line 432 of file NuAnaAlg.cxx.

434  {
435  // Just loop over the photons Trajectory points in momentum space
436  // and wait for a change. If it gets through the whole trajectory,
437  // the conversion point must be the end of the trajectory.
438  // When it changed, return the trajectory point immediately
439  // before to get the last point of the photon that's not converted.
440 
441  // Going to be watching for changes in energy,
442  // so we'd better know the start energy.
443  double E = photon->E(0);
444  // std::cout << "Starting Energy is " << E << std::endl;
445  // std::cout << "Looking at photon 4 vector, momentum currently starts \n\t("
446  // << photon->Momentum().X() << ", "
447  // << photon->Momentum().Y() << ", "
448  // << photon->Momentum().Z() << ", "
449  // << photon->Momentum().T() << ") ";
450  // std::cout << " At position \n\t( "
451  // << photon->Position().X() << ", "
452  // << photon->Position().Y() << ", "
453  // << photon->Position().Z() << ", "
454  // << photon->Position().T() << ") " << std::endl;
455  //Catch some special cases first.
456  if (photon->NumberTrajectoryPoints() == 0) return;
457  if (photon->NumberTrajectoryPoints() == 1){
458  ConversionPos = photon->EndPosition();
459  ConversionMom = photon->EndMomentum();
460  return;
461  }
462 
463  for(unsigned int i = 1; i < photon->NumberTrajectoryPoints(); i ++){
464  //ok, check if the energy has changed.
465  // std::cout << "Looking at photon 4 vector, momentum currently is \n\t("
466  // << photon->Momentum(i).X() << ", "
467  // << photon->Momentum(i).Y() << ", "
468  // << photon->Momentum(i).Z() << ", "
469  // << photon->Momentum(i).T() << ") ";
470  // std::cout << " At position \n\t("
471  // << photon->Position(i).X() << ", "
472  // << photon->Position(i).Y() << ", "
473  // << photon->Position(i).Z() << ", "
474  // << photon->Position(i).T() << ") " << std::endl;
475  if (photon->E(i) != E) {
476  //then the energy is different. Must be scattering or something.
477  //set the conversion points and bail!
478  ConversionPos = photon->Position(i);
479  ConversionMom = photon->Momentum(i-1); //<- This i OK since i starts at 1.
480  // std::cout << "Conversion Pos is " << ConversionPos << std::endl;
481  return;
482  }
483  }
484 
485  // If we made it out here, the photon must have ended without changing energy.
486  // So send back end position, momentum.
487  ConversionPos = photon->EndPosition();
488  ConversionMom = photon->EndMomentum();
489 
490  return;
491  }
photon
process_name can override from command line with o or output photon
Definition: runPID.fcl:28
E
process_name E
Definition: rawtorootconvert.fcl:26
bool sbnd::NuAnaAlg::isInTPC ( const TVector3 &  v) const
private

Definition at line 423 of file NuAnaAlg.cxx.

423  {
424  if (v.X() > xhigh || v.X() < xlow) return false;
425  if (v.Y() > yhigh || v.Y() < ylow) return false;
426  if (v.Z() > zhigh || v.Z() < zlow) return false;
427  return true;
428  }
sbnd::NuAnaAlg::xhigh
double xhigh
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::xlow
double xlow
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::zhigh
double zhigh
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::ylow
double ylow
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::yhigh
double yhigh
Definition: NuAnaAlg.h:123
sbnd::NuAnaAlg::zlow
double zlow
Definition: NuAnaAlg.h:123
void sbnd::NuAnaAlg::pack3Vector ( const TVector3 &  input,
std::vector< float > &  output 
) const
private

Definition at line 680 of file NuAnaAlg.cxx.

681  {
682  if (output.size() != 3) output.resize(3);
683  output[0] = input.X();
684  output[1] = input.Y();
685  output[2] = input.Z();
686  return;
687  }
void sbnd::NuAnaAlg::pack4Vector ( const TLorentzVector &  input,
std::vector< float > &  output 
) const
private

Definition at line 671 of file NuAnaAlg.cxx.

672  {
673  if (output.size() != 4) output.resize(4);
674  output[0] = input.E();
675  output[1] = input.X();
676  output[2] = input.Y();
677  output[3] = input.Z();
678  return;
679  }
void sbnd::NuAnaAlg::packFluxInfo ( art::Ptr< simb::MCFlux >  flux,
int &  ptype,
int &  tptype,
int &  ndecay,
std::vector< float > &  neutVertexInWindow,
std::vector< float > &  ParentVertex,
std::vector< float > &  nuParentMomAtDecay,
std::vector< float > &  nuParentMomAtProd,
std::vector< float > &  nuParentMomTargetExit 
)

Definition at line 333 of file NuAnaAlg.cxx.

339  {
340 
341  ptype = flux -> fptype;
342  tptype = flux->ftptype;
343  ndecay = flux->fndecay;
344 
345  neutVertexInWindow.clear();
346  ParentVertex.clear();
347  nuParentMomAtDecay.clear();
348  nuParentMomAtProd.clear();
349  nuParentMomTargetExit.clear();
350 
351  neutVertexInWindow.push_back(flux -> fgenx);
352  neutVertexInWindow.push_back(flux -> fgeny);
353  neutVertexInWindow.push_back(flux -> fgenz);
354  ParentVertex.push_back(flux -> fvx);
355  ParentVertex.push_back(flux -> fvy);
356  ParentVertex.push_back(flux -> fvz);
357  nuParentMomAtDecay.push_back(flux -> fpdpx);
358  nuParentMomAtDecay.push_back(flux -> fpdpy);
359  nuParentMomAtDecay.push_back(flux -> fpdpz);
360  nuParentMomAtProd.push_back(flux -> fppdxdz);
361  nuParentMomAtProd.push_back(flux -> fppdydz);
362  nuParentMomAtProd.push_back(flux -> fpppz);
363  nuParentMomTargetExit.push_back(flux -> ftpx);
364  nuParentMomTargetExit.push_back(flux -> ftpy);
365  nuParentMomTargetExit.push_back(flux -> ftpz);
366 
367  return;
368  }
void sbnd::NuAnaAlg::packFluxWeight ( art::Ptr< simb::MCFlux >  flux,
std::vector< std::vector< float >> &   
)

Definition at line 712 of file NuAnaAlg.cxx.

714  {
715 
716  std::cerr << "You are asking to pack flux weights but the version "
717  << "of nutools you use does not appear to support it.\n";
718  return;
719  }
cerr
BEGIN_PROLOG could also be cerr
Definition: messageservice.fcl:10
void sbnd::NuAnaAlg::packGenieInfo ( art::Ptr< simb::MCTruth >  truth,
std::vector< int > &  GeniePDG,
std::vector< std::vector< float >> &  GenieMomentum,
std::vector< std::string > &  GenieProc,
int &  NPi0FinalState,
int &  NGamma,
int &  NChargedPions 
)

Definition at line 371 of file NuAnaAlg.cxx.

377  {
378  int i = 0;
379  std::vector<float> tempMomentum;
380  tempMomentum.resize(4);
381  while( i < truth -> NParticles()){
382  auto part = truth -> GetParticle(i);
383  if (part.StatusCode() == 1){
384  GeniePDG.push_back(part.PdgCode());
385  GenieMomentum.push_back(tempMomentum);
386  pack4Vector(part.Momentum(),GenieMomentum.back());
387  GenieProc.push_back(part.Process());
388  if (part.PdgCode() == 22) NGamma ++;
389  if (part.PdgCode() == 111) NPi0FinalState ++;
390  if (abs(part.PdgCode()) == 211) NChargedPions ++;
391  }
392  i++;
393  }
394  }
sbnd::NuAnaAlg::pack4Vector
void pack4Vector(const TLorentzVector &input, std::vector< float > &output) const
Definition: NuAnaAlg.cxx:671
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
void sbnd::NuAnaAlg::packLarg4Info ( art::Handle< std::vector< simb::MCParticle > >  mclarg4,
int  ,
int  ,
int  ,
int  ,
std::vector< std::vector< float > > &  leptonPos,
std::vector< std::vector< float > > &  leptonMom,
std::vector< std::vector< float > > &  p1PhotonConversionPos,
std::vector< std::vector< float > > &  p1PhotonConversionMom,
std::vector< std::vector< float > > &  p2PhotonConversionPos,
std::vector< std::vector< float > > &  p2PhotonConversionMom,
std::vector< std::vector< float > > &  miscPhotonConversionPos,
std::vector< std::vector< float > > &  miscPhotonConversionMom,
std::vector< std::vector< float > > &  pionPos,
std::vector< std::vector< float > > &  pionMom,
std::vector< std::vector< std::vector< float > > > &  chargedPionPos,
std::vector< std::vector< std::vector< float > > > &  chargedPionMom,
std::vector< int > &  chargePionSign 
)

Definition at line 493 of file NuAnaAlg.cxx.

508  {
509 
510  // This function loops over the larg4 object and packs all the info
511  // up into the vectors provided (which are written to ntuple)
512  // We know from the genie function how many pi0 and gammas there are already
513  std::vector<float> empty4Vector;
514  empty4Vector.resize(4);
515 
516  // prepare the pi0 and photon vectors:
517  p1PhotonConversionPos.reserve(NPi0FinalState);
518  p1PhotonConversionMom.reserve(NPi0FinalState);
519  p2PhotonConversionPos.reserve(NPi0FinalState);
520  p2PhotonConversionMom.reserve(NPi0FinalState);
521  pionPos.reserve(NPi0FinalState);
522  pionMom.reserve(NPi0FinalState);
523  miscPhotonConversionPos.reserve(NGamma);
524  miscPhotonConversionMom.reserve(NGamma);
525 
526  // some counting variables to ensure everything is done correctly
527  int nPrimaryLepton(0); // should be == 1 at the end
528  int nPrimaryGamma(0); // should be == NGamma at the end
529  int nPrimaryPi0(0); // should be == NPi0FinalState
530 
531  chargedPionPos.resize(NChargedPions);
532  chargedPionMom.resize(NChargedPions);
533 
534 
535  // loop over the particles, extending the pion vectors as it goes.
536  for (unsigned int i = 0; i < mclarg4 -> size(); i ++){
537  art::Ptr<simb::MCParticle> particle(mclarg4,i);
538 
539 
540 
541  if (particle -> Mother() == 0 ){ // then this is a primary
542  // std::cout << "On particle " << particle -> TrackId()
543  // << " with PDG " << particle -> PdgCode() << std::endl;
544 
545  // For older files, set the nPrimaryLepton up since it didn't track neutrinos
546  if (isCC == 0) nPrimaryLepton ++;
547 
548  if (abs(particle -> PdgCode()) == 11 ||
549  abs(particle -> PdgCode()) == 12 ||
550  abs(particle -> PdgCode()) == 13 ||
551  abs(particle -> PdgCode()) == 14 )
552  {
553 
554  // then this is definitely the lepton.
555  nPrimaryLepton ++;
556  int nTrajectoryPoints = particle -> NumberTrajectoryPoints();
557  leptonPos.reserve(nTrajectoryPoints);
558  leptonMom.reserve(nTrajectoryPoints);
559  for (int traj_point = 0; traj_point < nTrajectoryPoints; ++traj_point){
560  leptonPos.push_back(empty4Vector);
561  pack4Vector(particle -> Position(traj_point),leptonPos.back());
562  leptonMom.push_back(empty4Vector);
563  pack4Vector(particle -> Momentum(traj_point),leptonMom.back());
564  if (!isInTPC(particle -> Position(traj_point).Vect())) break;
565  }
566  } // end lepton if block
567 
568  if (particle -> PdgCode() == 22) // misc gamma
569  {
570  nPrimaryGamma ++;
571  TLorentzVector conversionPoint, conversionMom;
572  GetPhotonConversionInfo(particle, conversionPoint,conversionMom);
573  miscPhotonConversionPos.push_back(empty4Vector);
574  miscPhotonConversionMom.push_back(empty4Vector);
575  pack4Vector(conversionPoint,miscPhotonConversionPos.back());
576  pack4Vector(conversionMom, miscPhotonConversionMom.back());
577  }
578 
579  if (particle -> PdgCode() == 111) // neutral pion
580  {
581  pionPos.push_back(empty4Vector);
582  pionMom.push_back(empty4Vector);
583  pack4Vector(particle->EndPosition(),pionPos.back());
584  pack4Vector(particle->Momentum(),pionMom.back());
585  nPrimaryPi0 ++;
586  // get the conversion point of each photon
587  if ( particle -> NumberDaughters() == 2){
588  // Get the info:
589  TLorentzVector conversionPoint, conversionMom;
590  art::Ptr<simb::MCParticle> daughter0
591  = getParticleByID(mclarg4,particle->Daughter(0));
592  GetPhotonConversionInfo(daughter0, conversionPoint,conversionMom);
593  // Pack it into the vectors
594  p1PhotonConversionPos.push_back(empty4Vector);
595  p1PhotonConversionMom.push_back(empty4Vector);
596  pack4Vector(conversionPoint,p1PhotonConversionPos.back());
597  pack4Vector(conversionMom, p1PhotonConversionMom.back());
598 
599  // get photon2:
600  art::Ptr<simb::MCParticle> daughter1
601  = getParticleByID(mclarg4,particle->Daughter(1));
602  GetPhotonConversionInfo(daughter1, conversionPoint,conversionMom);
603  // Pack it into the vectors:
604  p2PhotonConversionPos.push_back(empty4Vector);
605  p2PhotonConversionMom.push_back(empty4Vector);
606  pack4Vector(conversionPoint,p2PhotonConversionPos.back());
607  pack4Vector(conversionMom, p2PhotonConversionMom.back());
608  }
609  else{ // this is the "dalitz" decay
610  // Only take the photon, but be sure to find it
611  TLorentzVector conversionPoint, conversionMom;
612  p2PhotonConversionPos.push_back(empty4Vector);
613  p2PhotonConversionMom.push_back(empty4Vector);
614  for (int daughter = 0; daughter < 3; daughter ++ ){
615  art::Ptr<simb::MCParticle> daughterParticle
616  = getParticleByID(mclarg4,particle->Daughter(daughter));
617  if (daughterParticle->PdgCode() == 22) {
618  GetPhotonConversionInfo(daughterParticle,
619  conversionPoint,
620  conversionMom);
621  p1PhotonConversionPos.push_back(empty4Vector);
622  p1PhotonConversionMom.push_back(empty4Vector);
623  pack4Vector(conversionPoint,p1PhotonConversionPos.back());
624  pack4Vector(conversionMom, p1PhotonConversionMom.back());
625  break;
626  }
627  }
628  }
629 
630  }
631 
632  if (abs(particle -> PdgCode()) == 211) // charged pion
633  {
634  chargedPionSign.push_back(particle ->PdgCode() / 211);
635  int nTrajectoryPoints = particle -> NumberTrajectoryPoints();
636  unsigned int index = chargedPionSign.size() - 1;
637  chargedPionPos[index].reserve(nTrajectoryPoints);
638  chargedPionMom[index].reserve(nTrajectoryPoints);
639  for (int traj_point = 0; traj_point < nTrajectoryPoints; ++ traj_point){
640  chargedPionPos[index].push_back(empty4Vector);
641  pack4Vector(particle -> Position(traj_point),
642  chargedPionPos[index].back());
643  chargedPionMom[index].push_back(empty4Vector);
644  pack4Vector(particle -> Momentum(traj_point),
645  chargedPionMom[index].back());
646  }
647  }
648  } // end of if primary
649  } // end of loop over particles
650 
651  // if (nPrimaryLepton < 1){
652  // std::cerr << "Major problems here ... nPrimaryLepton\n"
653  // << " Should be " << 1
654  // << " but is " << nPrimaryLepton << "\n";
655  // exit(-1);
656  // }
657  if (nPrimaryGamma != NGamma){
658  std::cerr << "Major problems here ... nPrimaryGamma\n"
659  << " Should be " << NGamma
660  << " but is " << nPrimaryGamma << "\n";
661  exit(-1);
662  }
663  if (nPrimaryPi0 != NPi0FinalState){
664  std::cerr << "Major problems here ... nPrimaryPi0\n"
665  << " Should be " << NPi0FinalState
666  << " but is " << nPrimaryPi0 << "\n";
667  exit(-1);
668  }
669  } // end of packLarg4Info
sbnd::NuAnaAlg::pack4Vector
void pack4Vector(const TLorentzVector &input, std::vector< float > &output) const
Definition: NuAnaAlg.cxx:671
sbnd::NuAnaAlg::GetPhotonConversionInfo
void GetPhotonConversionInfo(art::Ptr< simb::MCParticle > photon, TLorentzVector &ConversionPos, TLorentzVector &ConversionMom)
Definition: NuAnaAlg.cxx:432
sbnd::NuAnaAlg::isInTPC
bool isInTPC(const TVector3 &) const
Definition: NuAnaAlg.cxx:423
cerr
BEGIN_PROLOG could also be cerr
Definition: messageservice.fcl:10
exit
exit
Definition: generate_icarus_geometry_splitwire.pl:430
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
sbnd::NuAnaAlg::getParticleByID
art::Ptr< simb::MCParticle > getParticleByID(art::Handle< std::vector< simb::MCParticle > > &mclistLARG4, int) const
Definition: NuAnaAlg.cxx:396
void sbnd::NuAnaAlg::packNeutrinoInfo ( simb::MCNeutrino *  neutrino,
int &  nuchan,
int &  inno,
double &  enugen,
int &  isCC,
int &  mode,
double &  thetaLep,
double &  phiLep,
double &  Elep,
std::vector< float > &  neutMom,
std::vector< float > &  vertex 
)

Definition at line 287 of file NuAnaAlg.cxx.

297  {
298  // Fill in all the variables that need to be filled!
299 
300  // prepare the vectors:
301  neutMom.clear();
302  neutMom.resize(4);
303  vertex.clear();
304  vertex.resize(3);
305 
306  // 1000 is the offset value used for NUANCE
307  nuchan = (neutrino->InteractionType() - 1000);
308  inno = neutrino->Nu().PdgCode();
309  enugen = neutrino->Nu().E();
310  // Is it a CC or NC interaction?
311  if (neutrino->CCNC()==simb::kCC) isCC = 1;
312  else if (neutrino->CCNC()==simb::kNC) isCC = 0;
313  mode = neutrino->Mode();
314 
315  simb::MCParticle lepton = neutrino->Lepton();
316  // thetaLep is the angle the lepton makes with the beam
317  double tempNumerator = lepton.Px()*lepton.Px() + lepton.Py()*lepton.Py();
318  thetaLep = atan(sqrt(tempNumerator)/lepton.Pz());
319  phiLep = atan(lepton.Px()/lepton.Py());
320  Elep = lepton.E();
321 
322  vertex[0] = neutrino->Nu().Position().X();
323  vertex[1] = neutrino->Nu().Position().Y();
324  vertex[2] = neutrino->Nu().Position().Z();
325  neutMom[0] = neutrino->Nu().Momentum().E();
326  neutMom[1] = neutrino->Nu().Momentum().X();
327  neutMom[2] = neutrino->Nu().Momentum().Y();
328  neutMom[3] = neutrino->Nu().Momentum().Z();
329 
330 
331  }
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:38
mode
const char mode
Definition: noise_ana.cxx:20
ana::Current::kNC
Neutral-current interactions.
Definition: IPrediction.h:39
void sbnd::NuAnaAlg::parseWeights ( const std::vector< std::string > &  string_weights,
std::vector< reweight > &  enum_weights 
)

Definition at line 229 of file NuAnaAlg.cxx.

230  {
231 
232  enum_weights.reserve(string_weights.size());
233  for( auto & s : string_weights){
234  if (s == "NCEL") enum_weights.push_back(kNCEL);
235  else if (s == "QEMA") enum_weights.push_back(kQEMA);
236  else if (s == "QEVec") enum_weights.push_back(kQEVec);
237  else if (s == "ResGanged") enum_weights.push_back(kResGanged);
238  else if (s == "CCRes") enum_weights.push_back(kCCRes);
239  else if (s == "NCRes") enum_weights.push_back(kNCRes);
240  else if (s == "Coh") enum_weights.push_back(kCoh);
241  else if (s == "NonResRvp1pi") enum_weights.push_back(kNonResRvp1pi);
242  else if (s == "NonResRvbarp1pi") enum_weights.push_back(kNonResRvbarp1pi);
243  else if (s == "NonResRvp2pi") enum_weights.push_back(kNonResRvp2pi);
244  else if (s == "NonResRvbarp2pi") enum_weights.push_back(kNonResRvbarp2pi);
245  else if (s == "ResDecay") enum_weights.push_back(kResDecay);
246  else if (s == "NC") enum_weights.push_back(kNC);
247  else if (s == "DIS") enum_weights.push_back(kDIS);
248  else if (s == "DISnucl") enum_weights.push_back(kDISnucl);
249  else if (s == "AGKY") enum_weights.push_back(kAGKY);
250  }
251 
252 
253 
254  }
s
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
unsigned int sbnd::NuAnaAlg::prepareSigmas ( int  NWeights,
unsigned int  RandSeed,
std::vector< std::vector< float > > &  reweightingSigmas 
)

Definition at line 211 of file NuAnaAlg.cxx.

214  {
215 
216  TRandom rand;
217  rand.SetSeed(RandSeed);
218 
219  reweightingSigmas.resize(kNReWeights);
220  for (unsigned int i = 0; i < reweightingSigmas.size(); ++i)
221  {
222  reweightingSigmas[i].resize(NWeights);
223  for (int j = 0; j < NWeights; j ++)
224  reweightingSigmas[i][j] = rand.Gaus(0,1);
225  }
226  return rand.GetSeed();
227  }

Member Data Documentation

std::vector<std::vector<rwgt::NuReweight *> > sbnd::NuAnaAlg::reweightVector
private

Definition at line 120 of file NuAnaAlg.h.

double sbnd::NuAnaAlg::xhigh
private

Definition at line 123 of file NuAnaAlg.h.

double sbnd::NuAnaAlg::xlow
private

Definition at line 123 of file NuAnaAlg.h.

double sbnd::NuAnaAlg::yhigh
private

Definition at line 123 of file NuAnaAlg.h.

double sbnd::NuAnaAlg::ylow
private

Definition at line 123 of file NuAnaAlg.h.

double sbnd::NuAnaAlg::zhigh
private

Definition at line 123 of file NuAnaAlg.h.

double sbnd::NuAnaAlg::zlow
private

Definition at line 123 of file NuAnaAlg.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5