yangyong/combineICEE/combineCalibConstantEndcapv1.C

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

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

TRandom3 *rgen_;

#include "usefullcode.cc"

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