EMCal/EMCALPositionScanCodes/FillEnergy.C

0001 #include <cmath>
0002 #include <TFile.h>
0003 #include <TString.h>
0004 #include <TLine.h>
0005 #include <TTree.h>
0006 #include <TLatex.h>
0007 #include <TGraphErrors.h>
0008 #include <cassert>
0009 #include "SaveCanvas.C"
0010 #include "SetOKStyle.C"
0011 #include <iostream>
0012 #include <fstream>
0013 #include "TROOT.h"
0014 #include "TH1.h"
0015 #include "TTree.h"
0016 using namespace std;
0017
0018
0019
0020 using std::cout;
0021 using std::endl;
0022 #endif
0023
0024
0025 void FillEnergy()
0026 {
0027   gSystem->Load("libg4eval.so");
0028   gSystem->Load("libqa_modules.so");
0029   gSystem->Load("libPrototype3.so");
0030  gStyle->SetOptFit(0);
0031    gStyle->SetOptStat(0);
0032
0033 double inteval = 1;
0034
0035 double Xmin = 170;
0036 double Xmax = 310.0;
0037
0038 double Ymin = 60;
0039 double Ymax = 245.0;
0040
0041 int XBins =(Xmax - Xmin)/inteval;
0042 int YBins = (Ymax - Ymin)/inteval;
0043
0044
0045
0046 int Ini = 3543;
0047 int Final = 3579;
0048
0049  int step = 1;
0050
0051 int num = (Final - Ini)/step;
0052
0053 cout << "num = " << num << endl;
0054
0055 const int N = 8;
0056
0057 const int NTotal = N * N;
0058
0059 double x;
0060     double y;
0061 char Filename[512];
0062 /*
0063
0064 double Horizontal[N];
0065 double Vertical[N];
0066 double Energy[N];
0067
0068 */
0069 char inputfile[512];
0070
0071 int j;
0072
0073 double xmax = 0;
0074 double ymax = 0;
0075
0076 double Energy;
0077 double Energy_err;
0078 double Width;
0079 double Width_err;
0080 double Center;
0081 double Center_err;
0082
0083 int index;
0084
0085
0086 //
0087
0088 for(int i =0; i < num; i ++)
0089 {
0090 index = 0;
0091
0092  gStyle->SetOptFit(0);
0093    gStyle->SetOptStat(0);
0094
0095 TCanvas *c1 = new TCanvas("c1", "c1",0,0,800,600);
0096 c1->cd();
0097
0098
0099 cout << " i = " << i << endl;
0100
0101
0102
0103 j = i * step + Ini;
0104
0105 cout << " j = " << j << endl;
0106
0107 if(j != 3416 && j != 3418 && j != 3420 && j != 3421 && j != 3426)
0108 {
0109
0110 //sprintf(inputfile,"/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/Production_0128_WithHCalCalib/beam_0000%d-0000_DSTReader.root",j);
0111
0112 sprintf(inputfile,"/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2017/Production_0130_WithEMCalCalib/beam_0000%d-0000_DSTReader.root",j);
0113
0114
0115 cout << "Infile name = " << inputfile << endl;
0116
0117 ifstream ifile(inputfile);
0118
0119
0120
0121 if(ifile){
0122
0123 TFile *fin = new TFile(inputfile);
0124
0125
0126
0127 TTree *t = (TTree *)fin->Get("T");
0128
0129   t->SetAlias("C2_Inner_e", "1*abs(TOWER_RAW_C2[2].energy)");
0130   t->SetAlias("C2_Outer_e", "1*abs(TOWER_RAW_C2[3].energy)");
0131 t->SetAlias("Average_column", "Sum$(TOWER_CALIB_CEMC.get_column() * TOWER_CALIB_CEMC.get_energy())/Sum$(TOWER_CALIB_CEMC.get_energy())");
0132   t->SetAlias("Average_HODO_HORIZONTAL", "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))");
0133   t->SetAlias("Valid_HODO_HORIZONTAL", "Sum$(abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) > 0");
0134   t->SetAlias("No_Triger_VETO", "Sum$(abs(TOWER_RAW_TRIGGER_VETO.energy)>15)==0");
0135  t->SetAlias("Valid_HODO_VERTICAL", "Sum$(abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) > 0");
0136   t->SetAlias("C2_Sum_e", "C2_Inner_e + C2_Outer_e");
0137   t->SetAlias("Average_HODO_VERTICAL","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))");
0138 //return;
0139   t->SetAlias("Energy_Sum_CEMC", "1*Sum$(TOWER_CALIB_CEMC.get_energy())");
0140
0141 //Cut on Events//
0142
0143 // TCut event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO && C2_Sum_e > 200 && fmod(Entry$,10)==0";
0144 TCut event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO && C2_Sum_e > 200";
0145
0146   t->Draw(">>EventList", event_sel);
0147
0148   TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
0149
0150 //cout << "Cut is " << event_sel.data() << endl;
0151
0152   cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected" << endl;
0153
0154   t->SetEventList(elist);
0155
0156
0157 TH1D *h5= new TH1D("h5","",200, -1, 10);
0158
0159
0160 /*
0161 t->Draw("beam_2CH_mm >> h1");
0162 t->Draw("beam_2CV_mm >> h2");
0163
0164
0165 t->Draw("Average_column >> h5");
0166
0167 t->Draw("Average_HODO_HORIZONTAL>>h3");
0168
0169 t->Draw("Average_HODO_VERTICAL>>h4");
0170
0171
0172 */
0173 sprintf(Filename,"Hisfiles130/His%d.root",i);
0174
0175 TFile *fout = new TFile(Filename,"RECREATE");
0176
0177   t->SetAlias("XPos", "beam_2CH_mm - 5* int( Average_HODO_HORIZONTAL + 0.5)");
0178   t->SetAlias("YPos", "beam_2CV_mm + 5* int(Average_HODO_VERTICAL + 0.5)");
0179
0180
0181
0182
0183 //cout <<"Start the scan..."<<endl;
0184 //t->Scan("XPos:YPos:Energy_Sum_CEMC");
0185 TH3F *Energyhis = new TH3F("Energyhis","",200,0,20,YBins,Ymin,Ymax,XBins,Xmin,Xmax);
0186
0187 // t->Draw("XPos:YPos:Energy_Sum_CEMC>>Energyhis(200,0,20,200,69.7,191.0,200,212,302.0)");
0188
0189  t->Draw("XPos:YPos:Energy_Sum_CEMC>>Energyhis");
0190
0191
0192 Energyhis->Write();
0193 //fout->Print();
0194
0195 /*
0196
0197 new TCanvas();
0198 t->Draw("XPos>>hXPos");
0199
0200 new TCanvas();
0201 t->Draw("YPos>>hYPos");
0202
0203 new TCanvas();
0204 t->Draw("Energy_Sum_CEMC>>hEnergy_CEMC");
0205
0206  fout->Write();
0207     fout->Print();
0208 */
0209 //t->Project("Energyhis","XPos:YPos:Energy_Sum_CEMC");
0210
0211 //TF1 *f1 = new TF1("f1","[0]*TMath::Exp(-(x - [1] )*(x - [1] )/(2 * [2] * [2] ))",7,12);
0212
0213
0214 //t->Draw("beam_2CV_mm >> h2");
0215
0216 //t->Draw("TOWER_CALIB_CEMC.energy >> h3");
0217
0218 //return;
0219 //
0220 //
0221 /*
0222 cout <<"n_eve = "<<n_eve<<", hist = "<<Energyhis->GetSum()<<endl;
0223 gDirectory->Print();
0224
0225 TCanvas *c22 = new TCanvas("c22", "c22",0,0,800,600);
0226
0227 c22->cd();
0228
0229 Energyhis->Draw();
0230
0231     c22->Update();
0232
0233   c22->SaveAs("Energy.png");
0234
0235
0236 */
0237
0238
0239 /*
0240   Energyhis->Fit("gaus");
0241
0242 c1->Update();
0243
0244   c1->SaveAs("AVHis.png");
0245
0246
0247 Energy = gaus->GetParameter(0);
0248 Energy_err = gaus->GetParError(0);
0249 Center = gaus->GetParameter(1);
0250 Center_err = gaus->GetParError(1);
0251 Width = gaus->GetParameter(2);
0252 Width_err = gaus->GetParError(2);
0253
0254 */
0255 //cout << "Energy =  "  << Energy << endl;
0256 /*
0257 for(int j = 0; j < N; j++)
0258 {
0259
0260 for(int k = 0; k < N; k++)
0261 {
0262
0263 x = h1->GetMean()-0.5*j;
0264 y = h2->GetMean()-0.5*k;
0265
0266
0267 //E[i] = h3->GetMean();
0268
0269 cout << "x = " <<  x << endl;
0270
0271 cout << "y = " <<  y << endl;
0272
0273
0274 Energyhis->Fill(x,y,Energy);
0275
0276 index = index + 1;
0277
0278 cout << "index = " << index << endl;
0279
0280 }
0281
0282 }
0283 */
0284
0285
0286 }
0287
0288 if (!ifile) {
0289
0290     i = i + 1;
0291
0292 cout << "File  " << i << "does not exist" << endl;
0293
0294 }
0295
0296
0297 }
0298
0299 }
0300
0301 /*
0302 TCanvas *c11 = new TCanvas("c11", "c11",0,0,800,600);
0303 c11->cd();
0304    gStyle->SetEndErrorSize(0.01);
0305    c11->SetFillColor(10);
0306    c11->SetBorderMode(0);
0307    c11->SetBorderSize(2);
0308    c11->SetFrameFillColor(0);
0309    c11->SetFrameBorderMode(0);
0310
0311    c11->SetGridx();
0312    c11->SetGridy();
0313
0314    c11->SetLeftMargin(0.13);
0315    c11->SetBottomMargin(0.13);
0316    c11->SetTopMargin(0.02);
0317    c11->SetRightMargin(0.06);
0318
0319    double x1 = 0;
0320    double x2 = xmax*1.2;
0321    double y1 = 0;
0322    double y2 = ymax*1.2;
0323
0324    TH1D *d0 = new TH1D("d0","",100,x1,x2);
0325    d0->SetMinimum(y1);
0326    d0->SetMaximum(y2);
0327    d0->GetXaxis()->SetNdivisions(208);
0328    d0->GetXaxis()->SetTitle("x");
0329    d0->GetXaxis()->SetTitleOffset(0.9);
0330    d0->GetXaxis()->SetTitleSize(0.06);
0331    d0->GetXaxis()->SetLabelOffset(0.01);
0332    d0->GetXaxis()->SetLabelSize(0.045);
0333    d0->GetXaxis()->SetLabelFont(42);
0334    d0->GetXaxis()->SetTitleFont(42);
0335    d0->GetYaxis()->SetNdivisions(205);
0336    d0->GetYaxis()->SetTitle("Energy");
0337    d0->GetYaxis()->SetTitleOffset(1.0);
0338    d0->GetYaxis()->SetTitleSize(0.06);
0339    d0->GetYaxis()->SetLabelOffset(0.005);
0340    d0->GetYaxis()->SetLabelSize(0.045);
0341    d0->GetYaxis()->SetLabelFont(42);
0342    d0->GetYaxis()->SetTitleFont(42);
0343    d0->SetLineWidth(2);
0344    d0->Draw();
0345
0346    TLine *l1 = new TLine(x1,y1,x2,y1);
0347    l1->SetLineWidth(3);
0348    l1->Draw("same");
0349    TLine *l2 = new TLine(x1,y2,x2,y2);
0350    l2->SetLineWidth(3);
0351    l2->Draw("same");
0352    TLine *l3 = new TLine(x1,y1,x1,y2);
0353    l3->SetLineWidth(3);
0354    l3->Draw("same");
0355    TLine *l4 = new TLine(x2,y1,x2,y2);
0356    l4->SetLineWidth(3);
0357    l4->Draw("same");
0358
0359 TGraph *gr = new TGraphErrors(N, x, y);
0360
0361
0362
0363
0364     gr->SetMarkerSize(1.1);
0365    gr->SetMarkerStyle(20);
0366    gr->SetLineWidth(2);
0367     gr->SetName("Energy vs x");
0368
0369     gr->Draw("p");
0370
0371     c11->Update();
0372
0373   c11->SaveAs("AVCH.png");
0374 */
0375
0376
0377
0378 }