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