Merge pull request #138 from giamman/new-mip-treatment-3

New treatment for MIPs in the calorimeters
cms-sw · Jul 18, 2013 · e8d9e5e · e8d9e5e
2 parents 1a44ed3 + 8fabd7d
commit e8d9e5e
Show file tree

Hide file tree

Showing 6 changed files with 120 additions and 6 deletions.
diff --git 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
@@ -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)
     )
-)
+)
diff --git 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
@@ -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<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
@@ -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");
@@ -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;
-}
+}
diff --git 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
@@ -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<<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;