ITS/UPGRADE/KMCDetector.cxx

   1 #include "KMCDetector.h"
   2 #include <TMath.h>
   3 #include <TH1F.h>
   4 #include <TProfile.h>
   5 #include <TF1.h>
   6 #include <TMatrixD.h>
   7 #include <TGraph.h>
   8 #include <TAxis.h>
   9 #include <TFormula.h>
  10 #include <TCanvas.h>
  11 #include <TEllipse.h>
  12 #include <TText.h>
  13 #include <TRandom3.h>
  14 #include <TGraphErrors.h>
  15 #include <TStopwatch.h>
  16
  17 /***********************************************************
  18
  19 Fast Simulation tool for Inner Tracker Systems
  20
  21 ***********************************************************/
  22
  23
  24 #define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector
  25
  26 #define Luminosity    1.e27       // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )
  27 #define SigmaD        6.0         // Size of the interaction diamond (cm) (LHC = 6.0 cm)
  28 #define dNdEtaMinB    1//950//660//950           // Multiplicity per unit Eta  (AuAu MinBias = 170, Central = 700)
  29 // #define dNdEtaCent    2300//15000 //1600//2300        // Multiplicity per unit Eta  (LHC at 5.5 TeV not known)
  30
  31 #define CrossSectionMinB         8    // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)
  32 #define AcceptanceOfTpcAndSi     1 //1//0.60 //0.35  // Assumed geometric acceptance (efficiency) of the TPC and Si detectors
  33 #define UPCBackgroundMultiplier  1.0   // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)
  34 #define OtherBackground          0.0   // Increase multiplicity in detector (0.0 to 1.0 * minBias)  (eg 0.0)
  35 #define EfficiencySearchFlag     2     // Define search method:
  36                                        // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0.
  37
  38 #define PionMass                 0.139  // Mass of the Pion
  39 #define KaonMass                 0.498  // Mass of the Kaon
  40 #define D0Mass                   1.865  // Mass of the D0
  41
  42 //TMatrixD *probKomb; // table for efficiency kombinatorics
  43
  44 ClassImp(KMCProbe)
  45
  46 Int_t    KMCProbe::fgNITSLayers = 0;
  47 Double_t KMCProbe::fgMissingHitPenalty = 2.;
  48 //__________________________________________________________________________
  49 KMCProbe& KMCProbe::operator=(const KMCProbe& src)
  50 {
  51   if (this!=&src) {
  52     AliExternalTrackParam::operator=(src);
  53     fMass = src.fMass;
  54     fChi2 = src.fChi2;
  55     fHits = src.fHits;
  56     fFakes = src.fFakes;
  57     fNHits = src.fNHits;
  58     fNHitsITS = src.fNHitsITS;
  59     fNHitsITSFake = src.fNHitsITSFake;
  60     fInnLrCheck     = src.fInnLrCheck;
  61     for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];
  62   }
  63   return *this;
  64 }
  65
  66 //__________________________________________________________________________
  67 KMCProbe::KMCProbe(KMCProbe& src)
  68   : AliExternalTrackParam(src),
  69     fMass(src.fMass),
  70     fChi2(src.fChi2),
  71     fHits(src.fHits),
  72     fFakes(src.fFakes),
  73     fNHits(src.fNHits),
  74     fNHitsITS(src.fNHitsITS),
  75     fNHitsITSFake(src.fNHitsITSFake),
  76     fInnLrCheck(src.fInnLrCheck)
  77 {
  78   for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];
  79 }
  80
  81 //__________________________________________________________________________
  82 void KMCProbe::ResetCovMat()
  83 {
  84   // reset errors
  85   double *trCov  = (double*)GetCovariance();
  86   double *trPars = (double*)GetParameter();
  87   const double kLargeErr2Coord = 50*50;
  88   const double kLargeErr2Dir = 0.6*0.6;
  89   const double kLargeErr2PtI = 0.5*0.5;
  90   for (int ic=15;ic--;) trCov[ic] = 0.;
  91   trCov[kY2]   = trCov[kZ2]   = kLargeErr2Coord;
  92   trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
  93   trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
  94   //
  95 }
  96
  97 //__________________________________________________________________________
  98 void KMCProbe::Print(Option_t* option) const
  99 {
 100   printf("M=%.3f Chi2=%7.2f (Norm:%6.2f|%d) Hits: Total:%d ITS:%d ITSFakes:%d | Y:%+8.4f Z: %+8.4f |",
 101          fMass,fChi2,GetNormChi2(kTRUE),fInnLrCheck,fNHits,fNHitsITS,fNHitsITSFake, GetY(),GetZ());
 102   for (int i=0;i<fgNITSLayers;i++) {
 103     if (!IsHit(i)) printf(".");
 104     else printf("%c",IsHitFake(i) ? '-':'+');
 105   }
 106   printf("|%s\n",IsKilled() ? " KILLED":"");
 107   TString opt = option;
 108   if (!opt.IsNull()) AliExternalTrackParam::Print(option);
 109 }
 110
 111 //__________________________________________________________________________
 112 Int_t KMCProbe::Compare(const TObject* obj) const
 113 {
 114   // compare to tracks
 115   const KMCProbe* trc = (KMCProbe*) obj;
 116   if (trc->IsKilled()) {
 117     if (IsKilled()) return 0;
 118     return -1;
 119   }
 120   else if (IsKilled()) return 1;
 121   double chi2a = GetNormChi2(kTRUE);
 122   double chi2b = trc->GetNormChi2(kTRUE);
 123   if (chi2a<chi2b) return  -1;
 124   if (chi2a>chi2b) return   1;
 125   return 0;
 126 }
 127
 128 //__________________________________________________________________________
 129 Bool_t KMCProbe::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const
 130 {
 131   // Get local X of the track position estimated at the radius lab radius r.
 132   // The track curvature is accounted exactly
 133   //
 134   // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
 135   // 0  - take the intersection closest to the current track position
 136   // >0 - go along the track (increasing fX)
 137   // <0 - go backward (decreasing fX)
 138   //
 139   // special case of R=0
 140   if (r<kAlmost0) {x=0; return kTRUE;}
 141
 142   const double* pars = GetParameter();
 143   const Double_t &fy=pars[0], &sn = pars[2];
 144   //
 145   double fx = GetX();
 146   double crv = GetC(bz);
 147   if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track
 148     if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
 149       double det = (r-fx)*(r+fx);
 150       if (det<0) return kFALSE;     // does not reach raduis r
 151       x = fx;
 152       if (dir==0) return kTRUE;
 153       det = TMath::Sqrt(det);
 154       if (dir>0) {                       // along the track direction
 155         if (sn>0) {if (fy>det)  return kFALSE;} // track is along Y axis and above the circle
 156         else      {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
 157       }
 158       else if(dir>0) {                                    // agains track direction
 159         if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
 160         else if (fy>det)  return kFALSE;        // track is against Y axis
 161       }
 162     }
 163     else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
 164       double det = (r-fy)*(r+fy);
 165       if (det<0) return kFALSE;     // does not reach raduis r
 166       det = TMath::Sqrt(det);
 167       if (!dir) {
 168         x = fx>0  ? det : -det;    // choose the solution requiring the smalest step
 169         return kTRUE;
 170       }
 171       else if (dir>0) {                    // along the track direction
 172         if      (fx > det) return kFALSE;  // current point is in on the right from the circle
 173         else if (fx <-det) x = -det;       // on the left
 174         else               x =  det;       // within the circle
 175       }
 176       else {                               // against the track direction
 177         if      (fx <-det) return kFALSE;
 178         else if (fx > det) x =  det;
 179         else               x = -det;
 180       }
 181     }
 182     else {                                 // general case of straight line
 183       double cs = TMath::Sqrt((1-sn)*(1+sn));
 184       double xsyc = fx*sn-fy*cs;
 185       double det = (r-xsyc)*(r+xsyc);
 186       if (det<0) return kFALSE;    // does not reach raduis r
 187       det = TMath::Sqrt(det);
 188       double xcys = fx*cs+fy*sn;
 189       double t = -xcys;
 190       if (dir==0) t += t>0 ? -det:det;  // chose the solution requiring the smalest step
 191       else if (dir>0) {                 // go in increasing fX direction. ( t+-det > 0)
 192         if (t>=-det) t += -det;         // take minimal step giving t>0
 193         else return kFALSE;             // both solutions have negative t
 194       }
 195       else {                            // go in increasing fx direction. (t+-det < 0)
 196         if (t<det) t -= det;            // take minimal step giving t<0
 197         else return kFALSE;             // both solutions have positive t
 198       }
 199       x = fx + cs*t;
 200     }
 201   }
 202   else {                                 // helix
 203     // get center of the track circle
 204     double tR = 1./crv;   // track radius (for the moment signed)
 205     double cs = TMath::Sqrt((1-sn)*(1+sn));
 206     double x0 = fx - sn*tR;
 207     double y0 = fy + cs*tR;
 208     double r0 = TMath::Sqrt(x0*x0+y0*y0);
 209     //    printf("Xc:%+e Yc:%+e\n",x0,y0);
 210     //
 211     if (r0<=kAlmost0) {
 212       AliDebug(2,Form("r0 = %f",r0));
 213       return kFALSE;
 214     }            // the track is concentric to circle
 215     tR = TMath::Abs(tR);
 216     double tR2r0 = tR/r0;
 217     double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
 218     double det = (1.-g)*(1.+g);
 219     if (det<0) {
 220       AliDebug(2,Form("g=%f tR=%f r0=%f\n",g,tR, r0));
 221       return kFALSE;
 222     }         // does not reach raduis r
 223     det = TMath::Sqrt(det);
 224     //
 225     // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
 226     // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
 227     // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
 228     //
 229     double tmp = 1.+g*tR2r0;
 230     x = x0*tmp;
 231     double y = y0*tmp;
 232     if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
 233       double dfx = tR2r0*TMath::Abs(y0)*det;
 234       double dfy = tR2r0*x0*TMath::Sign(det,y0);
 235       if (dir==0) {                    // chose the one which corresponds to smallest step
 236         double delta = (x-fx)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
 237         if (delta<0) x += dfx;
 238         else         x -= dfx;
 239       }
 240       else if (dir>0) {  // along track direction: x must be > fx
 241         x -= dfx; // try the smallest step (dfx is positive)
 242         if (x<fx && (x+=dfx+dfx)<fx) return kFALSE;
 243       }
 244       else { // backward: x must be < fx
 245         x += dfx; // try the smallest step (dfx is positive)
 246         if (x>fx && (x-=dfx+dfx)>fx) return kFALSE;
 247       }
 248     }
 249     else { // special case: track touching the circle just in 1 point
 250       if ( (dir>0&&x<fx) || (dir<0&&x>fx) ) return kFALSE;
 251     }
 252   }
 253   //
 254   return kTRUE;
 255 }
 256
 257 //____________________________________
 258 Bool_t KMCProbe::PropagateToR(double r, double b, int dir)
 259 {
 260   // go to radius R
 261   //
 262   double xR = 0;
 263   double rr = r*r;
 264   int iter = 0;
 265   const double kTiny = 1e-4;
 266   while(1) {
 267     if (!GetXatLabR(r ,xR, b, dir)) {
 268       //      printf("Track with pt=%f cannot reach radius %f\n",Pt(),r);
 269       //      Print("l");
 270       return kFALSE;
 271     }
 272
 273     if (!PropagateTo(xR, b)) {
 274       if (AliLog::GetGlobalDebugLevel()>2) {
 275         printf("Failed to propagate to X=%f for R=%f\n",xR,r);
 276         Print("l");
 277       }
 278       return kFALSE;
 279     }
 280     double rcurr2 = xR*xR + GetY()*GetY();
 281     if (TMath::Abs(rcurr2-rr)<kTiny || rr<kAlmost0) return kTRUE;
 282     //
 283     // two radii correspond to this X...
 284     double pos[3]; GetXYZ(pos);
 285     double phi = TMath::ATan2(pos[1],pos[0]);
 286     if (!Rotate(phi)) {
 287       if (AliLog::GetGlobalDebugLevel()>2) {
 288         printf("Failed to rotate to %f to propagate to R=%f\n",phi,r);
 289         Print("l");
 290       }
 291       return kFALSE;
 292     }
 293     if (++iter>8) {
 294       if (AliLog::GetGlobalDebugLevel()>2) {
 295         printf("Failed to propagate to R=%f after %d steps\n",r,iter);
 296         Print("l");
 297       }
 298       return kFALSE;
 299     }
 300   }
 301   return kTRUE;
 302 }
 303
 304
 305 //__________________________________________________________________________
 306 Bool_t KMCProbe::CorrectForMeanMaterial(const KMCLayer* lr, Bool_t inward)
 307 {
 308   //  printf("before at r=%.1f p=%.4f\n",lr->fR, P());
 309   if (AliExternalTrackParam::CorrectForMeanMaterial(lr->fx2X0, inward ? lr->fXRho : -lr->fXRho, GetMass() , kTRUE)) {
 310     //  printf("after  at r=%.1f p=%.4f\n",lr->fR, P());
 311     return kTRUE;
 312   }
 313   AliDebug(2,Form("Failed to apply material correction, X/X0=%.4f", lr->fx2X0));
 314   if (AliLog::GetGlobalDebugLevel()>1) Print();
 315   return kFALSE;
 316 }
 317
 318 /////////////////////////////////////////////////////////////////////////////
 319 ClassImp(KMCCluster)
 320
 321 //_________________________________________________________________________
 322 KMCCluster::KMCCluster(KMCCluster &src)
 323 : TObject(src),
 324   fY(src.fY),fZ(src.fZ),fX(src.fX),fPhi(src.fPhi)
 325 {}
 326
 327 //__________________________________________________________________________
 328 KMCCluster& KMCCluster::operator=(const KMCCluster& src)
 329 {
 330   if (this!=&src) {
 331     TObject::operator=(src);
 332     fY = src.fY;
 333     fZ = src.fZ;
 334     fX = src.fX;
 335     fPhi = src.fPhi;
 336   }
 337   return *this;
 338 }
 339
 340 //_________________________________________________________________________
 341 void KMCCluster::Print(Option_t *) const
 342 {
 343   printf(" Local YZ = (%3.4lf,%3.4lf) | X=%3.4lf  phi: %+.3f %s\n",fY,fZ,fX,fPhi,IsKilled()?"Killed":"");
 344 }
 345
 346 /////////////////////////////////////////////////////////////////////////////
 347 ClassImp(KMCLayer)
 348
 349 Double_t KMCLayer::fgDefEff = 1.0;
 350 //__________________________________________________________________________
 351 KMCLayer::KMCLayer(char *name) :
 352   TNamed(name,name),fR(0),fx2X0(0),fPhiRes(0),fZRes(0),fEff(0),fIsDead(kFALSE),fType(-1),fActiveID(-1),fSig2EstD(999),fSig2EstZ(999),
 353   fClCorr(),fClMC(),fClBg("KMCCluster",5), fTrCorr(), fTrMC("KMCProbe",5)
 354 {
 355   Reset();
 356 }
 357
 358 //__________________________________________________________________________
 359 void KMCLayer::Reset()
 360 {
 361   fTrCorr.Reset();
 362   fClCorr.Reset();
 363   ResetMC();
 364   fSig2EstD = fSig2EstZ = 999;
 365   //
 366 }
 367
 368 //__________________________________________________________________________
 369 KMCProbe* KMCLayer::AddMCTrack(KMCProbe* src)
 370 {
 371   int ntr = GetNMCTracks();
 372   KMCProbe* prb = 0;
 373   if (src) prb = new(fTrMC[ntr]) KMCProbe(*src);
 374   else     prb = new(fTrMC[ntr]) KMCProbe();
 375   if (!IsDead()) prb->ResetHit(GetActiveID());
 376   return prb;
 377 }
 378
 379 //__________________________________________________________________________
 380 void KMCLayer::Print(Option_t *opt) const
 381 {
 382   printf("Lr%3d(A%3d) %10s R=%5.1f X2X0=%.3f XRho=%.3f SigY=%.4f SigZ=%.4f Eff:%4.2f\n",GetUniqueID(),fActiveID,GetName(), fR, fx2X0,fXRho,fPhiRes,fZRes,fEff);
 383   TString opts = opt; opts.ToLower();
 384   if (opts.Contains("c")) {
 385     printf("Clusters: MC: %+7.4f:%+7.4f Ideal: %+7.4f:%+7.4f  NBgCl: %3d NTrMC: %4d\n",fClMC.fY,fClMC.fZ, fClCorr.fY,fClCorr.fZ, GetNBgClusters(),GetNMCTracks());
 386   }
 387 }
 388
 389 /////////////////////////////////////////////////////////////////////////////
 390 Double_t KMCDetector::fgVtxConstraint[2]={-1,-1};
 391
 392 ClassImp(KMCDetector)
 393 KMCDetector::KMCDetector() :
 394 TNamed("test_detector","detector"),
 395   fNLayers(0),
 396   fNActiveLayers(0),
 397   fNActiveITSLayers(0),
 398   fLastActiveLayer(-1),
 399   fLastActiveITSLayer(-1),
 400   fLastActiveLayerTracked(-1),
 401   fBFieldG(5.),
 402   fLhcUPCscale(1.0),
 403   fIntegrationTime(0.02), // in ms
 404   fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence
 405   fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle
 406   fParticleMass(0.140),    // Standard: pion mass
 407   fMaxRadiusSlowDet(10.),
 408   fAtLeastCorr(-1),     // if -1, then correct hit on all ITS layers
 409   fAtLeastFake(1),       // if at least x fakes, track is considered fake ...
 410   fdNdEtaCent(2300),
 411   fMaxChi2Cl(25.),
 412   fMaxNormChi2NDF(5.),
 413   fMinITSHits(4),
 414 //
 415   fMaxChi2ClSQ(4.), // precalulated internally
 416   fMaxSeedToPropagate(300),
 417 //
 418   fUseBackground(kFALSE),
 419   fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),
 420   fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),
 421   fDensFactorEta(1.),
 422 //
 423   fUpdCalls(0),
 424   fHMCLrResidRPhi(0),
 425   fHMCLrResidZ(0),
 426   fHMCLrChi2(0),
 427   fPattITS(0)
 428 {
 429   //
 430   // default constructor
 431   //
 432   //  fLayers = new TObjArray();
 433   RequireMaxChi2Cl(fMaxChi2Cl); // just to precalulate default square
 434 }
 435
 436 KMCDetector::KMCDetector(char *name, char *title)
 437   : TNamed(name,title),
 438     fNLayers(0),
 439     fNActiveLayers(0),
 440     fNActiveITSLayers(0),
 441     fLastActiveLayer(-1),
 442     fLastActiveITSLayer(-1),
 443     fLastActiveLayerTracked(-1),
 444     fBFieldG(5.0),
 445     fLhcUPCscale(1.0),
 446     fIntegrationTime(0.02),  // in ms
 447     fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence
 448     fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle
 449     fParticleMass(0.140),     // Standard: pion mass
 450     fMaxRadiusSlowDet(10.),
 451     fAtLeastCorr(-1),     // if -1, then correct hit on all ITS layers
 452     fAtLeastFake(1),       // if at least x fakes, track is considered fake ...
 453     fdNdEtaCent(2300),
 454     fMaxChi2Cl(9.),
 455     fMaxNormChi2NDF(5.),
 456     fMinITSHits(4),
 457     //
 458     fMaxChi2ClSQ(3.), // precalulated internally
 459     fMaxSeedToPropagate(50),
 460     //
 461     fUseBackground(kFALSE),
 462     fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),
 463     fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),
 464     fDensFactorEta(1.),
 465     //
 466     fUpdCalls(0),
 467     fHMCLrResidRPhi(0),
 468     fHMCLrResidZ(0),
 469     fHMCLrChi2(0),
 470     fPattITS(0)
 471 {
 472   //
 473   // default constructor, that set the name and title
 474   //
 475   //  fLayers = new TObjArray();
 476 }
 477 KMCDetector::~KMCDetector() { //
 478   // virtual destructor
 479   //
 480   //  delete fLayers;
 481 }
 482
 483 void KMCDetector::InitMCWatchHistos()
 484 {
 485   // init utility histos used for MC tuning
 486   enum {kNBinRes=1000};
 487   const double kMaxResidRPhi=1.0,kMaxResidZ=1.0,kMaxChi2=50;
 488   int nlr = fNActiveITSLayers;
 489   TString nam = "mc_residrhi";
 490   TString tit = "MC $Delta Cl-Tr R#phi";
 491   fHMCLrResidRPhi = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidRPhi,nlr,-0.5,nlr-0.5);
 492   fHMCLrResidRPhi->GetXaxis()->SetTitle("cl-tr #Delta r#phi");
 493   fHMCLrResidRPhi->Sumw2();
 494     //
 495   nam = "mc_residz";
 496   tit = "MC $Delta Cl-Tr Z";
 497   fHMCLrResidZ = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidZ,nlr,-0.5,nlr-0.5);
 498   fHMCLrResidZ->GetXaxis()->SetTitle("cl-tr #Delta z");
 499   fHMCLrResidZ->Sumw2();
 500     //
 501   nam = "mc_chi2";
 502   tit = "MC #chi^{2} Cl-Tr Z";
 503   fHMCLrChi2 = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxChi2,nlr,-0.5,nlr-0.5);
 504   fHMCLrChi2->GetXaxis()->SetTitle("cl-tr #chi^{2}");
 505   fHMCLrChi2->Sumw2();
 506   //
 507   SetBit(kUtilHisto);
 508 }
 509
 510
 511 void KMCDetector::AddLayer(char *name, Float_t radius, Float_t x2X0, Float_t xrho, Float_t phiRes, Float_t zRes, Float_t eff) {
 512   //
 513   // Add additional layer to the list of layers (ordered by radius)
 514   //
 515
 516   KMCLayer *newLayer = (KMCLayer*) fLayers.FindObject(name);
 517
 518   if (!newLayer) {
 519     newLayer = new KMCLayer(name);
 520     newLayer->fR = radius;
 521     newLayer->fx2X0 = x2X0;
 522     newLayer->fXRho  = xrho;
 523     newLayer->fPhiRes = phiRes;
 524     newLayer->fZRes = zRes;
 525     eff = TMath::Min(1.f,eff);
 526     newLayer->fEff = eff <0 ? KMCLayer::GetDefEff() : eff;
 527     newLayer->fActiveID = -2;
 528     TString lname = name;
 529     newLayer->fType = KMCLayer::kTypeNA;
 530     if      (lname.Contains("tpc")) newLayer->fType = KMCLayer::kTPC;
 531     else if (lname.Contains("its")) newLayer->fType = KMCLayer::kITS;
 532     if (lname.Contains("vertex")) newLayer->SetBit(KMCLayer::kBitVertex);
 533     //
 534     if (newLayer->fType==KMCLayer::kTypeNA) printf("Attention: the layer %s has undefined type\n",name);
 535     //
 536     newLayer->fIsDead =  (newLayer->fPhiRes==RIDICULOUS && newLayer->fZRes==RIDICULOUS);
 537     //
 538     if (fLayers.GetEntries()==0)
 539       fLayers.Add(newLayer);
 540     else {
 541       //
 542       for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
 543         KMCLayer *l = (KMCLayer*)fLayers.At(i);
 544         if (radius<l->fR) { fLayers.AddBefore(l,newLayer); break; }
 545           if (radius>l->fR && (i+1)==fLayers.GetEntries() ) fLayers.Add(newLayer); // even bigger then last one
 546       }
 547       //
 548     }
 549     //
 550     ClassifyLayers();
 551     //
 552   } else {
 553     printf("Layer with the name %s does already exist\n",name);
 554   }
 555 }
 556
 557 //____________________________________________________________
 558 void KMCDetector::ClassifyLayers()
 559 {
 560   // assign active Id's, etc
 561   fLastActiveLayer = -1;
 562   fLastActiveITSLayer = -1;
 563   fNActiveLayers = 0;
 564   fNActiveITSLayers = 0;
 565   //
 566   int nl = GetNLayers();
 567   for (int il=0;il<nl;il++) {
 568     KMCLayer* lr = GetLayer(il);
 569     lr->SetUniqueID(il);
 570     if (!lr->IsDead()) {
 571       fLastActiveLayer = il;
 572       lr->fActiveID = fNActiveLayers++;
 573       if (lr->IsITS()) {
 574         fLastActiveITSLayer = il;
 575         fNActiveITSLayers++;
 576       }
 577     }
 578   }
 579   //
 580   KMCProbe::SetNITSLayers(fNActiveITSLayers);
 581 }
 582
 583 //____________________________________________________________
 584 void KMCDetector::KillLayer(char *name) {
 585   //
 586   // Marks layer as dead. Contribution only by Material Budget
 587   //
 588
 589   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 590   if (!tmp)
 591     printf("Layer %s not found - cannot mark as dead\n",name);
 592   else {
 593      tmp->fPhiRes = 999999;
 594      tmp->fZRes = 999999;
 595      tmp->fIsDead = kTRUE;
 596      ClassifyLayers();
 597   }
 598 }
 599
 600 void KMCDetector::SetRadius(char *name, Float_t radius) {
 601   //
 602   // Set layer radius [cm]
 603   //
 604   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 605   //
 606   if (!tmp) {
 607     printf("Layer %s not found - cannot set radius\n",name);
 608   } else {
 609
 610     Float_t tmpRadL  = tmp->fx2X0;
 611     Float_t tmpPhiRes = tmp->fPhiRes;
 612     Float_t tmpZRes = tmp->fZRes;
 613     Float_t tmpXRho = tmp->fXRho;
 614     RemoveLayer(name); // so that the ordering is correct
 615     AddLayer(name,radius,tmpRadL,tmpXRho,tmpPhiRes,tmpZRes);
 616   }
 617 }
 618
 619 Float_t KMCDetector::GetRadius(char *name) {
 620   //
 621   // Return layer radius [cm]
 622   //
 623
 624   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 625   if (!tmp)
 626     printf("Layer %s not found - cannot get radius\n",name);
 627   else
 628     return tmp->fR;
 629
 630   return 0;
 631 }
 632
 633 void KMCDetector::SetRadiationLength(char *name, Float_t x2X0) {
 634   //
 635   // Set layer material [cm]
 636   //
 637
 638   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 639   if (!tmp)
 640     printf("Layer %s not found - cannot set layer material\n",name);
 641   else {
 642     tmp->fx2X0 = x2X0;
 643   }
 644 }
 645
 646 Float_t KMCDetector::GetRadiationLength(char *name) {
 647   //
 648   // Return layer radius [cm]
 649   //
 650
 651   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 652   if (!tmp)
 653     printf("Layer %s not found - cannot get layer material\n",name);
 654   else
 655     return tmp->fx2X0;
 656
 657   return 0;
 658
 659 }
 660
 661 void KMCDetector::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
 662   //
 663   // Set layer resolution in [cm]
 664   //
 665
 666   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 667   if (!tmp)
 668     printf("Layer %s not found - cannot set resolution\n",name);
 669   else {
 670     tmp->fPhiRes = phiRes;
 671     tmp->fZRes = zRes;
 672     tmp->fIsDead = (zRes==RIDICULOUS && phiRes==RIDICULOUS);
 673     ClassifyLayers();
 674   }
 675 }
 676
 677 Float_t KMCDetector::GetResolution(char *name, Int_t axis) {
 678   //
 679   // Return layer resolution in [cm]
 680   // axis = 0: resolution in rphi
 681   // axis = 1: resolution in z
 682   //
 683
 684   KMCLayer *tmp = GetLayer(name);
 685   if (!tmp)
 686     printf("Layer %s not found - cannot get resolution\n",name);
 687   else {
 688     if (axis==0) return tmp->fPhiRes;
 689     if (axis==1) return tmp->fZRes;
 690     printf("error: axis must be either 0 or 1 (rphi or z axis)\n");
 691   }
 692   return 0;
 693 }
 694
 695 void KMCDetector::SetLayerEfficiency(char *name, Float_t eff) {
 696   //
 697   // Set layer efficnecy (prop that his is missed within this layer)
 698   //
 699
 700   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 701   if (!tmp)
 702     printf("Layer %s not found - cannot set layer efficiency\n",name);
 703   else {
 704     tmp->fEff = eff;
 705   }
 706 }
 707
 708 Float_t KMCDetector::GetLayerEfficiency(char *name) {
 709   //
 710   // Get layer efficnecy (prop that his is missed within this layer)
 711   //
 712
 713   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 714   if (!tmp)
 715     printf("Layer %s not found - cannot get layer efficneicy\n",name);
 716   else
 717     return tmp->fEff;
 718
 719   return 0;
 720
 721 }
 722
 723 void KMCDetector::RemoveLayer(char *name) {
 724   //
 725   // Removes a layer from the list
 726   //
 727
 728   KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
 729   if (!tmp)
 730     printf("Layer %s not found - cannot remove it\n",name);
 731   else {
 732     fLayers.Remove(tmp);
 733     ClassifyLayers();
 734   }
 735 }
 736
 737
 738 void KMCDetector::PrintLayout() {
 739   //
 740   // Prints the detector layout
 741   //
 742
 743   printf("Detector %s: \"%s\"\n",GetName(),GetTitle());
 744
 745   if (fLayers.GetEntries()>0)
 746     printf("  Name \t\t r [cm] \t  X0 \t  phi & z res [um]\n");
 747
 748   KMCLayer *tmp = 0;
 749   for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
 750     tmp = (KMCLayer*)fLayers.At(i);
 751
 752     // don't print all the tpc layers
 753     TString name(tmp->GetName());
 754     if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
 755
 756     printf("%d. %s \t %03.2f   \t%1.4f\t  ",i,
 757            tmp->GetName(), tmp->fR, tmp->fx2X0);
 758     if (tmp->fPhiRes==RIDICULOUS)
 759       printf("  -  ");
 760     else
 761       printf("%3.0f   ",tmp->fPhiRes*10000);
 762     if (tmp->fZRes==RIDICULOUS)
 763       printf("  -\n");
 764     else
 765       printf("%3.0f\n",tmp->fZRes*10000);
 766   }
 767 }
 768
 769 void KMCDetector::PlotLayout(Int_t plotDead) {
 770   //
 771   // Plots the detector layout in Front view
 772   //
 773
 774   Double_t x0=0, y0=0;
 775
 776   TGraphErrors *gr = new TGraphErrors();
 777   gr->SetPoint(0,0,0);
 778   KMCLayer *lastLayer = (KMCLayer*)fLayers.At(fLayers.GetEntries()-1);  Double_t maxRad = lastLayer->fR;
 779   gr->SetPointError(0,maxRad,maxRad);
 780   gr->Draw("APE");
 781
 782
 783   KMCLayer *tmp = 0;
 784   for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {
 785     tmp = (KMCLayer*)fLayers.At(i);
 786
 787
 788     Double_t txtpos = tmp->fR;
 789     if ((tmp->IsDead())) txtpos*=-1; //
 790     TText *txt = new TText(x0,txtpos,tmp->GetName());
 791     txt->SetTextSizePixels(5); txt->SetTextAlign(21);
 792     if (!tmp->IsDead() || plotDead) txt->Draw();
 793
 794     TEllipse *layEl = new TEllipse(x0,y0,tmp->fR);
 795     //  layEl->SetFillColor(5);
 796     layEl->SetFillStyle(5001);
 797     layEl->SetLineStyle(tmp->IsDead()+1); // dashed if not active
 798     layEl->SetLineColor(4);
 799     TString name(tmp->GetName());
 800     if (!tmp->IsDead()) layEl->SetLineWidth(2);
 801     if (name.Contains("tpc") )  layEl->SetLineColor(29);
 802
 803     if (!tmp->IsDead() || plotDead) layEl->Draw();
 804
 805   }
 806
 807 }
 808
 809
 810
 811 void KMCDetector::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {
 812   //
 813   // Emulates the TPC
 814   //
 815   // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow
 816
 817
 818   AddLayer((char*)"IFCtpc",   77.8,0.01367, 0); // Inner Field cage (RS: set correct xrho for eloss)
 819
 820   // % Radiation Lengths ... Average per TPC row  (i.e. total/159 )
 821   Float_t x2X0PerRow = 0.000036;
 822
 823   Float_t tpcInnerRadialPitch  =    0.75 ;    // cm
 824   Float_t tpcMiddleRadialPitch =    1.0  ;    // cm
 825   Float_t tpcOuterRadialPitch  =    1.5  ;    // cm
 826   //  Float_t tpcInnerPadWidth     =    0.4  ;    // cm
 827   //  Float_t tpcMiddlePadWidth    =    0.6   ;   // cm
 828   //  Float_t tpcOuterPadWidth     =    0.6   ;   // cm
 829   Float_t innerRows            =   63 ;
 830   Float_t middleRows           =   64  ;
 831   Float_t outerRows            =   32  ;
 832   Float_t tpcRows            =   (innerRows + middleRows + outerRows) ;
 833   Float_t rowOneRadius         =   85.2  ;    // cm
 834   Float_t row64Radius          =  135.1  ;    // cm
 835   Float_t row128Radius         =  199.2  ;    // cm
 836
 837   for ( Int_t k = 0 ; k < tpcRows ; k++ ) {
 838
 839     Float_t rowRadius =0;
 840     if (k<innerRows)
 841       rowRadius =  rowOneRadius + k*tpcInnerRadialPitch ;
 842     else if ( k>=innerRows && k<(innerRows+middleRows) )
 843       rowRadius =  row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;
 844     else if (k>=(innerRows+middleRows) && k<tpcRows )
 845       rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;
 846
 847     if ( k%skip == 0 )
 848       AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0, phiResMean,zResMean);
 849     else
 850       AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0); // non "active" row
 851
 852
 853   }
 854
 855 }
 856
 857 void KMCDetector::RemoveTPC() {
 858
 859   // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)
 860   KMCLayer *tmp = 0;
 861   for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
 862     tmp = (KMCLayer*)fLayers.At(i);
 863     TString name(tmp->GetName());
 864     if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }
 865   }
 866   RemoveLayer((char*)"IFC");
 867
 868 }
 869
 870
 871 Double_t KMCDetector::ThetaMCS ( Double_t mass, Double_t x2X0, Double_t momentum ) const
 872 {
 873   //
 874   // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)
 875   //
 876
 877   Double_t beta  =  momentum / TMath::Sqrt(momentum*momentum+mass*mass)  ;
 878   Double_t theta =  0.0 ;    // Momentum and mass in GeV
 879   // if ( RadLength > 0 ) theta  =  0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );
 880   if ( x2X0 > 0 ) theta  =  0.0136 * TMath::Sqrt(x2X0) / ( beta * momentum ) * (1+0.038*TMath::Log(x2X0)) ;
 881   return (theta) ;
 882 }
 883
 884
 885 Double_t KMCDetector::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
 886 {
 887   // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm
 888   // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm
 889   // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius
 890   Double_t sx, sy, goodHit ;
 891   sx = 2 * TMath::Pi() *  searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;
 892   sy = 2 * TMath::Pi() *  searchRadiusZ    * searchRadiusZ    * HitDensity(radius) ;
 893   goodHit =  TMath::Sqrt(1./((1+sx)*(1+sy)))  ;
 894   return ( goodHit ) ;
 895 }
 896
 897
 898 Double_t KMCDetector::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
 899 {
 900   // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm
 901   // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
 902   Double_t sx, goodHit ;
 903   sx = 2 * TMath::Pi() *  searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;
 904   goodHit =  1./(1+sx) ;
 905   return ( goodHit ) ;
 906 }
 907
 908 Double_t KMCDetector::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
 909 {
 910   // Based on work by Ruben Shahoyen
 911   // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
 912   // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account
 913   // Following is correct for 2 DOF
 914
 915   Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
 916   Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
 917   Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c));
 918   return ( goodHit ) ;
 919 }
 920
 921 Double_t KMCDetector::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
 922 {
 923   // Based on work by Ruben Shahoyen
 924   // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)
 925
 926   Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
 927   Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
 928   Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));
 929   return ( nullHit ) ;
 930 }
 931
 932 Double_t KMCDetector::HitDensity ( Double_t radius )
 933 {
 934   // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.
 935   // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results
 936   // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.
 937   // Note that this function assumes we are working in CM and CM**2 [not meters].
 938   // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.
 939
 940   //  Double_t MaxRadiusSlowDet = 0.1; //?   // Maximum radius for slow detectors.  Fast detectors
 941                                         // and only fast detectors reside outside this radius.
 942   Double_t arealDensity = 0 ;
 943   if (radius<0.01) return 0;
 944
 945   if ( radius >= fMaxRadiusSlowDet )
 946     {
 947       arealDensity  = OneEventHitDensity(fdNdEtaCent,radius)  ; // Fast detectors see central collision density (only)
 948       arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius)  ;  // Increase density due to background
 949     }
 950
 951   if (radius < fMaxRadiusSlowDet )
 952     { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.
 953       arealDensity  = OneEventHitDensity(fdNdEtaCent,radius)
 954                     + IntegratedHitDensity(dNdEtaMinB,radius)
 955                     + UpcHitDensity(radius) ;
 956       arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;
 957       // Increase density due to background
 958     }
 959
 960   return ( arealDensity ) ;
 961 }
 962
 963
 964 double KMCDetector::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const
 965 {
 966   // This is for one event at the vertex.  No smearing.
 967   double den   = multiplicity / (2.*TMath::Pi()*radius*radius) * fDensFactorEta ; // 2 eta ?
 968   // note: surface of sphere is  '4*pi*r^2'
 969   //       surface of cylinder is '2*pi*r* h'
 970   return den ;
 971 }
 972
 973
 974 double KMCDetector::IntegratedHitDensity(Double_t multiplicity, Double_t radius)
 975 {
 976   // The integral of minBias events smeared over a gaussian vertex distribution.
 977   // Based on work by Yan Lu 12/20/2006, all radii in centimeters.
 978
 979   Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;
 980   Double_t den   = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;
 981
 982   // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.
 983   double dens1 = OneEventHitDensity(multiplicity,radius);
 984   if ( den < dens1 )  den = dens1;
 985
 986   return den ;
 987 }
 988
 989
 990 double KMCDetector::UpcHitDensity(Double_t radius)
 991 {
 992   // QED electrons ...
 993
 994   Double_t mUPCelectrons ;                                 ;
 995   //  mUPCelectrons =  fLhcUPCscale * (1.23 - radius/6.5)      ;  // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC
 996   mUPCelectrons = fLhcUPCscale*5456/(radius*radius)/dNdEtaMinB;      // Fit to 'Rossegger,Sadovsky'-Alice simulation
 997   if ( mUPCelectrons < 0 ) mUPCelectrons =  0.0             ;  // UPC electrons fall off quickly and don't go to large R
 998   mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ;  // UPCs increase Mulitiplicty ~ proportional to MinBias rate
 999   mUPCelectrons *= UPCBackgroundMultiplier                 ;  // Allow for an external multiplier (eg 0-1) to turn off UPC
1000
1001   return mUPCelectrons ;
1002 }
1003
1004
1005 double KMCDetector::Dist(double z, double r)
1006 {
1007   // Convolute dEta/dZ  distribution with assumed Gaussian of vertex z distribution
1008   // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm
1009   // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.
1010   Int_t    index  =  1     ;     // Start weight at 1 for Simpsons rule integration
1011   Int_t    nsteps =  301   ;     // NSteps must be odd for Simpson's rule to work
1012   double   dist   =  0.0   ;
1013   double   dz0    =  ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1)  ;  //cm
1014   double    z0    =  0.0   ;     //cm
1015   for(int i=0; i<nsteps; i++){
1016     if ( i == nsteps-1 ) index = 1 ;
1017     z0 = -4*SigmaD + i*dz0 ;
1018     dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) *
1019       (1/sqrt((z-z0)*(z-z0) + r*r)) ;
1020     if ( index != 4 ) index = 4; else index = 2 ;
1021   }
1022   return dist;
1023 }
1024
1025 #define  PZero   0.861  // Momentum of back to back decay particles in the CM frame
1026 #define  EPiZero 0.872  // Energy of the pion from a D0 decay at rest
1027 #define  EKZero  0.993  // Energy of the Kaon from a D0 decay at rest
1028
1029 Double_t KMCDetector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][20] ) const {
1030   // Math from Ron Longacre.  Note hardwired energy to bin conversion for PtK and PtPi.
1031
1032   Double_t const1  =  pt / D0Mass ;
1033   Double_t const2  =  TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;
1034   Double_t sum, ptPi, ptK ;
1035   Double_t effp, effk ;
1036
1037   sum = 0.0 ;
1038   for ( Int_t k = 0 ; k < 360 ; k++ )   {
1039
1040     Double_t theta = k * TMath::Pi() / 180. ;
1041
1042     ptPi = TMath::Sqrt(
1043                        PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
1044                        const1*const1*EPiZero*EPiZero -
1045                        2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +
1046                        PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
1047                        ) ;
1048
1049     ptK = TMath::Sqrt(
1050                       PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
1051                       const1*const1*EKZero*EKZero +
1052                       2*PZero*TMath::Cos(theta)*const2*const1*EKZero +
1053                       PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
1054                       ) ;
1055
1056     // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays
1057     Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ;
1058     Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ;
1059
1060     if ( pionindex >= 20 ) pionindex = 399 ;
1061     if ( pionindex >= 0 )   effp = corrEfficiency[0][pionindex] ;
1062     if ( pionindex <  0 )   effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd
1063     if ( effp < 0 )         effp = 0 ;
1064
1065     if ( kaonindex >= 20 ) kaonindex = 399 ;
1066     if ( kaonindex >= 0 )   effk = corrEfficiency[1][kaonindex] ;
1067     if ( kaonindex <  0 )   effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd
1068     if ( effk < 0 )         effk = 0 ;
1069
1070     // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.
1071
1072     sum += effp * effk ;
1073
1074   }
1075
1076   Double_t mean =sum/360;
1077   return mean ;
1078
1079 }
1080
1081 KMCProbe* KMCDetector::PrepareKalmanTrack(double pt, double lambda, double mass, int charge, double phi, double x,double y, double z)
1082 {
1083   // Prepare trackable Kalman track at the farthest position
1084   //
1085   // Set track parameters
1086   // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis
1087   fProbe.Reset();
1088   fProbe.SetMass(mass);
1089   KMCProbe* probe = new KMCProbe(fProbe);
1090   double *trPars = (double*)probe->GetParameter();
1091   double *trCov  = (double*)probe->GetCovariance();
1092   double xyz[3] = {x,y,z};
1093   probe->Global2LocalPosition(xyz,phi);
1094   probe->Set(xyz[0],phi,trPars,trCov);
1095   trPars[KMCProbe::kY] = xyz[1];
1096   trPars[KMCProbe::kZ] = xyz[2];
1097   trPars[KMCProbe::kSnp] = 0;                       //            track along X axis at the vertex
1098   trPars[KMCProbe::kTgl] = TMath::Tan(lambda);                // dip
1099   trPars[KMCProbe::kPtI] = charge/pt;               //            q/pt
1100   //
1101   // put tiny errors to propagate to the outer-most radius
1102   trCov[KMCProbe::kY2] = trCov[KMCProbe::kZ2] = trCov[KMCProbe::kSnp2] = trCov[KMCProbe::kTgl2] = trCov[KMCProbe::kPtI2] = 1e-20;
1103   fProbe = *probe;  // store original track
1104   //
1105   // propagate to last layer
1106   fLastActiveLayerTracked = 0;
1107   for (Int_t j=0; j<=fLastActiveLayer; j++) {
1108     KMCLayer* lr = GetLayer(j);
1109     lr->Reset();
1110     //
1111     if (!PropagateToLayer(probe,lr,1)) break;
1112     if (!probe->CorrectForMeanMaterial(lr, kFALSE)) break;
1113     //
1114     lr->fClCorr.Set(probe->GetY(),probe->GetZ(), probe->GetX(), probe->GetAlpha());
1115     if (!lr->IsDead()) fLastActiveLayerTracked = j;
1116   }
1117   probe->ResetCovMat();// reset cov.matrix
1118   printf("Last active layer trracked: %d (out of %d)\n",fLastActiveLayerTracked,fLastActiveLayer);
1119   //
1120   return probe;
1121 }
1122
1123
1124
1125 TGraph * KMCDetector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
1126   //
1127   // returns the momentum resolution
1128   //
1129
1130   TGraph *graph = new TGraph(20, fTransMomenta, fMomentumRes);
1131   graph->SetTitle("Momentum Resolution .vs. Pt" ) ;
1132   //  graph->GetXaxis()->SetRangeUser(0.,5.0) ;
1133   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1134   graph->GetXaxis()->CenterTitle();
1135   graph->GetXaxis()->SetNoExponent(1) ;
1136   graph->GetXaxis()->SetMoreLogLabels(1) ;
1137   graph->GetYaxis()->SetTitle("Momentum Resolution (%)") ;
1138   graph->GetYaxis()->CenterTitle();
1139
1140   graph->SetMaximum(20) ;
1141   graph->SetMinimum(0.1) ;
1142   graph->SetLineColor(color);
1143   graph->SetMarkerColor(color);
1144   graph->SetLineWidth(linewidth);
1145
1146   return graph;
1147
1148 }
1149
1150 TGraph * KMCDetector::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {
1151
1152   // Returns the pointing resolution
1153   // axis = 0 ... rphi pointing resolution
1154   // axis = 1 ... z pointing resolution
1155   //
1156
1157   TGraph * graph =  0;
1158
1159   if (axis==0) {
1160     graph = new TGraph ( 20, fTransMomenta, fResolutionRPhi ) ;
1161     graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;
1162     graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
1163   } else {
1164     graph =  new TGraph ( 20, fTransMomenta, fResolutionZ ) ;
1165     graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
1166     graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;
1167   }
1168
1169   graph->SetMinimum(1) ;
1170   graph->SetMaximum(300.1) ;
1171   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1172   graph->GetXaxis()->CenterTitle();
1173   graph->GetXaxis()->SetNoExponent(1) ;
1174   graph->GetXaxis()->SetMoreLogLabels(1) ;
1175   graph->GetYaxis()->CenterTitle();
1176
1177   graph->SetLineWidth(linewidth);
1178   graph->SetLineColor(color);
1179   graph->SetMarkerColor(color);
1180
1181   return graph;
1182
1183 }
1184
1185
1186 TGraph * KMCDetector::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {
1187   //
1188   // returns the Pointing resolution (accoring to Telescope equation)
1189   // axis =0 ... in rphi
1190   // axis =1 ... in z
1191   //
1192
1193   Double_t resolution[20];
1194
1195   Double_t layerResolution[2];
1196   Double_t layerRadius[2];
1197   Double_t layerThickness[2];
1198
1199   Int_t count =0; // search two first active layers
1200   printf("Telescope equation for layers:  ");
1201   for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
1202     KMCLayer *l = (KMCLayer*)fLayers.At(i);
1203     if (!l->IsDead() && l->fR>0) {
1204       layerRadius[count]     = l->fR;
1205       layerThickness[count]  = l->fx2X0;
1206       if (axis==0) {
1207         layerResolution[count] = l->fPhiRes;
1208       } else {
1209         layerResolution[count] = l->fZRes;
1210       }
1211       printf("%s, ",l->GetName());
1212       count++;
1213     }
1214     if (count>=2) break;
1215   }
1216   printf("\n");
1217
1218   Double_t pt, momentum, thickness,aMCS ;
1219   Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1220
1221   for ( Int_t i = 0 ; i < 20 ; i++ ) {
1222     // Reference data as if first two layers were acting all alone
1223     pt  =  fTransMomenta[i]  ;
1224     momentum = pt / TMath::Cos(lambda)   ;  // Total momentum
1225     resolution[i] =  layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1]
1226       +  layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;
1227     resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;
1228     thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;
1229     aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;
1230     resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;
1231     resolution[i] =  TMath::Sqrt(resolution[i]) * 10000.0 ;  // result in microns
1232   }
1233
1234
1235
1236   TGraph* graph = new TGraph ( 20, fTransMomenta, resolution ) ;
1237
1238   if (axis==0) {
1239     graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;
1240     graph->GetYaxis()->SetTitle("RPhi Pointing Resolution (#mum) ") ;
1241   } else {
1242     graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
1243     graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum) ") ;
1244   }
1245   graph->SetMinimum(1) ;
1246   graph->SetMaximum(300.1) ;
1247   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1248   graph->GetXaxis()->CenterTitle();
1249   graph->GetXaxis()->SetNoExponent(1) ;
1250   graph->GetXaxis()->SetMoreLogLabels(1) ;
1251   graph->GetYaxis()->CenterTitle();
1252
1253   graph->SetLineColor(color);
1254   graph->SetMarkerColor(color);
1255   graph->SetLineStyle(kDashed);
1256   graph->SetLineWidth(linewidth);
1257
1258   return graph;
1259
1260 }
1261
1262 TGraph * KMCDetector::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {
1263   //
1264   // particle = 0 ... choosen particle (setted particleMass)
1265   // particle = 1 ... Pion
1266   // particle = 2 ... Kaon
1267   // particle = 3 ... D0
1268   //
1269   Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1270
1271   Double_t particleEfficiency[20]; // with chosen particle mass
1272   Double_t kaonEfficiency[20], pionEfficiency[20], d0efficiency[20];
1273   Double_t partEfficiency[2][20];
1274
1275   if (particle != 0) {
1276     // resulting Pion and Kaon efficiency scaled with overall efficiency
1277     Double_t doNotDecayFactor;
1278     for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
1279
1280       for ( Int_t j = 0 ; j < 20 ; j++ ) {
1281         // JT Test Let the kaon decay.  If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
1282         Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda)           ;  // Total momentum at average rapidity
1283         if ( massloop == 1 ) { // KAON
1284           doNotDecayFactor  = TMath::Exp(-130/(371*momentum/KaonMass)) ;  // Decay length for kaon is 371 cm.
1285           kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
1286         } else { // PION
1287           doNotDecayFactor = 1.0 ;
1288           pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
1289         }
1290         partEfficiency[0][j] = pionEfficiency[j];
1291         partEfficiency[1][j] = kaonEfficiency[j];
1292       }
1293     }
1294
1295     // resulting estimate of the D0 efficiency
1296     for ( Int_t j = 0 ; j < 20 ; j++ ) {
1297       d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);
1298     }
1299   } else {
1300     for ( Int_t j = 0 ; j < 20 ; j++ ) {
1301       particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;
1302       // NOTE: Decay factor (see kaon) should be included to be realiable
1303     }
1304   }
1305
1306   for ( Int_t j = 0 ; j < 20 ; j++ ) {
1307     pionEfficiency[j]     *= 100;
1308     kaonEfficiency[j]     *= 100;
1309     d0efficiency[j]       *= 100;
1310     particleEfficiency[j] *= 100;
1311   }
1312
1313   TGraph * graph =  0;
1314   if (particle==0) {
1315     graph = new TGraph ( 20, fTransMomenta, particleEfficiency ) ; // choosen mass
1316     graph->SetLineWidth(1);
1317   }  else if (particle==1) {
1318     graph = new TGraph ( 20, fTransMomenta, pionEfficiency ) ;
1319     graph->SetLineWidth(1);
1320   }  else if (particle ==2) {
1321     graph = new TGraph ( 20, fTransMomenta, kaonEfficiency ) ;
1322     graph->SetLineWidth(1);
1323   }  else if (particle ==3) {
1324     graph = new TGraph ( 20, fTransMomenta, d0efficiency ) ;
1325     graph->SetLineStyle(kDashed);
1326   } else
1327     return 0;
1328
1329   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1330   graph->GetXaxis()->CenterTitle();
1331   graph->GetXaxis()->SetNoExponent(1) ;
1332   graph->GetXaxis()->SetMoreLogLabels(1) ;
1333   graph->GetYaxis()->SetTitle("Efficiency (%)") ;
1334   graph->GetYaxis()->CenterTitle();
1335
1336   graph->SetMinimum(0.01) ;
1337   graph->SetMaximum(100)  ;
1338
1339   graph->SetLineColor(color);
1340   graph->SetMarkerColor(color);
1341   graph->SetLineWidth(linewidth);
1342
1343   return graph;
1344 }
1345
1346 TGraph * KMCDetector::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) {
1347   //
1348   // particle = 0 ... choosen particle (setted particleMass)
1349   // particle = 1 ... Pion
1350   // particle = 2 ... Kaon
1351   //
1352
1353   Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
1354
1355   Double_t particleFake[20]; // with chosen particle mass
1356   Double_t kaonFake[20], pionFake[20];
1357   Double_t partFake[2][20];
1358
1359   if (particle != 0) {
1360     // resulting Pion and Kaon efficiency scaled with overall efficiency
1361     Double_t doNotDecayFactor;
1362     for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
1363
1364       for ( Int_t j = 0 ; j < 20 ; j++ ) {
1365         // JT Test Let the kaon decay.  If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
1366         Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda)           ;  // Total momentum at average rapidity
1367         if ( massloop == 1 ) { // KAON
1368           doNotDecayFactor  = TMath::Exp(-130/(371*momentum/KaonMass)) ;  // Decay length for kaon is 371 cm.
1369           kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;
1370         } else { // PION
1371           pionFake[j] = fFake[0][j] ;
1372         }
1373         partFake[0][j] = pionFake[j];
1374         partFake[1][j] = kaonFake[j];
1375       }
1376     }
1377
1378   } else {
1379     for ( Int_t j = 0 ; j < 20 ; j++ ) {
1380       particleFake[j] = fFake[2][j];
1381       // NOTE: Decay factor (see kaon) should be included to be realiable
1382     }
1383   }
1384
1385   for ( Int_t j = 0 ; j < 20 ; j++ ) {
1386     pionFake[j]     *= 100;
1387     kaonFake[j]     *= 100;
1388     particleFake[j] *= 100;
1389   }
1390
1391   TGraph * graph =  0;
1392   if (particle==0) {
1393     graph = new TGraph ( 20, fTransMomenta, particleFake ) ; // choosen mass
1394     graph->SetLineWidth(1);
1395   }  else if (particle==1) {
1396     graph = new TGraph ( 20, fTransMomenta, pionFake ) ;
1397     graph->SetLineWidth(1);
1398   }  else if (particle ==2) {
1399     graph = new TGraph ( 20, fTransMomenta, kaonFake ) ;
1400     graph->SetLineWidth(1);
1401   }
1402
1403   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1404   graph->GetXaxis()->CenterTitle();
1405   graph->GetXaxis()->SetNoExponent(1) ;
1406   graph->GetXaxis()->SetMoreLogLabels(1) ;
1407   graph->GetYaxis()->SetTitle("Fake (%)") ;
1408   graph->GetYaxis()->CenterTitle();
1409
1410   graph->SetMinimum(0.01) ;
1411   graph->SetMaximum(100)  ;
1412
1413   graph->SetLineColor(color);
1414   graph->SetMarkerColor(color);
1415   graph->SetLineWidth(linewidth);
1416
1417   return graph;
1418 }
1419
1420
1421 TGraph* KMCDetector::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {
1422   //
1423   // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)
1424   // mode 0: impact parameter (convolution of pointing and vertex resolution)
1425   // mode 1: pointing resolution
1426   // mode 2: vtx resolution
1427
1428
1429   TGraph *graph = new TGraph();
1430
1431   //  TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?
1432   TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); //
1433   TFormula vtxResZ("vtxResZ","600/(x+6)+10"); //
1434
1435   TGraph *trackRes = GetGraphPointingResolution(axis,1);
1436   Double_t *pt = trackRes->GetX();
1437   Double_t *trRes = trackRes->GetY();
1438   for (Int_t ip =0; ip<trackRes->GetN(); ip++) {
1439     Double_t vtxRes = 0;
1440     if (axis==0)
1441       vtxRes = vtxResRPhi.Eval(pt[ip]);
1442     else
1443       vtxRes = vtxResZ.Eval(pt[ip]);
1444
1445     if (mode==0)
1446       graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));
1447     else if (mode ==1)
1448       graph->SetPoint(ip,pt[ip],trRes[ip]);
1449     else
1450       graph->SetPoint(ip,pt[ip],vtxRes);
1451   }
1452
1453   graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;
1454   graph->GetYaxis()->SetTitle("d_{0} r#phi resolution (#mum)") ;
1455
1456   graph->SetMinimum(1) ;
1457   graph->SetMaximum(300.1) ;
1458   graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
1459   graph->GetXaxis()->CenterTitle();
1460   graph->GetXaxis()->SetNoExponent(1) ;
1461   graph->GetXaxis()->SetMoreLogLabels(1) ;
1462   graph->GetYaxis()->CenterTitle();
1463
1464   graph->SetLineColor(color);
1465   graph->SetMarkerColor(color);
1466   graph->SetLineWidth(linewidth);
1467
1468   return graph;
1469
1470 }
1471
1472 TGraph* KMCDetector::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
1473   //
1474   // returns graph according to the number
1475   //
1476   switch(number) {
1477   case 1:
1478     return GetGraphPointingResolution(0,color, linewidth); // dr
1479   case 2:
1480     return GetGraphPointingResolution(1,color, linewidth); // dz
1481   case 3:
1482     return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele
1483   case 4:
1484     return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele
1485   case 5:
1486     return GetGraphMomentumResolution(color, linewidth); // pt resolution
1487   case 10:
1488     return GetGraphRecoEfficiency(0, color, linewidth);  // tracked particle
1489   case 11:
1490     return GetGraphRecoEfficiency(1, color, linewidth);  // eff. pion
1491   case 12:
1492     return GetGraphRecoEfficiency(2, color, linewidth);  // eff. kaon
1493   case 13:
1494     return GetGraphRecoEfficiency(3, color, linewidth);  // eff. D0
1495   case 15:
1496     return GetGraphRecoFakes(0, color, linewidth);  // Fake tracked particle
1497   case 16:
1498     return GetGraphRecoFakes(1, color, linewidth);  // Fake pion
1499   case 17:
1500     return GetGraphRecoFakes(2, color, linewidth);  // Fake kaon
1501   default:
1502     printf(" Error: chosen graph number not valid\n");
1503   }
1504   return 0;
1505
1506 }
1507
1508 void KMCDetector::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon, int setVer) {
1509
1510   // All New configuration with X0 = 0.3 and resolution = 4 microns
1511
1512   AddLayer((char*)"bpipe_its",2.0,0.0022, 0.092); // beam pipe, 0.5 mm Be
1513   AddLayer((char*)"vertex_its",     0,     0); // dummy vertex for matrix calculation
1514   if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {
1515     printf("vertex will work as constraint: %.4f %.4f\n",fgVtxConstraint[0],fgVtxConstraint[1]);
1516     SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);
1517   }
1518   //
1519   // new ideal Pixel properties?
1520   Double_t x0     = 0.0050;
1521   Double_t resRPhi = 0.0006;
1522   Double_t resZ   = 0.0006;
1523   Double_t xrho = 0.0116;  // assume 0.5mm of silicon for eloss
1524   //
1525   if (flagMon) {
1526     x0 *= 3./5.;
1527     xrho *=3./5.;
1528     resRPhi = 0.0004;
1529     resZ   = 0.0004;
1530   }
1531
1532   double sclD = 1;
1533   double sclZ = 1;
1534
1535   // default all new
1536   if (setVer<=0) {
1537     AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1538     AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ);
1539     AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ);
1540     AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ);
1541     AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ);
1542     AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ);
1543     AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ);
1544   }
1545   else if (setVer==1) {
1546     AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1547     AddLayer((char*)"its2",  2.8 ,  x0, xrho, resRPhi, resZ);
1548     AddLayer((char*)"its3",  3.6 ,  x0, xrho, resRPhi, resZ);
1549     AddLayer((char*)"its4", 20.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1550     AddLayer((char*)"its5", 22.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1551     AddLayer((char*)"its6", 43.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1552     AddLayer((char*)"its7", 43.6 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1553     //
1554     /*
1555     UInt_t patt[] = {
1556       KMCTrackSummary::Bits(1,1,1),
1557       KMCTrackSummary::Bits(0,0,0,1,1),
1558       KMCTrackSummary::Bits(0,0,0,0,0,1,1)
1559     };
1560     RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
1561     */
1562   }
1563   else if (setVer==2) {
1564     AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1565     AddLayer((char*)"its1a", 2.8 ,  x0, xrho, resRPhi, resZ);
1566     AddLayer((char*)"its2",  3.6 ,  x0, xrho, resRPhi, resZ);
1567     AddLayer((char*)"its2a", 4.2 ,  x0, xrho, resRPhi, resZ);
1568     AddLayer((char*)"its3", 20.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1569     AddLayer((char*)"its4", 22.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1570     AddLayer((char*)"its5", 33.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1571     AddLayer((char*)"its6", 43.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1572     AddLayer((char*)"its7", 43.6 ,  x0, xrho, resRPhi*sclD, resZ*sclZ);
1573     //
1574     /*
1575     UInt_t patt[] = {
1576       KMCTrackSummary::Bits(1,1,1,1),
1577       KMCTrackSummary::Bits(0,0,0,0,1,1),
1578       KMCTrackSummary::Bits(0,0,0,0,0,0,1,1,1)
1579     };
1580     RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
1581     */
1582   }
1583    /*
1584   // last 2 layers strips
1585   double resRPStr=0.0020, resZStr = 0.0830;
1586   AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1587   AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ);
1588   AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ);
1589   AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ);
1590   AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ);
1591   AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPStr, resZStr);
1592   AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPStr, resZStr);
1593    */
1594   //*/
1595   /*
1596   // resolution scaled with R as res(2.2) * R/2.2
1597   AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1598   AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ);
1599   AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ);
1600   AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ);
1601   AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ);
1602   AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ);
1603   AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ);
1604   KMCLayer* lr0 = 0;
1605   for (int i=0;i<GetNActiveITSLayers();i++) {
1606     KMCLayer *lr = GetActiveLayer(i);
1607     if (lr->IsVertex()) continue;
1608     if (lr0==0) {
1609       lr0=lr;
1610       //  continue;
1611     }
1612     double scl = 5*lr->GetRadius()/lr0->GetRadius();
1613     SetResolution((char*)lr->GetName(), resRPhi*scl, resZ*scl*4);
1614   }
1615   */
1616   /*
1617   // 1st 2 layers double
1618   AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1619   AddLayer((char*)"its1a",  2.8 ,  x0, xrho, resRPhi, resZ);
1620
1621   AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ);
1622   AddLayer((char*)"its2a",  4.2 ,  x0, xrho, resRPhi, resZ);
1623
1624   AddLayer((char*)"its3a",  6.4 ,  x0, xrho, resRPhi, resZ);
1625   AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ);
1626
1627   AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ);
1628   AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ);
1629   AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ);
1630   AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ);
1631   */
1632   /*
1633   // last 4 layers doubled
1634   AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1635   AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ);
1636   AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ);
1637
1638   AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ);
1639   AddLayer((char*)"its4a", 12.8 ,  x0, xrho, resRPhi, resZ);
1640
1641   AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ);
1642   AddLayer((char*)"its5a", 23.9 ,  x0, xrho, resRPhi, resZ);
1643
1644   AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ);
1645   AddLayer((char*)"its6a", 40.0 ,  x0, xrho, resRPhi, resZ);
1646
1647   AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ);
1648   AddLayer((char*)"its7a", 43.4 ,  x0, xrho, resRPhi, resZ);
1649   */
1650
1651   /* //last 3 lr together
1652   AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ);
1653   AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ);
1654   AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ);
1655   AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ);
1656   AddLayer((char*)"its5", 42.2 ,  x0, xrho, resRPhi, resZ);
1657   AddLayer((char*)"its6", 42.6 ,  x0, xrho, resRPhi, resZ);
1658   AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ);
1659   */
1660   if (flagTPC) {
1661     AddTPC(0.1,0.1);                        // TPC
1662   }
1663   PrintLayout();
1664 }
1665
1666 void KMCDetector::MakeAliceCurrent(Bool_t flagTPC, Int_t AlignResiduals) {
1667
1668   // Numbers taken from
1669   // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks
1670   // number for misalingment: private communication with Andrea Dainese
1671
1672   //  1.48e-01, 2.48e-01,2.57e-01, 1.34e-01, 3.34e-01,3.50e-01, 2.22e-01, 2.38e-01,2.25e-01};
1673
1674   AddLayer((char*)"vertex_its",     0,     0); // dummy vertex for matrix calculation
1675   if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {
1676     printf("vertex will work as constraint: %.4f %.4f",fgVtxConstraint[0],fgVtxConstraint[1]);
1677     SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);
1678   }
1679   //
1680   AddLayer((char*)"bpipe_its",2.94,0.0022, 1.48e-01); // beam pipe
1681   AddLayer((char*)"tshld1_its",11.5,0.0065, 1.34e-01); // Thermal shield  // 1.3% /2
1682   AddLayer((char*)"tshld2_its",31.0,0.0108, 2.22e-01); // Thermal shield  // 1.3% /2
1683   //
1684   if (flagTPC) {
1685     AddTPC(0.1,0.1);                        // TPC
1686   }
1687   // Adding the ITS - current configuration
1688
1689   if (AlignResiduals==0) {
1690     AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, 0.0012, 0.0130);
1691     AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, 0.0012, 0.0130);
1692     AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, 0.0035, 0.0025);
1693     AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, 0.0035, 0.0025);
1694     AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, 0.0020, 0.0830);
1695     AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, 0.0020, 0.0830);
1696     /*
1697     UInt_t patt[] = {
1698       KMCTrackSummary::Bits(1,1),
1699       KMCTrackSummary::Bits(0,0,1,1),
1700       KMCTrackSummary::Bits(0,0,0,0,1,1)
1701     };
1702     RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
1703     */
1704   } else if (AlignResiduals==1) {
1705
1706     // tracking errors ...
1707     // (Additional systematic errors due to misalignments) ...
1708     // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020);  // [cm]
1709     // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
1710
1711     AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010),
1712              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1713     AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),
1714              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1715     AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),
1716              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1717     AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),
1718              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1719     AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
1720              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1721     AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
1722              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1723
1724   } else if (AlignResiduals==2) {
1725
1726
1727     // tracking errors ... PLUS ... module misalignment
1728
1729     // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020);  // [cm]
1730     // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
1731
1732     //  the ITS modules are misalignment with small gaussian smearings with
1733     //  sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
1734
1735     AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008),
1736              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1737     AddLayer((char*)"spd2_ots", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008),
1738              TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
1739     AddLayer((char*)"sdd1_ots",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
1740              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1741     AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
1742              TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
1743     AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
1744              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1745     AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
1746              TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
1747
1748   } else {
1749
1750       //  the ITS modules are misalignment with small gaussian smearings with
1751       //  sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
1752       //  unknown in Z ????
1753
1754     AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
1755              TMath::Sqrt(0.0130*0.0130+0.000*0.000));
1756     AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
1757              TMath::Sqrt(0.0130*0.0130+0.000*0.000));
1758     AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
1759              TMath::Sqrt(0.0025*0.0025+0.000*0.000));
1760     AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
1761              TMath::Sqrt(0.0025*0.0025+0.000*0.000));
1762     AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
1763              TMath::Sqrt(0.0830*0.0830+0.000*0.000));
1764     AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
1765              TMath::Sqrt(0.0830*0.0830+0.000*0.000));
1766   }
1767
1768 }
1769
1770
1771 void KMCDetector::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {
1772   //
1773   // Produces the standard performace plots
1774   //
1775
1776   if (!add) {
1777
1778     TCanvas *c1 = new TCanvas("c1","c1");//,100,100,500,500);
1779     c1->Divide(2,2);
1780
1781     c1->cd(1);  gPad->SetGridx();   gPad->SetGridy();
1782     gPad->SetLogx();
1783     TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth);
1784     eff->SetTitle("Efficiencies");
1785     eff->Draw("AL");
1786     if (!onlyPionEff) {
1787       GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
1788       GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
1789     }
1790     c1->cd(2); gPad->SetGridx();   gPad->SetGridy();
1791     gPad->SetLogy();  gPad->SetLogx();
1792     GetGraphMomentumResolution(color,linewidth)->Draw("AL");
1793
1794     c1->cd(3); gPad->SetGridx();   gPad->SetGridy();
1795     gPad->SetLogx();
1796     GetGraphPointingResolution(0,color,linewidth)->Draw("AL");
1797
1798     c1->cd(4); gPad->SetGridx();   gPad->SetGridy();
1799     gPad->SetLogx();
1800     GetGraphPointingResolution(1,color,linewidth)->Draw("AL");
1801
1802   } else {
1803
1804     TVirtualPad *c1 = gPad->GetMother();
1805
1806     c1->cd(1);
1807     GetGraphRecoEfficiency(1,color,linewidth)->Draw("L");
1808     if (!onlyPionEff) {
1809       GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
1810       GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
1811     }
1812     c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("L");
1813
1814     c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("L");
1815
1816     c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("L");
1817
1818   }
1819
1820 }
1821
1822 void KMCDetector::ApplyMS(KMCProbe* trc, double x2X0) const
1823 {
1824   // simulate random modification of track params due to the MS
1825   if (x2X0<=0) return;
1826   double alpha = trc->GetAlpha(); // store original alpha
1827   double mass = trc->GetMass();
1828   //
1829   double snp = trc->GetSnp();
1830   double dip = trc->GetTgl();
1831   Double_t angle=TMath::Sqrt((1.+ dip*dip)/((1-snp)*(1.+snp)));
1832   x2X0 *= angle;
1833   //
1834   static double covCorr[15],covDum[21]={0};
1835   static double mom[3],pos[3];
1836   double *cov = (double*) trc->GetCovariance();
1837   memcpy(covCorr,cov,15*sizeof(double));
1838   trc->GetXYZ(pos);
1839   trc->GetPxPyPz(mom);
1840   double pt2 = mom[0]*mom[0]+mom[1]*mom[1];
1841   double pt = TMath::Sqrt(pt2);
1842   double ptot2 = pt2 + mom[2]*mom[2];
1843   double ptot  = TMath::Sqrt(ptot2);
1844   double beta = ptot/TMath::Sqrt(ptot2 + mass*mass);
1845   double sigth = TMath::Sqrt(x2X0)*0.014/(ptot*beta);
1846   //
1847   // a la geant
1848   double phiSC = gRandom->Rndm()*TMath::Pi();
1849   double thtSC = gRandom->Gaus(0,1.4142*sigth);
1850   //  printf("MS phi: %+.5f tht: %+.5f\n",phiSC,thtSC);
1851   double sn = TMath::Sin(thtSC);
1852   double dx = sn*TMath::Sin(phiSC);
1853   double dy = sn*TMath::Cos(phiSC);
1854   double dz = TMath::Cos(thtSC);
1855   double v[3];
1856   //  printf("Before: %+.3e %+.3e %+.3e | MS: %+.3e %+.3e\n",mom[0],mom[1],mom[2],thtSC,phiSC);
1857   for (int i=3;i--;) mom[i] /= ptot;
1858   double vmm = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
1859   if (!IsZero(pt)) {
1860     double pd1 = mom[0]/vmm;
1861     double pd2 = mom[1]/vmm;
1862     v[0] = pd1*mom[2]*dx - pd2*dy + mom[0]*dz;
1863     v[1] = pd2*mom[2]*dx + pd1*dy + mom[1]*dz;
1864     v[2] = -vmm*dx                + mom[2]*dz;
1865   }
1866   else {
1867     v[0] = dx;
1868     v[1] = dy;
1869     v[2] = dz*TMath::Sign(1.,mom[2]);
1870   }
1871   double nrm = TMath::Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1872   //  printf("before :%+e %+e %+e  || %+e %+e %+e %+e\n",mom[0],mom[1],mom[2],  sigth, x2X0, pt, beta);
1873   //  trc->Print();
1874   // direction cosines -> p
1875   for (int i=3;i--;) mom[i] = ptot*v[i]/nrm;
1876   //  printf("After : %+.3e %+.3e %+.3e\n",mom[0],mom[1],mom[2]);
1877   trc->Set(pos,mom,covDum,trc->Charge());
1878   //
1879   trc->Rotate(alpha);
1880   memcpy(cov,covCorr,15*sizeof(double));
1881   //
1882 }
1883
1884 //________________________________________________________________________________
1885 Bool_t KMCDetector::TransportKalmanTrackWithMS(KMCProbe *probTr, int maxLr) const
1886 {
1887   // Transport track till layer maxLr, applying random MS
1888   //
1889   for (Int_t j=0; j<maxLr; j++) {
1890     KMCLayer* lr0 = (KMCLayer*)fLayers.At(j);
1891     KMCLayer* lr  = (KMCLayer*)fLayers.At(j+1);
1892     if (!lr0->IsVertex()) {
1893       ApplyMS(probTr,lr0->fx2X0); // apply MS
1894       if (!probTr->CorrectForMeanMaterial(lr0,kFALSE)) return kFALSE;
1895     }
1896     //
1897     if (!PropagateToLayer(probTr,lr,1)) return kFALSE;
1898     // store randomized cluster local coordinates and phi
1899     lr->ResetBgClusters();
1900     double rz,ry;
1901     gRandom->Rannor(rz,ry);
1902     lr->GetMCCluster()->Set(probTr->GetY()+ry*lr->GetPhiRes(),probTr->GetZ()+rz*lr->GetZRes(),
1903                                probTr->GetX(), probTr->GetAlpha() );
1904     //
1905   }
1906   //
1907   return kTRUE;
1908 }
1909
1910 //________________________________________________________________________________
1911 Bool_t KMCDetector::SolveSingleTrack(Double_t mass, Double_t pt, Double_t eta, TObjArray* sumArr,int nMC, int offset)
1912 {
1913   // analityc and fullMC (nMC trials) evaluaion of tracks with given kinematics.
1914   // the results are filled in KMCTrackSummary objects provided via summArr array
1915   //
1916   int progressP = 10;//0; // print progress in percents
1917   //
1918   progressP = int(nMC*0.01*progressP);
1919
1920   if (!SolveSingleTrackViaKalman(mass,pt,eta)) return kFALSE;
1921   //
1922   // Store non-updated track errors of inward propagated seed >>>>>>>>
1923   int maxLr = fLastActiveITSLayer + offset;
1924   if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;
1925   KMCProbe probeTmp = fProbe; // original probe at vertex
1926   KMCLayer* lr = 0;
1927   for (Int_t j=1; j<=maxLr; j++) {
1928     lr = GetLayer(j);
1929     //    printf("Here0: %d\n",j);
1930     if (!PropagateToLayer(&probeTmp,lr,1)) return 0;
1931     if (j!=maxLr) if (!probeTmp.CorrectForMeanMaterial(lr, kFALSE)) return 0;
1932     //    printf("Prelim. Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(probeTmp.GetSigmaY2()),TMath::Sqrt(probeTmp.GetSigmaZ2()));
1933   }
1934   for (Int_t j=maxLr; j>0; j--) {
1935     lr = GetLayer(j);
1936     //    printf("Here1: %d\n",j);
1937     if (j!=maxLr) if (!PropagateToLayer(&probeTmp,lr,-1)) return 0;
1938     lr->fSig2EstD = probeTmp.GetSigmaY2();
1939     lr->fSig2EstZ = probeTmp.GetSigmaZ2();
1940     //    probeTmp.Print("l");
1941     printf("Natural Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(lr->fSig2EstD),TMath::Sqrt(lr->fSig2EstZ));
1942     if (!probeTmp.CorrectForMeanMaterial(lr, kTRUE)) return 0;
1943   }
1944   // Store non-updated track errors of inward propagated seed <<<<<<<<
1945   //
1946   int nsm = sumArr ? sumArr->GetEntriesFast() : 0;
1947   KMCLayer* vtx = GetLayer(0);
1948   //
1949   for (int i=0;i<nsm;i++) {
1950     KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(i);
1951     if (!tsm) continue;
1952     tsm->SetRefProbe( GetProbeTrack() ); // attach reference track (generated)
1953     tsm->SetAnProbe( vtx->GetAnProbe() ); // attach analitycal solution
1954   }
1955   //
1956   TStopwatch sw;
1957   sw.Start();
1958   for (int it=0;it<nMC;it++) {
1959     printf("ev: %d\n",it);
1960     SolveSingleTrackViaKalmanMC(offset);
1961     KMCProbe* trc = vtx->GetWinnerMCTrack();
1962     vtx->GetMCTracks()->Print();
1963     if (progressP==1 || (progressP>0 &&  (it%progressP)==0)) {
1964       printf("%d%% done |",it*100/nMC);
1965       sw.Stop(); sw.Print(); sw.Start(kFALSE);
1966     }
1967     for (int ism=nsm;ism--;) { // account the track in each of summaries
1968       KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(ism);
1969       if (!tsm) continue;
1970       tsm->AddUpdCalls(GetUpdCalls());
1971       tsm->AddTrack(trc);
1972     }
1973   }
1974   //
1975   sw.Stop();
1976   printf("Total time: "); sw.Print();
1977   return kTRUE;
1978 }
1979
1980 //________________________________________________________________________________
1981 Int_t KMCDetector::GetLayerID(int actID) const
1982 {
1983   // find physical layer id from active id
1984   if (actID<0 || actID>fNActiveLayers) return -1;
1985   int start = actID<fNActiveITSLayers ? fLastActiveITSLayer : fLastActiveLayer;
1986   for (int i=start+1;i--;) {
1987     if (GetLayer(i)->GetActiveID()==actID) return i;
1988   }
1989   return -1;
1990 }
1991
1992 //________________________________________________________________________________
1993 KMCProbe* KMCDetector::KalmanSmooth(int actLr, int actMin,int actMax) const
1994 {
1995   // estimate kalman smoothed track params at given active lr
1996   // from fit at layers actMin:actMax (excluding actLr)
1997   // SolveSingleTrackViaKalman must have been called before
1998   //
1999   if (actMin>actMax) swap(actMin,actMax);
2000   if (actMax>=fNActiveLayers) actMax = fNActiveLayers-1;
2001   int nlrfit = actMax-actMin;
2002   if (actLr>=actMin && actLr<=actMax) nlrfit-=1;
2003   if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}
2004   static KMCProbe iwd,owd;
2005   //
2006   // find phisical layer id's
2007   int pLr  = GetLayerID(actLr);
2008   int pMin = GetLayerID(actMin);
2009   int pMax = GetLayerID(actMax);
2010   //
2011   //  printf(">>> %d %d %d\n",pLr, pMin,pMax);
2012   Bool_t useIwd=kFALSE, useOwd=kFALSE;
2013   if (pLr<pMax) { // need inward piece
2014     iwd = GetLayer(pMax)->fTrCorr;
2015     iwd.ResetCovMat();
2016     iwd.GetHitsPatt() = 0;
2017     for (int i=pMax;i>=pLr;i--) {
2018       KMCLayer* lr = GetLayer(i);
2019       //      printf("IWD %d\n",i);
2020       if (!lr->IsDead() && i!=pLr && i>=pMin) if (!UpdateTrack(&iwd,lr,&lr->fClCorr))  return 0;
2021       if (i!=pLr) {
2022         if (!iwd.CorrectForMeanMaterial(lr,kTRUE)) return 0; // correct for materials of this layer
2023         if (!PropagateToLayer(&iwd,GetLayer(i-1),-1)) return 0;      // propagate to next layer
2024       }
2025       //  printf("IWD%d:  ",i); iwd.Print("l");
2026     }
2027     useIwd = kTRUE;
2028   }
2029   if (pLr>pMin) { // need outward piece
2030     owd = GetLayer(pMin)->fTrCorr;
2031     owd.ResetCovMat();
2032     owd.GetHitsPatt() = 0;
2033     for (int i=pMin;i<=pLr;i++) {
2034       KMCLayer* lr = GetLayer(i);
2035       //      printf("OWD %d\n",i);
2036       if (!lr->IsDead() && i!=pLr && i<=pMax) if (!UpdateTrack(&owd,lr,&lr->fClCorr))  return 0;
2037       if (i!=pLr) {
2038         if (!owd.CorrectForMeanMaterial(lr,0)) return 0; // correct for materials of this layer
2039         if (!PropagateToLayer(&owd,GetLayer(i+1), 1)) return 0;      // propagate to next layer
2040       }
2041       //      printf("OWD%d:  ",i); owd.Print("l");
2042     }
2043     useOwd = kTRUE;
2044   }
2045   //
2046   // was this extrapolation outside the fit range?
2047   if (!useIwd) return (KMCProbe*)&owd;
2048   if (!useOwd) return (KMCProbe*)&iwd;
2049   //
2050   // weight both tracks
2051   if (!iwd.Propagate(owd.GetAlpha(),owd.GetX(),fBFieldG)) return 0;
2052   double meas[2] = {owd.GetY(),owd.GetZ()};
2053   double measErr2[3] = {owd.GetSigmaY2(), owd.GetSigmaZY(), owd.GetSigmaZ2()};
2054   //  printf("Weighting\n");
2055   //  owd.Print("l");
2056   //  iwd.Print("l");
2057   if (!iwd.Update(meas,measErr2)) return 0;
2058   iwd.GetHitsPatt() |= owd.GetHitsPatt();
2059
2060   //  printf("->\n");
2061   //  iwd.Print("l");
2062
2063   return (KMCProbe*)&iwd;
2064   //
2065 }
2066
2067 //________________________________________________________________________________
2068 KMCProbe* KMCDetector::KalmanSmoothFull(int actLr, int actMin,int actMax) const
2069 {
2070   // estimate kalman smoothed track params at given active lr
2071   // from fit at layers actMin:actMax (excluding actLr)
2072   // SolveSingleTrackViaKalman must have been called before
2073   //
2074   static TClonesArray prediction("KMCProbe",10);
2075   static TClonesArray update("KMCProbe",10);
2076   static KMCProbe res;
2077   //
2078   if (actMin>actMax) swap(actMin,actMax);
2079   int nlrfit = actMax-actMin;
2080   if (actLr>=actMin && actLr<=actMax) nlrfit-=1;
2081   if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}
2082   //
2083   // find phisical layer id's
2084   int pLr  = GetLayerID(actLr);
2085   int pMin = GetLayerID(actMin);
2086   int pMax = GetLayerID(actMax);
2087   //
2088   int dir=0,dirInt=0;
2089   if      (pLr<=pMin) dir=-1; // inward extrapolation
2090   else if (pLr>=pMax) dir= 1; // outward extrapolation
2091   else if (actMax-actLr >= actLr-actMin) dirInt = -1; // inward  interpolation (the test point is closer to inner layer)
2092   else    dirInt = 1;                                 // outward interpolation (the test point is closer to outer layer)
2093   //
2094   if (dir!=0) { // no sens to do smoothing: simple Kalman filtering extrapolation
2095     int start = dir<0 ? pMax : pMin;
2096     res = GetLayer(start)->fTrCorr;
2097     res.ResetCovMat();
2098     KMCLayer* lr = 0;
2099     for (int i=(dir<0?pMax:pMin); i!=pLr; i+=dir) { // track till nearest layer to pLr
2100       lr = GetLayer(i);
2101       if (!lr->IsDead() && !(i<pMin ||i>pMax)) if (!UpdateTrack(&res,lr,&lr->fClCorr))  return 0; // update only with layers in fit range
2102       if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE))   return 0; // correct for materials of this layer
2103       if (!PropagateToLayer(&res,GetLayer(i+dir),dir))            return 0; // propagate to next layer
2104     }
2105     if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE))   return 0; // correct for materials of this nearest layer
2106     if (!PropagateToLayer(&res,GetLayer(pLr), dir)) return 0; // propagate to test layer
2107     return (KMCProbe*)&res;
2108   }
2109   //
2110   // too bad, need to do real filtering
2111   //
2112   int start = dirInt<0 ? pMax : pMin;
2113   int stop  = dirInt<0 ? pMin-1 : pMax+1;
2114   res = GetLayer(start)->fTrCorr;
2115   res.ResetCovMat();
2116   KMCLayer* lr = 0;
2117   int count = 0;
2118   for (int i=start; i!=stop; i+=dirInt) { // track in full range, storing updates and predictions
2119     new(prediction[count]) KMCProbe(res);
2120     lr = GetLayer(i);
2121     if (!lr->IsDead() && i!=pLr) if (!UpdateTrack(&res,lr,&lr->fClCorr))  return 0; // update only with layers in fit range
2122     new(update[count]) KMCProbe(res);
2123     if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE))   return 0; // correct for materials of this layer
2124     if (!PropagateToLayer(&res,GetLayer(i+dir),dir))            return 0; // propagate to next layer
2125     count++;
2126   }
2127   return (KMCProbe*)&res;
2128   //
2129 }
2130
2131 //________________________________________________________________________________
2132 Bool_t KMCDetector::SolveSingleTrackViaKalman(Double_t mass, Double_t pt, Double_t eta)
2133 {
2134   // analytical estimate of tracking resolutions
2135   //  fProbe.SetUseLogTermMS(kTRUE);
2136   //
2137   if (fMinITSHits>fNActiveITSLayers) {fMinITSHits = fNActiveITSLayers; printf("Redefined request of min N ITS hits to %d\n",fMinITSHits);}
2138   if (TMath::Abs(eta)<1e-3) fDensFactorEta = 1.;
2139   else {
2140     fDensFactorEta = TMath::Tan( 2.*TMath::ATan(TMath::Exp(-TMath::Abs(eta))) );
2141     fDensFactorEta = 1./TMath::Sqrt( 1. + 1./fDensFactorEta/fDensFactorEta);
2142   }
2143   double lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-eta));
2144   KMCProbe* probe = PrepareKalmanTrack(pt,lambda,mass,-1);
2145   if (!probe) return kFALSE;
2146   //
2147   KMCLayer *lr = 0;
2148   //
2149   //
2150   // Start the track fitting --------------------------------------------------------
2151   //
2152   // Back-propagate the covariance matrix along the track.
2153   // Kalman loop over the layers
2154   //
2155   KMCProbe* currTr = 0;
2156   lr = (KMCLayer*)fLayers.At(fLastActiveLayerTracked);
2157   lr->fTrCorr = *probe;
2158   delete probe; // rethink...
2159   //
2160   for (Int_t j=fLastActiveLayerTracked; j--; ) {  // Layer loop
2161     //
2162     KMCLayer *lrP = lr;
2163     lr = (KMCLayer*)fLayers.At(j);
2164     //
2165     lr->fTrCorr = lrP->fTrCorr;
2166     currTr = &lr->fTrCorr;
2167     currTr->ResetHit(lrP->GetActiveID());
2168     //
2169     // if there was a measurement on prev layer, update the track
2170     if (!lrP->IsDead()) { // include measurement
2171       KMCCluster cl(currTr->GetY(),currTr->GetZ(), currTr->GetX(), currTr->GetAlpha());
2172       if (!UpdateTrack(currTr,lrP,&cl))  return kFALSE;
2173     }
2174     if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) return kFALSE; // correct for materials of this layer
2175     if (!PropagateToLayer(currTr,lr,-1)) return kFALSE;      // propagate to current layer
2176     //
2177   } // end loop over layers
2178   //
2179   return kTRUE;
2180 }
2181
2182 //____________________________________________________________
2183 Bool_t KMCDetector::SolveSingleTrackViaKalmanMC(int offset)
2184 {
2185   // MC estimate of tracking resolutions/effiencies. Requires that the SolveSingleTrackViaKalman
2186   // was called before, since it uses data filled by this method
2187   //
2188   // The MC tracking will be done starting from fLastActiveITSLayer + offset (before analytical estimate will be used)
2189   //
2190   // At this point, the fProbe contains the track params generated at vertex.
2191   // Clone it and propagate to target layer to generate hit positions affected by MS
2192   //
2193   fUpdCalls = 0.;
2194   KMCProbe *currTrP=0,*currTr=0;
2195   int maxLr = fLastActiveITSLayer + offset;
2196   if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;
2197   ResetMCTracks(maxLr);
2198   KMCLayer* lr = (KMCLayer*)fLayers.At(maxLr);
2199   currTr = lr->AddMCTrack(&fProbe); // start with original track at vertex
2200   //
2201   if (!TransportKalmanTrackWithMS(currTr, maxLr)) return kFALSE; // transport it to outermost layer where full MC is done
2202   //
2203   if (fLastActiveITSLayer<fLastActiveLayerTracked) { // prolongation from TPC
2204     // start from correct track propagated from above till maxLr
2205     double *covMS = (double*)currTr->GetCovariance();
2206     const double *covIdeal =lr->fTrCorr.GetCovariance();
2207     for (int i=15;i--;) covMS[i] = covIdeal[i];
2208   }
2209   else { // ITS SA: randomize the starting point
2210     //    double *pars = (double*)currTr->GetParameter();
2211     //    pars[0] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaY2()));
2212     //    pars[1] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaZ2()));
2213     //
2214     currTr->ResetCovMat();
2215     /*
2216     double *trCov  = (double*)currTr->GetCovariance();
2217     double *trPars = (double*)currTr->GetParameter();
2218     const double kLargeErr2PtI = 0.3*0.3;
2219     trCov[14] = TMath::Max(trCov[14],kLargeErr2PtI*trPars[4]*trPars[4]);
2220     */
2221   }
2222   //
2223   for (Int_t j=maxLr; j--; ) {  // Layer loop
2224     //
2225     KMCLayer *lrP = lr;
2226     lr = (KMCLayer*)fLayers.At(j);
2227     int ntPrev = lrP->GetNMCTracks();
2228     //
2229     if (lrP->IsDead()) { // for passive layer just propagate the copy of all tracks of prev layer >>>
2230       for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
2231         currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
2232         currTr = lr->AddMCTrack( currTrP );
2233         if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer
2234         if (!PropagateToLayer(currTr,lr,-1))      {currTr->Kill(); continue;} // propagate to current layer
2235       }
2236       continue;
2237     } // treatment of dead layer <<<
2238     //
2239     if (lrP->IsTPC()) { // we don't consider bg hits in TPC, just update with MC cluster
2240       for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
2241         currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
2242         currTr = lr->AddMCTrack( currTrP );
2243         if (!UpdateTrack(currTr, lrP, lrP->GetMCCluster(), kTRUE)) {currTr->Kill(); continue;} // update with correct MC cl.
2244         if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer
2245         if (!PropagateToLayer(currTr,lr,-1))      {currTr->Kill(); continue;} // propagate to current layer
2246       }
2247       continue;
2248     } // treatment of ideal (TPC?) layer <<<
2249     //
2250     // active layer under eff. study (ITS?): propagate copy of every track to MC cluster frame (to have them all in the same frame)
2251     // and calculate the limits of bg generation
2252     KMCCluster* clMC = lrP->GetMCCluster();
2253     if (lrP->GetLayerEff()<gRandom->Rndm()) clMC->Kill(); // simulate inefficiency
2254     ResetSearchLimits();
2255     int nseeds = 0;
2256     for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
2257       currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
2258       currTr = lr->AddMCTrack( currTrP );
2259       if (!currTr->PropagateToCluster(clMC,fBFieldG)) {currTr->Kill(); continue;} // go to MC cluster
2260       if ( !(currTr->GetNITSHits()>0 && currTr->GetNITSHits()==currTr->GetNFakeITSHits()) ) UpdateSearchLimits(currTr, lrP); // RS
2261       nseeds++;
2262     }
2263     //
2264     //    printf("%3d seeds\n",nseeds);
2265     if (fUseBackground && lrP->IsITS()) GenBgClusters(lrP); //  generate background hits
2266     //
2267     ntPrev = lr->GetNMCTracks();
2268     for (int itr=ntPrev;itr--;) { // loop over all tracks PROPAGATED from previous layer to clusters frame on previous layer
2269       currTrP = lr->GetMCTrack(itr); // this is a seed from prev layer. The new clusters are attached to its copies, the seed itself
2270                                      // will be propagated w/o cluster update if it does not violate requested "reconstructed" track settings
2271       if (currTrP->IsKilled()) continue;
2272       //printf("Check    %d %p %d\n",itr,currTrP,currTrP->GetUniqueID()); currTrP->Print();
2273       CheckTrackProlongations(currTrP, lr, lrP);
2274       if (NeedToKill(currTrP)) currTrP->Kill(); // kill track which was not updated at lrP
2275       //currTrP->Kill(); // kill track which was not updated at lrP
2276     }
2277     //
2278     lr->GetMCTracks()->Sort();
2279     int ntTot = lr->GetNMCTracks(); // propagate max amount of allowed tracks to current layer
2280     if (ntTot>fMaxSeedToPropagate && fMaxSeedToPropagate>0) {
2281       for (int itr=ntTot;itr>=fMaxSeedToPropagate;itr--)  lr->GetMCTracks()->RemoveAt(itr);
2282       ntTot = fMaxSeedToPropagate;
2283     }
2284     //
2285     for (int itr=ntTot;itr--;) {
2286       currTr = lr->GetMCTrack(itr);
2287       if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill();continue;} // correct for materials of prev. layer
2288       if (!PropagateToLayer(currTr,lr,-1))      {currTr->Kill();continue;} // propagate to current layer
2289     }
2290     AliDebug(1,Form("Got %d tracks on layer %s",ntTot,lr->GetName()));
2291     //    lr->GetMCTracks()->Print();
2292     //
2293   } // end loop over layers
2294   //
2295   // do we use vertex constraint?
2296   KMCLayer *vtx = GetLayer(0);
2297   if (!vtx->IsDead() && vtx->IsITS()) {
2298     int ntr = vtx->GetNMCTracks();
2299     for (int itr=0;itr<ntr;itr++) {
2300       currTr = vtx->GetMCTrack(itr);
2301       if (currTr->IsKilled()) continue;
2302       KMCCluster* clv = vtx->GetMCCluster();
2303       double meas[2] = {clv->GetY(),clv->GetZ()};
2304       double measErr2[3] = {vtx->fPhiRes*vtx->fPhiRes,0,vtx->fZRes*vtx->fZRes};
2305       double chi2v = currTr->GetPredictedChi2(meas,measErr2);
2306       currTr->AddHit(vtx->GetActiveID(), chi2v, -1);
2307       currTr->SetInnerLrChecked(vtx->GetActiveID());
2308       if (NeedToKill(currTr)) currTr->Kill();
2309       // if (vtx->IsITS()) {if (!UpdateTrack(currTr, vtx, vtx->GetMCCluster(), kFALSE)) {currTr->Kill();continue;}}
2310     }
2311   }
2312   EliminateUnrelated();
2313
2314   return kTRUE;
2315 }
2316
2317 //____________________________________________________________________________
2318 Bool_t KMCDetector::PropagateToLayer(KMCProbe* trc, KMCLayer* lr, int dir) const
2319 {
2320   // bring the track to layer and rotat to frame normal to its surface
2321   if (!trc->PropagateToR(lr->fR,fBFieldG, dir)) return kFALSE;
2322   //
2323   // rotate to frame with X axis normal to the surface (defined by ideal track)
2324   if (!lr->IsVertex()) {
2325     double pos[3];
2326     trc->GetXYZ(pos);  // lab position
2327     double phi = TMath::ATan2(pos[1],pos[0]);
2328     if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;
2329     if (!trc->Rotate(phi)) {
2330       AliDebug(2,Form("Failed to rotate to the frame (phi:%+.3f) of layer at %.2f at XYZ: %+.3f %+.3f %+.3f",
2331                       phi,lr->fR,pos[0],pos[1],pos[2]));
2332       if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");
2333       return kFALSE;
2334     }
2335   }
2336   //
2337   return kTRUE;
2338 }
2339
2340 //____________________________________________________________________________
2341 Bool_t KMCDetector::UpdateTrack(KMCProbe* trc, KMCLayer* lr, KMCCluster* cl, Bool_t goToCluster) const
2342 {
2343   // update track with measured cluster
2344   // propagate to cluster
2345   double meas[2] = {cl->GetY(),cl->GetZ()}; // ideal cluster coordinate
2346   double measErr2[3] = {lr->fPhiRes*lr->fPhiRes,0,lr->fZRes*lr->fZRes};
2347   //
2348   if (goToCluster) if (!trc->PropagateToCluster(cl,fBFieldG)) return kFALSE; // track was not propagated to cluster frame
2349   //
2350   double chi2 = trc->GetPredictedChi2(meas,measErr2);
2351   //  if (chi2>fMaxChi2Cl) return kTRUE; // chi2 is too large
2352   //
2353   if (!trc->Update(meas,measErr2)) {
2354     AliDebug(2,Form("layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",
2355                     lr->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));
2356     if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");
2357     return kFALSE;
2358   }
2359   trc->AddHit(lr->GetActiveID(), chi2);
2360   //
2361   return kTRUE;
2362 }
2363
2364 //____________________________________________________________________________
2365 Int_t KMCDetector::GenBgClusters(KMCLayer* lr)
2366 {
2367   // Generate fake clusters in precalculated RPhi,Z range
2368   if (fNBgLimits<1) return 0; // limits were not set - no track was prolongated
2369   //
2370   // Fix search limits to avoid seeds which will anyway point very far from the vertex
2371   double tolY = TMath::Sqrt(lr->fSig2EstD)*fMaxChi2ClSQ;
2372   double tolZ = TMath::Sqrt(lr->fSig2EstZ)*fMaxChi2ClSQ;
2373
2374   //  printf("Before: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);
2375   if (fBgYMin < lr->fClCorr.fY-tolY) fBgYMin = lr->fClCorr.fY-tolY;
2376   if (fBgYMax > lr->fClCorr.fY+tolY) fBgYMax = lr->fClCorr.fY+tolY;
2377   if (fBgZMin < lr->fClCorr.fZ-tolZ) fBgZMin = lr->fClCorr.fZ-tolZ;
2378   if (fBgZMax > lr->fClCorr.fZ+tolZ) fBgZMax = lr->fClCorr.fZ+tolZ;
2379   //printf("After: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);
2380   //
2381   double dy = fBgYMax - fBgYMin;
2382   double dz = fBgZMax - fBgZMin;
2383   double surf = dy*dz;               // surface of generation
2384   if (surf<0) return 0;
2385   double poissProb = surf*HitDensity(lr->fR)*lr->GetLayerEff();
2386   AliDebug(2,Form("Bg for Lr %s (r=%.2f) : Density %.2f on surface %.2e [%+.4f : %+.4f][%+.4f %+.4f]",
2387                   lr->GetName(),lr->fR,HitDensity(lr->fR),surf,fBgYMin,fBgYMax,fBgZMin,fBgZMax));
2388   int nFakesGen = gRandom->Poisson( poissProb ); // preliminary number of extra clusters to test
2389   KMCCluster *refCl = lr->GetMCCluster();
2390   double sig2y = lr->GetPhiRes()*lr->GetPhiRes();
2391   double sig2z = lr->GetZRes()*lr->GetZRes();
2392   for (int ic=nFakesGen;ic--;) {
2393     double y = fBgYMin+dy*gRandom->Rndm();
2394     double z = fBgZMin+dz*gRandom->Rndm();
2395     double dfy = y-refCl->GetY();
2396     double dfz = z-refCl->GetZ();
2397     double dist = (dfy*dfy)/sig2y + (dfz*dfz)/sig2z;
2398     if (dist<4) continue; // avoid overlap with MC cluster
2399     lr->AddBgCluster(y, z, refCl->GetX(), refCl->GetPhi());
2400   }
2401   AliDebug(2,Form("Added %6d noise clusters on lr %s (poisson Prob=%8.2f for surface %.2e) DY:%7.4f DZ: %7.4f",
2402                    lr->GetNBgClusters(),lr->GetName(),poissProb,surf,dy,dz));
2403   return nFakesGen;
2404   //
2405 }
2406
2407 //____________________________________________________________________________
2408 void KMCDetector::UpdateSearchLimits(KMCProbe* probe, KMCLayer* lr)
2409 {
2410   // define the search window for track on layer (where the bg hist will be generated)
2411   static double *currYMin = fBgYMinTr.GetArray();
2412   static double *currYMax = fBgYMaxTr.GetArray();
2413   static double *currZMin = fBgZMinTr.GetArray();
2414   static double *currZMax = fBgZMaxTr.GetArray();
2415   //
2416   double sizeY = probe->GetSigmaY2(), sizeZ = probe->GetSigmaZ2();
2417   //  /*
2418   if (probe->GetNITSHits()<3) sizeY = 10*lr->fSig2EstD;
2419   if (probe->GetNITSHits()<2) sizeZ = 10*lr->fSig2EstZ;
2420   sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY+lr->fPhiRes*lr->fPhiRes); // max deviation in rphi to accept
2421   sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ+lr->fZRes*lr->fZRes); // max deviation in dz to accept
2422   //  */
2423   //
2424   /*
2425   if (probe->GetNITSHits()<3) sizeY = 1./(1./lr->fSig2EstD + 1./sizeY);
2426   if (probe->GetNITSHits()<2) sizeZ = 1./(1./lr->fSig2EstZ + 1./sizeZ);
2427   sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY); // max deviation in rphi to accept
2428   sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ); // max deviation in dz to accept
2429   */
2430   //
2431   //  if (sizeY>2) sizeY=2;
2432   //  if (sizeZ>2) sizeZ=2;
2433   //  printf("Sizes at %s: %.5f %.5f\n",lr->GetName(), sizeY,sizeZ);
2434   //
2435   if (fNBgLimits>=fBgYMinTr.GetSize()) { // expand arrays, update pointers
2436     fBgYMinTr.Set(2*(fNBgLimits+1));
2437     fBgYMaxTr.Set(2*(fNBgLimits+1));
2438     fBgZMinTr.Set(2*(fNBgLimits+1));
2439     fBgZMaxTr.Set(2*(fNBgLimits+1));
2440     currYMin = fBgYMinTr.GetArray();
2441     currYMax = fBgYMaxTr.GetArray();
2442     currZMin = fBgZMinTr.GetArray();
2443     currZMax = fBgZMaxTr.GetArray();
2444   }
2445   if (fBgYMin > (currYMin[fNBgLimits]=probe->GetY()-sizeY) ) fBgYMin = currYMin[fNBgLimits];
2446   if (fBgYMax < (currYMax[fNBgLimits]=probe->GetY()+sizeY) ) fBgYMax = currYMax[fNBgLimits];
2447   if (fBgZMin > (currZMin[fNBgLimits]=probe->GetZ()-sizeZ) ) fBgZMin = currZMin[fNBgLimits];
2448   if (fBgZMax < (currZMax[fNBgLimits]=probe->GetZ()+sizeZ) ) fBgZMax = currZMax[fNBgLimits];
2449   if (AliLog::GetGlobalDebugLevel()>=2) {
2450     probe->Print("l");
2451     AliInfo(Form("Seed%3d Lr %s limits for y:%+8.4f z:%+8.4f [%+.4f : %+.4f][%+.4f %+.4f]",fNBgLimits,lr->GetName(),probe->GetY(),probe->GetZ(),currYMin[fNBgLimits],currYMax[fNBgLimits],currZMin[fNBgLimits],currZMax[fNBgLimits]));
2452     AliInfo(Form("Global Limits Lr %s                            [%+.4f : %+.4f][%+.4f %+.4f]",lr->GetName(),fBgYMin,fBgYMax,fBgZMin,fBgZMax));
2453     AliInfo(Form("MC Cluster: %+.4f : %+.4f",lr->fClMC.fY, lr->fClMC.fZ));
2454   }
2455   probe->SetUniqueID(fNBgLimits++);
2456   //
2457   if (lr->IsITS() && probe->GetNFakeITSHits()==0) {
2458     if (fHMCLrResidRPhi) fHMCLrResidRPhi->Fill(probe->GetY() - lr->GetMCCluster()->GetY(), lr->GetActiveID());
2459     if (fHMCLrResidZ)    fHMCLrResidZ->Fill(probe->GetZ() - lr->GetMCCluster()->GetZ(),lr->GetActiveID());
2460   }
2461   //
2462 }
2463
2464 //____________________________________________________________________________
2465 void KMCDetector::CheckTrackProlongations(KMCProbe *probe, KMCLayer* lr, KMCLayer* lrP)
2466 {
2467   // explore prolongation of probe from lrP to lr with all possible clusters of lrP
2468   // the probe is already brought to clusters frame
2469   int nCl = lrP->GetNBgClusters();
2470   double measErr2[3] = {lrP->fPhiRes*lrP->fPhiRes,0,lrP->fZRes*lrP->fZRes};
2471   double meas[2] = {0,0};
2472   UInt_t tmpID = probe->GetUniqueID();
2473   double yMin = fBgYMinTr[tmpID];
2474   double yMax = fBgYMaxTr[tmpID];
2475   double zMin = fBgZMinTr[tmpID];
2476   double zMax = fBgZMaxTr[tmpID];
2477   //
2478   probe->SetInnerLrChecked(lrP->GetActiveID());
2479   for (int icl=-1;icl<nCl;icl++) {
2480     KMCCluster* cl = icl<0 ? lrP->GetMCCluster() : lrP->GetBgCluster(icl);  // -1 is for true MC cluster
2481     if (cl->IsKilled()) {
2482       if (AliLog::GetGlobalDebugLevel()>1) {printf("Skip cluster %d ",icl); cl->Print();}
2483       continue;
2484     }
2485     double y = cl->GetY();
2486     double z = cl->GetZ();
2487     AliDebug(2,Form("Check seed%d against cl#%d out of %d at layer %s | y:%+8.4f z:%+8.4f [%+.4f:%+.4f]  [%+.4f:%+.4f]",tmpID,icl,nCl,lrP->GetName(),y,z,yMin,yMax,zMin,zMax));
2488     if (AliLog::GetGlobalDebugLevel()>0) {
2489       if (icl==-1 && probe->GetNFakeITSHits()==0) {
2490         meas[0] = y; meas[1] = z;
2491         double chi2a = probe->GetPredictedChi2(meas,measErr2);
2492         if (chi2a>fMaxChi2Cl || (y<yMin || y>yMax) || (z<zMin || z>zMax)) {
2493           probe->Print();
2494           printf("Loosing good point (y:%+8.4f z:%+8.4f) on lr %s: chi2: %.2f  | dy:%+8.4f dz:%+8.4f [%+.4f:%+.4f]  [%+.4f:%+.4f] |x: %.2f %.2f | phi: %.2f %.2f\n",
2495                  y,z,lrP->GetName(),chi2a,y-probe->GetY(),z-probe->GetZ(),yMin,yMax,zMin,zMax, probe->GetX(), cl->GetX(), probe->GetAlpha(), cl->GetPhi());
2496         }
2497       }
2498     }
2499     if (y<yMin || y>yMax) continue; // preliminary check on Y
2500     if (z<zMin || z>zMax) continue; // preliminary check on Z
2501     meas[0] = y; meas[1] = z;
2502     double chi2 = probe->GetPredictedChi2(meas,measErr2);
2503     if (fHMCLrChi2 && probe->GetNFakeITSHits()==0 && icl==-1) fHMCLrChi2->Fill(chi2,lrP->GetActiveID());
2504     AliDebug(2,Form("Seed-to-cluster chi2 = Chi2=%.2f",chi2));
2505     if (chi2>fMaxChi2Cl) continue;
2506     //
2507     // update track copy
2508     KMCProbe* newTr = lr->AddMCTrack( probe );
2509     fUpdCalls++;
2510     if (!newTr->Update(meas,measErr2)) {
2511       AliDebug(2,Form("Layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",
2512                       lrP->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));
2513       if (AliLog::GetGlobalDebugLevel()>1) newTr->Print("l");
2514       newTr->Kill();
2515       continue;
2516     }
2517     newTr->AddHit(lrP->GetActiveID(), chi2, icl);
2518     if (AliLog::GetGlobalDebugLevel()>1) {
2519       AliInfo("Cloned updated track is:");
2520       newTr->Print();
2521     }
2522     if (NeedToKill(newTr)) newTr->Kill();
2523   }
2524   //
2525 }
2526
2527 //____________________________________________________________________________
2528 void KMCDetector::ResetMCTracks(Int_t maxLr)
2529 {
2530   int nl = GetNLayers();
2531   if (maxLr<0 || maxLr>=nl) maxLr = nl-1;
2532   for (int i=maxLr+1;i--;) GetLayer(i)->ResetMCTracks();
2533 }
2534
2535 //____________________________________________________________________________
2536 Bool_t KMCDetector::NeedToKill(KMCProbe* probe) const
2537 {
2538   // check if the seed at given layer (last one where update was tried)
2539   // still has chances to be reconstructed
2540   const Bool_t kModeKillMiss = kFALSE;
2541   //
2542   Bool_t kill = kFALSE;
2543   while (1) {
2544     int il = probe->GetInnerLayerChecked();
2545     int nITS = probe->GetNITSHits();
2546     int nITSMax = nITS + il; // maximum it can have
2547     if (nITSMax<fMinITSHits) {
2548       kill = kTRUE;
2549       break;
2550     }    // has no chance to collect enough ITS hits
2551     //
2552     int ngr = fPattITS.GetSize();
2553     if (ngr>0) { // check pattern
2554       UInt_t patt = probe->GetHitsPatt();
2555       // complete the layers not checked yet
2556       for (int i=il;i--;) patt |= (0x1<<i);
2557       for (int ig=ngr;ig--;)
2558         if (!(((UInt_t)fPattITS[ig]) & patt)) {
2559           kill = kTRUE;
2560           break;
2561         }
2562       //
2563     }
2564     //
2565     if (nITS>2) {  // check if smallest possible norm chi2/ndf is acceptable
2566       double chi2min = probe->GetChi2();
2567       if (kModeKillMiss) {
2568         int nMiss = fNActiveITSLayers - probe->GetInnerLayerChecked() - nITS; // layers already missed
2569         chi2min = nMiss*probe->GetMissingHitPenalty();
2570       }
2571       chi2min /= ((nITSMax<<1)-KMCProbe::kNDOF);
2572       if (chi2min>fMaxNormChi2NDF) {
2573         kill = kTRUE;
2574         break;
2575       }
2576     }
2577     //
2578     // loose vertex constraint
2579     double dst;
2580     if (nITS>=2) {
2581       probe->GetZAt(0,fBFieldG,dst);
2582       //printf("Zd (F%d): %f\n",probe->GetNFakeITSHits(),dst);
2583       if (TMath::Abs(dst)>10.) {
2584         kill = kTRUE;
2585         break;
2586       }
2587     }
2588     if (nITS>=3) {
2589       probe->GetYAt(0,fBFieldG,dst);
2590       //printf("Dd (F%d): %f\n",probe->GetNFakeITSHits(),dst);
2591       if (TMath::Abs(dst)>10.) {
2592         kill = kTRUE;
2593         break;
2594       }
2595     }
2596     //
2597     break;
2598   }
2599   if (kill && AliLog::GetGlobalDebugLevel()>1 && probe->GetNFakeITSHits()==0) {
2600     printf("Killing good seed, last upd layer was %d\n",probe->GetInnerLayerChecked());
2601     probe->Print("l");
2602   }
2603   return kill;
2604 }
2605
2606 //_____________________________________________________________________
2607 void KMCDetector::EliminateUnrelated()
2608 {
2609   // kill useless tracks
2610   KMCLayer* lr = GetLayer(0);
2611   int ntr = lr->GetNMCTracks();
2612   int nval = 0;
2613   for (int itr=0;itr<ntr;itr++) {
2614     KMCProbe* probe = lr->GetMCTrack(itr);
2615     if (probe->IsKilled()) continue;
2616     if (probe->GetNITSHits()-probe->GetNFakeITSHits()<1) {probe->Kill(); continue;}
2617     nval++;
2618   }
2619   lr->GetMCTracks()->Sort();
2620   const int kDump = 0;
2621   if (kDump>0) {
2622     printf("Valid %d out of %d\n",nval, ntr);
2623     ntr = ntr>kDump ? kDump:0;
2624     for (int itr=0;itr<ntr;itr++) {
2625       lr->GetMCTrack(itr)->Print();
2626     }
2627   }
2628 }
2629
2630 //_____________________________________________________________________
2631 void KMCDetector::RequirePattern(UInt_t *patt, int groups)
2632 {
2633   if (groups<1) {fPattITS.Set(0); return;}
2634   fPattITS.Set(groups);
2635   for (int i=0;i<groups;i++) fPattITS[i] = patt[i];
2636 }
2637
2638 //_____________________________________________________________________
2639 void KMCDetector::CalcHardSearchLimits(double dzv)
2640 {
2641   //
2642   TArrayD zlims;
2643   zlims.Set(fNActiveITSLayers);
2644   for (int il0=0;il0<fNActiveITSLayers;il0++) {
2645     KMCLayer* lr0 = GetActiveLayer(il0);
2646     double angZTol = dzv/lr0->GetRadius();
2647     for (int il1=0;il1<fNActiveITSLayers;il1++) {
2648       if (il1==il0) continue;
2649       KMCLayer* lr1 = GetActiveLayer(il1);
2650       double ztol = angZTol*TMath::Abs(lr0->GetRadius() - lr1->GetRadius());
2651       if (ztol>zlims[il1]) zlims[il1] = ztol;
2652     }
2653   }
2654   //
2655   for (int il=0;il<fNActiveITSLayers;il++) printf("ZTol%d: %8.4f\n",il,zlims[il]);
2656 }
2657
2658 //_______________________________________________________
2659 double KMCDetector::PropagateBack(KMCProbe* trc)
2660 {
2661   static KMCProbe bwd;
2662   bwd = *trc;
2663   bwd.ResetCovMat();
2664   static double measErr2[3] = {0,0,0};
2665   static double meas[2] = {0,0};
2666   int icl = 0;
2667   double chi2Tot = 0;
2668   for (int il=1;il<=fLastActiveITSLayer;il++) {
2669     KMCLayer* lr = GetLayer(il);
2670     if (!PropagateToLayer(&bwd,lr,1)) return -1;
2671     int aID = lr->GetActiveID();
2672     if (aID>-1 && (icl=bwd.fClID[aID])>=-1) {
2673       KMCCluster* clMC =  icl<0 ? lr->GetMCCluster() : lr->GetBgCluster(icl);
2674       if (!bwd.PropagateToCluster(clMC,fBFieldG)) return -1;
2675       meas[0] = clMC->GetY(); meas[1] = clMC->GetZ();
2676       measErr2[0] = lr->fPhiRes*lr->fPhiRes;
2677       measErr2[2] = lr->fZRes*lr->fZRes;
2678       double chi2a = bwd.GetPredictedChi2(meas,measErr2);
2679       chi2Tot += chi2a;
2680       printf("Chis %d (cl%+3d): t2c: %6.3f tot: %6.3f\n",aID,icl,chi2a, chi2Tot);
2681       bwd.Update(meas,measErr2);
2682       bwd.AddHit(aID, chi2a, icl);
2683     }
2684     if (!bwd.CorrectForMeanMaterial(lr,kFALSE)) return -1;
2685   }
2686   return chi2Tot;
2687 }
2688
2689
2690 //////////////////////////////////////////////////////////////////////////////////////////////
2691 //--------------------------------------------------------------------------------------------
2692 ClassImp(KMCTrackSummary)
2693
2694 Int_t KMCTrackSummary::fgSumCounter = 0;
2695
2696 //____________________________________________________________________________
2697 KMCTrackSummary::KMCTrackSummary(const char* name,const char* title, int nlr) :
2698   TNamed(name,name), fNITSLayers(nlr),
2699   fPattITS(0),fPattITSFakeExcl(0),fPattITSCorr(0),
2700   fHMCChi2(0),fHMCSigDCARPhi(0),fHMCSigDCAZ(0),fHMCSigPt(0),fHMCNCl(0),fHMCFakePatt(0),
2701   fCountAcc(0),fCountTot(0),fUpdCalls(0),
2702   fRefProbe(),fAnProbe()
2703 {
2704   // create summary structure for specific track conditions
2705   //
2706   SetMinMaxClITS();
2707   SetMinMaxClITSFake();
2708   SetMinMaxClITSCorr();
2709   //
2710   if (name==0 && title==0 && nlr==0) return;  // default dummy contructor
2711   //
2712   enum {kNBinRes=1000};
2713   const double kMaxChi2(10),kMaxResRPhi=0.1, kMaxResZ=0.1,kMaxResPt=0.3;
2714   //
2715   TString strN = name, strT = title;
2716   if (strN.IsNull()) strN = Form("TrSum%d",fgSumCounter);
2717   if (strT.IsNull()) strT = strN;
2718   fgSumCounter++;
2719   //
2720   if (fNITSLayers<1) {
2721     fNITSLayers=KMCProbe::GetNITSLayers();
2722     if (fNITSLayers<1) {AliError("N ITS layer is not provided and not available from KMCProbe::GetNITSLayers()"); exit(1);}
2723     AliInfo(Form("%s No N Layers provided, set %d",strN.Data(),fNITSLayers));
2724   }
2725   nlr = fNITSLayers;
2726   //
2727   TString nam,tit;
2728   //
2729   nam = Form("%s_mc_chi2",strN.Data());
2730   tit = Form("%s MC #chi^{2}",strT.Data());
2731   fHMCChi2 = new TH1F(nam.Data(),tit.Data(),kNBinRes,0,kMaxChi2);
2732   fHMCChi2->GetXaxis()->SetTitle("#chi^{2}");
2733   fHMCChi2->Sumw2();
2734   //
2735   nam = Form("%s_mc_DCArphi",strN.Data());
2736   tit = Form("%s MC #DeltaDCA r#phi",strT.Data());
2737   fHMCSigDCARPhi = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResRPhi,kMaxResRPhi);
2738   fHMCSigDCARPhi->GetXaxis()->SetTitle("#Delta r#phi");
2739   fHMCSigDCARPhi->Sumw2();
2740   //
2741   nam = Form("%s_mc_DCAz",strN.Data());
2742   tit = Form("%s MC #DeltaDCA Z",strT.Data());
2743   fHMCSigDCAZ = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResZ,kMaxResZ);
2744   fHMCSigDCAZ->GetXaxis()->SetTitle("#Delta Z");
2745   fHMCSigDCAZ->Sumw2();
2746   //
2747   nam = Form("%s_mc_pt",strN.Data());
2748   tit = Form("%s MC $Deltap_{T}/p_{T}",strT.Data());
2749   fHMCSigPt = new TH1F(nam.Data(),tit.Data(),2*kNBinRes,-kMaxResPt,kMaxResPt);
2750   fHMCSigPt->GetXaxis()->SetTitle("#Deltap_{T}/p_{T}");
2751   fHMCSigPt->Sumw2();
2752   //
2753   nam = Form("%s_n_cl",strN.Data());
2754   tit = Form("%s N Clusters",strT.Data());
2755   fHMCNCl = new TH2F(nam.Data(),tit.Data(),nlr,1,nlr+1,nlr,1,nlr+1);
2756   fHMCNCl->GetXaxis()->SetTitle("N Clusters Tot");
2757   fHMCNCl->GetYaxis()->SetTitle("N Clusters Fake");
2758   fHMCNCl->Sumw2();
2759   //
2760   nam = Form("%s_fake_cl",strN.Data());
2761   tit = Form("%s Fake Clusters",strT.Data());
2762   fHMCFakePatt = new TH1F(nam.Data(),tit.Data(),nlr,-0.5,nlr-0.5);
2763   fHMCFakePatt->GetXaxis()->SetTitle("Fake clusters pattern");
2764   fHMCFakePatt->Sumw2();
2765   //
2766 }
2767
2768 //____________________________________________________________________________
2769 KMCTrackSummary::~KMCTrackSummary()
2770 {
2771   if (fHMCChi2)         delete fHMCChi2;
2772   if (fHMCSigDCARPhi)   delete fHMCSigDCARPhi;
2773   if (fHMCSigDCAZ)      delete fHMCSigDCAZ;
2774   if (fHMCSigPt)        delete fHMCSigPt;
2775   if (fHMCNCl)          delete fHMCNCl;
2776   if (fHMCFakePatt)     delete fHMCFakePatt;
2777 }
2778
2779
2780 //____________________________________________________________________________
2781 void KMCTrackSummary::Add(const KMCTrackSummary* src, double scl)
2782 {
2783   if (fHMCChi2)         fHMCChi2->Add(src->fHMCChi2,scl);
2784   if (fHMCSigDCARPhi)   fHMCSigDCARPhi->Add(src->fHMCSigDCARPhi,scl);
2785   if (fHMCSigDCAZ)      fHMCSigDCAZ->Add(src->fHMCSigDCAZ,scl);
2786   if (fHMCSigPt)        fHMCSigPt->Add(src->fHMCSigPt,scl);
2787   if (fHMCNCl)          fHMCNCl->Add(src->fHMCNCl,scl);
2788   if (fHMCFakePatt)     fHMCFakePatt->Add(src->fHMCFakePatt,scl);
2789   fCountAcc += src->fCountAcc;
2790   fUpdCalls += src->fUpdCalls;
2791 }
2792
2793 //____________________________________________________________________________
2794 void KMCTrackSummary::PrependTNamed(TNamed* obj, const char* nm, const char* tit)
2795 {
2796   // prepend  name, title by this prefix
2797   TString name = obj->GetName(),title = obj->GetTitle();
2798   name.Prepend(nm); title.Prepend(tit);
2799   obj->SetNameTitle(name.Data(),title.Data());
2800   //
2801 }
2802
2803 //____________________________________________________________________________
2804 void KMCTrackSummary::SetNamePrefix(const char* pref)
2805 {
2806   // prepend all names by this prefix
2807   TString nmT,nmP = pref;
2808   if (nmP.IsNull()) return;
2809   nmP = Form("%s_",pref);
2810   nmT = Form("%s ",pref);
2811   PrependTNamed(this, nmP.Data(), nmT.Data());
2812   PrependTNamed(fHMCChi2,       nmP.Data(), nmT.Data());
2813   PrependTNamed(fHMCSigDCARPhi, nmP.Data(), nmT.Data());
2814   PrependTNamed(fHMCSigDCAZ,    nmP.Data(), nmT.Data());
2815   PrependTNamed(fHMCSigPt,      nmP.Data(), nmT.Data());
2816   PrependTNamed(fHMCNCl,        nmP.Data(), nmT.Data());
2817   PrependTNamed(fHMCFakePatt,   nmP.Data(), nmT.Data());
2818   //
2819 }
2820
2821 //____________________________________________________________________________
2822 Bool_t KMCTrackSummary::CheckTrack(KMCProbe* trc)
2823 {
2824   // check if the track satisfies to selection conditions
2825   fCountTot++;
2826   if (!trc) return kFALSE;
2827   if (OutOfRange(trc->GetNITSHits(), fMinMaxClITS)) return kFALSE;
2828   if (OutOfRange(trc->GetNFakeITSHits(), fMinMaxClITSFake)) return kFALSE;
2829   if (OutOfRange(trc->GetNITSHits()-trc->GetNFakeITSHits(), fMinMaxClITSCorr)) return kFALSE;
2830   //
2831   // check layer patterns if requested
2832   UInt_t patt  = trc->GetHitsPatt();
2833   UInt_t pattF = trc->GetFakesPatt();
2834   UInt_t pattC = patt&(~pattF);    // correct hits
2835   // is there at least one hit in each of requested layer groups?
2836   for (int ip=fPattITS.GetSize();ip--;)     if (!CheckPattern(patt,fPattITS[ip])) return kFALSE;
2837   // is there at least one fake in any of requested layer groups?
2838   for (int ip=fPattITSFakeExcl.GetSize();ip--;) if (CheckPattern(pattF,fPattITSFakeExcl[ip])) return kFALSE;
2839   // is there at least one correct hit in each of requested layer groups?
2840   for (int ip=fPattITSCorr.GetSize();ip--;) if (!CheckPattern(pattC,fPattITSCorr[ip])) return kFALSE;
2841   //
2842   fCountAcc++;
2843   return kTRUE;
2844 }
2845
2846 //____________________________________________________________________________
2847 void KMCTrackSummary::AddTrack(KMCProbe* trc)
2848 {
2849   // fill track info
2850   if (!CheckTrack(trc)) return;
2851   //
2852   fHMCChi2->Fill(trc->GetNormChi2(kFALSE));
2853   fHMCSigDCARPhi->Fill(trc->GetY());
2854   fHMCSigDCAZ->Fill(trc->GetZ());
2855   if (fRefProbe.Pt()>0) fHMCSigPt->Fill( trc->Pt()/fRefProbe.Pt()-1.);
2856   //  printf("Pts: %.3f %.3f -> %.3f\n",trc->Pt(),fRefProbe.Pt(),trc->Pt()/fRefProbe.Pt()-1.);
2857   //  trc->Print("l");
2858   //  fRefProbe.Print("l");
2859   fHMCNCl->Fill(trc->GetNITSHits(),trc->GetNFakeITSHits());
2860   for (int i=trc->GetNITSLayers();i--;) if (trc->IsHitFake(i)) fHMCFakePatt->Fill(i);
2861   //
2862 }
2863
2864 //____________________________________________________________________________
2865 UInt_t KMCTrackSummary::Bits(Bool_t l0, Bool_t l1, Bool_t l2, Bool_t l3, Bool_t l4, Bool_t l5, Bool_t l6, Bool_t l7,Bool_t  l8, Bool_t l9,
2866                           Bool_t l10,Bool_t l11,Bool_t l12,Bool_t l13,Bool_t l14,Bool_t l15,Bool_t l16,Bool_t l17,Bool_t l18,Bool_t l19,
2867                           Bool_t l20,Bool_t l21,Bool_t l22,Bool_t l23,Bool_t l24,Bool_t l25,Bool_t l26,Bool_t l27,Bool_t l28,Bool_t l29,
2868                           Bool_t l30,Bool_t l31)
2869 {
2870   // create corresponding bit pattern
2871   UInt_t patt = 0;
2872   if (l0 ) patt |= (0x1<<0);    if (l1 ) patt |= (0x1<<1);   if (l2 ) patt |= (0x1<<2);
2873   if (l3 ) patt |= (0x1<<3);    if (l4 ) patt |= (0x1<<4);   if (l5 ) patt |= (0x1<<5);
2874   if (l6 ) patt |= (0x1<<6);    if (l7 ) patt |= (0x1<<7);   if (l8 ) patt |= (0x1<<8);
2875   if (l9 ) patt |= (0x1<<9);    if (l10) patt |= (0x1<<10);  if (l11) patt |= (0x1<<11);
2876   if (l12) patt |= (0x1<<12);   if (l13) patt |= (0x1<<13);  if (l14) patt |= (0x1<<14);
2877   if (l15) patt |= (0x1<<15);   if (l16) patt |= (0x1<<16);  if (l17) patt |= (0x1<<17);
2878   if (l18) patt |= (0x1<<18);   if (l19) patt |= (0x1<<19);  if (l20) patt |= (0x1<<20);
2879   if (l21) patt |= (0x1<<21);   if (l22) patt |= (0x1<<22);  if (l23) patt |= (0x1<<23);
2880   if (l24) patt |= (0x1<<24);   if (l25) patt |= (0x1<<25);  if (l26) patt |= (0x1<<26);
2881   if (l27) patt |= (0x1<<27);   if (l28) patt |= (0x1<<28);  if (l29) patt |= (0x1<<29);
2882   if (l30) patt |= (0x1<<30);   if (l31) patt |= (0x1<<31);
2883   return patt;
2884 }
2885
2886 //__________________________________________________________________________
2887 void KMCTrackSummary::Print(Option_t* ) const
2888 {
2889   printf("%s: summary for track M=%5.3f pT: %6.3f eta: %.2f\n", GetName(),
2890          fRefProbe.GetMass(),fRefProbe.Pt(), fRefProbe.Eta());
2891   printf("Cuts on NCl ITS: Tot: %2d - %2d Fake: %2d - %2d Corr: %2d - %2d\n",
2892          fMinMaxClITS[0],fMinMaxClITS[1],
2893          fMinMaxClITSFake[0],fMinMaxClITSFake[1],
2894          fMinMaxClITSCorr[0],fMinMaxClITSCorr[1]);
2895   //
2896   int nlr = fNITSLayers;
2897   if (fPattITS.GetSize()) {
2898     printf("Require at least 1 hit in groups: ");
2899     printf("Hits obligatory in groups: ");
2900     for (int i=fPattITS.GetSize();i--;) {
2901       UInt_t pat = (UInt_t)fPattITS[i];
2902       printf("[");
2903       for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
2904       printf("] ");
2905     }
2906     printf("\n");
2907   }
2908   //
2909   if (fPattITSFakeExcl.GetSize()) {
2910     printf("Fake hit forbidden in groups    : ");
2911     for (int i=fPattITSFakeExcl.GetSize();i--;) {
2912       UInt_t pat = (UInt_t)fPattITSFakeExcl[i];
2913       printf("[");
2914       for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
2915       printf("] ");
2916     }
2917     printf("\n");
2918   }
2919   //
2920   if (fPattITSCorr.GetSize()) {
2921     printf("Correct hit obligatory in groups: ");
2922     for (int i=fPattITSCorr.GetSize();i--;) {
2923       UInt_t pat = (UInt_t)fPattITSCorr[i];
2924       printf("[");
2925       for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
2926       printf("] ");
2927     }
2928     printf("\n");
2929   }
2930   //
2931   printf("Entries: Tot: %.4e Acc: %.4e -> Eff: %.4f+-%.4f %.2e KMC updates\n",fCountTot,fCountAcc,GetEff(),GetEffErr(),GetUpdCalls());
2932 }
2933