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