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 "AliTRDmatrix.h"
  45 #include "AliTRDsim.h"
  46 #include "AliTRDv1.h"
  47 #include "AliMC.h"
  48
  49 ClassImp(AliTRDv1)
  50
  51 //_____________________________________________________________________________
  52 AliTRDv1::AliTRDv1():AliTRD()
  53 {
  54   //
  55   // Default constructor
  56   //
  57
  58   fSensSelect      =  0;
  59   fSensPlane       = -1;
  60   fSensChamber     = -1;
  61   fSensSector      = -1;
  62   fSensSectorRange =  0;
  63
  64   fDeltaE          = NULL;
  65   fTR              = NULL;
  66
  67 }
  68
  69 //_____________________________________________________________________________
  70 AliTRDv1::AliTRDv1(const char *name, const char *title)
  71          :AliTRD(name, title)
  72 {
  73   //
  74   // Standard constructor for Transition Radiation Detector version 1
  75   //
  76
  77   fSensSelect      =  0;
  78   fSensPlane       = -1;
  79   fSensChamber     = -1;
  80   fSensSector      = -1;
  81   fSensSectorRange =  0;
  82
  83   fDeltaE          = NULL;
  84   fTR              = NULL;
  85
  86   SetBufferSize(128000);
  87
  88 }
  89
  90 //_____________________________________________________________________________
  91 AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
  92 {
  93   //
  94   // Copy constructor
  95   //
  96
  97   ((AliTRDv1 &) trd).Copy(*this);
  98
  99 }
 100
 101 //_____________________________________________________________________________
 102 AliTRDv1::~AliTRDv1()
 103 {
 104   //
 105   // AliTRDv1 destructor
 106   //
 107
 108   if (fDeltaE) delete fDeltaE;
 109   if (fTR)     delete fTR;
 110
 111 }
 112
 113 //_____________________________________________________________________________
 114 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
 115 {
 116   //
 117   // Assignment operator
 118   //
 119
 120   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
 121   return *this;
 122
 123 }
 124
 125 //_____________________________________________________________________________
 126 void AliTRDv1::Copy(TObject &trd)
 127 {
 128   //
 129   // Copy function
 130   //
 131
 132   ((AliTRDv1 &) trd).fSensSelect      = fSensSelect;
 133   ((AliTRDv1 &) trd).fSensPlane       = fSensPlane;
 134   ((AliTRDv1 &) trd).fSensChamber     = fSensChamber;
 135   ((AliTRDv1 &) trd).fSensSector      = fSensSector;
 136   ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
 137
 138   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
 139   fTR->Copy(*((AliTRDv1 &) trd).fTR);
 140
 141 }
 142
 143 //_____________________________________________________________________________
 144 void AliTRDv1::CreateGeometry()
 145 {
 146   //
 147   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 148   // This version covers the full azimuth.
 149   //
 150
 151   // Check that FRAME is there otherwise we have no place where to put the TRD
 152   AliModule* frame = gAlice->GetModule("FRAME");
 153   if (!frame) return;
 154
 155   // Define the chambers
 156   AliTRD::CreateGeometry();
 157
 158 }
 159
 160 //_____________________________________________________________________________
 161 void AliTRDv1::CreateMaterials()
 162 {
 163   //
 164   // Create materials for the Transition Radiation Detector version 1
 165   //
 166
 167   AliTRD::CreateMaterials();
 168
 169 }
 170
 171 //_____________________________________________________________________________
 172 void AliTRDv1::CreateTRhit(Int_t det)
 173 {
 174   //
 175   // Creates an electron cluster from a TR photon.
 176   // The photon is assumed to be created a the end of the radiator. The
 177   // distance after which it deposits its energy takes into account the
 178   // absorbtion of the entrance window and of the gas mixture in drift
 179   // volume.
 180   //
 181
 182   // PDG code electron
 183   const Int_t   kPdgElectron = 11;
 184
 185   // Ionization energy
 186   const Float_t kWion        = 22.04;
 187
 188   // Maximum number of TR photons per track
 189   const Int_t   kNTR         = 50;
 190
 191   TLorentzVector mom, pos;
 192
 193   // Create TR at the entrance of the chamber
 194   if (gMC->IsTrackEntering()) {
 195
 196     // Create TR only for electrons
 197     Int_t iPdg = gMC->TrackPid();
 198     if (TMath::Abs(iPdg) != kPdgElectron) return;
 199
 200     Float_t eTR[kNTR];
 201     Int_t   nTR;
 202
 203     // Create TR photons
 204     gMC->TrackMomentum(mom);
 205     Float_t pTot = mom.Rho();
 206     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
 207     if (nTR > kNTR) {
 208       printf("AliTRDv1::CreateTRhit -- ");
 209       printf("Boundary error: nTR = %d, kNTR = %d\n",nTR,kNTR);
 210       exit(1);
 211     }
 212
 213     // Loop through the TR photons
 214     for (Int_t iTR = 0; iTR < nTR; iTR++) {
 215
 216       Float_t energyMeV = eTR[iTR] * 0.001;
 217       Float_t energyeV  = eTR[iTR] * 1000.0;
 218       Float_t absLength = 0;
 219       Float_t sigma     = 0;
 220
 221       // Take the absorbtion in the entrance window into account
 222       Double_t muMy = fTR->GetMuMy(energyMeV);
 223       sigma = muMy * fFoilDensity;
 224       if (sigma > 0.0) {
 225         absLength = gRandom->Exp(1.0/sigma);
 226         if (absLength < AliTRDgeometry::MyThick()) continue;
 227       }
 228       else {
 229         continue;
 230       }
 231
 232       // The absorbtion cross sections in the drift gas
 233       if (fGasMix == 1) {
 234         // Gas-mixture (Xe/CO2)
 235         Double_t muXe = fTR->GetMuXe(energyMeV);
 236         Double_t muCO = fTR->GetMuCO(energyMeV);
 237         sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
 238       }
 239       else {
 240         // Gas-mixture (Xe/Isobutane)
 241         Double_t muXe = fTR->GetMuXe(energyMeV);
 242         Double_t muBu = fTR->GetMuBu(energyMeV);
 243         sigma = (0.97 * muXe + 0.03 * muBu) * fGasDensity * fTR->GetTemp();
 244       }
 245
 246       // The distance after which the energy of the TR photon
 247       // is deposited.
 248       if (sigma > 0.0) {
 249         absLength = gRandom->Exp(1.0/sigma);
 250         if (absLength > AliTRDgeometry::DrThick()) continue;
 251       }
 252       else {
 253         continue;
 254       }
 255
 256       // The position of the absorbtion
 257       Float_t posHit[3];
 258       gMC->TrackPosition(pos);
 259       posHit[0] = pos[0] + mom[0] / pTot * absLength;
 260       posHit[1] = pos[1] + mom[1] / pTot * absLength;
 261       posHit[2] = pos[2] + mom[2] / pTot * absLength;
 262
 263       // Create the charge
 264       Int_t q = ((Int_t) (energyeV / kWion));
 265
 266       // Add the hit to the array. TR photon hits are marked
 267       // by negative charge
 268       AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
 269
 270     }
 271
 272   }
 273
 274 }
 275
 276 //_____________________________________________________________________________
 277 void AliTRDv1::Init()
 278 {
 279   //
 280   // Initialise Transition Radiation Detector after geometry has been built.
 281   //
 282
 283   AliTRD::Init();
 284
 285   if(fDebug) printf("%s: Slow simulator\n",ClassName());
 286   if (fSensSelect) {
 287     if (fSensPlane   >= 0)
 288       printf("          Only plane %d is sensitive\n",fSensPlane);
 289     if (fSensChamber >= 0)
 290       printf("          Only chamber %d is sensitive\n",fSensChamber);
 291     if (fSensSector  >= 0) {
 292       Int_t sens1  = fSensSector;
 293       Int_t sens2  = fSensSector + fSensSectorRange;
 294             sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 295                    * AliTRDgeometry::Nsect();
 296       printf("          Only sectors %d - %d are sensitive\n",sens1,sens2-1);
 297     }
 298   }
 299   if (fTR)
 300     printf("%s: TR simulation on\n",ClassName());
 301   else
 302     printf("%s: TR simulation off\n",ClassName());
 303   printf("\n");
 304
 305   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 306   const Float_t kPoti = 12.1;
 307   // Maximum energy (50 keV);
 308   const Float_t kEend = 50000.0;
 309   // Ermilova distribution for the delta-ray spectrum
 310   Float_t poti = TMath::Log(kPoti);
 311   Float_t eEnd = TMath::Log(kEend);
 312   fDeltaE = new TF1("deltae",Ermilova,poti,eEnd,0);
 313
 314   if(fDebug) {
 315     printf("%s: ",ClassName());
 316     for (Int_t i = 0; i < 80; i++) printf("*");
 317     printf("\n");
 318   }
 319
 320 }
 321
 322 //_____________________________________________________________________________
 323 AliTRDsim *AliTRDv1::CreateTR()
 324 {
 325   //
 326   // Enables the simulation of TR
 327   //
 328
 329   fTR = new AliTRDsim();
 330   return fTR;
 331
 332 }
 333
 334 //_____________________________________________________________________________
 335 void AliTRDv1::SetSensPlane(Int_t iplane)
 336 {
 337   //
 338   // Defines the hit-sensitive plane (0-5)
 339   //
 340
 341   if ((iplane < 0) || (iplane > 5)) {
 342     printf("Wrong input value: %d\n",iplane);
 343     printf("Use standard setting\n");
 344     fSensPlane  = -1;
 345     fSensSelect =  0;
 346     return;
 347   }
 348
 349   fSensSelect = 1;
 350   fSensPlane  = iplane;
 351
 352 }
 353
 354 //_____________________________________________________________________________
 355 void AliTRDv1::SetSensChamber(Int_t ichamber)
 356 {
 357   //
 358   // Defines the hit-sensitive chamber (0-4)
 359   //
 360
 361   if ((ichamber < 0) || (ichamber > 4)) {
 362     printf("Wrong input value: %d\n",ichamber);
 363     printf("Use standard setting\n");
 364     fSensChamber = -1;
 365     fSensSelect  =  0;
 366     return;
 367   }
 368
 369   fSensSelect  = 1;
 370   fSensChamber = ichamber;
 371
 372 }
 373
 374 //_____________________________________________________________________________
 375 void AliTRDv1::SetSensSector(Int_t isector)
 376 {
 377   //
 378   // Defines the hit-sensitive sector (0-17)
 379   //
 380
 381   SetSensSector(isector,1);
 382
 383 }
 384
 385 //_____________________________________________________________________________
 386 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
 387 {
 388   //
 389   // Defines a range of hit-sensitive sectors. The range is defined by
 390   // <isector> (0-17) as the starting point and <nsector> as the number
 391   // of sectors to be included.
 392   //
 393
 394   if ((isector < 0) || (isector > 17)) {
 395     printf("Wrong input value <isector>: %d\n",isector);
 396     printf("Use standard setting\n");
 397     fSensSector      = -1;
 398     fSensSectorRange =  0;
 399     fSensSelect      =  0;
 400     return;
 401   }
 402
 403   if ((nsector < 1) || (nsector > 18)) {
 404     printf("Wrong input value <nsector>: %d\n",nsector);
 405     printf("Use standard setting\n");
 406     fSensSector      = -1;
 407     fSensSectorRange =  0;
 408     fSensSelect      =  0;
 409     return;
 410   }
 411
 412   fSensSelect      = 1;
 413   fSensSector      = isector;
 414   fSensSectorRange = nsector;
 415
 416 }
 417
 418 //_____________________________________________________________________________
 419 void AliTRDv1::StepManager()
 420 {
 421   //
 422   // Slow simulator. Every charged track produces electron cluster as hits
 423   // along its path across the drift volume. The step size is set acording
 424   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 425   // a spectrum taken from Ermilova et al.
 426   //
 427
 428   Int_t    pla = 0;
 429   Int_t    cha = 0;
 430   Int_t    sec = 0;
 431   Int_t    det = 0;
 432   Int_t    iPdg;
 433   Int_t    qTot;
 434
 435   Float_t  hits[3];
 436   Double_t  random[1];
 437   Float_t  charge;
 438   Float_t  aMass;
 439
 440   Double_t pTot = 0;
 441   Double_t eDelta;
 442   Double_t betaGamma, pp;
 443   Double_t stepSize;
 444
 445   Bool_t   drRegion = kFALSE;
 446   Bool_t   amRegion = kFALSE;
 447
 448   TString  cIdCurrent;
 449   TString  cIdSensDr = "J";
 450   TString  cIdSensAm = "K";
 451   Char_t   cIdChamber[3];
 452            cIdChamber[2] = 0;
 453
 454   TLorentzVector pos, mom;
 455
 456   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 457   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 458   const Int_t    kNdetsec     = kNplan * kNcham;
 459
 460   const Double_t kBig         = 1.0E+12;
 461
 462   // Ionization energy
 463   const Float_t  kWion        = 22.04;
 464   // Maximum momentum for e+ e- g
 465   const Float_t  kPTotMaxEl   = 0.002;
 466   // Minimum energy for the step size adjustment
 467   const Float_t  kEkinMinStep = 1.0e-5;
 468   // Plateau value of the energy-loss for electron in xenon
 469   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 470   //const Double_t kPlateau = 1.70;
 471   // the averaged value (26/3/99)
 472   const Float_t  kPlateau     = 1.55;
 473   // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 474   const Float_t  kPrim        = 48.0;
 475   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 476   const Float_t  kPoti        = 12.1;
 477
 478   // PDG code electron
 479   const Int_t    kPdgElectron = 11;
 480
 481   // Set the maximum step size to a very large number for all
 482   // neutral particles and those outside the driftvolume
 483   gMC->SetMaxStep(kBig);
 484
 485   // Use only charged tracks
 486   if (( gMC->TrackCharge()       ) &&
 487       (!gMC->IsTrackStop()       ) &&
 488       (!gMC->IsTrackDisappeared())) {
 489
 490     // Inside a sensitive volume?
 491     drRegion = kFALSE;
 492     amRegion = kFALSE;
 493     cIdCurrent = gMC->CurrentVolName();
 494     if (cIdSensDr == cIdCurrent[1]) {
 495       drRegion = kTRUE;
 496     }
 497     if (cIdSensAm == cIdCurrent[1]) {
 498       amRegion = kTRUE;
 499     }
 500     if (drRegion || amRegion) {
 501
 502       // The hit coordinates and charge
 503       gMC->TrackPosition(pos);
 504       hits[0] = pos[0];
 505       hits[1] = pos[1];
 506       hits[2] = pos[2];
 507
 508       // The sector number (0 - 17)
 509       // The numbering goes clockwise and starts at y = 0
 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 }