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