First V0 MC Analysis from H.Ricaud
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizer.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /* $Id$ */
18
19 ///////////////////////////////////////////////////////////////////////////////
20 //                                                                           //
21 // TRD cluster finder                                                        //
22 //                                                                           //
23 ///////////////////////////////////////////////////////////////////////////////
24
25 #include <TF1.h>
26 #include <TTree.h>
27 #include <TH1.h>
28 #include <TFile.h>
29 #include <TObjArray.h>
30
31 #include "AliRunLoader.h"
32 #include "AliLoader.h"
33 #include "AliRawReader.h"
34 #include "AliLog.h"
35 #include "AliAlignObj.h"
36
37 #include "AliTRDclusterizer.h"
38 #include "AliTRDcluster.h"
39 #include "AliTRDgeometry.h"
40 #include "AliTRDdataArrayF.h"
41 #include "AliTRDdataArrayI.h"
42 #include "AliTRDdigitsManager.h"
43 #include "AliTRDrawData.h"
44 #include "AliTRDcalibDB.h"
45 #include "AliTRDSimParam.h"
46 #include "AliTRDRecParam.h"
47 #include "AliTRDCommonParam.h"
48 #include "AliTRDtransform.h"
49 #include "AliTRDSignalIndex.h"
50 #include "AliTRDRawStream.h"
51 #include "AliTRDRawStreamV2.h"
52 #include "AliTRDfeeParam.h"
53
54 #include "Cal/AliTRDCalROC.h"
55 #include "Cal/AliTRDCalDet.h"
56
57 ClassImp(AliTRDclusterizer)
58
59 //_____________________________________________________________________________
60 AliTRDclusterizer::AliTRDclusterizer()
61   :TNamed()
62   ,fRunLoader(NULL)
63   ,fClusterTree(NULL)
64   ,fRecPoints(NULL)
65   ,fDigitsManager(NULL)
66   ,fAddLabels(kTRUE)
67   ,fRawVersion(2)
68   ,fIndexesOut(NULL)
69   ,fIndexesMaxima(NULL)
70   ,fTransform(NULL)
71 {
72   //
73   // AliTRDclusterizer default constructor
74   //
75
76   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
77
78 }
79
80 //_____________________________________________________________________________
81 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title)
82   :TNamed(name,title)
83   ,fRunLoader(NULL)
84   ,fClusterTree(NULL)
85   ,fRecPoints(NULL)
86   ,fDigitsManager(new AliTRDdigitsManager())
87   ,fAddLabels(kTRUE)
88   ,fRawVersion(2)
89   ,fIndexesOut(NULL)
90   ,fIndexesMaxima(NULL)
91   ,fTransform(new AliTRDtransform(0))
92 {
93   //
94   // AliTRDclusterizer constructor
95   //
96
97   fDigitsManager->CreateArrays();
98
99   fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
100
101 }
102
103 //_____________________________________________________________________________
104 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
105   :TNamed(c)
106   ,fRunLoader(NULL)
107   ,fClusterTree(NULL)
108   ,fRecPoints(NULL)
109   ,fDigitsManager(NULL)
110   ,fAddLabels(kTRUE)
111   ,fRawVersion(2)
112   ,fIndexesOut(NULL)
113   ,fIndexesMaxima(NULL)
114   ,fTransform(NULL)
115 {
116   //
117   // AliTRDclusterizer copy constructor
118   //
119
120 }
121
122 //_____________________________________________________________________________
123 AliTRDclusterizer::~AliTRDclusterizer()
124 {
125   //
126   // AliTRDclusterizer destructor
127   //
128
129   if (fRecPoints) 
130     {
131       fRecPoints->Delete();
132       delete fRecPoints;
133     }
134
135   if (fDigitsManager) 
136     {
137       delete fDigitsManager;
138       fDigitsManager = NULL;
139     }
140
141   if (fIndexesOut)
142     {
143       delete fIndexesOut;
144       fIndexesOut    = NULL;
145     }
146
147   if (fIndexesMaxima)
148     {
149       delete fIndexesMaxima;
150       fIndexesMaxima = NULL;
151     }
152
153   if (fTransform)
154     {
155       delete fTransform;
156       fTransform     = NULL;
157     }
158
159 }
160
161 //_____________________________________________________________________________
162 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
163 {
164   //
165   // Assignment operator
166   //
167
168   if (this != &c) 
169     {
170       ((AliTRDclusterizer &) c).Copy(*this);
171     }
172   return *this;
173
174 }
175
176 //_____________________________________________________________________________
177 void AliTRDclusterizer::Copy(TObject &c) const
178 {
179   //
180   // Copy function
181   //
182
183   ((AliTRDclusterizer &) c).fClusterTree   = NULL;
184   ((AliTRDclusterizer &) c).fRecPoints     = NULL;  
185   ((AliTRDclusterizer &) c).fDigitsManager = NULL;
186   ((AliTRDclusterizer &) c).fAddLabels     = fAddLabels;
187   ((AliTRDclusterizer &) c).fRawVersion    = fRawVersion;
188   ((AliTRDclusterizer &) c).fIndexesOut    = NULL;
189   ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
190   ((AliTRDclusterizer &) c).fTransform     = NULL;
191
192 }
193
194 //_____________________________________________________________________________
195 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
196 {
197   //
198   // Opens the AliROOT file. Output and input are in the same file
199   //
200
201   TString evfoldname = AliConfig::GetDefaultEventFolderName();
202   fRunLoader         = AliRunLoader::GetRunLoader(evfoldname);
203
204   if (!fRunLoader) {
205     fRunLoader = AliRunLoader::Open(name);
206   }
207
208   if (!fRunLoader) {
209     AliError(Form("Can not open session for file %s.",name));
210     return kFALSE;
211   }
212
213   OpenInput(nEvent);
214   OpenOutput();
215
216   return kTRUE;
217
218 }
219
220 //_____________________________________________________________________________
221 Bool_t AliTRDclusterizer::OpenOutput()
222 {
223   //
224   // Open the output file
225   //
226
227   TObjArray *ioArray = 0;
228
229   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
230   loader->MakeTree("R");
231
232   fClusterTree = loader->TreeR();
233   fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
234
235   return kTRUE;
236
237 }
238
239 //_____________________________________________________________________________
240 Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
241 {
242   //
243   // Connect the output tree
244   //
245
246   TObjArray *ioArray = 0;
247
248   fClusterTree = clusterTree;
249   fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
250
251   return kTRUE;
252
253 }
254
255 //_____________________________________________________________________________
256 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
257 {
258   //
259   // Opens a ROOT-file with TRD-hits and reads in the digits-tree
260   //
261
262   // Import the Trees for the event nEvent in the file
263   fRunLoader->GetEvent(nEvent);
264   
265   return kTRUE;
266
267 }
268
269 //_____________________________________________________________________________
270 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
271 {
272   //
273   // Fills TRDcluster branch in the tree with the clusters 
274   // found in detector = det. For det=-1 writes the tree. 
275   //
276
277   if ((det <                      -1) || 
278       (det >= AliTRDgeometry::Ndet())) {
279     AliError(Form("Unexpected detector index %d.\n",det));
280     return kFALSE;
281   }
282  
283   TBranch *branch = fClusterTree->GetBranch("TRDcluster");
284   if (!branch) {
285     TObjArray *ioArray = 0;
286     branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
287   }
288
289   if ((det >=                      0) && 
290       (det <  AliTRDgeometry::Ndet())) {
291
292     Int_t nRecPoints = RecPoints()->GetEntriesFast();
293     TObjArray *detRecPoints = new TObjArray(400);
294
295     for (Int_t i = 0; i < nRecPoints; i++) {
296       AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
297       if (det == c->GetDetector()) {
298         detRecPoints->AddLast(c);
299       }
300       else {
301         AliError(Form("Attempt to write a cluster with unexpected detector index: got=%d expected=%d\n"
302                      ,c->GetDetector()
303                      ,det));
304       }
305     }
306
307     branch->SetAddress(&detRecPoints);
308     fClusterTree->Fill();
309
310     delete detRecPoints;
311     
312     return kTRUE;
313
314   }
315
316   if (det == -1) {
317
318     AliInfo(Form("Writing the cluster tree %s for event %d."
319                 ,fClusterTree->GetName(),fRunLoader->GetEventNumber()));
320
321     if (fRecPoints) {
322
323       branch->SetAddress(&fRecPoints);
324
325       AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
326       loader->WriteRecPoints("OVERWRITE");
327   
328     }
329     else {
330
331       AliError("Cluster tree does not exist. Cannot write clusters.\n");
332       return kFALSE;
333
334     }
335
336     return kTRUE;  
337
338   }
339
340   AliError(Form("Unexpected detector index %d.\n",det));
341  
342   return kFALSE;  
343   
344 }
345
346 //_____________________________________________________________________________
347 void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
348 {
349   // 
350   // Reset the helper indexes
351   //
352
353   if (fIndexesOut)
354     {
355       // carefull here - we assume that only row number may change - most probable
356       if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
357         fIndexesOut->ResetContent();
358       else
359         fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
360                                            , indexesIn->GetNcol()
361                                            , indexesIn->GetNtime());
362     }
363   else
364     {
365       fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
366                                         , indexesIn->GetNcol() 
367                                         , indexesIn->GetNtime());
368     }
369   
370   if (fIndexesMaxima)
371     {
372       // carefull here - we assume that only row number may change - most probable
373       if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow()) 
374         {
375           fIndexesMaxima->ResetContent();
376         }
377       else
378         {
379           fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
380                                                 , indexesIn->GetNcol()
381                                                 , indexesIn->GetNtime());
382         }
383     }
384   else
385     {
386       fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
387                                            , indexesIn->GetNcol()
388                                            , indexesIn->GetNtime());
389     }
390
391 }
392
393 //_____________________________________________________________________________
394 Bool_t AliTRDclusterizer::ReadDigits()
395 {
396   //
397   // Reads the digits arrays from the input aliroot file
398   //
399
400   if (!fRunLoader) {
401     AliError("No run loader available");
402     return kFALSE;
403   }
404
405   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
406   if (!loader->TreeD()) {
407     loader->LoadDigits();
408   }
409
410   // Read in the digit arrays
411   return (fDigitsManager->ReadDigits(loader->TreeD()));
412
413 }
414
415 //_____________________________________________________________________________
416 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
417 {
418   //
419   // Reads the digits arrays from the input tree
420   //
421
422   // Read in the digit arrays
423   return (fDigitsManager->ReadDigits(digitsTree));
424
425 }
426
427 //_____________________________________________________________________________
428 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
429 {
430   //
431   // Reads the digits arrays from the ddl file
432   //
433
434   AliTRDrawData raw;
435   fDigitsManager = raw.Raw2Digits(rawReader);
436
437   return kTRUE;
438
439 }
440
441 //_____________________________________________________________________________
442 Bool_t AliTRDclusterizer::MakeClusters()
443 {
444   //
445   // Creates clusters from digits
446   //
447
448   // Propagate info from the digits manager
449   if (fAddLabels == kTRUE)
450     {
451       fAddLabels = fDigitsManager->UsesDictionaries();
452     }
453
454   Bool_t fReturn = kTRUE;
455   for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
456     {
457
458       AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(i);      
459       // This is to take care of switched off super modules
460       if (digitsIn->GetNtime() == 0) 
461         {
462           continue;
463         }
464       digitsIn->Expand();
465       AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
466       if (indexes->IsAllocated() == kFALSE)
467         {
468           fDigitsManager->BuildIndexes(i);
469         }
470
471       Bool_t fR = kFALSE;
472       if (indexes->HasEntry())
473         {
474           if (fAddLabels)
475             {
476               for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++) 
477                 {
478                   AliTRDdataArrayI *tracksIn = 0;
479                   tracksIn = fDigitsManager->GetDictionary(i,iDict);
480                   tracksIn->Expand();
481                 }
482             }
483           fR = MakeClusters(i);
484           fReturn = fR && fReturn;
485         }
486
487       if (fR == kFALSE)
488         {
489           WriteClusters(i);
490           ResetRecPoints();
491         }
492
493       // No compress just remove
494       fDigitsManager->RemoveDigits(i);
495       fDigitsManager->RemoveDictionaries(i);      
496       fDigitsManager->ClearIndexes(i);
497
498     }
499
500   return fReturn;
501
502 }
503
504 //_____________________________________________________________________________
505 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
506 {
507   //
508   // Creates clusters from raw data
509   //
510
511   AliTRDdataArrayI *digits = 0;
512   AliTRDdataArrayI *track0 = 0;
513   AliTRDdataArrayI *track1 = 0;
514   AliTRDdataArrayI *track2 = 0; 
515
516   AliTRDSignalIndex *indexes = 0;
517
518   // Create the digits manager
519   if (!fDigitsManager)
520     {
521       fDigitsManager = new AliTRDdigitsManager();
522       fDigitsManager->CreateArrays();
523     }
524
525   AliTRDRawStreamV2 input(rawReader);
526   input.SetRawVersion( fRawVersion );
527   input.Init();
528
529   AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
530
531   // Loop through the digits
532   Int_t lastdet = -1;
533   Int_t det     =  0;
534   Int_t it      =  0;
535   while (input.Next()) 
536     {
537
538       det = input.GetDet();
539
540       if (det != lastdet) 
541         {
542         
543           if (lastdet != -1)
544             {
545               digits = fDigitsManager->GetDigits(lastdet);
546               Bool_t iclusterBranch = kFALSE;
547               if (indexes->HasEntry())
548                 iclusterBranch = MakeClusters(lastdet);
549               if (iclusterBranch == kFALSE)
550                 {
551                   WriteClusters(lastdet);
552                   ResetRecPoints();
553                 }
554             }
555
556           if (digits)
557             {
558               fDigitsManager->RemoveDigits(lastdet);
559               fDigitsManager->RemoveDictionaries(lastdet);
560               fDigitsManager->ClearIndexes(lastdet);
561             }
562
563           lastdet = det;
564
565           // Add a container for the digits of this detector
566           digits = fDigitsManager->GetDigits(det);
567           track0 = fDigitsManager->GetDictionary(det,0);
568           track1 = fDigitsManager->GetDictionary(det,1);
569           track2 = fDigitsManager->GetDictionary(det,2);
570
571           // Allocate memory space for the digits buffer
572           if (digits->GetNtime() == 0) 
573             {
574               //AliDebug(5, Form("Alloc digits for det %d", det));
575               digits->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
576               track0->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
577               track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
578               track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
579             }
580           
581           indexes = fDigitsManager->GetIndexes(det);
582           indexes->SetSM(input.GetSM());
583           indexes->SetStack(input.GetStack());
584           indexes->SetLayer(input.GetLayer());
585           indexes->SetDetNumber(det);
586           if (indexes->IsAllocated() == kFALSE)
587             {
588               indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins());
589             }
590
591         }
592       
593       for (it = 0; it < 3; it++)
594         {
595           if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() )
596             {
597               if (input.GetSignals()[it] > 0)
598                 {
599                   digits->SetDataUnchecked(input.GetRow(), input.GetCol(),
600                                            input.GetTimeBin() + it, input.GetSignals()[it]);
601
602                   indexes->AddIndexTBin(input.GetRow(), input.GetCol(),
603                                         input.GetTimeBin() + it);
604                   track0->SetDataUnchecked(input.GetRow(), input.GetCol(),
605                                            input.GetTimeBin() + it, 0);
606                   track1->SetDataUnchecked(input.GetRow(), input.GetCol(),
607                                            input.GetTimeBin() + it, 0);
608                   track2->SetDataUnchecked(input.GetRow(), input.GetCol(),
609                                            input.GetTimeBin() + it, 0);
610                 }
611             }
612         }
613
614     }
615
616   if (lastdet != -1)
617     {
618       Bool_t iclusterBranch = kFALSE;
619       if (indexes->HasEntry()) 
620         {
621           iclusterBranch = MakeClusters(lastdet);
622         }
623       if (iclusterBranch == kFALSE)
624         {
625           WriteClusters(lastdet);
626           ResetRecPoints();
627         }
628       //MakeClusters(lastdet);
629       if (digits)
630         {
631           fDigitsManager->RemoveDigits(lastdet);
632           fDigitsManager->RemoveDictionaries(lastdet);
633           fDigitsManager->ClearIndexes(lastdet);
634         }
635     }
636
637   delete fDigitsManager;
638   fDigitsManager = NULL;
639   return kTRUE;
640
641 }
642
643 //_____________________________________________________________________________
644 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
645 {
646   //
647   // Creates clusters from raw data
648   //
649
650   // Create the digits manager
651   if (!fDigitsManager)
652     {
653       fDigitsManager = new AliTRDdigitsManager();
654       fDigitsManager->CreateArrays();
655     }
656
657   fDigitsManager->SetUseDictionaries(fAddLabels);
658
659   AliTRDRawStreamV2 input(rawReader);
660   input.SetRawVersion( fRawVersion );
661   input.Init();
662
663   AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
664   
665   Int_t det    = 0;
666   while ((det = input.NextChamber(fDigitsManager)) >= 0)
667     {
668       Bool_t iclusterBranch = kFALSE;
669       if (fDigitsManager->GetIndexes(det)->HasEntry())
670         {
671           iclusterBranch = MakeClusters(det);
672         }
673       if (iclusterBranch == kFALSE)
674         {
675           WriteClusters(det);
676           ResetRecPoints();
677         }
678       fDigitsManager->RemoveDigits(det);
679       fDigitsManager->RemoveDictionaries(det);      
680       fDigitsManager->ClearIndexes(det);
681     }
682
683   delete fDigitsManager;
684   fDigitsManager = NULL;
685   return kTRUE;
686
687 }
688
689 //_____________________________________________________________________________
690 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
691 {
692   //
693   // Generates the cluster.
694   //
695
696   // Get the digits
697   //   digits should be expanded beforehand!
698   //   digitsIn->Expand();
699   AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(det);      
700
701   // This is to take care of switched off super modules
702   if (digitsIn->GetNtime() == 0) 
703     {
704       return kFALSE;
705     }
706
707   AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
708   if (indexesIn->IsAllocated() == kFALSE)
709     {
710       AliError("Indexes do not exist!");
711       return kFALSE;      
712     }
713     
714   AliTRDcalibDB  *calibration = AliTRDcalibDB::Instance();
715   if (!calibration) 
716     {
717       AliFatal("No AliTRDcalibDB instance available\n");
718       return kFALSE;  
719     }
720
721   AliTRDRecParam *recParam    = AliTRDRecParam::Instance();
722   if (!recParam) 
723     {
724       AliError("No AliTRDRecParam instance available\n");
725       return kFALSE;  
726     }
727
728   // ADC thresholds
729   // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
730   Float_t ADCthreshold   = 0; 
731
732   // Threshold value for the maximum
733   Float_t maxThresh      = recParam->GetClusMaxThresh();
734   // Threshold value for the digit signal
735   Float_t sigThresh      = recParam->GetClusSigThresh();
736
737   // Iteration limit for unfolding procedure
738   const Float_t kEpsilon = 0.01;             
739   const Int_t   kNclus   = 3;  
740   const Int_t   kNsig    = 5;
741
742   Int_t    iUnfold       = 0;  
743   Double_t ratioLeft     = 1.0;
744   Double_t ratioRight    = 1.0;
745
746   Double_t padSignal[kNsig];   
747   Double_t clusterSignal[kNclus];
748
749   Int_t icham = indexesIn->GetChamber();
750   Int_t iplan = indexesIn->GetPlane();
751   Int_t isect = indexesIn->GetSM();
752
753   // Start clustering in the chamber
754
755   Int_t idet  = AliTRDgeometry::GetDetector(iplan,icham,isect);
756   if (idet != det)
757     {
758       AliError("Strange Detector number Missmatch!");
759       return kFALSE;
760     }
761
762   // TRD space point transformation
763   fTransform->SetDetector(det);
764
765   Int_t    ilayer  = AliGeomManager::kTRD1 + iplan;
766   Int_t    imodule = icham + AliTRDgeometry::Ncham() * isect;
767   UShort_t volid   = AliGeomManager::LayerToVolUID(ilayer,imodule); 
768
769   Int_t nColMax    = digitsIn->GetNcol();
770   Int_t nTimeTotal = digitsIn->GetNtime();
771
772   // Detector wise calibration object for the gain factors
773   const AliTRDCalDet *calGainFactorDet      = calibration->GetGainFactorDet();
774   // Calibration object with pad wise values for the gain factors
775   AliTRDCalROC       *calGainFactorROC      = calibration->GetGainFactorROC(idet);
776   // Calibration value for chamber wise gain factor
777   Float_t             calGainFactorDetValue = calGainFactorDet->GetValue(idet);
778
779   Int_t nClusters = 0;
780
781   AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow()
782                                                     ,digitsIn->GetNcol()
783                                                     ,digitsIn->GetNtime()); 
784
785   ResetHelperIndexes(indexesIn);
786
787   // Apply the gain and the tail cancelation via digital filter
788   TailCancelation(digitsIn
789                  ,digitsOut  
790                  ,indexesIn
791                  ,fIndexesOut
792                  ,nTimeTotal
793                  ,ADCthreshold
794                  ,calGainFactorROC
795                  ,calGainFactorDetValue);       
796         
797   Int_t row  = 0;
798   Int_t col  = 0;
799   Int_t time = 0;
800   Int_t iPad = 0;
801     
802   fIndexesOut->ResetCounters();
803   while (fIndexesOut->NextRCTbinIndex(row, col, time))
804     {
805
806       Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
807             
808       // Look for the maximum
809       if (signalM >= maxThresh) 
810         {
811                 
812           if (col + 1 >= nColMax || col-1 < 0)
813             continue;
814
815           Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
816           Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
817
818           if ((TMath::Abs(signalL) <= signalM) && 
819               (TMath::Abs(signalR) <  signalM)) 
820             {
821               if ((TMath::Abs(signalL) >= sigThresh) ||
822                   (TMath::Abs(signalR) >= sigThresh)) 
823                 {
824                   // Maximum found, mark the position by a negative signal
825                   digitsOut->SetDataUnchecked(row,col,time,-signalM);
826                   fIndexesMaxima->AddIndexTBin(row,col,time);
827                 }
828             }
829
830         }
831
832     }
833                
834   // The index to the first cluster of a given ROC
835   Int_t firstClusterROC = -1;
836   // The number of cluster in a given ROC
837   Int_t nClusterROC     =  0;
838
839   // Now check the maxima and calculate the cluster position
840   fIndexesMaxima->ResetCounters();
841   while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) 
842     {
843
844       // Maximum found ?             
845       if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) 
846         {
847
848           for (iPad = 0; iPad < kNclus; iPad++) 
849             {
850               Int_t iPadCol = col - 1 + iPad;
851               clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
852             }
853
854           // Count the number of pads in the cluster
855           Int_t nPadCount = 0;
856           Int_t ii;
857           // Look to the left
858           ii = 0;
859           while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii  ,time)) >= sigThresh) 
860             {
861               nPadCount++;
862               ii++;
863               if (col-ii   <        0) break;
864             }
865           // Look to the right
866           ii = 0;
867           while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) 
868             {
869               nPadCount++;
870               ii++;
871               if (col+ii+1 >= nColMax) break;
872             }
873           nClusters++;
874
875           // Look for 5 pad cluster with minimum in the middle
876           Bool_t fivePadCluster = kFALSE;
877           if (col < (nColMax - 3)) 
878             {
879               if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) 
880                 {
881                   fivePadCluster = kTRUE;
882                 }
883               if ((fivePadCluster) && (col < (nColMax - 5))) 
884                 {
885                   if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) 
886                     {
887                       fivePadCluster = kFALSE;
888                     }
889                 }
890               if ((fivePadCluster) && (col >             1)) 
891                 {
892                   if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) 
893                     {
894                       fivePadCluster = kFALSE;
895                     }
896                 }
897             }
898
899           // 5 pad cluster
900           // Modify the signal of the overlapping pad for the left part 
901           // of the cluster which remains from a previous unfolding
902           if (iUnfold) 
903             {
904               clusterSignal[0] *= ratioLeft;
905               iUnfold = 0;
906             }
907
908           // Unfold the 5 pad cluster
909           if (fivePadCluster) 
910             {
911               for (iPad = 0; iPad < kNsig; iPad++) 
912                 {
913                   padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
914                                                                           ,col-1+iPad
915                                                                           ,time));
916                 }
917               // Unfold the two maxima and set the signal on 
918               // the overlapping pad to the ratio
919               ratioRight        = Unfold(kEpsilon,iplan,padSignal);
920               ratioLeft         = 1.0 - ratioRight; 
921               clusterSignal[2] *= ratioRight;
922               iUnfold = 1;
923             }
924
925           // The position of the cluster in COL direction relative to the center pad (pad units)
926           Double_t clusterPosCol = 0.0;
927           if (recParam->LUTOn()) 
928             {
929               // Calculate the position of the cluster by using the
930               // lookup table method
931               clusterPosCol = recParam->LUTposition(iplan
932                                                    ,clusterSignal[0]
933                                                    ,clusterSignal[1]
934                                                    ,clusterSignal[2]);
935             }
936           else 
937             {
938               // Calculate the position of the cluster by using the
939               // center of gravity method
940               for (Int_t i = 0; i < kNsig; i++) 
941                 {
942                   padSignal[i] = 0.0;
943                 }
944               padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col  ,time)); // Central pad
945               padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left    pad
946               padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right   pad
947               if ((col >           2) && 
948                   (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) 
949                 {
950                   padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
951                 }
952               if ((col < nColMax - 3) &&
953                   (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])) 
954                 {
955                   padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
956                 }  
957               clusterPosCol = GetCOG(padSignal);
958             }
959
960           // Store the amplitudes of the pads in the cluster for later analysis
961           Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
962           for (Int_t jPad = col-3; jPad <= col+3; jPad++) 
963             {
964               if ((jPad <          0) || 
965                   (jPad >= nColMax-1)) 
966                 {
967                   continue;
968                 }
969               signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
970             }
971
972           // Transform the local cluster coordinates into calibrated 
973           // space point positions defined in the local tracking system.
974           // Here the calibration for T0, Vdrift and ExB is applied as well.
975           Double_t clusterXYZ[6];
976           clusterXYZ[0] = clusterPosCol;
977           clusterXYZ[1] = clusterSignal[0];
978           clusterXYZ[2] = clusterSignal[1];
979           clusterXYZ[3] = clusterSignal[2];
980           clusterXYZ[4] = 0.0;
981           clusterXYZ[5] = 0.0;
982           Int_t    clusterRCT[3];
983           clusterRCT[0] = row;
984           clusterRCT[1] = col;
985           clusterRCT[2] = 0;
986
987           if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),0)) {
988
989             // Add the cluster to the output array
990             // The track indices will be stored later 
991             Float_t clusterPos[3];
992             clusterPos[0] = clusterXYZ[0];
993             clusterPos[1] = clusterXYZ[1];
994             clusterPos[2] = clusterXYZ[2];
995             Float_t clusterSig[2];
996             clusterSig[0] = clusterXYZ[4];
997             clusterSig[1] = clusterXYZ[5];
998             Double_t clusterCharge  = clusterXYZ[3];
999             Char_t   clusterTimeBin = ((Char_t) clusterRCT[2]);
1000             AliTRDcluster *cluster = new AliTRDcluster(idet
1001                                                       ,clusterCharge
1002                                                       ,clusterPos
1003                                                       ,clusterSig
1004                                                       ,0x0
1005                                                       ,((Char_t) nPadCount)
1006                                                       ,signals
1007                                                       ,((UChar_t) col)
1008                                                       ,((UChar_t) row)
1009                                                       ,((UChar_t) time)
1010                                                       ,clusterTimeBin
1011                                                       ,clusterPosCol
1012                                                       ,volid);
1013
1014             // Temporarily store the row, column and time bin of the center pad
1015             // Used to later on assign the track indices
1016             cluster->SetLabel( row,0);
1017             cluster->SetLabel( col,1);
1018             cluster->SetLabel(time,2);
1019
1020             RecPoints()->Add(cluster);
1021
1022             // Store the index of the first cluster in the current ROC
1023             if (firstClusterROC < 0) 
1024               {
1025                 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
1026               }
1027
1028             // Count the number of cluster in the current ROC
1029             nClusterROC++;
1030
1031           } // if: Transform ok ?
1032
1033         } // if: Maximum found ?
1034
1035     }
1036
1037   delete digitsOut;
1038
1039   if (fAddLabels) 
1040     {
1041       AddLabels(idet, firstClusterROC, nClusterROC);
1042     }
1043
1044   // Write the cluster and reset the array
1045   WriteClusters(idet);
1046   ResetRecPoints();
1047
1048   return kTRUE;
1049
1050 }
1051
1052 //_____________________________________________________________________________
1053 Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
1054 {
1055   //
1056   // Add the track indices to the found clusters
1057   //
1058   
1059   const Int_t   kNclus  = 3;  
1060   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1061   const Int_t   kNtrack = kNdict * kNclus;
1062
1063   Int_t iClusterROC = 0;
1064
1065   Int_t row  = 0;
1066   Int_t col  = 0;
1067   Int_t time = 0;
1068   Int_t iPad = 0;
1069
1070   // Temporary array to collect the track indices
1071   Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1072
1073   // Loop through the dictionary arrays one-by-one
1074   // to keep memory consumption low
1075   AliTRDdataArrayI *tracksIn = 0;
1076   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1077
1078     // tracksIn should be expanded beforehand!
1079     tracksIn = fDigitsManager->GetDictionary(idet,iDict);
1080
1081     // Loop though the clusters found in this ROC
1082     for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1083  
1084       AliTRDcluster *cluster = (AliTRDcluster *)
1085         RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1086       row  = cluster->GetLabel(0);
1087       col  = cluster->GetLabel(1);
1088       time = cluster->GetLabel(2);
1089
1090       for (iPad = 0; iPad < kNclus; iPad++) {
1091         Int_t iPadCol = col - 1 + iPad;
1092         Int_t index   = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1093         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;     
1094       }
1095
1096     }
1097
1098   }
1099
1100   // Copy the track indices into the cluster
1101   // Loop though the clusters found in this ROC
1102   for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1103  
1104     AliTRDcluster *cluster = (AliTRDcluster *)
1105       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1106     cluster->SetLabel(-9999,0);
1107     cluster->SetLabel(-9999,1);
1108     cluster->SetLabel(-9999,2);
1109   
1110     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1111
1112   }
1113
1114   delete [] idxTracks;
1115
1116   return kTRUE;
1117
1118 }
1119
1120 //_____________________________________________________________________________
1121 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5])
1122 {
1123   //
1124   // Get COG position
1125   // Used for clusters with more than 3 pads - where LUT not applicable
1126   //
1127
1128   Double_t sum = signal[0]
1129                + signal[1]
1130                + signal[2] 
1131                + signal[3]
1132                + signal[4];
1133
1134   Double_t res = (0.0 * (-signal[0] + signal[4])
1135                       + (-signal[1] + signal[3])) / sum;
1136
1137   return res;             
1138
1139 }
1140
1141 //_____________________________________________________________________________
1142 Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
1143 {
1144   //
1145   // Method to unfold neighbouring maxima.
1146   // The charge ratio on the overlapping pad is calculated
1147   // until there is no more change within the range given by eps.
1148   // The resulting ratio is then returned to the calling method.
1149   //
1150
1151   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1152   if (!calibration) {
1153     AliError("No AliTRDcalibDB instance available\n");
1154     return kFALSE;  
1155   }
1156   
1157   Int_t   irc                = 0;
1158   Int_t   itStep             = 0;                 // Count iteration steps
1159
1160   Double_t ratio             = 0.5;               // Start value for ratio
1161   Double_t prevRatio         = 0.0;               // Store previous ratio
1162
1163   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1164   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1165   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1166
1167   // Start the iteration
1168   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1169
1170     itStep++;
1171     prevRatio = ratio;
1172
1173     // Cluster position according to charge ratio
1174     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
1175                       / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1176     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
1177                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1178
1179     // Set cluster charge ratio
1180     irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
1181     Double_t ampLeft  = padSignal[1] / newSignal[1];
1182     irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
1183     Double_t ampRight = padSignal[3] / newSignal[1];
1184
1185     // Apply pad response to parameters
1186     irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
1187     irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
1188
1189     // Calculate new overlapping ratio
1190     ratio = TMath::Min((Double_t) 1.0
1191                       ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1192
1193   }
1194
1195   return ratio;
1196
1197 }
1198
1199 //_____________________________________________________________________________
1200 void AliTRDclusterizer::TailCancelation(AliTRDdataArrayI *digitsIn
1201                                       , AliTRDdataArrayF *digitsOut
1202                                       , AliTRDSignalIndex *indexesIn
1203                                       , AliTRDSignalIndex *indexesOut
1204                                       , Int_t nTimeTotal
1205                                       , Float_t ADCthreshold
1206                                       , AliTRDCalROC *calGainFactorROC
1207                                       , Float_t calGainFactorDetValue)
1208 {
1209   //
1210   // Applies the tail cancelation and gain factors: 
1211   // Transform digitsIn to digitsOut
1212   //
1213
1214   Int_t iRow  = 0;
1215   Int_t iCol  = 0;
1216   Int_t iTime = 0;
1217
1218   AliTRDRecParam *recParam = AliTRDRecParam::Instance();
1219   if (!recParam) 
1220     {
1221       AliError("No AliTRDRecParam instance available\n");
1222       return;
1223     }
1224
1225   Double_t *inADC  = new Double_t[nTimeTotal];  // ADC data before tail cancellation
1226   Double_t *outADC = new Double_t[nTimeTotal];  // ADC data after tail cancellation
1227   indexesIn->ResetCounters();
1228   while (indexesIn->NextRCIndex(iRow, iCol))
1229     {
1230       Float_t  calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1231       Double_t gain                  = calGainFactorDetValue 
1232                                      * calGainFactorROCValue;
1233
1234       for (iTime = 0; iTime < nTimeTotal; iTime++) 
1235         {         
1236           // Apply gain gain factor
1237           inADC[iTime]   = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1238           inADC[iTime]  /= gain;
1239           outADC[iTime]  = inADC[iTime];
1240         }
1241
1242       // Apply the tail cancelation via the digital filter
1243       if (recParam->TCOn()) 
1244         {
1245           DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp());
1246         }
1247
1248       indexesIn->ResetTbinCounter();
1249       while (indexesIn->NextTbinIndex(iTime))
1250         {
1251           // Store the amplitude of the digit if above threshold
1252           if (outADC[iTime] > ADCthreshold) 
1253             {
1254               digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1255               indexesOut->AddIndexTBin(iRow,iCol,iTime);
1256             }     
1257         } // while itime
1258
1259     } // while irow icol
1260   
1261   delete [] inADC;
1262   delete [] outADC;
1263
1264   return;
1265
1266 }
1267
1268 //_____________________________________________________________________________
1269 void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1270                                  , Int_t n, Int_t nexp) 
1271 {
1272   //
1273   // Tail cancellation by deconvolution for PASA v4 TRF
1274   //
1275
1276   Double_t rates[2];
1277   Double_t coefficients[2];
1278
1279   // Initialization (coefficient = alpha, rates = lambda)
1280   Double_t r1 = 1.0;
1281   Double_t r2 = 1.0;
1282   Double_t c1 = 0.5;
1283   Double_t c2 = 0.5;
1284
1285   if (nexp == 1) {   // 1 Exponentials
1286     r1 = 1.156;
1287     r2 = 0.130;
1288     c1 = 0.066;
1289     c2 = 0.000;
1290   }
1291   if (nexp == 2) {   // 2 Exponentials
1292     r1 = 1.156;
1293     r2 = 0.130;
1294     c1 = 0.114;
1295     c2 = 0.624;
1296   }
1297
1298   coefficients[0] = c1;
1299   coefficients[1] = c2;
1300
1301   Double_t dt = 0.1;
1302
1303   rates[0] = TMath::Exp(-dt/(r1));
1304   rates[1] = TMath::Exp(-dt/(r2));
1305   
1306   Int_t i = 0;
1307   Int_t k = 0;
1308
1309   Double_t reminder[2];
1310   Double_t correction = 0.0;
1311   Double_t result     = 0.0;
1312
1313   // Attention: computation order is important
1314   for (k = 0; k < nexp; k++) {
1315     reminder[k] = 0.0;
1316   }
1317
1318   for (i = 0; i < n; i++) {
1319
1320     result    = (source[i] - correction);    // No rescaling
1321     target[i] = result;
1322
1323     for (k = 0; k < nexp; k++) {
1324       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1325     }
1326
1327     correction = 0.0;
1328     for (k = 0; k < nexp; k++) {
1329       correction += reminder[k];
1330     }
1331
1332   }
1333
1334 }
1335
1336 //_____________________________________________________________________________
1337 void AliTRDclusterizer::ResetRecPoints() 
1338 {
1339   //
1340   // Resets the list of rec points
1341   //
1342
1343   if (fRecPoints) {
1344     fRecPoints->Delete();
1345   }
1346
1347 }
1348
1349 //_____________________________________________________________________________
1350 TObjArray *AliTRDclusterizer::RecPoints() 
1351 {
1352   //
1353   // Returns the list of rec points
1354   //
1355
1356   if (!fRecPoints) {
1357     fRecPoints = new TObjArray(400);
1358   }
1359  
1360   return fRecPoints;
1361
1362 }