analysis_INTT/src/plotRecoVtx.cxx

0001 #include <TCanvas.h>
0002 #include <TCut.h>
0003 #include <TFile.h>
0004 #include <TH1F.h>
0005 #include <TH2F.h>
0006 #include <TMath.h>
0007 #include <TObjString.h>
0008 #include <TRandom3.h>
0009 #include <TTree.h>
0010 #include <TTreeIndex.h>
0011
0012 #include <fstream>
0013 #include <iostream>
0014 #include <sstream>
0015 #include <string>
0016 #include <vector>
0017
0018 #include "Utilities.h"
0019
0020 using namespace std;
0021
0022 int main(int argc, char *argv[])
0023 {
0024     if (argc != 4)
0025     {
0026         std::cout << "Usage: ./plotRecoVtx [isdata] [infilename (ntuple)] [outfilename (histogram)]" << std::endl;
0027         return 0;
0028     }
0029
0030     for (int i = 0; i < argc; i++)
0031     {
0032         std::cout << "argv[" << i << "] = " << argv[i] << std::endl;
0033     }
0034
0035     bool isdata = atoi(argv[1]);
0036     TString infilename = TString(argv[2]);
0037     TString outfilename = TString(argv[3]);
0038
0039     TH2F *hM_INTTVtxZ_MBDVtxZ_Inclusive = new TH2F("hM_INTTVtxZ_MBDVtxZ_Inclusive", "hM_INTTVtxZ_MBDVtxZ_Inclusive", 160, -40, 40, 160, -40, 40);
0040     TH1F *hM_INTTVtxZ_Inclusive = new TH1F("hM_INTTVtxZ_Inclusive", "hM_INTTVtxZ_Inclusive", 160, -40, 40);
0041     TH1F *hM_INTTVtxZ_MBDAsymLe1 = new TH1F("hM_INTTVtxZ_MBDAsymLe1", "hM_INTTVtxZ_MBDAsymLe1", 160, -40, 40);
0042     TH1F *hM_INTTVtxZ_MBDAsymLe1_VtxZm10to10 = new TH1F("hM_INTTVtxZ_MBDAsymLe1_VtxZm10to10", "hM_INTTVtxZ_MBDAsymLe1_VtxZm10to10", 100, -10, 10);
0043     TH1F *hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10 = new TH1F("hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10", "hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10", 100, -10, 10);
0044     // for reweighting
0045     TH1F *hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse = new TH1F("hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse", "hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse", 60, -30, 30);
0046     TH1F *hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse = new TH1F("hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse", "hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse", 20, -10, 10);
0047     TH1F *hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse = new TH1F("hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse", "hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse", 60, -30, 30);
0048     TH1F *hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse = new TH1F("hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse", "hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse", 20, -10, 10);
0049     TH2F *hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive = new TH2F("hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive", "hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive", 160, -40, 40, 200, -1, 1);
0050     TH2F *hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10 = new TH2F("hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10", "hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10", 200, -10, 10, 200, -1, 1);
0051     TH2F *hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1 = new TH2F("hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1", "hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1", 160, -40, 40, 160, -40, 40);
0052     vector<float> centrality_cut = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100};
0053     vector<TH1F *> hM_INTTVtxZ_Centrality;
0054     vector<TH1F *> hM_MBDVtxZ_Centrality;
0055     vector<TH1F *> hM_INTTVtxZ_Centrality_MBDAsymLe1;
0056     vector<TH2F *> hM_INTTVtxZ_MBDVtxZ_Centrality;
0057     vector<TH1F *> hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10;
0058     vector<TH2F *> hM_INTTVtxZ_MBDAsymm_Centrality_MBDAsymLe1_VtxZm10to10;
0059     for (int i = 0; i < centrality_cut.size() - 1; i++)
0060     {
0061         hM_INTTVtxZ_Centrality.push_back(new TH1F(Form("hM_INTTVtxZ_Centrality_%dto%d", (int)centrality_cut[i], (int)centrality_cut[i + 1]), Form("hM_INTTVtxZ_Centrality_%dto%d", (int)centrality_cut[i], (int)centrality_cut[i + 1]), 160, -40, 40));
0062         hM_MBDVtxZ_Centrality.push_back(new TH1F(Form("hM_MBDVtxZ_Centrality_%dto%d", (int)centrality_cut[i], (int)centrality_cut[i + 1]), Form("hM_MBDVtxZ_Centrality_%dto%d", (int)centrality_cut[i], (int)centrality_cut[i + 1]), 160, -40, 40));
0063         hM_INTTVtxZ_Centrality_MBDAsymLe1.push_back(new TH1F(Form("hM_INTTVtxZ_Centrality_%dto%d_MBDAsymLe1", (int)centrality_cut[i], (int)centrality_cut[i + 1]), Form("hM_INTTVtxZ_Centrality_%dto%d_MBDAsymLe1", (int)centrality_cut[i], (int)centrality_cut[i + 1]), 160, -40, 40));
0064         hM_INTTVtxZ_MBDVtxZ_Centrality.push_back(new TH2F(Form("hM_INTTVtxZ_MBDVtxZ_Centrality_%dto%d", (int)centrality_cut[i], (int)centrality_cut[i + 1]), Form("hM_INTTVtxZ_MBDVtxZ_Centrality_%dto%d", (int)centrality_cut[i], (int)centrality_cut[i + 1]), 160, -40, 40, 160, -40, 40));
0065         hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10.push_back(new TH1F(Form("hM_INTTVtxZ_Centrality_%dto%d_MBDAsymLe1_VtxZm10to10", (int)centrality_cut[i], (int)centrality_cut[i + 1]), Form("hM_INTTVtxZ_Centrality_%dto%d_MBDAsymLe1_VtxZm10to10", (int)centrality_cut[i], (int)centrality_cut[i + 1]), 100, -10, 10));
0066         hM_INTTVtxZ_MBDAsymm_Centrality_MBDAsymLe1_VtxZm10to10.push_back(new TH2F(Form("hM_INTTVtxZ_MBDAsymm_Centrality_%dto%d_MBDAsymLe1_VtxZm10to10", (int)centrality_cut[i], (int)centrality_cut[i + 1]), Form("hM_INTTVtxZ_MBDAsymm_Centrality_%dto%d_MBDAsymLe1_VtxZm10to10", (int)centrality_cut[i], (int)centrality_cut[i + 1]), 100, -10, 10, 200, -1, 1));
0067     }
0068
0069     TFile *f = new TFile(infilename, "READ");
0070     TTree *t = (TTree *)f->Get("minitree");
0071     int event;
0072     vector<int> *firedTriggers = 0;
0073     bool is_min_bias;
0074     int NClusLayer1, NTruthVtx;
0075     float TruthPV_x, TruthPV_y, TruthPV_z, Centrality_bimp, Centrality_impactparam, MBD_centrality, mbd_charge_sum, mbd_south_charge_sum, mbd_north_charge_sum, mbd_charge_asymm, mbd_z_vtx, PV_x, PV_y, PV_z;
0076     t->SetBranchAddress("event", &event);
0077     if (isdata)
0078     {
0079         t->SetBranchAddress("firedTriggers", &firedTriggers);
0080     }
0081     t->SetBranchAddress("is_min_bias", &is_min_bias);
0082     t->SetBranchAddress("NClusLayer1", &NClusLayer1);
0083     t->SetBranchAddress("NTruthVtx", &NTruthVtx);
0084     t->SetBranchAddress("TruthPV_x", &TruthPV_x);
0085     t->SetBranchAddress("TruthPV_y", &TruthPV_y);
0086     t->SetBranchAddress("TruthPV_z", &TruthPV_z);
0087     t->SetBranchAddress("Centrality_bimp", &Centrality_bimp);
0088     t->SetBranchAddress("Centrality_impactparam", &Centrality_impactparam);
0089     t->SetBranchAddress("MBD_centrality", &MBD_centrality);
0090     t->SetBranchAddress("mbd_charge_sum", &mbd_charge_sum);
0091     t->SetBranchAddress("mbd_south_charge_sum", &mbd_south_charge_sum);
0092     t->SetBranchAddress("mbd_north_charge_sum", &mbd_north_charge_sum);
0093     t->SetBranchAddress("mbd_charge_asymm", &mbd_charge_asymm);
0094     t->SetBranchAddress("mbd_z_vtx", &mbd_z_vtx);
0095     t->SetBranchAddress("PV_x", &PV_x);
0096     t->SetBranchAddress("PV_y", &PV_y);
0097     t->SetBranchAddress("PV_z", &PV_z);
0098
0099     for (int ev = 0; ev < t->GetEntriesFast(); ev++)
0100     {
0101         t->GetEntry(ev);
0102
0103         // set up selection criteria
0104         bool goodMbdAsymm = (fabs(mbd_charge_asymm) <= 1.0);
0105         bool validMbdVtx = (mbd_z_vtx == mbd_z_vtx);
0106         bool MbdNScharnge = (mbd_north_charge_sum > 0 || mbd_south_charge_sum > 0);
0107         bool MbdZvtxRange_wide = (mbd_z_vtx >= -30 && mbd_z_vtx <= 30);
0108         bool MbdZvtxRange_narrow = (mbd_z_vtx >= -10 && mbd_z_vtx <= 10);
0109         bool firedTrig10 = (isdata) ? (find(firedTriggers->begin(), firedTriggers->end(), 10) != firedTriggers->end()) : true; // MBD N&S >= 2
0110         bool firedTrig12 = (isdata) ? (find(firedTriggers->begin(), firedTriggers->end(), 12) != firedTriggers->end()) : true; // MBD N&S >= 2, vtx < 10 cm
0111         bool firedTrig13 = (isdata) ? (find(firedTriggers->begin(), firedTriggers->end(), 13) != firedTriggers->end()) : true; // MBD N&S >= 2, vtx < 30 cm
0112         bool goodCentrality = (MBD_centrality >= 0 && MBD_centrality <= 70);
0113         // bool goodCentrality = true;
0114         bool InttZvtxRange_wide = (PV_z >= -30 && PV_z <= 30);
0115         bool InttZvtxRange_narrow = (PV_z >= -10 && PV_z <= 10);
0116
0117         // Histograms - Inclusive, wide z vtx range
0118         if (is_min_bias && validMbdVtx && MbdNScharnge && firedTrig10 && MbdZvtxRange_wide)
0119         {
0120             hM_INTTVtxZ_MBDVtxZ_Inclusive->Fill(PV_z, mbd_z_vtx);
0121             hM_INTTVtxZ_Inclusive->Fill(PV_z);
0122
0123             if (goodCentrality)
0124             {
0125                 hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive->Fill(PV_z, mbd_charge_asymm);
0126             }
0127         }
0128
0129         if (is_min_bias && firedTrig10 && goodMbdAsymm)
0130         {
0131             hM_INTTVtxZ_MBDAsymLe1->Fill(PV_z);
0132
0133             if (InttZvtxRange_narrow)
0134             {
0135                 hM_INTTVtxZ_MBDAsymLe1_VtxZm10to10->Fill(PV_z);
0136                 if (goodCentrality)
0137                 {
0138                     hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10->Fill(PV_z);
0139                     hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse->Fill(PV_z); // for reweighting
0140                     hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse->Fill(mbd_z_vtx); // for reweighting
0141                     hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10->Fill(PV_z, mbd_charge_asymm);
0142                 }
0143             }
0144
0145             if (InttZvtxRange_wide)
0146             {
0147                 if (goodCentrality)
0148                 {
0149                     hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse->Fill(PV_z); // for reweighting
0150                     hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse->Fill(mbd_z_vtx);
0151                 }
0152             }
0153
0154             if (goodCentrality)
0155             {
0156                 hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1->Fill(PV_z, mbd_z_vtx);
0157             }
0158         }
0159
0160         for (int i = 0; i < centrality_cut.size() - 1; i++)
0161         {
0162             if (is_min_bias && firedTrig10 && MBD_centrality >= centrality_cut[i] && MBD_centrality < centrality_cut[i + 1])
0163             {
0164                 hM_INTTVtxZ_Centrality[i]->Fill(PV_z);
0165                 hM_MBDVtxZ_Centrality[i]->Fill(mbd_z_vtx);
0166                 hM_INTTVtxZ_MBDVtxZ_Centrality[i]->Fill(PV_z, mbd_z_vtx);
0167
0168                 if (goodMbdAsymm)
0169                 {
0170                     hM_INTTVtxZ_Centrality_MBDAsymLe1[i]->Fill(PV_z);
0171                     if (InttZvtxRange_narrow)
0172                     {
0173                         hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10[i]->Fill(PV_z);
0174                         hM_INTTVtxZ_MBDAsymm_Centrality_MBDAsymLe1_VtxZm10to10[i]->Fill(PV_z, mbd_charge_asymm);
0175                     }
0176                 }
0177             }
0178         }
0179     }
0180
0181     TFile *fout = new TFile(outfilename, "RECREATE");
0182     fout->cd();
0183     hM_INTTVtxZ_MBDVtxZ_Inclusive->Write();
0184     hM_INTTVtxZ_Inclusive->Write();
0185     hM_INTTVtxZ_MBDAsymLe1->Write();
0186     hM_INTTVtxZ_MBDAsymLe1_VtxZm10to10->Write();
0187     hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10->Write();
0188     hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse->Write(); // for reweighting
0189     hM_INTTVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse->Write(); // for reweighting
0190     hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm30to30_coarse->Write(); // for reweighting
0191     hM_MBDVtxZ_Centrality0to70_MBDAsymLe1_VtxZm10to10_coarse->Write(); // for reweighting
0192     hM_INTTVtxZ_MBDAsymm_Centrality0to70_MBDAsymLe1_VtxZm10to10->Write();
0193     hM_INTTVtxZ_MBDAsymm_Centrality0to70_Inclusive->Write();
0194     hM_INTTVtxZ_MBDVtxZ_Centrality0to70_MBDAsymLe1->Write();
0195     for (int i = 0; i < centrality_cut.size() - 1; i++)
0196     {
0197         hM_INTTVtxZ_Centrality[i]->Write();
0198         hM_MBDVtxZ_Centrality[i]->Write();
0199         hM_INTTVtxZ_MBDVtxZ_Centrality[i]->Write();
0200         hM_INTTVtxZ_Centrality_MBDAsymLe1[i]->Write();
0201         hM_INTTVtxZ_Centrality_MBDAsymLe1_VtxZm10to10[i]->Write();
0202         hM_INTTVtxZ_MBDAsymm_Centrality_MBDAsymLe1_VtxZm10to10[i]->Write();
0203     }
0204     fout->Close();
0205 }