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