File indexing completed on 2024-04-06 11:57:53
0001
0002
0003
0004
0005
0006
0007
0008 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TShapeAnalysis.h"
0009 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/TFParams.h"
0010
0011 #include <iostream>
0012 #include <cmath>
0013 #include <ctime>
0014 #include <cassert>
0015
0016 #include "TFile.h"
0017 #include "TTree.h"
0018 #include "TBranch.h"
0019
0020
0021
0022 using namespace std;
0023
0024
0025 TShapeAnalysis::TShapeAnalysis(double alpha0, double beta0, double width0, double chi20) {
0026 init(alpha0, beta0, width0, chi20);
0027 }
0028
0029 TShapeAnalysis::TShapeAnalysis(TTree *tAB, double alpha0, double beta0, double width0, double chi20) {
0030 init(tAB, alpha0, beta0, width0, chi20);
0031 }
0032
0033
0034 TShapeAnalysis::~TShapeAnalysis() {}
0035
0036 void TShapeAnalysis::init(double alpha0, double beta0, double width0, double chi20) {
0037 tABinit = nullptr;
0038 nchsel = fNchsel;
0039 for (int cris = 0; cris < fNchsel; cris++) {
0040 index[cris] = -1;
0041 npass[cris] = 0;
0042 npassok[cris] = 0.;
0043
0044 alpha_val[cris] = alpha0;
0045 beta_val[cris] = beta0;
0046 width_val[cris] = width0;
0047 chi2_val[cris] = chi20;
0048 flag_val[cris] = 0;
0049
0050 alpha_init[cris] = alpha0;
0051 beta_init[cris] = beta0;
0052 width_init[cris] = width0;
0053 chi2_init[cris] = chi20;
0054 flag_init[cris] = 0;
0055
0056 phi_init[cris] = 0;
0057 eta_init[cris] = 0;
0058 side_init[cris] = 0;
0059 dcc_init[cris] = 0;
0060 tower_init[cris] = 0;
0061 ch_init[cris] = 0;
0062 assignChannel(cris, cris);
0063 }
0064 }
0065
0066 void TShapeAnalysis::init(TTree *tAB, double alpha0, double beta0, double width0, double chi20) {
0067 init(alpha0, beta0, width0, chi20);
0068
0069 tABinit = tAB->CloneTree();
0070
0071
0072 Int_t sidei;
0073 Int_t iphii;
0074 Int_t ietai;
0075 Int_t dccIDi;
0076 Int_t towerIDi;
0077 Int_t channelIDi;
0078 Double_t alphai;
0079 Double_t betai;
0080 Double_t widthi;
0081 Double_t chi2i;
0082 Int_t flagi;
0083
0084
0085 TBranch *b_iphi;
0086 TBranch *b_ieta;
0087 TBranch *b_side;
0088 TBranch *b_dccID;
0089 TBranch *b_towerID;
0090 TBranch *b_channelID;
0091 TBranch *b_alpha;
0092 TBranch *b_beta;
0093 TBranch *b_width;
0094 TBranch *b_chi2;
0095 TBranch *b_flag;
0096
0097 if (tABinit) {
0098 tABinit->SetBranchAddress("iphi", &iphii, &b_iphi);
0099 tABinit->SetBranchAddress("ieta", &ietai, &b_ieta);
0100 tABinit->SetBranchAddress("side", &sidei, &b_side);
0101 tABinit->SetBranchAddress("dccID", &dccIDi, &b_dccID);
0102 tABinit->SetBranchAddress("towerID", &towerIDi, &b_towerID);
0103 tABinit->SetBranchAddress("channelID", &channelIDi, &b_channelID);
0104 tABinit->SetBranchAddress("alpha", &alphai, &b_alpha);
0105 tABinit->SetBranchAddress("beta", &betai, &b_beta);
0106 tABinit->SetBranchAddress("width", &widthi, &b_width);
0107 tABinit->SetBranchAddress("chi2", &chi2i, &b_chi2);
0108 tABinit->SetBranchAddress("flag", &flagi, &b_flag);
0109
0110 nchsel = tABinit->GetEntries();
0111 assert(nchsel <= fNchsel);
0112
0113 for (int cris = 0; cris < nchsel; cris++) {
0114 tABinit->GetEntry(cris);
0115
0116
0117
0118 putalphaVal(cris, alphai);
0119 putchi2Val(cris, chi2i);
0120 putbetaVal(cris, betai);
0121 putwidthVal(cris, widthi);
0122 putflagVal(cris, flagi);
0123
0124 putalphaInit(cris, alphai);
0125 putchi2Init(cris, chi2i);
0126 putbetaInit(cris, betai);
0127 putwidthInit(cris, widthi);
0128 putflagInit(cris, flagi);
0129 putetaInit(cris, ietai);
0130 putphiInit(cris, iphii);
0131 }
0132 }
0133 }
0134
0135 void TShapeAnalysis::set_const(int ns, int ns1, int ns2, int ps, int nevtmax, double noise_val, double chi2_cut) {
0136 nsamplecristal = ns;
0137 presample = ps;
0138 sampbmax = ns1;
0139 sampamax = ns2;
0140 nevt = nevtmax;
0141 noise = noise_val;
0142 chi2cut = chi2_cut;
0143 }
0144
0145 void TShapeAnalysis::set_presample(int ps) { presample = ps; }
0146 void TShapeAnalysis::set_nch(int nch) {
0147 assert(nch <= fNchsel);
0148 if (tABinit)
0149 assert(nch == nchsel);
0150 nchsel = nch;
0151 }
0152 void TShapeAnalysis::assignChannel(int n, int ch) {
0153 if (n >= nchsel)
0154 printf(" number of channels exceed maximum allowed\n");
0155
0156 index[n] = ch;
0157 }
0158
0159 void TShapeAnalysis::putDateStart(long int timecur) { timestart = timecur; }
0160
0161 void TShapeAnalysis::putDateStop(long int timecur) { timestop = timecur; }
0162
0163 void TShapeAnalysis::getDateStart() {
0164 time_t t, timecur;
0165 timecur = time(&t);
0166 timestart = ((long int)timecur);
0167 }
0168
0169 void TShapeAnalysis::getDateStop() {
0170 time_t t, timecur;
0171 timecur = time(&t);
0172 timestop = ((long int)timecur);
0173 }
0174
0175 void TShapeAnalysis::putAllVals(int ch, double *sampl, int ieta, int iphi, int dcc, int side, int tower, int chid) {
0176 dcc_init[ch] = dcc;
0177 tower_init[ch] = side;
0178 ch_init[ch] = chid;
0179 side_init[ch] = side;
0180 eta_init[ch] = ieta;
0181 phi_init[ch] = iphi;
0182 putAllVals(ch, sampl, ieta, iphi);
0183 }
0184
0185 void TShapeAnalysis::putAllVals(int ch, double *sampl, int ieta, int iphi) {
0186 int i, k;
0187 int n = -1;
0188 for (i = 0; i < nchsel; i++)
0189 if (index[i] == ch)
0190 n = i;
0191
0192 if (n >= 0) {
0193 if (npass[n] < nevt) {
0194 for (k = 0; k < nsamplecristal; k++) {
0195 rawsglu[n][npass[n]][k] = sampl[k];
0196 }
0197
0198 npass[n]++;
0199 }
0200 } else {
0201 printf("no index found for ch=%d\n", ch);
0202 }
0203 }
0204
0205 void TShapeAnalysis::computeShape(string namefile, TTree *tAB) {
0206 double tm_atmax[200];
0207 double parout[3];
0208
0209 double chi2_all = 0.;
0210
0211 double **dbi;
0212 dbi = new double *[200];
0213 for (int k = 0; k < 200; k++)
0214 dbi[k] = new double[2];
0215
0216 double **signalu;
0217 signalu = new double *[200];
0218 for (int k = 0; k < 200; k++)
0219 signalu[k] = new double[10];
0220
0221 TFParams *pjf = new TFParams();
0222
0223 for (int i = 0; i < nchsel; i++) {
0224 if (index[i] >= 0) {
0225 if (npass[i] <= 10) {
0226 putalphaVal(i, alpha_init[i]);
0227 putbetaVal(i, beta_init[i]);
0228 putwidthVal(i, width_init[i]);
0229 putchi2Val(i, chi2_init[i]);
0230 putflagVal(i, 0);
0231
0232 } else {
0233 pjf->set_const(nsamplecristal, sampbmax, sampamax, alpha_init[i], beta_init[i], npass[i]);
0234
0235 for (int pass = 0; pass < npass[i]; pass++) {
0236 double ped = 0;
0237 for (int k = 0; k < presample; k++) {
0238 ped += rawsglu[i][pass][k];
0239 }
0240 ped /= double(presample);
0241
0242 for (int k = 0; k < nsamplecristal; k++) {
0243 signalu[pass][k] = rawsglu[i][pass][k] - ped;
0244 }
0245 }
0246
0247 int debug = 0;
0248 chi2_all = pjf->fitpj(signalu, &parout[0], dbi, noise, debug);
0249
0250 if (parout[0] >= 0.0 && parout[1] >= 0.0 && chi2_all <= chi2cut && chi2_all > 0.0) {
0251 putalphaVal(i, parout[0]);
0252 putbetaVal(i, parout[1]);
0253 putchi2Val(i, chi2_all);
0254 putflagVal(i, 1);
0255
0256 } else {
0257 putalphaVal(i, alpha_init[i]);
0258 putbetaVal(i, beta_init[i]);
0259 putwidthVal(i, width_init[i]);
0260 putchi2Val(i, chi2_init[i]);
0261 putflagVal(i, 0);
0262 }
0263
0264 for (int kj = 0; kj < npass[i]; kj++) {
0265 tm_atmax[kj] = dbi[kj][1];
0266 }
0267 computetmaxVal(i, &tm_atmax[0]);
0268 }
0269 }
0270 }
0271
0272 if (tAB)
0273 tABinit = tAB->CloneTree();
0274
0275
0276 Int_t sidei;
0277 Int_t iphii;
0278 Int_t ietai;
0279 Int_t dccIDi;
0280 Int_t towerIDi;
0281 Int_t channelIDi;
0282 Double_t alphai;
0283 Double_t betai;
0284 Double_t widthi;
0285 Double_t chi2i;
0286 Int_t flagi;
0287
0288
0289 TBranch *b_iphi;
0290 TBranch *b_ieta;
0291 TBranch *b_side;
0292 TBranch *b_dccID;
0293 TBranch *b_towerID;
0294 TBranch *b_channelID;
0295 TBranch *b_alpha;
0296 TBranch *b_beta;
0297 TBranch *b_width;
0298 TBranch *b_chi2;
0299 TBranch *b_flag;
0300
0301 if (tABinit) {
0302 tABinit->SetBranchAddress("iphi", &iphii, &b_iphi);
0303 tABinit->SetBranchAddress("ieta", &ietai, &b_ieta);
0304 tABinit->SetBranchAddress("side", &sidei, &b_side);
0305 tABinit->SetBranchAddress("dccID", &dccIDi, &b_dccID);
0306 tABinit->SetBranchAddress("towerID", &towerIDi, &b_towerID);
0307 tABinit->SetBranchAddress("channelID", &channelIDi, &b_channelID);
0308 tABinit->SetBranchAddress("alpha", &alphai, &b_alpha);
0309 tABinit->SetBranchAddress("beta", &betai, &b_beta);
0310 tABinit->SetBranchAddress("width", &widthi, &b_width);
0311 tABinit->SetBranchAddress("chi2", &chi2i, &b_chi2);
0312 tABinit->SetBranchAddress("flag", &flagi, &b_flag);
0313 }
0314
0315 TFile *fABout = new TFile(namefile.c_str(), "RECREATE");
0316 tABout = new TTree("ABCol0", "ABCol0");
0317
0318
0319 Int_t side;
0320 Int_t iphi;
0321 Int_t ieta;
0322 Int_t dccID;
0323 Int_t towerID;
0324 Int_t channelID;
0325 Double_t alpha;
0326 Double_t beta;
0327 Double_t width;
0328 Double_t chi2;
0329 Int_t flag;
0330
0331 tABout->Branch("iphi", &iphi, "iphi/I");
0332 tABout->Branch("ieta", &ieta, "ieta/I");
0333 tABout->Branch("side", &side, "side/I");
0334 tABout->Branch("dccID", &dccID, "dccID/I");
0335 tABout->Branch("towerID", &towerID, "towerID/I");
0336 tABout->Branch("channelID", &channelID, "channelID/I");
0337 tABout->Branch("alpha", &alpha, "alpha/D");
0338 tABout->Branch("beta", &beta, "beta/D");
0339 tABout->Branch("width", &width, "width/D");
0340 tABout->Branch("chi2", &chi2, "chi2/D");
0341 tABout->Branch("flag", &flag, "flag/I");
0342
0343 tABout->SetBranchAddress("ieta", &ieta);
0344 tABout->SetBranchAddress("iphi", &iphi);
0345 tABout->SetBranchAddress("side", &side);
0346 tABout->SetBranchAddress("dccID", &dccID);
0347 tABout->SetBranchAddress("towerID", &towerID);
0348 tABout->SetBranchAddress("channelID", &channelID);
0349 tABout->SetBranchAddress("alpha", &alpha);
0350 tABout->SetBranchAddress("beta", &beta);
0351 tABout->SetBranchAddress("width", &width);
0352 tABout->SetBranchAddress("chi2", &chi2);
0353 tABout->SetBranchAddress("flag", &flag);
0354
0355 for (int i = 0; i < nchsel; i++) {
0356 if (tABinit) {
0357 tABinit->GetEntry(i);
0358 iphi = iphii;
0359 ieta = ietai;
0360 side = sidei;
0361 dccID = dccIDi;
0362 towerID = towerIDi;
0363 channelID = channelIDi;
0364
0365 } else {
0366 iphi = phi_init[i];
0367 ieta = eta_init[i];
0368 side = side_init[i];
0369 dccID = dcc_init[i];
0370 towerID = tower_init[i];
0371 channelID = ch_init[i];
0372 }
0373
0374 alpha = alpha_val[i];
0375 beta = beta_val[i];
0376 width = width_val[i];
0377 chi2 = chi2_val[i];
0378 flag = flag_val[i];
0379
0380 tABout->Fill();
0381 }
0382
0383 tABout->Write();
0384 fABout->Close();
0385
0386 delete pjf;
0387 }
0388
0389 void TShapeAnalysis::computetmaxVal(int i, double *tm_val) {
0390 double tm_mean = 0;
0391
0392 double tm = 0.;
0393 for (int k = 0; k < npass[i] - 1; k++) {
0394 if (1. < tm_val[k] && tm_val[k] < 10.) {
0395 npassok[i]++;
0396 tm += tm_val[k];
0397 }
0398 }
0399 if (npassok[i] <= 0) {
0400 tm_mean = 0.;
0401 } else {
0402 for (int k = 0; k < npass[i] - 1; k++) {
0403 if (1. < tm_val[k] && tm_val[k] < 10.) {
0404 tm_mean = tm / npassok[i];
0405 }
0406 }
0407 }
0408
0409 putwidthVal(i, tm_mean);
0410 }
0411
0412 void TShapeAnalysis::putalphaVal(int n, double val) { alpha_val[n] = val; }
0413
0414 void TShapeAnalysis::putchi2Val(int n, double val) { chi2_val[n] = val; }
0415 void TShapeAnalysis::putbetaVal(int n, double val) { beta_val[n] = val; }
0416
0417 void TShapeAnalysis::putwidthVal(int n, double val) { width_val[n] = val; }
0418
0419 void TShapeAnalysis::putflagVal(int n, int val) { flag_val[n] = val; }
0420
0421 void TShapeAnalysis::putalphaInit(int n, double val) { alpha_init[n] = val; }
0422
0423 void TShapeAnalysis::putchi2Init(int n, double val) { chi2_init[n] = val; }
0424 void TShapeAnalysis::putbetaInit(int n, double val) { beta_init[n] = val; }
0425
0426 void TShapeAnalysis::putwidthInit(int n, double val) { width_init[n] = val; }
0427
0428 void TShapeAnalysis::putetaInit(int n, int val) { eta_init[n] = val; }
0429
0430 void TShapeAnalysis::putphiInit(int n, int val) { phi_init[n] = val; }
0431
0432 void TShapeAnalysis::putflagInit(int n, int val) { flag_init[n] = val; }
0433 std::vector<double> TShapeAnalysis::getVals(int n) {
0434 std::vector<double> v;
0435
0436 v.push_back(alpha_val[n]);
0437 v.push_back(beta_val[n]);
0438 v.push_back(width_val[n]);
0439 v.push_back(chi2_val[n]);
0440 v.push_back(flag_val[n]);
0441
0442 return v;
0443 }
0444 std::vector<double> TShapeAnalysis::getInitVals(int n) {
0445 std::vector<double> v;
0446
0447 v.push_back(alpha_init[n]);
0448 v.push_back(beta_init[n]);
0449 v.push_back(width_init[n]);
0450 v.push_back(chi2_init[n]);
0451 v.push_back(flag_init[n]);
0452
0453 return v;
0454 }
0455
0456 void TShapeAnalysis::printshapeData(int gRunNumber) {
0457 FILE *fd;
0458 int nev;
0459 sprintf(filename, "runABW%d.pedestal", gRunNumber);
0460 fd = fopen(filename, "w");
0461 if (fd == nullptr)
0462 printf("Error while opening file : %s\n", filename);
0463
0464 for (int i = 0; i < nchsel; i++) {
0465 if (index[i] >= 0) {
0466 nev = (int)npassok[i];
0467 double trise = alpha_val[i] * beta_val[i];
0468 fprintf(fd,
0469 "%d %d 1 %ld %ld %f %f %f %f\n",
0470 index[i],
0471 nev,
0472 timestart,
0473 timestop,
0474 alpha_val[i],
0475 beta_val[i],
0476 trise,
0477 width_val[i]);
0478 }
0479 }
0480 int iret = fclose(fd);
0481 printf(" Closing file : %d\n", iret);
0482 }
