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