PWGLF/SPECTRA/IdentifiedHighPt/macros/my_functions.C

   1 #include <TMath.h>
   2 #include <TH2.h>
   3 #include <TProfile.h>
   4
   5 #include <iostream>
   6
   7 using namespace std;
   8
   9 const Double_t eMass  = 0.000511;               //electron mass
  10 const Double_t piMass = 0.13957;               //pion mass
  11 const Double_t kMass  = 0.493676999999999977;  //kaon mass
  12 const Double_t pMass  = 0.938271999999999995;  //proton mass
  13
  14
  15 Double_t fitf3G(Double_t* xx, Double_t* par);
  16 Double_t FitFunc(Double_t* x, Double_t *par);
  17 Double_t SigmaFunc(Double_t* x, Double_t *par);
  18
  19 TH2D* hDeDxVsP = 0;
  20 TProfile* hMeanP = 0;
  21 Double_t fixMIP     = 50.0;
  22 Double_t fixPlateau = 75.0;
  23
  24
  25 class DeDxFitInfo : public TObject
  26 {
  27  public:
  28
  29   DeDxFitInfo();
  30   void Print(Option_t* option="") const;
  31
  32   Double_t MIP;
  33   Double_t plateau;
  34
  35   Int_t optionDeDx;
  36   Int_t nDeDxPar;
  37   Double_t parDeDx[8];
  38
  39   Int_t optionSigma;
  40   Int_t nSigmaPar;
  41   Double_t parSigma[8];
  42
  43   TString calibFileName;
  44
  45   ClassDef(DeDxFitInfo, 1);    // Help class
  46 };
  47
  48 //_____________________________________________________________________________
  49 ClassImp(DeDxFitInfo)
  50
  51 DeDxFitInfo::DeDxFitInfo():
  52 TObject(),
  53   MIP(0),
  54   plateau(0),
  55   optionDeDx(-1),
  56   nDeDxPar(-1),
  57   optionSigma(-1),
  58   nSigmaPar(-1),
  59   calibFileName("")
  60 {
  61   // default constructor
  62   for(Int_t i = 0; i < 8; i++) {
  63     parDeDx[i]  = 0;
  64     parSigma[i] = 0;
  65   }
  66 }
  67
  68 //_________________________________________________________
  69 void DeDxFitInfo::Print(Option_t* option) const
  70 {
  71   if(option)
  72     cout << "Option: " << option << endl;
  73
  74   cout << ClassName() << " : " << GetName() << endl
  75        << "MIP: " << MIP << endl
  76        << "Plateau: " << plateau << endl
  77        << "OptionDeDx: " << optionDeDx << endl
  78        << "nDeDxPar: " << nDeDxPar << endl;
  79   for(Int_t i = 0; i < nDeDxPar; i++) {
  80
  81     cout << "parDeDx[" << i << "] = " << parDeDx[i] << endl;
  82   }
  83   cout << "OptionSigma: " << optionSigma << endl
  84        << "nSigmaPar: " << nSigmaPar << endl;
  85   for(Int_t i = 0; i < nSigmaPar; i++) {
  86
  87     cout << "parSigma[" << i << "] = " << parSigma[i] << endl;
  88   }
  89
  90   if(calibFileName.IsNull()) {
  91     cout << "No eta calibration file." << endl;
  92   } else {
  93     cout << "Eta calibration file: " << calibFileName.Data() << endl;
  94   }
  95 }
  96
  97 //______________________________________________________________________________
  98 Double_t fitf3G(Double_t* xx, Double_t* par)
  99 {
 100   //
 101   // Could speed up fit by forcing it to use <p>. In that way the parameters
 102   // could be amde statis cand only changed when going to a new p bin
 103   //
 104   Double_t p = xx[0];
 105   Double_t dedx = xx[1];
 106
 107   const Int_t bin = hDeDxVsP->GetXaxis()->FindBin(p);
 108
 109   if(hMeanP) {
 110     //    cout << "p before: " << p;
 111     p = hMeanP->GetBinContent(bin);
 112     //    cout << ", p after: " << p << endl;
 113   }
 114   const Int_t binStart = Int_t(par[0]);
 115   const Int_t binStop  = Int_t(par[1]);
 116
 117   if(bin<binStart || bin>binStop) {
 118
 119     cout << "Error: bin " << bin << " not inside inteval [" << binStart
 120          << "; " << binStop << "]" << endl;
 121     return 0;
 122   }
 123
 124   const Int_t nParDeDx = Int_t(par[2]);
 125
 126   Double_t* parDeDx = &par[3];
 127
 128   Int_t offset = 4 + nParDeDx; // binStart + binStop + nParDeDx + optionDedx + nParDeDx parameters
 129
 130   const Int_t nParSigma = Int_t(par[offset]);
 131   offset += 1; // nParSigma
 132
 133   Double_t* parSigma = &par[offset];
 134   offset += 1 + nParSigma; // optionSigma + nParSigma parameters
 135
 136   Double_t piMean = FitFunc(&p, parDeDx);
 137   Double_t pKeff = p*piMass/kMass; // corresponding p of a pion with same dE/dx
 138   Double_t kMean  = FitFunc(&pKeff, parDeDx);
 139   Double_t pPeff = p*piMass/pMass; // corresponding p of a pion with same dE/dx
 140   Double_t pMean  = FitFunc(&pPeff, parDeDx);
 141
 142   const Double_t piSigma = SigmaFunc(&piMean, parSigma);
 143   const Double_t kSigma  = SigmaFunc(&kMean,  parSigma);
 144   const Double_t pSigma  = SigmaFunc(&pMean,  parSigma);
 145
 146
 147   const Int_t j = bin - binStart;
 148   const Double_t piYield   = par[j * 3 + offset + 0];
 149   const Double_t kYield    = par[j * 3 + offset + 1];
 150   const Double_t pYield    = par[j * 3 + offset + 2];
 151
 152   return piYield* TMath::Gaus(dedx, piMean, piSigma, kTRUE)
 153     +    kYield * TMath::Gaus(dedx, kMean,  kSigma,  kTRUE)
 154     +    pYield * TMath::Gaus(dedx, pMean,  pSigma,  kTRUE);
 155 }
 156
 157
 158 //______________________________________________________________________________
 159 Double_t FitFunc(Double_t* x, Double_t *par)
 160 {
 161   static const Double_t bgMIP    = 0.5/piMass;
 162   static const Double_t beta2MIP = bgMIP*bgMIP / (1.0+bgMIP*bgMIP);
 163   //  static const Double_t betapowMIP = TMath;
 164   static const Double_t logMIP   = TMath::Log(1+bgMIP);
 165
 166
 167
 168   Int_t option = TMath::Nint(par[0]);
 169   Int_t specie = option;
 170   option = option%10;
 171   specie -= option;
 172   specie /= 10;
 173
 174
 175   Double_t bg = 0;
 176   switch (specie) {
 177
 178   case 0: // pion
 179     bg = x[0]/piMass;
 180     break;
 181   case 1: // kaon
 182     bg = x[0]/kMass;
 183     break;
 184   case 2: // proton
 185     bg = x[0]/pMass;
 186     break;
 187   case 3: // electron
 188     bg = x[0]/eMass;
 189     break;
 190   case 4: // just use bg
 191     bg = x[0];
 192     break;
 193   default:
 194     cout << "Error in FitFunc: specie " << specie << " not supported!!!!!" << endl;
 195     return 0;
 196     break;
 197   }
 198
 199
 200   if(bg > 10000.0)
 201     bg = 10000.0;
 202
 203   const Double_t beta2 = bg*bg / (1.0+bg*bg);
 204
 205   switch (option) {
 206
 207   case 1: // standard parametrisation
 208     {
 209       /*
 210         c0/beta^2 + c1 * log (1+x)
 211        */
 212       const Double_t c0 = par[1];
 213       const Double_t c1 = par[2];
 214
 215       const Double_t value = c0/beta2 + c1*TMath::Log(1+bg);
 216       return value;
 217     }
 218     break;
 219   case 2: // fix the dE/dx to 50 at 0.5 GeV/c
 220     {
 221       const Double_t c1 = par[1];
 222       const Double_t c0 = (fixMIP-par[1]*logMIP) * beta2MIP;
 223
 224       const Double_t value = c0/beta2 + c1*TMath::Log(1+bg);
 225       return value;
 226     }
 227     break;
 228   case 3: // fix the dE/dx to 50 at 0.5 GeV/c and the plateau to 75
 229     {
 230       /*
 231         a/beta^2 + b/c*log( (1+x)^c / (1 + d*(1+x)^c) )
 232
 233         Assymptotic behavior:
 234
 235         1) Small bg (and d small so that d*(1+x)^c << 1)
 236
 237         a/beta^2 + b * log (1+x)
 238
 239         So this is the same beavior as the standard expression.
 240
 241         2) Large bg where d*(1+x)^c >> 1
 242         a - b/c*log(d) = plateau
 243         -> d = exp(c*(a-plateau)/b)
 244
 245        */
 246       const Double_t b = par[1];
 247       const Double_t a = (fixMIP-par[1]*logMIP) * beta2MIP;
 248       const Double_t c = par[2];
 249       const Double_t d = TMath::Exp(c*(a-fixPlateau)/b);
 250
 251       //      cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
 252
 253       const Double_t powbg = TMath::Power(1.0+bg, c);
 254
 255       const Double_t value = a/beta2 + b/c*TMath::Log(powbg/(1.0 + d*powbg));
 256       return value;
 257     }
 258     break;
 259   case 4: // fix the dE/dx to 50 at 0.5 GeV/c and the plateau to 75
 260     {
 261       /*
 262         a/beta^2 + b/c*log( (1+x)^c / (1 + d*(1+x)^c) )
 263
 264         Assymptotic behavior:
 265
 266         1) Small bg (and d small so that d*(1+x)^c << 1)
 267
 268         a/beta^2 + b * log (1+x)
 269
 270         So this is the same beavior as the standard expression.
 271
 272         2) Large bg where d*(1+x)^c >> 1
 273         a - b/c*log(d) = plateau
 274         -> d = exp(c*(a-plateau)/b)
 275
 276        */
 277       const Double_t a = par[1];
 278       const Double_t b = par[2];
 279       const Double_t c = par[3];
 280       const Double_t d = TMath::Exp(c*(a-fixPlateau)/b);
 281
 282       //      cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
 283
 284       const Double_t powbg = TMath::Power(1.0+bg, c);
 285
 286       const Double_t value = a/beta2 + b/c*TMath::Log(powbg/(1.0 + d*powbg));
 287       return value;
 288     }
 289     break;
 290   case 5: // fix the dE/dx to 50 at 0.5 GeV/c and the plateau to 75
 291     {
 292       /*
 293         a/beta^2 + b/c*log( (1+x)^c / (1 + d*(1+x)^c) )
 294
 295         Assymptotic behavior:
 296
 297         1) Small bg (and d small so that d*(1+x)^c << 1)
 298
 299         a/beta^2 + b * log (1+x)
 300
 301         So this is the same beavior as the standard expression.
 302
 303         2) Large bg where d*(1+x)^c >> 1
 304         a - b/c*log(d) = plateau
 305         -> d = exp(c*(a-plateau)/b)
 306
 307        */
 308       const Double_t a = par[1];
 309       const Double_t b = par[2];
 310       const Double_t c = par[3];
 311       const Double_t d = TMath::Exp(c*(a-fixPlateau)/b);
 312       const Double_t e = par[4];
 313
 314       //      cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
 315
 316       const Double_t powbg = TMath::Power(1.0+bg, c);
 317
 318       const Double_t value = a/TMath::Power(beta2,e) + b/c*TMath::Log(powbg/(1.0 + d*powbg));
 319       return value;
 320     }
 321     break;
 322   case 6:
 323     {
 324       /*
 325         a/beta^(e/2) + b/c*log( (1+x)^c / (1 + d*(1+x)^c) )
 326
 327         Assymptotic behavior:
 328
 329         1) Small bg (and d small so that d*(1+x)^c << 1)
 330
 331         a/beta^(e/2) + b * log (1+x)
 332
 333         So this is the same beavior as the standard expression.
 334
 335         2) Large bg where d*(1+x)^c >> 1
 336         a - b/c*log(d) = plateau
 337         -> d = exp(c*(a-plateau)/b)
 338
 339         In this version we have 2 plateaus!!!!
 340         Plateau 1 = electrons!
 341
 342        */
 343       if(specie==3)
 344         return fixPlateau;
 345
 346       const Double_t a = par[1];
 347       const Double_t b = par[2];
 348       const Double_t c = par[3];
 349       const Double_t d = TMath::Exp(c*(a-par[5])/b);
 350       const Double_t e = par[4];
 351
 352       //      cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
 353
 354       const Double_t powbg = TMath::Power(1.0+bg, c);
 355
 356       const Double_t value = a/TMath::Power(beta2,e) + b/c*TMath::Log(powbg/(1.0 + d*powbg));
 357       return value;
 358     }
 359     break;
 360   case 7:
 361     {
 362       /*
 363         a/beta^(d/2) - b*log( c + 1.0/(1.0+x) )
 364
 365         Assymptotic behavior:
 366
 367         1) Small bg (and d small so that d*(1+x)^c << 1)
 368
 369         a/beta^(d/2) - b * log (1+x)
 370
 371         So this is the same beavior as the standard expression.
 372
 373         2) Large bg where c << 1
 374         -b*log(c) = plateau-a
 375         -> c = exp((a-plateau)/b)
 376
 377        */
 378       const Double_t a = par[1];
 379       const Double_t b = par[2];
 380       const Double_t c = TMath::Exp((a-fixPlateau)/b);
 381       const Double_t d = par[3];
 382
 383       //      cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
 384
 385       const Double_t value = a/TMath::Power(beta2,d) - b*TMath::Log(c + 1.0/(1.0+bg));
 386       return value;
 387     }
 388   case 8:
 389     {
 390       /*
 391         a/beta^(d/2) - b*log( c + 1.0/(1.0+x^e) )
 392
 393         Assymptotic behavior:
 394
 395         1) Small bg (and d small so that d*(1+x)^c << 1)
 396
 397         a/beta^(d/2) - b * log (1+x^e)
 398
 399         So this is the same beavior as the standard expression.
 400
 401         2) Large bg where c << 1
 402         -b*log(c) = plateau-a
 403         -> c = exp((a-plateau)/b)
 404
 405        */
 406       const Double_t a = par[1];
 407       const Double_t b = par[2];
 408       const Double_t c = TMath::Exp((a-fixPlateau)/b);
 409       const Double_t d = par[3];
 410       const Double_t e = par[4];
 411
 412       //      cout << bg << ": " << a << ", " << b << ", " << c << ", " << d << endl;
 413
 414       const Double_t value = a/TMath::Power(beta2,d) - b*TMath::Log(c + 1.0/(1.0+bg) + e/(1.0+TMath::Power(bg, 2)));
 415       return value;
 416     }
 417     break;
 418     // case 3: // fix the dE/dx to 50 at 0.5 GeV/c + powerlaw
 419
 420     //   static const bgMIP    = 0.5/piMass;
 421     //   static const beta2MIP = bgMIP*bgMIP / (1.0+bgMIP*bgMIP);
 422     //   static const logMIP   = TMath::Log(1+bgMIP);
 423
 424     //   const Double_t c1 = par[1];
 425     //   const Double_t c2 = par[2]; // expect it to be 0.75 from Bichsel (beta**-1.5 instead of -2)
 426     //   const Double_t c0 = (50.0-par[1]*logMIP) * beta2MIP;
 427
 428     //   const Double_t value = TMathh::Power(c0/beta2, c2) + c1*TMath::Log(1+bg);
 429     //   return value;
 430
 431   default:
 432     break;
 433   }
 434
 435   cout << "Error in FitFunc: option " << option << " not supported!!!!!" << endl;
 436   return 0;
 437 }
 438
 439 //______________________________________________________________________________
 440 Double_t SigmaFunc(Double_t* x, Double_t *par)
 441 {
 442   Int_t option = Int_t(par[0]);
 443
 444   switch (option) {
 445   case 1: // fixed sigma
 446     return par[1];
 447     break;
 448   case 2: // relative sigma
 449     return par[1]*x[0];
 450     break;
 451   case 3: // relative sigma + extrapolation
 452     return (par[1] + (x[0]-fixMIP)*par[2])*x[0];
 453     break;
 454   case 4: // relative sigma with dE/dx to some power close to 1
 455     return par[1]*TMath::Power(x[0], par[2]);
 456     break;
 457   case 5: // relative sigma with dE/dx to some power close to 1
 458     return TMath::Sqrt(par[1]*par[1]*x[0]*x[0] + par[2]*par[2]);
 459     break;
 460   case 6: // relative sigma with dE/dx to some power close to 1
 461     return par[1]*x[0]*TMath::Power(x[0]/50.0, par[2]);
 462     break;
 463   case 7: // 1/x^some power + constant
 464     return x[0]*(par[2]*TMath::Power(x[0], -par[3]) + par[1]);
 465     break;
 466   case 8: // for fitting relative sigma
 467     return par[2]*TMath::Power(x[0], -par[3]) + par[1];
 468     break;
 469   case 9: // for fitting relative sigma
 470     return par[1]+par[2]*x[0]+par[3]*x[0]*x[0];
 471     break;
 472   case 10: // for fitting relative sigma
 473     return par[1]+par[2]*x[0]+par[3]*x[0]*x[0]-par[4]/(x[0]*x[0]*x[0])-par[5]/(x[0]*x[0]);
 474     break;
 475   case 11: // for fitting relative sigma
 476     return x[0]*(par[1]+par[2]*x[0]+par[3]*x[0]*x[0]-par[4]/(x[0]*x[0]*x[0])-par[5]/(x[0]*x[0]));
 477     break;
 478   case 12: // for fitting relative sigma
 479     return par[1]+par[2]*x[0]+par[3]*x[0]*x[0];
 480     break;
 481   case 13: // for fitting relative sigma
 482     return x[0]*(par[1]+par[2]*x[0]+par[3]*x[0]*x[0]);
 483     break;
 484   case 14: // for fitting relative sigma
 485     return par[1]+par[2]*x[0];
 486     break;
 487   case 15: // for fitting relative sigma
 488     return x[0]*(par[1]+par[2]*x[0]);
 489     break;
 490
 491
 492
 493
 494
 495
 496   default:
 497     break;
 498   }
 499
 500   cout << "Error in SigmaFunc: option " << option << " not supported!!!!!" << endl;
 501   return 0;
 502 }