Line data Source code
1 : /**************************************************************************
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 : **************************************************************************/
15 :
16 : /* $Id$ */
17 :
18 : ///////////////////////////////////////////////////////////////////////////////
19 : // //
20 : // TRD cluster finder //
21 : // //
22 : ///////////////////////////////////////////////////////////////////////////////
23 :
24 : #include <TClonesArray.h>
25 : #include <TObjArray.h>
26 :
27 : #include "AliRunLoader.h"
28 : #include "AliLoader.h"
29 : #include "AliTreeLoader.h"
30 : #include "AliAlignObj.h"
31 :
32 : #include "AliTRDclusterizer.h"
33 : #include "AliTRDcluster.h"
34 : #include "AliTRDReconstructor.h"
35 : #include "AliTRDgeometry.h"
36 : #include "AliTRDarrayDictionary.h"
37 : #include "AliTRDarrayADC.h"
38 : #include "AliTRDdigitsManager.h"
39 : #include "AliTRDdigitsParam.h"
40 : #include "AliTRDrawData.h"
41 : #include "AliTRDcalibDB.h"
42 : #include "AliTRDtransform.h"
43 : #include "AliTRDSignalIndex.h"
44 : #include "AliTRDrawStream.h"
45 : #include "AliTRDfeeParam.h"
46 : #include "AliTRDtrackletWord.h"
47 : #include "AliTRDtrackletMCM.h"
48 : #include "AliTRDtrackGTU.h"
49 : #include "AliESDTrdTrack.h"
50 :
51 : #include "TTreeStream.h"
52 :
53 : #include "AliTRDCalROC.h"
54 : #include "AliTRDCalDet.h"
55 : #include "AliTRDCalSingleChamberStatus.h"
56 : #include "AliTRDCalOnlineGainTableROC.h"
57 :
58 48 : ClassImp(AliTRDclusterizer)
59 :
60 : //_____________________________________________________________________________
61 : AliTRDclusterizer::AliTRDclusterizer(const AliTRDReconstructor *const rec)
62 0 : :TNamed()
63 0 : ,fReconstructor(rec)
64 0 : ,fRunLoader(NULL)
65 0 : ,fClusterTree(NULL)
66 0 : ,fRecPoints(NULL)
67 0 : ,fTracklets(NULL)
68 0 : ,fTracks(NULL)
69 0 : ,fDigitsManager(new AliTRDdigitsManager())
70 0 : ,fRawVersion(2)
71 0 : ,fTransform(new AliTRDtransform(0))
72 0 : ,fDigits(NULL)
73 0 : ,fDigitsRaw(NULL)
74 0 : ,fIndexes(NULL)
75 0 : ,fMaxThresh(0)
76 0 : ,fMaxThreshTest(0)
77 0 : ,fSigThresh(0)
78 0 : ,fMinMaxCutSigma(0)
79 0 : ,fMinLeftRightCutSigma(0)
80 0 : ,fLayer(0)
81 0 : ,fDet(0)
82 0 : ,fVolid(0)
83 0 : ,fColMax(0)
84 0 : ,fTimeTotal(0)
85 0 : ,fTimeBinsDCS(-999)
86 0 : ,fRun(-1)
87 0 : ,fCalGainFactorROC(NULL)
88 0 : ,fCalGainFactorDetValue(0)
89 0 : ,fCalNoiseROC(NULL)
90 0 : ,fCalNoiseDetValue(0)
91 0 : ,fCalPadStatusROC(NULL)
92 0 : ,fCalOnlGainROC(NULL)
93 0 : ,fClusterROC(0)
94 0 : ,firstClusterROC(0)
95 0 : ,fNoOfClusters(0)
96 0 : ,fBaseline(0)
97 0 : ,fRawStream(NULL)
98 0 : ,fTrgFlags()
99 0 : {
100 : //
101 : // AliTRDclusterizer default constructor
102 : //
103 :
104 0 : SetBit(kLabels, kTRUE);
105 0 : SetBit(knewDM, kFALSE);
106 :
107 0 : fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
108 :
109 : // Initialize debug stream
110 0 : if(fReconstructor){
111 0 : if(fReconstructor->GetRecoParam()->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 1){
112 0 : TDirectory *savedir = gDirectory;
113 : //fgGetDebugStream = new TTreeSRedirector("TRD.ClusterizerDebug.root");
114 0 : savedir->cd();
115 0 : }
116 : }
117 :
118 0 : }
119 :
120 : //_____________________________________________________________________________
121 : AliTRDclusterizer::AliTRDclusterizer(const Text_t *name
122 : , const Text_t *title
123 : , const AliTRDReconstructor *const rec)
124 2 : :TNamed(name,title)
125 2 : ,fReconstructor(rec)
126 2 : ,fRunLoader(NULL)
127 2 : ,fClusterTree(NULL)
128 2 : ,fRecPoints(NULL)
129 2 : ,fTracklets(NULL)
130 2 : ,fTracks(NULL)
131 6 : ,fDigitsManager(new AliTRDdigitsManager())
132 2 : ,fRawVersion(2)
133 6 : ,fTransform(new AliTRDtransform(0))
134 2 : ,fDigits(NULL)
135 2 : ,fDigitsRaw(NULL)
136 2 : ,fIndexes(NULL)
137 2 : ,fMaxThresh(0)
138 2 : ,fMaxThreshTest(0)
139 2 : ,fSigThresh(0)
140 2 : ,fMinMaxCutSigma(0)
141 2 : ,fMinLeftRightCutSigma(0)
142 2 : ,fLayer(0)
143 2 : ,fDet(0)
144 2 : ,fVolid(0)
145 2 : ,fColMax(0)
146 2 : ,fTimeTotal(0)
147 2 : ,fTimeBinsDCS(-999)
148 2 : ,fRun(-1)
149 2 : ,fCalGainFactorROC(NULL)
150 2 : ,fCalGainFactorDetValue(0)
151 2 : ,fCalNoiseROC(NULL)
152 2 : ,fCalNoiseDetValue(0)
153 2 : ,fCalPadStatusROC(NULL)
154 2 : ,fCalOnlGainROC(NULL)
155 2 : ,fClusterROC(0)
156 2 : ,firstClusterROC(0)
157 2 : ,fNoOfClusters(0)
158 2 : ,fBaseline(0)
159 2 : ,fRawStream(NULL)
160 2 : ,fTrgFlags()
161 10 : {
162 : //
163 : // AliTRDclusterizer constructor
164 : //
165 :
166 2 : SetBit(kLabels, kTRUE);
167 2 : SetBit(knewDM, kFALSE);
168 :
169 2 : fDigitsManager->CreateArrays();
170 :
171 4 : fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
172 :
173 : //FillLUT();
174 :
175 4 : }
176 :
177 : //_____________________________________________________________________________
178 : AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
179 0 : :TNamed(c)
180 0 : ,fReconstructor(c.fReconstructor)
181 0 : ,fRunLoader(NULL)
182 0 : ,fClusterTree(NULL)
183 0 : ,fRecPoints(NULL)
184 0 : ,fTracklets(NULL)
185 0 : ,fTracks(NULL)
186 0 : ,fDigitsManager(NULL)
187 0 : ,fRawVersion(2)
188 0 : ,fTransform(NULL)
189 0 : ,fDigits(NULL)
190 0 : ,fDigitsRaw(NULL)
191 0 : ,fIndexes(NULL)
192 0 : ,fMaxThresh(0)
193 0 : ,fMaxThreshTest(0)
194 0 : ,fSigThresh(0)
195 0 : ,fMinMaxCutSigma(0)
196 0 : ,fMinLeftRightCutSigma(0)
197 0 : ,fLayer(0)
198 0 : ,fDet(0)
199 0 : ,fVolid(0)
200 0 : ,fColMax(0)
201 0 : ,fTimeTotal(0)
202 0 : ,fTimeBinsDCS(-999)
203 0 : ,fRun(-1)
204 0 : ,fCalGainFactorROC(NULL)
205 0 : ,fCalGainFactorDetValue(0)
206 0 : ,fCalNoiseROC(NULL)
207 0 : ,fCalNoiseDetValue(0)
208 0 : ,fCalPadStatusROC(NULL)
209 0 : ,fCalOnlGainROC(NULL)
210 0 : ,fClusterROC(0)
211 0 : ,firstClusterROC(0)
212 0 : ,fNoOfClusters(0)
213 0 : ,fBaseline(0)
214 0 : ,fRawStream(NULL)
215 0 : ,fTrgFlags()
216 0 : {
217 : //
218 : // AliTRDclusterizer copy constructor
219 : //
220 :
221 0 : SetBit(kLabels, kTRUE);
222 0 : SetBit(knewDM, kFALSE);
223 :
224 : //FillLUT();
225 :
226 0 : }
227 :
228 : //_____________________________________________________________________________
229 : AliTRDclusterizer::~AliTRDclusterizer()
230 12 : {
231 : //
232 : // AliTRDclusterizer destructor
233 : //
234 :
235 2 : if (fDigitsManager) {
236 4 : delete fDigitsManager;
237 2 : fDigitsManager = NULL;
238 2 : }
239 :
240 2 : if (fDigitsRaw) {
241 4 : delete fDigitsRaw;
242 2 : fDigitsRaw = NULL;
243 2 : }
244 :
245 2 : if (fTransform){
246 4 : delete fTransform;
247 2 : fTransform = NULL;
248 2 : }
249 :
250 2 : if (fRawStream){
251 2 : delete fRawStream;
252 1 : fRawStream = NULL;
253 1 : }
254 6 : }
255 :
256 : //_____________________________________________________________________________
257 : AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
258 : {
259 : //
260 : // Assignment operator
261 : //
262 :
263 0 : if (this != &c)
264 : {
265 0 : ((AliTRDclusterizer &) c).Copy(*this);
266 0 : }
267 :
268 0 : return *this;
269 :
270 : }
271 :
272 : //_____________________________________________________________________________
273 : void AliTRDclusterizer::Copy(TObject &c) const
274 : {
275 : //
276 : // Copy function
277 : //
278 :
279 0 : ((AliTRDclusterizer &) c).fClusterTree = NULL;
280 0 : ((AliTRDclusterizer &) c).fRecPoints = NULL;
281 0 : ((AliTRDclusterizer &) c).fDigitsManager = NULL;
282 0 : ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
283 0 : ((AliTRDclusterizer &) c).fTransform = NULL;
284 0 : ((AliTRDclusterizer &) c).fDigits = NULL;
285 0 : ((AliTRDclusterizer &) c).fDigitsRaw = NULL;
286 0 : ((AliTRDclusterizer &) c).fIndexes = NULL;
287 0 : ((AliTRDclusterizer &) c).fMaxThresh = 0;
288 0 : ((AliTRDclusterizer &) c).fMaxThreshTest = 0;
289 0 : ((AliTRDclusterizer &) c).fSigThresh = 0;
290 0 : ((AliTRDclusterizer &) c).fMinMaxCutSigma= 0;
291 0 : ((AliTRDclusterizer &) c).fMinLeftRightCutSigma = 0;
292 0 : ((AliTRDclusterizer &) c).fLayer = 0;
293 0 : ((AliTRDclusterizer &) c).fDet = 0;
294 0 : ((AliTRDclusterizer &) c).fVolid = 0;
295 0 : ((AliTRDclusterizer &) c).fColMax = 0;
296 0 : ((AliTRDclusterizer &) c).fTimeTotal = 0;
297 0 : ((AliTRDclusterizer &) c).fTimeBinsDCS = -999;
298 0 : ((AliTRDclusterizer &) c).fRun = -1;
299 0 : ((AliTRDclusterizer &) c).fCalGainFactorROC = NULL;
300 0 : ((AliTRDclusterizer &) c).fCalGainFactorDetValue = 0;
301 0 : ((AliTRDclusterizer &) c).fCalNoiseROC = NULL;
302 0 : ((AliTRDclusterizer &) c).fCalNoiseDetValue = 0;
303 0 : ((AliTRDclusterizer &) c).fCalPadStatusROC = NULL;
304 0 : ((AliTRDclusterizer &) c).fClusterROC = 0;
305 0 : ((AliTRDclusterizer &) c).firstClusterROC= 0;
306 0 : ((AliTRDclusterizer &) c).fNoOfClusters = 0;
307 0 : ((AliTRDclusterizer &) c).fBaseline = 0;
308 0 : ((AliTRDclusterizer &) c).fRawStream = NULL;
309 :
310 0 : }
311 :
312 : //_____________________________________________________________________________
313 : Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
314 : {
315 : //
316 : // Opens the AliROOT file. Output and input are in the same file
317 : //
318 :
319 0 : TString evfoldname = AliConfig::GetDefaultEventFolderName();
320 0 : fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
321 :
322 0 : if (!fRunLoader) {
323 0 : fRunLoader = AliRunLoader::Open(name);
324 0 : }
325 :
326 0 : if (!fRunLoader) {
327 0 : AliError(Form("Can not open session for file %s.",name));
328 0 : return kFALSE;
329 : }
330 :
331 0 : OpenInput(nEvent);
332 0 : OpenOutput();
333 :
334 0 : return kTRUE;
335 :
336 0 : }
337 :
338 : //_____________________________________________________________________________
339 : Bool_t AliTRDclusterizer::OpenOutput()
340 : {
341 : //
342 : // Open the output file
343 : //
344 :
345 0 : if (!fReconstructor->IsWritingClusters()) return kTRUE;
346 :
347 0 : TObjArray *ioArray = 0x0;
348 :
349 0 : AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
350 0 : loader->MakeTree("R");
351 :
352 0 : fClusterTree = loader->TreeR();
353 0 : fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
354 :
355 : return kTRUE;
356 :
357 0 : }
358 :
359 : //_____________________________________________________________________________
360 : Bool_t AliTRDclusterizer::OpenOutput(TTree *const clusterTree)
361 : {
362 : //
363 : // Connect the output tree
364 : //
365 :
366 : // clusters writing
367 16 : if (fReconstructor->IsWritingClusters()){
368 8 : TObjArray *ioArray = 0x0;
369 8 : fClusterTree = clusterTree;
370 8 : fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
371 8 : }
372 8 : return kTRUE;
373 : }
374 :
375 : //_____________________________________________________________________________
376 : Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
377 : {
378 : //
379 : // Opens a ROOT-file with TRD-hits and reads in the digits-tree
380 : //
381 :
382 : // Import the Trees for the event nEvent in the file
383 0 : fRunLoader->GetEvent(nEvent);
384 :
385 0 : return kTRUE;
386 :
387 : }
388 :
389 : //_____________________________________________________________________________
390 : Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
391 : {
392 : //
393 : // Fills TRDcluster branch in the tree with the clusters
394 : // found in detector = det. For det=-1 writes the tree.
395 : //
396 :
397 24 : if ((det < -1) ||
398 8 : (det >= AliTRDgeometry::Ndet())) {
399 0 : AliError(Form("Unexpected detector index %d.\n",det));
400 0 : return kFALSE;
401 : }
402 8 : Int_t nRecPoints = RecPoints()->GetEntriesFast();
403 8 : if(!nRecPoints) return kTRUE;
404 :
405 8 : const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
406 8 : TObjArray arrClTmp(400), *ioArray = &arrClTmp; // RS don't trash the heap
407 : //
408 8 : TBranch *branch = fClusterTree->GetBranch("TRDcluster");
409 8 : if (!branch) {
410 0 : fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
411 8 : } else branch->SetAddress(&ioArray);
412 :
413 : AliTRDcluster *c(NULL);
414 8 : if(det >= 0){
415 0 : for (Int_t i = 0; i < nRecPoints; i++) {
416 0 : if(!(c = (AliTRDcluster *) RecPoints()->UncheckedAt(i))) continue;
417 0 : if(det != c->GetDetector()) continue;
418 0 : ioArray->AddLast(c);
419 : }
420 0 : fClusterTree->Fill();
421 0 : ioArray->Clear();
422 : }
423 : else {
424 16 : if(!(c = (AliTRDcluster*)RecPoints()->UncheckedAt(0))){
425 0 : AliError("Missing first cluster.");
426 0 : return kFALSE;
427 : }
428 8 : Int_t detOld(c->GetDetector()), nw(0);
429 8 : ioArray->AddLast(c);
430 36904 : for (Int_t i(1); i<nRecPoints; i++) {
431 36888 : if(!(c = (AliTRDcluster *) RecPoints()->UncheckedAt(i))) continue;
432 : // Check on total cluster charge
433 18444 : if (c->GetQ() < recoParam->GetClusterQmin()) continue;
434 18444 : if(c->GetDetector() != detOld){
435 682 : nw += ioArray->GetEntriesFast();
436 : // fill & clear previously detector set of clusters
437 341 : fClusterTree->Fill();
438 341 : ioArray->Clear();
439 341 : detOld = c->GetDetector();
440 341 : }
441 18444 : ioArray->AddLast(c);
442 : }
443 16 : if(ioArray->GetEntriesFast()){
444 16 : nw += ioArray->GetEntriesFast();
445 : // fill & clear last detector set of clusters (if any)
446 8 : fClusterTree->Fill();
447 8 : ioArray->Clear();
448 : }
449 40 : AliDebug(2, Form("Clusters FOUND[%d] WRITTEN[%d] STATUS[%s]", nRecPoints, nw, nw==nRecPoints?"OK":"FAILED"));
450 8 : if(nw!=nRecPoints) AliWarning(Form("Clusters FOUND[%d] WRITTEN[%d]", nRecPoints, nw));
451 : }
452 :
453 8 : return kTRUE;
454 16 : }
455 :
456 : //_____________________________________________________________________________
457 : Bool_t AliTRDclusterizer::ReadDigits()
458 : {
459 : //
460 : // Reads the digits arrays from the input aliroot file
461 : //
462 :
463 0 : if (!fRunLoader) {
464 0 : AliError("No run loader available");
465 0 : return kFALSE;
466 : }
467 :
468 0 : AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
469 0 : if (!loader->TreeD()) {
470 0 : loader->LoadDigits();
471 0 : }
472 :
473 : // Read in the digit arrays
474 0 : return (fDigitsManager->ReadDigits(loader->TreeD()));
475 :
476 0 : }
477 :
478 : //_____________________________________________________________________________
479 : Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
480 : {
481 : //
482 : // Reads the digits arrays from the input tree
483 : //
484 :
485 : // Read in the digit arrays
486 8 : return (fDigitsManager->ReadDigits(digitsTree));
487 :
488 : }
489 :
490 : //_____________________________________________________________________________
491 : Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
492 : {
493 : //
494 : // Reads the digits arrays from the ddl file
495 : //
496 :
497 0 : AliTRDrawData raw;
498 0 : fDigitsManager = raw.Raw2Digits(rawReader);
499 :
500 : return kTRUE;
501 :
502 0 : }
503 :
504 : Bool_t AliTRDclusterizer::ReadTracklets()
505 : {
506 : //
507 : // Reads simulated tracklets from the input aliroot file
508 : //
509 :
510 8 : AliRunLoader *runLoader = AliRunLoader::Instance();
511 4 : if (!runLoader) {
512 0 : AliError("No run loader available");
513 0 : return kFALSE;
514 : }
515 :
516 4 : AliLoader* loader = runLoader->GetLoader("TRDLoader");
517 :
518 4 : AliDataLoader *trackletLoader = loader->GetDataLoader("tracklets");
519 4 : if (!trackletLoader) {
520 0 : return kFALSE;
521 : }
522 4 : trackletLoader->Load();
523 :
524 : Bool_t loaded = kFALSE;
525 : // look for simulated tracklets
526 4 : TTree *trackletTree = trackletLoader->Tree();
527 :
528 4 : if (trackletTree) {
529 4 : TBranch *trklbranch = trackletTree->GetBranch("mcmtrklbranch");
530 8 : TClonesArray *trklArray = TrackletsArray("AliTRDtrackletMCM");
531 4 : if (trklbranch && trklArray) {
532 4 : AliTRDtrackletMCM *trkl = 0x0;
533 4 : trklbranch->SetAddress(&trkl);
534 4 : Int_t nTracklets = trklbranch->GetEntries();
535 728 : for (Int_t iTracklet = 0; iTracklet < nTracklets; iTracklet++) {
536 360 : trklbranch->GetEntry(iTracklet);
537 360 : new ((*trklArray)[trklArray->GetEntries()]) AliTRDtrackletMCM(*trkl);
538 : }
539 : loaded = kTRUE;
540 4 : }
541 4 : }
542 : else {
543 : // if no simulated tracklets found, look for raw tracklets
544 0 : AliTreeLoader *treeLoader = (AliTreeLoader*) trackletLoader->GetBaseLoader("tracklets-raw");
545 0 : trackletTree = treeLoader ? treeLoader->Load(), treeLoader->Tree() : 0x0;
546 :
547 0 : if (trackletTree) {
548 0 : TClonesArray *trklArray = TrackletsArray("AliTRDtrackletWord");
549 :
550 0 : Int_t hc;
551 0 : TClonesArray *ar = 0x0;
552 0 : trackletTree->SetBranchAddress("hc", &hc);
553 0 : trackletTree->SetBranchAddress("trkl", &ar);
554 :
555 0 : Int_t nEntries = trackletTree->GetEntries();
556 0 : for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
557 0 : trackletTree->GetEntry(iEntry);
558 0 : Int_t nTracklets = ar->GetEntriesFast();
559 0 : AliDebug(2, Form("%i tracklets in HC %i", nTracklets, hc));
560 0 : for (Int_t iTracklet = 0; iTracklet < nTracklets; iTracklet++) {
561 0 : AliTRDtrackletWord *trklWord = (AliTRDtrackletWord*) (*ar)[iTracklet];
562 0 : new ((*trklArray)[trklArray->GetEntries()]) AliTRDtrackletWord(trklWord->GetTrackletWord(), hc);
563 : }
564 : }
565 : loaded = kTRUE;
566 0 : }
567 : }
568 :
569 4 : trackletLoader->UnloadAll();
570 4 : trackletLoader->CloseFile();
571 :
572 4 : return loaded;
573 4 : }
574 :
575 : Bool_t AliTRDclusterizer::ReadTracks()
576 : {
577 : //
578 : // Reads simulated GTU tracks from the input aliroot file
579 : //
580 :
581 8 : AliRunLoader *runLoader = AliRunLoader::Instance();
582 :
583 4 : if (!runLoader) {
584 0 : AliError("No run loader available");
585 0 : return kFALSE;
586 : }
587 :
588 4 : AliLoader* loader = runLoader->GetLoader("TRDLoader");
589 4 : if (!loader) {
590 0 : return kFALSE;
591 : }
592 :
593 4 : AliDataLoader *trackLoader = loader->GetDataLoader("gtutracks");
594 4 : if (!trackLoader) {
595 0 : return kFALSE;
596 : }
597 :
598 : Bool_t loaded = kFALSE;
599 :
600 4 : trackLoader->Load();
601 :
602 4 : TTree *trackTree = trackLoader->Tree();
603 4 : if (trackTree) {
604 4 : TClonesArray *trackArray = TracksArray();
605 4 : AliTRDtrackGTU *trk = 0x0;
606 4 : trackTree->SetBranchAddress("TRDtrackGTU", &trk);
607 38 : for (Int_t iTrack = 0; iTrack < trackTree->GetEntries(); iTrack++) {
608 15 : trackTree->GetEntry(iTrack);
609 30 : new ((*trackArray)[trackArray->GetEntries()]) AliESDTrdTrack(*(trk->CreateTrdTrack()));
610 : }
611 : loaded = kTRUE;
612 4 : }
613 :
614 4 : trackLoader->UnloadAll();
615 4 : trackLoader->CloseFile();
616 :
617 4 : return loaded;
618 4 : }
619 :
620 : //_____________________________________________________________________________
621 : Bool_t AliTRDclusterizer::MakeClusters()
622 : {
623 : //
624 : // Creates clusters from digits
625 : //
626 :
627 : // Propagate info from the digits manager
628 8 : if (TestBit(kLabels)){
629 4 : SetBit(kLabels, fDigitsManager->UsesDictionaries());
630 4 : }
631 :
632 : Bool_t fReturn = kTRUE;
633 4328 : for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
634 :
635 2160 : AliTRDarrayADC *digitsIn = (AliTRDarrayADC*) fDigitsManager->GetDigits(i); //mod
636 : // This is to take care of switched off super modules
637 4132 : if (!digitsIn->HasData()) continue;
638 188 : digitsIn->Expand();
639 188 : digitsIn->DeleteNegatives(); // Restore digits array to >=0 values
640 188 : AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
641 : // if (indexes->IsAllocated() == kFALSE){ // A.B.
642 188 : fDigitsManager->BuildIndexes(i);
643 : // }
644 :
645 : Bool_t fR(kFALSE);
646 188 : if (indexes->HasEntry()){
647 188 : if (TestBit(kLabels)){
648 : Int_t nDict(0);
649 1504 : for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
650 : AliTRDarrayDictionary *tracksIn(NULL); //mod
651 564 : tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(i,iDict); //mod
652 : // This is to take care of data reconstruction
653 564 : if (!tracksIn->GetDim()) continue;
654 564 : tracksIn->Expand(); nDict++;
655 564 : }
656 188 : if(!nDict){
657 0 : AliDebug(1, "MC labels not available. Switch them off.");
658 0 : SetUseLabels(kFALSE);
659 0 : }
660 188 : }
661 188 : fR = MakeClusters(i);
662 188 : fReturn = fR && fReturn;
663 188 : }
664 :
665 : //if (fR == kFALSE){
666 : // if(IsWritingClusters()) WriteClusters(i);
667 : // ResetRecPoints();
668 : //}
669 :
670 : // Clear arrays of this chamber, to prepare for next event
671 188 : fDigitsManager->ClearArrays(i);
672 188 : }
673 :
674 8 : if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
675 :
676 68 : AliInfo(Form("Found :: clusters[%d] tracklets[%d] tracks[%d]",
677 : RecPoints()?RecPoints()->GetEntriesFast():0,
678 : TrackletsArray()?TrackletsArray()->GetEntriesFast():0,
679 : TracksArray()?TracksArray()->GetEntriesFast():0));
680 :
681 4 : return fReturn;
682 :
683 0 : }
684 :
685 : //_____________________________________________________________________________
686 : Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
687 : {
688 : //
689 : // Creates clusters from raw data
690 : //
691 0 : return Raw2ClustersChamber(rawReader);
692 :
693 : }
694 :
695 : //_____________________________________________________________________________
696 : Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
697 : {
698 : //
699 : // Creates clusters from raw data
700 : //
701 : // if first call or run number was changed, update cached value ot timebinsDCS
702 8 : AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance();
703 4 : if (!calibration) {
704 0 : AliFatal("No AliTRDcalibDB instance available\n");
705 0 : return kFALSE;
706 : }
707 :
708 7 : if (fTimeBinsDCS==-999 || fRun!=(int)calibration->GetRun()) {
709 1 : fRun = calibration->GetRun();
710 1 : fTimeBinsDCS = calibration->GetNumberOfTimeBinsDCS();
711 5 : AliInfoF("Set number of DCS time bins to %d for run %d",fTimeBinsDCS,fRun);
712 1 : }
713 :
714 : // Create the digits manager
715 4 : if (!fDigitsManager){
716 1 : SetBit(knewDM, kTRUE);
717 2 : fDigitsManager = new AliTRDdigitsManager(kTRUE);
718 1 : fDigitsManager->CreateArrays();
719 1 : }
720 :
721 4 : fDigitsManager->SetUseDictionaries(TestBit(kLabels));
722 :
723 4 : if(!fRawStream)
724 4 : fRawStream = new AliTRDrawStream(rawReader);
725 : else
726 2 : fRawStream->SetReader(rawReader);
727 :
728 : //if(fReconstructor->IsHLT()){
729 4 : fRawStream->DisableErrorStorage();
730 : //}
731 :
732 : // register tracklet array for output
733 8 : fRawStream->SetTrackletArray(TrackletsArray("AliTRDtrackletWord"));
734 4 : fRawStream->SetTrackArray(TracksArray());
735 :
736 : UInt_t det = 0;
737 2168 : while ((det = fRawStream->NextChamber(fDigitsManager)) < AliTRDgeometry::kNdet){
738 2160 : if (fDigitsManager->GetIndexes(det)->HasEntry()){
739 188 : MakeClusters(det);
740 188 : fDigitsManager->ClearArrays(det);
741 188 : }
742 : }
743 :
744 152 : for (Int_t iSector = 0; iSector < AliTRDgeometry::kNsector; iSector++) {
745 72 : fTrgFlags[iSector] = fRawStream->GetTriggerFlags(iSector);
746 : }
747 :
748 8 : if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
749 :
750 4 : if(!TestBit(knewDM)){
751 2 : delete fDigitsManager;
752 1 : fDigitsManager = NULL;
753 2 : delete fRawStream;
754 1 : fRawStream = NULL;
755 1 : }
756 :
757 68 : AliInfo(Form("Found :: clusters[%d] tracklets[%d] tracks[%d]",
758 : RecPoints()?RecPoints()->GetEntriesFast():0,
759 : TrackletsArray()?TrackletsArray()->GetEntriesFast():0,
760 : TracksArray()?TracksArray()->GetEntriesFast():0));
761 : return kTRUE;
762 :
763 4 : }
764 :
765 : //_____________________________________________________________________________
766 : UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
767 : {
768 : //
769 : // Check if a pad is masked
770 : //
771 :
772 : UChar_t status = 0;
773 :
774 0 : if(signal>0 && TESTBIT(signal, 10)){
775 0 : CLRBIT(signal, 10);
776 0 : for(int ibit=0; ibit<4; ibit++){
777 0 : if(TESTBIT(signal, 11+ibit)){
778 0 : SETBIT(status, ibit);
779 0 : CLRBIT(signal, 11+ibit);
780 0 : }
781 : }
782 0 : }
783 0 : return status;
784 : }
785 :
786 : //_____________________________________________________________________________
787 : void AliTRDclusterizer::SetPadStatus(const UChar_t status, UChar_t &out) const {
788 : //
789 : // Set the pad status into out
790 : // First three bits are needed for the position encoding
791 : //
792 0 : out |= status << 3;
793 0 : }
794 :
795 : //_____________________________________________________________________________
796 : UChar_t AliTRDclusterizer::GetPadStatus(UChar_t encoding) const {
797 : //
798 : // return the staus encoding of the corrupted pad
799 : //
800 0 : return static_cast<UChar_t>(encoding >> 3);
801 : }
802 :
803 : //_____________________________________________________________________________
804 : Int_t AliTRDclusterizer::GetCorruption(UChar_t encoding) const {
805 : //
806 : // Return the position of the corruption
807 : //
808 38116 : return encoding & 7;
809 : }
810 :
811 : //_____________________________________________________________________________
812 : Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
813 : {
814 : //
815 : // Generates the cluster.
816 : //
817 :
818 : // Get the digits
819 376 : fDigits = (AliTRDarrayADC *) fDigitsManager->GetDigits(det); //mod
820 376 : fBaseline = fDigitsManager->GetDigitsParam()->GetADCbaseline(det);
821 :
822 : // This is to take care of switched off super modules
823 376 : if (!fDigits->HasData()) return kFALSE;
824 :
825 376 : fIndexes = fDigitsManager->GetIndexes(det);
826 376 : if (fIndexes->IsAllocated() == kFALSE) {
827 0 : AliError("Indexes do not exist!");
828 0 : return kFALSE;
829 : }
830 :
831 376 : AliTRDcalibDB* const calibration = AliTRDcalibDB::Instance();
832 376 : if (!calibration) {
833 0 : AliFatal("No AliTRDcalibDB instance available\n");
834 0 : return kFALSE;
835 : }
836 :
837 751 : if (fTimeBinsDCS==-999 || fRun!=(int)calibration->GetRun()) {
838 1 : fRun = calibration->GetRun();
839 1 : fTimeBinsDCS = calibration->GetNumberOfTimeBinsDCS();
840 5 : AliInfoF("Set number of DCS time bins to %d for run %d",fTimeBinsDCS,fRun);
841 1 : }
842 :
843 :
844 376 : if (!fReconstructor){
845 0 : AliError("Reconstructor not set\n");
846 0 : return kFALSE;
847 : }
848 :
849 376 : const AliTRDrecoParam *const recoParam = fReconstructor->GetRecoParam();
850 :
851 376 : fMaxThresh = recoParam->GetClusMaxThresh();
852 376 : fMaxThreshTest = (recoParam->GetClusMaxThresh()/2+fBaseline);
853 376 : fSigThresh = recoParam->GetClusSigThresh();
854 376 : fMinMaxCutSigma = recoParam->GetMinMaxCutSigma();
855 376 : fMinLeftRightCutSigma = recoParam->GetMinLeftRightCutSigma();
856 376 : const Int_t iEveryNTB = recoParam->GetRecEveryNTB();
857 :
858 376 : Int_t istack = fIndexes->GetStack();
859 376 : fLayer = fIndexes->GetLayer();
860 376 : Int_t isector = fIndexes->GetSM();
861 :
862 : // Start clustering in the chamber
863 :
864 376 : fDet = AliTRDgeometry::GetDetector(fLayer,istack,isector);
865 376 : if (fDet != det) {
866 0 : AliError(Form("Detector number missmatch! Request[%03d] RAW[%03d]", det, fDet));
867 0 : return kFALSE;
868 : }
869 :
870 1128 : AliDebug(2, Form("Det[%d] @ Sec[%d] Stk[%d] Ly[%d]", fDet, isector, istack, fLayer));
871 :
872 : // TRD space point transformation
873 376 : fTransform->SetDetector(det);
874 :
875 376 : Int_t iGeoLayer = AliGeomManager::kTRD1 + fLayer;
876 376 : Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
877 376 : fVolid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
878 :
879 376 : fColMax = fDigits->GetNcol();
880 376 : fTimeTotal = fDigitsManager->GetDigitsParam()->GetNTimeBins(det);
881 :
882 : // Check consistency between Geometry and raw data
883 376 : AliTRDpadPlane *pp(fTransform->GetPadPlane());
884 376 : Int_t ncols(pp->GetNcols()), nrows(pp->GetNrows());
885 376 : if(ncols != fColMax) AliDebug(1, Form("N cols missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, ncols, fColMax));
886 682 : if(nrows != fDigits->GetNrow()) AliDebug(1, Form("N rows missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, nrows, fDigits->GetNrow()));
887 376 : if(ncols != fIndexes->GetNcol()) AliDebug(1, Form("N cols missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, ncols, fIndexes->GetNcol()));
888 682 : if(nrows != fIndexes->GetNrow()) AliDebug(1, Form("N rows missmatch in Digits for Det[%3d] :: Geom[%3d] RAW[%3d]", fDet, nrows, fIndexes->GetNrow()));
889 :
890 : // Check consistency between OCDB and raw data
891 752 : if(fReconstructor->IsHLT()){
892 376 : if((fTimeBinsDCS > -1) && (fTimeTotal != fTimeBinsDCS)){
893 0 : AliWarning(Form("Number of timebins does not match OCDB value (RAW[%d] OCDB[%d]), using raw value"
894 : ,fTimeTotal,fTimeBinsDCS));
895 0 : }
896 : }else{
897 376 : if(fTimeBinsDCS == -1){
898 0 : AliDebug(1, "Undefined number of timebins in OCDB, using value from raw data.");
899 0 : if(!(fTimeTotal>0)){
900 0 : AliError(Form("Number of timebins in raw data is negative, skipping chamber[%3d]!", fDet));
901 0 : return kFALSE;
902 : }
903 376 : }else if(fTimeBinsDCS == -2){
904 0 : AliError("Mixed number of timebins in OCDB, no reconstruction of TRD data!");
905 0 : return kFALSE;
906 376 : }else if(fTimeTotal != fTimeBinsDCS){
907 0 : AliError(Form("Number of timebins in raw data does not match OCDB value (RAW[%d] OCDB[%d]), skipping chamber[%3d]!"
908 : ,fTimeTotal,fTimeBinsDCS, fDet));
909 0 : return kFALSE;
910 : }
911 : }
912 1128 : AliDebug(1, Form("Using %2d number of timebins for Det[%03d].", fTimeTotal, fDet));
913 :
914 : // Detector wise calibration object for the gain factors
915 376 : const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
916 : // Calibration object with pad wise values for the gain factors
917 376 : fCalGainFactorROC = calibration->GetGainFactorROC(fDet);
918 : // Calibration value for chamber wise gain factor
919 376 : fCalGainFactorDetValue = calGainFactorDet->GetValue(fDet);
920 :
921 : // Detector wise calibration object for the noise
922 376 : const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
923 : // Calibration object with pad wise values for the noise
924 376 : fCalNoiseROC = calibration->GetNoiseROC(fDet);
925 : // Calibration value for chamber wise noise
926 376 : fCalNoiseDetValue = calNoiseDet->GetValue(fDet);
927 :
928 : // Calibration object with the pad status
929 376 : fCalPadStatusROC = calibration->GetPadStatusROC(fDet);
930 : // Calibration object of the online gain
931 376 : fCalOnlGainROC = 0x0;
932 376 : if (calibration->HasOnlineFilterGain()) {
933 0 : fCalOnlGainROC = calibration->GetOnlineGainTableROC(fDet);
934 0 : }
935 :
936 376 : firstClusterROC = -1;
937 376 : fClusterROC = 0;
938 :
939 376 : SetBit(kLUT, recoParam->UseLUT());
940 376 : SetBit(kGAUS, recoParam->UseGAUS());
941 :
942 : // Apply the gain and the tail cancelation via digital filter
943 : // Use the configuration from the DCS to find out whether online
944 : // tail cancellation was applied
945 752 : if ((!calibration->HasOnlineTailCancellation()) &&
946 376 : (recoParam->UseTailCancelation())) {
947 : // save a copy of raw data
948 376 : if(TestBit(kRawSignal)){
949 376 : if(fDigitsRaw){
950 374 : fDigitsRaw->~AliTRDarrayADC();
951 374 : new(fDigitsRaw) AliTRDarrayADC(*fDigits);
952 4 : } else fDigitsRaw = new AliTRDarrayADC(*fDigits);
953 : }
954 376 : TailCancelation(recoParam);
955 376 : }
956 :
957 376 : MaxStruct curr, last;
958 376 : Int_t nMaximas = 0, nCorrupted = 0;
959 :
960 : // Here the clusterfining is happening
961 :
962 21056 : for(curr.time = 0; curr.time < fTimeTotal; curr.time+=iEveryNTB){
963 190134 : while(fIndexes->NextRCIndex(curr.row, curr.col)){
964 243025 : if(fDigits->GetData(curr.row, curr.col, curr.time) > fMaxThreshTest && IsMaximum(curr, curr.padStatus, &curr.signals[0])){
965 19058 : if(last.row>-1){
966 20459 : if(curr.col==last.col+2 && curr.row==last.row && curr.time==last.time) FivePadCluster(last, curr);
967 18682 : CreateCluster(last);
968 18682 : }
969 19058 : last=curr; curr.fivePad=kFALSE;
970 19058 : }
971 : }
972 : }
973 752 : if(last.row>-1) CreateCluster(last);
974 :
975 376 : if(recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 2 && fReconstructor->IsDebugStreaming()){
976 0 : TTreeSRedirector* fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
977 0 : (*fDebugStream) << "MakeClusters"
978 0 : << "Detector=" << det
979 0 : << "NMaxima=" << nMaximas
980 0 : << "NClusters=" << fClusterROC
981 0 : << "NCorrupted=" << nCorrupted
982 0 : << "\n";
983 0 : }
984 564 : if (TestBit(kLabels)) AddLabels();
985 :
986 : return kTRUE;
987 :
988 752 : }
989 :
990 : //_____________________________________________________________________________
991 : Bool_t AliTRDclusterizer::IsMaximum(const MaxStruct &Max, UChar_t &padStatus, Float_t *const Signals)
992 : {
993 : //
994 : // Returns true if this row,col,time combination is a maximum.
995 : // Gives back the padStatus and the signals of the center pad and the two neighbouring pads.
996 : //
997 :
998 126086 : Float_t gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col,Max.row);
999 126086 : Float_t onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col) : 1;
1000 63043 : Signals[1] = (fDigits->GetData(Max.row, Max.col, Max.time) - fBaseline) /(onlcf * gain) + 0.5f;
1001 72601 : if(Signals[1] <= fMaxThresh) return kFALSE;
1002 :
1003 107112 : if(Max.col < 1 || Max.col + 1 >= fColMax) return kFALSE;
1004 :
1005 53242 : Float_t noiseMiddleThresh = fMinMaxCutSigma*fCalNoiseDetValue*fCalNoiseROC->GetValue(Max.col, Max.row);
1006 54838 : if (Signals[1] <= noiseMiddleThresh) return kFALSE;
1007 :
1008 : Char_t status[3]={
1009 51646 : fCalPadStatusROC->GetStatus(Max.col-1, Max.row)
1010 51646 : ,fCalPadStatusROC->GetStatus(Max.col, Max.row)
1011 51646 : ,fCalPadStatusROC->GetStatus(Max.col+1, Max.row)
1012 : };
1013 :
1014 : Short_t signal(0);
1015 51646 : if((signal = fDigits->GetData(Max.row, Max.col-1, Max.time))){
1016 51616 : gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col-1,Max.row);
1017 103232 : onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col-1) : 1;
1018 51616 : Signals[0] = (signal - fBaseline) /( onlcf * gain) + 0.5f;
1019 51646 : } else Signals[0] = 0.;
1020 51646 : if((signal = fDigits->GetData(Max.row, Max.col+1, Max.time))){
1021 51578 : gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(Max.col+1,Max.row);
1022 103156 : onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,Max.col+1) : 1;
1023 51578 : Signals[2] = (signal - fBaseline) /( onlcf * gain) + 0.5f;
1024 51646 : } else Signals[2] = 0.;
1025 :
1026 51646 : if(!(status[0] | status[1] | status[2])) {//all pads are good
1027 86924 : if ((Signals[2] <= Signals[1]) && (Signals[0] < Signals[1])) {
1028 25485 : if ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh)) {
1029 19386 : if(Signals[0]<0) Signals[0]=0.;
1030 19460 : if(Signals[2]<0) Signals[2]=0.;
1031 19078 : Float_t noiseSumThresh = fMinLeftRightCutSigma * fCalNoiseDetValue
1032 19078 : * fCalNoiseROC->GetValue(Max.col, Max.row);
1033 19098 : if ((Signals[2]+Signals[0]+Signals[1]) <= noiseSumThresh) return kFALSE;
1034 19058 : padStatus = 0;
1035 19058 : return kTRUE;
1036 : }
1037 : }
1038 : } else { // at least one of the pads is bad, and reject candidates with more than 1 problematic pad
1039 0 : if(Signals[0]<0)Signals[0]=0;
1040 0 : if(Signals[2]<0)Signals[2]=0;
1041 0 : if (status[2] && (!(status[0] || status[1])) && Signals[1] > Signals[0] && Signals[0] > fSigThresh) {
1042 0 : Signals[2]=0;
1043 0 : SetPadStatus(status[2], padStatus);
1044 0 : return kTRUE;
1045 : }
1046 0 : else if (status[0] && (!(status[1] || status[2])) && Signals[1] >= Signals[2] && Signals[2] > fSigThresh) {
1047 0 : Signals[0]=0;
1048 0 : SetPadStatus(status[0], padStatus);
1049 0 : return kTRUE;
1050 : }
1051 0 : else if (status[1] && (!(status[0] || status[2])) && ((Signals[2] > fSigThresh) || (Signals[0] > fSigThresh))) {
1052 0 : Signals[1] = fMaxThresh;
1053 0 : SetPadStatus(status[1], padStatus);
1054 0 : return kTRUE;
1055 : }
1056 : }
1057 32568 : return kFALSE;
1058 63043 : }
1059 :
1060 : //_____________________________________________________________________________
1061 : Bool_t AliTRDclusterizer::FivePadCluster(MaxStruct &ThisMax, MaxStruct &NeighbourMax)
1062 : {
1063 : //
1064 : // Look for 5 pad cluster with minimum in the middle
1065 : // Gives back the ratio
1066 : //
1067 :
1068 1118 : if (ThisMax.col >= fColMax - 3) return kFALSE;
1069 : Float_t gain;
1070 : Float_t onlcf;
1071 559 : if (ThisMax.col < fColMax - 5){
1072 525 : gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col+4,ThisMax.row);
1073 1050 : onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(ThisMax.row,ThisMax.col+4) : 1;
1074 525 : if (fDigits->GetData(ThisMax.row, ThisMax.col+4, ThisMax.time) - fBaseline >= fSigThresh * gain * onlcf)
1075 214 : return kFALSE;
1076 : }
1077 345 : if (ThisMax.col > 1) {
1078 336 : gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(ThisMax.col-2,ThisMax.row);
1079 672 : onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(ThisMax.row,ThisMax.col-2) : 1;
1080 336 : if (fDigits->GetData(ThisMax.row, ThisMax.col-2, ThisMax.time) - fBaseline >= fSigThresh * gain * onlcf)
1081 76 : return kFALSE;
1082 : }
1083 :
1084 : const Float_t kEpsilon = 0.01;
1085 807 : Double_t padSignal[5] = {ThisMax.signals[0], ThisMax.signals[1], ThisMax.signals[2],
1086 538 : NeighbourMax.signals[1], NeighbourMax.signals[2]};
1087 :
1088 : // Unfold the two maxima and set the signal on
1089 : // the overlapping pad to the ratio
1090 269 : Float_t ratio = Unfold(kEpsilon,fLayer,padSignal);
1091 269 : ThisMax.signals[2] = ThisMax.signals[2]*ratio + 0.5f;
1092 269 : NeighbourMax.signals[0] = NeighbourMax.signals[0]*(1-ratio) + 0.5f;
1093 269 : ThisMax.fivePad=kTRUE;
1094 269 : NeighbourMax.fivePad=kTRUE;
1095 : return kTRUE;
1096 :
1097 828 : }
1098 :
1099 : //_____________________________________________________________________________
1100 : void AliTRDclusterizer::CreateCluster(const MaxStruct &Max)
1101 : {
1102 : //
1103 : // Creates a cluster at the given position and saves it in RecPoints
1104 : //
1105 :
1106 38116 : Int_t nPadCount = 1;
1107 19058 : Short_t signals[7] = { 0, 0, (Short_t)Max.signals[0], (Short_t)Max.signals[1], (Short_t)Max.signals[2], 0, 0 };
1108 38116 : if(!fReconstructor->IsHLT()) CalcAdditionalInfo(Max, signals, nPadCount);
1109 :
1110 19058 : AliTRDcluster cluster(fDet, ((UChar_t) Max.col), ((UChar_t) Max.row), ((UChar_t) Max.time), signals, fVolid);
1111 19058 : cluster.SetNPads(nPadCount);
1112 19058 : cluster.SetQ(Max.signals[0]+Max.signals[1]+Max.signals[2]);
1113 38116 : if(TestBit(kLUT)) cluster.SetRPhiMethod(AliTRDcluster::kLUT);
1114 0 : else if(TestBit(kGAUS)) cluster.SetRPhiMethod(AliTRDcluster::kGAUS);
1115 0 : else cluster.SetRPhiMethod(AliTRDcluster::kCOG);
1116 :
1117 19058 : cluster.SetFivePad(Max.fivePad);
1118 : // set pads status for the cluster
1119 19058 : UChar_t maskPosition = GetCorruption(Max.padStatus);
1120 19058 : if (maskPosition) {
1121 0 : cluster.SetPadMaskedPosition(maskPosition);
1122 0 : cluster.SetPadMaskedStatus(GetPadStatus(Max.padStatus));
1123 : }
1124 19058 : cluster.SetXcorr(fReconstructor->UseClusterRadialCorrection());
1125 :
1126 : // Transform the local cluster coordinates into calibrated
1127 : // space point positions defined in the local tracking system.
1128 : // Here the calibration for T0, Vdrift and ExB is applied as well.
1129 57780 : if(!TestBit(kSkipTrafo)) if(!fTransform->Transform(&cluster)) return;
1130 : // Store raw signals in cluster. This MUST be called after position reconstruction !
1131 : // Xianguo Lu and Alex Bercuci 19.03.2012
1132 36904 : if(TestBit(kRawSignal) && fDigitsRaw){
1133 : Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
1134 18452 : Short_t rawSignal[7] = {0};
1135 295232 : for(Int_t ipad(Max.col-3), iRawId(0); ipad<=Max.col+3; ipad++, iRawId++){
1136 257804 : if(ipad<0 || ipad>=fColMax) continue;
1137 256560 : if(!fCalOnlGainROC){
1138 128280 : rawSignal[iRawId] = fDigitsRaw->GetData(Max.row, ipad, Max.time);
1139 128280 : continue;
1140 : }
1141 : // Deconvolute online gain calibration when available
1142 : // Alex Bercuci 27.04.2012
1143 128280 : tmp = (fDigitsRaw->GetData(Max.row, ipad, Max.time) - fBaseline)/fCalOnlGainROC->GetGainCorrectionFactor(Max.row, ipad) + 0.5f;
1144 0 : rawSignal[iRawId] = (Short_t)TMath::Min(tmp, kMaxShortVal);
1145 0 : }
1146 18452 : cluster.SetSignals(rawSignal, kTRUE);
1147 18452 : }
1148 : // Temporarily store the Max.Row, column and time bin of the center pad
1149 : // Used to later on assign the track indices
1150 18452 : cluster.SetLabel(Max.row, 0);
1151 18452 : cluster.SetLabel(Max.col, 1);
1152 18452 : cluster.SetLabel(Max.time,2);
1153 :
1154 : //needed for HLT reconstruction
1155 18452 : AddClusterToArray(&cluster);
1156 :
1157 : // Store the index of the first cluster in the current ROC
1158 18801 : if (firstClusterROC < 0) firstClusterROC = fNoOfClusters;
1159 :
1160 18452 : fNoOfClusters++;
1161 18452 : fClusterROC++;
1162 37510 : }
1163 :
1164 : //_____________________________________________________________________________
1165 : void AliTRDclusterizer::CalcAdditionalInfo(const MaxStruct &Max, Short_t *const signals, Int_t &nPadCount)
1166 : {
1167 : // Calculate number of pads/cluster and
1168 : // ADC signals at position 0, 1, 5 and 6
1169 :
1170 : Float_t tmp(0.), kMaxShortVal(32767.); // protect against data overflow due to wrong gain calibration
1171 : Float_t gain(1.); Float_t onlcf(1.); Short_t signal(0);
1172 : // Store the amplitudes of the pads in the cluster for later analysis
1173 : // and check whether one of these pads is masked in the database
1174 38116 : signals[3]=Max.signals[1];
1175 : Int_t ipad(1), jpad(0);
1176 : // Look to the right
1177 86760 : while((jpad = Max.col-ipad)){
1178 42980 : if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1179 42824 : gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1180 85648 : onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,jpad) : 1;
1181 42824 : tmp = (signal - fBaseline) / (onlcf * gain) + 0.5f;
1182 42824 : signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1183 42824 : if(signal<fSigThresh) break; // signal under threshold
1184 24322 : nPadCount++;
1185 44211 : if(ipad<=3) signals[3 - ipad] = signal;
1186 24322 : ipad++;
1187 : }
1188 : ipad=1;
1189 : // Look to the left
1190 87548 : while((jpad = Max.col+ipad)<fColMax){
1191 43524 : if(!(signal = fDigits->GetData(Max.row, jpad, Max.time))) break; // empty digit !
1192 43343 : gain = fCalGainFactorDetValue * fCalGainFactorROC->GetValue(jpad, Max.row);
1193 86686 : onlcf = fCalOnlGainROC ? fCalOnlGainROC->GetGainCorrectionFactor(Max.row,jpad) : 1;
1194 43343 : tmp = (signal - fBaseline) / (onlcf * gain) + 0.5f;
1195 43343 : signal = (Short_t)TMath::Min(tmp, kMaxShortVal);
1196 43343 : if(signal<fSigThresh) break; // signal under threshold
1197 24716 : nPadCount++;
1198 44888 : if(ipad<=3) signals[3 + ipad] = signal;
1199 24716 : ipad++;
1200 : }
1201 :
1202 57174 : AliDebug(4, Form("Signals[%3d %3d %3d %3d %3d %3d %3d] Npads[%d]."
1203 : , signals[0], signals[1], signals[2], signals[3], signals[4], signals[5], signals[6], nPadCount));
1204 19058 : }
1205 :
1206 : //_____________________________________________________________________________
1207 : void AliTRDclusterizer::AddClusterToArray(AliTRDcluster* cluster)
1208 : {
1209 : //
1210 : // Add a cluster to the array
1211 : //
1212 :
1213 36904 : Int_t n = RecPoints()->GetEntriesFast();
1214 18452 : if(n!=fNoOfClusters)AliError(Form("fNoOfClusters != RecPoints()->GetEntriesFast %i != %i \n", fNoOfClusters, n));
1215 18452 : new((*RecPoints())[n]) AliTRDcluster(*cluster);
1216 18452 : }
1217 :
1218 : //_____________________________________________________________________________
1219 : Bool_t AliTRDclusterizer::AddLabels()
1220 : {
1221 : //
1222 : // Add the track indices to the found clusters
1223 : //
1224 :
1225 : const Int_t kNclus = 3;
1226 : const Int_t kNdict = AliTRDdigitsManager::kNDict;
1227 : const Int_t kNtrack = kNdict * kNclus;
1228 :
1229 : Int_t iClusterROC = 0;
1230 :
1231 : Int_t row = 0;
1232 : Int_t col = 0;
1233 : Int_t time = 0;
1234 : Int_t iPad = 0;
1235 :
1236 : // Temporary array to collect the track indices
1237 376 : Int_t *idxTracks = new Int_t[kNtrack*fClusterROC];
1238 :
1239 : // Loop through the dictionary arrays one-by-one
1240 : // to keep memory consumption low
1241 : AliTRDarrayDictionary *tracksIn(NULL); //mod
1242 1504 : for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1243 :
1244 : // tracksIn should be expanded beforehand!
1245 564 : tracksIn = (AliTRDarrayDictionary *) fDigitsManager->GetDictionary(fDet,iDict);
1246 :
1247 : // Loop though the clusters found in this ROC
1248 58332 : for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1249 :
1250 28602 : AliTRDcluster *cluster = (AliTRDcluster *)
1251 28602 : RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1252 28602 : row = cluster->GetLabel(0);
1253 28602 : col = cluster->GetLabel(1);
1254 28602 : time = cluster->GetLabel(2);
1255 :
1256 228816 : for (iPad = 0; iPad < kNclus; iPad++) {
1257 85806 : Int_t iPadCol = col - 1 + iPad;
1258 85806 : Int_t index = tracksIn->GetData(row,iPadCol,time); //Modification of -1 in Track
1259 85806 : idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1260 : }
1261 :
1262 : }
1263 :
1264 : }
1265 :
1266 : // Copy the track indices into the cluster
1267 : // Loop though the clusters found in this ROC
1268 19444 : for (iClusterROC = 0; iClusterROC < fClusterROC; iClusterROC++) {
1269 :
1270 9534 : AliTRDcluster *cluster = (AliTRDcluster *)
1271 9534 : RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1272 9534 : cluster->SetLabel(-9999,0);
1273 9534 : cluster->SetLabel(-9999,1);
1274 9534 : cluster->SetLabel(-9999,2);
1275 :
1276 9534 : cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1277 :
1278 : }
1279 :
1280 376 : delete [] idxTracks;
1281 :
1282 188 : return kTRUE;
1283 :
1284 : }
1285 :
1286 : //_____________________________________________________________________________
1287 : Float_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, const Double_t *const padSignal) const
1288 : {
1289 : //
1290 : // Method to unfold neighbouring maxima.
1291 : // The charge ratio on the overlapping pad is calculated
1292 : // until there is no more change within the range given by eps.
1293 : // The resulting ratio is then returned to the calling method.
1294 : //
1295 :
1296 538 : AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1297 269 : if (!calibration) {
1298 0 : AliError("No AliTRDcalibDB instance available\n");
1299 0 : return kFALSE;
1300 : }
1301 :
1302 : Int_t irc = 0;
1303 : Int_t itStep = 0; // Count iteration steps
1304 :
1305 : Double_t ratio = 0.5; // Start value for ratio
1306 : Double_t prevRatio = 0.0; // Store previous ratio
1307 :
1308 269 : Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1309 269 : Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1310 269 : Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1311 :
1312 : // Start the iteration
1313 2692 : while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1314 :
1315 1077 : itStep++;
1316 : prevRatio = ratio;
1317 :
1318 : // Cluster position according to charge ratio
1319 1077 : Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1320 1077 : / (padSignal[0] + padSignal[1] + ratio * padSignal[2]);
1321 1077 : Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1322 1077 : / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1323 :
1324 : // Set cluster charge ratio
1325 1077 : irc = calibration->PadResponse(1.0, maxLeft, layer, newSignal);
1326 1077 : Double_t ampLeft = padSignal[1] / newSignal[1];
1327 1077 : irc = calibration->PadResponse(1.0, maxRight, layer, newSignal);
1328 1077 : Double_t ampRight = padSignal[3] / newSignal[1];
1329 :
1330 : // Apply pad response to parameters
1331 1077 : irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1332 1077 : irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1333 :
1334 : // Calculate new overlapping ratio
1335 : // Coverity
1336 1077 : if (irc != 0) {
1337 1077 : ratio = TMath::Min((Double_t) 1.0
1338 1077 : ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1339 1077 : }
1340 :
1341 : }
1342 :
1343 269 : return ratio;
1344 :
1345 538 : }
1346 :
1347 : //_____________________________________________________________________________
1348 : void AliTRDclusterizer::TailCancelation(const AliTRDrecoParam* const recoParam)
1349 : {
1350 : //
1351 : // Applies the tail cancelation
1352 : //
1353 :
1354 752 : Int_t nexp = recoParam->GetTCnexp();
1355 376 : if(!nexp) return;
1356 :
1357 376 : Int_t iRow = 0;
1358 376 : Int_t iCol = 0;
1359 : Int_t iTime = 0;
1360 :
1361 376 : TTreeSRedirector *fDebugStream = fReconstructor->GetDebugStream(AliTRDrecoParam::kClusterizer);
1362 752 : Bool_t debugStreaming = recoParam->GetStreamLevel(AliTRDrecoParam::kClusterizer) > 7 && fReconstructor->IsDebugStreaming();
1363 14084 : while(fIndexes->NextRCIndex(iRow, iCol))
1364 : {
1365 : // if corrupted then don't make the tail cancallation
1366 6666 : if (fCalPadStatusROC->GetStatus(iCol, iRow)) continue;
1367 :
1368 6666 : if(debugStreaming){
1369 0 : for (iTime = 0; iTime < fTimeTotal; iTime++)
1370 0 : (*fDebugStream) << "TailCancellation"
1371 0 : << "col=" << iCol
1372 0 : << "row=" << iRow
1373 0 : << "\n";
1374 : }
1375 :
1376 : // Apply the tail cancelation via the digital filter
1377 : //DeConvExp(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1378 6666 : ApplyTCTM(fDigits->GetDataAddress(iRow,iCol),fTimeTotal,nexp);
1379 : } // while irow icol
1380 :
1381 : return;
1382 :
1383 752 : }
1384 :
1385 :
1386 : //_____________________________________________________________________________
1387 : void AliTRDclusterizer::ApplyTCTM(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1388 : {
1389 : //
1390 : // Steer tail cancellation
1391 : //
1392 :
1393 :
1394 19998 : switch(nexp) {
1395 : case 1:
1396 : case 2:
1397 6666 : DeConvExp(arr,nTime,nexp);
1398 6666 : break;
1399 : case -1:
1400 0 : ConvExp(arr,nTime);
1401 0 : break;
1402 : case -2:
1403 0 : DeConvExp(arr,nTime,1);
1404 0 : ConvExp(arr,nTime);
1405 0 : break;
1406 : default:
1407 : break;
1408 : }
1409 6666 : }
1410 :
1411 :
1412 : //_____________________________________________________________________________
1413 : void AliTRDclusterizer::ConvExp(Short_t *const arr, const Int_t nTime)
1414 : {
1415 : //
1416 : // Tail maker
1417 : //
1418 :
1419 : // Initialization (coefficient = alpha, rates = lambda)
1420 : Float_t slope = 1.0;
1421 : Float_t coeff = 0.5;
1422 : Float_t rate;
1423 :
1424 0 : Double_t par[4];
1425 0 : fReconstructor->GetRecoParam()->GetTCParams(par);
1426 0 : slope = par[1];
1427 0 : coeff = par[3];
1428 :
1429 : Double_t dt = 0.1;
1430 :
1431 0 : rate = TMath::Exp(-dt/(slope));
1432 :
1433 : Float_t reminder = .0;
1434 : Float_t correction = 0.0;
1435 : Float_t result = 0.0;
1436 :
1437 0 : for (int i = nTime-1; i >= 0; i--) {
1438 :
1439 0 : result = arr[i] + correction - fBaseline; // No rescaling
1440 0 : arr[i] = (Short_t)(result + fBaseline + 0.5f);
1441 :
1442 : correction = 0.0;
1443 :
1444 0 : correction += reminder = rate * (reminder + coeff * result);
1445 : }
1446 0 : }
1447 :
1448 :
1449 : //_____________________________________________________________________________
1450 : void AliTRDclusterizer::DeConvExp(Short_t *const arr, const Int_t nTime, const Int_t nexp)
1451 : {
1452 : //
1453 : // Tail cancellation by deconvolution for PASA v4 TRF
1454 : //
1455 :
1456 13332 : Float_t rates[2];
1457 6666 : Float_t coefficients[2];
1458 :
1459 : // Initialization (coefficient = alpha, rates = lambda)
1460 : Float_t r1 = 1.0;
1461 : Float_t r2 = 1.0;
1462 : Float_t c1 = 0.5;
1463 : Float_t c2 = 0.5;
1464 :
1465 6666 : if (nexp == 1) { // 1 Exponentials
1466 : r1 = 1.156;
1467 : r2 = 0.130;
1468 : c1 = 0.066;
1469 : c2 = 0.000;
1470 6666 : }
1471 6666 : if (nexp == 2) { // 2 Exponentials
1472 0 : Double_t par[4];
1473 0 : fReconstructor->GetRecoParam()->GetTCParams(par);
1474 0 : r1 = par[0];//1.156;
1475 0 : r2 = par[1];//0.130;
1476 0 : c1 = par[2];//0.114;
1477 0 : c2 = par[3];//0.624;
1478 0 : }
1479 :
1480 6666 : coefficients[0] = c1;
1481 6666 : coefficients[1] = c2;
1482 :
1483 : Double_t dt = 0.1;
1484 :
1485 6666 : rates[0] = TMath::Exp(-dt/(r1));
1486 13332 : rates[1] = (nexp == 1) ? .0 : TMath::Exp(-dt/(r2));
1487 :
1488 6666 : Float_t reminder[2] = { .0, .0 };
1489 : Float_t correction = 0.0;
1490 : Float_t result = 0.0;
1491 :
1492 373296 : for (int i = 0; i < nTime; i++) {
1493 :
1494 179982 : result = arr[i] - correction - fBaseline; // No rescaling
1495 179982 : arr[i] = (Short_t)(result + fBaseline + 0.5f);
1496 :
1497 : correction = 0.0;
1498 1079892 : for (int k = 0; k < 2; k++) {
1499 359964 : correction += reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1500 : }
1501 : }
1502 6666 : }
1503 :
1504 : //_____________________________________________________________________________
1505 : TClonesArray *AliTRDclusterizer::RecPoints()
1506 : {
1507 : //
1508 : // Returns the list of rec points
1509 : //
1510 :
1511 187032 : fRecPoints = AliTRDReconstructor::GetClusters();
1512 93516 : if (!fRecPoints) AliError("Missing cluster array");
1513 93516 : return fRecPoints;
1514 : }
1515 :
1516 : //_____________________________________________________________________________
1517 : TClonesArray *AliTRDclusterizer::TrackletsArray(const TString &trkltype)
1518 : {
1519 : //
1520 : // Returns the array of on-line tracklets
1521 : //
1522 48 : fTracklets = AliTRDReconstructor::GetTracklets(trkltype.Data());
1523 24 : if (!fTracklets) AliError("Missing online tracklets array");
1524 24 : return fTracklets;
1525 : }
1526 :
1527 : //_____________________________________________________________________________
1528 : TClonesArray* AliTRDclusterizer::TracksArray()
1529 : {
1530 : // return array of GTU tracks (create TClonesArray if necessary)
1531 :
1532 48 : fTracks = AliTRDReconstructor::GetTracks();
1533 24 : if (!fTracks) AliError("Missing online tracks array");
1534 24 : return fTracks;
1535 : }
|
