Subversion Repositories f9daq

Rev

Rev 71 | Rev 84 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 71 Rev 73
Line 58... Line 58...
58 
 
 58 
 
59 
void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}
 59 
void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}
60 
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
 60 
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
61 
*/
 61 
 */
62 
//-----------------------------------------------------------------------------
 62 
//-----------------------------------------------------------------------------
63 
void Init(double rsys_scale = 7.0)
 63 
void Init(double rsys_scale = 10.0)
64 
{      
 64 
{      
65 
        RTSetStyle(gStyle);
 65 
  RTSetStyle(gStyle);
66 
        gStyle->SetOptStat(10);
 66 
  gStyle->SetOptStat(10);
67 
       
 67 
 
68 
        if(show_3d) {
 68 
  if(show_3d) {
69 
                double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
 69 
      double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
70 
                double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
 70 
      double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
71 
               
 71 
 
72 
                c3dview = (TCanvas*)gROOT->FindObject("c3dview");
 72 
      c3dview = (TCanvas*)gROOT->FindObject("c3dview");
73 
                if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);      c3dview->SetFillColor(0);}
 73 
      if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);
- 
 
 74 
      c3dview->SetFillColor(0);}
74 
                else c3dview->Clear();
 75 
      else c3dview->Clear();
75 
                       
 76 
 
76 
                TView3D *view = new TView3D(1, rsys0, rsys1);
 77 
      TView3D *view = new TView3D(1, rsys0, rsys1);
77 
                view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
 78 
      view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
78 
                //view->SetPerspective();
 79 
      //view->SetPerspective();
Line 175... Line 176...
175 
}
 176 
}
176 
//-----------------------------------------------------------------------------
 177 
//-----------------------------------------------------------------------------
177 
       
 178 
 
178 
//-----------------------------------------------------------------------------
 179 
//-----------------------------------------------------------------------------
179 
// en zarek
 180 
// en zarek
- 
 
 181 
double Single(CDetector *detector, DetectorParameters& parameters,
180 
double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
 182 
TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
181 
{
 183 
{
182 
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
 - 
 
183 
        theta = Pi()*theta/180.0;
 184 
  theta = Pi()*theta/180.0;
184 
        if(theta < 1e-6) theta = 1e-6;
 185 
  if(theta < 1e-6) theta = 1e-6;
185 
        phi = phi*Pi()/180.0;
 186 
  phi = phi*Pi()/180.0;
186 
        if(phi < 1e-6) phi = 1e-6;
 187 
  if(phi < 1e-6) phi = 1e-6;
187 
       
 188 
 
188 
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 189 
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
189 
       
 190 
 
190 
        //detector->Init();
 191 
  //detector->Init();
191 
        if(show_3d) detector->Draw(draw_width);
 192 
  if(show_3d) detector->Draw(draw_width);
192 
       
 193 
 
193 
        CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
 194 
  CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
194 
                          TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
 195 
      TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi),
- 
 
 196 
      TMath::Sin(theta)*TMath::Cos(phi));
- 
 
 197 
  // Set z-polarization == vertical
- 
 
 198 
  TVector3 polarization(0,0,1);
- 
 
 199 
  polarization.RotateY(-theta);
- 
 
 200 
  polarization.RotateZ(phi);
- 
 
 201 
  if (polarization.Dot(ray0->GetK()) > 1e-5) printf("ERROR: pol not perep\n");
- 
 
 202 
  ray0->setPolarization(polarization);
195 
        CRay *ray1 = new CRay;
 203 
  CRay *ray1 = new CRay;
196 
       
 204 
 
197 
        detector->Propagate(*ray0, ray1, show_3d);
 205 
  detector->Propagate(*ray0, ray1, show_3d);
- 
 
 206 
 
- 
 
 207 
  CRay *incidentPolarization = new CRay;
- 
 
 208 
  incidentPolarization->SetColor(kGreen);
- 
 
 209 
  TVector3 drawPosition = ray0->GetR();
- 
 
 210 
  drawPosition.SetX(drawPosition.X() - 5);
- 
 
 211 
  incidentPolarization->Set(drawPosition, polarization);
- 
 
 212 
  incidentPolarization->DrawS(drawPosition.X(), 1);
- 
 
 213 
 
- 
 
 214 
  TVector3 outPolarization = ray1->GetP();
- 
 
 215 
  drawPosition = ray1->GetR();
- 
 
 216 
  drawPosition.SetX(drawPosition.X() + 5);
- 
 
 217 
  CRay* rayPol = new CRay(drawPosition, outPolarization);
- 
 
 218 
  rayPol->SetColor(kBlack);
- 
 
 219 
  rayPol->DrawS(drawPosition.X(), 1);
198 
       
 220 
 
199 
        delete ray0;
 221 
  delete ray0;
200 
        delete ray1;
 222 
  delete ray1;
201 
       
 223 
 
202 
        return (detector->GetHActive())->GetEntries() / (double)(1);
 224 
  return (detector->GetHActive())->GetEntries() / (double)(1);
203 
}
 225 
}
204 
//-----------------------------------------------------------------------------
 226 
