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