EWCosmics_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: EWCosmics
3 // Module Type: analyzer
4 // File: EWCosmics_module.cc
5 ////////////////////////////////////////////////////////////////////////
6 
17 #include "fhiclcpp/ParameterSet.h"
19 #include "Geometry/Geometry.h"
20 #include "RecoBase/Track.h"
21 #include "RecoBase/FitSum.h"
22 
23 #include <iostream>
24 #include <cassert>
25 #include <vector>
26 #include <fstream>
27 #include <string>
28 #include <map>
29 #include "TMath.h"
30 #include "TF1.h"
31 #include "TH2D.h"
32 #include "TGraph.h"
33 #include "TTree.h"
34 #include "TGeoManager.h"
35 
36 #include "RecoBase/RecoHit.h"
37 #include "RecoBase/Cluster.h"
38 #include "MCCheater/BackTracker.h"
39 #include "Simulation/FLSHit.h"
40 #include "Simulation/FLSHitList.h"
42 
43 #include "RecoBase/FilterList.h"
44 #include "RecoBase/Cluster.h"
45 #include "RecoBase/Vertex.h"
46 #include "RecoBase/Prong.h"
47 #include "Utilities/AssociationUtil.h"
49 
50 #include "CLHEP/Random/RandFlat.h"
52 
53 namespace red {
54  class MatterProbe;
55  class EWCosmics;
56  class Energystimatter;
58  class EnergystimatterVar;
59 }
60 
62  public:
64  virtual double EforX(double)=0;
65 };
66 
68  public:
70  void SetS(double);
71  double EforX(double);
72  private:
73  double S;
74 };
75 
77  S=alpha;
78 }
79 
81  return S*x;
82 }
83 
85  public:
87  private:
88 };
89 
91  public:
93  double mattersum(TVector3 startp, TVector3 startd);
94  void SetGeoManager(TGeoManager* gm){fGeoManager=gm;}
95  void SetMap();
96  void ECalc(TVector3 startp, TVector3 startd);
97  double GetDelE();
98  private:
99  TGeoManager* fGeoManager = nullptr;
100  std::map<std::string,TGraph> tab,tabT;
101  double deltaE;
102 };
103 
104 double red::MatterProbe::mattersum(TVector3 startp, TVector3 startd){
105  fGeoManager->InitTrack(startp.X(),startp.Y(),startp.Z(),-startd.X(),-startd.Y(),-startd.Z());
106  double step,dens,sum=0;
107  TGeoNode *nextnode;
108  nextnode = fGeoManager->FindNextBoundaryAndStep();
109  while(nextnode){
110  step = fGeoManager->GetStep();
111  if (step>1E10) break;
112 
113  //check that we are not in the detector volume anymore
114  TString path=fGeoManager->GetPath();
115  if(!path.Contains("/vDetEnclosure_")){
116  TGeoVolume *vol;
117  vol = fGeoManager->GetCurrentVolume();
118  TGeoMaterial *mat;
119  mat = vol->GetMaterial();
120  dens = mat->GetDensity();
121  sum += step*dens;
122  }
123  nextnode = fGeoManager->FindNextBoundaryAndStep();
124  }
125  return sum;
126 }
127 
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);
144  for(auto it : tab){
145  tabT[it.first]=TGraph(it.second.GetN(),it.second.GetY(),it.second.GetX());
146  }
147 
148 }
149 
150 void red::MatterProbe::ECalc(TVector3 startp, TVector3 startd){
151  TGraph *tg,*tgT;
153  tg =tfs2->make<TGraph>();
154  tgT=tfs2->make<TGraph>();
155 
156  fGeoManager->InitTrack(startp.X(),startp.Y(),startp.Z(),-startd.X(),-startd.Y(),-startd.Z());
157  double step,dens,cX;
158  TGeoNode *nextnode;
159  nextnode = fGeoManager->FindNextBoundaryAndStep();
160  double dE=0;
161  while(nextnode){
162  step = fGeoManager->GetStep();
163  if (step>1E10) break;
164 
165  //check that we are not in the detector volume anymore
166  TString path=fGeoManager->GetPath();
167  if(!path.Contains("/vDetEnclosure_")){
168 
169  TGeoVolume *vol;
170  vol = fGeoManager->GetCurrentVolume();
171  TGeoMaterial *mat;
172  mat = vol->GetMaterial();
173  dens= mat->GetDensity();
174  cX = step*dens;
175 
176  std::string matname=mat->GetName();
177  auto it=tab.find(matname);
178  if(it!=tab.end()){
179  *tg=TGraph(it->second);
180  cX+=tg->Eval(dE);
181  }
182  auto itT=tabT.find(matname);
183  if(itT!=tabT.end()){
184  *tgT=TGraph(itT->second);
185  dE=tgT->Eval(cX);
186  }
187  }
188  nextnode = fGeoManager->FindNextBoundaryAndStep();
189  }
190  deltaE=dE/1000.;
191 }
192 
194  return deltaE;
195 }
196 
198  public:
199  explicit EWCosmics(fhicl::ParameterSet const & p);
200  EWCosmics(EWCosmics const &) = delete;
201  EWCosmics(EWCosmics &&) = delete;
202  EWCosmics & operator = (EWCosmics const &) = delete;
203  EWCosmics & operator = (EWCosmics &&) = delete;
204 
205  void analyze(art::Event const & e) override;
206 
207  void beginJob() override;
208  void endJob() override;
209 
210  double distToBorder(const TVector3& p);
211 
212  void FirstWay(TVector3, double&, double&);
213  void SecondWay(const TVector3&, double&, double&);
214  private:
217 
225  double theta,phi,matter,Esurfv,Esurfc,th_d,ph_d,Ecal,EBPF,Erec;
226  unsigned long Ne=0; // Event number
227  TVector3 point,start,stop,startdir;
228  TH2D *hEarth;
229  TTree *tMu;
230 
231  const double phitoN=0;//26;
232  const double FDmaxZ=5962.;//5959.08;//geo->DetLength();
233  const double FDminZ=0.;
234  const double FDmaxX=782.95;//geo->DetHalfWidth();
235  const double FDminX=-FDmaxX;
236  const double FDmaxY=774.7;//geo->DetHalfHeight();
237  const double FDminY=-FDmaxY;
238 };
239 
241  EDAnalyzer(p) // ,
242  // More initializers here.
243 {
244  std::cout<<"EWCosmicslyzer starting"<<std::endl;
245 
246  fSliceLabel=p.get<std::string>("SliceLabel");
247  fTrackLabel=p.get<std::string>("TrackLabel");
248  fFitSumLabel=p.get<std::string>("FitSumLabel");
249  fVertexLabel=p.get<std::string>("VertexLabel");
250  fClusterLabel=p.get<std::string>("ClusterLabel");
251  fClusterInstance=p.get<std::string>("ClusterInstance");
252 }
253 
256  // hEarth=tfs->make<TH2D>("hEarth","Cosmic tracks; #phi;#theta",360,0,360,90,0,90);
257  tMu=tfs->make<TTree>("tMu","Stopped cosmic mu tracks");
258  //tMu->Branch("theta",&theta,"theta/D");
259  //tMu->Branch("phi",&phi,"phi/D");
260  tMu->Branch("th_d",&th_d,"th_d/D");
261  tMu->Branch("ph_d",&ph_d,"ph_d/D");
262  tMu->Branch("start",&start);
263  tMu->Branch("stop",&stop);
264  tMu->Branch("Ne",&Ne,"Ne/l");
265  tMu->Branch("Esurfv",&Esurfv, "Esurfv/D");
266  //tMu->Branch("Esurfc",&Esurfc, "Esurfc/D");
267  //tMu->Branch("Ecal",&Ecal, "Ecal/D");
268  tMu->Branch("EBPF",&EBPF, "EBPF/D");
269 // tMu->Branch("Erec",&Erec, "Erec/D");
270 }
271 
273  std::cout<<"bye-bye!"<<std::endl;
274 }
275 
276 double red::EWCosmics::distToBorder(const TVector3& p){
277  double dx=std::max(std::min(p.X()-FDminX,FDmaxX-p.X()),0.);
278  double dy=std::max(std::min(p.Y()-FDminY,FDmaxY-p.Y()),0.);
279  double dz=std::max(std::min(p.Z()-FDminZ,FDmaxZ-p.Z()),0.);
280  return std::min(dz,std::min(dx,dy));
281 }
282 
283 void red::EWCosmics::FirstWay(TVector3 v, double& T, double& P){
284  //Okay kid, this is where it gets complicated
285  //We have Z pointing along the detector, X pointing left, and Y pointing up
286  //We want X pointing along the detector, and Z pointing down
287  //Let's rotate the vector
288  v.RotateX(-TMath::PiOver2()); //now for vector new Y is old Z, Z is opposite old Y
289  v.RotateZ(-TMath::PiOver2()); //now for vector the newest Y is opposite the first X. Phi is going clockwise
290  double theta=v.Theta()*TMath::RadToDeg();
291  //theta=180-theta; //because Z pointing down, but
292  //we are interested in angles where CR comes FROM, not where TO, so
293  //theta=180-theta; //again. Just ignore both operations
294  double phi=v.Phi()*TMath::RadToDeg();
295  //We are interested in angles where CR comes FROM, not where TO
296  phi=-phi;
297  //Phi is going clockwise. 0 is correspond particles coming along the detector from the North, 90 more or less from East
298  phi-=phitoN; //now 0-particles are strict from the North, 90 from East, 270 from West
299  //But azimuth should be counted from the South
300  phi+=180; //now 180 is correspond particles coming from the North, 270 from East, 90(+360) from West. Oh, shut up
301  while(phi< 0){phi+=360;}
302  while(phi>=360){phi-=360;}
303  T=theta;
304  P=phi;
305 }
306 
307 void red::EWCosmics::SecondWay(const TVector3& v, double& T, double& P){
308  //Les just substitute axes in formulas
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();
311  //We are interested in angles where CR comes FROM, not where TO
312  theta=180-theta;
313  phi=-phi;
314  //Phi is going clockwise. 0 is correspond particles coming along the detector from the North, 90 more or less from East
315  phi-=phitoN; //now 0-particles are strict from the North, 90 from East, 270 from West
316  //But azimuth should be counted from the South
317  phi+=180; //now 180 is correspond particles coming from the North, 270 from East, 90(+360) from West. Oh, shut up
318  while(phi< 0){phi+=360;}
319  while(phi>=360){phi-=360;}
320  T=theta;
321  P=phi;
322 }
323 
325  Ne = e.id().event();//event number
326 
328  auto geomngr = geo->ROOTGeoManager();
329 
330  fMatterProber.SetGeoManager(geomngr);
331  //fEMC.SetS(0.002);
332 
334 
336  e.getByLabel(fSliceLabel, slice);
337  if((!slice.failedToGet())&&(!slice->empty())){
338  art::FindManyP<rb::Track> trackh(slice, e, fTrackLabel);
340 
341  for(size_t i = 0; i < trackh.size(); ++i) {
342  // Find the breakpoint tracks associated to this slice
343  std::vector< art::Ptr<rb::Track> > tracks = trackh.at(i);
344  // Only look at slices with one track
345  if (tracks.size()!=1)continue;
346 
347  const class rb::Track t=*tracks.at(0);
348 
349  start =t.Start();
350  stop =t.Stop();
351  startdir=t.Dir();
352 
353  if((red::EWCosmics::distToBorder(stop)<100)||(red::EWCosmics::distToBorder(start)>50)||(t.NXCell()<5)||(t.NYCell()<5)){continue;}
355  if(theta>90)continue;
356 
357  th_d=startdir.Theta()*TMath::RadToDeg();
358  ph_d=startdir.Phi()*TMath::RadToDeg();
359  while(ph_d< 0){ph_d+=360;}
360  while(ph_d>=360){ph_d-=360;}
361 
364  matter=overcalc->traceBack(t);
365  std::cout<<"Matter = "<<matter<<std::endl;
366  //Esurfc=fEMC.EforX(matter);
367 
368  art::FindManyP<rb::FitSum> kininf(tracks, e, fFitSumLabel);
369  std::vector<art::Ptr<rb::FitSum>> kinin = kininf.at(0);
370  art::Ptr<rb::FitSum> kin = kinin.at(0); //track kinetic information
371  TLorentzVector mom = kin->FourMom();
372  EBPF = mom.E();
373 
374 // Erec=Esurfv+EBPF;
375  tMu->Fill();}
376  }
377 }
378 
T max(const caf::Proxy< T > &a, T b)
back track the reconstruction to the simulation
const double phitoN
double traceBack(const rb::Track &t) const
set< int >::iterator it
#define S(x, n)
red::MatterProbe fMatterProber
std::string fFitSumLabel
void analyze(art::Event const &e) override
const char * p
Definition: xmltok.h:285
TGeoManager * ROOTGeoManager() const
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
red::EnergystimatterConst fEMC
std::string tracksLabel
double distToBorder(const TVector3 &p)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
double dE
void SetGeoManager(TGeoManager *gm)
const double FDmaxY
void ECalc(TVector3 startp, TVector3 startd)
#define P(a, b, c, d, e, x)
const double FDminY
std::string fSliceLabel
const double FDmaxX
double dy[NP][NC]
double dx[NP][NC]
std::string fClusterLabel
void endJob() override
void beginJob()
double dz[NP][NC]
Encapsulate the geometry of one entire detector (near, far, ndos)
const TLorentzVector FourMom() const
Definition: FitSum.h:72
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void FirstWay(TVector3, double &, double &)
Vertex location in position and time.
EWCosmics(fhicl::ParameterSet const &p)
std::map< std::string, TGraph > tabT
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
const std::string path
Definition: plot_BEN.C:43
const double FDminX
EventNumber_t event() const
Definition: EventID.h:116
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
std::string fClusterInstance
void beginJob() override
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
void SecondWay(const TVector3 &, double &, double &)
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
const double FDmaxZ
double T
Definition: Xdiff_gwt.C:5
const double FDminZ
std::string fVertexLabel
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
Helper for AttenCurve.
Definition: Path.h:10
double mattersum(TVector3 startp, TVector3 startd)
Definition: fwd.h:28
TGeoManager * gm
Definition: make_fe_box.C:4
EventID id() const
Definition: Event.h:56
std::string fTrackLabel
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
Eigen::MatrixXd mat
unsigned long Ne
virtual double EforX(double)=0