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