yangyong/combineICEE/combineCalibConstantEndcapv1.C

#include "rootheader.h"
#include "testSelection_cluster.h"

#include "common_functions.cc"
#include "foldEndcap.cc"

TRandom3 *rgen_;

#include "usefullcode.cc"


float ccalibpretag[2][101][101];

float interCalib_preCalib[170][360];
float interCalibEndcap_preCalib[2][101][101];


bool isTestBeamSM(int iSM){
    

  //if( iSM == 1 || iSM == 2  || iSM == 11 || iSM == 15  || iSM == 21 || iSM == 24) return true; 
  
  if( iSM == 1  || iSM == 11 || iSM == 15  || iSM == 24) return true; 
  
  
  else return false; 
    
}


#include "effSigma.C"

#include "getCrystaldeadflag.cc"

#include "getCalibConstants.cc"

#include "gausfit.cc"


float fitWind = 3; 

void fitHistogram(TH1F *h1, double res[]){
  
  float mean = h1->GetMean();
  float meanErr = h1->GetMeanError();
  float rmsEff = 100 * h1->GetRMS();
  float rmsEffErr = 100 * h1->GetRMSError();
  float rmsGaus = 100 * h1->GetRMS();
  float rmsGausErr = 100 * h1->GetRMSError();
  if( h1->Integral()>50){
    double resd[10];
    effSigma(h1,resd);
    rmsEff = resd[2]*100;
    float resf[20];
    
    fitgauswind2(h1,fitWind,fitWind,resf);
    //fitgauswindRefit(h1,resf);
    rmsGaus = resf[6]*100;
    rmsGausErr = resf[3] * 100;
    rmsEffErr = resf[3] * 100;
    mean = resf[0];
    meanErr = resf[1];
  }
  res[0] = mean; 
  res[1] = meanErr; 
  res[2] = rmsEff;
  res[3] = rmsEffErr; 
  res[4] = rmsGaus; 
  res[5] = rmsGausErr; 
}


float ic[50][2][101][101]; //at most 10 sets

int ndeadflag_ic[50][2][101][101];

float icwt_period[25][2][101][101];  ///pi0&eta combined of each period


