PWGLF/FORWARD/analysis2/corrs/PDFLandau.C

   1 #include <TF1.h>
   2 #include <TH1.h>
   3 #include <TMath.h>
   4 #include <cmath>
   5
   6 struct PDF
   7 {
   8
   9   static double cheb(double v, double* p, double* q, int n) {
  10     double num = 0;
  11     double den = 0;
  12     for (int i = n-1; i >= 0; i--) {
  13       num = (p[i] + num)*v;
  14       den = (q[i] + den)*v;
  15     }
  16     return num/den;
  17   }
  18   static double f1(double v) {
  19     static double a1[3] = {0.04166666667,-0.01996527778, 0.02709538966};
  20
  21     double u = std::exp(v+1);
  22     if (u < 1e-10) return 0;
  23
  24     double ue = std::exp(-1/u);
  25     double us = std::sqrt(u);
  26     return 0.3989422803*(ue/us)*(1+(a1[0]+(a1[1]+a1[2]*u)*u)*u);
  27   }
  28   static double f2(double v) {
  29     static double p1[] = {0.4259894875,  -0.1249762550,
  30                           0.03984243700, -0.006298287635,
  31                           0.001511162253};
  32     static double q1[] = {1.0,           -0.3388260629,
  33                           0.09594393323, -0.01608042283,
  34                           0.003778942063};
  35
  36     double u = std::exp(-v-1);
  37     return std::exp(-u)*std::sqrt(u)*cheb(v,p1,q1,5);
  38   }
  39   static double f3(double v) {
  40     static double p2[] = {0.1788541609,   0.1173957403,
  41                           0.01488850518, -0.001394989411,
  42                           0.0001283617211};
  43     static double q2[] = {1.0,            0.7428795082,
  44                           0.3153932961,   0.06694219548,
  45                           0.008790609714};
  46     return cheb(v,p2,q2,5);
  47   }
  48   static double f4(double v) {
  49     static double p3[] = {0.1788544503,   0.09359161662,
  50                           0.006325387654, 0.00006611667319,
  51                           -0.000002031049101};
  52     static double q3[] = {1.0,            0.6097809921,
  53                           0.2560616665,   0.04746722384,
  54                           0.006957301675};
  55     return cheb(v,p3,q3,5);
  56   }
  57   static double f5(double v) {
  58     static double p4[] = {0.9874054407, 118.6723273,
  59                           849.2794360, -743.7792444,
  60                           427.0262186};
  61     static double q4[] = {1.0,           106.8615961,
  62                           337.6496214,   2016.712389,
  63                           1597.063511};
  64     double u   = 1/v;
  65     return u*u*cheb(u,p4,q4,5);
  66   }
  67   static double f6(double v) {
  68     static double p5[] = {1.003675074,   167.5702434,
  69                            4789.711289, 21217.86767,
  70                            -22324.94910};
  71     static double q5[] = {1.0,          156.9424537,
  72                            3745.310488,  9834.698876,
  73                            66924.28357};
  74     double u   = 1/v;
  75     return u*u*cheb(u,p5,q5,5);
  76   }
  77   static double f7(double v) {
  78     static double p6[] = {1.000827619,  664.9143136,
  79                           62972.92665,  475554.6998,
  80                           -5743609.109};
  81     static double q6[] = {1.0,          651.4101098,
  82                           56974.73333,  165917.4725,
  83                           -2815759.939};
  84     double u   = 1/v;
  85     return u*u*cheb(u,p6,q6,5);
  86   }
  87   static double f8(double v) {
  88     static double a2[2] = {-1.845568670,-4.284640743};
  89     double u = 1/(v-v*std::log(v)/(v+1));
  90     return u*u*(1+(a2[0]+a2[1]*u)*u);
  91   }
  92
  93   static Double_t wrap(Double_t* px, Double_t* pp) {
  94     Int_t    i = pp[0];
  95     Double_t v = px[0];
  96     switch (i) {
  97     case 1: return f1(v);
  98     case 2: return f2(v);
  99     case 3: return f3(v);
 100     case 4: return f4(v);
 101     case 5: return f5(v);
 102     case 6: return f6(v);
 103     case 7: return f7(v);
 104     case 8: return f8(v);
 105     }
 106     return 0;
 107   }
 108
 109   static TF1* getComp(Int_t i) {
 110     Double_t l = -10;
 111     Double_t h = -5.5;
 112     switch (i) {
 113     case 2: l =   -5.5; h =   -1; break;
 114     case 3: l =   -1;   h =    1; break;
 115     case 4: l =    1;   h =    5; break;
 116     case 5: l =    5;   h =   12; break;
 117     case 6: l =   12;   h =   50; break;
 118     case 7: l =   50;   h =  300; break;
 119     case 8: l =  300;   h = 1000; break;
 120     default: i = 1; break;
 121     }
 122
 123     l *= (l < 0 ? 1.5 : 0.4);
 124     h *= (h < 0 ? 0.5 : 1.5);
 125
 126     TF1* f = new TF1(Form("f%d", i), &PDF::wrap, l, h, 1);
 127     f->SetParameter(0, i);
 128     f->SetLineColor(i);
 129     return f;
 130   }
 131   static void draw() {
 132     TH1* h = new TH1F("h","h", 1000, -10, 1000);
 133     h->Draw();
 134
 135     TF1* l = new TF1("l", "landau", -10, 1000);
 136     l->SetParameters(1,0,1);
 137     l->SetLineStyle(2);
 138     l->Draw("same");
 139
 140     Double_t max = 0;
 141     for (Int_t i = 1; i <= 8; i++) {
 142       TF1* f = getComp(i);
 143       f->Draw("same");
 144       max = TMath::Max(f->GetMaximum(),max);
 145     }
 146     h->SetMaximum(1.1*max);
 147   }
 148 };
 149
 150
 151