Belle II Software  release-05-02-19
OverlapResidualsModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Filippo Dattola *
7  * *
8  * Module to study VXD hits from overlapping sensors of a same VXD layer. *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 
12 #include <tracking/modules/VXDOverlaps/OverlapResidualsModule.h>
13 #include <framework/datastore/StoreArray.h>
14 #include <framework/gearbox/Unit.h>
15 #include <framework/core/Environment.h>
16 #include <svd/dataobjects/SVDCluster.h>
17 #include <svd/dataobjects/SVDTrueHit.h>
18 #include <vxd/dataobjects/VxdID.h>
19 #include <vxd/geometry/GeoCache.h>
20 #include <vxd/geometry/SensorInfoBase.h>
21 #include <tracking/dataobjects/RecoTrack.h>
22 #include <tracking/dataobjects/RecoHitInformation.h>
23 #include <genfit/TrackPoint.h>
24 #include <TVector3.h>
25 #include <TDirectory.h>
26 #include <math.h>
27 #include <iostream>
28 
29 using namespace Belle2;
30 using namespace std;
31 
32 
33 //-----------------------------------------------------------------
34 // Register the Module
35 //-----------------------------------------------------------------
36 REG_MODULE(OverlapResiduals)
37 
38 //-----------------------------------------------------------------
39 // Implementation
40 //-----------------------------------------------------------------
41 
42 OverlapResidualsModule::OverlapResidualsModule() : HistoModule()
43 {
44  //Set module properties
45  setDescription("The module studies consecutive hits in overlapping sensors of a same VXD layer, and the differences of their residuals, to monitor the detector alignment.");
46  //Parameter to take only specific RecoTracks as input
47  addParam("recoTracksStoreArrayName", m_recoTracksStoreArrayName, "StoreArray name of the input and output RecoTracks.",
48  m_recoTracksStoreArrayName);
49  //Parameter to produce TTrees storing global information on VXD overlaps
50  addParam("ExpertLevel", m_ExpertLevel,
51  "If set, enables the production of TTrees containing low level information on PXD and SVD overlaps", false);
52 }
53 
54 void OverlapResidualsModule::initialize()
55 {
56  StoreArray<RecoTrack> recoTracks(m_recoTracksStoreArrayName);
57  recoTracks.isOptional();
58  //Register histograms (calls back defineHisto)
59  REG_HISTOGRAM
60 }
61 
62 void OverlapResidualsModule::defineHisto()
63 {
64  //Create directory to store monitoring histograms
65  TDirectory* oldDir = gDirectory;
66  TDirectory* ResDir = oldDir->mkdir("Monitoring_VXDOverlaps");
67  ResDir->cd();
68  //Define histograms of residuals differences
69  h_U_DeltaRes = new TH1F("h_U_DeltaRes", "Histrogram of residual difference #Delta res_{u} for overlapping hits", 100, -1000, 1000);
70  h_U_DeltaRes->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
71  h_U_DeltaRes->GetYaxis()->SetTitle("counts");
72  h_V_DeltaRes = new TH1F("h_V_DeltaRes", "Histrogram of residual difference #Delta res_{v} for overlapping hits", 100, -1000, 1000);
73  h_V_DeltaRes->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
74  h_V_DeltaRes->GetYaxis()->SetTitle("counts");
75  h_U_DeltaRes_PXD = new TH1F("h_U_DeltaRes_PXD", "Histrogram of residual difference #Delta res_{u} for overlapping PXD hits", 100,
76  -1000, 1000);
77  h_U_DeltaRes_PXD->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
78  h_U_DeltaRes_PXD->GetYaxis()->SetTitle("counts");
79  h_U_DeltaRes_PXD_Lyr1 = new TH1F("h_U_DeltaRes_PXD_Lyr1",
80  "Layer1: histrogram of residual difference #Delta res_{u} for overlapping PXD hits", 100, -1000, 1000);
81  h_U_DeltaRes_PXD_Lyr1->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
82  h_U_DeltaRes_PXD_Lyr1->GetYaxis()->SetTitle("counts");
83  h_U_DeltaRes_PXD_Lyr2 = new TH1F("h_U_DeltaRes_PXD_Lyr2",
84  "Layer 2: hstrogram of residual difference #Delta res_{u} for overlapping PXD hits", 100, -1000, 1000);
85  h_U_DeltaRes_PXD_Lyr2->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
86  h_U_DeltaRes_PXD_Lyr2->GetYaxis()->SetTitle("counts");
87  h_V_DeltaRes_PXD = new TH1F("h_V_DeltaRes_PXD", "Histrogram of residual difference #Delta res_{v} for overlapping PXD hits", 100,
88  -1000, 1000);
89  h_V_DeltaRes_PXD->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
90  h_V_DeltaRes_PXD->GetYaxis()->SetTitle("counts");
91  h_V_DeltaRes_PXD_Lyr1 = new TH1F("h_V_DeltaRes_PXD_Lyr1",
92  "Layer 1: histrogram of residual difference #Delta res_{v} for overlapping PXD hits", 100, -1000, 1000);
93  h_V_DeltaRes_PXD_Lyr1->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
94  h_V_DeltaRes_PXD_Lyr1->GetYaxis()->SetTitle("counts");
95  h_V_DeltaRes_PXD_Lyr2 = new TH1F("h_V_DeltaRes_PXD_Lyr2",
96  "Layer 2: histrogram of residual difference #Delta res_{v} for overlapping PXD hits", 100, -1000, 1000);
97  h_V_DeltaRes_PXD_Lyr2->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
98  h_V_DeltaRes_PXD_Lyr2->GetYaxis()->SetTitle("counts");
99  h_U_DeltaRes_SVD = new TH1F("h_U_DeltaRes_SVD", "Histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100,
100  -1000, 1000);
101  h_U_DeltaRes_SVD->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
102  h_U_DeltaRes_SVD->GetYaxis()->SetTitle("counts");
103  h_U_DeltaRes_SVD_Lyr3 = new TH1F("h_U_DeltaRes_SVD_Lyr3",
104  "Layer 3: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
105  h_U_DeltaRes_SVD_Lyr3->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
106  h_U_DeltaRes_SVD_Lyr3->GetYaxis()->SetTitle("counts");
107  h_U_DeltaRes_SVD_Lyr4 = new TH1F("h_U_DeltaRes_SVD_Lyr4",
108  "Layer 4: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
109  h_U_DeltaRes_SVD_Lyr4->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
110  h_U_DeltaRes_SVD_Lyr4->GetYaxis()->SetTitle("counts");
111  h_U_DeltaRes_SVD_Lyr5 = new TH1F("h_U_DeltaRes_SVD_Lyr5",
112  "Layer 5: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
113  h_U_DeltaRes_SVD_Lyr5->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
114  h_U_DeltaRes_SVD_Lyr5->GetYaxis()->SetTitle("counts");
115  h_U_DeltaRes_SVD_Lyr6 = new TH1F("h_U_DeltaRes_SVD_Lyr6",
116  "Layer 6: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
117  h_U_DeltaRes_SVD_Lyr6->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
118  h_U_DeltaRes_SVD_Lyr6->GetYaxis()->SetTitle("counts");
119  h_V_DeltaRes_SVD = new TH1F("h_V_DeltaRes_SVD", "Histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100,
120  -1000, 1000);
121  h_V_DeltaRes_SVD->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
122  h_V_DeltaRes_SVD->GetYaxis()->SetTitle("counts");
123  h_V_DeltaRes_SVD_Lyr3 = new TH1F("h_V_DeltaRes_SVD_Lyr3",
124  "Layer 3: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
125  h_V_DeltaRes_SVD_Lyr3->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
126  h_V_DeltaRes_SVD_Lyr3->GetYaxis()->SetTitle("counts");
127  h_V_DeltaRes_SVD_Lyr4 = new TH1F("h_V_DeltaRes_SVD_Lyr4",
128  "Layer 4: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
129  h_V_DeltaRes_SVD_Lyr4->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
130  h_V_DeltaRes_SVD_Lyr4->GetYaxis()->SetTitle("counts");
131  h_V_DeltaRes_SVD_Lyr5 = new TH1F("h_V_DeltaRes_SVD_Lyr5",
132  "Layer 5: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
133  h_V_DeltaRes_SVD_Lyr5->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
134  h_V_DeltaRes_SVD_Lyr5->GetYaxis()->SetTitle("counts");
135  h_V_DeltaRes_SVD_Lyr6 = new TH1F("h_V_DeltaRes_SVD_Lyr6",
136  "Layer 6: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
137  h_V_DeltaRes_SVD_Lyr6->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
138  h_V_DeltaRes_SVD_Lyr6->GetYaxis()->SetTitle("counts");
139  h_SVDstrips_Mult = new TH1F("h_SVDstrips_Mult", "SVD strips multipicity for SVD clusters in overlapping sensors", 15, 0.5, 15.5);
140  h_SVDstrips_Mult->GetXaxis()->SetTitle("N. of SVD strips contributing to the cluster");
141  h_SVDstrips_Mult->GetYaxis()->SetTitle("counts");
142  //Define 2D histograms: difference of u-residuals vs phi of VXD overlaps for each layer (1 to 6)
143  h_DeltaResUPhi_Lyr1 = new TH2F("h_DeltaResUPhi_Lyr1", "Layer 1: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
144  h_DeltaResUPhi_Lyr1->GetXaxis()->SetTitle("#phi(rad)");
145  h_DeltaResUPhi_Lyr1->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
146  h_DeltaResUPhi_Lyr2 = new TH2F("h_DeltaResUPhi_Lyr2", "Layer 2: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
147  h_DeltaResUPhi_Lyr2->GetXaxis()->SetTitle("#phi(rad)");
148  h_DeltaResUPhi_Lyr2->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
149  h_DeltaResUPhi_Lyr3 = new TH2F("h_DeltaResUPhi_Lyr3", "Layer 3: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
150  h_DeltaResUPhi_Lyr3->GetXaxis()->SetTitle("#phi(rad)");
151  h_DeltaResUPhi_Lyr3->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
152  h_DeltaResUPhi_Lyr4 = new TH2F("h_DeltaResUPhi_Lyr4", "Layer 4: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
153  h_DeltaResUPhi_Lyr4->GetXaxis()->SetTitle("#phi(rad)");
154  h_DeltaResUPhi_Lyr4->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
155  h_DeltaResUPhi_Lyr5 = new TH2F("h_DeltaResUPhi_Lyr5", "Layer 5: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
156  h_DeltaResUPhi_Lyr5->GetXaxis()->SetTitle("#phi(rad)");
157  h_DeltaResUPhi_Lyr5->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
158  h_DeltaResUPhi_Lyr6 = new TH2F("h_DeltaResUPhi_Lyr6", "Layer 6: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
159  h_DeltaResUPhi_Lyr6->GetXaxis()->SetTitle("#phi(rad)");
160  h_DeltaResUPhi_Lyr6->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
161  //Define 2D histograms: difference of u-residuals vs z of VXD overlaps for each layer (1 to 6)
162  h_DeltaResUz_Lyr1 = new TH2F("h_DeltaResUz_Lyr1", "Layer 1: #Delta res_{u} vs z", 100, -4, 8, 100, -200, 200);
163  h_DeltaResUz_Lyr1->GetXaxis()->SetTitle("z (cm)");
164  h_DeltaResUz_Lyr1->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
165  h_DeltaResUz_Lyr2 = new TH2F("h_DeltaResUz_Lyr2", "Layer 2: #Delta res_{u} vs z", 100, -10, 15, 100, -200, 200);
166  h_DeltaResUz_Lyr2->GetXaxis()->SetTitle("z (cm)");
167  h_DeltaResUz_Lyr2->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
168  h_DeltaResUz_Lyr3 = new TH2F("h_DeltaResUz_Lyr3", "Layer 3: #Delta res_{u} vs z", 250, -15, 20, 100, -1000, 1000);
169  h_DeltaResUz_Lyr3->GetXaxis()->SetTitle("z (cm)");
170  h_DeltaResUz_Lyr3->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
171  h_DeltaResUz_Lyr4 = new TH2F("h_DeltaResUz_Lyr4", "Layer 4: #Delta res_{u} vs z", 250, -20, 25, 100, -1000, 1000);
172  h_DeltaResUz_Lyr4->GetXaxis()->SetTitle("z (cm)");
173  h_DeltaResUz_Lyr4->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
174  h_DeltaResUz_Lyr5 = new TH2F("h_DeltaResUz_Lyr5", "Layer 5: #Delta res_{u} vs z", 250, -25, 35, 100, -1000, 1000);
175  h_DeltaResUz_Lyr5->GetXaxis()->SetTitle("z (cm)");
176  h_DeltaResUz_Lyr5->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
177  h_DeltaResUz_Lyr6 = new TH2F("h_DeltaResUz_Lyr6", "Layer 6: #Delta res_{u} vs z", 250, -30, 45, 100, -1000, 1000);
178  h_DeltaResUz_Lyr6->GetXaxis()->SetTitle("z (cm)");
179  h_DeltaResUz_Lyr6->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
180  //Define 2D histograms: difference of u-residuals vs phi of VXD overlaps for each layer (1 to 6)
181  h_DeltaResVPhi_Lyr1 = new TH2F("h_DeltaResVPhi_Lyr1", "Layer 1: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
182  h_DeltaResVPhi_Lyr1->GetXaxis()->SetTitle("#phi(rad)");
183  h_DeltaResVPhi_Lyr1->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
184  h_DeltaResVPhi_Lyr2 = new TH2F("h_DeltaResVPhi_Lyr2", "Layer 2: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
185  h_DeltaResVPhi_Lyr2->GetXaxis()->SetTitle("#phi(rad)");
186  h_DeltaResVPhi_Lyr2->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
187  h_DeltaResVPhi_Lyr3 = new TH2F("h_DeltaResVPhi_Lyr3", "Layer 3: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
188  h_DeltaResVPhi_Lyr3->GetXaxis()->SetTitle("#phi(rad)");
189  h_DeltaResVPhi_Lyr3->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
190  h_DeltaResVPhi_Lyr4 = new TH2F("h_DeltaResVPhi_Lyr4", "Layer 4: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
191  h_DeltaResVPhi_Lyr4->GetXaxis()->SetTitle("#phi(rad)");
192  h_DeltaResVPhi_Lyr4->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
193  h_DeltaResVPhi_Lyr5 = new TH2F("h_DeltaResVPhi_Lyr5", "Layer 5: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
194  h_DeltaResVPhi_Lyr5->GetXaxis()->SetTitle("#phi(rad)");
195  h_DeltaResVPhi_Lyr5->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
196  h_DeltaResVPhi_Lyr6 = new TH2F("h_DeltaResVPhi_Lyr6", "Layer 6: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
197  h_DeltaResVPhi_Lyr6->GetXaxis()->SetTitle("#phi(rad)");
198  h_DeltaResVPhi_Lyr6->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
199  //Define 2D histograms: difference of v-residuals vs z of VXD overlaps for each layer (1 to 6)
200  h_DeltaResVz_Lyr1 = new TH2F("h_DeltaResVz_Lyr1", "Layer 1: #Delta res_{v} vs z", 100, -4, 8, 100, -200, 200);
201  h_DeltaResVz_Lyr1->GetXaxis()->SetTitle("z (cm)");
202  h_DeltaResVz_Lyr1->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
203  h_DeltaResVz_Lyr2 = new TH2F("h_DeltaResVz_Lyr2", "Layer 2: #Delta res_{v} vs z", 100, -10, 15, 100, -200, 200);
204  h_DeltaResVz_Lyr2->GetXaxis()->SetTitle("z (cm)");
205  h_DeltaResVz_Lyr2->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
206  h_DeltaResVz_Lyr3 = new TH2F("h_DeltaResVz_Lyr3", "Layer 3: #Delta res_{v} vs z", 250, -15, 20, 100, -1000, 1000);
207  h_DeltaResVz_Lyr3->GetXaxis()->SetTitle("z (cm)");
208  h_DeltaResVz_Lyr3->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
209  h_DeltaResVz_Lyr4 = new TH2F("h_DeltaResVz_Lyr4", "Layer 4: #Delta res_{v} vs z", 250, -20, 25, 100, -1000, 1000);
210  h_DeltaResVz_Lyr4->GetXaxis()->SetTitle("z (cm)");
211  h_DeltaResVz_Lyr4->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
212  h_DeltaResVz_Lyr5 = new TH2F("h_DeltaResVz_Lyr5", "Layer 5: #Delta res_{v} vs z", 250, -25, 35, 100, -1000, 1000);
213  h_DeltaResVz_Lyr5->GetXaxis()->SetTitle("z (cm)");
214  h_DeltaResVz_Lyr5->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
215  h_DeltaResVz_Lyr6 = new TH2F("h_DeltaResVz_Lyr6", "Layer 6: #Delta res_{v} vs z", 250, -30, 45, 100, -1000, 1000);
216  h_DeltaResVz_Lyr6->GetXaxis()->SetTitle("z (cm)");
217  h_DeltaResVz_Lyr6->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
218  //Restricting to SVD clusters sizes
219  for (int i = 1; i < 5; i++) {
220  //The name is the product of cluster sizes for 2 consecutive hits (maximum size considered is 2)
221  TString h_name_U = "h_U_Cl1Cl2_" + std::to_string(i);
222  TString h_name_V = "h_V_Cl1Cl2_" + std::to_string(i);
223  TString title_U = "#Delta res_{u}: SVDClusterSize_1 x SVDClusterSize_2 = " + std::to_string(i);
224  TString title_V = "#Delta res_{v}: SVDClusterSize_1 x SVDClusterSize_2 = " + std::to_string(i);
225  h_U_Cl1Cl2_DeltaRes[i] = new TH1F(h_name_U, title_U, 100, -1000, 1000);
226  h_U_Cl1Cl2_DeltaRes[i]->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
227  h_U_Cl1Cl2_DeltaRes[i]->GetYaxis()->SetTitle("counts");
228  h_V_Cl1Cl2_DeltaRes[i] = new TH1F(h_name_V, title_V, 100, -1000, 1000);
229  h_V_Cl1Cl2_DeltaRes[i]->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
230  h_V_Cl1Cl2_DeltaRes[i]->GetYaxis()->SetTitle("counts");
231  }
232 
233  //Special case of ExpertLevel option enabled
234  if (m_ExpertLevel) {
235  //Create directory to store PXD and SVD hitmaps for overlapping hits
236  TDirectory* HMDir = oldDir->mkdir("HitMaps_VXDOverlaps");
237  HMDir->cd();
238  //Define 2D sensor hit-maps for reconstructed hits
239  for (int i = 1; i <= 5; i++) {
240  for (int j = 1; j <= 16; j++) {
241  TString h_name = "h_6" + std::to_string(j) + std::to_string(i);
242  TString title = "Layer:Ladder:Sensor = 6:" + std::to_string(j) + ":" + std::to_string(i);
243  h_Lyr6[j][i] = new TH2F(h_name, title, 100, -2.88, 2.88, 100, -6.14, 6.14);
244  h_Lyr6[j][i]->GetXaxis()->SetTitle("u (cm)");
245  h_Lyr6[j][i]->GetYaxis()->SetTitle("v (cm)");
246  }
247  }
248  for (int i = 1; i <= 4; i++) {
249  for (int j = 1; j <= 12; j++) {
250  TString h_name = "h_5" + std::to_string(j) + std::to_string(i);
251  TString title = "Layer:Ladder:Sensor = 5:" + std::to_string(j) + ":" + std::to_string(i);
252  h_Lyr5[j][i] = new TH2F(h_name, title, 100, -2.88, 2.88, 100, -6.14, 6.14);
253  h_Lyr5[j][i]->GetXaxis()->SetTitle("u (cm)");
254  h_Lyr5[j][i]->GetYaxis()->SetTitle("v (cm)");
255  }
256  }
257  for (int i = 1; i <= 3; i++) {
258  for (int j = 1; j <= 10; j++) {
259  TString h_name = "h_4" + std::to_string(j) + std::to_string(i);
260  TString title = "Layer:Ladder:Sensor = 4:" + std::to_string(j) + ":" + std::to_string(i);
261  h_Lyr4[j][i] = new TH2F(h_name, title, 100, -2.88, 2.88, 100, -6.14, 6.14);
262  h_Lyr4[j][i]->GetXaxis()->SetTitle("u (cm)");
263  h_Lyr4[j][i]->GetYaxis()->SetTitle("v (cm)");
264  }
265  }
266  for (int i = 1; i <= 2; i++) {
267  for (int j = 1; j <= 7; j++) {
268  TString h_name = "h_3" + std::to_string(j) + std::to_string(i);
269  TString title = "Layer:Ladder:Sensor = 3:" + std::to_string(j) + ":" + std::to_string(i);
270  h_Lyr3[j][i] = new TH2F(h_name, title, 100, -1.92, 1.92, 100, -6.14, 6.14);
271  h_Lyr3[j][i]->GetXaxis()->SetTitle("u (cm)");
272  h_Lyr3[j][i]->GetYaxis()->SetTitle("v (cm)");
273  }
274  }
275  for (int i = 1; i <= 2; i++) {
276  for (int j = 1; j <= 12; j++) {
277  TString h_name = "h_2" + std::to_string(j) + std::to_string(i);
278  TString title = "Layer:Ladder:Sensor = 2:" + std::to_string(j) + ":" + std::to_string(i);
279  h_Lyr2[j][i] = new TH2F(h_name, title, 100, -0.625, 0.625, 100, -3.072, 3.072);
280  h_Lyr2[j][i]->GetXaxis()->SetTitle("u (cm)");
281  h_Lyr2[j][i]->GetYaxis()->SetTitle("v (cm)");
282  }
283  }
284  for (int i = 1; i <= 2; i++) {
285  for (int j = 1; j <= 8; j++) {
286  TString h_name = "h_1" + std::to_string(j) + std::to_string(i);
287  TString title = "Layer:Ladder:Sensor = 1:" + std::to_string(j) + ":" + std::to_string(i);
288  h_Lyr1[j][i] = new TH2F(h_name, title, 100, -0.625, 0.625, 100, -2.24, 2.24);
289  h_Lyr1[j][i]->GetXaxis()->SetTitle("u (cm)");
290  h_Lyr1[j][i]->GetYaxis()->SetTitle("v (cm)");
291  }
292  }
293  //Create directory to store PXD and SVD trees
294  TDirectory* TreeDir = oldDir->mkdir("Trees_VXDOverlaps");
295  TreeDir->cd();
296  //Tree for PXD
297  t_PXD = new TTree("t_PXD", "Tree for PXD overlaps");
298  t_PXD->Branch("deltaResU_PXD", &deltaResU_PXD, "deltaResU_PXD/F");
299  t_PXD->Branch("intResU_PXD", &intResU_PXD, "intResU_PXD/F");
300  t_PXD->Branch("intResV_PXD", &intResV_PXD, "intResV_PXD/F");
301  t_PXD->Branch("intU_PXD", &intU_PXD, "intU_PXD/F");
302  t_PXD->Branch("intV_PXD", &intV_PXD, "intV_PXD/F");
303  t_PXD->Branch("intPhi_PXD", &intPhi_PXD, "intPhi_PXD/F");
304  t_PXD->Branch("intZ_PXD", &intZ_PXD, "intZ_PXD/F");
305  t_PXD->Branch("intLayer_PXD", &intLayer_PXD, "intLayer_PXD/i");
306  t_PXD->Branch("intLadder_PXD", &intLadder_PXD, "intLadder_PXD/i");
307  t_PXD->Branch("intSensor_PXD", &intSensor_PXD, "intSensor_PXD/i");
308  t_PXD->Branch("extResU_PXD", &extResU_PXD, "extResU_PXD/F");
309  t_PXD->Branch("extResV_PXD", &extResV_PXD, "extResV_PXD/F");
310  t_PXD->Branch("extU_PXD", &extU_PXD, "extU_PXD/F");
311  t_PXD->Branch("extV_PXD", &extV_PXD, "extV_PXD/F");
312  t_PXD->Branch("extPhi_PXD", &extPhi_PXD, "extPhi_PXD/F");
313  t_PXD->Branch("extZ_PXD", &extZ_PXD, "extZ_PXD/F");
314  t_PXD->Branch("extLayer_PXD", &extLayer_PXD, "extLayer_PXD/i");
315  t_PXD->Branch("extLadder_PXD", &extLadder_PXD, "extLadder_PXD/i");
316  t_PXD->Branch("extSensor_PXD", &extSensor_PXD, "extSensor_PXD/i");
317  //Tree for SVD u overlapping clusters
318  t_SVD_U = new TTree("t_SVD_U", "Tree for SVD u-overlaps");
319  t_SVD_U->Branch("deltaRes", &deltaRes_SVD_U, "deltaResU/F");
320  t_SVD_U->Branch("intRes", &intRes_SVD_U, "intRes/F");
321  t_SVD_U->Branch("intClPos", &intClPos_SVD_U, "intClPos/F");
322  t_SVD_U->Branch("intClPosErr", &intClPosErr_SVD_U, "intClPosErr/F");
323  t_SVD_U->Branch("intTruePos", &intTruePos_SVD_U, "intTruePos/F");
324  t_SVD_U->Branch("intClPhi", &intClPhi_SVD_U, "intClPhi/F");
325  t_SVD_U->Branch("intClZ", &intClZ_SVD_U, "intClZ/F");
326  t_SVD_U->Branch("intTrkPos", &intTrkPos_SVD_U, "intTrkPos/F");
327  t_SVD_U->Branch("intTrkPosErr", &intTrkPosErr_SVD_U, "intTrkPosErr/F");
328  t_SVD_U->Branch("intTrkQoP", &intTrkQoP_SVD_U, "intTrkQoP/F");
329  t_SVD_U->Branch("intTrkPrime", &intTrkPrime_SVD_U, "intTrkPrime/F");
330  t_SVD_U->Branch("intLayer", &intLayer_SVD_U, "intLayer/i");
331  t_SVD_U->Branch("intLadder", &intLadder_SVD_U, "intLadder/i");
332  t_SVD_U->Branch("intSensor", &intSensor_SVD_U, "intSensor/i");
333  t_SVD_U->Branch("intSize", &intSize_SVD_U, "intSize/i");
334  t_SVD_U->Branch("deltaRes", &deltaRes_SVD_U, "deltaResU/F");
335  t_SVD_U->Branch("extRes", &extRes_SVD_U, "extRes/F");
336  t_SVD_U->Branch("extClPos", &extClPos_SVD_U, "extClPos/F");
337  t_SVD_U->Branch("extClPosErr", &extClPosErr_SVD_U, "extClPosErr/F");
338  t_SVD_U->Branch("extTruePos", &extTruePos_SVD_U, "extTruePos/F");
339  t_SVD_U->Branch("extClPhi", &extClPhi_SVD_U, "extClPhi/F");
340  t_SVD_U->Branch("extClZ", &extClZ_SVD_U, "extClZ/F");
341  t_SVD_U->Branch("extTrkPos", &extTrkPos_SVD_U, "extTrkPos/F");
342  t_SVD_U->Branch("extTrkPosErr", &extTrkPosErr_SVD_U, "extTrkPosErr/F");
343  t_SVD_U->Branch("extTrkQoP", &extTrkQoP_SVD_U, "extTrkQoP/F");
344  t_SVD_U->Branch("extTrkPrime", &extTrkPrime_SVD_U, "extTrkPrime/F");
345  t_SVD_U->Branch("extLayer", &extLayer_SVD_U, "extLayer/i");
346  t_SVD_U->Branch("extLadder", &extLadder_SVD_U, "extLadder/i");
347  t_SVD_U->Branch("extSensor", &extSensor_SVD_U, "extSensor/i");
348  t_SVD_U->Branch("extSize", &extSize_SVD_U, "extSize/i");
349  //Tree for SVD v overlapping clusters
350  t_SVD_V = new TTree("t_SVD_V", "Tree for SVD v-overlaps");
351  t_SVD_V->Branch("deltaRes", &deltaRes_SVD_V, "deltaResV/F");
352  t_SVD_V->Branch("intRes", &intRes_SVD_V, "intRes/F");
353  t_SVD_V->Branch("intClPos", &intClPos_SVD_V, "intClPos/F");
354  t_SVD_V->Branch("intClPosErr", &intClPosErr_SVD_V, "intClPosErr/F");
355  t_SVD_V->Branch("intTruePos", &intTruePos_SVD_V, "intTruePos/F");
356  t_SVD_V->Branch("intClPhi", &intClPhi_SVD_V, "intClPhi/F");
357  t_SVD_V->Branch("intClZ", &intClZ_SVD_V, "intClZ/F");
358  t_SVD_V->Branch("intTrkPos", &intTrkPos_SVD_V, "intTrkPos/F");
359  t_SVD_V->Branch("intTrkPosErr", &intTrkPosErr_SVD_V, "intTrkPosErr/F");
360  t_SVD_V->Branch("intTrkQoP", &intTrkQoP_SVD_V, "intTrkQoP/F");
361  t_SVD_V->Branch("intTrkPrime", &intTrkPrime_SVD_V, "intTrkPrime/F");
362  t_SVD_V->Branch("intLayer", &intLayer_SVD_V, "intLayer/i");
363  t_SVD_V->Branch("intLadder", &intLadder_SVD_V, "intLadder/i");
364  t_SVD_V->Branch("intSensor", &intSensor_SVD_V, "intSensor/i");
365  t_SVD_V->Branch("intSize", &intSize_SVD_V, "intSize/i");
366  t_SVD_V->Branch("deltaRes", &deltaRes_SVD_V, "deltaResV/F");
367  t_SVD_V->Branch("extRes", &extRes_SVD_V, "extRes/F");
368  t_SVD_V->Branch("extClPos", &extClPos_SVD_V, "extClPos/F");
369  t_SVD_V->Branch("extClPosErr", &extClPosErr_SVD_V, "extClPosErr/F");
370  t_SVD_V->Branch("extTruePos", &extTruePos_SVD_V, "extTruePos/F");
371  t_SVD_V->Branch("extClPhi", &extClPhi_SVD_V, "extClPhi/F");
372  t_SVD_V->Branch("extClZ", &extClZ_SVD_V, "extClZ/F");
373  t_SVD_V->Branch("extTrkPos", &extTrkPos_SVD_V, "extTrkPos/F");
374  t_SVD_V->Branch("extTrkPosErr", &extTrkPosErr_SVD_V, "extTrkPosErr/F");
375  t_SVD_V->Branch("extTrkQoP", &extTrkQoP_SVD_V, "extTrkQoP/F");
376  t_SVD_V->Branch("extTrkPrime", &extTrkPrime_SVD_V, "extTrkPrime/F");
377  t_SVD_V->Branch("extLayer", &extLayer_SVD_V, "extLayer/i");
378  t_SVD_V->Branch("extLadder", &extLadder_SVD_V, "extLadder/i");
379  t_SVD_V->Branch("extSensor", &extSensor_SVD_V, "extSensor/i");
380  t_SVD_V->Branch("extSize", &extSize_SVD_V, "extSize/i");
381  }
382  //Go back to the default output directory
383  oldDir->cd();
384 }
385 
386 void OverlapResidualsModule::event()
387 {
388  bool isMC = Environment::Instance().isMC();
389 
390  static VXD::GeoCache& geo = VXD::GeoCache::getInstance();
391  StoreArray<RecoTrack> recoTracks(m_recoTracksStoreArrayName);
392  for (const auto& trk : recoTracks) {
393  if (! trk.wasFitSuccessful()) {
394  continue;
395  }
396  const vector<PXDCluster* > pxdClusters = trk.getPXDHitList();
397  const vector<SVDCluster* > svdClusters = trk.getSVDHitList();
398  B2DEBUG(40, "FITTED TRACK: NUMBER OF PXD HITS = " << pxdClusters.size() << " NUMBER OF SVD HITS = " << svdClusters.size());
399  //LOOKING FOR 2 CONSECUTIVE PXD HITS IN OVERLAPPING MODULES OF A SAME LAYER
400  for (unsigned int i = 0; i < pxdClusters.size(); i++) {
401  const PXDCluster* pxd_1 = pxdClusters[i];
402  const RecoHitInformation* infoPXD_1 = trk.getRecoHitInformation(pxd_1);
403  if (!infoPXD_1) {
404  continue;
405  }
406  const auto* hitTrackPoint_1 = trk.getCreatedTrackPoint(infoPXD_1);
407  const auto* fittedResult_1 = hitTrackPoint_1->getFitterInfo();
408  if (!fittedResult_1) {
409  continue;
410  }
411  const VxdID pxd_id_1 = pxd_1->getSensorID();
412  const unsigned short pxd_Layer_1 = pxd_id_1.getLayerNumber();
413  const unsigned short pxd_Ladder_1 = pxd_id_1.getLadderNumber();
414  const unsigned short pxd_Sensor_1 = pxd_id_1.getSensorNumber();
415  for (unsigned int l = i + 1; l < pxdClusters.size(); l++) {
416  const PXDCluster* pxd_2 = pxdClusters[l];
417  const RecoHitInformation* infoPXD_2 = trk.getRecoHitInformation(pxd_2);
418  if (!infoPXD_2) {
419  continue;
420  }
421  const auto* hitTrackPoint_2 = trk.getCreatedTrackPoint(infoPXD_2);
422  const auto* fittedResult_2 = hitTrackPoint_2->getFitterInfo();
423  if (!fittedResult_2) {
424  continue;
425  }
426  const VxdID pxd_id_2 = pxd_2->getSensorID();
427  const unsigned short pxd_Layer_2 = pxd_id_2.getLayerNumber();
428  const unsigned short pxd_Ladder_2 = pxd_id_2.getLadderNumber();
429  const unsigned short pxd_Sensor_2 = pxd_id_2.getSensorNumber();
430  if (pxd_Layer_1 == pxd_Layer_2 && ((pxd_Ladder_2 - pxd_Ladder_1) == -1.0 || (pxd_Ladder_2 - pxd_Ladder_1) == 7.0
431  || (pxd_Ladder_2 - pxd_Ladder_1) == 11.0)) {
432  B2DEBUG(40, " ============= 2 hits in a PXD overlap ============= ");
433  const TVectorD resUnBias_PXD_1 = fittedResult_1->getResidual(0, false).getState();
434  const TVectorD resUnBias_PXD_2 = fittedResult_2->getResidual(0, false).getState();
435  const float res_U_1 = resUnBias_PXD_1.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
436  const float res_V_1 = resUnBias_PXD_1.GetMatrixArray()[1] * Unit::convertValueToUnit(1.0, "um");
437  const double res_U_2 = resUnBias_PXD_2.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
438  const double res_V_2 = resUnBias_PXD_2.GetMatrixArray()[1] * Unit::convertValueToUnit(1.0, "um");
439  const float over_U_PXD = res_U_2 - res_U_1;
440  const float over_V_PXD = res_V_2 - res_V_1;
441  const TVector3 pxdLocal_1(pxd_1->getU(), pxd_1->getV(), 0.);
442  const TVector3 pxdLocal_2(pxd_2->getU(), pxd_2->getV(), 0.);
443  const VXD::SensorInfoBase& pxdSensor_1 = geo.get(pxd_id_1);
444  const VXD::SensorInfoBase& pxdSensor_2 = geo.get(pxd_id_2);
445  const TVector3& pxdGlobal_1 = pxdSensor_1.pointToGlobal(pxdLocal_1);
446  const TVector3& pxdGlobal_2 = pxdSensor_2.pointToGlobal(pxdLocal_2);
447  double pxdPhi_1 = atan2(pxdGlobal_1(1), pxdGlobal_1(0));
448  double pxdPhi_2 = atan2(pxdGlobal_2(1), pxdGlobal_2(0));
449  double pxdZ_1 = pxdGlobal_1(2);
450  double pxdZ_2 = pxdGlobal_2(2);
451  B2DEBUG(40, "PXD: difference of residuals " << over_U_PXD << " " << over_V_PXD);
452  //Fill PXD tree for overlaps if required by the user
453  if (m_ExpertLevel) {
454  deltaResU_PXD = over_U_PXD;
455  intResU_PXD = res_U_1;
456  intResV_PXD = res_V_1;
457  intU_PXD = pxd_1->getU();
458  intV_PXD = pxd_1->getV();
459  intPhi_PXD = pxdPhi_1;
460  intZ_PXD = pxdZ_1;
461  intLayer_PXD = pxd_Layer_1;
462  intLadder_PXD = pxd_Ladder_1;
463  intSensor_PXD = pxd_Sensor_1;
464  extResU_PXD = res_U_2;
465  extResV_PXD = res_V_2;
466  extU_PXD = pxd_2->getU();
467  extV_PXD = pxd_2->getV();
468  extPhi_PXD = pxdPhi_2;
469  extZ_PXD = pxdZ_2;
470  extLayer_PXD = pxd_Layer_2;
471  extLadder_PXD = pxd_Ladder_2;
472  extSensor_PXD = pxd_Sensor_2;
473  t_PXD->Fill();
474  }
475  //Fill histograms of residuals differences with PXD clusters
476  h_U_DeltaRes->Fill(over_U_PXD);
477  h_V_DeltaRes->Fill(over_V_PXD);
478  h_U_DeltaRes_PXD->Fill(over_U_PXD);
479  h_V_DeltaRes_PXD->Fill(over_V_PXD);
480  //Fill sensor hit-maps and 2D histograms with PXD clusters
481  if (pxd_Layer_1 == 1 && pxd_Layer_2 == 1) {
482  h_Lyr1[pxd_Ladder_1][pxd_Sensor_1]->Fill(pxd_1->getU(), pxd_1->getV());
483  h_Lyr1[pxd_Ladder_2][pxd_Sensor_2]->Fill(pxd_2->getU(), pxd_2->getV());
484  h_U_DeltaRes_PXD_Lyr1->Fill(over_U_PXD);
485  h_V_DeltaRes_PXD_Lyr1->Fill(over_V_PXD);
486  h_DeltaResUPhi_Lyr1->Fill(pxdPhi_1, over_U_PXD);
487  h_DeltaResUPhi_Lyr1->Fill(pxdPhi_2, over_U_PXD);
488  h_DeltaResVPhi_Lyr1->Fill(pxdPhi_1, over_V_PXD);
489  h_DeltaResVPhi_Lyr1->Fill(pxdPhi_2, over_V_PXD);
490  h_DeltaResVz_Lyr1->Fill(pxdZ_1, over_V_PXD);
491  h_DeltaResVz_Lyr1->Fill(pxdZ_2, over_V_PXD);
492  h_DeltaResUz_Lyr1->Fill(pxdZ_1, over_U_PXD);
493  h_DeltaResUz_Lyr1->Fill(pxdZ_2, over_U_PXD);
494  }
495  if (pxd_Layer_1 == 2 && pxd_Layer_2 == 2) {
496  h_Lyr2[pxd_Ladder_1][pxd_Sensor_1]->Fill(pxd_1->getU(), pxd_1->getV());
497  h_Lyr2[pxd_Ladder_2][pxd_Sensor_2]->Fill(pxd_2->getU(), pxd_2->getV());
498  h_U_DeltaRes_PXD_Lyr2->Fill(over_U_PXD);
499  h_V_DeltaRes_PXD_Lyr2->Fill(over_V_PXD);
500  h_DeltaResUPhi_Lyr2->Fill(pxdPhi_1, over_U_PXD);
501  h_DeltaResUPhi_Lyr2->Fill(pxdPhi_2, over_U_PXD);
502  h_DeltaResVPhi_Lyr2->Fill(pxdPhi_1, over_V_PXD);
503  h_DeltaResVPhi_Lyr2->Fill(pxdPhi_2, over_V_PXD);
504  h_DeltaResVz_Lyr2->Fill(pxdZ_1, over_V_PXD);
505  h_DeltaResVz_Lyr2->Fill(pxdZ_2, over_V_PXD);
506  h_DeltaResUz_Lyr2->Fill(pxdZ_1, over_U_PXD);
507  h_DeltaResUz_Lyr2->Fill(pxdZ_2, over_U_PXD);
508  }
509  }
510  }
511 
512  }
513 
514  //LOOKING FOR 2 CONSECUTIVE SVD HITS IN OVERLAPPING MODULES OF A SAME LAYER
515  for (unsigned int i = 0; i < svdClusters.size(); i++) {
516  const SVDCluster* svd_1 = svdClusters[i];
517 
518  //get true hits, used only if isMC
519  RelationVector<SVDTrueHit> trueHit_1 = DataStore::getRelationsWithObj<SVDTrueHit>(svd_1);
520 
521  const RecoHitInformation* infoSVD_1 = trk.getRecoHitInformation(svd_1);
522  if (!infoSVD_1) {
523  continue;
524  }
525  const auto* hitTrackPoint_1 = trk.getCreatedTrackPoint(infoSVD_1);
526  const auto* fittedResult_1 = hitTrackPoint_1->getFitterInfo();
527  if (!fittedResult_1) {
528  continue;
529  }
530  const VxdID svd_id_1 = svd_1->getSensorID();
531  const unsigned short svd_Layer_1 = svd_id_1.getLayerNumber();
532  const unsigned short svd_Ladder_1 = svd_id_1.getLadderNumber();
533  const unsigned short svd_Sensor_1 = svd_id_1.getSensorNumber();
534  for (unsigned int l = i + 1; l < svdClusters.size(); l++) {
535  const SVDCluster* svd_2 = svdClusters[l];
536 
537  //get true hits, used only if isMC
538  RelationVector<SVDTrueHit> trueHit_2 = DataStore::getRelationsWithObj<SVDTrueHit>(svd_2);
539 
540  const RecoHitInformation* infoSVD_2 = trk.getRecoHitInformation(svd_2);
541  if (!infoSVD_2) {
542  continue;
543  }
544  const auto* hitTrackPoint_2 = trk.getCreatedTrackPoint(infoSVD_2);
545  const auto* fittedResult_2 = hitTrackPoint_2->getFitterInfo();
546  if (!fittedResult_2) {
547  continue;
548  }
549  const VxdID svd_id_2 = svd_2->getSensorID();
550  const unsigned short svd_Layer_2 = svd_id_2.getLayerNumber();
551  const unsigned short svd_Ladder_2 = svd_id_2.getLadderNumber();
552  const unsigned short svd_Sensor_2 = svd_id_2.getSensorNumber();
553  if (svd_Layer_1 == svd_Layer_2 && ((svd_Ladder_2 - svd_Ladder_1) == 1.0 || (svd_Ladder_2 - svd_Ladder_1) == -6.0
554  || (svd_Ladder_2 - svd_Ladder_1) == -9.0 || (svd_Ladder_2 - svd_Ladder_1) == -11.0 || (svd_Ladder_2 - svd_Ladder_1) == -15.0)) {
555  B2DEBUG(40, " ============= 2 hits in a SVD overlap ============= ");
556  const TVectorD resUnBias_SVD_1 = fittedResult_1->getResidual(0, false).getState();
557  const TVectorD resUnBias_SVD_2 = fittedResult_2->getResidual(0, false).getState();
558  genfit::MeasuredStateOnPlane state_1 = trk.getMeasuredStateOnPlaneFromRecoHit(infoSVD_1);
559  const TVectorD& svd_predIntersect_1 = state_1.getState();
560  const TMatrixDSym& covMatrix_1 = state_1.getCov();
561  genfit::MeasuredStateOnPlane state_2 = trk.getMeasuredStateOnPlaneFromRecoHit(infoSVD_2);
562  const TVectorD& svd_predIntersect_2 = state_2.getState();
563  const TMatrixDSym& covMatrix_2 = state_2.getCov();
564  //Restricting to consecutive SVD u-clusters
565  if (svd_1->isUCluster() == true && svd_2->isUCluster() == true) {
566  const int strips_1 = svd_1->getSize();
567  h_SVDstrips_Mult->Fill(strips_1);
568  const int strips_2 = svd_2->getSize();
569  h_SVDstrips_Mult->Fill(strips_2);
570  const double res_U_1 = resUnBias_SVD_1.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
571  const double res_U_2 = resUnBias_SVD_2.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
572  const float over_U_SVD = res_U_2 - res_U_1;
573  const TVector3 svdLocal_1(svd_1->getPosition(), svd_predIntersect_1[4], 0.);
574  const TVector3 svdLocal_2(svd_2->getPosition(), svd_predIntersect_2[4], 0.);
575  const VXD::SensorInfoBase& svdSensor_1 = geo.get(svd_id_1);
576  const VXD::SensorInfoBase& svdSensor_2 = geo.get(svd_id_2);
577  const TVector3& svdGlobal_1 = svdSensor_1.pointToGlobal(svdLocal_1);
578  const TVector3& svdGlobal_2 = svdSensor_2.pointToGlobal(svdLocal_2);
579  double svdPhi_1 = atan2(svdGlobal_1(1), svdGlobal_1(0));
580  double svdPhi_2 = atan2(svdGlobal_2(1), svdGlobal_2(0));
581  double svdZ_1 = svdGlobal_1(2);
582  double svdZ_2 = svdGlobal_2(2);
583  B2DEBUG(40, "SVD: difference of u-residuals =========> " << over_U_SVD);
584  //Fill SVD tree for u-overlaps if required by the user
585  if (m_ExpertLevel) {
586  deltaRes_SVD_U = over_U_SVD;
587  intRes_SVD_U = res_U_1;
588  intClPos_SVD_U = svd_1->getPosition();
589  intClPosErr_SVD_U = svd_1->getPositionSigma();
590  if (isMC && trueHit_1.size() > 0)
591  intTruePos_SVD_U = trueHit_1[0]->getU();
592  else
593  intTruePos_SVD_U = -99;
594  intClPhi_SVD_U = svdPhi_1;
595  intClZ_SVD_U = svdZ_1;
596  intTrkPos_SVD_U = svd_predIntersect_1[3];
597  intTrkPosErr_SVD_U = sqrt(covMatrix_1[3][3]);
598  intTrkQoP_SVD_U = svd_predIntersect_1[0];
599  intTrkPrime_SVD_U = svd_predIntersect_1[1];
600  intLayer_SVD_U = svd_Layer_1;
601  intLadder_SVD_U = svd_Ladder_1;
602  intSensor_SVD_U = svd_Sensor_1;
603  intSize_SVD_U = strips_1;
604  extRes_SVD_U = res_U_2;
605  extClPos_SVD_U = svd_2->getPosition();
606  extClPosErr_SVD_U = svd_2->getPositionSigma();
607  if (isMC && trueHit_2.size() > 0)
608  extTruePos_SVD_U = trueHit_2[0]->getU();
609  else
610  extTruePos_SVD_U = -99;
611  extClPhi_SVD_U = svdPhi_2;
612  extClZ_SVD_U = svdZ_2;
613  extTrkPos_SVD_U = svd_predIntersect_2[3];
614  extTrkPosErr_SVD_U = sqrt(covMatrix_2[3][3]);
615  extTrkQoP_SVD_U = svd_predIntersect_2[0];
616  extTrkPrime_SVD_U = svd_predIntersect_2[1];
617  extLayer_SVD_U = svd_Layer_2;
618  extLadder_SVD_U = svd_Ladder_2;
619  extSensor_SVD_U = svd_Sensor_2;
620  extSize_SVD_U = strips_2;
621  t_SVD_U->Fill();
622  }
623  //Fill histograms of residuals differences with SVD u clusters
624  h_U_DeltaRes->Fill(over_U_SVD);
625  h_U_DeltaRes_SVD->Fill(over_U_SVD);
626  if (svd_1->getSize() < 3 && svd_2->getSize() < 3) { //Consider only clusters with 3 or less strips involved
627  h_U_Cl1Cl2_DeltaRes[svd_1->getSize() * svd_2->getSize()]->Fill(over_U_SVD);
628  }
629  //Fill sensor hit-maps and 2D histograms with SVD u clusters
630  if (svd_Layer_1 == 3 && svd_Layer_2 == 3) {
631  h_Lyr3[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
632  h_Lyr3[svd_Ladder_2][svd_Sensor_2]->Fill(svd_2->getPosition(), 0.);
633  h_U_DeltaRes_SVD_Lyr3->Fill(over_U_SVD);
634  h_DeltaResUPhi_Lyr3->Fill(svdPhi_1, over_U_SVD);
635  h_DeltaResUPhi_Lyr3->Fill(svdPhi_2, over_U_SVD);
636  h_DeltaResUz_Lyr3->Fill(svdZ_1, over_U_SVD);
637  h_DeltaResUz_Lyr3->Fill(svdZ_2, over_U_SVD);
638  }
639  if (svd_Layer_1 == 4 && svd_Layer_2 == 4) {
640  h_Lyr4[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
641  h_Lyr4[svd_Ladder_2][svd_Sensor_2]->Fill(svd_2->getPosition(), 0.);
642  h_U_DeltaRes_SVD_Lyr4->Fill(over_U_SVD);
643  h_DeltaResUPhi_Lyr4->Fill(svdPhi_1, over_U_SVD);
644  h_DeltaResUPhi_Lyr4->Fill(svdPhi_2, over_U_SVD);
645  h_DeltaResUz_Lyr4->Fill(svdZ_1, over_U_SVD);
646  h_DeltaResUz_Lyr4->Fill(svdZ_2, over_U_SVD);
647  }
648  if (svd_Layer_1 == 5 && svd_Layer_2 == 5) {
649  h_Lyr5[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
650  h_Lyr5[svd_Ladder_2][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
651  h_U_DeltaRes_SVD_Lyr5->Fill(over_U_SVD);
652  h_DeltaResUPhi_Lyr5->Fill(svdPhi_1, over_U_SVD);
653  h_DeltaResUPhi_Lyr5->Fill(svdPhi_2, over_U_SVD);
654  h_DeltaResUz_Lyr5->Fill(svdZ_1, over_U_SVD);
655  h_DeltaResUz_Lyr5->Fill(svdZ_2, over_U_SVD);
656  }
657  if (svd_Layer_1 == 6 && svd_Layer_2 == 6) {
658  h_Lyr6[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
659  h_Lyr6[svd_Ladder_2][svd_Sensor_2]->Fill(svd_2->getPosition(), 0.);
660  h_U_DeltaRes_SVD_Lyr6->Fill(over_U_SVD);
661  h_DeltaResUPhi_Lyr6->Fill(svdPhi_1, over_U_SVD);
662  h_DeltaResUPhi_Lyr6->Fill(svdPhi_2, over_U_SVD);
663  h_DeltaResUz_Lyr6->Fill(svdZ_1, over_U_SVD);
664  h_DeltaResUz_Lyr6->Fill(svdZ_2, over_U_SVD);
665  }
666  }
667  //Restricting to consecutive SVD v-clusters
668  if (svd_1->isUCluster() != true && svd_2->isUCluster() != true) {
669  const int strips_1 = svd_1->getSize();
670  h_SVDstrips_Mult->Fill(strips_1);
671  const int strips_2 = svd_2->getSize();
672  h_SVDstrips_Mult->Fill(strips_2);
673  const double res_V_1 = resUnBias_SVD_1.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
674  const double res_V_2 = resUnBias_SVD_2.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
675  const float over_V_SVD = res_V_2 - res_V_1;
676  const TVector3 svdLocal_1(svd_predIntersect_1[3], svd_1->getPosition(), 0.);
677  const TVector3 svdLocal_2(svd_predIntersect_2[3], svd_2->getPosition(), 0.);
678  const VXD::SensorInfoBase& svdSensor_1 = geo.get(svd_id_1);
679  const VXD::SensorInfoBase& svdSensor_2 = geo.get(svd_id_2);
680  const TVector3& svdGlobal_1 = svdSensor_1.pointToGlobal(svdLocal_1);
681  const TVector3& svdGlobal_2 = svdSensor_2.pointToGlobal(svdLocal_2);
682  double svdPhi_1 = atan2(svdGlobal_1(1), svdGlobal_1(0));
683  double svdPhi_2 = atan2(svdGlobal_2(1), svdGlobal_2(0));
684  double svdZ_1 = svdGlobal_1(2);
685  double svdZ_2 = svdGlobal_2(2);
686  B2DEBUG(40, "SVD: difference of v-residuals =========> " << over_V_SVD);
687  //Fill SVD tree for v-overlaps if required by the user
688  if (m_ExpertLevel) {
689  deltaRes_SVD_V = over_V_SVD;
690  intRes_SVD_V = res_V_1;
691  intClPos_SVD_V = svd_1->getPosition();
692  intClPosErr_SVD_V = svd_1->getPositionSigma();
693  if (isMC && trueHit_1.size() > 0)
694  intTruePos_SVD_V = trueHit_1[0]->getV();
695  else
696  intTruePos_SVD_V = -99;
697  intClPhi_SVD_V = svdPhi_1;
698  intClZ_SVD_V = svdZ_1;
699  intTrkPos_SVD_V = svd_predIntersect_1[4];
700  intTrkPosErr_SVD_V = sqrt(covMatrix_1[4][4]);
701  intTrkQoP_SVD_V = svd_predIntersect_1[0];
702  intTrkPrime_SVD_V = svd_predIntersect_1[2];
703  intLayer_SVD_V = svd_Layer_1;
704  intLadder_SVD_V = svd_Ladder_1;
705  intSensor_SVD_V = svd_Sensor_1;
706  intSize_SVD_V = strips_1;
707  extRes_SVD_V = res_V_2;
708  extClPos_SVD_V = svd_2->getPosition();
709  extClPosErr_SVD_V = svd_2->getPositionSigma();
710  if (isMC && trueHit_2.size() > 0)
711  extTruePos_SVD_V = trueHit_2[0]->getV();
712  else
713  extTruePos_SVD_V = -99;
714  extClPhi_SVD_V = svdPhi_2;
715  extClZ_SVD_V = svdZ_2;
716  extTrkPos_SVD_V = svd_predIntersect_2[4];
717  extTrkPosErr_SVD_V = sqrt(covMatrix_2[4][4]);
718  extTrkQoP_SVD_V = svd_predIntersect_2[0];
719  extTrkPrime_SVD_V = svd_predIntersect_2[2];
720  extLayer_SVD_V = svd_Layer_2;
721  extLadder_SVD_V = svd_Ladder_2;
722  extSensor_SVD_V = svd_Sensor_2;
723  extSize_SVD_V = strips_2;
724  t_SVD_V->Fill();
725  }
726  //Fill histograms of residuals differences with SVD v clusters
727  h_V_DeltaRes->Fill(over_V_SVD);
728  h_V_DeltaRes_SVD->Fill(over_V_SVD);
729  if (svd_1->getSize() < 3 && svd_2->getSize() < 3) { //Consider only clusters with 3 or less strips involved
730  h_V_Cl1Cl2_DeltaRes[svd_1->getSize() * svd_2->getSize()]->Fill(over_V_SVD);
731  }
732  //Fill sensor hit-maps and 2D histograms with SVD v clusters
733  if (svd_Layer_1 == 3 && svd_Layer_2 == 3) {
734  h_Lyr3[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
735  h_Lyr3[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
736  h_V_DeltaRes_SVD_Lyr3->Fill(over_V_SVD);
737  h_DeltaResVPhi_Lyr3->Fill(svdPhi_1, over_V_SVD);
738  h_DeltaResVPhi_Lyr3->Fill(svdPhi_2, over_V_SVD);
739  h_DeltaResVz_Lyr3->Fill(svdZ_1, over_V_SVD);
740  h_DeltaResVz_Lyr3->Fill(svdZ_2, over_V_SVD);
741  }
742  if (svd_Layer_1 == 4 && svd_Layer_2 == 4) {
743  h_Lyr4[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
744  h_Lyr4[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
745  h_V_DeltaRes_SVD_Lyr4->Fill(over_V_SVD);
746  h_DeltaResVPhi_Lyr4->Fill(svdPhi_1, over_V_SVD);
747  h_DeltaResVPhi_Lyr4->Fill(svdPhi_2, over_V_SVD);
748  h_DeltaResVz_Lyr4->Fill(svdZ_1, over_V_SVD);
749  h_DeltaResVz_Lyr4->Fill(svdZ_2, over_V_SVD);
750  }
751  if (svd_Layer_1 == 5 && svd_Layer_2 == 5) {
752  h_Lyr5[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
753  h_Lyr5[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
754  h_V_DeltaRes_SVD_Lyr5->Fill(over_V_SVD);
755  h_DeltaResVPhi_Lyr5->Fill(svdPhi_1, over_V_SVD);
756  h_DeltaResVPhi_Lyr5->Fill(svdPhi_2, over_V_SVD);
757  h_DeltaResVz_Lyr5->Fill(svdZ_1, over_V_SVD);
758  h_DeltaResVz_Lyr5->Fill(svdZ_2, over_V_SVD);
759  }
760  if (svd_Layer_1 == 6 && svd_Layer_2 == 6) {
761  h_Lyr6[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
762  h_Lyr6[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
763  h_V_DeltaRes_SVD_Lyr6->Fill(over_V_SVD);
764  h_DeltaResVPhi_Lyr6->Fill(svdPhi_1, over_V_SVD);
765  h_DeltaResVPhi_Lyr6->Fill(svdPhi_2, over_V_SVD);
766  h_DeltaResVz_Lyr6->Fill(svdZ_1, over_V_SVD);
767  h_DeltaResVz_Lyr6->Fill(svdZ_2, over_V_SVD);
768  }
769  }
770  }
771  }
772  }
773  }
774 }
775 
776 
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:98
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
Belle2::VXD::GeoCache::get
static const SensorInfoBase & get(Belle2::VxdID id)
Return a reference to the SensorInfo of a given SensorID.
Definition: GeoCache.h:141
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::PXDCluster::getU
float getU() const
Get u coordinate of hit position.
Definition: PXDCluster.h:136
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::VXD::SensorInfoBase
Base class to provide Sensor Information for PXD and SVD.
Definition: SensorInfoBase.h:40
Belle2::OverlapResidualsModule::defineHisto
void defineHisto() override
Histogram definitions such as TH1(), TH2(), TNtuple(), TTree()....
Definition: OverlapResidualsModule.cc:62
Belle2::Environment::isMC
bool isMC() const
Do we have generated, not real data?
Definition: Environment.cc:57
Belle2::OverlapResidualsModule::event
void event() override
Compute the difference of coordinate residuals between two hits in overlapping sensors of a same VXD ...
Definition: OverlapResidualsModule.cc:386
Belle2::SVDCluster::getPositionSigma
float getPositionSigma() const
Get the error of the reconstructed hit coordinate.
Definition: SVDCluster.h:137
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:215
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::RecoHitInformation
This class stores additional information to every CDC/SVD/PXD hit stored in a RecoTrack.
Definition: RecoHitInformation.h:48
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::Unit::convertValueToUnit
static double convertValueToUnit(double value, const std::string &unitString)
Converts a floating point value from the standard framework unit to the given unit.
Definition: UnitConst.cc:146
Belle2::VXD::SensorInfoBase::pointToGlobal
TVector3 pointToGlobal(const TVector3 &local, bool reco=false) const
Convert a point from local to global coordinates.
Definition: SensorInfoBase.h:373
Belle2::SVDCluster::isUCluster
bool isUCluster() const
Get the direction of strips.
Definition: SVDCluster.h:118
Belle2::OverlapResidualsModule::initialize
void initialize() override
Register input and output data.
Definition: OverlapResidualsModule.cc:54
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::OverlapResidualsModule
The module studies VXD hits from overlapping sensors of a same VXD layer.
Definition: OverlapResidualsModule.h:46
Belle2::SVDCluster::getSensorID
VxdID getSensorID() const
Get the sensor ID.
Definition: SVDCluster.h:110
Belle2::SVDCluster::getPosition
float getPosition(double v=0) const
Get the coordinate of reconstructed hit.
Definition: SVDCluster.h:125
Belle2::PXDCluster::getSensorID
VxdID getSensorID() const
Get the sensor ID.
Definition: PXDCluster.h:131
Belle2::SVDCluster::getSize
unsigned short getSize() const
Get cluster size.
Definition: SVDCluster.h:162
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
Belle2::VXD::GeoCache
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
Definition: GeoCache.h:41
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::Environment::Instance
static Environment & Instance()
Static method to get a reference to the Environment instance.
Definition: Environment.cc:31
Belle2::PXDCluster::getV
float getV() const
Get v coordinate of hit position.
Definition: PXDCluster.h:141
Belle2::HistoModule
HistoModule.h is supposed to be used instead of Module.h for the modules with histogram definitions t...
Definition: HistoModule.h:29