* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
+
///////////////////////////////////////////////////////////////////////////////
// //
// TRD MCM (Multi Chip Module) simulator //
+// which simulates the TRAP processing after the AD-conversion. //
+// The relevant parameters (i.e. configuration settings of the TRAP) //
+// are taken from AliTRDtrapConfig. //
// //
///////////////////////////////////////////////////////////////////////////////
-/* $Id$ */
-
-/*
-
- New release on 2007/08/17
-
-AliTRDmcmSim is now stably working and zero suppression function seems ok.
-From now, the default version of raw data is set to 3 in AliTRDfeeParam.
-
-The following internal parameters were abolished because it is useless and
-made trouble:
-
- fColOfADCbeg
- fColOfADCend
-
-GetCol member was modified accordingly.
-
-New member function DumpData was prepared for diagnostics.
-
-ZSMapping member function was debugged. It was causing crash due to
-wrong indexing in 1 dimensional numbering. Also code was shaped up better.
-
-*/
-
-/*Semi-final version of TRD raw data simulation code with zero suppression (ZS)
-similar to TRD FEE. ZS is realized by the class group:
-
- AliTRDfeeParam
- AliTRDmcmSim
- AliTRDrawData
-
-AliTRDfeeParam has been modified to have more parameters like raw data
-production version and so on. AliTRDmcmSim is new class and this is the core
-of MCM (PASA+TRAP) simulator. It has still very simple function and it will be
-another project to improve this to make it closer to the reall FEE.
-AliTRDrawData has been modified to use new class AliTRDmcmSim.
-
-These modifications were tested on Aug. 02 HEAD version that code itself
-compiles. I'm sure there must be still bugs and we need testing by as many as
-possible persons now. Especially it seems HLT part is impacted by problems
-because some parameters were moved from AliTRDrawData to AliTRDfeeParam (like
-fRawVersion disappeared from AliTRDrawData).
-
-In TRD definition, we have now 4 raw data versions.
-
- 0 very old offline version (by Bogdan)
- 1 test version (no zero suppression)
- 2 final version (no zero suppression)
- 3 test version (with zero suppression)
-
-The default is still set to 2 in AliTRDfeeParam::fgkRAWversion and it uses
-previously existing codes. If you set this to 3, AliTRDrawData changes behavior
-to use AliTRDmcmSim with ZS.
-
-Plan is after we make sure it works stably, we delete AliTRDmcm which is obsolete.
-However it still take time because tracklet part is not yet touched.
-The default raw version is 2.
-
- Ken Oyama
-*/
-
-// if no histo is drawn, these are obsolete
-#include <TH1.h>
-#include <TCanvas.h>
-
-// only needed if I/O of tracklets is activated
-#include <TObject.h>
-#include <TFile.h>
+#include <iostream>
+#include <iomanip>
+
+#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 <TTree.h>
-#include <TSystem.h>
-
-#include <fstream>
-
-#include <TMath.h>
#include "AliLog.h"
+#include "AliRunLoader.h"
+#include "AliLoader.h"
-#include "AliTRDmcmSim.h"
#include "AliTRDfeeParam.h"
-#include "AliTRDSimParam.h"
-#include "AliTRDCommonParam.h"
-#include "AliTRDgeometry.h"
-#include "AliTRDcalibDB.h"
+#include "AliTRDtrapConfig.h"
#include "AliTRDdigitsManager.h"
#include "AliTRDarrayADC.h"
-// additional for new tail filter and/or tracklet
-#include "AliTRDtrapAlu.h"
-#include "AliTRDpadPlane.h"
+#include "AliTRDarrayDictionary.h"
#include "AliTRDtrackletMCM.h"
-
-#include "AliRun.h"
-#include "AliLoader.h"
+#include "AliTRDmcmSim.h"
ClassImp(AliTRDmcmSim)
-//_____________________________________________________________________________
-AliTRDmcmSim::AliTRDmcmSim() :TObject()
- ,fInitialized(kFALSE)
- ,fNextEvent(-1)
- ,fMaxTracklets(-1)
- ,fChaId(-1)
- ,fSector(-1)
- ,fStack(-1)
- ,fLayer(-1)
- ,fRobPos(-1)
- ,fMcmPos(-1)
- ,fNADC(-1)
- ,fNTimeBin(-1)
- ,fRow (-1)
- ,fADCR(NULL)
- ,fADCF(NULL)
- ,fADCT(NULL)
- ,fPosLUT(NULL)
- ,fMCMT(NULL)
- ,fZSM(NULL)
- ,fZSM1Dim(NULL)
- ,fFeeParam(NULL)
- ,fSimParam(NULL)
- ,fCal(NULL)
- ,fGeo(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.
-}
-//_____________________________________________________________________________
-AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
- :TObject(m)
- ,fInitialized(kFALSE)
- ,fNextEvent(-1)
- ,fMaxTracklets(-1)
- ,fChaId(-1)
- ,fSector(-1)
- ,fStack(-1)
- ,fLayer(-1)
- ,fRobPos(-1)
- ,fMcmPos(-1)
- ,fNADC(-1)
- ,fNTimeBin(-1)
- ,fRow(-1)
- ,fADCR(NULL)
- ,fADCF(NULL)
- ,fADCT(NULL)
- ,fPosLUT(NULL)
- ,fMCMT(NULL)
- ,fZSM(NULL)
- ,fZSM1Dim(NULL)
- ,fFeeParam(NULL)
- ,fSimParam(NULL)
- ,fCal(NULL)
- ,fGeo(NULL)
-
-{
- //
- // AliTRDmcmSim copy constructor
- //
+ for (Int_t iDict = 0; iDict < 3; iDict++)
+ fDict[iDict] = 0x0;
- // By default, nothing is initialized.
- // It is necessary to issue Init before use.
+ fFitPtr[0] = 0;
+ fFitPtr[1] = 0;
+ fFitPtr[2] = 0;
+ fFitPtr[3] = 0;
}
-//_____________________________________________________________________________
-AliTRDmcmSim::~AliTRDmcmSim()
+AliTRDmcmSim::~AliTRDmcmSim()
{
//
// AliTRDmcmSim destructor
//
- if( fADCR != NULL ) {
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- delete [] fADCR[iadc];
- delete [] fADCF[iadc];
- delete [] fADCT[iadc];
- delete [] fZSM [iadc];
+ if(fInitialized) {
+ for( Int_t iAdc = 0 ; iAdc < fgkNADC; iAdc++ ) {
+ delete [] fADCR[iAdc];
+ delete [] fADCF[iAdc];
}
delete [] fADCR;
delete [] fADCF;
- delete [] fADCT;
- delete [] fZSM;
- delete [] fZSM1Dim;
- }
-
- if(fInitialized){
- delete [] fPosLUT;
+ delete [] fZSMap;
delete [] fMCMT;
- }
- delete fGeo;
+ delete [] fPedAcc;
+ delete [] fGainCounterA;
+ delete [] fGainCounterB;
+ delete [] fTailAmplLong;
+ delete [] fTailAmplShort;
+ delete [] fFitReg;
+ fTrackletArray->Delete();
+ delete fTrackletArray;
+ }
}
-//_____________________________________________________________________________
-AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m)
+void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ )
{
//
- // Assignment operator
+ // Initialize the class with new MCM position information
+ // memory is allocated in the first initialization
//
- if (this != &m) {
- ((AliTRDmcmSim &) m).Copy(*this);
+ if (!fInitialized) {
+ fFeeParam = AliTRDfeeParam::Instance();
+ fTrapConfig = AliTRDtrapConfig::Instance();
}
- return *this;
-}
+ fDetector = det;
+ fRobPos = robPos;
+ fMcmPos = mcmPos;
+ fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
-//_____________________________________________________________________________
-void AliTRDmcmSim::Copy(TObject &m) const
-{
- //
- // Copy function
- //
- ((AliTRDmcmSim &) m).fNextEvent = 0; //new
- ((AliTRDmcmSim &) m).fMaxTracklets = 0; //new
- ((AliTRDmcmSim &) m).fInitialized = 0;
- ((AliTRDmcmSim &) m).fChaId = 0;
- ((AliTRDmcmSim &) m).fSector = 0;
- ((AliTRDmcmSim &) m).fStack = 0;
- ((AliTRDmcmSim &) m).fLayer = 0;
- ((AliTRDmcmSim &) m).fRobPos = 0;
- ((AliTRDmcmSim &) m).fMcmPos = 0;
- ((AliTRDmcmSim &) m).fNADC = 0;
- ((AliTRDmcmSim &) m).fNTimeBin = 0;
- ((AliTRDmcmSim &) m).fRow = 0;
- ((AliTRDmcmSim &) m).fADCR = 0;
- ((AliTRDmcmSim &) m).fADCF = 0;
- ((AliTRDmcmSim &) m).fADCT = 0; //new
- ((AliTRDmcmSim &) m).fPosLUT = 0; //new
- ((AliTRDmcmSim &) m).fMCMT = 0; //new
- ((AliTRDmcmSim &) m).fZSM = 0;
- ((AliTRDmcmSim &) m).fZSM1Dim = 0;
- ((AliTRDmcmSim &) m).fFeeParam = 0;
- ((AliTRDmcmSim &) m).fSimParam = 0;
- ((AliTRDmcmSim &) m).fCal = 0;
- ((AliTRDmcmSim &) m).fGeo = 0;
+ if (!fInitialized) {
+ 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[fgkNADC]; // accumulator for pedestal filter
+ fTailAmplLong = new UShort_t[fgkNADC];
+ fTailAmplShort = new UShort_t[fgkNADC];
+
+ // tracklet calculation
+ fFitReg = new FitReg_t[fgkNADC];
+ fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fgkMaxTracklets);
+
+ fMCMT = new UInt_t[fgkMaxTracklets];
+ }
+
+ fInitialized = kTRUE;
-//_____________________________________________________________________________
+ Reset();
+}
-//void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
-void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos, Bool_t newEvent = kFALSE ) // only for readout tree (new event)
+void AliTRDmcmSim::Reset()
{
- //
- // Initialize the class with new geometry information
- // fADC array will be reused with filled by zero
- //
-
- fNextEvent = 0;
- fFeeParam = AliTRDfeeParam::Instance();
- fSimParam = AliTRDSimParam::Instance();
- fCal = AliTRDcalibDB::Instance();
- fGeo = new AliTRDgeometry();
- fChaId = chaId;
- fSector = fGeo->GetSector( fChaId );
- fStack = fGeo->GetStack( fChaId );
- fLayer = fGeo->GetLayer( fChaId );
- fRobPos = robPos;
- fMcmPos = mcmPos;
- fNADC = fFeeParam->GetNadcMcm();
- fNTimeBin = fCal->GetNumberOfTimeBins();
- fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
+ // Resets the data values and internal filter registers
+ // by re-initialising them
- fMaxTracklets = fFeeParam->GetMaxNrOfTracklets();
+ 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;
+ }
+ fZSMap[iAdc] = -1; // Default unread, low active bit mask
+ fGainCounterA[iAdc] = 0;
+ fGainCounterB[iAdc] = 0;
+ }
-
- if (newEvent == kTRUE) {
- fNextEvent = 1;
+ 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();
+}
+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)
- // Allocate ADC data memory if not yet done
- if( fADCR == NULL ) {
- fADCR = new Int_t *[fNADC];
- fADCF = new Int_t *[fNADC];
- fADCT = new Int_t *[fNADC]; //new
- fZSM = new Int_t *[fNADC];
- fZSM1Dim = new Int_t [fNADC];
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- fADCR[iadc] = new Int_t[fNTimeBin];
- fADCF[iadc] = new Int_t[fNTimeBin];
- fADCT[iadc] = new Int_t[fNTimeBin]; //new
- fZSM [iadc] = new Int_t[fNTimeBin];
- }
+ 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];
}
+}
- // Initialize ADC data
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
- fADCR[iadc][it] = 0;
- fADCF[iadc][it] = 0;
- fADCT[iadc][it] = -1; //new
- fZSM [iadc][it] = 1; // Default unread = 1
- }
- fZSM1Dim[iadc] = 1; // Default unread = 1
+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.
+ // This method is meant for rare execution, e.g. in the visualization. When called
+ // frequently use SetData(...) instead.
+
+ Init(det, rob, mcm);
+
+ if (!runloader) {
+ AliError("No Runloader given");
+ return kFALSE;
}
-
- //new:
- fPosLUT = new Int_t[128];
- for(Int_t i = 0; i<128; i++){
- fPosLUT[i] = 0;
+
+ AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
+ if (!trdLoader) {
+ AliError("Could not get TRDLoader");
+ return kFALSE;
}
-
- fMCMT = new UInt_t[fMaxTracklets];
- for(Int_t i = 0; i < fMaxTracklets; i++) {
- fMCMT[i] = 0;
+
+ 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()) {
+ digits->Expand();
+
+ if (fNTimeBin != digits->GetNtime()) {
+ AliWarning(Form("Changing no. of timebins from %i to %i", fNTimeBin, digits->GetNtime()));
+ SetNTimebins(digits->GetNtime());
+ }
+
+ SetData(digits);
}
+ else
+ retval = kFALSE;
+ delete digMgr;
- fInitialized = 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.
+ // 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
+ // 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
+ // 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("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
+ 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);
+
+ if (inputTail == 0)
+ valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
+ else if (inputTail == 1)
+ valuet = FilterTailNextSample(1, valuep);
+ else
+ valuet = FilterTailNextSample(1, valueg);
+
+ hfp->SetBinContent(i, valuep >> 2);
+ hfg->SetBinContent(i, valueg >> 2);
+ hft->SetBinContent(i, valuet >> 2);
+ }
+
+ TCanvas *c = new TCanvas;
+ c->Divide(2,2);
+ c->cd(1);
+ h->Draw();
+ c->cd(2);
+ hfp->Draw();
+ c->cd(3);
+ hfg->Draw();
+ c->cd(4);
+ 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::SetPosLUT() {
- Double_t iHi = (Double_t)fCal->GetPRFhi();
- Double_t iLo = (Double_t)fCal->GetPRFlo();
- Int_t nBin = fCal->GetPRFbin();
- Int_t iOff = fLayer * nBin;
- Int_t kNlayer = fGeo->Nlayer();
-
- Float_t *sPRFsmp = new Float_t[nBin*kNlayer];
- Double_t *sPRFlayer = new Double_t[nBin];
-
-
- for(Int_t i = 0; i<nBin*kNlayer; i++){
-
- //printf("%f\n",fCal->GetSampledPRF()[i]);
- sPRFsmp[i] = fCal->GetSampledPRF()[i];
-
- }
+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 following characters:
+ // R - prints raw ADC data
+ // F - prints filtered data
+ // H - prints detected hits
+ // T - prints found tracklets
+ // The later stages are only meaningful after the corresponding calculations
+ // have been performed.
+
+ if ( !CheckInitialized() )
+ return;
- Double_t sWidth = (iHi-iLo)/((Double_t) nBin);
- Int_t sPad = (Int_t) (1.0/sWidth);
-
- // get the PRF for actual layer (interpolated to ibin data-points; 61 measured)
- for(Int_t iBin = 0; iBin < nBin; iBin++){
- sPRFlayer[iBin] = (Double_t)sPRFsmp[iOff+iBin];
- }
+ printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
- Int_t bin0 = (Int_t)(-iLo / sWidth - 0.5); // bin-nr. for pad-position 0
-
- Int_t bin1 = (Int_t)((Double_t)(0.5 - iLo) / sWidth - 0.5); // bin-nr. for pad-position 0.5
- bin1 = bin1 + 1;
- bin0 = bin0 + 1; //avoid negative values in aYest (start right of symmetry center)
- while (bin0-sPad<0) {
- bin0 = bin0 + 1;
+ TString opt = option;
+ if (opt.Contains("R") || opt.Contains("F")) {
+ std::cout << *this;
}
- while (bin1+sPad>=nBin) {
- bin1 = bin1 - 1;
+
+ 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].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos);
+ }
}
-
- Double_t* aYest = new Double_t[bin1-bin0+1];
-
- /*TH1F* hist1 = new TH1F("h1","yest(y)",128,0,0.5);
- TH1F* hist2 = new TH1F("h2","y(yest)",128,0,0.5);
- TH1F* hist3 = new TH1F("h3","y(yest)-yest",128,0,0.5);
- TH1F* hist4 = new TH1F("h4","y(yest)-yest,discrete",128,0,0.5);
-
- TCanvas *c1 = new TCanvas("c1","c1",800,1000);
- hist1->Draw();
- TCanvas *c2 = new TCanvas("c2","c2",800,1000);
- hist2->Draw();
- TCanvas *c3 = new TCanvas("c3","c3",800,1000);
- hist3->Draw();
- TCanvas *c4 = new TCanvas("c4","c4",800,1000);
- hist4->Draw();*/
-
- for(Int_t iBin = bin0; iBin <= bin1; iBin++){
- aYest[iBin-bin0] = 0.5*(sPRFlayer[iBin-sPad] - sPRFlayer[iBin+sPad])/(sPRFlayer[iBin]); // estimated position from PRF; between 0 and 1
- //Double_t position = ((Double_t)(iBin)+0.5)*sWidth+iLo;
- // hist1->Fill(position,aYest[iBin-bin0]);
+
+ if (opt.Contains("T")) {
+ printf("Tracklets:\n");
+ for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntriesFast(); iTrkl++) {
+ printf("tracklet %i: 0x%08x\n", iTrkl, ((AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl])->GetTrackletWord());
+ }
}
-
+}
+
+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 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
+ // T - plot tracklets
+
+ if( !CheckInitialized() )
+ return;
+ TString opt = option;
- Double_t aY[128]; // reversed function
+ TH2F *hist = new TH2F("mcmdata", Form("Data of MCM %i on ROB %i in detector %i", \
+ fMcmPos, fRobPos, fDetector), \
+ fgkNADC, -0.5, fgkNADC-.5, fNTimeBin, -.5, fNTimeBin-.5);
+ hist->GetXaxis()->SetTitle("ADC Channel");
+ hist->GetYaxis()->SetTitle("Timebin");
+ hist->SetStats(kFALSE);
- AliTRDtrapAlu a;
- a.Init(1,8,0,31);
-
- for(Int_t j = 0; j<128; j++) { // loop over all Yest; LUT has 128 entries;
- Double_t yest = ((Double_t)j)/256;
-
- Int_t iBin = 0;
- while (yest>aYest[iBin] && iBin<(bin1-bin0)) {
- iBin = iBin+1;
- }
- if((iBin == bin1 - bin0)&&(yest>aYest[iBin])) {
- aY[j] = 0.5; // yest too big
- //hist2->Fill(yest,aY[j]);
-
+ if (opt.Contains("R")) {
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
+ hist->SetBinContent(iAdc+1, iTimeBin+1, fADCR[iAdc][iTimeBin] >> fgkAddDigits);
+ }
}
- else {
- Int_t bin_d = iBin + bin0 - 1;
- Int_t bin_u = iBin + bin0;
- Double_t y_d = ((Double_t)bin_d + 0.5)*sWidth + iLo; // lower y
- Double_t y_u = ((Double_t)bin_u + 0.5)*sWidth + iLo; // upper y
- Double_t yest_d = aYest[iBin-1]; // lower estimated y
- Double_t yest_u = aYest[iBin]; // upper estimated y
-
- aY[j] = ((yest-yest_d)/(yest_u-yest_d))*(y_u-y_d) + y_d;
- //hist2->Fill(yest,aY[j]);
-
+ }
+ else {
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
+ hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
+ }
}
- aY[j] = aY[j] - yest;
- //hist3->Fill(yest,aY[j]);
- // formatting
- a.AssignDouble(aY[j]);
- //a.WriteWord();
- fPosLUT[j] = a.GetValue(); // 1+8Bit value;128 entries;LUT is steered by abs(Q(i+1)-Q(i-1))/Q(i)=COG and gives the correction to COG/2
- //hist4->Fill(yest,fPosLUT[j]);
-
}
-
-
-
- delete [] sPRFsmp;
- delete [] sPRFlayer;
- delete [] aYest;
-
-}
+ hist->Draw("colz");
+ if (opt.Contains("H")) {
+ TGraph *grHits = new TGraph();
+ for (Int_t iHit = 0; iHit < fNHits; iHit++) {
+ grHits->SetPoint(iHit,
+ fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256.,
+ fHits[iHit].fTimebin);
+ }
+ grHits->Draw("*");
+ }
-//_____________________________________________________________________________
-Int_t* AliTRDmcmSim::GetPosLUT(){
- return fPosLUT;
-}
+ if (opt.Contains("T")) {
+ TLine *trklLines = new TLine[4];
+ for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) {
+ AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
+ 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);
+ trklLines[iTrkl].Draw();
+ }
+ }
+}
-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
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]);
+ 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
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;
+ fADCR[adc][it] = data << fgkAddDigits;
+ fADCF[adc][it] = data << fgkAddDigits;
}
-//_____________________________________________________________________________
-void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
+void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager * const digitsManager)
{
- //
- // Store ADC data into array of raw data
- //
+ // Set the ADC data from an AliTRDarrayADC
- if( !CheckInitialized() ) return;
-
- if( iadc < 0 || iadc >= fNADC ) {
- //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
+ if( !CheckInitialized() )
return;
- }
- for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
- fADCR[iadc][it] = fSimParam->GetADCbaseline();
+ 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", fDetector, iDict));
+ fDict[iDict] = 0x0;
+ }
+ }
}
-}
-//_____________________________________________________________________________
-Int_t AliTRDmcmSim::GetCol( Int_t iadc )
-{
- //
- // Return column id of the pad for the given ADC channel
- //
+ if (fNTimeBin != adcArray->GetNtime())
+ SetNTimebins(adcArray->GetNtime());
- if( !CheckInitialized() ) return -1;
+ Int_t offset = (fMcmPos % 4 + 1) * 21 + (fRobPos % 2) * 84 - 1;
- return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
+ 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);
+ }
+ }
+ }
}
-//_____________________________________________________________________________
-Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
+void AliTRDmcmSim::SetDataByPad(AliTRDarrayADC* const adcArray, AliTRDdigitsManager * const digitsManager)
{
- //
- // Produce raw data stream from this MCM and put in buf
- // Returns number of words filled, or negative value
- // with -1 * number of overflowed words
- //
+ // Set the ADC data from an AliTRDarrayADC
+ // (by pad, to be used during initial reading in simulation)
- UInt_t x;
- UInt_t iEv = 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;
+ if( !CheckInitialized() )
+ return;
- if( !CheckInitialized() ) return 0;
+ 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( fFeeParam->GetRAWstoreRaw() ) {
- adc = fADCR;
- } else {
- adc = fADCF;
- }
+ if (fDict[iDict] == newDict)
+ continue;
- // Produce MCM header
- x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
- if (nw < maxSize) {
- buf[nw++] = x;
- }
- else {
- of++;
- }
+ fDict[iDict] = newDict;
+ fDict[iDict]->Expand();
+ }
+ else {
+ fDict[iDict] = newDict;
+ if (fDict[iDict])
+ fDict[iDict]->Expand();
+ }
- // Produce ADC mask
- 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);
+ // 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", fDetector, iDict));
+ fDict[iDict] = 0x0;
}
}
- if (nw < maxSize) {
- buf[nw++] = x;
- }
- else {
- of++;
- }
}
- // Produce ADC data. 3 timebins are packed into one 32 bits word
- // In this version, different ADC channel will NOT share the same word
+ if (fNTimeBin != adcArray->GetNtime())
+ SetNTimebins(adcArray->GetNtime());
- UInt_t aa=0, a1=0, a2=0, a3=0;
+ Int_t offset = (fMcmPos % 4 + 1) * 18 + (fRobPos % 2) * 72 + 1;
- for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
- if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed
- aa = !(iAdc & 1) + 2;
- for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
- a1 = ((iT ) < fNTimeBin ) ? adc[iAdc][iT ] : 0;
- a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
- a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
- x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
- if (nw < maxSize) {
- buf[nw++] = x;
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ 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 {
- of++;
+ fZSMap[iAdc] = 0;
+ fADCR[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits);
+ fADCF[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits);
}
}
}
+}
- if( of != 0 ) return -of; else return nw;
+void AliTRDmcmSim::SetDataPedestal( Int_t adc )
+{
+ //
+ // Store ADC data into array of raw data
+ //
+
+ if( !CheckInitialized() )
+ return;
+
+ if( adc < 0 || adc >= fgkNADC ) {
+ return;
+ }
+
+ for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
+ fADCR[adc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits);
+ fADCF[adc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP, fDetector, fRobPos, fMcmPos) + (fgAddBaseline << fgkAddDigits);
+ }
+}
+
+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;
+
+ Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adc);
+ if (col < 0 || col >= fFeeParam->GetNcol())
+ return -1;
+ else
+ return col;
}
-//_____________________________________________________________________________
-Int_t AliTRDmcmSim::ProduceRawStreamV2( 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 iEv = 0;
+ 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 : xrrr mmmm eeee eeee eeee eeee eeee 1100
- // x : 0 before , 1 since 10.2007
- // r : Readout board position (Alice numbering)
- // m : MCM posi
- // e : Event counter from 1
- //x = (1<<31) | ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
- 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
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
- aa = !(iAdc & 1) + 2; // 3 for the even ADC channel , 2 for the odd ADC channel
+ 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 ] : 0;
- a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
- a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
+ 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) {
- buf[nw++] = x;
- //printf("%08X ",x);
+ if (nw < bufSize) {
+ buf[nw++] = x;
}
else {
- of++;
+ of++;
}
}
}
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;
}
-
-//_____________________________________________________________________________
void AliTRDmcmSim::Filter()
{
//
- // Apply digital 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
+ // sequentially for parameter tuning.
//
- if( !CheckInitialized() ) return;
+ if( !CheckInitialized() )
+ return;
- // Initialize filtered data array with raw data
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
- fADCF[iadc][it] = fADCR[iadc][it];
- }
- }
+ // Apply filters sequentially. Bypass is handled by filters
+ // 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.
+
+ // Non-linearity filter not implemented.
+ FilterPedestal();
+ FilterGain();
+ FilterTail();
+ // Crosstalk filter not implemented.
+}
- // Then apply fileters one by one to filtered data array
- if( fFeeParam->IsPFon() ) FilterPedestal();
- if( fFeeParam->IsGFon() ) FilterGain();
- if( fFeeParam->IsTFon() ) FilterTail();
+void AliTRDmcmSim::FilterPedestalInit(Int_t baseline)
+{
+ // Initializes the pedestal filter assuming that the input has
+ // been constant for a long time (compared to the time constant).
+
+ UShort_t fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC, fDetector, fRobPos, fMcmPos); // 0..3, 0 - fastest, 3 - slowest
+
+ for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++)
+ fPedAcc[iAdc] = (baseline << 2) * (1 << fgkFPshifts[fptc]);
}
-//_____________________________________________________________________________
-void AliTRDmcmSim::FilterPedestal()
+UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
{
+ // Returns the output of the pedestal filter given the input value.
+ // The output depends on the internal registers and, thus, the
+ // history of the filter.
-
- //
- // Apply pedestal filter
- //
+ 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;
- Int_t ap = fSimParam->GetADCbaseline(); // ADC instrinsic pedestal
- Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
- //Int_t tc = fFeeParam->GetPFtimeConstant(); // this makes no sense yet
+ inpAdd = value + fpnp;
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
- fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
- }
+ 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
}
-}
-//_____________________________________________________________________________
-void AliTRDmcmSim::FilterGain()
-{
- //
- // Apply gain filter (not implemented)
- // Later it will be implemented because gain digital filiter will
- // increase noise level.
- //
+ if (fpby == 0)
+ return value;
+ if (inpAdd <= accumulatorShifted)
+ return 0;
+ else
+ {
+ inpAdd = inpAdd - accumulatorShifted;
+ if (inpAdd > 0xFFF)
+ return 0xFFF;
+ else
+ return inpAdd;
+ }
}
-//_____________________________________________________________________________
-void AliTRDmcmSim::FilterTail()
+void AliTRDmcmSim::FilterPedestal()
{
//
- // Apply exponential tail filter (Bogdan's version)
+ // 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
+ // the input has been stable for a sufficiently long time.
- Double_t *dtarg = new Double_t[fNTimeBin];
- Int_t *itarg = new Int_t[fNTimeBin];
- Int_t nexp = fFeeParam->GetTFnExp();
- Int_t tftype = fFeeParam->GetTFtype();
-
- switch( tftype ) {
-
- case 0: // Exponential Filter Analog Bogdan
- for (Int_t iCol = 0; iCol < fNADC; iCol++) {
- FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
- for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
- fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
- }
- }
- break;
-
- case 1: // Exponential filter digital Bogdan
- for (Int_t iCol = 0; iCol < fNADC; iCol++) {
- FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
- for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
- fADCF[iCol][iTime] = itarg[iTime];
- }
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
+ fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
}
- break;
-
- case 2: // Exponential filter Marian special
- for (Int_t iCol = 0; iCol < fNADC; iCol++) {
- FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
- for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
- fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
- }
- }
- break;
-
- //new
- case 3: // Exponential filter using AliTRDtrapAlu class
- for (Int_t iCol = 0; iCol < fNADC; iCol++) {
- FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp);
- for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
- fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data
- fADCT[iCol][iTime] = itarg[iTime]; // 12bits; to be used for tracklet; tracklet will have own container;
- }
- }
- break;
-
-
- default:
- AliError(Form("Invalid filter type %d ! \n", tftype ));
- break;
}
+}
- delete [] dtarg;
- delete [] itarg;
+void AliTRDmcmSim::FilterGainInit()
+{
+ // Initializes the gain filter. In this case, only threshold
+ // counters are reset.
+ 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;
+ }
+}
+
+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
+ // history of the filter.
+
+ 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 corr; // corrected value
+
+ value &= 0xFFF;
+ corr = (value * fgf) >> 11;
+ corr = corr > 0xfff ? 0xfff : corr;
+ corr = AddUintClipping(corr, fga, 12);
+
+ // Update threshold counters
+ // not really useful as they are cleared with every new event
+ if (!((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF)))
+ // stop when full
+ {
+ if (corr >= fgtb)
+ fGainCounterB[adc]++;
+ else if (corr >= fgta)
+ fGainCounterA[adc]++;
+ }
+
+ 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 < fgkNADC; iAdc++) {
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
+ }
+ }
+}
+
+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
+ // sufficiently long time.
+
+ // exponents and weight calculated from configuration
+ 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);
+ 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;
+
+ Float_t kt, ql, qs;
+ UShort_t aout;
+
+ 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 < 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));
+ }
+}
+
+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
+ // history of the filter.
+
+ // exponents and weight calculated from configuration
+ 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
+
+ // intermediate signals
+ UInt_t aDiff;
+ UInt_t alInpv;
+ UShort_t aQ;
+ 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
+ return aDiff;
+}
+
+void AliTRDmcmSim::FilterTail()
+{
+ // Apply tail cancellation filter to all data.
+
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
+ fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
+ }
+ }
}
-//_____________________________________________________________________________
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
//
- Int_t eBIS = fFeeParam->GetEBsglIndThr(); // TRAP default = 0x4 (Tis=4)
- Int_t eBIT = fFeeParam->GetEBsumIndThr(); // TRAP default = 0x28 (Tit=40)
- Int_t eBIL = fFeeParam->GetEBindLUT(); // TRAP default = 0xf0
- // (lookup table accept (I2,I1,I0)=(111)
- // or (110) or (101) or (100))
- Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // 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);
- if( !CheckInitialized() ) return;
+ Int_t **adc = fADCF;
- for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
+ for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++)
+ fZSMap[iAdc] = -1;
- // Get ADC data currently in filter buffer
- Int_t ap = fADCF[iadc-1][it] - ep; // previous
- Int_t ac = fADCF[iadc ][it] - ep; // current
- Int_t an = fADCF[iadc+1][it] - ep; // next
+ for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
+ 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)
- // 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
+ // ----- first channel -----
+ iAdc = 0;
- 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)
+ ap = 0; // previous
+ ac = adc[iAdc ][it]; // current
+ an = adc[iAdc+1][it]; // next
- fZSM[iadc][it] &= d;
- if( eBIN == 0 ) { // turn on neighboring ADCs
- fZSM[iadc-1][it] &= d;
- fZSM[iadc+1][it] &= d;
- }
+ mask = ( ac >= ap && ac >= an ) ? 0 : 0x1; // peak center detection
+ mask += ( ap + ac + an > eBIT ) ? 0 : 0x2; // cluster
+ mask += ( ac > eBIS ) ? 0 : 0x4; // absolute large peak
+
+ supp = (eBIL >> mask) & 1;
+ fZSMap[iAdc] &= ~((1-supp) << it);
+ if( eBIN == 0 ) { // neighbour sensitivity
+ fZSMap[iAdc+1] &= ~((1-supp) << it);
}
- }
- // do 1 dim projection
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
- fZSM1Dim[iadc] &= fZSM[iadc][it];
+ // ----- last channel -----
+ iAdc = fgkNADC - 1;
+
+ ap = adc[iAdc-1][it]; // previous
+ ac = adc[iAdc ][it]; // current
+ an = 0; // next
+
+ mask = ( ac >= ap && ac >= an ) ? 0 : 0x1; // peak center detection
+ mask += ( ap + ac + an > eBIT ) ? 0 : 0x2; // cluster
+ mask += ( ac > eBIS ) ? 0 : 0x4; // absolute large peak
+
+ supp = (eBIL >> mask) & 1;
+
+ fZSMap[iAdc] &= ~((1-supp) << it);
+ if( eBIN == 0 ) { // neighbour sensitivity
+ fZSMap[iAdc-1] &= ~((1-supp) << it);
}
+
+ // ----- 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[])
+{
+ // 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, 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].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].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::DumpData( char *f, char *target )
+void AliTRDmcmSim::CalcFitreg()
{
- //
- // 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
- //
+ // Preprocessing.
+ // Detect the hits and fill the fit registers.
+ // Requires 12-bit data from fADCF which means Filter()
+ // has to be called before even if all filters are bypassed.
+
+ //??? to be clarified:
+ UInt_t adcMask = 0xffffffff;
+
+ UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
+ Short_t ypos, fromLeft, fromRight, found;
+ 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, 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, fDetector, fRobPos, fMcmPos);
+
+ // reset the fit registers
+ fNHits = 0;
+ for (adcch = 0; adcch < fgkNADC-2; adcch++) // due to border channels
+ {
+ 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 qTotal array
+ // in case of not hit store 0 there.
+ for (adcch = 0; adcch < fgkNADC-2; adcch++) {
+ if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
+ {
+ adcLeft = fADCF[adcch ][timebin];
+ adcCentral = fADCF[adcch+1][timebin];
+ adcRight = fADCF[adcch+2][timebin];
+ if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY, 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!!!
+ qtotTemp = adcLeft + adcCentral + adcRight;
+ if ( (hitQual) &&
+ (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT, fDetector, fRobPos, fMcmPos)) &&
+ (adcLeft <= adcCentral) &&
+ (adcCentral > adcRight) )
+ qTotal[adcch] = qtotTemp;
+ else
+ qTotal[adcch] = 0;
+ }
+ else
+ qTotal[adcch] = 0; //jkl
+ if (qTotal[adcch] != 0)
+ AliDebug(10,Form("ch %2d qTotal %5d",adcch, qTotal[adcch]));
+ }
+
+ fromLeft = -1;
+ adcch = 0;
+ found = 0;
+ marked[4] = 19; // invalid channel
+ marked[5] = 19; // invalid channel
+ qTotal[19] = 0;
+ while ((adcch < 16) && (found < 3))
+ {
+ if (qTotal[adcch] > 0)
+ {
+ fromLeft = adcch;
+ marked[2*found+1]=adcch;
+ found++;
+ }
+ adcch++;
+ }
+
+ fromRight = -1;
+ adcch = 18;
+ found = 0;
+ while ((adcch > 2) && (found < 3))
+ {
+ if (qTotal[adcch] > 0)
+ {
+ marked[2*found]=adcch;
+ found++;
+ fromRight = adcch;
+ }
+ adcch--;
+ }
- UInt_t tempbuf[1024];
+ AliDebug(10,Form("Fromleft=%d, Fromright=%d",fromLeft, fromRight));
+ // here mask the hit candidates in the middle, if any
+ if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
+ for (adcch = fromLeft+1; adcch < fromRight; adcch++)
+ qTotal[adcch] = 0;
- if( !CheckInitialized() ) return;
+ found = 0;
+ for (adcch = 0; adcch < 19; adcch++)
+ if (qTotal[adcch] > 0) found++;
+ // NOT READY
- 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",
- fChaId, fSector, fStack, fLayer, 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");
+ if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
+ {
+ if (marked[4] == marked[5]) marked[5] = 19;
+ for (found=0; found<6; found++)
+ {
+ qMarked[found] = qTotal[marked[found]] >> 4;
+ AliDebug(10,Form("ch_%d qTotal %d qTotals %d",marked[found],qTotal[marked[found]],qMarked[found]));
}
- 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");
+
+ 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],
+ &worse1, &worse2);
+ // Now mask the two channels with the smallest charge
+ if (worse1 < 19)
+ {
+ qTotal[worse1] = 0;
+ AliDebug(10,Form("Kill ch %d\n",worse1));
}
- 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]);
+ if (worse2 < 19)
+ {
+ qTotal[worse2] = 0;
+ AliDebug(10,Form("Kill ch %d\n",worse2));
}
}
- }
-}
-//_____________________________________________________________________________
-void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
- , Int_t n, Int_t nexp)
-{
- //
- // Exponential filter "analog"
- // source will not be changed
- //
+ for (adcch = 0; adcch < 19; adcch++) {
+ if (qTotal[adcch] > 0) // the channel is marked for processing
+ {
+ 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 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
+ // 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 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]];
+ }
+ }
- Int_t i = 0;
- Int_t k = 0;
- Double_t reminder[2];
- Double_t correction;
- Double_t result;
- Double_t rates[2];
- Double_t coefficients[2];
-
- // Initialize (coefficient = alpha, rates = lambda)
- // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
-
- Double_t r1 = (Double_t)fFeeParam->GetTFr1();
- Double_t r2 = (Double_t)fFeeParam->GetTFr2();
- Double_t c1 = (Double_t)fFeeParam->GetTFc1();
- Double_t c2 = (Double_t)fFeeParam->GetTFc2();
-
- coefficients[0] = c1;
- coefficients[1] = c2;
-
- Double_t dt = 0.1;
- rates[0] = TMath::Exp(-dt/(r1));
- rates[1] = TMath::Exp(-dt/(r2));
-
- // Attention: computation order is important
- correction = 0.0;
- for (k = 0; k < nexp; k++) {
- reminder[k] = 0.0;
- }
-
- for (i = 0; i < n; i++) {
-
- result = ((Double_t)source[i] - correction); // no rescaling
- target[i] = result;
-
- for (k = 0; k < nexp; k++) {
- reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
+ // add the hit to the fitregister
+ AddHitToFitreg(adcch, timebin, qTotal[adcch] >> fgkAddDigits, ypos, mcLabel);
+ }
}
-
- correction = 0.0;
- for (k = 0; k < nexp; k++) {
- correction += reminder[k];
+ }
+
+ 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::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
- , Int_t nexp)
+void AliTRDmcmSim::TrackletSelection()
{
- //
- // Exponential filter "digital"
- // source will not be changed
- //
-
- Int_t i = 0;
- Int_t fAlphaL = 0;
- Int_t fAlphaS = 0;
- Int_t fTailPed = 0;
- Int_t iAlphaL = 0;
- Int_t iAlphaS = 0;
-
- // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
- // initialize (coefficient = alpha, rates = lambda)
-
- Double_t dt = 0.1;
- Double_t r1 = (Double_t)fFeeParam->GetTFr1();
- Double_t r2 = (Double_t)fFeeParam->GetTFr2();
- Double_t c1 = (Double_t)fFeeParam->GetTFc1();
- Double_t c2 = (Double_t)fFeeParam->GetTFc2();
-
- Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
- Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
- Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
- Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
-
- if (nexp == 1) {
- fAlphaL = (Int_t) (c1 * 2048.0);
- iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
- }
- if (nexp == 2) {
- fAlphaL = (Int_t) (c1 * 2048.0);
- fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
- iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
- iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
- }
-
- Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
- Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
- Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
- Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
-
- Int_t h1;
- Int_t h2;
- Int_t rem1;
- Int_t rem2;
- Int_t correction;
- Int_t result;
- Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
-
- Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
- rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
- rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
-
- // further initialization
- if ((rem1 + rem2) > 0x0FFF) {
- correction = 0x0FFF;
- }
- else {
- correction = (rem1 + rem2) & 0x0FFF;
- }
-
- fTailPed = iFactor - correction;
-
- for (i = 0; i < n; i++) {
-
- result = (source[i] - correction);
- if (result < 0) { // Too much undershoot
- result = 0;
+ // Select up to 4 tracklet candidates from the fit registers
+ // and assign them to the CPUs.
+
+ 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 (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)))
+ {
+ 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<ntracks;i++)
+ AliDebug(10,Form("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]));
+
+ if (ntracks > 4)
+ {
+ // primitive sorting according to the number of hits
+ for (j = 0; j < (ntracks-1); j++)
+ {
+ for (i = j+1; i < ntracks; i++)
+ {
+ if ( (trackletCand[j][1] < trackletCand[i][1]) ||
+ ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
+ {
+ // swap j & i
+ 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
- target[i] = result;
-
- h1 = (rem1 + ((iAlphaL * result) >> 11));
- if (h1 > 0x0FFF) {
- h1 = 0x0FFF;
- }
- else {
- h1 &= 0x0FFF;
+ // 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 (trackletCand[j][0] < trackletCand[i][0])
+ {
+ // swap j & i
+ 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] = 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)
+ 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]));
+}
- h2 = (rem2 + ((iAlphaS * result) >> 11));
- if (h2 > 0x0FFF) {
- h2 = 0x0FFF;
- }
- else {
- h2 &= 0x0FFF;
+void AliTRDmcmSim::FitTracklet()
+{
+ // 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 rndAdd = 0;
+ 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;
+
+ Int_t ndriftDp = 5; // decimal places for drift time
+ Long64_t shift = ((Long64_t) 1 << 32);
+
+ // 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, 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
+ 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,
+// 0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
+// 0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
+
+ for (Int_t cpu = 0; cpu < 4; cpu++) {
+ if (fFitPtr[cpu] == 31)
+ {
+ fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
}
-
- rem1 = (iLambdaL * h1 ) >> 11;
- rem2 = (iLambdaS * h2 ) >> 11;
-
- if ((rem1 + rem2) > 0x0FFF) {
- correction = 0x0FFF;
- }
- else {
- correction = (rem1 + rem2) & 0x0FFF;
+ else
+ {
+ fit0 = &fFitReg[fFitPtr[cpu] ];
+ fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
+
+ mult = 1;
+ mult = mult << (32 + decPlaces);
+ mult = -mult;
+
+ // Merging
+ nHits = fit0->fNhits + fit1->fNhits; // number of hits
+ sumX = fit0->fSumX + fit1->fSumX;
+ sumX2 = fit0->fSumX2 + fit1->fSumX2;
+ denom = nHits * sumX2 - sumX * sumX;
+
+ mult = mult / denom; // exactly like in the TRAP program
+ q0 = fit0->fQ0 + fit1->fQ0;
+ q1 = fit0->fQ1 + fit1->fQ1;
+ sumY = fit0->fSumY + fit1->fSumY + 256*fit1->fNhits;
+ sumXY = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
+ sumY2 = fit0->fSumY2 + fit1->fSumY2 + 512*fit1->fSumY + 256*256*fit1->fNhits;
+
+ slope = nHits*sumXY - sumX * sumY;
+ 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;
+ 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));
+
+ 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
+ {
+ 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)
+ AliWarning("Overflow in offset");
+ offset = offset & 0x1FFF; // 13 bit
+
+ 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] = (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)));
+ }
}
-
}
-
}
-//_____________________________________________________________________________
-void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
- , Int_t n)
+void AliTRDmcmSim::Tracklet()
{
- //
- // Exponential filter (M. Ivanov)
- // source will not be changed
- //
-
- Int_t i = 0;
- Double_t sig1[100];
- Double_t sig2[100];
- Double_t sig3[100];
+ // Run the tracklet calculation by calling sequentially:
+ // CalcFitreg(); TrackletSelection(); FitTracklet()
+ // and store the tracklets
- for (i = 0; i < n; i++) {
- sig1[i] = (Double_t)source[i];
+ if (!fInitialized) {
+ AliError("Called uninitialized! Nothing done!");
+ return;
}
- Float_t dt = 0.1;
- Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
- Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
-
- FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
- FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
-
- for (i = 0; i < n; i++) {
- target[i] = sig3[i];
- }
+ fTrackletArray->Delete();
+ CalcFitreg();
+ if (fNHits == 0)
+ return;
+ TrackletSelection();
+ FitTracklet();
}
-//______________________________________________________________________________
-void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
- , Double_t lambda, Int_t n)
+Bool_t AliTRDmcmSim::StoreTracklets()
{
- //
- // Special filter (M. Ivanov)
- //
-
- Int_t i = 0;
- Double_t l = TMath::Exp(-lambda*0.5);
- Double_t in[1000];
- Double_t out[1000];
+ // store the found tracklets via the loader
- // Initialize in[] and out[] goes 0 ... 2*n+19
- for (i = 0; i < n*2+20; i++) {
- in[i] = out[i] = 0;
- }
+ if (fTrackletArray->GetEntriesFast() == 0)
+ return kTRUE;
- // in[] goes 0, 1
- in[0] = ampin[0];
- in[1] = (ampin[0] + ampin[1]) * 0.5;
-
- // Add charge to the end
- for (i = 0; i < 22; i++) {
- in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
+ AliRunLoader *rl = AliRunLoader::Instance();
+ AliDataLoader *dl = 0x0;
+ if (rl)
+ dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
+ if (!dl) {
+ AliError("Could not get the tracklets data loader!");
+ return kFALSE;
}
- // Use arithmetic mean
- for (i = 1; i < n-1; i++) {
- in[2*i] = ampin[i]; // in[] goes 2, 3, ... , 2*n-4, 2*n-3
- in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
+ TTree *trackletTree = dl->Tree();
+ if (!trackletTree) {
+ dl->MakeTree();
+ trackletTree = dl->Tree();
}
- Double_t temp;
- out[2*n] = in[2*n];
- temp = 0;
- for (i = 2*n; i >= 0; i--) {
- out[i] = in[i] + temp;
- temp = l*(temp+in[i]);
- }
+ AliTRDtrackletMCM *trkl = 0x0;
+ TBranch *trkbranch = trackletTree->GetBranch(fTrklBranchName.Data());
+ if (!trkbranch)
+ trkbranch = trackletTree->Branch(fTrklBranchName.Data(), "AliTRDtrackletMCM", &trkl, 32000);
- for (i = 0; i < n; i++){
- //ampout[i] = out[2*i+1]; // org
- ampout[i] = out[2*i];
+ for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
+ trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
+ trkbranch->SetAddress(&trkl);
+ trkbranch->Fill();
}
+ return kTRUE;
}
-//______________________________________________________________________________
-void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
- , Double_t norm, Double_t lambda
- , Int_t n)
+void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
{
- //
- // Special filter (M. Ivanov)
- //
+ // write back the processed data configured by EBSF
+ // EBSF = 1: unfiltered data; EBSF = 0: filtered data
+ // zero-suppressed valued are written as -1 to digits
- Int_t i = 0;
+ if( !CheckInitialized() )
+ return;
- Double_t l = TMath::Exp(-lambda*0.5);
- Double_t k = l*(1.0 - norm*lambda*0.5);
- Double_t in[1000];
- Double_t out[1000];
+ Int_t offset = (fMcmPos % 4 + 1) * 21 + (fRobPos % 2) * 84 - 1;
- // Initialize in[] and out[] goes 0 ... 2*n+19
- for (i = 0; i < n*2+20; i++) {
- in[i] = out[i] = 0;
+ if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF, fDetector, fRobPos, fMcmPos) != 0) // store unfiltered data
+ {
+ for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
+ if (~fZSMap[iAdc] == 0) {
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; 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 = 0; iAdc < fgkNADC; iAdc++) {
+ if (~fZSMap[iAdc] != 0) {
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, (fADCF[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline);
+ }
+ }
+ else {
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, -1);
+ }
+ }
+ }
}
+}
- // in[] goes 0, 1
- in[0] = ampin[0];
- in[1] = (ampin[0]+ampin[1])*0.5;
- // Add charge to the end
- for (i =-2; i < 22; i++) {
- // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
- in[2*(n-1)+i] = ampin[n-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;
+}
- for (i = 1; i < n-2; i++) {
- // in[] goes 2, 3, ... , 2*n-6, 2*n-5
- in[2*i] = ampin[i];
- in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
- //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
- }
- Double_t temp;
- out[0] = in[0];
- temp = in[0];
- for (i = 1; i <= 2*n; i++) {
- out[i] = in[i] + (k-l)*temp;
- temp = in[i] + k *temp;
- }
- for (i = 0; i < n; i++) {
- //ampout[i] = out[2*i+1]; // org
- //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
- ampout[i] = TMath::Max(out[2*i],0.0);
+// help functions, to be cleaned up
+
+UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
+{
+ //
+ // This function adds a and b (unsigned) and clips to
+ // the specified number of bits.
+ //
+
+ UInt_t sum = a + b;
+ if (nbits < 32)
+ {
+ UInt_t maxv = (1 << nbits) - 1;;
+ if (sum > maxv)
+ sum = maxv;
+ }
+ else
+ {
+ if ((sum < a) || (sum < b))
+ sum = 0xFFFFFFFF;
}
+ return sum;
}
+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
-//_____________________________________________________________________________________
-//the following filter uses AliTRDtrapAlu-class
-
-void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
- //static Int_t count = 0;
-
- Double_t dt = 0.1;
- Double_t r1 = (Double_t)fFeeParam->GetTFr1();
- Double_t r2 = (Double_t)fFeeParam->GetTFr2();
- Double_t c1 = (Double_t)fFeeParam->GetTFc1();
- Double_t c2 = (Double_t)fFeeParam->GetTFc2();
-
- nexp = 1;
-
- //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
- //parameters need to be adjusted
- AliTRDtrapAlu lambdaL;
- AliTRDtrapAlu lambdaS;
- AliTRDtrapAlu alphaL;
- AliTRDtrapAlu alphaS;
-
- AliTRDtrapAlu correction;
- AliTRDtrapAlu result;
- AliTRDtrapAlu bufL;
- AliTRDtrapAlu bufS;
-
- AliTRDtrapAlu bSource;
-
- lambdaL.Init(1,11);
- lambdaS.Init(1,11);
- alphaL.Init(1,11);
- alphaS.Init(1,11);
-
- //count=count+1;
-
- lambdaL.AssignDouble(TMath::Exp(-dt/r1));
- lambdaS.AssignDouble(TMath::Exp(-dt/r2));
- alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
- alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
-
- //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
- correction.Init(10,2);
- result.Init(10,2);
- bufL.Init(10,2);
- bufS.Init(10,2);
- bSource.Init(10,2);
-
- for(Int_t i = 0; i < n; i++) {
- bSource.AssignInt(source[i]);
- result = bSource - correction; // subtraction can produce an underflow
- if(result.GetSign() == kTRUE) {
- result.AssignInt(0);
+ if (val1i > val2i)
+ {
+ *idx1o = idx1i;
+ *idx2o = idx2i;
+ *val1o = val1i;
+ *val2o = val2i;
+ }
+ else
+ {
+ *idx1o = idx2i;
+ *idx2o = idx1i;
+ *val1o = val2i;
+ *val2o = val1i;
}
-
- //target[i] = result.GetValuePre(); // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
-
- target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2
-
- //printf("target-Wert zur Zeit %d : %d",i,target[i]);
- //printf("\n");
-
- bufL = bufL + (result * alphaL);
- bufL = bufL * lambdaL;
-
- bufS = bufS + (result * alphaS);
- bufS = bufS * lambdaS; // eventually this should look like:
- // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
-
- correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
- }
-
-
}
+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)
+{
+ // 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;
+ switch(sel)
+ {
+ case 6 : // 1 > 2 > 3 => 1 2 3
+ case 0 : // 1 = 2 = 3 => 1 2 3 : in this case doesn't matter, but so is in hardware!
+ *idx1o = idx1i;
+ *idx2o = idx2i;
+ *idx3o = idx3i;
+ *val1o = val1i;
+ *val2o = val2i;
+ *val3o = val3i;
+ break;
+
+ case 4 : // 1 > 2, 2 <= 3, 3 <= 1 => 1 3 2
+ *idx1o = idx1i;
+ *idx2o = idx3i;
+ *idx3o = idx2i;
+ *val1o = val1i;
+ *val2o = val3i;
+ *val3o = val2i;
+ break;
+
+ case 2 : // 1 <= 2, 2 > 3, 3 <= 1 => 2 1 3
+ *idx1o = idx2i;
+ *idx2o = idx1i;
+ *idx3o = idx3i;
+ *val1o = val2i;
+ *val2o = val1i;
+ *val3o = val3i;
+ break;
+
+ case 3 : // 1 <= 2, 2 > 3, 3 > 1 => 2 3 1
+ *idx1o = idx2i;
+ *idx2o = idx3i;
+ *idx3o = idx1i;
+ *val1o = val2i;
+ *val2o = val3i;
+ *val3o = val1i;
+ break;
+
+ case 1 : // 1 <= 2, 2 <= 3, 3 > 1 => 3 2 1
+ *idx1o = idx3i;
+ *idx2o = idx2i;
+ *idx3o = idx1i;
+ *val1o = val3i;
+ *val2o = val2i;
+ *val3o = val1i;
+ break;
+
+ case 5 : // 1 > 2, 2 <= 3, 3 > 1 => 3 1 2
+ *idx1o = idx3i;
+ *idx2o = idx1i;
+ *idx3o = idx2i;
+ *val1o = val3i;
+ *val2o = val1i;
+ *val3o = val2i;
+ break;
+
+ default: // the rest should NEVER happen!
+ AliError("ERROR in Sort3!!!\n");
+ break;
+ }
+}
+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
+ 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,
+ val1o, &val21s, &val23as);
+ Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
+ idx2o, &idx22s, &idx23bs,
+ val2o, &val22s, &val23bs);
-// in order to use the Tracklets, please first
-// -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
-// -- set AliTRDfeeParam::fgkTFtype to 3, in order to use the new tail cancellation filter
-// currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
-// -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with information about the Tracklet position (MCM, channel number)
+ Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
-// The code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format
+ Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
+ idx3o, idx4o, &dummy,
+ val3o, val4o, &dummy);
-void AliTRDmcmSim::Tracklet(){
- // tracklet calculation
- // if you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
-
- if(!CheckInitialized()){ return; }
-
- Bool_t filtered = kTRUE;
-
-
-
- AliTRDtrapAlu data;
- data.Init(10,2);
- if(fADCT[0][0]==-1){ // check if filter was applied
- filtered = kFALSE;
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
- data.AssignInt(fADCR[iadc][iT]);
- fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called, this is done correctly
- }
- }
-
- }
-
- // the online ordering of mcm's is reverse to the TRAP-manual-ordering! reverse fADCT (to be consistent to TRAP), then do all calculations
- // reverse fADCT:
- Int_t** rev0 = new Int_t *[fNADC];
- Int_t** rev1 = new Int_t *[fNADC];
-
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- rev0[iadc] = new Int_t[fNTimeBin];
- rev1[iadc] = new Int_t[fNTimeBin];
- for( Int_t iT = 0; iT < fNTimeBin; iT++) {
- if( iadc <= fNADC-iadc-1 ) {
- rev0[iadc][iT] = fADCT[fNADC-iadc-1][iT];
- rev1[iadc][iT] = fADCT[iadc][iT];
- fADCT[iadc][iT] = rev0[iadc][iT];
- }
- else {
- rev0[iadc][iT] = rev1[fNADC-iadc-1][iT];
- fADCT[iadc][iT] = rev0[iadc][iT];
- }
- }
- }
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
- delete[] rev0[iadc];
- delete[] rev1[iadc];
- }
-
- delete[] rev0;
- delete[] rev1;
-
- rev0 = NULL;
- rev1 = NULL;
-
- // get the filtered pedestal; supports only electronic tail-cancellation filter
- AliTRDtrapAlu filPed;
- Int_t ep = 0;
- Int_t *ieffped = new Int_t[fNTimeBin];
- for(Int_t iT = 0; iT < fNTimeBin; iT++){
- ieffped[iT] = ep;
- }
-
- if( filtered == kTRUE ) {
- if( fFeeParam->IsPFon() ){
- ep = fFeeParam->GetPFeffectPedestal();
- }
- Int_t nexp = fFeeParam->GetTFnExp();
- Int_t *isource = new Int_t[fNTimeBin];
- filPed.Init(10,2);
- filPed.AssignInt(ep);
- Int_t epf = filPed.GetValue();
- for(Int_t iT = 0; iT < fNTimeBin; iT++){
- isource[iT] = ep;
- ieffped[iT] = epf;
- }
-
- if( fFeeParam->IsTFon() ) {
- FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
- }
-
- delete[] isource;
- }
-
- //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
- //naming follows conventions in TRAP-manual
-
-
- Bool_t bVBY = kTRUE; // cluster-verification bypass
-
- Double_t cQTParam = 0; // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
- AliTRDtrapAlu cQTAlu;
- cQTAlu.Init(1,10,0,63);
- cQTAlu.AssignDouble(cQTParam);
- Int_t cQT = cQTAlu.GetValue();
-
- // linear fit
- Int_t tFS = fFeeParam->GetLinearFitStart(); // linear fit start
- Int_t tFE = fFeeParam->GetLinearFitEnd(); // linear fit stop
-
- // charge accumulators
- Int_t tQS0 = fFeeParam->GetQacc0Start(); // start-time for charge-accumulator 0
- Int_t tQE0 = fFeeParam->GetQacc0End(); // stop-time for charge-accumulator 0
- Int_t tQS1 = fFeeParam->GetQacc1Start(); // start-time for charge-accumulator 1
- Int_t tQE1 = fFeeParam->GetQacc1End(); // stop-time for charge-accumulator 1
- // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
-
- Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
- AliTRDtrapAlu cTHAlu;
- cTHAlu.Init(12,2);
- cTHAlu.AssignDouble(cTHParam);
- Int_t cTH = cTHAlu.GetValue(); // cTH used for comparison
-
- struct List_t {
- List_t *next;
- Int_t iadc;
- Int_t value;
- };
-
- List_t selection[7]; // list with 7 elements
- List_t *list = NULL;
- List_t *listLeft = NULL;
-
- Int_t* qsum = new Int_t[fNADC];
-
- // fit sums
- AliTRDtrapAlu qsumAlu;
- qsumAlu.Init(12,2); // charge sum will be 12+2 bits
- AliTRDtrapAlu dCOGAlu;
- dCOGAlu.Init(1,7,0,127); // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
- AliTRDtrapAlu yrawAlu;
- yrawAlu.Init(1,8,-1,255);
- AliTRDtrapAlu yAlu;
- yAlu.Init(1,16,-1,0xFF00); // only first 8 past-comma bits filled;additional 8 bits for accuracy;maximum 1 - 2^-8; sign is given by + or -
- AliTRDtrapAlu xAlu;
- xAlu.Init(5,8); // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
- AliTRDtrapAlu xxAlu;
- xxAlu.Init(10,0);
- AliTRDtrapAlu yyAlu;
- yyAlu.Init(1,16,0,0xFFFF); // maximum is 2^16-1; 16Bit for past-commas
- AliTRDtrapAlu xyAlu;
- xyAlu.Init(6,8);
- AliTRDtrapAlu XAlu;
- XAlu.Init(9,0);
- AliTRDtrapAlu XXAlu;
- XXAlu.Init(14,0);
- AliTRDtrapAlu YAlu;
- YAlu.Init(5,8); // 14 bit, 1 is sign-bit; therefore only 13 bit
- AliTRDtrapAlu YYAlu;
- YYAlu.Init(5,16);
- AliTRDtrapAlu XYAlu;
- XYAlu.Init(8,8); // 17 bit, 1 is sign-bit; therefore only 16 bit
- AliTRDtrapAlu qtruncAlu;
- qtruncAlu.Init(12,0);
- AliTRDtrapAlu QT0Alu;
- QT0Alu.Init(15,0);
- AliTRDtrapAlu QT1Alu;
- QT1Alu.Init(16,0);
-
- AliTRDtrapAlu oneAlu;
- oneAlu.Init(1,8);
-
-
- AliTRDtrapAlu inverseNAlu;
- inverseNAlu.Init(1,8); // simulates the LUT for 1/N
- AliTRDtrapAlu MeanChargeAlu; // mean charge in ADC counts
- MeanChargeAlu.Init(8,0);
- AliTRDtrapAlu TotalChargeAlu;
- TotalChargeAlu.Init(17,8);
- //nr of post comma bits should be the same for inverseN and TotalCharge
-
-
- SetPosLUT(); // initialize the position correction LUT for this MCM;
-
-
- // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
- Int_t *X = new Int_t[fNADC-2];
- Int_t *XX = new Int_t[fNADC-2];
- Int_t *Y = new Int_t[fNADC-2];
- Int_t *YY = new Int_t[fNADC-2];
- Int_t *XY = new Int_t[fNADC-2];
- Int_t *N = new Int_t[fNADC-2];
- Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
- Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
-
- for (Int_t iCol = 0; iCol < fNADC-2; iCol++) {
-
- // initialize fit-sums
- X[iCol] = 0;
- XX[iCol] = 0;
- Y[iCol] = 0;
- YY[iCol] = 0;
- XY[iCol] = 0;
- N[iCol] = 0;
- QT0[iCol] = 0;
- QT1[iCol] = 0;
- }
-
-
- filPed.Init(7,2); // convert filtered pedestal into 7+2Bits
-
- for(Int_t iT = 0; iT < fNTimeBin; iT++){
-
- if(iT<tFS || iT>=tFE) continue; // linear fit yes/no?
-
- // reset
- Int_t portChannel[4] = {-1,-1,-1,-1};
- Int_t clusterCharge[4] = {0,0,0,0};
- Int_t leftCharge[4] = {0,0,0,0};
- Int_t centerCharge[4] = {0,0,0,0};
- Int_t rightCharge[4] = {0,0,0,0};
-
- Int_t mark = 0;
-
- filPed.AssignFormatted(ieffped[iT]); // no size-checking when using AssignFormatted; ieffped>=0
- filPed = filPed; // this checks the size
-
- ieffped[iT] = filPed.GetValue();
-
- for(Int_t i = 0; i<7; i++){
- selection[i].next = NULL;
- selection[i].iadc = -1; // value of -1: invalid adc
- selection[i].value = 0;
-
- }
- // selection[0] is starting list-element; just for pointing
-
- // loop over inner adc's
- for (Int_t iCol = 1; iCol < fNADC-1; iCol++) {
-
- Int_t left = fADCT[iCol-1][iT];
- Int_t center = fADCT[iCol][iT];
- Int_t right = fADCT[iCol+1][iT];
-
- Int_t sum = left + center + right; // cluster charge sum
- qsumAlu.AssignFormatted(sum);
- qsumAlu = qsumAlu; // size-checking; redundant
-
- qsum[iCol] = qsumAlu.GetValue();
-
- //hit detection and masking
- if(center>=left){
- if(center>right){
- if(qsum[iCol]>=(cTH + 3*ieffped[iT])){ // effective pedestal of all three channels must be added to cTH(+20); this is not parallel to TRAP manual; maybe cTH has to be adjusted in fFeeParam; therefore channels are not yet reduced by their pedestal
- mark |= 1; // marker
- }
- }
- }
- mark = mark<<1;
- }
- mark = mark>>1;
-
-
- // get selection of 6 adc's and sort,starting with greatest values
-
- //read three from right side and sort (primitive sorting algorithm)
- Int_t i = 0; // adc number
- Int_t j = 1; // selection number
- while(i<fNADC-2 && j<=3){
- i = i + 1;
- if( ((mark>>(i-1)) & 1) == 1) {
- selection[j].iadc = fNADC-1-i;
- selection[j].value = qsum[fNADC-1-i]>>6; // for hit-selection only the first 8 out of the 14 Bits are used for comparison
-
- // insert into sorted list
- listLeft = &selection[0];
- list = listLeft->next;
-
- if(list!=NULL) {
- while((list->next != NULL) && (selection[j].value <= list->value)){
- listLeft = list;
- list = list->next;
- }
-
- if(selection[j].value<=list->value){
- selection[j].next = list->next;
- list->next = &selection[j];
- }
- else {
- listLeft->next = &selection[j];
- selection[j].next = list;
- }
- }
- else{
- listLeft->next = &selection[j];
- selection[j].next = list;
- }
-
- j = j + 1;
- }
- }
+}
+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
- // read three from left side
- Int_t k = fNADC-2;
- while(k>i && j<=6) {
- if( ((mark>>(k-1)) & 1) == 1) {
- selection[j].iadc = fNADC-1-k;
- selection[j].value = qsum[fNADC-1-k]>>6;
-
- listLeft = &selection[0];
- list = listLeft->next;
-
- if(list!=NULL){
- while((list->next != NULL) && (selection[j].value <= list->value)){
- listLeft = list;
- list = list->next;
- }
-
- if(selection[j].value<=list->value){
- selection[j].next = list->next;
- list->next = &selection[j];
- }
- else {
- listLeft->next = &selection[j];
- selection[j].next = list;
- }
- }
- else{
- listLeft->next = &selection[j];
- selection[j].next = list;
- }
+ UShort_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
+ UShort_t val21s, val22s, val23s;
+ UShort_t idx23as, idx23bs;
+ UShort_t val23as, val23bs;
- j = j + 1;
- }
- k = k - 1;
- }
+ Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
+ &dummy1, &idx21s, &idx23as,
+ &dummy2, &val21s, &val23as);
- // get the four with greatest charge-sum
- list = &selection[0];
- for(i = 0; i<4; i++){
- if(list->next == NULL) continue;
- list = list->next;
- if(list->iadc == -1) continue;
- Int_t adc = list->iadc; // channel number with selected hit
-
- // the following arrays contain the four chosen channels in 1 time-bin
- portChannel[i] = adc;
- clusterCharge[i] = qsum[adc];
- leftCharge[i] = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
- centerCharge[i] = fADCT[adc][iT] - ieffped[iT];
- rightCharge[i] = fADCT[adc+1][iT] - ieffped[iT];
- }
+ Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
+ &dummy1, &idx22s, &idx23bs,
+ &dummy2, &val22s, &val23bs);
- // arithmetic unit
-
- // cluster verification
- if(!bVBY){
- for(i = 0; i<4; i++){
- Int_t lr = leftCharge[i]*rightCharge[i]*1024;
- Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
- if (lr>=cc){
- portChannel[i] = -1; // set to invalid address
- clusterCharge[i] = 0;
- }
- }
- }
+ Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
- // fit-sums of valid channels
- // local hit position
- for(i = 0; i<4; i++){
- if (centerCharge[i] == 0) {
- portChannel[i] = -1;
- }// prevent division by 0
-
- if (portChannel[i] == -1) continue;
-
- Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
-
- Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
- dCOG = (sign == kFALSE) ? dCOG : -dCOG; // AssignDouble doesn't allow for signed doubles
- dCOGAlu.AssignDouble(dCOG);
- Int_t iLUTpos = dCOGAlu.GetValue(); // steers position in LUT
-
- dCOG = dCOG/2;
- yrawAlu.AssignDouble(dCOG);
- Int_t iCOG = yrawAlu.GetValue();
- Int_t y = iCOG + fPosLUT[iLUTpos % 128]; // local position in pad-units
- yrawAlu.AssignFormatted(y); // 0<y<1
- yAlu = yrawAlu; // convert to 16 past-comma bits
-
- if(sign == kTRUE) yAlu.SetSign(-1); // buffer width of 9 bits; sign on real (not estimated) position
- xAlu.AssignInt(iT); // buffer width of 5 bits
-
-
- xxAlu = xAlu * xAlu; // buffer width of 10 bits -> fulfilled by x*x
-
- yyAlu = yAlu * yAlu; // buffer width of 16 bits
-
- xyAlu = xAlu * yAlu; // buffer width of 14 bits
-
- Int_t adc = portChannel[i]-1; // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
-
- // calculate fit-sums recursively
- // interpretation of their bit-length is given as comment
-
- // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
-
- XAlu.AssignFormatted(X[adc]);
- XAlu = XAlu + xAlu; // buffer width of 9 bits
- X[adc] = XAlu.GetValue();
-
- XXAlu.AssignFormatted(XX[adc]);
- XXAlu = XXAlu + xxAlu; // buffer width of 14 bits
- XX[adc] = XXAlu.GetValue();
-
- if (Y[adc] < 0) {
- YAlu.AssignFormatted(-Y[adc]); // make sure that only positive values are assigned; sign-setting must be done by hand
- YAlu.SetSign(-1);
- }
- else {
- YAlu.AssignFormatted(Y[adc]);
- YAlu.SetSign(1);
- }
-
- YAlu = YAlu + yAlu; // buffer width of 14 bits (8 past-comma);
- Y[adc] = YAlu.GetSignedValue();
-
- YYAlu.AssignFormatted(YY[adc]);
- YYAlu = YYAlu + yyAlu; // buffer width of 21 bits (16 past-comma)
- YY[adc] = YYAlu.GetValue();
-
- if (XY[adc] < 0) {
- XYAlu.AssignFormatted(-XY[adc]);
- XYAlu.SetSign(-1);
- }
- else {
- XYAlu.AssignFormatted(XY[adc]);
- XYAlu.SetSign(1);
- }
+ Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
+ &dummy1, &dummy2, idx6o,
+ &dummy3, &dummy4, &dummy5);
+}
- XYAlu = XYAlu + xyAlu; // buffer allows 17 bits (8 past-comma)
- XY[adc] = XYAlu.GetSignedValue();
-
- N[adc] = N[adc] + 1;
-
-
- // accumulated charge
- qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
- qtruncAlu = qsumAlu;
-
- if(iT>=tQS0 && iT<=tQE0){
- QT0Alu.AssignFormatted(QT0[adc]);
- QT0Alu = QT0Alu + qtruncAlu;
- QT0[adc] = QT0Alu.GetValue();
- //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
- }
-
- if(iT>=tQS1 && iT<=tQE1){
- QT1Alu.AssignFormatted(QT1[adc]);
- QT1Alu = QT1Alu + qtruncAlu;
- QT1[adc] = QT1Alu.GetValue();
- //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
- }
- }// i
-
- // remapping is done!!
-
- }//iT
-
-
-
- // tracklet-assembly
-
- // put into AliTRDfeeParam and take care that values are in proper range
- const Int_t cTCL = 1; // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??)
- const Int_t cTCT = 8; // joint number of hits; 8<=TCT<=31; note that according to TRAP manual this number cannot be lower than 8; however it should be adjustable to the number of hits in the fit time range (40%)
-
- Int_t mPair = 0; // marker for possible tracklet pairs
- Int_t* hitSum = new Int_t[fNADC-3];
- // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18;
-
- // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
- for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
- hitSum[iCol] = N[iCol] + N[iCol+1];
- if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
- mPair |= 1; // mark as possible channel-pair
-
- }
- mPair = mPair<<1;
- }
- mPair = mPair>>1;
-
- List_t* selectPair = new List_t[fNADC-2]; // list with 18 elements (0..18) containing the left channel-nr and hit sums
- // selectPair[18] is starting list-element just for pointing
- for(Int_t k = 0; k<fNADC-2; k++){
- selectPair[k].next = NULL;
- selectPair[k].iadc = -1; // invalid adc
- selectPair[k].value = 0;
-
- }
- list = NULL;
- listLeft = NULL;
-
- // read marker and sort according to hit-sum
-
- Int_t adcL = 0; // left adc-channel-number (remapped)
- Int_t selNr = 0; // current number in list
-
- // insert marked channels into list and sort according to hit-sum
- while(adcL < fNADC-3 && selNr < fNADC-3){
-
- if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
- selectPair[selNr].iadc = adcL;
- selectPair[selNr].value = hitSum[adcL];
-
- listLeft = &selectPair[fNADC-3];
- list = listLeft->next;
-
- if(list!=NULL) {
- while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
- listLeft = list;
- list = list->next;
- }
-
- if(selectPair[selNr].value <= list->value){
- selectPair[selNr].next = list->next;
- list->next = &selectPair[selNr];
- }
- else {
- listLeft->next = &selectPair[selNr];
- selectPair[selNr].next = list;
- }
-
- }
- else{
- listLeft->next = &selectPair[selNr];
- selectPair[selNr].next = list;
- }
-
- selNr = selNr + 1;
- }
- adcL = adcL + 1;
- }
-
- //select up to 4 channels with maximum number of hits
- Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
- Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
- list = &selectPair[fNADC-3];
-
- for (Int_t i = 0; i<4; i++) {
- if(list->next == NULL) continue;
- list = list->next;
- if(list->iadc == -1) continue;
- lpairChannel[i] = list->iadc; // channel number with selected hit
- rpairChannel[i] = lpairChannel[i]+1;
- }
-
- // avoid submission of double tracklets
- for (Int_t i = 3; i>0; i--) {
- for (Int_t j = i-1; j>-1; j--) {
- if(lpairChannel[i] == rpairChannel[j]) {
- lpairChannel[i] = -1;
- rpairChannel[i] = -1;
- break;
- }
- /* if(rpairChannel[i] == lpairChannel[j]) {
- lpairChannel[i] = -1;
- rpairChannel[i] = -1;
- break;
- }*/
- }
- }
-
- // merging of the fit-sums of the remainig channels
- // assume same data-word-width as for fit-sums for 1 channel
- // relative scales!
- Int_t mADC[4];
- Int_t mN[4];
- Int_t mQT0[4];
- Int_t mQT1[4];
- Int_t mX[4];
- Int_t mXX[4];
- Int_t mY[4];
- Int_t mYY[4];
- Int_t mXY[4];
- Int_t mOffset[4];
- Int_t mSlope[4];
- Int_t mMeanCharge[4];
- Int_t inverseN = 0;
- Double_t invN = 0;
- Int_t one = 0;
-
- for (Int_t i = 0; i<4; i++){
- mADC[i] = -1; // set to invalid number
- mN[i] = 0;
- mQT0[i] = 0;
- mQT1[i] = 0;
- mX[i] = 0;
- mXX[i] = 0;
- mY[i] = 0;
- mYY[i] = 0;
- mXY[i] = 0;
- mOffset[i] = 0;
- mSlope[i] = 0;
- mMeanCharge[i] = 0;
- }
-
- oneAlu.AssignInt(1);
- one = oneAlu.GetValue(); // one with 8 past comma bits
-
- for (Int_t i = 0; i<4; i++){
-
-
- mADC[i] = lpairChannel[i]; // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
- // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
- Int_t madc = mADC[i];
- if (madc == -1) continue;
-
- YAlu.AssignInt(N[rpairChannel[i]]);
- Int_t wpad = YAlu.GetValue(); // enlarge hit counter of right channel by 8 past-comma bits; YAlu can have 5 pre-comma bits (values up to 63); hit counter<=nr of time bins (24)
-
- mN[i] = hitSum[madc];
-
- // don't merge fit sums in case of a stand-alone tracklet (consisting of only 1 channel); in that case only left channel makes up the fit sums
- if (N[madc+1] == 0) {
- mQT0[i] = QT0[madc];
- mQT1[i] = QT1[madc];
-
- }
- else {
-
- // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
-
- mQT0[i] = QT0[madc] + QT0[madc+1];
- QT0Alu.AssignFormatted(mQT0[i]);
- QT0Alu = QT0Alu; // size-check
- mQT0[i] = QT0Alu.GetValue(); // write back
-
- mQT1[i] = QT1[madc] + QT1[madc+1];
- QT1Alu.AssignFormatted(mQT1[i]);
- QT1Alu = QT1Alu;
- mQT1[i] = QT1Alu.GetValue();
- }
-
- // calculate the mean charge in adc values; later to be replaced by electron likelihood
- mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
- mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
- // simulate LUT for 1/N; LUT is fed with the double-accurate pre-calculated value of 1/N; accuracy of entries has to be adjusted to real TRAP
- invN = 1.0/(mN[i]);
- inverseNAlu.AssignDouble(invN);
- inverseN = inverseNAlu.GetValue();
- mMeanCharge[i] = mMeanCharge[i] * inverseN; // now to be interpreted with 8 past-comma bits
- TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
- TotalChargeAlu = TotalChargeAlu;
- MeanChargeAlu = TotalChargeAlu;
- mMeanCharge[i] = MeanChargeAlu.GetValue();
-
- // this check is not necessary; it is just for efficiency reasons
- if (N[madc+1] == 0) {
- mX[i] = X[madc];
- mXX[i] = XX[madc];
- mY[i] = Y[madc];
- mXY[i] = XY[madc];
- mYY[i] = YY[madc];
- }
- else {
-
- mX[i] = X[madc] + X[madc+1];
- XAlu.AssignFormatted(mX[i]);
- XAlu = XAlu;
- mX[i] = XAlu.GetValue();
-
- mXX[i] = XX[madc] + XX[madc+1];
- XXAlu.AssignFormatted(mXX[i]);
- XXAlu = XXAlu;
- mXX[i] = XXAlu.GetValue();
-
-
- mY[i] = Y[madc] + Y[madc+1] + wpad;
- if (mY[i] < 0) {
- YAlu.AssignFormatted(-mY[i]);
- YAlu.SetSign(-1);
- }
- else {
- YAlu.AssignFormatted(mY[i]);
- YAlu.SetSign(1);
- }
- YAlu = YAlu;
- mY[i] = YAlu.GetSignedValue();
-
- mXY[i] = XY[madc] + XY[madc+1] + X[madc+1]*one; // multiplication by one to maintain the data format
-
- if (mXY[i] < 0) {
- XYAlu.AssignFormatted(-mXY[i]);
- XYAlu.SetSign(-1);
- }
- else {
- XYAlu.AssignFormatted(mXY[i]);
- XYAlu.SetSign(1);
- }
- XYAlu = XYAlu;
- mXY[i] = XYAlu.GetSignedValue();
-
- mYY[i] = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
- if (mYY[i] < 0) {
- YYAlu.AssignFormatted(-mYY[i]);
- YYAlu.SetSign(-1);
- }
- else {
- YYAlu.AssignFormatted(mYY[i]);
- YYAlu.SetSign(1);
- }
-
- YYAlu = YYAlu;
- mYY[i] = YYAlu.GetSignedValue();
- }
-
- }
-
- // calculation of offset and slope from the merged fit-sums;
- // YY is needed for some error measure only; still to be done
- // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
-
- // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- // !!important note: the offset is calculated from hits in the time bin range between tFS and tFE; it corresponds to the value at the height of the time bin tFS which does NOT need to correspond to the upper side of the drift !!
- // !!volume (cathode wire plane). The offset cannot be rescaled as long as it is unknown which is the first time bin that contains hits from the drift region and thus to which distance from the cathode plane tFS corresponds. !!
- // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below). !!
- // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used. !!
- // !!The corrections to the offset (e.g. no ExB correction applied as offset is supposed to be on top of drift region; however not at anode wire, so some inclination of drifting clusters due to Lorentz angle exists) are only !!
- // !!valid (in approximation) if tFS is close to the beginning of the drift region. !!
- // !!The slope however can be converted to a deflection length between electrode and cathode wire plane as it is clear that the drift region is sampled 20 times !!
- // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- // which formats should be chosen?
- AliTRDtrapAlu denomAlu;
- denomAlu.Init(20,8);
- AliTRDtrapAlu numAlu;
- numAlu.Init(20,8);
- // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
- // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
-
- for (Int_t i = 0; i<4; i++) {
- if (mADC[i] == -1) continue;
-
- Int_t num0 = (mN[i]*mXX[i]-mX[i]*mX[i]);
- if (num0 < 0) {
- denomAlu.AssignInt(-num0); // num0 does not have to be interpreted as having past-comma bits -> AssignInt
- denomAlu.SetSign(-1);
- }
- else {
- denomAlu.AssignInt(num0);
- denomAlu.SetSign(1);
- }
-
- Int_t num1 = mN[i]*mXY[i] - mX[i]*mY[i];
- if (num1 < 0) {
- numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
- numAlu.SetSign(-1);
- }
- else {
- numAlu.AssignFormatted(num1);
- numAlu.SetSign(1);
- }
- numAlu = numAlu/denomAlu;
- mSlope[i] = numAlu.GetSignedValue();
-
- Int_t num2 = mXX[i]*mY[i] - mX[i]*mXY[i];
-
- if (num2 < 0) {
- numAlu.AssignFormatted(-num2);
- numAlu.SetSign(-1);
- }
- else {
- numAlu.AssignFormatted(num2);
- numAlu.SetSign(1);
- }
-
- numAlu = numAlu/denomAlu;
-
-
- mOffset[i] = numAlu.GetSignedValue();
- numAlu.SetSign(1);
- denomAlu.SetSign(1);
-
-
- //numAlu.AssignInt(mADC[i]+1); // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)
- numAlu.AssignDouble((Double_t)mADC[i] + 1.5); // numAlu has enough pre-comma place for that; correct trafo, best values
- mOffset[i] = mOffset[i] + numAlu.GetValue(); // transform offset to a coord.system relative to chip; +1 to avoid neg. values
-
- // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
- // reverse adc-counting to be again in line with the online mapping
- mADC[i] = fNADC - 4 - mADC[i]; // fNADC-4-mADC[i]: 0..17; remapping necessary;
- mADC[i] = mADC[i] + 2;
- // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
- }
+// ----- I/O implementation -----
- // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected;
- // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
-
- // transform parameters to the local coordinate-system of a stack (used by GTU)
- AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
-
- Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
- //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
-
- // difference between width of inner and outer pads of a row is not accounted for;
-
- Double_t magField = 0.4; // z-component of magnetic field in Tesla; adjust to current simulation!!; magnetic field can hardly be evaluated for the position of each mcm
- Double_t eCharge = 0.3; // unit charge in (GeV/c)/m*T
- Double_t ptMin = 2.3; // minimum transverse momentum (GeV/c); to be adjusted(?)
-
- Double_t granularityOffset = 0.160; // granularity for offset in mm
- Double_t granularitySlope = 0.140; // granularity for slope in mm
-
- // get the coordinates in SM-system; parameters:
-
- Double_t zPos = (padPlane->GetRowPos(fRow))*10.0; // z-position of the MCM; fRow is counted on a chamber; SM consists of 5
- // zPos is position of pad-borders;
- Double_t zOffset = 0.0;
- if ( fRow == 0 || fRow == 15 ) {
- zOffset = padPlane->GetLengthOPad();
- }
- else {
- zOffset = padPlane->GetLengthIPad();
- }
- zOffset = (-1.0) * zOffset/2.0;
- // turn zPos to be z-coordinate at middle of pad-row
- zPos = zPos + zOffset;
-
-
- Double_t xPos = 0.0; // x-position of the upper border of the drift-chamber of actual layer
- Int_t icol = 0; // column-number of adc-channel
- Double_t yPos[4]; // y-position of the pad to which ADC is connected
- Double_t dx = 30.0; // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
- Double_t freqSample = fFeeParam->GetSamplingFrequency(); // retrieve the sampling frequency (10.019750 MHz)
- Double_t vdrift = fCal->GetVdriftAverage(fChaId); // averaged drift velocity for this detector (1.500000 cm/us)
- Int_t nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
- Int_t nrOfAmplTimeBins = 2; // the number of time bins between anode wire and cathode wires in ampl.region (3.5mm)(guess)(suppose v_drift+3.5cm/us there=>all clusters arrive at anode wire within one time bin (100ns))
- Int_t nrOfOffsetCorrTimeBins = tFS - nrOfAmplTimeBins - 1; // -1 is to be conservative; offset correction will not remove the shift but is supposed to improve it; if tFS = 5, 2 drift time bins before tFS are assumed
- if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed
- Double_t lorTan = AliTRDCommonParam::Instance()->GetOmegaTau(vdrift); // tan of the Lorentz-angle for this detector; could be evaluated and set as a parameter for each mcm
- //Double_t lorAngle = 7.0; // Lorentz-angle in degrees
- Double_t tiltAngle = padPlane->GetTiltingAngle(); // sign-respecting tilting angle of pads in actual layer
- Double_t tiltTan = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
- //Double_t lorTan = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
-
- Double_t alphaMax[4]; // maximum deflection from the direction to the primary vertex; granularity of hit pads
- Double_t slopeMin[4]; // local limits for the deflection
- Double_t slopeMax[4];
- Int_t mslopeMin[4]; // in granularity units; to be compared to mSlope[i]
- Int_t mslopeMax[4];
-
-
- // x coord. of upper side of drift chambers in local SM-system (in mm)
- // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
- // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
- switch(fLayer)
- {
- case 0:
- xPos = 3003.0;
- break;
- case 1:
- xPos = 3129.0;
- break;
- case 2:
- xPos = 3255.0;
- break;
- case 3:
- xPos = 3381.0;
- break;
- case 4:
- xPos = 3507.0;
- break;
- case 5:
- xPos = 3633.0;
- break;
- }
-
- // calculation of offset-correction n:
-
- Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));
-
- nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
- if (nCorrectOffset < 0) {
- numAlu.AssignInt(-nCorrectOffset);
- numAlu.SetSign(-1);
- }
- else {
- numAlu.AssignInt(nCorrectOffset);
- numAlu.SetSign(1);
- }
- nCorrectOffset = numAlu.GetSignedValue();
+ostream& AliTRDmcmSim::Text(ostream& os)
+{
+ // manipulator to activate output in text format (default)
- // the Lorentz correction to the offset
- Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
-
+ os.iword(fgkFormatIndex) = 0;
+ return os;
+}
- lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
-
- if(lorCorrectOffset < 0) {
- numAlu.AssignDouble(-lorCorrectOffset);
- numAlu.SetSign(-1);
- }
- else{
- numAlu.AssignDouble(lorCorrectOffset);
- numAlu.SetSign(1);
- }
-
- Int_t mlorCorrectOffset = numAlu.GetSignedValue();
-
-
- Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
-
- // calculation of slope-correction
-
- // this is only true for tracks coming (approx.) from primary vertex
- // everything is evaluated for a tracklet covering the whole drift chamber
- Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
- // Double_t cCorrectSlope = zPos/xPos*dx*tiltTan/granularitySlope;
- // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
-
- if (cCorrectSlope < 0) {
- numAlu.AssignDouble(-cCorrectSlope);
- numAlu.SetSign(-1);
- }
- else {
- numAlu.AssignDouble(cCorrectSlope);
- numAlu.SetSign(1);
- }
- cCorrectSlope = numAlu.GetSignedValue();
-
- // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
- // different pad-width of outer pads of a pad-plane not taken into account
- // note that the fit was only done in the range tFS to tFE, however this range does not need to cover the whole drift region (neither start nor end of it)
- // however the tracklets are supposed to be a fit in the drift region thus the linear function is stretched to fit the drift region of 30 mm
-
-
- Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope; // >= 0.0
-
- AliTRDtrapAlu correctAlu;
- correctAlu.Init(20,8);
-
- AliTRDtrapAlu offsetAlu;
- offsetAlu.Init(13,0,-0x1000,0x0FFF); // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
-
- AliTRDtrapAlu slopeAlu;
- slopeAlu.Init(7,0,-0x40,0x3F); // 7 bit-word; 2-complement (1 sign-bit);
-
- for (Int_t i = 0; i<4; i++) {
-
- if (mADC[i] == -1) continue;
-
- icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
- yPos[i] = (padPlane->GetColPos(icol))*10.0;
-
-
- // offset:
-
- correctAlu.AssignDouble(mCorrectOffset); // done because max. accuracy is 8 bit
- mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
- mOffset[i] = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset));
- //mOffset[i] = mOffset[i]*(-1); // adjust to direction of y-axes in online simulation
-
- if (mOffset[i] < 0) {
- numAlu.AssignFormatted(-mOffset[i]);
- numAlu.SetSign(-1);
- }
- else {
- numAlu.AssignFormatted(mOffset[i]);
- numAlu.SetSign(1);
- }
+ostream& AliTRDmcmSim::Cfdat(ostream& os)
+{
+ // manipulator to activate output in CFDAT format
+ // to send to the FEE via SCSN
- offsetAlu = numAlu;
- mOffset[i] = offsetAlu.GetSignedValue();
+ os.iword(fgkFormatIndex) = 1;
+ return os;
+}
-
- // slope:
-
- correctAlu.AssignDouble(mCorrectSlope);
- mCorrectSlope = correctAlu.GetValueWhole();
-
- mSlope[i] = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
+ostream& AliTRDmcmSim::Raw(ostream& os)
+{
+ // manipulator to activate output as raw data dump
- if (mSlope[i] < 0) {
- numAlu.AssignFormatted(-mSlope[i]);
- numAlu.SetSign(-1);
- }
- else {
- numAlu.AssignFormatted(mSlope[i]);
- numAlu.SetSign(1);
- }
+ os.iword(fgkFormatIndex) = 2;
+ return os;
+}
- slopeAlu = numAlu; // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
- mSlope[i] = slopeAlu.GetSignedValue();
-
- // local (LTU) limits for the deflection
- // ATan returns angles in radian
- alphaMax[i] = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
- slopeMin[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
- slopeMax[i] = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
-
- if (slopeMin[i] < 0) {
- slopeAlu.AssignDouble(-slopeMin[i]);
- slopeAlu.SetSign(-1);
- }
- else {
- slopeAlu.AssignDouble(slopeMin[i]);
- slopeAlu.SetSign(1);
- }
- mslopeMin[i] = slopeAlu.GetSignedValue(); // the borders should lie inside the range of mSlope -> usage of slopeAlu again
-
- if (slopeMax[i] < 0) {
- slopeAlu.AssignDouble(-slopeMax[i]);
- slopeAlu.SetSign(-1);
- }
- else {
- slopeAlu.AssignDouble(slopeMax[i]);
- slopeAlu.SetSign(1);
+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;
}
- mslopeMax[i] = slopeAlu.GetSignedValue();
}
- // suppress submission of tracks with low stiffness
- // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber)
-
- // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
- // up to now they are sorted according to maximum hit sum
- // is the sorting really done in the TRAP-chip?
-
- Int_t order[4] = {-1,-1,-1,-1};
- Int_t wordnr = 0; // number of tracklet-words
-
- for(Int_t j = 0; j < fMaxTracklets; j++) {
- //if( mADC[j] == -1) continue;
- if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
- wordnr++;
- if( wordnr-1 == 0) {
- order[0] = j;
- continue;
- }
- // wordnr-1>0, wordnr-1<4
- order[wordnr-1] = j;
- for( Int_t k = 0; k < wordnr-1; k++) {
- if( mOffset[j] < mOffset[order[k]] ) {
- for( Int_t l = wordnr-1; l > k; l-- ) {
- order[l] = order[l-1];
- }
- order[k] = j;
- break;
- }
-
+ // ----- 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;
}
- }
-
- // fill the bit-words in ascending order and without gaps
- UInt_t bitWord[4] = {0,0,0,0}; // attention: unsigned int to have real 32 bits (no 2-complement)
- for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
- //Bool_t rem1 = kTRUE;
-
- Int_t i = order[j];
- //bit-word is 2-complement and therefore without sign
- bitWord[j] = 1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
- //printf("\n");
- UInt_t shift = 0;
- UInt_t shift2 = 0;
-
-
-
-
- /*printf("mean charge: %d\n",mMeanCharge[i]);
- printf("row: %d\n",fRow);
- printf("slope: %d\n",mSlope[i]);
- printf("pad position: %d\n",mOffset[i]);
- printf("channel: %d\n",mADC[i]);*/
-
- // electron probability (currently not implemented; the mean charge is just scaled)
- shift = (UInt_t)mMeanCharge[i];
- for(Int_t iBit = 0; iBit < 8; iBit++) {
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= (shift>>(7-iBit))&1;
- //printf("0");
+ os << std::endl;
}
+ }
- // pad row
- shift = (UInt_t)fRow;
- for(Int_t iBit = 0; iBit < 4; iBit++) {
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= (shift>>(3-iBit))&1;
- //printf("%d", (fRow>>(3-iBit))&1);
- }
-
- // deflection length
- if(mSlope[i] < 0) {
- shift = (UInt_t)(-mSlope[i]);
- // shift2 is 2-complement of shift
- shift2 = 1;
- for(Int_t iBit = 1; iBit < 7; iBit++) {
- shift2 = shift2<<1;
- shift2 |= (1- (((shift)>>(6-iBit))&1) );
- //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
- }
- shift2 = shift2 + 1;
- //printf("1");
- for(Int_t iBit = 0; iBit < 7; iBit++) {
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= (shift2>>(6-iBit))&1;
- //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
- }
- }
- else {
- shift = (UInt_t)(mSlope[i]);
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= 0;
- //printf("0");
- for(Int_t iBit = 1; iBit < 7; iBit++) {
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= (shift>>(6-iBit))&1;
- //printf("%d",(mSlope[i]>>(6-iBit))&1);
- }
- }
+ // ----- raw data ouptut -----
+ else if (os.iword(AliTRDmcmSim::fgkFormatIndex) == 2) {
+ Int_t bufSize = 300;
+ UInt_t *buf = new UInt_t[bufSize];
- // pad position
- if(mOffset[i] < 0) {
- shift = (UInt_t)(-mOffset[i]);
- shift2 = 1;
- for(Int_t iBit = 1; iBit < 13; iBit++) {
- shift2 = shift2<<1;
- shift2 |= (1-(((shift)>>(12-iBit))&1));
- //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
- }
- shift2 = shift2 + 1;
- //printf("1");
- for(Int_t iBit = 0; iBit < 13; iBit++) {
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= (shift2>>(12-iBit))&1;
- //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
- }
- }
- else {
- shift = (UInt_t)mOffset[i];
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= 0;
- //printf("0");
- for(Int_t iBit = 1; iBit < 13; iBit++) {
- bitWord[j] = bitWord[j]<<1;
- bitWord[j] |= (shift>>(12-iBit))&1;
- //printf("%d",(mOffset[i]>>(12-iBit))&1);
- }
- }
+ 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;
-
- //printf("bitWord: %u\n",bitWord[j]);
- //printf("adc: %d\n",mADC[i]);
- fMCMT[j] = bitWord[j];
- }
-
- //printf("\n");
-
-
- delete [] qsum;
- delete [] ieffped;
-
- delete [] X;
- delete [] XX;
- delete [] Y;
- delete [] YY;
- delete [] XY;
- delete [] N;
- delete [] QT0;
- delete [] QT1;
-
- delete [] hitSum;
- delete [] selectPair;
-
- delete padPlane;
-
-//if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
-
- if (!fFeeParam->GetMCTrackletOutput())
- return;
-
- AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
- AliLog::SetFileOutput("../log/tracklet.log");
-
- // testing for wordnr in order to speed up the simulation
- if (wordnr == 0)
- return;
-
- UInt_t *trackletWord = new UInt_t[fMaxTracklets];
- Int_t *adcChannel = new Int_t[fMaxTracklets];
- Int_t *trackRef = new Int_t[fMaxTracklets];
-
- Int_t u = 0;
-
- AliTRDdigitsManager *digman = new AliTRDdigitsManager();
- digman->ReadDigits(AliRunLoader::Instance()->GetLoader("TRDLoader")->TreeD());
- digman->SetUseDictionaries(kTRUE);
- AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
-
- for (Int_t j = 0; j < fMaxTracklets; j++) {
- Int_t i = order[j];
- trackletWord[j] = 0;
- adcChannel[j] = -1;
- if (bitWord[j]!=0) {
- trackletWord[u] = bitWord[j];
- adcChannel[u] = mADC[i]; // mapping onto the original adc-array to be in line with the digits-adc-ordering (21 channels in total on 1 mcm, 18 belonging to pads); mADC[i] should be >-1 in case bitWord[i]>0
-
-// Finding label of MC track
- TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
- Int_t track[3];
- Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
- Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
- Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
- Int_t det = 30 * fSector + 6 * fStack + fLayer;
- for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
- track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
- track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
- track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
- hTrkRef->Fill(track[0]);
- if (track[1] != track[0] && track[1] != -1)
- hTrkRef->Fill(track[1]);
- if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
- hTrkRef->Fill(track[2]);
- if (padcol_ngb >= 0) {
- track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
- track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
- track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
- hTrkRef->Fill(track[0]);
- if (track[1] != track[0] && track[1] != -1)
- hTrkRef->Fill(track[1]);
- if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
- hTrkRef->Fill(track[2]);
- }
- }
- trackRef[u] = hTrkRef->GetMaximumBin() - 1;
- delete hTrkRef;
- u = u + 1;
- }
+ delete [] buf;
}
- AliDataLoader *dl = AliRunLoader::Instance()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
- if (!dl) {
- AliError("Could not get the tracklets data loader!");
- }
else {
- TTree *trackletTree = dl->Tree();
- if (!trackletTree)
- dl->MakeTree();
- trackletTree = dl->Tree();
-
- AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM();
- TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
- if (!trkbranch)
- trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
- trkbranch->SetAddress(&trkl);
-
- for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
- if (trackletWord[iTracklet] == 0)
- continue;
- trkl->SetTrackletWord(trackletWord[iTracklet]);
- trkl->SetDetector(30*fSector + 6*fStack + fLayer);
- trkl->SetROB(fRobPos);
- trkl->SetMCM(fMcmPos);
- trkl->SetLabel(trackRef[iTracklet]);
- trackletTree->Fill();
- }
- delete trkl;
- dl->WriteData("OVERWRITE");
+ os << "unknown format set" << std::endl;
}
- delete [] trackletWord;
- delete [] adcChannel;
- delete [] trackRef;
- delete digman;
-
- // to be done:
- // error measure for quality of fit (not necessarily needed for the trigger)
- // cluster quality threshold (not yet set)
- // electron probability
+ return os;
}
-//_____________________________________________________________________________________
-void AliTRDmcmSim::GeneratefZSM1Dim()
-{
- //
- // Generate the array fZSM1Dim necessary
- // for the method ProduceRawStream
- //
- // Fill the mapping
- // Supressed zeros indicated by -1 in digits array
- for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
- {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ )
- {
-
- if(fADCF[iadc][it]==-1) // If is a supressed value
- {
- fZSM[iadc][it]=1;
- }
- else // Not suppressed
- {
- fZSM[iadc][it]=0;
- }
- }
- }
- // Make the 1 dim projection
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
- {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ )
- {
- fZSM1Dim[iadc] &= fZSM[iadc][it];
- }
- }
-}
-//_______________________________________________________________________________________
-void AliTRDmcmSim::CopyArrays()
+void AliTRDmcmSim::PrintFitRegXml(ostream& os) const
{
- //
- // Initialize filtered data array with raw data
- // Method added for internal consistency
- //
-
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
- {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ )
- {
- fADCF[iadc][it] = fADCR[iadc][it];
+ // 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 << "<nginject>" << std::endl;
+ os << "<ack roc=\""<< fDetector << "\" cmndid=\"0\">" << std::endl;
+ os << "<dmem-readout>" << std::endl;
+ os << "<d det=\"" << fDetector << "\">" << std::endl;
+ os << " <ro-board rob=\"" << fRobPos << "\">" << std::endl;
+ os << " <m mcm=\"" << fMcmPos << "\">" << std::endl;
+
+ for(int cpu=0; cpu<4; cpu++) {
+ os << " <c cpu=\"" << cpu << "\">" << std::endl;
+ if(fFitPtr[cpu] != 31) {
+ for(int adcch=fFitPtr[cpu]; adcch<fFitPtr[cpu]+2; adcch++) {
+ os << " <ch chnr=\"" << adcch << "\">"<< std::endl;
+ os << " <hits>" << fFitReg[adcch].fNhits << "</hits>"<< std::endl;
+ os << " <q0>" << fFitReg[adcch].fQ0/4 << "</q0>"<< std::endl; // divided by 4 because in simulation we have 2 additional decimal places
+ os << " <q1>" << fFitReg[adcch].fQ1/4 << "</q1>"<< std::endl; // in the output
+ os << " <sumx>" << fFitReg[adcch].fSumX << "</sumx>"<< std::endl;
+ os << " <sumxsq>" << fFitReg[adcch].fSumX2 << "</sumxsq>"<< std::endl;
+ os << " <sumy>" << fFitReg[adcch].fSumY << "</sumy>"<< std::endl;
+ os << " <sumysq>" << fFitReg[adcch].fSumY2 << "</sumysq>"<< std::endl;
+ os << " <sumxy>" << fFitReg[adcch].fSumXY << "</sumxy>"<< std::endl;
+ os << " </ch>" << std::endl;
+ }
}
- }
+ os << " </c>" << std::endl;
+ }
+ os << " </m>" << std::endl;
+ os << " </ro-board>" << std::endl;
+ os << "</d>" << std::endl;
+ os << "</dmem-readout>" << std::endl;
+ os << "</ack>" << std::endl;
+ os << "</nginject>" << std::endl;
+ }
}
-//_______________________________________________________________________________________
-void AliTRDmcmSim::StartfastZS(Int_t pads, Int_t timebins)
+
+
+void AliTRDmcmSim::PrintTrackletsXml(ostream& os) const
{
- //
- // Initialize just the necessary elements to perform
- // the zero suppression in the digitizer
- //
-
- fFeeParam = AliTRDfeeParam::Instance();
- fSimParam = AliTRDSimParam::Instance();
- fNADC = pads;
- fNTimeBin = timebins;
+ // print tracklets in XML format
+
+ os << "<nginject>" << std::endl;
+ os << "<ack roc=\""<< fDetector << "\" cmndid=\"0\">" << std::endl;
+ os << "<dmem-readout>" << std::endl;
+ os << "<d det=\"" << fDetector << "\">" << std::endl;
+ os << " <ro-board rob=\"" << fRobPos << "\">" << std::endl;
+ os << " <m mcm=\"" << fMcmPos << "\">" << 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 ) ;
- if( fADCR == NULL )
- {
- fADCR = new Int_t *[fNADC];
- fADCF = new Int_t *[fNADC];
- fADCT = new Int_t *[fNADC];
- fZSM = new Int_t *[fNADC];
- fZSM1Dim = new Int_t [fNADC];
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
- {
- fADCR[iadc] = new Int_t[fNTimeBin];
- fADCF[iadc] = new Int_t[fNTimeBin];
- fADCT[iadc] = new Int_t[fNTimeBin];
- fZSM [iadc] = new Int_t[fNTimeBin];
}
- }
+ os << " <trk> <pid>" << pid << "</pid>" << " <padrow>" << padrow << "</padrow>"
+ << " <slope>" << slope << "</slope>" << " <offset>" << offset << "</offset>" << "</trk>" << std::endl;
+ }
+
+ os << " </m>" << std::endl;
+ os << " </ro-board>" << std::endl;
+ os << "</d>" << std::endl;
+ os << "</dmem-readout>" << std::endl;
+ os << "</ack>" << std::endl;
+ os << "</nginject>" << std::endl;
+}
- for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
- {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ )
- {
- fADCR[iadc][it] = 0;
- fADCF[iadc][it] = 0;
- fADCT[iadc][it] = -1;
- fZSM [iadc][it] = 1;
- }
- fZSM1Dim[iadc] = 1;
- }
-
- fInitialized = kTRUE;
+
+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::FlagDigitsArray(AliTRDarrayADC *tempdigs, Int_t valrow)
+
+
+void AliTRDmcmSim::PrintAdcDatXml(ostream& os) const
{
- //
- // Modify the digits array to flag suppressed values
- //
+ // print ADC data in XML format
+
+ os << "<nginject>" << std::endl;
+ os << "<ack roc=\""<< fDetector << "\" cmndid=\"0\">" << std::endl;
+ os << "<dmem-readout>" << std::endl;
+ os << "<d det=\"" << fDetector << "\">" << std::endl;
+ os << " <ro-board rob=\"" << fRobPos << "\">" << std::endl;
+ os << " <m mcm=\"" << fMcmPos << "\">" << std::endl;
+
+ for(Int_t iChannel = 0; iChannel < fgkNADC; iChannel++) {
+ os << " <ch chnr=\"" << iChannel << "\">" << std::endl;
+ for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
+ os << "<tb>" << fADCF[iChannel][iTimeBin]/4 << "</tb>";
+ }
+ os << " </ch>" << std::endl;
+ }
+
+ os << " </m>" << std::endl;
+ os << " </ro-board>" << std::endl;
+ os << "</d>" << std::endl;
+ os << "</dmem-readout>" << std::endl;
+ os << "</ack>" << std::endl;
+ os << "</nginject>" << std::endl;
+}
- for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
- {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ )
- {
- if(fZSM[iadc][it]==1)
- {
- tempdigs->SetData(valrow,iadc,it,-1);
- }
- }
- }
+
+
+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::RestoreZeros()
+
+
+void AliTRDmcmSim::PrintPidLutHuman()
{
- //
- // Restore the zero-suppressed values (set as -1) to the value 0
- //
+ // print PID LUT in human readable format
- for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ )
- {
- for( Int_t it = 0 ; it < fNTimeBin ; it++ )
- {
-
- if(fADCF[iadc][it]==-1) //if is a supressed zero, reset to zero
- {
- fADCF[iadc][it]=0;
- fADCR[iadc][it]=0;
- }
- }
- }
+ 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;
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff