Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:21 

0001 #define SDMULT 1
0002 #define CANVASW 800
0003 #define CANVASH 700
0004 
0005 #include <TCanvas.h>
0006 #include <TColor.h>
0007 #include <TCut.h>
0008 #include <TExec.h>
0009 #include <TF1.h>
0010 #include <TFile.h>
0011 #include <TGraphAsymmErrors.h>
0012 #include <TH1F.h>
0013 #include <TH2F.h>
0014 #include <TH3F.h>
0015 #include <TLatex.h>
0016 #include <TLegend.h>
0017 #include <TLine.h>
0018 #include <TMath.h>
0019 #include <TObjString.h>
0020 #include <TROOT.h>
0021 #include <TRandom3.h>
0022 #include <TStyle.h>
0023 #include <TTree.h>
0024 #include <TTreeIndex.h>
0025 
0026 #include <fstream>
0027 #include <iostream>
0028 #include <sstream>
0029 #include <string>
0030 #include <vector>
0031 
0032 #include "acceptance.h"
0033 
0034 #include "/sphenix/user/hjheng/TrackletAna/analysis/plot/sPHENIXStyle/sPhenixStyle.C"
0035 
0036 using namespace std;
0037 
0038 void calccorr(const TString infilename,              //
0039               int CentLow = -1,                      //
0040               int CentHigh = 10,                     //
0041               float PvzMin = -30.,                   //
0042               float PvzMax = -10.,                   //
0043               bool applyc = false,                   //
0044               bool applyg = false,                   //
0045               bool applym = false,                   //
0046               const TString estag = "null",          //
0047               TString aselstr = "",                  //
0048               const TString corfiletag = "baseline", //
0049               const TString outfiletag = "baseline", //
0050               bool debug = true                      //
0051 )
0052 {
0053     SetsPhenixStyle();
0054     gStyle->SetOptTitle(0);
0055     gStyle->SetOptStat(0);
0056     gStyle->SetPalette(kBird);
0057 
0058     TString outbasedir = TString::Format("./plot/corrections/%s", outfiletag.Data());
0059     TString corfiledir = TString::Format("./plot/corrections/%s", corfiletag.Data());
0060     TString zvtxmin_str = TString::Format("%.1f", PvzMin).ReplaceAll(".", "p").ReplaceAll("-", "m");
0061     TString zvtxmax_str = TString::Format("%.1f", PvzMax).ReplaceAll(".", "p").ReplaceAll("-", "m");
0062     TString asel_tag = (aselstr.EqualTo("") || aselstr.EqualTo("none")) ? "noasel" : aselstr.ReplaceAll(" ", "_").ReplaceAll("&&", "AND").ReplaceAll("||", "OR").ReplaceAll("(", "").ReplaceAll(")", "").ReplaceAll("==", "eq").ReplaceAll(">", "gt").ReplaceAll("<", "lt").ReplaceAll(">=", "geq").ReplaceAll("<=", "leq").ReplaceAll("!", "not");
0063     TString outdirtag = TString::Format("Centrality%dto%d_Zvtx%sto%s_%s", CentLow, CentHigh, zvtxmin_str.Data(), zvtxmax_str.Data(), asel_tag.Data());
0064 
0065     std::cout << "[INFO] Output directory tag: " << outdirtag << std::endl;
0066 
0067     TString outdir = Form("%s/%s", outbasedir.Data(), outdirtag.Data());
0068     system(Form("mkdir -p %s", outdir.Data()));
0069 
0070     bool IsData = (infilename.Contains("Data")) ? true : false;
0071 
0072     TString legtitle_head = (IsData) ? "Data" : "Simulation, closure test";
0073     TString legtitle_centrality = TString::Format("Centrality %d-%d", CentLow, CentHigh) + "%";
0074     TString legtitle_zvtx = TString::Format("Z_{vtx} [%.1f, %.1f] cm", PvzMin, PvzMax);
0075     TString legtitle = legtitle_head + TString::Format(": centrality %d-%d", CentLow, CentHigh) + "%, " + legtitle_zvtx;
0076 
0077     TFile *finput = new TFile(infilename, "read");
0078     TTree *tinput = (TTree *)finput->Get("minitree");
0079 
0080     TFile *fcorr = 0;
0081     TFile *faccep = 0;
0082     TFile *fes = 0;
0083 
0084     if (applyc)
0085     {
0086         cout << "[INFO] Applying correction factors" << endl;
0087         fcorr = new TFile(Form("%s/%s/correction_hists.root", corfiledir.Data(), outdirtag.Data()), "READ");
0088         if (!fcorr || fcorr->IsZombie())
0089         {
0090             cout << "[ERROR] No correction file found - exit" << endl;
0091             exit(1);
0092         }
0093 
0094         if (applyg)
0095         {
0096             faccep = new TFile("./plot/corrections/GeoAccepCorr_Run54280.root", "READ");
0097             // faccep = new TFile("./plot/corrections/GeoAccepCorr_Run54280_coarseEta.root", "READ");
0098             if (!faccep || faccep->IsZombie())
0099             {
0100                 cout << "[ERROR] No geometric correction file found - exit" << endl;
0101                 exit(1);
0102             }
0103             cout << "[INFO] Applying geometric correction" << endl;
0104         }
0105     }
0106 
0107     if (applym)
0108     {
0109         cout << "[INFO] Applying external acceptance maps" << endl;
0110     }
0111 
0112     if (estag.EqualTo("")) // for the single-diffractive event fraction
0113     {
0114         // fes = new TFile(Form("output/correction-%s-%i.root", estag, type));
0115         if (!fes || fes->IsZombie())
0116         {
0117             cout << "[ERROR] No external event selection correction file found - exit" << endl;
0118             exit(1);
0119         }
0120         cout << "[INFO] Applying external event selection corrections" << endl;
0121     }
0122 
0123     TCut asel = (aselstr.EqualTo("")) ? "1" : aselstr.Data();                                                                 // Additional event selection (for example, BCO ranges for cross check)
0124     TCut vsel = Form("(PV_z<=%f && PV_z>=%f)", PvzMax, PvzMin);                                                               // Vertex selection
0125     TCut ssel = "1";                                                                                                          // Signal range selection (tracklets)
0126     TCut csel = Form("(MBD_centrality>%f && MBD_centrality<=%f)", static_cast<float>(CentLow), static_cast<float>(CentHigh)); // Centrality selection
0127     TCut osel = (IsData) ? "firedTrig10_MBDSNgeq2==1 && is_min_bias==1" : "is_min_bias==1";                                   // Offline event selection
0128     TCut psel = (IsData) ? "1" : "1";                                                                                         // Single diffractive process (simulation) - not used for now
0129     TCut esel = asel && vsel && csel && osel;
0130     TCut gsel = asel && vsel && csel && psel && ("is_min_bias==1");
0131 
0132     // alpha factor criteria
0133     float alpha_min = 0.0;
0134     float alpha_max = 3.6;
0135     float alpha_pull = 5.0;
0136 
0137     TH1::SetDefaultSumw2();
0138 
0139     /* Setup bins for correction histograms */
0140 #define INCLUDE_VZ_BINS
0141 #define INCLUDE_ETA_BINS
0142 #define INCLUDE_MULT_BINS
0143 #include "bins.h"
0144 
0145     /* Setup histograms */
0146     TH3F *h3WEhadron = new TH3F("h3WEhadron", "h3WEhadron", neta, etab, nmult, multb, nvz, vzb); // Gen-hadron - passing event selection
0147     TH3F *h3WGhadron = new TH3F("h3WGhadron", "h3WGhadron", neta, etab, nmult, multb, nvz, vzb); // Gen-hadron - passing gen-level selection
0148     TH3F *h3WEtruth;
0149 
0150     TH3F *h3WEraw = new TH3F("h3WEraw", "h3WEraw", neta, etab, nmult, multb, nvz, vzb);    // raw reco-tracklets, before alpha correction, passing esel (vz cut, trigger)
0151     TH3F *h3WEcorr = new TH3F("h3WEcorr", "h3WEcorr", neta, etab, nmult, multb, nvz, vzb); // after alpha correction; corr = raw * alpha
0152     /* Acceptance & efficiency correction - alpha */
0153     TH3F *h3alpha = new TH3F("h3alpha", "h3alpha", neta, etab, nmult, multb, nvz, vzb); // h3WEhadron (gen-hadron passing evt selection) / h3WEraw
0154     TH3F *h3alphafinal = new TH3F("h3alphafinal", "h3alphafinal", neta, etab, nmult, multb, nvz, vzb);
0155     /* Vertex distribution: reweighting */
0156     TH1F *h1WEvz = new TH1F("h1WEvz", "h1WEvz", nvz, vzb);                           // v_z distribution passing event selection
0157     TH2F *h2WEvzmult = new TH2F("h2WEvzmult", "h2WEvzmult", nvz, vzb, nmult, multb); // v_z v.s tracklet multiplicity distribution passing event selection
0158     /* Acceptance region */
0159     TH2F *h2amapxev_prealpha; // Before alpha correction is applied, for checking
0160     TH2F *h2amapxev = new TH2F("h2amapxev", "h2amapxev", neta, etab, nvz, vzb);
0161     TH3F *h3amapxemv = new TH3F("h3amapxemv", "h3amapxemv", neta, etab, nmult, multb, nvz, vzb);
0162     /* Trigger/event selection efficiency */
0163     TH1F *h1teff = new TH1F("h1teff", "h1teff", nmult, multb);
0164     TH1F *h1WGOXteff = new TH1F("h1WGOXteff", "h1WGOXteff", nmult, multb); // passing (gsel && osel)
0165     TH1F *h1WGXteff = new TH1F("h1WGXteff", "h1WGXteff", nmult, multb);    // passing gsel
0166     /* Single diffractive event fraction */
0167     TH1F *h1sdf = new TH1F("h1sdf", "h1sdf", nmult, multb);
0168     TH1F *h1WEsdf = new TH1F("h1WEsdf", "h1WEsdf", nmult, multb);
0169     TH1F *h1WENGsdf = new TH1F("h1WENGsdf", "h1WENGsdf", nmult, multb);
0170     /* Empty correction */
0171     TH1F *h1empty = new TH1F("h1empty", "h1empty", neta, etab);
0172     /* Alpha factor in 1D and 2D (for checking) */
0173     TH1F *h1alpha[neta][nvz];
0174     TF1 *falpha[neta][nvz];
0175     TH2F *h2alpha_multb[nmult];
0176     for (int i = 0; i < neta; i++)
0177     {
0178         for (int j = 0; j < nvz; j++)
0179         {
0180             h1alpha[i][j] = new TH1F(Form("alpha_etabin%i_vzbin%i", i, j), Form("alpha_etabin%i_vzbin%i", i, j), nmult, multb);
0181             falpha[i][j] = new TF1(Form("falpha_etabin%i_vzbin%i", i, j), "[0]+[1]/(x+[2])+[3]/(x*x)", 1, 12000);
0182         }
0183     }
0184     for (int k = 0; k < nmult; k++)
0185     {
0186         h2alpha_multb[k] = new TH2F(Form("alphafinal2D_multbin%d", k), Form("alphafinal2D_multbin%d", k), neta, etab, nvz, vzb);
0187     }
0188     /*------------------------------------------------------------------------------------------------------*/
0189 
0190     if (applyc)
0191     {
0192         delete h2amapxev;
0193         delete h1teff;
0194         delete h1sdf;
0195         delete h1empty;
0196 
0197         h2amapxev = (TH2F *)fcorr->Get("h2amapxev")->Clone();
0198         h1teff = (TH1F *)fcorr->Get("h1teff")->Clone();
0199         h1sdf = (TH1F *)fcorr->Get("h1sdf")->Clone();
0200         h1empty = (TH1F *)fcorr->Get("h1empty")->Clone();
0201 
0202         for (int i = 0; i < neta; i++)
0203         {
0204             for (int j = 0; j < nvz; j++)
0205             {
0206                 delete h1alpha[i][j];
0207                 delete falpha[i][j];
0208 
0209                 h1alpha[i][j] = (TH1F *)fcorr->Get(Form("alpha_etabin%i_vzbin%i", i, j));
0210                 falpha[i][j] = (TF1 *)fcorr->Get(Form("falpha_etabin%i_vzbin%i", i, j));
0211             }
0212         }
0213     }
0214 
0215     /*------------------------------------------------------------------------------------------------------*/
0216     TH2F *haccepmc = 0;
0217     TH2F *haccepdata = 0;
0218     TH2F *haccepratio = 0;
0219     // TH2F *hgaccep = 0;
0220 
0221     if (applyg)
0222     {
0223         haccepmc = (TH2F *)faccep->Get("hM_sim_coarse");
0224         haccepdata = (TH2F *)faccep->Get("hM_data_coarse");
0225         haccepratio = (TH2F *)faccep->Get("hM_ratio");
0226         // hgaccep = (TH2F *)haccepmc->Clone("hgaccep");
0227         // hgaccep->Divide(haccepdata);
0228     }
0229     /*------------------------------------------------------------------------------------------------------*/
0230     /* nevents: for normalization */
0231     // Number of events passing event selection and gen-level selection
0232     TH1F *h1WEGevent = new TH1F("h1WEGevent", "", nvz, vzb);
0233     int nWEGentry = tinput->Draw("PV_z>>h1WEGevent", "vtxzwei" * (esel && gsel), "goff");
0234     float nWEGevent = h1WEGevent->Integral(0, h1WEGevent->GetNbinsX() + 1);
0235     // Number of events passing gen-level selection
0236     TH1F *h1WGevent = new TH1F("h1WGevent", "", nvz, vzb);
0237     tinput->Draw("PV_z>>h1WGevent", "vtxzwei" * (gsel), "goff");
0238     float nWGevent = h1WGevent->Integral(0, h1WGevent->GetNbinsX() + 1);
0239 
0240     cout << "[INFO] weighted events: " << nWEGevent << ", entries: " << nWEGentry << endl;
0241     cout << "[INFO] weighted gen-level events: " << nWGevent << endl;
0242     if (nWEGentry < 1)
0243     {
0244         cout << "[ERROR] No events selected - stopping" << endl;
0245         return;
0246     }
0247 
0248     /* set acceptance maps */
0249     cout << "[INFO] Setup the acceptance maps" << endl;
0250     tinput->Project("h1WEvz", "PV_z", "vtxzwei" * (esel));
0251     tinput->Project("h2WEvzmult", "NRecotkl_Raw:PV_z", "vtxzwei" * (esel));
0252 
0253     const int *amap = 0;
0254     if (applym)
0255     {
0256         amap = ext_accep_map();
0257     }
0258 
0259     for (int i = 1; i <= neta; i++)
0260     {
0261         for (int j = 1; j <= nvz; j++)
0262         {
0263             if (!applyc || h2amapxev->GetBinContent(i, j) != 0)
0264             {
0265                 if (applym && !amap[(nvz - j) * neta + i - 1])
0266                 {
0267                     h2amapxev->SetBinContent(i, j, 0);
0268                     h2amapxev->SetBinError(i, j, 0);
0269                     continue;
0270                 }
0271 
0272                 h2amapxev->SetBinContent(i, j, h1WEvz->GetBinContent(j)); // vertex distribution
0273                 h2amapxev->SetBinError(i, j, 0);
0274                 for (int k = 1; k <= nmult; k++)
0275                     h3amapxemv->SetBinContent(i, k, j, h2WEvzmult->GetBinContent(j, k));
0276             }
0277         }
0278     }
0279 
0280     h2amapxev_prealpha = (TH2F *)h2amapxev->Clone("h2amapxev_prealpha");
0281 
0282     /* vertex distribution */
0283     h1WEvz->Scale(1. / h1WEvz->GetEntries());
0284     h1WEvz->Fit("gaus");
0285 
0286     /* generator-level hadrons */
0287     tinput->Project("h3WEhadron", "PV_z:NRecotkl_Raw:GenHadron_eta", "vtxzwei" * (esel && "abs(GenHadron_eta)<4"));
0288     tinput->Project("h3WGhadron", "PV_z:NRecotkl_Raw:GenHadron_eta", "vtxzwei" * (gsel && "abs(GenHadron_eta)<4"));
0289     // print out number of entries
0290     cout << "[INFO] Entries in h3WGhadron: " << h3WGhadron->GetEntries() << endl;
0291 
0292     h3WEtruth = (TH3F *)h3WEhadron->Clone("h3WEtruth");
0293 
0294     /* reconstructed tracklets */
0295     tinput->Project("h3WEraw", "PV_z:NRecotkl_Raw:recotklraw_eta", "vtxzwei" * (ssel && esel));
0296 
0297     /* calculate alpha corrections */
0298     if (!applyc)
0299     {
0300         cout << "[INFO] Calculate the alpha correction" << endl;
0301         for (int x = 1; x <= neta; x++)
0302         {
0303             for (int z = 1; z <= nvz; z++)
0304             {
0305                 if (h2amapxev->GetBinContent(x, z) == 0)
0306                 {
0307                     if (debug)
0308                         printf("   # acceptance region: eta: %2i, vz: %2i          is 0\n", x, z);
0309                     continue;
0310                 }
0311 
0312                 int count = 0;
0313 
0314                 for (int y = 1; y <= nmult; y++)
0315                 {
0316                     h1alpha[x - 1][z - 1]->SetBinContent(y, 0);
0317                     h1alpha[x - 1][z - 1]->SetBinError(y, 0);
0318 
0319                     if (h3WEraw->GetBinContent(x, y, z) != 0 && h3WEhadron->GetBinContent(x, y, z) != 0)
0320                     {
0321                         double raw = h3WEraw->GetBinContent(x, y, z);
0322                         double rawerr = h3WEraw->GetBinError(x, y, z);
0323                         double truth = h3WEhadron->GetBinContent(x, y, z);
0324                         double trutherr = h3WEhadron->GetBinError(x, y, z);
0325 
0326                         double alpha = truth / raw;
0327                         double alphaerr = alpha * sqrt(rawerr / raw * rawerr / raw + trutherr / truth * trutherr / truth);
0328                         if (debug)
0329                             printf("   ^ alpha calculation: eta: %2i, vz: %2i, ntl: %2i, alpha: %5.3f [%5.3f], alpha/alphaerr: %5.3f, raw/truth: %9.2f/%9.2f\n", x, z, y, alpha, alphaerr, (alpha / alphaerr), raw, truth);
0330 
0331                         if (alpha > alpha_min && ((alpha / alphaerr > alpha_pull && alpha < 3) || (alpha < 2)))
0332                         {
0333                             h3alpha->SetBinContent(x, y, z, alpha);
0334                             h3alpha->SetBinError(x, y, z, alphaerr);
0335                             h1alpha[x - 1][z - 1]->SetBinContent(y, alpha);
0336                             h1alpha[x - 1][z - 1]->SetBinError(y, alphaerr);
0337                         }
0338                         else
0339                         {
0340                             if (debug)
0341                                 printf("     ! alpha not used\n");
0342                         }
0343 
0344                         ++count;
0345                     }
0346                 }
0347 
0348                 if (CentHigh - CentLow > 95 && count < nmult / 2)
0349                 {
0350                     if (debug)
0351                         printf("   # ! acceptance map:  eta: %2i, vz: %2i          set to 0 (bad acceptance/statistics)\n", x, z);
0352                     h2amapxev->SetBinContent(x, z, 0);
0353                 }
0354             }
0355         }
0356 
0357         for (int j = 0; j < nvz; j++)
0358         {
0359             for (int i = 0; i < neta; i++)
0360             {
0361                 falpha[i][j]->SetParameter(0, 0.84);
0362                 falpha[i][j]->SetParLimits(1, 1, 256);
0363                 falpha[i][j]->SetParLimits(2, 0, 512);
0364                 falpha[i][j]->SetParLimits(3, -64, 128);
0365                 falpha[i][j]->SetLineWidth(1);
0366                 falpha[i][j]->SetLineColor(46);
0367                 h1alpha[i][j]->Fit(falpha[i][j], "M Q", "", 8, 7000);
0368                 h1alpha[i][j]->Fit(falpha[i][j], "M E Q", "", 8, 7000);
0369                 falpha[i][j] = h1alpha[i][j]->GetFunction(Form("falpha_%i_%i", i, j));
0370             }
0371         }
0372 
0373         /* Trigger & offline selection */
0374         cout << "[INFO] Calculate the trigger & offline selection efficiency" << endl;
0375         tinput->Project("h1WGOXteff", "NRecotkl_Raw", "vtxzwei" * (gsel && osel));
0376         tinput->Project("h1WGXteff", "NRecotkl_Raw", "vtxzwei" * (gsel));
0377         std::cout << "h1WGOXteff: " << h1WGOXteff->Integral(-1, -1) << " h1WGXteff: " << h1WGXteff->Integral(-1, -1) << " Trigger&offline selection efficiency = " << h1WGOXteff->Integral(-1, -1) / h1WGXteff->Integral(-1, -1) << std::endl;
0378         h1teff = (TH1F *)h1WGOXteff->Clone("h1teff");
0379         h1teff->Divide(h1WGXteff);
0380 
0381         /* Single diffractive event fraction (complement of gen selection) */
0382         cout << "[INFO] Calculate the single diffractive event fraction" << endl;
0383         tinput->Project("h1WENGsdf", "NRecotkl_Raw", "vtxzwei" * (esel && !gsel));
0384         tinput->Project("h1WEsdf", "NRecotkl_Raw", "vtxzwei" * (esel));
0385         h1sdf = (TH1F *)h1WENGsdf->Clone("h1sdf");
0386         h1sdf->Divide(h1WEsdf);
0387 
0388         /* vertexing efficiency */
0389     }
0390 
0391     /* external single-diffractive fraction */
0392     if (fes)
0393     {
0394         delete h1sdf;
0395         // delete h1teff;
0396         delete h1empty;
0397         h1sdf = (TH1F *)fes->Get("h1sdf")->Clone();
0398         // h1teff = (TH1F *)fes->Get("h1teff")->Clone();
0399         h1empty = (TH1F *)fes->Get("h1empty")->Clone();
0400     }
0401 
0402     /* Apply alpha correction */
0403     cout << "-------------------------------------------------------------" << endl << "[INFO] Apply the alpha correction" << endl;
0404     for (int x = 1; x <= neta; x++)
0405     {
0406         for (int z = 1; z <= nvz; z++)
0407         {
0408             if (h2amapxev->GetBinContent(x, z) == 0)
0409                 continue;
0410 
0411             for (int y = 1; y <= nmult; y++)
0412             {
0413                 double raw = h3WEraw->GetBinContent(x, y, z);
0414                 double rawerr = h3WEraw->GetBinError(x, y, z);
0415 
0416                 double alpha = h1alpha[x - 1][z - 1]->GetBinContent(y);
0417                 double alphaerr = h1alpha[x - 1][z - 1]->GetBinError(y);
0418 
0419                 if (debug)
0420                     printf("   ^ alpha application: eta: %2i, vz: %2i, ntl: %2i, alpha: %5.3f [%5.3f], raw: %9.2f\n", x, z, y, alpha, alphaerr, raw);
0421 
0422                 if (alpha == 0 && falpha[x - 1][z - 1] != 0)
0423                 {
0424                     alpha = falpha[x - 1][z - 1]->Eval((multb[y] + multb[y - 1]) / 2);
0425                     if (debug)
0426                         printf("     # alpha from histogram is 0 -> check fit value: %.3f\n", alpha);
0427                 }
0428 
0429                 if (alpha == 0)
0430                 {
0431                     for (int k = 0; k < y; k++)
0432                     {
0433                         alpha = h1alpha[x - 1][z - 1]->GetBinContent(y - k);
0434                         alphaerr = h1alpha[x - 1][z - 1]->GetBinError(y - k);
0435                         if (debug)
0436                             printf("     # check other bin: %.3f, bin: %2i\n", alpha, y - k);
0437                         if (alpha != 0)
0438                             break;
0439                     }
0440                 }
0441 
0442                 if (alpha <= alpha_min || alpha > alpha_max)
0443                 {
0444                     if (debug)
0445                         printf("     ! invalid value: %.3f, set to 1\n", alpha);
0446                     alpha = 1;
0447                 }
0448 
0449                 if (applyg)
0450                 {
0451                     double gaccepdata = haccepdata->GetBinContent(x, z);
0452                     double gaccepmc = haccepmc->GetBinContent(x, z);
0453                     double gaccepratio = haccepratio->GetBinContent(x, z);
0454                     double gaccepratio_err = haccepratio->GetBinError(x, z);
0455 
0456                     if (gaccepdata && gaccepmc)
0457                     {
0458                         alpha = alpha * gaccepratio;
0459                         // alphaerr = sqrt(pow(alphaerr, 2) + pow(alpha * gaccepratio_err, 2));
0460                         alphaerr = alphaerr * gaccepratio;
0461                         if (debug)
0462                             printf("     & apply geo accep: %.3f", gaccepmc / gaccepdata);
0463                     }
0464                     else
0465                     {
0466                         if (debug)
0467                             printf("     ! geo accep error: eta: %2i, vz: %2i\n", x, z);
0468                     }
0469                 }
0470 
0471                 if (debug)
0472                     printf("     + alpha applied: [ %.3f ]\n", alpha);
0473 
0474                 h3alphafinal->SetBinContent(x, y, z, alpha);
0475                 h3alphafinal->SetBinError(x, y, z, alphaerr);
0476                 h2alpha_multb[y - 1]->SetBinContent(x, z, alpha);
0477                 h2alpha_multb[y - 1]->SetBinError(x, z, alphaerr);
0478 
0479                 double ncorr = raw * alpha;
0480                 // double ncorrerr = rawerr * alpha;
0481                 double ncorrerr = sqrt(pow(rawerr * alpha, 2) + pow(raw * alphaerr, 2)); // statistical (+) systematic uncertainties
0482 
0483                 h3WEcorr->SetBinContent(x, y, z, ncorr);
0484                 h3WEcorr->SetBinError(x, y, z, ncorrerr);
0485             }
0486 
0487             bool valid = false;
0488             for (int y = 1; y <= nmult; y++)
0489             {
0490                 if (h3alphafinal->GetBinContent(x, y, z) != 1)
0491                     valid = true;
0492             }
0493 
0494             if (!valid)
0495             {
0496                 if (debug)
0497                     printf("   # ! acceptance map:  eta: %2i, vz: %2i          set to 0 (invalid alpha)\n", x, z);
0498                 h2amapxev->SetBinContent(x, z, 0);
0499             }
0500         }
0501 
0502         for (int z = 1; z <= nvz; z++)
0503         {
0504             if (h2amapxev->GetBinContent(x, z) == 0)
0505             {
0506                 for (int y = 1; y <= nmult; y++)
0507                 {
0508                     h3alphafinal->SetBinContent(x, y, z, 0);
0509                     h3alphafinal->SetBinError(x, y, z, 0);
0510                     h2alpha_multb[y - 1]->SetBinContent(x, z, 0);
0511                     h2alpha_multb[y - 1]->SetBinError(x, z, 0);
0512 
0513                     h3WEcorr->SetBinContent(x, y, z, 0);
0514                     h3WEcorr->SetBinError(x, y, z, 0);
0515                 }
0516             }
0517         }
0518     }
0519 
0520     /* 2d corrected, raw tracklets - project to eta and v_z*/
0521     TH2D *h2WEcorr = (TH2D *)h3WEcorr->Project3D("zx");
0522     h2WEcorr->SetName("h2WEcorr");
0523     TH2D *h2WEraw = (TH2D *)h3WEraw->Project3D("zx");
0524     h2WEraw->SetName("h2WEraw");
0525 
0526     /* project 1d acceptance */
0527     TH1F *h1accep2xe_original = (TH1F *)h2amapxev->ProjectionX();
0528     h1accep2xe_original->SetName("h1accep2xe_original");
0529     h1accep2xe_original->Scale(1. / h1accep2xe_original->GetMaximum());
0530     TH1F *h1accep2xe = (TH1F *)h2amapxev->ProjectionX();
0531     h1accep2xe->SetName("h1accep2xe");
0532     h1accep2xe->Scale(1. / h1accep2xe->GetMaximum());
0533     // loop over the bins to check the acceptance. If the bin content is less than 0.8 (8%), set the bin content to 0 (exclude from the analysis)
0534     for (int i = 1; i <= h1accep2xe->GetNbinsX(); i++)
0535     {
0536         if (h1accep2xe->GetBinContent(i) < 0.8)
0537         {
0538             h1accep2xe->SetBinContent(i, 0);
0539             h1accep2xe->SetBinError(i, 0);
0540         }
0541     }
0542 
0543     /* truth (hadrons within acceptance) */
0544     TH1F *h1WEtruth = (TH1F *)h3WEtruth->Project3D("x");
0545     h1WEtruth->SetName("h1WEtruth");
0546     h1WEtruth->Scale(1. / nWEGevent, "width");
0547     h1WEtruth->Divide(h1accep2xe);
0548 
0549     /* reconstructed tracklets */
0550     TH1F *h1WEraw = (TH1F *)h3WEraw->Project3D("x");
0551     h1WEraw->SetName("h1WEraw");
0552     h1WEraw->Scale(1. / nWEGevent, "width");
0553     // h1WEraw->Divide(h1accep2xe);
0554 
0555     /* reconstructed tracklets + acceptance correction */
0556     TH1F *h1WErawacc = (TH1F *)h3WEraw->Project3D("x");
0557     h1WErawacc->SetName("h1WErawacc");
0558     h1WErawacc->Scale(1. / nWEGevent, "width");
0559     h1WErawacc->Divide(h1accep2xe);
0560 
0561     /* reconstructed tracklets + acceptance correction + alpha correction */
0562     TH1F *h1WEcorr = (TH1F *)h3WEcorr->Project3D("x");
0563     h1WEcorr->SetName("h1WEcorr");
0564     // print out the bin content of the h1WEcorr histogram and h1accep2xe
0565     if (debug)
0566     {
0567         for (int i = 1; i <= h1WEcorr->GetNbinsX(); i++)
0568         {
0569             printf("   ^ after alpha correction: bin %2i, h1WEcorr: %.3f, h1accep2xe: %.3f\n", i, h1WEcorr->GetBinContent(i), h1accep2xe->GetBinContent(i));
0570         }
0571     }
0572     h1WEcorr->Scale(1. / nWEGevent, "width");
0573     h1WEcorr->Divide(h1accep2xe);
0574 
0575     /* generator-level hadrons */
0576     TH1F *h1WGhadron = (TH1F *)h3WGhadron->Project3D("x");
0577     h1WGhadron->SetName("h1WGhadron");
0578     h1WGhadron->Scale(1. / nWGevent, "width");
0579 
0580     TH2F *h2WGhadron = (TH2F *)h3WGhadron->Project3D("zx");
0581     h2WGhadron->SetName("h2WGhadron");
0582 
0583     TH1F *h1WGhadronxm = (TH1F *)h3WGhadron->Project3D("y");
0584     h1WGhadronxm->SetName("h1WGhadronxm");
0585     h1WGhadronxm->Scale(1. / nWGevent, "width");
0586 
0587     /*------------------------------------------------------------------------------------------------------*/
0588     /* calculate final results */
0589     cout << "[INFO] Calculate final results with trigger and SDF corrections" << endl;
0590     TH2F *h2WEtcorr = new TH2F("h2WEtcorr", "", neta, etab, nmult, multb);
0591     TH2F *h2WEttruth = new TH2F("h2WEttruth", "", neta, etab, nmult, multb);
0592 
0593     for (int x = 1; x <= neta; x++)
0594     {
0595         for (int y = 1; y <= nmult; y++)
0596         {
0597             double sum = 0, sumerr = 0;
0598             double mcsum = 0, mcsumerr = 0;
0599 
0600             for (int z = 1; z <= nvz; z++)
0601             {
0602                 if (h2amapxev->GetBinContent(x, z) != 0)
0603                 {
0604                     sum += h3WEcorr->GetBinContent(x, y, z);
0605                     double err = h3WEcorr->GetBinError(x, y, z);
0606                     sumerr = sqrt(sumerr * sumerr + err * err);
0607 
0608                     mcsum += h3WEtruth->GetBinContent(x, y, z);
0609                     double mcerr = h3WEtruth->GetBinError(x, y, z);
0610                     mcsumerr = sqrt(mcsumerr * mcsumerr + mcerr * mcerr);
0611                 }
0612             }
0613 
0614             h2WEtcorr->SetBinContent(x, y, sum);
0615             h2WEtcorr->SetBinError(x, y, sumerr);
0616             h2WEttruth->SetBinContent(x, y, mcsum);
0617             h2WEttruth->SetBinError(x, y, mcsumerr);
0618             if (debug)
0619                 printf("   ^ before appling trigger correction: eta: %2i, mult: %2i, h2WEtcorr bin content: %.3f +- %.3f\n", x, y, h2WEtcorr->GetBinContent(x, y), h2WEtcorr->GetBinError(x, y));
0620         }
0621     }
0622 
0623     /* Apply trigger, sd fraction corrections */
0624     for (int y = 1; y <= nmult; y++)
0625     {
0626         double sdfrac = h1sdf->GetBinContent(y);
0627         double trigeff = h1teff->GetBinContent(y);
0628 
0629         double totalc = (1 - sdfrac * SDMULT) / trigeff; 
0630 
0631         for (int x = 1; x <= neta; x++)
0632         {
0633             if (trigeff != 0)
0634             {
0635                 h2WEtcorr->SetBinContent(x, y, h2WEtcorr->GetBinContent(x, y) * totalc);
0636                 h2WEtcorr->SetBinError(x, y, h2WEtcorr->GetBinError(x, y) * totalc);
0637                 if (debug)
0638                     printf("   ^ apply trigger correction: eta: %2i, mult: %2i, trigeff: %.3f, sdfrac: %.3f, totalc: %.3f, h2WEtcorr bin content: %.3f +- %.3f\n", x, y, trigeff, sdfrac, totalc, h2WEtcorr->GetBinContent(x, y), h2WEtcorr->GetBinError(x, y));
0639                 h2WEttruth->SetBinContent(x, y, h2WEttruth->GetBinContent(x, y) * totalc);
0640                 h2WEttruth->SetBinError(x, y, h2WEttruth->GetBinError(x, y) * totalc);
0641             }
0642 
0643             for (int z = 1; z <= nvz; z++)
0644             {
0645                 if (h2amapxev->GetBinContent(x, z) != 0)
0646                 {
0647                     if (trigeff != 0)
0648                         h3amapxemv->SetBinContent(x, y, z, h3amapxemv->GetBinContent(x, y, z) * totalc);
0649                 }
0650                 else
0651                 {
0652                     h3amapxemv->SetBinContent(x, y, z, 0);
0653                 }
0654             }
0655         }
0656     }
0657 
0658     /* project 2d acceptances */
0659     TH1F *h1accep3xe = (TH1F *)h3amapxemv->Project3D("x");
0660     h1accep3xe->SetName("h1accep3xe");
0661     TH1F *h1accep3xe_norm = (TH1F *)h3amapxemv->Project3D("x");
0662     h1accep3xe_norm->SetName("h1accep3xe_norm");
0663     h1accep3xe_norm->Scale(1. / h1accep3xe_norm->GetMaximum());
0664     for (int i = 1; i <= h1accep3xe_norm->GetNbinsX(); i++)
0665     {
0666         if (h1accep3xe_norm->GetBinContent(i) < 0.8)
0667         {
0668             h1accep3xe->SetBinContent(i, 0);
0669             h1accep3xe->SetBinError(i, 0);
0670         }
0671     }
0672 
0673     TH1F *h1accep3xm = (TH1F *)h3amapxemv->Project3D("y");
0674     h1accep3xm->SetName("h1accep3xm");
0675 
0676     TH1F *h1WEtcorr = (TH1F *)h2WEtcorr->ProjectionX();
0677     h1WEtcorr->SetName("h1WEtcorr");
0678     // print out the bin content of the h1WEtcorr histogram and h1accep3xe
0679     if (debug)
0680     {
0681         for (int i = 1; i <= h1WEtcorr->GetNbinsX(); i++)
0682         {
0683             printf("   ^ after trigger correction: bin %2i, h1WEtcorr: %.3f, h1accep3xe: %.3f\n", i, h1WEtcorr->GetBinContent(i), h1accep3xe->GetBinContent(i));
0684         }
0685     }
0686     h1WEtcorr->Scale(1., "width");
0687     h1WEtcorr->Divide(h1accep3xe);
0688 
0689     /* Calculate/apply empty correction */
0690     if (!applyc && !fes)
0691     {
0692         cout << "[INFO] Calculate and apply the empty correction" << endl;
0693         h1empty = (TH1F *)h1WGhadron->Clone("h1empty");
0694         h1empty->Divide(h1WEtcorr);
0695     }
0696 
0697     // The eta bin with an empty correction larger than 1.2 or less than 0.8 should be excluded
0698     for (int i = 1; i <= h1empty->GetNbinsX(); i++)
0699     {
0700         if (h1empty->GetBinContent(i) > 1.2 || h1empty->GetBinContent(i) < 0.8)
0701         {
0702             h1empty->SetBinContent(i, 0);
0703             h1empty->SetBinError(i, 0);
0704         }
0705     }
0706 
0707     cout << "[INFO] Apply the empty correction" << endl;
0708     TH1F *h1WEprefinal = (TH1F *)h1WEtcorr->Clone("h1WEprefinal");
0709     h1WEprefinal->Multiply(h1empty);
0710 
0711     TH1F *h1WEfinal = (TH1F *)h1WEprefinal->Clone("h1WEfinal");
0712 
0713     /*------------------------------------------------------------------------------------------------------*/
0714     /* Make diagonostic plots */
0715     cout << "[INFO] Making checking plots: intermediate steps" << endl;
0716     TGraphAsymmErrors *trigeff = new TGraphAsymmErrors();
0717     trigeff->Divide(h1WGOXteff, h1WGXteff, "cl=0.683 b(1,1) mode");
0718     TGraphAsymmErrors *sdfrac = new TGraphAsymmErrors();
0719     sdfrac->Divide(h1WENGsdf, h1WEsdf, "cl=0.683 b(1,1) mode");
0720 
0721     TCanvas *ccheck = new TCanvas("ccheck", "ccheck", 600, 600);
0722     gPad->SetTopMargin(0.1);
0723     gPad->SetRightMargin(0.15);
0724     gPad->SetLeftMargin(0.15);
0725     TLatex *t = new TLatex();
0726     t->SetTextAlign(23);
0727 
0728     /* Vertex distribution (passing esel)*/
0729     ccheck->cd();
0730     gPad->SetRightMargin(0.09);
0731     h1WEvz->GetXaxis()->SetTitle("vtx_{Z} [cm]");
0732     h1WEvz->GetYaxis()->SetTitle("Entries");
0733     h1WEvz->GetXaxis()->SetRangeUser(-12, 12);
0734     h1WEvz->GetYaxis()->SetRangeUser(0, h1WEvz->GetMaximum() * 1.2);
0735     h1WEvz->GetYaxis()->SetTitleOffset(1.4);
0736     h1WEvz->SetLineColor(1);
0737     h1WEvz->Draw("histtext");
0738     t->DrawLatexNDC(0.5, 0.97, h1WEvz->GetName());
0739     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h1WEvz->GetName()));
0740     ccheck->Clear();
0741 
0742     /* Vertex distribution (passing esel) in 2D (eta,v_z), before applying alpha */
0743     ccheck->cd();
0744     gPad->SetRightMargin(0.17);
0745     h2amapxev_prealpha->GetXaxis()->SetTitle("Tracklet #eta");
0746     h2amapxev_prealpha->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0747     h2amapxev_prealpha->Draw("colz");
0748     t->DrawLatexNDC(0.5, 0.97, h2amapxev_prealpha->GetName());
0749     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h2amapxev_prealpha->GetName()));
0750     ccheck->Clear();
0751 
0752     /* vz v.s tracklet multiplicity distribution passing esel; for 3D acceptance map */
0753     ccheck->cd();
0754     gPad->SetLeftMargin(0.18);
0755     gPad->SetLogy(0);
0756     h2WEvzmult->GetYaxis()->SetMoreLogLabels();
0757     h2WEvzmult->GetXaxis()->SetTitle("vtx_{Z} [cm]");
0758     h2WEvzmult->GetYaxis()->SetTitle("Multiplicity");
0759     h2WEvzmult->GetYaxis()->SetTitleOffset(1.8);
0760     h2WEvzmult->Draw("colz");
0761     t->DrawLatexNDC(0.5, 0.97, h2WEvzmult->GetName());
0762     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h2WEvzmult->GetName()));
0763     ccheck->Clear();
0764 
0765     /* Trigger efficiency */
0766     ccheck->cd();
0767     gPad->SetRightMargin(0.1);
0768     gPad->SetLogx(0);
0769     gPad->SetLogy(0);
0770     // trigeff->GetXaxis()->SetMoreLogLabels();
0771     trigeff->GetXaxis()->SetMaxDigits(2);
0772     trigeff->GetXaxis()->SetTitle("Multiplicity");
0773     trigeff->GetYaxis()->SetTitle("Efficiency");
0774     trigeff->GetYaxis()->SetTitleOffset(1.4);
0775     trigeff->SetMinimum(0);
0776     trigeff->SetMaximum(1.2);
0777     trigeff->SetMarkerStyle(20);
0778     trigeff->SetMarkerColor(1);
0779     trigeff->SetLineColor(1);
0780     trigeff->Draw("ALPE");
0781     t->DrawLatexNDC(0.5, 0.97, "Offline event selection efficiency");
0782     ccheck->SaveAs(Form("%s/TriggerEfficiency.pdf", outdir.Data(), trigeff->GetName()));
0783     ccheck->Clear();
0784 
0785     /* Single-diffractive fraction */
0786     ccheck->cd();
0787     gPad->SetRightMargin(0.1);
0788     gPad->SetLogx(0);
0789     gPad->SetLogy(0);
0790     // sdfrac->GetXaxis()->SetMoreLogLabels();
0791     sdfrac->GetXaxis()->SetMaxDigits(2);
0792     sdfrac->GetXaxis()->SetTitle("Multiplicity");
0793     sdfrac->GetYaxis()->SetTitle("SD event fraction");
0794     sdfrac->GetYaxis()->SetTitleOffset(1.4);
0795     sdfrac->SetMinimum(0);
0796     sdfrac->SetMaximum(1.2);
0797     sdfrac->SetMarkerStyle(20);
0798     sdfrac->SetMarkerColor(1);
0799     sdfrac->SetLineColor(1);
0800     sdfrac->Draw("ALPE");
0801     t->DrawLatexNDC(0.5, 0.97, "SD event fraction");
0802     ccheck->SaveAs(Form("%s/SingleDiffractiveFraction.pdf", outdir.Data(), sdfrac->GetName()));
0803     ccheck->Clear();
0804 
0805     /* Acceptance, after applying alpha */
0806     ccheck->cd();
0807     gPad->SetRightMargin(0.18);
0808     gPad->SetLogx(0);
0809     gPad->SetLogy(0);
0810     h2amapxev->GetXaxis()->SetTitle("Tracklet #eta");
0811     h2amapxev->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0812     h2amapxev->Draw("colz");
0813     t->DrawLatexNDC(0.5, 0.97, Form("%s_postalpha", h2amapxev->GetName()));
0814     ccheck->SaveAs(Form("%s/%s_postalpha.pdf", outdir.Data(), h2amapxev->GetName()));
0815     ccheck->Clear();
0816 
0817     /* Project the 2D acceptance to 1D*/
0818     ccheck->cd();
0819     gPad->SetRightMargin(0.09);
0820     // gStyle->SetPaintTextFormat("1.3f");
0821     h1accep2xe_original->GetXaxis()->SetTitle("Tracklet #eta");
0822     h1accep2xe_original->GetYaxis()->SetTitle("A.U");
0823     h1accep2xe_original->GetYaxis()->SetRangeUser(0, h1accep2xe_original->GetMaximum() * 1.7);
0824     h1accep2xe_original->SetLineColor(1);
0825     h1accep2xe_original->Draw("histtext");
0826     h1accep2xe->SetLineColor(2);
0827     h1accep2xe->Draw("hist same");
0828     t->DrawLatexNDC(0.5, 0.97, legtitle_centrality.Data());
0829     TLegend *leg = new TLegend(0.3, 0.77, 0.85, 0.87);
0830     leg->AddEntry(h1accep2xe_original, "Acceptance correction (Original)", "l");
0831     leg->AddEntry(h1accep2xe, "Acceptance correction (>0.8)", "l");
0832     leg->SetTextSize(0.035);
0833     leg->SetBorderSize(0);
0834     leg->SetFillStyle(0);
0835     leg->Draw();
0836     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h1accep2xe->GetName()));
0837     ccheck->Clear();
0838 
0839     ccheck->cd();
0840     gPad->SetRightMargin(0.09);
0841     gPad->SetLogx(0);
0842     gPad->SetLogy(0);
0843     h1accep3xe->GetXaxis()->SetTitle("Tracklet #eta");
0844     h1accep3xe->GetYaxis()->SetTitle("A.U");
0845     h1accep3xe->GetYaxis()->SetRangeUser(0, h1accep3xe->GetMaximum() * 1.2);
0846     h1accep3xe->SetLineColor(1);
0847     h1accep3xe->Draw("histtext");
0848     t->DrawLatexNDC(0.5, 0.97, legtitle_centrality.Data());
0849     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h1accep3xe->GetName()));
0850     ccheck->Clear();
0851 
0852     /* Raw tracklets, eta v.s vz*/
0853     ccheck->cd();
0854     gPad->SetRightMargin(0.18);
0855     gPad->SetLeftMargin(0.18);
0856     gPad->SetLogx(0);
0857     gPad->SetLogy(0);
0858     h2WEraw->GetXaxis()->SetTitle("Tracklet #eta");
0859     h2WEraw->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0860     h2WEraw->GetYaxis()->SetTitleOffset(1.8);
0861     h2WEraw->Draw("colz");
0862     t->DrawLatexNDC(0.5, 0.97, "Reco-tracklets before corrections (h2WEraw)");
0863     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h2WEraw->GetName()));
0864     ccheck->Clear();
0865 
0866     /* Gen hadron; eta v.s vz */
0867     ccheck->cd();
0868     gPad->SetRightMargin(0.18);
0869     gPad->SetLeftMargin(0.18);
0870     gPad->SetLogx(0);
0871     gPad->SetLogy(0);
0872     h2WGhadron->GetXaxis()->SetTitle("#eta");
0873     h2WGhadron->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0874     h2WGhadron->GetYaxis()->SetTitleOffset(1.8);
0875     h2WGhadron->Draw("colz");
0876     t->DrawLatexNDC(0.5, 0.97, "Generated hadrons (h2WGhadron)");
0877     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h2WGhadron->GetName()));
0878 
0879     /* Reco-tracklet after alpha, trigger, and SDF correction in (eta, tracklet multiplicity) */
0880     ccheck->cd();
0881     gPad->SetRightMargin(0.18);
0882     gPad->SetLeftMargin(0.18);
0883     gPad->SetLogx(0);
0884     gPad->SetLogy(0);
0885     h2WEtcorr->GetYaxis()->SetMoreLogLabels();
0886     h2WEtcorr->GetXaxis()->SetTitle("Tracklet #eta");
0887     h2WEtcorr->GetYaxis()->SetTitle("Multiplicity");
0888     h2WEtcorr->GetYaxis()->SetTitleOffset(1.8);
0889     h2WEtcorr->Draw("colz");
0890     t->DrawLatexNDC(0.5, 0.97, h2WEtcorr->GetName());
0891     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h2WEtcorr->GetName()));
0892     ccheck->Clear();
0893 
0894     /* Empty event correction */
0895     ccheck->cd();
0896     gPad->SetRightMargin(0.09);
0897     gPad->SetLogx(0);
0898     gPad->SetLogy(0);
0899     h1empty->GetXaxis()->SetTitle("Tracklet #eta");
0900     h1empty->GetYaxis()->SetTitle("Ratio");
0901     h1empty->GetYaxis()->SetRangeUser(0, h1empty->GetMaximum() * 1.2);
0902     h1empty->SetLineColor(1);
0903     h1empty->Draw("PE1");
0904     t->DrawLatexNDC(0.5, 0.97, "Empty correction");
0905     ccheck->SaveAs(Form("%s/%s.pdf", outdir.Data(), h1empty->GetName()));
0906     ccheck->Clear();
0907 
0908     /* Draw 1D alpha and fits */
0909     // cout << "[INFO] Making checking plots: 1D alpha and fits" << endl;
0910     // TLatex *t1 = new TLatex();
0911     // t1->SetTextAlign(23);
0912     // TCanvas *cfalphavz = new TCanvas("cfalphavz", "", 2000, 1600);
0913     // cfalphavz->Divide(5, 4);
0914 
0915     // for (int x = 1; x <= neta; x++)
0916     // {
0917     //     cfalphavz->Clear("d");
0918     //     for (int z = 1; z <= nvz; z++)
0919     //     {
0920     //         cfalphavz->cd(z);
0921     //         gPad->SetLogx(0);
0922     //         gPad->SetTickx();
0923     //         gPad->SetTicky();
0924     //         h1alpha[x - 1][z - 1]->GetYaxis()->SetRangeUser(0, h1alpha[x - 1][z - 1]->GetMaximum() * 2.5);
0925     //         h1alpha[x - 1][z - 1]->SetMarkerStyle(20);
0926     //         h1alpha[x - 1][z - 1]->SetMarkerSize(1);
0927     //         h1alpha[x - 1][z - 1]->SetMarkerColor(1);
0928     //         h1alpha[x - 1][z - 1]->SetLineColor(1);
0929     //         h1alpha[x - 1][z - 1]->Draw();
0930 
0931     //         TLine *l = new TLine(multb[0], 1, multb[nmult], 1);
0932     //         l->SetLineColorAlpha(38, 0.7);
0933     //         l->Draw("same");
0934 
0935     //         t1->DrawLatexNDC(0.5, 1.0, Form("%.1f < #eta < %.1f, %.1f < vtx_{Z} < %.1f", etab[x - 1], etab[x], vzb[z - 1], vzb[z]));
0936     //     }
0937 
0938     //     cfalphavz->SaveAs(Form("%s/alpha1Dfit-etabin%i.pdf", outdir.Data(), x - 1));
0939     // }
0940 
0941     // TCanvas *cfalphaeta = new TCanvas("cfalphaeta", "", 2400, 2400);
0942     // cfalphaeta->Divide(6, 6);
0943 
0944     // for (int z = 1; z <= nvz; z++)
0945     // {
0946     //     cfalphaeta->Clear("d");
0947     //     for (int x = 1; x <= neta; x++)
0948     //     {
0949     //         cfalphaeta->cd(x);
0950     //         gPad->SetLogx(0);
0951     //         gPad->SetTickx();
0952     //         gPad->SetTicky();
0953     //         h1alpha[x - 1][z - 1]->SetMarkerStyle(20);
0954     //         h1alpha[x - 1][z - 1]->SetMarkerSize(1);
0955     //         h1alpha[x - 1][z - 1]->SetMarkerColor(1);
0956     //         h1alpha[x - 1][z - 1]->SetLineColor(1);
0957     //         h1alpha[x - 1][z - 1]->Draw();
0958 
0959     //         TLine *l = new TLine(multb[0], 1, multb[nmult], 1);
0960     //         l->SetLineColorAlpha(38, 0.7);
0961     //         l->Draw("same");
0962 
0963     //         t1->DrawLatexNDC(0.5, 1.0, Form("%.1f < #eta < %.1f, %.1f < vtx_{Z} < %.1f", etab[x - 1], etab[x], vzb[z - 1], vzb[z]));
0964     //     }
0965 
0966     //     cfalphaeta->SaveAs(Form("%s/alpha1Dfit-vzbin%i.pdf", outdir.Data(), z));
0967     // }
0968 
0969     /* Draw 2D alpha (eta, vz), inclusive tracklet multiplicity */
0970     cout << "[INFO] Making checking plots: 2D alpha (eta, vz), inclusive tracklet multiplicity" << endl;
0971     TCanvas *calpha = new TCanvas("calpha", "calpha", CANVASW, CANVASH);
0972     gPad->SetRightMargin(0.18);
0973     gPad->SetLeftMargin(0.13);
0974     calpha->cd();
0975     TH2D *h2alphafinal = (TH2D *)h2WEcorr->Clone("h2alphafinal");
0976     h2alphafinal->Divide(h2WEraw);
0977     h2alphafinal->GetXaxis()->SetTitle("#eta");
0978     h2alphafinal->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0979     h2alphafinal->GetYaxis()->SetTitleOffset(1.3);
0980     // h2alphafinal->GetYaxis()->SetRangeUser(-12, 12);
0981     h2alphafinal->GetZaxis()->SetTitle("#alpha (#eta, vtx_{Z})");
0982     h2alphafinal->GetZaxis()->SetTitleOffset(1.3);
0983     h2alphafinal->GetZaxis()->SetRangeUser(alpha_min, alpha_max);
0984     gStyle->SetPaintTextFormat("1.3f");
0985     h2alphafinal->SetMarkerSize(0.6);
0986     h2alphafinal->SetContour(1000);
0987     h2alphafinal->Draw("colztext45");
0988     // TLegend *la = new TLegend(0.18, 0.18, 0.4, 0.3);
0989     // la->SetHeader(legtitle_centrality.Data(), "l");
0990     // la->Draw();
0991     TText *tt = new TText();
0992     // right and bottom adjusted
0993     tt->SetTextAlign(31);
0994     tt->SetTextSize(0.04);
0995     tt->DrawTextNDC(1 - gPad->GetRightMargin(), (1 - gPad->GetTopMargin()) + 0.01, legtitle_centrality.Data());
0996     calpha->SaveAs(Form("%s/alpha2D_inclTklMult.pdf", outdir.Data()));
0997 
0998     for (int k = 0; k < nmult; k++)
0999     {
1000         calpha->Clear();
1001         calpha->cd();
1002         h2alpha_multb[k]->GetXaxis()->SetTitle("#eta");
1003         h2alpha_multb[k]->GetYaxis()->SetTitle("vtx_{Z} [cm]");
1004         h2alpha_multb[k]->GetYaxis()->SetTitleOffset(1.3);
1005         h2alpha_multb[k]->GetZaxis()->SetTitle("#alpha (#eta, vtx_{Z})");
1006         h2alpha_multb[k]->GetZaxis()->SetTitleOffset(1.3);
1007         h2alpha_multb[k]->GetZaxis()->SetRangeUser(alpha_min, alpha_max);
1008         gStyle->SetPaintTextFormat("1.3f");
1009         h2alpha_multb[k]->SetMarkerSize(0.6);
1010         h2alpha_multb[k]->SetContour(1000);
1011         h2alpha_multb[k]->Draw("colztext45");
1012         // la->Draw();
1013         tt->DrawTextNDC(1 - gPad->GetRightMargin(), (1 - gPad->GetTopMargin()) + 0.01, legtitle_centrality.Data());
1014         calpha->SaveAs(Form("%s/alpha2D_TklMultb%d.pdf", outdir.Data(), k));
1015     }
1016 
1017     /* Analysis stages - Closure test*/
1018     // vector<TH1F *> vechist = {h1WGhadron, h1WEraw, h1WErawacc, h1WEcorr, h1WEtcorr, h1WEfinal};
1019     // vector<TH1F *> vechist = {h1WGhadron, h1WEraw, h1WErawacc, h1WEtcorr, h1WEfinal};
1020     vector<TH1F *> vechist = {h1WGhadron, h1WEraw, h1WErawacc, h1WEfinal};
1021     // vector<TH1F *> vechist = {h1WGhadron, h1WEraw, h1WEfinal};
1022     // vector<const char *> vcolor{"#20262E", "#7209b7", "#D04848", "#1C82AD", "#1F8A70", "#FC7300"};
1023     // vector<const char *> vcolor{"#D04848", "#1C82AD", "#1F8A70", "#FC7300", "#20262E"};
1024     vector<const char *> vcolor{"#D04848", "#1C82AD", "#1F8A70", "#20262E"};
1025     // vector<const char *> vcolor{"#D04848", "#1C82AD", "#20262E"};
1026     // vector<int> vlwidth = {3, 3, 3, 3, 2, 2};
1027     // vector<int> vlwidth = {3, 3, 3, 3, 3};
1028     vector<int> vlwidth = {3, 3, 3, 3};
1029     // vector<int> vlwidth = {3, 3, 3};
1030     // vector<int> vstyle = {1, 1, 1, 1, 2, 1};
1031     // vector<int> vstyle = {1, 1, 1, 1, 1};
1032     vector<int> vstyle = {1, 1, 1, 1};
1033     // vector<int> vstyle = {1, 1, 1};
1034     TCanvas *cstage = new TCanvas("cstage", "cstage", CANVASW, CANVASH);
1035     gPad->SetRightMargin(0.05);
1036     gPad->SetTopMargin(0.05);
1037     gPad->SetTickx();
1038     gPad->SetTicky();
1039     // cstage->SetLogy();
1040     TH1F *hframe = new TH1F("hframe", "", 1, etamin, etamax);
1041     float ymax = -1;
1042     for (size_t i = 0; i < vechist.size(); i++)
1043     {
1044         if (vechist[i]->GetMaximum() > ymax)
1045             ymax = vechist[i]->GetMaximum();
1046     }
1047     hframe->GetYaxis()->SetRangeUser(0, ymax * 1.6);
1048     // hframe->GetYaxis()->SetRangeUser(0.1, ymax * 150);
1049     hframe->GetXaxis()->SetTitle("#eta");
1050     hframe->GetYaxis()->SetTitle("dN_{ch}/d#eta");
1051     hframe->GetYaxis()->SetTitleOffset(1.4);
1052     hframe->Draw();
1053     for (size_t i = 0; i < vechist.size(); i++)
1054     {
1055         vechist[i]->SetLineColor(TColor::GetColor(vcolor[i]));
1056         vechist[i]->SetLineWidth(vlwidth[i]);
1057         vechist[i]->SetLineStyle(vstyle[i]);
1058         if (i == 0)
1059             vechist[i]->Draw("hist same");
1060         else
1061         {
1062             vechist[i]->SetFillColorAlpha(TColor::GetColor(vcolor[i]), 0.6);
1063             vechist[i]->SetMarkerStyle(54);
1064             // vechist[i]->SetMarkerSize(0.6);
1065             vechist[i]->SetMarkerColor(TColor::GetColor(vcolor[i]));
1066             vechist[i]->Draw("E5 same");
1067         }
1068     }
1069 
1070     TLegend *l1 = new TLegend(0.18, 0.68, 0.5, 0.92);
1071     l1->SetHeader(legtitle.Data(), "l");
1072     l1->AddEntry(h1WGhadron, "Truth", "l");
1073     l1->AddEntry(h1WEraw, "Reco-tracklets before corrections", "pl");
1074     l1->AddEntry(h1WErawacc, "Corrected for acceptance", "pl");
1075     // l1->AddEntry(h1WEcorr, "Corrected for accep. & #alpha", "pl");
1076     // l1->AddEntry(h1WEtcorr, "Corrected for accep. & #alpha, evt. sel.", "pl");
1077     l1->AddEntry(h1WEfinal, "After corrections", "pl");
1078     l1->SetTextSize(0.03);
1079     l1->SetFillStyle(0);
1080     l1->SetBorderSize(0);
1081     l1->Draw();
1082     cstage->Draw();
1083     cstage->SaveAs(Form("%s/Closure.pdf", outdir.Data()));
1084 
1085     /* Write histograms to file*/
1086     TFile *outf = new TFile(Form("%s/correction_hists.root", outdir.Data()), "recreate");
1087     outf->cd();
1088     h2amapxev->Write();
1089     h1WGhadron->Write();
1090     h1WEraw->Write();
1091     h1WEcorr->Write();
1092     h1WEtcorr->Write();
1093     h1WEfinal->Write(); // Final result, reco-tracklet after alpha, trg&sdf, empty event corrections
1094     trigeff->Write();
1095     sdfrac->Write();
1096     h1teff->Write();
1097     h1WGOXteff->Write();
1098     h1WGXteff->Write();
1099     h1sdf->Write();
1100     h1WEsdf->Write();
1101     h1WENGsdf->Write();
1102     h1empty->Write();
1103     for (int i = 0; i < neta; i++)
1104     {
1105         for (int j = 0; j < nvz; j++)
1106         {
1107             if (falpha[i][j])
1108                 falpha[i][j]->Write();
1109             if (h1alpha[i][j])
1110                 h1alpha[i][j]->Write();
1111         }
1112     }
1113 
1114     for (int k = 0; k < nmult; k++)
1115         h2alpha_multb[k]->Write();
1116 
1117     outf->Write("", TObject::kOverwrite);
1118     outf->Close();
1119 }
1120 
1121 int main(int argc, char *argv[])
1122 {
1123     if (argc != 14)
1124     {
1125         std::cout << "Usage: ./Corrections [infile] [CentLow] [CentHigh] [pvzmin] [pvzmax] [applyc] [applyg] [applym] [estag] [aselstring] [correctionfiletag] [outfilepath] [debug]" << std::endl;
1126         return 0;
1127     }
1128 
1129     for (int i = 0; i < argc; i++)
1130     {
1131         std::cout << "argv[" << i << "] = " << argv[i] << std::endl;
1132     }
1133 
1134     const TString input = TString(argv[1]);
1135     int CentLow = TString(argv[2]).Atoi();
1136     int CentHigh = TString(argv[3]).Atoi();
1137     float PVzMin = TString(argv[4]).Atof();
1138     float PVzMax = TString(argv[5]).Atof();
1139     bool applyc = (TString(argv[6]).Atoi() == 1) ? true : false;
1140     bool applyg = (TString(argv[7]).Atoi() == 1) ? true : false; // geometric difference; applied on data only
1141     bool applym = (TString(argv[8]).Atoi() == 1) ? true : false; // geometric difference map; applied on data only
1142     const TString estag = TString(argv[9]);
1143     const TString asel_string = TString(argv[10]);
1144     const TString correctionfiletag = TString(argv[11]);
1145     const TString outfilepath = TString(argv[12]);
1146     bool debug = (TString(argv[13]).Atoi() == 1) ? true : false;
1147 
1148     gStyle->SetPaintTextFormat();
1149 
1150     const TString aselstring = (asel_string == "null") ? "" : asel_string;
1151     calccorr(input, CentLow, CentHigh, PVzMin, PVzMax, applyc, applyg, applym, estag, aselstring, correctionfiletag, outfilepath, debug);
1152 
1153     return 0;
1154 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!