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