fix stupid bug which prevented clusterizer to write clusters from last detector to
[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   
620   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
621
622   AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast())); 
623
624   return fReturn;
625
626 }
627
628 //_____________________________________________________________________________
629 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
630 {
631   //
632   // Creates clusters from raw data
633   //
634
635   return Raw2ClustersChamber(rawReader);
636
637 }
638
639 //_____________________________________________________________________________
640 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
641 {
642   //
643   // Creates clusters from raw data
644   //
645
646   // Create the digits manager
647   if (!fDigitsManager){
648     SetBit(knewDM, kTRUE);
649     fDigitsManager = new AliTRDdigitsManager(kTRUE);
650     fDigitsManager->CreateArrays();
651   }
652
653   fDigitsManager->SetUseDictionaries(TestBit(kLabels));
654
655   // tracklet container for raw tracklet writing
656   if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
657     // maximum tracklets for one HC
658     const Int_t kTrackletChmb=256;
659     fTrackletContainer = new UInt_t *[2];
660     fTrackletContainer[0] = new UInt_t[kTrackletChmb]; 
661     fTrackletContainer[1] = new UInt_t[kTrackletChmb]; 
662   }
663
664   if(!fRawStream)
665     fRawStream = AliTRDrawStreamBase::GetRawStream(rawReader);
666   else
667     fRawStream->SetReader(rawReader);
668
669   if(fReconstructor->IsHLT())
670     fRawStream->SetSharedPadReadout(kFALSE);
671
672   AliInfo(Form("Stream version: %s", fRawStream->IsA()->GetName()));
673   
674   Int_t det    = 0;
675   while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
676     Bool_t iclusterBranch = kFALSE;
677     if (fDigitsManager->GetIndexes(det)->HasEntry()){
678       iclusterBranch = MakeClusters(det);
679     }
680
681     fDigitsManager->ClearArrays(det);
682
683     if (!fReconstructor->IsWritingTracklets()) continue;
684     if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
685   }
686   
687   if (fTrackletContainer){
688     delete [] fTrackletContainer[0];
689     delete [] fTrackletContainer[1];
690     delete [] fTrackletContainer;
691     fTrackletContainer = NULL;
692   }
693
694   if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
695
696   if(!TestBit(knewDM)){
697     delete fDigitsManager;
698     fDigitsManager = NULL;
699   }
700
701   AliInfo(Form("Number of found clusters : %d", fNoOfClusters)); 
702   return kTRUE;
703
704 }
705
706 //_____________________________________________________________________________
707 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
708 {
709   //
710   // Check if a pad is masked
711   //
712
713   UChar_t status = 0;
714
715   if(signal>0 && TESTBIT(signal, 10)){
716     CLRBIT(signal, 10);
717     for(int ibit=0; ibit<4; ibit++){
718       if(TESTBIT(signal, 11+ibit)){
719         SETBIT(status, ibit);
720         CLRBIT(signal, 11+ibit);
721       } 
722     }
723   }
724   return status;
725 }
726
727 //_____________________________________________________________________________
728 void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
729   //
730   // Set the pad status into out
731   // First three bits are needed for the position encoding
732   //
733   out |= status << 3;
734 }
735
736 //_____________________________________________________________________________
737 UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
738   //
739   // return the staus encoding of the corrupted pad
740   //
741   return static_cast<UChar_t>(encoding >> 3);
742 }
743
744 //_____________________________________________________________________________
745 Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
746   //
747   // Return the position of the corruption
748   //
749   return encoding & 7;
750 }
751
752 //_____________________________________________________________________________
753 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
754 {
755   //
756   // Generates the cluster.
757   //
758
759   // Get the digits
760   fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
761   fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline();
762   
763   // This is to take care of switched off super modules
764   if (!fDigits->HasData()) return kFALSE;
765
766   fIndexes = fDigitsManager->GetIndexes(det);
767   if (fIndexes->IsAllocated() == kFALSE) {
768     AliError("Indexes do not exist!");
769     return kFALSE;      
770   }
771
772   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
773   if (!calibration) {
774     AliFatal("No AliTRDcalibDB instance available\n");
775     return kFALSE;  
776   }
777
778   if (!fReconstructor){
779     AliError("Reconstructor not set\n");
780     return kFALSE;
781   }
782
783
784   fMaxThresh            = fReconstructor->GetRecoParam()->GetClusMaxThresh();
785   fSigThresh            = fReconstructor->GetRecoParam()->GetClusSigThresh();
786   fMinMaxCutSigma       = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
787   fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
788
789   Int_t istack  = fIndexes->GetStack();
790   fLayer  = fIndexes->GetLayer();
791   Int_t isector = fIndexes->GetSM();
792
793   // Start clustering in the chamber
794
795   fDet  = AliTRDgeometry::GetDetector(fLayer,istack,isector);
796   if (fDet != det) {
797     AliError("Strange Detector number Missmatch!");
798     return kFALSE;
799   }
800
801   AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
802
803   // TRD space point transformation
804   fTransform->SetDetector(det);
805
806   Int_t    iGeoLayer  = AliGeomManager::kTRD1 + fLayer;
807   Int_t    iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
808   fVolid      = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule); 
809
810   if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
811     AddTrackletsToArray();
812
813   fColMax    = fDigits->GetNcol();
814   //Int_t nRowMax    = fDigits->GetNrow();
815   fTimeTotal = fDigits->GetNtime();
816
817   // Detector wise calibration object for the gain factors
818   const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
819   // Calibration object with pad wise values for the gain factors
820   fCalGainFactorROC      = calibration->GetGainFactorROC(fDet);
821   // Calibration value for chamber wise gain factor
822   fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
823
824   // Detector wise calibration object for the noise
825   const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
826   // Calibration object with pad wise values for the noise
827   fCalNoiseROC           = calibration->GetNoiseROC(fDet);
828   // Calibration value for chamber wise noise
829   fCalNoiseDetValue      = calNoiseDet->GetValue(fDet);
830   
831   // Calibration object with the pad status
832   fCalPadStatusROC       = calibration->GetPadStatusROC(fDet);
833
834   SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
835   SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
836   SetBit(kHLT, fReconstructor->IsHLT());
837   
838   firstClusterROC = -1;
839   fClusterROC     =  0;
840
841   // Apply the gain and the tail cancelation via digital filter
842   if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
843
844   MaxStruct curr, last;
845   Int_t nMaximas = 0, nCorrupted = 0;
846
847   // Here the clusterfining is happening
848   
849   for(curr.time = 0; curr.time < fTimeTotal; curr.time++){
850     while(fIndexes->NextRCIndex(curr.row, curr.col)){
851       if(IsMaximum(curr, curr.padStatus, &curr.signals[0])){
852         if(last.row>-1){
853           if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
854           CreateCluster(last);
855         }
856         last=curr; curr.fivePad=kFALSE;
857       }
858     }
859   }
860   if(last.row>-1) CreateCluster(last);
861
862   if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
863     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
864     (*fDebugStream) << "MakeClusters"
865       << "Detector="   << det
866       << "NMaxima="    << nMaximas
867       << "NClusters="  << fClusterROC
868       << "NCorrupted=" << nCorrupted
869       << "\n";
870   }
871   if (TestBit(kLabels)) AddLabels();
872
873   return kTRUE;
874
875 }
876
877 //_____________________________________________________________________________
878 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals) 
879 {
880   //
881   // Returns true if this row,col,time combination is a maximum. 
882   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
883   //
884
885   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
886   Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
887   if(Signals[1] < fMaxThresh) return kFALSE;
888
889   Float_t  noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
890   if (Signals[1] < noiseMiddleThresh) return kFALSE;
891
892   if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
893
894   UChar_t status[3]={
895     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
896    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
897    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
898   };
899
900   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
901   Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
902   gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
903   Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
904
905   if(!(status[0] | status[1] | status[2])) {//all pads are good
906     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
907       if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
908         if(Signals[0]<0)Signals[0]=0;
909         if(Signals[2]<0)Signals[2]=0;
910         Float_t  noiseSumThresh = fMinLeftRightCutSigma
911           * fCalNoiseDetValue
912           * fCalNoiseROC->GetValue(Max.col, Max.row);
913         if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
914         padStatus = 0;
915         return kTRUE;
916       }
917     }
918   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
919     if(Signals[0]<0)Signals[0]=0;
920     if(Signals[2]<0)Signals[2]=0;
921     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) { 
922       Signals[2]=0;
923       SetPadStatus(status[2], padStatus);
924       return kTRUE;
925     } 
926     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
927       Signals[0]=0;
928       SetPadStatus(status[0], padStatus);
929       return kTRUE;
930     }
931     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
932       Signals[1] = (Short_t)(fMaxThresh + 0.5f);
933       SetPadStatus(status[1], padStatus);
934       return kTRUE;
935     }
936   }
937   return kFALSE;
938 }
939
940 //_____________________________________________________________________________
941 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
942 {
943   //
944   // Look for 5 pad cluster with minimum in the middle
945   // Gives back the ratio
946   //
947   
948   if (ThisMax.col >= fColMax - 3) return kFALSE;
949   Float_t gain;
950   if (ThisMax.col < fColMax - 5){
951     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
952     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
953       return kFALSE;
954   }
955   if (ThisMax.col > 1) {
956     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
957     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
958       return kFALSE;
959   }
960   
961   const Float_t kEpsilon = 0.01;
962   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
963       NeighbourMax.signals[1], NeighbourMax.signals[2]};
964   
965   // Unfold the two maxima and set the signal on 
966   // the overlapping pad to the ratio
967   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
968   ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
969   NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
970   ThisMax.fivePad=kTRUE;
971   NeighbourMax.fivePad=kTRUE;
972   return kTRUE;
973
974 }
975
976 //_____________________________________________________________________________
977 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
978 {
979   //
980   // Creates a cluster at the given position and saves it in fRecPoints
981   //
982
983   Int_t nPadCount = 1;
984   Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
985   if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
986
987   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
988   cluster.SetNPads(nPadCount);
989   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
990   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
991   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
992
993   cluster.SetFivePad(Max.fivePad);
994   // set pads status for the cluster
995   UChar_t maskPosition = GetCorruption(Max.padStatus);
996   if (maskPosition) { 
997     cluster.SetPadMaskedPosition(maskPosition);
998     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
999   }
1000
1001   // Transform the local cluster coordinates into calibrated 
1002   // space point positions defined in the local tracking system.
1003   // Here the calibration for T0, Vdrift and ExB is applied as well.
1004   if(!fTransform->Transform(&cluster)) return;
1005   // Temporarily store the Max.Row, column and time bin of the center pad
1006   // Used to later on assign the track indices
1007   cluster.SetLabel(Max.row, 0);
1008   cluster.SetLabel(Max.col, 1);
1009   cluster.SetLabel(Max.time,2);
1010
1011   //needed for HLT reconstruction
1012   AddClusterToArray(&cluster); 
1013
1014   // Store the index of the first cluster in the current ROC
1015   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1016   
1017   fNoOfClusters++;
1018   fClusterROC++;
1019 }
1020
1021 //_____________________________________________________________________________
1022 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1023 {
1024   // Look to the right
1025   Int_t ii = 1;
1026   while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
1027     nPadCount++;
1028     ii++;
1029     if (Max.col < ii) break;
1030   }
1031   // Look to the left
1032   ii = 1;
1033   while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
1034     nPadCount++;
1035     ii++;
1036     if (Max.col+ii >= fColMax) break;
1037   }
1038
1039   // Store the amplitudes of the pads in the cluster for later analysis
1040   // and check whether one of these pads is masked in the database
1041   signals[2]=Max.signals[0];
1042   signals[3]=Max.signals[1];
1043   signals[4]=Max.signals[2];
1044   Float_t gain;
1045   for(Int_t i = 0; i<2; i++)
1046     {
1047       if(Max.col+i >= 3){
1048         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1049         signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1050       }
1051       if(Max.col+3-i < fColMax){
1052         gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1053         signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1054       }
1055     }
1056   /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1057     if ((jPad >= 0) && (jPad < fColMax))
1058       signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1059       }*/
1060 }
1061
1062 //_____________________________________________________________________________
1063 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1064 {
1065   //
1066   // Add a cluster to the array
1067   //
1068
1069   Int_t n = RecPoints()->GetEntriesFast();
1070   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1071   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1072 }
1073
1074 //_____________________________________________________________________________
1075 void AliTRDclusterizer::AddTrackletsToArray()
1076 {
1077   //
1078   // Add the online tracklets of this chamber to the array
1079   //
1080
1081   UInt_t* trackletword;
1082   for(Int_t side=0; side<2; side++)
1083     {
1084       Int_t trkl=0;
1085       trackletword=fTrackletContainer[side];
1086       while(trackletword[trkl]>0){
1087         Int_t n = TrackletsArray()->GetEntriesFast();
1088         AliTRDtrackletWord tmp(trackletword[trkl]);
1089         new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
1090         trkl++;
1091       }
1092     }
1093 }
1094
1095 //_____________________________________________________________________________
1096 Bool_t AliTRDclusterizer::AddLabels()
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*fClusterROC];
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(fDet,iDict);
1123
1124     // Loop though the clusters found in this ROC
1125     for (iClusterROC = 0; iClusterROC < fClusterROC; 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 < fClusterROC; 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 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1165 {
1166   //
1167   // Method to unfold neighbouring maxima.
1168   // The charge ratio on the overlapping pad is calculated
1169   // until there is no more change within the range given by eps.
1170   // The resulting ratio is then returned to the calling method.
1171   //
1172
1173   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1174   if (!calibration) {
1175     AliError("No AliTRDcalibDB instance available\n");
1176     return kFALSE;  
1177   }
1178   
1179   Int_t   irc                = 0;
1180   Int_t   itStep             = 0;                 // Count iteration steps
1181
1182   Double_t ratio             = 0.5;               // Start value for ratio
1183   Double_t prevRatio         = 0.0;               // Store previous ratio
1184
1185   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1186   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1187   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1188
1189   // Start the iteration
1190   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1191
1192     itStep++;
1193     prevRatio = ratio;
1194
1195     // Cluster position according to charge ratio
1196     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
1197                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1198     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
1199                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1200
1201     // Set cluster charge ratio
1202     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1203     Double_t ampLeft  = padSignal[1] / newSignal[1];
1204     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1205     Double_t ampRight = padSignal[3] / newSignal[1];
1206
1207     // Apply pad response to parameters
1208     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1209     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1210
1211     // Calculate new overlapping ratio
1212     ratio = TMath::Min((Double_t) 1.0
1213                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1214
1215   }
1216
1217   return ratio;
1218
1219 }
1220
1221 //_____________________________________________________________________________
1222 void AliTRDclusterizer::TailCancelation()
1223 {
1224   //
1225   // Applies the tail cancelation and gain factors: 
1226   // Transform fDigits to fDigits
1227   //
1228
1229   Int_t iRow  = 0;
1230   Int_t iCol  = 0;
1231   Int_t iTime = 0;
1232
1233   Double_t *inADC = new Double_t[fTimeTotal];  // ADC data before tail cancellation
1234   Double_t *outADC = new Double_t[fTimeTotal];  // ADC data after tail cancellation
1235
1236   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1237   Bool_t debugStreaming = fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1238   while(fIndexes->NextRCIndex(iRow, iCol))
1239     {
1240       Bool_t corrupted = kFALSE;
1241       if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1242
1243       // Save data into the temporary processing array and substract the baseline,
1244       // since DeConvExp does not expect a baseline
1245       for (iTime = 0; iTime < fTimeTotal; iTime++) 
1246         inADC[iTime]   = fDigits->GetData(iRow,iCol,iTime)-fBaseline;
1247           
1248       if(debugStreaming){
1249         for (iTime = 0; iTime < fTimeTotal; iTime++) 
1250           (*fDebugStream) << "TailCancellation"
1251                           << "col="  << iCol
1252                           << "row="  << iRow
1253                           << "time=" << iTime
1254                           << "inADC=" << inADC[iTime]
1255                           << "outADC=" << outADC[iTime]
1256                           << "corrupted=" << corrupted
1257                           << "\n";
1258       }
1259       
1260       if (!corrupted)
1261         {
1262           // Apply the tail cancelation via the digital filter
1263           // (only for non-coorupted pads)
1264           DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1265         }
1266       else memcpy(&outADC[0],&inADC[0],fTimeTotal*sizeof(inADC[0]));
1267
1268       // Save tailcancalled data and add the baseline
1269       for(iTime = 0; iTime < fTimeTotal; iTime++)
1270         fDigits->SetData(iRow,iCol,iTime,(Short_t)(outADC[iTime] + fBaseline + 0.5));
1271       
1272     } // while irow icol
1273
1274   delete [] inADC;
1275   delete [] outADC;
1276
1277   return;
1278
1279 }
1280
1281 //_____________________________________________________________________________
1282 void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1283                                   ,const Int_t n, const Int_t nexp) 
1284 {
1285   //
1286   // Tail cancellation by deconvolution for PASA v4 TRF
1287   //
1288
1289   Double_t rates[2];
1290   Double_t coefficients[2];
1291
1292   // Initialization (coefficient = alpha, rates = lambda)
1293   Double_t r1 = 1.0;
1294   Double_t r2 = 1.0;
1295   Double_t c1 = 0.5;
1296   Double_t c2 = 0.5;
1297
1298   if (nexp == 1) {   // 1 Exponentials
1299     r1 = 1.156;
1300     r2 = 0.130;
1301     c1 = 0.066;
1302     c2 = 0.000;
1303   }
1304   if (nexp == 2) {   // 2 Exponentials
1305     Double_t par[4];
1306     fReconstructor->GetRecoParam()->GetTCParams(par);
1307     r1 = par[0];//1.156;
1308     r2 = par[1];//0.130;
1309     c1 = par[2];//0.114;
1310     c2 = par[3];//0.624;
1311   }
1312
1313   coefficients[0] = c1;
1314   coefficients[1] = c2;
1315
1316   Double_t dt = 0.1;
1317
1318   rates[0] = TMath::Exp(-dt/(r1));
1319   rates[1] = TMath::Exp(-dt/(r2));
1320   
1321   Int_t i = 0;
1322   Int_t k = 0;
1323
1324   Double_t reminder[2];
1325   Double_t correction = 0.0;
1326   Double_t result     = 0.0;
1327
1328   // Attention: computation order is important
1329   for (k = 0; k < nexp; k++) {
1330     reminder[k] = 0.0;
1331   }
1332
1333   for (i = 0; i < n; i++) {
1334
1335     result    = (source[i] - correction);    // No rescaling
1336     target[i] = result;
1337
1338     for (k = 0; k < nexp; k++) {
1339       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1340     }
1341
1342     correction = 0.0;
1343     for (k = 0; k < nexp; k++) {
1344       correction += reminder[k];
1345     }
1346
1347   }
1348
1349 }
1350
1351 //_____________________________________________________________________________
1352 void AliTRDclusterizer::ResetRecPoints() 
1353 {
1354   //
1355   // Resets the list of rec points
1356   //
1357
1358   if (fRecPoints) {
1359     fRecPoints->Delete();
1360     delete fRecPoints;
1361   }
1362 }
1363
1364 //_____________________________________________________________________________
1365 TClonesArray *AliTRDclusterizer::RecPoints() 
1366 {
1367   //
1368   // Returns the list of rec points
1369   //
1370
1371   if (!fRecPoints) {
1372     if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1373       // determine number of clusters which has to be allocated
1374       Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1375
1376       fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1377     }
1378     //SetClustersOwner(kTRUE);
1379     AliTRDReconstructor::SetClusters(0x0);
1380   }
1381   return fRecPoints;
1382
1383 }
1384
1385 //_____________________________________________________________________________
1386 TClonesArray *AliTRDclusterizer::TrackletsArray() 
1387 {
1388   //
1389   // Returns the list of rec points
1390   //
1391
1392   if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
1393     fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
1394     //SetClustersOwner(kTRUE);
1395     //AliTRDReconstructor::SetTracklets(0x0);
1396   }
1397   return fTracklets;
1398
1399 }
1400