EMCal/EMCALmacros/DataAna.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 <iostream>
0010 #include <fstream>
0011 #include "TROOT.h"
0012 #include "TH1.h"
0013 #include "TTree.h"
0014 using namespace std;
0015
0016
0017
0018 using std::cout;
0019 using std::endl;
0020 #endif
0021
0022
0023
0024 void DataAna(){
0025  gStyle->SetOptFit(0);
0026    gStyle->SetOptStat(0);
0027
0028 TFile *fin = new TFile("Interpolated2.root");
0029
0030 const int NBins = 10;
0031 int bins = 5;
0032 double step = 1;
0033
0034 const int All = 7;
0035 double inteval = 12.5;
0036 double xinitial = -188;
0037 //double x[All];
0038 double yinitial = -77;
0039 double x;
0040 char HistoName[512];
0041 //Set Tower Vertical Lines//
0042
0043 double x1 = 205;
0044 double x2 = 228;
0045 double x3 = 253;
0046 double x4 = 277;
0047 double x5 = 301;
0048 char outname[512];
0049 /*
0050 for(int i = 0; i < All; i ++)
0051 {
0052
0053
0054
0055
0056
0057 }
0058
0059 */
0060
0061 double Xmin = -340;
0062 double Xmax = -170;
0063
0064 double Ymax = -60;
0065 double Ymin = -245.0;
0066
0067 int binlow;
0068
0069 int binhigh;
0070
0071 double Center;
0072
0073 double xlow;
0074 double xhigh;
0075
0076 double ylow;
0077 double yhigh;
0078 int ybinlow;
0079 int ybinhigh;
0080 double y;
0081
0082 int ycent;
0083
0084 double mean;
0085
0086 int YBins;
0087
0088 int N;
0089
0090 int Ybinmax;
0091
0092 int Ybinmin;
0093
0094 char Title[512];
0095
0096 for(int j = 0; j < All; j++)
0097 {
0098 TCanvas *c22 = new TCanvas("c22", "c22",0,0,800,600);
0099
0100 c22->cd();
0101
0102     TFile *fin = new TFile("Interpolated2.root");
0103     x = xinitial - inteval * j;
0104
0105
0106 //cout << "x = " << x << endl;
0107 //cout << "j  =  "  << j << endl;
0108
0109 Center = Inter->GetXaxis()->FindBin(x);
0110
0111
0112
0113 //double x[NBins];
0114 //int xbin[NBins];
0115
0116 cout << "OK 0" << endl;
0117
0118  xlow = x - bins * step;
0119  xhigh = x + bins * step;
0120
0121  binlow = Inter->GetXaxis()->FindBin(xlow);
0122
0123  binhigh = Inter->GetXaxis()->FindBin(xhigh);
0124
0125  YBins = (Ymax - Ymin)/step;
0126
0127 cout << "binlow =  " << binlow << endl;
0128 cout << "binhigh =  " << binhigh << endl;
0129
0130
0131 TH1D *h1 = new TH1D("h1","",YBins,Ymin,Ymax);
0132
0133
0134 Inter->ProjectionY("h1",binlow,binhigh);
0135
0136 //cout << "mean = " << Inter->GetMean() << endl;
0137
0138
0139
0140 //h1->Draw();
0141
0142
0143  Ybinmin = h1->GetXaxis()->FindBin(Ymin);
0144
0145
0146  Ybinmax = h1->GetXaxis()->FindBin(Ymax);
0147 cout << " OK 1" << endl;
0148
0149  N = (Ybinmax - Ybinmin)/NBins;
0150 /*
0151 for(int i = 0; i < N; i++)
0152 {
0153
0154     y = Ymin + (i+0.5) * NBins;
0155
0156 cout << "y = " <<  y << endl;
0157 ylow = Ymin + i * NBins;
0158
0159 yhigh = Ymin + (i+1)*NBins;
0160
0161 ybinlow = h1->GetXaxis()->FindBin(ylow);
0162 ybinhigh = h1->GetXaxis()->FindBin(yhigh);
0163
0164 ycent = h2->GetXaxis()->FindBin(y);
0165
0166
0167 mean = h1->Integral(ybinlow,ybinhigh)/NBins/(ybinhigh - ybinlow);
0168
0169 cout << "Difference  = " << ybinhigh - ybinlow << endl;
0170 cout << "y Bin = " << ycent << endl;
0171 cout << "Average is = " << mean << endl;
0172
0173 h2->SetBinContent(ycent,mean);
0174
0175 }
0176 */
0177
0178
0179 h1->Scale(0.1);
0180 //h1->Draw();
0181 //h2 = h1;
0182 //h2->Draw();
0183
0184 sprintf(HistoName,"Result%d.root",j);
0185
0186   TFile *fout = new TFile(HistoName,"RECREATE");
0187   h1->GetXaxis()->SetTitle("y (mm)");
0188   h1->GetYaxis()->SetTitle("Average Energy (GeV)");
0189 int k = 7- (j+1)/2;
0190 int l = 6- (j+1)/2;
0191 int m = 7 - (j+2)/2;
0192
0193
0194 if(j == 1 || j == 3 || j == 5) sprintf(Title,"10 X Bin Average Energy vs Vertical Position Tower Between %d and %d",k,l);
0195
0196 if(j == 0 || j == 2 || j == 4 || j == 6 ) sprintf(Title,"10 X Bin Average Energy vs Vertical Position Tower %d",m);
0197
0198   h1->SetTitle(Title);
0199
0200   h1->Write();
0201
0202
0203
0204 //return;
0205
0206 h1->Draw("same");
0207
0208 sprintf(outname,"Result/Energy vs 10-Bin Average Y %d.png",j);
0209
0210
0211
0212 c22->Update();
0213
0214 c22->SaveAs(outname);
0215
0216 //cout << " OK 2" << endl;
0217 }
0218
0219 int HAll = All + 4;
0220
0221 for(int j = 0; j < HAll; j++)
0222 {
0223 TCanvas *c22 = new TCanvas("c22", "c22",0,0,800,600);
0224
0225 c22->cd();
0226
0227     TFile *fin = new TFile("Interpolated2.root");
0228     y = yinitial - inteval * j;
0229
0230
0231 //cout << "x = " << x << endl;
0232 //cout << "j  =  "  << j << endl;
0233
0234 Center = Inter->GetXaxis()->FindBin(x);
0235
0236
0237
0238 //double x[NBins];
0239 //int xbin[NBins];
0240
0241 //cout << "OK 0" << endl;
0242
0243  ylow = y - bins * step;
0244  yhigh = y + bins * step;
0245
0246 ybinlow = Inter->GetYaxis()->FindBin(ylow);
0247
0248  ybinhigh = Inter->GetYaxis()->FindBin(yhigh);
0249
0250
0251
0252  YBins = (Ymax - Ymin)/step;
0253
0254 cout << "ybinlow =  " <<ybinlow << endl;
0255 cout << "ybinhigh =  " << ybinhigh << endl;
0256
0257
0258 TH1D *h1 = new TH1D("h1","",YBins,Ymin,Ymax);
0259
0260
0261 Inter->ProjectionX("h1",ybinlow,ybinhigh);
0262
0263 //cout << "mean = " << Inter->GetMean() << endl;
0264
0265
0266
0267 //h1->Draw();
0268
0269
0270  Ybinmin = h1->GetXaxis()->FindBin(Ymin);
0271
0272
0273  Ybinmax = h1->GetXaxis()->FindBin(Ymax);
0274 cout << " OK 1" << endl;
0275
0276  N = (Ybinmax - Ybinmin)/NBins;
0277 /*
0278 for(int i = 0; i < N; i++)
0279 {
0280
0281     y = Ymin + (i+0.5) * NBins;
0282
0283 cout << "y = " <<  y << endl;
0284 ylow = Ymin + i * NBins;
0285
0286 yhigh = Ymin + (i+1)*NBins;
0287
0288 ybinlow = h1->GetXaxis()->FindBin(ylow);
0289 ybinhigh = h1->GetXaxis()->FindBin(yhigh);
0290
0291 ycent = h2->GetXaxis()->FindBin(y);
0292
0293
0294 mean = h1->Integral(ybinlow,ybinhigh)/NBins/(ybinhigh - ybinlow);
0295
0296 cout << "Difference  = " << ybinhigh - ybinlow << endl;
0297 cout << "y Bin = " << ycent << endl;
0298 cout << "Average is = " << mean << endl;
0299
0300 h2->SetBinContent(ycent,mean);
0301
0302 }
0303 */
0304
0305
0306 h1->Scale(0.1);
0307 //h1->Draw();
0308 //h2 = h1;
0309 //h2->Draw();
0310
0311 sprintf(HistoName,"Result%d.root",j);
0312
0313   TFile *fout = new TFile(HistoName,"RECREATE");
0314   h1->GetXaxis()->SetTitle("x (mm)");
0315   h1->GetYaxis()->SetTitle("Average Energy (GeV)");
0316 int k = 7- (j+1)/2;
0317 int l = 6- (j+1)/2;
0318 int m = 7 - (j+2)/2;
0319
0320
0321 if(j == 1 || j == 3 || j == 5 || j == 7 || j == 9) sprintf(Title,"10 X Bin Average Energy vs Horizontal Position Tower Between %d and %d",k,l);
0322
0323 if(j == 0 || j == 2 || j == 4 || j == 6 || j == 8 || j == 10 ) sprintf(Title,"10 X Bin Average Energy vs Horizontal Position Tower %d",m);
0324
0325   h1->SetTitle(Title);
0326
0327   h1->Write();
0328
0329
0330
0331 //return;
0332
0333 h1->Draw("same");
0334
0335 sprintf(outname,"Result/Energy vs 10-Bin Average X %d.png",j);
0336
0337
0338
0339 c22->Update();
0340
0341 c22->SaveAs(outname);
0342
0343 cout << " OK 2" << endl;
0344 }
0345
0346
0347
0348 }