radcorr/RadCorrCalculator.h at main · clarenc3/radcorr · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#ifndef _RADCORRCALC_H_SEEN_
#define _RADCORRCALC_H_SEEN_

#include <iostream>
#include <iomanip>

#include "TGraph.h"
#include "TFile.h"

// Just a handy enum
enum NuType {kNuMu = 0, kNuMuBar = 1, kNuE = 2, kNuEBar = 3};

class RadCorrCalc {
  public:
    // Constructor and destructor
    RadCorrCalc();

    void SetupNuMu();
    void SetupNuE();
    void SetupGraphs();

    // Just delete the graphs
    ~RadCorrCalc() {
      // Delete the numu(bar) graphs
      for (int i = 0; i < kNuMuBar+1; ++i) {
        for (int j = 0; j < nEnuNuMu; ++j) {
          if (GraphsNuMu[i][j] != NULL) delete GraphsNuMu[i][j];
        }
        if (GraphsNuMu[i] != NULL) delete[] GraphsNuMu[i];
      }

      // Delete the nue(bar) graphs
      for (int i = 0; i < kNuMuBar+1; ++i) {
        for (int j = 0; j < nEnuNuE; ++j) {
          if (GraphsNuE[i][j] != NULL) delete GraphsNuE[i][j];
        }
        if (GraphsNuE[i] != NULL) delete[] GraphsNuE[i];
      }
    }

    // Calculate a weight using a fixed Enu, Q2 and neutrino type
    double CalcWeight(double Enu, double Q2);

    // Old method, mostly for documentation. Produces some discontinuities in the interpolation
    double CalcWeightOld(double Enu, double Q2);

    // Setters
    void SetLeptonMass(double input) { leptonmass = input; }; // In GeV

    // Set the neutrino type, and the lepton mass
    void SetNuType(NuType nu) {
      nutype = nu;
      if (nu == kNuMu || nu == kNuMuBar) {
        leptonmass = 0.10566;
        nEnu = nEnuNuMu;
        EnuRange = EnuRangeNuMu;
      } else if (nu == kNuE || nu == kNuEBar) {
        leptonmass = 0.511e-3;
        nEnu = nEnuNuE;
        EnuRange = EnuRangeNuE;
      } else {
        std::cerr << "Calculator does not support neutrino type " << nu << std::endl;
        std::cerr << "0 = numu, 1 = numubar, 2 = nue, 3 = nuebar" << std::endl;
        throw;
      }
    };

    void SetLinearInterp() { drawcmd = ""; };
    void SetSplineInterp() { drawcmd = "S"; };

    // Getters
    std::string GetSplineInterp() { return drawcmd; };
    double* GetEnuRange() { return EnuRange; };
    double* GetEnuRangeNuMu() { return EnuRangeNuMu; };
    double* GetEnuRangeNuE() { return EnuRangeNuE; };
    int GetNEnu() { return nEnu; };
    int GetNEnuNuMu() { return nEnuNuMu; };
    int GetNEnuNuE() { return nEnuNuE; };

    // Get the maximum Q2 for a given Enu and muon mass
    // Get the max Q2 for a given Enu to check interpolaton
    // CCQE only
    double GetQ2max(double Enu) {
      // Nucleon mass
      const double Mn = 0.93956542052;
      const double Mp = 0.93827208816;
      const double M = (Mn+Mp)/2.;
      const double M2 = M*M;
      const double leptonmass2 = leptonmass*leptonmass;
      const double Enu2 = Enu*Enu;

      double val = -(M+Enu)*leptonmass2+2*M*Enu2+Enu*sqrt(pow(2*M*Enu-leptonmass2, 2)-4*leptonmass2*M2);
      val /= (M+2*Enu);
      return val;
    }

  private:
    // Input TGraphs, first index is numu or numubar, second index is the fixed enu
    TGraph ***GraphsNuMu;
    TGraph ***GraphsNuE;
    // Have different number of pre-calculated points for numu and nue, so keep this pointing to the relevant one
    TGraph ***Graphs;

    // Enu Range that the inputs come in
    double EnuRangeNuMu[20];
    double EnuRangeNuE[9];
    double *EnuRange;

    // Just a hard-coded check to see if all the graphs are there
    int nEnuNuMu;
    int nEnuNuE;
    int nEnu;

    // Lepton mass in GeV
    double leptonmass;
    NuType nutype;

    // Interpolation method, spline or linear
    std::string drawcmd;

    bool CheckSetup(double &Enu, double &Q2);

    bool Checked;

    // Q2 tolerance, used in old calculation for checking Q2 range
    double Q2tol;
};

#endif