Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New mip treatment 2 #113

Closed
wants to merge 9 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions FastSimulation/Calorimetry/interface/HCALResponse.h
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,9 @@ class HCALResponse
// mip = 2 means "mean" response regardless actual mip
double responseHCAL(int _mip, double energy, double eta, int partype);

//Get the energy and eta dependent mip fraction
double getMIPfraction(double energy, double eta);

// legacy methods using simple formulae
double getHCALEnergyResponse(double e, int hit);

Expand Down Expand Up @@ -101,6 +104,7 @@ class HCALResponse
// muon histos
// indices: responseMU[energy][eta][bin]
vec3 responseMU;
vec3 mipfraction;

// Famos random engine
const RandomEngine* random;
Expand Down
53 changes: 51 additions & 2 deletions FastSimulation/Calorimetry/python/HcalResponse_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -896,7 +896,56 @@
108.11042, 116.28395, 120.52682, 116.59945, 134.65273, 125.06872, 127.37372, 131.68186, 135.85193, 140.72941, 139.99588, 146.59417, 129.46715, 141.59213, 151.45253, 3.90576, 65.59708, 45.40182, 16.60240, 1.01002, #1000
113.89983, 122.30338, 120.67639, 118.65816, 129.41775, 135.38712, 126.18339, 141.89882, 144.88250, 126.33036, 11.82460, 14.81691, 1.01000, 2.47610, 3.38382, 3.54591, 2.34644, 8.70430, 166.11478, 1.01000 ] #3000
),