//-----------------------------------------------------------------------------
205 
// zarki, razporejeni v mrezi
 227 
// zarki, razporejeni v mrezi
206 
double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
 228 
double Grid(CDetector *detector, DetectorParameters& parameters,
- 
 
 229 
int NN = 10, double theta = 0.0)
207 
{
 230 
{
208 
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
 231 
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
209 
        theta = Pi()*theta/180.0;
 232 
  theta = Pi()*theta/180.0;
210 
        if(theta < 1e-6) theta = 1e-6;
 233 
  if(theta < 1e-6) theta = 1e-6;
211 
       
 234 
 
212 
        //detector->Init();
 235 
  //detector->Init();
213 
        if(show_3d) detector->Draw(draw_width);
 236 
  if(show_3d) detector->Draw(draw_width);
214 
 
 237 
 
215 
        const double b = parameters.getB();
 238 
  const double b = parameters.getB();
216 
       
 239 
 
217 
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 240 
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
218 
       
 241 
 
219 
        for(int i = 0; i < NN; i++) {
 242 
  for(int i = 0; i < NN; i++) {
220 
                for(int j = 0; j < NN; j++) {
 243 
      for(int j = 0; j < NN; j++) {
221 
                CRay *ray0 = new CRay(CENTER.x() - offset,
 244 
          CRay *ray0 = new CRay(CENTER.x() - offset,
222 
                              0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
 245 
              0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
223 
                              0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
 246 
              0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
224 
                                          TMath::Cos(theta),
 247 
              TMath::Cos(theta),
225 
                                          0.0,
 248 
              0.0,
226 
                                          TMath::Sin(theta));
 249 
              TMath::Sin(theta));
- 
 
 250 
          TVector3 polarization(0, 1, 0);
- 
 
 251 
          polarization.RotateY(-theta);
- 
 
 252 
          ray0->setPolarization(polarization);
227 
                CRay *ray1 = new CRay;
 253 
          CRay *ray1 = new CRay;
228 
                if(i == (NN/2))
 254 
          if(i == (NN/2))
229 
                        detector->Propagate(*ray0, ray1, show_3d);
 255 
            detector->Propagate(*ray0, ray1, show_3d);
230 
                else
 256 
          else
231 
                        detector->Propagate(*ray0, ray1, 0);
 257 
            detector->Propagate(*ray0, ray1, 0);
Line 245... Line 271...
245 
        return acceptance;
 271 
  return acceptance;
246 
}
 272 
}
247 
//-----------------------------------------------------------------------------
 273 
//-----------------------------------------------------------------------------
248 
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
 274 
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
249 
// vsi pod kotom (theta, phi)
 275 
// vsi pod kotom (theta, phi)
250 
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phi, int show_rays)
 276 
double RandYZ(CDetector *detector, DetectorParameters& parameters,
- 
 
 277 
int NN, double theta, double phi, int show_rays)
251 
{
 278 
{
252 
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
 279 
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
253 
        theta = theta*3.14159265358979312/180.0;
 280 
  theta = theta*3.14159265358979312/180.0;
254 
        if(theta < MARGIN) theta = MARGIN;
 281 
  if(theta < MARGIN) theta = MARGIN;
255 
        phi = phi*3.14159265358979312/180.0;
 282 
  phi = phi*3.14159265358979312/180.0;
256 
        if(phi < MARGIN) phi = MARGIN;
 283 
  if(phi < MARGIN) phi = MARGIN;
257 
       
 284 
 
258 
        TDatime now;
 285 
  TDatime now;
Line 265... Line 292...
265 
        double SiPM = parameters.getA();
 292 
  double SiPM = parameters.getA();
266 
        double    M = parameters.getM();
 293 
  double    M = parameters.getM();
267 
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 294 
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
268 
       
 295 
 
269 
        for(int i = 0; i < NN; i++) {
 296 
  for(int i = 0; i < NN; i++) {
270 
               
 297 
 
271 
                double rx = CENTER.x() - offset;
 298 
      double rx = CENTER.x() - offset;
272 
                double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
 299 
      double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
273 
                double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
 300 
      double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
274 
                               
 301 
 
275 
                double rl = TMath::Cos(theta);
 302 
      double rl = TMath::Cos(theta);
276 
                double rm = TMath::Sin(theta)*TMath::Sin(phi);
 303 
      double rm = TMath::Sin(theta)*TMath::Sin(phi);
277 
                double rn = TMath::Sin(theta)*TMath::Cos(phi); 
 304 
      double rn = TMath::Sin(theta)*TMath::Cos(phi);
278 
               
 305 
 
279 
                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
 306 
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
- 
 
 307 
      TVector3 polarization(0, 0, 1);
- 
 
 308 
      polarization.RotateY(-theta);
- 
 
 309 
      polarization.RotateZ(phi);
- 
 
 310 
      ray0->setPolarization(polarization);
280 
                CRay *ray1 = new CRay;
 311 
      CRay *ray1 = new CRay;
281 
               
 312 
 
282 
                if(i < show_rays)
 313 
      if(i < show_rays)
283 
                        detector->Propagate(*ray0, ray1, show_3d);
 314 
        detector->Propagate(*ray0, ray1, show_3d);
284 
                else
 315 
      else
285 
                        detector->Propagate(*ray0, ray1, 0);
 316 
        detector->Propagate(*ray0, ray1, 0);
286 
                //delete ray0;
 317 
      delete ray1;
287 
                //delete ray1; 
318 
      delete ray0;
288 
        }
 319 
  }
289 
       
 320 
 
290 
        double acceptance = 0.0;
 321 
  double acceptance = 0.0;
291 
        /*
 322 
  /*
292 
        if( !(parameters.getPlateOn()) )
 323 
        if( !(parameters.getPlateOn()) )
293 
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
 324 
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
294 
  else
 325 
  else
Line 302... Line 333...
302 
// = nakljucni vstopni polozaj in kot
 333 
// = nakljucni vstopni polozaj in kot
303 
// NN - number of rays to be simulated with angles theta distributed uniformly
 334 
// NN - number of rays to be simulated with angles theta distributed uniformly
304 
//      in cos theta and phi uniformly:
 335 
//      in cos theta and phi uniformly:
305 
//      theta [0, theta]
 336 
//      theta [0, theta]
306 
//      phi [0,360]
 337 
//      phi [0,360]
- 
 
 338 
double RandIso(CDetector *detector, DetectorParameters& parameters,
307 
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
 339 
int NN = 1e3, double thetaMin=0.0, double theta = 30.0, int show_rays = 30, int show_rand = 0)
308 
{
 340 
{
309 
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
 341 
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
310 
        double pi = 3.14159265358979312;
 342 
  double pi = 3.14159265358979312;
311 
        theta = theta*3.14159265358979312/180.0;
 343 
  theta = theta*3.14159265358979312/180.0;
312 
        if(theta < MARGIN) theta = MARGIN;
 344 
  if(theta < MARGIN) theta = MARGIN;
Line 317... Line 349...
317 
        double rx, ry, rz, rl, rm, rn;
 349 
  double rx, ry, rz, rl, rm, rn;
318 
        double rphi, rtheta;
 350 
  double rphi, rtheta;
319 
        //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
 351 
  //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
320 
        double theta_min_rad = TMath::Cos(theta);
 352 
  double theta_min_rad = TMath::Cos(theta);
321 
        double theta_max_rad = TMath::Cos(thetaMin);
 353 
  double theta_max_rad = TMath::Cos(thetaMin);
322 
 
 354 
 
323 
          TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
 355 
  TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
324 
          hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
 356 
  hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
325 
          hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    
 357 
  hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);
326 
          htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
 358 
  htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
327 
          htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     
 359 
  htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);
328 
          hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
 360 
  hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
329 
          hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       
 361 
  hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);
Line 341... Line 373...
341 
        double SiPM = parameters.getA();
 373 
  double SiPM = parameters.getA();
342 
        double M = parameters.getM();
 374 
  double M = parameters.getM();
343 
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
 375 
  double offset = (parameters._plateOn ? parameters._plateWidth : 0);
344 
       
 376 
 
345 
        for(int i = 0; i < NN; i++) {
 377 
  for(int i = 0; i < NN; i++) {
346 
               
 378 
 
347 
                rx = CENTER.x() - offset;
 379 
      rx = CENTER.x() - offset;
348 
                ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
 380 
      ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
349 
                rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
 381 
      rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
350 
                               
 382 
 
351 
                rphi = rand.Uniform(0.0, 2.0*pi);
 383 
      rphi = rand.Uniform(0.0, 2.0*pi);
Line 353... Line 385...
353 
                //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
 385 
      //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
354 
                rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
 386 
      rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
355 
               
 387 
 
356 
                rl = TMath::Cos(rtheta);
 388 
      rl = TMath::Cos(rtheta);
357 
                rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
 389 
      rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
358 
                rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      
 390 
      rn = TMath::Sin(rtheta)*TMath::Cos(rphi);
359 
               
 391 
 
360 
                if(show_rand) {
 392 
      if(show_rand) {
361 
                        htheta->Fill(rtheta);
 393 
          htheta->Fill(rtheta);
362 
                        hcostheta->Fill( TMath::Cos(rtheta) );
 394 
          hcostheta->Fill( TMath::Cos(rtheta) );
363 
                        hphi->Fill(rphi);
 395 
          hphi->Fill(rphi);
364 
                        hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
 396 
          hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
Line 367... Line 399...
367 
                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
 399 
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
368 
                CRay *ray1 = new CRay;
 400 
      CRay *ray1 = new CRay;
369 
               
 401 
 
370 
                if(i < show_rays) {
 402 
      if(i < show_rays) {
371 
                        detector->Propagate(*ray0, ray1, show_3d);
 403 
          detector->Propagate(*ray0, ray1, show_3d);
372 
                        }
 404 
      }
373 
                else {
 405 
      else {
374 
                        detector->Propagate(*ray0, ray1, 0);
 406 
          detector->Propagate(*ray0, ray1, 0);
375 
                        }
 407 
      }
376 
               
 408 
 
377 
          delete ray0;
 409 
      delete ray0;
378 
          delete ray1; 
 410 
      delete ray1;
379 
        }
 411 
  }
380 
       
 412 
 
381 
        if(show_rand) {
 413 
  if(show_rand) {
Line 388... Line 420...
388 
                c2rand->cd(3); htheta->Draw();
 420 
      c2rand->cd(3); htheta->Draw();
389 
                c2rand->cd(5); hcostheta->Draw();
 421 
      c2rand->cd(5); hcostheta->Draw();
390 
                c2rand->cd(2); hl->Draw();
 422 
      c2rand->cd(2); hl->Draw();
391 
                c2rand->cd(4); hm->Draw();
 423 
      c2rand->cd(4); hm->Draw();
392 
                c2rand->cd(6); hn->Draw();
 424 
      c2rand->cd(6); hn->Draw();
393 
        }
 425 
  }
394 
       
 426 
 
395 
        double acceptance = 0.0;
 427 
  double acceptance = 0.0;
396 
        /*
 428 
  /*
397 
        if( !(parameters.getPlateOn()) )
 429 
        if( !(parameters.getPlateOn()) )
398 
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
 430 
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
Line 403... Line 435...
403 
        return acceptance;
 435 
  return acceptance;
404 
}
 436 
}
405 
 
 437 
 
406 
// Beamtest distribution
 438 
// Beamtest distribution
407 
// with fixed theta and phi in interval phiMin, phiMax
 439 
// with fixed theta and phi in interval phiMin, phiMax
- 
 
 440 
double beamtest(CDetector *detector, DetectorParameters& parameters,
408 
double beamtest(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
 441 
int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
409 
{
 442 
{
410 
        double pi = 3.14159265358979312;
 443 
  double pi = 3.14159265358979312;
411 
        theta *= pi/180.0;
 444 
  theta *= pi/180.0;
412 
        if(theta < MARGIN) theta = MARGIN;
 445 
  if(theta < MARGIN) theta = MARGIN;
413 
        phiMin *= pi/180.0;
 446 
  phiMin *= pi/180.0;
Line 416... Line 449...
416 
        TDatime now;
 449 
  TDatime now;
417 
        TRandom rand;
 450 
  TRandom rand;
418 
        rand.SetSeed(now.Get());
 451 
  rand.SetSeed(now.Get());
419 
        double rx, ry, rz, rl, rm, rn;
 452 
  double rx, ry, rz, rl, rm, rn;
420 
        double rphi;
 453 
  double rphi;
421 
        //double rtheta;
 - 
 
422 
        //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
 - 
 
423 
        //double theta_min_rad = TMath::Cos(theta);
 - 
 
424 
 
 454 
 
425 
        TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
 455 
  TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
426 
        hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
 456 
  hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
427 
        hphi = new TH1F("hphi", "hphi", 100, -pi, pi);  
 457 
  hphi = new TH1F("hphi", "hphi", 100, -pi, pi);
428 
        htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
 458 
  htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
Line 459... Line 489...
459 
               
 489 
 
460 
                if(show_rand) {
 490 
      if(show_rand) {
461 
                        //htheta->Fill(rtheta);
 491 
          //htheta->Fill(rtheta);
462 
                        //hcostheta->Fill( TMath::Cos(rtheta) );
 492 
          //hcostheta->Fill( TMath::Cos(rtheta) );
463 
                        hphi->Fill(rphi);
 493 
          hphi->Fill(rphi);
- 
 
 494 
          hl->Fill(rl);
- 
 
 495 
          hm->Fill(rm);
464 
                        hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
 496 
          hn->Fill(rn);
465 
                }
 497 
      }
466 
               
 498 
 
467 
                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
 499 
      CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
- 
 
 500 
      // inital polarizaton
- 
 
 501 
      TVector3 polarization(0, 0, 1);
- 
 
 502 
      polarization.RotateY(-theta);
- 
 
 503 
      polarization.RotateX(rphi);
- 
 
 504 
      ray0->setPolarization(polarization);
468 
                CRay *ray1 = new CRay;
 505 
      CRay *ray1 = new CRay;
469 
               
 506 
 
470 
                if(i < show_rays) {
 507 
      if(i < show_rays) {
471 
                        detector->Propagate(*ray0, ray1, show_3d);
 508 
          detector->Propagate(*ray0, ray1, show_3d);
472 
                        }
 509 
      }