bug fix from Julien - thanks
[u/mrichter/AliRoot.git] / 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 }