f_mip_HB = cms.vdouble( *[
0.97421, 0.97097, 0.97074, 0.97488, 0.97791, 0.97398, 0.97588, 0.98007, 0.98398, 0.98396, 0.98807, 0.99079, 0.99170, 0.99238, 0.99494, 0.99607, #1
0.53277, 0.51659, 0.52109, 0.52950, 0.54883, 0.54117, 0.54443, 0.57778, 0.57320, 0.58880, 0.62213, 0.67456, 0.68939, 0.70992, 0.77665, 0.81718, #2
0.43487, 0.41283, 0.40378, 0.40406, 0.42833, 0.41668, 0.41240, 0.43882, 0.42568, 0.42633, 0.43602, 0.47671, 0.48651, 0.50125, 0.56429, 0.64965, #3
0.39616, 0.37508, 0.36708, 0.36680, 0.37725, 0.36048, 0.36680, 0.37988, 0.35196, 0.34160, 0.33828, 0.35790, 0.35327, 0.35615, 0.43074, 0.49421, #5
0.37106, 0.35794, 0.33717, 0.34122, 0.35064, 0.33553, 0.32955, 0.32902, 0.30803, 0.29422, 0.27866, 0.27681, 0.24911, 0.25587, 0.31394, 0.34733, #9
0.36939, 0.35728, 0.34662, 0.33784, 0.34244, 0.33774, 0.33343, 0.34176, 0.30816, 0.28949, 0.27744, 0.26656, 0.25620, 0.24061, 0.28918, 0.33004, #11
0.36506, 0.35876, 0.34990, 0.34575, 0.35501, 0.33175, 0.32858, 0.33924, 0.30566, 0.28754, 0.27834, 0.27160, 0.23347, 0.23038, 0.28632, 0.30970, #15
0.37756, 0.37268, 0.35152, 0.33697, 0.35208, 0.33986, 0.33675, 0.33959, 0.30494, 0.28411, 0.28194, 0.26983, 0.23786, 0.23273, 0.28151, 0.31958, #20
0.37558, 0.36670, 0.35471, 0.34682, 0.35926, 0.33765, 0.32953, 0.33973, 0.30641, 0.28765, 0.26840, 0.26171, 0.22832, 0.22373, 0.27078, 0.30927, #30
0.36568, 0.34685, 0.35615, 0.32804, 0.34532, 0.33953, 0.32270, 0.32192, 0.29795, 0.27939, 0.26075, 0.25055, 0.21651, 0.21725, 0.27376, 0.30446, #50
0.34547, 0.33897, 0.32693, 0.31872, 0.33163, 0.30426, 0.30069, 0.30882, 0.27316, 0.25439, 0.23999, 0.23656, 0.18286, 0.18518, 0.26812, 0.30828, #100
0.33340, 0.31901, 0.30205, 0.30504, 0.31130, 0.28607, 0.28387, 0.29463, 0.25310, 0.23397, 0.21072, 0.20689, 0.15714, 0.16812, 0.25404, 0.28392, #150
0.30253, 0.30323, 0.27515, 0.27576, 0.28682, 0.26440, 0.24795, 0.26074, 0.21725, 0.20939, 0.19017, 0.17910, 0.12227, 0.14179, 0.23668, 0.27608, #225
0.29090, 0.26856, 0.25165, 0.25174, 0.27146, 0.23324, 0.22383, 0.24070, 0.20317, 0.18744, 0.17008, 0.15445, 0.10284, 0.13657, 0.22832, 0.26915, #300
0.18384, 0.17105, 0.16635, 0.15558, 0.16177, 0.14182, 0.13942, 0.14725, 0.11509, 0.10264, 0.08775, 0.07937, 0.04351, 0.07630, 0.18804, 0.20948, #1000
0.09650, 0.08967, 0.08041, 0.08125, 0.08452, 0.06989, 0.06312, 0.07437, 0.05247, 0.04502, 0.04377, 0.02963, 0.01353, 0.03486, 0.13261, 0.15558 ] #3000
),
f_mip_HE = cms.vdouble( *[
0.99429, 0.99470, 0.99554, 0.99474, 0.99459, 0.99390, 0.99364, 0.99441, 0.99423, 0.99397, 0.99372, 0.99215, 0.99618, 0.99857, #1
0.71084, 0.71267, 0.71142, 0.71727, 0.71976, 0.70605, 0.71374, 0.69655, 0.72856, 0.73079, 0.72505, 0.72142, 0.80199, 0.94369, #2
0.55383, 0.54833, 0.54737, 0.55910, 0.56216, 0.56889, 0.55886, 0.55588, 0.58548, 0.58822, 0.58324, 0.58123, 0.62950, 0.85925, #3
0.41426, 0.41792, 0.41450, 0.42887, 0.44435, 0.44461, 0.45633, 0.45664, 0.47544, 0.48763, 0.49746, 0.48061, 0.49873, 0.74229, #5
0.28631, 0.29634, 0.30339, 0.31450, 0.32598, 0.33113, 0.33803, 0.34820, 0.36268, 0.37627, 0.38924, 0.40029, 0.38658, 0.61738, #9
0.26654, 0.26989, 0.27516, 0.28648, 0.29743, 0.31362, 0.31887, 0.33636, 0.33665, 0.34946, 0.36885, 0.38871, 0.36999, 0.60707, #11
0.24445, 0.26080, 0.25045, 0.27080, 0.26225, 0.28372, 0.28250, 0.30529, 0.29972, 0.30712, 0.33278, 0.36247, 0.33486, 0.56842, #15
0.25251, 0.25010, 0.24310, 0.25404, 0.24810, 0.27070, 0.27268, 0.29224, 0.28615, 0.28841, 0.30734, 0.32993, 0.32243, 0.53050, #20
0.24741, 0.25035, 0.24194, 0.25060, 0.24644, 0.26411, 0.26911, 0.27443, 0.26398, 0.27609, 0.28997, 0.33629, 0.30714, 0.51338, #30
0.24610, 0.25224, 0.23687, 0.25137, 0.25578, 0.25856, 0.25814, 0.27202, 0.26798, 0.26745, 0.28856, 0.31921, 0.29920, 0.49866, #50
0.23329, 0.23228, 0.22574, 0.24405, 0.25053, 0.25262, 0.26126, 0.26116, 0.25810, 0.25921, 0.28177, 0.31705, 0.29397, 0.48651, #100
0.22955, 0.22464, 0.21474, 0.21835, 0.23654, 0.23484, 0.23957, 0.26140, 0.24523, 0.25540, 0.26628, 0.29354, 0.28267, 0.48171, #150
0.21477, 0.20715, 0.19498, 0.20039, 0.20960, 0.22549, 0.22204, 0.24198, 0.22433, 0.23217, 0.25065, 0.28216, 0.26425, 0.46886, #225
0.19793, 0.19363, 0.18807, 0.19255, 0.19163, 0.20793, 0.21452, 0.23016, 0.21682, 0.22363, 0.24070, 0.27456, 0.25035, 0.46689, #300
0.14020, 0.12920, 0.12511, 0.13196, 0.13438, 0.13828, 0.15417, 0.16181, 0.15755, 0.15575, 0.17157, 0.21012, 0.18866, 0.42957, #1000
0.07904, 0.06727, 0.06703, 0.06885, 0.07529, 0.07981, 0.08565, 0.09368, 0.09294, 0.09604, 0.11844, 0.14221, 0.12614, 0.38910 ] #3000
),
f_mip_HF = cms.vdouble( *[
0.88658, 0.94872, 0.98165, 0.99411, 0.99960, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #15
0.86531, 0.93844, 0.97610, 0.99275, 0.99980, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #20
0.84301, 0.91987, 0.97340, 0.99266, 0.99990, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #30
0.82148, 0.91559, 0.96855, 0.98924, 0.99950, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #50
0.81098, 0.90218, 0.96255, 0.98793, 0.99890, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #100
0.80433, 0.89087, 0.96120, 0.98545, 0.99950, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #150
0.79887, 0.88446, 0.95396, 0.98387, 0.99889, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #225
0.78883, 0.87948, 0.95116, 0.97972, 0.99850, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #300
0.76272, 0.85972, 0.93617, 0.97497, 0.99858, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, #1000
0.72955, 0.83651, 0.92511, 0.96447, 0.99733, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000 ] #3000
),

