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