26 #include "TLorentzVector.h" 91 const int TrackEva::GAMMAID = 22;
92 const int TrackEva::PIZEROID = 111;
93 const int TrackEva::PIPLUSID = 211;
94 const int TrackEva::PIMINUSID = -211;
95 const int TrackEva::NUEID = 12;
96 const int TrackEva::MUONID = 13;
97 const int TrackEva::NUMUID = 14;
98 const int TrackEva::NUTAUID = 16;
99 const int TrackEva::PROTONID = 2212;
137 MUON_MC = tfs->
make<TNtuple>(
"MUON_MC",
"muon_mc_info",
138 "energy:trackLength:Dir_Cos:ID");
141 "energy:trackLength:Dir_Cos:ID");
143 MUON_RC = tfs->
make<TNtuple>(
"MUON_RC",
"muon_rc_info",
144 "tracklength_diff:energy_diff:Angle:TL_MC:E_MC:Dir_Cos:ID");
147 "tracklength_diff:energy_diff:Angle:TL_MC:E_MC:Dir_Cos:ID");
150 track_Nb = tfs->
make<TH1I>(
"track_Nb",
"Total Reconstructed Track-Number",
153 proton_Nb = tfs->
make<TH1I>(
"proton_Nb",
"Total Proton Number",
155 proton_Tr = tfs->
make<TH1I>(
"proton_Tr",
"Reco Proton Track Number",
158 muon_Nb = tfs->
make<TH1I>(
"muon_Nb",
"Total Muon Number",
160 muon_Tr = tfs->
make<TH1I>(
"muon_Tr",
"Reco Muon Track Number",
167 std::vector<double> Mu_E, P_E, LMC_mu, LMC_p;
168 std::vector<TVector3> P_Momenta, Mu_Momenta;
170 TVector3 momentum_P, momentum_Mu;
183 int NbTrack = trackhandle->size();
188 std::vector<int>
MuonID, ProtonID;
195 momentum_P.SetXYZ(p->
Px(),
198 P_Momenta.push_back(momentum_P);
199 P_E.push_back(p->
E());
200 ProtonID.push_back(p->
TrackId());
204 momentum_Mu.SetXYZ(p->
Px(),
207 Mu_Momenta.push_back(momentum_Mu);
208 Mu_E.push_back( p->
E());
209 MuonID.push_back(p->
TrackId());
214 for (
unsigned int n=0;
n!=P_E.size();
n++){
217 P_Momenta[n].
Z()/P_Momenta[n].
Mag(),
222 for (
unsigned int n=0;
n!=Mu_E.size();++
n){
225 Mu_Momenta[n].
Z()/Mu_Momenta[n].
Mag(),
240 for (
int i=0;
i!=NbTrack;
i++){
242 int cell_Nb = track->
NCell();
245 for (
unsigned int j=0;
j!=MuonID.size();
j++){
248 for (
int idx = 0;
idx != cell_Nb;
idx++){
250 bool real_hit=
false;
254 for (
unsigned int pts=0; pts!=trackides.size(); pts++){
255 if (trackides[pts].
trackID==MuonID[
j]) real_hit=
true;
258 if (real_hit) cell_count++;
260 if ((cell_count+0.0)/(cell_Nb+0.0)>=
fMatching){
265 Mu_Momenta[
j].Angle(track->
Dir()),
268 Mu_Momenta[j].
Z()/Mu_Momenta[
j].Mag(),
276 for (
unsigned int j=0;
j<ProtonID.size();
j++){
281 bool real_hit=
false;
283 for (
unsigned int pts=0; pts!=trackides.size(); pts++){
284 if (trackides[pts].
trackID==ProtonID[
j]) real_hit=
true;
287 if (real_hit) cell_count++;
289 if ((cell_count+0.0)/(cell_Nb+0.0)>
fMatching){
294 P_Momenta[
j].Angle(track->
Dir()),
297 P_Momenta[j].
Z()/P_Momenta[
j].Mag(),
318 TLorentzVector pos1 =
it->first;
319 TLorentzVector pos2 = (
it+1)->first;
322 if(pos2.Z()>9.96604 && pos2.Z()<
fEdgeZ){
323 length+=(pos1.Vect()-pos2.Vect()).
Mag();
double E(const int i=0) const
This is the version 2, works for multi-track event.
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
static const int PIZEROID
double Py(const int i=0) const
fvar< T > fabs(const fvar< T > &x)
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
const simb::MCTrajectory & Trajectory() const
list_type::const_iterator const_iterator
double Px(const int i=0) const
static const int PROTONID
double Tracklength(const simb::MCTrajectory &trajectory)
void reconfigure(fhicl::ParameterSet const &pset)
DEFINE_ART_MODULE(TestTMapFile)
virtual double TotalLength() const
Length (cm) of all the track segments.
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
list_type::const_iterator const_iterator
T get(std::string const &key) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
void analyze(const art::Event &evt)
This class describes a particle created in the detector Monte Carlo simulation.
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
static const int PIPLUSID
T * make(ARGS...args) const
float MuonID(const caf::SRProxy *sr)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
double Pz(const int i=0) const
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
static const int PIMINUSID