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