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