RICHAnalysis_macros/macros_fun4all/G4_FGEM_EIC.C

0001 // $Id: G4_FGEM_fsPHENIX.C,v 1.2 2014/01/22 01:44:13 jinhuang Exp $
0002
0003 /*!
0004  * \file G4_FGEM_fsPHENIX.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision: 1.2 $
0008  * \date $Date: 2014/01/22 01:44:13 $
0009  */
0010
0011 using namespace std;
0012
0013 void
0014 FGEM_Init()
0015 {
0016
0017 }
0018
0019 void
0020 FGEMSetup(PHG4Reco* g4Reco, const int N_Sector = 8, //
0021     const double min_eta = 1.45 //
0022     )
0023 {
0024
0025   const double tilt = .1;
0026
0027   string name;
0028   double etamax;
0029   double etamin;
0030   double zpos;
0031   PHG4SectorSubsystem *gem;
0032
0033   make_GEM_station("FGEM_0", g4Reco, 17, 1.01, 2.7, N_Sector);
0034   make_GEM_station("FGEM_1", g4Reco, 62, 2.15, 4.0, N_Sector);
0035
0036   ///////////////////////////////////////////////////////////////////////////
0037
0038   name = "FGEM_2";
0039   etamax = 4;
0040   etamin = min_eta;
0041   zpos = 1.2e2;
0042
0043   gem = new PHG4SectorSubsystem(name.c_str());
0044
0045   gem->get_geometry().set_normal_polar_angle(tilt);
0046   gem->get_geometry().set_normal_start(
0047       zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
0048   gem->get_geometry().set_min_polar_angle(
0049       PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
0050   gem->get_geometry().set_max_polar_angle(
0051       PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamin));
0052   gem->get_geometry().set_max_polar_edge(
0053       PHG4Sector::Sector_Geometry::FlatEdge());
0054   gem->get_geometry().set_material("G4_METHANE");
0055   gem->get_geometry().set_N_Sector(N_Sector);
0056   gem->OverlapCheck(overlapcheck);
0057   AddLayers_MiniTPCDrift(gem);
0058   gem->get_geometry().AddLayers_HBD_GEM();
0059   g4Reco->registerSubsystem(gem);
0060
0061   ///////////////////////////////////////////////////////////////////////////
0062
0063   name = "FGEM_3";
0064   etamax = 4;
0065   etamin = min_eta;
0066   zpos = 1.6e2;
0067   gem = new PHG4SectorSubsystem(name.c_str());
0068
0069   gem->SuperDetector(name);
0070   gem->get_geometry().set_normal_polar_angle(tilt);
0071   gem->get_geometry().set_normal_start(
0072       zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
0073   gem->get_geometry().set_min_polar_angle(
0074       PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
0075   gem->get_geometry().set_max_polar_angle(
0076       PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
0077   gem->get_geometry().set_max_polar_edge(
0078       PHG4Sector::Sector_Geometry::FlatEdge());
0079   gem->get_geometry().set_material("G4_METHANE");
0080   gem->get_geometry().set_N_Sector(N_Sector);
0081   gem->OverlapCheck(overlapcheck);
0082   AddLayers_MiniTPCDrift(gem);
0083   gem->get_geometry().AddLayers_HBD_GEM();
0084   g4Reco->registerSubsystem(gem);
0085
0086   gem = new PHG4SectorSubsystem(name + "_LowerEta");
0087   gem->SuperDetector(name);
0088
0089   zpos = zpos
0090       - (zpos * sin(tilt)
0091           + zpos * cos(tilt)
0092               * tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2) - tilt))
0093           * sin(tilt);
0094
0095   gem->get_geometry().set_normal_polar_angle(
0096       (PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta)
0097           + PHG4Sector::Sector_Geometry::eta_to_polar_angle(2)) / 2);
0098   gem->get_geometry().set_normal_start(
0099       zpos * PHG4Sector::Sector_Geometry::Unit_cm(),
0100       PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
0101   gem->get_geometry().set_min_polar_angle(
0102       PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
0103   gem->get_geometry().set_max_polar_angle(
0104       PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta));
0105   gem->get_geometry().set_material("G4_METHANE");
0106   gem->get_geometry().set_N_Sector(N_Sector);
0107   gem->get_geometry().set_min_polar_edge(
0108       PHG4Sector::Sector_Geometry::FlatEdge());
0109
0110   AddLayers_MiniTPCDrift(gem);
0111   gem->get_geometry().AddLayers_HBD_GEM();
0112   gem->OverlapCheck(overlapcheck);
0113   g4Reco->registerSubsystem(gem);
0114
0115   ///////////////////////////////////////////////////////////////////////////
0116
0117   name = "FGEM_4";
0118   etamax = 4;
0119   etamin = min_eta;
0120   zpos = 2.75e2;
0121   gem = new PHG4SectorSubsystem(name.c_str());
0122
0123   gem->SuperDetector(name);
0124   gem->get_geometry().set_normal_polar_angle(tilt);
0125   gem->get_geometry().set_normal_start(
0126       zpos * PHG4Sector::Sector_Geometry::Unit_cm(), 0);
0127   gem->get_geometry().set_min_polar_angle(
0128       PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
0129   gem->get_geometry().set_max_polar_angle(
0130       PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
0131   gem->get_geometry().set_max_polar_edge(
0132       PHG4Sector::Sector_Geometry::FlatEdge());
0133   gem->get_geometry().set_material("G4_METHANE");
0134   gem->get_geometry().set_N_Sector(N_Sector);
0135   gem->OverlapCheck(overlapcheck);
0136   AddLayers_MiniTPCDrift(gem);
0137   gem->get_geometry().AddLayers_HBD_GEM();
0138   g4Reco->registerSubsystem(gem);
0139
0140   zpos = zpos
0141       - (zpos * sin(tilt)
0142           + zpos * cos(tilt)
0143               * tan(PHG4Sector::Sector_Geometry::eta_to_polar_angle(2) - tilt))
0144           * sin(tilt);
0145
0146   gem = new PHG4SectorSubsystem(name + "_LowerEta");
0147   gem->SuperDetector(name);
0148
0149   gem->get_geometry().set_normal_polar_angle(
0150       (PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta)
0151           + PHG4Sector::Sector_Geometry::eta_to_polar_angle(2)) / 2);
0152   gem->get_geometry().set_normal_start(
0153       zpos * PHG4Sector::Sector_Geometry::Unit_cm(),
0154       PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
0155   gem->get_geometry().set_min_polar_angle(
0156       PHG4Sector::Sector_Geometry::eta_to_polar_angle(2));
0157   gem->get_geometry().set_max_polar_angle(
0158       PHG4Sector::Sector_Geometry::eta_to_polar_angle(min_eta));
0159   gem->get_geometry().set_material("G4_METHANE");
0160   gem->get_geometry().set_N_Sector(N_Sector);
0161   gem->get_geometry().set_min_polar_edge(
0162       PHG4Sector::Sector_Geometry::FlatEdge());
0163
0164   AddLayers_MiniTPCDrift(gem);
0165   gem->get_geometry().AddLayers_HBD_GEM();
0166   gem->OverlapCheck(overlapcheck);
0167   g4Reco->registerSubsystem(gem);
0168
0169   ///////////////////////////////////////////////////////////////////////////
0170
0171 }
0172
0173 //! Add drift layers to mini TPC
0174 void
0175 AddLayers_MiniTPCDrift(PHG4SectorSubsystem *gem)
0176 {
0177   assert(gem);
0178
0179   const double cm = PHG4Sector::Sector_Geometry::Unit_cm();
0180   const double mm = .1 * cm;
0181   const double um = 1e-3 * mm;
0182
0183 //  const int N_Layers = 70; // used for mini-drift TPC timing digitalization
0184   const int N_Layers = 1; // simplified setup
0185   const double thickness = 2 * cm;
0186
0187   gem->get_geometry().AddLayer("EntranceWindow", "G4_MYLAR", 25 * um, false,
0188       100);
0189   gem->get_geometry().AddLayer("Cathode", "G4_GRAPHITE", 10 * um, false, 100);
0190
0191   for (int d = 1; d <= N_Layers; d++)
0192     {
0193       stringstream s;
0194       s << "DriftLayer_";
0195       s << d;
0196
0197       gem->get_geometry().AddLayer(s.str(), "G4_METHANE", thickness / N_Layers,
0198           true);
0199
0200     }
0201 }
0202
0203 int
0204 make_GEM_station(string name, PHG4Reco* g4Reco, double zpos, double etamin,
0205     double etamax,  const int N_Sector = 8)
0206 {
0207
0208 //  cout
0209 //      << "make_GEM_station - GEM construction with PHG4SectorSubsystem - make_GEM_station_EdgeReadout  of "
0210 //      << name << endl;
0211
0212   double polar_angle = 0;
0213
0214   if (zpos < 0)
0215     {
0216       zpos = -zpos;
0217       polar_angle = TMath::Pi();
0218
0219     }
0220   if (etamax < etamin)
0221     {
0222       double t = etamax;
0223       etamax = etamin;
0224       etamin = t;
0225     }
0226
0227   PHG4SectorSubsystem *gem;
0228   gem = new PHG4SectorSubsystem(name.c_str());
0229
0230   gem->SuperDetector(name);
0231
0232   gem->get_geometry().set_normal_polar_angle(polar_angle);
0233   gem->get_geometry().set_normal_start(
0234       zpos * PHG4Sector::Sector_Geometry::Unit_cm());
0235   gem->get_geometry().set_min_polar_angle(
0236       PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamax));
0237   gem->get_geometry().set_max_polar_angle(
0238       PHG4Sector::Sector_Geometry::eta_to_polar_angle(etamin));
0239   gem->get_geometry().set_max_polar_edge(
0240       PHG4Sector::Sector_Geometry::FlatEdge());
0241   gem->get_geometry().set_min_polar_edge(
0242       PHG4Sector::Sector_Geometry::FlatEdge());
0243   gem->get_geometry().set_N_Sector(N_Sector);
0244   gem->get_geometry().set_material("G4_METHANE");
0245   gem->OverlapCheck(overlapcheck);
0246
0247   AddLayers_MiniTPCDrift(gem);
0248   gem->get_geometry().AddLayers_HBD_GEM();
0249   g4Reco->registerSubsystem(gem);
0250
0251 }
0252
0253 void FGEM_FastSim_Reco(int verbosity = 0) {
0254
0255   //---------------
0256   // Load libraries
0257   //---------------
0258
0259   gSystem->Load("libfun4all.so");
0260   gSystem->Load("libg4hough.so");
0261
0262   //---------------
0263   // Fun4All server
0264   //---------------
0265
0266   Fun4AllServer *se = Fun4AllServer::instance();
0267
0268   PHG4TrackFastSim* kalman = new PHG4TrackFastSim("PHG4TrackFastSim");
0269   kalman->Verbosity(0);
0270
0271   kalman->set_use_vertex_in_fitting(true);
0272   kalman->set_vertex_xy_resolution(50E-4);
0273   kalman->set_vertex_z_resolution(50E-4);
0274
0275   kalman->set_detector_type(PHG4TrackFastSim::Vertical_Plane); // Vertical_Plane, Cylinder
0276   kalman->set_phi_resolution(50E-4);
0277   kalman->set_r_resolution(1.);
0278
0279   //kalman->set_mag_field_file_name("/phenix/upgrades/decadal/fieldmaps/fsPHENIX.2d.root");
0280   //kalman->set_mag_field_re_scaling_factor(1.);
0281
0282   kalman->set_pat_rec_hit_finding_eff(1.);
0283   kalman->set_pat_rec_noise_prob(0.);
0284
0285   std::string phg4hits_names[] = {"G4HIT_FGEM_0","G4HIT_FGEM_1","G4HIT_FGEM_2","G4HIT_FGEM_3","G4HIT_FGEM_4"};
0286   kalman->set_phg4hits_names(phg4hits_names, 5);
0287   kalman->set_sub_top_node_name("SVTX");
0288   kalman->set_trackmap_out_name("SvtxTrackMap");
0289
0290   // Saved track states (projections)
0291   std::string state_names[] = {"RICH"};
0292   kalman->set_state_names(state_names, 1);
0293
0294   kalman->set_fit_alg_name("KalmanFitterRefTrack");//
0295   kalman->set_primary_assumption_pid(13);
0296   kalman->set_do_evt_display(false);
0297
0298   se->registerSubsystem(kalman);
0299
0300 }
0301