34 #include "TGeoManager.h" 47 #include "Utilities/AssociationUtil.h" 64 virtual double EforX(
double)=0;
93 double mattersum(TVector3 startp, TVector3 startd);
96 void ECalc(TVector3 startp, TVector3 startd);
99 TGeoManager* fGeoManager =
nullptr;
100 std::map<std::string,TGraph> tab,
tabT;
105 fGeoManager->InitTrack(startp.X(),startp.Y(),startp.Z(),-startd.X(),-startd.Y(),-startd.Z());
108 nextnode = fGeoManager->FindNextBoundaryAndStep();
110 step = fGeoManager->GetStep();
111 if (step>1E10)
break;
114 TString
path=fGeoManager->GetPath();
115 if(!path.Contains(
"/vDetEnclosure_")){
117 vol = fGeoManager->GetCurrentVolume();
119 mat = vol->GetMaterial();
120 dens = mat->GetDensity();
123 nextnode = fGeoManager->FindNextBoundaryAndStep();
129 const char*
temp=
"%*lg %lg %*lg %*lg %*lg %*lg %*lg %*lg %lg %*lg %*lg";
130 tab[
"Vacuum"] =TGraph(
"EWCosmicAnalyzer/materials/muE_air_dry_1_atm.txt",temp);
131 tab[
"Air"] =TGraph(
"EWCosmicAnalyzer/materials/muE_air_dry_1_atm.txt",temp);
132 tab[
"Water"] =TGraph(
"EWCosmicAnalyzer/materials/muE_water_liquid.txt",temp);
133 tab[
"ShotRock"] =TGraph(
"EWCosmicAnalyzer/materials/muE_standard_rock.txt",temp);
134 tab[
"Granite"] =TGraph(
"EWCosmicAnalyzer/materials/muE_standard_rock.txt",temp);
135 tab[
"Dirt"] =TGraph(
"EWCosmicAnalyzer/materials/muE_standard_rock.txt",temp);
136 tab[
"Concrete"] =TGraph(
"EWCosmicAnalyzer/materials/muE_shielding_concrete.txt",temp);
137 tab[
"BariteRock"] =TGraph(
"EWCosmicAnalyzer/materials/muE_barium_sulfate.txt",temp);
138 tab[
"Scintillator"]=TGraph(
"EWCosmicAnalyzer/materials/muE_paraffin.txt",temp);
139 tab[
"PVC"] =TGraph(
"EWCosmicAnalyzer/materials/muE_polyvinylchloride_PVC.txt",temp);
140 tab[
"WLSFiber"] =TGraph(
"EWCosmicAnalyzer/materials/muE_polyvinylchloride_PVC.txt",temp);
141 tab[
"Glue"] =TGraph(
"EWCosmicAnalyzer/materials/muE_carbon_amorphous_C.txt",temp);
142 tab[
"Steel"] =TGraph(
"EWCosmicAnalyzer/materials/muE_iron_Fe.txt",temp);
143 tab[
"PivoterSteel"]=TGraph(
"EWCosmicAnalyzer/materials/muE_iron_Fe.txt",temp);
145 tabT[
it.first]=TGraph(
it.second.GetN(),
it.second.GetY(),
it.second.GetX());
153 tg =tfs2->
make<TGraph>();
154 tgT=tfs2->
make<TGraph>();
156 fGeoManager->InitTrack(startp.X(),startp.Y(),startp.Z(),-startd.X(),-startd.Y(),-startd.Z());
159 nextnode = fGeoManager->FindNextBoundaryAndStep();
162 step = fGeoManager->GetStep();
163 if (step>1E10)
break;
166 TString
path=fGeoManager->GetPath();
167 if(!path.Contains(
"/vDetEnclosure_")){
170 vol = fGeoManager->GetCurrentVolume();
172 mat = vol->GetMaterial();
173 dens= mat->GetDensity();
177 auto it=tab.find(matname);
179 *tg=TGraph(
it->second);
182 auto itT=tabT.find(matname);
184 *tgT=TGraph(itT->second);
188 nextnode = fGeoManager->FindNextBoundaryAndStep();
208 void endJob()
override;
210 double distToBorder(
const TVector3& p);
212 void FirstWay(TVector3,
double&,
double&);
213 void SecondWay(
const TVector3&,
double&,
double&);
225 double theta,phi,matter,Esurfv,Esurfc,th_d,ph_d,Ecal,EBPF,Erec;
231 const double phitoN=0;
232 const double FDmaxZ=5962.;
233 const double FDminZ=0.;
234 const double FDmaxX=782.95;
235 const double FDminX=-FDmaxX;
236 const double FDmaxY=774.7;
237 const double FDminY=-FDmaxY;
257 tMu=tfs->
make<TTree>(
"tMu",
"Stopped cosmic mu tracks");
260 tMu->Branch(
"th_d",&
th_d,
"th_d/D");
261 tMu->Branch(
"ph_d",&
ph_d,
"ph_d/D");
264 tMu->Branch(
"Ne",&
Ne,
"Ne/l");
265 tMu->Branch(
"Esurfv",&
Esurfv,
"Esurfv/D");
268 tMu->Branch(
"EBPF",&
EBPF,
"EBPF/D");
288 v.RotateX(-TMath::PiOver2());
289 v.RotateZ(-TMath::PiOver2());
290 double theta=v.Theta()*TMath::RadToDeg();
294 double phi=v.Phi()*TMath::RadToDeg();
301 while(phi< 0){phi+=360;}
302 while(phi>=360){phi-=360;}
309 double theta=TMath::ACos(v.Y()/
TMath::Sqrt(v.X()*v.X()+v.Y()*v.Y()+v.Z()*v.Z()))*TMath::RadToDeg();
310 double phi=TMath::ATan2(-v.X(),v.Z())*TMath::RadToDeg();
318 while(phi< 0){phi+=360;}
319 while(phi>=360){phi-=360;}
341 for(
size_t i = 0;
i < trackh.size(); ++
i) {
343 std::vector< art::Ptr<rb::Track> > tracks = trackh.at(
i);
345 if (tracks.size()!=1)
continue;
355 if(
theta>90)
continue;
369 std::vector<art::Ptr<rb::FitSum>> kinin = kininf.at(0);
371 TLorentzVector mom = kin->
FourMom();
T max(const caf::Proxy< T > &a, T b)
back track the reconstruction to the simulation
double traceBack(const rb::Track &t) const
red::MatterProbe fMatterProber
void analyze(art::Event const &e) override
TGeoManager * ROOTGeoManager() const
A rb::Prong with full reconstructed trajectory.
red::EnergystimatterConst fEMC
double distToBorder(const TVector3 &p)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
void SetGeoManager(TGeoManager *gm)
void ECalc(TVector3 startp, TVector3 startd)
#define P(a, b, c, d, e, x)
std::string fClusterLabel
Encapsulate the geometry of one entire detector (near, far, ndos)
const TLorentzVector FourMom() const
T get(std::string const &key) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
void FirstWay(TVector3, double &, double &)
Vertex location in position and time.
EWCosmics(fhicl::ParameterSet const &p)
std::map< std::string, TGraph > tabT
unsigned int NYCell() const
Number of cells in the y-view.
EventNumber_t event() const
T * make(ARGS...args) const
std::string fClusterInstance
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.
void SecondWay(const TVector3 &, double &, double &)
unsigned int NXCell() const
Number of cells in the x-view.
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
double mattersum(TVector3 startp, TVector3 startd)
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual double EforX(double)=0