yangyong/combineICEB/getCalibConstants.cc



void scaleMeanToUnit(float C[170][360]){

  
  double mean = 0; 
  int n = 0; 
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      
      if( C[j][k] >0){
        mean += C[j][k];
        n ++; 
      } 
    }
  }
  mean /= n;
  cout<<"mean " << mean<<" "<<n<<endl; 

  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      if(C[j][k] >0){
        C[j][k] /= mean; 
      }
    }
  }
  
}


// void scaleMeanToUnitEndcap(float C[2][101][101]){

  
//   double mean = 0; 
//   int n = 0; 
//   for(int j=0; j< 170; j++){
//     for(int k=0; k< 360; k++){
      
//       if( C[j][k] >0){
//      mean += C[j][k];
//      n ++; 
//       }      
//     }
//   }
//   mean /= n;
//   cout<<"mean " << mean<<" "<<n<<endl; 

//   for(int j=0; j< 170; j++){
//     for(int k=0; k< 360; k++){
//       if(C[j][k] >0){
//      C[j][k] /= mean; 
//       }
//     }
//   }
  
// }


void NormIetaToUnitTestBeamSMOnly(float C[170][360]){
  
  float mean_ieta[170] = {0}; 
  int nmean_ieta[170] = {0};

  float mean_ietaTT[170] = {0}; 
  int nmean_ietaTT[170] = {0};
  

  for(int j=0; j< 170; j++){
    
    int ieta = j-85; 
    if( ieta >=0) ieta += 1; 
    
    for(int k=0; k< 360; k++){
            
      int iphi = k; 
      if( k==0) iphi = 360; 
      
      int iSM = (iphi-1)/20+1;
      if( ieta<0) iSM += 18;
      int smTB = isTestBeamSM(iSM);
      if( smTB ==0) continue; 
      
      //if( C[j][k] >0 ){
      if( C[j][k] >0.5 && C[j][k] < 1.5 ){
        mean_ieta[j] += C[j][k];
        nmean_ieta[j] ++; 

        mean_ietaTT[j/5] += C[j][k];
        nmean_ietaTT[j/5] ++; 
        
      }
      
    }
  }
  
  
  for(int j=0; j< 170; j++){
    mean_ieta[j] /= nmean_ieta[j];
    //cout<<" mean_ieta[" << j<<"] "<<  mean_ieta[j] <<endl; 
  }
   
//   for(int j=0; j< 34; j++){
//     mean_ietaTT[j] /= nmean_ietaTT[j];
//     cout<<" mean_ietaTT[" << j<<"] "<<  mean_ietaTT[j] <<endl; 
//   }
  
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      
      if(  C[j][k] >0 ){
        C[j][k] /=  mean_ieta[j] ; 
      }
    }
  }
  
}


void NormIetaAbsToUnitTestBeamSMOnly(float C[170][360]){
  
  
  float mean_ieta[170] = {0}; 
  int nmean_ieta[170] = {0};

  float mean_ietaTT[170] = {0}; 
  int nmean_ietaTT[170] = {0};
  

  for(int j=0; j< 170; j++){
    
    int ieta = j-85; 
    if( ieta >=0) ieta += 1; 
    
    for(int k=0; k< 360; k++){
            
      int iphi = k; 
      if( k==0) iphi = 360; 
      
      int iSM = (iphi-1)/20+1;
      if( ieta<0) iSM += 18;
      int smTB = isTestBeamSM(iSM);
      if( smTB ==0) continue; 
      
      int ietaAbs = abs(ieta);
      
      //if( C[j][k] >0 ){
      if( C[j][k] >0.5 && C[j][k] < 1.5 ){
        mean_ieta[ietaAbs-1] += C[j][k];
        nmean_ieta[ietaAbs-1] ++; 
        
        mean_ietaTT[j/5] += C[j][k];
        nmean_ietaTT[j/5] ++; 
        
      }
    }
  }
  
  
  for(int j=0; j< 85; j++){
    mean_ieta[j] /= nmean_ieta[j];
    //cout<<" mean_ietaAbs[" << j<<"] "<<  mean_ieta[j] <<endl; 
  }
  
//   for(int j=0; j< 34; j++){
//     mean_ietaTT[j] /= nmean_ietaTT[j];
//     cout<<" mean_ietaTT[" << j<<"] "<<  mean_ietaTT[j] <<endl; 
//   }
  
  for(int j=0; j< 170; j++){
    
    int ieta = j-85; 
    if( ieta >=0) ieta += 1; 
    
    for(int k=0; k< 360; k++){
      
      int ietaAbs = abs(ieta);
      
      if(  C[j][k] >0 ){
                
        C[j][k] /=  mean_ieta[ietaAbs-1] ; 
        
      }
    }
  }
  
}


void NormCrystalDeadFlag_v1(float C[170][360], int ndeadflagietaiphi[170][360]){
  
  float mean_deadflag[20] = {0};
  float nmean_deadflag[20] = {0};
                       
  int ndeadCrystals = 0; 
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      if( C[j][k] <0){
        ndeadCrystals++; 
        continue; 
      }
      int ndeadflag = ndeadflagietaiphi[j][k];
      if( ndeadflag >=0){
        mean_deadflag[ndeadflag] += C[j][k];
        nmean_deadflag[ndeadflag] ++; 
      }
    }
  }
  cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl; 
  for(int j=0; j<20; j++){
    if( nmean_deadflag[j] >0){
      mean_deadflag[j] /= nmean_deadflag[j];
      cout<<" mean_deadflag "<<j<<"  "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;  
    }
  }
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      if( C[j][k] <0) continue;
      int ndeadflag = ndeadflagietaiphi[j][k];
      if( ndeadflag >=0){
        C[j][k] /= mean_deadflag[ndeadflag];
      }
    }
  }
  
}

void NormCrystalDeadFlag(float C[170][360]){
  
  float mean_deadflag[20] = {0};
  float nmean_deadflag[20] = {0};
                       
  int ndeadCrystals = 0; 
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      if( C[j][k] <0){
        ndeadCrystals++; 
        continue; 
      }
      int ndeadflag = ndeadflag_ietaiphi[j][k];
      if( ndeadflag >=0){
        mean_deadflag[ndeadflag] += C[j][k];
        nmean_deadflag[ndeadflag] ++; 
      }
    }
  }
  cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl; 
  for(int j=0; j<20; j++){
    if( nmean_deadflag[j] >0){
      mean_deadflag[j] /= nmean_deadflag[j];
      cout<<" mean_deadflag "<<j<<"  "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;  
    }
  }
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      if( C[j][k] <0) continue;
      int ndeadflag = ndeadflag_ietaiphi[j][k];
      if( ndeadflag >=0){
        C[j][k] /= mean_deadflag[ndeadflag];
      }
    }
  }
  
}


