ITS/UPGRADE/Detector.cxx

   1 #include "Detector.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 /***********************************************************
  13
  14 Fast Simulation tool for Inner Tracker Systems
  15
  16 original code of using the billoir technique was developed
  17 for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov
  18 http://rnc.lbl.gov/~jhthomas
  19
  20 Changes by S. Rossegger -> see header file
  21
  22 ***********************************************************/
  23
  24
  25 #define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector
  26
  27 #define Luminosity    1.e27       // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )
  28 #define SigmaD        6.0         // Size of the interaction diamond (cm) (LHC = 6.0 cm)
  29 #define dNdEtaMinB    950         // Multiplicity per unit Eta  (AuAu MinBias = 170, Central = 700)
  30 #define dNdEtaCent    2300        // Multiplicity per unit Eta  (LHC at 5.5 TeV not known)
  31
  32 #define CrossSectionMinB         8    // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)
  33 #define AcceptanceOfTpcAndSi     1//0.60 //0.35  // Assumed geometric acceptance (efficiency) of the TPC and Si detectors
  34 #define UPCBackgroundMultiplier  1.0   // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)
  35 #define OtherBackground          0.0   // Increase multiplicity in detector (0.0 to 1.0 * minBias)  (eg 0.0)
  36 #define EfficiencySearchFlag     2     // Define search method:
  37                                        // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0.
  38
  39 #define PionMass                 0.139  // Mass of the Pion
  40 #define KaonMass                 0.498  // Mass of the Kaon
  41 #define D0Mass                   1.865  // Mass of the D0
  42
  43
  44
  45 class CylLayer : public TNamed {
  46 public:
  47
  48   CylLayer(char *name) : TNamed(name,name) {}
  49
  50   Float_t GetRadius()   const {return radius;}
  51   Float_t GetRadL()     const {return radL;}
  52   Float_t GetPhiRes()   const {return phiRes;}
  53   Float_t GetZRes()     const {return zRes;}
  54   Float_t GetLayerEff() const {return eff;}
  55
  56   //  void Print() {printf("  r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }
  57   Float_t radius; Float_t radL; Float_t phiRes; Float_t zRes;
  58   Float_t eff;
  59   Bool_t isDead;
  60
  61  ClassDef(CylLayer,1);
  62 };
  63
  64 ClassImp(Detector)
  65 Detector::Detector()
  66   : TNamed("test_detector","detector"),
  67     fNumberOfLayers(0),
  68     fNumberOfActiveLayers(0),
  69     fBField(0.5),
  70     fLhcUPCscale(2.5),
  71     fIntegrationTime(0.02), // in ms
  72     fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence
  73     fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle
  74     fParticleMass(0.140),    // Standard: pion mass
  75     fMaxRadiusSlowDet(10.)
  76 {
  77   //
  78   // default constructor
  79   //
  80   //  fLayers = new TObjArray();
  81
  82 }
  83
  84 Detector::Detector(char *name, char *title)
  85   : TNamed(name,title),
  86     fNumberOfLayers(0),
  87     fNumberOfActiveLayers(0),
  88     fBField(0.5),
  89     fLhcUPCscale(2.5),
  90     fIntegrationTime(0.02),  // in ms
  91     fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence
  92     fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle
  93     fParticleMass(0.140),     // Standard: pion mass
  94     fMaxRadiusSlowDet(10.)
  95 {
  96   //
  97   // default constructor, that set the name and title
  98   //
  99   //  fLayers = new TObjArray();
 100 }
 101 Detector::~Detector() { //
 102   // virtual destructor
 103   //
 104   //  delete fLayers;
 105 }
 106
 107 void Detector::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes, Float_t zRes, Float_t eff) {
 108   //
 109   // Add additional layer to the list of layers (ordered by radius)
 110   //
 111
 112   CylLayer *newLayer = (CylLayer*) fLayers.FindObject(name);
 113
 114   if (!newLayer) {
 115     newLayer = new CylLayer(name);
 116     newLayer->radius = radius;
 117     newLayer->radL = radL;
 118     newLayer->phiRes = phiRes;
 119     newLayer->zRes = zRes;
 120     newLayer->eff = eff;
 121
 122     if (newLayer->zRes==RIDICULOUS && newLayer->zRes==RIDICULOUS)
 123       newLayer->isDead = kTRUE;
 124     else
 125       newLayer->isDead = kFALSE;
 126
 127     if (fLayers.GetEntries()==0)
 128       fLayers.Add(newLayer);
 129     else {
 130
 131       for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
 132         CylLayer *l = (CylLayer*)fLayers.At(i);
 133         if (radius<l->radius) {
 134           fLayers.AddBefore(l,newLayer);
 135           break;
 136         }
 137         if (radius>l->radius && (i+1)==fLayers.GetEntries() ) {
 138           // even bigger then last one
 139           fLayers.Add(newLayer);
 140         }
 141       }
 142
 143     }
 144     fNumberOfLayers += 1;
 145     if (!(newLayer->isDead)) fNumberOfActiveLayers += 1;
 146
 147
 148   } else {
 149     printf("Layer with the name %s does already exist\n",name);
 150   }
 151
 152
 153 }
 154
 155 void Detector::KillLayer(char *name) {
 156   //
 157   // Marks layer as dead. Contribution only by Material Budget
 158   //
 159
 160   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 161   if (!tmp)
 162     printf("Layer %s not found - cannot mark as dead\n",name);
 163   else {
 164      tmp->phiRes = 999999;
 165      tmp->zRes = 999999;
 166      if (!(tmp->isDead)) {
 167        tmp->isDead = kTRUE;
 168        fNumberOfActiveLayers -= 1;
 169      }
 170   }
 171 }
 172
 173 void Detector::SetRadius(char *name, Float_t radius) {
 174   //
 175   // Set layer radius [cm]
 176   //
 177
 178   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 179
 180
 181   if (!tmp) {
 182     printf("Layer %s not found - cannot set radius\n",name);
 183   } else {
 184
 185     Float_t tmpRadL  = tmp->radL;
 186     Float_t tmpPhiRes = tmp->phiRes;
 187     Float_t tmpZRes = tmp->zRes;
 188
 189     RemoveLayer(name); // so that the ordering is correct
 190     AddLayer(name,radius,tmpRadL,tmpPhiRes,tmpZRes);
 191   }
 192 }
 193
 194 Float_t Detector::GetRadius(char *name) {
 195   //
 196   // Return layer radius [cm]
 197   //
 198
 199   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 200   if (!tmp)
 201     printf("Layer %s not found - cannot get radius\n",name);
 202   else
 203     return tmp->radius;
 204
 205   return 0;
 206 }
 207
 208 void Detector::SetRadiationLength(char *name, Float_t radL) {
 209   //
 210   // Set layer material [cm]
 211   //
 212
 213   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 214   if (!tmp)
 215     printf("Layer %s not found - cannot set layer material\n",name);
 216   else {
 217     tmp->radL = radL;
 218   }
 219 }
 220
 221 Float_t Detector::GetRadiationLength(char *name) {
 222   //
 223   // Return layer radius [cm]
 224   //
 225
 226   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 227   if (!tmp)
 228     printf("Layer %s not found - cannot get layer material\n",name);
 229   else
 230     return tmp->radL;
 231
 232   return 0;
 233
 234 }
 235
 236 void Detector::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
 237   //
 238   // Set layer resolution in [cm]
 239   //
 240
 241   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 242   if (!tmp)
 243     printf("Layer %s not found - cannot set resolution\n",name);
 244   else {
 245
 246     Bool_t wasDead = tmp->isDead;
 247
 248     tmp->phiRes = phiRes;
 249     tmp->zRes = zRes;
 250
 251     if (zRes==RIDICULOUS && phiRes==RIDICULOUS) {
 252       tmp->isDead = kTRUE;
 253       if (!wasDead) fNumberOfActiveLayers -= 1;
 254     } else {
 255       tmp->isDead = kFALSE;
 256       if (wasDead) fNumberOfActiveLayers += 1;
 257     }
 258
 259
 260   }
 261 }
 262
 263 Float_t Detector::GetResolution(char *name, Int_t axis) {
 264   //
 265   // Return layer resolution in [cm]
 266   // axis = 0: resolution in rphi
 267   // axis = 1: resolution in z
 268   //
 269
 270   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 271   if (!tmp)
 272     printf("Layer %s not found - cannot get resolution\n",name);
 273   else {
 274     if (axis==0) return tmp->phiRes;
 275     if (axis==1) return tmp->zRes;
 276     printf("error: axis must be either 0 or 1 (rphi or z axis)\n");
 277   }
 278   return 0;
 279 }
 280
 281 void Detector::SetLayerEfficiency(char *name, Float_t eff) {
 282   //
 283   // Set layer efficnecy (prop that his is missed within this layer)
 284   //
 285
 286   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 287   if (!tmp)
 288     printf("Layer %s not found - cannot set layer efficiency\n",name);
 289   else {
 290     tmp->eff = eff;
 291   }
 292 }
 293
 294 Float_t Detector::GetLayerEfficiency(char *name) {
 295   //
 296   // Get layer efficnecy (prop that his is missed within this layer)
 297   //
 298
 299   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 300   if (!tmp)
 301     printf("Layer %s not found - cannot get layer efficneicy\n",name);
 302   else
 303     return tmp->eff;
 304
 305   return 0;
 306
 307 }
 308
 309 void Detector::RemoveLayer(char *name) {
 310   //
 311   // Removes a layer from the list
 312   //
 313
 314   CylLayer *tmp = (CylLayer*) fLayers.FindObject(name);
 315   if (!tmp)
 316     printf("Layer %s not found - cannot remove it\n",name);
 317   else {
 318     Bool_t wasDead = tmp->isDead;
 319     fLayers.Remove(tmp);
 320     fNumberOfLayers -= 1;
 321     if (!wasDead) fNumberOfActiveLayers -= 1;
 322   }
 323 }
 324
 325
 326 void Detector::PrintLayout() {
 327   //
 328   // Prints the detector layout
 329   //
 330
 331   printf("Detector %s: \"%s\"\n",GetName(),GetTitle());
 332
 333   if (fLayers.GetEntries()>0)
 334     printf("  Name \t\t r [cm] \t  X0 \t  phi & z res [um]\n");
 335
 336   CylLayer *tmp = 0;
 337   for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
 338     tmp = (CylLayer*)fLayers.At(i);
 339
 340     // don't print all the tpc layers
 341     TString name(tmp->GetName());
 342     if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
 343
 344     printf("%d. %s \t %03.2f   \t%1.4f\t  ",i,
 345            tmp->GetName(), tmp->radius, tmp->radL);
 346     if (tmp->phiRes==RIDICULOUS)
 347       printf("  -  ");
 348     else
 349       printf("%3.0f   ",tmp->phiRes*10000);
 350     if (tmp->zRes==RIDICULOUS)
 351       printf("  -\n");
 352     else
 353       printf("%3.0f\n",tmp->zRes*10000);
 354   }
 355 }
 356
 357 void Detector::PlotLayout(Int_t plotDead) {
 358   //
 359   // Plots the detector layout in Front view
 360   //
 361
 362   Double_t x0=0, y0=0;
 363
 364   TGraphErrors *gr = new TGraphErrors();
 365   gr->SetPoint(0,0,0);
 366   CylLayer *lastLayer = (CylLayer*)fLayers.At(fLayers.GetEntries()-1);  Double_t maxRad = lastLayer->radius;
 367   gr->SetPointError(0,maxRad,maxRad);
 368   gr->Draw("APE");
 369
 370
 371   CylLayer *tmp = 0;
 372   for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {
 373     tmp = (CylLayer*)fLayers.At(i);
 374
 375
 376     Double_t txtpos = tmp->radius;
 377     if ((tmp->isDead)) txtpos*=-1; //
 378     TText *txt = new TText(x0,txtpos,tmp->GetName());
 379     txt->SetTextSizePixels(5); txt->SetTextAlign(21);
 380     if (!tmp->isDead || plotDead) txt->Draw();
 381
 382     TEllipse *layEl = new TEllipse(x0,y0,tmp->radius);
 383     //  layEl->SetFillColor(5);
 384     layEl->SetFillStyle(5001);
 385     layEl->SetLineStyle(tmp->isDead+1); // dashed if not active
 386     layEl->SetLineColor(4);
 387     TString name(tmp->GetName());
 388     if (!tmp->isDead) layEl->SetLineWidth(2);
 389     if (name.Contains("tpc") )  layEl->SetLineColor(29);
 390
 391     if (!tmp->isDead || plotDead) layEl->Draw();
 392
 393   }
 394
 395 }
 396
 397
 398
 399 void Detector::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {
 400   //
 401   // Emulates the TPC
 402   //
 403   // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow
 404
 405
 406   AddLayer((char*)"IFC",   77.8,0.01367); // Inner Field cage
 407
 408   // % Radiation Lengths ... Average per TPC row  (i.e. total/159 )
 409   Float_t radLPerRow = 0.000036;
 410
 411   Float_t tpcInnerRadialPitch  =    0.75 ;    // cm
 412   Float_t tpcMiddleRadialPitch =    1.0  ;    // cm
 413   Float_t tpcOuterRadialPitch  =    1.5  ;    // cm
 414   //  Float_t tpcInnerPadWidth     =    0.4  ;    // cm
 415   //  Float_t tpcMiddlePadWidth    =    0.6   ;   // cm
 416   //  Float_t tpcOuterPadWidth     =    0.6   ;   // cm
 417   Float_t innerRows            =   63 ;
 418   Float_t middleRows           =   64  ;
 419   Float_t outerRows            =   32  ;
 420   Float_t tpcRows            =   (innerRows + middleRows + outerRows) ;
 421   Float_t rowOneRadius         =   85.2  ;    // cm
 422   Float_t row64Radius          =  135.1  ;    // cm
 423   Float_t row128Radius         =  199.2  ;    // cm
 424
 425   for ( Int_t k = 0 ; k < tpcRows ; k++ ) {
 426
 427     Float_t rowRadius =0;
 428     if (k<innerRows)
 429       rowRadius =  rowOneRadius + k*tpcInnerRadialPitch ;
 430     else if ( k>=innerRows && k<(innerRows+middleRows) )
 431       rowRadius =  row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;
 432     else if (k>=(innerRows+middleRows) && k<tpcRows )
 433       rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;
 434
 435     if ( k%skip == 0 )
 436       AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow,phiResMean,zResMean);
 437     else
 438       AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow); // non "active" row
 439
 440
 441   }
 442
 443 }
 444
 445 void Detector::RemoveTPC() {
 446
 447   // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)
 448   CylLayer *tmp = 0;
 449   for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
 450     tmp = (CylLayer*)fLayers.At(i);
 451     TString name(tmp->GetName());
 452     if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }
 453   }
 454   RemoveLayer((char*)"IFC");
 455
 456 }
 457
 458
 459 Double_t Detector::ThetaMCS ( Double_t mass, Double_t radLength, Double_t momentum ) const
 460 {
 461   //
 462   // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)
 463   //
 464
 465   Double_t beta  =  momentum / TMath::Sqrt(momentum*momentum+mass*mass)  ;
 466   Double_t theta =  0.0 ;    // Momentum and mass in GeV
 467   // if ( RadLength > 0 ) theta  =  0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );
 468   if ( radLength > 0 ) theta  =  0.0136 * TMath::Sqrt(radLength) / ( beta * momentum ) * (1+0.038*TMath::Log(radLength)) ;
 469   return (theta) ;
 470 }
 471
 472
 473 Double_t Detector::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
 474 {
 475   // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm
 476   // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm
 477   // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius
 478   Double_t sx, sy, goodHit ;
 479   sx = 2 * TMath::Pi() *  searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;
 480   sy = 2 * TMath::Pi() *  searchRadiusZ    * searchRadiusZ    * HitDensity(radius) ;
 481   goodHit =  TMath::Sqrt(1./((1+sx)*(1+sy)))  ;
 482   return ( goodHit ) ;
 483 }
 484
 485
 486 Double_t Detector::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
 487 {
 488   // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm
 489   // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
 490   Double_t sx, goodHit ;
 491   sx = 2 * TMath::Pi() *  searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;
 492   goodHit =  1./(1+sx) ;
 493   return ( goodHit ) ;
 494 }
 495
 496 Double_t Detector::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
 497 {
 498   // Based on work by Ruben Shahoyen
 499   // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
 500   // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account
 501   // Following is correct for 2 DOF
 502
 503   Double_t gamma = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
 504   Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
 505   Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*gamma));
 506
 507   return ( goodHit ) ;
 508 }
 509
 510
 511 Double_t Detector::HitDensity ( Double_t radius )
 512 {
 513   // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.
 514   // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results
 515   // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.
 516   // Note that this function assumes we are working in CM and CM**2 [not meters].
 517   // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.
 518
 519   //  Double_t MaxRadiusSlowDet = 0.1; //?   // Maximum radius for slow detectors.  Fast detectors
 520                                         // and only fast detectors reside outside this radius.
 521   Double_t arealDensity = 0 ;
 522
 523   if ( radius > fMaxRadiusSlowDet )
 524     {
 525       arealDensity  = OneEventHitDensity(dNdEtaCent,radius)  ; // Fast detectors see central collision density (only)
 526       arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius)  ;  // Increase density due to background
 527     }
 528
 529   if (radius < fMaxRadiusSlowDet )
 530     { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.
 531       arealDensity  = OneEventHitDensity(dNdEtaCent,radius)
 532                     + IntegratedHitDensity(dNdEtaMinB,radius)
 533                     + UpcHitDensity(radius) ;
 534       arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;
 535       // Increase density due to background
 536     }
 537
 538   return ( arealDensity ) ;
 539 }
 540
 541
 542 double Detector::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const
 543 {
 544   // This is for one event at the vertex.  No smearing.
 545   double den   = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?
 546   // note: surface of sphere is  '4*pi*r^2'
 547   //       surface of cylinder is '2*pi*r* h'
 548   return den ;
 549 }
 550
 551
 552 double Detector::IntegratedHitDensity(Double_t multiplicity, Double_t radius)
 553 {
 554   // The integral of minBias events smeared over a gaussian vertex distribution.
 555   // Based on work by Yan Lu 12/20/2006, all radii in centimeters.
 556
 557   Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;
 558   Double_t den   = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;
 559
 560   // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.
 561   if ( den < OneEventHitDensity(multiplicity,radius) )  den = OneEventHitDensity(multiplicity,radius) ;
 562
 563   return den ;
 564 }
 565
 566
 567 double Detector::UpcHitDensity(Double_t radius)
 568 {
 569   // QED electrons ...
 570
 571   Double_t mUPCelectrons ;                                 ;
 572   mUPCelectrons =  fLhcUPCscale * (1.23 - radius/6.5)      ;  // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC
 573   if ( mUPCelectrons < 0 ) mUPCelectrons =  0.0             ;  // UPC electrons fall off quickly and don't go to large R
 574   mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ;  // UPCs increase Mulitiplicty ~ proportional to MinBias rate
 575   mUPCelectrons *= UPCBackgroundMultiplier                 ;  // Allow for an external multiplier (eg 0-1) to turn off UPC
 576
 577   return mUPCelectrons ;
 578 }
 579
 580
 581 double Detector::Dist(double z, double r)
 582 {
 583   // Convolute dEta/dZ  distribution with assumed Gaussian of vertex z distribution
 584   // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm
 585   // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.
 586   Int_t    index  =  1     ;     // Start weight at 1 for Simpsons rule integration
 587   Int_t    nsteps =  301   ;     // NSteps must be odd for Simpson's rule to work
 588   double   dist   =  0.0   ;
 589   double   dz0    =  ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1)  ;  //cm
 590   double    z0    =  0.0   ;     //cm
 591   for(int i=0; i<nsteps; i++){
 592     if ( i == nsteps-1 ) index = 1 ;
 593     z0 = -4*SigmaD + i*dz0 ;
 594     dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) *
 595       (1/sqrt((z-z0)*(z-z0) + r*r)) ;
 596     if ( index != 4 ) index = 4; else index = 2 ;
 597   }
 598   return dist;
 599 }
 600
 601 #define  PZero   0.861  // Momentum of back to back decay particles in the CM frame
 602 #define  EPiZero 0.872  // Energy of the pion from a D0 decay at rest
 603 #define  EKZero  0.993  // Energy of the Kaon from a D0 decay at rest
 604
 605 Double_t Detector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const {
 606   // Math from Ron Longacre.  Note hardwired energy to bin conversion for PtK and PtPi.
 607
 608   Double_t const1  =  pt / D0Mass ;
 609   Double_t const2  =  TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;
 610   Double_t sum, ptPi, ptK ;
 611   Double_t effp, effk ;
 612
 613   sum = 0.0 ;
 614   for ( Int_t k = 0 ; k < 360 ; k++ )   {
 615
 616     Double_t theta = k * TMath::Pi() / 180. ;
 617
 618     ptPi = TMath::Sqrt(
 619                        PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
 620                        const1*const1*EPiZero*EPiZero -
 621                        2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +
 622                        PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
 623                        ) ;
 624
 625     ptK = TMath::Sqrt(
 626                       PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
 627                       const1*const1*EKZero*EKZero +
 628                       2*PZero*TMath::Cos(theta)*const2*const1*EKZero +
 629                       PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
 630                       ) ;
 631
 632     // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays
 633     Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
 634     Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
 635
 636     if ( pionindex >= 400 ) pionindex = 399 ;
 637     if ( pionindex >= 0 )   effp = corrEfficiency[0][pionindex] ;
 638     if ( pionindex <  0 )   effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd
 639     if ( effp < 0 )         effp = 0 ;
 640
 641     if ( kaonindex >= 400 ) kaonindex = 399 ;
 642     if ( kaonindex >= 0 )   effk = corrEfficiency[1][kaonindex] ;
 643     if ( kaonindex <  0 )   effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd
 644     if ( effk < 0 )         effk = 0 ;
 645
 646     // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.
 647
 648     sum += effp * effk ;
 649
 650   }
 651
 652   Double_t mean =sum/360;
 653   return mean ;
 654
 655 }
 656
 657
 658 void Detector::SolveDOFminusOneAverage() {
 659   //
 660   // Short study to address "# layers-1 efficiencies"
 661   // saves the means in the according arrays
 662   // Note: Obviously, does not work for the Telescope equation
 663   //
 664
 665   Double_t fMomentumResM[400], fResolutionRPhiM[400], fResolutionZM[400];
 666   Double_t efficiencyM[3][400];
 667   for (Int_t i=0; i<400; i++) {
 668     fMomentumResM[i] = 0;   // Momentum resolution
 669     fResolutionRPhiM[i] = 0;   // Resolution in R
 670     fResolutionZM[i] = 0; // Resolution in Z
 671     for (Int_t part=0; part<3; part++)
 672       efficiencyM[part][i] = 0; // efficiencies
 673   }
 674
 675   // loop over active layers in ITS (remove 1 by 1)
 676   Int_t nITSLayers = 0;
 677   CylLayer *layer =0;
 678   for (Int_t j=0; j<(fLayers.GetEntries()-1); j++) {
 679     layer = (CylLayer*)fLayers.At(j);
 680     TString name(layer->GetName());
 681     if (name.Contains("tpc")) continue;
 682     if (!(layer->isDead))  {
 683
 684       nITSLayers++;
 685       printf("Kill Layer %s\n",name.Data());
 686       Double_t rRes = GetResolution((char*)name.Data(),0);
 687       Double_t zRes = GetResolution((char*)name.Data(),1);
 688       KillLayer((char*)name.Data());
 689       //   PrintLayout();
 690       SolveViaBilloir(1,0);
 691
 692       // produce sum for the mean calculation
 693       for (Int_t i=0; i<400; i++) {
 694         fMomentumResM[i] += fMomentumRes[i];   // Momentum resolution
 695         fResolutionRPhiM[i] += fResolutionRPhi[i];   // Resolution in R
 696         fResolutionZM[i] += fResolutionZ[i]; // Resolution in Z
 697         for (Int_t part=0; part<3; part++)
 698           efficiencyM[part][i] += fEfficiency[part][i]; // efficiencies
 699       }
 700
 701       // "Restore" layer ...
 702       SetResolution((char*)name.Data(),rRes,zRes);
 703
 704     }
 705   }
 706
 707   // save means in "std. Arrays"
 708   for (Int_t i=0; i<400; i++) {
 709     fMomentumRes[i] = fMomentumResM[i]/nITSLayers;   // Momentum resolution
 710     fResolutionRPhi[i] = fResolutionRPhiM[i]/nITSLayers;   // Resolution in R
 711     fResolutionZ[i] = fResolutionZM[i]/nITSLayers; // Resolution in Z
 712     for (Int_t part=0; part<3; part++)
 713       fEfficiency[part][i] = efficiencyM[part][i]/nITSLayers; // efficiencies
 714   }
 715
 716
 717 }
 718
 719 void Detector::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t meanPt) {
 720   //
 721   // Solves the current geometry with the Billoir technique
 722   // ( see P. Billoir, Nucl. Instr. and Meth. 225 (1984), p. 352. )
 723   //
 724
 725   Int_t nPt = 400;
 726   // Clean up ......
 727   for (Int_t i=0; i<kMaxNumberOfDetectors; i++) {
 728     for (Int_t j=0; j<nPt; j++) {
 729       fDetPointRes[i][j]  = RIDICULOUS;
 730       fDetPointZRes[i][j] = RIDICULOUS;
 731       fTransMomenta[i] =0;
 732       fMomentumRes[i] =0;
 733       fResolutionRPhi[i] =0;
 734     }
 735   }
 736
 737   if (!allPt) { // not the whole pt range -> allows a faster minimization at a defined 'meanpt'
 738     nPt = 3;
 739   }
 740
 741
 742
 743   // Calculate track parameters using Billoirs method of matrices
 744
 745   Double_t pt, pz, lambda, momentum, rho, deltaPoverP  ;
 746   Double_t layerThickness, aMCS, mmm, charge ;
 747   Double_t vF1, vF2, vF3, vG1, vG2, vG3 ;
 748   Double_t dx, sigma2, sigma2Z ;
 749   Double_t mass[3] ;
 750   Int_t printOnce = 1 ;
 751
 752   mass[0] = PionMass ; mass[1] = KaonMass ;  // Loop twice for the D0;  first pi then k
 753
 754   mass[2] = fParticleMass;  // third loop
 755
 756   Int_t mStart =0;
 757   if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode
 758
 759   for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ )  {
 760
 761     // PseudoRapidity OK, used as an angle
 762     lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity))  ;
 763
 764     // Track moves along the x axis and deviations are in the y and z directions.
 765     Double_t y, z, a, b, d, e ;
 766     // a = py/px = tan phi, b = pz/px = tan lambda / cos phi, d = 0.3*Ze/Abs(p)
 767
 768     for ( Int_t i = 0 ; i < nPt ; i++ ) { // pt loop
 769
 770       CylLayer *last = (CylLayer*) fLayers.At((fLayers.GetEntries()-1));
 771
 772       // Convert to Meters, Tesla, and GeV
 773
 774       // Starting values based on radius of outermost layer ... log10 steps to ~20 GeV
 775       Double_t bigRad = last->radius/100 /2 ;   // min. pt which the algorithm below could handle
 776       if (bigRad<0.61) bigRad=0.61; // -> min pt around 100 MeV for Bz=0.5T (don't overdo it ... ;-) )
 777       fTransMomenta[i] =  ( 0.3*bigRad*TMath::Abs(fBField) ) - 0.08 + TMath::Power(10,2.3*i/nPt) / 10.0 ;
 778       if (!allPt) { // just 3 points around meanPt
 779         fTransMomenta[i] = meanPt-0.1+i*0.1;
 780       }
 781
 782
 783       // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis
 784       // These are the EndPoint values for y, z, a, b, and d
 785       pt  =  fTransMomenta[i]                       ;  // GeV/c
 786       rho =  pt / TMath::Abs( 0.3 * fBField );  // Radius of curvature of the track (meters)
 787       momentum = pt / TMath::Cos(lambda)     ;  // Total momentum
 788       charge   = -1                          ;  // Assume an electron
 789       pz  =  pt * TMath::Tan(lambda)         ;  // GeV/c
 790       d   =  0.3 * charge / momentum         ;  // Its an electrons so q = -1, which makes d a negative number
 791
 792       //      printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d);
 793
 794       // Matrices for Billoir
 795       TMatrixD bigA(2,2); TMatrixD littleA(2,2);
 796       TMatrixD mIstar(5,5); TMatrixD mD(5,5); TMatrixD mDT(5,5); TMatrixD mA(5,5); TMatrixD mM(5,5); TMatrixD mWork(5,5) ;
 797       TMatrixD mCovStart(5,5);
 798       bigA.Zero();littleA.Zero();mD.Zero();mDT.Zero();mA.Zero();mM.Zero();mWork.Zero();
 799
 800       // Set Detector-Efficiency Storage area to unity
 801       fEfficiency[massloop][i] = 1.0 ;
 802       // Back-propagate the covariance matrix along the track.  I is the inverse of the covariance matrix.
 803       // Start with small variations to a straight line 'information' matrix.  Tuning this matrix is
 804       // equivalent to 'starting the recursion' in Billoirs paper.  It is a tricky business.  You need test
 805       // data and many layers for the matrices to stabilize.  Do not believe results for a small number of layers.
 806       // In our case, always include the TPC and this matrix will be well conditioned before it gets to the Si.
 807
 808
 809       mCovStart.UnitMatrix(); // start with unity (-> sigma estimate of 1cm)
 810
 811       // track parametrization: (y,z,a,b,1/p) with a=py/px=tan(phi) and b = pz/px = tan(lambda)/cos(phi)
 812       mCovStart(0,0) = (2*2);     // 2cm variance in y
 813       mCovStart(1,1) = (2*2);     // 2cm variance in z
 814       mCovStart(2,2) = (0.1*0.1); // approx 0.1rad in angles phi and lambda
 815       mCovStart(3,3) = (1*1);     // approx 0.1rad in angle phi and lambda
 816       mCovStart(4,4) = (1*1);     // 1 GeV error in pt
 817
 818
 819       mIstar = mCovStart; // start with above covariance Matrix
 820       mIstar.Invert();
 821
 822       // mIstar.UnitMatrix(); // start with unity (-> sigma estimate of ~ 1cm)
 823
 824       // 'A' is the "angle" matrix for aMCS at each step.
 825       // 'D' is the "distance" matrix at each step.  It propagates the particle backwards along the track.
 826       // 'M' is the "measurement" matrix. It represents the measurement resolution at each step.
 827
 828       CylLayer *layer = 0;
 829       CylLayer *nextlayer = 0;
 830
 831       // find last "active layer" - start tracking at the last active layer
 832       Int_t lastActiveLayer = 0;
 833       for (Int_t j=(fLayers.GetEntries()-1); j>0; j--) {
 834         layer = (CylLayer*)fLayers.At(j);
 835         if (!(layer->isDead)) { // is alive
 836           lastActiveLayer = j;
 837           break;
 838         }
 839       }
 840
 841       for (Int_t j=lastActiveLayer; j>0; j--) {  // Layer loop
 842
 843         layer = (CylLayer*)fLayers.At(j);
 844         nextlayer = (CylLayer*)fLayers.At(j-1);
 845
 846         // Convert to Meters, Tesla, and GeV
 847         Float_t radius = layer->radius/100;
 848         Float_t radLength = layer->radL;
 849
 850         Float_t nextRadius = nextlayer->radius /100;
 851         Float_t nextPhiResolution = nextlayer->phiRes /100;
 852         Float_t nextZResolution = nextlayer->zRes /100;
 853
 854         //      if ( radius >= rho*TMath::Cos(TMath::Pi()/4) )    // (meters) protect against too low pt
 855         //        { // printf("pt lower bound is too low for this algorithm to take this layer into account (r=%lf)\n",radius);
 856         //          continue ; }
 857
 858         // Jims approximation of the "ideal position", breaks down when rho~radius because of the hyperbola approximation
 859         //      y   =  rho - TMath::Sqrt(rho*rho-radius*radius)  ; // These are 'ideal' locations and
 860         //      a   =  radius / ( rho - y )                      ; // not propagated values which should
 861         //      b   =  rho * TMath::Tan(lambda) / ( rho - y )    ; // be done if we had data. But we don't.
 862         //      z   =  rho * TMath::Tan(lambda) * TMath::ATan(a) ;
 863
 864         // Helix format -> exact intersection at the layer radii
 865         y=0;z=0; // not used
 866         a = radius*radius/TMath::Sqrt(4*radius*radius*rho*rho - TMath::Power(radius,4)); // = Tan phi = y_c/x_c (in cartesian coordinates)
 867         b = TMath::Tan(lambda)/TMath::Cos(TMath::ATan(a));
 868
 869         e   =  TMath::Sqrt( 1 + a*a + b*b )              ;
 870
 871         layerThickness = radLength / TMath::Sin(TMath::Pi()/2 - lambda)  ; // X0 in lambda direction
 872
 873         aMCS =  ThetaMCS(mass[massloop], layerThickness, momentum)           ;
 874         mmm =  ( 1 + a*a + b*b ) ;
 875         bigA(0,0) = mmm*(1+a*a) ; bigA(0,1) = mmm*a*b     ; littleA(0,0) = mIstar(2,2)  ; littleA(0,1) = mIstar(2,3) ;
 876         bigA(1,0) = mmm*a*b     ; bigA(1,1) = mmm*(1+b*b) ; littleA(1,0) = mIstar(3,2)  ; littleA(1,1) = mIstar(3,3) ;
 877         littleA = bigA.Invert() + aMCS*aMCS*littleA ;
 878         littleA = littleA.Invert() ;
 879         mA(0,0) = 0.0 ;  mA(0,1) = 0.0  ;  mA(0,2) = 0.0          ;  mA(0,3) = 0.0          ;  mA(0,4) = 0.0  ;
 880         mA(1,0) = 0.0 ;  mA(1,1) = 0.0  ;  mA(1,2) = 0.0          ;  mA(1,3) = 0.0          ;  mA(1,4) = 0.0  ;
 881         mA(2,0) = 0.0 ;  mA(2,1) = 0.0  ;  mA(2,2) = littleA(0,0) ;  mA(2,3) = littleA(0,1) ;  mA(2,4) = 0.0  ;
 882         mA(3,0) = 0.0 ;  mA(3,1) = 0.0  ;  mA(3,2) = littleA(1,0) ;  mA(3,3) = littleA(1,1) ;  mA(3,4) = 0.0  ;
 883         mA(4,0) = 0.0 ;  mA(4,1) = 0.0  ;  mA(4,2) = 0.0          ;  mA(4,3) = 0.0          ;  mA(4,4) = 0.0  ;
 884         mIstar = mIstar - aMCS*aMCS*mIstar*mA*mIstar ;
 885         dx     = ( radius - nextRadius )         ;        // Billoir works with absolute magnitude of distance
 886         vF3     = e * ( -1 * ( 1 + a*a ) * fBField )   ;        // Assume fBField is on Z axis, only.
 887         vF1     = d * ( a*vF3/(e*e) - 2*e*a*fBField )   ;
 888         vF2     = d * ( b*vF3/(e*e) )                  ;
 889         vG3     = e * ( -1 * a * b * fBField )         ;
 890         vG1     = d * ( a*vG3/(e*e) + e*b*fBField )     ;
 891         vG2     = d * ( b*vG3/(e*e) - e*a*fBField )     ;
 892         mD(0,0) = 1.0 ; mD(0,1) = 0.0 ; mD(0,2) = -1*dx+vF1*dx*dx/2 ; mD(0,3) = vF2*dx*dx/2  ;  mD(0,4) = vF3*dx*dx/2 ;
 893         mD(1,0) = 0.0 ; mD(1,1) = 1.0 ; mD(1,2) = vG1*dx*dx/2 ;  mD(1,3) = -1*dx+vG2*dx*dx/2 ;  mD(1,4) = vG3*dx*dx/2 ;
 894         mD(2,0) = 0.0 ; mD(2,1) = 0.0 ; mD(2,2) = 1-vF1*dx    ;  mD(2,3) = -1*vF2*dx         ;  mD(2,4) = -1*vF3*dx   ;
 895         mD(3,0) = 0.0 ; mD(3,1) = 0.0 ; mD(3,2) = -1*vG1*dx   ;  mD(3,3) = 1-vG2*dx          ;  mD(3,4) = -1*vG3*dx   ;
 896         mD(4,0) = 0.0 ; mD(4,1) = 0.0 ; mD(4,2) = 0.0        ;  mD(4,3) = 0.0              ;  mD(4,4) = 1.0        ;
 897         mDT.Transpose(mD) ;
 898         mIstar  = mDT.Invert() * mIstar * mD.Invert() ;
 899         // Prepare to save Detector efficiencies
 900         mWork = mIstar  ;          // Working copy of matrix
 901         mWork.Invert() ;          // Invert the Matrix to recover the convariance matrix
 902         fDetPointRes [j-1][i]     =  TMath::Sqrt( mWork(0,0) )  ;     // result in meters
 903         fDetPointZRes[j-1][i]     =  TMath::Sqrt( mWork(1,1) )  ;     // result in meters
 904         // End save
 905
 906         sigma2  = ( nextPhiResolution  * nextPhiResolution  )   ;
 907         sigma2Z = ( nextZResolution    * nextZResolution )   ;
 908         mM(0,0) = 1/sigma2; mM(1,1) = 1/sigma2Z  ; // setting detector resolution matrix
 909         mIstar = mIstar + mM ;
 910
 911       }
 912
 913       // Pattern recognition is done .... save values like vertex resolution etc.
 914
 915       // Invert the Matrix to recover the convariance matrix
 916       mIstar.Invert() ;
 917       // Convert the Convariance matrix parameters into physical quantities
 918       // The results are propogated to the previous point but *do not* include the measurement at that point.
 919       deltaPoverP    =  TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3  ;  // Absolute magnitude so ignore charge
 920       fMomentumRes[i]      =  100.* TMath::Abs( deltaPoverP )             ;  // results in percent
 921       fResolutionRPhi[i]      =  TMath::Sqrt( mIstar(0,0) ) * 1.e6            ;  // result in microns
 922       fResolutionZ[i]     =  TMath::Sqrt( mIstar(1,1) ) * 1.e6            ;  // result in microns
 923       //      equivalent[i]  =  TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i])           ;  // Equivalent circular radius
 924
 925
 926       if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
 927         printf("Number of active layers: %d\n",fNumberOfActiveLayers) ;
 928         printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000);
 929         printf("Name   Radius Thickness PointResOn PointResOnZ  DetRes  DetResZ  Density Efficiency\n") ;
 930         //      printOnce =0;
 931       }
 932
 933       // print out and efficiency calculation
 934       for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {  // Layer loop
 935
 936         layer = (CylLayer*)fLayers.At(j);
 937
 938         // Convert to Meters, Tesla, and GeV
 939         Float_t radius = layer->radius /100;
 940         Float_t phiRes = layer->phiRes /100;
 941         Float_t zRes = layer->zRes /100;
 942         Float_t radLength = layer->radL;
 943         Float_t leff = layer->eff; // basic layer efficiency
 944         Bool_t isDead = layer->isDead;
 945
 946
 947         if ( (!isDead && radLength >0) )  {
 948             Double_t rphiError  =  TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
 949                                                 phiRes * phiRes ) * 100.  ; // work in cm
 950             Double_t zError     =  TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
 951                                                 zRes * zRes ) * 100.  ; // work in cm
 952
 953             Double_t layerEfficiency = 0;
 954             if ( EfficiencySearchFlag == 0 )
 955               layerEfficiency =  ProbGoodHit( radius*100, rphiError , zError  ) ;
 956             else if ( EfficiencySearchFlag == 1 )
 957               layerEfficiency =  ProbGoodChiSqHit( radius*100, rphiError , zError  ) ;
 958             else if ( EfficiencySearchFlag == 2 )
 959               layerEfficiency =  ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError  ) ;
 960
 961             TString name(layer->GetName());
 962             if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
 963
 964             if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
 965
 966
 967               printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
 968                      layer->GetName(), radius*100, radLength,
 969                      fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
 970                      phiRes*1.e6, zRes*1.e6,
 971                      HitDensity(radius*100)) ;
 972               if (!name.Contains("tpc"))
 973                 printf("%10.3f\n", layerEfficiency);
 974               else
 975                 printf("        -  \n");
 976             }
 977
 978             if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency;
 979
 980         }
 981         /*
 982         // vertex print
 983         if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) {
 984         printf("%s:\t -----    ----- %10.0f %11.0f \n", layer->GetName(),fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6);
 985         }
 986         */
 987       }
 988       if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
 989         if (fNumberOfActiveLayers >=15) printOnce = 0 ;
 990         printf("\n")  ;
 991       }
 992
 993
 994
 995
 996       if (fNumberOfActiveLayers <15 ) {
 997
 998
 999
1000         // BACKWORD TRACKING +++++++++++++++++
1001         // number of layers is quite low ... efficiency calculation was probably nonsense
1002         // Tracking outward (backword) to get reliable efficiencies from "smoothed estimates"
1003
1004         // For below, see paper, NIM A262 (1987) p.444, eqs.12.
1005         // Equivalently, one can simply combine the forward and backward estimates. Assuming
1006         // pf,Cf and pb,Cb as extrapolated position estimates and errors from fwd and bwd passes one can
1007         // use a weighted estimate Cw = (Cf^-1 + Cb^-1)^-1,  pw = Cw (pf Cf^-1 + pb Cb^-1).
1008         // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything.
1009
1010         Bool_t doLikeAliRoot = kFALSE; // don't do the "combined info" but do like in Aliroot
1011
1012         if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1013           printf("- Numbers of active layer is low (%d):\n    -> \"outward\" fitting done as well to get reliable eff.estimates\n",
1014                  fNumberOfActiveLayers);
1015         }
1016
1017         // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam)
1018         //      mIstar.UnitMatrix(); // start with unity
1019         if (doLikeAliRoot) {
1020           Double_t factor = 10;
1021           TMatrixD diagElements(5,1);
1022           for (Int_t ip =0; ip<5; ip++) diagElements(ip,0) = mIstar(ip,ip)*factor;
1023           mIstar.Zero();
1024           for (Int_t ip =0; ip<5; ip++) mIstar(ip,ip) = diagElements(ip,0);
1025         } else {
1026           // Complete RESET of the covariance matrix
1027           mIstar = mCovStart;
1028           // starts with above covariance Matrix and uses "combinded resolutions" (all infos), see below
1029           // this is the "honest" approach without biasing the results
1030         }
1031         mIstar.Invert() ; // undo the revert from above
1032
1033         // Clean up and storing of "forward estimates"
1034         Double_t detPointResForw[kMaxNumberOfDetectors][400], detPointZResForw[kMaxNumberOfDetectors][400] ;
1035         for (Int_t k=0; k<kMaxNumberOfDetectors; k++) {
1036           for (Int_t l=0; l<nPt; l++) {
1037             detPointResForw[k][l]  = fDetPointRes[k][l];
1038             if (!doLikeAliRoot) fDetPointRes[k][l]  = RIDICULOUS;
1039             detPointZResForw[k][l] = fDetPointZRes[k][l];
1040             if (!doLikeAliRoot) fDetPointZRes[k][l] = RIDICULOUS;
1041           }
1042         }
1043
1044         // find first "active layer" - start tracking at the first active layer
1045         Int_t firstActiveLayer = 0;
1046         for (Int_t j=0; j<(fLayers.GetEntries()-1); j++) {
1047           layer = (CylLayer*)fLayers.At(j);
1048           if (!(layer->isDead)) { // is alive
1049             firstActiveLayer = j;
1050             break;
1051           }
1052         }
1053
1054         for (Int_t j=firstActiveLayer; j<(fLayers.GetEntries()-1); j++) {  // Layer loop
1055
1056           layer = (CylLayer*)fLayers.At(j);
1057           nextlayer = (CylLayer*)fLayers.At(j+1);
1058
1059           // Convert to Meters, Tesla, and GeV
1060           Float_t radius = layer->radius/100;
1061           Float_t radLength = layer->radL;
1062
1063           Float_t nextRadius = nextlayer->radius /100;
1064           Float_t nextPhiResolution = nextlayer->phiRes /100;
1065           Float_t nextZResolution = nextlayer->zRes /100;
1066
1067           y   =  rho - TMath::Sqrt(rho*rho-radius*radius)  ; // These are 'ideal' locations and
1068           a   =  radius / ( rho - y )                      ; // not propagated values which should
1069           b   =  rho * TMath::Tan(lambda) / ( rho - y )    ; // be done if we had data. But we don't.
1070           z   =  rho * TMath::Tan(lambda) * TMath::ATan(a) ;
1071           e   =  TMath::Sqrt( 1 + a*a + b*b )              ;
1072
1073           layerThickness = radLength / TMath::Sin(TMath::Pi()/2 - lambda)  ;
1074
1075           aMCS =  ThetaMCS(mass[massloop], layerThickness, momentum)           ;
1076           mmm =  ( 1 + a*a + b*b ) ;
1077           bigA(0,0) = mmm*(1+a*a) ; bigA(0,1) = mmm*a*b     ; littleA(0,0) = mIstar(2,2)  ; littleA(0,1) = mIstar(2,3) ;
1078           bigA(1,0) = mmm*a*b     ; bigA(1,1) = mmm*(1+b*b) ; littleA(1,0) = mIstar(3,2)  ; littleA(1,1) = mIstar(3,3) ;
1079           littleA = bigA.Invert() + aMCS*aMCS*littleA ;
1080           littleA = littleA.Invert() ;
1081           mA(0,0) = 0.0 ;  mA(0,1) = 0.0  ;  mA(0,2) = 0.0          ;  mA(0,3) = 0.0          ;  mA(0,4) = 0.0  ;
1082           mA(1,0) = 0.0 ;  mA(1,1) = 0.0  ;  mA(1,2) = 0.0          ;  mA(1,3) = 0.0          ;  mA(1,4) = 0.0  ;
1083           mA(2,0) = 0.0 ;  mA(2,1) = 0.0  ;  mA(2,2) = littleA(0,0) ;  mA(2,3) = littleA(0,1) ;  mA(2,4) = 0.0  ;
1084           mA(3,0) = 0.0 ;  mA(3,1) = 0.0  ;  mA(3,2) = littleA(1,0) ;  mA(3,3) = littleA(1,1) ;  mA(3,4) = 0.0  ;
1085           mA(4,0) = 0.0 ;  mA(4,1) = 0.0  ;  mA(4,2) = 0.0          ;  mA(4,3) = 0.0          ;  mA(4,4) = 0.0  ;
1086           mIstar = mIstar - aMCS*aMCS*mIstar*mA*mIstar ;
1087           dx     = ( radius - nextRadius )         ;        // Billoir works with absolute magnitude of distance
1088           vF3     = e * ( -1 * ( 1 + a*a ) * fBField )   ;        // Assume fBField is on Z axis, only.
1089           vF1     = d * ( a*vF3/(e*e) - 2*e*a*fBField )   ;
1090           vF2     = d * ( b*vF3/(e*e) )                  ;
1091           vG3     = e * ( -1 * a * b * fBField )         ;
1092           vG1     = d * ( a*vG3/(e*e) + e*b*fBField )     ;
1093           vG2     = d * ( b*vG3/(e*e) - e*a*fBField )     ;
1094           mD(0,0) = 1.0 ; mD(0,1) = 0.0 ; mD(0,2) = -1*dx+vF1*dx*dx/2 ; mD(0,3) = vF2*dx*dx/2  ;  mD(0,4) = vF3*dx*dx/2 ;
1095           mD(1,0) = 0.0 ; mD(1,1) = 1.0 ; mD(1,2) = vG1*dx*dx/2 ;  mD(1,3) = -1*dx+vG2*dx*dx/2 ;  mD(1,4) = vG3*dx*dx/2 ;
1096           mD(2,0) = 0.0 ; mD(2,1) = 0.0 ; mD(2,2) = 1-vF1*dx    ;  mD(2,3) = -1*vF2*dx         ;  mD(2,4) = -1*vF3*dx   ;
1097           mD(3,0) = 0.0 ; mD(3,1) = 0.0 ; mD(3,2) = -1*vG1*dx   ;  mD(3,3) = 1-vG2*dx          ;  mD(3,4) = -1*vG3*dx   ;
1098           mD(4,0) = 0.0 ; mD(4,1) = 0.0 ; mD(4,2) = 0.0        ;  mD(4,3) = 0.0              ;  mD(4,4) = 1.0        ;
1099           mDT.Transpose(mD) ;
1100           mIstar  = mDT.Invert() * mIstar * mD.Invert() ;
1101           // Prepare to save Detector efficiencies
1102           mWork = mIstar  ;          // Working copy of matrix
1103           mWork.Invert() ;          // Invert the Matrix to recover the convariance matrix
1104           fDetPointRes [j+1][i]     =  TMath::Sqrt( mWork(0,0) )  ;     // result in meters
1105           fDetPointZRes[j+1][i]     =  TMath::Sqrt( mWork(1,1) )  ;     // result in meters
1106           // End save
1107           sigma2  = ( nextPhiResolution  * nextPhiResolution  )   ;
1108           sigma2Z = ( nextZResolution    * nextZResolution )   ;
1109           mM(0,0) = 1/sigma2; mM(1,1) = 1/sigma2Z  ;
1110           mIstar = mIstar + mM ;
1111         }
1112
1113         // values below NOT REALIABLE -> they do not point to the vertex but outwards !!!!!!!
1114         // ++++++++++++++
1115         // also update the values for the track position ??????
1116         /*
1117         // Pattern recognition is done .... save values like vertex resolution etc.
1118
1119         // Invert the Matrix to recover the convariance matrix
1120         mIstar.Invert() ;
1121         // Convert the Convariance matrix parameters into physical quantities
1122         // The results are propogated to the previous point but *do not* include the measurement at that point.
1123         deltaPoverP    =  TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3  ;  // Absolute magnitude so ignore charge
1124         fMomentumRes[i]      =  100.* TMath::Abs( deltaPoverP )             ;  // results in percent
1125         fResolutionRPhi[i]      =  TMath::Sqrt( mIstar(0,0) ) * 1.e6            ;  // result in microns
1126         fResolutionZ[i]     =  TMath::Sqrt( mIstar(1,1) ) * 1.e6            ;  // result in microns
1127         //      equivalent[i]  =  TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i])           ;  // Equivalent circular radius
1128         */
1129
1130         // Weighted combination of the forward and backward estimates
1131         if (!doLikeAliRoot) {
1132           for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {
1133             fDetPointRes[j][i] = 1/(1/detPointResForw[j][i] + 1/fDetPointRes[j][i]);
1134             fDetPointZRes[j][i] = 1/(1/detPointZResForw[j][i] + 1/fDetPointZRes[j][i]);
1135           }
1136         }
1137         // Set Detector-Efficiency Storage area to unity
1138         fEfficiency[massloop][i] = 1.0 ;
1139
1140         // print out and efficiency calculation
1141         for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {  // Layer loop
1142
1143           layer = (CylLayer*)fLayers.At(j);
1144
1145           // Convert to Meters, Tesla, and GeV
1146           Float_t radius = layer->radius /100;
1147           Float_t phiRes = layer->phiRes /100;
1148           Float_t zRes = layer->zRes /100;
1149           Float_t radLength = layer->radL;
1150           Float_t leff = layer->eff;
1151           Bool_t isDead = layer->isDead;
1152
1153
1154           if ( (!isDead && radLength >0) )  {
1155             Double_t rphiError  =  TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
1156                                                 phiRes * phiRes ) * 100.  ; // work in cm
1157             Double_t zError     =  TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
1158                                                 zRes * zRes ) * 100.  ; // work in cm
1159
1160             Double_t layerEfficiency = 0;
1161             if ( EfficiencySearchFlag == 0 )
1162               layerEfficiency =  ProbGoodHit( radius*100, rphiError , zError  ) ;
1163             else if ( EfficiencySearchFlag == 1 )
1164               layerEfficiency =  ProbGoodChiSqHit( radius*100, rphiError , zError  ) ;
1165             else if ( EfficiencySearchFlag == 2 )
1166               layerEfficiency =  ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError  ) ;
1167
1168             TString name(layer->GetName());
1169             if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
1170
1171             if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1172
1173
1174               printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
1175                      layer->GetName(), radius*100, radLength,
1176                      fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
1177                      phiRes*1.e6, zRes*1.e6,
1178                      HitDensity(radius*100)) ;
1179               if (!name.Contains("tpc"))
1180                 printf("%10.3f\n", layerEfficiency);
1181               else
1182                 printf("        -  \n");
1183             }
1184
1185             if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency;
1186
1187           }
1188         }
1189         if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
1190           printOnce = 0 ;
1191           printf("\n")  ;
1192         }
1193       }
1194     } // mass loop
1195   } // pt loop
1196
1197 }
1198
1199
1200 TGraph * Detector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
1201   //
1202   // returns the momentum resolution
1203   //
1204
1205   TGraph *graph = new TGraph(400, fTransMomenta, fMomentumRes);
1206   graph->SetTitle("Momentum Resolution .vs. Pt" ) ;
1207   //  graph->GetXaxis()->SetRangeUser(0.,5.0) ;
1208   graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ;
1209   graph->GetXaxis()->CenterTitle();
1210   graph->GetXaxis()->SetNoExponent(1) ;
1211   graph->GetXaxis()->SetMoreLogLabels(1) ;
1212   graph->GetYaxis()->SetTitle("Momentum Resolution [%]") ;
1213   graph->GetYaxis()->CenterTitle();
1214
1215   graph->SetMaximum(20) ;
1216   graph->SetMinimum(0.1) ;
1217   graph->SetLineColor(color);
1218   graph->SetMarkerColor(color);
1219   graph->SetLineWidth(linewidth);
1220
1221   return graph;
1222
1223 }
1224
1225 TGraph * Detector::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {
1226
1227   // Returns the pointing resolution
1228   // axis = 0 ... rphi pointing resolution
1229   // axis = 1 ... z pointing resolution
1230   //
1231
1232   TGraph * graph =  0;
1233
1234   if (axis==0) {
1235     graph = new TGraph ( 400, fTransMomenta, fResolutionRPhi ) ;
1236     graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;
1237     graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution [#mum]") ;
1238   } else {
1239     graph =  new TGraph ( 400, fTransMomenta, fResolutionZ ) ;
1240     graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
1241     graph->GetYaxis()->SetTitle("Z Pointing Resolution [#mum]") ;
1242   }
1243
1244   graph->SetMinimum(1) ;
1245   graph->SetMaximum(300.1) ;
1246   graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ;
1247   graph->GetXaxis()->CenterTitle();
1248   graph->GetXaxis()->SetNoExponent(1) ;
1249   graph->GetXaxis()->SetMoreLogLabels(1) ;
1250   graph->GetYaxis()->CenterTitle();
1251
1252   graph->SetLineWidth(linewidth);
1253   graph->SetLineColor(color);
1254   graph->SetMarkerColor(color);
1255
1256   return graph;
1257
1258 }
1259
1260
1261 TGraph * Detector::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {
1262   //
1263   // returns the Pointing resolution (accoring to Telescope equation)
1264   // axis =0 ... in rphi
1265   // axis =1 ... in z
1266   //
1267
1268   Double_t resolution[400];
1269
1270   Double_t layerResolution[2];
1271   Double_t layerRadius[2];
1272   Double_t layerThickness[2];
1273
1274   Int_t count =0; // search two first active layers
1275   printf("Telescope equation for layers:  ");
1276   for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
1277     CylLayer *l = (CylLayer*)fLayers.At(i);
1278     if (!l->isDead && l->radius>0) {
1279       layerRadius[count]     = l->radius;
1280       layerThickness[count]  = l->radL;
1281       if (axis==0) {
1282         layerResolution[count] = l->phiRes;
1283       } else {
1284         layerResolution[count] = l->zRes;
1285       }
1286       printf("%s, ",l->GetName());
1287       count++;
1288     }
1289     if (count>=2) break;
1290   }
1291   printf("\n");
1292
1293   Double_t pt, momentum, thickness,aMCS ;
1294   Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1295
1296   for ( Int_t i = 0 ; i < 400 ; i++ ) {
1297     // Reference data as if first two layers were acting all alone
1298     pt  =  fTransMomenta[i]  ;
1299     momentum = pt / TMath::Cos(lambda)   ;  // Total momentum
1300     resolution[i] =  layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1]
1301       +  layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;
1302     resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;
1303     thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;
1304     aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;
1305     resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;
1306     resolution[i] =  TMath::Sqrt(resolution[i]) * 10000.0 ;  // result in microns
1307   }
1308
1309
1310
1311   TGraph* graph = new TGraph ( 400, fTransMomenta, resolution ) ;
1312
1313   if (axis==0) {
1314     graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;
1315     graph->GetYaxis()->SetTitle("RPhi Pointing Resolution [#mum] ") ;
1316   } else {
1317     graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
1318     graph->GetYaxis()->SetTitle("Z Pointing Resolution [#mum] ") ;
1319   }
1320   graph->SetMinimum(1) ;
1321   graph->SetMaximum(300.1) ;
1322   graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ;
1323   graph->GetXaxis()->CenterTitle();
1324   graph->GetXaxis()->SetNoExponent(1) ;
1325   graph->GetXaxis()->SetMoreLogLabels(1) ;
1326   graph->GetYaxis()->CenterTitle();
1327
1328   graph->SetLineColor(color);
1329   graph->SetMarkerColor(color);
1330   graph->SetLineStyle(kDashed);
1331   graph->SetLineWidth(linewidth);
1332
1333   return graph;
1334
1335 }
1336
1337 TGraph * Detector::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {
1338   //
1339   // particle = 0 ... choosen particle (setted particleMass)
1340   // particle = 1 ... Pion
1341   // particle = 2 ... Kaon
1342   // particle = 3 ... D0
1343   //
1344   Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1345
1346   Double_t particleEfficiency[400]; // with chosen particle mass
1347   Double_t kaonEfficiency[400], pionEfficiency[400], d0efficiency[400];
1348   Double_t partEfficiency[2][400];
1349
1350   if (particle != 0) {
1351     // resulting Pion and Kaon efficiency scaled with overall efficiency
1352     Double_t doNotDecayFactor;
1353     for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
1354
1355       for ( Int_t j = 0 ; j < 400 ; j++ ) {
1356         // JT Test Let the kaon decay.  If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
1357         Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda)           ;  // Total momentum at average rapidity
1358         if ( massloop == 1 ) { // KAON
1359           doNotDecayFactor  = TMath::Exp(-130/(371*momentum/KaonMass)) ;  // Decay length for kaon is 371 cm.
1360           kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
1361         } else { // PION
1362           doNotDecayFactor = 1.0 ;
1363           pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
1364         }
1365         partEfficiency[0][j] = pionEfficiency[j];
1366         partEfficiency[1][j] = kaonEfficiency[j];
1367       }
1368     }
1369
1370     // resulting estimate of the D0 efficiency
1371     for ( Int_t j = 0 ; j < 400 ; j++ ) {
1372       d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);
1373     }
1374   } else {
1375     for ( Int_t j = 0 ; j < 400 ; j++ ) {
1376       particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;
1377       // NOTE: Decay factor (see kaon) should be included to be realiable
1378     }
1379   }
1380
1381   TGraph * graph =  0;
1382   if (particle==0) {
1383     graph = new TGraph ( 400, fTransMomenta, particleEfficiency ) ; // choosen mass
1384     graph->SetLineWidth(1);
1385   }  else if (particle==1) {
1386     graph = new TGraph ( 400, fTransMomenta, pionEfficiency ) ;
1387     graph->SetLineWidth(1);
1388   }  else if (particle ==2) {
1389     graph = new TGraph ( 400, fTransMomenta, kaonEfficiency ) ;
1390     graph->SetLineWidth(1);
1391   }  else if (particle ==3) {
1392     graph = new TGraph ( 400, fTransMomenta, d0efficiency ) ;
1393     graph->SetLineStyle(kDashed);
1394   } else
1395     return 0;
1396
1397   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1398   graph->GetXaxis()->CenterTitle();
1399   graph->GetXaxis()->SetNoExponent(1) ;
1400   graph->GetXaxis()->SetMoreLogLabels(1) ;
1401   graph->GetYaxis()->SetTitle("Efficiency (arbitrary units)") ;
1402   graph->GetYaxis()->CenterTitle();
1403
1404   graph->SetMinimum(0.01) ;
1405   graph->SetMaximum(1.0)  ;
1406
1407   graph->SetLineColor(color);
1408   graph->SetMarkerColor(color);
1409   graph->SetLineWidth(linewidth);
1410
1411   return graph;
1412 }
1413
1414
1415 TGraph* Detector::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {
1416   //
1417   // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)
1418   // mode 0: impact parameter (convolution of pointing and vertex resolution)
1419   // mode 1: pointing resolution
1420   // mode 2: vtx resolution
1421
1422
1423   TGraph *graph = new TGraph();
1424
1425   //  TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?
1426   TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); //
1427   TFormula vtxResZ("vtxResZ","600/(x+6)+10"); //
1428
1429   TGraph *trackRes = GetGraphPointingResolution(axis,1);
1430   Double_t *pt = trackRes->GetX();
1431   Double_t *trRes = trackRes->GetY();
1432   for (Int_t ip =0; ip<trackRes->GetN(); ip++) {
1433     Double_t vtxRes = 0;
1434     if (axis==0)
1435       vtxRes = vtxResRPhi.Eval(pt[ip]);
1436     else
1437       vtxRes = vtxResZ.Eval(pt[ip]);
1438
1439     if (mode==0)
1440       graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));
1441     else if (mode ==1)
1442       graph->SetPoint(ip,pt[ip],trRes[ip]);
1443     else
1444       graph->SetPoint(ip,pt[ip],vtxRes);
1445   }
1446
1447   graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;
1448   graph->GetYaxis()->SetTitle("d_{0} r#phi resolution [#mum]") ;
1449
1450   graph->SetMinimum(1) ;
1451   graph->SetMaximum(300.1) ;
1452   graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ;
1453   graph->GetXaxis()->CenterTitle();
1454   graph->GetXaxis()->SetNoExponent(1) ;
1455   graph->GetXaxis()->SetMoreLogLabels(1) ;
1456   graph->GetYaxis()->CenterTitle();
1457
1458   graph->SetLineColor(color);
1459   graph->SetMarkerColor(color);
1460   graph->SetLineWidth(linewidth);
1461
1462   return graph;
1463
1464 }
1465
1466 TGraph* Detector::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
1467   //
1468   // returns graph according to the number
1469   //
1470   switch(number) {
1471   case 1:
1472     return GetGraphPointingResolution(0,color, linewidth); // dr
1473   case 2:
1474     return GetGraphPointingResolution(1,color, linewidth); // dz
1475   case 3:
1476     return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele
1477   case 4:
1478     return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele
1479   case 5:
1480     return GetGraphMomentumResolution(color, linewidth); // pt resolution
1481   case 10:
1482     return GetGraphRecoEfficiency(0, color, linewidth);  // tracked particle
1483   case 11:
1484     return GetGraphRecoEfficiency(1, color, linewidth);  // eff. pion
1485   case 12:
1486     return GetGraphRecoEfficiency(2, color, linewidth);  // eff. kaon
1487   case 13:
1488     return GetGraphRecoEfficiency(3, color, linewidth);  // eff. D0
1489   default:
1490     printf(" Error: chosen graph number not valid\n");
1491   }
1492   return 0;
1493
1494 }
1495
1496
1497 void Detector::MakeAliceCurrent(Int_t AlignResiduals, Bool_t flagTPC) {
1498
1499   // Numbers taken from
1500   // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks
1501   // number for misalingment: private communication with Andrea Dainese
1502
1503   AddLayer((char*)"bpipe",2.94,0.0022); // beam pipe
1504   AddLayer((char*)"vertex",     0,     0); // dummy vertex for matrix calculation
1505   AddLayer((char*)"tshld1",11.5,0.0065); // Thermal shield  // 1.3% /2
1506   AddLayer((char*)"tshld2",31.0,0.0065); // Thermal shield  // 1.3% /2
1507
1508
1509   if (flagTPC) {
1510     AddTPC(0.1,0.1);                        // TPC
1511   }
1512   // Adding the ITS - current configuration
1513
1514   if (AlignResiduals==0) {
1515
1516     AddLayer((char*)"spd1", 3.9, 0.0114, 0.0012, 0.0130);
1517     AddLayer((char*)"spd2", 7.6, 0.0114, 0.0012, 0.0130);
1518     AddLayer((char*)"sdd1",15.0, 0.0113, 0.0035, 0.0025);
1519     AddLayer((char*)"sdd2",23.9, 0.0126, 0.0035, 0.0025);
1520     AddLayer((char*)"ssd1",38.0, 0.0083, 0.0020, 0.0830);
1521     AddLayer((char*)"ssd2",43.0, 0.0086, 0.0020, 0.0830);
1522
1523   } else if (AlignResiduals==1) {
1524
1525     // tracking errors ...
1526     // (Additional systematic errors due to misalignments) ...
1527     // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020);  // [cm]
1528     // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
1529
1530     AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010),
1531              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1532     AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),
1533              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1534     AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
1535              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1536     AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
1537              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1538     AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
1539              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1540     AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
1541              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1542
1543   } else if (AlignResiduals==2) {
1544
1545     // tracking errors ... PLUS ... module misalignment
1546
1547     // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020);  // [cm]
1548     // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
1549
1550     //  the ITS modules are misalignment with small gaussian smearings with
1551     //  sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
1552
1553     AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008),
1554              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1555     AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008),
1556              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1557     AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
1558              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1559     AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
1560              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1561     AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
1562              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1563     AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
1564              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1565
1566   } else {
1567
1568       //  the ITS modules are misalignment with small gaussian smearings with
1569       //  sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
1570       //  unknown in Z ????
1571
1572     AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
1573              TMath::Sqrt(0.0130*0.0130+0.000*0.000));
1574     AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
1575              TMath::Sqrt(0.0130*0.0130+0.000*0.000));
1576     AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
1577              TMath::Sqrt(0.0025*0.0025+0.000*0.000));
1578     AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
1579              TMath::Sqrt(0.0025*0.0025+0.000*0.000));
1580     AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
1581              TMath::Sqrt(0.0830*0.0830+0.000*0.000));
1582     AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
1583              TMath::Sqrt(0.0830*0.0830+0.000*0.000));
1584
1585
1586   }
1587
1588 }
1589
1590
1591 void Detector::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {
1592   //
1593   // Produces the standard performace plots
1594   //
1595
1596   if (!add) {
1597
1598     TCanvas *c1 = new TCanvas("c1","c1");
1599     c1->Divide(2,2);
1600
1601     c1->cd(1);  gPad->SetGridx();   gPad->SetGridy();
1602     gPad->SetLogx();
1603     TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth);
1604     eff->SetTitle("Efficiencies");
1605     eff->Draw("AC");
1606     if (!onlyPionEff) {
1607       GetGraphRecoEfficiency(2,color,linewidth)->Draw("C");
1608       GetGraphRecoEfficiency(3,color,linewidth)->Draw("C");
1609     }
1610     c1->cd(2); gPad->SetGridx();   gPad->SetGridy();
1611     gPad->SetLogy();  gPad->SetLogx();
1612     GetGraphMomentumResolution(color,linewidth)->Draw("AC");
1613
1614     c1->cd(3); gPad->SetGridx();   gPad->SetGridy();
1615     gPad->SetLogx();
1616     GetGraphPointingResolution(0,color,linewidth)->Draw("AC");
1617
1618     c1->cd(4); gPad->SetGridx();   gPad->SetGridy();
1619     gPad->SetLogx();
1620     GetGraphPointingResolution(1,color,linewidth)->Draw("AC");
1621
1622   } else {
1623
1624     TVirtualPad *c1 = gPad->GetMother();
1625
1626     c1->cd(1);
1627     GetGraphRecoEfficiency(1,color,linewidth)->Draw("C");
1628     if (!onlyPionEff) {
1629       GetGraphRecoEfficiency(2,color,linewidth)->Draw("C");
1630       GetGraphRecoEfficiency(3,color,linewidth)->Draw("C");
1631     }
1632     c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("C");
1633
1634     c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("C");
1635
1636     c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("C");
1637
1638   }
1639
1640 }