macros/leptoquarks_fun4all/G4Setup_fsPHENIX.C

0001
0002 double no_overlapp = 0.0001; // added to radii to avoid overlapping volumes
0003 bool overlapcheck = false; // set to true if you want to check for overlaps
0004
0005 void G4Init(bool do_svtx = true,
0006         bool do_cemc = true,
0007         bool do_hcalin = true,
0008         bool do_magnet = true,
0009         bool do_hcalout = true,
0010             bool do_pipe = true,
0011             bool do_FGEM = true,
0012         bool do_FEMC = true,
0013         bool do_FHCAL = true,
0014             int n_TPC_layers = 40) {
0015
0016   // load detector/material macros and execute Init() function
0017
0018   if (do_pipe)
0019     {
0020       gROOT->LoadMacro("G4_Pipe.C");
0021       PipeInit();
0022     }
0023   if (do_svtx)
0024     {
0025       gROOT->LoadMacro("G4_Svtx_maps_ladders+intt_ladders+tpc_KalmanPatRec.C");
0026       SvtxInit(n_TPC_layers);
0027     }
0028
0029   if (do_cemc)
0030     {
0031       gROOT->LoadMacro("G4_CEmc_Spacal.C");
0032       CEmcInit(72); // make it 2*2*2*3*3 so we can try other combinations
0033     }
0034
0035   if (do_hcalin)
0036     {
0037       gROOT->LoadMacro("G4_HcalIn_ref.C");
0038       HCalInnerInit();
0039     }
0040
0041   if (do_magnet)
0042     {
0043       gROOT->LoadMacro("G4_Magnet.C");
0044       MagnetInit();
0045     }
0046   if (do_hcalout)
0047     {
0048       gROOT->LoadMacro("G4_HcalOut_ref.C");
0049       HCalOuterInit();
0050     }
0051
0052   if (do_FGEM)
0053     {
0054       gROOT->LoadMacro("G4_FGEM_fsPHENIX.C");
0055       FGEM_Init();
0056     }
0057
0058   if (do_FEMC)
0059     {
0060       gROOT->LoadMacro("G4_FEMC.C");
0061       FEMCInit();
0062     }
0063
0064   if (do_FHCAL)
0065     {
0066       gROOT->LoadMacro("G4_FHCAL.C");
0067       FHCALInit();
0068     }
0069 }
0070
0071
0072 int G4Setup(const int absorberactive = 0,
0073         const string &field ="1.5",
0074         const EDecayType decayType = TPythia6Decayer::kAll,
0075         const bool do_svtx = true,
0076         const bool do_cemc = true,
0077         const bool do_hcalin = true,
0078         const bool do_magnet = true,
0079         const bool do_hcalout = true,
0080         const bool do_pipe = true,
0081         const bool do_FGEM = true,
0082         const bool do_FEMC = false,
0083         const bool do_FHCAL = false,
0084             const float magfield_rescale = 1.0) {
0085
0086   //---------------
0087   // Load libraries
0088   //---------------
0089
0090   gSystem->Load("libg4detectors.so");
0091   gSystem->Load("libg4testbench.so");
0092
0093   //---------------
0094   // Fun4All server
0095   //---------------
0096
0097   Fun4AllServer *se = Fun4AllServer::instance();
0098
0099   // read-in HepMC events to Geant4 if there is any
0100   HepMCNodeReader *hr = new HepMCNodeReader();
0101   se->registerSubsystem(hr);
0102
0103   PHG4Reco* g4Reco = new PHG4Reco();
0104   g4Reco->save_DST_geometry(true); //Save geometry from Geant4 to DST
0105   g4Reco->set_rapidity_coverage(1.1); // according to drawings
0106 // uncomment to set QGSP_BERT_HP physics list for productions
0107 // (default is QGSP_BERT for speed)
0108   //  g4Reco->SetPhysicsList("QGSP_BERT_HP");
0109
0110   if (decayType != TPythia6Decayer::kAll) {
0111     g4Reco->set_force_decay(decayType);
0112   }
0113
0114   double fieldstrength;
0115   istringstream stringline(field);
0116   stringline >> fieldstrength;
0117   if (stringline.fail()) { // conversion to double fails -> we have a string
0118
0119     if (field.find("sPHENIX.root") != string::npos) {
0120       g4Reco->set_field_map(field, 1);
0121     } else {
0122       g4Reco->set_field_map(field, 2);
0123     }
0124   } else {
0125     g4Reco->set_field(fieldstrength); // use const soleniodal field
0126   }
0127   g4Reco->set_field_rescale(magfield_rescale);
0128
0129   double radius = 0.;
0130
0131   //----------------------------------------
0132   // PIPE
0133   if (do_pipe) radius = Pipe(g4Reco, radius, absorberactive);
0134
0135   //----------------------------------------
0136   // SVTX
0137   if (do_svtx) radius = Svtx(g4Reco, radius, absorberactive);
0138
0139   //----------------------------------------
0140   // CEMC
0141 //
0142   if (do_cemc) radius = CEmc(g4Reco, radius, 8, absorberactive);
0143 //  if (do_cemc) radius = CEmc_Vis(g4Reco, radius, 8, absorberactive);// for visualization substructure of SPACAL, slow to render
0144
0145   //----------------------------------------
0146   // HCALIN
0147
0148   if (do_hcalin) radius = HCalInner(g4Reco, radius, 4, absorberactive);
0149
0150   //----------------------------------------
0151   // MAGNET
0152
0153   if (do_magnet) radius = Magnet(g4Reco, radius, 0, absorberactive);
0154
0155   //----------------------------------------
0156   // HCALOUT
0157
0158   if (do_hcalout) radius = HCalOuter(g4Reco, radius, 4, absorberactive);
0159
0160   //----------------------------------------
0161   // Forward tracking
0162
0163   if ( do_FGEM )
0164     FGEMSetup(g4Reco);
0165
0166   //----------------------------------------
0167   // FEMC
0168
0169   if ( do_FEMC )
0170     FEMCSetup(g4Reco, absorberactive);
0171
0172   //----------------------------------------
0173   // FHCAL
0174
0175   if ( do_FHCAL )
0176     FHCALSetup(g4Reco, absorberactive);
0177
0178   // sPHENIX forward flux return(s)
0179   PHG4CylinderSubsystem *flux_return_plus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
0180   flux_return_plus->set_int_param("lengthviarapidity",0);
0181   flux_return_plus->set_double_param("length",10.2);
0182   flux_return_plus->set_double_param("radius",2.1);
0183   flux_return_plus->set_double_param("thickness",263.5-5.0);
0184   flux_return_plus->set_double_param("place_z",335.9);
0185   flux_return_plus->set_string_param("material","G4_Fe");
0186   flux_return_plus->SetActive(false);
0187   flux_return_plus->SuperDetector("FLUXRET_ETA_PLUS");
0188   flux_return_plus->OverlapCheck(overlapcheck);
0189   g4Reco->registerSubsystem(flux_return_plus);
0190
0191   PHG4CylinderSubsystem *flux_return_minus = new PHG4CylinderSubsystem("FWDFLUXRET", 0);
0192   flux_return_minus->set_int_param("lengthviarapidity",0);
0193   flux_return_minus->set_double_param("length",10.2);
0194   flux_return_minus->set_double_param("radius",2.1);
0195   flux_return_minus->set_double_param("place_z",-335.9);
0196   flux_return_minus->set_double_param("thickness",263.5-5.0);
0197   flux_return_minus->set_string_param("material","G4_Fe");
0198   flux_return_minus->SetActive(false);
0199   flux_return_minus->SuperDetector("FLUXRET_ETA_MINUS");
0200   flux_return_minus->OverlapCheck(overlapcheck);
0201   g4Reco->registerSubsystem(flux_return_minus);
0202
0203   //----------------------------------------
0204   // piston magnet
0205   make_piston("magpiston", g4Reco);
0206
0207   //----------------------------------------
0208   // BLACKHOLE
0209
0210   // swallow all particles coming out of the backend of sPHENIX
0211   PHG4CylinderSubsystem *blackhole = new PHG4CylinderSubsystem("BH", 1);
0212 blackhole->set_double_param("radius",radius + 10); // add 10 cm
0213
0214   blackhole->set_int_param("lengthviarapidity",0);
0215   blackhole->set_double_param("length",g4Reco->GetWorldSizeZ() - no_overlapp); // make it cover the world in length
0216   blackhole->BlackHole();
0217   blackhole->set_double_param("thickness",0.1); // it needs some thickness
0218   blackhole->SetActive(); // always see what leaks out
0219   blackhole->OverlapCheck(overlapcheck);
0220   g4Reco->registerSubsystem(blackhole);
0221
0222   //----------------------------------------
0223   // FORWARD BLACKHOLEs
0224   // +Z
0225   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_PLUS", 1);
0226   blackhole->SuperDetector("BH_FORWARD_PLUS");
0227   blackhole->set_double_param("radius",0); // add 10 cm
0228   blackhole->set_int_param("lengthviarapidity",0);
0229   blackhole->set_double_param("length",0.1); // make it cover the world in length
0230   blackhole->set_double_param("place_z",g4Reco->GetWorldSizeZ()/2. - 0.1  - no_overlapp);
0231   blackhole->BlackHole();
0232   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0233   blackhole->SetActive(); // always see what leaks out
0234   blackhole->OverlapCheck(overlapcheck);
0235   g4Reco->registerSubsystem(blackhole);
0236
0237   blackhole = new PHG4CylinderSubsystem("BH_FORWARD_NEG", 1);
0238   blackhole->SuperDetector("BH_FORWARD_NEG");
0239   blackhole->set_double_param("radius",0); // add 10 cm
0240   blackhole->set_int_param("lengthviarapidity",0);
0241   blackhole->set_double_param("length",0.1); // make it cover the world in length
0242   blackhole->set_double_param("place_z", - g4Reco->GetWorldSizeZ()/2. +0.1  + no_overlapp);
0243   blackhole->BlackHole();
0244   blackhole->set_double_param("thickness",radius - no_overlapp); // it needs some thickness
0245   blackhole->SetActive(); // always see what leaks out
0246   blackhole->OverlapCheck(overlapcheck);
0247   g4Reco->registerSubsystem(blackhole);
0248
0249   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0250   g4Reco->registerSubsystem(truth);
0251   se->registerSubsystem( g4Reco );
0252 }
0253
0254
0255 void ShowerCompress(int verbosity = 0) {
0256
0257   gSystem->Load("libfun4all.so");
0258   gSystem->Load("libg4eval.so");
0259
0260   Fun4AllServer *se = Fun4AllServer::instance();
0261
0262   PHG4DstCompressReco* compress = new PHG4DstCompressReco("PHG4DstCompressReco");
0263   compress->AddHitContainer("G4HIT_PIPE");
0264   compress->AddHitContainer("G4HIT_SVTXSUPPORT");
0265   compress->AddHitContainer("G4HIT_CEMC_ELECTRONICS");
0266   compress->AddHitContainer("G4HIT_CEMC");
0267   compress->AddHitContainer("G4HIT_ABSORBER_CEMC");
0268   compress->AddHitContainer("G4HIT_CEMC_SPT");
0269   compress->AddHitContainer("G4HIT_ABSORBER_HCALIN");
0270   compress->AddHitContainer("G4HIT_HCALIN");
0271   compress->AddHitContainer("G4HIT_HCALIN_SPT");
0272   compress->AddHitContainer("G4HIT_MAGNET");
0273   compress->AddHitContainer("G4HIT_ABSORBER_HCALOUT");
0274   compress->AddHitContainer("G4HIT_HCALOUT");
0275   compress->AddHitContainer("G4HIT_BH_1");
0276   compress->AddHitContainer("G4HIT_BH_FORWARD_PLUS");
0277   compress->AddHitContainer("G4HIT_BH_FORWARD_NEG");
0278   compress->AddCellContainer("G4CELL_CEMC");
0279   compress->AddCellContainer("G4CELL_HCALIN");
0280   compress->AddCellContainer("G4CELL_HCALOUT");
0281   compress->AddTowerContainer("TOWER_SIM_CEMC");
0282   compress->AddTowerContainer("TOWER_RAW_CEMC");
0283   compress->AddTowerContainer("TOWER_CALIB_CEMC");
0284   compress->AddTowerContainer("TOWER_SIM_HCALIN");
0285   compress->AddTowerContainer("TOWER_RAW_HCALIN");
0286   compress->AddTowerContainer("TOWER_CALIB_HCALIN");
0287   compress->AddTowerContainer("TOWER_SIM_HCALOUT");
0288   compress->AddTowerContainer("TOWER_RAW_HCALOUT");
0289   compress->AddTowerContainer("TOWER_CALIB_HCALOUT");
0290
0291   compress->AddHitContainer("G4HIT_FEMC");
0292   compress->AddHitContainer("G4HIT_ABSORBER_FEMC");
0293   compress->AddHitContainer("G4HIT_FHCAL");
0294   compress->AddHitContainer("G4HIT_ABSORBER_FHCAL");
0295   compress->AddCellContainer("G4CELL_FEMC");
0296   compress->AddCellContainer("G4CELL_FHCAL");
0297   compress->AddTowerContainer("TOWER_SIM_FEMC");
0298   compress->AddTowerContainer("TOWER_RAW_FEMC");
0299   compress->AddTowerContainer("TOWER_CALIB_FEMC");
0300   compress->AddTowerContainer("TOWER_SIM_FHCAL");
0301   compress->AddTowerContainer("TOWER_RAW_FHCAL");
0302   compress->AddTowerContainer("TOWER_CALIB_FHCAL");
0303
0304   se->registerSubsystem(compress);
0305
0306   return;
0307 }
0308
0309 void DstCompress(Fun4AllDstOutputManager* out) {
0310   if (out) {
0311     out->StripNode("G4HIT_PIPE");
0312     out->StripNode("G4HIT_SVTXSUPPORT");
0313     out->StripNode("G4HIT_CEMC_ELECTRONICS");
0314     out->StripNode("G4HIT_CEMC");
0315     out->StripNode("G4HIT_ABSORBER_CEMC");
0316     out->StripNode("G4HIT_CEMC_SPT");
0317     out->StripNode("G4HIT_ABSORBER_HCALIN");
0318     out->StripNode("G4HIT_HCALIN");
0319     out->StripNode("G4HIT_HCALIN_SPT");
0320     out->StripNode("G4HIT_MAGNET");
0321     out->StripNode("G4HIT_ABSORBER_HCALOUT");
0322     out->StripNode("G4HIT_HCALOUT");
0323     out->StripNode("G4HIT_BH_1");
0324     out->StripNode("G4HIT_BH_FORWARD_PLUS");
0325     out->StripNode("G4HIT_BH_FORWARD_NEG");
0326     out->StripNode("G4CELL_CEMC");
0327     out->StripNode("G4CELL_HCALIN");
0328     out->StripNode("G4CELL_HCALOUT");
0329
0330     out->StripNode("G4HIT_FEMC");
0331     out->StripNode("G4HIT_ABSORBER_FEMC");
0332     out->StripNode("G4HIT_FHCAL");
0333     out->StripNode("G4HIT_ABSORBER_FHCAL");
0334     out->StripNode("G4CELL_FEMC");
0335     out->StripNode("G4CELL_FHCAL");
0336   }
0337 }
0338
0339 int
0340 make_piston(string name, PHG4Reco* g4Reco)
0341 {
0342   double be_pipe_radius    = 2.16;   // 2.16 cm based on spec sheet
0343   double be_pipe_thickness = 0.0760; // 760 um based on spec sheet
0344   double be_pipe_length    = 80.0;   // +/- 40 cm
0345
0346   double al_pipe_radius    = 2.16;   // same as Be pipe
0347   double al_pipe_thickness = 0.1600; // 1.6 mm based on spec
0348   double al_pipe_length    = 88.3;   // extension beyond +/- 40 cm
0349
0350   const double zpos0 = al_pipe_length + be_pipe_length * 0.5; // first large GEM station
0351   const double zpos1 = 305 - 20; // front of forward ECal/MPC
0352   const double zpos2 = 335.9 - 10.2 / 2.; // front of the forward field endcap
0353   const double calorimeter_hole_diamater = 9.92331 *2; // side length of the middle hole of MPC that can hold the piston. Also the max diameter of the piston in that region
0354
0355   const double beampipe_radius = be_pipe_radius;
0356
0357   // teeth cone section specific
0358   const double number_of_wteeth = 100;
0359   const double teeth_thickness = 0.3504 * 2; //2 X0
0360   const double eta_inner = -log(tan(atan((beampipe_radius + 0.1) / zpos0) / 2));
0361   const double eta_outter = 4.2;
0362   const double eta_teeth_outter = 4.05;
0363   double pos = zpos0 + (zpos1 - zpos0) / 2;
0364 //  cout << "MAGNETIC PISTON:" << eta_inner << " " << eta_outter << " " << pos
0365 //      << endl;
0366
0367   PHG4ConeSubsystem *magpiston = new PHG4ConeSubsystem("Piston", 0);
0368   magpiston->SetZlength((zpos1 - zpos0) / 2);
0369   magpiston->SetPlaceZ((zpos1 + zpos0) / 2);
0370   magpiston->SetR1(beampipe_radius,
0371       tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(eta_outter)) * zpos0);
0372   magpiston->SetR2(beampipe_radius,
0373       tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(eta_outter)) * zpos1);
0374   magpiston->SetMaterial("G4_Fe");
0375   magpiston->OverlapCheck(overlapcheck);
0376   g4Reco->registerSubsystem(magpiston);
0377
0378 //  PHG4ConeSubsystem *magpiston = new PHG4ConeSubsystem(name.c_str(), 1);
0379 //  magpiston->SetZlength((zpos1 - zpos0) / 2);
0380 //  magpiston->SetPlaceZ(pos);
0381 //  magpiston->Set_eta_range(eta_outter, eta_inner);
0382 //  magpiston->SetMaterial("G4_Fe");
0383 //  magpiston->SuperDetector(name.c_str());
0384 //  magpiston->SetActive(false);
0385 //  g4Reco->registerSubsystem(magpiston);
0386
0387   pos = zpos0 + 1.0 + teeth_thickness / 2;
0388   for (int i = 0; i < number_of_wteeth; i++)
0389     {
0390       stringstream s;
0391       s << name;
0392       s << "_teeth_";
0393       s << i;
0394
0395       magpiston = new PHG4ConeSubsystem(s.str(), i);
0396       magpiston->SuperDetector(name);
0397       magpiston->SetZlength(teeth_thickness / 2);
0398       magpiston->SetPlaceZ(pos);
0399       magpiston->SetR1(
0400           //
0401           tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(eta_outter - .01))
0402               * (pos - teeth_thickness / 2), //
0403           tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(eta_teeth_outter))
0404               * (pos - teeth_thickness / 2) //
0405               );
0406       magpiston->SetR2(
0407           //
0408           tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(eta_outter - .01))
0409               * (pos + teeth_thickness / 2), //
0410           tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(eta_outter - .01))
0411               * (pos + teeth_thickness / 2) + .1 //
0412               );
0413       magpiston->SetMaterial("G4_W");
0414       magpiston->SuperDetector(name.c_str());
0415       magpiston->SetActive(false);
0416       magpiston->OverlapCheck(overlapcheck);
0417       g4Reco->registerSubsystem(magpiston);
0418       pos += ((zpos1 - zpos0 - 10) / number_of_wteeth);
0419     }
0420
0421   // last piece connect to the field return
0422   PHG4CylinderSubsystem *magpiston2 = new PHG4CylinderSubsystem("Piston_EndSection", 0);
0423   magpiston2->set_int_param("lengthviarapidity",0);
0424   magpiston2->set_double_param("length",zpos2 - zpos1);
0425   magpiston2->set_double_param("place_z", (zpos2 + zpos1) / 2.);
0426   magpiston2->set_double_param("radius",beampipe_radius);
0427   magpiston2->set_double_param("thickness",calorimeter_hole_diamater / 2. - beampipe_radius);
0428   magpiston2->set_string_param("material","G4_Fe");
0429   magpiston2->SuperDetector(name);
0430   magpiston2->SetActive(false);
0431   magpiston2->OverlapCheck(overlapcheck);
0432   g4Reco->registerSubsystem(magpiston2);
0433
0434   return 0;
0435 }
0436