Faster version of baseline correction by Theodor with additional protection against...
[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   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod     
732   
733   // This is to take care of switched off super modules
734   if (!fDigits->HasData()) return kFALSE;
735
736   fIndexes = fDigitsManager->GetIndexes(det);
737   if (fIndexes->IsAllocated() == kFALSE) {
738     AliError("Indexes do not exist!");
739     return kFALSE;      
740   }
741
742   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
743   if (!calibration) {
744     AliFatal("No AliTRDcalibDB instance available\n");
745     return kFALSE;  
746   }
747
748   if (!fReconstructor){
749     AliError("Reconstructor not set\n");
750     return kFALSE;
751   }
752
753   fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
754   fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
755   fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
756   fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
757
758   Int_t istack  = fIndexes->GetStack();
759   fLayer  = fIndexes->GetLayer();
760   Int_t isector = fIndexes->GetSM();
761
762   // Start clustering in the chamber
763
764   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
765   if (fDet != det) {
766     AliError("Strange Detector number Missmatch!");
767     return kFALSE;
768   }
769
770   // TRD space point transformation
771   fTransform->SetDetector(det);
772
773   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
774   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
775   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule); 
776
777   if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
778     AddTrackletsToArray();
779
780   fColMax    = fDigits->GetNcol();
781   //Int_t nRowMax    = fDigits->GetNrow();
782   fTimeTotal = fDigits->GetNtime();
783
784   // Detector wise calibration object for the gain factors
785   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
786   // Calibration object with pad wise values for the gain factors
787   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
788   // Calibration value for chamber wise gain factor
789   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
790
791   // Detector wise calibration object for the noise
792   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
793   // Calibration object with pad wise values for the noise
794   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
795   // Calibration value for chamber wise noise
796   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
797   
798   // Calibration object with the pad status
799   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
800
801   SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
802   SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
803   SetBit(kHLT, fReconstructor->IsHLT());
804   
805   firstClusterROC = -1;
806   fClusterROC     =  0;
807
808   // Apply the gain and the tail cancelation via digital filter
809   if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
810
811   MaxStruct curr, last;
812   Int_t nMaximas = 0, nCorrupted = 0;
813
814   // Here the clusterfining is happening
815   
816   for(curr.Time = 0; curr.Time < fTimeTotal; curr.Time++){
817     while(fIndexes->NextRCIndex(curr.Row, curr.Col)){
818       //printf("\nCHECK r[%2d] c[%3d] t[%d]\n", curr.Row, curr.Col, curr.Time);
819       if(IsMaximum(curr, curr.padStatus, &curr.Signals[0])){
820         //printf("\tMAX s[%d %d %d]\n", curr.Signals[0], curr.Signals[1], curr.Signals[2]);
821         if(last.Row>-1){
822           if(curr.Time==last.Time && curr.Row==last.Row && curr.Col==last.Col+2) FivePadCluster(last, curr);
823           CreateCluster(last);
824         }
825         last=curr; curr.FivePad=kFALSE;
826       }
827       //printf("\t--- s[%d %d %d]\n", curr.Signals[0], curr.Signals[1], curr.Signals[2]);
828     }
829   }
830   if(last.Row>-1) CreateCluster(last);
831
832   if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
833     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
834     (*fDebugStream) << "MakeClusters"
835       << "Detector="   << det
836       << "NMaxima="    << nMaximas
837       << "NClusters="  << fClusterROC
838       << "NCorrupted=" << nCorrupted
839       << "\n";
840   }
841   if (TestBit(kLabels)) AddLabels();
842
843   return kTRUE;
844
845 }
846
847 //_____________________________________________________________________________
848 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals) 
849 {
850   //
851   // Returns true if this row,col,time combination is a maximum. 
852   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
853   //
854
855   Signals[1] = fDigits->GetData(Max.Row, Max.Col, Max.Time);
856   if(Signals[1] < fMaxThresh) return kFALSE;
857
858   Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.Col, Max.Row);
859   if (Signals[1] < noiseMiddleThresh) return kFALSE;
860
861   if (Max.Col + 1 >= fColMax || Max.Col < 1) return kFALSE;
862
863   UChar_t status[3]={
864     fCalPadStatusROC->GetStatus(Max.Col-1, Max.Row)
865    ,fCalPadStatusROC->GetStatus(Max.Col,   Max.Row)
866    ,fCalPadStatusROC->GetStatus(Max.Col+1, Max.Row)
867   };
868
869   Signals[0] = fDigits->GetData(Max.Row, Max.Col-1, Max.Time);
870   Signals[2] = fDigits->GetData(Max.Row, Max.Col+1, Max.Time);  
871
872   if(!(status[0] | status[1] | status[2])) {//all pads are good
873     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
874       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
875         Float_t  noiseSumThresh = fMinLeftRightCutSigma
876           * fCalNoiseDetValue
877           * fCalNoiseROC->GetValue(Max.Col, Max.Row);
878         if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
879         padStatus = 0;
880         return kTRUE;
881       }
882     }
883   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
884     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) { 
885       Signals[2]=0;
886       SetPadStatus(status[2], padStatus);
887       return kTRUE;
888     } 
889     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
890       Signals[0]=0;
891       SetPadStatus(status[0], padStatus);
892       return kTRUE;
893     }
894     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
895       Signals[1]=TMath::Nint(fMaxThresh);
896       SetPadStatus(status[1], padStatus);
897       return kTRUE;
898     }
899   }
900   return kFALSE;
901 }
902
903 //_____________________________________________________________________________
904 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
905 {
906   //
907   // Look for 5 pad cluster with minimum in the middle
908   // Gives back the ratio
909   //
910   if (ThisMax.Col >= fColMax - 3) return kFALSE;
911   if (ThisMax.Col < fColMax - 5){
912     if (fDigits->GetData(ThisMax.Row, ThisMax.Col+4, ThisMax.Time) >= fSigThresh)
913       return kFALSE;
914   }
915   if (ThisMax.Col > 1) {
916     if (fDigits->GetData(ThisMax.Row, ThisMax.Col-2, ThisMax.Time) >= fSigThresh)
917       return kFALSE;
918   }
919   
920   const Float_t kEpsilon = 0.01;
921   Double_t padSignal[5] = {ThisMax.Signals[0], ThisMax.Signals[1], ThisMax.Signals[2],
922       NeighbourMax.Signals[1], NeighbourMax.Signals[2]};
923   
924   // Unfold the two maxima and set the signal on 
925   // the overlapping pad to the ratio
926   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
927   ThisMax.Signals[2] = TMath::Nint(ThisMax.Signals[2]*ratio);
928   NeighbourMax.Signals[0] = TMath::Nint(NeighbourMax.Signals[0]*(1-ratio));
929   ThisMax.FivePad=kTRUE;
930   NeighbourMax.FivePad=kTRUE;
931   return kTRUE;
932
933 }
934
935 //_____________________________________________________________________________
936 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
937 {
938   //
939   // Creates a cluster at the given position and saves it in fRecPoints
940   //
941
942   Int_t nPadCount = 1;
943   Short_t signals[7] = { 0, 0, Max.Signals[0], Max.Signals[1], Max.Signals[2], 0, 0 };
944   if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
945
946   AliTRDcluster cluster(fDet, ((UChar_t) Max.Col), ((UChar_t) Max.Row), ((UChar_t) Max.Time), signals, fVolid);
947   cluster.SetNPads(nPadCount);
948   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
949   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
950   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
951
952   cluster.SetFivePad(Max.FivePad);
953   // set pads status for the cluster
954   UChar_t maskPosition = GetCorruption(Max.padStatus);
955   if (maskPosition) { 
956     cluster.SetPadMaskedPosition(maskPosition);
957     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
958   }
959
960   // Transform the local cluster coordinates into calibrated 
961   // space point positions defined in the local tracking system.
962   // Here the calibration for T0, Vdrift and ExB is applied as well.
963   if(!fTransform->Transform(&cluster)) return;
964   // Temporarily store the Max.Row, column and time bin of the center pad
965   // Used to later on assign the track indices
966   cluster.SetLabel(Max.Row, 0);
967   cluster.SetLabel(Max.Col, 1);
968   cluster.SetLabel(Max.Time,2);
969
970   //needed for HLT reconstruction
971   AddClusterToArray(&cluster); 
972
973   // Store the index of the first cluster in the current ROC
974   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
975   
976   fNoOfClusters++;
977   fClusterROC++;
978 }
979
980 //_____________________________________________________________________________
981 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
982 {
983   // Look to the right
984   Int_t ii = 1;
985   while (fDigits->GetData(Max.Row, Max.Col-ii, Max.Time) >= fSigThresh) {
986     nPadCount++;
987     ii++;
988     if (Max.Col < ii) break;
989   }
990   // Look to the left
991   ii = 1;
992   while (fDigits->GetData(Max.Row, Max.Col+ii, Max.Time) >= fSigThresh) {
993     nPadCount++;
994     ii++;
995     if (Max.Col+ii >= fColMax) break;
996   }
997
998   // Store the amplitudes of the pads in the cluster for later analysis
999   // and check whether one of these pads is masked in the database
1000   signals[2]=Max.Signals[0];
1001   signals[3]=Max.Signals[1];
1002   signals[4]=Max.Signals[2];
1003   for(Int_t i = 0; i<2; i++)
1004     {
1005       if(Max.Col+i >= 3)
1006         signals[i] = fDigits->GetData(Max.Row, Max.Col-3+i, Max.Time);
1007       if(Max.Col+3-i < fColMax)
1008         signals[6-i] = fDigits->GetData(Max.Row, Max.Col+3-i, Max.Time);
1009     }
1010   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1011     if ((jPad >= 0) && (jPad < fColMax))
1012       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1013       }*/
1014 }
1015
1016 //_____________________________________________________________________________
1017 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1018 {
1019   //
1020   // Add a cluster to the array
1021   //
1022
1023   Int_t n = RecPoints()->GetEntriesFast();
1024   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1025   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1026 }
1027
1028 //_____________________________________________________________________________
1029 void AliTRDclusterizer::AddTrackletsToArray()
1030 {
1031   //
1032   // Add the online tracklets of this chamber to the array
1033   //
1034
1035   UInt_t* trackletword;
1036   for(Int_t side=0; side<2; side++)
1037     {
1038       Int_t trkl=0;
1039       trackletword=fTrackletContainer[side];
1040       while(trackletword[trkl]>0){
1041         Int_t n = TrackletsArray()->GetEntriesFast();
1042         AliTRDtrackletWord tmp(trackletword[trkl]);
1043         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1044         trkl++;
1045       }
1046     }
1047 }
1048
1049 //_____________________________________________________________________________
1050 Bool_t AliTRDclusterizer::AddLabels()
1051 {
1052   //
1053   // Add the track indices to the found clusters
1054   //
1055   
1056   const Int_t   kNclus  = 3;  
1057   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1058   const Int_t   kNtrack = kNdict * kNclus;
1059
1060   Int_t iClusterROC = 0;
1061
1062   Int_t row  = 0;
1063   Int_t col  = 0;
1064   Int_t time = 0;
1065   Int_t iPad = 0;
1066
1067   // Temporary array to collect the track indices
1068   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1069
1070   // Loop through the dictionary arrays one-by-one
1071   // to keep memory consumption low
1072   AliTRDarrayDictionary *tracksIn = 0;  //mod
1073   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1074
1075     // tracksIn should be expanded beforehand!
1076     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1077
1078     // Loop though the clusters found in this ROC
1079     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1080
1081       AliTRDcluster *cluster = (AliTRDcluster *)
1082                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1083       row  = cluster->GetLabel(0);
1084       col  = cluster->GetLabel(1);
1085       time = cluster->GetLabel(2);
1086
1087       for (iPad = 0; iPad < kNclus; iPad++) {
1088         Int_t iPadCol = col - 1 + iPad;
1089         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1090         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;     
1091       }
1092
1093     }
1094
1095   }
1096
1097   // Copy the track indices into the cluster
1098   // Loop though the clusters found in this ROC
1099   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1100
1101     AliTRDcluster *cluster = (AliTRDcluster *)
1102       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1103     cluster->SetLabel(-9999,0);
1104     cluster->SetLabel(-9999,1);
1105     cluster->SetLabel(-9999,2);
1106   
1107     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1108
1109   }
1110
1111   delete [] idxTracks;
1112
1113   return kTRUE;
1114
1115 }
1116
1117 //_____________________________________________________________________________
1118 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1119 {
1120   //
1121   // Method to unfold neighbouring maxima.
1122   // The charge ratio on the overlapping pad is calculated
1123   // until there is no more change within the range given by eps.
1124   // The resulting ratio is then returned to the calling method.
1125   //
1126
1127   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1128   if (!calibration) {
1129     AliError("No AliTRDcalibDB instance available\n");
1130     return kFALSE;  
1131   }
1132   
1133   Int_t   irc                = 0;
1134   Int_t   itStep             = 0;                 // Count iteration steps
1135
1136   Double_t ratio             = 0.5;               // Start value for ratio
1137   Double_t prevRatio         = 0.0;               // Store previous ratio
1138
1139   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1140   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1141   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1142
1143   // Start the iteration
1144   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1145
1146     itStep++;
1147     prevRatio = ratio;
1148
1149     // Cluster position according to charge ratio
1150     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
1151                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1152     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
1153                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1154
1155     // Set cluster charge ratio
1156     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1157     Double_t ampLeft  = padSignal[1] / newSignal[1];
1158     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1159     Double_t ampRight = padSignal[3] / newSignal[1];
1160
1161     // Apply pad response to parameters
1162     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1163     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1164
1165     // Calculate new overlapping ratio
1166     ratio = TMath::Min((Double_t) 1.0
1167                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1168
1169   }
1170
1171   return ratio;
1172
1173 }
1174
1175 //_____________________________________________________________________________
1176 void AliTRDclusterizer::TailCancelation()
1177 {
1178   //
1179   // Applies the tail cancelation and gain factors: 
1180   // Transform fDigits to fDigits
1181   //
1182
1183   Int_t iRow  = 0;
1184   Int_t iCol  = 0;
1185   Int_t iTime = 0;
1186
1187   Double_t *inADC = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1188   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1189
1190   fIndexes->ResetCounters();
1191   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1192   while(fIndexes->NextRCIndex(iRow, iCol))
1193     {
1194       Float_t  fCalGainFactorROCValue = fCalGainFactorROC->GetValue(iCol,iRow);
1195       Double_t gain                  = fCalGainFactorDetValue 
1196                                     * fCalGainFactorROCValue;
1197
1198       Bool_t corrupted = kFALSE;
1199       for (iTime = 0; iTime < fTimeTotal; iTime++) 
1200         {         
1201           // Apply gain gain factor
1202           inADC[iTime]   = fDigits->GetData(iRow,iCol,iTime);
1203           if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1204           inADC[iTime]  /= gain;
1205           outADC[iTime]  = inADC[iTime];
1206           if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming()){
1207             (*fDebugStream) << "TailCancellation"
1208                               << "col="  << iCol
1209                               << "row="  << iRow
1210                               << "time=" << iTime
1211                               << "inADC=" << inADC[iTime]
1212                               << "gain=" << gain
1213                               << "outADC=" << outADC[iTime]
1214                               << "corrupted=" << corrupted
1215                               << "\n";
1216           }
1217         }
1218       if (!corrupted)
1219         {
1220           // Apply the tail cancelation via the digital filter
1221           // (only for non-coorupted pads)
1222           DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1223         }
1224
1225       for(iTime = 0; iTime < fTimeTotal; iTime++)//while (fIndexes->NextTbinIndex(iTime))
1226         {
1227           // Store the amplitude of the digit if above threshold
1228           if (outADC[iTime] > 0)
1229             fDigits->SetData(iRow,iCol,iTime,TMath::Nint(outADC[iTime]));
1230           else
1231             fDigits->SetData(iRow,iCol,iTime,0);
1232         } // while itime
1233
1234     } // while irow icol
1235
1236   delete [] inADC;
1237   delete [] outADC;
1238
1239   return;
1240
1241 }
1242
1243 //_____________________________________________________________________________
1244 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1245                                   ,const Int_t n, const Int_t nexp) 
1246 {
1247   //
1248   // Tail cancellation by deconvolution for PASA v4 TRF
1249   //
1250
1251   Double_t rates[2];
1252   Double_t coefficients[2];
1253
1254   // Initialization (coefficient = alpha, rates = lambda)
1255   Double_t r1 = 1.0;
1256   Double_t r2 = 1.0;
1257   Double_t c1 = 0.5;
1258   Double_t c2 = 0.5;
1259
1260   if (nexp == 1) {   // 1 Exponentials
1261     r1 = 1.156;
1262     r2 = 0.130;
1263     c1 = 0.066;
1264     c2 = 0.000;
1265   }
1266   if (nexp == 2) {   // 2 Exponentials
1267     Double_t par[4];
1268     fReconstructor->GetRecoParam()->GetTCParams(par);
1269     r1 = par[0];//1.156;
1270     r2 = par[1];//0.130;
1271     c1 = par[2];//0.114;
1272     c2 = par[3];//0.624;
1273   }
1274
1275   coefficients[0] = c1;
1276   coefficients[1] = c2;
1277
1278   Double_t dt = 0.1;
1279
1280   rates[0] = TMath::Exp(-dt/(r1));
1281   rates[1] = TMath::Exp(-dt/(r2));
1282   
1283   Int_t i = 0;
1284   Int_t k = 0;
1285
1286   Double_t reminder[2];
1287   Double_t correction = 0.0;
1288   Double_t result     = 0.0;
1289
1290   // Attention: computation order is important
1291   for (k = 0; k < nexp; k++) {
1292     reminder[k] = 0.0;
1293   }
1294
1295   for (i = 0; i < n; i++) {
1296
1297     result    = (source[i] - correction);    // No rescaling
1298     target[i] = result;
1299
1300     for (k = 0; k < nexp; k++) {
1301       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1302     }
1303
1304     correction = 0.0;
1305     for (k = 0; k < nexp; k++) {
1306       correction += reminder[k];
1307     }
1308
1309   }
1310
1311 }
1312
1313 //_____________________________________________________________________________
1314 void AliTRDclusterizer::ResetRecPoints() 
1315 {
1316   //
1317   // Resets the list of rec points
1318   //
1319
1320   if (fRecPoints) {
1321     fRecPoints->Delete();
1322     delete fRecPoints;
1323   }
1324 }
1325
1326 //_____________________________________________________________________________
1327 TClonesArray *AliTRDclusterizer::RecPoints() 
1328 {
1329   //
1330   // Returns the list of rec points
1331   //
1332
1333   if (!fRecPoints) {
1334     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1335       // determine number of clusters which has to be allocated
1336       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1337
1338       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1339     }
1340     //SetClustersOwner(kTRUE);
1341     AliTRDReconstructor::SetClusters(0x0);
1342   }
1343   return fRecPoints;
1344
1345 }
1346
1347 //_____________________________________________________________________________
1348 TClonesArray *AliTRDclusterizer::TrackletsArray() 
1349 {
1350   //
1351   // Returns the list of rec points
1352   //
1353
1354   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1355     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1356     //SetClustersOwner(kTRUE);
1357     //AliTRDReconstructor::SetTracklets(0x0);
1358   }
1359   return fTracklets;
1360
1361 }
1362