Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-08-06 22:36:16

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 //  This provides a provision of executing CalibMonitor, CalibProperties,
0004 //  CalibTree, or CalibSplit in batch operation.
0005 //
0006 //  Usage:
0007 //  ./calibMain.exe <mode> <other parameters depending on mode>
0008 //  mode = 0:CalibMonitor, 1:CalibProperties, 2:CalibTree; 3: CalibSplit
0009 //
0010 //  Other parameters for CalibMonitor:
0011 //  <InputFile> <HistogramFile> <Flag> <DirectoryName> <Prefix> <PUcorr>
0012 //  <Truncate> <Nmax> <datamc> <numb> <usegen> <scale> <usescale> <etalo>
0013 //  <etahi> <runlo> <runhi> <phimin> <phimax> <zside> <nvxlo> <nvxhi>
0014 //  <exclude> <etamax> <append> <all> <corrfile> <rcorfile> <dupfile>
0015 //  <rbxfile> <comfile> <outfile>
0016 //
0017 //  Other parameters for CalibProperties:
0018 //  <InputFile> <HistogramFile> <Flag> <DirectoryName> <Prefix> <PUcorr>
0019 //  <Truncate> <Nmax> <datamc> <usegen> <scale> <usescale> <etalo> <etahi>
0020 //  <runlo> <runhi> <phimin> <phimax> <zside> <nvxlo> <nvxhi> <exclude>
0021 //  <etamax> <append> <all> <corrfile> <rcorfile> <dupfile> <rbxfile>
0022 //
0023 //  Other parameters for CalibTree:
0024 //  <InputFile> <OutputFile> <Flag> <DirectoryName> <Prefix> <PUcorr>
0025 //  <Truncate> <Nmax> <maxIter> <corrfile> <applyl1> <l1cut> <useiter>
0026 //  <useweight> <usemean> <nmin> <inverse> <ratmin> <ratmax> <ietamax>
0027 //  <ietatrack> <sysmode> <rcorform> <usegen> <runlo> <runhi> <phimin>
0028 //  <phimax> <zside> <nvxlo> <nvxhi> <exclude> <higheta> <fraction>
0029 //  <writehisto> <debug> <rcorfile> <dupfile> <rbxfile> <treename>
0030 //
0031 //  Other parameters for CalibSplit:
0032 //  <InputFile> <HistogramFile> <Flag> <DirectoryName> <Prefix> <PUcorr>
0033 //  <Truncate> <Nmax> <pmin> <pmax> <runMin> <runMax> <debug>
0034 //
0035 //
0036 ////////////////////////////////////////////////////////////////////////////////
0037 
0038 #include <TROOT.h>
0039 #include <TChain.h>
0040 #include <TFile.h>
0041 #include <TH1D.h>
0042 #include <TH2F.h>
0043 #include <TProfile.h>
0044 #include <TFitResult.h>
0045 #include <TFitResultPtr.h>
0046 #include <TStyle.h>
0047 #include <TCanvas.h>
0048 #include <TLegend.h>
0049 #include <TPaveStats.h>
0050 #include <TPaveText.h>
0051 
0052 #include <algorithm>
0053 #include <iomanip>
0054 #include <iostream>
0055 #include <fstream>
0056 #include <sstream>
0057 #include <map>
0058 #include <vector>
0059 #include <string>
0060 
0061 void unpackDetId(unsigned int, int&, int&, int&, int&, int&);
0062 
0063 #include "CalibMonitor.C"
0064 #include "CalibPlotProperties.C"
0065 #include "CalibTree.C"
0066 
0067 int main(Int_t argc, Char_t* argv[]) {
0068   if (argc < 10) {
0069     std::cerr << "Please give N arguments \n"
0070               << "Mode (0 CalibMonitor; 1 CalibProperties; 2 CalibTree; 3 CalibSplit)\n"
0071               << "Input File Name\n"
0072               << "Output File Name(ROOT)\n"
0073               << "Flag\n"
0074               << "Directory Name\n"
0075               << "Prefix\n"
0076               << "PUcorr\n"
0077               << "Truncate\n"
0078               << "Nmax\n"
0079               << " .... Other parameters depending on mode\n"
0080               << std::endl;
0081     return -1;
0082   }
0083 
0084   int mode = std::atoi(argv[1]);
0085   const char* infile = argv[2];
0086   std::string histfile(argv[3]);
0087   int flag = std::atoi(argv[4]);
0088   const char* dirname = argv[5];
0089   std::string prefix(argv[6]);
0090   int pucorr = std::atoi(argv[7]);
0091   int truncate = std::atoi(argv[8]);
0092   Long64_t nmax = static_cast<Long64_t>(std::atoi(argv[9]));
0093 
0094   if (mode == 0) {
0095     // CalibMonitor
0096     bool datamc = (argc > 10) ? (std::atoi(argv[10]) > 0) : true;
0097     int numb = (argc > 11) ? std::atoi(argv[11]) : 50;
0098     bool usegen = (argc > 12) ? (std::atoi(argv[12]) > 0) : false;
0099     double scale = (argc > 13) ? std::atof(argv[13]) : 1.0;
0100     int usescale = (argc > 14) ? std::atoi(argv[14]) : 0;
0101     int etalo = (argc > 15) ? std::atoi(argv[15]) : 0;
0102     int etahi = (argc > 16) ? std::atoi(argv[16]) : 30;
0103     int runlo = (argc > 17) ? std::atoi(argv[17]) : 0;
0104     int runhi = (argc > 18) ? std::atoi(argv[18]) : 99999999;
0105     int phimin = (argc > 19) ? std::atoi(argv[19]) : 1;
0106     int phimax = (argc > 20) ? std::atoi(argv[20]) : 72;
0107     int zside = (argc > 21) ? std::atoi(argv[21]) : 1;
0108     int nvxlo = (argc > 22) ? std::atoi(argv[22]) : 0;
0109     int nvxhi = (argc > 23) ? std::atoi(argv[23]) : 1000;
0110     bool exclude = (argc > 24) ? (std::atoi(argv[24]) > 0) : false;
0111     bool etamax = (argc > 25) ? (std::atoi(argv[25]) > 0) : false;
0112     bool debug = (argc > 26) ? (std::atoi(argv[26]) > 0) : false;
0113     bool append = (argc > 27) ? (std::atoi(argv[27]) > 0) : true;
0114     bool all = (argc > 28) ? (std::atoi(argv[28]) > 0) : true;
0115     const char* corrfile = (argc > 29) ? argv[29] : "";
0116     const char* rcorfile = (argc > 30) ? argv[30] : "";
0117     const char* dupfile = (argc > 31) ? argv[31] : "";
0118     const char* rbxfile = (argc > 32) ? argv[32] : "";
0119     const char* comfile = (argc > 33) ? argv[33] : "";
0120     const char* outfile = (argc > 34) ? argv[34] : "";
0121     if (strcmp(corrfile, "junk.txt") == 0)
0122       corrfile = "";
0123     if (strcmp(rcorfile, "junk.txt") == 0)
0124       rcorfile = "";
0125     if (strcmp(dupfile, "junk.txt") == 0)
0126       dupfile = "";
0127     if (strcmp(comfile, "junk.txt") == 0)
0128       comfile = "";
0129     if (strcmp(rbxfile, "junk.txt") == 0)
0130       rbxfile = "";
0131 
0132     std::cout << "Execute CalibMonitor with infile:" << infile << " dirName: " << dirname << " dupFile: " << dupfile
0133               << " comFile:" << comfile << " outFile:" << outfile << " prefix: " << prefix << " corrFile: " << corrfile
0134               << " rcoFile:" << rcorfile << " puCorr:" << pucorr << " Flag:" << flag << " Numb:" << numb
0135               << " dataMC:" << datamc << " truncate:" << truncate << " useGen:" << usegen << " Scale:" << scale
0136               << " useScale:" << usescale << " etaRange:" << etalo << ":" << etahi << " runRange:" << runlo << ":"
0137               << runhi << " phiRange:" << phimin << ":" << phimax << " zside:" << zside << " nvxRange:" << nvxlo << ":"
0138               << nvxhi << " rbxFile:" << rbxfile << " exclude:" << exclude << " etaMax:" << etamax
0139               << " histFile:" << histfile << " append:" << append << " all:" << all << " nmax:" << nmax << std::endl;
0140 
0141     CalibMonitor c1(infile,
0142                     dirname,
0143                     dupfile,
0144                     comfile,
0145                     outfile,
0146                     prefix,
0147                     corrfile,
0148                     rcorfile,
0149                     pucorr,
0150                     flag,
0151                     numb,
0152                     datamc,
0153                     truncate,
0154                     usegen,
0155                     scale,
0156                     usescale,
0157                     etalo,
0158                     etahi,
0159                     runlo,
0160                     runhi,
0161                     phimin,
0162                     phimax,
0163                     zside,
0164                     nvxlo,
0165                     nvxhi,
0166                     rbxfile,
0167                     exclude,
0168                     etamax);
0169     c1.Loop(nmax, debug);
0170     c1.savePlot(histfile, append, all);
0171   } else if (mode == 1) {
0172     // CalibPlotProperties
0173     bool datamc = (argc > 10) ? (std::atoi(argv[10]) > 0) : true;
0174     bool usegen = (argc > 11) ? (std::atoi(argv[11]) > 0) : false;
0175     double scale = (argc > 12) ? std::atof(argv[12]) : 1.0;
0176     int usescale = (argc > 13) ? std::atoi(argv[13]) : 0;
0177     int etalo = (argc > 14) ? std::atoi(argv[14]) : 0;
0178     int etahi = (argc > 15) ? std::atoi(argv[15]) : 30;
0179     int runlo = (argc > 16) ? std::atoi(argv[16]) : 0;
0180     int runhi = (argc > 17) ? std::atoi(argv[17]) : 99999999;
0181     int phimin = (argc > 18) ? std::atoi(argv[18]) : 1;
0182     int phimax = (argc > 19) ? std::atoi(argv[19]) : 72;
0183     int zside = (argc > 20) ? std::atoi(argv[20]) : 1;
0184     int nvxlo = (argc > 21) ? std::atoi(argv[21]) : 0;
0185     int nvxhi = (argc > 22) ? std::atoi(argv[22]) : 1000;
0186     bool exclude = (argc > 23) ? (std::atoi(argv[23]) > 0) : false;
0187     bool etamax = (argc > 24) ? (std::atoi(argv[24]) > 0) : false;
0188     bool append = (argc > 25) ? (std::atoi(argv[25]) > 0) : true;
0189     bool all = (argc > 26) ? (std::atoi(argv[26]) > 0) : true;
0190     const char* corrfile = (argc > 27) ? argv[27] : "";
0191     const char* rcorfile = (argc > 28) ? argv[28] : "";
0192     const char* dupfile = (argc > 29) ? argv[29] : "";
0193     const char* rbxfile = (argc > 30) ? argv[30] : "";
0194     if (strcmp(corrfile, "junk.txt") == 0)
0195       corrfile = "";
0196     if (strcmp(rcorfile, "junk.txt") == 0)
0197       rcorfile = "";
0198     if (strcmp(dupfile, "junk.txt") == 0)
0199       dupfile = "";
0200     if (strcmp(rbxfile, "junk.txt") == 0)
0201       rbxfile = "";
0202     bool debug(false);
0203 
0204     CalibPlotProperties c1(infile,
0205                            dirname,
0206                            dupfile,
0207                            prefix,
0208                            corrfile,
0209                            rcorfile,
0210                            pucorr,
0211                            flag,
0212                            datamc,
0213                            truncate,
0214                            usegen,
0215                            scale,
0216                            usescale,
0217                            etalo,
0218                            etahi,
0219                            runlo,
0220                            runhi,
0221                            phimin,
0222                            phimax,
0223                            zside,
0224                            nvxlo,
0225                            nvxhi,
0226                            rbxfile,
0227                            exclude,
0228                            etamax);
0229     c1.Loop(nmax);
0230     c1.savePlot(histfile, append, all, debug);
0231   } else if (mode == 2) {
0232     // CalibTree
0233     int maxIter = (argc > 10) ? std::atoi(argv[10]) : 30;
0234     const char* corrfile = (argc > 11) ? argv[11] : "";
0235     int applyl1 = (argc > 12) ? std::atoi(argv[12]) : 1;
0236     double l1cut = (argc > 13) ? std::atof(argv[13]) : 0.5;
0237     bool useiter = (argc > 14) ? (std::atoi(argv[14]) > 0) : true;
0238     bool useweight = (argc > 15) ? (std::atoi(argv[15]) > 0) : true;
0239     bool usemean = (argc > 16) ? (std::atoi(argv[16]) > 0) : false;
0240     int nmin = (argc > 17) ? std::atoi(argv[17]) : 0;
0241     bool inverse = (argc > 18) ? (std::atoi(argv[18]) > 0) : true;
0242     double ratmin = (argc > 19) ? std::atof(argv[19]) : 0.25;
0243     double ratmax = (argc > 20) ? std::atof(argv[20]) : 3.0;
0244     int ietamax = (argc > 21) ? std::atoi(argv[21]) : 25;
0245     int ietatrack = (argc > 22) ? std::atoi(argv[22]) : -1;
0246     int sysmode = (argc > 23) ? std::atoi(argv[23]) : -1;
0247     int rcorform = (argc > 24) ? std::atoi(argv[24]) : 0;
0248     bool usegen = (argc > 25) ? (std::atoi(argv[25]) > 0) : false;
0249     int runlo = (argc > 26) ? std::atoi(argv[26]) : 0;
0250     int runhi = (argc > 27) ? std::atoi(argv[27]) : 99999999;
0251     int phimin = (argc > 28) ? std::atoi(argv[28]) : 1;
0252     int phimax = (argc > 29) ? std::atoi(argv[29]) : 72;
0253     int zside = (argc > 30) ? std::atoi(argv[30]) : 0;
0254     int nvxlo = (argc > 31) ? std::atoi(argv[31]) : 0;
0255     int nvxhi = (argc > 32) ? std::atoi(argv[32]) : 1000;
0256     bool exclude = (argc > 33) ? (std::atoi(argv[33]) > 0) : false;
0257     int higheta = (argc > 34) ? std::atoi(argv[34]) : 1;
0258     double fraction = (argc > 35) ? std::atof(argv[35]) : 1.0;
0259     bool writehisto = (argc > 36) ? (std::atoi(argv[36]) > 0) : false;
0260     double pmin = (argc > 37) ? std::atof(argv[37]) : 40.0;
0261     double pmax = (argc > 38) ? std::atof(argv[38]) : 60.0;
0262     bool debug = (argc > 39) ? (std::atoi(argv[39]) > 0) : false;
0263     const char* rcorfile = (argc > 40) ? argv[40] : "";
0264     const char* dupfile = (argc > 41) ? argv[41] : "";
0265     const char* rbxfile = (argc > 42) ? argv[42] : "";
0266     const char* treename = (argc > 43) ? argv[43] : "CalibTree";
0267     if (strcmp(rcorfile, "junk.txt") == 0)
0268       rcorfile = "";
0269     if (strcmp(dupfile, "junk.txt") == 0)
0270       dupfile = "";
0271     if (strcmp(rbxfile, "junk.txt") == 0)
0272       rbxfile = "";
0273 
0274     char name[500];
0275     sprintf(name, "%s/%s", dirname, treename);
0276     TChain* chain = new TChain(name);
0277     std::cout << "Create a chain for " << name << " from " << infile << std::endl;
0278 
0279     if (!fillChain(chain, infile)) {
0280       std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0281     } else {
0282       std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0283       Long64_t nentryTot = chain->GetEntries();
0284       Long64_t nentries = (fraction > 0.01 && fraction < 0.99) ? (Long64_t)(fraction * nentryTot) : nentryTot;
0285       static const int maxIterMax = 100;
0286       if (maxIter > maxIterMax)
0287         maxIter = maxIterMax;
0288       std::cout << "Tree " << name << " " << chain << " in directory " << dirname << " from file " << infile
0289                 << " with nentries (tracks): " << nentries << std::endl;
0290       unsigned int k(0), kmax(maxIter);
0291       std::cout << "Proceed using CalibTree with dupFile:" << dupfile << " rcorFile:" << rcorfile
0292                 << " trunCate:" << truncate << " useIter:" << useiter << " useMean:" << usemean << " runRange:" << runlo
0293                 << ":" << runhi << " phiRange:" << phimin << ":" << phimax << " zSide:" << zside
0294                 << " nvxRange:" << nvxlo << ":" << nvxhi << " sysMode:" << sysmode << " rbxFile:" << rbxfile
0295                 << " puCorr:" << pucorr << " rcorForm:" << rcorform << " useGen:" << usegen << " exclude:" << exclude
0296                 << " highEta:" << higheta << " pRange:" << pmin << ":" << pmax << std::endl;
0297       CalibTree t(dupfile,
0298                   rcorfile,
0299                   truncate,
0300                   useiter,
0301                   usemean,
0302                   runlo,
0303                   runhi,
0304                   phimin,
0305                   phimax,
0306                   zside,
0307                   nvxlo,
0308                   nvxhi,
0309                   sysmode,
0310                   rbxfile,
0311                   pucorr,
0312                   rcorform,
0313                   usegen,
0314                   exclude,
0315                   higheta,
0316                   pmin,
0317                   pmax,
0318                   chain);
0319       t.h_pbyE = new TH1D("pbyE", "pbyE", 100, -1.0, 9.0);
0320       t.h_Ebyp_bfr = new TProfile("Ebyp_bfr", "Ebyp_bfr", 60, -30, 30, 0, 10);
0321       t.h_Ebyp_aftr = new TProfile("Ebyp_aftr", "Ebyp_aftr", 60, -30, 30, 0, 10);
0322       t.h_cvg = new TH1D("Cvg0", "Convergence", kmax, 0, kmax);
0323       t.h_cvg->SetMarkerStyle(7);
0324       t.h_cvg->SetMarkerSize(5.0);
0325 
0326       TFile* fout = new TFile(histfile.c_str(), "RECREATE");
0327       std::cout << "Output file: " << histfile << " opened in recreate mode" << std::endl;
0328       fout->cd();
0329 
0330       double cvgs[maxIterMax], itrs[maxIterMax];
0331       t.getDetId(fraction, ietatrack, debug, nmax);
0332 
0333       for (; k <= kmax; ++k) {
0334         std::cout << "Calling Loop() " << k << "th time" << std::endl;
0335         double cvg = t.Loop(k,
0336                             fout,
0337                             useweight,
0338                             nmin,
0339                             inverse,
0340                             ratmin,
0341                             ratmax,
0342                             ietamax,
0343                             ietatrack,
0344                             applyl1,
0345                             l1cut,
0346                             k == kmax,
0347                             fraction,
0348                             writehisto,
0349                             debug,
0350                             nmax);
0351         itrs[k] = k;
0352         cvgs[k] = cvg;
0353         if (cvg < 0.00001)
0354           break;
0355       }
0356 
0357       t.writeCorrFactor(corrfile, ietamax);
0358 
0359       fout->cd();
0360       TGraph* g_cvg;
0361       g_cvg = new TGraph(k, itrs, cvgs);
0362       g_cvg->SetMarkerStyle(7);
0363       g_cvg->SetMarkerSize(5.0);
0364       g_cvg->Draw("AP");
0365       g_cvg->Write("Cvg");
0366       std::cout << "Finish looping after " << k << " iterations" << std::endl;
0367       t.makeplots(ratmin, ratmax, ietamax, useweight, fraction, debug, nmax);
0368       fout->Close();
0369     }
0370   } else {
0371     // CalibSplit
0372     double pmin = (argc > 10) ? std::atof(argv[10]) : 40.0;
0373     double pmax = (argc > 11) ? std::atof(argv[11]) : 60.0;
0374     int runMin = (argc > 12) ? std::atoi(argv[12]) : -1;
0375     int runMax = (argc > 13) ? std::atoi(argv[13]) : -1;
0376     bool debug = (argc > 14) ? (std::atoi(argv[14]) > 0) : false;
0377     CalibSplit c1(infile, dirname, histfile, pmin, pmax, runMin, runMax, debug);
0378     c1.Loop(nmax);
0379   }
0380   return 0;
0381 }