d3sigma_calc.C

Go to the documentation of this file.

// *********************************************************************
//  Xsec (3D) of single kaon production - Martti Nirkko (28th Nov 2014)
//  Compile and run in terminal:     root -l -q d3sigma_calc.C+g
// *********************************************************************

#include "code/singlekaon_xsec.cxx"
#include <TMath.h>
#include <TString.h>

// Compile using:   root -l -b -q d3sigma_calc.C+g
void d3sigma_calc() {
  
  // ***************************
  // *   STEERING PARAMETERS   *
  // ***************************
  const int verbose = 0;                // debugging mode on/off
  const int COMP = 10;                  // time complexity (runs with O(n^3))
  const int pLeptonIsUsed = 1;          // use momentum instead of kinetic energy
  const int thetaIsUsed = 1;            // use angle instead of cosine
  const int varPhiKqIsUsed = 1;         // use phi_kq as function of E_nu
  const double pi = TMath::Pi();
  
  // NOTE: All results for a changing "pleptonIsUsed" and "thetaIsUsed" give the same results within
  //       about 1%. Here, both are set to 1 because this is closest to what the Fortran code does.
  //       To use this in comparison, it is important to use the same parameters as d4sigma_calc.C !
  
  const int type = 2;                   // lepton:   1=electron, 2=muon, 3=tau
  const int reac = 3;                   // reaction: 1=NN, 2=NP, 3=PP
  
  // Phi_kq will run with energy from phikqmin to phikqmax (can be constant)
  double phikqmin, phikqmax, FRAC;      // used to set value of phi_kq [pi]
  if (varPhiKqIsUsed) {
    phikqmin = 0.50*pi;
    phikqmax = 0.55*pi;
  } else {
    printf("Please enter phi_kq as a fraction of pi: ");
    scanf("%lf", &FRAC);
    printf("Phi_kq set to fixed value %.2lf*pi...\n", FRAC);
    phikqmin = FRAC*pi;
    phikqmax = FRAC*pi;
  }
  
  double Mlep = 0.;                     // lepton mass [GeV]
  if      (type==1) Mlep = 0.510998928e-3;
  else if (type==2) Mlep = 0.1056583715;
  else if (type==3) Mlep = 1.77682;
  else {std::cout<<"ERROR: Invalid type!"<<std::endl; return;}
  
  double Mkaon = 0.;                    // kaon mass [GeV]
  if      (reac==1) Mkaon = 0.493677;
  else if (reac==2) Mkaon = 0.497614;
  else if (reac==3) Mkaon = 0.493677;
  else {std::cout<<"ERROR: Invalid reaction!"<<std::endl; return;}
  
  double Enu = 0.;                      // neutrino energy [GeV]
  singlekaon_xsec *inst = new singlekaon_xsec();
  inst->init(Enu, type, reac);
  
  double Emin = inst->threshold;
  double Emax = 3.1;
  int Esteps = (int)(20.0*(Emax-Emin));
  
  // Initialisation of variables
  int i,j,k,itenu;
  const int nsteps2 = 10*COMP, nsteps1 = 10*COMP, nsteps3 = 10*COMP;
  double varmin1, varmin2, varmin3, varmax1, varmax2, varmax3, temp;
  double binsvar2[nsteps2+1], binsvar1[nsteps1+1], binsvar3[nsteps3+1];
  double varvals2[nsteps2],   varvals1[nsteps1],   varvals3[nsteps3];
  
  double Tlep,  Elep;                   // LEPTON ENERGY [GeV]
  double Tkaon, Ekaon;                  // KAON ENERGY [GeV]
  double theta;                         // LEPTON ANGLE [rad] - FIXED VALUE IN 3D CASE !!!
  double phikq;                         // KAON ANGLE [rad]
  
  double diff3sigma;
  double diff2sigma;
  double diff1sigma;
  double sigma;                         // CROSS SECTION [cm^2]
  
  std::cout << std::endl;
  std::cout << "  E [GeV]\txsec [cm^2]" << std::endl;
  std::cout << "-----------------------------" << std::endl;
  
  std::string fpath = "./data/";
  std::string fname;
  if (varPhiKqIsUsed)   fname = "d3sigma_phi=var.out";
  else                  fname = Form("d3sigma_phi=%.2lfpi.out",FRAC);
  FILE *outfile = fopen((fpath+fname).c_str(),"w");
  fprintf(outfile, "#  E [GeV]\txsec [cm^2]\n");
  fprintf(outfile, "#-----------------------------\n");
  
  for (itenu=0; itenu<=Esteps; itenu++) {
    
    // Initialisation of energy
    Enu = Emin + itenu*0.05;
    inst->init(Enu, type, reac);
    
    // Value of phi_kq is a linear function of energy (3D case)
    phikq = phikqmin + (double)itenu/Esteps*(phikqmax-phikqmin);
    
    // Integration ranges
    varmin1 = 0.; varmax1 = Enu-Mkaon-Mlep;     // Tkaon
    // varmin2, varmax2 filled in loop below
    if (thetaIsUsed) {                          // theta
      varmin3 = 0.;
      varmax3 = pi;
    } else {                                    // cos(theta)
      varmin3 = -1.0;
      varmax3 = 1.0;
    }
    
    // Calculate edges of bins
    for (i=0; i<=nsteps1; i++) {
      binsvar1[i] = varmin1 + i*(varmax1-varmin1)/nsteps1;
    }
    // binsvar2 filled in loop below
    for (k=0; k<=nsteps3; k++) {
      if (thetaIsUsed) {
        temp = varmax3 - k*(varmax3-varmin3)/nsteps3;   // theta
        binsvar3[k] = cos(temp);                        // cos(theta)
      } else {
        binsvar3[k] = varmin3 + k*(varmax3-varmin3)/nsteps3;
      }
      if (verbose && !itenu && k)
        std::cout << "Bin #" << k << " has cos(theta) = " << 0.5*(binsvar3[k]+binsvar3[k-1])
                  << ",\t width = " << binsvar3[k]-binsvar3[k-1] << std::endl;
    }
    
    // Calculate central bin values (average of left and right edge of bin)
    for (i=0; i<nsteps1; i++) varvals1[i] = 0.5*(binsvar1[i+1]+binsvar1[i]);
    // binsvar2 filled in loop below
    for (k=0; k<nsteps3; k++) varvals3[k] = 0.5*(binsvar3[k+1]+binsvar3[k]);
    
    sigma = 0.;  // reset integral over Tkaon
    
    // **********************
    //  Integrate over Tkaon
    // **********************
    for (i=0; i<nsteps1; i++) {
      if (binsvar1[i+1]==binsvar1[i]) continue;
      Tkaon = varvals1[i];
      Ekaon = Tkaon+Mkaon;
      
      // ----------------------------------------------
      // Define binning for second integration variable
      // ----------------------------------------------
      if (pLeptonIsUsed) {                        // plep
        varmin2 = 0.;
        varmax2 = sqrt((Enu-Ekaon)*(Enu-Ekaon)-Mlep*Mlep);
      } else {                                    // Tlep
        varmin2 = 0.;
        varmax2 = Enu-Ekaon-Mlep;
      }
      for (j=0; j<=nsteps2; j++) {
        temp = varmin2 + j*(varmax2-varmin2)/nsteps2;   // plep OR Tlep
        if (pLeptonIsUsed) {
          binsvar2[j] = sqrt(temp*temp+Mlep*Mlep)-Mlep; // plep -> Tlep
        } else {
          binsvar2[j] = temp;                           // Tlep
        }
        if (verbose && !pLeptonIsUsed && !itenu && j)
          std::cout << "Bin edge #" << j << " has Tlep = " << binsvar2[j] 
                    << ",\t width = " << binsvar2[j]-binsvar2[j-1] << std::endl;
      }
      for (j=0; j<nsteps2; j++) varvals2[j] = 0.5*(binsvar2[j+1]+binsvar2[j]);
      // ----------------------------------------------
      
      diff1sigma = 0.;  // reset integral over plep
      
      // *******************************
      //  Integrate over plep (or Tlep)
      // *******************************
      for (j=0; j<nsteps2; j++) {
        if (binsvar2[j+1]==binsvar2[j]) continue;
        Tlep = varvals2[j];
        Elep = Tlep+Mlep;
        
        diff2sigma = 0.;  // reset integral over theta
        
        // ************************************
        //  Integrate over theta (or costheta)
        // ************************************
        for (k=0; k<nsteps3; k++) {
          if (binsvar3[k+1]==binsvar3[k]) continue;
          theta = acos(varvals3[k]);
          
          // Calculate 3D differential xsec
          diff3sigma = inst->diffxsec(Tlep,Tkaon,theta,phikq);
          diff3sigma *= pow(2*pi,2);
          
          // This Jacobian no longer needed, since binning is adjusted to cos(theta)
          // diff3sigma *= sin(theta);
          
          diff2sigma += diff3sigma*(binsvar3[k+1]-binsvar3[k]);
        }
        
        // Multiplication with Jacobian for transformation plep->Tlep (it is E/p)
        diff2sigma *= Elep/sqrt(Elep*Elep-Mlep*Mlep);
        
        diff1sigma += diff2sigma*(binsvar2[j+1]-binsvar2[j]);
      }
      
      sigma += diff1sigma*(binsvar1[i+1]-binsvar1[i]);
    }
    
    // Total cross section for this energy
    std::cout << Form("%12.6lf",Enu) << "\t" << sigma << std::endl;
    fprintf(outfile, "%lf\t%e\n",Enu,sigma);
  }
  
  fclose(outfile);
  std::cout << std::endl << "Output written to file: " << fpath << fname << std::endl << std::endl;
}