yangyong/combineICEE/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 = "interCalib_GR09_R_V6A_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 = "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};
  int 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 && C[iz][j][k] >0.7 && C[iz][j][k] <1.3 ){
          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.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
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 = "interCalib_GR09_R_V6A_barrel.txt";
774			cout<<"READING fro inter-calib file "<<file_input<<endl;
775			ifstream inputcc(file_input,ios::in);
776			if (inputcc.fail()){
777			cout<<"error open file barrel.. " << file_input<<endl;
778			exit(1);
779			}
780
781			int eta;
782			int phi;
783			float c = 1;
784			float c1 = 1;
785
786			int n = 0;
787
788			double mean = 0;
789
790
791			while (inputcc.good()){
792
793			inputcc >> eta >> phi >>c1;
794
795			// cout<<phi<<" "<<eta<<" "<<c<<endl;
796
797			convxtalid(phi,eta);
798
799			eta += 85;
800
801			///iCalConst[eta][phi] = c;
802
803			//iCalConst_data[eta][phi] = c1 * 0.038012/ 0.038; ///Now this is exactly the same as the inter-calibration used in data.
804
805			interCalib_preCalib[eta][phi] = c1 ; /// now OK iN MC startupv8p.
806
807
808			n++;
809			// cout<<phi<<" "<<eta<<" "<<c<<endl;
810			mean += c1;
811
812			/// hh_imcEB->Fill(c);
813			if(n >= 360*170) break;
814
815			}
816
817
818			mean /= 61200;
819
820			cout<<n<<" constants read. "<< mean <<endl;
821
822
823			}
824
825
826
827
828
829
830
831
832			///read initi calibration constants for each crystal
833			void readInterCalibEndcap_GR09_V8(){
834			char *file_input = new char[500];
835			file_input = "interCalib_GR09_H_V6OFF_endcap.txt";
836
837			cout<<"READING fro inter-calib file "<<file_input<<endl;
838			ifstream inputcc(file_input,ios::in);
839			if (inputcc.fail()){
840			cout<<"error open file barrel.. " << file_input<<endl;
841			exit(1);
842			}
843
844			int eta;
845			int phi;
846			int iz;
847			float c = 1;
848			int n = 0;
849
850
851			while (inputcc.good()){
852
853			inputcc >> iz >> eta >> phi >> c;
854
855			int izz = iz < 0? 0:1;
856
857			interCalibEndcap_preCalib[izz][eta][phi] = c;
858			validRecHitEndCap[izz][eta][phi] = 1;
859
860			n ++;
861
862			}
863
864			cout<<n<<" constants read. "<<endl;
865
866			}
867
868
869			void getInterCalibEndcapv1(const char *inputfile,float C[2][101][101]){
870
871
872
873			ifstream inputcc(inputfile);
874			if (inputcc.fail()) {
875			cout<<"error "<< inputfile<<" can not be opened."<<endl;
876			exit(1);
877			}
878			int eta;
879			int phi;
880			int iz;
881			float c = 1;
882			int n = 0;
883
884			double mean = 0;
885
886			while (inputcc.good()){
887
888			inputcc >> iz >> eta >> phi >> c;
889
890			int izz = iz < 0? 0:1;
891
892			C[izz][eta][phi] = c;
893			mean += c;
894
895			n ++;
896
897			if( n>= 14648){
898			break;
899			}
900
901			}
902
903			cout<<n<<" constants read. "<< mean /n <<endl;
904
905			}
906
907
908
909			void getInterCalibEndcap(const char *inputfile, float C[2][101][101]){
910
911
912			ifstream inputcc(inputfile);
913			if (inputcc.fail()) {
914			cout<<"error "<< inputfile<<" can not be opened."<<endl;
915			exit(1);
916			}
917
918			int ix,iy,iz;
919			float cc;
920
921			float ccMeaniring[50] ;
922			int nccMeaniring[50];
923
924
925			// for(int j=0; j<2;j++){
926			for(int k=0; k< kEndcEtaRings; k++){
927			ccMeaniring[k] = 0;
928			nccMeaniring[k] = 0;
929			}
930			// }
931
932			int nread = 0;
933			while( inputcc.good()){
934			inputcc >> iz >> ix >> iy >> cc;
935			C[iz][ix][iy] = cc;
936			int iring = iRingEndCap(2*iz-1,ix,iy); ///input -1/1,
937			if( cc > 0){
938			ccMeaniring[iring] += cc;
939			nccMeaniring[iring] += 1;
940			}
941
942			nread ++;
943			if( iz==1 && ix == 100 && iy == 60) break;
944			}
945
946			cout<<"nb of crystals endcap " << nread <<endl;
947
948
949			//for(int j=0; j<2;j++){
950			for(int k=0; k< kEndcEtaRings; k++){
951			ccMeaniring[k] /= nccMeaniring[k];
952			// }
953			}
954
955			double meanall = 0;
956			int nall = 0;
957			for(int iz=0; iz<2; iz++){
958			for(int j=0; j<101; j++){
959			for(int k=0; k< 101; k++){
960			if( validRecHitEndCap[iz][j][k] <1) continue;
961			int iring = iRingEndCap(2*iz-1,j,k); ///input -1/1,
962			if( C[iz][j][k] >0) {
963			C[iz][j][k] /= ccMeaniring[iring];
964			meanall += C[iz][j][k];
965			nall ++;
966			}
967
968			}
969			}
970			}
971
972			meanall /= nall;
973			cout<<" corrfactoriZiXiY read " << meanall <<" "<< nall <<endl;
974
975
976			}
977
978
979
980			void NormCrystalDeadFlagEndcap_v1(float C[2][101][101],int ndeadflag[2][101][101]){
981
982			float mean_deadflag[20] = {0};
983			int nmean_deadflag[20] = {0};
984
985			int ndeadCrystals = 0;
986
987
988			for(int iz =0; iz < 2; iz++){
989			for(int j=0; j< 101; j++){
990			for(int k=0; k<101; k++){
991
992			if( validRecHitEndCap[iz][j][k] ==0) continue;
993			if( C[iz][j][k] <0){
994			continue;
995			}
996			int ndead = ndeadflag[iz][j][k];
997			if( ndead >=0 && C[iz][j][k] >0.7 && C[iz][j][k] <1.3 ){
998			mean_deadflag[ndead] += C[iz][j][k];
999			nmean_deadflag[ndead] ++;
1000			}
1001
1002			}
1003			}
1004			}
1005			cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
1006			for(int j=0; j<20; j++){
1007			if( nmean_deadflag[j] >0){
1008			mean_deadflag[j] /= nmean_deadflag[j];
1009			cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
1010			}
1011			}
1012
1013
1014			for(int iz =0; iz < 2; iz++){
1015			for(int j=0; j< 101; j++){
1016			for(int k=0; k<101; k++){
1017
1018			if( validRecHitEndCap[iz][j][k] ==0) continue;
1019			if( C[iz][j][k] <0){
1020			continue;
1021			}
1022
1023			int ndead = ndeadflag[iz][j][k];
1024			if( ndead >=0){
1025			C[iz][j][k] /= mean_deadflag[ndead];
1026			}
1027
1028			}
1029			}
1030			}
1031
1032
1033
1034			}
1035
1036
1037
1038
1039
1040
1041			void NormCrystalDeadFlagEndcap(float C[2][101][101]){
1042
1043			float mean_deadflag[20] = {0};
1044			float nmean_deadflag[20] = {0};
1045
1046			int ndeadCrystals = 0;
1047
1048
1049			for(int iz =0; iz < 2; iz++){
1050			for(int j=0; j< 101; j++){
1051			for(int k=0; k<101; k++){
1052
1053			if( validRecHitEndCap[iz][j][k] ==0) continue;
1054			if( C[iz][j][k] <0){
1055			continue;
1056			}
1057			int ndeadflag = ndeadflag_endcap[iz][j][k];
1058			if( ndeadflag >=0){
1059			mean_deadflag[ndeadflag] += C[iz][j][k];
1060			nmean_deadflag[ndeadflag] ++;
1061			}
1062
1063			}
1064			}
1065			}
1066			cout<<" NormCrystalDeadFlag ndeadCrystals: " << ndeadCrystals <<endl;
1067			for(int j=0; j<20; j++){
1068			if( nmean_deadflag[j] >0){
1069			mean_deadflag[j] /= nmean_deadflag[j];
1070			cout<<" mean_deadflag "<<j<<" "<< mean_deadflag[j]<< " "<< nmean_deadflag[j] <<endl;
1071			}
1072			}
1073
1074
1075			for(int iz =0; iz < 2; iz++){
1076			for(int j=0; j< 101; j++){
1077			for(int k=0; k<101; k++){
1078
1079			if( validRecHitEndCap[iz][j][k] ==0) continue;
1080			if( C[iz][j][k] <0){
1081			continue;
1082			}
1083
1084			int ndeadflag = ndeadflag_endcap[iz][j][k];
1085			if( ndeadflag >=0){
1086			C[iz][j][k] /= mean_deadflag[ndeadflag];
1087			}
1088
1089			}
1090			}
1091			}
1092
1093
1094
1095			}
1096
1097
1098
1099			void scaleMeanToUnitEndcap(float C[2][101][101]){
1100
1101
1102			double mean = 0;
1103			int n = 0;
1104
1105
1106			for(int iz =0; iz < 2; iz++){
1107			for(int j=0; j< 101; j++){
1108			for(int k=0; k<101; k++){
1109			if( validRecHitEndCap[iz][j][k] ==0) continue;
1110			if( C[iz][j][k] <0){
1111			continue;
1112			}
1113			mean += C[iz][j][k];
1114			n ++;
1115			}
1116			}
1117			}
1118
1119			mean /= n;
1120			cout<<"mean " << mean<<" "<<n<<endl;
1121
1122
1123			for(int iz =0; iz < 2; iz++){
1124			for(int j=0; j< 101; j++){
1125			for(int k=0; k<101; k++){
1126			if( validRecHitEndCap[iz][j][k] ==0) continue;
1127			if( C[iz][j][k] <0){
1128			continue;
1129			}
1130
1131			C[iz][j][k] /= mean;
1132			}
1133			}
1134			}
1135
1136			}