#define HEX_LADO 354.49/2.0 #define HEX_ALTU 307.00/2.0 #define HEX_VERT 66.4 #define VIDRO_NHSEG 100 #define VIDRO_NVSEG 30 //#define PMT_RAIO 103.5 #define PMT_RAIO 10. #define PMT_HEIGHT 82.9 #define COS30 cos(3.141592654/6.0) #define STEP_Y HEX_ALTU*2.0 #define STEP_X STEP_Y*COS30 #define DET_W1 1152.10 #define DET_W2 1381.5 #define DET_H 900.0 + 101.5 #define OFF_X (DET_W1/2.0 - HEX_LADO) #define OFF_Y1 (DET_W2/2.0 - HEX_ALTU) #define OFF_Y2 (DET_W2/2.0 - STEP_Y) void teste7(Double_t theta0=45., Double_t phi0=120.) { // cria mundo --------------------------------------------------------- TGeoMaterial *mat = new TGeoMaterial("Agua",0,0,0); TGeoMedium *med = new TGeoMedium ("Agua",1,mat); TGeoMaterial *matglass = new TGeoMaterial("Vidro",0,0,0); TGeoMedium *medglass= new TGeoMedium ("Vidro",1,mat); // caixa do detector -------------------------------------------------- TGeoVolume *det = gGeoManager->MakeBox("Detector", med, DET_W1/2.0, DET_W2/2.0 + HEX_VERT, DET_H/2.0 + HEX_VERT); // espelho inferior --------------------------------------------------- TGeoShape *sbox = new TGeoBBox("mbox", DET_W1/2.0, DET_W2/2.0, HEX_VERT/2.0); TGeoShape *scone = new TGeoCone("cone", HEX_VERT/2.0*1.1, 0.0, PMT_RAIO, 0.0, HEX_LADO); TGeoTranslation *t1 = new TGeoTranslation("t1", -OFF_X, -OFF_Y1, 0); TGeoSubtraction *s1 = new TGeoSubtraction(sbox, scone, NULL, t1); TGeoCompositeShape *c1 = new TGeoCompositeShape("c1", s1); TGeoTranslation *t2 = new TGeoTranslation("t2", -OFF_X, -OFF_Y1+STEP_Y, 0); TGeoSubtraction *s2 = new TGeoSubtraction(c1, scone, NULL, t2); TGeoCompositeShape *c2 = new TGeoCompositeShape("c2", s2); TGeoTranslation *t3 = new TGeoTranslation("t3", -OFF_X, -OFF_Y1+2.0*STEP_Y, 0); TGeoSubtraction *s3 = new TGeoSubtraction(c2, scone, NULL, t3); TGeoCompositeShape *c3 = new TGeoCompositeShape("c3", s3); TGeoTranslation *t4 = new TGeoTranslation("t4", -OFF_X, -OFF_Y1+3.0*STEP_Y, 0); TGeoSubtraction *s4 = new TGeoSubtraction(c3, scone, NULL, t4); TGeoCompositeShape *c4 = new TGeoCompositeShape("c4", s4); TGeoTranslation *t5 = new TGeoTranslation("t5", -OFF_X+STEP_X, -OFF_Y2, 0); TGeoSubtraction *s5 = new TGeoSubtraction(c4, scone, NULL, t5); TGeoCompositeShape *c5 = new TGeoCompositeShape("c5", s5); TGeoTranslation *t6 = new TGeoTranslation("t6", -OFF_X+STEP_X, -OFF_Y2+STEP_Y, 0); TGeoSubtraction *s6 = new TGeoSubtraction(c5, scone, NULL, t6); TGeoCompositeShape *c6 = new TGeoCompositeShape("c6", s6); TGeoTranslation *t7 = new TGeoTranslation("t7", -OFF_X+STEP_X, -OFF_Y2+2.0*STEP_Y, 0); TGeoSubtraction *s7 = new TGeoSubtraction(c6, scone, NULL, t7); TGeoCompositeShape *c7 = new TGeoCompositeShape("c7", s7); TGeoTranslation *t8 = new TGeoTranslation("t8", -OFF_X+STEP_X, -OFF_Y2+3.0*STEP_Y, 0); TGeoSubtraction *s8 = new TGeoSubtraction(c7, scone, NULL, t8); TGeoCompositeShape *c8 = new TGeoCompositeShape("c8", s8); TGeoTranslation *t9 = new TGeoTranslation("t9", -OFF_X+2.0*STEP_X, -OFF_Y1, 0); TGeoSubtraction *s9 = new TGeoSubtraction(c8, scone, NULL, t9); TGeoCompositeShape *c9 = new TGeoCompositeShape("c9", s9); TGeoTranslation *t10 = new TGeoTranslation("t10", -OFF_X+2.0*STEP_X, -OFF_Y1+STEP_Y, 0); TGeoSubtraction *s10 = new TGeoSubtraction(c9, scone, NULL, t10); TGeoCompositeShape *c10 = new TGeoCompositeShape("c10", s10); TGeoTranslation *t11 = new TGeoTranslation("t11", -OFF_X+2.0*STEP_X, -OFF_Y1+2.0*STEP_Y, 0); TGeoSubtraction *s11 = new TGeoSubtraction(c10, scone, NULL, t11); TGeoCompositeShape *c11 = new TGeoCompositeShape("c11", s11); TGeoTranslation *t12 = new TGeoTranslation("t12", -OFF_X+2.0*STEP_X, -OFF_Y1+3.0*STEP_Y, 0); TGeoSubtraction *s12 = new TGeoSubtraction(c11, scone, NULL, t12); TGeoCompositeShape *c12 = new TGeoCompositeShape("c12", s12); TGeoTranslation *t13 = new TGeoTranslation("t13", -OFF_X+3.0*STEP_X, -OFF_Y2, 0); TGeoSubtraction *s13 = new TGeoSubtraction(c12, scone, NULL, t13); TGeoCompositeShape *c13 = new TGeoCompositeShape("c13", s13); TGeoTranslation *t14 = new TGeoTranslation("t14", -OFF_X+3.0*STEP_X, -OFF_Y2+STEP_Y, 0); TGeoSubtraction *s14 = new TGeoSubtraction(c13, scone, NULL, t14); TGeoCompositeShape *c14 = new TGeoCompositeShape("c14", s14); TGeoTranslation *t15 = new TGeoTranslation("t15", -OFF_X+3.0*STEP_X, -OFF_Y2+2.0*STEP_Y, 0); TGeoSubtraction *s15 = new TGeoSubtraction(c14, scone, NULL, t15); TGeoCompositeShape *c15 = new TGeoCompositeShape("c15", s15); TGeoTranslation *t16 = new TGeoTranslation("t16", -OFF_X+3.0*STEP_X, -OFF_Y2+3.0*STEP_Y, 0); TGeoSubtraction *s16 = new TGeoSubtraction(c15, scone, NULL, t16); TGeoCompositeShape *c16 = new TGeoCompositeShape("c16", s16); TGeoVolume *v16 = new TGeoVolume ("v16", c16, medglass); v16->SetLineColor(kGray); TGeoTranslation *t17 = new TGeoTranslation("t17", 0.0, 0.0, -(DET_H/2.0 + HEX_VERT/2.0)); TGeoRotation *r1 = new TGeoRotation("r1", 0.0, 180.0, 0.0); TGeoCombiTrans *ct1 = new TGeoCombiTrans("ct1", 0.0, 0.0, DET_H/2.0 + HEX_VERT/2.0, r1); det->AddNode(v16, 1, t17); det->AddNode(v16, 2, ct1); // espelho lateral ---------------------------------------------------- TGeoBBox *slbox = new TGeoBBox("lbox", DET_W1/2.0, DET_H/2.0 + HEX_VERT, HEX_VERT/2.0); TGeoTranslation *t18 = new TGeoTranslation("t18", HEX_ALTU, 0, 0); TGeoSubtraction *s18 = new TGeoSubtraction(slbox, scone, NULL, t18); TGeoCompositeShape *c18 = new TGeoCompositeShape("c18", s18); TGeoTranslation *t19 = new TGeoTranslation("t19", -HEX_ALTU, 0, 0); TGeoSubtraction *s19 = new TGeoSubtraction(c18, scone, NULL, t19); TGeoCompositeShape *c19 = new TGeoCompositeShape("c19", s19); TGeoTranslation *t20 = new TGeoTranslation("t20", 0, STEP_X, 0); TGeoSubtraction *s20 = new TGeoSubtraction(c19, scone, NULL, t20); TGeoCompositeShape *c20 = new TGeoCompositeShape("c20", s20); TGeoTranslation *t21 = new TGeoTranslation("t21", 0, -STEP_X, 0); TGeoSubtraction *s21 = new TGeoSubtraction(c20, scone, NULL, t21); TGeoCompositeShape *c21 = new TGeoCompositeShape("c21", s21); TGeoVolume *v21 = new TGeoVolume ("v21", c21, medglass); v21->SetLineColor(kGray); TGeoRotation *r22 = new TGeoRotation("r22", 0.0, 90.0, 0.0); TGeoCombiTrans *ct22 = new TGeoCombiTrans("ct22", 0.0, DET_W2/2.0 + HEX_VERT/2.0, 0.0, r22); det->AddNode(v21, 1, ct22); TGeoTranslation *t23 = new TGeoTranslation("t23", STEP_X, 0, 0); TGeoSubtraction *s23 = new TGeoSubtraction(slbox, scone, NULL, t23); TGeoCompositeShape *c23 = new TGeoCompositeShape("c23", s23); TGeoTranslation *t24 = new TGeoTranslation("t24", -STEP_X, 0, 0); TGeoSubtraction *s24 = new TGeoSubtraction(c23, scone, NULL, t24); TGeoCompositeShape *c24 = new TGeoCompositeShape("c24", s24); TGeoTranslation *t25 = new TGeoTranslation("t25", 0, HEX_ALTU, 0); TGeoSubtraction *s25 = new TGeoSubtraction(c24, scone, NULL, t25); TGeoCompositeShape *c25 = new TGeoCompositeShape("c25", s25); TGeoTranslation *t26 = new TGeoTranslation("t26", 0, -HEX_ALTU, 0); TGeoSubtraction *s26 = new TGeoSubtraction(c25, scone, NULL, t26); TGeoCompositeShape *c26 = new TGeoCompositeShape("c26", s26); TGeoVolume *v26 = new TGeoVolume ("v26", c26, medglass); v26->SetLineColor(kGray); TGeoRotation *r27 = new TGeoRotation("r27", 0.0, -90.0, 0.0); TGeoCombiTrans *ct27 = new TGeoCombiTrans("ct27", 0.0, -DET_W2/2.0 - HEX_VERT/2.0, 0.0, r27); det->AddNode(v26, 1, ct27); // fecha a geometria -------------------------------------------------- gGeoManager->SetTopVolume(det); gGeoManager->CloseGeometry(); gGeoManager->SetTopVisible(); det->SetTransparency(50); det->Draw("ogl"); // Inicializa uma trajetoria aleatoriamente --------------------------- Double_t *pos; pos = new Double_t[3]; Double_t *dir; dir = new Double_t[3]; TGeoNode *node = 0; // Reflete photon ate encontrar um vidro de PMT ----------------------- Int_t nphotons = 100; Double_t m, cosa; Double_t *n; Int_t N = -1; Double_t comp_tot = 0.0; TString start_path; for (Int_t i = 0; i < nphotons; i++) { node = 0; // Generate photon while (node != gGeoManager->GetTopNode()) { pos[0] = gRandom->Uniform(-0.5, 0.5); pos[1] = gRandom->Uniform(-0.5, 0.5); pos[2] = gRandom->Uniform(-0.5, 0.5); // Make sure you generate in the top volume node = gGeoManager->FindNode(pos[0],pos[1],pos[2]); } start_path = gGeoManager->GetPath(); // Double_t theta = TMath::ACos((1.-2.*gRandom->Rndm())); // Double_t phi = TMath::TwoPi()*gRandom->Rndm(); Double_t theta = theta0*TMath::DegToRad(); Double_t phi = phi0*TMath::DegToRad(); dir[0] = TMath::Sin(theta)*TMath::Cos(phi); dir[1] = TMath::Sin(theta)*TMath::Sin(phi); dir[2] = TMath::Cos(theta); // inicializa track ------------------------------------------- TPolyLine3D *pl = new TPolyLine3D(10); pl->SetLineColor(kRed); pl->SetPoint(0, pos[0], pos[1], pos[2]); comp_tot = 0; N = 0; while (1) { // propaga track ate proximo objeto --------------------------- gGeoManager->SetCurrentDirection(dir); printf("dir:%f,%f,%f\n", dir[0], dir[1], dir[2]); node = gGeoManager->FindNextBoundaryAndStep(); n = gGeoManager->FindNormalFast(); comp_tot += gGeoManager->GetStep(); const Double_t *crtpos = gGeoManager->GetCurrentPoint(); memcpy(pos, crtpos, 3*sizeof(Double_t)); pl->SetNextPoint(pos[0], pos[1], pos[2]); if (!node || node->GetVolume()->GetMedium() != medglass) break; // Reflection N++; TPolyLine3D *norm = new TPolyLine3D(2); norm->SetLineColor(kGreen); norm->SetPoint(0, pos[0], pos[1], pos[2]); norm->SetPoint(1, pos[0] - n[0]*50.0, pos[1] - n[1]*50.0, pos[2] - n[2]*50.0); norm->Draw(); // computa proxima direcao de propagacao ---------------------- cosa = dir[0]*n[0]+dir[1]*n[1]+dir[2]*n[2]; dir[0] = dir[0] - 2.0*cosa*n[0]; dir[1] = dir[1] - 2.0*cosa*n[1]; dir[2] = dir[2] - 2.0*cosa*n[2]; TMath::Normalize(dir); gGeoManager->cd(start_path); } pl->Draw(); cout << "Track " << i << ": Numero de reflexoes = " << N << endl; cout << "Comprimento total = " << comp_tot/1000.0 << " metros." << endl; } }