]> 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
- put baseline subtraction into the clusterizer and remove all other
[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 <TClonesArray.h>
25#include <TObjArray.h>
26
27#include "AliRunLoader.h"
28#include "AliLoader.h"
29#include "AliAlignObj.h"
30
31#include "AliTRDclusterizer.h"
32#include "AliTRDcluster.h"
33#include "AliTRDReconstructor.h"
34#include "AliTRDgeometry.h"
35#include "AliTRDarrayDictionary.h"
36#include "AliTRDarrayADC.h"
37#include "AliTRDdigitsManager.h"
38#include "AliTRDdigitsParam.h"
39#include "AliTRDrawData.h"
40#include "AliTRDcalibDB.h"
41#include "AliTRDtransform.h"
42#include "AliTRDSignalIndex.h"
43#include "AliTRDrawStreamBase.h"
44#include "AliTRDfeeParam.h"
45#include "AliTRDtrackletWord.h"
46
47#include "TTreeStream.h"
48
49#include "Cal/AliTRDCalROC.h"
50#include "Cal/AliTRDCalDet.h"
51#include "Cal/AliTRDCalSingleChamberStatus.h"
52
53ClassImp(AliTRDclusterizer)
54
55//_____________________________________________________________________________
56AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
57 :TNamed()
58 ,fReconstructor(rec)
59 ,fRunLoader(NULL)
60 ,fClusterTree(NULL)
61 ,fRecPoints(NULL)
62 ,fTracklets(NULL)
63 ,fTrackletTree(NULL)
64 ,fDigitsManager(new AliTRDdigitsManager())
65 ,fTrackletContainer(NULL)
66 ,fRawVersion(2)
67 ,fTransform(new AliTRDtransform(0))
68 ,fDigits(NULL)
69 ,fIndexes(NULL)
70 ,fMaxThresh(0)
71 ,fSigThresh(0)
72 ,fMinMaxCutSigma(0)
73 ,fMinLeftRightCutSigma(0)
74 ,fLayer(0)
75 ,fDet(0)
76 ,fVolid(0)
77 ,fColMax(0)
78 ,fTimeTotal(0)
79 ,fCalGainFactorROC(NULL)
80 ,fCalGainFactorDetValue(0)
81 ,fCalNoiseROC(NULL)
82 ,fCalNoiseDetValue(0)
83 ,fCalPadStatusROC(NULL)
84 ,fClusterROC(0)
85 ,firstClusterROC(0)
86 ,fNoOfClusters(0)
87 ,fBaseline(0)
88{
89 //
90 // AliTRDclusterizer default constructor
91 //
92
93 SetBit(kLabels, kTRUE);
94 SetBit(knewDM, kFALSE);
95
96 AliTRDcalibDB *trd = 0x0;
97 if (!(trd = AliTRDcalibDB::Instance())) {
98 AliFatal("Could not get calibration object");
99 }
100
101 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
102
103 // Initialize debug stream
104 if(fReconstructor){
105 if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 1){
106 TDirectory *savedir = gDirectory;
107 //fgGetDebugStream = new TTreeSRedirector("TRD.ClusterizerDebug.root");
108 savedir->cd();
109 }
110 }
111
112}
113
114//_____________________________________________________________________________
115AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, const AliTRDReconstructor *const rec)
116 :TNamed(name,title)
117 ,fReconstructor(rec)
118 ,fRunLoader(NULL)
119 ,fClusterTree(NULL)
120 ,fRecPoints(NULL)
121 ,fTracklets(NULL)
122 ,fTrackletTree(NULL)
123 ,fDigitsManager(new AliTRDdigitsManager())
124 ,fTrackletContainer(NULL)
125 ,fRawVersion(2)
126 ,fTransform(new AliTRDtransform(0))
127 ,fDigits(NULL)
128 ,fIndexes(NULL)
129 ,fMaxThresh(0)
130 ,fSigThresh(0)
131 ,fMinMaxCutSigma(0)
132 ,fMinLeftRightCutSigma(0)
133 ,fLayer(0)
134 ,fDet(0)
135 ,fVolid(0)
136 ,fColMax(0)
137 ,fTimeTotal(0)
138 ,fCalGainFactorROC(NULL)
139 ,fCalGainFactorDetValue(0)
140 ,fCalNoiseROC(NULL)
141 ,fCalNoiseDetValue(0)
142 ,fCalPadStatusROC(NULL)
143 ,fClusterROC(0)
144 ,firstClusterROC(0)
145 ,fNoOfClusters(0)
146 ,fBaseline(0)
147{
148 //
149 // AliTRDclusterizer constructor
150 //
151
152 SetBit(kLabels, kTRUE);
153 SetBit(knewDM, kFALSE);
154
155 AliTRDcalibDB *trd = 0x0;
156 if (!(trd = AliTRDcalibDB::Instance())) {
157 AliFatal("Could not get calibration object");
158 }
159
160 fDigitsManager->CreateArrays();
161
162 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
163
164 //FillLUT();
165
166}
167
168//_____________________________________________________________________________
169AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
170 :TNamed(c)
171 ,fReconstructor(c.fReconstructor)
172 ,fRunLoader(NULL)
173 ,fClusterTree(NULL)
174 ,fRecPoints(NULL)
175 ,fTracklets(NULL)
176 ,fTrackletTree(NULL)
177 ,fDigitsManager(NULL)
178 ,fTrackletContainer(NULL)
179 ,fRawVersion(2)
180 ,fTransform(NULL)
181 ,fDigits(NULL)
182 ,fIndexes(NULL)
183 ,fMaxThresh(0)
184 ,fSigThresh(0)
185 ,fMinMaxCutSigma(0)
186 ,fMinLeftRightCutSigma(0)
187 ,fLayer(0)
188 ,fDet(0)
189 ,fVolid(0)
190 ,fColMax(0)
191 ,fTimeTotal(0)
192 ,fCalGainFactorROC(NULL)
193 ,fCalGainFactorDetValue(0)
194 ,fCalNoiseROC(NULL)
195 ,fCalNoiseDetValue(0)
196 ,fCalPadStatusROC(NULL)
197 ,fClusterROC(0)
198 ,firstClusterROC(0)
199 ,fNoOfClusters(0)
200 ,fBaseline(0)
201{
202 //
203 // AliTRDclusterizer copy constructor
204 //
205
206 SetBit(kLabels, kTRUE);
207 SetBit(knewDM, kFALSE);
208
209 //FillLUT();
210
211}
212
213//_____________________________________________________________________________
214AliTRDclusterizer::~AliTRDclusterizer()
215{
216 //
217 // AliTRDclusterizer destructor
218 //
219
220 if (fRecPoints/* && IsClustersOwner()*/){
221 fRecPoints->Delete();
222 delete fRecPoints;
223 }
224
225 if (fTracklets){
226 fTracklets->Delete();
227 delete fTracklets;
228 }
229
230 if (fDigitsManager) {
231 delete fDigitsManager;
232 fDigitsManager = NULL;
233 }
234
235 if (fTrackletContainer){
236 delete [] fTrackletContainer[0];
237 delete [] fTrackletContainer[1];
238 delete [] fTrackletContainer;
239 fTrackletContainer = NULL;
240 }
241
242 if (fTransform){
243 delete fTransform;
244 fTransform = NULL;
245 }
246
247}
248
249//_____________________________________________________________________________
250AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
251{
252 //
253 // Assignment operator
254 //
255
256 if (this != &c)
257 {
258 ((AliTRDclusterizer &) c).Copy(*this);
259 }
260
261 return *this;
262
263}
264
265//_____________________________________________________________________________
266void AliTRDclusterizer::Copy(TObject &c) const
267{
268 //
269 // Copy function
270 //
271
272 ((AliTRDclusterizer &) c).fClusterTree = NULL;
273 ((AliTRDclusterizer &) c).fRecPoints = NULL;
274 ((AliTRDclusterizer &) c).fTrackletTree = NULL;
275 ((AliTRDclusterizer &) c).fDigitsManager = NULL;
276 ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
277 ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
278 ((AliTRDclusterizer &) c).fTransform = NULL;
279 ((AliTRDclusterizer &) c).fDigits = NULL;
280 ((AliTRDclusterizer &) c).fIndexes = NULL;
281 ((AliTRDclusterizer &) c).fMaxThresh = 0;
282 ((AliTRDclusterizer &) c).fSigThresh = 0;
283 ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0;
284 ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0;
285 ((AliTRDclusterizer &) c).fLayer = 0;
286 ((AliTRDclusterizer &) c).fDet = 0;
287 ((AliTRDclusterizer &) c).fVolid = 0;
288 ((AliTRDclusterizer &) c).fColMax = 0;
289 ((AliTRDclusterizer &) c).fTimeTotal = 0;
290 ((AliTRDclusterizer &) c).fCalGainFactorROC = NULL;
291 ((AliTRDclusterizer &) c).fCalGainFactorDetValue = 0;
292 ((AliTRDclusterizer &) c).fCalNoiseROC = NULL;
293 ((AliTRDclusterizer &) c).fCalNoiseDetValue = 0;
294 ((AliTRDclusterizer &) c).fCalPadStatusROC = NULL;
295 ((AliTRDclusterizer &) c).fClusterROC = 0;
296 ((AliTRDclusterizer &) c).firstClusterROC= 0;
297 ((AliTRDclusterizer &) c).fNoOfClusters = 0;
298 ((AliTRDclusterizer &) c).fBaseline = 0;
299
300}
301
302//_____________________________________________________________________________
303Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
304{
305 //
306 // Opens the AliROOT file. Output and input are in the same file
307 //
308
309 TString evfoldname = AliConfig::GetDefaultEventFolderName();
310 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
311
312 if (!fRunLoader) {
313 fRunLoader = AliRunLoader::Open(name);
314 }
315
316 if (!fRunLoader) {
317 AliError(Form("Can not open session for file %s.",name));
318 return kFALSE;
319 }
320
321 OpenInput(nEvent);
322 OpenOutput();
323
324 return kTRUE;
325
326}
327
328//_____________________________________________________________________________
329Bool_t AliTRDclusterizer::OpenOutput()
330{
331 //
332 // Open the output file
333 //
334
335 if (!fReconstructor->IsWritingClusters()) return kTRUE;
336
337 TObjArray *ioArray = 0x0;
338
339 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
340 loader->MakeTree("R");
341
342 fClusterTree = loader->TreeR();
343 fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
344
345 return kTRUE;
346
347}
348
349//_____________________________________________________________________________
350Bool_t AliTRDclusterizer::OpenOutput(TTree *const clusterTree)
351{
352 //
353 // Connect the output tree
354 //
355
356 // clusters writing
357 if (fReconstructor->IsWritingClusters()){
358 TObjArray *ioArray = 0x0;
359 fClusterTree = clusterTree;
360 fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
361 }
362 return kTRUE;
363}
364
365//_____________________________________________________________________________
366Bool_t AliTRDclusterizer::OpenTrackletOutput()
367{
368 //
369 // Tracklet writing
370 //
371
372 if (fReconstructor->IsWritingTracklets()){
373 TString evfoldname = AliConfig::GetDefaultEventFolderName();
374 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
375
376 if (!fRunLoader) {
377 fRunLoader = AliRunLoader::Open("galice.root");
378 }
379 if (!fRunLoader) {
380 AliError(Form("Can not open session for file galice.root."));
381 return kFALSE;
382 }
383
384 UInt_t **leaves = new UInt_t *[2];
385 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
386 if (!dl) {
387 AliError("Could not get the tracklets data loader!");
388 dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
389 fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
390 }
391 fTrackletTree = dl->Tree();
392 if (!fTrackletTree)
393 {
394 dl->MakeTree();
395 fTrackletTree = dl->Tree();
396 }
397 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
398 if (!trkbranch)
399 fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
400 }
401
402 return kTRUE;
403}
404
405//_____________________________________________________________________________
406Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
407{
408 //
409 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
410 //
411
412 // Import the Trees for the event nEvent in the file
413 fRunLoader->GetEvent(nEvent);
414
415 return kTRUE;
416
417}
418
419//_____________________________________________________________________________
420Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
421{
422 //
423 // Fills TRDcluster branch in the tree with the clusters
424 // found in detector = det. For det=-1 writes the tree.
425 //
426
427 if ((det < -1) ||
428 (det >= AliTRDgeometry::Ndet())) {
429 AliError(Form("Unexpected detector index %d.\n",det));
430 return kFALSE;
431 }
432
433 TObjArray *ioArray = new TObjArray(400);
434 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
435 if (!branch) {
436 branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
437 } else branch->SetAddress(&ioArray);
438
439 Int_t nRecPoints = RecPoints()->GetEntriesFast();
440 if(det >= 0){
441 for (Int_t i = 0; i < nRecPoints; i++) {
442 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
443 if(det != c->GetDetector()) continue;
444 ioArray->AddLast(c);
445 }
446 fClusterTree->Fill();
447 } else {
448
449 Int_t detOld = -1;
450 for (Int_t i = 0; i < nRecPoints; i++) {
451 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
452 if(c->GetDetector() != detOld){
453 fClusterTree->Fill();
454 ioArray->Clear();
455 detOld = c->GetDetector();
456 }
457 ioArray->AddLast(c);
458 }
459 }
460 delete ioArray;
461
462 return kTRUE;
463
464}
465
466//_____________________________________________________________________________
467Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
468{
469 //
470 // Write the raw data tracklets into seperate file
471 //
472
473 UInt_t **leaves = new UInt_t *[2];
474 for (Int_t i=0; i<2 ;i++){
475 leaves[i] = new UInt_t[258];
476 leaves[i][0] = det; // det
477 leaves[i][1] = i; // side
478 memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
479 }
480
481 if (!fTrackletTree){
482 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
483 dl->MakeTree();
484 fTrackletTree = dl->Tree();
485 }
486
487 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
488 if (!trkbranch) {
489 trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
490 }
491
492 for (Int_t i=0; i<2; i++){
493 if (leaves[i][2]>0) {
494 trkbranch->SetAddress(leaves[i]);
495 fTrackletTree->Fill();
496 }
497 }
498
499 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
500 dl->WriteData("OVERWRITE");
501 //dl->Unload();
502 delete [] leaves;
503
504 return kTRUE;
505
506}
507
508//_____________________________________________________________________________
509Bool_t AliTRDclusterizer::ReadDigits()
510{
511 //
512 // Reads the digits arrays from the input aliroot file
513 //
514
515 if (!fRunLoader) {
516 AliError("No run loader available");
517 return kFALSE;
518 }
519
520 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
521 if (!loader->TreeD()) {
522 loader->LoadDigits();
523 }
524
525 // Read in the digit arrays
526 return (fDigitsManager->ReadDigits(loader->TreeD()));
527
528}
529
530//_____________________________________________________________________________
531Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
532{
533 //
534 // Reads the digits arrays from the input tree
535 //
536
537 // Read in the digit arrays
538 return (fDigitsManager->ReadDigits(digitsTree));
539
540}
541
542//_____________________________________________________________________________
543Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
544{
545 //
546 // Reads the digits arrays from the ddl file
547 //
548
549 AliTRDrawData raw;
550 fDigitsManager = raw.Raw2Digits(rawReader);
551
552 return kTRUE;
553
554}
555
556//_____________________________________________________________________________
557Bool_t AliTRDclusterizer::MakeClusters()
558{
559 //
560 // Creates clusters from digits
561 //
562
563 // Propagate info from the digits manager
564 if (TestBit(kLabels)){
565 SetBit(kLabels, fDigitsManager->UsesDictionaries());
566 }
567
568 Bool_t fReturn = kTRUE;
569 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
570
571 AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
572 // This is to take care of switched off super modules
573 if (!digitsIn->HasData()) continue;
574 digitsIn->Expand();
575 digitsIn->DeleteNegatives(); // Restore digits array to >=0 values
576 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
577 if (indexes->IsAllocated() == kFALSE){
578 fDigitsManager->BuildIndexes(i);
579 }
580
581 Bool_t fR = kFALSE;
582 if (indexes->HasEntry()){
583 if (TestBit(kLabels)){
584 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
585 AliTRDarrayDictionary *tracksIn = 0; //mod
586 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict); //mod
587 tracksIn->Expand();
588 }
589 }
590 fR = MakeClusters(i);
591 fReturn = fR && fReturn;
592 }
593
594 //if (fR == kFALSE){
595 // if(IsWritingClusters()) WriteClusters(i);
596 // ResetRecPoints();
597 //}
598
599 // No compress just remove
600 fDigitsManager->RemoveDigits(i);
601 fDigitsManager->RemoveDictionaries(i);
602 fDigitsManager->ClearIndexes(i);
603 }
604
605 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
606
607 AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
608
609 return fReturn;
610
611}
612
613//_____________________________________________________________________________
614Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
615{
616 //
617 // Creates clusters from raw data
618 //
619
620 return Raw2ClustersChamber(rawReader);
621
622}
623
624//_____________________________________________________________________________
625Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
626{
627 //
628 // Creates clusters from raw data
629 //
630
631 // Create the digits manager
632 if (!fDigitsManager){
633 SetBit(knewDM, kTRUE);
634 fDigitsManager = new AliTRDdigitsManager(kTRUE);
635 fDigitsManager->CreateArrays();
636 }
637
638 fDigitsManager->SetUseDictionaries(TestBit(kLabels));
639
640 // tracklet container for raw tracklet writing
641 if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
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 AliTRDrawStreamBase *input = AliTRDrawStreamBase::GetRawStream(rawReader);
650 if(fReconstructor->IsHLT())
651 input->SetSharedPadReadout(kFALSE);
652
653 AliInfo(Form("Stream version: %s", input->IsA()->GetName()));
654
655 Int_t det = 0;
656 while ((det = input->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
657 Bool_t iclusterBranch = kFALSE;
658 if (fDigitsManager->GetIndexes(det)->HasEntry()){
659 iclusterBranch = MakeClusters(det);
660 }
661
662 fDigitsManager->ClearArrays(det);
663
664 if (!fReconstructor->IsWritingTracklets()) continue;
665 if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
666 }
667
668 if (fTrackletContainer){
669 delete [] fTrackletContainer[0];
670 delete [] fTrackletContainer[1];
671 delete [] fTrackletContainer;
672 fTrackletContainer = NULL;
673 }
674
675 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
676
677 if(!TestBit(knewDM)){
678 delete fDigitsManager;
679 fDigitsManager = NULL;
680 }
681
682 delete input;
683 input = NULL;
684
685 AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
686 return kTRUE;
687
688}
689
690//_____________________________________________________________________________
691UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
692{
693 //
694 // Check if a pad is masked
695 //
696
697 UChar_t status = 0;
698
699 if(signal>0 && TESTBIT(signal, 10)){
700 CLRBIT(signal, 10);
701 for(int ibit=0; ibit<4; ibit++){
702 if(TESTBIT(signal, 11+ibit)){
703 SETBIT(status, ibit);
704 CLRBIT(signal, 11+ibit);
705 }
706 }
707 }
708 return status;
709}
710
711//_____________________________________________________________________________
712void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
713 //
714 // Set the pad status into out
715 // First three bits are needed for the position encoding
716 //
717 out |= status << 3;
718}
719
720//_____________________________________________________________________________
721UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
722 //
723 // return the staus encoding of the corrupted pad
724 //
725 return static_cast<UChar_t>(encoding >> 3);
726}
727
728//_____________________________________________________________________________
729Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
730 //
731 // Return the position of the corruption
732 //
733 return encoding & 7;
734}
735
736//_____________________________________________________________________________
737Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
738{
739 //
740 // Generates the cluster.
741 //
742
743 // Get the digits
744 fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
745 fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline();
746
747 // This is to take care of switched off super modules
748 if (!fDigits->HasData()) return kFALSE;
749
750 fIndexes = fDigitsManager->GetIndexes(det);
751 if (fIndexes->IsAllocated() == kFALSE) {
752 AliError("Indexes do not exist!");
753 return kFALSE;
754 }
755
756 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
757 if (!calibration) {
758 AliFatal("No AliTRDcalibDB instance available\n");
759 return kFALSE;
760 }
761
762 if (!fReconstructor){
763 AliError("Reconstructor not set\n");
764 return kFALSE;
765 }
766
767 fMaxThresh = fReconstructor->GetRecoParam()->GetClusMaxThresh();
768 fSigThresh = fReconstructor->GetRecoParam()->GetClusSigThresh();
769 fMinMaxCutSigma = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
770 fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
771
772 Int_t istack = fIndexes->GetStack();
773 fLayer = fIndexes->GetLayer();
774 Int_t isector = fIndexes->GetSM();
775
776 // Start clustering in the chamber
777
778 fDet = AliTRDgeometry::GetDetector(fLayer,istack,isector);
779 if (fDet != det) {
780 AliError("Strange Detector number Missmatch!");
781 return kFALSE;
782 }
783
784 // TRD space point transformation
785 fTransform->SetDetector(det);
786
787 Int_t iGeoLayer = AliGeomManager::kTRD1 + fLayer;
788 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
789 fVolid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
790
791 if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
792 AddTrackletsToArray();
793
794 fColMax = fDigits->GetNcol();
795 //Int_t nRowMax = fDigits->GetNrow();
796 fTimeTotal = fDigits->GetNtime();
797
798 // Detector wise calibration object for the gain factors
799 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
800 // Calibration object with pad wise values for the gain factors
801 fCalGainFactorROC = calibration->GetGainFactorROC(fDet);
802 // Calibration value for chamber wise gain factor
803 fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
804
805 // Detector wise calibration object for the noise
806 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
807 // Calibration object with pad wise values for the noise
808 fCalNoiseROC = calibration->GetNoiseROC(fDet);
809 // Calibration value for chamber wise noise
810 fCalNoiseDetValue = calNoiseDet->GetValue(fDet);
811
812 // Calibration object with the pad status
813 fCalPadStatusROC = calibration->GetPadStatusROC(fDet);
814
815 SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
816 SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
817 SetBit(kHLT, fReconstructor->IsHLT());
818
819 firstClusterROC = -1;
820 fClusterROC = 0;
821
822 // Apply the gain and the tail cancelation via digital filter
823 if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
824
825 MaxStruct curr, last;
826 Int_t nMaximas = 0, nCorrupted = 0;
827
828 // Here the clusterfining is happening
829
830 for(curr.time = 0; curr.time < fTimeTotal; curr.time++){
831 while(fIndexes->NextRCIndex(curr.row, curr.col)){
832 if(IsMaximum(curr, curr.padStatus, &curr.signals[0])){
833 if(last.row>-1){
834 if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
835 CreateCluster(last);
836 }
837 last=curr; curr.fivePad=kFALSE;
838 }
839 }
840 }
841 if(last.row>-1) CreateCluster(last);
842
843 if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
844 TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
845 (*fDebugStream) << "MakeClusters"
846 << "Detector=" << det
847 << "NMaxima=" << nMaximas
848 << "NClusters=" << fClusterROC
849 << "NCorrupted=" << nCorrupted
850 << "\n";
851 }
852 if (TestBit(kLabels)) AddLabels();
853
854 return kTRUE;
855
856}
857
858//_____________________________________________________________________________
859Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
860{
861 //
862 // Returns true if this row,col,time combination is a maximum.
863 // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
864 //
865
866 Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
867 Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
868 if(Signals[1] < fMaxThresh) return kFALSE;
869
870 Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
871 if (Signals[1] < noiseMiddleThresh) return kFALSE;
872
873 if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
874
875 UChar_t status[3]={
876 fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
877 ,fCalPadStatusROC->GetStatus(Max.col, Max.row)
878 ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
879 };
880
881gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
882 Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
883 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
884 Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
885
886 if(!(status[0] | status[1] | status[2])) {//all pads are good
887 if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) {
888 if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
889 Float_t noiseSumThresh = fMinLeftRightCutSigma
890 * fCalNoiseDetValue
891 * fCalNoiseROC->GetValue(Max.col, Max.row);
892 if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
893 padStatus = 0;
894 return kTRUE;
895 }
896 }
897 } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
898 if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
899 Signals[2]=0;
900 SetPadStatus(status[2], padStatus);
901 return kTRUE;
902 }
903 else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
904 Signals[0]=0;
905 SetPadStatus(status[0], padStatus);
906 return kTRUE;
907 }
908 else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
909 Signals[1] = (Short_t)(fMaxThresh + 0.5f);
910 SetPadStatus(status[1], padStatus);
911 return kTRUE;
912 }
913 }
914 return kFALSE;
915}
916
917//_____________________________________________________________________________
918Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
919{
920 //
921 // Look for 5 pad cluster with minimum in the middle
922 // Gives back the ratio
923 //
924
925 if (ThisMax.col >= fColMax - 3) return kFALSE;
926 Float_t gain;
927 if (ThisMax.col < fColMax - 5){
928 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
929 if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
930 return kFALSE;
931 }
932 if (ThisMax.col > 1) {
933 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
934 if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
935 return kFALSE;
936 }
937
938 const Float_t kEpsilon = 0.01;
939 Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
940 NeighbourMax.signals[1], NeighbourMax.signals[2]};
941
942 // Unfold the two maxima and set the signal on
943 // the overlapping pad to the ratio
944 Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
945 ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
946 NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
947 ThisMax.fivePad=kTRUE;
948 NeighbourMax.fivePad=kTRUE;
949 return kTRUE;
950
951}
952
953//_____________________________________________________________________________
954void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
955{
956 //
957 // Creates a cluster at the given position and saves it in fRecPoints
958 //
959
960 Int_t nPadCount = 1;
961 Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
962 if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
963
964 AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
965 cluster.SetNPads(nPadCount);
966 if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
967 else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
968 else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
969
970 cluster.SetFivePad(Max.fivePad);
971 // set pads status for the cluster
972 UChar_t maskPosition = GetCorruption(Max.padStatus);
973 if (maskPosition) {
974 cluster.SetPadMaskedPosition(maskPosition);
975 cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
976 }
977
978 // Transform the local cluster coordinates into calibrated
979 // space point positions defined in the local tracking system.
980 // Here the calibration for T0, Vdrift and ExB is applied as well.
981 if(!fTransform->Transform(&cluster)) return;
982 // Temporarily store the Max.Row, column and time bin of the center pad
983 // Used to later on assign the track indices
984 cluster.SetLabel(Max.row, 0);
985 cluster.SetLabel(Max.col, 1);
986 cluster.SetLabel(Max.time,2);
987
988 //needed for HLT reconstruction
989 AddClusterToArray(&cluster);
990
991 // Store the index of the first cluster in the current ROC
992 if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
993
994 fNoOfClusters++;
995 fClusterROC++;
996}
997
998//_____________________________________________________________________________
999void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1000{
1001 // Look to the right
1002 Int_t ii = 1;
1003 while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
1004 nPadCount++;
1005 ii++;
1006 if (Max.col < ii) break;
1007 }
1008 // Look to the left
1009 ii = 1;
1010 while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
1011 nPadCount++;
1012 ii++;
1013 if (Max.col+ii >= fColMax) break;
1014 }
1015
1016 // Store the amplitudes of the pads in the cluster for later analysis
1017 // and check whether one of these pads is masked in the database
1018 signals[2]=Max.signals[0];
1019 signals[3]=Max.signals[1];
1020 signals[4]=Max.signals[2];
1021 Float_t gain;
1022 for(Int_t i = 0; i<2; i++)
1023 {
1024 if(Max.col+i >= 3){
1025 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1026 signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1027 }
1028 if(Max.col+3-i < fColMax){
1029 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1030 signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1031 }
1032 }
1033 /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1034 if ((jPad >= 0) && (jPad < fColMax))
1035 signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1036 }*/
1037}
1038
1039//_____________________________________________________________________________
1040void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1041{
1042 //
1043 // Add a cluster to the array
1044 //
1045
1046 Int_t n = RecPoints()->GetEntriesFast();
1047 if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1048 new((*RecPoints())[n]) AliTRDcluster(*cluster);
1049}
1050
1051//_____________________________________________________________________________
1052void AliTRDclusterizer::AddTrackletsToArray()
1053{
1054 //
1055 // Add the online tracklets of this chamber to the array
1056 //
1057
1058 UInt_t* trackletword;
1059 for(Int_t side=0; side<2; side++)
1060 {
1061 Int_t trkl=0;
1062 trackletword=fTrackletContainer[side];
1063 while(trackletword[trkl]>0){
1064 Int_t n = TrackletsArray()->GetEntriesFast();
1065 AliTRDtrackletWord tmp(trackletword[trkl]);
1066 new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1067 trkl++;
1068 }
1069 }
1070}
1071
1072//_____________________________________________________________________________
1073Bool_t AliTRDclusterizer::AddLabels()
1074{
1075 //
1076 // Add the track indices to the found clusters
1077 //
1078
1079 const Int_t kNclus = 3;
1080 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1081 const Int_t kNtrack = kNdict * kNclus;
1082
1083 Int_t iClusterROC = 0;
1084
1085 Int_t row = 0;
1086 Int_t col = 0;
1087 Int_t time = 0;
1088 Int_t iPad = 0;
1089
1090 // Temporary array to collect the track indices
1091 Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1092
1093 // Loop through the dictionary arrays one-by-one
1094 // to keep memory consumption low
1095 AliTRDarrayDictionary *tracksIn = 0; //mod
1096 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1097
1098 // tracksIn should be expanded beforehand!
1099 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1100
1101 // Loop though the clusters found in this ROC
1102 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1103
1104 AliTRDcluster *cluster = (AliTRDcluster *)
1105 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1106 row = cluster->GetLabel(0);
1107 col = cluster->GetLabel(1);
1108 time = cluster->GetLabel(2);
1109
1110 for (iPad = 0; iPad < kNclus; iPad++) {
1111 Int_t iPadCol = col - 1 + iPad;
1112 Int_t index = tracksIn->GetData(row,iPadCol,time); //Modification of -1 in Track
1113 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1114 }
1115
1116 }
1117
1118 }
1119
1120 // Copy the track indices into the cluster
1121 // Loop though the clusters found in this ROC
1122 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1123
1124 AliTRDcluster *cluster = (AliTRDcluster *)
1125 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1126 cluster->SetLabel(-9999,0);
1127 cluster->SetLabel(-9999,1);
1128 cluster->SetLabel(-9999,2);
1129
1130 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1131
1132 }
1133
1134 delete [] idxTracks;
1135
1136 return kTRUE;
1137
1138}
1139
1140//_____________________________________________________________________________
1141Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1142{
1143 //
1144 // Method to unfold neighbouring maxima.
1145 // The charge ratio on the overlapping pad is calculated
1146 // until there is no more change within the range given by eps.
1147 // The resulting ratio is then returned to the calling method.
1148 //
1149
1150 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1151 if (!calibration) {
1152 AliError("No AliTRDcalibDB instance available\n");
1153 return kFALSE;
1154 }
1155
1156 Int_t irc = 0;
1157 Int_t itStep = 0; // Count iteration steps
1158
1159 Double_t ratio = 0.5; // Start value for ratio
1160 Double_t prevRatio = 0.0; // Store previous ratio
1161
1162 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1163 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1164 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1165
1166 // Start the iteration
1167 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1168
1169 itStep++;
1170 prevRatio = ratio;
1171
1172 // Cluster position according to charge ratio
1173 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1174 / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1175 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1176 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1177
1178 // Set cluster charge ratio
1179 irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1180 Double_t ampLeft = padSignal[1] / newSignal[1];
1181 irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1182 Double_t ampRight = padSignal[3] / newSignal[1];
1183
1184 // Apply pad response to parameters
1185 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1186 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1187
1188 // Calculate new overlapping ratio
1189 ratio = TMath::Min((Double_t) 1.0
1190 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1191
1192 }
1193
1194 return ratio;
1195
1196}
1197
1198//_____________________________________________________________________________
1199void AliTRDclusterizer::TailCancelation()
1200{
1201 //
1202 // Applies the tail cancelation and gain factors:
1203 // Transform fDigits to fDigits
1204 //
1205
1206 Int_t iRow = 0;
1207 Int_t iCol = 0;
1208 Int_t iTime = 0;
1209
1210 Double_t *inADC = new Double_t[fTimeTotal]; // ADC data before tail cancellation
1211 Double_t *outADC = new Double_t[fTimeTotal]; // ADC data after tail cancellation
1212
1213 fIndexes->ResetCounters();
1214 TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1215 while(fIndexes->NextRCIndex(iRow, iCol))
1216 {
1217 Bool_t corrupted = kFALSE;
1218 for (iTime = 0; iTime < fTimeTotal; iTime++)
1219 {
1220 // Apply gain gain factor
1221 inADC[iTime] = fDigits->GetData(iRow,iCol,iTime);
1222 if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1223 outADC[iTime] = inADC[iTime];
1224 if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming()){
1225 (*fDebugStream) << "TailCancellation"
1226 << "col=" << iCol
1227 << "row=" << iRow
1228 << "time=" << iTime
1229 << "inADC=" << inADC[iTime]
1230 << "outADC=" << outADC[iTime]
1231 << "corrupted=" << corrupted
1232 << "\n";
1233 }
1234 }
1235 if (!corrupted)
1236 {
1237 // Apply the tail cancelation via the digital filter
1238 // (only for non-coorupted pads)
1239 DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1240 }
1241
1242 for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime))
1243 {
1244 // Store the amplitude of the digit if above threshold
1245 if (outADC[iTime] > 0)
1246 fDigits->SetData(iRow,iCol,iTime,TMath::Nint(outADC[iTime]));
1247 else
1248 fDigits->SetData(iRow,iCol,iTime,0);
1249 } // while itime
1250
1251 } // while irow icol
1252
1253 delete [] inADC;
1254 delete [] outADC;
1255
1256 return;
1257
1258}
1259
1260//_____________________________________________________________________________
1261void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1262 ,const Int_t n, const Int_t nexp)
1263{
1264 //
1265 // Tail cancellation by deconvolution for PASA v4 TRF
1266 //
1267
1268 Double_t rates[2];
1269 Double_t coefficients[2];
1270
1271 // Initialization (coefficient = alpha, rates = lambda)
1272 Double_t r1 = 1.0;
1273 Double_t r2 = 1.0;
1274 Double_t c1 = 0.5;
1275 Double_t c2 = 0.5;
1276
1277 if (nexp == 1) { // 1 Exponentials
1278 r1 = 1.156;
1279 r2 = 0.130;
1280 c1 = 0.066;
1281 c2 = 0.000;
1282 }
1283 if (nexp == 2) { // 2 Exponentials
1284 Double_t par[4];
1285 fReconstructor->GetRecoParam()->GetTCParams(par);
1286 r1 = par[0];//1.156;
1287 r2 = par[1];//0.130;
1288 c1 = par[2];//0.114;
1289 c2 = par[3];//0.624;
1290 }
1291
1292 coefficients[0] = c1;
1293 coefficients[1] = c2;
1294
1295 Double_t dt = 0.1;
1296
1297 rates[0] = TMath::Exp(-dt/(r1));
1298 rates[1] = TMath::Exp(-dt/(r2));
1299
1300 Int_t i = 0;
1301 Int_t k = 0;
1302
1303 Double_t reminder[2];
1304 Double_t correction = 0.0;
1305 Double_t result = 0.0;
1306
1307 // Attention: computation order is important
1308 for (k = 0; k < nexp; k++) {
1309 reminder[k] = 0.0;
1310 }
1311
1312 for (i = 0; i < n; i++) {
1313
1314 result = (source[i] - correction); // No rescaling
1315 target[i] = result;
1316
1317 for (k = 0; k < nexp; k++) {
1318 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1319 }
1320
1321 correction = 0.0;
1322 for (k = 0; k < nexp; k++) {
1323 correction += reminder[k];
1324 }
1325
1326 }
1327
1328}
1329
1330//_____________________________________________________________________________
1331void AliTRDclusterizer::ResetRecPoints()
1332{
1333 //
1334 // Resets the list of rec points
1335 //
1336
1337 if (fRecPoints) {
1338 fRecPoints->Delete();
1339 delete fRecPoints;
1340 }
1341}
1342
1343//_____________________________________________________________________________
1344TClonesArray *AliTRDclusterizer::RecPoints()
1345{
1346 //
1347 // Returns the list of rec points
1348 //
1349
1350 if (!fRecPoints) {
1351 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1352 // determine number of clusters which has to be allocated
1353 Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1354
1355 fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1356 }
1357 //SetClustersOwner(kTRUE);
1358 AliTRDReconstructor::SetClusters(0x0);
1359 }
1360 return fRecPoints;
1361
1362}
1363
1364//_____________________________________________________________________________
1365TClonesArray *AliTRDclusterizer::TrackletsArray()
1366{
1367 //
1368 // Returns the list of rec points
1369 //
1370
1371 if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1372 fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1373 //SetClustersOwner(kTRUE);
1374 //AliTRDReconstructor::SetTracklets(0x0);
1375 }
1376 return fTracklets;
1377
1378}
1379