Fork me on GitHub

Changeset 71 in svn for trunk/src


Ignore:
Timestamp:
Dec 3, 2008, 2:47:52 AM (16 years ago)
Author:
Xavier Rouby
Message:

iEta et iPhi. Verification non complete.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SmearUtil.cc

    r69 r71  
    5252ELG_Ncen =        0.25 ;                 // N term for central ECAL
    5353ELG_Ccen =        0.0055 ;                // C term for central ECAL
    54 ELG_Sfwd =        1.5;                    // S term for forward ECAL
    55 ELG_Nfwd =        0.0;                    // N term for central ECAL
    56 ELG_Cfwd =        0.06;                   // C term for forward ECAL
     54ELG_Cfwd =        0.107 ;                 // S term for forward ECAL
     55ELG_Sfwd =        2.084 ;                 // C term for forward ECAL
     56ELG_Nfwd =        0.0   ;                 // N term for central ECAL
    5757
    5858HAD_Shcal =       1.5 ;                    // S term for central HCAL // hadronic calorimeter
    … …  
    104104M_MAXITERATIONS     =   100;
    105105
     106// Define Calorimeter Towers
     107NTOWERS = 40;
     108
     109const float tower_eta_edges[41] = {
     110    0., 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566,
     111    1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.868, 2.950, 3.125, 3.300, 3.475, 3.650, 3.825, 4.000, 4.175,
     112    4.350, 4.525,  4.700, 5.000}; // temporary object
     113TOWER_ETA_EDGES = new float[NTOWERS+1];
     114for(unsigned int i=0; i<NTOWERS+1; i++) TOWER_ETA_EDGES[i] = tower_eta_edges[i];
     115
     116const float tower_dphi[40] = {
     117     5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 10,
     118    10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; // temporary object
     119TOWER_DPHI = new float[NTOWERS];
     120for(unsigned int i=0; i<NTOWERS; i++) TOWER_DPHI[i] = tower_dphi[i];
     121
    106122}
    107123
    … …  
    122138    curstring.str(temp_string);
    123139    string varname;
    124     float value;
     140    float value; int ivalue;
    125141   
    126142    if(strstr(temp_string.c_str(),"#")) { }
    … …  
    129145    else if(strstr(temp_string.c_str(),"MAX_CALO_FWD")){curstring >> varname >> value; MAX_CALO_FWD = value;}
    130146    else if(strstr(temp_string.c_str(),"MAX_MU")){curstring >> varname >> value; MAX_MU = value;}
    131     else if(strstr(temp_string.c_str(),"TRACKING_RADIUS")){curstring >> varname >> value; TRACKING_RADIUS = (int)value;}
    132     else if(strstr(temp_string.c_str(),"TRACKING_LENGTH")){curstring >> varname >> value; TRACKING_LENGTH = (int)value;}
     147    else if(strstr(temp_string.c_str(),"TRACKING_RADIUS")){curstring >> varname >> ivalue; TRACKING_RADIUS = ivalue;}
     148    else if(strstr(temp_string.c_str(),"TRACKING_LENGTH")){curstring >> varname >> ivalue; TRACKING_LENGTH = ivalue;}
    133149    else if(strstr(temp_string.c_str(),"BFIELD_X")){curstring >> varname >> value; BFIELD_X = value;}
    134150    else if(strstr(temp_string.c_str(),"BFIELD_Y")){curstring >> varname >> value; BFIELD_Y = value;}
    … …  
    151167    else if(strstr(temp_string.c_str(),"PT_TRACK_TAU")){curstring >> varname >> value; PT_TRACK_TAU = value;}
    152168    else if(strstr(temp_string.c_str(),"PT_TRACKS_MIN")){curstring >> varname >> value; PT_TRACKS_MIN = value;}
    153     else if(strstr(temp_string.c_str(),"TAGGING_B")){curstring >> varname >> value; TAGGING_B = (int)value;}
    154     else if(strstr(temp_string.c_str(),"MISTAGGING_C")){curstring >> varname >> value; MISTAGGING_C = (int)value;}
    155     else if(strstr(temp_string.c_str(),"MISTAGGING_L")){curstring >> varname >> value; MISTAGGING_L = (int)value;}
     169    else if(strstr(temp_string.c_str(),"TAGGING_B")){curstring >> varname >> ivalue; TAGGING_B = ivalue;}
     170    else if(strstr(temp_string.c_str(),"MISTAGGING_C")){curstring >> varname >> ivalue; MISTAGGING_C = ivalue;}
     171    else if(strstr(temp_string.c_str(),"MISTAGGING_L")){curstring >> varname >> ivalue; MISTAGGING_L = ivalue;}
    156172    else if(strstr(temp_string.c_str(),"CONERADIUS")){curstring >> varname >> value; CONERADIUS = value;}
    157     else if(strstr(temp_string.c_str(),"JETALGO")){curstring >> varname >> value; JETALGO = (int)value;}
    158     else if(strstr(temp_string.c_str(),"TRACKING_EFF")){curstring >> varname >> value; TRACKING_EFF = (int)value;}
     173    else if(strstr(temp_string.c_str(),"JETALGO")){curstring >> varname >> ivalue; JETALGO = ivalue;}
     174    else if(strstr(temp_string.c_str(),"TRACKING_EFF")){curstring >> varname >> ivalue; TRACKING_EFF = ivalue;}
    159175    else if(strstr(temp_string.c_str(),"ELEC_pt")){curstring >> varname >> value; ELEC_pt = value;}
    160176    else if(strstr(temp_string.c_str(),"MUON_pt")){curstring >> varname >> value; MUON_pt = value;}
    161177    else if(strstr(temp_string.c_str(),"JET_pt")){curstring >> varname >> value; JET_pt = value;}
    162178    else if(strstr(temp_string.c_str(),"TAUJET_pt")){curstring >> varname >> value; TAUJET_pt = value;}
     179    else if(strstr(temp_string.c_str(),"NTOWERS")){curstring >> varname >> ivalue; NTOWERS = ivalue;}
     180    else if(strstr(temp_string.c_str(),"TOWER_ETA_EDGES")){
     181        curstring >> varname;   for(unsigned int i=0; i<NTOWERS+1; i++) {curstring >> value; TOWER_ETA_EDGES[i] = value;} }
     182    else if(strstr(temp_string.c_str(),"TOWER_DPHI")){
     183        curstring >> varname;   for(unsigned int i=0; i<NTOWERS; i++) {curstring >> value; TOWER_DPHI[i] = value;} }
     184   
    163185
    164186  }
    … …  
    182204 PROTOJET_PTMIN = 0.0;
    183205
    184  DOTRIGGER=1;
    185206
    186207}
    187208
    188209void RESOLution::Logfile(string LogName) {
     210//void RESOLution::Logfile(string outputfilename) {
    189211
    190212  ofstream f_out(LogName.c_str());
    … …  
    266288        << left << setw(5) <<BFIELD_Z <<""<< right << setw(10)<<"*"<<"\n";
    267289  f_out<<"*                                                                    *"<<"\n";
     290
     291
     292  f_out<<"*                                                                    *"<<"\n";
     293  f_out<<"#********************                                                *"<<"\n";
     294  f_out<<"# Calorimetric Towers                                                *"<<"\n";
     295  f_out<<"#********************                                                *"<<"\n";
     296  f_out << left << setw(55) <<"* Number of calorimetric towers in eta, for eta>0: "<<""
     297        << left << setw(5) << NTOWERS <<""<< right << setw(10)<<"*"<<"\n";
     298  f_out << left << setw(55) <<"* Tower edges in eta, for eta>0: "<<"" << right << setw(15)<<"*"<<"\n";
     299          f_out << "*   ";
     300          for (unsigned int i=0; i<NTOWERS+1; i++) {
     301                f_out << left << setw(7) << TOWER_ETA_EDGES[i];
     302                if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "*   ";
     303          }
     304          for (unsigned int i=(NTOWERS+1)%9; i<9; i++) f_out << left << setw(7) << "";
     305          f_out << right << setw(3)<<"*"<<"\n";
     306  f_out << left << setw(55) <<"* Tower sizes in phi, for eta>0 [degree]:"<<"" << right << setw(15)<<"*"<<"\n";
     307          f_out << "*   ";
     308          for (unsigned int i=0; i<NTOWERS; i++) {
     309                f_out << left << setw(7) << TOWER_DPHI[i];
     310                if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "*   ";
     311          }
     312          for (unsigned int i=(NTOWERS)%9; i<9; i++) f_out << left << setw(7) << "";
     313          f_out << right << setw(3)<<"*"<<"\n";
     314  f_out<<"*                                                                    *"<<"\n";
     315
    268316  f_out<<"#************************************                                *"<<"\n";
    269317  f_out<<"# Electromagnetic smearing parameters                                *"<<"\n";
    … …  
    446494  float energy = electron.E();  // before smearing
    447495  float energyS = 0.0;          // after smearing  // \sigma/E = C + N/E + S/\sqrt{E}
    448   float eta=fabs(electron.Eta());
    449  
    450   if(eta < MAX_TRACKER) { // if the electron is inside the tracker
     496 
     497  if(fabs(electron.Eta()) < MAX_TRACKER) { // if the electron is inside the tracker
    451498    energyS = gRandom->Gaus(energy, sqrt(
    452499                                         pow(ELG_Ncen,2) +
    … …  
    454501                                         pow(ELG_Scen*sqrt(energy),2) ));
    455502  }
    456   if(eta > MAX_TRACKER && eta < MAX_CALO_FWD){
     503  if(fabs(electron.Eta()) > MAX_TRACKER && fabs(electron.Eta()) < MAX_CALO_FWD){
    457504    energyS = gRandom->Gaus(energy, sqrt(
    458505                                         pow(ELG_Nfwd,2) +
    … …  
    609656
    610657
     658  //********** returns a segmented value for eta and phi, for calo towers *****
     659void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
     660   iEta = -100;
     661   int index=-100;
     662   for (unsigned int i=1; i< NTOWERS+1; i++) {
     663        if(fabs(eta)>TOWER_ETA_EDGES[i-1] && fabs(eta)<TOWER_ETA_EDGES[i]) {
     664                iEta = (eta>0) ? TOWER_ETA_EDGES[i-1] : -TOWER_ETA_EDGES[i];
     665                index = i-1;
     666                //cout << setw(15) << left << eta << "\t" << iEta << endl;
     667                break;
     668        }
     669   }
     670   if(index==-100) return;
     671   iPhi = -100;
     672   float dphi = TOWER_DPHI[index]*PI/180.;
     673   for (unsigned int i=1; i < 360/TOWER_DPHI[index]; i++ ) {
     674        float low = -PI+(i-1)*dphi;
     675        float high= low+dphi;
     676        if(phi > low && phi < high ){
     677                iPhi = low;
     678                break;
     679        }
     680   }
     681   if (phi > PI-dphi) iPhi = PI-dphi;
     682}
     683
     684
    611685//**************************** Returns the delta Phi ****************************
    612686float DeltaPhi(const float phi1, const float phi2) {
Note: See TracChangeset for help on using the changeset viewer.