]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterizer.cxx
Bug fix in AliESDHeader::AddTriggerIR.
[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
17e0e535 670 if(fReconstructor->IsHLT()){
8dbc94c7 671 fRawStream->SetSharedPadReadout(kFALSE);
17e0e535 672 fRawStream->SetNoErrorWarning();
673 }
3fe61b77 674
17e0e535 675 AliDebug(1,Form("Stream version: %s", fRawStream->IsA()->GetName()));
3fe61b77 676
677 Int_t det = 0;
8dbc94c7 678 while ((det = fRawStream->NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
66f6bfd9 679 Bool_t iclusterBranch = kFALSE;
680 if (fDigitsManager->GetIndexes(det)->HasEntry()){
486c2339 681 iclusterBranch = MakeClusters(det);
3fe61b77 682 }
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;
703 }
dfbb4bb9 704
3d0c7d6d 705 AliInfo(Form("Number of found clusters : %d", fNoOfClusters));
66f6bfd9 706 return kTRUE;
eb52b657 707
3fe61b77 708}
709
fa7427d0 710//_____________________________________________________________________________
711UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
712{
023b669c 713 //
714 // Check if a pad is masked
715 //
716
317b5644 717 UChar_t status = 0;
718
719 if(signal>0 && TESTBIT(signal, 10)){
720 CLRBIT(signal, 10);
721 for(int ibit=0; ibit<4; ibit++){
722 if(TESTBIT(signal, 11+ibit)){
723 SETBIT(status, ibit);
724 CLRBIT(signal, 11+ibit);
725 }
726 }
727 }
728 return status;
fa7427d0 729}
730
6d9e79cc 731//_____________________________________________________________________________
e7295a3a 732void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
6d9e79cc 733 //
734 // Set the pad status into out
735 // First three bits are needed for the position encoding
736 //
064d808d 737 out |= status << 3;
6d9e79cc 738}
739
740//_____________________________________________________________________________
064d808d 741UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
6d9e79cc 742 //
743 // return the staus encoding of the corrupted pad
744 //
745 return static_cast<UChar_t>(encoding >> 3);
746}
747
748//_____________________________________________________________________________
064d808d 749Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
6d9e79cc 750 //
751 // Return the position of the corruption
752 //
753 return encoding & 7;
754}
755
3fe61b77 756//_____________________________________________________________________________
757Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
758{
759 //
760 // Generates the cluster.
761 //
762
763 // Get the digits
615f0826 764 fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
765 fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline();
f5375dcb 766
3fe61b77 767 // This is to take care of switched off super modules
b72f4eaf 768 if (!fDigits->HasData()) return kFALSE;
3fe61b77 769
064d808d 770 fIndexes = fDigitsManager->GetIndexes(det);
b72f4eaf 771 if (fIndexes->IsAllocated() == kFALSE) {
772 AliError("Indexes do not exist!");
773 return kFALSE;
774 }
064d808d 775
3fe61b77 776 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
b72f4eaf 777 if (!calibration) {
778 AliFatal("No AliTRDcalibDB instance available\n");
779 return kFALSE;
780 }
3fe61b77 781
3a039a31 782 if (!fReconstructor){
783 AliError("Reconstructor not set\n");
784 return kFALSE;
785 }
c5f589b9 786
064d808d 787 fMaxThresh = fReconstructor->GetRecoParam()->GetClusMaxThresh();
788 fSigThresh = fReconstructor->GetRecoParam()->GetClusSigThresh();
789 fMinMaxCutSigma = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
790 fMinLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
3fe61b77 791
064d808d 792 Int_t istack = fIndexes->GetStack();
793 fLayer = fIndexes->GetLayer();
794 Int_t isector = fIndexes->GetSM();
3fe61b77 795
796 // Start clustering in the chamber
797
064d808d 798 fDet = AliTRDgeometry::GetDetector(fLayer,istack,isector);
799 if (fDet != det) {
828c6f80 800 AliError("Strange Detector number Missmatch!");
b65e5048 801 return kFALSE;
802 }
3fe61b77 803
839f6b9b 804 AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
805
3fe61b77 806 // TRD space point transformation
807 fTransform->SetDetector(det);
808
064d808d 809 Int_t iGeoLayer = AliGeomManager::kTRD1 + fLayer;
053767a4 810 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
064d808d 811 fVolid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
3fe61b77 812
a5b99acd 813 if(fReconstructor->IsProcessingTracklets() && fTrackletContainer)
814 AddTrackletsToArray();
815
1b3a719f 816 fColMax = fDigits->GetNcol();
817 //Int_t nRowMax = fDigits->GetNrow();
818 fTimeTotal = fDigits->GetNtime();
3fe61b77 819
69a850a0 820 // Check consistency between OCDB and raw data
4e6823c2 821 Int_t nTimeOCDB = calibration->GetNumberOfTimeBinsDCS();
822 if ((nTimeOCDB > -1) &&
823 (fTimeTotal != nTimeOCDB)) {
828c6f80 824 AliError(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d])"
825 ,fTimeTotal,calibration->GetNumberOfTimeBinsDCS()));
69a850a0 826 }
827
3fe61b77 828 // Detector wise calibration object for the gain factors
064d808d 829 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
3fe61b77 830 // Calibration object with pad wise values for the gain factors
064d808d 831 fCalGainFactorROC = calibration->GetGainFactorROC(fDet);
3fe61b77 832 // Calibration value for chamber wise gain factor
064d808d 833 fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
0e09df31 834
df83a620 835 // Detector wise calibration object for the noise
064d808d 836 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
df83a620 837 // Calibration object with pad wise values for the noise
064d808d 838 fCalNoiseROC = calibration->GetNoiseROC(fDet);
df83a620 839 // Calibration value for chamber wise noise
064d808d 840 fCalNoiseDetValue = calNoiseDet->GetValue(fDet);
1b3a719f 841
842 // Calibration object with the pad status
843 fCalPadStatusROC = calibration->GetPadStatusROC(fDet);
6b0d4ad5 844
a2fbb6ec 845 SetBit(kLUT, fReconstructor->GetRecoParam()->UseLUT());
846 SetBit(kGAUS, fReconstructor->GetRecoParam()->UseGAUS());
b72f4eaf 847 SetBit(kHLT, fReconstructor->IsHLT());
6b0d4ad5 848
064d808d 849 firstClusterROC = -1;
850 fClusterROC = 0;
3fe61b77 851
b72f4eaf 852 // Apply the gain and the tail cancelation via digital filter
a2fbb6ec 853 if(fReconstructor->GetRecoParam()->UseTailCancelation()) TailCancelation();
023b669c 854
486c2339 855 MaxStruct curr, last;
064d808d 856 Int_t nMaximas = 0, nCorrupted = 0;
f5375dcb 857
064d808d 858 // Here the clusterfining is happening
859
615f0826 860 for(curr.time = 0; curr.time < fTimeTotal; curr.time++){
861 while(fIndexes->NextRCIndex(curr.row, curr.col)){
862 if(IsMaximum(curr, curr.padStatus, &curr.signals[0])){
863 if(last.row>-1){
864 if(curr.time==last.time && curr.row==last.row && curr.col==last.col+2) FivePadCluster(last, curr);
b72f4eaf 865 CreateCluster(last);
866 }
615f0826 867 last=curr; curr.fivePad=kFALSE;
b72f4eaf 868 }
b72f4eaf 869 }
870 }
615f0826 871 if(last.row>-1) CreateCluster(last);
df83a620 872
a2fbb6ec 873 if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
874 TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
064d808d 875 (*fDebugStream) << "MakeClusters"
b72f4eaf 876 << "Detector=" << det
877 << "NMaxima=" << nMaximas
878 << "NClusters=" << fClusterROC
879 << "NCorrupted=" << nCorrupted
880 << "\n";
df83a620 881 }
b72f4eaf 882 if (TestBit(kLabels)) AddLabels();
f5375dcb 883
064d808d 884 return kTRUE;
f5375dcb 885
064d808d 886}
3fe61b77 887
064d808d 888//_____________________________________________________________________________
1b3a719f 889Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Short_t *const Signals)
064d808d 890{
891 //
892 // Returns true if this row,col,time combination is a maximum.
893 // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
894 //
c96d03ba 895
615f0826 896 Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
897 Signals[1] = (Short_t)((fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) / gain + 0.5f);
c96d03ba 898 if(Signals[1] < fMaxThresh) return kFALSE;
3fe61b77 899
615f0826 900 Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
c96d03ba 901 if (Signals[1] < noiseMiddleThresh) return kFALSE;
f5375dcb 902
615f0826 903 if (Max.col + 1 >= fColMax || Max.col < 1) return kFALSE;
3fe61b77 904
b72f4eaf 905 UChar_t status[3]={
615f0826 906 fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
907 ,fCalPadStatusROC->GetStatus(Max.col, Max.row)
908 ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
b72f4eaf 909 };
486c2339 910
615f0826 911 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
19dd8eb9 912 Signals[0] = (Short_t)((fDigits->GetData(Max.row, Max.col-1, Max.time) - fBaseline) / gain + 0.5f);
913 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
914 Signals[2] = (Short_t)((fDigits->GetData(Max.row, Max.col+1, Max.time) - fBaseline) / gain + 0.5f);
c96d03ba 915
486c2339 916 if(!(status[0] | status[1] | status[2])) {//all pads are good
c96d03ba 917 if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) {
918 if ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh)) {
5a8db426 919 if(Signals[0]<0)Signals[0]=0;
920 if(Signals[2]<0)Signals[2]=0;
b72f4eaf 921 Float_t noiseSumThresh = fMinLeftRightCutSigma
922 * fCalNoiseDetValue
615f0826 923 * fCalNoiseROC->GetValue(Max.col, Max.row);
c96d03ba 924 if ((Signals[2]+Signals[0]+Signals[1]) < noiseSumThresh) return kFALSE;
b72f4eaf 925 padStatus = 0;
926 return kTRUE;
eb52b657 927 }
064d808d 928 }
b72f4eaf 929 } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
5a8db426 930 if(Signals[0]<0)Signals[0]=0;
931 if(Signals[2]<0)Signals[2]=0;
c96d03ba 932 if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] >= fSigThresh) {
933 Signals[2]=0;
486c2339 934 SetPadStatus(status[2], padStatus);
935 return kTRUE;
936 }
c96d03ba 937 else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] >= fSigThresh) {
938 Signals[0]=0;
486c2339 939 SetPadStatus(status[0], padStatus);
940 return kTRUE;
941 }
c96d03ba 942 else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] >= fSigThresh) || (Signals[0] >= fSigThresh))) {
615f0826 943 Signals[1] = (Short_t)(fMaxThresh + 0.5f);
486c2339 944 SetPadStatus(status[1], padStatus);
064d808d 945 return kTRUE;
946 }
947 }
948 return kFALSE;
949}
3fe61b77 950
064d808d 951//_____________________________________________________________________________
1b3a719f 952Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
064d808d 953{
954 //
955 // Look for 5 pad cluster with minimum in the middle
956 // Gives back the ratio
957 //
615f0826 958
959 if (ThisMax.col >= fColMax - 3) return kFALSE;
960 Float_t gain;
961 if (ThisMax.col < fColMax - 5){
962 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
963 if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 964 return kFALSE;
064d808d 965 }
615f0826 966 if (ThisMax.col > 1) {
967 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
968 if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain)
486c2339 969 return kFALSE;
970 }
971
486c2339 972 const Float_t kEpsilon = 0.01;
615f0826 973 Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
974 NeighbourMax.signals[1], NeighbourMax.signals[2]};
486c2339 975
976 // Unfold the two maxima and set the signal on
977 // the overlapping pad to the ratio
1b3a719f 978 Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
615f0826 979 ThisMax.signals[2] = (Short_t)(ThisMax.signals[2]*ratio + 0.5f);
980 NeighbourMax.signals[0] = (Short_t)(NeighbourMax.signals[0]*(1-ratio) + 0.5f);
981 ThisMax.fivePad=kTRUE;
982 NeighbourMax.fivePad=kTRUE;
486c2339 983 return kTRUE;
c96d03ba 984
064d808d 985}
eb52b657 986
064d808d 987//_____________________________________________________________________________
486c2339 988void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
064d808d 989{
990 //
3d0c7d6d 991 // Creates a cluster at the given position and saves it in fRecPoints
064d808d 992 //
eb52b657 993
c96d03ba 994 Int_t nPadCount = 1;
615f0826 995 Short_t signals[7] = { 0, 0, Max.signals[0], Max.signals[1], Max.signals[2], 0, 0 };
c96d03ba 996 if(!TestBit(kHLT)) CalcAdditionalInfo(Max, signals, nPadCount);
997
615f0826 998 AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
c96d03ba 999 cluster.SetNPads(nPadCount);
b72f4eaf 1000 if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1001 else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1002 else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
3fe61b77 1003
615f0826 1004 cluster.SetFivePad(Max.fivePad);
b72f4eaf 1005 // set pads status for the cluster
1006 UChar_t maskPosition = GetCorruption(Max.padStatus);
1007 if (maskPosition) {
1008 cluster.SetPadMaskedPosition(maskPosition);
1009 cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1010 }
064d808d 1011
1012 // Transform the local cluster coordinates into calibrated
1013 // space point positions defined in the local tracking system.
1014 // Here the calibration for T0, Vdrift and ExB is applied as well.
b72f4eaf 1015 if(!fTransform->Transform(&cluster)) return;
1016 // Temporarily store the Max.Row, column and time bin of the center pad
1017 // Used to later on assign the track indices
615f0826 1018 cluster.SetLabel(Max.row, 0);
1019 cluster.SetLabel(Max.col, 1);
1020 cluster.SetLabel(Max.time,2);
c96d03ba 1021
b72f4eaf 1022 //needed for HLT reconstruction
1023 AddClusterToArray(&cluster);
1024
1025 // Store the index of the first cluster in the current ROC
1026 if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1027
1028 fNoOfClusters++;
1029 fClusterROC++;
3fe61b77 1030}
1031
3d0c7d6d 1032//_____________________________________________________________________________
1033void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1034{
1035 // Look to the right
1036 Int_t ii = 1;
615f0826 1037 while (fDigits->GetData(Max.row, Max.col-ii, Max.time) >= fSigThresh) {
3d0c7d6d 1038 nPadCount++;
1039 ii++;
615f0826 1040 if (Max.col < ii) break;
3d0c7d6d 1041 }
1042 // Look to the left
1043 ii = 1;
615f0826 1044 while (fDigits->GetData(Max.row, Max.col+ii, Max.time) >= fSigThresh) {
3d0c7d6d 1045 nPadCount++;
1046 ii++;
615f0826 1047 if (Max.col+ii >= fColMax) break;
3d0c7d6d 1048 }
1049
1050 // Store the amplitudes of the pads in the cluster for later analysis
1051 // and check whether one of these pads is masked in the database
615f0826 1052 signals[2]=Max.signals[0];
1053 signals[3]=Max.signals[1];
1054 signals[4]=Max.signals[2];
1055 Float_t gain;
3d0c7d6d 1056 for(Int_t i = 0; i<2; i++)
1057 {
615f0826 1058 if(Max.col+i >= 3){
1059 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-3+i,Max.row);
1060 signals[i] = (Short_t)((fDigits->GetData(Max.row, Max.col-3+i, Max.time) - fBaseline) / gain + 0.5f);
1061 }
1062 if(Max.col+3-i < fColMax){
1063 gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+3-i,Max.row);
1064 signals[6-i] = (Short_t)((fDigits->GetData(Max.row, Max.col+3-i, Max.time) - fBaseline) / gain + 0.5f);
1065 }
3d0c7d6d 1066 }
1067 /*for (Int_t jPad = Max.Col-3; jPad <= Max.Col+3; jPad++) {
1068 if ((jPad >= 0) && (jPad < fColMax))
1069 signals[jPad-Max.Col+3] = TMath::Nint(fDigits->GetData(Max.Row,jPad,Max.Time));
1070 }*/
1071}
1072
1073//_____________________________________________________________________________
e7295a3a 1074void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
3d0c7d6d 1075{
1076 //
1077 // Add a cluster to the array
1078 //
1079
1080 Int_t n = RecPoints()->GetEntriesFast();
1081 if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1082 new((*RecPoints())[n]) AliTRDcluster(*cluster);
1083}
1084
a5b99acd 1085//_____________________________________________________________________________
1086void AliTRDclusterizer::AddTrackletsToArray()
1087{
1088 //
1089 // Add the online tracklets of this chamber to the array
1090 //
1091
1092 UInt_t* trackletword;
1093 for(Int_t side=0; side<2; side++)
1094 {
1095 Int_t trkl=0;
1096 trackletword=fTrackletContainer[side];
97be9934 1097 while(trackletword[trkl]>0){
a5b99acd 1098 Int_t n = TrackletsArray()->GetEntriesFast();
35943b1e 1099 AliTRDtrackletWord tmp(trackletword[trkl]);
1100 new((*TrackletsArray())[n]) AliTRDcluster(&tmp,fDet,fVolid);
a5b99acd 1101 trkl++;
97be9934 1102 }
a5b99acd 1103 }
1104}
1105
3fe61b77 1106//_____________________________________________________________________________
6cd9a21f 1107Bool_t AliTRDclusterizer::AddLabels()
3551db50 1108{
1109 //
3fe61b77 1110 // Add the track indices to the found clusters
3551db50 1111 //
1112
3fe61b77 1113 const Int_t kNclus = 3;
1114 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1115 const Int_t kNtrack = kNdict * kNclus;
1116
1117 Int_t iClusterROC = 0;
1118
1119 Int_t row = 0;
1120 Int_t col = 0;
1121 Int_t time = 0;
1122 Int_t iPad = 0;
1123
1124 // Temporary array to collect the track indices
6cd9a21f 1125 Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
3fe61b77 1126
1127 // Loop through the dictionary arrays one-by-one
1128 // to keep memory consumption low
b65e5048 1129 AliTRDarrayDictionary *tracksIn = 0; //mod
3fe61b77 1130 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1131
1132 // tracksIn should be expanded beforehand!
6cd9a21f 1133 tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
3fe61b77 1134
1135 // Loop though the clusters found in this ROC
6cd9a21f 1136 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1137
3fe61b77 1138 AliTRDcluster *cluster = (AliTRDcluster *)
b65e5048 1139 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
3fe61b77 1140 row = cluster->GetLabel(0);
1141 col = cluster->GetLabel(1);
1142 time = cluster->GetLabel(2);
1143
1144 for (iPad = 0; iPad < kNclus; iPad++) {
b65e5048 1145 Int_t iPadCol = col - 1 + iPad;
1146 Int_t index = tracksIn->GetData(row,iPadCol,time); //Modification of -1 in Track
1147 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
3fe61b77 1148 }
1149
1150 }
1151
1152 }
1153
1154 // Copy the track indices into the cluster
1155 // Loop though the clusters found in this ROC
6cd9a21f 1156 for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
317b5644 1157
3fe61b77 1158 AliTRDcluster *cluster = (AliTRDcluster *)
1159 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1160 cluster->SetLabel(-9999,0);
1161 cluster->SetLabel(-9999,1);
1162 cluster->SetLabel(-9999,2);
1163
1164 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1165
1166 }
1167
1168 delete [] idxTracks;
1169
1170 return kTRUE;
1171
1172}
1173
3fe61b77 1174//_____________________________________________________________________________
e7295a3a 1175Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
3fe61b77 1176{
1177 //
1178 // Method to unfold neighbouring maxima.
1179 // The charge ratio on the overlapping pad is calculated
1180 // until there is no more change within the range given by eps.
1181 // The resulting ratio is then returned to the calling method.
1182 //
1183
1184 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 1185 if (!calibration) {
3fe61b77 1186 AliError("No AliTRDcalibDB instance available\n");
1187 return kFALSE;
6d50f529 1188 }
3fe61b77 1189
1190 Int_t irc = 0;
1191 Int_t itStep = 0; // Count iteration steps
1192
1193 Double_t ratio = 0.5; // Start value for ratio
1194 Double_t prevRatio = 0.0; // Store previous ratio
1195
1196 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1197 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1198 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1199
1200 // Start the iteration
1201 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1202
1203 itStep++;
1204 prevRatio = ratio;
1205
1206 // Cluster position according to charge ratio
1207 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
064d808d 1208 / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
3fe61b77 1209 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1210 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1211
1212 // Set cluster charge ratio
064d808d 1213 irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
3fe61b77 1214 Double_t ampLeft = padSignal[1] / newSignal[1];
064d808d 1215 irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
3fe61b77 1216 Double_t ampRight = padSignal[3] / newSignal[1];
1217
1218 // Apply pad response to parameters
053767a4 1219 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1220 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
3fe61b77 1221
1222 // Calculate new overlapping ratio
1223 ratio = TMath::Min((Double_t) 1.0
1224 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1225
b43a3e17 1226 }
88719a08 1227
3fe61b77 1228 return ratio;
1229
1230}
88719a08 1231
3fe61b77 1232//_____________________________________________________________________________
064d808d 1233void AliTRDclusterizer::TailCancelation()
3fe61b77 1234{
1235 //
1236 // Applies the tail cancelation and gain factors:
1b3a719f 1237 // Transform fDigits to fDigits
3fe61b77 1238 //
1239
1240 Int_t iRow = 0;
1241 Int_t iCol = 0;
1242 Int_t iTime = 0;
1243
1b3a719f 1244 Double_t *inADC = new Double_t[fTimeTotal]; // ADC data before tail cancellation
064d808d 1245 Double_t *outADC = new Double_t[fTimeTotal]; // ADC data after tail cancellation
1b3a719f 1246
a2fbb6ec 1247 TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
cb5522da 1248 Bool_t debugStreaming = fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
064d808d 1249 while(fIndexes->NextRCIndex(iRow, iCol))
3fe61b77 1250 {
f5375dcb 1251 Bool_t corrupted = kFALSE;
5a8db426 1252 if (fCalPadStatusROC->GetStatus(iCol, iRow)) corrupted = kTRUE;
1253
cb5522da 1254 // Save data into the temporary processing array and substract the baseline,
1255 // since DeConvExp does not expect a baseline
064d808d 1256 for (iTime = 0; iTime < fTimeTotal; iTime++)
cb5522da 1257 inADC[iTime] = fDigits->GetData(iRow,iCol,iTime)-fBaseline;
5a8db426 1258
cb5522da 1259 if(debugStreaming){
1260 for (iTime = 0; iTime < fTimeTotal; iTime++)
5a8db426 1261 (*fDebugStream) << "TailCancellation"
1262 << "col=" << iCol
1263 << "row=" << iRow
1264 << "time=" << iTime
1265 << "inADC=" << inADC[iTime]
1266 << "outADC=" << outADC[iTime]
1267 << "corrupted=" << corrupted
1268 << "\n";
cb5522da 1269 }
5a8db426 1270
b65e5048 1271 if (!corrupted)
1272 {
1273 // Apply the tail cancelation via the digital filter
1274 // (only for non-coorupted pads)
1b3a719f 1275 DeConvExp(&inADC[0],&outADC[0],fTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
b65e5048 1276 }
5a8db426 1277 else memcpy(&outADC[0],&inADC[0],fTimeTotal*sizeof(inADC[0]));
3fe61b77 1278
cb5522da 1279 // Save tailcancalled data and add the baseline
5a8db426 1280 for(iTime = 0; iTime < fTimeTotal; iTime++)
cb5522da 1281 fDigits->SetData(iRow,iCol,iTime,(Short_t)(outADC[iTime] + fBaseline + 0.5));
5a8db426 1282
3fe61b77 1283 } // while irow icol
1b3a719f 1284
3fe61b77 1285 delete [] inADC;
1286 delete [] outADC;
1287
1288 return;
1289
1290}
1291
1292//_____________________________________________________________________________
064d808d 1293void AliTRDclusterizer::DeConvExp(const Double_t *const source, Double_t *const target
1294 ,const Int_t n, const Int_t nexp)
3fe61b77 1295{
1296 //
1297 // Tail cancellation by deconvolution for PASA v4 TRF
1298 //
1299
1300 Double_t rates[2];
1301 Double_t coefficients[2];
1302
1303 // Initialization (coefficient = alpha, rates = lambda)
1304 Double_t r1 = 1.0;
1305 Double_t r2 = 1.0;
1306 Double_t c1 = 0.5;
1307 Double_t c2 = 0.5;
1308
1309 if (nexp == 1) { // 1 Exponentials
1310 r1 = 1.156;
1311 r2 = 0.130;
1312 c1 = 0.066;
1313 c2 = 0.000;
1314 }
1315 if (nexp == 2) { // 2 Exponentials
181c7f7e 1316 Double_t par[4];
a2fbb6ec 1317 fReconstructor->GetRecoParam()->GetTCParams(par);
181c7f7e 1318 r1 = par[0];//1.156;
1319 r2 = par[1];//0.130;
1320 c1 = par[2];//0.114;
1321 c2 = par[3];//0.624;
3fe61b77 1322 }
1323
1324 coefficients[0] = c1;
1325 coefficients[1] = c2;
1326
1327 Double_t dt = 0.1;
1328
1329 rates[0] = TMath::Exp(-dt/(r1));
1330 rates[1] = TMath::Exp(-dt/(r2));
1331
1332 Int_t i = 0;
1333 Int_t k = 0;
1334
1335 Double_t reminder[2];
1336 Double_t correction = 0.0;
1337 Double_t result = 0.0;
1338
1339 // Attention: computation order is important
1340 for (k = 0; k < nexp; k++) {
1341 reminder[k] = 0.0;
1342 }
1343
1344 for (i = 0; i < n; i++) {
1345
1346 result = (source[i] - correction); // No rescaling
1347 target[i] = result;
1348
1349 for (k = 0; k < nexp; k++) {
1350 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1351 }
1352
1353 correction = 0.0;
1354 for (k = 0; k < nexp; k++) {
1355 correction += reminder[k];
1356 }
1357
1358 }
6d50f529 1359
1360}
1361
1362//_____________________________________________________________________________
1363void AliTRDclusterizer::ResetRecPoints()
1364{
1365 //
1366 // Resets the list of rec points
1367 //
1368
1369 if (fRecPoints) {
1370 fRecPoints->Delete();
48f8adf3 1371 delete fRecPoints;
6d50f529 1372 }
6d50f529 1373}
1374
1375//_____________________________________________________________________________
66f6bfd9 1376TClonesArray *AliTRDclusterizer::RecPoints()
6d50f529 1377{
1378 //
1379 // Returns the list of rec points
1380 //
1381
1382 if (!fRecPoints) {
48f8adf3 1383 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
8ae98148 1384 // determine number of clusters which has to be allocated
1385 Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
8ae98148 1386
1387 fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
48f8adf3 1388 }
1389 //SetClustersOwner(kTRUE);
1390 AliTRDReconstructor::SetClusters(0x0);
6d50f529 1391 }
6d50f529 1392 return fRecPoints;
eb52b657 1393
3551db50 1394}
fc546d21 1395
a5b99acd 1396//_____________________________________________________________________________
1397TClonesArray *AliTRDclusterizer::TrackletsArray()
1398{
1399 //
1400 // Returns the list of rec points
1401 //
1402
1403 if (!fTracklets && fReconstructor->IsProcessingTracklets()) {
94ac01ef 1404 fTracklets = new TClonesArray("AliTRDcluster", 2*MAXTRACKLETSPERHC);
a5b99acd 1405 //SetClustersOwner(kTRUE);
1406 //AliTRDReconstructor::SetTracklets(0x0);
1407 }
1408 return fTracklets;
1409
1410}
1411