source: pandaroot/trunk/PndTools/AnalysisTools/PndRhoTupleQA.cxx @ 30295

Revision 30135, 42.2 KB checked in by kgoetzen, 9 months ago (diff)

Updated PndRhoTupleQA/PndSimpleCombinerTask to store reco detector info. Added PndSimpleNtuple?.

Line 
1
2#include "PndRhoTupleQA.h"
3
4#include "RhoCandidate.h"
5#include "RhoCandList.h"
6#include "RhoTuple.h"
7#include "RhoFitterBase.h"
8#include "PndEventShape.h"
9#include "PndAnalysis.h"
10#include "PndPidCandidate.h"
11#include "PndVtxPoca.h"
12#include "PndVtxPRG.h"
13#include "TVector.h"
14#include <iostream>
15
16PndRhoTupleQA::PndRhoTupleQA(PndAnalysis *ana, double pbarmom)
17{
18  fAnalysis = ana;
19  fVtxPoca = new PndVtxPoca();
20  double mp=0.938272;
21  fIniP4.SetXYZT(0,0,pbarmom, sqrt(pbarmom*pbarmom+mp*mp)+mp);
22}
23// -------------------------------------------------------------------------
24
25PndRhoTupleQA::~PndRhoTupleQA()
26{
27  delete fVtxPoca;
28}
29// -------------------------------------------------------------------------
30
31void PndRhoTupleQA::qaESPidMult(TString pre, PndEventShape *evsh, double prob, double pmin, RhoTuple *n)
32{
33  n->Column(pre+"npide", (Int_t) evsh->MultElectronPminLab(prob,pmin),          0 );
34  n->Column(pre+"npidmu",(Int_t) evsh->MultMuonPminLab(prob,pmin),                      0 );
35  n->Column(pre+"npidpi",(Int_t) evsh->MultPionPminLab(prob,pmin),                      0 );
36  n->Column(pre+"npidk", (Int_t) evsh->MultKaonPminLab(prob,pmin),                      0 );
37  n->Column(pre+"npidp", (Int_t) evsh->MultProtonPminLab(prob,pmin),                    0 );
38}
39
40// -------------------------------------------------------------------------
41
42void PndRhoTupleQA::qaESMult(TString pre, PndEventShape *evsh, RhoTuple *n)
43{
44  n->Column(pre+"npart", (Int_t)        evsh->NParticles(),     0 );
45  n->Column(pre+"nneut", (Int_t)                evsh->NNeutral(),       0 );
46  n->Column(pre+"nchrg", (Int_t)        evsh->NCharged(),       0 );
47}
48
49// -------------------------------------------------------------------------
50
51void PndRhoTupleQA::qaESSum(TString pre, PndEventShape *evsh, RhoTuple *n)
52{
53  n->Column(pre+"sumpc",        (Float_t)  evsh->ChrgPSumCms(),         0.0f );
54  n->Column(pre+"sumpcl", (Float_t)  evsh->ChrgPSumLab(),               0.0f );
55  n->Column(pre+"sumen",        (Float_t)  evsh->NeutESumCms(),         0.0f );
56  n->Column(pre+"sumenl", (Float_t)  evsh->NeutESumLab(),               0.0f );
57
58  n->Column(pre+"sumpt",        (Float_t)  evsh->PtSumCms(),                    0.0f );
59  n->Column(pre+"sumptl",  (Float_t) evsh->PtSumLab(),                  0.0f );
60  n->Column(pre+"sumptc", (Float_t)  evsh->ChrgPtSumCms(),              0.0f );
61  n->Column(pre+"sumptcl",(Float_t)  evsh->ChrgPtSumLab(),              0.0f );
62  n->Column(pre+"sumetn",  (Float_t) evsh->NeutEtSumCms(),              0.0f );
63  n->Column(pre+"sumetnl", (Float_t) evsh->NeutEtSumLab(),              0.0f );
64}
65
66// -------------------------------------------------------------------------
67
68void PndRhoTupleQA::qaESMinMax(TString pre, PndEventShape *evsh, RhoTuple *n)
69{
70  n->Column(pre+"pmax", (Float_t)         evsh->PmaxCms(),      0.0f );
71  n->Column(pre+"ptmax",(Float_t)         evsh->Ptmax() ,       0.0f );
72  n->Column(pre+"pmaxl",(Float_t)         evsh->PmaxLab() ,     0.0f );
73  n->Column(pre+"pmin", (Float_t)         evsh->PminCms(),      0.0f );
74  n->Column(pre+"ptmin",(Float_t)   evsh->Ptmin(),      0.0f );
75  n->Column(pre+"pminl",(Float_t)         evsh->PminLab(),      0.0f );
76}
77
78// -------------------------------------------------------------------------
79
80void PndRhoTupleQA::qaESEventVars(TString pre, PndEventShape *evsh, RhoTuple *n)
81{
82  n->Column(pre+"prapmax",(Float_t) evsh->PRapmax(),            0.0f );
83  n->Column(pre+"sph", (Float_t)          evsh->Sphericity(),   0.0f );
84  n->Column(pre+"apl", (Float_t)          evsh->Aplanarity(),   0.0f );
85  n->Column(pre+"pla", (Float_t)          evsh->Planarity(),    0.0f );
86  n->Column(pre+"thr",(Float_t)   evsh->Thrust(),               0.0f );
87  n->Column(pre+"cir", (Float_t)          evsh->Circularity(),  0.0f );
88
89  n->Column(pre+"fw1", (Float_t)          evsh->FoxWolfMomR(1), 0.0f );
90  n->Column(pre+"fw2", (Float_t)          evsh->FoxWolfMomR(2), 0.0f );
91  n->Column(pre+"fw3", (Float_t)          evsh->FoxWolfMomR(3), 0.0f );
92  n->Column(pre+"fw4", (Float_t)          evsh->FoxWolfMomR(4), 0.0f );
93  n->Column(pre+"fw5", (Float_t)          evsh->FoxWolfMomR(5), 0.0f );
94}
95
96// -------------------------------------------------------------------------
97
98void PndRhoTupleQA::qaEventShape(TString pre, PndEventShape *evsh, RhoTuple *n)
99{
100  if (n==0) return;
101
102  // basic multiplicities
103  qaESMult(pre,  evsh, n);
104
105  // PID multiplicities
106  qaESPidMult(pre+"l",  evsh, 0.25, 0.0, n);
107  qaESPidMult(pre+"l1", evsh, 0.25, 1.0, n);
108  qaESPidMult(pre+"t",  evsh, 0.5,  0.0, n);
109  qaESPidMult(pre+"t1", evsh, 0.5,  1.0, n);
110  qaESPidMult(pre+"vt", evsh, 0.9,  0.0, n);
111
112  // event vars like thrust, sphericity etc
113  qaESEventVars(pre, evsh, n);
114
115  // standard sums over pt, p, E for all, neutral, charged in lab, cms
116  qaESSum(pre, evsh, n);
117
118  // standard min, max values of p, pt in lab, cms
119  qaESMinMax(pre, evsh, n);
120
121  // Multiplicities with min momemtum cut (cms)
122  n->Column(pre+"np05", (Int_t)   evsh->MultPminCms(0.5),       0 );
123  n->Column(pre+"np10", (Int_t)   evsh->MultPminCms(1.0),       0 );
124  n->Column(pre+"np20", (Int_t)   evsh->MultPminCms(2.0),       0 );
125  n->Column(pre+"np30", (Int_t)   evsh->MultPminCms(3.0),       0 );
126  n->Column(pre+"np40", (Int_t)   evsh->MultPminCms(4.0),       0 );
127  n->Column(pre+"np50", (Int_t)   evsh->MultPminCms(5.0),       0 );
128
129  // Multiplicities with min momemtum cut (lab)
130  n->Column(pre+"np05l", (Int_t)  evsh->MultPminLab(0.5),       0 );
131  n->Column(pre+"np10l", (Int_t)  evsh->MultPminLab(1.0),       0 );
132  n->Column(pre+"np20l", (Int_t)  evsh->MultPminLab(2.0),       0 );
133  n->Column(pre+"np30l", (Int_t)  evsh->MultPminLab(3.0),       0 );
134  n->Column(pre+"np40l", (Int_t)  evsh->MultPminLab(4.0),       0 );
135  n->Column(pre+"np50l",        (Int_t)  evsh->MultPminLab(5.0),        0 );
136
137  // Multiplicities with min p_t cut (cms)
138  n->Column(pre+"npt05",        (Int_t)  evsh->MultPtminCms(0.5),       0 );
139  n->Column(pre+"npt10",        (Int_t)  evsh->MultPtminCms(1.0),       0 );
140  n->Column(pre+"npt15",        (Int_t)  evsh->MultPtminCms(1.5),       0 );
141  n->Column(pre+"npt20",        (Int_t)  evsh->MultPtminCms(2.0),       0 );
142  n->Column(pre+"npt25",        (Int_t)  evsh->MultPtminCms(2.5),       0 );
143  n->Column(pre+"npt30",        (Int_t)  evsh->MultPtminCms(3.0),       0 );
144
145  // Neutral multiplicities with min energy cut (lab)
146  n->Column(pre+"nne003l",  (Int_t) evsh->MultNeutEminLab(0.03),        0 );
147  n->Column(pre+"nne005l",  (Int_t) evsh->MultNeutEminLab(0.05),        0 );
148  n->Column(pre+"nne01l",   (Int_t) evsh->MultNeutEminLab(0.1), 0 );
149  n->Column(pre+"nne05l",   (Int_t) evsh->MultNeutEminLab(0.5), 0 );
150
151  // Charged multiplicities with min momentum cut (lab)
152  n->Column(pre+"ncp005l", (Int_t)        evsh->MultChrgPminLab(0.05),  0 );
153  n->Column(pre+"ncp01l", (Int_t)         evsh->MultChrgPminLab(0.1),   0 );
154  n->Column(pre+"ncp02l", (Int_t)         evsh->MultChrgPminLab(0.2),   0 );
155  n->Column(pre+"ncp05l", (Int_t)         evsh->MultChrgPminLab(0.5),   0 );
156  n->Column(pre+"ncp10l", (Int_t)         evsh->MultChrgPminLab(1.0),   0 );
157
158  // Charged multiplicities with min momentum cut (cms)
159  n->Column(pre+"ncp005", (Int_t)         evsh->MultChrgPminCms(0.05),  0 );
160  n->Column(pre+"ncp01", (Int_t)          evsh->MultChrgPminCms(0.1),   0 );
161  n->Column(pre+"ncp02", (Int_t)          evsh->MultChrgPminCms(0.2),   0 );
162  n->Column(pre+"ncp05", (Int_t)          evsh->MultChrgPminCms(0.5),   0 );
163  n->Column(pre+"ncp10", (Int_t)          evsh->MultChrgPminCms(1.0),   0 );
164
165  // Sum of p_t  with min momentum cut (cms)
166  n->Column(pre+"sumpt05", (Float_t) evsh->SumPtminCms(0.5),    0.0f );
167  n->Column(pre+"sumpt10", (Float_t) evsh->SumPtminCms(1.0),    0.0f );
168
169  // Sum of charged momenta with min momentum cut (cms)
170  n->Column(pre+"sumpc05", (Float_t) evsh->SumChrgPminCms(0.5) ,0.0f );
171  n->Column(pre+"sumpc10", (Float_t) evsh->SumChrgPminCms(1.0),0.0f );
172
173  // Sum of charged momenta with min momentum cut (lab)
174  n->Column(pre+"sumpc05l", (Float_t)evsh->SumChrgPminLab(0.5) ,0.0f );
175  n->Column(pre+"sumpc10l",(Float_t) evsh->SumChrgPminLab(1.0) ,0.0f );
176
177  // Sum of neutral energy with min energy cut (cms)
178  n->Column(pre+"sumen05",  (Float_t)evsh->SumNeutEminCms(0.5) ,0.0f );
179  n->Column(pre+"sumen10", (Float_t) evsh->SumNeutEminCms(1.0) ,0.0f );
180
181  // Sum of neutral energy with min energy cut (lab)
182  n->Column(pre+"sumen05l",(Float_t) evsh->SumNeutEminLab(0.5) ,0.0f );
183  n->Column(pre+"sumen10l",(Float_t) evsh->SumNeutEminLab(1.0) ,0.0f );
184
185  // detector specific variables
186  n->Column(pre+"detemcsum",(Float_t) evsh->DetEmcSum()  ,0.0f );
187  n->Column(pre+"detemcmax",(Float_t) evsh->DetEmcMax()  ,0.0f );
188}
189// -------------------------------------------------------------------------
190void PndRhoTupleQA::qaEventShapeShort(TString pre, PndEventShape *evsh, RhoTuple *n)
191{
192  if (n==0) return;
193  // *** vars for PID multiplicity
194  //int ne, nmu, npi, nk, np;
195
196  // basic multiplicities
197  qaESMult(pre,  evsh, n);
198
199  // PID multiplicities
200  qaESPidMult(pre+"l",  evsh, 0.25, 0.0, n);
201  qaESPidMult(pre+"l1", evsh, 0.25, 1.0, n);
202
203  // event vars like thrust, sphericity etc
204  qaESEventVars(pre, evsh, n);
205
206  // standard sums over pt, p, E for all, neutral, charged in lab, cms
207  qaESSum(pre, evsh, n);
208
209  // standard min, max values of p, pt in lab, cms
210  qaESMinMax(pre, evsh, n);
211
212  // some multiplicities with min momentum cut (lab, cms)
213  n->Column(pre+"np10",  (Int_t)  evsh->MultPminCms(1.0),               0 );
214  n->Column(pre+"npt10", (Int_t)  evsh->MultPtminCms(1.0),      0 );
215  n->Column(pre+"ncp10l", (Int_t) evsh->MultChrgPminLab(1.0), 0 );
216  n->Column(pre+"nne10l",(Int_t)  evsh->MultNeutEminLab(1.0), 0 );
217
218  // sum of charged momenta with min momentum cut
219  n->Column(pre+"sumpc05", (Float_t) evsh->SumChrgPminCms(0.5) ,0.0f );
220}
221// -------------------------------------------------------------------------
222// *** store QA for PocaVtx
223void PndRhoTupleQA::qaPoca(TString pre, RhoCandidate *c, RhoTuple *n)
224{
225  if (n==0) return;
226
227  // *** simple vtx finder
228  TVector3 vtx, altvtx, primvtx;
229  double qavtx = fVtxPoca->GetPocaVtx(vtx, c);
230
231  // *** determine poca of rest of tracks
232  RhoCandList l;
233  fAnalysis->FillList(l, "Charged");
234  if (fVtxPoca->GetPocaVtx(primvtx, l)>998.) primvtx.SetXYZ(0.,0.,0.);
235
236  l.RemoveFamily(c);
237
238  if (l.GetLength()>1) fVtxPoca->GetPocaVtx(altvtx, l);
239  else altvtx = primvtx;
240
241  double dist=999.;
242  if (altvtx.Mag()>0.) dist = (vtx-altvtx).Mag();
243
244  // *** store QA info
245  n->Column(pre+"pocvx",  (Float_t) vtx.X(),   0.0f);
246  n->Column(pre+"pocvy",  (Float_t) vtx.Y(),   0.0f);
247  n->Column(pre+"pocvz",  (Float_t) vtx.Z(),   0.0f);
248  n->Column(pre+"altvx",  (Float_t) altvtx.X(),0.0f);
249  n->Column(pre+"altvy",  (Float_t) altvtx.Y(),0.0f);
250  n->Column(pre+"altvz",  (Float_t) altvtx.Z(),0.0f);
251  n->Column(pre+"pocmag", (Float_t) vtx.Mag(), 0.0f);
252  n->Column(pre+"pocqa",  (Float_t) qavtx,     0.0f);
253  n->Column(pre+"pocdist",(Float_t) dist,      0.0f);
254  n->Column(pre+"pocctau",(Float_t) (dist*c->M()/c->P()), 0.0f); // decay length
255}
256
257// -------------------------------------------------------------------------
258// *** store QA for PRG Vtx
259void PndRhoTupleQA::qaPRG(TString pre, RhoCandidate *c, RhoTuple *n)
260{
261  if (n==0) return;
262
263  PndVtxPRG vtxPRG(c);
264  vtxPRG.SetSilent();
265
266  // *** PRG vtx finder
267  TVector3 vtx(0,0,0);
268  TMatrixD vmatrix(3,3);
269  double qavtx = vtxPRG.FitVertexFast(vtx, vmatrix, true, 2);//iteration 2=default
270
271  // *** store QA info
272  n->Column(pre+"prgvx",  (Float_t) vtx.X(), 0.0f);
273  n->Column(pre+"prgvy",  (Float_t) vtx.Y(), 0.0f);
274  n->Column(pre+"prgvz",  (Float_t) vtx.Z(), 0.0f);
275  n->Column(pre+"prgqa",  (Float_t) qavtx,   0.0f);
276}
277
278// -------------------------------------------------------------------------
279// *** store QA for composite particles
280
281void PndRhoTupleQA::qaComp(TString pre, RhoCandidate *c, RhoTuple *n, bool covs, bool pulls)
282{
283  if (n==0) return;
284
285  // what kind of particle?
286  //int pdg = c->PdgCode(); //[R.K. 01/2017] unused variable
287
288  // special composite particle?
289  // pi0 or eta?
290  //if (pdg == 111 || pdg == 221)
291  //{
292    //qaPi0(pre, c, n);
293    //return;
294  //}
295  //// K_S?
296  //if (pdg == 310)
297  //{
298    //qaKs0(pre, c, n);
299    //return;
300  //}
301
302  // how many daughters?
303  int nd = c->NDaughters();
304
305  // truth match already done?
306  RhoCandidate *truth = c->GetMcTruth();
307
308  // if not, try one
309  if (truth == 0 && fAnalysis != 0)
310  {
311    fAnalysis->McTruthMatch(c);
312    truth = c->GetMcTruth();
313  }
314
315  bool mct = false;
316  if (truth!=0)
317  {
318    if (nd>0) mct = true;
319    else mct = (truth->PdgCode()==c->PdgCode());
320  }
321  // store cand info in lab and cms
322  qaCand(pre,   c,      n);
323  qaP4Cms(pre, c->P4(), n);
324  if(covs)qaP4Cov(pre,  c,      n);
325  if(pulls) qaPull(pre,c,n, (0==truth));
326  n->Column(pre+"mct",  (Float_t) mct, 0.0f);
327
328  // for mass difference e.g. D* -> D pi decays
329  if (nd>=2) n->Column(pre+"mdif", (Float_t) (c->M() - c->Daughter(0)->M()), 0.0f);
330
331  // cand is final state particle
332  if (nd==0)
333  {
334    // charged particle -> store PID info
335    if ( fabs(c->Charge())>0.1 ) qaPid(pre,     c,      n);
336
337    // and pdg code from truth
338    Float_t trpdg = 0.;
339    if (truth) trpdg = truth->PdgCode();
340    n->Column(pre+"trpdg", (Float_t) trpdg, 0.0f);
341  }
342  // cand is composite
343  else
344  {
345    // if 2 daughters -> store decay angle etc
346    if (nd==2) qa2Body(pre, c, n);
347    // if 3 daughters -> dalitz plot vars
348    if (nd==3) qaDalitz(pre, c, n);
349
350    // counter for charged final states
351    int nchrgfs = 0;
352
353    for (int i=0; i<nd; ++i)
354    {
355      RhoCandidate *dau = c->Daughter(i);
356      TString name=TString::Format("%sd%d",pre.Data(),i);
357
358      // count charged final states for possible vertexing later
359      if (dau->NDaughters()==0 && fabs(dau->Charge())>0.01) nchrgfs++;
360
361      // recursive call of qaComp
362      qaComp(name, dau, n, covs, pulls);
363    }
364    // only charged final state daughters -> Vtx info with PndVtxPoca
365    if (nchrgfs > 1)
366    {
367      qaPoca(pre, c, n);
368      qaVtx(pre, c, n);
369    }
370  }
371}
372
373// -------------------------------------------------------------------------
374
375void PndRhoTupleQA::qaPi0(TString pre, RhoCandidate *c, RhoTuple *n)
376{
377  if (n==0) return;
378
379  qaCand(pre, c, n);
380
381  RhoCandidate *d0 = c->Daughter(0);
382  RhoCandidate *d1 = c->Daughter(1);
383  double ang = d0->P3().Angle(d1->P3());
384
385  if (fAnalysis) fAnalysis->McTruthMatch(c);
386  RhoCandidate *truth = c->GetMcTruth();
387
388  qaCand(pre, c, n);
389  n->Column(pre+"oang",(Float_t) ang,           0.0f);
390
391  qaCand(pre+"d0", d0, n);
392  qaEmc(pre+"d0", d0, n);
393
394  qaCand(pre+"d1", d1, n);
395  qaEmc(pre+"d1", d1, n);
396
397  if (truth!=0)
398  {
399    qaCand("t"+pre, truth, n);
400    n->Column(pre+"mct", 1.0f, 0.0f);
401  }
402  else
403  {
404    qaCand("t"+pre, NULL , n, true);
405    n->Column(pre+"mct", 0.0f, 0.0f);
406  }
407}
408// -------------------------------------------------------------------------
409
410void PndRhoTupleQA::qaKs0(TString pre, RhoCandidate *c, RhoTuple *n)
411{
412  if (n==0) return;
413
414  RhoCandidate *d0 = c->Daughter(0);
415  RhoCandidate *d1 = c->Daughter(1);
416  double ang = d0->P3().Angle(d1->P3());
417
418  if (fAnalysis!=0) fAnalysis->McTruthMatch(c);
419  RhoCandidate *truth = c->GetMcTruth();
420
421  qaCand(pre, c, n);
422  qaPoca(pre, c, n);
423  qaVtx(pre, c, n);
424
425  n->Column(pre+"oang",(Float_t) ang,           0.0f);
426
427  qaCand(pre+"d0", d0, n);
428  qaPid(pre+"d0",  d0, n);
429
430  qaCand(pre+"d1", d1, n);
431  qaPid(pre+"d1",  d1, n);
432
433  if (truth!=0)
434  {
435    TVector3 vdist = truth->Pos() - truth->Daughter(0)->Pos();
436    Float_t dist = vdist.Mag();
437    Float_t ctau = dist * truth->M() / truth->P();
438
439    qaCand("t"+pre, truth, n);
440    n->Column(pre+"mct", 1.0f, 0.0f);
441    n->Column("t"+pre+"dist",    (Float_t) dist , 0.0f);
442    n->Column("t"+pre+"ctau", (Float_t) ctau , 0.0f);
443  }
444  else
445  {
446    qaCand("t"+pre, NULL , n, true);
447    n->Column(pre+"mct", 0.0f, 0.0f);
448    n->Column("t"+pre+"dist", (Float_t) -999. , 0.0f);
449    n->Column("t"+pre+"ctau", (Float_t) -999. , 0.0f);
450  }
451}
452
453// -------------------------------------------------------------------------
454
455void PndRhoTupleQA::qaP4(TString pre, TLorentzVector c, RhoTuple *n, bool skip)
456{
457  if (n==0) return;
458
459  if (!skip)
460  {
461    n->Column(pre+"px",  (Float_t) c.Px(),     0.0f );
462    n->Column(pre+"py",  (Float_t) c.Py(),     0.0f );
463    n->Column(pre+"pz",  (Float_t) c.Pz(),     0.0f );
464    n->Column(pre+"e",   (Float_t) c.E(),      0.0f );
465    n->Column(pre+"p",   (Float_t) c.P(),      0.0f );
466    n->Column(pre+"tht", (Float_t) c.Theta(),  0.0f );
467    n->Column(pre+"phi", (Float_t) c.Phi(),    0.0f );
468    n->Column(pre+"pt",  (Float_t) c.Pt(),     0.0f );
469    n->Column(pre+"m",   (Float_t) c.M(),      0.0f );
470  }
471  else
472  {
473    n->Column(pre+"px",  (Float_t) -999.,                       0.0f );
474    n->Column(pre+"py",  (Float_t) -999.,                       0.0f );
475    n->Column(pre+"pz",  (Float_t) -999.,                       0.0f );
476    n->Column(pre+"e",   (Float_t) -999.,                       0.0f );
477    n->Column(pre+"p",   (Float_t) -999.,                       0.0f );
478    n->Column(pre+"tht", (Float_t) -999.,                       0.0f );
479    n->Column(pre+"phi", (Float_t) -999.,                       0.0f );
480    n->Column(pre+"pt",  (Float_t) -999.,                       0.0f );
481    n->Column(pre+"m",   (Float_t) -999.,                       0.0f );
482  }
483}
484// -------------------------------------------------------------------------
485
486void PndRhoTupleQA::qaP4Cms(TString pre, TLorentzVector c, RhoTuple *n, bool skip)
487{
488  if (n==0) return;
489
490  if (!skip)
491  {
492    c.Boost(-fIniP4.BoostVector());
493
494    n->Column(pre+"pxcm",  (Float_t) c.Px(),     0.0f );
495    n->Column(pre+"pycm",  (Float_t) c.Py(),     0.0f );
496    n->Column(pre+"pzcm",  (Float_t) c.Pz(),     0.0f );
497    n->Column(pre+"ecm",   (Float_t) c.E(),      0.0f );
498    n->Column(pre+"pcm",   (Float_t) c.P(),      0.0f );
499    n->Column(pre+"thtcm", (Float_t) c.Theta(),  0.0f );
500    n->Column(pre+"phicm", (Float_t) c.Phi(),    0.0f );
501  }
502  else
503  {
504    n->Column(pre+"pxcm",  (Float_t) -999.,                     0.0f );
505    n->Column(pre+"pycm",  (Float_t) -999.,                     0.0f );
506    n->Column(pre+"pzcm",  (Float_t) -999.,                     0.0f );
507    n->Column(pre+"ecm",   (Float_t) -999.,                     0.0f );
508    n->Column(pre+"pcm",   (Float_t) -999.,                     0.0f );
509    n->Column(pre+"thtcm", (Float_t) -999.,                     0.0f );
510    n->Column(pre+"phicm", (Float_t) -999.,                     0.0f );
511  }
512}
513
514// -------------------------------------------------------------------------
515
516void PndRhoTupleQA::qaP4Cov(TString pre, RhoCandidate *c, RhoTuple *n, bool skip)
517{
518  if (n==0) return;
519
520  RhoError cov = c->P4Cov();
521
522  if (!skip)
523  {
524    n->Column(pre+"covpxpx", (Float_t) cov(0,0), 0.0f);
525    n->Column(pre+"covpxpy", (Float_t) cov(0,1), 0.0f);
526    n->Column(pre+"covpxpz", (Float_t) cov(0,2), 0.0f);
527    n->Column(pre+"covpxe",  (Float_t) cov(0,3), 0.0f);
528    n->Column(pre+"covpypy", (Float_t) cov(1,1), 0.0f);
529    n->Column(pre+"covpypz", (Float_t) cov(1,2), 0.0f);
530    n->Column(pre+"covpye",  (Float_t) cov(1,3), 0.0f);
531    n->Column(pre+"covpzpz", (Float_t) cov(2,2), 0.0f);
532    n->Column(pre+"covpze",  (Float_t) cov(2,3), 0.0f);
533    n->Column(pre+"covee",   (Float_t) cov(3,3), 0.0f);
534  }
535  else
536  {
537    n->Column(pre+"covpxpx", (Float_t) -999., 0.0f);
538    n->Column(pre+"covpxpy", (Float_t) -999., 0.0f);
539    n->Column(pre+"covpxpz", (Float_t) -999., 0.0f);
540    n->Column(pre+"covpxe",  (Float_t) -999., 0.0f);
541    n->Column(pre+"covpypy", (Float_t) -999., 0.0f);
542    n->Column(pre+"covpypz", (Float_t) -999., 0.0f);
543    n->Column(pre+"covpyee", (Float_t) -999., 0.0f);
544    n->Column(pre+"covpzpz", (Float_t) -999., 0.0f);
545    n->Column(pre+"covpze",  (Float_t) -999., 0.0f);
546    n->Column(pre+"covee",   (Float_t) -999., 0.0f);
547  }
548}
549
550// -------------------------------------------------------------------------
551
552void PndRhoTupleQA::qaCand(TString pre, RhoCandidate *cc, RhoTuple *n, bool skip)
553{
554  if (n==0) return;
555
556  if (!skip)
557  {
558    TLorentzVector c=cc->P4();
559    TVector3       p=cc->Pos();
560
561    qaP4(pre, c, n);
562    qaPos(pre, p, n);
563    n->Column(pre+"chg", (Float_t) cc->Charge(),        0.0f );
564    n->Column(pre+"pdg", (Float_t) cc->PdgCode(),       0.0f );
565  }
566  else
567  {
568    TLorentzVector dummy;
569    TVector3 dummypos;
570    qaP4(pre, dummy, n, true);
571    qaPos(pre, dummypos, n, true);
572    n->Column(pre+"chg", (Float_t) -999.,                       0.0f );
573    n->Column(pre+"pdg", (Float_t) -999.,                       0.0f );
574  }
575}
576// -------------------------------------------------------------------------
577
578void PndRhoTupleQA::qaPos(TString pre, TVector3 p, RhoTuple *n, bool skip)
579{
580  if (n==0) return;
581
582  if (!skip)
583  {
584    n->Column(pre+"x",   (Float_t) p.X(),       0.0f );
585    n->Column(pre+"y",   (Float_t) p.Y(),                       0.0f );
586    n->Column(pre+"z",   (Float_t) p.Z(),                       0.0f );
587    n->Column(pre+"l",   (Float_t) p.Mag(),             0.0f );
588  }
589  else
590  {
591    n->Column(pre+"x",   (Float_t) -999.,                       0.0f );
592    n->Column(pre+"y",   (Float_t) -999.,                       0.0f );
593    n->Column(pre+"z",   (Float_t) -999.,                       0.0f );
594    n->Column(pre+"l",   (Float_t) -999.,                       0.0f );
595  }
596}
597// -------------------------------------------------------------------------
598
599void PndRhoTupleQA::qaPull(TString pre, RhoCandidate *c, RhoTuple *n, bool skip)
600{
601  if (!skip)
602  {
603    RhoCandidate *mct = c->GetMcTruth();
604    if(mct) {
605      TLorentzVector difp4  = c->P4()  - mct->P4();
606      TVector3       difpos = c->Pos() - mct->Pos();
607      RhoError       covp4  = c->P4Cov();
608      RhoError       covpos = c->PosCov();
609      n->Column(pre+"pullpx",   (Float_t) (difp4.Px()/sqrt(covp4(0,0))),               0.0f );
610      n->Column(pre+"pullpy",   (Float_t) (difp4.Py()/sqrt(covp4(1,1))),               0.0f );
611      n->Column(pre+"pullpz",   (Float_t) (difp4.Pz()/sqrt(covp4(2,2))),               0.0f );
612      n->Column(pre+"pulle",    (Float_t) (difp4.E()/sqrt(covp4(3,3))),                0.0f );
613      n->Column(pre+"pullx",    (Float_t) (difpos.X()/sqrt(covpos(0,0))),              0.0f );
614      n->Column(pre+"pully",    (Float_t) (difpos.Y()/sqrt(covpos(1,1))),              0.0f );
615      n->Column(pre+"pullz",    (Float_t) (difpos.Z()/sqrt(covpos(2,2))),              0.0f );
616      n->Column(pre+"pullpxpy", (Float_t) (sqrt(difp4.Px()*difp4.Py()/covp4(0,1))),    0.0f );
617      n->Column(pre+"pullpxpz", (Float_t) (sqrt(difp4.Px()*difp4.Pz()/covp4(0,2))),    0.0f );
618      n->Column(pre+"pullpypz", (Float_t) (sqrt(difp4.Py()*difp4.Pz()/covp4(1,2))),    0.0f );
619      n->Column(pre+"pullpxe",  (Float_t) (sqrt(difp4.Px()*difp4.E()/covp4(0,3))),     0.0f );
620      n->Column(pre+"pullpye",  (Float_t) (sqrt(difp4.Py()*difp4.E()/covp4(1,3))),     0.0f );
621      n->Column(pre+"pullpze",  (Float_t) (sqrt(difp4.Pz()*difp4.E()/covp4(2,3))),     0.0f );
622      n->Column(pre+"pullxy",   (Float_t) (sqrt(difpos.X()*difpos.Y()/covpos(0,1))),   0.0f );
623      n->Column(pre+"pullxz",   (Float_t) (sqrt(difpos.Y()*difpos.Z()/covpos(0,2))),   0.0f );
624      n->Column(pre+"pullyz",   (Float_t) (sqrt(difpos.Y()*difpos.Z()/covpos(1,2))),   0.0f );
625    } else {skip=true;}
626  }
627 
628  if (skip)
629  {
630      n->Column(pre+"pullpx",   (Float_t) -999.,        0.0f );
631      n->Column(pre+"pullpy",   (Float_t) -999.,        0.0f );
632      n->Column(pre+"pullpz",   (Float_t) -999.,        0.0f );
633      n->Column(pre+"pulle",    (Float_t) -999.,        0.0f );
634      n->Column(pre+"pullx",    (Float_t) -999.,        0.0f );
635      n->Column(pre+"pully",    (Float_t) -999.,        0.0f );
636      n->Column(pre+"pullz",    (Float_t) -999.,        0.0f );
637      n->Column(pre+"pullpxpy", (Float_t) -999.,        0.0f );
638      n->Column(pre+"pullpxpz", (Float_t) -999.,        0.0f );
639      n->Column(pre+"pullpypz", (Float_t) -999.,        0.0f );
640      n->Column(pre+"pullpxe",  (Float_t) -999.,        0.0f );
641      n->Column(pre+"pullpye",  (Float_t) -999.,        0.0f );
642      n->Column(pre+"pullpze",  (Float_t) -999.,        0.0f );
643      n->Column(pre+"pullxy",   (Float_t) -999.,        0.0f );
644      n->Column(pre+"pullxz",   (Float_t) -999.,        0.0f );
645      n->Column(pre+"pullyz",   (Float_t) -999.,        0.0f );
646 }
647}
648// -------------------------------------------------------------------------
649
650void PndRhoTupleQA::qaTrk(TString pre, RhoCandidate *c, RhoTuple *n)
651{
652  if (n==0) return;
653
654  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
655
656  if (mic)
657  {
658    n->Column(pre+"trkdof",     (Int_t)   mic->GetDegreesOfFreedom(),   0 );
659    n->Column(pre+"trkstat",    (Int_t)   mic->GetFitStatus(),                  0 );
660    n->Column(pre+"trkchi2",    (Float_t) mic->GetChiSquared(),                 0.0f );
661    n->Column(pre+"trkidx",     (Int_t)   mic->GetTrackIndex(),                 0 );
662    n->Column(pre+"trkbranch",  (Int_t)   mic->GetTrackBranch() ,               0 );
663  }
664}
665
666// -------------------------------------------------------------------------
667
668void PndRhoTupleQA::qaPid(TString pre, RhoCandidate *c, RhoTuple *n)
669{
670  if (n==0) return;
671
672  // pinf[0]...pinf[4]: cache P values, pinf[5] = maximum
673  double pinf[6] = {-1000.,-1000.,-1000.,-1000.,-1000.,-1000.};
674  // index of particle type (e=0 ... p=4) with maximum P
675  int bestidx = -1;
676
677  if (fabs(c->Charge())>0.1)
678  {
679    for (int i=0;i<5;++i)
680    {
681      pinf[i] = c->GetPidInfo(i);
682
683      if (pinf[i]>pinf[5])
684      {
685        pinf[5] = pinf[i];
686        bestidx = i;
687      }
688    }
689  }
690
691  n->Column(pre+"pide",         (Float_t) pinf[0],              0.0f );
692  n->Column(pre+"pidmu",        (Float_t) pinf[1],              0.0f );
693  n->Column(pre+"pidpi",        (Float_t) pinf[2],              0.0f );
694  n->Column(pre+"pidk",         (Float_t) pinf[3],              0.0f );
695  n->Column(pre+"pidp",         (Float_t) pinf[4],              0.0f );
696
697  n->Column(pre+"pidmax",       (Float_t) pinf[5],              0.0f );
698  n->Column(pre+"pidbest",(Float_t) bestidx,            0.0f );
699
700}
701
702// -------------------------------------------------------------------------
703void PndRhoTupleQA::qa2Body(TString pre, RhoCandidate *c, RhoTuple *n)
704{
705  if (n==0) return;
706
707  if (c->NDaughters()!=2) return;
708
709  RhoCandidate *d0 = c->Daughter(0);
710  RhoCandidate *d1 = c->Daughter(1);
711
712  // opening angle lab
713  double oang = d0->P3().Angle(d1->P3());
714
715  // decay angle
716  TLorentzVector d_cms = d0->P4();
717  d_cms.Boost(-(c->P4().BoostVector()));
718  Float_t dec  = d_cms.Vect().Angle(c->P3());
719  Float_t cdec = cos(dec);
720
721  n->Column(pre+"oang",          (Float_t) oang,                0.0f );
722  n->Column(pre+"decang",  (Float_t) dec,               0.0f );
723  n->Column(pre+"cdecang", (Float_t) cdec,              0.0f );
724}
725
726
727// -------------------------------------------------------------------------
728void PndRhoTupleQA::qaDalitz(TString pre, RhoCandidate *c, RhoTuple *n)
729{
730  if (n==0) return;
731
732  if (c->NDaughters()!=3) return;
733
734  TLorentzVector l01 = c->Daughter(0)->P4()+c->Daughter(1)->P4();
735  TLorentzVector l12 = c->Daughter(1)->P4()+c->Daughter(2)->P4();
736  TLorentzVector l02 = c->Daughter(2)->P4()+c->Daughter(0)->P4();
737
738  n->Column(pre+"m01", (Float_t) l01.M(),               0.0f );
739  n->Column(pre+"m12", (Float_t) l12.M(),               0.0f );
740  n->Column(pre+"m02", (Float_t) l02.M(),               0.0f );
741
742  n->Column(pre+"dal01", (Float_t) l01.M2(),            0.0f );
743  n->Column(pre+"dal12", (Float_t) l12.M2(),            0.0f );
744  n->Column(pre+"dal02", (Float_t) l02.M2(),            0.0f );
745}
746
747
748// -------------------------------------------------------------------------
749
750void PndRhoTupleQA::qaVtx(TString pre, RhoCandidate *c, RhoTuple *n)
751{
752  if (n==0) return;
753
754        RhoCandidate *d = c->Daughter(0);
755  TVector3 v = c->DecayVtx();
756  if (d)
757  {
758    if (v.X()==0 && v.Y()==0 && v.Z()==0)
759      v = d->Pos();
760    //          TLorentzVector d_cms = d->P4();
761    //          d_cms.Boost(-(c->P4().BoostVector()));
762
763    Float_t ctau = v.Mag()*c->M()/c->P();
764    /*          Float_t dec  = d_cms.Vect().Angle(c->P3());
765     Float_t cdec = cos(dec);*/
766
767    // if primary Vertex available, compute ctau relative to that one
768    Float_t ctaud = -999.;
769    TVector3 primvtx;
770
771    // *** determine poca of all charged tracks as primary vertex
772    RhoCandList l;
773    fAnalysis->FillList(l, "Charged");
774    if (fVtxPoca->GetPocaVtx(primvtx, l)>998.) primvtx.SetXYZ(0.,0.,0.);
775
776    if (primvtx.Mag()>0) ctaud = (v-primvtx).Mag()*c->M()/c->P();
777
778    n->Column(pre+"vx",         (Float_t) v.X(),                0.0f );
779    n->Column(pre+"vy",         (Float_t) v.Y(),                0.0f );
780    n->Column(pre+"vz",         (Float_t) v.Z(),                0.0f );
781    n->Column(pre+"len",        (Float_t) v.Mag(),              0.0f );
782    n->Column(pre+"ctau",       (Float_t) ctau,                 0.0f );
783    n->Column(pre+"ctaud",      (Float_t) ctaud,                0.0f );
784    /*          n->Column(pre+"decang", (Float_t) dec,                  0.0f );
785     n->Column(pre+"cdecang",(Float_t) cdec,                    0.0f );*/
786  }
787  else
788  {
789    n->Column(pre+"vx",         (Float_t) -999.0,               0.0f );
790    n->Column(pre+"vy",         (Float_t) -999.0,               0.0f );
791    n->Column(pre+"vz",         (Float_t) -999.0,               0.0f );
792    n->Column(pre+"len",        (Float_t) -999.0,               0.0f );
793    n->Column(pre+"ctau",       (Float_t) -999.0,               0.0f );
794    n->Column(pre+"ctaud",      (Float_t) -999.0,               0.0f );
795    /*          n->Column(pre+"decang", (Float_t) -999.0,               0.0f );
796     n->Column(pre+"cdecang",(Float_t) -999.0,          0.0f );*/
797  }
798}
799
800// -------------------------------------------------------------------------
801
802void PndRhoTupleQA::qaEmc(TString pre, RhoCandidate *c, RhoTuple *n)
803{
804  if (n==0) return;
805
806  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
807
808  if (mic)
809  {
810    n->Column(pre+"emceraw",  (Float_t) mic->GetEmcRawEnergy(),                 0.0f );
811    n->Column(pre+"emcecal",  (Float_t) mic->GetEmcCalEnergy(),                 0.0f );
812    n->Column(pre+"emcqa",    (Float_t) mic->GetEmcQuality(),                   0.0f );
813    n->Column(pre+"emcnx",    (Int_t)   mic->GetEmcNumberOfCrystals(),  0 );
814    n->Column(pre+"emcnb",    (Int_t)   mic->GetEmcNumberOfBumps(),     0 );
815    n->Column(pre+"emcz20",   (Float_t) mic->GetEmcClusterZ20(),                0.0f );
816    n->Column(pre+"emcz53",   (Float_t) mic->GetEmcClusterZ53(),                0.0f );
817    n->Column(pre+"emclat",   (Float_t) mic->GetEmcClusterLat(),                0.0f );
818    n->Column(pre+"emce1",    (Float_t) mic->GetEmcClusterE1(),                 0.0f );
819    n->Column(pre+"emce9",    (Float_t) mic->GetEmcClusterE9(),                 0.0f );
820    n->Column(pre+"emce25",   (Float_t) mic->GetEmcClusterE25(),                0.0f );
821    n->Column(pre+"emcmod",   (Int_t)   mic->GetEmcModule(),                    0 );
822    n->Column(pre+"emcidx",   (Int_t)   mic->GetEmcIndex(),                     0 );
823  }
824}
825// -------------------------------------------------------------------------
826
827void PndRhoTupleQA::qaGem(TString pre, RhoCandidate *c, RhoTuple *n)
828{
829  if (n==0) return;
830
831  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
832
833  if (mic)
834  {
835    n->Column(pre+"gemnhits",  (Int_t) mic->GetGemHits(),                       0 );
836  }
837}
838
839
840// -------------------------------------------------------------------------
841
842void PndRhoTupleQA::qaMvd(TString pre, RhoCandidate *c, RhoTuple *n)
843{
844  if (n==0) return;
845
846  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
847
848  if (mic)
849  {
850    n->Column(pre+"mvddedx",  (Float_t) mic->GetMvdDEDX(),                      0.0f );
851    n->Column(pre+"mvdhits",  (Int_t)   mic->GetMvdHits(),                      0 );
852  }
853}
854
855// -------------------------------------------------------------------------
856
857void PndRhoTupleQA::qaStt(TString pre, RhoCandidate *c, RhoTuple *n)
858{
859  if (n==0) return;
860
861  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
862
863  if (mic)
864  {
865    n->Column(pre+"sttdedx",  (Float_t) mic->GetSttMeanDEDX(),          0.0f );
866    n->Column(pre+"stthits",  (Int_t)   mic->GetSttHits(),                      0 );
867  }
868}
869
870// -------------------------------------------------------------------------
871
872void PndRhoTupleQA::qaDrc(TString pre, RhoCandidate *c, RhoTuple *n)
873{
874  if (n==0) return;
875
876  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
877
878  if (mic)
879  {
880    n->Column(pre+"drcthtc",    (Float_t) mic->GetDrcThetaC(),                  0.0f );
881    n->Column(pre+"drcdthtc",   (Float_t) mic->GetDrcThetaCErr(),               0.0f );
882    n->Column(pre+"drcqa",      (Float_t) mic->GetDrcQuality(),                 0.0f );
883    n->Column(pre+"drcnphot",   (Int_t)   mic->GetDrcNumberOfPhotons(), 0 );
884    n->Column(pre+"drcidx",     (Int_t)   mic->GetDrcIndex(),                   0 );
885  }
886}
887
888// -------------------------------------------------------------------------
889
890void PndRhoTupleQA::qaDsc(TString pre, RhoCandidate *c, RhoTuple *n)
891{
892  if (n==0) return;
893
894  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
895
896  if (mic)
897  {
898    n->Column(pre+"dscthtc",    (Float_t) mic->GetDiscThetaC(),                 0.0f );
899    n->Column(pre+"dscdthtc",   (Float_t) mic->GetDiscThetaCErr(),              0.0f );
900    n->Column(pre+"dscqa",      (Float_t) mic->GetDiscQuality(),                0.0f );
901    n->Column(pre+"dscnphot",   (Int_t)   mic->GetDiscNumberOfPhotons(),0 );
902    n->Column(pre+"dscidx",     (Int_t)   mic->GetDiscIndex(),                  0 );
903  }
904}
905
906// -------------------------------------------------------------------------
907
908void PndRhoTupleQA::qaRich(TString pre, RhoCandidate *c, RhoTuple *n)
909{
910  if (n==0) return;
911
912  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
913
914  if (mic)
915  {
916    n->Column(pre+"richthtc",   (Float_t) mic->GetRichThetaC(),                 0.0f );
917    n->Column(pre+"richdthtc",  (Float_t) mic->GetRichThetaCErr(),              0.0f );
918    n->Column(pre+"richqa",     (Float_t) mic->GetRichQuality(),                0.0f );
919    n->Column(pre+"richnphot",  (Int_t)   mic->GetRichNumberOfPhotons(),0 );
920    n->Column(pre+"richidx",    (Int_t)    mic->GetRichIndex(),                 0 );
921  }
922}
923
924// -------------------------------------------------------------------------
925
926void PndRhoTupleQA::qaMuo(TString pre, RhoCandidate *c, RhoTuple *n)
927{
928  if (n==0) return;
929
930  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
931
932  if (mic)
933  {
934    n->Column(pre+"muonlay",    (Int_t)   mic->GetMuoNumberOfLayers(),  0 );;
935    n->Column(pre+"muoprob",    (Float_t) mic->GetMuoProbability(),             0.0f );;
936    n->Column(pre+"muoqa",      (Float_t) mic->GetMuoQuality() ,                0.0f ); ;
937    n->Column(pre+"muoiron",    (Float_t) mic->GetMuoIron() ,                   0.0f );  ;
938    n->Column(pre+"muopin",     (Float_t) mic->GetMuoMomentumIn(),              0.0f );;
939    n->Column(pre+"muomod",     (Int_t)   mic->GetMuoModule(),                  0 );  ;
940    n->Column(pre+"muohits",    (Int_t)   mic->GetMuoHits(),                    0 );  ;
941    n->Column(pre+"muoidx",     (Int_t)   mic->GetMuoIndex(),                   0 ); ;
942  }
943}
944
945// -------------------------------------------------------------------------
946
947void PndRhoTupleQA::qaTof(TString pre, RhoCandidate *c, RhoTuple *n)
948{
949  if (n==0) return;
950
951  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
952
953  if (mic)
954  {
955    n->Column(pre+"toftime",    (Float_t) mic->GetTofStopTime(),        0.0f );
956    n->Column(pre+"tofm2",      (Float_t) mic->GetTofM2(),                      0.0f );
957    n->Column(pre+"toflen",     (Float_t) mic->GetTofTrackLength(),     0.0f );
958    n->Column(pre+"tofqa",      (Float_t) mic->GetTofQuality(),         0.0f );
959    n->Column(pre+"tofidx",     (Int_t)   mic->GetTofIndex(),           0 );
960    n->Column(pre+"tofbeta",    (Float_t) mic->GetTofBeta(),            0.0f );
961  }
962}
963
964// -------------------------------------------------------------------------
965void PndRhoTupleQA::qaRecoShort(TString pre, RhoCandidate *c, RhoTuple *n)
966{
967  if (n==0) return;
968
969  PndPidCandidate *mic = (PndPidCandidate*)c->GetRecoCandidate();
970
971  if (mic)
972  {
973    n->Column(pre+"emcecal",  (Float_t) mic->GetEmcCalEnergy(),                 0.0f );
974    n->Column(pre+"emcnx",    (Int_t)   mic->GetEmcNumberOfCrystals(),  0 );
975    n->Column(pre+"emcnb",    (Int_t)   mic->GetEmcNumberOfBumps(),     0 );
976
977    n->Column(pre+"gemnhits",  (Int_t) mic->GetGemHits(),                           0 );
978
979    n->Column(pre+"mvddedx",  (Float_t) mic->GetMvdDEDX(),                          0.0f );
980    n->Column(pre+"mvdhits",  (Int_t)   mic->GetMvdHits(),                          0 );
981       
982    n->Column(pre+"sttdedx",  (Float_t) mic->GetSttMeanDEDX(),              0.0f );
983    n->Column(pre+"stthits",  (Int_t)   mic->GetSttHits(),                          0 );
984       
985    n->Column(pre+"drcthtc",    (Float_t) mic->GetDrcThetaC(),                  0.0f );
986    n->Column(pre+"drcnphot",   (Int_t)   mic->GetDrcNumberOfPhotons(), 0 );
987
988    n->Column(pre+"dscthtc",    (Float_t) mic->GetDiscThetaC(),                 0.0f );
989    n->Column(pre+"dscnphot",   (Int_t)   mic->GetDiscNumberOfPhotons(),0 );
990
991    n->Column(pre+"richthtc",   (Float_t) mic->GetRichThetaC(),                 0.0f );
992    n->Column(pre+"richnphot",  (Int_t)   mic->GetRichNumberOfPhotons(),0 );
993
994    n->Column(pre+"muonlay",    (Int_t)   mic->GetMuoNumberOfLayers(),  0 );;
995    n->Column(pre+"muoiron",    (Float_t) mic->GetMuoIron() ,                   0.0f );  ;
996
997    n->Column(pre+"tofm2",      (Float_t) mic->GetTofM2(),                          0.0f );
998    n->Column(pre+"tofbeta",    (Float_t) mic->GetTofBeta(),                0.0f );
999  }
1000}
1001
1002// -------------------------------------------------------------------------
1003void PndRhoTupleQA::qaRecoShortTree(TString pre, RhoCandidate *c, RhoTuple *n)
1004{
1005  if (n==0) return;
1006  if (c==0) return;
1007  int nd = c->NDaughters();
1008  if (nd>0)  {
1009    for (int i=0; i<nd; ++i)
1010    {
1011      RhoCandidate *dau = c->Daughter(i);
1012      TString name=TString::Format("%sd%d",pre.Data(),i);
1013      qaRecoShortTree(name,dau,n);
1014    }
1015  } else {
1016    qaRecoShort(pre,c,n);
1017  }
1018}
1019
1020
1021// -------------------------------------------------------------------------
1022void PndRhoTupleQA::qaRecoFull(TString pre, RhoCandidate *c, RhoTuple *n)
1023{
1024  if (n==0) return;
1025
1026  qaEmc( pre, c, n );
1027  qaMvd( pre, c, n );
1028  qaStt( pre, c, n );
1029  qaGem( pre, c, n );
1030  qaDrc( pre, c, n );
1031  qaDsc( pre, c, n );
1032  qaRich( pre, c, n );
1033  qaTof( pre, c, n );
1034  qaMuo( pre, c, n );
1035  qaTrk( pre, c, n);
1036}
1037
1038// -------------------------------------------------------------------------
1039void PndRhoTupleQA::qaRecoFullTree(TString pre, RhoCandidate *c, RhoTuple *n)
1040{
1041  if (n==0) return;
1042  if (c==0) return;
1043  int nd = c->NDaughters();
1044  if (nd>0)  {
1045    for (int i=0; i<nd; ++i)
1046    {
1047      RhoCandidate *dau = c->Daughter(i);
1048      TString name=TString::Format("%sd%d",pre.Data(),i);
1049      qaRecoFullTree(name,dau,n);
1050    }
1051  } else {
1052    qaRecoFull(pre,c,n);
1053  }
1054}
1055
1056// -------------------------------------------------------------------------
1057void PndRhoTupleQA::qaMc(TString pre, RhoCandidate *c, RhoTuple *n, bool skip)
1058{
1059  if (n==0) return;
1060  if (c==0) return;
1061
1062  RhoCandidate *mct=c->GetMcTruth();
1063  if (mct){
1064          qaCand(pre+"mc",mct,n,skip);
1065  }
1066}
1067
1068// -------------------------------------------------------------------------
1069void PndRhoTupleQA::qaMcList(RhoTuple *n, int max)
1070{
1071  RhoCandList mc;
1072  fAnalysis->FillList(mc,"McTruth");
1073  qaMcList("",mc,n,max);
1074}
1075
1076// -------------------------------------------------------------------------
1077void PndRhoTupleQA::qaMcList(TString pre, RhoCandList &l, RhoTuple *n, int max)
1078{
1079  if (n==0) return;
1080
1081  int npart = l.GetLength();
1082  if (npart>max) npart=max;
1083
1084  TVector vpart(npart), vpdg(npart), vmoth(npart),
1085  vp(npart),    vmass(npart), vndau(npart),
1086  vpx(npart),   vpy(npart),  vpz(npart), ve(npart),
1087  vtht(npart),  vphi(npart),
1088  vx(npart),    vy(npart),   vz(npart);
1089
1090  for (int j=0;j<npart;++j)
1091  {
1092    RhoCandidate *moth = l[j]->TheMother();
1093
1094    vpart(j) = j;
1095    vpdg(j)  = l[j]->PdgCode();
1096    vmoth(j) = (moth!=0x0) ? moth->GetTrackNumber() : -1;
1097    vndau(j) = l[j]->NDaughters();
1098
1099    vmass(j) = l[j]->Mass();
1100    vp(j)    = l[j]->P();
1101
1102    vpx(j)   = l[j]->Px();
1103    vpy(j)   = l[j]->Py();
1104    vpz(j)   = l[j]->Pz();
1105    ve(j)    = l[j]->E();
1106
1107    vtht(j)  = l[j]->P3().Theta();
1108    vphi(j)  = l[j]->P3().Phi();
1109
1110    vx(j)    = l[j]->Pos().X();
1111    vy(j)    = l[j]->Pos().Y();
1112    vz(j)    = l[j]->Pos().Z();
1113  }
1114
1115  n->Column(pre+"npart", (Int_t) npart);
1116
1117  n->Column(pre+"part",  vpart,  pre+"npart");
1118  n->Column(pre+"pdg",   vpdg,   pre+"npart");
1119  n->Column(pre+"moth",  vmoth,  pre+"npart");
1120  n->Column(pre+"ndau",  vndau,  pre+"npart");
1121
1122  n->Column(pre+"m",     vmass,  pre+"npart");
1123  n->Column(pre+"p",     vp,     pre+"npart");
1124
1125  n->Column(pre+"px",    vpx,    pre+"npart");
1126  n->Column(pre+"py",    vpy,    pre+"npart");
1127  n->Column(pre+"pz",    vpz,    pre+"npart");
1128  n->Column(pre+"e",     ve,     pre+"npart");
1129
1130  n->Column(pre+"tht",   vtht,   pre+"npart");
1131  n->Column(pre+"phi",   vphi,   pre+"npart");
1132
1133  n->Column(pre+"x",     vx,     pre+"npart");
1134  n->Column(pre+"y",     vy,     pre+"npart");
1135  n->Column(pre+"z",     vz,     pre+"npart");
1136}
1137
1138void PndRhoTupleQA::qaMcDiff(TString pre, RhoCandidate *c, RhoTuple *n, bool skip)
1139{
1140  if (n==0) return;
1141  if (c==0) return;
1142
1143  RhoCandidate *mct=c->GetMcTruth();
1144  if (!skip && mct)
1145  {
1146    TLorentzVector p4=c->P4();
1147    TLorentzVector mcp4=mct->P4();
1148    TLorentzVector diff = p4 - mcp4;
1149    TVector3 v=c->GetPosition();
1150    TVector3 mcv=mct->GetPosition();
1151    TVector3 vdiff= v - mcv;
1152    TMatrixD cov7 = c->Cov7();
1153
1154    n->Column(pre+"mcdiffvx",     (Float_t) vdiff.x() ,       0.0f );
1155    n->Column(pre+"mcdiffvy",     (Float_t) vdiff.y() ,       0.0f );
1156    n->Column(pre+"mcdiffvz",     (Float_t) vdiff.z() ,       0.0f );
1157    n->Column(pre+"mcdiffpx",     (Float_t) diff.Px(),        0.0f );
1158    n->Column(pre+"mcdiffpy",     (Float_t) diff.Py(),        0.0f );
1159    n->Column(pre+"mcdiffpz",     (Float_t) diff.Pz(),        0.0f );
1160    n->Column(pre+"mcdiffe",      (Float_t) diff.E(),         0.0f );
1161
1162    n->Column(pre+"mcdiffp",      (Float_t) (p4.P()-mcp4.P()),         0.0f );
1163    n->Column(pre+"mcdifftht",    (Float_t) (p4.Theta()-mcp4.Theta()), 0.0f );
1164    n->Column(pre+"mcdiffphi",    (Float_t) (p4.Phi()-mcp4.Phi()),     0.0f );
1165
1166    n->Column(pre+"mcpullvx",     (Float_t) ( vdiff.x()/sqrt(cov7(0,0)) ),       0.0f );
1167    n->Column(pre+"mcpullvy",     (Float_t) ( vdiff.y()/sqrt(cov7(1,1)) ),       0.0f );
1168    n->Column(pre+"mcpullvz",     (Float_t) ( vdiff.z()/sqrt(cov7(2,2)) ),       0.0f );
1169    n->Column(pre+"mcpullpx",     (Float_t) ( diff.Px()/sqrt(cov7(3,3)) ),       0.0f );
1170    n->Column(pre+"mcpullpy",     (Float_t) ( diff.Py()/sqrt(cov7(4,4)) ),       0.0f );
1171    n->Column(pre+"mcpullpz",     (Float_t) ( diff.Pz()/sqrt(cov7(5,5)) ),       0.0f );
1172    n->Column(pre+"mcpulle",      (Float_t) ( diff.E() /sqrt(cov7(6,6)) ),       0.0f );
1173
1174    n->Column(pre+"mcerrvx",     (Float_t) sqrt(cov7(0,0)),       0.0f );
1175    n->Column(pre+"mcerrvy",     (Float_t) sqrt(cov7(1,1)),       0.0f );
1176    n->Column(pre+"mcerrvz",     (Float_t) sqrt(cov7(2,2)),       0.0f );
1177    n->Column(pre+"mcerrpx",     (Float_t) sqrt(cov7(3,3)),       0.0f );
1178    n->Column(pre+"mcerrpy",     (Float_t) sqrt(cov7(4,4)),       0.0f );
1179    n->Column(pre+"mcerrpz",     (Float_t) sqrt(cov7(5,5)),       0.0f );
1180    n->Column(pre+"mcerre",      (Float_t) sqrt(cov7(6,6)),       0.0f );
1181  }
1182  else
1183  {
1184    TLorentzVector dummy;
1185    n->Column(pre+"mcdiffvx",     (Float_t) -999.0,       0.0f );
1186    n->Column(pre+"mcdiffvy",     (Float_t) -999.0,       0.0f );
1187    n->Column(pre+"mcdiffvz",     (Float_t) -999.0,       0.0f );
1188    n->Column(pre+"mcdiffpx",     (Float_t) -999.0,       0.0f );
1189    n->Column(pre+"mcdiffpy",     (Float_t) -999.0,       0.0f );
1190    n->Column(pre+"mcdiffpz",     (Float_t) -999.0,       0.0f );
1191    n->Column(pre+"mcdiffe",      (Float_t) -999.0,       0.0f );
1192    n->Column(pre+"mcdiffp",      (Float_t) -999.0,       0.0f );
1193    n->Column(pre+"mcdifftht",    (Float_t) -999.0,       0.0f );
1194    n->Column(pre+"mcdiffphi",    (Float_t) -999.0,       0.0f );
1195
1196    n->Column(pre+"mcpullvx",     (Float_t) -999.0,       0.0f );
1197    n->Column(pre+"mcpullvy",     (Float_t) -999.0,       0.0f );
1198    n->Column(pre+"mcpullvz",     (Float_t) -999.0,       0.0f );
1199    n->Column(pre+"mcpullpx",     (Float_t) -999.0,       0.0f );
1200    n->Column(pre+"mcpullpy",     (Float_t) -999.0,       0.0f );
1201    n->Column(pre+"mcpullpz",     (Float_t) -999.0,       0.0f );
1202    n->Column(pre+"mcpulle",      (Float_t) -999.0,       0.0f );
1203
1204    n->Column(pre+"mcerrvx",     (Float_t) -999.0,       0.0f );
1205    n->Column(pre+"mcerrvy",     (Float_t) -999.0,       0.0f );
1206    n->Column(pre+"mcerrvz",     (Float_t) -999.0,       0.0f );
1207    n->Column(pre+"mcerrpx",     (Float_t) -999.0,       0.0f );
1208    n->Column(pre+"mcerrpy",     (Float_t) -999.0,       0.0f );
1209    n->Column(pre+"mcerrpz",     (Float_t) -999.0,       0.0f );
1210    n->Column(pre+"mcerre",      (Float_t) -999.0,       0.0f );
1211  }
1212}
1213
1214void PndRhoTupleQA::qaFitter(TString pre, RhoFitterBase* fitter, RhoTuple* n, bool skip)
1215{
1216        if(!skip){
1217                n->Column(pre+"chisq" , (Float_t) fitter->GetChi2(),       0.0f);
1218                n->Column(pre+"ndf" ,   (Float_t) fitter->GetNdf(),        0.0f);
1219                n->Column(pre+"prob" ,  (Float_t) fitter->GetProb(),       0.0f);
1220        } else {
1221                n->Column(pre+"chisq" , (Float_t) -999.0,       0.0f);
1222                n->Column(pre+"ndf" ,   (Float_t) -999.0,       0.0f);
1223                n->Column(pre+"prob" ,  (Float_t) -999.0,       0.0f);
1224        }
1225}
1226
Note: See TracBrowser for help on using the repository browser.