ITS/UPGRADE/testITSU/DetectorK.cxx

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