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