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