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       // Not fully consistent to new corrdinate schema!!!
 511       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 512       if (phi < 90.)
 513         phi = phi + 270.;
 514       else
 515         phi = phi -  90.;
 516       sec = ((Int_t) (phi / 20));
 517
 518       // The plane and chamber number
 519       cIdChamber[0] = cIdCurrent[2];
 520       cIdChamber[1] = cIdCurrent[3];
 521       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 522       cha = ((Int_t) idChamber / kNplan);
 523       pla = ((Int_t) idChamber % kNplan);
 524
 525       // Check on selected volumes
 526       Int_t addthishit = 1;
 527       if (fSensSelect) {
 528         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 529         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 530         if (fSensSector  >= 0) {
 531           Int_t sens1  = fSensSector;
 532           Int_t sens2  = fSensSector + fSensSectorRange;
 533                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 534                        * AliTRDgeometry::Nsect();
 535           if (sens1 < sens2) {
 536             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 537           }
 538           else {
 539             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 540           }
 541         }
 542       }
 543
 544       // Add this hit
 545       if (addthishit) {
 546
 547         // The detector number
 548         det = fGeometry->GetDetector(pla,cha,sec);
 549
 550         // Special hits and TR photons only in the drift region
 551         if (drRegion) {
 552
 553           // Create a track reference at the entrance and
 554           // exit of each chamber that contain the
 555           // momentum components of the particle
 556           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
 557             gMC->TrackMomentum(mom);
 558             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 559           }
 560
 561           // Create the hits from TR photons
 562           if (fTR) CreateTRhit(det);
 563
 564         }
 565
 566         // Calculate the energy of the delta-electrons
 567         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 568         eDelta = TMath::Max(eDelta,0.0);
 569
 570         // The number of secondary electrons created
 571         qTot = ((Int_t) (eDelta / kWion) + 1);
 572
 573         // Create a new dEdx hit
 574         if (drRegion) {
 575           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,kTRUE);
 576         }
 577         else {
 578           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,kFALSE);
 579         }
 580
 581         // Calculate the maximum step size for the next tracking step
 582         // Produce only one hit if Ekin is below cutoff
 583         aMass = gMC->TrackMass();
 584         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 585
 586           // The energy loss according to Bethe Bloch
 587           iPdg  = TMath::Abs(gMC->TrackPid());
 588           if ( (iPdg != kPdgElectron) ||
 589               ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
 590             gMC->TrackMomentum(mom);
 591             pTot      = mom.Rho();
 592             betaGamma = pTot / aMass;
 593             pp        = kPrim * BetheBloch(betaGamma);
 594             // Take charge > 1 into account
 595             charge = gMC->TrackCharge();
 596             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 597           }
 598           // Electrons above 20 Mev/c are at the plateau
 599           else {
 600             pp = kPrim * kPlateau;
 601           }
 602
 603           if (pp > 0) {
 604             do
 605             gMC->GetRandom()->RndmArray(1, random);
 606             while ((random[0] == 1.) || (random[0] == 0.));
 607             stepSize = - TMath::Log(random[0]) / pp;
 608             gMC->SetMaxStep(stepSize);
 609           }
 610
 611         }
 612
 613       }
 614
 615     }
 616
 617   }
 618
 619 }
 620
 621 //_____________________________________________________________________________
 622 Double_t AliTRDv1::BetheBloch(Double_t bg)
 623 {
 624   //
 625   // Parametrization of the Bethe-Bloch-curve
 626   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
 627   //
 628
 629   // This parameters have been adjusted to averaged values from GEANT
 630   const Double_t kP1 = 7.17960e-02;
 631   const Double_t kP2 = 8.54196;
 632   const Double_t kP3 = 1.38065e-06;
 633   const Double_t kP4 = 5.30972;
 634   const Double_t kP5 = 2.83798;
 635
 636   // This parameters have been adjusted to Xe-data found in:
 637   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 638   //const Double_t kP1 = 0.76176E-1;
 639   //const Double_t kP2 = 10.632;
 640   //const Double_t kP3 = 3.17983E-6;
 641   //const Double_t kP4 = 1.8631;
 642   //const Double_t kP5 = 1.9479;
 643
 644   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
 645   const Double_t kBgMin = 0.8;
 646   const Double_t kBBMax = 6.83298;
 647   //const Double_t kBgMin = 0.6;
 648   //const Double_t kBBMax = 17.2809;
 649   //const Double_t kBgMin = 0.4;
 650   //const Double_t kBBMax = 82.0;
 651
 652   if (bg > kBgMin) {
 653     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
 654     Double_t aa = TMath::Power(yy,kP4);
 655     Double_t bb = TMath::Power((1./bg),kP5);
 656              bb = TMath::Log(kP3 + bb);
 657     return ((kP2 - aa - bb)*kP1 / aa);
 658   }
 659   else {
 660     return kBBMax;
 661   }
 662
 663 }
 664
 665 //_____________________________________________________________________________
 666 Double_t Ermilova(Double_t *x, Double_t *)
 667 {
 668   //
 669   // Calculates the delta-ray energy distribution according to Ermilova.
 670   // Logarithmic scale !
 671   //
 672
 673   Double_t energy;
 674   Double_t dpos;
 675   Double_t dnde;
 676
 677   Int_t    pos1, pos2;
 678
 679   const Int_t kNv = 31;
 680
 681   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
 682                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
 683                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
 684                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
 685                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
 686                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
 687                      ,11.5129 };
 688
 689   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
 690                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
 691                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
 692                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
 693                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
 694                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
 695                      ,  0.004 };
 696
 697   energy = x[0];
 698
 699   // Find the position
 700   pos1 = pos2 = 0;
 701   dpos = 0;
 702   do {
 703     dpos = energy - vxe[pos2++];
 704   }
 705   while (dpos > 0);
 706   pos2--;
 707   if (pos2 > kNv) pos2 = kNv - 1;
 708   pos1 = pos2 - 1;
 709
 710   // Differentiate between the sampling points
 711   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
 712
 713   return dnde;
 714
 715 }