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