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