void SetSMScale(float C[170][360],TH1F *hhtmp){
  float mean_ism[36] = {0};
  float nmean_ism[36]= {0};

  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){

      int ieta = j-85;
      if( j >=0) ieta += 1;
      int iphi = k;
      if( k==0) iphi = 360;
      int iSM = (iphi-1)/20+1;
      if( ieta<0) iSM += 18;
      
      if(  C[j][k] >0 ){
        mean_ism[iSM-1 ]  += C[j][k];
        nmean_ism[iSM-1] ++;
      }
    }
  }

  for(int j=0; j<36; j++){
    mean_ism[j] /=  nmean_ism[j];
    ///    cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl;
  }
  
  
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){

      int ieta = j-85;
      if( j >=0) ieta += 1;
      int iphi = k;
      if( k==0) iphi = 360;
      int iSM = (iphi-1)/20+1;
      if( ieta<0) iSM += 18;
      
      if(  C[j][k] >0 ){
        C[j][k] *= hhtmp->GetBinContent(iSM);
      }
    }
  }
  
  
}


void NormSMScaleToUnit(float C[170][360]){
  float mean_ism[36] = {0}; 
  float nmean_ism[36]= {0};
  
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){
      
      int ieta = j-85; 
      if( j >=0) ieta += 1; 
      int iphi = k; 
      if( k==0) iphi = 360; 
      int iSM = (iphi-1)/20+1;
      if( ieta<0) iSM += 18;
      
      
      if(  C[j][k] >0 ){
        mean_ism[iSM-1 ]  += C[j][k]; 
        nmean_ism[iSM-1] ++; 
      }
    }
  }
  
  for(int j=0; j<36; j++){
    mean_ism[j] /=  nmean_ism[j];
    //cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl; 
  }
  
  for(int j=0; j< 170; j++){
    for(int k=0; k< 360; k++){

      int ieta = j-85; 
      if( j >=0) ieta += 1; 
      int iphi = k; 
      if( k==0) iphi = 360; 
      int iSM = (iphi-1)/20+1;
      if( ieta<0) iSM += 18;
      
      if(  C[j][k] >0 ){
        C[j][k] /=  mean_ism[iSM-1];
      }
      
    }
  }
  
}


void readInterCalibConstEBSimple(const char *input,float C[170][360]){
  
  ifstream txtin1(input,ios::in);
  if (txtin1.fail()){
    cout<<"error open file barrel.. " << input<<endl;
    exit(1);
  }
  
  
  int eta;
  int phi;
  
  float cc; 

  int n =0; 
  while(txtin1.good()){
    
    
    txtin1 >> eta >> phi >> cc; 

    int iphi = phi;
    int ieta = eta;

    // cout<< eta <<" "<< phi <<" "<< cc <<endl; 
    
    convxtalid(phi,eta);
    eta += 85;
    C[eta][phi] = cc; 
    
    n++; 
    if( n>= 61200) break; 
    
  }
  
}


void readInterCalibConstEBSimplev1(const char *input,float C[170][360]){
  
  ifstream txtin1(input,ios::in);
  if (txtin1.fail()){
    cout<<"error open file barrel.. " << input<<endl;
    exit(1);
  }
  
  
  int eta;
  int phi;
  
  float cc; 
  double mean = 0; 
  int ngood = 0; 
  int n =0; 
  while(txtin1.good()){
    
    txtin1 >> eta >> phi >> cc; 

    int iphi = phi;
    int ieta = eta;
    
    if( !(eta>=0 && eta <= 169 && phi >=0 && phi<=359)) {
      cout<<"wrong eta/phi " << eta<<" "<<phi <<endl; 
      exit(1);
    }

    C[eta][phi] = cc; 
    if( cc >0){
      mean += cc; 
      ngood ++; 
    }
    n++; 
    if( n>= 61200) break; 
  }
  mean /= ngood; 
  cout<<"read ngood "<< ngood <<" mean: "<< mean<<endl; 
  
  for(int j=0; j<170; j++){
    for(int k=0; k<360;k++){
      if(C[j][k] >0) {
        C[j][k] /= mean;
      }
    }
  }
  
}


void readInterCalibBarre_Phiv1new(){
  
  char *file_input = new char[500];
  sprintf(file_input,"phisym_EcalIntercalib-160M-corrected.txt");
  
  sprintf(file_input,"EcalIntercalibConstants_Run2010A_corr.xml.txt");
  
  cout<<"READING fro inter-calib file "<<file_input<<endl;
  ifstream inputcc(file_input,ios::in);
  if (inputcc.fail()){
    cout<<"error open file barrel.. " << file_input<<endl;
    exit(1);
  }
  
  float c1; 
  for(int ieta=-85; ieta <= 85; ieta++){
    if( ieta==0) continue; 
    for(int iphi=1; iphi <= 360; iphi++){
      inputcc >> c1; 
      
      int eta = ieta; 
      int phi = iphi; 
      convxtalid(phi,eta);
      
      if( c1 <=0){
        c1 =1; 
      }
      
      if( c1 == 1 ){
        c1 = - 1; 
      }
      if( c1 <0.6 || c1 >2.5){
        c1 = -1; 
      }
      
      
      CphiCorr[eta+85][phi] = c1; 
    }
  }
    
}


void readInterCalibBarre_Phi2011a(){
  
  char *file_input = new char[500];

  
  sprintf(file_input,"EcalIntercalib-2011A-corrected.xmlv1");
  
  cout<<"READING fro inter-calib file "<<file_input<<endl;
  ifstream inputcc(file_input,ios::in);
  if (inputcc.fail()){
    cout<<"error open file barrel.. " << file_input<<endl;
    exit(1);
  }
  
  float c1; 
  for(int ieta=-85; ieta <= 85; ieta++){
    if( ieta==0) continue; 
    for(int iphi=1; iphi <= 360; iphi++){
      inputcc >> c1; 
      
      int eta = ieta; 
      int phi = iphi; 
      convxtalid(phi,eta);
      
      if( c1 <=0){
        c1 =1; 
      }
      
      if( c1 == 1 ){
        c1 = - 1; 
      }
      if( c1 <0.6 || c1 >2.5){
        c1 = -1; 
      }
      
      
      CphiCorr[eta+85][phi] = c1; 
    }
  }
  

}


