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 ///////////////////////////////////////////////////////////////////////////////
29 #include <TObjArray.h>
31 #include "AliRunLoader.h"
32 #include "AliLoader.h"
33 #include "AliRawReader.h"
35 #include "AliAlignObj.h"
37 #include "AliTRDclusterizer.h"
38 #include "AliTRDcluster.h"
39 #include "AliTRDReconstructor.h"
40 #include "AliTRDgeometry.h"
41 #include "AliTRDdataArrayF.h"
42 #include "AliTRDdataArrayI.h"
43 #include "AliTRDdataArrayS.h"
44 #include "AliTRDdataArrayDigits.h"
45 #include "AliTRDdigitsManager.h"
46 #include "AliTRDrawData.h"
47 #include "AliTRDcalibDB.h"
48 #include "AliTRDrecoParam.h"
49 #include "AliTRDCommonParam.h"
50 #include "AliTRDtransform.h"
51 #include "AliTRDSignalIndex.h"
52 #include "AliTRDrawStreamBase.h"
53 #include "AliTRDfeeParam.h"
55 #include "Cal/AliTRDCalROC.h"
56 #include "Cal/AliTRDCalDet.h"
57 #include "Cal/AliTRDCalSingleChamberStatus.h"
59 ClassImp(AliTRDclusterizer)
61 //_____________________________________________________________________________
62 AliTRDclusterizer::AliTRDclusterizer()
72 ,fTransform(new AliTRDtransform(0))
77 // AliTRDclusterizer default constructor
80 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
84 //_____________________________________________________________________________
85 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title)
90 ,fDigitsManager(new AliTRDdigitsManager())
95 ,fTransform(new AliTRDtransform(0))
100 // AliTRDclusterizer constructor
103 fDigitsManager->CreateArrays();
105 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
111 //_____________________________________________________________________________
112 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
117 ,fDigitsManager(NULL)
121 ,fIndexesMaxima(NULL)
127 // AliTRDclusterizer copy constructor
134 //_____________________________________________________________________________
135 AliTRDclusterizer::~AliTRDclusterizer()
138 // AliTRDclusterizer destructor
141 printf("----> AliTRDclusterizer::dtor\n");
145 fRecPoints->Delete();
151 delete fDigitsManager;
152 fDigitsManager = NULL;
163 delete fIndexesMaxima;
164 fIndexesMaxima = NULL;
181 //_____________________________________________________________________________
182 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
185 // Assignment operator
190 ((AliTRDclusterizer &) c).Copy(*this);
196 //_____________________________________________________________________________
197 void AliTRDclusterizer::Copy(TObject &c) const
203 ((AliTRDclusterizer &) c).fClusterTree = NULL;
204 ((AliTRDclusterizer &) c).fRecPoints = NULL;
205 ((AliTRDclusterizer &) c).fDigitsManager = NULL;
206 ((AliTRDclusterizer &) c).fAddLabels = fAddLabels;
207 ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
208 ((AliTRDclusterizer &) c).fIndexesOut = NULL;
209 ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
210 ((AliTRDclusterizer &) c).fTransform = NULL;
211 ((AliTRDclusterizer &) c).fLUTbin = 0;
212 ((AliTRDclusterizer &) c).fLUT = NULL;
216 //_____________________________________________________________________________
217 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
220 // Opens the AliROOT file. Output and input are in the same file
223 TString evfoldname = AliConfig::GetDefaultEventFolderName();
224 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
227 fRunLoader = AliRunLoader::Open(name);
231 AliError(Form("Can not open session for file %s.",name));
242 //_____________________________________________________________________________
243 Bool_t AliTRDclusterizer::OpenOutput()
246 // Open the output file
249 TObjArray *ioArray = 0;
251 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
252 loader->MakeTree("R");
254 fClusterTree = loader->TreeR();
255 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
261 //_____________________________________________________________________________
262 Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
265 // Connect the output tree
268 TObjArray *ioArray = 0;
270 fClusterTree = clusterTree;
271 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
277 //_____________________________________________________________________________
278 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
281 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
284 // Import the Trees for the event nEvent in the file
285 fRunLoader->GetEvent(nEvent);
291 //_____________________________________________________________________________
292 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
295 // Fills TRDcluster branch in the tree with the clusters
296 // found in detector = det. For det=-1 writes the tree.
300 (det >= AliTRDgeometry::Ndet())) {
301 AliError(Form("Unexpected detector index %d.\n",det));
305 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
307 TObjArray *ioArray = 0;
308 branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
312 (det < AliTRDgeometry::Ndet())) {
314 Int_t nRecPoints = RecPoints()->GetEntriesFast();
315 TObjArray *detRecPoints = new TObjArray(400);
317 for (Int_t i = 0; i < nRecPoints; i++) {
318 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
319 if (det == c->GetDetector()) {
320 detRecPoints->AddLast(c);
323 AliError(Form("Attempt to write a cluster with unexpected detector index: got=%d expected=%d\n"
329 branch->SetAddress(&detRecPoints);
330 fClusterTree->Fill();
340 AliInfo(Form("Writing the cluster tree %s for event %d."
341 ,fClusterTree->GetName(),fRunLoader->GetEventNumber()));
345 branch->SetAddress(&fRecPoints);
347 AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
348 loader->WriteRecPoints("OVERWRITE");
353 AliError("Cluster tree does not exist. Cannot write clusters.\n");
362 AliError(Form("Unexpected detector index %d.\n",det));
368 //_____________________________________________________________________________
369 void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
372 // Reset the helper indexes
377 // carefull here - we assume that only row number may change - most probable
378 if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
379 fIndexesOut->ResetContent();
381 fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
382 , indexesIn->GetNcol()
383 , indexesIn->GetNtime());
387 fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
388 , indexesIn->GetNcol()
389 , indexesIn->GetNtime());
394 // carefull here - we assume that only row number may change - most probable
395 if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
397 fIndexesMaxima->ResetContent();
401 fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
402 , indexesIn->GetNcol()
403 , indexesIn->GetNtime());
408 fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
409 , indexesIn->GetNcol()
410 , indexesIn->GetNtime());
415 //_____________________________________________________________________________
416 Bool_t AliTRDclusterizer::ReadDigits()
419 // Reads the digits arrays from the input aliroot file
423 AliError("No run loader available");
427 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
428 if (!loader->TreeD()) {
429 loader->LoadDigits();
432 // Read in the digit arrays
433 return (fDigitsManager->ReadDigits(loader->TreeD()));
437 //_____________________________________________________________________________
438 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
441 // Reads the digits arrays from the input tree
444 // Read in the digit arrays
445 return (fDigitsManager->ReadDigits(digitsTree));
449 //_____________________________________________________________________________
450 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
453 // Reads the digits arrays from the ddl file
457 fDigitsManager = raw.Raw2Digits(rawReader);
463 //_____________________________________________________________________________
464 Bool_t AliTRDclusterizer::MakeClusters()
467 // Creates clusters from digits
470 // Propagate info from the digits manager
471 if (fAddLabels == kTRUE)
473 fAddLabels = fDigitsManager->UsesDictionaries();
476 Bool_t fReturn = kTRUE;
477 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
480 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits*) fDigitsManager->GetDigits(i);
481 // This is to take care of switched off super modules
482 if (!digitsIn->HasData())
487 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
488 if (indexes->IsAllocated() == kFALSE)
490 fDigitsManager->BuildIndexes(i);
494 if (indexes->HasEntry())
498 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
500 AliTRDdataArrayI *tracksIn = 0;
501 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(i,iDict);
505 fR = MakeClusters(i);
506 fReturn = fR && fReturn;
515 // No compress just remove
516 fDigitsManager->RemoveDigits(i);
517 fDigitsManager->RemoveDictionaries(i);
518 fDigitsManager->ClearIndexes(i);
526 //_____________________________________________________________________________
527 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
530 // Creates clusters from raw data
533 return Raw2ClustersChamber(rawReader);
537 //_____________________________________________________________________________
538 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
541 // Creates clusters from raw data
544 // Create the digits manager
547 fDigitsManager = new AliTRDdigitsManager();
548 fDigitsManager->CreateArrays();
551 fDigitsManager->SetUseDictionaries(fAddLabels);
553 AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
554 AliTRDrawStreamBase &input = *pinput;
556 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
559 while ((det = input.NextChamber(fDigitsManager)) >= 0)
561 Bool_t iclusterBranch = kFALSE;
562 if (fDigitsManager->GetIndexes(det)->HasEntry())
564 iclusterBranch = MakeClusters(det);
566 if (iclusterBranch == kFALSE)
571 fDigitsManager->RemoveDigits(det);
572 fDigitsManager->RemoveDictionaries(det);
573 fDigitsManager->ClearIndexes(det);
576 delete fDigitsManager;
577 fDigitsManager = NULL;
585 //_____________________________________________________________________________
586 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
589 // check if pad is masked
590 if(signal>0 && TESTBIT(signal, 10)){
592 for(int ibit=0; ibit<4; ibit++){
593 if(TESTBIT(signal, 11+ibit)){
594 SETBIT(status, ibit);
595 CLRBIT(signal, 11+ibit);
602 //_____________________________________________________________________________
603 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
606 // Generates the cluster.
610 // digits should be expanded beforehand!
611 // digitsIn->Expand();
612 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits *) fDigitsManager->GetDigits(det);
614 // This is to take care of switched off super modules
615 if (!digitsIn->HasData())
620 AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
621 if (indexesIn->IsAllocated() == kFALSE)
623 AliError("Indexes do not exist!");
627 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
630 AliFatal("No AliTRDcalibDB instance available\n");
635 // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
636 Float_t adcThreshold = 0;
638 if (!AliTRDReconstructor::RecoParam())
640 AliError("RecoParam does not exist\n");
644 // Threshold value for the maximum
645 Float_t maxThresh = AliTRDReconstructor::RecoParam()->GetClusMaxThresh();
646 // Threshold value for the digit signal
647 Float_t sigThresh = AliTRDReconstructor::RecoParam()->GetClusSigThresh();
649 // Iteration limit for unfolding procedure
650 const Float_t kEpsilon = 0.01;
651 const Int_t kNclus = 3;
652 const Int_t kNsig = 5;
655 Double_t ratioLeft = 1.0;
656 Double_t ratioRight = 1.0;
658 Double_t padSignal[kNsig];
659 Double_t clusterSignal[kNclus];
661 Int_t icham = indexesIn->GetChamber();
662 Int_t iplan = indexesIn->GetPlane();
663 Int_t isect = indexesIn->GetSM();
665 // Start clustering in the chamber
667 Int_t idet = AliTRDgeometry::GetDetector(iplan,icham,isect);
670 AliError("Strange Detector number Missmatch!");
674 // TRD space point transformation
675 fTransform->SetDetector(det);
677 Int_t ilayer = AliGeomManager::kTRD1 + iplan;
678 Int_t imodule = icham + AliTRDgeometry::Ncham() * isect;
679 UShort_t volid = AliGeomManager::LayerToVolUID(ilayer,imodule);
681 Int_t nColMax = digitsIn->GetNcol();
682 Int_t nRowMax = digitsIn->GetNrow();
683 Int_t nTimeTotal = digitsIn->GetNtime();
685 // Detector wise calibration object for the gain factors
686 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
687 // Calibration object with pad wise values for the gain factors
688 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
689 // Calibration value for chamber wise gain factor
690 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
694 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(nRowMax, nColMax, nTimeTotal);
695 AliTRDdataArrayS padStatus(nRowMax, nColMax, nTimeTotal);
697 ResetHelperIndexes(indexesIn);
699 // Apply the gain and the tail cancelation via digital filter
700 TailCancelation(digitsIn
707 ,calGainFactorDetValue);
714 UChar_t status[3]={0, 0, 0}, ipos = 0;
715 fIndexesOut->ResetCounters();
716 while (fIndexesOut->NextRCTbinIndex(row, col, time)){
717 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
718 status[1] = digitsIn->GetPadStatus(row,col,time);
719 if(status[1]) SETBIT(ipos, AliTRDcluster::kMaskedCenter);
721 // Look for the maximum
722 if (signalM >= maxThresh) {
723 if (col + 1 >= nColMax || col-1 < 0) continue;
725 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
726 status[0] = digitsIn->GetPadStatus(row,col+1,time);
727 if(status[0]) SETBIT(ipos, AliTRDcluster::kMaskedLeft);
729 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
730 status[2] = digitsIn->GetPadStatus(row,col-1,time);
731 if(status[2]) SETBIT(ipos, AliTRDcluster::kMaskedRight);
733 // reject candidates with more than 1 problematic pad
734 if(ipos == 3 || ipos > 4) continue;
736 if(!status[1]){ // good central pad
737 if(!ipos){ // all pads are OK
738 if ((signalL <= signalM) && (signalR < signalM)) {
739 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
740 // Maximum found, mark the position by a negative signal
741 digitsOut->SetDataUnchecked(row,col,time,-signalM);
742 fIndexesMaxima->AddIndexTBin(row,col,time);
743 padStatus.SetDataUnchecked(row, col, time, ipos);
746 } else { // one of the neighbouring pads are bad
747 if(status[0] && signalR < signalM && signalR >= sigThresh){
748 digitsOut->SetDataUnchecked(row,col,time,-signalM);
749 digitsOut->SetDataUnchecked(row, col, time+1, 0.);
750 fIndexesMaxima->AddIndexTBin(row,col,time);
751 padStatus.SetDataUnchecked(row, col, time, ipos);
752 } else if(status[2] && signalL <= signalM && signalL >= sigThresh){
753 digitsOut->SetDataUnchecked(row,col,time,-signalM);
754 digitsOut->SetDataUnchecked(row, col, time-1, 0.);
755 fIndexesMaxima->AddIndexTBin(row,col,time);
756 padStatus.SetDataUnchecked(row, col, time, ipos);
759 } else { // wrong maximum pad
760 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
761 // Maximum found, mark the position by a negative signal
762 digitsOut->SetDataUnchecked(row,col,time,-maxThresh);
763 fIndexesMaxima->AddIndexTBin(row,col,time);
764 padStatus.SetDataUnchecked(row, col, time, ipos);
770 // The index to the first cluster of a given ROC
771 Int_t firstClusterROC = -1;
772 // The number of cluster in a given ROC
773 Int_t nClusterROC = 0;
775 // Now check the maxima and calculate the cluster position
776 fIndexesMaxima->ResetCounters();
777 while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) {
780 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
781 for (iPad = 0; iPad < kNclus; iPad++) {
782 Int_t iPadCol = col - 1 + iPad;
783 clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
786 // Count the number of pads in the cluster
791 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) {
794 if (col-ii < 0) break;
798 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh){
801 if (col+ii+1 >= nColMax) break;
805 // Look for 5 pad cluster with minimum in the middle
806 Bool_t fivePadCluster = kFALSE;
807 if (col < (nColMax - 3)){
808 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
809 fivePadCluster = kTRUE;
811 if ((fivePadCluster) && (col < (nColMax - 5))) {
812 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh){
813 fivePadCluster = kFALSE;
816 if ((fivePadCluster) && (col > 1)){
817 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh){
818 fivePadCluster = kFALSE;
824 // Modify the signal of the overlapping pad for the left part
825 // of the cluster which remains from a previous unfolding
827 clusterSignal[0] *= ratioLeft;
831 // Unfold the 5 pad cluster
833 for (iPad = 0; iPad < kNsig; iPad++) {
834 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
838 // Unfold the two maxima and set the signal on
839 // the overlapping pad to the ratio
840 ratioRight = Unfold(kEpsilon,iplan,padSignal);
841 ratioLeft = 1.0 - ratioRight;
842 clusterSignal[2] *= ratioRight;
846 // The position of the cluster in COL direction relative to the center pad (pad units)
847 Double_t clusterPosCol = 0.0;
848 if (AliTRDReconstructor::RecoParam()->LUTOn()) {
849 // Calculate the position of the cluster by using the
850 // lookup table method
851 clusterPosCol = LUTposition(iplan,clusterSignal[0]
855 // Calculate the position of the cluster by using the
856 // center of gravity method
857 for (Int_t i = 0; i < kNsig; i++) {
860 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
861 padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left pad
862 padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right pad
865 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
866 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
868 if ((col < nColMax - 3) &&
869 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])){
870 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
872 clusterPosCol = GetCOG(padSignal);
875 // Store the amplitudes of the pads in the cluster for later analysis
876 // and check whether one of these pads is masked in the database
877 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
878 for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
880 (jPad >= nColMax-1)) {
883 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
886 // Transform the local cluster coordinates into calibrated
887 // space point positions defined in the local tracking system.
888 // Here the calibration for T0, Vdrift and ExB is applied as well.
889 Double_t clusterXYZ[6];
890 clusterXYZ[0] = clusterPosCol;
891 clusterXYZ[1] = clusterSignal[0];
892 clusterXYZ[2] = clusterSignal[1];
893 clusterXYZ[3] = clusterSignal[2];
902 if (fTransform->Transform(clusterXYZ, clusterRCT, ((UInt_t) time), out, 0)) {
904 // Add the cluster to the output array
905 // The track indices will be stored later
906 Float_t clusterPos[3];
907 clusterPos[0] = clusterXYZ[0];
908 clusterPos[1] = clusterXYZ[1];
909 clusterPos[2] = clusterXYZ[2];
910 Float_t clusterSig[2];
911 clusterSig[0] = clusterXYZ[4];
912 clusterSig[1] = clusterXYZ[5];
913 Double_t clusterCharge = clusterXYZ[3];
914 Char_t clusterTimeBin = ((Char_t) clusterRCT[2]);
915 AliTRDcluster *cluster = new AliTRDcluster(idet
920 ,((Char_t) nPadCount)
928 cluster->SetInChamber(!out);
929 UChar_t maskPosition = padStatus.GetDataUnchecked(row, col, time);
931 cluster->SetPadMaskedPosition(maskPosition);
933 if(maskPosition & AliTRDcluster::kMaskedLeft) cluster->SetPadMaskedStatus(status[0]);
934 else if(maskPosition & AliTRDcluster::kMaskedCenter) cluster->SetPadMaskedStatus(status[1]);
935 else cluster->SetPadMaskedStatus(status[2]);
938 // Temporarily store the row, column and time bin of the center pad
939 // Used to later on assign the track indices
940 cluster->SetLabel( row,0);
941 cluster->SetLabel( col,1);
942 cluster->SetLabel(time,2);
944 RecPoints()->Add(cluster);
946 // Store the index of the first cluster in the current ROC
947 if (firstClusterROC < 0) {
948 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
951 // Count the number of cluster in the current ROC
954 } // if: Transform ok ?
955 } // if: Maximum found ?
960 if (fAddLabels) AddLabels(idet, firstClusterROC, nClusterROC);
962 // Write the cluster and reset the array
970 //_____________________________________________________________________________
971 Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
974 // Add the track indices to the found clusters
977 const Int_t kNclus = 3;
978 const Int_t kNdict = AliTRDdigitsManager::kNDict;
979 const Int_t kNtrack = kNdict * kNclus;
981 Int_t iClusterROC = 0;
988 // Temporary array to collect the track indices
989 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
991 // Loop through the dictionary arrays one-by-one
992 // to keep memory consumption low
993 AliTRDdataArrayI *tracksIn = 0;
994 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
996 // tracksIn should be expanded beforehand!
997 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
999 // Loop though the clusters found in this ROC
1000 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1002 AliTRDcluster *cluster = (AliTRDcluster *)
1003 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1004 row = cluster->GetLabel(0);
1005 col = cluster->GetLabel(1);
1006 time = cluster->GetLabel(2);
1008 for (iPad = 0; iPad < kNclus; iPad++) {
1009 Int_t iPadCol = col - 1 + iPad;
1010 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1011 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1018 // Copy the track indices into the cluster
1019 // Loop though the clusters found in this ROC
1020 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1022 AliTRDcluster *cluster = (AliTRDcluster *)
1023 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1024 cluster->SetLabel(-9999,0);
1025 cluster->SetLabel(-9999,1);
1026 cluster->SetLabel(-9999,2);
1028 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1032 delete [] idxTracks;
1038 //_____________________________________________________________________________
1039 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
1043 // Used for clusters with more than 3 pads - where LUT not applicable
1046 Double_t sum = signal[0]
1052 Double_t res = (0.0 * (-signal[0] + signal[4])
1053 + (-signal[1] + signal[3])) / sum;
1059 //_____________________________________________________________________________
1060 Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
1063 // Method to unfold neighbouring maxima.
1064 // The charge ratio on the overlapping pad is calculated
1065 // until there is no more change within the range given by eps.
1066 // The resulting ratio is then returned to the calling method.
1069 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1071 AliError("No AliTRDcalibDB instance available\n");
1076 Int_t itStep = 0; // Count iteration steps
1078 Double_t ratio = 0.5; // Start value for ratio
1079 Double_t prevRatio = 0.0; // Store previous ratio
1081 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1082 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1083 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1085 // Start the iteration
1086 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1091 // Cluster position according to charge ratio
1092 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1093 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1094 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1095 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1097 // Set cluster charge ratio
1098 irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
1099 Double_t ampLeft = padSignal[1] / newSignal[1];
1100 irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
1101 Double_t ampRight = padSignal[3] / newSignal[1];
1103 // Apply pad response to parameters
1104 irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
1105 irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
1107 // Calculate new overlapping ratio
1108 ratio = TMath::Min((Double_t) 1.0
1109 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1117 //_____________________________________________________________________________
1118 void AliTRDclusterizer::TailCancelation(AliTRDdataArrayDigits *digitsIn
1119 , AliTRDdataArrayF *digitsOut
1120 , AliTRDSignalIndex *indexesIn
1121 , AliTRDSignalIndex *indexesOut
1123 , Float_t adcThreshold
1124 , AliTRDCalROC *calGainFactorROC
1125 , Float_t calGainFactorDetValue)
1128 // Applies the tail cancelation and gain factors:
1129 // Transform digitsIn to digitsOut
1136 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
1137 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
1138 indexesIn->ResetCounters();
1139 while (indexesIn->NextRCIndex(iRow, iCol))
1141 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1142 Double_t gain = calGainFactorDetValue
1143 * calGainFactorROCValue;
1145 Bool_t corrupted = kFALSE;
1146 for (iTime = 0; iTime < nTimeTotal; iTime++)
1148 // Apply gain gain factor
1149 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1150 if(digitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
1151 inADC[iTime] /= gain;
1152 outADC[iTime] = inADC[iTime];
1156 // Apply the tail cancelation via the digital filter
1157 // (only for non-coorupted pads)
1158 if (AliTRDReconstructor::RecoParam()->TCOn())
1160 DeConvExp(inADC,outADC,nTimeTotal,AliTRDReconstructor::RecoParam()->GetTCnexp());
1164 indexesIn->ResetTbinCounter();
1165 while (indexesIn->NextTbinIndex(iTime))
1167 // Store the amplitude of the digit if above threshold
1168 if (outADC[iTime] > adcThreshold)
1170 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1171 indexesOut->AddIndexTBin(iRow,iCol,iTime);
1175 } // while irow icol
1184 //_____________________________________________________________________________
1185 void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1186 , Int_t n, Int_t nexp)
1189 // Tail cancellation by deconvolution for PASA v4 TRF
1193 Double_t coefficients[2];
1195 // Initialization (coefficient = alpha, rates = lambda)
1201 if (nexp == 1) { // 1 Exponentials
1207 if (nexp == 2) { // 2 Exponentials
1214 coefficients[0] = c1;
1215 coefficients[1] = c2;
1219 rates[0] = TMath::Exp(-dt/(r1));
1220 rates[1] = TMath::Exp(-dt/(r2));
1225 Double_t reminder[2];
1226 Double_t correction = 0.0;
1227 Double_t result = 0.0;
1229 // Attention: computation order is important
1230 for (k = 0; k < nexp; k++) {
1234 for (i = 0; i < n; i++) {
1236 result = (source[i] - correction); // No rescaling
1239 for (k = 0; k < nexp; k++) {
1240 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1244 for (k = 0; k < nexp; k++) {
1245 correction += reminder[k];
1252 //_____________________________________________________________________________
1253 void AliTRDclusterizer::ResetRecPoints()
1256 // Resets the list of rec points
1260 fRecPoints->Delete();
1265 //_____________________________________________________________________________
1266 TObjArray *AliTRDclusterizer::RecPoints()
1269 // Returns the list of rec points
1273 fRecPoints = new TObjArray(400);
1280 //_____________________________________________________________________________
1281 void AliTRDclusterizer::FillLUT()
1287 const Int_t kNlut = 128;
1289 fLUTbin = AliTRDgeometry::kNplan * kNlut;
1291 // The lookup table from Bogdan
1292 Float_t lut[AliTRDgeometry::kNplan][kNlut] = {
1294 0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1295 0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1296 0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1297 0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1298 0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1299 0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1300 0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1301 0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1302 0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1303 0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1304 0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1305 0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1306 0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1307 0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1308 0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1309 0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1312 0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1313 0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1314 0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1315 0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1316 0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1317 0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1318 0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1319 0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1320 0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1321 0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1322 0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1323 0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1324 0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1325 0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1326 0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1327 0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1330 0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1331 0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1332 0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1333 0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1334 0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1335 0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1336 0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1337 0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1338 0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1339 0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1340 0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1341 0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1342 0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1343 0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1344 0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1345 0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1348 0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1349 0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1350 0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1351 0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1352 0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1353 0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1354 0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1355 0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1356 0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1357 0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1358 0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1359 0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1360 0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1361 0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1362 0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1363 0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1366 0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1367 0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1368 0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1369 0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1370 0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1371 0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1372 0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1373 0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1374 0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1375 0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1376 0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1377 0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1378 0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1379 0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1380 0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1381 0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1384 0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1385 0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1386 0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1387 0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1388 0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1389 0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1390 0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1391 0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1392 0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1393 0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1394 0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1395 0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1396 0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1397 0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1398 0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1399 0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1406 fLUT = new Double_t[fLUTbin];
1408 for (Int_t iplan = 0; iplan < AliTRDgeometry::kNplan; iplan++) {
1409 for (Int_t ilut = 0; ilut < kNlut; ilut++ ) {
1410 fLUT[iplan*kNlut+ilut] = lut[iplan][ilut];
1416 //_____________________________________________________________________________
1417 Double_t AliTRDclusterizer::LUTposition(Int_t iplane, Double_t ampL
1418 , Double_t ampC, Double_t ampR) const
1421 // Calculates the cluster position using the lookup table.
1422 // Method provided by Bogdan Vulpescu.
1425 const Int_t kNlut = 128;
1436 Double_t xMin[AliTRDgeometry::kNplan] = { 0.006492, 0.006377, 0.006258
1437 , 0.006144, 0.006030, 0.005980 };
1438 Double_t xMax[AliTRDgeometry::kNplan] = { 0.960351, 0.965870, 0.970445
1439 , 0.974352, 0.977667, 0.996101 };
1442 x = (ampL - ampR) / ampC;
1445 else if (ampL < ampR) {
1446 x = (ampR - ampL) / ampC;
1452 xmin = xMin[iplane] + 0.000005;
1453 xmax = xMax[iplane] - 0.000005;
1454 xwid = (xmax - xmin) / 127.0;
1459 else if (x > xmax) {
1460 pos = side * 0.5000;
1463 ix = (Int_t) ((x - xmin) / xwid);
1464 pos = side * fLUT[iplane*kNlut+ix];