Several changes; mainly for BH and adding HepMC pointer to Particle

doc/wup.tex

@@ -15,6 +15,53 @@
 \bigskip
 \centerline{\large May 2016}
 
+\section{Optimisation of layer weights}
+Each event $e$ gives some values $d_{e,i}$ of the deposited energy in 
+layer $i$; these
+can be in any units, e.g. MIPs.
+Assume these are to be multiplied by some constant coefficients $a_i$ (which
+are approximately the integrated dE/dx values if the $d_{e,i}$ are in MIPs)
+to give the estimate of the incoming EM photon or electron energy. 
+Hence, the energy
+estimation for event $e$ is
+\begin{equation} 
+E_e = \sum_i a_i d_{e,i}
+\end{equation} 
+To find the optimal coefficients, then we need to know the truth energy per
+event $T_e$. For a given set of coefficients, the RMS$^2$ of the energy 
+around the truth value is given by
+\begin{equation} 
+\mathrm{RMS}^2 = \frac{1}{N} \sum_e (E_e - T_e)^2
+= \frac{1}{N} \sum_e \left(\sum_i a_i d_{e,i} - T_e\right)^2
+\end{equation} 
+This can be thought of as similar to a chi-squared; we want to minimise this
+expression. If all the $a_i$ are considered as independent parameters (so 28 for
+the EE only), then explicitly
+\begin{equation} 
+\frac{\partial \mathrm{RMS}^2}{\partial a_j}
+= \frac{1}{N} \sum_e 2d_{e,j} \left(\sum_i a_i d_{e,i} - T_e\right)
+= \frac{2}{N} \sum_i a_i \left(\sum_e d_{e,j} d_{e,i} \right)
+- \frac{2}{N} \sum_e d_{e,j} T_e 
+\end{equation} 
+Hence, for the minimum, we require
+\begin{equation} 
+\sum_i \frac{\sum_e d_{e,j} d_{e,i}}{N} a_i
+= \frac{\sum_e d_{e,j} T_e}{N} 
+\end{equation} 
+Writing in matrix notation with $M$ and $v$ defined as
+\begin{equation} 
+M_{ji} = \frac{\sum_e d_{e,j} d_{e,i}}{N},\qquad
+v_j = \frac{\sum_e d_{e,j} T_e}{N}
+\end{equation} 
+then the requirement is
+\begin{equation} 
+M a = v\qquad\mathrm{so}\qquad a = M^{-1}v
+\end{equation}
+Inverting the large matrix $M$ is required to give the solution for the
+optimal $a_i$.
+Note, $M$ is similar (but not identical) to the error matrix of the $d_i$.
+
+
 \section{Units}
 Keeping quantities to 16-bit integers.
 

interface/AnalysisPhoton.hh

@@ -63,7 +63,6 @@ public:
 
       fNum=1.0;
     }
