X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TRD%2FAliTRDmcmSim.cxx;h=162a623eafe3b0e16b10165acb80a125299a73dc;hp=a6a5fd4900b5e99fc1056650e7ff80a03e75bdc2;hb=9eee9fe98421cc964099854d967785b3c73c44f0;hpb=c8b1590d38895c0471accffb8a0bf536e4aba304 diff --git a/TRD/AliTRDmcmSim.cxx b/TRD/AliTRDmcmSim.cxx index a6a5fd4900b..162a623eafe 100644 --- a/TRD/AliTRDmcmSim.cxx +++ b/TRD/AliTRDmcmSim.cxx @@ -13,166 +13,162 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ + /////////////////////////////////////////////////////////////////////////////// // // // TRD MCM (Multi Chip Module) simulator // -// which simulated the TRAP processing after the AD-conversion // -// The relevant parameters (i.e. configuration registers of the TRAP // -// configuration are taken from AliTRDtrapConfig. // +// which simulates the TRAP processing after the AD-conversion. // +// The relevant parameters (i.e. configuration settings of the TRAP) // +// are taken from AliTRDtrapConfig. // // // /////////////////////////////////////////////////////////////////////////////// -#include // needed for raw data dump +#include +#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "TCanvas.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TGraph.h" +#include "TLine.h" +#include "TRandom.h" +#include "TClonesArray.h" +#include "TMath.h" +#include #include "AliLog.h" -#include "AliRun.h" #include "AliRunLoader.h" #include "AliLoader.h" -#include "AliTRDdigit.h" #include "AliTRDfeeParam.h" -#include "AliTRDtrapConfig.h" -#include "AliTRDSimParam.h" -#include "AliTRDgeometry.h" #include "AliTRDcalibDB.h" +#include "AliTRDtrapConfig.h" #include "AliTRDdigitsManager.h" #include "AliTRDarrayADC.h" #include "AliTRDarrayDictionary.h" -#include "AliTRDpadPlane.h" #include "AliTRDtrackletMCM.h" #include "AliTRDmcmSim.h" - -#include "AliMagF.h" -#include "TGeoGlobalMagField.h" +#include "TTreeStream.h" ClassImp(AliTRDmcmSim) Bool_t AliTRDmcmSim::fgApplyCut = kTRUE; - -//_____________________________________________________________________________ -AliTRDmcmSim::AliTRDmcmSim() : TObject() - ,fInitialized(kFALSE) - ,fMaxTracklets(-1) - ,fDetector(-1) - ,fRobPos(-1) - ,fMcmPos(-1) - ,fRow (-1) - ,fNADC(-1) - ,fNTimeBin(-1) - ,fADCR(NULL) - ,fADCF(NULL) - ,fMCMT(NULL) - ,fTrackletArray(NULL) - ,fZSM(NULL) - ,fZSM1Dim(NULL) - ,fFeeParam(NULL) - ,fTrapConfig(NULL) - ,fSimParam(NULL) - ,fCommonParam(NULL) - ,fCal(NULL) - ,fGeo(NULL) - ,fDigitsManager(NULL) - ,fPedAcc(NULL) - ,fGainCounterA(NULL) - ,fGainCounterB(NULL) - ,fTailAmplLong(NULL) - ,fTailAmplShort(NULL) - ,fNHits(0) - ,fFitReg(NULL) +Int_t AliTRDmcmSim::fgAddBaseline = 0; +Bool_t AliTRDmcmSim::fgStoreClusters = kFALSE; + +const Int_t AliTRDmcmSim::fgkFormatIndex = std::ios_base::xalloc(); + +const UShort_t AliTRDmcmSim::fgkFPshifts[4] = {11, 14, 17, 21}; + + +AliTRDmcmSim::AliTRDmcmSim() : + TObject(), + fInitialized(kFALSE), + fDetector(-1), + fRobPos(-1), + fMcmPos(-1), + fRow (-1), + fNTimeBin(-1), + fADCR(NULL), + fADCF(NULL), + fMCMT(NULL), + fTrackletArray(NULL), + fZSMap(NULL), + fTrklBranchName("mcmtrklbranch"), + fFeeParam(NULL), + fTrapConfig(NULL), + fDigitsManager(NULL), + fPedAcc(NULL), + fGainCounterA(NULL), + fGainCounterB(NULL), + fTailAmplLong(NULL), + fTailAmplShort(NULL), + fNHits(0), + fFitReg(NULL), + fDebugStream(0x0) { // // AliTRDmcmSim default constructor // By default, nothing is initialized. // It is necessary to issue Init before use. + + for (Int_t iDict = 0; iDict < 3; iDict++) + fDict[iDict] = 0x0; + + fFitPtr[0] = 0; + fFitPtr[1] = 0; + fFitPtr[2] = 0; + fFitPtr[3] = 0; } -AliTRDmcmSim::~AliTRDmcmSim() +AliTRDmcmSim::~AliTRDmcmSim() { // // AliTRDmcmSim destructor // if(fInitialized) { - for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { - delete [] fADCR[iadc]; - delete [] fADCF[iadc]; - delete [] fZSM [iadc]; + for( Int_t iAdc = 0 ; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++ ) { + delete [] fADCR[iAdc]; + delete [] fADCF[iAdc]; } delete [] fADCR; delete [] fADCF; - delete [] fZSM; - delete [] fZSM1Dim; + delete [] fZSMap; delete [] fMCMT; - + delete [] fPedAcc; delete [] fGainCounterA; delete [] fGainCounterB; delete [] fTailAmplLong; delete [] fTailAmplShort; delete [] fFitReg; - + fTrackletArray->Delete(); delete fTrackletArray; - delete fGeo; } } -void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ ) +void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ ) { // - // Initialize the class with new geometry information - // fADC array will be reused with filled by zero + // Initialize the class with new MCM position information + // memory is allocated in the first initialization // - + if (!fInitialized) { fFeeParam = AliTRDfeeParam::Instance(); - fTrapConfig = AliTRDtrapConfig::Instance(); - fSimParam = AliTRDSimParam::Instance(); - fCommonParam = AliTRDCommonParam::Instance(); - fCal = AliTRDcalibDB::Instance(); - fGeo = new AliTRDgeometry(); + fTrapConfig = AliTRDcalibDB::Instance()->GetTrapConfig(); } fDetector = det; fRobPos = robPos; fMcmPos = mcmPos; - fNADC = fFeeParam->GetNadcMcm(); - fNTimeBin = fCal->GetNumberOfTimeBins(); fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos ); - fMaxTracklets = fFeeParam->GetMaxNrOfTracklets(); - + if (!fInitialized) { - fADCR = new Int_t *[fNADC]; - fADCF = new Int_t *[fNADC]; - fZSM = new Int_t *[fNADC]; - fZSM1Dim = new Int_t [fNADC]; - fGainCounterA = new UInt_t[fNADC]; - fGainCounterB = new UInt_t[fNADC]; - for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { - fADCR[iadc] = new Int_t[fNTimeBin]; - fADCF[iadc] = new Int_t[fNTimeBin]; - fZSM [iadc] = new Int_t[fNTimeBin]; + fADCR = new Int_t *[AliTRDfeeParam::GetNadcMcm()]; + fADCF = new Int_t *[AliTRDfeeParam::GetNadcMcm()]; + fZSMap = new Int_t [AliTRDfeeParam::GetNadcMcm()]; + fGainCounterA = new UInt_t[AliTRDfeeParam::GetNadcMcm()]; + fGainCounterB = new UInt_t[AliTRDfeeParam::GetNadcMcm()]; + fNTimeBin = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kC13CPUA, fDetector, fRobPos, fMcmPos); + for( Int_t iAdc = 0 ; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++ ) { + fADCR[iAdc] = new Int_t[fNTimeBin]; + fADCF[iAdc] = new Int_t[fNTimeBin]; } - + // filter registers - fPedAcc = new UInt_t[fNADC]; // accumulator for pedestal filter - fTailAmplLong = new UShort_t[fNADC]; - fTailAmplShort = new UShort_t[fNADC]; - + fPedAcc = new UInt_t[AliTRDfeeParam::GetNadcMcm()]; // accumulator for pedestal filter + fTailAmplLong = new UShort_t[AliTRDfeeParam::GetNadcMcm()]; + fTailAmplShort = new UShort_t[AliTRDfeeParam::GetNadcMcm()]; + // tracklet calculation - fFitReg = new FitReg_t[fNADC]; - fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets); - - fMCMT = new UInt_t[fMaxTracklets]; + fFitReg = new FitReg_t[AliTRDfeeParam::GetNadcMcm()]; + fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fgkMaxTracklets); + + fMCMT = new UInt_t[fgkMaxTracklets]; } fInitialized = kTRUE; @@ -185,29 +181,53 @@ void AliTRDmcmSim::Reset() // Resets the data values and internal filter registers // by re-initialising them - for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { + if( !CheckInitialized() ) + return; + + for( Int_t iAdc = 0 ; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++ ) { for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { - fADCR[iadc][it] = 0; - fADCF[iadc][it] = 0; - fZSM [iadc][it] = 1; // Default unread = 1 + fADCR[iAdc][it] = 0; + fADCF[iAdc][it] = 0; } - fZSM1Dim[iadc] = 1; // Default unread = 1 - fGainCounterA[iadc] = 0; - fGainCounterB[iadc] = 0; + fZSMap[iAdc] = -1; // Default unread, low active bit mask + fGainCounterA[iAdc] = 0; + fGainCounterB[iAdc] = 0; } - - for(Int_t i = 0; i < fMaxTracklets; i++) { + + for(Int_t i = 0; i < fgkMaxTracklets; i++) { fMCMT[i] = 0; } - + + for (Int_t iDict = 0; iDict < 3; iDict++) + fDict[iDict] = 0x0; + FilterPedestalInit(); FilterGainInit(); - FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active + FilterTailInit(); +} + +void AliTRDmcmSim::SetNTimebins(Int_t ntimebins) +{ + // Reallocate memory if a change in the number of timebins + // is needed (should not be the case for real data) + + if( !CheckInitialized() ) + return; + + fNTimeBin = ntimebins; + for( Int_t iAdc = 0 ; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++ ) { + delete [] fADCR[iAdc]; + delete [] fADCF[iAdc]; + fADCR[iAdc] = new Int_t[fNTimeBin]; + fADCF[iAdc] = new Int_t[fNTimeBin]; + } } -Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm) +Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm) { - // loads the ADC data as obtained from the digitsManager for the specified MCM + // loads the ADC data as obtained from the digitsManager for the specified MCM. + // This method is meant for rare execution, e.g. in the visualization. When called + // frequently use SetData(...) instead. Init(det, rob, mcm); @@ -222,6 +242,7 @@ Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob return kFALSE; } + Bool_t retval = kTRUE; trdLoader->LoadDigits(); fDigitsManager = 0x0; AliTRDdigitsManager *digMgr = new AliTRDdigitsManager(); @@ -229,91 +250,121 @@ Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob digMgr->CreateArrays(); digMgr->ReadDigits(trdLoader->TreeD()); AliTRDarrayADC *digits = (AliTRDarrayADC*) digMgr->GetDigits(det); - if (!digits->HasData()) - return kFALSE; - digits->Expand(); - - Int_t padrow = fFeeParam->GetPadRowFromMCM(rob, mcm); - Int_t padcol = 0; - for (Int_t ch = 0; ch < fNADC; ch++) { - padcol = GetCol(ch); - for (Int_t tb = 0; tb < fNTimeBin; tb++) { - if (padcol < 0) { - fADCR[ch][tb] = 0; - fADCF[ch][tb] = 0; - } - else { - if (digits->GetData(padrow,padcol, tb) < 0) { - fADCR[ch][tb] = 0; - fADCF[ch][tb] = 0; - } - else { - fADCR[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits; - fADCF[ch][tb] = digits->GetData(padrow, padcol, tb) << fgkAddDigits; - } - } + if (digits->HasData()) { + digits->Expand(); + + if (fNTimeBin != digits->GetNtime()) { + AliWarning(Form("Changing no. of timebins from %i to %i", fNTimeBin, digits->GetNtime())); + SetNTimebins(digits->GetNtime()); } + + SetData(digits); } + else + retval = kFALSE; + delete digMgr; - return kTRUE; + return retval; } void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail) { - // This function can be used to test the filters. + // This function can be used to test the filters. // It feeds nsamples of ADC values with a gaussian distribution specified by mean and sigma. // The filter chain implemented here consists of: // Pedestal -> Gain -> Tail - // With inputGain and inputTail the input to the gain and tail filter, respectively, - // can be chosen where + // With inputGain and inputTail the input to the gain and tail filter, respectively, + // can be chosen where // 0: noise input // 1: pedestal output // 2: gain output - // The input has to be chosen from a stage before. - // The filter behaviour is controlled by the TRAP parameters from AliTRDtrapConfig in the + // The input has to be chosen from a stage before. + // The filter behaviour is controlled by the TRAP parameters from AliTRDtrapConfig in the // same way as in normal simulation. // The functions produces four histograms with the values at the different stages. + if( !CheckInitialized() ) + return; + + TString nameInputGain; + TString nameInputTail; + + switch (inputGain) { + case 0: + nameInputGain = "Noise"; + break; + + case 1: + nameInputGain = "Pedestal"; + break; + + default: + AliError("Undefined input to tail cancellation filter"); + return; + } + + switch (inputTail) { + case 0: + nameInputTail = "Noise"; + break; + + case 1: + nameInputTail = "Pedestal"; + break; + + case 2: + nameInputTail = "Gain"; + break; + + default: + AliError("Undefined input to tail cancellation filter"); + return; + } + TH1F *h = new TH1F("noise", "Gaussian Noise;sample;ADC count", nsamples, 0, nsamples); - TH1F *hfp = new TH1F("pedf", "Noise #rightarrow Pedestal filter;sample;ADC count", nsamples, 0, nsamples); - TH1F *hfg = new TH1F("pedg", "Pedestal #rightarrow Gain;sample;ADC count", nsamples, 0, nsamples); - TH1F *hft = new TH1F("pedt", "Gain #rightarrow Tail;sample;ADC count", nsamples, 0, nsamples); + TH1F *hfp = new TH1F("ped", "Noise #rightarrow Pedestal filter;sample;ADC count", nsamples, 0, nsamples); + TH1F *hfg = new TH1F("gain", + (nameInputGain + "#rightarrow Gain;sample;ADC count").Data(), + nsamples, 0, nsamples); + TH1F *hft = new TH1F("tail", + (nameInputTail + "#rightarrow Tail;sample;ADC count").Data(), + nsamples, 0, nsamples); h->SetStats(kFALSE); hfp->SetStats(kFALSE); hfg->SetStats(kFALSE); hft->SetStats(kFALSE); - + Int_t value; // ADC count with noise (10 bit) Int_t valuep; // pedestal filter output (12 bit) Int_t valueg; // gain filter output (12 bit) Int_t valuet; // tail filter value (12 bit) - + for (Int_t i = 0; i < nsamples; i++) { - value = (Int_t) gRandom->Gaus(mean, sigma); // generate noise with gaussian distribution + value = (Int_t) gRandom->Gaus(mean, sigma); // generate noise with gaussian distribution h->SetBinContent(i, value); valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2); - + if (inputGain == 0) valueg = FilterGainNextSample(1, ((Int_t) value) << 2); - else - valueg = FilterGainNextSample(1, valuep); - + else + valueg = FilterGainNextSample(1, valuep); + if (inputTail == 0) valuet = FilterTailNextSample(1, ((Int_t) value) << 2); else if (inputTail == 1) - valuet = FilterTailNextSample(1, valuep); + valuet = FilterTailNextSample(1, valuep); else - valuet = FilterTailNextSample(1, valueg); + valuet = FilterTailNextSample(1, valueg); hfp->SetBinContent(i, valuep >> 2); hfg->SetBinContent(i, valueg >> 2); hft->SetBinContent(i, valuet >> 2); } - TCanvas *c = new TCanvas; + TCanvas *c = new TCanvas; c->Divide(2,2); c->cd(1); h->Draw(); @@ -325,51 +376,38 @@ void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inpu hft->Draw(); } -Bool_t AliTRDmcmSim::CheckInitialized() +Bool_t AliTRDmcmSim::CheckInitialized() const { // // Check whether object is initialized // - if( ! fInitialized ) { - AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called.")); - } + if( ! fInitialized ) + AliError(Form ("AliTRDmcmSim is not initialized but function other than Init() is called.")); + return fInitialized; } void AliTRDmcmSim::Print(Option_t* const option) const { // Prints the data stored and/or calculated for this MCM. - // The output is controlled by option which can be a sequence of any of + // The output is controlled by option which can be a sequence of any of // the following characters: // R - prints raw ADC data - // F - prints filtered data + // F - prints filtered data // H - prints detected hits // T - prints found tracklets - // The later stages are only useful when the corresponding calculations + // The later stages are only meaningful after the corresponding calculations // have been performed. + if ( !CheckInitialized() ) + return; + printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector); TString opt = option; - if (opt.Contains("R")) { - printf("Raw ADC data (10 bit):\n"); - for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) { - printf("%5i", fADCR[iChannel][iTimeBin] >> fgkAddDigits); - } - printf("\n"); - } - } - - if (opt.Contains("F")) { - printf("Filtered data (12 bit):\n"); - for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) { - printf("%5i", fADCF[iChannel][iTimeBin]); - } - printf("\n"); - } + if (opt.Contains("R") || opt.Contains("F")) { + std::cout << *this; } if (opt.Contains("H")) { @@ -388,36 +426,39 @@ void AliTRDmcmSim::Print(Option_t* const option) const } } -void AliTRDmcmSim::Draw(Option_t* const option) +void AliTRDmcmSim::Draw(Option_t* const option) { // Plots the data stored in a 2-dim. timebin vs. ADC channel plot. - // The option selects what data is plotted and can be a sequence of + // The option selects what data is plotted and can be a sequence of // the following characters: // R - plot raw data (default) // F - plot filtered data (meaningless if R is specified) // In addition to the ADC values: - // H - plot hits + // H - plot hits // T - plot tracklets + if( !CheckInitialized() ) + return; + TString opt = option; TH2F *hist = new TH2F("mcmdata", Form("Data of MCM %i on ROB %i in detector %i", \ fMcmPos, fRobPos, fDetector), \ - fNADC, -0.5, fNADC-.5, fNTimeBin, -.5, fNTimeBin-.5); + AliTRDfeeParam::GetNadcMcm(), -0.5, AliTRDfeeParam::GetNadcMcm()-.5, fNTimeBin, -.5, fNTimeBin-.5); hist->GetXaxis()->SetTitle("ADC Channel"); hist->GetYaxis()->SetTitle("Timebin"); hist->SetStats(kFALSE); if (opt.Contains("R")) { for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { hist->SetBinContent(iAdc+1, iTimeBin+1, fADCR[iAdc][iTimeBin] >> fgkAddDigits); } } } else { for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits); } } @@ -427,8 +468,8 @@ void AliTRDmcmSim::Draw(Option_t* const option) if (opt.Contains("H")) { TGraph *grHits = new TGraph(); for (Int_t iHit = 0; iHit < fNHits; iHit++) { - grHits->SetPoint(iHit, - fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256., + grHits->SetPoint(iHit, + fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256., fHits[iHit].fTimebin); } grHits->Draw("*"); @@ -437,13 +478,21 @@ void AliTRDmcmSim::Draw(Option_t* const option) if (opt.Contains("T")) { TLine *trklLines = new TLine[4]; for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) { - AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector); AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl]; - Float_t offset = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19)) + 19 * pp->GetWidthIPad(); - trklLines[iTrkl].SetX1((offset - trkl->GetY()) / pp->GetWidthIPad()); - trklLines[iTrkl].SetY1(0); - trklLines[iTrkl].SetX2((offset - (trkl->GetY() + ((Float_t) trkl->GetdY())*140e-4)) / pp->GetWidthIPad()); - trklLines[iTrkl].SetY2(fNTimeBin - 1); + Float_t padWidth = 0.635 + 0.03 * (fDetector % 6); + Float_t offset = padWidth/256. * ((((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - ((18*4*2 - 18*2 - 3) << 7)); // revert adding offset in FitTracklet + Int_t ndrift = fTrapConfig->GetDmemUnsigned(fgkDmemAddrNdrift, fDetector, fRobPos, fMcmPos) >> 5; + Float_t slope = 0; + if (ndrift) + slope = trkl->GetdY() * 140e-4 / ndrift; + + Int_t t0 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS, fDetector, fRobPos, fMcmPos); + Int_t t1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE, fDetector, fRobPos, fMcmPos); + + trklLines[iTrkl].SetX1((offset - (trkl->GetY() - slope * t0)) / padWidth); // ??? sign? + trklLines[iTrkl].SetY1(t0); + trklLines[iTrkl].SetX2((offset - (trkl->GetY() - slope * t1)) / padWidth); // ??? sign? + trklLines[iTrkl].SetY2(t1); trklLines[iTrkl].SetLineColor(2); trklLines[iTrkl].SetLineWidth(2); printf("Tracklet %i: y = %f, dy = %f, offset = %f\n", iTrkl, trkl->GetY(), (trkl->GetdY() * 140e-4), offset); @@ -452,7 +501,7 @@ void AliTRDmcmSim::Draw(Option_t* const option) } } -void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc ) +void AliTRDmcmSim::SetData( Int_t adc, const Int_t* const data ) { // // Store ADC data into array of raw data @@ -460,18 +509,18 @@ void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc ) if( !CheckInitialized() ) return; - if( iadc < 0 || iadc >= fNADC ) { - //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1)); + if( adc < 0 || adc >= AliTRDfeeParam::GetNadcMcm() ) { + AliError(Form ("Error: ADC %i is out of range (0 .. %d).", adc, AliTRDfeeParam::GetNadcMcm()-1)); return; } - for( int it = 0 ; it < fNTimeBin ; it++ ) { - fADCR[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits; - fADCF[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits; + for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { + fADCR[adc][it] = (Int_t) (data[it]) << fgkAddDigits; + fADCF[adc][it] = (Int_t) (data[it]) << fgkAddDigits; } } -void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc ) +void AliTRDmcmSim::SetData( Int_t adc, Int_t it, Int_t data ) { // // Store ADC data into array of raw data @@ -479,149 +528,246 @@ void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc ) if( !CheckInitialized() ) return; - if( iadc < 0 || iadc >= fNADC ) { - //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1)); + if( adc < 0 || adc >= AliTRDfeeParam::GetNadcMcm() ) { + AliError(Form ("Error: ADC %i is out of range (0 .. %d).", adc, AliTRDfeeParam::GetNadcMcm()-1)); return; } - fADCR[iadc][it] = adc << fgkAddDigits; - fADCF[iadc][it] = adc << fgkAddDigits; + fADCR[adc][it] = data << fgkAddDigits; + fADCF[adc][it] = data << fgkAddDigits; } -void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager) +void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager * const digitsManager) { // Set the ADC data from an AliTRDarrayADC - if (!fInitialized) { - AliError("Called uninitialized! Nothing done!"); + if( !CheckInitialized() ) return; - } fDigitsManager = digitsManager; + if (fDigitsManager) { + for (Int_t iDict = 0; iDict < 3; iDict++) { + AliTRDarrayDictionary *newDict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict); + if (fDict[iDict] != 0x0 && newDict != 0x0) { - Int_t firstAdc = 0; - Int_t lastAdc = fNADC-1; + if (fDict[iDict] == newDict) + continue; - while (GetCol(firstAdc) < 0) { - for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - fADCR[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; - fADCF[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; + fDict[iDict] = newDict; + if(fDict[iDict]->GetDim() != 0) + fDict[iDict]->Expand(); + } + else { + fDict[iDict] = newDict; + if (fDict[iDict] && (fDict[iDict]->GetDim() != 0) ) + fDict[iDict]->Expand(); + } + + // If there is no data, set dictionary to zero to avoid crashes + if (fDict[iDict]->GetDim() == 0) { + // AliError(Form("Dictionary %i of det. %i has dim. 0", iDict, fDetector)); + fDict[iDict] = 0x0; + } } - firstAdc++; } - while (GetCol(lastAdc) < 0) { - for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - fADCR[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; - fADCF[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; + if (fNTimeBin != adcArray->GetNtime()) + SetNTimebins(adcArray->GetNtime()); + + Int_t offset = (fMcmPos % 4 + 1) * 21 + (fRobPos % 2) * 84 - 1; + + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { + Int_t value = adcArray->GetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin); + // treat 0 as suppressed, + // this is not correct but reported like that from arrayADC + if (value <= 0 || (offset - iAdc < 1) || (offset - iAdc > 165)) { + fADCR[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits); + fADCF[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits); + } + else { + fZSMap[iAdc] = 0; + fADCR[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits); + fADCF[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits); + } + } + } +} + +void AliTRDmcmSim::SetDataByPad(const AliTRDarrayADC* const adcArray, AliTRDdigitsManager * const digitsManager) +{ + // Set the ADC data from an AliTRDarrayADC + // (by pad, to be used during initial reading in simulation) + + if( !CheckInitialized() ) + return; + + fDigitsManager = digitsManager; + if (fDigitsManager) { + for (Int_t iDict = 0; iDict < 3; iDict++) { + AliTRDarrayDictionary *newDict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict); + if (fDict[iDict] != 0x0 && newDict != 0x0) { + + if (fDict[iDict] == newDict) + continue; + + fDict[iDict] = newDict; + fDict[iDict]->Expand(); + } + else { + fDict[iDict] = newDict; + if (fDict[iDict]) + fDict[iDict]->Expand(); + } + + // If there is no data, set dictionary to zero to avoid crashes + if (fDict[iDict]->GetDim() == 0) { + AliError(Form("Dictionary %i of det. %i has dim. 0", iDict, fDetector)); + fDict[iDict] = 0x0; + } } - lastAdc--; } + if (fNTimeBin != adcArray->GetNtime()) + SetNTimebins(adcArray->GetNtime()); + + Int_t offset = (fMcmPos % 4 + 1) * 18 + (fRobPos % 2) * 72 + 1; + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) { - Int_t value = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin); - if (value < 0) { - fADCR[iAdc][iTimeBin] = 0; - fADCF[iAdc][iTimeBin] = 0; + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { + Int_t value = -1; + Int_t pad = offset - iAdc; + if (pad > -1 && pad < 144) + value = adcArray->GetData(GetRow(), offset - iAdc, iTimeBin); + // Int_t value = adcArray->GetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin); + if (value < 0 || (offset - iAdc < 1) || (offset - iAdc > 165)) { + fADCR[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits); + fADCF[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits); } else { - fADCR[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits; - fADCF[iAdc][iTimeBin] = adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin) << fgkAddDigits; + fZSMap[iAdc] = 0; + fADCR[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits); + fADCF[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits); } } } } -void AliTRDmcmSim::SetDataPedestal( Int_t iadc ) +void AliTRDmcmSim::SetDataPedestal( Int_t adc ) { // // Store ADC data into array of raw data // - if( !CheckInitialized() ) return; + if( !CheckInitialized() ) + return; - if( iadc < 0 || iadc >= fNADC ) { - //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1)); + if( adc < 0 || adc >= AliTRDfeeParam::GetNadcMcm() ) { return; } for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { - fADCR[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits; - fADCF[iadc][it] = fSimParam->GetADCbaseline() << fgkAddDigits; + fADCR[adc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits); + fADCF[adc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits); } } -Int_t AliTRDmcmSim::GetCol( Int_t iadc ) +Bool_t AliTRDmcmSim::GetHit(Int_t index, Int_t &channel, Int_t &timebin, Int_t &qtot, Int_t &ypos, Float_t &y, Int_t &label) const +{ + // retrieve the MC hit information (not available in TRAP hardware) + + if (index < 0 || index >= fNHits) + return kFALSE; + + channel = fHits[index].fChannel; + timebin = fHits[index].fTimebin; + qtot = fHits[index].fQtot; + ypos = fHits[index].fYpos; + y = (Float_t) ((((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - ((18*4*2 - 18*2 - 1) << 7) - + (channel << 8) - ypos) + * (0.635 + 0.03 * (fDetector % 6)) + / 256.0; + label = fHits[index].fLabel[0]; + + return kTRUE; +} + +Int_t AliTRDmcmSim::GetCol( Int_t adc ) { // // Return column id of the pad for the given ADC channel // - if( !CheckInitialized() ) + if( !CheckInitialized() ) return -1; - Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc); - if (col < 0 || col >= fFeeParam->GetNcol()) + Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adc); + if (col < 0 || col >= fFeeParam->GetNcol()) return -1; - else + else return col; } -Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv) +Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t bufSize, UInt_t iEv) const { // // Produce raw data stream from this MCM and put in buf - // Returns number of words filled, or negative value + // Returns number of words filled, or negative value // with -1 * number of overflowed words // + if( !CheckInitialized() ) + return 0; + UInt_t x; + UInt_t mcmHeader = 0; + UInt_t adcMask = 0; Int_t nw = 0; // Number of written words Int_t of = 0; // Number of overflowed words Int_t rawVer = fFeeParam->GetRAWversion(); Int_t **adc; Int_t nActiveADC = 0; // number of activated ADC bits in a word - if( !CheckInitialized() ) return 0; + if( !CheckInitialized() ) + return 0; - if( fFeeParam->GetRAWstoreRaw() ) { + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF, fDetector, fRobPos, fMcmPos) != 0) // store unfiltered data adc = fADCR; - } else { + else adc = fADCF; - } - - // Produce MCM header - x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC; - - if (nw < maxSize) { - buf[nw++] = x; - //printf("\nMCM header: %X ",x); - } - else { - of++; - } // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100 // n : unused , c : ADC count, m : selected ADCs - if( rawVer >= 3 ) { - x = 0; - for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) { - if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed - x = x | (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc - nActiveADC++; // number of 1 in mmm....m + if( rawVer >= 3 && + (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kC15CPUA, fDetector, fRobPos, fMcmPos) & (1 << 13))) { // check for zs flag in TRAP configuration + for( Int_t iAdc = 0 ; iAdc < AliTRDfeeParam::GetNadcMcm() ; iAdc++ ) { + if( ~fZSMap[iAdc] != 0 ) { // 0 means not suppressed + adcMask |= (1 << (iAdc+4) ); // last 4 digit reserved for 1100=0xc + nActiveADC++; // number of 1 in mmm....m } } - x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC; // nn = 01, ccccc are inverted, 0xc=1100 - //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x ); - if (nw < maxSize) { - buf[nw++] = x; - //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC); - } - else { + if ((nActiveADC == 0) && + (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kC15CPUA, fDetector, fRobPos, fMcmPos) & (1 << 8))) // check for DEH flag in TRAP configuration + return 0; + + // assemble adc mask word + adcMask |= (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC; // nn = 01, ccccc are inverted, 0xc=1100 + } + + // MCM header + mcmHeader = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC; + if (nw < bufSize) + buf[nw++] = mcmHeader; + else + of++; + + // ADC mask + if( adcMask != 0 ) { + if (nw < bufSize) + buf[nw++] = adcMask; + else of++; - } } // Produce ADC data. 3 timebins are packed into one 32 bits word @@ -630,16 +776,15 @@ Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv) UInt_t aa=0, a1=0, a2=0, a3=0; for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) { - if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed + if( rawVer>= 3 && ~fZSMap[iAdc] == 0 ) continue; // Zero Suppression, 0 means not suppressed aa = !(iAdc & 1) + 2; for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) { a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] >> fgkAddDigits : 0; a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0; a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0; x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa; - if (nw < maxSize) { + if (nw < bufSize) { buf[nw++] = x; - //printf("%08X ",x); } else { of++; @@ -650,29 +795,30 @@ Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv) if( of != 0 ) return -of; else return nw; } -Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize ) +Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t bufSize ) { // // Produce tracklet data stream from this MCM and put in buf - // Returns number of words filled, or negative value + // Returns number of words filled, or negative value // with -1 * number of overflowed words // + if( !CheckInitialized() ) + return 0; + Int_t nw = 0; // Number of written words Int_t of = 0; // Number of overflowed words - - if( !CheckInitialized() ) return 0; - // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM + // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM // fMCMT is filled continuously until no more tracklet words available for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) { - if (nw < maxSize) + if (nw < bufSize) buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord(); - else + else of++; } - + if( of != 0 ) return -of; else return nw; } @@ -681,22 +827,20 @@ void AliTRDmcmSim::Filter() // // Filter the raw ADC values. The active filter stages and their // parameters are taken from AliTRDtrapConfig. - // The raw data is stored separate from the filtered data. Thus, - // it is possible to run the filters on a set of raw values + // The raw data is stored separate from the filtered data. Thus, + // it is possible to run the filters on a set of raw values // sequentially for parameter tuning. // - if( !CheckInitialized() ) { - AliError("got called before initialization! Nothing done!"); + if( !CheckInitialized() ) return; - } // Apply filters sequentially. Bypass is handled by filters - // since counters and internal registers may be updated even + // since counters and internal registers may be updated even // if the filter is bypassed. - // The first filter takes the data from fADCR and - // outputs to fADCF. - + // The first filter takes the data from fADCR and + // outputs to fADCF. + // Non-linearity filter not implemented. FilterPedestal(); FilterGain(); @@ -704,54 +848,51 @@ void AliTRDmcmSim::Filter() // Crosstalk filter not implemented. } -void AliTRDmcmSim::FilterPedestalInit() +void AliTRDmcmSim::FilterPedestalInit(Int_t baseline) { - // Initializes the pedestal filter assuming that the input has + // Initializes the pedestal filter assuming that the input has // been constant for a long time (compared to the time constant). -// UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output - UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest - UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to take shifts from? + UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC, fDetector, fRobPos, fMcmPos); // 0..3, 0 - fastest, 3 - slowest - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) - fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output - UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest - UShort_t fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 the bypass, active low - UShort_t shifts[4] = {11, 14, 17, 21}; //??? where to come from + UShort_t fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP, fDetector, fRobPos, fMcmPos); // 0..511 -> 0..127.75, pedestal at the output + UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC, fDetector, fRobPos, fMcmPos); // 0..3, 0 - fastest, 3 - slowest + UShort_t fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY, fDetector, fRobPos, fMcmPos); // 0..1 bypass, active low UShort_t accumulatorShifted; Int_t correction; UShort_t inpAdd; - - inpAdd = value + fpnp; - if (fpby == 0) //??? before or after update of accumulator - return value; + inpAdd = value + fpnp; - accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF; // 10 bits + accumulatorShifted = (fPedAcc[adc] >> fgkFPshifts[fptc]) & 0x3FF; // 10 bits if (timebin == 0) // the accumulator is disabled in the drift time { correction = (value & 0x3FF) - accumulatorShifted; fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF; // 31 bits } - + + if (fpby == 0) + return value; + if (inpAdd <= accumulatorShifted) return 0; else { inpAdd = inpAdd - accumulatorShifted; - if (inpAdd > 0xFFF) + if (inpAdd > 0xFFF) return 0xFFF; - else + else return inpAdd; } } @@ -761,14 +902,14 @@ void AliTRDmcmSim::FilterPedestal() // // Apply pedestal filter // - // As the first filter in the chain it reads data from fADCR - // and outputs to fADCF. - // It has only an effect if previous samples have been fed to - // find the pedestal. Currently, the simulation assumes that + // As the first filter in the chain it reads data from fADCR + // and outputs to fADCF. + // It has only an effect if previous samples have been fed to + // find the pedestal. Currently, the simulation assumes that // the input has been stable for a sufficiently long time. for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]); } } @@ -776,11 +917,11 @@ void AliTRDmcmSim::FilterPedestal() void AliTRDmcmSim::FilterGainInit() { - // Initializes the gain filter. In this case, only threshold + // Initializes the gain filter. In this case, only threshold // counters are reset. - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) { - // these are counters which in hardware continue + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { + // these are counters which in hardware continue // until maximum or reset fGainCounterA[iAdc] = 0; fGainCounterB[iAdc] = 0; @@ -791,42 +932,48 @@ UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value) { // Apply the gain filter to the given value. // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX - // The output depends on the internal registers and, thus, the + // The output depends on the internal registers and, thus, the // history of the filter. - UShort_t fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low - UShort_t fgf = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff); - UShort_t fga = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40; - UShort_t fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20; - UShort_t fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060; + UShort_t fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY, fDetector, fRobPos, fMcmPos); // bypass, active low + UShort_t fgf = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc), fDetector, fRobPos, fMcmPos); // 0x700 + (0 & 0x1ff); + UShort_t fga = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc), fDetector, fRobPos, fMcmPos); // 40; + UShort_t fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA, fDetector, fRobPos, fMcmPos); // 20; + UShort_t fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB, fDetector, fRobPos, fMcmPos); // 2060; - UInt_t tmp; + UInt_t fgfExtended = 0x700 + fgf; // The corr factor which is finally applied has to be extended by 0x700 (hex) or 0.875 (dec) + // because fgf=0 correspons to 0.875 and fgf=511 correspons to 1.125 - 2^(-11) + // (see TRAP User Manual for details) - value &= 0xFFF; - tmp = (value * fgf) >> 11; - if (tmp > 0xFFF) tmp = 0xFFF; + UInt_t corr; // corrected value - if (fgby == 1) - value = AddUintClipping(tmp, fga, 12); + value &= 0xFFF; + corr = (value * fgfExtended) >> 11; + corr = corr > 0xfff ? 0xfff : corr; + corr = AddUintClipping(corr, fga, 12); - // Update threshold counters + // Update threshold counters // not really useful as they are cleared with every new event - if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF)) + if (!((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))) + // stop when full { - if (value >= fgtb) + if (corr >= fgtb) fGainCounterB[adc]++; - else if (value >= fgta) + else if (corr >= fgta) fGainCounterA[adc]++; } - return value; + if (fgby == 1) + return corr; + else + return value; } void AliTRDmcmSim::FilterGain() { // Read data from fADCF and apply gain filter. - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]); } @@ -835,14 +982,14 @@ void AliTRDmcmSim::FilterGain() void AliTRDmcmSim::FilterTailInit(Int_t baseline) { - // Initializes the tail filter assuming that the input has - // been at the baseline value (configured by FTFP) for a + // Initializes the tail filter assuming that the input has + // been at the baseline value (configured by FTFP) for a // sufficiently long time. // exponents and weight calculated from configuration - UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component - UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier - UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier + UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL, fDetector, fRobPos, fMcmPos); // the weight of the long component + UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL, fDetector, fRobPos, fMcmPos) & 0x1FF); // the multiplier + UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS, fDetector, fRobPos, fMcmPos) & 0x1FF); // the multiplier Float_t lambdaL = lambdaLong * 1.0 / (1 << 11); Float_t lambdaS = lambdaShort * 1.0 / (1 << 11); @@ -854,13 +1001,22 @@ void AliTRDmcmSim::FilterTailInit(Int_t baseline) Float_t kt, ql, qs; UShort_t aout; - - kt = kdc * baseline; - aout = baseline - (UShort_t) kt; + + if (baseline < 0) + baseline = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP, fDetector, fRobPos, fMcmPos); + ql = lambdaL * (1 - lambdaS) * alphaL; qs = lambdaS * (1 - lambdaL) * (1 - alphaL); - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { + Int_t value = baseline & 0xFFF; + Int_t corr = (value * fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + iAdc), fDetector, fRobPos, fMcmPos)) >> 11; + corr = corr > 0xfff ? 0xfff : corr; + corr = AddUintClipping(corr, fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + iAdc), fDetector, fRobPos, fMcmPos), 12); + + kt = kdc * baseline; + aout = baseline - (UShort_t) kt; + fTailAmplLong[iAdc] = (UShort_t) (aout * ql / (ql + qs)); fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs)); } @@ -868,67 +1024,58 @@ void AliTRDmcmSim::FilterTailInit(Int_t baseline) UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value) { - // Returns the output of the tail filter for the given input value. - // The output depends on the internal registers and, thus, the + // Returns the output of the tail filter for the given input value. + // The output depends on the internal registers and, thus, the // history of the filter. // exponents and weight calculated from configuration - UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component - UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier - UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier - - Float_t lambdaL = lambdaLong * 1.0 / (1 << 11); - Float_t lambdaS = lambdaShort * 1.0 / (1 << 11); - Float_t alphaL = alphaLong * 1.0 / (1 << 11); - Float_t qup, qdn; - qup = (1 - lambdaL) * (1 - lambdaS); - qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL); -// Float_t kdc = qup/qdn; + UShort_t alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL, fDetector, fRobPos, fMcmPos); // the weight of the long component + UShort_t lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL, fDetector, fRobPos, fMcmPos) & 0x1FF); // the multiplier of the long component + UShort_t lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS, fDetector, fRobPos, fMcmPos) & 0x1FF); // the multiplier of the short component - UInt_t aDiff; - UInt_t alInpv; + // intermediate signals + UInt_t aDiff; + UInt_t alInpv; UShort_t aQ; - UInt_t tmp; - + UInt_t tmp; + UShort_t inpVolt = value & 0xFFF; // 12 bits - - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low + + // add the present generator outputs + aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12); + + // calculate the difference between the input and the generated signal + if (inpVolt > aQ) + aDiff = inpVolt - aQ; + else + aDiff = 0; + + // the inputs to the two generators, weighted + alInpv = (aDiff * alphaLong) >> 11; + + // the new values of the registers, used next time + // long component + tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12); + tmp = (tmp * lambdaLong) >> 11; + fTailAmplLong[adc] = tmp & 0xFFF; + // short component + tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12); + tmp = (tmp * lambdaShort) >> 11; + fTailAmplShort[adc] = tmp & 0xFFF; + + // the output of the filter + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY, fDetector, fRobPos, fMcmPos) == 0) // bypass mode, active low return value; else - { - // add the present generator outputs - aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12); - - // calculate the difference between the input the generated signal - if (inpVolt > aQ) - aDiff = inpVolt - aQ; - else - aDiff = 0; - - // the inputs to the two generators, weighted - alInpv = (aDiff * alphaLong) >> 11; - - // the new values of the registers, used next time - // long component - tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12); - tmp = (tmp * lambdaLong) >> 11; - fTailAmplLong[adc] = tmp & 0xFFF; - // short component - tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12); - tmp = (tmp * lambdaShort) >> 11; - fTailAmplShort[adc] = tmp & 0xFFF; - - // the output of the filter return aDiff; - } } void AliTRDmcmSim::FilterTail() { - // Apply tail filter + // Apply tail cancellation filter to all data. for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]); } } @@ -938,146 +1085,108 @@ void AliTRDmcmSim::ZSMapping() { // // Zero Suppression Mapping implemented in TRAP chip + // only implemented for up to 30 timebins // // See detail TRAP manual "Data Indication" section: // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf // - //??? values should come from TRAPconfig - Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); // TRAP default = 0x4 (Tis=4) - Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); // TRAP default = 0x28 (Tit=40) - Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); // TRAP default = 0xf0 - // (lookup table accept (I2,I1,I0)=(111) - // or (110) or (101) or (100)) - Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity) - Int_t ep = 0; // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); //??? really subtracted here + if( !CheckInitialized() ) + return; + + Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS, fDetector, fRobPos, fMcmPos); + Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT, fDetector, fRobPos, fMcmPos); + Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL, fDetector, fRobPos, fMcmPos); + Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN, fDetector, fRobPos, fMcmPos); Int_t **adc = fADCF; - if( !CheckInitialized() ) { - AliError("got called uninitialized! Nothing done!"); - return; - } + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) + fZSMap[iAdc] = -1; for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { - for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) { - - // Get ADC data currently in filter buffer - Int_t ap = adc[iadc-1][it] - ep; // previous - Int_t ac = adc[iadc ][it] - ep; // current - Int_t an = adc[iadc+1][it] - ep; // next - - // evaluate three conditions - Int_t i0 = ( ac >= ap && ac >= an ) ? 0 : 1; // peak center detection - Int_t i1 = ( ap + ac + an > eBIT ) ? 0 : 1; // cluster - Int_t i2 = ( ac > eBIS ) ? 0 : 1; // absolute large peak - - Int_t i = i2 * 4 + i1 * 2 + i0; // Bit position in lookup table - Int_t d = (eBIL >> i) & 1; // Looking up (here d=0 means true - // and d=1 means false according to TRAP manual) - - fZSM[iadc][it] &= d; - if( eBIN == 0 ) { // turn on neighboring ADCs - fZSM[iadc-1][it] &= d; - fZSM[iadc+1][it] &= d; - } + Int_t iAdc; // current ADC channel + Int_t ap; + Int_t ac; + Int_t an; + Int_t mask; + Int_t supp; // suppression of the current channel (low active) + + // ----- first channel ----- + iAdc = 0; + + ap = 0; // previous + ac = adc[iAdc ][it]; // current + an = adc[iAdc+1][it]; // next + + mask = ( ac >= ap && ac >= an ) ? 0 : 0x1; // peak center detection + mask += ( ap + ac + an > eBIT ) ? 0 : 0x2; // cluster + mask += ( ac > eBIS ) ? 0 : 0x4; // absolute large peak + + supp = (eBIL >> mask) & 1; + + fZSMap[iAdc] &= ~((1-supp) << it); + if( eBIN == 0 ) { // neighbour sensitivity + fZSMap[iAdc+1] &= ~((1-supp) << it); } - } - // do 1 dim projection - for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { - for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { - fZSM1Dim[iadc] &= fZSM[iadc][it]; + // ----- last channel ----- + iAdc = AliTRDfeeParam::GetNadcMcm() - 1; + + ap = adc[iAdc-1][it]; // previous + ac = adc[iAdc ][it]; // current + an = 0; // next + + mask = ( ac >= ap && ac >= an ) ? 0 : 0x1; // peak center detection + mask += ( ap + ac + an > eBIT ) ? 0 : 0x2; // cluster + mask += ( ac > eBIS ) ? 0 : 0x4; // absolute large peak + + supp = (eBIL >> mask) & 1; + + fZSMap[iAdc] &= ~((1-supp) << it); + if( eBIN == 0 ) { // neighbour sensitivity + fZSMap[iAdc-1] &= ~((1-supp) << it); } - } -} -void AliTRDmcmSim::DumpData( char *f, char *target ) -{ - // - // Dump data stored (for debugging). - // target should contain one or multiple of the following characters - // R for raw data - // F for filtered data - // Z for zero suppression map - // S Raw dat astream - // other characters are simply ignored - // + // ----- middle channels ----- + for( iAdc = 1 ; iAdc < AliTRDfeeParam::GetNadcMcm()-1; iAdc++ ) { + ap = adc[iAdc-1][it]; // previous + ac = adc[iAdc ][it]; // current + an = adc[iAdc+1][it]; // next - UInt_t tempbuf[1024]; + mask = ( ac >= ap && ac >= an ) ? 0 : 0x1; // peak center detection + mask += ( ap + ac + an > eBIT ) ? 0 : 0x2; // cluster + mask += ( ac > eBIS ) ? 0 : 0x4; // absolute large peak - if( !CheckInitialized() ) return; + supp = (eBIL >> mask) & 1; - std::ofstream of( f, std::ios::out | std::ios::app ); - of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n", - fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector), - fGeo->GetSector(fDetector), fRobPos, fMcmPos ); - - for( int t=0 ; target[t] != 0 ; t++ ) { - switch( target[t] ) { - case 'R' : - case 'r' : - of << Form("fADCR (raw ADC data)\n"); - for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { - of << Form(" ADC %02d: ", iadc); - for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { - of << Form("% 4d", fADCR[iadc][it]); - } - of << Form("\n"); - } - break; - case 'F' : - case 'f' : - of << Form("fADCF (filtered ADC data)\n"); - for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { - of << Form(" ADC %02d: ", iadc); - for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { - of << Form("% 4d", fADCF[iadc][it]); - } - of << Form("\n"); - } - break; - case 'Z' : - case 'z' : - of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n"); - for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) { - of << Form(" ADC %02d: ", iadc); - if( fZSM1Dim[iadc] == 0 ) { of << " R " ; } else { of << " . "; } // R:read .:suppressed - for( Int_t it = 0 ; it < fNTimeBin ; it++ ) { - if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed - } - of << Form("\n"); - } - break; - case 'S' : - case 's' : - Int_t s = ProduceRawStream( tempbuf, 1024 ); - of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion()); - of << Form(" address data\n"); - for( int i = 0 ; i < s ; i++ ) { - of << Form(" %04x %08x\n", i, tempbuf[i]); + fZSMap[iAdc] &= ~((1-supp) << it); + if( eBIN == 0 ) { // neighbour sensitivity + fZSMap[iAdc-1] &= ~((1-supp) << it); + fZSMap[iAdc+1] &= ~((1-supp) << it); } } + } } -void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label) +void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label[]) { - // Add the given hit to the fit register which is lateron used for - // the tracklet calculation. - // In addition to the fit sums in the fit register MC information + // Add the given hit to the fit register which is lateron used for + // the tracklet calculation. + // In addition to the fit sums in the fit register MC information // is stored. - if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) && - (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0))) + if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0, fDetector, fRobPos, fMcmPos)) && + (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0, fDetector, fRobPos, fMcmPos))) fFitReg[adc].fQ0 += qtot; - - if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) && - (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))) + + if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1, fDetector, fRobPos, fMcmPos)) && + (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1, fDetector, fRobPos, fMcmPos))) fFitReg[adc].fQ1 += qtot; - - if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) && - (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE))) + + if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS, fDetector, fRobPos, fMcmPos) ) && + (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE, fDetector, fRobPos, fMcmPos))) { fFitReg[adc].fSumX += timebin; fFitReg[adc].fSumX2 += timebin*timebin; @@ -1085,6 +1194,9 @@ void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Sh fFitReg[adc].fSumY += ypos; fFitReg[adc].fSumY2 += ypos*ypos; fFitReg[adc].fSumXY += timebin*ypos; + AliDebug(10, Form("fitreg[%2i] in timebin %2i: X=%i, X2=%i, N=%i, Y=%i, Y2=%i, XY=%i, Q0=%i, Q1=%i", + adc, timebin, fFitReg[adc].fSumX, fFitReg[adc].fSumX2, fFitReg[adc].fNhits, + fFitReg[adc].fSumY, fFitReg[adc].fSumY2, fFitReg[adc].fSumXY, fFitReg[adc].fQ0, fFitReg[adc].fQ1)); } // register hits (MC info) @@ -1092,45 +1204,56 @@ void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Sh fHits[fNHits].fQtot = qtot; fHits[fNHits].fYpos = ypos; fHits[fNHits].fTimebin = timebin; - fHits[fNHits].fLabel = label; + fHits[fNHits].fLabel[0] = label[0]; + fHits[fNHits].fLabel[1] = label[1]; + fHits[fNHits].fLabel[2] = label[2]; fNHits++; } -void AliTRDmcmSim::CalcFitreg() +void AliTRDmcmSim::CalcFitreg() { // Preprocessing. // Detect the hits and fill the fit registers. - // Requires 12-bit data from fADCF which means Filter() + // Requires 12-bit data from fADCF which means Filter() // has to be called before even if all filters are bypassed. - //??? - // TRAP parameters: - const UShort_t lutPos[128] = { // move later to some other file - 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, - 16, 16, 16, 17, 17, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26, - 27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 27, 27, 27, 27, 26, - 26, 26, 26, 25, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 20, 20, 19, 18, 18, 17, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 7}; - //??? to be clarified: UInt_t adcMask = 0xffffffff; - - UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp; + + Bool_t hitQual; + Int_t adcLeft, adcCentral, adcRight; + UShort_t timebin, adcch, timebin1, timebin2, qtotTemp; Short_t ypos, fromLeft, fromRight, found; - UShort_t qTotal[19]; // the last is dummy + UShort_t qTotal[19+1]; // the last is dummy UShort_t marked[6], qMarked[6], worse1, worse2; - - timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS); - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) - < timebin1) - timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0); - timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE); - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1) - > timebin2) - timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1); + + if (fgStoreClusters) { + timebin1 = 0; + timebin2 = fNTimeBin; + } + else { + // find first timebin to be looked at + timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS, fDetector, fRobPos, fMcmPos); + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0, fDetector, fRobPos, fMcmPos) + < timebin1) + timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0, fDetector, fRobPos, fMcmPos); + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1, fDetector, fRobPos, fMcmPos) + < timebin1) + timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1, fDetector, fRobPos, fMcmPos); + + // find last timebin to be looked at + timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE, fDetector, fRobPos, fMcmPos); + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0, fDetector, fRobPos, fMcmPos) + > timebin2) + timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0, fDetector, fRobPos, fMcmPos); + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1, fDetector, fRobPos, fMcmPos) + > timebin2) + timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1, fDetector, fRobPos, fMcmPos); + } // reset the fit registers - fNHits = 0; - for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels + fNHits = 0; + for (adcch = 0; adcch < AliTRDfeeParam::GetNadcMcm()-2; adcch++) // due to border channels { fFitReg[adcch].fNhits = 0; fFitReg[adcch].fQ0 = 0; @@ -1141,35 +1264,54 @@ void AliTRDmcmSim::CalcFitreg() fFitReg[adcch].fSumY2 = 0; fFitReg[adcch].fSumXY = 0; } - + for (timebin = timebin1; timebin < timebin2; timebin++) { // first find the hit candidates and store the total cluster charge in qTotal array // in case of not hit store 0 there. - for (adcch = 0; adcch < fNADC-2; adcch++) { + for (adcch = 0; adcch < AliTRDfeeParam::GetNadcMcm()-2; adcch++) { if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS { adcLeft = fADCF[adcch ][timebin]; adcCentral = fADCF[adcch+1][timebin]; adcRight = fADCF[adcch+2][timebin]; - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1) - hitQual = ( (adcLeft * adcRight) < - (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) ); - else - hitQual = 1; + + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY, fDetector, fRobPos, fMcmPos) == 0) { + // bypass the cluster verification + hitQual = kTRUE; + } + else { + hitQual = ( (adcLeft * adcRight) < + ((fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT, fDetector, fRobPos, fMcmPos) * adcCentral*adcCentral) >> 10) ); + if (hitQual) + AliDebug(5, Form("cluster quality cut passed with %3i, %3i, %3i - threshold %3i -> %i", + adcLeft, adcCentral, adcRight, + fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT, fDetector, fRobPos, fMcmPos), + fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT, fDetector, fRobPos, fMcmPos) * adcCentral*adcCentral)); + } + // The accumulated charge is with the pedestal!!! qtotTemp = adcLeft + adcCentral + adcRight; + if ((fDebugStream) && (qtotTemp > 130)) { + (*fDebugStream) << "testtree" + << "qtot=" << qtotTemp + << "qleft=" << adcLeft + << "qcent=" << adcCentral + << "qright=" << adcRight + << "\n"; + } if ( (hitQual) && - (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) && + (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT, fDetector, fRobPos, fMcmPos)) && (adcLeft <= adcCentral) && (adcCentral > adcRight) ) qTotal[adcch] = qtotTemp; else qTotal[adcch] = 0; - //printf("ch %2d qTotal %5d\n",adcch, qTotal[adcch]); } else qTotal[adcch] = 0; //jkl + if (qTotal[adcch] != 0) + AliDebug(10,Form("ch %2d qTotal %5d",adcch, qTotal[adcch])); } fromLeft = -1; @@ -1188,7 +1330,7 @@ void AliTRDmcmSim::CalcFitreg() } adcch++; } - + fromRight = -1; adcch = 18; found = 0; @@ -1203,12 +1345,12 @@ void AliTRDmcmSim::CalcFitreg() adcch--; } - //printf("Fromleft=%d, Fromright=%d\n",fromLeft, fromRight); + AliDebug(10,Form("Fromleft=%d, Fromright=%d",fromLeft, fromRight)); // here mask the hit candidates in the middle, if any if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight)) for (adcch = fromLeft+1; adcch < fromRight; adcch++) qTotal[adcch] = 0; - + found = 0; for (adcch = 0; adcch < 19; adcch++) if (qTotal[adcch] > 0) found++; @@ -1220,10 +1362,9 @@ void AliTRDmcmSim::CalcFitreg() for (found=0; found<6; found++) { qMarked[found] = qTotal[marked[found]] >> 4; - //printf("ch_%d qTotal %d qTotals %d |",marked[found],qTotal[marked[found]],qMarked[found]); + AliDebug(10,Form("ch_%d qTotal %d qTotals %d",marked[found],qTotal[marked[found]],qMarked[found])); } - //printf("\n"); - + Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5], qMarked[0], qMarked[3], @@ -1236,15 +1377,15 @@ void AliTRDmcmSim::CalcFitreg() if (worse1 < 19) { qTotal[worse1] = 0; - //printf("Kill ch %d\n",worse1); + AliDebug(10,Form("Kill ch %d\n",worse1)); } if (worse2 < 19) { qTotal[worse2] = 0; - //printf("Kill ch %d\n",worse2); + AliDebug(10,Form("Kill ch %d\n",worse2)); } } - + for (adcch = 0; adcch < 19; adcch++) { if (qTotal[adcch] > 0) // the channel is marked for processing { @@ -1253,11 +1394,11 @@ void AliTRDmcmSim::CalcFitreg() adcRight = fADCF[adcch+2][timebin]; // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels! // subtract the pedestal TPFP, clipping instead of wrapping - - Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP); -// printf("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n", -// timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP, -// fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)); + + Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP, fDetector, fRobPos, fMcmPos); + AliDebug(10, Form("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n", + timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP, + fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT, fDetector, fRobPos, fMcmPos))); if (adcLeft < regTPFP) adcLeft = 0; else adcLeft -= regTPFP; if (adcCentral < regTPFP) adcCentral = 0; else adcCentral -= regTPFP; @@ -1267,71 +1408,79 @@ void AliTRDmcmSim::CalcFitreg() // checking for adcCentral != 0 (in case of "bad" configuration) if (adcCentral == 0) continue; - ypos = 128*(adcLeft - adcRight) / adcCentral; + ypos = 128*(adcRight - adcLeft) / adcCentral; if (ypos < 0) ypos = -ypos; - // make the correction using the LUT - ypos = ypos + lutPos[ypos & 0x7F]; + // make the correction using the position LUT + ypos = ypos + fTrapConfig->GetTrapReg((AliTRDtrapConfig::TrapReg_t) (AliTRDtrapConfig::kTPL00 + (ypos & 0x7F)), + fDetector, fRobPos, fMcmPos); if (adcLeft > adcRight) ypos = -ypos; - // label calculation - Int_t mcLabel = -1; + // label calculation (up to 3) + Int_t mcLabel[] = {-1, -1, -1}; if (fDigitsManager) { - Int_t label[9] = { 0 }; // up to 9 different labels possible - Int_t count[9] = { 0 }; - Int_t maxIdx = -1; - Int_t maxCount = 0; + const Int_t maxLabels = 9; + Int_t label[maxLabels] = { 0 }; // up to 9 different labels possible + Int_t count[maxLabels] = { 0 }; Int_t nLabels = 0; - Int_t padcol[3]; + Int_t padcol[3]; padcol[0] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch); padcol[1] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+1); padcol[2] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+2); Int_t padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos); for (Int_t iDict = 0; iDict < 3; iDict++) { - if (!fDigitsManager->UsesDictionaries() || fDigitsManager->GetDictionary(fDetector, iDict) == 0) { - AliError("Cannot get dictionary"); - continue; - } - AliTRDarrayDictionary *dict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict); - if (dict->GetDim() == 0) { - AliError("Dictionary has dim. 0"); + if (!fDict[iDict]) continue; - } - dict->Expand(); for (Int_t iPad = 0; iPad < 3; iPad++) { - if (padcol[iPad] < 0) + if (padcol[iPad] < 0) continue; - Int_t currLabel = dict->GetData(padrow, padcol[iPad], timebin); //fDigitsManager->GetTrack(iDict, padrow, padcol, timebin, fDetector); -// printf("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin); + Int_t currLabel = fDict[iDict]->GetData(padrow, padcol[iPad], timebin); + AliDebug(10, Form("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin)); for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) { if (currLabel == label[iLabel]) { count[iLabel]++; - if (count[iLabel] > maxCount) { - maxCount = count[iLabel]; - maxIdx = iLabel; - } - currLabel = 0; + currLabel = -1; break; } - } - if (currLabel > 0) { - label[nLabels++] = currLabel; + } + if (currLabel >= 0) { + label[nLabels] = currLabel; + count[nLabels] = 1; + nLabels++; } } } - if (maxIdx >= 0) - mcLabel = label[maxIdx]; + Int_t index[2*maxLabels]; + TMath::Sort(maxLabels, count, index); + for (Int_t i = 0; i < 3; i++) { + if (count[index[i]] <= 0) + break; + mcLabel[i] = label[index[i]]; + } } // add the hit to the fitregister - AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, mcLabel); + AddHitToFitreg(adcch, timebin, qTotal[adcch] >> fgkAddDigits, ypos, mcLabel); } } } + + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { + if (fFitReg[iAdc].fNhits != 0) { + AliDebug(2, Form("fitreg[%i]: nHits = %i, sumX = %i, sumY = %i, sumX2 = %i, sumY2 = %i, sumXY = %i", iAdc, + fFitReg[iAdc].fNhits, + fFitReg[iAdc].fSumX, + fFitReg[iAdc].fSumY, + fFitReg[iAdc].fSumX2, + fFitReg[iAdc].fSumY2, + fFitReg[iAdc].fSumXY + )); + } + } } -void AliTRDmcmSim::TrackletSelection() +void AliTRDmcmSim::TrackletSelection() { - // Select up to 4 tracklet candidates from the fit registers + // Select up to 4 tracklet candidates from the fit registers // and assign them to the CPUs. UShort_t adcIdx, i, j, ntracks, tmp; @@ -1339,18 +1488,19 @@ void AliTRDmcmSim::TrackletSelection() ntracks = 0; for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs - if ( (fFitReg[adcIdx].fNhits - >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) && + if ( (fFitReg[adcIdx].fNhits + >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL, fDetector, fRobPos, fMcmPos)) && (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits - >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT))) + >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT, fDetector, fRobPos, fMcmPos))) { trackletCand[ntracks][0] = adcIdx; trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits; - //printf("%d %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]); + AliDebug(10,Form("%d %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1])); ntracks++; }; - // for (i=0; i 4) { @@ -1375,7 +1525,7 @@ void AliTRDmcmSim::TrackletSelection() ntracks = 4; // cut the rest, 4 is the max } // else is not necessary to sort - + // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP for (j = 0; j < (ntracks-1); j++) { @@ -1397,67 +1547,49 @@ void AliTRDmcmSim::TrackletSelection() fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i] for (i = ntracks; i < 4; i++) // CPUs without tracklets fFitPtr[i] = 31; // pointer to the left channel with tracklet for CPU[i] = 31 (invalid) -// printf("found %i tracklet candidates\n", ntracks); -// for (i = 0; i < 4; i++) -// printf("fitPtr[%i]: %i\n", i, fFitPtr[i]); + AliDebug(10,Form("found %i tracklet candidates\n", ntracks)); + for (i = 0; i < 4; i++) + AliDebug(10,Form("fitPtr[%i]: %i\n", i, fFitPtr[i])); } void AliTRDmcmSim::FitTracklet() { - // Perform the actual tracklet fit based on the fit sums - // which have been filled in the fit registers. + // Perform the actual tracklet fit based on the fit sums + // which have been filled in the fit registers. // parameters in fitred.asm (fit program) - Int_t decPlaces = 5; Int_t rndAdd = 0; - if (decPlaces > 1) + Int_t decPlaces = 5; // must be larger than 1 or change the following code + // if (decPlaces > 1) rndAdd = (1 << (decPlaces-1)) + 1; - else if (decPlaces == 1) - rndAdd = 1; + // else if (decPlaces == 1) + // rndAdd = 1; - // should come from trapConfig (DMEM) - AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector); + Int_t ndriftDp = 5; // decimal places for drift time Long64_t shift = ((Long64_t) 1 << 32); - UInt_t scaleY = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4))); - UInt_t scaleD = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4))); - Float_t scaleSlope = (256 / pp->GetWidthIPad()) * (1 << decPlaces); -// printf("scaleSlope: %f \n", scaleSlope); - int padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos); - int yoffs = (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19) - fFeeParam->GetNcol()/2) << (8 + decPlaces); - int ndrift = 20; //??? value in simulation? - Int_t deflCorr = -1 * (Int_t) (TMath::Tan(fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector))) * fGeo->CdrHght() * scaleSlope); // -370; - Int_t tiltCorr = -1 * (Int_t) (pp->GetRowPos(padrow) / fGeo->GetTime0(fDetector % 6) * fGeo->CdrHght() * scaleSlope * - TMath::Tan(pp->GetTiltingAngle() / 180. * TMath::Pi())); -// printf("vdrift av.: %f\n", fCal->GetVdriftAverage(fDetector)); -// printf("chamber height: %f\n", fGeo->CdrHght()); -// printf("omega tau: %f\n", fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector))); -// printf("deflection correction: %i\n", deflCorr); - Float_t ptcut = 2.3; - AliMagF* fld = (AliMagF *) TGeoGlobalMagField::Instance()->GetField(); - Double_t bz = 0; - if (fld) { - bz = 0.1 * fld->SolenoidField(); // kGauss -> Tesla - } -// printf("Bz: %f\n", bz); - Float_t x0 = fGeo->GetTime0(fDetector % 6); - Float_t y0 = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10)); - Float_t alphaMax = TMath::ASin( (TMath::Sqrt(TMath::Power(x0/100., 2) + TMath::Power(y0/100., 2)) * - 0.3 * TMath::Abs(bz) ) / (2 * ptcut)); -// printf("alpha max: %f\n", alphaMax * 180/TMath::Pi()); - Int_t minslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) + alphaMax) * scaleSlope); - Int_t maxslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) - alphaMax) * scaleSlope); -// printf("min y-defl: %i\n", minslope); -// printf("max y-defl: %i\n", maxslope); + + // calculated in fitred.asm + Int_t padrow = ((fRobPos >> 1) << 2) | (fMcmPos >> 2); + Int_t yoffs = (((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - + ((18*4*2 - 18*2 - 1) << 7); + yoffs = yoffs << decPlaces; // holds position of ADC channel 1 + Int_t layer = fDetector % 6; + UInt_t scaleY = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 160.0e-4) * shift); + UInt_t scaleD = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 140.0e-4) * shift); + + Int_t deflCorr = (Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrDeflCorr, fDetector, fRobPos, fMcmPos); + Int_t ndrift = (Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrNdrift, fDetector, fRobPos, fMcmPos); // local variables for calculation Long64_t mult, temp, denom; //??? - UInt_t q0, q1, qTotal; // charges in the two windows and total charge + UInt_t q0, q1, pid; // charges in the two windows and total charge UShort_t nHits; // number of hits Int_t slope, offset; // slope and offset of the tracklet Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers - //int32_t SumY2; // not used in the current TRAP program + Int_t sumY2; // not used in the current TRAP program, now used for error calculation (simulation only) + Float_t fitError, fitSlope, fitOffset; FitReg_t *fit0, *fit1; // pointers to relevant fit registers - + // const uint32_t OneDivN[32] = { // 2**31/N : exactly like in the TRAP, the simple division here gives the same result! // 0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490, // 0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880, @@ -1467,7 +1599,7 @@ void AliTRDmcmSim::FitTracklet() for (Int_t cpu = 0; cpu < 4; cpu++) { if (fFitPtr[cpu] == 31) { - fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker(); + fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker(); } else { @@ -1482,113 +1614,189 @@ void AliTRDmcmSim::FitTracklet() nHits = fit0->fNhits + fit1->fNhits; // number of hits sumX = fit0->fSumX + fit1->fSumX; sumX2 = fit0->fSumX2 + fit1->fSumX2; - denom = nHits*sumX2 - sumX*sumX; + denom = ((Long64_t) nHits)*((Long64_t) sumX2) - ((Long64_t) sumX)*((Long64_t) sumX); mult = mult / denom; // exactly like in the TRAP program q0 = fit0->fQ0 + fit1->fQ0; q1 = fit0->fQ1 + fit1->fQ1; sumY = fit0->fSumY + fit1->fSumY + 256*fit1->fNhits; sumXY = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX; + sumY2 = fit0->fSumY2 + fit1->fSumY2 + 512*fit1->fSumY + 256*256*fit1->fNhits; slope = nHits*sumXY - sumX * sumY; -// printf("slope from fitreg: %i\n", slope); offset = sumX2*sumY - sumX * sumXY; temp = mult * slope; slope = temp >> 32; // take the upper 32 bits + slope = -slope; temp = mult * offset; offset = temp >> 32; // take the upper 32 bits - offset = offset + yoffs + (18 << (8 + decPlaces)); -// printf("slope: %i, slope * ndrift: %i, deflCorr: %i, tiltCorr: %i\n", slope, slope * ndrift, deflCorr, tiltCorr); - slope = slope * ndrift + deflCorr + tiltCorr; + offset = offset + yoffs; + AliDebug(10, Form("slope = %i, slope * ndrift = %i, deflCorr: %i", + slope, slope * ndrift, deflCorr)); + slope = ((slope * ndrift) >> ndriftDp) + deflCorr; offset = offset - (fFitPtr[cpu] << (8 + decPlaces)); - -// printf("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i ", fDetector, fRobPos, fMcmPos, slope, minslope, maxslope); + + temp = slope; + temp = temp * scaleD; + slope = (temp >> 32); + temp = offset; + temp = temp * scaleY; + offset = (temp >> 32); + + // rounding, like in the TRAP + slope = (slope + rndAdd) >> decPlaces; + offset = (offset + rndAdd) >> decPlaces; + + AliDebug(5, Form("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i", + fDetector, fRobPos, fMcmPos, slope, + (Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrDeflCutStart + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos), + (Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrDeflCutStart + 1 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos))); + + AliDebug(5, Form("Fit sums: x = %i, X = %i, y = %i, Y = %i, Z = %i, q0 = %i, q1 = %i", + sumX, sumX2, sumY, sumY2, sumXY, q0, q1)); + + fitSlope = (Float_t) (nHits * sumXY - sumX * sumY) / (nHits * sumX2 - sumX*sumX); + + fitOffset = (Float_t) (sumX2 * sumY - sumX * sumXY) / (nHits * sumX2 - sumX*sumX); + + Float_t sx = (Float_t) sumX; + Float_t sx2 = (Float_t) sumX2; + Float_t sy = (Float_t) sumY; + Float_t sy2 = (Float_t) sumY2; + Float_t sxy = (Float_t) sumXY; + fitError = sy2 - (sx2 * sy*sy - 2 * sx * sxy * sy + nHits * sxy*sxy) / (nHits * sx2 - sx*sx); + //fitError = (Float_t) sumY2 - (Float_t) (sumY*sumY) / nHits - fitSlope * ((Float_t) (sumXY - sumX*sumY) / nHits); + Bool_t rejected = kFALSE; - if (GetApplyCut() && ((slope < minslope) || (slope > maxslope))) + // deflection range table from DMEM + if ((slope < ((Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrDeflCutStart + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos))) || + (slope > ((Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrDeflCutStart + 1 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)))) rejected = kTRUE; - if (rejected) + + if (rejected && GetApplyCut()) { -// printf("rejected\n"); fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker(); } else { -// printf("accepted\n"); - temp = slope; - temp = temp * scaleD; - slope = (temp >> 32); -// printf("slope after scaling: %i\n", slope); - - temp = offset; - temp = temp * scaleY; - offset = (temp >> 32); - - // rounding, like in the TRAP - slope = (slope + rndAdd) >> decPlaces; -// printf("slope after shifting: %i\n", slope); - offset = (offset + rndAdd) >> decPlaces; - - if (slope > 63) { // wrapping in TRAP! - AliError(Form("Overflow in slope: %i, tracklet discarded!", slope)); - fMCMT[cpu] = 0x10001000; - continue; - } - else if (slope < -64) { - AliError(Form("Underflow in slope: %i, tracklet discarded!", slope)); + if (slope > 63 || slope < -64) { // wrapping in TRAP! + AliDebug(1,Form("Overflow in slope: %i, tracklet discarded!", slope)); fMCMT[cpu] = 0x10001000; continue; } - else { - slope = slope & 0x7F; // 7 bit - } -// printf("slope after clipping: 0x%02x\n", slope); - if (offset > 0xfff || offset < -0xfff) + slope = slope & 0x7F; // 7 bit + + if (offset > 0xfff || offset < -0xfff) AliWarning("Overflow in offset"); offset = offset & 0x1FFF; // 13 bit - qTotal = (q1 / nHits) >> 1; - if (qTotal > 0xff) - AliWarning("Overflow in charge"); - qTotal = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program + pid = GetPID(q0, q1); - // assemble and store the tracklet word - fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset; + if (pid > 0xff) + AliWarning("Overflow in PID"); + pid = pid & 0xFF; // 8 bit, exactly like in the TRAP program - // calculate MC label - Int_t mcLabel = -1; - if (fDigitsManager) { - Int_t label[30] = {0}; // up to 30 different labels possible - Int_t count[30] = {0}; - Int_t maxIdx = -1; - Int_t maxCount = 0; - Int_t nLabels = 0; - for (Int_t iHit = 0; iHit < fNHits; iHit++) { - if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) || - (fHits[iHit].fChannel - fFitPtr[cpu] > 1)) - continue; - Int_t currLabel = fHits[iHit].fLabel; - for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) { - if (currLabel == label[iLabel]) { - count[iLabel]++; - if (count[iLabel] > maxCount) { - maxCount = count[iLabel]; - maxIdx = iLabel; - } - currLabel = 0; - break; - } - } - if (currLabel > 0) { - label[nLabels++] = currLabel; - } - } - if (maxIdx >= 0) - mcLabel = label[maxIdx]; + // assemble and store the tracklet word + fMCMT[cpu] = (pid << 24) | (padrow << 20) | (slope << 13) | offset; + + // calculate number of hits and MC label + Int_t mcLabel[] = { -1, -1, -1}; + Int_t nHits0 = 0; + Int_t nHits1 = 0; + + const Int_t maxLabels = 30; + Int_t label[maxLabels] = {0}; // up to 30 different labels possible + Int_t count[maxLabels] = {0}; + Int_t nLabels = 0; + + for (Int_t iHit = 0; iHit < fNHits; iHit++) { + if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) || + (fHits[iHit].fChannel - fFitPtr[cpu] > 1)) + continue; + + // counting contributing hits + if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0, fDetector, fRobPos, fMcmPos) && + fHits[iHit].fTimebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0, fDetector, fRobPos, fMcmPos)) + nHits0++; + if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1, fDetector, fRobPos, fMcmPos) && + fHits[iHit].fTimebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1, fDetector, fRobPos, fMcmPos)) + nHits1++; + + // label calculation only if there is a digitsmanager to get the labels from + if (fDigitsManager) { + for (Int_t i = 0; i < 3; i++) { + Int_t currLabel = fHits[iHit].fLabel[i]; + for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) { + if (currLabel == label[iLabel]) { + count[iLabel]++; + currLabel = -1; + break; + } + } + if (currLabel >= 0 && nLabels < maxLabels) { + label[nLabels] = currLabel; + count[nLabels]++; + nLabels++; + } + } + } + + if (fDigitsManager) { + Int_t index[2*maxLabels]; + TMath::Sort(maxLabels, count, index); + for (Int_t i = 0; i < 3; i++) { + if (count[index[i]] <= 0) + break; + mcLabel[i] = label[index[i]]; + } + } } new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos); ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel); + + + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits(fit0->fNhits + fit1->fNhits); + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits0(nHits0); + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits1(nHits1); + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0); + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1); + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetSlope(fitSlope); + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetOffset(fitOffset); + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetError(TMath::Sqrt(TMath::Abs(fitError)/nHits)); + + // store cluster information (if requested) + if (fgStoreClusters) { + Float_t *res = new Float_t[fNTimeBin]; + Float_t *qtot = new Float_t[fNTimeBin]; + for (Int_t iTimebin = 0; iTimebin < fNTimeBin; ++iTimebin) { + res[iTimebin] = 0; + qtot[iTimebin] = 0; + } + for (Int_t iHit = 0; iHit < fNHits; iHit++) { + Int_t timebin = fHits[iHit].fTimebin; + + // check if hit contributes + if (fHits[iHit].fChannel == fFitPtr[cpu]) { + res[timebin] = fHits[iHit].fYpos - (fitSlope * timebin + fitOffset); + qtot[timebin] = fHits[iHit].fQtot; + } + else if (fHits[iHit].fChannel == fFitPtr[cpu] + 1) { + res[timebin] = fHits[iHit].fYpos + 256 - (fitSlope * timebin + fitOffset); + qtot[timebin] = fHits[iHit].fQtot; + } + } + ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetClusters(res, qtot, fNTimeBin); + delete [] res; + delete [] qtot; + } + + if (fitError < 0) + AliError(Form("Strange fit error: %f from Sx: %i, Sy: %i, Sxy: %i, Sx2: %i, Sy2: %i, nHits: %i", + fitError, sumX, sumY, sumXY, sumX2, sumY2, nHits)); + AliDebug(3, Form("fit slope: %f, offset: %f, error: %f", + fitSlope, fitOffset, TMath::Sqrt(TMath::Abs(fitError)/nHits))); } } } @@ -1598,7 +1806,7 @@ void AliTRDmcmSim::Tracklet() { // Run the tracklet calculation by calling sequentially: // CalcFitreg(); TrackletSelection(); FitTracklet() - // and store the tracklets + // and store the tracklets if (!fInitialized) { AliError("Called uninitialized! Nothing done!"); @@ -1614,9 +1822,11 @@ void AliTRDmcmSim::Tracklet() FitTracklet(); } -Bool_t AliTRDmcmSim::StoreTracklets() +Bool_t AliTRDmcmSim::StoreTracklets() { - if (fTrackletArray->GetEntriesFast() == 0) + // store the found tracklets via the loader + + if (fTrackletArray->GetEntriesFast() == 0) return kTRUE; AliRunLoader *rl = AliRunLoader::Instance(); @@ -1633,19 +1843,17 @@ Bool_t AliTRDmcmSim::StoreTracklets() dl->MakeTree(); trackletTree = dl->Tree(); } - + AliTRDtrackletMCM *trkl = 0x0; - TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch"); + TBranch *trkbranch = trackletTree->GetBranch(fTrklBranchName.Data()); if (!trkbranch) - trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000); - + trkbranch = trackletTree->Branch(fTrklBranchName.Data(), "AliTRDtrackletMCM", &trkl, 32000); + for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) { trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]); trkbranch->SetAddress(&trkl); -// printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord()); trkbranch->Fill(); } - dl->WriteData("OVERWRITE"); return kTRUE; } @@ -1656,68 +1864,136 @@ void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits) // EBSF = 1: unfiltered data; EBSF = 0: filtered data // zero-suppressed valued are written as -1 to digits - if (!fInitialized) { - AliError("Called uninitialized! Nothing done!"); + if( !CheckInitialized() ) return; - } - Int_t firstAdc = 0; - Int_t lastAdc = fNADC - 1; + Int_t offset = (fMcmPos % 4 + 1) * 21 + (fRobPos % 2) * 84 - 1; - while (GetCol(firstAdc) < 0) - firstAdc++; - - while (GetCol(lastAdc) < 0) - lastAdc--; - - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF, fDetector, fRobPos, fMcmPos) != 0) // store unfiltered data { - for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) { - if (fZSM1Dim[iAdc] == 1) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { + if (~fZSMap[iAdc] == 0) { for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1); -// printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin, -// digits->GetData(GetRow(), GetCol(iAdc), iTimeBin)); + digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, -1); + } + } + else if (iAdc < 2 || iAdc == 20) { + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { + digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, (fADCR[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline); } } } } else { - for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) { - if (fZSM1Dim[iAdc] == 0) { + for (Int_t iAdc = 0; iAdc < AliTRDfeeParam::GetNadcMcm(); iAdc++) { + if (~fZSMap[iAdc] != 0) { for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, fADCF[iAdc][iTimeBin] >> fgkAddDigits); + digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, (fADCF[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline); } } else { for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1); -// printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin); + digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, -1); } } } } } + +// ****************************** +// PID section +// +// Memory area for the LUT: 0xC100 to 0xC3FF +// +// The addresses for the parameters (the order is optimized for maximum calculation speed in the MCMs): +// 0xC028: cor1 +// 0xC029: nBins(sF) +// 0xC02A: cor0 +// 0xC02B: TableLength +// Defined in AliTRDtrapConfig.h +// +// The algorithm implemented in the TRAP program of the MCMs (Venelin Angelov) +// 1) set the read pointer to the beginning of the Parameters in DMEM +// 2) shift right the FitReg with the Q0 + (Q1 << 16) to get Q1 +// 3) read cor1 with rpointer++ +// 4) start cor1*Q1 +// 5) read nBins with rpointer++ +// 6) start nBins*cor1*Q1 +// 7) read cor0 with rpointer++ +// 8) swap hi-low parts in FitReg, now is Q1 + (Q0 << 16) +// 9) shift right to get Q0 +// 10) start cor0*Q0 +// 11) read TableLength +// 12) compare cor0*Q0 with nBins +// 13) if >=, clip cor0*Q0 to nBins-1 +// 14) add cor0*Q0 to nBins*cor1*Q1 +// 15) compare the result with TableLength +// 16) if >=, clip to TableLength-1 +// 17) read from the LUT 8 bits + + +Int_t AliTRDmcmSim::GetPID(Int_t q0, Int_t q1) +{ + // return PID calculated from charges accumulated in two time windows + + ULong64_t addrQ0; + ULong64_t addr; + + UInt_t nBinsQ0 = fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTnbins, fDetector, fRobPos, fMcmPos); // number of bins in q0 / 4 !! + UInt_t pidTotalSize = fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTLength, fDetector, fRobPos, fMcmPos); + if(nBinsQ0==0 || pidTotalSize==0) // make sure we don't run into trouble if the value for Q0 is not configured + return 0; // Q1 not configured is ok for 1D LUT + + ULong_t corrQ0 = fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTcor0, fDetector, fRobPos, fMcmPos); + ULong_t corrQ1 = fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTcor1, fDetector, fRobPos, fMcmPos); + if(corrQ0==0) // make sure we don't run into trouble if one of the values is not configured + return 0; + + addrQ0 = corrQ0; + addrQ0 = (((addrQ0*q0)>>16)>>16); // because addrQ0 = (q0 * corrQ0) >> 32; does not work for unknown reasons + + if(addrQ0 >= nBinsQ0) { // check for overflow + AliDebug(5,Form("Overflow in q0: %llu/4 is bigger then %u", addrQ0, nBinsQ0)); + addrQ0 = nBinsQ0 -1; + } + + addr = corrQ1; + addr = (((addr*q1)>>16)>>16); + addr = addrQ0 + nBinsQ0*addr; // because addr = addrQ0 + nBinsQ0* (((corrQ1*q1)>>32); does not work + + if(addr >= pidTotalSize) { + AliDebug(5,Form("Overflow in q1. Address %llu/4 is bigger then %u", addr, pidTotalSize)); + addr = pidTotalSize -1; + } + + // For a LUT with 11 input and 8 output bits, the first memory address is set to LUT[0] | (LUT[1] << 8) | (LUT[2] << 16) | (LUT[3] << 24) + // and so on + UInt_t result = fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTStart+(addr/4), fDetector, fRobPos, fMcmPos); + return (result>>((addr%4)*8)) & 0xFF; +} + + + // help functions, to be cleaned up UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const { - // - // This function adds a and b (unsigned) and clips to - // the specified number of bits. - // + // + // This function adds a and b (unsigned) and clips to + // the specified number of bits. + // UInt_t sum = a + b; if (nbits < 32) { UInt_t maxv = (1 << nbits) - 1;; - if (sum > maxv) + if (sum > maxv) sum = maxv; } else { - if ((sum < a) || (sum < b)) + if ((sum < a) || (sum < b)) sum = 0xFFFFFFFF; } return sum; @@ -1725,8 +2001,8 @@ UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const void AliTRDmcmSim::Sort2(UShort_t idx1i, UShort_t idx2i, \ UShort_t val1i, UShort_t val2i, \ - UShort_t *idx1o, UShort_t *idx2o, \ - UShort_t *val1o, UShort_t *val2o) const + UShort_t * const idx1o, UShort_t * const idx2o, \ + UShort_t * const val1o, UShort_t * const val2o) const { // sorting for tracklet selection @@ -1748,18 +2024,17 @@ void AliTRDmcmSim::Sort2(UShort_t idx1i, UShort_t idx2i, \ void AliTRDmcmSim::Sort3(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, \ UShort_t val1i, UShort_t val2i, UShort_t val3i, \ - UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, \ - UShort_t *val1o, UShort_t *val2o, UShort_t *val3o) + UShort_t * const idx1o, UShort_t * const idx2o, UShort_t * const idx3o, \ + UShort_t * const val1o, UShort_t * const val2o, UShort_t * const val3o) { // sorting for tracklet selection - int sel; + Int_t sel; if (val1i > val2i) sel=4; else sel=0; if (val2i > val3i) sel=sel + 2; if (val3i > val1i) sel=sel + 1; - //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel); switch(sel) { case 6 : // 1 > 2 > 3 => 1 2 3 @@ -1821,13 +2096,12 @@ void AliTRDmcmSim::Sort3(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, \ AliError("ERROR in Sort3!!!\n"); break; } -// printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o); } void AliTRDmcmSim::Sort6To4(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, UShort_t idx4i, UShort_t idx5i, UShort_t idx6i, \ UShort_t val1i, UShort_t val2i, UShort_t val3i, UShort_t val4i, UShort_t val5i, UShort_t val6i, \ - UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, UShort_t *idx4o, \ - UShort_t *val1o, UShort_t *val2o, UShort_t *val3o, UShort_t *val4o) + UShort_t * const idx1o, UShort_t * const idx2o, UShort_t * const idx3o, UShort_t * const idx4o, \ + UShort_t * const val1o, UShort_t * const val2o, UShort_t * const val3o, UShort_t * const val4o) { // sorting for tracklet selection @@ -1854,7 +2128,7 @@ void AliTRDmcmSim::Sort6To4(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, U void AliTRDmcmSim::Sort6To2Worst(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, UShort_t idx4i, UShort_t idx5i, UShort_t idx6i, \ UShort_t val1i, UShort_t val2i, UShort_t val3i, UShort_t val4i, UShort_t val5i, UShort_t val6i, \ - UShort_t *idx5o, UShort_t *idx6o) + UShort_t * const idx5o, UShort_t * const idx6o) { // sorting for tracklet selection @@ -1876,8 +2150,485 @@ void AliTRDmcmSim::Sort6To2Worst(UShort_t idx1i, UShort_t idx2i, UShort_t idx Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s, &dummy1, &dummy2, idx6o, &dummy3, &dummy4, &dummy5); -// printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n", -// idx21s, idx23as, idx22s, idx23bs, *idx5o, *idx6o); } +// ----- I/O implementation ----- + +ostream& AliTRDmcmSim::Text(ostream& os) +{ + // manipulator to activate output in text format (default) + + os.iword(fgkFormatIndex) = 0; + return os; +} + +ostream& AliTRDmcmSim::Cfdat(ostream& os) +{ + // manipulator to activate output in CFDAT format + // to send to the FEE via SCSN + + os.iword(fgkFormatIndex) = 1; + return os; +} + +ostream& AliTRDmcmSim::Raw(ostream& os) +{ + // manipulator to activate output as raw data dump + + os.iword(fgkFormatIndex) = 2; + return os; +} + +ostream& operator<<(ostream& os, const AliTRDmcmSim& mcm) +{ + // output implementation + + // no output for non-initialized MCM + if (!mcm.CheckInitialized()) + return os; + + // ----- human-readable output ----- + if (os.iword(AliTRDmcmSim::fgkFormatIndex) == 0) { + + os << "MCM " << mcm.fMcmPos << " on ROB " << mcm.fRobPos << + " in detector " << mcm.fDetector << std::endl; + + os << "----- Unfiltered ADC data (10 bit) -----" << std::endl; + os << "ch "; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) + os << std::setw(5) << iChannel; + os << std::endl; + for (Int_t iTimeBin = 0; iTimeBin < mcm.fNTimeBin; iTimeBin++) { + os << "tb " << std::setw(2) << iTimeBin << ":"; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) { + os << std::setw(5) << (mcm.fADCR[iChannel][iTimeBin] >> mcm.fgkAddDigits); + } + os << std::endl; + } + + os << "----- Filtered ADC data (10+2 bit) -----" << std::endl; + os << "ch "; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) + os << std::setw(4) << iChannel + << ((~mcm.fZSMap[iChannel] != 0) ? "!" : " "); + os << std::endl; + for (Int_t iTimeBin = 0; iTimeBin < mcm.fNTimeBin; iTimeBin++) { + os << "tb " << std::setw(2) << iTimeBin << ":"; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) { + os << std::setw(4) << (mcm.fADCF[iChannel][iTimeBin]) + << (((mcm.fZSMap[iChannel] & (1 << iTimeBin)) == 0) ? "!" : " "); + } + os << std::endl; + } + } + + // ----- CFDAT output ----- + else if(os.iword(AliTRDmcmSim::fgkFormatIndex) == 1) { + Int_t dest = 127; + Int_t addrOffset = 0x2000; + Int_t addrStep = 0x80; + + for (Int_t iTimeBin = 0; iTimeBin < mcm.fNTimeBin; iTimeBin++) { + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) { + os << std::setw(5) << 10 + << std::setw(5) << addrOffset + iChannel * addrStep + iTimeBin + << std::setw(5) << (mcm.fADCF[iChannel][iTimeBin]) + << std::setw(5) << dest << std::endl; + } + os << std::endl; + } + } + + // ----- raw data ouptut ----- + else if (os.iword(AliTRDmcmSim::fgkFormatIndex) == 2) { + Int_t bufSize = 300; + UInt_t *buf = new UInt_t[bufSize]; + + Int_t bufLength = mcm.ProduceRawStream(&buf[0], bufSize); + + for (Int_t i = 0; i < bufLength; i++) + std::cout << "0x" << std::hex << buf[i] << std::dec << std::endl; + + delete [] buf; + } + + else { + os << "unknown format set" << std::endl; + } + + return os; +} + + +void AliTRDmcmSim::PrintFitRegXml(ostream& os) const +{ + // print fit registres in XML format + + bool tracklet=false; + + for (Int_t cpu = 0; cpu < 4; cpu++) { + if(fFitPtr[cpu] != 31) + tracklet=true; + } + + if(tracklet==true) { + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << " " << std::endl; + os << " " << std::endl; + + for(int cpu=0; cpu<4; cpu++) { + os << " " << std::endl; + if(fFitPtr[cpu] != 31) { + for(int adcch=fFitPtr[cpu]; adcch"<< std::endl; + os << " " << fFitReg[adcch].fNhits << ""<< std::endl; + os << " " << fFitReg[adcch].fQ0 << ""<< std::endl; + os << " " << fFitReg[adcch].fQ1 << ""<< std::endl; + os << " " << fFitReg[adcch].fSumX << ""<< std::endl; + os << " " << fFitReg[adcch].fSumX2 << ""<< std::endl; + os << " " << fFitReg[adcch].fSumY << ""<< std::endl; + os << " " << fFitReg[adcch].fSumY2 << ""<< std::endl; + os << " " << fFitReg[adcch].fSumXY << ""<< std::endl; + os << " " << std::endl; + } + } + os << " " << std::endl; + } + os << " " << std::endl; + os << " " << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + } +} + + +void AliTRDmcmSim::PrintTrackletsXml(ostream& os) const +{ + // print tracklets in XML format + + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << " " << std::endl; + os << " " << std::endl; + + Int_t pid, padrow, slope, offset; + for(Int_t cpu=0; cpu<4; cpu++) { + if(fMCMT[cpu] == 0x10001000) { + pid=-1; + padrow=-1; + slope=-1; + offset=-1; + } + else { + pid = (fMCMT[cpu] & 0xFF000000) >> 24; + padrow = (fMCMT[cpu] & 0xF00000 ) >> 20; + slope = (fMCMT[cpu] & 0xFE000 ) >> 13; + offset = (fMCMT[cpu] & 0x1FFF ) ; + + } + os << " " << pid << "" << " " << padrow << "" + << " " << slope << "" << " " << offset << "" << "" << std::endl; + } + + os << " " << std::endl; + os << " " << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; +} + + +void AliTRDmcmSim::PrintAdcDatTxt(ostream& os) const +{ + // print ADC data in text format (suitable as Modelsim stimuli) + + os << "# MCM " << fMcmPos << " on ROB " << fRobPos << + " in detector " << fDetector << std::endl; + + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); ++iChannel) { + os << std::setw(5) << (fADCR[iChannel][iTimeBin] >> fgkAddDigits); + } + os << std::endl; + } +} + + +void AliTRDmcmSim::PrintAdcDatHuman(ostream& os) const +{ + // print ADC data in human-readable format + + os << "MCM " << fMcmPos << " on ROB " << fRobPos << + " in detector " << fDetector << std::endl; + + os << "----- Unfiltered ADC data (10 bit) -----" << std::endl; + os << "ch "; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) + os << std::setw(5) << iChannel; + os << std::endl; + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { + os << "tb " << std::setw(2) << iTimeBin << ":"; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) { + os << std::setw(5) << (fADCR[iChannel][iTimeBin] >> fgkAddDigits); + } + os << std::endl; + } + + os << "----- Filtered ADC data (10+2 bit) -----" << std::endl; + os << "ch "; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) + os << std::setw(4) << iChannel + << ((~fZSMap[iChannel] != 0) ? "!" : " "); + os << std::endl; + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { + os << "tb " << std::setw(2) << iTimeBin << ":"; + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) { + os << std::setw(4) << (fADCF[iChannel][iTimeBin]) + << (((fZSMap[iChannel] & (1 << iTimeBin)) == 0) ? "!" : " "); + } + os << std::endl; + } +} + + +void AliTRDmcmSim::PrintAdcDatXml(ostream& os) const +{ + // print ADC data in XML format + + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << " " << std::endl; + os << " " << std::endl; + + for(Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) { + os << " " << std::endl; + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { + os << "" << fADCF[iChannel][iTimeBin]/4 << ""; + } + os << " " << std::endl; + } + + os << " " << std::endl; + os << " " << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; + os << "" << std::endl; +} + + + +void AliTRDmcmSim::PrintAdcDatDatx(ostream& os, Bool_t broadcast, Int_t timeBinOffset) const +{ + // print ADC data in datx format (to send to FEE) + + fTrapConfig->PrintDatx(os, 2602, 1, 0, 127); // command to enable the ADC clock - necessary to write ADC values to MCM + os << std::endl; + + Int_t addrOffset = 0x2000; + Int_t addrStep = 0x80; + Int_t addrOffsetEBSIA = 0x20; + + for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { + for (Int_t iChannel = 0; iChannel < AliTRDfeeParam::GetNadcMcm(); iChannel++) { + if ((iTimeBin < timeBinOffset) || (iTimeBin >= fNTimeBin+timeBinOffset)) { + if(broadcast==kFALSE) + fTrapConfig->PrintDatx(os, addrOffset+iChannel*addrStep+addrOffsetEBSIA+iTimeBin, 10, GetRobPos(), GetMcmPos()); + else + fTrapConfig->PrintDatx(os, addrOffset+iChannel*addrStep+addrOffsetEBSIA+iTimeBin, 10, 0, 127); + } + else { + if(broadcast==kFALSE) + fTrapConfig->PrintDatx(os, addrOffset+iChannel*addrStep+addrOffsetEBSIA+iTimeBin, (fADCF[iChannel][iTimeBin-timeBinOffset]/4), GetRobPos(), GetMcmPos()); + else + fTrapConfig->PrintDatx(os, addrOffset+iChannel*addrStep+addrOffsetEBSIA+iTimeBin, (fADCF[iChannel][iTimeBin-timeBinOffset]/4), 0, 127); + } + } + os << std::endl; + } +} + + +void AliTRDmcmSim::PrintPidLutHuman() +{ + // print PID LUT in human readable format + + UInt_t result; + + UInt_t addrEnd = fgkDmemAddrLUTStart + fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTLength, fDetector, fRobPos, fMcmPos)/4; // /4 because each addr contains 4 values + UInt_t nBinsQ0 = fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTnbins, fDetector, fRobPos, fMcmPos); + + std::cout << "nBinsQ0: " << nBinsQ0 << std::endl; + std::cout << "LUT table length: " << fTrapConfig->GetDmemUnsigned(fgkDmemAddrLUTLength, fDetector, fRobPos, fMcmPos) << std::endl; + + if (nBinsQ0>0) { + for(UInt_t addr=fgkDmemAddrLUTStart; addr< addrEnd; addr++) { + result = fTrapConfig->GetDmemUnsigned(addr, fDetector, fRobPos, fMcmPos); + std::cout << addr << " # x: " << ((addr-fgkDmemAddrLUTStart)%((nBinsQ0)/4))*4 << ", y: " <<(addr-fgkDmemAddrLUTStart)/(nBinsQ0/4) + << " # " <<((result>>0)&0xFF) + << " | " << ((result>>8)&0xFF) + << " | " << ((result>>16)&0xFF) + << " | " << ((result>>24)&0xFF) << std::endl; + } + } +} + + +Bool_t AliTRDmcmSim::ReadPackedConfig(AliTRDtrapConfig *cfg, Int_t hc, UInt_t *data, Int_t size) +{ + // Read the packed configuration from the passed memory block + // + // To be used to retrieve the TRAP configuration from the + // configuration as sent in the raw data. + + AliDebugClass(1, "Reading packed configuration"); + + Int_t det = hc/2; + + Int_t idx = 0; + Int_t err = 0; + Int_t step, bwidth, nwords, exitFlag, bitcnt; + + UShort_t caddr; + UInt_t dat, msk, header, dataHi; + + while (idx < size && *data != 0x00000000) { + + Int_t rob = (*data >> 28) & 0x7; + Int_t mcm = (*data >> 24) & 0xf; + + AliDebugClass(1, Form("Config of det. %3i MCM %i:%02i (0x%08x)", det, rob, mcm, *data)); + data++; + + while (idx < size && *data != 0x00000000) { + + header = *data; + data++; + idx++; + + AliDebugClass(5, Form("read: 0x%08x", header)); + + if (header & 0x01) // single data + { + dat = (header >> 2) & 0xFFFF; // 16 bit data + caddr = (header >> 18) & 0x3FFF; // 14 bit address + + if (caddr != 0x1FFF) // temp!!! because the end marker was wrong + { + if (header & 0x02) // check if > 16 bits + { + dataHi = *data; + AliDebugClass(5, Form("read: 0x%08x", dataHi)); + data++; + idx++; + err += ((dataHi ^ (dat | 1)) & 0xFFFF) != 0; + dat = (dataHi & 0xFFFF0000) | dat; + } + AliDebugClass(5, Form("addr=0x%04x (%s) data=0x%08x\n", caddr, cfg->GetRegName(cfg->GetRegByAddress(caddr)), dat)); + if ( ! cfg->Poke(caddr, dat, det, rob, mcm) ) + AliDebugClass(5, Form("(single-write): non-existing address 0x%04x containing 0x%08x\n", caddr, header)); + if (idx > size) + { + AliDebugClass(5, Form("(single-write): no more data, missing end marker\n")); + return -err; + } + } + else + { + AliDebugClass(5, Form("(single-write): address 0x%04x => old endmarker?\n", caddr)); + return err; + } + } + + else // block of data + { + step = (header >> 1) & 0x0003; + bwidth = ((header >> 3) & 0x001F) + 1; + nwords = (header >> 8) & 0x00FF; + caddr = (header >> 16) & 0xFFFF; + exitFlag = (step == 0) || (step == 3) || (nwords == 0); + + if (exitFlag) + break; + + switch (bwidth) + { + case 15: + case 10: + case 7: + case 6: + case 5: + { + msk = (1 << bwidth) - 1; + bitcnt = 0; + while (nwords > 0) + { + nwords--; + bitcnt -= bwidth; + if (bitcnt < 0) + { + header = *data; + AliDebugClass(5, Form("read 0x%08x", header)); + data++; + idx++; + err += (header & 1); + header = header >> 1; + bitcnt = 31 - bwidth; + } + AliDebugClass(5, Form("addr=0x%04x (%s) data=0x%08x\n", caddr, cfg->GetRegName(cfg->GetRegByAddress(caddr)), header & msk)); + if ( ! cfg->Poke(caddr, header & msk, det, rob, mcm) ) + AliDebugClass(5, Form("(single-write): non-existing address 0x%04x containing 0x%08x\n", caddr, header)); + + caddr += step; + header = header >> bwidth; + if (idx >= size) + { + AliDebugClass(5, Form("(block-write): no end marker! %d words read\n", idx)); + return -err; + } + } + break; + } // end case 5-15 + case 31: + { + while (nwords > 0) + { + header = *data; + AliDebugClass(5, Form("read 0x%08x", header)); + data++; + idx++; + nwords--; + err += (header & 1); + + AliDebugClass(5, Form("addr=0x%04x (%s) data=0x%08x", caddr, cfg->GetRegName(cfg->GetRegByAddress(caddr)), header >> 1)); + if ( ! cfg->Poke(caddr, header >> 1, det, rob, mcm) ) + AliDebugClass(5, Form("(single-write): non-existing address 0x%04x containing 0x%08x\n", caddr, header)); + + caddr += step; + if (idx >= size) + { + AliDebugClass(5, Form("no end marker! %d words read", idx)); + return -err; + } + } + break; + } + default: return err; + } // end switch + } // end block case + } + } // end while + AliDebugClass(5, Form("no end marker! %d words read", idx)); + return -err; // only if the max length of the block reached! +}