File indexing completed on 2022-09-22 23:03:21
0001 #include <TROOT.h>
0002 #include <TChain.h>
0003 #include <TFile.h>
0004 #include <TH1D.h>
0005 #include <TH2F.h>
0006 #include <TProfile.h>
0007 #include <TFitResult.h>
0008 #include <TFitResultPtr.h>
0009 #include <TStyle.h>
0010 #include <TCanvas.h>
0011 #include <TLegend.h>
0012 #include <TPaveStats.h>
0013 #include <TPaveText.h>
0014
0015 #include <algorithm>
0016 #include <iomanip>
0017 #include <iostream>
0018 #include <fstream>
0019 #include <sstream>
0020 #include <map>
0021 #include <vector>
0022 #include <string>
0023
0024 void unpackDetId(unsigned int, int&, int&, int&, int&, int&);
0025
0026 #include "CalibMonitor.C"
0027 #include "CalibPlotProperties.C"
0028 #include "CalibTree.C"
0029
0030 int main(Int_t argc, Char_t* argv[]) {
0031 if (argc < 10) {
0032 std::cerr << "Please give N arguments \n"
0033 << "Mode (0 CalibMonitor; 1 CalibProperties; 2 CalibTree)\n"
0034 << "Input File Name\n"
0035 << "Output File Name(ROOT)\n"
0036 << "Correction File Name\n"
0037 << "Directory Name\n"
0038 << "Duplicate File Name\n"
0039 << "Prefix\n"
0040 << "PUcorr\n"
0041 << "Truncate\n"
0042 << "Nmax\n"
0043 << " .... Other parameters depending on mode\n"
0044 << std::endl;
0045 return -1;
0046 }
0047
0048 int mode = std::atoi(argv[1]);
0049 const char* infile = argv[2];
0050 std::string histfile(argv[3]);
0051 int flag = std::atoi(argv[4]);
0052 const char* dirname = argv[5];
0053 std::string prefix(argv[6]);
0054 int pucorr = std::atoi(argv[7]);
0055 int truncate = std::atoi(argv[8]);
0056 Long64_t nmax = static_cast<Long64_t>(std::atoi(argv[9]));
0057
0058 if (mode == 0) {
0059
0060 bool datamc = (argc > 10) ? (std::atoi(argv[10]) < 1) : true;
0061 int numb = (argc > 11) ? std::atoi(argv[11]) : 50;
0062 bool usegen = (argc > 12) ? (std::atoi(argv[12]) < 1) : false;
0063 double scale = (argc > 13) ? std::atof(argv[13]) : 1.0;
0064 int usescale = (argc > 14) ? std::atoi(argv[14]) : 0;
0065 int etalo = (argc > 15) ? std::atoi(argv[15]) : 0;
0066 int etahi = (argc > 16) ? std::atoi(argv[16]) : 30;
0067 int runlo = (argc > 17) ? std::atoi(argv[17]) : 0;
0068 int runhi = (argc > 18) ? std::atoi(argv[18]) : 99999999;
0069 int phimin = (argc > 19) ? std::atoi(argv[19]) : 1;
0070 int phimax = (argc > 20) ? std::atoi(argv[20]) : 72;
0071 int zside = (argc > 21) ? std::atoi(argv[21]) : 1;
0072 int nvxlo = (argc > 22) ? std::atoi(argv[22]) : 0;
0073 int nvxhi = (argc > 23) ? std::atoi(argv[23]) : 1000;
0074 int rbx = (argc > 24) ? std::atoi(argv[24]) : 0;
0075 bool exclude = (argc > 25) ? (std::atoi(argv[25]) > 0) : false;
0076 bool etamax = (argc > 26) ? (std::atoi(argv[26]) > 0) : false;
0077 const char* corrfile = (argc > 27) ? argv[27] : "";
0078 const char* dupfile = (argc > 28) ? argv[28] : "";
0079 const char* comfile = (argc > 29) ? argv[29] : "";
0080 const char* outfile = (argc > 30) ? argv[30] : "";
0081 const char* rcorfile = (argc > 31) ? argv[31] : "";
0082 bool append = (argc > 32) ? (std::atoi(argv[32]) > 0) : true;
0083 bool all = (argc > 33) ? (std::atoi(argv[33]) > 0) : true;
0084 CalibMonitor c1(infile,
0085 dirname,
0086 dupfile,
0087 comfile,
0088 outfile,
0089 prefix,
0090 corrfile,
0091 rcorfile,
0092 pucorr,
0093 flag,
0094 numb,
0095 datamc,
0096 truncate,
0097 usegen,
0098 scale,
0099 usescale,
0100 etalo,
0101 etahi,
0102 runlo,
0103 runhi,
0104 phimin,
0105 phimax,
0106 zside,
0107 nvxlo,
0108 nvxhi,
0109 rbx,
0110 exclude,
0111 etamax);
0112 c1.Loop(nmax);
0113 c1.savePlot(histfile, append, all);
0114 } else if (mode == 1) {
0115
0116 bool datamc = (argc > 10) ? (std::atoi(argv[10]) < 1) : true;
0117 bool usegen = (argc > 11) ? (std::atoi(argv[11]) < 1) : false;
0118 double scale = (argc > 12) ? std::atof(argv[12]) : 1.0;
0119 int usescale = (argc > 13) ? std::atoi(argv[13]) : 0;
0120 int etalo = (argc > 14) ? std::atoi(argv[14]) : 0;
0121 int etahi = (argc > 15) ? std::atoi(argv[15]) : 30;
0122 int runlo = (argc > 16) ? std::atoi(argv[16]) : 0;
0123 int runhi = (argc > 17) ? std::atoi(argv[17]) : 99999999;
0124 int phimin = (argc > 18) ? std::atoi(argv[18]) : 1;
0125 int phimax = (argc > 19) ? std::atoi(argv[19]) : 72;
0126 int zside = (argc > 20) ? std::atoi(argv[20]) : 1;
0127 int nvxlo = (argc > 21) ? std::atoi(argv[21]) : 0;
0128 int nvxhi = (argc > 22) ? std::atoi(argv[22]) : 1000;
0129 int rbx = (argc > 23) ? std::atoi(argv[23]) : 0;
0130 bool exclude = (argc > 24) ? (std::atoi(argv[24]) > 0) : false;
0131 bool etamax = (argc > 25) ? (std::atoi(argv[25]) > 0) : false;
0132 const char* corrfile = (argc > 26) ? argv[26] : "";
0133 const char* dupfile = (argc > 27) ? argv[27] : "";
0134 const char* rcorfile = (argc > 28) ? argv[28] : "";
0135 bool append = (argc > 29) ? (std::atoi(argv[29]) > 0) : true;
0136 bool all = (argc > 30) ? (std::atoi(argv[30]) > 0) : true;
0137 bool debug(false);
0138 CalibPlotProperties c1(infile,
0139 dirname,
0140 dupfile,
0141 prefix,
0142 corrfile,
0143 rcorfile,
0144 pucorr,
0145 flag,
0146 datamc,
0147 truncate,
0148 usegen,
0149 scale,
0150 usescale,
0151 etalo,
0152 etahi,
0153 runlo,
0154 runhi,
0155 phimin,
0156 phimax,
0157 zside,
0158 nvxlo,
0159 nvxhi,
0160 rbx,
0161 exclude,
0162 etamax);
0163 c1.Loop(nmax);
0164 c1.savePlot(histfile, append, all, debug);
0165 } else {
0166
0167 int maxIter = (argc > 10) ? std::atoi(argv[10]) : 30;
0168 const char* corrfile = (argc > 11) ? argv[11] : "";
0169 int applyl1 = (argc > 12) ? std::atoi(argv[12]) : 1;
0170 double l1cut = (argc > 13) ? std::atof(argv[13]) : 0.5;
0171 bool useiter = (argc > 14) ? (std::atoi(argv[14]) < 1) : true;
0172 bool useweight = (argc > 15) ? (std::atoi(argv[15]) < 1) : true;
0173 bool usemean = (argc > 16) ? (std::atoi(argv[16]) < 1) : false;
0174 int nmin = (argc > 17) ? std::atoi(argv[17]) : 0;
0175 bool inverse = (argc > 18) ? (std::atoi(argv[18]) < 1) : true;
0176 double ratmin = (argc > 19) ? std::atof(argv[19]) : 0.25;
0177 double ratmax = (argc > 20) ? std::atof(argv[20]) : 3.0;
0178 int ietamax = (argc > 21) ? std::atoi(argv[21]) : 25;
0179 int ietatrack = (argc > 22) ? std::atoi(argv[22]) : -1;
0180 int sysmode = (argc > 23) ? std::atoi(argv[23]) : -1;
0181 int rcorform = (argc > 24) ? std::atoi(argv[24]) : 0;
0182 bool usegen = (argc > 25) ? (std::atoi(argv[25]) < 1) : false;
0183 int runlo = (argc > 26) ? std::atoi(argv[26]) : 0;
0184 int runhi = (argc > 27) ? std::atoi(argv[27]) : 99999999;
0185 int phimin = (argc > 28) ? std::atoi(argv[28]) : 1;
0186 int phimax = (argc > 29) ? std::atoi(argv[29]) : 72;
0187 int zside = (argc > 30) ? std::atoi(argv[30]) : 0;
0188 int nvxlo = (argc > 31) ? std::atoi(argv[31]) : 0;
0189 int nvxhi = (argc > 32) ? std::atoi(argv[32]) : 1000;
0190 int rbx = (argc > 33) ? std::atoi(argv[33]) : 0;
0191 bool exclude = (argc > 34) ? (std::atoi(argv[34]) > 0) : false;
0192 int higheta = (argc > 35) ? std::atoi(argv[35]) : 1;
0193 double fraction = (argc > 36) ? std::atof(argv[36]) : 1.0;
0194 bool writehisto = (argc > 37) ? (std::atoi(argv[37]) > 0) : false;
0195 bool debug = (argc > 38) ? (std::atoi(argv[38]) > 0) : false;
0196 const char* treename = (argc > 39) ? argv[39] : "CalibTree";
0197 const char* dupfile = (argc > 40) ? argv[40] : "";
0198 const char* rcorfile = (argc > 41) ? argv[41] : "";
0199
0200 char name[500];
0201 sprintf(name, "%s/%s", dirname, treename);
0202 TChain* chain = new TChain(name);
0203 std::cout << "Create a chain for " << name << " from " << infile << std::endl;
0204 if (!fillChain(chain, infile)) {
0205 std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0206 } else {
0207 std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0208 Long64_t nentryTot = chain->GetEntries();
0209 Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0210 static const int maxIterMax = 100;
0211 if (maxIter > maxIterMax)
0212 maxIter = maxIterMax;
0213 std::cout << "Tree " << name << " " << chain << " in directory " << dirname << " from file " << infile
0214 << " with nentries (tracks): " << nentries << std::endl;
0215 unsigned int k(0), kmax(maxIter);
0216 CalibTree t(dupfile,
0217 rcorfile,
0218 truncate,
0219 useiter,
0220 usemean,
0221 runlo,
0222 runhi,
0223 phimin,
0224 phimax,
0225 zside,
0226 nvxlo,
0227 nvxhi,
0228 sysmode,
0229 rbx,
0230 pucorr,
0231 rcorform,
0232 usegen,
0233 exclude,
0234 higheta,
0235 chain);
0236 t.h_pbyE = new TH1D("pbyE", "pbyE", 100, -1.0, 9.0);
0237 t.h_Ebyp_bfr = new TProfile("Ebyp_bfr", "Ebyp_bfr", 60, -30, 30, 0, 10);
0238 t.h_Ebyp_aftr = new TProfile("Ebyp_aftr", "Ebyp_aftr", 60, -30, 30, 0, 10);
0239 t.h_cvg = new TH1D("Cvg0", "Convergence", kmax, 0, kmax);
0240 t.h_cvg->SetMarkerStyle(7);
0241 t.h_cvg->SetMarkerSize(5.0);
0242
0243 TFile* fout = new TFile(histfile.c_str(), "RECREATE");
0244 std::cout << "Output file: " << histfile << " opened in recreate mode" << std::endl;
0245 fout->cd();
0246
0247 double cvgs[maxIterMax], itrs[maxIterMax];
0248 t.getDetId(fraction, ietatrack, debug, nmax);
0249
0250 for (; k <= kmax; ++k) {
0251 std::cout << "Calling Loop() " << k << "th time" << std::endl;
0252 double cvg = t.Loop(k,
0253 fout,
0254 useweight,
0255 nmin,
0256 inverse,
0257 ratmin,
0258 ratmax,
0259 ietamax,
0260 ietatrack,
0261 applyl1,
0262 l1cut,
0263 k == kmax,
0264 fraction,
0265 writehisto,
0266 debug,
0267 nmax);
0268 itrs[k] = k;
0269 cvgs[k] = cvg;
0270 if (cvg < 0.00001)
0271 break;
0272 }
0273
0274 t.writeCorrFactor(corrfile, ietamax);
0275
0276 fout->cd();
0277 TGraph* g_cvg;
0278 g_cvg = new TGraph(k, itrs, cvgs);
0279 g_cvg->SetMarkerStyle(7);
0280 g_cvg->SetMarkerSize(5.0);
0281 g_cvg->Draw("AP");
0282 g_cvg->Write("Cvg");
0283 std::cout << "Finish looping after " << k << " iterations" << std::endl;
0284 t.makeplots(ratmin, ratmax, ietamax, useweight, fraction, debug, nmax);
0285 fout->Close();
0286 }
0287 }
0288
0289 return 0;
0290 }
