ITS/UPGRADE/DetectorK.cxx

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