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