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