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