EMCal/EMCALmacros/DataAna2.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 DataAna2(){
0025  gStyle->SetOptFit(0);
0026    gStyle->SetOptStat(0);
0027
0028 TFile *fin = new TFile("Interpolated.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 = -300;
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 if(0){
0097
0098 for(int j = 0; j < All; j++)
0099 {
0100 TCanvas *c22 = new TCanvas("c22", "c22",0,0,800,600);
0101
0102 c22->cd();
0103
0104     TFile *fin = new TFile("Interpolated2.root");
0105     x = xinitial - inteval * j;
0106
0107
0108 //cout << "x = " << x << endl;
0109 //cout << "j  =  "  << j << endl;
0110
0111 Center = Inter->GetXaxis()->FindBin(x);
0112
0113
0114
0115 //double x[NBins];
0116 //int xbin[NBins];
0117
0118 cout << "OK 0" << endl;
0119
0120  xlow = x - bins * step;
0121  xhigh = x + bins * step;
0122
0123  binlow = Inter->GetXaxis()->FindBin(xlow);
0124
0125  binhigh = Inter->GetXaxis()->FindBin(xhigh);
0126
0127  YBins = (Ymax - Ymin)/step;
0128
0129 cout << "binlow =  " << binlow << endl;
0130 cout << "binhigh =  " << binhigh << endl;
0131
0132
0133 TH1D *h1 = new TH1D("h1","",YBins,Ymin,Ymax);
0134
0135
0136 Inter->ProjectionY("h1",binlow,binhigh);
0137
0138 //cout << "mean = " << Inter->GetMean() << endl;
0139
0140
0141
0142 //h1->Draw();
0143
0144
0145  Ybinmin = h1->GetXaxis()->FindBin(Ymin);
0146
0147
0148  Ybinmax = h1->GetXaxis()->FindBin(Ymax);
0149 cout << " OK 1" << endl;
0150
0151  N = (Ybinmax - Ybinmin)/NBins;
0152 /*
0153 for(int i = 0; i < N; i++)
0154 {
0155
0156     y = Ymin + (i+0.5) * NBins;
0157
0158 cout << "y = " <<  y << endl;
0159 ylow = Ymin + i * NBins;
0160
0161 yhigh = Ymin + (i+1)*NBins;
0162
0163 ybinlow = h1->GetXaxis()->FindBin(ylow);
0164 ybinhigh = h1->GetXaxis()->FindBin(yhigh);
0165
0166 ycent = h2->GetXaxis()->FindBin(y);
0167
0168
0169 mean = h1->Integral(ybinlow,ybinhigh)/NBins/(ybinhigh - ybinlow);
0170
0171 cout << "Difference  = " << ybinhigh - ybinlow << endl;
0172 cout << "y Bin = " << ycent << endl;
0173 cout << "Average is = " << mean << endl;
0174
0175 h2->SetBinContent(ycent,mean);
0176
0177 }
0178 */
0179
0180
0181 h1->Scale(0.1);
0182 //h1->Draw();
0183 //h2 = h1;
0184 //h2->Draw();
0185
0186 sprintf(HistoName,"Result%d.root",j);
0187
0188   TFile *fout = new TFile(HistoName,"RECREATE");
0189   h1->GetXaxis()->SetTitle("y (mm)");
0190   h1->GetYaxis()->SetTitle("Average Energy (GeV)");
0191 int k = 7- (j+1)/2;
0192 int l = 6- (j+1)/2;
0193 int m = 7 - (j+2)/2;
0194
0195
0196 if(j == 1 || j == 3 || j == 5) sprintf(Title,"10 X Bin Average Energy vs Vertical Position Tower Between %d and %d",k,l);
0197
0198 if(j == 0 || j == 2 || j == 4 || j == 6 ) sprintf(Title,"10 X Bin Average Energy vs Vertical Position Tower %d",m);
0199
0200   h1->SetTitle(Title);
0201
0202   h1->Write();
0203
0204
0205
0206 //return;
0207
0208 h1->Draw("same");
0209
0210 sprintf(outname,"Result/Energy vs 10-Bin Average Y %d.png",j);
0211
0212
0213
0214 c22->Update();
0215
0216 c22->SaveAs(outname);
0217
0218 //cout << " OK 2" << endl;
0219 }
0220
0221 }
0222 int HAll = All + 4;
0223
0224 for(int j = 0; j < HAll; j++)
0225 {
0226 TCanvas *c22 = new TCanvas("c22", "c22",0,0,800,600);
0227
0228 c22->cd();
0229
0230     TFile *fin = new TFile("Interpolated2.root");
0231     y = yinitial - inteval * j;
0232
0233
0234 //cout << "x = " << x << endl;
0235 //cout << "j  =  "  << j << endl;
0236
0237 Center = Inter->GetXaxis()->FindBin(x);
0238
0239
0240
0241 //double x[NBins];
0242 //int xbin[NBins];
0243
0244 //cout << "OK 0" << endl;
0245
0246  ylow = y - bins * step;
0247  yhigh = y + bins * step;
0248
0249 ybinlow = Inter->GetYaxis()->FindBin(ylow);
0250
0251  ybinhigh = Inter->GetYaxis()->FindBin(yhigh);
0252
0253
0254
0255  YBins = (Ymax - Ymin)/step;
0256
0257 cout << "ybinlow =  " <<ybinlow << endl;
0258 cout << "ybinhigh =  " << ybinhigh << endl;
0259
0260
0261 TH1D *h1 = new TH1D("h1","",YBins,Ymin,Ymax);
0262
0263
0264 Inter->ProjectionX("h1",ybinlow,ybinhigh);
0265
0266 //cout << "mean = " << Inter->GetMean() << endl;
0267
0268
0269
0270 //h1->Draw();
0271
0272
0273  Ybinmin = h1->GetXaxis()->FindBin(Ymin);
0274
0275
0276  Ybinmax = h1->GetXaxis()->FindBin(Ymax);
0277 cout << " OK 1" << endl;
0278
0279  N = (Ybinmax - Ybinmin)/NBins;
0280 /*
0281 for(int i = 0; i < N; i++)
0282 {
0283
0284     y = Ymin + (i+0.5) * NBins;
0285
0286 cout << "y = " <<  y << endl;
0287 ylow = Ymin + i * NBins;
0288
0289 yhigh = Ymin + (i+1)*NBins;
0290
0291 ybinlow = h1->GetXaxis()->FindBin(ylow);
0292 ybinhigh = h1->GetXaxis()->FindBin(yhigh);
0293
0294 ycent = h2->GetXaxis()->FindBin(y);
0295
0296
0297 mean = h1->Integral(ybinlow,ybinhigh)/NBins/(ybinhigh - ybinlow);
0298
0299 cout << "Difference  = " << ybinhigh - ybinlow << endl;
0300 cout << "y Bin = " << ycent << endl;
0301 cout << "Average is = " << mean << endl;
0302
0303 h2->SetBinContent(ycent,mean);
0304
0305 }
0306 */
0307
0308
0309 h1->Scale(0.1);
0310 //h1->Draw();
0311 //h2 = h1;
0312 //h2->Draw();
0313
0314 sprintf(HistoName,"Result%d.root",j);
0315
0316   TFile *fout = new TFile(HistoName,"RECREATE");
0317   h1->GetXaxis()->SetTitle("x (mm)");
0318   h1->GetYaxis()->SetTitle("Average Energy (GeV)");
0319 int k = 7- (j+1)/2;
0320 int l = 6- (j+1)/2;
0321 int m = 7 - (j+2)/2;
0322
0323
0324 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);
0325
0326 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);
0327
0328   h1->SetTitle(Title);
0329
0330   h1->Write();
0331
0332
0333
0334 //return;
0335
0336 h1->Draw("same");
0337
0338 sprintf(outname,"Result2/Energy vs 10-Bin Average X %d.png",j);
0339
0340
0341
0342 c22->Update();
0343
0344 c22->SaveAs(outname);
0345
0346 cout << " OK 2" << endl;
0347 }
0348
0349
0350
0351 }