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 "AliLog.h"
  42 #include "AliMC.h"
  43 #include "AliRun.h"
  44 #include "AliTRDgeometry.h"
  45 #include "AliTRDhit.h"
  46 #include "AliTRDsim.h"
  47 #include "AliTRDv1.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   fDeltaG            = NULL;
  66   fTR                = NULL;
  67   fTRon              = kFALSE;
  68
  69   fStepSize          = 0.1;
  70   fTypeOfStepManager = 1;
  71
  72 }
  73
  74 //_____________________________________________________________________________
  75 AliTRDv1::AliTRDv1(const char *name, const char *title)
  76          :AliTRD(name, title)
  77 {
  78   //
  79   // Standard constructor for Transition Radiation Detector version 1
  80   //
  81
  82   fSensSelect        =  0;
  83   fSensPlane         = -1;
  84   fSensChamber       = -1;
  85   fSensSector        = -1;
  86   fSensSectorRange   =  0;
  87
  88   fDeltaE            = NULL;
  89   fDeltaG            = NULL;
  90   fTR                = NULL;
  91   fTRon              = kTRUE;
  92   fStepSize          = 0.1;
  93   fTypeOfStepManager = 1;
  94
  95   SetBufferSize(128000);
  96
  97 }
  98
  99 //_____________________________________________________________________________
 100 AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
 101 {
 102   //
 103   // Copy constructor
 104   //
 105
 106   ((AliTRDv1 &) trd).Copy(*this);
 107
 108 }
 109
 110 //_____________________________________________________________________________
 111 AliTRDv1::~AliTRDv1()
 112 {
 113   //
 114   // AliTRDv1 destructor
 115   //
 116
 117   if (fDeltaE) delete fDeltaE;
 118   if (fDeltaG) delete fDeltaG;
 119   if (fTR)     delete fTR;
 120
 121 }
 122
 123 //_____________________________________________________________________________
 124 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
 125 {
 126   //
 127   // Assignment operator
 128   //
 129
 130   if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
 131   return *this;
 132
 133 }
 134
 135 //_____________________________________________________________________________
 136 void AliTRDv1::Copy(TObject &trd) const
 137 {
 138         printf("void AliTRDv1::Copy(TObject &trd) const\n");
 139   //
 140   // Copy function
 141   //
 142
 143   ((AliTRDv1 &) trd).fSensSelect        = fSensSelect;
 144   ((AliTRDv1 &) trd).fSensPlane         = fSensPlane;
 145   ((AliTRDv1 &) trd).fSensChamber       = fSensChamber;
 146   ((AliTRDv1 &) trd).fSensSector        = fSensSector;
 147   ((AliTRDv1 &) trd).fSensSectorRange   = fSensSectorRange;
 148
 149   ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
 150   ((AliTRDv1 &) trd).fStepSize          = fStepSize;
 151
 152   ((AliTRDv1 &) trd).fTRon              = fTRon;
 153
 154   fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
 155   fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
 156   fTR->Copy(*((AliTRDv1 &) trd).fTR);
 157
 158 }
 159
 160 //_____________________________________________________________________________
 161 void AliTRDv1::CreateGeometry()
 162 {
 163   //
 164   // Create the GEANT geometry for the Transition Radiation Detector - Version 1
 165   // This version covers the full azimuth.
 166   //
 167
 168   // Check that FRAME is there otherwise we have no place where to put the TRD
 169   AliModule* frame = gAlice->GetModule("FRAME");
 170   if (!frame) return;
 171
 172   // Define the chambers
 173   AliTRD::CreateGeometry();
 174
 175 }
 176
 177 //_____________________________________________________________________________
 178 void AliTRDv1::CreateMaterials()
 179 {
 180   //
 181   // Create materials for the Transition Radiation Detector version 1
 182   //
 183
 184   AliTRD::CreateMaterials();
 185
 186 }
 187
 188 //_____________________________________________________________________________
 189 void AliTRDv1::CreateTRhit(Int_t det)
 190 {
 191   //
 192   // Creates an electron cluster from a TR photon.
 193   // The photon is assumed to be created a the end of the radiator. The
 194   // distance after which it deposits its energy takes into account the
 195   // absorbtion of the entrance window and of the gas mixture in drift
 196   // volume.
 197   //
 198
 199   // PDG code electron
 200   const Int_t   kPdgElectron = 11;
 201
 202   // Ionization energy
 203   const Float_t kWion        = 23.53;
 204
 205   // Maximum number of TR photons per track
 206   const Int_t   kNTR         = 50;
 207
 208   TLorentzVector mom, pos;
 209
 210   // Create TR at the entrance of the chamber
 211   if (gMC->IsTrackEntering()) {
 212
 213     // Create TR only for electrons
 214     Int_t iPdg = gMC->TrackPid();
 215     if (TMath::Abs(iPdg) != kPdgElectron) return;
 216
 217     Float_t eTR[kNTR];
 218     Int_t   nTR;
 219
 220     // Create TR photons
 221     gMC->TrackMomentum(mom);
 222     Float_t pTot = mom.Rho();
 223     fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
 224     if (nTR > kNTR) {
 225       AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
 226     }
 227
 228     // Loop through the TR photons
 229     for (Int_t iTR = 0; iTR < nTR; iTR++) {
 230
 231       Float_t energyMeV = eTR[iTR] * 0.001;
 232       Float_t energyeV  = eTR[iTR] * 1000.0;
 233       Float_t absLength = 0;
 234       Float_t sigma     = 0;
 235
 236       // Take the absorbtion in the entrance window into account
 237       Double_t muMy = fTR->GetMuMy(energyMeV);
 238       sigma = muMy * fFoilDensity;
 239       if (sigma > 0.0) {
 240         absLength = gRandom->Exp(1.0/sigma);
 241         if (absLength < AliTRDgeometry::MyThick()) continue;
 242       }
 243       else {
 244         continue;
 245       }
 246
 247       // The absorbtion cross sections in the drift gas
 248       // Gas-mixture (Xe/CO2)
 249       Double_t muXe = fTR->GetMuXe(energyMeV);
 250       Double_t muCO = fTR->GetMuCO(energyMeV);
 251       sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
 252
 253       // The distance after which the energy of the TR photon
 254       // is deposited.
 255       if (sigma > 0.0) {
 256         absLength = gRandom->Exp(1.0/sigma);
 257         if (absLength > (AliTRDgeometry::DrThick()
 258                        + AliTRDgeometry::AmThick())) {
 259           continue;
 260         }
 261       }
 262       else {
 263         continue;
 264       }
 265
 266       // The position of the absorbtion
 267       Float_t posHit[3];
 268       gMC->TrackPosition(pos);
 269       posHit[0] = pos[0] + mom[0] / pTot * absLength;
 270       posHit[1] = pos[1] + mom[1] / pTot * absLength;
 271       posHit[2] = pos[2] + mom[2] / pTot * absLength;
 272
 273       // Create the charge
 274       Int_t q = ((Int_t) (energyeV / kWion));
 275
 276       // Add the hit to the array. TR photon hits are marked
 277       // by negative charge
 278      AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
 279
 280     }
 281
 282   }
 283
 284 }
 285
 286 //_____________________________________________________________________________
 287 void AliTRDv1::Init()
 288 {
 289   //
 290   // Initialise Transition Radiation Detector after geometry has been built.
 291   //
 292
 293   AliTRD::Init();
 294
 295   AliDebug(1,"Slow simulator\n");
 296   if (fSensSelect) {
 297     if (fSensPlane   >= 0)
 298       AliInfo(Form("Only plane %d is sensitive"));
 299     if (fSensChamber >= 0)
 300       AliInfo(Form("Only chamber %d is sensitive",fSensChamber));
 301     if (fSensSector  >= 0) {
 302       Int_t sens1  = fSensSector;
 303       Int_t sens2  = fSensSector + fSensSectorRange;
 304             sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 305                    * AliTRDgeometry::Nsect();
 306             AliInfo(Form("Only sectors %d - %d are sensitive\n",sens1,sens2-1));
 307     }
 308   }
 309
 310   // Switch on TR simulation as default
 311   if (!fTRon) {
 312     AliInfo("TR simulation off");
 313   }
 314   else {
 315     fTR = new AliTRDsim();
 316   }
 317
 318   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 319   const Float_t kPoti = 12.1;
 320   // Maximum energy (50 keV);
 321   const Float_t kEend = 50000.0;
 322   // Ermilova distribution for the delta-ray spectrum
 323   Float_t poti = TMath::Log(kPoti);
 324   Float_t eEnd = TMath::Log(kEend);
 325
 326   // Ermilova distribution for the delta-ray spectrum
 327   fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
 328
 329   // Geant3 distribution for the delta-ray spectrum
 330   fDeltaG = new TF1("deltag",IntSpecGeant,2.421257,28.536469,0);
 331
 332   AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
 333
 334 }
 335
 336 //_____________________________________________________________________________
 337 void AliTRDv1::SetSensPlane(Int_t iplane)
 338 {
 339   //
 340   // Defines the hit-sensitive plane (0-5)
 341   //
 342
 343   if ((iplane < 0) || (iplane > 5)) {
 344     AliWarning(Form("Wrong input value:%d",iplane));
 345     AliWarning("Use standard setting");
 346     fSensPlane  = -1;
 347     fSensSelect =  0;
 348     return;
 349   }
 350
 351   fSensSelect = 1;
 352   fSensPlane  = iplane;
 353
 354 }
 355
 356 //_____________________________________________________________________________
 357 void AliTRDv1::SetSensChamber(Int_t ichamber)
 358 {
 359   //
 360   // Defines the hit-sensitive chamber (0-4)
 361   //
 362
 363   if ((ichamber < 0) || (ichamber > 4)) {
 364     AliWarning(Form("Wrong input value: %d",ichamber));
 365     AliWarning("Use standard setting");
 366     fSensChamber = -1;
 367     fSensSelect  =  0;
 368     return;
 369   }
 370
 371   fSensSelect  = 1;
 372   fSensChamber = ichamber;
 373
 374 }
 375
 376 //_____________________________________________________________________________
 377 void AliTRDv1::SetSensSector(Int_t isector)
 378 {
 379   //
 380   // Defines the hit-sensitive sector (0-17)
 381   //
 382
 383   SetSensSector(isector,1);
 384
 385 }
 386
 387 //_____________________________________________________________________________
 388 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
 389 {
 390   //
 391   // Defines a range of hit-sensitive sectors. The range is defined by
 392   // <isector> (0-17) as the starting point and <nsector> as the number
 393   // of sectors to be included.
 394   //
 395
 396   if ((isector < 0) || (isector > 17)) {
 397     AliWarning(Form("Wrong input value <isector>: %d",isector));
 398     AliWarning("Use standard setting");
 399     fSensSector      = -1;
 400     fSensSectorRange =  0;
 401     fSensSelect      =  0;
 402     return;
 403   }
 404
 405   if ((nsector < 1) || (nsector > 18)) {
 406     AliWarning(Form("Wrong input value <nsector>: %d",nsector));
 407     AliWarning("Use standard setting");
 408     fSensSector      = -1;
 409     fSensSectorRange =  0;
 410     fSensSelect      =  0;
 411     return;
 412   }
 413
 414   fSensSelect      = 1;
 415   fSensSector      = isector;
 416   fSensSectorRange = nsector;
 417
 418 }
 419
 420 //_____________________________________________________________________________
 421 void AliTRDv1::StepManager()
 422 {
 423   //
 424   // Slow simulator. Every charged track produces electron cluster as hits
 425   // along its path across the drift volume.
 426   //
 427
 428   switch (fTypeOfStepManager) {
 429     case 0  : StepManagerErmilova();  break;  // 0 is Ermilova
 430     case 1  : StepManagerGeant();     break;  // 1 is Geant
 431     case 2  : StepManagerFixedStep(); break;  // 2 is fixed step
 432     default : AliWarning("Not a valid Step Manager.");
 433   }
 434
 435 }
 436
 437 //_____________________________________________________________________________
 438 void AliTRDv1::SelectStepManager(Int_t t)
 439 {
 440   //
 441   // Selects a step manager type:
 442   //   0 - Ermilova
 443   //   1 - Geant3
 444   //   2 - Fixed step size
 445   //
 446
 447 /*  if (t == 1) {
 448     AliWarning("Sorry, Geant parametrization step manager is not implemented yet. Please ask K.Oyama for detail.");
 449   }
 450 */
 451   fTypeOfStepManager = t;
 452   AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
 453
 454 }
 455
 456 //_____________________________________________________________________________
 457 void AliTRDv1::StepManagerGeant()
 458 {
 459   //
 460   // Slow simulator. Every charged track produces electron cluster as hits
 461   // along its path across the drift volume. The step size is set acording
 462   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 463   // a spectrum taken from Geant3.
 464   //
 465   Int_t    pla = 0;
 466   Int_t    cha = 0;
 467   Int_t    sec = 0;
 468   Int_t    det = 0;
 469   Int_t    iPdg;
 470   Int_t    qTot;
 471
 472   Float_t  hits[3];
 473   Float_t  charge;
 474   Float_t  aMass;
 475
 476   Double_t pTot = 0;
 477   Double_t eDelta;
 478   Double_t betaGamma, pp;
 479   Double_t stepSize=0;
 480
 481   Bool_t   drRegion = kFALSE;
 482   Bool_t   amRegion = kFALSE;
 483
 484   TString  cIdCurrent;
 485   TString  cIdSensDr = "J";
 486   TString  cIdSensAm = "K";
 487   Char_t   cIdChamber[3];
 488            cIdChamber[2] = 0;
 489
 490   TLorentzVector pos, mom;
 491
 492   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 493   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 494   const Int_t    kNdetsec     = kNplan * kNcham;
 495
 496   const Double_t kBig         = 1.0E+12; // Infinitely big
 497   const Float_t  kWion        = 23.53;   // Ionization energy
 498   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 499
 500   // Minimum energy for the step size adjustment
 501   const Float_t  kEkinMinStep = 1.0e-5;
 502   // energy threshold for production of delta electrons
 503         const Float_t  kECut = 1.0e4;
 504         // Parameters entering the parametrized range for delta electrons
 505         const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
 506         // Gas density -> To be made user adjustable !
 507         const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
 508
 509   // Plateau value of the energy-loss for electron in xenon
 510   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 511   //const Double_t kPlateau = 1.70;
 512   // the averaged value (26/3/99)
 513   const Float_t  kPlateau     = 1.55;
 514
 515   const Float_t  kPrim        = 19.34;  // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 516   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 517   const Float_t  kPoti        = 12.1;
 518
 519   const Int_t    kPdgElectron = 11;  // PDG code electron
 520
 521   // Set the maximum step size to a very large number for all
 522   // neutral particles and those outside the driftvolume
 523   gMC->SetMaxStep(kBig);
 524
 525   // Use only charged tracks
 526   if (( gMC->TrackCharge()       ) &&
 527       (!gMC->IsTrackStop()       ) &&
 528       (!gMC->IsTrackDisappeared())) {
 529
 530     // Inside a sensitive volume?
 531     drRegion = kFALSE;
 532     amRegion = kFALSE;
 533     cIdCurrent = gMC->CurrentVolName();
 534     if (cIdSensDr == cIdCurrent[1]) {
 535       drRegion = kTRUE;
 536     }
 537     if (cIdSensAm == cIdCurrent[1]) {
 538       amRegion = kTRUE;
 539     }
 540     if (drRegion || amRegion) {
 541
 542       // The hit coordinates and charge
 543       gMC->TrackPosition(pos);
 544       hits[0] = pos[0];
 545       hits[1] = pos[1];
 546       hits[2] = pos[2];
 547
 548       // The sector number (0 - 17)
 549       // The numbering goes clockwise and starts at y = 0
 550       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 551       if (phi < 90.)
 552         phi = phi + 270.;
 553       else
 554         phi = phi -  90.;
 555       sec = ((Int_t) (phi / 20));
 556
 557       // The plane and chamber number
 558       cIdChamber[0] = cIdCurrent[2];
 559       cIdChamber[1] = cIdCurrent[3];
 560       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 561       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 562       pla = ((Int_t) idChamber % kNplan);
 563
 564       // Check on selected volumes
 565       Int_t addthishit = 1;
 566       if (fSensSelect) {
 567         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 568         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 569         if (fSensSector  >= 0) {
 570           Int_t sens1  = fSensSector;
 571           Int_t sens2  = fSensSector + fSensSectorRange;
 572                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 573                        * AliTRDgeometry::Nsect();
 574           if (sens1 < sens2) {
 575             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 576                                 }
 577           else {
 578             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 579                                 }
 580                                 }
 581       }
 582
 583       // Add this hit
 584       if (addthishit) {
 585
 586         // The detector number
 587         det = fGeometry->GetDetector(pla,cha,sec);
 588
 589         // Special hits only in the drift region
 590         if (drRegion) {
 591           // Create a track reference at the entrance and
 592           // exit of each chamber that contain the
 593                 // momentum components of the particle
 594           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
 595             gMC->TrackMomentum(mom);
 596             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 597           }
 598
 599                                         if (gMC->IsTrackEntering() && !gMC->IsNewTrack()) {
 600                                                 // determine if hit belong to primary track
 601                                                 fPrimaryTrackPid=gAlice->GetMCApp()->GetCurrentTrackNumber();
 602                                                 //determine track length when entering the detector
 603                                                 fTrackLength0=gMC->TrackLength();
 604                                         }
 605
 606                                         // Create the hits from TR photons
 607           if (fTR) CreateTRhit(det);
 608                                 }
 609
 610                                 // Calculate the energy of the delta-electrons
 611                                 // modified by Alex Bercuci (A.Bercuci@gsi.de) on 26.01.06
 612                                 // take into account correlation with the underlying GEANT tracking
 613                                 // mechanism. see
 614         // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
 615
 616                                 // determine the most significant process (last on the processes list)
 617                                 // which caused this hit
 618         TArrayI processes;
 619         gMC->StepProcesses(processes);
 620         int nofprocesses=processes.GetSize(), pid;
 621                                 if(!nofprocesses) pid=0;
 622                                 else pid=       processes[nofprocesses-1];
 623
 624                                 // generate Edep according to GEANT parametrisation
 625                                 eDelta =TMath::Exp(fDeltaG->GetRandom()) - kPoti;
 626         eDelta=TMath::Max(eDelta,0.0);
 627                                 float pr_range=0.;
 628                                 float range=gMC->TrackLength()-fTrackLength0;
 629                                 // merge GEANT tracker information with localy cooked one
 630                                 if(gAlice->GetMCApp()->GetCurrentTrackNumber()==fPrimaryTrackPid) {
 631 //                                      printf("primary pid=%d eDelta=%f\n",pid,eDelta);
 632                                         if(pid==27){
 633                                                 if(eDelta>=kECut){
 634                                                         pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
 635                                 if(pr_range>=(3.7-range)) eDelta*=.1;
 636                                                 }
 637                                         } else if(pid==1){
 638                                                 if(eDelta<kECut) eDelta*=.5;
 639                                                 else {
 640                                                         pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
 641                                 if(pr_range>=((AliTRDgeometry::DrThick()
 642                        + AliTRDgeometry::AmThick())-range)) eDelta*=.05;
 643                                                         else eDelta*=.5;
 644                                                 }
 645                                         } else eDelta=0.;
 646                                 } else eDelta=0.;
 647
 648         // Generate the electron cluster size
 649         if(eDelta==0.) qTot=0;
 650                                 else qTot = ((Int_t) (eDelta / kWion) + 1);
 651                                 // Create a new dEdx hit
 652         AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot, drRegion);
 653
 654         // Calculate the maximum step size for the next tracking step
 655         // Produce only one hit if Ekin is below cutoff
 656         aMass = gMC->TrackMass();
 657         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 658
 659           // The energy loss according to Bethe Bloch
 660           iPdg  = TMath::Abs(gMC->TrackPid());
 661           if ( (iPdg != kPdgElectron) ||
 662                                 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
 663             gMC->TrackMomentum(mom);
 664             pTot      = mom.Rho();
 665             betaGamma = pTot / aMass;
 666             pp        = BetheBlochGeant(betaGamma);
 667                         // Take charge > 1 into account
 668             charge = gMC->TrackCharge();
 669             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 670           } else { // Electrons above 20 Mev/c are at the plateau
 671                                                 pp = kPrim * kPlateau;
 672           }
 673
 674           Int_t nsteps = 0;
 675           do {nsteps = gRandom->Poisson(pp);} while(!nsteps);
 676           stepSize = 1./nsteps;
 677           gMC->SetMaxStep(stepSize);
 678         }
 679       }
 680     }
 681   }
 682 }
 683
 684 //_____________________________________________________________________________
 685 void AliTRDv1::StepManagerErmilova()
 686 {
 687   //
 688   // Slow simulator. Every charged track produces electron cluster as hits
 689   // along its path across the drift volume. The step size is set acording
 690   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 691   // a spectrum taken from Ermilova et al.
 692   //
 693
 694   Int_t    pla = 0;
 695   Int_t    cha = 0;
 696   Int_t    sec = 0;
 697   Int_t    det = 0;
 698   Int_t    iPdg;
 699   Int_t    qTot;
 700
 701   Float_t  hits[3];
 702   Double_t random[1];
 703   Float_t  charge;
 704   Float_t  aMass;
 705
 706   Double_t pTot = 0;
 707   Double_t eDelta;
 708   Double_t betaGamma, pp;
 709   Double_t stepSize;
 710
 711   Bool_t   drRegion = kFALSE;
 712   Bool_t   amRegion = kFALSE;
 713
 714   TString  cIdCurrent;
 715   TString  cIdSensDr = "J";
 716   TString  cIdSensAm = "K";
 717   Char_t   cIdChamber[3];
 718            cIdChamber[2] = 0;
 719
 720   TLorentzVector pos, mom;
 721
 722   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 723   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 724   const Int_t    kNdetsec     = kNplan * kNcham;
 725
 726   const Double_t kBig         = 1.0E+12; // Infinitely big
 727   const Float_t  kWion        = 23.53;   // Ionization energy
 728   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 729
 730   // energy threshold for production of delta electrons
 731   //const Float_t  kECut = 1.0e4;
 732   // Parameters entering the parametrized range for delta electrons
 733   //const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
 734   // Gas density -> To be made user adjustable !
 735   //const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
 736
 737   // Minimum energy for the step size adjustment
 738   const Float_t  kEkinMinStep = 1.0e-5;
 739
 740   // Plateau value of the energy-loss for electron in xenon
 741   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 742   //const Double_t kPlateau = 1.70;
 743   // the averaged value (26/3/99)
 744   const Float_t  kPlateau     = 1.55;
 745
 746   const Float_t  kPrim        = 48.0;  // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 747   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 748   const Float_t  kPoti        = 12.1;
 749
 750   const Int_t    kPdgElectron = 11;  // PDG code electron
 751
 752   // Set the maximum step size to a very large number for all
 753   // neutral particles and those outside the driftvolume
 754   gMC->SetMaxStep(kBig);
 755
 756   // Use only charged tracks
 757   if (( gMC->TrackCharge()       ) &&
 758       (!gMC->IsTrackStop()       ) &&
 759       (!gMC->IsTrackDisappeared())) {
 760
 761     // Inside a sensitive volume?
 762     drRegion = kFALSE;
 763     amRegion = kFALSE;
 764     cIdCurrent = gMC->CurrentVolName();
 765     if (cIdSensDr == cIdCurrent[1]) {
 766       drRegion = kTRUE;
 767     }
 768     if (cIdSensAm == cIdCurrent[1]) {
 769       amRegion = kTRUE;
 770     }
 771     if (drRegion || amRegion) {
 772
 773       // The hit coordinates and charge
 774       gMC->TrackPosition(pos);
 775       hits[0] = pos[0];
 776       hits[1] = pos[1];
 777       hits[2] = pos[2];
 778
 779       // The sector number (0 - 17)
 780       // The numbering goes clockwise and starts at y = 0
 781       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 782       if (phi < 90.)
 783         phi = phi + 270.;
 784       else
 785         phi = phi -  90.;
 786       sec = ((Int_t) (phi / 20));
 787
 788       // The plane and chamber number
 789       cIdChamber[0] = cIdCurrent[2];
 790       cIdChamber[1] = cIdCurrent[3];
 791       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 792       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 793       pla = ((Int_t) idChamber % kNplan);
 794
 795       // Check on selected volumes
 796       Int_t addthishit = 1;
 797       if (fSensSelect) {
 798         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 799         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 800         if (fSensSector  >= 0) {
 801           Int_t sens1  = fSensSector;
 802           Int_t sens2  = fSensSector + fSensSectorRange;
 803                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 804                        * AliTRDgeometry::Nsect();
 805           if (sens1 < sens2) {
 806             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 807                                 }
 808           else {
 809             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 810                                 }
 811                                 }
 812       }
 813
 814       // Add this hit
 815       if (addthishit) {
 816
 817         // The detector number
 818         det = fGeometry->GetDetector(pla,cha,sec);
 819
 820         // Special hits only in the drift region
 821         if (drRegion) {
 822
 823           // Create a track reference at the entrance and
 824           // exit of each chamber that contain the
 825           // momentum components of the particle
 826           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
 827             gMC->TrackMomentum(mom);
 828             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 829           }
 830           // Create the hits from TR photons
 831           if (fTR) CreateTRhit(det);
 832                                 }
 833
 834
 835         // Calculate the energy of the delta-electrons
 836         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 837         eDelta = TMath::Max(eDelta,0.0);
 838         // Generate the electron cluster size
 839         if(eDelta==0.) qTot=0;
 840                                 else qTot = ((Int_t) (eDelta / kWion) + 1);
 841
 842                                 // Create a new dEdx hit
 843         if (drRegion) {
 844           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 845                 ,det,hits,qTot, kTRUE);
 846                                 }
 847         else {
 848           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 849                 ,det,hits,qTot,kFALSE);
 850                                 }
 851
 852         // Calculate the maximum step size for the next tracking step
 853         // Produce only one hit if Ekin is below cutoff
 854         aMass = gMC->TrackMass();
 855         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 856
 857           // The energy loss according to Bethe Bloch
 858           iPdg  = TMath::Abs(gMC->TrackPid());
 859           if ( (iPdg != kPdgElectron) ||
 860                                 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
 861             gMC->TrackMomentum(mom);
 862             pTot      = mom.Rho();
 863             betaGamma = pTot / aMass;
 864             pp        = kPrim * BetheBloch(betaGamma);
 865             // Take charge > 1 into account
 866             charge = gMC->TrackCharge();
 867             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 868           } else { // Electrons above 20 Mev/c are at the plateau
 869                                                 pp = kPrim * kPlateau;
 870           }
 871
 872           if (pp > 0) {
 873             do
 874             gMC->GetRandom()->RndmArray(1, random);
 875             while ((random[0] == 1.) || (random[0] == 0.));
 876             stepSize = - TMath::Log(random[0]) / pp;
 877             gMC->SetMaxStep(stepSize);
 878                                 }
 879                                 }
 880       }
 881     }
 882   }
 883 }
 884
 885 //_____________________________________________________________________________
 886 void AliTRDv1::StepManagerFixedStep()
 887 {
 888   //
 889   // Slow simulator. Every charged track produces electron cluster as hits
 890   // along its path across the drift volume. The step size is fixed in
 891   // this version of the step manager.
 892   //
 893
 894   Int_t    pla = 0;
 895   Int_t    cha = 0;
 896   Int_t    sec = 0;
 897   Int_t    det = 0;
 898   Int_t    qTot;
 899
 900   Float_t  hits[3];
 901   Double_t eDep;
 902
 903   Bool_t   drRegion = kFALSE;
 904   Bool_t   amRegion = kFALSE;
 905
 906   TString  cIdCurrent;
 907   TString  cIdSensDr = "J";
 908   TString  cIdSensAm = "K";
 909   Char_t   cIdChamber[3];
 910   cIdChamber[2] = 0;
 911
 912   TLorentzVector pos, mom;
 913
 914   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 915   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 916   const Int_t    kNdetsec     = kNplan * kNcham;
 917
 918   const Double_t kBig         = 1.0E+12;
 919
 920   const Float_t  kWion        = 23.53;   // Ionization energy
 921   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
 922
 923   // Set the maximum step size to a very large number for all
 924   // neutral particles and those outside the driftvolume
 925   gMC->SetMaxStep(kBig);
 926
 927   // If not charged track or already stopped or disappeared, just return.
 928   if ((!gMC->TrackCharge()) ||
 929         gMC->IsTrackStop()  ||
 930         gMC->IsTrackDisappeared()) return;
 931
 932   // Inside a sensitive volume?
 933   cIdCurrent = gMC->CurrentVolName();
 934
 935   if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
 936   if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
 937
 938   if ((!drRegion) && (!amRegion)) return;
 939
 940   // The hit coordinates and charge
 941   gMC->TrackPosition(pos);
 942   hits[0] = pos[0];
 943   hits[1] = pos[1];
 944   hits[2] = pos[2];
 945
 946   // The sector number (0 - 17)
 947   // The numbering goes clockwise and starts at y = 0
 948   Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 949   if (phi < 90.) phi += 270.;
 950   else           phi -=  90.;
 951   sec = ((Int_t) (phi / 20.));
 952
 953   // The plane and chamber number
 954   cIdChamber[0] = cIdCurrent[2];
 955   cIdChamber[1] = cIdCurrent[3];
 956   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 957   cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 958   pla = ((Int_t) idChamber % kNplan);
 959
 960   // Check on selected volumes
 961   Int_t addthishit = 1;
 962   if(fSensSelect) {
 963     if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 964     if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 965     if (fSensSector  >= 0) {
 966       Int_t sens1  = fSensSector;
 967       Int_t sens2  = fSensSector + fSensSectorRange;
 968       sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
 969       if (sens1 < sens2) {
 970         if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 971       }
 972       else {
 973         if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 974       }
 975     }
 976   }
 977
 978   if (!addthishit) return;
 979
 980   det = fGeometry->GetDetector(pla,cha,sec);  // The detector number
 981
 982   Int_t trkStat = 0;  // 0: InFlight 1:Entering 2:Exiting
 983
 984   // Special hits only in the drift region
 985   if (drRegion) {
 986
 987     // Create a track reference at the entrance and exit of each
 988     // chamber that contain the momentum components of the particle
 989
 990     if (gMC->IsTrackEntering()) {
 991       gMC->TrackMomentum(mom);
 992       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 993       trkStat = 1;
 994     }
 995     if (gMC->IsTrackExiting()) {
 996       gMC->TrackMomentum(mom);
 997       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 998       trkStat = 2;
 999     }
1000
1001     // Create the hits from TR photons
1002     if (fTR) CreateTRhit(det);
1003
1004   }
1005
1006   // Calculate the charge according to GEANT Edep
1007   // Create a new dEdx hit
1008   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1009   qTot = (Int_t) (eDep / kWion);
1010   AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1011         ,det,hits,qTot,drRegion);
1012
1013   // Set Maximum Step Size
1014   // Produce only one hit if Ekin is below cutoff
1015   if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
1016   gMC->SetMaxStep(fStepSize);
1017
1018 }
1019
1020 //_____________________________________________________________________________
1021 Double_t AliTRDv1::BetheBloch(Double_t bg)
1022 {
1023   //
1024   // Parametrization of the Bethe-Bloch-curve
1025   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1026   //
1027
1028   // This parameters have been adjusted to averaged values from GEANT
1029   const Double_t kP1 = 7.17960e-02;
1030   const Double_t kP2 = 8.54196;
1031   const Double_t kP3 = 1.38065e-06;
1032   const Double_t kP4 = 5.30972;
1033   const Double_t kP5 = 2.83798;
1034
1035   // This parameters have been adjusted to Xe-data found in:
1036   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
1037   //const Double_t kP1 = 0.76176E-1;
1038   //const Double_t kP2 = 10.632;
1039   //const Double_t kP3 = 3.17983E-6;
1040   //const Double_t kP4 = 1.8631;
1041   //const Double_t kP5 = 1.9479;
1042
1043   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1044   const Double_t kBgMin = 0.8;
1045   const Double_t kBBMax = 6.83298;
1046   //const Double_t kBgMin = 0.6;
1047   //const Double_t kBBMax = 17.2809;
1048   //const Double_t kBgMin = 0.4;
1049   //const Double_t kBBMax = 82.0;
1050
1051   if (bg > kBgMin) {
1052     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
1053     Double_t aa = TMath::Power(yy,kP4);
1054     Double_t bb = TMath::Power((1./bg),kP5);
1055              bb = TMath::Log(kP3 + bb);
1056     return ((kP2 - aa - bb)*kP1 / aa);
1057   }
1058   else {
1059     return kBBMax;
1060   }
1061
1062 }
1063
1064 //_____________________________________________________________________________
1065 Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
1066 {
1067   //
1068   // Return dN/dx (number of primary collisions per centimeter)
1069   // for given beta*gamma factor.
1070   //
1071   // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1072   // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1073   // This must be used as a set with IntSpecGeant.
1074   //
1075
1076   Double_t arr_g[20] = {
1077     1.100000,   1.200000,    1.300000,    1.500000,
1078     1.800000,   2.000000,    2.500000,    3.000000,
1079     4.000000,   7.000000,    10.000000,   20.000000,
1080     40.000000,  70.000000,   100.000000,  300.000000,
1081     600.000000, 1000.000000, 3000.000000, 10000.000000 };
1082
1083   Double_t arr_nc[20] = {
1084     75.009056,   45.508083,   35.299252,   27.116327,
1085     22.734999,   21.411915,   19.934095,   19.449375,
1086     19.344431,   20.185553,   21.027925,   22.912676,
1087     24.933352,   26.504053,   27.387468,   29.566597,
1088     30.353779,   30.787134,   31.129285,   31.157350 };
1089
1090   // betagamma to gamma
1091   Double_t g = TMath::Sqrt( 1. + bg*bg );
1092
1093   // Find the index just before the point we need.
1094   int i;
1095   for( i = 0 ; i < 18 ; i++ )
1096     if( arr_g[i] < g && arr_g[i+1] > g )
1097       break;
1098
1099   // Simple interpolation.
1100   Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
1101                  (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
1102
1103   return pp; //arr_nc[8];
1104
1105 }
1106
1107 //_____________________________________________________________________________
1108 void AliTRDv1::Stepping()
1109 {
1110 // Stepping info
1111 // ---
1112
1113    cout << "X(cm)    "
1114        << "Y(cm)    "
1115        << "Z(cm)  "
1116        << "KinE(MeV)   "
1117        << "dE(MeV) "
1118        << "Step(cm) "
1119        << "TrackL(cm) "
1120        << "Volume  "
1121        << "Process "
1122        << endl;
1123
1124    // Position
1125     //
1126     Double_t x, y, z;
1127     gMC->TrackPosition(x, y, z);
1128     cout << setw(8) << setprecision(3) << x << " "
1129          << setw(8) << setprecision(3) << y << " "
1130          << setw(8) << setprecision(3) << z << "  ";
1131
1132     // Kinetic energy
1133     //
1134     Double_t px, py, pz, etot;
1135     gMC->TrackMomentum(px, py, pz, etot);
1136     Double_t ekin = etot - gMC->TrackMass();
1137     cout << setw(9) << setprecision(4) << ekin*1e03 << " ";
1138
1139     // Energy deposit
1140     //
1141     cout << setw(9) << setprecision(4) << gMC->Edep()*1e03 << " ";
1142
1143     // Step length
1144     //
1145     cout << setw(8) << setprecision(3) << gMC->TrackStep() << " ";
1146
1147     // Track length
1148     //
1149     cout << setw(8) << setprecision(3) << gMC->TrackLength() << "     ";
1150
1151     // Volume
1152     //
1153     if (gMC->CurrentVolName() != 0)
1154       cout << setw(4) << gMC->CurrentVolName() << "  ";
1155     else
1156       cout << setw(4) << "None"  << "  ";
1157
1158     // Process
1159     //
1160     TArrayI processes;
1161     Int_t nofProcesses = gMC->StepProcesses(processes);
1162     for(int ip=0;ip<nofProcesses; ip++)
1163       cout << TMCProcessName[processes[ip]]<<" / ";
1164
1165     cout << endl;
1166 }
1167
1168
1169 //_____________________________________________________________________________
1170 Double_t Ermilova(Double_t *x, Double_t *)
1171 {
1172   //
1173   // Calculates the delta-ray energy distribution according to Ermilova.
1174   // Logarithmic scale !
1175   //
1176
1177   Double_t energy;
1178   Double_t dpos;
1179   Double_t dnde;
1180
1181   Int_t    pos1, pos2;
1182
1183   const Int_t kNv = 31;
1184
1185   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
1186                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
1187                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
1188                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
1189                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
1190                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
1191                      ,11.5129 };
1192
1193   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
1194                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
1195                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
1196                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
1197                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
1198                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
1199                      ,  0.004 };
1200
1201   energy = x[0];
1202
1203   // Find the position
1204   pos1 = pos2 = 0;
1205   dpos = 0;
1206   do {
1207     dpos = energy - vxe[pos2++];
1208   }
1209   while (dpos > 0);
1210   pos2--;
1211   if (pos2 > kNv) pos2 = kNv - 1;
1212   pos1 = pos2 - 1;
1213
1214   // Differentiate between the sampling points
1215   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1216
1217   return dnde;
1218
1219 }
1220
1221 //_____________________________________________________________________________
1222 Double_t IntSpecGeant(Double_t *x, Double_t *)
1223 {
1224   //
1225   // Integrated spectrum from Geant3
1226   //
1227
1228   const Int_t npts = 83;
1229   Double_t arre[npts] = {
1230     2.421257,  2.483278,  2.534301,  2.592230,
1231     2.672067,  2.813299,  3.015059,  3.216819,
1232     3.418579,  3.620338,  3.868209,  3.920198,
1233     3.978284,  4.063923,  4.186264,  4.308605,
1234     4.430946,  4.553288,  4.724261,  4.837736,
1235     4.999842,  5.161949,  5.324056,  5.486163,
1236     5.679688,  5.752998,  5.857728,  5.962457,
1237     6.067185,  6.171914,  6.315653,  6.393674,
1238     6.471694,  6.539689,  6.597658,  6.655627,
1239     6.710957,  6.763648,  6.816338,  6.876198,
1240     6.943227,  7.010257,  7.106285,  7.252151,
1241     7.460531,  7.668911,  7.877290,  8.085670,
1242     8.302979,  8.353585,  8.413120,  8.483500,
1243     8.541030,  8.592857,  8.668865,  8.820485,
1244     9.037086,  9.253686,  9.470286,  9.686887,
1245     9.930838,  9.994655, 10.085822, 10.176990,
1246     10.268158, 10.359325, 10.503614, 10.627565,
1247     10.804637, 10.981709, 11.158781, 11.335854,
1248     11.593397, 11.781165, 12.049404, 12.317644,
1249     12.585884, 12.854123, 14.278421, 16.975889,
1250     20.829416, 24.682943, 28.536469
1251   };
1252   /*
1253   Double_t arrdndx[npts] = {
1254     19.344431, 18.664679, 18.136106, 17.567745,
1255     16.836426, 15.677382, 14.281277, 13.140237,
1256     12.207677, 11.445510, 10.697049, 10.562296,
1257     10.414673, 10.182341,  9.775256,  9.172330,
1258     8.240271,  6.898587,  4.808303,  3.889751,
1259     3.345288,  3.093431,  2.897347,  2.692470,
1260     2.436222,  2.340029,  2.208579,  2.086489,
1261     1.975535,  1.876519,  1.759626,  1.705024,
1262     1.656374,  1.502638,  1.330566,  1.200697,
1263     1.101168,  1.019323,  0.943867,  0.851951,
1264     0.755229,  0.671576,  0.570675,  0.449672,
1265     0.326722,  0.244225,  0.188225,  0.149608,
1266     0.121529,  0.116289,  0.110636,  0.103490,
1267     0.096147,  0.089191,  0.079780,  0.063927,
1268     0.047642,  0.036341,  0.028250,  0.022285,
1269     0.017291,  0.016211,  0.014802,  0.013533,
1270     0.012388,  0.011352,  0.009803,  0.008537,
1271     0.007039,  0.005829,  0.004843,  0.004034,
1272     0.003101,  0.002564,  0.001956,  0.001494,
1273     0.001142,  0.000873,  0.000210,  0.000014,
1274     0.000000,  0.000000,  0.000000
1275   };
1276   */
1277   // Differentiate
1278   //  dnde = (arrdndx[i-1] - arrdndx[i]) / (arre[i] - arre[i-1]);
1279
1280   Double_t arrdnde[npts] = {
1281     10.960000, 10.960000, 10.359500,  9.811340,
1282     9.1601500,  8.206670,  6.919630,  5.655430,
1283     4.6221300,  3.777610,  3.019560,  2.591950,
1284     2.5414600,  2.712920,  3.327460,  4.928240,
1285     7.6185300, 10.966700, 12.225800,  8.094750,
1286     3.3586900,  1.553650,  1.209600,  1.263840,
1287     1.3241100,  1.312140,  1.255130,  1.165770,
1288     1.0594500,  0.945450,  0.813231,  0.699837,
1289     0.6235580,  2.260990,  2.968350,  2.240320,
1290     1.7988300,  1.553300,  1.432070,  1.535520,
1291     1.4429900,  1.247990,  1.050750,  0.829549,
1292     0.5900280,  0.395897,  0.268741,  0.185320,
1293     0.1292120,  0.103545,  0.0949525, 0.101535,
1294     0.1276380,  0.134216,  0.123816,  0.104557,
1295     0.0751843,  0.0521745, 0.0373546, 0.0275391,
1296     0.0204713,  0.0169234, 0.0154552, 0.0139194,
1297     0.0125592,  0.0113638, 0.0107354, 0.0102137,
1298     0.00845984, 0.00683338, 0.00556836, 0.00456874,
1299     0.0036227,  0.00285991, 0.00226664, 0.00172234,
1300     0.00131226, 0.00100284, 0.000465492, 7.26607e-05,
1301     3.63304e-06, 0.0000000, 0.0000000
1302   };
1303
1304   Int_t i;
1305   Double_t energy = x[0];
1306
1307   for( i = 0 ; i < npts ; i++ )
1308     if( energy < arre[i] ) break;
1309
1310   if( i == 0 )
1311     AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");
1312
1313    return arrdnde[i];
1314
1315 }