StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
electron_master.C
1 //copied from /star/u/staremc/2009/pp200/StRoot/StEmcPool/StEmcOfflineCalibrationMaker/macros/electron_master.C
2 
3 #include <iostream>
4 #include <fstream>
5 using namespace std;
6 
7 #include "TFile.h"
8 #include "TH1.h"
9 #include "TH2.h"
10 #include "TF1.h"
11 #include "TMath.h"
12 #include "TPostScript.h"
13 #include "TCanvas.h"
14 #include "TStyle.h"
15 #include "TLatex.h"
16 #include "TTree.h"
17 #include "TLine.h"
18 
19 #include "CalibrationHelperFunctions.h"
20 #include "CalibrationHelperFunctions.cxx"
21 
22 void drawTower(TH1* h, TF1* f, int id, CalibrationHelperFunctions* helper);
23 double fit_function(double *x, double *par);
24 double fit_function2(double *x, double *par);
25 double background_only_fit(double *x, double *par);
26 int lookup_crate(float,float);
27 
28 //void electron_master(const char* file_list="checkP10ic/electrons.list",const char* output="checkP10ic/electroncalib.root", const char* dbDate="2009-04-10 00:11:00", const char* geantfile="geant_func.root", const char* gfname="checkP10ic/mip.gains", const char* ngname="checkP10ic/electron.gains", const char* rgname="matrel.gains",const char* ffname="checkP10ic/electron.fits"){
29 void electron_master(const char* file_list="checkP11id/electrons.list",const char* output="checkP11id/electroncalib2.root", const char* dbDate="2009-04-10 00:11:00", const char* gfname="checkP11id/mip.gains", const char* ngname="checkP11id/electron2.gains",const char* ffname="checkP11id/electron2.fits"){
30  //**********************************************//
31  //Load Libraries //
32  //**********************************************//
33  gROOT->Macro("LoadLogger.C");
34  gROOT->Macro("loadMuDst.C");
35  gSystem->Load("StTpcDb");
36  gSystem->Load("StDaqLib");
37  gSystem->Load("StDetectorDbMaker");
38  gSystem->Load("St_db_Maker");
39  gSystem->Load("StDbUtilities");
40  gSystem->Load("StEmcRawMaker");
41  gSystem->Load("StMcEvent");
42  gSystem->Load("StMcEventMaker");//***
43  gSystem->Load("StEmcSimulatorMaker");//***
44  gSystem->Load("StEmcADCtoEMaker");
45  gSystem->Load("StEpcMaker");
46  gSystem->Load("StDbBroker");
47  gSystem->Load("StEEmcUtil");
48  gSystem->Load("StAssociationMaker");
49  gSystem->Load("StEmcTriggerMaker");
50  gSystem->Load("StTriggerUtilities");
51  gSystem->Load("StEmcOfflineCalibrationMaker");
52 
53  cout<<"input filelist: "<<file_list<<endl;
54  cout<<"output: "<<output<<endl;
55  cout<<"db Date: "<<dbDate<<endl;
56 
57 
58  //**********************************************//
59  //Initialize Stuff //
60  //**********************************************//
61 
62  CalibrationHelperFunctions *helper = new CalibrationHelperFunctions();
63 
64  StBemcTablesWriter *bemctables = new StBemcTablesWriter();
65  bemctables->loadTables(dbDate,"sim");
66  StEmcDecoder* decoder = bemctables->getDecoder();
67 
68 
69  //chain all input files together
70  char file[300];
71  TChain* tree = new TChain("skimTree");
72  ifstream filelist(file_list);
73  while(1){
74  filelist >> file;
75  if(!filelist.good()) break;
76  cout<<file<<endl;
77  tree->Add(file);
78  }
79 
80  const int ntowers = 4800;
81  const int nrings = 40;
82  const int ncrates = 30;
83  float gains[ntowers];
84  int status[ntowers];
85  float peaks[ntowers];
86  float peakerr[ntowers];
87  float gainerr[ntowers];
88 
89  ifstream gainfile(gfname);
90  while(1){
91  int id,stat;
92  float peak,err,gain,eta,theta;
93  gainfile >> id >> peak >> err >> stat;
94  if(!gainfile.good())break;
95  eta = helper->getEta(id);
96  theta = helper->getTheta(id);
97  gain = 0.264*(1+0.056*eta*eta)/(sin(theta)*peak);
98  peaks[id-1] = peak;
99  gains[id-1] = gain;
100  if(status[id-1] != 1 || stat != 1) status[id-1] = stat;
101  peakerr[id-1] = err;
102  //cout<<id<<" "<<gain;
103  }
104 
105 
106  TFile* geant_file = new TFile("geant_fits.root","READ");
107  TF1* geant_fits[20];
108  for(int i = 0; i < 20; i++){
109  TString fname = "fit_";
110  fname += i;
111  geant_fits[i] = (TF1*)geant_file->Get(fname);
112  }
113 
114 
115  TFile outfile(output,"RECREATE");
116 
117  TH1 *crate_histo[ncrates];
118  TF1 *crate_fit[ncrates];
119  TH1 *electron_histo[ntowers];
120  TF1 *fit[ntowers];
121  TH1 *prs_histo[ntowers];
122  TH1 *ring_histo[nrings];
123  TH1 *ring2_histo[nrings/2];
124  TF1 *ring2fit[nrings/2];
125 
126  TH1 *etacrate_histo[ncrates*20];
127  TH1 *etacrate_histo_p[ncrates*20];
128  TH1 *etacrate_histo_n[ncrates*20];
129 
130  TH1 *bghisto[21];
131  TH2F* drvsep = new TH2F("drvsep","drvsep",60,0,3.0,100,0.0,0.03);
132 
133  TF1 *ringfit[nrings];
134  TH2F* ring_pve[nrings];
135  TH2F* jan_pve[6];
136  float eta;
137  int etaindex;
138  TH1F* ringprec = new TH1F("ringprec","",40,-1.0,1.0);
139  TH1F* ringprec2 = new TH1F("ringprec2","",40,-1.0,1.0);
140  ringprec->SetYTitle("E/p");
141  ringprec->SetXTitle("#eta");
142  ringprec2->SetYTitle("E/p");
143  ringprec2->SetXTitle("#eta");
144  TH1F* ring2prec = new TH1F("ring2prec","",nrings/2,-1.0,1.0);
145  TH1F* ring2prec2 = new TH1F("ring2prec2","",nrings/2,-1.0,1.0);
146  ring2prec->SetYTitle("E/p");
147  ring2prec->SetXTitle("#eta");
148  ring2prec2->SetYTitle("E/p");
149  ring2prec2->SetXTitle("#eta");
150  double ew[nrings];
151 
152  TH2F* trgcheck = new TH2F("trgcheck","trgcheck",2,-0.5,1.5,3,-0.5,2.5);
153 
154  TH1F* crateprec = new TH1F("crateprec","",30,0.5,30.5);
155  crateprec->SetXTitle("Crate");
156  crateprec->SetYTitle("E/p");
157 
158  TH2F *energyleak = new TH2F("energyleak","",20,0.0,0.03,20,0.0,1.0);
159  TH2F *findbg = new TH2F("findbg","",20,0.0,0.03,30,0.0,5.0);
160  TH1F *energymean = new TH1F("energymean","",20,0.0,0.03);
161  TH1F *leakmean = new TH1F("leakmean","",20,0.0,0.03);
162  energyleak->SetXTitle("#DeltaR");
163  energyleak->SetYTitle("leaked energy / total energy");
164  findbg->SetXTitle("#DeltaR");
165  findbg->SetYTitle("Total energy / track momentum");
166  TH2F *tevsp = new TH2F("tevsp","",50,0.0,20.0,50,0.0,30.0);
167  TH1F *pmean = new TH1F("pmean","",20,0.0,15.0);
168  tevsp->SetXTitle("track momentum (GeV)");
169  tevsp->SetYTitle("Total Energy (GeV)");
170  TH2F *tevspcent = new TH2F("tevspcent","",20,0.0,15.0,20,0.0,15.0);
171  tevspcent->SetXTitle("track momentum (GeV)");
172  tevspcent->SetYTitle("Energy in central tower (GeV)");
173  TH1F *cmean = new TH1F("cmain","",20,0.0,15.0);
174  TH2F *sistertracks = new TH2F("sistertracks","",20,0.0,8.0,20,0.0,8.0);
175  sistertracks->SetXTitle("Track momentum (GeV)");
176  sistertracks->SetYTitle("Neighbor momentum (GeV)");
177  TH2F* dEdxvsp = new TH2F("dEdxvsp","",100,0.15,1.3,100,-5.7,-5.);
178  dEdxvsp->SetXTitle("Log(p)");
179  dEdxvsp->SetYTitle("Log(dE/dx)");
180  TH2F* dEdxvsp_east = new TH2F("dEdxvsp_east","",100,0.15,1.3,100,-5.7,-5.0);
181  dEdxvsp_east->SetXTitle("Log(p)");
182  dEdxvsp_east->SetYTitle("Log(dE/dx)");
183  TH2F* dEdxvsp_west = new TH2F("dEdxvsp_west","",100,0.15,1.3,100,-5.7,-5.0);
184  dEdxvsp_west->SetXTitle("Log(p)");
185  dEdxvsp_west->SetYTitle("Log(dE/dx)");
186  TH2F* energyleak2 = new TH2F("energyleak2","",20,0.0,0.03,20,0.0,1.0);
187  TH1F* energymean2 = new TH1F("energymean2","",20,0.0,0.03);
188  TH1F* towermult = new TH1F("towermult","",9,0.0,9.0);
189  energyleak2->SetXTitle("#DeltaR");
190  energyleak2->SetYTitle("leaked energy/total energy");
191  towermult->SetXTitle("Neighbors with energy");
192  TH2F* multvsp = new TH2F("multvsp","",20,0.0,20.0,9,0.0,9.0);
193  multvsp->SetXTitle("Track momentum (GeV)");
194  multvsp->SetYTitle("Neighbors with energy");
195  TH1F* multmean = new TH1F("multmean","",20,0.0,20.0);
196  TH2F* tep3 = new TH2F("tep3","2 < p < 3",20,0.0,0.03,20,0.0,4.0);
197  TH2F* tep5 = new TH2F("tep5","3 < p < 5",20,0.0,0.03,20,0.0,4.0);
198  TH2F* tep10 = new TH2F("tep10","5 < p < 10",20,0.0,0.03,20,0.0,4.0);
199  tep3->SetXTitle("DeltaR");
200  tep3->SetYTitle("Total energy / track momentum");
201  tep5->SetXTitle("DeltaR");
202  tep5->SetYTitle("Total energy / track momentum");
203  tep10->SetXTitle("DeltaR");
204  tep10->SetYTitle("Total energy / track momentum");
205  TH1F* tep3mean = new TH1F("tep3mean","",20,0,0.03);
206  TH1F* tep5mean = new TH1F("tep5mean","",20,0,0.03);
207  TH1F* tep10mean = new TH1F("tep10mean","",20,0,0.03);
208  TH1F* multen = new TH1F("multen","",40,0.0,1.0);
209  multen->SetXTitle("Energy in neighbor towers (GeV)");
210  TH1F* east_histo = new TH1F("east_histo","Electron E/p in East",40,0.0,2.0);
211  TH1F* west_histo = new TH1F("west_histo","Electron E/p in West",40,0.0,2.0);
212  TH1F* all_histo = new TH1F("all_histo","Electron E/p",40,0.0,2.0);
213  TH1F* notrg = new TH1F("notrg","Electron E/p",40,0.0,2.0);
214  TH2F* pvsep = new TH2F("pvsep","Electron p vs E/p",120,0,3.0,20,0,20.0);
215  pvsep->SetYTitle("p (Gev)");
216  pvsep->SetXTitle("E/p");
217  TH2F* pvsep0 = new TH2F("pvsep0","Electron p vs E/p",120,0,3.0,20,0,20.0);
218  pvsep0->SetYTitle("p (Gev)");
219  pvsep0->SetXTitle("E/p");
220  TH2F* evsep = new TH2F("evsep","Electron E vs E/p",120,0,3.0,20,0,20.0);
221  evsep->SetYTitle("E (GeV)");
222  evsep->SetXTitle("E/p");
223 
224  TH1F* bsmde = new TH1F("bsmde","BSMDE ADC TOT",1500,0,1500);
225  TH1F* bsmdp = new TH1F("bsmdp","BSMDP ADC TOT",1500,0,1500);
226 
227  TH1F* bsmde_central = new TH1F("bsmde_central","BSMDE ADC TOT",100,0,1500);
228  TH1F* bsmde_mid = new TH1F("bsmde_mid","BSMDE ADC TOT",100,0,1500);
229  TH1F* bsmde_outer = new TH1F("bsmde_outer","BSMDE ADC TOT",100,0,1500);
230 
231  TH2F* bsmdep = new TH2F("bsmdep","BSMDE v BSMDP",100,0,1500,100,0,1500);
232 
233  TH2F* bsmdevp = new TH2F("bsmdevp","BSMDE v p",100,1.0,15.0,100,0,1500);
234  TH2F* bsmdpvp = new TH2F("bsmdpvp","BSMDP v p",100,1.0,15.0,100,0,1500);
235  TH2F* bsmdevep = new TH2F("bsmdevep","BSMDE v E/p",100,0.0,2.0,100,0,1500);
236  TH2F* bsmdpvep = new TH2F("bsmdpvep","BSMDP v E/p",100,0.0,2.0,100,0,1500);
237 
238  TH2F* bsmdeve = new TH2F("bsmdeve","BSMDE v E",100,1.0,30.0,100,0,1500);
239  TH2F* bsmdpve = new TH2F("bsmdpve","BSMDP v E",100,1.0,30.0,100,0,1500);
240 
241  TH1F* httrig = new TH1F("httrig","HT Trigger",5,-0.5,4.5);
242 
243  TH1F* pplus = new TH1F("pplus","e+ p",100,0,20);
244  TH1F* pminus = new TH1F("pminus","e- p",100,0,20);
245 
246  TH1F* posep = new TH1F("posep","e+ E/p",60,0,3.0);
247  TH1F* negep = new TH1F("negep","e- E/p",60,0,3.0);
248 
249 
250  //create the tower histograms
251  char name1[100];
252 
253  for(int i = 0; i < ncrates; i++){
254  for(int j = 0; j < 20; j++){
255  TString ecname,ecnamep,ecnamen;
256  ecname += "etacrate_";
257  int cr = i+1;
258  int et = j;
259  ecname += cr;
260  ecname += "_";
261  ecname += et;
262  ecnamep += ecname;
263  ecnamen += ecname;
264  ecnamep.ReplaceAll("te_","te_p_");
265  ecnamen.ReplaceAll("te_","te_n_");
266  etacrate_histo[i*20+j] = new TH1F(ecname.Data(),ecname.Data(),60,0.0,3.0);
267  etacrate_histo_p[i*20+j] = new TH1F(ecnamep.Data(),ecnamep.Data(),60,0.0,3.0);
268  etacrate_histo_n[i*20+j] = new TH1F(ecnamen.Data(),ecnamen.Data(),60,0.0,3.0);
269 
270  }
271  }
272 
273  for(int i=0; i<ntowers; i++){
274  char nameeh[100];
275  sprintf(nameeh,"electron_histo_%i",i+1);
276  electron_histo[i] = new TH1D(nameeh,"",60,0.,3.0);
277  electron_histo[i]->SetXTitle("E/p");
278  //sprintf(name1,"prs_histo_%i",i+1);
279  //prs_histo[i] = new TH1D(name1,"",60,0.,500.);
280  //prs_histo[i]->SetXTitle("ADC");
281  }
282 
283  for(int i = 0; i < 21; i++){
284  TString namebg;
285  namebg += "bg_";
286  namebg += i+25;
287  bghisto[i] = new TH1F(namebg.Data(),namebg.Data(),60,0,3.0);
288  }
289 
290  for(int i = 0; i < ncrates; i++){
291  char name[100];
292  char title[100];
293  sprintf(name,"crate_%i",i+1);
294  sprintf(title,"E/p for Crate %i",i+1);
295  crate_histo[i] = new TH1F(name,title,60,0.,3.0);
296  crate_histo[i]->SetXTitle("E/p");
297  }
298 
299  for(int i = 0; i < nrings/2; i++){
300  char name[100];
301  sprintf(name,"ring2_histo_%i",i);
302  ring2_histo[i] = new TH1F(name,"",60,0.,3.0);
303  ring2_histo[i]->SetXTitle("E/p");
304  }
305 
306  char name2[100];
307  for(int i=0; i<nrings;i++)
308  {
309  sprintf(name2,"ring_histo_%i",i);
310  //ring_histo[i] = new TH1D(name2,"",30,0.,140.0);
311  //ring_histo[i]->SetXTitle("ADC / GeV Sin(#theta)");
312  ring_histo[i] = new TH1D(name2,"",60,0.,3.0);
313  ring_histo[i]->SetXTitle("E/p");
314  char namerpve[100];
315  sprintf(namerpve,"ring_pve_%i",i);
316  ring_pve[i] = new TH2F(namerpve,"",20,0,20.0,20,0,20.0);
317  ring_pve[i]->SetXTitle("E (GeV)");
318  ring_pve[i]->SetYTitle("p (GeV)");
319  }
320 
321  for(int i = 0; i < 6; i++){
322  char jname[100];
323  sprintf(jname,"jan_pve_%i",i);
324  jan_pve[i] = new TH2F(jname,"",120,0,3.0,20,0,20.0);
325  jan_pve[i]->SetXTitle("E/p");
326  jan_pve[i]->SetYTitle("p (GeV)");
327  }
328  //global graphics functions
329  gStyle->SetOptStat("oue");
330  gStyle->SetOptFit(111);
331  gStyle->SetCanvasColor(10);
332  gStyle->SetCanvasBorderMode(0);
333  gStyle->SetOptTitle(0);
334  gStyle->SetPalette(1);
335  gStyle->SetStatColor(0);
336 
337  StEmcOfflineCalibrationTrack* track = new StEmcOfflineCalibrationTrack();
338  StEmcOfflineCalibrationCluster* cluster = new StEmcOfflineCalibrationCluster();
339  tree->SetBranchAddress("clusters",&cluster);
340 
341  //**********************************************//
342  //Loop Over Tracks, Fill Histograms //
343  //**********************************************//
344 
345  int nentries = tree->GetEntries();
346  cout<<nentries<<endl;
347  int ngoodhit = 0;
348  int nplt10 = 0;
349  int nnosis = 0;
350  int nfinal = 0;
351  int nbsmdgood = 0;
352  int nnottrig = 0;
353  int nfidu = 0;
354  int nenterexit = 0;
355  for(int j=0; j<nentries; j++){
356  tree->GetEntry(j);
357  track = &(cluster->centralTrack);
358  TClonesArray *tracks = cluster->tracks;
359  if(j % 500000 == 0) cout<<"reading "<<j<<" of "<<nentries<<endl;
360 
361  httrig->Fill((float)track->htTrig);
362 
363  if(track->charge > 0)pplus->Fill(track->p);
364  if(track->charge < 0)pminus->Fill(track->p);
365 
366  int bsmdeadctot = 0;
367  int bsmdpadctot = 0;
368  for(int i = 0; i < 11; i++){
369  if(track->smde_adc[i] > track->smde_pedestal[i])bsmdeadctot += track->smde_adc[i] - track->smde_pedestal[i];
370  if(track->smdp_adc[i] > track->smdp_pedestal[i])bsmdpadctot += track->smdp_adc[i] - track->smdp_pedestal[i];
371  }
372 
373  double dR = TMath::Sqrt(track->deta*track->deta + track->dphi*track->dphi);
374 
375 
376  double scaled_adc = (track->tower_adc[0] - track->tower_pedestal[0]) / track->p;
377 
378  int index = track->tower_id[0];
379 
380  //figure out eta and etaindex
381  eta = helper->getEta(index);
382  if(TMath::Abs(eta) > 0.968) eta += 0.005 * TMath::Abs(eta)/eta;
383  etaindex = ((TMath::Nint(eta * 1000.0) + 25)/50 + 19);
384  int geantetaindex = ((TMath::Nint(fabs(eta) * 1000.0) + 25)/50 - 1);
385 
386  double geant_scale = geant_fits[geantetaindex]->Eval(dR);
387  scaled_adc /= geant_scale;
388  //double geant_scale = geant_fit->Eval(dR);
389  //scaled_adc *= geant_scale;
390  //cout<<scaled_adc<<endl;
391 
392  //now rescale dR for last ring to make cuts work
393  if(geantetaindex == 19)dR *= 0.025/0.017;
394 
395  //double tgain = bemctables->calib(1,track->tower_id[0])*gains[index-1];
396  double tgain = gains[index-1];
397 
398 
399  if((track->tower_adc[0] - track->tower_pedestal[0]) < 2.5 * track->tower_pedestal_rms[0])continue;
400  ngoodhit++;
401  dEdxvsp->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
402  if(track->tower_id[0] <= 2400)dEdxvsp_west->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
403  if(track->tower_id[0] > 2400)dEdxvsp_east->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
404 
405  trgcheck->Fill(track->nonhtTrig,track->htTrig);
406 
407  if(track->tower_id[0] != track->tower_id_exit)continue;
408  nenterexit++;
409 
410  if(status[index-1]!=1)continue;
411 
412  pvsep0->Fill(scaled_adc*tgain,track->p);
413  if(track->p > 10)continue;
414  nplt10++;
415 
416  //if(track->p < 1.5)continue;
417  //if(track->p < 3.0)continue;
418  if(track->htTrig == 2 && track->nonhtTrig == 0)continue;
419  nnottrig++;
420  //change the fiducial cut to a square with diagonal = 0.06 in deta, dphi space
421  float squarefid = 0.02;//0.03/TMath::Sqrt(2.0);
422  //if(TMath::Abs(track->deta) > squarefid || TMath::Abs(track->dphi) > squarefid)continue;
423 
424  //calculate geant scaled, pedestal subtracted adc
425  //if(dR > 0.0125)continue;
426  int numsis = tracks->GetEntries();
427  float totalbtow = 0;
428  float maxEt = 0;
429  int maxId = -1;
430  for(int i = 0; i < 9; i++){
431  if(track->tower_adc[i] - track->tower_pedestal[i] < 0)continue;
432  float theta = helper->getTheta(track->tower_id[i]);
433  float nextEt = (track->tower_adc[i] - track->tower_pedestal[i]) * bemctables->calib(1,track->tower_id[i])*sin(theta);
434  totalbtow += nextEt;
435  if(nextEt > maxEt){
436  maxEt = nextEt;
437  maxId = i;
438  }
439  }
440 
441  if(track->dEdx > 3.5e-6 && track->dEdx <5.0e-6 && numsis ==0 &&maxId == 0)drvsep->Fill(scaled_adc*tgain,dR);
442  if(dR > 0.02)continue;
443  nfidu++;
444  //if(track->p > 6.0)continue;
445  //if(track->p > 15)continue;
446  //cout<<track->dEdx<<endl;
447 
448 
449  //if(track->dEdx*1000000 > 4.5 || track->dEdx*1000000 < 3.5)continue;
450 
451  //cout<<track->htTrig<<endl;
452 
453  for(int i = 0; i < 21; i++){
454  if(numsis > 0)break;
455  if(maxId != 0)break;
456  if(track->dEdx*1e7 > 25 + i && track->dEdx*1e7 < 26+i)bghisto[i]->Fill(scaled_adc*tgain);
457  }
458 
459  //if(track->dEdx < 3.5e-6 || track->dEdx > 5.0e-6)continue;
460  if(track->dEdx < 3.5e-6 || track->dEdx > 4.5e-6)continue;
461  //if((bsmdeadctot > -1 && bsmdeadctot < 50) && (bsmdpadctot < 50 && bsmdpadctot > -1))continue;
462  nbsmdgood++;
463  //if(bsmdeadctot < 500) continue;
464 
465  //if(bsmdeadctot < 84.*track->p)continue;
466  //if(bsmdeadctot > 200.*track->p + 1500)continue;
467 
468  //if(bsmdpadctot < 800)continue;
469 
470  if(numsis > 0)continue;
471 
472  nnosis++;
473  if(maxId != 0) continue;
474 
475 
476  nfinal++;
477  if(track->htTrig!=2)notrg->Fill(scaled_adc*tgain);
478 
479  //if(!track->nonhtTrig)continue;
480  //if(track->nHits < 25)continue;
481  //if(status[index-1]==1)ring_histo[etaindex]->Fill(scaled_adc*gains[index-1]);
482 
483  //scaled_adc = totalbtow/track->p;
484  //tgain = 1.0;
485 
486  float phi = helper->getPhi(index);
487  int crate,sequence;
488  decoder->GetCrateFromTowerId(index,crate,sequence);
489  etacrate_histo[(crate-1)*20+geantetaindex]->Fill(scaled_adc*tgain);
490  if(track->charge > 0)etacrate_histo_p[(crate-1)*20+geantetaindex]->Fill(scaled_adc*tgain);
491  if(track->charge < 0)etacrate_histo_n[(crate-1)*20+geantetaindex]->Fill(scaled_adc*tgain);
492 
493  electron_histo[index-1]->Fill(scaled_adc*tgain);
494  ring_histo[etaindex]->Fill(scaled_adc*tgain);
495  if(etaindex == 0 || etaindex == 39){
496  //cout<<etaindex<<" "<<tgain<<" "<<track->p<<" "<<scaled_adc*tgain*track->p<<" "<<scaled_adc*tgain<<endl;
497  }
498  //cout<<index<<" "<<gains[index-1]<<" "<<scaled_adc*track->p<<" "<<track->p<<" "<<status[index-1]<<endl;
499  ring_pve[etaindex]->Fill(scaled_adc*tgain*track->p,track->p);
500 
501  dEdxvsp->Fill(TMath::Log10(track->p),TMath::Log10(track->dEdx));
502 
503  float abseta = TMath::Abs(eta);
504  if(abseta > 0.95){
505  jan_pve[5]->Fill(scaled_adc*tgain,track->p);
506  }else if(abseta > 0.9){
507  jan_pve[4]->Fill(scaled_adc*tgain,track->p);
508  }else if(abseta > 0.6){
509  jan_pve[3]->Fill(scaled_adc*tgain,track->p);
510  }else if(abseta > 0.3){
511  jan_pve[2]->Fill(scaled_adc*tgain,track->p);
512  }else if(abseta > 0.05){
513  jan_pve[1]->Fill(scaled_adc*tgain,track->p);
514  }else{
515  jan_pve[0]->Fill(scaled_adc*tgain,track->p);
516  }
517 
518  all_histo->Fill(scaled_adc*tgain);
519  pvsep->Fill(scaled_adc*tgain,track->p);
520  evsep->Fill(scaled_adc*tgain,track->p*scaled_adc*tgain);
521  tevsp->Fill(track->p,scaled_adc*tgain*track->p);
522 
523  if(track->charge > 0)posep->Fill(scaled_adc*tgain);
524  if(track->charge < 0)negep->Fill(scaled_adc*tgain);
525 
526  //if(scaled_adc*tgain < 0.7 || scaled_adc*tgain > 5.0)continue;
527  bsmde->Fill(bsmdeadctot);
528  bsmdp->Fill(bsmdpadctot);
529  bsmdep->Fill(bsmdeadctot,bsmdpadctot);
530 
531  if(abseta > 0.6){
532  bsmde_outer->Fill(bsmdeadctot);
533  }else if(abseta > 0.3){
534  bsmde_mid->Fill(bsmdeadctot);
535  }else{
536  bsmde_central->Fill(bsmdeadctot);
537  }
538 
539  bsmdevp->Fill(track->p,bsmdeadctot);
540  bsmdpvp->Fill(track->p,bsmdpadctot);
541  bsmdevep->Fill(scaled_adc*tgain,bsmdeadctot);
542  bsmdpvep->Fill(scaled_adc*tgain,bsmdpadctot);
543 
544  bsmdeve->Fill(scaled_adc*tgain*track->p,bsmdeadctot);
545  bsmdpve->Fill(scaled_adc*tgain*track->p,bsmdpadctot);
546 
547  if(dR > 0.015)continue;
548  if(track->tower_id[0] <= 2400){
549  west_histo->Fill(scaled_adc*tgain);
550  }
551  if(track->tower_id[0] > 2400){
552  east_histo->Fill(scaled_adc*tgain);
553  }
554 
555  }
556  cout<<"processed electron tree"<<endl;
557  cout<<"ngoodhit: "<<ngoodhit<<endl;
558  cout<<"nenterexit: "<<nenterexit<<endl;
559  cout<<"n not trig: "<<nnottrig<<endl;
560  cout<<"n p < 10: "<<nplt10<<endl;
561  cout<<"n in fidu: "<<nfidu<<endl;
562  cout<<"nbsmdhit: "<<nbsmdgood<<endl;
563  cout<<"n no sis: "<<nnosis<<endl;
564  cout<<"n final: "<<nfinal<<endl;
565 
566  //double ew[21];
567  for(int h=0;h<21;h++){
568  TH1D* projection = energyleak->ProjectionY("projection",h,h);
569  float mean = projection->GetMean();
570  energymean->SetBinContent(h,mean);
571  TH1D* projection1 = findbg->ProjectionY("projection1",h,h);
572  float mean1 = projection1->GetMean();
573  leakmean->SetBinContent(h,mean1);
574  TH1D* projection2 = tevsp->ProjectionY("projection2",h,h);
575  float mean2 = projection2->GetMean();
576  pmean->SetBinContent(h,mean2);
577  //ew[h] = projection2->GetRMS();
578  TH1D* projection3 = tevspcent->ProjectionY("projection3",h,h);
579  float mean3 = projection3->GetMean();
580  cmean->SetBinContent(h,mean3);
581  TH1D* projection4 = energyleak2->ProjectionY("projection4",h,h);
582  float mean4 = projection4->GetMean();
583  energymean2->SetBinContent(h,mean4);
584  TH1D* projection5 = multvsp->ProjectionY("projection5",h,h);
585  float mean5 = projection5->GetMean();
586  multmean->SetBinContent(h,mean5);
587  TH1D* projection6 = tep3->ProjectionY("projection6",h,h);
588  float mean6 = projection6->GetMean();
589  tep3mean->SetBinContent(h,mean6);
590  TH1D* projection7 = tep3->ProjectionY("projection7",h,h);
591  float mean7 = projection7->GetMean();
592  tep5mean->SetBinContent(h,mean7);
593  TH1D* projection8 = tep3->ProjectionY("projection8",h,h);
594  float mean8 = projection8->GetMean();
595  tep10mean->SetBinContent(h,mean8);
596  }
597 
598  TF1* fitleak = new TF1("fitleak","[0]",0,0.03);
599  fitleak->SetLineWidth(0.1);
600  leakmean->Fit(fitleak,"rq");
601 
602 
603  //**********************************************//
604  //Fit Tower Histograms //
605  //**********************************************//
606  /*
607  for(int i=0; i<ntowers; i++){
608 
609  if(i%600 == 0) cout<<"fitting tower "<<i+1<<" of "<<ntowers<<endl;
610 
611  sprintf(name,"fit_%i",i+1);
612 
613  //this fit is for the electron tree
614  fit[i] = new TF1(name,fit_function,0.,140.,6);
615  fit[i]->SetParameter(1,65.);
616  fit[i]->SetParameter(2,10.);
617  fit[i]->SetParameter(3,10.); //relative height of peak to bg
618  fit[i]->SetParameter(4,10.);
619  fit[i]->SetParameter(5,3.);
620  fit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma");
621 
622  fit[i]->SetLineColor(kGreen);
623  fit[i]->SetLineWidth(0.6);
624 
625  electron_histo[i]->Fit(fit[i],"rq");
626  }
627  */
628  //**********************************************//
629  //Fit Ring Histograms //
630  //**********************************************//
631 
632  for(int i=0; i<nrings; i++){
633 
634  //cout<<"fitting ring "<<i+1<<" of "<<nrings<<endl;
635 
636  sprintf(name,"ring_fit_%i",i);
637  /*
638  ringfit[i] = new TF1(name,fit_function,0.,140.,6);
639  ringfit[i]->SetParameter(1,1.);
640  ringfit[i]->SetParameter(2,0.2);
641  ringfit[i]->SetParameter(3,1.5); //relative height of peak to bg
642  ringfit[i]->SetParameter(4,0.15);
643  ringfit[i]->SetParameter(5,0.8);
644  ringfit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma");
645  */
646  //TF1* ffff = new TF1("ffff","expo(0) + gaus(2)",0.4,1.7);
647 
648  ring_histo[i]->Sumw2();
649  //ring_histo[i]->Rebin(3);
650 
651  ringfit[i] = new TF1(name,"pol1(0) + gaus(2)");
652  ringfit[i]->SetParLimits(0,0,10.0*ring_histo[i]->GetBinContent(1));
653  ringfit[i]->SetParLimits(1,-10000,0);
654  ringfit[i]->SetParLimits(2,0,10.0*ring_histo[i]->GetMaximum());
655  ringfit[i]->SetParLimits(3,0,10);
656  ringfit[i]->SetParameter(0,ring_histo[i]->GetBinContent(1));
657  ringfit[i]->SetParameter(1,-ring_histo[i]->GetBinContent(1)/6.0);
658  ringfit[i]->SetParameter(2,ring_histo[i]->GetMaximum());
659  ringfit[i]->SetParameter(3,0.95);
660  ringfit[i]->SetParameter(4,0.15);
661  ringfit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
662  /*
663  ringfit[i] = new TF1(name,fit_function2,0.1,2.5,8);
664  ringfit[i]->SetParameter(1,1.);
665  ringfit[i]->SetParameter(2,0.2);
666  ringfit[i]->SetParameter(3,1.5); //relative height of peak to bg
667  ringfit[i]->SetParameter(4,0.25);
668  ringfit[i]->SetParameter(5,0.15);
669  ringfit[i]->SetParameter(7,0.8);
670  ringfit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma","Bg2 constant","Bg2 decay");
671 
672  */
673  ringfit[i]->SetLineColor(kBlue);
674  ringfit[i]->SetLineWidth(0.6);
675 
676  ring_histo[i]->Fit(ringfit[i],"rql","",0.2,1.7);
677 
678  ringprec->SetBinContent(i+1,(ringfit[i]->GetParameter(3)));
679  ringprec->SetBinError(i+1,ringfit[i]->GetParameter(4));
680  ringprec2->SetBinContent(i+1,(ringfit[i]->GetParameter(3)));
681  ringprec2->SetBinError(i+1,ringfit[i]->GetParError(3));
682  //ew[i] = 4066/(60*(fit[i]->GetParameter(2))*(fit[i]->GetParameter(2)));
683 
684  float mean = ringfit[i]->GetParameter(3);
685  float merr = ringfit[i]->GetParError(3);
686  cout<<"ring "<<i<<" "<<mean<<" "<<merr/mean<<" "<<ring_histo[i]->GetEntries()<<endl;
687  }
688 
689  for(int i = 0; i < nrings/2; i++){
690  ring2_histo[i]->Add(ring_histo[2*i]);
691  ring2_histo[i]->Add(ring_histo[2*i+1]);
692  sprintf(name,"ring2_fit_%i",i);
693  ring2fit[i] = new TF1(name,"pol1(0) + gaus(2)",0.3,1.7);
694 
695  ring2_histo[i]->Rebin(3);
696 
697  ring2fit[i]->SetParLimits(0,0,10.0*ring2_histo[i]->GetBinContent(1));
698  ring2fit[i]->SetParLimits(1,-10000,0);
699  ring2fit[i]->SetParLimits(2,0,10.0*ring2_histo[i]->GetMaximum());
700  ring2fit[i]->SetParLimits(3,0,10);
701  ring2fit[i]->SetParLimits(4,0.17,0.175);
702  ring2fit[i]->SetParameter(0,ring2_histo[i]->GetBinContent(1));
703  ring2fit[i]->SetParameter(1,-ring2_histo[i]->GetBinContent(1)/6.0);
704  ring2fit[i]->SetParameter(2,ring2_histo[i]->GetMaximum());
705  ring2fit[i]->SetParameter(3,0.95);
706  ring2fit[i]->SetParameter(4,0.11245);
707  ring2fit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
708 
709  ring2_histo[i]->Fit(ring2fit[i],"rql","",0.3,1.7);
710 
711  ring2prec->SetBinContent(i+1,(ring2fit[i]->GetParameter(3)));
712  ring2prec->SetBinError(i+1,ring2fit[i]->GetParameter(4));
713  ring2prec2->SetBinContent(i+1,(ring2fit[i]->GetParameter(3)));
714  ring2prec2->SetBinError(i+1,ring2fit[i]->GetParError(3));
715 
716  cout<<"ring2 "<<i<<" "<<ring2fit[i]->GetParameter(3)<<" "<<ring2fit[i]->GetParError(3)<<endl;
717  }
718 
719  for(int i = 0; i < ntowers; i++){
720  char name[100];
721  sprintf(name,"electron_fit_%i",i+1);
722  /*
723  fit[i] = new TF1(name,"pol1(0)+gaus(2)");
724  fit[i]->SetParLimits(0,0,10.0*electron_histo[i]->GetBinContent(1));
725  fit[i]->SetParLimits(1,-10000,0);
726  fit[i]->SetParLimits(2,0,10.0*electron_histo[i]->GetMaximum());
727  fit[i]->SetParLimits(3,0,10);
728  fit[i]->SetParameter(0,electron_histo[i]->GetBinContent(1));
729  fit[i]->SetParameter(1,-electron_histo[i]->GetBinContent(1)/6.0);
730  fit[i]->SetParameter(2,electron_histo[i]->GetMaximum());
731  fit[i]->SetParameter(3,0.95);
732  fit[i]->SetParameter(4,0.15);
733  fit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
734  */
735  fit[i] = new TF1(name,"pol0(0)+gaus(1)");
736  fit[i]->SetParLimits(0,0,10.0*electron_histo[i]->GetBinContent(1));
737  fit[i]->SetParLimits(1,0,10.0*electron_histo[i]->GetMaximum());
738  fit[i]->SetParLimits(2,0,10);
739  fit[i]->SetParameter(0,electron_histo[i]->GetBinContent(1));
740  fit[i]->SetParameter(1,electron_histo[i]->GetMaximum());
741  fit[i]->SetParameter(2,0.95);
742  fit[i]->SetParameter(3,0.15);
743  fit[i]->SetParNames("constant1","constant2","Mean","Sigma");
744 
745  electron_histo[i]->Fit(fit[i],"rql","",0.3,1.7);
746  }
747 
748  ofstream fitfile(ffname);
749 
750  TF1* etacrate_fit[ncrates*20];
751  for(int ii = 0; ii < ncrates; ii++){
752  for(int j = 0; j < 20; j++){
753  TString ecname;
754  ecname += "fit";
755  int cr = ii+1;
756  int et = j;
757  int i = ii*20 + j;
758  //ecname += cr;
759  //ecname += "_";
760  //ecname += et;
761  etacrate_fit[i] = new TF1(ecname.Data(),"pol1(0) + gaus(2)",0.25,1.6);
762  //etacrate_fit[i]->SetParLimits(0,0,10.0*etacrate_histo[i]->GetBinContent(1));
763  etacrate_histo[i]->Rebin();
764  etacrate_fit[i]->SetParLimits(1,-10000,0);
765  etacrate_fit[i]->SetParLimits(2,0,10.0*etacrate_histo[i]->GetMaximum());
766  etacrate_fit[i]->SetParLimits(3,0,10);
767  etacrate_fit[i]->SetParameter(0,etacrate_histo[i]->GetBinContent(2));
768  etacrate_fit[i]->SetParameter(1,-etacrate_histo[i]->GetBinContent(2)/3.0);
769  etacrate_fit[i]->SetParameter(2,etacrate_histo[i]->GetMaximum());
770  etacrate_fit[i]->SetParameter(3,0.96134);
771  etacrate_fit[i]->SetParameter(4,0.141123);
772  etacrate_fit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
773 
774  etacrate_histo[i]->Fit(etacrate_fit[i],"rql","",0.25,1.5);
775  etacrate_histo_p[i]->Fit(etacrate_fit[i],"rql","",0.25,1.5);
776  etacrate_histo_n[i]->Fit(etacrate_fit[i],"rql","",0.25,1.5);
777  fitfile << i << " " << etacrate_histo[i]->GetFunction("fit")->GetParameter(3) << " " << etacrate_histo_p[i]->GetFunction("fit")->GetParameter(3) << " " << etacrate_histo_n[i]->GetFunction("fit")->GetParameter(3) << endl;
778 
779  }
780  }
781 
782  ofstream newgain(ngname);
783 
784 
785  float gains2[ntowers];
786  float gerr2[ntowers];
787  for(int i = 0; i < ntowers; i++){
788  float eta = helper->getEta(i+1);
789  int crate,sequence;
790  decoder->GetCrateFromTowerId(i+1,crate,sequence);
791  int geantetaindex = ((TMath::Nint(fabs(eta) * 1000.0) + 25)/50 - 1);
792  TString ecname;
793  ecname += "fit";
794  //int cr = ii+1;
795  //int et = j;
796  //int i = ii*20 + j;
797  //ecname += crate;
798  //ecname += "_";
799  //ecname += geantetaindex;
800  if(TMath::Abs(eta) > 0.968) eta += 0.005 * TMath::Abs(eta)/eta;
801  int etaindex = ((TMath::Nint(eta * 1000.0) + 25)/50 + 19);
802  //float adjust = ring2fit[(int)etaindex/2]->GetParameter(3);
803  //cout<<etaindex<<" "<<(int)etaindex/2<<" "<<adjust<<endl;
804  //float adjust = fit[i]->GetParameter(3);
805  float ng = 0;
806  float ne = 0;
807  if(status[i] == 1){
808  //float og = bemctables->calib(1,i+1)*gains[i];
809  float og = gains[i];
810  //float aerr = ring2fit[(int)etaindex/2]->GetParError(3);
811  float adjust = etacrate_histo[(crate-1)*20+geantetaindex]->GetFunction("fit")->GetParameter(3);
812  float aerr = etacrate_histo[(crate-1)*20+geantetaindex]->GetFunction("fit")->GetParError(3);
813  if(geantetaindex == 19){
814  adjust = ringfit[etaindex]->GetParameter(3);
815  aerr = ringfit[etaindex]->GetParError(3);
816  }
817  ng = og/adjust;
818  float gerr = peakerr[i]*gains[i]/peaks[i];
819  ne = sqrt(pow(og*aerr/(adjust*adjust),2) + pow(gerr/adjust,2));
820  }
821  newgain << i+1 << " " << ng << " " << ne << " " << status[i] << endl;
822  gains2[i] = ng;
823  gerr2[i] = ne;
824 
825 
826  }
827 
828  newgain.close();
829 
831  //Using new gains regenerate by crate//
833 
834 
835  /*
836  for(int j=0; j<nentries; j++){
837  tree->GetEntry(j);
838  track = &(cluster->centralTrack);
839  TClonesArray *tracks = cluster->tracks;
840 
841  int bsmdeadctot = 0;
842  int bsmdpadctot = 0;
843  for(int i = 0; i < 11; i++){
844  if(track->smde_adc[i] > track->smde_pedestal[i])bsmdeadctot += track->smde_adc[i] - track->smde_pedestal[i];
845  if(track->smdp_adc[i] > track->smdp_pedestal[i])bsmdpadctot += track->smdp_adc[i] - track->smdp_pedestal[i];
846  }
847 
848  double dR = TMath::Sqrt(track->deta*track->deta + track->dphi*track->dphi);
849 
850 
851  double scaled_adc = (track->tower_adc[0] - track->tower_pedestal[0]) / track->p;
852  double geant_scale = geant_fit->Eval(dR);
853  scaled_adc *= geant_scale;
854  //cout<<scaled_adc<<endl;
855  int index = track->tower_id[0];
856 
857  //figure out eta and etaindex
858  eta = helper->getEta(index);
859  if(TMath::Abs(eta) > 0.968) eta += 0.005 * TMath::Abs(eta)/eta;
860  etaindex = ((TMath::Nint(eta * 1000.0) + 25)/50 + 19);
861 
862  double tgain = bemctables->calib(1,track->tower_id[0])*gains[index-1]*gains2[index-1];
863  //double tgain = gains[index-1];
864 
865 
866  if((track->tower_adc[0] - track->tower_pedestal[0]) < 2.5 * track->tower_pedestal_rms[0])continue;
867 
868  if(track->tower_id[0] != track->tower_id_exit)continue;
869 
870  if(track->htTrig == 2)continue;
871  if(track->p > 6)continue;
872 
873  if(dR > 0.025)continue;
874 
875  if(track->dEdx < 3.4e-6)continue;
876  if((bsmdeadctot > -1 && bsmdeadctot < 50) && (bsmdpadctot < 50 && bsmdpadctot > -1))continue;
877 
878  int numsis = tracks->GetEntries();
879  if(numsis > 0)continue;
880 
881  float totalbtow = 0;
882  float maxEt = 0;
883  int maxId = -1;
884  for(int i = 0; i < 9; i++){
885  if(track->tower_adc[i] - track->tower_pedestal[i] < 0)continue;
886  float theta = helper->getTheta(track->tower_id[i]);
887  float nextEt = (track->tower_adc[i] - track->tower_pedestal[i]) * bemctables->calib(1,track->tower_id[i])*sin(theta);
888  totalbtow += nextEt;
889  if(nextEt > maxEt){
890  maxEt = nextEt;
891  maxId = i;
892  }
893  }
894  if(maxId != 0) continue;
895 
896  eta = helper->getEta(index);
897  float phi = helper->getPhi(index);
898  int crate = lookup_crate(eta,phi);
899  if(status[index-1]==1)crate_histo[crate-1]->Fill(scaled_adc*tgain);
900  }
901 
902  for(int i = 0; i < ncrates; i++){
903  sprintf(name,"crate_fit_%i",i);
904  crate_histo[i]->Sumw2();
905  crate_histo[i]->Rebin(4);
906 
907  crate_fit[i] = new TF1(name,"pol1(0) + gaus(2)",0.3,1.7);
908  crate_fit[i]->SetParLimits(0,0,10.0*crate_histo[i]->GetBinContent(1));
909  crate_fit[i]->SetParLimits(1,-10000,0);
910  crate_fit[i]->SetParLimits(2,0,10.0*crate_histo[i]->GetMaximum());
911  crate_fit[i]->SetParLimits(3,0,10);
912  crate_fit[i]->SetParameter(0,crate_histo[i]->GetBinContent(1));
913  crate_fit[i]->SetParameter(1,-crate_histo[i]->GetBinContent(1)/6.0);
914  crate_fit[i]->SetParameter(2,-crate_histo[i]->GetMaximum());
915  crate_fit[i]->SetParameter(3,0.929);
916  crate_fit[i]->SetParameter(4,0.156);
917  crate_fit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
918  crate_fit[i]->SetLineColor(kBlue);
919  crate_fit[i]->SetLineWidth(0.6);
920 
921 
922  crate_histo[i]->Fit(crate_fit[i],"rql","",0.3,1.8);
923  float mean = crate_histo[i]->GetFunction(name)->GetParameter(3);
924  float merr = crate_histo[i]->GetFunction(name)->GetParError(3);
925  crateprec->SetBinContent(i+1,mean);
926  crateprec->SetBinError(i+1,merr);
927  cout<<"crate "<<i+1<<" "<<mean<<" "<<merr/mean<<endl;
928  }
929 
930  ofstream newgain(ngname);
931 
932  float gains3[ntowers];
933  float gerr3[ntowers];
934  for(int i = 0; i < ntowers; i++){
935  float eta = helper->getEta(i+1);
936  float phi = helper->getPhi(i+1);
937  int crate = lookup_crate(eta,phi);
938  float adjust = crate_fit[crate-1]->GetParameter(3);
939  float og = bemctables->calib(1,i)*gains[i]*gains2[i];
940  float ng = og;
941  float aerr = crate_fit[crate-1]->GetParError(3);
942  float ne = sqrt(pow(gains[i]*gerr2[i],2) + pow(gains2[i]*gainerr[i],2));
943  if(fabs(adjust-1)/aerr > 1.5){
944  ne = sqrt(pow(ne/(adjust),2)+pow(og*aerr/(adjust*adjust),2));
945  ng /= adjust;
946  }
947  newgain << i+1 << " " << ng << " " << ne << " " << status[i] << endl;
948  gains3[i] = ng;
949  gerr3[i] = ne;
950  }
951 
952  newgain.close();
953  */
954  outfile.Write();
955  outfile.Close();
956 
957 }
958 
959 //constrained double-Gaussian to fit the background
960 double background_only_fit(double *x, double *par){
961  double par3 = par[0]/1.5;
962  double par4 = par[1] + 10.;
963  double par5 = par[2] * 6.5;
964 
965  double fitval = 0;
966  if(par[2] != 0){
967  double arg1 = (x[0]-par[1])/par[2];
968  double arg2 = (x[0]-par4)/par5;
969  fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2);
970  }
971 
972  return fitval;
973 }
974 
975 double fit_function(double *x, double *par){
976  //6-parameter fit includes
977  //3 param electron Gaussian
978  //par3 = relative height of peak/bg ~ 10
979  //par4 = mean of main bg Gaussian ~ 10
980  //par5 = width of main bg Gaussian ~ 3
981 
982  double par3 = 0;
983  if(par[3] != 0)par3 = par[0] / par[3];
984  double par6 = par3/1.5;
985  double par7 = par[4] + 10.;
986  double par8 = par[5] * 6.5;
987 
988  double fitval = 0;
989  if(par[2] != 0 && par[5] != 0){
990  double arg1 = (x[0]-par[1])/par[2];
991  double arg2 = (x[0]-par[4])/par[5];
992  double arg3 = (x[0]-par7)/par8;
993  fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2) + par6*TMath::Exp(-0.5*arg3*arg3);
994  }
995 
996  return fitval;
997 }
998 
999 double fit_function2(double *x, double *par){
1000  //6-parameter fit includes
1001  //3 param electron Gaussian
1002  //par3 = relative height of peak/bg ~ 10
1003  //par4 = mean of main bg Gaussian ~ 10
1004  //par5 = width of main bg Gaussian ~ 3
1005 
1006  double par3 = 0;
1007  if(par[3] != 0)par3 = par[0] / par[3];
1008  double par6 = par[6];
1009  double par7 = par[4] + 10.;
1010  double par8 = par[7];
1011 
1012  double fitval = 0;
1013  if(par[2] != 0 && par[5] != 0 && par8 != 0){
1014  double arg1 = (x[0]-par[1])/par[2];
1015  double arg2 = (x[0]-par[4])/par[5];
1016  double arg3 = (x[0])/par8;
1017  fitval = par[0]*TMath::Exp(-0.5*arg1*arg1) + par3*TMath::Exp(-0.5*arg2*arg2) + par6*TMath::Exp(-arg3);
1018  }
1019  return fitval;
1020 }
1021 
1022 int lookup_crate(float eta,float phi)
1023 {
1024  if (eta < 0){
1025  if (phi < -2.72) return 5;
1026  else if (phi < -2.30) return 6;
1027  else if (phi < -1.88) return 7;
1028  else if (phi < -1.46) return 8;
1029  else if( phi < -1.04) return 9;
1030  else if( phi < -0.62) return 10;
1031  else if( phi < -0.20) return 11;
1032  else if( phi < 0.22) return 12;
1033  else if( phi < 0.64) return 13;
1034  else if( phi < 1.06) return 14;
1035  else if( phi < 1.48) return 15;
1036  else if( phi < 1.90) return 1;
1037  else if( phi < 2.32) return 2;
1038  else if( phi < 2.74) return 3;
1039  else return 4;
1040  }else{
1041  if (phi < -2.72) return 20;
1042  else if( phi < -2.30) return 21;
1043  else if( phi < -1.88) return 22;
1044  else if( phi < -1.46) return 23;
1045  else if( phi < -1.04) return 24;
1046  else if( phi < -0.62) return 25;
1047  else if( phi < -0.20) return 26;
1048  else if( phi < 0.22) return 27;
1049  else if( phi < 0.64) return 28;
1050  else if( phi < 1.06) return 29;
1051  else if( phi < 1.48) return 30;
1052  else if( phi < 1.90) return 16;
1053  else if( phi < 2.32) return 17;
1054  else if( phi < 2.74) return 18;
1055  else return 19;
1056  }
1057 }
StBemcTablesWriter::loadTables
void loadTables(const char *sqlTime, const char *flavor="ofl")
load directly from DB, no Maker needed
Definition: StBemcTablesWriter.cxx:61
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
StBemcTables::getDecoder
StEmcDecoder * getDecoder()
Return pointer to decoder.
Definition: StBemcTables.h:137
Generated by   doxygen 1.8.5