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