]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterizer.cxx
Some fixes for the mode w/o auto branch loading
[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
66f6bfd9 24#include <TClonesArray.h>
6d50f529 25#include <TObjArray.h>
f7336fa3 26
88cb7938 27#include "AliRunLoader.h"
28#include "AliLoader.h"
cc26f39c 29#include "AliTreeLoader.h"
3fe61b77 30#include "AliAlignObj.h"
88cb7938 31
f7336fa3 32#include "AliTRDclusterizer.h"
3e1a3ad8 33#include "AliTRDcluster.h"
fc546d21 34#include "AliTRDReconstructor.h"
793ff80c 35#include "AliTRDgeometry.h"
b65e5048 36#include "AliTRDarrayDictionary.h"
37#include "AliTRDarrayADC.h"
3fe61b77 38#include "AliTRDdigitsManager.h"
e570c2a5 39#include "AliTRDdigitsParam.h"
3fe61b77 40#include "AliTRDrawData.h"
3551db50 41#include "AliTRDcalibDB.h"
3fe61b77 42#include "AliTRDtransform.h"
43#include "AliTRDSignalIndex.h"
8dcbd181 44#include "AliTRDrawStream.h"
3fe61b77 45#include "AliTRDfeeParam.h"
a5b99acd 46#include "AliTRDtrackletWord.h"
3fe61b77 47
29f95561 48#include "TTreeStream.h"
49
3fe61b77 50#include "Cal/AliTRDCalROC.h"
51#include "Cal/AliTRDCalDet.h"
0e09df31 52#include "Cal/AliTRDCalSingleChamberStatus.h"
f7336fa3 53
54ClassImp(AliTRDclusterizer)
55
56//_____________________________________________________________________________
486c2339 57AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
6d50f529 58 :TNamed()
3a039a31 59 ,fReconstructor(rec)
6d50f529 60 ,fRunLoader(NULL)
61 ,fClusterTree(NULL)
62 ,fRecPoints(NULL)
a5b99acd 63 ,fTracklets(NULL)
9c7c9ec1 64 ,fTrackletTree(NULL)
3d0c7d6d 65 ,fDigitsManager(new AliTRDdigitsManager())
9c7c9ec1 66 ,fTrackletContainer(NULL)
3fe61b77 67 ,fRawVersion(2)
dfbb4bb9 68 ,fTransform(new AliTRDtransform(0))
1b3a719f 69 ,fDigits(NULL)
064d808d 70 ,fIndexes(NULL)
064d808d 71 ,fMaxThresh(0)
1ab4da8c 72 ,fMaxThreshTest(0)
064d808d 73 ,fSigThresh(0)
74 ,fMinMaxCutSigma(0)
75 ,fMinLeftRightCutSigma(0)
76 ,fLayer(0)
77 ,fDet(0)
78 ,fVolid(0)
79 ,fColMax(0)
80 ,fTimeTotal(0)
81 ,fCalGainFactorROC(NULL)
82 ,fCalGainFactorDetValue(0)
83 ,fCalNoiseROC(NULL)
84 ,fCalNoiseDetValue(0)
1b3a719f 85 ,fCalPadStatusROC(NULL)
064d808d 86 ,fClusterROC(0)
87 ,firstClusterROC(0)
3d0c7d6d 88 ,fNoOfClusters(0)
615f0826 89 ,fBaseline(0)
8dbc94c7 90 ,fRawStream(NULL)
f7336fa3 91{
92 //
93 // AliTRDclusterizer default constructor
94 //
95
b72f4eaf 96 SetBit(kLabels, kTRUE);
8cbe253a 97 SetBit(knewDM, kFALSE);
3d0c7d6d 98
3fe61b77 99 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
100
3a039a31 101 // Initialize debug stream
102 if(fReconstructor){
a2fbb6ec 103 if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 1){
3a039a31 104 TDirectory *savedir = gDirectory;
29f95561 105 //fgGetDebugStream = new TTreeSRedirector("TRD.ClusterizerDebug.root");
3a039a31 106 savedir->cd();
107 }
108 }
eb52b657 109
f7336fa3 110}
111
112//_____________________________________________________________________________
a987273c 113AliTRDclusterizer::AliTRDclusterizer(const Text_t *name
114 , const Text_t *title
115 , const AliTRDReconstructor *const rec)
6d50f529 116 :TNamed(name,title)
3a039a31 117 ,fReconstructor(rec)
6d50f529 118 ,fRunLoader(NULL)
119 ,fClusterTree(NULL)
120 ,fRecPoints(NULL)
a5b99acd 121 ,fTracklets(NULL)
9c7c9ec1 122 ,fTrackletTree(NULL)
3d0c7d6d 123 ,fDigitsManager(new AliTRDdigitsManager())
9c7c9ec1 124 ,fTrackletContainer(NULL)
3fe61b77 125 ,fRawVersion(2)
3fe61b77 126 ,fTransform(new AliTRDtransform(0))
1b3a719f 127 ,fDigits(NULL)
064d808d 128 ,fIndexes(NULL)
064d808d 129 ,fMaxThresh(0)
1ab4da8c 130 ,fMaxThreshTest(0)
064d808d 131 ,fSigThresh(0)
132 ,fMinMaxCutSigma(0)
133 ,fMinLeftRightCutSigma(0)
134 ,fLayer(0)
135 ,fDet(0)
136 ,fVolid(0)
137 ,fColMax(0)
138 ,fTimeTotal(0)
139 ,fCalGainFactorROC(NULL)
140 ,fCalGainFactorDetValue(0)
141 ,fCalNoiseROC(NULL)
142 ,fCalNoiseDetValue(0)
1b3a719f 143 ,fCalPadStatusROC(NULL)
064d808d 144 ,fClusterROC(0)
145 ,firstClusterROC(0)
3d0c7d6d 146 ,fNoOfClusters(0)
615f0826 147 ,fBaseline(0)
8dbc94c7 148 ,fRawStream(NULL)
f7336fa3 149{
150 //
6d50f529 151 // AliTRDclusterizer constructor
f7336fa3 152 //
153
b72f4eaf 154 SetBit(kLabels, kTRUE);
8cbe253a 155 SetBit(knewDM, kFALSE);
3d0c7d6d 156
3fe61b77 157 fDigitsManager->CreateArrays();
158
159 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
160
b72f4eaf 161 //FillLUT();
023b669c 162
f7336fa3 163}
164
8230f242 165//_____________________________________________________________________________
6d50f529 166AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
167 :TNamed(c)
3a039a31 168 ,fReconstructor(c.fReconstructor)
6d50f529 169 ,fRunLoader(NULL)
170 ,fClusterTree(NULL)
171 ,fRecPoints(NULL)
a5b99acd 172 ,fTracklets(NULL)
9c7c9ec1 173 ,fTrackletTree(NULL)
3fe61b77 174 ,fDigitsManager(NULL)
66f6bfd9 175 ,fTrackletContainer(NULL)
3fe61b77 176 ,fRawVersion(2)
3fe61b77 177 ,fTransform(NULL)
1b3a719f 178 ,fDigits(NULL)
064d808d 179 ,fIndexes(NULL)
064d808d 180 ,fMaxThresh(0)
1ab4da8c 181 ,fMaxThreshTest(0)
064d808d 182 ,fSigThresh(0)
183 ,fMinMaxCutSigma(0)
184 ,fMinLeftRightCutSigma(0)
185 ,fLayer(0)
186 ,fDet(0)
187 ,fVolid(0)
188 ,fColMax(0)
189 ,fTimeTotal(0)
190 ,fCalGainFactorROC(NULL)
191 ,fCalGainFactorDetValue(0)
192 ,fCalNoiseROC(NULL)
193 ,fCalNoiseDetValue(0)
1b3a719f 194 ,fCalPadStatusROC(NULL)
064d808d 195 ,fClusterROC(0)
196 ,firstClusterROC(0)
3d0c7d6d 197 ,fNoOfClusters(0)
615f0826 198 ,fBaseline(0)
8dbc94c7 199 ,fRawStream(NULL)
8230f242 200{
201 //
202 // AliTRDclusterizer copy constructor
203 //
204
b72f4eaf 205 SetBit(kLabels, kTRUE);
8cbe253a 206 SetBit(knewDM, kFALSE);
3d0c7d6d 207
b72f4eaf 208 //FillLUT();
fc546d21 209
8230f242 210}
211
f7336fa3 212//_____________________________________________________________________________
213AliTRDclusterizer::~AliTRDclusterizer()
214{
8230f242 215 //
216 // AliTRDclusterizer destructor
217 //
f7336fa3 218
48f8adf3 219 if (fRecPoints/* && IsClustersOwner()*/){
66f6bfd9 220 fRecPoints->Delete();
221 delete fRecPoints;
222 }
3fe61b77 223
a5b99acd 224 if (fTracklets){
225 fTracklets->Delete();
226 delete fTracklets;
227 }
228
66f6bfd9 229 if (fDigitsManager) {
230 delete fDigitsManager;
231 fDigitsManager = NULL;
232 }
9c7c9ec1 233
66f6bfd9 234 if (fTrackletContainer){
8cbe253a 235 delete [] fTrackletContainer[0];
236 delete [] fTrackletContainer[1];
237 delete [] fTrackletContainer;
66f6bfd9 238 fTrackletContainer = NULL;
239 }
3fe61b77 240
66f6bfd9 241 if (fTransform){
242 delete fTransform;
8dbc94c7 243 fTransform = NULL;
244 }
245
246 if (fRawStream){
247 delete fRawStream;
248 fRawStream = NULL;
66f6bfd9 249 }
fc546d21 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//_____________________________________________________________________________
0ab5ff49 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).fRawVersion = fRawVersion;
3fe61b77 282 ((AliTRDclusterizer &) c).fTransform = NULL;
1b3a719f 283 ((AliTRDclusterizer &) c).fDigits = NULL;
064d808d 284 ((AliTRDclusterizer &) c).fIndexes = NULL;
064d808d 285 ((AliTRDclusterizer &) c).fMaxThresh = 0;
1ab4da8c 286 ((AliTRDclusterizer &) c).fMaxThreshTest = 0;
064d808d 287 ((AliTRDclusterizer &) c).fSigThresh = 0;
288 ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0;
289 ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0;
290 ((AliTRDclusterizer &) c).fLayer = 0;
291 ((AliTRDclusterizer &) c).fDet = 0;
292 ((AliTRDclusterizer &) c).fVolid = 0;
293 ((AliTRDclusterizer &) c).fColMax = 0;
294 ((AliTRDclusterizer &) c).fTimeTotal = 0;
295 ((AliTRDclusterizer &) c).fCalGainFactorROC = NULL;
296 ((AliTRDclusterizer &) c).fCalGainFactorDetValue = 0;
297 ((AliTRDclusterizer &) c).fCalNoiseROC = NULL;
298 ((AliTRDclusterizer &) c).fCalNoiseDetValue = 0;
1b3a719f 299 ((AliTRDclusterizer &) c).fCalPadStatusROC = NULL;
064d808d 300 ((AliTRDclusterizer &) c).fClusterROC = 0;
301 ((AliTRDclusterizer &) c).firstClusterROC= 0;
3d0c7d6d 302 ((AliTRDclusterizer &) c).fNoOfClusters = 0;
615f0826 303 ((AliTRDclusterizer &) c).fBaseline = 0;
8dbc94c7 304 ((AliTRDclusterizer &) c).fRawStream = NULL;
615f0826 305
8230f242 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//_____________________________________________________________________________
e7295a3a 356Bool_t AliTRDclusterizer::OpenOutput(TTree *const clusterTree)
25ca55ce 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 }
11d0be11 368 return kTRUE;
369}
370
3e1a3ad8 371//_____________________________________________________________________________
88cb7938 372Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
f7336fa3 373{
374 //
375 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
376 //
377
f7336fa3 378 // Import the Trees for the event nEvent in the file
bdbb05bb 379 fRunLoader->GetEvent(nEvent);
88cb7938 380
f7336fa3 381 return kTRUE;
382
383}
384
385//_____________________________________________________________________________
793ff80c 386Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
f7336fa3 387{
388 //
3e1a3ad8 389 // Fills TRDcluster branch in the tree with the clusters
793ff80c 390 // found in detector = det. For det=-1 writes the tree.
a3c76cdc 391 //
793ff80c 392
6d50f529 393 if ((det < -1) ||
394 (det >= AliTRDgeometry::Ndet())) {
395 AliError(Form("Unexpected detector index %d.\n",det));
3e1a3ad8 396 return kFALSE;
793ff80c 397 }
317b5644 398
66f6bfd9 399 TObjArray *ioArray = new TObjArray(400);
3e1a3ad8 400 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
401 if (!branch) {
97b70082 402 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
66f6bfd9 403 } else branch->SetAddress(&ioArray);
66f6bfd9 404
405 Int_t nRecPoints = RecPoints()->GetEntriesFast();
406 if(det >= 0){
3e1a3ad8 407 for (Int_t i = 0; i < nRecPoints; i++) {
bdbb05bb 408 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
66f6bfd9 409 if(det != c->GetDetector()) continue;
410 ioArray->AddLast(c);
793ff80c 411 }
3e1a3ad8 412 fClusterTree->Fill();
839f6b9b 413 ioArray->Clear();
66f6bfd9 414 } else {
839f6b9b 415 Int_t detOld = -1, nw(0);
66f6bfd9 416 for (Int_t i = 0; i < nRecPoints; i++) {
417 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
418 if(c->GetDetector() != detOld){
839f6b9b 419 nw += ioArray->GetEntriesFast();
66f6bfd9 420 fClusterTree->Fill();
421 ioArray->Clear();
422 detOld = c->GetDetector();
423 }
424 ioArray->AddLast(c);
a6dd11e9 425 }
839f6b9b 426 if(ioArray->GetEntriesFast()){
427 nw += ioArray->GetEntriesFast();
428 fClusterTree->Fill();
429 ioArray->Clear();
430 }
431 AliDebug(2, Form("Clusters FOUND[%d] WRITTEN[%d] STATUS[%s]", nRecPoints, nw, nw==nRecPoints?"OK":"FAILED"));
793ff80c 432 }
66f6bfd9 433 delete ioArray;
6d50f529 434
66f6bfd9 435 return kTRUE;
f7336fa3 436}
793ff80c 437
9c7c9ec1 438//_____________________________________________________________________________
439Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
440{
eb52b657 441 //
442 // Write the raw data tracklets into seperate file
443 //
9c7c9ec1 444
445 UInt_t **leaves = new UInt_t *[2];
446 for (Int_t i=0; i<2 ;i++){
317b5644 447 leaves[i] = new UInt_t[258];
448 leaves[i][0] = det; // det
449 leaves[i][1] = i; // side
450 memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
9c7c9ec1 451 }
452
9c7c9ec1 453 if (!fTrackletTree){
454 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
455 dl->MakeTree();
456 fTrackletTree = dl->Tree();
457 }
458
459 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
460 if (!trkbranch) {
461 trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
462 }
463
464 for (Int_t i=0; i<2; i++){
317b5644 465 if (leaves[i][2]>0) {
466 trkbranch->SetAddress(leaves[i]);
467 fTrackletTree->Fill();
468 }
9c7c9ec1 469 }
470
f2ab58d4 471// AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets"); //jkl: wrong here
472// dl->WriteData("OVERWRITE"); //jkl: wrong here
9c7c9ec1 473 //dl->Unload();
474 delete [] leaves;
475
eb52b657 476 return kTRUE;
477
9c7c9ec1 478}
479
3fe61b77 480//_____________________________________________________________________________
481Bool_t AliTRDclusterizer::ReadDigits()
482{
483 //
484 // Reads the digits arrays from the input aliroot file
485 //
486
487 if (!fRunLoader) {
488 AliError("No run loader available");
489 return kFALSE;
490 }
491
492 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
493 if (!loader->TreeD()) {
494 loader->LoadDigits();
495 }
496
497 // Read in the digit arrays
498 return (fDigitsManager->ReadDigits(loader->TreeD()));
499
500}
501
502//_____________________________________________________________________________
503Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
504{
505 //
506 // Reads the digits arrays from the input tree
507 //
508
509 // Read in the digit arrays
510 return (fDigitsManager->ReadDigits(digitsTree));
511
512}
513
514//_____________________________________________________________________________
515Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
516{
517 //
518 // Reads the digits arrays from the ddl file
519 //
520
521 AliTRDrawData raw;
522 fDigitsManager = raw.Raw2Digits(rawReader);
523
524 return kTRUE;
525
526}
527
528//_____________________________________________________________________________
529Bool_t AliTRDclusterizer::MakeClusters()
530{
531 //
532 // Creates clusters from digits
533 //
534
535 // Propagate info from the digits manager
b72f4eaf 536 if (TestBit(kLabels)){
537 SetBit(kLabels, fDigitsManager->UsesDictionaries());
66f6bfd9 538 }
539
3fe61b77 540 Bool_t fReturn = kTRUE;
66f6bfd9 541 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
542
b65e5048 543 AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
66f6bfd9 544 // This is to take care of switched off super modules
545 if (!digitsIn->HasData()) continue;
546 digitsIn->Expand();
0d64b05f 547 digitsIn->DeleteNegatives(); // Restore digits array to >=0 values
66f6bfd9 548 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
549 if (indexes->IsAllocated() == kFALSE){
550 fDigitsManager->BuildIndexes(i);
551 }
552
553 Bool_t fR = kFALSE;
554 if (indexes->HasEntry()){
b72f4eaf 555 if (TestBit(kLabels)){
66f6bfd9 556 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
b65e5048 557 AliTRDarrayDictionary *tracksIn = 0; //mod
558 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict); //mod
66f6bfd9 559 tracksIn->Expand();
3fe61b77 560 }
66f6bfd9 561 }
562 fR = MakeClusters(i);
563 fReturn = fR && fReturn;
3fe61b77 564 }
66f6bfd9 565
566 //if (fR == kFALSE){
567 // if(IsWritingClusters()) WriteClusters(i);
568 // ResetRecPoints();
569 //}
570
5127281e 571 // Clear arrays of this chamber, to prepare for next event
572 fDigitsManager->ClearArrays(i);
66f6bfd9 573 }
574
575 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
3fe61b77 576
66f6bfd9 577 AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
3fe61b77 578
66f6bfd9 579 return fReturn;
eb52b657 580
3fe61b77 581}
582
583//_____________________________________________________________________________
584Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
585{
586 //
587 // Creates clusters from raw data
588 //
589
fc546d21 590 return Raw2ClustersChamber(rawReader);
3fe61b77 591
3fe61b77 592}
593
594//_____________________________________________________________________________
595Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
596{
597 //
598 // Creates clusters from raw data
599 //
600
601 // Create the digits manager
66f6bfd9 602 if (!fDigitsManager){
8cbe253a 603 SetBit(knewDM, kTRUE);
3d0c7d6d 604 fDigitsManager = new AliTRDdigitsManager(kTRUE);
66f6bfd9 605 fDigitsManager->CreateArrays();
606 }
3fe61b77 607
b72f4eaf 608 fDigitsManager->SetUseDictionaries(TestBit(kLabels));
3fe61b77 609
f2ab58d4 610 // ----- preparing tracklet output -----
611 if (fReconstructor->IsWritingTracklets()) {
612 AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
613 if (!trklLoader) {
cc26f39c 614 //AliInfo("Could not get the tracklets data loader, adding it now!");
f2ab58d4 615 trklLoader = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
616 AliRunLoader::Instance()->GetLoader("TRDLoader")->AddDataLoader(trklLoader);
617 }
cc26f39c 618 AliTreeLoader *trklTreeLoader = dynamic_cast<AliTreeLoader*> (trklLoader->GetBaseLoader("tracklets-raw"));
619 if (!trklTreeLoader) {
620 trklTreeLoader = new AliTreeLoader("tracklets-raw", trklLoader);
621 trklLoader->AddBaseLoader(trklTreeLoader);
622 }
623 if (!trklTreeLoader->Tree())
624 trklTreeLoader->MakeTree();
f2ab58d4 625 }
626
317b5644 627 // tracklet container for raw tracklet writing
a5b99acd 628 if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
317b5644 629 // maximum tracklets for one HC
630 const Int_t kTrackletChmb=256;
631 fTrackletContainer = new UInt_t *[2];
632 fTrackletContainer[0] = new UInt_t[kTrackletChmb];
633 fTrackletContainer[1] = new UInt_t[kTrackletChmb];
f2ab58d4 634 memset(fTrackletContainer[0], 0, kTrackletChmb*sizeof(UInt_t)); //jkl
635 memset(fTrackletContainer[1], 0, kTrackletChmb*sizeof(UInt_t)); //jkl
317b5644 636 }
637
8dbc94c7 638 if(!fRawStream)
8df1f8f5 639 fRawStream = new AliTRDrawStream(rawReader);
8dbc94c7 640 else
641 fRawStream->SetReader(rawReader);
642
d03bfe2b 643 if(fReconstructor->IsHLT()){
8dcbd181 644 if(fRawStream->InheritsFrom(AliTRDrawStream::Class()))
645 ((AliTRDrawStream*)fRawStream)->DisableErrorStorage();
646 else{
647 fRawStream->SetSharedPadReadout(kFALSE);
648 fRawStream->SetNoErrorWarning();
649 }
17e0e535 650 }
3fe61b77 651
17e0e535 652 AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
3fe61b77 653
b14f2845 654 UInt_t det = 0;
655 while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) < AliTRDgeometry::kNdet){
dc4457f0 656 if (fDigitsManager->GetIndexes(det)->HasEntry()){
c6f7c6cb 657 MakeClusters(det);
dc4457f0 658 fDigitsManager->ClearArrays(det);
659 }
660
317b5644 661 if (!fReconstructor->IsWritingTracklets()) continue;
662 if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
9c7c9ec1 663 }
486c2339 664
f2ab58d4 665 if (fReconstructor->IsWritingTracklets()) {
666 if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
cc26f39c 667 if (trklLoader) {
668 if (AliTreeLoader *trklTreeLoader = (AliTreeLoader*) trklLoader->GetBaseLoader("tracklets-raw"))
669 trklTreeLoader->WriteData("OVERWRITE");
670 trklLoader->UnloadAll();
671 }
f2ab58d4 672 }
673 }
674
a5b99acd 675 if (fTrackletContainer){
317b5644 676 delete [] fTrackletContainer[0];
677 delete [] fTrackletContainer[1];
678 delete [] fTrackletContainer;
679 fTrackletContainer = NULL;
680 }
681
66f6bfd9 682 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
3fe61b77 683
8cbe253a 684 if(!TestBit(knewDM)){
685 delete fDigitsManager;
686 fDigitsManager = NULL;
ad808511 687 delete fRawStream;
688 fRawStream = NULL;
8cbe253a 689 }
dfbb4bb9 690
3d0c7d6d 691 AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
66f6bfd9 692 return kTRUE;
eb52b657 693
3fe61b77 694}
695
fa7427d0 696//_____________________________________________________________________________
697UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
698{
023b669c 699 //
700 // Check if a pad is masked
701 //
702
317b5644 703 UChar_t status = 0;
704
705 if(signal>0 && TESTBIT(signal, 10)){
706 CLRBIT(signal, 10);
707 for(int ibit=0; ibit<4; ibit++){
708 if(TESTBIT(signal, 11+ibit)){
709 SETBIT(status, ibit);
710 CLRBIT(signal, 11+ibit);
711 }
712 }
713 }
714 return status;
fa7427d0 715}
716
6d9e79cc 717//_____________________________________________________________________________
e7295a3a 718void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
6d9e79cc 719 //
720 // Set the pad status into out
721 // First three bits are needed for the position encoding
722 //
064d808d 723 out |= status << 3;
6d9e79cc 724}
725
726//_____________________________________________________________________________
064d808d 727UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
6d9e79cc 728 //
729 // return the staus encoding of the corrupted pad
730 //
731 return static_cast<UChar_t>(encoding >> 3);
732}
733
734//_____________________________________________________________________________
064d808d 735Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
6d9e79cc 736 //
737 // Return the position of the corruption
738 //
739 return encoding & 7;
740}
741
3fe61b77 742//_____________________________________________________________________________
743Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
744{
745 //
746 // Generates the cluster.
747 //
748
749 // Get the digits
615f0826 750 fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
a0446ff1 751 fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline(det);
f5375dcb 752
3fe61b77 753 // This is to take care of switched off super modules
b72f4eaf 754 if (!fDigits->HasData()) return kFALSE;
3fe61b77 755
064d808d 756 fIndexes = fDigitsManager->GetIndexes(det);
b72f4eaf 757 if (fIndexes->IsAllocated() == kFALSE) {
758 AliError("Indexes do not exist!");
759 return kFALSE;
760 }
064d808d 761
fc0882f3 762 AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance();
b72f4eaf 763 if (!calibration) {
764 AliFatal("No AliTRDcalibDB instance available\n");
765 return kFALSE;
766 }
3fe61b77 767
3a039a31 768 if (!fReconstructor){
769 AliError("Reconstructor not set\n");
770 return kFALSE;
771 }
c5f589b9 772
c6f7c6cb 773 const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
fc0882f3 774
1ab4da8c 775 fMaxThresh = (Short_t)recoParam->GetClusMaxThresh();
7a72d9ad 776 fMaxThreshTest = (Short_t)(recoParam->GetClusMaxThresh()/2+fBaseline);
1ab4da8c 777 fSigThresh = (Short_t)recoParam->GetClusSigThresh();
fc0882f3 778 fMinMaxCutSigma = recoParam->GetMinMaxCutSigma();
779 fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma();
8c499dbf 780 const Int_t iEveryNTB = recoParam->GetRecEveryNTB();
3fe61b77 781
064d808d 782 Int_t istack = fIndexes->GetStack();
783 fLayer = fIndexes->GetLayer();
784 Int_t isector = fIndexes->GetSM();
3fe61b77 785
786 // Start clustering in the chamber
787
064d808d 788 fDet = AliTRDgeometry::GetDetector(fLayer,istack,isector);
789 if (fDet != det) {
828c6f80 790 AliError("Strange Detector number Missmatch!");
b65e5048 791 return kFALSE;
792 }
3fe61b77 793
839f6b9b 794 AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
795
3fe61b77 796 // TRD space point transformation
797 fTransform->SetDetector(det);
798
064d808d 799 Int_t iGeoLayer = AliGeomManager::kTRD1 + fLayer;
053767a4 800 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
064d808d 801 fVolid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
3fe61b77 802
a5b99acd 803 if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
804 AddTrackletsToArray();
805
1b3a719f 806 fColMax = fDigits->GetNcol();
8b8d2999 807 fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
3fe61b77 808
69a850a0 809 // Check consistency between OCDB and raw data
4e6823c2 810 Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
d03bfe2b 811 if(fReconstructor->IsHLT()){
c6f7c6cb 812 if((nTimeOCDB > -1) && (fTimeTotal != nTimeOCDB)){
813 AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value"
814 ,fTimeTotal,nTimeOCDB));
815 }
816 }else{
817 if(nTimeOCDB == -1){
818 AliWarning("Undefined number of timebins in OCDB, using value from raw data.");
819 if(!fTimeTotal>0){
820 AliError("Number of timebins in raw data is negative, skipping chamber!");
821 return kFALSE;
822 }
823 }else if(nTimeOCDB == -2){
824 AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!");
825 return kFALSE;
826 }else if(fTimeTotal != nTimeOCDB){
827 AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber!"
828 ,fTimeTotal,nTimeOCDB));
829 return kFALSE;
830 }
69a850a0 831 }
832
3fe61b77 833 // Detector wise calibration object for the gain factors
064d808d 834 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
3fe61b77 835 // Calibration object with pad wise values for the gain factors
064d808d 836 fCalGainFactorROC = calibration->GetGainFactorROC(fDet);
3fe61b77 837 // Calibration value for chamber wise gain factor
064d808d 838 fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
0e09df31 839
df83a620 840 // Detector wise calibration object for the noise
064d808d 841 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
df83a620 842 // Calibration object with pad wise values for the noise
064d808d 843 fCalNoiseROC = calibration->GetNoiseROC(fDet);
df83a620 844 // Calibration value for chamber wise noise
064d808d 845 fCalNoiseDetValue = calNoiseDet->GetValue(fDet);
1b3a719f 846
847 // Calibration object with the pad status
848 fCalPadStatusROC = calibration->GetPadStatusROC(fDet);
6b0d4ad5 849
064d808d 850 firstClusterROC = -1;
851 fClusterROC = 0;
3fe61b77 852
c6f7c6cb 853 SetBit(kLUT, recoParam->UseLUT());
854 SetBit(kGAUS, recoParam->UseGAUS());
855
b72f4eaf 856 // Apply the gain and the tail cancelation via digital filter
5f5f179e 857 // Use the configuration from the DCS to find out whether online
858 // tail cancellation was applied
859 if(!calibration->HasOnlineTailCancellation()) TailCancelation(recoParam);
023b669c 860
486c2339 861 MaxStruct curr, last;
064d808d 862 Int_t nMaximas = 0, nCorrupted = 0;
f5375dcb 863
064d808d 864 // Here the clusterfining is happening
865
8c499dbf 866 for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
615f0826 867 while(fIndexes->NextRCIndex(curr.row, curr.col)){
1ab4da8c 868 if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
615f0826 869 if(last.row>-1){
634a1625 870 if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
b72f4eaf 871 CreateCluster(last);
872 }
615f0826 873 last=curr; curr.fivePad=kFALSE;
b72f4eaf 874 }
b72f4eaf 875 }
876 }
615f0826 877 if(last.row>-1) CreateCluster(last);
df83a620 878
fc0882f3 879 if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
a2fbb6ec 880 TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
064d808d 881 (*fDebugStream) << "MakeClusters"
b72f4eaf 882 << "Detector=" << det
883 << "NMaxima=" << nMaximas
884 << "NClusters=" << fClusterROC
885 << "NCorrupted=" << nCorrupted
886 << "\n";
df83a620 887 }
b72f4eaf 888 if (TestBit(kLabels)) AddLabels();
f5375dcb 889
064d808d 890 return kTRUE;
f5375dcb 891
064d808d 892}
3fe61b77 893
064d808d 894//_____________________________________________________________________________
1b3a719f 895Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
064d808d 896{
897 //
898 // Returns true if this row,col,time combination is a maximum.
899 // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
900 //
c96d03ba 901
615f0826 902 Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
903 Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
1ab4da8c 904 if(Signals[1] <= fMaxThresh) return kFALSE;
3fe61b77 905
634a1625 906 if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
f5375dcb 907
634a1625 908 Short_t noiseMiddleThresh = (Short_t)(fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row));
909 if (Signals[1] <= noiseMiddleThresh) return kFALSE;
3fe61b77 910
b72f4eaf 911 UChar_t status[3]={
615f0826 912 fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
913 ,fCalPadStatusROC->GetStatus(Max.col, Max.row)
914 ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
b72f4eaf 915 };
486c2339 916
c66dc420 917 Short_t signal(0);
918 if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
919 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
920 Signals[0] = (Short_t)((signal - fBaseline) / gain + 0.5f);
921 } else Signals[0] = 0;
922 if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
923 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
924 Signals[2] = (Short_t)((signal - fBaseline) / gain + 0.5f);
925 } else Signals[2] = 0;
c96d03ba 926
486c2339 927 if(!(status[0] | status[1] | status[2])) {//all pads are good
c96d03ba 928 if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) {
1ab4da8c 929 if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
5a8db426 930 if(Signals[0]<0)Signals[0]=0;
931 if(Signals[2]<0)Signals[2]=0;
1ab4da8c 932 Short_t noiseSumThresh = (Short_t)(fMinLeftRightCutSigma * fCalNoiseDetValue
933 * fCalNoiseROC->GetValue(Max.col, Max.row));
934 if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
b72f4eaf 935 padStatus = 0;
936 return kTRUE;
eb52b657 937 }
064d808d 938 }
b72f4eaf 939 } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
5a8db426 940 if(Signals[0]<0)Signals[0]=0;
941 if(Signals[2]<0)Signals[2]=0;
1ab4da8c 942 if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
c96d03ba 943 Signals[2]=0;
486c2339 944 SetPadStatus(status[2], padStatus);
945 return kTRUE;
946 }
1ab4da8c 947 else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
c96d03ba 948 Signals[0]=0;
486c2339 949 SetPadStatus(status[0], padStatus);
950 return kTRUE;
951 }
1ab4da8c 952 else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
953 Signals[1] = fMaxThresh;
486c2339 954 SetPadStatus(status[1], padStatus);
064d808d 955 return kTRUE;
956 }
957 }
958 return kFALSE;
959}
3fe61b77 960
064d808d 961//_____________________________________________________________________________
1b3a719f 962Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
064d808d 963{
964 //
965 // Look for 5 pad cluster with minimum in the middle
966 // Gives back the ratio
967 //
615f0826 968
969 if (ThisMax.col >= fColMax - 3) return kFALSE;
970 Float_t gain;
971 if (ThisMax.col < fColMax - 5){
972 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
973 if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 974 return kFALSE;
064d808d 975 }
615f0826 976 if (ThisMax.col > 1) {
977 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
978 if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 979 return kFALSE;
980 }
981
486c2339 982 const Float_t kEpsilon = 0.01;
615f0826 983 Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
984 NeighbourMax.signals[1], NeighbourMax.signals[2]};
486c2339 985
986 // Unfold the two maxima and set the signal on
987 // the overlapping pad to the ratio
1b3a719f 988 Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
615f0826 989 ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
990 NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
991 ThisMax.fivePad=kTRUE;
992 NeighbourMax.fivePad=kTRUE;
486c2339 993 return kTRUE;
c96d03ba 994
064d808d 995}
eb52b657 996
064d808d 997//_____________________________________________________________________________
486c2339 998void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
064d808d 999{
1000 //
3d0c7d6d 1001 // Creates a cluster at the given position and saves it in fRecPoints
064d808d 1002 //
eb52b657 1003
c96d03ba 1004 Int_t nPadCount = 1;
615f0826 1005 Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
d03bfe2b 1006 if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
c96d03ba 1007
615f0826 1008 AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
c96d03ba 1009 cluster.SetNPads(nPadCount);
b72f4eaf 1010 if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1011 else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1012 else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
3fe61b77 1013
615f0826 1014 cluster.SetFivePad(Max.fivePad);
b72f4eaf 1015 // set pads status for the cluster
1016 UChar_t maskPosition = GetCorruption(Max.padStatus);
1017 if (maskPosition) {
1018 cluster.SetPadMaskedPosition(maskPosition);
1019 cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1020 }
8fc736d7 1021 cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
064d808d 1022
1023 // Transform the local cluster coordinates into calibrated
1024 // space point positions defined in the local tracking system.
1025 // Here the calibration for T0, Vdrift and ExB is applied as well.
d03bfe2b 1026 if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
82ddb093 1027
b72f4eaf 1028 // Temporarily store the Max.Row, column and time bin of the center pad
1029 // Used to later on assign the track indices
615f0826 1030 cluster.SetLabel(Max.row, 0);
1031 cluster.SetLabel(Max.col, 1);
1032 cluster.SetLabel(Max.time,2);
c96d03ba 1033
b72f4eaf 1034 //needed for HLT reconstruction
8c499dbf 1035 AddClusterToArray(&cluster);
b72f4eaf 1036
1037 // Store the index of the first cluster in the current ROC
1038 if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1039
1040 fNoOfClusters++;
1041 fClusterROC++;
3fe61b77 1042}
1043
3d0c7d6d 1044//_____________________________________________________________________________
1045void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1046{
c66dc420 1047// Calculate number of pads/cluster and
1048// ADC signals at position 0, 1, 5 and 6
1049
1050 Float_t gain(1.); Short_t signal(0.);
1051 // Store the amplitudes of the pads in the cluster for later analysis
1052 // and check whether one of these pads is masked in the database
1053 signals[3]=Max.signals[1];
1054 Int_t ipad(1), jpad(0);
3d0c7d6d 1055 // Look to the right
c66dc420 1056 while((jpad = Max.col-ipad)){
1057 if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1058 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1059 signal = (Short_t)((signal - fBaseline) / gain + 0.5f);
1060 if(signal<fSigThresh) break; // signal under threshold
3d0c7d6d 1061 nPadCount++;
c66dc420 1062 if(ipad<=3) signals[3 - ipad] = signal;
1063 ipad++;
3d0c7d6d 1064 }
c66dc420 1065 ipad=1;
3d0c7d6d 1066 // Look to the left
c66dc420 1067 while((jpad = Max.col+ipad)<fColMax){
1068 if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1069 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1070 signal = (Short_t)((signal - fBaseline) / gain + 0.5f);
1071 if(signal<fSigThresh) break; // signal under threshold
3d0c7d6d 1072 nPadCount++;
c66dc420 1073 if(ipad<=3) signals[3 + ipad] = signal;
1074 ipad++;
3d0c7d6d 1075 }
1076
c66dc420 1077 AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1078 , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
3d0c7d6d 1079}
1080
1081//_____________________________________________________________________________
e7295a3a 1082void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
3d0c7d6d 1083{
1084 //
1085 // Add a cluster to the array
1086 //
1087
1088 Int_t n = RecPoints()->GetEntriesFast();
1089 if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1090 new((*RecPoints())[n]) AliTRDcluster(*cluster);
1091}
1092
a5b99acd 1093//_____________________________________________________________________________
1094void AliTRDclusterizer::AddTrackletsToArray()
1095{
1096 //
1097 // Add the online tracklets of this chamber to the array
1098 //
1099
1100 UInt_t* trackletword;
1101 for(Int_t side=0; side<2; side++)
1102 {
1103 Int_t trkl=0;
1104 trackletword=fTrackletContainer[side];
97be9934 1105 while(trackletword[trkl]>0){
a5b99acd 1106 Int_t n = TrackletsArray()->GetEntriesFast();
35943b1e 1107 AliTRDtrackletWord tmp(trackletword[trkl]);
1108 new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
a5b99acd 1109 trkl++;
97be9934 1110 }
a5b99acd 1111 }
1112}
1113
3fe61b77 1114//_____________________________________________________________________________
6cd9a21f 1115Bool_t AliTRDclusterizer::AddLabels()
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
6cd9a21f 1133 Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
3fe61b77 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!
6cd9a21f 1141 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
3fe61b77 1142
1143 // Loop though the clusters found in this ROC
6cd9a21f 1144 for (iClusterROC = 0; iClusterROC < fClusterROC; 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
6cd9a21f 1164 for (iClusterROC = 0; iClusterROC < fClusterROC; 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
3fe61b77 1182//_____________________________________________________________________________
e7295a3a 1183Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
3fe61b77 1184{
1185 //
1186 // Method to unfold neighbouring maxima.
1187 // The charge ratio on the overlapping pad is calculated
1188 // until there is no more change within the range given by eps.
1189 // The resulting ratio is then returned to the calling method.
1190 //
1191
1192 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1193 if (!calibration) {
3fe61b77 1194 AliError("No AliTRDcalibDB instance available\n");
1195 return kFALSE;
6d50f529 1196 }
3fe61b77 1197
1198 Int_t irc = 0;
1199 Int_t itStep = 0; // Count iteration steps
1200
1201 Double_t ratio = 0.5; // Start value for ratio
1202 Double_t prevRatio = 0.0; // Store previous ratio
1203
1204 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1205 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1206 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1207
1208 // Start the iteration
1209 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1210
1211 itStep++;
1212 prevRatio = ratio;
1213
1214 // Cluster position according to charge ratio
1215 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
064d808d 1216 / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
3fe61b77 1217 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1218 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1219
1220 // Set cluster charge ratio
064d808d 1221 irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
3fe61b77 1222 Double_t ampLeft = padSignal[1] / newSignal[1];
064d808d 1223 irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
3fe61b77 1224 Double_t ampRight = padSignal[3] / newSignal[1];
1225
1226 // Apply pad response to parameters
053767a4 1227 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1228 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
3fe61b77 1229
1230 // Calculate new overlapping ratio
1231 ratio = TMath::Min((Double_t) 1.0
1232 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1233
b43a3e17 1234 }
88719a08 1235
3fe61b77 1236 return ratio;
1237
1238}
88719a08 1239
3fe61b77 1240//_____________________________________________________________________________
fc0882f3 1241void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
3fe61b77 1242{
1243 //
fc0882f3 1244 // Applies the tail cancelation
3fe61b77 1245 //
1246
1247 Int_t iRow = 0;
1248 Int_t iCol = 0;
1249 Int_t iTime = 0;
1250
a2fbb6ec 1251 TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
fc0882f3 1252 Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1253 Int_t nexp = recoParam->GetTCnexp();
064d808d 1254 while(fIndexes->NextRCIndex(iRow, iCol))
3fe61b77 1255 {
664030bc 1256 // if corrupted then don't make the tail cancallation
1257 if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
5a8db426 1258
cb5522da 1259 if(debugStreaming){
1ab4da8c 1260 for (iTime = 0; iTime < fTimeTotal; iTime++)
1261 (*fDebugStream) << "TailCancellation"
1262 << "col=" << iCol
1263 << "row=" << iRow
1264 << "\n";
cb5522da 1265 }
5a8db426 1266
664030bc 1267 // Apply the tail cancelation via the digital filter
1ab4da8c 1268 DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
5a8db426 1269
3fe61b77 1270 } // while irow icol
1b3a719f 1271
3fe61b77 1272 return;
1273
1274}
1275
1276//_____________________________________________________________________________
1ab4da8c 1277void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
3fe61b77 1278{
1279 //
1280 // Tail cancellation by deconvolution for PASA v4 TRF
1281 //
1282
664030bc 1283 Float_t rates[2];
1284 Float_t coefficients[2];
3fe61b77 1285
1286 // Initialization (coefficient = alpha, rates = lambda)
664030bc 1287 Float_t r1 = 1.0;
1288 Float_t r2 = 1.0;
1289 Float_t c1 = 0.5;
1290 Float_t c2 = 0.5;
3fe61b77 1291
1292 if (nexp == 1) { // 1 Exponentials
1293 r1 = 1.156;
1294 r2 = 0.130;
1295 c1 = 0.066;
1296 c2 = 0.000;
1297 }
1298 if (nexp == 2) { // 2 Exponentials
181c7f7e 1299 Double_t par[4];
a2fbb6ec 1300 fReconstructor->GetRecoParam()->GetTCParams(par);
181c7f7e 1301 r1 = par[0];//1.156;
1302 r2 = par[1];//0.130;
1303 c1 = par[2];//0.114;
1304 c2 = par[3];//0.624;
3fe61b77 1305 }
1306
1307 coefficients[0] = c1;
1308 coefficients[1] = c2;
1309
fc0882f3 1310 Double_t dt = 0.1;
3fe61b77 1311
1312 rates[0] = TMath::Exp(-dt/(r1));
634a1625 1313 rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
3fe61b77 1314
634a1625 1315 Float_t reminder[2] = { .0, .0 };
664030bc 1316 Float_t correction = 0.0;
1317 Float_t result = 0.0;
3fe61b77 1318
634a1625 1319 for (int i = 0; i < nTime; i++) {
3fe61b77 1320
1ab4da8c 1321 result = arr[i] - correction - fBaseline; // No rescaling
1322 arr[i] = (Short_t)(result + fBaseline + 0.5f);
3fe61b77 1323
3fe61b77 1324 correction = 0.0;
634a1625 1325 for (int k = 0; k < 2; k++) {
1326 correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
3fe61b77 1327 }
1328
1329 }
6d50f529 1330
1331}
1332
1333//_____________________________________________________________________________
1334void AliTRDclusterizer::ResetRecPoints()
1335{
1336 //
1337 // Resets the list of rec points
1338 //
1339
1340 if (fRecPoints) {
da0a4b87 1341 fRecPoints->Clear();
1342 fNoOfClusters = 0;
1343 // delete fRecPoints;
6d50f529 1344 }
6d50f529 1345}
1346
1347//_____________________________________________________________________________
66f6bfd9 1348TClonesArray *AliTRDclusterizer::RecPoints()
6d50f529 1349{
1350 //
1351 // Returns the list of rec points
1352 //
1353
1354 if (!fRecPoints) {
48f8adf3 1355 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
8ae98148 1356 // determine number of clusters which has to be allocated
1357 Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
8ae98148 1358
1359 fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
48f8adf3 1360 }
1361 //SetClustersOwner(kTRUE);
1362 AliTRDReconstructor::SetClusters(0x0);
6d50f529 1363 }
6d50f529 1364 return fRecPoints;
eb52b657 1365
3551db50 1366}
fc546d21 1367
a5b99acd 1368//_____________________________________________________________________________
1369TClonesArray *AliTRDclusterizer::TrackletsArray()
1370{
1371 //
1372 // Returns the list of rec points
1373 //
1374
1375 if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
94ac01ef 1376 fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
a5b99acd 1377 //SetClustersOwner(kTRUE);
1378 //AliTRDReconstructor::SetTracklets(0x0);
1379 }
1380 return fTracklets;
1381
1382}
1383