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