void readInterCalibBarre_Phiv1(){
    
  char *file_input = new char[500];
  sprintf(file_input,"phisym_EcalIntercalib-160M-corrected.txt");
  
  
  cout<<"READING fro inter-calib file "<<file_input<<endl;
  ifstream inputcc(file_input,ios::in);
  if (inputcc.fail()){
    cout<<"error open file barrel.. " << file_input<<endl;
    exit(1);
  }
  
  int eta; 
  int phi; 
  float c = 1; 
  float c1 = 1; 
  
  int n = 0; 

  eta  = -85; 
  phi = 0; 
  
  float mean = 0; 
  
  int cureta = -85;
  
  int ngood = 0; 
  
  
  while (inputcc.good()){

    inputcc >> c1; 

    if(c1==1) c1 = -1; 
    
    
    if(  n% 360 ==0 && n>0 ) {
      eta++; 
      if(eta==0) eta = 1;  
    }
    
    if( eta == cureta){
      phi ++; 
    }else{
      cureta = eta; 
      phi = 1; 
      
    }
    
    
    n++; 
    
    
    int iSM = (phi-1)/20+1;
    if( eta<0) iSM += 18;
        
    
    int ieta = eta; 
    int iphi = phi; 
    
    
    convxtalid(iphi,ieta);
    
    ///txtout<<eta<<" "<<phi<<endl; 
    
    CphiCorr[ieta+85][iphi] = c1;  
    
    
    if(c1 >0 ){
      mean += c1; 
      ngood ++; 
    }
    
    
    if(n >= 360*170) break; 
    
  }
  
  
  mean /=  ngood; 
  cout<<n<<" phi  constants read. "<< mean <<endl; 
  
  
//   for(int j=0; j< 170; j++){
//     mean_ieta[j] /= nmean_ieta[j];
//     txtout<<"mean_ieta_phisym "<< j<<" "<< mean_ieta[j] <<" "<< nmean_ieta[j] <<endl; 
//   }
  
  
//   for(int j=0; j< 170; j++){
//     for(int k=0; k< 360; k++){
      
//       if(  CphiCorr[j][k] >0 ){
//      CphiCorr[j][k] /=  mean_ieta[j] ; 
        
//       }
//       Cphi[j][k] =   CphiCorr[j][k] *  iCalConst_data[j][k];
//     }
//   }
  
//   NormSMScaleToUnit(CphiCorr);
  
}


///read initi calibration constants for each crystal
void readInterCalibBarre_BSv1(){
  
  
  char *file_input = new char[500];
  sprintf(file_input,"IC_barrel_splash09_AllRuns.txt");
  
  
  cout<<"READING fro inter-calib bsfile "<<file_input<<endl;
  txtout<<"READING fro inter-calib bsfile "<<file_input<<endl;
  ifstream inputcc(file_input,ios::in);
  if (inputcc.fail()){
    cout<<"error open file barrel.. " << file_input<<endl;
    exit(1);
  }
  
  int eta; 
  int phi; 
  float c = 1; 
  float c1 = 1; 
  
  int n = 0; 
  int ngood = 0; 
  

  float mean = 0; 
   
  
  float mean_ieta[170] = {0}; 
  int nmean_ieta[170] = {0};
  

  while (inputcc.good()){

    inputcc >> eta >> phi >>c1; 
    
    //  cout<<phi<<" "<<eta<<" "<<c<<endl;
    
    convxtalid(phi,eta);
    
    eta += 85; 
    
    if( c1 == 1){
      c1 = -1; 
    }
    
    if( c1 <0.6 || c1 >2.5){
      c1 = -1; 
    }
    
    
    ///iCalConstPi0[eta][phi] = c1 ; 
    
    CBSCorr[eta][phi] = c1;  
    
    if(c1 >0 ){
      mean += c1; 
      ngood ++;
    }

    if(c1 >0 ){
      mean_ieta[eta] += c1; 
      nmean_ieta[eta] ++; 
    }
    
    n++; 
    
    /// hh_imcEB->Fill(c);
    if(n >= 360*170) break; 
    
  }
  
  
  mean /= ngood; 
  
  
  cout<<n<<" bs  constants read. "<< mean << endl;
  
    
//   for(int j=0; j< 170; j++){
//     for(int k=0; k< 360; k++){
      
//       if(  CBSCorr[j][k] >0 ){
//      CBSCorr[j][k] /=  mean_ieta[j];
//       }
      
//       CBS[j][k] = CBSCorr[j][k] * iCalConst_data[j][k];
      
      
//     }
//   }
//   NormSMScaleToUnit(CBSCorr);
  
  
}


///read initi calibration constants for each crystal
void readInterCalibBarre_GR09_V8(){


  char *file_input = new char[500];
  //file_input ="/data/yangyong/data/p
  file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_R_V6A_barrel.txt";
  
  ///sprintf(file_input,"/uscms_data/d2/yongy/calib/backup/interCalib_STARTUP3X_V8P_barrel_mc.txt");  ///now the same used in data
  
  //sprintf(file_input,"/uscms/home/yongy/work/interCalib_STARTUP3X_V8P_barrel.txt");
  
  //  sprintf(file_input,"/uscms/home/yongy/work/interCalib_GR_R_36X_V11A_Barrel.txt");

  cout<<"READING fro inter-calib file "<<file_input<<endl;
  ifstream inputcc(file_input,ios::in);
  if (inputcc.fail()){
    cout<<"error open file barrel.. " << file_input<<endl;
    exit(1);
  }
  
  int eta; 
  int phi; 
  float c = 1; 
  float c1 = 1; 
  
  int n = 0; 
  
  double mean = 0; 
  

  while (inputcc.good()){

    inputcc >> eta >> phi >>c1; 
    
    //  cout<<phi<<" "<<eta<<" "<<c<<endl;
    
    convxtalid(phi,eta);
    
    eta += 85; 
    
    ///iCalConst[eta][phi] = c; 
    
    //iCalConst_data[eta][phi] = c1 * 0.038012/ 0.038;  ///Now this is exactly the same as the inter-calibration used in data. 
    
    interCalib_preCalib[eta][phi] = c1 ; /// now OK iN MC startupv8p. 
    
    
    n++; 
    //  cout<<phi<<" "<<eta<<" "<<c<<endl;
    mean += c1; 
    
    /// hh_imcEB->Fill(c);
    if(n >= 360*170) break; 
    
  }
  
 
  mean /= 61200; 
    
  cout<<n<<" constants read. "<< mean <<endl; 
  
  
}


///read initi calibration constants for each crystal
void readInterCalibEndcap_GR09_V8(){
  char *file_input = new char[500];
  //file_input ="/data/yangyong/data/p
  ///file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_R_V6A_barrel.txt";
  
  file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_H_V6OFF_endcap.txt";
    
  cout<<"READING fro inter-calib file "<<file_input<<endl;
  ifstream inputcc(file_input,ios::in);
  if (inputcc.fail()){
    cout<<"error open file barrel.. " << file_input<<endl;
    exit(1);
  }
  
  int eta; 
  int phi; 
  int iz; 
  float c = 1; 
  int n = 0; 
  
  
  while (inputcc.good()){
    
    inputcc >> iz >> eta >> phi  >> c; 
    
    int izz = iz < 0? 0:1; 

    interCalibEndcap_preCalib[izz][eta][phi] = c; 
    validRecHitEndCap[izz][eta][phi] = 1; 
    
    n ++; 
    
  }
  
  cout<<n<<" constants read. "<<endl;
  
}


void getInterCalibEndcapv1(const char *inputfile,float C[2][101][101]){

  
  ifstream inputcc(inputfile);
  if (inputcc.fail()) {
    cout<<"error "<< inputfile<<" can not be opened."<<endl; 
    exit(1);
  }
  int eta; 
  int phi; 
  int iz; 
  float c = 1; 
  int n = 0; 
  
  double mean = 0; 
  
  while (inputcc.good()){

    inputcc >> iz >> eta >> phi  >> c; 
    
    int izz = iz < 0? 0:1; 

    C[izz][eta][phi] = c; 
    mean += c; 
    
    n ++; 

    if( n>= 14648){
      break; 
    }
    
  }
  
  cout<<n<<" constants read. "<< mean /n <<endl;
  
}


