ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/yangyong/combineICEB/combineCalibConstantv2.C
Revision: 1.2
Committed: Tue Apr 10 19:41:08 2012 UTC (13 years, 1 month ago) by yangyong
Content type: text/plain
Branch: MAIN
CVS Tags: V2012combv1, V2011combv1, pi0etaeb_laser20111122
Changes since 1.1: +1667 -0 lines
Log Message: 
LaserTag:
EcalLaserAPDPNRatios_data_20111122_158851_180363

and alpha Tag ( for Endcap Only)
EcalLaserAlphas_lto420-620_progr_data_20111122

File Contents

# User Rev Content
1 yangyong 1.2 #include "rootheader.h"
2     #include "testSelection_cluster.h"
3    
4     #include "common_functions.cc"
5     #include "foldEndcap.cc"
6    
7     TRandom3 *rgen_;
8    
9     #include "usefullcode.cc"
10    
11     float interCalib_preCalib[170][360];
12     float interCalib_preCalibTrue[170][360];
13     float interCalibEndcap_preCalib[2][101][101];
14    
15     void copyConstant(float c1[170][360], float c2[170][360]){
16    
17     for(int j=0; j<170; j++){
18     for(int k=0; k<360; k++){
19     c2[j][k] = c1[j][k];
20     }
21     }
22    
23     }
24    
25     void isAtEcalBarrelModuleCracks(int ieta, int iphi, bool &etaCracks, bool &phiCracks){
26     int absieta = abs(ieta);
27    
28     etaCracks = false;
29     phiCracks = false;
30    
31     if( absieta == 25 || absieta == 26 || absieta == 45 || absieta == 46 || absieta == 65 || absieta == 66) etaCracks = true;
32     if( iphi %20 == 1 || iphi %20 == 0) phiCracks = true;
33    
34    
35     }
36    
37     float ic[50][170][360]; ///at most 20 set
38    
39     float icv1[50][170][360];
40    
41     int ndeadflagietaiphi_ic[50][170][360];
42    
43     float icwt_period[50][170][360]; //pi0&eta combined for each period
44    
45    
46    
47     bool isTestBeamSM(int iSM){
48    
49    
50    
51     //if( iSM == 1 || iSM == 2 || iSM == 11 || iSM == 15 || iSM == 21 || iSM == 24) return true;
52    
53     if( iSM == 1 || iSM == 11 || iSM == 15 || iSM == 24) return true;
54    
55    
56     else return false;
57    
58     }
59    
60    
61     #include "effSigma.C"
62    
63     #include "getCrystaldeadflag.cc"
64    
65     #include "getCalibConstants.cc"
66    
67     #include "gausfit.cc"
68    
69     float corrGR10_P_V10_over_GR09_R_V6A[170][360];
70     float corrGR_R_311_V1A_over_GR09_R_V6A[170][360];
71     float corrGR_R_311_V1A_over_GR10_P_V10[170][360];
72    
73    
74     float CphiCorrin[170][360];
75     float CBSCorrin[170][360];
76    
77     float fitWind = 3;
78     float sigmaTB = 0.55;
79     //float sigmaTB = 0.44;
80     //float sigmaTB = 0.95;
81    
82     void fitHistogram(TH1F *h1, double res[]){
83    
84     float mean = h1->GetMean();
85     float meanErr = h1->GetMeanError();
86     float rmsEff = 100 * h1->GetRMS();
87     float rmsEffErr = 100 * h1->GetRMSError();
88     float rmsGaus = 100 * h1->GetRMS();
89     float rmsGausErr = 100 * h1->GetRMSError();
90     if( h1->Integral()>50){
91     double resd[10];
92     effSigma(h1,resd);
93     rmsEff = resd[2]*100;
94     float resf[20];
95    
96     fitgauswind2(h1,fitWind,fitWind,resf);
97     //fitgauswindRefit(h1,resf);
98     rmsGaus = resf[6]*100;
99     rmsGausErr = resf[3] * 100;
100     rmsEffErr = resf[3] * 100;
101     mean = resf[0];
102     meanErr = resf[1];
103     }
104     res[0] = mean;
105     res[1] = meanErr;
106     res[2] = rmsEff;
107     res[3] = rmsEffErr;
108     res[4] = rmsGaus;
109     res[5] = rmsGausErr;
110     }
111    
112    
113    
114     float cpizv1[170][360];
115     float cpizv2[170][360];
116     float cpizv1in[170][360];
117     float cpizv2in[170][360];
118    
119     float cetav1[170][360];
120     float cetav2[170][360];
121     float cetav1in[170][360];
122     float cetav2in[170][360];
123    
124    
125     float cpizoutput[170][360];
126     float cetaoutput[170][360];
127     float cpizetaoutput[170][360];
128    
129     float cpizetaphibsoutput[170][360];
130     float cpizetaphibsoutputv1[170][360];
131    
132     float cpizetaphibsprecaliboutput[170][360];
133     float cpizetaphibsprecaliboutputFinal[170][360];
134     float cpizetaphibsprecaliboutputv1[170][360];
135    
136    
137     float cphibsoutput[170][360];
138     float cphibsoutputv1[170][360];
139    
140    
141     float cpizoutputv1[170][360];
142     float cetaoutputv1[170][360];
143     float cpizetaoutputv1[170][360];
144    
145    
146     void rescaleConstTo(float Cin[170][360], float Corr[170][360]){
147    
148     float mean = 0;
149     int nmean = 0;
150     for(int j=0; j< 170; j++){
151     for(int k=0; k< 360; k++){
152     if( Cin[j][k]> 0){
153     Cin[j][k] /= Corr[j][k];
154     nmean ++;
155     mean += Cin[j][k];
156     }
157     }
158     }
159     mean /= nmean;
160     for(int j=0; j< 170; j++){
161     for(int k=0; k< 360; k++){
162     if( Cin[j][k]> 0){
163     Cin[j][k] /= mean;
164     }
165     }
166     }
167     }
168    
169     ///Cprecaib * c1 == c_new * c2
170     ///c2 = c1 / ( c_new / cprecalib)
171    
172     void rescaleConstTov1(float Cin[170][360], float Corr[170][360]){
173    
174     float mean = 0;
175     int nmean = 0;
176     for(int j=0; j< 170; j++){
177     for(int k=0; k< 360; k++){
178     if( Cin[j][k]> 0){
179     Cin[j][k] *= Corr[j][k];
180     nmean ++;
181     mean += Cin[j][k];
182     }
183     }
184     }
185     mean /= nmean;
186     for(int j=0; j< 170; j++){
187     for(int k=0; k< 360; k++){
188     if( Cin[j][k]> 0){
189     Cin[j][k] /= mean;
190     }
191     }
192     }
193     }
194    
195    
196     ///code to combined different pi0/eta EB IC for each laserTag
197    
198    
199     void combineCalibConstantv2(){
200    
201    
202    
203     readInterCalibConstEBSimple("interCalibEB_GR_R_42_V13.txt",interCalib_preCalib); ///2010IC
204    
205    
206     ///pre-calibration constants
207     readInterCalibConstEBSimple("interCalib_GR09_H_V6OFF_barrel.txt",interCalib_preCalibTrue);
208    
209    
210     map<int, string> smScaleFiles;
211    
212     //7 Sets IC
213     /// 160404to166923
214     /// 166946to167913
215     /// 169985to172630
216     /// 172635to172791
217     /// 175832to176929
218     /// 176933to177878
219     /// 178003to180252
220    
221    
222     //SM-scale files Pi0
223     smScaleFiles[0] = "calibres/deriveCalibConst.testCalibv1.dflag60.pe1.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
224     smScaleFiles[1] = "calibres/deriveCalibConst.testCalibv1.dflag61.pe1.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
225     smScaleFiles[2] = "calibres/deriveCalibConst.testCalibv1.dflag62.pe1.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
226     smScaleFiles[3] = "calibres/deriveCalibConst.testCalibv1.dflag63.pe1.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
227     smScaleFiles[4] = "calibres/deriveCalibConst.testCalibv1.dflag72.pe1.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
228     smScaleFiles[5] = "calibres/deriveCalibConst.testCalibv1.dflag73.pe1.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
229     smScaleFiles[6] = "calibres/deriveCalibConst.testCalibv1.dflag64.pe1.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
230    
231     //SM-scale files Eta
232     smScaleFiles[7] = "calibres/deriveCalibConst.testCalibv1.dflag250.pe2.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
233     smScaleFiles[8] = "calibres/deriveCalibConst.testCalibv1.dflag251.pe2.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
234     smScaleFiles[9] = "calibres/deriveCalibConst.testCalibv1.dflag252.pe2.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
235     smScaleFiles[10] = "calibres/deriveCalibConst.testCalibv1.dflag253.pe2.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
236     smScaleFiles[11] = "calibres/deriveCalibConst.testCalibv1.dflag227.pe2.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
237     smScaleFiles[12] = "calibres/deriveCalibConst.testCalibv1.dflag228.pe2.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
238     smScaleFiles[13] = "calibres/deriveCalibConst.testCalibv1.dflag254.pe2.cut0.rmOvlap0.step1.method1.corrEta-1.corrPhi-1.corrSM0.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root";
239    
240    
241     // number of ICs
242     int nIC = 7;
243    
244    
245    
246     map<int,string> icFiles;
247    
248     //IC Pi0
249     icFiles[0] = "calibres/deriveCalibConst.testCalibv1.dflag60.pe1.cut0.rmOvlap0.step25.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
250     icFiles[1] = "calibres/deriveCalibConst.testCalibv1.dflag61.pe1.cut0.rmOvlap0.step25.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
251     icFiles[2] = "calibres/deriveCalibConst.testCalibv1.dflag62.pe1.cut0.rmOvlap0.step25.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
252     icFiles[3] = "calibres/deriveCalibConst.testCalibv1.dflag63.pe1.cut0.rmOvlap0.step25.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
253     icFiles[4] = "calibres/deriveCalibConst.testCalibv1.dflag72.pe1.cut0.rmOvlap0.step25.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
254     icFiles[5] = "calibres/deriveCalibConst.testCalibv1.dflag73.pe1.cut0.rmOvlap0.step25.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
255     icFiles[6] = "calibres/deriveCalibConst.testCalibv1.dflag64.pe1.cut0.rmOvlap0.step25.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
256    
257     //IC Eta
258     icFiles[7] = "calibres/deriveCalibConst.testCalibv1.dflag250.pe2.cut0.rmOvlap0.step45.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
259     icFiles[8] = "calibres/deriveCalibConst.testCalibv1.dflag251.pe2.cut0.rmOvlap0.step45.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
260     icFiles[9] = "calibres/deriveCalibConst.testCalibv1.dflag252.pe2.cut0.rmOvlap0.step45.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
261     icFiles[10] = "calibres/deriveCalibConst.testCalibv1.dflag253.pe2.cut0.rmOvlap0.step45.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
262     icFiles[11] = "calibres/deriveCalibConst.testCalibv1.dflag227.pe2.cut0.rmOvlap0.step45.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
263     icFiles[12] = "calibres/deriveCalibConst.testCalibv1.dflag228.pe2.cut0.rmOvlap0.step45.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
264     icFiles[13] = "calibres/deriveCalibConst.testCalibv1.dflag254.pe2.cut0.rmOvlap0.step45.method2.corrEta21.corrPhi11.corrSM1.corrDead1.precalib1.vtx1.encorr0.evtNot-1.trig0.txt";
265    
266    
267     //crystal dead flag
268     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag60.txt",ndeadflagietaiphi_ic[0]);
269     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag61.txt",ndeadflagietaiphi_ic[1]);
270     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag62.txt",ndeadflagietaiphi_ic[2]);
271     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag63.txt",ndeadflagietaiphi_ic[3]);
272     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag72.txt",ndeadflagietaiphi_ic[4]);
273     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag73.txt",ndeadflagietaiphi_ic[5]);
274     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag64.txt",ndeadflagietaiphi_ic[6]);
275    
276     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag250.txt",ndeadflagietaiphi_ic[7]);
277     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag251.txt",ndeadflagietaiphi_ic[8]);
278     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag252.txt",ndeadflagietaiphi_ic[9]);
279     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag253.txt",ndeadflagietaiphi_ic[10]);
280     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag227.txt",ndeadflagietaiphi_ic[11]);
281     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag228.txt",ndeadflagietaiphi_ic[12]);
282     getCrystaldeadflagBarrel_v1("deadflag/crystal_deadflag_eb_dflag254.txt",ndeadflagietaiphi_ic[13]);
283    
284    
285     ofstream txtoutTocheck("combinedCalibConstantv2.txt");
286    
287    
288     int nConstantSet = int(icFiles.size());
289    
290     cout<<" nConstantSet " << nConstantSet <<endl;
291    
292    
293     map<int,TH1F*> hh_smscales;
294     float cc[170][360];
295    
296    
297     int nMaxICSet = 50;
298     if( nConstantSet > nMaxICSet){
299     cout<<"more than " << nMaxICSet<<" "<< nConstantSet <<endl;
300     return;
301     }
302    
303    
304    
305     float icwt[170][360];
306     float icwtv1[170][360];
307    
308    
309    
310     // float icwt1[170][360];
311     // float icwt2[170][360];
312    
313    
314     /// ieta -21 to -15 iph 280 to 290
315     bool fixTempTT = true; ///this is a fix-me by hand, not elegant, be careful!
316    
317     // ieta -30 to -25, iphi 165 to 171
318     // ieta -81 to -75, iphi 300 to 310
319    
320    
321     for(int n=0; n< nConstantSet; n++){
322     string filename = smScaleFiles[n];
323     TFile *ff = new TFile(filename.c_str(),"read");
324     TH1F *hhtmp = (TH1F*)ff->Get("hh_corr_smv2");
325     hh_smscales[n] = hhtmp;
326    
327     filename = icFiles[n];
328    
329     readInterCalibConstEBSimplev1(filename.c_str(),cc);
330    
331    
332     if( fixTempTT){
333    
334     float meanTemp = 0;
335     int nTemp = 0;
336     float meanTemp2 =0;
337     int nTemp2 = 0;
338    
339    
340     float meanTemp3=0;
341     int nTemp3 = 0;
342     float meanTemp4=0;
343     int nTemp4 = 0;
344    
345    
346     float meanTemp5=0;
347     int nTemp5 = 0;
348     float meanTemp6=0;
349     int nTemp6 = 0;
350    
351    
352     float meanTemp7=0;
353     int nTemp7 = 0;
354     float meanTemp8=0;
355     int nTemp8 = 0;
356    
357    
358     for(int j=0; j< 170; j++){
359     for(int k=0; k< 360; k++){
360     int ieta = j-85;
361     if( ieta >=0) ieta += 1;
362     int iphi = k;
363     if( k==0) iphi = 360;
364    
365    
366     ////<<<<<<<
367     if( (ieta == 35 && iphi >= 311 && iphi <= 315) ||
368     (ieta == 41 && iphi >= 311 && iphi <= 315) ||
369     (ieta >= 36 && ieta <= 40 && iphi == 310 ) ||
370     (ieta >= 36 && ieta <= 40 && iphi == 316 )
371     ){
372     if( cc[j][k] >0){
373     meanTemp7 += cc[j][k];
374     nTemp7 ++;
375     }
376     }
377     if( (ieta== 35 && iphi == 310) ||
378     (ieta== 35 && iphi == 316) ||
379     (ieta== 41 && iphi == 310) ||
380     (ieta== 41 && iphi == 316) ){
381     if( cc[j][k] >0){
382     meanTemp8 += cc[j][k];
383     nTemp8 ++;
384     }
385     }
386     ////>>>>>>>
387    
388    
389     ////<<<<<<<
390     if( (ieta == -81 && iphi >= 301 && iphi <= 310) ||
391     (ieta == -75 && iphi >= 301 && iphi <= 310) ||
392     (ieta >= -80 && ieta <= -76 && iphi == 300 ) ||
393     (ieta >= -80 && ieta <= -76 && iphi == 311 )
394     ){
395     if( cc[j][k] >0){
396     meanTemp5 += cc[j][k];
397     nTemp5 ++;
398     }
399     }
400     if( (ieta== -81 && iphi == 300) ||
401     (ieta== -81 && iphi == 311) ||
402     (ieta== -75 && iphi == 300) ||
403     (ieta== -75 && iphi == 311) ){
404     if( cc[j][k] >0){
405     meanTemp6 += cc[j][k];
406     nTemp6 ++;
407     }
408     }
409     ////>>>>>>>
410    
411    
412     ////<<<<<<<
413     if( (ieta == -31 && iphi >= 166 && iphi <= 170) ||
414     (ieta == -25 && iphi >= 166 && iphi <= 170) ||
415     (ieta >= -30 && ieta <= -26 && iphi == 165 ) ||
416     (ieta >= -30 && ieta <= -26 && iphi == 171 )
417     ){
418     if( cc[j][k] >0){
419     meanTemp3 += cc[j][k];
420     nTemp3 ++;
421     }
422     }
423     if( (ieta== -31 && iphi == 165) ||
424     (ieta== -31 && iphi == 171) ||
425     (ieta== -25 && iphi == 165) ||
426     (ieta== -25 && iphi == 171) ){
427     if( cc[j][k] >0){
428     meanTemp4 += cc[j][k];
429     nTemp4 ++;
430     }
431     }
432     ////>>>>>>>
433    
434    
435     ////<<<<<<<
436     if( (ieta == -21 && iphi >= 281 && iphi <= 290) ||
437     (ieta == -15 && iphi >= 281 && iphi <= 290) ||
438     (ieta >= -20 && ieta <= -16 && iphi == 280 ) ||
439     (ieta >= -20 && ieta <= -16 && iphi == 291 )
440     ){
441     if( cc[j][k] >0){
442     meanTemp += cc[j][k];
443     nTemp ++;
444     }
445    
446     }
447    
448     if( (ieta== -21 && iphi == 291) ||
449     (ieta== -21 && iphi == 280) ||
450     (ieta== -15 && iphi == 280) ||
451     (ieta== -15 && iphi == 291) ){
452    
453     if( cc[j][k] > 0){
454     meanTemp2 += cc[j][k];
455     nTemp2 ++;
456     }
457     }
458     ////>>>>>>>
459    
460    
461     }
462     }
463    
464    
465    
466     if( n==2 || n==3 || n==9 || n==10){
467     meanTemp5 /= nTemp5;
468     meanTemp6 /= nTemp6;
469     cout<<" meanTemp5 6 " << meanTemp5 <<" "<<meanTemp6 <<endl;
470    
471     for(int j=0; j< 170; j++){
472     for(int k=0; k< 360; k++){
473     int ieta = j-85;
474     if( ieta >=0) ieta += 1;
475     int iphi = k;
476     if( k==0) iphi = 360;
477    
478     if( (ieta == -81 && iphi >= 301 && iphi <= 310) ||
479     (ieta == -75 && iphi >= 301 && iphi <= 310) ||
480     (ieta >= -80 && ieta <= -76 && iphi == 300 ) ||
481     (ieta >= -80 && ieta <= -76 && iphi == 311 )
482     ){
483     if( cc[j][k] >0){
484     cc[j][k] /= meanTemp5;
485     }
486     }
487     if( (ieta== -81 && iphi == 300) ||
488     (ieta== -81 && iphi == 311) ||
489     (ieta== -75 && iphi == 300) ||
490     (ieta== -75 && iphi == 311) ){
491     if( cc[j][k] >0){
492     cc[j][k] /= meanTemp6;
493     }
494     }
495     }
496     }
497     }
498    
499    
500     if( n==0 || n==1 || n==2 || n==3 || n ==4 || n==5 || n==6 || n==7 || n==8 || n==9 || n==10|| n== 11 || n== 12 || n== 13){
501     meanTemp3 /= nTemp3;
502     meanTemp4 /= nTemp4;
503     cout<<" meanTemp3 4 " << meanTemp3 <<" "<<meanTemp4 <<endl;
504    
505     for(int j=0; j< 170; j++){
506     for(int k=0; k< 360; k++){
507     int ieta = j-85;
508     if( ieta >=0) ieta += 1;
509     int iphi = k;
510     if( k==0) iphi = 360;
511    
512     if( (ieta == -31 && iphi >= 166 && iphi <= 170) ||
513     (ieta == -25 && iphi >= 166 && iphi <= 170) ||
514     (ieta >= -30 && ieta <= -26 && iphi == 165 ) ||
515     (ieta >= -30 && ieta <= -26 && iphi == 171 )
516     ){
517     if( cc[j][k] >0){
518     cc[j][k] /= meanTemp3;
519     }
520     }
521     if( (ieta== -31 && iphi == 165) ||
522     (ieta== -31 && iphi == 171) ||
523     (ieta== -25 && iphi == 165) ||
524     (ieta== -25 && iphi == 171) ){
525     if( cc[j][k] >0){
526     cc[j][k] /= meanTemp4;
527     }
528     }
529     }
530     }
531     }
532    
533    
534    
535     if( n==4 || n==5 || n== 11 || n==12){ ///September
536     meanTemp /= nTemp;
537     meanTemp2 /=nTemp2;
538     meanTemp7 /= nTemp7;
539     meanTemp8 /= nTemp8;
540     cout<<" meanTemp " << meanTemp <<" "<<meanTemp2 <<" "<< meanTemp7 <<" "<< meanTemp8<<endl;
541    
542    
543     for(int j=0; j< 170; j++){
544     for(int k=0; k< 360; k++){
545     int ieta = j-85;
546     if( ieta >=0) ieta += 1;
547     int iphi = k;
548     if( k==0) iphi = 360;
549    
550     if( (ieta == -21 && iphi >= 281 && iphi <= 290) ||
551     (ieta == -15 && iphi >= 281 && iphi <= 290) ||
552     (ieta >= -20 && ieta <= -16 && iphi == 280 ) ||
553     (ieta >= -20 && ieta <= -16 && iphi == 291 )
554     ){
555     if( cc[j][k] >0){
556     cc[j][k] /= meanTemp;
557     }
558     }
559    
560     if( (ieta== -21 && iphi == 291) ||
561     (ieta== -21 && iphi == 280) ||
562     (ieta== -15 && iphi == 280) ||
563     (ieta== -15 && iphi == 291) ){
564     if( cc[j][k] > 0){
565     cc[j][k] /= meanTemp2;
566     }
567     }
568    
569     if( (ieta == 35 && iphi >= 311 && iphi <= 315) ||
570     (ieta == 41 && iphi >= 311 && iphi <= 315) ||
571     (ieta >= 36 && ieta <= 40 && iphi == 310 ) ||
572     (ieta >= 36 && ieta <= 40 && iphi == 316 )
573     ){
574     if( cc[j][k] >0){
575     cc[j][k] /= meanTemp7;
576     }
577     }
578     if( (ieta== 35 && iphi == 310) ||
579     (ieta== 35 && iphi == 316) ||
580     (ieta== 41 && iphi == 310) ||
581     (ieta== 41 && iphi == 316) ){
582     if( cc[j][k] >0){
583     cc[j][k] /= meanTemp8;
584     }
585     }
586    
587     }
588     }
589    
590     }
591    
592     }
593    
594    
595    
596     for(int j=0; j< 170; j++){
597     for(int k=0; k< 360; k++){
598     ic[n][j][k] = cc[j][k];
599     }
600     }
601    
602    
603     NormIetaAbsToUnitTestBeamSMOnly(ic[n]);
604    
605     if(n==0) cout<<"check ic[0][1] 1 " << ic[n][0][1]<<endl;
606    
607     SetSMScale(ic[n], hh_smscales[n]);
608    
609     if(n==0) cout<<"check ic[0][1] 2 " << ic[n][0][1]<<endl;
610    
611     NormCrystalDeadFlag_v1(ic[n],ndeadflagietaiphi_ic[n]);
612    
613     if(n==0) cout<<"check ic[0][1] 3 " << ic[n][0][1]<<endl;
614    
615    
616     copyConstant(ic[n],icv1[n]);
617    
618     /////for estimating precision
619     NormSMScaleToUnit(icv1[n]);
620     }
621    
622     //return;
623    
624    
625     TFile *fnew = new TFile("combineCalibConstantv2.root","recreate");
626    
627    
628     TH1F *hh_csm[51][36][3]; //for each SM, all/centra/outer
629    
630    
631     //|ieta|<=45 , all, removing eta phi boundaries, at eta moduaries, at phi boduaris
632     TH1F *hh_csmtb_ietaMod12[50][5];
633     for(int n=0; n< nConstantSet; n++){
634     for(int k=0; k<5; k++){
635     TString histname = TString(Form("hh_csmtb_ietaMod12_%d_%d",n,k));
636     hh_csmtb_ietaMod12[n][k] = new TH1F(histname,histname,500,0,2);
637     }
638     }
639    
640    
641     TH1F *hh_c_ieta[50][170];
642     TH1F *hh_c_ietaAbs[50][85];
643    
644    
645    
646    
647     TH1F *hh_csmtb_ietaTT[50][34];
648     TH1F *hh_csmtb_ietaTTAbs[50][17];
649     TH1F *hh_csmco_ietaTTAbs[50][17];
650     TH1F *hh_csmco_ietaTT[50][34];
651    
652     TH1F *hh_csmall_ietaTTAbs[50][17];
653    
654    
655    
656     TH1F *hh_wtavg_csmtb_ietaTT[34];
657     TH1F *hh_wtavg_csmtb_ietaTTAbs[17];
658     TH1F *hh_wtavg_csmco_ietaTT[34];
659     TH1F *hh_wtavg_csmco_ietaTTAbs[17];
660    
661     TH1F *hh_wtavg_csmall_ietaTTAbs[17];
662    
663    
664     TH1F *hh_res_csmtbietaTT[50][3];
665     TH1F *hh_res_csmtbietaTTAbs[50][3];
666    
667     TH1F *hh_res_csmcoietaTT[50][3];
668     TH1F *hh_res_csmcoietaTTAbs[50][3];
669    
670     TH1F *hh_res_csmallietaTTAbs[50][3];
671    
672    
673    
674     TH1F *hh_res_wtavg_csmtbietaTT[3]; //mean/rmsEff/rmsGaus
675    
676     TH1F *hh_res_wtavg_csmtbietaTTAbs[3]; //mean/rmsEff/rmsGaus
677    
678     TH1F *hh_res_wtavg_csmallietaTTAbs[3]; //mean/rmsEff/rmsGaus
679    
680    
681     TH1F *hh_res_wtavg_csmcoietaTT[3]; //mean/rmsEff/rmsGaus
682     TH1F *hh_res_wtavg_csmcoietaTTAbs[3]; //mean/rmsEff/rmsGaus
683    
684     TH1F *hh_diff_csmtb_ietaTTAbs[50][50][17];
685     TH1F *hh_diff_csmtb_ietaTT[50][50][34];
686     TH1F *hh_diff_csmco_ietaTTAbs[50][50][17];
687     TH1F *hh_diff_csmco_ietaTT[50][50][34];
688     TH1F *hh_diff_csmall_ietaTT[50][50][34];
689     TH1F *hh_diff_csmall_ietaTTAbs[50][50][17];
690    
691     TH1F *hh_diff_csm[50][50][36][3]; //by each SM, all/central/outer
692    
693    
694     TH1F *hh_res_csm[51][3][3];
695    
696     TH2F *hh2_c[51];
697     TH2F *hh2_c_period[51];
698    
699    
700     float xbinLow[35];
701     for(int j=0; j< 35; j++){
702     float xbin = -85 + 5*j;
703     xbinLow[j] = xbin;
704     }
705     float xbinLow1[35];
706     for(int j=0; j<= 17 ; j++){
707     float xbin = 5 * j;
708     xbinLow1[j] = xbin+1;
709     }
710    
711    
712     TH1F *hh_c_ietaAbs_period[50][85];
713     TH1F *hh_c_ietaTTAbs_period[50][17];
714     TH1F *hh_c_smtbietaTTAbs_period[50][17];
715    
716    
717     for(int n=0; n< nIC+1; n++){
718     for(int k=0; k< 85; k++){
719     TString histname = TString(Form("hh_c_ietaAbs_period_%d_%d",n,k));
720     hh_c_ietaAbs_period[n][k] = new TH1F(histname,histname,500,0,2);
721     }
722     for(int k=0; k< 17; k++){
723     TString histname = TString(Form("hh_c_ietaTTAbs_period_%d_%d",n,k));
724     hh_c_ietaTTAbs_period[n][k] = new TH1F(histname,histname,500,0,2);
725     histname = TString(Form("hh_c_smtbietaTTAbs_period_%d_%d",n,k));
726     hh_c_smtbietaTTAbs_period[n][k] = new TH1F(histname,histname,500,0,2);
727     }
728     }
729    
730    
731    
732    
733     TH1F *hh_res_cietaAbs_period[50][3];
734     TH1F *hh_res_cietaTTAbs_period[50][3];
735     TH1F *hh_res_csmtbietaTTAbs_period[50][3];
736     for(int n=0; n< nIC+1; n++){
737     for(int k=0;k<3; k++){
738     TString histname = TString(Form("hh_res_cietaAbs_period_%d_%d",n,k));
739     hh_res_cietaAbs_period[n][k] = new TH1F(histname,histname,85,1,86);
740     histname = TString(Form("hh_res_cietaTTAbs_period_%d_%d",n,k));
741     hh_res_cietaTTAbs_period[n][k] = new TH1F(histname,histname,17,xbinLow1);
742     histname = TString(Form("hh_res_csmtbietaTTAbs_period_%d_%d",n,k));
743     hh_res_csmtbietaTTAbs_period[n][k] = new TH1F(histname,histname,17,xbinLow1);
744     }
745     }
746    
747    
748    
749     for(int n=0; n< nIC; n++){
750     TString histname = TString(Form("hh2_c_period_%d",n));
751     hh2_c_period[n] = new TH2F(histname,histname,171,-85,86,360,1,361);
752    
753     }
754    
755    
756     for(int n=0; n< nConstantSet+1; n++){
757     TString histname = TString(Form("hh2_c_%d",n));
758     hh2_c[n] = new TH2F(histname,histname,171,-85,86,360,1,361);
759    
760     }
761    
762     for(int n=0; n< nConstantSet+1; n++){
763     for(int j=0;j<36;j++){
764     for(int k=0;k<3; k++){
765     TString histname = TString(Form("hh_csm_%d_%d_%d",n,j,k));
766     hh_csm[n][j][k] = new TH1F(histname,histname,500,0,2);
767     }
768     }
769    
770     for(int j=0;j<3; j++){
771     for(int k=0;k<3; k++){
772     TString histname = TString(Form("hh_res_csm_%d_%d_%d",n,j,k));
773     hh_res_csm[n][j][k] = new TH1F(histname,histname,36,1,37);
774     }
775     }
776     }
777    
778    
779    
780     for(int n=0; n< nConstantSet; n++){
781    
782    
783     for(int j=0; j<17; j++){
784     TString histname = TString(Form("hh_csmtb_ietaTTAbs_%d_%d",n,j));
785     hh_csmtb_ietaTTAbs[n][j] = new TH1F(histname,histname,500,0,2);
786     histname = TString(Form("hh_csmco_ietaTTAbs_%d_%d",n,j));
787     hh_csmco_ietaTTAbs[n][j] = new TH1F(histname,histname,500,0,2);
788    
789     histname = TString(Form("hh_csmall_ietaTTAbs_%d_%d",n,j));
790     hh_csmall_ietaTTAbs[n][j] = new TH1F(histname,histname,500,0,2);
791    
792     }
793    
794    
795     for(int j=0; j<34; j++){
796     TString histname = TString(Form("hh_csmtb_ietaTT_%d_%d",n,j));
797     hh_csmtb_ietaTT[n][j] = new TH1F(histname,histname,500,0,2);
798     histname = TString(Form("hh_csmco_ietaTT_%d_%d",n,j));
799     hh_csmco_ietaTT[n][j] = new TH1F(histname,histname,500,0,2);
800     }
801     }
802    
803     for(int j=0; j<34; j++){
804     TString histname = TString(Form("hh_wtavg_csmtb_ietaTT_%d",j));
805     hh_wtavg_csmtb_ietaTT[j] = new TH1F(histname,histname,500,0,2);
806     histname = TString(Form("hh_wtavg_csmco_ietaTT_%d",j));
807     hh_wtavg_csmco_ietaTT[j] = new TH1F(histname,histname,500,0,2);
808     }
809     for(int j=0; j<17; j++){
810     TString histname = TString(Form("hh_wtavg_csmtb_ietaTTAbs_%d",j));
811     hh_wtavg_csmtb_ietaTTAbs[j] = new TH1F(histname,histname,500,0,2);
812     histname = TString(Form("hh_wtavg_csmco_ietaTTAbs_%d",j));
813     hh_wtavg_csmco_ietaTTAbs[j] = new TH1F(histname,histname,500,0,2);
814    
815     histname = TString(Form("hh_wtavg_csmall_ietaTTAbs_%d",j));
816     hh_wtavg_csmall_ietaTTAbs[j] = new TH1F(histname,histname,500,0,2);
817    
818     }
819    
820    
821     TH1F *hh_c_deadflag[50][20];
822    
823     TH1F *hh_c_deadflag_period[50][20];
824     TH1F *hh_cc_deadflag_period[50][20];
825     TH1F *hh_cc1_deadflag_period[50][20];
826    
827     for(int j=0; j< nConstantSet; j++){
828     for(int k=0; k<20; k++){
829     TString histname = TString(Form("hh_c_deadflag_%d_%d",j,k));
830     hh_c_deadflag[j][k] = new TH1F(histname,histname,500,0,2);
831    
832     histname = TString(Form("hh_c_deadflag_period_%d_%d",j,k));
833     hh_c_deadflag_period[j][k] = new TH1F(histname,histname,500,0,2);
834    
835     histname = TString(Form("hh_cc_deadflag_period_%d_%d",j,k));
836     hh_cc_deadflag_period[j][k] = new TH1F(histname,histname,500,0,2);
837     histname = TString(Form("hh_cc1_deadflag_period_%d_%d",j,k));
838     hh_cc1_deadflag_period[j][k] = new TH1F(histname,histname,500,0,2);
839    
840     }
841     }
842    
843    
844    
845    
846     for(int j=0; j< nConstantSet; j++){
847     for(int k=0; k<3; k++){
848    
849     TString histname = TString(Form("hh_res_csmtbietaTT_%d_%d",j,k));
850     hh_res_csmtbietaTT[j][k] = new TH1F(histname,histname,34,xbinLow);
851     histname = TString(Form("hh_res_csmtbietaTTAbs_%d_%d",j,k));
852     hh_res_csmtbietaTTAbs[j][k] = new TH1F(histname,histname,17,xbinLow1);
853    
854     histname = TString(Form("hh_res_csmcoietaTT_%d_%d",j,k));
855     hh_res_csmcoietaTT[j][k] = new TH1F(histname,histname,34,xbinLow);
856     histname = TString(Form("hh_res_csmcoietaTTAbs_%d_%d",j,k));
857     hh_res_csmcoietaTTAbs[j][k] = new TH1F(histname,histname,17,xbinLow1);
858    
859     }
860     }
861    
862     for(int k=0; k<3; k++){
863     TString histname = TString(Form("hh_res_wtavg_csmtbietaTT_%d",k));
864     hh_res_wtavg_csmtbietaTT[k] = new TH1F(histname,histname,34,xbinLow);
865     histname = TString(Form("hh_res_wtavg_csmtbietaTTAbs_%d",k));
866     hh_res_wtavg_csmtbietaTTAbs[k] = new TH1F(histname,histname,17,xbinLow1);
867    
868     histname = TString(Form("hh_res_wtavg_csmallietaTTAbs_%d",k));
869     hh_res_wtavg_csmallietaTTAbs[k] = new TH1F(histname,histname,17,xbinLow1);
870    
871     histname = TString(Form("hh_res_wtavg_csmcoietaTT_%d",k));
872     hh_res_wtavg_csmcoietaTT[k] = new TH1F(histname,histname,34,xbinLow);
873     histname = TString(Form("hh_res_wtavg_csmcoietaTTAbs_%d",k));
874     hh_res_wtavg_csmcoietaTTAbs[k] = new TH1F(histname,histname,17,xbinLow1);
875     }
876     TH2F *hh2_diff[50][50];
877    
878    
879     TH1F *hh_res_diff_csmtbietaTTAbs[50][50][3];
880     TH1F *hh_res_diff_csmcoietaTTAbs[50][50][3];
881     TH1F *hh_res_diff_csmtbietaTT[50][50][3];
882     TH1F *hh_res_diff_csmcoietaTT[50][50][3];
883     TH1F *hh_res_diff_csmallietaTT[50][50][3];
884     TH1F *hh_res_diff_csmallietaTTAbs[50][50][3];
885    
886     for(int n=0; n< nConstantSet; n++){
887     for(int k=n+1; k< nConstantSet; k++){
888     TString histname = TString (Form("hh2_diff_%dand%d",n,k));
889     hh2_diff[n][k] = new TH2F(histname,histname,171,-85,86,360,1,361);
890    
891     for(int j=0; j<36; j++){
892     for(int m=0; m<3; m++){
893     histname = TString ( Form("hh_diff_csm_%dand%d_%d_%d",n,k,j,m));
894     hh_diff_csm[n][k][j][m] = new TH1F(histname,histname,50,-0.1,0.1);
895     }
896     }
897    
898     for(int j=0; j<17; j++){
899     histname = TString (Form("hh_diff_csmtb_ietaTTAbs_%dand%d_%d",n,k,j));
900     hh_diff_csmtb_ietaTTAbs[n][k][j] = new TH1F(histname,histname,500,-0.1,0.1);
901     histname = TString (Form("hh_diff_csmco_ietaTTAbs_%dand%d_%d",n,k,j));
902     hh_diff_csmco_ietaTTAbs[n][k][j] = new TH1F(histname,histname,500,-0.1,0.1);
903     histname = TString (Form("hh_diff_csmall_ietaTTAbs_%dand%d_%d",n,k,j));
904     hh_diff_csmall_ietaTTAbs[n][k][j] = new TH1F(histname,histname,500,-0.1,0.1);
905    
906     }
907     for(int j=0; j<34; j++){
908     histname = TString (Form("hh_diff_csmtb_ietaTT_%dand%d_%d",n,k,j));
909     hh_diff_csmtb_ietaTT[n][k][j] = new TH1F(histname,histname,500,-0.1,0.1);
910     histname = TString (Form("hh_diff_csmco_ietaTT_%dand%d_%d",n,k,j));
911     hh_diff_csmco_ietaTT[n][k][j] = new TH1F(histname,histname,500,-0.1,0.1);
912    
913     histname = TString (Form("hh_diff_csmall_ietaTT_%dand%d_%d",n,k,j));
914     hh_diff_csmall_ietaTT[n][k][j] = new TH1F(histname,histname,500,-0.1,0.1);
915    
916     }
917     for(int j=0;j<3; j++){
918     histname = TString(Form("hh_res_diff_csmtbietaTTAbs_%dand%d_%d",n,k,j));
919     hh_res_diff_csmtbietaTTAbs[n][k][j] = new TH1F(histname,histname,17,xbinLow1);
920     histname = TString(Form("hh_res_diff_csmcoietaTTAbs_%dand%d_%d",n,k,j));
921     hh_res_diff_csmcoietaTTAbs[n][k][j] = new TH1F(histname,histname,17,xbinLow1);
922     histname = TString(Form("hh_res_diff_csmtbietaTT_%dand%d_%d",n,k,j));
923     hh_res_diff_csmtbietaTT[n][k][j] = new TH1F(histname,histname,34,xbinLow);
924     histname = TString(Form("hh_res_diff_csmcoietaTT_%dand%d_%d",n,k,j));
925     hh_res_diff_csmcoietaTT[n][k][j] = new TH1F(histname,histname,34,xbinLow);
926    
927     histname = TString(Form("hh_res_diff_csmallietaTT_%dand%d_%d",n,k,j));
928     hh_res_diff_csmallietaTT[n][k][j] = new TH1F(histname,histname,34,xbinLow);
929     histname = TString(Form("hh_res_diff_csmallietaTTAbs_%dand%d_%d",n,k,j));
930     hh_res_diff_csmallietaTTAbs[n][k][j] = new TH1F(histname,histname,17,xbinLow1);
931     }
932    
933     }
934     }
935    
936    
937     TH2F *hh2_largeICdiff[50][50];
938     for(int n=0; n< nIC; n++){
939     for(int k=n+1; k< nIC; k++){
940     TString histname = TString (Form("hh2_largeICdiff_%dand%d",n,k));
941     hh2_largeICdiff[n][k] =new TH2F(histname,histname,171,-85,86,360,1,361);
942     }
943     }
944     TH2F *hh2_largeICdiff_all =new TH2F("hh2_largeICdiff_all","hh2_largeICdiff_all",171,-85,86,360,1,361);
945    
946     vector<string> inputfileStat;
947    
948     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag60.pe1.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
949     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag61.pe1.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
950     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag62.pe1.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
951     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag63.pe1.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
952     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag72.pe1.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
953     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag73.pe1.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
954     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag64.pe1.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
955    
956     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag250.pe2.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
957     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag251.pe2.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
958     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag252.pe2.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
959     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag253.pe2.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
960     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag227.pe2.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
961     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag228.pe2.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
962     inputfileStat.push_back("calibres/deriveCalibConst.testCalibv1.dflag254.pe2.cut0.rmOvlap0.step11.method1.corrEta21.corrPhi10.corrSM1.corrDead-1.precalib1.vtx1.encorr0.evtNot-1.trig0.root");
963    
964    
965     TH1F *hh_res_ietaTTAbs[50][5];////[4] means the stat error.
966    
967     for(int n=0; n< int(inputfileStat.size()); n++){
968     for(int k=0;k<5;k++){
969     string filename = string(Form("hh_res_ietaTTAbs_%d_%d",n,k));
970     hh_res_ietaTTAbs[n][k] =new TH1F(filename.c_str(),filename.c_str(),17,xbinLow1);
971     }
972     }
973    
974     TH1F *hh_statErr_ietaTTAbs_period[50]; //pi0&eta combined for each period
975     for(int j=0; j< nIC; j++){
976     string filename = string(Form("hh_statErr_ietaTTAbs_period_%d",j));
977     hh_statErr_ietaTTAbs_period[j]=new TH1F(filename.c_str(),filename.c_str(),17,xbinLow1);
978     }
979    
980    
981     ///using MC-based forumula + 0.5/1 % sys.
982     // 7.4/resolution * 17 /sqrt(N) * sqrt( 1+ 1.8/sob) + 0.5/1%
983    
984    
985     for(int j=0; j< int(inputfileStat.size());j++){
986     string filename = inputfileStat[j];
987     TFile *f1 = new TFile(filename.c_str(),"read");
988     for(int k=0;k<4;k++){
989     string histname = string (Form("hh_res_ietaTTAbs_%d",k));
990     TH1F *hhtmp = (TH1F*)f1->Get(histname.c_str());
991     hh_res_ietaTTAbs[j][k]->Add(hhtmp);
992     }
993     }
994    
995     for(int j=0; j< int(inputfileStat.size());j++){
996    
997     for(int n=0; n< 17; n++){
998     float npiz = hh_res_ietaTTAbs[j][1]->GetBinContent(n+1)*0.5/(360*10) ;
999     float sob = hh_res_ietaTTAbs[j][2]->GetBinContent(n+1);
1000     float reso = hh_res_ietaTTAbs[j][3]->GetBinContent(n+1)/hh_res_ietaTTAbs[j][0]->GetBinContent(n+1) * 100;
1001     float statErr = reso/7.4 * 17/sqrt(npiz) * sqrt( 1+ 1.8/sob);
1002     hh_res_ietaTTAbs[j][4]->SetBinContent(n+1,statErr);
1003     }
1004     }
1005    
1006     ///the combined IC from all ICs
1007     ofstream txtout("interCalibConstants.combinedPi0EtaAllPeriod.EcalBarrel.txt",ios::out);
1008    
1009    
1010     ofstream txtout_period[50]; //pi0eta combined for each period
1011     for(int j=0; j< nIC; j++){
1012     string filename = string(Form("interCalibConstants.combinedPi0EtaPeriod%d.EcalBarrel.txt",j));
1013     txtout_period[j].open(filename.c_str(),ios::out);
1014     }
1015    
1016     ofstream txtout_periodv1[50]; //pi0eta combined for each period / averaged all
1017     for(int j=0; j< nIC; j++){
1018     string filename = string(Form("interCalibConstants.combinedPi0EtaPeriod%d.EcalBarrel.txtv1",j));
1019     txtout_periodv1[j].open(filename.c_str(),ios::out);
1020     }
1021    
1022    
1023     cout<<"fill" <<endl;
1024    
1025    
1026    
1027    
1028     for(int n=0; n< nConstantSet; n++){
1029    
1030     for(int j=0; j<170; j++){
1031    
1032     int ieta = j-85;
1033     if( ieta >=0) ieta += 1;
1034    
1035     for(int k=0; k<360; k++){
1036     int iphi = k;
1037     if( k==0) iphi = 360;
1038     int ietaTTAbs = (abs(ieta)-1)/5;
1039     int iSM = (iphi-1)/20+1;
1040     if( ieta<0) iSM += 18;
1041     int smTB = isTestBeamSM(iSM);
1042     int beta = ieta+85+1;
1043     int bphi = iphi;
1044     float c = icv1[n][j][k]; ///sm normalized to unit
1045     if( c>0){
1046    
1047     int ndflag = -1;
1048     ndflag = ndeadflagietaiphi_ic[n][j][k];
1049    
1050     if( ndflag <0){
1051     cout<<"wrong dflag? " << ndflag <<" "<<n<<" "<<j<<" "<<k<<endl;
1052     return;
1053     }
1054     hh_c_deadflag[n][ndflag]->Fill(c);
1055    
1056    
1057     hh2_c[n]->SetBinContent(beta,bphi,ic[n][j][k]);
1058    
1059     int ietaflag[3] = {abs(ieta)<25,abs(ieta)>=25 && abs(ieta)<60,abs(ieta)>=60};
1060     for(int jj=0;jj<3; jj++){
1061     if( ietaflag[jj]==0) continue;
1062     hh_csm[n][iSM-1][jj]->Fill(c);
1063     }
1064    
1065     if( smTB){
1066     hh_csmtb_ietaTTAbs[n][ietaTTAbs]->Fill(c);
1067     hh_csmtb_ietaTT[n][j/5]->Fill(c);
1068    
1069     bool etaCracks;
1070     bool phiCracks;
1071     isAtEcalBarrelModuleCracks(ieta,iphi,etaCracks,phiCracks);
1072     int crackFlag[5] = {1,!etaCracks && !phiCracks, etaCracks, phiCracks, etaCracks || phiCracks};
1073     for(int t =0; t<5; t++){
1074     if(crackFlag[t]==0) continue;
1075     hh_csmtb_ietaMod12[n][t]->Fill(c);
1076     }
1077    
1078     }else{
1079     hh_csmco_ietaTTAbs[n][ietaTTAbs]->Fill(c);
1080     hh_csmco_ietaTT[n][j/5]->Fill(c);
1081     }
1082     }else{
1083     hh2_c[n]->SetBinContent(beta,bphi,-1);
1084     }
1085    
1086     }
1087    
1088     }
1089     }
1090    
1091    
1092     for(int j=0; j<170; j++){
1093    
1094     int ieta = j-85;
1095     if( ieta >=0) ieta += 1;
1096    
1097     for(int k=0; k<360; k++){
1098     int iphi = k;
1099     if( k==0) iphi = 360;
1100     int ietaTTAbs = (abs(ieta)-1)/5;
1101     int iSM = (iphi-1)/20+1;
1102     if( ieta<0) iSM += 18;
1103     int smTB = isTestBeamSM(iSM);
1104     int beta = ieta+85+1;
1105     int bphi = iphi;
1106    
1107     for(int n1 =0; n1 < nConstantSet; n1++){
1108     for(int n2 =n1+1; n2 < nConstantSet; n2++){
1109     if(ic[n1][j][k] >0 && ic[n2][j][k] > 0){
1110    
1111     float diff = ic[n1][j][k] - ic[n2][j][k];
1112     float average = 0.5*(ic[n1][j][k] + ic[n2][j][k]);
1113     diff /= average;
1114    
1115     hh2_diff[n1][n2] ->SetBinContent(beta,bphi, diff);
1116    
1117     if(smTB){
1118     hh_diff_csmtb_ietaTTAbs[n1][n2][ietaTTAbs]->Fill(diff);
1119     hh_diff_csmtb_ietaTT[n1][n2][j/5]->Fill(diff);
1120     }else{
1121     hh_diff_csmco_ietaTTAbs[n1][n2][ietaTTAbs]->Fill(diff);
1122     hh_diff_csmco_ietaTT[n1][n2][j/5]->Fill(diff);
1123     }
1124     hh_diff_csmall_ietaTT[n1][n2][j/5]->Fill(diff);
1125     hh_diff_csmall_ietaTTAbs[n1][n2][ietaTTAbs]->Fill(diff);
1126    
1127     }else{
1128     hh2_diff[n1][n2] ->SetBinContent(beta,bphi,-1);
1129     }
1130     }
1131     }
1132     }
1133     }
1134    
1135    
1136    
1137     cout<<"fitting "<<endl;
1138     double resfit[10];
1139    
1140    
1141    
1142     for(int n1 =0; n1 < nConstantSet; n1++){
1143     for(int n2 =n1+1; n2 < nConstantSet; n2++){
1144    
1145     for(int j=0;j<34; j++){
1146     fitHistogram(hh_diff_csmtb_ietaTT[n1][n2][j],resfit);
1147     for(int n=0;n<3; n++){
1148     hh_res_diff_csmtbietaTT[n1][n2][n]->SetBinContent(j+1,resfit[2*n]);
1149     hh_res_diff_csmtbietaTT[n1][n2][n]->SetBinError(j+1,resfit[2*n+1]);
1150     }
1151     fitHistogram(hh_diff_csmco_ietaTT[n1][n2][j],resfit);
1152     for(int n=0;n<3; n++){
1153     hh_res_diff_csmcoietaTT[n1][n2][n]->SetBinContent(j+1,resfit[2*n]);
1154     hh_res_diff_csmcoietaTT[n1][n2][n]->SetBinError(j+1,resfit[2*n+1]);
1155     }
1156     fitHistogram(hh_diff_csmall_ietaTT[n1][n2][j],resfit);
1157     for(int n=0;n<3; n++){
1158     hh_res_diff_csmallietaTT[n1][n2][n]->SetBinContent(j+1,resfit[2*n]);
1159     hh_res_diff_csmallietaTT[n1][n2][n]->SetBinError(j+1,resfit[2*n+1]);
1160     }
1161    
1162     }
1163     for(int j=0;j<17; j++){
1164     fitHistogram(hh_diff_csmtb_ietaTTAbs[n1][n2][j],resfit);
1165     for(int n=0;n<3; n++){
1166     hh_res_diff_csmtbietaTTAbs[n1][n2][n]->SetBinContent(j+1,resfit[2*n]);
1167     hh_res_diff_csmtbietaTTAbs[n1][n2][n]->SetBinError(j+1,resfit[2*n+1]);
1168     }
1169     fitHistogram(hh_diff_csmco_ietaTTAbs[n1][n2][j],resfit);
1170     for(int n=0;n<3; n++){
1171     hh_res_diff_csmcoietaTTAbs[n1][n2][n]->SetBinContent(j+1,resfit[2*n]);
1172     hh_res_diff_csmcoietaTTAbs[n1][n2][n]->SetBinError(j+1,resfit[2*n+1]);
1173     }
1174     fitHistogram(hh_diff_csmall_ietaTTAbs[n1][n2][j],resfit);
1175     for(int n=0;n<3; n++){
1176     hh_res_diff_csmallietaTTAbs[n1][n2][n]->SetBinContent(j+1,resfit[2*n]);
1177     hh_res_diff_csmallietaTTAbs[n1][n2][n]->SetBinError(j+1,resfit[2*n+1]);
1178     }
1179    
1180     }
1181    
1182     }
1183     }
1184    
1185    
1186     for(int k=0;k<nConstantSet; k++){
1187    
1188    
1189     for(int j=0;j<34; j++){
1190     fitHistogram(hh_csmtb_ietaTT[k][j],resfit);
1191     for(int n=0;n<3; n++){
1192     hh_res_csmtbietaTT[k][n]->SetBinContent(j+1,resfit[2*n]);
1193     hh_res_csmtbietaTT[k][n]->SetBinError(j+1,resfit[2*n+1]);
1194     }
1195     }
1196     for(int j=0;j<17; j++){
1197     fitHistogram(hh_csmtb_ietaTTAbs[k][j],resfit);
1198     for(int n=0;n<3; n++){
1199     hh_res_csmtbietaTTAbs[k][n]->SetBinContent(j+1,resfit[2*n]);
1200     hh_res_csmtbietaTTAbs[k][n]->SetBinError(j+1,resfit[2*n+1]);
1201     }
1202     }
1203    
1204     for(int j=0;j<34; j++){
1205     fitHistogram(hh_csmco_ietaTT[k][j],resfit);
1206     for(int n=0;n<3; n++){
1207     hh_res_csmcoietaTT[k][n]->SetBinContent(j+1,resfit[2*n]);
1208     hh_res_csmcoietaTT[k][n]->SetBinError(j+1,resfit[2*n+1]);
1209     }
1210     }
1211     for(int j=0;j<17; j++){
1212     fitHistogram(hh_csmco_ietaTTAbs[k][j],resfit);
1213     for(int n=0;n<3; n++){
1214     hh_res_csmcoietaTTAbs[k][n]->SetBinContent(j+1,resfit[2*n]);
1215     hh_res_csmcoietaTTAbs[k][n]->SetBinError(j+1,resfit[2*n+1]);
1216     }
1217     }
1218    
1219     }
1220    
1221     cout<<"combining " <<endl;
1222    
1223    
1224     float wtSumC_period[50] = {0};
1225     float wtSumS_period[50] = {0};
1226    
1227    
1228    
1229     for(int j=0; j<170; j++){
1230     int ieta = j-85;
1231     if( ieta >=0) ieta += 1;
1232     for(int k=0; k<360; k++){
1233     int iphi = k;
1234     if( k==0) iphi = 360;
1235     int ietaTTAbs = (abs(ieta)-1)/5;
1236     int iSM = (iphi-1)/20+1;
1237     if( ieta<0) iSM += 18;
1238     int smTB = isTestBeamSM(iSM);
1239     int beta = ieta+85+1;
1240     int bphi = iphi;
1241    
1242    
1243     float wtSumC = 0;
1244     float wtSumS = 0;
1245    
1246     for(int n=0; n< nIC; n++){ ///for each period
1247     wtSumC_period[n] = 0;
1248     wtSumS_period[n] = 0;
1249     }
1250    
1251    
1252     for(int n=0; n< nConstantSet; n++){
1253     float c = ic[n][j][k];
1254    
1255    
1256    
1257    
1258     //float sigma = hh_res_csmtbietaTTAbs[n][2]->GetBinContent(ietaTTAbs+1);
1259     //if( sigma > sigmaTB){
1260     //sigma = sqrt( sigma * sigma - sigmaTB * sigmaTB);
1261     //}
1262     float statErr = 0;
1263     float sysErr = 0.5;
1264     if(abs(ieta)>=60) sysErr = 1;
1265    
1266     // if(n>=0 && n<=1){
1267     // statErr = hh_res_diff_csmallietaTTAbs[0][1][2]->GetBinContent(ietaTTAbs+1)/sqrt(2);
1268     // }else if( n>=2 && n<=3){
1269     // statErr = hh_res_diff_csmallietaTTAbs[2][3][2]->GetBinContent(ietaTTAbs+1)/sqrt(2);
1270     // }
1271    
1272    
1273     //now use MC-predicted precision
1274     statErr = hh_res_ietaTTAbs[n][4]->GetBinContent(ietaTTAbs+1);
1275    
1276     float sigma = sqrt( statErr * statErr + sysErr * sysErr);
1277    
1278     //if( ieta==1) cout<<"sigma " << ieta<<" n"<< n<<" "<< sigma <<" "<< statErr <<" "<< sysErr <<endl;
1279    
1280     if( c > 0){
1281    
1282     float tmp1 = c/ ( sigma * sigma);
1283     float tmp2 = 1/(sigma * sigma);
1284     int nperiod = n% nIC;
1285    
1286     wtSumC_period[nperiod] += tmp1;
1287     wtSumS_period[nperiod] += tmp2;
1288    
1289    
1290     ////// Set by hand, for some crystals'IC, not used for combination
1291     if( n==4 || n==11 || n==5 || n== 12){ //// fix by hand September 's IC in a region ( deadflag txt file was wrong )
1292     //set by hand those to be -1 for combination, including +2 crystals next to the "dead region"
1293     if( (ieta >=-22 && ieta <= -12) &&
1294     (iphi >=279 && iphi <= 292 )
1295     ){
1296     c = -1;
1297     }
1298     }
1299     if( n==4 || n==11 || n==5 || n== 12){ ///one more TT
1300     if( (ieta >=35 && ieta <= 41) &&
1301     (iphi >=310 && iphi <= 316 )
1302     ){
1303     c = -1;
1304     }
1305     }
1306     if( n==0 || n==1 || n==2 || n==3 || n ==4 || n==6 || n==7 || n==8 || n==9 || n==10|| n== 11 || n== 13){ //using only 5 and 12 for combination
1307     if( ieta >= -31 && ieta <= -25 && iphi >= 165 && iphi <= 171){
1308     c = -1;
1309     }
1310     }
1311     if( n==2 || n==3 || n==9 || n==10){
1312     if( ieta >= -81 && ieta <= -75 && iphi >= 300 && iphi <= 311){
1313     c = -1;
1314     }
1315     }
1316    
1317     if( c> 0){ //all combined
1318     wtSumC += tmp1;
1319     wtSumS += tmp2;
1320     }
1321     }
1322     }
1323    
1324    
1325     if( wtSumC > 0){ //combined all
1326     icwt[j][k] = wtSumC / wtSumS;
1327     hh2_c[nConstantSet]->SetBinContent(beta,bphi,icwt[j][k]);
1328     }else{
1329     icwt[j][k] = -1;
1330     hh2_c[nConstantSet]->SetBinContent(beta,bphi,-1);
1331     }
1332    
1333     for(int n=0; n< nIC; n++){ ///for each period
1334     if( wtSumC_period[n] > 0){
1335     icwt_period[n][j][k] = wtSumC_period[n] / wtSumS_period[n];
1336     hh2_c_period[n]->SetBinContent(beta,bphi,icwt_period[n][j][k]);
1337     }else{
1338     icwt_period[n][j][k] = -1;
1339     hh2_c_period[n]->SetBinContent(beta,bphi,-1);
1340     }
1341     }
1342    
1343     }
1344     }
1345    
1346    
1347    
1348     float statErr_allCombined[85] = {0} ;
1349    
1350     for(int b=1; b<= hh_res_ietaTTAbs[0][4]->GetNbinsX(); b++){
1351     float sumStatErr2 = 0;
1352     for(int n=0; n< nConstantSet; n++){
1353     float statErr = hh_res_ietaTTAbs[n][4]->GetBinContent(b);
1354     sumStatErr2 += 1./ ( statErr * statErr );
1355     }
1356     float statErr = sqrt( 1./ sumStatErr2 );
1357     statErr_allCombined[b-1] = statErr;
1358    
1359     cout<<" statErr_allCombined " << b <<" "<< statErr_allCombined[b-1] << " % "<<endl;
1360    
1361     }
1362    
1363    
1364     ///Stat error for each period
1365     for(int j=0; j< nIC; j++){
1366     for(int b=1; b<= hh_res_ietaTTAbs[j][4]->GetNbinsX();b++){
1367     float statErrPi0 = hh_res_ietaTTAbs[j][4]->GetBinContent(b);
1368     float statErrEta = hh_res_ietaTTAbs[j+nIC][4]->GetBinContent(b);
1369     float statErr = 1/(statErrPi0*statErrPi0) + 1./(statErrEta*statErrEta);
1370     statErr = sqrt(1/statErr);
1371     hh_statErr_ietaTTAbs_period[j]->SetBinContent(b,statErr);
1372     }
1373     }
1374    
1375    
1376     txtoutTocheck <<"checkme " <<endl;
1377    
1378     for(int j=0; j<170; j++){
1379    
1380     int ieta = j-85;
1381     if( ieta >=0) ieta += 1;
1382    
1383     for(int k=0; k<360; k++){
1384     int iphi = k;
1385     if( k==0) iphi = 360;
1386     int beta = ieta+85+1;
1387     int bphi = iphi;
1388     int ietaTTAbs = (abs(ieta)-1)/5;
1389    
1390     bool largeIC = false;
1391    
1392     for(int n=0; n< nIC; n++){ ///for each period
1393     if( icwt_period[n][j][k] > 0 && (icwt_period[n][j][k] >1.2 || icwt_period[n][j][k] < 0.8)){
1394     largeIC = true;
1395     }
1396     }
1397     if(largeIC){
1398     txtoutTocheck<<"largeIC "<< j<<" "<<k <<" "<<ieta <<" "<<iphi<<" " ;
1399     for(int n1 =0; n1 < nIC ; n1++){
1400     txtoutTocheck<<icwt_period[n1][j][k]<<" ";
1401     }
1402     txtoutTocheck<<endl;
1403     }
1404    
1405     bool checkICdiff = false;
1406     for(int n1=0; n1< nIC; n1++){ ///for each period
1407     for(int n2=n1+1;n2< nIC; n2++){ ///for each period
1408     if( icwt_period[n1][j][k]>0 && icwt_period[n2][j][k]>0){
1409     float reldiff = fabs( icwt_period[n1][j][k] - icwt_period[n2][j][k] )/ ( 0.5* (icwt_period[n1][j][k] + icwt_period[n2][j][k]));
1410    
1411     float statErr1 = hh_statErr_ietaTTAbs_period[n1]->GetBinContent(ietaTTAbs+1)/100;
1412     float statErr2 = hh_statErr_ietaTTAbs_period[n2]->GetBinContent(ietaTTAbs+1)/100;
1413    
1414     ////assuming same sys.
1415     float diff_statErr = sqrt( statErr1*statErr1 + statErr2*statErr2);
1416    
1417     checkICdiff = false;
1418     // if( abs(ieta)<60 && reldiff >0.02){
1419     // checkICdiff = true;
1420     // }
1421     // if( abs(ieta)>=60 && reldiff >0.1){
1422     // checkICdiff = true;
1423     // }
1424     checkICdiff = reldiff > 3 * diff_statErr ;
1425     if(checkICdiff){
1426     hh2_largeICdiff[n1][n2]->SetBinContent(beta,bphi,1);
1427     }else{
1428     hh2_largeICdiff[n1][n2]->SetBinContent(beta,bphi,-1);
1429     }
1430    
1431     }else{
1432     hh2_largeICdiff[n1][n2]->SetBinContent(beta,bphi,-1);
1433     }
1434     }
1435     }
1436     if(checkICdiff){
1437     hh2_largeICdiff_all->SetBinContent(beta,bphi,1);
1438     txtoutTocheck<<"checkICdiff "<< j<<" "<<k <<" "<<ieta <<" "<<iphi<<" " ;
1439     for(int n1 =0; n1 < nIC ; n1++){
1440     txtoutTocheck<<icwt_period[n1][j][k]<<" ";
1441     }
1442     txtoutTocheck<<endl;
1443     }else{
1444     hh2_largeICdiff_all->SetBinContent(beta,bphi,-1);
1445     }
1446    
1447     }
1448     }
1449    
1450    
1451    
1452    
1453     scaleMeanToUnit(icwt);
1454     ///for estimating precision of combined IC
1455     copyConstant(icwt,icwtv1);
1456     NormSMScaleToUnit(icwtv1);
1457    
1458    
1459     for(int n=0; n< nIC; n++){ ///for each period
1460     scaleMeanToUnit(icwt_period[n]);
1461     }
1462    
1463    
1464     for(int j=0; j<170; j++){
1465     int ieta = j-85;
1466     if( ieta >=0) ieta += 1;
1467     for(int k=0; k<360; k++){
1468     int iphi = k;
1469     if( k==0) iphi = 360;
1470     int ietaTTAbs = (abs(ieta)-1)/5;
1471     int iSM = (iphi-1)/20+1;
1472     if( ieta<0) iSM += 18;
1473     int smTB = isTestBeamSM(iSM);
1474    
1475    
1476     float sysErr = 0.5;
1477     if(abs(ieta)>=60) sysErr = 1;
1478    
1479     float icErr = sysErr;
1480     float statErr = statErr_allCombined[ietaTTAbs];
1481     icErr = sqrt(sysErr *sysErr + statErr * statErr);
1482     icErr /=100;
1483    
1484    
1485     if( icwt[j][k] > 0){
1486     txtout<<ieta<<" "<<iphi<<" "<< icwt[j][k]*interCalib_preCalib[j][k]<<" "<< icErr* icwt[j][k]*interCalib_preCalib[j][k]<<endl;
1487     //txtout<<ieta<<" "<<iphi<<" "<< icwt[j][k]*interCalib_preCalib[j][k]<<endl;
1488     }else{
1489     txtout<<ieta<<" "<<iphi<<" "<< -1 <<" "<< 999 <<endl;
1490     //txtout<<ieta<<" "<<iphi<<" "<< interCalib_preCalib[j][k]<<endl;
1491    
1492     }
1493    
1494     float c = icwtv1[j][k]; ///sm normalized to unit
1495    
1496    
1497     int absieta = abs(ieta);
1498    
1499     for(int n=0; n< nIC+1; n++){ ///for each period and all combined
1500    
1501     if(n< nIC){
1502     if( icwt_period[n][j][k] > 0){
1503     hh_c_ietaAbs_period[n][absieta-1]->Fill(icwt_period[n][j][k]);
1504     hh_c_ietaTTAbs_period[n][ietaTTAbs]->Fill(icwt_period[n][j][k]);
1505    
1506     if( smTB){
1507     hh_c_smtbietaTTAbs_period[n][ietaTTAbs]->Fill(icwt_period[n][j][k]);
1508     }
1509     int deadflag = ndeadflagietaiphi_ic[n][j][k];
1510     if( deadflag<0){
1511     cout<<"wrong deadlfag !!! n " << n <<" "<<endl;
1512     exit(1);
1513     }
1514     hh_c_deadflag_period[n][deadflag]->Fill(icwt_period[n][j][k]);
1515     hh_cc_deadflag_period[n][deadflag]->Fill(icwt_period[n][j][k] * interCalib_preCalibTrue[j][k]);
1516     hh_cc1_deadflag_period[n][deadflag]->Fill(interCalib_preCalibTrue[j][k]);
1517     if(deadflag>0){
1518     hh_c_deadflag_period[n][19]->Fill(icwt_period[n][j][k]);
1519     hh_cc_deadflag_period[n][19]->Fill( interCalib_preCalib[j][k] * icwt_period[n][j][k]);
1520     hh_cc1_deadflag_period[n][19]->Fill( interCalib_preCalibTrue[j][k]);
1521     }
1522    
1523     }
1524     }else{
1525     if( icwt[j][k] > 0){
1526     hh_c_ietaAbs_period[n][absieta-1]->Fill(icwt[j][k]);
1527     hh_c_ietaTTAbs_period[n][ietaTTAbs]->Fill(icwt[j][k]);
1528    
1529     int deadflag = ndeadflagietaiphi_ic[0][j][k];
1530     if( deadflag<0){
1531     cout<<"wrong deadlfag !!! comb " << n <<" "<<endl;
1532     exit(1);
1533     }
1534     hh_c_deadflag_period[n][deadflag]->Fill(icwt[j][k]);
1535     if(deadflag>0){
1536     hh_c_deadflag_period[n][19]->Fill(icwt[j][k]);
1537     }
1538    
1539     if( smTB){
1540     hh_c_smtbietaTTAbs_period[n][ietaTTAbs]->Fill(icwt[j][k]);
1541     }
1542    
1543     }
1544     }
1545     }
1546    
1547    
1548     for(int n=0; n< nIC; n++){ ///for each period
1549     statErr = hh_statErr_ietaTTAbs_period[n]->GetBinContent(ietaTTAbs+1);
1550     icErr = sqrt( statErr * statErr + sysErr * sysErr);
1551     icErr /= 100;
1552    
1553     if( icwt_period[n][j][k] > 0){
1554     txtout_period[n]<<ieta<<" "<<iphi<<" "<< icwt_period[n][j][k]*interCalib_preCalib[j][k]<<" "<< icErr*icwt_period[n][j][k]*interCalib_preCalib[j][k] <<endl;
1555     //txtout_period[n]<<ieta<<" "<<iphi<<" "<< icwt_period[n][j][k]*interCalib_preCalib[j][k]<<endl;
1556     }else{
1557     txtout_period[n]<<ieta<<" "<<iphi<<" "<< -1<<" "<< 999 <<endl;
1558     //txtout_period[n]<<ieta<<" "<<iphi<<" "<< interCalib_preCalib[j][k]<<" "<< 999 <<endl;
1559    
1560     }
1561     }
1562    
1563    
1564    
1565     if( c>0){
1566    
1567     int ietaflag[3] = {abs(ieta)<25,abs(ieta)>=25 && abs(ieta)<60,abs(ieta)>=60};
1568     for(int jj=0;jj<3; jj++){
1569     if( ietaflag[jj]==0) continue;
1570     hh_csm[nConstantSet][iSM-1][jj]->Fill(c);
1571     }
1572    
1573     if( smTB){
1574     hh_wtavg_csmtb_ietaTT[j/5]->Fill(c);
1575     hh_wtavg_csmtb_ietaTTAbs[ietaTTAbs]->Fill(c);
1576     }else{
1577     hh_wtavg_csmco_ietaTT[j/5]->Fill(c);
1578     hh_wtavg_csmco_ietaTTAbs[ietaTTAbs]->Fill(c);
1579     }
1580    
1581     ///all sm
1582     hh_wtavg_csmall_ietaTTAbs[ietaTTAbs]->Fill(c);
1583    
1584     }
1585    
1586     }
1587     }
1588    
1589     cout <<" fitting combined " <<endl;
1590    
1591    
1592     for(int n=0; n< nConstantSet+1; n++){
1593     for(int j=0; j<36; j++){
1594    
1595     for(int k=0;k<3;k++){
1596     fitHistogram(hh_csm[n][j][k],resfit);
1597     for(int jj=0;jj<3; jj++){
1598     hh_res_csm[n][k][jj]->SetBinContent(j+1,resfit[2*jj]);
1599     hh_res_csm[n][k][jj]->SetBinError(j+1,resfit[2*jj+1]);
1600     }
1601     }
1602     }
1603     }
1604    
1605    
1606    
1607     for(int n1=0; n1< nIC+1; n1++){ ///for each period and all combined
1608     for(int j=0; j< 85; j++){
1609     fitHistogram(hh_c_ietaAbs_period[n1][j],resfit);
1610     for(int n=0;n<3; n++){
1611     hh_res_cietaAbs_period[n1][n]->SetBinContent(j+1,resfit[2*n]);
1612     hh_res_cietaAbs_period[n1][n]->SetBinError(j+1,resfit[2*n+1]);
1613     }
1614     }
1615     for(int j=0; j< 17; j++){
1616     fitHistogram(hh_c_ietaTTAbs_period[n1][j],resfit);
1617     for(int n=0;n<3; n++){
1618     hh_res_cietaTTAbs_period[n1][n]->SetBinContent(j+1,resfit[2*n]);
1619     hh_res_cietaTTAbs_period[n1][n]->SetBinError(j+1,resfit[2*n+1]);
1620     }
1621     fitHistogram(hh_c_smtbietaTTAbs_period[n1][j],resfit);
1622     for(int n=0;n<3; n++){
1623     hh_res_csmtbietaTTAbs_period[n1][n]->SetBinContent(j+1,resfit[2*n]);
1624     hh_res_csmtbietaTTAbs_period[n1][n]->SetBinError(j+1,resfit[2*n+1]);
1625     }
1626     }
1627     }
1628    
1629    
1630     for(int j=0;j<34; j++){
1631     fitHistogram(hh_wtavg_csmtb_ietaTT[j],resfit);
1632     for(int n=0;n<3; n++){
1633     hh_res_wtavg_csmtbietaTT[n]->SetBinContent(j+1,resfit[2*n]);
1634     hh_res_wtavg_csmtbietaTT[n]->SetBinError(j+1,resfit[2*n+1]);
1635     }
1636     fitHistogram(hh_wtavg_csmco_ietaTT[j],resfit);
1637     for(int n=0;n<3; n++){
1638     hh_res_wtavg_csmcoietaTT[n]->SetBinContent(j+1,resfit[2*n]);
1639     hh_res_wtavg_csmcoietaTT[n]->SetBinError(j+1,resfit[2*n+1]);
1640     }
1641     }
1642     for(int j=0;j<17; j++){
1643     fitHistogram(hh_wtavg_csmtb_ietaTTAbs[j],resfit);
1644     for(int n=0;n<3; n++){
1645     hh_res_wtavg_csmtbietaTTAbs[n]->SetBinContent(j+1,resfit[2*n]);
1646     hh_res_wtavg_csmtbietaTTAbs[n]->SetBinError(j+1,resfit[2*n+1]);
1647     }
1648    
1649    
1650     fitHistogram(hh_wtavg_csmall_ietaTTAbs[j],resfit);
1651     for(int n=0;n<3; n++){
1652     hh_res_wtavg_csmallietaTTAbs[n]->SetBinContent(j+1,resfit[2*n]);
1653     hh_res_wtavg_csmallietaTTAbs[n]->SetBinError(j+1,resfit[2*n+1]);
1654     }
1655    
1656     fitHistogram(hh_wtavg_csmco_ietaTTAbs[j],resfit);
1657     for(int n=0;n<3; n++){
1658     hh_res_wtavg_csmcoietaTTAbs[n]->SetBinContent(j+1,resfit[2*n]);
1659     hh_res_wtavg_csmcoietaTTAbs[n]->SetBinError(j+1,resfit[2*n+1]);
1660     }
1661     }
1662    
1663    
1664     fnew->Write();
1665     fnew->Close();
1666    
1667     }