SiStrip/macros/noise_analysis.C

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 /*void initMap(map<int, map<double, Entry> > &toinit){
0065   map<double, Entry> dummy;
0066   for(int i=0; i<4; i++)
0067     toinit.insert(make_pair<int, map<double, Entry> >(i, dummy));
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     //cout << ipoint << " " << it->first << " " << it->second.mean() << endl;
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     //cout << regions[i] << endl;
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         // cout << "time to flush? " << !((apv == prev_apv_ && det == prev_det_) || prev_apv_ == 0) <<
0110         //  " apv: " << apv << " prev_apv: " << prev_apv_ << " det: " << det << " prev_det: " << prev_det_ << endl;
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     //cout << "Monitor2D::flush" << endl;
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     //cout << "Filler::add(" << prefix_ << ", " << op_mode_ << ")"<<endl;
0164     EntryMapIT found = emap_[idx].find(length);
0165     if (found == emap_[idx].end()) {
0166       //cout << "adding new entry in map"<<endl;
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   //region, strip length, Entries
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   //book branches (I know, hand-made, I hate it)
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   //tree->SetBranchAddress("ring/i", &ring);
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       //on last event force flushing
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       //Get Region ID
0322       size_t region_ID = 0;  //diagonal by default
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   //TText* clone_info = (TText*) info->Clone("DBTags");
0384   //clone_info->
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     // HMap* hmap = fill_gain.hmap();
0393     // for(int i=0; i<4; i++){
0394     //  for(HMapIT it = hmap[i].begin(); it != hmap[i].end(); ++it){
0395     //      cout << "saving " << it->second.GetName() << endl;
0396     //      it->second.Write();
0397     //  }
0398     // }
0399     // fill_gain.save_harray();
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 }