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