]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/Detector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing compilation warnings
[u/mrichter/AliRoot.git] / ITS / UPGRADE / Detector.cxx
CommitLineData
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
14Fast Simulation tool for Inner Tracker Systems
15
16original code of using the billoir technique was developed
17for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov
18http://rnc.lbl.gov/~jhthomas
19
c40d3fcb 20Changes 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
45class CylLayer : public TNamed {
46public:
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
64ClassImp(Detector)
65Detector::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
84Detector::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}
101Detector::~Detector() { //
102 // virtual destructor
103 //
104 // delete fLayers;
105}
106
c40d3fcb 107void 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
155void 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
173void 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
194Float_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
208void 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
221Float_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
236void 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
263Float_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 281void 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
294Float_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 309void 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
326void 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 357void 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 399void 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
445void 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 459Double_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 473Double_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 486Double_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 496Double_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
511Double_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 542double 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 552double 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 567double 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
581double 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 605Double_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 658void 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 719void 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
1200TGraph * 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
1225TGraph * 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
1261TGraph * 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
1337TGraph * 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
1415TGraph* 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
1466TGraph* 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
1497void 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
1591void 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}