Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:34:31

0001 {
0002 gROOT->Reset();
0003 gROOT->SetStyle("Plain");
0004 
0005 gStyle->SetOptStat(1111);
0006 gStyle->SetOptFit(111);
0007 
0008    Float_t plsignal[50][73][5][5],plnoise[50][73][5][5];
0009    Float_t minsignal[50][73][5][5],minnoise[50][73][5][5];
0010    Float_t excluded_min[50][73][5][5];
0011    
0012    FILE *Out1 = fopen("coefficients_8.9mln.txt", "w+");  
0013    
0014     for(Int_t k=1; k<50; k++) 
0015     {
0016     for(Int_t i=1; i<73; i++) 
0017     {  
0018     for(Int_t j=1; j<5; j++) 
0019     {  
0020     for(Int_t l=1; l<5; l++) 
0021     {  
0022       excluded_min[k][i][j][l] = 0.;
0023     } // l  
0024     } // j
0025     } // i
0026     } // k 
0027    
0028    std::string line;
0029    
0030    std::ifstream in21( "HB_exclusion.txt" );
0031 
0032    while( std::getline( in21, line)){
0033    istringstream linestream(line);
0034 
0035    int eta,phi,dep;
0036    linestream>>dep>>eta>>phi;
0037     cout<<" Eta="<<eta<<endl;
0038    if( eta > 0 )
0039    {
0040      cout<<" eta > 0 "<<endl;
0041    } else
0042    {
0043      excluded_min[abs(eta)][phi][dep][1] = 1;
0044    }
0045   }
0046 
0047 
0048    cout<<" Read "<<endl;
0049    
0050    std::ifstream in20( "var_noise_8.9mln.txt" );
0051    
0052    while( std::getline( in20, line)){
0053    istringstream linestream(line);
0054    Float_t var,err;
0055    int subd,eta,phi,dep;
0056    linestream>>subd>>eta>>phi>>dep>>var;
0057    
0058    
0059    if( eta > 0 )
0060    {
0061      plnoise[eta][phi][dep][subd] = var;
0062    } else
0063    {
0064       minnoise[abs(eta)][phi][dep][subd] = var;
0065    }
0066    }
0067    cout<<" End of noise read "<<endl;
0068    
0069    std::ifstream in21( "var_minbias_8.9mln.txt" );
0070    
0071    while( std::getline( in21, line)){
0072    istringstream linestream(line);
0073    Float_t var,err;
0074    int subd,eta,phi,dep;
0075    linestream>>subd>>eta>>phi>>dep>>var>>err;
0076    if( eta > 0 )
0077    {
0078      plsignal[eta][phi][dep][subd] = var;
0079    } else
0080    {
0081       minsignal[abs(eta)][phi][dep][subd] = var;
0082    }
0083    }
0084    
0085    cout<<" End of signal read "<<endl;
0086 // 
0087 // Choose depth
0088 //   
0089    
0090    Float_t plmean[50][5][5]; 
0091    Float_t minmean[50][5][5];
0092    Int_t numchan[50][5][5];
0093    
0094    for(Int_t k=1; k<5; k++) 
0095    {
0096     for(Int_t idep0=1; idep0<5; idep0++) 
0097    {  
0098     for(Int_t i=1; i<42; i++) {
0099     
0100      plmean[i][idep0][k] = 0.; 
0101      minmean[i][idep0][k] = 0.;
0102 
0103      Int_t nch1 = 0; 
0104      Int_t nchp = 0; 
0105      Int_t nchm = 0;
0106        
0107      for(Int_t j=1; j<73; j++){
0108         if(minsignal[i][j][idep0][k]>0.) nch1++;
0109         if(plsignal[i][j][idep0][k]>0.) nchp++;
0110         if(minsignal[i][j][idep0][k]>0.) nchm++;
0111         plmean[i][idep0][k] = plmean[i][idep0][k] + plsignal[i][j][idep0][k] - plnoise[i][j][idep0][k];
0112         minmean[i][idep0][k] = minmean[i][idep0][k] + minsignal[i][j][idep0][k] - minnoise[i][j][idep0][k];
0113       } // j
0114           
0115     numchan[i][idep0][k] = nch1;    
0116     if(nch1 == 0) continue;      
0117     plmean[i][idep0][k] = plmean[i][idep0][k]/nchp;
0118     minmean[i][idep0][k] = minmean[i][idep0][k]/nchm;
0119     cout<<" k, idep0, i, nch1= "<<k<<" "<<idep0<<" "<<i<<" "<<nch1<<endl;
0120   // Do not calibrate HO
0121     Float_t err = 0.00001;
0122     
0123     if( k != 3 ) {
0124     for(Int_t j=1;j<73;j++) {
0125        if( plsignal[i][j][idep0][k] > 0.) {       
0126        Float_t tt0; 
0127        if( plsignal[i][j][idep0][k] - plnoise[i][j][idep0][k]>0.0000001 && plmean[i][idep0][k]>0.) {
0128          tt0 = sqrt(plmean[i][idep0][k]/(plsignal[i][j][idep0][k] - plnoise[i][j][idep0][k]));
0129        } else { 
0130      tt0 = 1.;
0131        }
0132          fprintf(Out1,"%d %d %d %d %.5f %.5f\n",k,idep0,i,j,tt0,err);       
0133        } // plnoise
0134     } // j
0135     for(Int_t j=1;j<73;j++) {   
0136        if( minsignal[i][j][idep0][k] > 0. ) {       
0137        Float_t tt0;        
0138        if(minsignal[i][j][idep0][k] - minnoise[i][j][idep0][k]>0.&& minmean[i][idep0][k]>0. ) {
0139           tt0 = sqrt(minmean[i][idep0][k]/(minsignal[i][j][idep0][k] - minnoise[i][j][idep0][k]));
0140         } else {  
0141           tt0 = 1.;
0142        }      
0143        Int_t ieta = -1*i;
0144        fprintf(Out1,"%d %d %d %d %.5f %.5f\n",k,idep0,ieta,j,tt0,err);
0145        } // minnoise
0146      } // j
0147    } // TMP  
0148     } else { 
0149         for(Int_t j=1;j<73;j++){   
0150            Float_t tt0=1.; 
0151            Int_t ieta = i;
0152            fprintf(Out1,"%d %d %d %d %.5f %.5f\n",k,idep0,ieta,j,tt0,err);
0153         }// j
0154          
0155         for(Int_t j=1;j<73;j++){   
0156            Float_t tt0=1.; 
0157            Int_t ieta = -1*i;
0158            fprintf(Out1,"%d %d %d %d %.5f %.5f\n",k,idep0,ieta,j,tt0,err);
0159         }// j
0160      } // k = 3
0161       cout<<" End "<<endl;
0162    } // i
0163    } // idep0
0164    } // k   
0165     
0166    fclose(Out1);
0167 }