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