]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/DetectorK.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
commit to clean dev branch
[u/mrichter/AliRoot.git] / ITS / UPGRADE / DetectorK.cxx
CommitLineData
84356182 1#include "DetectorK.h"
2#include <TMath.h>
3#include <TMatrixD.h>
4#include <TGraph.h>
5#include <TAxis.h>
6#include <TFormula.h>
7#include <TCanvas.h>
8#include <TEllipse.h>
9#include <TText.h>
10#include <TGraphErrors.h>
11
12#include "AliExternalTrackParam.h"
13
14/***********************************************************
15
16Fast Simulation tool for Inner Tracker Systems
17
18original code of using the billoir technique was developed
19for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov
20http://rnc.lbl.gov/~jhthomas
21
22Changes by S. Rossegger -> see header file
23
24***********************************************************/
25
26
27#define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector
28
29#define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )
30#define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm)
31#define dNdEtaMinB 1//950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700)
32// #define dNdEtaCent 2300//15000 //1600//2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known)
33
34#define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)
35#define AcceptanceOfTpcAndSi 1 //1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors
36#define UPCBackgroundMultiplier 1.0 // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)
37#define OtherBackground 0.0 // Increase multiplicity in detector (0.0 to 1.0 * minBias) (eg 0.0)
38#define EfficiencySearchFlag 2 // Define search method:
39 // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0.
40
41#define PionMass 0.139 // Mass of the Pion
42#define KaonMass 0.498 // Mass of the Kaon
43#define D0Mass 1.865 // Mass of the D0
44
45//TMatrixD *probKomb; // table for efficiency kombinatorics
46
47
48class CylLayerK : public TNamed {
49public:
50
51 CylLayerK(char *name) : TNamed(name,name) {}
52
53 Float_t GetRadius() const {return radius;}
54 Float_t GetRadL() const {return radL;}
55 Float_t GetPhiRes() const {return phiRes;}
56 Float_t GetZRes() const {return zRes;}
57 Float_t GetLayerEff() const {return eff;}
58
59 // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }
60 Float_t radius; Float_t radL; Float_t phiRes; Float_t zRes;
61 Float_t eff;
62 Bool_t isDead;
63
64 ClassDef(CylLayerK,1);
65};
66
67
68class ForwardLayer : public TNamed {
69public:
70 ForwardLayer(char *name) : TNamed(name,name) {}
71
72 Float_t GetZ() const {return zPos;}
73 Float_t GetXRes() const {return xRes;}
74 Float_t GetYRes() const {return yRes;}
75 Float_t GetThickness() const {return thickness;}
76 Float_t Getdensity() const {return density;}
77 Float_t GetLayerEff() const {return eff;}
78
79 // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }
80 Float_t zPos; Float_t xRes; Float_t yRes;
81 Float_t radL;
82 Float_t thickness;
83 Float_t density;
84 Float_t eff;
85 Bool_t isDead;
86
87 ClassDef(ForwardLayer,1);
88};
89
90
91ClassImp(DetectorK)
92DetectorK::DetectorK()
93 : TNamed("test_detector","detector"),
94 fNumberOfLayers(0),
95 fNumberOfActiveLayers(0),
96 fNumberOfActiveITSLayers(0),
97 fBField(0.5),
98 fLhcUPCscale(1.0),
99 fIntegrationTime(0.02), // in ms
100 fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
101 fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
102 fParticleMass(0.140), // Standard: pion mass
103 fMaxRadiusSlowDet(10.),
104 fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
105 fAtLeastFake(1), // if at least x fakes, track is considered fake ...
106 fMaxSeedRadius(50000),
107 fptScale(10.),
108 fdNdEtaCent(2300)
109{
110 //
111 // default constructor
112 //
113 // fLayers = new TObjArray();
114
115}
116
117DetectorK::DetectorK(char *name, char *title)
118 : TNamed(name,title),
119 fNumberOfLayers(0),
120 fNumberOfActiveLayers(0),
121 fNumberOfActiveITSLayers(0),
122 fBField(0.5),
123 fLhcUPCscale(1.0),
124 fIntegrationTime(0.02), // in ms
125 fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
126 fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
127 fParticleMass(0.140), // Standard: pion mass
128 fMaxRadiusSlowDet(10.),
129 fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
130 fAtLeastFake(1), // if at least x fakes, track is considered fake ...
131 fMaxSeedRadius(50000),
132 fptScale(10.),
133 fdNdEtaCent(2200)
134{
135 //
136 // default constructor, that set the name and title
137 //
138 // fLayers = new TObjArray();
139}
140DetectorK::~DetectorK() { //
141 // virtual destructor
142 //
143 // delete fLayers;
144}
145
146void DetectorK::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes, Float_t zRes, Float_t eff) {
147 //
148 // Add additional layer to the list of layers (ordered by radius)
149 //
150
151 CylLayerK *newLayer = (CylLayerK*) fLayers.FindObject(name);
152
153 if (!newLayer) {
154 newLayer = new CylLayerK(name);
155 newLayer->radius = radius;
156 newLayer->radL = radL;
157 newLayer->phiRes = phiRes;
158 newLayer->zRes = zRes;
159 newLayer->eff = eff;
160
161 if (newLayer->zRes==RIDICULOUS && newLayer->zRes==RIDICULOUS)
162 newLayer->isDead = kTRUE;
163 else
164 newLayer->isDead = kFALSE;
165
166 if (fLayers.GetEntries()==0)
167 fLayers.Add(newLayer);
168 else {
169
170 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
171 CylLayerK *l = (CylLayerK*)fLayers.At(i);
172 if (radius<l->radius) {
173 fLayers.AddBefore(l,newLayer);
174 break;
175 }
176 if (radius>l->radius && (i+1)==fLayers.GetEntries() ) {
177 // even bigger then last one
178 fLayers.Add(newLayer);
179 }
180 }
181
182 }
183 fNumberOfLayers += 1;
184 if (!(newLayer->isDead)) {
185 fNumberOfActiveLayers += 1;
186 TString lname(newLayer->GetName());
187 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;
188 }
189
190
191 } else {
192 printf("Layer with the name %s does already exist\n",name);
193 }
194
195
196}
197
198void DetectorK::KillLayer(char *name) {
199 //
200 // Marks layer as dead. Contribution only by Material Budget
201 //
202
203 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
204 if (!tmp)
205 printf("Layer %s not found - cannot mark as dead\n",name);
206 else {
207 tmp->phiRes = 999999;
208 tmp->zRes = 999999;
209 if (!(tmp->isDead)) {
210 tmp->isDead = kTRUE;
211 fNumberOfActiveLayers -= 1;
212 TString lname(tmp->GetName());
213 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
214 }
215 }
216}
217
218void DetectorK::SetRadius(char *name, Float_t radius) {
219 //
220 // Set layer radius [cm]
221 //
222
223 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
224
225
226 if (!tmp) {
227 printf("Layer %s not found - cannot set radius\n",name);
228 } else {
229
230 Float_t tmpRadL = tmp->radL;
231 Float_t tmpPhiRes = tmp->phiRes;
232 Float_t tmpZRes = tmp->zRes;
233
234 RemoveLayer(name); // so that the ordering is correct
235 AddLayer(name,radius,tmpRadL,tmpPhiRes,tmpZRes);
236 }
237}
238
239Float_t DetectorK::GetRadius(char *name) {
240 //
241 // Return layer radius [cm]
242 //
243
244 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
245 if (!tmp)
246 printf("Layer %s not found - cannot get radius\n",name);
247 else
248 return tmp->radius;
249
250 return 0;
251}
252
253void DetectorK::SetRadiationLength(char *name, Float_t radL) {
254 //
255 // Set layer material [cm]
256 //
257
258 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
259 if (!tmp)
260 printf("Layer %s not found - cannot set layer material\n",name);
261 else {
262 tmp->radL = radL;
263 }
264}
265
266Float_t DetectorK::GetRadiationLength(char *name) {
267 //
268 // Return layer radius [cm]
269 //
270
271 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
272 if (!tmp)
273 printf("Layer %s not found - cannot get layer material\n",name);
274 else
275 return tmp->radL;
276
277 return 0;
278
279}
280
281void DetectorK::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
282 //
283 // Set layer resolution in [cm]
284 //
285
286 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
287 if (!tmp)
288 printf("Layer %s not found - cannot set resolution\n",name);
289 else {
290
291 Bool_t wasDead = tmp->isDead;
292
293 tmp->phiRes = phiRes;
294 tmp->zRes = zRes;
295 TString lname(tmp->GetName());
296
297 if (zRes==RIDICULOUS && phiRes==RIDICULOUS) {
298 tmp->isDead = kTRUE;
299 if (!wasDead) {
300 fNumberOfActiveLayers -= 1;
301 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
302 }
303 } else {
304 tmp->isDead = kFALSE;
305 if (wasDead) {
306 fNumberOfActiveLayers += 1;
307 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;
308 }
309 }
310
311
312 }
313}
314
315Float_t DetectorK::GetResolution(char *name, Int_t axis) {
316 //
317 // Return layer resolution in [cm]
318 // axis = 0: resolution in rphi
319 // axis = 1: resolution in z
320 //
321
322 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
323 if (!tmp)
324 printf("Layer %s not found - cannot get resolution\n",name);
325 else {
326 if (axis==0) return tmp->phiRes;
327 if (axis==1) return tmp->zRes;
328 printf("error: axis must be either 0 or 1 (rphi or z axis)\n");
329 }
330 return 0;
331}
332
333void DetectorK::SetLayerEfficiency(char *name, Float_t eff) {
334 //
335 // Set layer efficnecy (prop that his is missed within this layer)
336 //
337
338 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
339 if (!tmp)
340 printf("Layer %s not found - cannot set layer efficiency\n",name);
341 else {
342 tmp->eff = eff;
343 }
344}
345
346Float_t DetectorK::GetLayerEfficiency(char *name) {
347 //
348 // Get layer efficnecy (prop that his is missed within this layer)
349 //
350
351 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
352 if (!tmp)
353 printf("Layer %s not found - cannot get layer efficneicy\n",name);
354 else
355 return tmp->eff;
356
357 return 0;
358
359}
360
361void DetectorK::RemoveLayer(char *name) {
362 //
363 // Removes a layer from the list
364 //
365
366 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
367 if (!tmp)
368 printf("Layer %s not found - cannot remove it\n",name);
369 else {
370 Bool_t wasDead = tmp->isDead;
371 fLayers.Remove(tmp);
372 fNumberOfLayers -= 1;
373 if (!wasDead) {
374 fNumberOfActiveLayers -= 1;
375 TString lname(tmp->GetName());
376 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
377
378 }
379 }
380}
381
382
383void DetectorK::PrintLayout() {
384 //
385 // Prints the detector layout
386 //
387
388 printf("Detector %s: \"%s\"\n",GetName(),GetTitle());
389
390 if (fLayers.GetEntries()>0)
391 printf(" Name \t\t r [cm] \t X0 \t phi & z res [um] layerEff \n");
392
393 CylLayerK *tmp = 0;
394 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
395 tmp = (CylLayerK*)fLayers.At(i);
396
397 // don't print all the tpc layers
398 TString name(tmp->GetName());
399 if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
400
401 printf("%d. %s \t %03.2f \t%1.4f\t ",i,
402 tmp->GetName(), tmp->radius, tmp->radL);
403 if (tmp->phiRes==RIDICULOUS)
404 printf(" - ");
405 else
406 printf("%3.0f ",tmp->phiRes*10000);
407 if (tmp->zRes==RIDICULOUS)
408 printf(" -");
409 else
410 printf("%3.0f",tmp->zRes*10000);
411
412 if (tmp->zRes==RIDICULOUS)
413 printf("\t -\n");
414 else
415 printf("\t%0.2f\n",tmp->eff);
416
417 }
418}
419
420void DetectorK::PlotLayout(Int_t plotDead) {
421 //
422 // Plots the detector layout in Front view
423 //
424
425 Double_t x0=0, y0=0;
426
427 TGraphErrors *gr = new TGraphErrors();
428 gr->SetPoint(0,0,0);
429 CylLayerK *lastLayer = (CylLayerK*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->radius;
430 gr->SetPointError(0,maxRad,maxRad);
431 gr->Draw("APE");
432
433
434 CylLayerK *tmp = 0;
435 for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {
436 tmp = (CylLayerK*)fLayers.At(i);
437
438
439 Double_t txtpos = tmp->radius;
440 if ((tmp->isDead)) txtpos*=-1; //
441 TText *txt = new TText(x0,txtpos,tmp->GetName());
442 txt->SetTextSizePixels(5); txt->SetTextAlign(21);
443 if (!tmp->isDead || plotDead) txt->Draw();
444
445 TEllipse *layEl = new TEllipse(x0,y0,tmp->radius);
446 // layEl->SetFillColor(5);
447 layEl->SetFillStyle(5001);
448 layEl->SetLineStyle(tmp->isDead+1); // dashed if not active
449 layEl->SetLineColor(4);
450 TString name(tmp->GetName());
451 if (!tmp->isDead) layEl->SetLineWidth(2);
452 if (name.Contains("tpc") ) layEl->SetLineColor(29);
453
454 if (!tmp->isDead || plotDead) layEl->Draw();
455
456 }
457
458}
459
460
461
462void DetectorK::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {
463 //
464 // Emulates the TPC
465 //
466 // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow
467
468
469 AddLayer((char*)"IFC", 77.8,0.01367); // Inner Field cage
470 AddLayer((char*)"OFC", 254.0,0.01367); // Outer Field cage
471
472 // % Radiation Lengths ... Average per TPC row (i.e. total/159 )
473 Float_t radLBoubdary = 0.0165;
474 Float_t rBoundary = 70.0; // cm
475
476 Float_t radLPerRow = 0.000036;
477
478 Float_t tpcInnerRadialPitch = 0.75 ; // cm
479 Float_t tpcMiddleRadialPitch = 1.0 ; // cm
480 Float_t tpcOuterRadialPitch = 1.5 ; // cm
481 // Float_t tpcInnerPadWidth = 0.4 ; // cm
482 // Float_t tpcMiddlePadWidth = 0.6 ; // cm
483 // Float_t tpcOuterPadWidth = 0.6 ; // cm
484 Float_t innerRows = 63 ;
485 Float_t middleRows = 64 ;
486 Float_t outerRows = 32 ;
487 Float_t tpcRows = (innerRows + middleRows + outerRows) ;
488 Float_t rowOneRadius = 85.2 ; // cm
489 Float_t row64Radius = 135.1 ; // cm
490 Float_t row128Radius = 199.2 ; // cm
491
492 // add boundary between ITS and TPC
493 AddLayer("tpc_boundary",rBoundary,radLBoubdary); // dummy errors
494
495 for ( Int_t k = 0 ; k < tpcRows ; k++ ) {
496
497 Float_t rowRadius =0;
498 if (k<innerRows)
499 rowRadius = rowOneRadius + k*tpcInnerRadialPitch ;
500 else if ( k>=innerRows && k<(innerRows+middleRows) )
501 rowRadius = row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;
502 else if (k>=(innerRows+middleRows) && k<tpcRows )
503 rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;
504
505 if ( k%skip == 0 )
506 AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow,phiResMean,zResMean);
507 else
508 AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow); // non "active" row
509
510
511 }
512
513}
514
515void DetectorK::RemoveTPC() {
516
517 // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)
518 CylLayerK *tmp = 0;
519 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
520 tmp = (CylLayerK*)fLayers.At(i);
521 TString name(tmp->GetName());
522 if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }
523 }
524 RemoveLayer((char*)"IFC");
525
526}
527
528
529Double_t DetectorK::ThetaMCS ( Double_t mass, Double_t radLength, Double_t momentum ) const
530{
531 //
532 // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)
533 //
534
535 Double_t beta = momentum / TMath::Sqrt(momentum*momentum+mass*mass) ;
536 Double_t theta = 0.0 ; // Momentum and mass in GeV
537 // if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );
538 if ( radLength > 0 ) theta = 0.0136 * TMath::Sqrt(radLength) / ( beta * momentum ) * (1+0.038*TMath::Log(radLength)) ;
539 return (theta) ;
540}
541
542
543Double_t DetectorK::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
544{
545 // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm
546 // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm
547 // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius
548 Double_t sx, sy, goodHit ;
549 sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;
550 sy = 2 * TMath::Pi() * searchRadiusZ * searchRadiusZ * HitDensity(radius) ;
551 goodHit = TMath::Sqrt(1./((1+sx)*(1+sy))) ;
552 return ( goodHit ) ;
553}
554
555
556Double_t DetectorK::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
557{
558 // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm
559 // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
560 Double_t sx, goodHit ;
561 sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;
562 goodHit = 1./(1+sx) ;
563 return ( goodHit ) ;
564}
565
566Double_t DetectorK::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
567{
568 // Based on work by Ruben Shahoyen
569 // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
570 // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account
571 // Following is correct for 2 DOF
572
573 Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
574 Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
575 Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c));
576 return ( goodHit ) ;
577}
578
579Double_t DetectorK::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
580{
581 // Based on work by Ruben Shahoyen
582 // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)
583
584 Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
585 Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
586 Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));
587 return ( nullHit ) ;
588}
589
590Double_t DetectorK::HitDensity ( Double_t radius )
591{
592 // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.
593 // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results
594 // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.
595 // Note that this function assumes we are working in CM and CM**2 [not meters].
596 // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.
597
598 // Double_t MaxRadiusSlowDet = 0.1; //? // Maximum radius for slow detectors. Fast detectors
599 // and only fast detectors reside outside this radius.
600 Double_t arealDensity = 0 ;
601
602 if ( radius > fMaxRadiusSlowDet )
603 {
604 arealDensity = OneEventHitDensity(fdNdEtaCent,radius) ; // Fast detectors see central collision density (only)
605 arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius) ; // Increase density due to background
606 }
607
608 if (radius < fMaxRadiusSlowDet )
609 { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.
610 arealDensity = OneEventHitDensity(fdNdEtaCent,radius)
611 + IntegratedHitDensity(dNdEtaMinB,radius)
612 + UpcHitDensity(radius) ;
613 arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;
614 // Increase density due to background
615 }
616
617 return ( arealDensity ) ;
618}
619
620
621double DetectorK::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const
622{
623 // This is for one event at the vertex. No smearing.
624
625 double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?
626 double tg = TMath::Tan(2*TMath::ATan(TMath::Exp(-fAvgRapidity)));
627 den = den/TMath::Sqrt(1 + 1/(tg*tg));
628
629 // double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?
630 // note: surface of sphere is '4*pi*r^2'
631 // surface of cylinder is '2*pi*r* h'
632
633
634
635 return den ;
636}
637
638
639double DetectorK::IntegratedHitDensity(Double_t multiplicity, Double_t radius)
640{
641 // The integral of minBias events smeared over a gaussian vertex distribution.
642 // Based on work by Yan Lu 12/20/2006, all radii in centimeters.
643
644 Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;
645 Double_t den = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;
646
647 // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.
648 if ( den < OneEventHitDensity(multiplicity,radius) ) den = OneEventHitDensity(multiplicity,radius) ;
649
650 return den ;
651}
652
653
654double DetectorK::UpcHitDensity(Double_t radius)
655{
656 // QED electrons ...
657
658 Double_t mUPCelectrons ; ;
659 // mUPCelectrons = fLhcUPCscale * (1.23 - radius/6.5) ; // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC
660 mUPCelectrons = fLhcUPCscale*5456/(radius*radius)/dNdEtaMinB; // Fit to 'Rossegger,Sadovsky'-Alice simulation
661 if ( mUPCelectrons < 0 ) mUPCelectrons = 0.0 ; // UPC electrons fall off quickly and don't go to large R
662 mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ; // UPCs increase Mulitiplicty ~ proportional to MinBias rate
663 mUPCelectrons *= UPCBackgroundMultiplier ; // Allow for an external multiplier (eg 0-1) to turn off UPC
664
665 return mUPCelectrons ;
666}
667
668
669double DetectorK::Dist(double z, double r)
670{
671 // Convolute dEta/dZ distribution with assumed Gaussian of vertex z distribution
672 // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm
673 // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.
674 Int_t index = 1 ; // Start weight at 1 for Simpsons rule integration
675 Int_t nsteps = 301 ; // NSteps must be odd for Simpson's rule to work
676 double dist = 0.0 ;
677 double dz0 = ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1) ; //cm
678 double z0 = 0.0 ; //cm
679 for(int i=0; i<nsteps; i++){
680 if ( i == nsteps-1 ) index = 1 ;
681 z0 = -4*SigmaD + i*dz0 ;
682 dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) *
683 (1/sqrt((z-z0)*(z-z0) + r*r)) ;
684 if ( index != 4 ) index = 4; else index = 2 ;
685 }
686 return dist;
687}
688
689#define PZero 0.861 // Momentum of back to back decay particles in the CM frame
690#define EPiZero 0.872 // Energy of the pion from a D0 decay at rest
691#define EKZero 0.993 // Energy of the Kaon from a D0 decay at rest
692
693Double_t DetectorK::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const {
694 // Math from Ron Longacre. Note hardwired energy to bin conversion for PtK and PtPi.
695
696 Double_t const1 = pt / D0Mass ;
697 Double_t const2 = TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;
698 Double_t sum, ptPi, ptK ;
699 Double_t effp, effk ;
700
701 sum = 0.0 ;
702 for ( Int_t k = 0 ; k < 360 ; k++ ) {
703
704 Double_t theta = k * TMath::Pi() / 180. ;
705
706 ptPi = TMath::Sqrt(
707 PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
708 const1*const1*EPiZero*EPiZero -
709 2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +
710 PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
711 ) ;
712
713 ptK = TMath::Sqrt(
714 PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
715 const1*const1*EKZero*EKZero +
716 2*PZero*TMath::Cos(theta)*const2*const1*EKZero +
717 PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
718 ) ;
719
720 // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays
721 Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
722 Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
723
724 if ( pionindex >= kNptBins ) pionindex = 399 ;
725 if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ;
726 if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd
727 if ( effp < 0 ) effp = 0 ;
728
729 if ( kaonindex >= kNptBins ) kaonindex = 399 ;
730 if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ;
731 if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd
732 if ( effk < 0 ) effk = 0 ;
733
734 // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.
735
736 sum += effp * effk ;
737
738 }
739
740 Double_t mean =sum/360;
741 return mean ;
742
743}
744
745
746
747void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t meanPt, char* detLayer) {
748 //
749 // Solves the current geometry with the Billoir technique
750 // ( see P. Billoir, Nucl. Instr. and Meth. 225 (1984), p. 352. )
751 // ABOVE IS OBSOLETE -> NOW, its uses the Aliroot Kalman technique
752 //
753 const float kTrackingMargin = 0.1;
754
755 static AliExternalTrackParam probTr; // track to propagate
756 probTr.SetUseLogTermMS(kTRUE);
757
758
759 Int_t nPt = kNptBins;
760 // Clean up ......
761 for (Int_t i=0; i<kMaxNumberOfDetectors; i++) {
762 for (Int_t j=0; j<nPt; j++) {
763 fDetPointRes[i][j] = RIDICULOUS;
764 fDetPointZRes[i][j] = RIDICULOUS;
765 fTransMomenta[i] =0;
766 fMomentumRes[i] =0;
767 fResolutionRPhi[i] =0;
768 }
769 }
770
771 if (!allPt) { // not the whole pt range -> allows a faster minimization at a defined 'meanpt'
772 nPt = 3;
773 }
774
775
776 // Calculate track parameters using Billoirs method of matrices
777
778 Double_t pt,tgl, pz, lambda, deltaPoverP ;
779 Double_t charge ;
780 Double_t mass[3] ;
781 Int_t printOnce = 1 ;
782
783 mass[0] = PionMass ; mass[1] = KaonMass ; // Loop twice for the D0; first pi then k
784
785 mass[2] = fParticleMass; // third loop
786
787 Int_t mStart =0;
788 if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode
789
790
791
792
793 // Prepare Probability Kombinations
794 Int_t nLayer = fNumberOfActiveITSLayers;
795 Int_t base = 3; // null, fake, correct
796
797 Int_t komb = (Int_t) TMath::Power(base,nLayer);
798
799 TMatrixD probLay(base,fNumberOfActiveITSLayers);
800 TMatrixD probKomb(komb,nLayer);
801 for (Int_t num=0; num<komb; num++) {
802 for (Int_t l=nLayer; l--;) {
803 Int_t pow = ((Int_t)TMath::Power(base,l+1));
804 probKomb(num,nLayer-1-l)=(num%pow)/((Int_t)TMath::Power(base,l));
805 }
806 }
807
808 TString detLayerStr(detLayer);
809 CylLayerK *theLayer = (CylLayerK*) fLayers.FindObject(detLayer);
810 if (!theLayer && detLayerStr.IsNull()!=1){
811 printf("Error: Layer with the name \"%s\" not found -> no detailed infos possible\n",detLayer);
812 return;
813 }
814 for (Int_t i=0; i<fLayers.GetEntries();i++) {
815 CylLayerK *l = (CylLayerK*) fLayers.At(i);
816 if (detLayerStr.CompareTo(l->GetName())==0) { // is the same
817 kDetLayer=i;
818 break;
819 }
820 }
821
822
823 for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ ) {
824
825 // PseudoRapidity OK, used as an angle
826 lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)) ;
827
828
829 for ( Int_t i = 0 ; i < nPt ; i++ ) { // pt loop
830
831 CylLayerK *last = (CylLayerK*) fLayers.At((fLayers.GetEntries()-1));
832
833 // Starting values based on radius of outermost layer ... log10 steps to ~20 GeV
834 Double_t bigRad = last->radius/2 ; // min. pt which the algorithm below could handle
835 // if (bigRad<61) bigRad=61; // -> min pt around 100 MeV for Bz=0.5T (don't overdo it ... ;-) )
836 fTransMomenta[i] = ( 0.3*bigRad*TMath::Abs(fBField)*1e-2 ) - 0.08 - (1./fptScale-0.1) + TMath::Power(10,2.3*i/nPt) / fptScale ;
837 if (!allPt) { // just 3 points around meanPt
838 fTransMomenta[i] = meanPt-0.001+(Double_t)(i)*0.001;
839 }
840
841 // New from here ................
842
843 // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis
844 // These are the EndPoint values for y, z, a, b, and d
845 double bGauss = fBField*10; // field in kgauss
846 pt = fTransMomenta[i]; // GeV/c
847 tgl = TMath::Tan(lambda); // dip
848 charge = -1; // Assume an electron
849 pz = pt * TMath::Tan(lambda) ; // GeV/
850 enum {kY,kZ,kSnp,kTgl,kPtI}; // track parameter aliases
851 enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2}; // cov.matrix aliases
852 //
853 probTr.Reset();
854 double *trPars = (double*)probTr.GetParameter();
855 double *trCov = (double*)probTr.GetCovariance();
856 trPars[kY] = 0; // start from Y = 0
857 trPars[kZ] = 0; // Z = 0
858 trPars[kSnp] = 0; // track along X axis at the vertex
859 trPars[kTgl] = TMath::Tan(lambda); // dip
860 trPars[kPtI] = charge/pt; // q/pt
861 //
862 // put tiny errors to propagate to the outer radius
863 trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9;
864 if (!PropagateToR(&probTr,last->radius+ kTrackingMargin,bGauss,1)) continue;
865 //if (!probTr.PropagateTo(last->radius,bGauss)) continue;
866 // reset cov.matrix
867 const double kLargeErr2Coord = 5*5;
868 const double kLargeErr2Dir = 0.7*0.7;
869 const double kLargeErr2PtI = 30.5*30.5;
870 for (int ic=15;ic--;) trCov[ic] = 0.;
871 trCov[kY2] = trCov[kZ2] = kLargeErr2Coord;
872 trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
873 trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
874 probTr.CheckCovariance();
875 //
876 // printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d);
877
878 // Set Detector-Efficiency Storage area to unity
879 fEfficiency[massloop][i] = 1.0 ;
880 //
881 // Back-propagate the covariance matrix along the track.
882
883 CylLayerK *layer = 0;
884
885 // find last "active layer" - start tracking at the last active layer
886 Int_t lastActiveLayer = 0;
887 for (Int_t j=fLayers.GetEntries(); j--;) {
888 layer = (CylLayerK*)fLayers.At(j);
889 if (!(layer->isDead)) { // is alive
890 lastActiveLayer = j;
891 break;
892 }
893 }
894 // probTr.Print();
895 for (Int_t j=lastActiveLayer+1; j--;) { // Layer loop
896
897 layer = (CylLayerK*)fLayers.At(j);
898
899 if (layer->radius>fMaxSeedRadius) continue; // no seeding beyond this radius
900
901 TString name(layer->GetName());
902 Bool_t isVertex = name.Contains("vertex");
903 //
904 if (!PropagateToR(&probTr,layer->radius,bGauss,-1)) exit(1);
905 // if (!probTr.PropagateTo(last->radius,bGauss)) exit(1); //
906 // rotate to frame with X axis normal to the surface
907 if (!isVertex) {
908 double pos[3];
909 probTr.GetXYZ(pos); // lab position
910 double phi = TMath::ATan2(pos[1],pos[0]);
911 if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
912 if (!probTr.Rotate(phi)) {
913 printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f (pt=%+.3f)\n",
914 phi,layer->radius,pos[0],pos[1],pos[2],pt);
915
916 probTr.Print();
917 exit(1);
918 }
919 }
920 // save resolutions at this layer
921 fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters
922 fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters
923 //printf(">> L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);
924 // End save
925 //
926 if (isVertex) continue;
927 //
928 // create fake measurement with the errors assigned to the layer
929 // account for the measurement there
930 double meas[2] = {probTr.GetY(),probTr.GetZ()};
931 double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};
932 //
933
934 if (!probTr.Update(meas,measErr2)) {
935 printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",
936 meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);
937 probTr.Print();
938 exit(1);
939 }
940
941 // correct for materials of this layer
942 // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param
943 if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {
944 printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);
945 probTr.Print();
946 exit(1);
947 }
948 //
949 }
950
951 // Pattern recognition is done .... save values like vertex resolution etc.
952
953 // Convert the Convariance matrix parameters into physical quantities
954 // The results are propogated to the previous point but *do not* include the measurement at that point.
955 // deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/probTr.Get1P(); // Absolute magnitude so ignore charge
956 deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/TMath::Abs(probTr.GetSigned1Pt());
957 fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ); // results in percent
958 fResolutionRPhi[i] = TMath::Sqrt( probTr.GetSigmaY2() ) * 1.e4; // result in microns
959 fResolutionZ[i] = TMath::Sqrt( probTr.GetSigmaZ2() ) * 1.e4; // result in microns
960 // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius
961
962
963 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
964 printf("Number of active layers: %d\n",fNumberOfActiveLayers) ;
965 if (fAtLeastCorr != -1) printf("Number of combinatorics for probabilities: %d\n",komb);
966 printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000);
967 printf("Name Radius Thickness PointResOn PointResOnZ DetRes DetResZ Density Efficiency\n") ;
968 // printOnce =0;
969 }
970
971 // print out and efficiency calculation
972 Int_t iLayActive=0;
973 for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop
974
975 layer = (CylLayerK*)fLayers.At(j);
976
977 // Convert to Meters, Tesla, and GeV
978 Float_t radius = layer->radius /100;
979 Float_t phiRes = layer->phiRes /100;
980 Float_t zRes = layer->zRes /100;
981 Float_t radLength = layer->radL;
982 Float_t leff = layer->eff; // basic layer efficiency
983 Bool_t isDead = layer->isDead;
984
985
986 if ( (!isDead && radLength >0) ) {
987
988 Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
989 phiRes * phiRes ) * 100. ; // work in cm
990 Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
991 zRes * zRes ) * 100. ; // work in cm
992
993
994 Double_t layerEfficiency = 0;
995 if ( EfficiencySearchFlag == 0 )
996 layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;
997 else if ( EfficiencySearchFlag == 1 )
998 layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;
999 else if ( EfficiencySearchFlag == 2 )
1000 layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;
1001
1002 TString name(layer->GetName());
1003 if (!name.Contains("tpc")) {
1004 probLay(2,iLayActive)= layerEfficiency ; // Pcorr
1005 probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull
1006 probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake
1007 iLayActive++;
1008 }
1009 if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
1010
1011 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1012
1013
1014 printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
1015 layer->GetName(), radius*100, radLength,
1016 fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
1017 phiRes*1.e6, zRes*1.e6,
1018 HitDensity(radius*100)) ;
1019 if (!name.Contains("tpc"))
1020 printf("%10.3f\n", layerEfficiency);
1021 else
1022 printf(" - \n");
1023 }
1024
1025 if (!name.Contains("tpc") && !name.Contains("trd")) fEfficiency[massloop][i] *= layerEfficiency;
1026
1027
1028 }
1029
1030 if (fAtLeastCorr != -1) {
1031 // Calculate probabilities from Kombinatorics tree ...
1032 Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);
1033 fEfficiency[massloop][i] = probs[0]; // efficiency
1034 fFake[massloop][i] = probs[1]; // fake
1035 }
1036
1037 /*
1038 // vertex print
1039 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) {
1040 printf("%s:\t ----- ----- %10.0f %11.0f \n", layer->GetName(),fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6);
1041 }
1042 */
1043 }
1044 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1045 if (fNumberOfActiveLayers >=1500) printOnce = 0 ;
1046 printf("\n") ;
1047 }
1048
1049
1050
1051
1052 if (fNumberOfActiveLayers <1500 ) {
1053
1054
1055
1056 // BACKWORD TRACKING +++++++++++++++++
1057 // number of layers is quite low ... efficiency calculation was probably nonsense
1058 // Tracking outward (backword) to get reliable efficiencies from "smoothed estimates"
1059
1060 // For below, see paper, NIM A262 (1987) p.444, eqs.12.
1061 // Equivalently, one can simply combine the forward and backward estimates. Assuming
1062 // pf,Cf and pb,Cb as extrapolated position estimates and errors from fwd and bwd passes one can
1063 // use a weighted estimate Cw = (Cf^-1 + Cb^-1)^-1, pw = Cw (pf Cf^-1 + pb Cb^-1).
1064 // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything.
1065
1066 Bool_t doLikeAliRoot = 0; // don't do the "combined info" but do like in Aliroot
1067
1068 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1069 printf("- Numbers of active layer is low (%d):\n -> \"outward\" fitting done as well to get reliable eff.estimates\n",
1070 fNumberOfActiveLayers);
1071 }
1072
1073 // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam)
1074 // mIstar.UnitMatrix(); // start with unity
1075 if (doLikeAliRoot) {
1076 probTr.ResetCovariance(10);
1077 } else {
1078 // cannot do complete reset, set to very large errors
1079 for (int ic=15;ic--;) trCov[ic] = 0.;
1080 trCov[kY2] = trCov[kZ2] = kLargeErr2Coord;
1081 trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
1082 trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
1083 probTr.CheckCovariance();
1084 // cout<<pt<<": "<<kLargeErr2Coord<<" "<<kLargeErr2Dir<<" "<<kLargeErr2PtI*trPars[kPtI]*trPars[kPtI]<<endl;
1085 }
1086 // Clean up and storing of "forward estimates"
1087 Double_t detPointResForw[kMaxNumberOfDetectors][kNptBins], detPointZResForw[kMaxNumberOfDetectors][kNptBins] ;
1088 for (Int_t k=0; k<kMaxNumberOfDetectors; k++) {
1089 for (Int_t l=0; l<nPt; l++) {
1090 detPointResForw[k][l] = fDetPointRes[k][l];
1091 if (!doLikeAliRoot) fDetPointRes[k][l] = RIDICULOUS;
1092 detPointZResForw[k][l] = fDetPointZRes[k][l];
1093 if (!doLikeAliRoot) fDetPointZRes[k][l] = RIDICULOUS;
1094 }
1095 }
1096
1097 // find first "active layer" - start tracking at the first active layer
1098 Int_t firstActiveLayer = 0;
1099 for (Int_t j=0; j<fLayers.GetEntries(); j++) {
1100 layer = (CylLayerK*)fLayers.At(j);
1101 if (!(layer->isDead)) { // is alive
1102 firstActiveLayer = j;
1103 break;
1104 }
1105 }
1106 probTr.Rotate(0);
1107 for (Int_t j=firstActiveLayer; j<(fLayers.GetEntries()); j++) { // Layer loop
1108
1109 layer = (CylLayerK*)fLayers.At(j);
1110 // CylLayerK *nextlayer = (CylLayerK*)fLayers.At(j+1);
1111
1112 TString name(layer->GetName());
1113 Bool_t isVertex = name.Contains("vertex");
1114 if (!PropagateToR(&probTr, layer->radius,bGauss,1)) exit(1);
1115 //if (!probTr.PropagateTo(last->radius,bGauss)) exit(1);
1116 if (!isVertex) {
1117 // rotate to frame with X axis normal to the surface
1118 double pos[3];
1119 probTr.GetXYZ(pos); // lab position
1120 double phi = TMath::ATan2(pos[1],pos[0]);
1121 if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
1122 if (!probTr.Rotate(phi)) {
1123 printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f (pt=%+.3f)\n",
1124 phi,layer->radius,pos[0],pos[1],pos[2],pt);
1125
1126 probTr.Print();
1127 exit(1);
1128 }
1129 }
1130 //
1131 fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters
1132 fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters
1133 //
1134 //printf("<< L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);
1135 // create fake measurement with the errors assigned to the layer
1136 // account for the measurement there
1137 if (isVertex) continue;
1138 double meas[2] = {probTr.GetY(),probTr.GetZ()};
1139 double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};
1140 //
1141 if (!probTr.Update(meas,measErr2)) {
1142 printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",
1143 meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);
1144 probTr.Print();
1145 exit(1);
1146 }
1147 // correct for materials of this layer
1148 // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param
1149 if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {
1150 printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);
1151 probTr.Print();
1152 exit(1);
1153 }
1154 }
1155
1156 // values below NOT REALIABLE -> they do not point to the vertex but outwards !!!!!!!
1157 // ++++++++++++++
1158 // also update the values for the track position ??????
1159 /*
1160 // Pattern recognition is done .... save values like vertex resolution etc.
1161
1162 // Invert the Matrix to recover the convariance matrix
1163 mIstar.Invert() ;
1164 // Convert the Convariance matrix parameters into physical quantities
1165 // The results are propogated to the previous point but *do not* include the measurement at that point.
1166 deltaPoverP = TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3 ; // Absolute magnitude so ignore charge
1167 fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ) ; // results in percent
1168 fResolutionRPhi[i] = TMath::Sqrt( mIstar(0,0) ) * 1.e6 ; // result in microns
1169 fResolutionZ[i] = TMath::Sqrt( mIstar(1,1) ) * 1.e6 ; // result in microns
1170 // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius
1171 */
1172
1173 // deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/TMath::Abs(probTr.GetSigned1Pt());
1174 // fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ); // results in percent
1175
1176 // Weighted combination of the forward and backward estimates
1177 if (!doLikeAliRoot) {
1178 for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {
1179 fDetPointRes[j][i] = detPointResForw[j][i]*fDetPointRes[j][i]/TMath::Sqrt((detPointResForw[j][i]*detPointResForw[j][i]) + (fDetPointRes[j][i]*fDetPointRes[j][i]));
1180 fDetPointZRes[j][i] = detPointZResForw[j][i]*fDetPointZRes[j][i]/TMath::Sqrt((detPointZResForw[j][i]*detPointZResForw[j][i]) + (fDetPointZRes[j][i]*fDetPointZRes[j][i]));
1181 }
1182 }
1183 // Set Detector-Efficiency Storage area to unity
1184 fEfficiency[massloop][i] = 1.0 ;
1185
1186 // print out and efficiency calculation
1187 iLayActive=0;
1188 for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop
1189
1190 layer = (CylLayerK*)fLayers.At(j);
1191
1192 // Convert to Meters, Tesla, and GeV
1193 Float_t radius = layer->radius /100;
1194 Float_t phiRes = layer->phiRes /100;
1195 Float_t zRes = layer->zRes /100;
1196 Float_t radLength = layer->radL;
1197 Float_t leff = layer->eff;
1198 Bool_t isDead = layer->isDead;
1199
1200
1201 Double_t layerEfficiency = 0;
1202 if ( (!isDead && radLength >0) ) {
1203 Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
1204 phiRes * phiRes ) * 100. ; // work in cm
1205 Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
1206 zRes * zRes ) * 100. ; // work in cm
1207 if ( EfficiencySearchFlag == 0 )
1208 layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;
1209 else if ( EfficiencySearchFlag == 1 )
1210 layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;
1211 else if ( EfficiencySearchFlag == 2 )
1212 layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;
1213
1214 TString name(layer->GetName());
1215 if (!name.Contains("tpc")) {
1216 probLay(2,iLayActive)= layerEfficiency ; // Pcorr
1217 probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull
1218 probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake
1219 iLayActive++;
1220 }
1221 if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
1222
1223 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1224
1225
1226 printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
1227 layer->GetName(), radius*100, radLength,
1228 fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
1229 phiRes*1.e6, zRes*1.e6,
1230 HitDensity(radius*100)) ;
1231 if (!name.Contains("tpc"))
1232 printf("%10.3f\n", layerEfficiency);
1233 else
1234 printf(" - \n");
1235 }
1236
1237 if (massloop==2 && j==kDetLayer) { // copy layer specific performances
1238 fEfficProlongLay[i] = layerEfficiency;
1239 }
1240
1241 if (!name.Contains("tpc") && !name.Contains("trd")) fEfficiency[massloop][i] *= layerEfficiency;
1242
1243
1244
1245 }
1246 if (fAtLeastCorr != -1) {
1247 // Calculate probabilities from Kombinatorics tree ...
1248 Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);
1249 fEfficiency[massloop][i] = probs[0]; // efficiency
1250 fFake[massloop][i] = probs[1]; // fake
1251 }
1252 }
1253 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1254 printOnce = 0 ;
1255 printf("\n") ;
1256 }
1257 }
1258
1259 if (massloop==2) { // copy layer specific performances
1260 fResolutionRPhiLay[i] = fDetPointRes[kDetLayer][i];
1261 fResolutionZLay[i] = fDetPointZRes[kDetLayer][i];
1262 }
1263
1264 } // pt loop
1265
1266
1267
1268 } // mass loop
1269
1270 probTr.SetUseLogTermMS(kFALSE); // Reset of MS term usage to avoid problems since its static
1271
1272
1273
1274}
1275
1276
1277TGraph * DetectorK::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
1278 //
1279 // returns the momentum resolution
1280 //
1281
1282 TGraph *graph = new TGraph(kNptBins, fTransMomenta, fMomentumRes);
1283 graph->SetTitle("Momentum Resolution .vs. Pt" ) ;
1284 // graph->GetXaxis()->SetRangeUser(0.,5.0) ;
1285 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1286 graph->GetXaxis()->CenterTitle();
1287 graph->GetXaxis()->SetNoExponent(1) ;
1288 graph->GetXaxis()->SetMoreLogLabels(1) ;
1289 graph->GetYaxis()->SetTitle("Momentum Resolution (%)") ;
1290 graph->GetYaxis()->CenterTitle();
1291
1292 graph->SetMaximum(20);
1293 graph->SetMinimum(0.1);
1294 graph->SetLineColor(color);
1295 graph->SetMarkerColor(color);
1296 graph->SetLineWidth(linewidth);
1297
1298 return graph;
1299
1300}
1301
1302TGraph * DetectorK::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {
1303
1304 // Returns the pointing resolution
1305 // axis = 0 ... rphi pointing resolution
1306 // axis = 1 ... z pointing resolution
1307 //
1308
1309 TGraph * graph = 0;
1310
1311 if (axis==0) {
1312 graph = new TGraph ( kNptBins, fTransMomenta, fResolutionRPhi ) ;
1313 graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;
1314 graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
1315 } else {
1316 graph = new TGraph ( kNptBins, fTransMomenta, fResolutionZ ) ;
1317 graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
1318 graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;
1319 }
1320
1321 graph->SetMinimum(1) ;
1322 graph->SetMaximum(1000.1) ;
1323 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1324 graph->GetXaxis()->CenterTitle();
1325 graph->GetXaxis()->SetNoExponent(1) ;
1326 graph->GetXaxis()->SetMoreLogLabels(1) ;
1327 graph->GetYaxis()->CenterTitle();
1328
1329 graph->SetLineWidth(linewidth);
1330 graph->SetLineColor(color);
1331 graph->SetMarkerColor(color);
1332
1333 return graph;
1334
1335}
1336
1337TGraph * DetectorK::GetGraphLayerInfo(Int_t plot, Int_t color, Int_t linewidth) {
1338
1339 // Returns the pointing resolution
1340 // plot = 0 ... rphi pointing resolution
1341 // plot = 1 ... z pointing resolution
1342 // plot = 2 ... prolongation efficiency (outwards)
1343 //
1344
1345
1346 Double_t fDet[kNptBins];
1347 for ( Int_t i = 0 ; i < kNptBins ; i++ ) { // pt loop
1348 if (plot==0)
1349 fDet[i] = fResolutionRPhiLay[i]*1e6; // in microns
1350 else if (plot==1)
1351 fDet[i] = fResolutionZLay[i]*1e6; // in microns
1352 else
1353 fDet[i] = fEfficProlongLay[i]*100; // in percent
1354 }
1355
1356 CylLayerK *l = (CylLayerK*) fLayers.At(kDetLayer);
1357 TGraph * graph = 0;
1358 graph = new TGraph ( kNptBins, fTransMomenta, fDet ) ;
1359 if (plot==0) {
1360 graph->SetTitle(Form("R-#phi Pointing Resolution onto layer \"%s\"",(char*)l->GetName()) );
1361 graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
1362 } else if (plot==1){
1363 graph->SetTitle(Form("Z Pointing Resolution onto layer \"%s\"",(char*)l->GetName()) ) ;
1364 graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;
1365 } else {
1366 graph->SetTitle(Form("Prolongation efficiency onto layer \"%s\"",(char*)l->GetName()) ) ;
1367 graph->GetYaxis()->SetTitle("Prolongation efficiency (%)") ;
1368 graph->SetMinimum(0);
1369 graph->SetMaximum(100);
1370 }
1371
1372 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1373 graph->GetXaxis()->CenterTitle();
1374 graph->GetXaxis()->SetNoExponent(1) ;
1375 graph->GetXaxis()->SetMoreLogLabels(1) ;
1376 graph->GetYaxis()->CenterTitle();
1377
1378 graph->SetLineWidth(linewidth);
1379 graph->SetLineColor(color);
1380 graph->SetMarkerColor(color);
1381
1382 return graph;
1383
1384}
1385
1386
1387
1388TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {
1389 //
1390 // returns the Pointing resolution (accoring to Telescope equation)
1391 // axis =0 ... in rphi
1392 // axis =1 ... in z
1393 //
1394
1395 Double_t resolution[kNptBins];
1396
1397 Double_t layerResolution[2];
1398 Double_t layerRadius[2];
1399 Double_t layerThickness[2];
1400
1401 Int_t count =0; // search two first active layers
1402 printf("Telescope equation for layers: ");
1403 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
1404 CylLayerK *l = (CylLayerK*)fLayers.At(i);
1405 if (!l->isDead && l->radius>0) {
1406 layerRadius[count] = l->radius;
1407 layerThickness[count] = l->radL;
1408 if (axis==0) {
1409 layerResolution[count] = l->phiRes;
1410 } else {
1411 layerResolution[count] = l->zRes;
1412 }
1413 printf("%s, ",l->GetName());
1414 count++;
1415 }
1416 if (count>=2) break;
1417 }
1418 printf("\n");
1419
1420 Double_t pt, momentum, thickness,aMCS ;
1421 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1422
1423 for ( Int_t i = 0 ; i < kNptBins ; i++ ) {
1424 // Reference data as if first two layers were acting all alone
1425 pt = fTransMomenta[i] ;
1426 momentum = pt / TMath::Cos(lambda) ; // Total momentum
1427 resolution[i] = layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1]
1428 + layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;
1429 resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;
1430 thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;
1431 aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;
1432 resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;
1433 resolution[i] = TMath::Sqrt(resolution[i]) * 10000.0 ; // result in microns
1434 }
1435
1436
1437
1438 TGraph* graph = new TGraph ( kNptBins, fTransMomenta, resolution ) ;
1439
1440 if (axis==0) {
1441 graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;
1442 graph->GetYaxis()->SetTitle("RPhi Pointing Resolution (#mum) ") ;
1443 } else {
1444 graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
1445 graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum) ") ;
1446 }
1447 graph->SetMinimum(1) ;
1448 graph->SetMaximum(300.1) ;
1449 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1450 graph->GetXaxis()->CenterTitle();
1451 graph->GetXaxis()->SetNoExponent(1) ;
1452 graph->GetXaxis()->SetMoreLogLabels(1) ;
1453 graph->GetYaxis()->CenterTitle();
1454
1455 graph->SetLineColor(color);
1456 graph->SetMarkerColor(color);
1457 graph->SetLineStyle(kDashed);
1458 graph->SetLineWidth(linewidth);
1459
1460 return graph;
1461
1462}
1463
1464TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {
1465 //
1466 // particle = 0 ... choosen particle (setted particleMass)
1467 // particle = 1 ... Pion
1468 // particle = 2 ... Kaon
1469 // particle = 3 ... D0
1470 //
1471 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1472
1473 Double_t particleEfficiency[kNptBins]; // with chosen particle mass
1474 Double_t kaonEfficiency[kNptBins], pionEfficiency[kNptBins], d0efficiency[kNptBins];
1475 Double_t partEfficiency[2][400];
1476
1477 if (particle != 0) {
1478 // resulting Pion and Kaon efficiency scaled with overall efficiency
1479 Double_t doNotDecayFactor;
1480 for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
1481
1482 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1483 // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
1484 Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
1485 if ( massloop == 1 ) { // KAON
1486 doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
1487 kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
1488 } else { // PION
1489 doNotDecayFactor = 1.0 ;
1490 pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
1491 }
1492 partEfficiency[0][j] = pionEfficiency[j];
1493 partEfficiency[1][j] = kaonEfficiency[j];
1494 }
1495 }
1496
1497 // resulting estimate of the D0 efficiency
1498 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1499 d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);
1500 }
1501 } else {
1502 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1503 particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;
1504 // NOTE: Decay factor (see kaon) should be included to be realiable
1505 }
1506 }
1507
1508 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1509 pionEfficiency[j] *= 100;
1510 kaonEfficiency[j] *= 100;
1511 d0efficiency[j] *= 100;
1512 particleEfficiency[j] *= 100;
1513 }
1514
1515 TGraph * graph = 0;
1516 if (particle==0) {
1517 graph = new TGraph ( kNptBins, fTransMomenta, particleEfficiency ) ; // choosen mass
1518 graph->SetLineWidth(1);
1519 } else if (particle==1) {
1520 graph = new TGraph ( kNptBins, fTransMomenta, pionEfficiency ) ;
1521 graph->SetLineWidth(1);
1522 } else if (particle ==2) {
1523 graph = new TGraph ( kNptBins, fTransMomenta, kaonEfficiency ) ;
1524 graph->SetLineWidth(1);
1525 } else if (particle ==3) {
1526 graph = new TGraph ( kNptBins, fTransMomenta, d0efficiency ) ;
1527 graph->SetLineStyle(kDashed);
1528 } else
1529 return 0;
1530
1531 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1532 graph->GetXaxis()->CenterTitle();
1533 graph->GetXaxis()->SetNoExponent(1) ;
1534 graph->GetXaxis()->SetMoreLogLabels(1) ;
1535 graph->GetYaxis()->SetTitle("Efficiency (%)") ;
1536 graph->GetYaxis()->CenterTitle();
1537
1538 graph->SetMinimum(0.01) ;
1539 graph->SetMaximum(100) ;
1540
1541 graph->SetLineColor(color);
1542 graph->SetMarkerColor(color);
1543 graph->SetLineWidth(linewidth);
1544
1545 return graph;
1546}
1547
1548TGraph * DetectorK::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) {
1549 //
1550 // particle = 0 ... choosen particle (setted particleMass)
1551 // particle = 1 ... Pion
1552 // particle = 2 ... Kaon
1553 //
1554
1555 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1556
1557 Double_t particleFake[kNptBins]; // with chosen particle mass
1558 Double_t kaonFake[kNptBins], pionFake[kNptBins];
1559 Double_t partFake[2][kNptBins];
1560
1561 if (particle != 0) {
1562 // resulting Pion and Kaon efficiency scaled with overall efficiency
1563 Double_t doNotDecayFactor;
1564 for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
1565
1566 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1567 // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
1568 Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
1569 if ( massloop == 1 ) { // KAON
1570 doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
1571 kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;
1572 } else { // PION
1573 pionFake[j] = fFake[0][j] ;
1574 }
1575 partFake[0][j] = pionFake[j];
1576 partFake[1][j] = kaonFake[j];
1577 }
1578 }
1579
1580 } else {
1581 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1582 particleFake[j] = fFake[2][j];
1583 // NOTE: Decay factor (see kaon) should be included to be realiable
1584 }
1585 }
1586
1587 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1588 pionFake[j] *= 100;
1589 kaonFake[j] *= 100;
1590 particleFake[j] *= 100;
1591 }
1592
1593 TGraph * graph = 0;
1594 if (particle==0) {
1595 graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass
1596 graph->SetLineWidth(1);
1597 } else if (particle==1) {
1598 graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;
1599 graph->SetLineWidth(1);
1600 } else if (particle ==2) {
1601 graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;
1602 graph->SetLineWidth(1);
1603 }
1604
1605 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1606 graph->GetXaxis()->CenterTitle();
1607 graph->GetXaxis()->SetNoExponent(1) ;
1608 graph->GetXaxis()->SetMoreLogLabels(1) ;
1609 graph->GetYaxis()->SetTitle("Fake (%)") ;
1610 graph->GetYaxis()->CenterTitle();
1611
1612 graph->SetMinimum(0.01) ;
1613 graph->SetMaximum(100) ;
1614
1615 graph->SetLineColor(color);
1616 graph->SetMarkerColor(color);
1617 graph->SetLineWidth(linewidth);
1618
1619 return graph;
1620}
1621TGraph * DetectorK::GetGraphRecoPurity(Int_t particle,Int_t color, Int_t linewidth) {
1622 //
1623 // particle = 0 ... choosen particle (setted particleMass)
1624 // particle = 1 ... Pion
1625 // particle = 2 ... Kaon
1626 //
1627
1628 // Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1629
1630 Double_t particleFake[kNptBins]; // with chosen particle mass
1631 Double_t kaonFake[kNptBins], pionFake[kNptBins];
1632 // Double_t partFake[2][kNptBins];
1633
1634 if (particle != 0) {
1635 cout <<" not implemented"<<endl;
1636
1637 } else {
1638 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1639 particleFake[j] = fFake[2][j];
1640 // NOTE: Decay factor (see kaon) should be included to be realiable
1641 }
1642 }
1643
1644 // Get Purity
1645 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
1646 pionFake[j] = (1-pionFake[j])*100;
1647 kaonFake[j] = (1-kaonFake[j])*100;
1648 particleFake[j] = (1-particleFake[j])*100;
1649 }
1650
1651 TGraph * graph = 0;
1652 if (particle==0) {
1653 graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass
1654 graph->SetLineWidth(1);
1655 } else if (particle==1) {
1656 graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;
1657 graph->SetLineWidth(1);
1658 } else if (particle ==2) {
1659 graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;
1660 graph->SetLineWidth(1);
1661 }
1662
1663 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1664 graph->GetXaxis()->CenterTitle();
1665 graph->GetXaxis()->SetNoExponent(1) ;
1666 graph->GetXaxis()->SetMoreLogLabels(1) ;
1667 graph->GetYaxis()->SetTitle("Purity (%)") ;
1668 graph->GetYaxis()->CenterTitle();
1669
1670 graph->SetMinimum(0.01) ;
1671 graph->SetMaximum(100) ;
1672
1673 graph->SetLineColor(color);
1674 graph->SetMarkerColor(color);
1675 graph->SetLineWidth(linewidth);
1676
1677 return graph;
1678}
1679
1680
1681TGraph* DetectorK::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {
1682 //
1683 // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)
1684 // mode 0: impact parameter (convolution of pointing and vertex resolution)
1685 // mode 1: pointing resolution
1686 // mode 2: vtx resolution
1687
1688
1689 TGraph *graph = new TGraph();
1690
1691 // TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?
1692 TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); //
1693 TFormula vtxResZ("vtxResZ","600/(x+6)+10"); //
1694
1695 TGraph *trackRes = GetGraphPointingResolution(axis,1);
1696 Double_t *pt = trackRes->GetX();
1697 Double_t *trRes = trackRes->GetY();
1698 for (Int_t ip =0; ip<trackRes->GetN(); ip++) {
1699 Double_t vtxRes = 0;
1700 if (axis==0)
1701 vtxRes = vtxResRPhi.Eval(pt[ip]);
1702 else
1703 vtxRes = vtxResZ.Eval(pt[ip]);
1704
1705 if (mode==0)
1706 graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));
1707 else if (mode ==1)
1708 graph->SetPoint(ip,pt[ip],trRes[ip]);
1709 else
1710 graph->SetPoint(ip,pt[ip],vtxRes);
1711 }
1712
1713 graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;
1714 graph->GetYaxis()->SetTitle("d_{0} r#phi resolution (#mum)") ;
1715
1716 graph->SetMinimum(1) ;
1717 graph->SetMaximum(300.1) ;
1718 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1719 graph->GetXaxis()->CenterTitle();
1720 graph->GetXaxis()->SetNoExponent(1) ;
1721 graph->GetXaxis()->SetMoreLogLabels(1) ;
1722 graph->GetYaxis()->CenterTitle();
1723
1724 graph->SetLineColor(color);
1725 graph->SetMarkerColor(color);
1726 graph->SetLineWidth(linewidth);
1727
1728 return graph;
1729
1730}
1731
1732TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
1733 //
1734 // returns graph according to the number
1735 //
1736 switch(number) {
1737 case 1:
1738 return GetGraphPointingResolution(0,color, linewidth); // dr
1739 case 2:
1740 return GetGraphPointingResolution(1,color, linewidth); // dz
1741 case 3:
1742 return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele
1743 case 4:
1744 return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele
1745 case 5:
1746 return GetGraphMomentumResolution(color, linewidth); // pt resolution
1747 case 10:
1748 return GetGraphRecoEfficiency(0, color, linewidth); // tracked particle
1749 case 11:
1750 return GetGraphRecoEfficiency(1, color, linewidth); // eff. pion
1751 case 12:
1752 return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon
1753 case 13:
1754 return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0
1755 case 15:
1756 return GetGraphRecoFakes(0, color, linewidth); // Fake tracked particle
1757 case 16:
1758 return GetGraphRecoFakes(1, color, linewidth); // Fake pion
1759 case 17:
1760 return GetGraphRecoFakes(2, color, linewidth); // Fake kaon
1761 default:
1762 printf(" Error: chosen graph number not valid\n");
1763 }
1764 return 0;
1765
1766}
1767
1768void DetectorK::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon) {
1769
1770 // All New configuration with X0 = 0.3 and resolution = 4 microns
1771
1772 AddLayer((char*)"bpipe",2.0,0.0022); // beam pipe
1773 AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation
1774
1775 // new ideal Pixel properties?
1776 Double_t x0 = 0.0050;
1777 Double_t resRPhi = 0.0006;
1778 Double_t resZ = 0.0006;
1779
1780 if (flagMon) {
1781 x0 = 0.0030;
1782 resRPhi = 0.0004;
1783 resZ = 0.0004;
1784 }
1785
1786 AddLayer((char*)"ddd1", 2.2 , x0, resRPhi, resZ);
1787 AddLayer((char*)"ddd2", 2.8 , x0, resRPhi, resZ);
1788 AddLayer((char*)"ddd3", 3.6 , x0, resRPhi, resZ);
1789 AddLayer((char*)"ddd4", 20.0 , x0, resRPhi, resZ);
1790 AddLayer((char*)"ddd5", 22.0 , x0, resRPhi, resZ);
1791 AddLayer((char*)"ddd6", 41.0 , x0, resRPhi, resZ);
1792 AddLayer((char*)"ddd7", 43.0 , x0, resRPhi, resZ);
1793
1794 if (flagTPC) {
1795 AddTPC(0.1,0.1); // TPC
1796 }
1797}
1798
1799void DetectorK::MakeAliceCurrent(Int_t AlignResiduals, Bool_t flagTPC) {
1800
1801 // Numbers taken from
1802 // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks
1803 // number for misalingment: private communication with Andrea Dainese
1804
1805 AddLayer((char*)"bpipe",2.94,0.0022); // beam pipe
1806 AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation
1807 AddLayer((char*)"tshld1",11.5,0.0065); // Thermal shield // 1.3% /2
1808 AddLayer((char*)"tshld2",31.0,0.0065); // Thermal shield // 1.3% /2
1809
1810
1811 if (flagTPC) {
1812 AddTPC(0.1,0.1); // TPC
1813 }
1814 // Adding the ITS - current configuration
1815
1816 if (AlignResiduals==0) {
1817
1818 AddLayer((char*)"spd1", 3.9, 0.0114, 0.0012, 0.0130);
1819 AddLayer((char*)"spd2", 7.6, 0.0114, 0.0012, 0.0130);
1820 AddLayer((char*)"sdd1",15.0, 0.0113, 0.0035, 0.0025);
1821 AddLayer((char*)"sdd2",23.9, 0.0126, 0.0035, 0.0025);
1822 AddLayer((char*)"ssd1",38.0, 0.0083, 0.0020, 0.0830);
1823 AddLayer((char*)"ssd2",43.0, 0.0086, 0.0020, 0.0830);
1824
1825 } else if (AlignResiduals==1) {
1826
1827 // tracking errors ...
1828 // (Additional systematic errors due to misalignments) ...
1829 // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
1830 // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
1831
1832 AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010),
1833 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1834 AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),
1835 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1836 AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
1837 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1838 AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
1839 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1840 AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
1841 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1842 AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
1843 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1844
1845 } else if (AlignResiduals==2) {
1846
2c4acef0 1847 // tracking errors ... PLUS ... chip misalignment
84356182 1848
1849 // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
1850 // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
1851
2c4acef0 1852 // the ITS chips are misalignment with small gaussian smearings with
84356182 1853 // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
1854
1855 AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008),
1856 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1857 AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008),
1858 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1859 AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
1860 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1861 AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
1862 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1863 AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
1864 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1865 AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
1866 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1867
1868 } else {
1869
2c4acef0 1870 // the ITS chips are misalignment with small gaussian smearings with
84356182 1871 // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
1872 // unknown in Z ????
1873
1874 AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
1875 TMath::Sqrt(0.0130*0.0130+0.000*0.000));
1876 AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
1877 TMath::Sqrt(0.0130*0.0130+0.000*0.000));
1878 AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
1879 TMath::Sqrt(0.0025*0.0025+0.000*0.000));
1880 AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
1881 TMath::Sqrt(0.0025*0.0025+0.000*0.000));
1882 AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
1883 TMath::Sqrt(0.0830*0.0830+0.000*0.000));
1884 AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
1885 TMath::Sqrt(0.0830*0.0830+0.000*0.000));
1886
1887
1888 }
1889
1890}
1891
1892
1893void DetectorK::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {
1894 //
1895 // Produces the standard performace plots
1896 //
1897 TGraph *eff,*momRes,*pointRes;
1898 if (!add) {
1899
1900 TCanvas *c1 = new TCanvas("c1","c1");//,100,100,500,500);
1901 c1->Divide(2,2);
1902
1903 c1->cd(1); gPad->SetGridx(); gPad->SetGridy();
1904 gPad->SetLogx();
1905 eff = GetGraphRecoEfficiency(1,color,linewidth);
1906 eff->SetName(Form("grEff%d",1));
1907 eff->SetTitle("Efficiencies");
1908 eff->Draw("AL");
1909 eff->SetMaximum(110);
1910 if (!onlyPionEff) {
1911 eff = GetGraphRecoEfficiency(2,color,linewidth);
1912 eff->SetName(Form("grEff%d",2));
1913 eff->Draw("L");
1914 eff = GetGraphRecoEfficiency(3,color,linewidth);
1915 eff->SetName(Form("grEff%d",3));
1916 eff->Draw("L");
1917 }
1918 c1->cd(2); gPad->SetGridx(); gPad->SetGridy();
1919 gPad->SetLogy(); gPad->SetLogx();
1920 momRes = GetGraphMomentumResolution(color,linewidth);
1921 momRes->SetName(Form("grMomRes%d",1));
1922 momRes->Draw("AL");
1923
1924 c1->cd(3); gPad->SetGridx(); gPad->SetGridy();
1925 gPad->SetLogx();
1926 pointRes = GetGraphPointingResolution(0,color,linewidth);
1927 pointRes->SetName(Form("pointRRes%d",0));
1928 pointRes->Draw("AL");
1929 //
1930 c1->cd(4); gPad->SetGridx(); gPad->SetGridy();
1931 gPad->SetLogx();
1932 pointRes = GetGraphPointingResolution(1,color,linewidth);
1933 pointRes->SetName(Form("pointZRes%d",0));
1934 pointRes->Draw("AL");
1935
1936 } else {
1937
1938 TVirtualPad *c1 = gPad->GetMother();
1939
1940 c1->cd(1);
1941 eff = GetGraphRecoEfficiency(1,color,linewidth);
1942 eff->SetName(Form("grEff%dadd",1));
1943 eff->Draw("L");
1944 if (!onlyPionEff) {
1945 eff = GetGraphRecoEfficiency(2,color,linewidth);
1946 eff->SetName(Form("grEff%dadd",2));
1947 eff->Draw("L");
1948 eff = GetGraphRecoEfficiency(3,color,linewidth);
1949 eff->SetName(Form("grEff%dadd",3));
1950 eff->Draw("L");
1951 }
1952 c1->cd(2);
1953 momRes = GetGraphMomentumResolution(color,linewidth);
1954 momRes->SetName(Form("grMomRes%dadd",1));
1955 momRes->Draw("L");
1956
1957 c1->cd(3);
1958 pointRes = GetGraphPointingResolution(0,color,linewidth);
1959 pointRes->SetName(Form("pointRRes%dadd",0));
1960 pointRes->Draw("L");
1961
1962 c1->cd(4);
1963 pointRes = GetGraphPointingResolution(1,color,linewidth);
1964 pointRes->SetName(Form("pointZRes%dadd",0));
1965 pointRes->Draw("L");
1966
1967 }
1968
1969}
1970
1971
1972Bool_t DetectorK::GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, Double_t bz, Int_t dir)
1973{
1974 // Get local X of the track position estimated at the radius lab radius r.
1975 // The track curvature is accounted exactly
1976 //
1977 // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
1978 // 0 - take the intersection closest to the current track position
1979 // >0 - go along the track (increasing fX)
1980 // <0 - go backward (decreasing fX)
1981 //
1982 // special case of R=0
1983 if (r<kAlmost0) {x=0; return kTRUE;}
1984
1985 const double* pars = tr->GetParameter();
1986 const Double_t &fy=pars[0], &sn = pars[2];
1987 //
1988 double fx = tr->GetX();
1989 double crv = tr->GetC(bz);
1990 if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track
1991 if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
1992 double det = (r-fx)*(r+fx);
1993 if (det<0) return kFALSE; // does not reach raduis r
1994 x = fx;
1995 if (dir==0) return kTRUE;
1996 det = TMath::Sqrt(det);
1997 if (dir>0) { // along the track direction
1998 if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle
1999 else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
2000 }
2001 else if(dir>0) { // agains track direction
2002 if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
2003 else if (fy>det) return kFALSE; // track is against Y axis
2004 }
2005 }
2006 else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
2007 double det = (r-fy)*(r+fy);
2008 if (det<0) return kFALSE; // does not reach raduis r
2009 det = TMath::Sqrt(det);
2010 if (!dir) {
2011 x = fx>0 ? det : -det; // choose the solution requiring the smalest step
2012 return kTRUE;
2013 }
2014 else if (dir>0) { // along the track direction
2015 if (fx > det) return kFALSE; // current point is in on the right from the circle
2016 else if (fx <-det) x = -det; // on the left
2017 else x = det; // within the circle
2018 }
2019 else { // against the track direction
2020 if (fx <-det) return kFALSE;
2021 else if (fx > det) x = det;
2022 else x = -det;
2023 }
2024 }
2025 else { // general case of straight line
2026 double cs = TMath::Sqrt((1-sn)*(1+sn));
2027 double xsyc = fx*sn-fy*cs;
2028 double det = (r-xsyc)*(r+xsyc);
2029 if (det<0) return kFALSE; // does not reach raduis r
2030 det = TMath::Sqrt(det);
2031 double xcys = fx*cs+fy*sn;
2032 double t = -xcys;
2033 if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step
2034 else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)
2035 if (t>=-det) t += -det; // take minimal step giving t>0
2036 else return kFALSE; // both solutions have negative t
2037 }
2038 else { // go in increasing fx direction. (t+-det < 0)
2039 if (t<det) t -= det; // take minimal step giving t<0
2040 else return kFALSE; // both solutions have positive t
2041 }
2042 x = fx + cs*t;
2043 }
2044 }
2045 else { // helix
2046 // get center of the track circle
2047 double tR = 1./crv; // track radius (for the moment signed)
2048 double cs = TMath::Sqrt((1-sn)*(1+sn));
2049 double x0 = fx - sn*tR;
2050 double y0 = fy + cs*tR;
2051 double r0 = TMath::Sqrt(x0*x0+y0*y0);
2052 // printf("Xc:%+e Yc:%+e tR:%e r0:%e\n",x0,y0,tR,r0);
2053 //
2054 if (r0<=kAlmost0) return kFALSE; // the track is concentric to circle
2055 tR = TMath::Abs(tR);
2056 double tR2r0 = tR/r0;
2057 double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
2058 double det = (1.-g)*(1.+g);
2059 if (det<0) return kFALSE; // does not reach raduis r
2060 det = TMath::Sqrt(det);
2061 //
2062 // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
2063 // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
2064 // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
2065 //
2066 double tmp = 1.+g*tR2r0;
2067 x = x0*tmp;
2068 double y = y0*tmp;
2069 if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
2070 double dfx = tR2r0*TMath::Abs(y0)*det;
2071 double dfy = tR2r0*x0*TMath::Sign(det,y0);
2072 if (dir==0) { // chose the one which corresponds to smallest step
2073 double delta = (x-fx)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
2074 if (delta<0) x += dfx;
2075 else x -= dfx;
2076 }
2077 else if (dir>0) { // along track direction: x must be > fx
2078 x -= dfx; // try the smallest step (dfx is positive)
2079 if (x<fx && (x+=dfx+dfx)<fx) return kFALSE;
2080 }
2081 else { // backward: x must be < fx
2082 x += dfx; // try the smallest step (dfx is positive)
2083 if (x>fx && (x-=dfx+dfx)>fx) return kFALSE;
2084 }
2085 }
2086 else { // special case: track touching the circle just in 1 point
2087 if ( (dir>0&&x<fx) || (dir<0&&x>fx) ) return kFALSE;
2088 }
2089 }
2090 //
2091 return kTRUE;
2092}
2093
2094
2095
2096Double_t* DetectorK::PrepareEffFakeKombinations(TMatrixD *probKomb, TMatrixD *probLay) {
2097
2098 if (!probLay) {
2099 printf("Error: Layer tracking efficiencies not set \n");
2100 return 0;
2101 }
2102
2103 TMatrixD &tProbKomb = *probKomb;
2104 TMatrixD &tProbLay = *probLay;
2105
2106
2107 // Int_t base = tProbLay.GetNcols(); // 3? null, fake, correct
2108 Int_t nLayer = tProbKomb.GetNcols(); // nlayer? - number of ITS layers
2109 Int_t komb = tProbKomb.GetNrows(); // 3^nlayer? - number of kombinations
2110
2111 // Fill probabilities
2112
2113 Double_t probEff =0;
2114 Double_t probFake =0;
2115 for (Int_t num=0; num<komb; num++) {
2116 Int_t flCorr=0, flFake=0, flNull=0;
2117 for (Int_t l=0; l<nLayer; l++) {
2118 if (tProbKomb(num,l)==0)
2119 flNull++;
2120 else if (tProbKomb(num,l)==1)
2121 flFake++;
2122 else if (tProbKomb(num,l)==2)
2123 flCorr++;
2124 else
2125 printf("Error: unexpected values in combinatorics table\n");
2126 }
2127
2128 Int_t fkAtLeastCorr = fAtLeastCorr;
2129 if (fAtLeastCorr == -1) fkAtLeastCorr = nLayer; // all hits are "correct"
2130
2131 if (flCorr>=fkAtLeastCorr && flFake==0) { // at least correct but zero fake
2132 Double_t probEffLayer = 1;
2133 for (Int_t l=0; l<nLayer; l++) {
2134 probEffLayer *= tProbLay((Int_t)tProbKomb(num,l),l);
2135 // cout<<a(num,l)<<" ";
2136 }
2137 // cout<<endl;
2138 probEff+=probEffLayer;
2139 }
2140
2141 if (flFake>=fAtLeastFake) {
2142 Double_t probFakeLayer = 1;
2143 for (Int_t l=0; l<nLayer; l++) {
2144 probFakeLayer *= tProbLay((Int_t)tProbKomb(num,l),l);
2145 // cout<<a(num,l)<<" ";
2146 }
2147 // cout<<endl;
2148 probFake+=probFakeLayer;
2149 }
2150
2151 }
20273357 2152 Double_t *probs = new Double_t[2];
84356182 2153 probs[0] = probEff; probs[1] = probFake;
2154 return probs;
2155
2156}
2157
2158//____________________________________
2159Bool_t DetectorK::PropagateToR(AliExternalTrackParam* trc, double r, double b, int dir)
2160{
2161 // go to radius R
2162 //
2163 double xR = 0;
2164 double rr = r*r;
2165 int iter = 0;
2166 const double kTiny = 1e-6;
2167 while(1) {
2168 if (!GetXatLabR(trc, r ,xR, b, dir)) {
2169 printf("Track with pt=%f cannot reach radius %f\n",trc->Pt(),r);
2170 trc->Print();
2171 return kFALSE;
2172 }
2173 if (!trc->PropagateTo(xR, b)) {printf("Failed to propagate to X=%f for R=%f\n",xR,r); trc->Print(); return kFALSE;}
2174 double rcurr2 = xR*xR + trc->GetY()*trc->GetY();
2175 if (TMath::Abs(rcurr2-rr)<kTiny || rr<kTiny) return kTRUE;
2176 //
2177 // two radii correspond to this X...
2178 double pos[3]; trc->GetXYZ(pos);
2179 double phi = TMath::ATan2(pos[1],pos[0]); //TMath::ASin( trc->GetSnp() );
2180 if (!trc->Rotate(phi)) {printf("Failed to rotate to %f to propagate to R=%f\n",phi,r); trc->Print(); return kFALSE;}
2181 if (++iter>8) {printf("Failed to propagate to R=%f after %d steps\n",r,iter); trc->Print(); return kFALSE;}
2182 }
2183 return kTRUE;
2184}
2185