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