EMCAL/SMcalib/WriteNewBias.C

   1 /*
   2 */
   3 const int debug = 0; // 0=quiet; 1=print some parameter info
   4
   5 const int kMaxHV = 395; // default max voltage limit; could possibly be relaxed in future
   6 const int kBreakDownMargin = 5; // at least 5 V away from breakDown voltage..
   7
   8 static const int fgkEmCalRows = 24; // number of rows per module for EMCAL
   9 static const int fgkEmCalCols = 48; // number of columns per module for EMCAL
  10
  11 Float_t previousVoltage[fgkEmCalCols][fgkEmCalRows];
  12 Float_t gainFactor[fgkEmCalCols][fgkEmCalRows];
  13 Float_t biasVoltage[fgkEmCalCols][fgkEmCalRows];
  14
  15 //____________________________________________________________________
  16 void WriteNewBias(const char * inputDBName, const char * inputMapName,
  17                   const char * previousVoltageFileName, const char * gainFactorFileName,
  18                   const char * outputFileName)
  19 {
  20
  21   ofstream outputFile(outputFileName);
  22   ReadFiles(previousVoltageFileName, gainFactorFileName);
  23
  24   SetBiasVoltage(inputDBName, inputMapName);
  25
  26   for (int icol=0; icol<fgkEmCalCols; icol++) {
  27     for (int irow=0; irow<fgkEmCalRows; irow++) {
  28       outputFile << icol << " " << irow << " " << biasVoltage[icol][irow] << endl;
  29     }
  30   }
  31
  32   outputFile.close();
  33 }
  34
  35 //____________________________________________________________________
  36 void ReadFiles(const char * previousVoltageFileName,
  37                const char * gainFactorFileName)
  38 {
  39   ifstream previousVoltageFile(previousVoltageFileName);
  40   ifstream gainFactorFile(gainFactorFileName);
  41
  42   int icolp, irowp;
  43   int icolg, irowg;
  44   Float_t gFactor;
  45   Float_t pVoltage;
  46
  47   for (int icol=0; icol<fgkEmCalCols; icol++) {
  48     for (int irow=0; irow<fgkEmCalRows; irow++) {
  49       previousVoltageFile >> icolp >> irowp >> pVoltage;
  50       gainFactorFile >> icolg >> irowg >> gFactor;
  51
  52       previousVoltage[icolp][irowp] = pVoltage;
  53       gainFactor[icolg][irowg] = gFactor;
  54     }
  55   }
  56
  57   previousVoltageFile.close();
  58   gainFactorFile.close();
  59
  60   return;
  61 }
  62
  63 //____________________________________________________________________
  64 void SetBiasVoltage(const char * inputDBName, const char * inputMapName)
  65 {
  66   gSystem->Load("AliEMCALCalibAPD_cxx");
  67   AliEMCALCalibAPD *calibAPD = new AliEMCALCalibAPD();
  68
  69   calibAPD->ReadCalibAPDInfo(10000, inputDBName);
  70   int fNCalibAPD = calibAPD->GetNCalibAPD();
  71   AliEMCALCalibAPD::AliEMCALCalibAPDData * fCalib = calibAPD->GetCalibAPDData();
  72
  73
  74   gSystem->Load("AliEMCALMapAPD_cxx");
  75   AliEMCALMapAPD *mapAPD = new AliEMCALMapAPD();
  76
  77   int nSM = 1;
  78   mapAPD->ReadMapAPDInfo(nSM, inputMapName);
  79   AliEMCALMapAPD::AliEMCALSuperModuleMapAPD * fMap = mapAPD->GetSuperModuleData();
  80
  81   int nFound = 0;
  82   for (int icol=0; icol<fgkEmCalCols; icol++) {
  83     for (int irow=0; irow<fgkEmCalRows; irow++) {
  84
  85       int apdMap = fMap[0].fAPDNum[icol][irow]; // 0 = nSM - 1
  86       int i = 0;
  87       int apdCalib = -1;
  88       while (i<fNCalibAPD && apdMap!=apdCalib) {
  89         apdCalib = fCalib[i].fAPDNum;
  90         i++;
  91       }
  92
  93       if (apdCalib == apdMap) { // found!
  94         i--; // go back to what we dound
  95
  96         // estimate what the new/target HV should be
  97         biasVoltage[icol][irow] = CalculateTargetHV(previousVoltage[icol][irow],
  98                                                     gainFactor[icol][irow],
  99                                                     fCalib[i].fPar[0],
 100                                                     fCalib[i].fPar[1],
 101                                                     fCalib[i].fPar[2],
 102                                                     fCalib[i].fBreakDown);
 103         nFound++;
 104       }
 105       else { // no calib info, just use old settings
 106         biasVoltage[icol][irow] = previousVoltage[icol][irow];
 107       }
 108
 109     }
 110   }
 111
 112   cout << " found " << nFound << " matches " << endl;
 113   return;
 114 }
 115
 116 //____________________________________________________________________
 117 Float_t CalculateTargetHV(Float_t initialHV, Float_t gainChange,
 118                           Float_t par0, Float_t par1, Float_t par2, Int_t breakDown)
 119 {
 120   if (debug) { printf("parameters p0:%g p1:%g p2:%g\n", par0, par1, par2); }
 121
 122   // figure out what new HV should be,
 123   // if we want to adjust the gain by some factor
 124   Float_t initialGain = par0 + par1 * exp(par2*initialHV);
 125   Float_t newGain = initialGain * gainChange; // = par0 + par1 * exp(par2*newHV);
 126
 127   if (debug) { printf("initialGain:%g newGain:%g\n", initialGain, newGain); }
 128
 129   Float_t newHV = -1;
 130   if ( par1>0 && par2>0 ) {
 131     newHV = log ( (newGain - par0)/par1 ) / par2;
 132   }
 133
 134   // check results before returning..
 135   if (newHV < 0) {
 136    // conversion failed:  let's just keep the old custom value then
 137     newHV = initialHV;
 138   }
 139   if (newHV>kMaxHV) {
 140     // we reached a too high voltage: let's keep the max then
 141    newHV = kMaxHV;
 142   }
 143
 144   // in case we increase the kMaxHV limit some time in the future we could
 145   // also enter get close to the Hamamatsu breakdown limit - let's avoid that
 146   if (newHV>breakDown) {
 147     newHV = breakDown - kBreakDownMargin;
 148   }
 149
 150   return newHV;
 151 }