// $Id$ // // File: JEventProcessor_pi0_res.cc // Created: Tue Nov 1 22:48:06 EDT 2011 // Creator: davidl (on Darwin Amelia.local 9.8.0 i386) // #include "JEventProcessor_pi0_res.h" using namespace jana; #include using namespace std; #include // Routine used to create our JEventProcessor #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_pi0_res()); } } // "C" //------------------ // JEventProcessor_pi0_res (Constructor) //------------------ JEventProcessor_pi0_res::JEventProcessor_pi0_res() { // Splines used for theta dependence of tdiff resolution double zrev_all[]={12.6509, 97.862, 149.661, 283.255, 341.992}; double sigma_334[]={1.17205, 1.17967, 1, 0.405325, 0.362295}; double sigma_1234[]={1.35234, 1.28165, 1, 0.606042, 0.555303}; sp334 = new TSpline3("334", &zrev_all[0], &sigma_334[0], 5); sp1234 = new TSpline3("1234", &zrev_all[0], &sigma_1234[0], 5); } //------------------ // ~JEventProcessor_pi0_res (Destructor) //------------------ JEventProcessor_pi0_res::~JEventProcessor_pi0_res() { } //------------------ // init //------------------ jerror_t JEventProcessor_pi0_res::init(void) { // This is called once at program startup return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_pi0_res::brun(JEventLoop *eventLoop, int runnumber) { sigma_theta_scale_bcal = 1.0; sigma_phi_scale_bcal = 1.0; sigma_E_scale_bcal = 1.0; sigma_theta_scale_fcal = 1.0; sigma_phi_scale_fcal = 1.0; sigma_E_scale_fcal = 1.0; use_334 = true; // true=334 scheme false=1234 scheme gPARMS->SetDefaultParameter("SIGMA_THETA_SCALE_BCAL", sigma_theta_scale_bcal); gPARMS->SetDefaultParameter("SIGMA_PHI_SCALE_BCAL", sigma_phi_scale_bcal); gPARMS->SetDefaultParameter("SIGMA_E_SCALE_BCAL", sigma_E_scale_bcal); gPARMS->SetDefaultParameter("SIGMA_THETA_SCALE_FCAL", sigma_theta_scale_fcal); gPARMS->SetDefaultParameter("SIGMA_PHI_SCALE_FCAL", sigma_phi_scale_fcal); gPARMS->SetDefaultParameter("SIGMA_E_SCALE_FCAL", sigma_E_scale_fcal); gPARMS->SetDefaultParameter("USE_334", use_334, "Use the 334 BCAL segmentation scheme. Otherwise use 1234 scheme"); pi0mass2 = new TH1D("pi0mass2", "", 100, -0.1, 0.2); pi0mass = new TH1D("pi0mass", "", 100, 0.0, 0.5); pi0mass2_bcal = (TH1D*)pi0mass2->Clone("pi0mass2_bcal"); pi0mass_bcal = (TH1D*)pi0mass->Clone("pi0mass_bcal"); pi0mass2_fcal = (TH1D*)pi0mass2->Clone("pi0mass2_fcal"); pi0mass_fcal = (TH1D*)pi0mass->Clone("pi0mass_fcal"); pi0mass2_bcal_and_fcal = (TH1D*)pi0mass2->Clone("pi0mass2_bcal_and_fcal"); pi0mass_bcal_and_fcal = (TH1D*)pi0mass->Clone("pi0mass_bcal_and_fcal"); if(use_334){ jout< throwns; loop->Get(throwns); // Loop over thrown particles and extract gammas // storing smeared 4-vectors in separate lists // for BCAL and FCAL vector gammas_bcal; vector gammas_fcal; for(unsigned int i=0; itype != 1)continue; DVector3 mom = thrown->momentum(); double theta = mom.Theta(); double theta_degrees = theta*180.0/3.14159265; // BCAL photons if(theta_degrees>11.0 && theta_degrees<=110.0){ double Egen = thrown->energy(); double sigma_theta = GetThetaSigmaBCAL(mom, Egen); double sigma_phi = GetPhiSigmaBCAL(mom, Egen); double sigma_E = GetESigmaBCAL(mom, Egen); // Opportunity to scale errors sigma_theta *= sigma_theta_scale_bcal; sigma_phi *= sigma_phi_scale_bcal; sigma_E *= sigma_E_scale_bcal; double theta = rnd.Gaus(mom.Theta(), sigma_theta); double phi = rnd.Gaus(mom.Phi(), sigma_phi); double E = rnd.Gaus(Egen, sigma_E); double px = E*sin(theta)*cos(phi); double py = E*sin(theta)*sin(phi); double pz = E*cos(theta); TLorentzVector gamma(px, py, pz, E); gammas_bcal.push_back(gamma); } // FCAL photons double phi = mom.Phi(); double L = 650.0 - 65.0; // rough distance to FCAL from center of target in z double r = L*sin(theta); double x = r*cos(phi); double y = r*cos(phi); if(theta_degrees<=11.0 && (fabs(x)>10.0 || fabs(y)>10.0)){ double Egen = thrown->energy(); double sigma_theta = GetThetaSigmaFCAL(mom, Egen); double sigma_phi = GetPhiSigmaFCAL(mom, Egen); double sigma_E = GetESigmaFCAL(mom, Egen); sigma_theta *= sigma_theta_scale_fcal; sigma_phi *= sigma_phi_scale_fcal; sigma_E *= sigma_E_scale_fcal; double theta = rnd.Gaus(mom.Theta(), sigma_theta); double phi = rnd.Gaus(mom.Phi(), sigma_phi); double E = rnd.Gaus(Egen, sigma_E); double px = E*sin(theta)*cos(phi); double py = E*sin(theta)*sin(phi); double pz = E*cos(theta); TLorentzVector gamma(px, py, pz, E); gammas_fcal.push_back(gamma); } } // BCAL only for(unsigned int i=0; iFill(pi0.M2()); if(pi0.M2()>=0.0)pi0mass->Fill(pi0.M()); pi0mass2_bcal->Fill(pi0.M2()); if(pi0.M2()>=0.0)pi0mass_bcal->Fill(pi0.M()); } } // FCAL only for(unsigned int i=0; iFill(pi0.M2()); if(pi0.M2()>=0.0)pi0mass->Fill(pi0.M()); pi0mass2_fcal->Fill(pi0.M2()); if(pi0.M2()>=0.0)pi0mass_fcal->Fill(pi0.M()); } } // One photon in BCAL, one in FCAL for(unsigned int i=0; iFill(pi0.M2()); if(pi0.M2()>=0.0)pi0mass->Fill(pi0.M()); pi0mass2_bcal_and_fcal->Fill(pi0.M2()); if(pi0.M2()>=0.0)pi0mass_bcal_and_fcal->Fill(pi0.M()); } } return NOERROR; } //------------------ // GetThetaSigmaBCAL //------------------ double JEventProcessor_pi0_res::GetThetaSigmaBCAL(DVector3 &mom, double E) { if(use_334){ double R_granularity = (3.0+3.0)/2.0; // average summed cell size in R for inner layers return GetThetaSigmaBCALParms(mom, E, 50.2, 78.8, R_granularity); }else{ double R_granularity = (2.0+4.0+6.0)/3.0; // average summed cell size in R for inner layers return GetThetaSigmaBCALParms(mom, E, 45.6, 40.3, R_granularity); } } //------------------ // GetThetaSigmaBCALParms //------------------ double JEventProcessor_pi0_res::GetThetaSigmaBCALParms(DVector3 &mom, double E, double p0, double p1, double R_granularity) { // Calculate z from downstream end. It was done this way since the // spline did not appear to behave as nicely when the x values of the // knots were decreasing. (Go figure.) double theta = mom.Theta(); double z = 70.0/tan(theta) + 65.0; double zrev = 407.0 - z; // Calculate tdiff resolution at 20 degrees for this energy double sigma_334_20degrees = sqrt(pow(p0, 2.0)*E + pow(p1, 2.0)); // Get ratio of tdiff res at this angle to res at 20 degrees // by evaluating spline double scale = sp334->Eval(zrev); // Calcuate uncertainty in reconstructed z position double sigma_tdiff = scale*sigma_334_20degrees; // in ps double ceff = 16.75; double sigma_z = sigma_tdiff*ceff/1000.0; // factor of 1000 converts to ns // Calculate uncertainty in reconstructed R position // This should be based on the granularity in R. // Assume a factor of 5 better than the granularity. // Include 1/sqrt(E) energy dependence with a 2mm floor term double sigma_R_floor = 0.2; double sigma_R_stochastic = R_granularity/5.0; double sigma_R = sqrt(sigma_R_stochastic*sigma_R_stochastic/E + sigma_R_floor*sigma_R_floor); // Combine errors in R and z to get error in theta. z -= 65.0; // shift z so target center is at z=0 double R = 65.0 + 12.5*sin(theta); // this is an estimate double dtheta_dR = z/(z*z + R*R); double dtheta_dZ = -R/(z*z + R*R); double sigma_tot2 = pow(dtheta_dR*sigma_R, 2.0) + pow(dtheta_dZ*sigma_z, 2.0); double sigma_tot = sqrt(sigma_tot2); return sigma_tot; } #if 0 //------------------ // GetThetaSigma334 //------------------ double JEventProcessor_pi0_res::GetThetaSigma334(DVector3 &mom, double E) { // Calculate z from downstream end. It was done this way since the // spline did not appear to behave as nicely when the x values of the // knots were decreasing. (Go figure.) double theta = mom.Theta(); double z = 70.0/tan(theta) + 65.0; double zrev = 407.0 - z; // Calculate tdiff resolution at 20 degrees for this energy double p0_334 = 50.2; double p1_334 = 78.8; double sigma_334_20degrees = sqrt(pow(p0_334, 2.0)*E + pow(p1_334, 2.0)); // Get ratio of tdiff res at this angle to res at 20 degrees // by evaluating spline double scale = sp334->Eval(zrev); // Calcuate uncertainty in reconstructed z position double sigma_tdiff = scale*sigma_334_20degrees; // in ps double ceff = 16.75; double sigma_z = sigma_tdiff*ceff/1000.0; // factor of 1000 converts to ns // Calculate uncertainty in reconstructed R position // This should be based on the granularity in R. // For the 334 scheme, the granularity is // 3*2cm = 6cm for the inner cells. Again assume a factor // of 5 better than the granularity. // Include 1/sqrt(E) energy dependence with a 2mm floor term double sigma_R_floor = 0.2; double sigma_R_stochastic = (3.0+3.0)/2.0/5.0; double sigma_R = sqrt(sigma_R_stochastic*sigma_R_stochastic/E + sigma_R_floor*sigma_R_floor); // Combine errors in R and z to get error in theta. z -= 65.0; // shift z so target center is at z=0 double R = 65.0 + 12.5*sin(theta); // this is an estimate double dtheta_dR = z/(z*z + R*R); double dtheta_dZ = -R/(z*z + R*R); double sigma_tot2 = pow(dtheta_dR*sigma_R, 2.0) + pow(dtheta_dZ*sigma_z, 2.0); double sigma_tot = sqrt(sigma_tot2); return sigma_tot; } //------------------ // GetThetaSigma1234 //------------------ double JEventProcessor_pi0_res::GetThetaSigma1234(DVector3 &mom, double E) { // Calculate z from downstream end. It was done this way since the // spline did not appear to behave as nicely when the x values of the // knots were decreasing. (Go figure.) double theta = mom.Theta(); double z = 70.0/tan(theta) + 65.0; double zrev = 407.0 - z; // Calculate tdiff resolution at 20 degrees for this energy double p0_1234 = 45.6; double p1_1234 = 40.3; double sigma_1234_20degrees = sqrt(pow(p0_1234, 2.0)*E + pow(p1_1234, 2.0)); // Get ratio of tdiff res at this angle to res at 20 degrees // by evaluating spline double scale = sp1234->Eval(zrev); // Calcuate uncertainty in reconstructed z position double sigma_tdiff = scale*sigma_1234_20degrees; // in ps double ceff = 16.75; double sigma_z = sigma_tdiff*ceff/1000.0; // factor of 1000 converts to ns // Calculate uncertainty in reconstructed R position // This should be based on the granularity in R. // For the 1234 scheme, the granularity could be considered // 4cm for the inner cells since that is the average summed // cell size. Again assume a factor // of 5 better than the granularity. // Include 1/sqrt(E) energy dependence with a 2mm floor term double sigma_R_floor = 0.2; double sigma_R_stochastic = (2.0+4.0+6.0)/3.0/5.0; double sigma_R = sqrt(sigma_R_stochastic*sigma_R_stochastic/E + sigma_R_floor*sigma_R_floor); // Combine errors in R and z to get error in theta. z -= 65.0; // shift z so target center is at z=0 double R = 65.0 + 12.5*sin(theta); // this is an estimate double dtheta_dR = z/(z*z + R*R); double dtheta_dZ = -R/(z*z + R*R); double sigma_tot2 = pow(dtheta_dR*sigma_R, 2.0) + pow(dtheta_dZ*sigma_z, 2.0); double sigma_tot = sqrt(sigma_tot2); return sigma_tot; } #endif //------------------ // GetPhiSigmaBCAL //------------------ double JEventProcessor_pi0_res::GetPhiSigmaBCAL(DVector3 &mom, double E) { // There are 48*4 columns in BCAL that separate it in the phi direction. // Calculate delta of one segment in radians double delta_phi = 2.0*3.14159/(48.0*4.0); // Assume a reconstructed phi resolution that is 5 times smaller than // granularity for 1GeV particles. Also assume a 6mrad floor term // (see PDG pg 320) double sigma_stochastic = delta_phi/5.0/sqrt(E); double sigma_floor =0.006; double sigma = sqrt(sigma_stochastic*sigma_stochastic + sigma_floor*sigma_floor); return sigma; } //------------------ // GetESigmaBCAL //------------------ double JEventProcessor_pi0_res::GetESigmaBCAL(DVector3 &mom, double E) { return sqrt(pow(0.057, 2.0)*E + pow(0.023*E, 2.0)); } //------------------ // GetThetaSigmaFCAL //------------------ double JEventProcessor_pi0_res::GetThetaSigmaFCAL(DVector3 &mom, double E) { // From Alex D's Mar. 15, 2008 note double sigma_r = 0.64/sqrt(E); // in cm double L = 650.0 - 65.0; // estimate center of target to FCAL front face double r = L*sin(mom.Theta()); // distance from beamline in plane of FCAL double sigma_theta = sigma_r * L/(L*L + r*r); return sigma_theta; } //------------------ // GetPhiSigmaFCAL //------------------ double JEventProcessor_pi0_res::GetPhiSigmaFCAL(DVector3 &mom, double E) { // From Alex D's Mar. 15, 2008 note double sigma_r = 0.64/sqrt(E); // in cm double L = 650.0 - 65.0; // estimate center of target to FCAL front face double r = L*sin(mom.Theta()); // distance from beamline in plane of FCAL double sigma_phi = atan(sigma_r/r); return sigma_phi; } //------------------ // GetESigmaFCAL //------------------ double JEventProcessor_pi0_res::GetESigmaFCAL(DVector3 &mom, double E) { // These numbers from GlueX-doc-1093 fig.3. It appears that the floor // term was not added in quadrature based on reading the resolution // off the image. Thus, we do not add in quadrature here. Note that this // is much better than previously used resolutions, including the numbers // in Alex Dzierba's note from Mar. 15, 2008 that had a 7.3% stochastic // term and a 3.6% floor term (added in quadrature). return 0.056*sqrt(E) + 0.006*E; } //------------------ // erun //------------------ jerror_t JEventProcessor_pi0_res::erun(void) { // This is called whenever the run number changes, before it is // changed to give you a chance to clean up before processing // events from the next run number. return NOERROR; } //------------------ // fini //------------------ jerror_t JEventProcessor_pi0_res::fini(void) { // Called before program exit after event processing is finished. return NOERROR; }