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