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