]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterizer.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
QA update by Sylwester
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizer.cxx
CommitLineData
3fe61b77 1
f7336fa3 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
88cb7938 17/* $Id$ */
f7336fa3 18
19///////////////////////////////////////////////////////////////////////////////
20// //
3fe61b77 21// TRD cluster finder //
f7336fa3 22// //
23///////////////////////////////////////////////////////////////////////////////
24
3fe61b77 25#include <TF1.h>
94de3818 26#include <TTree.h>
3fe61b77 27#include <TH1.h>
793ff80c 28#include <TFile.h>
6d50f529 29#include <TObjArray.h>
f7336fa3 30
88cb7938 31#include "AliRunLoader.h"
32#include "AliLoader.h"
3fe61b77 33#include "AliRawReader.h"
6d50f529 34#include "AliLog.h"
3fe61b77 35#include "AliAlignObj.h"
88cb7938 36
f7336fa3 37#include "AliTRDclusterizer.h"
3e1a3ad8 38#include "AliTRDcluster.h"
fc546d21 39#include "AliTRDReconstructor.h"
793ff80c 40#include "AliTRDgeometry.h"
3fe61b77 41#include "AliTRDdataArrayF.h"
42#include "AliTRDdataArrayI.h"
625f5260 43#include "AliTRDdataArrayS.h"
fb44eb8e 44#include "AliTRDdataArrayDigits.h"
3fe61b77 45#include "AliTRDdigitsManager.h"
46#include "AliTRDrawData.h"
3551db50 47#include "AliTRDcalibDB.h"
fc546d21 48#include "AliTRDrecoParam.h"
b43a3e17 49#include "AliTRDCommonParam.h"
3fe61b77 50#include "AliTRDtransform.h"
51#include "AliTRDSignalIndex.h"
dfbb4bb9 52#include "AliTRDrawStreamBase.h"
3fe61b77 53#include "AliTRDfeeParam.h"
54
55#include "Cal/AliTRDCalROC.h"
56#include "Cal/AliTRDCalDet.h"
0e09df31 57#include "Cal/AliTRDCalSingleChamberStatus.h"
f7336fa3 58
59ClassImp(AliTRDclusterizer)
60
61//_____________________________________________________________________________
6d50f529 62AliTRDclusterizer::AliTRDclusterizer()
63 :TNamed()
64 ,fRunLoader(NULL)
65 ,fClusterTree(NULL)
66 ,fRecPoints(NULL)
3fe61b77 67 ,fDigitsManager(NULL)
68 ,fAddLabels(kTRUE)
69 ,fRawVersion(2)
70 ,fIndexesOut(NULL)
71 ,fIndexesMaxima(NULL)
dfbb4bb9 72 ,fTransform(new AliTRDtransform(0))
fc546d21 73 ,fLUTbin(0)
74 ,fLUT(NULL)
f7336fa3 75{
76 //
77 // AliTRDclusterizer default constructor
78 //
79
a7ac01d2 80 AliTRDcalibDB *trd = 0x0;
81 if (!(trd = AliTRDcalibDB::Instance())) {
82 AliFatal("Could not get calibration object");
83 }
84
3fe61b77 85 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
86
a7ac01d2 87 // retrive reco params
88 AliTRDrecoParam *rec = 0x0;
89 if (!(rec = AliTRDReconstructor::RecoParam())){
90 if(!(rec = trd->GetRecoParam(0))){
91 AliInfo("Using default RecoParams = LowFluxParam.");
92 rec = AliTRDrecoParam::GetLowFluxParam();
93 }
94 AliTRDReconstructor::SetRecoParam(rec);
95 }
f7336fa3 96}
97
98//_____________________________________________________________________________
3fe61b77 99AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title)
6d50f529 100 :TNamed(name,title)
101 ,fRunLoader(NULL)
102 ,fClusterTree(NULL)
103 ,fRecPoints(NULL)
3fe61b77 104 ,fDigitsManager(new AliTRDdigitsManager())
105 ,fAddLabels(kTRUE)
106 ,fRawVersion(2)
107 ,fIndexesOut(NULL)
108 ,fIndexesMaxima(NULL)
109 ,fTransform(new AliTRDtransform(0))
fc546d21 110 ,fLUTbin(0)
111 ,fLUT(NULL)
f7336fa3 112{
113 //
6d50f529 114 // AliTRDclusterizer constructor
f7336fa3 115 //
116
a7ac01d2 117 AliTRDcalibDB *trd = 0x0;
118 if (!(trd = AliTRDcalibDB::Instance())) {
119 AliFatal("Could not get calibration object");
120 }
121
122 // retrive reco params
123 AliTRDrecoParam *rec = 0x0;
124 if (!(rec = AliTRDReconstructor::RecoParam())){
125 if(!(rec = trd->GetRecoParam(0))){
126 AliInfo("Using default RecoParams = LowFluxParam.");
127 rec = AliTRDrecoParam::GetLowFluxParam();
128 }
129 AliTRDReconstructor::SetRecoParam(rec);
130 }
131
3fe61b77 132 fDigitsManager->CreateArrays();
133
134 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
135
fc546d21 136 FillLUT();
f7336fa3 137}
138
8230f242 139//_____________________________________________________________________________
6d50f529 140AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
141 :TNamed(c)
142 ,fRunLoader(NULL)
143 ,fClusterTree(NULL)
144 ,fRecPoints(NULL)
3fe61b77 145 ,fDigitsManager(NULL)
146 ,fAddLabels(kTRUE)
147 ,fRawVersion(2)
148 ,fIndexesOut(NULL)
149 ,fIndexesMaxima(NULL)
150 ,fTransform(NULL)
fc546d21 151 ,fLUTbin(0)
152 ,fLUT(0)
8230f242 153{
154 //
155 // AliTRDclusterizer copy constructor
156 //
157
fc546d21 158 FillLUT();
159
8230f242 160}
161
f7336fa3 162//_____________________________________________________________________________
163AliTRDclusterizer::~AliTRDclusterizer()
164{
8230f242 165 //
166 // AliTRDclusterizer destructor
167 //
f7336fa3 168
3fe61b77 169 if (fRecPoints)
170 {
171 fRecPoints->Delete();
172 delete fRecPoints;
173 }
174
175 if (fDigitsManager)
176 {
177 delete fDigitsManager;
178 fDigitsManager = NULL;
179 }
180
181 if (fIndexesOut)
182 {
183 delete fIndexesOut;
184 fIndexesOut = NULL;
185 }
186
187 if (fIndexesMaxima)
188 {
189 delete fIndexesMaxima;
190 fIndexesMaxima = NULL;
191 }
192
193 if (fTransform)
194 {
195 delete fTransform;
196 fTransform = NULL;
197 }
6d50f529 198
fc546d21 199 if (fLUT)
200 {
201 delete [] fLUT;
202 fLUT = NULL;
203 }
204
f7336fa3 205}
206
8230f242 207//_____________________________________________________________________________
dd9a6ee3 208AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
209{
210 //
211 // Assignment operator
212 //
213
3fe61b77 214 if (this != &c)
215 {
216 ((AliTRDclusterizer &) c).Copy(*this);
217 }
053767a4 218
dd9a6ee3 219 return *this;
220
221}
222
223//_____________________________________________________________________________
e0d47c25 224void AliTRDclusterizer::Copy(TObject &c) const
8230f242 225{
226 //
227 // Copy function
228 //
229
3fe61b77 230 ((AliTRDclusterizer &) c).fClusterTree = NULL;
231 ((AliTRDclusterizer &) c).fRecPoints = NULL;
232 ((AliTRDclusterizer &) c).fDigitsManager = NULL;
233 ((AliTRDclusterizer &) c).fAddLabels = fAddLabels;
234 ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
235 ((AliTRDclusterizer &) c).fIndexesOut = NULL;
236 ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
237 ((AliTRDclusterizer &) c).fTransform = NULL;
fc546d21 238 ((AliTRDclusterizer &) c).fLUTbin = 0;
239 ((AliTRDclusterizer &) c).fLUT = NULL;
8230f242 240
241}
242
f7336fa3 243//_____________________________________________________________________________
3e1a3ad8 244Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
f7336fa3 245{
246 //
3e1a3ad8 247 // Opens the AliROOT file. Output and input are in the same file
f7336fa3 248 //
6d50f529 249
e191bb57 250 TString evfoldname = AliConfig::GetDefaultEventFolderName();
6d50f529 251 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
252
253 if (!fRunLoader) {
19dd5b2f 254 fRunLoader = AliRunLoader::Open(name);
6d50f529 255 }
256
257 if (!fRunLoader) {
258 AliError(Form("Can not open session for file %s.",name));
259 return kFALSE;
260 }
88cb7938 261
262 OpenInput(nEvent);
263 OpenOutput();
6d50f529 264
3e1a3ad8 265 return kTRUE;
f7336fa3 266
6d50f529 267}
3e1a3ad8 268
269//_____________________________________________________________________________
88cb7938 270Bool_t AliTRDclusterizer::OpenOutput()
3e1a3ad8 271{
272 //
273 // Open the output file
274 //
275
3e1a3ad8 276 TObjArray *ioArray = 0;
88cb7938 277
278 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
279 loader->MakeTree("R");
6d50f529 280
88cb7938 281 fClusterTree = loader->TreeR();
365d0374 282 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
3e1a3ad8 283
3e1a3ad8 284 return kTRUE;
285
286}
287
25ca55ce 288//_____________________________________________________________________________
289Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
290{
291 //
292 // Connect the output tree
293 //
294
295 TObjArray *ioArray = 0;
296
297 fClusterTree = clusterTree;
298 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
299
300 return kTRUE;
301
302}
303
3e1a3ad8 304//_____________________________________________________________________________
88cb7938 305Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
f7336fa3 306{
307 //
308 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
309 //
310
f7336fa3 311 // Import the Trees for the event nEvent in the file
bdbb05bb 312 fRunLoader->GetEvent(nEvent);
88cb7938 313
f7336fa3 314 return kTRUE;
315
316}
317
318//_____________________________________________________________________________
793ff80c 319Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
f7336fa3 320{
321 //
3e1a3ad8 322 // Fills TRDcluster branch in the tree with the clusters
793ff80c 323 // found in detector = det. For det=-1 writes the tree.
a3c76cdc 324 //
793ff80c 325
6d50f529 326 if ((det < -1) ||
327 (det >= AliTRDgeometry::Ndet())) {
328 AliError(Form("Unexpected detector index %d.\n",det));
3e1a3ad8 329 return kFALSE;
793ff80c 330 }
3e1a3ad8 331
3e1a3ad8 332 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
333 if (!branch) {
793ff80c 334 TObjArray *ioArray = 0;
3e1a3ad8 335 branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
793ff80c 336 }
337
6d50f529 338 if ((det >= 0) &&
339 (det < AliTRDgeometry::Ndet())) {
793ff80c 340
bdbb05bb 341 Int_t nRecPoints = RecPoints()->GetEntriesFast();
3e1a3ad8 342 TObjArray *detRecPoints = new TObjArray(400);
343
344 for (Int_t i = 0; i < nRecPoints; i++) {
bdbb05bb 345 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
3e1a3ad8 346 if (det == c->GetDetector()) {
347 detRecPoints->AddLast(c);
348 }
349 else {
3fe61b77 350 AliError(Form("Attempt to write a cluster with unexpected detector index: got=%d expected=%d\n"
351 ,c->GetDetector()
352 ,det));
3e1a3ad8 353 }
793ff80c 354 }
355
3e1a3ad8 356 branch->SetAddress(&detRecPoints);
357 fClusterTree->Fill();
358
d9b8978b 359 delete detRecPoints;
ecf39416 360
793ff80c 361 return kTRUE;
3e1a3ad8 362
793ff80c 363 }
364
365 if (det == -1) {
366
6d50f529 367 AliInfo(Form("Writing the cluster tree %s for event %d."
368 ,fClusterTree->GetName(),fRunLoader->GetEventNumber()));
369
a6dd11e9 370 if (fRecPoints) {
371
372 branch->SetAddress(&fRecPoints);
373
374 AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
375 loader->WriteRecPoints("OVERWRITE");
c630aafd 376
a6dd11e9 377 }
378 else {
379
380 AliError("Cluster tree does not exist. Cannot write clusters.\n");
381 return kFALSE;
382
383 }
384
793ff80c 385 return kTRUE;
3e1a3ad8 386
793ff80c 387 }
6d50f529 388
389 AliError(Form("Unexpected detector index %d.\n",det));
3e1a3ad8 390
793ff80c 391 return kFALSE;
88cb7938 392
f7336fa3 393}
793ff80c 394
3551db50 395//_____________________________________________________________________________
3fe61b77 396void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
397{
398 //
399 // Reset the helper indexes
400 //
401
402 if (fIndexesOut)
403 {
404 // carefull here - we assume that only row number may change - most probable
405 if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
406 fIndexesOut->ResetContent();
407 else
408 fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
409 , indexesIn->GetNcol()
410 , indexesIn->GetNtime());
411 }
412 else
413 {
414 fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
415 , indexesIn->GetNcol()
416 , indexesIn->GetNtime());
417 }
418
419 if (fIndexesMaxima)
420 {
421 // carefull here - we assume that only row number may change - most probable
422 if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
423 {
424 fIndexesMaxima->ResetContent();
425 }
426 else
427 {
428 fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
429 , indexesIn->GetNcol()
430 , indexesIn->GetNtime());
431 }
432 }
433 else
434 {
435 fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
436 , indexesIn->GetNcol()
437 , indexesIn->GetNtime());
438 }
439
440}
441
442//_____________________________________________________________________________
443Bool_t AliTRDclusterizer::ReadDigits()
444{
445 //
446 // Reads the digits arrays from the input aliroot file
447 //
448
449 if (!fRunLoader) {
450 AliError("No run loader available");
451 return kFALSE;
452 }
453
454 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
455 if (!loader->TreeD()) {
456 loader->LoadDigits();
457 }
458
459 // Read in the digit arrays
460 return (fDigitsManager->ReadDigits(loader->TreeD()));
461
462}
463
464//_____________________________________________________________________________
465Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
466{
467 //
468 // Reads the digits arrays from the input tree
469 //
470
471 // Read in the digit arrays
472 return (fDigitsManager->ReadDigits(digitsTree));
473
474}
475
476//_____________________________________________________________________________
477Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
478{
479 //
480 // Reads the digits arrays from the ddl file
481 //
482
483 AliTRDrawData raw;
484 fDigitsManager = raw.Raw2Digits(rawReader);
485
486 return kTRUE;
487
488}
489
490//_____________________________________________________________________________
491Bool_t AliTRDclusterizer::MakeClusters()
492{
493 //
494 // Creates clusters from digits
495 //
496
497 // Propagate info from the digits manager
498 if (fAddLabels == kTRUE)
499 {
500 fAddLabels = fDigitsManager->UsesDictionaries();
501 }
502
503 Bool_t fReturn = kTRUE;
504 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
505 {
506
fb44eb8e 507 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits*) fDigitsManager->GetDigits(i);
3fe61b77 508 // This is to take care of switched off super modules
625f5260 509 if (!digitsIn->HasData())
3fe61b77 510 {
511 continue;
512 }
513 digitsIn->Expand();
514 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
515 if (indexes->IsAllocated() == kFALSE)
516 {
517 fDigitsManager->BuildIndexes(i);
518 }
519
520 Bool_t fR = kFALSE;
521 if (indexes->HasEntry())
522 {
523 if (fAddLabels)
524 {
525 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
526 {
527 AliTRDdataArrayI *tracksIn = 0;
625f5260 528 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(i,iDict);
3fe61b77 529 tracksIn->Expand();
530 }
531 }
532 fR = MakeClusters(i);
533 fReturn = fR && fReturn;
534 }
535
536 if (fR == kFALSE)
537 {
538 WriteClusters(i);
539 ResetRecPoints();
540 }
541
542 // No compress just remove
543 fDigitsManager->RemoveDigits(i);
544 fDigitsManager->RemoveDictionaries(i);
545 fDigitsManager->ClearIndexes(i);
546
547 }
548
549 return fReturn;
550
551}
552
553//_____________________________________________________________________________
554Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
555{
556 //
557 // Creates clusters from raw data
558 //
559
fc546d21 560 return Raw2ClustersChamber(rawReader);
3fe61b77 561
3fe61b77 562}
563
564//_____________________________________________________________________________
565Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
566{
567 //
568 // Creates clusters from raw data
569 //
570
571 // Create the digits manager
572 if (!fDigitsManager)
573 {
574 fDigitsManager = new AliTRDdigitsManager();
575 fDigitsManager->CreateArrays();
576 }
577
578 fDigitsManager->SetUseDictionaries(fAddLabels);
579
dfbb4bb9 580 AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
581 AliTRDrawStreamBase &input = *pinput;
3fe61b77 582
583 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
584
585 Int_t det = 0;
586 while ((det = input.NextChamber(fDigitsManager)) >= 0)
587 {
588 Bool_t iclusterBranch = kFALSE;
589 if (fDigitsManager->GetIndexes(det)->HasEntry())
590 {
591 iclusterBranch = MakeClusters(det);
592 }
593 if (iclusterBranch == kFALSE)
594 {
595 WriteClusters(det);
596 ResetRecPoints();
597 }
598 fDigitsManager->RemoveDigits(det);
599 fDigitsManager->RemoveDictionaries(det);
600 fDigitsManager->ClearIndexes(det);
601 }
602
603 delete fDigitsManager;
604 fDigitsManager = NULL;
dfbb4bb9 605
606 delete pinput;
607 pinput = NULL;
3fe61b77 608 return kTRUE;
609
610}
611
fa7427d0 612//_____________________________________________________________________________
613UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
614{
615 UChar_t status = 0;
616 // check if pad is masked
617 if(signal>0 && TESTBIT(signal, 10)){
618 CLRBIT(signal, 10);
619 for(int ibit=0; ibit<4; ibit++){
620 if(TESTBIT(signal, 11+ibit)){
621 SETBIT(status, ibit);
622 CLRBIT(signal, 11+ibit);
623 }
624 }
625 }
626 return status;
627}
628
3fe61b77 629//_____________________________________________________________________________
630Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
631{
632 //
633 // Generates the cluster.
634 //
635
636 // Get the digits
637 // digits should be expanded beforehand!
638 // digitsIn->Expand();
fb44eb8e 639 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits *) fDigitsManager->GetDigits(det);
f5375dcb 640
3fe61b77 641 // This is to take care of switched off super modules
625f5260 642 if (!digitsIn->HasData())
3fe61b77 643 {
644 return kFALSE;
645 }
646
647 AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
648 if (indexesIn->IsAllocated() == kFALSE)
649 {
650 AliError("Indexes do not exist!");
651 return kFALSE;
652 }
653
654 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
655 if (!calibration)
656 {
657 AliFatal("No AliTRDcalibDB instance available\n");
658 return kFALSE;
659 }
660
3fe61b77 661 // ADC thresholds
662 // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
acc49af9 663 Float_t adcThreshold = 0;
3fe61b77 664
c5f589b9 665 if (!AliTRDReconstructor::RecoParam())
666 {
667 AliError("RecoParam does not exist\n");
668 return kFALSE;
669 }
670
3fe61b77 671 // Threshold value for the maximum
fc546d21 672 Float_t maxThresh = AliTRDReconstructor::RecoParam()->GetClusMaxThresh();
3fe61b77 673 // Threshold value for the digit signal
fc546d21 674 Float_t sigThresh = AliTRDReconstructor::RecoParam()->GetClusSigThresh();
3fe61b77 675
df83a620 676 // Threshold value for the maximum ( cut noise)
677 Float_t minMaxCutSigma = AliTRDReconstructor::RecoParam()->GetMinMaxCutSigma();
678 // Threshold value for the sum pad ( cut noise)
679 Float_t minLeftRightCutSigma = AliTRDReconstructor::RecoParam()->GetMinLeftRightCutSigma();
680
3fe61b77 681 // Iteration limit for unfolding procedure
682 const Float_t kEpsilon = 0.01;
683 const Int_t kNclus = 3;
684 const Int_t kNsig = 5;
685
686 Int_t iUnfold = 0;
687 Double_t ratioLeft = 1.0;
688 Double_t ratioRight = 1.0;
689
690 Double_t padSignal[kNsig];
691 Double_t clusterSignal[kNclus];
692
053767a4 693 Int_t istack = indexesIn->GetStack();
694 Int_t ilayer = indexesIn->GetLayer();
695 Int_t isector = indexesIn->GetSM();
3fe61b77 696
697 // Start clustering in the chamber
698
053767a4 699 Int_t idet = AliTRDgeometry::GetDetector(ilayer,istack,isector);
3fe61b77 700 if (idet != det)
701 {
702 AliError("Strange Detector number Missmatch!");
703 return kFALSE;
704 }
705
706 // TRD space point transformation
707 fTransform->SetDetector(det);
708
053767a4 709 Int_t iGeoLayer = AliGeomManager::kTRD1 + ilayer;
710 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
711 UShort_t volid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
3fe61b77 712
713 Int_t nColMax = digitsIn->GetNcol();
fa7427d0 714 Int_t nRowMax = digitsIn->GetNrow();
3fe61b77 715 Int_t nTimeTotal = digitsIn->GetNtime();
716
717 // Detector wise calibration object for the gain factors
0e09df31 718 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
3fe61b77 719 // Calibration object with pad wise values for the gain factors
0e09df31 720 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
3fe61b77 721 // Calibration value for chamber wise gain factor
0e09df31 722 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
723
df83a620 724
725 // Detector wise calibration object for the noise
726 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
727 // Calibration object with pad wise values for the noise
728 AliTRDCalROC *calNoiseROC = calibration->GetNoiseROC(idet);
729 // Calibration value for chamber wise noise
730 Float_t calNoiseDetValue = calNoiseDet->GetValue(idet);
731
3fe61b77 732 Int_t nClusters = 0;
733
fa7427d0 734 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(nRowMax, nColMax, nTimeTotal);
735 AliTRDdataArrayS padStatus(nRowMax, nColMax, nTimeTotal);
3fe61b77 736
737 ResetHelperIndexes(indexesIn);
738
739 // Apply the gain and the tail cancelation via digital filter
740 TailCancelation(digitsIn
f5375dcb 741 ,digitsOut
742 ,indexesIn
743 ,fIndexesOut
744 ,nTimeTotal
745 ,adcThreshold
746 ,calGainFactorROC
747 ,calGainFactorDetValue);
3fe61b77 748
749 Int_t row = 0;
750 Int_t col = 0;
751 Int_t time = 0;
752 Int_t iPad = 0;
753
f5375dcb 754 UChar_t status[3]={0, 0, 0}, ipos = 0;
3fe61b77 755 fIndexesOut->ResetCounters();
fa7427d0 756 while (fIndexesOut->NextRCTbinIndex(row, col, time)){
f5375dcb 757 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
758 status[1] = digitsIn->GetPadStatus(row,col,time);
759 if(status[1]) SETBIT(ipos, AliTRDcluster::kMaskedCenter);
760
df83a620 761 if(signalM < maxThresh) continue;
762
763 Float_t noiseMiddleThresh = minMaxCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
764 if (signalM < noiseMiddleThresh) continue;
765
766 if (col + 1 >= nColMax || col-1 < 0) continue;
f5375dcb 767
df83a620 768 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
769 status[0] = digitsIn->GetPadStatus(row,col+1,time);
770 if(status[0]) SETBIT(ipos, AliTRDcluster::kMaskedLeft);
f5375dcb 771
df83a620 772 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
773 status[2] = digitsIn->GetPadStatus(row,col-1,time);
774 if(status[2]) SETBIT(ipos, AliTRDcluster::kMaskedRight);
f5375dcb 775
df83a620 776 // reject candidates with more than 1 problematic pad
777 if(ipos == 3 || ipos > 4) continue;
f5375dcb 778
df83a620 779 if(!status[1]){ // good central pad
780 if(!ipos){ // all pads are OK
781 if ((signalL <= signalM) && (signalR < signalM)) {
782 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
783 Float_t noiseSumThresh = minLeftRightCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
784 if((signalL+signalR+signalM) >= noiseSumThresh){
785 // Maximum found, mark the position by a negative signal
786 digitsOut->SetDataUnchecked(row,col,time,-signalM);
787 fIndexesMaxima->AddIndexTBin(row,col,time);
788 padStatus.SetDataUnchecked(row, col, time, ipos);
789 }
790 }
791 }
792 } else { // one of the neighbouring pads are bad
793 if(status[0] && signalR < signalM && signalR >= sigThresh){
794 digitsOut->SetDataUnchecked(row,col,time,-signalM);
795 digitsOut->SetDataUnchecked(row, col, time+1, 0.);
796 fIndexesMaxima->AddIndexTBin(row,col,time);
797 padStatus.SetDataUnchecked(row, col, time, ipos);
798 } else if(status[2] && signalL <= signalM && signalL >= sigThresh){
799 digitsOut->SetDataUnchecked(row,col,time,-signalM);
800 digitsOut->SetDataUnchecked(row, col, time-1, 0.);
f5375dcb 801 fIndexesMaxima->AddIndexTBin(row,col,time);
802 padStatus.SetDataUnchecked(row, col, time, ipos);
df83a620 803 }
804 }
805 } else { // wrong maximum pad
806 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
807 // Maximum found, mark the position by a negative signal
808 digitsOut->SetDataUnchecked(row,col,time,-maxThresh);
809 fIndexesMaxima->AddIndexTBin(row,col,time);
810 padStatus.SetDataUnchecked(row, col, time, ipos);
f5375dcb 811 }
812 }
df83a620 813 }
f5375dcb 814
815 // The index to the first cluster of a given ROC
df83a620 816 Int_t firstClusterROC = -1;
f5375dcb 817 // The number of cluster in a given ROC
818 Int_t nClusterROC = 0;
819
820 // Now check the maxima and calculate the cluster position
821 fIndexesMaxima->ResetCounters();
822 while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) {
3fe61b77 823
824 // Maximum found ?
f5375dcb 825 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
826 for (iPad = 0; iPad < kNclus; iPad++) {
827 Int_t iPadCol = col - 1 + iPad;
828 clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
829 }
3fe61b77 830
f5375dcb 831 // Count the number of pads in the cluster
832 Int_t nPadCount = 0;
833 Int_t ii;
834 // Look to the left
835 ii = 0;
836 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) {
837 nPadCount++;
838 ii++;
839 if (col-ii < 0) break;
840 }
841 // Look to the right
842 ii = 0;
843 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh){
844 nPadCount++;
845 ii++;
846 if (col+ii+1 >= nColMax) break;
847 }
848 nClusters++;
849
850 // Look for 5 pad cluster with minimum in the middle
851 Bool_t fivePadCluster = kFALSE;
852 if (col < (nColMax - 3)){
853 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
854 fivePadCluster = kTRUE;
855 }
856 if ((fivePadCluster) && (col < (nColMax - 5))) {
857 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh){
858 fivePadCluster = kFALSE;
859 }
860 }
861 if ((fivePadCluster) && (col > 1)){
862 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh){
863 fivePadCluster = kFALSE;
864 }
865 }
866 }
3fe61b77 867
f5375dcb 868 // 5 pad cluster
869 // Modify the signal of the overlapping pad for the left part
870 // of the cluster which remains from a previous unfolding
871 if (iUnfold) {
872 clusterSignal[0] *= ratioLeft;
873 iUnfold = 0;
874 }
3fe61b77 875
f5375dcb 876 // Unfold the 5 pad cluster
877 if (fivePadCluster){
878 for (iPad = 0; iPad < kNsig; iPad++) {
879 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
880 ,col-1+iPad
881 ,time));
882 }
883 // Unfold the two maxima and set the signal on
884 // the overlapping pad to the ratio
053767a4 885 ratioRight = Unfold(kEpsilon,ilayer,padSignal);
f5375dcb 886 ratioLeft = 1.0 - ratioRight;
887 clusterSignal[2] *= ratioRight;
888 iUnfold = 1;
889 }
3fe61b77 890
f5375dcb 891 // The position of the cluster in COL direction relative to the center pad (pad units)
892 Double_t clusterPosCol = 0.0;
a7ac01d2 893 if (AliTRDReconstructor::RecoParam()->IsLUT()) {
f5375dcb 894 // Calculate the position of the cluster by using the
895 // lookup table method
053767a4 896 clusterPosCol = LUTposition(ilayer,clusterSignal[0]
f5375dcb 897 ,clusterSignal[1]
898 ,clusterSignal[2]);
899 } else {
900 // Calculate the position of the cluster by using the
901 // center of gravity method
902 for (Int_t i = 0; i < kNsig; i++) {
903 padSignal[i] = 0.0;
904 }
905 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
906 padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left pad
907 padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right pad
908
909 if ((col > 2) &&
910 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
911 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
912 }
913 if ((col < nColMax - 3) &&
914 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])){
915 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
916 }
917 clusterPosCol = GetCOG(padSignal);
918 }
3fe61b77 919
f5375dcb 920 // Store the amplitudes of the pads in the cluster for later analysis
921 // and check whether one of these pads is masked in the database
922 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
923 for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
924 if ((jPad < 0) ||
925 (jPad >= nColMax-1)) {
926 continue;
927 }
928 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
929 }
3fe61b77 930
f5375dcb 931 // Transform the local cluster coordinates into calibrated
932 // space point positions defined in the local tracking system.
933 // Here the calibration for T0, Vdrift and ExB is applied as well.
934 Double_t clusterXYZ[6];
935 clusterXYZ[0] = clusterPosCol;
936 clusterXYZ[1] = clusterSignal[0];
937 clusterXYZ[2] = clusterSignal[1];
938 clusterXYZ[3] = clusterSignal[2];
939 clusterXYZ[4] = 0.0;
940 clusterXYZ[5] = 0.0;
941 Int_t clusterRCT[3];
942 clusterRCT[0] = row;
943 clusterRCT[1] = col;
944 clusterRCT[2] = 0;
945
946 Bool_t out = kTRUE;
947 if (fTransform->Transform(clusterXYZ, clusterRCT, ((UInt_t) time), out, 0)) {
948
949 // Add the cluster to the output array
950 // The track indices will be stored later
951 Float_t clusterPos[3];
952 clusterPos[0] = clusterXYZ[0];
953 clusterPos[1] = clusterXYZ[1];
954 clusterPos[2] = clusterXYZ[2];
955 Float_t clusterSig[2];
956 clusterSig[0] = clusterXYZ[4];
957 clusterSig[1] = clusterXYZ[5];
958 Double_t clusterCharge = clusterXYZ[3];
959 Char_t clusterTimeBin = ((Char_t) clusterRCT[2]);
960 AliTRDcluster *cluster = new AliTRDcluster(idet
961 ,clusterCharge
962 ,clusterPos
963 ,clusterSig
964 ,0x0
965 ,((Char_t) nPadCount)
966 ,signals
967 ,((UChar_t) col)
968 ,((UChar_t) row)
969 ,((UChar_t) time)
970 ,clusterTimeBin
971 ,clusterPosCol
972 ,volid);
973 cluster->SetInChamber(!out);
974 UChar_t maskPosition = padStatus.GetDataUnchecked(row, col, time);
975 if(maskPosition){
976 cluster->SetPadMaskedPosition(maskPosition);
977
978 if(maskPosition & AliTRDcluster::kMaskedLeft) cluster->SetPadMaskedStatus(status[0]);
979 else if(maskPosition & AliTRDcluster::kMaskedCenter) cluster->SetPadMaskedStatus(status[1]);
980 else cluster->SetPadMaskedStatus(status[2]);
981 }
3fe61b77 982
f5375dcb 983 // Temporarily store the row, column and time bin of the center pad
984 // Used to later on assign the track indices
985 cluster->SetLabel( row,0);
986 cluster->SetLabel( col,1);
987 cluster->SetLabel(time,2);
988
989 RecPoints()->Add(cluster);
3fe61b77 990
f5375dcb 991 // Store the index of the first cluster in the current ROC
992 if (firstClusterROC < 0) {
993 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
994 }
995
996 // Count the number of cluster in the current ROC
997 nClusterROC++;
998
999 } // if: Transform ok ?
1000 } // if: Maximum found ?
1001 }
3fe61b77 1002
1003 delete digitsOut;
1004
f5375dcb 1005 if (fAddLabels) AddLabels(idet, firstClusterROC, nClusterROC);
3fe61b77 1006
1007 // Write the cluster and reset the array
1008 WriteClusters(idet);
1009 ResetRecPoints();
1010
1011 return kTRUE;
1012
1013}
1014
1015//_____________________________________________________________________________
1016Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
3551db50 1017{
1018 //
3fe61b77 1019 // Add the track indices to the found clusters
3551db50 1020 //
1021
3fe61b77 1022 const Int_t kNclus = 3;
1023 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1024 const Int_t kNtrack = kNdict * kNclus;
1025
1026 Int_t iClusterROC = 0;
1027
1028 Int_t row = 0;
1029 Int_t col = 0;
1030 Int_t time = 0;
1031 Int_t iPad = 0;
1032
1033 // Temporary array to collect the track indices
1034 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1035
1036 // Loop through the dictionary arrays one-by-one
1037 // to keep memory consumption low
1038 AliTRDdataArrayI *tracksIn = 0;
1039 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1040
1041 // tracksIn should be expanded beforehand!
625f5260 1042 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
3fe61b77 1043
1044 // Loop though the clusters found in this ROC
1045 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1046
1047 AliTRDcluster *cluster = (AliTRDcluster *)
1048 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1049 row = cluster->GetLabel(0);
1050 col = cluster->GetLabel(1);
1051 time = cluster->GetLabel(2);
1052
1053 for (iPad = 0; iPad < kNclus; iPad++) {
1054 Int_t iPadCol = col - 1 + iPad;
1055 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1056 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1057 }
1058
1059 }
1060
1061 }
1062
1063 // Copy the track indices into the cluster
1064 // Loop though the clusters found in this ROC
1065 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1066
1067 AliTRDcluster *cluster = (AliTRDcluster *)
1068 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1069 cluster->SetLabel(-9999,0);
1070 cluster->SetLabel(-9999,1);
1071 cluster->SetLabel(-9999,2);
1072
1073 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1074
1075 }
1076
1077 delete [] idxTracks;
1078
1079 return kTRUE;
1080
1081}
1082
1083//_____________________________________________________________________________
acc49af9 1084Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
3fe61b77 1085{
1086 //
1087 // Get COG position
1088 // Used for clusters with more than 3 pads - where LUT not applicable
1089 //
1090
1091 Double_t sum = signal[0]
1092 + signal[1]
1093 + signal[2]
1094 + signal[3]
1095 + signal[4];
1096
1097 Double_t res = (0.0 * (-signal[0] + signal[4])
1098 + (-signal[1] + signal[3])) / sum;
1099
1100 return res;
1101
1102}
1103
1104//_____________________________________________________________________________
053767a4 1105Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal)
3fe61b77 1106{
1107 //
1108 // Method to unfold neighbouring maxima.
1109 // The charge ratio on the overlapping pad is calculated
1110 // until there is no more change within the range given by eps.
1111 // The resulting ratio is then returned to the calling method.
1112 //
1113
1114 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1115 if (!calibration) {
3fe61b77 1116 AliError("No AliTRDcalibDB instance available\n");
1117 return kFALSE;
6d50f529 1118 }
3fe61b77 1119
1120 Int_t irc = 0;
1121 Int_t itStep = 0; // Count iteration steps
1122
1123 Double_t ratio = 0.5; // Start value for ratio
1124 Double_t prevRatio = 0.0; // Store previous ratio
1125
1126 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1127 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1128 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1129
1130 // Start the iteration
1131 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1132
1133 itStep++;
1134 prevRatio = ratio;
1135
1136 // Cluster position according to charge ratio
1137 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1138 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1139 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1140 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1141
1142 // Set cluster charge ratio
053767a4 1143 irc = calibration->PadResponse(1.0,maxLeft ,layer,newSignal);
3fe61b77 1144 Double_t ampLeft = padSignal[1] / newSignal[1];
053767a4 1145 irc = calibration->PadResponse(1.0,maxRight,layer,newSignal);
3fe61b77 1146 Double_t ampRight = padSignal[3] / newSignal[1];
1147
1148 // Apply pad response to parameters
053767a4 1149 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1150 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
3fe61b77 1151
1152 // Calculate new overlapping ratio
1153 ratio = TMath::Min((Double_t) 1.0
1154 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1155
b43a3e17 1156 }
88719a08 1157
3fe61b77 1158 return ratio;
1159
1160}
88719a08 1161
3fe61b77 1162//_____________________________________________________________________________
fb44eb8e 1163void AliTRDclusterizer::TailCancelation(AliTRDdataArrayDigits *digitsIn
3fe61b77 1164 , AliTRDdataArrayF *digitsOut
1165 , AliTRDSignalIndex *indexesIn
1166 , AliTRDSignalIndex *indexesOut
1167 , Int_t nTimeTotal
acc49af9 1168 , Float_t adcThreshold
3fe61b77 1169 , AliTRDCalROC *calGainFactorROC
1170 , Float_t calGainFactorDetValue)
1171{
1172 //
1173 // Applies the tail cancelation and gain factors:
1174 // Transform digitsIn to digitsOut
1175 //
1176
1177 Int_t iRow = 0;
1178 Int_t iCol = 0;
1179 Int_t iTime = 0;
1180
3fe61b77 1181 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
1182 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
1183 indexesIn->ResetCounters();
1184 while (indexesIn->NextRCIndex(iRow, iCol))
1185 {
1186 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1187 Double_t gain = calGainFactorDetValue
1188 * calGainFactorROCValue;
1189
f5375dcb 1190 Bool_t corrupted = kFALSE;
3fe61b77 1191 for (iTime = 0; iTime < nTimeTotal; iTime++)
1192 {
1193 // Apply gain gain factor
1194 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
f5375dcb 1195 if(digitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
3fe61b77 1196 inADC[iTime] /= gain;
1197 outADC[iTime] = inADC[iTime];
1198 }
f5375dcb 1199 if(!corrupted)
1200 {
1201 // Apply the tail cancelation via the digital filter
1202 // (only for non-coorupted pads)
a7ac01d2 1203 if (AliTRDReconstructor::RecoParam()->IsTailCancelation())
f5375dcb 1204 {
1205 DeConvExp(inADC,outADC,nTimeTotal,AliTRDReconstructor::RecoParam()->GetTCnexp());
1206 }
1207 }
3fe61b77 1208
1209 indexesIn->ResetTbinCounter();
1210 while (indexesIn->NextTbinIndex(iTime))
1211 {
1212 // Store the amplitude of the digit if above threshold
acc49af9 1213 if (outADC[iTime] > adcThreshold)
3fe61b77 1214 {
1215 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1216 indexesOut->AddIndexTBin(iRow,iCol,iTime);
1217 }
1218 } // while itime
1219
1220 } // while irow icol
1221
1222 delete [] inADC;
1223 delete [] outADC;
1224
1225 return;
1226
1227}
1228
1229//_____________________________________________________________________________
1230void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1231 , Int_t n, Int_t nexp)
1232{
1233 //
1234 // Tail cancellation by deconvolution for PASA v4 TRF
1235 //
1236
1237 Double_t rates[2];
1238 Double_t coefficients[2];
1239
1240 // Initialization (coefficient = alpha, rates = lambda)
1241 Double_t r1 = 1.0;
1242 Double_t r2 = 1.0;
1243 Double_t c1 = 0.5;
1244 Double_t c2 = 0.5;
1245
1246 if (nexp == 1) { // 1 Exponentials
1247 r1 = 1.156;
1248 r2 = 0.130;
1249 c1 = 0.066;
1250 c2 = 0.000;
1251 }
1252 if (nexp == 2) { // 2 Exponentials
1253 r1 = 1.156;
1254 r2 = 0.130;
1255 c1 = 0.114;
1256 c2 = 0.624;
1257 }
1258
1259 coefficients[0] = c1;
1260 coefficients[1] = c2;
1261
1262 Double_t dt = 0.1;
1263
1264 rates[0] = TMath::Exp(-dt/(r1));
1265 rates[1] = TMath::Exp(-dt/(r2));
1266
1267 Int_t i = 0;
1268 Int_t k = 0;
1269
1270 Double_t reminder[2];
1271 Double_t correction = 0.0;
1272 Double_t result = 0.0;
1273
1274 // Attention: computation order is important
1275 for (k = 0; k < nexp; k++) {
1276 reminder[k] = 0.0;
1277 }
1278
1279 for (i = 0; i < n; i++) {
1280
1281 result = (source[i] - correction); // No rescaling
1282 target[i] = result;
1283
1284 for (k = 0; k < nexp; k++) {
1285 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1286 }
1287
1288 correction = 0.0;
1289 for (k = 0; k < nexp; k++) {
1290 correction += reminder[k];
1291 }
1292
1293 }
6d50f529 1294
1295}
1296
1297//_____________________________________________________________________________
1298void AliTRDclusterizer::ResetRecPoints()
1299{
1300 //
1301 // Resets the list of rec points
1302 //
1303
1304 if (fRecPoints) {
1305 fRecPoints->Delete();
1306 }
1307
1308}
1309
1310//_____________________________________________________________________________
34eaaa7e 1311TObjArray *AliTRDclusterizer::RecPoints()
6d50f529 1312{
1313 //
1314 // Returns the list of rec points
1315 //
1316
1317 if (!fRecPoints) {
1318 fRecPoints = new TObjArray(400);
1319 }
1320
1321 return fRecPoints;
1322
3551db50 1323}
fc546d21 1324
1325//_____________________________________________________________________________
1326void AliTRDclusterizer::FillLUT()
1327{
1328 //
1329 // Create the LUT
1330 //
1331
1332 const Int_t kNlut = 128;
1333
053767a4 1334 fLUTbin = AliTRDgeometry::kNlayer * kNlut;
fc546d21 1335
1336 // The lookup table from Bogdan
053767a4 1337 Float_t lut[AliTRDgeometry::kNlayer][kNlut] = {
fc546d21 1338 {
1339 0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1340 0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1341 0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1342 0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1343 0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1344 0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1345 0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1346 0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1347 0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1348 0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1349 0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1350 0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1351 0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1352 0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1353 0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1354 0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1355 },
1356 {
1357 0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1358 0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1359 0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1360 0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1361 0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1362 0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1363 0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1364 0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1365 0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1366 0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1367 0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1368 0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1369 0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1370 0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1371 0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1372 0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1373 },
1374 {
1375 0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1376 0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1377 0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1378 0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1379 0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1380 0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1381 0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1382 0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1383 0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1384 0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1385 0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1386 0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1387 0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1388 0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1389 0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1390 0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1391 },
1392 {
1393 0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1394 0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1395 0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1396 0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1397 0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1398 0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1399 0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1400 0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1401 0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1402 0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1403 0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1404 0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1405 0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1406 0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1407 0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1408 0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1409 },
1410 {
1411 0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1412 0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1413 0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1414 0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1415 0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1416 0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1417 0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1418 0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1419 0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1420 0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1421 0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1422 0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1423 0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1424 0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1425 0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1426 0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1427 },
1428 {
1429 0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1430 0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1431 0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1432 0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1433 0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1434 0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1435 0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1436 0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1437 0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1438 0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1439 0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1440 0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1441 0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1442 0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1443 0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1444 0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1445 }
1446 };
1447
1448 if (fLUT) {
1449 delete [] fLUT;
1450 }
1451 fLUT = new Double_t[fLUTbin];
1452
053767a4 1453 for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
fc546d21 1454 for (Int_t ilut = 0; ilut < kNlut; ilut++ ) {
053767a4 1455 fLUT[ilayer*kNlut+ilut] = lut[ilayer][ilut];
fc546d21 1456 }
1457 }
1458
1459}
1460
1461//_____________________________________________________________________________
053767a4 1462Double_t AliTRDclusterizer::LUTposition(Int_t ilayer, Double_t ampL
fc546d21 1463 , Double_t ampC, Double_t ampR) const
1464{
1465 //
1466 // Calculates the cluster position using the lookup table.
1467 // Method provided by Bogdan Vulpescu.
1468 //
1469
1470 const Int_t kNlut = 128;
1471
1472 Double_t pos;
1473 Double_t x = 0.0;
1474 Double_t xmin;
1475 Double_t xmax;
1476 Double_t xwid;
1477
1478 Int_t side = 0;
1479 Int_t ix;
1480
053767a4 1481 Double_t xMin[AliTRDgeometry::kNlayer] = { 0.006492, 0.006377, 0.006258
1482 , 0.006144, 0.006030, 0.005980 };
1483 Double_t xMax[AliTRDgeometry::kNlayer] = { 0.960351, 0.965870, 0.970445
1484 , 0.974352, 0.977667, 0.996101 };
fc546d21 1485
1486 if (ampL > ampR) {
1487 x = (ampL - ampR) / ampC;
1488 side = -1;
1489 }
1490 else if (ampL < ampR) {
1491 x = (ampR - ampL) / ampC;
1492 side = +1;
1493 }
1494
1495 if (ampL != ampR) {
1496
053767a4 1497 xmin = xMin[ilayer] + 0.000005;
1498 xmax = xMax[ilayer] - 0.000005;
fc546d21 1499 xwid = (xmax - xmin) / 127.0;
1500
1501 if (x < xmin) {
1502 pos = 0.0000;
1503 }
1504 else if (x > xmax) {
1505 pos = side * 0.5000;
1506 }
1507 else {
1508 ix = (Int_t) ((x - xmin) / xwid);
053767a4 1509 pos = side * fLUT[ilayer*kNlut+ix];
fc546d21 1510 }
1511
1512 }
1513 else {
1514
1515 pos = 0.0;
1516
1517 }
1518
1519 return pos;
1520
1521}