- // 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];
- }
-
- // 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;
- }
- }
- }
-
- // 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);
- }
-
- 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
-
- 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;
-
- 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;
- }
-
-
- YAlu.AssignInt(1);
- Int_t wpad = YAlu.GetValue(); // 1 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;
- 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();
-
- 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]*wpad;
-
- 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]*wpad + wpad*wpad;
- 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 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)
- 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:
-
- 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();
- 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));
- 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);
- }
-
- offsetAlu = numAlu;
- mOffset[i] = offsetAlu.GetSignedValue();
-
-
- // slope:
-
- correctAlu.AssignDouble(mCorrectSlope);
- mCorrectSlope = correctAlu.GetValueWhole();
-
- mSlope[i] = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
-
- if (mSlope[i] < 0) {
- numAlu.AssignFormatted(-mSlope[i]);
- numAlu.SetSign(-1);
- }
- else {
- numAlu.AssignFormatted(mSlope[i]);
- numAlu.SetSign(1);
- }
-
- 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);
- }
- 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;
- }
-
- }
- }
-
- // 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");
- }
-
- // 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);
- }
- }
-
- // 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);
- }
- }
-
-
-
- //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;
-
-
-
-
-
- // output-part; creates an independent output .root folder if uncommented;
-
- // structure: in the current directory a root-file called "TRD_readout_tree.root" is stored with subdirectories SMxx/sx (supermodule, stack);
- // in each of these subdirectories 6 trees according to layers are saved, called lx;
- // whenever a mcm of that layer had a bit-word!=0, a branch containing an array with 4 (possibly some valued 0) elements is added;
- // branch-name: mcmxxxwd;
- // another branch contains the channel-number (mcmxxxch)
-
- /*
- AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
- AliLog::SetFileOutput("../log/tracklet.log");
-
-
- UInt_t* trackletWord;
- Int_t* adcChannel;
-
- Int_t u = 0;
-
- // testing for wordnr in order to speed up the simulation
-
- if (wordnr == 0 && fNextEvent == 0) {
- return;
- }
-
-
- Int_t mcmNr = fRobPos * (fGeo->MCMmax()) + fMcmPos;
-
- Char_t* SMName = new Char_t[4];
- Char_t* stackName = new Char_t[2];
- Char_t* layerName = new Char_t[2];
-
- Char_t* treeName = new Char_t[2];
- Char_t* treeTitle = new Char_t[8];
- Char_t* branchNameWd = new Char_t[8]; // mcmxxxwd bit word
- Char_t* branchNameCh = new Char_t[8]; // mcmxxxch channel
-
- Char_t* dirName = NULL;
- Char_t* treeFile = NULL;
- Char_t* evFile = NULL;
- Char_t* curDir = new Char_t[255];
-
- for (Int_t i = 0; i<255; i++) {
- curDir[i]='n';
- }
- sprintf(curDir,"%s",gSystem->BaseName(gSystem->WorkingDirectory()));
- Int_t rawEvent = 0;
- Int_t nrPos = 3;
-
-
- while(curDir[nrPos]!='n'){
-
-
- switch(curDir[nrPos]) {
- case '0':
- rawEvent = rawEvent*10 + 0;
- break;
- case '1':
- rawEvent = rawEvent*10 + 1;
- break;
- case '2':
- rawEvent = rawEvent*10 + 2;
- break;
- case '3':
- rawEvent = rawEvent*10 + 3;
- break;
- case '4':
- rawEvent = rawEvent*10 + 4;
- break;
- case '5':
- rawEvent = rawEvent*10 + 5;
- break;
- case '6':
- rawEvent = rawEvent*10 + 6;
- break;
- case '7':
- rawEvent = rawEvent*10 + 7;
- break;
- case '8':
- rawEvent = rawEvent*10 + 8;
- break;
- case '9':
- rawEvent = rawEvent*10 + 9;
- break;
-
- }
- nrPos = nrPos + 1;
- }
- delete [] curDir;
-
- //if (!gSystem->ChangeDirectory("../TRD_Tracklet")) {
- // gSystem->MakeDirectory("../TRD_Tracklet");
- //gSystem->ChangeDirectory("../TRD_Tracklet");
- //}
-
- gSystem->ChangeDirectory("..");
-
-
- TFile *f = new TFile("TRD_readout_tree.root","update");
- TTree *tree = NULL;
- TBranch *branch = NULL;
- TBranch *branchCh = NULL;
-
- Int_t iEventNr = 0;
- Int_t dignr = 10; // nr of digits of a integer
- Int_t space = 1; // additional char-space
-
-
- evFile = new Char_t[2+space];
- sprintf(evFile,"ev%d",iEventNr);
-
-
- while(f->cd(evFile)){
- iEventNr = iEventNr + 1;
- if (iEventNr/dignr > 0) {
- dignr = dignr * 10;
- space = space + 1;
- }
- delete [] evFile;
- evFile = NULL;
- evFile = new Char_t[2+space];
- sprintf(evFile,"ev%d",iEventNr);
- }
-
- if(iEventNr == rawEvent) { fNextEvent = 1; } // new event
-
- if (fNextEvent == 1) {
- fNextEvent = 0;
- // turn to head directory
- f->mkdir(evFile);
- f->cd(evFile);
- // create all subdirectories and trees in case of new event
-
-
- for (Int_t iSector = 0; iSector < 18; iSector++) {
-
- if (iSector < 10) {
- sprintf(SMName,"SM0%d",iSector);
- }
- else {
- sprintf(SMName,"SM%d",iSector);
- }
-
- for (Int_t iStack = 0; iStack < 5; iStack++) {
- sprintf(stackName,"s%d",iStack);
-
- f->cd(evFile);
- if (iStack == 0) {
- gDirectory->mkdir(SMName);
- }
- gDirectory->cd(SMName);
- gDirectory->mkdir(stackName);
- gDirectory->cd(stackName);
-
- for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
- sprintf(layerName,"l%d",iLayer);
- sprintf(treeName,"%s",layerName);
- sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
- tree = new TTree(treeName,treeTitle);
- tree->Write("",TObject::kOverwrite);
- delete tree;
- tree = NULL;
- }
- }
- }
-
-
- }
-
-
- else {
- iEventNr = iEventNr - 1;
- dignr = dignr/10;
- if (iEventNr/dignr == 0) space = space - 1;
- delete [] evFile;
- evFile = NULL;
- evFile = new Char_t[2+space];
- sprintf(evFile,"ev%d",iEventNr);
- }
-
- if (wordnr == 0) {
- delete [] SMName;
- delete [] stackName;
- delete [] layerName;
- delete [] treeName;
- delete [] treeTitle;
- delete [] branchNameWd;
- delete [] branchNameCh;
- delete [] evFile;
- f->Close();
- dirName = new Char_t[6+space];
- //sprintf(dirName,"../raw%d",iEventNr);
- sprintf(dirName,"raw%d",iEventNr);
- gSystem->ChangeDirectory(dirName);
- delete [] dirName;
- return;
- }
-
- dirName = new Char_t[6+space];
- //sprintf(dirName,"../raw%d",iEventNr);
- sprintf(dirName,"raw%d",iEventNr);
-
- f->cd(evFile);
-
- if (fSector < 10) {
- sprintf(SMName,"SM0%d",fSector);
- }
- else {
- sprintf(SMName,"SM%d",fSector);
- }
- sprintf(stackName,"s%d",fStack);
- sprintf(layerName,"l%d",fLayer);
- sprintf(treeName,"%s",layerName);
- sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
-
- treeFile = new Char_t[13+space];
- sprintf(treeFile,"%s/%s/%s/%s",evFile,SMName,stackName,treeName);
- delete [] evFile;
- evFile = NULL;
- gDirectory->cd(SMName);
- gDirectory->cd(stackName);
- tree = (TTree*)f->Get(treeFile);
- delete [] treeFile;
- treeFile = NULL;
-
-
- //make branch with number of words and fill
-
- if (mcmNr < 10) {
- sprintf(branchNameWd,"mcm00%dwd",mcmNr);
- sprintf(branchNameCh,"mcm00%dch",mcmNr);
- }
- else {
- if (mcmNr < 100) {
- sprintf(branchNameWd,"mcm0%dwd",mcmNr);
- sprintf(branchNameCh,"mcm0%dch",mcmNr);
- }
- else {
- sprintf(branchNameWd,"mcm%dwd",mcmNr);
- sprintf(branchNameCh,"mcm%dch",mcmNr);
- }
- }
-
-
-
- // fill the tracklet word; here wordnr > 0
-
- trackletWord = new UInt_t[fMaxTracklets];
- adcChannel = new Int_t[fMaxTracklets];
-
- 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
-
- //fMCMT[u] = bitWord[j];
- u = u + 1;
- }
- }
-
- branch = tree->GetBranch(branchNameWd);
- if(!branch) {
- //make branch and fill
- branch = tree->Branch(branchNameWd,trackletWord,"trackletWord[4]/i"); // 32 bit unsigned integer
- branch->Fill();
- }
-
- branchCh = tree->GetBranch(branchNameCh);
- if(!branchCh) {
- //make branch and fill
- branchCh = tree->Branch(branchNameCh,adcChannel,"adcChannel[4]/i"); // 32 bit unsigned integer
- branchCh->Fill();
- }
-
-
- tree->Write("",TObject::kOverwrite);
-
-
- delete [] SMName;
- delete [] stackName;
- delete [] layerName;
- delete [] treeName;
- delete [] treeTitle;
- delete [] branchNameWd;
- delete [] branchNameCh;
- delete [] trackletWord;
- delete [] adcChannel;
-
- f->Close();
- gSystem->ChangeDirectory(dirName);
- delete [] dirName;
- */
-
-
- // to be done:
- // error measure for quality of fit (not necessarily needed for the trigger)
- // cluster quality threshold (not yet set)
- // electron probability
-
-
-