- 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
- }
-
- // 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 = fCal->GetOmegaTau(vdrift,magField); // 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: