diff --git a/FastSimulation/Calorimetry/interface/HCALResponse.h b/FastSimulation/Calorimetry/interface/HCALResponse.h index 6af7188344230..179580e337062 100644 --- a/FastSimulation/Calorimetry/interface/HCALResponse.h +++ b/FastSimulation/Calorimetry/interface/HCALResponse.h @@ -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); @@ -101,6 +104,7 @@ class HCALResponse // muon histos // indices: responseMU[energy][eta][bin] vec3 responseMU; + vec3 mipfraction; // Famos random engine const RandomEngine* random; diff --git a/FastSimulation/Calorimetry/python/HcalResponse_cfi.py b/FastSimulation/Calorimetry/python/HcalResponse_cfi.py index 3b1db2ea42857..20560c8d0ab8a 100644 --- a/FastSimulation/Calorimetry/python/HcalResponse_cfi.py +++ b/FastSimulation/Calorimetry/python/HcalResponse_cfi.py @@ -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), @@ -986,4 +1035,4 @@ eResponsePlateauHB = cms.double(0.95), ElectronForwardResolution_Noise = cms.double(0.0) ) -) \ No newline at end of file +) diff --git a/FastSimulation/Calorimetry/src/CalorimetryManager.cc b/FastSimulation/Calorimetry/src/CalorimetryManager.cc index af881f7bf6f16..ff7b1593bab4a 100644 --- a/FastSimulation/Calorimetry/src/CalorimetryManager.cc +++ b/FastSimulation/Calorimetry/src/CalorimetryManager.cc @@ -1,3 +1,4 @@ +//updated by Reza Goldouzian //Framework headers #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -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; //=========================================================================== @@ -704,6 +709,7 @@ void CalorimetryManager::HDShowerSimulation(const FSimTrack& myTrack){//, &myHcalHitMaker, onECAL, eGen, + pmip, dbe); status = theShower.compute(); mip = theShower.getmip(); diff --git a/FastSimulation/Calorimetry/src/HCALResponse.cc b/FastSimulation/Calorimetry/src/HCALResponse.cc index 1575e77491c3d..01f616fea1785 100644 --- a/FastSimulation/Calorimetry/src/HCALResponse.cc +++ b/FastSimulation/Calorimetry/src/HCALResponse.cc @@ -1,3 +1,4 @@ +//updated by Reza Goldouzian //FastSimulation headers #include "FastSimulation/Calorimetry/interface/HCALResponse.h" #include "FastSimulation/Utilities/interface/RandomEngine.h" @@ -75,7 +76,7 @@ HCALResponse::HCALResponse(const edm::ParameterSet& pset, const RandomEngine* en parNames = pset.getParameter >("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 @@ -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(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("muStep"); @@ -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; @@ -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]; @@ -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; -} \ No newline at end of file +} diff --git a/FastSimulation/ShowerDevelopment/interface/HDShower.h b/FastSimulation/ShowerDevelopment/interface/HDShower.h index 91cc15448e8bc..8217970d170ea 100644 --- a/FastSimulation/ShowerDevelopment/interface/HDShower.h +++ b/FastSimulation/ShowerDevelopment/interface/HDShower.h @@ -38,6 +38,7 @@ class HDShower HcalHitMaker* myHcalHitMaker, int onECAL, double epart, + double pmip, DQMStore * const dbeIn); int getmip() {return mip;} diff --git a/FastSimulation/ShowerDevelopment/src/HDShower.cc b/FastSimulation/ShowerDevelopment/src/HDShower.cc index 162c5c53552eb..6cd77e1b31c05 100644 --- a/FastSimulation/ShowerDevelopment/src/HDShower.cc +++ b/FastSimulation/ShowerDevelopment/src/HDShower.cc @@ -1,3 +1,4 @@ +//updated by Reza Goldouzian //FastSimulation Headers #include "FastSimulation/ShowerDevelopment/interface/HDShower.h" //#include "FastSimulation/Utilities/interface/Histos.h" @@ -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) { @@ -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 = 0" << std::endl;