Warning, /GeneratorInterface/CosmicMuonGenerator/test/geo_el is written in an unsupported language. File is not indexed.
0001 {
0002 cout << "macro geo_el: check geometry (and energy loss)" << endl;
0003 gROOT->Reset();
0004 gROOT->LoadMacro("GeneratorInterface/CosmicMuonGenerator/src/SingleParticleEvent.cc");
0005
0006 double PlugVx = PlugOnShaftVx;
0007 double PlugVz = PlugOnShaftVz;
0008 //double PlugVz = -33000.; //[mm]
0009
0010 int colVec[6] = {0000,2222, 8888, 1111, 9999, 4444};
0011
0012
0013 gStyle->SetPalette(6,colVec);
0014
0015
0016 //////////////////////////////////////////////////////////////////////////////////////
0017 TCanvas *C1 = new TCanvas("densities","Densities",10,10,250,250);
0018 C1->SetBottomMargin(0.15);
0019 C1->Divide(1,1);
0020 C1->cd(1); //gPad->SetTicks(1,1);
0021 TH2F* col = new TH2F("col","col",5, 0.55, 5.55, 5, 0.5, 5.5);
0022 col->Fill(1, 1, 3);
0023 col->Fill(1, 2, 4);
0024 col->Fill(1, 3, 5);
0025 col->Fill(1, 4, 2);
0026 col->Fill(1, 5, 1);
0027 col->SetMinimum(0);
0028 col->SetMaximum(6);
0029 col->Draw("COLAH");
0030
0031 TLatex* line = new TLatex(1.7, 5.0, "Plug #rho=6.3 g/cm^{3}");
0032 line->SetTextSize(0.075);
0033 line->Draw();
0034 TLatex* line = new TLatex(1.7, 4.0, "Concrete #rho=2.5 g/cm^{3}");
0035 line->SetTextSize(0.075);
0036 line->Draw();
0037 TLatex* line = new TLatex(1.7, 3.0, "Rock #rho=2.2 g/cm^{3}");
0038 line->SetTextSize(0.075);
0039 line->Draw();
0040 TLatex* line = new TLatex(1.7, 2.0, "Clay #rho=1.5 g/cm^{3}");
0041 line->SetTextSize(0.075);
0042 line->Draw();
0043 TLatex* line = new TLatex(1.7, 1.0, "Air #rho=0.0 g/cm^{3}");
0044 line->SetTextSize(0.075);
0045 line->Draw();
0046
0047
0048 gPad->Update();
0049
0050
0051
0052 // scan geometry
0053 //////////////////////////////////////////////////////////////////////////////////////
0054 TH2D* zyv = new TH2D("zyv","Z-Y view at X=0",500,-65.,35.,600,-20.,100.);
0055 double vx = 0.;
0056 for (int iz=0; iz<500; ++iz){
0057 double vz = double(iz)*0.2 - 64.9;
0058 for (int iy=0; iy<600; ++iy){
0059 double vy = double(iy)*0.2 - 19.9;
0060 zyv->Fill(float(vz),float(vy), inMat(vx*1000.,vy*1000.,vz*1000., PlugVx, PlugVz));
0061
0062 }
0063 }
0064 TCanvas *C1 = new TCanvas("zy_view","Z-Y view at X=0",50,50,440,550);
0065 C1->SetBottomMargin(0.15);
0066 C1->Divide(1,1);
0067 C1->cd(1); gPad->SetTicks(1,1);
0068 zyv->SetMinimum(0);
0069 zyv->SetMaximum(6);
0070 zyv->Draw("COL");
0071 zyv->GetXaxis()->SetLabelSize(0.035);
0072 zyv->GetXaxis()->SetTitleSize(0.035);
0073 zyv->GetXaxis()->SetTitleOffset(1.5);
0074 zyv->GetXaxis()->SetTitle("Z [m]");
0075 zyv->GetXaxis()->SetLabelOffset(0.001);
0076 zyv->GetYaxis()->SetLabelSize(0.035);
0077 zyv->GetYaxis()->SetTitleSize(0.035);
0078 zyv->GetYaxis()->SetTitleOffset(2.0);
0079 zyv->GetYaxis()->SetTitle("Y [m]");
0080 gPad->Update();
0081
0082 //////////////////////////////////////////////////////////////////////////////////////
0083 TH2D* xyv1 = new TH2D("xyv1","X-Y view at Z=-38",250,-20.,30.,600,-20.,100.);
0084 double vz = -38.227;
0085 for (int ix=0; ix<250; ++ix){
0086 double vx = double(ix)*0.2 - 19.9;
0087 for (int iy=0; iy<600; ++iy){
0088 double vy = double(iy)*0.2 - 19.9;
0089 xyv1->Fill(float(vx),float(vy),inMat(vx*1000.,vy*1000.,vz*1000., PlugVx, PlugVz));
0090 }
0091 }
0092 TCanvas *C2 = new TCanvas("xy_view1","X-Y view at Z=-38",500,50,220,550);
0093 C2->SetBottomMargin(0.15);
0094 C2->SetLeftMargin(0.15);
0095 C2->Divide(1,1);
0096 C2->cd(1); gPad->SetTicks(1,1);
0097 gPad->SetLeftMargin(0.15);
0098 gPad->SetLeftMargin(0.2);
0099 xyv1->SetMinimum(0);
0100 xyv1->SetMaximum(6);
0101 xyv1->Draw("COL");
0102 xyv1->GetXaxis()->SetLabelSize(0.07);
0103 xyv1->GetXaxis()->SetTitleSize(0.055);
0104 xyv1->GetXaxis()->SetTitleOffset(0.8);
0105 xyv1->GetXaxis()->SetTitle("X [m]");
0106 xyv1->GetXaxis()->SetLabelOffset(-0.022);
0107 xyv1->GetYaxis()->SetLabelSize(0.07);
0108 xyv1->GetYaxis()->SetTitleSize(0.055);
0109 xyv1->GetYaxis()->SetTitleOffset(1.9);
0110 xyv1->GetYaxis()->SetTitle("Y [m]");
0111 gPad->Update();
0112
0113 //////////////////////////////////////////////////////////////////////////////////////
0114 TH2D* xyv2 = new TH2D("xyv2","X-Y view at Z=15",300,-20.,40.,600,-20.,100.);
0115 double vz = 15.;
0116 for (int ix=0; ix<300; ++ix){
0117 double vx = double(ix)*0.2 - 19.9;
0118 for (int iy=0; iy<600; ++iy){
0119 double vy = double(iy)*0.2 - 19.9;
0120 xyv2->Fill(float(vx),float(vy),inMat(vx*1000.,vy*1000.,vz*1000., PlugVx, PlugVz));
0121 }
0122 }
0123 TCanvas *C3 = new TCanvas("xy_view2","X-Y view at Z=14",730,50,270,550);
0124 C3->SetBottomMargin(0.15);
0125 C3->Divide(1,1);
0126 C3->cd(1); gPad->SetTicks(1,1);
0127 gPad->SetLeftMargin(0.2);
0128 xyv2->SetMinimum(0);
0129 xyv2->SetMaximum(6);
0130 xyv2->Draw("COL");
0131 xyv2->GetXaxis()->SetLabelSize(0.06);
0132 xyv2->GetXaxis()->SetTitleSize(0.055);
0133 xyv2->GetXaxis()->SetTitleOffset(0.7);
0134 xyv2->GetXaxis()->SetTitle("X [m]");
0135 xyv2->GetXaxis()->SetLabelOffset(-0.022);
0136 xyv2->GetYaxis()->SetLabelSize(0.06);
0137 xyv2->GetYaxis()->SetTitleSize(0.055);
0138 xyv2->GetYaxis()->SetTitleOffset(1.7);
0139 xyv2->GetYaxis()->SetTitle("Y [m]");
0140 gPad->Update();
0141
0142 //////////////////////////////////////////////////////////////////////////////////////
0143 TH2D* xzv1 = new TH2D("xzv1","X-Z view at Y=0",350,-20.,50.,675,-65.,70.);
0144 double vy = 0.;
0145 for (int ix=0; ix<350; ++ix){
0146 double vx = double(ix)*0.2 - 19.9;
0147 for (int iz=0; iz<675; ++iz){
0148 double vz = double(iz)*0.2 - 64.9;
0149 xzv1->Fill(float(vx),float(vz),inMat(vx*1000.,vy*1000.,vz*1000., PlugVx, PlugVz));
0150 }
0151 }
0152 TCanvas *C4 = new TCanvas("xz_view1","X-Z view at Y=0",150,500,250,480);
0153 C4->SetBottomMargin(0.15);
0154 //C4->SetLeftMargin(3.15);
0155 C4->Divide(1,1);
0156 C4->cd(1); gPad->SetTicks(1,1);
0157 gPad->SetLeftMargin(0.2);
0158
0159 xzv1->SetMinimum(0); //needed to get right color interval
0160 xzv1->SetMaximum(6); //needed to get right color interval
0161 xzv1->Draw("COL");
0162 xzv1->GetXaxis()->SetLabelSize(0.07);
0163 xzv1->GetXaxis()->SetTitleSize(0.055);
0164 xzv1->GetXaxis()->SetTitleOffset(0.8);
0165 xzv1->GetXaxis()->SetTitle("X [m]");
0166 xzv1->GetXaxis()->SetLabelOffset(-0.022);
0167 xzv1->GetYaxis()->SetLabelSize(0.07);
0168 xzv1->GetYaxis()->SetTitleSize(0.055);
0169 xzv1->GetYaxis()->SetTitleOffset(1.7);
0170 xzv1->GetYaxis()->SetTitle("Z [m]");
0171 gPad->Update();
0172
0173 //////////////////////////////////////////////////////////////////////////////////////
0174 TH2D* xzv2 = new TH2D("xzv2","X-Z view at Y=50",275,-20.,35.,500,-70.,30.);
0175 double vy = 50.;
0176 for (int ix=0; ix<275; ++ix){
0177 double vx = double(ix)*0.2 - 19.9;
0178 for (int iz=0; iz<500; ++iz){
0179 double vz = double(iz)*0.2 - 69.9;
0180 xzv2->Fill(float(vx),float(vz),inMat(vx*1000.,vy*1000.,vz*1000., PlugVx, PlugVz));
0181 }
0182 }
0183 TCanvas *C5 = new TCanvas("xz_view2","X-Z view at Y=50",360,500,230,460);
0184 C5->SetBottomMargin(0.15);
0185 C5->Divide(1,1);
0186 C5->cd(1); gPad->SetTicks(1,1);
0187 xzv2->SetMinimum(0); //needed to get right color interval
0188 xzv2->SetMaximum(6); //needed to get right color interval
0189 xzv2->Draw("COL");
0190 xzv2->GetXaxis()->SetLabelSize(0.055);
0191 xzv2->GetXaxis()->SetTitleSize(0.055);
0192 xzv2->GetXaxis()->SetTitleOffset(0.4);
0193 xzv2->GetXaxis()->SetTitle("X [m]");
0194 xzv2->GetXaxis()->SetLabelOffset(-0.022);
0195 xzv2->GetYaxis()->SetLabelSize(0.055);
0196 xzv2->GetYaxis()->SetTitleSize(0.055);
0197 xzv2->GetYaxis()->SetTitleOffset(1.3);
0198 xzv2->GetYaxis()->SetTitle("Z [m]");
0199 gPad->Update();
0200
0201 //////////////////////////////////////////////////////////////////////////////////////
0202 TH2D* xzv3 = new TH2D("xzv3","X-Z view at Y=surf.",275,-20.,35.,500,-50.,50.);
0203 double vy = 88.87;
0204 for (int ix=0; ix<275; ++ix){
0205 double vx = double(ix)*0.2 - 19.9;
0206 for (int iz=0; iz<500; ++iz){
0207 double vz = double(iz)*0.2 - 49.9;
0208 xzv3->Fill(float(vx),float(vz),inMat(vx*1000.,vy*1000.,vz*1000., PlugVx, PlugVz));
0209 }
0210 }
0211 TCanvas *C6 = new TCanvas("xz_view3","X-Z view at Y=surf.",600,500,250,480);
0212 C6->SetBottomMargin(0.15);
0213 C6->Divide(1,1);
0214 C6->cd(1); gPad->SetTicks(1,1);
0215 xzv3->SetMinimum(1.);
0216 xzv3->SetMinimum(0.); //needed to get right color interval
0217 xzv3->SetMaximum(5.5); //needed to get right color interval
0218 xzv3->Draw("COL");
0219 xzv3->GetXaxis()->SetLabelSize(0.055);
0220 xzv3->GetXaxis()->SetTitleSize(0.055);
0221 xzv3->GetXaxis()->SetTitleOffset(0.4);
0222 xzv3->GetXaxis()->SetTitle("X [m]");
0223 xzv3->GetXaxis()->SetLabelOffset(-0.022);
0224 xzv3->GetYaxis()->SetLabelSize(0.055);
0225 xzv3->GetYaxis()->SetTitleSize(0.055);
0226 xzv3->GetYaxis()->SetTitleOffset(1.3);
0227 xzv3->GetYaxis()->SetTitle("Z [m]");
0228 gPad->Update();
0229
0230 //////////////////////////////////////////////////////////////////////////////////////
0231 TH2D* xzv4 = new TH2D("xzv3","X-Z view at Y=surf.+1m",200,-20.,20.,450,-45.,45.);
0232 double vy = 88.87+1.00;
0233 for (int ix=0; ix<200; ++ix){
0234 double vx = double(ix)*0.2 - 19.9;
0235 for (int iz=0; iz<450; ++iz){
0236 double vz = double(iz)*0.2 - 44.9;
0237 xzv4->Fill(float(vx),float(vz),inMat(vx*1000.,vy*1000.,vz*1000., PlugVx, PlugVz));
0238 }
0239 }
0240 TCanvas *C7 = new TCanvas("xz_view4","X-Z view at Y=surf.",860,500,180,410);
0241 C7->SetBottomMargin(0.15);
0242 C7->Divide(1,1);
0243 C7->cd(1); gPad->SetTicks(1,1);
0244 xzv4->SetMinimum(0.); //needed to get right color interval
0245 xzv4->SetMaximum(5.5); //needed to get right color interval
0246 xzv4->Draw("COL");
0247 xzv4->GetXaxis()->SetLabelSize(0.055);
0248 xzv4->GetXaxis()->SetTitleSize(0.055);
0249 xzv4->GetXaxis()->SetTitleOffset(0.4);
0250 xzv4->GetXaxis()->SetTitle("X [m]");
0251 xzv4->GetXaxis()->SetLabelOffset(-0.022);
0252 xzv4->GetYaxis()->SetLabelSize(0.055);
0253 xzv4->GetYaxis()->SetTitleSize(0.055);
0254 xzv4->GetYaxis()->SetTitleOffset(1.3);
0255 xzv4->GetYaxis()->SetTitle("Z [m]");
0256 gPad->Update();
0257
0258 //////////////////////////////////////////////////////////////////////////////////////
0259 /*
0260 // check energy loss
0261 SingleParticleEvent TestEvt;
0262 int id = 13;
0263 double E = 100.;
0264 double Px = 0.;
0265 double Py = -sqrt(E*E - MuonMass*MuonMass);
0266 double Pz = 0.;
0267 double Vy = SurfaceOfEarth;
0268 double T0 = 0.;
0269
0270 TH2D* exz = new TH2D("exz","X-Z view of E underground",80,-8.,8.,75,-15.,15.);
0271 for (int ix=0; ix<80; ++ix){
0272 double Vx = double(ix)/5. - 7.9;
0273 for (int iz=0; iz<75; ++iz){
0274 double Vz = double(iz)/2.5 - 14.8;
0275 TestEvt.create(id, Px, Py, Pz, E, MuonMass, (Vx*1000.), Vy, (Vz*1000.), T0);
0276 TestEvt.propagate(1.0, 8000., 15000., false, false);
0277 if (TestEvt.hitTarget()) exz->Fill(float(Vx),float(Vz),float(TestEvt.e()));
0278 }
0279 }
0280 TCanvas *C7 = new TCanvas("el_xz","X-Z of E underground",111,111,400,400);
0281 C7->SetBottomMargin(0.15);
0282 C7->Divide(1,1);
0283 C7->cd(1); gPad->SetTicks(1,1);
0284 exz->Draw("LEGO");
0285 exz->GetXaxis()->SetLabelSize(0.055);
0286 exz->GetXaxis()->SetTitleSize(0.055);
0287 exz->GetXaxis()->SetTitleOffset(1.1);
0288 exz->GetXaxis()->SetTitle("X [m]");
0289 exz->GetYaxis()->SetLabelSize(0.055);
0290 exz->GetYaxis()->SetTitleSize(0.055);
0291 exz->GetYaxis()->SetTitleOffset(1.5);
0292 exz->GetYaxis()->SetTitle("Z [m]");
0293 exz->GetZaxis()->SetLabelSize(0.055);
0294 exz->GetZaxis()->SetTitleSize(0.055);
0295 exz->GetZaxis()->SetTitleOffset(1.2);
0296 exz->GetZaxis()->SetTitle("E_{underground} [GeV]");
0297 gPad->Update();
0298
0299 //////////////////////////////////////////////////////////////////////////////////////
0300 TH2D* eth = new TH2D("eth","E_surf. vs. theta",100,0.,1000.,90,0.,90.);
0301 for (int it=0; it<90; ++it){
0302 double Theta = double(it) + 0.5; // [deg]
0303 double ThRad = Theta*Deg2Rad; // [rad]
0304 double Vz = SurfaceOfEarth*tan(ThRad);
0305 for (int ie=0; ie<100; ++ie){
0306 E = double(ie)*10. + 5.;
0307 double P = sqrt(E*E - MuonMass*MuonMass);
0308 Py = -P*cos(ThRad);
0309 Pz = -P*sin(ThRad);
0310 TestEvt.create(id, Px, Py, Pz, E, MuonMass, 0., Vy, Vz, T0);
0311 TestEvt.propagate(1.0, 8000., 15000., false, false);
0312 if (TestEvt.hitTarget()) eth->Fill(float(E),float(Theta));
0313 }
0314 }
0315 TCanvas *C8 = new TCanvas("e_the","E surface vesrus theta",222,222,400,400);
0316 C8->SetBottomMargin(0.15);
0317 C8->Divide(1,1);
0318 C8->cd(1); gPad->SetTicks(1,1);
0319 eth->Draw("COL");
0320 eth->GetXaxis()->SetLabelSize(0.05);
0321 eth->GetXaxis()->SetTitleSize(0.05);
0322 eth->GetXaxis()->SetTitleOffset(0.9);
0323 eth->GetXaxis()->SetTitle("E_{surface} [GeV]");
0324 eth->GetYaxis()->SetLabelSize(0.05);
0325 eth->GetYaxis()->SetTitleSize(0.05);
0326 eth->GetYaxis()->SetTitleOffset(0.9);
0327 eth->GetYaxis()->SetTitle("#vartheta [#circ]");
0328 gPad->Update();
0329 */
0330 }