EMCal/macros/getMeanRMS.C

0001 #include <cmath>
0002 #include <vector>
0003 #include <TFile.h>
0004 #include <TString.h>
0005 #include <TLine.h>
0006 #include <TGraph.h>
0007 #include <TLatex.h>
0008 #include <TAxis.h>
0009 #include <TGraphErrors.h>
0010 #include <TMultiGraph.h>
0011 #include <TLegend.h>
0012 #include <TCanvas.h>
0013 #include <TF1.h>
0014 #include <TH1.h>
0015 #include <TMath.h>
0016 #include <cassert>
0017 #include <iostream>
0018 #include <fstream>
0019
0020 #include "SaveCanvas.C"
0021
0022 Double_t
0023 langaufun(Double_t *x, Double_t *par)
0024 {
0025
0026   //Fit parameters:
0027   //par[0]=Width (scale) parameter of Landau density
0028   //par[1]=Most Probable (MP, location) parameter of Landau density
0029   //par[2]=Total area (integral -inf to inf, normalization constant)
0030   //par[3]=Width (sigma) of convoluted Gaussian function
0031   //
0032   //In the Landau distribution (represented by the CERNLIB approximation),
0033   //the maximum is located at x=-0.22278298 with the location parameter=0.
0034   //This shift is corrected within this function, so that the actual
0035   //maximum is identical to the MP parameter.
0036
0037   // Numeric constants
0038   Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
0039   Double_t mpshift = -0.22278298; // Landau maximum location
0040
0041   // Control constants
0042   Double_t np = 100.0; // number of convolution steps
0043   Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
0044
0045   // Variables
0046   Double_t xx;
0047   Double_t mpc;
0048   Double_t fland;
0049   Double_t sum = 0.0;
0050   Double_t xlow, xupp;
0051   Double_t step;
0052   Double_t i;
0053
0054   // MP shift correction
0055   mpc = par[1] - mpshift * par[0];
0056
0057   // Range of convolution integral
0058   xlow = x[0] - sc * par[3];
0059   xupp = x[0] + sc * par[3];
0060
0061   step = (xupp - xlow) / np;
0062
0063   // Convolution integral of Landau and Gaussian by sum
0064   for (i = 1.0; i <= np / 2; i++)
0065     {
0066       xx = xlow + (i - .5) * step;
0067
0068       fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0069       sum += fland * TMath::Gaus(x[0], xx, par[3]);
0070
0071       xx = xupp - (i - .5) * step;
0072       fland = TMath::Landau(xx, mpc, par[0]) / par[0];
0073       sum += fland * TMath::Gaus(x[0], xx, par[3]);
0074     }
0075
0076   return (par[2] * step * sum * invsq2pi / par[3]);
0077 }
0078
0079 void
0080 getMeanRMS(void)
0081 {
0082   ofstream Output;
0083   Output.open("meanRMS_v1.txt");
0084
0085   char name[20000];
0086
0087   int runList[100] =
0088     { 2181, 2182, 2184, 2185, 2186, 2187, 2188, 2189, 2190, 2192 };
0089   double vert[10] =
0090     { 22, 42, 62, 82, 102, 122, 142, 152, 162, 172 };
0091   double mu[8][100];
0092   double rms[8][100];
0093   double err[8][100];
0094   double vt[8][100];
0095   double up[8][100];
0096   double mip[80] =
0097     { 43.4997, 45.6349, 40.8449, 47.7902, 48.4663, 39.2068, 48.3755, 83.1102 };
0098
0099   for (int j = 0; j < 8; j++)
0100     {
0101       sprintf(name,
0102           "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/Production_1101_EMCalSet2_HCalPR12/beam_0000%d-0000_DSTReader.root_DrawPrototype2MIPAnalysis_EMCDistribution_RAW_Valid_HODO_MIP_Col2_PedestalOther_TimingGood.root",
0103           runList[j]);
0104       TFile *_file0 = TFile::Open(name);
0105
0106       TCanvas *c1 = new TCanvas(Form("getMeanRMS_run%d", runList[j]),
0107           Form("getMeanRMS_run%d", runList[j]), 1800, 950);
0108       c1->Divide(4, 2);
0109
0110       for (int i = 0; i < 8; i++)
0111         {
0112           double min = 0;
0113           double max = 150;
0114           //double max=400;
0115           //if(runList[j]==2418 || runList[j]==2417)max=400;
0116           //if(runList[j]==2551)max=200;
0117           //if(runList[j]==2560)max=400;
0118           /*TF1 *fitFunc = new TF1("fitFunc","[0]*TMath::Landau(x,[1],[2])+[3]*TMath::Gaus(x,[4],[5])",min,max);
0119            //TF1 *fitFunc = new TF1("fitFunc","[0]*TMath::Landau(x,[1],[2])",min,max);
0120            //TF1 *fitFunc = new TF1("fitFunc","[0]*TMath::Gaus(x,[1],[2])",min,max);
0121            fitFunc->SetParameter(0,100);
0122            fitFunc->SetParameter(1,10);
0123            fitFunc->SetParameter(2,10);
0124            fitFunc->SetParameter(3,10);
0125            fitFunc->SetParameter(4,10);
0126            fitFunc->SetParameter(5,10);*/
0127
0128           c1->cd(i + 1);
0129           c1->cd(i + 1)->SetLogy();
0130           sprintf(name, "hEnergy_ieta2_iphi%d", i);
0131           TH1F *h = (TH1F*) gDirectory->Get(name)->Clone();
0132           assert(h);
0133           h->Draw();
0134
0135           TF1 *ffit = new TF1(TString(h->GetName()) + "_langaufun", langaufun,
0136               0, 500, 4);
0137           ffit->SetParameters(6, 50, 3000, 15);
0138           ffit->SetParNames("Width", "MP", "Area", "GSigma");
0139           ffit->SetParLimits(0,0,20);
0140           ffit->SetParLimits(1,0,200);
0141           ffit->SetParLimits(2,0,20000);
0142 //          ffit->SetParLimits(3,10,20);
0143           ffit->SetParLimits(3,1,30);
0144
0145           h->Fit(ffit, "RBM");
0146
0147           c1->Update();
0148
0149 //          Output << col[j] << "\t" << i << "\t" << fitFunc->GetMaximumX(20, 500)
0150 //              << endl;
0151           Output << runList[j] << "\t" << vert[j] << "\t" << i << "\t"
0152               << ffit->GetParameter(1) << "\t" << ffit->GetParError(1) << "\t"
0153               << h->GetEntries() << endl;
0154 //          mu[i][j]=h->GetMean();
0155           mu[i][j] = ffit->GetParameter(1);
0156 //          h->Fit("landau");
0157 //          mu[i][j] = landau->GetParameter(1);
0158 //          rms[i][j] = h->GetRMS();
0159           err[i][j] = ffit->GetParError(1);
0160           vt[i][j] = vert[j]; //+i*.1;
0161 //          up[i][j] =
0162 //              (h->GetMean() + 1.281 * h->GetRMS() / sqrt(h->GetEntries()))
0163 //                  / mip[i];
0164 //          //cout<<up[i][j]<<endl;
0165 //          Output << runList[j] << "\t" << vert[j] << "\t" << i << "\t"
0166 //              << h->GetMean() << "\t" << h->GetRMS() << "\t" << h->GetEntries()
0167 //              << "\t" << up[i][j] << endl;
0168
0169           //if(vert[j]<160 || vert[j]>170)vt[i][j]=vert[j]+i*.4;
0170
0171         }
0172 //      sprintf(name, "meanRMS_%d_v2.root", runList[j]);
0173 //      c1->Print(name);
0174       SaveCanvas(c1, TString(name) + TString(c1->GetName()), kTRUE);
0175     }
0176
0177 //  double avg[7];
0178 //  for (int j = 0; j < 7; j++) //veritcal
0179 //    {
0180 //      double sum = 0.;
0181 //      for (int i = 0; i < 8; i++) //tower
0182 //        {
0183 //          sum = sum + mu[i][j];
0184 //        }
0185 //      avg[j] = sum / 8.;
0186 //    }
0187
0188 //  TCanvas c4("c4");
0189 //  for (int i = 0; i < 8; i++)
0190 //    {
0191 //      TGraphErrors *gr = new TGraphErrors(7, vert, avg, 0, 0);
0192 //      gr->SetMarkerStyle(20);
0193 //      //gr->SetMaximum(27);
0194 //      //gr->SetMinimum(23);
0195 //      gr->SetTitle("");
0196 //      gr->GetXaxis()->SetTitle("vertical position (mm)");
0197 //      gr->GetYaxis()->SetTitle("<MPV ADC>");
0198 //      gr->Draw("ap");
0199 //    }
0200
0201   TLegend * leg = new TLegend(0.91, 0.51, 0.99, 0.89, NULL, "brNDC");
0202
0203   TCanvas * c1 = new TCanvas(Form("getMeanRMS_Summary"),
0204       Form("getMeanRMS_Summary"), 1800, 950);
0205   gPad->DrawFrame(0, 0, 200, 200);
0206   for (int i = 0; i < 8; i++)
0207     {
0208       TGraphErrors *gr = new TGraphErrors(7, vt[i], mu[i], 0, err[i]);
0209       //TGraphErrors *gr = new TGraphErrors(11,vt[i],mu[i],0,err[i]);
0210       gr->SetMarkerStyle(20);
0211       gr->SetMarkerColor(i + 1);
0212       //gr->SetMaximum(32);
0213       //gr->SetMinimum(-14);
0214       gr->SetMaximum(60);
0215       gr->SetMinimum(10);
0216       gr->SetTitle("");
0217       gr->GetXaxis()->SetTitle("vertical position (mm)");
0218       gr->GetYaxis()->SetTitle("MPV ADC");
0219 //      if (i == 0)
0220 //        gr->Draw("ap");
0221 //      if (i > 0)
0222       gr->Draw("p");
0223       //gr->Fit("pol1");
0224       //gr->GetXaxis()->TAxis::SetNdivisions(505);
0225       //gr->GetXaxis()->SetLimits(160.6,166);
0226       //gr->GetXaxis()->SetLimits(130,190);
0227       sprintf(name, "row %d", i);
0228       leg->AddEntry(gr, name, "p");
0229
0230 //      gr->Draw("p");
0231     }
0232   leg->Draw();
0233   SaveCanvas(c1, TString(c1->GetName()), kTRUE);
0234
0235 //  leg2 = new TLegend(0.54, 0.35, 0.62, 0.73, NULL, "brNDC");
0236 //  TCanvas c3("c3");
0237 //  for (int i = 0; i < 8; i++)
0238 //    {
0239 //      TGraphErrors *gr2 = new TGraphErrors(11, vert, up[i], 0, 0);
0240 //      //TGraphErrors *gr2 = new TGraphErrors(11,vert,mu[i],0,err[i]);
0241 //      gr2->SetMarkerStyle(20);
0242 //      gr2->SetMarkerColor(i + 1);
0243 //      //gr2->SetMaximum(32);
0244 //      //gr2->SetMinimum(-14);
0245 //      gr2->SetMaximum(.4);
0246 //      gr2->SetMinimum(-.04);
0247 //      gr2->SetTitle("");
0248 //      gr2->GetXaxis()->SetTitle("vertical position (mm)");
0249 //      gr2->GetYaxis()->SetTitle("Cherenkov/MIP");
0250 //      if (i == 0)
0251 //        gr2->Draw("ap");
0252 //      if (i > 0)
0253 //        gr2->Draw("p");
0254 //      //gr2->GetXaxis()->TAxis::SetNdivisions(505);
0255 //      //gr2->GetXaxis()->SetLimits(160.6,166);
0256 //      gr2->GetXaxis()->SetLimits(130, 190);
0257 //      sprintf(name, "row %d", i);
0258 //      leg2->AddEntry(gr2, name, "p");
0259 //    }
0260 //  leg2->Draw();
0261
0262 }
0263