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