Make some calculations optional for HLT
[u/mrichter/AliRoot.git] / TRD / AliTRDv1.cxx
CommitLineData
4c039060 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
88cb7938 16/* $Id$ */
4c039060 17
030b4415 18////////////////////////////////////////////////////////////////////////////
19// //
20// Transition Radiation Detector version 1 -- slow simulator //
21// //
22////////////////////////////////////////////////////////////////////////////
fe4da5cc 23
769257f4 24#include <stdlib.h>
25
1819f4bb 26#include <TLorentzVector.h>
88cb7938 27#include <TMath.h>
28#include <TRandom.h>
29#include <TVector.h>
30#include <TVirtualMC.h>
f57bb418 31#include <TGeoManager.h>
268f57b1 32#include <TGeoMatrix.h>
170c35f1 33#include <TGeoPhysicalNode.h>
fe4da5cc 34
d3f347ff 35#include "AliConst.h"
45160b1f 36#include "AliLog.h"
e6add757 37#include "AliTrackReference.h"
45160b1f 38#include "AliMC.h"
88cb7938 39#include "AliRun.h"
02cb65d6 40#include "AliGeomManager.h"
030b4415 41
88cb7938 42#include "AliTRDgeometry.h"
a076fc2f 43#include "AliTRDCommonParam.h"
793ff80c 44#include "AliTRDhit.h"
cb2f9e9b 45#include "AliTRDsimTR.h"
88cb7938 46#include "AliTRDv1.h"
851d3db9 47
fe4da5cc 48ClassImp(AliTRDv1)
8230f242 49
50//_____________________________________________________________________________
030b4415 51AliTRDv1::AliTRDv1()
52 :AliTRD()
53 ,fTRon(kFALSE)
54 ,fTR(NULL)
030b4415 55 ,fStepSize(0)
f2979d08 56 ,fWion(0)
8230f242 57{
58 //
59 // Default constructor
60 //
61
8230f242 62}
63
fe4da5cc 64//_____________________________________________________________________________
65AliTRDv1::AliTRDv1(const char *name, const char *title)
030b4415 66 :AliTRD(name,title)
67 ,fTRon(kTRUE)
68 ,fTR(NULL)
030b4415 69 ,fStepSize(0.1)
f2979d08 70 ,fWion(0)
fe4da5cc 71{
72 //
851d3db9 73 // Standard constructor for Transition Radiation Detector version 1
fe4da5cc 74 //
82bbf98a 75
5c7f4665 76 SetBufferSize(128000);
77
a076fc2f 78 if (AliTRDCommonParam::Instance()->IsXenon()) {
f2979d08 79 fWion = 23.53; // Ionization energy XeCO2 (85/15)
80 }
a076fc2f 81 else if (AliTRDCommonParam::Instance()->IsArgon()) {
f2979d08 82 fWion = 27.21; // Ionization energy ArCO2 (82/18)
83 }
84 else {
85 AliFatal("Wrong gas mixture");
86 exit(1);
87 }
88
5c7f4665 89}
90
91//_____________________________________________________________________________
92AliTRDv1::~AliTRDv1()
93{
dd9a6ee3 94 //
95 // AliTRDv1 destructor
96 //
82bbf98a 97
030b4415 98 if (fTR) {
99 delete fTR;
100 fTR = 0;
101 }
82bbf98a 102
fe4da5cc 103}
104
dd9a6ee3 105//_____________________________________________________________________________
f57bb418 106void AliTRDv1::AddAlignableVolumes() const
107{
108 //
109 // Create entries for alignable volumes associating the symbolic volume
110 // name with the corresponding volume path. Needs to be syncronized with
111 // eventual changes in the geometry.
112 //
113
114 TString volPath;
115 TString symName;
116
92cd97ad 117 TString vpStr = "ALIC_1/B077_1/BSEGMO";
118 TString vpApp1 = "_1/BTRD";
119 TString vpApp2 = "_1";
120 TString vpApp3a = "/UTR1_1/UTS1_1/UTI1_1/UT";
121 TString vpApp3b = "/UTR2_1/UTS2_1/UTI2_1/UT";
122 TString vpApp3c = "/UTR3_1/UTS3_1/UTI3_1/UT";
f57bb418 123
92cd97ad 124 TString snStr = "TRD/sm";
125 TString snApp1 = "/st";
126 TString snApp2 = "/pl";
f57bb418 127
128 //
129 // The super modules
130 // The symbolic names are: TRD/sm00
131 // ...
132 // TRD/sm17
133 //
053767a4 134 for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
f57bb418 135
136 volPath = vpStr;
053767a4 137 volPath += isector;
f57bb418 138 volPath += vpApp1;
053767a4 139 volPath += isector;
f57bb418 140 volPath += vpApp2;
141
142 symName = snStr;
053767a4 143 symName += Form("%02d",isector);
f57bb418 144
145 gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data());
146
147 }
148
149 //
150 // The readout chambers
151 // The symbolic names are: TRD/sm00/st0/pl0
152 // ...
153 // TRD/sm17/st4/pl5
154 //
02cb65d6 155 AliGeomManager::ELayerID idTRD1 = AliGeomManager::kTRD1;
156 Int_t layer, modUID;
157
053767a4 158 for (Int_t isector = 0; isector < AliTRDgeometry::Nsector(); isector++) {
8bf0cd64 159
053767a4 160 if (fGeometry->GetSMstatus(isector) == 0) continue;
8bf0cd64 161
053767a4 162 for (Int_t istack = 0; istack < AliTRDgeometry::Nstack(); istack++) {
163 for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) {
f57bb418 164
053767a4 165 layer = idTRD1 + ilayer;
166 modUID = AliGeomManager::LayerToVolUIDSafe(layer,isector*5+istack);
167
168 Int_t idet = AliTRDgeometry::GetDetectorSec(ilayer,istack);
f57bb418 169
170 volPath = vpStr;
053767a4 171 volPath += isector;
f57bb418 172 volPath += vpApp1;
053767a4 173 volPath += isector;
f57bb418 174 volPath += vpApp2;
053767a4 175 switch (isector) {
92cd97ad 176 case 13:
177 case 14:
178 case 15:
053767a4 179 if (istack == 2) {
92cd97ad 180 continue;
181 }
182 volPath += vpApp3c;
183 break;
184 case 11:
185 case 12:
186 volPath += vpApp3b;
187 break;
188 default:
189 volPath += vpApp3a;
190 };
f57bb418 191 volPath += Form("%02d",idet);
192 volPath += vpApp2;
193
194 symName = snStr;
053767a4 195 symName += Form("%02d",isector);
f57bb418 196 symName += snApp1;
053767a4 197 symName += istack;
f57bb418 198 symName += snApp2;
053767a4 199 symName += ilayer;
f57bb418 200
51a5f1d0 201 TGeoPNEntry *alignableEntry =
02cb65d6 202 gGeoManager->SetAlignableEntry(symName.Data(),volPath.Data(),modUID);
f57bb418 203
170c35f1 204 // Add the tracking to local matrix following the TPC example
51a5f1d0 205 if (alignableEntry) {
02cb65d6 206 TGeoHMatrix *globMatrix = alignableEntry->GetGlobalOrig();
053767a4 207 Double_t sectorAngle = 20.0 * (isector % 18) + 10.0;
51a5f1d0 208 TGeoHMatrix *t2lMatrix = new TGeoHMatrix();
209 t2lMatrix->RotateZ(sectorAngle);
210 t2lMatrix->MultiplyLeft(&(globMatrix->Inverse()));
211 alignableEntry->SetMatrix(t2lMatrix);
212 }
213 else {
214 AliError(Form("Alignable entry %s is not valid!",symName.Data()));
215 }
8bf0cd64 216
f57bb418 217 }
218 }
219 }
220
221}
222
223//_____________________________________________________________________________
fe4da5cc 224void AliTRDv1::CreateGeometry()
225{
226 //
851d3db9 227 // Create the GEANT geometry for the Transition Radiation Detector - Version 1
5c7f4665 228 // This version covers the full azimuth.
d3f347ff 229 //
230
82bbf98a 231 // Check that FRAME is there otherwise we have no place where to put the TRD
8230f242 232 AliModule* frame = gAlice->GetModule("FRAME");
030b4415 233 if (!frame) {
234 AliError("TRD needs FRAME to be present\n");
235 return;
236 }
d3f347ff 237
82bbf98a 238 // Define the chambers
239 AliTRD::CreateGeometry();
d3f347ff 240
fe4da5cc 241}
242
243//_____________________________________________________________________________
244void AliTRDv1::CreateMaterials()
245{
246 //
851d3db9 247 // Create materials for the Transition Radiation Detector version 1
fe4da5cc 248 //
82bbf98a 249
d3f347ff 250 AliTRD::CreateMaterials();
82bbf98a 251
fe4da5cc 252}
253
254//_____________________________________________________________________________
793ff80c 255void AliTRDv1::CreateTRhit(Int_t det)
256{
257 //
258 // Creates an electron cluster from a TR photon.
259 // The photon is assumed to be created a the end of the radiator. The
260 // distance after which it deposits its energy takes into account the
261 // absorbtion of the entrance window and of the gas mixture in drift
262 // volume.
263 //
264
793ff80c 265 // Maximum number of TR photons per track
266 const Int_t kNTR = 50;
267
030b4415 268 TLorentzVector mom;
269 TLorentzVector pos;
793ff80c 270
ce0d6231 271 Float_t eTR[kNTR];
272 Int_t nTR;
793ff80c 273
ce0d6231 274 // Create TR photons
275 gMC->TrackMomentum(mom);
276 Float_t pTot = mom.Rho();
277 fTR->CreatePhotons(11,pTot,nTR,eTR);
278 if (nTR > kNTR) {
279 AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
280 }
f73816f5 281
ce0d6231 282 // Loop through the TR photons
283 for (Int_t iTR = 0; iTR < nTR; iTR++) {
793ff80c 284
ce0d6231 285 Float_t energyMeV = eTR[iTR] * 0.001;
286 Float_t energyeV = eTR[iTR] * 1000.0;
287 Float_t absLength = 0.0;
288 Float_t sigma = 0.0;
793ff80c 289
ce0d6231 290 // Take the absorbtion in the entrance window into account
291 Double_t muMy = fTR->GetMuMy(energyMeV);
292 sigma = muMy * fFoilDensity;
293 if (sigma > 0.0) {
294 absLength = gRandom->Exp(1.0/sigma);
295 if (absLength < AliTRDgeometry::MyThick()) {
842287f2 296 continue;
297 }
ce0d6231 298 }
299 else {
300 continue;
301 }
793ff80c 302
ce0d6231 303 // The absorbtion cross sections in the drift gas
304 // Gas-mixture (Xe/CO2)
f2979d08 305 Double_t muNo = 0.0;
a076fc2f 306 if (AliTRDCommonParam::Instance()->IsXenon()) {
f2979d08 307 muNo = fTR->GetMuXe(energyMeV);
308 }
a076fc2f 309 else if (AliTRDCommonParam::Instance()->IsArgon()) {
f2979d08 310 muNo = fTR->GetMuAr(energyMeV);
311 }
ce0d6231 312 Double_t muCO = fTR->GetMuCO(energyMeV);
f2979d08 313 sigma = (fGasNobleFraction * muNo + (1.0 - fGasNobleFraction) * muCO)
314 * fGasDensity
315 * fTR->GetTemp();
ce0d6231 316
317 // The distance after which the energy of the TR photon
318 // is deposited.
319 if (sigma > 0.0) {
320 absLength = gRandom->Exp(1.0/sigma);
321 if (absLength > (AliTRDgeometry::DrThick()
322 + AliTRDgeometry::AmThick())) {
842287f2 323 continue;
324 }
ce0d6231 325 }
326 else {
327 continue;
328 }
793ff80c 329
ce0d6231 330 // The position of the absorbtion
331 Float_t posHit[3];
332 gMC->TrackPosition(pos);
333 posHit[0] = pos[0] + mom[0] / pTot * absLength;
334 posHit[1] = pos[1] + mom[1] / pTot * absLength;
335 posHit[2] = pos[2] + mom[2] / pTot * absLength;
793ff80c 336
ce0d6231 337 // Create the charge
f2979d08 338 Int_t q = ((Int_t) (energyeV / fWion));
793ff80c 339
ce0d6231 340 // Add the hit to the array. TR photon hits are marked
341 // by negative charge
342 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
343 ,det
344 ,posHit
345 ,-q
25ca55ce 346 ,gMC->TrackTime()*1.0e06
d4c6453d 347 ,kTRUE);
793ff80c 348
349 }
350
351}
352
353//_____________________________________________________________________________
5c7f4665 354void AliTRDv1::Init()
355{
356 //
357 // Initialise Transition Radiation Detector after geometry has been built.
5c7f4665 358 //
359
360 AliTRD::Init();
361
45160b1f 362 AliDebug(1,"Slow simulator\n");
bd0f8685 363
364 // Switch on TR simulation as default
365 if (!fTRon) {
45160b1f 366 AliInfo("TR simulation off");
bd0f8685 367 }
368 else {
cb2f9e9b 369 fTR = new AliTRDsimTR();
bd0f8685 370 }
5c7f4665 371
45160b1f 372 AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
5c7f4665 373
fe4da5cc 374}
375
376//_____________________________________________________________________________
5c7f4665 377void AliTRDv1::StepManager()
378{
379 //
a328fff9 380 // Slow simulator. Every charged track produces electron cluster as hits
381 // along its path across the drift volume. The step size is fixed in
382 // this version of the step manager.
383 //
f2979d08 384 // Works for Xe/CO2 as well as Ar/CO2
385 //
a328fff9 386
ce0d6231 387 // PDG code electron
388 const Int_t kPdgElectron = 11;
389
053767a4 390 Int_t layer = 0;
391 Int_t stack = 0;
392 Int_t sector = 0;
393 Int_t det = 0;
a328fff9 394 Int_t qTot;
395
396 Float_t hits[3];
397 Double_t eDep;
398
399 Bool_t drRegion = kFALSE;
400 Bool_t amRegion = kFALSE;
401
2c8bf4aa 402 TString cIdPath;
403 Char_t cIdSector[3];
404 cIdSector[2] = 0;
405
a328fff9 406 TString cIdCurrent;
407 TString cIdSensDr = "J";
408 TString cIdSensAm = "K";
409 Char_t cIdChamber[3];
2c8bf4aa 410 cIdChamber[2] = 0;
a328fff9 411
030b4415 412 TLorentzVector pos;
413 TLorentzVector mom;
a328fff9 414
053767a4 415 const Int_t kNlayer = AliTRDgeometry::Nlayer();
416 const Int_t kNstack = AliTRDgeometry::Nstack();
417 const Int_t kNdetsec = kNlayer * kNstack;
a328fff9 418
030b4415 419 const Double_t kBig = 1.0e+12;
a328fff9 420 const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment
421
422 // Set the maximum step size to a very large number for all
423 // neutral particles and those outside the driftvolume
424 gMC->SetMaxStep(kBig);
425
426 // If not charged track or already stopped or disappeared, just return.
427 if ((!gMC->TrackCharge()) ||
ce0d6231 428 gMC->IsTrackDisappeared()) {
429 return;
430 }
a328fff9 431
432 // Inside a sensitive volume?
433 cIdCurrent = gMC->CurrentVolName();
434
ce0d6231 435 if (cIdSensDr == cIdCurrent[1]) {
436 drRegion = kTRUE;
437 }
438 if (cIdSensAm == cIdCurrent[1]) {
439 amRegion = kTRUE;
440 }
a328fff9 441
030b4415 442 if ((!drRegion) &&
443 (!amRegion)) {
444 return;
445 }
a328fff9 446
447 // The hit coordinates and charge
448 gMC->TrackPosition(pos);
449 hits[0] = pos[0];
450 hits[1] = pos[1];
451 hits[2] = pos[2];
452
2c8bf4aa 453 // The sector number (0 - 17), according to standard coordinate system
454 cIdPath = gGeoManager->GetPath();
455 cIdSector[0] = cIdPath[21];
456 cIdSector[1] = cIdPath[22];
053767a4 457 sector = atoi(cIdSector);
a328fff9 458
459 // The plane and chamber number
030b4415 460 cIdChamber[0] = cIdCurrent[2];
461 cIdChamber[1] = cIdCurrent[3];
a328fff9 462 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
053767a4 463 stack = ((Int_t) idChamber / kNlayer);
464 layer = ((Int_t) idChamber % kNlayer);
e0d47c25 465
030b4415 466 // The detector number
053767a4 467 det = fGeometry->GetDetector(layer,stack,sector);
030b4415 468
25ca55ce 469 // 0: InFlight 1:Entering 2:Exiting
030b4415 470 Int_t trkStat = 0;
a328fff9 471
472 // Special hits only in the drift region
ce0d6231 473 if ((drRegion) &&
474 (gMC->IsTrackEntering())) {
a328fff9 475
ce0d6231 476 // Create a track reference at the entrance of each
477 // chamber that contains the momentum components of the particle
478 gMC->TrackMomentum(mom);
e6add757 479 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
ce0d6231 480 trkStat = 1;
a328fff9 481
ce0d6231 482 // Create the hits from TR photons if electron/positron is
483 // entering the drift volume
484 if ((fTR) &&
485 (TMath::Abs(gMC->TrackPid()) == kPdgElectron)) {
f2979d08 486 CreateTRhit(det);
030b4415 487 }
a328fff9 488
489 }
ce0d6231 490 else if ((amRegion) &&
491 (gMC->IsTrackExiting())) {
492
493 // Create a track reference at the exit of each
494 // chamber that contains the momentum components of the particle
495 gMC->TrackMomentum(mom);
e6add757 496 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kTRD);
ce0d6231 497 trkStat = 2;
498
499 }
a328fff9 500
501 // Calculate the charge according to GEANT Edep
502 // Create a new dEdx hit
503 eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
f2979d08 504 qTot = (Int_t) (eDep / fWion);
ce0d6231 505 if ((qTot) ||
506 (trkStat)) {
507 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
508 ,det
509 ,hits
510 ,qTot
25ca55ce 511 ,gMC->TrackTime()*1.0e06
ce0d6231 512 ,drRegion);
513 }
a328fff9 514
515 // Set Maximum Step Size
516 // Produce only one hit if Ekin is below cutoff
030b4415 517 if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) {
518 return;
519 }
a328fff9 520 gMC->SetMaxStep(fStepSize);
521
522}