macros/leptoquarks_fun4all/G4_FGEM_EIC.C

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