TOF + macros to each detector folder
[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){ // A.B.
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       (recoParam->UseTailCancelation())) {
919     // save a copy of raw data
920     if(TestBit(kRawSignal)){
921       if(fDigitsRaw){
922         fDigitsRaw->~AliTRDarrayADC();
923         new(fDigitsRaw) AliTRDarrayADC(*fDigits);
924       } else fDigitsRaw = new AliTRDarrayADC(*fDigits);
925     }
926     TailCancelation(recoParam);
927   }
928
929   MaxStruct curr, last;
930   Int_t nMaximas = 0, nCorrupted = 0;
931
932   // Here the clusterfining is happening
933   
934   for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
935     while(fIndexes->NextRCIndex(curr.row, curr.col)){
936       if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
937         if(last.row>-1){
938           if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
939           CreateCluster(last);
940         }
941         last=curr; curr.fivePad=kFALSE;
942       }
943     }
944   }
945   if(last.row>-1) CreateCluster(last);
946
947   if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
948     TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
949     (*fDebugStream) << "MakeClusters"
950       << "Detector="   << det
951       << "NMaxima="    << nMaximas
952       << "NClusters="  << fClusterROC
953       << "NCorrupted=" << nCorrupted
954       << "\n";
955   }
956   if (TestBit(kLabels)) AddLabels();
957
958   return kTRUE;
959
960 }
961
962 //_____________________________________________________________________________
963 Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Float_t *const Signals)
964 {
965   //
966   // Returns true if this row,col,time combination is a maximum. 
967   // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
968   //
969
970   Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
971   Float_t onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col) : 1;
972   Signals[1] = (fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) /(onlcf * gain) + 0.5f;
973   if(Signals[1] <= fMaxThresh) return kFALSE;
974
975   if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
976
977   Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
978   if (Signals[1] <= noiseMiddleThresh) return kFALSE;  
979
980   Char_t status[3]={
981     fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
982    ,fCalPadStatusROC->GetStatus(Max.col,   Max.row)
983    ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
984   };
985
986   Short_t signal(0);
987   if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
988     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
989     onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col-1) : 1;
990     Signals[0] = (signal - fBaseline) /( onlcf * gain) + 0.5f;
991   } else Signals[0] = 0.;
992   if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
993     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
994     onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col+1) : 1;
995     Signals[2] = (signal - fBaseline) /( onlcf *  gain) + 0.5f;
996   } else Signals[2] = 0.;
997
998   if(!(status[0] | status[1] | status[2])) {//all pads are good
999     if ((Signals[2] <= Signals[1]) && (Signals[0] <  Signals[1])) {
1000       if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
1001         if(Signals[0]<0) Signals[0]=0.;
1002         if(Signals[2]<0) Signals[2]=0.;
1003         Float_t noiseSumThresh = fMinLeftRightCutSigma * fCalNoiseDetValue
1004                                * fCalNoiseROC->GetValue(Max.col, Max.row);
1005         if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
1006         padStatus = 0;
1007         return kTRUE;
1008       }
1009     }
1010   } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
1011     if(Signals[0]<0)Signals[0]=0;
1012     if(Signals[2]<0)Signals[2]=0;
1013     if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) { 
1014       Signals[2]=0;
1015       SetPadStatus(status[2], padStatus);
1016       return kTRUE;
1017     } 
1018     else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
1019       Signals[0]=0;
1020       SetPadStatus(status[0], padStatus);
1021       return kTRUE;
1022     }
1023     else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
1024       Signals[1] = fMaxThresh;
1025       SetPadStatus(status[1], padStatus);
1026       return kTRUE;
1027     }
1028   }
1029   return kFALSE;
1030 }
1031
1032 //_____________________________________________________________________________
1033 Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
1034 {
1035   //
1036   // Look for 5 pad cluster with minimum in the middle
1037   // Gives back the ratio
1038   //
1039   
1040   if (ThisMax.col >= fColMax - 3) return kFALSE;
1041   Float_t gain;
1042   Float_t onlcf;
1043   if (ThisMax.col < fColMax - 5){
1044     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
1045     onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(ThisMax.row,ThisMax.col+4) : 1;
1046     if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain * onlcf)
1047       return kFALSE;
1048   }
1049   if (ThisMax.col > 1) {
1050     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
1051     onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(ThisMax.row,ThisMax.col-2) : 1;
1052     if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain * onlcf)
1053       return kFALSE;
1054   }
1055   
1056   const Float_t kEpsilon = 0.01;
1057   Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
1058       NeighbourMax.signals[1], NeighbourMax.signals[2]};
1059   
1060   // Unfold the two maxima and set the signal on 
1061   // the overlapping pad to the ratio
1062   Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
1063   ThisMax.signals[2] = ThisMax.signals[2]*ratio + 0.5f;
1064   NeighbourMax.signals[0] = NeighbourMax.signals[0]*(1-ratio) + 0.5f;
1065   ThisMax.fivePad=kTRUE;
1066   NeighbourMax.fivePad=kTRUE;
1067   return kTRUE;
1068
1069 }
1070
1071 //_____________________________________________________________________________
1072 void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
1073 {
1074   //
1075   // Creates a cluster at the given position and saves it in RecPoints
1076   //
1077
1078   Int_t nPadCount = 1;
1079   Short_t signals[7] = { 0, 0, (Short_t)Max.signals[0], (Short_t)Max.signals[1], (Short_t)Max.signals[2], 0, 0 };
1080   if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
1081
1082   AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
1083   cluster.SetNPads(nPadCount);
1084   cluster.SetQ(Max.signals[0]+Max.signals[1]+Max.signals[2]);
1085   if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1086   else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1087   else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
1088
1089   cluster.SetFivePad(Max.fivePad);
1090   // set pads status for the cluster
1091   UChar_t maskPosition = GetCorruption(Max.padStatus);
1092   if (maskPosition) { 
1093     cluster.SetPadMaskedPosition(maskPosition);
1094     cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1095   }
1096   cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
1097
1098   // Transform the local cluster coordinates into calibrated 
1099   // space point positions defined in the local tracking system.
1100   // Here the calibration for T0, Vdrift and ExB is applied as well.
1101   if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
1102   // Store raw signals in cluster. This MUST be called after position reconstruction !
1103   // Xianguo Lu and Alex Bercuci 19.03.2012
1104   if(TestBit(kRawSignal) && fDigitsRaw){
1105     Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
1106     Short_t rawSignal[7] = {0};
1107     for(Int_t ipad(Max.col-3), iRawId(0); ipad<=Max.col+3; ipad++, iRawId++){
1108       if(ipad<0 || ipad>=fColMax) continue;
1109       if(!fCalOnlGainROC){
1110         rawSignal[iRawId] = fDigitsRaw->GetData(Max.row, ipad, Max.time);
1111         continue;
1112       }
1113       // Deconvolute online gain calibration when available
1114       // Alex Bercuci 27.04.2012
1115       tmp = (fDigitsRaw->GetData(Max.row, ipad, Max.time) - fBaseline)/fCalOnlGainROC->GetGainCorrectionFactor(Max.row, ipad) + 0.5f;
1116       rawSignal[iRawId] = (Short_t)TMath::Min(tmp, kMaxShortVal);
1117     }
1118     cluster.SetSignals(rawSignal, kTRUE);
1119   }
1120   // Temporarily store the Max.Row, column and time bin of the center pad
1121   // Used to later on assign the track indices
1122   cluster.SetLabel(Max.row, 0);
1123   cluster.SetLabel(Max.col, 1);
1124   cluster.SetLabel(Max.time,2);
1125
1126   //needed for HLT reconstruction
1127   AddClusterToArray(&cluster);
1128
1129   // Store the index of the first cluster in the current ROC
1130   if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1131   
1132   fNoOfClusters++;
1133   fClusterROC++;
1134 }
1135
1136 //_____________________________________________________________________________
1137 void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1138 {
1139 // Calculate number of pads/cluster and
1140 // ADC signals at position 0, 1, 5 and 6
1141
1142   Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
1143   Float_t gain(1.); Float_t onlcf(1.); Short_t signal(0);
1144   // Store the amplitudes of the pads in the cluster for later analysis
1145   // and check whether one of these pads is masked in the database
1146   signals[3]=Max.signals[1];
1147   Int_t ipad(1), jpad(0);
1148   // Look to the right
1149   while((jpad = Max.col-ipad)){
1150     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1151     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1152     onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,jpad) : 1;
1153     tmp = (signal - fBaseline) / (onlcf * gain) + 0.5f;
1154     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1155     if(signal<fSigThresh) break; // signal under threshold
1156     nPadCount++;
1157     if(ipad<=3) signals[3 - ipad] = signal;
1158     ipad++;
1159   }
1160   ipad=1;
1161   // Look to the left
1162   while((jpad = Max.col+ipad)<fColMax){
1163     if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1164     gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1165     onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,jpad) : 1;
1166     tmp = (signal - fBaseline) / (onlcf * gain) + 0.5f;
1167     signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1168     if(signal<fSigThresh) break; // signal under threshold
1169     nPadCount++;
1170     if(ipad<=3) signals[3 + ipad] = signal;
1171     ipad++;
1172   }
1173
1174   AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1175     , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
1176 }
1177
1178 //_____________________________________________________________________________
1179 void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1180 {
1181   //
1182   // Add a cluster to the array
1183   //
1184
1185   Int_t n = RecPoints()->GetEntriesFast();
1186   if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1187   new((*RecPoints())[n]) AliTRDcluster(*cluster);
1188 }
1189
1190 //_____________________________________________________________________________
1191 Bool_t AliTRDclusterizer::AddLabels()
1192 {
1193   //
1194   // Add the track indices to the found clusters
1195   //
1196   
1197   const Int_t   kNclus  = 3;  
1198   const Int_t   kNdict  = AliTRDdigitsManager::kNDict;
1199   const Int_t   kNtrack = kNdict * kNclus;
1200
1201   Int_t iClusterROC = 0;
1202
1203   Int_t row  = 0;
1204   Int_t col  = 0;
1205   Int_t time = 0;
1206   Int_t iPad = 0;
1207
1208   // Temporary array to collect the track indices
1209   Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1210
1211   // Loop through the dictionary arrays one-by-one
1212   // to keep memory consumption low
1213   AliTRDarrayDictionary *tracksIn(NULL);  //mod
1214   for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1215
1216     // tracksIn should be expanded beforehand!
1217     tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1218
1219     // Loop though the clusters found in this ROC
1220     for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1221
1222       AliTRDcluster *cluster = (AliTRDcluster *)
1223                                RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1224       row  = cluster->GetLabel(0);
1225       col  = cluster->GetLabel(1);
1226       time = cluster->GetLabel(2);
1227
1228       for (iPad = 0; iPad < kNclus; iPad++) {
1229         Int_t iPadCol = col - 1 + iPad;
1230         Int_t index   = tracksIn->GetData(row,iPadCol,time);  //Modification of -1 in Track
1231         idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;     
1232       }
1233
1234     }
1235
1236   }
1237
1238   // Copy the track indices into the cluster
1239   // Loop though the clusters found in this ROC
1240   for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1241
1242     AliTRDcluster *cluster = (AliTRDcluster *)
1243       RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1244     cluster->SetLabel(-9999,0);
1245     cluster->SetLabel(-9999,1);
1246     cluster->SetLabel(-9999,2);
1247   
1248     cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1249
1250   }
1251
1252   delete [] idxTracks;
1253
1254   return kTRUE;
1255
1256 }
1257
1258 //_____________________________________________________________________________
1259 Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1260 {
1261   //
1262   // Method to unfold neighbouring maxima.
1263   // The charge ratio on the overlapping pad is calculated
1264   // until there is no more change within the range given by eps.
1265   // The resulting ratio is then returned to the calling method.
1266   //
1267
1268   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1269   if (!calibration) {
1270     AliError("No AliTRDcalibDB instance available\n");
1271     return kFALSE;  
1272   }
1273   
1274   Int_t   irc                = 0;
1275   Int_t   itStep             = 0;                 // Count iteration steps
1276
1277   Double_t ratio             = 0.5;               // Start value for ratio
1278   Double_t prevRatio         = 0.0;               // Store previous ratio
1279
1280   Double_t newLeftSignal[3]  = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1281   Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1282   Double_t newSignal[3]      = { 0.0, 0.0, 0.0 };
1283
1284   // Start the iteration
1285   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1286
1287     itStep++;
1288     prevRatio = ratio;
1289
1290     // Cluster position according to charge ratio
1291     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
1292                       / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1293     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
1294                       / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1295
1296     // Set cluster charge ratio
1297     irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1298     Double_t ampLeft  = padSignal[1] / newSignal[1];
1299     irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1300     Double_t ampRight = padSignal[3] / newSignal[1];
1301
1302     // Apply pad response to parameters
1303     irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1304     irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1305
1306     // Calculate new overlapping ratio
1307     // Coverity
1308     if (irc != 0) {
1309       ratio = TMath::Min((Double_t) 1.0
1310                         ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1311     }
1312
1313   }
1314
1315   return ratio;
1316
1317 }
1318
1319 //_____________________________________________________________________________
1320 void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1321 {
1322   //
1323   // Applies the tail cancelation
1324   //
1325
1326   Int_t nexp = recoParam->GetTCnexp();
1327   if(!nexp) return;
1328   
1329   Int_t iRow  = 0;
1330   Int_t iCol  = 0;
1331   Int_t iTime = 0;
1332
1333   TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1334   Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1335   while(fIndexes->NextRCIndex(iRow, iCol))
1336     {
1337       // if corrupted then don't make the tail cancallation
1338       if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1339
1340       if(debugStreaming){
1341         for (iTime = 0; iTime < fTimeTotal; iTime++) 
1342           (*fDebugStream) << "TailCancellation"
1343                           << "col="  << iCol
1344                           << "row="  << iRow
1345                           << "\n";
1346       }
1347       
1348       // Apply the tail cancelation via the digital filter
1349       //DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1350       ApplyTCTM(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1351     } // while irow icol
1352
1353   return;
1354
1355 }
1356
1357
1358 //_____________________________________________________________________________
1359 void AliTRDclusterizer::ApplyTCTM(Short_t *const arr, const Int_t nTime, const Int_t nexp) 
1360 {
1361   //
1362   // Steer tail cancellation
1363   //
1364
1365
1366   switch(nexp) {
1367   case 1:
1368   case 2:
1369     DeConvExp(arr,nTime,nexp);
1370     break;
1371   case -1:
1372     ConvExp(arr,nTime);
1373     break;
1374   case -2:
1375     DeConvExp(arr,nTime,1);
1376     ConvExp(arr,nTime);
1377     break;
1378   default:
1379     break;
1380   }
1381 }
1382
1383
1384 //_____________________________________________________________________________
1385 void AliTRDclusterizer::ConvExp(Short_t *const arr, const Int_t nTime)
1386 {
1387   //
1388   // Tail maker
1389   //
1390
1391   // Initialization (coefficient = alpha, rates = lambda)
1392   Float_t slope = 1.0;
1393   Float_t coeff = 0.5;
1394   Float_t rate;
1395
1396   Double_t par[4];
1397   fReconstructor->GetRecoParam()->GetTCParams(par);
1398   slope = par[1];
1399   coeff = par[3];  
1400
1401   Double_t dt = 0.1;
1402
1403   rate = TMath::Exp(-dt/(slope));
1404    
1405   Float_t reminder =  .0;
1406   Float_t correction = 0.0;
1407   Float_t result     = 0.0;
1408
1409   for (int i = nTime-1; i >= 0; i--) {
1410
1411     result = arr[i] + correction - fBaseline;    // No rescaling
1412     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1413
1414     correction = 0.0;
1415     
1416     correction += reminder = rate * (reminder + coeff * result);
1417   }
1418 }
1419
1420
1421 //_____________________________________________________________________________
1422 void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1423 {
1424   //
1425   // Tail cancellation by deconvolution for PASA v4 TRF
1426   //
1427
1428   Float_t rates[2];
1429   Float_t coefficients[2];
1430
1431   // Initialization (coefficient = alpha, rates = lambda)
1432   Float_t r1 = 1.0;
1433   Float_t r2 = 1.0;
1434   Float_t c1 = 0.5;
1435   Float_t c2 = 0.5;
1436
1437   if (nexp == 1) {   // 1 Exponentials
1438     r1 = 1.156;
1439     r2 = 0.130;
1440     c1 = 0.066;
1441     c2 = 0.000;
1442   }
1443   if (nexp == 2) {   // 2 Exponentials
1444     Double_t par[4];
1445     fReconstructor->GetRecoParam()->GetTCParams(par);
1446     r1 = par[0];//1.156;
1447     r2 = par[1];//0.130;
1448     c1 = par[2];//0.114;
1449     c2 = par[3];//0.624;
1450   }
1451
1452   coefficients[0] = c1;
1453   coefficients[1] = c2;
1454
1455   Double_t dt = 0.1;
1456
1457   rates[0] = TMath::Exp(-dt/(r1));
1458   rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
1459   
1460   Float_t reminder[2] = { .0, .0 };
1461   Float_t correction = 0.0;
1462   Float_t result     = 0.0;
1463
1464   for (int i = 0; i < nTime; i++) {
1465
1466     result = arr[i] - correction - fBaseline;    // No rescaling
1467     arr[i] = (Short_t)(result + fBaseline + 0.5f);
1468
1469     correction = 0.0;
1470     for (int k = 0; k < 2; k++) {
1471       correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1472     }
1473   }
1474 }
1475
1476 //_____________________________________________________________________________
1477 TClonesArray *AliTRDclusterizer::RecPoints() 
1478 {
1479   //
1480   // Returns the list of rec points
1481   //
1482
1483   fRecPoints = AliTRDReconstructor::GetClusters();
1484   if (!fRecPoints) AliError("Missing cluster array");
1485   return fRecPoints;
1486 }
1487
1488 //_____________________________________________________________________________
1489 TClonesArray *AliTRDclusterizer::TrackletsArray(const TString &trkltype)
1490 {
1491   //
1492   // Returns the array of on-line tracklets
1493   //
1494   fTracklets = AliTRDReconstructor::GetTracklets(trkltype.Data());
1495   if (!fTracklets) AliError("Missing online tracklets array");
1496   return fTracklets;
1497 }
1498
1499 //_____________________________________________________________________________
1500 TClonesArray* AliTRDclusterizer::TracksArray()
1501 {
1502   // return array of GTU tracks (create TClonesArray if necessary)
1503
1504   fTracks = AliTRDReconstructor::GetTracks();
1505   if (!fTracks) AliError("Missing online tracks array");
1506   return fTracks;
1507 }