35 gROOT->LoadMacro(
"$VMCWORKDIR/macro/run/Tools.C");
38 gRandom->SetSeed(ranseed);
41 fMcCands =
new TClonesArray(
"RhoCandidate");
42 fCands =
new TClonesArray(
"RhoCandidate");
50 double min=0,
max=4.5,mdf=2.0;
51 double min2=0.13,max2=0.15;
52 double min3=1.,max3=1.5;
58 double angdiffran=0.05;
60 TList* histolist =
new TList();
63 TH1F *hVtxDiffX=
new TH1F(
"hvtxDiffx",
"Vertex Smearing;(x_{mc}-x)/cm",bins,-vrange,vrange);
64 TH1F *hVtxDiffY=
new TH1F(
"hvtxDiffy",
"Vertex Smearing;(y_{mc}-y)/cm",bins,-vrange,vrange);
65 TH1F *hVtxDiffZ=
new TH1F(
"hvtxDiffz",
"Vertex Smearing;(z_{mc}-z)/cm",bins,-vrange,vrange);
67 TH1F *hVtxDiffPX=
new TH1F(
"hvtxdiffpx",
"Momentum reco;(p_{x,mc}-p_{x} / GeV/c)",bins,-momrange,momrange);
68 TH1F *hVtxDiffPY=
new TH1F(
"hvtxdiffpy",
"Momentum reco;(p_{y,mc}-p_{y} / GeV/c)",bins,-momrange,momrange);
69 TH1F *hVtxDiffPZ=
new TH1F(
"hvtxdiffpz",
"Momentum reco;(p_{z,mc}-p_{z}) / GeV/c",bins,-momrange,momrange);
70 TH1F *hVtxDiffE=
new TH1F(
"hvtxdiffe",
"Energy reco;(E_{mc}-E) / GeV",bins,-momrange,momrange);
72 TH1F *hVtxPullPX=
new TH1F(
"hvtxpullpx",
"Momentum pull distribution reco;(p_{x,mc}-p_{x})/#sigma_{p_{x}}",bins,-pullrange,pullrange);
73 TH1F *hVtxPullPY=
new TH1F(
"hvtxpullpy",
"Momentum pull distribution reco;(p_{y,mc}-p_{y})/#sigma_{p_{y}}",bins,-pullrange,pullrange);
74 TH1F *hVtxPullPZ=
new TH1F(
"hvtxpullpz",
"Momentum pull distribution reco;(p_{z,mc}-p_{z})/#sigma_{p_{z}}",bins,-pullrange,pullrange);
75 TH1F *hVtxPullE=
new TH1F(
"hvtxpulle",
"Energy pull distribution reco;(E_{mc}-E)/#sigma_{E}",bins,-pullrange,pullrange);
77 TH1F *hVtxDiffThe=
new TH1F(
"hvtxdiffthe",
"Momentum reco;(#Theta_{mc}-#Theta}) / GeV/c",bins,-angdiffran,angdiffran);
78 TH1F *hVtxDiffPhi=
new TH1F(
"hvtxdiffphi",
"Momentum reco;(#Phi_{mc}-#Phi}) / GeV/c",bins,-angdiffran,angdiffran);
82 TH1F *hVtxPocaX=
new TH1F(
"hvtxPocax",
"Vertex POCA;(x_{mc}-x)/cm",bins,-vrange,vrange);
83 TH1F *hVtxPocaY=
new TH1F(
"hvtxPocay",
"Vertex POCA;(y_{mc}-y)/cm",bins,-vrange,vrange);
84 TH1F *hVtxPocaZ=
new TH1F(
"hvtxPocaz",
"Vertex POCA;(z_{mc}-z)/cm",bins,-vrange,vrange);
86 TH2F *hVtxPocaXY=
new TH2F(
"hvtxPocaxy",
"Vertex POCA X-Y view;(x_{mc}-x)/#mum;(y_{mc}-y)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
87 TH2F *hVtxPocaRZ=
new TH2F(
"hvtxPocarz",
"Vertex POCA R-Z view;(z_{mc}-z)/#mum;(r_{mc}-r)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
89 TH1F *hVtxPullPocaX=
new TH1F(
"hvtxpullpocax",
"Vertex pull distribution POCA;(x_{mc}-x)/#sigma_{x}",bins,-pullrange,pullrange);
90 TH1F *hVtxPullPocaY=
new TH1F(
"hvtxpullpocay",
"Vertex pull distribution POCA;(y_{mc}-y)/#sigma_{y}",bins,-pullrange,pullrange);
91 TH1F *hVtxPullPocaZ=
new TH1F(
"hvtxpullpocaz",
"Vertex pull distribution POCA;(z_{mc}-z)/#sigma_{z}",bins,-pullrange,pullrange);
93 TH1F *hVtxPocas=
new TH1F(
"hvtxpocas",
"Vertex POCA D value;D/cm",bins,0.,vrange);
94 TH1F *hVtxPocaEmpty=
new TH1F(
"hvtxcpocaempty",
"",0,0.,1.);
98 TH1F *hVtxFastX=
new TH1F(
"hvtxFastx",
"Vertex fast fit;(x_{mc}-x)/cm",bins,-vrange,vrange);
99 TH1F *hVtxFastY=
new TH1F(
"hvtxTasty",
"Vertex fast fit;(y_{mc}-y)/cm",bins,-vrange,vrange);
100 TH1F *hVtxFastZ=
new TH1F(
"hvtxFastz",
"Vertex fast fit;(z_{mc}-z)/cm",bins,-vrange,vrange);
102 TH2F *hVtxFastXY=
new TH2F(
"hvtxfastxy",
"Vertex fast fit X-Y view;(x_{mc}-x)/#mum;(y_{mc}-y)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
103 TH2F *hVtxFastRZ=
new TH2F(
"hvtxfastrz",
"Vertex fast fit R-Z view;(z_{mc}-z)/#mum;(r_{mc}-r)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
105 TH1F *hVtxErrFastX=
new TH1F(
"hvtxerrfastx",
"Vertex error distribution fast;#sigma_{x}",bins,0.,vrange);
106 TH1F *hVtxErrFastY=
new TH1F(
"hvtxerrfasty",
"Vertex error distribution fast;#sigma_{y}",bins,0.,vrange);
107 TH1F *hVtxErrFastZ=
new TH1F(
"hvtxerrfastz",
"Vertex error distribution fast;#sigma_{z}",bins,0.,vrange);
108 TH1F *hVtxPullFastX=
new TH1F(
"hvtxpullfastx",
"Vertex pull distribution fast;(x_{mc}-x)/#sigma_{x}",bins,-pullrange,pullrange);
109 TH1F *hVtxPullFastY=
new TH1F(
"hvtxpullfasty",
"Vertex pull distribution fast;(y_{mc}-y)/#sigma_{y}",bins,-pullrange,pullrange);
110 TH1F *hVtxPullFastZ=
new TH1F(
"hvtxpullfastz",
"Vertex pull distribution fast;(z_{mc}-z)/#sigma_{z}",bins,-pullrange,pullrange);
112 TH1F *hVtxChi2Fast=
new TH1F(
"hvtxchi2fast",
"Vertex Fast #Chi^{2};#Chi^{2}",bins,0,maxchi);
113 TH1F *hVtxChiProbFast=
new TH1F(
"hvtxchipropfast",
"P(#Chi^{2}) Fast;p",bins,0.,1.);
116 TH1F *hVtxFitX=
new TH1F(
"hvtxFitx",
"Vertex full fit;(x_{mc}-x)/cm",bins,-vrange,vrange);
117 TH1F *hVtxFitY=
new TH1F(
"hvtxFity",
"Vertex full fit;(y_{mc}-y)/cm",bins,-vrange,vrange);
118 TH1F *hVtxFitZ=
new TH1F(
"hvtxFitz",
"Vertex full fit;(z_{mc}-z)/cm",bins,-vrange,vrange);
120 TH2F *hVtxFitXY=
new TH2F(
"hvtxfitxy",
"Vertex full fit X-Y view;(x_{mc}-x)/#mum;(y_{mc}-y)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
121 TH2F *hVtxFitRZ=
new TH2F(
"hvtxfitrz",
"Vertex full fit R-Z view;(z_{mc}-z)/#mum;(r_{mc}-r)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
123 TH1F *hVtxErrFitX=
new TH1F(
"hvtxerrfitx",
"Vertex error distribution full fit;#sigma_{x}",bins,0.,vrange);
124 TH1F *hVtxErrFitY=
new TH1F(
"hvtxerrfity",
"Vertex error distribution full fit;#sigma_{y}",bins,0.,vrange);
125 TH1F *hVtxErrFitZ=
new TH1F(
"hvtxerrfitz",
"Vertex error distribution full fit;#sigma_{z}",bins,0.,vrange);
126 TH1F *hVtxPullFitX=
new TH1F(
"hvtxpullfitx",
"Vertex pull distribution full fit;(x_{mc}-x)/#sigma_{x}",bins,-pullrange,pullrange);
127 TH1F *hVtxPullFitY=
new TH1F(
"hvtxpullfity",
"Vertex pull distribution full fit;(y_{mc}-y)/#sigma_{y}",bins,-pullrange,pullrange);
128 TH1F *hVtxPullFitZ=
new TH1F(
"hvtxpullfitz",
"Vertex pull distribution full fit;(z_{mc}-z)/#sigma_{z}",bins,-pullrange,pullrange);
129 TH1F *hVtxDiffFitPX=
new TH1F(
"hvtxdifffitpx",
"Momentum full fit;(p_{x,mc}-p_{x} / GeV/c)",bins,-momrange,momrange);
130 TH1F *hVtxDiffFitPY=
new TH1F(
"hvtxdifffitpy",
"Momentum full fit;(p_{y,mc}-p_{y} / GeV/c)",bins,-momrange,momrange);
131 TH1F *hVtxDiffFitPZ=
new TH1F(
"hvtxdifffitpz",
"Momentum full fit;(p_{z,mc}-p_{z}) / GeV/c",bins,-momrange,momrange);
132 TH1F *hVtxDiffFitE=
new TH1F(
"hvtxdifffite",
"Energy full fit;(E_{mc}-E) / GeV",bins,-momrange,momrange);
133 TH1F *hVtxDiffFitThe=
new TH1F(
"hvtxdifffitthe",
"Momentum full fit;(#Theta_{mc}-#Theta}) / GeV/c",bins,-angdiffran,angdiffran);
134 TH1F *hVtxDiffFitPhi=
new TH1F(
"hvtxdifffitphi",
"Momentum full fit;(#Phi_{mc}-#Phi}) / GeV/c",bins,-angdiffran,angdiffran);
136 TH1F *hVtxPullFitPX=
new TH1F(
"hvtxpullfitpx",
"Momentum pull distribution full fit;(p_{x,mc}-p_{x})/#sigma_{p_{x}}",bins,-pullrange,pullrange);
137 TH1F *hVtxPullFitPY=
new TH1F(
"hvtxpullfitpy",
"Momentum pull distribution full fit;(p_{y,mc}-p_{y})/#sigma_{p_{y}}",bins,-pullrange,pullrange);
138 TH1F *hVtxPullFitPZ=
new TH1F(
"hvtxpullfitpz",
"Momentum pull distribution full fit;(p_{z,mc}-p_{z})/#sigma_{p_{z}}",bins,-pullrange,pullrange);
139 TH1F *hVtxPullFitE=
new TH1F(
"hvtxpullfite",
"Energy pull distribution full fit;(E_{mc}-E)/#sigma_{E}",bins,-pullrange,pullrange);
141 TH1F *hVtxChi2Fit=
new TH1F(
"hvtxchi2fit",
"Vertex Fit #Chi^{2};#Chi^{2}",bins,0,maxchi);
142 TH1F *hVtxChiProbFit=
new TH1F(
"hvtxchipropfit",
"P(#Chi^{2}) Fit;p",bins,0.,1.);
147 TH1F *hVtxKinX=
new TH1F(
"hvtxKinx",
"Vertex full kin;(x_{mc}-x)/cm",bins,-vrange,vrange);
148 TH1F *hVtxKinY=
new TH1F(
"hvtxKiny",
"Vertex full kin;(y_{mc}-y)/cm",bins,-vrange,vrange);
149 TH1F *hVtxKinZ=
new TH1F(
"hvtxKinz",
"Vertex full kin;(z_{mc}-z)/cm",bins,-vrange,vrange);
151 TH2F *hVtxKinXY=
new TH2F(
"hvtxkinxy",
"Vertex full kin X-Y view;(x_{mc}-x)/#mum;(y_{mc}-y)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
152 TH2F *hVtxKinRZ=
new TH2F(
"hvtxkinrz",
"Vertex full kin R-Z view;(z_{mc}-z)/#mum;(r_{mc}-r)/cm",bins,-vrange,vrange,bins,-vrange,vrange);
154 TH1F *hVtxErrKinX=
new TH1F(
"hvtxerrkinx",
"Vertex error distribution full kin;#sigma_{x}",bins,0.,vrange);
155 TH1F *hVtxErrKinY=
new TH1F(
"hvtxerrkiny",
"Vertex error distribution full kin;#sigma_{y}",bins,0.,vrange);
156 TH1F *hVtxErrKinZ=
new TH1F(
"hvtxerrkinz",
"Vertex error distribution full kin;#sigma_{z}",bins,0.,vrange);
157 TH1F *hVtxPullKinX=
new TH1F(
"hvtxpullkinx",
"Vertex pull distribution full kin;(x_{mc}-x)/#sigma_{x}",bins,-pullrange,pullrange);
158 TH1F *hVtxPullKinY=
new TH1F(
"hvtxpullkiny",
"Vertex pull distribution full kin;(y_{mc}-y)/#sigma_{y}",bins,-pullrange,pullrange);
159 TH1F *hVtxPullKinZ=
new TH1F(
"hvtxpullkinz",
"Vertex pull distribution full kin;(z_{mc}-z)/#sigma_{z}",bins,-pullrange,pullrange);
160 TH1F *hVtxDiffKinPX=
new TH1F(
"hvtxdiffkinpx",
"Momentum full kin;(p_{x,mc}-p_{x} / GeV/c)",bins,-momrange,momrange);
161 TH1F *hVtxDiffKinPY=
new TH1F(
"hvtxdiffkinpy",
"Momentum full kin;(p_{y,mc}-p_{y} / GeV/c)",bins,-momrange,momrange);
162 TH1F *hVtxDiffKinPZ=
new TH1F(
"hvtxdiffkinpz",
"Momentum full kin;(p_{z,mc}-p_{z}) / GeV/c",bins,-momrange,momrange);
163 TH1F *hVtxDiffKinE=
new TH1F(
"hvtxdiffkine",
"Energy full kin;(E_{mc}-E) / GeV",bins,-momrange,momrange);
164 TH1F *hVtxDiffKinThe=
new TH1F(
"hvtxdiffkinthe",
"Momentum full kin;(#Theta_{mc}-#Theta}) / GeV/c",bins,-angdiffran,angdiffran);
165 TH1F *hVtxDiffKinPhi=
new TH1F(
"hvtxdiffkinphi",
"Momentum full kin;(#Phi_{mc}-#Phi}) / GeV/c",bins,-angdiffran,angdiffran);
167 TH1F *hVtxPullKinPX=
new TH1F(
"hvtxpullkinpx",
"Momentum pull distribution full kin;(p_{x,mc}-p_{x})/#sigma_{p_{x}}",bins,-pullrange,pullrange);
168 TH1F *hVtxPullKinPY=
new TH1F(
"hvtxpullkinpy",
"Momentum pull distribution full kin;(p_{y,mc}-p_{y})/#sigma_{p_{y}}",bins,-pullrange,pullrange);
169 TH1F *hVtxPullKinPZ=
new TH1F(
"hvtxpullkinpz",
"Momentum pull distribution full kin;(p_{z,mc}-p_{z})/#sigma_{p_{z}}",bins,-pullrange,pullrange);
170 TH1F *hVtxPullKinE=
new TH1F(
"hvtxpullkine",
"Energy pull distribution full kin;(E_{mc}-E)/#sigma_{E}",bins,-pullrange,pullrange);
172 TH1F *hVtxChi2Kin=
new TH1F(
"hvtxchi2kin",
"Vertex Kin #Chi^{2};#Chi^{2}",bins,0,maxchi);
173 TH1F *hVtxChiProbKin=
new TH1F(
"hvtxchipropkin",
"P(#Chi^{2}) Kin;p",bins,0.,1.);
176 TH1F *hPrgPull0 =
new TH1F(
"hprgpull0",
"PRG Pull Distribution (0)",bins,-pullrange,pullrange);
177 TH1F *hPrgPull1 =
new TH1F(
"hprgpull1",
"PRG Pull Distribution (1)",bins,-pullrange,pullrange);
178 TH1F *hPrgPull2 =
new TH1F(
"hprgpull2",
"PRG Pull Distribution (2)",bins,-pullrange,pullrange);
179 TH1F *hPrgPull3 =
new TH1F(
"hprgpull3",
"PRG Pull Distribution (3)",bins,-pullrange,pullrange);
180 TH1F *hPrgPull4 =
new TH1F(
"hprgpull4",
"PRG Pull Distribution (4)",bins,-pullrange,pullrange);
182 histolist->Add(hVtxDiffX);
183 histolist->Add(hVtxDiffY);
184 histolist->Add(hVtxDiffZ);
186 histolist->Add(hVtxPocaX);
187 histolist->Add(hVtxPocaY);
188 histolist->Add(hVtxPocaZ);
190 histolist->Add(hVtxFastX);
191 histolist->Add(hVtxFastY);
192 histolist->Add(hVtxFastZ);
194 histolist->Add(hVtxFitX);
195 histolist->Add(hVtxFitY);
196 histolist->Add(hVtxFitZ);
198 histolist->Add(hVtxKinX);
199 histolist->Add(hVtxKinY);
200 histolist->Add(hVtxKinZ);
202 histolist->Add(hVtxErrFastX);
203 histolist->Add(hVtxErrFastY);
204 histolist->Add(hVtxErrFastZ);
206 histolist->Add(hVtxErrFitX);
207 histolist->Add(hVtxErrFitY);
208 histolist->Add(hVtxErrFitZ);
210 histolist->Add(hVtxErrKinX);
211 histolist->Add(hVtxErrKinY);
212 histolist->Add(hVtxErrKinZ);
214 histolist->Add(hVtxPullPocaX);
215 histolist->Add(hVtxPullPocaY);
216 histolist->Add(hVtxPullPocaZ);
218 histolist->Add(hVtxPullFastX);
219 histolist->Add(hVtxPullFastY);
220 histolist->Add(hVtxPullFastZ);
222 histolist->Add(hVtxPullFitX);
223 histolist->Add(hVtxPullFitY);
224 histolist->Add(hVtxPullFitZ);
226 histolist->Add(hVtxPullKinX);
227 histolist->Add(hVtxPullKinY);
228 histolist->Add(hVtxPullKinZ);
230 histolist->Add(hVtxDiffPX);
231 histolist->Add(hVtxDiffPY);
232 histolist->Add(hVtxDiffPZ);
234 histolist->Add(hVtxPullPX);
235 histolist->Add(hVtxPullPY);
236 histolist->Add(hVtxPullPZ);
238 histolist->Add(hVtxDiffFitPX);
239 histolist->Add(hVtxDiffFitPY);
240 histolist->Add(hVtxDiffFitPZ);
242 histolist->Add(hVtxPullFitPX);
243 histolist->Add(hVtxPullFitPY);
244 histolist->Add(hVtxPullFitPZ);
246 histolist->Add(hVtxDiffKinPX);
247 histolist->Add(hVtxDiffKinPY);
248 histolist->Add(hVtxDiffKinPZ);
250 histolist->Add(hVtxPullKinPX);
251 histolist->Add(hVtxPullKinPY);
252 histolist->Add(hVtxPullKinPZ);
254 histolist->Add(hVtxDiffE);
255 histolist->Add(hVtxDiffFitE);
256 histolist->Add(hVtxDiffKinE);
257 histolist->Add(hVtxPullE);
258 histolist->Add(hVtxPullFitE);
259 histolist->Add(hVtxPullKinE);
261 histolist->Add(hVtxDiffThe);
262 histolist->Add(hVtxDiffPhi);
264 histolist->Add(hVtxDiffFitThe);
265 histolist->Add(hVtxDiffFitPhi);
267 histolist->Add(hVtxDiffKinThe);
268 histolist->Add(hVtxDiffKinPhi);
270 histolist->Add(hVtxPocas);
271 histolist->Add(hVtxChi2Fast);
272 histolist->Add(hVtxChi2Fit);
273 histolist->Add(hVtxChi2Kin);
274 histolist->Add(hVtxChiProbFast);
275 histolist->Add(hVtxChiProbFit);
276 histolist->Add(hVtxChiProbKin);
278 histolist->Add(hVtxPocaXY);
279 histolist->Add(hVtxFastXY);
280 histolist->Add(hVtxFitXY);
281 histolist->Add(hVtxKinXY);
282 histolist->Add(hVtxPocaRZ);
283 histolist->Add(hVtxFastRZ);
284 histolist->Add(hVtxFitRZ);
285 histolist->Add(hVtxKinRZ);
287 histolist->Add(hPrgPull0);
288 histolist->Add(hPrgPull1);
289 histolist->Add(hPrgPull2);
290 histolist->Add(hPrgPull3);
291 histolist->Add(hPrgPull4);
295 TVector3 vertexMC,vertexPoc,vertexFast,vertexFit,vertexKin,vertexDiff,vertexDiffFast,vertexDiffFit,vertexDiffKin;
297 const TVector3 nullpunkt(0.,0.,0.);
298 double rnd1=0.,rnd2=0.,rnd3=0.;
304 Float_t helixparams[5];
308 for(
int i=0;
i<nevt;
i++)
310 printf(
"Event %i of %i. (seed=%i)\n",
i,nevt,ranseed);
316 if (
fCands->GetEntriesFast()<2) {
317 cout<<
"Warning: too few particles created."<<endl;
323 if(
fCands->GetEntriesFast()==3)
326 combiCand = aCand->
Combine(bCand,cCand);
327 combiCand1 = aCand->
Combine(bCand,cCand);
328 combiCand2 = aCand->
Combine(bCand,cCand);
329 }
else if(
fCands->GetEntriesFast()==4)
333 combiCand = aCand->
Combine(bCand,cCand,dCand);
334 combiCand1 = aCand->
Combine(bCand,cCand,dCand);
335 combiCand2 = aCand->
Combine(bCand,cCand,dCand);
337 combiCand = aCand->
Combine(bCand);
338 combiCand1 = aCand->
Combine(bCand);
339 combiCand2 = aCand->
Combine(bCand);
342 if(
fVerbose>1) cout<<
"aCand, bCand, comboCand, aPos, bPos, combipos"<<endl;
345 if(
fVerbose>1) cout<<*combiCand<<endl;
347 if(
fVerbose>1) bCand->GetPosition().Print();
357 vertexDiff=mccand->
Pos();vertexDiff-=tmpcand->
Pos();
358 hVtxDiffX->Fill(vertexDiff.X());
359 hVtxDiffY->Fill(vertexDiff.Y());
360 hVtxDiffZ->Fill(vertexDiff.Z());
365 amomcov=tmpcand->
P4Cov();
366 hVtxDiffPX->Fill(momdiff.X());
367 hVtxDiffPY->Fill(momdiff.Y());
368 hVtxDiffPZ->Fill(momdiff.Z());
369 hVtxDiffE->Fill(momdiff.E());
371 hVtxPullPX->Fill(momdiff.X()/
sqrt(amomcov[0][0]));
372 hVtxPullPY->Fill(momdiff.Y()/
sqrt(amomcov[1][1]));
373 hVtxPullPZ->Fill(momdiff.Z()/
sqrt(amomcov[2][2]));
374 hVtxPullE->Fill(momdiff.E()/
sqrt(amomcov[3][3]));
376 hVtxDiffThe->Fill(mommc.Theta()-mom.Theta());
377 hVtxDiffPhi->Fill(mommc.Phi()-mom.Phi());
380 Bool_t prgcheckmc = PndAnalysisCalcTools::P7toPRG(mccand->Pos(),mccand->P4(),mccand->Charge(),mccand->Cov7(),nullpunkt, mcparams, helixCov, jacobian, kFALSE);
381 Bool_t prgchecktc = PndAnalysisCalcTools::P7toPRG(tmpcand->
Pos(),tmpcand->
P4(),tmpcand->
Charge(),tmpcand->
Cov7(),nullpunkt, helixparams, helixCov, jacobian, kFALSE);
384 hPrgPull0->Fill((mcparams[0]-helixparams[0])/
sqrt(helixCov[0][0]));
385 hPrgPull1->Fill((mcparams[1]-helixparams[1])/
sqrt(helixCov[1][1]));
386 hPrgPull2->Fill((mcparams[2]-helixparams[2])/
sqrt(helixCov[2][2]));
387 hPrgPull3->Fill((mcparams[3]-helixparams[3])/
sqrt(helixCov[3][3]));
388 hPrgPull4->Fill((mcparams[4]-helixparams[4])/
sqrt(helixCov[4][4]));
394 if(
fVerbose>0) cout<<
"poca "<<flush;
395 double dist = vPoca.
GetPocaVtx(vertexPoc,combiCand);
397 vertexDiff=vertexMC - vertexPoc;
402 hVtxPocaX->Fill(vertexDiff.X());
403 hVtxPocaY->Fill(vertexDiff.Y());
404 hVtxPocaZ->Fill(vertexDiff.Z());
405 hVtxPocaXY->Fill(vertexDiff.X(),vertexDiff.Y());
406 hVtxPocaRZ->Fill(vertexDiff.Z(),vertexMC.Perp()-vertexPoc.Perp());
407 hVtxPocas->Fill(dist);
408 hVtxPullPocaX->Fill(vertexDiff.X()/dist);
409 hVtxPullPocaY->Fill(vertexDiff.Y()/dist);
410 hVtxPullPocaZ->Fill(vertexDiff.Z()/dist);
417 vertexFast=nullpunkt;
419 double chiq=vFastter.
GetChi2();
420 int ndf = vFastter.
GetNdf();
421 double prob = vFastter.
GetProb();
425 hVtxChi2Fast->Fill(chiq);
426 hVtxChiProbFast->Fill(prob);
428 vertexDiffFast=vertexMC - vertexFast;
431 if(
fVerbose>0) vertexDiffFast.Print();
433 hVtxFastX->Fill(vertexDiffFast.X());
434 hVtxFastY->Fill(vertexDiffFast.Y());
435 hVtxFastZ->Fill(vertexDiffFast.Z());
436 hVtxFastXY->Fill(vertexDiffFast.X(),vertexDiffFast.Y());
437 hVtxFastRZ->Fill(vertexDiffFast.Z(),vertexMC.Perp()-vertexFast.Perp());
439 hVtxErrFastX->Fill(
sqrt(covV[0][0]));
440 hVtxErrFastY->Fill(
sqrt(covV[1][1]));
441 hVtxErrFastZ->Fill(
sqrt(covV[2][2]));
442 hVtxPullFastX->Fill(vertexDiffFast.X()/
sqrt(covV[0][0]));
443 hVtxPullFastY->Fill(vertexDiffFast.Y()/
sqrt(covV[1][1]));
444 hVtxPullFastZ->Fill(vertexDiffFast.Z()/
sqrt(covV[2][2]));
454 bool check = vFitter.
Fit();
455 double chiqfit = vFitter.
GetChi2();
456 int ndffit = vFitter.
GetNdf();
457 double probfit = vFitter.
GetProb();
459 if(
fVerbose>0) cout<<
" - chiq = "<<chiqfit<<
" \t"<<flush;
if(
fVerbose>1)cout<<endl;
463 hVtxChi2Fit->Fill(chiqfit);
464 hVtxChiProbFit->Fill(probfit);
466 vertexFit=combiCand1->GetFit()->DecayVtx();
467 covVfit=combiCand1->GetFit()->DecayVtx().CovMatrix();
468 vertexDiffFit=vertexMC - vertexFit;
471 if(
fVerbose>1) vertexDiffFit.Print();
473 hVtxFitX->Fill(vertexDiffFit.X());
474 hVtxFitY->Fill(vertexDiffFit.Y());
475 hVtxFitZ->Fill(vertexDiffFit.Z());
476 hVtxFitXY->Fill(vertexDiffFit.X(),vertexDiffFit.Y());
477 hVtxFitRZ->Fill(vertexDiffFit.Z(),vertexMC.Perp()-vertexFit.Perp());
479 hVtxErrFitX->Fill(
sqrt(covVfit[0][0]));
480 hVtxErrFitY->Fill(
sqrt(covVfit[1][1]));
481 hVtxErrFitZ->Fill(
sqrt(covVfit[2][2]));
482 hVtxPullFitX->Fill(vertexDiffFit.X()/
sqrt(covVfit[0][0]));
483 hVtxPullFitY->Fill(vertexDiffFit.Y()/
sqrt(covVfit[1][1]));
484 hVtxPullFitZ->Fill(vertexDiffFit.Z()/
sqrt(covVfit[2][2]));
486 for(
int k=0;k<combiCand1->NDaughters();k++)
490 tmpcand = tmpcand->
GetFit();
494 amomcov=tmpcand->
P4Cov();
495 hVtxDiffFitPX->Fill(momdiff.X());
496 hVtxDiffFitPY->Fill(momdiff.Y());
497 hVtxDiffFitPZ->Fill(momdiff.Z());
498 hVtxDiffFitE->Fill(momdiff.E());
500 hVtxPullFitPX->Fill(momdiff.X()/
sqrt(amomcov[0][0]));
501 hVtxPullFitPY->Fill(momdiff.Y()/
sqrt(amomcov[1][1]));
502 hVtxPullFitPZ->Fill(momdiff.Z()/
sqrt(amomcov[2][2]));
503 hVtxPullFitE->Fill(momdiff.E()/
sqrt(amomcov[3][3]));
505 hVtxDiffFitThe->Fill(mommc.Theta()-mom.Theta());
506 hVtxDiffFitPhi->Fill(mommc.Phi()-mom.Phi());
514 bool test = kFitter.
Fit();
515 double chiqkin = kFitter.
GetChi2();
516 int ndfkin = kFitter.
GetNdf();
517 double probkin = kFitter.
GetProb();
519 if(
fVerbose>0) cout<<
" - chiq = "<<chiqkin<<
" \t"<<flush;
if(
fVerbose>1)cout<<endl;
523 vertexKin=combiCand2->GetFit()->DecayVtx();
524 aposcov=combiCand2->GetFit()->DecayVtx().CovMatrix();
525 hVtxChi2Kin->Fill(chiqkin);
526 hVtxChiProbKin->Fill(probkin);
528 vertexDiffKin=vertexMC - vertexKin;
531 if(
fVerbose>0) vertexDiffKin.Print();
533 hVtxKinX->Fill(vertexDiffKin.X());
534 hVtxKinY->Fill(vertexDiffKin.Y());
535 hVtxKinZ->Fill(vertexDiffKin.Z());
536 hVtxKinXY->Fill(vertexDiffKin.X(),vertexDiffKin.Y());
537 hVtxKinRZ->Fill(vertexDiffKin.Z(),vertexMC.Perp()-vertexKin.Perp());
539 hVtxErrKinX->Fill(
sqrt(aposcov[0][0]));
540 hVtxErrKinY->Fill(
sqrt(aposcov[1][1]));
541 hVtxErrKinZ->Fill(
sqrt(aposcov[2][2]));
542 hVtxPullKinX->Fill(vertexDiffKin.X()/
sqrt(aposcov[0][0]));
543 hVtxPullKinY->Fill(vertexDiffKin.Y()/
sqrt(aposcov[1][1]));
544 hVtxPullKinZ->Fill(vertexDiffKin.Z()/
sqrt(aposcov[2][2]));
546 for(
int k=0;k<combiCand2->NDaughters();k++)
550 tmpcand = tmpcand->
GetFit();
554 amomcov=tmpcand->
P4Cov();
555 hVtxDiffKinPX->Fill(momdiff.X());
556 hVtxDiffKinPY->Fill(momdiff.Y());
557 hVtxDiffKinPZ->Fill(momdiff.Z());
558 hVtxDiffKinE->Fill(momdiff.E());
560 hVtxPullKinPX->Fill(momdiff.X()/
sqrt(amomcov[0][0]));
561 hVtxPullKinPY->Fill(momdiff.Y()/
sqrt(amomcov[1][1]));
562 hVtxPullKinPZ->Fill(momdiff.Z()/
sqrt(amomcov[2][2]));
563 hVtxPullKinE->Fill(momdiff.E()/
sqrt(amomcov[3][3]));
565 hVtxDiffKinThe->Fill(mommc.Theta()-mom.Theta());
566 hVtxDiffKinPhi->Fill(mommc.Phi()-mom.Phi());
571 cout<<
" ---------- calculations done, writing histograms."<<endl;
576 TString hfile=
"Data/HistoVertexing";
577 if(seed>=0){hfile+=
"-";hfile+=ranseed;}
579 TFile* histofile =
new TFile(hfile.Data(),
"RECREATE");
581 TListIter iter(histolist);
583 while( tmph=(TH1*)iter() ) tmph->Write();
598 double rnd1=widx*2*(0.5-gRandom->Rndm());
599 double rnd2=widy*2*(0.5-gRandom->Rndm());
600 double rnd3=widz*2*(0.5-gRandom->Rndm());
601 TVector3 vertex(rnd1,rnd2,rnd3);
602 if(
fVerbose>2) cout<<
"RollVertexBox: ("<<vertex.x()<<
","<<vertex.y()<<
","<<vertex.z()<<
")"<<endl;
608 double rnd1=gRandom->Gaus(vertex.x(),
fSigVx);
609 double rnd2=gRandom->Gaus(vertex.y(),
fSigVy);
610 double rnd3=gRandom->Gaus(vertex.z(),
fSigVz);
611 if(
fVerbose>2) cout<<
"SmearVertex (B): ("<<vertex.x()<<
","<<vertex.y()<<
","<<vertex.z()<<
")"<<endl;
612 vertex.SetXYZ(rnd1,rnd2,rnd3);
613 if(
fVerbose>2) cout<<
"SmearVertex (A): ("<<vertex.x()<<
","<<vertex.y()<<
","<<vertex.z()<<
")"<<endl;
622 double theta = vtx.Theta() + dtheta*2*(0.5-gRandom->Rndm());
624 double pt = ptmin + (ptmax-ptmin)*gRandom->Rndm();
626 momentum.SetPtThetaPhi(pt,theta,phi);
627 if(
fVerbose>2) cout<<
"RollMomentumBox: ("<<momentum.x()<<
","<<momentum.y()<<
","<<momentum.z()<<
")"<<endl;
633 double rnd1=gRandom->Gaus(momentum.Px(),
fSigPx);
634 double rnd2=gRandom->Gaus(momentum.Py(),
fSigPy);
635 double rnd3=gRandom->Gaus(momentum.Pz(),
fSigPz);
636 momentum.SetXYZ(rnd1,rnd2,rnd3);
642 if(
fVerbose>1) cout<<
" --------------------------------------------------- "<<endl;
643 if(
fVerbose>0) cout<<
"Poor man tracks."<<endl;
655 for (
int iTr=0; iTr<
fNumTrk; iTr++)
657 fPDG = TDatabasePDG::Instance()->GetParticle(
fPID);
660 int sizeMc =
fMcCands->GetEntriesFast();
675 int size =
fCands->GetEntriesFast();
676 TVector3 vertexMeasure =
fVertex;
678 TVector3 momMeasure = momentum;
700 if(
fVerbose>1) cout<<
" ---------\n pmc, p, pmcpos, ppos, fVertex"<<endl;
701 if(
fVerbose>1) std::cout<<*pmc<<std::endl;
702 if(
fVerbose>1) std::cout<<*p<<std::endl;
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
RhoCandidate * Combine(RhoCandidate *c)
SmearMomentum(TVector3 &momentum)
void SetPos(const TVector3 &pos)
friend F32vec4 sqrt(const F32vec4 &a)
void SmearVertex(TVector3 &vertex)
RhoCandidate * Daughter(Int_t n)
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
TVector3 RollVertexBox(double widx, double widy, double widz)
TVector3 fVertex(0., 0., 0.)
TString pt(TString pts, TString exts="px py pz")
TVector3 GetPosition() const
Double_t GetPocaVtx(TVector3 &vertex, RhoCandidate *composite)
TLorentzVector P4() const
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
RhoCandidate * GetMcTruth() const
double FitVertexFast(TVector3 &vtx, TMatrixD &cov, bool skipcov=false, int niter=2)
void SetMcTruth(RhoCandidate *mct)
RhoCandidate * GetFit() const
int poormantracks(int nevt=250, int laut=0, int seed=-1)
drchit SetVerbose(iVerbose)
TVector3 RollMomentumBox(const TVector3 &vtx, double dtheta, double ptmin, double ptmax)
void SetCov7(const TMatrixD &cov7)
TMatrixT< double > TMatrixD