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