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