macros/leptoquarks_fun4all/G4_Svtx.C

0001 // reference SVTX macro used for MIE projections
0002
0003 int Min_si_layer = 0;
0004 int Max_si_layer = 6;
0005
0006 void SvtxInit(int verbosity = 0)
0007 {
0008   Min_si_layer = 0;
0009   Max_si_layer = 6;
0010 }
0011
0012 double Svtx(PHG4Reco* g4Reco, double radius,
0013         const int absorberactive = 0,
0014         int verbosity = 0) {
0015
0016   float svtx_inner_radius = 2.71;
0017
0018   if (radius > svtx_inner_radius) {
0019     cout << "inconsistency: radius: " << radius
0020      << " larger than SVTX inner radius: " << svtx_inner_radius << endl;
0021     gSystem->Exit(-1);
0022   }
0023
0024   //---------------
0025   // Load libraries
0026   //---------------
0027
0028   gSystem->Load("libg4detectors.so");
0029   gSystem->Load("libg4testbench.so");
0030
0031   PHG4CylinderSubsystem *cyl;
0032
0033   //======================================================================================================
0034   // The thicknesses from Yasuyuki on June 12, 2014 are as follows:
0035   // For Si 1mm = 1.07% X_0
0036   // For Cu 1mm = 6.96% X_0
0037   // The thickness of the tracking layers is:
0038   // Pixels:         1.3% X_0  (0.21% sensor +  1.07% support)  sensor = 200 mc Si, support = 154 mc Cu
0039   // Stripixel:      5%   X_0  (0.67% sensor + 4.3% support)    sensor = 624 mc Si, support = 618 mc Cu
0040   // Outer strips:   2%   X_0 (conservative)  (0.34% sensor + 1.66% support) sensor = 320 mc Si, support = 238 mc Cu
0041   //=======================================================================================================
0042
0043   double si_thickness[7] = {0.02, 0.02, 0.032, 0.032, 0.032, 0.032, 0.032};
0044   double svxrad[7] = {svtx_inner_radius, 4.63, 9.5, 10.5, 44.5, 45.5, 80.0}; // provides 98 MeV Upsilon resolution
0045   // Thicknesses (in % X_0) of 1.3,1.3,2.7/2,2.7/2,2.0/2,2.0/2,2.0 - YA's most conservative case
0046   double support_thickness[7] = {0.0154, 0.0154, 0.0338/2.0, 0.0338/2.0, 0.0238/2.0, 0.0238/2.0, 0.0238};
0047   double length[7] = {20., 20., -1, -1., - 1., - 1., -1}; // -1 use eta coverage to determine length
0048
0049   // here is our silicon:
0050   double inner_radius = radius;
0051   for (int ilayer = Min_si_layer; ilayer <= Max_si_layer; ilayer++)
0052     {
0053       cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0054       radius = svxrad[ilayer];
0055       // protect against installing layer with radius < inner radius from argument
0056       if (radius < inner_radius)
0057     {
0058       cout << "current radius " << radius << " smaller than inner radius "
0059            << inner_radius << endl;
0060       gSystem->Exit(-1);
0061     }
0062       cyl->set_double_param("radius",radius);
0063       if (length[ilayer] > 0)
0064         {
0065           cyl->set_int_param("lengthviarapidity",0);
0066           cyl->set_double_param("length",length[ilayer]);
0067         }
0068       else
0069     {
0070           cyl->set_int_param("lengthviarapidity",1);
0071     }
0072       cyl->set_string_param("material","G4_Si");
0073       cyl->set_double_param("thickness",si_thickness[ilayer]);
0074       cyl->SetActive();
0075       cyl->SuperDetector("SVTX");
0076       g4Reco->registerSubsystem( cyl );
0077
0078       radius += si_thickness[ilayer] + no_overlapp;
0079       cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", ilayer);
0080       cyl->set_double_param("radius",radius);
0081       if (length[ilayer] > 0)
0082         {
0083           cyl->set_int_param("lengthviarapidity",0);
0084           cyl->set_double_param("length",length[ilayer]);
0085         }
0086       else
0087     {
0088           cyl->set_int_param("lengthviarapidity",1);
0089     }
0090       cyl->set_string_param("material","G4_Cu");
0091       cyl->set_double_param("thickness",support_thickness[ilayer]);
0092       if (absorberactive)  cyl->SetActive();
0093       cyl->SuperDetector("SVTXSUPPORT");
0094       g4Reco->registerSubsystem( cyl );
0095     }
0096   if (ilayer != (Max_si_layer+1)) // coming out of the loop, layer is layer+1
0097     {
0098       cout << "layer number mismatch for Max_si_layer, Max_si_layer "
0099            << Max_si_layer << " should be " << ilayer << endl;
0100       gSystem->Exit(-1);
0101     }
0102   radius += support_thickness[Max_si_layer] + no_overlapp;
0103
0104   if (verbosity > 0) {
0105     cout << "============================ G4_Svtx.C::Svtx() ============================" << endl;
0106     cout << " SVTX Material Description:" << endl;
0107     for (int ilayer = Min_si_layer; ilayer <= Max_si_layer; ilayer++) {
0108       cout << "  layer " << ilayer
0109        << "  radius " << svxrad[ilayer]
0110        << "  zlength " << length[ilayer]
0111        << "  thickness (Si) " << si_thickness[ilayer]
0112        << "  support thickness (Cu) " << support_thickness[ilayer]
0113        << endl;
0114     }
0115     cout << "===========================================================================" << endl;
0116   }
0117   return radius;
0118 }
0119
0120 void Svtx_Cells(int verbosity = 0)
0121 {
0122   // runs the cellularization of the energy deposits (g4hits)
0123   // into detector hits (g4cells)
0124
0125   //---------------
0126   // Load libraries
0127   //---------------
0128
0129   gSystem->Load("libfun4all.so");
0130   gSystem->Load("libg4detectors.so");
0131
0132   //---------------
0133   // Fun4All server
0134   //---------------
0135
0136   Fun4AllServer *se = Fun4AllServer::instance();
0137
0138   //-----------
0139   // SVTX cells
0140   //-----------
0141
0142   // The pattern recognition layers (4 & 6) are set at 2 mm in Z and 240 microns
0143   // pitch to make the area the same as the long strips
0144   // Layers 3, 5 and 7 are long strips parallel to the beam line
0145
0146   // 60 micron X strips, 240 micron pattern reco strips
0147   double svxcellsizex[7] = {0.0050, 0.0050, 0.0060, 0.0240, 0.0060, 0.0240, 0.0060};
0148
0149   // 8 mm tracking strips, 2 mm pattern reco strips
0150   double svxcellsizey[7] = {0.0425, 0.0425, 0.8000, 0.2000, 0.8000, 0.2000, 0.8000};
0151
0152   PHG4CylinderCellReco *svtx_cells = new PHG4CylinderCellReco();
0153   svtx_cells->Detector("SVTX");
0154   svtx_cells->Verbosity(verbosity);
0155   int idx = 0;
0156   for (int i = Min_si_layer; i <= Max_si_layer; ++i) {
0157     svtx_cells->cellsize(i, svxcellsizex[idx], svxcellsizey[idx]);
0158     ++idx;
0159   }
0160   se->registerSubsystem(svtx_cells);
0161
0162   return;
0163 }
0164
0165 void Svtx_Reco(int verbosity = 0)
0166 {
0167   // reconstructs the MIE Svtx v2 detector (7 layers)
0168
0169   // requires Svtx setup and Svtx cell routines
0170   // prefers calorimeter reconstruction prior (once projections are working)
0171
0172   //---------------
0173   // Load libraries
0174   //---------------
0175
0176   gSystem->Load("libfun4all.so");
0177   gSystem->Load("libg4hough.so");
0178
0179   //---------------
0180   // Fun4All server
0181   //---------------
0182
0183   Fun4AllServer *se = Fun4AllServer::instance();
0184
0185   //----------------------------------
0186   // Digitize the cell energy into ADC
0187   //----------------------------------
0188   // defaults to 8-bit ADC with MIP at 0.25% dynamic range
0189   PHG4SvtxDigitizer* digi = new PHG4SvtxDigitizer();
0190   digi->Verbosity(verbosity);
0191   digi->set_adc_scale(0, 255, 1.0e-6); // 1.0 keV / bit
0192   digi->set_adc_scale(1, 255, 1.0e-6); // 1.0 keV / bit
0193   digi->set_adc_scale(2, 255, 1.6e-6); // 1.6 keV / bit
0194   digi->set_adc_scale(3, 255, 1.6e-6); // 1.6 keV / bit
0195   digi->set_adc_scale(4, 255, 1.6e-6); // 1.6 keV / bit
0196   digi->set_adc_scale(5, 255, 1.6e-6); // 1.6 keV / bit
0197   digi->set_adc_scale(6, 255, 1.6e-6); // 1.6 keV / bit
0198   se->registerSubsystem( digi );
0199
0200   //-------------------------------------
0201   // Apply Live Area Inefficiency to Hits
0202   //-------------------------------------
0203   // defaults to 1.0 (fully active)
0204   PHG4SvtxDeadArea* deadarea = new PHG4SvtxDeadArea();
0205   deadarea->Verbosity(verbosity);
0206   // deadarea->set_hit_efficiency(0,0.90);
0207   // deadarea->set_hit_efficiency(1,0.90);
0208   // deadarea->set_hit_efficiency(2,0.98);
0209   // deadarea->set_hit_efficiency(3,0.98);
0210   // deadarea->set_hit_efficiency(4,0.98);
0211   // deadarea->set_hit_efficiency(5,0.98);
0212   // deadarea->set_hit_efficiency(6,0.98);
0213   se->registerSubsystem( deadarea );
0214
0215   //-----------------------------
0216   // Apply MIP thresholds to Hits
0217   //-----------------------------
0218   PHG4SvtxThresholds* thresholds = new PHG4SvtxThresholds();
0219   thresholds->Verbosity(verbosity);
0220   thresholds->set_threshold(0,0.33);
0221   thresholds->set_threshold(1,0.33);
0222   thresholds->set_threshold(2,0.33);
0223   thresholds->set_threshold(3,0.33);
0224   thresholds->set_threshold(4,0.33);
0225   thresholds->set_threshold(5,0.33);
0226   thresholds->set_threshold(6,0.33);
0227   thresholds->set_use_thickness_mip(0, true);
0228   se->registerSubsystem( thresholds );
0229
0230   //-------------
0231   // Cluster Hits
0232   //-------------
0233   // needs to have clusters hold hit ids not cell ids
0234   // will require changes to evaluation
0235   PHG4SvtxClusterizer* clusterizer = new PHG4SvtxClusterizer();
0236   clusterizer->Verbosity(verbosity);
0237   clusterizer->set_threshold(0.33);
0238   clusterizer->set_z_clustering(2, false);
0239   clusterizer->set_z_clustering(3, false);
0240   clusterizer->set_z_clustering(4, false);
0241   clusterizer->set_z_clustering(5, false);
0242   clusterizer->set_z_clustering(6, false);
0243   // clusterizer->set_energy_weighting(2, true);
0244   // clusterizer->set_energy_weighting(3, true);
0245   // clusterizer->set_energy_weighting(4, true);
0246   // clusterizer->set_energy_weighting(5, true);
0247   // clusterizer->set_energy_weighting(6, true);
0248   se->registerSubsystem( clusterizer );
0249
0250   //---------------------
0251   // Track reconstruction
0252   //---------------------
0253   PHG4HoughTransform* hough = new PHG4HoughTransform(7,7);
0254   hough->Verbosity(verbosity);
0255   hough->set_mag_field(1.4);
0256   hough->set_material(0, 0.013);
0257   hough->set_material(1, 0.013);
0258   hough->set_material(2, 0.013);
0259   hough->set_material(3, 0.013);
0260   hough->set_material(4, 0.010);
0261   hough->set_material(5, 0.010);
0262   hough->set_material(6, 0.020);
0263   hough->setPtRescaleFactor(0.9972);
0264   hough->set_chi2_cut_init(3.0);
0265   hough->set_chi2_cut_full(3.0);
0266   hough->set_ca_chi2_cut(3.0);
0267   hough->setCutOnDCA(true);
0268   hough->setDCACut(0.1);
0269   hough->setDCAZCut(0.1);
0270   hough->setRejectGhosts(false);
0271   hough->setRemoveHits(false);
0272   se->registerSubsystem( hough );
0273
0274   //---------------------
0275   // Ghost rejection
0276   //---------------------
0277   // needs updates to merge split tracks when possible
0278   PHG4TrackGhostRejection* rejection = new PHG4TrackGhostRejection(7);
0279   rejection->Verbosity(verbosity);
0280   rejection->set_max_shared_hits(3);
0281   se->registerSubsystem( rejection );
0282
0283   //------------------------
0284   // Final Track Refitting
0285   //------------------------
0286   // PHG4TrackKalmanFitter *kalman = new PHG4TrackKalmanFitter
0287   // we need a module to redo the Kalman fit with G4 material and real mag field
0288   // to update the track container
0289
0290   //------------------------
0291   // Primary Track Refitting
0292   //------------------------
0293   // PHG4TrackKalmanFitter *kalman = new PHG4TrackKalmanFitter
0294   // we need a module to redo the Kalman fit including the vertex position
0295   // and create a separate stream of output tracks
0296
0297   //------------------
0298   // Track Projections
0299   //------------------
0300   PHG4SvtxTrackProjection* projection = new PHG4SvtxTrackProjection();
0301   projection->Verbosity(verbosity);
0302   se->registerSubsystem( projection );
0303
0304   //----------------------
0305   // Beam Spot Calculation
0306   //----------------------
0307   PHG4SvtxBeamSpotReco* beamspot = new PHG4SvtxBeamSpotReco();
0308   beamspot->Verbosity(verbosity);
0309   se->registerSubsystem( beamspot );
0310
0311   return;
0312 }
0313
0314 void G4_Svtx_Reco()
0315 {
0316   cout << "\033[31;1m"
0317        << "Warning: G4_Svtx_Reco() was moved to G4_Svtx.C and renamed to Svtx_Reco(), please update macros"
0318        << "\033[0m" << endl;
0319   Svtx_Reco();
0320
0321   return;
0322 }
0323
0324 void Svtx_Eval(std::string outputfile, int verbosity = 0)
0325 {
0326   //---------------
0327   // Load libraries
0328   //---------------
0329
0330   gSystem->Load("libfun4all.so");
0331   gSystem->Load("libg4detectors.so");
0332   gSystem->Load("libg4hough.so");
0333   gSystem->Load("libg4eval.so");
0334
0335   //---------------
0336   // Fun4All server
0337   //---------------
0338
0339   Fun4AllServer *se = Fun4AllServer::instance();
0340
0341   //----------------
0342   // SVTX evaluation
0343   //----------------
0344
0345   SubsysReco* eval = new SvtxEvaluator("SVTXEVALUATOR", outputfile.c_str());
0346   eval->Verbosity(verbosity);
0347   se->registerSubsystem( eval );
0348
0349   return;
0350 }