2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
19 ///////////////////////////////////////////////////////////////////////////////
21 // TRD cluster finder //
23 ///////////////////////////////////////////////////////////////////////////////
30 #include "AliRunLoader.h"
31 #include "AliLoader.h"
32 #include "AliRawReader.h"
34 #include "AliAlignObj.h"
36 #include "AliTRDclusterizerV2.h"
37 #include "AliTRDgeometry.h"
38 #include "AliTRDdataArrayF.h"
39 #include "AliTRDdataArrayI.h"
40 #include "AliTRDdigitsManager.h"
41 #include "AliTRDpadPlane.h"
42 #include "AliTRDrawData.h"
43 #include "AliTRDcalibDB.h"
44 #include "AliTRDSimParam.h"
45 #include "AliTRDRecParam.h"
46 #include "AliTRDcluster.h"
48 #include "Cal/AliTRDCalROC.h"
49 #include "Cal/AliTRDCalDet.h"
51 #include "AliTRDSignalIndex.h"
52 #include "AliTRDRawStream.h"
53 #include "AliTRDRawStreamV2.h"
55 #include "AliTRDfeeParam.h"
57 ClassImp(AliTRDclusterizerV2)
59 //_____________________________________________________________________________
60 AliTRDclusterizerV2::AliTRDclusterizerV2()
70 // AliTRDclusterizerV2 default constructor
73 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
76 //_____________________________________________________________________________
77 AliTRDclusterizerV2::AliTRDclusterizerV2(const Text_t *name, const Text_t *title)
78 :AliTRDclusterizer(name,title)
79 ,fDigitsManager(new AliTRDdigitsManager())
87 // AliTRDclusterizerV2 constructor
90 fDigitsManager->CreateArrays();
91 fGeometry = new AliTRDgeometry;
93 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
96 //_____________________________________________________________________________
97 AliTRDclusterizerV2::AliTRDclusterizerV2(const AliTRDclusterizerV2 &c)
104 ,fIndexesMaxima(NULL)
107 // AliTRDclusterizerV2 copy constructor
112 //_____________________________________________________________________________
113 AliTRDclusterizerV2::~AliTRDclusterizerV2()
116 // AliTRDclusterizerV2 destructor
119 if (fDigitsManager) {
120 delete fDigitsManager;
121 fDigitsManager = NULL;
138 delete fIndexesMaxima;
139 fIndexesMaxima = NULL;
143 //_____________________________________________________________________________
144 AliTRDclusterizerV2 &AliTRDclusterizerV2::operator=(const AliTRDclusterizerV2 &c)
147 // Assignment operator
150 if (this != &c) ((AliTRDclusterizerV2 &) c).Copy(*this);
155 //_____________________________________________________________________________
156 void AliTRDclusterizerV2::Copy(TObject &c) const
162 ((AliTRDclusterizerV2 &) c).fDigitsManager = NULL;
163 ((AliTRDclusterizerV2 &) c).fGeometry = NULL;
164 ((AliTRDclusterizerV2 &) c).fAddLabels = fAddLabels;
165 ((AliTRDclusterizerV2 &) c).fRawVersion = fRawVersion;
166 ((AliTRDclusterizerV2 &) c).fIndexesOut = NULL;
167 ((AliTRDclusterizerV2 &) c).fIndexesMaxima = NULL;
169 AliTRDclusterizer::Copy(c);
173 //_____________________________________________________________________________
174 void AliTRDclusterizerV2::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
177 // Reset the helper indexes
181 // carefull here - we assume that only row number may change - most probable
182 if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
183 fIndexesOut->ResetContent();
185 fIndexesOut->ResetContentConditional(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
189 fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
194 // carefull here - we assume that only row number may change - most probable
195 if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
196 fIndexesMaxima->ResetContent();
198 fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
202 fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow(), indexesIn->GetNcol(), indexesIn->GetNtime());
206 //_____________________________________________________________________________
207 Bool_t AliTRDclusterizerV2::ReadDigits()
210 // Reads the digits arrays from the input aliroot file
214 AliError("No run loader available");
218 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
219 if (!loader->TreeD()) {
220 loader->LoadDigits();
223 // Read in the digit arrays
224 return (fDigitsManager->ReadDigits(loader->TreeD()));
228 //_____________________________________________________________________________
229 Bool_t AliTRDclusterizerV2::ReadDigits(TTree *digitsTree)
232 // Reads the digits arrays from the input tree
235 // Read in the digit arrays
236 return (fDigitsManager->ReadDigits(digitsTree));
240 //_____________________________________________________________________________
241 Bool_t AliTRDclusterizerV2::ReadDigits(AliRawReader *rawReader)
244 // Reads the digits arrays from the ddl file
248 fDigitsManager = raw.Raw2Digits(rawReader);
254 //_____________________________________________________________________________
255 Bool_t AliTRDclusterizerV2::MakeClusters()
258 // Creates clusters from digits
261 //propagate info from the digits manager
262 if (fAddLabels == kTRUE)
263 fAddLabels = fDigitsManager->UsesDictionaries();
265 Bool_t fReturn = kTRUE;
266 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
268 AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(i);
269 // This is to take care of switched off super modules
270 if (digitsIn->GetNtime() == 0) {
273 //AliInfo(Form("digitsIn->Expand() 0x%x", digitsIn));
275 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
276 if (indexes->IsAllocated() == kFALSE)
278 fDigitsManager->BuildIndexes(i);
282 if (indexes->HasEntry())
286 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
288 AliTRDdataArrayI *tracksIn = 0;
289 tracksIn = fDigitsManager->GetDictionary(i,iDict);
293 fR = MakeClusters(i);
294 fReturn = fR && fReturn;
302 //digitsIn->Compress(1,0);
303 // no compress just remove
304 fDigitsManager->RemoveDigits(i);
305 fDigitsManager->RemoveDictionaries(i);
306 fDigitsManager->ClearIndexes(i);
313 //_____________________________________________________________________________
314 Bool_t AliTRDclusterizerV2::Raw2Clusters(AliRawReader *rawReader)
317 // Creates clusters from raw data
320 //AliDebug(2, "Read all");
322 AliTRDdataArrayI *digits = 0;
323 AliTRDdataArrayI *track0 = 0;
324 AliTRDdataArrayI *track1 = 0;
325 AliTRDdataArrayI *track2 = 0;
327 AliTRDSignalIndex *indexes = 0;
328 // Create the digits manager
331 fDigitsManager = new AliTRDdigitsManager();
332 fDigitsManager->CreateArrays();
335 //AliTRDRawStream input(rawReader);
336 AliTRDRawStreamV2 input(rawReader);
337 input.SetRawVersion( fRawVersion );
340 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
342 // Loop through the digits
349 det = input.GetDet();
355 digits = fDigitsManager->GetDigits(lastdet);
356 Bool_t iclusterBranch = kFALSE;
357 if (indexes->HasEntry())
358 iclusterBranch = MakeClusters(lastdet);
359 if (iclusterBranch == kFALSE)
361 WriteClusters(lastdet);
368 fDigitsManager->RemoveDigits(lastdet);
369 fDigitsManager->RemoveDictionaries(lastdet);
370 fDigitsManager->ClearIndexes(lastdet);
375 // Add a container for the digits of this detector
376 digits = fDigitsManager->GetDigits(det);
377 track0 = fDigitsManager->GetDictionary(det,0);
378 track1 = fDigitsManager->GetDictionary(det,1);
379 track2 = fDigitsManager->GetDictionary(det,2);
381 // Allocate memory space for the digits buffer
382 if (digits->GetNtime() == 0)
384 //AliDebug(5, Form("Alloc digits for det %d", det));
385 digits->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
386 track0->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
387 track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
388 track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins());
391 indexes = fDigitsManager->GetIndexes(det);
392 indexes->SetSM(input.GetSM());
393 indexes->SetStack(input.GetStack());
394 indexes->SetLayer(input.GetLayer());
395 indexes->SetDetNumber(det);
396 if (indexes->IsAllocated() == kFALSE)
397 indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins());
400 for (it = 0; it < 3; it++)
402 if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() )
404 if (input.GetSignals()[it] > 0)
406 digits->SetDataUnchecked(input.GetRow(), input.GetCol(),
407 input.GetTimeBin() + it, input.GetSignals()[it]);
409 indexes->AddIndexTBin(input.GetRow(), input.GetCol(),
410 input.GetTimeBin() + it);
411 track0->SetDataUnchecked(input.GetRow(), input.GetCol(),
412 input.GetTimeBin() + it, 0);
413 track1->SetDataUnchecked(input.GetRow(), input.GetCol(),
414 input.GetTimeBin() + it, 0);
415 track2->SetDataUnchecked(input.GetRow(), input.GetCol(),
416 input.GetTimeBin() + it, 0);
424 Bool_t iclusterBranch = kFALSE;
425 if (indexes->HasEntry())
426 iclusterBranch = MakeClusters(lastdet);
427 if (iclusterBranch == kFALSE)
429 WriteClusters(lastdet);
432 //MakeClusters(lastdet);
435 fDigitsManager->RemoveDigits(lastdet);
436 fDigitsManager->RemoveDictionaries(lastdet);
437 fDigitsManager->ClearIndexes(lastdet);
441 delete fDigitsManager;
442 fDigitsManager = NULL;
446 //_____________________________________________________________________________
447 Bool_t AliTRDclusterizerV2::Raw2ClustersChamber(AliRawReader *rawReader)
450 // Creates clusters from raw data
453 //AliDebug(1, "Raw2ClustersChamber");
455 // Create the digits manager
456 // Create the digits manager
459 fDigitsManager = new AliTRDdigitsManager();
460 fDigitsManager->CreateArrays();
463 fDigitsManager->SetUseDictionaries(fAddLabels);
465 //AliTRDRawStream input(rawReader);
466 AliTRDRawStreamV2 input(rawReader);
467 input.SetRawVersion( fRawVersion );
470 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
473 while ((det = input.NextChamber(fDigitsManager)) >= 0)
475 Bool_t iclusterBranch = kFALSE;
476 if (fDigitsManager->GetIndexes(det)->HasEntry())
477 iclusterBranch = MakeClusters(det);
478 if (iclusterBranch == kFALSE)
483 fDigitsManager->RemoveDigits(det);
484 fDigitsManager->RemoveDictionaries(det);
485 fDigitsManager->ClearIndexes(det);
488 delete fDigitsManager;
489 fDigitsManager = NULL;
493 //_____________________________________________________________________________
494 Bool_t AliTRDclusterizerV2::MakeClusters(Int_t det)
497 // Generates the cluster.
501 // digits should be expanded beforehand!
502 // digitsIn->Expand();
503 AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(det);
504 // This is to take care of switched off super modules
505 if (digitsIn->GetNtime() == 0)
507 //AliDebug(5, Form("digitsIn->GetNtime() == 0 [%d]", det));
511 AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
512 if (indexesIn->IsAllocated() == kFALSE)
514 AliError("Indexes do not exist!");
518 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
520 AliFatal("No AliTRDcalibDB instance available\n");
524 AliTRDSimParam *simParam = AliTRDSimParam::Instance();
526 AliError("No AliTRDSimParam instance available\n");
530 AliTRDRecParam *recParam = AliTRDRecParam::Instance();
532 AliError("No AliTRDRecParam instance available\n");
537 // Float_t ADCthreshold = simParam->GetADCthreshold();
538 Float_t ADCthreshold = 0; // There is no ADC threshold anymore, and simParam should not be used ni clusterizer. KO
540 // Threshold value for the maximum
541 Float_t maxThresh = recParam->GetClusMaxThresh();
542 // Threshold value for the digit signal
543 Float_t sigThresh = recParam->GetClusSigThresh();
545 // Detector wise calibration object for t0
546 const AliTRDCalDet *calT0Det = calibration->GetT0Det();
547 // Detector wise calibration object for the gain factors
548 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
550 // Iteration limit for unfolding procedure
551 const Float_t kEpsilon = 0.01;
552 const Int_t kNclus = 3;
553 const Int_t kNsig = 5;
556 Double_t ratioLeft = 1.0;
557 Double_t ratioRight = 1.0;
559 Double_t padSignal[kNsig];
560 Double_t clusterSignal[kNclus];
561 Double_t clusterPads[kNclus];
563 Int_t icham = indexesIn->GetChamber();
564 Int_t iplan = indexesIn->GetPlane();
565 Int_t isect = indexesIn->GetSM();
567 // Start clustering in the chamber
569 Int_t idet = fGeometry->GetDetector(iplan,icham,isect);
572 AliError("Strange Detector number Missmatch!");
576 Int_t ilayer = AliGeomManager::kTRD1 + iplan;
577 Int_t imodule = icham + AliTRDgeometry::Ncham() * isect;
578 UShort_t volid = AliGeomManager::LayerToVolUID(ilayer,imodule);
580 //Int_t nRowMax = digitsIn->GetNrow();
581 Int_t nColMax = digitsIn->GetNcol();
582 Int_t nTimeTotal = digitsIn->GetNtime();
584 AliTRDpadPlane *padPlane = fGeometry->GetPadPlane(iplan,icham);
586 // Calibration object with pad wise values for t0
587 AliTRDCalROC *calT0ROC = calibration->GetT0ROC(idet);
588 // Calibration object with pad wise values for the gain factors
589 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
590 // Calibration value for chamber wise t0
591 Float_t calT0DetValue = calT0Det->GetValue(idet);
592 // Calibration value for chamber wise gain factor
593 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
597 // Apply the gain and the tail cancelation via digital filter
598 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow()
600 ,digitsIn->GetNtime());
602 // AliInfo(Form("nrows %d cols %d time %d",
603 // digitsIn->GetNrow()
604 // ,digitsIn->GetNcol()
605 // ,digitsIn->GetNtime()));
607 ResetHelperIndexes(indexesIn);
616 ,calGainFactorDetValue);
623 fIndexesOut->ResetCounters();
624 while (fIndexesOut->NextRCTbinIndex(row, col, time))
626 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
628 // Look for the maximum
629 if (signalM >= maxThresh)
632 if (col + 1 >= nColMax || col-1 < 0)
635 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1 ,time));
637 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
639 if ((TMath::Abs(signalL) <= signalM) &&
640 (TMath::Abs(signalR) < signalM))
642 if ((TMath::Abs(signalL) >= sigThresh) ||
643 (TMath::Abs(signalR) >= sigThresh))
645 // Maximum found, mark the position by a negative signal
646 digitsOut->SetDataUnchecked(row,col,time,-signalM);
647 fIndexesMaxima->AddIndexTBin(row,col,time);
653 // The index to the first cluster of a given ROC
654 Int_t firstClusterROC = -1;
655 // The number of cluster in a given ROC
656 Int_t nClusterROC = 0;
658 // Now check the maxima and calculate the cluster position
659 fIndexesMaxima->ResetCounters();
660 while (fIndexesMaxima->NextRCTbinIndex(row, col, time))
663 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
665 for (iPad = 0; iPad < kNclus; iPad++) {
666 Int_t iPadCol = col - 1 + iPad;
667 clusterSignal[iPad] =
668 TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
671 // Count the number of pads in the cluster
676 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) {
679 if (col-ii < 0) break;
683 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) {
686 if (col+ii+1 >= nColMax) break;
690 // Look for 5 pad cluster with minimum in the middle
691 Bool_t fivePadCluster = kFALSE;
692 if (col < (nColMax - 3)) {
693 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
694 fivePadCluster = kTRUE;
696 if ((fivePadCluster) && (col < (nColMax - 5))) {
697 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) {
698 fivePadCluster = kFALSE;
701 if ((fivePadCluster) && (col > 1)) {
702 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) {
703 fivePadCluster = kFALSE;
709 // Modify the signal of the overlapping pad for the left part
710 // of the cluster which remains from a previous unfolding
712 clusterSignal[0] *= ratioLeft;
716 // Unfold the 5 pad cluster
717 if (fivePadCluster) {
718 for (iPad = 0; iPad < kNsig; iPad++) {
719 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
723 // Unfold the two maxima and set the signal on
724 // the overlapping pad to the ratio
725 ratioRight = Unfold(kEpsilon,iplan,padSignal);
726 ratioLeft = 1.0 - ratioRight;
727 clusterSignal[2] *= ratioRight;
731 Double_t clusterCharge = clusterSignal[0]
735 // The position of the cluster
736 clusterPads[0] = row + 0.5;
737 // Take the shift of the additional time bins into account
738 clusterPads[2] = time + 0.5;
740 if (recParam->LUTOn()) {
741 // Calculate the position of the cluster by using the
742 // lookup table method
743 clusterPads[1] = recParam->LUTposition(iplan,clusterSignal[0]
748 // Calculate the position of the cluster by using the
749 // center of gravity method
750 for (Int_t i = 0; i < kNsig; i++) {
753 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
754 padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left pad
755 padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right pad
757 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
758 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
760 if ((col < nColMax - 3) &&
761 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])) {
762 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
764 clusterPads[1] = GetCOG(padSignal);
767 Double_t q0 = clusterSignal[0];
768 Double_t q1 = clusterSignal[1];
769 Double_t q2 = clusterSignal[2];
770 Double_t clusterSigmaY2 = (q1 * (q0 + q2) + 4.0 * q0 * q2)
771 / (clusterCharge*clusterCharge);
774 // Calculate the position and the error
777 // Correct for t0 (sum of chamber and pad wise values !!!)
778 Float_t calT0ROCValue = calT0ROC->GetValue(col,row);
779 Char_t clusterTimeBin = ((Char_t) TMath::Nint(time - (calT0DetValue + calT0ROCValue)));
780 Double_t colSize = padPlane->GetColSize(col);
781 Double_t rowSize = padPlane->GetRowSize(row);
783 Float_t clusterPos[3];
784 clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1] + 0.5) * colSize;
785 clusterPos[1] = padPlane->GetRowPos(row) - 0.5 * rowSize;
786 clusterPos[2] = CalcXposFromTimebin(clusterPads[2],idet,col,row);
787 Float_t clusterSig[2];
788 clusterSig[0] = (clusterSigmaY2 + 1.0/12.0) * colSize*colSize;
789 clusterSig[1] = rowSize * rowSize / 12.0;
791 // Store the amplitudes of the pads in the cluster for later analysis
792 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
793 for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
795 (jPad >= nColMax-1)) {
798 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
801 // Add the cluster to the output array
802 // The track indices will be stored later
803 AliTRDcluster *cluster = new AliTRDcluster(idet
808 ,((Char_t) nPadCount)
814 // Temporarily store the row, column and time bin of the center pad
815 // Used to later on assign the track indices
816 cluster->SetLabel( row,0);
817 cluster->SetLabel( col,1);
818 cluster->SetLabel(time,2);
819 RecPoints()->Add(cluster);
821 // Store the index of the first cluster in the current ROC
822 if (firstClusterROC < 0) {
823 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
825 // Count the number of cluster in the current ROC
828 } // if: Maximum found ?
832 //delete fIndexesOut;
833 //delete fIndexesMaxima;
836 AddLabels(idet, firstClusterROC, nClusterROC);
838 // Write the cluster and reset the array
846 //_____________________________________________________________________________
847 Bool_t AliTRDclusterizerV2::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
850 // Add the track indices to the found clusters
853 const Int_t kNclus = 3;
854 const Int_t kNdict = AliTRDdigitsManager::kNDict;
855 const Int_t kNtrack = kNdict * kNclus;
857 Int_t iClusterROC = 0;
864 // Temporary array to collect the track indices
865 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
867 // Loop through the dictionary arrays one-by-one
868 // to keep memory consumption low
869 AliTRDdataArrayI *tracksIn = 0;
870 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
872 tracksIn = fDigitsManager->GetDictionary(idet,iDict);
873 // tracksIn should be expanded beforehand!
875 // Loop though the clusters found in this ROC
876 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
878 AliTRDcluster *cluster = (AliTRDcluster *)
879 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
880 row = cluster->GetLabel(0);
881 col = cluster->GetLabel(1);
882 time = cluster->GetLabel(2);
884 for (iPad = 0; iPad < kNclus; iPad++) {
885 Int_t iPadCol = col - 1 + iPad;
886 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
887 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
892 // Compress the arrays
893 // no do not compress - we will delete them when we are done with the detector
894 //tracksIn->Compress(1,0);
898 // Copy the track indices into the cluster
899 // Loop though the clusters found in this ROC
900 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
902 AliTRDcluster *cluster = (AliTRDcluster *)
903 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
904 cluster->SetLabel(-9999,0);
905 cluster->SetLabel(-9999,1);
906 cluster->SetLabel(-9999,2);
908 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
917 //_____________________________________________________________________________
918 Double_t AliTRDclusterizerV2::GetCOG(Double_t signal[5])
922 // Used for clusters with more than 3 pads - where LUT not applicable
925 Double_t sum = signal[0]
931 Double_t res = (0.0 * (-signal[0] + signal[4])
932 + (-signal[1] + signal[3])) / sum;
938 //_____________________________________________________________________________
939 Double_t AliTRDclusterizerV2::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
942 // Method to unfold neighbouring maxima.
943 // The charge ratio on the overlapping pad is calculated
944 // until there is no more change within the range given by eps.
945 // The resulting ratio is then returned to the calling method.
948 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
950 AliError("No AliTRDcalibDB instance available\n");
955 Int_t itStep = 0; // Count iteration steps
957 Double_t ratio = 0.5; // Start value for ratio
958 Double_t prevRatio = 0.0; // Store previous ratio
960 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
961 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
962 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
964 // Start the iteration
965 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
970 // Cluster position according to charge ratio
971 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
972 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
973 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
974 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
976 // Set cluster charge ratio
977 irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
978 Double_t ampLeft = padSignal[1] / newSignal[1];
979 irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
980 Double_t ampRight = padSignal[3] / newSignal[1];
982 // Apply pad response to parameters
983 irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
984 irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
986 // Calculate new overlapping ratio
987 ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] /
988 (newLeftSignal[2] + newRightSignal[0]));
996 //_____________________________________________________________________________
997 void AliTRDclusterizerV2::Transform(AliTRDdataArrayI *digitsIn
998 , AliTRDdataArrayF *digitsOut
999 , AliTRDSignalIndex *indexesIn
1000 , AliTRDSignalIndex *indexesOut
1002 , Float_t ADCthreshold
1003 , AliTRDCalROC *calGainFactorROC
1004 , Float_t calGainFactorDetValue)
1007 // Apply gain factor
1008 // Apply tail cancelation: Transform digitsIn to digitsOut
1015 AliTRDRecParam *recParam = AliTRDRecParam::Instance();
1017 AliError("No AliTRDRecParam instance available\n");
1021 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
1022 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
1023 indexesIn->ResetCounters();
1024 while (indexesIn->NextRCIndex(iRow, iCol))
1026 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1027 Double_t gain = calGainFactorDetValue
1028 * calGainFactorROCValue;
1030 for (iTime = 0; iTime < nTimeTotal; iTime++)
1035 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1036 inADC[iTime] /= gain;
1037 outADC[iTime] = inADC[iTime];
1040 // Apply the tail cancelation via the digital filter
1041 if (recParam->TCOn()) {
1042 DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp());
1045 indexesIn->ResetTbinCounter();
1046 while (indexesIn->NextTbinIndex(iTime))
1048 // Store the amplitude of the digit if above threshold
1049 if (outADC[iTime] > ADCthreshold)
1051 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1052 //AliDebug(5, Form("add index %d", indexesIn->GetDetNumber()));
1053 indexesOut->AddIndexTBin(iRow,iCol,iTime);
1061 //AliDebug(5, Form("Stop %d", indexesIn->GetDetNumber()));
1067 //_____________________________________________________________________________
1068 void AliTRDclusterizerV2::DeConvExp(Double_t *source, Double_t *target
1069 , Int_t n, Int_t nexp)
1072 // Tail cancellation by deconvolution for PASA v4 TRF
1076 Double_t coefficients[2];
1078 // Initialization (coefficient = alpha, rates = lambda)
1084 if (nexp == 1) { // 1 Exponentials
1090 if (nexp == 2) { // 2 Exponentials
1097 coefficients[0] = C1;
1098 coefficients[1] = C2;
1102 rates[0] = TMath::Exp(-Dt/(R1));
1103 rates[1] = TMath::Exp(-Dt/(R2));
1108 Double_t reminder[2];
1109 Double_t correction;
1112 // Attention: computation order is important
1114 for (k = 0; k < nexp; k++) {
1117 for (i = 0; i < n; i++) {
1118 result = (source[i] - correction); // No rescaling
1121 for (k = 0; k < nexp; k++) {
1122 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1125 for (k = 0; k < nexp; k++) {
1126 correction += reminder[k];