EICMCAnalysis/analysis_eicsmear/eic_sphenix.C

0001 /**
0002    ROOT Macro using eic-smear library to 'smear' energy and momentum of MonteCarlo generated
0003    particles to reflect the performance of an EIC detector based on sPHENIX.
0004
0005    See link for documentation of eic-smear package:
0006
0007    https://wiki.bnl.gov/eic/index.php/Smearing
0008 */
0009
0010 /** Convert pseudorapidity eta to polar angle theta */
0011 float eta2theta( double eta )
0012 {
0013   float theta = 2.0 * atan( exp( -1 * eta ) );
0014   return theta;
0015 }
0016
0017 /** Build EIC detector
0018  */
0019 Smear::Detector BuildEicSphenix() {
0020
0021   gSystem->Load("libeicsmear");
0022
0023   /* CALORIMETER
0024    *
0025    * Calorimeter resolution usually given as sigma_E/E = const% + stocastic%/Sqrt{E}
0026    * EIC Smear needs absolute sigma: sigma_E = Sqrt{const*const*E*E + stoc*stoc*E}
0027    *
0028    * Genre == 1 (third argument) means only photons and electrons are smeared
0029    *
0030    * Accept.Charge( kAll ) seems not to be working properly
0031    */
0032
0033   /* Create "electron-going" (backward) electromagnetic calorimeter (PbWO4)*/
0034   Smear::Acceptance::Zone zone_eemc( eta2theta( -1.55 ), eta2theta( -4.00 ) );
0035
0036   Smear::Device eemcE(Smear::kE,  "sqrt(0.01*0.01*E*E + 0.025*0.025*E)");
0037   eemcE.Accept.SetGenre(Smear::kElectromagnetic);
0038   eemcE.Accept.AddZone(zone_eemc);
0039
0040   /* Create "central rapidity" electromagnetic calorimeter (W-SciFi) */
0041   Smear::Acceptance::Zone zone_cemc( eta2theta(  1.24 ), eta2theta( -1.55 ) );
0042
0043   Smear::Device cemcE(Smear::kE,  "sqrt(0.05*0.05*E*E + 0.16*0.16*E)");
0044   cemcE.Accept.SetGenre(Smear::kElectromagnetic);
0045   cemcE.Accept.AddZone(zone_cemc);
0046
0047   /* Create "hadron-going" (forward) electromagnetic calorimeter (PbScint) */
0048   Smear::Acceptance::Zone zone_femc( eta2theta(  4.00 ), eta2theta(  1.24 ) );
0049
0050   Smear::Device femcE(Smear::kE,  "sqrt(0.02*0.02*E*E + 0.08*0.08*E)");
0051   femcE.Accept.SetGenre(Smear::kElectromagnetic);
0052   femcE.Accept.AddZone(zone_femc);
0053
0054   /* Create "central rapidity" hadron calorimeter (inner+outer) (Fe Scint + Steel Scint) */
0055   Smear::Acceptance::Zone zone_chcal( eta2theta(  1.10 ), eta2theta( -1.10 ) );
0056
0057   Smear::Device chcalE(Smear::kE,  "sqrt(0.12*0.12*E*E + 0.81*0.81*E)");
0058   chcalE.Accept.SetGenre(Smear::kHadronic);
0059   chcalE.Accept.AddZone(zone_chcal);
0060
0061   /* Create "hadron-going" (forward) hadron calorimeter (Fe Scint) */
0062   Smear::Acceptance::Zone zone_fhcal( eta2theta(  4.00 ), eta2theta(  1.24 ) );
0063
0064   Smear::Device fhcalE(Smear::kE,  "sqrt(0.0*0.0*E*E + 0.70*0.70*E)");
0065   fhcalE.Accept.SetGenre(Smear::kHadronic);
0066   fhcalE.Accept.AddZone(zone_fhcal);
0067
0068   /* TRACKING SYSTEM */
0069   /* Create our tracking capabilities, by a combination of mometum, theta and phi Devices.
0070    * The momentum parametrization (a*p + b) gives sigma_P/P in percent.
0071    * So Multiply through by P and divide by 100 to get absolute sigma_P
0072    * Theta and Phi parametrizations give absolute sigma in miliradians
0073    *
0074    * Note: eic-smear only saves smeared parameters for 'smeared' particle, i.e. if you want to
0075    * save e.g. the 'true' energy for particles measured with the tracker, you need to smear the
0076    * energy with '0' smearing for those particles.
0077    */
0078
0079   /* Create tracking system. */
0080   Smear::Acceptance::Zone zone_tracking( eta2theta( 4 ), eta2theta( -4.00 ) );
0081
0082   /* Try parametrization based on sPHENIX Geant4 simulation. The fit to the simulations only covers pseudorapidity -2.5 to + 2.5, so be extra careful with the parametrization outside that range. */
0083   Smear::Device trackingMomentum(Smear::kP, "P * sqrt( (7.82e-3 + 4.39e-4*(-log(tan(theta/2.0)))**2 + 7.55e-4*(-log(tan(theta/2.0)))**4)**2 + ( (1.44e-3 + -5.35e-4*(-log(tan(theta/2.0)))**2 + -6.73e-5*(-log(tan(theta/2.0)))**3 + 1.37e-4*(-log(tan(theta/2.0)))**4) * P )**2 )");
0084   trackingMomentum.Accept.SetGenre(Smear::kAll);
0085   trackingMomentum.Accept.SetCharge(Smear::kCharged);
0086   trackingMomentum.Accept.AddZone(zone_tracking);
0087
0088   Smear::Device trackingTheta(Smear::kTheta, "0");
0089   trackingTheta.Accept.SetGenre(Smear::kAll);
0090   trackingTheta.Accept.SetCharge(Smear::kCharged);
0091   trackingTheta.Accept.AddZone(zone_tracking);
0092
0093   Smear::Device trackingPhi(Smear::kPhi, "0");
0094   trackingPhi.Accept.SetGenre(Smear::kAll);
0095   trackingPhi.Accept.SetCharge(Smear::kCharged);
0096   trackingPhi.Accept.AddZone(zone_tracking);
0097
0098   /* @TODO Below numbers are for tracking with BeAST. We need to get the parametrization for EIC-sPHENIX!!! */
0099   //  Smear::Device trackingMomentum(Smear::kP, "(P*P*(0.0182031 + 0.00921047*pow((-log(tan(theta/2.0))), 2) - 0.00291243*pow((-log(tan(theta/2.0))), 4) + 0.000264353*pow((-log(tan(theta/2.0))), 6)) + P*(0.209681 + 0.275144*pow((-log(tan(theta/2.0))), 2) - 0.0436536*pow((-log(tan(theta/2.0))), 4) + 0.00367412*pow((-log(tan(theta/2.0))), 6)))*0.01");
0100   //  trackingMomentum.Accept.SetCharge(Smear::kCharged);
0101   //  trackingMomentum.Accept.AddZone(zone_tracking);
0102
0103   //  Smear::Device trackingTheta(Smear::kTheta, "((1.0/(1.0*P))*(0.752935 + 0.280370*pow((-log(tan(theta/2.0))), 2) - 0.0359713*pow((-log(tan(theta/2.0))), 4) + 0.00200623*pow((-log(tan(theta/2.0))), 6)) + 0.0282315 - 0.00998623*pow((-log(tan(theta/2.0))), 2) + 0.00117487*pow((-log(tan(theta/2.0))), 4) - 0.0000443918*pow((-log(tan(theta/2.0))), 6))*0.001");
0104   //  trackingTheta.Accept.SetCharge(Smear::kCharged);
0105   //  trackingTheta.Accept.AddZone(zone_tracking);
0106
0107   //  Smear::Device trackingPhi(Smear::kPhi, "((1.0/(1.0*P))*(0.743977 + 0.753393*pow((-log(tan(theta/2.0))), 2) + 0.0634184*pow((-log(tan(theta/2.0))), 4) + 0.0128001*pow((-log(tan(theta/2.0))), 6)) + 0.0308753 + 0.0480770*pow((-log(tan(theta/2.0))), 2) - 0.0129859*pow((-log(tan(theta/2.0))), 4) + 0.00109374*pow((-log(tan(theta/2.0))), 6))*0.001");
0108   //  trackingPhi.Accept.SetCharge(Smear::kCharged);
0109   //  trackingPhi.Accept.AddZone(zone_tracking);
0110
0111   /* Create mRICH detector parameterization */
0112
0113   //dRICH parameterization is available for studies, but the current design in the 2018 Detector Design Study only used h-side mRICH, e-side mRICH, DIRC and gas RICH, so the dRICH is commented out here
0114   /*Smear::Acceptance::Zone zone_dRICH( eta2theta( 3.95 ), eta2theta( 1.24 ));
0115   Smear::ParticleID dRICH_KPi("PIDMatrixFiles/dRICH_KPiPIDMatrix.dat");
0116   Smear::ParticleID dRICH_ePi("PIDMatrixFiles/dRICH_ePiPIDMatrix.dat");
0117   dRICH_KPi.Accept.AddZone(zone_dRICH);
0118   dRICH_ePi.Accept.AddZone(zone_dRICH);*/
0119
0120   Smear::Acceptance::Zone zone_gasRICH( eta2theta( 3.95 ), eta2theta( 1.24 ));
0121   Smear::ParticleID gasRICH_KPi("PIDMatrixFiles/gasRICH_KPiPIDMatrix.dat");
0122   Smear::ParticleID gasRICH_ePi("PIDMatrixFiles/gasRICH_ePiPIDMatrix.dat");
0123   gasRICH_KPi.Accept.AddZone(zone_gasRICH);
0124   gasRICH_ePi.Accept.AddZone(zone_gasRICH);
0125
0126   Smear::Acceptance::Zone zone_hside_mRICH( eta2theta( 1.85 ), eta2theta( 1.10 ));
0127   Smear::ParticleID hside_mRICH_KPi("PIDMatrixFiles/mRICH_KPiPIDMatrix.dat");
0128   Smear::ParticleID hside_mRICH_ePi("PIDMatrixFiles/mRICH_ePiPIDMatrix.dat");
0129   hside_mRICH_KPi.Accept.AddZone(zone_hside_mRICH);
0130   hside_mRICH_ePi.Accept.AddZone(zone_hside_mRICH);
0131
0132   Smear::Acceptance::Zone zone_DIRC( eta2theta( 1.24 ), eta2theta( -1.4 ));
0133   Smear::ParticleID DIRC("PIDMatrixFiles/DIRCPIDMatrix.dat");
0134   DIRC.Accept.AddZone(zone_DIRC);
0135
0136   Smear::Acceptance::Zone zone_eside_mRICH( eta2theta( -1.4 ), eta2theta( -3.9 ));
0137   Smear::ParticleID eside_mRICH_KPi("PIDMatrixFiles/mRICH_KPiPIDMatrix.dat");
0138   Smear::ParticleID eside_mRICH_ePi("PIDMatrixFiles/mRICH_ePiPIDMatrix.dat");
0139   eside_mRICH_KPi.Accept.AddZone(zone_eside_mRICH);
0140   eside_mRICH_ePi.Accept.AddZone(zone_eside_mRICH);
0141
0142   /* Create a DETECTOR and add the devices
0143    */
0144   Smear::Detector det;
0145
0146   det.AddDevice(eemcE);
0147   det.AddDevice(cemcE);
0148   det.AddDevice(femcE);
0149
0150   det.AddDevice(chcalE);
0151   det.AddDevice(fhcalE);
0152
0153   det.AddDevice(trackingMomentum);
0154   det.AddDevice(trackingTheta);
0155   det.AddDevice(trackingPhi);
0156
0157   //det.AddDevice(dRICH_KPi);
0158   //det.AddDevice(dRICH_ePi);
0159   det.AddDevice(gasRICH_KPi);
0160   det.AddDevice(gasRICH_ePi);
0161   det.AddDevice(hside_mRICH_KPi);
0162   det.AddDevice(hside_mRICH_ePi);
0163   det.AddDevice(DIRC);
0164   det.AddDevice(eside_mRICH_KPi);
0165   det.AddDevice(eside_mRICH_ePi);
0166
0167   det.SetEventKinematicsCalculator("NM JB DA"); // The detector will calculate event kinematics from smeared values
0168
0169   return det;
0170 }