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