Fork me on GitHub

Changeset a5bbe8a in git for modules


Ignore:
Timestamp:
Dec 21, 2014, 7:48:30 PM (10 years ago)
Author:
Pavel Demin <pavel.demin@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
d4b9697
Parents:
e5767b57
Message:

fix probability calculation in IdentificationMap

Location:
modules
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • modules/IdentificationMap.cc

    re5767b57 ra5bbe8a  
    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
     
    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    pdgCodeIn = candidate->PID;
    155146    charge = candidate->Charge;
    156147
    157    // cout<<"------------ New Candidate ------------"<<endl;
    158    // cout<<candidate->PID<<"   "<<pt<<","<<eta<<","<<phi<<endl;
     148    // first check that PID of this particle is specified in the map
     149    // otherwise, look for PID = 0
    159150
    160     P = 1.0;
     151    itEfficiencyMap = fEfficiencyMap.find(pdgCodeIn);
    161152
    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);
     153    range = fEfficiencyMap.equal_range(pdgCodeIn);
     154    if(range.first == range.second) range = fEfficiencyMap.equal_range(-pdgCodeIn);
    168155    if(range.first == range.second) range = fEfficiencyMap.equal_range(0);
    169156
    170     //loop over submap for this pid
    171     for (TMisIDMap::iterator it=range.first; it!=range.second; ++it)
     157    r = gRandom->Uniform();
     158    total = 0.0;
     159
     160    // loop over sub-map for this PID
     161    for(TMisIDMap::iterator it = range.first; it != range.second; ++it)
    172162    {
     163      formula = (it->second).second;
     164      pdgCodeOut = (it->second).first;
    173165
    174       formula = (it->second).second;
    175       pdgOut = (it->second).first;
     166      p = formula->Eval(pt, eta);
    176167
    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
     168      if(total <= r && r < total + p)
    184169      {
    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        }
     170        // change PID of particle
     171        if(pdgCodeOut != 0) candidate->PID = charge*pdgCodeOut;
     172        fOutputArray->Add(candidate);
     173        break;
    196174      }
    197175
     176      total += p;
    198177    }
    199 
    200    }
    201 
     178  }
    202179}
    203180
  • modules/IdentificationMap.h

    re5767b57 ra5bbe8a  
    4949private:
    5050
    51   typedef std::multimap< Int_t, std::pair<Int_t , DelphesFormula * > > TMisIDMap; //!
     51  typedef std::multimap< Int_t, std::pair< Int_t, DelphesFormula * > > TMisIDMap; //!
    5252
    53   TMisIDMap fEfficiencyMap;
     53  TMisIDMap fEfficiencyMap; //!
    5454
    5555  TIterator *fItInputArray; //!
Note: See TracChangeset for help on using the changeset viewer.