]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterizer.cxx
Use configuration from the DCS to switch on/off offline tail cancellation dependent...
[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
5f5f179e 858 // Use the configuration from the DCS to find out whether online
859 // tail cancellation was applied
860 if(!calibration->HasOnlineTailCancellation()) TailCancelation(recoParam);
023b669c 861
486c2339 862 MaxStruct curr, last;
064d808d 863 Int_t nMaximas = 0, nCorrupted = 0;
f5375dcb 864
064d808d 865 // Here the clusterfining is happening
866
8c499dbf 867 for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
615f0826 868 while(fIndexes->NextRCIndex(curr.row, curr.col)){
1ab4da8c 869 if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
615f0826 870 if(last.row>-1){
634a1625 871 if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
b72f4eaf 872 CreateCluster(last);
873 }
615f0826 874 last=curr; curr.fivePad=kFALSE;
b72f4eaf 875 }
b72f4eaf 876 }
877 }
615f0826 878 if(last.row>-1) CreateCluster(last);
df83a620 879
fc0882f3 880 if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
a2fbb6ec 881 TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
064d808d 882 (*fDebugStream) << "MakeClusters"
b72f4eaf 883 << "Detector=" << det
884 << "NMaxima=" << nMaximas
885 << "NClusters=" << fClusterROC
886 << "NCorrupted=" << nCorrupted
887 << "\n";
df83a620 888 }
b72f4eaf 889 if (TestBit(kLabels)) AddLabels();
f5375dcb 890
064d808d 891 return kTRUE;
f5375dcb 892
064d808d 893}
3fe61b77 894
064d808d 895//_____________________________________________________________________________
1b3a719f 896Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
064d808d 897{
898 //
899 // Returns true if this row,col,time combination is a maximum.
900 // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
901 //
c96d03ba 902
615f0826 903 Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
904 Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
1ab4da8c 905 if(Signals[1] <= fMaxThresh) return kFALSE;
3fe61b77 906
634a1625 907 if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
f5375dcb 908
634a1625 909 Short_t noiseMiddleThresh = (Short_t)(fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row));
910 if (Signals[1] <= noiseMiddleThresh) return kFALSE;
3fe61b77 911
b72f4eaf 912 UChar_t status[3]={
615f0826 913 fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
914 ,fCalPadStatusROC->GetStatus(Max.col, Max.row)
915 ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
b72f4eaf 916 };
486c2339 917
c66dc420 918 Short_t signal(0);
919 if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
920 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
921 Signals[0] = (Short_t)((signal - fBaseline) / gain + 0.5f);
922 } else Signals[0] = 0;
923 if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
924 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
925 Signals[2] = (Short_t)((signal - fBaseline) / gain + 0.5f);
926 } else Signals[2] = 0;
c96d03ba 927
486c2339 928 if(!(status[0] | status[1] | status[2])) {//all pads are good
c96d03ba 929 if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) {
1ab4da8c 930 if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
5a8db426 931 if(Signals[0]<0)Signals[0]=0;
932 if(Signals[2]<0)Signals[2]=0;
1ab4da8c 933 Short_t noiseSumThresh = (Short_t)(fMinLeftRightCutSigma * fCalNoiseDetValue
934 * fCalNoiseROC->GetValue(Max.col, Max.row));
935 if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
b72f4eaf 936 padStatus = 0;
937 return kTRUE;
eb52b657 938 }
064d808d 939 }
b72f4eaf 940 } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
5a8db426 941 if(Signals[0]<0)Signals[0]=0;
942 if(Signals[2]<0)Signals[2]=0;
1ab4da8c 943 if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
c96d03ba 944 Signals[2]=0;
486c2339 945 SetPadStatus(status[2], padStatus);
946 return kTRUE;
947 }
1ab4da8c 948 else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
c96d03ba 949 Signals[0]=0;
486c2339 950 SetPadStatus(status[0], padStatus);
951 return kTRUE;
952 }
1ab4da8c 953 else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
954 Signals[1] = fMaxThresh;
486c2339 955 SetPadStatus(status[1], padStatus);
064d808d 956 return kTRUE;
957 }
958 }
959 return kFALSE;
960}
3fe61b77 961
064d808d 962//_____________________________________________________________________________
1b3a719f 963Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
064d808d 964{
965 //
966 // Look for 5 pad cluster with minimum in the middle
967 // Gives back the ratio
968 //
615f0826 969
970 if (ThisMax.col >= fColMax - 3) return kFALSE;
971 Float_t gain;
972 if (ThisMax.col < fColMax - 5){
973 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
974 if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 975 return kFALSE;
064d808d 976 }
615f0826 977 if (ThisMax.col > 1) {
978 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
979 if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 980 return kFALSE;
981 }
982
486c2339 983 const Float_t kEpsilon = 0.01;
615f0826 984 Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
985 NeighbourMax.signals[1], NeighbourMax.signals[2]};
486c2339 986
987 // Unfold the two maxima and set the signal on
988 // the overlapping pad to the ratio
1b3a719f 989 Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
615f0826 990 ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
991 NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
992 ThisMax.fivePad=kTRUE;
993 NeighbourMax.fivePad=kTRUE;
486c2339 994 return kTRUE;
c96d03ba 995
064d808d 996}
eb52b657 997
064d808d 998//_____________________________________________________________________________
486c2339 999void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
064d808d 1000{
1001 //
3d0c7d6d 1002 // Creates a cluster at the given position and saves it in fRecPoints
064d808d 1003 //
eb52b657 1004
c96d03ba 1005 Int_t nPadCount = 1;
615f0826 1006 Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
d03bfe2b 1007 if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
c96d03ba 1008
615f0826 1009 AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
c96d03ba 1010 cluster.SetNPads(nPadCount);
b72f4eaf 1011 if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1012 else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1013 else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
3fe61b77 1014
615f0826 1015 cluster.SetFivePad(Max.fivePad);
b72f4eaf 1016 // set pads status for the cluster
1017 UChar_t maskPosition = GetCorruption(Max.padStatus);
1018 if (maskPosition) {
1019 cluster.SetPadMaskedPosition(maskPosition);
1020 cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1021 }
8fc736d7 1022 cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
064d808d 1023
1024 // Transform the local cluster coordinates into calibrated
1025 // space point positions defined in the local tracking system.
1026 // Here the calibration for T0, Vdrift and ExB is applied as well.
d03bfe2b 1027 if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
82ddb093 1028
b72f4eaf 1029 // Temporarily store the Max.Row, column and time bin of the center pad
1030 // Used to later on assign the track indices
615f0826 1031 cluster.SetLabel(Max.row, 0);
1032 cluster.SetLabel(Max.col, 1);
1033 cluster.SetLabel(Max.time,2);
c96d03ba 1034
b72f4eaf 1035 //needed for HLT reconstruction
8c499dbf 1036 AddClusterToArray(&cluster);
b72f4eaf 1037
1038 // Store the index of the first cluster in the current ROC
1039 if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1040
1041 fNoOfClusters++;
1042 fClusterROC++;
3fe61b77 1043}
1044
3d0c7d6d 1045//_____________________________________________________________________________
1046void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1047{
c66dc420 1048// Calculate number of pads/cluster and
1049// ADC signals at position 0, 1, 5 and 6
1050
1051 Float_t gain(1.); Short_t signal(0.);
1052 // Store the amplitudes of the pads in the cluster for later analysis
1053 // and check whether one of these pads is masked in the database
1054 signals[3]=Max.signals[1];
1055 Int_t ipad(1), jpad(0);
3d0c7d6d 1056 // Look to the right
c66dc420 1057 while((jpad = Max.col-ipad)){
1058 if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1059 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1060 signal = (Short_t)((signal - fBaseline) / gain + 0.5f);
1061 if(signal<fSigThresh) break; // signal under threshold
3d0c7d6d 1062 nPadCount++;
c66dc420 1063 if(ipad<=3) signals[3 - ipad] = signal;
1064 ipad++;
3d0c7d6d 1065 }
c66dc420 1066 ipad=1;
3d0c7d6d 1067 // Look to the left
c66dc420 1068 while((jpad = Max.col+ipad)<fColMax){
1069 if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1070 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1071 signal = (Short_t)((signal - fBaseline) / gain + 0.5f);
1072 if(signal<fSigThresh) break; // signal under threshold
3d0c7d6d 1073 nPadCount++;
c66dc420 1074 if(ipad<=3) signals[3 + ipad] = signal;
1075 ipad++;
3d0c7d6d 1076 }
1077
c66dc420 1078 AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1079 , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
3d0c7d6d 1080}
1081
1082//_____________________________________________________________________________
e7295a3a 1083void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
3d0c7d6d 1084{
1085 //
1086 // Add a cluster to the array
1087 //
1088
1089 Int_t n = RecPoints()->GetEntriesFast();
1090 if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1091 new((*RecPoints())[n]) AliTRDcluster(*cluster);
1092}
1093
a5b99acd 1094//_____________________________________________________________________________
1095void AliTRDclusterizer::AddTrackletsToArray()
1096{
1097 //
1098 // Add the online tracklets of this chamber to the array
1099 //
1100
1101 UInt_t* trackletword;
1102 for(Int_t side=0; side<2; side++)
1103 {
1104 Int_t trkl=0;
1105 trackletword=fTrackletContainer[side];
97be9934 1106 while(trackletword[trkl]>0){
a5b99acd 1107 Int_t n = TrackletsArray()->GetEntriesFast();
35943b1e 1108 AliTRDtrackletWord tmp(trackletword[trkl]);
1109 new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
a5b99acd 1110 trkl++;
97be9934 1111 }
a5b99acd 1112 }
1113}
1114
3fe61b77 1115//_____________________________________________________________________________
6cd9a21f 1116Bool_t AliTRDclusterizer::AddLabels()
3551db50 1117{
1118 //
3fe61b77 1119 // Add the track indices to the found clusters
3551db50 1120 //
1121
3fe61b77 1122 const Int_t kNclus = 3;
1123 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1124 const Int_t kNtrack = kNdict * kNclus;
1125
1126 Int_t iClusterROC = 0;
1127
1128 Int_t row = 0;
1129 Int_t col = 0;
1130 Int_t time = 0;
1131 Int_t iPad = 0;
1132
1133 // Temporary array to collect the track indices
6cd9a21f 1134 Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
3fe61b77 1135
1136 // Loop through the dictionary arrays one-by-one
1137 // to keep memory consumption low
b65e5048 1138 AliTRDarrayDictionary *tracksIn = 0; //mod
3fe61b77 1139 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1140
1141 // tracksIn should be expanded beforehand!
6cd9a21f 1142 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
3fe61b77 1143
1144 // Loop though the clusters found in this ROC
6cd9a21f 1145 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1146
3fe61b77 1147 AliTRDcluster *cluster = (AliTRDcluster *)
b65e5048 1148 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
3fe61b77 1149 row = cluster->GetLabel(0);
1150 col = cluster->GetLabel(1);
1151 time = cluster->GetLabel(2);
1152
1153 for (iPad = 0; iPad < kNclus; iPad++) {
b65e5048 1154 Int_t iPadCol = col - 1 + iPad;
1155 Int_t index = tracksIn->GetData(row,iPadCol,time); //Modification of -1 in Track
1156 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
3fe61b77 1157 }
1158
1159 }
1160
1161 }
1162
1163 // Copy the track indices into the cluster
1164 // Loop though the clusters found in this ROC
6cd9a21f 1165 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1166
3fe61b77 1167 AliTRDcluster *cluster = (AliTRDcluster *)
1168 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1169 cluster->SetLabel(-9999,0);
1170 cluster->SetLabel(-9999,1);
1171 cluster->SetLabel(-9999,2);
1172
1173 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1174
1175 }
1176
1177 delete [] idxTracks;
1178
1179 return kTRUE;
1180
1181}
1182
3fe61b77 1183//_____________________________________________________________________________
e7295a3a 1184Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
3fe61b77 1185{
1186 //
1187 // Method to unfold neighbouring maxima.
1188 // The charge ratio on the overlapping pad is calculated
1189 // until there is no more change within the range given by eps.
1190 // The resulting ratio is then returned to the calling method.
1191 //
1192
1193 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1194 if (!calibration) {
3fe61b77 1195 AliError("No AliTRDcalibDB instance available\n");
1196 return kFALSE;
6d50f529 1197 }
3fe61b77 1198
1199 Int_t irc = 0;
1200 Int_t itStep = 0; // Count iteration steps
1201
1202 Double_t ratio = 0.5; // Start value for ratio
1203 Double_t prevRatio = 0.0; // Store previous ratio
1204
1205 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1206 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1207 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1208
1209 // Start the iteration
1210 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1211
1212 itStep++;
1213 prevRatio = ratio;
1214
1215 // Cluster position according to charge ratio
1216 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
064d808d 1217 / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
3fe61b77 1218 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1219 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1220
1221 // Set cluster charge ratio
064d808d 1222 irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
3fe61b77 1223 Double_t ampLeft = padSignal[1] / newSignal[1];
064d808d 1224 irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
3fe61b77 1225 Double_t ampRight = padSignal[3] / newSignal[1];
1226
1227 // Apply pad response to parameters
053767a4 1228 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1229 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
3fe61b77 1230
1231 // Calculate new overlapping ratio
1232 ratio = TMath::Min((Double_t) 1.0
1233 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1234
b43a3e17 1235 }
88719a08 1236
3fe61b77 1237 return ratio;
1238
1239}
88719a08 1240
3fe61b77 1241//_____________________________________________________________________________
fc0882f3 1242void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
3fe61b77 1243{
1244 //
fc0882f3 1245 // Applies the tail cancelation
3fe61b77 1246 //
1247
1248 Int_t iRow = 0;
1249 Int_t iCol = 0;
1250 Int_t iTime = 0;
1251
a2fbb6ec 1252 TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
fc0882f3 1253 Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1254 Int_t nexp = recoParam->GetTCnexp();
064d808d 1255 while(fIndexes->NextRCIndex(iRow, iCol))
3fe61b77 1256 {
664030bc 1257 // if corrupted then don't make the tail cancallation
1258 if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
5a8db426 1259
cb5522da 1260 if(debugStreaming){
1ab4da8c 1261 for (iTime = 0; iTime < fTimeTotal; iTime++)
1262 (*fDebugStream) << "TailCancellation"
1263 << "col=" << iCol
1264 << "row=" << iRow
1265 << "\n";
cb5522da 1266 }
5a8db426 1267
664030bc 1268 // Apply the tail cancelation via the digital filter
1ab4da8c 1269 DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
5a8db426 1270
3fe61b77 1271 } // while irow icol
1b3a719f 1272
3fe61b77 1273 return;
1274
1275}
1276
1277//_____________________________________________________________________________
1ab4da8c 1278void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
3fe61b77 1279{
1280 //
1281 // Tail cancellation by deconvolution for PASA v4 TRF
1282 //
1283
664030bc 1284 Float_t rates[2];
1285 Float_t coefficients[2];
3fe61b77 1286
1287 // Initialization (coefficient = alpha, rates = lambda)
664030bc 1288 Float_t r1 = 1.0;
1289 Float_t r2 = 1.0;
1290 Float_t c1 = 0.5;
1291 Float_t c2 = 0.5;
3fe61b77 1292
1293 if (nexp == 1) { // 1 Exponentials
1294 r1 = 1.156;
1295 r2 = 0.130;
1296 c1 = 0.066;
1297 c2 = 0.000;
1298 }
1299 if (nexp == 2) { // 2 Exponentials
181c7f7e 1300 Double_t par[4];
a2fbb6ec 1301 fReconstructor->GetRecoParam()->GetTCParams(par);
181c7f7e 1302 r1 = par[0];//1.156;
1303 r2 = par[1];//0.130;
1304 c1 = par[2];//0.114;
1305 c2 = par[3];//0.624;
3fe61b77 1306 }
1307
1308 coefficients[0] = c1;
1309 coefficients[1] = c2;
1310
fc0882f3 1311 Double_t dt = 0.1;
3fe61b77 1312
1313 rates[0] = TMath::Exp(-dt/(r1));
634a1625 1314 rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
3fe61b77 1315
634a1625 1316 Float_t reminder[2] = { .0, .0 };
664030bc 1317 Float_t correction = 0.0;
1318 Float_t result = 0.0;
3fe61b77 1319
634a1625 1320 for (int i = 0; i < nTime; i++) {
3fe61b77 1321
1ab4da8c 1322 result = arr[i] - correction - fBaseline; // No rescaling
1323 arr[i] = (Short_t)(result + fBaseline + 0.5f);
3fe61b77 1324
3fe61b77 1325 correction = 0.0;
634a1625 1326 for (int k = 0; k < 2; k++) {
1327 correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
3fe61b77 1328 }
1329
1330 }
6d50f529 1331
1332}
1333
1334//_____________________________________________________________________________
1335void AliTRDclusterizer::ResetRecPoints()
1336{
1337 //
1338 // Resets the list of rec points
1339 //
1340
1341 if (fRecPoints) {
da0a4b87 1342 fRecPoints->Clear();
1343 fNoOfClusters = 0;
1344 // delete fRecPoints;
6d50f529 1345 }
6d50f529 1346}
1347
1348//_____________________________________________________________________________
66f6bfd9 1349TClonesArray *AliTRDclusterizer::RecPoints()
6d50f529 1350{
1351 //
1352 // Returns the list of rec points
1353 //
1354
1355 if (!fRecPoints) {
48f8adf3 1356 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
8ae98148 1357 // determine number of clusters which has to be allocated
1358 Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
8ae98148 1359
1360 fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
48f8adf3 1361 }
1362 //SetClustersOwner(kTRUE);
1363 AliTRDReconstructor::SetClusters(0x0);
6d50f529 1364 }
6d50f529 1365 return fRecPoints;
eb52b657 1366
3551db50 1367}
fc546d21 1368
a5b99acd 1369//_____________________________________________________________________________
1370TClonesArray *AliTRDclusterizer::TrackletsArray()
1371{
1372 //
1373 // Returns the list of rec points
1374 //
1375
1376 if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
94ac01ef 1377 fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
a5b99acd 1378 //SetClustersOwner(kTRUE);
1379 //AliTRDReconstructor::SetTracklets(0x0);
1380 }
1381 return fTracklets;
1382
1383}
1384