]> 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
Added OutsideThreshold function. Corrections to TTimeStamp syntax. Set startTime...
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizer.cxx
CommitLineData
3fe61b77 1
f7336fa3 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
88cb7938 17/* $Id$ */
f7336fa3 18
19///////////////////////////////////////////////////////////////////////////////
20// //
3fe61b77 21// TRD cluster finder //
f7336fa3 22// //
23///////////////////////////////////////////////////////////////////////////////
24
3fe61b77 25#include <TF1.h>
94de3818 26#include <TTree.h>
3fe61b77 27#include <TH1.h>
793ff80c 28#include <TFile.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"
793ff80c 39#include "AliTRDgeometry.h"
3fe61b77 40#include "AliTRDdataArrayF.h"
41#include "AliTRDdataArrayI.h"
625f5260 42#include "AliTRDdataArrayS.h"
3fe61b77 43#include "AliTRDdigitsManager.h"
44#include "AliTRDrawData.h"
3551db50 45#include "AliTRDcalibDB.h"
3fe61b77 46#include "AliTRDSimParam.h"
47#include "AliTRDRecParam.h"
b43a3e17 48#include "AliTRDCommonParam.h"
3fe61b77 49#include "AliTRDtransform.h"
50#include "AliTRDSignalIndex.h"
dfbb4bb9 51#include "AliTRDrawStreamBase.h"
52// #include "AliTRDRawStream.h"
53// #include "AliTRDRawStreamV2.h"
3fe61b77 54#include "AliTRDfeeParam.h"
55
56#include "Cal/AliTRDCalROC.h"
57#include "Cal/AliTRDCalDet.h"
f7336fa3 58
59ClassImp(AliTRDclusterizer)
60
61//_____________________________________________________________________________
6d50f529 62AliTRDclusterizer::AliTRDclusterizer()
63 :TNamed()
64 ,fRunLoader(NULL)
65 ,fClusterTree(NULL)
66 ,fRecPoints(NULL)
3fe61b77 67 ,fDigitsManager(NULL)
68 ,fAddLabels(kTRUE)
69 ,fRawVersion(2)
70 ,fIndexesOut(NULL)
71 ,fIndexesMaxima(NULL)
dfbb4bb9 72 ,fTransform(new AliTRDtransform(0))
f7336fa3 73{
74 //
75 // AliTRDclusterizer default constructor
76 //
77
3fe61b77 78 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
79
f7336fa3 80}
81
82//_____________________________________________________________________________
3fe61b77 83AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title)
6d50f529 84 :TNamed(name,title)
85 ,fRunLoader(NULL)
86 ,fClusterTree(NULL)
87 ,fRecPoints(NULL)
3fe61b77 88 ,fDigitsManager(new AliTRDdigitsManager())
89 ,fAddLabels(kTRUE)
90 ,fRawVersion(2)
91 ,fIndexesOut(NULL)
92 ,fIndexesMaxima(NULL)
93 ,fTransform(new AliTRDtransform(0))
f7336fa3 94{
95 //
6d50f529 96 // AliTRDclusterizer constructor
f7336fa3 97 //
98
3fe61b77 99 fDigitsManager->CreateArrays();
100
101 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
102
f7336fa3 103}
104
8230f242 105//_____________________________________________________________________________
6d50f529 106AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
107 :TNamed(c)
108 ,fRunLoader(NULL)
109 ,fClusterTree(NULL)
110 ,fRecPoints(NULL)
3fe61b77 111 ,fDigitsManager(NULL)
112 ,fAddLabels(kTRUE)
113 ,fRawVersion(2)
114 ,fIndexesOut(NULL)
115 ,fIndexesMaxima(NULL)
116 ,fTransform(NULL)
8230f242 117{
118 //
119 // AliTRDclusterizer copy constructor
120 //
121
8230f242 122}
123
f7336fa3 124//_____________________________________________________________________________
125AliTRDclusterizer::~AliTRDclusterizer()
126{
8230f242 127 //
128 // AliTRDclusterizer destructor
129 //
f7336fa3 130
3fe61b77 131 if (fRecPoints)
132 {
133 fRecPoints->Delete();
134 delete fRecPoints;
135 }
136
137 if (fDigitsManager)
138 {
139 delete fDigitsManager;
140 fDigitsManager = NULL;
141 }
142
143 if (fIndexesOut)
144 {
145 delete fIndexesOut;
146 fIndexesOut = NULL;
147 }
148
149 if (fIndexesMaxima)
150 {
151 delete fIndexesMaxima;
152 fIndexesMaxima = NULL;
153 }
154
155 if (fTransform)
156 {
157 delete fTransform;
158 fTransform = NULL;
159 }
6d50f529 160
f7336fa3 161}
162
8230f242 163//_____________________________________________________________________________
dd9a6ee3 164AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
165{
166 //
167 // Assignment operator
168 //
169
3fe61b77 170 if (this != &c)
171 {
172 ((AliTRDclusterizer &) c).Copy(*this);
173 }
dd9a6ee3 174 return *this;
175
176}
177
178//_____________________________________________________________________________
e0d47c25 179void AliTRDclusterizer::Copy(TObject &c) const
8230f242 180{
181 //
182 // Copy function
183 //
184
3fe61b77 185 ((AliTRDclusterizer &) c).fClusterTree = NULL;
186 ((AliTRDclusterizer &) c).fRecPoints = NULL;
187 ((AliTRDclusterizer &) c).fDigitsManager = NULL;
188 ((AliTRDclusterizer &) c).fAddLabels = fAddLabels;
189 ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
190 ((AliTRDclusterizer &) c).fIndexesOut = NULL;
191 ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
192 ((AliTRDclusterizer &) c).fTransform = NULL;
8230f242 193
194}
195
f7336fa3 196//_____________________________________________________________________________
3e1a3ad8 197Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
f7336fa3 198{
199 //
3e1a3ad8 200 // Opens the AliROOT file. Output and input are in the same file
f7336fa3 201 //
6d50f529 202
e191bb57 203 TString evfoldname = AliConfig::GetDefaultEventFolderName();
6d50f529 204 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
205
206 if (!fRunLoader) {
19dd5b2f 207 fRunLoader = AliRunLoader::Open(name);
6d50f529 208 }
209
210 if (!fRunLoader) {
211 AliError(Form("Can not open session for file %s.",name));
212 return kFALSE;
213 }
88cb7938 214
215 OpenInput(nEvent);
216 OpenOutput();
6d50f529 217
3e1a3ad8 218 return kTRUE;
f7336fa3 219
6d50f529 220}
3e1a3ad8 221
222//_____________________________________________________________________________
88cb7938 223Bool_t AliTRDclusterizer::OpenOutput()
3e1a3ad8 224{
225 //
226 // Open the output file
227 //
228
3e1a3ad8 229 TObjArray *ioArray = 0;
88cb7938 230
231 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
232 loader->MakeTree("R");
6d50f529 233
88cb7938 234 fClusterTree = loader->TreeR();
365d0374 235 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
3e1a3ad8 236
3e1a3ad8 237 return kTRUE;
238
239}
240
25ca55ce 241//_____________________________________________________________________________
242Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
243{
244 //
245 // Connect the output tree
246 //
247
248 TObjArray *ioArray = 0;
249
250 fClusterTree = clusterTree;
251 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
252
253 return kTRUE;
254
255}
256
3e1a3ad8 257//_____________________________________________________________________________
88cb7938 258Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
f7336fa3 259{
260 //
261 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
262 //
263
f7336fa3 264 // Import the Trees for the event nEvent in the file
bdbb05bb 265 fRunLoader->GetEvent(nEvent);
88cb7938 266
f7336fa3 267 return kTRUE;
268
269}
270
271//_____________________________________________________________________________
793ff80c 272Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
f7336fa3 273{
274 //
3e1a3ad8 275 // Fills TRDcluster branch in the tree with the clusters
793ff80c 276 // found in detector = det. For det=-1 writes the tree.
a3c76cdc 277 //
793ff80c 278
6d50f529 279 if ((det < -1) ||
280 (det >= AliTRDgeometry::Ndet())) {
281 AliError(Form("Unexpected detector index %d.\n",det));
3e1a3ad8 282 return kFALSE;
793ff80c 283 }
3e1a3ad8 284
3e1a3ad8 285 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
286 if (!branch) {
793ff80c 287 TObjArray *ioArray = 0;
3e1a3ad8 288 branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
793ff80c 289 }
290
6d50f529 291 if ((det >= 0) &&
292 (det < AliTRDgeometry::Ndet())) {
793ff80c 293
bdbb05bb 294 Int_t nRecPoints = RecPoints()->GetEntriesFast();
3e1a3ad8 295 TObjArray *detRecPoints = new TObjArray(400);
296
297 for (Int_t i = 0; i < nRecPoints; i++) {
bdbb05bb 298 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
3e1a3ad8 299 if (det == c->GetDetector()) {
300 detRecPoints->AddLast(c);
301 }
302 else {
3fe61b77 303 AliError(Form("Attempt to write a cluster with unexpected detector index: got=%d expected=%d\n"
304 ,c->GetDetector()
305 ,det));
3e1a3ad8 306 }
793ff80c 307 }
308
3e1a3ad8 309 branch->SetAddress(&detRecPoints);
310 fClusterTree->Fill();
311
d9b8978b 312 delete detRecPoints;
ecf39416 313
793ff80c 314 return kTRUE;
3e1a3ad8 315
793ff80c 316 }
317
318 if (det == -1) {
319
6d50f529 320 AliInfo(Form("Writing the cluster tree %s for event %d."
321 ,fClusterTree->GetName(),fRunLoader->GetEventNumber()));
322
a6dd11e9 323 if (fRecPoints) {
324
325 branch->SetAddress(&fRecPoints);
326
327 AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
328 loader->WriteRecPoints("OVERWRITE");
c630aafd 329
a6dd11e9 330 }
331 else {
332
333 AliError("Cluster tree does not exist. Cannot write clusters.\n");
334 return kFALSE;
335
336 }
337
793ff80c 338 return kTRUE;
3e1a3ad8 339
793ff80c 340 }
6d50f529 341
342 AliError(Form("Unexpected detector index %d.\n",det));
3e1a3ad8 343
793ff80c 344 return kFALSE;
88cb7938 345
f7336fa3 346}
793ff80c 347
3551db50 348//_____________________________________________________________________________
3fe61b77 349void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
350{
351 //
352 // Reset the helper indexes
353 //
354
355 if (fIndexesOut)
356 {
357 // carefull here - we assume that only row number may change - most probable
358 if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
359 fIndexesOut->ResetContent();
360 else
361 fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
362 , indexesIn->GetNcol()
363 , indexesIn->GetNtime());
364 }
365 else
366 {
367 fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
368 , indexesIn->GetNcol()
369 , indexesIn->GetNtime());
370 }
371
372 if (fIndexesMaxima)
373 {
374 // carefull here - we assume that only row number may change - most probable
375 if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
376 {
377 fIndexesMaxima->ResetContent();
378 }
379 else
380 {
381 fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
382 , indexesIn->GetNcol()
383 , indexesIn->GetNtime());
384 }
385 }
386 else
387 {
388 fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
389 , indexesIn->GetNcol()
390 , indexesIn->GetNtime());
391 }
392
393}
394
395//_____________________________________________________________________________
396Bool_t AliTRDclusterizer::ReadDigits()
397{
398 //
399 // Reads the digits arrays from the input aliroot file
400 //
401
402 if (!fRunLoader) {
403 AliError("No run loader available");
404 return kFALSE;
405 }
406
407 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
408 if (!loader->TreeD()) {
409 loader->LoadDigits();
410 }
411
412 // Read in the digit arrays
413 return (fDigitsManager->ReadDigits(loader->TreeD()));
414
415}
416
417//_____________________________________________________________________________
418Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
419{
420 //
421 // Reads the digits arrays from the input tree
422 //
423
424 // Read in the digit arrays
425 return (fDigitsManager->ReadDigits(digitsTree));
426
427}
428
429//_____________________________________________________________________________
430Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
431{
432 //
433 // Reads the digits arrays from the ddl file
434 //
435
436 AliTRDrawData raw;
437 fDigitsManager = raw.Raw2Digits(rawReader);
438
439 return kTRUE;
440
441}
442
443//_____________________________________________________________________________
444Bool_t AliTRDclusterizer::MakeClusters()
445{
446 //
447 // Creates clusters from digits
448 //
449
450 // Propagate info from the digits manager
451 if (fAddLabels == kTRUE)
452 {
453 fAddLabels = fDigitsManager->UsesDictionaries();
454 }
455
456 Bool_t fReturn = kTRUE;
457 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
458 {
459
625f5260 460 AliTRDdataArrayS *digitsIn = (AliTRDdataArrayS *) fDigitsManager->GetDigits(i);
3fe61b77 461 // This is to take care of switched off super modules
625f5260 462 if (!digitsIn->HasData())
3fe61b77 463 {
464 continue;
465 }
466 digitsIn->Expand();
467 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
468 if (indexes->IsAllocated() == kFALSE)
469 {
470 fDigitsManager->BuildIndexes(i);
471 }
472
473 Bool_t fR = kFALSE;
474 if (indexes->HasEntry())
475 {
476 if (fAddLabels)
477 {
478 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
479 {
480 AliTRDdataArrayI *tracksIn = 0;
625f5260 481 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(i,iDict);
3fe61b77 482 tracksIn->Expand();
483 }
484 }
485 fR = MakeClusters(i);
486 fReturn = fR && fReturn;
487 }
488
489 if (fR == kFALSE)
490 {
491 WriteClusters(i);
492 ResetRecPoints();
493 }
494
495 // No compress just remove
496 fDigitsManager->RemoveDigits(i);
497 fDigitsManager->RemoveDictionaries(i);
498 fDigitsManager->ClearIndexes(i);
499
500 }
501
502 return fReturn;
503
504}
505
506//_____________________________________________________________________________
507Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
508{
509 //
510 // Creates clusters from raw data
511 //
512
dfbb4bb9 513// AliTRDdataArrayS *digits = 0;
514// AliTRDdataArrayI *track0 = 0;
515// AliTRDdataArrayI *track1 = 0;
516// AliTRDdataArrayI *track2 = 0;
3fe61b77 517
dfbb4bb9 518// AliTRDSignalIndex *indexes = 0;
3fe61b77 519
dfbb4bb9 520// // Create the digits manager
521// if (!fDigitsManager)
522// {
523// fDigitsManager = new AliTRDdigitsManager();
524// fDigitsManager->CreateArrays();
525// }
3fe61b77 526
dfbb4bb9 527// // AliTRDRawStreamV2 input(rawReader);
528// // input.SetRawVersion( fRawVersion );
529// // input.Init();
530// AliTRDrawStreamBase &input = *AliTRDrawStreamBase::GetRawStream(rawReader);
3fe61b77 531
dfbb4bb9 532// AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
3fe61b77 533
dfbb4bb9 534// // Loop through the digits
535// Int_t lastdet = -1;
536// Int_t det = 0;
537// Int_t it = 0;
538// while (input.Next())
539// {
3fe61b77 540
dfbb4bb9 541// det = input.GetDet();
3fe61b77 542
dfbb4bb9 543// if (det != lastdet)
544// {
3fe61b77 545
dfbb4bb9 546// if (lastdet != -1)
547// {
548// digits = (AliTRDdataArrayS *) fDigitsManager->GetDigits(lastdet);
549// Bool_t iclusterBranch = kFALSE;
550// if (indexes->HasEntry())
551// iclusterBranch = MakeClusters(lastdet);
552// if (iclusterBranch == kFALSE)
553// {
554// WriteClusters(lastdet);
555// ResetRecPoints();
556// }
557// }
558
559// if (digits)
560// {
561// fDigitsManager->RemoveDigits(lastdet);
562// fDigitsManager->RemoveDictionaries(lastdet);
563// fDigitsManager->ClearIndexes(lastdet);
564// }
565
566// lastdet = det;
567
568// // Add a container for the digits of this detector
569// digits = (AliTRDdataArrayS *) fDigitsManager->GetDigits(det);
570// track0 = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(det,0);
571// track1 = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(det,1);
572// track2 = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(det,2);
573
574// // Allocate memory space for the digits buffer
575// if (!digits->HasData())
576// {
577// //AliDebug(5, Form("Alloc digits for det %d", det));
578// digits->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
579// track0->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
580// track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
581// track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
582// }
3fe61b77 583
dfbb4bb9 584// indexes = fDigitsManager->GetIndexes(det);
585// indexes->SetSM(input.GetSM());
586// indexes->SetStack(input.GetStack());
587// indexes->SetLayer(input.GetLayer());
588// indexes->SetDetNumber(det);
589// if (indexes->IsAllocated() == kFALSE)
590// {
591// indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins());
592// }
593
594// }
3fe61b77 595
dfbb4bb9 596// for (it = 0; it < 3; it++)
597// {
598// if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() )
599// {
600// if (input.GetSignals()[it] > 0)
601// {
602// digits->SetDataUnchecked(input.GetRow(), input.GetCol(),
603// input.GetTimeBin() + it, input.GetSignals()[it]);
604
605// indexes->AddIndexTBin(input.GetRow(), input.GetCol(),
606// input.GetTimeBin() + it);
607// track0->SetDataUnchecked(input.GetRow(), input.GetCol(),
608// input.GetTimeBin() + it, 0);
609// track1->SetDataUnchecked(input.GetRow(), input.GetCol(),
610// input.GetTimeBin() + it, 0);
611// track2->SetDataUnchecked(input.GetRow(), input.GetCol(),
612// input.GetTimeBin() + it, 0);
613// }
614// }
615// }
616
617// }
618
619// if (lastdet != -1)
620// {
621// Bool_t iclusterBranch = kFALSE;
622// if (indexes->HasEntry())
623// {
624// iclusterBranch = MakeClusters(lastdet);
625// }
626// if (iclusterBranch == kFALSE)
627// {
628// WriteClusters(lastdet);
629// ResetRecPoints();
630// }
631// //MakeClusters(lastdet);
632// if (digits)
633// {
634// fDigitsManager->RemoveDigits(lastdet);
635// fDigitsManager->RemoveDictionaries(lastdet);
636// fDigitsManager->ClearIndexes(lastdet);
637// }
638// }
639
640// delete fDigitsManager;
641// fDigitsManager = NULL;
642// return kTRUE;
643
644return Raw2ClustersChamber(rawReader);
3fe61b77 645}
646
647//_____________________________________________________________________________
648Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
649{
650 //
651 // Creates clusters from raw data
652 //
653
654 // Create the digits manager
655 if (!fDigitsManager)
656 {
657 fDigitsManager = new AliTRDdigitsManager();
658 fDigitsManager->CreateArrays();
659 }
660
661 fDigitsManager->SetUseDictionaries(fAddLabels);
662
dfbb4bb9 663// AliTRDRawStreamV2 input(rawReader);
664// input.SetRawVersion( fRawVersion );
665// input.Init();
666 AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
667 AliTRDrawStreamBase &input = *pinput;
3fe61b77 668
669 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
670
671 Int_t det = 0;
672 while ((det = input.NextChamber(fDigitsManager)) >= 0)
673 {
674 Bool_t iclusterBranch = kFALSE;
675 if (fDigitsManager->GetIndexes(det)->HasEntry())
676 {
677 iclusterBranch = MakeClusters(det);
678 }
679 if (iclusterBranch == kFALSE)
680 {
681 WriteClusters(det);
682 ResetRecPoints();
683 }
684 fDigitsManager->RemoveDigits(det);
685 fDigitsManager->RemoveDictionaries(det);
686 fDigitsManager->ClearIndexes(det);
687 }
688
689 delete fDigitsManager;
690 fDigitsManager = NULL;
dfbb4bb9 691
692 delete pinput;
693 pinput = NULL;
3fe61b77 694 return kTRUE;
695
696}
697
698//_____________________________________________________________________________
699Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
700{
701 //
702 // Generates the cluster.
703 //
704
705 // Get the digits
706 // digits should be expanded beforehand!
707 // digitsIn->Expand();
625f5260 708 AliTRDdataArrayS *digitsIn = (AliTRDdataArrayS *) fDigitsManager->GetDigits(det);
3fe61b77 709
710 // This is to take care of switched off super modules
625f5260 711 if (!digitsIn->HasData())
3fe61b77 712 {
713 return kFALSE;
714 }
715
716 AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
717 if (indexesIn->IsAllocated() == kFALSE)
718 {
719 AliError("Indexes do not exist!");
720 return kFALSE;
721 }
722
723 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
724 if (!calibration)
725 {
726 AliFatal("No AliTRDcalibDB instance available\n");
727 return kFALSE;
728 }
729
730 AliTRDRecParam *recParam = AliTRDRecParam::Instance();
731 if (!recParam)
732 {
733 AliError("No AliTRDRecParam instance available\n");
734 return kFALSE;
735 }
736
737 // ADC thresholds
738 // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
acc49af9 739 Float_t adcThreshold = 0;
3fe61b77 740
741 // Threshold value for the maximum
742 Float_t maxThresh = recParam->GetClusMaxThresh();
743 // Threshold value for the digit signal
744 Float_t sigThresh = recParam->GetClusSigThresh();
745
746 // Iteration limit for unfolding procedure
747 const Float_t kEpsilon = 0.01;
748 const Int_t kNclus = 3;
749 const Int_t kNsig = 5;
750
751 Int_t iUnfold = 0;
752 Double_t ratioLeft = 1.0;
753 Double_t ratioRight = 1.0;
754
755 Double_t padSignal[kNsig];
756 Double_t clusterSignal[kNclus];
757
758 Int_t icham = indexesIn->GetChamber();
759 Int_t iplan = indexesIn->GetPlane();
760 Int_t isect = indexesIn->GetSM();
761
762 // Start clustering in the chamber
763
764 Int_t idet = AliTRDgeometry::GetDetector(iplan,icham,isect);
765 if (idet != det)
766 {
767 AliError("Strange Detector number Missmatch!");
768 return kFALSE;
769 }
770
771 // TRD space point transformation
772 fTransform->SetDetector(det);
773
774 Int_t ilayer = AliGeomManager::kTRD1 + iplan;
775 Int_t imodule = icham + AliTRDgeometry::Ncham() * isect;
776 UShort_t volid = AliGeomManager::LayerToVolUID(ilayer,imodule);
777
778 Int_t nColMax = digitsIn->GetNcol();
779 Int_t nTimeTotal = digitsIn->GetNtime();
780
781 // Detector wise calibration object for the gain factors
782 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
783 // Calibration object with pad wise values for the gain factors
784 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
785 // Calibration value for chamber wise gain factor
786 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
787
788 Int_t nClusters = 0;
789
790 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow()
791 ,digitsIn->GetNcol()
792 ,digitsIn->GetNtime());
793
794 ResetHelperIndexes(indexesIn);
795
796 // Apply the gain and the tail cancelation via digital filter
797 TailCancelation(digitsIn
798 ,digitsOut
799 ,indexesIn
800 ,fIndexesOut
801 ,nTimeTotal
acc49af9 802 ,adcThreshold
3fe61b77 803 ,calGainFactorROC
804 ,calGainFactorDetValue);
805
806 Int_t row = 0;
807 Int_t col = 0;
808 Int_t time = 0;
809 Int_t iPad = 0;
810
811 fIndexesOut->ResetCounters();
812 while (fIndexesOut->NextRCTbinIndex(row, col, time))
813 {
814
815 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
816
817 // Look for the maximum
818 if (signalM >= maxThresh)
819 {
820
821 if (col + 1 >= nColMax || col-1 < 0)
822 continue;
823
824 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
825 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
826
827 if ((TMath::Abs(signalL) <= signalM) &&
828 (TMath::Abs(signalR) < signalM))
829 {
830 if ((TMath::Abs(signalL) >= sigThresh) ||
831 (TMath::Abs(signalR) >= sigThresh))
832 {
833 // Maximum found, mark the position by a negative signal
834 digitsOut->SetDataUnchecked(row,col,time,-signalM);
835 fIndexesMaxima->AddIndexTBin(row,col,time);
836 }
837 }
838
839 }
840
841 }
842
843 // The index to the first cluster of a given ROC
844 Int_t firstClusterROC = -1;
845 // The number of cluster in a given ROC
846 Int_t nClusterROC = 0;
847
848 // Now check the maxima and calculate the cluster position
849 fIndexesMaxima->ResetCounters();
850 while (fIndexesMaxima->NextRCTbinIndex(row, col, time))
851 {
852
853 // Maximum found ?
854 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0)
855 {
856
857 for (iPad = 0; iPad < kNclus; iPad++)
858 {
859 Int_t iPadCol = col - 1 + iPad;
860 clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
861 }
862
863 // Count the number of pads in the cluster
864 Int_t nPadCount = 0;
865 Int_t ii;
866 // Look to the left
867 ii = 0;
868 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh)
869 {
870 nPadCount++;
871 ii++;
872 if (col-ii < 0) break;
873 }
874 // Look to the right
875 ii = 0;
876 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh)
877 {
878 nPadCount++;
879 ii++;
880 if (col+ii+1 >= nColMax) break;
881 }
882 nClusters++;
883
884 // Look for 5 pad cluster with minimum in the middle
885 Bool_t fivePadCluster = kFALSE;
886 if (col < (nColMax - 3))
887 {
888 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0)
889 {
890 fivePadCluster = kTRUE;
891 }
892 if ((fivePadCluster) && (col < (nColMax - 5)))
893 {
894 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh)
895 {
896 fivePadCluster = kFALSE;
897 }
898 }
899 if ((fivePadCluster) && (col > 1))
900 {
901 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh)
902 {
903 fivePadCluster = kFALSE;
904 }
905 }
906 }
907
908 // 5 pad cluster
909 // Modify the signal of the overlapping pad for the left part
910 // of the cluster which remains from a previous unfolding
911 if (iUnfold)
912 {
913 clusterSignal[0] *= ratioLeft;
914 iUnfold = 0;
915 }
916
917 // Unfold the 5 pad cluster
918 if (fivePadCluster)
919 {
920 for (iPad = 0; iPad < kNsig; iPad++)
921 {
922 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
923 ,col-1+iPad
924 ,time));
925 }
926 // Unfold the two maxima and set the signal on
927 // the overlapping pad to the ratio
928 ratioRight = Unfold(kEpsilon,iplan,padSignal);
929 ratioLeft = 1.0 - ratioRight;
930 clusterSignal[2] *= ratioRight;
931 iUnfold = 1;
932 }
933
934 // The position of the cluster in COL direction relative to the center pad (pad units)
935 Double_t clusterPosCol = 0.0;
936 if (recParam->LUTOn())
937 {
938 // Calculate the position of the cluster by using the
939 // lookup table method
940 clusterPosCol = recParam->LUTposition(iplan
941 ,clusterSignal[0]
942 ,clusterSignal[1]
943 ,clusterSignal[2]);
944 }
945 else
946 {
947 // Calculate the position of the cluster by using the
948 // center of gravity method
949 for (Int_t i = 0; i < kNsig; i++)
950 {
951 padSignal[i] = 0.0;
952 }
953 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
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
956 if ((col > 2) &&
957 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1]))
958 {
959 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
960 }
961 if ((col < nColMax - 3) &&
962 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3]))
963 {
964 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
965 }
966 clusterPosCol = GetCOG(padSignal);
967 }
968
969 // Store the amplitudes of the pads in the cluster for later analysis
970 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
971 for (Int_t jPad = col-3; jPad <= col+3; jPad++)
972 {
973 if ((jPad < 0) ||
974 (jPad >= nColMax-1))
975 {
976 continue;
977 }
978 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
979 }
980
981 // Transform the local cluster coordinates into calibrated
982 // space point positions defined in the local tracking system.
983 // Here the calibration for T0, Vdrift and ExB is applied as well.
984 Double_t clusterXYZ[6];
985 clusterXYZ[0] = clusterPosCol;
986 clusterXYZ[1] = clusterSignal[0];
987 clusterXYZ[2] = clusterSignal[1];
988 clusterXYZ[3] = clusterSignal[2];
989 clusterXYZ[4] = 0.0;
990 clusterXYZ[5] = 0.0;
991 Int_t clusterRCT[3];
992 clusterRCT[0] = row;
993 clusterRCT[1] = col;
994 clusterRCT[2] = 0;
bcb6fb78 995
996 Bool_t out = kTRUE;
997 if (fTransform->Transform(clusterXYZ, clusterRCT, ((UInt_t) time), out, 0)) {
9bf8c575 998
999 // Add the cluster to the output array
1000 // The track indices will be stored later
1001 Float_t clusterPos[3];
1002 clusterPos[0] = clusterXYZ[0];
1003 clusterPos[1] = clusterXYZ[1];
1004 clusterPos[2] = clusterXYZ[2];
1005 Float_t clusterSig[2];
1006 clusterSig[0] = clusterXYZ[4];
1007 clusterSig[1] = clusterXYZ[5];
1008 Double_t clusterCharge = clusterXYZ[3];
1009 Char_t clusterTimeBin = ((Char_t) clusterRCT[2]);
1010 AliTRDcluster *cluster = new AliTRDcluster(idet
1011 ,clusterCharge
1012 ,clusterPos
1013 ,clusterSig
1014 ,0x0
1015 ,((Char_t) nPadCount)
1016 ,signals
1017 ,((UChar_t) col)
1018 ,((UChar_t) row)
1019 ,((UChar_t) time)
1020 ,clusterTimeBin
1021 ,clusterPosCol
1022 ,volid);
bcb6fb78 1023 cluster->SetInChamber(!out);
1024
9bf8c575 1025 // Temporarily store the row, column and time bin of the center pad
1026 // Used to later on assign the track indices
1027 cluster->SetLabel( row,0);
1028 cluster->SetLabel( col,1);
1029 cluster->SetLabel(time,2);
1030
1031 RecPoints()->Add(cluster);
1032
1033 // Store the index of the first cluster in the current ROC
1034 if (firstClusterROC < 0)
1035 {
1036 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
1037 }
1038
1039 // Count the number of cluster in the current ROC
1040 nClusterROC++;
1041
1042 } // if: Transform ok ?
3fe61b77 1043
1044 } // if: Maximum found ?
1045
1046 }
1047
1048 delete digitsOut;
1049
1050 if (fAddLabels)
1051 {
1052 AddLabels(idet, firstClusterROC, nClusterROC);
1053 }
1054
1055 // Write the cluster and reset the array
1056 WriteClusters(idet);
1057 ResetRecPoints();
1058
1059 return kTRUE;
1060
1061}
1062
1063//_____________________________________________________________________________
1064Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
3551db50 1065{
1066 //
3fe61b77 1067 // Add the track indices to the found clusters
3551db50 1068 //
1069
3fe61b77 1070 const Int_t kNclus = 3;
1071 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1072 const Int_t kNtrack = kNdict * kNclus;
1073
1074 Int_t iClusterROC = 0;
1075
1076 Int_t row = 0;
1077 Int_t col = 0;
1078 Int_t time = 0;
1079 Int_t iPad = 0;
1080
1081 // Temporary array to collect the track indices
1082 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1083
1084 // Loop through the dictionary arrays one-by-one
1085 // to keep memory consumption low
1086 AliTRDdataArrayI *tracksIn = 0;
1087 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1088
1089 // tracksIn should be expanded beforehand!
625f5260 1090 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
3fe61b77 1091
1092 // Loop though the clusters found in this ROC
1093 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1094
1095 AliTRDcluster *cluster = (AliTRDcluster *)
1096 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1097 row = cluster->GetLabel(0);
1098 col = cluster->GetLabel(1);
1099 time = cluster->GetLabel(2);
1100
1101 for (iPad = 0; iPad < kNclus; iPad++) {
1102 Int_t iPadCol = col - 1 + iPad;
1103 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1104 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1105 }
1106
1107 }
1108
1109 }
1110
1111 // Copy the track indices into the cluster
1112 // Loop though the clusters found in this ROC
1113 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1114
1115 AliTRDcluster *cluster = (AliTRDcluster *)
1116 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1117 cluster->SetLabel(-9999,0);
1118 cluster->SetLabel(-9999,1);
1119 cluster->SetLabel(-9999,2);
1120
1121 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1122
1123 }
1124
1125 delete [] idxTracks;
1126
1127 return kTRUE;
1128
1129}
1130
1131//_____________________________________________________________________________
acc49af9 1132Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
3fe61b77 1133{
1134 //
1135 // Get COG position
1136 // Used for clusters with more than 3 pads - where LUT not applicable
1137 //
1138
1139 Double_t sum = signal[0]
1140 + signal[1]
1141 + signal[2]
1142 + signal[3]
1143 + signal[4];
1144
1145 Double_t res = (0.0 * (-signal[0] + signal[4])
1146 + (-signal[1] + signal[3])) / sum;
1147
1148 return res;
1149
1150}
1151
1152//_____________________________________________________________________________
1153Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
1154{
1155 //
1156 // Method to unfold neighbouring maxima.
1157 // The charge ratio on the overlapping pad is calculated
1158 // until there is no more change within the range given by eps.
1159 // The resulting ratio is then returned to the calling method.
1160 //
1161
1162 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1163 if (!calibration) {
3fe61b77 1164 AliError("No AliTRDcalibDB instance available\n");
1165 return kFALSE;
6d50f529 1166 }
3fe61b77 1167
1168 Int_t irc = 0;
1169 Int_t itStep = 0; // Count iteration steps
1170
1171 Double_t ratio = 0.5; // Start value for ratio
1172 Double_t prevRatio = 0.0; // Store previous ratio
1173
1174 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1175 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1176 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1177
1178 // Start the iteration
1179 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1180
1181 itStep++;
1182 prevRatio = ratio;
1183
1184 // Cluster position according to charge ratio
1185 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1186 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1187 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1188 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1189
1190 // Set cluster charge ratio
1191 irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
1192 Double_t ampLeft = padSignal[1] / newSignal[1];
1193 irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
1194 Double_t ampRight = padSignal[3] / newSignal[1];
1195
1196 // Apply pad response to parameters
1197 irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
1198 irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
1199
1200 // Calculate new overlapping ratio
1201 ratio = TMath::Min((Double_t) 1.0
1202 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1203
b43a3e17 1204 }
88719a08 1205
3fe61b77 1206 return ratio;
1207
1208}
88719a08 1209
3fe61b77 1210//_____________________________________________________________________________
625f5260 1211void AliTRDclusterizer::TailCancelation(AliTRDdataArrayS *digitsIn
3fe61b77 1212 , AliTRDdataArrayF *digitsOut
1213 , AliTRDSignalIndex *indexesIn
1214 , AliTRDSignalIndex *indexesOut
1215 , Int_t nTimeTotal
acc49af9 1216 , Float_t adcThreshold
3fe61b77 1217 , AliTRDCalROC *calGainFactorROC
1218 , Float_t calGainFactorDetValue)
1219{
1220 //
1221 // Applies the tail cancelation and gain factors:
1222 // Transform digitsIn to digitsOut
1223 //
1224
1225 Int_t iRow = 0;
1226 Int_t iCol = 0;
1227 Int_t iTime = 0;
1228
1229 AliTRDRecParam *recParam = AliTRDRecParam::Instance();
1230 if (!recParam)
1231 {
1232 AliError("No AliTRDRecParam instance available\n");
1233 return;
1234 }
88719a08 1235
3fe61b77 1236 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
1237 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
1238 indexesIn->ResetCounters();
1239 while (indexesIn->NextRCIndex(iRow, iCol))
1240 {
1241 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1242 Double_t gain = calGainFactorDetValue
1243 * calGainFactorROCValue;
1244
1245 for (iTime = 0; iTime < nTimeTotal; iTime++)
1246 {
1247 // Apply gain gain factor
1248 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1249 inADC[iTime] /= gain;
1250 outADC[iTime] = inADC[iTime];
1251 }
1252
1253 // Apply the tail cancelation via the digital filter
1254 if (recParam->TCOn())
1255 {
1256 DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp());
1257 }
1258
1259 indexesIn->ResetTbinCounter();
1260 while (indexesIn->NextTbinIndex(iTime))
1261 {
1262 // Store the amplitude of the digit if above threshold
acc49af9 1263 if (outADC[iTime] > adcThreshold)
3fe61b77 1264 {
1265 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1266 indexesOut->AddIndexTBin(iRow,iCol,iTime);
1267 }
1268 } // while itime
1269
1270 } // while irow icol
1271
1272 delete [] inADC;
1273 delete [] outADC;
1274
1275 return;
1276
1277}
1278
1279//_____________________________________________________________________________
1280void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1281 , Int_t n, Int_t nexp)
1282{
1283 //
1284 // Tail cancellation by deconvolution for PASA v4 TRF
1285 //
1286
1287 Double_t rates[2];
1288 Double_t coefficients[2];
1289
1290 // Initialization (coefficient = alpha, rates = lambda)
1291 Double_t r1 = 1.0;
1292 Double_t r2 = 1.0;
1293 Double_t c1 = 0.5;
1294 Double_t c2 = 0.5;
1295
1296 if (nexp == 1) { // 1 Exponentials
1297 r1 = 1.156;
1298 r2 = 0.130;
1299 c1 = 0.066;
1300 c2 = 0.000;
1301 }
1302 if (nexp == 2) { // 2 Exponentials
1303 r1 = 1.156;
1304 r2 = 0.130;
1305 c1 = 0.114;
1306 c2 = 0.624;
1307 }
1308
1309 coefficients[0] = c1;
1310 coefficients[1] = c2;
1311
1312 Double_t dt = 0.1;
1313
1314 rates[0] = TMath::Exp(-dt/(r1));
1315 rates[1] = TMath::Exp(-dt/(r2));
1316
1317 Int_t i = 0;
1318 Int_t k = 0;
1319
1320 Double_t reminder[2];
1321 Double_t correction = 0.0;
1322 Double_t result = 0.0;
1323
1324 // Attention: computation order is important
1325 for (k = 0; k < nexp; k++) {
1326 reminder[k] = 0.0;
1327 }
1328
1329 for (i = 0; i < n; i++) {
1330
1331 result = (source[i] - correction); // No rescaling
1332 target[i] = result;
1333
1334 for (k = 0; k < nexp; k++) {
1335 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1336 }
1337
1338 correction = 0.0;
1339 for (k = 0; k < nexp; k++) {
1340 correction += reminder[k];
1341 }
1342
1343 }
6d50f529 1344
1345}
1346
1347//_____________________________________________________________________________
1348void AliTRDclusterizer::ResetRecPoints()
1349{
1350 //
1351 // Resets the list of rec points
1352 //
1353
1354 if (fRecPoints) {
1355 fRecPoints->Delete();
1356 }
1357
1358}
1359
1360//_____________________________________________________________________________
34eaaa7e 1361TObjArray *AliTRDclusterizer::RecPoints()
6d50f529 1362{
1363 //
1364 // Returns the list of rec points
1365 //
1366
1367 if (!fRecPoints) {
1368 fRecPoints = new TObjArray(400);
1369 }
1370
1371 return fRecPoints;
1372
3551db50 1373}