macros/leptoquarks_fun4all/G4_Svtx_MAPScyl_ITTcyl_TPC.C

0001 #include "G4_TPC.C"
0002
0003 const int n_ib_layer = 3;   // number of maps inner barrel layers
0004 const int n_intt_layer = 4; // number of int. tracker layers. Make this number 0 to use MAPS + TPC only.
0005 const int n_gas_layer = 40; // number of TPC layers
0006
0007 int Min_si_layer = 0;
0008 int Max_si_layer = n_ib_layer + n_intt_layer + n_gas_layer;
0009
0010 void SvtxInit(int verbosity = 0)
0011 {
0012   Min_si_layer = 0;
0013   Max_si_layer = n_ib_layer + n_intt_layer + n_gas_layer;
0014 }
0015
0016 double Svtx(PHG4Reco* g4Reco, double radius,
0017       const int absorberactive = 0,
0018       int verbosity = 0) {
0019
0020   float svtx_inner_radius = 2.3;
0021
0022   if (radius > svtx_inner_radius) {
0023     cout << "inconsistency: radius: " << radius
0024    << " larger than SVTX inner radius: " << svtx_inner_radius << endl;
0025     gSystem->Exit(-1);
0026   }
0027
0028   PHG4CylinderSubsystem *cyl;
0029
0030   // silicon layers ------------------------------------------------------------
0031
0032   // inner barrel
0033
0034   double ib_si_thickness[3] = {0.0050, 0.0050, 0.0050};
0035   double ib_rad[3] = {svtx_inner_radius, 3.2, 3.9};
0036   double ib_support_thickness[3] = {0.0036, 0.0036, 0.0036};
0037   double ib_length[3] = {27.0, 27.0, 27.0};
0038
0039   for (int ilayer=0;ilayer<n_ib_layer;++ilayer) {
0040     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0041     cyl->Verbosity(verbosity);
0042     radius = ib_rad[ilayer];
0043     cyl->set_double_param("radius",radius);
0044     //cyl->set_int_param("lengthviarapidity",0);
0045     cyl->set_double_param("length",ib_length[ilayer]);
0046     cyl->set_string_param("material","G4_Si");
0047     cyl->set_double_param("thickness",ib_si_thickness[ilayer]);
0048     cyl->SetActive();
0049     cyl->SuperDetector("SVTX");
0050     g4Reco->registerSubsystem( cyl );
0051
0052     radius += ib_si_thickness[ilayer] + no_overlapp;
0053
0054     cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", ilayer);
0055     cyl->Verbosity(verbosity);
0056     cyl->set_double_param("radius",radius);
0057     //cyl->set_int_param("lengthviarapidity",1);
0058     cyl->set_double_param("length",ib_length[ilayer]);
0059     cyl->set_string_param("material","G4_Cu");
0060     cyl->set_double_param("thickness",ib_support_thickness[ilayer]);
0061     cyl->SuperDetector("SVTXSUPPORT");
0062     g4Reco->registerSubsystem( cyl );
0063
0064     cout << "Added inner barrel layer with radius " << ib_rad[ilayer]
0065      << " si thickness " << ib_si_thickness[ilayer]
0066      << " support thickness " << ib_support_thickness[ilayer]
0067      << " length " << ib_length[ilayer]
0068      << endl;
0069   }
0070
0071   // intermediate tracker
0072   // parameters from RIKEN
0073   double intt_si_thickness[4] = {0.0120, 0.0120, 0.0120,0.0120};
0074   double intt_rad[4] = { 6.0, 8.0, 10.0, 12.0};
0075   // 120 microns of silicon is 0.13% of X_0, so to get 1% total we need 0.87% more in the Cu
0076   double multiplier = 0.87;  // how many times 1% do you want?
0077   double apercent = 0.0144;  // Cu thickness in cm corresponding to 1% X_0
0078   double intt_support_thickness[4] = {apercent*multiplier, apercent*multiplier, apercent*multiplier, apercent*multiplier};
0079   double intt_length[4] = {50.0, 50.0, 50.0, 50.0};
0080
0081   for (int ilayer=n_ib_layer;ilayer<n_intt_layer+n_ib_layer;++ilayer) {
0082     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0083     cyl->Verbosity(verbosity);
0084     radius = intt_rad[ilayer-n_ib_layer];
0085     cyl->set_double_param("radius",radius);
0086     cyl->set_int_param("lengthviarapidity",1);
0087     //cyl->set_double_param("length",intt_length[ilayer-n_ib_layer]);
0088     cyl->set_string_param("material","G4_Si");
0089     cyl->set_double_param("thickness",intt_si_thickness[ilayer-n_ib_layer]);
0090     cyl->SetActive();
0091     cyl->SuperDetector("SVTX");
0092     g4Reco->registerSubsystem( cyl );
0093
0094     radius += intt_si_thickness[ilayer-n_ib_layer] + no_overlapp;
0095
0096     cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", ilayer);
0097     cyl->Verbosity(verbosity);
0098     cyl->set_double_param("radius",radius);
0099     cyl->set_int_param("lengthviarapidity",1);
0100     //cyl->set_double_param("length", intt_length[ilayer-n_ib_layer]);
0101     cyl->set_string_param("material","G4_Cu");
0102     cyl->set_double_param("thickness",intt_support_thickness[ilayer-n_ib_layer]);
0103     cyl->SuperDetector("SVTXSUPPORT");
0104     g4Reco->registerSubsystem( cyl );
0105
0106     cout << "Added intermediate tracker layer with radius " << intt_rad[ilayer-n_ib_layer]
0107      << " si thickness " << intt_si_thickness[ilayer-n_ib_layer]
0108      << " support thickness " << intt_support_thickness[ilayer-n_ib_layer]
0109      << " length " << intt_length[ilayer-n_ib_layer]
0110      << endl;
0111   }
0112
0113   //TPC
0114   radius = AddTPC2G4Geo(g4Reco,radius,0);
0115
0116   return radius;
0117 }
0118
0119 void Svtx_Cells(int verbosity = 0)
0120 {
0121   // runs the cellularization of the energy deposits (g4hits)
0122   // into detector hits (g4cells)
0123
0124   //---------------
0125   // Load libraries
0126   //---------------
0127
0128   gSystem->Load("libfun4all.so");
0129   gSystem->Load("libg4detectors.so");
0130
0131   //---------------
0132   // Fun4All server
0133   //---------------
0134
0135   Fun4AllServer *se = Fun4AllServer::instance();
0136
0137   //-----------
0138   // cells
0139   //-----------
0140
0141   // inner barrel
0142   double svxcellsizex[3] = {0.0020, 0.0020, 0.0020};
0143   double svxcellsizey[3] = {0.0020, 0.0020, 0.0020};
0144
0145   // intermediate tracker
0146   double intt_cellsizex[4] = { 0.0080, 0.0080, 0.0080, 0.0080}; // cm
0147   double intt_cellsizey[4] = { 1.2, 1.2, 1.2, 1.2}; // cm
0148
0149   PHG4CylinderCellTPCReco *svtx_cells = new PHG4CylinderCellTPCReco(n_ib_layer+n_intt_layer);
0150   svtx_cells->Verbosity(1);
0151   svtx_cells->Detector("SVTX");
0152
0153   for (int i=0;i<n_ib_layer;++i) {
0154     svtx_cells->cellsize(i, svxcellsizex[i], svxcellsizey[i]);
0155     svtx_cells->set_timing_window(i, -2000.0, +2000.0);
0156   }
0157   for (int i=n_ib_layer;i<n_ib_layer+n_intt_layer;++i) {
0158     svtx_cells->cellsize(i, intt_cellsizex[i-n_ib_layer], intt_cellsizey[i-n_ib_layer]);
0159     svtx_cells->set_timing_window(i, -50.0, +50.0);
0160   }
0161
0162   // TPC
0163   AddTPC2CellReco(svtx_cells);
0164
0165   se->registerSubsystem(svtx_cells);
0166   return;
0167 }
0168
0169 void Svtx_Reco(int verbosity = 0)
0170 {
0171   //---------------
0172   // Load libraries
0173   //---------------
0174
0175   gSystem->Load("libfun4all.so");
0176   gSystem->Load("libg4hough.so");
0177
0178   //---------------
0179   // Fun4All server
0180   //---------------
0181
0182   Fun4AllServer *se = Fun4AllServer::instance();
0183
0184   //----------------------------------
0185   // Digitize the cell energy into ADC
0186   //----------------------------------
0187   PHG4SvtxDigitizer* digi = new PHG4SvtxDigitizer();
0188   digi->Verbosity(0);
0189   for (int i=0;i<n_ib_layer+n_intt_layer;++i) {
0190     digi->set_adc_scale(i, 255, 1.0e-6);
0191   }
0192   se->registerSubsystem( digi );
0193
0194   //-------------------------------------
0195   // Apply Live Area Inefficiency to Hits
0196   //-------------------------------------
0197   // defaults to 1.0 (fully active)
0198
0199   PHG4SvtxDeadArea* deadarea = new PHG4SvtxDeadArea();
0200   deadarea->Verbosity(verbosity);
0201   for(int i=0;i<n_ib_layer;i++)
0202     deadarea->set_hit_efficiency(i,0.99);
0203
0204   for(int i=n_ib_layer;i<n_ib_layer+n_intt_layer;i++)
0205     deadarea->set_hit_efficiency(i-n_ib_layer,0.99);
0206
0207   se->registerSubsystem( deadarea );
0208
0209   //-----------------------------
0210   // Apply MIP thresholds to Hits
0211   //-----------------------------
0212
0213   PHG4SvtxThresholds* thresholds = new PHG4SvtxThresholds();
0214   thresholds->Verbosity(verbosity);
0215   for(int i=0;i<n_ib_layer;i++)
0216     thresholds->set_threshold(i,0.25);  // reduce to 0.1 for increased efficiency
0217
0218   for(int i=n_ib_layer;i<n_ib_layer+n_intt_layer;i++)
0219     thresholds->set_threshold(i,0.25);
0220
0221   se->registerSubsystem( thresholds );
0222
0223   //-------------
0224   // Cluster Hits
0225   //-------------
0226
0227   PHG4SvtxClusterizer* clusterizer = new PHG4SvtxClusterizer("PHG4SvtxClusterizer",Min_si_layer,n_ib_layer+n_intt_layer-1);
0228   clusterizer->set_threshold(0.25);  // reduced from 0.5, should be same as cell threshold, since many hits are single cell
0229   se->registerSubsystem( clusterizer );
0230
0231   //---------------------
0232   // Track reconstruction
0233   //---------------------
0234   PHG4HoughTransformTPC* hough = new PHG4HoughTransformTPC(Max_si_layer,Max_si_layer-30);
0235   hough->set_mag_field(1.4);
0236   hough->setPtRescaleFactor(1.00/0.993892);
0237   hough->set_use_vertex(true);
0238   hough->setRemoveHits(true);
0239   hough->setRejectGhosts(true);
0240   hough->set_min_pT(0.2);
0241   hough->set_chi2_cut_full( 2.0 );
0242   hough->set_chi2_cut_init( 2.0 );
0243
0244   hough->setBinScale(1.0);
0245   hough->setZBinScale(1.0);
0246
0247   hough->Verbosity(verbosity);
0248
0249   double mat_scale = 1.0;
0250   // maps inner barrel, total of 0.3% of X_0 per layer
0251   for(int i=0;i<n_ib_layer;i++)
0252     hough->set_material(i, mat_scale*0.003);
0253   // intermediate tracker, total 1% of X_0 pair layer
0254   for(int i=n_ib_layer;i<n_ib_layer+n_intt_layer;i++)
0255     hough->set_material(i-n_ib_layer, mat_scale*0.010);
0256
0257   // silicon
0258   for (int i=0;i<n_ib_layer+n_intt_layer;++i) {
0259     hough->setVoteErrorScale(i, 1.0);
0260   }
0261   for (int i=0;i<n_ib_layer+n_intt_layer;++i) {
0262     hough->setFitErrorScale(i, 1./sqrt(12.));
0263   }
0264
0265   //TPC
0266   AddTPC2Reco(digi,hough);
0267
0268   se->registerSubsystem( hough );
0269
0270   //-----------------------
0271   // Momentum Recalibration
0272   //-----------------------
0273   TF1 *corr = new TF1("corr","1.0/(1+0.00908642+5.91337e-05*x+-1.87201e-05*x*x+-3.31928e-06*x*x*x+1.03004e-07*x*x*x*x+-1.05111e-09*x*x*x*x*x)",0.0,40.0);
0274   PHG4SvtxMomentumRecal* recal = new PHG4SvtxMomentumRecal("PHG4SvtxMomentumRecal",corr);
0275   se->registerSubsystem(recal);
0276
0277   //------------------
0278   // Track Projections
0279   //------------------
0280   PHG4SvtxTrackProjection* projection = new PHG4SvtxTrackProjection();
0281   projection->Verbosity(verbosity);
0282   se->registerSubsystem( projection );
0283
0284   //----------------------
0285   // Beam Spot Calculation
0286   //----------------------
0287   PHG4SvtxBeamSpotReco* beamspot = new PHG4SvtxBeamSpotReco();
0288   beamspot->Verbosity(verbosity);
0289   se->registerSubsystem( beamspot );
0290
0291   return;
0292 }
0293
0294 void G4_Svtx_Reco()
0295 {
0296   cout << "\033[31;1m"
0297        << "Warning: G4_Svtx_Reco() was moved to G4_Svtx.C and renamed to Svtx_Reco(), please update macros"
0298        << "\033[0m" << endl;
0299   Svtx_Reco();
0300
0301   return;
0302 }
0303
0304 void Svtx_Eval(std::string outputfile, int verbosity = 0)
0305 {
0306   //---------------
0307   // Load libraries
0308   //---------------
0309
0310   gSystem->Load("libfun4all.so");
0311   gSystem->Load("libg4detectors.so");
0312   gSystem->Load("libg4hough.so");
0313   gSystem->Load("libg4eval.so");
0314
0315   //---------------
0316   // Fun4All server
0317   //---------------
0318
0319   Fun4AllServer *se = Fun4AllServer::instance();
0320
0321   //----------------
0322   // SVTX evaluation
0323   //----------------
0324
0325   SvtxEvaluator* eval = new SvtxEvaluator("SVTXEVALUATOR", outputfile.c_str());
0326   eval->do_cluster_eval(true);   // make cluster ntuple
0327   eval->do_g4hit_eval(false);     // make g4hit ntuple
0328   eval->do_hit_eval(false);         // make hit ntuple
0329   eval->do_gpoint_eval(false);
0330   //eval->scan_for_embedded(true);  // evaluator will only collect embedded tracks - it will also ignore decay tracks from embedded particles!
0331   eval->scan_for_embedded(false); // evaluator takes all tracks
0332   eval->Verbosity(verbosity);
0333   se->registerSubsystem( eval );
0334
0335   // MomentumEvaluator* eval = new MomentumEvaluator(outputfile.c_str(),0.2,0.4,Max_si_layer,2,Max_si_layer-4,10.,80.);
0336   // se->registerSubsystem( eval );
0337
0338   return;
0339 }