Back to home page

Project CMSSW displayed by LXR

 
 

    


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   // event variables
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   // tower variables
0024   unsigned int nbOfTowers;  //max 4032 EB+EE
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 //! branch addresses settings
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   // Paga: to be confirmed, specially in EE:
0103   return 0.0174 * std::abs(ietaTower);
0104 }
0105 
0106 ///////  Main program /////////
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   // book the histograms
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   // Chain the trees:
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   // Main loop over entries
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     // trigger selection if any
0230     bool keep(false);
0231     if (algobits.empty())
0232       keep = true;  // keep all events when no trigger selection
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     // loop on towers
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       // Fill TP spctrum
0270       TP->Fill(tp);
0271       TPEmul->Fill(emul[ref]);
0272       TPEmulMax->Fill(maxOfTPEmul);
0273       TPspectrumMap3D->Fill(iphi, ieta, tp);
0274 
0275       // Fill TP occupancy
0276       if (tp > occupancyCut)
0277         occupancyTP->Fill(iphi, ieta);
0278       if (emul[ref] > occupancyCut)
0279         occupancyTPEmul->Fill(iphi, ieta);
0280 
0281       // Fill TP-Emulator matching
0282       // comparison is meaningful when:
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       // Fill TTF mismatch
0310       if ((ttf == 1 || ttf == 3) && nbXtals != 25)
0311         ttfMismatch->Fill(iphi, ieta);
0312 
0313     }  // end loop towers
0314 
0315   }  // endloop entries
0316 
0317   ///////////////////////
0318   // Format & write histos
0319   ///////////////////////
0320 
0321   // 1. TP Spectrum
0322   TProfile2D* TPspectrumMap = TPspectrumMap3D->Project3DProfile("yx");
0323   TPspectrumMap->SetName("TPspectrumMap");
0324 
0325   // 2. TP Timing
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;                                                  // empty cell
0349       TPMatchEmul2D->SetBinContent(binx, biny, float(maxBinz) - 2.);  //z must be in [-1,5]
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 }