TRD/AliTRDtrapConfigHandler.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 ////////////////////////////////////////////////////////////////////////////
  18 //                                                                        //
  19 //  MCM configuraton handler                                              //
  20 //                                                                        //
  21 //  Author: U. Westerhoff                                                 //
  22 //                                                                        //
  23 ////////////////////////////////////////////////////////////////////////////
  24
  25
  26
  27 #include "AliTRDtrapConfigHandler.h"
  28
  29 #include <iostream>
  30 #include <iomanip>
  31
  32 #include "AliTRDpadPlane.h"
  33 #include "AliTRDtrapConfig.h"
  34 #include "AliTRDcalibDB.h"
  35 #include "AliTRDCommonParam.h"
  36 #include "AliLog.h"
  37
  38 #include "AliTRDarrayDictionary.h"
  39
  40 #include "AliMagF.h"
  41 #include "TGeoGlobalMagField.h"
  42 #include "AliMagWrapCheb.h"
  43 #include "AliTRDgeometry.h"
  44
  45 #include "TMath.h"
  46 #include "TGeoMatrix.h"
  47
  48 using namespace std;
  49
  50 ClassImp(AliTRDtrapConfigHandler)
  51
  52
  53 AliTRDtrapConfigHandler::AliTRDtrapConfigHandler() :
  54      fGeo(NULL)
  55      , fDet(0)
  56      , fBField(0)
  57      , fOmegaTau(1)
  58      , fPtMin(0)
  59      , fNTimebins(0)
  60      , fScaleQ0(1)
  61      , fScaleQ1(1)
  62      , fPidTracklengthCorr(kFALSE)
  63      , fTiltCorr(kTRUE)
  64 {
  65    fGeo = new AliTRDgeometry;
  66 }
  67
  68
  69 AliTRDtrapConfigHandler::~AliTRDtrapConfigHandler()
  70 {
  71    delete fGeo;
  72 }
  73
  74 void AliTRDtrapConfigHandler::ResetMCMs()
  75 {
  76    //
  77    // Reset all MCM registers and DMEM
  78    //
  79
  80    AliTRDtrapConfig *cfg = AliTRDtrapConfig::Instance();
  81    cfg->ResetRegs();
  82    cfg->ResetDmem();
  83 }
  84
  85
  86 Int_t AliTRDtrapConfigHandler::LoadConfig(TString filename, Int_t det)
  87 {
  88    //
  89   // load a TRAP configuration from a file
  90    // The file format is the format created by the standalone
  91    // command coder scc / show_cfdat but without the running number
  92    // scc /show_cfdat adds as the first column
  93    // which are two tools to inspect/export configurations from wingDB
  94    //
  95
  96
  97    fDet = det;
  98    Int_t ignoredLines=0;
  99    Int_t ignoredCmds=0;
 100    Int_t readLines=0;
 101
 102
 103    AliTRDtrapConfig *cfg = AliTRDtrapConfig::Instance();
 104
 105    AliDebug(5, Form("Processing file %s", filename.Data()));
 106    std::ifstream infile;
 107    infile.open(filename.Data(), std::ifstream::in);
 108    if (!infile.is_open()) {
 109     AliError(Form("Can not open MCM configuration file %s", filename.Data()));
 110     return kFALSE;
 111    }
 112
 113    UInt_t cmd;
 114    Int_t extali, addr, data;
 115
 116    while(infile.good()) {
 117       cmd=999;
 118       extali=-1;
 119       addr=-1;
 120       data=-1;
 121       infile >> std::skipws >> cmd >> addr >> data >> extali;
 122       //      std::cout << "no: " << no << ", cmd " << cmd << ", extali " << extali << ", addr " << addr << ", data " << data <<  endl;
 123
 124       if(cmd!=999 && extali!=-1 && addr != -1 && data!= -1 && extali!=-1) {
 125          if(cmd==fgkScsnCmdWrite)
 126             cfg->AddValues(det, cmd, extali, addr, data);
 127          else if(cmd == fgkScsnLTUparam)
 128             ProcessLTUparam(extali, addr, data);
 129          else
 130             ignoredCmds++;
 131
 132          readLines++;
 133       }
 134       else if(!infile.eof() && !infile.good()) {
 135          infile.clear();
 136          infile.ignore(256, '\n');
 137          ignoredLines++;
 138       }
 139
 140       if(!infile.eof())
 141          infile.clear();
 142    }
 143
 144    infile.close();
 145
 146    AliDebug(5, Form("Ignored lines: %i, ignored cmds: %i", ignoredLines, ignoredCmds));
 147
 148
 149    if(ignoredLines>readLines)
 150       AliError(Form("More than 50 %% of the input file could not be processed. Perhaps you should check the input file %s", filename.Data()));
 151
 152
 153    return kTRUE;
 154 }
 155
 156
 157
 158 void AliTRDtrapConfigHandler::ProcessLTUparam(Int_t dest, Int_t addr, UInt_t data)
 159 {
 160    //
 161    // Process the LTU parameters and stores them in internal class variables
 162    // or transfer the stored values to AliTRDtrapConfig, depending on the dest parameter
 163    //
 164
 165    switch (dest) {
 166
 167    case 0: // set the parameters in AliTRDtrapConfig
 168       ConfigureDyCorr();
 169       ConfigureDRange(); // deflection range
 170       ConfigureNTimebins();  // timebins in the drift region
 171       ConfigurePIDcorr();  // scaling parameters for the PID
 172       break;
 173
 174    case 1: // set variables
 175       switch (addr) {
 176
 177       case 0: fPtMin =  float(data) / 1000.; break; // pt_min in GeV/c (*1000)
 178       case 1: fBField = float(data) / 1000. ; break; // B in T (*1000)
 179       case 2: fOmegaTau = float(data) / 1.0E6  ; break; // omega*tau
 180       case 3: fNTimebins = data; break;
 181         // ntimbins: drift time (for 3 cm) in timebins (5 add. bin. digits)
 182       case 4: fScaleQ0 = data; break;
 183       case 5: fScaleQ1 = data; break;
 184       case 6: fPidTracklengthCorr = (bool) data; break;
 185       case 7: fTiltCorr = (bool) data; break;
 186       }
 187       break;
 188
 189    default:
 190       AliError(Form("dest %i not implemented", dest));
 191    }
 192
 193 }
 194
 195
 196 void AliTRDtrapConfigHandler::ConfigureNTimebins()
 197 {
 198    //
 199    // Set timebins in the drift region
 200    //
 201    AliTRDtrapConfig::Instance()->AddValues(fDet, fgkScsnCmdWrite, 127, AliTRDtrapConfig::fgkDmemAddrNdrift, fNTimebins);
 202 }
 203
 204
 205
 206 void AliTRDtrapConfigHandler::ConfigureDyCorr()
 207 {
 208    //
 209    //  Deflection length correction
 210    //  due to Lorentz angle and tilted pad correction
 211    //  This correction is in units of padwidth / (256*32)
 212    //
 213
 214    Int_t nRobs;
 215    if(fGeo->GetStack(fDet) == 2)
 216       nRobs=6;
 217    else
 218       nRobs=8;
 219
 220    Double_t dyLorentz = - fOmegaTau * fGeo->CdrHght();    // CdrHght: Height of the drift region
 221    Double_t globalPos[3];
 222    Double_t tiltingAngle = fGeo->GetPadPlane(fDet)->GetTiltingAngle();
 223    Double_t scalePad = 256. * 32.;
 224    Int_t dyCorrInt=0;
 225
 226     for (Int_t r = 0; r < nRobs; r++) {
 227         for (Int_t m = 0; m < 16; m++) {
 228            if(GetPadPosNonRot(r, m, 9, globalPos)==0) {
 229               double dyTilt = ( fGeo->CdrHght() *  TMath::Tan(tiltingAngle * TMath::DegToRad()) * globalPos[2]/globalPos[0]);
 230
 231               double dyCorr;
 232               if(fTiltCorr==kTRUE)
 233                  dyCorr = dyLorentz + dyTilt;
 234               else
 235                  dyCorr = dyTilt;
 236
 237
 238               dyCorrInt =  TMath::Nint(dyCorr * scalePad /  fGeo->GetPadPlane(fDet)->GetWidthIPad());  // The correction is in units of 1/256 of the
 239                                                                                                        // pad width, including 5 binary digits
 240            }
 241            else {
 242               AliError("No transformation matrix available");
 243               dyCorrInt=0;
 244            }
 245            Int_t dest =  1<<10 | r<<7 | m;
 246            AliTRDtrapConfig::Instance()->AddValues(fDet, fgkScsnCmdWrite, dest, AliTRDtrapConfig::fgkDmemAddrDeflCorr, dyCorrInt);
 247         }
 248     }
 249 }
 250
 251
 252
 253
 254
 255 void AliTRDtrapConfigHandler::ConfigureDRange()
 256 {
 257    //
 258    // deflection range LUT
 259    // range calculated according to B-field (in T) and pt_min (in GeV/c)
 260    // if pt_min < 0.1 GeV/c the maximal allowed range for the tracklet
 261    // deflection (-64..63) is used
 262    //
 263
 264    static const int x=0;
 265    static const int y=1;
 266    static const Double_t dyBin = 140e-6;
 267
 268    Int_t dyMinInt=fgkDyMinCut;
 269    Int_t dyMaxInt=fgkDyMaxCut;
 270    Double_t mcmPos[3];
 271
 272    Int_t nRobs=-1;
 273    if(fGeo->GetStack(fDet) == 2)
 274       nRobs=6;
 275    else
 276       nRobs=8;
 277
 278    for (Int_t r = 0; r < nRobs; r++) {
 279       for (Int_t m = 0; m < 16; m++) {
 280          for (Int_t c = 0; c < 18; c++) {
 281
 282             if(fPtMin<0.1) {
 283                dyMinInt=fgkDyMinCut;
 284                dyMaxInt=fgkDyMaxCut;
 285             }
 286             else {
 287                if(GetPadPosNonRot(r, m, c, mcmPos)==0) {
 288                   Double_t radius = fPtMin/(0.3*fBField);
 289
 290                   double vertexPos[2] = {0,0};
 291
 292                   Double_t distanceX = (vertexPos[x]-mcmPos[x]) / 100.; // cm -> m
 293                   Double_t distanceY = (vertexPos[y]-mcmPos[y]) / 100.; // cm -> m
 294
 295                   Double_t maxDeflTemp = (TMath::Sqrt( Square(distanceX) + Square(distanceY)) / 2) / radius;
 296                   Double_t localPhi = TMath::ATan2(distanceY, distanceX);
 297
 298                   Double_t maxDeflAngle=0;
 299                   if(maxDeflTemp < 1. ) {
 300                      maxDeflAngle = TMath::ASin(maxDeflTemp);
 301                      Double_t dyMin = fGeo->CdrHght()/100. * TMath::Tan(localPhi - maxDeflAngle);  // CdrHght: Height of the drift region in cm
 302                      Double_t dyMax = fGeo->CdrHght()/100. * TMath::Tan(localPhi + maxDeflAngle);
 303
 304                      dyMinInt = Int_t (dyMin / dyBin);
 305                      dyMaxInt = Int_t (dyMax / dyBin);
 306
 307                      if(dyMinInt < fgkDyMinCut)
 308                         dyMinInt = fgkDyMinCut;
 309                      if(dyMaxInt > fgkDyMaxCut)
 310                         dyMaxInt = fgkDyMaxCut;
 311                   }
 312                   else {
 313                      dyMinInt = fgkDyMinCut;
 314                      dyMaxInt = fgkDyMaxCut;
 315                   }
 316
 317 //                cout << "maxdefl: " << maxDeflAngle << ", localPhi " << localPhi << endl;
 318 //                cout << "r " << r << ", m" << m << ", c " << c << ", min angle: " << localPhi-maxDeflAngle << ", max: " << localPhi+maxDeflAngle
 319 //                     << ", min int: " << dyMinInt << ", max int: " << dyMaxInt << endl;
 320                }
 321                else {
 322                   AliError("No geometry model loaded\n");
 323                }
 324             }
 325
 326             Int_t dest =  1<<10 | r<<7 | m;
 327             Int_t lutAddr = AliTRDtrapConfig::fgkDmemAddrDeflCutStart + 2*c;
 328             AliTRDtrapConfig::Instance()->AddValues(fDet, fgkScsnCmdWrite, dest, lutAddr+0, dyMinInt);
 329             AliTRDtrapConfig::Instance()->AddValues(fDet, fgkScsnCmdWrite, dest, lutAddr+1, dyMaxInt);
 330          }
 331       }
 332    }
 333 }
 334
 335 void AliTRDtrapConfigHandler::PrintGeoTest()
 336 {
 337    //
 338    // Prints some information about the geometry. Only for debugging
 339    //
 340
 341    Double_t mcmPos[3];
 342    int sm=0;
 343    //   for(int sm=0; sm<6; sm++) {
 344    for(int stack=0; stack<5; stack++) {
 345       for(int layer=0; layer<6; layer++) {
 346
 347          fDet = sm*30+stack*6+layer;
 348          for (Int_t r = 0; r < 6; r++) {
 349             for (Int_t m = 0; m < 16; m++) {
 350                for (Int_t c = 7; c < 8; c++) {
 351                   GetPadPosNonRot(r, m, c, mcmPos);
 352                   cout << stack << ";" << layer << ";" << r << ";" << m
 353                        << ";" << mcmPos[0] << ";" << mcmPos[1] << ";" << mcmPos[2] << endl;
 354                }
 355             }
 356          }
 357       }
 358    }
 359    // }
 360 }
 361
 362
 363
 364 void AliTRDtrapConfigHandler::ConfigurePIDcorr()
 365 {
 366    //
 367    // Calculate the MCM individual correction factors for the PID
 368    // and transfer them to AliTRDtrapConfig
 369    //
 370
 371    static const Int_t addrLUTcor0 = AliTRDtrapConfig::fgkDmemAddrLUTcor0;
 372    static const Int_t addrLUTcor1 = AliTRDtrapConfig::fgkDmemAddrLUTcor1;
 373
 374    UInt_t cor0;
 375    UInt_t cor1;
 376
 377    Double_t globalPos[3];
 378
 379    Int_t nRobs=-1;
 380    if(fGeo->GetStack(fDet) == 2)
 381       nRobs=6;
 382    else
 383       nRobs=8;
 384
 385    for (Int_t r=0; r<nRobs; r++) {
 386       for(Int_t m=0; m<16; m++) {
 387
 388          if(GetPadPosNonRot(r, m, 9, globalPos)==0) {
 389             Double_t elongation = TMath::Abs(TMath::Sqrt(globalPos[0]*globalPos[0] + globalPos[1]*globalPos[1] + globalPos[2]*globalPos[2]) / globalPos[0]);
 390
 391             if(fPidTracklengthCorr==kFALSE) {
 392                cor0 = fScaleQ0;
 393                cor1 = fScaleQ1;
 394             }
 395             else {
 396                cor0 = Int_t ((1.0*fScaleQ0* (1/elongation) ));
 397                cor1 = Int_t ((1.0*fScaleQ1* (1/elongation) ));
 398             }
 399
 400             Int_t dest =  1<<10 | r<<7 | m;
 401             AliTRDtrapConfig::Instance()->AddValues(fDet, fgkScsnCmdWrite, dest, addrLUTcor0, cor0);
 402             AliTRDtrapConfig::Instance()->AddValues(fDet, fgkScsnCmdWrite, dest, addrLUTcor1, cor1);
 403          }
 404          else {
 405             AliError("No transformation matrix available");
 406          }
 407       }
 408    }
 409 }
 410
 411
 412
 413 Int_t AliTRDtrapConfigHandler::GetPadPosNonRot(Int_t rob, Int_t mcm, Int_t channel, Double_t trackCoor[3])
 414 {
 415    //
 416    // Calcutate the gobal coordinates for an mcm channel in the supermodule at position -0.5
 417    //
 418
 419    Int_t stack = fGeo->GetStack(fDet);
 420    Int_t layer = fGeo->GetLayer(fDet);
 421
 422    AliTRDpadPlane *plane = fGeo->GetPadPlane(layer, stack);
 423    if(plane==NULL) {
 424       AliError(Form("stack %i, layer %i, det %i", stack, layer, fDet));
 425       return 1;
 426    }
 427
 428    Double_t locYZ[2];
 429    Double_t loc[3];
 430
 431    Double_t radialPos = fGeo->AnodePos()-0.83; //
 432    //   Double_t radialPos = 300.65 + 12.60 * layer; // cm // taken from libTRD/geometry/geometry.cc, probably not the final value
 433
 434    GetLocalPadPos(plane, rob, mcm, channel, locYZ);
 435
 436    loc[0] = radialPos;
 437    loc[1] = locYZ[0];
 438    loc[2] = locYZ[1];
 439
 440    //transform from loc[3] - coordinates of the pad
 441    // Go to tracking coordinates
 442
 443    TGeoHMatrix *fMatrix  = fGeo->GetClusterMatrix(fDet);
 444    if(fMatrix==NULL) {
 445       AliError(Form("stack %i, layer %i, det %i", stack, layer, fDet));
 446       return 2;
 447    }
 448    fMatrix->LocalToMaster(loc, trackCoor);
 449    return 0;
 450 }
 451
 452
 453 void AliTRDtrapConfigHandler::GetLocalPadPos(AliTRDpadPlane *plane, Int_t rob, Int_t mcm, Int_t channel, Double_t result[2])
 454 {
 455    //
 456    // calculate the local coordinates for an mcm channel
 457    //
 458
 459    Double_t localY, localZ;
 460
 461     Int_t padCol;
 462     if(rob%2 == 0) //side a
 463        padCol  = (mcm % fgkMCMperROBCol) * fgkPadsPerMCM + channel;
 464     else
 465        padCol  = (mcm % fgkMCMperROBCol) * fgkPadsPerMCM + (plane->GetNcols()/2) + channel;
 466
 467     Int_t padRow = ((Int_t) floor(rob/2.0)) * fgkMCMperROBRow + ((Int_t) floor(mcm/4));
 468
 469     if(padCol<0 || padCol>= plane->GetNcols())
 470        AliError(Form("Invalid pad col: %i\n", padCol));
 471
 472     if(padRow<0 || padRow>= plane->GetNrows())
 473        AliError(Form("Invalid pad row: %i\n", padRow));
 474
 475     if(padCol+1 == plane->GetNcols()) // last pad
 476        localY = plane->GetColPos(padCol) + (plane->GetColEnd()-plane->GetColPos(padCol))/2;
 477     else
 478        localY = plane->GetColPos(padCol) + (plane->GetColPos(padCol+1)-plane->GetColPos(padCol))/2;
 479
 480     if(padRow+1 == plane->GetNrows())
 481        localZ = plane->GetRowPosROC(padRow) + (plane->GetRowEndROC() - plane->GetRowPosROC(padRow))/2;
 482     else
 483        localZ = plane->GetRowPosROC(padRow) + (plane->GetRowPosROC(padRow+1) - plane->GetRowPosROC(padRow))/2;
 484
 485     //    std::cout << "pad col " << padCol << ", pad row " << padRow << std::endl;
 486     //    std::cout << "pos y (col) " << localY << ", pos z (row) " << localZ << std::endl;
 487
 488     result[0]=localY;
 489     result[1]=localZ;
 490 }
 491
 492
 493 Double_t AliTRDtrapConfigHandler::Square(Double_t val)
 494 {
 495    //
 496    // calculate the square of the argument
 497    //
 498
 499    return val*val;
 500 }