#muon values
muStep = cms.double(0.25),
maxMUe = cms.int32(4),
Expand Down Expand Up @@ -986,4 +1035,4 @@
eResponsePlateauHB = cms.double(0.95),
ElectronForwardResolution_Noise = cms.double(0.0)
)
)
)
8 changes: 7 additions & 1 deletion FastSimulation/Calorimetry/src/CalorimetryManager.cc
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
//updated by Reza Goldouzian
//Framework headers
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
Expand Down Expand Up @@ -608,7 +609,11 @@ void CalorimetryManager::HDShowerSimulation(const FSimTrack& myTrack){//,
double eint = moment.e();
double eGen = myTrack.hcalEntrance().e();

double emeas = 0.;
double emeas = 0.;
double pmip= myHDResponse_->getMIPfraction(eGen, pathEta);
// std::cout << " CalorimetryManager onHcal " << pmip << std::endl;


//double emeas = -0.000001;

//===========================================================================
Expand Down Expand Up @@ -704,6 +709,7 @@ void CalorimetryManager::HDShowerSimulation(const FSimTrack& myTrack){//,
&myHcalHitMaker,
onECAL,
eGen,
pmip,
dbe);
status = theShower.compute();
mip = theShower.getmip();
Expand Down
52 changes: 49 additions & 3 deletions FastSimulation/Calorimetry/src/HCALResponse.cc
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
//updated by Reza Goldouzian
//FastSimulation headers
#include "FastSimulation/Calorimetry/interface/HCALResponse.h"
#include "FastSimulation/Utilities/interface/RandomEngine.h"
Expand Down Expand Up @@ -75,7 +76,7 @@ HCALResponse::HCALResponse(const edm::ParameterSet& pset, const RandomEngine* en
parNames = pset.getParameter<std::vector<std::string> >("parNames");
std::string detNames[] = {"_HB","_HE","_HF"};
std::string mipNames[] = {"_mip","_nomip",""};

std::string fraction="f";
//setup parameters (5D vector)
parameters = vec5(nPar,vec4(3,vec3(3)));
for(int p = 0; p < nPar; p++){ //loop over parameters
Expand All @@ -101,6 +102,24 @@ HCALResponse::HCALResponse(const edm::ParameterSet& pset, const RandomEngine* en
}
}

//MIP fraction fill in 3d vector
////--------------------------------------------------------------------
mipfraction = vec3(3);
for(int d = 0; d < 3; d++){ //loop over dets: HB, HE, HF
//get from python
std::string mipname = fraction + mipNames[0] + detNames[d] ;
vec1 tmp1 = pset.getParameter<vec1>(mipname);
mipfraction[d].resize(maxHDe[d]);
for(int i = 0; i < maxHDe[d]; i++){ //loop over energy for det d
//resize vector for eta range of det d
mipfraction[d][i].resize(maxHDetas[d]);
for(int j = 0; j < maxHDetas[d]; j++){ //loop over eta for det d
//fill in parameters vector from python
mipfraction[d][i][j]= tmp1[i*maxHDetas[d] + j];
}
}
}

// MUON probability histos for bin size = 0.25 GeV (0-10 GeV, 40 bins)
//--------------------------------------------------------------------
muStep = pset.getParameter<double>("muStep");
Expand Down Expand Up @@ -171,7 +190,33 @@ HCALResponse::HCALResponse(const edm::ParameterSet& pset, const RandomEngine* en

}


double HCALResponse::getMIPfraction(double energy, double eta){
int ieta = abs((int)(eta / etaStep)) ;
int ie = -1;
//check eta and det
int det = getDet(ieta);
int deta = ieta - HDeta[det];
if(deta >= maxHDetas[det]) deta = maxHDetas[det] - 1;
else if(deta < 0 ) deta = 0;
//find energy range
for (int i = 0; i < maxHDe[det]; i++) {
if(energy < eGridHD[det][i]) {
if(i == 0) return mipfraction [det][0][deta]; // less than minimal - the first value is used instead of extrapolating
else ie = i-1;
break;
}
}
if(ie == -1) return mipfraction [det][maxHDe[det]][deta]; // more than maximal - the last value is used instead of extrapolating
double y1, y2;
double x1 = eGridHD[det][ie];
double x2 = eGridHD[det][ie+1];
y1=mipfraction[det][ie][deta];
y2=mipfraction[det][ie+1][deta];
double mean = 0;
mean=(y1*(x2-energy) + y2*(energy-x1))/(x2-x1);
return mean;
}

double HCALResponse::responseHCAL(int _mip, double energy, double eta, int partype){
int ieta = abs((int)(eta / etaStep)) ;
int ie = -1;
Expand Down Expand Up @@ -316,6 +361,7 @@ double HCALResponse::interMU(double e, int ie, int ieta) {

double HCALResponse::interHD(int mip, double e, int ie, int ieta, int det) {
double y1, y2;
if(det==2) mip=2; //ignore mip status for HF

double x1 = eGridHD[det][ie];
double x2 = eGridHD[det][ie+1];
Expand Down Expand Up @@ -459,4 +505,4 @@ double HCALResponse::cballShootNoNegative(double mu, double sigma, double aL, do
//else give up on re-trying, otherwise too much time can be lost before emeas comes out positive

return out;
}
}
1 change: 1 addition & 0 deletions FastSimulation/ShowerDevelopment/interface/HDShower.h
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ class HDShower
HcalHitMaker* myHcalHitMaker,
int onECAL,
double epart,
double pmip,
DQMStore * const dbeIn);

int getmip() {return mip;}
Expand Down
8 changes: 8 additions & 0 deletions FastSimulation/ShowerDevelopment/src/HDShower.cc
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
//updated by Reza Goldouzian
//FastSimulation Headers
#include "FastSimulation/ShowerDevelopment/interface/HDShower.h"
//#include "FastSimulation/Utilities/interface/Histos.h"
Expand Down Expand Up @@ -36,12 +37,14 @@ HDShower::HDShower(const RandomEngine* engine,
HcalHitMaker* myHcalHitMaker,
int onECAL,
double epart,
double pmip,
DQMStore * const dbeIn)
: theParam(myParam),
theGrid(myGrid),
theHcalHitMaker(myHcalHitMaker),
onEcal(onECAL),
e(epart),
// pmip(pmip),
random(engine),
dbe(dbeIn)
{
Expand Down Expand Up @@ -240,6 +243,11 @@ HDShower::HDShower(const RandomEngine* engine,
// if no HCAL material behind - force to deposit in ECAL
double maxDepth = depthToHCAL + depthHCAL - 1.1 * depthStep;
double depthStart = std::log(1./random->flatShoot()); // starting point lambda unts
//std::cout<<"generated depth "<<depthStart<<std::endl;
if (pmip==1) {depthStart=depthToHCAL;}
else {depthStart=depthStart*0.9*depthECAL/std::log(1./pmip);}
// std::cout<<"modified depth "<<depthStart<<std::endl;


if(e < emin) {
if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : e <emin -> depthStart = 0" << std::endl;
Expand Down
2 changes: 2 additions & 0 deletions Tutorial/Test/bin/BuildFile.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
<bin file="*.cc" name="git-tutorial-test">
</bin>