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