Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CDC_gains.C
Go to the documentation of this file.
1 
2 
3 void CDC_gains(int EXIT_EARLY=0) {
4 
5  // set EXIT_EARLY to 1 to fit the sum histograms only, 0 for the individual straw gains
6 
7 
8  // Fits histograms from CDC_amp plugin to estimate CDC gains
9  // writes digi_scales/ascale to cdc_new_ascale.txt
10  // writes wire_gains to cdc_new_wiregains.txt
11  // writes fitted histograms to cdc_amphistos.root
12  //
13  // wire_gains are gain correction constants for individual straws, with respect to sum histogram
14  // fit landau to histogram for sum of all straws and histograms for each individual straw
15  // set wire_gains to scale individual fit mpvs to the mpv for the sum of all straws histogram.
16 
17  // digi_scales/ascale is the overall gain correction factor for the chamber
18  // Reference values of IDEALMPV and ASCALE are chosen from a high gain (low pressure) run 11621
19  // The new ascale is ASCALE*IDEALMPV/sum fit mpv
20  //
21  // (For the first round of gain calibrations, ASCALE was chosen to match the simulation charge histogram.
22  // For the next round, ASCALE was chosen to match the dedx histograms from the first round).
23 
24  // CDC_amp provides 2D histograms of pulse amplitude vs straw
25  // The histograms are projected here into 1D histograms, one per straw
26  // For Landau-like fits, need tracked hits with a narrow range of theta
27 
28  // NSJ 18 Nov 2016
29 
30  // Need 80+ run files to get good fits for individual straws. 2 run files are enough for sum.
31 
32 
33  // NSJ 24 Apr 2018
34  // Updated in 2018 to use hits within the first 100ns only. This avoids the gain loss at later drift times which is probably caused by oxygen in the chamber.
35 
36 
37 
38  // reference values
39 
40  // const float IDEALMPV=32.385; //mpv for tracked hits with restricted z,theta from low pressure run 11621
41 
42  // const float IDEALMPV=37.9649; //mpv for tracked hits ***at 0-100ns*** with restricted z,theta from low pressure run 11621
43 
44  // NSJ 24 Aug 2018 modified to use 30 degree tracks instead of 90 degree tracks.
45  // Still use the restricted time range 0-100ns
46  // Both of these give a higher amplitude peak which is easier to fit for high pressure runs
47 
48  const float IDEALMPV=88.5315*2.02/2.59; //mpv for hits in 0-100ns with theta 28-32 degrees, z 52-78cm, from low pressure run 11621, adjusted to put dE/dx (1.5GeV) at 2.02 keV/cm to match geant
49 
50  const float ASCALE=0.176; //ascale for tracked ztheta hits from 011621
51 
52  const bool REBIN_HISTOS=kFALSE;
53 
54 
55  // which theta histos to use for fits
56 
57  int theta = 30;
58  int theta_low = 28;
59  int theta_high = 32;
60 
61 
62  // get histograms
63 
64  TDirectory *fmain = (TDirectory*)gDirectory->FindObjectAny("CDC_amp");
65  if (!fmain) printf("Cannot find directory CDC_amp\n");
66  if (!fmain) return;
67  fmain->cd();
68 
69 
70  TH1I *attsum = (TH1I*)fmain->Get(Form("a%i_100ns",theta));
71  if (!attsum) printf("Cannot find histogram a%i_100ns\n",theta);
72  if (!attsum) return;
73 
74  // get sum amplitude histo to find readout range
75 
76  TH1I *asum_all= (TH1I*)fmain->Get("a");
77  if (!asum_all) printf("Cannot find amplitude (sum over all straws) histogram a\n");
78  if (!asum_all) return;
79 
80  // find amplitude range
81  int SCALE_UP = 1;
82  double highcounts = asum_all->FindLastBinAbove(0);
83  if (highcounts > 512) SCALE_UP = 8; // this is full range 0-4095
84 
85  if (SCALE_UP==1) cout << "Amp range 0 to 511" << endl;
86  if (SCALE_UP==8) cout << "Amp range 0 to 4095" << endl;
87 
88  // fit range constants
89 
90  int AL;
91  int ANL;
92  int AH;
93 
94  int INCREMENT_FITSTART;
95 
96 
97  if (SCALE_UP==1) {
98 
99  AL = 20; // lower limit attsum
100 
101  ANL=14; // lower limit attsum indiv straw fits
102 
103  AH=400; //amp fit upper limit landau
104 
105  INCREMENT_FITSTART=1;
106 
107  } else {
108 
109  AL = 130; // lower limit attsum
110 
111  ANL=110; // lower limit attsum indiv straw fits
112 
113  AH=3000; //amp fit upper limit landau
114 
115  INCREMENT_FITSTART=8;
116 
117  }
118 
119  const int MINCOUNTS=1000; //counts required to fit histogram
120 
121 
122  gStyle->SetOptFit(1);
123  gStyle->SetFuncWidth(1);
124 
125  TF1 *f = new TF1("f","landau");
126  f->SetLineColor(6);
127 
128  int a_fitstat,q_fitstat;
129 
130  float newascale = 0;
131 
132  float thismpv = 0;
133 
134 
135  new TCanvas;
136  printf("\nattsum fit, tracked hits at 0-100ns, restricted z and theta:\n");
137  f->SetRange(AL,AH);
138  a_fitstat = attsum->Fit(f,"RW");
139 
140  if (!a_fitstat) {
141  thismpv = f->GetParameter(1);
142  printf("\n This mpv: %.3f ideal mpv: %.3f\n",thismpv,IDEALMPV*SCALE_UP);
143  newascale = ASCALE*IDEALMPV*SCALE_UP/thismpv;
144  printf("\nnew digi_scales/ascale should be %.3f\n\n",newascale);
145  }
146 
147  if (a_fitstat) printf("\nSum histogram fit error, exiting script\n");
148  if (a_fitstat) return;
149 
150 
151  FILE *outfile = fopen("cdc_new_ascale.txt","w");
152  fprintf(outfile,"%.3f 0.8\n",newascale); //0.8 is the time scale, for in the 2nd column of ccdb's CDC/digi_scales table
153  fclose(outfile);
154 
155 
156  TFile *hfile = new TFile("cdc_amphistos.root","RECREATE");
157 
158  hfile->cd();
159  attsum->Write();
160 
161  if (EXIT_EARLY==1) hfile->Write();
162  if (EXIT_EARLY==1) hfile->Close();
163  if (EXIT_EARLY==1) return;
164 
165 
166 
167  TTree *attstats = new TTree("attstats","fit stats for tracked hits with restricted z and theta");
168 
169  int n,ring; //straw number (1-3522), ring number (1-28)
170 
171  int a_n; //count
172  double a_mean;
173  double a_c; //const
174  double a_mpv;
175  double a_mpverr;
176  double a_sig;
177  double a_chisq;
178 
179  int fitlowerlimit;
180 
181  attstats->Branch("n",&n,"n/I");
182  attstats->Branch("hits",&a_n,"hits/I");
183  attstats->Branch("mean",&a_mean,"mean/D");
184  attstats->Branch("c",&a_c,"c/D");
185  attstats->Branch("mpv",&a_mpv,"mpv/D");
186  attstats->Branch("mpverr",&a_mpverr,"mpverr/D");
187  attstats->Branch("sig",&a_sig,"sig/D");
188  attstats->Branch("chisq",&a_chisq,"chisq/D");
189  attstats->Branch("fitstat",&a_fitstat,"fitstat/I");
190  attstats->Branch("fitlowerlimit",&fitlowerlimit,"fitlowerlimit/I");
191 
192  int i;
193 
194  TH1D *ahisto;
195  TH2I *anhisto;
196  char htitle[300];
197 
198 
199  int badfit, nofit; //counters
200  int bincont, lastbin, nbins, lastbintocheck; // to flag low-gain preamp channels
201 
202  float wiregain;
203 
204 
205 
206  hfile->cd();
207 
208  TDirectory *hmain = gDirectory;
209 
210  outfile = fopen("cdc_new_wiregains.txt","w");
211 
212  new TCanvas;
213 
214  hmain->cd();
215 
216  gDirectory->mkdir("theta")->cd();
217 
218  TDirectory *hsub = gDirectory;
219 
220  fmain->cd();
221 
222  anhisto = (TH2I*)fmain->Get(Form("an%i_100ns",theta));
223 
224 
225  badfit = 0;
226  nofit = 0;
227 
228 
229  for (i=1; i<3523; i++) {
230  // for (i=1; i<10; i++) {
231 
232  fitlowerlimit = ANL;
233 
234  f->SetRange(fitlowerlimit,AH);
235 
236  a_c = 0;
237  a_mpv = 0;
238  a_mpverr = 0;
239  a_sig = 0;
240  a_chisq = 0;
241  a_fitstat = 0;
242 
243  n = i;
244 
245  ahisto = anhisto->ProjectionY(Form("a[%i]",i),i+1,i+1); //straw 1 is in bin 2
246 
247  if (REBIN_HISTOS) ahisto->Rebin(4);
248 
249 
250  sprintf(htitle,"Amplitude, tracked hits, z=52 to 78 cm, theta=%i to %i degrees, straw %i; amplitude-pedestal",theta_low,theta_high,i);
251 
252  ahisto->SetTitle(htitle);
253 
254  a_n = ahisto->GetEntries();
255  a_mean = ahisto->GetMean();
256 
257  // cout << "fitting straw " << i << endl;
258 
259  a_fitstat = -1;
260  if (a_n > MINCOUNTS) {
261  a_fitstat = ahisto->Fit(f,"QW","",fitlowerlimit,AH); //landau LL is for low stats
262  }
263 
264 
265  // if fit did not converge, increase lower limit
266 
267  while (a_fitstat==4 && fitlowerlimit<200) {
268  fitlowerlimit += INCREMENT_FITSTART;
269  a_fitstat = ahisto->Fit(f,"W","",fitlowerlimit,AH);
270  }
271 
272  if (a_fitstat==0) {
273 
274  a_c = f->GetParameter(0);
275  a_mpv = f->GetParameter(1);
276  a_mpverr = f->GetParError(1);
277  a_sig = f->GetParameter(2);
278 
279  if (f->GetNDF()>0) a_chisq = f->GetChisquare()/f->GetNDF();
280 
281  } else if (a_fitstat==4){
282  printf("straw %i unconverged fit \n",i);
283 
284  } else if (a_fitstat>0){
285  printf("straw %i fitstatus %i fit MPV %f\n",i,a_fitstat,f->GetParameter(1));
286 
287  } else if (a_fitstat<0){
288  printf("straw %i not enough counts \n",i);
289  nofit++;
290 
291  }
292 
293 
294 
295 
296  if (a_fitstat>0) badfit++;
297 
298 
299  if (a_mpv<0) printf("straw %i fit MPV %f\n",i,a_mpv);
300  if (a_mpv>0 && a_mpv<fitlowerlimit) printf("straw %i MPV %f below fit lower limit\n",i,a_mpv);
301  if (a_mpv>0 && a_mean>3.3*a_mpv) printf("straw %i mean %f MPV %f\n",i,a_mean,a_mpv);
302 
303 
304  attstats->Fill();
305 
306  wiregain=0;
307  if (a_mpv>0) wiregain = thismpv/a_mpv;
308  fprintf(outfile,"%.3f\n",wiregain);
309 
310 
311  hsub->cd();
312 
313  ahisto->Write();
314 
315  }
316 
317  printf("\nfitstatus 4 count: %i \n",badfit);
318  printf("no fit count: %i \n",nofit);
319 
320 
321 
322  fclose(outfile);
323 
324  hfile->cd();
325 
326  attstats->Write();
327 
328 
329 }
sprintf(text,"Post KinFit Cut")
void CDC_gains(int EXIT_EARLY=0)
Definition: CDC_gains.C:3
TF1 * f
Definition: FitGains.C:21
TFile * outfile
Definition: tw_corr.C:46
printf("string=%s", string)