127 produces< std::vector<cheat::SimHit> >();
128 produces< std::vector<cheat::SimTrack> >();
133 mf::LogInfo(
"MCCheater") <<
" MCCheater::MCCheater()\n";
136 mf::LogInfo(
"MCCheater") <<
" ****************************************\n" 137 <<
" * MCCheater Version Number 4 *\n" 138 <<
" ****************************************\n";
159 fMuE = tfs->
make<TH1F>(
"fMuE",
";True Muon Energy (GeV)", 50, 0., 10.);
160 fElE = tfs->
make<TH1F>(
"fElE",
";True Electron Energy (GeV)", 50, 0., 10.);
161 felL = tfs->
make<TH1F>(
"felL",
";length(cm)", 60, 0., 1200.);
162 fmuL = tfs->
make<TH1F>(
"fmuL",
";length(cm)", 60, 0., 1200.);
163 ftotalE = tfs->
make<TH1F>(
"ftotalE",
";E of SimHit (MeV)", 100, 0., 50.);
164 ftotEXV = tfs->
make<TH1F>(
"ftotEXV",
";E of Xview SimHit (MeV)", 100, 0., 50.);
165 ftotEYV = tfs->
make<TH1F>(
"ftotEYV",
";E of Yview SimHit (MeV)", 100, 0., 50.);
166 fdcosX = tfs->
make<TH1F>(
"fdcosX",
";Dir cos X", 100, -1., 1.);
167 fdcosY = tfs->
make<TH1F>(
"fdcosY",
";Dir cos Y", 100, -1., 1.);
168 fdcosZ = tfs->
make<TH1F>(
"fdcosZ",
";Dir cos Z", 100, -1., 1.);
169 fthreshE = tfs->
make<TH1F>(
"fthreshE",
";E of threshold SimHit (MeV)", 100, 0., 50.);
170 fsimhitE = tfs->
make<TH1F>(
"fsimhitE",
";Total SimHit Energy per Event (GeV)", 200, 0., 20.);
172 ftlen = tfs->
make<TH1F>(
"ftlen",
";Track Length (cm)", 250, 0., 500.);
174 fFLSE = tfs->
make<TH1F>(
"fFLSE",
";Total FLSHit Energy per Event (GeV)", 200, 0., 20.);
175 fnumsh = tfs->
make<TH1F>(
"fnumsh",
";Number of SimHits per Event", 50, -0.5, 999.5);
176 fnplanes = tfs->
make<TH1F>(
"fnplanes",
";Number of Planes in Track", 100, 0., 200.);
182 fNumFLS = tfs->
make<TH1F>(
"fNumFLS",
";Number of FLSHits in Event", 100, -0.5, 1999.5);
183 fnpart = tfs->
make<TH1F>(
"fnpart",
";Number of Particles in Event", 50, -0.5, 49.5);
184 fparte = tfs->
make<TH1F>(
"fparte",
";Total Energy of Particles in Event (GeV)", 200, 0., 20.);
185 fntrajpts = tfs->
make<TH1F>(
"fntrajpts",
";Number of Trajectory Points per Particle", 50, 0., 50.);
190 fcorx = tfs->
make<TH1F>(
"fcorx",
";X Coord of xview plane", 120, -810., 810.);
191 fcory = tfs->
make<TH1F>(
"fcory",
";Y Coord of yview plane", 120, -810., 810.);
192 fclhit = tfs->
make<TH1F>(
"fclhit",
";Penetration Depth with SimHits", 125, 0., 250.);
193 fzstart = tfs->
make<TH1F>(
"fzstart",
";Plane of particle start", 100, -1.0, 199.0);
204 mf::LogInfo(
"MCCheater") <<
"MCCheater::produce() " <<
" New Reco Event\n";
250 std::map<int, const sim::Particle*> trpartmap;
251 std::map<int, int> trpidmap;
253 unsigned int NPart = 0;
269 if (
abs(pid)==13)
fMuE->Fill(p.
E());
270 if (
abs(pid)==11)
fElE->Fill(p.
E());
273 if (
abs(pid)==12 ||
abs(pid)==14)
continue;
278 mf::LogInfo(
"MCCheater") <<
"Final particle count after cuts in particle list = " << pcount <<
"\n";
290 if (hitlist->empty()){
291 mf::LogError(
"MCCheater") <<
"No FLSHit Lists - bummer\n";
298 std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidpcpmap;
314 std::map<int, std::vector<const sim::FLSHit*> > trhitmap;
317 std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> > pcpflsmap;
319 for (
unsigned int i=0;
i<hitlist->size(); ++
i)
325 const std::vector<sim::FLSHit>&
hits = hlist->
fHits;
327 nFLShits += hits.size();
329 for (
unsigned int j=0;
j<hits.size(); ++
j)
331 int tid = hits[
j].GetTrackID();
350 if (smoth > -2) smoth = 1;
351 if (tmoth==-2)
mf::LogInfo(
"MCCheater") <<
" THIS TrackID NOT IN PARTICLE LIST = " << tid <<
" " << hits[
i].GetPDG() <<
"\n";
353 std::pair<int, int> mgmpair(tmoth, tgmoth);
355 unsigned int planeid = hits[
j].GetPlaneID();
356 unsigned int cellid = hits[
j].GetCellID();
358 std::pair<unsigned int, unsigned int> pcpair(planeid, cellid);
360 std::pair<unsigned int, unsigned int> tppair(tid, tpid);
363 xyzloc[0] = 0.5*(hits[
j].GetEntryX()+hits[
j].GetExitX());
364 xyzloc[1] = 0.5*(hits[
j].GetEntryY()+hits[
j].GetExitY());
365 xyzloc[2] = 0.5*(hits[
j].GetEntryZ()+hits[
j].GetExitZ());
369 trhitmap[tid].push_back(& hits[
j]);
370 pcpflsmap[pcpair].push_back(& hits[j]);
371 double end = hits[
j].GetEdep();
379 if (end>0.000000000001)
381 tidpcpmap[tid].push_back(pcpair);
393 fFLSE->Fill(flshitE);
422 std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidupcpmap;
423 for (std::map<
unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidpcpmap.begin(); itm != tidpcpmap.end(); ++itm)
425 unsigned int trackid = itm->first;
426 std::vector<std::pair<unsigned int, unsigned int> > pcarray = itm->second;
432 std::sort(pcarray.begin(), pcarray.end());
433 std::vector<std::pair<unsigned int, unsigned int> >::iterator endloc;
434 endloc = std::unique(pcarray.begin(), pcarray.end());
435 pcarray.erase(endloc, pcarray.end());
436 for (
unsigned int i=0;
i<pcarray.size(); ++
i)
438 tidupcpmap[trackid].push_back(pcarray[
i]);
446 unsigned int NumSH = 0;
447 std::vector<cheat::SimHit> simhits;
448 for (std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> >::iterator itm = pcpflsmap.begin(); itm != pcpflsmap.end(); ++itm)
452 std::pair<unsigned int, unsigned int> pcp = itm->first;
453 std::vector<const sim::FLSHit*> fh = itm->second;
454 unsigned int hsize = fh.size();
458 unsigned int planen = pcp.first;
460 unsigned int celln = pcp.second;
468 std::vector<double> earr;
469 std::vector<double> carr;
470 std::vector<int> tarr;
471 std::vector<int> parr;
473 std::vector<std::pair<int, int> > tparr;
474 std::vector<std::pair<int, double> > tearr;
475 std::vector<std::pair<int, double> > tcarr;
496 for (
unsigned int i = 0;
i<hsize; ++
i)
498 earr.push_back(fh[
i]->GetEdep());
499 std::pair<int, double> tep(fh[
i]->GetTrackID(), fh[
i]->GetEdep());
500 tearr.push_back(tep);
501 carr.push_back(fh[
i]->GetEntryT());
502 std::pair<int, double> ttp(fh[
i]->GetTrackID(), fh[
i]->GetEntryT());
503 tcarr.push_back(ttp);
505 tarr.push_back(fh[
i]->GetTrackID());
507 for (std::map<int, int>::iterator itm = trpidmap.begin(); itm != trpidmap.end(); ++itm)
509 int trid = itm->first;
510 if (trid == fh[
i]->GetTrackID())
512 parr.push_back(itm->second);
513 std::pair<int, int>
tpp(fh[
i]->GetTrackID(), itm->second);
514 tparr.push_back(tpp);
517 ene += fh[
i]->GetEdep();
525 double zexo = fh[0]->GetExitZ();
526 double zeno = fh[0]->GetEntryZ();
527 loc0[2] = 0.5*(zexo+zeno);
528 double xexo = fh[0]->GetExitX();
529 double xeno = fh[0]->GetEntryX();
530 loc0[0] = 0.5*(xexo+xeno);
531 double yexo = fh[0]->GetExitY();
532 double yeno = fh[0]->GetEntryY();
533 loc0[1] = 0.5*(yexo+yeno);
536 if (view==
geo::kX) coory = wor0[1];
537 if (view==
geo::kY) coorx = wor0[0];
542 hit.
SetPE(ene*5000.);
560 std::vector<int>::iterator startloc(tarr.begin());
561 std::vector<int>::iterator endloc(tarr.end());
562 int value(*startloc);
564 while (++startloc < endloc)
566 endloc = std::remove(startloc, endloc, value);
569 tarr.erase(endloc,tarr.end());
574 std::vector<int> mparr (1,parr[0]);
576 std::vector<double> mcarr (1,carr[0]);
579 for (
unsigned int i=0;
i<hsize; ++
i)
583 std::vector<double> mearr (1,esum);
590 if (hsize==tarr.size())
602 hit.
SetNP(tarr.size());
606 std::vector<int> pvec(tarr.size());
607 std::vector<double> evec(tarr.size());
608 std::vector<double> tvec(tarr.size());
609 for (
unsigned int i=0;
i<tarr.size(); ++
i)
611 for (
unsigned int j=0;
j<tparr.size(); ++
j)
613 double lowt = 9999999999.;
614 if (tparr[
j].first == tarr[
i])
616 pvec[
i] = tparr[
j].second;
617 evec[
i] += tearr[
j].second;
620 tvec[
i] = tcarr[
j].second;
621 lowt = tcarr[
j].second;
635 if( hit.
fTotE > 0.00000005)
638 simhits.push_back(hit);
646 std::cout <<
" Total E (GeV) in SimHits = " << Etotev <<
"\n";
649 std::cout <<
"Size of SimHit List = " << simhits.size() <<
"\n";
651 for (
unsigned int i=0;
i<simhits.size(); ++
i)
673 std::map<unsigned int, std::vector<cheat::SimHit*> > simtrhitmap;
674 std::vector<cheat::SimTrack> simtracks;
675 for (std::map<
unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidupcpmap.begin(); itm != tidupcpmap.end(); ++itm)
677 unsigned int tid = itm->first;
678 unsigned int upcsize = itm->second.size();
689 bool foundStart =
false;
690 unsigned int startpoint = 0;
692 bool foundEnd =
false;
693 unsigned int endpoint = 0;
694 std::vector<TLorentzVector> trpoints;
697 std::map<int, const sim::Particle*>::const_iterator
it;
698 for (it = trpartmap.begin(); it != trpartmap.end(); ++
it)
700 unsigned int trid = it->first;
708 unsigned int trajPoint = 0;
709 while(!foundStart && trajPoint<p->NumberTrajectoryPoints()){
712 scoord.SetXYZ(start.X(),start.Y(),start.Z());
714 startpoint = trajPoint;
722 endpoint = startpoint;
725 while(!foundEnd && trajPoint>=startpoint){
728 ecoord.SetXYZ(end.X(),end.Y(),end.Z());
730 endpoint = trajPoint;
738 if(foundStart && foundEnd && endpoint>startpoint){
739 for(
unsigned int trpoint = startpoint+1; trpoint <= endpoint;++trpoint){
742 trpoints.push_back(p->
Position(trpoint));
749 if(foundStart && foundEnd && endpoint>=startpoint){
756 for(
unsigned int trpoint = 0; trpoint < trpoints.size();++trpoint){
757 TVector3
pos(trpoints[trpoint].
X(),trpoints[trpoint].
Y(),trpoints[trpoint].
Z());
763 double dist =
sqrt((ecoord.X()-trpoints[trpoint].X())*(ecoord.X()-trpoints[trpoint].X())+
764 (ecoord.Y()-trpoints[trpoint].Y())*(ecoord.Y()-trpoints[trpoint].Y())+
765 (ecoord.Z()-trpoints[trpoint].Z())*(ecoord.Z()-trpoints[trpoint].Z()));
766 TVector3
dir((ecoord.X()-trpoints[trpoint].X())/dist,(ecoord.Y()-trpoints[trpoint].Y())/dist,(ecoord.Z()-trpoints[trpoint].Z())/dist);
772 simtracks.push_back(strack);
777 for (
unsigned int i=0;
i<simhits.size(); ++
i)
779 bool ismatch =
false;
782 for (
unsigned int j=0;
j<h.
fNP; ++
j)
784 if (h.
fTrID[
j]==(
int) tid) ismatch =
true;
786 if (ismatch) simtrhitmap[tid].push_back(& h);
792 std::unique_ptr<std::vector<cheat::SimTrack> > simtrackcol(
new std::vector<cheat::SimTrack>(simtracks) );
793 evt.
put(std::move(simtrackcol));
923 std::unique_ptr<std::vector<cheat::SimHit> > simhitcol(
new std::vector<cheat::SimHit>(simhits) );
924 std::cout <<
" Wrote SimHits to event of size = " << simhits.size() <<
"\n";
925 evt.
put(std::move(simhitcol));
double E(const int i=0) const
unsigned int NumberTrajectoryPoints() const
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
virtual void SetStart(TVector3 start)
std::vector< sim::FLSHit > fHits
void produce(art::Event &evt)
void GetCenter(double *xyz, double localz=0.0) const
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void LocalToWorld(const double *local, double *world) const
const CellGeo * Cell(int icell) const
list_type::const_iterator const_iterator
Vertical planes which measure X.
void SetChannel(uint32_t iChan)
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void SetTotE(double tote)
void SetTrackID(int trackid)
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
void SetView(geo::View_t view)
void SetPlane(unsigned short plane)
ProductID put(std::unique_ptr< PROD > &&product)
void SetCell(unsigned short cell)
View_t View() const
Which coordinate does this plane measure.
MCCheater(fhicl::ParameterSet const &pset)
const XML_Char int const XML_Char * value
void SetNP(unsigned short np)
void SetCoorZ(double coorz)
virtual void SetDir(TVector3 dir)
void AppendTrajectoryPoint(TVector3 pt)
void SetPID(std::vector< int > PID)
void SetTim(std::vector< double > Tim)
std::string fMCTruthListLabel
double DetHalfHeight() const
void SetTrID(std::vector< int > TrID)
void SetCoorX(double coorx)
void SetEnD(std::vector< double > EnD)
T * make(ARGS...args) const
code to link reconstructed objects back to the MC truth information
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::string fFLSHitListLabel
double fTotE
Total Energy dep from all particles.
void SetCoorY(double coory)
unsigned short fNP
size of vectors
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< int > fTrID
vector of sim track IDs
int EveId(const int trackID) const