Back to home page

Project CMSSW displayed by LXR

 
 

    


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 }