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