10 #define CosmicTrends_h 27 #include "TStopwatch.h" 83 unsigned int const& qualityLevel,
84 std::map<
unsigned short, std::set<unsigned short> >
const& planesAndCells,
127 int numTriCellTracks = 0;
131 for(
size_t t = 0;
t < th->size(); ++
t){
137 std::map<unsigned short, std::set<unsigned short> > planesAndCells;
141 double useableCells[2] = {0.};
142 double totalCells[2] = {0.};
144 for(
size_t h = 0;
h < track.
NCell(); ++
h)
147 useableCells[0] = 0.;
148 useableCells[1] = 0.;
152 for(
auto const& pacitr : planesAndCells){
154 unsigned short plane = pacitr.first;
157 std::set<unsigned short> cellSet = pacitr.second;
168 for(
auto citr : cellSet ){
169 if( cellSet.find(citr - 1) != cellSet.end() &&
170 cellSet.find(citr + 1) != cellSet.end() ){
194 double usefulFrac = 0.;
195 if(totalCells[0] + totalCells[1] > 0)
196 usefulFrac = (useableCells[0] + useableCells[1])/(totalCells[0] + totalCells[1]);
202 if(totalCells[0] > 0){
207 if(totalCells[1] > 0){
244 double maxCells = (yCells > xCells) ? yCells : xCells;
249 if(q ==
kTriCell) qualString =
"TriCell";
252 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"stopperCountX"+qualString).c_str(),
254 (
int)totalPlanes, 0., totalPlanes,
255 (
int)maxCells, 0., maxCells));
258 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"totalCountX"+qualString).c_str(),
260 (
int)totalPlanes, 0., totalPlanes,
261 (
int)maxCells, 0., maxCells));
264 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"stopperCountY"+qualString).c_str(),
266 (
int)totalPlanes, 0., totalPlanes,
267 (
int)maxCells, 0., maxCells));
269 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"totalCountY"+qualString).c_str(),
271 (
int)totalPlanes, 0., totalPlanes,
272 (
int)maxCells, 0., maxCells));
275 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"cosThetaVsAz"+qualString).c_str(),
276 ";cos(#theta);#phi (degrees)",
281 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"cosThetaVsAzTru"+qualString).c_str(),
282 ";True cos(#theta);True #phi (degrees)",
287 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"cosThetaVsPlanes"+qualString).c_str(),
288 ";cos(#theta);Planes Crossed",
290 (
int)totalPlanes, 0., totalPlanes));
293 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"azimuthVsPlanes"+qualString).c_str(),
294 ";#phi (degrees);Planes Crossed",
296 (
int)totalPlanes, 0., totalPlanes));
299 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"deltaAzVsAz"+qualString).c_str(),
300 ";#phi (degrees);#Delta#phi (degrees)",
305 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"deltaZenVsZen"+qualString).c_str(),
306 ";cos(#theta);#Deltacos(#theta)",
311 fOneDHists[q].push_back(tfs->
make<TH1D>((
"fractionUsableX"+qualString).c_str(),
312 ";Usable Fraction;Number of Muons",
316 fOneDHists[q].push_back(tfs->
make<TH1D>((
"fractionUsableY"+qualString).c_str(),
317 ";Usable Fraction;Number of Muons",
321 fOneDHists[q].push_back(tfs->
make<TH1D>((
"recoCosTheta"+qualString).c_str(),
322 ";cos(#theta);Number of Muons",
326 fOneDHists[q].push_back(tfs->
make<TH1D>((
"recoAzimuth"+qualString).c_str(),
327 ";#phi (degrees);Number of Muons",
331 fOneDHists[q].push_back(tfs->
make<TH1D>((
"planesCrossed"+qualString).c_str(),
332 ";Planes Crossed;Number of Muons",
333 (
int)totalPlanes*0.2, 0., (
int)totalPlanes));
336 fOneDHists[q].push_back(tfs->
make<TH1D>((
"dcosZ"+qualString).c_str(),
337 ";dz/ds;Number of Tracks",
341 fOneDHists[q].push_back(tfs->
make<TH1D>((
"dcosX"+qualString).c_str(),
342 ";dx/ds;Number of Tracks",
346 fOneDHists[q].push_back(tfs->
make<TH1D>((
"dcosY"+qualString).c_str(),
347 ";dy/ds;Number of TriCell Tracks",
351 fOneDHists[q].push_back(tfs->
make<TH1D>((
"numTracks"+qualString).c_str(),
352 ";Number of TriCell Tracks;Events",
356 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"totalHitCountX"+qualString).c_str(),
358 (
int)totalPlanes, 0., totalPlanes,
359 (
int)maxCells, 0., maxCells));
361 fTwoDHists[q].push_back(tfs->
make<TH2D>((
"totalHitCountY"+qualString).c_str(),
363 (
int)totalPlanes, 0., totalPlanes,
364 (
int)maxCells, 0., maxCells));
377 fMinPlanes = p.
get<
unsigned int >(
"MinimumPlanesPerTrack", 8);
386 unsigned int const& qualityLevel,
387 std::map<
unsigned short, std::set<unsigned short> >
const& planesAndCells,
395 TVector3 startDir = track.
Dir();
396 double cosTheta = -startDir.y()/startDir.Mag();
397 double azimuth =
std::atan2(startDir.x()/startDir.Mag(), -startDir.z()/startDir.Mag()) * 180./TMath::Pi();
398 if(azimuth < 0) azimuth += 360.;
417 LOG_WARNING(
"CosmicTrends") <<
"no simulated particle matches current track, skip it";
419 float trueCosTheta = -truth->
Py()/truth->
P();
420 float trueAzimuth =
std::atan2(truth->
Px()/truth->
P(), -truth->
Pz()/truth->
P()) * 180./TMath::Pi();
421 if(trueAzimuth < 0) trueAzimuth += 360.;
429 TVector3 startPoint(track.
Start());
430 TVector3 endPoint(track.
Stop());
431 bool stopper =
false;
438 for(
auto const& pacitr : planesAndCells){
440 unsigned short plane = pacitr.first;
443 std::set<unsigned short> cellSet = pacitr.second;
448 for(
auto citr : cellSet ){
back track the reconstruction to the simulation
double fContainmentCut
#cm in from the detector boundary to be called a stopper
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
virtual void analyze(art::Event const &e)
double Py(const int i=0) const
unsigned short Plane() const
art::ServiceHandle< geo::Geometry > fGeo
std::map< unsigned int, std::vector< TH2D * > > fTwoDHists
all 2D hists to make
Vertical planes which measure X.
unsigned int Ncells() const
Number of cells in this plane.
double Px(const int i=0) const
A rb::Prong with full reconstructed trajectory.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
std::map< unsigned int, std::vector< TH1D * > > fOneDHists
all 1D hists to make
unsigned short Cell() const
View_t View() const
Which coordinate does this plane measure.
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
double P(const int i=0) const
T get(std::string const &key) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
double fTriCellFrac
minimum fraction of tri-cells in a track
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
#define LOG_WARNING(category)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
T * make(ARGS...args) const
double DetHalfWidth() const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::string fTrackLabel
label of module creating the tracks we use
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
void FillHistograms(rb::Track const &track, unsigned int const &qualityLevel, std::map< unsigned short, std::set< unsigned short > > const &planesAndCells, bool mc)
double Pz(const int i=0) const
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
virtual void reconfigure(fhicl::ParameterSet const &p)
unsigned int NPlanes() const
TVector3 Stop() const
Position of the final trajectory point.
virtual void beginRun(art::Run const &)
Encapsulate the geometry of one entire detector (near, far, ndos)
CosmicTrends(fhicl::ParameterSet const &p)
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.