LightFlavorRatios/kinematic_acceptance/CalculateKinematicAcceptance.C

0001 #include <phool/PHRandomSeed.h>
0002 #include <gsl/gsl_rng.h>
0003
0004 namespace Settings
0005 {
0006
0007   const std::string mother_name = "#Lambda";
0008   const float mother_mass = 1.115;
0009
0010   const std::string daughter1_name = "proton";
0011   const float daughter1_mass = .938;
0012
0013   const std::string daughter2_name = "#pi";
0014   const float daughter2_mass = .1395;
0015
0016 /*
0017   const std::string mother_name = "K_{S}^{0}";
0018   const float mother_mass = 0.4976;
0019
0020   const std::string daughter1_name = "#pi^{+}";
0021   const float daughter1_mass = .1395;
0022
0023   const std::string daughter2_name = "#pi^{-}";
0024   const float daughter2_mass = .1395;
0025 */
0026
0027   const float mother_pt_low = 0.;
0028   const float mother_pt_high = 4.;
0029   const int mother_pt_nbins = 400;
0030
0031   const float daughter_pt_low = 0.;
0032   const float daughter_pt_high = 4.;
0033   const int daughter_pt_nbins = 400;
0034
0035   const int n_mother_samples = 200;
0036   const int n_daughter_samples = 200;
0037
0038   const float min_daughter_pt = 0.160;
0039
0040   const bool debug = false;
0041 };
0042
0043 float get_daughter_pt_CM(const float mother_mass, const float daughter1_mass, const float daughter2_mass)
0044 {
0045   const float daughter1_gamma_CM = (pow(mother_mass,2.)+pow(daughter1_mass,2.)-pow(daughter2_mass,2.))/(2.*mother_mass*daughter1_mass);
0046   return sqrt(pow(daughter1_gamma_CM,2.)-1)*daughter1_mass;
0047 }
0048
0049 std::string bin_name(const int bin_number)
0050 {
0051   const float bin_pt_low = Settings::mother_pt_low+((float)bin_number)/Settings::mother_pt_nbins * (Settings::mother_pt_high-Settings::mother_pt_low);
0052   const float bin_pt_high = bin_pt_low + (1./Settings::mother_pt_nbins) * (Settings::mother_pt_high-Settings::mother_pt_low);
0053   return std::to_string(bin_pt_low) + " <= mother pT <= " + std::to_string(bin_pt_high);
0054 }
0055
0056 void CalculateKinematicAcceptance()
0057 {
0058   TFile* fout = new TFile("kinematic_acceptance.root","RECREATE");
0059
0060   TTree* outtree = new TTree("candidates","Candidates");
0061
0062   int mother_ptbin;
0063   float mother_mass = Settings::mother_mass;
0064   float mother_pt;
0065   float mother_gamma;
0066   float mother_beta;
0067   float mother_E;
0068
0069   float theta_xy;
0070   float theta_rz;
0071   float daughter_p_CM;
0072   float daughter_px_CM;
0073   float daughter_py_CM;
0074   float daughter_pz_CM;
0075
0076   float daughter1_mass = Settings::daughter1_mass;
0077   float daughter1_gamma_CM;
0078   float daughter1_E_CM;
0079   float daughter1_px_lab;
0080   float daughter1_py_lab;
0081   float daughter1_pt_lab;
0082
0083   float daughter2_mass = Settings::daughter2_mass;
0084   float daughter2_gamma_CM;
0085   float daughter2_E_CM;
0086   float daughter2_px_lab;
0087   float daughter2_py_lab;
0088   float daughter2_pt_lab;
0089
0090   outtree->Branch("mother_ptbin",&mother_ptbin);
0091   outtree->Branch("mother_mass",&mother_mass);
0092   outtree->Branch("mother_pt",&mother_pt);
0093   outtree->Branch("mother_gamma",&mother_gamma);
0094   outtree->Branch("mother_beta",&mother_beta);
0095   outtree->Branch("mother_E",&mother_E);
0096
0097   outtree->Branch("theta_xy",&theta_xy);
0098   outtree->Branch("theta_rz",&theta_rz);
0099   outtree->Branch("daughter_p_CM",&daughter_p_CM);
0100   outtree->Branch("daughter_px_CM",&daughter_px_CM);
0101   outtree->Branch("daughter_py_CM",&daughter_py_CM);
0102
0103   outtree->Branch("daughter1_mass",&daughter1_mass);
0104   outtree->Branch("daughter1_gamma_CM",&daughter1_gamma_CM);
0105   outtree->Branch("daughter1_E_CM",&daughter1_E_CM);
0106   outtree->Branch("daughter1_px_lab",&daughter1_px_lab);
0107   outtree->Branch("daughter1_py_lab",&daughter1_py_lab);
0108   outtree->Branch("daughter1_pt_lab",&daughter1_pt_lab);
0109
0110   outtree->Branch("daughter2_mass",&daughter2_mass);
0111   outtree->Branch("daughter2_gamma_CM",&daughter2_gamma_CM);
0112   outtree->Branch("daughter2_E_CM",&daughter2_E_CM);
0113   outtree->Branch("daughter2_px_lab",&daughter2_px_lab);
0114   outtree->Branch("daughter2_py_lab",&daughter2_py_lab);
0115   outtree->Branch("daughter2_pt_lab",&daughter2_pt_lab);
0116
0117   TH1F* acceptance = new TH1F("acceptance",(Settings::mother_name+" kinematic acceptance;"+Settings::mother_name+" pT [GeV];acceptance").c_str(),
0118                               Settings::mother_pt_nbins,Settings::mother_pt_low,Settings::mother_pt_high);
0119   TH1F* pt1 = new TH1F("pt1",(Settings::daughter1_name+" pT").c_str(),
0120                               Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high);
0121   TH1F* pt2 = new TH1F("pt2",(Settings::daughter2_name+" pT").c_str(),
0122                               Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high);
0123   TH2F* pt1pt2 = new TH2F("pt1pt2",("("+Settings::daughter1_name+", "+Settings::daughter2_name+") pT").c_str(),
0124                               Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high,
0125                               Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high);
0126
0127   std::vector<TH1F*> pt1_vs_motherpt;
0128   std::vector<TH1F*> pt2_vs_motherpt;
0129   std::vector<TH2F*> pt1pt2_vs_motherpt;
0130
0131   for(int i=0; i<Settings::mother_pt_nbins; i++)
0132   {
0133     std::string name_suffix = "_bin"+std::to_string(i);
0134     std::string title_suffix = bin_name(i);
0135     pt1_vs_motherpt.push_back(new TH1F(("pt1_vs_motherpt"+name_suffix).c_str(),(Settings::daughter1_name+" pT, "+title_suffix).c_str(),
0136                                         Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high));
0137     pt2_vs_motherpt.push_back(new TH1F(("pt2_vs_motherpt"+name_suffix).c_str(),(Settings::daughter2_name+" pT, "+title_suffix).c_str(),
0138                                         Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high));
0139     pt1pt2_vs_motherpt.push_back(new TH2F(("pt1pt2_vs_motherpt"+name_suffix).c_str(),("("+Settings::daughter1_name+","+Settings::daughter2_name+") pT, "+title_suffix).c_str(),
0140                                         Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high,
0141                                         Settings::daughter_pt_nbins,Settings::daughter_pt_low,Settings::daughter_pt_high));
0142   }
0143
0144   const uint seed = PHRandomSeed();
0145   std::unique_ptr<gsl_rng> m_rng;
0146   m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
0147   gsl_rng_set(m_rng.get(), seed);
0148
0149   for(int impt = 0; impt<Settings::mother_pt_nbins; impt++)
0150   {
0151     mother_ptbin = impt;
0152
0153     const float bin_pt_low = acceptance->GetBinLowEdge(impt+1);
0154     const float bin_pt_high = bin_pt_low + acceptance->GetBinWidth(impt+1);
0155
0156     std::cout << "mother pt bin " << impt << std::endl;
0157
0158     for(int ims=0; ims<Settings::n_mother_samples; ims++)
0159     {
0160       mother_pt = bin_pt_low + (bin_pt_high-bin_pt_low) * gsl_rng_uniform_pos(m_rng.get());
0161
0162       mother_gamma = sqrt(1.+pow(mother_pt/mother_mass,2.));
0163       mother_beta = mother_pt/(mother_gamma*mother_mass);
0164       mother_E = mother_gamma*mother_mass;
0165
0166       if(Settings::debug)
0167       {
0168         std::cout << "pt in (" << bin_pt_low << ", " << bin_pt_high << ")" << std::endl;
0169         std::cout << "mother pt " << mother_pt << " E " << mother_E << " gamma " << mother_gamma << " beta " << mother_beta << std::endl;
0170       }
0171
0172       for(int is=0; is<Settings::n_daughter_samples; is++)
0173       {
0174         // pick random decay axis (equivalent to random point on sphere) in a way that respects angle dependence of solid-angle element
0175         theta_xy = 2.*M_PI*gsl_rng_uniform_pos(m_rng.get());
0176         theta_rz = acos(2*gsl_rng_uniform_pos(m_rng.get())-1.);
0177
0178         daughter1_gamma_CM = (pow(Settings::mother_mass,2.)+pow(Settings::daughter1_mass,2.)-pow(Settings::daughter2_mass,2.))/(2.*Settings::mother_mass*Settings::daughter1_mass);
0179         daughter2_gamma_CM = (pow(Settings::mother_mass,2.)+pow(Settings::daughter2_mass,2.)-pow(Settings::daughter1_mass,2.))/(2.*Settings::mother_mass*Settings::daughter2_mass);
0180
0181         daughter1_E_CM = daughter1_gamma_CM * Settings::daughter1_mass;
0182         daughter2_E_CM = daughter2_gamma_CM * Settings::daughter2_mass;
0183
0184         daughter_p_CM = sqrt(pow(daughter1_gamma_CM,2.)-1.) * Settings::daughter1_mass;
0185         daughter_px_CM = daughter_p_CM * cos(theta_xy) * cos(theta_rz);
0186         daughter_py_CM = daughter_p_CM * sin(theta_xy) * cos(theta_rz);
0187         daughter_pz_CM = daughter_p_CM * sin(theta_rz);
0188
0189         if(gsl_rng_uniform_pos(m_rng.get())>0.5)
0190         {
0191           daughter1_px_lab = mother_gamma * (daughter_px_CM + mother_beta*daughter1_E_CM);
0192           daughter2_px_lab = mother_gamma * (-daughter_px_CM + mother_beta*daughter2_E_CM);
0193
0194           daughter1_py_lab = daughter_py_CM;
0195           daughter2_py_lab = -daughter_py_CM;
0196         }
0197         else
0198         {
0199           daughter1_px_lab = mother_gamma * (-daughter_px_CM + mother_beta*daughter1_E_CM);
0200           daughter2_px_lab = mother_gamma * (daughter_px_CM + mother_beta*daughter2_E_CM);
0201
0202           daughter1_py_lab = -daughter_py_CM;
0203           daughter2_py_lab = daughter_py_CM;
0204         }
0205
0206         daughter1_pt_lab = sqrt(pow(daughter1_px_lab,2.)+pow(daughter_py_CM,2.));
0207         daughter2_pt_lab = sqrt(pow(daughter2_px_lab,2.)+pow(daughter_py_CM,2.));
0208
0209         if(Settings::debug)
0210         {
0211           std::cout << "sample " << is << std::endl;
0212           std::cout << "theta_xy " << theta_xy << " theta_rz " << theta_rz << std::endl;
0213           std::cout << Settings::daughter1_name << " gammaCM " << daughter1_gamma_CM << " E_CM " << daughter1_E_CM << " p_CM " << daughter_p_CM << " pt_lab " << daughter1_pt_lab << std::endl;
0214           std::cout << Settings::daughter2_name << " gammaCM " << daughter2_gamma_CM << " E_CM " << daughter2_E_CM << " p_CM " << daughter_p_CM << " pt_lab " << daughter2_pt_lab << std::endl;
0215         }
0216
0217         outtree->Fill();
0218
0219         pt1->Fill(daughter1_pt_lab);
0220         pt2->Fill(daughter2_pt_lab);
0221         pt1pt2->Fill(daughter1_pt_lab,daughter2_pt_lab);
0222
0223         pt1_vs_motherpt[impt]->Fill(daughter1_pt_lab);
0224         pt2_vs_motherpt[impt]->Fill(daughter2_pt_lab);
0225         pt1pt2_vs_motherpt[impt]->Fill(daughter1_pt_lab,daughter2_pt_lab);
0226
0227         if(daughter1_pt_lab>Settings::min_daughter_pt && daughter2_pt_lab>Settings::min_daughter_pt)
0228         {
0229           acceptance->Fill(mother_pt);
0230         }
0231       }
0232     }
0233   }
0234   acceptance->Scale(1./(Settings::n_mother_samples*Settings::n_daughter_samples));
0235   acceptance->SetMinimum(0.);
0236
0237   outtree->Write();
0238
0239   acceptance->Write();
0240   pt1->Write();
0241   pt2->Write();
0242   pt1pt2->Write();
0243   for(int i=0; i<Settings::mother_pt_nbins; i++)
0244   {
0245     pt1_vs_motherpt[i]->Write();
0246     pt2_vs_motherpt[i]->Write();
0247     pt1pt2_vs_motherpt[i]->Write();
0248   }
0249 }