X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TRD%2FAliTRDmcmSim.cxx;h=7aefe9ab531451a1e85267795d093199343a1876;hp=4f6f0907d8f3dc79f53a0c3dc0ee3f153840017c;hb=64eb00cbccc0b83ca964e23ddd2ecc5af20c6260;hpb=b0a41e803818ef232a76f440857103fc49d25f42 diff --git a/TRD/AliTRDmcmSim.cxx b/TRD/AliTRDmcmSim.cxx index 4f6f0907d8f..7aefe9ab531 100644 --- a/TRD/AliTRDmcmSim.cxx +++ b/TRD/AliTRDmcmSim.cxx @@ -13,157 +13,159 @@ * 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 "AliTRDdigitsManager.h" #include "AliTRDarrayADC.h" -#include "AliTRDpadPlane.h" +#include "AliTRDarrayDictionary.h" #include "AliTRDtrackletMCM.h" #include "AliTRDmcmSim.h" ClassImp(AliTRDmcmSim) -//_____________________________________________________________________________ -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) - ,fCal(NULL) - ,fGeo(NULL) - ,fPedAcc(NULL) - ,fGainCounterA(NULL) - ,fGainCounterB(NULL) - ,fTailAmplLong(NULL) - ,fTailAmplShort(NULL) - ,fNHits(0) - ,fFitReg(NULL) +Bool_t AliTRDmcmSim::fgApplyCut = kTRUE; +Int_t AliTRDmcmSim::fgAddBaseline = 0; + +const Int_t AliTRDmcmSim::fgkFormatIndex = std::ios_base::xalloc(); + +const Int_t AliTRDmcmSim::fgkNADC = AliTRDfeeParam::GetNadcMcm(); +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) { // // 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 < fgkNADC; 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(); - fCal = AliTRDcalibDB::Instance(); - fGeo = new AliTRDgeometry(); } 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 *[fgkNADC]; + fADCF = new Int_t *[fgkNADC]; + fZSMap = new Int_t [fgkNADC]; + fGainCounterA = new UInt_t[fgkNADC]; + fGainCounterB = new UInt_t[fgkNADC]; + fNTimeBin = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kC13CPUA, fDetector, fRobPos, fMcmPos); + for( Int_t iAdc = 0 ; iAdc < fgkNADC; 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[fgkNADC]; // accumulator for pedestal filter + fTailAmplLong = new UShort_t[fgkNADC]; + fTailAmplShort = new UShort_t[fgkNADC]; + // tracklet calculation - fFitReg = new FitReg_t[fNADC]; - fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets); - - fMCMT = new UInt_t[fMaxTracklets]; + fFitReg = new FitReg_t[fgkNADC]; + fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fgkMaxTracklets); + + fMCMT = new UInt_t[fgkMaxTracklets]; } fInitialized = kTRUE; @@ -176,32 +178,55 @@ 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 < fgkNADC; 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 < fgkNADC; iAdc++ ) { + delete [] fADCR[iAdc]; + delete [] fADCF[iAdc]; + fADCR[iAdc] = new Int_t[fNTimeBin]; + fADCF[iAdc] = new Int_t[fNTimeBin]; + } } -Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader *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. - if (!CheckInitialized()) - Init(det, rob, mcm); + Init(det, rob, mcm); if (!runloader) { AliError("No Runloader given"); @@ -214,98 +239,129 @@ Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader *runloader, Int_t det, Int_t rob, Int_ return kFALSE; } + Bool_t retval = kTRUE; trdLoader->LoadDigits(); + fDigitsManager = 0x0; AliTRDdigitsManager *digMgr = new AliTRDdigitsManager(); digMgr->SetSDigits(0); 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++) { - for (Int_t tb = 0; tb < fNTimeBin; tb++) { - padcol = fFeeParam->GetPadColFromADC(rob, mcm, ch); - 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); } - digMgr->RemoveDigits(det); + 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(); @@ -317,58 +373,45 @@ 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* option) const +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("U")) { - 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")) { printf("Found %i hits:\n", fNHits); for (Int_t iHit = 0; iHit < fNHits; iHit++) { printf("Hit %3i in timebin %2i, ADC %2i has charge %3i and position %3i\n", - iHit, fHits[iHit].timebin, fHits[iHit].channel, fHits[iHit].qtot, fHits[iHit].ypos); + iHit, fHits[iHit].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos); } } @@ -380,36 +423,39 @@ void AliTRDmcmSim::Print(Option_t* option) const } } -void AliTRDmcmSim::Draw(Option_t* 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); + fgkNADC, -0.5, fgkNADC-.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 < fgkNADC; 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 < fgkNADC; iAdc++) { hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits); } } @@ -419,25 +465,29 @@ void AliTRDmcmSim::Draw(Option_t* option) if (opt.Contains("H")) { TGraph *grHits = new TGraph(); for (Int_t iHit = 0; iHit < fNHits; iHit++) { - grHits->SetPoint(iHit, - fHits[iHit].channel + 1 + fHits[iHit].ypos/256., - fHits[iHit].timebin); + grHits->SetPoint(iHit, + fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256., + fHits[iHit].fTimebin); } grHits->Draw("*"); } if (opt.Contains("T")) { TLine *trklLines = new TLine[4]; - for (Int_t iTrkl = 0; iTrkl < 4; iTrkl++) { - if (fMCMT[iTrkl] == 0x10001000) - break; - AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector); + for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) { 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(AliTRDtrapConfig::fgkDmemAddrNdrift, fDetector, fRobPos, fMcmPos) >> 5; + Float_t 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); @@ -446,7 +496,7 @@ void AliTRDmcmSim::Draw(Option_t* option) } } -void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc ) +void AliTRDmcmSim::SetData( Int_t adc, Int_t* const data ) { // // Store ADC data into array of raw data @@ -454,18 +504,18 @@ void AliTRDmcmSim::SetData( Int_t iadc, 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 >= fgkNADC ) { + AliError(Form ("Error: ADC %i is out of range (0 .. %d).", adc, fgkNADC-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 @@ -473,140 +523,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 >= fgkNADC ) { + AliError(Form ("Error: ADC %i is out of range (0 .. %d).", adc, fgkNADC-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 *adcArray) +void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager * const digitsManager) { - if (!fInitialized) { - AliError("Called uninitialized! Aborting!"); + // Set the ADC data from an AliTRDarrayADC + + 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; + 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; + } + } } - Int_t firstAdc = 0; - Int_t lastAdc = fNADC-1; + if (fNTimeBin != adcArray->GetNtime()) + SetNTimebins(adcArray->GetNtime()); - if (GetCol(firstAdc) > 143) { - for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - fADCR[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; - fADCF[firstAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; + 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 < fgkNADC; 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); + } } - firstAdc = 1; } +} - if (GetCol(lastAdc) < 0) { - for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { - fADCR[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; - fADCF[lastAdc][iTimeBin] = fSimParam->GetADCbaseline() << fgkAddDigits; +void AliTRDmcmSim::SetDataByPad(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 = fNADC - 2; } + 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 < fgkNADC; 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 >= fgkNADC ) { 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() ) return -1; + if( !CheckInitialized() ) + return -1; - return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc); + Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adc); + if (col < 0 || col >= fFeeParam->GetNcol()) + return -1; + 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 < fgkNADC ; 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 @@ -615,16 +771,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++; @@ -635,35 +790,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 // - UInt_t x; + 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 - Int_t wd = 0; - while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){ - x = fMCMT[wd]; - if (nw < maxSize) { - buf[nw++] = x; - } - else { - of++; - } - wd++; + for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) { + if (nw < bufSize) + buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord(); + else + of++; } - + if( of != 0 ) return -of; else return nw; } @@ -672,22 +822,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(); @@ -695,54 +843,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; } } @@ -752,14 +897,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 < fgkNADC; iAdc++) { fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]); } } @@ -767,11 +912,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 < fgkNADC; iAdc++) { + // these are counters which in hardware continue // until maximum or reset fGainCounterA[iAdc] = 0; fGainCounterB[iAdc] = 0; @@ -782,42 +927,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 < fgkNADC; iAdc++) { for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) { fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]); } @@ -826,14 +977,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); @@ -845,13 +996,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 < fgkNADC; 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)); } @@ -859,67 +1019,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 + 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 - 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; - - UInt_t aDiff; - UInt_t alInpv; + // intermediate signals + UInt_t aDiff; + UInt_t alInpv; UShort_t aQ; - UInt_t tmp; - - UShort_t inp_volt = value & 0xFFF; // 12 bits - - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low + UInt_t tmp; + + UShort_t inpVolt = value & 0xFFF; // 12 bits + + // 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 (inp_volt > aQ) - aDiff = inp_volt - 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 < fgkNADC; iAdc++) { fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]); } } @@ -929,239 +1080,198 @@ 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 = AliTRDfeeParam::GetPFeffectPedestal(); + 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! Aborting!"); - return; - } + for (Int_t iAdc = 0; iAdc < fgkNADC; 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) - // 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]; + // ----- 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); } - } -} + // ----- last channel ----- + iAdc = fgkNADC - 1; -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 - // + ap = adc[iAdc-1][it]; // previous + ac = adc[iAdc ][it]; // current + an = 0; // 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); + } + + // ----- middle channels ----- + for( iAdc = 1 ; iAdc < fgkNADC-1; iAdc++ ) { + ap = adc[iAdc-1][it]; // 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); + 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))) - fFitReg[adc].Q0 += qtot; - - if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) && - (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))) - fFitReg[adc].Q1 += qtot; - - if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) && - (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE))) + 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, fDetector, fRobPos, fMcmPos)) && + (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1, fDetector, fRobPos, fMcmPos))) + fFitReg[adc].fQ1 += qtot; + + if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS, fDetector, fRobPos, fMcmPos) ) && + (timebin < fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE, fDetector, fRobPos, fMcmPos))) { - fFitReg[adc].SumX += timebin; - fFitReg[adc].SumX2 += timebin*timebin; - fFitReg[adc].Nhits++; - fFitReg[adc].SumY += ypos; - fFitReg[adc].SumY2 += ypos*ypos; - fFitReg[adc].SumXY += timebin*ypos; + fFitReg[adc].fSumX += timebin; + fFitReg[adc].fSumX2 += timebin*timebin; + fFitReg[adc].fNhits++; + 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) - fHits[fNHits].channel = adc; - fHits[fNHits].qtot = qtot; - fHits[fNHits].ypos = ypos; - fHits[fNHits].timebin = timebin; - fHits[fNHits].label = label; + fHits[fNHits].fChannel = adc; + fHits[fNHits].fQtot = qtot; + fHits[fNHits].fYpos = ypos; + fHits[fNHits].fTimebin = timebin; + 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 uint16_t LUT_POS[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 adc_mask = 0xfffff; - - UShort_t timebin, adcch, adc_left, adc_cent, adc_right, hit_qual, timebin1, timebin2, Qtot_tmp; + UInt_t adcMask = 0xffffffff; + + UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp; Short_t ypos, fromLeft, fromRight, found; - UShort_t Qtot[19]; // the last is dummy - UShort_t marked[6], Qmarked[6], worse1, worse2; - - timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS); - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) + UShort_t qTotal[19+1]; // the last is dummy + UShort_t marked[6], qMarked[6], worse1, worse2; + + timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS, fDetector, fRobPos, fMcmPos); + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0, fDetector, fRobPos, fMcmPos) < timebin1) - timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0); - timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE); - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1) + timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0, fDetector, fRobPos, fMcmPos); + timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE, fDetector, fRobPos, fMcmPos); + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1, fDetector, fRobPos, fMcmPos) > timebin2) - timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1); + 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 < fgkNADC-2; adcch++) // due to border channels { - fFitReg[adcch].Nhits = 0; - fFitReg[adcch].Q0 = 0; - fFitReg[adcch].Q1 = 0; - fFitReg[adcch].SumX = 0; - fFitReg[adcch].SumY = 0; - fFitReg[adcch].SumX2 = 0; - fFitReg[adcch].SumY2 = 0; - fFitReg[adcch].SumXY = 0; + fFitReg[adcch].fNhits = 0; + fFitReg[adcch].fQ0 = 0; + fFitReg[adcch].fQ1 = 0; + fFitReg[adcch].fSumX = 0; + fFitReg[adcch].fSumY = 0; + fFitReg[adcch].fSumX2 = 0; + 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 Qtot array + // 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++) { - if ( ( (adc_mask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS + for (adcch = 0; adcch < fgkNADC-2; adcch++) { + if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS { - adc_left = fADCF[adcch ][timebin]; - adc_cent = fADCF[adcch+1][timebin]; - adc_right = fADCF[adcch+2][timebin]; - if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1) - hit_qual = ( (adc_left * adc_right) < - (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adc_cent) ); - else - hit_qual = 1; + adcLeft = fADCF[adcch ][timebin]; + adcCentral = fADCF[adcch+1][timebin]; + adcRight = fADCF[adcch+2][timebin]; + if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY, fDetector, fRobPos, fMcmPos) == 1) + hitQual = ( (adcLeft * adcRight) < + (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT, fDetector, fRobPos, fMcmPos) * adcCentral) ); + else + hitQual = 1; // The accumulated charge is with the pedestal!!! - Qtot_tmp = adc_left + adc_cent + adc_right; - if ( (hit_qual) && - (Qtot_tmp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) && - (adc_left <= adc_cent) && - (adc_cent > adc_right) ) - Qtot[adcch] = Qtot_tmp; + qtotTemp = adcLeft + adcCentral + adcRight; + if ( (hitQual) && + (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT, fDetector, fRobPos, fMcmPos)) && + (adcLeft <= adcCentral) && + (adcCentral > adcRight) ) + qTotal[adcch] = qtotTemp; else - Qtot[adcch] = 0; - //printf("ch %2d Qtot %5d\n",adcch, Qtot[adcch]); + qTotal[adcch] = 0; } else - Qtot[adcch] = 0; //jkl + qTotal[adcch] = 0; //jkl + if (qTotal[adcch] != 0) + AliDebug(10,Form("ch %2d qTotal %5d",adcch, qTotal[adcch])); } fromLeft = -1; @@ -1169,10 +1279,10 @@ void AliTRDmcmSim::CalcFitreg() found = 0; marked[4] = 19; // invalid channel marked[5] = 19; // invalid channel - Qtot[19] = 0; + qTotal[19] = 0; while ((adcch < 16) && (found < 3)) { - if (Qtot[adcch] > 0) + if (qTotal[adcch] > 0) { fromLeft = adcch; marked[2*found+1]=adcch; @@ -1180,13 +1290,13 @@ void AliTRDmcmSim::CalcFitreg() } adcch++; } - + fromRight = -1; adcch = 18; found = 0; while ((adcch > 2) && (found < 3)) { - if (Qtot[adcch] > 0) + if (qTotal[adcch] > 0) { marked[2*found]=adcch; found++; @@ -1195,15 +1305,15 @@ 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++) - Qtot[adcch] = 0; - + qTotal[adcch] = 0; + found = 0; for (adcch = 0; adcch < 19; adcch++) - if (Qtot[adcch] > 0) found++; + if (qTotal[adcch] > 0) found++; // NOT READY if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates! @@ -1211,83 +1321,146 @@ void AliTRDmcmSim::CalcFitreg() if (marked[4] == marked[5]) marked[5] = 19; for (found=0; found<6; found++) { - Qmarked[found] = Qtot[marked[found]] >> 4; - //printf("ch_%d Qtot %d Qtots %d |",marked[found],Qtot[marked[found]],Qmarked[found]); + qMarked[found] = qTotal[marked[found]] >> 4; + 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], - Qmarked[4], - Qmarked[1], - Qmarked[2], - Qmarked[5], + qMarked[0], + qMarked[3], + qMarked[4], + qMarked[1], + qMarked[2], + qMarked[5], &worse1, &worse2); // Now mask the two channels with the smallest charge if (worse1 < 19) { - Qtot[worse1] = 0; - //printf("Kill ch %d\n",worse1); + qTotal[worse1] = 0; + AliDebug(10,Form("Kill ch %d\n",worse1)); } if (worse2 < 19) { - Qtot[worse2] = 0; - //printf("Kill ch %d\n",worse2); + qTotal[worse2] = 0; + AliDebug(10,Form("Kill ch %d\n",worse2)); } } - + for (adcch = 0; adcch < 19; adcch++) { - if (Qtot[adcch] > 0) // the channel is marked for processing + if (qTotal[adcch] > 0) // the channel is marked for processing { - adc_left = fADCF[adcch ][timebin]; - adc_cent = fADCF[adcch+1][timebin]; - adc_right = fADCF[adcch+2][timebin]; + adcLeft = fADCF[adcch ][timebin]; + adcCentral = fADCF[adcch+1][timebin]; + 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 TPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP); -// printf("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, TPFP=%d, TPHT=%d\n", -// timebin, adcch, adcch+1, adcch+2, adc_left, adc_cent, adc_right, TPFP, -// fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)); - - if (adc_left < TPFP) adc_left = 0; else adc_left -= TPFP; - if (adc_cent < TPFP) adc_cent = 0; else adc_cent -= TPFP; - if (adc_right < TPFP) adc_right = 0; else adc_right -= TPFP; + + 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; + if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP; + // Calculate the center of gravity - ypos = 128*(adc_left - adc_right) / adc_cent; + // checking for adcCentral != 0 (in case of "bad" configuration) + if (adcCentral == 0) + continue; + ypos = 128*(adcRight - adcLeft) / adcCentral; if (ypos < 0) ypos = -ypos; - // make the correction using the LUT - ypos = ypos + LUT_POS[ypos & 0x7F]; - if (adc_left > adc_right) ypos = -ypos; - AddHitToFitreg(adcch, timebin, Qtot[adcch], ypos, -1); + // 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 (up to 3) + Int_t mcLabel[] = {-1, -1, -1}; + if (fDigitsManager) { + 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]; + 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 (!fDict[iDict]) + continue; + for (Int_t iPad = 0; iPad < 3; iPad++) { + if (padcol[iPad] < 0) + continue; + 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]++; + currLabel = -1; + break; + } + } + if (currLabel >= 0) { + label[nLabels] = currLabel; + count[nLabels] = 1; + nLabels++; + } + } + } + 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] >> fgkAddDigits, ypos, mcLabel); } } } + + for (Int_t iAdc = 0; iAdc < fgkNADC; 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 adc_idx, i, j, ntracks, tmp; - UShort_t track_p[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates + UShort_t adcIdx, i, j, ntracks, tmp; + UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates ntracks = 0; - for (adc_idx = 0; adc_idx < 18; adc_idx++) // ADCs - if ( (fFitReg[adc_idx].Nhits - >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) && - (fFitReg[adc_idx].Nhits+fFitReg[adc_idx+1].Nhits - >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT))) + for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs + if ( (fFitReg[adcIdx].fNhits + >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL, fDetector, fRobPos, fMcmPos)) && + (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits + >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT, fDetector, fRobPos, fMcmPos))) { - track_p[ntracks][0] = adc_idx; - track_p[ntracks][1] = fFitReg[adc_idx].Nhits+fFitReg[adc_idx+1].Nhits; - //printf("%d %2d %4d\n", ntracks, track_p[ntracks][0], track_p[ntracks][1]); + trackletCand[ntracks][0] = adcIdx; + trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits; + AliDebug(10,Form("%d %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1])); ntracks++; }; - // for (i=0; i 4) { @@ -1296,81 +1469,87 @@ void AliTRDmcmSim::TrackletSelection() { for (i = j+1; i < ntracks; i++) { - if ( (track_p[j][1] < track_p[i][1]) || - ( (track_p[j][1] == track_p[i][1]) && (track_p[j][0] < track_p[i][0]) ) ) + if ( (trackletCand[j][1] < trackletCand[i][1]) || + ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) ) { // swap j & i - tmp = track_p[j][1]; - track_p[j][1] = track_p[i][1]; - track_p[i][1] = tmp; - tmp = track_p[j][0]; - track_p[j][0] = track_p[i][0]; - track_p[i][0] = tmp; + tmp = trackletCand[j][1]; + trackletCand[j][1] = trackletCand[i][1]; + trackletCand[i][1] = tmp; + tmp = trackletCand[j][0]; + trackletCand[j][0] = trackletCand[i][0]; + trackletCand[i][0] = tmp; } } } 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++) { for (i = j+1; i < ntracks; i++) { - if (track_p[j][0] < track_p[i][0]) + if (trackletCand[j][0] < trackletCand[i][0]) { // swap j & i - tmp = track_p[j][1]; - track_p[j][1] = track_p[i][1]; - track_p[i][1] = tmp; - tmp = track_p[j][0]; - track_p[j][0] = track_p[i][0]; - track_p[i][0] = tmp; + tmp = trackletCand[j][1]; + trackletCand[j][1] = trackletCand[i][1]; + trackletCand[i][1] = tmp; + tmp = trackletCand[j][0]; + trackletCand[j][0] = trackletCand[i][0]; + trackletCand[i][0] = tmp; } } } for (i = 0; i < ntracks; i++) // CPUs with tracklets. - fFitPtr[i] = track_p[i][0]; // pointer to the left channel with tracklet for CPU[i] + 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); + 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 scale_y = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4))); - UInt_t scale_d = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4))); - 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 defl_cor = 0; // -370; - int minslope = -10000; // no pt-cut so far - int maxslope = 10000; // no pt-cut so far + + // 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(AliTRDtrapConfig::fgkDmemAddrDeflCorr, fDetector, fRobPos, fMcmPos); + Int_t ndrift = (Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::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, @@ -1380,7 +1559,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 { @@ -1392,62 +1571,181 @@ void AliTRDmcmSim::FitTracklet() mult = -mult; // Merging - nHits = fit0->Nhits + fit1->Nhits; // number of hits - sumX = fit0->SumX + fit1->SumX; - sumX2 = fit0->SumX2 + fit1->SumX2; - denom = nHits*sumX2 - sumX*sumX; + nHits = fit0->fNhits + fit1->fNhits; // number of hits + sumX = fit0->fSumX + fit1->fSumX; + sumX2 = fit0->fSumX2 + fit1->fSumX2; + denom = ((Long64_t) nHits)*((Long64_t) sumX2) - ((Long64_t) sumX)*((Long64_t) sumX); mult = mult / denom; // exactly like in the TRAP program - q0 = fit0->Q0 + fit1->Q0; - q1 = fit0->Q1 + fit1->Q1; - sumY = fit0->SumY + fit1->SumY + 256*fit1->Nhits; - sumXY = fit0->SumXY + fit1->SumXY + 256*fit1->SumX; + 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; 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)); - slope = slope * ndrift + defl_cor; + 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)); - - if ((slope < minslope) || (slope > 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(AliTRDtrapConfig::fgkDmemAddrDeflCutStart + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos), + (Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCutStart + 1 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos))); + + AliDebug(5, Form("Fit sums: x = %i, X = %i, y = %i, Y = %i, Z = %i", + sumX, sumX2, sumY, sumY2, sumXY)); + + 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; + // deflection range table from DMEM + if ((slope < ((Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCutStart + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos))) || + (slope > ((Int_t) fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrDeflCutStart + 1 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)))) + rejected = kTRUE; + + if (rejected && GetApplyCut()) { fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker(); } else { - temp = slope; - temp = temp * scale_d; - slope = (temp >> 32); - - temp = offset; - temp = temp * scale_y; - offset = (temp >> 32); - - // rounding, like in the TRAP - slope = (slope + rndAdd) >> decPlaces; - offset = (offset + rndAdd) >> decPlaces; - - if (slope > 0x3f || slope < -0x3f) - AliWarning("Overflow in slope"); + if (slope > 63 || slope < -64) { // wrapping in TRAP! + AliDebug(1,Form("Overflow in slope: %i, tracklet discarded!", slope)); + fMCMT[cpu] = 0x10001000; + continue; + } + slope = slope & 0x7F; // 7 bit - if (offset > 0xfff || offset < 0xfff) + 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); + + if (pid > 0xff) + AliWarning("Overflow in PID"); + pid = pid & 0xFF; // 8 bit, exactly like in the TRAP program // assemble and store the tracklet word - fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset; - new ((*fTrackletArray)[cpu]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos); + fMCMT[cpu] = (pid << 24) | (padrow << 20) | (slope << 13) | offset; + + // calculate MC label + Int_t mcLabel[] = { -1, -1, -1}; + Int_t nHits0 = 0; + Int_t nHits1 = 0; + if (fDigitsManager) { + 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++; + + 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++; + } + } + } + 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)); + +// // cluster information +// Float_t *res = new Float_t[nHits]; +// Float_t *qtot = new Float_t[nHits]; +// Int_t nCls = 0; +// for (Int_t iHit = 0; iHit < fNHits; iHit++) { +// // check if hit contributes +// if (fHits[iHit].fChannel == fFitPtr[cpu]) { +// res[nCls] = fHits[iHit].fYpos - (fitSlope * fHits[iHit].fTimebin + fitOffset); +// qtot[nCls] = fHits[iHit].fQtot; +// nCls++; +// } +// else if (fHits[iHit].fChannel == fFitPtr[cpu] + 1) { +// res[nCls] = fHits[iHit].fYpos + 256 - (fitSlope * fHits[iHit].fTimebin + fitOffset); +// qtot[nCls] = fHits[iHit].fQtot; +// nCls++; +// } +// } +// ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetClusters(res, qtot, nCls); +// 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))); } } } @@ -1455,16 +1753,30 @@ void AliTRDmcmSim::FitTracklet() void AliTRDmcmSim::Tracklet() { + // Run the tracklet calculation by calling sequentially: + // CalcFitreg(); TrackletSelection(); FitTracklet() + // and store the tracklets + if (!fInitialized) { - AliError("Called uninitialized! Aborting!"); + AliError("Called uninitialized! Nothing done!"); return; } fTrackletArray->Delete(); CalcFitreg(); + if (fNHits == 0) + return; TrackletSelection(); FitTracklet(); +} + +Bool_t AliTRDmcmSim::StoreTracklets() +{ + // store the found tracklets via the loader + + if (fTrackletArray->GetEntriesFast() == 0) + return kTRUE; AliRunLoader *rl = AliRunLoader::Instance(); AliDataLoader *dl = 0x0; @@ -1472,27 +1784,27 @@ void AliTRDmcmSim::Tracklet() dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets"); if (!dl) { AliError("Could not get the tracklets data loader!"); + return kFALSE; } - else { - TTree *trackletTree = dl->Tree(); - if (!trackletTree) - dl->MakeTree(); + + TTree *trackletTree = dl->Tree(); + if (!trackletTree) { + dl->MakeTree(); trackletTree = dl->Tree(); + } - AliTRDtrackletMCM *trkl = 0x0; - TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch"); - if (!trkbranch) - trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000); -// trkbranch = trackletTree->Branch("mcmtrklbranch", &fTrackletArray, 32000, 2); - - 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"); + AliTRDtrackletMCM *trkl = 0x0; + TBranch *trkbranch = trackletTree->GetBranch(fTrklBranchName.Data()); + if (!trkbranch) + trkbranch = trackletTree->Branch(fTrklBranchName.Data(), "AliTRDtrackletMCM", &trkl, 32000); + + for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) { + trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]); + trkbranch->SetAddress(&trkl); + trkbranch->Fill(); } + + return kTRUE; } void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits) @@ -1501,78 +1813,147 @@ 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! Aborting!"); + if( !CheckInitialized() ) return; - } - - Int_t firstAdc = 0; - Int_t lastAdc = fNADC - 1; - if (GetCol(firstAdc) > 143) - firstAdc = 1; + Int_t offset = (fMcmPos % 4 + 1) * 21 + (fRobPos % 2) * 84 - 1; - if (GetCol(lastAdc) < 0) - lastAdc = fNADC - 2; - - 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 < fgkNADC; 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 < fgkNADC; 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(AliTRDtrapConfig::fgkDmemAddrLUTnbins); // number of bins in q0 / 4 !! + UInt_t pidTotalSize = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTLength); + 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(AliTRDtrapConfig::fgkDmemAddrLUTcor0, fDetector, fRobPos, fMcmPos); + ULong_t corrQ1 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::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(AliTRDtrapConfig::fgkDmemAddrLUTStart+(addr/4)); + 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) +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; } -void AliTRDmcmSim::Sort2(uint16_t idx1i, uint16_t idx2i, \ - uint16_t val1i, uint16_t val2i, \ - uint16_t *idx1o, uint16_t *idx2o, \ - uint16_t *val1o, uint16_t *val2o) +void AliTRDmcmSim::Sort2(UShort_t idx1i, UShort_t idx2i, \ + UShort_t val1i, UShort_t val2i, \ + UShort_t * const idx1o, UShort_t * const idx2o, \ + UShort_t * const val1o, UShort_t * const val2o) const { + // sorting for tracklet selection if (val1i > val2i) { @@ -1590,18 +1971,19 @@ void AliTRDmcmSim::Sort2(uint16_t idx1i, uint16_t idx2i, \ } } -void AliTRDmcmSim::Sort3(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, \ - uint16_t val1i, uint16_t val2i, uint16_t val3i, \ - uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, \ - uint16_t *val1o, uint16_t *val2o, uint16_t *val3o) +void AliTRDmcmSim::Sort3(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, \ + UShort_t val1i, UShort_t val2i, UShort_t val3i, \ + UShort_t * const idx1o, UShort_t * const idx2o, UShort_t * const idx3o, \ + UShort_t * const val1o, UShort_t * const val2o, UShort_t * const val3o) { - int sel; + // sorting for tracklet selection + + 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 @@ -1660,22 +2042,22 @@ void AliTRDmcmSim::Sort3(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, \ break; default: // the rest should NEVER happen! - printf("ERROR in Sort3!!!\n"); + 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(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \ - uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \ - uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, uint16_t *idx4o, \ - uint16_t *val1o, uint16_t *val2o, uint16_t *val3o, uint16_t *val4o) +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 * 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 - uint16_t idx21s, idx22s, idx23s, dummy; - uint16_t val21s, val22s, val23s; - uint16_t idx23as, idx23bs; - uint16_t val23as, val23bs; + UShort_t idx21s, idx22s, idx23s, dummy; + UShort_t val21s, val22s, val23s; + UShort_t idx23as, idx23bs; + UShort_t val23as, val23bs; Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i, idx1o, &idx21s, &idx23as, @@ -1693,15 +2075,16 @@ void AliTRDmcmSim::Sort6To4(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, u } -void AliTRDmcmSim::Sort6To2Worst(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \ - uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \ - uint16_t *idx5o, uint16_t *idx6o) +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 * const idx5o, UShort_t * const idx6o) { + // sorting for tracklet selection - uint16_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5; - uint16_t val21s, val22s, val23s; - uint16_t idx23as, idx23bs; - uint16_t val23as, val23bs; + UShort_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5; + UShort_t val21s, val22s, val23s; + UShort_t idx23as, idx23bs; + UShort_t val23as, val23bs; Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i, &dummy1, &idx21s, &idx23as, @@ -1716,7 +2099,318 @@ void AliTRDmcmSim::Sort6To2Worst(uint16_t idx1i, uint16_t idx2i, uint16_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 < mcm.fgkNADC; 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 < mcm.fgkNADC; 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 < mcm.fgkNADC; 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 < mcm.fgkNADC; 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 < mcm.fgkNADC; 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/4 << ""<< std::endl; // divided by 4 because in simulation we have 2 additional decimal places + os << " " << fFitReg[adcch].fQ1/4 << ""<< std::endl; // in the output + 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::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 < fgkNADC; 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 < fgkNADC; 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 < fgkNADC; 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 < fgkNADC; 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 < fgkNADC; 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 < fgkNADC; 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 = AliTRDtrapConfig::fgkDmemAddrLUTStart + fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTLength)/4; // /4 because each addr contains 4 values + UInt_t nBinsQ0 = fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTnbins); + + std::cout << "nBinsQ0: " << nBinsQ0 << std::endl; + std::cout << "LUT table length: " << fTrapConfig->GetDmemUnsigned(AliTRDtrapConfig::fgkDmemAddrLUTLength) << std::endl; + + for(UInt_t addr=AliTRDtrapConfig::fgkDmemAddrLUTStart; addr< addrEnd; addr++) { + result = fTrapConfig->GetDmemUnsigned(addr); + std::cout << addr << " # x: " << ((addr-AliTRDtrapConfig::fgkDmemAddrLUTStart)%((nBinsQ0)/4))*4 << ", y: " <<(addr-AliTRDtrapConfig::fgkDmemAddrLUTStart)/(nBinsQ0/4) + << " # " <<((result>>0)&0xFF) + << " | " << ((result>>8)&0xFF) + << " | " << ((result>>16)&0xFF) + << " | " << ((result>>24)&0xFF) << std::endl; + } +}