+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(fgkDmemAddrDeflCorr, fDetector, fRobPos, fMcmPos);
+ Int_t ndrift = (Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrNdrift, fDetector, fRobPos, fMcmPos);
+
+ // local variables for calculation
+ Long64_t mult, temp, denom; //???
+ UInt_t q0, q1, 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();
+ }
+ 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 = ((Long64_t) nHits)*((Long64_t) sumX2) - ((Long64_t) sumX)*((Long64_t) sumX);
+
+ mult = mult / denom; // exactly like in the TRAP program
+ q0 = fit0->fQ0 + fit1->fQ0;
+ q1 = fit0->fQ1 + fit1->fQ1;
+ sumY = fit0->fSumY + fit1->fSumY + 256*fit1->fNhits;
+ sumXY = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
+ sumY2 = fit0->fSumY2 + fit1->fSumY2 + 512*fit1->fSumY + 256*256*fit1->fNhits;
+
+ slope = nHits*sumXY - sumX * sumY;
+ 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(fgkDmemAddrDeflCutStart + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos),
+ (Int_t) fTrapConfig->GetDmemUnsigned(fgkDmemAddrDeflCutStart + 1 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)));
+
+ AliDebug(5, Form("Fit sums: x = %i, X = %i, y = %i, Y = %i, Z = %i",
+ 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(fgkDmemAddrDeflCutStart + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos))) ||
+ (slope > ((Int_t) fTrapConfig->GetDmemUnsigned(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));
+
+ // store cluster information (if requested)
+ if (fgStoreClusters) {
+ Float_t *res = new Float_t[fNTimeBin];
+ Float_t *qtot = new Float_t[fNTimeBin];
+ for (Int_t iTimebin = 0; iTimebin < fNTimeBin; ++iTimebin) {
+ res[iTimebin] = 0;
+ qtot[iTimebin] = 0;
+ }
+ for (Int_t iHit = 0; iHit < fNHits; iHit++) {
+ Int_t timebin = fHits[iHit].fTimebin;
+
+ // check if hit contributes
+ if (fHits[iHit].fChannel == fFitPtr[cpu]) {
+ res[timebin] = fHits[iHit].fYpos - (fitSlope * timebin + fitOffset);
+ qtot[timebin] = fHits[iHit].fQtot;
+ }
+ else if (fHits[iHit].fChannel == fFitPtr[cpu] + 1) {
+ res[timebin] = fHits[iHit].fYpos + 256 - (fitSlope * timebin + fitOffset);
+ qtot[timebin] = fHits[iHit].fQtot;
+ }
+ }
+ ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetClusters(res, qtot, fNTimeBin);
+ delete [] res;
+ delete [] qtot;
+ }
+
+ if (fitError < 0)
+ AliError(Form("Strange fit error: %f from Sx: %i, Sy: %i, Sxy: %i, Sx2: %i, Sy2: %i, nHits: %i",
+ fitError, sumX, sumY, sumXY, sumX2, sumY2, nHits));
+ AliDebug(3, Form("fit slope: %f, offset: %f, error: %f",
+ fitSlope, fitOffset, TMath::Sqrt(TMath::Abs(fitError)/nHits)));
+ }
+ }
+ }
+}