delphes/scripts/dloop.C

0001
0002 //
0003 // root -l dloop.C
0004 //
0005
0006 #include <dirent.h>
0007
0008 #define _AEROGEL_
0009
0010 void dloop(const char *fname = 0)
0011 {
0012 #ifdef _AEROGEL_
0013   auto delphes = new DelphesConfigRICH("dualRICH_aerogel");
0014   delphes->SetRefractiveIndex(1.0190);
0015   const char *dname = "EVALUATION.A";
0016 #else
0017   auto delphes = new DelphesConfigRICH("dualRICH_c2f6");
0018   delphes->SetRefractiveIndex(1.00076);
0019   const char *dname = "EVALUATION.G";
0020 #endif
0021
0022   // Define particle mass hypotheses in ascending mass order; yes, there is no
0023   // reason to overcomplicate things;
0024   delphes->AddMassHypothesis("pi+");
0025   delphes->AddMassHypothesis("K+");
0026   delphes->AddMassHypothesis("proton");
0027
0028   // Imitate tracker resolution;
0029   delphes->SetAdditionalSmearing(0.50);
0030
0031   // Loop through all .root file in this directory; the files presumably contain a single
0032   // tree with a single 'th' variable;
0033
0034   struct dirent **namelist;
0035   int n = scandir(dname, &namelist, 0, alphasort);
0036
0037     //while((curr_file = readdir(curr_dir))) {
0038   for(unsigned iq=0; iq<n; iq++) {
0039     if (!strcmp(namelist[iq]->d_name, ".") || !strcmp(namelist[iq]->d_name, "..")) continue;
0040     {
0041       unsigned len = strlen(namelist[iq]->d_name);
0042       if (len < 5 || strcmp(namelist[iq]->d_name + len - 5, ".root")) continue;
0043     }
0044
0045     if (fname && strcmp(fname, namelist[iq]->d_name)) continue;
0046
0047     auto ifdata = new TFile((std::string(dname) + "/" + namelist[iq]->d_name).c_str());
0048     if (!ifdata) {
0049       printf("failed to open '%s'\n", namelist[iq]->d_name);
0050       exit(0);
0051     } //if
0052     TTree *it = dynamic_cast<TTree*>(ifdata->Get("t"));
0053     if (!it) {
0054       printf("input file '%s' does not have \"t\" tree\n", namelist[iq]->d_name);
0055       exit(0);
0056     } //if
0057
0058     int pdg;
0059     float emin, emax, pmin, pmax;
0060     const char *format = "drich-data.%d..%f-%f..%f-%f..juggler.evaluation.root";
0061     sscanf(namelist[iq]->d_name, format, &pdg, &emin, &emax, &pmin, &pmax);
0062
0063     double par[100]; memset(par, 0x00, sizeof(par));
0064
0065     const char *var = "1000*th";
0066     TH1D *dth1 = new TH1D("dth1", "dth1", 50, -10.0, 10.0);
0067     it->Project("dth1", var);
0068     TF1 *fq1 = new TF1("fq1", "gaus(0)",    -10.0, 10.0);
0069     par[0] = 100.0; par[1] = 0.0; par[2] = 1.0;
0070     fq1->SetParameters(par);
0071     dth1->Fit("fq1","R");
0072     fq1->GetParameters(par);
0073
0074     double mean2 = par[1], sigma2 = fabs(par[2]), min2 = mean2-5*sigma2, max2 = mean2+5*sigma2;
0075
0076     TH1D *dth2 = new TH1D("dth2", "dth2", 50, min2, max2);
0077     it->Project("dth2", var);
0078     TF1 *fq2 = new TF1("fq2", "gaus(0)", min2, max2);
0079     fq2->SetParameters(par);
0080     dth2->Fit("fq2","R");
0081     fq2->GetParameters(par);
0082
0083     double mean3 = par[1], sigma3 = fabs(par[2]), min3 = mean3-5*sigma3, max3 = mean3+5*sigma3;
0084
0085     TH1D *dth3 = new TH1D("dth3", "dth3", 50, min3, max3);
0086     it->Project("dth3", var);
0087     TF1 *fq3 = new TF1("fq3", "gaus(0)", min3, max3);
0088     fq3->SetParameters(par);
0089     dth3->Fit("fq3","R");
0090
0091     printf("@@@ \"%s\" -> %6.2f\n", namelist[iq]->d_name, fabs(par[2]));
0092
0093     // Shamelessly use the fact that 211/321/2212 are in alphabetic order, and
0094     // alphasort takes care about ordering; then all eta ranges and all momentum
0095     // ranges will be imported in order, and since 211 files are always present
0096     // the whole flakey system will seemingly work;
0097     auto erange = delphes->AddEtaRange    (emin, emax);
0098     auto mrange = erange->GetMomentumRange(pmin, pmax);
0099     // FIXME: should check the return code;
0100     {
0101       bool ret = delphes->StoreSigmaEntry(mrange, pdg, fabs(par[2]));
0102       if (!ret) {
0103     printf("Failed to insert a sigma entry\n");
0104     exit(0);
0105       } //if
0106     }
0107   } //for iq
0108
0109   delphes->AddZeroSigmaEntries();
0110   delphes->Print();
0111   delphes->WriteTcl();
0112
0113   if (!fname) exit(0);
0114 } // dconfig()