File indexing completed on 2024-04-06 12:03:21
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
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
0022 static const double rangeConversion = 7.;
0023
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
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
0041 TH2D *h2d0=0,*h2d1,*h2d2,*h2d3,*h2d4,*h2d5,*h2d6,*h2d7,*h2d8,*h2d9,
0042 *h2d10,*h2d11,*h2d12,*h2d13,*h2d14,*h2d15;
0043
0044 TH2D *h2dp0, *h2dp1, *h2dp2, *h2dp3, *h2dchis;
0045
0046
0047
0048
0049
0050
0051 Double_t Fitfcn( Double_t *x, Double_t *par) {
0052
0053
0054
0055
0056 Double_t y = par[3] + par[2] * TMath::TanH(par[0]*x[0]-par[1]);
0057 return y;
0058 }
0059
0060
0061 Double_t FitLinear( Double_t *x, Double_t *par) {
0062
0063
0064
0065
0066 Double_t y = par[0] + par[1] * x[0];
0067 return y;
0068 }
0069
0070
0071 const char *Fitfcn() {
0072
0073
0074 return "par[3] + par[2] * TMath::TanH(par[0]*x[0] - par[1])";
0075 }
0076
0077 const char *FitLinear() {
0078
0079
0080 return "par[0] + par[1] * x[0]";
0081 }
0082
0083
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
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
0172
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
0179 for (int i = 0; i < 2*vcalSteps; i++) {
0180 xErr[i] = 10.;
0181 yErr[i] = 2.;
0182 }
0183
0184
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
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
0203 void Fit(bool linear) {
0204 bool verbose = false;
0205
0206 if (graph) delete graph;
0207
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
0219 if(!linear) phFit->SetRange(ymin, ymax);
0220 else phLinear->SetRange(ymin, ymax);
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232 if (n < 7 || linear) {
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 {
0241
0242 for (int i = 0; i < nFitParams; i++) phFit->ReleaseParameter(i);
0243
0244
0245
0246
0247
0248
0249
0250
0251
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
0261
0262
0263
0264
0265
0266
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
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
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
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
0324 for (int iCol = 0; iCol < 52; iCol++) {
0325
0326 printf("col %i ", iCol);
0327 for (int iRow = 0; iRow < 80; iRow++) {
0328 printf(".");fflush(stdout);
0329
0330 n = 0;
0331 for (int i = 0; i < 2*vcalSteps; i++) {
0332 fscanf(inputFile, "%s", string);
0333
0334 if (!linear) {
0335
0336 if (strcmp(string, "N/A") == 0);
0337 else {
0338 ph[i] = atoi(string);
0339
0340
0341 x[n] = (double) transformADC(ph[i]);
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 {
0368
0369
0370 if ((strcmp(string, "N/A") == 0) || (i < 1) || (i > 2*vcalSteps - 3));
0371 else {
0372 ph[i] = atoi(string);
0373
0374 x[n] = (double) transformADC(ph[i]);
0375 y[n] = vcalLow[i];
0376
0377 n++;
0378 }
0379 }
0380 }
0381 fscanf(inputFile, "%s %2i %2i", string, &a, &b);
0382
0383
0384 if (n != 0) {
0385 Fit(linear);
0386
0387 if(!linear) chiSquare = phFit->GetChisquare()/phFit->GetNDF();
0388 else chiSquare = phLinear->GetChisquare()/phLinear->GetNDF();
0389 histoChi->Fill(chiSquare);
0390 if (chiSquare > maxChiSquare) {
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
0402
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
0431
0432 ;
0433 }
0434 }
0435
0436
0437
0438 } else {
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
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
0483 x[n] = (double) transformADC(ph[i]);
0484 y[n] = vcalLow[i];
0485 n++;
0486 }
0487 } else {
0488
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
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);
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
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 }