ITS/UPGRADE/DetectorK.cxx

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