TRD/AliTRDv1.cxx

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