void combinCalibConstantEndcapv1(){


  cout.setf(ios::fixed,ios::floatfield);
  cout.precision(8);
  

  for(int j=0; j<2; j++){
    for(int x =0; x<101; x++){
      for(int y=0; y<101; y++){
        validRecHitEndCap[j][x][y] = 0;
      }
    }
  }
  
  
  readInterCalibEndcap_GR09_V8();

  getInterCalibEndcapv1("interCalibEE_GR_P_V39.txt",ccalibpretag);
  
  
  get_xyzEBrechits();
  setEtaRingBoundaryEndcap();
  
  map<int,string> icFiles; 
  
  
  icFiles[0] = "calibres/deriveCalibConst.dflag2.pe1.step2.iter50.txt";
  icFiles[1] = "calibres/deriveCalibConst.dflag3.pe2.step2.iter50.txt";
  
  int nIC = 1; 
  
  getCrystaldeadflagEndcap_v1("deadflag/crystal_deadflag_ee_dflag2.txt",ndeadflag_ic[0]);
  getCrystaldeadflagEndcap_v1("deadflag/crystal_deadflag_ee_dflag3.txt",ndeadflag_ic[1]);
  
  
  int nConstantSet = int(icFiles.size());
  
  if( nConstantSet >50){
    cout<<"at most 50 IC " << nConstantSet <<endl; 
    return; 
  }
  
  
  float cc[2][101][101];
    
  TFile *fnew = new TFile("combineCalibConstantEndcapv1.root","recreate");
  
    
  TH1F *hh_cc_ietaring[50][3][40]; //file,ee-/ee+/both, ring, 
  TH1F *hh_res_cc_ietaring[50][3][3];
  
  //wted average cc
  TH1F *hh_ccwtavg_ietaring[3][40];
  TH1F *hh_res_ccwtavg_ietaring[3][3];

  TH1F *hh_ccwtavg_ietaring_period[25][3][40];
  TH1F *hh_res_ccwtavg_ietaring_period[25][3][3]; //of each period
  
  for(int n=0; n< nIC; n++){
    for(int j=0;j<3; j++){
      for(int k=0;k<40;k++){
        TString filename = TString (Form("hh_ccwtavg_ietaring_period_%d_%d_%d",n,j,k));
        hh_ccwtavg_ietaring_period[n][j][k] = new TH1F(filename,filename,200,0,4);
      }
    }
    for(int j=0;j<3; j++){
      for(int k=0;k<3;k++){
        TString filename = TString (Form("hh_res_ccwtavg_ietaring_period_%d_%d_%d",n,j,k));
        hh_res_ccwtavg_ietaring_period[n][j][k] = new TH1F(filename,filename,40,0,40);
      }
    }
  }
      
  
  TH2F *hh2_cc_xy[51][2];
  for(int j=0; j< nConstantSet+1; j++){
    for(int k=0;k<2; k++){
      TString filename = TString (Form("hh2_cc_xy_%d_%d",j,k));
      hh2_cc_xy[j][k] = new TH2F(filename,filename,100,1,101,100,1,101);
    }
  }
  
  TH2F *hh_corr_cc[50][10];
  for(int n=0; n<nConstantSet; n++){
    for(int k=n+1; k<nConstantSet; k++){
      TString filename = TString (Form("hh_corr_cc_%dand%d",n,k));
      hh_corr_cc[n][k] = new TH2F(filename,filename,400,0,2,400,0,2);
    }
  }
  
  
  for(int n=0; n<nConstantSet; n++){
    for(int j=0;j<3; j++){
      for(int k=0;k<40;k++){
        TString filename = TString (Form("hh_cc_ietaring_%d_%d_%d",n,j,k));
        hh_cc_ietaring[n][j][k] = new TH1F(filename,filename,200,0,4);
      }
    }
    for(int j=0;j<3; j++){
      for(int k=0;k<3;k++){
        TString filename = TString (Form("hh_res_cc_ietaring_%d_%d_%d",n,j,k));
        hh_res_cc_ietaring[n][j][k] = new TH1F(filename,filename,40,0,40);
      }
    }
  }
    
  for(int j=0;j<3; j++){
    for(int k=0;k<40;k++){
      TString filename = TString (Form("hh_ccwtavg_ietaring_%d_%d",j,k));
      hh_ccwtavg_ietaring[j][k] = new TH1F(filename,filename,200,0,4);
    }
  }
  for(int j=0;j<3; j++){
    for(int k=0;k<3;k++){
      TString filename = TString (Form("hh_res_ccwtavg_ietaring_%d_%d",j,k));
      hh_res_ccwtavg_ietaring[j][k] = new TH1F(filename,filename,40,0,40);
    }
  }

  
  //(c1 - c2) / ( average) 
  TH1F *hh_diff_cc_ietaring[50][50][3][40];
  TH1F *hh_res_diff_cc_ietaring[50][50][3][3];
  TH2F *hh2_diff_cc[50][50][2];
  
  
  for(int n=0; n<nConstantSet; n++){
    for(int k=n+1; k<nConstantSet; k++){
      for(int j=0;j<2;j++){
        TString filename = TString (Form("hh2_diff_cc_%dand%d_%d",n,k,j));
        hh2_diff_cc[n][k][j] = new TH2F(filename,filename,100,1,101,100,1,101);
      }
    }
  }
  
  
  for(int n=0; n<nConstantSet; n++){
    for(int k=n+1; k<nConstantSet; k++){

      for(int m=0; m<3; m++){
        for(int j=0;j<40;j++){
          TString filename = TString (Form("hh_diff_cc_ietaring_%dand%d_%d_%d",n,k,m,j));
          hh_diff_cc_ietaring[n][k][m][j] = new TH1F(filename,filename,100,-0.5,0.5);
        }
      }
      
      for(int m=0; m<3; m++){
        for(int j=0;j<3;j++){
          TString filename = TString (Form("hh_res_diff_cc_ietaring_%dand%d_%d_%d",n,k,m,j));
          hh_res_diff_cc_ietaring[n][k][m][j] = new TH1F(filename,filename,40,0,40);
        }
      }
    }
  }
  
  
  TH1F *hh_c_deadflag[51][30]; 
  TH1F *hh_c_deadflag_period[51][30]; 
  for(int j=0; j< nConstantSet+1; j++){
    for(int k=0; k<30; k++){
      TString histname = TString(Form("hh_c_deadflag_%d_%d",j,k));
      hh_c_deadflag[j][k] = new TH1F(histname,histname,500,0,2);

      histname = TString(Form("hh_c_deadflag_period_%d_%d",j,k));
      hh_c_deadflag_period[j][k] = new TH1F(histname,histname,500,0,2);
      
    }
  }
  
  
  ofstream txtoutcheck("combineCalibConstantEndcapv1.txt",ios::out);
    
  
  for(int n=0; n< nConstantSet; n++){
    string  filename = icFiles[n];


    getInterCalibEndcap(filename.c_str(),cc);


    NormCrystalDeadFlagEndcap_v1(cc,ndeadflag_ic[n]);
    
    
    scaleMeanToUnitEndcap(cc);
    for(int iz=0; iz<2; iz++){
      for(int j=0; j<101; j++){
        for(int k=0; k< 101; k++){
          ic[n][iz][j][k] = cc[iz][j][k];

          if( ndeadflag_ic[n][iz][j][k] < 0  && ic[n][iz][j][k] >0){
            cout<<"warning (can be ignored) dead crystal  " <<ic[n][iz][j][k] <<endl; 
            ic[n][iz][j][k] = -1; 
          }
          
        }
      }
    }
    
    
  }
  

  cout<<" nConstantSet " << nConstantSet <<endl; 
  

  for(int n=0; n< nConstantSet; n++){
    
     for(int iz=0; iz<2; iz++){
      for(int j=0; j<101; j++){
        for(int k=0; k< 101; k++){
          
          if( validRecHitEndCap[iz][j][k] <1) {
            continue; 
          }
          
          int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,  
          float c = ic[n][iz][j][k];
          if(c >0){
            hh_cc_ietaring[n][iz][iring]->Fill(c);
            hh_cc_ietaring[n][2][iring]->Fill(c);
            hh2_cc_xy[n][iz]->SetBinContent(j,k,c);
          }
        }
      }
     }
  }
  
  
  for(int n=0; n< nConstantSet; n++){
    for(int m=n+1; m < nConstantSet; m++){
      
      for(int iz=0; iz<2; iz++){
        for(int j=0; j<101; j++){
          for(int k=0; k< 101; k++){
            
            if( validRecHitEndCap[iz][j][k] <1) {
              hh2_diff_cc[n][m][iz]->SetBinContent(j,k,-1);
              continue; 
            }
            int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,  
            float c1 = ic[n][iz][j][k];
            float c2 = ic[m][iz][j][k];
            if(c1 >0 && c2 > 0 ){
              hh_corr_cc[n][m]->Fill(c1,c2);
              hh2_diff_cc[n][m][iz]->SetBinContent(j,k,c1-c2);
              hh_diff_cc_ietaring[n][m][iz][iring]->Fill( (c1 -c2)/(0.5*(c1+c2)));
              hh_diff_cc_ietaring[n][m][2][iring]->Fill( (c1 -c2)/(0.5*(c1+c2)));
              
            }else{
              hh2_diff_cc[n][m][iz]->SetBinContent(j,k,-1);
            }
          }
        }
      }
    }
  }
  
  
  cout<<"fit " <<endl; 
  
  
  double resfit[20];

  
  for(int n=0; n< nConstantSet; n++){
    for(int m=n+1; m < nConstantSet; m++){
      for(int k=0; k< kEndcEtaRings; k++){
        
        for(int iz=0; iz<3; iz++){
          fitHistogram(hh_diff_cc_ietaring[n][m][iz][k],resfit);
          
          for(int j=0; j<3; j++){
            hh_res_diff_cc_ietaring[n][m][iz][j]->SetBinContent(k+1,resfit[2*j]);
            hh_res_diff_cc_ietaring[n][m][iz][j]->SetBinError(k+1,resfit[2*j+1]);
          }
        }
      }
    }
  }
  
  
  for(int n=0; n< nConstantSet; n++){
    for(int j=0; j<3; j++){
      for(int k=0; k< kEndcEtaRings; k++){
        
        fitHistogram(hh_cc_ietaring[n][j][k],resfit);
        for(int m=0;m<3; m++){
          hh_res_cc_ietaring[n][j][m]->SetBinContent(k+1,resfit[2*m]);
          hh_res_cc_ietaring[n][j][m]->SetBinError(k+1,resfit[2*m+1]);
        }
        
      }
    }
  }
  
  
  float sigmaSys = 2; ///precision of precalib+LC 
  
  
  cout<<"combine " <<endl; 
  float icwt[2][101][101];
  
  
  float wtSumC_period[25] = {0};
  float wtSumS_period[25] = {0};
    

  for(int iz=0; iz<2; iz++){
    for(int j=0; j<101; j++){
      for(int k=0; k< 101; k++){
          
        if( validRecHitEndCap[iz][j][k] <1) {
          continue; 
        }
        int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,  
        
        float wtSumC = 0; 
        float wtSumS = 0; 
        for(int n=0; n< nIC; n++){ ///for each period
          wtSumC_period[n] = 0; 
          wtSumS_period[n] = 0; 
        }
        
        for(int n=0; n< nConstantSet; n++){
          float c = ic[n][iz][j][k]; 


          float sigma = hh_res_cc_ietaring[n][2][2]->GetBinContent(iring+1);
          sigma = sqrt( sigma * sigma + sigmaSys * sigmaSys);
          
          if( c > 0){
            float tmp1 = c/ ( sigma * sigma);
            float tmp2 = 1/(sigma * sigma);
            
            wtSumC += tmp1; 
            wtSumS += tmp2; 
            
            int nperiod = n% nIC; 
            wtSumC_period[nperiod] += tmp1; 
            wtSumS_period[nperiod] += tmp2; 
            
            int deadflag = ndeadflag_ic[n][iz][j][k];
            if(deadflag>=0){
              hh_c_deadflag[n][deadflag]->Fill( c);
            }

          }
        }
        
        if( wtSumC > 0){
          icwt[iz][j][k] = wtSumC / wtSumS; 


          hh_ccwtavg_ietaring[iz][iring] ->Fill( icwt[iz][j][k] );
          hh_ccwtavg_ietaring[2][iring] ->Fill( icwt[iz][j][k] );
          
        }else{
          icwt[iz][j][k] = -1;
        }
        for(int n=0; n< nIC; n++){ ///for each period
          if( wtSumC_period[n] > 0){
            icwt_period[n][iz][j][k] = wtSumC_period[n] / wtSumS_period[n]; 
            
            hh_ccwtavg_ietaring_period[n][iz][iring]->Fill( icwt_period[n][iz][j][k] );
            hh_ccwtavg_ietaring_period[n][2][iring]->Fill( icwt_period[n][iz][j][k] );
            
          }else{
            icwt_period[n][iz][j][k] = -1; 
          }
        }
        
      }

    }
  }
  
  
  for(int j=0; j<3; j++){
    for(int k=0; k< kEndcEtaRings; k++){
      fitHistogram(hh_ccwtavg_ietaring[j][k],resfit);
      for(int m=0;m<3; m++){
        hh_res_ccwtavg_ietaring[j][m]->SetBinContent(k+1,resfit[2*m]);
        hh_res_ccwtavg_ietaring[j][m]->SetBinError(k+1,resfit[2*m+1]);
      }
    }
  }
    
  for(int n=0; n< nIC; n++){ ///for each period
    for(int j=0; j<3; j++){
      for(int k=0; k< kEndcEtaRings; k++){
        fitHistogram(hh_ccwtavg_ietaring_period[n][j][k],resfit);
        for(int m=0;m<3; m++){
          hh_res_ccwtavg_ietaring_period[n][j][m]->SetBinContent(k+1,resfit[2*m]);
          hh_res_ccwtavg_ietaring_period[n][j][m]->SetBinError(k+1,resfit[2*m+1]);
        }
      }
    }
  }
  
  
  scaleMeanToUnitEndcap(icwt);
    
  ofstream txtout("interCalibConstants.combined.EcalEndcap.txt",ios::out);
  
  
  ofstream txtout_period[50]; //pi0eta combined for each period
  for(int j=0; j< nIC; j++){
    string filename = string(Form("interCalibConstants.combinedPi0EtaPeriod%d.EcalEndcap.txt",j));
    txtout_period[j].open(filename.c_str(),ios::out);
  }
  
  
  cout<<"print out final " <<endl; 
  for(int iz=0; iz<2; iz++){
    for(int j=0; j<101; j++){
      for(int k=0; k< 101; k++){
        if( validRecHitEndCap[iz][j][k] <1) continue;

        
        int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,  
        
        for(int n=0; n< nIC; n++){ ///for each period
          
          float sigmaC = hh_res_ccwtavg_ietaring_period[n][2][2]->GetBinContent(iring+1);
          
          float cErr = sqrt( sigmaC *sigmaC + sigmaSys * sigmaSys);
          cErr /=100;
          
          float c = icwt_period[n][iz][j][k];
          if( c >0){
            
            int deadflag1 = ndeadflag_ic[n%nIC][iz][j][k];
            int deadflag2 = ndeadflag_ic[nIC+n%nIC][iz][j][k];
            if(deadflag1<0 && deadflag2<0){
              cout<<"wrong deadflag ! " << n << " "<<iz <<" "<< j<<" "<<k <<endl; 
              return; 
            }
            if(deadflag1>0 ){
              hh_c_deadflag_period[n][deadflag1]->Fill( c);
              hh_c_deadflag_period[n][19]->Fill( c);
            }
            
          }
          
          if( c > 0){
            txtout_period[n]<<j<<" "<<k<<" "<< 2*iz-1<<" "<< c*ccalibpretag[iz][j][k]<<" "<< cErr * c*ccalibpretag[iz][j][k] <<endl; 
            
          }else{        
            txtout_period[n]<<j<<" "<<k<<" "<< 2*iz-1<<" "<<-1 <<" "<< 999 <<endl; 
          }
        }
        
        
        float c = icwt[iz][j][k];
        if( c>0){
          hh2_cc_xy[nConstantSet][iz]->SetBinContent(j,k,c);
          
          float sigmaC = hh_res_ccwtavg_ietaring[2][2]->GetBinContent(iring+1);
          
          float cErr = sqrt( sigmaC *sigmaC + sigmaSys * sigmaSys);
          cErr /=100;
          txtout<<j<<" "<<k<<" "<< 2*iz-1<<" "<< c * ccalibpretag[iz][j][k] <<" "<< cErr * c * ccalibpretag[iz][j][k] <<endl; 
        }else{
          txtout<<j<<" "<<k<<" "<< 2*iz-1<<" "<<-1 <<" "<< 999 <<endl; 
        }
      }
    }
  }
  
  
  fnew->Write();
  fnew->Close();
  
  
}
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:	+55 -82 lines
Log Message:	* empty log message *
#	Content
1	#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
12
13	float ccalibpretag[2][101][101];
14
15	float interCalib_preCalib[170][360];
16	float interCalibEndcap_preCalib[2][101][101];
17
18
19
20	bool isTestBeamSM(int iSM){
21
22
23
24	//if( iSM == 1 \|\| iSM == 2 \|\| iSM == 11 \|\| iSM == 15 \|\| iSM == 21 \|\| iSM == 24) return true;
25
26	if( iSM == 1 \|\| iSM == 11 \|\| iSM == 15 \|\| iSM == 24) return true;
27
28
29	else return false;
30
31	}
32
33
34	#include "effSigma.C"
35
36	#include "getCrystaldeadflag.cc"
37
38	#include "getCalibConstants.cc"
39
40	#include "gausfit.cc"
41
42
43	float fitWind = 3;
44
45	void fitHistogram(TH1F *h1, double res[]){
46
47	float mean = h1->GetMean();
48	float meanErr = h1->GetMeanError();
49	float rmsEff = 100 * h1->GetRMS();
50	float rmsEffErr = 100 * h1->GetRMSError();
51	float rmsGaus = 100 * h1->GetRMS();
52	float rmsGausErr = 100 * h1->GetRMSError();
53	if( h1->Integral()>50){
54	double resd[10];
55	effSigma(h1,resd);
56	rmsEff = resd[2]*100;
57	float resf[20];
58
59	fitgauswind2(h1,fitWind,fitWind,resf);
60	//fitgauswindRefit(h1,resf);
61	rmsGaus = resf[6]*100;
62	rmsGausErr = resf[3] * 100;
63	rmsEffErr = resf[3] * 100;
64	mean = resf[0];
65	meanErr = resf[1];
66	}
67	res[0] = mean;
68	res[1] = meanErr;
69	res[2] = rmsEff;
70	res[3] = rmsEffErr;
71	res[4] = rmsGaus;
72	res[5] = rmsGausErr;
73	}
74
75
76	float ic[50][2][101][101]; //at most 10 sets
77
78	int ndeadflag_ic[50][2][101][101];
79
80	float icwt_period[25][2][101][101]; ///pi0&eta combined of each period
81
82
83	void combinCalibConstantEndcapv1(){
84
85
86	cout.setf(ios::fixed,ios::floatfield);
87	cout.precision(8);
88
89
90	for(int j=0; j<2; j++){
91	for(int x =0; x<101; x++){
92	for(int y=0; y<101; y++){
93	validRecHitEndCap[j][x][y] = 0;
94	}
95	}
96	}
97
98
99	readInterCalibEndcap_GR09_V8();
100
101	getInterCalibEndcapv1("interCalibEE_GR_P_V39.txt",ccalibpretag);
102
103
104	get_xyzEBrechits();
105	setEtaRingBoundaryEndcap();
106
107	map<int,string> icFiles;
108
109
110	icFiles[0] = "calibres/deriveCalibConst.dflag2.pe1.step2.iter50.txt";
111	icFiles[1] = "calibres/deriveCalibConst.dflag3.pe2.step2.iter50.txt";
112
113	int nIC = 1;
114
115	getCrystaldeadflagEndcap_v1("deadflag/crystal_deadflag_ee_dflag2.txt",ndeadflag_ic[0]);
116	getCrystaldeadflagEndcap_v1("deadflag/crystal_deadflag_ee_dflag3.txt",ndeadflag_ic[1]);
117
118
119	int nConstantSet = int(icFiles.size());
120
121	if( nConstantSet >50){
122	cout<<"at most 50 IC " << nConstantSet <<endl;
123	return;
124	}
125
126
127
128	float cc[2][101][101];
129
130	TFile *fnew = new TFile("combineCalibConstantEndcapv1.root","recreate");
131
132
133
134	TH1F *hh_cc_ietaring[50][3][40]; //file,ee-/ee+/both, ring,
135	TH1F *hh_res_cc_ietaring[50][3][3];
136
137	//wted average cc
138	TH1F *hh_ccwtavg_ietaring[3][40];
139	TH1F *hh_res_ccwtavg_ietaring[3][3];
140
141	TH1F *hh_ccwtavg_ietaring_period[25][3][40];
142	TH1F *hh_res_ccwtavg_ietaring_period[25][3][3]; //of each period
143
144	for(int n=0; n< nIC; n++){
145	for(int j=0;j<3; j++){
146	for(int k=0;k<40;k++){
147	TString filename = TString (Form("hh_ccwtavg_ietaring_period_%d_%d_%d",n,j,k));
148	hh_ccwtavg_ietaring_period[n][j][k] = new TH1F(filename,filename,200,0,4);
149	}
150	}
151	for(int j=0;j<3; j++){
152	for(int k=0;k<3;k++){
153	TString filename = TString (Form("hh_res_ccwtavg_ietaring_period_%d_%d_%d",n,j,k));
154	hh_res_ccwtavg_ietaring_period[n][j][k] = new TH1F(filename,filename,40,0,40);
155	}
156	}
157	}
158
159
160	TH2F *hh2_cc_xy[51][2];
161	for(int j=0; j< nConstantSet+1; j++){
162	for(int k=0;k<2; k++){
163	TString filename = TString (Form("hh2_cc_xy_%d_%d",j,k));
164	hh2_cc_xy[j][k] = new TH2F(filename,filename,100,1,101,100,1,101);
165	}
166	}
167
168	TH2F *hh_corr_cc[50][10];
169	for(int n=0; n<nConstantSet; n++){
170	for(int k=n+1; k<nConstantSet; k++){
171	TString filename = TString (Form("hh_corr_cc_%dand%d",n,k));
172	hh_corr_cc[n][k] = new TH2F(filename,filename,400,0,2,400,0,2);
173	}
174	}
175
176
177	for(int n=0; n<nConstantSet; n++){
178	for(int j=0;j<3; j++){
179	for(int k=0;k<40;k++){
180	TString filename = TString (Form("hh_cc_ietaring_%d_%d_%d",n,j,k));
181	hh_cc_ietaring[n][j][k] = new TH1F(filename,filename,200,0,4);
182	}
183	}
184	for(int j=0;j<3; j++){
185	for(int k=0;k<3;k++){
186	TString filename = TString (Form("hh_res_cc_ietaring_%d_%d_%d",n,j,k));
187	hh_res_cc_ietaring[n][j][k] = new TH1F(filename,filename,40,0,40);
188	}
189	}
190	}
191
192	for(int j=0;j<3; j++){
193	for(int k=0;k<40;k++){
194	TString filename = TString (Form("hh_ccwtavg_ietaring_%d_%d",j,k));
195	hh_ccwtavg_ietaring[j][k] = new TH1F(filename,filename,200,0,4);
196	}
197	}
198	for(int j=0;j<3; j++){
199	for(int k=0;k<3;k++){
200	TString filename = TString (Form("hh_res_ccwtavg_ietaring_%d_%d",j,k));
201	hh_res_ccwtavg_ietaring[j][k] = new TH1F(filename,filename,40,0,40);
202	}
203	}
204
205
206	//(c1 - c2) / ( average)
207	TH1F *hh_diff_cc_ietaring[50][50][3][40];
208	TH1F *hh_res_diff_cc_ietaring[50][50][3][3];
209	TH2F *hh2_diff_cc[50][50][2];
210
211
212	for(int n=0; n<nConstantSet; n++){
213	for(int k=n+1; k<nConstantSet; k++){
214	for(int j=0;j<2;j++){
215	TString filename = TString (Form("hh2_diff_cc_%dand%d_%d",n,k,j));
216	hh2_diff_cc[n][k][j] = new TH2F(filename,filename,100,1,101,100,1,101);
217	}
218	}
219	}
220
221
222	for(int n=0; n<nConstantSet; n++){
223	for(int k=n+1; k<nConstantSet; k++){
224
225	for(int m=0; m<3; m++){
226	for(int j=0;j<40;j++){
227	TString filename = TString (Form("hh_diff_cc_ietaring_%dand%d_%d_%d",n,k,m,j));
228	hh_diff_cc_ietaring[n][k][m][j] = new TH1F(filename,filename,100,-0.5,0.5);
229	}
230	}
231
232	for(int m=0; m<3; m++){
233	for(int j=0;j<3;j++){
234	TString filename = TString (Form("hh_res_diff_cc_ietaring_%dand%d_%d_%d",n,k,m,j));
235	hh_res_diff_cc_ietaring[n][k][m][j] = new TH1F(filename,filename,40,0,40);
236	}
237	}
238	}
239	}
240
241
242
243	TH1F *hh_c_deadflag[51][30];
244	TH1F *hh_c_deadflag_period[51][30];
245	for(int j=0; j< nConstantSet+1; j++){
246	for(int k=0; k<30; k++){
247	TString histname = TString(Form("hh_c_deadflag_%d_%d",j,k));
248	hh_c_deadflag[j][k] = new TH1F(histname,histname,500,0,2);
249
250	histname = TString(Form("hh_c_deadflag_period_%d_%d",j,k));
251	hh_c_deadflag_period[j][k] = new TH1F(histname,histname,500,0,2);
252
253	}
254	}
255
256
257
258
259	ofstream txtoutcheck("combineCalibConstantEndcapv1.txt",ios::out);
260
261
262	for(int n=0; n< nConstantSet; n++){
263	string filename = icFiles[n];
264
265
266	getInterCalibEndcap(filename.c_str(),cc);
267
268
269
270	NormCrystalDeadFlagEndcap_v1(cc,ndeadflag_ic[n]);
271
272
273
274	scaleMeanToUnitEndcap(cc);
275	for(int iz=0; iz<2; iz++){
276	for(int j=0; j<101; j++){
277	for(int k=0; k< 101; k++){
278	ic[n][iz][j][k] = cc[iz][j][k];
279
280	if( ndeadflag_ic[n][iz][j][k] < 0 && ic[n][iz][j][k] >0){
281	cout<<"warning (can be ignored) dead crystal " <<ic[n][iz][j][k] <<endl;
282	ic[n][iz][j][k] = -1;
283	}
284
285	}
286	}
287	}
288
289
290
291	}
292
293
294	cout<<" nConstantSet " << nConstantSet <<endl;
295
296
297	for(int n=0; n< nConstantSet; n++){
298
299	for(int iz=0; iz<2; iz++){
300	for(int j=0; j<101; j++){
301	for(int k=0; k< 101; k++){
302
303	if( validRecHitEndCap[iz][j][k] <1) {
304	continue;
305	}
306
307	int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
308	float c = ic[n][iz][j][k];
309	if(c >0){
310	hh_cc_ietaring[n][iz][iring]->Fill(c);
311	hh_cc_ietaring[n][2][iring]->Fill(c);
312	hh2_cc_xy[n][iz]->SetBinContent(j,k,c);
313	}
314	}
315	}
316	}
317	}
318
319
320	for(int n=0; n< nConstantSet; n++){
321	for(int m=n+1; m < nConstantSet; m++){
322
323	for(int iz=0; iz<2; iz++){
324	for(int j=0; j<101; j++){
325	for(int k=0; k< 101; k++){
326
327	if( validRecHitEndCap[iz][j][k] <1) {
328	hh2_diff_cc[n][m][iz]->SetBinContent(j,k,-1);
329	continue;
330	}
331	int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
332	float c1 = ic[n][iz][j][k];
333	float c2 = ic[m][iz][j][k];
334	if(c1 >0 && c2 > 0 ){
335	hh_corr_cc[n][m]->Fill(c1,c2);
336	hh2_diff_cc[n][m][iz]->SetBinContent(j,k,c1-c2);
337	hh_diff_cc_ietaring[n][m][iz][iring]->Fill( (c1 -c2)/(0.5*(c1+c2)));
338	hh_diff_cc_ietaring[n][m][2][iring]->Fill( (c1 -c2)/(0.5*(c1+c2)));
339
340	}else{
341	hh2_diff_cc[n][m][iz]->SetBinContent(j,k,-1);
342	}
343	}
344	}
345	}
346	}
347	}
348
349
350
351	cout<<"fit " <<endl;
352
353
354
355	double resfit[20];
356
357
358
359
360	for(int n=0; n< nConstantSet; n++){
361	for(int m=n+1; m < nConstantSet; m++){
362	for(int k=0; k< kEndcEtaRings; k++){
363
364	for(int iz=0; iz<3; iz++){
365	fitHistogram(hh_diff_cc_ietaring[n][m][iz][k],resfit);
366
367	for(int j=0; j<3; j++){
368	hh_res_diff_cc_ietaring[n][m][iz][j]->SetBinContent(k+1,resfit[2*j]);
369	hh_res_diff_cc_ietaring[n][m][iz][j]->SetBinError(k+1,resfit[2*j+1]);
370	}
371	}
372	}
373	}
374	}
375
376
377	for(int n=0; n< nConstantSet; n++){
378	for(int j=0; j<3; j++){
379	for(int k=0; k< kEndcEtaRings; k++){
380
381	fitHistogram(hh_cc_ietaring[n][j][k],resfit);
382	for(int m=0;m<3; m++){
383	hh_res_cc_ietaring[n][j][m]->SetBinContent(k+1,resfit[2*m]);
384	hh_res_cc_ietaring[n][j][m]->SetBinError(k+1,resfit[2*m+1]);
385	}
386
387	}
388	}
389	}
390
391
392
393	float sigmaSys = 2; ///precision of precalib+LC
394
395
396	cout<<"combine " <<endl;
397	float icwt[2][101][101];
398
399
400	float wtSumC_period[25] = {0};
401	float wtSumS_period[25] = {0};
402
403
404	for(int iz=0; iz<2; iz++){
405	for(int j=0; j<101; j++){
406	for(int k=0; k< 101; k++){
407
408	if( validRecHitEndCap[iz][j][k] <1) {
409	continue;
410	}
411	int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
412
413	float wtSumC = 0;
414	float wtSumS = 0;
415	for(int n=0; n< nIC; n++){ ///for each period
416	wtSumC_period[n] = 0;
417	wtSumS_period[n] = 0;
418	}
419
420	for(int n=0; n< nConstantSet; n++){
421	float c = ic[n][iz][j][k];
422
423
424	float sigma = hh_res_cc_ietaring[n][2][2]->GetBinContent(iring+1);
425	sigma = sqrt( sigma * sigma + sigmaSys * sigmaSys);
426
427	if( c > 0){
428	float tmp1 = c/ ( sigma * sigma);
429	float tmp2 = 1/(sigma * sigma);
430
431	wtSumC += tmp1;
432	wtSumS += tmp2;
433
434	int nperiod = n% nIC;
435	wtSumC_period[nperiod] += tmp1;
436	wtSumS_period[nperiod] += tmp2;
437
438	int deadflag = ndeadflag_ic[n][iz][j][k];
439	if(deadflag>=0){
440	hh_c_deadflag[n][deadflag]->Fill( c);
441	}
442
443	}
444	}
445
446	if( wtSumC > 0){
447	icwt[iz][j][k] = wtSumC / wtSumS;
448
449
450	hh_ccwtavg_ietaring[iz][iring] ->Fill( icwt[iz][j][k] );
451	hh_ccwtavg_ietaring[2][iring] ->Fill( icwt[iz][j][k] );
452
453	}else{
454	icwt[iz][j][k] = -1;
455	}
456	for(int n=0; n< nIC; n++){ ///for each period
457	if( wtSumC_period[n] > 0){
458	icwt_period[n][iz][j][k] = wtSumC_period[n] / wtSumS_period[n];
459
460	hh_ccwtavg_ietaring_period[n][iz][iring]->Fill( icwt_period[n][iz][j][k] );
461	hh_ccwtavg_ietaring_period[n][2][iring]->Fill( icwt_period[n][iz][j][k] );
462
463	}else{
464	icwt_period[n][iz][j][k] = -1;
465	}
466	}
467
468	}
469
470	}
471	}
472
473
474	for(int j=0; j<3; j++){
475	for(int k=0; k< kEndcEtaRings; k++){
476	fitHistogram(hh_ccwtavg_ietaring[j][k],resfit);
477	for(int m=0;m<3; m++){
478	hh_res_ccwtavg_ietaring[j][m]->SetBinContent(k+1,resfit[2*m]);
479	hh_res_ccwtavg_ietaring[j][m]->SetBinError(k+1,resfit[2*m+1]);
480	}
481	}
482	}
483
484	for(int n=0; n< nIC; n++){ ///for each period
485	for(int j=0; j<3; j++){
486	for(int k=0; k< kEndcEtaRings; k++){
487	fitHistogram(hh_ccwtavg_ietaring_period[n][j][k],resfit);
488	for(int m=0;m<3; m++){
489	hh_res_ccwtavg_ietaring_period[n][j][m]->SetBinContent(k+1,resfit[2*m]);
490	hh_res_ccwtavg_ietaring_period[n][j][m]->SetBinError(k+1,resfit[2*m+1]);
491	}
492	}
493	}
494	}
495
496
497	scaleMeanToUnitEndcap(icwt);
498
499	ofstream txtout("interCalibConstants.combined.EcalEndcap.txt",ios::out);
500
501
502	ofstream txtout_period[50]; //pi0eta combined for each period
503	for(int j=0; j< nIC; j++){
504	string filename = string(Form("interCalibConstants.combinedPi0EtaPeriod%d.EcalEndcap.txt",j));
505	txtout_period[j].open(filename.c_str(),ios::out);
506	}
507
508
509	cout<<"print out final " <<endl;
510	for(int iz=0; iz<2; iz++){
511	for(int j=0; j<101; j++){
512	for(int k=0; k< 101; k++){
513	if( validRecHitEndCap[iz][j][k] <1) continue;
514
515
516	int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
517
518	for(int n=0; n< nIC; n++){ ///for each period
519
520	float sigmaC = hh_res_ccwtavg_ietaring_period[n][2][2]->GetBinContent(iring+1);
521
522	float cErr = sqrt( sigmaC sigmaC + sigmaSys sigmaSys);
523	cErr /=100;
524
525	float c = icwt_period[n][iz][j][k];
526	if( c >0){
527
528	int deadflag1 = ndeadflag_ic[n%nIC][iz][j][k];
529	int deadflag2 = ndeadflag_ic[nIC+n%nIC][iz][j][k];
530	if(deadflag1<0 && deadflag2<0){
531	cout<<"wrong deadflag ! " << n << " "<<iz <<" "<< j<<" "<<k <<endl;
532	return;
533	}
534	if(deadflag1>0 ){
535	hh_c_deadflag_period[n][deadflag1]->Fill( c);
536	hh_c_deadflag_period[n][19]->Fill( c);
537	}
538
539	}
540
541	if( c > 0){
542	txtout_period[n]<<j<<" "<<k<<" "<< 2iz-1<<" "<< cccalibpretag[iz][j][k]<<" "<< cErr * c*ccalibpretag[iz][j][k] <<endl;
543
544	}else{
545	txtout_period[n]<<j<<" "<<k<<" "<< 2*iz-1<<" "<<-1 <<" "<< 999 <<endl;
546	}
547	}
548
549
550	float c = icwt[iz][j][k];
551	if( c>0){
552	hh2_cc_xy[nConstantSet][iz]->SetBinContent(j,k,c);
553
554	float sigmaC = hh_res_ccwtavg_ietaring[2][2]->GetBinContent(iring+1);
555
556	float cErr = sqrt( sigmaC sigmaC + sigmaSys sigmaSys);
557	cErr /=100;
558	txtout<<j<<" "<<k<<" "<< 2iz-1<<" "<< c ccalibpretag[iz][j][k] <<" "<< cErr * c * ccalibpretag[iz][j][k] <<endl;
559	}else{
560	txtout<<j<<" "<<k<<" "<< 2*iz-1<<" "<<-1 <<" "<< 999 <<endl;
561	}
562	}
563	}
564	}
565
566
567	fnew->Write();
568	fnew->Close();
569
570
571
572	}