TRD/AliTRDv1.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Transition Radiation Detector version 1 -- slow simulator                //
  21 //                                                                           //
  22 //Begin_Html
  23 /*
  24 <img src="picts/AliTRDfullClass.gif">
  25 */
  26 //End_Html
  27 //                                                                           //
  28 //                                                                           //
  29 ///////////////////////////////////////////////////////////////////////////////
  30
  31 #include <stdlib.h>
  32
  33 #include <TF1.h>
  34 #include <TLorentzVector.h>
  35 #include <TMath.h>
  36 #include <TRandom.h>
  37 #include <TVector.h>
  38 #include <TVirtualMC.h>
  39
  40 #include "AliConst.h"
  41 #include "AliRun.h"
  42 #include "AliTRDgeometry.h"
  43 #include "AliTRDhit.h"
  44 #include "AliTRDsim.h"
  45 #include "AliTRDv1.h"
  46 #include "AliMC.h"
  47
  48 ClassImp(AliTRDv1)
  49
  50 //_____________________________________________________________________________
  51 AliTRDv1::AliTRDv1():AliTRD()
  52 {
  53   //
  54   // Default constructor
  55   //
  56
  57   fSensSelect      =  0;
  58   fSensPlane       = -1;
  59   fSensChamber     = -1;
  60   fSensSector      = -1;
  61   fSensSectorRange =  0;
  62
  63   fDeltaE          = NULL;
  64   fTR              = NULL;
  65
  66 }
  67
  68 //_____________________________________________________________________________
  69 AliTRDv1::AliTRDv1(const char *name, const char *title)
  70          :AliTRD(name, title)
  71 {
  72   //
  73   // Standard constructor for Transition Radiation Detector version 1
  74   //
  75
  76   fSensSelect      =  0;
  77   fSensPlane       = -1;
  78   fSensChamber     = -1;
  79   fSensSector      = -1;
  80   fSensSectorRange =  0;
  81
  82   fDeltaE          = NULL;
  83   fTR              = NULL;
  84
  85   SetBufferSize(128000);
  86
  87 }
  88
  89 //_____________________________________________________________________________
  90 AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
  91 {
  92   //
  93   // Copy constructor
  94   //
  95
  96   ((AliTRDv1 &) trd).Copy(*this);
  97
  98 }
  99
 100 //_____________________________________________________________________________
 101 AliTRDv1::~AliTRDv1()
 102 {
 103   //
 104   // AliTRDv1 destructor
 105   //
 106
 107   if (fDeltaE) delete fDeltaE;
 108   if (fTR)     delete fTR;
 109
 110 }
 111
 112 //_____________________________________________________________________________
 113 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
 114 {
 115   //
 116   // Assignment operator
 117   //
 118
 119   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
 120   return *this;
 121
 122 }
 123
 124 //_____________________________________________________________________________
 125 void AliTRDv1::Copy(TObject &trd)
 126 {
 127   //
 128   // Copy function
 129   //
 130
 131   ((AliTRDv1 &) trd).fSensSelect      = fSensSelect;
 132   ((AliTRDv1 &) trd).fSensPlane       = fSensPlane;
 133   ((AliTRDv1 &) trd).fSensChamber     = fSensChamber;
 134   ((AliTRDv1 &) trd).fSensSector      = fSensSector;
 135   ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
 136
 137   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
 138   fTR->Copy(*((AliTRDv1 &) trd).fTR);
 139
 140 }
 141
 142 //_____________________________________________________________________________
 143 void AliTRDv1::CreateGeometry()
 144 {
 145   //
 146   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 147   // This version covers the full azimuth.
 148   //
 149
 150   // Check that FRAME is there otherwise we have no place where to put the TRD
 151   AliModule* frame = gAlice->GetModule("FRAME");
 152   if (!frame) return;
 153
 154   // Define the chambers
 155   AliTRD::CreateGeometry();
 156
 157 }
 158
 159 //_____________________________________________________________________________
 160 void AliTRDv1::CreateMaterials()
 161 {
 162   //
 163   // Create materials for the Transition Radiation Detector version 1
 164   //
 165
 166   AliTRD::CreateMaterials();
 167
 168 }
 169
 170 //_____________________________________________________________________________
 171 void AliTRDv1::CreateTRhit(Int_t det)
 172 {
 173   //
 174   // Creates an electron cluster from a TR photon.
 175   // The photon is assumed to be created a the end of the radiator. The
 176   // distance after which it deposits its energy takes into account the
 177   // absorbtion of the entrance window and of the gas mixture in drift
 178   // volume.
 179   //
 180
 181   // PDG code electron
 182   const Int_t   kPdgElectron = 11;
 183
 184   // Ionization energy
 185   const Float_t kWion        = 22.04;
 186
 187   // Maximum number of TR photons per track
 188   const Int_t   kNTR         = 50;
 189
 190   TLorentzVector mom, pos;
 191
 192   // Create TR at the entrance of the chamber
 193   if (gMC->IsTrackEntering()) {
 194
 195     // Create TR only for electrons
 196     Int_t iPdg = gMC->TrackPid();
 197     if (TMath::Abs(iPdg) != kPdgElectron) return;
 198
 199     Float_t eTR[kNTR];
 200     Int_t   nTR;
 201
 202     // Create TR photons
 203     gMC->TrackMomentum(mom);
 204     Float_t pTot = mom.Rho();
 205     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
 206     if (nTR > kNTR) {
 207       printf("AliTRDv1::CreateTRhit -- ");
 208       printf("Boundary error: nTR = %d, kNTR = %d\n",nTR,kNTR);
 209       exit(1);
 210     }
 211
 212     // Loop through the TR photons
 213     for (Int_t iTR = 0; iTR < nTR; iTR++) {
 214
 215       Float_t energyMeV = eTR[iTR] * 0.001;
 216       Float_t energyeV  = eTR[iTR] * 1000.0;
 217       Float_t absLength = 0;
 218       Float_t sigma     = 0;
 219
 220       // Take the absorbtion in the entrance window into account
 221       Double_t muMy = fTR->GetMuMy(energyMeV);
 222       sigma = muMy * fFoilDensity;
 223       if (sigma > 0.0) {
 224         absLength = gRandom->Exp(1.0/sigma);
 225         if (absLength < AliTRDgeometry::MyThick()) continue;
 226       }
 227       else {
 228         continue;
 229       }
 230
 231       // The absorbtion cross sections in the drift gas
 232       if (fGasMix == 1) {
 233         // Gas-mixture (Xe/CO2)
 234         Double_t muXe = fTR->GetMuXe(energyMeV);
 235         Double_t muCO = fTR->GetMuCO(energyMeV);
 236         sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
 237       }
 238       else {
 239         // Gas-mixture (Xe/Isobutane)
 240         Double_t muXe = fTR->GetMuXe(energyMeV);
 241         Double_t muBu = fTR->GetMuBu(energyMeV);
 242         sigma = (0.97 * muXe + 0.03 * muBu) * fGasDensity * fTR->GetTemp();
 243       }
 244
 245       // The distance after which the energy of the TR photon
 246       // is deposited.
 247       if (sigma > 0.0) {
 248         absLength = gRandom->Exp(1.0/sigma);
 249         if (absLength > AliTRDgeometry::DrThick()) continue;
 250       }
 251       else {
 252         continue;
 253       }
 254
 255       // The position of the absorbtion
 256       Float_t posHit[3];
 257       gMC->TrackPosition(pos);
 258       posHit[0] = pos[0] + mom[0] / pTot * absLength;
 259       posHit[1] = pos[1] + mom[1] / pTot * absLength;
 260       posHit[2] = pos[2] + mom[2] / pTot * absLength;
 261
 262       // Create the charge
 263       Int_t q = ((Int_t) (energyeV / kWion));
 264
 265       // Add the hit to the array. TR photon hits are marked
 266       // by negative charge
 267       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
 268
 269     }
 270
 271   }
 272
 273 }
 274
 275 //_____________________________________________________________________________
 276 void AliTRDv1::Init()
 277 {
 278   //
 279   // Initialise Transition Radiation Detector after geometry has been built.
 280   //
 281
 282   AliTRD::Init();
 283
 284   if(fDebug) printf("%s: Slow simulator\n",ClassName());
 285   if (fSensSelect) {
 286     if (fSensPlane   >= 0)
 287       printf("          Only plane %d is sensitive\n",fSensPlane);
 288     if (fSensChamber >= 0)
 289       printf("          Only chamber %d is sensitive\n",fSensChamber);
 290     if (fSensSector  >= 0) {
 291       Int_t sens1  = fSensSector;
 292       Int_t sens2  = fSensSector + fSensSectorRange;
 293             sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 294                    * AliTRDgeometry::Nsect();
 295       printf("          Only sectors %d - %d are sensitive\n",sens1,sens2-1);
 296     }
 297   }
 298   if (fTR)
 299     printf("%s: TR simulation on\n",ClassName());
 300   else
 301     printf("%s: TR simulation off\n",ClassName());
 302   printf("\n");
 303
 304   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 305   const Float_t kPoti = 12.1;
 306   // Maximum energy (50 keV);
 307   const Float_t kEend = 50000.0;
 308   // Ermilova distribution for the delta-ray spectrum
 309   Float_t poti = TMath::Log(kPoti);
 310   Float_t eEnd = TMath::Log(kEend);
 311   fDeltaE = new TF1("deltae",Ermilova,poti,eEnd,0);
 312
 313   if(fDebug) {
 314     printf("%s: ",ClassName());
 315     for (Int_t i = 0; i < 80; i++) printf("*");
 316     printf("\n");
 317   }
 318
 319 }
 320
 321 //_____________________________________________________________________________
 322 AliTRDsim *AliTRDv1::CreateTR()
 323 {
 324   //
 325   // Enables the simulation of TR
 326   //
 327
 328   fTR = new AliTRDsim();
 329   return fTR;
 330
 331 }
 332
 333 //_____________________________________________________________________________
 334 void AliTRDv1::SetSensPlane(Int_t iplane)
 335 {
 336   //
 337   // Defines the hit-sensitive plane (0-5)
 338   //
 339
 340   if ((iplane < 0) || (iplane > 5)) {
 341     printf("Wrong input value: %d\n",iplane);
 342     printf("Use standard setting\n");
 343     fSensPlane  = -1;
 344     fSensSelect =  0;
 345     return;
 346   }
 347
 348   fSensSelect = 1;
 349   fSensPlane  = iplane;
 350
 351 }
 352
 353 //_____________________________________________________________________________
 354 void AliTRDv1::SetSensChamber(Int_t ichamber)
 355 {
 356   //
 357   // Defines the hit-sensitive chamber (0-4)
 358   //
 359
 360   if ((ichamber < 0) || (ichamber > 4)) {
 361     printf("Wrong input value: %d\n",ichamber);
 362     printf("Use standard setting\n");
 363     fSensChamber = -1;
 364     fSensSelect  =  0;
 365     return;
 366   }
 367
 368   fSensSelect  = 1;
 369   fSensChamber = ichamber;
 370
 371 }
 372
 373 //_____________________________________________________________________________
 374 void AliTRDv1::SetSensSector(Int_t isector)
 375 {
 376   //
 377   // Defines the hit-sensitive sector (0-17)
 378   //
 379
 380   SetSensSector(isector,1);
 381
 382 }
 383
 384 //_____________________________________________________________________________
 385 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
 386 {
 387   //
 388   // Defines a range of hit-sensitive sectors. The range is defined by
 389   // <isector> (0-17) as the starting point and <nsector> as the number
 390   // of sectors to be included.
 391   //
 392
 393   if ((isector < 0) || (isector > 17)) {
 394     printf("Wrong input value <isector>: %d\n",isector);
 395     printf("Use standard setting\n");
 396     fSensSector      = -1;
 397     fSensSectorRange =  0;
 398     fSensSelect      =  0;
 399     return;
 400   }
 401
 402   if ((nsector < 1) || (nsector > 18)) {
 403     printf("Wrong input value <nsector>: %d\n",nsector);
 404     printf("Use standard setting\n");
 405     fSensSector      = -1;
 406     fSensSectorRange =  0;
 407     fSensSelect      =  0;
 408     return;
 409   }
 410
 411   fSensSelect      = 1;
 412   fSensSector      = isector;
 413   fSensSectorRange = nsector;
 414
 415 }
 416
 417 //_____________________________________________________________________________
 418 void AliTRDv1::StepManager()
 419 {
 420   //
 421   // Slow simulator. Every charged track produces electron cluster as hits
 422   // along its path across the drift volume. The step size is set acording
 423   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 424   // a spectrum taken from Ermilova et al.
 425   //
 426
 427   Int_t    pla = 0;
 428   Int_t    cha = 0;
 429   Int_t    sec = 0;
 430   Int_t    det = 0;
 431   Int_t    iPdg;
 432   Int_t    qTot;
 433
 434   Float_t  hits[3];
 435   Double_t  random[1];
 436   Float_t  charge;
 437   Float_t  aMass;
 438
 439   Double_t pTot = 0;
 440   Double_t eDelta;
 441   Double_t betaGamma, pp;
 442   Double_t stepSize;
 443
 444   Bool_t   drRegion = kFALSE;
 445   Bool_t   amRegion = kFALSE;
 446
 447   TString  cIdCurrent;
 448   TString  cIdSensDr = "J";
 449   TString  cIdSensAm = "K";
 450   Char_t   cIdChamber[3];
 451            cIdChamber[2] = 0;
 452
 453   TLorentzVector pos, mom;
 454
 455   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 456   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 457   const Int_t    kNdetsec     = kNplan * kNcham;
 458
 459   const Double_t kBig         = 1.0E+12;
 460
 461   // Ionization energy
 462   const Float_t  kWion        = 22.04;
 463   // Maximum momentum for e+ e- g
 464   const Float_t  kPTotMaxEl   = 0.002;
 465   // Minimum energy for the step size adjustment
 466   const Float_t  kEkinMinStep = 1.0e-5;
 467   // Plateau value of the energy-loss for electron in xenon
 468   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 469   //const Double_t kPlateau = 1.70;
 470   // the averaged value (26/3/99)
 471   const Float_t  kPlateau     = 1.55;
 472   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 473   const Float_t  kPrim        = 48.0;
 474   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 475   const Float_t  kPoti        = 12.1;
 476
 477   // PDG code electron
 478   const Int_t    kPdgElectron = 11;
 479
 480   // Set the maximum step size to a very large number for all
 481   // neutral particles and those outside the driftvolume
 482   gMC->SetMaxStep(kBig);
 483
 484   // Use only charged tracks
 485   if (( gMC->TrackCharge()       ) &&
 486       (!gMC->IsTrackStop()       ) &&
 487       (!gMC->IsTrackDisappeared())) {
 488
 489     // Inside a sensitive volume?
 490     drRegion = kFALSE;
 491     amRegion = kFALSE;
 492     cIdCurrent = gMC->CurrentVolName();
 493     if (cIdSensDr == cIdCurrent[1]) {
 494       drRegion = kTRUE;
 495     }
 496     if (cIdSensAm == cIdCurrent[1]) {
 497       amRegion = kTRUE;
 498     }
 499     if (drRegion || amRegion) {
 500
 501       // The hit coordinates and charge
 502       gMC->TrackPosition(pos);
 503       hits[0] = pos[0];
 504       hits[1] = pos[1];
 505       hits[2] = pos[2];
 506
 507       // The sector number (0 - 17)
 508       // The numbering goes clockwise and starts at y = 0
 509       // Not fully consistent to new corrdinate schema!!!
 510       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 511       if (phi < 90.)
 512         phi = phi + 270.;
 513       else
 514         phi = phi -  90.;
 515       sec = ((Int_t) (phi / 20));
 516
 517       // The plane and chamber number
 518       cIdChamber[0] = cIdCurrent[2];
 519       cIdChamber[1] = cIdCurrent[3];
 520       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 521       cha = ((Int_t) idChamber / kNplan);
 522       pla = ((Int_t) idChamber % kNplan);
 523
 524       // Check on selected volumes
 525       Int_t addthishit = 1;
 526       if (fSensSelect) {
 527         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 528         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 529         if (fSensSector  >= 0) {
 530           Int_t sens1  = fSensSector;
 531           Int_t sens2  = fSensSector + fSensSectorRange;
 532                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 533                        * AliTRDgeometry::Nsect();
 534           if (sens1 < sens2) {
 535             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 536           }
 537           else {
 538             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 539           }
 540         }
 541       }
 542
 543       // Add this hit
 544       if (addthishit) {
 545
 546         // The detector number
 547         det = fGeometry->GetDetector(pla,cha,sec);
 548
 549         // Special hits and TR photons only in the drift region
 550         if (drRegion) {
 551
 552           // Create a track reference at the entrance and
 553           // exit of each chamber that contain the
 554           // momentum components of the particle
 555           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
 556             gMC->TrackMomentum(mom);
 557             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 558           }
 559
 560           // Create the hits from TR photons
 561           if (fTR) CreateTRhit(det);
 562
 563         }
 564
 565         // Calculate the energy of the delta-electrons
 566         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 567         eDelta = TMath::Max(eDelta,0.0);
 568
 569         // The number of secondary electrons created
 570         qTot = ((Int_t) (eDelta / kWion) + 1);
 571
 572         // Create a new dEdx hit
 573         if (drRegion) {
 574           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,kTRUE);
 575         }
 576         else {
 577           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,kFALSE);
 578         }
 579
 580         // Calculate the maximum step size for the next tracking step
 581         // Produce only one hit if Ekin is below cutoff
 582         aMass = gMC->TrackMass();
 583         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 584
 585           // The energy loss according to Bethe Bloch
 586           iPdg  = TMath::Abs(gMC->TrackPid());
 587           if ( (iPdg != kPdgElectron) ||
 588               ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
 589             gMC->TrackMomentum(mom);
 590             pTot      = mom.Rho();
 591             betaGamma = pTot / aMass;
 592             pp        = kPrim * BetheBloch(betaGamma);
 593             // Take charge > 1 into account
 594             charge = gMC->TrackCharge();
 595             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 596           }
 597           // Electrons above 20 Mev/c are at the plateau
 598           else {
 599             pp = kPrim * kPlateau;
 600           }
 601
 602           if (pp > 0) {
 603             do
 604             gMC->GetRandom()->RndmArray(1, random);
 605             while ((random[0] == 1.) || (random[0] == 0.));
 606             stepSize = - TMath::Log(random[0]) / pp;
 607             gMC->SetMaxStep(stepSize);
 608           }
 609
 610         }
 611
 612       }
 613
 614     }
 615
 616   }
 617
 618 }
 619
 620 //_____________________________________________________________________________
 621 Double_t AliTRDv1::BetheBloch(Double_t bg)
 622 {
 623   //
 624   // Parametrization of the Bethe-Bloch-curve
 625   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
 626   //
 627
 628   // This parameters have been adjusted to averaged values from GEANT
 629   const Double_t kP1 = 7.17960e-02;
 630   const Double_t kP2 = 8.54196;
 631   const Double_t kP3 = 1.38065e-06;
 632   const Double_t kP4 = 5.30972;
 633   const Double_t kP5 = 2.83798;
 634
 635   // This parameters have been adjusted to Xe-data found in:
 636   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 637   //const Double_t kP1 = 0.76176E-1;
 638   //const Double_t kP2 = 10.632;
 639   //const Double_t kP3 = 3.17983E-6;
 640   //const Double_t kP4 = 1.8631;
 641   //const Double_t kP5 = 1.9479;
 642
 643   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
 644   const Double_t kBgMin = 0.8;
 645   const Double_t kBBMax = 6.83298;
 646   //const Double_t kBgMin = 0.6;
 647   //const Double_t kBBMax = 17.2809;
 648   //const Double_t kBgMin = 0.4;
 649   //const Double_t kBBMax = 82.0;
 650
 651   if (bg > kBgMin) {
 652     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
 653     Double_t aa = TMath::Power(yy,kP4);
 654     Double_t bb = TMath::Power((1./bg),kP5);
 655              bb = TMath::Log(kP3 + bb);
 656     return ((kP2 - aa - bb)*kP1 / aa);
 657   }
 658   else {
 659     return kBBMax;
 660   }
 661
 662 }
 663
 664 //_____________________________________________________________________________
 665 Double_t Ermilova(Double_t *x, Double_t *)
 666 {
 667   //
 668   // Calculates the delta-ray energy distribution according to Ermilova.
 669   // Logarithmic scale !
 670   //
 671
 672   Double_t energy;
 673   Double_t dpos;
 674   Double_t dnde;
 675
 676   Int_t    pos1, pos2;
 677
 678   const Int_t kNv = 31;
 679
 680   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
 681                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
 682                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
 683                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
 684                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
 685                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
 686                      ,11.5129 };
 687
 688   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
 689                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
 690                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
 691                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
 692                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
 693                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
 694                      ,  0.004 };
 695
 696   energy = x[0];
 697
 698   // Find the position
 699   pos1 = pos2 = 0;
 700   dpos = 0;
 701   do {
 702     dpos = energy - vxe[pos2++];
 703   }
 704   while (dpos > 0);
 705   pos2--;
 706   if (pos2 > kNv) pos2 = kNv - 1;
 707   pos1 = pos2 - 1;
 708
 709   // Differentiate between the sampling points
 710   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
 711
 712   return dnde;
 713
 714 }