]>
Commit | Line | Data |
---|---|---|
ca21baaa | 1 | |
2 | /************************************************************************** | |
3 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
4 | * * | |
5 | * Author: The ALICE Off-line Project. * | |
6 | * Contributors are mentioned in the code where appropriate. * | |
7 | * * | |
8 | * Permission to use, copy, modify and distribute this software and its * | |
9 | * documentation strictly for non-commercial purposes is hereby granted * | |
10 | * without fee, provided that the above copyright notice appears in all * | |
11 | * copies and that both the copyright notice and this permission notice * | |
12 | * appear in the supporting documentation. The authors make no claims * | |
13 | * about the suitability of this software for any purpose. It is * | |
14 | * provided "as is" without express or implied warranty. * | |
15 | **************************************************************************/ | |
16 | ||
17 | /* $Id$ */ | |
18 | ||
19 | /////////////////////////////////////////////////////////////////////////////// | |
20 | // // | |
21 | // TRD cluster finder // | |
22 | // // | |
23 | /////////////////////////////////////////////////////////////////////////////// | |
24 | ||
25 | #include <TF1.h> | |
26 | #include <TTree.h> | |
27 | #include <TH1.h> | |
28 | #include <TFile.h> | |
29 | ||
30 | #include "AliRunLoader.h" | |
31 | #include "AliLoader.h" | |
32 | #include "AliRawReader.h" | |
33 | #include "AliLog.h" | |
34 | #include "AliAlignObj.h" | |
35 | ||
36 | #include "AliTRDclusterizerV2.h" | |
37 | #include "AliTRDgeometry.h" | |
38 | #include "AliTRDdataArrayF.h" | |
39 | #include "AliTRDdataArrayI.h" | |
40 | #include "AliTRDdigitsManager.h" | |
ca21baaa | 41 | #include "AliTRDrawData.h" |
42 | #include "AliTRDcalibDB.h" | |
43 | #include "AliTRDSimParam.h" | |
44 | #include "AliTRDRecParam.h" | |
45 | #include "AliTRDcluster.h" | |
af26ce80 | 46 | #include "AliTRDtransform.h" |
ca21baaa | 47 | #include "AliTRDSignalIndex.h" |
48 | #include "AliTRDRawStream.h" | |
ecf39416 | 49 | #include "AliTRDRawStreamV2.h" |
ecf39416 | 50 | #include "AliTRDfeeParam.h" |
ca21baaa | 51 | |
af26ce80 | 52 | #include "Cal/AliTRDCalROC.h" |
53 | #include "Cal/AliTRDCalDet.h" | |
54 | ||
ca21baaa | 55 | ClassImp(AliTRDclusterizerV2) |
56 | ||
57 | //_____________________________________________________________________________ | |
58 | AliTRDclusterizerV2::AliTRDclusterizerV2() | |
59 | :AliTRDclusterizer() | |
60 | ,fDigitsManager(NULL) | |
ca21baaa | 61 | ,fAddLabels(kTRUE) |
f582272a | 62 | ,fRawVersion(2) |
001be664 | 63 | ,fIndexesOut(NULL) |
64 | ,fIndexesMaxima(NULL) | |
af26ce80 | 65 | ,fTransform(NULL) |
ca21baaa | 66 | { |
67 | // | |
68 | // AliTRDclusterizerV2 default constructor | |
69 | // | |
70 | ||
ecf39416 | 71 | fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion(); |
af26ce80 | 72 | |
ca21baaa | 73 | } |
74 | ||
75 | //_____________________________________________________________________________ | |
76 | AliTRDclusterizerV2::AliTRDclusterizerV2(const Text_t *name, const Text_t *title) | |
77 | :AliTRDclusterizer(name,title) | |
78 | ,fDigitsManager(new AliTRDdigitsManager()) | |
ca21baaa | 79 | ,fAddLabels(kTRUE) |
f582272a | 80 | ,fRawVersion(2) |
001be664 | 81 | ,fIndexesOut(NULL) |
82 | ,fIndexesMaxima(NULL) | |
af26ce80 | 83 | ,fTransform(new AliTRDtransform(0)) |
ca21baaa | 84 | { |
85 | // | |
86 | // AliTRDclusterizerV2 constructor | |
87 | // | |
88 | ||
89 | fDigitsManager->CreateArrays(); | |
ecf39416 | 90 | |
91 | fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion(); | |
af26ce80 | 92 | |
ca21baaa | 93 | } |
94 | ||
95 | //_____________________________________________________________________________ | |
96 | AliTRDclusterizerV2::AliTRDclusterizerV2(const AliTRDclusterizerV2 &c) | |
97 | :AliTRDclusterizer(c) | |
98 | ,fDigitsManager(NULL) | |
ca21baaa | 99 | ,fAddLabels(kTRUE) |
f582272a | 100 | ,fRawVersion(2) |
001be664 | 101 | ,fIndexesOut(NULL) |
102 | ,fIndexesMaxima(NULL) | |
af26ce80 | 103 | ,fTransform(NULL) |
ca21baaa | 104 | { |
105 | // | |
106 | // AliTRDclusterizerV2 copy constructor | |
107 | // | |
108 | ||
109 | } | |
110 | ||
111 | //_____________________________________________________________________________ | |
112 | AliTRDclusterizerV2::~AliTRDclusterizerV2() | |
113 | { | |
114 | // | |
115 | // AliTRDclusterizerV2 destructor | |
116 | // | |
117 | ||
af26ce80 | 118 | if (fDigitsManager) |
ca21baaa | 119 | { |
af26ce80 | 120 | delete fDigitsManager; |
121 | fDigitsManager = NULL; | |
ca21baaa | 122 | } |
001be664 | 123 | |
124 | if (fIndexesOut) | |
125 | { | |
126 | delete fIndexesOut; | |
af26ce80 | 127 | fIndexesOut = NULL; |
001be664 | 128 | } |
129 | ||
130 | if (fIndexesMaxima) | |
131 | { | |
132 | delete fIndexesMaxima; | |
133 | fIndexesMaxima = NULL; | |
134 | } | |
af26ce80 | 135 | |
136 | if (fTransform) | |
137 | { | |
138 | delete fTransform; | |
139 | fTransform = NULL; | |
140 | } | |
141 | ||
ca21baaa | 142 | } |
143 | ||
144 | //_____________________________________________________________________________ | |
145 | AliTRDclusterizerV2 &AliTRDclusterizerV2::operator=(const AliTRDclusterizerV2 &c) | |
146 | { | |
147 | // | |
148 | // Assignment operator | |
149 | // | |
150 | ||
151 | if (this != &c) ((AliTRDclusterizerV2 &) c).Copy(*this); | |
152 | return *this; | |
153 | ||
154 | } | |
155 | ||
156 | //_____________________________________________________________________________ | |
157 | void AliTRDclusterizerV2::Copy(TObject &c) const | |
158 | { | |
159 | // | |
160 | // Copy function | |
161 | // | |
162 | ||
001be664 | 163 | ((AliTRDclusterizerV2 &) c).fDigitsManager = NULL; |
af26ce80 | 164 | ((AliTRDclusterizerV2 &) c).fAddLabels = fAddLabels; |
165 | ((AliTRDclusterizerV2 &) c).fRawVersion = fRawVersion; | |
166 | ((AliTRDclusterizerV2 &) c).fIndexesOut = NULL; | |
001be664 | 167 | ((AliTRDclusterizerV2 &) c).fIndexesMaxima = NULL; |
af26ce80 | 168 | ((AliTRDclusterizerV2 &) c).fTransform = NULL; |
001be664 | 169 | |
ca21baaa | 170 | AliTRDclusterizer::Copy(c); |
171 | ||
172 | } | |
173 | ||
001be664 | 174 | //_____________________________________________________________________________ |
175 | void AliTRDclusterizerV2::ResetHelperIndexes(AliTRDSignalIndex *indexesIn) | |
176 | { | |
177 | // | |
178 | // Reset the helper indexes | |
179 | // | |
af26ce80 | 180 | |
001be664 | 181 | if (fIndexesOut) |
182 | { | |
183 | // carefull here - we assume that only row number may change - most probable | |
184 | if (indexesIn->GetNrow() <= fIndexesOut->GetNrow()) | |
185 | fIndexesOut->ResetContent(); | |
186 | else | |
af26ce80 | 187 | fIndexesOut->ResetContentConditional(indexesIn->GetNrow() |
188 | , indexesIn->GetNcol() | |
189 | , indexesIn->GetNtime()); | |
001be664 | 190 | } |
191 | else | |
192 | { | |
af26ce80 | 193 | fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow() |
194 | , indexesIn->GetNcol() | |
195 | , indexesIn->GetNtime()); | |
001be664 | 196 | } |
197 | ||
198 | if (fIndexesMaxima) | |
199 | { | |
200 | // carefull here - we assume that only row number may change - most probable | |
af26ce80 | 201 | if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow()) |
202 | { | |
203 | fIndexesMaxima->ResetContent(); | |
204 | } | |
001be664 | 205 | else |
af26ce80 | 206 | { |
207 | fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow() | |
208 | , indexesIn->GetNcol() | |
209 | , indexesIn->GetNtime()); | |
210 | } | |
001be664 | 211 | } |
212 | else | |
213 | { | |
af26ce80 | 214 | fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow() |
215 | , indexesIn->GetNcol() | |
216 | , indexesIn->GetNtime()); | |
001be664 | 217 | } |
af26ce80 | 218 | |
001be664 | 219 | } |
220 | ||
ca21baaa | 221 | //_____________________________________________________________________________ |
222 | Bool_t AliTRDclusterizerV2::ReadDigits() | |
223 | { | |
224 | // | |
225 | // Reads the digits arrays from the input aliroot file | |
226 | // | |
227 | ||
228 | if (!fRunLoader) { | |
229 | AliError("No run loader available"); | |
230 | return kFALSE; | |
231 | } | |
232 | ||
233 | AliLoader* loader = fRunLoader->GetLoader("TRDLoader"); | |
234 | if (!loader->TreeD()) { | |
235 | loader->LoadDigits(); | |
236 | } | |
237 | ||
238 | // Read in the digit arrays | |
239 | return (fDigitsManager->ReadDigits(loader->TreeD())); | |
240 | ||
241 | } | |
242 | ||
243 | //_____________________________________________________________________________ | |
244 | Bool_t AliTRDclusterizerV2::ReadDigits(TTree *digitsTree) | |
245 | { | |
246 | // | |
247 | // Reads the digits arrays from the input tree | |
248 | // | |
249 | ||
250 | // Read in the digit arrays | |
251 | return (fDigitsManager->ReadDigits(digitsTree)); | |
252 | ||
253 | } | |
254 | ||
255 | //_____________________________________________________________________________ | |
256 | Bool_t AliTRDclusterizerV2::ReadDigits(AliRawReader *rawReader) | |
257 | { | |
258 | // | |
259 | // Reads the digits arrays from the ddl file | |
260 | // | |
261 | ||
262 | AliTRDrawData raw; | |
263 | fDigitsManager = raw.Raw2Digits(rawReader); | |
264 | ||
265 | return kTRUE; | |
266 | ||
267 | } | |
268 | ||
269 | //_____________________________________________________________________________ | |
270 | Bool_t AliTRDclusterizerV2::MakeClusters() | |
271 | { | |
272 | // | |
273 | // Creates clusters from digits | |
274 | // | |
275 | ||
af26ce80 | 276 | // Propagate info from the digits manager |
001be664 | 277 | if (fAddLabels == kTRUE) |
af26ce80 | 278 | { |
279 | fAddLabels = fDigitsManager->UsesDictionaries(); | |
280 | } | |
281 | ||
ca21baaa | 282 | Bool_t fReturn = kTRUE; |
283 | for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++) | |
284 | { | |
af26ce80 | 285 | |
ca21baaa | 286 | AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(i); |
287 | // This is to take care of switched off super modules | |
af26ce80 | 288 | if (digitsIn->GetNtime() == 0) |
289 | { | |
290 | continue; | |
291 | } | |
ca21baaa | 292 | digitsIn->Expand(); |
293 | AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i); | |
294 | if (indexes->IsAllocated() == kFALSE) | |
295 | { | |
296 | fDigitsManager->BuildIndexes(i); | |
297 | } | |
298 | ||
ecf39416 | 299 | Bool_t fR = kFALSE; |
ca21baaa | 300 | if (indexes->HasEntry()) |
301 | { | |
302 | if (fAddLabels) | |
303 | { | |
304 | for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++) | |
305 | { | |
306 | AliTRDdataArrayI *tracksIn = 0; | |
307 | tracksIn = fDigitsManager->GetDictionary(i,iDict); | |
308 | tracksIn->Expand(); | |
309 | } | |
310 | } | |
ecf39416 | 311 | fR = MakeClusters(i); |
ca21baaa | 312 | fReturn = fR && fReturn; |
313 | } | |
ecf39416 | 314 | |
315 | if (fR == kFALSE) | |
316 | { | |
317 | WriteClusters(i); | |
318 | ResetRecPoints(); | |
319 | } | |
af26ce80 | 320 | |
321 | // No compress just remove | |
ca21baaa | 322 | fDigitsManager->RemoveDigits(i); |
323 | fDigitsManager->RemoveDictionaries(i); | |
324 | fDigitsManager->ClearIndexes(i); | |
af26ce80 | 325 | |
ca21baaa | 326 | } |
001be664 | 327 | |
ca21baaa | 328 | return fReturn; |
ca21baaa | 329 | |
af26ce80 | 330 | } |
ca21baaa | 331 | |
332 | //_____________________________________________________________________________ | |
333 | Bool_t AliTRDclusterizerV2::Raw2Clusters(AliRawReader *rawReader) | |
334 | { | |
335 | // | |
336 | // Creates clusters from raw data | |
337 | // | |
338 | ||
ca21baaa | 339 | AliTRDdataArrayI *digits = 0; |
340 | AliTRDdataArrayI *track0 = 0; | |
341 | AliTRDdataArrayI *track1 = 0; | |
342 | AliTRDdataArrayI *track2 = 0; | |
343 | ||
344 | AliTRDSignalIndex *indexes = 0; | |
af26ce80 | 345 | |
ca21baaa | 346 | // Create the digits manager |
347 | if (!fDigitsManager) | |
348 | { | |
349 | fDigitsManager = new AliTRDdigitsManager(); | |
350 | fDigitsManager->CreateArrays(); | |
351 | } | |
ca21baaa | 352 | |
ecf39416 | 353 | AliTRDRawStreamV2 input(rawReader); |
f582272a | 354 | input.SetRawVersion( fRawVersion ); |
ca21baaa | 355 | input.Init(); |
356 | ||
ecf39416 | 357 | AliInfo(Form("Stream version: %s", input.IsA()->GetName())); |
358 | ||
ca21baaa | 359 | // Loop through the digits |
360 | Int_t lastdet = -1; | |
af26ce80 | 361 | Int_t det = 0; |
362 | Int_t it = 0; | |
ca21baaa | 363 | while (input.Next()) |
364 | { | |
365 | ||
af26ce80 | 366 | det = input.GetDet(); |
ca21baaa | 367 | |
368 | if (det != lastdet) | |
af26ce80 | 369 | { |
370 | ||
ca21baaa | 371 | if (lastdet != -1) |
372 | { | |
373 | digits = fDigitsManager->GetDigits(lastdet); | |
ecf39416 | 374 | Bool_t iclusterBranch = kFALSE; |
ca21baaa | 375 | if (indexes->HasEntry()) |
ecf39416 | 376 | iclusterBranch = MakeClusters(lastdet); |
377 | if (iclusterBranch == kFALSE) | |
378 | { | |
379 | WriteClusters(lastdet); | |
380 | ResetRecPoints(); | |
381 | } | |
ca21baaa | 382 | } |
383 | ||
384 | if (digits) | |
385 | { | |
386 | fDigitsManager->RemoveDigits(lastdet); | |
387 | fDigitsManager->RemoveDictionaries(lastdet); | |
388 | fDigitsManager->ClearIndexes(lastdet); | |
389 | } | |
390 | ||
391 | lastdet = det; | |
392 | ||
393 | // Add a container for the digits of this detector | |
394 | digits = fDigitsManager->GetDigits(det); | |
395 | track0 = fDigitsManager->GetDictionary(det,0); | |
396 | track1 = fDigitsManager->GetDictionary(det,1); | |
397 | track2 = fDigitsManager->GetDictionary(det,2); | |
398 | ||
399 | // Allocate memory space for the digits buffer | |
400 | if (digits->GetNtime() == 0) | |
401 | { | |
402 | //AliDebug(5, Form("Alloc digits for det %d", det)); | |
403 | digits->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); | |
404 | track0->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); | |
405 | track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); | |
406 | track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); | |
407 | } | |
408 | ||
409 | indexes = fDigitsManager->GetIndexes(det); | |
410 | indexes->SetSM(input.GetSM()); | |
411 | indexes->SetStack(input.GetStack()); | |
412 | indexes->SetLayer(input.GetLayer()); | |
413 | indexes->SetDetNumber(det); | |
414 | if (indexes->IsAllocated() == kFALSE) | |
af26ce80 | 415 | { |
416 | indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins()); | |
417 | } | |
418 | ||
ca21baaa | 419 | } |
420 | ||
421 | for (it = 0; it < 3; it++) | |
422 | { | |
423 | if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() ) | |
424 | { | |
425 | if (input.GetSignals()[it] > 0) | |
426 | { | |
427 | digits->SetDataUnchecked(input.GetRow(), input.GetCol(), | |
428 | input.GetTimeBin() + it, input.GetSignals()[it]); | |
429 | ||
430 | indexes->AddIndexTBin(input.GetRow(), input.GetCol(), | |
431 | input.GetTimeBin() + it); | |
432 | track0->SetDataUnchecked(input.GetRow(), input.GetCol(), | |
433 | input.GetTimeBin() + it, 0); | |
434 | track1->SetDataUnchecked(input.GetRow(), input.GetCol(), | |
435 | input.GetTimeBin() + it, 0); | |
436 | track2->SetDataUnchecked(input.GetRow(), input.GetCol(), | |
437 | input.GetTimeBin() + it, 0); | |
438 | } | |
439 | } | |
440 | } | |
af26ce80 | 441 | |
442 | } | |
ca21baaa | 443 | |
444 | if (lastdet != -1) | |
445 | { | |
ecf39416 | 446 | Bool_t iclusterBranch = kFALSE; |
af26ce80 | 447 | if (indexes->HasEntry()) |
448 | { | |
449 | iclusterBranch = MakeClusters(lastdet); | |
450 | } | |
ecf39416 | 451 | if (iclusterBranch == kFALSE) |
452 | { | |
453 | WriteClusters(lastdet); | |
454 | ResetRecPoints(); | |
455 | } | |
456 | //MakeClusters(lastdet); | |
ca21baaa | 457 | if (digits) |
458 | { | |
459 | fDigitsManager->RemoveDigits(lastdet); | |
460 | fDigitsManager->RemoveDictionaries(lastdet); | |
461 | fDigitsManager->ClearIndexes(lastdet); | |
462 | } | |
463 | } | |
464 | ||
001be664 | 465 | delete fDigitsManager; |
466 | fDigitsManager = NULL; | |
ca21baaa | 467 | return kTRUE; |
af26ce80 | 468 | |
ca21baaa | 469 | } |
470 | ||
471 | //_____________________________________________________________________________ | |
472 | Bool_t AliTRDclusterizerV2::Raw2ClustersChamber(AliRawReader *rawReader) | |
473 | { | |
474 | // | |
475 | // Creates clusters from raw data | |
476 | // | |
477 | ||
ca21baaa | 478 | // Create the digits manager |
479 | if (!fDigitsManager) | |
480 | { | |
481 | fDigitsManager = new AliTRDdigitsManager(); | |
482 | fDigitsManager->CreateArrays(); | |
483 | } | |
ca21baaa | 484 | |
001be664 | 485 | fDigitsManager->SetUseDictionaries(fAddLabels); |
ca21baaa | 486 | |
ecf39416 | 487 | AliTRDRawStreamV2 input(rawReader); |
f582272a | 488 | input.SetRawVersion( fRawVersion ); |
ca21baaa | 489 | input.Init(); |
ecf39416 | 490 | |
491 | AliInfo(Form("Stream version: %s", input.IsA()->GetName())); | |
ca21baaa | 492 | |
493 | Int_t det = 0; | |
494 | while ((det = input.NextChamber(fDigitsManager)) >= 0) | |
495 | { | |
ecf39416 | 496 | Bool_t iclusterBranch = kFALSE; |
ca21baaa | 497 | if (fDigitsManager->GetIndexes(det)->HasEntry()) |
af26ce80 | 498 | { |
499 | iclusterBranch = MakeClusters(det); | |
500 | } | |
ecf39416 | 501 | if (iclusterBranch == kFALSE) |
502 | { | |
503 | WriteClusters(det); | |
504 | ResetRecPoints(); | |
505 | } | |
ca21baaa | 506 | fDigitsManager->RemoveDigits(det); |
507 | fDigitsManager->RemoveDictionaries(det); | |
508 | fDigitsManager->ClearIndexes(det); | |
509 | } | |
510 | ||
001be664 | 511 | delete fDigitsManager; |
512 | fDigitsManager = NULL; | |
ca21baaa | 513 | return kTRUE; |
af26ce80 | 514 | |
ca21baaa | 515 | } |
516 | ||
517 | //_____________________________________________________________________________ | |
518 | Bool_t AliTRDclusterizerV2::MakeClusters(Int_t det) | |
519 | { | |
520 | // | |
521 | // Generates the cluster. | |
522 | // | |
523 | ||
524 | // Get the digits | |
525 | // digits should be expanded beforehand! | |
526 | // digitsIn->Expand(); | |
527 | AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(det); | |
af26ce80 | 528 | |
ca21baaa | 529 | // This is to take care of switched off super modules |
530 | if (digitsIn->GetNtime() == 0) | |
531 | { | |
ca21baaa | 532 | return kFALSE; |
533 | } | |
534 | ||
535 | AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det); | |
536 | if (indexesIn->IsAllocated() == kFALSE) | |
537 | { | |
538 | AliError("Indexes do not exist!"); | |
539 | return kFALSE; | |
540 | } | |
541 | ||
af26ce80 | 542 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); |
543 | if (!calibration) | |
544 | { | |
545 | AliFatal("No AliTRDcalibDB instance available\n"); | |
546 | return kFALSE; | |
547 | } | |
548 | ||
549 | AliTRDRecParam *recParam = AliTRDRecParam::Instance(); | |
550 | if (!recParam) | |
551 | { | |
552 | AliError("No AliTRDRecParam instance available\n"); | |
553 | return kFALSE; | |
554 | } | |
ca21baaa | 555 | |
556 | // ADC thresholds | |
af26ce80 | 557 | // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO |
558 | Float_t ADCthreshold = 0; | |
ecf39416 | 559 | |
ca21baaa | 560 | // Threshold value for the maximum |
561 | Float_t maxThresh = recParam->GetClusMaxThresh(); | |
562 | // Threshold value for the digit signal | |
563 | Float_t sigThresh = recParam->GetClusSigThresh(); | |
564 | ||
ca21baaa | 565 | // Iteration limit for unfolding procedure |
566 | const Float_t kEpsilon = 0.01; | |
567 | const Int_t kNclus = 3; | |
568 | const Int_t kNsig = 5; | |
569 | ||
570 | Int_t iUnfold = 0; | |
571 | Double_t ratioLeft = 1.0; | |
572 | Double_t ratioRight = 1.0; | |
573 | ||
574 | Double_t padSignal[kNsig]; | |
575 | Double_t clusterSignal[kNclus]; | |
ca21baaa | 576 | |
577 | Int_t icham = indexesIn->GetChamber(); | |
578 | Int_t iplan = indexesIn->GetPlane(); | |
579 | Int_t isect = indexesIn->GetSM(); | |
580 | ||
581 | // Start clustering in the chamber | |
582 | ||
af26ce80 | 583 | Int_t idet = AliTRDgeometry::GetDetector(iplan,icham,isect); |
ca21baaa | 584 | if (idet != det) |
585 | { | |
586 | AliError("Strange Detector number Missmatch!"); | |
587 | return kFALSE; | |
588 | } | |
589 | ||
af26ce80 | 590 | // TRD space point transformation |
591 | fTransform->SetDetector(det); | |
592 | ||
ca21baaa | 593 | Int_t ilayer = AliGeomManager::kTRD1 + iplan; |
594 | Int_t imodule = icham + AliTRDgeometry::Ncham() * isect; | |
595 | UShort_t volid = AliGeomManager::LayerToVolUID(ilayer,imodule); | |
596 | ||
af26ce80 | 597 | Int_t nColMax = digitsIn->GetNcol(); |
ca21baaa | 598 | Int_t nTimeTotal = digitsIn->GetNtime(); |
599 | ||
af26ce80 | 600 | // Detector wise calibration object for the gain factors |
601 | const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet(); | |
ca21baaa | 602 | // Calibration object with pad wise values for the gain factors |
af26ce80 | 603 | AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet); |
ca21baaa | 604 | // Calibration value for chamber wise gain factor |
af26ce80 | 605 | Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet); |
ca21baaa | 606 | |
af26ce80 | 607 | Int_t nClusters = 0; |
ca21baaa | 608 | |
ca21baaa | 609 | AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow() |
af26ce80 | 610 | ,digitsIn->GetNcol() |
611 | ,digitsIn->GetNtime()); | |
ca21baaa | 612 | |
001be664 | 613 | ResetHelperIndexes(indexesIn); |
ca21baaa | 614 | |
af26ce80 | 615 | // Apply the gain and the tail cancelation via digital filter |
616 | TailCancelation(digitsIn | |
617 | ,digitsOut | |
618 | ,indexesIn | |
619 | ,fIndexesOut | |
620 | ,nTimeTotal | |
621 | ,ADCthreshold | |
622 | ,calGainFactorROC | |
623 | ,calGainFactorDetValue); | |
ca21baaa | 624 | |
af26ce80 | 625 | Int_t row = 0; |
626 | Int_t col = 0; | |
627 | Int_t time = 0; | |
628 | Int_t iPad = 0; | |
ca21baaa | 629 | |
001be664 | 630 | fIndexesOut->ResetCounters(); |
631 | while (fIndexesOut->NextRCTbinIndex(row, col, time)) | |
ca21baaa | 632 | { |
af26ce80 | 633 | |
ca21baaa | 634 | Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time)); |
635 | ||
636 | // Look for the maximum | |
637 | if (signalM >= maxThresh) | |
638 | { | |
639 | ||
640 | if (col + 1 >= nColMax || col-1 < 0) | |
641 | continue; | |
642 | ||
af26ce80 | 643 | Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); |
ca21baaa | 644 | Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); |
645 | ||
646 | if ((TMath::Abs(signalL) <= signalM) && | |
647 | (TMath::Abs(signalR) < signalM)) | |
648 | { | |
649 | if ((TMath::Abs(signalL) >= sigThresh) || | |
650 | (TMath::Abs(signalR) >= sigThresh)) | |
651 | { | |
652 | // Maximum found, mark the position by a negative signal | |
653 | digitsOut->SetDataUnchecked(row,col,time,-signalM); | |
001be664 | 654 | fIndexesMaxima->AddIndexTBin(row,col,time); |
ca21baaa | 655 | } |
af26ce80 | 656 | } |
657 | ||
658 | } | |
659 | ||
ca21baaa | 660 | } |
661 | ||
662 | // The index to the first cluster of a given ROC | |
663 | Int_t firstClusterROC = -1; | |
664 | // The number of cluster in a given ROC | |
665 | Int_t nClusterROC = 0; | |
666 | ||
667 | // Now check the maxima and calculate the cluster position | |
001be664 | 668 | fIndexesMaxima->ResetCounters(); |
669 | while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) | |
ca21baaa | 670 | { |
af26ce80 | 671 | |
ca21baaa | 672 | // Maximum found ? |
af26ce80 | 673 | if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) |
674 | { | |
ca21baaa | 675 | |
af26ce80 | 676 | for (iPad = 0; iPad < kNclus; iPad++) |
677 | { | |
678 | Int_t iPadCol = col - 1 + iPad; | |
679 | clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time)); | |
680 | } | |
ca21baaa | 681 | |
af26ce80 | 682 | // Count the number of pads in the cluster |
683 | Int_t nPadCount = 0; | |
684 | Int_t ii; | |
685 | // Look to the left | |
686 | ii = 0; | |
687 | while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) | |
688 | { | |
689 | nPadCount++; | |
690 | ii++; | |
691 | if (col-ii < 0) break; | |
ca21baaa | 692 | } |
af26ce80 | 693 | // Look to the right |
694 | ii = 0; | |
695 | while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) | |
696 | { | |
697 | nPadCount++; | |
698 | ii++; | |
699 | if (col+ii+1 >= nColMax) break; | |
700 | } | |
701 | nClusters++; | |
702 | ||
703 | // Look for 5 pad cluster with minimum in the middle | |
704 | Bool_t fivePadCluster = kFALSE; | |
705 | if (col < (nColMax - 3)) | |
706 | { | |
707 | if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) | |
708 | { | |
709 | fivePadCluster = kTRUE; | |
710 | } | |
711 | if ((fivePadCluster) && (col < (nColMax - 5))) | |
712 | { | |
713 | if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) | |
714 | { | |
715 | fivePadCluster = kFALSE; | |
716 | } | |
717 | } | |
718 | if ((fivePadCluster) && (col > 1)) | |
719 | { | |
720 | if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) | |
721 | { | |
722 | fivePadCluster = kFALSE; | |
723 | } | |
724 | } | |
ca21baaa | 725 | } |
ca21baaa | 726 | |
af26ce80 | 727 | // 5 pad cluster |
728 | // Modify the signal of the overlapping pad for the left part | |
729 | // of the cluster which remains from a previous unfolding | |
730 | if (iUnfold) | |
731 | { | |
732 | clusterSignal[0] *= ratioLeft; | |
733 | iUnfold = 0; | |
734 | } | |
ca21baaa | 735 | |
af26ce80 | 736 | // Unfold the 5 pad cluster |
737 | if (fivePadCluster) | |
738 | { | |
739 | for (iPad = 0; iPad < kNsig; iPad++) | |
740 | { | |
741 | padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row | |
742 | ,col-1+iPad | |
743 | ,time)); | |
744 | } | |
745 | // Unfold the two maxima and set the signal on | |
746 | // the overlapping pad to the ratio | |
747 | ratioRight = Unfold(kEpsilon,iplan,padSignal); | |
748 | ratioLeft = 1.0 - ratioRight; | |
749 | clusterSignal[2] *= ratioRight; | |
750 | iUnfold = 1; | |
751 | } | |
ca21baaa | 752 | |
af26ce80 | 753 | // The position of the cluster in COL direction relative to the center pad (pad units) |
754 | Double_t clusterPosCol = 0.0; | |
755 | if (recParam->LUTOn()) | |
756 | { | |
757 | // Calculate the position of the cluster by using the | |
758 | // lookup table method | |
759 | clusterPosCol = recParam->LUTposition(iplan | |
760 | ,clusterSignal[0] | |
761 | ,clusterSignal[1] | |
762 | ,clusterSignal[2]); | |
763 | } | |
764 | else | |
765 | { | |
766 | // Calculate the position of the cluster by using the | |
767 | // center of gravity method | |
768 | for (Int_t i = 0; i < kNsig; i++) | |
769 | { | |
770 | padSignal[i] = 0.0; | |
771 | } | |
772 | padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad | |
773 | padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left pad | |
774 | padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right pad | |
775 | if ((col > 2) && | |
776 | (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) | |
777 | { | |
778 | padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)); | |
779 | } | |
780 | if ((col < nColMax - 3) && | |
781 | (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])) | |
782 | { | |
783 | padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)); | |
784 | } | |
785 | clusterPosCol = GetCOG(padSignal); | |
786 | } | |
ca21baaa | 787 | |
af26ce80 | 788 | // Store the amplitudes of the pads in the cluster for later analysis |
789 | Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 }; | |
790 | for (Int_t jPad = col-3; jPad <= col+3; jPad++) | |
791 | { | |
792 | if ((jPad < 0) || | |
793 | (jPad >= nColMax-1)) | |
794 | { | |
795 | continue; | |
796 | } | |
797 | signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time))); | |
798 | } | |
ca21baaa | 799 | |
af26ce80 | 800 | // Transform the local cluster coordinates into calibrated |
801 | // space point positions defined in the local tracking system. | |
802 | // Here the calibration for T0, Vdrift and ExB is applied as well. | |
803 | Double_t clusterXYZ[6]; | |
804 | clusterXYZ[0] = clusterPosCol; | |
805 | clusterXYZ[1] = clusterSignal[0]; | |
806 | clusterXYZ[2] = clusterSignal[1]; | |
807 | clusterXYZ[3] = clusterSignal[2]; | |
808 | clusterXYZ[4] = 0.0; | |
809 | clusterXYZ[5] = 0.0; | |
810 | Int_t clusterRCT[3]; | |
811 | clusterRCT[0] = row; | |
812 | clusterRCT[1] = col; | |
813 | clusterRCT[2] = 0; | |
814 | fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),0); | |
815 | ||
816 | // Add the cluster to the output array | |
817 | // The track indices will be stored later | |
818 | Float_t clusterPos[3]; | |
819 | clusterPos[0] = clusterXYZ[0]; | |
820 | clusterPos[1] = clusterXYZ[1]; | |
821 | clusterPos[2] = clusterXYZ[2]; | |
822 | Float_t clusterSig[2]; | |
823 | clusterSig[0] = clusterXYZ[4]; | |
824 | clusterSig[1] = clusterXYZ[5]; | |
825 | Double_t clusterCharge = clusterXYZ[3]; | |
826 | Char_t clusterTimeBin = ((Char_t) clusterRCT[2]); | |
827 | AliTRDcluster *cluster = new AliTRDcluster(idet | |
828 | ,clusterCharge | |
829 | ,clusterPos | |
830 | ,clusterSig | |
831 | ,0x0 | |
832 | ,((Char_t) nPadCount) | |
833 | ,signals | |
834 | ,((UChar_t) col) | |
835 | ,((UChar_t) row) | |
836 | ,((UChar_t) time) | |
837 | ,clusterTimeBin | |
838 | ,clusterPosCol | |
839 | ,volid); | |
840 | ||
841 | // Temporarily store the row, column and time bin of the center pad | |
842 | // Used to later on assign the track indices | |
843 | cluster->SetLabel( row,0); | |
844 | cluster->SetLabel( col,1); | |
845 | cluster->SetLabel(time,2); | |
846 | ||
847 | RecPoints()->Add(cluster); | |
848 | ||
849 | // Store the index of the first cluster in the current ROC | |
850 | if (firstClusterROC < 0) | |
851 | { | |
852 | firstClusterROC = RecPoints()->GetEntriesFast() - 1; | |
853 | } | |
ca21baaa | 854 | |
af26ce80 | 855 | // Count the number of cluster in the current ROC |
856 | nClusterROC++; | |
857 | ||
858 | } // if: Maximum found ? | |
ca21baaa | 859 | |
860 | } | |
af26ce80 | 861 | |
ca21baaa | 862 | delete digitsOut; |
ca21baaa | 863 | |
af26ce80 | 864 | if (fAddLabels) |
865 | { | |
866 | AddLabels(idet, firstClusterROC, nClusterROC); | |
867 | } | |
ca21baaa | 868 | |
869 | // Write the cluster and reset the array | |
870 | WriteClusters(idet); | |
871 | ResetRecPoints(); | |
872 | ||
873 | return kTRUE; | |
874 | ||
875 | } | |
876 | ||
877 | //_____________________________________________________________________________ | |
878 | Bool_t AliTRDclusterizerV2::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC) | |
879 | { | |
880 | // | |
881 | // Add the track indices to the found clusters | |
882 | // | |
883 | ||
af26ce80 | 884 | const Int_t kNclus = 3; |
885 | const Int_t kNdict = AliTRDdigitsManager::kNDict; | |
886 | const Int_t kNtrack = kNdict * kNclus; | |
ca21baaa | 887 | |
af26ce80 | 888 | Int_t iClusterROC = 0; |
ca21baaa | 889 | |
af26ce80 | 890 | Int_t row = 0; |
891 | Int_t col = 0; | |
892 | Int_t time = 0; | |
893 | Int_t iPad = 0; | |
ca21baaa | 894 | |
895 | // Temporary array to collect the track indices | |
896 | Int_t *idxTracks = new Int_t[kNtrack*nClusterROC]; | |
897 | ||
898 | // Loop through the dictionary arrays one-by-one | |
899 | // to keep memory consumption low | |
900 | AliTRDdataArrayI *tracksIn = 0; | |
901 | for (Int_t iDict = 0; iDict < kNdict; iDict++) { | |
902 | ||
ca21baaa | 903 | // tracksIn should be expanded beforehand! |
af26ce80 | 904 | tracksIn = fDigitsManager->GetDictionary(idet,iDict); |
ca21baaa | 905 | |
906 | // Loop though the clusters found in this ROC | |
907 | for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) { | |
908 | ||
909 | AliTRDcluster *cluster = (AliTRDcluster *) | |
910 | RecPoints()->UncheckedAt(firstClusterROC+iClusterROC); | |
911 | row = cluster->GetLabel(0); | |
912 | col = cluster->GetLabel(1); | |
913 | time = cluster->GetLabel(2); | |
914 | ||
915 | for (iPad = 0; iPad < kNclus; iPad++) { | |
916 | Int_t iPadCol = col - 1 + iPad; | |
917 | Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1; | |
918 | idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index; | |
919 | } | |
920 | ||
921 | } | |
922 | ||
ca21baaa | 923 | } |
924 | ||
925 | // Copy the track indices into the cluster | |
926 | // Loop though the clusters found in this ROC | |
927 | for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) { | |
928 | ||
929 | AliTRDcluster *cluster = (AliTRDcluster *) | |
930 | RecPoints()->UncheckedAt(firstClusterROC+iClusterROC); | |
931 | cluster->SetLabel(-9999,0); | |
932 | cluster->SetLabel(-9999,1); | |
933 | cluster->SetLabel(-9999,2); | |
934 | ||
935 | cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]); | |
936 | ||
937 | } | |
938 | ||
939 | delete [] idxTracks; | |
940 | ||
941 | return kTRUE; | |
af26ce80 | 942 | |
ca21baaa | 943 | } |
944 | ||
945 | //_____________________________________________________________________________ | |
946 | Double_t AliTRDclusterizerV2::GetCOG(Double_t signal[5]) | |
947 | { | |
948 | // | |
949 | // Get COG position | |
950 | // Used for clusters with more than 3 pads - where LUT not applicable | |
951 | // | |
952 | ||
953 | Double_t sum = signal[0] | |
954 | + signal[1] | |
955 | + signal[2] | |
956 | + signal[3] | |
957 | + signal[4]; | |
958 | ||
959 | Double_t res = (0.0 * (-signal[0] + signal[4]) | |
960 | + (-signal[1] + signal[3])) / sum; | |
961 | ||
962 | return res; | |
963 | ||
964 | } | |
965 | ||
966 | //_____________________________________________________________________________ | |
967 | Double_t AliTRDclusterizerV2::Unfold(Double_t eps, Int_t plane, Double_t *padSignal) | |
968 | { | |
969 | // | |
970 | // Method to unfold neighbouring maxima. | |
971 | // The charge ratio on the overlapping pad is calculated | |
972 | // until there is no more change within the range given by eps. | |
973 | // The resulting ratio is then returned to the calling method. | |
974 | // | |
975 | ||
976 | AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); | |
977 | if (!calibration) { | |
978 | AliError("No AliTRDcalibDB instance available\n"); | |
979 | return kFALSE; | |
980 | } | |
981 | ||
982 | Int_t irc = 0; | |
983 | Int_t itStep = 0; // Count iteration steps | |
984 | ||
985 | Double_t ratio = 0.5; // Start value for ratio | |
986 | Double_t prevRatio = 0.0; // Store previous ratio | |
987 | ||
988 | Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal | |
989 | Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal | |
990 | Double_t newSignal[3] = { 0.0, 0.0, 0.0 }; | |
991 | ||
992 | // Start the iteration | |
993 | while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { | |
994 | ||
995 | itStep++; | |
996 | prevRatio = ratio; | |
997 | ||
998 | // Cluster position according to charge ratio | |
999 | Double_t maxLeft = (ratio*padSignal[2] - padSignal[0]) | |
1000 | / (padSignal[0] + padSignal[1] + ratio*padSignal[2]); | |
1001 | Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) | |
1002 | / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]); | |
1003 | ||
1004 | // Set cluster charge ratio | |
1005 | irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal); | |
1006 | Double_t ampLeft = padSignal[1] / newSignal[1]; | |
1007 | irc = calibration->PadResponse(1.0,maxRight,plane,newSignal); | |
1008 | Double_t ampRight = padSignal[3] / newSignal[1]; | |
1009 | ||
1010 | // Apply pad response to parameters | |
1011 | irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal ); | |
1012 | irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal); | |
1013 | ||
1014 | // Calculate new overlapping ratio | |
af26ce80 | 1015 | ratio = TMath::Min((Double_t) 1.0 |
1016 | ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0])); | |
ca21baaa | 1017 | |
1018 | } | |
1019 | ||
1020 | return ratio; | |
1021 | ||
1022 | } | |
1023 | ||
1024 | //_____________________________________________________________________________ | |
af26ce80 | 1025 | void AliTRDclusterizerV2::TailCancelation(AliTRDdataArrayI *digitsIn |
1026 | , AliTRDdataArrayF *digitsOut | |
1027 | , AliTRDSignalIndex *indexesIn | |
1028 | , AliTRDSignalIndex *indexesOut | |
1029 | , Int_t nTimeTotal | |
1030 | , Float_t ADCthreshold | |
1031 | , AliTRDCalROC *calGainFactorROC | |
1032 | , Float_t calGainFactorDetValue) | |
ca21baaa | 1033 | { |
1034 | // | |
af26ce80 | 1035 | // Applies the tail cancelation and gain factors: |
1036 | // Transform digitsIn to digitsOut | |
ca21baaa | 1037 | // |
1038 | ||
1039 | Int_t iRow = 0; | |
1040 | Int_t iCol = 0; | |
1041 | Int_t iTime = 0; | |
1042 | ||
1043 | AliTRDRecParam *recParam = AliTRDRecParam::Instance(); | |
af26ce80 | 1044 | if (!recParam) |
1045 | { | |
1046 | AliError("No AliTRDRecParam instance available\n"); | |
1047 | return; | |
1048 | } | |
ca21baaa | 1049 | |
1050 | Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation | |
1051 | Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation | |
1052 | indexesIn->ResetCounters(); | |
1053 | while (indexesIn->NextRCIndex(iRow, iCol)) | |
1054 | { | |
1055 | Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow); | |
1056 | Double_t gain = calGainFactorDetValue | |
1057 | * calGainFactorROCValue; | |
1058 | ||
1059 | for (iTime = 0; iTime < nTimeTotal; iTime++) | |
001be664 | 1060 | { |
af26ce80 | 1061 | // Apply gain gain factor |
ca21baaa | 1062 | inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime); |
1063 | inADC[iTime] /= gain; | |
1064 | outADC[iTime] = inADC[iTime]; | |
1065 | } | |
1066 | ||
1067 | // Apply the tail cancelation via the digital filter | |
af26ce80 | 1068 | if (recParam->TCOn()) |
1069 | { | |
1070 | DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp()); | |
1071 | } | |
ecf39416 | 1072 | |
ca21baaa | 1073 | indexesIn->ResetTbinCounter(); |
1074 | while (indexesIn->NextTbinIndex(iTime)) | |
1075 | { | |
1076 | // Store the amplitude of the digit if above threshold | |
1077 | if (outADC[iTime] > ADCthreshold) | |
1078 | { | |
1079 | digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]); | |
ca21baaa | 1080 | indexesOut->AddIndexTBin(iRow,iCol,iTime); |
1081 | } | |
1082 | } //while itime | |
af26ce80 | 1083 | |
ca21baaa | 1084 | }//while irow icol |
1085 | ||
1086 | delete [] inADC; | |
1087 | delete [] outADC; | |
1088 | ||
ca21baaa | 1089 | return; |
1090 | ||
1091 | } | |
1092 | ||
1093 | //_____________________________________________________________________________ | |
1094 | void AliTRDclusterizerV2::DeConvExp(Double_t *source, Double_t *target | |
1095 | , Int_t n, Int_t nexp) | |
1096 | { | |
1097 | // | |
1098 | // Tail cancellation by deconvolution for PASA v4 TRF | |
1099 | // | |
1100 | ||
1101 | Double_t rates[2]; | |
1102 | Double_t coefficients[2]; | |
1103 | ||
1104 | // Initialization (coefficient = alpha, rates = lambda) | |
af26ce80 | 1105 | Double_t r1 = 1.0; |
1106 | Double_t r2 = 1.0; | |
1107 | Double_t c1 = 0.5; | |
1108 | Double_t c2 = 0.5; | |
ca21baaa | 1109 | |
1110 | if (nexp == 1) { // 1 Exponentials | |
af26ce80 | 1111 | r1 = 1.156; |
1112 | r2 = 0.130; | |
1113 | c1 = 0.066; | |
1114 | c2 = 0.000; | |
ca21baaa | 1115 | } |
1116 | if (nexp == 2) { // 2 Exponentials | |
af26ce80 | 1117 | r1 = 1.156; |
1118 | r2 = 0.130; | |
1119 | c1 = 0.114; | |
1120 | c2 = 0.624; | |
ca21baaa | 1121 | } |
1122 | ||
af26ce80 | 1123 | coefficients[0] = c1; |
1124 | coefficients[1] = c2; | |
ca21baaa | 1125 | |
af26ce80 | 1126 | Double_t dt = 0.1; |
ca21baaa | 1127 | |
af26ce80 | 1128 | rates[0] = TMath::Exp(-dt/(r1)); |
1129 | rates[1] = TMath::Exp(-dt/(r2)); | |
ca21baaa | 1130 | |
1131 | Int_t i = 0; | |
1132 | Int_t k = 0; | |
1133 | ||
1134 | Double_t reminder[2]; | |
af26ce80 | 1135 | Double_t correction = 0.0; |
1136 | Double_t result = 0.0; | |
ca21baaa | 1137 | |
1138 | // Attention: computation order is important | |
ca21baaa | 1139 | for (k = 0; k < nexp; k++) { |
1140 | reminder[k] = 0.0; | |
1141 | } | |
af26ce80 | 1142 | |
ca21baaa | 1143 | for (i = 0; i < n; i++) { |
af26ce80 | 1144 | |
ca21baaa | 1145 | result = (source[i] - correction); // No rescaling |
1146 | target[i] = result; | |
1147 | ||
1148 | for (k = 0; k < nexp; k++) { | |
1149 | reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result); | |
1150 | } | |
af26ce80 | 1151 | |
ca21baaa | 1152 | correction = 0.0; |
1153 | for (k = 0; k < nexp; k++) { | |
1154 | correction += reminder[k]; | |
1155 | } | |
af26ce80 | 1156 | |
ca21baaa | 1157 | } |
1158 | ||
1159 | } |