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