10 #include "TDatabasePDG.h" 13 #include "TParticle.h" 14 #include "TPolyMarker.h" 15 #include "TPolyMarker3D.h" 16 #include "TPolyLine.h" 17 #include "TPolyLine3D.h" 44 void writeErrMsg(
const char*
fcn,
48 <<
" failed with message:\n" 72 if ( (drawopt->
fText &
79 std::vector<simb::MCTruth> mctruth;
85 const unsigned int mctextmax = 2000;
86 bool cappedtext =
false;
90 for (
unsigned int i = 0;
i < mctruth.size(); ++
i) {
91 if(
i > 0) mctext +=
", ";
100 for (
int j=0;
j<mctruth[
i].NParticles(); ++
j) {
106 const unsigned int bufsize = 1024;
109 snprintf(buff, bufsize,
"#color[%d]{%s #scale[0.75]{[%.1f#scale[0.5]{GeV/c}]}}",
115 snprintf(buff, bufsize,
"#color[%d]{%s}",
120 if(mctext.size() + incoming.size() + 3 + strlen(buff) > mctextmax){
124 if (!firstin) incoming +=
" + ";
129 if(mctext.size() + outgoing.size() + 3 + strlen(buff) > mctextmax){
133 if (!firstout) outgoing +=
" + ";
138 if (origin==
"" && incoming==
"") {
139 if(mctext.size() + outgoing.size() > mctextmax){
146 const int intType = mctruth[
i].GetNeutrino().InteractionType();
148 const std::string newtext = origin+incoming+
" #rightarrow "+outgoing+suffix;
149 if(mctext.size() + newtext.size() > mctextmax){
155 if(cappedtext)
break;
159 if(mctext.empty()) mctext =
"Too many particles";
160 else mctext +=
"... and more!";
165 latex.SetTextSize(0.15);
185 const TDatabasePDG*
db = TDatabasePDG::Instance();
186 const TParticlePDG* defn = db->GetParticle(part->
PdgCode());
188 ss << defn->GetName();
190 else if(part->
PdgCode() > 1000000000){
191 const int strange = (part->
PdgCode()/10000000)%10;
192 const int element = (part->
PdgCode()/10000)%1000;
193 const int A = (part->
PdgCode()/10)%1000;
194 const int I = part->
PdgCode()%10;
196 case 1: ss <<
"H";
break;
197 case 2: ss <<
"He";
break;
198 case 3: ss <<
"Li";
break;
199 case 4: ss <<
"Be";
break;
200 case 5: ss <<
"B";
break;
201 case 6: ss <<
"C";
break;
202 case 7: ss <<
"N";
break;
203 case 8: ss <<
"O";
break;
204 case 9: ss <<
"F";
break;
205 case 10: ss <<
"Ne";
break;
206 case 11: ss <<
"Na";
break;
207 case 12: ss <<
"Mg";
break;
208 case 13: ss <<
"Al";
break;
209 case 14: ss <<
"Si";
break;
210 case 15: ss <<
"P";
break;
211 case 16: ss <<
"S";
break;
212 case 17: ss <<
"Cl";
break;
213 case 18: ss <<
"Ar";
break;
214 case 19: ss <<
"K";
break;
215 case 20: ss <<
"Ca";
break;
216 case 21: ss <<
"Sc";
break;
217 case 22: ss <<
"Ti";
break;
218 case 23: ss <<
"V";
break;
219 case 24: ss <<
"Cr";
break;
220 case 25: ss <<
"Mn";
break;
221 case 26: ss <<
"Fe";
break;
222 case 27: ss <<
"Co";
break;
223 case 49: ss <<
"In";
break;
224 case 50: ss <<
"Sn";
break;
225 case 51: ss <<
"Sb";
break;
226 case 52: ss <<
"Te";
break;
227 case 53: ss <<
"I";
break;
228 case 54: ss <<
"Xe";
break;
229 case 55: ss <<
"Cs";
break;
230 case 56: ss <<
"Ba";
break;
231 case 57: ss <<
"La";
break;
232 default: ss << element <<
"-";
239 if(I != 0) ss <<
"*";
241 const char* Lambda =
"\u039b";
242 for(
int i = 0;
i < strange;
i++) ss << Lambda;
250 const double ke = part->
E()-part->
Mass();
253 sprintf(Estr,
"%.3lf GeV", ke);
255 sprintf(Estr,
"%.3lf MeV", 1e3*ke);
257 sprintf(Estr,
"%.3lf keV", 1e6*ke);
259 ss <<
"\t(" << Estr <<
")\t[" << part->
Process() <<
"]";
261 ss <<
" at (" << part->
Vx() <<
", " << part->
Vy() <<
", " << part->
Vz() <<
"; " << part->
T() <<
") ";
264 double ptot=part->
P();
265 if (showVtxs) ss<<
",";
267 ss <<
" p-hat (" << part->
Px()/ptot <<
", " << part->
Py()/ptot <<
", " << part->
Pz()/ptot <<
") ";
275 if(depth >= depthLimit && depthLimit > 0)
return;
277 const char* kHoriz =
"\u2500";
278 const char* kVert =
"\u2502";
279 const char* kCorner =
"\u2514";
280 const char* kTee =
"\u251c";
284 for(
int n = 0;
n < N; ++
n){
289 if(childIt == pnav.
end())
continue;
292 ss << prefix << kCorner << kHoriz << kHoriz;
294 prefix+
" ", depth+1);
297 ss << prefix << kTee << kHoriz << kHoriz;
299 prefix+kVert+
" ", depth+1);
312 if ((drawopt->
fText &
323 std::stringstream
ss;
324 ss<<
"Long MCTruth text ( from sim::Particle )"<<
std::endl;
325 ss<<
" Particle ( KineticEnergy ) [ Process ] ";
326 if (showVtxs) ss<<
" at (Vx, Vy, Vz, T)";
327 if (showDirs) ss<<
" with direction (Px, Py, Pz)";
332 ss <<
"\u2500\u2500";
345 if ( (drawopt->
fDraw &
349 std::vector<double>
x;
350 std::vector<double>
y;
351 std::vector<double>
z;
357 for (i=0; i<x.size(); ++
i) {
358 if(xzview) xzview->
AddMarker(z[i], x[i], gold, star, sz);
359 if(yzview) yzview->
AddMarker(z[i], y[i], gold, star, sz);
369 if ( (drawopt->
fDraw &
373 std::vector<double>
x;
374 std::vector<double>
y;
375 std::vector<double>
z;
384 for (i=0; i<x.size(); ++
i) {
385 pm.SetPoint(i, x[i], y[i], z[i]);
396 if ( (drawopt->
fDraw &
402 std::vector<simb::MCTruth> mctruth;
406 for (
unsigned int i=0;
i<mctruth.size(); ++
i) {
407 for (
int j=0;
j<mctruth[
i].NParticles(); ++
j) {
413 double r = p.
P()*100.0;
420 double x1 = x0 + r * p.
Px()/p.
P();
421 double y1 = y0 + r * p.
Py()/p.
P();
422 double z1 = z0 + r * p.
Pz()/p.
P();
424 TLine&
l = xzview->
AddLine(z0, x0, z1, x1);
428 TLine&
l = yzview->
AddLine(z0, y0, z1, y1);
441 if ( (drawopt->
fDraw &
447 std::vector<simb::MCTruth> mctruth;
451 for (
unsigned int i=0;
i<mctruth.size(); ++
i) {
452 for (
int j=0;
j<mctruth[
i].NParticles(); ++
j) {
458 double r = p.
P()*100.0;
465 double x1 = x0 + r * p.
Px()/p.
P();
466 double y1 = y0 + r * p.
Py()/p.
P();
467 double z1 = z0 + r * p.
Pz()/p.
P();
473 l.SetPoint(0, x0, y0, z0);
474 l.SetPoint(1, x1, y1, z1);
487 if ( (drawopt->
fDraw &
493 std::vector<simb::MCTruth> mctruth;
508 if (
abs(pdg) == 12 ||
abs(pdg) == 14 ||
abs(pdg) == 16)
continue;
516 bool draw_neutrals = (drawopt->
fDraw &
519 bool draw_gammas = (drawopt->
fDraw &
521 if (!draw_neutrals && s==kDotted)
continue;
522 if (!draw_gammas && s==kDashed)
continue;
524 TPolyLine* xz = xzview ? &xzview->
AddPolyLine(N, c, w, s) : 0;
525 TPolyLine* yz = yzview ? &yzview->
AddPolyLine(N, c, w, s) : 0;
526 TPolyLine3D* xyz = view3D ? &view3D->
AddPolyLine3D(N, c, w, s) : 0;
528 for(
int n = 0;
n < N; ++
n){
530 if(xz) xz->SetPoint(
n, vec.Z(), vec.X());
531 if(yz) yz->SetPoint(
n, vec.Z(), vec.Y());
532 if(xyz) xyz->SetPoint(
n, vec.X(), vec.Y(), vec.Z());
535 if(xz) xz->SetBit(kCannotPick);
536 if(yz) yz->SetBit(kCannotPick);
537 if(xyz) xyz->SetBit(kCannotPick);
562 if ( (drawopt->
fDraw &
573 std::vector<sim::FLSHitList> flshits;
577 std::map<int, TPolyMarker*> pms_x, pms_y;
580 for (
unsigned int i=0;
i<flshits.size(); ++
i) {
581 for (
unsigned int j=0;
j<flshits[
i].fHits.size(); ++
j) {
598 l.SetPoint(
i, xyzwor[2], xyzwor[
view]);
608 const double sz = 0.5;
611 TPolyMarker* pm = pms_x[
c];
616 pm->SetNextPoint(xyzwor[2], xyzwor[0]);
619 TPolyMarker* pm = pms_y[
c];
624 pm->SetNextPoint(xyzwor[2], xyzwor[1]);
628 if(xzview) xzview->
AddMarker(xyzwor[2],xyzwor[0],c,kOpenCircle,sz);
629 if(yzview) yzview->
AddMarker(xyzwor[2],xyzwor[1],c,kOpenCircle,sz);
644 if ( (drawopt->
fDraw &
657 std::vector<sim::FLSHitList> flshits;
659 for (
unsigned int i=0;
i<flshits.size(); ++
i) {
660 for (
unsigned int j=0;
j<flshits[
i].fHits.size(); ++
j) {
678 p.SetPoint(0, xyzwor[0], xyzwor[1], xyzwor[2]);
682 p.SetPoint(0, xyzwor[0], xyzwor[1], xyzwor[2]);
694 std::vector<simb::MCTruth>& mctruth)
702 std::vector<const simb::MCTruth*> truths;
708 return mctruth.size();
715 std::vector<sim::FLSHitList>& flshits)
723 std::vector<sim::FLSHitList>
temp(fcol->size());
724 for(
size_t f = 0;
f < fcol->size(); ++
f){
725 flshits.push_back(fcol->at(
f));
729 writeErrMsg(
"GetFLSHits", e);
732 return flshits.size();
738 std::vector<double>&
x,
739 std::vector<double>&
y,
740 std::vector<double>&
z)
750 std::vector<simb::MCTruth> mctruth;
752 for (
unsigned int i=0;
i<mctruth.size(); ++
i) {
755 for (
int j=0;
j<mctruth[
i].NParticles(); ++
j) {
768 for (
unsigned int k=0; k<x.size(); ++k) {
769 in_list = ((x0 == x[k]) && (y0 == y[k]) && (z0 == z[k]));
772 if (in_list==
false) {
843 double& zmin,
double&
zmax,
844 const std::set<geo::OfflineChan>& hmap)
852 xmin = +1e10; ymin = +1e10; zmin = +1e10;
853 xmax = -1e10; ymax = -1e10; zmax = -1e10;
856 std::vector<sim::FLSHitList> lists;
859 for(
unsigned int i = 0;
i < lists.size(); ++
i){
860 const std::vector<sim::FLSHit>&
hits = lists[
i].fHits;
862 for(
unsigned int j = 0;
j < hits.size(); ++
j){
867 hit.
GetCellID())) == hmap.end())
continue;
873 if(view ==
geo::kX && xyz[0] < xmin) xmin = xyz[0];
874 if(view ==
geo::kX && xyz[0] > xmax) xmax = xyz[0];
875 if(view ==
geo::kY && xyz[1] < ymin) ymin = xyz[1];
876 if(view ==
geo::kY && xyz[1] > ymax) ymax = xyz[1];
877 if(xyz[2] < zmin) zmin = xyz[2];
878 if(xyz[2] > zmax) zmax = xyz[2];
887 double& tmin,
double& tmax)
897 std::vector<sim::FLSHitList> lists;
900 for(
unsigned int i = 0;
i < lists.size(); ++
i){
901 const std::vector<sim::FLSHit>&
hits = lists[
i].fHits;
903 for(
unsigned int j = 0;
j < hits.size(); ++
j){
double E(const int i=0) const
neutral current quasi-elastic
static const int kTEXT_SHORT
float GetEntryX() const
Entry point of the particle (position, time and energy)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
unsigned int NumberTrajectoryPoints() const
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
resonant charged current, nubar p -> nubar n pi+
double Py(const int i=0) const
void MCTruthLongText(const art::Event &evt, evdb::View2D *view)
std::map< std::string, double > xmax
void FLSHit2D(const art::Event &evt, evdb::View2D *xzview, evdb::View2D *yzview)
resonant neutral current, nu p -> nu p pi0
charged current deep inelastic scatter
int fDraw
Which MC Truth objects to draw.
int GetPlaneID() const
Plane ID.
void HiLite(int trkId, bool hlt=true)
resonant charged current, nubar p -> l+ p pi-
Float_t y1[n_points_granero]
TPolyLine & AddPolyLine(int n, int c, int w, int s)
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
static const int kTEXT_LONG
Float_t x1[n_points_granero]
int GetCellID() const
Cell ID.
int GetNPoints() const
Number of end points in FLSHit (could be Geant end points)
TLine & AddLine(double x1, double y1, double x2, double y2)
resonant charged current, nubar p -> l+ n pi0
list_type::const_iterator const_iterator
Vertical planes which measure X.
double Px(const int i=0) const
::xsd::cxx::tree::exception< char > exception
A collection of drawable 2-D objects.
A single unit of energy deposition in the liquid scintillator.
static const int kDRAW_NEUTRALS
charged current quasi-elastic
int fFLSHitStyle
How to render FLS hits.
const PlaneGeo * Plane(unsigned int i) const
size_t GetProducts(const art::Event &evt, const std::string &which, const std::string &instance, std::vector< const T * > &prods)
std::string Process() const
void FromPDG(TAttLine *line, int pdgcode)
float GetExitX() const
Exit point of the particle (position, time and energy)
void MCTruthVectors2D(const art::Event &evt, evdb::View2D *xzview, evdb::View2D *yzview)
Global drawing options that apply to all displays.
int NumberDaughters() const
void GetTimeLimits(const art::Event *evt, double &tmin, double &tmax)
static const int kDRAW_HITS
Horizontal planes which measure Y.
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
resonant charged current, nubar p -> nubar p pi0
int LineWidthFromPDG(int pdgcode)
static const int kFLSHIT_AS_LINE
int Daughter(const int i) const
resonant charged current, nu n -> l- p pi0
TPolyLine3D & AddPolyLine3D(int n, int c, int w, int s)
int fText
What text to draw?
void PrintParticleAndOffspring(std::stringstream &ss, const sim::Particle *part, std::string prefix, int depth) const
Helper function for MCTruthLongText.
const sim::Particle * Primary(const int) const
resonant neutral current, nu n -> nu n pi0
std::vector< std::string > fMCTruthModules
Modules to load the data from.
resonant charged current, nu n -> l- n pi+
void MCTruthTrajectories2D(const art::Event &evt, evdb::View2D *xzview, evdb::View2D *yzview)
resonant charged current, nubar n -> nubar p pi-
float GetY(const int ipoint) const
float GetXAverage(const int step) const
Get X-average for the step. This is in local coordinates.
double P(const int i=0) const
TPolyMarker & AddPolyMarker(int n, int c, int st, double sz)
int ColorFromPDG(int pdgcode)
int fTextIncludeDirections
double T(const int i=0) const
void MCTruthVertices2D(const art::Event &evt, evdb::View2D *xzview, evdb::View2D *yzview)
int GetVertexPoints(const art::Event &evt, std::vector< double > &x, std::vector< double > &y, std::vector< double > &z)
TLatex & AddLatex(double x, double y, const char *text)
float fFLSHitThresh
Threshold to apply to FLS hits (MeV)
void GetLimits(const art::Event *evt, double &xmin, double &xmax, double &ymin, double &ymax, double &zmin, double &zmax, const std::set< geo::OfflineChan > &hmap={})
void MCTruthShortText(const art::Event &evt, evdb::View2D *view)
charged current deep inelastic scatter
resonant charged current, nu p -> l- p pi+
void MCTruthTrajectories3D(const art::Event &evt, evdb::View3D *view)
void MCTruthVectors3D(const art::Event &evt, evdb::View3D *view)
double Vx(const int i=0) const
int GetMCTruth(const art::Event &evt, std::vector< simb::MCTruth > &mctruth)
std::string ShortInteractionSuffix(int iType) const
charged current coherent pion
resonant neutral current, nu n -> nu p pi-
int GetTrackID() const
Track ID.
static const int kDRAW_VERTEX
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
float GetZ(const int ipoint) const
static const char * LatexName(int pdgcode)
Convert PDG code to a latex string (root-style)
void MCTruthVertices3D(const art::Event &evt, evdb::View3D *view)
double Pz(const int i=0) const
Render the objects from the Simulation package.
resonant charged current, nubar n -> l+ n pi-
void MCTruthTrajectoriesAnyD(const art::Event &evt, evdb::View2D *xzview, evdb::View2D *yzview, evdb::View3D *view3D)
Helper function for MCTruthTrajectories[2|3]D.
double Vz(const int i=0) const
assert(nhit_max >=nhit_nbins)
resonant charged current, nubar n -> nubar n pi0
int LineStyleFromPDG(int pdgcode)
A collection of 3D drawable objects.
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
std::vector< std::string > fFLSHitListModules
FLSHitLists here.
Event generator information.
float GetZAverage(const int step) const
Get Z-average for the step. This is in local coordinates.
float GetYAverage(const int step) const
Get Y-average for the step. This is in local coordinates.
float GetX(const int ipoint) const
Get a point of the particle in the FLSHit (position, time and energy)
TMarker & AddMarker(double x, double y, int c, int st, double sz)
static const int kDRAW_GAMMAS
Simple object representing a (plane, cell) pair.
iterator find(const key_type &key)
int GetFLSHits(const art::Event &evt, const char *which, std::vector< sim::FLSHitList > &flshits)
std::map< int, bool > fHighlite
static const int kDRAW_TRAJECT
void FLSHit3D(const art::Event &evt, evdb::View3D *xzview)
static const int kDRAW_VECTORS
double Vy(const int i=0) const
resonant neutral current, nu p -> nu p pi+
Encapsulate the geometry of one entire detector (near, far, ndos)
TPolyMarker3D & AddPolyMarker3D(int n, int c, int st, double sz)
int NumberOfPrimaries() const