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("libPrototype4.so");
0028
0029   // Let's take a look at run 668, which is -24 GeV/c secondary beam centered on EMCal tower 18.
0030   // More runs are produced in data production: https://wiki.bnl.gov/sPHENIX/index.php/2018_calorimeter_beam_test/Data_Production_and_Analysis#Production_output
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           "/sphenix/data/data03/phnxreco/sphenix/t1044/production/production_0306/beam_00000668-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[0].energy + TOWER_CALIB_C2[1].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
0055   T->SetAlias("Average_HODO_VERTICAL",
0056               "Sum$(TOWER_CALIB_HODO_VERTICAL.towerid * (abs(TOWER_CALIB_HODO_VERTICAL.energy)>0.5) * abs(TOWER_CALIB_HODO_VERTICAL.energy))/Sum$((abs(TOWER_CALIB_HODO_VERTICAL.energy)>0.5) * abs(TOWER_CALIB_HODO_VERTICAL.energy))");
0057   T->SetAlias("Valid_HODO_VERTICAL", "Sum$((TOWER_CALIB_HODO_VERTICAL.energy)>0.5) > 0");
0058
0059   T->SetAlias("Average_HODO_HORIZONTAL",
0060               "Sum$(TOWER_CALIB_HODO_HORIZONTAL.towerid * (abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>0.5) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))/Sum$((abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>0.5) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))");
0061   T->SetAlias("Valid_HODO_HORIZONTAL", "Sum$((TOWER_CALIB_HODO_HORIZONTAL.energy)>0.5) > 0");
0062
0063   T->SetAlias("Hodoscope_h3_v4",
0064       "abs(Average_HODO_HORIZONTAL-3)<.5 && abs(Average_HODO_VERTICAL-4)<.5");
0065
0066   // Nothing in trigger veto counter
0067   T->SetAlias("No_Triger_VETO",
0068               "Sum$((TOWER_CALIB_TRIGGER_VETO.energy)>0.2)==0");
0069
0070   // EMCal, HCal_in and HCal_out energy sums
0071   T->SetAlias("Energy_Sum_CEMC", "1*Sum$(TOWER_CALIB_CEMC.get_energy())");
0072   T->SetAlias("Energy_Sum_HCALIN",
0073       "1*Sum$(TOWER_CALIB_LG_HCALIN.get_energy())");
0074   T->SetAlias("Energy_Sum_HCALOUT",
0075       "1*Sum$(TOWER_CALIB_LG_HCALOUT.get_energy())");
0076
0077   // Good event
0078   T->SetAlias("Valid_Hodoscope_Clean_VETO_Counter",
0079       "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO");
0080   T->SetAlias("CherenkovCut", "C2_Sum_e>1000");
0081
0082   // Now make the plots
0083   TText * t;
0084   TCanvas *c1 = new TCanvas("ExampleAnalysisDSTReader",
0085       "ExampleAnalysisDSTReader", 1200, 1000);
0086   c1->Divide(2, 2);
0087   int idx = 1;
0088   TPad * p;
0089
0090   c1->Update();
0091   p = (TPad *) c1->cd(idx++);
0092   p->SetLogy();
0093
0094   T->Draw("C2_Sum_e>>hC2_Sum_e(1000,-1000,17000)",
0095       "Valid_Hodoscope_Clean_VETO_Counter", "");
0096   T->Draw("C2_Sum_e>>hC2_Sum_e_c(1000,-1000,17000)",
0097       "Valid_Hodoscope_Clean_VETO_Counter && CherenkovCut", "same");
0098
0099   hC2_Sum_e->SetTitle(
0100       "Cherenkov C2 Sum, inclusive (blue) and with C2 cut (green);Cherenkov C2 Sum (ADC)");
0101 //  hC2_Sum_e_c->SetLineColor(kGreen + 2);
0102   hC2_Sum_e_c->SetFillColor(kGreen);
0103
0104   c1->Update();
0105   p = (TPad *) c1->cd(idx++);
0106   p->SetLogy();
0107
0108   T->Draw("Energy_Sum_CEMC>>hEnergy_Sum_CEMC(400,0,40)",
0109       "Valid_Hodoscope_Clean_VETO_Counter", "");
0110   T->Draw("Energy_Sum_CEMC>>hEnergy_Sum_CEMC_C(400,0,40)",
0111       "Valid_Hodoscope_Clean_VETO_Counter && CherenkovCut", "same");
0112   hEnergy_Sum_CEMC->SetTitle(
0113       "EMCal Sum, inclusive (black) and with C2 cut (green);EMCal Calibrated Energy Sum (GeV)");
0114 //  hEnergy_Sum_CEMC_C->SetLineColor(kGreen + 2);
0115   hEnergy_Sum_CEMC_C->SetFillColor(kGreen);
0116
0117   c1->Update();
0118   p = (TPad *) c1->cd(idx++);
0119   p->SetLogy();
0120
0121   T->Draw(
0122       "Energy_Sum_CEMC + Energy_Sum_HCALIN + Energy_Sum_HCALOUT>>hEnergy_Sum(400,0,40)",
0123       "Valid_Hodoscope_Clean_VETO_Counter", "");
0124   hEnergy_Sum->SetTitle(
0125       "All three calorimeter energy Sum;Energy Sum (GeV)");
0126
0127   c1->Update();
0128   p = (TPad *) c1->cd(idx++);
0129   p->SetLogy();
0130
0131   TH1 * hEnergy_Sum_CEMC_v2h3 = new TH1F("hEnergy_Sum_CEMC_v2h3",
0132       ";EMCal Calibrated Energy Sum (GeV);Count / bin", 100, 0, 40);
0133
0134   T->Draw("Energy_Sum_CEMC>>hEnergy_Sum_CEMC_v2h3",
0135       "Valid_Hodoscope_Clean_VETO_Counter && CherenkovCut && Hodoscope_h3_v4",
0136       "");
0137
0138   TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus",
0139       hEnergy_Sum_CEMC_v2h3->GetMean() - 2 * hEnergy_Sum_CEMC_v2h3->GetRMS(),
0140       hEnergy_Sum_CEMC_v2h3->GetMean() + 2 * hEnergy_Sum_CEMC_v2h3->GetRMS());
0141   fgaus_g->SetParameters(1,
0142       hEnergy_Sum_CEMC_v2h3->GetMean() - 2 * hEnergy_Sum_CEMC_v2h3->GetRMS(),
0143       hEnergy_Sum_CEMC_v2h3->GetMean() + 2 * hEnergy_Sum_CEMC_v2h3->GetRMS());
0144   hEnergy_Sum_CEMC_v2h3->Fit(fgaus_g, "MR0N");
0145
0146   TF1 * fgaus = new TF1("fgaus_LG", "gaus",
0147       fgaus_g->GetParameter(1) - 3 * fgaus_g->GetParameter(2),
0148       fgaus_g->GetParameter(1) + 3 * fgaus_g->GetParameter(2));
0149   fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
0150       fgaus_g->GetParameter(2));
0151   hEnergy_Sum_CEMC_v2h3->Fit(fgaus, "MR");
0152
0153   hEnergy_Sum_CEMC_v2h3->SetTitle(
0154       Form("EMCal energy (1x1 hodoscope H=3 V=2), #DeltaE/<E> = %.1f%%",
0155           100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
0156
0157 }
0158