void getInterCalibEndcap(const char *inputfile, float C[2][101][101]){
  
  
  ifstream inputcc(inputfile);
  if (inputcc.fail()) {
    cout<<"error "<< inputfile<<" can not be opened."<<endl; 
    exit(1);
  }
  
  int ix,iy,iz; 
  float cc; 

  float ccMeaniring[50] ;
  int nccMeaniring[50];
  

//  for(int j=0; j<2;j++){
  for(int k=0; k< kEndcEtaRings; k++){
    ccMeaniring[k] = 0; 
    nccMeaniring[k] = 0; 
  }
  //  }
  
  int nread = 0; 
  while( inputcc.good()){
    inputcc >> iz >> ix >> iy >> cc; 
    C[iz][ix][iy] = cc; 
    int iring = iRingEndCap(2*iz-1,ix,iy); ///input -1/1, 
    if( cc > 0){
      ccMeaniring[iring] += cc; 
      nccMeaniring[iring] += 1; 
    }
    
    nread ++; 
    if( iz==1 && ix == 100 && iy == 60) break; 
  }
  
  cout<<"nb of crystals endcap " << nread <<endl; 
  
  
  //for(int j=0; j<2;j++){
  for(int k=0; k< kEndcEtaRings; k++){
    ccMeaniring[k] /= nccMeaniring[k];
    //  }
  }
  
  double meanall = 0; 
  int nall = 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,                      
        if( C[iz][j][k] >0) {
          C[iz][j][k] /= ccMeaniring[iring];
          meanall += C[iz][j][k];
          nall ++; 
        }
        
      }
    }
  }
  
  meanall /= nall; 
  cout<<" corrfactoriZiXiY read " << meanall <<" "<< nall <<endl;
  
  
}


void NormCrystalDeadFlagEndcap_v1(float C[2][101][101],int ndeadflag[2][101][101]){
  
  float mean_deadflag[20] = {0};
  float nmean_deadflag[20] = {0};
                       
  int ndeadCrystals = 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] ==0) continue; 
        if( C[iz][j][k] <0){
          continue; 
        }
        int ndead = ndeadflag[iz][j][k];
        if( ndead >=0){
          mean_deadflag[ndead] += C[iz][j][k];
          nmean_deadflag[ndead] ++; 
        }
        
      }
    }
  }
  cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl; 
  for(int j=0; j<20; j++){
    if( nmean_deadflag[j] >0){
      mean_deadflag[j] /= nmean_deadflag[j];
      cout<<" mean_deadflag "<<j<<"  "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<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] ==0) continue; 
        if( C[iz][j][k] <0){
          continue; 
        }
        
        int ndead = ndeadflag[iz][j][k];
        if( ndead >=0){
          C[iz][j][k] /= mean_deadflag[ndead];
        }
        
      }
    }
  }
  
  
}


void NormCrystalDeadFlagEndcap(float C[2][101][101]){
  
  float mean_deadflag[20] = {0};
  float nmean_deadflag[20] = {0};
                       
  int ndeadCrystals = 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] ==0) continue; 
        if( C[iz][j][k] <0){
          continue; 
        }
        int ndeadflag = ndeadflag_endcap[iz][j][k];
        if( ndeadflag >=0){
          mean_deadflag[ndeadflag] += C[iz][j][k];
          nmean_deadflag[ndeadflag] ++; 
        }
        
      }
    }
  }
  cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl; 
  for(int j=0; j<20; j++){
    if( nmean_deadflag[j] >0){
      mean_deadflag[j] /= nmean_deadflag[j];
      cout<<" mean_deadflag "<<j<<"  "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<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] ==0) continue; 
        if( C[iz][j][k] <0){
          continue; 
        }
        
        int ndeadflag = ndeadflag_endcap[iz][j][k];
        if( ndeadflag >=0){
          C[iz][j][k] /= mean_deadflag[ndeadflag];
        }
        
      }
    }
  }
  
  
}


