]>
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$ | |
90f8d287 | 18 | Revision 1.15 1999/11/02 17:20:19 fca |
19 | initialise nbytes before using it | |
20 | ||
036da493 | 21 | Revision 1.14 1999/11/02 17:15:54 fca |
22 | Correct ansi scoping not accepted by HP compilers | |
23 | ||
0549c011 | 24 | Revision 1.13 1999/11/02 17:14:51 fca |
25 | Correct ansi scoping not accepted by HP compilers | |
26 | ||
9c767df4 | 27 | Revision 1.12 1999/11/02 16:35:56 fca |
28 | New version of TRD introduced | |
29 | ||
5c7f4665 | 30 | Revision 1.11 1999/11/01 20:41:51 fca |
31 | Added protections against using the wrong version of FRAME | |
32 | ||
ab76897d | 33 | Revision 1.10 1999/09/29 09:24:35 fca |
34 | Introduction of the Copyright and cvs Log | |
35 | ||
4c039060 | 36 | */ |
37 | ||
fe4da5cc | 38 | /////////////////////////////////////////////////////////////////////////////// |
39 | // // | |
5c7f4665 | 40 | // Transition Radiation Detector version 2 -- slow simulator // |
fe4da5cc | 41 | // // |
42 | //Begin_Html | |
43 | /* | |
5c7f4665 | 44 | <img src="picts/AliTRDfullClass.gif"> |
fe4da5cc | 45 | */ |
46 | //End_Html | |
47 | // // | |
48 | // // | |
49 | /////////////////////////////////////////////////////////////////////////////// | |
50 | ||
51 | #include <TMath.h> | |
fe4da5cc | 52 | #include <TVector.h> |
5c7f4665 | 53 | #include <TRandom.h> |
fe4da5cc | 54 | |
fe4da5cc | 55 | #include "AliTRDv1.h" |
5c7f4665 | 56 | #include "AliTRDmatrix.h" |
fe4da5cc | 57 | #include "AliRun.h" |
fe4da5cc | 58 | #include "AliMC.h" |
d3f347ff | 59 | #include "AliConst.h" |
5c7f4665 | 60 | |
fe4da5cc | 61 | ClassImp(AliTRDv1) |
62 | ||
63 | //_____________________________________________________________________________ | |
64 | AliTRDv1::AliTRDv1(const char *name, const char *title) | |
65 | :AliTRD(name, title) | |
66 | { | |
67 | // | |
5c7f4665 | 68 | // Standard constructor for Transition Radiation Detector version 2 |
fe4da5cc | 69 | // |
82bbf98a | 70 | |
5c7f4665 | 71 | fIdSens = 0; |
72 | ||
73 | fIdChamber1 = 0; | |
74 | fIdChamber2 = 0; | |
75 | fIdChamber3 = 0; | |
76 | ||
77 | fSensSelect = 0; | |
78 | fSensPlane = 0; | |
79 | fSensChamber = 0; | |
80 | fSensSector = 0; | |
82bbf98a | 81 | |
5c7f4665 | 82 | fGasGain = 0; |
83 | fNoise = 0; | |
84 | fChipGain = 0; | |
85 | fADCoutRange = 0; | |
86 | fADCinRange = 0; | |
87 | fADCthreshold = 0; | |
88 | ||
89 | fDiffusionT = 0; | |
90 | fDiffusionL = 0; | |
91 | ||
92 | fClusMaxThresh = 0; | |
93 | fClusSigThresh = 0; | |
94 | fClusMethod = 0; | |
95 | ||
96 | fDeltaE = NULL; | |
97 | ||
98 | SetBufferSize(128000); | |
99 | ||
100 | } | |
101 | ||
102 | //_____________________________________________________________________________ | |
103 | AliTRDv1::~AliTRDv1() | |
104 | { | |
82bbf98a | 105 | |
5c7f4665 | 106 | if (fDeltaE) delete fDeltaE; |
82bbf98a | 107 | |
fe4da5cc | 108 | } |
109 | ||
110 | //_____________________________________________________________________________ | |
111 | void AliTRDv1::CreateGeometry() | |
112 | { | |
113 | // | |
5c7f4665 | 114 | // Create the GEANT geometry for the Transition Radiation Detector - Version 2 |
115 | // This version covers the full azimuth. | |
d3f347ff | 116 | // |
117 | ||
82bbf98a | 118 | // Check that FRAME is there otherwise we have no place where to put the TRD |
119 | AliModule* FRAME = gAlice->GetModule("FRAME"); | |
120 | if (!FRAME) return; | |
d3f347ff | 121 | |
82bbf98a | 122 | // Define the chambers |
123 | AliTRD::CreateGeometry(); | |
d3f347ff | 124 | |
fe4da5cc | 125 | } |
126 | ||
127 | //_____________________________________________________________________________ | |
128 | void AliTRDv1::CreateMaterials() | |
129 | { | |
130 | // | |
5c7f4665 | 131 | // Create materials for the Transition Radiation Detector version 2 |
fe4da5cc | 132 | // |
82bbf98a | 133 | |
d3f347ff | 134 | AliTRD::CreateMaterials(); |
82bbf98a | 135 | |
fe4da5cc | 136 | } |
137 | ||
138 | //_____________________________________________________________________________ | |
5c7f4665 | 139 | void AliTRDv1::Diffusion(Float_t driftlength, Float_t *xyz) |
fe4da5cc | 140 | { |
141 | // | |
5c7f4665 | 142 | // Applies the diffusion smearing to the position of a single electron |
fe4da5cc | 143 | // |
82bbf98a | 144 | |
5c7f4665 | 145 | if ((driftlength > 0) && |
146 | (driftlength < kDrThick)) { | |
147 | Float_t driftSqrt = TMath::Sqrt(driftlength); | |
148 | Float_t sigmaT = driftSqrt * fDiffusionT; | |
149 | Float_t sigmaL = driftSqrt * fDiffusionL; | |
150 | xyz[0] = gRandom->Gaus(xyz[0], sigmaL); | |
151 | xyz[1] = gRandom->Gaus(xyz[1], sigmaT); | |
152 | xyz[2] = gRandom->Gaus(xyz[2], sigmaT); | |
153 | } | |
154 | else { | |
155 | xyz[0] = 0.0; | |
156 | xyz[1] = 0.0; | |
157 | xyz[2] = 0.0; | |
158 | } | |
ab76897d | 159 | |
5c7f4665 | 160 | } |
82bbf98a | 161 | |
5c7f4665 | 162 | //_____________________________________________________________________________ |
163 | void AliTRDv1::Hits2Digits() | |
164 | { | |
165 | // | |
166 | // Creates TRD digits from hits. This procedure includes the following: | |
167 | // - Diffusion | |
168 | // - Gas gain including fluctuations | |
169 | // - Pad-response (simple Gaussian approximation) | |
170 | // - Electronics noise | |
171 | // - Electronics gain | |
172 | // - Digitization | |
173 | // - ADC threshold | |
174 | // The corresponding parameter can be adjusted via the various Set-functions. | |
175 | // If these parameters are not explicitly set, default values are used (see | |
176 | // Init-function). | |
177 | // To produce digits from a root-file with TRD-hits use the | |
178 | // slowDigitsCreate.C macro. | |
ab76897d | 179 | // |
5c7f4665 | 180 | |
181 | printf("AliTRDv1::Hits2Digits -- Start creating digits\n"); | |
182 | ||
183 | /////////////////////////////////////////////////////////////// | |
184 | // Parameter | |
185 | /////////////////////////////////////////////////////////////// | |
186 | ||
187 | // Converts number of electrons to fC | |
188 | const Float_t el2fC = 1.602E-19 * 1.0E15; | |
189 | ||
190 | /////////////////////////////////////////////////////////////// | |
191 | ||
192 | Int_t nBytes = 0; | |
193 | ||
9c767df4 | 194 | Int_t iRow; |
195 | ||
5c7f4665 | 196 | AliTRDhit *TRDhit; |
197 | ||
198 | // Get the pointer to the hit tree | |
199 | TTree *HitTree = gAlice->TreeH(); | |
200 | // Get the pointer to the digits tree | |
201 | TTree *DigitsTree = gAlice->TreeD(); | |
202 | ||
203 | // Get the number of entries in the hit tree | |
204 | // (Number of primary particles creating a hit somewhere) | |
205 | Int_t nTrack = (Int_t) HitTree->GetEntries(); | |
206 | ||
207 | Int_t chamBeg = 0; | |
208 | Int_t chamEnd = kNcham; | |
209 | if (fSensChamber) chamEnd = chamBeg = fSensChamber; | |
210 | Int_t planBeg = 0; | |
211 | Int_t planEnd = kNplan; | |
212 | if (fSensPlane) planEnd = planBeg = fSensPlane; | |
213 | Int_t sectBeg = 0; | |
214 | Int_t sectEnd = kNsect; | |
215 | if (fSensSector) sectEnd = sectBeg = fSensSector; | |
216 | ||
217 | // Loop through all the chambers | |
218 | for (Int_t icham = chamBeg; icham < chamEnd; icham++) { | |
219 | for (Int_t iplan = planBeg; iplan < planEnd; iplan++) { | |
220 | for (Int_t isect = sectBeg; isect < sectEnd; isect++) { | |
221 | ||
222 | Int_t nDigits = 0; | |
223 | ||
224 | printf("AliTRDv1::Hits2Digits -- Digitizing chamber %d, plane %d, sector %d\n" | |
225 | ,icham+1,iplan+1,isect+1); | |
226 | ||
227 | // Create a detector matrix to keep the signal and track numbers | |
228 | AliTRDmatrix *matrix = new AliTRDmatrix(fRowMax[iplan][icham][isect] | |
229 | ,fColMax[iplan] | |
230 | ,fTimeMax | |
231 | ,isect+1,icham+1,iplan+1); | |
232 | ||
233 | // Loop through all entries in the tree | |
234 | for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) { | |
235 | ||
236 | gAlice->ResetHits(); | |
237 | nBytes += HitTree->GetEvent(iTrack); | |
238 | ||
239 | // Get the number of hits in the TRD created by this particle | |
240 | Int_t nHit = fHits->GetEntriesFast(); | |
241 | ||
242 | // Loop through the TRD hits | |
243 | for (Int_t iHit = 0; iHit < nHit; iHit++) { | |
244 | ||
245 | if (!(TRDhit = (AliTRDhit *) fHits->UncheckedAt(iHit))) | |
246 | continue; | |
247 | ||
248 | Float_t x = TRDhit->fX; | |
249 | Float_t y = TRDhit->fY; | |
250 | Float_t z = TRDhit->fZ; | |
251 | Float_t q = TRDhit->fQ; | |
252 | Int_t track = TRDhit->fTrack; | |
253 | Int_t plane = TRDhit->fPlane; | |
254 | Int_t sector = TRDhit->fSector; | |
255 | Int_t chamber = TRDhit->fChamber; | |
256 | ||
257 | if ((sector != isect+1) || | |
258 | (plane != iplan+1) || | |
259 | (chamber != icham+1)) | |
260 | continue; | |
261 | ||
262 | // Rotate the sectors on top of each other | |
263 | Float_t phi = 2.0 * kPI / (Float_t) kNsect | |
264 | * ((Float_t) sector - 0.5); | |
265 | Float_t xRot = -x * TMath::Cos(phi) + y * TMath::Sin(phi); | |
266 | Float_t yRot = x * TMath::Sin(phi) + y * TMath::Cos(phi); | |
267 | Float_t zRot = z; | |
268 | ||
269 | // The hit position in pad coordinates (center pad) | |
270 | // The pad row (z-direction) | |
271 | Int_t rowH = (Int_t) ((zRot - fRow0[iplan][icham][isect]) / fRowPadSize); | |
272 | // The pad column (rphi-direction) | |
273 | Int_t colH = (Int_t) ((yRot - fCol0[iplan] ) / fColPadSize); | |
274 | // The time bucket | |
275 | Int_t timeH = (Int_t) ((xRot - fTime0[iplan] ) / fTimeBinSize); | |
276 | ||
277 | // Array to sum up the signal in a box surrounding the | |
278 | // hit postition | |
279 | const Int_t timeBox = 5; | |
280 | const Int_t colBox = 7; | |
281 | const Int_t rowBox = 5; | |
282 | Float_t signalSum[rowBox][colBox][timeBox]; | |
9c767df4 | 283 | for (iRow = 0; iRow < rowBox; iRow++ ) { |
5c7f4665 | 284 | for (Int_t iCol = 0; iCol < colBox; iCol++ ) { |
285 | for (Int_t iTime = 0; iTime < timeBox; iTime++) { | |
286 | signalSum[iRow][iCol][iTime] = 0; | |
287 | } | |
288 | } | |
289 | } | |
290 | ||
291 | // Loop over all electrons of this hit | |
292 | Int_t nEl = (Int_t) q; | |
293 | for (Int_t iEl = 0; iEl < nEl; iEl++) { | |
294 | ||
295 | // Apply the diffusion smearing | |
296 | Float_t driftlength = xRot - fTime0[iplan]; | |
297 | Float_t xyz[3]; | |
298 | xyz[0] = xRot; | |
299 | xyz[1] = yRot; | |
300 | xyz[2] = zRot; | |
301 | Diffusion(driftlength,xyz); | |
302 | ||
303 | // At this point absorption effects that depend on the | |
304 | // driftlength could be taken into account. | |
305 | ||
306 | // The electron position and the distance to the hit position | |
307 | // in pad units | |
308 | // The pad row (z-direction) | |
309 | Int_t rowE = (Int_t) ((xyz[2] - fRow0[iplan][icham][isect]) / fRowPadSize); | |
310 | Int_t rowD = rowH - rowE; | |
311 | // The pad column (rphi-direction) | |
312 | Int_t colE = (Int_t) ((xyz[1] - fCol0[iplan] ) / fColPadSize); | |
313 | Int_t colD = colH - colE; | |
314 | // The time bucket | |
315 | Int_t timeE = (Int_t) ((xyz[0] - fTime0[iplan] ) / fTimeBinSize); | |
316 | Int_t timeD = timeH - timeE; | |
317 | ||
318 | // Apply the gas gain including fluctuations | |
319 | Int_t signal = (Int_t) (-fGasGain * TMath::Log(gRandom->Rndm())); | |
320 | ||
321 | // The distance of the electron to the center of the pad | |
322 | // in units of pad width | |
323 | Float_t dist = (xyz[1] - fCol0[iplan] - (colE + 0.5) * fColPadSize) | |
324 | / fColPadSize; | |
325 | ||
326 | // Sum up the signal in the different pixels | |
327 | // and apply the pad response | |
328 | Int_t rowIdx = rowD + (Int_t) ( rowBox / 2); | |
329 | Int_t colIdx = colD + (Int_t) ( colBox / 2); | |
330 | Int_t timeIdx = timeD + (Int_t) (timeBox / 2); | |
331 | signalSum[rowIdx][colIdx-1][timeIdx] += PadResponse(dist-1.) * signal; | |
332 | signalSum[rowIdx][colIdx ][timeIdx] += PadResponse(dist ) * signal; | |
333 | signalSum[rowIdx][colIdx+1][timeIdx] += PadResponse(dist+1.) * signal; | |
334 | ||
335 | } | |
336 | ||
337 | // Add the padcluster to the detector matrix | |
9c767df4 | 338 | for (iRow = 0; iRow < rowBox; iRow++ ) { |
5c7f4665 | 339 | for (Int_t iCol = 0; iCol < colBox; iCol++ ) { |
340 | for (Int_t iTime = 0; iTime < timeBox; iTime++) { | |
341 | ||
342 | Int_t rowB = rowH + iRow - (Int_t) ( rowBox / 2); | |
343 | Int_t colB = colH + iCol - (Int_t) ( colBox / 2); | |
344 | Int_t timeB = timeH + iTime - (Int_t) (timeBox / 2); | |
345 | ||
346 | Float_t signalB = signalSum[iRow][iCol][iTime]; | |
347 | if (signalB > 0.0) { | |
348 | matrix->AddSignal(rowB,colB,timeB,signalB); | |
349 | if (!(matrix->AddTrack(rowB,colB,timeB,track))) | |
350 | printf(" More than three tracks in a pixel!\n"); | |
351 | } | |
352 | ||
353 | } | |
354 | } | |
355 | } | |
356 | ||
357 | } | |
358 | ||
359 | } | |
360 | ||
361 | // Create the hits for this chamber | |
362 | for (Int_t iRow = 0; iRow < fRowMax[iplan][icham][isect]; iRow++ ) { | |
363 | for (Int_t iCol = 0; iCol < fColMax[iplan] ; iCol++ ) { | |
364 | for (Int_t iTime = 0; iTime < fTimeMax ; iTime++) { | |
365 | ||
366 | Float_t signalAmp = matrix->GetSignal(iRow,iCol,iTime); | |
367 | ||
368 | // Add the noise | |
369 | signalAmp = TMath::Max(gRandom->Gaus(signalAmp,fNoise),(Float_t) 0.0); | |
370 | // Convert to fC | |
371 | signalAmp *= el2fC; | |
372 | // Convert to mV | |
373 | signalAmp *= fChipGain; | |
374 | // Convert to ADC counts | |
375 | Int_t adc = (Int_t) (signalAmp * (fADCoutRange / fADCinRange)); | |
376 | ||
377 | // Apply threshold on ADC value | |
378 | if (adc > fADCthreshold) { | |
379 | ||
380 | Int_t trackSave[3]; | |
381 | for (Int_t ii = 0; ii < 3; ii++) { | |
382 | trackSave[ii] = matrix->GetTrack(iRow,iCol,iTime,ii); | |
383 | } | |
384 | ||
385 | Int_t digits[7]; | |
386 | digits[0] = matrix->GetSector(); | |
387 | digits[1] = matrix->GetChamber(); | |
388 | digits[2] = matrix->GetPlane(); | |
389 | digits[3] = iRow; | |
390 | digits[4] = iCol; | |
391 | digits[5] = iTime; | |
392 | digits[6] = adc; | |
393 | ||
394 | // Add this digit to the TClonesArray | |
395 | AddDigit(trackSave,digits); | |
396 | nDigits++; | |
397 | ||
398 | } | |
399 | ||
400 | } | |
401 | } | |
402 | } | |
403 | ||
404 | printf("AliTRDv1::Hits2Digits -- Number of digits found: %d\n",nDigits); | |
405 | ||
406 | // Clean up | |
407 | delete matrix; | |
408 | ||
409 | } | |
ab76897d | 410 | } |
5c7f4665 | 411 | } |
ab76897d | 412 | |
5c7f4665 | 413 | // Fill the digits-tree |
414 | printf("AliTRDv1::Hits2Digits -- Fill the digits tree\n"); | |
415 | DigitsTree->Fill(); | |
416 | ||
417 | } | |
418 | ||
419 | //_____________________________________________________________________________ | |
420 | void AliTRDv1::Digits2Clusters() | |
421 | { | |
422 | ||
423 | // | |
424 | // Method to convert AliTRDdigits created by AliTRDv1::Hits2Digits() | |
425 | // into AliTRDclusters | |
426 | // To produce cluster from a root-file with TRD-digits use the | |
427 | // slowClusterCreate.C macro. | |
428 | // | |
9c767df4 | 429 | |
0549c011 | 430 | Int_t row; |
5c7f4665 | 431 | |
432 | printf("AliTRDv1::Digits2Clusters -- Start creating clusters\n"); | |
433 | ||
434 | AliTRDdigit *TRDdigit; | |
435 | TClonesArray *TRDDigits; | |
436 | ||
437 | // Parameters | |
438 | Float_t maxThresh = fClusMaxThresh; // threshold value for maximum | |
439 | Float_t signalThresh = fClusSigThresh; // threshold value for digit signal | |
440 | Int_t clusteringMethod = fClusMethod; // clustering method option (for testing) | |
441 | ||
442 | const Float_t epsilon = 0.01; // iteration limit for unfolding procedure | |
443 | ||
444 | // Get the pointer to the digits tree | |
445 | TTree *DigitTree = gAlice->TreeD(); | |
446 | // Get the pointer to the cluster tree | |
447 | TTree *ClusterTree = gAlice->TreeD(); | |
448 | ||
449 | // Get the pointer to the digits container | |
450 | TRDDigits = Digits(); | |
451 | ||
452 | Int_t chamBeg = 0; | |
453 | Int_t chamEnd = kNcham; | |
454 | if (fSensChamber) chamEnd = chamBeg = fSensChamber; | |
455 | Int_t planBeg = 0; | |
456 | Int_t planEnd = kNplan; | |
457 | if (fSensPlane) planEnd = planBeg = fSensPlane; | |
458 | Int_t sectBeg = 0; | |
459 | Int_t sectEnd = kNsect; | |
460 | if (fSensSector) sectEnd = sectBeg = fSensSector; | |
461 | ||
462 | // Import the digit tree | |
463 | gAlice->ResetDigits(); | |
036da493 | 464 | Int_t nbytes=0; |
5c7f4665 | 465 | nbytes += DigitTree->GetEvent(1); |
466 | ||
467 | // Get the number of digits in the detector | |
468 | Int_t nTRDDigits = TRDDigits->GetEntriesFast(); | |
469 | ||
470 | // *** Start clustering *** in every chamber | |
471 | for (Int_t icham = chamBeg; icham < chamEnd; icham++) { | |
472 | for (Int_t iplan = planBeg; iplan < planEnd; iplan++) { | |
473 | for (Int_t isect = sectBeg; isect < sectEnd; isect++) { | |
474 | ||
475 | Int_t nClusters = 0; | |
476 | printf("AliTRDv1::Digits2Clusters -- Finding clusters in chamber %d, plane %d, sector %d\n" | |
477 | ,icham+1,iplan+1,isect+1); | |
478 | ||
479 | // Create a detector matrix to keep maxima | |
480 | AliTRDmatrix *digitMatrix = new AliTRDmatrix(fRowMax[iplan][icham][isect] | |
481 | ,fColMax[iplan] | |
482 | ,fTimeMax,isect+1 | |
483 | ,icham+1,iplan+1); | |
484 | // Create a matrix to contain maximum flags | |
485 | AliTRDmatrix *maximaMatrix = new AliTRDmatrix(fRowMax[iplan][icham][isect] | |
486 | ,fColMax[iplan] | |
487 | ,fTimeMax | |
488 | ,isect+1,icham+1,iplan+1); | |
489 | ||
490 | // Loop through all TRD digits | |
491 | for (Int_t iTRDDigits = 0; iTRDDigits < nTRDDigits; iTRDDigits++) { | |
492 | ||
493 | // Get the information for this digit | |
494 | TRDdigit = (AliTRDdigit*) TRDDigits->UncheckedAt(iTRDDigits); | |
495 | Int_t signal = TRDdigit->fSignal; | |
496 | Int_t sector = TRDdigit->fSector; | |
497 | Int_t chamber = TRDdigit->fChamber; | |
498 | Int_t plane = TRDdigit->fPlane; | |
499 | Int_t row = TRDdigit->fRow; | |
500 | Int_t col = TRDdigit->fCol; | |
501 | Int_t time = TRDdigit->fTime; | |
502 | ||
503 | Int_t track[3]; | |
504 | for (Int_t iTrack = 0; iTrack < 3; iTrack++) { | |
505 | track[iTrack] = TRDdigit->AliDigit::fTracks[iTrack]; | |
506 | } | |
507 | ||
508 | if ((sector != isect+1) || | |
509 | (plane != iplan+1) || | |
510 | (chamber != icham+1)) | |
511 | continue; | |
512 | ||
513 | // Fill the detector matrix | |
514 | if (signal > signalThresh) { | |
515 | digitMatrix->SetSignal(row,col,time,signal); | |
516 | for (Int_t iTrack = 0; iTrack < 3; iTrack++) { | |
517 | if (track[iTrack] > 0) { | |
518 | digitMatrix->AddTrack(row,col,time,track[iTrack]); | |
519 | } | |
520 | } | |
521 | } | |
522 | ||
523 | } | |
524 | ||
525 | // Loop chamber and find maxima in digitMatrix | |
9c767df4 | 526 | for (row = 0; row < fRowMax[iplan][icham][isect]; row++) { |
5c7f4665 | 527 | for (Int_t col = 1; col < fColMax[iplan] ; col++) { |
528 | for (Int_t time = 0; time < fTimeMax ; time++) { | |
529 | ||
530 | if (digitMatrix->GetSignal(row,col,time) | |
531 | < digitMatrix->GetSignal(row,col - 1,time)) { | |
532 | // really maximum? | |
533 | if (col > 1) { | |
534 | if (digitMatrix->GetSignal(row,col - 2,time) | |
535 | < digitMatrix->GetSignal(row,col - 1,time)) { | |
536 | // yes, so set maximum flag | |
537 | maximaMatrix->SetSignal(row,col - 1,time,1); | |
538 | } | |
539 | else maximaMatrix->SetSignal(row,col - 1,time,0); | |
540 | } | |
541 | } | |
542 | ||
543 | } // time | |
544 | } // col | |
545 | } // row | |
546 | ||
547 | // now check maxima and calculate cluster position | |
9c767df4 | 548 | for (row = 0; row < fRowMax[iplan][icham][isect]; row++) { |
5c7f4665 | 549 | for (Int_t col = 1; col < fColMax[iplan] ; col++) { |
550 | for (Int_t time = 0; time < fTimeMax ; time++) { | |
551 | ||
552 | if ((maximaMatrix->GetSignal(row,col,time) > 0) | |
553 | && (digitMatrix->GetSignal(row,col,time) > maxThresh)) { | |
554 | ||
555 | Int_t clusters[5] = {0}; // cluster-object data | |
556 | ||
557 | Float_t ratio = 0; // ratio resulting from unfolding | |
558 | Float_t padSignal[5] = {0}; // signals on max and neighbouring pads | |
559 | Float_t clusterSignal[3] = {0}; // signals from cluster | |
560 | Float_t clusterPos[3] = {0}; // cluster in ALICE refFrame coords | |
561 | Float_t clusterPads[6] = {0}; // cluster pad info | |
562 | ||
563 | // setting values | |
564 | clusters[0] = isect+1; // = isect ???? | |
565 | clusters[1] = icham+1; // = ichamber ???? | |
566 | clusters[2] = iplan+1; // = iplane ???? | |
567 | clusters[3] = time; | |
568 | ||
569 | clusterPads[0] = icham+1; | |
570 | clusterPads[1] = isect+1; | |
571 | clusterPads[2] = iplan+1; | |
572 | ||
573 | for (Int_t iPad = 0; iPad < 3; iPad++) { | |
574 | clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time); | |
575 | } | |
576 | ||
577 | // neighbouring maximum on right side? | |
578 | if (col < fColMax[iplan] - 2) { | |
579 | if (maximaMatrix->GetSignal(row,col + 2,time) > 0) { | |
580 | for (Int_t iPad = 0; iPad < 5; iPad++) { | |
581 | padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time); | |
582 | } | |
583 | ||
584 | // unfold: | |
585 | ratio = Unfold(epsilon, padSignal); | |
586 | ||
587 | // set signal on overlapping pad to ratio | |
588 | clusterSignal[2] *= ratio; | |
589 | } | |
590 | } | |
591 | ||
592 | switch (clusteringMethod) { | |
593 | case 1: | |
594 | // method 1: simply center of mass | |
595 | clusterPads[3] = row + 0.5; | |
596 | clusterPads[4] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) / | |
597 | (clusterSignal[1] + clusterSignal[2] + clusterSignal[3]); | |
598 | clusterPads[5] = time + 0.5; | |
599 | ||
600 | nClusters++; | |
601 | break; | |
602 | case 2: | |
603 | // method 2: integral gauss fit on 3 pads | |
604 | TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads" | |
605 | , 5, -1.5, 3.5); | |
606 | for (Int_t iCol = -1; iCol <= 3; iCol++) { | |
607 | if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1; | |
608 | hPadCharges->Fill(iCol, clusterSignal[iCol]); | |
609 | } | |
610 | hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5); | |
611 | TF1 *fPadChargeFit = hPadCharges->GetFunction("gaus"); | |
612 | Double_t colMean = fPadChargeFit->GetParameter(1); | |
613 | ||
614 | clusterPads[3] = row + 0.5; | |
615 | clusterPads[4] = col - 1.5 + colMean; | |
616 | clusterPads[5] = time + 0.5; | |
617 | ||
618 | delete hPadCharges; | |
619 | ||
620 | nClusters++; | |
621 | break; | |
622 | } | |
623 | ||
624 | Float_t clusterCharge = clusterSignal[0] | |
625 | + clusterSignal[1] | |
626 | + clusterSignal[2]; | |
627 | clusters[4] = (Int_t)clusterCharge; | |
628 | ||
629 | Int_t trackSave[3]; | |
630 | for (Int_t iTrack = 0; iTrack < 3; iTrack++) { | |
631 | trackSave[iTrack] = digitMatrix->GetTrack(row,col,time,iTrack); | |
632 | } | |
633 | ||
634 | // Calculate cluster position in ALICE refFrame coords | |
635 | // and set array clusterPos to calculated values | |
636 | Pads2XYZ(clusterPads, clusterPos); | |
637 | ||
638 | // Add cluster to reconstruction tree | |
639 | AddCluster(trackSave,clusters,clusterPos); | |
640 | ||
641 | } | |
82bbf98a | 642 | |
5c7f4665 | 643 | } // time |
644 | } // col | |
645 | } // row | |
646 | ||
647 | printf("AliTRDv1::Digits2Clusters -- Number of clusters found: %d\n",nClusters); | |
648 | ||
649 | delete digitMatrix; | |
650 | delete maximaMatrix; | |
651 | ||
652 | } // isect | |
653 | } // iplan | |
654 | } // icham | |
655 | ||
656 | // Fill the cluster-tree | |
657 | printf("AliTRDv1::Digits2Clusters -- Total number of clusters found: %d\n" | |
658 | ,fClusters->GetEntries()); | |
659 | printf("AliTRDv1::Digits2Clusters -- Fill the cluster tree\n"); | |
660 | ClusterTree->Fill(); | |
661 | ||
662 | } | |
663 | ||
664 | //_____________________________________________________________________________ | |
665 | void AliTRDv1::Init() | |
666 | { | |
667 | // | |
668 | // Initialise Transition Radiation Detector after geometry has been built. | |
669 | // Includes the default settings of all parameter for the digitization. | |
670 | // | |
671 | ||
672 | AliTRD::Init(); | |
673 | ||
674 | printf(" Slow simulator\n"); | |
675 | if (fSensPlane) | |
676 | printf(" Only plane %d is sensitive\n",fSensPlane); | |
677 | if (fSensChamber) | |
678 | printf(" Only chamber %d is sensitive\n",fSensChamber); | |
679 | if (fSensSector) | |
680 | printf(" Only sector %d is sensitive\n",fSensSector); | |
681 | ||
682 | // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2) | |
683 | const Float_t kPoti = 12.1; | |
684 | // Maximum energy (50 keV); | |
685 | const Float_t kEend = 50000.0; | |
686 | // Ermilova distribution for the delta-ray spectrum | |
687 | Float_t Poti = TMath::Log(kPoti); | |
688 | Float_t Eend = TMath::Log(kEend); | |
689 | fDeltaE = new TF1("deltae",Ermilova,Poti,Eend,0); | |
690 | ||
691 | // Identifier of the sensitive volume (drift region) | |
692 | fIdSens = gMC->VolId("UL05"); | |
82bbf98a | 693 | |
694 | // Identifier of the TRD-driftchambers | |
695 | fIdChamber1 = gMC->VolId("UCIO"); | |
696 | fIdChamber2 = gMC->VolId("UCIM"); | |
697 | fIdChamber3 = gMC->VolId("UCII"); | |
698 | ||
5c7f4665 | 699 | // The default parameter for the digitization |
700 | if (!(fGasGain)) fGasGain = 2.0E3; | |
701 | if (!(fNoise)) fNoise = 3000.; | |
702 | if (!(fChipGain)) fChipGain = 10.; | |
703 | if (!(fADCoutRange)) fADCoutRange = 255.; | |
704 | if (!(fADCinRange)) fADCinRange = 2000.; | |
705 | if (!(fADCthreshold)) fADCthreshold = 1; | |
706 | ||
707 | // Transverse and longitudinal diffusion coefficients (Xe/Isobutane) | |
708 | if (!(fDiffusionT)) fDiffusionT = 0.060; | |
709 | if (!(fDiffusionL)) fDiffusionL = 0.017; | |
710 | ||
711 | // The default parameter for the clustering | |
712 | if (!(fClusMaxThresh)) fClusMaxThresh = 5.0; | |
713 | if (!(fClusSigThresh)) fClusSigThresh = 2.0; | |
714 | if (!(fClusMethod)) fClusMethod = 1; | |
715 | ||
716 | for (Int_t i = 0; i < 80; i++) printf("*"); | |
717 | printf("\n"); | |
718 | ||
fe4da5cc | 719 | } |
720 | ||
721 | //_____________________________________________________________________________ | |
5c7f4665 | 722 | Float_t AliTRDv1::PadResponse(Float_t x) |
fe4da5cc | 723 | { |
724 | // | |
5c7f4665 | 725 | // The pad response for the chevron pads. |
726 | // We use a simple Gaussian approximation which should be good | |
727 | // enough for our purpose. | |
fe4da5cc | 728 | // |
d3f347ff | 729 | |
5c7f4665 | 730 | // The parameters for the response function |
731 | const Float_t aa = 0.8872; | |
732 | const Float_t bb = -0.00573; | |
733 | const Float_t cc = 0.454; | |
734 | const Float_t cc2 = cc*cc; | |
735 | ||
736 | Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2))); | |
737 | ||
738 | return (pr); | |
739 | ||
740 | } | |
741 | ||
742 | //_____________________________________________________________________________ | |
743 | void AliTRDv1::SetSensPlane(Int_t iplane) | |
744 | { | |
745 | // | |
746 | // Defines the hit-sensitive plane (1-6) | |
747 | // | |
82bbf98a | 748 | |
5c7f4665 | 749 | if ((iplane < 0) || (iplane > 6)) { |
750 | printf("Wrong input value: %d\n",iplane); | |
751 | printf("Use standard setting\n"); | |
752 | fSensPlane = 0; | |
753 | fSensSelect = 0; | |
754 | return; | |
755 | } | |
82bbf98a | 756 | |
5c7f4665 | 757 | fSensSelect = 1; |
758 | fSensPlane = iplane; | |
82bbf98a | 759 | |
5c7f4665 | 760 | } |
761 | ||
762 | //_____________________________________________________________________________ | |
763 | void AliTRDv1::SetSensChamber(Int_t ichamber) | |
764 | { | |
765 | // | |
766 | // Defines the hit-sensitive chamber (1-5) | |
767 | // | |
768 | ||
769 | if ((ichamber < 0) || (ichamber > 5)) { | |
770 | printf("Wrong input value: %d\n",ichamber); | |
771 | printf("Use standard setting\n"); | |
772 | fSensChamber = 0; | |
773 | fSensSelect = 0; | |
774 | return; | |
775 | } | |
776 | ||
777 | fSensSelect = 1; | |
778 | fSensChamber = ichamber; | |
779 | ||
780 | } | |
781 | ||
782 | //_____________________________________________________________________________ | |
783 | void AliTRDv1::SetSensSector(Int_t isector) | |
784 | { | |
785 | // | |
786 | // Defines the hit-sensitive sector (1-18) | |
787 | // | |
788 | ||
789 | if ((isector < 0) || (isector > 18)) { | |
790 | printf("Wrong input value: %d\n",isector); | |
791 | printf("Use standard setting\n"); | |
792 | fSensSector = 0; | |
793 | fSensSelect = 0; | |
794 | return; | |
795 | } | |
796 | ||
797 | fSensSelect = 1; | |
798 | fSensSector = isector; | |
799 | ||
800 | } | |
801 | ||
802 | //_____________________________________________________________________________ | |
803 | void AliTRDv1::StepManager() | |
804 | { | |
805 | // | |
806 | // Called at every step in the Transition Radiation Detector version 2. | |
807 | // Slow simulator. Every charged track produces electron cluster as hits | |
808 | // along its path across the drift volume. The step size is set acording | |
809 | // to Bethe-Bloch. The energy distribution of the delta electrons follows | |
810 | // a spectrum taken from Ermilova et al. | |
811 | // | |
812 | ||
813 | Int_t iIdSens, icSens; | |
814 | Int_t iIdSpace, icSpace; | |
815 | Int_t iIdChamber, icChamber; | |
816 | Int_t vol[3]; | |
817 | Int_t iPid; | |
818 | ||
819 | Float_t hits[4]; | |
820 | Float_t random[1]; | |
821 | Float_t charge; | |
822 | Float_t aMass; | |
823 | ||
824 | Double_t pTot; | |
825 | Double_t qTot; | |
826 | Double_t eDelta; | |
827 | Double_t betaGamma, pp; | |
828 | ||
829 | TLorentzVector pos, mom; | |
82bbf98a | 830 | TClonesArray &lhits = *fHits; |
831 | ||
5c7f4665 | 832 | const Double_t kBig = 1.0E+12; |
833 | ||
834 | // Ionization energy | |
835 | const Float_t kWion = 22.04; | |
836 | // Maximum energy for e+ e- g for the step-size calculation | |
837 | const Float_t kPTotMax = 0.002; | |
838 | // Plateau value of the energy-loss for electron in xenon | |
839 | // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253 | |
840 | //const Double_t kPlateau = 1.70; | |
841 | // the averaged value (26/3/99) | |
842 | const Float_t kPlateau = 1.55; | |
843 | // dN1/dx|min for the gas mixture (90% Xe + 10% CO2) | |
844 | const Float_t kPrim = 48.0; | |
845 | // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2) | |
846 | const Float_t kPoti = 12.1; | |
847 | ||
848 | // Set the maximum step size to a very large number for all | |
849 | // neutral particles and those outside the driftvolume | |
850 | gMC->SetMaxStep(kBig); | |
851 | ||
852 | // Use only charged tracks | |
853 | if (( gMC->TrackCharge() ) && | |
854 | (!gMC->IsTrackStop() ) && | |
855 | (!gMC->IsTrackDisappeared())) { | |
fe4da5cc | 856 | |
5c7f4665 | 857 | // Inside a sensitive volume? |
82bbf98a | 858 | iIdSens = gMC->CurrentVolID(icSens); |
859 | if (iIdSens == fIdSens) { | |
860 | ||
82bbf98a | 861 | iIdSpace = gMC->CurrentVolOffID(4,icSpace ); |
862 | iIdChamber = gMC->CurrentVolOffID(1,icChamber); | |
fe4da5cc | 863 | |
5c7f4665 | 864 | // Calculate the energy of the delta-electrons |
865 | eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti; | |
866 | eDelta = TMath::Max(eDelta,0.0); | |
867 | ||
868 | // The number of secondary electrons created | |
869 | qTot = (Double_t) ((Int_t) (eDelta / kWion) + 1); | |
870 | ||
871 | // The hit coordinates and charge | |
872 | gMC->TrackPosition(pos); | |
873 | hits[0] = pos[0]; | |
874 | hits[1] = pos[1]; | |
875 | hits[2] = pos[2]; | |
876 | hits[3] = qTot; | |
877 | ||
fe4da5cc | 878 | // The sector number |
90f8d287 | 879 | Float_t phi = pos[1] != 0 ? kRaddeg*TMath::ATan2(pos[0],pos[1]) : (pos[0] > 0 ? 180. : 0.); |
5c7f4665 | 880 | vol[0] = ((Int_t) (phi / 20)) + 1; |
82bbf98a | 881 | |
d3f347ff | 882 | // The chamber number |
883 | // 1: outer left | |
82bbf98a | 884 | // 2: middle left |
d3f347ff | 885 | // 3: inner |
82bbf98a | 886 | // 4: middle right |
d3f347ff | 887 | // 5: outer right |
82bbf98a | 888 | if (iIdChamber == fIdChamber1) |
889 | vol[1] = (hits[2] < 0 ? 1 : 5); | |
890 | else if (iIdChamber == fIdChamber2) | |
891 | vol[1] = (hits[2] < 0 ? 2 : 4); | |
892 | else if (iIdChamber == fIdChamber3) | |
d3f347ff | 893 | vol[1] = 3; |
82bbf98a | 894 | |
fe4da5cc | 895 | // The plane number |
82bbf98a | 896 | vol[2] = icChamber - TMath::Nint((Float_t) (icChamber / 7)) * 6; |
897 | ||
5c7f4665 | 898 | // Check on selected volumes |
899 | Int_t addthishit = 1; | |
900 | if (fSensSelect) { | |
901 | if ((fSensPlane) && (vol[2] != fSensPlane )) addthishit = 0; | |
902 | if ((fSensChamber) && (vol[1] != fSensChamber)) addthishit = 0; | |
903 | if ((fSensSector) && (vol[0] != fSensSector )) addthishit = 0; | |
904 | } | |
905 | ||
906 | // Add this hit | |
907 | if (addthishit) { | |
908 | ||
909 | new(lhits[fNhits++]) AliTRDhit(fIshunt,gAlice->CurrentTrack(),vol,hits); | |
910 | ||
911 | // The energy loss according to Bethe Bloch | |
912 | gMC->TrackMomentum(mom); | |
913 | pTot = mom.Rho(); | |
914 | iPid = gMC->TrackPid(); | |
915 | if ( (iPid > 3) || | |
916 | ((iPid <= 3) && (pTot < kPTotMax))) { | |
917 | aMass = gMC->TrackMass(); | |
918 | betaGamma = pTot / aMass; | |
919 | pp = kPrim * BetheBloch(betaGamma); | |
920 | // Take charge > 1 into account | |
921 | charge = gMC->TrackCharge(); | |
922 | if (TMath::Abs(charge) > 1) pp = pp * charge*charge; | |
923 | } | |
924 | // Electrons above 20 Mev/c are at the plateau | |
925 | else { | |
926 | pp = kPrim * kPlateau; | |
927 | } | |
928 | ||
929 | // Calculate the maximum step size for the next tracking step | |
930 | if (pp > 0) { | |
931 | do | |
932 | gMC->Rndm(random,1); | |
933 | while ((random[0] == 1.) || (random[0] == 0.)); | |
934 | gMC->SetMaxStep( - TMath::Log(random[0]) / pp); | |
935 | } | |
936 | ||
937 | } | |
938 | else { | |
939 | // set step size to maximal value | |
940 | gMC->SetMaxStep(kBig); | |
941 | } | |
d3f347ff | 942 | |
943 | } | |
944 | ||
5c7f4665 | 945 | } |
946 | ||
947 | } | |
948 | ||
949 | //_____________________________________________________________________________ | |
950 | Double_t AliTRDv1::BetheBloch(Double_t bg) | |
951 | { | |
952 | // | |
953 | // Parametrization of the Bethe-Bloch-curve | |
954 | // The parametrization is the same as for the TPC and is taken from Lehrhaus. | |
955 | // | |
956 | ||
957 | // This parameters have been adjusted to averaged values from GEANT | |
958 | const Double_t kP1 = 7.17960e-02; | |
959 | const Double_t kP2 = 8.54196; | |
960 | const Double_t kP3 = 1.38065e-06; | |
961 | const Double_t kP4 = 5.30972; | |
962 | const Double_t kP5 = 2.83798; | |
963 | ||
964 | // This parameters have been adjusted to Xe-data found in: | |
965 | // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253 | |
966 | //const Double_t kP1 = 0.76176E-1; | |
967 | //const Double_t kP2 = 10.632; | |
968 | //const Double_t kP3 = 3.17983E-6; | |
969 | //const Double_t kP4 = 1.8631; | |
970 | //const Double_t kP5 = 1.9479; | |
971 | ||
972 | if (bg > 0) { | |
973 | Double_t yy = bg / TMath::Sqrt(1. + bg*bg); | |
974 | Double_t aa = TMath::Power(yy,kP4); | |
975 | Double_t bb = TMath::Power((1./bg),kP5); | |
976 | bb = TMath::Log(kP3 + bb); | |
977 | return ((kP2 - aa - bb)*kP1 / aa); | |
978 | } | |
979 | else | |
980 | return 0; | |
d3f347ff | 981 | |
fe4da5cc | 982 | } |
5c7f4665 | 983 | |
984 | //_____________________________________________________________________________ | |
985 | Double_t Ermilova(Double_t *x, Double_t *) | |
986 | { | |
987 | // | |
988 | // Calculates the delta-ray energy distribution according to Ermilova. | |
989 | // Logarithmic scale ! | |
990 | // | |
991 | ||
992 | Double_t energy; | |
993 | Double_t dpos; | |
994 | Double_t dnde; | |
995 | ||
996 | Int_t pos1, pos2; | |
997 | ||
998 | const Int_t nV = 31; | |
999 | ||
1000 | Float_t vxe[nV] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120 | |
1001 | , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052 | |
1002 | , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915 | |
1003 | , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009 | |
1004 | , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103 | |
1005 | , 9.4727, 9.9035,10.3735,10.5966,10.8198 | |
1006 | ,11.5129 }; | |
1007 | ||
1008 | Float_t vye[nV] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0 | |
1009 | , 20.9 , 20.8 , 20.0 , 16.0 , 11.0 | |
1010 | , 8.0 , 6.0 , 5.2 , 4.6 , 4.0 | |
1011 | , 3.5 , 3.0 , 1.4 , 0.67 , 0.44 | |
1012 | , 0.3 , 0.18 , 0.12 , 0.08 , 0.056 | |
1013 | , 0.04 , 0.023, 0.015, 0.011, 0.01 | |
1014 | , 0.004 }; | |
1015 | ||
1016 | energy = x[0]; | |
1017 | ||
1018 | // Find the position | |
1019 | pos1 = pos2 = 0; | |
1020 | dpos = 0; | |
1021 | do { | |
1022 | dpos = energy - vxe[pos2++]; | |
1023 | } | |
1024 | while (dpos > 0); | |
1025 | pos2--; | |
1026 | if (pos2 > nV) pos2 = nV; | |
1027 | pos1 = pos2 - 1; | |
1028 | ||
1029 | // Differentiate between the sampling points | |
1030 | dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]); | |
1031 | ||
1032 | return dnde; | |
1033 | ||
1034 | } | |
1035 | ||
1036 | //_____________________________________________________________________________ | |
1037 | void AliTRDv1::Pads2XYZ(Float_t *pads, Float_t *pos) | |
1038 | { | |
1039 | // Method to convert pad coordinates (row,col,time) | |
1040 | // into ALICE reference frame coordinates (x,y,z) | |
1041 | ||
1042 | Int_t chamber = (Int_t) pads[0]; // chamber info (1-5) | |
1043 | Int_t sector = (Int_t) pads[1]; // sector info (1-18) | |
1044 | Int_t plane = (Int_t) pads[2]; // plane info (1-6) | |
1045 | ||
1046 | Int_t icham = chamber - 1; // chamber info (0-4) | |
1047 | Int_t isect = sector - 1; // sector info (0-17) | |
1048 | Int_t iplan = plane - 1; // plane info (0-5) | |
1049 | ||
1050 | Float_t padRow = pads[3]; // Pad Row position | |
1051 | Float_t padCol = pads[4]; // Pad Column position | |
1052 | Float_t timeSlice = pads[5]; // Time "position" | |
1053 | ||
1054 | // calculate (x,y) position in rotated chamber | |
1055 | Float_t yRot = fCol0[iplan] + padCol * fColPadSize; | |
1056 | Float_t xRot = fTime0[iplan] + timeSlice * fTimeBinSize; | |
1057 | // calculate z-position: | |
1058 | Float_t z = fRow0[iplan][icham][isect] + padRow * fRowPadSize; | |
1059 | ||
1060 | /** | |
1061 | rotate chamber back to original position | |
1062 | 1. mirror at y-axis, 2. rotate back to position (-phi) | |
1063 | / cos(phi) -sin(phi) \ / -1 0 \ / -cos(phi) -sin(phi) \ | |
1064 | \ sin(phi) cos(phi) / * \ 0 1 / = \ -sin(phi) cos(phi) / | |
1065 | **/ | |
1066 | //Float_t phi = 2*kPI / kNsect * ((Float_t) sector - 0.5); | |
1067 | //Float_t x = -xRot * TMath::Cos(phi) - yRot * TMath::Sin(phi); | |
1068 | //Float_t y = -xRot * TMath::Sin(phi) + yRot * TMath::Cos(phi); | |
1069 | Float_t phi = 2*kPI / kNsect * ((Float_t) sector - 0.5); | |
1070 | Float_t x = -xRot * TMath::Cos(phi) + yRot * TMath::Sin(phi); | |
1071 | Float_t y = xRot * TMath::Sin(phi) + yRot * TMath::Cos(phi); | |
1072 | ||
1073 | // Setting values | |
1074 | pos[0] = x; | |
1075 | pos[1] = y; | |
1076 | pos[2] = z; | |
1077 | ||
1078 | } | |
1079 | ||
1080 | //_____________________________________________________________________________ | |
1081 | Float_t AliTRDv1::Unfold(Float_t eps, Float_t* padSignal) | |
1082 | { | |
1083 | // Method to unfold neighbouring maxima. | |
1084 | // The charge ratio on the overlapping pad is calculated | |
1085 | // until there is no more change within the range given by eps. | |
1086 | // The resulting ratio is then returned to the calling method. | |
1087 | ||
1088 | Int_t itStep = 0; // count iteration steps | |
1089 | ||
1090 | Float_t ratio = 0.5; // start value for ratio | |
1091 | Float_t prevRatio = 0; // store previous ratio | |
1092 | ||
1093 | Float_t newLeftSignal[3] = {0}; // array to store left cluster signal | |
1094 | Float_t newRightSignal[3] = {0}; // array to store right cluster signal | |
1095 | ||
1096 | // start iteration: | |
1097 | while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { | |
1098 | ||
1099 | itStep++; | |
1100 | prevRatio = ratio; | |
1101 | ||
1102 | // cluster position according to charge ratio | |
1103 | Float_t maxLeft = (ratio*padSignal[2] - padSignal[0]) / | |
1104 | (padSignal[0] + padSignal[1] + ratio*padSignal[2]); | |
1105 | Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) / | |
1106 | ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]); | |
1107 | ||
1108 | // set cluster charge ratio | |
1109 | Float_t ampLeft = padSignal[1]; | |
1110 | Float_t ampRight = padSignal[3]; | |
1111 | ||
1112 | // apply pad response to parameters | |
1113 | newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft); | |
1114 | newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft); | |
1115 | newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft); | |
1116 | ||
1117 | newRightSignal[0] = ampRight*PadResponse(-1 - maxRight); | |
1118 | newRightSignal[1] = ampRight*PadResponse( 0 - maxRight); | |
1119 | newRightSignal[2] = ampRight*PadResponse( 1 - maxRight); | |
1120 | ||
1121 | // calculate new overlapping ratio | |
1122 | ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]); | |
1123 | ||
1124 | } | |
1125 | ||
1126 | return ratio; | |
1127 | ||
1128 | } | |
1129 |