21 #include "Utilities/AssociationUtil.h" 210 DedxSample = tfs->
make<TTree>(
"fDedxSample",
"Variables for Creating Dedx Sample Histograms");
255 fPOT = tfs->
make<TH1D>(
"POT",
";;Total POT",1, 0.0, 1.0);
295 for(
size_t h = 0;
h < hithdl->size(); ++
h){
306 for(
unsigned int i = 0;
i<slicecol->size();++
i){
318 for(
unsigned int iSlice = 0; iSlice<slicelist.
size(); ++iSlice){
321 if(slicelist[iSlice]->IsNoise()){
continue; }
324 std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
326 if(tracks.size() == 0){
continue; }
329 if(tracks.size() != 1){
continue; }
336 if(!tracks[bestidx]->Is3D()) is3D =
false;
337 if (!is3D){
continue; }
361 bool isStopper =
true;
362 if (distEndToEdge < 50) isStopper =
false;
365 if (!isStopper)
continue;
368 fzStart = tracks[bestidx]->Start().Z();
369 fzEnd = tracks[bestidx]->Stop().Z();
372 fxStart = tracks[bestidx]->Start().X();
373 fyStart = tracks[bestidx]->Start().Y();
374 fxEnd = tracks[bestidx]->Stop().X();
375 fyEnd = tracks[bestidx]->Stop().Y();
377 fxDir = tracks[bestidx]->Dir().X();
378 fyDir = tracks[bestidx]->Dir().Y();
379 fzDir = tracks[bestidx]->Dir().Z();
381 double trklen = tracks[bestidx]->TotalLength();
383 unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();
386 std::vector<double> lengths(ntraj);
389 for(
int itrj = ntraj-1; itrj >= 0; --itrj){
390 TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
391 if(itrj == (
int)ntraj-1) lengths[itrj] = 0.0;
393 TVector3 lasttrajpoint = tracks[bestidx]->TrajectoryPoint(itrj+1);
394 lengths[itrj] = (trajpoint-lasttrajpoint).
Mag() + lengths[itrj+1];
398 unsigned int minPlane = tracks[bestidx]->MinPlane();
399 unsigned int maxPlane = tracks[bestidx]->MaxPlane();
402 double hitCosx[tracks[bestidx]->NXCell()];
403 double hitCosy[tracks[bestidx]->NYCell()];
405 for(
unsigned int ihit = 0; ihit<tracks[bestidx]->NXCell(); ++ihit){
406 if(tracks[bestidx]->XCell(ihit)->Plane() > maxXPlane){ maxXPlane = tracks[bestidx]->XCell(ihit)->Plane(); }
410 double dz = geom->
Plane(tracks[bestidx]->XCell(ihit)->Plane())->
Cell(0)->
HalfD();
412 double minz = xyzplane[2]-
dz;
413 double maxz = xyzplane[2]+
dz;
414 unsigned int currtrj = 0;
416 for(
unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
417 if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
419 TVector3 direction = tracks[bestidx]->Dir();
420 if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
421 direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
423 direction = direction.Unit();
424 hitCosx[ihit] =
sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
427 for(
unsigned int ihit = 0; ihit<tracks[bestidx]->NYCell(); ++ihit){
428 if(tracks[bestidx]->YCell(ihit)->Plane() > maxYPlane){ maxYPlane = tracks[bestidx]->YCell(ihit)->Plane(); }
432 double dz = geom->
Plane(tracks[bestidx]->YCell(ihit)->Plane())->
Cell(0)->
HalfD();
434 double minz = xyzplane[2]-
dz;
435 double maxz = xyzplane[2]+
dz;
436 unsigned int currtrj = 0;
438 for(
unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
439 if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
441 TVector3 direction = tracks[bestidx]->Dir();
442 if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
443 direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
445 direction = direction.Unit();
446 hitCosy[ihit] =
sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
452 std::vector<double> dedxs;
453 std::vector<double> des;
454 std::vector<double> pecorrs;
455 std::vector<double> mips;
456 std::vector<double> dxs;
457 std::vector<double> deads;
458 std::vector<double> xdedxs;
459 std::vector<unsigned int> dedxpls;
461 std::vector<double> pes;
462 std::vector<double> adcs;
463 std::vector<double> xhits;
464 std::vector<double> yhits;
465 std::vector<double> zhits;
466 std::vector<double> thits;
467 std::vector<int> views;
468 std::vector<int> cells;
469 std::vector<int> planes;
471 std::vector<double> mcdes;
472 std::vector<double> mcdebs;
473 std::vector<double> mcdxs;
475 for(
unsigned int iPlane = minPlane;iPlane<=
maxPlane;++iPlane){
478 double planePECorr = 0;
492 double mcPlaneEBirks = 0;
493 double mcPathLength = 0;
494 double totcosinview = 0;
495 double avgCosOtherView = 1.0;
496 bool useOtherView =
true;
502 double minz = xyzplane[2]-
dz;
503 double maxz = xyzplane[2]+
dz;
507 double nhitsother = 0;
512 for(
unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
513 if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){matchpt = itraj;}
519 double xlength = trklen - lengths[matchpt];
522 if(minPlane == maxPlane){xlength = trklen/2;}
523 if(tracks[bestidx]->Stop().Y()>tracks[bestidx]->Start().Y()){xlength = (trklen-xlength);}
527 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
528 if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane && tracks[bestidx]->XCell(iHit)->View() !=
geo::kY){
530 planeADC+=chit->
ADC();
531 planePE +=chit->
PE();
533 plane = chit->
Plane();
535 if(minPlane == maxPlane){
536 rhit = cal->
MakeRecoHit(*(tracks[bestidx]->XCell(iHit)), 0.25);
539 planeGeV+=rhit.
GeV();
540 planePECorr+=rhit.
PECorr();
541 planeMIP+=rhit.
MIP();
546 totcosinview+=hitCosx[iHit];
548 if(tracks[bestidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->XCell(iHit)->Cell();}
549 if(tracks[bestidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->XCell(iHit)->Cell();}
555 const std::vector< sim::FLSHit > flsHits =
bt->
HitToFLSHit(chit);
556 for(
int j = 0;
j < (
int)flsHits.size();
j++){
557 mcPlaneE += flsHits[
j].GetEdep();
558 mcPlaneEBirks += flsHits[
j].GetEdepBirks();
559 mcPathLength += flsHits[
j].GetTotalPathLength();
569 double nhitsbefore = 0;
570 double nhitsafter = 0;
571 double totCosBefore = 0;
572 double totCosAfter = 0;
573 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
574 if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane-1){
576 totCosBefore+=hitCosy[iHit];
578 else if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane+1){
580 totCosAfter+=hitCosy[iHit];
583 if(nhitsbefore > 0 && nhitsafter > 0){
584 avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
585 nhitsother = (nhitsbefore+nhitsafter)/2;
587 else if(nhitsbefore > 0){
588 avgCosOtherView = totCosBefore/nhitsbefore;
589 nhitsother = nhitsbefore;
591 else if(nhitsafter > 0){
592 avgCosOtherView = totCosAfter/nhitsafter;
593 nhitsother = nhitsafter;
595 else{ useOtherView =
false; }
600 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
601 if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane && tracks[bestidx]->YCell(iHit)->View() !=
geo::kX 602 && tracks[bestidx]->YCell(iHit)->View() !=
geo::kXorY){
604 planeADC+=chit->
ADC();
605 planePE +=chit->
PE();
607 plane = chit->
Plane();
610 if(minPlane == maxPlane){
611 rhit = cal->
MakeRecoHit(*(tracks[bestidx]->YCell(iHit)), 0.25);
614 planeGeV+=rhit.
GeV();
615 planePECorr+=rhit.
PECorr();
616 planeMIP+=rhit.
MIP();
621 totcosinview+=hitCosy[iHit];
623 if(tracks[bestidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->YCell(iHit)->Cell();}
624 if(tracks[bestidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->YCell(iHit)->Cell();}
630 const std::vector< sim::FLSHit > flsHits =
bt->
HitToFLSHit(chit);
631 for(
int j = 0;
j < (
int)flsHits.size();
j++){
632 mcPlaneE += flsHits[
j].GetEdep();
633 mcPlaneEBirks += flsHits[
j].GetEdepBirks();
634 mcPathLength += flsHits[
j].GetTotalPathLength();
644 double nhitsbefore = 0;
645 double nhitsafter = 0;
646 double totCosBefore = 0;
647 double totCosAfter = 0;
648 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
649 if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane-1){
651 totCosBefore+=hitCosx[iHit];
653 else if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane+1){
655 totCosAfter+=hitCosx[iHit];
658 if(nhitsbefore > 0 && nhitsafter > 0){
659 avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
660 nhitsother = (nhitsbefore+nhitsafter)/2;
662 else if(nhitsbefore > 0){
663 avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
664 nhitsother = nhitsbefore;
666 else if(nhitsafter > 0){
667 avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
668 nhitsother = nhitsbefore;
670 else{ useOtherView =
false; }
680 double avgCos = TMath::Abs(totcosinview/nhits);
689 double liveLength = 0;
690 for(
int i = minCell;
i<maxCell+1;++
i){
694 if(
i == minCell ||
i == maxCell){ liveLength+=c->
HalfW(); }
695 else{ liveLength+=(2.0*c->
HalfW() ); }
698 double deadLength = 0;
699 if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
700 double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
701 double planePathLength =
TMath::Sqrt(planewidth*planewidth + planeheight*planeheight);
704 if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }
706 double deltaz = planewidth;
707 double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
708 if(avgCos == 0 || avgCos!=avgCos){
712 if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
716 double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
720 if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
721 deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
723 planePathLength =
TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
726 if(minPlane == maxPlane){ planePathLength = liveLength; }
728 xdedxs.push_back(xlength);
730 assert(planeGeV/(planePathLength)>=0);
731 dedxs.push_back(planeGeV/(planePathLength));
732 des.push_back(planeGeV);
733 pecorrs.push_back(planePECorr);
734 mips.push_back(planeMIP);
735 dxs.push_back(planePathLength);
736 deads.push_back(deadLength);
737 xhits.push_back(planeX/nhits);
738 yhits.push_back(planeY/nhits);
739 zhits.push_back(planeZ/nhits);
740 thits.push_back(planeT/nhits);
741 views.push_back(planeView);
742 planes.push_back(plane);
743 cells.push_back(cell);
744 pes.push_back(planePE);
745 adcs.push_back(planeADC);
747 mcdes.push_back(mcPlaneE);
748 mcdebs.push_back(mcPlaneEBirks);
749 mcdxs.push_back(mcPathLength);
758 if (!fndmnyNumuEnergy.isValid()){
continue; }
759 std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
760 if (energyVec.empty()){
continue; }
761 numuE = *energyVec[0];
764 if (numuE.
E() == -5){
continue; }
769 double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy());
773 int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell());
779 const std::vector<const sim::Particle*> test_mu =
bt->
HitsToParticle(tracks[bestidx]->AllCells());
780 if (!test_mu.empty()) {
781 fpdg = test_mu[0]->PdgCode();
792 for(
int idedx = 0; idedx<
fnDedx; ++idedx){
801 fxHit = xhits[idedx];
802 fyHit = yhits[idedx];
803 fzHit = zhits[idedx];
804 ftHit = thits[idedx];
805 fView = views[idedx];
807 fCell = cells[idedx];
811 fmcDe = mcdes[idedx];
813 fmcDx = mcdxs[idedx];
818 fDistDedx = (1.0 -xdedxs[idedx]/
trklen)*tracks[bestidx]->TotalLength();
842 bool contained =
false;
852 TVector3 dirbwk = -trk->
Dir();
856 int sliceNHit = slice->
NCell();
862 int trackNHit = trk->
NCell();
866 if((sliceNHit > 20)&&(trackNHit > 5)){
871 if(fwdcell > 4) contained =
true;
size_t NTrajectoryPoints() const
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void GetCenter(double *xyz, double localz=0.0) const
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
std::string fSpillQualLabel
art::ServiceHandle< photrans::FiberBrightness > fb
unsigned short Plane() const
double DistToEdgeXY(TVector3 vertex)
const CellGeo * Cell(int icell) const
Vertical planes which measure X.
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
std::string fNumiBeamLabel
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
Horizontal planes which measure Y.
art::ServiceHandle< nova::dbi::RunHistoryService > frh
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
unsigned short Cell() const
Track finder for cosmic rays.
std::string fGeneratorLabel
art::ServiceHandle< cheat::BackTracker > bt
View_t View() const
Which coordinate does this plane measure.
Far Detector at Ash River, MN.
double DistToEdgeZ(TVector3 vertex)
void push_back(Ptr< U > const &p)
TH1D * fPOT
Histogram of POT passing spill cuts.
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
virtual TVector3 Dir() const
Unit vector describing prong direction.
Geometry information for a single readout plane.
TH1D * fTrueSlicesPassingSelectionCuts
std::string fNumuEnergyLabel
Collect Geo headers and supply basic geometry functions.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
EDAnalyzer(Table< Config > const &config)
std::string fSliceLabel
Where to find reconstructed slices.
std::string fHitModuleLabel
TH1D * fSlicesPassingSelectionCuts
art::ServiceHandle< geo::LiveGeometry > flivegeom
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
std::string fTrackLabel
Where to find recondtructed tracks.
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
art::ServiceHandle< geo::Geometry > fgeom
assert(nhit_max >=nhit_nbins)
void beginRun(const art::Run &run)
ReMIdDedxFD(fhicl::ParameterSet const &pset)
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
Encapsulate the cell geometry.
bool fPassHadEOnMuonTrack
std::string fPidLabel
Where to find the pids.
void analyze(art::Event const &evt)
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.