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