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        = 22.57;
 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        = 22.57;   // 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           stepSize = 1./gRandom->Poisson(pp);
 679                                         gMC->SetMaxStep(stepSize);
 680                                 }
 681       }
 682     }
 683   }
 684 }
 685
 686 //_____________________________________________________________________________
 687 void AliTRDv1::StepManagerErmilova()
 688 {
 689   //
 690   // Slow simulator. Every charged track produces electron cluster as hits
 691   // along its path across the drift volume. The step size is set acording
 692   // to Bethe-Bloch. The energy distribution of the delta electrons follows
 693   // a spectrum taken from Ermilova et al.
 694   //
 695
 696   Int_t    pla = 0;
 697   Int_t    cha = 0;
 698   Int_t    sec = 0;
 699   Int_t    det = 0;
 700   Int_t    iPdg;
 701   Int_t    qTot;
 702
 703   Float_t  hits[3];
 704   Double_t random[1];
 705   Float_t  charge;
 706   Float_t  aMass;
 707
 708   Double_t pTot = 0;
 709   Double_t eDelta;
 710   Double_t betaGamma, pp;
 711   Double_t stepSize;
 712
 713   Bool_t   drRegion = kFALSE;
 714   Bool_t   amRegion = kFALSE;
 715
 716   TString  cIdCurrent;
 717   TString  cIdSensDr = "J";
 718   TString  cIdSensAm = "K";
 719   Char_t   cIdChamber[3];
 720            cIdChamber[2] = 0;
 721
 722   TLorentzVector pos, mom;
 723
 724   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 725   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 726   const Int_t    kNdetsec     = kNplan * kNcham;
 727
 728   const Double_t kBig         = 1.0E+12; // Infinitely big
 729   const Float_t  kWion        = 22.57;   // Ionization energy
 730   const Float_t  kPTotMaxEl   = 0.002;   // Maximum momentum for e+ e- g
 731
 732   // energy threshold for production of delta electrons
 733   //const Float_t  kECut = 1.0e4;
 734   // Parameters entering the parametrized range for delta electrons
 735   //const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
 736   // Gas density -> To be made user adjustable !
 737   //const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
 738
 739   // Minimum energy for the step size adjustment
 740   const Float_t  kEkinMinStep = 1.0e-5;
 741
 742   // Plateau value of the energy-loss for electron in xenon
 743   // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
 744   //const Double_t kPlateau = 1.70;
 745   // the averaged value (26/3/99)
 746   const Float_t  kPlateau     = 1.55;
 747
 748   const Float_t  kPrim        = 48.0;  // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
 749   // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
 750   const Float_t  kPoti        = 12.1;
 751
 752   const Int_t    kPdgElectron = 11;  // PDG code electron
 753
 754   // Set the maximum step size to a very large number for all
 755   // neutral particles and those outside the driftvolume
 756   gMC->SetMaxStep(kBig);
 757
 758   // Use only charged tracks
 759   if (( gMC->TrackCharge()       ) &&
 760       (!gMC->IsTrackStop()       ) &&
 761       (!gMC->IsTrackDisappeared())) {
 762
 763     // Inside a sensitive volume?
 764     drRegion = kFALSE;
 765     amRegion = kFALSE;
 766     cIdCurrent = gMC->CurrentVolName();
 767     if (cIdSensDr == cIdCurrent[1]) {
 768       drRegion = kTRUE;
 769     }
 770     if (cIdSensAm == cIdCurrent[1]) {
 771       amRegion = kTRUE;
 772     }
 773     if (drRegion || amRegion) {
 774
 775       // The hit coordinates and charge
 776       gMC->TrackPosition(pos);
 777       hits[0] = pos[0];
 778       hits[1] = pos[1];
 779       hits[2] = pos[2];
 780
 781       // The sector number (0 - 17)
 782       // The numbering goes clockwise and starts at y = 0
 783       Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 784       if (phi < 90.)
 785         phi = phi + 270.;
 786       else
 787         phi = phi -  90.;
 788       sec = ((Int_t) (phi / 20));
 789
 790       // The plane and chamber number
 791       cIdChamber[0] = cIdCurrent[2];
 792       cIdChamber[1] = cIdCurrent[3];
 793       Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 794       cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 795       pla = ((Int_t) idChamber % kNplan);
 796
 797       // Check on selected volumes
 798       Int_t addthishit = 1;
 799       if (fSensSelect) {
 800         if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 801         if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 802         if (fSensSector  >= 0) {
 803           Int_t sens1  = fSensSector;
 804           Int_t sens2  = fSensSector + fSensSectorRange;
 805                 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
 806                        * AliTRDgeometry::Nsect();
 807           if (sens1 < sens2) {
 808             if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 809                                 }
 810           else {
 811             if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 812                                 }
 813                                 }
 814       }
 815
 816       // Add this hit
 817       if (addthishit) {
 818
 819         // The detector number
 820         det = fGeometry->GetDetector(pla,cha,sec);
 821
 822         // Special hits only in the drift region
 823         if (drRegion) {
 824
 825           // Create a track reference at the entrance and
 826           // exit of each chamber that contain the
 827           // momentum components of the particle
 828           if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
 829             gMC->TrackMomentum(mom);
 830             AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 831           }
 832           // Create the hits from TR photons
 833           if (fTR) CreateTRhit(det);
 834                                 }
 835
 836
 837         // Calculate the energy of the delta-electrons
 838         eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
 839         eDelta = TMath::Max(eDelta,0.0);
 840         // Generate the electron cluster size
 841         if(eDelta==0.) qTot=0;
 842                                 else qTot = ((Int_t) (eDelta / kWion) + 1);
 843
 844                                 // Create a new dEdx hit
 845         if (drRegion) {
 846           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 847                 ,det,hits,qTot, kTRUE);
 848                                 }
 849         else {
 850           AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
 851                 ,det,hits,qTot,kFALSE);
 852                                 }
 853
 854         // Calculate the maximum step size for the next tracking step
 855         // Produce only one hit if Ekin is below cutoff
 856         aMass = gMC->TrackMass();
 857         if ((gMC->Etot() - aMass) > kEkinMinStep) {
 858
 859           // The energy loss according to Bethe Bloch
 860           iPdg  = TMath::Abs(gMC->TrackPid());
 861           if ( (iPdg != kPdgElectron) ||
 862                                 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
 863             gMC->TrackMomentum(mom);
 864             pTot      = mom.Rho();
 865             betaGamma = pTot / aMass;
 866             pp        = kPrim * BetheBloch(betaGamma);
 867             // Take charge > 1 into account
 868             charge = gMC->TrackCharge();
 869             if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
 870           } else { // Electrons above 20 Mev/c are at the plateau
 871                                                 pp = kPrim * kPlateau;
 872           }
 873
 874           if (pp > 0) {
 875             do
 876             gMC->GetRandom()->RndmArray(1, random);
 877             while ((random[0] == 1.) || (random[0] == 0.));
 878             stepSize = - TMath::Log(random[0]) / pp;
 879             gMC->SetMaxStep(stepSize);
 880                                 }
 881                                 }
 882       }
 883     }
 884   }
 885 }
 886
 887 //_____________________________________________________________________________
 888 void AliTRDv1::StepManagerFixedStep()
 889 {
 890   //
 891   // Slow simulator. Every charged track produces electron cluster as hits
 892   // along its path across the drift volume. The step size is fixed in
 893   // this version of the step manager.
 894   //
 895
 896   Int_t    pla = 0;
 897   Int_t    cha = 0;
 898   Int_t    sec = 0;
 899   Int_t    det = 0;
 900   Int_t    qTot;
 901
 902   Float_t  hits[3];
 903   Double_t eDep;
 904
 905   Bool_t   drRegion = kFALSE;
 906   Bool_t   amRegion = kFALSE;
 907
 908   TString  cIdCurrent;
 909   TString  cIdSensDr = "J";
 910   TString  cIdSensAm = "K";
 911   Char_t   cIdChamber[3];
 912   cIdChamber[2] = 0;
 913
 914   TLorentzVector pos, mom;
 915
 916   const Int_t    kNplan       = AliTRDgeometry::Nplan();
 917   const Int_t    kNcham       = AliTRDgeometry::Ncham();
 918   const Int_t    kNdetsec     = kNplan * kNcham;
 919
 920   const Double_t kBig         = 1.0E+12;
 921
 922   const Float_t  kWion        = 22.57;   // Ionization energy
 923   const Float_t  kEkinMinStep = 1.0e-5;  // Minimum energy for the step size adjustment
 924
 925   // Set the maximum step size to a very large number for all
 926   // neutral particles and those outside the driftvolume
 927   gMC->SetMaxStep(kBig);
 928
 929   // If not charged track or already stopped or disappeared, just return.
 930   if ((!gMC->TrackCharge()) ||
 931         gMC->IsTrackStop()  ||
 932         gMC->IsTrackDisappeared()) return;
 933
 934   // Inside a sensitive volume?
 935   cIdCurrent = gMC->CurrentVolName();
 936
 937   if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
 938   if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
 939
 940   if ((!drRegion) && (!amRegion)) return;
 941
 942   // The hit coordinates and charge
 943   gMC->TrackPosition(pos);
 944   hits[0] = pos[0];
 945   hits[1] = pos[1];
 946   hits[2] = pos[2];
 947
 948   // The sector number (0 - 17)
 949   // The numbering goes clockwise and starts at y = 0
 950   Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
 951   if (phi < 90.) phi += 270.;
 952   else           phi -=  90.;
 953   sec = ((Int_t) (phi / 20.));
 954
 955   // The plane and chamber number
 956   cIdChamber[0] = cIdCurrent[2];
 957   cIdChamber[1] = cIdCurrent[3];
 958   Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
 959   cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
 960   pla = ((Int_t) idChamber % kNplan);
 961
 962   // Check on selected volumes
 963   Int_t addthishit = 1;
 964   if(fSensSelect) {
 965     if ((fSensPlane   >= 0) && (pla != fSensPlane  )) addthishit = 0;
 966     if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
 967     if (fSensSector  >= 0) {
 968       Int_t sens1  = fSensSector;
 969       Int_t sens2  = fSensSector + fSensSectorRange;
 970       sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
 971       if (sens1 < sens2) {
 972         if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
 973       }
 974       else {
 975         if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
 976       }
 977     }
 978   }
 979
 980   if (!addthishit) return;
 981
 982   det = fGeometry->GetDetector(pla,cha,sec);  // The detector number
 983
 984   Int_t trkStat = 0;  // 0: InFlight 1:Entering 2:Exiting
 985
 986   // Special hits only in the drift region
 987   if (drRegion) {
 988
 989     // Create a track reference at the entrance and exit of each
 990     // chamber that contain the momentum components of the particle
 991
 992     if (gMC->IsTrackEntering()) {
 993       gMC->TrackMomentum(mom);
 994       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
 995       trkStat = 1;
 996     }
 997     if (gMC->IsTrackExiting()) {
 998       gMC->TrackMomentum(mom);
 999       AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1000       trkStat = 2;
1001     }
1002
1003     // Create the hits from TR photons
1004     if (fTR) CreateTRhit(det);
1005
1006   }
1007
1008   // Calculate the charge according to GEANT Edep
1009   // Create a new dEdx hit
1010   eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1011   qTot = (Int_t) (eDep / kWion);
1012   AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1013         ,det,hits,qTot,drRegion);
1014
1015   // Set Maximum Step Size
1016   // Produce only one hit if Ekin is below cutoff
1017   if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
1018   gMC->SetMaxStep(fStepSize);
1019
1020 }
1021
1022 //_____________________________________________________________________________
1023 Double_t AliTRDv1::BetheBloch(Double_t bg)
1024 {
1025   //
1026   // Parametrization of the Bethe-Bloch-curve
1027   // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1028   //
1029
1030   // This parameters have been adjusted to averaged values from GEANT
1031   const Double_t kP1 = 7.17960e-02;
1032   const Double_t kP2 = 8.54196;
1033   const Double_t kP3 = 1.38065e-06;
1034   const Double_t kP4 = 5.30972;
1035   const Double_t kP5 = 2.83798;
1036
1037   // This parameters have been adjusted to Xe-data found in:
1038   // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
1039   //const Double_t kP1 = 0.76176E-1;
1040   //const Double_t kP2 = 10.632;
1041   //const Double_t kP3 = 3.17983E-6;
1042   //const Double_t kP4 = 1.8631;
1043   //const Double_t kP5 = 1.9479;
1044
1045   // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1046   const Double_t kBgMin = 0.8;
1047   const Double_t kBBMax = 6.83298;
1048   //const Double_t kBgMin = 0.6;
1049   //const Double_t kBBMax = 17.2809;
1050   //const Double_t kBgMin = 0.4;
1051   //const Double_t kBBMax = 82.0;
1052
1053   if (bg > kBgMin) {
1054     Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
1055     Double_t aa = TMath::Power(yy,kP4);
1056     Double_t bb = TMath::Power((1./bg),kP5);
1057              bb = TMath::Log(kP3 + bb);
1058     return ((kP2 - aa - bb)*kP1 / aa);
1059   }
1060   else {
1061     return kBBMax;
1062   }
1063
1064 }
1065
1066 //_____________________________________________________________________________
1067 Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
1068 {
1069   //
1070   // Return dN/dx (number of primary collisions per centimeter)
1071   // for given beta*gamma factor.
1072   //
1073   // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1074   // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1075   // This must be used as a set with IntSpecGeant.
1076   //
1077
1078   Double_t arr_g[20] = {
1079     1.100000,   1.200000,    1.300000,    1.500000,
1080     1.800000,   2.000000,    2.500000,    3.000000,
1081     4.000000,   7.000000,    10.000000,   20.000000,
1082     40.000000,  70.000000,   100.000000,  300.000000,
1083     600.000000, 1000.000000, 3000.000000, 10000.000000 };
1084
1085   Double_t arr_nc[20] = {
1086     75.009056,   45.508083,   35.299252,   27.116327,
1087     22.734999,   21.411915,   19.934095,   19.449375,
1088     19.344431,   20.185553,   21.027925,   22.912676,
1089     24.933352,   26.504053,   27.387468,   29.566597,
1090     30.353779,   30.787134,   31.129285,   31.157350 };
1091
1092   // betagamma to gamma
1093   Double_t g = TMath::Sqrt( 1. + bg*bg );
1094
1095   // Find the index just before the point we need.
1096   int i;
1097   for( i = 0 ; i < 18 ; i++ )
1098     if( arr_g[i] < g && arr_g[i+1] > g )
1099       break;
1100
1101   // Simple interpolation.
1102   Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
1103                  (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
1104
1105   return pp; //arr_nc[8];
1106
1107 }
1108
1109 //_____________________________________________________________________________
1110 void AliTRDv1::Stepping()
1111 {
1112 // Stepping info
1113 // ---
1114
1115    cout << "X(cm)    "
1116        << "Y(cm)    "
1117        << "Z(cm)  "
1118        << "KinE(MeV)   "
1119        << "dE(MeV) "
1120        << "Step(cm) "
1121        << "TrackL(cm) "
1122        << "Volume  "
1123        << "Process "
1124        << endl;
1125
1126    // Position
1127     //
1128     Double_t x, y, z;
1129     gMC->TrackPosition(x, y, z);
1130     cout << setw(8) << setprecision(3) << x << " "
1131          << setw(8) << setprecision(3) << y << " "
1132          << setw(8) << setprecision(3) << z << "  ";
1133
1134     // Kinetic energy
1135     //
1136     Double_t px, py, pz, etot;
1137     gMC->TrackMomentum(px, py, pz, etot);
1138     Double_t ekin = etot - gMC->TrackMass();
1139     cout << setw(9) << setprecision(4) << ekin*1e03 << " ";
1140
1141     // Energy deposit
1142     //
1143     cout << setw(9) << setprecision(4) << gMC->Edep()*1e03 << " ";
1144
1145     // Step length
1146     //
1147     cout << setw(8) << setprecision(3) << gMC->TrackStep() << " ";
1148
1149     // Track length
1150     //
1151     cout << setw(8) << setprecision(3) << gMC->TrackLength() << "     ";
1152
1153     // Volume
1154     //
1155     if (gMC->CurrentVolName() != 0)
1156       cout << setw(4) << gMC->CurrentVolName() << "  ";
1157     else
1158       cout << setw(4) << "None"  << "  ";
1159
1160     // Process
1161     //
1162     TArrayI processes;
1163     Int_t nofProcesses = gMC->StepProcesses(processes);
1164     for(int ip=0;ip<nofProcesses; ip++)
1165       cout << TMCProcessName[processes[ip]]<<" / ";
1166
1167     cout << endl;
1168 }
1169
1170
1171 //_____________________________________________________________________________
1172 Double_t Ermilova(Double_t *x, Double_t *)
1173 {
1174   //
1175   // Calculates the delta-ray energy distribution according to Ermilova.
1176   // Logarithmic scale !
1177   //
1178
1179   Double_t energy;
1180   Double_t dpos;
1181   Double_t dnde;
1182
1183   Int_t    pos1, pos2;
1184
1185   const Int_t kNv = 31;
1186
1187   Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
1188                      , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
1189                      , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
1190                      , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
1191                      , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
1192                      , 9.4727, 9.9035,10.3735,10.5966,10.8198
1193                      ,11.5129 };
1194
1195   Float_t vye[kNv] = { 80.0  , 31.0  , 23.3  , 21.1  , 21.0
1196                      , 20.9  , 20.8  , 20.0  , 16.0  , 11.0
1197                      ,  8.0  ,  6.0  ,  5.2  ,  4.6  ,  4.0
1198                      ,  3.5  ,  3.0  ,  1.4  ,  0.67 ,  0.44
1199                      ,  0.3  ,  0.18 ,  0.12 ,  0.08 ,  0.056
1200                      ,  0.04 ,  0.023,  0.015,  0.011,  0.01
1201                      ,  0.004 };
1202
1203   energy = x[0];
1204
1205   // Find the position
1206   pos1 = pos2 = 0;
1207   dpos = 0;
1208   do {
1209     dpos = energy - vxe[pos2++];
1210   }
1211   while (dpos > 0);
1212   pos2--;
1213   if (pos2 > kNv) pos2 = kNv - 1;
1214   pos1 = pos2 - 1;
1215
1216   // Differentiate between the sampling points
1217   dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1218
1219   return dnde;
1220
1221 }
1222
1223 //_____________________________________________________________________________
1224 Double_t IntSpecGeant(Double_t *x, Double_t *)
1225 {
1226   //
1227   // Integrated spectrum from Geant3
1228   //
1229
1230   const Int_t n_pts = 83;
1231   Double_t arr_e[n_pts] = {
1232     2.421257,  2.483278,  2.534301,  2.592230,
1233     2.672067,  2.813299,  3.015059,  3.216819,
1234     3.418579,  3.620338,  3.868209,  3.920198,
1235     3.978284,  4.063923,  4.186264,  4.308605,
1236     4.430946,  4.553288,  4.724261,  4.837736,
1237     4.999842,  5.161949,  5.324056,  5.486163,
1238     5.679688,  5.752998,  5.857728,  5.962457,
1239     6.067185,  6.171914,  6.315653,  6.393674,
1240     6.471694,  6.539689,  6.597658,  6.655627,
1241     6.710957,  6.763648,  6.816338,  6.876198,
1242     6.943227,  7.010257,  7.106285,  7.252151,
1243     7.460531,  7.668911,  7.877290,  8.085670,
1244     8.302979,  8.353585,  8.413120,  8.483500,
1245     8.541030,  8.592857,  8.668865,  8.820485,
1246     9.037086,  9.253686,  9.470286,  9.686887,
1247     9.930838,  9.994655, 10.085822, 10.176990,
1248     10.268158, 10.359325, 10.503614, 10.627565,
1249     10.804637, 10.981709, 11.158781, 11.335854,
1250     11.593397, 11.781165, 12.049404, 12.317644,
1251     12.585884, 12.854123, 14.278421, 16.975889,
1252     20.829416, 24.682943, 28.536469
1253   };
1254   Double_t arr_dndx[n_pts] = {
1255     19.344431, 18.664679, 18.136106, 17.567745,
1256     16.836426, 15.677382, 14.281277, 13.140237,
1257     12.207677, 11.445510, 10.697049, 10.562296,
1258     10.414673, 10.182341,  9.775256,  9.172330,
1259     8.240271,  6.898587,  4.808303,  3.889751,
1260     3.345288,  3.093431,  2.897347,  2.692470,
1261     2.436222,  2.340029,  2.208579,  2.086489,
1262     1.975535,  1.876519,  1.759626,  1.705024,
1263     1.656374,  1.502638,  1.330566,  1.200697,
1264     1.101168,  1.019323,  0.943867,  0.851951,
1265     0.755229,  0.671576,  0.570675,  0.449672,
1266     0.326722,  0.244225,  0.188225,  0.149608,
1267     0.121529,  0.116289,  0.110636,  0.103490,
1268     0.096147,  0.089191,  0.079780,  0.063927,
1269     0.047642,  0.036341,  0.028250,  0.022285,
1270     0.017291,  0.016211,  0.014802,  0.013533,
1271     0.012388,  0.011352,  0.009803,  0.008537,
1272     0.007039,  0.005829,  0.004843,  0.004034,
1273     0.003101,  0.002564,  0.001956,  0.001494,
1274     0.001142,  0.000873,  0.000210,  0.000014,
1275     0.000000,  0.000000,  0.000000
1276   };
1277
1278   Int_t i;
1279   Double_t energy = x[0];
1280   Double_t dnde;
1281
1282   for( i = 0 ; i < n_pts ; i++ )
1283     if( energy < arr_e[i] ) break;
1284
1285   if( i == 0 )
1286     AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");
1287
1288   // Differentiate
1289   dnde = (arr_dndx[i-1] - arr_dndx[i]) / (arr_e[i] - arr_e[i-1]);
1290
1291   return dnde;
1292
1293 }