macros/leptoquarks_fun4all/G4_Svtx_maps_ladders%2Bintt_ladders%2Btpc_KalmanPatRec.C

0001 #include <vector>
0002
0003 // if true, refit tracks with primary vertex included in track fit  - good for analysis of prompt tracks only
0004 // Adds second node to node tree, keeps original track node undisturbed
0005 // Adds second evaluator to process refitted tracks and outputs separate ntuples
0006 bool use_primary_vertex = false;
0007
0008 const int n_maps_layer = 3;  // must be 0-3, setting it to zero removes MVTX completely, n < 3 gives the first n layers
0009 const int n_intt_layer = 4;  // must be 0-4, setting this to zero will remove the INTT completely, n < 4 gives you the first n layers
0010
0011 int n_tpc_layer_inner = 8;
0012 double tpc_layer_thick_inner = 1.25;
0013 int tpc_layer_rphi_count_inner = 1920;
0014
0015 // make inner layers a factor of 2 thinner
0016 /*
0017 int n_tpc_layer_inner = 8 * 2;
0018 double tpc_layer_thick_inner = 1.25 / 2.0;
0019 int tpc_layer_rphi_count_inner = 1920 / 2;
0020 */
0021
0022 int n_tpc_layer_mid = 16;
0023 double tpc_layer_thick_mid = 1.25;
0024 int tpc_layer_rphi_count_mid = 1536;
0025
0026 int n_tpc_layer_outer = 16;
0027 double tpc_layer_thick_outer = 1.125; // outer later reach from 60-78 cm (instead of 80 cm), that leads to radial thickness of 1.125 cm
0028 int tpc_layer_rphi_count_outer = 2304;
0029
0030 int n_gas_layer = n_tpc_layer_inner + n_tpc_layer_mid + n_tpc_layer_outer;
0031
0032 double inner_cage_radius = 20.;
0033 double inner_readout_radius = 30.;
0034
0035 // TPC gas parameters
0036 // These are set for a variety of gas choices...
0037 //==============================================
0038 enum TPC_Gas
0039 {
0040   Ne2K_100,
0041   Ne2K_400,
0042   NeCF4_100,
0043   NeCF4_300,
0044   NeCF4_400,
0045   ByHand
0046 };
0047 TPC_Gas ether = TPC_Gas::NeCF4_400;
0048 //TPC_Gas ether = TPC_Gas::ByHand;
0049
0050 // Data on gasses @20 C and 760 Torr from the following source:
0051 // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
0052 double Ne_dEdx = 1.56;    // keV/cm
0053 double CF4_dEdx = 7.00;   // keV/cm
0054 double iBut_dEdx = 5.93;  // keV/cm
0055
0056 double Ne_NPrimary = 12;    // Number/cm
0057 double CF4_NPrimary = 51;   // Number/cm
0058 double iBut_NPrimary = 84;  // Number/cm
0059
0060 double Ne_NTotal = 43;     // Number/cm
0061 double CF4_NTotal = 100;   // Number/cm
0062 double iBut_NTotal = 195;  // Number/cm
0063
0064 // TPC Performance Parameter (applies extra smear to mimic the avalanche):
0065 double TPC_SigmaT = 0.03;  // 0.03 means 300 microns...Prakhar Garg Simulation...desire measurement...
0066
0067 // to be overwritten...
0068 double TPCDriftVelocity;
0069 double TPC_Trans_Diffusion;
0070 double TPC_Long_Diffusion;
0071 double TPC_dEdx;
0072 double TPC_NPri;
0073 double TPC_NTot;
0074 double TPC_ElectronsPerKeV;
0075
0076 // TPC readout shaping time and ADC clock parameters
0077 // these set the Z size of the TPC cells
0078 // These need to be set in the init since some of them require the drift velocity...
0079 //=======================================
0080 double TPCADCClock;
0081 double TPCShapingRMSLead;
0082 double TPCShapingRMSTail;
0083 double tpc_cell_z;
0084 double TPC_SmearRPhi;
0085 double TPC_SmearZ;
0086
0087 int Max_si_layer;
0088
0089 void SvtxInit(int n_TPC_layers = 40, int verbosity = 0)
0090 {
0091   Max_si_layer = n_maps_layer + n_intt_layer + n_gas_layer;
0092
0093   switch (ether)
0094   {
0095   // https://www.phenix.bnl.gov/WWW/p/draft/prakhar/tpc/HTML_Gas_Linear/Ne_CF4_IC4H10_95_3_2.html
0096   case TPC_Gas::Ne2K_100:
0097   {
0098     if (verbosity)
0099       cout << "Gas Choice:  TPC_Gas::Ne2K_100" << endl;
0100     TPCDriftVelocity = 3.2 / 1000.0;  // cm/ns
0101     TPC_Trans_Diffusion = 0.0065;     // cm/SQRT(cm)
0102     TPC_Long_Diffusion = 0.0300;      // cm/SQRT(cm)
0103     TPC_dEdx = 0.95 * Ne_dEdx + 0.03 * CF4_dEdx + 0.02 * iBut_dEdx;
0104     TPC_NPri = 0.95 * Ne_NPrimary + 0.03 * CF4_NPrimary + 0.02 * iBut_NPrimary;
0105     TPC_NTot = 0.95 * Ne_NTotal + 0.03 * CF4_NTotal + 0.02 * iBut_NTotal;
0106     break;
0107   }
0108   case TPC_Gas::Ne2K_400:
0109   {
0110     if (verbosity)
0111       cout << "Gas Choice:  TPC_Gas::Ne2K_400" << endl;
0112     TPCDriftVelocity = 5.5 / 1000.0;  // cm/ns
0113     TPC_Trans_Diffusion = 0.0120;     // cm/SQRT(cm)
0114     TPC_Long_Diffusion = 0.0175;      // cm/SQRT(cm)
0115     TPC_dEdx = 0.95 * Ne_dEdx + 0.03 * CF4_dEdx + 0.02 * iBut_dEdx;
0116     TPC_NPri = 0.95 * Ne_NPrimary + 0.03 * CF4_NPrimary + 0.02 * iBut_NPrimary;
0117     TPC_NTot = 0.95 * Ne_NTotal + 0.03 * CF4_NTotal + 0.02 * iBut_NTotal;
0118     break;
0119   }
0120   // https://www.phenix.bnl.gov/WWW/p/draft/prakhar/tpc/HTML_Gas_Linear/Ne_CF4_90_10.html
0121   case TPC_Gas::NeCF4_100:
0122   {
0123     if (verbosity)
0124       cout << "Gas Choice:  TPC_Gas::NeCF4_100" << endl;
0125     TPCDriftVelocity = 4.0 / 1000.0;  // cm/ns
0126     TPC_Trans_Diffusion = 0.0045;     // cm/SQRT(cm)
0127     TPC_Long_Diffusion = 0.0270;      // cm/SQRT(cm)
0128     TPC_dEdx = 0.90 * Ne_dEdx + 0.10 * CF4_dEdx;
0129     TPC_NPri = 0.90 * Ne_NPrimary + 0.10 * CF4_NPrimary;
0130     TPC_NTot = 0.90 * Ne_NTotal + 0.10 * CF4_NTotal;
0131     break;
0132   }
0133   case TPC_Gas::NeCF4_300:
0134   {
0135     if (verbosity)
0136       cout << "Gas Choice:  TPC_Gas::NeCF4_300" << endl;
0137     TPCDriftVelocity = 7.0 / 1000.0;  // cm/ns
0138     TPC_Trans_Diffusion = 0.0052;     // cm/SQRT(cm)
0139     TPC_Long_Diffusion = 0.0170;      // cm/SQRT(cm)
0140     TPC_dEdx = 0.90 * Ne_dEdx + 0.10 * CF4_dEdx;
0141     TPC_NPri = 0.90 * Ne_NPrimary + 0.10 * CF4_NPrimary;
0142     TPC_NTot = 0.90 * Ne_NTotal + 0.10 * CF4_NTotal;
0143     break;
0144   }
0145   case TPC_Gas::NeCF4_400:
0146   {
0147     if (verbosity)
0148       cout << "Gas Choice:  TPC_Gas::NeCF4_400" << endl;
0149     TPCDriftVelocity = 8.0 / 1000.0;  // cm/ns
0150     TPC_Trans_Diffusion = 0.0060;     // cm/SQRT(cm)
0151     TPC_Long_Diffusion = 0.0150;      // cm/SQRT(cm)
0152     TPC_dEdx = 0.90 * Ne_dEdx + 0.10 * CF4_dEdx;
0153     TPC_NPri = 0.90 * Ne_NPrimary + 0.10 * CF4_NPrimary;
0154     TPC_NTot = 0.90 * Ne_NTotal + 0.10 * CF4_NTotal;
0155     break;
0156   }
0157   case TPC_Gas::ByHand:
0158   {
0159     if (verbosity)
0160       cout << "Gas Choice:  TPC_Gas::ByHand" << endl;
0161     TPCDriftVelocity = 6.0 / 1000.0;  // cm/ns
0162     TPC_Trans_Diffusion = 0.0130;     // cm/SQRT(cm)
0163     TPC_Long_Diffusion = 0.0130;      // cm/SQRT(cm)
0164     TPC_ElectronsPerKeV = 28.0;
0165     TPC_dEdx = 0.90 * Ne_dEdx + 0.10 * CF4_dEdx;
0166     TPC_NPri = 0.90 * Ne_NPrimary + 0.10 * CF4_NPrimary;
0167     TPC_NTot = TPC_ElectronsPerKeV * TPC_dEdx;
0168     break;
0169   }
0170   default:  // defaults to NeCF4_400
0171   {
0172     if (verbosity)
0173       cout << "Gas Choice Undefined...using TPC_Gas::NeCF4_400" << endl;
0174     TPCDriftVelocity = 8.0 / 1000.0;  // cm/ns
0175     TPC_Trans_Diffusion = 0.0060;     // cm/SQRT(cm)
0176     TPC_Long_Diffusion = 0.0150;      // cm/SQRT(cm)
0177     TPC_dEdx = 0.90 * Ne_dEdx + 0.10 * CF4_dEdx;
0178     TPC_NPri = 0.90 * Ne_NPrimary + 0.10 * CF4_NPrimary;
0179     TPC_NTot = 0.90 * Ne_NTotal + 0.10 * CF4_NTotal;
0180     break;
0181   }
0182   }
0183
0184   TPC_ElectronsPerKeV = TPC_NTot / TPC_dEdx;
0185
0186   // TPC readout shaping time and ADC clock parameters
0187   // these set the Z size of the TPC cells
0188   //=======================================
0189   // TPCShapingRMSLead = 32.0;  // ns, rising RMS equivalent of shaping amplifier for 80 ns SAMPA
0190   // TPCShapingRMSTail = 48.0;  // ns, falling RMS equivalent of shaping amplifier for 80 ns SAMPA
0191   TPCADCClock = 53.0;                           // ns, corresponds to an ADC clock rate of 18.8 MHz
0192   TPCShapingRMSLead = 16.0;                     // ns, rising RMS equivalent of shaping amplifier for 40 ns SAMPA
0193   TPCShapingRMSTail = 24.0;                     // ns, falling RMS equivalent of shaping amplifier for 40 ns SAMPA
0194                                                 // TPCADCClock = 27.0;  // ns, corresponds to an ADC clock rate of 18.8 MHz * 2
0195   tpc_cell_z = TPCADCClock * TPCDriftVelocity;  // cm
0196
0197   //  TKH does not understand the physical origin of these parameters.
0198   //  however, their impact seems quite small...
0199   TPC_SmearRPhi = 0.10;
0200   TPC_SmearZ = 0.15;
0201 }
0202
0203 double Svtx(PHG4Reco* g4Reco, double radius,
0204             const int absorberactive = 0,
0205             int verbosity = 0)
0206 {
0207   if (n_maps_layer > 0)
0208   {
0209     bool maps_overlapcheck = false;  // set to true if you want to check for overlaps
0210
0211     /*
0212     The numbers used in the macro below are from the xml file dump of ITS.gdml
0213     As a sanity check, I got numbers  from Walt Sondheim's drawings, sent by email June 20, 2017:
0214     OD of Be beam pipe is 41.53 mm, ID is 40 mm
0215     Layer 0: radius 23.44 mm to sensor center, tilt from normal to radial vector:  17.37 degrees (0.303 rad), spacing btw sensor centers: 30 deg, arc spacing 12.27 mm
0216     Layer 1: radius 31.54 mm to sensor center, ttilt from normal to radial vector:  17.53 degrees (0.306 rad), spacing btw sensor centers: 22.5 deg, arc spacing 12.38 mm
0217     Layer 2: radius 39.29 to sensor center, tilt from normal to radial vector: 17.02 degrees (0.297 rad), spacing btw sensor centers: 18.0 deg, arc spacing 12.34 mm
0218     These are in reasonable agreement with the numbers I extracted from the gdml file, which are what we use below.
0219     These use a spacing in arc length of 12.37 mm and a tilt of 0.304 for all of the first three layers
0220       */
0221
0222     // MAPS inner barrel layers
0223     //======================================================
0224
0225     double maps_layer_radius[3] = {23.4, 31.5, 39.3};  // mm  - precise numbers from ITS.gdml
0226     //double maps_layer_radius[3] = {24.9, 33.0, 40.8};   // mm  - precise numbers from ITS.gdml + 1.5 mm for greater clearance from beam pipe
0227
0228     // type 1 = inner barrel stave, 2 = middle barrel stave, 3 = outer barrel stave
0229     // we use only type 0 here
0230     int stave_type[3] = {0, 0, 0};
0231     int staves_in_layer[3] = {12, 16, 20};       // Number of staves per layer in sPHENIX MVTX
0232     double phi_tilt[3] = {0.304, 0.304, 0.304};  // radians, from the gdml file, 0.304 radians is 17.4 degrees
0233
0234     for (int ilayer = 0; ilayer < n_maps_layer; ilayer++)
0235     {
0236       if (verbosity)
0237         cout << "Create Maps layer " << ilayer << " with radius " << maps_layer_radius[ilayer] << " mm, stave type " << stave_type[ilayer]
0238              << " pixel size 30 x 30 microns "
0239              << " active pixel thickness 0.0018 microns" << endl;
0240
0241       PHG4MapsSubsystem* lyr = new PHG4MapsSubsystem("MAPS", ilayer, stave_type[ilayer]);
0242       lyr->Verbosity(verbosity);
0243
0244       lyr->set_double_param("layer_nominal_radius", maps_layer_radius[ilayer]);  // thickness in cm
0245       lyr->set_int_param("N_staves", staves_in_layer[ilayer]);                   // uses fixed number of staves regardless of radius, if set. Otherwise, calculates optimum number of staves
0246
0247       // The cell size is used only during pixilization of sensor hits, but it is convemient to set it now because the geometry object needs it
0248       lyr->set_double_param("pixel_x", 0.0030);          // pitch in cm
0249       lyr->set_double_param("pixel_z", 0.0030);          // length in cm
0250       lyr->set_double_param("pixel_thickness", 0.0018);  // thickness in cm
0251       lyr->set_double_param("phitilt", phi_tilt[ilayer]);
0252
0253       lyr->set_int_param("active", 1);
0254       lyr->OverlapCheck(maps_overlapcheck);
0255
0256       lyr->set_string_param("stave_geometry_file",
0257                             string(getenv("CALIBRATIONROOT")) + string("/Tracking/geometry/ALICE_ITS_tgeo.gdml"));
0258
0259       g4Reco->registerSubsystem(lyr);
0260
0261       radius = maps_layer_radius[ilayer];
0262     }
0263   }
0264
0265   if (n_intt_layer > 0)
0266   {
0267     //-------------------
0268     // INTT ladders
0269     //-------------------
0270
0271     bool intt_overlapcheck = false;  // set to true if you want to check for overlaps
0272
0273     // instantiate the Silicon tracker subsystem and register it
0274     // We make one instance of PHG4TrackerSubsystem for all four layers of tracker
0275     // dimensions are in mm, angles are in radians
0276
0277     // PHG4SiliconTrackerSubsystem creates the detetor layer using PHG4SiliconTrackerDetector
0278     // and instantiates the appropriate PHG4SteppingAction
0279     const double intt_radius_max = 140.;  // including stagger radius (mm)
0280
0281     // The length of vpair is used to determine the number of layers
0282     std::vector<std::pair<int, int>> vpair;  // (sphxlayer, inttlayer)
0283     for (int i = 0; i < n_intt_layer; i++)
0284     {
0285       // We want the sPHENIX layer numbers for the INTT to be from n_maps_layer to n_maps_layer+n_intt_layer - 1
0286       vpair.push_back(std::make_pair(n_maps_layer + i, i));  // sphxlayer=n_maps_layer+i corresponding to inttlayer=i
0287       if (verbosity) cout << "Create strip tracker layer " << vpair[i].second << " as  sphenix layer  " << vpair[i].first << endl;
0288     }
0289     PHG4SiliconTrackerSubsystem* sitrack = new PHG4SiliconTrackerSubsystem("SILICON_TRACKER", vpair);
0290     sitrack->Verbosity(verbosity);
0291     sitrack->SetActive(1);
0292     sitrack->OverlapCheck(intt_overlapcheck);
0293     g4Reco->registerSubsystem(sitrack);
0294
0295     // outer radius marker (translation back to cm)
0296     radius = intt_radius_max * 0.1;
0297   }
0298
0299   // time projection chamber layers --------------------------------------------
0300
0301   PHG4CylinderSubsystem* cyl;
0302
0303   radius = inner_cage_radius;
0304
0305   double cage_length = 211.0;  // From TPC group, gives eta = 1.1 at 78 cm
0306   double n_rad_length_cage = 1.13e-02;
0307   double cage_thickness = 28.6 * n_rad_length_cage;  // Kapton X_0 = 28.6 cm  // mocks up Kapton + carbon fiber structure
0308
0309   // inner field cage
0310   cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", n_maps_layer + n_intt_layer);
0311   cyl->set_double_param("radius", radius);
0312   cyl->set_int_param("lengthviarapidity", 0);
0313   cyl->set_double_param("length", cage_length);
0314   cyl->set_string_param("material", "G4_KAPTON");
0315   cyl->set_double_param("thickness", cage_thickness);
0316   cyl->SuperDetector("SVTXSUPPORT");
0317   cyl->Verbosity(0);
0318   g4Reco->registerSubsystem(cyl);
0319
0320   radius += cage_thickness;
0321
0322   double inner_readout_radius = 30.;
0323   if (inner_readout_radius < radius) inner_readout_radius = radius;
0324
0325   string tpcgas = "sPHENIX_TPC_Gas";  //  Ne(90%) CF4(10%) - defined in g4main/PHG4Reco.cc
0326
0327   // Layer of inert TPC gas from 20-30 cm
0328   if (inner_readout_radius - radius > 0)
0329   {
0330     cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", n_maps_layer + n_intt_layer + 1);
0331     cyl->set_double_param("radius", radius);
0332     cyl->set_int_param("lengthviarapidity", 0);
0333     cyl->set_double_param("length", cage_length);
0334     cyl->set_string_param("material", tpcgas.c_str());
0335     cyl->set_double_param("thickness", inner_readout_radius - radius);
0336     cyl->SuperDetector("SVTXSUPPORT");
0337     g4Reco->registerSubsystem(cyl);
0338   }
0339
0340   radius = inner_readout_radius;
0341
0342   double outer_radius = 78.;
0343   //int npoints = Max_si_layer - n_maps_layer - n_intt_layer;
0344   //double delta_radius =  ( outer_radius - inner_readout_radius )/( (double)npoints );
0345
0346   // Active layers of the TPC from 30-40 cm (inner layers)
0347
0348   for (int ilayer = n_maps_layer + n_intt_layer; ilayer < (n_maps_layer + n_intt_layer + n_tpc_layer_inner); ++ilayer)
0349   {
0350     if (verbosity)
0351       cout << "Create TPC gas layer " << ilayer << " with inner radius " << radius << " cm "
0352            << " thickness " << tpc_layer_thick_inner - 0.01 << " length " << cage_length << endl;
0353
0354     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0355     cyl->set_double_param("radius", radius);
0356     cyl->set_int_param("lengthviarapidity", 0);
0357     cyl->set_double_param("length", cage_length);
0358     cyl->set_string_param("material", tpcgas.c_str());
0359     cyl->set_double_param("thickness", tpc_layer_thick_inner - 0.01);
0360     cyl->SetActive();
0361     cyl->SuperDetector("SVTX");
0362     g4Reco->registerSubsystem(cyl);
0363
0364     radius += tpc_layer_thick_inner;
0365   }
0366
0367   // Active layers of the TPC from 40-60 cm (mid layers)
0368
0369   for (int ilayer = n_maps_layer + n_intt_layer + n_tpc_layer_inner; ilayer < (n_maps_layer + n_intt_layer + n_tpc_layer_inner + n_tpc_layer_mid); ++ilayer)
0370   {
0371     if (verbosity)
0372       cout << "Create TPC gas layer " << ilayer << " with inner radius " << radius << " cm "
0373            << " thickness " << tpc_layer_thick_mid - 0.01 << " length " << cage_length << endl;
0374
0375     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0376     cyl->set_double_param("radius", radius);
0377     cyl->set_int_param("lengthviarapidity", 0);
0378     cyl->set_double_param("length", cage_length);
0379     cyl->set_string_param("material", tpcgas.c_str());
0380     cyl->set_double_param("thickness", tpc_layer_thick_mid - 0.01);
0381     cyl->SetActive();
0382     cyl->SuperDetector("SVTX");
0383     g4Reco->registerSubsystem(cyl);
0384
0385     radius += tpc_layer_thick_mid;
0386   }
0387
0388   // Active layers of the TPC from 60-80 cm (outer layers)
0389
0390   for (int ilayer = n_maps_layer + n_intt_layer + n_tpc_layer_inner + n_tpc_layer_mid; ilayer < (n_maps_layer + n_intt_layer + n_tpc_layer_inner + n_tpc_layer_mid + n_tpc_layer_outer); ++ilayer)
0391   {
0392     if (verbosity)
0393       cout << "Create TPC gas layer " << ilayer << " with inner radius " << radius << " cm "
0394            << " thickness " << tpc_layer_thick_outer - 0.01 << " length " << cage_length << endl;
0395
0396     cyl = new PHG4CylinderSubsystem("SVTX", ilayer);
0397     cyl->set_double_param("radius", radius);
0398     cyl->set_int_param("lengthviarapidity", 0);
0399     cyl->set_double_param("length", cage_length);
0400     cyl->set_string_param("material", tpcgas.c_str());
0401     cyl->set_double_param("thickness", tpc_layer_thick_outer - 0.01);
0402     cyl->SetActive();
0403     cyl->SuperDetector("SVTX");
0404     g4Reco->registerSubsystem(cyl);
0405
0406     radius += tpc_layer_thick_outer;
0407   }
0408
0409   // outer field cage
0410   cyl = new PHG4CylinderSubsystem("SVTXSUPPORT", n_maps_layer + n_intt_layer + n_gas_layer);
0411   cyl->set_double_param("radius", radius);
0412   cyl->set_int_param("lengthviarapidity", 0);
0413   cyl->set_double_param("length", cage_length);
0414   cyl->set_string_param("material", "G4_KAPTON");
0415   cyl->set_double_param("thickness", cage_thickness);  // Kapton X_0 = 28.6 cm
0416   cyl->SuperDetector("SVTXSUPPORT");
0417   g4Reco->registerSubsystem(cyl);
0418
0419   radius += cage_thickness;
0420
0421   return radius;
0422 }
0423
0424 void Svtx_Cells(int verbosity = 0)
0425 {
0426   // runs the cellularization of the energy deposits (g4hits)
0427   // into detector hits (g4cells)
0428
0429   //---------------
0430   // Load libraries
0431   //---------------
0432
0433   gSystem->Load("libfun4all.so");
0434   gSystem->Load("libg4detectors.so");
0435
0436   //---------------
0437   // Fun4All server
0438   //---------------
0439
0440   Fun4AllServer* se = Fun4AllServer::instance();
0441
0442   //-----------
0443   // SVTX cells
0444   //-----------
0445
0446   if (verbosity)
0447   {
0448     cout << "  TPC Drift Velocity: " << TPCDriftVelocity << " cm/nsec" << endl;
0449     cout << "  TPC Transverse Diffusion: " << TPC_Trans_Diffusion << " cm/SQRT(cm)" << endl;
0450     cout << "  TPC Longitudinal Diffusion: " << TPC_Long_Diffusion << " cm/SQRT(cm)" << endl;
0451     cout << "  TPC dE/dx: " << TPC_dEdx << " keV/cm" << endl;
0452     cout << "  TPC N Primary: " << TPC_NPri << " electrons/cm" << endl;
0453     cout << "  TPC N Total: " << TPC_NTot << " electrons/cm" << endl;
0454     cout << "  TPC Electrons Per keV: " << TPC_ElectronsPerKeV << " electrons/keV" << endl;
0455     cout << "  TPC ADC Clock: " << TPCADCClock << " nsec" << endl;
0456     cout << "  TPC ADC Rate: " << 1000.0 / TPCADCClock << " MHZ" << endl;
0457     cout << "  TPC Shaping Lead: " << TPCShapingRMSLead << " nsec" << endl;
0458     cout << "  TPC Shaping Tail: " << TPCShapingRMSTail << " nsec" << endl;
0459     cout << "  TPC z cell " << tpc_cell_z << " cm" << endl;
0460     cout << "  TPC Smear R-Phi " << TPC_SmearRPhi << " cm" << endl;
0461     cout << "  TPC Smear Z " << TPC_SmearZ << " cm" << endl;
0462   }
0463
0464   if (n_maps_layer > 0)
0465   {
0466     // MAPS cells
0467     PHG4MapsCellReco* maps_cells = new PHG4MapsCellReco("MAPS");
0468     maps_cells->Verbosity(verbosity);
0469     for (int ilayer = 0; ilayer < n_maps_layer; ilayer++)
0470     {
0471       maps_cells->set_timing_window(ilayer, -2000, 2000);
0472     }
0473     se->registerSubsystem(maps_cells);
0474   }
0475
0476   if (n_intt_layer > 0)
0477   {
0478     // INTT cells
0479     PHG4SiliconTrackerCellReco* reco = new PHG4SiliconTrackerCellReco("SILICON_TRACKER");
0480     // The timing windows are hard-coded in the INTT ladder model
0481     reco->Verbosity(verbosity);
0482     se->registerSubsystem(reco);
0483   }
0484
0485   // TPC cell sizes are defined at top of macro, this is for backward compatibility with old hits files
0486   if (n_gas_layer == 60)
0487   {
0488     TPCDriftVelocity = 6.0 / 1000.0;  // cm/ns
0489     tpc_cell_x = 0.12;
0490     tpc_cell_y = 0.17;
0491   }
0492
0493   // Main switch for TPC distortion
0494   const bool do_tpc_distortion = true;
0495   PHG4TPCSpaceChargeDistortion* tpc_distortion = NULL;
0496   if (do_tpc_distortion)
0497   {
0498     if (inner_cage_radius != 20. && inner_cage_radius != 30.)
0499     {
0500       cout << "Svtx_Cells - Fatal Error - TPC distortion required that "
0501               "inner_cage_radius is either 20 or 30 cm."
0502            << endl;
0503       exit(3);
0504     }
0505
0506     string TPC_distortion_file =
0507         string(getenv("CALIBRATIONROOT")) +
0508         Form("/Tracking/TPC/SpaceChargeDistortion/TPCCAGE_20_78_211_2.root");
0509     PHG4TPCSpaceChargeDistortion* tpc_distortion =
0510         new PHG4TPCSpaceChargeDistortion(TPC_distortion_file);
0511     //tpc_distortion -> setAccuracy(0); // option to over write default  factors
0512     //tpc_distortion -> setPrecision(0.001); // option to over write default  factors      // default is 0.001
0513   }
0514
0515   PHG4CylinderCellTPCReco* svtx_cells = new PHG4CylinderCellTPCReco(n_maps_layer + n_intt_layer);
0516   svtx_cells->Detector("SVTX");
0517   svtx_cells->setDistortion(tpc_distortion);
0518   if (n_gas_layer != 60)
0519   {
0520     svtx_cells->setDiffusionT(TPC_Trans_Diffusion);
0521     svtx_cells->setDiffusionL(TPC_Long_Diffusion);
0522     svtx_cells->setSigmaT(TPC_SigmaT);
0523
0524     svtx_cells->setShapingRMSLead(TPCShapingRMSLead * TPCDriftVelocity);
0525     svtx_cells->setShapingRMSTail(TPCShapingRMSTail * TPCDriftVelocity);
0526     // Expected cluster resolutions:
0527     //    r-phi: diffusion + GEM smearing = 750 microns, assume resolution is 20% of that => 150 microns
0528     //    Z:  amplifier shaping time (RMS 32 ns, 48 ns) and drift vel of 3 cm/microsec gives smearing of 3 x (32+48/2 = 1.2 mm, assume resolution is 20% of that => 240 microns
0529     svtx_cells->setSmearRPhi(TPC_SmearRPhi);  // additional random displacement of cloud positions wrt hits to give expected cluster resolution of 150 microns for charge at membrane
0530     svtx_cells->setSmearZ(TPC_SmearZ);        // additional random displacement of cloud positions wrt hits to give expected cluster rsolution of 240 microns for charge at membrane
0531   }
0532   else
0533   {
0534     // 60 layer tune
0535     svtx_cells->setDiffusionT(0.0120);
0536     svtx_cells->setDiffusionL(0.0120);
0537     svtx_cells->setSmearRPhi(0.09);  // additional smearing of cluster positions
0538     svtx_cells->setSmearZ(0.06);     // additional smearing of cluster positions
0539   }
0540   svtx_cells->set_drift_velocity(TPCDriftVelocity);
0541   svtx_cells->setHalfLength(105.5);
0542   svtx_cells->setElectronsPerKeV(TPC_ElectronsPerKeV);
0543   svtx_cells->Verbosity(0);
0544
0545   // The maps cell size is set when the detector is constructed because it is needed by the geometry object
0546   // The INTT ladder cell size is set in the detector construction code
0547
0548   // set cylinder cell TPC cell sizes
0549   //======================
0550
0551   double tpc_timing_window = 105.5 / TPCDriftVelocity;  // half length in cm / Vd in cm/ns => ns
0552
0553   // inner layers
0554   double radius_layer = inner_readout_radius + tpc_layer_thick_inner / 2.0;
0555   int counter_layer = 0;
0556   for (int i = n_maps_layer + n_intt_layer; i < n_maps_layer + n_intt_layer + n_tpc_layer_inner; i++)
0557   {
0558     // this calculates the radius at the middle of the layer
0559     if (counter_layer > 0)
0560       radius_layer += tpc_layer_thick_inner;
0561     double tpc_cell_rphi = 2 * TMath::Pi() * radius_layer / (double) tpc_layer_rphi_count_inner;
0562     svtx_cells->cellsize(i, tpc_cell_rphi, tpc_cell_z);
0563     svtx_cells->set_timing_window(i, -tpc_timing_window, +tpc_timing_window);
0564     if (verbosity)
0565       cout << "TPC cells inner: layer " << i << " center radius " << radius_layer << " tpc_cell_rphi " << tpc_cell_rphi << " tpc_cell_z " << tpc_cell_z << endl;
0566     counter_layer++;
0567   }
0568   radius_layer += tpc_layer_thick_inner / 2.0;
0569
0570   // mid layers
0571   radius_layer += tpc_layer_thick_mid / 2.0;
0572   counter_layer = 0;
0573   for (int i = n_maps_layer + n_intt_layer + n_tpc_layer_inner; i < n_maps_layer + n_intt_layer + n_tpc_layer_inner + n_tpc_layer_mid; i++)
0574   {
0575     if (counter_layer > 0)
0576       radius_layer += tpc_layer_thick_mid;
0577     double tpc_cell_rphi = 2 * TMath::Pi() * radius_layer / (double) tpc_layer_rphi_count_mid;
0578     svtx_cells->cellsize(i, tpc_cell_rphi, tpc_cell_z);
0579     svtx_cells->set_timing_window(i, -tpc_timing_window, +tpc_timing_window);
0580     if (verbosity)
0581       cout << "TPC cells mid: layer " << i << " center radius " << radius_layer << " tpc_cell_rphi " << tpc_cell_rphi << " tpc_cell_z " << tpc_cell_z << endl;
0582     counter_layer++;
0583   }
0584   radius_layer += tpc_layer_thick_mid / 2.0;
0585
0586   // outer layers
0587   radius_layer += tpc_layer_thick_outer / 2.0;
0588   counter_layer = 0;
0589   for (int i = n_maps_layer + n_intt_layer + n_tpc_layer_inner + n_tpc_layer_mid; i < n_maps_layer + n_intt_layer + n_tpc_layer_inner + n_tpc_layer_mid + n_tpc_layer_outer; i++)
0590   {
0591     if (counter_layer > 0)
0592       radius_layer += tpc_layer_thick_outer;
0593     double tpc_cell_rphi = 2 * TMath::Pi() * radius_layer / (double) tpc_layer_rphi_count_outer;
0594     svtx_cells->cellsize(i, tpc_cell_rphi, tpc_cell_z);
0595     svtx_cells->set_timing_window(i, -tpc_timing_window, +tpc_timing_window);
0596     if (verbosity)
0597       cout << "TPC cells outer: layer " << i << " center radius " << radius_layer << " tpc_cell_rphi " << tpc_cell_rphi << " tpc_cell_z " << tpc_cell_z << endl;
0598     counter_layer++;
0599   }
0600
0601   se->registerSubsystem(svtx_cells);
0602
0603   return;
0604 }
0605
0606 void Svtx_Reco(int verbosity = 0)
0607 {
0608   //---------------
0609   // Load libraries
0610   //---------------
0611
0612   gSystem->Load("libfun4all.so");
0613   gSystem->Load("libg4hough.so");
0614
0615   //---------------
0616   // Fun4All server
0617   //---------------
0618
0619   Fun4AllServer* se = Fun4AllServer::instance();
0620
0621   //----------------------------------
0622   // Digitize the cell energy into ADC
0623   //----------------------------------
0624   PHG4SvtxDigitizer* digi = new PHG4SvtxDigitizer();
0625   digi->Verbosity(0);
0626   for (int i = 0; i < n_maps_layer; ++i)
0627   {
0628     digi->set_adc_scale(i, 255, 0.4e-6);  // reduced by a factor of 2.5 when going from maps thickess of 50 microns to 18 microns
0629   }
0630
0631   if (n_intt_layer > 0)
0632   {
0633     // INTT
0634     std::vector<double> userrange;  // 3-bit ADC threshold relative to the mip_e at each layer.
0635     // these should be used for the INTT
0636     userrange.push_back(0.05);
0637     userrange.push_back(0.10);
0638     userrange.push_back(0.15);
0639     userrange.push_back(0.20);
0640     userrange.push_back(0.25);
0641     userrange.push_back(0.30);
0642     userrange.push_back(0.35);
0643     userrange.push_back(0.40);
0644
0645     PHG4SiliconTrackerDigitizer* digiintt = new PHG4SiliconTrackerDigitizer();
0646     digiintt->Verbosity(verbosity);
0647     for (int i = 0; i < n_intt_layer; i++)
0648     {
0649       digiintt->set_adc_scale(n_maps_layer + i, userrange);
0650     }
0651     se->registerSubsystem(digiintt);
0652   }
0653
0654   // TPC layers
0655   for (int i = n_maps_layer + n_intt_layer; i < Max_si_layer; ++i)
0656   {
0657     digi->set_adc_scale(i, 30000, 1.0);
0658   }
0659   se->registerSubsystem(digi);
0660
0661   //-------------------------------------
0662   // Apply Live Area Inefficiency to Hits
0663   //-------------------------------------
0664   // defaults to 1.0 (fully active)
0665
0666   PHG4SvtxDeadArea* deadarea = new PHG4SvtxDeadArea();
0667
0668   for (int i = 0; i < n_maps_layer; i++)
0669   {
0670     deadarea->Verbosity(verbosity);
0671     //deadarea->set_hit_efficiency(i,0.99);
0672     deadarea->set_hit_efficiency(i, 1.0);
0673   }
0674   for (int i = n_maps_layer; i < n_maps_layer + n_intt_layer; i++)
0675   {
0676     //deadarea->set_hit_efficiency(i,0.99);
0677     deadarea->set_hit_efficiency(i, 1.0);
0678   }
0679   se->registerSubsystem(deadarea);
0680
0681   //-----------------------------
0682   // Apply MIP thresholds to Hits
0683   //-----------------------------
0684
0685   PHG4SvtxThresholds* thresholds = new PHG4SvtxThresholds();
0686   thresholds->Verbosity(verbosity);
0687
0688   // maps
0689   for (int i = 0; i < n_maps_layer; i++)
0690   {
0691     // reduced by x2.5 when going from cylinder maps with 50 microns thickness to actual maps with 18 microns thickness
0692     // Note the non-use of set_using_thickness here, this is so that the shortest dimension of the cell sets the mip energy loss
0693     thresholds->set_threshold(i, 0.1);
0694   }
0695   // INTT
0696   for (int i = n_maps_layer; i < n_maps_layer + n_intt_layer; i++)
0697   {
0698     thresholds->set_threshold(i, 0.1);
0699     thresholds->set_use_thickness_mip(i, true);
0700   }
0701
0702   se->registerSubsystem(thresholds);
0703
0704   //-------------
0705   // Cluster Hits
0706   //-------------
0707
0708   PHG4SvtxClusterizer* clusterizer = new PHG4SvtxClusterizer("PHG4SvtxClusterizer", 0, n_maps_layer + n_intt_layer - 1);
0709   clusterizer->Verbosity(verbosity);
0710   // Reduced by 2 relative to the cylinder cell maps macro. I found this necessary to get full efficiency
0711   // Many hits in the present simulation are single cell hits, so it is not clear why the cluster threshold should be higher than the cell threshold
0712   clusterizer->set_threshold(0.1);  // fraction of a mip
0713   // no Z clustering for INTT layers (only)
0714   for (int i = n_maps_layer; i < n_maps_layer + n_intt_layer; i++)
0715   {
0716     clusterizer->set_z_clustering(i, false);
0717   }
0718
0719   se->registerSubsystem(clusterizer);
0720
0721   PHG4TPCClusterizer* tpcclusterizer = new PHG4TPCClusterizer();
0722   tpcclusterizer->Verbosity(0);
0723   tpcclusterizer->setRangeLayers(n_maps_layer + n_intt_layer, Max_si_layer);
0724   if (n_gas_layer == 40)
0725   {
0726     tpcclusterizer->setEnergyCut(12 /*15 adc*/);
0727     tpcclusterizer->setFitWindowSigmas(0.0160, 0.0160);  // should be changed when TPC cluster resolution changes
0728     tpcclusterizer->setFitWindowMax(4 /*rphibins*/, 6 /*zbins*/);
0729     tpcclusterizer->setFitEnergyThreshold(0.05 /*fraction*/);
0730   }
0731   else
0732   {
0733     // 60 layer tune
0734     tpcclusterizer->setEnergyCut(15 /*adc*/);
0735     tpcclusterizer->setFitWindowSigmas(0.0150, 0.0160);  // should be changed when TPC cluster resolution changes
0736     tpcclusterizer->setFitWindowMax(4 /*rphibins*/, 3 /*zbins*/);
0737     tpcclusterizer->setFitEnergyThreshold(0.05 /*fraction*/);
0738   }
0739   se->registerSubsystem(tpcclusterizer);
0740
0741   // This should be true for everything except testing!
0742   const bool use_kalman_pat_rec = true;
0743   if (use_kalman_pat_rec)
0744   {
0745     //---------------------
0746     // PHG4KalmanPatRec
0747     //---------------------
0748
0749     PHG4KalmanPatRec* kalman_pat_rec = new PHG4KalmanPatRec("PHG4KalmanPatRec", n_maps_layer, n_intt_layer, n_gas_layer);
0750     kalman_pat_rec->Verbosity(0);
0751     se->registerSubsystem(kalman_pat_rec);
0752   }
0753   else
0754   {
0755     //---------------------
0756     // Truth Pattern Recognition
0757     //---------------------
0758     PHG4TruthPatRec* pat_rec = new PHG4TruthPatRec();
0759     se->registerSubsystem(pat_rec);
0760   }
0761
0762   //---------------------
0763   // Kalman Filter
0764   //---------------------
0765
0766   PHG4TrackKalmanFitter* kalman = new PHG4TrackKalmanFitter();
0767   kalman->Verbosity(0);
0768   if (use_primary_vertex)
0769     kalman->set_fit_primary_tracks(true);  // include primary vertex in track fit if true
0770   se->registerSubsystem(kalman);
0771
0772   //------------------
0773   // Track Projections
0774   //------------------
0775   PHG4GenFitTrackProjection* projection = new PHG4GenFitTrackProjection();
0776   projection->Verbosity(verbosity);
0777   se->registerSubsystem(projection);
0778
0779   /*
0780   //----------------------
0781   // Beam Spot Calculation
0782   //----------------------
0783   PHG4SvtxBeamSpotReco* beamspot = new PHG4SvtxBeamSpotReco();
0784   beamspot->Verbosity(verbosity);
0785   se->registerSubsystem( beamspot );
0786   */
0787
0788   return;
0789 }
0790
0791 void G4_Svtx_Reco()
0792 {
0793   cout << "\033[31;1m"
0794        << "Warning: G4_Svtx_Reco() was moved to G4_Svtx.C and renamed to Svtx_Reco(), please update macros"
0795        << "\033[0m" << endl;
0796   Svtx_Reco();
0797
0798   return;
0799 }
0800
0801 void Svtx_Eval(std::string outputfile, int verbosity = 0)
0802 {
0803   //---------------
0804   // Load libraries
0805   //---------------
0806
0807   gSystem->Load("libfun4all.so");
0808   gSystem->Load("libg4detectors.so");
0809   gSystem->Load("libg4hough.so");
0810   gSystem->Load("libg4eval.so");
0811
0812   //---------------
0813   // Fun4All server
0814   //---------------
0815
0816   Fun4AllServer* se = Fun4AllServer::instance();
0817
0818   //----------------
0819   // SVTX evaluation
0820   //----------------
0821
0822   SvtxEvaluator* eval;
0823   eval = new SvtxEvaluator("SVTXEVALUATOR", outputfile.c_str());
0824   eval->do_cluster_eval(true);
0825   eval->do_g4hit_eval(true);
0826   eval->do_hit_eval(false);  // enable to see the hits that includes the chamber physics...
0827   eval->do_gpoint_eval(false);
0828   eval->scan_for_embedded(false);  // take all tracks if false - take only embedded tracks if true
0829   eval->Verbosity(verbosity);
0830   se->registerSubsystem(eval);
0831
0832   if (use_primary_vertex)
0833   {
0834     // make a second evaluator that records tracks fitted with primary vertex included
0835     // good for analysis of prompt tracks, particularly if MVTX is not present
0836     SvtxEvaluator* evalp;
0837     evalp = new SvtxEvaluator("SVTXEVALUATOR", string(outputfile.c_str()) + "_primary_eval.root", "PrimaryTrackMap");
0838     evalp->do_cluster_eval(true);
0839     evalp->do_g4hit_eval(true);
0840     evalp->do_hit_eval(false);
0841     evalp->do_gpoint_eval(false);
0842     evalp->scan_for_embedded(true);  // take all tracks if false - take only embedded tracks if true
0843     evalp->Verbosity(0);
0844     se->registerSubsystem(evalp);
0845   }
0846
0847   // MomentumEvaluator* eval = new MomentumEvaluator(outputfile.c_str(),0.2,0.4,Max_si_layer,2,Max_si_layer-4,10.,80.);
0848   // se->registerSubsystem( eval );
0849
0850   return;
0851 }