#include "AliTRDcalibDB.h"
#include "AliTRDdigitsManager.h"
#include "AliTRDarrayADC.h"
+#include "AliTRDarrayDictionary.h"
#include "AliTRDpadPlane.h"
#include "AliTRDtrackletMCM.h"
#include "AliTRDmcmSim.h"
+#include "AliMagF.h"
+#include "TGeoGlobalMagField.h"
+
ClassImp(AliTRDmcmSim)
+Bool_t AliTRDmcmSim::fgApplyCut = kTRUE;
+
//_____________________________________________________________________________
AliTRDmcmSim::AliTRDmcmSim() : TObject()
,fInitialized(kFALSE)
,fFeeParam(NULL)
,fTrapConfig(NULL)
,fSimParam(NULL)
+ ,fCommonParam(NULL)
,fCal(NULL)
,fGeo(NULL)
+ ,fDigitsManager(NULL)
,fPedAcc(NULL)
,fGainCounterA(NULL)
,fGainCounterB(NULL)
fFeeParam = AliTRDfeeParam::Instance();
fTrapConfig = AliTRDtrapConfig::Instance();
fSimParam = AliTRDSimParam::Instance();
+ fCommonParam = AliTRDCommonParam::Instance();
fCal = AliTRDcalibDB::Instance();
fGeo = new AliTRDgeometry();
}
{
// loads the ADC data as obtained from the digitsManager for the specified MCM
- if (!CheckInitialized())
- Init(det, rob, mcm);
+ Init(det, rob, mcm);
if (!runloader) {
AliError("No Runloader given");
}
trdLoader->LoadDigits();
+ fDigitsManager = 0x0;
AliTRDdigitsManager *digMgr = new AliTRDdigitsManager();
digMgr->SetSDigits(0);
digMgr->CreateArrays();
}
}
}
- digMgr->RemoveDigits(det);
delete digMgr;
return kTRUE;
if (opt.Contains("T")) {
TLine *trklLines = new TLine[4];
- for (Int_t iTrkl = 0; iTrkl < 4; iTrkl++) {
- if (fMCMT[iTrkl] == 0x10001000)
- break;
+ for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) {
AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
Float_t offset = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19)) + 19 * pp->GetWidthIPad();
fADCF[iadc][it] = adc << fgkAddDigits;
}
-void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray)
+void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
{
// Set the ADC data from an AliTRDarrayADC
return;
}
+ fDigitsManager = digitsManager;
+
Int_t firstAdc = 0;
Int_t lastAdc = fNADC-1;
// Return column id of the pad for the given ADC channel
//
- if( !CheckInitialized() ) return -1;
+ if( !CheckInitialized() )
+ return -1;
Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
if (col < 0 || col >= fFeeParam->GetNcol())
if (nw < maxSize) {
buf[nw++] = x;
- //printf("\nMCM header: %X ",x);
+ //printf("\nMCM header: %X ",x);
}
else {
of++;
if (nw < maxSize) {
buf[nw++] = x;
- //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
+ //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
}
else {
of++;
// with -1 * number of overflowed words
//
- UInt_t x;
Int_t nw = 0; // Number of written words
Int_t of = 0; // Number of overflowed words
// 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 < maxSize)
+ buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
+ else
+ of++;
}
if( of != 0 ) return -of; else return nw;
// (lookup table accept (I2,I1,I0)=(111)
// or (110) or (101) or (100))
Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity)
- Int_t ep = AliTRDfeeParam::GetPFeffectPedestal();
+ Int_t ep = 0; // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); //??? really subtracted here
Int_t **adc = fADCF;
fZSM1Dim[iadc] &= fZSM[iadc][it];
}
}
-
}
void AliTRDmcmSim::DumpData( char *f, char *target )
//???
// TRAP parameters:
- const uint16_t lutPos[128] = { // move later to some other file
+ const UShort_t lutPos[128] = { // move later to some other file
0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15,
16, 16, 16, 17, 17, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26,
27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 27, 27, 27, 27, 26,
26, 26, 26, 25, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 20, 20, 19, 18, 18, 17, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 7};
//??? to be clarified:
- UInt_t adcMask = 0xfffff;
+ UInt_t adcMask = 0xffffffff;
UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
Short_t ypos, fromLeft, fromRight, found;
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*(adcLeft - adcRight) / adcCentral;
if (ypos < 0) ypos = -ypos;
// make the correction using the LUT
ypos = ypos + lutPos[ypos & 0x7F];
if (adcLeft > adcRight) ypos = -ypos;
- AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, -1);
+
+ // label calculation
+ Int_t mcLabel = -1;
+ if (fDigitsManager) {
+ Int_t label[9] = { 0 }; // up to 9 different labels possible
+ Int_t count[9] = { 0 };
+ Int_t maxIdx = -1;
+ Int_t maxCount = 0;
+ 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 (!fDigitsManager->UsesDictionaries() || fDigitsManager->GetDictionary(fDetector, iDict) == 0) {
+ AliError("Cannot get dictionary");
+ continue;
+ }
+ AliTRDarrayDictionary *dict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict);
+ if (dict->GetDim() == 0) {
+ AliError("Dictionary has dim. 0");
+ continue;
+ }
+ dict->Expand();
+ for (Int_t iPad = 0; iPad < 3; iPad++) {
+ if (padcol[iPad] < 0)
+ continue;
+ Int_t currLabel = dict->GetData(padrow, padcol[iPad], timebin); //fDigitsManager->GetTrack(iDict, padrow, padcol, timebin, fDetector);
+// printf("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin);
+ for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
+ if (currLabel == label[iLabel]) {
+ count[iLabel]++;
+ if (count[iLabel] > maxCount) {
+ maxCount = count[iLabel];
+ maxIdx = iLabel;
+ }
+ currLabel = 0;
+ break;
+ }
+ }
+ if (currLabel > 0) {
+ label[nLabels++] = currLabel;
+ }
+ }
+ }
+ if (maxIdx >= 0)
+ mcLabel = label[maxIdx];
+ }
+
+ // add the hit to the fitregister
+ AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, mcLabel);
}
}
}
for (i = ntracks; i < 4; i++) // CPUs without tracklets
fFitPtr[i] = 31; // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
// printf("found %i tracklet candidates\n", ntracks);
+// for (i = 0; i < 4; i++)
+// printf("fitPtr[%i]: %i\n", i, fFitPtr[i]);
}
void AliTRDmcmSim::FitTracklet()
Long64_t shift = ((Long64_t) 1 << 32);
UInt_t scaleY = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4)));
UInt_t scaleD = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4)));
+ Float_t scaleSlope = (256 / pp->GetWidthIPad()) * (1 << decPlaces);
+// printf("scaleSlope: %f \n", scaleSlope);
int padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
int yoffs = (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19) - fFeeParam->GetNcol()/2) << (8 + decPlaces);
int ndrift = 20; //??? value in simulation?
- int deflCorr = 0; // -370;
- int minslope = -10000; // no pt-cut so far
- int maxslope = 10000; // no pt-cut so far
+ Int_t deflCorr = -1 * (Int_t) (TMath::Tan(fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector))) * fGeo->CdrHght() * scaleSlope); // -370;
+ Int_t tiltCorr = -1 * (Int_t) (pp->GetRowPos(padrow) / fGeo->GetTime0(fDetector % 6) * fGeo->CdrHght() * scaleSlope *
+ TMath::Tan(pp->GetTiltingAngle() / 180. * TMath::Pi()));
+// printf("vdrift av.: %f\n", fCal->GetVdriftAverage(fDetector));
+// printf("chamber height: %f\n", fGeo->CdrHght());
+// printf("omega tau: %f\n", fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector)));
+// printf("deflection correction: %i\n", deflCorr);
+ Float_t ptcut = 2.3;
+ AliMagF* fld = (AliMagF *) TGeoGlobalMagField::Instance()->GetField();
+ Double_t bz = 0;
+ if (fld) {
+ bz = 0.1 * fld->SolenoidField(); // kGauss -> Tesla
+ }
+// printf("Bz: %f\n", bz);
+ Float_t x0 = fGeo->GetTime0(fDetector % 6);
+ Float_t y0 = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10));
+ Float_t alphaMax = TMath::ASin( (TMath::Sqrt(TMath::Power(x0/100., 2) + TMath::Power(y0/100., 2)) *
+ 0.3 * TMath::Abs(bz) ) / (2 * ptcut));
+// printf("alpha max: %f\n", alphaMax * 180/TMath::Pi());
+ Int_t minslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) + alphaMax) * scaleSlope);
+ Int_t maxslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) - alphaMax) * scaleSlope);
+// printf("min y-defl: %i\n", minslope);
+// printf("max y-defl: %i\n", maxslope);
// local variables for calculation
Long64_t mult, temp, denom; //???
sumXY = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
slope = nHits*sumXY - sumX * sumY;
+// printf("slope from fitreg: %i\n", slope);
offset = sumX2*sumY - sumX * sumXY;
temp = mult * slope;
slope = temp >> 32; // take the upper 32 bits
offset = temp >> 32; // take the upper 32 bits
offset = offset + yoffs + (18 << (8 + decPlaces));
- slope = slope * ndrift + deflCorr;
+// printf("slope: %i, slope * ndrift: %i, deflCorr: %i, tiltCorr: %i\n", slope, slope * ndrift, deflCorr, tiltCorr);
+ slope = slope * ndrift + deflCorr + tiltCorr;
offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
- if ((slope < minslope) || (slope > maxslope))
+// printf("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i ", fDetector, fRobPos, fMcmPos, slope, minslope, maxslope);
+ Bool_t rejected = kFALSE;
+ if (GetApplyCut() && ((slope < minslope) || (slope > maxslope)))
+ rejected = kTRUE;
+ if (rejected)
{
+// printf("rejected\n");
fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
}
else
{
+// printf("accepted\n");
temp = slope;
temp = temp * scaleD;
slope = (temp >> 32);
-
+// printf("slope after scaling: %i\n", slope);
+
temp = offset;
temp = temp * scaleY;
offset = (temp >> 32);
// rounding, like in the TRAP
slope = (slope + rndAdd) >> decPlaces;
+// printf("slope after shifting: %i\n", slope);
offset = (offset + rndAdd) >> decPlaces;
- if (slope > 0x3f || slope < -0x3f)
- AliWarning("Overflow in slope");
- slope = slope & 0x7F; // 7 bit
+ if (slope > 63) { // wrapping in TRAP!
+ AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
+ fMCMT[cpu] = 0x10001000;
+ continue;
+ }
+ else if (slope < -64) {
+ AliError(Form("Underflow in slope: %i, tracklet discarded!", slope));
+ fMCMT[cpu] = 0x10001000;
+ continue;
+ }
+ else {
+ slope = slope & 0x7F; // 7 bit
+ }
+// printf("slope after clipping: 0x%02x\n", slope);
- if (offset > 0xfff || offset < 0xfff)
+ if (offset > 0xfff || offset < -0xfff)
AliWarning("Overflow in offset");
offset = offset & 0x1FFF; // 13 bit
// assemble and store the tracklet word
fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
- new ((*fTrackletArray)[cpu]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
+
+ // calculate MC label
+ Int_t mcLabel = -1;
+ if (fDigitsManager) {
+ Int_t label[30] = {0}; // up to 30 different labels possible
+ Int_t count[30] = {0};
+ Int_t maxIdx = -1;
+ Int_t maxCount = 0;
+ Int_t nLabels = 0;
+ for (Int_t iHit = 0; iHit < fNHits; iHit++) {
+ if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) ||
+ (fHits[iHit].fChannel - fFitPtr[cpu] > 1))
+ continue;
+ Int_t currLabel = fHits[iHit].fLabel;
+ for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
+ if (currLabel == label[iLabel]) {
+ count[iLabel]++;
+ if (count[iLabel] > maxCount) {
+ maxCount = count[iLabel];
+ maxIdx = iLabel;
+ }
+ currLabel = 0;
+ break;
+ }
+ }
+ if (currLabel > 0) {
+ label[nLabels++] = currLabel;
+ }
+ }
+ if (maxIdx >= 0)
+ mcLabel = label[maxIdx];
+ }
+ new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
+ ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
}
}
}
fTrackletArray->Delete();
CalcFitreg();
+ if (fNHits == 0)
+ return;
TrackletSelection();
FitTracklet();
+ if (fTrackletArray->GetEntriesFast() == 0)
+ return;
AliRunLoader *rl = AliRunLoader::Instance();
AliDataLoader *dl = 0x0;
}
else {
TTree *trackletTree = dl->Tree();
- if (!trackletTree)
+ if (!trackletTree) {
dl->MakeTree();
- trackletTree = dl->Tree();
+ trackletTree = dl->Tree();
+ }
AliTRDtrackletMCM *trkl = 0x0;
TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
if (!trkbranch)
trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
-// trkbranch = trackletTree->Branch("mcmtrklbranch", &fTrackletArray, 32000, 2);
for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
return sum;
}
-void AliTRDmcmSim::Sort2(uint16_t idx1i, uint16_t idx2i, \
- uint16_t val1i, uint16_t val2i, \
- uint16_t *idx1o, uint16_t *idx2o, \
- uint16_t *val1o, uint16_t *val2o) const
+void AliTRDmcmSim::Sort2(UShort_t idx1i, UShort_t idx2i, \
+ UShort_t val1i, UShort_t val2i, \
+ UShort_t *idx1o, UShort_t *idx2o, \
+ UShort_t *val1o, UShort_t *val2o) const
{
// sorting for tracklet selection
}
}
-void AliTRDmcmSim::Sort3(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, \
- uint16_t val1i, uint16_t val2i, uint16_t val3i, \
- uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, \
- uint16_t *val1o, uint16_t *val2o, uint16_t *val3o)
+void AliTRDmcmSim::Sort3(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, \
+ UShort_t val1i, UShort_t val2i, UShort_t val3i, \
+ UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, \
+ UShort_t *val1o, UShort_t *val2o, UShort_t *val3o)
{
// sorting for tracklet selection
break;
default: // the rest should NEVER happen!
- printf("ERROR in Sort3!!!\n");
+ AliError("ERROR in Sort3!!!\n");
break;
}
// printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
}
-void AliTRDmcmSim::Sort6To4(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
- uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
- uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, uint16_t *idx4o, \
- uint16_t *val1o, uint16_t *val2o, uint16_t *val3o, uint16_t *val4o)
+void AliTRDmcmSim::Sort6To4(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, UShort_t idx4i, UShort_t idx5i, UShort_t idx6i, \
+ UShort_t val1i, UShort_t val2i, UShort_t val3i, UShort_t val4i, UShort_t val5i, UShort_t val6i, \
+ UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, UShort_t *idx4o, \
+ UShort_t *val1o, UShort_t *val2o, UShort_t *val3o, UShort_t *val4o)
{
// sorting for tracklet selection
- uint16_t idx21s, idx22s, idx23s, dummy;
- uint16_t val21s, val22s, val23s;
- uint16_t idx23as, idx23bs;
- uint16_t val23as, val23bs;
+ UShort_t idx21s, idx22s, idx23s, dummy;
+ UShort_t val21s, val22s, val23s;
+ UShort_t idx23as, idx23bs;
+ UShort_t val23as, val23bs;
Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
idx1o, &idx21s, &idx23as,
}
-void AliTRDmcmSim::Sort6To2Worst(uint16_t idx1i, uint16_t idx2i, uint16_t idx3i, uint16_t idx4i, uint16_t idx5i, uint16_t idx6i, \
- uint16_t val1i, uint16_t val2i, uint16_t val3i, uint16_t val4i, uint16_t val5i, uint16_t val6i, \
- uint16_t *idx5o, uint16_t *idx6o)
+void AliTRDmcmSim::Sort6To2Worst(UShort_t idx1i, UShort_t idx2i, UShort_t idx3i, UShort_t idx4i, UShort_t idx5i, UShort_t idx6i, \
+ UShort_t val1i, UShort_t val2i, UShort_t val3i, UShort_t val4i, UShort_t val5i, UShort_t val6i, \
+ UShort_t *idx5o, UShort_t *idx6o)
{
// sorting for tracklet selection
- uint16_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
- uint16_t val21s, val22s, val23s;
- uint16_t idx23as, idx23bs;
- uint16_t val23as, val23bs;
+ UShort_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
+ UShort_t val21s, val22s, val23s;
+ UShort_t idx23as, idx23bs;
+ UShort_t val23as, val23bs;
Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
&dummy1, &idx21s, &idx23as,
// printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
// idx21s, idx23as, idx22s, idx23bs, *idx5o, *idx6o);
}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff