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 *
#	User	Rev	Content
1	yangyong	1.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	yangyong	1.2
12
13			float ccalibpretag[2][101][101];
14	yangyong	1.1
15			float interCalib_preCalib[170][360];
16			float interCalibEndcap_preCalib[2][101][101];
17
18
19	yangyong	1.2
20	yangyong	1.1	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	yangyong	1.2
85
86	yangyong	1.1	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	yangyong	1.2	getInterCalibEndcapv1("interCalibEE_GR_P_V39.txt",ccalibpretag);
102	yangyong	1.1
103
104			get_xyzEBrechits();
105			setEtaRingBoundaryEndcap();
106
107			map<int,string> icFiles;
108
109
110	yangyong	1.2	icFiles[0] = "calibres/deriveCalibConst.dflag2.pe1.step2.iter50.txt";
111			icFiles[1] = "calibres/deriveCalibConst.dflag3.pe2.step2.iter50.txt";
112	yangyong	1.1
113	yangyong	1.2	int nIC = 1;
114	yangyong	1.1
115	yangyong	1.2	getCrystaldeadflagEndcap_v1("deadflag/crystal_deadflag_ee_dflag2.txt",ndeadflag_ic[0]);
116			getCrystaldeadflagEndcap_v1("deadflag/crystal_deadflag_ee_dflag3.txt",ndeadflag_ic[1]);
117	yangyong	1.1
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	yangyong	1.2
133	yangyong	1.1
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	yangyong	1.2	TH1F *hh_c_deadflag[51][30];
244			TH1F *hh_c_deadflag_period[51][30];
245	yangyong	1.1	for(int j=0; j< nConstantSet+1; j++){
246	yangyong	1.2	for(int k=0; k<30; k++){
247	yangyong	1.1	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	yangyong	1.2
265
266	yangyong	1.1	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	yangyong	1.2
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	yangyong	1.1	}
286			}
287			}
288	yangyong	1.2
289
290
291	yangyong	1.1	}
292
293	yangyong	1.2
294			cout<<" nConstantSet " << nConstantSet <<endl;
295
296
297	yangyong	1.1	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	yangyong	1.2
307	yangyong	1.1	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	yangyong	1.2
393			float sigmaSys = 2; ///precision of precalib+LC
394	yangyong	1.1
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	yangyong	1.2	sigma = sqrt( sigma * sigma + sigmaSys * sigmaSys);
426
427	yangyong	1.1	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	yangyong	1.2	if(deadflag>=0){
440			hh_c_deadflag[n][deadflag]->Fill( c);
441	yangyong	1.1	}
442	yangyong	1.2
443	yangyong	1.1	}
444			}
445
446			if( wtSumC > 0){
447			icwt[iz][j][k] = wtSumC / wtSumS;
448
449	yangyong	1.2
450	yangyong	1.1	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	yangyong	1.2
508	yangyong	1.1
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	yangyong	1.2
522			float cErr = sqrt( sigmaC sigmaC + sigmaSys sigmaSys);
523	yangyong	1.1	cErr /=100;
524
525			float c = icwt_period[n][iz][j][k];
526			if( c >0){
527
528	yangyong	1.2	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	yangyong	1.1	}
534	yangyong	1.2	if(deadflag1>0 ){
535			hh_c_deadflag_period[n][deadflag1]->Fill( c);
536	yangyong	1.1	hh_c_deadflag_period[n][19]->Fill( c);
537			}
538	yangyong	1.2
539	yangyong	1.1	}
540
541			if( c > 0){
542	yangyong	1.2	txtout_period[n]<<j<<" "<<k<<" "<< 2iz-1<<" "<< cccalibpretag[iz][j][k]<<" "<< cErr * c*ccalibpretag[iz][j][k] <<endl;
543	yangyong	1.1
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	yangyong	1.2	float cErr = sqrt( sigmaC sigmaC + sigmaSys sigmaSys);
557	yangyong	1.1	cErr /=100;
558	yangyong	1.2	txtout<<j<<" "<<k<<" "<< 2iz-1<<" "<< c ccalibpretag[iz][j][k] <<" "<< cErr * c * ccalibpretag[iz][j][k] <<endl;
559	yangyong	1.1	}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	yangyong	1.2
572	yangyong	1.1	}