StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
electron_drawfits.C
1 /*
2  * electron_drawfits.C
3  * Original author: Matt Walker
4  * Update Author: J. Kevin Adkins, University of Kentucky
5  * Updates: Removed a lot of code. Much of it was commented
6  * out by Alice Ohlson when she made previous changes to the
7  * algorithm. Also, removed a lot of code which was unnecessary,
8  * and updated naming conventions for clearer reading.
9  *
10  * Most recent update: Feb 25, 2015
11  */
12 
13 #include <iostream>
14 #include <fstream>
15 using namespace std;
16 
17 #include "TFile.h"
18 #include "TH1.h"
19 #include "TH2.h"
20 #include "TF1.h"
21 #include "TMath.h"
22 #include "TPostScript.h"
23 #include "TCanvas.h"
24 #include "TStyle.h"
25 #include "TLatex.h"
26 #include "TTree.h"
27 #include "TLine.h"
28 
29 void electron_drawfits(const char* infile="CombinedElectrons.root", const char* mipGainFilename="mip.gains", const char* absoluteGainFilename="electron.gains",
30  const char* psName="electron.ps", const char* ffile="electron.fits")
31 {
32  cout << "Input File: " << infile << endl;
33  cout << "Plots File: " << psName << endl;
34 
35  // Load library and set up decoder & geometry
36  gROOT->Macro("loadMuDst.C");
37  StEmcDecoder *mDecoder = new StEmcDecoder();
38  StEmcGeom *mEmcGeom = StEmcGeom::instance("bemc");
39 
40  //**********************************************//
41  //Initialization //
42  //**********************************************//
43  const Int_t nTowers = 4800;
44  const Int_t nRings = 40;
45  const Int_t nCrates = 30;
46  const Int_t nCrateSlices = nCrates*20;
47  const Double_t pi = TMath::Pi();
48  Double_t mipGains[nTowers];
49  Int_t mipTowerStatus[nTowers];
50  Double_t mipAdc[nTowers];
51  Double_t mipAdcErr[nTowers];
52  Double_t ringPars[5], crateslicePars[5];
53 
54  ifstream mipGainFile(mipGainFilename);
55  while(1){
56  Int_t softId,towerStat;
57  Double_t mipAdcVal , mipAdcErrVal, mipGain;
58  mipGainFile >> softId >> mipAdcVal >> mipAdcErrVal >> towerStat;
59  if(!mipGainFile.good())break;
60  Float_t towerEta, towerTheta;
61  mEmcGeom->getEta(softId,towerEta);
62  mEmcGeom->getTheta(softId,towerTheta);
63  mipGain = 0.264*(1+0.056*towerEta*towerEta)/(sin(towerTheta)*mipAdcVal);
64  mipAdc[softId-1] = mipAdcVal;
65  mipGains[softId-1] = mipGain;
66  mipTowerStatus[softId-1] = towerStat;
67  mipAdcErr[softId-1] = mipAdcErrVal;
68  }
69  mipGainFile.close();
70 
71  TPostScript *ps = new TPostScript(psName);
72 
73  TFile *rootfile = new TFile(infile,"read");
74 
75  TH1 *crateslice_histo[nCrates][nRings/2];
76  TH1 *ring_histo[nRings];
77 
78  TF1 *ringFit[nRings], *ringExpoFit[nRings], *ringGausFit[nRings];
79  TF1 *cratesliceFit[nCrates][nRings/2], *cratesliceExpoFit[nCrates][nRings/2], *cratesliceGausFit[nCrates][nRings/2];
80 
81  Int_t ringIndex, sliceEtaIndex;
82 
83  // compare fit parameters
84  TH2F *hChi2EtaPhi = new TH2F("hChi2EtaPhi","Chi2/NDF",40,-1,1,120,-pi,pi);
85  TH2F *hNewGains = new TH2F("hNewGains","New Gains",40,-1.,1.,120,-pi,pi);
86 
87  char name[100], title[100];
88  for(Int_t iRing = 0; iRing < nRings; ++iRing){
89  Int_t ringId = iRing + 1;
90  sprintf(name,"ringHisto_%i",iRing+1);
91  ring_histo[iRing] = (TH1F*)rootfile->Get(name);
92  ring_histo[iRing]->SetTitle("");
93  }
94  for(Int_t iCrate = 0; iCrate < nCrates; ++iCrate){
95  Int_t crateId = iCrate + 1;
96  for(Int_t iRing = 0; iRing < nRings/2; ++iRing){
97  Int_t ringId = iRing + 1;
98  sprintf(name,"cratesliceHisto_%i_%i",crateId,ringId);
99  crateslice_histo[iCrate][iRing] = (TH1F*)rootfile->Get(name);
100  crateslice_histo[iCrate][iRing]->SetTitle("");
101  }
102  }
103 
104  //global graphics functions
105  gStyle->SetOptFit(111);
106  gStyle->SetCanvasColor(10);
107  gStyle->SetCanvasBorderMode(0);
108  gStyle->SetPalette(1);
109  gStyle->SetStatColor(0);
110  gStyle->SetOptDate(0);
111 
112  //**********************************************//
113  //Fit Ring Histograms //
114  //**********************************************//
115 
116  TCanvas *c = new TCanvas("c","",100,100,600.,800.);
117  ofstream fitfile(ffile);
118  for(Int_t iRing = 0; iRing < nRings; iRing++){
119  Int_t ringId = iRing + 1;
120  if(iRing%2==0){
121  c->Update();
122  ps->NewPage();
123  c->Clear();
124  c->Divide(1,2);
125  }
126  c->cd((iRing%2)+1);
127 
128  sprintf(name,"ringExpoFit_%i",ringId);
129  ringExpoFit[iRing] = new TF1(name,"expo",0.25,1.7);
130  sprintf(name,"ringGausFit_%i",ringId);
131  ringGausFit[iRing] = new TF1(name,"gaus",0.65,1.15);
132  sprintf(name,"ringFit_%i",ringId);
133  ringFit[iRing] = new TF1(name,"expo(0)+gaus(2)",0.25,1.7);
134 
135  ring_histo[iRing]->Fit(ringExpoFit[iRing],"RQ0");
136  ring_histo[iRing]->Fit(ringGausFit[iRing],"RQ0");
137  ringExpoFit[iRing]->GetParameters(&ringPars[0]);
138  ringGausFit[iRing]->GetParameters(&ringPars[2]);
139  ringFit[iRing]->SetParameters(ringPars);
140  ringFit[iRing]->SetParNames("ExpConst","ExpSlope","GausConst","GausMean","GausSigma");
141  ring_histo[iRing]->Fit(ringFit[iRing],"RQ0");
142 
143  ring_histo[iRing]->DrawCopy();
144  ringFit[iRing]->DrawCopy("same");
145 
146  Double_t ringMean = ringFit[iRing]->GetParameter(3);
147  Double_t ringErr = ringFit[iRing]->GetParError(3);
148 
149  fitfile << "Ring " << ringId << " " << ringMean << " " << ringErr << endl;
150  cout << "Ring " << ringId << " Fit: " << ringMean << " " << ringErr/ringMean << " " << ring_histo[iRing]->GetEntries() << endl;
151  }
152 
153  //**********************************************//
154  //Fit Crate Slice Histograms //
155  //**********************************************//
156  for(Int_t iCrate = 0; iCrate < nCrates; iCrate++){
157  Int_t crateId = iCrate + 1;
158  for(Int_t iRing = 0; iRing < nRings/2; iRing++){
159  Int_t ringId = iRing + 1;
160  if(iRing%10 == 0){
161  c->Update();
162  ps->NewPage();
163  c->Clear();
164  c->Divide(2,5);
165  }
166  c->cd((iRing%10)+1);
167 
168  sprintf(name,"cratesliceExpoFit_%i_%i",crateId,ringId);
169  cratesliceExpoFit[iCrate][iRing] = new TF1(name,"expo",0.25,1.7);
170  sprintf(name,"cratesliceGausFit_%i_%i",crateId,ringId);
171  cratesliceGausFit[iCrate][iRing] = new TF1(name,"gaus",0.65,1.15);
172  sprintf(name,"cratesliceFit_%i_%i",crateId,ringId);
173  cratesliceFit[iCrate][iRing] = new TF1(name,"expo(0)+gaus(2)",0.25,1.7);
174 
175  crateslice_histo[iCrate][iRing]->Fit(cratesliceExpoFit[iCrate][iRing],"RQ0");
176  crateslice_histo[iCrate][iRing]->Fit(cratesliceGausFit[iCrate][iRing],"RQ0");
177  cratesliceExpoFit[iCrate][iRing]->GetParameters(&crateslicePars[0]);
178  cratesliceGausFit[iCrate][iRing]->GetParameters(&crateslicePars[2]);
179  cratesliceFit[iCrate][iRing]->SetParameters(crateslicePars);
180  cratesliceFit[iCrate][iRing]->SetParNames("ExpConst","ExpSlope","GausConst","GausMean","GausSigma");
181  crateslice_histo[iCrate][iRing]->Fit(cratesliceFit[iCrate][iRing],"RQ0");
182 
183  crateslice_histo[iCrate][iRing]->DrawCopy();
184  cratesliceFit[iCrate][iRing]->DrawCopy("same");
185 
186  Double_t cratesliceMean = cratesliceFit[iCrate][iRing]->GetParameter(3);
187  Double_t cratesliceErr = cratesliceFit[iCrate][iRing]->GetParError(3);
188 
189  fitfile << "CrateSlice " << iCrate*20+iRing+1 << " " << cratesliceMean << " " << cratesliceErr << endl;
190  cout << "CrateSlice " << iCrate*20+iRing+1 << " Fit: " << cratesliceMean << " +/- " << cratesliceErr << endl;
191  }
192  }
193 
194  ofstream newGainFile(absoluteGainFilename);
195  Double_t absoluteGain, absoluteGainErr;
196  Double_t adjustVal, adjustErr;
197  Double_t mipGain, mipGainErr;
198  Int_t softId, crateId, sequence;
199  Float_t mTowerEta, mTowerPhi;
200  // Loop through all towers and calculate absoluate gains
201  for(Int_t iTow = 0; iTow < nTowers; iTow++){
202  absoluteGain = absoluteGainErr = 0.;
203  adjustVal = adjustErr = 0.;
204  mipGain = 0.;
205  crateId = sequence = -1;
206  mTowerEta = mTowerPhi = -1.;
207  softId = iTow + 1;
208  mDecoder->GetCrateFromTowerId(softId,crateId,sequence);
209  mEmcGeom->getEta(softId,mTowerEta);
210  mEmcGeom->getPhi(softId,mTowerPhi);
211  if (fabs(mTowerEta) > 0.965) // Outer tower Eta is 0.967, bump this up to 0.975 for calculating sliceEtaIndex correctly
212  mTowerEta += 0.008*fabs(mTowerEta)/mTowerEta;
213  sliceEtaIndex = ((TMath::Nint(fabs(mTowerEta) * 1000.0) + 25)/50 - 1);
214  ringIndex = ((TMath::Nint(mTowerEta * 1000.0) + 25)/50 + 19);
215  if(mipTowerStatus[iTow] == 1){ // Don't calculate if bad MIP status
216  mipGain = mipGains[iTow];
217  if (sliceEtaIndex <= 12){
218  adjustVal = cratesliceFit[crateId-1][sliceEtaIndex]->GetParameter(3);
219  adjustErr = cratesliceFit[crateId-1][sliceEtaIndex]->GetParError(3);
220  hChi2EtaPhi->Fill(mTowerEta,mTowerPhi,cratesliceFit[crateId-1][sliceEtaIndex]->GetChisquare()/(Double_t)cratesliceFit[crateId-1][sliceEtaIndex]->GetNDF());
221  }
222  else{ //sliceEtaIndex >= 13 - Use ring values on outer eta for combined statstics
223  adjustVal = ringFit[ringIndex]->GetParameter(3);
224  adjustErr = ringFit[ringIndex]->GetParError(3);
225  hChi2EtaPhi->Fill(mTowerEta,mTowerPhi,ringFit[ringIndex]->GetChisquare()/(Double_t)ringFit[ringIndex]->GetNDF());
226  }
227  absoluteGain = mipGain/adjustVal;
228  mipGainErr = mipAdcErr[iTow]*mipGains[iTow]/mipAdc[iTow];
229  absoluteGainErr = sqrt(pow(mipGain*adjustErr/(adjustVal*adjustVal),2) + pow(mipGainErr/adjustVal,2));
230  }
231  newGainFile << softId << " " << absoluteGain << " " << absoluteGainErr << " " << mipTowerStatus[iTow] << endl;
232  cout << "Relative Adjustment Value: " << adjustVal << " +/- " << adjustErr << endl;
233  cout << softId << " " << absoluteGain << " " << absoluteGainErr << " " << mipTowerStatus[iTow] << endl;
234  if (mipTowerStatus[iTow] == 1)
235  hNewGains->Fill(mTowerEta,mTowerPhi,absoluteGain);
236  }
237 
238  newGainFile.close();
239 
240  c->Update();
241  ps->NewPage();
242  c->Clear();
243  c->Divide(1,2);
244  c->cd(1);
245  hNewGains->SetTitle("Absolute Gains");
246  hNewGains->GetXaxis()->SetTitle("#eta");
247  hNewGains->GetXaxis()->CenterTitle();
248  hNewGains->GetYaxis()->SetTitle("#phi");
249  hNewGains->GetYaxis()->CenterTitle();
250  hNewGains->SetStats(kFALSE);
251  hNewGains->Draw("colz");
252 
253  c->cd(2);
254  hChi2EtaPhi->SetTitle("#chi^{2}/dof");
255  hChi2EtaPhi->GetXaxis()->SetTitle("#eta");
256  hChi2EtaPhi->GetXaxis()->CenterTitle();
257  hChi2EtaPhi->GetYaxis()->SetTitle("#phi");
258  hChi2EtaPhi->GetYaxis()->CenterTitle();
259  hChi2EtaPhi->SetStats(kFALSE);
260  hChi2EtaPhi->Draw("colz");
261 
262  ps->Close();
263 
264  // Delete pointers to decoder & geometry
265  if (mDecoder) delete mDecoder;
266  if (mEmcGeom) delete mEmcGeom;
267 }
StEmcDecoder::GetCrateFromTowerId
int GetCrateFromTowerId(int softId, int &crate, int &sequence) const
Get crate number and position in crate for Software Id.
Definition: StEmcDecoder.cxx:556
Generated by   doxygen 1.8.5