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 <TClonesArray.h>
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(AliTRDReconstructor *rec)
70 ,fTrackletContainer(NULL)
75 ,fTransform(new AliTRDtransform(0))
80 // AliTRDclusterizer default constructor
83 AliTRDcalibDB *trd = 0x0;
84 if (!(trd = AliTRDcalibDB::Instance())) {
85 AliFatal("Could not get calibration object");
88 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
90 // Initialize debug stream
92 if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 1){
93 TDirectory *savedir = gDirectory;
94 //fgDebugStreamer = new TTreeSRedirector("TRD.ClusterizerDebug.root");
101 //_____________________________________________________________________________
102 AliTRDclusterizer::AliTRDclusterizer(const Text_t *name, const Text_t *title, AliTRDReconstructor *rec)
109 ,fDigitsManager(new AliTRDdigitsManager())
110 ,fTrackletContainer(NULL)
114 ,fIndexesMaxima(NULL)
115 ,fTransform(new AliTRDtransform(0))
120 // AliTRDclusterizer constructor
123 AliTRDcalibDB *trd = 0x0;
124 if (!(trd = AliTRDcalibDB::Instance())) {
125 AliFatal("Could not get calibration object");
128 // Initialize debug stream
130 if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kClusterizer) > 1){
131 TDirectory *savedir = gDirectory;
132 //fgDebugStreamer = new TTreeSRedirector("TRD.ClusterizerDebug.root");
137 fDigitsManager->CreateArrays();
139 fRawVersion = AliTRDfeeParam::Instance()->GetRAWversion();
145 //_____________________________________________________________________________
146 AliTRDclusterizer::AliTRDclusterizer(const AliTRDclusterizer &c)
148 ,fReconstructor(c.fReconstructor)
153 ,fDigitsManager(NULL)
154 ,fTrackletContainer(NULL)
158 ,fIndexesMaxima(NULL)
164 // AliTRDclusterizer copy constructor
171 //_____________________________________________________________________________
172 AliTRDclusterizer::~AliTRDclusterizer()
175 // AliTRDclusterizer destructor
178 if (fRecPoints/* && IsClustersOwner()*/){
179 fRecPoints->Delete();
183 if (fDigitsManager) {
184 delete fDigitsManager;
185 fDigitsManager = NULL;
188 if (fTrackletContainer){
189 delete fTrackletContainer;
190 fTrackletContainer = NULL;
199 delete fIndexesMaxima;
200 fIndexesMaxima = NULL;
215 //_____________________________________________________________________________
216 AliTRDclusterizer &AliTRDclusterizer::operator=(const AliTRDclusterizer &c)
219 // Assignment operator
224 ((AliTRDclusterizer &) c).Copy(*this);
231 //_____________________________________________________________________________
232 void AliTRDclusterizer::Copy(TObject &c) const
238 ((AliTRDclusterizer &) c).fClusterTree = NULL;
239 ((AliTRDclusterizer &) c).fRecPoints = NULL;
240 ((AliTRDclusterizer &) c).fTrackletTree = NULL;
241 ((AliTRDclusterizer &) c).fDigitsManager = NULL;
242 ((AliTRDclusterizer &) c).fTrackletContainer = NULL;
243 ((AliTRDclusterizer &) c).fAddLabels = fAddLabels;
244 ((AliTRDclusterizer &) c).fRawVersion = fRawVersion;
245 ((AliTRDclusterizer &) c).fIndexesOut = NULL;
246 ((AliTRDclusterizer &) c).fIndexesMaxima = NULL;
247 ((AliTRDclusterizer &) c).fTransform = NULL;
248 ((AliTRDclusterizer &) c).fLUTbin = 0;
249 ((AliTRDclusterizer &) c).fLUT = NULL;
253 //_____________________________________________________________________________
254 Bool_t AliTRDclusterizer::Open(const Char_t *name, Int_t nEvent)
257 // Opens the AliROOT file. Output and input are in the same file
260 TString evfoldname = AliConfig::GetDefaultEventFolderName();
261 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
264 fRunLoader = AliRunLoader::Open(name);
268 AliError(Form("Can not open session for file %s.",name));
279 //_____________________________________________________________________________
280 Bool_t AliTRDclusterizer::OpenOutput()
283 // Open the output file
286 if (!fReconstructor->IsWritingClusters()) return kTRUE;
288 TObjArray *ioArray = 0x0;
290 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
291 loader->MakeTree("R");
293 fClusterTree = loader->TreeR();
294 fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
300 //_____________________________________________________________________________
301 Bool_t AliTRDclusterizer::OpenOutput(TTree *clusterTree)
304 // Connect the output tree
308 if (fReconstructor->IsWritingClusters()){
309 TObjArray *ioArray = 0x0;
310 fClusterTree = clusterTree;
311 fClusterTree->Branch("TRDcluster", "TObjArray", &ioArray, 32000, 0);
315 if (fReconstructor->IsWritingTracklets()){
316 TString evfoldname = AliConfig::GetDefaultEventFolderName();
317 fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
320 fRunLoader = AliRunLoader::Open("galice.root");
323 AliError(Form("Can not open session for file galice.root."));
327 UInt_t **leaves = new UInt_t *[2];
328 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
330 AliError("Could not get the tracklets data loader!");
331 dl = new AliDataLoader("TRD.Tracklets.root","tracklets", "tracklets");
332 fRunLoader->GetLoader("TRDLoader")->AddDataLoader(dl);
335 fTrackletTree = dl->Tree();
339 fTrackletTree = dl->Tree();
341 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
343 fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
351 //_____________________________________________________________________________
352 Bool_t AliTRDclusterizer::OpenInput(Int_t nEvent)
355 // Opens a ROOT-file with TRD-hits and reads in the digits-tree
358 // Import the Trees for the event nEvent in the file
359 fRunLoader->GetEvent(nEvent);
365 //_____________________________________________________________________________
366 Bool_t AliTRDclusterizer::WriteClusters(Int_t det)
369 // Fills TRDcluster branch in the tree with the clusters
370 // found in detector = det. For det=-1 writes the tree.
374 (det >= AliTRDgeometry::Ndet())) {
375 AliError(Form("Unexpected detector index %d.\n",det));
379 TObjArray *ioArray = new TObjArray(400);
380 TBranch *branch = fClusterTree->GetBranch("TRDcluster");
382 branch = fClusterTree->Branch("TRDcluster","TObjArray",&ioArray,32000,0);
383 } else branch->SetAddress(&ioArray);
385 Int_t nRecPoints = RecPoints()->GetEntriesFast();
387 for (Int_t i = 0; i < nRecPoints; i++) {
388 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
389 if(det != c->GetDetector()) continue;
392 fClusterTree->Fill();
396 for (Int_t i = 0; i < nRecPoints; i++) {
397 AliTRDcluster *c = (AliTRDcluster *) RecPoints()->UncheckedAt(i);
398 if(c->GetDetector() != detOld){
399 fClusterTree->Fill();
401 detOld = c->GetDetector();
412 //_____________________________________________________________________________
413 Bool_t AliTRDclusterizer::WriteTracklets(Int_t det)
416 // Write the raw data tracklets into seperate file
419 UInt_t **leaves = new UInt_t *[2];
420 for (Int_t i=0; i<2 ;i++){
421 leaves[i] = new UInt_t[258];
422 leaves[i][0] = det; // det
423 leaves[i][1] = i; // side
424 memcpy(leaves[i]+2, fTrackletContainer[i], sizeof(UInt_t) * 256);
428 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
430 fTrackletTree = dl->Tree();
433 TBranch *trkbranch = fTrackletTree->GetBranch("trkbranch");
435 trkbranch = fTrackletTree->Branch("trkbranch",leaves[0],"det/i:side/i:tracklets[256]/i");
438 for (Int_t i=0; i<2; i++){
439 if (leaves[i][2]>0) {
440 trkbranch->SetAddress(leaves[i]);
441 fTrackletTree->Fill();
445 AliDataLoader *dl = fRunLoader->GetLoader("TRDLoader")->GetDataLoader("tracklets");
446 dl->WriteData("OVERWRITE");
454 //_____________________________________________________________________________
455 void AliTRDclusterizer::ResetHelperIndexes(AliTRDSignalIndex *indexesIn)
458 // Reset the helper indexes
463 // carefull here - we assume that only row number may change - most probable
464 if (indexesIn->GetNrow() <= fIndexesOut->GetNrow())
465 fIndexesOut->ResetContent();
467 fIndexesOut->ResetContentConditional(indexesIn->GetNrow()
468 , indexesIn->GetNcol()
469 , indexesIn->GetNtime());
473 fIndexesOut = new AliTRDSignalIndex(indexesIn->GetNrow()
474 , indexesIn->GetNcol()
475 , indexesIn->GetNtime());
480 // carefull here - we assume that only row number may change - most probable
481 if (indexesIn->GetNrow() <= fIndexesMaxima->GetNrow())
483 fIndexesMaxima->ResetContent();
487 fIndexesMaxima->ResetContentConditional(indexesIn->GetNrow()
488 , indexesIn->GetNcol()
489 , indexesIn->GetNtime());
494 fIndexesMaxima = new AliTRDSignalIndex(indexesIn->GetNrow()
495 , indexesIn->GetNcol()
496 , indexesIn->GetNtime());
501 //_____________________________________________________________________________
502 Bool_t AliTRDclusterizer::ReadDigits()
505 // Reads the digits arrays from the input aliroot file
509 AliError("No run loader available");
513 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
514 if (!loader->TreeD()) {
515 loader->LoadDigits();
518 // Read in the digit arrays
519 return (fDigitsManager->ReadDigits(loader->TreeD()));
523 //_____________________________________________________________________________
524 Bool_t AliTRDclusterizer::ReadDigits(TTree *digitsTree)
527 // Reads the digits arrays from the input tree
530 // Read in the digit arrays
531 return (fDigitsManager->ReadDigits(digitsTree));
535 //_____________________________________________________________________________
536 Bool_t AliTRDclusterizer::ReadDigits(AliRawReader *rawReader)
539 // Reads the digits arrays from the ddl file
543 fDigitsManager = raw.Raw2Digits(rawReader);
549 //_____________________________________________________________________________
550 Bool_t AliTRDclusterizer::MakeClusters()
553 // Creates clusters from digits
556 // Propagate info from the digits manager
557 if (fAddLabels == kTRUE){
558 fAddLabels = fDigitsManager->UsesDictionaries();
561 Bool_t fReturn = kTRUE;
562 for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++){
564 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits*) fDigitsManager->GetDigits(i);
565 // This is to take care of switched off super modules
566 if (!digitsIn->HasData()) continue;
568 AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i);
569 if (indexes->IsAllocated() == kFALSE){
570 fDigitsManager->BuildIndexes(i);
574 if (indexes->HasEntry()){
576 for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++){
577 AliTRDdataArrayI *tracksIn = 0;
578 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(i,iDict);
582 fR = MakeClusters(i);
583 fReturn = fR && fReturn;
587 // if(IsWritingClusters()) WriteClusters(i);
591 // No compress just remove
592 fDigitsManager->RemoveDigits(i);
593 fDigitsManager->RemoveDictionaries(i);
594 fDigitsManager->ClearIndexes(i);
597 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
599 AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
605 //_____________________________________________________________________________
606 Bool_t AliTRDclusterizer::Raw2Clusters(AliRawReader *rawReader)
609 // Creates clusters from raw data
612 return Raw2ClustersChamber(rawReader);
616 //_____________________________________________________________________________
617 Bool_t AliTRDclusterizer::Raw2ClustersChamber(AliRawReader *rawReader)
620 // Creates clusters from raw data
623 // Create the digits manager
624 if (!fDigitsManager){
625 fDigitsManager = new AliTRDdigitsManager();
626 fDigitsManager->CreateArrays();
629 fDigitsManager->SetUseDictionaries(fAddLabels);
631 // tracklet container for raw tracklet writing
632 if (!fTrackletContainer && fReconstructor->IsWritingTracklets()) {
633 // maximum tracklets for one HC
634 const Int_t kTrackletChmb=256;
635 fTrackletContainer = new UInt_t *[2];
636 fTrackletContainer[0] = new UInt_t[kTrackletChmb];
637 fTrackletContainer[1] = new UInt_t[kTrackletChmb];
641 AliTRDrawStreamBase *pinput = AliTRDrawStreamBase::GetRawStream(rawReader);
642 AliTRDrawStreamBase &input = *pinput;
644 AliInfo(Form("Stream version: %s", input.IsA()->GetName()));
647 while ((det = input.NextChamber(fDigitsManager,fTrackletContainer)) >= 0){
648 Bool_t iclusterBranch = kFALSE;
649 if (fDigitsManager->GetIndexes(det)->HasEntry()){
650 iclusterBranch = MakeClusters(det);
653 fDigitsManager->RemoveDigits(det);
654 fDigitsManager->RemoveDictionaries(det);
655 fDigitsManager->ClearIndexes(det);
657 if (!fReconstructor->IsWritingTracklets()) continue;
658 if (*(fTrackletContainer[0]) > 0 || *(fTrackletContainer[1]) > 0) WriteTracklets(det);
660 if (fReconstructor->IsWritingTracklets()){
661 delete [] fTrackletContainer[0];
662 delete [] fTrackletContainer[1];
663 delete [] fTrackletContainer;
664 fTrackletContainer = NULL;
667 if(fReconstructor->IsWritingClusters()) WriteClusters(-1);
669 delete fDigitsManager;
670 fDigitsManager = NULL;
675 AliInfo(Form("Number of found clusters : %d", RecPoints()->GetEntriesFast()));
680 //_____________________________________________________________________________
681 UChar_t AliTRDclusterizer::GetStatus(Short_t &signal)
684 // Check if a pad is masked
689 if(signal>0 && TESTBIT(signal, 10)){
691 for(int ibit=0; ibit<4; ibit++){
692 if(TESTBIT(signal, 11+ibit)){
693 SETBIT(status, ibit);
694 CLRBIT(signal, 11+ibit);
701 //_____________________________________________________________________________
702 Bool_t AliTRDclusterizer::MakeClusters(Int_t det)
705 // Generates the cluster.
709 // digits should be expanded beforehand!
710 // digitsIn->Expand();
711 AliTRDdataArrayDigits *digitsIn = (AliTRDdataArrayDigits *) fDigitsManager->GetDigits(det);
713 // This is to take care of switched off super modules
714 if (!digitsIn->HasData())
719 AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det);
720 if (indexesIn->IsAllocated() == kFALSE)
722 AliError("Indexes do not exist!");
726 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
729 AliFatal("No AliTRDcalibDB instance available\n");
734 // There is no ADC threshold anymore, and simParam should not be used in clusterizer. KO
735 Float_t adcThreshold = 0;
737 if (!fReconstructor){
738 AliError("Reconstructor not set\n");
742 // Threshold value for the maximum
743 Float_t maxThresh = fReconstructor->GetRecoParam()->GetClusMaxThresh();
744 // Threshold value for the digit signal
745 Float_t sigThresh = fReconstructor->GetRecoParam()->GetClusSigThresh();
747 // Threshold value for the maximum ( cut noise)
748 Float_t minMaxCutSigma = fReconstructor->GetRecoParam()->GetMinMaxCutSigma();
749 // Threshold value for the sum pad ( cut noise)
750 Float_t minLeftRightCutSigma = fReconstructor->GetRecoParam()->GetMinLeftRightCutSigma();
752 // Iteration limit for unfolding procedure
753 const Float_t kEpsilon = 0.01;
754 const Int_t kNclus = 3;
755 const Int_t kNsig = 5;
758 Double_t ratioLeft = 1.0;
759 Double_t ratioRight = 1.0;
761 Double_t padSignal[kNsig];
762 Double_t clusterSignal[kNclus];
764 Int_t istack = indexesIn->GetStack();
765 Int_t ilayer = indexesIn->GetLayer();
766 Int_t isector = indexesIn->GetSM();
768 // Start clustering in the chamber
770 Int_t idet = AliTRDgeometry::GetDetector(ilayer,istack,isector);
773 AliError("Strange Detector number Missmatch!");
777 // TRD space point transformation
778 fTransform->SetDetector(det);
780 Int_t iGeoLayer = AliGeomManager::kTRD1 + ilayer;
781 Int_t iGeoModule = istack + AliTRDgeometry::Nstack() * isector;
782 UShort_t volid = AliGeomManager::LayerToVolUID(iGeoLayer,iGeoModule);
784 Int_t nColMax = digitsIn->GetNcol();
785 Int_t nRowMax = digitsIn->GetNrow();
786 Int_t nTimeTotal = digitsIn->GetNtime();
788 // Detector wise calibration object for the gain factors
789 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
790 // Calibration object with pad wise values for the gain factors
791 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
792 // Calibration value for chamber wise gain factor
793 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
795 // Detector wise calibration object for the noise
796 const AliTRDCalDet *calNoiseDet = calibration->GetNoiseDet();
797 // Calibration object with pad wise values for the noise
798 AliTRDCalROC *calNoiseROC = calibration->GetNoiseROC(idet);
799 // Calibration value for chamber wise noise
800 Float_t calNoiseDetValue = calNoiseDet->GetValue(idet);
804 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(nRowMax, nColMax, nTimeTotal);
805 AliTRDdataArrayS padStatus(nRowMax, nColMax, nTimeTotal);
807 ResetHelperIndexes(indexesIn);
809 // Apply the gain and the tail cancelation via digital filter
810 TailCancelation(digitsIn
817 ,calGainFactorDetValue);
824 UChar_t status[3]={0, 0, 0}, ipos = 0;
825 fIndexesOut->ResetCounters();
826 while (fIndexesOut->NextRCTbinIndex(row, col, time)) {
828 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time));
829 status[1] = digitsIn->GetPadStatus(row,col,time);
830 if(status[1]) SETBIT(ipos, AliTRDcluster::kMaskedCenter);
832 if(signalM < maxThresh) continue;
834 Float_t noiseMiddleThresh = minMaxCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
835 if (signalM < noiseMiddleThresh) continue;
837 if (col + 1 >= nColMax || col-1 < 0) continue;
839 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time));
840 status[0] = digitsIn->GetPadStatus(row,col+1,time);
841 if(status[0]) SETBIT(ipos, AliTRDcluster::kMaskedLeft);
843 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
844 status[2] = digitsIn->GetPadStatus(row,col-1,time);
845 if(status[2]) SETBIT(ipos, AliTRDcluster::kMaskedRight);
847 // reject candidates with more than 1 problematic pad
848 if(ipos == 3 || ipos > 4) continue;
850 if(!status[1]){ // good central pad
851 if(!ipos){ // all pads are OK
852 if ((signalL <= signalM) && (signalR < signalM)) {
853 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
854 Float_t noiseSumThresh = minLeftRightCutSigma*calNoiseDetValue*calNoiseROC->GetValue(col,row);
855 if((signalL+signalR+signalM) >= noiseSumThresh){
856 // Maximum found, mark the position by a negative signal
857 digitsOut->SetDataUnchecked(row,col,time,-signalM);
858 fIndexesMaxima->AddIndexTBin(row,col,time);
859 padStatus.SetDataUnchecked(row, col, time, ipos);
863 } else { // one of the neighbouring pads are bad
864 if(status[0] && signalR < signalM && signalR >= sigThresh){
865 digitsOut->SetDataUnchecked(row,col,time,-signalM);
866 digitsOut->SetDataUnchecked(row, col, time+1, 0.);
867 fIndexesMaxima->AddIndexTBin(row,col,time);
868 padStatus.SetDataUnchecked(row, col, time, ipos);
869 } else if(status[2] && signalL <= signalM && signalL >= sigThresh){
870 digitsOut->SetDataUnchecked(row,col,time,-signalM);
871 digitsOut->SetDataUnchecked(row, col, time-1, 0.);
872 fIndexesMaxima->AddIndexTBin(row,col,time);
873 padStatus.SetDataUnchecked(row, col, time, ipos);
876 } else { // wrong maximum pad
877 if ((signalL >= sigThresh) || (signalR >= sigThresh)) {
878 // Maximum found, mark the position by a negative signal
879 digitsOut->SetDataUnchecked(row,col,time,-maxThresh);
880 fIndexesMaxima->AddIndexTBin(row,col,time);
881 padStatus.SetDataUnchecked(row, col, time, ipos);
886 // The index to the first cluster of a given ROC
887 Int_t firstClusterROC = -1;
888 // The number of cluster in a given ROC
889 Int_t nClusterROC = 0;
891 // Now check the maxima and calculate the cluster position
892 fIndexesMaxima->ResetCounters();
893 while (fIndexesMaxima->NextRCTbinIndex(row, col, time)) {
896 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
898 for (iPad = 0; iPad < kNclus; iPad++) {
899 Int_t iPadCol = col - 1 + iPad;
900 clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
903 // Count the number of pads in the cluster
908 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) {
911 if (col-ii < 0) break;
915 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) {
918 if (col+ii+1 >= nColMax) break;
922 // Look for 5 pad cluster with minimum in the middle
923 Bool_t fivePadCluster = kFALSE;
924 if (col < (nColMax - 3)){
925 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
926 fivePadCluster = kTRUE;
928 if ((fivePadCluster) && (col < (nColMax - 5))) {
929 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) {
930 fivePadCluster = kFALSE;
933 if ((fivePadCluster) && (col > 1)) {
934 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) {
935 fivePadCluster = kFALSE;
941 // Modify the signal of the overlapping pad for the left part
942 // of the cluster which remains from a previous unfolding
944 clusterSignal[0] *= ratioLeft;
948 // Unfold the 5 pad cluster
950 for (iPad = 0; iPad < kNsig; iPad++) {
951 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
955 // Unfold the two maxima and set the signal on
956 // the overlapping pad to the ratio
957 ratioRight = Unfold(kEpsilon,ilayer,padSignal);
958 ratioLeft = 1.0 - ratioRight;
959 clusterSignal[2] *= ratioRight;
963 // The position of the cluster in COL direction relative to the center pad (pad units)
964 Double_t clusterPosCol = 0.0;
965 if (fReconstructor->GetRecoParam()->IsLUT()) {
966 // Calculate the position of the cluster by using the
967 // lookup table method
968 clusterPosCol = LUTposition(ilayer,clusterSignal[0]
973 // Calculate the position of the cluster by using the
974 // center of gravity method
975 for (Int_t i = 0; i < kNsig; i++) {
978 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
979 padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Left pad
980 padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Right pad
982 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
983 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
985 if ((col < nColMax - 3) &&
986 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])) {
987 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
989 clusterPosCol = GetCOG(padSignal);
992 // Store the amplitudes of the pads in the cluster for later analysis
993 // and check whether one of these pads is masked in the database
994 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
995 for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
997 (jPad >= nColMax-1)) {
1000 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
1003 // Transform the local cluster coordinates into calibrated
1004 // space point positions defined in the local tracking system.
1005 // Here the calibration for T0, Vdrift and ExB is applied as well.
1006 Double_t clusterXYZ[6];
1007 clusterXYZ[0] = clusterPosCol;
1008 clusterXYZ[1] = clusterSignal[2];
1009 clusterXYZ[2] = clusterSignal[1];
1010 clusterXYZ[3] = clusterSignal[0];
1011 clusterXYZ[4] = 0.0;
1012 clusterXYZ[5] = 0.0;
1013 Int_t clusterRCT[3];
1014 clusterRCT[0] = row;
1015 clusterRCT[1] = col;
1019 if (fTransform->Transform(clusterXYZ,clusterRCT,((UInt_t) time),out,0)) {
1021 // Add the cluster to the output array
1022 // The track indices will be stored later
1023 Float_t clusterPos[3];
1024 clusterPos[0] = clusterXYZ[0];
1025 clusterPos[1] = clusterXYZ[1];
1026 clusterPos[2] = clusterXYZ[2];
1027 Float_t clusterSig[2];
1028 clusterSig[0] = clusterXYZ[4];
1029 clusterSig[1] = clusterXYZ[5];
1030 Double_t clusterCharge = clusterXYZ[3];
1031 Char_t clusterTimeBin = ((Char_t) clusterRCT[2]);
1033 Int_t n = RecPoints()->GetEntriesFast();
1034 AliTRDcluster *cluster = new((*RecPoints())[n]) AliTRDcluster(idet
1039 ,((Char_t) nPadCount)
1047 cluster->SetInChamber(!out);
1049 UChar_t maskPosition = padStatus.GetDataUnchecked(row, col, time);
1051 cluster->SetPadMaskedPosition(maskPosition);
1052 if (maskPosition & AliTRDcluster::kMaskedLeft) {
1053 cluster->SetPadMaskedStatus(status[0]);
1055 else if (maskPosition & AliTRDcluster::kMaskedCenter) {
1056 cluster->SetPadMaskedStatus(status[1]);
1059 cluster->SetPadMaskedStatus(status[2]);
1063 // Temporarily store the row, column and time bin of the center pad
1064 // Used to later on assign the track indices
1065 cluster->SetLabel( row,0);
1066 cluster->SetLabel( col,1);
1067 cluster->SetLabel(time,2);
1069 // Store the index of the first cluster in the current ROC
1070 if (firstClusterROC < 0) {
1071 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
1074 // Count the number of cluster in the current ROC
1077 } // if: Transform ok ?
1079 } // if: Maximum found ?
1086 AddLabels(idet,firstClusterROC,nClusterROC);
1093 //_____________________________________________________________________________
1094 Bool_t AliTRDclusterizer::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC)
1097 // Add the track indices to the found clusters
1100 const Int_t kNclus = 3;
1101 const Int_t kNdict = AliTRDdigitsManager::kNDict;
1102 const Int_t kNtrack = kNdict * kNclus;
1104 Int_t iClusterROC = 0;
1111 // Temporary array to collect the track indices
1112 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
1114 // Loop through the dictionary arrays one-by-one
1115 // to keep memory consumption low
1116 AliTRDdataArrayI *tracksIn = 0;
1117 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
1119 // tracksIn should be expanded beforehand!
1120 tracksIn = (AliTRDdataArrayI *) fDigitsManager->GetDictionary(idet,iDict);
1122 // Loop though the clusters found in this ROC
1123 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1125 AliTRDcluster *cluster = (AliTRDcluster *)
1126 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1127 row = cluster->GetLabel(0);
1128 col = cluster->GetLabel(1);
1129 time = cluster->GetLabel(2);
1131 for (iPad = 0; iPad < kNclus; iPad++) {
1132 Int_t iPadCol = col - 1 + iPad;
1133 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
1134 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
1141 // Copy the track indices into the cluster
1142 // Loop though the clusters found in this ROC
1143 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
1145 AliTRDcluster *cluster = (AliTRDcluster *)
1146 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
1147 cluster->SetLabel(-9999,0);
1148 cluster->SetLabel(-9999,1);
1149 cluster->SetLabel(-9999,2);
1151 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
1155 delete [] idxTracks;
1161 //_____________________________________________________________________________
1162 Double_t AliTRDclusterizer::GetCOG(Double_t signal[5]) const
1166 // Used for clusters with more than 3 pads - where LUT not applicable
1169 Double_t sum = signal[0]
1176 // Go to 3 pad COG ????
1178 Double_t res = (0.0 * (-signal[0] + signal[4])
1179 + (-signal[1] + signal[3])) / sum;
1185 //_____________________________________________________________________________
1186 Double_t AliTRDclusterizer::Unfold(Double_t eps, Int_t layer, Double_t *padSignal)
1189 // Method to unfold neighbouring maxima.
1190 // The charge ratio on the overlapping pad is calculated
1191 // until there is no more change within the range given by eps.
1192 // The resulting ratio is then returned to the calling method.
1195 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
1197 AliError("No AliTRDcalibDB instance available\n");
1202 Int_t itStep = 0; // Count iteration steps
1204 Double_t ratio = 0.5; // Start value for ratio
1205 Double_t prevRatio = 0.0; // Store previous ratio
1207 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
1208 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
1209 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
1211 // Start the iteration
1212 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
1217 // Cluster position according to charge ratio
1218 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
1219 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
1220 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
1221 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
1223 // Set cluster charge ratio
1224 irc = calibration->PadResponse(1.0,maxLeft ,layer,newSignal);
1225 Double_t ampLeft = padSignal[1] / newSignal[1];
1226 irc = calibration->PadResponse(1.0,maxRight,layer,newSignal);
1227 Double_t ampRight = padSignal[3] / newSignal[1];
1229 // Apply pad response to parameters
1230 irc = calibration->PadResponse(ampLeft ,maxLeft ,layer,newLeftSignal );
1231 irc = calibration->PadResponse(ampRight,maxRight,layer,newRightSignal);
1233 // Calculate new overlapping ratio
1234 ratio = TMath::Min((Double_t) 1.0
1235 ,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]));
1243 //_____________________________________________________________________________
1244 void AliTRDclusterizer::TailCancelation(AliTRDdataArrayDigits *digitsIn
1245 , AliTRDdataArrayF *digitsOut
1246 , AliTRDSignalIndex *indexesIn
1247 , AliTRDSignalIndex *indexesOut
1249 , Float_t adcThreshold
1250 , AliTRDCalROC *calGainFactorROC
1251 , Float_t calGainFactorDetValue)
1254 // Applies the tail cancelation and gain factors:
1255 // Transform digitsIn to digitsOut
1262 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
1263 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
1264 indexesIn->ResetCounters();
1265 while (indexesIn->NextRCIndex(iRow, iCol))
1267 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
1268 Double_t gain = calGainFactorDetValue
1269 * calGainFactorROCValue;
1271 Bool_t corrupted = kFALSE;
1272 for (iTime = 0; iTime < nTimeTotal; iTime++)
1274 // Apply gain gain factor
1275 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
1276 if(digitsIn->GetPadStatus(iRow, iCol, iTime)) corrupted = kTRUE;
1277 inADC[iTime] /= gain;
1278 outADC[iTime] = inADC[iTime];
1282 // Apply the tail cancelation via the digital filter
1283 // (only for non-coorupted pads)
1284 if (fReconstructor->GetRecoParam() ->IsTailCancelation())
1286 DeConvExp(inADC,outADC,nTimeTotal,fReconstructor->GetRecoParam() ->GetTCnexp());
1290 indexesIn->ResetTbinCounter();
1291 while (indexesIn->NextTbinIndex(iTime))
1293 // Store the amplitude of the digit if above threshold
1294 if (outADC[iTime] > adcThreshold)
1296 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
1297 indexesOut->AddIndexTBin(iRow,iCol,iTime);
1301 } // while irow icol
1310 //_____________________________________________________________________________
1311 void AliTRDclusterizer::DeConvExp(Double_t *source, Double_t *target
1312 , Int_t n, Int_t nexp)
1315 // Tail cancellation by deconvolution for PASA v4 TRF
1319 Double_t coefficients[2];
1321 // Initialization (coefficient = alpha, rates = lambda)
1327 if (nexp == 1) { // 1 Exponentials
1333 if (nexp == 2) { // 2 Exponentials
1335 fReconstructor->GetTCParams(par);
1336 r1 = par[0];//1.156;
1337 r2 = par[1];//0.130;
1338 c1 = par[2];//0.114;
1339 c2 = par[3];//0.624;
1342 coefficients[0] = c1;
1343 coefficients[1] = c2;
1347 rates[0] = TMath::Exp(-dt/(r1));
1348 rates[1] = TMath::Exp(-dt/(r2));
1353 Double_t reminder[2];
1354 Double_t correction = 0.0;
1355 Double_t result = 0.0;
1357 // Attention: computation order is important
1358 for (k = 0; k < nexp; k++) {
1362 for (i = 0; i < n; i++) {
1364 result = (source[i] - correction); // No rescaling
1367 for (k = 0; k < nexp; k++) {
1368 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1372 for (k = 0; k < nexp; k++) {
1373 correction += reminder[k];
1380 //_____________________________________________________________________________
1381 void AliTRDclusterizer::ResetRecPoints()
1384 // Resets the list of rec points
1388 fRecPoints->Delete();
1393 //_____________________________________________________________________________
1394 TClonesArray *AliTRDclusterizer::RecPoints()
1397 // Returns the list of rec points
1401 if(!(fRecPoints = AliTRDReconstructor::GetClusters())){
1402 // determine number of clusters which has to be allocated
1403 Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
1404 if(fReconstructor->IsHLT()) nclusters /= AliTRDgeometry::kNsector;
1406 fRecPoints = new TClonesArray("AliTRDcluster", Int_t(nclusters));
1408 //SetClustersOwner(kTRUE);
1409 AliTRDReconstructor::SetClusters(0x0);
1415 //_____________________________________________________________________________
1416 void AliTRDclusterizer::FillLUT()
1422 const Int_t kNlut = 128;
1424 fLUTbin = AliTRDgeometry::kNlayer * kNlut;
1426 // The lookup table from Bogdan
1427 Float_t lut[AliTRDgeometry::kNlayer][kNlut] = {
1429 0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611,
1430 0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187,
1431 0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704,
1432 0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160,
1433 0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557,
1434 0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901,
1435 0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207,
1436 0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483,
1437 0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727,
1438 0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952,
1439 0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157,
1440 0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345,
1441 0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520,
1442 0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683,
1443 0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824,
1444 0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978
1447 0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642,
1448 0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241,
1449 0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770,
1450 0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229,
1451 0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627,
1452 0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968,
1453 0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266,
1454 0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536,
1455 0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772,
1456 0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991,
1457 0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187,
1458 0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369,
1459 0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538,
1460 0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694,
1461 0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832,
1462 0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1465 0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674,
1466 0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296,
1467 0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839,
1468 0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299,
1469 0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693,
1470 0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033,
1471 0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323,
1472 0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586,
1473 0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815,
1474 0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027,
1475 0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217,
1476 0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393,
1477 0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555,
1478 0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705,
1479 0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839,
1480 0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1483 0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708,
1484 0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356,
1485 0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910,
1486 0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371,
1487 0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759,
1488 0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095,
1489 0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378,
1490 0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633,
1491 0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857,
1492 0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061,
1493 0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244,
1494 0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414,
1495 0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571,
1496 0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715,
1497 0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846,
1498 0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1501 0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744,
1502 0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419,
1503 0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984,
1504 0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444,
1505 0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824,
1506 0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152,
1507 0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432,
1508 0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677,
1509 0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896,
1510 0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094,
1511 0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270,
1512 0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435,
1513 0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586,
1514 0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725,
1515 0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852,
1516 0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1519 0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792,
1520 0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505,
1521 0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081,
1522 0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541,
1523 0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917,
1524 0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235,
1525 0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511,
1526 0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748,
1527 0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962,
1528 0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152,
1529 0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324,
1530 0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483,
1531 0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629,
1532 0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759,
1533 0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859,
1534 0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1541 fLUT = new Double_t[fLUTbin];
1543 for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1544 for (Int_t ilut = 0; ilut < kNlut; ilut++ ) {
1545 fLUT[ilayer*kNlut+ilut] = lut[ilayer][ilut];
1551 //_____________________________________________________________________________
1552 Double_t AliTRDclusterizer::LUTposition(Int_t ilayer
1555 , Double_t ampR) const
1558 // Calculates the cluster position using the lookup table.
1559 // Method provided by Bogdan Vulpescu.
1562 const Int_t kNlut = 128;
1573 Double_t xMin[AliTRDgeometry::kNlayer] = { 0.006492, 0.006377, 0.006258
1574 , 0.006144, 0.006030, 0.005980 };
1575 Double_t xMax[AliTRDgeometry::kNlayer] = { 0.960351, 0.965870, 0.970445
1576 , 0.974352, 0.977667, 0.996101 };
1579 x = (ampL - ampR) / ampC;
1582 else if (ampL < ampR) {
1583 x = (ampR - ampL) / ampC;
1589 xmin = xMin[ilayer] + 0.000005;
1590 xmax = xMax[ilayer] - 0.000005;
1591 xwid = (xmax - xmin) / 127.0;
1596 else if (x > xmax) {
1597 pos = side * 0.5000;
1600 ix = (Int_t) ((x - xmin) / xwid);
1601 pos = side * fLUT[ilayer*kNlut+ix];