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