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