//g++ particleGun.cxx -o particleGun `root-config --glibs --cflags` #include #include #include #include #include #include"ConfigFile.h" #include "TRandom.h" #include "TApplication.h" #include "TCanvas.h" #include "TH1.h" #include "TVector3.h" #include "TF1.h" #include "TMath.h" #include "TPolyLine3D.h" std::map confI; std::map confD; ConfigFile cf; TRandom myRand(0); void readConfI(std::string var) { if(!(cf.readInto(confI[var] , var) )) { std::cerr << "Reading parameter " << var << " from conf file failed ->abort" << std::endl; throw; } } void readConfD(std::string var) { if(!(cf.readInto(confD[var] , var) )) { std::cerr << "Reading parameter " << var << " from conf file failed ->abort" << std::endl; throw; } } void getXY(double z,const TVector3& pos, const TVector3& mom,double& x,double& y){ /* X=P+t*M z-component: z=Pz+t*Mz <=> t=(z-Pz)/Mz */ double t=(z-pos.Z())/mom.Z(); x=pos.X()+t*mom.X(); y=pos.Y()+t*mom.Y(); } void getXZ(double y,const TVector3& pos, const TVector3& mom,double& x,double& z){ /* X=P+t*M y-component: y=Py+t*My <=> t=(y-Py)/My */ double t=(y-pos.Y())/mom.Y(); x=pos.X()+t*mom.X(); z=pos.Z()+t*mom.Z(); } void getYZ(double x,const TVector3& pos, const TVector3& mom,double& y,double& z){ /* X=P+t*M x-component: x=Px+t*Mx <=> t=(x-Px)/Mx */ double t=(x-pos.X())/mom.X(); y=pos.Y()+t*mom.Y(); z=pos.Z()+t*mom.Z(); } bool acceptance(const TVector3& pos, const TVector3& mom){ double x,y,z; //check for hit in bottom scintillator getXY(confD["BOTZ"],pos,mom,x,y); if(xconfD["BOTXMAX"]) return false; if(yconfD["BOTYMAX"]) return false; //check for penetration of upper and lower 'endcap' of readout volume getXY(confD["TPCZMIN"],pos,mom,x,y); if(x>confD["TPCXMIN"] && xconfD["TPCYMIN"] && yconfD["TPCXMIN"] && xconfD["TPCYMIN"] && yconfD["TPCYMIN"] && yconfD["TPCZMIN"] && zconfD["TPCYMIN"] && yconfD["TPCZMIN"] && zconfD["TPCXMIN"] && xconfD["TPCZMIN"] && zconfD["TPCXMIN"] && xconfD["TPCZMIN"] && zDraw(); TH1D *hTL = new TH1D("hTL","hTL",100,0.,13.); TH1I *hNCL = new TH1I("hNCL","hNCL",30,0,30); TH1D *hZ = new TH1D("hZ","hZ",100,0.,8.); static const double avDistance=.5; char buf[50]; sprintf(buf,"1/%f*exp(-x/%f)",avDistance,avDistance); TF1* distF = new TF1("distF",buf,0.,1.E2); for(int iev=0;ievconfD["TPCXMIN"]&&x1zconfD["TPCYMIN"]&&y1zconfD["TPCXMIN"]&&x2zconfD["TPCYMIN"]&&y2zconfD["TPCYMIN"]&&y1xconfD["TPCZMIN"]&&z1xconfD["TPCYMIN"]&&y2xconfD["TPCZMIN"]&&z2xconfD["TPCXMIN"]&&x1yconfD["TPCZMIN"]&&z1yconfD["TPCXMIN"]&&x2yconfD["TPCZMIN"]&&z2y vecs; TVector3 p1; if(p1x){p1.SetXYZ(confD["TPCXMIN"],y1x,z1x);vecs.push_back(p1);} if(p2x){p1.SetXYZ(confD["TPCXMAX"],y2x,z2x);vecs.push_back(p1);} if(p1y){p1.SetXYZ(x1y,confD["TPCYMIN"],z1y);vecs.push_back(p1);} if(p2y){p1.SetXYZ(x2y,confD["TPCYMAX"],z2y);vecs.push_back(p1);} if(p1z){p1.SetXYZ(x1z,y1z,confD["TPCZMIN"]);vecs.push_back(p1);} if(p2z){p1.SetXYZ(x2z,y2z,confD["TPCZMAX"]);vecs.push_back(p1);} assert(vecs.size()==2); trackLength=(vecs[0]-vecs[1]).Mag(); // std::cout << trackLength << std::endl << "===============" << std::endl; hTL->Fill(trackLength); if(vecs[0].z() i){ nClus=i+2; } } hNCL->Fill(nClus); if(nClus>1){ distance=trackLength/(nClus-1); std::cout << avDistance << " " << distance << std::endl; //TAG hZ->Fill(vecs[0].z()); TVector3 point=vecs[0]; TVector3 direction=mom; direction.Print(); TVector3 v=vecs[1]-vecs[0]; v.SetMag(1.); v.Print(); direction.SetMag(1.); for(int i=0;iFill(point.z()); } } } new TCanvas("c1","c1"); hTL->Draw(); new TCanvas("c2","c2"); // hZ->Draw(); for(int i=0;i<1000000;++i){ // hZ->Fill(distF->GetRandom()); double r = gRandom->Uniform(); } hZ->Draw(); new TCanvas("c3","c3"); hNCL->Draw(); //distF->Draw(); app->Run(); }