void scaleMeanToUnitEndcap(float C[2][101][101]){
  
  
  double mean = 0; 
  int n = 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] ==0) continue; 
        if( C[iz][j][k] <0){
          continue; 
        }
        mean += C[iz][j][k]; 
        n ++; 
      }
    }
  }
    
  mean /= n;
  cout<<"mean " << mean<<" "<<n<<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] ==0) continue; 
        if( C[iz][j][k] <0){
          continue; 
        }
        
        C[iz][j][k] /= mean; 
      }
    }
  }
  
}
Revision:	1.2
Committed:	Tue Apr 10 19:41:08 2012 UTC (13 years, 1 month ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V2012combv2b, V2012combv2a, V2012combv2, V2012combv1, V2011combv1, pi0etaeb_laser20111122
Changes since 1.1:	+1147 -0 lines
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.2
2
3			void scaleMeanToUnit(float C[170][360]){
4
5
6			double mean = 0;
7			int n = 0;
8			for(int j=0; j< 170; j++){
9			for(int k=0; k< 360; k++){
10
11			if( C[j][k] >0){
12			mean += C[j][k];
13			n ++;
14			}
15			}
16			}
17			mean /= n;
18			cout<<"mean " << mean<<" "<<n<<endl;
19
20			for(int j=0; j< 170; j++){
21			for(int k=0; k< 360; k++){
22			if(C[j][k] >0){
23			C[j][k] /= mean;
24			}
25			}
26			}
27
28			}
29
30
31
32
33			// void scaleMeanToUnitEndcap(float C[2][101][101]){
34
35
36			// double mean = 0;
37			// int n = 0;
38			// for(int j=0; j< 170; j++){
39			// for(int k=0; k< 360; k++){
40
41			// if( C[j][k] >0){
42			// mean += C[j][k];
43			// n ++;
44			// }
45			// }
46			// }
47			// mean /= n;
48			// cout<<"mean " << mean<<" "<<n<<endl;
49
50			// for(int j=0; j< 170; j++){
51			// for(int k=0; k< 360; k++){
52			// if(C[j][k] >0){
53			// C[j][k] /= mean;
54			// }
55			// }
56			// }
57
58			// }
59
60
61
62
63			void NormIetaToUnitTestBeamSMOnly(float C[170][360]){
64
65			float mean_ieta[170] = {0};
66			int nmean_ieta[170] = {0};
67
68			float mean_ietaTT[170] = {0};
69			int nmean_ietaTT[170] = {0};
70
71
72			for(int j=0; j< 170; j++){
73
74			int ieta = j-85;
75			if( ieta >=0) ieta += 1;
76
77			for(int k=0; k< 360; k++){
78
79			int iphi = k;
80			if( k==0) iphi = 360;
81
82			int iSM = (iphi-1)/20+1;
83			if( ieta<0) iSM += 18;
84			int smTB = isTestBeamSM(iSM);
85			if( smTB ==0) continue;
86
87			//if( C[j][k] >0 ){
88			if( C[j][k] >0.5 && C[j][k] < 1.5 ){
89			mean_ieta[j] += C[j][k];
90			nmean_ieta[j] ++;
91
92			mean_ietaTT[j/5] += C[j][k];
93			nmean_ietaTT[j/5] ++;
94
95			}
96
97			}
98			}
99
100
101			for(int j=0; j< 170; j++){
102			mean_ieta[j] /= nmean_ieta[j];
103			//cout<<" mean_ieta[" << j<<"] "<< mean_ieta[j] <<endl;
104			}
105
106			// for(int j=0; j< 34; j++){
107			// mean_ietaTT[j] /= nmean_ietaTT[j];
108			// cout<<" mean_ietaTT[" << j<<"] "<< mean_ietaTT[j] <<endl;
109			// }
110
111			for(int j=0; j< 170; j++){
112			for(int k=0; k< 360; k++){
113
114			if( C[j][k] >0 ){
115			C[j][k] /= mean_ieta[j] ;
116			}
117			}
118			}
119
120			}
121
122
123
124			void NormIetaAbsToUnitTestBeamSMOnly(float C[170][360]){
125
126
127			float mean_ieta[170] = {0};
128			int nmean_ieta[170] = {0};
129
130			float mean_ietaTT[170] = {0};
131			int nmean_ietaTT[170] = {0};
132
133
134			for(int j=0; j< 170; j++){
135
136			int ieta = j-85;
137			if( ieta >=0) ieta += 1;
138
139			for(int k=0; k< 360; k++){
140
141			int iphi = k;
142			if( k==0) iphi = 360;
143
144			int iSM = (iphi-1)/20+1;
145			if( ieta<0) iSM += 18;
146			int smTB = isTestBeamSM(iSM);
147			if( smTB ==0) continue;
148
149			int ietaAbs = abs(ieta);
150
151			//if( C[j][k] >0 ){
152			if( C[j][k] >0.5 && C[j][k] < 1.5 ){
153			mean_ieta[ietaAbs-1] += C[j][k];
154			nmean_ieta[ietaAbs-1] ++;
155
156			mean_ietaTT[j/5] += C[j][k];
157			nmean_ietaTT[j/5] ++;
158
159			}
160			}
161			}
162
163
164			for(int j=0; j< 85; j++){
165			mean_ieta[j] /= nmean_ieta[j];
166			//cout<<" mean_ietaAbs[" << j<<"] "<< mean_ieta[j] <<endl;
167			}
168
169			// for(int j=0; j< 34; j++){
170			// mean_ietaTT[j] /= nmean_ietaTT[j];
171			// cout<<" mean_ietaTT[" << j<<"] "<< mean_ietaTT[j] <<endl;
172			// }
173
174			for(int j=0; j< 170; j++){
175
176			int ieta = j-85;
177			if( ieta >=0) ieta += 1;
178
179			for(int k=0; k< 360; k++){
180
181			int ietaAbs = abs(ieta);
182
183			if( C[j][k] >0 ){
184
185			C[j][k] /= mean_ieta[ietaAbs-1] ;
186
187			}
188			}
189			}
190
191			}
192
193
194
195			void NormCrystalDeadFlag_v1(float C[170][360], int ndeadflagietaiphi[170][360]){
196
197			float mean_deadflag[20] = {0};
198			float nmean_deadflag[20] = {0};
199
200			int ndeadCrystals = 0;
201			for(int j=0; j< 170; j++){
202			for(int k=0; k< 360; k++){
203			if( C[j][k] <0){
204			ndeadCrystals++;
205			continue;
206			}
207			int ndeadflag = ndeadflagietaiphi[j][k];
208			if( ndeadflag >=0){
209			mean_deadflag[ndeadflag] += C[j][k];
210			nmean_deadflag[ndeadflag] ++;
211			}
212			}
213			}
214			cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
215			for(int j=0; j<20; j++){
216			if( nmean_deadflag[j] >0){
217			mean_deadflag[j] /= nmean_deadflag[j];
218			cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
219			}
220			}
221			for(int j=0; j< 170; j++){
222			for(int k=0; k< 360; k++){
223			if( C[j][k] <0) continue;
224			int ndeadflag = ndeadflagietaiphi[j][k];
225			if( ndeadflag >=0){
226			C[j][k] /= mean_deadflag[ndeadflag];
227			}
228			}
229			}
230
231			}
232
233			void NormCrystalDeadFlag(float C[170][360]){
234
235			float mean_deadflag[20] = {0};
236			float nmean_deadflag[20] = {0};
237
238			int ndeadCrystals = 0;
239			for(int j=0; j< 170; j++){
240			for(int k=0; k< 360; k++){
241			if( C[j][k] <0){
242			ndeadCrystals++;
243			continue;
244			}
245			int ndeadflag = ndeadflag_ietaiphi[j][k];
246			if( ndeadflag >=0){
247			mean_deadflag[ndeadflag] += C[j][k];
248			nmean_deadflag[ndeadflag] ++;
249			}
250			}
251			}
252			cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
253			for(int j=0; j<20; j++){
254			if( nmean_deadflag[j] >0){
255			mean_deadflag[j] /= nmean_deadflag[j];
256			cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
257			}
258			}
259			for(int j=0; j< 170; j++){
260			for(int k=0; k< 360; k++){
261			if( C[j][k] <0) continue;
262			int ndeadflag = ndeadflag_ietaiphi[j][k];
263			if( ndeadflag >=0){
264			C[j][k] /= mean_deadflag[ndeadflag];
265			}
266			}
267			}
268
269			}
270
271
272
273			void SetSMScale(float C[170][360],TH1F *hhtmp){
274			float mean_ism[36] = {0};
275			float nmean_ism[36]= {0};
276
277			for(int j=0; j< 170; j++){
278			for(int k=0; k< 360; k++){
279
280			int ieta = j-85;
281			if( j >=0) ieta += 1;
282			int iphi = k;
283			if( k==0) iphi = 360;
284			int iSM = (iphi-1)/20+1;
285			if( ieta<0) iSM += 18;
286
287			if( C[j][k] >0 ){
288			mean_ism[iSM-1 ] += C[j][k];
289			nmean_ism[iSM-1] ++;
290			}
291			}
292			}
293
294			for(int j=0; j<36; j++){
295			mean_ism[j] /= nmean_ism[j];
296			/// cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl;
297			}
298
299
300			for(int j=0; j< 170; j++){
301			for(int k=0; k< 360; k++){
302
303			int ieta = j-85;
304			if( j >=0) ieta += 1;
305			int iphi = k;
306			if( k==0) iphi = 360;
307			int iSM = (iphi-1)/20+1;
308			if( ieta<0) iSM += 18;
309
310			if( C[j][k] >0 ){
311			C[j][k] *= hhtmp->GetBinContent(iSM);
312			}
313			}
314			}
315
316
317			}
318
319
320
321
322
323
324
325			void NormSMScaleToUnit(float C[170][360]){
326			float mean_ism[36] = {0};
327			float nmean_ism[36]= {0};
328
329			for(int j=0; j< 170; j++){
330			for(int k=0; k< 360; k++){
331
332			int ieta = j-85;
333			if( j >=0) ieta += 1;
334			int iphi = k;
335			if( k==0) iphi = 360;
336			int iSM = (iphi-1)/20+1;
337			if( ieta<0) iSM += 18;
338
339
340			if( C[j][k] >0 ){
341			mean_ism[iSM-1 ] += C[j][k];
342			nmean_ism[iSM-1] ++;
343			}
344			}
345			}
346
347			for(int j=0; j<36; j++){
348			mean_ism[j] /= nmean_ism[j];
349			//cout<<"sm mean: " << j+1 <<" "<< mean_ism[j] <<" "<< nmean_ism[j]<<endl;
350			}
351
352			for(int j=0; j< 170; j++){
353			for(int k=0; k< 360; k++){
354
355			int ieta = j-85;
356			if( j >=0) ieta += 1;
357			int iphi = k;
358			if( k==0) iphi = 360;
359			int iSM = (iphi-1)/20+1;
360			if( ieta<0) iSM += 18;
361
362			if( C[j][k] >0 ){
363			C[j][k] /= mean_ism[iSM-1];
364			}
365
366			}
367			}
368
369			}
370
371
372
373			void readInterCalibConstEBSimple(const char *input,float C[170][360]){
374
375			ifstream txtin1(input,ios::in);
376			if (txtin1.fail()){
377			cout<<"error open file barrel.. " << input<<endl;
378			exit(1);
379			}
380
381
382			int eta;
383			int phi;
384
385			float cc;
386
387			int n =0;
388			while(txtin1.good()){
389
390
391			txtin1 >> eta >> phi >> cc;
392
393			int iphi = phi;
394			int ieta = eta;
395
396			// cout<< eta <<" "<< phi <<" "<< cc <<endl;
397
398			convxtalid(phi,eta);
399			eta += 85;
400			C[eta][phi] = cc;
401
402			n++;
403			if( n>= 61200) break;
404
405			}
406
407			}
408
409
410
411
412
413			void readInterCalibConstEBSimplev1(const char *input,float C[170][360]){
414
415			ifstream txtin1(input,ios::in);
416			if (txtin1.fail()){
417			cout<<"error open file barrel.. " << input<<endl;
418			exit(1);
419			}
420
421
422			int eta;
423			int phi;
424
425			float cc;
426			double mean = 0;
427			int ngood = 0;
428			int n =0;
429			while(txtin1.good()){
430
431			txtin1 >> eta >> phi >> cc;
432
433			int iphi = phi;
434			int ieta = eta;
435
436			if( !(eta>=0 && eta <= 169 && phi >=0 && phi<=359)) {
437			cout<<"wrong eta/phi " << eta<<" "<<phi <<endl;
438			exit(1);
439			}
440
441			C[eta][phi] = cc;
442			if( cc >0){
443			mean += cc;
444			ngood ++;
445			}
446			n++;
447			if( n>= 61200) break;
448			}
449			mean /= ngood;
450			cout<<"read ngood "<< ngood <<" mean: "<< mean<<endl;
451
452			for(int j=0; j<170; j++){
453			for(int k=0; k<360;k++){
454			if(C[j][k] >0) {
455			C[j][k] /= mean;
456			}
457			}
458			}
459
460			}
461
462
463
464			void readInterCalibBarre_Phiv1new(){
465
466			char *file_input = new char[500];
467			sprintf(file_input,"phisym_EcalIntercalib-160M-corrected.txt");
468
469			sprintf(file_input,"EcalIntercalibConstants_Run2010A_corr.xml.txt");
470
471			cout<<"READING fro inter-calib file "<<file_input<<endl;
472			ifstream inputcc(file_input,ios::in);
473			if (inputcc.fail()){
474			cout<<"error open file barrel.. " << file_input<<endl;
475			exit(1);
476			}
477
478			float c1;
479			for(int ieta=-85; ieta <= 85; ieta++){
480			if( ieta==0) continue;
481			for(int iphi=1; iphi <= 360; iphi++){
482			inputcc >> c1;
483
484			int eta = ieta;
485			int phi = iphi;
486			convxtalid(phi,eta);
487
488			if( c1 <=0){
489			c1 =1;
490			}
491
492			if( c1 == 1 ){
493			c1 = - 1;
494			}
495			if( c1 <0.6 \|\| c1 >2.5){
496			c1 = -1;
497			}
498
499
500			CphiCorr[eta+85][phi] = c1;
501			}
502			}
503
504			}
505
506
507			void readInterCalibBarre_Phi2011a(){
508
509			char *file_input = new char[500];
510
511
512			sprintf(file_input,"EcalIntercalib-2011A-corrected.xmlv1");
513
514			cout<<"READING fro inter-calib file "<<file_input<<endl;
515			ifstream inputcc(file_input,ios::in);
516			if (inputcc.fail()){
517			cout<<"error open file barrel.. " << file_input<<endl;
518			exit(1);
519			}
520
521			float c1;
522			for(int ieta=-85; ieta <= 85; ieta++){
523			if( ieta==0) continue;
524			for(int iphi=1; iphi <= 360; iphi++){
525			inputcc >> c1;
526
527			int eta = ieta;
528			int phi = iphi;
529			convxtalid(phi,eta);
530
531			if( c1 <=0){
532			c1 =1;
533			}
534
535			if( c1 == 1 ){
536			c1 = - 1;
537			}
538			if( c1 <0.6 \|\| c1 >2.5){
539			c1 = -1;
540			}
541
542
543			CphiCorr[eta+85][phi] = c1;
544			}
545			}
546
547
548
549			}
550
551
552
553
554			void readInterCalibBarre_Phiv1(){
555
556			char *file_input = new char[500];
557			sprintf(file_input,"phisym_EcalIntercalib-160M-corrected.txt");
558
559
560			cout<<"READING fro inter-calib file "<<file_input<<endl;
561			ifstream inputcc(file_input,ios::in);
562			if (inputcc.fail()){
563			cout<<"error open file barrel.. " << file_input<<endl;
564			exit(1);
565			}
566
567			int eta;
568			int phi;
569			float c = 1;
570			float c1 = 1;
571
572			int n = 0;
573
574			eta = -85;
575			phi = 0;
576
577			float mean = 0;
578
579			int cureta = -85;
580
581			int ngood = 0;
582
583
584
585			while (inputcc.good()){
586
587			inputcc >> c1;
588
589			if(c1==1) c1 = -1;
590
591
592			if( n% 360 ==0 && n>0 ) {
593			eta++;
594			if(eta==0) eta = 1;
595			}
596
597			if( eta == cureta){
598			phi ++;
599			}else{
600			cureta = eta;
601			phi = 1;
602
603			}
604
605
606			n++;
607
608
609			int iSM = (phi-1)/20+1;
610			if( eta<0) iSM += 18;
611
612
613			int ieta = eta;
614			int iphi = phi;
615
616
617			convxtalid(iphi,ieta);
618
619			///txtout<<eta<<" "<<phi<<endl;
620
621			CphiCorr[ieta+85][iphi] = c1;
622
623
624			if(c1 >0 ){
625			mean += c1;
626			ngood ++;
627			}
628
629
630			if(n >= 360*170) break;
631
632			}
633
634
635			mean /= ngood;
636			cout<<n<<" phi constants read. "<< mean <<endl;
637
638
639			// for(int j=0; j< 170; j++){
640			// mean_ieta[j] /= nmean_ieta[j];
641			// txtout<<"mean_ieta_phisym "<< j<<" "<< mean_ieta[j] <<" "<< nmean_ieta[j] <<endl;
642			// }
643
644
645			// for(int j=0; j< 170; j++){
646			// for(int k=0; k< 360; k++){
647
648			// if( CphiCorr[j][k] >0 ){
649			// CphiCorr[j][k] /= mean_ieta[j] ;
650
651			// }
652			// Cphi[j][k] = CphiCorr[j][k] * iCalConst_data[j][k];
653			// }
654			// }
655
656			// NormSMScaleToUnit(CphiCorr);
657
658			}
659
660
661			///read initi calibration constants for each crystal
662			void readInterCalibBarre_BSv1(){
663
664
665			char *file_input = new char[500];
666			sprintf(file_input,"IC_barrel_splash09_AllRuns.txt");
667
668
669			cout<<"READING fro inter-calib bsfile "<<file_input<<endl;
670			txtout<<"READING fro inter-calib bsfile "<<file_input<<endl;
671			ifstream inputcc(file_input,ios::in);
672			if (inputcc.fail()){
673			cout<<"error open file barrel.. " << file_input<<endl;
674			exit(1);
675			}
676
677			int eta;
678			int phi;
679			float c = 1;
680			float c1 = 1;
681
682			int n = 0;
683			int ngood = 0;
684
685
686			float mean = 0;
687
688
689
690
691			float mean_ieta[170] = {0};
692			int nmean_ieta[170] = {0};
693
694
695			while (inputcc.good()){
696
697			inputcc >> eta >> phi >>c1;
698
699			// cout<<phi<<" "<<eta<<" "<<c<<endl;
700
701			convxtalid(phi,eta);
702
703			eta += 85;
704
705			if( c1 == 1){
706			c1 = -1;
707			}
708
709			if( c1 <0.6 \|\| c1 >2.5){
710			c1 = -1;
711			}
712
713
714			///iCalConstPi0[eta][phi] = c1 ;
715
716			CBSCorr[eta][phi] = c1;
717
718			if(c1 >0 ){
719			mean += c1;
720			ngood ++;
721			}
722
723			if(c1 >0 ){
724			mean_ieta[eta] += c1;
725			nmean_ieta[eta] ++;
726			}
727
728			n++;
729
730			/// hh_imcEB->Fill(c);
731			if(n >= 360*170) break;
732
733			}
734
735
736
737			mean /= ngood;
738
739
740
741			cout<<n<<" bs constants read. "<< mean << endl;
742
743
744
745			// for(int j=0; j< 170; j++){
746			// for(int k=0; k< 360; k++){
747
748			// if( CBSCorr[j][k] >0 ){
749			// CBSCorr[j][k] /= mean_ieta[j];
750			// }
751
752			// CBS[j][k] = CBSCorr[j][k] * iCalConst_data[j][k];
753
754
755			// }
756			// }
757			// NormSMScaleToUnit(CBSCorr);
758
759
760
761			}
762
763
764
765
766
767			///read initi calibration constants for each crystal
768			void readInterCalibBarre_GR09_V8(){
769
770
771
772			char *file_input = new char[500];
773			//file_input ="/data/yangyong/data/p
774			file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_R_V6A_barrel.txt";
775
776			///sprintf(file_input,"/uscms_data/d2/yongy/calib/backup/interCalib_STARTUP3X_V8P_barrel_mc.txt"); ///now the same used in data
777
778			//sprintf(file_input,"/uscms/home/yongy/work/interCalib_STARTUP3X_V8P_barrel.txt");
779
780			// sprintf(file_input,"/uscms/home/yongy/work/interCalib_GR_R_36X_V11A_Barrel.txt");
781
782			cout<<"READING fro inter-calib file "<<file_input<<endl;
783			ifstream inputcc(file_input,ios::in);
784			if (inputcc.fail()){
785			cout<<"error open file barrel.. " << file_input<<endl;
786			exit(1);
787			}
788
789			int eta;
790			int phi;
791			float c = 1;
792			float c1 = 1;
793
794			int n = 0;
795
796			double mean = 0;
797
798
799			while (inputcc.good()){
800
801			inputcc >> eta >> phi >>c1;
802
803			// cout<<phi<<" "<<eta<<" "<<c<<endl;
804
805			convxtalid(phi,eta);
806
807			eta += 85;
808
809			///iCalConst[eta][phi] = c;
810
811			//iCalConst_data[eta][phi] = c1 * 0.038012/ 0.038; ///Now this is exactly the same as the inter-calibration used in data.
812
813			interCalib_preCalib[eta][phi] = c1 ; /// now OK iN MC startupv8p.
814
815
816			n++;
817			// cout<<phi<<" "<<eta<<" "<<c<<endl;
818			mean += c1;
819
820			/// hh_imcEB->Fill(c);
821			if(n >= 360*170) break;
822
823			}
824
825
826			mean /= 61200;
827
828			cout<<n<<" constants read. "<< mean <<endl;
829
830
831			}
832
833
834
835
836
837
838
839
840			///read initi calibration constants for each crystal
841			void readInterCalibEndcap_GR09_V8(){
842			char *file_input = new char[500];
843			//file_input ="/data/yangyong/data/p
844			///file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_R_V6A_barrel.txt";
845
846			file_input = "/afs/cern.ch/cms/cit/yongy/interCalib_GR09_H_V6OFF_endcap.txt";
847
848			cout<<"READING fro inter-calib file "<<file_input<<endl;
849			ifstream inputcc(file_input,ios::in);
850			if (inputcc.fail()){
851			cout<<"error open file barrel.. " << file_input<<endl;
852			exit(1);
853			}
854
855			int eta;
856			int phi;
857			int iz;
858			float c = 1;
859			int n = 0;
860
861
862			while (inputcc.good()){
863
864			inputcc >> iz >> eta >> phi >> c;
865
866			int izz = iz < 0? 0:1;
867
868			interCalibEndcap_preCalib[izz][eta][phi] = c;
869			validRecHitEndCap[izz][eta][phi] = 1;
870
871			n ++;
872
873			}
874
875			cout<<n<<" constants read. "<<endl;
876
877			}
878
879
880			void getInterCalibEndcapv1(const char *inputfile,float C[2][101][101]){
881
882
883
884			ifstream inputcc(inputfile);
885			if (inputcc.fail()) {
886			cout<<"error "<< inputfile<<" can not be opened."<<endl;
887			exit(1);
888			}
889			int eta;
890			int phi;
891			int iz;
892			float c = 1;
893			int n = 0;
894
895			double mean = 0;
896
897			while (inputcc.good()){
898
899			inputcc >> iz >> eta >> phi >> c;
900
901			int izz = iz < 0? 0:1;
902
903			C[izz][eta][phi] = c;
904			mean += c;
905
906			n ++;
907
908			if( n>= 14648){
909			break;
910			}
911
912			}
913
914			cout<<n<<" constants read. "<< mean /n <<endl;
915
916			}
917
918
919
920			void getInterCalibEndcap(const char *inputfile, float C[2][101][101]){
921
922
923			ifstream inputcc(inputfile);
924			if (inputcc.fail()) {
925			cout<<"error "<< inputfile<<" can not be opened."<<endl;
926			exit(1);
927			}
928
929			int ix,iy,iz;
930			float cc;
931
932			float ccMeaniring[50] ;
933			int nccMeaniring[50];
934
935
936			// for(int j=0; j<2;j++){
937			for(int k=0; k< kEndcEtaRings; k++){
938			ccMeaniring[k] = 0;
939			nccMeaniring[k] = 0;
940			}
941			// }
942
943			int nread = 0;
944			while( inputcc.good()){
945			inputcc >> iz >> ix >> iy >> cc;
946			C[iz][ix][iy] = cc;
947			int iring = iRingEndCap(2*iz-1,ix,iy); ///input -1/1,
948			if( cc > 0){
949			ccMeaniring[iring] += cc;
950			nccMeaniring[iring] += 1;
951			}
952
953			nread ++;
954			if( iz==1 && ix == 100 && iy == 60) break;
955			}
956
957			cout<<"nb of crystals endcap " << nread <<endl;
958
959
960			//for(int j=0; j<2;j++){
961			for(int k=0; k< kEndcEtaRings; k++){
962			ccMeaniring[k] /= nccMeaniring[k];
963			// }
964			}
965
966			double meanall = 0;
967			int nall = 0;
968			for(int iz=0; iz<2; iz++){
969			for(int j=0; j<101; j++){
970			for(int k=0; k< 101; k++){
971			if( validRecHitEndCap[iz][j][k] <1) continue;
972			int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
973			if( C[iz][j][k] >0) {
974			C[iz][j][k] /= ccMeaniring[iring];
975			meanall += C[iz][j][k];
976			nall ++;
977			}
978
979			}
980			}
981			}
982
983			meanall /= nall;
984			cout<<" corrfactoriZiXiY read " << meanall <<" "<< nall <<endl;
985
986
987			}
988
989
990
991			void NormCrystalDeadFlagEndcap_v1(float C[2][101][101],int ndeadflag[2][101][101]){
992
993			float mean_deadflag[20] = {0};
994			float nmean_deadflag[20] = {0};
995
996			int ndeadCrystals = 0;
997
998
999			for(int iz =0; iz < 2; iz++){
1000			for(int j=0; j< 101; j++){
1001			for(int k=0; k<101; k++){
1002
1003			if( validRecHitEndCap[iz][j][k] ==0) continue;
1004			if( C[iz][j][k] <0){
1005			continue;
1006			}
1007			int ndead = ndeadflag[iz][j][k];
1008			if( ndead >=0){
1009			mean_deadflag[ndead] += C[iz][j][k];
1010			nmean_deadflag[ndead] ++;
1011			}
1012
1013			}
1014			}
1015			}
1016			cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
1017			for(int j=0; j<20; j++){
1018			if( nmean_deadflag[j] >0){
1019			mean_deadflag[j] /= nmean_deadflag[j];
1020			cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
1021			}
1022			}
1023
1024
1025			for(int iz =0; iz < 2; iz++){
1026			for(int j=0; j< 101; j++){
1027			for(int k=0; k<101; k++){
1028
1029			if( validRecHitEndCap[iz][j][k] ==0) continue;
1030			if( C[iz][j][k] <0){
1031			continue;
1032			}
1033
1034			int ndead = ndeadflag[iz][j][k];
1035			if( ndead >=0){
1036			C[iz][j][k] /= mean_deadflag[ndead];
1037			}
1038
1039			}
1040			}
1041			}
1042
1043
1044
1045			}
1046
1047
1048
1049
1050
1051
1052			void NormCrystalDeadFlagEndcap(float C[2][101][101]){
1053
1054			float mean_deadflag[20] = {0};
1055			float nmean_deadflag[20] = {0};
1056
1057			int ndeadCrystals = 0;
1058
1059
1060			for(int iz =0; iz < 2; iz++){
1061			for(int j=0; j< 101; j++){
1062			for(int k=0; k<101; k++){
1063
1064			if( validRecHitEndCap[iz][j][k] ==0) continue;
1065			if( C[iz][j][k] <0){
1066			continue;
1067			}
1068			int ndeadflag = ndeadflag_endcap[iz][j][k];
1069			if( ndeadflag >=0){
1070			mean_deadflag[ndeadflag] += C[iz][j][k];
1071			nmean_deadflag[ndeadflag] ++;
1072			}
1073
1074			}
1075			}
1076			}
1077			cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
1078			for(int j=0; j<20; j++){
1079			if( nmean_deadflag[j] >0){
1080			mean_deadflag[j] /= nmean_deadflag[j];
1081			cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
1082			}
1083			}
1084
1085
1086			for(int iz =0; iz < 2; iz++){
1087			for(int j=0; j< 101; j++){
1088			for(int k=0; k<101; k++){
1089
1090			if( validRecHitEndCap[iz][j][k] ==0) continue;
1091			if( C[iz][j][k] <0){
1092			continue;
1093			}
1094
1095			int ndeadflag = ndeadflag_endcap[iz][j][k];
1096			if( ndeadflag >=0){
1097			C[iz][j][k] /= mean_deadflag[ndeadflag];
1098			}
1099
1100			}
1101			}
1102			}
1103
1104
1105
1106			}
1107
1108
1109
1110			void scaleMeanToUnitEndcap(float C[2][101][101]){
1111
1112
1113			double mean = 0;
1114			int n = 0;
1115
1116
1117			for(int iz =0; iz < 2; iz++){
1118			for(int j=0; j< 101; j++){
1119			for(int k=0; k<101; k++){
1120			if( validRecHitEndCap[iz][j][k] ==0) continue;
1121			if( C[iz][j][k] <0){
1122			continue;
1123			}
1124			mean += C[iz][j][k];
1125			n ++;
1126			}
1127			}
1128			}
1129
1130			mean /= n;
1131			cout<<"mean " << mean<<" "<<n<<endl;
1132
1133
1134			for(int iz =0; iz < 2; iz++){
1135			for(int j=0; j< 101; j++){
1136			for(int k=0; k<101; k++){
1137			if( validRecHitEndCap[iz][j][k] ==0) continue;
1138			if( C[iz][j][k] <0){
1139			continue;
1140			}
1141
1142			C[iz][j][k] /= mean;
1143			}
1144			}
1145			}
1146
1147			}