Belle II Software  release-06-02-00
OverlapResidualsModule.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <tracking/modules/VXDOverlaps/OverlapResidualsModule.h>
10 #include <framework/datastore/StoreArray.h>
11 #include <framework/gearbox/Unit.h>
12 #include <framework/core/Environment.h>
13 #include <svd/dataobjects/SVDCluster.h>
14 #include <svd/dataobjects/SVDRecoDigit.h>
15 #include <svd/dataobjects/SVDShaperDigit.h>
16 #include <svd/dataobjects/SVDTrueHit.h>
17 #include <vxd/dataobjects/VxdID.h>
18 #include <vxd/geometry/GeoCache.h>
19 #include <vxd/geometry/SensorInfoBase.h>
20 #include <tracking/dataobjects/RecoTrack.h>
21 #include <tracking/dataobjects/RecoHitInformation.h>
22 #include <genfit/TrackPoint.h>
23 #include <TVector3.h>
24 #include <TDirectory.h>
25 #include <math.h>
26 #include <iostream>
27 #include <algorithm>
28 #include <mdst/dataobjects/Track.h>
29 #include <mdst/dataobjects/TrackFitResult.h>
30 
31 using namespace Belle2;
32 using namespace std;
33 
34 
35 //-----------------------------------------------------------------
36 // Register the Module
37 //-----------------------------------------------------------------
38 REG_MODULE(OverlapResiduals)
39 
40 //-----------------------------------------------------------------
41 // Implementation
42 //-----------------------------------------------------------------
43 
44 OverlapResidualsModule::OverlapResidualsModule() : HistoModule()
45 {
46  //Set module properties
47  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.");
48  //Parameter to take only specific RecoTracks as input
49  addParam("recoTracksStoreArrayName", m_recoTracksStoreArrayName, "StoreArray name of the input and output RecoTracks.",
50  m_recoTracksStoreArrayName);
51  //Parameter to produce TTrees storing global information on VXD overlaps
52  addParam("ExpertLevel", m_ExpertLevel,
53  "If set, enables the production of TTrees containing low level information on PXD and SVD overlaps", false);
54 }
55 
56 void OverlapResidualsModule::initialize()
57 {
58  StoreArray<RecoTrack> recoTracks(m_recoTracksStoreArrayName);
59  recoTracks.isOptional();
60  //Register histograms (calls back defineHisto)
61  REG_HISTOGRAM
62 }
63 
64 void OverlapResidualsModule::defineHisto()
65 {
66  //Create directory to store monitoring histograms
67  TDirectory* oldDir = gDirectory;
68  TDirectory* ResDir = oldDir->mkdir("Monitoring_VXDOverlaps");
69  ResDir->cd();
70  //Define histograms of residuals differences
71  h_U_DeltaRes = new TH1F("h_U_DeltaRes", "Histrogram of residual difference #Delta res_{u} for overlapping hits", 100, -1000, 1000);
72  h_U_DeltaRes->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
73  h_U_DeltaRes->GetYaxis()->SetTitle("counts");
74  h_V_DeltaRes = new TH1F("h_V_DeltaRes", "Histrogram of residual difference #Delta res_{v} for overlapping hits", 100, -1000, 1000);
75  h_V_DeltaRes->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
76  h_V_DeltaRes->GetYaxis()->SetTitle("counts");
77  h_U_DeltaRes_PXD = new TH1F("h_U_DeltaRes_PXD", "Histrogram of residual difference #Delta res_{u} for overlapping PXD hits", 100,
78  -1000, 1000);
79  h_U_DeltaRes_PXD->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
80  h_U_DeltaRes_PXD->GetYaxis()->SetTitle("counts");
81  h_U_DeltaRes_PXD_Lyr1 = new TH1F("h_U_DeltaRes_PXD_Lyr1",
82  "Layer1: histrogram of residual difference #Delta res_{u} for overlapping PXD hits", 100, -1000, 1000);
83  h_U_DeltaRes_PXD_Lyr1->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
84  h_U_DeltaRes_PXD_Lyr1->GetYaxis()->SetTitle("counts");
85  h_U_DeltaRes_PXD_Lyr2 = new TH1F("h_U_DeltaRes_PXD_Lyr2",
86  "Layer 2: hstrogram of residual difference #Delta res_{u} for overlapping PXD hits", 100, -1000, 1000);
87  h_U_DeltaRes_PXD_Lyr2->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
88  h_U_DeltaRes_PXD_Lyr2->GetYaxis()->SetTitle("counts");
89  h_V_DeltaRes_PXD = new TH1F("h_V_DeltaRes_PXD", "Histrogram of residual difference #Delta res_{v} for overlapping PXD hits", 100,
90  -1000, 1000);
91  h_V_DeltaRes_PXD->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
92  h_V_DeltaRes_PXD->GetYaxis()->SetTitle("counts");
93  h_V_DeltaRes_PXD_Lyr1 = new TH1F("h_V_DeltaRes_PXD_Lyr1",
94  "Layer 1: histrogram of residual difference #Delta res_{v} for overlapping PXD hits", 100, -1000, 1000);
95  h_V_DeltaRes_PXD_Lyr1->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
96  h_V_DeltaRes_PXD_Lyr1->GetYaxis()->SetTitle("counts");
97  h_V_DeltaRes_PXD_Lyr2 = new TH1F("h_V_DeltaRes_PXD_Lyr2",
98  "Layer 2: histrogram of residual difference #Delta res_{v} for overlapping PXD hits", 100, -1000, 1000);
99  h_V_DeltaRes_PXD_Lyr2->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
100  h_V_DeltaRes_PXD_Lyr2->GetYaxis()->SetTitle("counts");
101  h_U_DeltaRes_SVD = new TH1F("h_U_DeltaRes_SVD", "Histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100,
102  -1000, 1000);
103  h_U_DeltaRes_SVD->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
104  h_U_DeltaRes_SVD->GetYaxis()->SetTitle("counts");
105  h_U_DeltaRes_SVD_Lyr3 = new TH1F("h_U_DeltaRes_SVD_Lyr3",
106  "Layer 3: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
107  h_U_DeltaRes_SVD_Lyr3->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
108  h_U_DeltaRes_SVD_Lyr3->GetYaxis()->SetTitle("counts");
109  h_U_DeltaRes_SVD_Lyr4 = new TH1F("h_U_DeltaRes_SVD_Lyr4",
110  "Layer 4: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
111  h_U_DeltaRes_SVD_Lyr4->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
112  h_U_DeltaRes_SVD_Lyr4->GetYaxis()->SetTitle("counts");
113  h_U_DeltaRes_SVD_Lyr5 = new TH1F("h_U_DeltaRes_SVD_Lyr5",
114  "Layer 5: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
115  h_U_DeltaRes_SVD_Lyr5->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
116  h_U_DeltaRes_SVD_Lyr5->GetYaxis()->SetTitle("counts");
117  h_U_DeltaRes_SVD_Lyr6 = new TH1F("h_U_DeltaRes_SVD_Lyr6",
118  "Layer 6: histrogram of residual difference #Delta res_{u} for overlapping SVD hits", 100, -1000, 1000);
119  h_U_DeltaRes_SVD_Lyr6->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
120  h_U_DeltaRes_SVD_Lyr6->GetYaxis()->SetTitle("counts");
121  h_V_DeltaRes_SVD = new TH1F("h_V_DeltaRes_SVD", "Histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100,
122  -1000, 1000);
123  h_V_DeltaRes_SVD->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
124  h_V_DeltaRes_SVD->GetYaxis()->SetTitle("counts");
125  h_V_DeltaRes_SVD_Lyr3 = new TH1F("h_V_DeltaRes_SVD_Lyr3",
126  "Layer 3: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
127  h_V_DeltaRes_SVD_Lyr3->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
128  h_V_DeltaRes_SVD_Lyr3->GetYaxis()->SetTitle("counts");
129  h_V_DeltaRes_SVD_Lyr4 = new TH1F("h_V_DeltaRes_SVD_Lyr4",
130  "Layer 4: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
131  h_V_DeltaRes_SVD_Lyr4->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
132  h_V_DeltaRes_SVD_Lyr4->GetYaxis()->SetTitle("counts");
133  h_V_DeltaRes_SVD_Lyr5 = new TH1F("h_V_DeltaRes_SVD_Lyr5",
134  "Layer 5: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
135  h_V_DeltaRes_SVD_Lyr5->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
136  h_V_DeltaRes_SVD_Lyr5->GetYaxis()->SetTitle("counts");
137  h_V_DeltaRes_SVD_Lyr6 = new TH1F("h_V_DeltaRes_SVD_Lyr6",
138  "Layer 6: histrogram of residual difference #Delta res_{v} for overlapping SVD hits", 100, -1000, 1000);
139  h_V_DeltaRes_SVD_Lyr6->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
140  h_V_DeltaRes_SVD_Lyr6->GetYaxis()->SetTitle("counts");
141  h_SVDstrips_Mult = new TH1F("h_SVDstrips_Mult", "SVD strips multipicity for SVD clusters in overlapping sensors", 15, 0.5, 15.5);
142  h_SVDstrips_Mult->GetXaxis()->SetTitle("N. of SVD strips contributing to the cluster");
143  h_SVDstrips_Mult->GetYaxis()->SetTitle("counts");
144  //Define 2D histograms: difference of u-residuals vs phi of VXD overlaps for each layer (1 to 6)
145  h_DeltaResUPhi_Lyr1 = new TH2F("h_DeltaResUPhi_Lyr1", "Layer 1: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
146  h_DeltaResUPhi_Lyr1->GetXaxis()->SetTitle("#phi(rad)");
147  h_DeltaResUPhi_Lyr1->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
148  h_DeltaResUPhi_Lyr2 = new TH2F("h_DeltaResUPhi_Lyr2", "Layer 2: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
149  h_DeltaResUPhi_Lyr2->GetXaxis()->SetTitle("#phi(rad)");
150  h_DeltaResUPhi_Lyr2->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
151  h_DeltaResUPhi_Lyr3 = new TH2F("h_DeltaResUPhi_Lyr3", "Layer 3: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
152  h_DeltaResUPhi_Lyr3->GetXaxis()->SetTitle("#phi(rad)");
153  h_DeltaResUPhi_Lyr3->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
154  h_DeltaResUPhi_Lyr4 = new TH2F("h_DeltaResUPhi_Lyr4", "Layer 4: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
155  h_DeltaResUPhi_Lyr4->GetXaxis()->SetTitle("#phi(rad)");
156  h_DeltaResUPhi_Lyr4->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
157  h_DeltaResUPhi_Lyr5 = new TH2F("h_DeltaResUPhi_Lyr5", "Layer 5: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
158  h_DeltaResUPhi_Lyr5->GetXaxis()->SetTitle("#phi(rad)");
159  h_DeltaResUPhi_Lyr5->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
160  h_DeltaResUPhi_Lyr6 = new TH2F("h_DeltaResUPhi_Lyr6", "Layer 6: #Delta res_{u} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
161  h_DeltaResUPhi_Lyr6->GetXaxis()->SetTitle("#phi(rad)");
162  h_DeltaResUPhi_Lyr6->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
163  //Define 2D histograms: difference of u-residuals vs z of VXD overlaps for each layer (1 to 6)
164  h_DeltaResUz_Lyr1 = new TH2F("h_DeltaResUz_Lyr1", "Layer 1: #Delta res_{u} vs z", 100, -4, 8, 100, -200, 200);
165  h_DeltaResUz_Lyr1->GetXaxis()->SetTitle("z (cm)");
166  h_DeltaResUz_Lyr1->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
167  h_DeltaResUz_Lyr2 = new TH2F("h_DeltaResUz_Lyr2", "Layer 2: #Delta res_{u} vs z", 100, -10, 15, 100, -200, 200);
168  h_DeltaResUz_Lyr2->GetXaxis()->SetTitle("z (cm)");
169  h_DeltaResUz_Lyr2->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
170  h_DeltaResUz_Lyr3 = new TH2F("h_DeltaResUz_Lyr3", "Layer 3: #Delta res_{u} vs z", 250, -15, 20, 100, -1000, 1000);
171  h_DeltaResUz_Lyr3->GetXaxis()->SetTitle("z (cm)");
172  h_DeltaResUz_Lyr3->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
173  h_DeltaResUz_Lyr4 = new TH2F("h_DeltaResUz_Lyr4", "Layer 4: #Delta res_{u} vs z", 250, -20, 25, 100, -1000, 1000);
174  h_DeltaResUz_Lyr4->GetXaxis()->SetTitle("z (cm)");
175  h_DeltaResUz_Lyr4->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
176  h_DeltaResUz_Lyr5 = new TH2F("h_DeltaResUz_Lyr5", "Layer 5: #Delta res_{u} vs z", 250, -25, 35, 100, -1000, 1000);
177  h_DeltaResUz_Lyr5->GetXaxis()->SetTitle("z (cm)");
178  h_DeltaResUz_Lyr5->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
179  h_DeltaResUz_Lyr6 = new TH2F("h_DeltaResUz_Lyr6", "Layer 6: #Delta res_{u} vs z", 250, -30, 45, 100, -1000, 1000);
180  h_DeltaResUz_Lyr6->GetXaxis()->SetTitle("z (cm)");
181  h_DeltaResUz_Lyr6->GetYaxis()->SetTitle("#Delta res_{u} (#mum)");
182  //Define 2D histograms: difference of u-residuals vs phi of VXD overlaps for each layer (1 to 6)
183  h_DeltaResVPhi_Lyr1 = new TH2F("h_DeltaResVPhi_Lyr1", "Layer 1: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
184  h_DeltaResVPhi_Lyr1->GetXaxis()->SetTitle("#phi(rad)");
185  h_DeltaResVPhi_Lyr1->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
186  h_DeltaResVPhi_Lyr2 = new TH2F("h_DeltaResVPhi_Lyr2", "Layer 2: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -200, 200);
187  h_DeltaResVPhi_Lyr2->GetXaxis()->SetTitle("#phi(rad)");
188  h_DeltaResVPhi_Lyr2->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
189  h_DeltaResVPhi_Lyr3 = new TH2F("h_DeltaResVPhi_Lyr3", "Layer 3: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
190  h_DeltaResVPhi_Lyr3->GetXaxis()->SetTitle("#phi(rad)");
191  h_DeltaResVPhi_Lyr3->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
192  h_DeltaResVPhi_Lyr4 = new TH2F("h_DeltaResVPhi_Lyr4", "Layer 4: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
193  h_DeltaResVPhi_Lyr4->GetXaxis()->SetTitle("#phi(rad)");
194  h_DeltaResVPhi_Lyr4->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
195  h_DeltaResVPhi_Lyr5 = new TH2F("h_DeltaResVPhi_Lyr5", "Layer 5: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
196  h_DeltaResVPhi_Lyr5->GetXaxis()->SetTitle("#phi(rad)");
197  h_DeltaResVPhi_Lyr5->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
198  h_DeltaResVPhi_Lyr6 = new TH2F("h_DeltaResVPhi_Lyr6", "Layer 6: #Delta res_{v} vs #phi", 200, -3.4, 3.4, 100, -1000, 1000);
199  h_DeltaResVPhi_Lyr6->GetXaxis()->SetTitle("#phi(rad)");
200  h_DeltaResVPhi_Lyr6->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
201  //Define 2D histograms: difference of v-residuals vs z of VXD overlaps for each layer (1 to 6)
202  h_DeltaResVz_Lyr1 = new TH2F("h_DeltaResVz_Lyr1", "Layer 1: #Delta res_{v} vs z", 100, -4, 8, 100, -200, 200);
203  h_DeltaResVz_Lyr1->GetXaxis()->SetTitle("z (cm)");
204  h_DeltaResVz_Lyr1->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
205  h_DeltaResVz_Lyr2 = new TH2F("h_DeltaResVz_Lyr2", "Layer 2: #Delta res_{v} vs z", 100, -10, 15, 100, -200, 200);
206  h_DeltaResVz_Lyr2->GetXaxis()->SetTitle("z (cm)");
207  h_DeltaResVz_Lyr2->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
208  h_DeltaResVz_Lyr3 = new TH2F("h_DeltaResVz_Lyr3", "Layer 3: #Delta res_{v} vs z", 250, -15, 20, 100, -1000, 1000);
209  h_DeltaResVz_Lyr3->GetXaxis()->SetTitle("z (cm)");
210  h_DeltaResVz_Lyr3->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
211  h_DeltaResVz_Lyr4 = new TH2F("h_DeltaResVz_Lyr4", "Layer 4: #Delta res_{v} vs z", 250, -20, 25, 100, -1000, 1000);
212  h_DeltaResVz_Lyr4->GetXaxis()->SetTitle("z (cm)");
213  h_DeltaResVz_Lyr4->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
214  h_DeltaResVz_Lyr5 = new TH2F("h_DeltaResVz_Lyr5", "Layer 5: #Delta res_{v} vs z", 250, -25, 35, 100, -1000, 1000);
215  h_DeltaResVz_Lyr5->GetXaxis()->SetTitle("z (cm)");
216  h_DeltaResVz_Lyr5->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
217  h_DeltaResVz_Lyr6 = new TH2F("h_DeltaResVz_Lyr6", "Layer 6: #Delta res_{v} vs z", 250, -30, 45, 100, -1000, 1000);
218  h_DeltaResVz_Lyr6->GetXaxis()->SetTitle("z (cm)");
219  h_DeltaResVz_Lyr6->GetYaxis()->SetTitle("#Delta res_{v} (#mum)");
220  //Restricting to SVD clusters sizes
221  for (int i = 1; i < 5; i++) {
222  //The name is the product of cluster sizes for 2 consecutive hits (maximum size considered is 2)
223  TString h_name_U = "h_U_Cl1Cl2_" + std::to_string(i);
224  TString h_name_V = "h_V_Cl1Cl2_" + std::to_string(i);
225  TString title_U = "#Delta res_{u}: SVDClusterSize_1 x SVDClusterSize_2 = " + std::to_string(i);
226  TString title_V = "#Delta res_{v}: SVDClusterSize_1 x SVDClusterSize_2 = " + std::to_string(i);
227  h_U_Cl1Cl2_DeltaRes[i] = new TH1F(h_name_U, title_U, 100, -1000, 1000);
228  h_U_Cl1Cl2_DeltaRes[i]->GetXaxis()->SetTitle("#Delta res_{u} (#mum)");
229  h_U_Cl1Cl2_DeltaRes[i]->GetYaxis()->SetTitle("counts");
230  h_V_Cl1Cl2_DeltaRes[i] = new TH1F(h_name_V, title_V, 100, -1000, 1000);
231  h_V_Cl1Cl2_DeltaRes[i]->GetXaxis()->SetTitle("#Delta res_{v} (#mum)");
232  h_V_Cl1Cl2_DeltaRes[i]->GetYaxis()->SetTitle("counts");
233  }
234 
235  //Special case of ExpertLevel option enabled
236  if (m_ExpertLevel) {
237  //Create directory to store PXD and SVD hitmaps for overlapping hits
238  TDirectory* HMDir = oldDir->mkdir("HitMaps_VXDOverlaps");
239  HMDir->cd();
240  //Define 2D sensor hit-maps for reconstructed hits
241  for (int i = 1; i <= 5; i++) {
242  for (int j = 1; j <= 16; j++) {
243  TString h_name = "h_6" + std::to_string(j) + std::to_string(i);
244  TString title = "Layer:Ladder:Sensor = 6:" + std::to_string(j) + ":" + std::to_string(i);
245  h_Lyr6[j][i] = new TH2F(h_name, title, 100, -2.88, 2.88, 100, -6.14, 6.14);
246  h_Lyr6[j][i]->GetXaxis()->SetTitle("u (cm)");
247  h_Lyr6[j][i]->GetYaxis()->SetTitle("v (cm)");
248  }
249  }
250  for (int i = 1; i <= 4; i++) {
251  for (int j = 1; j <= 12; j++) {
252  TString h_name = "h_5" + std::to_string(j) + std::to_string(i);
253  TString title = "Layer:Ladder:Sensor = 5:" + std::to_string(j) + ":" + std::to_string(i);
254  h_Lyr5[j][i] = new TH2F(h_name, title, 100, -2.88, 2.88, 100, -6.14, 6.14);
255  h_Lyr5[j][i]->GetXaxis()->SetTitle("u (cm)");
256  h_Lyr5[j][i]->GetYaxis()->SetTitle("v (cm)");
257  }
258  }
259  for (int i = 1; i <= 3; i++) {
260  for (int j = 1; j <= 10; j++) {
261  TString h_name = "h_4" + std::to_string(j) + std::to_string(i);
262  TString title = "Layer:Ladder:Sensor = 4:" + std::to_string(j) + ":" + std::to_string(i);
263  h_Lyr4[j][i] = new TH2F(h_name, title, 100, -2.88, 2.88, 100, -6.14, 6.14);
264  h_Lyr4[j][i]->GetXaxis()->SetTitle("u (cm)");
265  h_Lyr4[j][i]->GetYaxis()->SetTitle("v (cm)");
266  }
267  }
268  for (int i = 1; i <= 2; i++) {
269  for (int j = 1; j <= 7; j++) {
270  TString h_name = "h_3" + std::to_string(j) + std::to_string(i);
271  TString title = "Layer:Ladder:Sensor = 3:" + std::to_string(j) + ":" + std::to_string(i);
272  h_Lyr3[j][i] = new TH2F(h_name, title, 100, -1.92, 1.92, 100, -6.14, 6.14);
273  h_Lyr3[j][i]->GetXaxis()->SetTitle("u (cm)");
274  h_Lyr3[j][i]->GetYaxis()->SetTitle("v (cm)");
275  }
276  }
277  for (int i = 1; i <= 2; i++) {
278  for (int j = 1; j <= 12; j++) {
279  TString h_name = "h_2" + std::to_string(j) + std::to_string(i);
280  TString title = "Layer:Ladder:Sensor = 2:" + std::to_string(j) + ":" + std::to_string(i);
281  h_Lyr2[j][i] = new TH2F(h_name, title, 100, -0.625, 0.625, 100, -3.072, 3.072);
282  h_Lyr2[j][i]->GetXaxis()->SetTitle("u (cm)");
283  h_Lyr2[j][i]->GetYaxis()->SetTitle("v (cm)");
284  }
285  }
286  for (int i = 1; i <= 2; i++) {
287  for (int j = 1; j <= 8; j++) {
288  TString h_name = "h_1" + std::to_string(j) + std::to_string(i);
289  TString title = "Layer:Ladder:Sensor = 1:" + std::to_string(j) + ":" + std::to_string(i);
290  h_Lyr1[j][i] = new TH2F(h_name, title, 100, -0.625, 0.625, 100, -2.24, 2.24);
291  h_Lyr1[j][i]->GetXaxis()->SetTitle("u (cm)");
292  h_Lyr1[j][i]->GetYaxis()->SetTitle("v (cm)");
293  }
294  }
295  //Create directory to store PXD and SVD trees
296  TDirectory* TreeDir = oldDir->mkdir("Trees_VXDOverlaps");
297  TreeDir->cd();
298  //Tree for PXD
299  t_PXD = new TTree("t_PXD", "Tree for PXD overlaps");
300  t_PXD->Branch("deltaResU_PXD", &deltaResU_PXD, "deltaResU_PXD/F");
301  t_PXD->Branch("intResU_PXD", &intResU_PXD, "intResU_PXD/F");
302  t_PXD->Branch("intResV_PXD", &intResV_PXD, "intResV_PXD/F");
303  t_PXD->Branch("intU_PXD", &intU_PXD, "intU_PXD/F");
304  t_PXD->Branch("intV_PXD", &intV_PXD, "intV_PXD/F");
305  t_PXD->Branch("intPhi_PXD", &intPhi_PXD, "intPhi_PXD/F");
306  t_PXD->Branch("intZ_PXD", &intZ_PXD, "intZ_PXD/F");
307  t_PXD->Branch("intLayer_PXD", &intLayer_PXD, "intLayer_PXD/i");
308  t_PXD->Branch("intLadder_PXD", &intLadder_PXD, "intLadder_PXD/i");
309  t_PXD->Branch("intSensor_PXD", &intSensor_PXD, "intSensor_PXD/i");
310  t_PXD->Branch("extResU_PXD", &extResU_PXD, "extResU_PXD/F");
311  t_PXD->Branch("extResV_PXD", &extResV_PXD, "extResV_PXD/F");
312  t_PXD->Branch("extU_PXD", &extU_PXD, "extU_PXD/F");
313  t_PXD->Branch("extV_PXD", &extV_PXD, "extV_PXD/F");
314  t_PXD->Branch("extPhi_PXD", &extPhi_PXD, "extPhi_PXD/F");
315  t_PXD->Branch("extZ_PXD", &extZ_PXD, "extZ_PXD/F");
316  t_PXD->Branch("extLayer_PXD", &extLayer_PXD, "extLayer_PXD/i");
317  t_PXD->Branch("extLadder_PXD", &extLadder_PXD, "extLadder_PXD/i");
318  t_PXD->Branch("extSensor_PXD", &extSensor_PXD, "extSensor_PXD/i");
319  //Tree for SVD u overlapping clusters
320  t_SVD_U = new TTree("t_SVD_U", "Tree for SVD u-overlaps");
321  t_SVD_U->Branch("svdDeltaRes", &svdDeltaRes_U, "svdDeltaRes/F");
322  t_SVD_U->Branch("svdTrkPXDHits", &svdTrkPXDHits, "svdTrkPXDHits/i");
323  t_SVD_U->Branch("svdTrkSVDHits", &svdTrkSVDHits, "svdTrkSVDHits/i");
324  t_SVD_U->Branch("svdTrkCDCHits", &svdTrkCDCHits, "svdTrkCDCHits/i");
325  t_SVD_U->Branch("svdTrkd0", &svdTrkd0, "svdTrkd0/F");
326  t_SVD_U->Branch("svdTrkz0", &svdTrkz0, "svdTrkz0/F");
327  t_SVD_U->Branch("svdTrkpT", &svdTrkpT, "svdTrkpT/F");
328  t_SVD_U->Branch("svdTrkpCM", &svdTrkpCM, "svdTrkpCM/F");
329  // Internal ladder variables
330  t_SVD_U->Branch("svdClSNR_int", &svdClSNR_U_int, "svdClSNR_int/F");
331  t_SVD_U->Branch("svdClCharge_int", &svdClCharge_U_int, "svdClCharge_int/F");
332  t_SVD_U->Branch("svdStripCharge_int", &svdStripCharge_U_int);
333  t_SVD_U->Branch("svdStrip6Samples_int", &svdStrip6Samples_U_int);
334  t_SVD_U->Branch("svdClTime_int", &svdClTime_U_int, "svdClTime_int/F");
335  t_SVD_U->Branch("svdStripTime_int", &svdStripTime_U_int);
336  t_SVD_U->Branch("svdStripPosition_int", &svdStripPosition_U_int);
337  t_SVD_U->Branch("svdRes_int", &svdRes_U_int, "svdRes_int/F");
338  t_SVD_U->Branch("svdClIntStrPos_int", &svdClIntStrPos_U_int, "svdClIntStrPos_int/F");
339  t_SVD_U->Branch("svdClPos_int", &svdClPos_U_int, "svdClPos_int/F");
340  t_SVD_U->Branch("svdClPosErr_int", &svdClPosErr_U_int, "svdClPosErr_int/F");
341  t_SVD_U->Branch("svdTruePos_int", &svdTruePos_U_int, "svdTruePos_int/F");
342  t_SVD_U->Branch("svdClPhi_int", &svdClPhi_U_int, "svdClPhi_int/F");
343  t_SVD_U->Branch("svdClZ_int", &svdClZ_U_int, "svdClZ_int/F");
344  t_SVD_U->Branch("svdTrkTraversedLength_int", &svdTrkTraversedLength_U_int, "svdTrkTraversedLength_int/F");
345  t_SVD_U->Branch("svdTrkPos_int", &svdTrkPos_U_int, "svdTrkPos_int/F");
346  t_SVD_U->Branch("svdTrkPosOS_int", &svdTrkPosOS_U_int, "svdTrkPosOS_int/F");
347  t_SVD_U->Branch("svdTrkPosErr_int", &svdTrkPosErr_U_int, "svdTrkPosErr_int/F");
348  t_SVD_U->Branch("svdTrkPosErrOS_int", &svdTrkPosErrOS_U_int, "svdTrkPosErrOS_int/F");
349  t_SVD_U->Branch("svdTrkQoP_int", &svdTrkQoP_U_int, "svdTrkQoP_int/F");
350  t_SVD_U->Branch("svdTrkPrime_int", &svdTrkPrime_U_int, "svdTrkPrime_int/F");
351  t_SVD_U->Branch("svdTrkPrimeOS_int", &svdTrkPrimeOS_U_int, "svdTrkPrimeOS_int/F");
352  t_SVD_U->Branch("svdTrkPosUnbiased_int", &svdTrkPosUnbiased_U_int, "svdTrkPosUnbiased_int/F");
353  t_SVD_U->Branch("svdTrkPosErrUnbiased_int", &svdTrkPosErrUnbiased_U_int, "svdTrkPosErrUnbiased_int/F");
354  t_SVD_U->Branch("svdTrkQoPUnbiased_int", &svdTrkQoPUnbiased_U_int, "svdTrkQoPUnbiased_int/F");
355  t_SVD_U->Branch("svdTrkPrimeUnbiased_int", &svdTrkPrimeUnbiased_U_int, "svdTrkPrimeUnbiased_int/F");
356  t_SVD_U->Branch("svdLayer_int", &svdLayer_U_int, "svdLayer_int/i");
357  t_SVD_U->Branch("svdLadder_int", &svdLadder_U_int, "svdLadder_int/i");
358  t_SVD_U->Branch("svdSensor_int", &svdSensor_U_int, "svdSensor_int/i");
359  t_SVD_U->Branch("svdSize_int", &svdSize_U_int, "svdSize_int/i");
360  // External ladder variables
361  t_SVD_U->Branch("svdClSNR_ext", &svdClSNR_U_ext, "svdClSNR_ext/F");
362  t_SVD_U->Branch("svdClCharge_ext", &svdClCharge_U_ext, "svdClCharge_ext/F");
363  t_SVD_U->Branch("svdStripCharge_ext", &svdStripCharge_U_ext);
364  t_SVD_U->Branch("svdStrip6Samples_ext", &svdStrip6Samples_U_ext);
365  t_SVD_U->Branch("svdClTime_ext", &svdClTime_U_ext, "svdClTime_ext/F");
366  t_SVD_U->Branch("svdStripTime_ext", &svdStripTime_U_ext);
367  t_SVD_U->Branch("svdStripPosition_ext", &svdStripPosition_U_ext);
368  t_SVD_U->Branch("svdRes_ext", &svdRes_U_ext, "svdRes_ext/F");
369  t_SVD_U->Branch("svdClIntStrPos_ext", &svdClIntStrPos_U_ext, "svdClIntStrPos_ext/F");
370  t_SVD_U->Branch("svdClPos_ext", &svdClPos_U_ext, "svdClPos_ext/F");
371  t_SVD_U->Branch("svdClPosErr_ext", &svdClPosErr_U_ext, "svdClPosErr_ext/F");
372  t_SVD_U->Branch("svdTruePos_ext", &svdTruePos_U_ext, "svdTruePos_ext/F");
373  t_SVD_U->Branch("svdClPhi_ext", &svdClPhi_U_ext, "svdClPhi_ext/F");
374  t_SVD_U->Branch("svdClZ_ext", &svdClZ_U_ext, "svdClZ_ext/F");
375  t_SVD_U->Branch("svdTrkTraversedLength_ext", &svdTrkTraversedLength_U_ext, "svdTrkTraversedLength_ext/F");
376  t_SVD_U->Branch("svdTrkPos_ext", &svdTrkPos_U_ext, "svdTrkPos_ext/F");
377  t_SVD_U->Branch("svdTrkPosOS_ext", &svdTrkPosOS_U_ext, "svdTrkPosOS_ext/F");
378  t_SVD_U->Branch("svdTrkPosErr_ext", &svdTrkPosErr_U_ext, "svdTrkPosErr_ext/F");
379  t_SVD_U->Branch("svdTrkPosErrOS_ext", &svdTrkPosErrOS_U_ext, "svdTrkPosErrOS_ext/F");
380  t_SVD_U->Branch("svdTrkQoP_ext", &svdTrkQoP_U_ext, "svdTrkQoP_ext/F");
381  t_SVD_U->Branch("svdTrkPrime_ext", &svdTrkPrime_U_ext, "svdTrkPrime_ext/F");
382  t_SVD_U->Branch("svdTrkPrimeOS_ext", &svdTrkPrimeOS_U_ext, "svdTrkPrimeOS_ext/F");
383  t_SVD_U->Branch("svdTrkPosUnbiased_ext", &svdTrkPosUnbiased_U_ext, "svdTrkPosUnbiased_ext/F");
384  t_SVD_U->Branch("svdTrkPosErrUnbiased_ext", &svdTrkPosErrUnbiased_U_ext, "svdTrkPosErrUnbiased_ext/F");
385  t_SVD_U->Branch("svdTrkQoPUnbiased_ext", &svdTrkQoPUnbiased_U_ext, "svdTrkQoPUnbiased_ext/F");
386  t_SVD_U->Branch("svdTrkPrimeUnbiased_ext", &svdTrkPrimeUnbiased_U_ext, "svdTrkPrimeUnbiased_ext/F");
387  t_SVD_U->Branch("svdLayer_ext", &svdLayer_U_ext, "svdLayer_ext/i");
388  t_SVD_U->Branch("svdLadder_ext", &svdLadder_U_ext, "svdLadder_ext/i");
389  t_SVD_U->Branch("svdSensor_ext", &svdSensor_U_ext, "svdSensor_ext/i");
390  t_SVD_U->Branch("svdSize_ext", &svdSize_U_ext, "svdSize_ext/i");
391  //Tree for SVD v overlapping clusters
392  t_SVD_V = new TTree("t_SVD_V", "Tree for SVD v-overlaps");
393  t_SVD_V->Branch("svdDeltaRes", &svdDeltaRes_V, "svdDeltaRes/F");
394  t_SVD_V->Branch("svdTrkPXDHits", &svdTrkPXDHits, "svdTrkPXDHits/i");
395  t_SVD_V->Branch("svdTrkSVDHits", &svdTrkSVDHits, "svdTrkSVDHits/i");
396  t_SVD_V->Branch("svdTrkCDCHits", &svdTrkCDCHits, "svdTrkCDCHits/i");
397  t_SVD_V->Branch("svdTrkd0", &svdTrkd0, "svdTrkd0/F");
398  t_SVD_V->Branch("svdTrkz0", &svdTrkz0, "svdTrkz0/F");
399  t_SVD_V->Branch("svdTrkpT", &svdTrkpT, "svdTrkpT/F");
400  t_SVD_V->Branch("svdTrkpCM", &svdTrkpCM, "svdTrkpCM/F");
401  // Internal ladder variables
402  t_SVD_V->Branch("svdClSNR_int", &svdClSNR_V_int, "svdClSNR_int/F");
403  t_SVD_V->Branch("svdClCharge_int", &svdClCharge_V_int, "svdClCharge_int/F");
404  t_SVD_V->Branch("svdStripCharge_int", &svdStripCharge_V_int);
405  t_SVD_V->Branch("svdStrip6Samples_int", &svdStrip6Samples_V_int);
406  t_SVD_V->Branch("svdClTime_int", &svdClTime_V_int, "svdClTime_int/F");
407  t_SVD_V->Branch("svdStripTime_int", &svdStripTime_V_int);
408  t_SVD_V->Branch("svdStripPosition_int", &svdStripPosition_V_int);
409  t_SVD_V->Branch("svdRes_int", &svdRes_V_int, "svdRes_int/F");
410  t_SVD_V->Branch("svdClIntStrPos_int", &svdClIntStrPos_V_int, "svdClIntStrPos_int/F");
411  t_SVD_V->Branch("svdClPos_int", &svdClPos_V_int, "svdClPos_int/F");
412  t_SVD_V->Branch("svdClPosErr_int", &svdClPosErr_V_int, "svdClPosErr_int/F");
413  t_SVD_V->Branch("svdTruePos_int", &svdTruePos_V_int, "svdTruePos_int/F");
414  t_SVD_V->Branch("svdClPhi_int", &svdClPhi_V_int, "svdClPhi_int/F");
415  t_SVD_V->Branch("svdClZ_int", &svdClZ_V_int, "svdClZ_int/F");
416  t_SVD_V->Branch("svdTrkTraversedLength_int", &svdTrkTraversedLength_V_int, "svdTrkTraversedLength_int/F");
417  t_SVD_V->Branch("svdTrkPos_int", &svdTrkPos_V_int, "svdTrkPos_int/F");
418  t_SVD_V->Branch("svdTrkPosOS_int", &svdTrkPosOS_V_int, "svdTrkPosOS_int/F");
419  t_SVD_V->Branch("svdTrkPosErr_int", &svdTrkPosErr_V_int, "svdTrkPosErr_int/F");
420  t_SVD_V->Branch("svdTrkPosErrOS_int", &svdTrkPosErrOS_V_int, "svdTrkPosErrOS_int/F");
421  t_SVD_V->Branch("svdTrkQoP_int", &svdTrkQoP_V_int, "svdTrkQoP_int/F");
422  t_SVD_V->Branch("svdTrkPrime_int", &svdTrkPrime_V_int, "svdTrkPrime_int/F");
423  t_SVD_V->Branch("svdTrkPrimeOS_int", &svdTrkPrimeOS_V_int, "svdTrkPrimeOS_int/F");
424  t_SVD_V->Branch("svdTrkPosUnbiased_int", &svdTrkPosUnbiased_V_int, "svdTrkPosUnbiased_int/F");
425  t_SVD_V->Branch("svdTrkPosErrUnbiased_int", &svdTrkPosErrUnbiased_V_int, "svdTrkPosErrUnbiased_int/F");
426  t_SVD_V->Branch("svdTrkQoPUnbiased_int", &svdTrkQoPUnbiased_V_int, "svdTrkQoPUnbiased_int/F");
427  t_SVD_V->Branch("svdTrkPrimeUnbiased_int", &svdTrkPrimeUnbiased_V_int, "svdTrkPrimeUnbiased_int/F");
428  t_SVD_V->Branch("svdLayer_int", &svdLayer_V_int, "svdLayer_int/i");
429  t_SVD_V->Branch("svdLadder_int", &svdLadder_V_int, "svdLadder_int/i");
430  t_SVD_V->Branch("svdSensor_int", &svdSensor_V_int, "svdSensor_int/i");
431  t_SVD_V->Branch("svdSize_int", &svdSize_V_int, "svdSize_int/i");
432  // External ladder variables
433  t_SVD_V->Branch("svdClSNR_ext", &svdClSNR_V_ext, "svdClSNR_ext/F");
434  t_SVD_V->Branch("svdClCharge_ext", &svdClCharge_V_ext, "svdClCharge_ext/F");
435  t_SVD_V->Branch("svdStripCharge_ext", &svdStripCharge_V_ext);
436  t_SVD_V->Branch("svdStrip6Samples_ext", &svdStrip6Samples_V_ext);
437  t_SVD_V->Branch("svdClTime_ext", &svdClTime_V_ext, "svdClTime_ext/F");
438  t_SVD_V->Branch("svdStripTime_ext", &svdStripTime_V_ext);
439  t_SVD_V->Branch("svdStripPosition_ext", &svdStripPosition_V_ext);
440  t_SVD_V->Branch("svdRes_ext", &svdRes_V_ext, "svdRes_ext/F");
441  t_SVD_V->Branch("svdClIntStrPos_ext", &svdClIntStrPos_V_ext, "svdClIntStrPos_ext/F");
442  t_SVD_V->Branch("svdClPos_ext", &svdClPos_V_ext, "svdClPos_ext/F");
443  t_SVD_V->Branch("svdClPosErr_ext", &svdClPosErr_V_ext, "svdClPosErr_ext/F");
444  t_SVD_V->Branch("svdTruePos_ext", &svdTruePos_V_ext, "svdTruePos_ext/F");
445  t_SVD_V->Branch("svdClPhi_ext", &svdClPhi_V_ext, "svdClPhi_ext/F");
446  t_SVD_V->Branch("svdClZ_ext", &svdClZ_V_ext, "svdClZ_ext/F");
447  t_SVD_V->Branch("svdTrkTraversedLength_ext", &svdTrkTraversedLength_V_ext, "svdTrkTraversedLength_ext/F");
448  t_SVD_V->Branch("svdTrkPos_ext", &svdTrkPos_V_ext, "svdTrkPos_ext/F");
449  t_SVD_V->Branch("svdTrkPosOS_ext", &svdTrkPosOS_V_ext, "svdTrkPosOS_ext/F");
450  t_SVD_V->Branch("svdTrkPosErr_ext", &svdTrkPosErr_V_ext, "svdTrkPosErr_ext/F");
451  t_SVD_V->Branch("svdTrkPosErrOS_ext", &svdTrkPosErrOS_V_ext, "svdTrkPosErrOS_ext/F");
452  t_SVD_V->Branch("svdTrkQoP_ext", &svdTrkQoP_V_ext, "svdTrkQoP_ext/F");
453  t_SVD_V->Branch("svdTrkPrime_ext", &svdTrkPrime_V_ext, "svdTrkPrime_ext/F");
454  t_SVD_V->Branch("svdTrkPrimeOS_ext", &svdTrkPrimeOS_V_ext, "svdTrkPrimeOS_ext/F");
455  t_SVD_V->Branch("svdTrkPosUnbiased_ext", &svdTrkPosUnbiased_V_ext, "svdTrkPosUnbiased_ext/F");
456  t_SVD_V->Branch("svdTrkPosErrUnbiased_ext", &svdTrkPosErrUnbiased_V_ext, "svdTrkPosErrUnbiased_ext/F");
457  t_SVD_V->Branch("svdTrkQoPUnbiased_ext", &svdTrkQoPUnbiased_V_ext, "svdTrkQoPUnbiased_ext/F");
458  t_SVD_V->Branch("svdTrkPrimeUnbiased_ext", &svdTrkPrimeUnbiased_V_ext, "svdTrkPrimeUnbiased_ext/F");
459  t_SVD_V->Branch("svdLayer_ext", &svdLayer_V_ext, "svdLayer_ext/i");
460  t_SVD_V->Branch("svdLadder_ext", &svdLadder_V_ext, "svdLadder_ext/i");
461  t_SVD_V->Branch("svdSensor_ext", &svdSensor_V_ext, "svdSensor_ext/i");
462  t_SVD_V->Branch("svdSize_ext", &svdSize_V_ext, "svdSize_ext/i");
463  }
464  //Go back to the default output directory
465  oldDir->cd();
466 }
467 
468 void OverlapResidualsModule::event()
469 {
470  bool isMC = Environment::Instance().isMC();
471 
472  static VXD::GeoCache& geo = VXD::GeoCache::getInstance();
473  StoreArray<RecoTrack> recoTracks(m_recoTracksStoreArrayName);
474  for (const auto& trk : recoTracks) {
475  if (! trk.wasFitSuccessful()) {
476  continue;
477  }
478  const vector<PXDCluster* > pxdClusters = trk.getPXDHitList();
479  const vector<SVDCluster* > svdClusters = trk.getSVDHitList();
480  B2DEBUG(40, "FITTED TRACK: NUMBER OF PXD HITS = " << pxdClusters.size() << " NUMBER OF SVD HITS = " << svdClusters.size());
481  //LOOKING FOR 2 CONSECUTIVE PXD HITS IN OVERLAPPING MODULES OF A SAME LAYER
482  for (unsigned int i = 0; i < pxdClusters.size(); i++) {
483  const PXDCluster* pxd_1 = pxdClusters[i];
484  const RecoHitInformation* infoPXD_1 = trk.getRecoHitInformation(pxd_1);
485  if (!infoPXD_1) {
486  continue;
487  }
488  const auto* hitTrackPoint_1 = trk.getCreatedTrackPoint(infoPXD_1);
489  const auto* fittedResult_1 = hitTrackPoint_1->getFitterInfo();
490  if (!fittedResult_1) {
491  continue;
492  }
493  const VxdID pxd_id_1 = pxd_1->getSensorID();
494  const unsigned short pxd_Layer_1 = pxd_id_1.getLayerNumber();
495  const unsigned short pxd_Ladder_1 = pxd_id_1.getLadderNumber();
496  const unsigned short pxd_Sensor_1 = pxd_id_1.getSensorNumber();
497  for (unsigned int l = i + 1; l < pxdClusters.size(); l++) {
498  const PXDCluster* pxd_2 = pxdClusters[l];
499  const RecoHitInformation* infoPXD_2 = trk.getRecoHitInformation(pxd_2);
500  if (!infoPXD_2) {
501  continue;
502  }
503  const auto* hitTrackPoint_2 = trk.getCreatedTrackPoint(infoPXD_2);
504  const auto* fittedResult_2 = hitTrackPoint_2->getFitterInfo();
505  if (!fittedResult_2) {
506  continue;
507  }
508  const VxdID pxd_id_2 = pxd_2->getSensorID();
509  const unsigned short pxd_Layer_2 = pxd_id_2.getLayerNumber();
510  const unsigned short pxd_Ladder_2 = pxd_id_2.getLadderNumber();
511  const unsigned short pxd_Sensor_2 = pxd_id_2.getSensorNumber();
512  if (pxd_Layer_1 == pxd_Layer_2 && ((pxd_Ladder_2 - pxd_Ladder_1) == -1.0 || (pxd_Ladder_2 - pxd_Ladder_1) == 7.0
513  || (pxd_Ladder_2 - pxd_Ladder_1) == 11.0)) {
514  B2DEBUG(40, " ============= 2 hits in a PXD overlap ============= ");
515  const TVectorD resUnBias_PXD_1 = fittedResult_1->getResidual(0, false).getState();
516  const TVectorD resUnBias_PXD_2 = fittedResult_2->getResidual(0, false).getState();
517  const float res_U_1 = resUnBias_PXD_1.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
518  const float res_V_1 = resUnBias_PXD_1.GetMatrixArray()[1] * Unit::convertValueToUnit(1.0, "um");
519  const double res_U_2 = resUnBias_PXD_2.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
520  const double res_V_2 = resUnBias_PXD_2.GetMatrixArray()[1] * Unit::convertValueToUnit(1.0, "um");
521  const float over_U_PXD = res_U_2 - res_U_1;
522  const float over_V_PXD = res_V_2 - res_V_1;
523  const TVector3 pxdLocal_1(pxd_1->getU(), pxd_1->getV(), 0.);
524  const TVector3 pxdLocal_2(pxd_2->getU(), pxd_2->getV(), 0.);
525  const VXD::SensorInfoBase& pxdSensor_1 = geo.get(pxd_id_1);
526  const VXD::SensorInfoBase& pxdSensor_2 = geo.get(pxd_id_2);
527  const TVector3& pxdGlobal_1 = pxdSensor_1.pointToGlobal(pxdLocal_1);
528  const TVector3& pxdGlobal_2 = pxdSensor_2.pointToGlobal(pxdLocal_2);
529  double pxdPhi_1 = atan2(pxdGlobal_1(1), pxdGlobal_1(0));
530  double pxdPhi_2 = atan2(pxdGlobal_2(1), pxdGlobal_2(0));
531  double pxdZ_1 = pxdGlobal_1(2);
532  double pxdZ_2 = pxdGlobal_2(2);
533  B2DEBUG(40, "PXD: difference of residuals " << over_U_PXD << " " << over_V_PXD);
534  //Fill PXD tree for overlaps if required by the user
535  if (m_ExpertLevel) {
536  deltaResU_PXD = over_U_PXD;
537  intResU_PXD = res_U_1;
538  intResV_PXD = res_V_1;
539  intU_PXD = pxd_1->getU();
540  intV_PXD = pxd_1->getV();
541  intPhi_PXD = pxdPhi_1;
542  intZ_PXD = pxdZ_1;
543  intLayer_PXD = pxd_Layer_1;
544  intLadder_PXD = pxd_Ladder_1;
545  intSensor_PXD = pxd_Sensor_1;
546  extResU_PXD = res_U_2;
547  extResV_PXD = res_V_2;
548  extU_PXD = pxd_2->getU();
549  extV_PXD = pxd_2->getV();
550  extPhi_PXD = pxdPhi_2;
551  extZ_PXD = pxdZ_2;
552  extLayer_PXD = pxd_Layer_2;
553  extLadder_PXD = pxd_Ladder_2;
554  extSensor_PXD = pxd_Sensor_2;
555  t_PXD->Fill();
556  }
557  //Fill histograms of residuals differences with PXD clusters
558  h_U_DeltaRes->Fill(over_U_PXD);
559  h_V_DeltaRes->Fill(over_V_PXD);
560  h_U_DeltaRes_PXD->Fill(over_U_PXD);
561  h_V_DeltaRes_PXD->Fill(over_V_PXD);
562  //Fill sensor hit-maps and 2D histograms with PXD clusters
563  if (pxd_Layer_1 == 1 && pxd_Layer_2 == 1) {
564  h_Lyr1[pxd_Ladder_1][pxd_Sensor_1]->Fill(pxd_1->getU(), pxd_1->getV());
565  h_Lyr1[pxd_Ladder_2][pxd_Sensor_2]->Fill(pxd_2->getU(), pxd_2->getV());
566  h_U_DeltaRes_PXD_Lyr1->Fill(over_U_PXD);
567  h_V_DeltaRes_PXD_Lyr1->Fill(over_V_PXD);
568  h_DeltaResUPhi_Lyr1->Fill(pxdPhi_1, over_U_PXD);
569  h_DeltaResUPhi_Lyr1->Fill(pxdPhi_2, over_U_PXD);
570  h_DeltaResVPhi_Lyr1->Fill(pxdPhi_1, over_V_PXD);
571  h_DeltaResVPhi_Lyr1->Fill(pxdPhi_2, over_V_PXD);
572  h_DeltaResVz_Lyr1->Fill(pxdZ_1, over_V_PXD);
573  h_DeltaResVz_Lyr1->Fill(pxdZ_2, over_V_PXD);
574  h_DeltaResUz_Lyr1->Fill(pxdZ_1, over_U_PXD);
575  h_DeltaResUz_Lyr1->Fill(pxdZ_2, over_U_PXD);
576  }
577  if (pxd_Layer_1 == 2 && pxd_Layer_2 == 2) {
578  h_Lyr2[pxd_Ladder_1][pxd_Sensor_1]->Fill(pxd_1->getU(), pxd_1->getV());
579  h_Lyr2[pxd_Ladder_2][pxd_Sensor_2]->Fill(pxd_2->getU(), pxd_2->getV());
580  h_U_DeltaRes_PXD_Lyr2->Fill(over_U_PXD);
581  h_V_DeltaRes_PXD_Lyr2->Fill(over_V_PXD);
582  h_DeltaResUPhi_Lyr2->Fill(pxdPhi_1, over_U_PXD);
583  h_DeltaResUPhi_Lyr2->Fill(pxdPhi_2, over_U_PXD);
584  h_DeltaResVPhi_Lyr2->Fill(pxdPhi_1, over_V_PXD);
585  h_DeltaResVPhi_Lyr2->Fill(pxdPhi_2, over_V_PXD);
586  h_DeltaResVz_Lyr2->Fill(pxdZ_1, over_V_PXD);
587  h_DeltaResVz_Lyr2->Fill(pxdZ_2, over_V_PXD);
588  h_DeltaResUz_Lyr2->Fill(pxdZ_1, over_U_PXD);
589  h_DeltaResUz_Lyr2->Fill(pxdZ_2, over_U_PXD);
590  }
591  }
592  }
593 
594  }
595 
596  //LOOKING FOR 2 CONSECUTIVE SVD HITS IN OVERLAPPING MODULES OF A SAME LAYER
597  RelationVector<Track> theTK = DataStore::getRelationsWithObj<Track>(&trk);
598 
599  const TrackFitResult* tfr = theTK[0]->getTrackFitResultWithClosestMass(Const::pion);
600  if (tfr) {
601  svdTrkd0 = tfr->getD0();
602  svdTrkz0 = tfr->getZ0();
603  svdTrkpT = tfr->getMomentum().Perp();
604  TLorentzVector pStar = tfr->get4Momentum();
605  pStar.Boost(0, 0, 3. / 11);
606  svdTrkpCM = pStar.P();
607  }
608 
609  svdTrkPXDHits = (trk.getPXDHitList()).size();
610  svdTrkSVDHits = (trk.getSVDHitList()).size();
611  svdTrkCDCHits = (trk.getCDCHitList()).size();
612 
613  for (unsigned int i = 0; i < svdClusters.size(); i++) {
614  const SVDCluster* svd_1 = svdClusters[i];
615 
616  //get true hits, used only if isMC
617  RelationVector<SVDTrueHit> trueHit_1 = DataStore::getRelationsWithObj<SVDTrueHit>(svd_1);
618 
619  const RecoHitInformation* infoSVD_1 = trk.getRecoHitInformation(svd_1);
620  if (!infoSVD_1) {
621  continue;
622  }
623  const auto* hitTrackPoint_1 = trk.getCreatedTrackPoint(infoSVD_1);
624  const auto* fittedResult_1 = hitTrackPoint_1->getFitterInfo();
625  if (!fittedResult_1) {
626  continue;
627  }
628  const VxdID svd_id_1 = svd_1->getSensorID();
629  const unsigned short svd_Layer_1 = svd_id_1.getLayerNumber();
630  const unsigned short svd_Ladder_1 = svd_id_1.getLadderNumber();
631  const unsigned short svd_Sensor_1 = svd_id_1.getSensorNumber();
632  for (unsigned int l = i + 1; l < svdClusters.size(); l++) {
633  const SVDCluster* svd_2 = svdClusters[l];
634 
635  //get true hits, used only if isMC
636  RelationVector<SVDTrueHit> trueHit_2 = DataStore::getRelationsWithObj<SVDTrueHit>(svd_2);
637 
638  const RecoHitInformation* infoSVD_2 = trk.getRecoHitInformation(svd_2);
639  if (!infoSVD_2) {
640  continue;
641  }
642  const auto* hitTrackPoint_2 = trk.getCreatedTrackPoint(infoSVD_2);
643  const auto* fittedResult_2 = hitTrackPoint_2->getFitterInfo();
644  if (!fittedResult_2) {
645  continue;
646  }
647  const VxdID svd_id_2 = svd_2->getSensorID();
648  const unsigned short svd_Layer_2 = svd_id_2.getLayerNumber();
649  const unsigned short svd_Ladder_2 = svd_id_2.getLadderNumber();
650  const unsigned short svd_Sensor_2 = svd_id_2.getSensorNumber();
651  if (svd_Layer_1 == svd_Layer_2 && ((svd_Ladder_2 - svd_Ladder_1) == 1.0 || (svd_Ladder_2 - svd_Ladder_1) == -6.0
652  || (svd_Ladder_2 - svd_Ladder_1) == -9.0 || (svd_Ladder_2 - svd_Ladder_1) == -11.0 || (svd_Ladder_2 - svd_Ladder_1) == -15.0)) {
653  B2DEBUG(40, " ============= 2 hits in a SVD overlap ============= ");
654  const TVectorD resUnBias_SVD_1 = fittedResult_1->getResidual(0, false).getState();
655  const TVectorD resUnBias_SVD_2 = fittedResult_2->getResidual(0, false).getState();
656  genfit::MeasuredStateOnPlane state_unbiased_1 = fittedResult_1->getFittedState(false);
657  genfit::MeasuredStateOnPlane state_unbiased_2 = fittedResult_2->getFittedState(false);
658  const TVectorD& svd_predIntersect_unbiased_1 = state_unbiased_1.getState();
659  const TVectorD& svd_predIntersect_unbiased_2 = state_unbiased_2.getState();
660  const TMatrixDSym& covMatrix_unbiased_1 = state_unbiased_1.getCov();
661  const TMatrixDSym& covMatrix_unbiased_2 = state_unbiased_2.getCov();
662  genfit::MeasuredStateOnPlane state_1 = trk.getMeasuredStateOnPlaneFromRecoHit(infoSVD_1);
663  genfit::MeasuredStateOnPlane state_2 = trk.getMeasuredStateOnPlaneFromRecoHit(infoSVD_2);
664  const TVectorD& svd_predIntersect_1 = state_1.getState();
665  const TVectorD& svd_predIntersect_2 = state_2.getState();
666  const TMatrixDSym& covMatrix_1 = state_1.getCov();
667  const TMatrixDSym& covMatrix_2 = state_2.getCov();
668  //Restricting to consecutive SVD u-clusters
669  if (svd_1->isUCluster() == true && svd_2->isUCluster() == true) {
670  const int strips_1 = svd_1->getSize();
671  h_SVDstrips_Mult->Fill(strips_1);
672  const int strips_2 = svd_2->getSize();
673  h_SVDstrips_Mult->Fill(strips_2);
674  const double res_U_1 = resUnBias_SVD_1.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
675  const double res_U_2 = resUnBias_SVD_2.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
676  const float over_U_SVD = res_U_2 - res_U_1;
677  const TVector3 svdLocal_1(svd_1->getPosition(), svd_predIntersect_1[4], 0.);
678  const TVector3 svdLocal_2(svd_2->getPosition(), svd_predIntersect_2[4], 0.);
679  const VXD::SensorInfoBase& svdSensor_1 = geo.get(svd_id_1);
680  const VXD::SensorInfoBase& svdSensor_2 = geo.get(svd_id_2);
681  const TVector3& svdGlobal_1 = svdSensor_1.pointToGlobal(svdLocal_1);
682  const TVector3& svdGlobal_2 = svdSensor_2.pointToGlobal(svdLocal_2);
683  double svdPhi_1 = atan2(svdGlobal_1(1), svdGlobal_1(0));
684  double svdPhi_2 = atan2(svdGlobal_2(1), svdGlobal_2(0));
685  double svdZ_1 = svdGlobal_1(2);
686  double svdZ_2 = svdGlobal_2(2);
687  B2DEBUG(40, "SVD: difference of u-residuals =========> " << over_U_SVD);
688  //Fill SVD tree for u-overlaps if required by the user
689  if (m_ExpertLevel) {
690  svdDeltaRes_U = over_U_SVD;
691  svdRes_U_int = res_U_1;
692  svdClTime_U_int = svd_1->getClsTime();
693  svdClSNR_U_int = svd_1->getSNR();
694  svdClCharge_U_int = svd_1->getCharge();
695  svdClPosErr_U_int = svd_1->getPositionSigma();
696  if (isMC && trueHit_1.size() > 0)
697  svdTruePos_U_int = trueHit_1[0]->getU();
698  else
699  svdTruePos_U_int = -99;
700  svdClPhi_U_int = svdPhi_1;
701  svdClZ_U_int = svdZ_1;
702  svdTrkPos_U_int = svd_predIntersect_1[3];
703  svdTrkPosOS_U_int = svd_predIntersect_1[4];
704  svdTrkPosErr_U_int = sqrt(covMatrix_1[3][3]);
705  svdTrkPosErrOS_U_int = sqrt(covMatrix_1[4][4]);
706  svdTrkQoP_U_int = svd_predIntersect_1[0];
707  svdTrkPrime_U_int = svd_predIntersect_1[1];
708  svdTrkPrimeOS_U_int = svd_predIntersect_1[2];
709  svdTrkTraversedLength_U_int = svdSensor_1.getThickness() * sqrt(1 + svdTrkPrimeOS_U_int * svdTrkPrimeOS_U_int + svdTrkPrime_U_int *
710  svdTrkPrime_U_int);
711  svdTrkPosUnbiased_U_int = svd_predIntersect_unbiased_1[3];
712  svdClPos_U_int = svdRes_U_int / 1e4 + svdTrkPosUnbiased_U_int;
713  svdTrkPosErrUnbiased_U_int = sqrt(covMatrix_unbiased_1[3][3]);
714  svdTrkQoPUnbiased_U_int = svd_predIntersect_unbiased_1[0];
715  svdTrkPrimeUnbiased_U_int = svd_predIntersect_unbiased_1[1];
716  svdLayer_U_int = svd_Layer_1;
717  svdLadder_U_int = svd_Ladder_1;
718  svdSensor_U_int = svd_Sensor_1;
719  svdSize_U_int = strips_1;
720 
721  float pitch = 50e-4;
722  float halfLength = 1.92;
723  if (svdLayer_U_int > 3) {
724  pitch = 75e-4;
725  halfLength = 2.88;
726  }
727  svdClIntStrPos_U_int = fmod(svdClPos_U_int + halfLength, pitch) / pitch;
728 
729  svdStripCharge_U_int.clear();
730  svdStrip6Samples_U_int.clear();
731  svdStripTime_U_int.clear();
732  svdStripPosition_U_int.clear();
733  //retrieve relations and set strip charges and times
734  RelationVector<SVDRecoDigit> theRecoDigits_1 = DataStore::getRelationsWithObj<SVDRecoDigit>(svd_1);
735  if ((theRecoDigits_1.size() != svdSize_U_int) && (svdSize_U_int != 128)) //virtual cluster
736  B2ERROR(" Inconsistency with cluster size! # recoDigits = " << theRecoDigits_1.size() << " != " << svdSize_U_int <<
737  " cluster size");
738 
739  //skip clusters created beacuse of missing APV
740  if (svdSize_U_int < 128)
741  for (unsigned int d = 0; d < svdSize_U_int; d++) {
742  svdStripCharge_U_int.push_back(theRecoDigits_1[d]->getCharge());
743  SVDShaperDigit* ShaperDigit_1 = theRecoDigits_1[d]->getRelated<SVDShaperDigit>();
744  array<float, 6> Samples_1 = ShaperDigit_1->getSamples();
745  std::copy(std::begin(Samples_1), std::end(Samples_1), std::back_inserter(svdStrip6Samples_U_int));
746  svdStripTime_U_int.push_back(theRecoDigits_1[d]->getTime());
747  double misalignedStripPos = svdSensor_1.getUCellPosition(theRecoDigits_1[d]->getCellID());
748  //aligned strip pos = misaligned strip - ( misaligned cluster - aligned cluster)
749  svdStripPosition_U_int.push_back(misalignedStripPos - svd_1->getPosition() + svdClPos_U_int);
750  }
751 
752  svdRes_U_ext = res_U_2;
753  svdClTime_U_ext = svd_2->getClsTime();
754  svdClSNR_U_ext = svd_2->getSNR();
755  svdClCharge_U_ext = svd_2->getCharge();
756  svdClPosErr_U_ext = svd_2->getPositionSigma();
757  if (isMC && trueHit_2.size() > 0)
758  svdTruePos_U_ext = trueHit_2[0]->getU();
759  else
760  svdTruePos_U_ext = -99;
761  svdClPhi_U_ext = svdPhi_2;
762  svdClZ_U_ext = svdZ_2;
763  svdTrkPos_U_ext = svd_predIntersect_2[3];
764  svdTrkPosOS_U_ext = svd_predIntersect_2[4];
765  svdTrkPosErr_U_ext = sqrt(covMatrix_2[3][3]);
766  svdTrkPosErrOS_U_ext = sqrt(covMatrix_2[4][4]);
767  svdTrkQoP_U_ext = svd_predIntersect_2[0];
768  svdTrkPrime_U_ext = svd_predIntersect_2[1];
769  svdTrkPrimeOS_U_ext = svd_predIntersect_2[2];
770  svdTrkTraversedLength_U_ext = svdSensor_2.getThickness() * sqrt(1 + svdTrkPrimeOS_U_ext * svdTrkPrimeOS_U_ext + svdTrkPrime_U_ext *
771  svdTrkPrime_U_ext);
772  svdTrkPosUnbiased_U_ext = svd_predIntersect_unbiased_2[3];
773  svdClPos_U_ext = svdRes_U_ext / 1e4 + svdTrkPosUnbiased_U_ext;
774  svdTrkPosErrUnbiased_U_ext = sqrt(covMatrix_unbiased_2[3][3]);
775  svdTrkQoPUnbiased_U_ext = svd_predIntersect_unbiased_2[0];
776  svdTrkPrimeUnbiased_U_ext = svd_predIntersect_unbiased_2[1];
777  svdLayer_U_ext = svd_Layer_2;
778  svdLadder_U_ext = svd_Ladder_2;
779  svdSensor_U_ext = svd_Sensor_2;
780  svdSize_U_ext = strips_2;
781 
782  svdClIntStrPos_U_ext = fmod(svdClPos_U_ext + halfLength, pitch) / pitch;
783 
784  svdStripCharge_U_ext.clear();
785  svdStrip6Samples_U_ext.clear();
786  svdStripTime_U_ext.clear();
787  svdStripPosition_U_ext.clear();
788  //retrieve relations and set strip charges and times
789  RelationVector<SVDRecoDigit> theRecoDigits_2 = DataStore::getRelationsWithObj<SVDRecoDigit>(svd_2);
790  if ((theRecoDigits_2.size() != svdSize_U_ext) && (svdSize_U_ext != 128)) //virtual cluster
791  B2ERROR(" Inconsistency with cluster size! # recoDigits = " << theRecoDigits_2.size() << " != " << svdSize_U_ext <<
792  " cluster size");
793  //skip clusters created beacuse of missing APV
794  if (svdSize_U_ext < 128)
795  for (unsigned int d = 0; d < svdSize_U_ext; d++) {
796  svdStripCharge_U_ext.push_back(theRecoDigits_2[d]->getCharge());
797  SVDShaperDigit* ShaperDigit_2 = theRecoDigits_2[d]->getRelated<SVDShaperDigit>();
798  array<float, 6> Samples_2 = ShaperDigit_2->getSamples();
799  std::copy(std::begin(Samples_2), std::end(Samples_2), std::back_inserter(svdStrip6Samples_U_ext));
800  svdStripTime_U_ext.push_back(theRecoDigits_2[d]->getTime());
801  double misalignedStripPos = svdSensor_2.getUCellPosition(theRecoDigits_2[d]->getCellID());
802  //aligned strip pos = misaligned strip - ( misaligned cluster - aligned cluster)
803  svdStripPosition_U_ext.push_back(misalignedStripPos - svd_2->getPosition() + svdClPos_U_ext);
804  }
805  t_SVD_U->Fill();
806  }
807  //Fill histograms of residuals differences with SVD u clusters
808  h_U_DeltaRes->Fill(over_U_SVD);
809  h_U_DeltaRes_SVD->Fill(over_U_SVD);
810  if (svd_1->getSize() < 3 && svd_2->getSize() < 3) { //Consider only clusters with 3 or less strips involved
811  h_U_Cl1Cl2_DeltaRes[svd_1->getSize() * svd_2->getSize()]->Fill(over_U_SVD);
812  }
813  //Fill sensor hit-maps and 2D histograms with SVD u clusters
814  if (svd_Layer_1 == 3 && svd_Layer_2 == 3) {
815  h_Lyr3[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
816  h_Lyr3[svd_Ladder_2][svd_Sensor_2]->Fill(svd_2->getPosition(), 0.);
817  h_U_DeltaRes_SVD_Lyr3->Fill(over_U_SVD);
818  h_DeltaResUPhi_Lyr3->Fill(svdPhi_1, over_U_SVD);
819  h_DeltaResUPhi_Lyr3->Fill(svdPhi_2, over_U_SVD);
820  h_DeltaResUz_Lyr3->Fill(svdZ_1, over_U_SVD);
821  h_DeltaResUz_Lyr3->Fill(svdZ_2, over_U_SVD);
822  }
823  if (svd_Layer_1 == 4 && svd_Layer_2 == 4) {
824  h_Lyr4[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
825  h_Lyr4[svd_Ladder_2][svd_Sensor_2]->Fill(svd_2->getPosition(), 0.);
826  h_U_DeltaRes_SVD_Lyr4->Fill(over_U_SVD);
827  h_DeltaResUPhi_Lyr4->Fill(svdPhi_1, over_U_SVD);
828  h_DeltaResUPhi_Lyr4->Fill(svdPhi_2, over_U_SVD);
829  h_DeltaResUz_Lyr4->Fill(svdZ_1, over_U_SVD);
830  h_DeltaResUz_Lyr4->Fill(svdZ_2, over_U_SVD);
831  }
832  if (svd_Layer_1 == 5 && svd_Layer_2 == 5) {
833  h_Lyr5[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
834  h_Lyr5[svd_Ladder_2][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
835  h_U_DeltaRes_SVD_Lyr5->Fill(over_U_SVD);
836  h_DeltaResUPhi_Lyr5->Fill(svdPhi_1, over_U_SVD);
837  h_DeltaResUPhi_Lyr5->Fill(svdPhi_2, over_U_SVD);
838  h_DeltaResUz_Lyr5->Fill(svdZ_1, over_U_SVD);
839  h_DeltaResUz_Lyr5->Fill(svdZ_2, over_U_SVD);
840  }
841  if (svd_Layer_1 == 6 && svd_Layer_2 == 6) {
842  h_Lyr6[svd_Ladder_1][svd_Sensor_1]->Fill(svd_1->getPosition(), 0.);
843  h_Lyr6[svd_Ladder_2][svd_Sensor_2]->Fill(svd_2->getPosition(), 0.);
844  h_U_DeltaRes_SVD_Lyr6->Fill(over_U_SVD);
845  h_DeltaResUPhi_Lyr6->Fill(svdPhi_1, over_U_SVD);
846  h_DeltaResUPhi_Lyr6->Fill(svdPhi_2, over_U_SVD);
847  h_DeltaResUz_Lyr6->Fill(svdZ_1, over_U_SVD);
848  h_DeltaResUz_Lyr6->Fill(svdZ_2, over_U_SVD);
849  }
850  }
851  //Restricting to consecutive SVD v-clusters
852  if (svd_1->isUCluster() != true && svd_2->isUCluster() != true) {
853  const int strips_1 = svd_1->getSize();
854  h_SVDstrips_Mult->Fill(strips_1);
855  const int strips_2 = svd_2->getSize();
856  h_SVDstrips_Mult->Fill(strips_2);
857  const double res_V_1 = resUnBias_SVD_1.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
858  const double res_V_2 = resUnBias_SVD_2.GetMatrixArray()[0] * Unit::convertValueToUnit(1.0, "um");
859  const float over_V_SVD = res_V_2 - res_V_1;
860  const TVector3 svdLocal_1(svd_predIntersect_1[3], svd_1->getPosition(), 0.);
861  const TVector3 svdLocal_2(svd_predIntersect_2[3], svd_2->getPosition(), 0.);
862  const VXD::SensorInfoBase& svdSensor_1 = geo.get(svd_id_1);
863  const VXD::SensorInfoBase& svdSensor_2 = geo.get(svd_id_2);
864  const TVector3& svdGlobal_1 = svdSensor_1.pointToGlobal(svdLocal_1);
865  const TVector3& svdGlobal_2 = svdSensor_2.pointToGlobal(svdLocal_2);
866  double svdPhi_1 = atan2(svdGlobal_1(1), svdGlobal_1(0));
867  double svdPhi_2 = atan2(svdGlobal_2(1), svdGlobal_2(0));
868  double svdZ_1 = svdGlobal_1(2);
869  double svdZ_2 = svdGlobal_2(2);
870  B2DEBUG(40, "SVD: difference of v-residuals =========> " << over_V_SVD);
871  //Fill SVD tree for v-overlaps if required by the user
872  if (m_ExpertLevel) {
873  svdDeltaRes_V = over_V_SVD;
874  svdRes_V_int = res_V_1;
875  svdClTime_V_int = svd_1->getClsTime();
876  svdClSNR_V_int = svd_1->getSNR();
877  svdClCharge_V_int = svd_1->getCharge();
878  svdClPosErr_V_int = svd_1->getPositionSigma();
879  if (isMC && trueHit_1.size() > 0)
880  svdTruePos_V_int = trueHit_1[0]->getV();
881  else
882  svdTruePos_V_int = -99;
883  svdClPhi_V_int = svdPhi_1;
884  svdClZ_V_int = svdZ_1;
885  svdTrkPos_V_int = svd_predIntersect_1[4];
886  svdTrkPosOS_V_int = svd_predIntersect_1[3];
887  svdTrkPosErr_V_int = sqrt(covMatrix_1[4][4]);
888  svdTrkPosErrOS_V_int = sqrt(covMatrix_1[3][3]);
889  svdTrkQoP_V_int = svd_predIntersect_1[0];
890  svdTrkPrime_V_int = svd_predIntersect_1[2];
891  svdTrkPrimeOS_V_int = svd_predIntersect_1[1];
892  svdTrkTraversedLength_V_int = svdSensor_1.getThickness() * sqrt(1 + svdTrkPrimeOS_V_int * svdTrkPrimeOS_V_int + svdTrkPrime_V_int *
893  svdTrkPrime_V_int);
894  svdTrkPosUnbiased_V_int = svd_predIntersect_unbiased_1[4];
895  svdClPos_V_int = svdRes_V_int / 1e4 + svdTrkPosUnbiased_V_int;
896  svdTrkPosErrUnbiased_V_int = sqrt(covMatrix_unbiased_1[4][4]);
897  svdTrkQoPUnbiased_V_int = svd_predIntersect_unbiased_1[0];
898  svdTrkPrimeUnbiased_V_int = svd_predIntersect_unbiased_1[2];
899  svdLayer_V_int = svd_Layer_1;
900  svdLadder_V_int = svd_Ladder_1;
901  svdSensor_V_int = svd_Sensor_1;
902  svdSize_V_int = strips_1;
903 
904  float pitch = 160e-4;
905  float halfLength = 6.144;
906  if (svdLayer_V_int > 3) {
907  pitch = 240e-4;
908  }
909  svdClIntStrPos_V_int = fmod(svdClPos_V_int + halfLength, pitch) / pitch;
910 
911  svdStripCharge_V_int.clear();
912  svdStrip6Samples_V_int.clear();
913  svdStripTime_V_int.clear();
914  svdStripPosition_V_int.clear();
915  //retrieve relations and set strip charges and times
916  RelationVector<SVDRecoDigit> theRecoDigits_1 = DataStore::getRelationsWithObj<SVDRecoDigit>(svd_1);
917  if ((theRecoDigits_1.size() != svdSize_V_int) && (svdSize_V_int != 128)) //virtual cluster
918  B2ERROR(" Inconsistency with cluster size! # recoDigits = " << theRecoDigits_1.size() << " != " << svdSize_V_int <<
919  " cluster size");
920  //skip clusters created beacuse of missing APV
921  if (svdSize_V_int < 128)
922  for (unsigned int d = 0; d < svdSize_V_int; d++) {
923  svdStripCharge_V_int.push_back(theRecoDigits_1[d]->getCharge());
924  SVDShaperDigit* ShaperDigit_1 = theRecoDigits_1[d]->getRelated<SVDShaperDigit>();
925  array<float, 6> Samples_1 = ShaperDigit_1->getSamples();
926  std::copy(std::begin(Samples_1), std::end(Samples_1), std::back_inserter(svdStrip6Samples_V_int));
927  svdStripTime_V_int.push_back(theRecoDigits_1[d]->getTime());
928  double misalignedStripPos = svdSensor_1.getVCellPosition(theRecoDigits_1[d]->getCellID());
929  //aligned strip pos = misaligned strip - ( misaligned cluster - aligned cluster)
930  svdStripPosition_V_int.push_back(misalignedStripPos - svd_1->getPosition() + svdClPos_V_int);
931  }
932 
933  svdRes_V_ext = res_V_2;
934  svdClTime_V_ext = svd_2->getClsTime();
935  svdClSNR_V_ext = svd_2->getSNR();
936  svdClCharge_V_ext = svd_2->getCharge();
937  svdClPosErr_V_ext = svd_2->getPositionSigma();
938  if (isMC && trueHit_2.size() > 0)
939  svdTruePos_V_ext = trueHit_2[0]->getV();
940  else
941  svdTruePos_V_ext = -99;
942  svdClPhi_V_ext = svdPhi_2;
943  svdClZ_V_ext = svdZ_2;
944  svdTrkPos_V_ext = svd_predIntersect_2[4];
945  svdTrkPosOS_V_ext = svd_predIntersect_2[3];
946  svdTrkPosErr_V_ext = sqrt(covMatrix_2[4][4]);
947  svdTrkPosErrOS_V_ext = sqrt(covMatrix_1[3][3]);
948  svdTrkQoP_V_ext = svd_predIntersect_2[0];
949  svdTrkPrime_V_ext = svd_predIntersect_2[2];
950  svdTrkPrimeOS_V_ext = svd_predIntersect_2[1];
951  svdTrkTraversedLength_V_ext = svdSensor_2.getThickness() * sqrt(1 + svdTrkPrimeOS_V_ext * svdTrkPrimeOS_V_ext + svdTrkPrime_V_ext *
952  svdTrkPrime_V_ext);
953  svdTrkPosUnbiased_V_ext = svd_predIntersect_unbiased_2[4];
954  svdClPos_V_ext = svdRes_V_ext / 1e4 + svdTrkPosUnbiased_V_ext;
955  svdTrkPosErrUnbiased_V_ext = sqrt(covMatrix_unbiased_2[4][4]);
956  svdTrkQoPUnbiased_V_ext = svd_predIntersect_unbiased_2[0];
957  svdTrkPrimeUnbiased_V_ext = svd_predIntersect_unbiased_2[2];
958  svdLayer_V_ext = svd_Layer_2;
959  svdLadder_V_ext = svd_Ladder_2;
960  svdSensor_V_ext = svd_Sensor_2;
961  svdSize_V_ext = strips_2;
962 
963  svdClIntStrPos_V_ext = fmod(svdClPos_V_ext + halfLength, pitch) / pitch;
964 
965  svdStripCharge_V_ext.clear();
966  svdStrip6Samples_V_ext.clear();
967  svdStripTime_V_ext.clear();
968  svdStripPosition_V_ext.clear();
969  //retrieve relations and set strip charges and times
970  RelationVector<SVDRecoDigit> theRecoDigits_2 = DataStore::getRelationsWithObj<SVDRecoDigit>(svd_2);
971  if ((theRecoDigits_2.size() != svdSize_V_ext) && (svdSize_V_ext != 128)) //virtual cluster
972  B2ERROR(" Inconsistency with cluster size! # recoDigits = " << theRecoDigits_2.size() << " != " << svdSize_V_ext <<
973  " cluster size");
974  //skip clusters created beacuse of missing APV
975  if (svdSize_V_ext < 128)
976  for (unsigned int d = 0; d < svdSize_V_ext; d++) {
977  svdStripCharge_V_ext.push_back(theRecoDigits_2[d]->getCharge());
978  SVDShaperDigit* ShaperDigit_2 = theRecoDigits_2[d]->getRelated<SVDShaperDigit>();
979  array<float, 6> Samples_2 = ShaperDigit_2->getSamples();
980  std::copy(std::begin(Samples_2), std::end(Samples_2), std::back_inserter(svdStrip6Samples_V_ext));
981  svdStripTime_V_ext.push_back(theRecoDigits_2[d]->getTime());
982  double misalignedStripPos = svdSensor_2.getVCellPosition(theRecoDigits_2[d]->getCellID());
983  //aligned strip pos = misaligned strip - ( misaligned cluster - aligned cluster)
984  svdStripPosition_V_ext.push_back(misalignedStripPos - svd_2->getPosition() + svdClPos_V_ext);
985  }
986 
987  t_SVD_V->Fill();
988  }
989  //Fill histograms of residuals differences with SVD v clusters
990  h_V_DeltaRes->Fill(over_V_SVD);
991  h_V_DeltaRes_SVD->Fill(over_V_SVD);
992  if (svd_1->getSize() < 3 && svd_2->getSize() < 3) { //Consider only clusters with 3 or less strips involved
993  h_V_Cl1Cl2_DeltaRes[svd_1->getSize() * svd_2->getSize()]->Fill(over_V_SVD);
994  }
995  //Fill sensor hit-maps and 2D histograms with SVD v clusters
996  if (svd_Layer_1 == 3 && svd_Layer_2 == 3) {
997  h_Lyr3[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
998  h_Lyr3[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
999  h_V_DeltaRes_SVD_Lyr3->Fill(over_V_SVD);
1000  h_DeltaResVPhi_Lyr3->Fill(svdPhi_1, over_V_SVD);
1001  h_DeltaResVPhi_Lyr3->Fill(svdPhi_2, over_V_SVD);
1002  h_DeltaResVz_Lyr3->Fill(svdZ_1, over_V_SVD);
1003  h_DeltaResVz_Lyr3->Fill(svdZ_2, over_V_SVD);
1004  }
1005  if (svd_Layer_1 == 4 && svd_Layer_2 == 4) {
1006  h_Lyr4[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
1007  h_Lyr4[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
1008  h_V_DeltaRes_SVD_Lyr4->Fill(over_V_SVD);
1009  h_DeltaResVPhi_Lyr4->Fill(svdPhi_1, over_V_SVD);
1010  h_DeltaResVPhi_Lyr4->Fill(svdPhi_2, over_V_SVD);
1011  h_DeltaResVz_Lyr4->Fill(svdZ_1, over_V_SVD);
1012  h_DeltaResVz_Lyr4->Fill(svdZ_2, over_V_SVD);
1013  }
1014  if (svd_Layer_1 == 5 && svd_Layer_2 == 5) {
1015  h_Lyr5[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
1016  h_Lyr5[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
1017  h_V_DeltaRes_SVD_Lyr5->Fill(over_V_SVD);
1018  h_DeltaResVPhi_Lyr5->Fill(svdPhi_1, over_V_SVD);
1019  h_DeltaResVPhi_Lyr5->Fill(svdPhi_2, over_V_SVD);
1020  h_DeltaResVz_Lyr5->Fill(svdZ_1, over_V_SVD);
1021  h_DeltaResVz_Lyr5->Fill(svdZ_2, over_V_SVD);
1022  }
1023  if (svd_Layer_1 == 6 && svd_Layer_2 == 6) {
1024  h_Lyr6[svd_Ladder_1][svd_Sensor_1]->Fill(0., svd_1->getPosition());
1025  h_Lyr6[svd_Ladder_2][svd_Sensor_2]->Fill(0., svd_2->getPosition());
1026  h_V_DeltaRes_SVD_Lyr6->Fill(over_V_SVD);
1027  h_DeltaResVPhi_Lyr6->Fill(svdPhi_1, over_V_SVD);
1028  h_DeltaResVPhi_Lyr6->Fill(svdPhi_2, over_V_SVD);
1029  h_DeltaResVz_Lyr6->Fill(svdZ_1, over_V_SVD);
1030  h_DeltaResVz_Lyr6->Fill(svdZ_2, over_V_SVD);
1031  }
1032  }
1033  }
1034  }
1035  }
1036  }
1037 }
1038 
1039 
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:542
Belle2::Environment::isMC
bool isMC() const
Do we have generated, not real data?
Definition: Environment.cc:55
Belle2::Environment::Instance
static Environment & Instance()
Static method to get a reference to the Environment instance.
Definition: Environment.cc:29
Belle2::HistoModule
HistoModule.h is supposed to be used instead of Module.h for the modules with histogram definitions t...
Definition: HistoModule.h:29
Belle2::OverlapResidualsModule
The module studies VXD hits from overlapping sensors of a same VXD layer.
Definition: OverlapResidualsModule.h:34
Belle2::OverlapResidualsModule::initialize
void initialize() override
Register input and output data.
Definition: OverlapResidualsModule.cc:56
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:468
Belle2::OverlapResidualsModule::defineHisto
void defineHisto() override
Histogram definitions such as TH1(), TH2(), TNtuple(), TTree()....
Definition: OverlapResidualsModule.cc:64
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:30
Belle2::PXDCluster::getV
float getV() const
Get v coordinate of hit position.
Definition: PXDCluster.h:136
Belle2::PXDCluster::getSensorID
VxdID getSensorID() const
Get the sensor ID.
Definition: PXDCluster.h:126
Belle2::PXDCluster::getU
float getU() const
Get u coordinate of hit position.
Definition: PXDCluster.h:131
Belle2::RecoHitInformation
This class stores additional information to every CDC/SVD/PXD hit stored in a RecoTrack.
Definition: RecoHitInformation.h:38
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: RelationVector.h:67
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:88
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:28
Belle2::SVDCluster::getClsTime
float getClsTime() const
Get average of waveform maximum times of cluster strip signals.
Definition: SVDCluster.h:133
Belle2::SVDCluster::getSNR
float getSNR() const
Get cluster SNR.
Definition: SVDCluster.h:158
Belle2::SVDCluster::getSize
unsigned short getSize() const
Get cluster size.
Definition: SVDCluster.h:153
Belle2::SVDCluster::getCharge
float getCharge() const
Get collected charge.
Definition: SVDCluster.h:143
Belle2::SVDCluster::getSensorID
VxdID getSensorID() const
Get the sensor ID.
Definition: SVDCluster.h:101
Belle2::SVDCluster::isUCluster
bool isUCluster() const
Get the direction of strips.
Definition: SVDCluster.h:109
Belle2::SVDCluster::getPosition
float getPosition(double v=0) const
Get the coordinate of reconstructed hit.
Definition: SVDCluster.h:116
Belle2::SVDCluster::getPositionSigma
float getPositionSigma() const
Get the error of the reconstructed hit coordinate.
Definition: SVDCluster.h:128
Belle2::SVDShaperDigit
The SVD ShaperDigit class.
Definition: SVDShaperDigit.h:32
Belle2::SVDShaperDigit::getSamples
APVFloatSamples getSamples() const
Get array of samples.
Definition: SVDShaperDigit.h:119
Belle2::StoreAccessorBase::isOptional
bool isOptional(const std::string &name="")
Tell the DataStore about an optional input.
Definition: StoreAccessorBase.h:93
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:49
Belle2::TrackFitResult::get4Momentum
TLorentzVector get4Momentum() const
Getter for the 4Momentum at the closest approach of the track in the r/phi projection.
Definition: TrackFitResult.h:115
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:106
Belle2::TrackFitResult::getD0
double getD0() const
Getter for d0.
Definition: TrackFitResult.h:168
Belle2::TrackFitResult::getZ0
double getZ0() const
Getter for z0.
Definition: TrackFitResult.h:190
Belle2::VXD::GeoCache
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
Definition: GeoCache.h:39
Belle2::VXD::GeoCache::getInstance
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition: GeoCache.cc:213
Belle2::VXD::GeoCache::get
static const SensorInfoBase & get(Belle2::VxdID id)
Return a reference to the SensorInfo of a given SensorID.
Definition: GeoCache.h:139
Belle2::VXD::SensorInfoBase
Base class to provide Sensor Information for PXD and SVD.
Definition: SensorInfoBase.h:30
Belle2::VXD::SensorInfoBase::getVCellPosition
double getVCellPosition(int vID) const
Return the position of a specific strip/pixel in v direction.
Definition: SensorInfoBase.h:181
Belle2::VXD::SensorInfoBase::getUCellPosition
double getUCellPosition(int uID, int vID=-1) const
Return the position of a specific strip/pixel in u direction.
Definition: SensorInfoBase.h:169
Belle2::VXD::SensorInfoBase::pointToGlobal
TVector3 pointToGlobal(const TVector3 &local, bool reco=false) const
Convert a point from local to global coordinates.
Definition: SensorInfoBase.h:365
Belle2::VXD::SensorInfoBase::getThickness
double getThickness() const
Return the thickness of the sensor.
Definition: SensorInfoBase.h:111
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:100
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:98
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:96
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
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:139
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17