From ca21baaaf0120b0716661e8d9c037bd2151a8206 Mon Sep 17 00:00:00 2001 From: cblume Date: Mon, 30 Jul 2007 13:11:35 +0000 Subject: [PATCH] New fast cluster finder by Mateusz --- TRD/AliTRDRawStream.cxx | 358 +++++++++++- TRD/AliTRDRawStream.h | 4 + TRD/AliTRDReconstructor.cxx | 43 +- TRD/AliTRDReconstructor.h | 3 +- TRD/AliTRDSignalIndex.cxx | 449 +++++++++++++++ TRD/AliTRDSignalIndex.h | 159 ++++++ TRD/AliTRDclusterizerV2.cxx | 1040 +++++++++++++++++++++++++++++++++++ TRD/AliTRDclusterizerV2.h | 75 +++ TRD/AliTRDdigitsManager.cxx | 130 ++++- TRD/AliTRDdigitsManager.h | 16 +- TRD/AliTRDrawData.cxx | 33 +- TRD/TRDbaseLinkDef.h | 2 + TRD/TRDrecLinkDef.h | 1 + TRD/libTRDbase.pkg | 1 + TRD/libTRDrec.pkg | 1 + 15 files changed, 2265 insertions(+), 50 deletions(-) create mode 100644 TRD/AliTRDSignalIndex.cxx create mode 100644 TRD/AliTRDSignalIndex.h create mode 100644 TRD/AliTRDclusterizerV2.cxx create mode 100644 TRD/AliTRDclusterizerV2.h diff --git a/TRD/AliTRDRawStream.cxx b/TRD/AliTRDRawStream.cxx index 5b65761b290..e415481dbb0 100644 --- a/TRD/AliTRDRawStream.cxx +++ b/TRD/AliTRDRawStream.cxx @@ -34,6 +34,10 @@ #include "AliTRDgeometry.h" #include "AliTRDcalibDB.h" +#include "AliTRDdigitsManager.h" +#include "AliTRDdataArrayI.h" +#include "AliTRDSignalIndex.h" + ClassImp(AliTRDRawStream) //_____________________________________________________________________________ @@ -54,6 +58,7 @@ AliTRDRawStream::AliTRDRawStream() ,fROW(0) ,fCOL(0) ,fDET(0) + ,fLastDET(-1) ,fBCctr(0) ,fPTctr(0) ,fPTphase(0) @@ -132,6 +137,7 @@ AliTRDRawStream::AliTRDRawStream(AliRawReader *rawReader) ,fROW(0) ,fCOL(0) ,fDET(0) + ,fLastDET(-1) ,fBCctr(0) ,fPTctr(0) ,fPTphase(0) @@ -212,6 +218,7 @@ AliTRDRawStream::AliTRDRawStream(const AliTRDRawStream& stream) ,fROW(-1) ,fCOL(-1) ,fDET(0) + ,fLastDET(-1) ,fBCctr(-1) ,fPTctr(-1) ,fPTphase(-1) @@ -355,6 +362,7 @@ Int_t AliTRDRawStream::Init() fWordCtr = 0; fDET = 0; + fLastDET = -1; fRetVal = 0; fEqID = 0; fDataSize = 0; @@ -461,19 +469,27 @@ Bool_t AliTRDRawStream::Next() fNextStatus = fkNextHC; continue; } - if ( fADC > 1 && fADC < (Int_t)fGeo->ADCmax()-1 ) - { - // Write Digits - if ( fCOL >= 0 && fCOL < fColMax && fROW >= 0 && fROW < fRowMax ) - { // A real pad - fTB += 3; - return kTRUE; - } - } - else + if ( fRetVal == 1) { - fCOL = -1; + { // A real pad + fTB += 3; + return kTRUE; + } } + // following ifs have been moved to DEcodeDatawordV1V2 +// if ( fADC > 1 && fADC < (Int_t)fGeo->ADCmax()-1 ) +// { +// // Write Digits +// if ( fCOL >= 0 && fCOL < fColMax && fROW >= 0 && fROW < fRowMax ) +// { // A real pad +// fTB += 3; +// return kTRUE; +// } +// } +// else +// { +// fCOL = -1; +// } }// fkNextData continue; @@ -537,6 +553,7 @@ Bool_t AliTRDRawStream::Next() if ( (*fDataWord & 0x00000003) == 1 ) { // HC header DecodeHCheader(fTimeBinsCalib); // This is the new header! + fLastDET = fDET; fDET = fGeo->GetDetector(fLAYER, fSTACK, fSM); fRowMax = fGeo->GetRowMax(fLAYER,fSTACK,fSM); fColMax = fGeo->GetColMax(fROC); @@ -595,6 +612,276 @@ Bool_t AliTRDRawStream::Next() return kFALSE; } +//____________________________________________________________________________ +Int_t AliTRDRawStream::NextChamber(AliTRDdigitsManager *man) +{ + // + // Updates the next data word pointer + // + + AliTRDdataArrayI *digits = 0; + AliTRDdataArrayI *track0 = 0; + AliTRDdataArrayI *track1 = 0; + AliTRDdataArrayI *track2 = 0; + AliTRDSignalIndex *indexes = 0; + + if (fNextStatus == fkStart) + { + Init(); + } + + while (fNextStatus != fkStop) + { // !fkStop + NextData(); + if (fNextStatus == fkNextMCM || fNextStatus == fkNextData) + //while (fNextStatus == fkNextMCM || fNextStatus == fkNextData) + { + fHCdataCtr += 4; + + if( ((*fDataWord & 0x80000000) == 0x0) && ((*fDataWord & 0x0000000f) == 0xC) ) + { // MCM Header + DecodeMCMheader(); + if ( fMCM < 0 || fMCM > 15 || fROB < 0 || fROB > 7 ) + { + AliWarning("Wrong fMCM or fROB. Skip this data"); + fRawReader->AddMajorErrorLog(kWrongMCMorROB,Form("MCM=%d, ROB=%d",fMCM,fROB)); + fNextStatus = fkNextHC; + continue; + } + fTbSwitch = 3; // For first adc channel we expect: (*fDataWord & 3) = 3 + fTbSwitchCtr = 0; // + fADC = fTB = 0; // Reset Counter + fNextStatus = fkNextData; + +// NextData(); // if while loop! + continue; // if if + } + + if ( *fDataWord == kEndofrawdatamarker ) + { // End of half-chamber data, finished + fGTUctr1 = -1; + fNextStatus = fkNextHC; + // full chamber processed ? + if (fChamberDone[fDET] == 2) + { + return fDET; + } + else + { +// break; // if while loop + continue; // if if + } + } + + if (fNextStatus == fkNextData ) + { // MCM header is set, ADC data is valid. + + // Found some data. Decode it now: + fRetVal = DecodeDataWord(); + if ( fRetVal == 0 ) continue; + if ( fRetVal == -1 ) + { + fNextStatus = fkNextHC; + +// NextData(); // if while loop! +// break; //if while loop! + continue;// if if + } + + if ( fRetVal == 1) + { + { // A real pad + // here fill the data arrays + //return kTRUE; + for (Int_t it = 0; it < 3; it++) + { + if ( GetTimeBin() + it < GetNumberOfTimeBins() ) + { + if (GetSignals()[it] > 0) + { + digits->SetDataUnchecked(fROW, fCOL, fTB + it, fSig[it]); + indexes->AddIndexTBin(fROW, fCOL, fTB + it); + track0->SetDataUnchecked(fROW, fCOL, fTB + it, 0); + track1->SetDataUnchecked(fROW, fCOL, fTB + it, 0); + track2->SetDataUnchecked(fROW, fCOL, fTB + it, 0); + } + } // check the tbins range + } // for each tbin of current 3 + fTB += 3; + }// real pad + } // if fRetVal == 1 + + // following ifs have been moved to DEcodeDatawordV1V2 +// if ( fADC > 1 && fADC < (Int_t)fGeo->ADCmax()-1 ) +// { +// // Write Digits +// if ( fCOL >= 0 && fCOL < fColMax && fROW >= 0 && fROW < fRowMax ) +// { // A real pad +// fTB += 3; +// return kTRUE; +// } +// } +// else +// { +// fCOL = -1; +// } + }// fkNextData + +// NextData(); // if while loop! + continue; //if if + } //next mcm + + if ( fNextStatus == fkNextHC ) + { + // + // 1) Find end_of_tracklet_marker + // + // GTU Link Mask? + if ( (*fDataWord & 0xfffff000) == 0xe0000000 ) + { + DecodeGTUlinkMask(); + continue; + } + + // endoftrackletmarker? + if ( *fDataWord == kEndoftrackletmarker ) + { + AliDebug(3, "end-of-tracklet-marker found"); + fNextStatus = fkSeekNonEoTracklet; + continue; + } + else + { + // Tracklets found + AliDebug(3, "Tracklet found"); + DecodeTracklet(); + continue; + } + } //if next HC + + if (fNextStatus == fkSeekNonEoTracklet) + { + // + // 2) Look for non-end_of_tracklet_marker + // + //printf("Word %d: 0x%08x\n", fWordCtr, *fDataWord); + + if ( *fDataWord != kEndoftrackletmarker ) + { + fNextStatus = fkDecodeHC; + AliDebug(3, "NON end-of-tracklet-marker found"); + //// no do not continue - this should be the hcheader + } + else + { + //just go on and find the non-end_of_tracklet_marker + continue; + } + } + + if ( fNextStatus == fkDecodeHC ) + { + AliDebug(3, "Decode HC"); + + // + // 3) This Word must be Half Chamber Header + // + if ( (*fDataWord & 0x00000003) == 1 ) + { // HC header + DecodeHCheader(fTimeBinsCalib); // This is the new header! + fDET = fGeo->GetDetector(fLAYER, fSTACK, fSM); + fRowMax = fGeo->GetRowMax(fLAYER,fSTACK,fSM); + fColMax = fGeo->GetColMax(fROC); + + if (fLastDET != fDET) + { + AliDebug(4, "New DET!"); + // allocate stuff for the new det + digits = man->GetDigits(fDET); + track0 = man->GetDictionary(fDET,0); + track1 = man->GetDictionary(fDET,1); + track2 = man->GetDictionary(fDET,2); + + // Allocate memory space for the digits buffer + if (digits->GetNtime() == 0) + { + AliDebug(4, "Allocating digits"); + //AliDebug(5, Form("Alloc digits for det %d", det)); + digits->Allocate(fRowMax, fColMax, fTBins); + track0->Allocate(fRowMax, fColMax, fTBins); + track1->Allocate(fRowMax, fColMax, fTBins); + track2->Allocate(fRowMax, fColMax, fTBins); + } + + indexes = man->GetIndexes(fDET); + indexes->SetSM(fSM); + indexes->SetStack(fSTACK); + indexes->SetLayer(fLAYER); + indexes->SetDetNumber(fDET); + + if (indexes->IsAllocated() == kFALSE) + { + AliDebug(4, "Allocating indexes"); + indexes->Allocate(fRowMax, fColMax, fTBins); + } + fLastDET = fDET; + } + + fMCMHctr2 = 0; + fHCdataCtr = 0; + + fChamberDone[fDET]++; + fNextStatus = fkNextMCM; + AliDebug(3, "Decode HC OK"); + continue; + } //HC header + else + { + AliDebug(3, "Decode HC NOT OK"); + fNextStatus = fkNextSM; + continue; + } + } // if decode HC + + if (fNextStatus == fkNextSM) + { + + fDET = 0; + fRetVal = 0; + fEqID = 0; + fDataSize = 0; + fSizeOK = kFALSE; + + // After reading the first word check for size of this data and get Eq. ID + if ( fWordCtr == 1 ) + { + fDataSize = fRawReader->GetDataSize()/4; // Size of this payload in 32bit words + fEqID = fRawReader->GetEquipmentId(); // Get Equipment ID + if ( fDataSize > 0 ) fSizeOK = kTRUE; + } + + // GTU Link Mask? + if ( (*fDataWord & 0xfffff000) == 0xe0000000 ) + { + DecodeGTUlinkMask(); + fNextStatus = fkNextHC; + continue; + } + else + { + AliWarning(Form("Equipment %d: First data word is not GTU Link Mask!", fEqID)); + fRawReader->AddMajorErrorLog(kGTULinkMaskMissing,Form("Equipment %d",fEqID)); + fNextStatus = fkStop; + } + }// if nextSM + + } // not fkStop + + AliDebug(1, Form("That's all folks! %d", fSM)); + //return kFALSE; + return -1; +} + //============================================================================ // Decoding functions //============================================================================ @@ -895,8 +1182,16 @@ Int_t AliTRDRawStream::DecodeDataWordV1V2() // return -1 means break data loop // +// // check the content first! - something wrong with that... +// // Decode 32 bit data words with information from 3 time bins and copy the data +// fSig[0] = (*fDataWord & 0x00000ffc) >> 2; +// fSig[1] = (*fDataWord & 0x003ff000) >> 12; +// fSig[2] = (*fDataWord & 0xffc00000) >> 22; +// if (fSig[0] <= 0 && fSig[1] <= 0 && fSig[2] <= 0) +// return 0; + if ( (*fDataWord & 0x00000003) != 0x2 && (*fDataWord & 0x00000003) != 0x3) { - AliWarning(Form("Data %08x : Data Word ends neither with b11 nor b10", (Int_t)*fDataWord)); + //AliWarning(Form("Data %08x : Data Word ends neither with b11 nor b10", (Int_t)*fDataWord)); fRawReader->AddMinorErrorLog(kDataMaskError,Form("Data %08x", (Int_t)*fDataWord)); return -1; } @@ -912,14 +1207,14 @@ Int_t AliTRDRawStream::DecodeDataWordV1V2() // We have only timeTotal time bins if ( fTbSwitchCtr > fTimeWords ) { - AliWarning(Form("Data is strange. Already found %d words for this ADC channel", (Int_t)fTbSwitchCtr)); + //AliWarning(Form("Data is strange. Already found %d words for this ADC channel", (Int_t)fTbSwitchCtr)); fRawReader->AddMinorErrorLog(kADCNumberOverflow,Form("%d words", (Int_t)fTbSwitchCtr)); return 0; } // We have only 21 ADC channels. if ( fADC > (Int_t)fGeo->ADCmax()-1 ) { - AliWarning(Form("Data %08x : Data is strange. fADC is already %d", (Int_t)*fDataWord, (Int_t)fADC)); + //AliWarning(Form("Data %08x : Data is strange. fADC is already %d", (Int_t)*fDataWord, (Int_t)fADC)); fRawReader->AddMinorErrorLog(kADCChannelOverflow,Form("Data %08x : fADC=%d", (Int_t)*fDataWord, (Int_t)fADC)); return 0; } @@ -932,18 +1227,27 @@ Int_t AliTRDRawStream::DecodeDataWordV1V2() fCOL = fGeo->GetPadColFromADC(fROB, fMCM, fADC); // We have only 144 Pad Columns - if ( fCOL > fColMax-1 || fCOL < 0 ) { - AliWarning(Form("SM%d L%dS%d: Wrong Pad column (%d) fROB=%d, fSIDE=%d, fMCM=%02d", fSM, - fLAYER, fSTACK, fCOL, fROB, fSIDE, fMCM )); - fRawReader->AddMajorErrorLog(kWrongPadcolumn,Form("SM%d L%dS%d: column (%d) fROB=%d, fSIDE=%d, fMCM=%02d", fSM, - fLAYER, fSTACK, fCOL, fROB, fSIDE, fMCM )); - } - - // Decode 32 bit data words with information from 3 time bins and copy the data - fSig[0] = (*fDataWord & 0x00000ffc) >> 2; - fSig[1] = (*fDataWord & 0x003ff000) >> 12; - fSig[2] = (*fDataWord & 0xffc00000) >> 22; - + //if ( fCOL > fColMax-1 || fCOL < 0 ) { + if ( fCOL >= 0 && fCOL < fColMax && fROW >= 0 && fROW < fRowMax ) + { + // Decode 32 bit data words with information from 3 time bins and copy the data + fSig[0] = (*fDataWord & 0x00000ffc) >> 2; + fSig[1] = (*fDataWord & 0x003ff000) >> 12; + fSig[2] = (*fDataWord & 0xffc00000) >> 22; + + if (fSig[0] > 0 || fSig[1] > 0 || fSig[2] > 0) + return 1; + else + return 0; + } + else + { +// AliWarning(Form("SM%d L%dS%d: Wrong Pad column (%d) fROB=%d, fSIDE=%d, fMCM=%02d", fSM, +// fLAYER, fSTACK, fCOL, fROB, fSIDE, fMCM )); + fRawReader->AddMajorErrorLog(kWrongPadcolumn,Form("SM%d L%dS%d: column (%d) fROB=%d, fSIDE=%d, fMCM=%02d", fSM, + fLAYER, fSTACK, fCOL, fROB, fSIDE, fMCM )); + return 0; + } // Print data to screen: // Do NOT switch on for default production, it is VERY slow // AliDebug(5, Form("SM%d L%dS%d: ROB%d MCM=%d ADC=%d (ROW=%d COL=%d): Data %04d %04d %04d\n", @@ -953,7 +1257,7 @@ Int_t AliTRDRawStream::DecodeDataWordV1V2() else { fCOL = -1; - + return 0; } return 1; diff --git a/TRD/AliTRDRawStream.h b/TRD/AliTRDRawStream.h index 06eec72a33a..ed254f42382 100644 --- a/TRD/AliTRDRawStream.h +++ b/TRD/AliTRDRawStream.h @@ -15,6 +15,7 @@ class AliTRDgeometry; class AliRawReader; +class AliTRDdigitsManager; // Some constants: const UInt_t kEndoftrackletmarker = 0xAAAAAAAA; /*This marks the end of tracklet data words*/ @@ -30,6 +31,7 @@ class AliTRDRawStream: public TObject { virtual ~AliTRDRawStream(); virtual Bool_t Next(); // Read the next data + virtual Int_t NextChamber(AliTRDdigitsManager *man); // read next chamber data virtual Int_t Init(); // Init for the fRawVersion > 1 enum { kDDLOffset = 0x400 }; // Offset for DDL numbers @@ -69,6 +71,7 @@ class AliTRDRawStream: public TObject { Int_t GetRow() const { return fROW;} Int_t GetCol() const { return fCOL;} // Detector Pad coordinates Int_t GetDet() const { return fDET;} // helper + Int_t GetLastDet() const { return fLastDET;} // helper Int_t GetMaxRow() const { return fRowMax;} Int_t GetMaxCol() const { return fColMax;} Int_t GetNumberOfTimeBins() const {return fTBins;} @@ -90,6 +93,7 @@ class AliTRDRawStream: public TObject { Int_t fROW; // Detector Row coordinates Int_t fCOL; // Detector Pad coordinates Int_t fDET; // Current detector - version > 1 + Int_t fLastDET; // Previous detector - version > 1 Int_t fBCctr; // Counters from HC header (>=V2) Int_t fPTctr; // Counters from HC header (>=V2) diff --git a/TRD/AliTRDReconstructor.cxx b/TRD/AliTRDReconstructor.cxx index daaf7a06d92..58ea95564b6 100644 --- a/TRD/AliTRDReconstructor.cxx +++ b/TRD/AliTRDReconstructor.cxx @@ -31,6 +31,8 @@ #include "AliTRDReconstructor.h" #include "AliTRDclusterizerV1.h" +#include "AliTRDclusterizerV2.h" +//#include "AliTRDclusterizerV2xMP.h" #include "AliTRDtracker.h" #include "AliTRDpidESD.h" #include "AliTRDgtuTrack.h" @@ -50,6 +52,8 @@ void AliTRDReconstructor::ConvertDigits(AliRawReader *rawReader // Convert raw data digits into digit objects in a root tree // + AliInfo("Convert raw data digits into digit objects [RawReader -> Digit TTree]"); + AliTRDrawData rawData; rawReader->Reset(); rawReader->Select("TRD"); @@ -67,7 +71,7 @@ void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader // Reconstruct clusters // - AliInfo("Reconstruct TRD clusters from RAW data"); + AliInfo("Reconstruct TRD clusters from RAW data [RunLoader, RawReader]"); AliLoader *loader = runLoader->GetLoader("TRDLoader"); loader->LoadRecPoints("recreate"); @@ -79,12 +83,22 @@ void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader rawReader->Select("TRD"); for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) { + if (!rawReader->NextEvent()) break; - AliTRDclusterizerV1 clusterer("clusterer","TRD clusterizer"); + + // Old (slow) cluster finder + //AliTRDclusterizerV1 clusterer("clusterer","TRD clusterizer"); + //clusterer.Open(runLoader->GetFileName(),iEvent); + //clusterer.ReadDigits(rawReader); + //clusterer.MakeClusters(); + + // New (fast) cluster finder + AliTRDclusterizerV2 clusterer("clusterer","TRD clusterizer"); clusterer.Open(runLoader->GetFileName(),iEvent); - clusterer.ReadDigits(rawReader); - clusterer.MakeClusters(); + clusterer.Raw2ClustersChamber(rawReader); + clusterer.WriteClusters(-1); + } loader->UnloadRecPoints(); @@ -99,15 +113,22 @@ void AliTRDReconstructor::Reconstruct(AliRawReader *rawReader // Reconstruct clusters // - AliInfo("Reconstruct TRD clusters from RAW data"); + AliInfo("Reconstruct TRD clusters from RAW data [RawReader -> Cluster TTree]"); rawReader->Reset(); rawReader->Select("TRD"); - AliTRDclusterizerV1 clusterer("clusterer","TRD clusterizer"); + // Old (slow) cluster finder + //AliTRDclusterizerV1 clusterer("clusterer","TRD clusterizer"); + //clusterer.OpenOutput(clusterTree); + //clusterer.ReadDigits(rawReader); + //clusterer.MakeClusters(); + + // New (fast) cluster finder + AliTRDclusterizerV2 clusterer("clusterer","TRD clusterizer"); clusterer.OpenOutput(clusterTree); - clusterer.ReadDigits(rawReader); - clusterer.MakeClusters(); + clusterer.SetAddLabels(kFALSE); + clusterer.Raw2ClustersChamber(rawReader); } @@ -118,8 +139,10 @@ void AliTRDReconstructor::Reconstruct(TTree *digitsTree // // Reconstruct clusters // + AliInfo("Reconstruct TRD clusters from Digits [Digit TTree -> Cluster TTree]"); AliTRDclusterizerV1 clusterer("clusterer","TRD clusterizer"); + //AliTRDclusterizerV2 clusterer("clusterer","TRD clusterizer"); clusterer.OpenOutput(clusterTree); clusterer.ReadDigits(digitsTree); clusterer.MakeClusters(); @@ -133,6 +156,7 @@ void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader) const // Reconstruct clusters // + AliInfo("Reconstruct TRD clusters [AliRunLoader]"); AliLoader *loader = runLoader->GetLoader("TRDLoader"); loader->LoadRecPoints("recreate"); @@ -141,6 +165,7 @@ void AliTRDReconstructor::Reconstruct(AliRunLoader *runLoader) const for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) { AliTRDclusterizerV1 clusterer("clusterer","TRD clusterizer"); + //AliTRDclusterizerV2 clusterer("clusterer","TRD clusterizer"); clusterer.Open(runLoader->GetFileName(),iEvent); clusterer.ReadDigits(); clusterer.MakeClusters(); @@ -195,7 +220,7 @@ void AliTRDReconstructor::FillESD(AliRawReader* /*rawReader*/ //_____________________________________________________________________________ void AliTRDReconstructor::FillESD(TTree* /*digitsTree*/ , TTree* /*clusterTree*/ - , AliESDEvent* /*esd*/) const + , AliESDEvent* /*esd*/) const { // // Make PID diff --git a/TRD/AliTRDReconstructor.h b/TRD/AliTRDReconstructor.h index fefaa4eb272..624746f69a5 100644 --- a/TRD/AliTRDReconstructor.h +++ b/TRD/AliTRDReconstructor.h @@ -22,7 +22,8 @@ class AliTRDReconstructor: public AliReconstructor { AliTRDReconstructor():AliReconstructor() { }; virtual ~AliTRDReconstructor() { }; - virtual Bool_t HasDigitConversion() const { return kTRUE; }; + //virtual Bool_t HasDigitConversion() const { return kTRUE; }; + virtual Bool_t HasDigitConversion() const { return kFALSE; }; virtual void ConvertDigits(AliRawReader *rawReader, TTree *digitsTree) const; virtual Bool_t HasLocalReconstruction() const { return kTRUE; }; diff --git a/TRD/AliTRDSignalIndex.cxx b/TRD/AliTRDSignalIndex.cxx new file mode 100644 index 00000000000..3dcdd0a02c5 --- /dev/null +++ b/TRD/AliTRDSignalIndex.cxx @@ -0,0 +1,449 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// General container for data from TRD detector segments // +// Adapted from AliDigits, origin M.Ivanov // +// // +// Author: M. Ploskon (ploskon@ikf.uni-frankfurt.de) // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include "AliTRDSignalIndex.h" + +#include "TArrayI.h" +#include "AliLog.h" + +ClassImp(AliTRDSignalIndex) + +//_____________________________________________________________________________ +AliTRDSignalIndex::AliTRDSignalIndex() + : TObject() + , fDet(-1) + , fLayer(-1) + , fStack(-1) + , fSM(-1) + , fIndex(NULL) + , fPositionRow(0) + , fPositionCol(0) + , fPositionTbin(0) + , fLastRow(0) + , fLastCol(0) + , fLastTbin(0) + , fNrows(0) + , fNcols(0) + , fNtbins(0) + , fMaxLimit(0) + , fResetCounters(kTRUE) + , fHasEntry(kFALSE) +{ + // + // default contructor + // + ResetCounters(); +} + +//_____________________________________________________________________________ +AliTRDSignalIndex::AliTRDSignalIndex(Int_t nrow, Int_t ncol,Int_t ntime) + : TObject() + , fDet(-1) + , fLayer(-1) + , fStack(-1) + , fSM(-1) + , fIndex(NULL) + , fPositionRow(0) + , fPositionCol(0) + , fPositionTbin(0) + , fLastRow(0) + , fLastCol(0) + , fLastTbin(0) + , fNrows(0) + , fNcols(0) + , fNtbins(0) + , fMaxLimit(0) + , fResetCounters(kTRUE) + , fHasEntry(kFALSE) +{ + // not the default contructor... hmmm... + Allocate(nrow, ncol, ntime); +} + +//_____________________________________________________________________________ +AliTRDSignalIndex::~AliTRDSignalIndex() +{ + // + // Destructor + // + + delete fIndex; + fIndex = NULL; +} + +//_____________________________________________________________________________ +AliTRDSignalIndex::AliTRDSignalIndex(const AliTRDSignalIndex &a) + : TObject(a) + , fDet(a.fDet) + , fLayer(a.fLayer) + , fStack(a.fStack) + , fSM(a.fSM) + , fIndex(a.fIndex) + , fPositionRow(a.fPositionRow) + , fPositionCol(a.fPositionCol) + , fPositionTbin(a.fPositionTbin) + , fLastRow(a.fLastRow) + , fLastCol(a.fLastCol) + , fLastTbin(a.fLastTbin) + , fNrows(a.fNrows) + , fNcols(a.fNcols) + , fNtbins(a.fNtbins) + , fMaxLimit(a.fMaxLimit) + , fResetCounters(a.fResetCounters) + , fHasEntry(a.fHasEntry) +{ + // + // Copy constructor + // +} + +//_____________________________________________________________________________ +void AliTRDSignalIndex::Copy(TObject &a) const +{ + // + // Copy function + // + + ((AliTRDSignalIndex &)a).fDet = fDet; + ((AliTRDSignalIndex &)a).fLayer = fLayer; + ((AliTRDSignalIndex &)a).fStack = fStack; + ((AliTRDSignalIndex &)a).fSM = fSM; + ((AliTRDSignalIndex &)a).fIndex = fIndex; + ((AliTRDSignalIndex &)a).fPositionRow = fPositionRow; + ((AliTRDSignalIndex &)a).fPositionTbin = fPositionTbin; + ((AliTRDSignalIndex &)a).fLastRow = fLastRow; + ((AliTRDSignalIndex &)a).fLastCol = fLastCol; + ((AliTRDSignalIndex &)a).fLastTbin = fLastTbin; + ((AliTRDSignalIndex &)a).fNrows = fNrows; + ((AliTRDSignalIndex &)a).fNcols = fNcols; + ((AliTRDSignalIndex &)a).fNtbins = fNtbins; + ((AliTRDSignalIndex &)a).fMaxLimit = fMaxLimit; + ((AliTRDSignalIndex &)a).fResetCounters = fResetCounters; + ((AliTRDSignalIndex &)a).fHasEntry = fHasEntry; +} + +//_____________________________________________________________________________ +AliTRDSignalIndex& AliTRDSignalIndex::operator = (const AliTRDSignalIndex& a) +{ + // + // Assignment operator + // + + if (this != &a) ((AliTRDSignalIndex &) a).Copy(*this); + return *this; + +} +//_____________________________________________________________________________ +void AliTRDSignalIndex::Allocate(Int_t nrow, Int_t ncol,Int_t ntime) +{ + // + // create the arrays + // + + fNrows = nrow; + fNcols = ncol; + fNtbins = ntime; + + fMaxLimit = nrow * ncol * ntime + nrow * ncol * 2; + if (fIndex) + { + delete fIndex; + fIndex = NULL; + } + + fIndex = new TArrayI(fMaxLimit); + //fIndex->Set(fMaxLimit); + fIndex->Reset(-1); + ResetCounters(); + fHasEntry = kFALSE; +} + +//_____________________________________________________________________________ +void AliTRDSignalIndex::Reset() +{ + // + // Reset the array but keep the size - realloc + // + + fDet = -1; + fLayer = -1; + fStack = -1; + fSM = -1; + + // all will be lost + Allocate(fNrows, fNcols, fNtbins); +} + +//_____________________________________________________________________________ +void AliTRDSignalIndex::ClearAll() +{ + // + // Reset the values - clear all! + // + + fDet = -1; + fLayer = -1; + fStack = -1; + fSM = -1; + + fNrows = -1; + fNcols = -1; + fNtbins = -1; + + if (fIndex) + { + delete fIndex; + fIndex = NULL; + } + fIndex = new TArrayI(); + ResetCounters(); + fHasEntry = kFALSE; +} + +//_____________________________________________________________________________ +void AliTRDSignalIndex::AddIndexTBin(Int_t row, Int_t col, Int_t tbin) +{ + // store the index row-column-tbin as an interesting one + // the RC index is updated to!!! + // this is to be used in the TRD clusterizer! + + if (fPositionCol + fNtbins >= fMaxLimit) + { + AliError(Form("Out-of-limits fPositionCol + fNtbins %d. Limit is: %d", fPositionCol + fNtbins, fMaxLimit)); + return; + } + + if (row != fLastRow || col != fLastCol) + { + // new RC combination + if (fResetCounters == kTRUE) + { + fPositionRow = 0; + fPositionCol = 1; + + fResetCounters = kFALSE; + } + else + { + fPositionRow += fNtbins + 2; + fPositionCol += fNtbins + 2; + } + + fPositionTbin = 1; + + (*fIndex)[fPositionRow] = row; + (*fIndex)[fPositionCol] = col; + (*fIndex)[fPositionCol + fPositionTbin] = tbin; + + ++fPositionTbin; + //AliDebug(3, Form("fNRCindexed=%d", fNRCindexed)); + } + else + { + // same RCT combination ? + // if (fLastTbin == tbin) + // { + // AliWarning(Form("Same RCT? %d %d %d", row, col, tbin)); + // } + + (*fIndex)[fPositionCol + fPositionTbin] = tbin; + ++fPositionTbin; + } + + fLastRow = row; + fLastCol = col; + fLastTbin = tbin; + + fHasEntry = kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliTRDSignalIndex::NextRCIndex(Int_t &row, Int_t &col) +{ + // return the position (index in the data array) of the next available pad + + if (fPositionCol + fNtbins >= fMaxLimit) + { + //AliDebug(8, "Out of index range"); + return kFALSE; + } + + if (fResetCounters == kTRUE) + { + fPositionRow = 0; + fPositionCol = 1; + + fResetCounters = kFALSE; + AliDebug(2, "Reset counters"); + } + else + { + fPositionRow += fNtbins + 2; + fPositionCol += fNtbins + 2; + } + + //AliDebug(8, Form("Next RC %d", fPositionRow / (fNtbins + 2))); + + fPositionTbin = 1; + + row = (*fIndex)[fPositionRow]; + col = (*fIndex)[fPositionCol]; + + if (row > -1 && col > -1) + return kTRUE; + + return kFALSE; +} + +//_____________________________________________________________________________ +Bool_t AliTRDSignalIndex::NextRCTbinIndex(Int_t &row, Int_t &col, Int_t &tbin) +{ + // return the position (index in the data array) of the next available tbin + // within the current pad + +// if (fNRCcounter >= fNRCindexed) +// return kFALSE; + + if (fPositionCol + fNtbins >= fMaxLimit) + { + return kFALSE; + } + + if (NextTbinIndex(tbin)) + { + row = (*fIndex)[fPositionRow]; + col = (*fIndex)[fPositionCol]; + fResetCounters = kFALSE; + return kTRUE; + } + else + { + if (NextRCIndex(row, col)) + { + //return NextTbinIndex(tbin); + return NextRCTbinIndex(row, col, tbin); + } + } + + return kFALSE; +} + +//_____________________________________________________________________________ +Bool_t AliTRDSignalIndex::NextTbinIndex(Int_t &tbin) +{ + // return the position (index in the data array) of the next available tbin + // within the current pad + +// if (fNRCcounter >= fNRCindexed) +// return kFALSE; + + if (fPositionCol + fNtbins >= fMaxLimit || fPositionTbin > fNtbins) + { + return kFALSE; + } + + tbin = (*fIndex)[fPositionCol + fPositionTbin]; + + if (tbin > -1) + { + ++fPositionTbin; + return kTRUE; + } + + return kFALSE; +} + +// void AliTRDSignalIndex::Dump() +// { +// AliInfo("R C T..."); +// Int_t i = 0; +// Int_t rcok = 0; +// while ( i < fMaxLimit ) +// { +// if (i % (fNtbins + 2) == 0) +// { +// if ((*(*fIndex))[i] > -1) +// { +// rcok = 1; +// printf("\n RC : "); +// } +// else +// rcok = 0; +// } +// if (rcok) +// if ((*(*fIndex))[i] > -1) +// printf("[%d] = %d \t", i, (*(*fIndex))[i]); +// i++; +// } +// } + +// //_____________________________________________________________________________ +// void AliTRDSignalIndex::Dump() +// { +// // +// // Dump the data +// // + +// AliInfo("R C T..."); +// Int_t i = 0; +// Int_t rcok = 0; +// while ( i < fMaxLimit ) +// { +// if (i % (fNtbins + 2) == 0) +// { +// if ((*fIndex)[i] > -1) +// { +// rcok = 1; +// printf("\n RC : "); +// } +// else +// rcok = 0; +// } +// if (rcok) +// if ((*fIndex)[i] > -1) +// printf("[%d] = %d \t", i, (*fIndex)[i]); +// i++; +// } +// } + +// //_____________________________________________________________________________ +// void AliTRDSignalIndex::Dump2() +// { +// // +// // Dump the data +// // + +// AliInfo("R C T..."); +// Int_t ir, ic, it; +// ResetCounters(); +// while (NextRCIndex(ir, ic)) +// { +// printf("\nR %d C %d t : ", ir, ic); +// while (NextTbinIndex(it)) +// printf("%d ", it); +// } +// printf("\n"); +// } diff --git a/TRD/AliTRDSignalIndex.h b/TRD/AliTRDSignalIndex.h new file mode 100644 index 00000000000..255ced22b75 --- /dev/null +++ b/TRD/AliTRDSignalIndex.h @@ -0,0 +1,159 @@ +#ifndef AliTRDSIGNALINDEX_H +#define AliTRDSIGNALINDEX_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +#include "TObject.h" + +///////////////////////////////////////////////////////////// +// General container for data from TRD detector segments // +// Adapted from AliDigits, origin M.Ivanov // +///////////////////////////////////////////////////////////// + +//class TArrayI; +#include "TArrayI.h" + +class AliTRDSignalIndex : public TObject +{ + public: + + AliTRDSignalIndex(); + AliTRDSignalIndex(Int_t nrow, Int_t ncol,Int_t ntime); + AliTRDSignalIndex(const AliTRDSignalIndex &d); + virtual ~AliTRDSignalIndex(); // destructor + AliTRDSignalIndex &operator=(const AliTRDSignalIndex &d); + virtual void Copy(TObject &d) const; + virtual void Allocate(Int_t nrow, Int_t ncol,Int_t ntime); + virtual void Reset(); + + virtual void ResetCounters() + { + // reset the counters/iterators + fPositionRow = 0; + fPositionCol = fPositionRow + 1; + fPositionTbin = 1; + + fLastRow = -1; + fLastCol = -1; + fLastTbin = -1; + + fResetCounters = kTRUE; + } + + virtual void ResetTbinCounter() + { + // reset the time bin counter + + fPositionTbin = 1; + } + + virtual void AddIndexTBin(Int_t row, Int_t col, Int_t tbin); + + Bool_t NextRCIndex(Int_t &row, Int_t &col); // get the next pad (row and column) and return kTRUE on success + Bool_t NextRCTbinIndex(Int_t &row, Int_t &col, Int_t &tbin); // get the next timebin of a pad (row and column) and return kTRUE on success + Bool_t NextTbinIndex(Int_t &tbin); // get the next active timebin and return kTRUE on success + + Int_t GetCurrentRow() {return (*fIndex)[fPositionRow];} // current row + Int_t GetCurrentCol() {return (*fIndex)[fPositionCol];} // current col + Int_t GetCurtentTbin() {return (*fIndex)[fPositionCol + fPositionTbin];} //current tbin + + void ClearAll(); // clear the array, actually destroy and recreate w/o allocating + + //void Dump(); // printf content - one way of iterating demo + //void Dump2(); // printf content - another way of iterating demo + + //Bool_t IsAllocated() const {if (fIndex) return kTRUE; else return kFALSE;} + Bool_t IsAllocated() const + { + // return kTRUE if array allocated and there is no need to call allocate + if (!fIndex) + return kFALSE; + if (fIndex->GetSize() <= 0) + return kFALSE; + else return kTRUE; + } + + void SetSM(Int_t ix) + { + // Set which SM + fSM = ix; + }; + + void SetStack(Int_t ix) + { + // Set which stack + fStack = ix; + }; + + void SetChamber(Int_t ix) + { + // aka set stack + SetStack(ix); + } + + void SetLayer(Int_t ix) + { + // Set which layer + fLayer = ix; + }; + + void SetPlane(Int_t ix) + { + // aka set plane + SetLayer(ix); + } + + void SetDetNumber(Int_t ix) + { + // Set Det Number + fDet = ix; + } + + const Int_t GetDetNumber() {return fDet;} // Get Det number + const Int_t GetLayer() {return fLayer;} // Layer = Plane = position of the chamber in TRD + const Int_t GetPlane() {return fLayer;} // Layer = Plane = position of the chamber in TRD + const Int_t GetStack() {return fStack;} // Stack = Chameber = position of the chamber in TRD + const Int_t GetChamber() {return fStack;} // Stack = Chameber = position of the chamber in TRD + const Int_t GetSM() {return fSM;} // Super module of the TRD + + const Bool_t HasEntry() {return fHasEntry;} // Return status if has an entry + + TArrayI *GetArray() {return fIndex;} // Get the tarrayi pointer for god knows what reason + + private: + + Int_t fDet; // det number + Int_t fLayer; // aka plane - position in the full TRD + Int_t fStack; // aka chamber - position in the full TRD + Int_t fSM; // super module - position in the full TRD + + protected: + + TArrayI *fIndex; //! monitor active pads and tbins + + Int_t fPositionRow; // position in the index - jumps by 1 + 1 + fNtbins + Int_t fPositionCol; // position in the index - jumps by 1 + 1 + fNtbins + Int_t fPositionTbin; // position in the tbin - goes from 0 to fNtbins + + Int_t fLastRow; // to keep track what is the RC combination + Int_t fLastCol; // to keep track what is the RC combination + Int_t fLastTbin; // to keep track what is the Tbin - will catch if raw data bogus + + Int_t fNrows; // number of rows in the chamber + Int_t fNcols; // number of cols in the chamber + Int_t fNtbins; // number of tbins in the chamber + + Int_t fMaxLimit; // max number of things in the array = nrow * ncol * ntime + nrow * ncol * 2 + + Bool_t fResetCounters; // reset counter status + + Bool_t fHasEntry; // kTRUE flag if we have an entry + + ClassDef(AliTRDSignalIndex,1) // Data container for one TRD detector segment + +}; + +#endif diff --git a/TRD/AliTRDclusterizerV2.cxx b/TRD/AliTRDclusterizerV2.cxx new file mode 100644 index 00000000000..56c8fc6cc9f --- /dev/null +++ b/TRD/AliTRDclusterizerV2.cxx @@ -0,0 +1,1040 @@ + +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ + +/////////////////////////////////////////////////////////////////////////////// +// // +// TRD cluster finder // +// // +/////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include + +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliRawReader.h" +#include "AliLog.h" +#include "AliAlignObj.h" + +#include "AliTRDclusterizerV2.h" +#include "AliTRDgeometry.h" +#include "AliTRDdataArrayF.h" +#include "AliTRDdataArrayI.h" +#include "AliTRDdigitsManager.h" +#include "AliTRDpadPlane.h" +#include "AliTRDrawData.h" +#include "AliTRDcalibDB.h" +#include "AliTRDSimParam.h" +#include "AliTRDRecParam.h" +#include "AliTRDcluster.h" + +#include "Cal/AliTRDCalROC.h" +#include "Cal/AliTRDCalDet.h" + +#include "AliTRDSignalIndex.h" +#include "AliTRDRawStream.h" + +ClassImp(AliTRDclusterizerV2) + +//_____________________________________________________________________________ +AliTRDclusterizerV2::AliTRDclusterizerV2() + :AliTRDclusterizer() + ,fDigitsManager(NULL) + ,fGeometry(NULL) + ,fAddLabels(kTRUE) +{ + // + // AliTRDclusterizerV2 default constructor + // + +} + +//_____________________________________________________________________________ +AliTRDclusterizerV2::AliTRDclusterizerV2(const Text_t *name, const Text_t *title) + :AliTRDclusterizer(name,title) + ,fDigitsManager(new AliTRDdigitsManager()) + ,fGeometry(NULL) + ,fAddLabels(kTRUE) +{ + // + // AliTRDclusterizerV2 constructor + // + + fDigitsManager->CreateArrays(); + fGeometry = new AliTRDgeometry; +} + +//_____________________________________________________________________________ +AliTRDclusterizerV2::AliTRDclusterizerV2(const AliTRDclusterizerV2 &c) + :AliTRDclusterizer(c) + ,fDigitsManager(NULL) + ,fGeometry(NULL) + ,fAddLabels(kTRUE) +{ + // + // AliTRDclusterizerV2 copy constructor + // + +} + +//_____________________________________________________________________________ +AliTRDclusterizerV2::~AliTRDclusterizerV2() +{ + // + // AliTRDclusterizerV2 destructor + // + + if (fDigitsManager) { + delete fDigitsManager; + fDigitsManager = NULL; + } + + if (fGeometry) + { + delete fGeometry; + fGeometry = NULL; + } +} + +//_____________________________________________________________________________ +AliTRDclusterizerV2 &AliTRDclusterizerV2::operator=(const AliTRDclusterizerV2 &c) +{ + // + // Assignment operator + // + + if (this != &c) ((AliTRDclusterizerV2 &) c).Copy(*this); + return *this; + +} + +//_____________________________________________________________________________ +void AliTRDclusterizerV2::Copy(TObject &c) const +{ + // + // Copy function + // + + ((AliTRDclusterizerV2 &) c).fDigitsManager = 0; + + AliTRDclusterizer::Copy(c); + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::ReadDigits() +{ + // + // Reads the digits arrays from the input aliroot file + // + + if (!fRunLoader) { + AliError("No run loader available"); + return kFALSE; + } + + AliLoader* loader = fRunLoader->GetLoader("TRDLoader"); + if (!loader->TreeD()) { + loader->LoadDigits(); + } + + // Read in the digit arrays + return (fDigitsManager->ReadDigits(loader->TreeD())); + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::ReadDigits(TTree *digitsTree) +{ + // + // Reads the digits arrays from the input tree + // + + // Read in the digit arrays + return (fDigitsManager->ReadDigits(digitsTree)); + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::ReadDigits(AliRawReader *rawReader) +{ + // + // Reads the digits arrays from the ddl file + // + + AliTRDrawData raw; + fDigitsManager = raw.Raw2Digits(rawReader); + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::MakeClusters() +{ + // + // Creates clusters from digits + // + + Bool_t fReturn = kTRUE; + for (Int_t i = 0; i < AliTRDgeometry::kNdet; i++) + { + AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(i); + // This is to take care of switched off super modules + if (digitsIn->GetNtime() == 0) { + continue; + } + //AliInfo(Form("digitsIn->Expand() 0x%x", digitsIn)); + digitsIn->Expand(); + AliTRDSignalIndex* indexes = fDigitsManager->GetIndexes(i); + if (indexes->IsAllocated() == kFALSE) + { + fDigitsManager->BuildIndexes(i); + } + + if (indexes->HasEntry()) + { + if (fAddLabels) + { + for (Int_t iDict = 0; iDict < AliTRDdigitsManager::kNDict; iDict++) + { + AliTRDdataArrayI *tracksIn = 0; + tracksIn = fDigitsManager->GetDictionary(i,iDict); + tracksIn->Expand(); + } + } + Bool_t fR = MakeClusters(i); + fReturn = fR && fReturn; + } + //digitsIn->Compress(1,0); + // no compress just remove + fDigitsManager->RemoveDigits(i); + fDigitsManager->RemoveDictionaries(i); + fDigitsManager->ClearIndexes(i); + } + return fReturn; +} + + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::Raw2Clusters(AliRawReader *rawReader) +{ + // + // Creates clusters from raw data + // + + //AliDebug(2, "Read all"); + + AliTRDdataArrayI *digits = 0; + AliTRDdataArrayI *track0 = 0; + AliTRDdataArrayI *track1 = 0; + AliTRDdataArrayI *track2 = 0; + + AliTRDSignalIndex *indexes = 0; + // Create the digits manager + if (!fDigitsManager) + { + fDigitsManager = new AliTRDdigitsManager(); + fDigitsManager->CreateArrays(); + } + //delete fDigitsManager; + + AliTRDRawStream input(rawReader); + //input.SetRawVersion( fRawVersion ); + input.SetRawVersion( 2 ); + input.Init(); + + // Loop through the digits + Int_t lastdet = -1; + Int_t det = 0; + Int_t it = 0; + while (input.Next()) + { + + det = input.GetDet(); + + if (det != lastdet) + { + if (lastdet != -1) + { + digits = fDigitsManager->GetDigits(lastdet); +// if (digits->GetNtime() == 0) +// { +// AliDebug(5, Form("No Ntime? %d", lastdet)); +// } +// else +// { +// AliDebug(5, Form("Ntime ok %d", lastdet)); +// } + if (indexes->HasEntry()) + MakeClusters(lastdet); + } + + if (digits) + { + fDigitsManager->RemoveDigits(lastdet); + fDigitsManager->RemoveDictionaries(lastdet); + fDigitsManager->ClearIndexes(lastdet); + } + + lastdet = det; + + // Add a container for the digits of this detector + digits = fDigitsManager->GetDigits(det); + track0 = fDigitsManager->GetDictionary(det,0); + track1 = fDigitsManager->GetDictionary(det,1); + track2 = fDigitsManager->GetDictionary(det,2); + + // Allocate memory space for the digits buffer + if (digits->GetNtime() == 0) + { + //AliDebug(5, Form("Alloc digits for det %d", det)); + digits->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); + track0->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); + track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); + track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); + } + + indexes = fDigitsManager->GetIndexes(det); + indexes->SetSM(input.GetSM()); + indexes->SetStack(input.GetStack()); + indexes->SetLayer(input.GetLayer()); + indexes->SetDetNumber(det); + if (indexes->IsAllocated() == kFALSE) + indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins()); + } + + for (it = 0; it < 3; it++) + { + if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() ) + { + if (input.GetSignals()[it] > 0) + { + digits->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, input.GetSignals()[it]); + + indexes->AddIndexTBin(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it); + track0->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, 0); + track1->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, 0); + track2->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, 0); + } + } + } + } + + if (lastdet != -1) + { + MakeClusters(lastdet); + if (digits) + { + fDigitsManager->RemoveDigits(lastdet); + fDigitsManager->RemoveDictionaries(lastdet); + fDigitsManager->ClearIndexes(lastdet); + } + } + +// fDigitsManager->Delete(); +// delete fDigitsManager; +// fDigitsManager = NULL; + return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::Raw2ClustersChamber(AliRawReader *rawReader) +{ + // + // Creates clusters from raw data + // + + //AliDebug(1, "Raw2ClustersChamber"); + + // Create the digits manager + // Create the digits manager + if (!fDigitsManager) + { + fDigitsManager = new AliTRDdigitsManager(); + fDigitsManager->CreateArrays(); + } + //delete fDigitsManager; + + + AliTRDRawStream input(rawReader); + //input.SetRawVersion( fRawVersion ); + input.SetRawVersion( 2 ); + input.Init(); + + Int_t det = 0; + while ((det = input.NextChamber(fDigitsManager)) >= 0) + { + if (fDigitsManager->GetIndexes(det)->HasEntry()) + MakeClusters(det); + fDigitsManager->RemoveDigits(det); + fDigitsManager->RemoveDictionaries(det); + fDigitsManager->ClearIndexes(det); + } + +// delete fDigitsManager; +// fDigitsManager = NULL; + return kTRUE; +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::MakeClusters(Int_t det) +{ + // + // Generates the cluster. + // + + // Get the digits + // digits should be expanded beforehand! + // digitsIn->Expand(); + AliTRDdataArrayI *digitsIn = fDigitsManager->GetDigits(det); + // This is to take care of switched off super modules + if (digitsIn->GetNtime() == 0) + { + //AliDebug(5, Form("digitsIn->GetNtime() == 0 [%d]", det)); + return kFALSE; + } + + AliTRDSignalIndex *indexesIn = fDigitsManager->GetIndexes(det); + if (indexesIn->IsAllocated() == kFALSE) + { + AliError("Indexes do not exist!"); + return kFALSE; + } + + AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); + if (!calibration) { + AliFatal("No AliTRDcalibDB instance available\n"); + return kFALSE; + } + + AliTRDSimParam *simParam = AliTRDSimParam::Instance(); + if (!simParam) { + AliError("No AliTRDSimParam instance available\n"); + return kFALSE; + } + + AliTRDRecParam *recParam = AliTRDRecParam::Instance(); + if (!recParam) { + AliError("No AliTRDRecParam instance available\n"); + return kFALSE; + } + + // ADC thresholds + Float_t ADCthreshold = simParam->GetADCthreshold(); + // Threshold value for the maximum + Float_t maxThresh = recParam->GetClusMaxThresh(); + // Threshold value for the digit signal + Float_t sigThresh = recParam->GetClusSigThresh(); + + // Detector wise calibration object for t0 + const AliTRDCalDet *calT0Det = calibration->GetT0Det(); + // Detector wise calibration object for the gain factors + const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet(); + + // Iteration limit for unfolding procedure + const Float_t kEpsilon = 0.01; + const Int_t kNclus = 3; + const Int_t kNsig = 5; + + Int_t iUnfold = 0; + Double_t ratioLeft = 1.0; + Double_t ratioRight = 1.0; + + Double_t padSignal[kNsig]; + Double_t clusterSignal[kNclus]; + Double_t clusterPads[kNclus]; + + Int_t icham = indexesIn->GetChamber(); + Int_t iplan = indexesIn->GetPlane(); + Int_t isect = indexesIn->GetSM(); + + // Start clustering in the chamber + + Int_t idet = fGeometry->GetDetector(iplan,icham,isect); + if (idet != det) + { + AliError("Strange Detector number Missmatch!"); + return kFALSE; + } + + Int_t ilayer = AliGeomManager::kTRD1 + iplan; + Int_t imodule = icham + AliTRDgeometry::Ncham() * isect; + UShort_t volid = AliGeomManager::LayerToVolUID(ilayer,imodule); + + Int_t nRowMax = digitsIn->GetNrow(); + Int_t nColMax = digitsIn->GetNcol(); + Int_t nTimeTotal = digitsIn->GetNtime(); + + AliTRDpadPlane *padPlane = fGeometry->GetPadPlane(iplan,icham); + + // Calibration object with pad wise values for t0 + AliTRDCalROC *calT0ROC = calibration->GetT0ROC(idet); + // Calibration object with pad wise values for the gain factors + AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet); + // Calibration value for chamber wise t0 + Float_t calT0DetValue = calT0Det->GetValue(idet); + // Calibration value for chamber wise gain factor + Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet); + + Int_t nClusters = 0; + + // Apply the gain and the tail cancelation via digital filter + AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow() + ,digitsIn->GetNcol() + ,digitsIn->GetNtime()); + + AliTRDSignalIndex *indexesOut = new AliTRDSignalIndex(digitsIn->GetNrow() + ,digitsIn->GetNcol() + ,digitsIn->GetNtime()); + + AliTRDSignalIndex *indexesMaxima = new AliTRDSignalIndex(digitsIn->GetNrow() + ,digitsIn->GetNcol() + ,digitsIn->GetNtime()); + + + Transform(digitsIn + ,digitsOut + ,indexesIn + ,indexesOut + ,nRowMax,nColMax,nTimeTotal + ,ADCthreshold + ,calGainFactorROC + ,calGainFactorDetValue); + + Int_t row = 0; + Int_t col = 0; + Int_t time = 0; + Int_t iPad = 0; + + indexesOut->ResetCounters(); + while (indexesOut->NextRCTbinIndex(row, col, time)) + { + Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time)); + + // Look for the maximum + if (signalM >= maxThresh) + { + + if (col + 1 >= nColMax || col-1 < 0) + continue; + + Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1 ,time)); + + Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); + + if ((TMath::Abs(signalL) <= signalM) && + (TMath::Abs(signalR) < signalM)) + { + if ((TMath::Abs(signalL) >= sigThresh) || + (TMath::Abs(signalR) >= sigThresh)) + { + // Maximum found, mark the position by a negative signal + digitsOut->SetDataUnchecked(row,col,time,-signalM); + indexesMaxima->AddIndexTBin(row,col,time); + } + } + } + } + + // The index to the first cluster of a given ROC + Int_t firstClusterROC = -1; + // The number of cluster in a given ROC + Int_t nClusterROC = 0; + + // Now check the maxima and calculate the cluster position + indexesMaxima->ResetCounters(); + while (indexesMaxima->NextRCTbinIndex(row, col, time)) + { + // Maximum found ? + if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) { + + for (iPad = 0; iPad < kNclus; iPad++) { + Int_t iPadCol = col - 1 + iPad; + clusterSignal[iPad] = + TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time)); + } + + // Count the number of pads in the cluster + Int_t nPadCount = 0; + Int_t ii; + // Look to the left + ii = 0; + while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) { + nPadCount++; + ii++; + if (col-ii < 0) break; + } + // Look to the right + ii = 0; + while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) { + nPadCount++; + ii++; + if (col+ii+1 >= nColMax) break; + } + nClusters++; + + // Look for 5 pad cluster with minimum in the middle + Bool_t fivePadCluster = kFALSE; + if (col < (nColMax - 3)) { + if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) { + fivePadCluster = kTRUE; + } + if ((fivePadCluster) && (col < (nColMax - 5))) { + if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) { + fivePadCluster = kFALSE; + } + } + if ((fivePadCluster) && (col > 1)) { + if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) { + fivePadCluster = kFALSE; + } + } + } + + // 5 pad cluster + // Modify the signal of the overlapping pad for the left part + // of the cluster which remains from a previous unfolding + if (iUnfold) { + clusterSignal[0] *= ratioLeft; + iUnfold = 0; + } + + // Unfold the 5 pad cluster + if (fivePadCluster) { + for (iPad = 0; iPad < kNsig; iPad++) { + padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row + ,col-1+iPad + ,time)); + } + // Unfold the two maxima and set the signal on + // the overlapping pad to the ratio + ratioRight = Unfold(kEpsilon,iplan,padSignal); + ratioLeft = 1.0 - ratioRight; + clusterSignal[2] *= ratioRight; + iUnfold = 1; + } + + Double_t clusterCharge = clusterSignal[0] + + clusterSignal[1] + + clusterSignal[2]; + + // The position of the cluster + clusterPads[0] = row + 0.5; + // Take the shift of the additional time bins into account + clusterPads[2] = time + 0.5; + + if (recParam->LUTOn()) { + // Calculate the position of the cluster by using the + // lookup table method + clusterPads[1] = recParam->LUTposition(iplan,clusterSignal[0] + ,clusterSignal[1] + ,clusterSignal[2]); + } + else { + // Calculate the position of the cluster by using the + // center of gravity method + for (Int_t i = 0; i < kNsig; i++) { + padSignal[i] = 0.0; + } + padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad + padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left pad + padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right pad + if ((col > 2) && + (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) { + padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)); + } + if ((col < nColMax - 3) && + (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])) { + padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)); + } + clusterPads[1] = GetCOG(padSignal); + } + + Double_t q0 = clusterSignal[0]; + Double_t q1 = clusterSignal[1]; + Double_t q2 = clusterSignal[2]; + Double_t clusterSigmaY2 = (q1 * (q0 + q2) + 4.0 * q0 * q2) + / (clusterCharge*clusterCharge); + + // + // Calculate the position and the error + // + + // Correct for t0 (sum of chamber and pad wise values !!!) + Float_t calT0ROCValue = calT0ROC->GetValue(col,row); + Char_t clusterTimeBin = ((Char_t) TMath::Nint(time - (calT0DetValue + calT0ROCValue))); + Double_t colSize = padPlane->GetColSize(col); + Double_t rowSize = padPlane->GetRowSize(row); + + Float_t clusterPos[3]; + clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1] + 0.5) * colSize; + clusterPos[1] = padPlane->GetRowPos(row) - 0.5 * rowSize; + clusterPos[2] = CalcXposFromTimebin(clusterPads[2],idet,col,row); + Float_t clusterSig[2]; + clusterSig[0] = (clusterSigmaY2 + 1.0/12.0) * colSize*colSize; + clusterSig[1] = rowSize * rowSize / 12.0; + + // Store the amplitudes of the pads in the cluster for later analysis + Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 }; + for (Int_t jPad = col-3; jPad <= col+3; jPad++) { + if ((jPad < 0) || + (jPad >= nColMax-1)) { + continue; + } + signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time))); + } + + // Add the cluster to the output array + // The track indices will be stored later + AliTRDcluster *cluster = new AliTRDcluster(idet + ,clusterCharge + ,clusterPos + ,clusterSig + ,0x0 + ,((Char_t) nPadCount) + ,signals + ,((UChar_t) col) + ,clusterTimeBin + ,clusterPads[1] + ,volid); + // Temporarily store the row, column and time bin of the center pad + // Used to later on assign the track indices + cluster->SetLabel( row,0); + cluster->SetLabel( col,1); + cluster->SetLabel(time,2); + RecPoints()->Add(cluster); + + // Store the index of the first cluster in the current ROC + if (firstClusterROC < 0) { + firstClusterROC = RecPoints()->GetEntriesFast() - 1; + } + // Count the number of cluster in the current ROC + nClusterROC++; + + } // if: Maximum found ? + + } + delete digitsOut; + delete indexesOut; + delete indexesMaxima; + + if (fAddLabels) + AddLabels(idet, firstClusterROC, nClusterROC); + + // Write the cluster and reset the array + WriteClusters(idet); + ResetRecPoints(); + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV2::AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC) +{ + // + // Add the track indices to the found clusters + // + + const Int_t kNclus = 3; + const Int_t kNdict = AliTRDdigitsManager::kNDict; + const Int_t kNtrack = kNdict * kNclus; + + Int_t iClusterROC = 0; + + Int_t row = 0; + Int_t col = 0; + Int_t time = 0; + Int_t iPad = 0; + + // Temporary array to collect the track indices + Int_t *idxTracks = new Int_t[kNtrack*nClusterROC]; + + // Loop through the dictionary arrays one-by-one + // to keep memory consumption low + AliTRDdataArrayI *tracksIn = 0; + for (Int_t iDict = 0; iDict < kNdict; iDict++) { + + tracksIn = fDigitsManager->GetDictionary(idet,iDict); + // tracksIn should be expanded beforehand! + + // Loop though the clusters found in this ROC + for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) { + + AliTRDcluster *cluster = (AliTRDcluster *) + RecPoints()->UncheckedAt(firstClusterROC+iClusterROC); + row = cluster->GetLabel(0); + col = cluster->GetLabel(1); + time = cluster->GetLabel(2); + + for (iPad = 0; iPad < kNclus; iPad++) { + Int_t iPadCol = col - 1 + iPad; + Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1; + idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index; + } + + } + + // Compress the arrays + // no do not compress - we will delete them when we are done with the detector + //tracksIn->Compress(1,0); + + } + + // Copy the track indices into the cluster + // Loop though the clusters found in this ROC + for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) { + + AliTRDcluster *cluster = (AliTRDcluster *) + RecPoints()->UncheckedAt(firstClusterROC+iClusterROC); + cluster->SetLabel(-9999,0); + cluster->SetLabel(-9999,1); + cluster->SetLabel(-9999,2); + + cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]); + + } + + delete [] idxTracks; + + return kTRUE; +} + +//_____________________________________________________________________________ +Double_t AliTRDclusterizerV2::GetCOG(Double_t signal[5]) +{ + // + // Get COG position + // Used for clusters with more than 3 pads - where LUT not applicable + // + + Double_t sum = signal[0] + + signal[1] + + signal[2] + + signal[3] + + signal[4]; + + Double_t res = (0.0 * (-signal[0] + signal[4]) + + (-signal[1] + signal[3])) / sum; + + return res; + +} + +//_____________________________________________________________________________ +Double_t AliTRDclusterizerV2::Unfold(Double_t eps, Int_t plane, Double_t *padSignal) +{ + // + // Method to unfold neighbouring maxima. + // The charge ratio on the overlapping pad is calculated + // until there is no more change within the range given by eps. + // The resulting ratio is then returned to the calling method. + // + + AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); + if (!calibration) { + AliError("No AliTRDcalibDB instance available\n"); + return kFALSE; + } + + Int_t irc = 0; + Int_t itStep = 0; // Count iteration steps + + Double_t ratio = 0.5; // Start value for ratio + Double_t prevRatio = 0.0; // Store previous ratio + + Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal + Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal + Double_t newSignal[3] = { 0.0, 0.0, 0.0 }; + + // Start the iteration + while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { + + itStep++; + prevRatio = ratio; + + // Cluster position according to charge ratio + Double_t maxLeft = (ratio*padSignal[2] - padSignal[0]) + / (padSignal[0] + padSignal[1] + ratio*padSignal[2]); + Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) + / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]); + + // Set cluster charge ratio + irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal); + Double_t ampLeft = padSignal[1] / newSignal[1]; + irc = calibration->PadResponse(1.0,maxRight,plane,newSignal); + Double_t ampRight = padSignal[3] / newSignal[1]; + + // Apply pad response to parameters + irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal ); + irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal); + + // Calculate new overlapping ratio + ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] / + (newLeftSignal[2] + newRightSignal[0])); + + } + + return ratio; + +} + +//_____________________________________________________________________________ +void AliTRDclusterizerV2::Transform(AliTRDdataArrayI *digitsIn + , AliTRDdataArrayF *digitsOut + , AliTRDSignalIndex *indexesIn + , AliTRDSignalIndex *indexesOut + , Int_t /*nRowMax*/, Int_t /*nColMax*/ + , Int_t nTimeTotal + , Float_t ADCthreshold + , AliTRDCalROC *calGainFactorROC + , Float_t calGainFactorDetValue) +{ + // + // Apply gain factor + // Apply tail cancelation: Transform digitsIn to digitsOut + // + + Int_t iRow = 0; + Int_t iCol = 0; + Int_t iTime = 0; + + AliTRDRecParam *recParam = AliTRDRecParam::Instance(); + if (!recParam) { + AliError("No AliTRDRecParam instance available\n"); + return; + } + + Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation + Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation + indexesIn->ResetCounters(); + while (indexesIn->NextRCIndex(iRow, iCol)) + { + Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow); + Double_t gain = calGainFactorDetValue + * calGainFactorROCValue; + + for (iTime = 0; iTime < nTimeTotal; iTime++) + { + + // + // Add gain + // + inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime); + inADC[iTime] /= gain; + outADC[iTime] = inADC[iTime]; + } + + // Apply the tail cancelation via the digital filter + if (recParam->TCOn()) { + DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp()); + } + indexesIn->ResetTbinCounter(); + while (indexesIn->NextTbinIndex(iTime)) + { + // Store the amplitude of the digit if above threshold + if (outADC[iTime] > ADCthreshold) + { + digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]); + AliDebug(5, Form("add index %d", indexesIn->GetDetNumber())); + indexesOut->AddIndexTBin(iRow,iCol,iTime); + } + } //while itime + }//while irow icol + + delete [] inADC; + delete [] outADC; + + AliDebug(5, Form("Stop %d", indexesIn->GetDetNumber())); + + return; + +} + +//_____________________________________________________________________________ +void AliTRDclusterizerV2::DeConvExp(Double_t *source, Double_t *target + , Int_t n, Int_t nexp) +{ + // + // Tail cancellation by deconvolution for PASA v4 TRF + // + + Double_t rates[2]; + Double_t coefficients[2]; + + // Initialization (coefficient = alpha, rates = lambda) + Double_t R1 = 1.0; + Double_t R2 = 1.0; + Double_t C1 = 0.5; + Double_t C2 = 0.5; + + if (nexp == 1) { // 1 Exponentials + R1 = 1.156; + R2 = 0.130; + C1 = 0.066; + C2 = 0.000; + } + if (nexp == 2) { // 2 Exponentials + R1 = 1.156; + R2 = 0.130; + C1 = 0.114; + C2 = 0.624; + } + + coefficients[0] = C1; + coefficients[1] = C2; + + Double_t Dt = 0.1; + + rates[0] = TMath::Exp(-Dt/(R1)); + rates[1] = TMath::Exp(-Dt/(R2)); + + Int_t i = 0; + Int_t k = 0; + + Double_t reminder[2]; + Double_t correction; + Double_t result; + + // Attention: computation order is important + correction = 0.0; + for (k = 0; k < nexp; k++) { + reminder[k] = 0.0; + } + for (i = 0; i < n; i++) { + result = (source[i] - correction); // No rescaling + target[i] = result; + + for (k = 0; k < nexp; k++) { + reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result); + } + correction = 0.0; + for (k = 0; k < nexp; k++) { + correction += reminder[k]; + } + } + +} diff --git a/TRD/AliTRDclusterizerV2.h b/TRD/AliTRDclusterizerV2.h new file mode 100644 index 00000000000..9ede011f7aa --- /dev/null +++ b/TRD/AliTRDclusterizerV2.h @@ -0,0 +1,75 @@ +#ifndef ALITRDCLUSTERIZERV2_H +#define ALITRDCLUSTERIZERV2_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$ */ + +//////////////////////////////////////////////////////////////////////////// +// // +// TRD cluster finder // +// // +//////////////////////////////////////////////////////////////////////////// + +#include "AliTRDclusterizer.h" + +class AliTRDdataArrayI; +class AliTRDdataArrayF; +class AliTRDdigitsManager; +class AliTRDCalROC; +class AliRawReader; +class AliTRDSignalIndex; +class AliTRDgeometry; + +class AliTRDclusterizerV2 : public AliTRDclusterizer { + + public: + + AliTRDclusterizerV2(); + AliTRDclusterizerV2(const Text_t* name, const Text_t* title); + AliTRDclusterizerV2(const AliTRDclusterizerV2 &c); + virtual ~AliTRDclusterizerV2(); + AliTRDclusterizerV2 &operator=(const AliTRDclusterizerV2 &c); + + virtual void Copy(TObject &c) const; + virtual Bool_t Raw2Clusters(AliRawReader *rawReader); + virtual Bool_t Raw2ClustersChamber(AliRawReader *rawReader); + virtual Bool_t MakeClusters(); + virtual Bool_t MakeClusters(Int_t det); + virtual Bool_t ReadDigits(); + virtual Bool_t ReadDigits(AliRawReader *rawReader); + virtual Bool_t ReadDigits(TTree *digitsTree); + + virtual Bool_t AddLabels(Int_t idet, Int_t firstClusterROC, Int_t nClusterROC); + + virtual Bool_t SetAddLabels(Bool_t kset) { fAddLabels = kset; return fAddLabels;} + + protected: + + void DeConvExp(Double_t *source, Double_t *target + , Int_t nTimeTotal, Int_t nexp); + void Transform(AliTRDdataArrayI *digitsIn + , AliTRDdataArrayF *digitsOut + , AliTRDSignalIndex *indexesIn + , AliTRDSignalIndex *indexesOut + , Int_t nRowMax, Int_t nColMax, Int_t nTimeTotal + , Float_t ADCthreshold + , AliTRDCalROC *calGainFactorROC + , Float_t calGainFactorDetValue); +/* void Transform(AliTRDdataArrayI *digitsIn, AliTRDdataArrayF *digitsOut */ +/* , Int_t nRowMax, Int_t nColMax, Int_t nTimeTotal */ +/* , Float_t ADCthreshold */ +/* , AliTRDCalROC *calGainFactorROC */ +/* , Float_t calGainFactorDetValue); */ + virtual Double_t Unfold(Double_t eps, Int_t plane, Double_t *padSignal); + Double_t GetCOG(Double_t signal[5]); + + AliTRDdigitsManager *fDigitsManager; //! TRD digits manager + AliTRDgeometry *fGeometry; //! default TRD geometry + + Bool_t fAddLabels; // should clusters have MC labels? + ClassDef(AliTRDclusterizerV2,1) // TRD-Cluster finder, slow simulator + +}; + +#endif diff --git a/TRD/AliTRDdigitsManager.cxx b/TRD/AliTRDdigitsManager.cxx index 3e8fe6b7d42..7c759249e7d 100644 --- a/TRD/AliTRDdigitsManager.cxx +++ b/TRD/AliTRDdigitsManager.cxx @@ -37,6 +37,8 @@ #include "AliTRDdigit.h" #include "AliTRDgeometry.h" +#include "AliTRDSignalIndex.h" + ClassImp(AliTRDdigitsManager) //_____________________________________________________________________________ @@ -52,6 +54,7 @@ AliTRDdigitsManager::AliTRDdigitsManager() ,fDigits(0) ,fIsRaw(0) ,fSDigits(0) + ,fSignalIndexes(NULL) { // // Default constructor @@ -60,7 +63,10 @@ AliTRDdigitsManager::AliTRDdigitsManager() for (Int_t iDict = 0; iDict < kNDict; iDict++) { fDictionary[iDict] = NULL; } - + + //fSignalIndexes = new TList(); + fSignalIndexes = new TObjArray(AliTRDgeometry::Ndet()); + } //_____________________________________________________________________________ @@ -71,6 +77,7 @@ AliTRDdigitsManager::AliTRDdigitsManager(const AliTRDdigitsManager &m) ,fDigits(0) ,fIsRaw(m.fIsRaw) ,fSDigits(m.fSDigits) + ,fSignalIndexes(NULL) { // // AliTRDdigitsManager copy constructor @@ -97,6 +104,11 @@ AliTRDdigitsManager::~AliTRDdigitsManager() fDictionary[iDict] = NULL; } + delete fSignalIndexes; + fSignalIndexes = NULL; +// for (Int_t i = 0; i < AliTRDgeometry::Ndet(); i++) +// delete fSignalIndexes[i]; + } //_____________________________________________________________________________ @@ -121,7 +133,8 @@ void AliTRDdigitsManager::Copy(TObject &m) const ((AliTRDdigitsManager &) m).fIsRaw = fIsRaw; ((AliTRDdigitsManager &) m).fEvent = fEvent; ((AliTRDdigitsManager &) m).fSDigits = fSDigits; - + + ((AliTRDdigitsManager &) m).fSignalIndexes = fSignalIndexes; TObject::Copy(m); } @@ -140,6 +153,10 @@ void AliTRDdigitsManager::CreateArrays() ,AliTRDgeometry::Ndet()); } + for (Int_t i = 0; i < AliTRDgeometry::Ndet(); i++) + { + fSignalIndexes->AddLast(new AliTRDSignalIndex()); + } } //_____________________________________________________________________________ void AliTRDdigitsManager::ResetArrays() @@ -161,6 +178,11 @@ void AliTRDdigitsManager::ResetArrays() ,AliTRDgeometry::Ndet()); } + for (Int_t i = 0; i < AliTRDgeometry::Ndet(); i++) + { + AliTRDSignalIndex *idx = (AliTRDSignalIndex *)fSignalIndexes->At(i); + idx->Reset(); + } } //_____________________________________________________________________________ @@ -396,3 +418,107 @@ Int_t AliTRDdigitsManager::GetTrack(Int_t track, AliTRDdigit *Digit) const return GetTrack(track,row,col,time,det); } + +//_____________________________________________________________________________ +AliTRDSignalIndex *AliTRDdigitsManager::GetIndexes(Int_t det) +{ + // + // Returns indexes of active pads + // + + return (AliTRDSignalIndex*)fSignalIndexes->At(det); + +} + +//_____________________________________________________________________________ +void AliTRDdigitsManager::RemoveDigits(Int_t det) +{ + // + // Clear memory + // + + fDigits->ClearSegment(det); + +} + +//_____________________________________________________________________________ +void AliTRDdigitsManager::RemoveDictionaries(Int_t det) +{ + // + // Clear memory + // + + for (Int_t i = 0; i < kNDict; i++) { + fDictionary[i]->ClearSegment(det); + } + +} + +//_____________________________________________________________________________ +void AliTRDdigitsManager::ClearIndexes(Int_t det) +{ + // + // Clear memory + // + fSignalIndexes->At(det)->Clear(); +} + + +//_____________________________________________________________________________ +Bool_t AliTRDdigitsManager::BuildIndexes(Int_t det) +{ + // + // Build the list of indices + // + + Int_t nRows = 0; + Int_t nCols = 0; + Int_t nTbins = 0; + + AliTRDgeometry geom; + AliTRDdataArrayI *digits = GetDigits(det); + //digits should be expanded by now!!! + if (digits->GetNtime() > 0) + { + digits->Expand(); + nRows = digits->GetNrow(); + nCols = digits->GetNcol(); + nTbins = digits->GetNtime(); + + //AliInfo(Form("rows %d cols %d tbins %d", nRows, nCols, nTbins)); + + AliTRDSignalIndex* indexes = GetIndexes(det); + indexes->SetSM(geom.GetSector(det)); + indexes->SetChamber(geom.GetChamber(det)); + indexes->SetPlane(geom.GetPlane(det)); + indexes->SetDetNumber(det); + if (indexes->IsAllocated() == kFALSE) + { + indexes->Allocate(nRows, nCols, nTbins); + //AliInfo(Form("Allocating 0x%x %d", indexes->GetArray(), indexes->GetArray()->GetSize())); + } + for (Int_t ir = 0; ir < nRows; ir++) + { + for (Int_t ic = 0; ic < nCols; ic++) + { + for (Int_t it = 0; it < nTbins; it++) + { + //AliInfo(Form("row %d col %d tbin %d", ir, ic, it)); + + Int_t isig = digits->GetDataUnchecked(ir, ic, it); + if (isig > 0) + { + //AliInfo(Form("row %d col %d tbin %d", ir, ic, it)); + indexes->AddIndexTBin(ir, ic, it); + } + } //tbins + } //cols + } // rows + } // if GetNtime + else + { + return kFALSE; + } + return kTRUE; + +} diff --git a/TRD/AliTRDdigitsManager.h b/TRD/AliTRDdigitsManager.h index d4769735667..d94a0766f6e 100644 --- a/TRD/AliTRDdigitsManager.h +++ b/TRD/AliTRDdigitsManager.h @@ -4,7 +4,7 @@ /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * See cxx source for full Copyright notice */ -/* $Id: AliTRDdigitsManager.h,v */ +/* $Id$ */ ///////////////////////////////////////////////////////////// // Manages the TRD digits // @@ -18,6 +18,7 @@ class TTree; class AliTRDsegmentArray; class AliTRDdataArrayI; class AliTRDdigit; +class AliTRDSignalIndex; class AliTRDdigitsManager : public TObject { @@ -55,9 +56,18 @@ class AliTRDdigitsManager : public TObject { AliTRDdataArrayI *GetDigits(Int_t det) const; AliTRDdataArrayI *GetDictionary(Int_t det, Int_t i) const; + + void RemoveDigits(Int_t det); + void RemoveDictionaries(Int_t det); + void ClearIndexes(Int_t det); + Int_t GetTrack(Int_t track, AliTRDdigit *Digit) const; Short_t GetDigitAmp(Int_t row, Int_t col, Int_t time, Int_t det) const; + AliTRDSignalIndex *GetIndexes(Int_t det); + TObjArray *GetIndexes() {return fSignalIndexes;}; + virtual Bool_t BuildIndexes(Int_t det); + protected: static const Int_t fgkNDict; // Number of track dictionary arrays @@ -72,7 +82,9 @@ class AliTRDdigitsManager : public TObject { Bool_t fIsRaw; // Flag indicating raw digits Bool_t fSDigits; // Switch for the summable digits - ClassDef(AliTRDdigitsManager,5) // Manages the TRD digits + TObjArray *fSignalIndexes; // Provides access to the active pads and tbins + + ClassDef(AliTRDdigitsManager,6) // Manages the TRD digits }; diff --git a/TRD/AliTRDrawData.cxx b/TRD/AliTRDrawData.cxx index 1cb00bf26c8..d01bd90cebe 100644 --- a/TRD/AliTRDrawData.cxx +++ b/TRD/AliTRDrawData.cxx @@ -38,6 +38,7 @@ #include "AliTRDcalibDB.h" #include "AliFstream.h" +#include "AliTRDSignalIndex.h" ClassImp(AliTRDrawData) //_____________________________________________________________________________ @@ -420,6 +421,7 @@ AliTRDdigitsManager *AliTRDrawData::Raw2Digits(AliRawReader *rawReader) AliTRDdataArrayI *track1 = 0; AliTRDdataArrayI *track2 = 0; + AliTRDSignalIndex *indexes = 0; // Create the digits manager AliTRDdigitsManager* digitsManager = new AliTRDdigitsManager(); digitsManager->CreateArrays(); @@ -461,21 +463,34 @@ AliTRDdigitsManager *AliTRDrawData::Raw2Digits(AliRawReader *rawReader) track1->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); track2->Allocate(input.GetMaxRow(),input.GetMaxCol(), input.GetNumberOfTimeBins()); } + + indexes = digitsManager->GetIndexes(det); + indexes->SetSM(input.GetSM()); + indexes->SetStack(input.GetStack()); + indexes->SetLayer(input.GetLayer()); + indexes->SetDetNumber(det); + if (indexes->IsAllocated() == kFALSE) + indexes->Allocate(input.GetMaxRow(), input.GetMaxCol(), input.GetNumberOfTimeBins()); } for (it = 0; it < 3; it++) { if ( input.GetTimeBin() + it < input.GetNumberOfTimeBins() ) { - digits->SetDataUnchecked(input.GetRow(), input.GetCol(), - input.GetTimeBin() + it, input.GetSignals()[it]); - track0->SetDataUnchecked(input.GetRow(), input.GetCol(), - input.GetTimeBin() + it, 0); - track1->SetDataUnchecked(input.GetRow(), input.GetCol(), - input.GetTimeBin() + it, 0); - track2->SetDataUnchecked(input.GetRow(), input.GetCol(), - input.GetTimeBin() + it, 0); - + if (input.GetSignals()[it] > 0) + { + digits->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, input.GetSignals()[it]); + + indexes->AddIndexTBin(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it); + track0->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, 0); + track1->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, 0); + track2->SetDataUnchecked(input.GetRow(), input.GetCol(), + input.GetTimeBin() + it, 0); + } } } } diff --git a/TRD/TRDbaseLinkDef.h b/TRD/TRDbaseLinkDef.h index c466b1246f0..58222ed916f 100644 --- a/TRD/TRDbaseLinkDef.h +++ b/TRD/TRDbaseLinkDef.h @@ -19,6 +19,8 @@ #pragma link C++ class AliTRDsegmentArrayBase+; #pragma link C++ class AliTRDsegmentArray+; +#pragma link C++ class AliTRDSignalIndex+; + #pragma link C++ class AliTRDgeometry+; #pragma link C++ class AliTRDpadPlane+; diff --git a/TRD/TRDrecLinkDef.h b/TRD/TRDrecLinkDef.h index 1dda033e462..50238d33ca5 100644 --- a/TRD/TRDrecLinkDef.h +++ b/TRD/TRDrecLinkDef.h @@ -13,6 +13,7 @@ #pragma link C++ class AliTRDclusterizer+; #pragma link C++ class AliTRDclusterizerV1+; +#pragma link C++ class AliTRDclusterizerV2+; #pragma link C++ class AliTRDclusterCorrection+; diff --git a/TRD/libTRDbase.pkg b/TRD/libTRDbase.pkg index 9eb38c52e3c..2c6b5c738b5 100644 --- a/TRD/libTRDbase.pkg +++ b/TRD/libTRDbase.pkg @@ -6,6 +6,7 @@ SRCS= AliTRDarrayI.cxx \ AliTRDdataArrayF.cxx \ AliTRDsegmentArrayBase.cxx \ AliTRDsegmentArray.cxx \ + AliTRDSignalIndex.cxx \ AliTRDgeometry.cxx \ AliTRDdigit.cxx \ AliTRDdigitsManager.cxx \ diff --git a/TRD/libTRDrec.pkg b/TRD/libTRDrec.pkg index 7938cbde194..35fa0d56dba 100644 --- a/TRD/libTRDrec.pkg +++ b/TRD/libTRDrec.pkg @@ -2,6 +2,7 @@ SRCS= AliTRDcluster.cxx \ AliTRDclusterMI.cxx \ AliTRDclusterizer.cxx \ AliTRDclusterizerV1.cxx \ + AliTRDclusterizerV2.cxx \ AliTRDclusterCorrection.cxx \ AliTRDtracklet.cxx \ AliTRDtrack.cxx \ -- 2.43.0