ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/yangyong/combineICEB/getCalibConstants.cc
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: V2012combv2b, V2012combv2a, V2012combv2, V2012combv1, V2011combv1, pi0etaeb_laser20111122
Changes since 1.1: +1147 -0 lines
Log Message:
LaserTag:
EcalLaserAPDPNRatios_data_20111122_158851_180363

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

File Contents

# Content
1
2
3 void scaleMeanToUnit(float C[170][360]){
4
5
6 double mean = 0;
7 int n = 0;
8 for(int j=0; j< 170; j++){
9 for(int k=0; k< 360; k++){
10
11 if( C[j][k] >0){
12 mean += C[j][k];
13 n ++;
14 }
15 }
16 }
17 mean /= n;
18 cout<<"mean " << mean<<" "<<n<<endl;
19
20 for(int j=0; j< 170; j++){
21 for(int k=0; k< 360; k++){
22 if(C[j][k] >0){
23 C[j][k] /= mean;
24 }
25 }
26 }
27
28 }
29
30
31
32
33 // void scaleMeanToUnitEndcap(float C[2][101][101]){
34
35
36 // double mean = 0;
37 // int n = 0;
38 // for(int j=0; j< 170; j++){
39 // for(int k=0; k< 360; k++){
40
41 // if( C[j][k] >0){
42 // mean += C[j][k];
43 // n ++;
44 // }
45 // }
46 // }
47 // mean /= n;
48 // cout<<"mean " << mean<<" "<<n<<endl;
49
50 // for(int j=0; j< 170; j++){
51 // for(int k=0; k< 360; k++){
52 // if(C[j][k] >0){
53 // C[j][k] /= mean;
54 // }
55 // }
56 // }
57
58 // }
59
60
61
62
63 void NormIetaToUnitTestBeamSMOnly(float C[170][360]){
64
65 float mean_ieta[170] = {0};
66 int nmean_ieta[170] = {0};
67
68 float mean_ietaTT[170] = {0};
69 int nmean_ietaTT[170] = {0};
70
71
72 for(int j=0; j< 170; j++){
73
74 int ieta = j-85;
75 if( ieta >=0) ieta += 1;
76
77 for(int k=0; k< 360; k++){
78
79 int iphi = k;
80 if( k==0) iphi = 360;
81
82 int iSM = (iphi-1)/20+1;
83 if( ieta<0) iSM += 18;
84 int smTB = isTestBeamSM(iSM);
85 if( smTB ==0) continue;
86
87 //if( C[j][k] >0 ){
88 if( C[j][k] >0.5 && C[j][k] < 1.5 ){
89 mean_ieta[j] += C[j][k];
90 nmean_ieta[j] ++;
91
92 mean_ietaTT[j/5] += C[j][k];
93 nmean_ietaTT[j/5] ++;
94
95 }
96
97 }
98 }
99
100
101 for(int j=0; j< 170; j++){
102 mean_ieta[j] /= nmean_ieta[j];
103 //cout<<" mean_ieta[" << j<<"] "<< mean_ieta[j] <<endl;
104 }
105
106 // for(int j=0; j< 34; j++){
107 // mean_ietaTT[j] /= nmean_ietaTT[j];
108 // cout<<" mean_ietaTT[" << j<<"] "<< mean_ietaTT[j] <<endl;
109 // }
110
111 for(int j=0; j< 170; j++){
112 for(int k=0; k< 360; k++){
113
114 if( C[j][k] >0 ){
115 C[j][k] /= mean_ieta[j] ;
116 }
117 }
118 }
119
120 }
121
122
123
124 void NormIetaAbsToUnitTestBeamSMOnly(float C[170][360]){
125
126
127 float mean_ieta[170] = {0};
128 int nmean_ieta[170] = {0};
129
130 float mean_ietaTT[170] = {0};
131 int nmean_ietaTT[170] = {0};
132
133
134 for(int j=0; j< 170; j++){
135
136 int ieta = j-85;
137 if( ieta >=0) ieta += 1;
138
139 for(int k=0; k< 360; k++){
140
141 int iphi = k;
142 if( k==0) iphi = 360;
143
144 int iSM = (iphi-1)/20+1;
145 if( ieta<0) iSM += 18;
146 int smTB = isTestBeamSM(iSM);
147 if( smTB ==0) continue;
148
149 int ietaAbs = abs(ieta);
150
151 //if( C[j][k] >0 ){
152 if( C[j][k] >0.5 && C[j][k] < 1.5 ){
153 mean_ieta[ietaAbs-1] += C[j][k];
154 nmean_ieta[ietaAbs-1] ++;
155
156 mean_ietaTT[j/5] += C[j][k];
157 nmean_ietaTT[j/5] ++;
158
159 }
160 }
161 }
162
163
164 for(int j=0; j< 85; j++){
165 mean_ieta[j] /= nmean_ieta[j];
166 //cout<<" mean_ietaAbs[" << j<<"] "<< mean_ieta[j] <<endl;
167 }
168
169 // for(int j=0; j< 34; j++){
170 // mean_ietaTT[j] /= nmean_ietaTT[j];
171 // cout<<" mean_ietaTT[" << j<<"] "<< mean_ietaTT[j] <<endl;
172 // }
173
174 for(int j=0; j< 170; j++){
175
176 int ieta = j-85;
177 if( ieta >=0) ieta += 1;
178
179 for(int k=0; k< 360; k++){
180
181 int ietaAbs = abs(ieta);
182
183 if( C[j][k] >0 ){
184
185 C[j][k] /= mean_ieta[ietaAbs-1] ;
186
187 }
188 }
189 }
190
191 }
192
193
194
195 void NormCrystalDeadFlag_v1(float C[170][360], int ndeadflagietaiphi[170][360]){
196
197 float mean_deadflag[20] = {0};
198 float nmean_deadflag[20] = {0};
199
200 int ndeadCrystals = 0;
201 for(int j=0; j< 170; j++){
202 for(int k=0; k< 360; k++){
203 if( C[j][k] <0){
204 ndeadCrystals++;
205 continue;
206 }
207 int ndeadflag = ndeadflagietaiphi[j][k];
208 if( ndeadflag >=0){
209 mean_deadflag[ndeadflag] += C[j][k];
210 nmean_deadflag[ndeadflag] ++;
211 }
212 }
213 }
214 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
215 for(int j=0; j<20; j++){
216 if( nmean_deadflag[j] >0){
217 mean_deadflag[j] /= nmean_deadflag[j];
218 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
219 }
220 }
221 for(int j=0; j< 170; j++){
222 for(int k=0; k< 360; k++){
223 if( C[j][k] <0) continue;
224 int ndeadflag = ndeadflagietaiphi[j][k];
225 if( ndeadflag >=0){
226 C[j][k] /= mean_deadflag[ndeadflag];
227 }
228 }
229 }
230
231 }
232
233 void NormCrystalDeadFlag(float C[170][360]){
234
235 float mean_deadflag[20] = {0};
236 float nmean_deadflag[20] = {0};
237
238 int ndeadCrystals = 0;
239 for(int j=0; j< 170; j++){
240 for(int k=0; k< 360; k++){
241 if( C[j][k] <0){
242 ndeadCrystals++;
243 continue;
244 }
245 int ndeadflag = ndeadflag_ietaiphi[j][k];
246 if( ndeadflag >=0){
247 mean_deadflag[ndeadflag] += C[j][k];
248 nmean_deadflag[ndeadflag] ++;
249 }
250 }
251 }
252 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
253 for(int j=0; j<20; j++){
254 if( nmean_deadflag[j] >0){
255 mean_deadflag[j] /= nmean_deadflag[j];
256 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
257 }
258 }
259 for(int j=0; j< 170; j++){
260 for(int k=0; k< 360; k++){
261 if( C[j][k] <0) continue;
262 int ndeadflag = ndeadflag_ietaiphi[j][k];
263 if( ndeadflag >=0){
264 C[j][k] /= mean_deadflag[ndeadflag];
265 }
266 }
267 }
268
269 }
270
271
272
273 void SetSMScale(float C[170][360],TH1F *hhtmp){
274 float mean_ism[36] = {0};
275 float nmean_ism[36]= {0};
276
277 for(int j=0; j< 170; j++){
278 for(int k=0; k< 360; k++){
279
280 int ieta = j-85;
281 if( j >=0) ieta += 1;
282 int iphi = k;
283 if( k==0) iphi = 360;
284 int iSM = (iphi-1)/20+1;
285 if( ieta<0) iSM += 18;
286
287 if( C[j][k] >0 ){
288 mean_ism[iSM-1 ] += C[j][k];
289 nmean_ism[iSM-1] ++;
290 }
291 }
292 }
293
294 for(int j=0; j<36; j++){
295 mean_ism[j] /= nmean_ism[j];
296 /// cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl;
297 }
298
299
300 for(int j=0; j< 170; j++){
301 for(int k=0; k< 360; k++){
302
303 int ieta = j-85;
304 if( j >=0) ieta += 1;
305 int iphi = k;
306 if( k==0) iphi = 360;
307 int iSM = (iphi-1)/20+1;
308 if( ieta<0) iSM += 18;
309
310 if( C[j][k] >0 ){
311 C[j][k] *= hhtmp->GetBinContent(iSM);
312 }
313 }
314 }
315
316
317 }
318
319
320
321
322
323
324
325 void NormSMScaleToUnit(float C[170][360]){
326 float mean_ism[36] = {0};
327 float nmean_ism[36]= {0};
328
329 for(int j=0; j< 170; j++){
330 for(int k=0; k< 360; k++){
331
332 int ieta = j-85;
333 if( j >=0) ieta += 1;
334 int iphi = k;
335 if( k==0) iphi = 360;
336 int iSM = (iphi-1)/20+1;
337 if( ieta<0) iSM += 18;
338
339
340 if( C[j][k] >0 ){
341 mean_ism[iSM-1 ] += C[j][k];
342 nmean_ism[iSM-1] ++;
343 }
344 }
345 }
346
347 for(int j=0; j<36; j++){
348 mean_ism[j] /= nmean_ism[j];
349 //cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl;
350 }
351
352 for(int j=0; j< 170; j++){
353 for(int k=0; k< 360; k++){
354
355 int ieta = j-85;
356 if( j >=0) ieta += 1;
357 int iphi = k;
358 if( k==0) iphi = 360;
359 int iSM = (iphi-1)/20+1;
360 if( ieta<0) iSM += 18;
361
362 if( C[j][k] >0 ){
363 C[j][k] /= mean_ism[iSM-1];
364 }
365
366 }
367 }
368
369 }
370
371
372
373 void readInterCalibConstEBSimple(const char *input,float C[170][360]){
374
375 ifstream txtin1(input,ios::in);
376 if (txtin1.fail()){
377 cout<<"error open file barrel.. " << input<<endl;
378 exit(1);
379 }
380
381
382 int eta;
383 int phi;
384
385 float cc;
386
387 int n =0;
388 while(txtin1.good()){
389
390
391 txtin1 >> eta >> phi >> cc;
392
393 int iphi = phi;
394 int ieta = eta;
395
396 // cout<< eta <<" "<< phi <<" "<< cc <<endl;
397
398 convxtalid(phi,eta);
399 eta += 85;
400 C[eta][phi] = cc;
401
402 n++;
403 if( n>= 61200) break;
404
405 }
406
407 }
408
409
410
411
412
413 void readInterCalibConstEBSimplev1(const char *input,float C[170][360]){
414
415 ifstream txtin1(input,ios::in);
416 if (txtin1.fail()){
417 cout<<"error open file barrel.. " << input<<endl;
418 exit(1);
419 }
420
421
422 int eta;
423 int phi;
424
425 float cc;
426 double mean = 0;
427 int ngood = 0;
428 int n =0;
429 while(txtin1.good()){
430
431 txtin1 >> eta >> phi >> cc;
432
433 int iphi = phi;
434 int ieta = eta;
435
436 if( !(eta>=0 && eta <= 169 && phi >=0 && phi<=359)) {
437 cout<<"wrong eta/phi " << eta<<" "<<phi <<endl;
438 exit(1);
439 }
440
441 C[eta][phi] = cc;
442 if( cc >0){
443 mean += cc;
444 ngood ++;
445 }
446 n++;
447 if( n>= 61200) break;
448 }
449 mean /= ngood;
450 cout<<"read ngood "<< ngood <<" mean: "<< mean<<endl;
451
452 for(int j=0; j<170; j++){
453 for(int k=0; k<360;k++){
454 if(C[j][k] >0) {
455 C[j][k] /= mean;
456 }
457 }
458 }
459
460 }
461
462
463
464 void readInterCalibBarre_Phiv1new(){
465
466 char *file_input = new char[500];
467 sprintf(file_input,"phisym_EcalIntercalib-160M-corrected.txt");
468
469 sprintf(file_input,"EcalIntercalibConstants_Run2010A_corr.xml.txt");
470
471 cout<<"READING fro inter-calib file "<<file_input<<endl;
472 ifstream inputcc(file_input,ios::in);
473 if (inputcc.fail()){
474 cout<<"error open file barrel.. " << file_input<<endl;
475 exit(1);
476 }
477
478 float c1;
479 for(int ieta=-85; ieta <= 85; ieta++){
480 if( ieta==0) continue;
481 for(int iphi=1; iphi <= 360; iphi++){
482 inputcc >> c1;
483
484 int eta = ieta;
485 int phi = iphi;
486 convxtalid(phi,eta);
487
488 if( c1 <=0){
489 c1 =1;
490 }
491
492 if( c1 == 1 ){
493 c1 = - 1;
494 }
495 if( c1 <0.6 || c1 >2.5){
496 c1 = -1;
497 }
498
499
500 CphiCorr[eta+85][phi] = c1;
501 }
502 }
503
504 }
505
506
507 void readInterCalibBarre_Phi2011a(){
508
509 char *file_input = new char[500];
510
511
512 sprintf(file_input,"EcalIntercalib-2011A-corrected.xmlv1");
513
514 cout<<"READING fro inter-calib file "<<file_input<<endl;
515 ifstream inputcc(file_input,ios::in);
516 if (inputcc.fail()){
517 cout<<"error open file barrel.. " << file_input<<endl;
518 exit(1);
519 }
520
521 float c1;
522 for(int ieta=-85; ieta <= 85; ieta++){
523 if( ieta==0) continue;
524 for(int iphi=1; iphi <= 360; iphi++){
525 inputcc >> c1;
526
527 int eta = ieta;
528 int phi = iphi;
529 convxtalid(phi,eta);
530
531 if( c1 <=0){
532 c1 =1;
533 }
534
535 if( c1 == 1 ){
536 c1 = - 1;
537 }
538 if( c1 <0.6 || c1 >2.5){
539 c1 = -1;
540 }
541
542
543 CphiCorr[eta+85][phi] = c1;
544 }
545 }
546
547
548
549 }
550
551
552
553
554 void readInterCalibBarre_Phiv1(){
555
556 char *file_input = new char[500];
557 sprintf(file_input,"phisym_EcalIntercalib-160M-corrected.txt");
558
559
560 cout<<"READING fro inter-calib file "<<file_input<<endl;
561 ifstream inputcc(file_input,ios::in);
562 if (inputcc.fail()){
563 cout<<"error open file barrel.. " << file_input<<endl;
564 exit(1);
565 }
566
567 int eta;
568 int phi;
569 float c = 1;
570 float c1 = 1;
571
572 int n = 0;
573
574 eta = -85;
575 phi = 0;
576
577 float mean = 0;
578
579 int cureta = -85;
580
581 int ngood = 0;
582
583
584
585 while (inputcc.good()){
586
587 inputcc >> c1;
588
589 if(c1==1) c1 = -1;
590
591
592 if( n% 360 ==0 && n>0 ) {
593 eta++;
594 if(eta==0) eta = 1;
595 }
596
597 if( eta == cureta){
598 phi ++;
599 }else{
600 cureta = eta;
601 phi = 1;
602
603 }
604
605
606 n++;
607
608
609 int iSM = (phi-1)/20+1;
610 if( eta<0) iSM += 18;
611
612
613 int ieta = eta;
614 int iphi = phi;
615
616
617 convxtalid(iphi,ieta);
618
619 ///txtout<<eta<<" "<<phi<<endl;
620
621 CphiCorr[ieta+85][iphi] = c1;
622
623
624 if(c1 >0 ){
625 mean += c1;
626 ngood ++;
627 }
628
629
630 if(n >= 360*170) break;
631
632 }
633
634
635 mean /= ngood;
636 cout<<n<<" phi constants read. "<< mean <<endl;
637
638
639 // for(int j=0; j< 170; j++){
640 // mean_ieta[j] /= nmean_ieta[j];
641 // txtout<<"mean_ieta_phisym "<< j<<" "<< mean_ieta[j] <<" "<< nmean_ieta[j] <<endl;
642 // }
643
644
645 // for(int j=0; j< 170; j++){
646 // for(int k=0; k< 360; k++){
647
648 // if( CphiCorr[j][k] >0 ){
649 // CphiCorr[j][k] /= mean_ieta[j] ;
650
651 // }
652 // Cphi[j][k] = CphiCorr[j][k] * iCalConst_data[j][k];
653 // }
654 // }
655
656 // NormSMScaleToUnit(CphiCorr);
657
658 }
659
660
661 ///read initi calibration constants for each crystal
662 void readInterCalibBarre_BSv1(){
663
664
665 char *file_input = new char[500];
666 sprintf(file_input,"IC_barrel_splash09_AllRuns.txt");
667
668
669 cout<<"READING fro inter-calib bsfile "<<file_input<<endl;
670 txtout<<"READING fro inter-calib bsfile "<<file_input<<endl;
671 ifstream inputcc(file_input,ios::in);
672 if (inputcc.fail()){
673 cout<<"error open file barrel.. " << file_input<<endl;
674 exit(1);
675 }
676
677 int eta;
678 int phi;
679 float c = 1;
680 float c1 = 1;
681
682 int n = 0;
683 int ngood = 0;
684
685
686 float mean = 0;
687
688
689
690
691 float mean_ieta[170] = {0};
692 int nmean_ieta[170] = {0};
693
694
695 while (inputcc.good()){
696
697 inputcc >> eta >> phi >>c1;
698
699 // cout<<phi<<" "<<eta<<" "<<c<<endl;
700
701 convxtalid(phi,eta);
702
703 eta += 85;
704
705 if( c1 == 1){
706 c1 = -1;
707 }
708
709 if( c1 <0.6 || c1 >2.5){
710 c1 = -1;
711 }
712
713
714 ///iCalConstPi0[eta][phi] = c1 ;
715
716 CBSCorr[eta][phi] = c1;
717
718 if(c1 >0 ){
719 mean += c1;
720 ngood ++;
721 }
722
723 if(c1 >0 ){
724 mean_ieta[eta] += c1;
725 nmean_ieta[eta] ++;
726 }
727
728 n++;
729
730 /// hh_imcEB->Fill(c);
731 if(n >= 360*170) break;
732
733 }
734
735
736
737 mean /= ngood;
738
739
740
741 cout<<n<<" bs constants read. "<< mean << endl;
742
743
744
745 // for(int j=0; j< 170; j++){
746 // for(int k=0; k< 360; k++){
747
748 // if( CBSCorr[j][k] >0 ){
749 // CBSCorr[j][k] /= mean_ieta[j];
750 // }
751
752 // CBS[j][k] = CBSCorr[j][k] * iCalConst_data[j][k];
753
754
755 // }
756 // }
757 // NormSMScaleToUnit(CBSCorr);
758
759
760
761 }
762
763
764
765
766
767 ///read initi calibration constants for each crystal
768 void readInterCalibBarre_GR09_V8(){
769
770
771
772 char *file_input = new char[500];
773 //file_input ="/data/yangyong/data/p
774 file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_R_V6A_barrel.txt";
775
776 ///sprintf(file_input,"/uscms_data/d2/yongy/calib/backup/interCalib_STARTUP3X_V8P_barrel_mc.txt"); ///now the same used in data
777
778 //sprintf(file_input,"/uscms/home/yongy/work/interCalib_STARTUP3X_V8P_barrel.txt");
779
780 // sprintf(file_input,"/uscms/home/yongy/work/interCalib_GR_R_36X_V11A_Barrel.txt");
781
782 cout<<"READING fro inter-calib file "<<file_input<<endl;
783 ifstream inputcc(file_input,ios::in);
784 if (inputcc.fail()){
785 cout<<"error open file barrel.. " << file_input<<endl;
786 exit(1);
787 }
788
789 int eta;
790 int phi;
791 float c = 1;
792 float c1 = 1;
793
794 int n = 0;
795
796 double mean = 0;
797
798
799 while (inputcc.good()){
800
801 inputcc >> eta >> phi >>c1;
802
803 // cout<<phi<<" "<<eta<<" "<<c<<endl;
804
805 convxtalid(phi,eta);
806
807 eta += 85;
808
809 ///iCalConst[eta][phi] = c;
810
811 //iCalConst_data[eta][phi] = c1 * 0.038012/ 0.038; ///Now this is exactly the same as the inter-calibration used in data.
812
813 interCalib_preCalib[eta][phi] = c1 ; /// now OK iN MC startupv8p.
814
815
816 n++;
817 // cout<<phi<<" "<<eta<<" "<<c<<endl;
818 mean += c1;
819
820 /// hh_imcEB->Fill(c);
821 if(n >= 360*170) break;
822
823 }
824
825
826 mean /= 61200;
827
828 cout<<n<<" constants read. "<< mean <<endl;
829
830
831 }
832
833
834
835
836
837
838
839
840 ///read initi calibration constants for each crystal
841 void readInterCalibEndcap_GR09_V8(){
842 char *file_input = new char[500];
843 //file_input ="/data/yangyong/data/p
844 ///file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_R_V6A_barrel.txt";
845
846 file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_H_V6OFF_endcap.txt";
847
848 cout<<"READING fro inter-calib file "<<file_input<<endl;
849 ifstream inputcc(file_input,ios::in);
850 if (inputcc.fail()){
851 cout<<"error open file barrel.. " << file_input<<endl;
852 exit(1);
853 }
854
855 int eta;
856 int phi;
857 int iz;
858 float c = 1;
859 int n = 0;
860
861
862 while (inputcc.good()){
863
864 inputcc >> iz >> eta >> phi >> c;
865
866 int izz = iz < 0? 0:1;
867
868 interCalibEndcap_preCalib[izz][eta][phi] = c;
869 validRecHitEndCap[izz][eta][phi] = 1;
870
871 n ++;
872
873 }
874
875 cout<<n<<" constants read. "<<endl;
876
877 }
878
879
880 void getInterCalibEndcapv1(const char *inputfile,float C[2][101][101]){
881
882
883
884 ifstream inputcc(inputfile);
885 if (inputcc.fail()) {
886 cout<<"error "<< inputfile<<" can not be opened."<<endl;
887 exit(1);
888 }
889 int eta;
890 int phi;
891 int iz;
892 float c = 1;
893 int n = 0;
894
895 double mean = 0;
896
897 while (inputcc.good()){
898
899 inputcc >> iz >> eta >> phi >> c;
900
901 int izz = iz < 0? 0:1;
902
903 C[izz][eta][phi] = c;
904 mean += c;
905
906 n ++;
907
908 if( n>= 14648){
909 break;
910 }
911
912 }
913
914 cout<<n<<" constants read. "<< mean /n <<endl;
915
916 }
917
918
919
920 void getInterCalibEndcap(const char *inputfile, float C[2][101][101]){
921
922
923 ifstream inputcc(inputfile);
924 if (inputcc.fail()) {
925 cout<<"error "<< inputfile<<" can not be opened."<<endl;
926 exit(1);
927 }
928
929 int ix,iy,iz;
930 float cc;
931
932 float ccMeaniring[50] ;
933 int nccMeaniring[50];
934
935
936 // for(int j=0; j<2;j++){
937 for(int k=0; k< kEndcEtaRings; k++){
938 ccMeaniring[k] = 0;
939 nccMeaniring[k] = 0;
940 }
941 // }
942
943 int nread = 0;
944 while( inputcc.good()){
945 inputcc >> iz >> ix >> iy >> cc;
946 C[iz][ix][iy] = cc;
947 int iring = iRingEndCap(2*iz-1,ix,iy); ///input -1/1,
948 if( cc > 0){
949 ccMeaniring[iring] += cc;
950 nccMeaniring[iring] += 1;
951 }
952
953 nread ++;
954 if( iz==1 && ix == 100 && iy == 60) break;
955 }
956
957 cout<<"nb of crystals endcap " << nread <<endl;
958
959
960 //for(int j=0; j<2;j++){
961 for(int k=0; k< kEndcEtaRings; k++){
962 ccMeaniring[k] /= nccMeaniring[k];
963 // }
964 }
965
966 double meanall = 0;
967 int nall = 0;
968 for(int iz=0; iz<2; iz++){
969 for(int j=0; j<101; j++){
970 for(int k=0; k< 101; k++){
971 if( validRecHitEndCap[iz][j][k] <1) continue;
972 int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
973 if( C[iz][j][k] >0) {
974 C[iz][j][k] /= ccMeaniring[iring];
975 meanall += C[iz][j][k];
976 nall ++;
977 }
978
979 }
980 }
981 }
982
983 meanall /= nall;
984 cout<<" corrfactoriZiXiY read " << meanall <<" "<< nall <<endl;
985
986
987 }
988
989
990
991 void NormCrystalDeadFlagEndcap_v1(float C[2][101][101],int ndeadflag[2][101][101]){
992
993 float mean_deadflag[20] = {0};
994 float nmean_deadflag[20] = {0};
995
996 int ndeadCrystals = 0;
997
998
999 for(int iz =0; iz < 2; iz++){
1000 for(int j=0; j< 101; j++){
1001 for(int k=0; k<101; k++){
1002
1003 if( validRecHitEndCap[iz][j][k] ==0) continue;
1004 if( C[iz][j][k] <0){
1005 continue;
1006 }
1007 int ndead = ndeadflag[iz][j][k];
1008 if( ndead >=0){
1009 mean_deadflag[ndead] += C[iz][j][k];
1010 nmean_deadflag[ndead] ++;
1011 }
1012
1013 }
1014 }
1015 }
1016 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
1017 for(int j=0; j<20; j++){
1018 if( nmean_deadflag[j] >0){
1019 mean_deadflag[j] /= nmean_deadflag[j];
1020 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
1021 }
1022 }
1023
1024
1025 for(int iz =0; iz < 2; iz++){
1026 for(int j=0; j< 101; j++){
1027 for(int k=0; k<101; k++){
1028
1029 if( validRecHitEndCap[iz][j][k] ==0) continue;
1030 if( C[iz][j][k] <0){
1031 continue;
1032 }
1033
1034 int ndead = ndeadflag[iz][j][k];
1035 if( ndead >=0){
1036 C[iz][j][k] /= mean_deadflag[ndead];
1037 }
1038
1039 }
1040 }
1041 }
1042
1043
1044
1045 }
1046
1047
1048
1049
1050
1051
1052 void NormCrystalDeadFlagEndcap(float C[2][101][101]){
1053
1054 float mean_deadflag[20] = {0};
1055 float nmean_deadflag[20] = {0};
1056
1057 int ndeadCrystals = 0;
1058
1059
1060 for(int iz =0; iz < 2; iz++){
1061 for(int j=0; j< 101; j++){
1062 for(int k=0; k<101; k++){
1063
1064 if( validRecHitEndCap[iz][j][k] ==0) continue;
1065 if( C[iz][j][k] <0){
1066 continue;
1067 }
1068 int ndeadflag = ndeadflag_endcap[iz][j][k];
1069 if( ndeadflag >=0){
1070 mean_deadflag[ndeadflag] += C[iz][j][k];
1071 nmean_deadflag[ndeadflag] ++;
1072 }
1073
1074 }
1075 }
1076 }
1077 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
1078 for(int j=0; j<20; j++){
1079 if( nmean_deadflag[j] >0){
1080 mean_deadflag[j] /= nmean_deadflag[j];
1081 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
1082 }
1083 }
1084
1085
1086 for(int iz =0; iz < 2; iz++){
1087 for(int j=0; j< 101; j++){
1088 for(int k=0; k<101; k++){
1089
1090 if( validRecHitEndCap[iz][j][k] ==0) continue;
1091 if( C[iz][j][k] <0){
1092 continue;
1093 }
1094
1095 int ndeadflag = ndeadflag_endcap[iz][j][k];
1096 if( ndeadflag >=0){
1097 C[iz][j][k] /= mean_deadflag[ndeadflag];
1098 }
1099
1100 }
1101 }
1102 }
1103
1104
1105
1106 }
1107
1108
1109
1110 void scaleMeanToUnitEndcap(float C[2][101][101]){
1111
1112
1113 double mean = 0;
1114 int n = 0;
1115
1116
1117 for(int iz =0; iz < 2; iz++){
1118 for(int j=0; j< 101; j++){
1119 for(int k=0; k<101; k++){
1120 if( validRecHitEndCap[iz][j][k] ==0) continue;
1121 if( C[iz][j][k] <0){
1122 continue;
1123 }
1124 mean += C[iz][j][k];
1125 n ++;
1126 }
1127 }
1128 }
1129
1130 mean /= n;
1131 cout<<"mean " << mean<<" "<<n<<endl;
1132
1133
1134 for(int iz =0; iz < 2; iz++){
1135 for(int j=0; j< 101; j++){
1136 for(int k=0; k<101; k++){
1137 if( validRecHitEndCap[iz][j][k] ==0) continue;
1138 if( C[iz][j][k] <0){
1139 continue;
1140 }
1141
1142 C[iz][j][k] /= mean;
1143 }
1144 }
1145 }
1146
1147 }