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