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