43 #include "TProfile2D.h" 151 const double X0 = 14;
230 fLongitudinal = tfs->
make<TH3F>(
"fLongitudinal",
"longitudinal e- profile;t;electron energy (MeV);Q",
LBINS,
LMIN,
LMAX,
EBINS,
EMIN,
EMAX, 100, 0, 150000);
231 fTransverse = tfs->
make<TH3F>(
"fTransverse",
"transverse e- profile;dist (cm);electron energy (MeV);Q",
TBINS,
TMIN,
TMAX,
EBINS,
EMIN,
EMAX, 100, 0, 150000);
232 fTransverse_1 = tfs->
make<TH3F>(
"fTransverse_1",
"transverse e- profile [0 <= t < 1];dist (cm);electron energy (MeV);Q",
TBINS,
TMIN,
TMAX,
EBINS,
EMIN,
EMAX, 100, 0, 100000);
233 fTransverse_2 = tfs->
make<TH3F>(
"fTransverse_2",
"transverse e- profile [1 <= t < 2];dist (cm);electron energy (MeV);Q",
TBINS,
TMIN,
TMAX,
EBINS,
EMIN,
EMAX, 100, 0, 100000);
234 fTransverse_3 = tfs->
make<TH3F>(
"fTransverse_3",
"transverse e- profile [2 <= t < 3];dist (cm);electron energy (MeV);Q",
TBINS,
TMIN,
TMAX,
EBINS,
EMIN,
EMAX, 100, 0, 100000);
235 fTransverse_4 = tfs->
make<TH3F>(
"fTransverse_4",
"transverse e- profile [3 <= t < 4];dist (cm);electron energy (MeV);Q",
TBINS,
TMIN,
TMAX,
EBINS,
EMIN,
EMAX, 100, 0, 100000);
236 fTransverse_5 = tfs->
make<TH3F>(
"fTransverse_5",
"transverse e- profile [4 <= t < 5];dist (cm);electron energy (MeV);Q",
TBINS,
TMIN,
TMAX,
EBINS,
EMIN,
EMAX, 100, 0, 100000);
270 std::vector<art::Ptr<recob::Hit> > hitlist;
275 std::vector<art::Ptr<sim::SimChannel> > simchanlist;
280 std::vector<art::Ptr<recob::Shower> > showerlist;
285 std::vector<art::Ptr<simb::MCTruth> > mclist;
296 if (showerlist.size()) {
297 std::vector< art::Ptr<recob::Hit> > showerhits = shwfm.at(0);
324 auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
329 double shwVtxTime = detprop->ConvertXToTicks(shwvtx[0], collectionPlane);
330 double shwVtxWire = geom->WireCoordinate(shwvtx[1], shwvtx[2], collectionPlane);
332 double shwTwoTime = detprop->ConvertXToTicks(shwvtx[0]+shwdir[0], collectionPlane);
333 double shwTwoWire = geom->WireCoordinate(shwvtx[1]+shwdir[1], shwvtx[2]+shwdir[2], collectionPlane);
338 TH1F* ttemp_1 =
new TH1F(
"ttemp_1",
"ttemp_1",
TBINS,
TMIN,
TMAX);
339 TH1F* ttemp_2 =
new TH1F(
"ttemp_2",
"ttemp_2",
TBINS,
TMIN,
TMAX);
340 TH1F* ttemp_3 =
new TH1F(
"ttemp_3",
"ttemp_3",
TBINS,
TMIN,
TMAX);
341 TH1F* ttemp_4 =
new TH1F(
"ttemp_4",
"ttemp_4",
TBINS,
TMIN,
TMAX);
342 TH1F* ttemp_5 =
new TH1F(
"ttemp_5",
"ttemp_5",
TBINS,
TMIN,
TMAX);
344 for (
size_t i = 0; i < showerhits.size(); ++i) {
345 if (showerhits[i]->WireID().Plane != collectionPlane.Plane)
continue;
347 double wirePitch = geom->WirePitch(showerhits[i]->WireID());
348 double tickToDist = detprop->DriftVelocity(detprop->Efield(),detprop->Temperature());
349 tickToDist *= 1.e-3 * detprop->SamplingRate();
351 double xvtx = shwVtxTime * tickToDist;
352 double yvtx = shwVtxWire * wirePitch;
354 double xtwo = shwTwoTime * tickToDist;
355 double ytwo = shwTwoWire * wirePitch;
357 double xtwoorth = (ytwo - yvtx) + xvtx;
358 double ytwoorth = -(xtwo - xvtx) + yvtx;
360 double xhit = showerhits[i]->PeakTime() * tickToDist;
361 double yhit = showerhits[i]->WireID().Wire * wirePitch;
363 double ldist = std::abs((ytwoorth-yvtx)*xhit - (xtwoorth-xvtx)*yhit + xtwoorth*yvtx - ytwoorth*xvtx)/std::sqrt( pow((ytwoorth-yvtx), 2) + pow((xtwoorth-xvtx), 2) );
364 double tdist = ((ytwo-yvtx)*xhit - (xtwo-xvtx)*yhit + xtwo*yvtx - ytwo*xvtx)/std::sqrt( pow((ytwo-yvtx), 2) + pow((xtwo-xvtx), 2) );
366 double to3D = 1. / sqrt( pow(xvtx-xtwo,2) + pow(yvtx-ytwo,2) ) ;
374 ttemp->Fill(tdist, Q);
376 if (t < 1) ttemp_1->Fill(tdist, Q);
377 else if (t < 2) ttemp_2->Fill(tdist, Q);
378 else if (t < 3) ttemp_3->Fill(tdist, Q);
379 else if (t < 4) ttemp_4->Fill(tdist, Q);
380 else if (t < 5) ttemp_5->Fill(tdist, Q);
384 for (
int i = 0; i <
LBINS; ++i) {
385 if (ltemp->GetBinContent(i+1) == 0)
continue;
388 fLongitudinal->Fill(ltemp->GetBinCenter(i+1), elep, ltemp->GetBinContent(i+1));
391 for (
int i = 0; i <
TBINS; ++i) {
392 if (ttemp->GetBinContent(i+1) == 0)
continue;
395 fTransverse->Fill(ttemp->GetBinCenter(i+1), elep, ttemp->GetBinContent(i+1));
403 fTransverse_1->Fill(ttemp_1->GetBinCenter(i+1), elep, ttemp_1->GetBinContent(i+1));
404 fTransverse_2->Fill(ttemp_2->GetBinCenter(i+1), elep, ttemp_2->GetBinContent(i+1));
405 fTransverse_3->Fill(ttemp_3->GetBinCenter(i+1), elep, ttemp_3->GetBinContent(i+1));
406 fTransverse_4->Fill(ttemp_4->GetBinCenter(i+1), elep, ttemp_4->GetBinContent(i+1));
407 fTransverse_5->Fill(ttemp_5->GetBinCenter(i+1), elep, ttemp_5->GetBinContent(i+1));
419 auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
424 std::map<int,double> trkID_E;
429 TH1F* ttemp_1 =
new TH1F(
"ttemp_1",
"ttemp_1",
TBINS,
TMIN,
TMAX);
430 TH1F* ttemp_2 =
new TH1F(
"ttemp_2",
"ttemp_2",
TBINS,
TMIN,
TMAX);
431 TH1F* ttemp_3 =
new TH1F(
"ttemp_3",
"ttemp_3",
TBINS,
TMIN,
TMAX);
432 TH1F* ttemp_4 =
new TH1F(
"ttemp_4",
"ttemp_4",
TBINS,
TMIN,
TMAX);
433 TH1F* ttemp_5 =
new TH1F(
"ttemp_5",
"ttemp_5",
TBINS,
TMIN,
TMAX);
441 double shwvtxT = -999;
442 double shwvtxW = -999;
443 double shwtwoT = -999;
444 double shwtwoW = -999;
446 double shwvtxx = -999;
447 double shwvtxy = -999;
448 double shwtwox = -999;
449 double shwtwoy = -999;
450 double xtwoorth = -999;
451 double ytwoorth = -999;
453 double wirePitch = -999;
454 double tickToDist = -999;
456 bool foundParent =
false;
458 for (
size_t i = 0; i < allhits.size(); ++i) {
459 if (allhits[i]->WireID().Plane != collectionPlane.Plane)
continue;
464 for (
size_t j = 0; j < trackIDs.size(); ++j) {
466 if ( std::abs((piserv->
TrackIdToParticle_P(trackIDs[j].trackID))->PdgCode()) != 11)
continue;
467 if ( std::abs((piserv->
TrackIdToParticle_P(trackIDs[j].trackID))->Mother()) != 0)
continue;
478 shwvtxT = detprop->ConvertXToTicks(xvtx, collectionPlane);
479 shwvtxW = geom->WireCoordinate(yvtx, zvtx, collectionPlane);
481 shwtwoT = detprop->ConvertXToTicks(xtwo, collectionPlane);
482 shwtwoW = geom->WireCoordinate(ytwo, ztwo, collectionPlane);
484 wirePitch = geom->WirePitch(hit->
WireID());
485 tickToDist = detprop->DriftVelocity(detprop->Efield(),detprop->Temperature());
486 tickToDist *= 1.e-3 * detprop->SamplingRate();
488 shwvtxx = shwvtxT * tickToDist;
489 shwvtxy = shwvtxW * wirePitch;
491 shwtwox = shwtwoT * tickToDist;
492 shwtwoy = shwtwoW * wirePitch;
494 xtwoorth = (shwtwoy - shwvtxy) + shwvtxx;
495 ytwoorth = -(shwtwox - shwvtxx) + shwvtxy;
499 double xhit = hit->
PeakTime() * tickToDist;
502 double ldist = abs((ytwoorth-shwvtxy)*xhit - (xtwoorth-shwvtxx)*yhit + xtwoorth*shwvtxy - ytwoorth*shwvtxx)/sqrt( pow((ytwoorth-shwvtxy), 2) + pow((xtwoorth-shwvtxx), 2) );
503 double tdist = ((shwtwoy-shwvtxy)*xhit - (shwtwox-shwvtxx)*yhit + shwtwox*shwvtxy - shwtwoy*shwvtxx)/sqrt( pow((shwtwoy-shwvtxy), 2) + pow((shwtwox-shwvtxx), 2) );
505 double to3D = sqrt( pow(xvtx-xtwo , 2) + pow(yvtx-ytwo , 2) + pow(zvtx-ztwo , 2) ) / sqrt( pow(shwvtxx-shwtwox,2) + pow(shwvtxy-shwtwoy,2) ) ;
510 double energy = trackIDs[j].energy;
512 ltemp->Fill(t, energy);
513 ttemp->Fill(tdist, energy);
515 if (t < 1) ttemp_1->Fill(tdist, energy);
516 else if (t < 2) ttemp_2->Fill(tdist, energy);
517 else if (t < 3) ttemp_3->Fill(tdist, energy);
518 else if (t < 4) ttemp_4->Fill(tdist, energy);
519 else if (t < 5) ttemp_5->Fill(tdist, energy);
525 for (
int i = 0; i <
LBINS; ++i) {
526 if (ltemp->GetBinContent(i+1) == 0)
continue;
531 for (
int i = 0; i <
TBINS; ++i) {
532 if (ttemp->GetBinContent(i+1) == 0)
continue;
553 std::vector<sim::MCEnDep> alledep;
558 TH1F* ttemp_1 =
new TH1F(
"ttemp_1",
"ttemp_1",
TBINS,
TMIN,
TMAX);
559 TH1F* ttemp_2 =
new TH1F(
"ttemp_2",
"ttemp_2",
TBINS,
TMIN,
TMAX);
560 TH1F* ttemp_3 =
new TH1F(
"ttemp_3",
"ttemp_3",
TBINS,
TMIN,
TMAX);
561 TH1F* ttemp_4 =
new TH1F(
"ttemp_4",
"ttemp_4",
TBINS,
TMIN,
TMAX);
562 TH1F* ttemp_5 =
new TH1F(
"ttemp_5",
"ttemp_5",
TBINS,
TMIN,
TMAX);
565 for (
size_t i = 0; i < allchan.size(); ++i) {
570 for(
auto const& tdc_ide_pair : tdc_ide_map) {
571 for (
auto const& ide : tdc_ide_pair.second) {
580 alledep.push_back(edep);
587 double x0 = electron.
Vx();
588 double y0 = electron.
Vy();
589 double z0 = electron.
Vz();
591 double x2 = electron.
Px();
592 double y2 = electron.
Py();
593 double z2 = electron.
Pz();
595 TVector3 v0(x2, y2, z2);
598 for (
size_t i = 0; i < alledep.size(); ++i) {
599 double x = (double)alledep[i].Vertex()[0];
600 double y = (double)alledep[i].Vertex()[1];
601 double z = (double)alledep[i].Vertex()[2];
603 TVector3 v1(x-x0, y-y0, z-z0);
605 double ldist = v0.Dot(v1);
607 double tdist = (v0.Orthogonal()).Dot(v1);
609 double energy = alledep[i].Energy();
611 ltemp->Fill(t, energy);
612 ttemp->Fill(tdist, energy);
614 if (t < 1) ttemp_1->Fill(tdist, energy);
615 else if (t < 2) ttemp_2->Fill(tdist, energy);
616 else if (t < 3) ttemp_3->Fill(tdist, energy);
617 else if (t < 4) ttemp_4->Fill(tdist, energy);
618 else if (t < 5) ttemp_5->Fill(tdist, energy);
622 for (
int i = 0; i <
LBINS; ++i) {
623 if (ltemp->GetBinContent(i+1) == 0)
continue;
628 for (
int i = 0; i <
TBINS; ++i) {
629 if (ttemp->GetBinContent(i+1) == 0)
continue;
double E(const int i=0) const
const simb::MCParticle * TrackIdToParticle_P(int const &id)
TProfile * fTransverse4_photon
TCShowerTemplateMaker(fhicl::ParameterSet const &pset)
const simb::MCNeutrino & GetNeutrino() const
double Py(const int i=0) const
TProfile2D * fShowerProfileSimTrans2D_5
TProfile * fTransverse3_photon
TProfile * fLongitudinal_photon
geo::WireID WireID() const
Initial tdc tick for hit.
TProfile2D * fShowerProfileSimTrans2D_2
Declaration of signal hit object.
const simb::MCParticle & Nu() const
TProfile * fShowerProfileRecoTrans
std::string fHitModuleLabel
The data type to uniquely identify a Plane.
double Px(const int i=0) const
std::string fShowerModuleLabel
TProfile * fLongitudinal_electron
TProfile2D * fShowerProfileHitTrans2D_2
void SetTrackId(unsigned int id)
WireID_t Wire
Index of the wire within its plane.
TProfile * fTransverse4_electron
TProfile2D * fShowerProfileSimTrans2D_1
TProfile * fTransverse2_photon
TProfile2D * fShowerProfileHitTrans2D_4
TProfile * fShowerProfileSimLong
object containing MC flux information
TProfile2D * fShowerProfileHitTrans2D_1
Float_t y2[n_points_geant4]
TProfile * fShowerProfileSimTrans
TProfile * fTransverse1_photon
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
TProfile2D * fShowerProfileRecoTrans2D_4
const simb::MCParticle & Lepton() const
#define DEFINE_ART_MODULE(klass)
Provides recob::Track data product.
void showerProfile(std::vector< art::Ptr< recob::Hit > > showerhits, TVector3 shwvtx, TVector3 shwdir, double elep)
TProfile * fTransverse1_other
std::string fDigitModuleLabel
TProfile2D * fShowerProfileRecoTrans2D_5
void SetVertex(float x, float y, float z)
TProfile * fTransverse5_other
EDAnalyzer(Table< Config > const &config)
TProfile * fLongitudinal_other
Declaration of cluster object.
TProfile2D * fShowerProfileHitTrans2D_3
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
TProfile * fShowerProfileRecoLong
TProfile2D * fShowerProfileRecoTrans2D_2
void analyze(const art::Event &evt)
Detector simulation of raw signals on wires.
raw::ChannelID_t Channel() const
Returns the readout channel this object describes.
double Vx(const int i=0) const
float PeakTime() const
Time of the signal peak, in tick units.
TProfile2D * fShowerProfileSimLong2D
virtual ~TCShowerTemplateMaker()
T * make(ARGS...args) const
TProfile * fTransverse5_electron
TProfile2D * fShowerProfileSimTrans2D
Utility object to perform functions of association.
TProfile * fShowerProfileHitTrans
TProfile * fTransverse3_electron
const std::vector< sim::TrackIDE > HitToEveTrackIDEs(recob::Hit const &hit)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
TProfile * fShowerProfileHitLong
TProfile2D * fShowerProfileRecoTrans2D_3
TProfile2D * fShowerProfileHitTrans2D
double Pz(const int i=0) const
object containing MC truth information necessary for making RawDigits and doing back tracking ...
double Vz(const int i=0) const
TProfile2D * fShowerProfileSimTrans2D_4
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
TProfile * fTransverse4_other
calo::CalorimetryAlg fCalorimetryAlg
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
const simb::MCParticle * TrackIdToMotherParticle_P(int const &id)
TProfile2D * fShowerProfileHitLong2D
Float_t x2[n_points_geant4]
TProfile * fTransverse5_photon
TProfile2D * fShowerProfileRecoLong2D
double LifetimeCorrection(double time, double T0=0) const
TProfile2D * fShowerProfileRecoTrans2D_1
TProfile * fTransverse3_other
TProfile2D * fShowerProfileRecoTrans2D
TProfile * fTransverse2_electron
TProfile2D * fShowerProfileHitTrans2D_5
void showerProfileTrue(std::vector< art::Ptr< recob::Hit > > allhits, double elep)
double Vy(const int i=0) const
art framework interface to geometry description
void reconfigure(fhicl::ParameterSet const &pset)
TProfile2D * fShowerProfileSimTrans2D_3
TProfile * fTransverse1_electron
std::string fGenieGenModuleLabel
TProfile * fTransverse2_other