HcalCalibAlgos/macros/get_correction_8.9mln.C

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 }