]> 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
fix "cw" option, small memory leak due to cluster array and
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizer.cxx
CommitLineData
f7336fa3 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
88cb7938 16/* $Id$ */
f7336fa3 17
18///////////////////////////////////////////////////////////////////////////////
19// //
3fe61b77 20// TRD cluster finder //
f7336fa3 21// //
22///////////////////////////////////////////////////////////////////////////////
23
3fe61b77 24#include <TF1.h>
94de3818 25#include <TTree.h>
3fe61b77 26#include <TH1.h>
793ff80c 27#include <TFile.h>
66f6bfd9 28#include <TClonesArray.h>
6d50f529 29#include <TObjArray.h>
f7336fa3 30
88cb7938 31#include "AliRunLoader.h"
32#include "AliLoader.h"
3fe61b77 33#include "AliRawReader.h"
6d50f529 34#include "AliLog.h"
3fe61b77 35#include "AliAlignObj.h"
88cb7938 36
f7336fa3 37#include "AliTRDclusterizer.h"
3e1a3ad8 38#include "AliTRDcluster.h"
fc546d21 39#include "AliTRDReconstructor.h"
793ff80c 40#include "AliTRDgeometry.h"
3fe61b77 41#include "AliTRDdataArrayF.h"
42#include "AliTRDdataArrayI.h"
625f5260 43#include "AliTRDdataArrayS.h"
fb44eb8e 44#include "AliTRDdataArrayDigits.h"
3fe61b77 45#include "AliTRDdigitsManager.h"
46#include "AliTRDrawData.h"
3551db50 47#include "AliTRDcalibDB.h"
fc546d21 48#include "AliTRDrecoParam.h"
b43a3e17 49#include "AliTRDCommonParam.h"
3fe61b77 50#include "AliTRDtransform.h"
51#include "AliTRDSignalIndex.h"
dfbb4bb9 52#include "AliTRDrawStreamBase.h"
3fe61b77 53#include "AliTRDfeeParam.h"
54
55#include "Cal/AliTRDCalROC.h"
56#include "Cal/AliTRDCalDet.h"
0e09df31 57#include "Cal/AliTRDCalSingleChamberStatus.h"
f7336fa3 58
59ClassImp(AliTRDclusterizer)
60
61//_____________________________________________________________________________
3a039a31 62AliTRDclusterizer::AliTRDclusterizer(AliTRDReconstructor *rec)
6d50f529 63 :TNamed()
3a039a31 64 ,fReconstructor(rec)
6d50f529 65 ,fRunLoader(NULL)
66 ,fClusterTree(NULL)
67 ,fRecPoints(NULL)
9c7c9ec1 68 ,fTrackletTree(NULL)
3fe61b77 69 ,fDigitsManager(NULL)
9c7c9ec1 70 ,fTrackletContainer(NULL)
3fe61b77 71 ,fAddLabels(kTRUE)
72 ,fRawVersion(2)
73 ,fIndexesOut(NULL)
74 ,fIndexesMaxima(NULL)
dfbb4bb9 75 ,fTransform(new AliTRDtransform(0))
fc546d21 76 ,fLUTbin(0)
77 ,fLUT(NULL)
f7336fa3 78{
79 //
80 // AliTRDclusterizer default constructor
81 //
82
a7ac01d2 83 AliTRDcalibDB *trd = 0x0;
84 if (!(trd = AliTRDcalibDB::Instance())) {
85 AliFatal("Could not get calibration object");
86 }
87
3fe61b77 88 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
89
3a039a31 90 // Initialize debug stream
91 if(fReconstructor){
92 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 1){
93 TDirectory *savedir = gDirectory;
94 //fgDebugStreamer = new TTreeSRedirector("TRD.ClusterizerDebug.root");
95 savedir->cd();
96 }
97 }
f7336fa3 98}
99
100//_____________________________________________________________________________
3a039a31 101AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, AliTRDReconstructor *rec)
6d50f529 102 :TNamed(name,title)
3a039a31 103 ,fReconstructor(rec)
6d50f529 104 ,fRunLoader(NULL)
105 ,fClusterTree(NULL)
106 ,fRecPoints(NULL)
9c7c9ec1 107 ,fTrackletTree(NULL)
3fe61b77 108 ,fDigitsManager(new AliTRDdigitsManager())
9c7c9ec1 109 ,fTrackletContainer(NULL)
3fe61b77 110 ,fAddLabels(kTRUE)
111 ,fRawVersion(2)
112 ,fIndexesOut(NULL)
113 ,fIndexesMaxima(NULL)
114 ,fTransform(new AliTRDtransform(0))
fc546d21 115 ,fLUTbin(0)
116 ,fLUT(NULL)
f7336fa3 117{
118 //
6d50f529 119 // AliTRDclusterizer constructor
f7336fa3 120 //
121
a7ac01d2 122 AliTRDcalibDB *trd = 0x0;
123 if (!(trd = AliTRDcalibDB::Instance())) {
124 AliFatal("Could not get calibration object");
125 }
126
3a039a31 127 // Initialize debug stream
128 if(fReconstructor){
129 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 1){
130 TDirectory *savedir = gDirectory;
131 //fgDebugStreamer = new TTreeSRedirector("TRD.ClusterizerDebug.root");
132 savedir->cd();
133 }
134 }
135
3fe61b77 136 fDigitsManager->CreateArrays();
137
138 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
139
fc546d21 140 FillLUT();
023b669c 141
f7336fa3 142}
143
8230f242 144//_____________________________________________________________________________
6d50f529 145AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
146 :TNamed(c)
3a039a31 147 ,fReconstructor(c.fReconstructor)
6d50f529 148 ,fRunLoader(NULL)
149 ,fClusterTree(NULL)
150 ,fRecPoints(NULL)
9c7c9ec1 151 ,fTrackletTree(NULL)
3fe61b77 152 ,fDigitsManager(NULL)
66f6bfd9 153 ,fTrackletContainer(NULL)
3fe61b77 154 ,fAddLabels(kTRUE)
155 ,fRawVersion(2)
156 ,fIndexesOut(NULL)
157 ,fIndexesMaxima(NULL)
158 ,fTransform(NULL)
fc546d21 159 ,fLUTbin(0)
160 ,fLUT(0)
8230f242 161{
162 //
163 // AliTRDclusterizer copy constructor
164 //
165
fc546d21 166 FillLUT();
167
8230f242 168}
169
f7336fa3 170//_____________________________________________________________________________
171AliTRDclusterizer::~AliTRDclusterizer()
172{
8230f242 173 //
174 // AliTRDclusterizer destructor
175 //
f7336fa3 176
48f8adf3 177 if (fRecPoints/* && IsClustersOwner()*/){
66f6bfd9 178 fRecPoints->Delete();
179 delete fRecPoints;
180 }
3fe61b77 181
66f6bfd9 182 if (fDigitsManager) {
183 delete fDigitsManager;
184 fDigitsManager = NULL;
185 }
9c7c9ec1 186
66f6bfd9 187 if (fTrackletContainer){
188 delete fTrackletContainer;
189 fTrackletContainer = NULL;
190 }
3fe61b77 191
66f6bfd9 192 if (fIndexesOut){
193 delete fIndexesOut;
194 fIndexesOut = NULL;
195 }
3fe61b77 196
66f6bfd9 197 if (fIndexesMaxima){
198 delete fIndexesMaxima;
199 fIndexesMaxima = NULL;
200 }
6d50f529 201
66f6bfd9 202 if (fTransform){
203 delete fTransform;
204 fTransform = NULL;
205 }
fc546d21 206
66f6bfd9 207 if (fLUT) {
208 delete [] fLUT;
209 fLUT = NULL;
210 }
f7336fa3 211}
212
8230f242 213//_____________________________________________________________________________
dd9a6ee3 214AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
215{
216 //
217 // Assignment operator
218 //
219
3fe61b77 220 if (this != &c)
221 {
222 ((AliTRDclusterizer &) c).Copy(*this);
223 }
053767a4 224
dd9a6ee3 225 return *this;
226
227}
228
229//_____________________________________________________________________________
e0d47c25 230void AliTRDclusterizer::Copy(TObject &c) const
8230f242 231{
232 //
233 // Copy function
234 //
235
3fe61b77 236 ((AliTRDclusterizer &) c).fClusterTree = NULL;
237 ((AliTRDclusterizer &) c).fRecPoints = NULL;
9c7c9ec1 238 ((AliTRDclusterizer &) c).fTrackletTree = NULL;
3fe61b77 239 ((AliTRDclusterizer &) c).fDigitsManager = NULL;
9c7c9ec1 240 ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
3fe61b77 241 ((AliTRDclusterizer &) c).fAddLabels = fAddLabels;
242 ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
243 ((AliTRDclusterizer &) c).fIndexesOut = NULL;
244 ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
245 ((AliTRDclusterizer &) c).fTransform = NULL;
fc546d21 246 ((AliTRDclusterizer &) c).fLUTbin = 0;
247 ((AliTRDclusterizer &) c).fLUT = NULL;
8230f242 248
249}
250
f7336fa3 251//_____________________________________________________________________________
3e1a3ad8 252Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
f7336fa3 253{
254 //
3e1a3ad8 255 // Opens the AliROOT file. Output and input are in the same file
f7336fa3 256 //
6d50f529 257
e191bb57 258 TString evfoldname = AliConfig::GetDefaultEventFolderName();
6d50f529 259 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
260
261 if (!fRunLoader) {
19dd5b2f 262 fRunLoader = AliRunLoader::Open(name);
6d50f529 263 }
264
265 if (!fRunLoader) {
266 AliError(Form("Can not open session for file %s.",name));
267 return kFALSE;
268 }
88cb7938 269
270 OpenInput(nEvent);
271 OpenOutput();
6d50f529 272
3e1a3ad8 273 return kTRUE;
f7336fa3 274
6d50f529 275}
3e1a3ad8 276
277//_____________________________________________________________________________
88cb7938 278Bool_t AliTRDclusterizer::OpenOutput()
3e1a3ad8 279{
280 //
281 // Open the output file
282 //
283
66f6bfd9 284 if (!fReconstructor->IsWritingClusters()) return kTRUE;
285
286 TObjArray *ioArray = 0x0;
88cb7938 287
288 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
289 loader->MakeTree("R");
6d50f529 290
88cb7938 291 fClusterTree = loader->TreeR();
66f6bfd9 292 fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
3e1a3ad8 293
3e1a3ad8 294 return kTRUE;
295
296}
297
25ca55ce 298//_____________________________________________________________________________
299Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
300{
301 //
302 // Connect the output tree
303 //
304
66f6bfd9 305 // clusters writing
306 if (fReconstructor->IsWritingClusters()){
307 TObjArray *ioArray = 0x0;
308 fClusterTree = clusterTree;
309 fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
310 }
9c7c9ec1 311
312 // tracklet writing
3a039a31 313 if (fReconstructor->IsWritingTracklets()){
9c7c9ec1 314 TString evfoldname = AliConfig::GetDefaultEventFolderName();
315 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
316
317 if (!fRunLoader) {
318 fRunLoader = AliRunLoader::Open("galice.root");
319 }
320 if (!fRunLoader) {
321 AliError(Form("Can not open session for file galice.root."));
322 return kFALSE;
323 }
324
325 UInt_t **leaves = new UInt_t *[2];
326 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
327 if (!dl) {
328 AliError("Could not get the tracklets data loader!");
329 AliDataLoader *dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
330 fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
331 }
332 else {
333 fTrackletTree = dl->Tree();
334 if (!fTrackletTree)
335 {
336 dl->MakeTree();
337 fTrackletTree = dl->Tree();
338 }
339 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
340 if (!trkbranch)
341 fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
342 }
343 }
344
25ca55ce 345 return kTRUE;
346
347}
348
3e1a3ad8 349//_____________________________________________________________________________
88cb7938 350Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
f7336fa3 351{
352 //
353 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
354 //
355
f7336fa3 356 // Import the Trees for the event nEvent in the file
bdbb05bb 357 fRunLoader->GetEvent(nEvent);
88cb7938 358
f7336fa3 359 return kTRUE;
360
361}
362
363//_____________________________________________________________________________
793ff80c 364Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
f7336fa3 365{
366 //
3e1a3ad8 367 // Fills TRDcluster branch in the tree with the clusters
793ff80c 368 // found in detector = det. For det=-1 writes the tree.
a3c76cdc 369 //
793ff80c 370
6d50f529 371 if ((det < -1) ||
372 (det >= AliTRDgeometry::Ndet())) {
373 AliError(Form("Unexpected detector index %d.\n",det));
3e1a3ad8 374 return kFALSE;
793ff80c 375 }
3e1a3ad8 376
66f6bfd9 377 TObjArray *ioArray = new TObjArray(400);
3e1a3ad8 378 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
379 if (!branch) {
3e1a3ad8 380 branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
66f6bfd9 381 } else branch->SetAddress(&ioArray);
3e1a3ad8 382
66f6bfd9 383
384 Int_t nRecPoints = RecPoints()->GetEntriesFast();
385 if(det >= 0){
3e1a3ad8 386 for (Int_t i = 0; i < nRecPoints; i++) {
bdbb05bb 387 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
66f6bfd9 388 if(det != c->GetDetector()) continue;
389 ioArray->AddLast(c);
793ff80c 390 }
3e1a3ad8 391 fClusterTree->Fill();
66f6bfd9 392 } else {
393 //AliInfo(Form("Writing the cluster tree %s for event %d.", fClusterTree->GetName(), fRunLoader->GetEventNumber()));
3e1a3ad8 394
66f6bfd9 395 Int_t detOld = -1;
396 for (Int_t i = 0; i < nRecPoints; i++) {
397 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
398 if(c->GetDetector() != detOld){
399 fClusterTree->Fill();
400 ioArray->Clear();
401 detOld = c->GetDetector();
402 }
403 ioArray->AddLast(c);
a6dd11e9 404 }
793ff80c 405 }
66f6bfd9 406 delete ioArray;
6d50f529 407
66f6bfd9 408 return kTRUE;
f7336fa3 409}
793ff80c 410
9c7c9ec1 411//_____________________________________________________________________________
412Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
413{
414
415 UInt_t **leaves = new UInt_t *[2];
416 for (Int_t i=0; i<2 ;i++){
417 leaves[i] = new UInt_t[258];
418 leaves[i][0] = det; // det
419 leaves[i][1] = i; // side
420 memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
421 }
422
423
424 if (!fTrackletTree){
425 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
426 dl->MakeTree();
427 fTrackletTree = dl->Tree();
428 }
429
430 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
431 if (!trkbranch) {
432 trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
433 }
434
435 for (Int_t i=0; i<2; i++){
436 if (leaves[i][2]>0) {
437 trkbranch->SetAddress(leaves[i]);
438 fTrackletTree->Fill();
439 }
440 }
441
442 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
443 dl->WriteData("OVERWRITE");
444 //dl->Unload();
445 delete [] leaves;
446
447 return kTRUE;
448}
449
3551db50 450//_____________________________________________________________________________
3fe61b77 451void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
452{
453 //
454 // Reset the helper indexes
455 //
456
457 if (fIndexesOut)
458 {
459 // carefull here - we assume that only row number may change - most probable
460 if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
461 fIndexesOut->ResetContent();
462 else
463 fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
464 , indexesIn->GetNcol()
465 , indexesIn->GetNtime());
466 }
467 else
468 {
469 fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
470 , indexesIn->GetNcol()
471 , indexesIn->GetNtime());
472 }
473
474 if (fIndexesMaxima)
475 {
476 // carefull here - we assume that only row number may change - most probable
477 if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
478 {
479 fIndexesMaxima->ResetContent();
480 }
481 else
482 {
483 fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
484 , indexesIn->GetNcol()
485 , indexesIn->GetNtime());
486 }
487 }
488 else
489 {
490 fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
491 , indexesIn->GetNcol()
492 , indexesIn->GetNtime());
493 }
494
495}
496
497//_____________________________________________________________________________
498Bool_t AliTRDclusterizer::ReadDigits()
499{
500 //
501 // Reads the digits arrays from the input aliroot file
502 //
503
504 if (!fRunLoader) {
505 AliError("No run loader available");
506 return kFALSE;
507 }
508
509 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
510 if (!loader->TreeD()) {
511 loader->LoadDigits();
512 }
513
514 // Read in the digit arrays
515 return (fDigitsManager->ReadDigits(loader->TreeD()));
516
517}
518
519//_____________________________________________________________________________
520Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
521{
522 //
523 // Reads the digits arrays from the input tree
524 //
525
526 // Read in the digit arrays
527 return (fDigitsManager->ReadDigits(digitsTree));
528
529}
530
531//_____________________________________________________________________________
532Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
533{
534 //
535 // Reads the digits arrays from the ddl file
536 //
537
538 AliTRDrawData raw;
539 fDigitsManager = raw.Raw2Digits(rawReader);
540
541 return kTRUE;
542
543}
544
545//_____________________________________________________________________________
546Bool_t AliTRDclusterizer::MakeClusters()
547{
548 //
549 // Creates clusters from digits
550 //
551
552 // Propagate info from the digits manager
66f6bfd9 553 if (fAddLabels == kTRUE){
554 fAddLabels = fDigitsManager->UsesDictionaries();
555 }
556
3fe61b77 557 Bool_t fReturn = kTRUE;
66f6bfd9 558 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
559
560 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits*) fDigitsManager->GetDigits(i);
561 // This is to take care of switched off super modules
562 if (!digitsIn->HasData()) continue;
563 digitsIn->Expand();
564 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
565 if (indexes->IsAllocated() == kFALSE){
566 fDigitsManager->BuildIndexes(i);
567 }
568
569 Bool_t fR = kFALSE;
570 if (indexes->HasEntry()){
571 if (fAddLabels){
572 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
573 AliTRDdataArrayI *tracksIn = 0;
574 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(i,iDict);
575 tracksIn->Expand();
3fe61b77 576 }
66f6bfd9 577 }
578 fR = MakeClusters(i);
579 fReturn = fR && fReturn;
3fe61b77 580 }
66f6bfd9 581
582 //if (fR == kFALSE){
583 // if(IsWritingClusters()) WriteClusters(i);
584 // ResetRecPoints();
585 //}
586
587 // No compress just remove
588 fDigitsManager->RemoveDigits(i);
589 fDigitsManager->RemoveDictionaries(i);
590 fDigitsManager->ClearIndexes(i);
591 }
592
593 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
3fe61b77 594
66f6bfd9 595 AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
3fe61b77 596
66f6bfd9 597 return fReturn;
3fe61b77 598}
599
600//_____________________________________________________________________________
601Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
602{
603 //
604 // Creates clusters from raw data
605 //
606
fc546d21 607 return Raw2ClustersChamber(rawReader);
3fe61b77 608
3fe61b77 609}
610
611//_____________________________________________________________________________
612Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
613{
614 //
615 // Creates clusters from raw data
616 //
617
618 // Create the digits manager
66f6bfd9 619 if (!fDigitsManager){
620 fDigitsManager = new AliTRDdigitsManager();
621 fDigitsManager->CreateArrays();
622 }
3fe61b77 623
624 fDigitsManager->SetUseDictionaries(fAddLabels);
625
dfbb4bb9 626 AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
627 AliTRDrawStreamBase &input = *pinput;
3fe61b77 628
629 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
630
631 Int_t det = 0;
66f6bfd9 632 while ((det = input.NextChamber(fDigitsManager)) >= 0){
633 Bool_t iclusterBranch = kFALSE;
634 if (fDigitsManager->GetIndexes(det)->HasEntry()){
635 iclusterBranch = MakeClusters(det);
3fe61b77 636 }
66f6bfd9 637
638 fDigitsManager->RemoveDigits(det);
639 fDigitsManager->RemoveDictionaries(det);
640 fDigitsManager->ClearIndexes(det);
9c7c9ec1 641 }
66f6bfd9 642 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
3fe61b77 643
644 delete fDigitsManager;
645 fDigitsManager = NULL;
dfbb4bb9 646
647 delete pinput;
648 pinput = NULL;
3fe61b77 649
66f6bfd9 650 AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
651 return kTRUE;
3fe61b77 652}
653
fa7427d0 654//_____________________________________________________________________________
655UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
656{
023b669c 657 //
658 // Check if a pad is masked
659 //
660
fa7427d0 661 UChar_t status = 0;
023b669c 662
fa7427d0 663 if(signal>0 && TESTBIT(signal, 10)){
664 CLRBIT(signal, 10);
665 for(int ibit=0; ibit<4; ibit++){
666 if(TESTBIT(signal, 11+ibit)){
667 SETBIT(status, ibit);
668 CLRBIT(signal, 11+ibit);
669 }
670 }
671 }
672 return status;
673}
674
3fe61b77 675//_____________________________________________________________________________
676Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
677{
678 //
679 // Generates the cluster.
680 //
681
682 // Get the digits
683 // digits should be expanded beforehand!
684 // digitsIn->Expand();
fb44eb8e 685 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits *) fDigitsManager->GetDigits(det);
f5375dcb 686
3fe61b77 687 // This is to take care of switched off super modules
625f5260 688 if (!digitsIn->HasData())
3fe61b77 689 {
690 return kFALSE;
691 }
692
693 AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
694 if (indexesIn->IsAllocated() == kFALSE)
695 {
696 AliError("Indexes do not exist!");
697 return kFALSE;
698 }
699
700 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
701 if (!calibration)
702 {
703 AliFatal("No AliTRDcalibDB instance available\n");
704 return kFALSE;
705 }
706
3fe61b77 707 // ADC thresholds
708 // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
acc49af9 709 Float_t adcThreshold = 0;
3fe61b77 710
3a039a31 711 if (!fReconstructor){
712 AliError("Reconstructor not set\n");
713 return kFALSE;
714 }
c5f589b9 715
3fe61b77 716 // Threshold value for the maximum
3a039a31 717 Float_t maxThresh = fReconstructor->GetRecoParam() ->GetClusMaxThresh();
3fe61b77 718 // Threshold value for the digit signal
3a039a31 719 Float_t sigThresh = fReconstructor->GetRecoParam() ->GetClusSigThresh();
3fe61b77 720
df83a620 721 // Threshold value for the maximum ( cut noise)
3a039a31 722 Float_t minMaxCutSigma = fReconstructor->GetRecoParam() ->GetMinMaxCutSigma();
df83a620 723 // Threshold value for the sum pad ( cut noise)
3a039a31 724 Float_t minLeftRightCutSigma = fReconstructor->GetRecoParam() ->GetMinLeftRightCutSigma();
df83a620 725
3fe61b77 726 // Iteration limit for unfolding procedure
727 const Float_t kEpsilon = 0.01;
728 const Int_t kNclus = 3;
729 const Int_t kNsig = 5;
730
731 Int_t iUnfold = 0;
732 Double_t ratioLeft = 1.0;
733 Double_t ratioRight = 1.0;
734
735 Double_t padSignal[kNsig];
736 Double_t clusterSignal[kNclus];
737
053767a4 738 Int_t istack = indexesIn->GetStack();
739 Int_t ilayer = indexesIn->GetLayer();
740 Int_t isector = indexesIn->GetSM();
3fe61b77 741
742 // Start clustering in the chamber
743
053767a4 744 Int_t idet = AliTRDgeometry::GetDetector(ilayer,istack,isector);
3fe61b77 745 if (idet != det)
746 {
747 AliError("Strange Detector number Missmatch!");
748 return kFALSE;
749 }
750
751 // TRD space point transformation
752 fTransform->SetDetector(det);
753
053767a4 754 Int_t iGeoLayer = AliGeomManager::kTRD1 + ilayer;
755 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
756 UShort_t volid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
3fe61b77 757
758 Int_t nColMax = digitsIn->GetNcol();
fa7427d0 759 Int_t nRowMax = digitsIn->GetNrow();
3fe61b77 760 Int_t nTimeTotal = digitsIn->GetNtime();
761
762 // Detector wise calibration object for the gain factors
0e09df31 763 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
3fe61b77 764 // Calibration object with pad wise values for the gain factors
0e09df31 765 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
3fe61b77 766 // Calibration value for chamber wise gain factor
0e09df31 767 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
768
df83a620 769
770 // Detector wise calibration object for the noise
771 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
772 // Calibration object with pad wise values for the noise
773 AliTRDCalROC *calNoiseROC = calibration->GetNoiseROC(idet);
774 // Calibration value for chamber wise noise
775 Float_t calNoiseDetValue = calNoiseDet->GetValue(idet);
776
3fe61b77 777 Int_t nClusters = 0;
778
fa7427d0 779 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(nRowMax, nColMax, nTimeTotal);
780 AliTRDdataArrayS padStatus(nRowMax, nColMax, nTimeTotal);
3fe61b77 781
782 ResetHelperIndexes(indexesIn);
783
784 // Apply the gain and the tail cancelation via digital filter
785 TailCancelation(digitsIn
023b669c 786 ,digitsOut
787 ,indexesIn
788 ,fIndexesOut
789 ,nTimeTotal
790 ,adcThreshold
791 ,calGainFactorROC
792 ,calGainFactorDetValue);
3fe61b77 793
794 Int_t row = 0;
795 Int_t col = 0;
796 Int_t time = 0;
797 Int_t iPad = 0;
798
f5375dcb 799 UChar_t status[3]={0, 0, 0}, ipos = 0;
3fe61b77 800 fIndexesOut->ResetCounters();
fa7427d0 801 while (fIndexesOut->NextRCTbinIndex(row, col, time)){
023b669c 802
f5375dcb 803 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
804 status[1] = digitsIn->GetPadStatus(row,col,time);
805 if(status[1]) SETBIT(ipos, AliTRDcluster::kMaskedCenter);
806
df83a620 807 if(signalM < maxThresh) continue;
808
809 Float_t noiseMiddleThresh = minMaxCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
810 if (signalM < noiseMiddleThresh) continue;
811
812 if (col + 1 >= nColMax || col-1 < 0) continue;
f5375dcb 813
df83a620 814 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
815 status[0] = digitsIn->GetPadStatus(row,col+1,time);
816 if(status[0]) SETBIT(ipos, AliTRDcluster::kMaskedLeft);
f5375dcb 817
df83a620 818 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
819 status[2] = digitsIn->GetPadStatus(row,col-1,time);
820 if(status[2]) SETBIT(ipos, AliTRDcluster::kMaskedRight);
f5375dcb 821
df83a620 822 // reject candidates with more than 1 problematic pad
823 if(ipos == 3 || ipos > 4) continue;
f5375dcb 824
df83a620 825 if(!status[1]){ // good central pad
826 if(!ipos){ // all pads are OK
827 if ((signalL <= signalM) && (signalR < signalM)) {
828 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
829 Float_t noiseSumThresh = minLeftRightCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
830 if((signalL+signalR+signalM) >= noiseSumThresh){
831 // Maximum found, mark the position by a negative signal
832 digitsOut->SetDataUnchecked(row,col,time,-signalM);
833 fIndexesMaxima->AddIndexTBin(row,col,time);
834 padStatus.SetDataUnchecked(row, col, time, ipos);
835 }
836 }
837 }
838 } else { // one of the neighbouring pads are bad
839 if(status[0] && signalR < signalM && signalR >= sigThresh){
840 digitsOut->SetDataUnchecked(row,col,time,-signalM);
841 digitsOut->SetDataUnchecked(row, col, time+1, 0.);
842 fIndexesMaxima->AddIndexTBin(row,col,time);
843 padStatus.SetDataUnchecked(row, col, time, ipos);
844 } else if(status[2] && signalL <= signalM && signalL >= sigThresh){
845 digitsOut->SetDataUnchecked(row,col,time,-signalM);
846 digitsOut->SetDataUnchecked(row, col, time-1, 0.);
023b669c 847 fIndexesMaxima->AddIndexTBin(row,col,time);
848 padStatus.SetDataUnchecked(row, col, time, ipos);
df83a620 849 }
850 }
851 } else { // wrong maximum pad
852 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
853 // Maximum found, mark the position by a negative signal
854 digitsOut->SetDataUnchecked(row,col,time,-maxThresh);
855 fIndexesMaxima->AddIndexTBin(row,col,time);
856 padStatus.SetDataUnchecked(row, col, time, ipos);
f5375dcb 857 }
858 }
df83a620 859 }
f5375dcb 860
861 // The index to the first cluster of a given ROC
df83a620 862 Int_t firstClusterROC = -1;
f5375dcb 863 // The number of cluster in a given ROC
864 Int_t nClusterROC = 0;
865
866 // Now check the maxima and calculate the cluster position
867 fIndexesMaxima->ResetCounters();
868 while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) {
3fe61b77 869
870 // Maximum found ?
f5375dcb 871 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
023b669c 872
f5375dcb 873 for (iPad = 0; iPad < kNclus; iPad++) {
874 Int_t iPadCol = col - 1 + iPad;
875 clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
876 }
3fe61b77 877
f5375dcb 878 // Count the number of pads in the cluster
879 Int_t nPadCount = 0;
880 Int_t ii;
023b669c 881 // Look to the right
f5375dcb 882 ii = 0;
883 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) {
884 nPadCount++;
885 ii++;
886 if (col-ii < 0) break;
887 }
023b669c 888 // Look to the left
f5375dcb 889 ii = 0;
890 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh){
891 nPadCount++;
892 ii++;
893 if (col+ii+1 >= nColMax) break;
894 }
895 nClusters++;
896
897 // Look for 5 pad cluster with minimum in the middle
898 Bool_t fivePadCluster = kFALSE;
899 if (col < (nColMax - 3)){
900 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
901 fivePadCluster = kTRUE;
902 }
903 if ((fivePadCluster) && (col < (nColMax - 5))) {
904 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh){
905 fivePadCluster = kFALSE;
906 }
907 }
908 if ((fivePadCluster) && (col > 1)){
909 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh){
910 fivePadCluster = kFALSE;
911 }
912 }
913 }
3fe61b77 914
f5375dcb 915 // 5 pad cluster
916 // Modify the signal of the overlapping pad for the left part
917 // of the cluster which remains from a previous unfolding
918 if (iUnfold) {
919 clusterSignal[0] *= ratioLeft;
920 iUnfold = 0;
921 }
3fe61b77 922
f5375dcb 923 // Unfold the 5 pad cluster
924 if (fivePadCluster){
925 for (iPad = 0; iPad < kNsig; iPad++) {
926 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
023b669c 927 ,col-1+iPad
928 ,time));
f5375dcb 929 }
930 // Unfold the two maxima and set the signal on
931 // the overlapping pad to the ratio
053767a4 932 ratioRight = Unfold(kEpsilon,ilayer,padSignal);
f5375dcb 933 ratioLeft = 1.0 - ratioRight;
934 clusterSignal[2] *= ratioRight;
935 iUnfold = 1;
936 }
3fe61b77 937
f5375dcb 938 // The position of the cluster in COL direction relative to the center pad (pad units)
939 Double_t clusterPosCol = 0.0;
3a039a31 940 if (fReconstructor->GetRecoParam() ->IsLUT()) {
f5375dcb 941 // Calculate the position of the cluster by using the
942 // lookup table method
023b669c 943 clusterPosCol = LUTposition(ilayer,clusterSignal[2]
944 ,clusterSignal[1]
945 ,clusterSignal[0]);
946 }
947 else {
f5375dcb 948 // Calculate the position of the cluster by using the
949 // center of gravity method
950 for (Int_t i = 0; i < kNsig; i++) {
951 padSignal[i] = 0.0;
952 }
953 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
023b669c 954 padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Left pad
955 padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Right pad
f5375dcb 956 if ((col > 2) &&
957 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
023b669c 958 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
f5375dcb 959 }
960 if ((col < nColMax - 3) &&
961 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])){
023b669c 962 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
f5375dcb 963 }
964 clusterPosCol = GetCOG(padSignal);
965 }
3fe61b77 966
f5375dcb 967 // Store the amplitudes of the pads in the cluster for later analysis
968 // and check whether one of these pads is masked in the database
969 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
970 for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
971 if ((jPad < 0) ||
972 (jPad >= nColMax-1)) {
973 continue;
974 }
975 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
976 }
3fe61b77 977
f5375dcb 978 // Transform the local cluster coordinates into calibrated
979 // space point positions defined in the local tracking system.
980 // Here the calibration for T0, Vdrift and ExB is applied as well.
981 Double_t clusterXYZ[6];
982 clusterXYZ[0] = clusterPosCol;
023b669c 983 clusterXYZ[1] = clusterSignal[2];
f5375dcb 984 clusterXYZ[2] = clusterSignal[1];
023b669c 985 clusterXYZ[3] = clusterSignal[0];
f5375dcb 986 clusterXYZ[4] = 0.0;
987 clusterXYZ[5] = 0.0;
988 Int_t clusterRCT[3];
989 clusterRCT[0] = row;
990 clusterRCT[1] = col;
991 clusterRCT[2] = 0;
992
993 Bool_t out = kTRUE;
994 if (fTransform->Transform(clusterXYZ, clusterRCT, ((UInt_t) time), out, 0)) {
995
996 // Add the cluster to the output array
997 // The track indices will be stored later
998 Float_t clusterPos[3];
999 clusterPos[0] = clusterXYZ[0];
1000 clusterPos[1] = clusterXYZ[1];
1001 clusterPos[2] = clusterXYZ[2];
1002 Float_t clusterSig[2];
1003 clusterSig[0] = clusterXYZ[4];
1004 clusterSig[1] = clusterXYZ[5];
1005 Double_t clusterCharge = clusterXYZ[3];
1006 Char_t clusterTimeBin = ((Char_t) clusterRCT[2]);
66f6bfd9 1007
1008 Int_t n = RecPoints()->GetEntriesFast();
1009 AliTRDcluster *cluster = new((*RecPoints())[n]) AliTRDcluster(idet
1010 ,clusterCharge
1011 ,clusterPos
1012 ,clusterSig
1013 ,0x0
1014 ,((Char_t) nPadCount)
1015 ,signals
1016 ,((UChar_t) col)
1017 ,((UChar_t) row)
1018 ,((UChar_t) time)
1019 ,clusterTimeBin
1020 ,clusterPosCol
1021 ,volid);
f5375dcb 1022 cluster->SetInChamber(!out);
023b669c 1023
f5375dcb 1024 UChar_t maskPosition = padStatus.GetDataUnchecked(row, col, time);
023b669c 1025 if (maskPosition) {
f5375dcb 1026 cluster->SetPadMaskedPosition(maskPosition);
023b669c 1027 if (maskPosition & AliTRDcluster::kMaskedLeft) {
1028 cluster->SetPadMaskedStatus(status[0]);
1029 }
1030 else if (maskPosition & AliTRDcluster::kMaskedCenter) {
1031 cluster->SetPadMaskedStatus(status[1]);
1032 }
1033 else {
1034 cluster->SetPadMaskedStatus(status[2]);
1035 }
f5375dcb 1036 }
3fe61b77 1037
f5375dcb 1038 // Temporarily store the row, column and time bin of the center pad
1039 // Used to later on assign the track indices
1040 cluster->SetLabel( row,0);
1041 cluster->SetLabel( col,1);
1042 cluster->SetLabel(time,2);
1043
66f6bfd9 1044 //RecPoints()->Add(cluster);
3fe61b77 1045
f5375dcb 1046 // Store the index of the first cluster in the current ROC
1047 if (firstClusterROC < 0) {
1048 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
1049 }
1050
1051 // Count the number of cluster in the current ROC
1052 nClusterROC++;
1053
1054 } // if: Transform ok ?
1055 } // if: Maximum found ?
1056 }
3fe61b77 1057
1058 delete digitsOut;
1059
f5375dcb 1060 if (fAddLabels) AddLabels(idet, firstClusterROC, nClusterROC);
3fe61b77 1061
1062 // Write the cluster and reset the array
66f6bfd9 1063 //WriteClusters(idet);
1064 //ResetRecPoints();
3fe61b77 1065
1066 return kTRUE;
1067
1068}
1069
1070//_____________________________________________________________________________
1071Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
3551db50 1072{
1073 //
3fe61b77 1074 // Add the track indices to the found clusters
3551db50 1075 //
1076
3fe61b77 1077 const Int_t kNclus = 3;
1078 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1079 const Int_t kNtrack = kNdict * kNclus;
1080
1081 Int_t iClusterROC = 0;
1082
1083 Int_t row = 0;
1084 Int_t col = 0;
1085 Int_t time = 0;
1086 Int_t iPad = 0;
1087
1088 // Temporary array to collect the track indices
1089 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1090
1091 // Loop through the dictionary arrays one-by-one
1092 // to keep memory consumption low
1093 AliTRDdataArrayI *tracksIn = 0;
1094 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1095
1096 // tracksIn should be expanded beforehand!
625f5260 1097 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
3fe61b77 1098
1099 // Loop though the clusters found in this ROC
1100 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1101
1102 AliTRDcluster *cluster = (AliTRDcluster *)
1103 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1104 row = cluster->GetLabel(0);
1105 col = cluster->GetLabel(1);
1106 time = cluster->GetLabel(2);
1107
1108 for (iPad = 0; iPad < kNclus; iPad++) {
1109 Int_t iPadCol = col - 1 + iPad;
1110 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1111 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1112 }
1113
1114 }
1115
1116 }
1117
1118 // Copy the track indices into the cluster
1119 // Loop though the clusters found in this ROC
1120 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1121
1122 AliTRDcluster *cluster = (AliTRDcluster *)
1123 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1124 cluster->SetLabel(-9999,0);
1125 cluster->SetLabel(-9999,1);
1126 cluster->SetLabel(-9999,2);
1127
1128 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1129
1130 }
1131
1132 delete [] idxTracks;
1133
1134 return kTRUE;
1135
1136}
1137
1138//_____________________________________________________________________________
acc49af9 1139Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
3fe61b77 1140{
1141 //
1142 // Get COG position
1143 // Used for clusters with more than 3 pads - where LUT not applicable
1144 //
1145
1146 Double_t sum = signal[0]
1147 + signal[1]
1148 + signal[2]
1149 + signal[3]
1150 + signal[4];
1151
023b669c 1152 // ???????????? CBL
3fe61b77 1153 Double_t res = (0.0 * (-signal[0] + signal[4])
1154 + (-signal[1] + signal[3])) / sum;
1155
1156 return res;
1157
1158}
1159
1160//_____________________________________________________________________________
053767a4 1161Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal)
3fe61b77 1162{
1163 //
1164 // Method to unfold neighbouring maxima.
1165 // The charge ratio on the overlapping pad is calculated
1166 // until there is no more change within the range given by eps.
1167 // The resulting ratio is then returned to the calling method.
1168 //
1169
1170 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1171 if (!calibration) {
3fe61b77 1172 AliError("No AliTRDcalibDB instance available\n");
1173 return kFALSE;
6d50f529 1174 }
3fe61b77 1175
1176 Int_t irc = 0;
1177 Int_t itStep = 0; // Count iteration steps
1178
1179 Double_t ratio = 0.5; // Start value for ratio
1180 Double_t prevRatio = 0.0; // Store previous ratio
1181
1182 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1183 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1184 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1185
1186 // Start the iteration
1187 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1188
1189 itStep++;
1190 prevRatio = ratio;
1191
1192 // Cluster position according to charge ratio
1193 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1194 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1195 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1196 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1197
1198 // Set cluster charge ratio
053767a4 1199 irc = calibration->PadResponse(1.0,maxLeft ,layer,newSignal);
3fe61b77 1200 Double_t ampLeft = padSignal[1] / newSignal[1];
053767a4 1201 irc = calibration->PadResponse(1.0,maxRight,layer,newSignal);
3fe61b77 1202 Double_t ampRight = padSignal[3] / newSignal[1];
1203
1204 // Apply pad response to parameters
053767a4 1205 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1206 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
3fe61b77 1207
1208 // Calculate new overlapping ratio
1209 ratio = TMath::Min((Double_t) 1.0
1210 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1211
b43a3e17 1212 }
88719a08 1213
3fe61b77 1214 return ratio;
1215
1216}
88719a08 1217
3fe61b77 1218//_____________________________________________________________________________
fb44eb8e 1219void AliTRDclusterizer::TailCancelation(AliTRDdataArrayDigits *digitsIn
3fe61b77 1220 , AliTRDdataArrayF *digitsOut
1221 , AliTRDSignalIndex *indexesIn
1222 , AliTRDSignalIndex *indexesOut
1223 , Int_t nTimeTotal
acc49af9 1224 , Float_t adcThreshold
3fe61b77 1225 , AliTRDCalROC *calGainFactorROC
1226 , Float_t calGainFactorDetValue)
1227{
1228 //
1229 // Applies the tail cancelation and gain factors:
1230 // Transform digitsIn to digitsOut
1231 //
1232
1233 Int_t iRow = 0;
1234 Int_t iCol = 0;
1235 Int_t iTime = 0;
1236
3fe61b77 1237 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
1238 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
1239 indexesIn->ResetCounters();
1240 while (indexesIn->NextRCIndex(iRow, iCol))
1241 {
1242 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1243 Double_t gain = calGainFactorDetValue
1244 * calGainFactorROCValue;
1245
f5375dcb 1246 Bool_t corrupted = kFALSE;
3fe61b77 1247 for (iTime = 0; iTime < nTimeTotal; iTime++)
1248 {
1249 // Apply gain gain factor
1250 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
f5375dcb 1251 if(digitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
3fe61b77 1252 inADC[iTime] /= gain;
1253 outADC[iTime] = inADC[iTime];
1254 }
f5375dcb 1255 if(!corrupted)
1256 {
1257 // Apply the tail cancelation via the digital filter
1258 // (only for non-coorupted pads)
3a039a31 1259 if (fReconstructor->GetRecoParam() ->IsTailCancelation())
f5375dcb 1260 {
3a039a31 1261 DeConvExp(inADC,outADC,nTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
f5375dcb 1262 }
1263 }
3fe61b77 1264
1265 indexesIn->ResetTbinCounter();
1266 while (indexesIn->NextTbinIndex(iTime))
1267 {
1268 // Store the amplitude of the digit if above threshold
acc49af9 1269 if (outADC[iTime] > adcThreshold)
3fe61b77 1270 {
1271 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1272 indexesOut->AddIndexTBin(iRow,iCol,iTime);
1273 }
1274 } // while itime
1275
1276 } // while irow icol
1277
1278 delete [] inADC;
1279 delete [] outADC;
1280
1281 return;
1282
1283}
1284
1285//_____________________________________________________________________________
1286void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1287 , Int_t n, Int_t nexp)
1288{
1289 //
1290 // Tail cancellation by deconvolution for PASA v4 TRF
1291 //
1292
1293 Double_t rates[2];
1294 Double_t coefficients[2];
1295
1296 // Initialization (coefficient = alpha, rates = lambda)
1297 Double_t r1 = 1.0;
1298 Double_t r2 = 1.0;
1299 Double_t c1 = 0.5;
1300 Double_t c2 = 0.5;
1301
1302 if (nexp == 1) { // 1 Exponentials
1303 r1 = 1.156;
1304 r2 = 0.130;
1305 c1 = 0.066;
1306 c2 = 0.000;
1307 }
1308 if (nexp == 2) { // 2 Exponentials
181c7f7e 1309 Double_t par[4];
1310 fReconstructor->GetTCParams(par);
1311 r1 = par[0];//1.156;
1312 r2 = par[1];//0.130;
1313 c1 = par[2];//0.114;
1314 c2 = par[3];//0.624;
3fe61b77 1315 }
1316
1317 coefficients[0] = c1;
1318 coefficients[1] = c2;
1319
1320 Double_t dt = 0.1;
1321
1322 rates[0] = TMath::Exp(-dt/(r1));
1323 rates[1] = TMath::Exp(-dt/(r2));
1324
1325 Int_t i = 0;
1326 Int_t k = 0;
1327
1328 Double_t reminder[2];
1329 Double_t correction = 0.0;
1330 Double_t result = 0.0;
1331
1332 // Attention: computation order is important
1333 for (k = 0; k < nexp; k++) {
1334 reminder[k] = 0.0;
1335 }
1336
1337 for (i = 0; i < n; i++) {
1338
1339 result = (source[i] - correction); // No rescaling
1340 target[i] = result;
1341
1342 for (k = 0; k < nexp; k++) {
1343 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1344 }
1345
1346 correction = 0.0;
1347 for (k = 0; k < nexp; k++) {
1348 correction += reminder[k];
1349 }
1350
1351 }
6d50f529 1352
1353}
1354
1355//_____________________________________________________________________________
1356void AliTRDclusterizer::ResetRecPoints()
1357{
1358 //
1359 // Resets the list of rec points
1360 //
1361
1362 if (fRecPoints) {
1363 fRecPoints->Delete();
48f8adf3 1364 delete fRecPoints;
6d50f529 1365 }
6d50f529 1366}
1367
1368//_____________________________________________________________________________
66f6bfd9 1369TClonesArray *AliTRDclusterizer::RecPoints()
6d50f529 1370{
1371 //
1372 // Returns the list of rec points
1373 //
1374
1375 if (!fRecPoints) {
48f8adf3 1376 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1377 fRecPoints = new TClonesArray("AliTRDcluster", 400);
1378 }
1379 //SetClustersOwner(kTRUE);
1380 AliTRDReconstructor::SetClusters(0x0);
6d50f529 1381 }
6d50f529 1382 return fRecPoints;
3551db50 1383}
fc546d21 1384
1385//_____________________________________________________________________________
1386void AliTRDclusterizer::FillLUT()
1387{
1388 //
1389 // Create the LUT
1390 //
1391
1392 const Int_t kNlut = 128;
1393
053767a4 1394 fLUTbin = AliTRDgeometry::kNlayer * kNlut;
fc546d21 1395
1396 // The lookup table from Bogdan
053767a4 1397 Float_t lut[AliTRDgeometry::kNlayer][kNlut] = {
fc546d21 1398 {
1399 0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1400 0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1401 0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1402 0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1403 0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1404 0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1405 0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1406 0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1407 0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1408 0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1409 0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1410 0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1411 0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1412 0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1413 0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1414 0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1415 },
1416 {
1417 0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1418 0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1419 0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1420 0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1421 0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1422 0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1423 0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1424 0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1425 0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1426 0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1427 0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1428 0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1429 0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1430 0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1431 0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1432 0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1433 },
1434 {
1435 0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1436 0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1437 0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1438 0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1439 0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1440 0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1441 0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1442 0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1443 0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1444 0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1445 0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1446 0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1447 0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1448 0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1449 0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1450 0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1451 },
1452 {
1453 0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1454 0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1455 0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1456 0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1457 0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1458 0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1459 0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1460 0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1461 0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1462 0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1463 0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1464 0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1465 0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1466 0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1467 0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1468 0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1469 },
1470 {
1471 0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1472 0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1473 0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1474 0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1475 0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1476 0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1477 0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1478 0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1479 0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1480 0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1481 0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1482 0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1483 0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1484 0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1485 0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1486 0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1487 },
1488 {
1489 0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1490 0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1491 0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1492 0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1493 0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1494 0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1495 0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1496 0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1497 0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1498 0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1499 0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1500 0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1501 0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1502 0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1503 0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1504 0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1505 }
1506 };
1507
1508 if (fLUT) {
1509 delete [] fLUT;
1510 }
1511 fLUT = new Double_t[fLUTbin];
1512
053767a4 1513 for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
fc546d21 1514 for (Int_t ilut = 0; ilut < kNlut; ilut++ ) {
053767a4 1515 fLUT[ilayer*kNlut+ilut] = lut[ilayer][ilut];
fc546d21 1516 }
1517 }
1518
1519}
1520
1521//_____________________________________________________________________________
053767a4 1522Double_t AliTRDclusterizer::LUTposition(Int_t ilayer, Double_t ampL
fc546d21 1523 , Double_t ampC, Double_t ampR) const
1524{
1525 //
1526 // Calculates the cluster position using the lookup table.
1527 // Method provided by Bogdan Vulpescu.
1528 //
1529
1530 const Int_t kNlut = 128;
1531
1532 Double_t pos;
1533 Double_t x = 0.0;
1534 Double_t xmin;
1535 Double_t xmax;
1536 Double_t xwid;
1537
1538 Int_t side = 0;
1539 Int_t ix;
1540
053767a4 1541 Double_t xMin[AliTRDgeometry::kNlayer] = { 0.006492, 0.006377, 0.006258
1542 , 0.006144, 0.006030, 0.005980 };
1543 Double_t xMax[AliTRDgeometry::kNlayer] = { 0.960351, 0.965870, 0.970445
1544 , 0.974352, 0.977667, 0.996101 };
fc546d21 1545
1546 if (ampL > ampR) {
1547 x = (ampL - ampR) / ampC;
1548 side = -1;
1549 }
1550 else if (ampL < ampR) {
1551 x = (ampR - ampL) / ampC;
1552 side = +1;
1553 }
1554
1555 if (ampL != ampR) {
1556
053767a4 1557 xmin = xMin[ilayer] + 0.000005;
1558 xmax = xMax[ilayer] - 0.000005;
fc546d21 1559 xwid = (xmax - xmin) / 127.0;
1560
1561 if (x < xmin) {
1562 pos = 0.0000;
1563 }
1564 else if (x > xmax) {
1565 pos = side * 0.5000;
1566 }
1567 else {
1568 ix = (Int_t) ((x - xmin) / xwid);
053767a4 1569 pos = side * fLUT[ilayer*kNlut+ix];
fc546d21 1570 }
1571
1572 }
1573 else {
1574
1575 pos = 0.0;
1576
1577 }
1578
1579 return pos;
1580
1581}