tutorial/macros/ExampleAnalysisDSTReader.C

0001 // $Id: $
0002
0003 /*!
0004  * \file ExampleAnalysisDSTReader.C
0005  * \brief example plotting macro using DSTReader TTrees from the production output
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include <cmath>
0012 #include <TFile.h>
0013 #include <TString.h>
0014 #include <TLine.h>
0015 #include <TTree.h>
0016 #include <TLatex.h>
0017 #include <TGraphErrors.h>
0018 #include <cassert>
0019 using namespace std;
0020
0021 void
0022 ExampleAnalysisDSTReader(
0023
0024 )
0025 {
0026   // here, load the lib for prototype3 so we can use get_*() functions in T->Draw functions.
0027   gSystem->Load("libPrototype3.so");
0028
0029   // Let's take a look at run 3533, which is -8 GeV/c secondary beam centered on EMCal tower 45.
0030   // More runs are produced in data production: https://wiki.bnl.gov/sPHENIX/index.php/2017_calorimeter_beam_test/Data_Production_and_Analysis#Production_Information
0031   // This tutorial focuses on the *_DSTReader.root, which is interctive TTree file containing all tower and channels.
0032   TFile * file =
0033       TFile::Open(
0034           "/gpfs/mnt/gpfs02/sphenix/data/data01/t1044-2016a/production.2017/Production_0203_init/beam_00003533-0000_DSTReader.root");
0035
0036   // Now you have a TTree called T accessible in your command line. Nevertheless, let's get it to use programtically.
0037   assert(file);
0038   TTree* T = (TTree *) file->GetObjectChecked("T", "TTree");
0039   assert(T);
0040
0041   // To make life easier, TTree allow us to define Alias of formulas, which can be used as if TTree varialbles.
0042
0043   // Cherenkov C2 energy sum
0044   T->SetAlias("C2_Sum_e",
0045       "TOWER_CALIB_C2[2].energy + TOWER_CALIB_C2[3].energy");
0046
0047   // Average column and rows on EMCal
0048   T->SetAlias("Average_column",
0049       "Sum$(TOWER_CALIB_CEMC.get_column() * TOWER_CALIB_CEMC.get_energy())/Sum$(TOWER_CALIB_CEMC.get_energy())");
0050   T->SetAlias("Average_row",
0051       "Sum$(TOWER_CALIB_CEMC.get_row() * TOWER_CALIB_CEMC.get_energy())/Sum$(TOWER_CALIB_CEMC.get_energy())");
0052
0053   // Average hodoscope and whether there is a valid hit
0054   T->SetAlias("Average_HODO_VERTICAL",
0055       "Sum$(TOWER_CALIB_HODO_VERTICAL.towerid * (abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))/Sum$((abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))");
0056   T->SetAlias("Valid_HODO_VERTICAL",
0057       "Sum$(abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) > 0");
0058
0059   T->SetAlias("Average_HODO_HORIZONTAL",
0060       "Sum$(TOWER_CALIB_HODO_HORIZONTAL.towerid * (abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))/Sum$((abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))");
0061   T->SetAlias("Valid_HODO_HORIZONTAL",
0062       "Sum$(abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) > 0");
0063
0064   T->SetAlias("Hodoscope_h3_v2",
0065       "abs(Average_HODO_HORIZONTAL-3)<.5 && abs(Average_HODO_VERTICAL-2)<.5");
0066
0067   // Nothing in trigger veto counter
0068   T->SetAlias("No_Triger_VETO",
0069       "Sum$(abs(TOWER_RAW_TRIGGER_VETO.energy)>15)==0");
0070
0071   // EMCal, HCal_in and HCal_out energy sums
0072   T->SetAlias("Energy_Sum_CEMC", "1*Sum$(TOWER_CALIB_CEMC.get_energy())");
0073   T->SetAlias("Energy_Sum_HCALIN",
0074       "1*Sum$(TOWER_CALIB_LG_HCALIN.get_energy())");
0075   T->SetAlias("Energy_Sum_HCALOUT",
0076       "1*Sum$(TOWER_CALIB_LG_HCALOUT.get_energy())");
0077
0078   // Good event
0079   T->SetAlias("Valid_Hodoscope_Clean_VETO_Counter",
0080       "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO");
0081   T->SetAlias("CherenkovCut", "C2_Sum_e>200");
0082
0083   // Now make the plots
0084   TText * t;
0085   TCanvas *c1 = new TCanvas("ExampleAnalysisDSTReader",
0086       "ExampleAnalysisDSTReader", 1200, 1000);
0087   c1->Divide(2, 2);
0088   int idx = 1;
0089   TPad * p;
0090
0091   c1->Update();
0092   p = (TPad *) c1->cd(idx++);
0093   p->SetLogy();
0094
0095   T->Draw("C2_Sum_e>>hC2_Sum_e(1000,0,2000)",
0096       "Valid_Hodoscope_Clean_VETO_Counter", "");
0097   T->Draw("C2_Sum_e>>hC2_Sum_e_c(1000,0,2000)",
0098       "Valid_Hodoscope_Clean_VETO_Counter && CherenkovCut", "same");
0099
0100   hC2_Sum_e->SetTitle(
0101       "Cherenkov C2 Sum, inclusive (blue) and with C2 cut (green);Cherenkov C2 Sum (ADC)");
0102 //  hC2_Sum_e_c->SetLineColor(kGreen + 2);
0103   hC2_Sum_e_c->SetFillColor(kGreen);
0104
0105   c1->Update();
0106   p = (TPad *) c1->cd(idx++);
0107   p->SetLogy();
0108
0109   T->Draw("Energy_Sum_CEMC>>hEnergy_Sum_CEMC(400,0,20)",
0110       "Valid_Hodoscope_Clean_VETO_Counter", "");
0111   T->Draw("Energy_Sum_CEMC>>hEnergy_Sum_CEMC_C(400,0,20)",
0112       "Valid_Hodoscope_Clean_VETO_Counter && CherenkovCut", "same");
0113   hEnergy_Sum_CEMC->SetTitle(
0114       "EMCal Sum, inclusive (black) and with C2 cut (green);EMCal Calibrated Energy Sum (GeV)");
0115 //  hEnergy_Sum_CEMC_C->SetLineColor(kGreen + 2);
0116   hEnergy_Sum_CEMC_C->SetFillColor(kGreen);
0117
0118   c1->Update();
0119   p = (TPad *) c1->cd(idx++);
0120   p->SetLogy();
0121
0122   T->Draw(
0123       "Energy_Sum_CEMC + Energy_Sum_HCALIN + Energy_Sum_HCALOUT>>hEnergy_Sum(400,0,20)",
0124       "Valid_Hodoscope_Clean_VETO_Counter", "");
0125   hEnergy_Sum->SetTitle(
0126       "All three calorimeter energy Sum;Energy Sum (GeV)");
0127
0128   c1->Update();
0129   p = (TPad *) c1->cd(idx++);
0130   p->SetLogy();
0131
0132   TH1 * hEnergy_Sum_CEMC_v2h3 = new TH1F("hEnergy_Sum_CEMC_v2h3",
0133       ";EMCal Calibrated Energy Sum (GeV);Count / bin", 100, 0, 20);
0134
0135   T->Draw("Energy_Sum_CEMC>>hEnergy_Sum_CEMC_v2h3",
0136       "Valid_Hodoscope_Clean_VETO_Counter && CherenkovCut && Hodoscope_h3_v2",
0137       "");
0138
0139   TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus",
0140       hEnergy_Sum_CEMC_v2h3->GetMean() - 2 * hEnergy_Sum_CEMC_v2h3->GetRMS(),
0141       hEnergy_Sum_CEMC_v2h3->GetMean() + 2 * hEnergy_Sum_CEMC_v2h3->GetRMS());
0142   fgaus_g->SetParameters(1,
0143       hEnergy_Sum_CEMC_v2h3->GetMean() - 2 * hEnergy_Sum_CEMC_v2h3->GetRMS(),
0144       hEnergy_Sum_CEMC_v2h3->GetMean() + 2 * hEnergy_Sum_CEMC_v2h3->GetRMS());
0145   hEnergy_Sum_CEMC_v2h3->Fit(fgaus_g, "MR0N");
0146
0147   TF1 * fgaus = new TF1("fgaus_LG", "gaus",
0148       fgaus_g->GetParameter(1) - 3 * fgaus_g->GetParameter(2),
0149       fgaus_g->GetParameter(1) + 3 * fgaus_g->GetParameter(2));
0150   fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
0151       fgaus_g->GetParameter(2));
0152   hEnergy_Sum_CEMC_v2h3->Fit(fgaus, "MR");
0153
0154   hEnergy_Sum_CEMC_v2h3->SetTitle(
0155       Form("EMCal energy (1x1 hodoscope H=3 V=2), #DeltaE/<E> = %.1f%%",
0156           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0157
0158 }
0159