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