CerenkovResponse_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CerenkovResponse
3 // Module Type: analyzer
4 // File: CerenkovResponse_module.cc
5 //
6 // Generated at Tue Apr 9 23:39:58 2013 by Adam Aurisano using artmod
7 ////////////////////////////////////////////////////////////////////////
8 
19 
20 #include "Geometry/Geometry.h"
21 #include "Geometry/LiveGeometry.h"
22 
24 #include "GeometryObjects/Geo.h"
27 
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/Cluster.h"
31 #include "RecoBase/Track.h"
32 #include "Utilities/func/MathUtil.h"
37 #include "NovaDAQConventions/DAQConventions.h"
38 #include "cetlib/search_path.h"
39 #include "RawData/RawTrigger.h"
40 
41 #include "MCCheater/BackTracker.h"
42 #include "Simulation/FLSHit.h"
43 
44 #include <cmath>
45 #include <iostream>
46 #include <iomanip>
47 #include <map>
48 #include <unordered_map>
49 #include <vector>
50 #include <string>
51 #include <set>
52 #include <utility>
53 #include <algorithm>
54 
55 #include "TFile.h"
56 #include "TTree.h"
57 #include "TH1D.h"
58 #include "TH2F.h"
59 #include "TH2D.h"
60 #include "TH3D.h"
61 #include "TProfile.h"
62 #include "TString.h"
63 #include "TVector3.h"
64 
65 #include <sys/stat.h>
66 
67 #include <iostream>
68 #include <iomanip>
69 
70 class CerenkovResponse;
71 
73 public:
74  explicit CerenkovResponse(fhicl::ParameterSet const & p);
75  virtual ~CerenkovResponse();
76 
77  void analyze(art::Event const & e) override;
78 
79  void beginJob() override;
80  void beginRun(const art::Run& run) override;
81  void beginSubRun(const art::SubRun& sr) override;
82  void endSubRun(const art::SubRun& sr) override;
83  void endRun(const art::Run& run) override;
84 
85  void CreateHistos();
86 
87 private:
88  // Declare member data here.
89 
90  TTree* fCerenkovTree;
91 
92  Float_t _pathLength;
93  Float_t _pe;
94  Int_t _plane;
95  Int_t _cell;
96  Float_t _distFromEnd;
97  Float_t _trackLength;
98  Float_t _flightLength;
99  Bool_t _isX;
100  Bool_t _isProton;
101  Bool_t _isMuon;
102  Float_t _cerenkovPho;
103  Float_t _eBirks;
104  Float_t _w;
105 
107 
110 
111 };
112 
114  EDAnalyzer(p)
115 {}
116 
118 {}
119 
121 {
122  fFirstTime = true;
123 }
124 
126 {
127 }
128 
130 {
131 }
132 
133 
135 {
136 }
137 
139 {
140  if (!fFirstTime) return;
141 
144 
145  fCerenkovTree = tfs->make<TTree>("cerenkovTree","Cerenkov Information");
146  fCerenkovTree->Branch("pathLength", &_pathLength, "pathLength/F");
147  fCerenkovTree->Branch("pe", &_pe, "pe/F");
148  fCerenkovTree->Branch("distFromEnd", &_distFromEnd, "distFromEnd/F");
149  fCerenkovTree->Branch("trackLength", &_trackLength, "trackLength/F");
150  fCerenkovTree->Branch("flightLength",&_flightLength,"flightLength/F");
151  fCerenkovTree->Branch("isX", &_isX, "isX/B");
152  fCerenkovTree->Branch("isProton", &_isProton, "isProton/B");
153  fCerenkovTree->Branch("isMuons", &_isMuon, "isMuon/B");
154  fCerenkovTree->Branch("cerenkovPho", &_cerenkovPho, "cerenkovPho/F");
155  fCerenkovTree->Branch("eBirks", &_eBirks, "eBirks/F");
156  fCerenkovTree->Branch("w", &_w, "w/F");
157 
158 // fQualLabels.push_back("pathQualXY");
159  fQualLabel = "pathQualTraj";
160  fFirstTime = false;
161  return;
162 }
163 
165 {
166  CreateHistos();
167  return;
168 }
169 
171 {
172 
173  std::string trackLabel = "cosmictrack";
174  std::string pclistLabel = "pclist";
175 
177 
178  //if there are no hits, this event was probably filtered
179  //TODO: ??????????
181  evt.getByLabel("calhit", hits);
182 
183  if (hits.failedToGet())
184  {
185 // std::cerr << "Failed to get calhits" << std::endl;
186 // return;
187  }
188 
190 
192  evt.getByLabel(trackLabel, tracks);
193  if (tracks.failedToGet()) {
194  std::cerr << "Failed to get cosmic tracks" << std::endl;
195  return;
196  }
197  if (tracks->empty()) {
198  std::cerr << "Empty cosmic tracks container" << std::endl;
199  return;
200  }
201 
202  art::FindManyP<caldp::PCHit> fmpcxy(tracks, evt, art::InputTag(pclistLabel, fQualLabel.c_str()) );
203  if (!fmpcxy.isValid()) {
204  std::cerr << "No associated pchits for cosmic tracks!" << std::endl;
205  return;
206  }
207 
208  for(unsigned int i = 0; i < tracks->size(); ++i){
209  const rb::Track* track(&tracks->at(i));
210  //3D tracks
211  if(track->View() != geo::kXorY) continue;
212 
213  //Determine if tracks come from outside and stop
215 
216  double distStartToFront = livegeo->DistToFront(track->Start());
217  double distStartToTop = livegeo->DistToTop(track->Start());
218  double distStartToSide = livegeo->DistToEdgeX(track->Start());
219 
220  //physical starts
221 
222  if (distStartToFront < 0) continue;
223  if (distStartToTop < 0) continue;
224  if (distStartToSide < 10) continue;
225  // endpoint isn't contained
226 
227  double distEndToEdge = std::min(livegeo->DistToEdgeXY(track->Stop()), livegeo->DistToEdgeZ(track->Stop()));
228 
229  bool isStopper = true;
230  if (distEndToEdge < 50) isStopper = false;
231  if (!isStopper) {continue;}
232 
233  const std::vector< art::Ptr<caldp::PCHit> > pchits = fmpcxy.at(i);
234 
235  for (auto& pchit : pchits) {
236  _pathLength = pchit->Path();
237  _pe = pchit->PE();
238  _plane = pchit->Plane();
239  _cell = pchit->Cell();
240  _distFromEnd = track->TotalLength() - pchit->FlightLen();
241  _trackLength = track->TotalLength();
242  _flightLength = pchit->FlightLen();
243  _w = pchit->W();
244  _isProton = false;
245  _isMuon = false;
246  _cerenkovPho = 0;
247  _eBirks = 0;
248 
249  if (pchit->View() == geo::kX) _isX = true;
250  else _isX = false;
251 
252  bool isLowEff = false;
253  //remove potentially low efficiency cells and low oil cells
254  int modcell = _cell%32;
255  if ( modcell == 15 || modcell == 31) isLowEff = true;
256  if ( !_isX && modcell == 30) isLowEff = true;
257  bool inMC = false;
258  if (geo->DetId() == novadaq::cnv::kNEARDET) {
259  if (_plane >= int(geo->FirstPlaneInMuonCatcher())) inMC = true;
260  }
261  if (isLowEff || inMC) continue;
262 
263  if (fBT->HaveTruthInfo()) {
264  for (unsigned int iCell = 0; iCell < track->NCell(); ++iCell) {
265  art::Ptr<rb::CellHit> chit = track->Cell(iCell);
266  if (chit->Plane() == _plane && chit->Cell() == _cell) {
267  const sim::Particle* part = fBT->HitToParticle(chit);
268  if(!part) {std::cout<<(track->Cell(iCell))->PE()<<std::endl;continue;};
269 
270  if (part->PdgCode() == 13) _isMuon = true;
271  else if (part->PdgCode() == 2212) _isProton = true;
272  std::vector<sim::FLSHit> flsHits = fBT->HitToFLSHit(chit);
273  if(flsHits.empty()) continue;
274  for(auto& fls: flsHits ) {
275  _cerenkovPho += fls.GetNCerenkov();
276  _eBirks += fls.GetEdepBirks();
277  }
278  break;
279  }
280  }
281  }
282  fCerenkovTree->Fill();
283  }//done pclist and filling
284  }
285 }
286 
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
void analyze(art::Event const &e) override
X or Y views.
Definition: PlaneGeo.h:30
unsigned short Plane() const
Definition: CellHit.h:39
double DistToEdgeXY(TVector3 vertex)
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
OStream cerr
Definition: OStream.cxx:7
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
double DistToEdgeX(TVector3 vertex)
double DistToFront(TVector3 vertex)
DEFINE_ART_MODULE(TestTMapFile)
const sim::Particle * HitToParticle(art::Ptr< rb::CellHit > const &hit, bool quiet=false) const
Returns the sim::Particle object contributing the most light to the specified rb::CellHit.
Definition: Run.h:31
unsigned short Cell() const
Definition: CellHit.h:40
double DistToTop(TVector3 vertex)
TString part[npart]
Definition: Style.C:32
void hits()
Definition: readHits.C:15
void beginRun(const art::Run &run) override
double DistToEdgeZ(TVector3 vertex)
art::ServiceHandle< cheat::BackTracker > fBT
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
void endRun(const art::Run &run) override
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
caf::StandardRecord * sr
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
CerenkovResponse(fhicl::ParameterSet const &p)
void beginSubRun(const art::SubRun &sr) override
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
T * make(ARGS...args) const
void endSubRun(const art::SubRun &sr) override
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
Helper for AttenCurve.
Definition: Path.h:10
Float_t track
Definition: plot.C:35
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
enum BeamMode string