{ char *name = "test"; char *listing = "test"; char *upcom = "Electromagnetic shower, crystal with APD"; char *downcom = "Cascade inside crystal, displaced along Oz"; gROOT->ProcessLine(".x InitSLitrani.C(5,name,listing,upcom,downcom,kTRUE,kTRUE,kFALSE)"); // Colors const char *comTL = "Cascade is a gamma of 100 Gev"; const Color_t matVacuumColor = 0; const Color_t PbWO4Color = 9; const Color_t AluColor = 3; const Color_t GlueColor = 49; const Color_t TotAbsColor = 1; const Color_t SiliciumColor = 2; const Color_t CradleColor = 5; const Double_t IrrA = 0.0; //Value used for A when irrelevant ! const Double_t IrrZ = 0.0; //Value used for Z when irrelevant ! const Double_t IrrRho = 0.0; //Value used for Rho when irrelevant ! Int_t mediumindex = 1; gCleanOut->fMaxInfo = 1000; Short_t ns = -1; Bool_t ok = kFALSE; ns = 0; //____________________________________________________________________________ //____________________________________________________________________________ // // Building the geometry //____________________________________________________________________________ //____________________________________________________________________________ // TGeoManager *geom = new TGeoManager("setup","test of new Litrani"); //____________________________________________________________________________ // // Materials and media //____________________________________________________________________________ // // (1) Vacuum for "TOP" // TGeoMaterial *vacuum_mat = new TGeoMaterial("Vacuum",IrrA,IrrZ,IrrRho); TGeoMedium *vacuum = new TGeoMedium("Vacuum",mediumindex++,vacuum_mat); // // (2) Definition of PbWO4! // const Double_t PbWO4_radlen = 0.893; //radiation length of PbWO4 const Double_t PbWO4_intlen = 19.5; //interaction length of PbWO4 const Double_t PbWO4_rMoliere = 2.0; //Moliere radius of PbWO4 in cm TGeoElementTable *table = gGeoManager->GetElementTable(); TGeoElement *Pb = table->FindElement("LEAD"); TGeoElement *W = table->FindElement("TUNGSTEN"); TGeoElement *O = table->FindElement("OXYGEN"); TGeoMixture *pbwo4_mix = new TGeoMixture("PbWO4",3,8.28); pbwo4_mix->AddElement(Pb,1); pbwo4_mix->AddElement(W,1); pbwo4_mix->AddElement(O,4); pbwo4_mix->SetRadLen(-PbWO4_radlen,PbWO4_intlen); TLitMedium *pbwo4 = new TLitMedium("PbWO4",mediumindex++,pbwo4_mix,kFALSE,1.0,"AbsorptionLength_PbWO4"); const Double_t ns0 = 1.5861; //index 0 for the Sellmeier law for PbWO4 const Double_t ns1 = 1.1062; //index 1 for the Sellmeier law for PbWO4 const Double_t ls1 = 270.63; //1st resonant wavelength for the Sellmeier law for PbWO4 pbwo4->IsIsotropic(1,ns0,ns1,ls1); pbwo4->FindSpectrum("Spectrum_PbWO4La"); pbwo4->SetCascades(PbWO4_rMoliere); //only useful when TLitCascade used pbwo4->SetPhotMev(100.0); // // (3) Full pedantic definition of Silicium // TGeoElement *Si = table->FindElement("SILICON"); TGeoMaterial *silicium_mat = new TGeoMaterial("Silicium",Si,2.33); TLitMedium *silicium = new TLitMedium("Silicium",mediumindex++,silicium_mat,kTRUE,1.0, "AbsorptionLength_Silicium"); silicium->IsIsotropic("RefrIndex_Silicium"); // // (4) Definition of medium "glue" // TGeoMaterial *glue_mat = new TGeoMaterial("Glue",IrrA,IrrZ,IrrRho); switch (ns) { case 1: TLitMedium *glue = new TLitMedium("Glue",mediumindex++,glue_mat,kFALSE,1.0,10.0); break; default: TLitMedium *glue = new TLitMedium("Glue",mediumindex++,glue_mat,kFALSE,1.0, "AbsorptionLength_Meltmount1704"); break; } glue->IsIsotropic("RefrIndex_RTV3145"); //____________________________________________________________________________ // // Wrappings //____________________________________________________________________________ // // // (6) Definition of "TotAbsorbing": a totally absorbing wrapping // TGeoMaterial *totabs_mat = new TGeoMaterial("TotAbsorbing",IrrA,IrrZ,IrrRho); TLitMedium *totabsorbing = new TLitMedium("TotAbsorbing",mediumindex++,totabs_mat,-1.0); totabsorbing->SetAsWrapping(0.0,0.0,1.0,1.0,1.0); totabsorbing->IsIsotropic(1.0); // // (7) Definition of wrapping medium aluminium // const Double_t Air_RefrIndex = 1.0003; //Refractive index of air const Double_t aluminium_diffus = 0.01; const Double_t AluminiumSupplAbs = 0.1; TGeoMaterial *aluminium_mat = new TGeoMaterial("Aluminium",IrrA,IrrZ,IrrRho); TLitMedium *aluminium = new TLitMedium("Aluminium",mediumindex++,aluminium_mat,1.0); aluminium->SetAsWrapping(aluminium_diffus,"RIndexRev_Aluminium","IIndexRev_Aluminium",1.0,AluminiumSupplAbs); aluminium->IsIsotropic(Air_RefrIndex); // // (8) Definition of strontium // const Double_t rho_strontium = 2.64; const Double_t absl_strontium = 0.5; //Irrelevant. Chosen at random! const Double_t index_strontium = 2.5; //Irrelevant. Chosen at random! TGeoMaterial *sr_mat = new TGeoMaterial("Strontium",table->FindElement("STRONTIUM"),rho_strontium); TLitMedium *sr = new TLitMedium("Strontium",mediumindex++,sr_mat,kFALSE,1.0,absl_strontium); sr->IsIsotropic(index_strontium); //To avoid a warning in TLitMedium::NewWaveLength() gROOT->ProcessLine(".L PhotoEl_Strontium.C"); //Found in directory FitMacros, not in SplineFitDB.rdb TSplineFit *fitPESr = PhotoEl_Strontium(); sr->SetXSectnPE(fitPESr); //____________________________________________________________________________ // // Dimensions //____________________________________________________________________________ // // // Dimensions of crystal // const Double_t crys_dz = 10.0; //half length in z of crystal of PbWO4 const Double_t crys_AF = 8.0; //side left, small face |_ to OZ const Double_t crys_BF = 8.0; //side back, small face |_ to OZ const Double_t crys_CF = 8.0; //side right, small face |_ to OZ const Double_t crys_AR = 8.0; //side left, large face |_ to OZ const Double_t crys_BR = 8.0; //side back, large face |_ to OZ const Double_t crys_CR = 8.0; //-------------------------------Test------------------------------------ Double_t x0,x1,x2,x3,x4,x5,x6,x7; Double_t y0,y1,y2,y3,y4,y5,y6,y7; //small face |_ to OZ x0 = - crys_BF*0.5; y0 = - crys_AF*0.5; x1 = x0 + crys_BF; y1 = y0; x2 = x1; y2 = y1 + crys_CF; x3 = x0; y3 = y0 + crys_AF; //large face |_ to OZ x5 = x1; y5 = y1; x4 = x5 - crys_BR; y4 = y5; x6 = x5; y6 = y5 + crys_CR; x7 = x4; y7 = y4 + crys_AR; Double_t punkte[16]; punkte[0]=x0; punkte[1]=y0; punkte[6]=x1; punkte[7]=y1; punkte[4]=x2; punkte[5]=y2; punkte[2]=x3; punkte[3]=y3; punkte[8]=x4; punkte[9]=y4; punkte[14]=x5; punkte[15]=y5; punkte[12]=x6; punkte[13]=y6; punkte[10]=x7; punkte[11]=y7; //It is clockwise!!! Double_t *vertices=&punkte; // Vertices for Crystal const Double_t crys_dx = 4.0; // still necessary for Revetment const Double_t crys_dy = 4.0; // // Dimensions of APD // const Double_t apd_dx = 1.8; const Double_t apd_dy = 1.8; const Double_t apd_dz = 1.0; // // Dimensions for glue of APD // const Double_t glueapd_dx = apd_dx; const Double_t glueapd_dy = apd_dy; const Double_t glueapd_dz[2] = {0.05,0.5}; // // Width of wrapping is wid // const Double_t wid[2] = {0.1,0.5}; Double_t ws2 = 0.5*wid[ns]; // // Dimensions of TGEOBBox of alu around crystal "CRYSTALWS" // Double_t aluc_dx = crys_dx + wid[ns]; Double_t aluc_dy = crys_dy + wid[ns]; Double_t aluc_dz = crys_dz + wid[ns]; // // Start position of cascade and step size // // If inlocal is true, works in LOCAL COORDINATE SYSTEM, using 1st constructor of // TLitCascade. If false, works in WORLD COORDINATE SYSTEM using the 3rd constructor. // Both are valid and equivalent. Try both! // Bool_t inlocal = kFALSE; const Double_t epsz = 0.001; //To be just inside crystal const Double_t frspy = 0.4; //frspy == free space in y between crystals const Double_t frspz = 1.0; //frspz == free space in z below crystals const Double_t stepz = 1.0; //Starting point of cascade is displaced along Oz in steps stepz Double_t starty; if (inlocal) { starty = 0.0; //Start position of cascade in y in LCS of crystal startz = -crys_dz + epsz; //Start position of cascade in z in LCS of crystal const Double_t endz = 0.0; //End position of cascade in z in LCS of crystal } else { starty = -aluc_dy - 0.5*frspy; //Start position of cascade in y in WCS startz = -crys_dz + epsz; //Start position of cascade in z in WCS const Double_t endz = 0.0; //End position of cascade in z in WCS } Double_t posz = startz; Double_t arun = (endz-startz)/stepz; Int_t nrun = (Int_t)(arun+epsz); // // Dimensions of TGEOBBox of alu around APD and glue of APD [partially] "APDWS" // Double_t alua_dx = apd_dx + wid[ns]; Double_t alua_dy = apd_dy + wid[ns]; Double_t alua_dz = apd_dz + glueapd_dz[ns]; // // Dimensions of TotAbsorbing wrapping around APD and glue of APD // Double_t tot_dx = apd_dx + wid[ns]; Double_t tot_dy = apd_dy + wid[ns]; Double_t tot_dz = apd_dz + glueapd_dz[ns] + ws2; // // Dimensions of TOP // const Double_t top_dx = 10.0; const Double_t top_dy = 10.0; Double_t top_dz = wid[ns] + crys_dz + glueapd_dz[ns] + apd_dz + frspz + 2.0; //____________________________________________________________________________ // // Positionning (translations) //____________________________________________________________________________ // const Double_t t_alua_crys_z = aluc_dz + alua_dz; TGeoTranslation *t1 = new TGeoTranslation("t1",0.0,0.0,t_alua_crys_z); t1->RegisterYourself(); // t3 is for positionning the TGeoBBox of "totabsorbing" around APD and around glue of APD // with respect to the TGeoBBox of aluminium around crystal [or equivalently, with respect // to the composite shape "rev_shape"]. // const Double_t t_tot_crys_z = crys_dz + alua_dz + ws2; TGeoTranslation *t3 = new TGeoTranslation("t3",0.0,0.0,t_tot_crys_z); Double_t t_apd_tot_z = glueapd_dz[ns] - ws2; TGeoTranslation *t4 = new TGeoTranslation("t4",0.0,0.0,t_apd_tot_z); const Double_t t_glueapd_tot_z = -(apd_dz + ws2); TGeoTranslation *t5 = new TGeoTranslation("t5",0.0,0.0,t_glueapd_tot_z); Double_t t_rev = aluc_dy + 0.5*frspy; TGeoTranslation *t6 = new TGeoTranslation("t6",0.0,-t_rev,0.0); //____________________________________________________________________________ // // Volumes and nodes //____________________________________________________________________________ // // // top box containing everything and defining the WCS, world coordinate system // TGeoVolume *top = geom->MakeBox("TOP",vacuum,top_dx,top_dy,top_dz); geom->SetTopVolume(top); // // Composite shape aluminium // TGeoBBox *revcrys_shape = new TGeoBBox("CRYSTALWS",aluc_dx,aluc_dy,aluc_dz); TGeoBBox *revapd_shape = new TGeoBBox("APDWS",alua_dx,alua_dy,alua_dz); TGeoCompositeShape *rev_shape = new TGeoCompositeShape("REV","CRYSTALWS+APDWS:t1"); TGeoVolume *rev = new TGeoVolume("REV",rev_shape,aluminium); top->AddNode(rev,1,t6); // // TGeoBBox of "totabsorbing" around APD and glue of APD // TGeoVolume *tot = geom->MakeBox("TOT",totabsorbing,tot_dx,tot_dy,tot_dz); rev->AddNode(tot,1,t3); // // TGeoBBox APD of silicium // TGeoVolume *apd = geom->MakeBox("APD",silicium,apd_dx,apd_dy,apd_dz); tot->AddNode(apd,1,t4); // apd being a detector, booking of a TLitVolume is necessary. All characteristics of the // APD will then be given by a call to TLitVolume::SetAPD(). TLitVolume *lit_apd = new TLitVolume(apd); lit_apd->SetAPD("GainProfile_CMSAPD",3); // // TGeoBBox of glue of APD // TGeoVolume *glueapd = geom->MakeBox("GlueAPD",glue,glueapd_dx,glueapd_dy,glueapd_dz[ns]); tot->AddNode(glueapd,1,t5); //-------------------------------------------------------------------------------------------- //------------Here is the responsible part for the errors!!!!--------------------------------- //TGeoVolume *crystal = geom->MakeTrd1("CRYSTAL",pbwo4,4.0,2.0,4.0,10.0); //TGeoVolume *crystal = geom->MakeTrd2("CRYSTAL",pbwo4,3.0,4.0,3.0,4.0,10.0); TGeoVolume *crystal = geom->MakeArb8("CRYSTAL",pbwo4,10,vertices); //--- //vertices is a pointer to the array which has the vertices coordinates as content (defined in line 177). //-------------------------------------------------------------------------------------------- rev->AddNode(crystal,1); geom->CloseGeometry(); geom->CheckOverlaps(0.01); //____________________________________________________________________________ //____________________________________________________________________________ // // End of Building the geometry //____________________________________________________________________________ //____________________________________________________________________________ // //____________________________________________________________________________ // // Colors and drawing //____________________________________________________________________________ // top->SetVisibility(kFALSE); top->SetVisContainers(); top->SetLineColor(1); top->SetLineWidth(1); rev->SetLineColor(AluColor); rev->SetLineWidth(1); tot->SetVisibility(kFALSE); tot->SetLineColor(TotAbsColor); tot->SetLineWidth(1); apd->SetVisibility(kTRUE); apd->SetLineColor(SiliciumColor); apd->SetLineWidth(1); glueapd->SetVisibility(kTRUE); glueapd->SetLineColor(GlueColor); glueapd->SetLineWidth(1); crystal->SetVisibility(kTRUE); crystal->SetLineColor(PbWO4Color); crystal->SetLineWidth(1); gGeoManager->SetTopVisible(1); gGeoManager->SetVisLevel(4); TLit::Get()->BookCanvas(4); gTwoPad->SetStateOfPads(Pad1Only); gTwoPad->SetAllGrey(); gTwoPad->CommentTL(comTL); top->Draw(""); //____________________________________________________________________________ //____________________________________________________________________________ // // Generation of photons //____________________________________________________________________________ //____________________________________________________________________________ // Int_t krun; TLit::Get()->PrintPeriod(25000); TString path = "/TOP_1/REV_1/CRYSTAL_1"; TVector3 StartCasc(0.0,starty,posz); TVector3 StartCasc(0.0,0.0,0.0); TVector3 AxisCasc(0.0,0.0,1.0); TLitCascade *casc; if (inlocal) casc = new TLitCascade("gamma","gamma",path.Data(),"PbWO4",22,5,StartCasc,AxisCasc); else casc = new TLitCascade("gamma","gamma","PbWO4",StartCasc,AxisCasc,22,5); for (krun=1;krun<=1;krun++) { casc->Gen(krun,1,posz); posz += stepz; if (krun != nrun) { StartCasc.SetXYZ(0.0,starty,posz); if (inlocal) casc->NewPosition(path.Data(),StartCasc,AxisCasc); else casc->NewPosition(StartCasc,AxisCasc); } } gLitGp->SetTitle("position z of cascade"); gLitGp->Summary(); gTwoPad->SetStateOfPads(Pad1AndPad2); gTwoPad->ChangePad(); gLitGs->DoStat(); TLit::Get()->CloseFiles(); }