]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterizer.cxx
extra cut for clusters
[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()){
92305359 644 ((AliTRDrawStream*)fRawStream)->DisableErrorStorage();
17e0e535 645 }
3fe61b77 646
17e0e535 647 AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
3fe61b77 648
b14f2845 649 UInt_t det = 0;
650 while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) < AliTRDgeometry::kNdet){
dc4457f0 651 if (fDigitsManager->GetIndexes(det)->HasEntry()){
c6f7c6cb 652 MakeClusters(det);
dc4457f0 653 fDigitsManager->ClearArrays(det);
654 }
655
317b5644 656 if (!fReconstructor->IsWritingTracklets()) continue;
657 if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
9c7c9ec1 658 }
486c2339 659
f2ab58d4 660 if (fReconstructor->IsWritingTracklets()) {
661 if (AliDataLoader *trklLoader = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets")) {
cc26f39c 662 if (trklLoader) {
663 if (AliTreeLoader *trklTreeLoader = (AliTreeLoader*) trklLoader->GetBaseLoader("tracklets-raw"))
664 trklTreeLoader->WriteData("OVERWRITE");
665 trklLoader->UnloadAll();
666 }
f2ab58d4 667 }
668 }
669
a5b99acd 670 if (fTrackletContainer){
317b5644 671 delete [] fTrackletContainer[0];
672 delete [] fTrackletContainer[1];
673 delete [] fTrackletContainer;
674 fTrackletContainer = NULL;
675 }
676
66f6bfd9 677 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
3fe61b77 678
8cbe253a 679 if(!TestBit(knewDM)){
680 delete fDigitsManager;
681 fDigitsManager = NULL;
ad808511 682 delete fRawStream;
683 fRawStream = NULL;
8cbe253a 684 }
dfbb4bb9 685
3d0c7d6d 686 AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
66f6bfd9 687 return kTRUE;
eb52b657 688
3fe61b77 689}
690
fa7427d0 691//_____________________________________________________________________________
692UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
693{
023b669c 694 //
695 // Check if a pad is masked
696 //
697
317b5644 698 UChar_t status = 0;
699
700 if(signal>0 && TESTBIT(signal, 10)){
701 CLRBIT(signal, 10);
702 for(int ibit=0; ibit<4; ibit++){
703 if(TESTBIT(signal, 11+ibit)){
704 SETBIT(status, ibit);
705 CLRBIT(signal, 11+ibit);
706 }
707 }
708 }
709 return status;
fa7427d0 710}
711
6d9e79cc 712//_____________________________________________________________________________
e7295a3a 713void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
6d9e79cc 714 //
715 // Set the pad status into out
716 // First three bits are needed for the position encoding
717 //
064d808d 718 out |= status << 3;
6d9e79cc 719}
720
721//_____________________________________________________________________________
064d808d 722UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
6d9e79cc 723 //
724 // return the staus encoding of the corrupted pad
725 //
726 return static_cast<UChar_t>(encoding >> 3);
727}
728
729//_____________________________________________________________________________
064d808d 730Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
6d9e79cc 731 //
732 // Return the position of the corruption
733 //
734 return encoding & 7;
735}
736
3fe61b77 737//_____________________________________________________________________________
738Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
739{
740 //
741 // Generates the cluster.
742 //
743
744 // Get the digits
615f0826 745 fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
a0446ff1 746 fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline(det);
f5375dcb 747
3fe61b77 748 // This is to take care of switched off super modules
b72f4eaf 749 if (!fDigits->HasData()) return kFALSE;
3fe61b77 750
064d808d 751 fIndexes = fDigitsManager->GetIndexes(det);
b72f4eaf 752 if (fIndexes->IsAllocated() == kFALSE) {
753 AliError("Indexes do not exist!");
754 return kFALSE;
755 }
064d808d 756
fc0882f3 757 AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance();
b72f4eaf 758 if (!calibration) {
759 AliFatal("No AliTRDcalibDB instance available\n");
760 return kFALSE;
761 }
3fe61b77 762
3a039a31 763 if (!fReconstructor){
764 AliError("Reconstructor not set\n");
765 return kFALSE;
766 }
c5f589b9 767
c6f7c6cb 768 const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
fc0882f3 769
1ab4da8c 770 fMaxThresh = (Short_t)recoParam->GetClusMaxThresh();
7a72d9ad 771 fMaxThreshTest = (Short_t)(recoParam->GetClusMaxThresh()/2+fBaseline);
1ab4da8c 772 fSigThresh = (Short_t)recoParam->GetClusSigThresh();
fc0882f3 773 fMinMaxCutSigma = recoParam->GetMinMaxCutSigma();
774 fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma();
8c499dbf 775 const Int_t iEveryNTB = recoParam->GetRecEveryNTB();
3fe61b77 776
064d808d 777 Int_t istack = fIndexes->GetStack();
778 fLayer = fIndexes->GetLayer();
779 Int_t isector = fIndexes->GetSM();
3fe61b77 780
781 // Start clustering in the chamber
782
064d808d 783 fDet = AliTRDgeometry::GetDetector(fLayer,istack,isector);
784 if (fDet != det) {
828c6f80 785 AliError("Strange Detector number Missmatch!");
b65e5048 786 return kFALSE;
787 }
3fe61b77 788
839f6b9b 789 AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
790
3fe61b77 791 // TRD space point transformation
792 fTransform->SetDetector(det);
793
064d808d 794 Int_t iGeoLayer = AliGeomManager::kTRD1 + fLayer;
053767a4 795 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
064d808d 796 fVolid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
3fe61b77 797
a5b99acd 798 if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
799 AddTrackletsToArray();
800
1b3a719f 801 fColMax = fDigits->GetNcol();
8b8d2999 802 fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
3fe61b77 803
69a850a0 804 // Check consistency between OCDB and raw data
4e6823c2 805 Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
d03bfe2b 806 if(fReconstructor->IsHLT()){
c6f7c6cb 807 if((nTimeOCDB > -1) && (fTimeTotal != nTimeOCDB)){
808 AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value"
809 ,fTimeTotal,nTimeOCDB));
810 }
811 }else{
812 if(nTimeOCDB == -1){
813 AliWarning("Undefined number of timebins in OCDB, using value from raw data.");
814 if(!fTimeTotal>0){
815 AliError("Number of timebins in raw data is negative, skipping chamber!");
816 return kFALSE;
817 }
818 }else if(nTimeOCDB == -2){
819 AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!");
820 return kFALSE;
821 }else if(fTimeTotal != nTimeOCDB){
822 AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber!"
823 ,fTimeTotal,nTimeOCDB));
824 return kFALSE;
825 }
69a850a0 826 }
827
3fe61b77 828 // Detector wise calibration object for the gain factors
064d808d 829 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
3fe61b77 830 // Calibration object with pad wise values for the gain factors
064d808d 831 fCalGainFactorROC = calibration->GetGainFactorROC(fDet);
3fe61b77 832 // Calibration value for chamber wise gain factor
064d808d 833 fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
0e09df31 834
df83a620 835 // Detector wise calibration object for the noise
064d808d 836 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
df83a620 837 // Calibration object with pad wise values for the noise
064d808d 838 fCalNoiseROC = calibration->GetNoiseROC(fDet);
df83a620 839 // Calibration value for chamber wise noise
064d808d 840 fCalNoiseDetValue = calNoiseDet->GetValue(fDet);
1b3a719f 841
842 // Calibration object with the pad status
843 fCalPadStatusROC = calibration->GetPadStatusROC(fDet);
6b0d4ad5 844
064d808d 845 firstClusterROC = -1;
846 fClusterROC = 0;
3fe61b77 847
c6f7c6cb 848 SetBit(kLUT, recoParam->UseLUT());
849 SetBit(kGAUS, recoParam->UseGAUS());
850
b72f4eaf 851 // Apply the gain and the tail cancelation via digital filter
5f5f179e 852 // Use the configuration from the DCS to find out whether online
853 // tail cancellation was applied
854 if(!calibration->HasOnlineTailCancellation()) TailCancelation(recoParam);
023b669c 855
486c2339 856 MaxStruct curr, last;
064d808d 857 Int_t nMaximas = 0, nCorrupted = 0;
f5375dcb 858
064d808d 859 // Here the clusterfining is happening
860
8c499dbf 861 for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
615f0826 862 while(fIndexes->NextRCIndex(curr.row, curr.col)){
1ab4da8c 863 if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
615f0826 864 if(last.row>-1){
634a1625 865 if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
b72f4eaf 866 CreateCluster(last);
867 }
615f0826 868 last=curr; curr.fivePad=kFALSE;
b72f4eaf 869 }
b72f4eaf 870 }
871 }
615f0826 872 if(last.row>-1) CreateCluster(last);
df83a620 873
fc0882f3 874 if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
a2fbb6ec 875 TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
064d808d 876 (*fDebugStream) << "MakeClusters"
b72f4eaf 877 << "Detector=" << det
878 << "NMaxima=" << nMaximas
879 << "NClusters=" << fClusterROC
880 << "NCorrupted=" << nCorrupted
881 << "\n";
df83a620 882 }
b72f4eaf 883 if (TestBit(kLabels)) AddLabels();
f5375dcb 884
064d808d 885 return kTRUE;
f5375dcb 886
064d808d 887}
3fe61b77 888
064d808d 889//_____________________________________________________________________________
1b3a719f 890Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
064d808d 891{
892 //
893 // Returns true if this row,col,time combination is a maximum.
894 // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
895 //
c96d03ba 896
f05239f1 897 Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
615f0826 898 Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
f05239f1 899 tmp = (fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f;
900 Signals[1] = (Short_t)TMath::Min(tmp, kMaxShortVal);
1ab4da8c 901 if(Signals[1] <= fMaxThresh) return kFALSE;
3fe61b77 902
634a1625 903 if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
f5375dcb 904
634a1625 905 Short_t noiseMiddleThresh = (Short_t)(fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row));
906 if (Signals[1] <= noiseMiddleThresh) return kFALSE;
3fe61b77 907
b72f4eaf 908 UChar_t status[3]={
615f0826 909 fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
910 ,fCalPadStatusROC->GetStatus(Max.col, Max.row)
911 ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
b72f4eaf 912 };
486c2339 913
c66dc420 914 Short_t signal(0);
915 if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
916 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
f05239f1 917 tmp = (signal - fBaseline) / gain + 0.5f;
918 Signals[0] = (Short_t)TMath::Min(tmp, kMaxShortVal);
c66dc420 919 } else Signals[0] = 0;
920 if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
921 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
f05239f1 922 tmp = (signal - fBaseline) / gain + 0.5f;
923 Signals[2] = (Short_t)TMath::Min(tmp, kMaxShortVal);
c66dc420 924 } else Signals[2] = 0;
c96d03ba 925
486c2339 926 if(!(status[0] | status[1] | status[2])) {//all pads are good
c96d03ba 927 if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) {
1ab4da8c 928 if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
5a8db426 929 if(Signals[0]<0)Signals[0]=0;
930 if(Signals[2]<0)Signals[2]=0;
1ab4da8c 931 Short_t noiseSumThresh = (Short_t)(fMinLeftRightCutSigma * fCalNoiseDetValue
932 * fCalNoiseROC->GetValue(Max.col, Max.row));
933 if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
b72f4eaf 934 padStatus = 0;
935 return kTRUE;
eb52b657 936 }
064d808d 937 }
b72f4eaf 938 } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
5a8db426 939 if(Signals[0]<0)Signals[0]=0;
940 if(Signals[2]<0)Signals[2]=0;
1ab4da8c 941 if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
c96d03ba 942 Signals[2]=0;
486c2339 943 SetPadStatus(status[2], padStatus);
944 return kTRUE;
945 }
1ab4da8c 946 else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
c96d03ba 947 Signals[0]=0;
486c2339 948 SetPadStatus(status[0], padStatus);
949 return kTRUE;
950 }
1ab4da8c 951 else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
952 Signals[1] = fMaxThresh;
486c2339 953 SetPadStatus(status[1], padStatus);
064d808d 954 return kTRUE;
955 }
956 }
957 return kFALSE;
958}
3fe61b77 959
064d808d 960//_____________________________________________________________________________
1b3a719f 961Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
064d808d 962{
963 //
964 // Look for 5 pad cluster with minimum in the middle
965 // Gives back the ratio
966 //
615f0826 967
968 if (ThisMax.col >= fColMax - 3) return kFALSE;
969 Float_t gain;
970 if (ThisMax.col < fColMax - 5){
971 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
972 if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 973 return kFALSE;
064d808d 974 }
615f0826 975 if (ThisMax.col > 1) {
976 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
977 if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 978 return kFALSE;
979 }
980
486c2339 981 const Float_t kEpsilon = 0.01;
615f0826 982 Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
983 NeighbourMax.signals[1], NeighbourMax.signals[2]};
486c2339 984
985 // Unfold the two maxima and set the signal on
986 // the overlapping pad to the ratio
1b3a719f 987 Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
615f0826 988 ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
989 NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
990 ThisMax.fivePad=kTRUE;
991 NeighbourMax.fivePad=kTRUE;
486c2339 992 return kTRUE;
c96d03ba 993
064d808d 994}
eb52b657 995
064d808d 996//_____________________________________________________________________________
486c2339 997void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
064d808d 998{
999 //
3d0c7d6d 1000 // Creates a cluster at the given position and saves it in fRecPoints
064d808d 1001 //
eb52b657 1002
c96d03ba 1003 Int_t nPadCount = 1;
615f0826 1004 Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
d03bfe2b 1005 if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
c96d03ba 1006
615f0826 1007 AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
c96d03ba 1008 cluster.SetNPads(nPadCount);
b72f4eaf 1009 if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1010 else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1011 else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
3fe61b77 1012
615f0826 1013 cluster.SetFivePad(Max.fivePad);
b72f4eaf 1014 // set pads status for the cluster
1015 UChar_t maskPosition = GetCorruption(Max.padStatus);
1016 if (maskPosition) {
1017 cluster.SetPadMaskedPosition(maskPosition);
1018 cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1019 }
8fc736d7 1020 cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
064d808d 1021
1022 // Transform the local cluster coordinates into calibrated
1023 // space point positions defined in the local tracking system.
1024 // Here the calibration for T0, Vdrift and ExB is applied as well.
d03bfe2b 1025 if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
82ddb093 1026
b72f4eaf 1027 // Temporarily store the Max.Row, column and time bin of the center pad
1028 // Used to later on assign the track indices
615f0826 1029 cluster.SetLabel(Max.row, 0);
1030 cluster.SetLabel(Max.col, 1);
1031 cluster.SetLabel(Max.time,2);
c96d03ba 1032
b72f4eaf 1033 //needed for HLT reconstruction
8c499dbf 1034 AddClusterToArray(&cluster);
b72f4eaf 1035
1036 // Store the index of the first cluster in the current ROC
1037 if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1038
1039 fNoOfClusters++;
1040 fClusterROC++;
3fe61b77 1041}
1042
3d0c7d6d 1043//_____________________________________________________________________________
1044void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1045{
c66dc420 1046// Calculate number of pads/cluster and
1047// ADC signals at position 0, 1, 5 and 6
1048
f05239f1 1049 Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
1050 Float_t gain(1.); Short_t signal(0);
c66dc420 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);
f05239f1 1059 tmp = (signal - fBaseline) / gain + 0.5f;
1060 signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
c66dc420 1061 if(signal<fSigThresh) break; // signal under threshold
3d0c7d6d 1062 nPadCount++;
c66dc420 1063 if(ipad<=3) signals[3 - ipad] = signal;
1064 ipad++;
3d0c7d6d 1065 }
c66dc420 1066 ipad=1;
3d0c7d6d 1067 // Look to the left
c66dc420 1068 while((jpad = Max.col+ipad)<fColMax){
1069 if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1070 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
f05239f1 1071 tmp = (signal - fBaseline) / gain + 0.5f;
1072 signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
c66dc420 1073 if(signal<fSigThresh) break; // signal under threshold
3d0c7d6d 1074 nPadCount++;
c66dc420 1075 if(ipad<=3) signals[3 + ipad] = signal;
1076 ipad++;
3d0c7d6d 1077 }
1078
c66dc420 1079 AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1080 , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
3d0c7d6d 1081}
1082
1083//_____________________________________________________________________________
e7295a3a 1084void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
3d0c7d6d 1085{
1086 //
1087 // Add a cluster to the array
1088 //
1089
1090 Int_t n = RecPoints()->GetEntriesFast();
1091 if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1092 new((*RecPoints())[n]) AliTRDcluster(*cluster);
1093}
1094
a5b99acd 1095//_____________________________________________________________________________
1096void AliTRDclusterizer::AddTrackletsToArray()
1097{
1098 //
1099 // Add the online tracklets of this chamber to the array
1100 //
1101
1102 UInt_t* trackletword;
1103 for(Int_t side=0; side<2; side++)
1104 {
1105 Int_t trkl=0;
1106 trackletword=fTrackletContainer[side];
97be9934 1107 while(trackletword[trkl]>0){
a5b99acd 1108 Int_t n = TrackletsArray()->GetEntriesFast();
35943b1e 1109 AliTRDtrackletWord tmp(trackletword[trkl]);
1110 new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
a5b99acd 1111 trkl++;
97be9934 1112 }
a5b99acd 1113 }
1114}
1115
3fe61b77 1116//_____________________________________________________________________________
6cd9a21f 1117Bool_t AliTRDclusterizer::AddLabels()
3551db50 1118{
1119 //
3fe61b77 1120 // Add the track indices to the found clusters
3551db50 1121 //
1122
3fe61b77 1123 const Int_t kNclus = 3;
1124 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1125 const Int_t kNtrack = kNdict * kNclus;
1126
1127 Int_t iClusterROC = 0;
1128
1129 Int_t row = 0;
1130 Int_t col = 0;
1131 Int_t time = 0;
1132 Int_t iPad = 0;
1133
1134 // Temporary array to collect the track indices
6cd9a21f 1135 Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
3fe61b77 1136
1137 // Loop through the dictionary arrays one-by-one
1138 // to keep memory consumption low
b65e5048 1139 AliTRDarrayDictionary *tracksIn = 0; //mod
3fe61b77 1140 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1141
1142 // tracksIn should be expanded beforehand!
6cd9a21f 1143 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
3fe61b77 1144
1145 // Loop though the clusters found in this ROC
6cd9a21f 1146 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1147
3fe61b77 1148 AliTRDcluster *cluster = (AliTRDcluster *)
b65e5048 1149 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
3fe61b77 1150 row = cluster->GetLabel(0);
1151 col = cluster->GetLabel(1);
1152 time = cluster->GetLabel(2);
1153
1154 for (iPad = 0; iPad < kNclus; iPad++) {
b65e5048 1155 Int_t iPadCol = col - 1 + iPad;
1156 Int_t index = tracksIn->GetData(row,iPadCol,time); //Modification of -1 in Track
1157 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
3fe61b77 1158 }
1159
1160 }
1161
1162 }
1163
1164 // Copy the track indices into the cluster
1165 // Loop though the clusters found in this ROC
6cd9a21f 1166 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1167
3fe61b77 1168 AliTRDcluster *cluster = (AliTRDcluster *)
1169 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1170 cluster->SetLabel(-9999,0);
1171 cluster->SetLabel(-9999,1);
1172 cluster->SetLabel(-9999,2);
1173
1174 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1175
1176 }
1177
1178 delete [] idxTracks;
1179
1180 return kTRUE;
1181
1182}
1183
3fe61b77 1184//_____________________________________________________________________________
e7295a3a 1185Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
3fe61b77 1186{
1187 //
1188 // Method to unfold neighbouring maxima.
1189 // The charge ratio on the overlapping pad is calculated
1190 // until there is no more change within the range given by eps.
1191 // The resulting ratio is then returned to the calling method.
1192 //
1193
1194 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1195 if (!calibration) {
3fe61b77 1196 AliError("No AliTRDcalibDB instance available\n");
1197 return kFALSE;
6d50f529 1198 }
3fe61b77 1199
1200 Int_t irc = 0;
1201 Int_t itStep = 0; // Count iteration steps
1202
1203 Double_t ratio = 0.5; // Start value for ratio
1204 Double_t prevRatio = 0.0; // Store previous ratio
1205
1206 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1207 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1208 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1209
1210 // Start the iteration
1211 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1212
1213 itStep++;
1214 prevRatio = ratio;
1215
1216 // Cluster position according to charge ratio
1217 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
064d808d 1218 / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
3fe61b77 1219 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1220 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1221
1222 // Set cluster charge ratio
064d808d 1223 irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
3fe61b77 1224 Double_t ampLeft = padSignal[1] / newSignal[1];
064d808d 1225 irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
3fe61b77 1226 Double_t ampRight = padSignal[3] / newSignal[1];
1227
1228 // Apply pad response to parameters
053767a4 1229 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1230 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
3fe61b77 1231
1232 // Calculate new overlapping ratio
1233 ratio = TMath::Min((Double_t) 1.0
1234 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1235
b43a3e17 1236 }
88719a08 1237
3fe61b77 1238 return ratio;
1239
1240}
88719a08 1241
3fe61b77 1242//_____________________________________________________________________________
fc0882f3 1243void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
3fe61b77 1244{
1245 //
fc0882f3 1246 // Applies the tail cancelation
3fe61b77 1247 //
1248
f05239f1 1249 Int_t nexp = recoParam->GetTCnexp();
1250 if(!nexp) return;
1251
3fe61b77 1252 Int_t iRow = 0;
1253 Int_t iCol = 0;
1254 Int_t iTime = 0;
1255
a2fbb6ec 1256 TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
fc0882f3 1257 Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
064d808d 1258 while(fIndexes->NextRCIndex(iRow, iCol))
3fe61b77 1259 {
664030bc 1260 // if corrupted then don't make the tail cancallation
1261 if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
5a8db426 1262
cb5522da 1263 if(debugStreaming){
1ab4da8c 1264 for (iTime = 0; iTime < fTimeTotal; iTime++)
1265 (*fDebugStream) << "TailCancellation"
1266 << "col=" << iCol
1267 << "row=" << iRow
1268 << "\n";
cb5522da 1269 }
5a8db426 1270
664030bc 1271 // Apply the tail cancelation via the digital filter
1ab4da8c 1272 DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
5a8db426 1273
3fe61b77 1274 } // while irow icol
1b3a719f 1275
3fe61b77 1276 return;
1277
1278}
1279
1280//_____________________________________________________________________________
1ab4da8c 1281void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
3fe61b77 1282{
1283 //
1284 // Tail cancellation by deconvolution for PASA v4 TRF
1285 //
1286
664030bc 1287 Float_t rates[2];
1288 Float_t coefficients[2];
3fe61b77 1289
1290 // Initialization (coefficient = alpha, rates = lambda)
664030bc 1291 Float_t r1 = 1.0;
1292 Float_t r2 = 1.0;
1293 Float_t c1 = 0.5;
1294 Float_t c2 = 0.5;
3fe61b77 1295
1296 if (nexp == 1) { // 1 Exponentials
1297 r1 = 1.156;
1298 r2 = 0.130;
1299 c1 = 0.066;
1300 c2 = 0.000;
1301 }
1302 if (nexp == 2) { // 2 Exponentials
181c7f7e 1303 Double_t par[4];
a2fbb6ec 1304 fReconstructor->GetRecoParam()->GetTCParams(par);
181c7f7e 1305 r1 = par[0];//1.156;
1306 r2 = par[1];//0.130;
1307 c1 = par[2];//0.114;
1308 c2 = par[3];//0.624;
3fe61b77 1309 }
1310
1311 coefficients[0] = c1;
1312 coefficients[1] = c2;
1313
fc0882f3 1314 Double_t dt = 0.1;
3fe61b77 1315
1316 rates[0] = TMath::Exp(-dt/(r1));
634a1625 1317 rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
3fe61b77 1318
634a1625 1319 Float_t reminder[2] = { .0, .0 };
664030bc 1320 Float_t correction = 0.0;
1321 Float_t result = 0.0;
3fe61b77 1322
634a1625 1323 for (int i = 0; i < nTime; i++) {
3fe61b77 1324
1ab4da8c 1325 result = arr[i] - correction - fBaseline; // No rescaling
1326 arr[i] = (Short_t)(result + fBaseline + 0.5f);
3fe61b77 1327
3fe61b77 1328 correction = 0.0;
634a1625 1329 for (int k = 0; k < 2; k++) {
1330 correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
3fe61b77 1331 }
1332
1333 }
6d50f529 1334
1335}
1336
1337//_____________________________________________________________________________
1338void AliTRDclusterizer::ResetRecPoints()
1339{
1340 //
1341 // Resets the list of rec points
1342 //
1343
1344 if (fRecPoints) {
da0a4b87 1345 fRecPoints->Clear();
1346 fNoOfClusters = 0;
1347 // delete fRecPoints;
6d50f529 1348 }
6d50f529 1349}
1350
1351//_____________________________________________________________________________
66f6bfd9 1352TClonesArray *AliTRDclusterizer::RecPoints()
6d50f529 1353{
1354 //
1355 // Returns the list of rec points
1356 //
1357
1358 if (!fRecPoints) {
48f8adf3 1359 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
8ae98148 1360 // determine number of clusters which has to be allocated
1361 Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
8ae98148 1362
1363 fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
48f8adf3 1364 }
1365 //SetClustersOwner(kTRUE);
1366 AliTRDReconstructor::SetClusters(0x0);
6d50f529 1367 }
6d50f529 1368 return fRecPoints;
eb52b657 1369
3551db50 1370}
fc546d21 1371
a5b99acd 1372//_____________________________________________________________________________
1373TClonesArray *AliTRDclusterizer::TrackletsArray()
1374{
1375 //
1376 // Returns the list of rec points
1377 //
1378
1379 if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
94ac01ef 1380 fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
a5b99acd 1381 //SetClustersOwner(kTRUE);
1382 //AliTRDReconstructor::SetTracklets(0x0);
1383 }
1384 return fTracklets;
1385
1386}
1387