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