ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/yangyong/combineICEE/getCalibConstants.cc
Revision: 1.2
Committed: Thu Aug 30 11:03:29 2012 UTC (12 years, 8 months ago) by yangyong
Content type: text/plain
Branch: MAIN
CVS Tags: V2012comb, HEAD
Changes since 1.1: +0 -393 lines
Error occurred while calculating annotation data.
Log Message:
*** empty log message ***

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 NormIetaToUnitTestBeamSMOnly(float C[170][360]){
34
35 float mean_ieta[170] = {0};
36 int nmean_ieta[170] = {0};
37
38 float mean_ietaTT[170] = {0};
39 int nmean_ietaTT[170] = {0};
40
41
42 for(int j=0; j< 170; j++){
43
44 int ieta = j-85;
45 if( ieta >=0) ieta += 1;
46
47 for(int k=0; k< 360; k++){
48
49 int iphi = k;
50 if( k==0) iphi = 360;
51
52 int iSM = (iphi-1)/20+1;
53 if( ieta<0) iSM += 18;
54 int smTB = isTestBeamSM(iSM);
55 if( smTB ==0) continue;
56
57 //if( C[j][k] >0 ){
58 if( C[j][k] >0.5 && C[j][k] < 1.5 ){
59 mean_ieta[j] += C[j][k];
60 nmean_ieta[j] ++;
61
62 mean_ietaTT[j/5] += C[j][k];
63 nmean_ietaTT[j/5] ++;
64
65 }
66
67 }
68 }
69
70
71 for(int j=0; j< 170; j++){
72 mean_ieta[j] /= nmean_ieta[j];
73 //cout<<" mean_ieta[" << j<<"] "<< mean_ieta[j] <<endl;
74 }
75
76 // for(int j=0; j< 34; j++){
77 // mean_ietaTT[j] /= nmean_ietaTT[j];
78 // cout<<" mean_ietaTT[" << j<<"] "<< mean_ietaTT[j] <<endl;
79 // }
80
81 for(int j=0; j< 170; j++){
82 for(int k=0; k< 360; k++){
83
84 if( C[j][k] >0 ){
85 C[j][k] /= mean_ieta[j] ;
86 }
87 }
88 }
89
90 }
91
92
93
94 void NormIetaAbsToUnitTestBeamSMOnly(float C[170][360]){
95
96
97 float mean_ieta[170] = {0};
98 int nmean_ieta[170] = {0};
99
100 float mean_ietaTT[170] = {0};
101 int nmean_ietaTT[170] = {0};
102
103
104 for(int j=0; j< 170; j++){
105
106 int ieta = j-85;
107 if( ieta >=0) ieta += 1;
108
109 for(int k=0; k< 360; k++){
110
111 int iphi = k;
112 if( k==0) iphi = 360;
113
114 int iSM = (iphi-1)/20+1;
115 if( ieta<0) iSM += 18;
116 int smTB = isTestBeamSM(iSM);
117 if( smTB ==0) continue;
118
119 int ietaAbs = abs(ieta);
120
121 //if( C[j][k] >0 ){
122 if( C[j][k] >0.5 && C[j][k] < 1.5 ){
123 mean_ieta[ietaAbs-1] += C[j][k];
124 nmean_ieta[ietaAbs-1] ++;
125
126 mean_ietaTT[j/5] += C[j][k];
127 nmean_ietaTT[j/5] ++;
128
129 }
130 }
131 }
132
133
134 for(int j=0; j< 85; j++){
135 mean_ieta[j] /= nmean_ieta[j];
136 //cout<<" mean_ietaAbs[" << j<<"] "<< mean_ieta[j] <<endl;
137 }
138
139 // for(int j=0; j< 34; j++){
140 // mean_ietaTT[j] /= nmean_ietaTT[j];
141 // cout<<" mean_ietaTT[" << j<<"] "<< mean_ietaTT[j] <<endl;
142 // }
143
144 for(int j=0; j< 170; j++){
145
146 int ieta = j-85;
147 if( ieta >=0) ieta += 1;
148
149 for(int k=0; k< 360; k++){
150
151 int ietaAbs = abs(ieta);
152
153 if( C[j][k] >0 ){
154
155 C[j][k] /= mean_ieta[ietaAbs-1] ;
156
157 }
158 }
159 }
160
161 }
162
163
164
165 void NormCrystalDeadFlag_v1(float C[170][360], int ndeadflagietaiphi[170][360]){
166
167 float mean_deadflag[20] = {0};
168 float nmean_deadflag[20] = {0};
169
170 int ndeadCrystals = 0;
171 for(int j=0; j< 170; j++){
172 for(int k=0; k< 360; k++){
173 if( C[j][k] <0){
174 ndeadCrystals++;
175 continue;
176 }
177 int ndeadflag = ndeadflagietaiphi[j][k];
178 if( ndeadflag >=0){
179 mean_deadflag[ndeadflag] += C[j][k];
180 nmean_deadflag[ndeadflag] ++;
181 }
182 }
183 }
184 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
185 for(int j=0; j<20; j++){
186 if( nmean_deadflag[j] >0){
187 mean_deadflag[j] /= nmean_deadflag[j];
188 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
189 }
190 }
191 for(int j=0; j< 170; j++){
192 for(int k=0; k< 360; k++){
193 if( C[j][k] <0) continue;
194 int ndeadflag = ndeadflagietaiphi[j][k];
195 if( ndeadflag >=0){
196 C[j][k] /= mean_deadflag[ndeadflag];
197 }
198 }
199 }
200
201 }
202
203 void NormCrystalDeadFlag(float C[170][360]){
204
205 float mean_deadflag[20] = {0};
206 float nmean_deadflag[20] = {0};
207
208 int ndeadCrystals = 0;
209 for(int j=0; j< 170; j++){
210 for(int k=0; k< 360; k++){
211 if( C[j][k] <0){
212 ndeadCrystals++;
213 continue;
214 }
215 int ndeadflag = ndeadflag_ietaiphi[j][k];
216 if( ndeadflag >=0){
217 mean_deadflag[ndeadflag] += C[j][k];
218 nmean_deadflag[ndeadflag] ++;
219 }
220 }
221 }
222 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
223 for(int j=0; j<20; j++){
224 if( nmean_deadflag[j] >0){
225 mean_deadflag[j] /= nmean_deadflag[j];
226 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
227 }
228 }
229 for(int j=0; j< 170; j++){
230 for(int k=0; k< 360; k++){
231 if( C[j][k] <0) continue;
232 int ndeadflag = ndeadflag_ietaiphi[j][k];
233 if( ndeadflag >=0){
234 C[j][k] /= mean_deadflag[ndeadflag];
235 }
236 }
237 }
238
239 }
240
241
242
243 void SetSMScale(float C[170][360],TH1F *hhtmp){
244 float mean_ism[36] = {0};
245 float nmean_ism[36]= {0};
246
247 for(int j=0; j< 170; j++){
248 for(int k=0; k< 360; k++){
249
250 int ieta = j-85;
251 if( j >=0) ieta += 1;
252 int iphi = k;
253 if( k==0) iphi = 360;
254 int iSM = (iphi-1)/20+1;
255 if( ieta<0) iSM += 18;
256
257 if( C[j][k] >0 ){
258 mean_ism[iSM-1 ] += C[j][k];
259 nmean_ism[iSM-1] ++;
260 }
261 }
262 }
263
264 for(int j=0; j<36; j++){
265 mean_ism[j] /= nmean_ism[j];
266 /// cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl;
267 }
268
269
270 for(int j=0; j< 170; j++){
271 for(int k=0; k< 360; k++){
272
273 int ieta = j-85;
274 if( j >=0) ieta += 1;
275 int iphi = k;
276 if( k==0) iphi = 360;
277 int iSM = (iphi-1)/20+1;
278 if( ieta<0) iSM += 18;
279
280 if( C[j][k] >0 ){
281 C[j][k] *= hhtmp->GetBinContent(iSM);
282 }
283 }
284 }
285
286
287 }
288
289
290
291
292
293
294
295 void NormSMScaleToUnit(float C[170][360]){
296 float mean_ism[36] = {0};
297 float nmean_ism[36]= {0};
298
299 for(int j=0; j< 170; j++){
300 for(int k=0; k< 360; k++){
301
302 int ieta = j-85;
303 if( j >=0) ieta += 1;
304 int iphi = k;
305 if( k==0) iphi = 360;
306 int iSM = (iphi-1)/20+1;
307 if( ieta<0) iSM += 18;
308
309
310 if( C[j][k] >0 ){
311 mean_ism[iSM-1 ] += C[j][k];
312 nmean_ism[iSM-1] ++;
313 }
314 }
315 }
316
317 for(int j=0; j<36; j++){
318 mean_ism[j] /= nmean_ism[j];
319 //cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl;
320 }
321
322 for(int j=0; j< 170; j++){
323 for(int k=0; k< 360; k++){
324
325 int ieta = j-85;
326 if( j >=0) ieta += 1;
327 int iphi = k;
328 if( k==0) iphi = 360;
329 int iSM = (iphi-1)/20+1;
330 if( ieta<0) iSM += 18;
331
332 if( C[j][k] >0 ){
333 C[j][k] /= mean_ism[iSM-1];
334 }
335
336 }
337 }
338
339 }
340
341
342
343 void readInterCalibConstEBSimple(const char *input,float C[170][360]){
344
345 ifstream txtin1(input,ios::in);
346 if (txtin1.fail()){
347 cout<<"error open file barrel.. " << input<<endl;
348 exit(1);
349 }
350
351
352 int eta;
353 int phi;
354
355 float cc;
356
357 int n =0;
358 while(txtin1.good()){
359
360
361 txtin1 >> eta >> phi >> cc;
362
363 int iphi = phi;
364 int ieta = eta;
365
366 // cout<< eta <<" "<< phi <<" "<< cc <<endl;
367
368 convxtalid(phi,eta);
369 eta += 85;
370 C[eta][phi] = cc;
371
372 n++;
373 if( n>= 61200) break;
374
375 }
376
377 }
378
379
380
381
382
383 void readInterCalibConstEBSimplev1(const char *input,float C[170][360]){
384
385 ifstream txtin1(input,ios::in);
386 if (txtin1.fail()){
387 cout<<"error open file barrel.. " << input<<endl;
388 exit(1);
389 }
390
391
392 int eta;
393 int phi;
394
395 float cc;
396 double mean = 0;
397 int ngood = 0;
398 int n =0;
399 while(txtin1.good()){
400
401 txtin1 >> eta >> phi >> cc;
402
403 int iphi = phi;
404 int ieta = eta;
405
406 if( !(eta>=0 && eta <= 169 && phi >=0 && phi<=359)) {
407 cout<<"wrong eta/phi " << eta<<" "<<phi <<endl;
408 exit(1);
409 }
410
411 C[eta][phi] = cc;
412 if( cc >0){
413 mean += cc;
414 ngood ++;
415 }
416 n++;
417 if( n>= 61200) break;
418 }
419 mean /= ngood;
420 cout<<"read ngood "<< ngood <<" mean: "<< mean<<endl;
421
422 for(int j=0; j<170; j++){
423 for(int k=0; k<360;k++){
424 if(C[j][k] >0) {
425 C[j][k] /= mean;
426 }
427 }
428 }
429
430 }
431
432
433
434
435
436
437
438
439 ///read initi calibration constants for each crystal
440 void readInterCalibEndcap_GR09_V8(){
441 char *file_input = new char[500];
442 file_input = "interCalib_GR09_H_V6OFF_endcap.txt";
443
444 cout<<"READING fro inter-calib file "<<file_input<<endl;
445 ifstream inputcc(file_input,ios::in);
446 if (inputcc.fail()){
447 cout<<"error open file barrel.. " << file_input<<endl;
448 exit(1);
449 }
450
451 int eta;
452 int phi;
453 int iz;
454 float c = 1;
455 int n = 0;
456
457
458 while (inputcc.good()){
459
460 inputcc >> iz >> eta >> phi >> c;
461
462 int izz = iz < 0? 0:1;
463
464 interCalibEndcap_preCalib[izz][eta][phi] = c;
465 validRecHitEndCap[izz][eta][phi] = 1;
466
467 n ++;
468
469 }
470
471 cout<<n<<" constants read. "<<endl;
472
473 }
474
475
476 void getInterCalibEndcapv1(const char *inputfile,float C[2][101][101]){
477
478
479
480 ifstream inputcc(inputfile);
481 if (inputcc.fail()) {
482 cout<<"error "<< inputfile<<" can not be opened."<<endl;
483 exit(1);
484 }
485 int eta;
486 int phi;
487 int iz;
488 float c = 1;
489 int n = 0;
490
491 double mean = 0;
492
493 while (inputcc.good()){
494
495 inputcc >> iz >> eta >> phi >> c;
496
497 int izz = iz < 0? 0:1;
498
499 C[izz][eta][phi] = c;
500 mean += c;
501
502 n ++;
503
504 if( n>= 14648){
505 break;
506 }
507
508 }
509
510 cout<<n<<" constants read. "<< mean /n <<endl;
511
512 }
513
514
515
516 void getInterCalibEndcap(const char *inputfile, float C[2][101][101]){
517
518
519 ifstream inputcc(inputfile);
520 if (inputcc.fail()) {
521 cout<<"error "<< inputfile<<" can not be opened."<<endl;
522 exit(1);
523 }
524
525 int ix,iy,iz;
526 float cc;
527
528 float ccMeaniring[50] ;
529 int nccMeaniring[50];
530
531
532 // for(int j=0; j<2;j++){
533 for(int k=0; k< kEndcEtaRings; k++){
534 ccMeaniring[k] = 0;
535 nccMeaniring[k] = 0;
536 }
537 // }
538
539 int nread = 0;
540 while( inputcc.good()){
541 inputcc >> iz >> ix >> iy >> cc;
542 C[iz][ix][iy] = cc;
543 int iring = iRingEndCap(2*iz-1,ix,iy); ///input -1/1,
544 if( cc > 0){
545 ccMeaniring[iring] += cc;
546 nccMeaniring[iring] += 1;
547 }
548
549 nread ++;
550 if( iz==1 && ix == 100 && iy == 60) break;
551 }
552
553 cout<<"nb of crystals endcap " << nread <<endl;
554
555
556 //for(int j=0; j<2;j++){
557 for(int k=0; k< kEndcEtaRings; k++){
558 ccMeaniring[k] /= nccMeaniring[k];
559 // }
560 }
561
562 double meanall = 0;
563 int nall = 0;
564 for(int iz=0; iz<2; iz++){
565 for(int j=0; j<101; j++){
566 for(int k=0; k< 101; k++){
567 if( validRecHitEndCap[iz][j][k] <1) continue;
568 int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
569 if( C[iz][j][k] >0) {
570 C[iz][j][k] /= ccMeaniring[iring];
571 meanall += C[iz][j][k];
572 nall ++;
573 }
574
575 }
576 }
577 }
578
579 meanall /= nall;
580 cout<<" corrfactoriZiXiY read " << meanall <<" "<< nall <<endl;
581
582
583 }
584
585
586
587 void NormCrystalDeadFlagEndcap_v1(float C[2][101][101],int ndeadflag[2][101][101]){
588
589 float mean_deadflag[20] = {0};
590 int nmean_deadflag[20] = {0};
591
592 int ndeadCrystals = 0;
593
594
595 for(int iz =0; iz < 2; iz++){
596 for(int j=0; j< 101; j++){
597 for(int k=0; k<101; k++){
598
599 if( validRecHitEndCap[iz][j][k] ==0) continue;
600 if( C[iz][j][k] <0){
601 continue;
602 }
603 int ndead = ndeadflag[iz][j][k];
604 if( ndead >=0 && C[iz][j][k] >0.7 && C[iz][j][k] <1.3 ){
605 mean_deadflag[ndead] += C[iz][j][k];
606 nmean_deadflag[ndead] ++;
607 }
608
609 }
610 }
611 }
612 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
613 for(int j=0; j<20; j++){
614 if( nmean_deadflag[j] >0){
615 mean_deadflag[j] /= nmean_deadflag[j];
616 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
617 }
618 }
619
620
621 for(int iz =0; iz < 2; iz++){
622 for(int j=0; j< 101; j++){
623 for(int k=0; k<101; k++){
624
625 if( validRecHitEndCap[iz][j][k] ==0) continue;
626 if( C[iz][j][k] <0){
627 continue;
628 }
629
630 int ndead = ndeadflag[iz][j][k];
631 if( ndead >=0){
632 C[iz][j][k] /= mean_deadflag[ndead];
633 }
634
635 }
636 }
637 }
638
639
640
641 }
642
643
644
645
646
647
648 void NormCrystalDeadFlagEndcap(float C[2][101][101]){
649
650 float mean_deadflag[20] = {0};
651 float nmean_deadflag[20] = {0};
652
653 int ndeadCrystals = 0;
654
655
656 for(int iz =0; iz < 2; iz++){
657 for(int j=0; j< 101; j++){
658 for(int k=0; k<101; k++){
659
660 if( validRecHitEndCap[iz][j][k] ==0) continue;
661 if( C[iz][j][k] <0){
662 continue;
663 }
664 int ndeadflag = ndeadflag_endcap[iz][j][k];
665 if( ndeadflag >=0){
666 mean_deadflag[ndeadflag] += C[iz][j][k];
667 nmean_deadflag[ndeadflag] ++;
668 }
669
670 }
671 }
672 }
673 cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
674 for(int j=0; j<20; j++){
675 if( nmean_deadflag[j] >0){
676 mean_deadflag[j] /= nmean_deadflag[j];
677 cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
678 }
679 }
680
681
682 for(int iz =0; iz < 2; iz++){
683 for(int j=0; j< 101; j++){
684 for(int k=0; k<101; k++){
685
686 if( validRecHitEndCap[iz][j][k] ==0) continue;
687 if( C[iz][j][k] <0){
688 continue;
689 }
690
691 int ndeadflag = ndeadflag_endcap[iz][j][k];
692 if( ndeadflag >=0){
693 C[iz][j][k] /= mean_deadflag[ndeadflag];
694 }
695
696 }
697 }
698 }
699
700
701
702 }
703
704
705
706 void scaleMeanToUnitEndcap(float C[2][101][101]){
707
708
709 double mean = 0;
710 int n = 0;
711
712
713 for(int iz =0; iz < 2; iz++){
714 for(int j=0; j< 101; j++){
715 for(int k=0; k<101; k++){
716 if( validRecHitEndCap[iz][j][k] ==0) continue;
717 if( C[iz][j][k] <0){
718 continue;
719 }
720 mean += C[iz][j][k];
721 n ++;
722 }
723 }
724 }
725
726 mean /= n;
727 cout<<"mean " << mean<<" "<<n<<endl;
728
729
730 for(int iz =0; iz < 2; iz++){
731 for(int j=0; j< 101; j++){
732 for(int k=0; k<101; k++){
733 if( validRecHitEndCap[iz][j][k] ==0) continue;
734 if( C[iz][j][k] <0){
735 continue;
736 }
737
738 C[iz][j][k] /= mean;
739 }
740 }
741 }
742
743 }