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