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