Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:52:26

0001 // Macro to fit Peters cal files.
0002 // This uses A3 + A3*tanh(A0*x-A1) from Urs.
0003 // This fits ADC vs VCAL (x-axis=VCAL) - opposite to Peter
0004 // I have tried VCAL vs ADC with ATanH() but it crashes.
0005 // Compress ADC range to 0-255 range to simulate FED data.
0006 // Usage:
0007 // .L PHCalibration.C
0008 //          dir  roc pixel linear
0009 // FitCurve("mtb",0,20,20,false)  - one pixel
0010 //              dir   linear
0011 // FitAllCurves("mtb",false)      - all pixels
0012 // Plots()                 - show some plost
0013 
0014 #include <iostream>
0015 #include <TSystem.h>
0016 #include "TF1.h"
0017 #include "TGraph.h"
0018 #include "TH1D.h"
0019 #include "TFile.h"
0020 
0021 // Conversion between high and low range 
0022 static const double rangeConversion = 7.;  // 
0023 // Number of VCAL point per range
0024 static const int vcalSteps = 5;
0025 double vcal[2*vcalSteps], vcalLow[2*vcalSteps];
0026 double x[2*vcalSteps], y[2*vcalSteps], xErr[2*vcalSteps], yErr[2*vcalSteps];
0027 int n;
0028 
0029 // Number of parameters to fit
0030 static const int nFitParams = 4;
0031 
0032 const bool HISTO = true;
0033 
0034 TH1D *histoFit[nFitParams], *histoChi;
0035 TH1D *histoFit0[16],*histoFit1[16],*histoFit2[16],*histoFit3[16];
0036 
0037 TGraph *graph;
0038 TF1 *phFit, *phLinear;
0039 
0040 // Calibration data points per ROC
0041 TH2D *h2d0=0,*h2d1,*h2d2,*h2d3,*h2d4,*h2d5,*h2d6,*h2d7,*h2d8,*h2d9,
0042   *h2d10,*h2d11,*h2d12,*h2d13,*h2d14,*h2d15;
0043 // Fit distribution 
0044 TH2D *h2dp0, *h2dp1, *h2dp2, *h2dp3, *h2dchis;
0045 
0046 //const double xCut = TMath::Pi()/2. - 0.0005;
0047 //const double tanXCut = TMath::Tan(xCut);
0048 
0049 //===========================================================================
0050 // Function to fit
0051 Double_t Fitfcn( Double_t *x, Double_t *par) {
0052 
0053   //if (par[0]*x[0] - par[4] > xCut) return tanXCut + (x[0] - (xCut + par[4])/par[0])* 1e8;
0054   //Double_t y = TMath::Tan(par[0]*x[0] - par[4]) + par[1]*x[0]*x[0]*x[0] + par[5]*x[0]*x[0] + par[2]*x[0] + par[3];
0055 
0056   Double_t y = par[3] + par[2] * TMath::TanH(par[0]*x[0]-par[1]);
0057   return y;
0058 }
0059 //===========================================================================
0060 // Function to fit
0061 Double_t FitLinear( Double_t *x, Double_t *par) {
0062 
0063   //if (par[0]*x[0] - par[4] > xCut) return tanXCut + (x[0] - (xCut + par[4])/par[0])* 1e8;
0064   //Double_t y = TMath::Tan(par[0]*x[0] - par[4]) + par[1]*x[0]*x[0]*x[0] + par[5]*x[0]*x[0] + par[2]*x[0] + par[3];
0065 
0066   Double_t y = par[0] + par[1] * x[0];
0067   return y;
0068 }
0069 //==========================================================================
0070 //String with the function 
0071 const char *Fitfcn() {
0072   //    return "TMath::Exp(par[1]*x[0] - par[0]) + par[2]*x[0]*x[0]*x[0] + par[3]*x[0]*x[0] + par[4]*x[0] + par[5]";
0073  // return "TMath::Tan(par[0]*x[0] - par[4]) + par[1]*x[0]*x[0]*x[0] + par[5]*x[0]*x[0] + par[2]*x[0] + par[3]";
0074   return "par[3] + par[2] * TMath::TanH(par[0]*x[0] - par[1])";
0075 }
0076 //String with the function 
0077 const char *FitLinear() {
0078   //    return "TMath::Exp(par[1]*x[0] - par[0]) + par[2]*x[0]*x[0]*x[0] + par[3]*x[0]*x[0] + par[4]*x[0] + par[5]";
0079  // return "TMath::Tan(par[0]*x[0] - par[4]) + par[1]*x[0]*x[0]*x[0] + par[5]*x[0]*x[0] + par[2]*x[0] + par[3]";
0080   return "par[0] + par[1] * x[0]";
0081 }
0082 //===========================================================================
0083 // Transform input ADC data to the 0-255 range
0084 int transformADC(int adc) {
0085   const int offset = 330.;
0086   const int gain = 6;
0087   int temp = (adc + offset)/gain;
0088   if(temp<0 || temp>255) cout<<" adc wrong "<<temp<<" " <<adc<<endl;
0089   return temp;
0090 }
0091 //=========================================================================== 
0092 // Init: test histos, fit function, errors
0093 void Initialize(bool linear=false) {
0094   gROOT->SetStyle("Plain");
0095   gStyle->SetTitleBorderSize(0);
0096   gStyle->SetPalette(1,0);
0097   
0098   float p0max = 0.01, p1max=2.,p2max=1000.,p3max=1000.;
0099   float p0min = 0.;
0100   if(linear) {p0max = 100.; p0min=-100.;}
0101  
0102   histoFit[0] = new TH1D("histoFit0", "histoFit0", 100, p0min, p0max);
0103   histoFit[1] = new TH1D("histoFit1", "histoFit1", 200, 0.0, p1max);
0104   histoFit[2] = new TH1D("histoFit2", "histoFit2", 100, 0.0, p2max);
0105   histoFit[3] = new TH1D("histoFit3", "histoFit3", 100, 0.0, p3max);
0106   histoChi = new TH1D("histoChi", "histoChi", 1000, 0., 10.);
0107 
0108   char hiname[20];
0109   
0110   for(int i=0;i<16;i++) {
0111     sprintf(hiname, "histoFit0_%i", i);
0112     histoFit0[i] = new TH1D(hiname,hiname, 100, p0min, p0max);
0113     sprintf(hiname, "histoFit1_%i", i);
0114     histoFit1[i] = new TH1D(hiname,hiname, 200, 0.0, p1max);
0115     sprintf(hiname, "histoFit2_%i", i);
0116     histoFit2[i] = new TH1D(hiname,hiname, 100, 0.0, p2max);
0117     sprintf(hiname, "histoFit3_%i", i);
0118     histoFit3[i] = new TH1D(hiname,hiname, 100, 0.0, p3max);
0119   }  
0120 
0121   if(HISTO) {
0122     const float xdmin=0., xdmax=260.;
0123     const int xdbins=130;
0124     if(h2d0!=0) delete h2d0;
0125     h2d0 = new TH2D("h2d0", "h2d0", xdbins,xdmin,xdmax, 150, 0., 1500.);
0126     if(h2d1!=0) delete h2d1;
0127     h2d1 = new TH2D("h2d1", "h2d1", xdbins,xdmin,xdmax, 150, 0., 1500.);
0128     if(h2d2!=0) delete h2d2;
0129     h2d2 = new TH2D("h2d2", "h2d2", xdbins,xdmin,xdmax, 150, 0., 1500.);
0130     if(h2d3!=0) delete h2d3;
0131     h2d3 = new TH2D("h2d3", "h2d3", xdbins,xdmin,xdmax, 150, 0., 1500.);
0132     if(h2d4!=0) delete h2d4;
0133     h2d4 = new TH2D("h2d4", "h2d4", xdbins,xdmin,xdmax, 150, 0., 1500.);
0134     if(h2d5!=0) delete h2d5;
0135     h2d5 = new TH2D("h2d5", "h2d5", xdbins,xdmin,xdmax, 150, 0., 1500.);
0136     if(h2d6!=0) delete h2d6;
0137     h2d6 = new TH2D("h2d6", "h2d6", xdbins,xdmin,xdmax, 150, 0., 1500.);
0138     if(h2d7!=0) delete h2d7;
0139     h2d7 = new TH2D("h2d7", "h2d7", xdbins,xdmin,xdmax, 150, 0., 1500.);
0140     if(h2d8!=0) delete h2d8;
0141     h2d8 = new TH2D("h2d8", "h2d8", xdbins,xdmin,xdmax, 150, 0., 1500.);
0142     if(h2d9!=0) delete h2d9;
0143     h2d9 = new TH2D("h2d9", "h2d9", xdbins,xdmin,xdmax, 150, 0., 1500.);
0144     if(h2d10!=0) delete h2d10;
0145     h2d10 = new TH2D("h2d10", "h2d10", xdbins,xdmin,xdmax, 150, 0., 1500.);
0146     if(h2d11!=0) delete h2d11;
0147     h2d11 = new TH2D("h2d11", "h2d11", xdbins,xdmin,xdmax, 150, 0., 1500.);
0148     if(h2d12!=0) delete h2d12;
0149     h2d12 = new TH2D("h2d12", "h2d12", xdbins,xdmin,xdmax, 150, 0., 1500.);
0150     if(h2d13!=0) delete h2d13;
0151     h2d13 = new TH2D("h2d13", "h2d13", xdbins,xdmin,xdmax, 150, 0., 1500.);
0152     if(h2d14!=0) delete h2d14;
0153     h2d14 = new TH2D("h2d14", "h2d14", xdbins,xdmin,xdmax, 150, 0., 1500.);
0154     if(h2d15!=0) delete h2d15;
0155     h2d15 = new TH2D("h2d15", "h2d15", xdbins,xdmin,xdmax, 150, 0., 1500.);
0156 
0157     if(h2dp0!=0) delete h2dp0;
0158     h2dp0 = new TH2D("h2dp0", "h2dp0", 100, p0min,p0max, 16, 0., 16.);
0159     if(h2dp1!=0) delete h2dp1;
0160     h2dp1 = new TH2D("h2dp1", "h2dp1", 100, 0.,p1max, 16, 0., 16.);
0161     if(h2dp2!=0) delete h2dp2;
0162     h2dp2 = new TH2D("h2dp2", "h2dp2", 100, 0.,p2max, 16, 0., 16.);
0163     if(h2dp3!=0) delete h2dp3;
0164     h2dp3 = new TH2D("h2dp3", "h2dp3", 100, 0.,p3max, 16, 0., 16.);
0165 
0166     if(h2dchis!=0) delete h2dchis;
0167     h2dchis = new TH2D("h2dchis", "h2dchis", 200, 0., 2000., 16, 0., 16.);
0168 
0169   }
0170 
0171   // Init fit function
0172   //phFit = new TF1("phFit", Fitfcn, -400., 1200., nFitParams);
0173   phFit = new TF1("phFit", Fitfcn, 0., 1600., nFitParams);
0174   phLinear = new TF1("phLinear", FitLinear, 0., 1600., 2);
0175 
0176   phFit->SetNpx(1000);
0177   
0178   // Init errors
0179   for (int i = 0; i < 2*vcalSteps; i++) {
0180     xErr[i] = 10.;
0181     yErr[i] = 2.;
0182   }
0183   
0184   // VCAL points used for calibration
0185   vcal[0] = 50.;
0186   vcal[1] = 100.;
0187   vcal[2] = 150.;
0188   vcal[3] = 200.;
0189   vcal[4] = 250.;
0190   vcal[5] = 30.;
0191   vcal[6] = 50.;
0192   vcal[7] = 70.;
0193   vcal[8] = 90.;
0194   vcal[9] = 200.;
0195   // Convert always to low range
0196   for (int i = 0; i < 2*vcalSteps; i++) {
0197     vcalLow[i] = vcal[i];
0198     if (i > (vcalSteps - 1)) vcalLow[i]*=rangeConversion;
0199   }
0200 }
0201 //====================================================================
0202 // Fit a pixel, linear=false to the full Fit, = true for linear fit
0203 void Fit(bool linear) {
0204   bool verbose = false;
0205   
0206   if (graph) delete graph;
0207   //graph = new TGraphErrors(n, x, y, xErr, yErr);
0208   graph = new TGraphErrors(n, y, x, yErr, xErr);
0209   
0210   double xmax = 0., xmin = 9999.; 
0211   double ymax = 0., ymin = 9999.;
0212   for (int i = 0; i < n; i++) {
0213     if (x[i] < xmin) xmin = x[i];
0214     if (x[i] > xmax) xmax = x[i];
0215     if (y[i] < ymin) ymin = y[i];
0216     if (y[i] > ymax) ymax = y[i];
0217   }
0218   //phFit->SetRange(xmin, xmax);
0219   if(!linear) phFit->SetRange(ymin, ymax);
0220   else phLinear->SetRange(ymin, ymax);
0221 
0222   //int upperPoint = vcalSteps+2 - 1;
0223   //int lowerPoint = vcalSteps/3 - 1;
0224   //double slope;
0225   
0226   //if ( (x[upperPoint]-x[lowerPoint]) != 0 ) slope = (y[upperPoint]-y[lowerPoint])/(x[upperPoint]-x[lowerPoint]);
0227   //else slope = 0.5;
0228   
0229   //phFit->SetParameter(2, slope);
0230   //phFit->SetParameter(3, y[upperPoint] - slope*x[upperPoint]);
0231   
0232   if (n < 7 || linear) {   // Use Linear Fit
0233     for (int i = 0; i < 2; i++) phFit->ReleaseParameter(i);
0234     phFit->SetParameter(0, 0.);
0235     phFit->SetParameter(1, 1.);
0236     phFit->FixParameter(2, 0.);
0237     phFit->FixParameter(3, 0.);
0238     if(n<7 && !linear) cout<<" switch to linear "<<n<<endl;
0239 
0240   } else {    // Use full fit
0241 
0242     for (int i = 0; i < nFitParams; i++) phFit->ReleaseParameter(i);
0243 
0244     //double par0 = (TMath::Pi()/2. - 1.4) / x[n-1];
0245     //      printf("par0 %e\n", par0);
0246 
0247 //     phFit->SetParameter(0, par0);
0248 //     phFit->SetParameter(1, 5.e-7);
0249 //     phFit->FixParameter(4, -1.4);
0250 //     if (x[upperPoint] != 0.) phFit->SetParameter(5, (y[upperPoint] - (TMath::Tan(phFit->GetParameter(0)*x[upperPoint] - phFit->GetParameter(4)) + phFit->GetParameter(1)*x[upperPoint]*x[upperPoint]*x[upperPoint] + slope*x[upperPoint] + phFit->GetParameter(3)))/(x[upperPoint]*x[upperPoint]));
0251 //     else phFit->SetParameter(5, 0.);
0252 
0253      phFit->SetParameter(0, 0.003);
0254      phFit->SetParameter(1, 1.);
0255      phFit->SetParameter(2, 100.);
0256      phFit->SetParameter(3, 100.);
0257 
0258   }
0259   
0260   //    for (int i = 0; i < 6; i++) 
0261   //    {
0262 //      if (phFit->GetParameter(i)>0) phFit->SetParLimits(i, 0.8*phFit->GetParameter(i), 1.2*phFit->GetParameter(i));
0263 //      else phFit->SetParLimits(i, 1.2*phFit->GetParameter(i), 0.8*phFit->GetParameter(i));
0264 //  }
0265 //  phFit->SetParLimits(1, 0., 1.e-6);
0266 //  for (int i = 0; i < 6; i++) phFit->FixParameter(i, phFit->GetParameter(i));
0267             
0268   if (verbose) printf("par1 %e\n", (TMath::Pi()/2. - 1.4) / y[n-1]);
0269   if (verbose) printf("par4 %e\n", slope);
0270   if (verbose) printf("x %e y %e\n", x[upperPoint], y[upperPoint]);
0271   if (verbose) printf("x %e y %e\n", x[lowerPoint], y[lowerPoint]); 
0272   if (verbose) printf("par5 %e\n", y[upperPoint] - slope*x[upperPoint]);
0273 
0274   if(linear) {
0275     if (verbose) graph->Fit("phLinear", "R", "");
0276     else graph->Fit("phLinear", "RQ", "");
0277     for (int i = 0; i < 2; i++) 
0278       {histoFit[i]->Fill(phLinear->GetParameter(i));}
0279   } else {
0280     if (verbose) graph->Fit("phFit", "R", "");
0281     else graph->Fit("phFit", "RQ", "");
0282     for (int i = 0; i < nFitParams; i++) 
0283       {histoFit[i]->Fill(phFit->GetParameter(i));}
0284   }
0285 
0286 }
0287 //========================================================================
0288 //Fit all pixels in all ROCs
0289 void FitAllCurves(char *dirName, bool linear = false) {
0290   FILE *inputFile, *outputFile;
0291   char fname[1000], string[500];
0292   int ph[2*vcalSteps], a, b, maxRoc, maxCol, maxRow;
0293   double chiSquare, maxChiSquare = 0.;
0294   int countErrors=0;
0295 
0296   Initialize(linear);
0297   
0298   for (int chip = 0; chip < 16; chip++) {
0299     //for (int chip = 0; chip < 1; chip++) {
0300     printf("Fitting pulse height curves for chip %i\n", chip);
0301     
0302     sprintf(fname, "%s/phCalibration_C%i.dat", dirName, chip);
0303     inputFile = fopen(fname, "r");
0304     if (!inputFile) {
0305       printf("!!!!!!!!!  ----> PHCalibration: Could not open file %s to read pulse height calibration\n", fname);
0306       return;
0307     }
0308     
0309     for (int i = 0; i < 4; i++) fgets(string, 500, inputFile);
0310     
0311     //sprintf(fname, "%s/phCalibrationFit_C%i.dat", dirName, chip);
0312     sprintf(fname, "./phCalibrationFit_C%i.dat", chip);
0313     outputFile = fopen(fname, "w");
0314     if (!outputFile) {
0315       printf("!!!!!!!!!  ----> PHCalibration: Could not open file %s to write the fit results\n", fname);
0316       return;
0317     }
0318     fprintf(outputFile, "Parameters of the vcal vs. pulse height fits\n");
0319     if(!linear) fprintf(outputFile, "%s\n", Fitfcn());
0320     else fprintf(outputFile, "%s\n", FitLinear());
0321     fprintf(outputFile, "\n");
0322     
0323     // Loop over columns
0324     for (int iCol = 0; iCol < 52; iCol++) {
0325       //for (int iCol = 0; iCol < 1; iCol++) {
0326       printf("col %i ", iCol);
0327       for (int iRow = 0; iRow < 80; iRow++) { // Loop over rows
0328     printf(".");fflush(stdout);
0329     //printf("col %i row %i\n", iCol, iRow);
0330     n = 0;
0331     for (int i = 0; i < 2*vcalSteps; i++) { // Loop over VCALs
0332       fscanf(inputFile, "%s", string);
0333         
0334       if (!linear) {  // Full fit
0335 
0336         if (strcmp(string, "N/A") == 0);  // invalid point
0337         else { // valid data point
0338           ph[i] = atoi(string);  // ADC as interger
0339           //printf("ph %i vcal %.0f\n", ph[i], vcalLow[i]);
0340           //x[n] = (double)ph[i]; 
0341           x[n] = (double) transformADC(ph[i]); // transform 
0342           y[n] = vcalLow[i];
0343 
0344           if(HISTO) {
0345         if(chip==0) h2d0->Fill(x[n],y[n]);
0346         else if(chip==1) h2d1->Fill(x[n],y[n]);
0347         else if(chip==2) h2d2->Fill(x[n],y[n]);
0348         else if(chip==3) h2d3->Fill(x[n],y[n]);
0349         else if(chip==4) h2d4->Fill(x[n],y[n]);
0350         else if(chip==5) h2d5->Fill(x[n],y[n]);
0351         else if(chip==6) h2d6->Fill(x[n],y[n]);
0352         else if(chip==7) h2d7->Fill(x[n],y[n]);
0353         else if(chip==8) h2d8->Fill(x[n],y[n]);
0354         else if(chip==9) h2d9->Fill(x[n],y[n]);
0355         else if(chip==10) h2d10->Fill(x[n],y[n]);
0356         else if(chip==11) h2d11->Fill(x[n],y[n]);
0357         else if(chip==12) h2d12->Fill(x[n],y[n]);
0358         else if(chip==13) h2d13->Fill(x[n],y[n]);
0359         else if(chip==14) h2d14->Fill(x[n],y[n]);
0360         else if(chip==15) h2d15->Fill(x[n],y[n]);
0361           }
0362 
0363           n++;
0364 
0365         }
0366 
0367       } else {  // Linear
0368 
0369         //if ((strcmp(string, "N/A") == 0) || (i < 2) || (i > 2*vcalSteps - 2));
0370         if ((strcmp(string, "N/A") == 0) || (i < 1) || (i > 2*vcalSteps - 3));
0371         else {
0372           ph[i] = atoi(string);
0373           //x[n] = (double)ph[i];
0374           x[n] = (double) transformADC(ph[i]); // transform 
0375           y[n] = vcalLow[i];
0376           //printf("ph %i vcal %.0f %f\n", ph[i], vcalLow[i],x[n]);
0377           n++;
0378         }
0379       }
0380     }
0381     fscanf(inputFile, "%s %2i %2i", string, &a, &b);  //comment
0382     
0383     // Do the Fit
0384     if (n != 0) {
0385       Fit(linear);  // Fit
0386       // Check chisq
0387       if(!linear) chiSquare = phFit->GetChisquare()/phFit->GetNDF();
0388       else chiSquare = phLinear->GetChisquare()/phLinear->GetNDF();
0389       histoChi->Fill(chiSquare);
0390       if (chiSquare > maxChiSquare) { // Find chisq
0391         maxChiSquare = chiSquare;
0392         maxRoc = chip;
0393         maxCol = iCol;
0394         maxRow = iRow;
0395       }
0396       if (chiSquare > 6.) {
0397         countErrors++;
0398         cout<<"roc "<<chip<<" col "<<iCol<<" row "<<iRow<<" chisq "
0399         <<chiSquare<<" errors "<<countErrors<<endl;
0400       }
0401       //cout<<iCol<<" "<<iRow<<" " <<chiSquare<<endl;
0402       // Save results in a file 
0403       for (int i = 0; i < nFitParams; i++) {
0404         if(!linear) fprintf(outputFile, "%+e ", phFit->GetParameter(i));
0405         else if(i<2) fprintf(outputFile, "%+e ", phLinear->GetParameter(i));
0406       }
0407 
0408       if(HISTO) {
0409             h2dchis->Fill(chiSquare,float(chip));
0410 
0411         if(!linear) {
0412           h2dp0->Fill(phFit->GetParameter(0),float(chip));
0413           h2dp1->Fill(phFit->GetParameter(1),float(chip));
0414           h2dp2->Fill(phFit->GetParameter(2),float(chip));
0415           h2dp3->Fill(phFit->GetParameter(3),float(chip));
0416           
0417           histoFit0[chip]->Fill(phFit->GetParameter(0));
0418           histoFit1[chip]->Fill(phFit->GetParameter(1));
0419           histoFit2[chip]->Fill(phFit->GetParameter(2));
0420           histoFit3[chip]->Fill(phFit->GetParameter(3));
0421 
0422         } else {
0423 
0424           h2dp0->Fill(phLinear->GetParameter(0),float(chip));
0425           h2dp1->Fill(phLinear->GetParameter(1),float(chip));
0426           
0427           histoFit0[chip]->Fill(phLinear->GetParameter(0));
0428           histoFit1[chip]->Fill(phLinear->GetParameter(1));
0429 
0430           //cout<<phLinear->GetParameter(0)<<" "
0431           //  <<phLinear->GetParameter(1)<<endl;
0432 ;
0433         }
0434           }
0435 
0436 
0437 
0438     } else {  // there may be dead pixels
0439       for (int i = 0; i < nFitParams; i++) {
0440         if(!linear || i<2) fprintf(outputFile, "%+e ", 0.);}
0441     }
0442     fprintf(outputFile, "    Pix %2i %2i\n", iCol, iRow);
0443       }
0444       printf("\n");
0445     }
0446     fclose(inputFile);
0447     fclose(outputFile);
0448   }
0449   printf(" %i %i %i Max ChiSquare/NDF %e\n", maxRoc, maxCol, maxRow, maxChiSquare);
0450 }
0451 //=====================================================================
0452 //Fit one pixel only 
0453 void FitCurve(char *dirName, int chip, int col, int row, 
0454           bool linear = false) {
0455 
0456   FILE *inputFile, *outputFile;
0457   char fname[1000], string[1000];
0458   int ph[2*vcalSteps], a, b;
0459   double chiSquare;
0460   
0461   cout<<"0"<<endl;
0462   Initialize();
0463   
0464   sprintf(fname, "%s/phCalibration_C%i.dat", dirName, chip);
0465   inputFile = fopen(fname, "r");
0466   if (!inputFile) {
0467     printf("!!!!!!!!!  ----> PHCalibration: Could not open file %s to read pulse height calibration\n", fname);
0468     return;
0469   }
0470 
0471   for (int i = 0; i < 4; i++) fgets(string, 1000, inputFile);
0472   for (int i = 0; i < col*80+row; i++) fgets(string, 1000, inputFile);
0473 
0474   n = 0;
0475   for (int i = 0; i < 2*vcalSteps; i++) {
0476     fscanf(inputFile, "%s", string);
0477     if (!linear) {
0478       if (strcmp(string, "N/A") == 0);
0479       else {
0480     ph[i] = atoi(string);
0481     printf("ph %i vcal %.0f\n", ph[i], vcalLow[i]);
0482     //x[n] = (double)ph[i];
0483     x[n] = (double) transformADC(ph[i]);
0484     y[n] = vcalLow[i];
0485     n++;
0486       }
0487     } else {
0488       //if ((strcmp(string, "N/A") == 0) || (i < 2) || (i > 2*vcalSteps - 2));
0489       if ((strcmp(string, "N/A") == 0) || (i < 1) || (i > 2*vcalSteps - 3));
0490       else {
0491     ph[i] = atoi(string);
0492     printf("ph %i vcal %.0f\n", ph[i], vcalLow[i]);
0493     //x[n] = (double)ph[i];
0494     x[n] = (double) transformADC(ph[i]);
0495     y[n] = vcalLow[i];
0496     n++;
0497       }
0498     }
0499   }
0500   fscanf(inputFile, "%s %2i %2i", string, &a, &b);  //comment
0501   
0502   if (n != 0) {
0503       Fit(linear);
0504 
0505       if(!linear) chiSquare = phFit->GetChisquare()/phFit->GetNDF();
0506       else chiSquare = phLinear->GetChisquare()/phLinear->GetNDF();
0507       printf("chiSquare/NDF %e\n", chiSquare);
0508       graph->SetTitle("");
0509       graph->GetYaxis()->SetTitle("Pulse height (ADC units)");
0510       graph->GetYaxis()->SetRangeUser(0.,260.);
0511       graph->GetXaxis()->SetTitle("Vcal (DAC units)");
0512       graph->GetXaxis()->SetTitleOffset(1.2);
0513       graph->GetXaxis()->SetRangeUser(0., 1700.);
0514       graph->Draw("A*");
0515 
0516       cout<<" Fit parameter "<<endl;
0517       for (int i = 0; i < nFitParams; i++) {
0518     if(!linear) cout<<i<<" "<<phFit->GetParameter(i)<<endl;
0519     else if(i<2) cout<<i<<" "<<phLinear->GetParameter(i)<<endl;
0520       }
0521     }
0522   else printf("Error: No measured pulse height values for this pixel\n");
0523   
0524 }
0525 //=========================================================================
0526 // Plots
0527 void Plots() {
0528   TCanvas *c1 = new TCanvas();
0529   c1->Divide(3,2);
0530   
0531   c1->cd(1);
0532   histoFit[0]->Draw();
0533   
0534   c1->cd(2);
0535   histoFit[1]->Draw();
0536   
0537   c1->cd(3);
0538   histoFit[2]->Draw();
0539   
0540   c1->cd(4);
0541   histoFit[3]->Draw();
0542   
0543   c1->cd(6);
0544   histoChi->Draw();
0545 }