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