ITS/UPGRADE/DetectorK.cxx

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