20 #include "Utilities/AssociationUtil.h" 237 DedxSample = tfs->
make<TTree>(
"fDedxSample",
"Variables for Creating Dedx Sample Histograms");
297 fPOT = tfs->
make<TH1D>(
"POT",
";;Total POT",1, 0.0, 1.0);
302 fTrueProtonEFromTrackLength = tfs->
make<TH2D>(
"TrueProtonEFromTrackLength",
";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
348 mf::LogError(
"Susie") <<
"Spill Data not found for real data event";
355 if ((spillPot->
spillpot)*1e12 < 2e12)
return;
356 if ((spillPot->
hornI < -202)||(spillPot->
hornI > -198))
return;
357 if ((spillPot->
posx < -2.00)||(spillPot->
posx > 2.00))
return;
358 if ((spillPot->
posy < -2.00)||(spillPot->
posy > 2.00))
return;
359 if ((spillPot->
widthx < 0.57)||(spillPot->
widthx > 1.58))
return;
360 if ((spillPot->
widthy < 0.57)||(spillPot->
widthy > 1.58))
return;
370 mf::LogError(
"CAFMaker") <<
"Spill EventQuality not found for real data event";
379 float spillpot = spillPot->
spillpot;
390 for(
size_t h = 0;
h < hithdl->size(); ++
h){
404 for(
unsigned int i = 0;
i<slicecol->size();++
i){
422 for(
unsigned int iSlice = 0; iSlice<slicelist.
size(); ++iSlice){
425 if(slicelist[iSlice]->IsNoise()){
continue; }
428 std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
430 if(tracks.size() == 0){
continue; }
433 if(tracks.size() != 1){
continue; }
448 for(
unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
449 if(!tracks[itrk]->Is3D()){
455 if(pid->
Value() > bestpid){
456 bestpid = pid->
Value();
461 if (!all3D){
continue; }
476 if(bestidx == -1 || bestpid < 0.8){
continue; }
478 if((bestidx != 0) && (bestidx !=1)){
479 std::cout<<
"BUMMER! A weird track id for muon that shouldn't happen!!"<<bestidx<<
std::endl;
497 fzStart = tracks[bestidx]->Start().Z();
498 fzEnd = tracks[bestidx]->Stop().Z();
502 fxStart = tracks[bestidx]->Start().X();
503 fyStart = tracks[bestidx]->Start().Y();
504 fxEnd = tracks[bestidx]->Stop().X();
505 fyEnd = tracks[bestidx]->Stop().Y();
507 fxDir = tracks[bestidx]->Dir().X();
508 fyDir = tracks[bestidx]->Dir().Y();
509 fzDir = tracks[bestidx]->Dir().Z();
511 unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();
513 for(
int itrj = ntraj-1; itrj >= 0; --itrj){
514 TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
518 unsigned int minPlane = tracks[bestidx]->MinPlane();
519 unsigned int maxPlane = tracks[bestidx]->MaxPlane();
522 double hitCosx[tracks[bestidx]->NXCell()];
523 double hitCosy[tracks[bestidx]->NYCell()];
525 for(
unsigned int ihit = 0; ihit<tracks[bestidx]->NXCell(); ++ihit){
526 if(tracks[bestidx]->XCell(ihit)->Plane() > maxXPlane){ maxXPlane = tracks[bestidx]->XCell(ihit)->Plane(); }
530 double dz = geom->
Plane(tracks[bestidx]->XCell(ihit)->Plane())->
Cell(0)->
HalfD();
532 double minz = xyzplane[2]-
dz;
533 double maxz = xyzplane[2]+
dz;
534 unsigned int currtrj = 0;
536 for(
unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
537 if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
539 TVector3 direction = tracks[bestidx]->Dir();
540 if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
541 direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
543 direction = direction.Unit();
544 hitCosx[ihit] =
sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
549 for(
unsigned int ihit = 0; ihit<tracks[bestidx]->NYCell(); ++ihit){
550 if(tracks[bestidx]->YCell(ihit)->Plane() > maxYPlane){ maxYPlane = tracks[bestidx]->YCell(ihit)->Plane(); }
554 double dz = geom->
Plane(tracks[bestidx]->YCell(ihit)->Plane())->
Cell(0)->
HalfD();
556 double minz = xyzplane[2]-
dz;
557 double maxz = xyzplane[2]+
dz;
558 unsigned int currtrj = 0;
560 for(
unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
561 if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
563 TVector3 direction = tracks[bestidx]->Dir();
564 if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
565 direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
567 direction = direction.Unit();
568 hitCosy[ihit] =
sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
574 std::vector<double> dedxs;
575 std::vector<double> des;
576 std::vector<double> pecorrs;
577 std::vector<double> mips;
578 std::vector<double> dxs;
579 std::vector<double> deads;
580 std::vector<double> xdedxs;
581 std::vector<unsigned int> dedxpls;
583 std::vector<double> pes;
584 std::vector<double> adcs;
585 std::vector<double> xhits;
586 std::vector<double> yhits;
587 std::vector<double> zhits;
588 std::vector<double> thits;
589 std::vector<int> views;
590 std::vector<int> cells;
591 std::vector<int> planes;
593 std::vector<double> mcdes;
594 std::vector<double> mcdebs;
595 std::vector<double> mcdxs;
597 std::vector< std::vector<int> > mcells;
598 std::vector< std::vector<float> > madcs;
599 std::vector< std::vector<float> > mpes;
600 std::vector< std::vector<float> > mpecorrs;
602 for(
unsigned int iPlane = minPlane;iPlane<=
maxPlane;++iPlane){
604 double planePECorr = 0;
617 std::vector<int> mCell;
618 std::vector<float> mPE;
619 std::vector<float> mADC;
620 std::vector<float> mPECorr;
623 double mcPlaneEBirks = 0;
624 double mcPathLength = 0;
625 double totcosinview = 0;
626 double avgCosOtherView = 1.0;
627 bool useOtherView =
true;
637 double nhitsother = 0;
650 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
651 if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane && tracks[bestidx]->XCell(iHit)->View() !=
geo::kY){
653 planeADC+=chit->
ADC();
654 planePE +=chit->
PE();
656 plane = chit->
Plane();
659 if(minPlane == maxPlane){
660 rhit = cal->
MakeRecoHit(*(tracks[bestidx]->XCell(iHit)), 0.25);
663 planeGeV+=rhit.
GeV();
664 planePECorr+=rhit.
PECorr();
665 planeMIP+=rhit.
MIP();
670 totcosinview+=hitCosx[iHit];
671 mCell.push_back( chit->
Plane() );
672 mPE.push_back( chit->
PE() );
673 mADC.push_back( chit->
ADC() );
674 mPECorr.push_back(rhit.
PECorr() );
677 if(tracks[bestidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->XCell(iHit)->Cell();}
678 if(tracks[bestidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->XCell(iHit)->Cell();}
684 const std::vector< sim::FLSHit > flsHits =
bt->
HitToFLSHit(chit);
685 for(
int j = 0;
j < (
int)flsHits.size();
j++){
686 mcPlaneE += flsHits[
j].GetEdep();
687 mcPlaneEBirks += flsHits[
j].GetEdepBirks();
688 mcPathLength += flsHits[
j].GetTotalPathLength();
698 double nhitsbefore = 0;
699 double nhitsafter = 0;
700 double totCosBefore = 0;
701 double totCosAfter = 0;
702 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
703 if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane-1){
705 totCosBefore+=hitCosy[iHit];
707 else if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane+1){
709 totCosAfter+=hitCosy[iHit];
712 if(nhitsbefore > 0 && nhitsafter > 0){
713 avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
714 nhitsother = (nhitsbefore+nhitsafter)/2;
716 else if(nhitsbefore > 0){
717 avgCosOtherView = totCosBefore/nhitsbefore;
718 nhitsother = nhitsbefore;
720 else if(nhitsafter > 0){
721 avgCosOtherView = totCosAfter/nhitsafter;
722 nhitsother = nhitsafter;
724 else{ useOtherView =
false; }
729 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
730 if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane && tracks[bestidx]->YCell(iHit)->View() !=
geo::kX 731 && tracks[bestidx]->YCell(iHit)->View() !=
geo::kXorY){
733 planeADC+=chit->
ADC();
734 planePE +=chit->
PE();
736 plane = chit->
Plane();
739 if(minPlane == maxPlane){
740 rhit = cal->
MakeRecoHit(*(tracks[bestidx]->YCell(iHit)), 0.25);
743 planeGeV+=rhit.
GeV();
744 planePECorr+=rhit.
PECorr();
745 planeMIP+=rhit.
MIP();
750 totcosinview+=hitCosy[iHit];
752 mCell.push_back( chit->
Plane() );
753 mPE.push_back( chit->
PE() );
754 mADC.push_back( chit->
ADC() );
755 mPECorr.push_back(rhit.
PECorr() );
758 if(tracks[bestidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->YCell(iHit)->Cell();}
759 if(tracks[bestidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->YCell(iHit)->Cell();}
765 const std::vector< sim::FLSHit > flsHits =
bt->
HitToFLSHit(chit);
766 for(
int j = 0;
j < (
int)flsHits.size();
j++){
767 mcPlaneE += flsHits[
j].GetEdep();
768 mcPlaneEBirks += flsHits[
j].GetEdepBirks();
769 mcPathLength += flsHits[
j].GetTotalPathLength();
779 double nhitsbefore = 0;
780 double nhitsafter = 0;
781 double totCosBefore = 0;
782 double totCosAfter = 0;
783 for(
unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
784 if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane-1){
786 totCosBefore+=hitCosx[iHit];
788 else if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane+1){
790 totCosAfter+=hitCosx[iHit];
793 if(nhitsbefore > 0 && nhitsafter > 0){
794 avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
795 nhitsother = (nhitsbefore+nhitsafter)/2;
797 else if(nhitsbefore > 0){
798 avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
799 nhitsother = nhitsbefore;
801 else if(nhitsafter > 0){
802 avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
803 nhitsother = nhitsbefore;
805 else{ useOtherView =
false; }
815 double avgCos = TMath::Abs(totcosinview/nhits);
824 double liveLength = 0;
825 for(
int i = minCell;
i<maxCell+1;++
i){
829 if(
i == minCell ||
i == maxCell){ liveLength+=c->
HalfW(); }
830 else{ liveLength+=(2.0*c->
HalfW() ); }
833 double deadLength = 0;
834 if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
835 double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
836 double planePathLength =
TMath::Sqrt(planewidth*planewidth + planeheight*planeheight);
839 if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }
841 double deltaz = planewidth;
842 double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
843 if(avgCos == 0 || avgCos!=avgCos){
847 if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
851 double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
855 if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
856 deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
858 planePathLength =
TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
861 if(minPlane == maxPlane){ planePathLength = liveLength; }
863 assert(planeGeV/(planePathLength)>=0);
864 dedxs.push_back(planeGeV/(planePathLength));
865 des.push_back(planeGeV);
866 pecorrs.push_back(planePECorr);
867 mips.push_back(planeMIP);
868 dxs.push_back(planePathLength);
869 deads.push_back(deadLength);
870 xhits.push_back(planeX/nhits);
871 yhits.push_back(planeY/nhits);
872 zhits.push_back(planeZ/nhits);
873 thits.push_back(planeT/nhits);
874 views.push_back(planeView);
875 planes.push_back(plane);
876 cells.push_back(cell);
877 pes.push_back(planePE);
878 adcs.push_back(planeADC);
880 mcdes.push_back(mcPlaneE);
881 mcdebs.push_back(mcPlaneEBirks);
882 mcdxs.push_back(mcPathLength);
884 mcells.push_back(mCell);
885 madcs.push_back(mADC);
887 mpecorrs.push_back(mPECorr);
889 mcdes.push_back(mcPlaneE);
890 mcdebs.push_back(mcPlaneEBirks);
891 mcdxs.push_back(mcPathLength);
902 if (!fndmnyNumuEnergy.isValid()){
continue; }
903 std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
904 if (energyVec.empty()){
continue; }
905 numuE = *energyVec[0];
908 if (numuE.
E() == -5){
continue; }
915 if (!fndmnyQepid.isValid()){
continue; }
916 std::vector<art::Ptr<qeef::QePId> > qVec = fndmnyQepid.at(iSlice);
917 if (qVec.empty()){
continue; }
927 double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy());
931 int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell());
944 if (!funcReturn.empty()){
947 const std::vector<const sim::Particle*> test_mu =
bt->
HitsToParticle(tracks[bestidx]->AllCells());
948 if (!test_mu.empty()) {
949 fpdg = test_mu[0]->PdgCode();
968 for(
unsigned int idedx = 0; idedx<pid->
NDedx(); ++idedx){
980 fxHit = xhits[idedx];
981 fyHit = yhits[idedx];
982 fzHit = zhits[idedx];
983 ftHit = thits[idedx];
984 fView = views[idedx];
986 fCell = cells[idedx];
995 fmcDe = mcdes[idedx];
997 fmcDx = mcdxs[idedx];
1021 bool contained =
false;
1027 unsigned int slcncellsfromedge = 999999999;
1028 for(
unsigned int ihit = 0; ihit < slice->
NCell(); ++ihit){
1032 if(cellsfromedge < slcncellsfromedge){
1033 slcncellsfromedge = cellsfromedge;
1039 unsigned int firstpl = slice->
MinPlane();
1040 unsigned int lastpl = slice->
MaxPlane();
1042 TVector3 dirbwk = -trk->
Dir();
1046 int sliceNHit = slice->
NCell();
1049 double slcTime = slice->
MeanTNS();
1052 int trackNHit = trk->
NCell();
1056 if((sliceNHit > 20)&&(slcncellsfromedge > 1)&&(slcTime>50000)&&(slcTime<450000)&&(trackNHit > 5)){
1057 if(firstpl > 1 && lastpl < 212){
1062 if(fwdcell > 4 && bwkcell > 8 && trk->
Start().Z() < 1150){
1065 if(trk->
Stop().Z() < 1275){
1091 if( ypostranKirk < 55){
size_t NTrajectoryPoints() const
back track the reconstruction to the simulation
double Dedx(unsigned int i) const
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
TH1D * fSlicesPassingSelectionCuts
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.
bool ContainedEvent(art::Ptr< rb::Cluster > slice, art::Ptr< rb::Track > trk)
bool fTrueNumuCCWithProton
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
unsigned short Plane() const
const CellGeo * Cell(int icell) const
TH1D * fSuperTrueSlicesPassingSelectionCuts
Vertical planes which measure X.
unsigned int Ncells() const
Number of cells in this plane.
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
art::ServiceHandle< photrans::FiberBrightness > fb
void beginRun(const art::Run &run)
art::ServiceHandle< geo::Geometry > fgeom
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
art::ServiceHandle< cheat::BackTracker > bt
virtual TVector3 Start() const
int nmissingdcms
of missing DCMs
Horizontal planes which measure Y.
std::string fNumuEnergyLabel
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
bool DedxVertex(unsigned int i) const
unsigned int getBrightnessIndex(unsigned int plane, unsigned int cell) const
unsigned short Cell() const
Track finder for cosmic rays.
View_t View() const
Which coordinate does this plane measure.
bool fSuperTrueNumuCCWithProton
signed long long int deltaspilltimensec
double fracdcm3hits
fraction of DCM3 hits in horizontal modules
void push_back(Ptr< U > const &p)
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
int NMeasurementPlanes() const
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.
std::string fHitModuleLabel
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
TTree * DedxSample
Tree for holding dE/dx variables needed for making signal sample histograms.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
bool DedxUsed(unsigned int i) const
std::string fPidLabel
Where to find the pids.
std::vector< int > fMultiCell
EDAnalyzer(Table< Config > const &config)
TH1D * fEventsPassingSpillCuts
ReMIdDedxRock(fhicl::ParameterSet const &pset)
std::string fSpillQualLabel
float NDHadCalTran() const
unsigned int NDedx() const
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
std::vector< float > fMultiADC
std::string fGeneratorLabel
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
double DedxLocation(unsigned int i) const
art::ServiceHandle< geo::LiveGeometry > flivegeom
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.
bool fPassHadEOnMuonTrack
std::string fSliceLabel
Where to find reconstructed slices.
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
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?)
std::vector< float > fMultiPECorr
void analyze(art::Event const &evt)
assert(nhit_max >=nhit_nbins)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
art::ServiceHandle< nova::dbi::RunHistoryService > frh
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
std::string fNumiBeamLabel
Encapsulate the cell geometry.
TH1D * fPOT
Histogram of POT passing spill cuts.
TH1D * fTrueSlicesPassingSelectionCuts
TH2D * fSuperTrueProtonEFromTrackLength
TH2D * fTrueProtonEFromTrackLength
float NDHadCalCat() const
double spillpot
POT for spill normalized by 10^12.
std::vector< float > fMultiPE
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.