macro/loc_lib/loc_lib.h

0001 #ifndef LOC_LIB_H
0002 #define LOC_LIB_H
0003
0004 // local library for jet unfolding
0005
0006 #include "noiBinVec.h"
0007 #include "pAu_bins.h"
0008 #include "noi_fnc.h"
0009
0010
0011 TH1D* norm_binwidths(TH1D* hg) {
0012     auto ax = hg->GetXaxis();
0013     for (int i=1;i<=ax->GetNbins();++i){
0014         double _val = hg->GetBinContent(i);
0015         double _err = hg->GetBinError(i);
0016         double weight = 1./ax->GetBinWidth(i);
0017         hg->SetBinContent(i,_val*weight);
0018         hg->SetBinError  (i,_err*weight);
0019     }
0020     return hg;
0021 }
0022
0023    float calc_dphi(float a, float b) {
0024      float _ = a-b;
0025      while (_ <-M_PI) _ += 2*M_PI;
0026      while (_ > M_PI) _ -= 2*M_PI;
0027      return _;
0028    }
0029
0030 TH1D* rebin_TH1D(TH1D* hg_in, vector<int> i_bins, string _name="", bool delete_hg_in=true, bool norm_by_cnt=true, bool norm_by_width = true) {
0031     // input: bins numbers for left-hand most bin in each new bin. Example: 1, 2, 3, 4, 20, would have 5 bins, with bin 1-1, 2-2, 3-3, 4-20 in each of the new bins
0032     double cnt = hg_in->Integral();
0033     string name = (_name=="") ? noiUniqueName() : _name;
0034     auto ax = hg_in->GetXaxis();
0035     vector<double> edges;
0036     int nbins = i_bins.size()-1;
0037
0038     if (nbins<1) return hg_in;
0039     for (int i=0; i<=nbins; ++i) {
0040         edges.push_back(ax->GetBinLowEdge(i_bins[i]));
0041     }
0042     tuBinVec new_bins { edges };
0043     auto hg_out = new TH1D(name.c_str(), Form("%s;%s;%s",hg_in->GetTitle(),
0044             hg_in->GetXaxis()->GetTitle(), hg_in->GetYaxis()->GetTitle()),
0045             new_bins, new_bins);
0046     for (int i=0; i<nbins; ++i) {
0047         int i0 = i_bins[i];
0048         int i1 = i_bins[i+1];
0049         double sum=0;
0050         for (int k=i0;k<i1;++k) {
0051             sum += hg_in->GetBinContent(k);
0052         }
0053         hg_out->SetBinContent(i+1, sum);
0054         hg_out->SetBinError(i+1, pow(sum,0.5));
0055     }
0056     if (norm_by_width) {
0057         for (int i=1; i<=nbins; ++i) {
0058             auto c = hg_out->GetBinContent(i);
0059             auto e = hg_out->GetBinError  (i);
0060             auto w = hg_out->GetXaxis()->GetBinWidth(i);
0061             hg_out->SetBinContent(i, c/w);
0062             hg_out->SetBinError  (i, e/w);
0063         }
0064     }
0065     if (norm_by_cnt) {
0066         hg_out->Scale(1./cnt);
0067     }
0068     if (delete_hg_in) delete hg_in;
0069     return hg_out;
0070 }
0071
0072 void div_by_W (TH1D* hg, TH1D* w) {
0073     for (int i=1; i<=hg->GetNbinsX(); ++i) {
0074         double W = w->GetBinContent(i);
0075         if (W==0) {
0076             hg->SetBinContent(i, 0.);
0077             hg->SetBinError  (i, 0.);
0078         } else {
0079             hg->SetBinContent(i, hg->GetBinContent(i)/W);
0080             hg->SetBinError  (i, hg->GetBinError  (i)/W);
0081         }
0082     }
0083 }
0084
0085 void set_range(int axis, int bin0, int bin1, vector<THnSparseD*> vec_sparse, bool print=false) {
0086     bool first = true;
0087     for (auto& sparse : vec_sparse) {
0088         TAxis* ax = sparse->GetAxis(axis);
0089         sparse->GetAxis(axis)->SetRange(bin0, bin1);
0090         if (first && print) {
0091             first = false;
0092             cout << Form("Set range(%2i) to %7.2f-%7.2f",axis,ax->GetBinLowEdge(bin0),ax->GetBinUpEdge(bin1)) << endl;
0093         }
0094     }
0095 }
0096
0097 pair<double,double> sparse_integral(THnSparseD* sp) {
0098     auto hg = sp->Projection(0);
0099     double _val, _err;
0100     _val = hg->IntegralAndError(1, hg->GetNbinsX(), _err);
0101     delete hg;
0102     return {_val,_err};
0103 }
0104
0105 TH1D* sparse_proj(THnSparseD* sp, int i_ax, string name="", bool norm_bin_width=false) {
0106     auto hg = (TH1D*) sp->Projection(i_ax);
0107     if (name != "") hg->SetName(name.c_str());
0108     else              hg->SetName(noiUniqueName());
0109     if (norm_bin_width) norm_binwidths(hg);
0110     return hg;
0111 }
0112
0113 void set_track_cuts(vector<THnSparseD*> data) {
0114     set_range(_dca,     1,  5,  data);
0115     set_range(_nhitfit, 20, 48, data);
0116     set_range(_nhitrat, 6,  10, data);
0117 }
0118
0119 TH1D* emb_eff_per_zdcx_binnedby3(THnSparseD* reco, THnSparseD* truth, string name) {
0120     array<int,5> i0_zdcx { 5, 8,  11, 14, 17 }; // binning by 3 for zdcx
0121     array<int,5> i1_zdcx { 7, 10, 13, 16, 19 };
0122
0123     tuBinVec bins_5by3to20 {{ 5000., 5000., 5, 20000 }};
0124     //
0125     /* array<int,7> i0_zdcx { 4, 7, 10, 13, 16, 19, 22 }; */
0126     /* array<int,7> i1_zdcx { 4, 7, 10, 13, 16, 19, 22 }; */
0127
0128     /* tuBinVec bins_5by3to20 {{ 5000., 5000., 5, 20000 }}; */
0129     TH1D* hg_out = new TH1D(name.c_str(), ";ZDCx;Tracking Efficiency", bins_5by3to20, bins_5by3to20 );
0130
0131     for (int i=0; i<5; ++i) {
0132         set_range(_zdcx, i0_zdcx[i], i1_zdcx[i], {reco, truth});
0133         auto h_reco  = (TH1D*) reco->Projection(_pt);  h_reco ->SetName(noiUniqueName());
0134         auto h_truth = (TH1D*) truth->Projection(_pt); h_truth->SetName(noiUniqueName());
0135         double _val, _err;
0136         _val = h_reco->IntegralAndError(3,h_reco->GetNbinsX(),_err);
0137         cout << " i("<<i<<") val("<<_val<<" +/- "<<_err<<")   ----> " << _err/_val << endl;
0138         double norm = h_truth->Integral(3,h_truth->GetNbinsX());
0139         delete h_reco;
0140         delete h_truth;
0141         if (norm == 0) continue;
0142         _val /= norm;
0143         _err /= norm;
0144         hg_out->SetBinContent(i+1,_val);
0145         hg_out->SetBinError  (i+1,_err);
0146     }
0147     return hg_out;
0148 }
0149
0150 TH1D* hg_abs_max(vector<TH1D*> vdat) {
0151     if (vdat.size()==0) return nullptr;
0152     auto hg_new = (TH1D*) vdat[0]->Clone(noiUniqueName());
0153     for (auto bin=1; bin<=hg_new->GetNbinsX(); ++bin) {
0154         double val0 = abs(hg_new->GetBinContent(bin));
0155         for (auto ihg=0; ihg<vdat.size(); ++ihg) {
0156             double val1 = abs(vdat[ihg]->GetBinContent(bin));
0157             if (val1>val0) val0 = val1;
0158         }
0159         hg_new->SetBinContent(bin,val0);
0160         hg_new->SetBinError(bin,0);
0161     }
0162     return hg_new;
0163 }
0164
0165 TH1D* hg_abs_deltarat(TH1D* nom, TH1D* comp) {
0166     auto hg_new = (TH1D*) nom->Clone(noiUniqueName());
0167     hg_new->Reset();
0168     for (int i=1; i<=nom->GetNbinsX(); ++i) {
0169         if (nom->GetBinContent(i) != 0) {
0170             auto v_nom  = nom->GetBinContent(i);
0171             auto v_comp = comp->GetBinContent(i);
0172             hg_new->SetBinContent(i, TMath::Abs(v_comp-v_nom)/v_nom);
0173         }
0174     }
0175     return hg_new;
0176 }
0177
0178 TH1D* hg_abs_delta(TH1D* nom, TH1D* comp) {
0179     auto hg_new = (TH1D*) nom->Clone(noiUniqueName());
0180     hg_new->Reset();
0181     for (int i=1; i<=nom->GetNbinsX(); ++i) {
0182         /* if (nom->GetBinContent(i) != 0) { */
0183             auto v_nom  = nom->GetBinContent(i);
0184             auto v_comp = comp->GetBinContent(i);
0185             hg_new->SetBinContent(i, TMath::Abs(v_comp-v_nom));
0186         /* } */
0187     }
0188     return hg_new;
0189 }
0190
0191 TH1D* hg_const(TH1D* nom, double val, string name) {
0192     auto hg_new = (TH1D*) nom->Clone(name.c_str());
0193     hg_new->Reset();
0194     for (int i=1; i<=hg_new->GetNbinsX(); ++i) {
0195         hg_new->SetBinContent(i, val);
0196     }
0197     return hg_new;
0198 }
0199
0200 TH1D* hg_RSS(vector<TH1D*> vec, string name) {
0201     auto hg_new = (TH1D*) vec[0]->Clone(name.c_str());
0202     hg_new->Reset();
0203     for (auto& hg : vec) {
0204         for (int i=1; i<=hg_new->GetNbinsX(); ++i) {
0205             hg_new->SetBinContent(i, hg_new->GetBinContent(i) + TMath::Sq(hg->GetBinContent(i)));
0206         }
0207     }
0208     for (int i=1; i<=hg_new->GetNbinsX(); ++i) {
0209         hg_new->SetBinContent(i, TMath::Sqrt(hg_new->GetBinContent(i)));
0210     }
0211     return hg_new;
0212 }
0213
0214 TH1D* hg_syserr(TH1D* hg_nom, TH1D* hg_rss, string name) {
0215     auto hg_new = (TH1D*) hg_nom->Clone(name.c_str());
0216     for (int i=1; i<=hg_new->GetNbinsX(); ++i) {
0217         hg_new->SetBinError(i, hg_nom->GetBinContent(i) * hg_rss->GetBinContent(i));
0218         /* hg_new->SetBinError(i, .4); */
0219         /* cout << "i " << i << " -> " << hg_new->GetBinContent(i) << "  " << hg_new->GetBinError(i) << endl; */
0220     }
0221     return hg_new;
0222 }
0223
0224 /* int noi_geant05(int geantid) { */
0225 /*     switch (geantid) { */
0226 /*         case 8: return 0; */
0227 /*         case 9: return 1; */
0228 /*         case 11: return 2; */
0229 /*         case 12: return 3; */
0230 /*         case 14: return 4; */
0231 /*         case 15: return 5; */
0232 /*     } */
0233 /*     return -1; */
0234 /* }; */
0235
0236
0237
0238 // make tuBinVec avialable locally
0239
0240 #endif