-    /*    
     TMatrixDSym fitMatrix() const {
       TMatrixDSym fm(40);
       for(unsigned i(0);i<40;i++) {
@@ -74,6 +73,7 @@ public:
       return fm;
     }
 
+    /*    
     TMatrixDSym errorMatrix() const {
       TMatrixDSym em(40);
       for(unsigned i(0);i<40;i++) {

interface/Event.hh

@@ -177,15 +177,15 @@ public:
 
   void fillGenVertex(HepMC::GenVertex *vtx, bool signal=false) {
     if(vtx==0) return;
-    std::cout << "New vertex with barcode " << vtx->barcode() << std::endl;
+    //std::cout << "New vertex with barcode " << vtx->barcode() << std::endl;
 
     // Loop over all particles from the vertex
     for (HepMC::GenVertex::particles_out_const_iterator p(vtx->particles_out_const_begin());p!= vtx->particles_out_const_end();++p) {
 
       const HepMC::GenParticle &ptc(*(*p));
       int pdgId(ptc.pdg_id());
-      std::cout << "New particle with barcode " << ptc.barcode()
-		<< (abs(pdgId)==11?" ELECTRON!!!":"") << std::endl;
+      //std::cout << "New particle with barcode " << ptc.barcode()
+      //<< (abs(pdgId)==11?" ELECTRON!!!":"") << std::endl;
 	  
       HepMC::FourVector ptc4(ptc.momentum());
       unsigned e(ptc4.z()<0.0?0:1);
@@ -228,15 +228,15 @@ public:
 
   void fillGenVertexList(const HepMC::GenVertex *vtx, std::vector<Particle> &v) {
     if(vtx==0) return;
-    std::cout << "New vertex with barcode " << vtx->barcode() << std::endl;
+    //std::cout << "New vertex with barcode " << vtx->barcode() << std::endl;
     
     // Loop over all particles from the vertex
     for (HepMC::GenVertex::particles_out_const_iterator p(vtx->particles_out_const_begin());p!= vtx->particles_out_const_end();++p) {
       
       const HepMC::GenParticle &ptc(*(*p));
-      int pdgId(ptc.pdg_id());
-      std::cout << "New particle with barcode " << ptc.barcode()
-		<< (abs(pdgId)==11?" ELECTRON!!!":"") << std::endl;
+      //int pdgId(ptc.pdg_id());
+      //std::cout << "New particle with barcode " << ptc.barcode()
+      //<< (abs(pdgId)==11?" ELECTRON!!!":"") << std::endl;
       
       HepMC::FourVector ptc4(ptc.momentum());
       //unsigned e(ptc4.z()<0.0?0:1);
@@ -301,6 +301,7 @@ public:
 	      if(abs(pdgId)==15) {
 		std::cout << std::endl << "New tau " << std::endl;
 		vTau[e].push_back(Particle(ptc));
+		vTau[e][vTau[e].size()-1].print();
 		vTauDaughters[e].push_back(std::vector<Particle>());
 		fillGenVertexList(&vtx,vTauDaughters[e][vTauDaughters[e].size()-1]);	
 		std::cout << std::endl << "New tau daughters = " 

interface/Geometry.hh

@@ -233,7 +233,7 @@ public:
     }
 
     for(unsigned i(0);i<kNumberOfBHLayers;i++) {
-      fLayerZBh[i]=420.0+10.0*i; // NOT CORRECT!!!
+      fLayerZBh[i]=421.0+9.0*i; // NOT EXACTLY CORRECT!!!
     }
 
     std::ifstream fin;

interface/Particle.hh

@@ -68,10 +68,12 @@ public:
     return q;
   }
 
-  Particle() {
+  Particle() : pGenParticle(0) {
   }
   
   Particle(const HepMC::GenParticle &p) {
+    pGenParticle=&p;
+
     fPdgId=p.pdg_id();
 
     HepMC::GenVertex *v(p.production_vertex());
@@ -86,6 +88,10 @@ public:
   ~Particle() {
   }
 
+  const HepMC::GenParticle* genParticle() const {
+    return pGenParticle;
+  }
+
   Point position() const {
     return fPosition;
   }
@@ -108,6 +114,8 @@ public:
 
   Particle project(double z) const {
     Particle p;
+    p.pGenParticle=pGenParticle;
+
     p.fPdgId=fPdgId;
 
     double dz(z-fPosition.z());
@@ -142,12 +150,15 @@ public:
     std::cout << "Particle PDG Id = " << fPdgId << std::endl;
     std::cout << " Position ";fPosition.print();
     std::cout << " Momentum ";fMomentum.print();
+    std::cout << " ";pGenParticle->print();
   }
 
 
 private:
   static const double fBField;
 
+  const HepMC::GenParticle *pGenParticle;
+  
   int fPdgId;
   Point fPosition;
   Point fMomentum;

interface/SimDetBh.hh

@@ -21,10 +21,9 @@ public:
     assert(h<1024);
     assert(f<1024);
 
-    fDetId=(5<<29)+(t<<26)+(l<<21)+(e<<20)+(h<<10)+f;
+    fDetId=(5<<29)+(t<<26)+((l+7)<<21)+(e<<20)+((h+1)<<10)+(f+1);
   }
 
-
   ~SimDetBh() {
   }
   
@@ -53,7 +52,7 @@ public:
   }
 
   uint32_t eta() const {
-    return ((fDetId&0x0007fc00)>>10)-1;
+    return ((fDetId&0x000ffc00)>>10)-1;
   }
 
   uint32_t endcap() const {
@@ -93,7 +92,7 @@ public:
   }
 
   void setSignal() {
-    fDetId|=0x00080000;
+    fDetId|=0x10000000;
   }
 
   void flipEndcap() {

src/DetIdTest.cpp

@@ -86,11 +86,11 @@ int main(int argc, const char **argv) {
     }
 
     std::cout << " BH ";
-    for(unsigned l(0);l<32;l++) {
+    for(unsigned l(0);l<25;l++) {
       std::cout << "." << std::flush;
       for(unsigned t(0);t<4;t++) {
-	for(unsigned h(0);h<512;h++) {
-	  for(unsigned f(0);f<1024;f++) {
+	for(unsigned h(0);h<1023;h++) {
+	  for(unsigned f(0);f<1023;f++) {
 	    SimDetBh s(e,l,t,h,f);
 	    assert(s.endcap()==e);
 	    assert(s.layer()==l);

src/readSimHitsTest.cpp

@@ -222,7 +222,8 @@ int main(int argc, const char **argv) {
 
   TH1F *hBhSimNumber,*hBhSimDetector,*hBhSimSubdetector,
     *hBhSimEndcap,*hBhSimLayer,*hBhSimEta;
-  TH2F *hBhSimPhiVsEta,*hBhSimDepthVsLayer;
+  TH2F *hBhSimPhiVsEta,*hBhSimLayerVsEta,
+    *hBhSimDepthVsLayer,*hBhSimSelDepthVsLayer;
 
   hBhSimNumber=new TH1F("BhSimNumber","Sim ",100,0.0,5000.0);
   hBhSimDetector=new TH1F("BhSimDetector","Sim ",10,0.0,10.0);
@@ -233,8 +234,12 @@ int main(int argc, const char **argv) {
   
   hBhSimPhiVsEta=new TH2F("BhSimPhiVsEta","Sim ",
 			  100,0.0,100,360,0.0,360.0);
+  hBhSimLayerVsEta=new TH2F("BhSimLayerVsEta","Sim ",
+			  100,0.0,100,12,0.0,12.0);
   hBhSimDepthVsLayer=new TH2F("BhSimDepthVsLayer","Sim ",
-			  20,0.0,20,4,0.0,4.0);
+			  12,0.0,12,4,0.0,4.0);
+  hBhSimSelDepthVsLayer=new TH2F("BhSimSelDepthVsLayer","Sim ",
+				 12,0.0,12,4,0.0,4.0);
 
   ///////////////////////////////////////////////////////////////////
 
@@ -350,7 +355,9 @@ int main(int argc, const char **argv) {
 
       hBhSimEta->Fill(vEvtSimHbh[evt][i].eta());
       hBhSimPhiVsEta->Fill(vEvtSimHbh[evt][i].eta(),vEvtSimHbh[evt][i].phi());
+      hBhSimLayerVsEta->Fill(vEvtSimHbh[evt][i].eta(),vEvtSimHbh[evt][i].layer());
       hBhSimDepthVsLayer->Fill(vEvtSimHbh[evt][i].layer(),vEvtSimHbh[evt][i].depth());
+      if(vEvtSimHbh[evt][i].eta()>=16) hBhSimSelDepthVsLayer->Fill(vEvtSimHbh[evt][i].layer(),vEvtSimHbh[evt][i].depth());
     }
   }