Fork me on GitHub

Changeset d77b51d in git for modules/IdentificationMap.cc


Ignore:
Timestamp:
Sep 29, 2015, 2:08:10 PM (9 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
a98c7ef
Parents:
d870fc5 (diff), 06ec139 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge remote-tracking branch 'upstream/master'

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/IdentificationMap.cc

    rd870fc5 rd77b51d  
    6969void IdentificationMap::Init()
    7070{
    71   // read efficiency formula
    72 
    73 
    7471  TMisIDMap::iterator itEfficiencyMap;
    7572  ExRootConfParam param;
     
    8784    formula->Compile(param[i*3 + 2].GetString());
    8885    pdg = param[i*3].GetInt();
    89     fEfficiencyMap.insert(make_pair(pdg,make_pair(param[i*3 + 1].GetInt(),formula)));
    90 
    91    // cout<<param[i*3].GetInt()<<","<<param[i*3+1].GetInt()<<","<<param[i*3 + 2].GetString()<<endl;
    92 
     86    fEfficiencyMap.insert(make_pair(pdg, make_pair(param[i*3 + 1].GetInt(), formula)));
    9387  }
    9488
     
    10094    formula->Compile("1.0");
    10195
    102     fEfficiencyMap.insert(make_pair(0,make_pair(0,formula)));
     96    fEfficiencyMap.insert(make_pair(0, make_pair(0, formula)));
    10397  }
    10498
     
    126120    if(formula) delete formula;
    127121  }
    128 
    129122}
    130123
     
    134127{
    135128  Candidate *candidate;
    136   Double_t pt, eta, phi;
     129  Double_t pt, eta, phi, e;
    137130  TMisIDMap::iterator itEfficiencyMap;
    138131  pair <TMisIDMap::iterator, TMisIDMap::iterator> range;
    139132  DelphesFormula *formula;
    140   Int_t pdgIn, pdgOut, charge;
     133  Int_t pdgCodeIn, pdgCodeOut, charge;
    141134
    142   Double_t P, Pi;
    143 
    144  // cout<<"------------ New Event ------------"<<endl;
     135  Double_t p, r, total;
    145136
    146137  fItInputArray->Reset();
     
    152143    phi = candidatePosition.Phi();
    153144    pt = candidateMomentum.Pt();
    154     pdgIn = candidate->PID;
     145    e = candidateMomentum.E();
     146   
     147    pdgCodeIn = candidate->PID;
    155148    charge = candidate->Charge;
    156149
    157    // cout<<"------------ New Candidate ------------"<<endl;
    158    // cout<<candidate->PID<<"   "<<pt<<","<<eta<<","<<phi<<endl;
     150    // first check that PID of this particle is specified in the map
     151    // otherwise, look for PID = 0
    159152
    160     P = 1.0;
     153    itEfficiencyMap = fEfficiencyMap.find(pdgCodeIn);
    161154
    162     //first check that PID of this particle is specified in cfg, if not set look for PID=0
    163 
    164     itEfficiencyMap = fEfficiencyMap.find(pdgIn);
    165 
    166     range = fEfficiencyMap.equal_range(pdgIn);
    167     if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgIn);
     155    range = fEfficiencyMap.equal_range(pdgCodeIn);
     156    if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgCodeIn);
    168157    if(range.first == range.second) range = fEfficiencyMap.equal_range(0);
    169158
    170     //loop over submap for this pid
    171     for (TMisIDMap::iterator it=range.first; it!=range.second; ++it)
     159    r = gRandom->Uniform();
     160    total = 0.0;
     161
     162    // loop over sub-map for this PID
     163    for(TMisIDMap::iterator it = range.first; it != range.second; ++it)
    172164    {
     165      formula = (it->second).second;
     166      pdgCodeOut = (it->second).first;
    173167
    174       formula = (it->second).second;
    175       pdgOut = (it->second).first;
     168      p = formula->Eval(pt, eta, phi, e);
    176169
    177       Pi = formula->Eval(pt, eta);
    178 
    179       // check that sum of probabilities does not exceed 1.
    180       P = (P - Pi)/P;
    181 
    182       if( P < 0.0 ) continue;
    183       else
     170      if(total <= r && r < total + p)
    184171      {
    185 
    186        //randomly assign a PID to particle according to map
    187        Double_t rndm = gRandom->Uniform();
    188 
    189        if(rndm > P)
    190        {
    191          //change PID of particle
    192          if(pdgOut != 0) candidate->PID = charge*pdgOut;
    193          fOutputArray->Add(candidate);
    194          break;
    195        }
     172        // change PID of particle
     173        if(pdgCodeOut != 0) candidate->PID = charge*pdgCodeOut;
     174        fOutputArray->Add(candidate);
     175        break;
    196176      }
    197177
     178      total += p;
    198179    }
    199 
    200    }
    201 
     180  }
    202181}
    203182
Note: See TracChangeset for help on using the changeset viewer.