File indexing completed on 2024-04-06 12:03:27
0001 #include <iostream>
0002 #include "TTree.h"
0003 #include "TMath.h"
0004 #include <string>
0005 #include <map>
0006 #include <vector>
0007 #include "TFile.h"
0008 #include "TText.h"
0009 #include "TGraphErrors.h"
0010 #include "TH1F.h"
0011 #include "TH2F.h"
0012 #include "TUUID.h"
0013 #include <sstream>
0014 #include <fstream>
0015 #include <set>
0016 #include <algorithm>
0017
0018 using namespace std;
0019
0020 class Mask {
0021 public:
0022 Mask() : container_() {}
0023 void add(unsigned int id, int apv) { container_.insert(make_pair(id, apv)); }
0024 bool has(unsigned int id, int apv) { return container_.find(make_pair(id, apv)) != container_.end(); }
0025
0026 private:
0027 set<pair<unsigned int, int> > container_;
0028 };
0029
0030 class Entry {
0031 public:
0032 Entry() : entries(0), sum(0), sq_sum(0) {}
0033
0034 double mean() { return sum / entries; }
0035 double std_dev() {
0036 double tmean = mean();
0037 return TMath::Sqrt((sq_sum - entries * tmean * tmean) / (entries - 1));
0038 }
0039 double mean_rms() { return std_dev() / TMath::Sqrt(entries); }
0040
0041 void add(double val) {
0042 entries++;
0043 sum += val;
0044 sq_sum += val * val;
0045 }
0046
0047 void reset() {
0048 entries = 0;
0049 sum = 0;
0050 sq_sum = 0;
0051 }
0052
0053 private:
0054 long int entries;
0055 double sum, sq_sum;
0056 };
0057
0058 typedef map<double, Entry> EntryMap;
0059 typedef EntryMap::iterator EntryMapIT;
0060
0061 typedef map<double, TH1F> HMap;
0062 typedef HMap::iterator HMapIT;
0063
0064
0065
0066
0067
0068
0069
0070 void loadGraph(EntryMap& input_map, TGraphErrors* graph) {
0071 int ipoint = 0;
0072 for (EntryMapIT it = input_map.begin(); it != input_map.end(); ++it) {
0073
0074 graph->SetPoint(ipoint, it->first, it->second.mean());
0075 graph->SetPointError(ipoint, 0., it->second.std_dev());
0076 ipoint++;
0077 }
0078 }
0079
0080 TDirectory* makeGraphs(TFile* file, string dirname, EntryMap* input_map) {
0081 TDirectory* dir = file->mkdir(dirname.c_str());
0082 dir->cd();
0083 string regions[4] = {"TIB", "TID", "TOB", "TEC"};
0084
0085 for (int i = 0; i < 4; i++) {
0086 TGraphErrors* graph = new TGraphErrors();
0087 graph->SetName(regions[i].c_str());
0088
0089 loadGraph(input_map[i], graph);
0090 graph->Write();
0091 }
0092 return dir;
0093 }
0094
0095 enum OpMode { STRIP_BASED, APV_BASED, MODULE_BASED };
0096
0097 class Monitor2D {
0098 public:
0099 Monitor2D(OpMode mode, const char* name, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax)
0100 : entryx_(), entryy_(), mode_(mode), obj_(name, "", nbinsx, xmin, xmax, nbinsy, ymin, ymax) {}
0101
0102 Monitor2D() : entryx_(), entryy_(), mode_(OpMode::STRIP_BASED), obj_() {}
0103
0104 ~Monitor2D() {}
0105
0106 void Fill(int apv, int det, double vx, double vy) {
0107 switch (mode_) {
0108 case (OpMode::APV_BASED):
0109
0110
0111 if (!((apv == prev_apv_ && det == prev_det_) || prev_apv_ == 0)) {
0112 flush();
0113 }
0114 prev_apv_ = apv;
0115 prev_det_ = det;
0116 break;
0117 case (OpMode::MODULE_BASED):
0118 if (!(det == prev_det_ || prev_det_ == 0)) {
0119 flush();
0120 }
0121 prev_det_ = det;
0122 break;
0123 case (OpMode::STRIP_BASED):
0124 flush();
0125 break;
0126 }
0127 entryx_.add(vx);
0128 entryy_.add(vy);
0129 }
0130
0131 void flush() {
0132
0133 obj_.Fill(entryx_.mean(), entryy_.mean());
0134 entryx_.reset();
0135 entryy_.reset();
0136 }
0137
0138 TH2F& hist() {
0139 flush();
0140 return obj_;
0141 }
0142
0143 private:
0144 int prev_apv_ = 0, prev_det_ = 0;
0145 Entry entryx_, entryy_;
0146 OpMode mode_;
0147 TH2F obj_;
0148 };
0149
0150 class Filler {
0151 public:
0152 Filler(string prefix, double hmax = 20) : emap_(), hmap_(), prefix_(prefix), hmax_(hmax) {
0153 string names[] = {
0154 "UNKNOWN", "IB1", "IB2", "OB1", "OB2", "W1A", "W2A", "W3A", "W1B", "W2B", "W3B", "W4", "W5", "W6", "W7"};
0155 for (size_t i = 0; i < 15; i++) {
0156 harray_[i] = TH1F((prefix_ + names[i]).c_str(), "", 100, 0, hmax_);
0157 }
0158 }
0159
0160 ~Filler(){};
0161
0162 void add(int idx, double length, int type, double val) {
0163
0164 EntryMapIT found = emap_[idx].find(length);
0165 if (found == emap_[idx].end()) {
0166
0167 emap_[idx][length] = Entry();
0168 stringstream ss;
0169 ss << prefix_ << regions[idx] << "_" << length;
0170 hmap_[idx].insert(make_pair(length, TH1F(ss.str().c_str(), "", 100, 0, hmax_)));
0171 }
0172
0173 emap_[idx][length].add(val);
0174 hmap_[idx][length].Fill(val);
0175 harray_[type].Fill(val);
0176 }
0177
0178 EntryMap* emap() { return emap_; }
0179 HMap* hmap() { return hmap_; }
0180 void save_harray() {
0181 for (size_t i = 0; i < 15; i++)
0182 harray_[i].Write();
0183 }
0184
0185 private:
0186 EntryMap emap_[4];
0187 HMap hmap_[4];
0188 TH1F harray_[15];
0189 string prefix_;
0190 const string regions[4] = {"TIB", "TID", "TOB", "TEC"};
0191 double hmax_;
0192 };
0193
0194 class TkMap {
0195 public:
0196 TkMap(string name) : detid_(0), counts_(0), file_(name) {}
0197 ~TkMap() { file_.close(); }
0198
0199 void add(unsigned int det) {
0200 if (detid_ && detid_ != det)
0201 flush();
0202 detid_ = det;
0203 counts_++;
0204 }
0205
0206 void flush() {
0207 file_ << detid_ << " " << counts_ << endl;
0208 counts_ = 0;
0209 }
0210
0211 private:
0212 unsigned int detid_, counts_;
0213 ofstream file_;
0214 };
0215
0216 void analyze_noise(string input_file,
0217 string output_file,
0218 bool gsim_,
0219 bool g1_,
0220 bool gratio_,
0221 bool gain_ = false,
0222 Mask mask_ = Mask(),
0223 OpMode op_mode_ = OpMode::STRIP_BASED) {
0224
0225 cout << "Running opts: " << endl
0226 << " input: " << input_file << endl
0227 << " output: " << output_file << endl
0228 << " gsim: " << gsim_ << endl
0229 << " g1: " << g1_ << endl
0230 << " gratio " << gratio_ << endl
0231 << " gain " << gain_ << endl
0232 << " op_mode:" << op_mode_ << endl;
0233 string regions[4] = {"TIB", "TID", "TOB", "TEC"};
0234
0235 Filler fill_gsim("GSim_");
0236 Filler fill_g1("G1_");
0237 Filler fill_gratio("GRatio_");
0238 Filler fill_gain("GAIN_", 2);
0239
0240 cout << "Booking 2D plots" << endl;
0241 string det_types[] = {
0242 "UNKNOWN", "IB1", "IB2", "OB1", "OB2", "W1A", "W2A", "W3A", "W1B", "W2B", "W3B", "W4", "W5", "W6", "W7"};
0243 Monitor2D noise_vs_gain[15][6];
0244 string base("_noise_vs_gain");
0245 for (size_t i = 0; i < 15; i++) {
0246 for (size_t j = 0; j < 6; j++) {
0247 TUUID id;
0248 string idc = id.AsString();
0249 noise_vs_gain[i][j] = Monitor2D(op_mode_, idc.c_str(), 100, 0, 2, 124, 0, 31);
0250 }
0251 }
0252
0253 cout << "Booking Tracker maps" << endl;
0254 TkMap* tkmaps[5];
0255 string region_names[] = {"diagonal", "underflow", "below", "above", "overflow", "masked"};
0256 for (size_t j = 0; j < 5; j++) {
0257 tkmaps[j] = new TkMap(region_names[j + 1] + ".detlist");
0258 }
0259 cout << "Everything booked " << endl;
0260
0261 unsigned int detId, ring, istrip, dtype;
0262 Int_t layer;
0263 float noise, gsim, g1, g2, length;
0264 bool isTIB, isTOB, isTEC, isTID;
0265
0266 TFile* infile = TFile::Open(input_file.c_str());
0267 TTree* tree = (TTree*)infile->Get("treeDump/StripDBTree");
0268
0269
0270 tree->SetBranchAddress("detId", &detId);
0271 tree->SetBranchAddress("noise", &noise);
0272 tree->SetBranchAddress("istrip", &istrip);
0273 tree->SetBranchAddress("detType", &dtype);
0274 tree->SetBranchAddress("gsim", &gsim);
0275 tree->SetBranchAddress("g1", &g1);
0276 tree->SetBranchAddress("g2", &g2);
0277 tree->SetBranchAddress("layer", &layer);
0278
0279 tree->SetBranchAddress("length", &length);
0280 tree->SetBranchAddress("isTIB", &isTIB);
0281 tree->SetBranchAddress("isTOB", &isTOB);
0282 tree->SetBranchAddress("isTEC", &isTEC);
0283 tree->SetBranchAddress("isTID", &isTID);
0284
0285 unsigned long int entries = tree->GetEntries();
0286 int cent = entries / 10;
0287 TH1::AddDirectory(kFALSE);
0288
0289 unsigned int prev_det = 0, prev_apv = 0;
0290 int prev_subdet = -1, prev_type = -1;
0291 double prev_length = -1;
0292 Entry enoise, eg1, egsim, eg2;
0293 bool masked = false;
0294 for (unsigned long int ientry = 0; ientry <= entries; ientry++) {
0295 if (ientry < entries)
0296 tree->GetEntry(ientry);
0297 else {
0298
0299 detId = 0;
0300 istrip = prev_apv * 128 + 100;
0301 }
0302 if (ientry % cent == 0) {
0303 cout << "reading entry " << ientry << " of " << entries << " (" << float(ientry) / entries << ")" << endl;
0304 }
0305 unsigned int idx = 0;
0306
0307 bool flush = false;
0308 switch (op_mode_) {
0309 case (OpMode::APV_BASED):
0310 flush = (prev_det != 0 && prev_apv != istrip / 128);
0311 break;
0312 case (OpMode::MODULE_BASED):
0313 flush = (prev_det != 0 && prev_det != detId);
0314 break;
0315 case (OpMode::STRIP_BASED):
0316 flush = (ientry != 0);
0317 break;
0318 }
0319
0320 if (flush) {
0321
0322 size_t region_ID = 0;
0323 if (masked) {
0324 region_ID = 5;
0325 } else if (enoise.mean() < 1)
0326 region_ID = 1;
0327 else if (enoise.mean() > 30)
0328 region_ID = 4;
0329 else if (eg1.mean() > 0.2 && (enoise.mean() - 2.5 * eg1.mean()) < 0.5)
0330 region_ID = 2;
0331 else if (enoise.mean() > 8.333 * eg1.mean())
0332 region_ID = 3;
0333
0334 if (region_ID >= 1)
0335 tkmaps[region_ID - 1]->add(prev_det);
0336
0337 if (gain_) {
0338 fill_gain.add(prev_subdet, prev_length, prev_type, eg1.mean());
0339 noise_vs_gain[prev_type][region_ID].Fill(prev_apv, prev_det, eg1.mean(), enoise.mean());
0340 }
0341 if (gsim_) {
0342 fill_gsim.add(prev_subdet, prev_length, prev_type, enoise.mean() / egsim.mean());
0343 }
0344 if (g1_) {
0345 fill_g1.add(prev_subdet, prev_length, prev_type, enoise.mean() / eg1.mean());
0346 }
0347 if (gratio_) {
0348 fill_gratio.add(prev_subdet, prev_length, prev_type, (eg1.mean() * eg2.mean() / egsim.mean()) - 1);
0349 }
0350 enoise.reset();
0351 eg1.reset();
0352 egsim.reset();
0353 eg2.reset();
0354 }
0355
0356 masked = mask_.has(detId, istrip / 128);
0357 if (masked && op_mode_ != OpMode::APV_BASED && !gain_)
0358 continue;
0359 if (ientry < entries) {
0360 if (isTOB) {
0361 idx = 2;
0362 } else if (isTEC) {
0363 idx = 3;
0364 } else if (isTID) {
0365 idx = 1;
0366 }
0367
0368 enoise.add(std::min<float>(noise, 30.5));
0369 eg1.add(g1);
0370 egsim.add(gsim);
0371 eg2.add(g2);
0372
0373 prev_det = detId;
0374 prev_apv = istrip / 128;
0375 prev_subdet = idx;
0376 prev_length = length;
0377 prev_type = dtype;
0378 }
0379 }
0380 cout << "loop done" << endl;
0381 TText* info = (TText*)infile->Get("DBTags");
0382 cout << "Got DB Info" << endl;
0383
0384
0385
0386 cout << "Opening output: " << output_file << endl;
0387 TFile* outfile = TFile::Open(output_file.c_str(), "RECREATE");
0388
0389 if (gain_) {
0390 cout << "Saving Gain" << endl;
0391 TDirectory* dir = makeGraphs(outfile, "Gain", fill_gain.emap());
0392
0393
0394
0395
0396
0397
0398
0399
0400 for (size_t i = 0; i < 15; i++) {
0401 dir->mkdir(det_types[i].c_str())->cd();
0402 for (size_t j = 0; j < 6; j++) {
0403 noise_vs_gain[i][j].hist().SetName(region_names[j].c_str());
0404 noise_vs_gain[i][j].hist().Write();
0405 }
0406 }
0407 }
0408 if (gsim_) {
0409 cout << "Saving GSim" << endl;
0410 makeGraphs(outfile, "GSim", fill_gsim.emap());
0411 HMap* hmap = fill_gsim.hmap();
0412 fill_gsim.save_harray();
0413 for (int i = 0; i < 4; i++) {
0414 for (HMapIT it = hmap[i].begin(); it != hmap[i].end(); ++it) {
0415 cout << "saving " << it->second.GetName() << endl;
0416 it->second.Write();
0417 }
0418 }
0419 }
0420 if (g1_) {
0421 cout << "Saving G1" << endl;
0422 makeGraphs(outfile, "G1", fill_g1.emap());
0423 HMap* hmap = fill_g1.hmap();
0424 fill_g1.save_harray();
0425 for (int i = 0; i < 4; i++) {
0426 for (HMapIT it = hmap[i].begin(); it != hmap[i].end(); ++it) {
0427 cout << "saving " << it->second.GetName() << endl;
0428 it->second.Write();
0429 }
0430 }
0431 }
0432 if (gratio_) {
0433 cout << "Saving GRatio" << endl;
0434 makeGraphs(outfile, "GRatio", fill_gratio.emap());
0435 HMap* hmap = fill_gratio.hmap();
0436 fill_gratio.save_harray();
0437 for (int i = 0; i < 4; i++) {
0438 for (HMapIT it = hmap[i].begin(); it != hmap[i].end(); ++it) {
0439 cout << "saving " << it->second.GetName() << endl;
0440 it->second.Write();
0441 }
0442 }
0443 }
0444 outfile->Write();
0445 outfile->Close();
0446 infile->Close();
0447 for (size_t j = 0; j < 4; j++) {
0448 delete tkmaps[j];
0449 }
0450 }