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