]>
Commit | Line | Data |
---|---|---|
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 | ||
16 | /* | |
17 | $Log$ | |
ab76897d | 18 | Revision 1.14 1999/10/04 14:48:07 fca |
19 | Avoid warnings on non-ansi compiler HP-UX CC | |
20 | ||
7d83513f | 21 | Revision 1.13 1999/09/29 09:24:35 fca |
22 | Introduction of the Copyright and cvs Log | |
23 | ||
4c039060 | 24 | */ |
25 | ||
fe4da5cc | 26 | /////////////////////////////////////////////////////////////////////////////// |
27 | // // | |
28 | // Transition Radiation Detector version 2 -- detailed simulation // | |
29 | // // | |
30 | //Begin_Html | |
31 | /* | |
1439f98e | 32 | <img src="picts/AliTRDv2Class.gif"> |
fe4da5cc | 33 | */ |
34 | //End_Html | |
35 | // // | |
36 | // // | |
37 | /////////////////////////////////////////////////////////////////////////////// | |
38 | ||
ab76897d | 39 | #include <stdlib.h> |
40 | ||
fe4da5cc | 41 | #include <TMath.h> |
fe4da5cc | 42 | #include <TVector.h> |
99d5402e | 43 | #include <TRandom.h> |
fe4da5cc | 44 | |
fe4da5cc | 45 | #include "AliTRDv2.h" |
99d5402e | 46 | #include "AliTRDmatrix.h" |
fe4da5cc | 47 | #include "AliRun.h" |
fe4da5cc | 48 | #include "AliMC.h" |
d3f347ff | 49 | #include "AliConst.h" |
fe4da5cc | 50 | |
51 | ClassImp(AliTRDv2) | |
52 | ||
53 | //_____________________________________________________________________________ | |
54 | AliTRDv2::AliTRDv2(const char *name, const char *title) | |
55 | :AliTRD(name, title) | |
56 | { | |
57 | // | |
58 | // Standard constructor for Transition Radiation Detector version 2 | |
59 | // | |
82bbf98a | 60 | |
99d5402e | 61 | fIdSens = 0; |
82bbf98a | 62 | |
99d5402e | 63 | fIdSpace1 = 0; |
64 | fIdSpace2 = 0; | |
65 | fIdSpace3 = 0; | |
82bbf98a | 66 | |
99d5402e | 67 | fIdChamber1 = 0; |
68 | fIdChamber2 = 0; | |
69 | fIdChamber3 = 0; | |
82bbf98a | 70 | |
99d5402e | 71 | fSensSelect = 0; |
72 | fSensPlane = 0; | |
73 | fSensChamber = 0; | |
74 | fSensSector = 0; | |
82bbf98a | 75 | |
7d83513f | 76 | Int_t iplan; |
77 | ||
78 | for (iplan = 0; iplan < kNplan; iplan++) { | |
99d5402e | 79 | for (Int_t icham = 0; icham < kNcham; icham++) { |
80 | fRowMax[iplan][icham] = 0; | |
81 | } | |
82 | fColMax[iplan] = 0; | |
83 | } | |
84 | fTimeMax = 0; | |
85 | ||
86 | fRowPadSize = 0; | |
87 | fColPadSize = 0; | |
88 | fTimeBinSize = 0; | |
89 | ||
90 | fGasGain = 0; | |
91 | fNoise = 0; | |
92 | fChipGain = 0; | |
93 | fADCoutRange = 0; | |
94 | fADCinRange = 0; | |
95 | fADCthreshold = 0; | |
96 | ||
97 | fDiffusionT = 0; | |
98 | fDiffusionL = 0; | |
99 | ||
100 | fDeltaE = NULL; | |
d3f347ff | 101 | |
fe4da5cc | 102 | SetBufferSize(128000); |
82bbf98a | 103 | |
fe4da5cc | 104 | } |
d3f347ff | 105 | |
99d5402e | 106 | //_____________________________________________________________________________ |
d3f347ff | 107 | AliTRDv2::~AliTRDv2() |
108 | { | |
82bbf98a | 109 | |
99d5402e | 110 | if (fDeltaE) delete fDeltaE; |
82bbf98a | 111 | |
d3f347ff | 112 | } |
fe4da5cc | 113 | |
114 | //_____________________________________________________________________________ | |
115 | void AliTRDv2::CreateGeometry() | |
116 | { | |
117 | // | |
82bbf98a | 118 | // Create the GEANT geometry for the Transition Radiation Detector - Version 2 |
fe4da5cc | 119 | // This version covers the full azimuth. |
d3f347ff | 120 | // |
82bbf98a | 121 | // Author: Christoph Blume (C.Blume@gsi.de) 20/07/99 |
fe4da5cc | 122 | // |
d3f347ff | 123 | |
82bbf98a | 124 | Float_t xpos, ypos, zpos; |
d3f347ff | 125 | |
82bbf98a | 126 | // Check that FRAME is there otherwise we have no place where to put the TRD |
127 | AliModule* FRAME = gAlice->GetModule("FRAME"); | |
128 | if (!FRAME) return; | |
d3f347ff | 129 | |
82bbf98a | 130 | // Define the chambers |
131 | AliTRD::CreateGeometry(); | |
fe4da5cc | 132 | |
82bbf98a | 133 | // Position the the TRD-sectors in all TRD-volumes in the spaceframe |
134 | xpos = 0.; | |
135 | ypos = 0.; | |
136 | zpos = 0.; | |
137 | gMC->Gspos("TRD ",1,"BTR1",xpos,ypos,zpos,0,"ONLY"); | |
138 | gMC->Gspos("TRD ",2,"BTR2",xpos,ypos,zpos,0,"ONLY"); | |
139 | gMC->Gspos("TRD ",3,"BTR3",xpos,ypos,zpos,0,"ONLY"); | |
fe4da5cc | 140 | |
fe4da5cc | 141 | } |
142 | ||
143 | //_____________________________________________________________________________ | |
144 | void AliTRDv2::CreateMaterials() | |
145 | { | |
146 | // | |
147 | // Create materials for the Transition Radiation Detector version 2 | |
148 | // | |
82bbf98a | 149 | |
fe4da5cc | 150 | AliTRD::CreateMaterials(); |
82bbf98a | 151 | |
fe4da5cc | 152 | } |
153 | ||
154 | //_____________________________________________________________________________ | |
99d5402e | 155 | void AliTRDv2::Diffusion(Float_t driftlength, Float_t *xyz) |
fe4da5cc | 156 | { |
157 | // | |
99d5402e | 158 | // Applies the diffusion smearing to the position of a single electron |
fe4da5cc | 159 | // |
d3f347ff | 160 | |
99d5402e | 161 | if ((driftlength > 0) && |
162 | (driftlength < kDrThick)) { | |
163 | Float_t driftSqrt = TMath::Sqrt(driftlength); | |
164 | Float_t sigmaT = driftSqrt * fDiffusionT; | |
165 | Float_t sigmaL = driftSqrt * fDiffusionL; | |
166 | xyz[0] = gRandom->Gaus(xyz[0], sigmaL); | |
167 | xyz[1] = gRandom->Gaus(xyz[1], sigmaT); | |
168 | xyz[2] = gRandom->Gaus(xyz[2], sigmaT); | |
169 | } | |
170 | else { | |
171 | xyz[0] = 0.0; | |
172 | xyz[1] = 0.0; | |
173 | xyz[2] = 0.0; | |
174 | } | |
175 | ||
176 | } | |
177 | ||
178 | //_____________________________________________________________________________ | |
179 | void AliTRDv2::Hits2Digits() | |
180 | { | |
181 | // | |
182 | // Creates TRD digits from hits. This procedure includes the following: | |
183 | // - Diffusion | |
184 | // - Gas gain including fluctuations | |
185 | // - Pad-response (simple Gaussian approximation) | |
186 | // - Electronics noise | |
187 | // - Electronics gain | |
188 | // - Digitization | |
189 | // - ADC threshold | |
190 | // The corresponding parameter can be adjusted via the various Set-functions. | |
191 | // If these parameters are not explicitly set, default values are used (see | |
192 | // Init-function). | |
193 | // To produce digits from a galice.root file with TRD-hits use the | |
194 | // digitsCreate.C macro. | |
195 | // | |
196 | ||
197 | printf(" Start creating digits\n"); | |
198 | ||
199 | /////////////////////////////////////////////////////////////// | |
200 | // Parameter | |
201 | /////////////////////////////////////////////////////////////// | |
202 | ||
203 | // Converts number of electrons to fC | |
204 | const Float_t el2fC = 1.602E-19 * 1.0E15; | |
205 | ||
206 | /////////////////////////////////////////////////////////////// | |
207 | ||
208 | Int_t nBytes = 0; | |
209 | ||
210 | AliTRDhit *TRDhit; | |
211 | ||
7d83513f | 212 | Int_t iplan; |
213 | Int_t iRow; | |
214 | ||
99d5402e | 215 | // Position of pad 0,0,0 |
216 | // | |
217 | // chambers seen from the top: | |
218 | // +----------------------------+ | |
219 | // | | | |
220 | // | | ^ | |
221 | // | | rphi| | |
222 | // | | | | |
223 | // |0 | | | |
224 | // +----------------------------+ +------> | |
225 | // z | |
226 | // chambers seen from the side: ^ | |
227 | // +----------------------------+ time| | |
228 | // | | | | |
229 | // |0 | | | |
230 | // +----------------------------+ +------> | |
231 | // z | |
232 | // | |
233 | // The pad row (z-direction) | |
234 | Float_t row0[kNplan][kNcham]; | |
7d83513f | 235 | for (iplan = 0; iplan < kNplan; iplan++) { |
99d5402e | 236 | row0[iplan][0] = -fClengthI[iplan]/2. - fClengthM[iplan] - fClengthO[iplan] |
237 | + kCcthick; | |
238 | row0[iplan][1] = -fClengthI[iplan]/2. - fClengthM[iplan] | |
239 | + kCcthick; | |
240 | row0[iplan][2] = -fClengthI[iplan]/2. | |
241 | + kCcthick; | |
242 | row0[iplan][3] = fClengthI[iplan]/2. | |
243 | + kCcthick; | |
244 | row0[iplan][4] = fClengthI[iplan]/2. + fClengthM[iplan] | |
245 | + kCcthick; | |
246 | } | |
247 | // The pad column (rphi-direction) | |
248 | Float_t col0[kNplan]; | |
7d83513f | 249 | for (iplan = 0; iplan < kNplan; iplan++) { |
99d5402e | 250 | col0[iplan] = -fCwidth[iplan]/2. + kCcthick; |
251 | } | |
252 | // The time bucket | |
253 | Float_t time0[kNplan]; | |
7d83513f | 254 | for (iplan = 0; iplan < kNplan; iplan++) { |
99d5402e | 255 | time0[iplan] = kRmin + kCcframe/2. + kDrZpos - 0.5 * kDrThick |
256 | + iplan * (kCheight + kCspace); | |
257 | } | |
258 | ||
259 | // Get the pointer to the hit tree | |
260 | TTree *HitTree = gAlice->TreeH(); | |
261 | // Get the pointer to the digits tree | |
262 | TTree *DigitsTree = gAlice->TreeD(); | |
263 | ||
264 | // Get the number of entries in the hit tree | |
265 | // (Number of primary particles creating a hit somewhere) | |
266 | Int_t nTrack = (Int_t) HitTree->GetEntries(); | |
267 | ||
268 | Int_t chamBeg = 0; | |
269 | Int_t chamEnd = kNcham; | |
270 | if (fSensChamber) chamEnd = chamBeg = fSensChamber; | |
271 | Int_t planBeg = 0; | |
272 | Int_t planEnd = kNplan; | |
273 | if (fSensPlane) planEnd = planBeg = fSensPlane; | |
274 | Int_t sectBeg = 0; | |
275 | Int_t sectEnd = kNsect; | |
276 | if (fSensSector) sectEnd = sectBeg = fSensSector; | |
277 | ||
278 | // Loop through all the chambers | |
279 | for (Int_t icham = chamBeg; icham < chamEnd; icham++) { | |
7d83513f | 280 | for (iplan = planBeg; iplan < planEnd; iplan++) { |
99d5402e | 281 | for (Int_t isect = sectBeg; isect < sectEnd; isect++) { |
282 | ||
283 | printf(" Digitizing chamber %d, plane %d, sector %d\n" | |
284 | ,icham+1,iplan+1,isect+1); | |
285 | ||
286 | // Create a detector matrix to keep the signal and track numbers | |
287 | AliTRDmatrix *matrix = new AliTRDmatrix(fRowMax[iplan][icham] | |
288 | ,fColMax[iplan] | |
289 | ,fTimeMax | |
290 | ,isect+1,icham+1,iplan+1); | |
291 | ||
292 | // Loop through all entries in the tree | |
293 | for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) { | |
294 | ||
295 | gAlice->ResetHits(); | |
296 | nBytes += HitTree->GetEvent(iTrack); | |
297 | ||
298 | // Get the number of hits in the TRD created by this particle | |
299 | Int_t nHit = fHits->GetEntriesFast(); | |
300 | ||
301 | // Loop through the TRD hits | |
302 | for (Int_t iHit = 0; iHit < nHit; iHit++) { | |
303 | ||
304 | if (!(TRDhit = (AliTRDhit *) fHits->UncheckedAt(iHit))) | |
305 | continue; | |
306 | ||
307 | Float_t x = TRDhit->fX; | |
308 | Float_t y = TRDhit->fY; | |
309 | Float_t z = TRDhit->fZ; | |
310 | Float_t q = TRDhit->fQ; | |
311 | Int_t track = TRDhit->fTrack; | |
312 | Int_t plane = TRDhit->fPlane; | |
313 | Int_t sector = TRDhit->fSector; | |
314 | Int_t chamber = TRDhit->fChamber; | |
315 | ||
316 | if ((sector != isect+1) || | |
317 | (plane != iplan+1) || | |
318 | (chamber != icham+1)) | |
319 | continue; | |
320 | ||
321 | // Rotate the sectors on top of each other | |
322 | Float_t phi = 2.0 * kPI / (Float_t) kNsect | |
323 | * ((Float_t) sector - 0.5); | |
324 | Float_t xRot = -x * TMath::Cos(phi) + y * TMath::Sin(phi); | |
325 | Float_t yRot = x * TMath::Sin(phi) + y * TMath::Cos(phi); | |
326 | Float_t zRot = z; | |
327 | ||
328 | // The hit position in pad coordinates (center pad) | |
329 | // The pad row (z-direction) | |
330 | Int_t rowH = (Int_t) ((zRot - row0[iplan][icham]) / fRowPadSize); | |
331 | // The pad column (rphi-direction) | |
332 | Int_t colH = (Int_t) ((yRot - col0[iplan] ) / fColPadSize); | |
333 | // The time bucket | |
334 | Int_t timeH = (Int_t) ((xRot - time0[iplan] ) / fTimeBinSize); | |
335 | ||
336 | // Array to sum up the signal in a box surrounding the | |
337 | // hit postition | |
338 | const Int_t timeBox = 5; | |
339 | const Int_t colBox = 7; | |
340 | const Int_t rowBox = 5; | |
341 | Float_t signalSum[rowBox][colBox][timeBox]; | |
7d83513f | 342 | for (iRow = 0; iRow < rowBox; iRow++ ) { |
99d5402e | 343 | for (Int_t iCol = 0; iCol < colBox; iCol++ ) { |
344 | for (Int_t iTime = 0; iTime < timeBox; iTime++) { | |
345 | signalSum[iRow][iCol][iTime] = 0; | |
346 | } | |
347 | } | |
348 | } | |
349 | ||
350 | // Loop over all electrons of this hit | |
351 | Int_t nEl = (Int_t) q; | |
352 | for (Int_t iEl = 0; iEl < nEl; iEl++) { | |
353 | ||
354 | // Apply the diffusion smearing | |
355 | Float_t driftlength = xRot - time0[iplan]; | |
356 | Float_t xyz[3]; | |
357 | xyz[0] = xRot; | |
358 | xyz[1] = yRot; | |
359 | xyz[2] = zRot; | |
360 | Diffusion(driftlength,xyz); | |
361 | ||
362 | // At this point absorption effects that depend on the | |
363 | // driftlength could be taken into account. | |
364 | ||
365 | // The electron position and the distance to the hit position | |
366 | // in pad units | |
367 | // The pad row (z-direction) | |
368 | Int_t rowE = (Int_t) ((xyz[2] - row0[iplan][icham]) / fRowPadSize); | |
369 | Int_t rowD = rowH - rowE; | |
370 | // The pad column (rphi-direction) | |
371 | Int_t colE = (Int_t) ((xyz[1] - col0[iplan] ) / fColPadSize); | |
372 | Int_t colD = colH - colE; | |
373 | // The time bucket | |
374 | Int_t timeE = (Int_t) ((xyz[0] - time0[iplan] ) / fTimeBinSize); | |
375 | Int_t timeD = timeH - timeE; | |
376 | ||
377 | // Apply the gas gain including fluctuations | |
378 | Int_t signal = (Int_t) (-fGasGain * TMath::Log(gRandom->Rndm())); | |
379 | ||
380 | // The distance of the electron to the center of the pad | |
381 | // in units of pad width | |
382 | Float_t dist = (xyz[1] - col0[iplan] - (colE + 0.5) * fColPadSize) | |
383 | / fColPadSize; | |
384 | ||
385 | // Sum up the signal in the different pixels | |
386 | // and apply the pad response | |
387 | Int_t rowIdx = rowD + (Int_t) ( rowBox / 2); | |
388 | Int_t colIdx = colD + (Int_t) ( colBox / 2); | |
389 | Int_t timeIdx = timeD + (Int_t) (timeBox / 2); | |
390 | signalSum[rowIdx][colIdx-1][timeIdx] += PadResponse(dist-1.) * signal; | |
391 | signalSum[rowIdx][colIdx ][timeIdx] += PadResponse(dist ) * signal; | |
392 | signalSum[rowIdx][colIdx+1][timeIdx] += PadResponse(dist+1.) * signal; | |
393 | ||
394 | } | |
395 | ||
396 | // Add the padcluster to the detector matrix | |
7d83513f | 397 | for (iRow = 0; iRow < rowBox; iRow++ ) { |
99d5402e | 398 | for (Int_t iCol = 0; iCol < colBox; iCol++ ) { |
399 | for (Int_t iTime = 0; iTime < timeBox; iTime++) { | |
400 | ||
401 | Int_t rowB = rowH + iRow - (Int_t) ( rowBox / 2); | |
402 | Int_t colB = colH + iCol - (Int_t) ( colBox / 2); | |
403 | Int_t timeB = timeH + iTime - (Int_t) (timeBox / 2); | |
404 | ||
405 | Float_t signalB = signalSum[iRow][iCol][iTime]; | |
406 | if (signalB > 0.0) { | |
407 | matrix->AddSignal(rowB,colB,timeB,signalB); | |
408 | if (!(matrix->AddTrack(rowB,colB,timeB,track))) | |
409 | printf("More than three tracks in a pixel!\n"); | |
410 | } | |
411 | ||
412 | } | |
413 | } | |
414 | } | |
415 | ||
416 | } | |
417 | ||
418 | } | |
419 | ||
420 | // Create the hits for this chamber | |
421 | for (Int_t iRow = 0; iRow < fRowMax[iplan][icham]; iRow++ ) { | |
422 | for (Int_t iCol = 0; iCol < fColMax[iplan] ; iCol++ ) { | |
423 | for (Int_t iTime = 0; iTime < fTimeMax ; iTime++) { | |
424 | ||
425 | Float_t signalAmp = matrix->GetSignal(iRow,iCol,iTime); | |
426 | ||
427 | // Add the noise | |
7d83513f | 428 | signalAmp = TMath::Max(gRandom->Gaus(signalAmp,fNoise),(Float_t) 0.0); |
99d5402e | 429 | // Convert to fC |
430 | signalAmp *= el2fC; | |
431 | // Convert to mV | |
432 | signalAmp *= fChipGain; | |
433 | // Convert to ADC counts | |
434 | Int_t adc = (Int_t) (signalAmp * (fADCoutRange / fADCinRange)); | |
435 | ||
436 | // Apply threshold on ADC value | |
437 | if (adc > fADCthreshold) { | |
438 | ||
439 | Int_t trackSave[3]; | |
440 | for (Int_t ii = 0; ii < 3; ii++) { | |
441 | trackSave[ii] = matrix->GetTrack(iRow,iCol,iTime,ii); | |
442 | } | |
443 | ||
444 | Int_t digits[7]; | |
445 | digits[0] = matrix->GetSector(); | |
446 | digits[1] = matrix->GetChamber(); | |
447 | digits[2] = matrix->GetPlane(); | |
448 | digits[3] = iRow; | |
449 | digits[4] = iCol; | |
450 | digits[5] = iTime; | |
451 | digits[6] = adc; | |
452 | ||
453 | // Add this digit to the TClonesArray | |
454 | AddDigit(trackSave,digits); | |
455 | ||
456 | } | |
457 | ||
458 | } | |
459 | } | |
460 | } | |
461 | ||
462 | // Clean up | |
463 | delete matrix; | |
464 | ||
465 | } | |
466 | } | |
467 | } | |
468 | ||
469 | // Fill the digits-tree | |
470 | DigitsTree->Fill(); | |
471 | ||
472 | } | |
473 | ||
474 | //_____________________________________________________________________________ | |
475 | void AliTRDv2::Init() | |
476 | { | |
477 | // | |
478 | // Initialise Transition Radiation Detector after geometry has been built. | |
479 | // Includes the default settings of all parameter for the digitization. | |
480 | // | |
d3f347ff | 481 | |
ab76897d | 482 | printf("**************************************" |
483 | " TRD " | |
484 | "**************************************\n"); | |
485 | printf("\n Version 2 of TRD initialing, " | |
486 | "symmetric TRD\n\n"); | |
487 | ||
fe4da5cc | 488 | AliTRD::Init(); |
d3f347ff | 489 | |
ab76897d | 490 | |
491 | // | |
492 | // Check that FRAME is there otherwise we have no place where to | |
493 | // put TRD | |
494 | AliModule* FRAME=gAlice->GetModule("FRAME"); | |
495 | if(!FRAME) { | |
496 | Error("Ctor","TRD needs FRAME to be present\n"); | |
497 | exit(1); | |
498 | } else | |
499 | if(FRAME->IsVersion()!=1) { | |
500 | Error("Ctor","FRAME version 1 needed with this version of TRD\n"); | |
501 | exit(1); | |
502 | } | |
503 | ||
82bbf98a | 504 | if (fSensPlane) |
505 | printf(" Only plane %d is sensitive\n",fSensPlane); | |
506 | if (fSensChamber) | |
507 | printf(" Only chamber %d is sensitive\n",fSensChamber); | |
508 | if (fSensSector) | |
509 | printf(" Only sector %d is sensitive\n",fSensSector); | |
d3f347ff | 510 | |
99d5402e | 511 | for (Int_t i = 0; i < 80; i++) printf("*"); |
82bbf98a | 512 | printf("\n"); |
d3f347ff | 513 | |
99d5402e | 514 | // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2) |
515 | const Float_t kPoti = 12.1; | |
516 | // Maximum energy (50 keV); | |
517 | const Float_t kEend = 50000.0; | |
d3f347ff | 518 | // Ermilova distribution for the delta-ray spectrum |
82bbf98a | 519 | Float_t Poti = TMath::Log(kPoti); |
520 | Float_t Eend = TMath::Log(kEend); | |
d3f347ff | 521 | fDeltaE = new TF1("deltae",Ermilova,Poti,Eend,0); |
522 | ||
82bbf98a | 523 | // Identifier of the sensitive volume (drift region) |
524 | fIdSens = gMC->VolId("UL05"); | |
525 | ||
526 | // Identifier of the TRD-spaceframe volumina | |
527 | fIdSpace1 = gMC->VolId("B028"); | |
528 | fIdSpace2 = gMC->VolId("B029"); | |
529 | fIdSpace3 = gMC->VolId("B030"); | |
530 | ||
531 | // Identifier of the TRD-driftchambers | |
532 | fIdChamber1 = gMC->VolId("UCIO"); | |
533 | fIdChamber2 = gMC->VolId("UCIM"); | |
534 | fIdChamber3 = gMC->VolId("UCII"); | |
535 | ||
99d5402e | 536 | // The default pad dimensions |
537 | if (!(fRowPadSize)) fRowPadSize = 4.5; | |
538 | if (!(fColPadSize)) fColPadSize = 1.0; | |
539 | if (!(fTimeBinSize)) fTimeBinSize = 0.1; | |
540 | ||
541 | // The maximum number of pads | |
542 | for (Int_t iplan = 0; iplan < kNplan; iplan++) { | |
543 | // Rows | |
544 | fRowMax[iplan][0] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick) | |
545 | / fRowPadSize - 0.5); | |
546 | fRowMax[iplan][1] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick) | |
547 | / fRowPadSize - 0.5); | |
548 | fRowMax[iplan][2] = 1 + TMath::Nint((fClengthI[iplan] - 2. * kCcthick) | |
549 | / fRowPadSize - 0.5); | |
550 | fRowMax[iplan][3] = 1 + TMath::Nint((fClengthM[iplan] - 2. * kCcthick) | |
551 | / fRowPadSize - 0.5); | |
552 | fRowMax[iplan][4] = 1 + TMath::Nint((fClengthO[iplan] - 2. * kCcthick) | |
553 | / fRowPadSize - 0.5); | |
554 | // Columns | |
555 | fColMax[iplan] = 1 + TMath::Nint((fCwidth[iplan] - 2. * kCcthick) | |
556 | / fColPadSize - 0.5); | |
557 | } | |
558 | // Time buckets | |
559 | fTimeMax = 1 + TMath::Nint(kDrThick / fTimeBinSize - 0.5); | |
560 | ||
561 | // The default parameter for the digitization | |
562 | if (!(fGasGain)) fGasGain = 2.0E3; | |
563 | if (!(fNoise)) fNoise = 3000.; | |
564 | if (!(fChipGain)) fChipGain = 10.; | |
565 | if (!(fADCoutRange)) fADCoutRange = 255.; | |
566 | if (!(fADCinRange)) fADCinRange = 2000.; | |
567 | if (!(fADCthreshold)) fADCthreshold = 0; | |
568 | ||
569 | // Transverse and longitudinal diffusion coefficients (Xe/Isobutane) | |
570 | if (!(fDiffusionT)) fDiffusionT = 0.060; | |
571 | if (!(fDiffusionL)) fDiffusionL = 0.017; | |
572 | ||
ab76897d | 573 | printf("**************************************" |
574 | " TRD " | |
575 | "**************************************\n"); | |
99d5402e | 576 | } |
577 | ||
578 | //_____________________________________________________________________________ | |
579 | void AliTRDv2::MakeBranch(Option_t* option) | |
580 | { | |
581 | // | |
582 | // Create Tree branches for the TRD digits. | |
583 | // | |
584 | ||
585 | Int_t buffersize = 4000; | |
586 | Char_t branchname[10]; | |
587 | ||
588 | sprintf(branchname,"%s",GetName()); | |
589 | ||
590 | AliDetector::MakeBranch(option); | |
591 | ||
592 | Char_t *D = strstr(option,"D"); | |
593 | if (fDigits && gAlice->TreeD() && D) { | |
594 | gAlice->TreeD()->Branch(branchname,&fDigits, buffersize); | |
595 | printf("Making Branch %s for digits\n",branchname); | |
596 | } | |
597 | ||
598 | } | |
599 | ||
600 | //_____________________________________________________________________________ | |
601 | Float_t AliTRDv2::PadResponse(Float_t x) | |
602 | { | |
603 | // | |
604 | // The pad response for the chevron pads. | |
605 | // We use a simple Gaussian approximation which should be good | |
606 | // enough for our purpose. | |
607 | // | |
608 | ||
609 | // The parameters for the response function | |
610 | const Float_t aa = 0.8872; | |
611 | const Float_t bb = -0.00573; | |
612 | const Float_t cc = 0.454; | |
613 | const Float_t cc2 = cc*cc; | |
614 | ||
615 | Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2))); | |
616 | ||
617 | //TF1 *funPR = new TF1("funPR","[0]*([1]+exp(-x*x /(2.*[2])))",-3,3); | |
618 | //funPR->SetParameter(0,aa ); | |
619 | //funPR->SetParameter(1,bb ); | |
620 | //funPR->SetParameter(2,cc2); | |
621 | // | |
622 | //Float_t pr = funPR->Eval(distance); | |
623 | // | |
624 | //delete funPR; | |
625 | ||
626 | return (pr); | |
627 | ||
fe4da5cc | 628 | } |
629 | ||
630 | //_____________________________________________________________________________ | |
82bbf98a | 631 | void AliTRDv2::SetSensPlane(Int_t iplane) |
fe4da5cc | 632 | { |
633 | // | |
82bbf98a | 634 | // Defines the hit-sensitive plane (1-6) |
fe4da5cc | 635 | // |
d3f347ff | 636 | |
82bbf98a | 637 | if ((iplane < 0) || (iplane > 6)) { |
638 | printf("Wrong input value: %d\n",iplane); | |
639 | printf("Use standard setting\n"); | |
640 | fSensPlane = 0; | |
641 | fSensSelect = 0; | |
642 | return; | |
643 | } | |
d3f347ff | 644 | |
82bbf98a | 645 | fSensSelect = 1; |
646 | fSensPlane = iplane; | |
d3f347ff | 647 | |
82bbf98a | 648 | } |
d3f347ff | 649 | |
82bbf98a | 650 | //_____________________________________________________________________________ |
651 | void AliTRDv2::SetSensChamber(Int_t ichamber) | |
652 | { | |
653 | // | |
654 | // Defines the hit-sensitive chamber (1-5) | |
655 | // | |
656 | ||
657 | if ((ichamber < 0) || (ichamber > 5)) { | |
658 | printf("Wrong input value: %d\n",ichamber); | |
659 | printf("Use standard setting\n"); | |
660 | fSensChamber = 0; | |
661 | fSensSelect = 0; | |
662 | return; | |
663 | } | |
664 | ||
665 | fSensSelect = 1; | |
666 | fSensChamber = ichamber; | |
667 | ||
668 | } | |
669 | ||
670 | //_____________________________________________________________________________ | |
671 | void AliTRDv2::SetSensSector(Int_t isector) | |
672 | { | |
673 | // | |
674 | // Defines the hit-sensitive sector (1-18) | |
675 | // | |
676 | ||
677 | if ((isector < 0) || (isector > 18)) { | |
678 | printf("Wrong input value: %d\n",isector); | |
679 | printf("Use standard setting\n"); | |
680 | fSensSector = 0; | |
681 | fSensSelect = 0; | |
682 | return; | |
683 | } | |
684 | ||
685 | fSensSelect = 1; | |
686 | fSensSector = isector; | |
687 | ||
688 | } | |
689 | ||
690 | //_____________________________________________________________________________ | |
691 | void AliTRDv2::StepManager() | |
692 | { | |
693 | // | |
694 | // Called at every step in the Transition Radiation Detector version 2. | |
695 | // Slow simulator. Every charged track produces electron cluster as hits | |
696 | // along its path across the drift volume. The step size is set acording | |
697 | // to Bethe-Bloch. The energy distribution of the delta electrons follows | |
698 | // a spectrum taken from Ermilova et al. | |
699 | // | |
700 | ||
701 | Int_t iIdSens, icSens; | |
702 | Int_t iIdSpace, icSpace; | |
703 | Int_t iIdChamber, icChamber; | |
704 | Int_t vol[3]; | |
705 | Int_t iPid; | |
706 | ||
707 | Int_t secMap1[10] = { 3, 7, 8, 9, 10, 11, 2, 1, 18, 17 }; | |
708 | Int_t secMap2[ 5] = { 16, 15, 14, 13, 12 }; | |
709 | Int_t secMap3[ 3] = { 5, 6, 4 }; | |
710 | ||
711 | Float_t hits[4]; | |
712 | Float_t random[1]; | |
713 | Float_t charge; | |
714 | Float_t aMass; | |
0a6d8768 | 715 | |
82bbf98a | 716 | Double_t pTot; |
717 | Double_t qTot; | |
718 | Double_t eDelta; | |
719 | Double_t betaGamma, pp; | |
d3f347ff | 720 | |
82bbf98a | 721 | TLorentzVector pos, mom; |
fe4da5cc | 722 | TClonesArray &lhits = *fHits; |
d3f347ff | 723 | |
82bbf98a | 724 | const Double_t kBig = 1.0E+12; |
725 | ||
fe4da5cc | 726 | // Ionization energy |
d3f347ff | 727 | const Float_t kWion = 22.04; |
fe4da5cc | 728 | // Maximum energy for e+ e- g for the step-size calculation |
d3f347ff | 729 | const Float_t kPTotMax = 0.002; |
fe4da5cc | 730 | // Plateau value of the energy-loss for electron in xenon |
731 | // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253 | |
d3f347ff | 732 | //const Double_t kPlateau = 1.70; |
733 | // the averaged value (26/3/99) | |
734 | const Float_t kPlateau = 1.55; | |
fe4da5cc | 735 | // dN1/dx|min for the gas mixture (90% Xe + 10% CO2) |
d3f347ff | 736 | const Float_t kPrim = 48.0; |
737 | // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2) | |
738 | const Float_t kPoti = 12.1; | |
739 | ||
fe4da5cc | 740 | // Set the maximum step size to a very large number for all |
741 | // neutral particles and those outside the driftvolume | |
cfce8870 | 742 | gMC->SetMaxStep(kBig); |
d3f347ff | 743 | |
fe4da5cc | 744 | // Use only charged tracks |
82bbf98a | 745 | if (( gMC->TrackCharge() ) && |
746 | (!gMC->IsTrackStop() ) && | |
0a6d8768 | 747 | (!gMC->IsTrackDisappeared())) { |
d3f347ff | 748 | |
fe4da5cc | 749 | // Inside a sensitive volume? |
82bbf98a | 750 | iIdSens = gMC->CurrentVolID(icSens); |
751 | if (iIdSens == fIdSens) { | |
752 | ||
753 | iIdSpace = gMC->CurrentVolOffID(4,icSpace ); | |
754 | iIdChamber = gMC->CurrentVolOffID(1,icChamber); | |
d3f347ff | 755 | |
756 | // Calculate the energy of the delta-electrons | |
757 | eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti; | |
758 | eDelta = TMath::Max(eDelta,0.0); | |
759 | ||
760 | // The number of secondary electrons created | |
761 | qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1); | |
762 | ||
82bbf98a | 763 | // The hit coordinates and charge |
764 | gMC->TrackPosition(pos); | |
765 | hits[0] = pos[0]; | |
766 | hits[1] = pos[1]; | |
767 | hits[2] = pos[2]; | |
768 | hits[3] = qTot; | |
d3f347ff | 769 | |
82bbf98a | 770 | // The sector number |
771 | if (iIdSpace == fIdSpace1) | |
772 | vol[0] = secMap1[icSpace-1]; | |
773 | else if (iIdSpace == fIdSpace2) | |
774 | vol[0] = secMap2[icSpace-1]; | |
775 | else if (iIdSpace == fIdSpace3) | |
776 | vol[0] = secMap3[icSpace-1]; | |
777 | ||
778 | // The chamber number | |
d3f347ff | 779 | // 1: outer left |
82bbf98a | 780 | // 2: middle left |
d3f347ff | 781 | // 3: inner |
82bbf98a | 782 | // 4: middle right |
d3f347ff | 783 | // 5: outer right |
82bbf98a | 784 | if (iIdChamber == fIdChamber1) |
785 | vol[1] = (hits[2] < 0 ? 1 : 5); | |
786 | else if (iIdChamber == fIdChamber2) | |
787 | vol[1] = (hits[2] < 0 ? 2 : 4); | |
788 | else if (iIdChamber == fIdChamber3) | |
789 | vol[1] = 3; | |
d3f347ff | 790 | |
82bbf98a | 791 | // The plane number |
792 | vol[2] = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6; | |
fe4da5cc | 793 | |
d3f347ff | 794 | // Check on selected volumes |
795 | Int_t addthishit = 1; | |
796 | if (fSensSelect) { | |
797 | if ((fSensPlane) && (vol[2] != fSensPlane )) addthishit = 0; | |
798 | if ((fSensChamber) && (vol[1] != fSensChamber)) addthishit = 0; | |
799 | if ((fSensSector) && (vol[0] != fSensSector )) addthishit = 0; | |
800 | } | |
801 | ||
82bbf98a | 802 | // Add this hit |
d3f347ff | 803 | if (addthishit) { |
804 | ||
d3f347ff | 805 | new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits); |
806 | ||
807 | // The energy loss according to Bethe Bloch | |
cfce8870 | 808 | gMC->TrackMomentum(mom); |
0a6d8768 | 809 | pTot = mom.Rho(); |
cfce8870 | 810 | iPid = gMC->TrackPid(); |
d3f347ff | 811 | if ( (iPid > 3) || |
812 | ((iPid <= 3) && (pTot < kPTotMax))) { | |
cfce8870 | 813 | aMass = gMC->TrackMass(); |
d3f347ff | 814 | betaGamma = pTot / aMass; |
815 | pp = kPrim * BetheBloch(betaGamma); | |
816 | // Take charge > 1 into account | |
cfce8870 | 817 | charge = gMC->TrackCharge(); |
d3f347ff | 818 | if (TMath::Abs(charge) > 1) pp = pp * charge*charge; |
819 | } | |
820 | // Electrons above 20 Mev/c are at the plateau | |
821 | else { | |
822 | pp = kPrim * kPlateau; | |
823 | } | |
fe4da5cc | 824 | |
d3f347ff | 825 | // Calculate the maximum step size for the next tracking step |
826 | if (pp > 0) { | |
827 | do | |
cfce8870 | 828 | gMC->Rndm(random,1); |
d3f347ff | 829 | while ((random[0] == 1.) || (random[0] == 0.)); |
cfce8870 | 830 | gMC->SetMaxStep( - TMath::Log(random[0]) / pp); |
d3f347ff | 831 | } |
832 | ||
fe4da5cc | 833 | } |
fe4da5cc | 834 | else { |
d3f347ff | 835 | // set step size to maximal value |
cfce8870 | 836 | gMC->SetMaxStep(kBig); |
fe4da5cc | 837 | } |
d3f347ff | 838 | |
fe4da5cc | 839 | } |
d3f347ff | 840 | |
841 | } | |
842 | ||
843 | } | |
844 | ||
845 | //_____________________________________________________________________________ | |
846 | Double_t AliTRDv2::BetheBloch(Double_t bg) | |
847 | { | |
848 | // | |
849 | // Parametrization of the Bethe-Bloch-curve | |
850 | // The parametrization is the same as for the TPC and is taken from Lehrhaus. | |
851 | // | |
852 | ||
d3f347ff | 853 | // This parameters have been adjusted to averaged values from GEANT |
854 | const Double_t kP1 = 7.17960e-02; | |
855 | const Double_t kP2 = 8.54196; | |
856 | const Double_t kP3 = 1.38065e-06; | |
857 | const Double_t kP4 = 5.30972; | |
858 | const Double_t kP5 = 2.83798; | |
859 | ||
99d5402e | 860 | // This parameters have been adjusted to Xe-data found in: |
861 | // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253 | |
862 | //const Double_t kP1 = 0.76176E-1; | |
863 | //const Double_t kP2 = 10.632; | |
864 | //const Double_t kP3 = 3.17983E-6; | |
865 | //const Double_t kP4 = 1.8631; | |
866 | //const Double_t kP5 = 1.9479; | |
867 | ||
d3f347ff | 868 | if (bg > 0) { |
869 | Double_t yy = bg / TMath::Sqrt(1. + bg*bg); | |
870 | Double_t aa = TMath::Power(yy,kP4); | |
871 | Double_t bb = TMath::Power((1./bg),kP5); | |
872 | bb = TMath::Log(kP3 + bb); | |
873 | return ((kP2 - aa - bb)*kP1 / aa); | |
fe4da5cc | 874 | } |
d3f347ff | 875 | else |
876 | return 0; | |
877 | ||
878 | } | |
879 | ||
880 | //_____________________________________________________________________________ | |
6fe53707 | 881 | Double_t Ermilova(Double_t *x, Double_t *) |
d3f347ff | 882 | { |
883 | // | |
99d5402e | 884 | // Calculates the delta-ray energy distribution according to Ermilova. |
d3f347ff | 885 | // Logarithmic scale ! |
886 | // | |
887 | ||
888 | Double_t energy; | |
889 | Double_t dpos; | |
890 | Double_t dnde; | |
891 | ||
892 | Int_t pos1, pos2; | |
893 | ||
894 | const Int_t nV = 31; | |
895 | ||
896 | Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120 | |
897 | , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052 | |
898 | , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915 | |
899 | , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009 | |
900 | , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103 | |
901 | , 9.4727, 9.9035,10.3735,10.5966,10.8198 | |
902 | ,11.5129 }; | |
903 | ||
904 | Float_t vye[nV] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0 | |
905 | , 20.9 , 20.8 , 20.0 , 16.0 , 11.0 | |
906 | , 8.0 , 6.0 , 5.2 , 4.6 , 4.0 | |
907 | , 3.5 , 3.0 , 1.4 , 0.67 , 0.44 | |
908 | , 0.3 , 0.18 , 0.12 , 0.08 , 0.056 | |
909 | , 0.04 , 0.023, 0.015, 0.011, 0.01 | |
910 | , 0.004 }; | |
911 | ||
912 | energy = x[0]; | |
913 | ||
914 | // Find the position | |
915 | pos1 = pos2 = 0; | |
916 | dpos = 0; | |
917 | do { | |
918 | dpos = energy - vxe[pos2++]; | |
919 | } | |
920 | while (dpos > 0); | |
921 | pos2--; | |
922 | if (pos2 > nV) pos2 = nV; | |
923 | pos1 = pos2 - 1; | |
924 | ||
925 | // Differentiate between the sampling points | |
926 | dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]); | |
927 | ||
928 | return dnde; | |
929 | ||
fe4da5cc | 930 | } |