1 #if !defined(__CINT__) || defined(__MAKECINT__)
4 #include <TDirectory.h>
8 #include <TGraphErrors.h>
9 #include "KMCDetector.h"
12 // Process the output file of testDetKMC.C macro to prepare summary
13 // for the "summary class" icl (definition of acceptable track).
14 // See ProcessSummary comment for details.
15 // The names of the graphs are prefixed by `pref`.
16 // If outF is provided, the graphs will be stored there
19 TObjArray* PrepSummaryKMC(const char* sumf, int cls, const char* pref=0, const char* outF=0);
20 TObjArray* ProcessSummary(TObjArray* sums, int icl, const char* pref=0);
22 enum {kSigAD,kSigAZ,kSigAPt,kSigD,kSigZ,kSigPt,kEff,kUpd};
25 TObjArray* PrepSummaryKMC(const char* sumf, int cls, const char* pref, const char* outF)
27 if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
28 TFile* fl = TFile::Open(sumf);
29 if (!fl) {printf("No file %s\n",sumf); return 0;}
30 TObjArray* arrs = (TObjArray*)gDirectory->Get("trSum");
31 if (!arrs) {printf("No summary in file %s\n",sumf); return 0;}
33 TObjArray* sums = ProcessSummary(arrs,cls,pref);
37 TFile* flOut = TFile::Open(outF,"update");
38 if (!flOut) {printf("Failed to open output file %s\n",outF);}
40 flOut->WriteObject(sums,sums->GetName(),"kSingleKey");
43 printf("Stored array %s in %s\n",sums->GetName(), outF);
51 TObjArray* ProcessSummary(TObjArray* arrs, int icl, const char* pref)
53 // Process TObjArray (e.g. for set of pt bins) of TObjArray of KMCTrackSummary objects:
54 // pick the KMCTrackSummary for "summary class" icl (definition of acceptable track) and create
56 // These graphs are returned in new TObjArray
59 if (!gs) gs = new TF1("gs","gaus",-1,1);
61 int nb = arrs->GetEntriesFast();
62 TObjArray* sums = (TObjArray*) arrs->At(0);
63 int nclass = sums->GetEntriesFast();
64 if (icl>=nclass) {printf("summary set has %d classes only, %d requested\n",nclass,icl);return 0;}
66 KMCTrackSummary* sm = (KMCTrackSummary*)sums->At(icl);
70 h = sm->GetHMCSigDCARPhi(); // MC resolution in transverse DCA
71 TGraphErrors * grSigD = 0;
73 grSigD = new TGraphErrors(nb);
74 grSigD->SetName(Form("%s%s",prefs.Data(),h->GetName()));
75 grSigD->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
78 TGraphErrors * grSigZ = 0;
79 h = sm->GetHMCSigDCAZ(); // MC resolution in Z DCA
81 grSigZ = new TGraphErrors(nb);
82 grSigZ->SetName(Form("%s%s",prefs.Data(),h->GetName()));
83 grSigZ->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
86 TGraphErrors * grSigAD = 0; // anaitical estimate for resolution in transverse DCA
88 grSigAD = new TGraphErrors(nb);
89 grSigAD->SetName(Form("%s%s",prefs.Data(),"sigmaDan"));
90 grSigAD->SetTitle(Form("%s%s",prefs.Data(),"#sigmaD an"));
93 TGraphErrors * grSigAZ = 0; // anaitical estimate for resolution in Z DCA
95 grSigAZ = new TGraphErrors(nb);
96 grSigAZ->SetName(Form("%s%s",prefs.Data(),"sigmaZan"));
97 grSigAZ->SetTitle(Form("%s%s",prefs.Data(),"#sigmaZ an"));
101 TGraphErrors * grSigPt = 0; // MC res. in pt
103 grSigPt = new TGraphErrors(nb);
104 grSigPt->SetName(Form("%s%s",prefs.Data(),"sigmaPt"));
105 grSigPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt"));
108 TGraphErrors * grSigAPt = 0; // analitycal res. in pt
110 grSigAPt = new TGraphErrors(nb);
111 grSigAPt->SetName(Form("%s%s",prefs.Data(),"sigmaPtan"));
112 grSigAPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt an"));
116 TGraphErrors * grEff = 0; // MC efficiency
118 grEff = new TGraphErrors(nb);
119 grEff->SetName(Form("%s_rate",prefs.Data()));
120 grEff->SetTitle(Form("%s Rate",prefs.Data()));
123 TGraphErrors * grUpd = 0; // number of Kalman track updates
125 grUpd = new TGraphErrors(nb);
126 grUpd->SetName(Form("%s_updCalls",prefs.Data()));
127 grUpd->SetTitle(Form("%s Updates",prefs.Data()));
130 for (int ib=0;ib<nb;ib++) {
131 sums = (TObjArray*) arrs->At(ib);
132 sm = (KMCTrackSummary*)sums->At(icl);
133 KMCProbe& prbRef = sm->GetRefProbe();
134 KMCProbe& prbAn = sm->GetAnProbe();
136 double pt = prbRef.Pt();
139 grSigAD->SetPoint(ib, pt,prbAn.GetSigmaY2()>0 ? TMath::Sqrt(prbAn.GetSigmaY2()) : 0.);
143 grSigAZ->SetPoint(ib, pt,prbAn.GetSigmaZ2()>0 ? TMath::Sqrt(prbAn.GetSigmaZ2()) : 0.);
147 double pts = TMath::Sqrt(prbAn.GetSigma1Pt2());
148 grSigAPt->SetPoint(ib, pt,pts>0 ? pts*pt : 0.);
152 h = sm->GetHMCSigPt();
154 grSigPt->SetPoint(ib, pt, gs->GetParameter(2));
155 grSigPt->SetPointError(ib, 0, gs->GetParError(2));
159 h = sm->GetHMCSigDCARPhi();
161 grSigD->SetPoint(ib, pt,gs->GetParameter(2));
162 grSigD->SetPointError(ib, 0,gs->GetParError(2));
166 h = sm->GetHMCSigDCAZ();
168 grSigZ->SetPoint(ib, pt,gs->GetParameter(2));
169 grSigZ->SetPointError(ib, 0,gs->GetParError(2));
173 grEff->SetPoint(ib, pt,sm->GetEff());
174 grEff->SetPointError(ib, 0,sm->GetEffErr());
178 grUpd->SetPoint(ib, pt,sm->GetUpdCalls());
179 grUpd->SetPointError(ib, 0, 0);
183 TObjArray* dest = new TObjArray();
184 dest->AddAtAndExpand(grSigAD,kSigAD);
185 dest->AddAtAndExpand(grSigAZ,kSigAZ);
186 dest->AddAtAndExpand(grSigAPt,kSigAPt);
187 dest->AddAtAndExpand(grSigD,kSigD);
188 dest->AddAtAndExpand(grSigZ,kSigZ);
189 dest->AddAtAndExpand(grSigPt,kSigPt);
190 dest->AddAtAndExpand(grEff,kEff);
191 dest->AddAtAndExpand(grUpd,kUpd);
193 if (!prefs.IsNull()) dest->SetName(pref);