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