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