Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-07-14 01:08:24

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