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