Back to home page

Project CMSSW displayed by LXR

 
 

    


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     // CalibMonitor
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     // CalibPlotProperties
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     // CalibTree
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 }