StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
mip_ring_fitter.C
1 #include <iostream>
2 #include <fstream>
3 #include <set>
4 using namespace std;
5 
6 #include "CalibrationHelperFunctions.cxx"
7 
8 void mip_ring_fitter(const char* file_list="", const char* skimfile="mipskimfile.root")
9 {
10  gROOT->Macro("LoadLogger.C");
11  gROOT->Macro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
12  gSystem->Load("StTpcDb");
13  gSystem->Load("StDaqLib");
14  gSystem->Load("StDetectorDbMaker");
15  gSystem->Load("St_db_Maker");
16  gSystem->Load("StDbUtilities");
17  gSystem->Load("StEmcRawMaker");
18  gSystem->Load("StMcEvent");
19  gSystem->Load("StMcEventMaker");//***
20  gSystem->Load("StEmcSimulatorMaker");//***
21  gSystem->Load("StEmcADCtoEMaker");
22  gSystem->Load("StEpcMaker");
23  gSystem->Load("StDbBroker");
24  gSystem->Load("StEEmcUtil");
25  gSystem->Load("StAssociationMaker");
26  gSystem->Load("StEmcTriggerMaker");
27 
28  gSystem->Load("StEmcOfflineCalibrationMaker");
29 
30  const int ntowers=4800;
31  const int nrings=40;
32  const bool lookForSwaps=false;
33 
34  CalibrationHelperFunctions* helper = new CalibrationHelperFunctions();
35 
36  cout<<"input filelist: "<<file_list<<endl;
37  cout<<"histogram file: "<<skimfile<<endl;
38 
39  //chain all input files together
40  char file[300];
41  TChain* calib_tree = new TChain("calibTree");
42  ifstream filelist(file_list);
43  TFile *test_file;
44  while(1){
45  filelist >> file;
46  if(!filelist.good()) break;
47  cout<<file<<endl;
48  calib_tree->Add(file);
49  }
50 
51  StEmcOfflineCalibrationEvent* myEvent = new StEmcOfflineCalibrationEvent();
52  calib_tree->SetBranchAddress("event_branch",&myEvent);
53  StEmcOfflineCalibrationTrack* mip;
54 
55  //create the 4800 mip histograms
56  TH1* ring_histo[nrings];
57  char name[100];
58  for(int k=0; k<nrings; k++){
59  sprintf(name,"ring_histo_%i",k+1);
60  ring_histo[k] = new TH1D(name,name,250,-50.5,199.5);
61  }
62 
63  //keep track of all hit towers and exclude any with >1 track/tower
64  set<int> track_towers;
65  set<int> excluded_towers;
66 
67  unsigned int nentries = calib_tree->GetEntries();
68  for(unsigned int i=0; i<nentries; i++){
69  if(i%100000 == 0) cout<<"processing "<<i<<" of "<<nentries<<endl;
70 
71  track_towers.clear();
72  excluded_towers.clear();
73 
74  calib_tree->GetEntry(i);
75 
76  //we're only working with events from vertexIndex==0
77  if(TMath::Abs(myEvent->vz[0]) > 30.) continue;
78 
79  //do a quick loop over tracks to get excluded towers
80  for(int j=0; j<myEvent->tracks->GetEntries(); j++){
81  mip = (StEmcOfflineCalibrationTrack*)myEvent->tracks->At(j);
82  int id = mip->tower_id[0];
83 
84  if(track_towers.find(id) != track_towers.end()){
85  excluded_towers.insert(id);
86  }
87  else{
88  track_towers.insert(id);
89  }
90  }
91 
92  //select on runnumbers to look for stability
93  //if(myEvent->run < 7135000) continue;
94 
95  for(int j=0; j<myEvent->tracks->GetEntries(); j++){
96  mip = (StEmcOfflineCalibrationTrack*)myEvent->tracks->At(j);
97 
98 /* if(lookForSwaps){
99  for(int tower=1; tower<9; tower++){
100  double pedsub = mip->tower_adc[tower] - mip->tower_pedestal[tower];
101 
102  if(mip->p < 1.) continue;
103  if(mip->tower_status[tower] != 1) continue;
104  if(TMath::Abs(pedsub) < 1.5*mip->tower_pedestal_rms[tower]) continue;
105  if(mip->tower_id[0] != mip->tower_id_exit) continue;
106  if(mip->vertexIndex > 0) continue;
107 
108  int index = mip->tower_id[0];
109  mip_histo[index-1]->Fill(pedsub);
110  }
111  }
112  else
113  {
114  double pedsub = mip->tower_adc[0] - mip->tower_pedestal[0];
115 
116  if(mip->p < 1.) continue;
117  if(mip->tower_status[0] != 1) continue;
118  if(mip->highest_neighbor > 2.) continue;
119  if(TMath::Abs(pedsub) < 1.5*mip->tower_pedestal_rms[0]) continue;
120  if(mip->tower_id[0] != mip->tower_id_exit) continue;
121  if(mip->vertexIndex > 0) continue;
122  if(excluded_towers.find(mip->tower_id[0]) != excluded_towers.end()) continue;
123 
124 
125  int index = mip->tower_id[0];
126  mip_histo[index-1]->Fill(pedsub);
127  }*/
128 
129  //preshower histograms for Rory
130  //double pedsub = mip->preshower_adc[0] - mip->preshower_pedestal[0];
131  double pedsub = mip->tower_adc[0] - mip->tower_pedestal[0];
132  //double pedsub = mip->tower_adc[0];
133  if(mip->p < 1.)continue;
134  if(mip->tower_status[0] != 1)continue;
135  if(mip->highest_neighbor > 2.)continue;
136  if(TMath::Abs(pedsub) < 1.5*mip->tower_pedestal_rms[0])continue;
137  if(mip->tower_id[0] != mip->tower_id_exit)continue;
138  if(mip->vertexIndex > 0)continue;
139  if(excluded_towers.find(mip->tower_id[0]) != excluded_towers.end()) continue;
140 
141  //cout<<mip->tower_pedestal_rms[0]<<endl;
142  int index = mip->tower_id[0];
143  double eta = helper->getEta(index);
144  if(TMath::Abs(eta) > 0.968) eta += 0.005 * TMath::Abs(eta)/eta;
145  int etaindex = ((TMath::Nint(eta * 1000.0) + 25) / 50 + 19);
146  ring_histo[etaindex]->Fill(pedsub);
147 
148  }
149  }
150 
151  TFile* output_file = new TFile(skimfile,"RECREATE");
152  //for(int k=0; k<ntowers; k++) mip_histo[k]->Write();
153  for(int k=0; k<nrings;k++) ring_histo[k]->Write();
154  output_file->Close();
155  /*
156  //draw the histograms and their fits
157  TPostScript* ps = new TPostScript(postscript);
158  TCanvas* c = new TCanvas("c","",100,100,600.,800.);
159  int pad = 16;
160  for(int i=0; i<nrings; i++){
161  if(pad%15 == 1){
162  c->Update();
163  ps->NewPage();
164  c->Clear();
165  c->Divide(3,5);
166  pad = 1;
167  }
168  c->cd(pad);
169 
170  if(i%600 == 0) cout<<"fitting tower "<<i+1<<" of "<<ntowers<<endl;
171 
172  sprintf(name,"fit_%i",i+1);
173 
174  //this fit is for the electron tree
175 // fit[i] = new TF1(name,"gaus",0.,160.);
176 // fit[i] = new TF1(name,"gaus(0)+gaus(3)",0.,140.);
177 // fit[i]->SetParameter(1,65.);
178 // fit[i]->SetParameter(2,10.);
179 // fit[i]->SetParameter(4,10.);
180 // fit[i]->SetParameter(5,7.);
181  fit[i] = new TF1(name,fit_function,0.,140.,6);
182  fit[i]->SetParameter(1,65.);
183  fit[i]->SetParameter(2,10.);
184  fit[i]->SetParameter(3,10.); //relative height of peak to bg
185  fit[i]->SetParameter(4,10.);
186  fit[i]->SetParameter(5,3.);
187  fit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma");
188 
189  //this fit is for the hadron tree
190 // fit[i] = new TF1(name,"landau",0.,140.);
191 // fit[i] = new TF1(name,"gaus(0)+gaus(3)",0.,140.);
192 // fit[i] = new TF1(name,fitf,0.,140.,3);
193 // fit[i]->SetParameter(1,10.);
194 // fit[i]->SetParameter(2,3.);
195 // fit[i]->SetParameter(4,15.);
196 // fit[i]->SetParameter(5,25.);
197 
198  fit[i]->SetLineColor(kGreen);
199  fit[i]->SetLineWidth(0.6);
200 
201  ring_histo[i]->Fit(fit[i],"rq");
202 
203  drawTower(ring_histo[i],fit[i],i, helper);
204 
205  pad++;
206  }
207 
208  ps->Close();
209  */
210  //print gains to a file
211  // ofstream gains(gainfile);
212  //char line[500];
213  /* for(int i=0; i<nrings; i++){
214  float gain = 1./fit[i]->GetParameter(1);
215 // cout<<i+1<<" "<<gain<<" "<<helper->getTheta(i+1)<<" "<<TMath::Sin(helper->getTheta(i+1))<<endl;
216  gain /= TMath::Sin(helper->getTheta(i+1));
217  if(!helper->isGoodTower2006(i+1)) gain = 0.;
218  if(gain<0.) gain=0.;
219 
220  sprintf(line,"%-4i % 1.6f",i+1,gain);
221  gains<<line<<endl;
222  }*/
223 }
224 
225 void drawTower(TH1* h, TF1* f, int id, CalibrationHelperFunctions* helper){
226  //calculate a few quantities
227  double peak = f->GetParameter(1);
228  double histo_height = h->GetBinContent(h->GetMaximumBin());
229  if(histo_height == 0) histo_height = 1.;
230 
231  //histogram options
232  h->SetXTitle("ADC/gev*sin(#theta)");
233  h->Draw("esame");
234 
235  //draw a line through the location of the MIP peak
236  TLine *gaussian_peak = new TLine(peak,0.,peak,histo_height+15);
237  gaussian_peak->SetLineColor(kGreen);
238  gaussian_peak->SetLineWidth(2.0);
239  gaussian_peak->Draw("same");
240 
241  //write the tower number on the plot
242  char tower_title[100];
243  float eta = (float)((id - 20) * 2 + 1)/40;
244  sprintf(tower_title,"eta = %f",eta);
245  TLatex title_latex;
246  title_latex.SetTextSize(0.15);
247  if(!helper->isGoodTower2006(id)) title_latex.SetTextColor(kRed);
248  title_latex.DrawTextNDC(0.13,0.78,tower_title);
249 
250  //draw the peak Gaussian
251  TF1 *f2 = new TF1(tower_title,"gaus",0.,140.);
252  f2->FixParameter(0,f->GetParameter(0));
253  f2->FixParameter(1,f->GetParameter(1));
254  f2->FixParameter(2,f->GetParameter(2));
255  f2->SetLineWidth(0.6);
256  f2->Draw("same");
257 }
258 
259 //constrained double-Gaussian to fit the background
260 double background_only_fit(double *x, double *par){
261  double par3 = par[0]/1.5;
262  double par4 = par[1] + 10.;
263  double par5 = par[2] * 6.5;
264 
265  double fitval = 0;
266  if(par[2] != 0){
267  double arg1 = (x[0]-par[1])/par[2];
268  double arg2 = (x[0]-par4)/par5;
269  fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2);
270  }
271 
272  return fitval;
273 }
274 
275 double fit_function(double *x, double *par){
276  //6-parameter fit includes
277  //3 param electron Gaussian
278  //par3 = relative height of peak/bg ~ 10
279  //par4 = mean of main bg Gaussian ~ 10
280  //par5 = width of main bg Gaussian ~ 3
281 
282  double par3 = par[0] / par[3];
283  double par6 = par3/1.5;
284  double par7 = par[4] + 10.;
285  double par8 = par[5] * 6.5;
286 
287  double fitval = 0;
288  if(par[2] != 0){
289  double arg1 = (x[0]-par[1])/par[2];
290  double arg2 = (x[0]-par[4])/par[5];
291  double arg3 = (x[0]-par7)/par8;
292  fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2) + par6*TMath::Exp(-0.5*arg3*arg3);
293  }
294 
295  return fitval;
296 
297 
298 }
Generated by   doxygen 1.8.5