10 #define CosmicTrends_h 39 #include "TStopwatch.h" 48 class StopperThreshold;
73 void FillTree(
rb::Track const& track, std::set<TVector3,
bool(*)(TVector3
const&, TVector3
const&) >
const& trajPoints);
74 void FillHist(
rb::Track const& track, std::set<TVector3,
bool(*)(TVector3
const&, TVector3
const&) >
const& trajPoints);
75 void testPath(
rb::Track const& track, std::set<TVector3,
bool(*)(TVector3
const&, TVector3
const&) >
const& trajPoints);
77 int CalculateWall(
float x0,
float z0,
float dx,
float dz,
float xlow,
float xhigh,
float zlow,
float zhigh);
166 return a.Z() < b.Z();
187 tInfoTree = tfs->
make<TTree>(
"InfoTree",
"Shadow and Threshold variables");
219 float maxADCrate = 200.0;
222 int wBins = detHalfWidth/5;
224 TString sview[4] = {
"Xall",
"Yall",
"Xmip",
"Ymip"};
226 for(
int i = 0;
i < 4; ++
i){
229 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
232 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
235 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
238 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
241 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
244 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADCrate);
247 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,20);
250 wBins,-1,1,nADCbin,0,20);
253 wBins,-1,1,nADCbin,0,20);
256 wBins,-1,1,nADCbin,0,20);
259 5,-0.5,4.5,nADCbin,0,20);
262 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,20);
265 wBins, -detHalfWidth, detHalfWidth, nADCbin, 0, 5);
268 wBins, -detHalfWidth, detHalfWidth, nADCbin, 0, 5);
271 wBins,-detHalfWidth,detHalfWidth,nADCbin,0,maxADC);
274 nADCbin,0,maxADC,nADCbin,0,maxADC);
277 wBins,0,20.0,nADCbin,0,maxADCrate);
280 wBins,0,20.0,nADCbin,0,maxADC);
283 wBins,0,20.0,nADCbin,0,maxADC);
286 nADCbin,0,maxADCrate);
304 wBins,-detHalfWidth,detHalfWidth);
307 wBins,-detHalfWidth,detHalfWidth);
362 for(
size_t t = 0;
t < th->size(); ++
t){
364 if(!track.
Is3D())
continue;
365 std::set<TVector3, bool(*)(TVector3 const&, TVector3 const&)> trajPoints(
ltTVec);
367 if(trajPoints.size() == 0)
continue;
381 std::set<TVector3,
bool(*)(TVector3
const&, TVector3
const&) > & trajPoints)
387 double totalDist = 0.;
407 std::set<TVector3,
bool(*)(TVector3
const&, TVector3
const&)>
const& trajPoints)
409 float minZ = trajPoints.begin()->Z();
410 float maxZ = trajPoints.rbegin()->Z();
419 minZ -= cellHalfDepth;
420 maxZ += cellHalfDepth;
421 float xi = trajPoints.begin()->X();
422 float yi = trajPoints.begin()->Y();
423 float zi = trajPoints.begin()->Z();
424 float xf = trajPoints.rbegin()->X();
425 float yf = trajPoints.rbegin()->Y();
426 float zf = trajPoints.rbegin()->Z();
427 float dr =
sqrt((xi-xf)*(xi-xf)+(yi-yf)*(yi-yf)+(zi-zf)*(zi-zf));
429 float dx = (xf-xi)/dr;
430 float dy = (yf-yi)/dr;
431 float dz = (zf-zi)/dr;
434 std::pair<uint32_t, uint32_t>
pc;
436 std::set<double> planeZBounds;
439 if(boundaryMap.size() < 1)
return;
440 for(
size_t c = 0;
c < track.
NCell(); ++
c){
443 float fpoissonlambda = 0.0;
444 float ftruemev = 0.0;
445 float ftruepath = -10.0;
446 float ftruew = -9999.0;
452 fpoissonlambda += pe.PoissonLambda();
457 for(
size_t n = 0;
n < ides.size(); ++
n)
458 ftruemev += ides[
n].
energy*1000;
464 auto particle = parts.front();
465 std::vector<sim::FLSHit> flsHits =
fBT->
ParticleToFLSHit(particle->TrackId(), particle->PdgCode());
466 for(
auto fls : flsHits){
467 if(fls.GetPlaneID() == track.
Cell(
c)->
Plane() &&
468 fls.GetCellID() == track.
Cell(
c)->
Cell() ){
469 ftruepath = fls.GetTotalPathLength();
470 ftruew = fls.GetZAverage();
476 if((ftruemev == 0.0) || (fpoissonlambda == 0.0))
continue;
505 for(
unsigned int tp=0;tp<ntp;tp++){
524 for(
unsigned int tp=0;tp<ntp;tp++){
534 vdir = (vtp1 - vtp0).
Unit();
537 pc.first = ch->
Plane();
538 pc.second = ch->
Cell();
540 point = TVector3(rh.
X(), rh.
Y(), rh.
Z());
559 double fadc = ch->
ADC();
560 double fpe = ch->
PE();
562 double fpecorr = -5.0;
570 if((ch->
View() ==
geo::kX) && (vdir.X() != 0 || vdir.Z() != 0)){
573 (xyz[0]-cellwidth), (xyz[0]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
576 fw = (xyz[2]-vtp0.Z())*vdir.Y()/vdir.Z() + vtp0.Y();
579 fw = (xyz[2]-zi)*dy/dz + yi;
583 if((ch->
View() ==
geo::kY) && (vdir.Y() != 0 || vdir.Z() != 0)){
586 (xyz[1]-cellwidth), (xyz[1]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
589 fw = (xyz[2]-vtp0.Z())*vdir.X()/vdir.Z() + vtp0.X();
592 fw = (xyz[2]-zi)*dx/dz + xi;
596 if((fpath < 1.0) || (fpath > 25.0) || (
abs(fw) > detHalfWidth) )
continue;
615 float fpathverylong[2] = {6.0,4.5};
616 float fpathlong[2] = {6.0,4.5};
628 hrecodx[iview]->Fill(vdir.X());
629 hrecody[iview]->Fill(vdir.Y());
630 hrecodz[iview]->Fill(vdir.Z());
643 if(((ftruepath > 3.7) && (ftruepath <= 3.9))
655 if( (
abs(vdir.Y()) < 0.5) && (ftruew > -1000)){
656 htruew[iview]->Fill(ftruew);
661 if( fpath > fpathverylong[iview]){
666 if( (fw > 550.0) && (fpath > fpathlong[iview])){
679 if(
isstopper && xyz[2] > minZ && xyz[2] < maxZ){
684 fLvsdxTot[iview+2]->Fill(vdir.X(),fpath);
685 fLvsdyTot[iview+2]->Fill(vdir.Y(),fpath);
686 fLvsdzTot[iview+2]->Fill(vdir.Z(),fpath);
687 hrecodx[iview+2]->Fill(vdir.X());
688 hrecody[iview+2]->Fill(vdir.Y());
689 hrecodz[iview+2]->Fill(vdir.Z());
703 if(((ftruepath > 3.7) && (ftruepath <= 3.9))
714 if(
abs(vdir.Y()) < 0.5 && (ftruew > -1000.0)){
715 htruew[iview+2]->Fill(ftruew);
716 hrecow[iview+2]->Fill(fw);
719 if(fpath > fpathverylong[iview]){
724 if((fw > 550.0) && (fpath > fpathlong[iview])){
743 std::set<TVector3,
bool(*)(TVector3
const&, TVector3
const&)>
const& trajPoints)
745 float minZ = trajPoints.begin()->Z();
746 float maxZ = trajPoints.rbegin()->Z();
750 minZ -= cellHalfDepth;
751 maxZ += cellHalfDepth;
752 float xi = trajPoints.begin()->X();
753 float yi = trajPoints.begin()->Y();
754 float zi = trajPoints.begin()->Z();
755 float xf = trajPoints.rbegin()->X();
756 float yf = trajPoints.rbegin()->Y();
757 float zf = trajPoints.rbegin()->Z();
758 float dr =
sqrt((xi-xf)*(xi-xf)+(yi-yf)*(yi-yf)+(zi-zf)*(zi-zf));
760 float dx = (xf-xi)/dr;
761 float dy = (yf-yi)/dr;
762 float dz = (zf-zi)/dr;
765 std::pair<uint32_t, uint32_t>
pc;
767 std::set<double> planeZBounds;
770 if(boundaryMap.size() < 1)
return;
773 for(
size_t c = 0;
c < track.
NCell(); ++
c){
778 float fpoissonlambda = 0.0;
779 float ftruemev = 0.0;
780 float ftruepath = -10.0;
781 float ftruew = -9999.0;
787 fpoissonlambda += pe.PoissonLambda();
792 for(
size_t n = 0;
n < ides.size(); ++
n)
793 ftruemev += ides[
n].
energy*1000;
799 auto particle = parts.front();
800 std::vector<sim::FLSHit> flsHits =
fBT->
ParticleToFLSHit(particle->TrackId(), particle->PdgCode());
801 for(
auto fls : flsHits){
802 if(fls.GetPlaneID() == track.
Cell(
c)->
Plane() &&
803 fls.GetCellID() == track.
Cell(
c)->
Cell() ){
804 ftruepath = fls.GetTotalPathLength();
805 ftruew = fls.GetZAverage();
811 if((ftruemev == 0.0) || (fpoissonlambda == 0.0))
continue;
829 for(
unsigned int tp=0;tp<ntp;tp++){
848 for(
unsigned int tp=0;tp<ntp;tp++){
858 vdir = (vtp1 - vtp0).
Unit();
861 pc.first = ch->
Plane();
862 pc.second = ch->
Cell();
864 point = TVector3(rh.
X(), rh.
Y(), rh.
Z());
866 double fadc = ch->
ADC();
867 double fpe = ch->
PE();
869 double fpecorr = -5.0;
877 if((ch->
View() ==
geo::kX) && (vdir.X() != 0 || vdir.Z() != 0)){
880 (xyz[0]-cellwidth), (xyz[0]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
883 fw = (xyz[2]-vtp0.Z())*vdir.Y()/vdir.Z() + vtp0.Y();
886 fw = (xyz[2]-zi)*dy/dz + yi;
890 if((ch->
View() ==
geo::kY) && (vdir.Y() != 0 || vdir.Z() != 0)){
893 (xyz[1]-cellwidth), (xyz[1]+cellwidth), (xyz[2]-cellHalfDepth), (xyz[2]+cellHalfDepth));
896 fw = (xyz[2]-vtp0.Z())*vdir.X()/vdir.Z() + vtp0.X();
899 fw = (xyz[2]-zi)*dx/dz + xi;
903 if((fpath < 1.0) || (fpath > 25.0) || (
abs(fw) > detHalfWidth) )
continue;
914 float fpathverylong[2] = {6.0,4.5};
915 float fpathlong[2] = {6.0,4.5};
938 bool istruepath =
false;
944 bool isanglerestr =
false;
946 if( (
abs(vdir.Y()) < 0.5) && (ftruew > -1000)){
953 if( fpath > fpathverylong[iview]){
958 bool islongnear =
false;
960 if( (fw > 550.0) && (fpath > fpathlong[iview])){
965 bool isMIPstopper =
false;
967 if(
isstopper && xyz[2] > minZ && xyz[2] < maxZ){
983 std::set<TVector3,
bool(*)(TVector3
const&, TVector3
const&)>
const& trajPoints)
985 float minZ = trajPoints.begin()->Z();
986 float maxZ = trajPoints.rbegin()->Z();
989 minZ -= cellHalfDepth;
990 maxZ += cellHalfDepth;
992 std::pair<uint32_t, uint32_t>
pc;
994 std::set<double> planeZBounds;
997 if(boundaryMap.size() < 1)
return;
998 for(
size_t c = 0;
c < track.
NCell(); ++
c){
1000 pc.first = ch->
Plane();
1001 pc.second = ch->
Cell();
1003 point = TVector3(rh.
X(), rh.
Y(), rh.
Z());
1007 double xyz[3]= {0.};
1009 if(
isstopper && xyz[2] > minZ && xyz[2] < maxZ)
1019 float xlow,
float xhigh,
float zlow,
float zhigh)
1022 float lambda[4] = {0.0};
1023 float otherview[4] = {0.0};
1024 float sign[4] = {0.0};
1026 if((x0 - xlow)*(x0 - xhigh) < 0){
1029 else if((x0 - xlow)*(x0 - xhigh) == 0){
1037 if((z0 - zlow)*(z0 - zhigh) < 0){
1040 else if((z0 - zlow)*(z0 - zhigh) == 0){
1049 lambda[0] = (xlow - x0)/dx;
1050 lambda[1] = (xhigh - x0)/dx;
1051 lambda[2] = (zlow - z0)/dz;
1052 lambda[3] = (zhigh - z0)/dz;
1053 otherview[0] = lambda[0]*dz + z0;
1054 otherview[1] = lambda[1]*dz + z0;
1055 otherview[2] = lambda[2]*dx + x0;
1056 otherview[3] = lambda[3]*dx + x0;
1057 sign[0] = (otherview[0] - zlow)*(otherview[0] - zhigh);
1058 sign[1] = (otherview[1] - zlow)*(otherview[1] - zhigh);
1059 sign[2] = (otherview[2] - xlow)*(otherview[2] - xhigh);
1060 sign[3] = (otherview[3] - xlow)*(otherview[3] - xhigh);
1062 if((sign[0] > 0) && (sign[1] > 0) &&
1063 (sign[2] > 0) && (sign[3] > 0)){
1066 else if((sign[0] == 0) || (sign[1] == 0) || (sign[2] == 0) || (sign[3] == 0)){
1069 else if((sign[0] < 0) && (sign[1] < 0)){
1072 else if((sign[2] < 0) && (sign[3] < 0)){
TH2 * ftrueThredVsWTot[4]
T max(const caf::Proxy< T > &a, T b)
size_t NTrajectoryPoints() const
void FillTree(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &) > const &trajPoints)
void FillHist(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &) > const &trajPoints)
back track the reconstruction to the simulation
virtual void analyze(art::Event const &e)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
StopperThreshold(fhicl::ParameterSet const &p)
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
zBoundMap MakeZBoundaryMap(std::set< double > const &planeZBounds, std::vector< TVector3 > const &trajectory)
Return a map of the z position of each trajectory point on a track to the bounding positions of where...
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
unsigned short Plane() const
double DistToEdgeXY(TVector3 vertex)
std::string fTrackLabel
label of module creating the tracks used
void testPath(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &) > const &trajPoints)
const CellGeo * Cell(int icell) const
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
Vertical planes which measure X.
virtual ~StopperThreshold()
A rb::Prong with full reconstructed trajectory.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
art::ServiceHandle< geo::LiveGeometry > fLivegeo
unsigned short Cell() const
virtual void endRun(art::Run const &)
TH2F * fADCratevsYlongTot[4]
double DistToEdgeZ(TVector3 vertex)
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
T get(std::string const &key) const
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
TH2 * ftrueShadowVsYTot[4]
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
double fMaxDistFromEnd
maximum distance from end of track to use a point
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
void FindTrajectoryPoints(rb::Track const &track, std::set< TVector3, bool(*)(TVector3 const &, TVector3 const &)> &trajPoints)
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
T * make(ARGS...args) const
double DetHalfWidth() const
int16_t ADC(uint32_t i) const
TH2 * fLvsWnothreshTot[4]
virtual void beginRun(art::Run const &)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
bool ltTVec(TVector3 const &a, TVector3 const &b)
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
bool cellSortZ(art::Ptr< rb::CellHit > const &a, art::Ptr< rb::CellHit > const &b)
art::ServiceHandle< cheat::BackTracker > fBT
TVector3 Stop() const
Position of the final trajectory point.
void FindZBoundaries(std::set< double > &planeZBounds)
Find the boundaries in the z direction of planes in the detector.
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
double fMinDistFromEnd
minimum distance from end of track to use a point
virtual void reconfigure(fhicl::ParameterSet const &p)
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
TH2F * fADCratevsWlongTot[4]
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.
int CalculateWall(float x0, float z0, float dx, float dz, float xlow, float xhigh, float zlow, float zhigh)
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
Return the path length of a track in the cell in question.