1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // TRD cluster finder //
22 ///////////////////////////////////////////////////////////////////////////////
28 #include <TObjArray.h>
30 #include "AliRunLoader.h"
31 #include "AliLoader.h"
32 #include "AliRawReader.h"
34 #include "AliAlignObj.h"
36 #include "AliTRDclusterizer.h"
37 #include "AliTRDcluster.h"
38 #include "AliTRDReconstructor.h"
39 #include "AliTRDgeometry.h"
40 #include "AliTRDdataArrayF.h"
41 #include "AliTRDdataArrayI.h"
42 #include "AliTRDdataArrayS.h"
43 #include "AliTRDdataArrayDigits.h"
44 #include "AliTRDdigitsManager.h"
45 #include "AliTRDrawData.h"
46 #include "AliTRDcalibDB.h"
47 #include "AliTRDrecoParam.h"
48 #include "AliTRDCommonParam.h"
49 #include "AliTRDtransform.h"
50 #include "AliTRDSignalIndex.h"
51 #include "AliTRDrawStreamBase.h"
52 #include "AliTRDfeeParam.h"
54 #include "Cal/AliTRDCalROC.h"
55 #include "Cal/AliTRDCalDet.h"
56 #include "Cal/AliTRDCalSingleChamberStatus.h"
58 ClassImp(AliTRDclusterizer)
60 //_____________________________________________________________________________
61 AliTRDclusterizer::AliTRDclusterizer()
68 ,fTrackletContainer(NULL)
73 ,fTransform(new AliTRDtransform(0))
78 // AliTRDclusterizer default constructor
81 AliTRDcalibDB *trd = 0x0;
82 if (!(trd = AliTRDcalibDB::Instance())) {
83 AliFatal("Could not get calibration object");
86 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
88 // retrive reco params
89 AliTRDrecoParam *rec = 0x0;
90 if (!(rec = AliTRDReconstructor::RecoParam())){
91 if(!(rec = trd->GetRecoParam(0))){
92 AliInfo("Using default RecoParams = LowFluxParam.");
93 rec = AliTRDrecoParam::GetLowFluxParam();
95 AliTRDReconstructor::SetRecoParam(rec);
100 //_____________________________________________________________________________
101 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title)
107 ,fDigitsManager(new AliTRDdigitsManager())
108 ,fTrackletContainer(NULL)
112 ,fIndexesMaxima(NULL)
113 ,fTransform(new AliTRDtransform(0))
118 // AliTRDclusterizer constructor
121 AliTRDcalibDB *trd = 0x0;
122 if (!(trd = AliTRDcalibDB::Instance())) {
123 AliFatal("Could not get calibration object");
126 // retrive reco params
127 AliTRDrecoParam *rec = 0x0;
128 if (!(rec = AliTRDReconstructor::RecoParam())){
129 if(!(rec = trd->GetRecoParam(0))){
130 AliInfo("Using default RecoParams = LowFluxParam.");
131 rec = AliTRDrecoParam::GetLowFluxParam();
133 AliTRDReconstructor::SetRecoParam(rec);
136 fDigitsManager->CreateArrays();
138 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
144 //_____________________________________________________________________________
145 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
151 ,fTrackletContainer(NULL)
152 ,fDigitsManager(NULL)
156 ,fIndexesMaxima(NULL)
162 // AliTRDclusterizer copy constructor
169 //_____________________________________________________________________________
170 AliTRDclusterizer::~AliTRDclusterizer()
173 // AliTRDclusterizer destructor
178 fRecPoints->Delete();
184 delete fDigitsManager;
185 fDigitsManager = NULL;
188 if (fTrackletContainer)
190 delete fTrackletContainer;
191 fTrackletContainer = NULL;
202 delete fIndexesMaxima;
203 fIndexesMaxima = NULL;
220 //_____________________________________________________________________________
221 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
224 // Assignment operator
229 ((AliTRDclusterizer &) c).Copy(*this);
236 //_____________________________________________________________________________
237 void AliTRDclusterizer::Copy(TObject &c) const
243 ((AliTRDclusterizer &) c).fClusterTree = NULL;
244 ((AliTRDclusterizer &) c).fRecPoints = NULL;
245 ((AliTRDclusterizer &) c).fTrackletTree = NULL;
246 ((AliTRDclusterizer &) c).fDigitsManager = NULL;
247 ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
248 ((AliTRDclusterizer &) c).fAddLabels = fAddLabels;
249 ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
250 ((AliTRDclusterizer &) c).fIndexesOut = NULL;
251 ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
252 ((AliTRDclusterizer &) c).fTransform = NULL;
253 ((AliTRDclusterizer &) c).fLUTbin = 0;
254 ((AliTRDclusterizer &) c).fLUT = NULL;
258 //_____________________________________________________________________________
259 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
262 // Opens the AliROOT file. Output and input are in the same file
265 TString evfoldname = AliConfig::GetDefaultEventFolderName();
266 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
269 fRunLoader = AliRunLoader::Open(name);
273 AliError(Form("Can not open session for file %s.",name));
284 //_____________________________________________________________________________
285 Bool_t AliTRDclusterizer::OpenOutput()
288 // Open the output file
291 TObjArray *ioArray = 0;
293 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
294 loader->MakeTree("R");
296 fClusterTree = loader->TreeR();
297 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
303 //_____________________________________________________________________________
304 Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
307 // Connect the output tree
310 TObjArray *ioArray = 0;
312 fClusterTree = clusterTree;
313 fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
317 if (AliTRDReconstructor::RecoParam()->IsTrackletWriteEnabled()){
318 TString evfoldname = AliConfig::GetDefaultEventFolderName();
319 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
322 fRunLoader = AliRunLoader::Open("galice.root");
325 AliError(Form("Can not open session for file galice.root."));
329 UInt_t **leaves = new UInt_t *[2];
330 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
332 AliError("Could not get the tracklets data loader!");
333 AliDataLoader *dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
334 fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
337 fTrackletTree = dl->Tree();
341 fTrackletTree = dl->Tree();
343 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
345 fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
353 //_____________________________________________________________________________
354 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
357 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
360 // Import the Trees for the event nEvent in the file
361 fRunLoader->GetEvent(nEvent);
367 //_____________________________________________________________________________
368 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
371 // Fills TRDcluster branch in the tree with the clusters
372 // found in detector = det. For det=-1 writes the tree.
376 (det >= AliTRDgeometry::Ndet())) {
377 AliError(Form("Unexpected detector index %d.\n",det));
381 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
383 TObjArray *ioArray = 0;
384 branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
388 (det < AliTRDgeometry::Ndet())) {
390 Int_t nRecPoints = RecPoints()->GetEntriesFast();
391 TObjArray *detRecPoints = new TObjArray(400);
393 for (Int_t i = 0; i < nRecPoints; i++) {
394 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
395 if (det == c->GetDetector()) {
396 detRecPoints->AddLast(c);
399 AliError(Form("Attempt to write a cluster with unexpected detector index: got=%d expected=%d\n"
405 branch->SetAddress(&detRecPoints);
406 fClusterTree->Fill();
416 AliInfo(Form("Writing the cluster tree %s for event %d."
417 ,fClusterTree->GetName(),fRunLoader->GetEventNumber()));
421 branch->SetAddress(&fRecPoints);
423 AliLoader *loader = fRunLoader->GetLoader("TRDLoader");
424 loader->WriteRecPoints("OVERWRITE");
429 AliError("Cluster tree does not exist. Cannot write clusters.\n");
438 AliError(Form("Unexpected detector index %d.\n",det));
444 //_____________________________________________________________________________
445 Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
448 UInt_t **leaves = new UInt_t *[2];
449 for (Int_t i=0; i<2 ;i++){
450 leaves[i] = new UInt_t[258];
451 leaves[i][0] = det; // det
452 leaves[i][1] = i; // side
453 memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
458 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
460 fTrackletTree = dl->Tree();
463 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
465 trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
468 for (Int_t i=0; i<2; i++){
469 if (leaves[i][2]>0) {
470 trkbranch->SetAddress(leaves[i]);
471 fTrackletTree->Fill();
475 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
476 dl->WriteData("OVERWRITE");
483 //_____________________________________________________________________________
484 void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
487 // Reset the helper indexes
492 // carefull here - we assume that only row number may change - most probable
493 if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
494 fIndexesOut->ResetContent();
496 fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
497 , indexesIn->GetNcol()
498 , indexesIn->GetNtime());
502 fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
503 , indexesIn->GetNcol()
504 , indexesIn->GetNtime());
509 // carefull here - we assume that only row number may change - most probable
510 if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
512 fIndexesMaxima->ResetContent();
516 fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
517 , indexesIn->GetNcol()
518 , indexesIn->GetNtime());
523 fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
524 , indexesIn->GetNcol()
525 , indexesIn->GetNtime());
530 //_____________________________________________________________________________
531 Bool_t AliTRDclusterizer::ReadDigits()
534 // Reads the digits arrays from the input aliroot file
538 AliError("No run loader available");
542 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
543 if (!loader->TreeD()) {
544 loader->LoadDigits();
547 // Read in the digit arrays
548 return (fDigitsManager->ReadDigits(loader->TreeD()));
552 //_____________________________________________________________________________
553 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
556 // Reads the digits arrays from the input tree
559 // Read in the digit arrays
560 return (fDigitsManager->ReadDigits(digitsTree));
564 //_____________________________________________________________________________
565 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
568 // Reads the digits arrays from the ddl file
572 fDigitsManager = raw.Raw2Digits(rawReader);
578 //_____________________________________________________________________________
579 Bool_t AliTRDclusterizer::MakeClusters()
582 // Creates clusters from digits
585 // Propagate info from the digits manager
586 if (fAddLabels == kTRUE)
588 fAddLabels = fDigitsManager->UsesDictionaries();
591 Bool_t fReturn = kTRUE;
592 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++)
595 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits*) fDigitsManager->GetDigits(i);
596 // This is to take care of switched off super modules
597 if (!digitsIn->HasData())
602 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
603 if (indexes->IsAllocated() == kFALSE)
605 fDigitsManager->BuildIndexes(i);
609 if (indexes->HasEntry())
613 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++)
615 AliTRDdataArrayI *tracksIn = 0;
616 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(i,iDict);
620 fR = MakeClusters(i);
621 fReturn = fR && fReturn;
630 // No compress just remove
631 fDigitsManager->RemoveDigits(i);
632 fDigitsManager->RemoveDictionaries(i);
633 fDigitsManager->ClearIndexes(i);
641 //_____________________________________________________________________________
642 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
645 // Creates clusters from raw data
648 return Raw2ClustersChamber(rawReader);
652 //_____________________________________________________________________________
653 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
656 // Creates clusters from raw data
659 // Create the digits manager
662 fDigitsManager = new AliTRDdigitsManager();
663 fDigitsManager->CreateArrays();
666 fDigitsManager->SetUseDictionaries(fAddLabels);
668 // tracklet container for raw tracklet writing
669 if (!fTrackletContainer && AliTRDReconstructor::RecoParam()->IsTrackletWriteEnabled())
671 fTrackletContainer = new UInt_t *[2];
672 for (Int_t i=0; i<2 ;i++){
673 fTrackletContainer[i] = new UInt_t[256]; // maximum tracklets for one HC
677 AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
678 AliTRDrawStreamBase &input = *pinput;
680 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
683 while ((det = input.NextChamber(fDigitsManager,fTrackletContainer)) >= 0)
685 Bool_t iclusterBranch = kFALSE;
686 if (fDigitsManager->GetIndexes(det)->HasEntry())
688 iclusterBranch = MakeClusters(det);
690 if (iclusterBranch == kFALSE)
695 fDigitsManager->RemoveDigits(det);
696 fDigitsManager->RemoveDictionaries(det);
697 fDigitsManager->ClearIndexes(det);
699 if (!AliTRDReconstructor::RecoParam()-> IsTrackletWriteEnabled()) continue;
700 if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det); // if there is tracklet words in this det
703 if (AliTRDReconstructor::RecoParam()->IsTrackletWriteEnabled()){
704 delete [] fTrackletContainer;
705 fTrackletContainer = NULL;
708 delete fDigitsManager;
709 fDigitsManager = NULL;
717 //_____________________________________________________________________________
718 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
721 // Check if a pad is masked
726 if(signal>0 && TESTBIT(signal, 10)){
728 for(int ibit=0; ibit<4; ibit++){
729 if(TESTBIT(signal, 11+ibit)){
730 SETBIT(status, ibit);
731 CLRBIT(signal, 11+ibit);
738 //_____________________________________________________________________________
739 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
742 // Generates the cluster.
746 // digits should be expanded beforehand!
747 // digitsIn->Expand();
748 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits *) fDigitsManager->GetDigits(det);
750 // This is to take care of switched off super modules
751 if (!digitsIn->HasData())
756 AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
757 if (indexesIn->IsAllocated() == kFALSE)
759 AliError("Indexes do not exist!");
763 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
766 AliFatal("No AliTRDcalibDB instance available\n");
771 // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
772 Float_t adcThreshold = 0;
774 if (!AliTRDReconstructor::RecoParam())
776 AliError("RecoParam does not exist\n");
780 // Threshold value for the maximum
781 Float_t maxThresh = AliTRDReconstructor::RecoParam()->GetClusMaxThresh();
782 // Threshold value for the digit signal
783 Float_t sigThresh = AliTRDReconstructor::RecoParam()->GetClusSigThresh();
785 // Threshold value for the maximum ( cut noise)
786 Float_t minMaxCutSigma = AliTRDReconstructor::RecoParam()->GetMinMaxCutSigma();
787 // Threshold value for the sum pad ( cut noise)
788 Float_t minLeftRightCutSigma = AliTRDReconstructor::RecoParam()->GetMinLeftRightCutSigma();
790 // Iteration limit for unfolding procedure
791 const Float_t kEpsilon = 0.01;
792 const Int_t kNclus = 3;
793 const Int_t kNsig = 5;
796 Double_t ratioLeft = 1.0;
797 Double_t ratioRight = 1.0;
799 Double_t padSignal[kNsig];
800 Double_t clusterSignal[kNclus];
802 Int_t istack = indexesIn->GetStack();
803 Int_t ilayer = indexesIn->GetLayer();
804 Int_t isector = indexesIn->GetSM();
806 // Start clustering in the chamber
808 Int_t idet = AliTRDgeometry::GetDetector(ilayer,istack,isector);
811 AliError("Strange Detector number Missmatch!");
815 // TRD space point transformation
816 fTransform->SetDetector(det);
818 Int_t iGeoLayer = AliGeomManager::kTRD1 + ilayer;
819 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
820 UShort_t volid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
822 Int_t nColMax = digitsIn->GetNcol();
823 Int_t nRowMax = digitsIn->GetNrow();
824 Int_t nTimeTotal = digitsIn->GetNtime();
826 // Detector wise calibration object for the gain factors
827 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
828 // Calibration object with pad wise values for the gain factors
829 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
830 // Calibration value for chamber wise gain factor
831 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
834 // Detector wise calibration object for the noise
835 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
836 // Calibration object with pad wise values for the noise
837 AliTRDCalROC *calNoiseROC = calibration->GetNoiseROC(idet);
838 // Calibration value for chamber wise noise
839 Float_t calNoiseDetValue = calNoiseDet->GetValue(idet);
843 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(nRowMax, nColMax, nTimeTotal);
844 AliTRDdataArrayS padStatus(nRowMax, nColMax, nTimeTotal);
846 ResetHelperIndexes(indexesIn);
848 // Apply the gain and the tail cancelation via digital filter
849 TailCancelation(digitsIn
856 ,calGainFactorDetValue);
863 UChar_t status[3]={0, 0, 0}, ipos = 0;
864 fIndexesOut->ResetCounters();
865 while (fIndexesOut->NextRCTbinIndex(row, col, time)){
867 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
868 status[1] = digitsIn->GetPadStatus(row,col,time);
869 if(status[1]) SETBIT(ipos, AliTRDcluster::kMaskedCenter);
871 if(signalM < maxThresh) continue;
873 Float_t noiseMiddleThresh = minMaxCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
874 if (signalM < noiseMiddleThresh) continue;
876 if (col + 1 >= nColMax || col-1 < 0) continue;
878 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
879 status[0] = digitsIn->GetPadStatus(row,col+1,time);
880 if(status[0]) SETBIT(ipos, AliTRDcluster::kMaskedLeft);
882 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
883 status[2] = digitsIn->GetPadStatus(row,col-1,time);
884 if(status[2]) SETBIT(ipos, AliTRDcluster::kMaskedRight);
886 // reject candidates with more than 1 problematic pad
887 if(ipos == 3 || ipos > 4) continue;
889 if(!status[1]){ // good central pad
890 if(!ipos){ // all pads are OK
891 if ((signalL <= signalM) && (signalR < signalM)) {
892 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
893 Float_t noiseSumThresh = minLeftRightCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
894 if((signalL+signalR+signalM) >= noiseSumThresh){
895 // Maximum found, mark the position by a negative signal
896 digitsOut->SetDataUnchecked(row,col,time,-signalM);
897 fIndexesMaxima->AddIndexTBin(row,col,time);
898 padStatus.SetDataUnchecked(row, col, time, ipos);
902 } else { // one of the neighbouring pads are bad
903 if(status[0] && signalR < signalM && signalR >= sigThresh){
904 digitsOut->SetDataUnchecked(row,col,time,-signalM);
905 digitsOut->SetDataUnchecked(row, col, time+1, 0.);
906 fIndexesMaxima->AddIndexTBin(row,col,time);
907 padStatus.SetDataUnchecked(row, col, time, ipos);
908 } else if(status[2] && signalL <= signalM && signalL >= sigThresh){
909 digitsOut->SetDataUnchecked(row,col,time,-signalM);
910 digitsOut->SetDataUnchecked(row, col, time-1, 0.);
911 fIndexesMaxima->AddIndexTBin(row,col,time);
912 padStatus.SetDataUnchecked(row, col, time, ipos);
915 } else { // wrong maximum pad
916 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
917 // Maximum found, mark the position by a negative signal
918 digitsOut->SetDataUnchecked(row,col,time,-maxThresh);
919 fIndexesMaxima->AddIndexTBin(row,col,time);
920 padStatus.SetDataUnchecked(row, col, time, ipos);
925 // The index to the first cluster of a given ROC
926 Int_t firstClusterROC = -1;
927 // The number of cluster in a given ROC
928 Int_t nClusterROC = 0;
930 // Now check the maxima and calculate the cluster position
931 fIndexesMaxima->ResetCounters();
932 while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) {
935 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
937 for (iPad = 0; iPad < kNclus; iPad++) {
938 Int_t iPadCol = col - 1 + iPad;
939 clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
942 // Count the number of pads in the cluster
947 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) {
950 if (col-ii < 0) break;
954 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh){
957 if (col+ii+1 >= nColMax) break;
961 // Look for 5 pad cluster with minimum in the middle
962 Bool_t fivePadCluster = kFALSE;
963 if (col < (nColMax - 3)){
964 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
965 fivePadCluster = kTRUE;
967 if ((fivePadCluster) && (col < (nColMax - 5))) {
968 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh){
969 fivePadCluster = kFALSE;
972 if ((fivePadCluster) && (col > 1)){
973 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh){
974 fivePadCluster = kFALSE;
980 // Modify the signal of the overlapping pad for the left part
981 // of the cluster which remains from a previous unfolding
983 clusterSignal[0] *= ratioLeft;
987 // Unfold the 5 pad cluster
989 for (iPad = 0; iPad < kNsig; iPad++) {
990 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
994 // Unfold the two maxima and set the signal on
995 // the overlapping pad to the ratio
996 ratioRight = Unfold(kEpsilon,ilayer,padSignal);
997 ratioLeft = 1.0 - ratioRight;
998 clusterSignal[2] *= ratioRight;
1002 // The position of the cluster in COL direction relative to the center pad (pad units)
1003 Double_t clusterPosCol = 0.0;
1004 if (AliTRDReconstructor::RecoParam()->IsLUT()) {
1005 // Calculate the position of the cluster by using the
1006 // lookup table method
1007 clusterPosCol = LUTposition(ilayer,clusterSignal[2]
1012 // Calculate the position of the cluster by using the
1013 // center of gravity method
1014 for (Int_t i = 0; i < kNsig; i++) {
1017 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
1018 padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Left pad
1019 padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Right pad
1021 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
1022 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
1024 if ((col < nColMax - 3) &&
1025 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])){
1026 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
1028 clusterPosCol = GetCOG(padSignal);
1031 // Store the amplitudes of the pads in the cluster for later analysis
1032 // and check whether one of these pads is masked in the database
1033 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
1034 for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
1036 (jPad >= nColMax-1)) {
1039 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
1042 // Transform the local cluster coordinates into calibrated
1043 // space point positions defined in the local tracking system.
1044 // Here the calibration for T0, Vdrift and ExB is applied as well.
1045 Double_t clusterXYZ[6];
1046 clusterXYZ[0] = clusterPosCol;
1047 clusterXYZ[1] = clusterSignal[2];
1048 clusterXYZ[2] = clusterSignal[1];
1049 clusterXYZ[3] = clusterSignal[0];
1050 clusterXYZ[4] = 0.0;
1051 clusterXYZ[5] = 0.0;
1052 Int_t clusterRCT[3];
1053 clusterRCT[0] = row;
1054 clusterRCT[1] = col;
1058 if (fTransform->Transform(clusterXYZ, clusterRCT, ((UInt_t) time), out, 0)) {
1060 // Add the cluster to the output array
1061 // The track indices will be stored later
1062 Float_t clusterPos[3];
1063 clusterPos[0] = clusterXYZ[0];
1064 clusterPos[1] = clusterXYZ[1];
1065 clusterPos[2] = clusterXYZ[2];
1066 Float_t clusterSig[2];
1067 clusterSig[0] = clusterXYZ[4];
1068 clusterSig[1] = clusterXYZ[5];
1069 Double_t clusterCharge = clusterXYZ[3];
1070 Char_t clusterTimeBin = ((Char_t) clusterRCT[2]);
1071 AliTRDcluster *cluster = new AliTRDcluster(idet
1076 ,((Char_t) nPadCount)
1084 cluster->SetInChamber(!out);
1086 UChar_t maskPosition = padStatus.GetDataUnchecked(row, col, time);
1088 cluster->SetPadMaskedPosition(maskPosition);
1089 if (maskPosition & AliTRDcluster::kMaskedLeft) {
1090 cluster->SetPadMaskedStatus(status[0]);
1092 else if (maskPosition & AliTRDcluster::kMaskedCenter) {
1093 cluster->SetPadMaskedStatus(status[1]);
1096 cluster->SetPadMaskedStatus(status[2]);
1100 // Temporarily store the row, column and time bin of the center pad
1101 // Used to later on assign the track indices
1102 cluster->SetLabel( row,0);
1103 cluster->SetLabel( col,1);
1104 cluster->SetLabel(time,2);
1106 RecPoints()->Add(cluster);
1108 // Store the index of the first cluster in the current ROC
1109 if (firstClusterROC < 0) {
1110 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
1113 // Count the number of cluster in the current ROC
1116 } // if: Transform ok ?
1117 } // if: Maximum found ?
1122 if (fAddLabels) AddLabels(idet, firstClusterROC, nClusterROC);
1124 // Write the cluster and reset the array
1125 WriteClusters(idet);
1132 //_____________________________________________________________________________
1133 Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
1136 // Add the track indices to the found clusters
1139 const Int_t kNclus = 3;
1140 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1141 const Int_t kNtrack = kNdict * kNclus;
1143 Int_t iClusterROC = 0;
1150 // Temporary array to collect the track indices
1151 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1153 // Loop through the dictionary arrays one-by-one
1154 // to keep memory consumption low
1155 AliTRDdataArrayI *tracksIn = 0;
1156 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1158 // tracksIn should be expanded beforehand!
1159 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
1161 // Loop though the clusters found in this ROC
1162 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1164 AliTRDcluster *cluster = (AliTRDcluster *)
1165 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1166 row = cluster->GetLabel(0);
1167 col = cluster->GetLabel(1);
1168 time = cluster->GetLabel(2);
1170 for (iPad = 0; iPad < kNclus; iPad++) {
1171 Int_t iPadCol = col - 1 + iPad;
1172 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1173 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1180 // Copy the track indices into the cluster
1181 // Loop though the clusters found in this ROC
1182 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1184 AliTRDcluster *cluster = (AliTRDcluster *)
1185 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1186 cluster->SetLabel(-9999,0);
1187 cluster->SetLabel(-9999,1);
1188 cluster->SetLabel(-9999,2);
1190 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1194 delete [] idxTracks;
1200 //_____________________________________________________________________________
1201 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
1205 // Used for clusters with more than 3 pads - where LUT not applicable
1208 Double_t sum = signal[0]
1215 Double_t res = (0.0 * (-signal[0] + signal[4])
1216 + (-signal[1] + signal[3])) / sum;
1222 //_____________________________________________________________________________
1223 Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal)
1226 // Method to unfold neighbouring maxima.
1227 // The charge ratio on the overlapping pad is calculated
1228 // until there is no more change within the range given by eps.
1229 // The resulting ratio is then returned to the calling method.
1232 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1234 AliError("No AliTRDcalibDB instance available\n");
1239 Int_t itStep = 0; // Count iteration steps
1241 Double_t ratio = 0.5; // Start value for ratio
1242 Double_t prevRatio = 0.0; // Store previous ratio
1244 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1245 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1246 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1248 // Start the iteration
1249 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1254 // Cluster position according to charge ratio
1255 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1256 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1257 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1258 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1260 // Set cluster charge ratio
1261 irc = calibration->PadResponse(1.0,maxLeft ,layer,newSignal);
1262 Double_t ampLeft = padSignal[1] / newSignal[1];
1263 irc = calibration->PadResponse(1.0,maxRight,layer,newSignal);
1264 Double_t ampRight = padSignal[3] / newSignal[1];
1266 // Apply pad response to parameters
1267 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1268 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1270 // Calculate new overlapping ratio
1271 ratio = TMath::Min((Double_t) 1.0
1272 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1280 //_____________________________________________________________________________
1281 void AliTRDclusterizer::TailCancelation(AliTRDdataArrayDigits *digitsIn
1282 , AliTRDdataArrayF *digitsOut
1283 , AliTRDSignalIndex *indexesIn
1284 , AliTRDSignalIndex *indexesOut
1286 , Float_t adcThreshold
1287 , AliTRDCalROC *calGainFactorROC
1288 , Float_t calGainFactorDetValue)
1291 // Applies the tail cancelation and gain factors:
1292 // Transform digitsIn to digitsOut
1299 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
1300 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
1301 indexesIn->ResetCounters();
1302 while (indexesIn->NextRCIndex(iRow, iCol))
1304 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1305 Double_t gain = calGainFactorDetValue
1306 * calGainFactorROCValue;
1308 Bool_t corrupted = kFALSE;
1309 for (iTime = 0; iTime < nTimeTotal; iTime++)
1311 // Apply gain gain factor
1312 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1313 if(digitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
1314 inADC[iTime] /= gain;
1315 outADC[iTime] = inADC[iTime];
1319 // Apply the tail cancelation via the digital filter
1320 // (only for non-coorupted pads)
1321 if (AliTRDReconstructor::RecoParam()->IsTailCancelation())
1323 DeConvExp(inADC,outADC,nTimeTotal,AliTRDReconstructor::RecoParam()->GetTCnexp());
1327 indexesIn->ResetTbinCounter();
1328 while (indexesIn->NextTbinIndex(iTime))
1330 // Store the amplitude of the digit if above threshold
1331 if (outADC[iTime] > adcThreshold)
1333 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1334 indexesOut->AddIndexTBin(iRow,iCol,iTime);
1338 } // while irow icol
1347 //_____________________________________________________________________________
1348 void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1349 , Int_t n, Int_t nexp)
1352 // Tail cancellation by deconvolution for PASA v4 TRF
1356 Double_t coefficients[2];
1358 // Initialization (coefficient = alpha, rates = lambda)
1364 if (nexp == 1) { // 1 Exponentials
1370 if (nexp == 2) { // 2 Exponentials
1377 coefficients[0] = c1;
1378 coefficients[1] = c2;
1382 rates[0] = TMath::Exp(-dt/(r1));
1383 rates[1] = TMath::Exp(-dt/(r2));
1388 Double_t reminder[2];
1389 Double_t correction = 0.0;
1390 Double_t result = 0.0;
1392 // Attention: computation order is important
1393 for (k = 0; k < nexp; k++) {
1397 for (i = 0; i < n; i++) {
1399 result = (source[i] - correction); // No rescaling
1402 for (k = 0; k < nexp; k++) {
1403 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1407 for (k = 0; k < nexp; k++) {
1408 correction += reminder[k];
1415 //_____________________________________________________________________________
1416 void AliTRDclusterizer::ResetRecPoints()
1419 // Resets the list of rec points
1423 fRecPoints->Delete();
1428 //_____________________________________________________________________________
1429 TObjArray *AliTRDclusterizer::RecPoints()
1432 // Returns the list of rec points
1436 fRecPoints = new TObjArray(400);
1443 //_____________________________________________________________________________
1444 void AliTRDclusterizer::FillLUT()
1450 const Int_t kNlut = 128;
1452 fLUTbin = AliTRDgeometry::kNlayer * kNlut;
1454 // The lookup table from Bogdan
1455 Float_t lut[AliTRDgeometry::kNlayer][kNlut] = {
1457 0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1458 0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1459 0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1460 0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1461 0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1462 0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1463 0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1464 0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1465 0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1466 0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1467 0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1468 0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1469 0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1470 0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1471 0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1472 0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1475 0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1476 0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1477 0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1478 0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1479 0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1480 0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1481 0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1482 0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1483 0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1484 0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1485 0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1486 0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1487 0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1488 0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1489 0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1490 0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1493 0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1494 0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1495 0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1496 0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1497 0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1498 0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1499 0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1500 0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1501 0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1502 0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1503 0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1504 0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1505 0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1506 0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1507 0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1508 0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1511 0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1512 0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1513 0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1514 0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1515 0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1516 0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1517 0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1518 0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1519 0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1520 0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1521 0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1522 0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1523 0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1524 0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1525 0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1526 0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1529 0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1530 0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1531 0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1532 0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1533 0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1534 0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1535 0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1536 0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1537 0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1538 0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1539 0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1540 0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1541 0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1542 0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1543 0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1544 0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1547 0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1548 0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1549 0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1550 0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1551 0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1552 0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1553 0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1554 0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1555 0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1556 0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1557 0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1558 0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1559 0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1560 0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1561 0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1562 0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1569 fLUT = new Double_t[fLUTbin];
1571 for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1572 for (Int_t ilut = 0; ilut < kNlut; ilut++ ) {
1573 fLUT[ilayer*kNlut+ilut] = lut[ilayer][ilut];
1579 //_____________________________________________________________________________
1580 Double_t AliTRDclusterizer::LUTposition(Int_t ilayer, Double_t ampL
1581 , Double_t ampC, Double_t ampR) const
1584 // Calculates the cluster position using the lookup table.
1585 // Method provided by Bogdan Vulpescu.
1588 const Int_t kNlut = 128;
1599 Double_t xMin[AliTRDgeometry::kNlayer] = { 0.006492, 0.006377, 0.006258
1600 , 0.006144, 0.006030, 0.005980 };
1601 Double_t xMax[AliTRDgeometry::kNlayer] = { 0.960351, 0.965870, 0.970445
1602 , 0.974352, 0.977667, 0.996101 };
1605 x = (ampL - ampR) / ampC;
1608 else if (ampL < ampR) {
1609 x = (ampR - ampL) / ampC;
1615 xmin = xMin[ilayer] + 0.000005;
1616 xmax = xMax[ilayer] - 0.000005;
1617 xwid = (xmax - xmin) / 127.0;
1622 else if (x > xmax) {
1623 pos = side * 0.5000;
1626 ix = (Int_t) ((x - xmin) / xwid);
1627 pos = side * fLUT[ilayer*kNlut+ix];