]>
Commit | Line | Data |
---|---|---|
54fc064a | 1 | #include "Detector.h" |
2 | #include <TMath.h> | |
3 | #include <TMatrixD.h> | |
4 | #include <TGraph.h> | |
5 | #include <TAxis.h> | |
6 | #include <TFormula.h> | |
c40d3fcb | 7 | #include <TCanvas.h> |
8 | #include <TEllipse.h> | |
9 | #include <TText.h> | |
10 | #include <TGraphErrors.h> | |
54fc064a | 11 | |
12 | /*********************************************************** | |
13 | ||
14 | Fast Simulation tool for Inner Tracker Systems | |
15 | ||
16 | original code of using the billoir technique was developed | |
17 | for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov | |
18 | http://rnc.lbl.gov/~jhthomas | |
19 | ||
c40d3fcb | 20 | Changes by S. Rossegger -> see header file |
54fc064a | 21 | |
22 | ***********************************************************/ | |
23 | ||
24 | ||
c40d3fcb | 25 | #define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector |
54fc064a | 26 | |
27 | #define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 ) | |
28 | #define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm) | |
29 | #define dNdEtaMinB 950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700) | |
30 | #define dNdEtaCent 2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known) | |
31 | ||
c40d3fcb | 32 | #define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns) |
54fc064a | 33 | #define AcceptanceOfTpcAndSi 1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors |
34 | #define UPCBackgroundMultiplier 1.0 // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0) | |
35 | #define OtherBackground 0.0 // Increase multiplicity in detector (0.0 to 1.0 * minBias) (eg 0.0) | |
c40d3fcb | 36 | #define EfficiencySearchFlag 2 // Define search method: |
37 | // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0. | |
54fc064a | 38 | |
39 | #define PionMass 0.139 // Mass of the Pion | |
40 | #define KaonMass 0.498 // Mass of the Kaon | |
41 | #define D0Mass 1.865 // Mass of the D0 | |
42 | ||
43 | ||
44 | ||
45 | class CylLayer : public TNamed { | |
46 | public: | |
47 | ||
48 | CylLayer(char *name) : TNamed(name,name) {} | |
c40d3fcb | 49 | |
50 | Float_t GetRadius() const {return radius;} | |
51 | Float_t GetRadL() const {return radL;} | |
52 | Float_t GetPhiRes() const {return phiRes;} | |
53 | Float_t GetZRes() const {return zRes;} | |
54 | Float_t GetLayerEff() const {return eff;} | |
54fc064a | 55 | |
56 | // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); } | |
57 | Float_t radius; Float_t radL; Float_t phiRes; Float_t zRes; | |
c40d3fcb | 58 | Float_t eff; |
54fc064a | 59 | Bool_t isDead; |
60 | ||
61 | ClassDef(CylLayer,1); | |
62 | }; | |
63 | ||
64 | ClassImp(Detector) | |
65 | Detector::Detector() | |
66 | : TNamed("test_detector","detector"), | |
c40d3fcb | 67 | fNumberOfLayers(0), |
68 | fNumberOfActiveLayers(0), | |
54fc064a | 69 | fBField(0.5), |
70 | fLhcUPCscale(2.5), | |
71 | fIntegrationTime(0.02), // in ms | |
c40d3fcb | 72 | fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence |
54fc064a | 73 | fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle |
c40d3fcb | 74 | fParticleMass(0.140), // Standard: pion mass |
75 | fMaxRadiusSlowDet(10.) | |
54fc064a | 76 | { |
77 | // | |
78 | // default constructor | |
79 | // | |
80 | // fLayers = new TObjArray(); | |
c40d3fcb | 81 | |
54fc064a | 82 | } |
83 | ||
84 | Detector::Detector(char *name, char *title) | |
85 | : TNamed(name,title), | |
c40d3fcb | 86 | fNumberOfLayers(0), |
87 | fNumberOfActiveLayers(0), | |
54fc064a | 88 | fBField(0.5), |
89 | fLhcUPCscale(2.5), | |
90 | fIntegrationTime(0.02), // in ms | |
c40d3fcb | 91 | fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence |
54fc064a | 92 | fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle |
c40d3fcb | 93 | fParticleMass(0.140), // Standard: pion mass |
94 | fMaxRadiusSlowDet(10.) | |
54fc064a | 95 | { |
96 | // | |
97 | // default constructor, that set the name and title | |
98 | // | |
99 | // fLayers = new TObjArray(); | |
100 | } | |
101 | Detector::~Detector() { // | |
102 | // virtual destructor | |
103 | // | |
104 | // delete fLayers; | |
105 | } | |
106 | ||
c40d3fcb | 107 | void Detector::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes, Float_t zRes, Float_t eff) { |
108 | // | |
109 | // Add additional layer to the list of layers (ordered by radius) | |
110 | // | |
111 | ||
54fc064a | 112 | CylLayer *newLayer = (CylLayer*) fLayers.FindObject(name); |
113 | ||
114 | if (!newLayer) { | |
115 | newLayer = new CylLayer(name); | |
116 | newLayer->radius = radius; | |
117 | newLayer->radL = radL; | |
118 | newLayer->phiRes = phiRes; | |
119 | newLayer->zRes = zRes; | |
c40d3fcb | 120 | newLayer->eff = eff; |
54fc064a | 121 | |
122 | if (newLayer->zRes==RIDICULOUS && newLayer->zRes==RIDICULOUS) | |
123 | newLayer->isDead = kTRUE; | |
124 | else | |
125 | newLayer->isDead = kFALSE; | |
126 | ||
127 | if (fLayers.GetEntries()==0) | |
128 | fLayers.Add(newLayer); | |
129 | else { | |
130 | ||
131 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { | |
132 | CylLayer *l = (CylLayer*)fLayers.At(i); | |
133 | if (radius<l->radius) { | |
134 | fLayers.AddBefore(l,newLayer); | |
135 | break; | |
136 | } | |
137 | if (radius>l->radius && (i+1)==fLayers.GetEntries() ) { | |
138 | // even bigger then last one | |
139 | fLayers.Add(newLayer); | |
140 | } | |
141 | } | |
142 | ||
143 | } | |
c40d3fcb | 144 | fNumberOfLayers += 1; |
145 | if (!(newLayer->isDead)) fNumberOfActiveLayers += 1; | |
146 | ||
54fc064a | 147 | |
148 | } else { | |
149 | printf("Layer with the name %s does already exist\n",name); | |
150 | } | |
151 | ||
152 | ||
153 | } | |
154 | ||
155 | void Detector::KillLayer(char *name) { | |
156 | // | |
157 | // Marks layer as dead. Contribution only by Material Budget | |
158 | // | |
159 | ||
160 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
161 | if (!tmp) | |
162 | printf("Layer %s not found - cannot mark as dead\n",name); | |
163 | else { | |
c40d3fcb | 164 | tmp->phiRes = 999999; |
165 | tmp->zRes = 999999; | |
166 | if (!(tmp->isDead)) { | |
167 | tmp->isDead = kTRUE; | |
168 | fNumberOfActiveLayers -= 1; | |
169 | } | |
54fc064a | 170 | } |
171 | } | |
172 | ||
173 | void Detector::SetRadius(char *name, Float_t radius) { | |
174 | // | |
175 | // Set layer radius [cm] | |
176 | // | |
177 | ||
178 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
179 | ||
180 | ||
181 | if (!tmp) { | |
c40d3fcb | 182 | printf("Layer %s not found - cannot set radius\n",name); |
54fc064a | 183 | } else { |
184 | ||
185 | Float_t tmpRadL = tmp->radL; | |
186 | Float_t tmpPhiRes = tmp->phiRes; | |
187 | Float_t tmpZRes = tmp->zRes; | |
188 | ||
189 | RemoveLayer(name); // so that the ordering is correct | |
190 | AddLayer(name,radius,tmpRadL,tmpPhiRes,tmpZRes); | |
191 | } | |
192 | } | |
193 | ||
194 | Float_t Detector::GetRadius(char *name) { | |
195 | // | |
196 | // Return layer radius [cm] | |
197 | // | |
198 | ||
199 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
200 | if (!tmp) | |
c40d3fcb | 201 | printf("Layer %s not found - cannot get radius\n",name); |
54fc064a | 202 | else |
c40d3fcb | 203 | return tmp->radius; |
54fc064a | 204 | |
205 | return 0; | |
206 | } | |
207 | ||
208 | void Detector::SetRadiationLength(char *name, Float_t radL) { | |
209 | // | |
c40d3fcb | 210 | // Set layer material [cm] |
54fc064a | 211 | // |
212 | ||
213 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
214 | if (!tmp) | |
c40d3fcb | 215 | printf("Layer %s not found - cannot set layer material\n",name); |
54fc064a | 216 | else { |
217 | tmp->radL = radL; | |
218 | } | |
219 | } | |
220 | ||
221 | Float_t Detector::GetRadiationLength(char *name) { | |
222 | // | |
223 | // Return layer radius [cm] | |
224 | // | |
225 | ||
226 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
227 | if (!tmp) | |
c40d3fcb | 228 | printf("Layer %s not found - cannot get layer material\n",name); |
54fc064a | 229 | else |
230 | return tmp->radL; | |
231 | ||
232 | return 0; | |
233 | ||
234 | } | |
235 | ||
236 | void Detector::SetResolution(char *name, Float_t phiRes, Float_t zRes) { | |
237 | // | |
238 | // Set layer resolution in [cm] | |
239 | // | |
240 | ||
241 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
242 | if (!tmp) | |
243 | printf("Layer %s not found - cannot set resolution\n",name); | |
244 | else { | |
c40d3fcb | 245 | |
246 | Bool_t wasDead = tmp->isDead; | |
247 | ||
54fc064a | 248 | tmp->phiRes = phiRes; |
249 | tmp->zRes = zRes; | |
c40d3fcb | 250 | |
251 | if (zRes==RIDICULOUS && phiRes==RIDICULOUS) { | |
54fc064a | 252 | tmp->isDead = kTRUE; |
c40d3fcb | 253 | if (!wasDead) fNumberOfActiveLayers -= 1; |
254 | } else { | |
54fc064a | 255 | tmp->isDead = kFALSE; |
c40d3fcb | 256 | if (wasDead) fNumberOfActiveLayers += 1; |
257 | } | |
258 | ||
259 | ||
54fc064a | 260 | } |
261 | } | |
262 | ||
263 | Float_t Detector::GetResolution(char *name, Int_t axis) { | |
264 | // | |
265 | // Return layer resolution in [cm] | |
266 | // axis = 0: resolution in rphi | |
267 | // axis = 1: resolution in z | |
268 | // | |
269 | ||
270 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
271 | if (!tmp) | |
c40d3fcb | 272 | printf("Layer %s not found - cannot get resolution\n",name); |
54fc064a | 273 | else { |
274 | if (axis==0) return tmp->phiRes; | |
275 | if (axis==1) return tmp->zRes; | |
276 | printf("error: axis must be either 0 or 1 (rphi or z axis)\n"); | |
277 | } | |
278 | return 0; | |
279 | } | |
280 | ||
c40d3fcb | 281 | void Detector::SetLayerEfficiency(char *name, Float_t eff) { |
282 | // | |
283 | // Set layer efficnecy (prop that his is missed within this layer) | |
284 | // | |
285 | ||
286 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
287 | if (!tmp) | |
288 | printf("Layer %s not found - cannot set layer efficiency\n",name); | |
289 | else { | |
290 | tmp->eff = eff; | |
291 | } | |
292 | } | |
293 | ||
294 | Float_t Detector::GetLayerEfficiency(char *name) { | |
295 | // | |
296 | // Get layer efficnecy (prop that his is missed within this layer) | |
297 | // | |
298 | ||
299 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
300 | if (!tmp) | |
301 | printf("Layer %s not found - cannot get layer efficneicy\n",name); | |
302 | else | |
303 | return tmp->eff; | |
304 | ||
305 | return 0; | |
306 | ||
307 | } | |
308 | ||
54fc064a | 309 | void Detector::RemoveLayer(char *name) { |
310 | // | |
311 | // Removes a layer from the list | |
312 | // | |
313 | ||
314 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); | |
315 | if (!tmp) | |
316 | printf("Layer %s not found - cannot remove it\n",name); | |
317 | else { | |
c40d3fcb | 318 | Bool_t wasDead = tmp->isDead; |
54fc064a | 319 | fLayers.Remove(tmp); |
c40d3fcb | 320 | fNumberOfLayers -= 1; |
321 | if (!wasDead) fNumberOfActiveLayers -= 1; | |
54fc064a | 322 | } |
323 | } | |
324 | ||
325 | ||
326 | void Detector::PrintLayout() { | |
327 | // | |
328 | // Prints the detector layout | |
329 | // | |
330 | ||
331 | printf("Detector %s: \"%s\"\n",GetName(),GetTitle()); | |
332 | ||
333 | if (fLayers.GetEntries()>0) | |
334 | printf(" Name \t\t r [cm] \t X0 \t phi & z res [um]\n"); | |
335 | ||
336 | CylLayer *tmp = 0; | |
337 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { | |
338 | tmp = (CylLayer*)fLayers.At(i); | |
339 | ||
340 | // don't print all the tpc layers | |
341 | TString name(tmp->GetName()); | |
342 | if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue; | |
343 | ||
344 | printf("%d. %s \t %03.2f \t%1.4f\t ",i, | |
345 | tmp->GetName(), tmp->radius, tmp->radL); | |
346 | if (tmp->phiRes==RIDICULOUS) | |
347 | printf(" - "); | |
348 | else | |
349 | printf("%3.0f ",tmp->phiRes*10000); | |
350 | if (tmp->zRes==RIDICULOUS) | |
351 | printf(" -\n"); | |
352 | else | |
353 | printf("%3.0f\n",tmp->zRes*10000); | |
354 | } | |
355 | } | |
356 | ||
c40d3fcb | 357 | void Detector::PlotLayout(Int_t plotDead) { |
358 | // | |
359 | // Plots the detector layout in Front view | |
360 | // | |
361 | ||
362 | Double_t x0=0, y0=0; | |
363 | ||
364 | TGraphErrors *gr = new TGraphErrors(); | |
365 | gr->SetPoint(0,0,0); | |
366 | CylLayer *lastLayer = (CylLayer*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->radius; | |
367 | gr->SetPointError(0,maxRad,maxRad); | |
368 | gr->Draw("APE"); | |
369 | ||
370 | ||
371 | CylLayer *tmp = 0; | |
372 | for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) { | |
373 | tmp = (CylLayer*)fLayers.At(i); | |
374 | ||
375 | ||
376 | Double_t txtpos = tmp->radius; | |
377 | if ((tmp->isDead)) txtpos*=-1; // | |
378 | TText *txt = new TText(x0,txtpos,tmp->GetName()); | |
379 | txt->SetTextSizePixels(5); txt->SetTextAlign(21); | |
380 | if (!tmp->isDead || plotDead) txt->Draw(); | |
381 | ||
382 | TEllipse *layEl = new TEllipse(x0,y0,tmp->radius); | |
383 | // layEl->SetFillColor(5); | |
384 | layEl->SetFillStyle(5001); | |
385 | layEl->SetLineStyle(tmp->isDead+1); // dashed if not active | |
386 | layEl->SetLineColor(4); | |
387 | TString name(tmp->GetName()); | |
388 | if (!tmp->isDead) layEl->SetLineWidth(2); | |
389 | if (name.Contains("tpc") ) layEl->SetLineColor(29); | |
390 | ||
391 | if (!tmp->isDead || plotDead) layEl->Draw(); | |
392 | ||
393 | } | |
394 | ||
395 | } | |
396 | ||
397 | ||
398 | ||
54fc064a | 399 | void Detector::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) { |
400 | // | |
401 | // Emulates the TPC | |
402 | // | |
403 | // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow | |
c40d3fcb | 404 | |
405 | ||
406 | AddLayer((char*)"IFC", 77.8,0.01367); // Inner Field cage | |
54fc064a | 407 | |
408 | // % Radiation Lengths ... Average per TPC row (i.e. total/159 ) | |
409 | Float_t radLPerRow = 0.000036; | |
c40d3fcb | 410 | |
54fc064a | 411 | Float_t tpcInnerRadialPitch = 0.75 ; // cm |
412 | Float_t tpcMiddleRadialPitch = 1.0 ; // cm | |
413 | Float_t tpcOuterRadialPitch = 1.5 ; // cm | |
414 | // Float_t tpcInnerPadWidth = 0.4 ; // cm | |
415 | // Float_t tpcMiddlePadWidth = 0.6 ; // cm | |
416 | // Float_t tpcOuterPadWidth = 0.6 ; // cm | |
417 | Float_t innerRows = 63 ; | |
418 | Float_t middleRows = 64 ; | |
419 | Float_t outerRows = 32 ; | |
420 | Float_t tpcRows = (innerRows + middleRows + outerRows) ; | |
421 | Float_t rowOneRadius = 85.2 ; // cm | |
422 | Float_t row64Radius = 135.1 ; // cm | |
423 | Float_t row128Radius = 199.2 ; // cm | |
424 | ||
425 | for ( Int_t k = 0 ; k < tpcRows ; k++ ) { | |
426 | ||
427 | Float_t rowRadius =0; | |
428 | if (k<innerRows) | |
429 | rowRadius = rowOneRadius + k*tpcInnerRadialPitch ; | |
430 | else if ( k>=innerRows && k<(innerRows+middleRows) ) | |
431 | rowRadius = row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ; | |
432 | else if (k>=(innerRows+middleRows) && k<tpcRows ) | |
433 | rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ; | |
434 | ||
435 | if ( k%skip == 0 ) | |
436 | AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow,phiResMean,zResMean); | |
437 | else | |
438 | AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow); // non "active" row | |
439 | ||
440 | ||
441 | } | |
c40d3fcb | 442 | |
443 | } | |
444 | ||
445 | void Detector::RemoveTPC() { | |
54fc064a | 446 | |
447 | // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-) | |
448 | CylLayer *tmp = 0; | |
449 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { | |
450 | tmp = (CylLayer*)fLayers.At(i); | |
451 | TString name(tmp->GetName()); | |
c40d3fcb | 452 | if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; } |
54fc064a | 453 | } |
c40d3fcb | 454 | RemoveLayer((char*)"IFC"); |
455 | ||
54fc064a | 456 | } |
457 | ||
458 | ||
c40d3fcb | 459 | Double_t Detector::ThetaMCS ( Double_t mass, Double_t radLength, Double_t momentum ) const |
54fc064a | 460 | { |
c40d3fcb | 461 | // |
462 | // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14) | |
463 | // | |
464 | ||
54fc064a | 465 | Double_t beta = momentum / TMath::Sqrt(momentum*momentum+mass*mass) ; |
466 | Double_t theta = 0.0 ; // Momentum and mass in GeV | |
467 | // if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum ); | |
c40d3fcb | 468 | if ( radLength > 0 ) theta = 0.0136 * TMath::Sqrt(radLength) / ( beta * momentum ) * (1+0.038*TMath::Log(radLength)) ; |
54fc064a | 469 | return (theta) ; |
470 | } | |
471 | ||
472 | ||
c40d3fcb | 473 | Double_t Detector::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) |
54fc064a | 474 | { |
475 | // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm | |
476 | // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm | |
477 | // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius | |
c40d3fcb | 478 | Double_t sx, sy, goodHit ; |
479 | sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ; | |
480 | sy = 2 * TMath::Pi() * searchRadiusZ * searchRadiusZ * HitDensity(radius) ; | |
481 | goodHit = TMath::Sqrt(1./((1+sx)*(1+sy))) ; | |
482 | return ( goodHit ) ; | |
54fc064a | 483 | } |
484 | ||
485 | ||
c40d3fcb | 486 | Double_t Detector::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) |
54fc064a | 487 | { |
488 | // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm | |
489 | // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function | |
c40d3fcb | 490 | Double_t sx, goodHit ; |
491 | sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ; | |
492 | goodHit = 1./(1+sx) ; | |
493 | return ( goodHit ) ; | |
54fc064a | 494 | } |
495 | ||
c40d3fcb | 496 | Double_t Detector::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) |
497 | { | |
498 | // Based on work by Ruben Shahoyen | |
499 | // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function | |
500 | // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account | |
501 | // Following is correct for 2 DOF | |
502 | ||
503 | Double_t gamma = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level | |
504 | Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; | |
505 | Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*gamma)); | |
506 | ||
507 | return ( goodHit ) ; | |
508 | } | |
54fc064a | 509 | |
c40d3fcb | 510 | |
511 | Double_t Detector::HitDensity ( Double_t radius ) | |
54fc064a | 512 | { |
513 | // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor. | |
514 | // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results | |
515 | // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda. | |
516 | // Note that this function assumes we are working in CM and CM**2 [not meters]. | |
517 | // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2. | |
518 | ||
c40d3fcb | 519 | // Double_t MaxRadiusSlowDet = 0.1; //? // Maximum radius for slow detectors. Fast detectors |
54fc064a | 520 | // and only fast detectors reside outside this radius. |
c40d3fcb | 521 | Double_t arealDensity = 0 ; |
54fc064a | 522 | |
c40d3fcb | 523 | if ( radius > fMaxRadiusSlowDet ) |
54fc064a | 524 | { |
c40d3fcb | 525 | arealDensity = OneEventHitDensity(dNdEtaCent,radius) ; // Fast detectors see central collision density (only) |
526 | arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius) ; // Increase density due to background | |
54fc064a | 527 | } |
528 | ||
c40d3fcb | 529 | if (radius < fMaxRadiusSlowDet ) |
54fc064a | 530 | { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero. |
c40d3fcb | 531 | arealDensity = OneEventHitDensity(dNdEtaCent,radius) |
532 | + IntegratedHitDensity(dNdEtaMinB,radius) | |
533 | + UpcHitDensity(radius) ; | |
534 | arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ; | |
54fc064a | 535 | // Increase density due to background |
536 | } | |
537 | ||
c40d3fcb | 538 | return ( arealDensity ) ; |
54fc064a | 539 | } |
540 | ||
541 | ||
c40d3fcb | 542 | double Detector::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const |
54fc064a | 543 | { |
544 | // This is for one event at the vertex. No smearing. | |
c40d3fcb | 545 | double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ? |
54fc064a | 546 | // note: surface of sphere is '4*pi*r^2' |
547 | // surface of cylinder is '2*pi*r* h' | |
548 | return den ; | |
549 | } | |
550 | ||
551 | ||
c40d3fcb | 552 | double Detector::IntegratedHitDensity(Double_t multiplicity, Double_t radius) |
54fc064a | 553 | { |
554 | // The integral of minBias events smeared over a gaussian vertex distribution. | |
555 | // Based on work by Yan Lu 12/20/2006, all radii in centimeters. | |
556 | ||
c40d3fcb | 557 | Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ; |
558 | Double_t den = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ; | |
54fc064a | 559 | |
560 | // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex. | |
c40d3fcb | 561 | if ( den < OneEventHitDensity(multiplicity,radius) ) den = OneEventHitDensity(multiplicity,radius) ; |
54fc064a | 562 | |
563 | return den ; | |
564 | } | |
565 | ||
566 | ||
c40d3fcb | 567 | double Detector::UpcHitDensity(Double_t radius) |
54fc064a | 568 | { |
569 | // QED electrons ... | |
570 | ||
c40d3fcb | 571 | Double_t mUPCelectrons ; ; |
572 | mUPCelectrons = fLhcUPCscale * (1.23 - radius/6.5) ; // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC | |
573 | if ( mUPCelectrons < 0 ) mUPCelectrons = 0.0 ; // UPC electrons fall off quickly and don't go to large R | |
574 | mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ; // UPCs increase Mulitiplicty ~ proportional to MinBias rate | |
575 | mUPCelectrons *= UPCBackgroundMultiplier ; // Allow for an external multiplier (eg 0-1) to turn off UPC | |
54fc064a | 576 | |
c40d3fcb | 577 | return mUPCelectrons ; |
54fc064a | 578 | } |
579 | ||
580 | ||
581 | double Detector::Dist(double z, double r) | |
582 | { | |
583 | // Convolute dEta/dZ distribution with assumed Gaussian of vertex z distribution | |
584 | // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm | |
585 | // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters. | |
c40d3fcb | 586 | Int_t index = 1 ; // Start weight at 1 for Simpsons rule integration |
587 | Int_t nsteps = 301 ; // NSteps must be odd for Simpson's rule to work | |
54fc064a | 588 | double dist = 0.0 ; |
c40d3fcb | 589 | double dz0 = ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1) ; //cm |
54fc064a | 590 | double z0 = 0.0 ; //cm |
c40d3fcb | 591 | for(int i=0; i<nsteps; i++){ |
592 | if ( i == nsteps-1 ) index = 1 ; | |
54fc064a | 593 | z0 = -4*SigmaD + i*dz0 ; |
c40d3fcb | 594 | dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) * |
54fc064a | 595 | (1/sqrt((z-z0)*(z-z0) + r*r)) ; |
c40d3fcb | 596 | if ( index != 4 ) index = 4; else index = 2 ; |
54fc064a | 597 | } |
598 | return dist; | |
599 | } | |
600 | ||
601 | #define PZero 0.861 // Momentum of back to back decay particles in the CM frame | |
602 | #define EPiZero 0.872 // Energy of the pion from a D0 decay at rest | |
603 | #define EKZero 0.993 // Energy of the Kaon from a D0 decay at rest | |
604 | ||
c40d3fcb | 605 | Double_t Detector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const { |
54fc064a | 606 | // Math from Ron Longacre. Note hardwired energy to bin conversion for PtK and PtPi. |
607 | ||
c40d3fcb | 608 | Double_t const1 = pt / D0Mass ; |
609 | Double_t const2 = TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ; | |
610 | Double_t sum, ptPi, ptK ; | |
54fc064a | 611 | Double_t effp, effk ; |
612 | ||
613 | sum = 0.0 ; | |
614 | for ( Int_t k = 0 ; k < 360 ; k++ ) { | |
615 | ||
616 | Double_t theta = k * TMath::Pi() / 180. ; | |
617 | ||
c40d3fcb | 618 | ptPi = TMath::Sqrt( |
619 | PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 + | |
620 | const1*const1*EPiZero*EPiZero - | |
621 | 2*PZero*TMath::Cos(theta)*const2*const1*EPiZero + | |
54fc064a | 622 | PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta) |
623 | ) ; | |
624 | ||
c40d3fcb | 625 | ptK = TMath::Sqrt( |
626 | PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 + | |
627 | const1*const1*EKZero*EKZero + | |
628 | 2*PZero*TMath::Cos(theta)*const2*const1*EKZero + | |
54fc064a | 629 | PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta) |
630 | ) ; | |
631 | ||
632 | // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays | |
c40d3fcb | 633 | Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; |
634 | Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; | |
54fc064a | 635 | |
636 | if ( pionindex >= 400 ) pionindex = 399 ; | |
637 | if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ; | |
638 | if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd | |
639 | if ( effp < 0 ) effp = 0 ; | |
640 | ||
641 | if ( kaonindex >= 400 ) kaonindex = 399 ; | |
642 | if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ; | |
643 | if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd | |
644 | if ( effk < 0 ) effk = 0 ; | |
645 | ||
646 | // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here. | |
647 | ||
648 | sum += effp * effk ; | |
649 | ||
650 | } | |
651 | ||
652 | Double_t mean =sum/360; | |
c40d3fcb | 653 | return mean ; |
654 | ||
655 | } | |
54fc064a | 656 | |
657 | ||
c40d3fcb | 658 | void Detector::SolveDOFminusOneAverage() { |
659 | // | |
660 | // Short study to address "# layers-1 efficiencies" | |
661 | // saves the means in the according arrays | |
662 | // Note: Obviously, does not work for the Telescope equation | |
663 | // | |
664 | ||
665 | Double_t fMomentumResM[400], fResolutionRPhiM[400], fResolutionZM[400]; | |
666 | Double_t efficiencyM[3][400]; | |
667 | for (Int_t i=0; i<400; i++) { | |
668 | fMomentumResM[i] = 0; // Momentum resolution | |
669 | fResolutionRPhiM[i] = 0; // Resolution in R | |
670 | fResolutionZM[i] = 0; // Resolution in Z | |
671 | for (Int_t part=0; part<3; part++) | |
672 | efficiencyM[part][i] = 0; // efficiencies | |
673 | } | |
54fc064a | 674 | |
c40d3fcb | 675 | // loop over active layers in ITS (remove 1 by 1) |
676 | Int_t nITSLayers = 0; | |
677 | CylLayer *layer =0; | |
678 | for (Int_t j=0; j<(fLayers.GetEntries()-1); j++) { | |
679 | layer = (CylLayer*)fLayers.At(j); | |
680 | TString name(layer->GetName()); | |
681 | if (name.Contains("tpc")) continue; | |
682 | if (!(layer->isDead)) { | |
683 | ||
684 | nITSLayers++; | |
685 | printf("Kill Layer %s\n",name.Data()); | |
686 | Double_t rRes = GetResolution((char*)name.Data(),0); | |
687 | Double_t zRes = GetResolution((char*)name.Data(),1); | |
688 | KillLayer((char*)name.Data()); | |
689 | // PrintLayout(); | |
690 | SolveViaBilloir(1,0); | |
691 | ||
692 | // produce sum for the mean calculation | |
693 | for (Int_t i=0; i<400; i++) { | |
694 | fMomentumResM[i] += fMomentumRes[i]; // Momentum resolution | |
695 | fResolutionRPhiM[i] += fResolutionRPhi[i]; // Resolution in R | |
696 | fResolutionZM[i] += fResolutionZ[i]; // Resolution in Z | |
697 | for (Int_t part=0; part<3; part++) | |
698 | efficiencyM[part][i] += fEfficiency[part][i]; // efficiencies | |
699 | } | |
700 | ||
701 | // "Restore" layer ... | |
702 | SetResolution((char*)name.Data(),rRes,zRes); | |
703 | ||
704 | } | |
705 | } | |
54fc064a | 706 | |
c40d3fcb | 707 | // save means in "std. Arrays" |
708 | for (Int_t i=0; i<400; i++) { | |
709 | fMomentumRes[i] = fMomentumResM[i]/nITSLayers; // Momentum resolution | |
710 | fResolutionRPhi[i] = fResolutionRPhiM[i]/nITSLayers; // Resolution in R | |
711 | fResolutionZ[i] = fResolutionZM[i]/nITSLayers; // Resolution in Z | |
712 | for (Int_t part=0; part<3; part++) | |
713 | fEfficiency[part][i] = efficiencyM[part][i]/nITSLayers; // efficiencies | |
714 | } | |
715 | ||
716 | ||
54fc064a | 717 | } |
718 | ||
c40d3fcb | 719 | void Detector::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t meanPt) { |
720 | // | |
721 | // Solves the current geometry with the Billoir technique | |
722 | // ( see P. Billoir, Nucl. Instr. and Meth. 225 (1984), p. 352. ) | |
723 | // | |
54fc064a | 724 | |
c40d3fcb | 725 | Int_t nPt = 400; |
726 | // Clean up ...... | |
727 | for (Int_t i=0; i<kMaxNumberOfDetectors; i++) { | |
728 | for (Int_t j=0; j<nPt; j++) { | |
729 | fDetPointRes[i][j] = RIDICULOUS; | |
730 | fDetPointZRes[i][j] = RIDICULOUS; | |
731 | fTransMomenta[i] =0; | |
732 | fMomentumRes[i] =0; | |
733 | fResolutionRPhi[i] =0; | |
734 | } | |
735 | } | |
736 | ||
737 | if (!allPt) { // not the whole pt range -> allows a faster minimization at a defined 'meanpt' | |
738 | nPt = 3; | |
739 | } | |
54fc064a | 740 | |
741 | ||
54fc064a | 742 | |
743 | // Calculate track parameters using Billoirs method of matrices | |
744 | ||
c40d3fcb | 745 | Double_t pt, pz, lambda, momentum, rho, deltaPoverP ; |
746 | Double_t layerThickness, aMCS, mmm, charge ; | |
747 | Double_t vF1, vF2, vF3, vG1, vG2, vG3 ; | |
54fc064a | 748 | Double_t dx, sigma2, sigma2Z ; |
749 | Double_t mass[3] ; | |
c40d3fcb | 750 | Int_t printOnce = 1 ; |
54fc064a | 751 | |
752 | mass[0] = PionMass ; mass[1] = KaonMass ; // Loop twice for the D0; first pi then k | |
753 | ||
754 | mass[2] = fParticleMass; // third loop | |
755 | ||
756 | Int_t mStart =0; | |
757 | if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode | |
758 | ||
759 | for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ ) { | |
760 | ||
761 | // PseudoRapidity OK, used as an angle | |
762 | lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)) ; | |
763 | ||
764 | // Track moves along the x axis and deviations are in the y and z directions. | |
765 | Double_t y, z, a, b, d, e ; | |
766 | // a = py/px = tan phi, b = pz/px = tan lambda / cos phi, d = 0.3*Ze/Abs(p) | |
767 | ||
c40d3fcb | 768 | for ( Int_t i = 0 ; i < nPt ; i++ ) { // pt loop |
54fc064a | 769 | |
770 | CylLayer *last = (CylLayer*) fLayers.At((fLayers.GetEntries()-1)); | |
771 | ||
772 | // Convert to Meters, Tesla, and GeV | |
c40d3fcb | 773 | |
774 | // Starting values based on radius of outermost layer ... log10 steps to ~20 GeV | |
775 | Double_t bigRad = last->radius/100 /2 ; // min. pt which the algorithm below could handle | |
776 | if (bigRad<0.61) bigRad=0.61; // -> min pt around 100 MeV for Bz=0.5T (don't overdo it ... ;-) ) | |
777 | fTransMomenta[i] = ( 0.3*bigRad*TMath::Abs(fBField) ) - 0.08 + TMath::Power(10,2.3*i/nPt) / 10.0 ; | |
778 | if (!allPt) { // just 3 points around meanPt | |
779 | fTransMomenta[i] = meanPt-0.1+i*0.1; | |
780 | } | |
54fc064a | 781 | |
782 | ||
783 | // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis | |
784 | // These are the EndPoint values for y, z, a, b, and d | |
c40d3fcb | 785 | pt = fTransMomenta[i] ; // GeV/c |
54fc064a | 786 | rho = pt / TMath::Abs( 0.3 * fBField ); // Radius of curvature of the track (meters) |
787 | momentum = pt / TMath::Cos(lambda) ; // Total momentum | |
c40d3fcb | 788 | charge = -1 ; // Assume an electron |
54fc064a | 789 | pz = pt * TMath::Tan(lambda) ; // GeV/c |
c40d3fcb | 790 | d = 0.3 * charge / momentum ; // Its an electrons so q = -1, which makes d a negative number |
54fc064a | 791 | |
c40d3fcb | 792 | // printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d); |
54fc064a | 793 | |
794 | // Matrices for Billoir | |
c40d3fcb | 795 | TMatrixD bigA(2,2); TMatrixD littleA(2,2); |
796 | TMatrixD mIstar(5,5); TMatrixD mD(5,5); TMatrixD mDT(5,5); TMatrixD mA(5,5); TMatrixD mM(5,5); TMatrixD mWork(5,5) ; | |
797 | TMatrixD mCovStart(5,5); | |
798 | bigA.Zero();littleA.Zero();mD.Zero();mDT.Zero();mA.Zero();mM.Zero();mWork.Zero(); | |
799 | ||
54fc064a | 800 | // Set Detector-Efficiency Storage area to unity |
c40d3fcb | 801 | fEfficiency[massloop][i] = 1.0 ; |
54fc064a | 802 | // Back-propagate the covariance matrix along the track. I is the inverse of the covariance matrix. |
803 | // Start with small variations to a straight line 'information' matrix. Tuning this matrix is | |
804 | // equivalent to 'starting the recursion' in Billoirs paper. It is a tricky business. You need test | |
805 | // data and many layers for the matrices to stabilize. Do not believe results for a small number of layers. | |
806 | // In our case, always include the TPC and this matrix will be well conditioned before it gets to the Si. | |
807 | ||
54fc064a | 808 | |
c40d3fcb | 809 | mCovStart.UnitMatrix(); // start with unity (-> sigma estimate of 1cm) |
810 | ||
811 | // track parametrization: (y,z,a,b,1/p) with a=py/px=tan(phi) and b = pz/px = tan(lambda)/cos(phi) | |
812 | mCovStart(0,0) = (2*2); // 2cm variance in y | |
813 | mCovStart(1,1) = (2*2); // 2cm variance in z | |
814 | mCovStart(2,2) = (0.1*0.1); // approx 0.1rad in angles phi and lambda | |
815 | mCovStart(3,3) = (1*1); // approx 0.1rad in angle phi and lambda | |
816 | mCovStart(4,4) = (1*1); // 1 GeV error in pt | |
817 | ||
818 | ||
819 | mIstar = mCovStart; // start with above covariance Matrix | |
820 | mIstar.Invert(); | |
821 | ||
822 | // mIstar.UnitMatrix(); // start with unity (-> sigma estimate of ~ 1cm) | |
823 | ||
824 | // 'A' is the "angle" matrix for aMCS at each step. | |
54fc064a | 825 | // 'D' is the "distance" matrix at each step. It propagates the particle backwards along the track. |
826 | // 'M' is the "measurement" matrix. It represents the measurement resolution at each step. | |
827 | ||
828 | CylLayer *layer = 0; | |
829 | CylLayer *nextlayer = 0; | |
c40d3fcb | 830 | |
831 | // find last "active layer" - start tracking at the last active layer | |
832 | Int_t lastActiveLayer = 0; | |
833 | for (Int_t j=(fLayers.GetEntries()-1); j>0; j--) { | |
834 | layer = (CylLayer*)fLayers.At(j); | |
835 | if (!(layer->isDead)) { // is alive | |
836 | lastActiveLayer = j; | |
837 | break; | |
838 | } | |
839 | } | |
840 | ||
841 | for (Int_t j=lastActiveLayer; j>0; j--) { // Layer loop | |
54fc064a | 842 | |
843 | layer = (CylLayer*)fLayers.At(j); | |
844 | nextlayer = (CylLayer*)fLayers.At(j-1); | |
845 | ||
54fc064a | 846 | // Convert to Meters, Tesla, and GeV |
847 | Float_t radius = layer->radius/100; | |
54fc064a | 848 | Float_t radLength = layer->radL; |
849 | ||
850 | Float_t nextRadius = nextlayer->radius /100; | |
851 | Float_t nextPhiResolution = nextlayer->phiRes /100; | |
852 | Float_t nextZResolution = nextlayer->zRes /100; | |
c40d3fcb | 853 | |
854 | // if ( radius >= rho*TMath::Cos(TMath::Pi()/4) ) // (meters) protect against too low pt | |
855 | // { // printf("pt lower bound is too low for this algorithm to take this layer into account (r=%lf)\n",radius); | |
856 | // continue ; } | |
54fc064a | 857 | |
c40d3fcb | 858 | // Jims approximation of the "ideal position", breaks down when rho~radius because of the hyperbola approximation |
859 | // y = rho - TMath::Sqrt(rho*rho-radius*radius) ; // These are 'ideal' locations and | |
860 | // a = radius / ( rho - y ) ; // not propagated values which should | |
861 | // b = rho * TMath::Tan(lambda) / ( rho - y ) ; // be done if we had data. But we don't. | |
862 | // z = rho * TMath::Tan(lambda) * TMath::ATan(a) ; | |
863 | ||
864 | // Helix format -> exact intersection at the layer radii | |
865 | y=0;z=0; // not used | |
866 | a = radius*radius/TMath::Sqrt(4*radius*radius*rho*rho - TMath::Power(radius,4)); // = Tan phi = y_c/x_c (in cartesian coordinates) | |
867 | b = TMath::Tan(lambda)/TMath::Cos(TMath::ATan(a)); | |
868 | ||
54fc064a | 869 | e = TMath::Sqrt( 1 + a*a + b*b ) ; |
870 | ||
c40d3fcb | 871 | layerThickness = radLength / TMath::Sin(TMath::Pi()/2 - lambda) ; // X0 in lambda direction |
54fc064a | 872 | |
c40d3fcb | 873 | aMCS = ThetaMCS(mass[massloop], layerThickness, momentum) ; |
54fc064a | 874 | mmm = ( 1 + a*a + b*b ) ; |
c40d3fcb | 875 | bigA(0,0) = mmm*(1+a*a) ; bigA(0,1) = mmm*a*b ; littleA(0,0) = mIstar(2,2) ; littleA(0,1) = mIstar(2,3) ; |
876 | bigA(1,0) = mmm*a*b ; bigA(1,1) = mmm*(1+b*b) ; littleA(1,0) = mIstar(3,2) ; littleA(1,1) = mIstar(3,3) ; | |
877 | littleA = bigA.Invert() + aMCS*aMCS*littleA ; | |
878 | littleA = littleA.Invert() ; | |
879 | mA(0,0) = 0.0 ; mA(0,1) = 0.0 ; mA(0,2) = 0.0 ; mA(0,3) = 0.0 ; mA(0,4) = 0.0 ; | |
880 | mA(1,0) = 0.0 ; mA(1,1) = 0.0 ; mA(1,2) = 0.0 ; mA(1,3) = 0.0 ; mA(1,4) = 0.0 ; | |
881 | mA(2,0) = 0.0 ; mA(2,1) = 0.0 ; mA(2,2) = littleA(0,0) ; mA(2,3) = littleA(0,1) ; mA(2,4) = 0.0 ; | |
882 | mA(3,0) = 0.0 ; mA(3,1) = 0.0 ; mA(3,2) = littleA(1,0) ; mA(3,3) = littleA(1,1) ; mA(3,4) = 0.0 ; | |
883 | mA(4,0) = 0.0 ; mA(4,1) = 0.0 ; mA(4,2) = 0.0 ; mA(4,3) = 0.0 ; mA(4,4) = 0.0 ; | |
884 | mIstar = mIstar - aMCS*aMCS*mIstar*mA*mIstar ; | |
54fc064a | 885 | dx = ( radius - nextRadius ) ; // Billoir works with absolute magnitude of distance |
c40d3fcb | 886 | vF3 = e * ( -1 * ( 1 + a*a ) * fBField ) ; // Assume fBField is on Z axis, only. |
887 | vF1 = d * ( a*vF3/(e*e) - 2*e*a*fBField ) ; | |
888 | vF2 = d * ( b*vF3/(e*e) ) ; | |
889 | vG3 = e * ( -1 * a * b * fBField ) ; | |
890 | vG1 = d * ( a*vG3/(e*e) + e*b*fBField ) ; | |
891 | vG2 = d * ( b*vG3/(e*e) - e*a*fBField ) ; | |
892 | mD(0,0) = 1.0 ; mD(0,1) = 0.0 ; mD(0,2) = -1*dx+vF1*dx*dx/2 ; mD(0,3) = vF2*dx*dx/2 ; mD(0,4) = vF3*dx*dx/2 ; | |
893 | mD(1,0) = 0.0 ; mD(1,1) = 1.0 ; mD(1,2) = vG1*dx*dx/2 ; mD(1,3) = -1*dx+vG2*dx*dx/2 ; mD(1,4) = vG3*dx*dx/2 ; | |
894 | mD(2,0) = 0.0 ; mD(2,1) = 0.0 ; mD(2,2) = 1-vF1*dx ; mD(2,3) = -1*vF2*dx ; mD(2,4) = -1*vF3*dx ; | |
895 | mD(3,0) = 0.0 ; mD(3,1) = 0.0 ; mD(3,2) = -1*vG1*dx ; mD(3,3) = 1-vG2*dx ; mD(3,4) = -1*vG3*dx ; | |
896 | mD(4,0) = 0.0 ; mD(4,1) = 0.0 ; mD(4,2) = 0.0 ; mD(4,3) = 0.0 ; mD(4,4) = 1.0 ; | |
897 | mDT.Transpose(mD) ; | |
898 | mIstar = mDT.Invert() * mIstar * mD.Invert() ; | |
54fc064a | 899 | // Prepare to save Detector efficiencies |
c40d3fcb | 900 | mWork = mIstar ; // Working copy of matrix |
901 | mWork.Invert() ; // Invert the Matrix to recover the convariance matrix | |
902 | fDetPointRes [j-1][i] = TMath::Sqrt( mWork(0,0) ) ; // result in meters | |
903 | fDetPointZRes[j-1][i] = TMath::Sqrt( mWork(1,1) ) ; // result in meters | |
54fc064a | 904 | // End save |
c40d3fcb | 905 | |
54fc064a | 906 | sigma2 = ( nextPhiResolution * nextPhiResolution ) ; |
907 | sigma2Z = ( nextZResolution * nextZResolution ) ; | |
c40d3fcb | 908 | mM(0,0) = 1/sigma2; mM(1,1) = 1/sigma2Z ; // setting detector resolution matrix |
909 | mIstar = mIstar + mM ; | |
54fc064a | 910 | |
911 | } | |
912 | ||
c40d3fcb | 913 | // Pattern recognition is done .... save values like vertex resolution etc. |
914 | ||
54fc064a | 915 | // Invert the Matrix to recover the convariance matrix |
c40d3fcb | 916 | mIstar.Invert() ; |
54fc064a | 917 | // Convert the Convariance matrix parameters into physical quantities |
918 | // The results are propogated to the previous point but *do not* include the measurement at that point. | |
c40d3fcb | 919 | deltaPoverP = TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3 ; // Absolute magnitude so ignore charge |
920 | fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ) ; // results in percent | |
921 | fResolutionRPhi[i] = TMath::Sqrt( mIstar(0,0) ) * 1.e6 ; // result in microns | |
922 | fResolutionZ[i] = TMath::Sqrt( mIstar(1,1) ) * 1.e6 ; // result in microns | |
923 | // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius | |
924 | ||
925 | ||
926 | if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { | |
927 | printf("Number of active layers: %d\n",fNumberOfActiveLayers) ; | |
928 | printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000); | |
929 | printf("Name Radius Thickness PointResOn PointResOnZ DetRes DetResZ Density Efficiency\n") ; | |
930 | // printOnce =0; | |
54fc064a | 931 | } |
932 | ||
c40d3fcb | 933 | // print out and efficiency calculation |
934 | for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop | |
54fc064a | 935 | |
c40d3fcb | 936 | layer = (CylLayer*)fLayers.At(j); |
54fc064a | 937 | |
938 | // Convert to Meters, Tesla, and GeV | |
939 | Float_t radius = layer->radius /100; | |
940 | Float_t phiRes = layer->phiRes /100; | |
941 | Float_t zRes = layer->zRes /100; | |
942 | Float_t radLength = layer->radL; | |
c40d3fcb | 943 | Float_t leff = layer->eff; // basic layer efficiency |
54fc064a | 944 | Bool_t isDead = layer->isDead; |
c40d3fcb | 945 | |
54fc064a | 946 | |
c40d3fcb | 947 | if ( (!isDead && radLength >0) ) { |
948 | Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] + | |
54fc064a | 949 | phiRes * phiRes ) * 100. ; // work in cm |
c40d3fcb | 950 | Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] + |
54fc064a | 951 | zRes * zRes ) * 100. ; // work in cm |
c40d3fcb | 952 | |
953 | Double_t layerEfficiency = 0; | |
954 | if ( EfficiencySearchFlag == 0 ) | |
955 | layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ; | |
956 | else if ( EfficiencySearchFlag == 1 ) | |
957 | layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ; | |
958 | else if ( EfficiencySearchFlag == 2 ) | |
959 | layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; | |
54fc064a | 960 | |
c40d3fcb | 961 | TString name(layer->GetName()); |
962 | if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue; | |
963 | ||
964 | if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { | |
965 | ||
966 | ||
967 | printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ", | |
968 | layer->GetName(), radius*100, radLength, | |
969 | fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6, | |
54fc064a | 970 | phiRes*1.e6, zRes*1.e6, |
c40d3fcb | 971 | HitDensity(radius*100)) ; |
972 | if (!name.Contains("tpc")) | |
973 | printf("%10.3f\n", layerEfficiency); | |
974 | else | |
975 | printf(" - \n"); | |
54fc064a | 976 | } |
977 | ||
c40d3fcb | 978 | if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency; |
979 | ||
980 | } | |
981 | /* | |
982 | // vertex print | |
983 | if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) { | |
984 | printf("%s:\t ----- ----- %10.0f %11.0f \n", layer->GetName(),fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6); | |
54fc064a | 985 | } |
c40d3fcb | 986 | */ |
54fc064a | 987 | } |
c40d3fcb | 988 | if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { |
989 | if (fNumberOfActiveLayers >=15) printOnce = 0 ; | |
54fc064a | 990 | printf("\n") ; |
991 | } | |
992 | ||
c40d3fcb | 993 | |
994 | ||
995 | ||
996 | if (fNumberOfActiveLayers <15 ) { | |
997 | ||
998 | ||
999 | ||
1000 | // BACKWORD TRACKING +++++++++++++++++ | |
1001 | // number of layers is quite low ... efficiency calculation was probably nonsense | |
1002 | // Tracking outward (backword) to get reliable efficiencies from "smoothed estimates" | |
1003 | ||
1004 | // For below, see paper, NIM A262 (1987) p.444, eqs.12. | |
1005 | // Equivalently, one can simply combine the forward and backward estimates. Assuming | |
1006 | // pf,Cf and pb,Cb as extrapolated position estimates and errors from fwd and bwd passes one can | |
1007 | // use a weighted estimate Cw = (Cf^-1 + Cb^-1)^-1, pw = Cw (pf Cf^-1 + pb Cb^-1). | |
1008 | // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything. | |
1009 | ||
1010 | Bool_t doLikeAliRoot = kFALSE; // don't do the "combined info" but do like in Aliroot | |
1011 | ||
1012 | if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { | |
1013 | printf("- Numbers of active layer is low (%d):\n -> \"outward\" fitting done as well to get reliable eff.estimates\n", | |
1014 | fNumberOfActiveLayers); | |
1015 | } | |
1016 | ||
1017 | // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam) | |
1018 | // mIstar.UnitMatrix(); // start with unity | |
1019 | if (doLikeAliRoot) { | |
1020 | Double_t factor = 10; | |
1021 | TMatrixD diagElements(5,1); | |
1022 | for (Int_t ip =0; ip<5; ip++) diagElements(ip,0) = mIstar(ip,ip)*factor; | |
1023 | mIstar.Zero(); | |
1024 | for (Int_t ip =0; ip<5; ip++) mIstar(ip,ip) = diagElements(ip,0); | |
1025 | } else { | |
1026 | // Complete RESET of the covariance matrix | |
1027 | mIstar = mCovStart; | |
1028 | // starts with above covariance Matrix and uses "combinded resolutions" (all infos), see below | |
1029 | // this is the "honest" approach without biasing the results | |
1030 | } | |
1031 | mIstar.Invert() ; // undo the revert from above | |
1032 | ||
1033 | // Clean up and storing of "forward estimates" | |
1034 | Double_t detPointResForw[kMaxNumberOfDetectors][400], detPointZResForw[kMaxNumberOfDetectors][400] ; | |
1035 | for (Int_t k=0; k<kMaxNumberOfDetectors; k++) { | |
1036 | for (Int_t l=0; l<nPt; l++) { | |
1037 | detPointResForw[k][l] = fDetPointRes[k][l]; | |
1038 | if (!doLikeAliRoot) fDetPointRes[k][l] = RIDICULOUS; | |
1039 | detPointZResForw[k][l] = fDetPointZRes[k][l]; | |
1040 | if (!doLikeAliRoot) fDetPointZRes[k][l] = RIDICULOUS; | |
1041 | } | |
1042 | } | |
1043 | ||
1044 | // find first "active layer" - start tracking at the first active layer | |
1045 | Int_t firstActiveLayer = 0; | |
1046 | for (Int_t j=0; j<(fLayers.GetEntries()-1); j++) { | |
1047 | layer = (CylLayer*)fLayers.At(j); | |
1048 | if (!(layer->isDead)) { // is alive | |
1049 | firstActiveLayer = j; | |
1050 | break; | |
1051 | } | |
1052 | } | |
1053 | ||
1054 | for (Int_t j=firstActiveLayer; j<(fLayers.GetEntries()-1); j++) { // Layer loop | |
1055 | ||
1056 | layer = (CylLayer*)fLayers.At(j); | |
1057 | nextlayer = (CylLayer*)fLayers.At(j+1); | |
1058 | ||
1059 | // Convert to Meters, Tesla, and GeV | |
1060 | Float_t radius = layer->radius/100; | |
1061 | Float_t radLength = layer->radL; | |
1062 | ||
1063 | Float_t nextRadius = nextlayer->radius /100; | |
1064 | Float_t nextPhiResolution = nextlayer->phiRes /100; | |
1065 | Float_t nextZResolution = nextlayer->zRes /100; | |
1066 | ||
1067 | y = rho - TMath::Sqrt(rho*rho-radius*radius) ; // These are 'ideal' locations and | |
1068 | a = radius / ( rho - y ) ; // not propagated values which should | |
1069 | b = rho * TMath::Tan(lambda) / ( rho - y ) ; // be done if we had data. But we don't. | |
1070 | z = rho * TMath::Tan(lambda) * TMath::ATan(a) ; | |
1071 | e = TMath::Sqrt( 1 + a*a + b*b ) ; | |
1072 | ||
1073 | layerThickness = radLength / TMath::Sin(TMath::Pi()/2 - lambda) ; | |
1074 | ||
1075 | aMCS = ThetaMCS(mass[massloop], layerThickness, momentum) ; | |
1076 | mmm = ( 1 + a*a + b*b ) ; | |
1077 | bigA(0,0) = mmm*(1+a*a) ; bigA(0,1) = mmm*a*b ; littleA(0,0) = mIstar(2,2) ; littleA(0,1) = mIstar(2,3) ; | |
1078 | bigA(1,0) = mmm*a*b ; bigA(1,1) = mmm*(1+b*b) ; littleA(1,0) = mIstar(3,2) ; littleA(1,1) = mIstar(3,3) ; | |
1079 | littleA = bigA.Invert() + aMCS*aMCS*littleA ; | |
1080 | littleA = littleA.Invert() ; | |
1081 | mA(0,0) = 0.0 ; mA(0,1) = 0.0 ; mA(0,2) = 0.0 ; mA(0,3) = 0.0 ; mA(0,4) = 0.0 ; | |
1082 | mA(1,0) = 0.0 ; mA(1,1) = 0.0 ; mA(1,2) = 0.0 ; mA(1,3) = 0.0 ; mA(1,4) = 0.0 ; | |
1083 | mA(2,0) = 0.0 ; mA(2,1) = 0.0 ; mA(2,2) = littleA(0,0) ; mA(2,3) = littleA(0,1) ; mA(2,4) = 0.0 ; | |
1084 | mA(3,0) = 0.0 ; mA(3,1) = 0.0 ; mA(3,2) = littleA(1,0) ; mA(3,3) = littleA(1,1) ; mA(3,4) = 0.0 ; | |
1085 | mA(4,0) = 0.0 ; mA(4,1) = 0.0 ; mA(4,2) = 0.0 ; mA(4,3) = 0.0 ; mA(4,4) = 0.0 ; | |
1086 | mIstar = mIstar - aMCS*aMCS*mIstar*mA*mIstar ; | |
1087 | dx = ( radius - nextRadius ) ; // Billoir works with absolute magnitude of distance | |
1088 | vF3 = e * ( -1 * ( 1 + a*a ) * fBField ) ; // Assume fBField is on Z axis, only. | |
1089 | vF1 = d * ( a*vF3/(e*e) - 2*e*a*fBField ) ; | |
1090 | vF2 = d * ( b*vF3/(e*e) ) ; | |
1091 | vG3 = e * ( -1 * a * b * fBField ) ; | |
1092 | vG1 = d * ( a*vG3/(e*e) + e*b*fBField ) ; | |
1093 | vG2 = d * ( b*vG3/(e*e) - e*a*fBField ) ; | |
1094 | mD(0,0) = 1.0 ; mD(0,1) = 0.0 ; mD(0,2) = -1*dx+vF1*dx*dx/2 ; mD(0,3) = vF2*dx*dx/2 ; mD(0,4) = vF3*dx*dx/2 ; | |
1095 | mD(1,0) = 0.0 ; mD(1,1) = 1.0 ; mD(1,2) = vG1*dx*dx/2 ; mD(1,3) = -1*dx+vG2*dx*dx/2 ; mD(1,4) = vG3*dx*dx/2 ; | |
1096 | mD(2,0) = 0.0 ; mD(2,1) = 0.0 ; mD(2,2) = 1-vF1*dx ; mD(2,3) = -1*vF2*dx ; mD(2,4) = -1*vF3*dx ; | |
1097 | mD(3,0) = 0.0 ; mD(3,1) = 0.0 ; mD(3,2) = -1*vG1*dx ; mD(3,3) = 1-vG2*dx ; mD(3,4) = -1*vG3*dx ; | |
1098 | mD(4,0) = 0.0 ; mD(4,1) = 0.0 ; mD(4,2) = 0.0 ; mD(4,3) = 0.0 ; mD(4,4) = 1.0 ; | |
1099 | mDT.Transpose(mD) ; | |
1100 | mIstar = mDT.Invert() * mIstar * mD.Invert() ; | |
1101 | // Prepare to save Detector efficiencies | |
1102 | mWork = mIstar ; // Working copy of matrix | |
1103 | mWork.Invert() ; // Invert the Matrix to recover the convariance matrix | |
1104 | fDetPointRes [j+1][i] = TMath::Sqrt( mWork(0,0) ) ; // result in meters | |
1105 | fDetPointZRes[j+1][i] = TMath::Sqrt( mWork(1,1) ) ; // result in meters | |
1106 | // End save | |
1107 | sigma2 = ( nextPhiResolution * nextPhiResolution ) ; | |
1108 | sigma2Z = ( nextZResolution * nextZResolution ) ; | |
1109 | mM(0,0) = 1/sigma2; mM(1,1) = 1/sigma2Z ; | |
1110 | mIstar = mIstar + mM ; | |
1111 | } | |
1112 | ||
1113 | // values below NOT REALIABLE -> they do not point to the vertex but outwards !!!!!!! | |
1114 | // ++++++++++++++ | |
1115 | // also update the values for the track position ?????? | |
1116 | /* | |
1117 | // Pattern recognition is done .... save values like vertex resolution etc. | |
1118 | ||
1119 | // Invert the Matrix to recover the convariance matrix | |
1120 | mIstar.Invert() ; | |
1121 | // Convert the Convariance matrix parameters into physical quantities | |
1122 | // The results are propogated to the previous point but *do not* include the measurement at that point. | |
1123 | deltaPoverP = TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3 ; // Absolute magnitude so ignore charge | |
1124 | fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ) ; // results in percent | |
1125 | fResolutionRPhi[i] = TMath::Sqrt( mIstar(0,0) ) * 1.e6 ; // result in microns | |
1126 | fResolutionZ[i] = TMath::Sqrt( mIstar(1,1) ) * 1.e6 ; // result in microns | |
1127 | // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius | |
1128 | */ | |
1129 | ||
1130 | // Weighted combination of the forward and backward estimates | |
1131 | if (!doLikeAliRoot) { | |
1132 | for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { | |
1133 | fDetPointRes[j][i] = 1/(1/detPointResForw[j][i] + 1/fDetPointRes[j][i]); | |
1134 | fDetPointZRes[j][i] = 1/(1/detPointZResForw[j][i] + 1/fDetPointZRes[j][i]); | |
1135 | } | |
1136 | } | |
1137 | // Set Detector-Efficiency Storage area to unity | |
1138 | fEfficiency[massloop][i] = 1.0 ; | |
1139 | ||
1140 | // print out and efficiency calculation | |
1141 | for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop | |
1142 | ||
1143 | layer = (CylLayer*)fLayers.At(j); | |
1144 | ||
1145 | // Convert to Meters, Tesla, and GeV | |
1146 | Float_t radius = layer->radius /100; | |
1147 | Float_t phiRes = layer->phiRes /100; | |
1148 | Float_t zRes = layer->zRes /100; | |
1149 | Float_t radLength = layer->radL; | |
1150 | Float_t leff = layer->eff; | |
1151 | Bool_t isDead = layer->isDead; | |
1152 | ||
1153 | ||
1154 | if ( (!isDead && radLength >0) ) { | |
1155 | Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] + | |
1156 | phiRes * phiRes ) * 100. ; // work in cm | |
1157 | Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] + | |
1158 | zRes * zRes ) * 100. ; // work in cm | |
1159 | ||
1160 | Double_t layerEfficiency = 0; | |
1161 | if ( EfficiencySearchFlag == 0 ) | |
1162 | layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ; | |
1163 | else if ( EfficiencySearchFlag == 1 ) | |
1164 | layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ; | |
1165 | else if ( EfficiencySearchFlag == 2 ) | |
1166 | layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; | |
1167 | ||
1168 | TString name(layer->GetName()); | |
1169 | if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue; | |
1170 | ||
1171 | if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { | |
1172 | ||
1173 | ||
1174 | printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ", | |
1175 | layer->GetName(), radius*100, radLength, | |
1176 | fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6, | |
1177 | phiRes*1.e6, zRes*1.e6, | |
1178 | HitDensity(radius*100)) ; | |
1179 | if (!name.Contains("tpc")) | |
1180 | printf("%10.3f\n", layerEfficiency); | |
1181 | else | |
1182 | printf(" - \n"); | |
1183 | } | |
1184 | ||
1185 | if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency; | |
1186 | ||
1187 | } | |
1188 | } | |
1189 | if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { | |
1190 | printOnce = 0 ; | |
1191 | printf("\n") ; | |
1192 | } | |
1193 | } | |
54fc064a | 1194 | } // mass loop |
1195 | } // pt loop | |
1196 | ||
1197 | } | |
1198 | ||
1199 | ||
1200 | TGraph * Detector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) { | |
c40d3fcb | 1201 | // |
1202 | // returns the momentum resolution | |
1203 | // | |
54fc064a | 1204 | |
c40d3fcb | 1205 | TGraph *graph = new TGraph(400, fTransMomenta, fMomentumRes); |
54fc064a | 1206 | graph->SetTitle("Momentum Resolution .vs. Pt" ) ; |
c40d3fcb | 1207 | // graph->GetXaxis()->SetRangeUser(0.,5.0) ; |
54fc064a | 1208 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; |
1209 | graph->GetXaxis()->CenterTitle(); | |
c40d3fcb | 1210 | graph->GetXaxis()->SetNoExponent(1) ; |
1211 | graph->GetXaxis()->SetMoreLogLabels(1) ; | |
54fc064a | 1212 | graph->GetYaxis()->SetTitle("Momentum Resolution [%]") ; |
1213 | graph->GetYaxis()->CenterTitle(); | |
c40d3fcb | 1214 | |
1215 | graph->SetMaximum(20) ; | |
1216 | graph->SetMinimum(0.1) ; | |
54fc064a | 1217 | graph->SetLineColor(color); |
1218 | graph->SetMarkerColor(color); | |
1219 | graph->SetLineWidth(linewidth); | |
1220 | ||
1221 | return graph; | |
1222 | ||
1223 | } | |
1224 | ||
1225 | TGraph * Detector::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) { | |
1226 | ||
1227 | // Returns the pointing resolution | |
1228 | // axis = 0 ... rphi pointing resolution | |
1229 | // axis = 1 ... z pointing resolution | |
1230 | // | |
1231 | ||
1232 | TGraph * graph = 0; | |
1233 | ||
1234 | if (axis==0) { | |
c40d3fcb | 1235 | graph = new TGraph ( 400, fTransMomenta, fResolutionRPhi ) ; |
54fc064a | 1236 | graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ; |
1237 | graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution [#mum]") ; | |
1238 | } else { | |
c40d3fcb | 1239 | graph = new TGraph ( 400, fTransMomenta, fResolutionZ ) ; |
54fc064a | 1240 | graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ; |
1241 | graph->GetYaxis()->SetTitle("Z Pointing Resolution [#mum]") ; | |
1242 | } | |
1243 | ||
1244 | graph->SetMinimum(1) ; | |
1245 | graph->SetMaximum(300.1) ; | |
1246 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; | |
1247 | graph->GetXaxis()->CenterTitle(); | |
1248 | graph->GetXaxis()->SetNoExponent(1) ; | |
1249 | graph->GetXaxis()->SetMoreLogLabels(1) ; | |
1250 | graph->GetYaxis()->CenterTitle(); | |
1251 | ||
1252 | graph->SetLineWidth(linewidth); | |
1253 | graph->SetLineColor(color); | |
1254 | graph->SetMarkerColor(color); | |
1255 | ||
1256 | return graph; | |
1257 | ||
1258 | } | |
1259 | ||
1260 | ||
1261 | TGraph * Detector::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) { | |
1262 | // | |
1263 | // returns the Pointing resolution (accoring to Telescope equation) | |
1264 | // axis =0 ... in rphi | |
1265 | // axis =1 ... in z | |
1266 | // | |
1267 | ||
1268 | Double_t resolution[400]; | |
1269 | ||
1270 | Double_t layerResolution[2]; | |
1271 | Double_t layerRadius[2]; | |
1272 | Double_t layerThickness[2]; | |
1273 | ||
1274 | Int_t count =0; // search two first active layers | |
c40d3fcb | 1275 | printf("Telescope equation for layers: "); |
54fc064a | 1276 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { |
1277 | CylLayer *l = (CylLayer*)fLayers.At(i); | |
1278 | if (!l->isDead && l->radius>0) { | |
1279 | layerRadius[count] = l->radius; | |
1280 | layerThickness[count] = l->radL; | |
1281 | if (axis==0) { | |
1282 | layerResolution[count] = l->phiRes; | |
1283 | } else { | |
1284 | layerResolution[count] = l->zRes; | |
1285 | } | |
c40d3fcb | 1286 | printf("%s, ",l->GetName()); |
54fc064a | 1287 | count++; |
1288 | } | |
1289 | if (count>=2) break; | |
1290 | } | |
c40d3fcb | 1291 | printf("\n"); |
54fc064a | 1292 | |
c40d3fcb | 1293 | Double_t pt, momentum, thickness,aMCS ; |
54fc064a | 1294 | Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); |
1295 | ||
1296 | for ( Int_t i = 0 ; i < 400 ; i++ ) { | |
1297 | // Reference data as if first two layers were acting all alone | |
c40d3fcb | 1298 | pt = fTransMomenta[i] ; |
54fc064a | 1299 | momentum = pt / TMath::Cos(lambda) ; // Total momentum |
1300 | resolution[i] = layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1] | |
1301 | + layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ; | |
1302 | resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ; | |
1303 | thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ; | |
c40d3fcb | 1304 | aMCS = ThetaMCS(fParticleMass, thickness, momentum) ; |
1305 | resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ; | |
54fc064a | 1306 | resolution[i] = TMath::Sqrt(resolution[i]) * 10000.0 ; // result in microns |
1307 | } | |
1308 | ||
1309 | ||
1310 | ||
c40d3fcb | 1311 | TGraph* graph = new TGraph ( 400, fTransMomenta, resolution ) ; |
54fc064a | 1312 | |
1313 | if (axis==0) { | |
1314 | graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ; | |
1315 | graph->GetYaxis()->SetTitle("RPhi Pointing Resolution [#mum] ") ; | |
1316 | } else { | |
1317 | graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ; | |
1318 | graph->GetYaxis()->SetTitle("Z Pointing Resolution [#mum] ") ; | |
1319 | } | |
1320 | graph->SetMinimum(1) ; | |
1321 | graph->SetMaximum(300.1) ; | |
1322 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; | |
1323 | graph->GetXaxis()->CenterTitle(); | |
1324 | graph->GetXaxis()->SetNoExponent(1) ; | |
1325 | graph->GetXaxis()->SetMoreLogLabels(1) ; | |
1326 | graph->GetYaxis()->CenterTitle(); | |
1327 | ||
1328 | graph->SetLineColor(color); | |
1329 | graph->SetMarkerColor(color); | |
1330 | graph->SetLineStyle(kDashed); | |
1331 | graph->SetLineWidth(linewidth); | |
1332 | ||
1333 | return graph; | |
1334 | ||
1335 | } | |
1336 | ||
1337 | TGraph * Detector::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) { | |
1338 | // | |
1339 | // particle = 0 ... choosen particle (setted particleMass) | |
1340 | // particle = 1 ... Pion | |
1341 | // particle = 2 ... Kaon | |
1342 | // particle = 3 ... D0 | |
1343 | // | |
1344 | Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); | |
1345 | ||
1346 | Double_t particleEfficiency[400]; // with chosen particle mass | |
c40d3fcb | 1347 | Double_t kaonEfficiency[400], pionEfficiency[400], d0efficiency[400]; |
54fc064a | 1348 | Double_t partEfficiency[2][400]; |
1349 | ||
1350 | if (particle != 0) { | |
1351 | // resulting Pion and Kaon efficiency scaled with overall efficiency | |
1352 | Double_t doNotDecayFactor; | |
1353 | for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon | |
1354 | ||
1355 | for ( Int_t j = 0 ; j < 400 ; j++ ) { | |
1356 | // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm. | |
c40d3fcb | 1357 | Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity |
54fc064a | 1358 | if ( massloop == 1 ) { // KAON |
1359 | doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm. | |
c40d3fcb | 1360 | kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ; |
54fc064a | 1361 | } else { // PION |
1362 | doNotDecayFactor = 1.0 ; | |
c40d3fcb | 1363 | pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ; |
54fc064a | 1364 | } |
1365 | partEfficiency[0][j] = pionEfficiency[j]; | |
1366 | partEfficiency[1][j] = kaonEfficiency[j]; | |
1367 | } | |
1368 | } | |
1369 | ||
1370 | // resulting estimate of the D0 efficiency | |
1371 | for ( Int_t j = 0 ; j < 400 ; j++ ) { | |
c40d3fcb | 1372 | d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency); |
54fc064a | 1373 | } |
1374 | } else { | |
1375 | for ( Int_t j = 0 ; j < 400 ; j++ ) { | |
c40d3fcb | 1376 | particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi; |
54fc064a | 1377 | // NOTE: Decay factor (see kaon) should be included to be realiable |
1378 | } | |
1379 | } | |
1380 | ||
1381 | TGraph * graph = 0; | |
1382 | if (particle==0) { | |
c40d3fcb | 1383 | graph = new TGraph ( 400, fTransMomenta, particleEfficiency ) ; // choosen mass |
54fc064a | 1384 | graph->SetLineWidth(1); |
1385 | } else if (particle==1) { | |
c40d3fcb | 1386 | graph = new TGraph ( 400, fTransMomenta, pionEfficiency ) ; |
54fc064a | 1387 | graph->SetLineWidth(1); |
1388 | } else if (particle ==2) { | |
c40d3fcb | 1389 | graph = new TGraph ( 400, fTransMomenta, kaonEfficiency ) ; |
54fc064a | 1390 | graph->SetLineWidth(1); |
1391 | } else if (particle ==3) { | |
c40d3fcb | 1392 | graph = new TGraph ( 400, fTransMomenta, d0efficiency ) ; |
54fc064a | 1393 | graph->SetLineStyle(kDashed); |
1394 | } else | |
1395 | return 0; | |
1396 | ||
1397 | graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ; | |
1398 | graph->GetXaxis()->CenterTitle(); | |
1399 | graph->GetXaxis()->SetNoExponent(1) ; | |
1400 | graph->GetXaxis()->SetMoreLogLabels(1) ; | |
1401 | graph->GetYaxis()->SetTitle("Efficiency (arbitrary units)") ; | |
1402 | graph->GetYaxis()->CenterTitle(); | |
1403 | ||
1404 | graph->SetMinimum(0.01) ; | |
1405 | graph->SetMaximum(1.0) ; | |
1406 | ||
1407 | graph->SetLineColor(color); | |
1408 | graph->SetMarkerColor(color); | |
1409 | graph->SetLineWidth(linewidth); | |
1410 | ||
1411 | return graph; | |
1412 | } | |
1413 | ||
1414 | ||
1415 | TGraph* Detector::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) { | |
1416 | // | |
1417 | // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution) | |
1418 | // mode 0: impact parameter (convolution of pointing and vertex resolution) | |
1419 | // mode 1: pointing resolution | |
1420 | // mode 2: vtx resolution | |
1421 | ||
1422 | ||
1423 | TGraph *graph = new TGraph(); | |
1424 | ||
1425 | // TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ? | |
1426 | TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); // | |
c40d3fcb | 1427 | TFormula vtxResZ("vtxResZ","600/(x+6)+10"); // |
54fc064a | 1428 | |
1429 | TGraph *trackRes = GetGraphPointingResolution(axis,1); | |
1430 | Double_t *pt = trackRes->GetX(); | |
1431 | Double_t *trRes = trackRes->GetY(); | |
1432 | for (Int_t ip =0; ip<trackRes->GetN(); ip++) { | |
1433 | Double_t vtxRes = 0; | |
1434 | if (axis==0) | |
1435 | vtxRes = vtxResRPhi.Eval(pt[ip]); | |
1436 | else | |
1437 | vtxRes = vtxResZ.Eval(pt[ip]); | |
1438 | ||
1439 | if (mode==0) | |
1440 | graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip])); | |
1441 | else if (mode ==1) | |
1442 | graph->SetPoint(ip,pt[ip],trRes[ip]); | |
1443 | else | |
1444 | graph->SetPoint(ip,pt[ip],vtxRes); | |
1445 | } | |
1446 | ||
1447 | graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ; | |
1448 | graph->GetYaxis()->SetTitle("d_{0} r#phi resolution [#mum]") ; | |
1449 | ||
1450 | graph->SetMinimum(1) ; | |
1451 | graph->SetMaximum(300.1) ; | |
1452 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; | |
1453 | graph->GetXaxis()->CenterTitle(); | |
1454 | graph->GetXaxis()->SetNoExponent(1) ; | |
1455 | graph->GetXaxis()->SetMoreLogLabels(1) ; | |
1456 | graph->GetYaxis()->CenterTitle(); | |
1457 | ||
1458 | graph->SetLineColor(color); | |
1459 | graph->SetMarkerColor(color); | |
1460 | graph->SetLineWidth(linewidth); | |
1461 | ||
1462 | return graph; | |
1463 | ||
1464 | } | |
1465 | ||
1466 | TGraph* Detector::GetGraph(Int_t number, Int_t color, Int_t linewidth) { | |
1467 | // | |
1468 | // returns graph according to the number | |
1469 | // | |
1470 | switch(number) { | |
1471 | case 1: | |
1472 | return GetGraphPointingResolution(0,color, linewidth); // dr | |
1473 | case 2: | |
1474 | return GetGraphPointingResolution(1,color, linewidth); // dz | |
1475 | case 3: | |
1476 | return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele | |
1477 | case 4: | |
1478 | return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele | |
1479 | case 5: | |
1480 | return GetGraphMomentumResolution(color, linewidth); // pt resolution | |
c40d3fcb | 1481 | case 10: |
1482 | return GetGraphRecoEfficiency(0, color, linewidth); // tracked particle | |
54fc064a | 1483 | case 11: |
1484 | return GetGraphRecoEfficiency(1, color, linewidth); // eff. pion | |
1485 | case 12: | |
1486 | return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon | |
1487 | case 13: | |
1488 | return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0 | |
1489 | default: | |
1490 | printf(" Error: chosen graph number not valid\n"); | |
1491 | } | |
1492 | return 0; | |
1493 | ||
1494 | } | |
1495 | ||
1496 | ||
1497 | void Detector::MakeAliceCurrent(Int_t AlignResiduals, Bool_t flagTPC) { | |
1498 | ||
1499 | // Numbers taken from | |
1500 | // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks | |
1501 | // number for misalingment: private communication with Andrea Dainese | |
1502 | ||
1503 | AddLayer((char*)"bpipe",2.94,0.0022); // beam pipe | |
1504 | AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation | |
1505 | AddLayer((char*)"tshld1",11.5,0.0065); // Thermal shield // 1.3% /2 | |
1506 | AddLayer((char*)"tshld2",31.0,0.0065); // Thermal shield // 1.3% /2 | |
54fc064a | 1507 | |
54fc064a | 1508 | |
c40d3fcb | 1509 | if (flagTPC) { |
1510 | AddTPC(0.1,0.1); // TPC | |
1511 | } | |
54fc064a | 1512 | // Adding the ITS - current configuration |
1513 | ||
1514 | if (AlignResiduals==0) { | |
1515 | ||
c40d3fcb | 1516 | AddLayer((char*)"spd1", 3.9, 0.0114, 0.0012, 0.0130); |
1517 | AddLayer((char*)"spd2", 7.6, 0.0114, 0.0012, 0.0130); | |
54fc064a | 1518 | AddLayer((char*)"sdd1",15.0, 0.0113, 0.0035, 0.0025); |
1519 | AddLayer((char*)"sdd2",23.9, 0.0126, 0.0035, 0.0025); | |
1520 | AddLayer((char*)"ssd1",38.0, 0.0083, 0.0020, 0.0830); | |
1521 | AddLayer((char*)"ssd2",43.0, 0.0086, 0.0020, 0.0830); | |
1522 | ||
1523 | } else if (AlignResiduals==1) { | |
1524 | ||
1525 | // tracking errors ... | |
1526 | // (Additional systematic errors due to misalignments) ... | |
1527 | // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm] | |
1528 | // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000); | |
1529 | ||
1530 | AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010), | |
c40d3fcb | 1531 | TMath::Sqrt(0.0130*0.0130+0.0050*0.0050)); |
54fc064a | 1532 | AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030), |
c40d3fcb | 1533 | TMath::Sqrt(0.0130*0.0130+0.0050*0.0050)); |
54fc064a | 1534 | AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500), |
1535 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); | |
1536 | AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500), | |
1537 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); | |
1538 | AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020), | |
1539 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); | |
1540 | AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020), | |
1541 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); | |
1542 | ||
1543 | } else if (AlignResiduals==2) { | |
1544 | ||
1545 | // tracking errors ... PLUS ... module misalignment | |
1546 | ||
1547 | // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm] | |
1548 | // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000); | |
1549 | ||
1550 | // the ITS modules are misalignment with small gaussian smearings with | |
1551 | // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD | |
1552 | ||
1553 | AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008), | |
c40d3fcb | 1554 | TMath::Sqrt(0.0130*0.0130+0.0050*0.0050)); |
54fc064a | 1555 | AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008), |
c40d3fcb | 1556 | TMath::Sqrt(0.0130*0.0130+0.0050*0.0050)); |
54fc064a | 1557 | AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010), |
1558 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); | |
1559 | AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010), | |
1560 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); | |
1561 | AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010), | |
1562 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); | |
1563 | AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010), | |
1564 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); | |
1565 | ||
1566 | } else { | |
1567 | ||
1568 | // the ITS modules are misalignment with small gaussian smearings with | |
1569 | // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD | |
1570 | // unknown in Z ???? | |
1571 | ||
1572 | AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008), | |
c40d3fcb | 1573 | TMath::Sqrt(0.0130*0.0130+0.000*0.000)); |
54fc064a | 1574 | AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008), |
c40d3fcb | 1575 | TMath::Sqrt(0.0130*0.0130+0.000*0.000)); |
54fc064a | 1576 | AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010), |
1577 | TMath::Sqrt(0.0025*0.0025+0.000*0.000)); | |
1578 | AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010), | |
1579 | TMath::Sqrt(0.0025*0.0025+0.000*0.000)); | |
1580 | AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010), | |
1581 | TMath::Sqrt(0.0830*0.0830+0.000*0.000)); | |
1582 | AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010), | |
1583 | TMath::Sqrt(0.0830*0.0830+0.000*0.000)); | |
1584 | ||
1585 | ||
1586 | } | |
1587 | ||
1588 | } | |
c40d3fcb | 1589 | |
1590 | ||
1591 | void Detector::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) { | |
1592 | // | |
1593 | // Produces the standard performace plots | |
1594 | // | |
1595 | ||
1596 | if (!add) { | |
1597 | ||
1598 | TCanvas *c1 = new TCanvas("c1","c1"); | |
1599 | c1->Divide(2,2); | |
1600 | ||
1601 | c1->cd(1); gPad->SetGridx(); gPad->SetGridy(); | |
1602 | gPad->SetLogx(); | |
1603 | TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth); | |
1604 | eff->SetTitle("Efficiencies"); | |
1605 | eff->Draw("AC"); | |
1606 | if (!onlyPionEff) { | |
1607 | GetGraphRecoEfficiency(2,color,linewidth)->Draw("C"); | |
1608 | GetGraphRecoEfficiency(3,color,linewidth)->Draw("C"); | |
1609 | } | |
1610 | c1->cd(2); gPad->SetGridx(); gPad->SetGridy(); | |
1611 | gPad->SetLogy(); gPad->SetLogx(); | |
1612 | GetGraphMomentumResolution(color,linewidth)->Draw("AC"); | |
1613 | ||
1614 | c1->cd(3); gPad->SetGridx(); gPad->SetGridy(); | |
1615 | gPad->SetLogx(); | |
1616 | GetGraphPointingResolution(0,color,linewidth)->Draw("AC"); | |
1617 | ||
1618 | c1->cd(4); gPad->SetGridx(); gPad->SetGridy(); | |
1619 | gPad->SetLogx(); | |
1620 | GetGraphPointingResolution(1,color,linewidth)->Draw("AC"); | |
1621 | ||
1622 | } else { | |
1623 | ||
1624 | TVirtualPad *c1 = gPad->GetMother(); | |
1625 | ||
1626 | c1->cd(1); | |
1627 | GetGraphRecoEfficiency(1,color,linewidth)->Draw("C"); | |
1628 | if (!onlyPionEff) { | |
1629 | GetGraphRecoEfficiency(2,color,linewidth)->Draw("C"); | |
1630 | GetGraphRecoEfficiency(3,color,linewidth)->Draw("C"); | |
1631 | } | |
1632 | c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("C"); | |
1633 | ||
1634 | c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("C"); | |
1635 | ||
1636 | c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("C"); | |
1637 | ||
1638 | } | |
1639 | ||
1640 | } |