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