]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterizer.cxx
Some hick-up with common blocks solved.
[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//_____________________________________________________________________________
0ab5ff49 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();
839f6b9b 456 ioArray->Clear();
66f6bfd9 457 } else {
839f6b9b 458 Int_t detOld = -1, nw(0);
66f6bfd9 459 for (Int_t i = 0; i < nRecPoints; i++) {
460 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
461 if(c->GetDetector() != detOld){
839f6b9b 462 nw += ioArray->GetEntriesFast();
66f6bfd9 463 fClusterTree->Fill();
464 ioArray->Clear();
465 detOld = c->GetDetector();
466 }
467 ioArray->AddLast(c);
a6dd11e9 468 }
839f6b9b 469 if(ioArray->GetEntriesFast()){
470 nw += ioArray->GetEntriesFast();
471 fClusterTree->Fill();
472 ioArray->Clear();
473 }
474 AliDebug(2, Form("Clusters FOUND[%d] WRITTEN[%d] STATUS[%s]", nRecPoints, nw, nw==nRecPoints?"OK":"FAILED"));
793ff80c 475 }
66f6bfd9 476 delete ioArray;
6d50f529 477
66f6bfd9 478 return kTRUE;
f7336fa3 479}
793ff80c 480
9c7c9ec1 481//_____________________________________________________________________________
482Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
483{
eb52b657 484 //
485 // Write the raw data tracklets into seperate file
486 //
9c7c9ec1 487
488 UInt_t **leaves = new UInt_t *[2];
489 for (Int_t i=0; i<2 ;i++){
317b5644 490 leaves[i] = new UInt_t[258];
491 leaves[i][0] = det; // det
492 leaves[i][1] = i; // side
493 memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
9c7c9ec1 494 }
495
9c7c9ec1 496 if (!fTrackletTree){
497 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
498 dl->MakeTree();
499 fTrackletTree = dl->Tree();
500 }
501
502 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
503 if (!trkbranch) {
504 trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
505 }
506
507 for (Int_t i=0; i<2; i++){
317b5644 508 if (leaves[i][2]>0) {
509 trkbranch->SetAddress(leaves[i]);
510 fTrackletTree->Fill();
511 }
9c7c9ec1 512 }
513
514 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
515 dl->WriteData("OVERWRITE");
516 //dl->Unload();
517 delete [] leaves;
518
eb52b657 519 return kTRUE;
520
9c7c9ec1 521}
522
3fe61b77 523//_____________________________________________________________________________
524Bool_t AliTRDclusterizer::ReadDigits()
525{
526 //
527 // Reads the digits arrays from the input aliroot file
528 //
529
530 if (!fRunLoader) {
531 AliError("No run loader available");
532 return kFALSE;
533 }
534
535 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
536 if (!loader->TreeD()) {
537 loader->LoadDigits();
538 }
539
540 // Read in the digit arrays
541 return (fDigitsManager->ReadDigits(loader->TreeD()));
542
543}
544
545//_____________________________________________________________________________
546Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
547{
548 //
549 // Reads the digits arrays from the input tree
550 //
551
552 // Read in the digit arrays
553 return (fDigitsManager->ReadDigits(digitsTree));
554
555}
556
557//_____________________________________________________________________________
558Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
559{
560 //
561 // Reads the digits arrays from the ddl file
562 //
563
564 AliTRDrawData raw;
565 fDigitsManager = raw.Raw2Digits(rawReader);
566
567 return kTRUE;
568
569}
570
571//_____________________________________________________________________________
572Bool_t AliTRDclusterizer::MakeClusters()
573{
574 //
575 // Creates clusters from digits
576 //
577
578 // Propagate info from the digits manager
b72f4eaf 579 if (TestBit(kLabels)){
580 SetBit(kLabels, fDigitsManager->UsesDictionaries());
66f6bfd9 581 }
582
3fe61b77 583 Bool_t fReturn = kTRUE;
66f6bfd9 584 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
585
b65e5048 586 AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
66f6bfd9 587 // This is to take care of switched off super modules
588 if (!digitsIn->HasData()) continue;
589 digitsIn->Expand();
0d64b05f 590 digitsIn->DeleteNegatives(); // Restore digits array to >=0 values
66f6bfd9 591 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
592 if (indexes->IsAllocated() == kFALSE){
593 fDigitsManager->BuildIndexes(i);
594 }
595
596 Bool_t fR = kFALSE;
597 if (indexes->HasEntry()){
b72f4eaf 598 if (TestBit(kLabels)){
66f6bfd9 599 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
b65e5048 600 AliTRDarrayDictionary *tracksIn = 0; //mod
601 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict); //mod
66f6bfd9 602 tracksIn->Expand();
3fe61b77 603 }
66f6bfd9 604 }
605 fR = MakeClusters(i);
606 fReturn = fR && fReturn;
3fe61b77 607 }
66f6bfd9 608
609 //if (fR == kFALSE){
610 // if(IsWritingClusters()) WriteClusters(i);
611 // ResetRecPoints();
612 //}
613
614 // No compress just remove
615 fDigitsManager->RemoveDigits(i);
616 fDigitsManager->RemoveDictionaries(i);
617 fDigitsManager->ClearIndexes(i);
618 }
828c6f80 619 fReconstructor->SetDigitsParam(fDigitsManager->GetDigitsParam());
66f6bfd9 620
621 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
3fe61b77 622
66f6bfd9 623 AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
3fe61b77 624
66f6bfd9 625 return fReturn;
eb52b657 626
3fe61b77 627}
628
629//_____________________________________________________________________________
630Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
631{
632 //
633 // Creates clusters from raw data
634 //
635
fc546d21 636 return Raw2ClustersChamber(rawReader);
3fe61b77 637
3fe61b77 638}
639
640//_____________________________________________________________________________
641Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
642{
643 //
644 // Creates clusters from raw data
645 //
646
647 // Create the digits manager
66f6bfd9 648 if (!fDigitsManager){
8cbe253a 649 SetBit(knewDM, kTRUE);
3d0c7d6d 650 fDigitsManager = new AliTRDdigitsManager(kTRUE);
66f6bfd9 651 fDigitsManager->CreateArrays();
652 }
3fe61b77 653
b72f4eaf 654 fDigitsManager->SetUseDictionaries(TestBit(kLabels));
3fe61b77 655
317b5644 656 // tracklet container for raw tracklet writing
a5b99acd 657 if (!fTrackletContainer && ( fReconstructor->IsWritingTracklets() || fReconstructor->IsProcessingTracklets() )) {
317b5644 658 // maximum tracklets for one HC
659 const Int_t kTrackletChmb=256;
660 fTrackletContainer = new UInt_t *[2];
661 fTrackletContainer[0] = new UInt_t[kTrackletChmb];
662 fTrackletContainer[1] = new UInt_t[kTrackletChmb];
663 }
664
8dbc94c7 665 if(!fRawStream)
666 fRawStream = AliTRDrawStreamBase::GetRawStream(rawReader);
667 else
668 fRawStream->SetReader(rawReader);
669
17e0e535 670 if(fReconstructor->IsHLT()){
8dbc94c7 671 fRawStream->SetSharedPadReadout(kFALSE);
17e0e535 672 fRawStream->SetNoErrorWarning();
673 }
3fe61b77 674
17e0e535 675 AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
3fe61b77 676
677 Int_t det = 0;
8dbc94c7 678 while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
66f6bfd9 679 Bool_t iclusterBranch = kFALSE;
680 if (fDigitsManager->GetIndexes(det)->HasEntry()){
486c2339 681 iclusterBranch = MakeClusters(det);
3fe61b77 682 }
66f6bfd9 683
66f68968 684 fDigitsManager->ClearArrays(det);
685
317b5644 686 if (!fReconstructor->IsWritingTracklets()) continue;
687 if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
9c7c9ec1 688 }
828c6f80 689 fReconstructor->SetDigitsParam(fDigitsManager->GetDigitsParam());
486c2339 690
a5b99acd 691 if (fTrackletContainer){
317b5644 692 delete [] fTrackletContainer[0];
693 delete [] fTrackletContainer[1];
694 delete [] fTrackletContainer;
695 fTrackletContainer = NULL;
696 }
697
66f6bfd9 698 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
3fe61b77 699
8cbe253a 700 if(!TestBit(knewDM)){
701 delete fDigitsManager;
702 fDigitsManager = NULL;
ad808511 703 delete fRawStream;
704 fRawStream = NULL;
8cbe253a 705 }
dfbb4bb9 706
3d0c7d6d 707 AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
66f6bfd9 708 return kTRUE;
eb52b657 709
3fe61b77 710}
711
fa7427d0 712//_____________________________________________________________________________
713UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
714{
023b669c 715 //
716 // Check if a pad is masked
717 //
718
317b5644 719 UChar_t status = 0;
720
721 if(signal>0 && TESTBIT(signal, 10)){
722 CLRBIT(signal, 10);
723 for(int ibit=0; ibit<4; ibit++){
724 if(TESTBIT(signal, 11+ibit)){
725 SETBIT(status, ibit);
726 CLRBIT(signal, 11+ibit);
727 }
728 }
729 }
730 return status;
fa7427d0 731}
732
6d9e79cc 733//_____________________________________________________________________________
e7295a3a 734void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
6d9e79cc 735 //
736 // Set the pad status into out
737 // First three bits are needed for the position encoding
738 //
064d808d 739 out |= status << 3;
6d9e79cc 740}
741
742//_____________________________________________________________________________
064d808d 743UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
6d9e79cc 744 //
745 // return the staus encoding of the corrupted pad
746 //
747 return static_cast<UChar_t>(encoding >> 3);
748}
749
750//_____________________________________________________________________________
064d808d 751Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
6d9e79cc 752 //
753 // Return the position of the corruption
754 //
755 return encoding & 7;
756}
757
3fe61b77 758//_____________________________________________________________________________
759Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
760{
761 //
762 // Generates the cluster.
763 //
764
765 // Get the digits
615f0826 766 fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
767 fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline();
f5375dcb 768
3fe61b77 769 // This is to take care of switched off super modules
b72f4eaf 770 if (!fDigits->HasData()) return kFALSE;
3fe61b77 771
064d808d 772 fIndexes = fDigitsManager->GetIndexes(det);
b72f4eaf 773 if (fIndexes->IsAllocated() == kFALSE) {
774 AliError("Indexes do not exist!");
775 return kFALSE;
776 }
064d808d 777
3fe61b77 778 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
b72f4eaf 779 if (!calibration) {
780 AliFatal("No AliTRDcalibDB instance available\n");
781 return kFALSE;
782 }
3fe61b77 783
3a039a31 784 if (!fReconstructor){
785 AliError("Reconstructor not set\n");
786 return kFALSE;
787 }
c5f589b9 788
064d808d 789 fMaxThresh = fReconstructor->GetRecoParam()->GetClusMaxThresh();
790 fSigThresh = fReconstructor->GetRecoParam()->GetClusSigThresh();
791 fMinMaxCutSigma = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
792 fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
3fe61b77 793
064d808d 794 Int_t istack = fIndexes->GetStack();
795 fLayer = fIndexes->GetLayer();
796 Int_t isector = fIndexes->GetSM();
3fe61b77 797
798 // Start clustering in the chamber
799
064d808d 800 fDet = AliTRDgeometry::GetDetector(fLayer,istack,isector);
801 if (fDet != det) {
828c6f80 802 AliError("Strange Detector number Missmatch!");
b65e5048 803 return kFALSE;
804 }
3fe61b77 805
839f6b9b 806 AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
807
3fe61b77 808 // TRD space point transformation
809 fTransform->SetDetector(det);
810
064d808d 811 Int_t iGeoLayer = AliGeomManager::kTRD1 + fLayer;
053767a4 812 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
064d808d 813 fVolid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
3fe61b77 814
a5b99acd 815 if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
816 AddTrackletsToArray();
817
1b3a719f 818 fColMax = fDigits->GetNcol();
819 //Int_t nRowMax = fDigits->GetNrow();
820 fTimeTotal = fDigits->GetNtime();
3fe61b77 821
69a850a0 822 // Check consistency between OCDB and raw data
4e6823c2 823 Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
824 if ((nTimeOCDB > -1) &&
825 (fTimeTotal != nTimeOCDB)) {
828c6f80 826 AliError(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d])"
827 ,fTimeTotal,calibration->GetNumberOfTimeBinsDCS()));
69a850a0 828 }
829
3fe61b77 830 // Detector wise calibration object for the gain factors
064d808d 831 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
3fe61b77 832 // Calibration object with pad wise values for the gain factors
064d808d 833 fCalGainFactorROC = calibration->GetGainFactorROC(fDet);
3fe61b77 834 // Calibration value for chamber wise gain factor
064d808d 835 fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
0e09df31 836
df83a620 837 // Detector wise calibration object for the noise
064d808d 838 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
df83a620 839 // Calibration object with pad wise values for the noise
064d808d 840 fCalNoiseROC = calibration->GetNoiseROC(fDet);
df83a620 841 // Calibration value for chamber wise noise
064d808d 842 fCalNoiseDetValue = calNoiseDet->GetValue(fDet);
1b3a719f 843
844 // Calibration object with the pad status
845 fCalPadStatusROC = calibration->GetPadStatusROC(fDet);
6b0d4ad5 846
a2fbb6ec 847 SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
848 SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
b72f4eaf 849 SetBit(kHLT, fReconstructor->IsHLT());
6b0d4ad5 850
064d808d 851 firstClusterROC = -1;
852 fClusterROC = 0;
3fe61b77 853
b72f4eaf 854 // Apply the gain and the tail cancelation via digital filter
a2fbb6ec 855 if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
023b669c 856
486c2339 857 MaxStruct curr, last;
064d808d 858 Int_t nMaximas = 0, nCorrupted = 0;
f5375dcb 859
064d808d 860 // Here the clusterfining is happening
861
615f0826 862 for(curr.time = 0; curr.time < fTimeTotal; curr.time++){
863 while(fIndexes->NextRCIndex(curr.row, curr.col)){
864 if(IsMaximum(curr, curr.padStatus, &curr.signals[0])){
865 if(last.row>-1){
866 if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
b72f4eaf 867 CreateCluster(last);
868 }
615f0826 869 last=curr; curr.fivePad=kFALSE;
b72f4eaf 870 }
b72f4eaf 871 }
872 }
615f0826 873 if(last.row>-1) CreateCluster(last);
df83a620 874
a2fbb6ec 875 if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
876 TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
064d808d 877 (*fDebugStream) << "MakeClusters"
b72f4eaf 878 << "Detector=" << det
879 << "NMaxima=" << nMaximas
880 << "NClusters=" << fClusterROC
881 << "NCorrupted=" << nCorrupted
882 << "\n";
df83a620 883 }
b72f4eaf 884 if (TestBit(kLabels)) AddLabels();
f5375dcb 885
064d808d 886 return kTRUE;
f5375dcb 887
064d808d 888}
3fe61b77 889
064d808d 890//_____________________________________________________________________________
1b3a719f 891Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
064d808d 892{
893 //
894 // Returns true if this row,col,time combination is a maximum.
895 // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
896 //
c96d03ba 897
615f0826 898 Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
899 Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
c96d03ba 900 if(Signals[1] < fMaxThresh) return kFALSE;
3fe61b77 901
615f0826 902 Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
c96d03ba 903 if (Signals[1] < noiseMiddleThresh) return kFALSE;
f5375dcb 904
615f0826 905 if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
3fe61b77 906
b72f4eaf 907 UChar_t status[3]={
615f0826 908 fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
909 ,fCalPadStatusROC->GetStatus(Max.col, Max.row)
910 ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
b72f4eaf 911 };
486c2339 912
615f0826 913 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
19dd8eb9 914 Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
915 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
916 Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
c96d03ba 917
486c2339 918 if(!(status[0] | status[1] | status[2])) {//all pads are good
c96d03ba 919 if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) {
920 if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
5a8db426 921 if(Signals[0]<0)Signals[0]=0;
922 if(Signals[2]<0)Signals[2]=0;
b72f4eaf 923 Float_t noiseSumThresh = fMinLeftRightCutSigma
924 * fCalNoiseDetValue
615f0826 925 * fCalNoiseROC->GetValue(Max.col, Max.row);
c96d03ba 926 if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
b72f4eaf 927 padStatus = 0;
928 return kTRUE;
eb52b657 929 }
064d808d 930 }
b72f4eaf 931 } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
5a8db426 932 if(Signals[0]<0)Signals[0]=0;
933 if(Signals[2]<0)Signals[2]=0;
c96d03ba 934 if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
935 Signals[2]=0;
486c2339 936 SetPadStatus(status[2], padStatus);
937 return kTRUE;
938 }
c96d03ba 939 else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
940 Signals[0]=0;
486c2339 941 SetPadStatus(status[0], padStatus);
942 return kTRUE;
943 }
c96d03ba 944 else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
615f0826 945 Signals[1] = (Short_t)(fMaxThresh + 0.5f);
486c2339 946 SetPadStatus(status[1], padStatus);
064d808d 947 return kTRUE;
948 }
949 }
950 return kFALSE;
951}
3fe61b77 952
064d808d 953//_____________________________________________________________________________
1b3a719f 954Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
064d808d 955{
956 //
957 // Look for 5 pad cluster with minimum in the middle
958 // Gives back the ratio
959 //
615f0826 960
961 if (ThisMax.col >= fColMax - 3) return kFALSE;
962 Float_t gain;
963 if (ThisMax.col < fColMax - 5){
964 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
965 if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 966 return kFALSE;
064d808d 967 }
615f0826 968 if (ThisMax.col > 1) {
969 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
970 if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 971 return kFALSE;
972 }
973
486c2339 974 const Float_t kEpsilon = 0.01;
615f0826 975 Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
976 NeighbourMax.signals[1], NeighbourMax.signals[2]};
486c2339 977
978 // Unfold the two maxima and set the signal on
979 // the overlapping pad to the ratio
1b3a719f 980 Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
615f0826 981 ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
982 NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
983 ThisMax.fivePad=kTRUE;
984 NeighbourMax.fivePad=kTRUE;
486c2339 985 return kTRUE;
c96d03ba 986
064d808d 987}
eb52b657 988
064d808d 989//_____________________________________________________________________________
486c2339 990void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
064d808d 991{
992 //
3d0c7d6d 993 // Creates a cluster at the given position and saves it in fRecPoints
064d808d 994 //
eb52b657 995
c96d03ba 996 Int_t nPadCount = 1;
615f0826 997 Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
c96d03ba 998 if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
999
615f0826 1000 AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
c96d03ba 1001 cluster.SetNPads(nPadCount);
b72f4eaf 1002 if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1003 else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1004 else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
3fe61b77 1005
615f0826 1006 cluster.SetFivePad(Max.fivePad);
b72f4eaf 1007 // set pads status for the cluster
1008 UChar_t maskPosition = GetCorruption(Max.padStatus);
1009 if (maskPosition) {
1010 cluster.SetPadMaskedPosition(maskPosition);
1011 cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1012 }
064d808d 1013
1014 // Transform the local cluster coordinates into calibrated
1015 // space point positions defined in the local tracking system.
1016 // Here the calibration for T0, Vdrift and ExB is applied as well.
b72f4eaf 1017 if(!fTransform->Transform(&cluster)) return;
1018 // Temporarily store the Max.Row, column and time bin of the center pad
1019 // Used to later on assign the track indices
615f0826 1020 cluster.SetLabel(Max.row, 0);
1021 cluster.SetLabel(Max.col, 1);
1022 cluster.SetLabel(Max.time,2);
c96d03ba 1023
b72f4eaf 1024 //needed for HLT reconstruction
1025 AddClusterToArray(&cluster);
1026
1027 // Store the index of the first cluster in the current ROC
1028 if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1029
1030 fNoOfClusters++;
1031 fClusterROC++;
3fe61b77 1032}
1033
3d0c7d6d 1034//_____________________________________________________________________________
1035void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1036{
1037 // Look to the right
1038 Int_t ii = 1;
615f0826 1039 while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
3d0c7d6d 1040 nPadCount++;
1041 ii++;
615f0826 1042 if (Max.col < ii) break;
3d0c7d6d 1043 }
1044 // Look to the left
1045 ii = 1;
615f0826 1046 while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
3d0c7d6d 1047 nPadCount++;
1048 ii++;
615f0826 1049 if (Max.col+ii >= fColMax) break;
3d0c7d6d 1050 }
1051
1052 // Store the amplitudes of the pads in the cluster for later analysis
1053 // and check whether one of these pads is masked in the database
615f0826 1054 signals[2]=Max.signals[0];
1055 signals[3]=Max.signals[1];
1056 signals[4]=Max.signals[2];
1057 Float_t gain;
3d0c7d6d 1058 for(Int_t i = 0; i<2; i++)
1059 {
615f0826 1060 if(Max.col+i >= 3){
1061 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1062 signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1063 }
1064 if(Max.col+3-i < fColMax){
1065 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1066 signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1067 }
3d0c7d6d 1068 }
1069 /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1070 if ((jPad >= 0) && (jPad < fColMax))
1071 signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1072 }*/
1073}
1074
1075//_____________________________________________________________________________
e7295a3a 1076void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
3d0c7d6d 1077{
1078 //
1079 // Add a cluster to the array
1080 //
1081
1082 Int_t n = RecPoints()->GetEntriesFast();
1083 if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1084 new((*RecPoints())[n]) AliTRDcluster(*cluster);
1085}
1086
a5b99acd 1087//_____________________________________________________________________________
1088void AliTRDclusterizer::AddTrackletsToArray()
1089{
1090 //
1091 // Add the online tracklets of this chamber to the array
1092 //
1093
1094 UInt_t* trackletword;
1095 for(Int_t side=0; side<2; side++)
1096 {
1097 Int_t trkl=0;
1098 trackletword=fTrackletContainer[side];
97be9934 1099 while(trackletword[trkl]>0){
a5b99acd 1100 Int_t n = TrackletsArray()->GetEntriesFast();
35943b1e 1101 AliTRDtrackletWord tmp(trackletword[trkl]);
1102 new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
a5b99acd 1103 trkl++;
97be9934 1104 }
a5b99acd 1105 }
1106}
1107
3fe61b77 1108//_____________________________________________________________________________
6cd9a21f 1109Bool_t AliTRDclusterizer::AddLabels()
3551db50 1110{
1111 //
3fe61b77 1112 // Add the track indices to the found clusters
3551db50 1113 //
1114
3fe61b77 1115 const Int_t kNclus = 3;
1116 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1117 const Int_t kNtrack = kNdict * kNclus;
1118
1119 Int_t iClusterROC = 0;
1120
1121 Int_t row = 0;
1122 Int_t col = 0;
1123 Int_t time = 0;
1124 Int_t iPad = 0;
1125
1126 // Temporary array to collect the track indices
6cd9a21f 1127 Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
3fe61b77 1128
1129 // Loop through the dictionary arrays one-by-one
1130 // to keep memory consumption low
b65e5048 1131 AliTRDarrayDictionary *tracksIn = 0; //mod
3fe61b77 1132 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1133
1134 // tracksIn should be expanded beforehand!
6cd9a21f 1135 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
3fe61b77 1136
1137 // Loop though the clusters found in this ROC
6cd9a21f 1138 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1139
3fe61b77 1140 AliTRDcluster *cluster = (AliTRDcluster *)
b65e5048 1141 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
3fe61b77 1142 row = cluster->GetLabel(0);
1143 col = cluster->GetLabel(1);
1144 time = cluster->GetLabel(2);
1145
1146 for (iPad = 0; iPad < kNclus; iPad++) {
b65e5048 1147 Int_t iPadCol = col - 1 + iPad;
1148 Int_t index = tracksIn->GetData(row,iPadCol,time); //Modification of -1 in Track
1149 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
3fe61b77 1150 }
1151
1152 }
1153
1154 }
1155
1156 // Copy the track indices into the cluster
1157 // Loop though the clusters found in this ROC
6cd9a21f 1158 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1159
3fe61b77 1160 AliTRDcluster *cluster = (AliTRDcluster *)
1161 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1162 cluster->SetLabel(-9999,0);
1163 cluster->SetLabel(-9999,1);
1164 cluster->SetLabel(-9999,2);
1165
1166 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1167
1168 }
1169
1170 delete [] idxTracks;
1171
1172 return kTRUE;
1173
1174}
1175
3fe61b77 1176//_____________________________________________________________________________
e7295a3a 1177Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
3fe61b77 1178{
1179 //
1180 // Method to unfold neighbouring maxima.
1181 // The charge ratio on the overlapping pad is calculated
1182 // until there is no more change within the range given by eps.
1183 // The resulting ratio is then returned to the calling method.
1184 //
1185
1186 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1187 if (!calibration) {
3fe61b77 1188 AliError("No AliTRDcalibDB instance available\n");
1189 return kFALSE;
6d50f529 1190 }
3fe61b77 1191
1192 Int_t irc = 0;
1193 Int_t itStep = 0; // Count iteration steps
1194
1195 Double_t ratio = 0.5; // Start value for ratio
1196 Double_t prevRatio = 0.0; // Store previous ratio
1197
1198 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1199 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1200 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1201
1202 // Start the iteration
1203 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1204
1205 itStep++;
1206 prevRatio = ratio;
1207
1208 // Cluster position according to charge ratio
1209 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
064d808d 1210 / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
3fe61b77 1211 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1212 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1213
1214 // Set cluster charge ratio
064d808d 1215 irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
3fe61b77 1216 Double_t ampLeft = padSignal[1] / newSignal[1];
064d808d 1217 irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
3fe61b77 1218 Double_t ampRight = padSignal[3] / newSignal[1];
1219
1220 // Apply pad response to parameters
053767a4 1221 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1222 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
3fe61b77 1223
1224 // Calculate new overlapping ratio
1225 ratio = TMath::Min((Double_t) 1.0
1226 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1227
b43a3e17 1228 }
88719a08 1229
3fe61b77 1230 return ratio;
1231
1232}
88719a08 1233
3fe61b77 1234//_____________________________________________________________________________
064d808d 1235void AliTRDclusterizer::TailCancelation()
3fe61b77 1236{
1237 //
1238 // Applies the tail cancelation and gain factors:
1b3a719f 1239 // Transform fDigits to fDigits
3fe61b77 1240 //
1241
1242 Int_t iRow = 0;
1243 Int_t iCol = 0;
1244 Int_t iTime = 0;
1245
1b3a719f 1246 Double_t *inADC = new Double_t[fTimeTotal]; // ADC data before tail cancellation
064d808d 1247 Double_t *outADC = new Double_t[fTimeTotal]; // ADC data after tail cancellation
1b3a719f 1248
a2fbb6ec 1249 TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
cb5522da 1250 Bool_t debugStreaming = fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
064d808d 1251 while(fIndexes->NextRCIndex(iRow, iCol))
3fe61b77 1252 {
f5375dcb 1253 Bool_t corrupted = kFALSE;
5a8db426 1254 if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1255
cb5522da 1256 // Save data into the temporary processing array and substract the baseline,
1257 // since DeConvExp does not expect a baseline
064d808d 1258 for (iTime = 0; iTime < fTimeTotal; iTime++)
cb5522da 1259 inADC[iTime] = fDigits->GetData(iRow,iCol,iTime)-fBaseline;
5a8db426 1260
cb5522da 1261 if(debugStreaming){
1262 for (iTime = 0; iTime < fTimeTotal; iTime++)
5a8db426 1263 (*fDebugStream) << "TailCancellation"
1264 << "col=" << iCol
1265 << "row=" << iRow
1266 << "time=" << iTime
1267 << "inADC=" << inADC[iTime]
1268 << "outADC=" << outADC[iTime]
1269 << "corrupted=" << corrupted
1270 << "\n";
cb5522da 1271 }
5a8db426 1272
b65e5048 1273 if (!corrupted)
1274 {
1275 // Apply the tail cancelation via the digital filter
1276 // (only for non-coorupted pads)
1b3a719f 1277 DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
b65e5048 1278 }
5a8db426 1279 else memcpy(&outADC[0],&inADC[0],fTimeTotal*sizeof(inADC[0]));
3fe61b77 1280
cb5522da 1281 // Save tailcancalled data and add the baseline
5a8db426 1282 for(iTime = 0; iTime < fTimeTotal; iTime++)
cb5522da 1283 fDigits->SetData(iRow,iCol,iTime,(Short_t)(outADC[iTime] + fBaseline + 0.5));
5a8db426 1284
3fe61b77 1285 } // while irow icol
1b3a719f 1286
3fe61b77 1287 delete [] inADC;
1288 delete [] outADC;
1289
1290 return;
1291
1292}
1293
1294//_____________________________________________________________________________
064d808d 1295void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1296 ,const Int_t n, const Int_t nexp)
3fe61b77 1297{
1298 //
1299 // Tail cancellation by deconvolution for PASA v4 TRF
1300 //
1301
1302 Double_t rates[2];
1303 Double_t coefficients[2];
1304
1305 // Initialization (coefficient = alpha, rates = lambda)
1306 Double_t r1 = 1.0;
1307 Double_t r2 = 1.0;
1308 Double_t c1 = 0.5;
1309 Double_t c2 = 0.5;
1310
1311 if (nexp == 1) { // 1 Exponentials
1312 r1 = 1.156;
1313 r2 = 0.130;
1314 c1 = 0.066;
1315 c2 = 0.000;
1316 }
1317 if (nexp == 2) { // 2 Exponentials
181c7f7e 1318 Double_t par[4];
a2fbb6ec 1319 fReconstructor->GetRecoParam()->GetTCParams(par);
181c7f7e 1320 r1 = par[0];//1.156;
1321 r2 = par[1];//0.130;
1322 c1 = par[2];//0.114;
1323 c2 = par[3];//0.624;
3fe61b77 1324 }
1325
1326 coefficients[0] = c1;
1327 coefficients[1] = c2;
1328
1329 Double_t dt = 0.1;
1330
1331 rates[0] = TMath::Exp(-dt/(r1));
1332 rates[1] = TMath::Exp(-dt/(r2));
1333
1334 Int_t i = 0;
1335 Int_t k = 0;
1336
1337 Double_t reminder[2];
1338 Double_t correction = 0.0;
1339 Double_t result = 0.0;
1340
1341 // Attention: computation order is important
1342 for (k = 0; k < nexp; k++) {
1343 reminder[k] = 0.0;
1344 }
1345
1346 for (i = 0; i < n; i++) {
1347
1348 result = (source[i] - correction); // No rescaling
1349 target[i] = result;
1350
1351 for (k = 0; k < nexp; k++) {
1352 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1353 }
1354
1355 correction = 0.0;
1356 for (k = 0; k < nexp; k++) {
1357 correction += reminder[k];
1358 }
1359
1360 }
6d50f529 1361
1362}
1363
1364//_____________________________________________________________________________
1365void AliTRDclusterizer::ResetRecPoints()
1366{
1367 //
1368 // Resets the list of rec points
1369 //
1370
1371 if (fRecPoints) {
1372 fRecPoints->Delete();
48f8adf3 1373 delete fRecPoints;
6d50f529 1374 }
6d50f529 1375}
1376
1377//_____________________________________________________________________________
66f6bfd9 1378TClonesArray *AliTRDclusterizer::RecPoints()
6d50f529 1379{
1380 //
1381 // Returns the list of rec points
1382 //
1383
1384 if (!fRecPoints) {
48f8adf3 1385 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
8ae98148 1386 // determine number of clusters which has to be allocated
1387 Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
8ae98148 1388
1389 fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
48f8adf3 1390 }
1391 //SetClustersOwner(kTRUE);
1392 AliTRDReconstructor::SetClusters(0x0);
6d50f529 1393 }
6d50f529 1394 return fRecPoints;
eb52b657 1395
3551db50 1396}
fc546d21 1397
a5b99acd 1398//_____________________________________________________________________________
1399TClonesArray *AliTRDclusterizer::TrackletsArray()
1400{
1401 //
1402 // Returns the list of rec points
1403 //
1404
1405 if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
94ac01ef 1406 fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
a5b99acd 1407 //SetClustersOwner(kTRUE);
1408 //AliTRDReconstructor::SetTracklets(0x0);
1409 }
1410 return fTracklets;
1411
1412}
1413