File indexing completed on 2024-04-06 12:00:01
0001 #include <iostream>
0002 #include <string>
0003 #include <sstream>
0004 #include <cmath>
0005 #include <boost/tokenizer.hpp>
0006
0007 #include <TChain.h>
0008 #include <TFile.h>
0009 #include <TProfile2D.h>
0010 #include <TH1.h>
0011 #include <TH2.h>
0012 #include <TH3.h>
0013
0014 struct EcalTPGVariables {
0015
0016 unsigned int runNb;
0017 unsigned int evtNb;
0018 unsigned int bxNb;
0019 unsigned int orbitNb;
0020 unsigned int nbOfActiveTriggers;
0021 int activeTriggers[128];
0022
0023
0024 unsigned int nbOfTowers;
0025 int ieta[4032];
0026 int iphi[4032];
0027 int nbOfXtals[4032];
0028 int rawTPData[4032];
0029 int rawTPEmul1[4032];
0030 int rawTPEmul2[4032];
0031 int rawTPEmul3[4032];
0032 int rawTPEmul4[4032];
0033 int rawTPEmul5[4032];
0034 float eRec[4032];
0035 };
0036
0037
0038 void setBranchAddresses(TChain* chain, EcalTPGVariables& treeVars) {
0039 chain->SetBranchAddress("runNb", &treeVars.runNb);
0040 chain->SetBranchAddress("evtNb", &treeVars.evtNb);
0041 chain->SetBranchAddress("bxNb", &treeVars.bxNb);
0042 chain->SetBranchAddress("orbitNb", &treeVars.orbitNb);
0043 chain->SetBranchAddress("nbOfActiveTriggers", &treeVars.nbOfActiveTriggers);
0044 chain->SetBranchAddress("activeTriggers", treeVars.activeTriggers);
0045
0046 chain->SetBranchAddress("nbOfTowers", &treeVars.nbOfTowers);
0047 chain->SetBranchAddress("ieta", treeVars.ieta);
0048 chain->SetBranchAddress("iphi", treeVars.iphi);
0049 chain->SetBranchAddress("nbOfXtals", treeVars.nbOfXtals);
0050 chain->SetBranchAddress("rawTPData", treeVars.rawTPData);
0051 chain->SetBranchAddress("rawTPEmul1", treeVars.rawTPEmul1);
0052 chain->SetBranchAddress("rawTPEmul2", treeVars.rawTPEmul2);
0053 chain->SetBranchAddress("rawTPEmul3", treeVars.rawTPEmul3);
0054 chain->SetBranchAddress("rawTPEmul4", treeVars.rawTPEmul4);
0055 chain->SetBranchAddress("rawTPEmul5", treeVars.rawTPEmul5);
0056 chain->SetBranchAddress("eRec", treeVars.eRec);
0057 }
0058
0059 void printHelp() {
0060 std::cout << "\n Help" << std::endl;
0061 std::cout << " -h : display help" << std::endl;
0062 std::cout << " -i <input root files containing the TPG tree>" << std::endl;
0063 std::cout << " 2 possible formats when more than 1 file:" << std::endl;
0064 std::cout << " - files names separated by comma without blanks: file1.root,file2.root,file3.root" << std::endl;
0065 std::cout << " - files names using the wildcard \":\" file:1:3:.root" << std::endl;
0066 std::cout << " with this last format, the loop from file1.root to file3.root is performed (see examples below)"
0067 << std::endl;
0068 std::cout << " -d <directory containing the input root files, default=ignored>" << std::endl;
0069 std::cout << " -o <output root file name, default=histoTPG.root>" << std::endl;
0070 std::cout << " -v <verbosity level(int), default=0 (minimal)>" << std::endl;
0071 std::cout << " -l1 <L1 algo bits required. If several, use a comma with no blank. Default=select all>" << std::endl;
0072 std::cout << " Additional long options are:" << std::endl;
0073 std::cout << " --cutTPOccup <minimal value of TP to fill the occupancy plot, default=0>" << std::endl;
0074
0075 std::cout << "\n Example:" << std::endl;
0076 std::cout << "1) TPGTreeReder -o myfile -l1 16,46 -i file1,file2" << std::endl;
0077 std::cout << "will produce histo in file myfile, selecting only events triggered " << std::endl;
0078 std::cout << "by algo bit 16 and 46 and using as inputs the file1 and file2" << std::endl;
0079 std::cout << "2) TPGTreeReder -d rfio:/castor/cern.ch/user/p/paganini/67977 -i run67977_:1:18:.root" << std::endl;
0080 std::cout << "will produce histo in the default file histoTPG.root, selecting all events" << std::endl;
0081 std::cout << "and using as inputs the 18 files run67977_1.root to run67977_18.root" << std::endl;
0082 std::cout << "located in the castor directory /castor/cern.ch/user/p/paganini " << std::endl;
0083 }
0084
0085 int getEt(int val) { return (val & 0xff); }
0086
0087 unsigned int getFg(unsigned int val) { return ((val & 0x100) != 0); }
0088
0089 unsigned int getTtf(unsigned int val) { return ((val >> 9) & 0x7); }
0090
0091 std::vector<std::string> split(std::string msg, std::string separator) {
0092 boost::char_separator<char> sep(separator.c_str());
0093 boost::tokenizer<boost::char_separator<char> > tok(msg, sep);
0094 std::vector<std::string> token;
0095 for (boost::tokenizer<boost::char_separator<char> >::const_iterator i = tok.begin(); i != tok.end(); ++i) {
0096 token.push_back(std::string(*i));
0097 }
0098 return token;
0099 }
0100
0101 double getEta(int ietaTower) {
0102
0103 return 0.0174 * std::abs(ietaTower);
0104 }
0105
0106
0107
0108 int main(int argc, char** argv) {
0109 if (argc < 3) {
0110 printHelp();
0111 exit(1);
0112 }
0113
0114 std::string inputfiles, inputdir;
0115 std::string outputRootName = "histoTPG.root";
0116 int verbose = 0;
0117 int occupancyCut = 0;
0118 std::string l1algo;
0119
0120 bool ok(false);
0121 for (int i = 0; i < argc; i++) {
0122 if (argv[i] == std::string("-h")) {
0123 printHelp();
0124 exit(1);
0125 }
0126 if (argv[i] == std::string("-i") && argc > i + 1) {
0127 ok = true;
0128 inputfiles = argv[i + 1];
0129 }
0130 if (argv[i] == std::string("-d") && argc > i + 1)
0131 inputdir = argv[i + 1];
0132 if (argv[i] == std::string("-o") && argc > i + 1)
0133 outputRootName = argv[i + 1];
0134 if (argv[i] == std::string("-v") && argc > i + 1)
0135 verbose = atoi(argv[i + 1]);
0136 if (argv[i] == std::string("-l1") && argc > i + 1)
0137 l1algo = std::string(argv[i + 1]);
0138 if (argv[i] == std::string("--cutTPOccup") && argc > i + 1)
0139 occupancyCut = atoi(argv[i + 1]);
0140 }
0141 if (!ok) {
0142 std::cout << "No input files have been given: nothing to do!" << std::endl;
0143 printHelp();
0144 exit(1);
0145 }
0146
0147 std::vector<int> algobits;
0148 std::vector<std::string> algos = split(l1algo, ",");
0149 for (unsigned int i = 0; i < algos.size(); i++)
0150 algobits.push_back(atoi(algos[i].c_str()));
0151
0152 unsigned int ref = 2;
0153
0154
0155
0156
0157
0158 TH2F* occupancyTP = new TH2F("occupancyTP", "Occupancy TP data", 72, 1, 73, 38, -19, 19);
0159 occupancyTP->GetYaxis()->SetTitle("eta index");
0160 occupancyTP->GetXaxis()->SetTitle("phi index");
0161 TH2F* occupancyTPEmul = new TH2F("occupancyTPEmul", "Occupancy TP emulator", 72, 1, 73, 38, -19, 19);
0162 occupancyTPEmul->GetYaxis()->SetTitle("eta index");
0163 occupancyTPEmul->GetXaxis()->SetTitle("phi index");
0164
0165 TH1F* TP = new TH1F("TP", "TP", 256, 0., 256.);
0166 TP->GetXaxis()->SetTitle("TP (ADC)");
0167 TH1F* TPEmul = new TH1F("TPEmul", "TP Emulator", 256, 0., 256.);
0168 TPEmul->GetXaxis()->SetTitle("TP (ADC)");
0169 TH1F* TPEmulMax = new TH1F("TPEmulMax", "TP Emulator max", 256, 0., 256.);
0170 TPEmulMax->GetXaxis()->SetTitle("TP (ADC)");
0171 TH3F* TPspectrumMap3D = new TH3F("TPspectrumMap3D", "TP data spectrum map", 72, 1, 73, 38, -19, 19, 256, 0., 256.);
0172 TPspectrumMap3D->GetYaxis()->SetTitle("eta index");
0173 TPspectrumMap3D->GetXaxis()->SetTitle("phi index");
0174
0175 TH1F* TPMatchEmul = new TH1F("TPMatchEmul", "TP data matching Emulator", 7, -1., 6.);
0176 TH1F* TPEmulMaxIndex = new TH1F("TPEmulMaxIndex", "Index of the max TP from Emulator", 7, -1., 6.);
0177 TH3I* TPMatchEmul3D = new TH3I("TPMatchEmul3D", "TP data matching Emulator", 72, 1, 73, 38, -19, 19, 7, -1, 6);
0178 TPMatchEmul3D->GetYaxis()->SetTitle("eta index");
0179 TPMatchEmul3D->GetXaxis()->SetTitle("phi index");
0180
0181 TH2I* ttfMismatch = new TH2I("ttfMismatch", "TTF mismatch map", 72, 1, 73, 38, -19, 19);
0182 ttfMismatch->GetYaxis()->SetTitle("eta index");
0183 ttfMismatch->GetXaxis()->SetTitle("phi index");
0184
0185
0186
0187
0188
0189 TChain* chain = new TChain("EcalTPGAnalysis");
0190 std::vector<std::string> files;
0191 if (inputfiles.find(std::string(",")) != std::string::npos)
0192 files = split(inputfiles, ",");
0193 if (inputfiles.find(std::string(":")) != std::string::npos) {
0194 std::vector<std::string> filesbase = split(inputfiles, ":");
0195 if (filesbase.size() == 4) {
0196 int first = atoi(filesbase[1].c_str());
0197 int last = atoi(filesbase[2].c_str());
0198 for (int i = first; i <= last; i++) {
0199 std::stringstream name;
0200 name << filesbase[0] << i << filesbase[3];
0201 files.push_back(name.str());
0202 }
0203 }
0204 }
0205 for (unsigned int i = 0; i < files.size(); i++) {
0206 files[i] = inputdir + "/" + files[i];
0207 std::cout << "Input file: " << files[i] << std::endl;
0208 chain->Add(files[i].c_str());
0209 }
0210
0211 EcalTPGVariables treeVars;
0212 setBranchAddresses(chain, treeVars);
0213
0214 int nEntries = chain->GetEntries();
0215 std::cout << "Number of entries: " << nEntries << std::endl;
0216
0217
0218
0219
0220
0221 for (int entry = 0; entry < nEntries; ++entry) {
0222 chain->GetEntry(entry);
0223 if (entry % 1000 == 0)
0224 std::cout << "------> " << entry + 1 << " entries processed"
0225 << " <------\n";
0226 if (verbose > 0)
0227 std::cout << "Run=" << treeVars.runNb << " Evt=" << treeVars.runNb << std::endl;
0228
0229
0230 bool keep(false);
0231 if (algobits.empty())
0232 keep = true;
0233 for (unsigned int algo = 0; algo < algobits.size(); algo++)
0234 for (unsigned int ntrig = 0; ntrig < treeVars.nbOfActiveTriggers; ntrig++)
0235 if (algobits[algo] == treeVars.activeTriggers[ntrig])
0236 keep = true;
0237 if (!keep)
0238 continue;
0239
0240
0241 for (unsigned int tower = 0; tower < treeVars.nbOfTowers; tower++) {
0242 int tp = getEt(treeVars.rawTPData[tower]);
0243 int emul[5] = {getEt(treeVars.rawTPEmul1[tower]),
0244 getEt(treeVars.rawTPEmul2[tower]),
0245 getEt(treeVars.rawTPEmul3[tower]),
0246 getEt(treeVars.rawTPEmul4[tower]),
0247 getEt(treeVars.rawTPEmul5[tower])};
0248 int maxOfTPEmul = 0;
0249 int indexOfTPEmulMax = -1;
0250 for (int i = 0; i < 5; i++)
0251 if (emul[i] > maxOfTPEmul) {
0252 maxOfTPEmul = emul[i];
0253 indexOfTPEmulMax = i;
0254 }
0255 int ieta = treeVars.ieta[tower];
0256 int iphi = treeVars.iphi[tower];
0257 int nbXtals = treeVars.nbOfXtals[tower];
0258 int ttf = getTtf(treeVars.rawTPData[tower]);
0259
0260 if (verbose > 9 && (tp > 0 || maxOfTPEmul > 0)) {
0261 std::cout << "(phi,eta, Nbxtals)=" << std::dec << iphi << " " << ieta << " " << nbXtals << std::endl;
0262 std::cout << "Data Et, TTF: " << tp << " " << ttf << std::endl;
0263 std::cout << "Emulator: ";
0264 for (int i = 0; i < 5; i++)
0265 std::cout << emul[i] << " ";
0266 std::cout << std::endl;
0267 }
0268
0269
0270 TP->Fill(tp);
0271 TPEmul->Fill(emul[ref]);
0272 TPEmulMax->Fill(maxOfTPEmul);
0273 TPspectrumMap3D->Fill(iphi, ieta, tp);
0274
0275
0276 if (tp > occupancyCut)
0277 occupancyTP->Fill(iphi, ieta);
0278 if (emul[ref] > occupancyCut)
0279 occupancyTPEmul->Fill(iphi, ieta);
0280
0281
0282
0283 if (tp > 0 && nbXtals == 25) {
0284 bool match(false);
0285 for (int i = 0; i < 5; i++) {
0286 if (tp == emul[i]) {
0287 TPMatchEmul->Fill(i + 1);
0288 TPMatchEmul3D->Fill(iphi, ieta, i + 1);
0289 match = true;
0290 }
0291 }
0292 if (!match) {
0293 TPMatchEmul->Fill(-1);
0294 TPMatchEmul3D->Fill(iphi, ieta, -1);
0295 if (verbose > 5) {
0296 std::cout << "MISMATCH" << std::endl;
0297 std::cout << "(phi,eta, Nbxtals)=" << std::dec << iphi << " " << ieta << " " << nbXtals << std::endl;
0298 std::cout << "Data Et, TTF: " << tp << " " << ttf << std::endl;
0299 std::cout << "Emulator: ";
0300 for (int i = 0; i < 5; i++)
0301 std::cout << emul[i] << " ";
0302 std::cout << std::endl;
0303 }
0304 }
0305 }
0306 if (maxOfTPEmul > 0)
0307 TPEmulMaxIndex->Fill(indexOfTPEmulMax + 1);
0308
0309
0310 if ((ttf == 1 || ttf == 3) && nbXtals != 25)
0311 ttfMismatch->Fill(iphi, ieta);
0312
0313 }
0314
0315 }
0316
0317
0318
0319
0320
0321
0322 TProfile2D* TPspectrumMap = TPspectrumMap3D->Project3DProfile("yx");
0323 TPspectrumMap->SetName("TPspectrumMap");
0324
0325
0326 TH2F* TPMatchEmul2D = new TH2F("TPMatchEmul2D", "TP data matching Emulator", 72, 1, 73, 38, -19, 19);
0327 TH2F* TPMatchFraction2D =
0328 new TH2F("TPMatchFraction2D", "TP data: fraction of non-single timing", 72, 1, 73, 38, -19, 19);
0329 TPMatchEmul2D->GetYaxis()->SetTitle("eta index");
0330 TPMatchEmul2D->GetXaxis()->SetTitle("phi index");
0331 TPMatchEmul2D->GetZaxis()->SetRangeUser(-1, 6);
0332 TPMatchFraction2D->GetYaxis()->SetTitle("eta index");
0333 TPMatchFraction2D->GetXaxis()->SetTitle("phi index");
0334 for (int binx = 1; binx <= 72; binx++)
0335 for (int biny = 1; biny <= 38; biny++) {
0336 int maxBinz = 5;
0337 double maxCell = TPMatchEmul3D->GetBinContent(binx, biny, maxBinz);
0338 double totalCell(0);
0339 for (int binz = 1; binz <= 7; binz++) {
0340 double content = TPMatchEmul3D->GetBinContent(binx, biny, binz);
0341 if (content > maxCell) {
0342 maxCell = content;
0343 maxBinz = binz;
0344 }
0345 totalCell += content;
0346 }
0347 if (maxCell <= 0)
0348 maxBinz = 2;
0349 TPMatchEmul2D->SetBinContent(binx, biny, float(maxBinz) - 2.);
0350 double fraction = 0;
0351 if (totalCell > 0)
0352 fraction = 1. - maxCell / totalCell;
0353 TPMatchFraction2D->SetBinContent(binx, biny, fraction);
0354 if (totalCell > maxCell && verbose > 9) {
0355 std::cout << "--->" << std::endl;
0356 for (int binz = 1; binz <= 7; binz++) {
0357 std::cout << "(phi,eta, z): (" << TPMatchEmul3D->GetXaxis()->GetBinLowEdge(binx) << ", "
0358 << TPMatchEmul3D->GetYaxis()->GetBinLowEdge(biny) << ", "
0359 << TPMatchEmul3D->GetZaxis()->GetBinLowEdge(binz)
0360 << ") Content=" << TPMatchEmul3D->GetBinContent(binx, biny, binz)
0361 << ", erro=" << TPMatchEmul3D->GetBinContent(binx, biny, binz) << std::endl;
0362 }
0363 }
0364 }
0365
0366 TFile saving(outputRootName.c_str(), "recreate");
0367 saving.cd();
0368
0369 occupancyTP->Write();
0370 occupancyTPEmul->Write();
0371
0372 TP->Write();
0373 TPEmul->Write();
0374 TPEmulMax->Write();
0375 TPspectrumMap->Write();
0376
0377 TPMatchEmul->Write();
0378 TPMatchEmul3D->Write();
0379 TPEmulMaxIndex->Write();
0380 TPMatchEmul2D->Write();
0381 TPMatchFraction2D->Write();
0382
0383 ttfMismatch->Write();
0384
0385 saving.Close();
0386 delete chain;
0387
0388 return 0;
0389 }