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
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 yangyong 1.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     }